MetaD box size variation in NPT condition

[Alexandre Oliveira]

Dear all,

I'm been using metadynamics, in this case, transition-tempered metadynamics, to study permeations of solutes through a membrane. My two CVs are distances of different parts of the solute to the center of a membrane. Then I realized huge jumps in the box size in z coordinate, as you see in this plot:

I'm using plumed 2.5.7 and from the manual since version 2.0 that plumed is able to perform bias simulation in the NPT conditions. My plumed input is attached. Am I missing something? have I done anything wrong?  Should I use a more recent version? I can add UPPER_WALLS and LOWER_WALLS in the cell size but I'm concerned that I'm no longer in NPT conditions.

Best regards,

Alexandre Oliveira

[Giovanni Bussi]

Hi,

which MD code are you using?

If the jump is in the box size (and not in the solute coordinate) there is something very weird happening. How can the box have this big change without the sytem exploding?

If instead the jump is in the scaled coordinates (CV1.c or CV2.c) then it is a PBC problem.

Giovanni

--
You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/27e34530-d0b9-4ba0-b559-ec1c093da8b8n%40googlegroups.com.

[Alexandre Oliveira]

Hi Professor,

I'm using GROMACS 2019.4. The box doesn't explode because it expands in x,y coordinate. 

This happens when my two CVs are distances and the height of the gaussian is 1.2 kJ/mol at lower values it doesn't seem to happen. Everything seems fine with the solute coordinate:

With other CVs, like density of water in certain zones of the box, it happens also something to the box, but I always assumed that it was because of the type of CV, maybe it's the same issue. When I visualize the trajectory, the system compacts in z, but now that you mention the PBC, it does look like there is some correlation between the solute crossing the frontier of the system and the jump in the box size, but I'm speculating. Could it be that this is happening because I'm using: WHOLEMOLECULES ENTITY0=... ?

Best regards,

Alexandre Oliveira

[Alexandre Oliveira]

Hi,

I've tested with plumed 2.8, the same plumed input and I get the same jumps in the box:

I've check the plumed manual, I don't see any mistakes in the input. In the gromacs inputs, I also don't see any mistakes. I've run a normal unbias simulation without plumed, I didn't get any issue. With other parameters in plumed input like lower height of the gaussian I also didn't get any issue. Could it be some bug in the code? I would really appreciate any help.

Thank you in advance,

Best regards,

Alexandre Oliveira

[Giovanni Bussi]

Hi.

If the system is dense (if you have a solvent around) it should explode in this regime. There is no way to have a jump on the cell side, even if compensated by a jump leading to the same volume on the other cell side.

Anyway, if the time scale is such that these are not instantaneous jumps, but you are deforming the box gradually, then the explanation could be the following.

Using relative coordinates (SCALED_COMPONENTS), you can increment a scaled coordinate both incrementing the real coordinate and decrementing the cell side. 

In your input file, there is a line that could be affected by jumps across PBCs: b1 is computed on a set of atoms for which the PBC image might change. Maybe you have jumps in this variable leading to strong forces that then deform the box for the reason above.

Notice that if b1 jumps by a lattice vector there is no problem, since it is used in a scaled component only. But if only some of the atoms used to compute b1 jump, then b1 will jump by a fraction of a lattice vector, leading to problems.

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/e0206dfb-e7ae-4f01-8d6a-bedc18b3a741n%40googlegroups.com.

[Alexandre Oliveira]

Thank you for the reply.

ok, I understand the situation. You are right, b1 is in fact the membrane atoms, and some of the lipids jump across PBC. Since I want the distance between the solute and membrane COM, so I shouldn't be using the scaled_components or if I want to use scaled_components I should rebuild all the lipids in each step (with wholemolecules)?

Best regards,

Alexandre Oliveira

[Giovanni Bussi]

Yes, exactly.

Notice that you would be able to reconstruct correctly the center of mass even using only COM without NOPBC (and without WHOLEMOLECULES).

You can also use less atoms to define the membrane COM, to make the calculation faster.

Finally, I guess (but it depends on the box size relative to the membrane size) that the order in which you place the (subsets of) lipid atoms is not relevant. If the COM position is wrong along x and y, there will be no effect on CV1.c and CV2.c. I guess that just using atom numbering order is fine. However, you should double-check carefully the box size relative to the membrane size. The lipids are usually oriented, if you pick atoms close to the center of the membrane it is less likely that this problem will show up, but on the other hand you will have a less homogeneous distribution of the reaction force on the membrane. You might construct lists with some script, something like:

- go through the first lipid molecule, possible skipping atoms, from the lower to the upper extreme (in z direction)

- then go through the second lipid molecule in the opposite order

- then go through the third lipid molecule in the same order as in the first one.

- etc

Make sure you understand how PBC reconstruction works in PLUMED because it's easy to make mistakes...

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/08c0899e-c00c-4c22-8220-be3ec91c1a38n%40googlegroups.com.

[Alexandre Oliveira]

Thank you for the reply,

I've tested with different options, with WHOLEMOLECULES, WRAPAROUND, with or without the NOPBC option. The only situation where I didn't get the Box compression was when I removed the SCALED_COMPONENTS, but I still get a very high change in the box dimensions:

That makes me think that perhaps, with biased simulations, I will always have something like this and that this is the best I can get. 

I'm wondering if the height of the gaussian as 1.3 kJ/mol isn't too high? The translocation barrier is very height that's why I've used 1.3 kJ/mol. Since you always recommend using half KbT, for me made sense to use 1.3 Kj/mol but I see from the literature that the height of the gaussian used for permeations processes is always smaller like for example 1/100 kbT. I'm wondering if the height of the gaussian as 1.3 kJ/mol isn't too high and if it is the reason for these fluctuations in the box size?

Best regards,

Alexandre Oliveira