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