Distance Components and correct periodic behavior

[lucie delemotte]

Hi guys,

I'm trying to make sense of the correct way to account for periodicity when computing distance components running with gromacs.

I am using the latest development version of plumed patched to gromacs 4.6.3.

I would like to bias the z component of a distance between a real particle and a ghost particle. The particle belongs to a protein and the ghost is defined with respect to 3 atoms in the same protein.

I understand I should specify before anything else in my plumed input

WHOLEMOLECULES ENTITY0=1-2619 (1-2619 being the protein atoms) so that the protein is considered as a whole

I then define my ghost particle 

GHOST ATOMS=286,936,1707 COORDINATES=-0.424,0.659,-0.780 LABEL=G1

and my regular particle

A: GROUP ATOMS=2264

I then calculate the distance between the 2 and its 3 cartesian components with

DISTANCE ...

   LABEL=A_G1

   ATOMS=A,G1

   COMPONENTS

... DISTANCE

This does not give the right behavior, along the trajectory, the z component of the distance is sometimes more than half the size of the box (I did get a warning in the output and from the manual that the component will not have the right behavior so no surprise there). I understand from the manual I should be using:

dA_G1: COMBINE ARG=A_G1.z PERIODIC=-Z,+Z

However, using pressure control, the size of the simulation box is fluctuating making it impossible to define the value of the z axis of the box. 

It doesn't seem to me that gromacs has a keyword that allows to keep molecules whole either.

What do you guys recommend ? I would like to stick with NPT if possible.

Best,

Lucie

[Giovanni Bussi]

Hi Lucie,

On Wed, Jan 22, 2014 at 1:23 AM, lucie delemotte <lucie.d...@gmail.com> wrote:

Hi guys,

I'm trying to make sense of the correct way to account for periodicity when computing distance components running with gromacs.

I am using the latest development version of plumed patched to gromacs 4.6.3.

I would like to bias the z component of a distance between a real particle and a ghost particle. The particle belongs to a protein and the ghost is defined with respect to 3 atoms in the same protein.

I understand I should specify before anything else in my plumed input

WHOLEMOLECULES ENTITY0=1-2619 (1-2619 being the protein atoms) so that the protein is considered as a whole

Be careful: this could slow down your simulation. If you can reduce this list of atoms to a list such that

1. consecutive atoms are always closest than half box

2. atoms you are interested in are part of the list

you should have the same effect at a lower cost (e.g., you could use the ordered list of all the c_alpha, plus atoms 286,936,1707, 2264).

 

I then define my ghost particle 

GHOST ATOMS=286,936,1707 COORDINATES=-0.424,0.659,-0.780 LABEL=G1

Another warning here (I thought it was in the manual, but actually we forgot to include it... I will asap): GHOST_ATOMS is possibly the only CV for which virial contribution is not computed correctly, so you could have artifacts in NPT simulations. (technical explanation here: github.com/plumed/plumed2/issues/74). I am thinking at very tiny artifacts, but if you can manage to have an equivalent virtual atom using CENTER that would be rigorously correct.

 

and my regular particle

A: GROUP ATOMS=2264

I then calculate the distance between the 2 and its 3 cartesian components with

DISTANCE ...

   LABEL=A_G1

   ATOMS=A,G1

   COMPONENTS

... DISTANCE

This does not give the right behavior, along the trajectory, the z component of the distance is sometimes more than half the size of the box (I did get a warning in the output and from the manual that the component will not have the right behavior so no surprise there).

Another comment: it looks your regular particle (2264) is in the protein. So, are you sure you want the z component to be periodic? You might want to use the NOPBC flag: component will be not periodic and you can stop reading here.

Now, I image you want to include cell periodicity.

The only situations where A_G1.z could be outside -Z/2,+Z/2 is:

1. If cell is not orthorombic (this is complicated, see below).

2. If cell is not fixed (NPT), since you cannot define Z (as you noticed).

 

I understand from the manual I should be using:

dA_G1: COMBINE ARG=A_G1.z PERIODIC=-Z,+Z

This is ok for NVT (notice that Z here should be *half* box size). Notice that PERIODIC here is not meant at bringing the variable in the right domain (it should already be there). The reason to use it is that if you use e.g. metadynamics and add a hill at z close to plus half box, you want it to be felt also when z is close to minus half box.

 

However, using pressure control, the size of the simulation box is fluctuating making it impossible to define the value of the z axis of the box. 

It doesn't seem to me that gromacs has a keyword that allows to keep molecules whole either.

What do you guys recommend ? I would like to stick with NPT if possible.

The only solution with NPT (that I know) is to use scaled coordinates. This is on the todo list (github.com/plumed/plumed2/issues/63). Meanwhile, as a hack, you can "implement" it in your input using something like

cell: CELL

dist: DISTANCE ATOMS=A,B COMPONENTS

distz: MATHEVAL ARG=dist.z,cell.cz FUNC=x/y PERIODIC=-0.5,0.5

Notice that distz will be in scaled units, you can convert it back to nm using a "nominal" cell size (this is just to simplify interpretation of the value).

For non-orthoromic boxes, you need to scale using the full cell matrix. That should also be feasible with MATHEVAL, but is more complicated.

Let us know if this solution works for you!

Giovanni

 

Best,

Lucie

--
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 post to this group, send email to plumed...@googlegroups.com.
Visit this group at http://groups.google.com/group/plumed-users.
For more options, visit https://groups.google.com/groups/opt_out.

[lucie delemotte]

Thanks for your answer !

I think I understand better now:

If you specify wholemolecules and consider distance between atoms that are either within the wholemolecule block or virtual particles based on such atoms, you can specify NOPBC because the molecule is considered as a whole anyways. Is that correct ?

I have not tried the hack you proposed since it seems that a few tests with nopbc worked for me.

Thanks also for all the warnings, I think I now have all the elements to get to work !

Best,

--
Lucie