getting the correct cartesian component distance between two groups (or, is matheval necessary)?

[Chris Neale]

Dear Plumed users:

I am trying to set up a collective variable that will bring a sodium ion across an ion channel. I've been following the example for that on the plumed webpage ( http://plumed.github.io/doc-v2.2/user-doc/html/_d_i_s_t_a_n_c_e.html ), but I think that I am running into problems because by arguments to PERIODIC within a COMBINE command are functions rather than constants (since I'm running with pressure coupling and therefore a box dimension that is not constant). I'm also heavily influenced by the fact that the plumed page noted above contains the text "Cartesian components will not have the proper periodicity", which I presume to mean that if I have a 10-nm box and two groups are 1 nm apart, then plumed COMPONENTS might report 1-nm or 11-nm (even though I have seen the value be correct in a very short test) whereas SCALED_COMPONENTS is guaranteed to be correct.

I first tried this:

WHOLEMOLECULES ENTITY0=1-4857

protein: COM ATOMS=1-4857

sodium: GROUP ATOMS=5341

distz: DISTANCE ATOMS=protein,sodium COMPONENTS

cell: CELL

halfboxz: COMBINE ARG=cell.cz COEFFICIENTS=0.5 POWERS=1 PERIODIC=NO

neghalfboxz: COMBINE ARG=cell.cz COEFFICIENTS=-0.5 POWERS=1 PERIODIC=NO

dz: COMBINE ARG=distz.z PERIODIC=neghalfboxz,halfboxz

PRINT ARG=dz FILE=look.plumed.dz.xvg

Which gives the error:

PLUMED: ERROR in input to action COMBINE with label dz : could not convert period string neghalfboxz to real

###########

So I tried this:

WHOLEMOLECULES ENTITY0=1-4857

protein: COM ATOMS=1-4857

sodium: GROUP ATOMS=5341

dist: DISTANCE ATOMS=protein,sodium SCALED_COMPONENTS

cell: CELL

dz: MATHEVAL ARG=dist.c,cell.cz FUNC=x*y PERIODIC=NO

PRINT ARG=dz FILE=look.plumed.dz.xvg

Which gives the error:

PLUMED: I cannot understand line: MATHEVAL LABEL=dz ARG=dist.c,cell.cz FUNC=x*y PERIODIC=NO

and I'm having trouble compiling plumed with libmatheval, but I wonder if this is the only way forward?

Thank you,

Chris.

[Giovanni Bussi]

Hi Chris,

a few comments on your input:

1. WHOLEMOLECULES ENTITY0=1-4857

If you just need it to compute the following:

  protein: COM ATOMS=1-4857

then you can remove it. WHOLEMOLECULES is not required (plumed does it on the fly since v2.2 when you compute the COM).

2. COM with 5000 atoms

This will slow down your simulation. Replace it with something like 1-4857:20 or better pick a set of meaningful atoms (I suggest all calpha, you can do it with a script).

3. dz: COMBINE ARG=distz.z PERIODIC=neghalfboxz,halfboxz

Domain should be defined in terms of constants. The quick trick is to set it by hand to the its typical value (e.g. -5,+5 if z dimension is ~10nm). However, this will introduce discontinuities on the distance when the sodium is at a distance close to +5 or -5. If you are 100% sure this will never happen (e.g. since you are using a wall to limit the exploration of the ion) you can ignore this warning and proceed. Otherwise the only solution is to use SCALED_COMPONENTS.

4. dz: MATHEVAL ARG=dist.c,cell.cz FUNC=x*y PERIODIC=NO

This is correct if you want to compute the value of dz to monitor it. However, if you apply a force on it (e.g. with walls or metad) you will still be in trouble. You are forcing plumed to assume that dz is not periodic, but actually it has a periodicity of neghalfboxz,halfboxz, and non-constant periodicities are not supported.

If you want to use the z position as an argument to a bias, you should use directly dist.c, which has period -0.5,+0.5.

If you don't like working with unitless variables (e.g. in METAD you want to set SIGMA in nanometers) you can (approximately) convert it to nm with COMBINE:

dz: COMBINE ARG=dist.c COEFFICIENTS=XXX PERIODIC=-YYY,YYY

where XXX is a precomputed typical value of cell.cz and YYY half of it. In an NPT simulation the actual value is going to fluctuate by max 1%. Notice that no discontinuity will be introduced in this way. The only artefact is that you will think distance is say 0.5 nm and it will be 0.505. Not a big problem.

Hope this helps.

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 post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/de339878-6cd2-492a-97f1-d4acb020b43b%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Chris Neale]

Dear Giovanni:

thank you for the suggestions, I will certainly try this. However, the variance is more than you suspect. It would be less with regular pressure coupling, but with membranes it is common to use semi-isotropic pressure coupling in which pressure along z is coupled separately from pressure in the xy plane (to allow the area per lipid to vary). Below are some numbers from a 700 ns simulation of the protein in the membrane (no plumed, just standard simulation):

box z = 11.51 +/- 0.12 nm

min = 11.13 nm

max = 11.90 nm

range = 0.77 nm = 6.7 %

For metadynamics, I suspect that this 3.8 angstrom uncertainty in the position of the ion (half of the 7.7 angstrom variance in box height) is going to introduce some unwanted noise. However, I admit that this 3.8 angstrom value is the absolute largest such deviation and even then it's really more like +/- 1.9 angstroms about the average. I suspect that I will have to fix the z-coordinate, which is probably not such a big deal anyway.

Thanks again for the help with the functions, I really appreciate it (as always).

Chris.

[Chris Neale]

Dear Giovanni:

one more question:

My box is 11 nm tall and I'm looking only at the ion near the center of the channel. I'm going to use upper walls and lower walls to keep the ion within a 3 nm region near the center of the channel. I'm then going to use metadynamics to bias sampling along z within this range. I wonder if I interpret your previous comment correctly that in fact it is entirely correct in this case to simply use dist.z (i.e., COMPONENTS rather than SCALED_COMPONENTS). For the system that I describe, I can't see any reason that there should be any issue with periodicity unless it is possible that when the ion is e.g., in the exact center of the protein plumed will sometimes report d.z as 0 and other times will report d.z as +/- the box height in z.

Thank you,

Chris.

[Giovanni Bussi]

Hi,

that's correct. z component of DISTANCE is going to be computed with the correct pbc (unless you add NOPBC to DISTANCE), so with walls you should have no 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 post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/921640b7-7620-4179-aaa7-cd8d7a0afb4b%40googlegroups.com.

[Sofya Lushchekina]

Dear Giovanni,

> If you don't like working with unitless variables (e.g. in METAD you want to set SIGMA in nanometers) you can (approximately) convert it to nm with COMBINE:

I have a question about SIGMA in METAD with SCALED_COMPONENTS. Do I have to scale SIGMA too?

E.g. I have GRID_MIN=-5 GRID_MAX=5 SIGMA=0.1, do I convert it to GRID_MIN=-0.5 GRID_MAX=0.5 SIGMA=0.01 with SCALED_COMPONENTS?

Thank you in advance,

Sofya