[gmx-users] FEP
David Mobley
dmobley at gmail.com
Thu Mar 30 17:23:01 CEST 2006
P,
I am not sure what's going on; it could be any number of things.
However, I would probably start by trying to create a topology which
will work for free_energy=no (which corresponds to lambda=0 which
corresponds to no perturbation (at least in the case when soft core is
not used)). If your topology works properly then (which it should),
then there are a couple of other things I would suggest.
1) sc-alpha IMO is usually better at 0.5 than 1.5 for particle
annhilation, at least with sc-power=1.0 (the default). I've done some
checking on this and this seems to generally be the optimal choice of
parameters (which is also what Michael Shirts has found). You might
check if changing this helps. (In my experience sc-alpha=1.5 can lead
to really steep slopes in dv/dlambda near lambda=0 and lambda=1, which
could lead to large forces and some of the instabilities you seem to
be seeing).
2) There are a number of people who do exactly what you're doing and
turn off charges and van der Waals interactions at the same time, but
I don't think it's really correct to do so. First, you want a smooth
dv/dlambda curve, but if you adjust soft core parameters to give you a
relatively smooth curve for turning off van der Waals interactions, it
gives you a very un-smooth curve for turning off charge interactions,
meaning that if you do both of them at the same time, your curve won't
be very smooth. In principle, there can be worse problems, too: Soft
core allows atoms to start interpenetrating somewhere before lambda=1,
when the charges are still on to some extent. So now you have a Cl-
with very weak vdW interactions such that other atoms can overlap
with it -- atoms like, for example, the hydrogens on water (which have
the opposite charge and typically no vdW radii either). So now you're
putting a positively charged atom almost on top of a (somewhat)
negatively charged ion, and bad things happen (usually the simulation
blows up. You could describe this as sort of a fusion event... :) ).
My typical approach is to turn off charges first in a separate
calculation without using soft core. The electrostatic interactions
are usually a very smooth function of lambda, so this can usually be
done with only 5 lambda values or so. THEN use soft core to turn off
the vdW interactions. Again, there are two advantages to this:
(a) Although it IS harder to run two sets of calculations instead of
one, you can usually get by with fewer lambda values on each since now
both transformations are relatively smooth functions of lambda. So
it's two sets of calculations but you may be able to do them with the
same total number of lambda values.
(b) You avoid charge-charge overlap.
IMO your problem is more likely to be the first than the second, since
you are seeing problems at lambda=0 as well. Do let me know what you
find out, though, especially if it is the second problem (I need to
demonstrate the charge-charge overlap problem, but it is only sporadic
in everything I've tried, so I'd like to find a system where it's
reproducible).
Incidentally, are you also using PME for long range electrostatics?
Best wishes,
David Mobley
UCSF
On 3/30/06, P <pyt1 at op.pl> wrote:
>
> Hi all.
> I'm trying FEP In gromacs3.3. My system (3nm; 3nm; 3nm) consists
> of one Na, Cl, DUM ion and solvent molecules.
> During my free energy calculations I want to perturb:
> 1) Cl- into a DUM atom
> 2) DUM into Cl-
>
>
> I'm using ffgmx. I've added new atom DUM into ffgmx.atp:
> …..
> SD 32.06000 ; DMSO Sulphur
> OD 15.99940 ; DMSO Oxygen
> CD 15.03500 ; DMSO Carbon
> CHE 12.01100 ; HEME RING CARBON
> MNH3 0 ; Dummy mass in rigid tetraedrical NH3 group
> MW 0 ; Dummy mass in rigid tyrosine rings
> IW 0 ; Dummy particle in TIP4P etc.
> DUM 0 ; moj atom
> …..
>
> I've updated ffgmxnb.itp:
> ….
> OWT3 15.99940 0.000 A 0.24889E-02 0.24352E-05
> SD 32.06000 0.000 A 0.10561E-01 0.21499E-04
> OD 15.99940 0.000 A 0.22715E-02 0.75147E-06
> CD 15.03500 0.000 A 0.90507E-02 0.21758E-04
> CHE 12.01100 0.000 A 0.23402E-02 0.33740E-05
> MNH3 0 0.000 A 0.0 0.0
> MW 0 0.000 A 0.0 0.0
> DUM 0 0.000 A 0.0 0.0
> ….
>
> My dual topology for Cl and DUM atoms:
>
> [ moleculetype ]
> ; molname nrexcl
> Cl 1
>
> [ atoms ]
> ; id at type res nr residu name at name cg nr charge
> 1 DUM 1 DUM DUM 1 0 35.45300
> Cl -1 35.45300
>
>
> … and
>
> [ moleculetype ]
> ; molname nrexcl
> ClB 1
>
> [ atoms ]
> ; id at type res nr residu name at name cg nr charge
> 1 Cl 1 Cl Cl 1 -1 35.45300
> DUM 0 35.45300
>
>
> I did calculations for 20 lambda values 0.05, 0.10 , 0.15 ….
> 0.95 and for those values delta free energy looks reasonable.
> My problem is that I get error every time I try to run
> FEP for lambda = 1 or lambnda = 0.
> Am I doing something wrong. Please give me some advice.
> (I'm using PME, DUM and Cl atom are restrained, step 1fs)
>
>
> My ERROR message for lambda = 0 is :
> …
> ….
>
> Large VCM(group Cl): -38.08457, 8.92837, -19.97197, ekin-cm:
> 6.83899e+04
> Large VCM(group Na): 0.04709, -0.21100, -0.29647, ekin-cm:
> 1.54762e+00
>
> t = 0.002 ps: Water molecule starting at atom 61 can not be settled.
> Check for bad contacts and/or reduce the timestep.Wrote pdb files with
> previous and current coordinates
> Large VCM(group SOL): -0.01922, -0.00999, -0.02540, ekin-cm:
> 8.98334e+00
> Large VCM(group Cl): 4.35556, 2.33494, 5.86521, ekin-cm:
> 2.08547e+03
> Large VCM(group Na): 0.04590, -0.19717, -0.27716, ekin-cm:
> 1.35412e+00
> Large VCM(group SOL): -0.05937, 0.00317, -0.07617, ekin-cm:
> 7.52736e+01
> Large VCM(group Cl): 13.48634, -0.66233, 17.40376, ekin-cm:
> 1.72022e+04
> Large VCM(group Na): 0.04470, -0.18096, -0.25439, ekin-cm:
> 1.14329e+00
> Large VCM(group SOL): -0.00393, -0.03332, 0.01434, ekin-cm:
> 1.07341e+01
> Large VCM(group Cl): 0.88076, 7.63229, -3.18263, ekin-cm:
> 2.45181e+03
>
> t = 0.003 ps: Water molecule starting at atom 631 can not be settled.
> Check for bad contacts and/or reduce the timestep.Wrote pdb
> files with previous and current coordinates
> Large VCM(group Cl): 6.73354, 14.12955, -14.62614, ekin-cm:
> 1.62697e+04
>
>
> Thank you for your help
> ALL BEST
>
> _______________________________________________
> 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.
> Can't post? Read
> http://www.gromacs.org/mailing_lists/users.php
>
>
More information about the gromacs.org_gmx-users
mailing list