[gmx-users] problem calculating PMF with distant constraints

David van der Spoel spoel at xray.bmc.uu.se
Fri Dec 16 09:15:39 CET 2005

```>
>
> We used the same combination rule as the reference and also used the second type of combination rule.
>
> Below is one of the topology files. Similar results we got also with SPC, SPCE and tip5p, and also with different LJ particles.
>
> The locations of the minimums of the PMF are good but the depth at contact is higher than the state at large distances, where is should be lower.
>

I assume you use the normal free energy code for this, then. The output
you get is dH/dlambda, and since you have carefully tuned the A and B to
mean lambda really corresponds to distance in nm, the dH/dlambda has the
dimension of a force. So far so good.

The only thing I can think of is a PBC problem. Your box is roughly 3.5
nm I presume, is your LJ cut-off larger than 1.22 nm?

Have you tried the alternative of freezing the neon positions and
(without explicit use of free energy code) compute the average force
between the particles as a function of distance? It should give the same
result. If not then something may be wrong with the free energy code, in
which case you should submit a bugzilla.

>
>>
>>>
>>>>>In particular when the constraint distance was bigger than the cut-off
>>>>
>>>>interaction (so the two atoms don't `see' each other) and bigger than the
>>
>>range
>>
>>>>of the solvent induced interactions, the constraint force was not zero but
>>>>replusive. This can be related to what was mentioned in the literature before
>>>>[Chem. Phys. Lett. 196 297 (1992), JCP 100 9129(1994)] as a `dynamic'
>>>>contribution to the free energy when the system allowed to rotate (the
>>>>centripetal force). Thus, even when there are no forces on the atoms, but ones
>>>>allow distance constraint to rotate the atoms will experience a force in the
>>>>opposite direction to their center of mass, which, of course, should not be in
>>>>the PMF - is there such correction to this in Gromacs?
>>>>
>>>>correction.
>>>>
>>>>
>>>>
>>>>>Nevertheless, if one restrains the two atoms from moving along the x and y
>>>>
>>>>axis and sets the constraint distance along the z-axis it still doesn't give
>>
>>the
>>
>>>>
>>>>
>>>>>Did someone reported these issuses before
>>>>
>>>>?
>>>>I'm not completely sure that you used the right method. If you put a
>>>>constraint between the atoms they won't interact at all, unless you set
>>>>nrexcl to zero. An alternative method is to freeze the Neon atoms and
>>>>measure the force on them.
>>>
>>>
>>>nrexcl was set to zero and fn. type 2 for the constraint `bond' was used, so
>>
>>they should see each other.
>>"should" is not good enough. have you checked the exclusions in the tpr
>>file?
>
>
> this is the first few lines of the excl groups from the top file (atom 0 and atom 1 are Neon the rest is water):
>
>             excl[0][0..0]={0}
>             excl[1][1..1]={1}
>             excl[2][2..4]={2, 3, 4}
>             excl[3][5..7]={2, 3, 4}
>             excl[4][8..10]={2, 3, 4}
>
>
>>>>>This was checked with 3.2.1 and 3.1.4 versions, with a few different water
>>>>>models, PME and cut-offs for the electrostatic interactions as well.
>>>>>
>>>>>Is there a simple fix? The group I am working in is quite concerned about
>>>>
>>>>this.
>>>>If you have a well defined system where you think you are doing exactly
>>>>the same as in a publication, and you still can not reproduce it then
>>>>feel free to post a bug report on bugzilla.gromacs.org. I assume you are
>>>>reasonably sure that your reference values are correct.
>>
>>This comment still holds. For your PMF it could e.g. be a single point.
>>
>>
>>
>>--
>
>
>
>
> #include "ffoplsaa.itp"
>
> [ atomtypes ]
> ;name   mass      charge   ptype   sigma        epsilon
> LJP   20.1797     0.000     A     0.2749        0.295108
>
>
> [ moleculetype ]
> ; Name            nrexcl
> LJG                 0
>
> [ atoms ]
> ;nr  type  resnr residue  atom   cgnr  charge  mass
> 1    LJP     1     LJG    LJ1     1     0.0    20.1797
> 2    LJP     1     LJG    LJ2     1     0.0    20.1797
>
> [ constraints ]
> ; #at.   f.tp  length_A  Length_B
> 1   2     2     0.220    1.220
>
>
> [ moleculetype ]
> ; molname       nrexcl
> SOL             2
>
> [ atoms ]
> ; id    at type res nr  residu name     at name         cg nr   charge
> #ifdef _FF_OPLS
> 1     opls_111  1       SOL              OW             1       -0.834
> 2     opls_112  1       SOL             HW1             1        0.417
> 3     opls_112  1       SOL             HW2             1        0.417
> #else
> 1       OWT3    1       SOL              OW             1       -0.834
> 2       HW      1       SOL             HW1             1        0.417
> 3       HW      1       SOL             HW2             1        0.417
> #endif
>
> #ifdef FLEXIBLE
> [ bonds ]
> ; i     j       funct   length  force.c.
> 1       2       1       0.09572 502416.0 0.09572        502416.0
> 1       3       1       0.09572 502416.0 0.09572        502416.0
>
>
> [ angles ]
> ; i     j       k       funct   angle   force.c.
> 2       1       3       1       104.52  628.02  104.52  628.02
>
> #else
> [ settles ]
> ; i     j       funct   length
> 1       1       0.09572 0.15139
>
> [ exclusions ]
> 1       2       3
> 2       1       3
> 3       1       2
> #endif
>
>
> [ system ]
> ; Name
> LJG   SOL
>
> [ molecules ]
> ; Compound        #mols
> LJG                 1
> SOL                1200
> _______________________________________________
> gmx-users mailing list
> gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.

--
David.
________________________________________________________________________
David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,  	75124 Uppsala, Sweden
phone:	46 18 471 4205		fax: 46 18 511 755
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://xray.bmc.uu.se/~spoel
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

```