Jonathan,<br><br><div><span class="gmail_quote">On 11/11/05, <b class="gmail_sendername">Moore, Jonathan (J)</b> <<a href="mailto:JMoore2@dow.com" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">JMoore2@dow.com
</a>> wrote:</span><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
David and David,<br><br>Thanks for the suggestions. I'm not currently using soft cores for fear of<br>the soft core bug in 3.2.1.<br><a href="http://www.gromacs.org/pipermail/gmx-users/2004-August/011830.html" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">
http://www.gromacs.org/pipermail/gmx-users/2004-August/011830.html
</a><br><br>I haven't moved to 3.3 yet on any of my production machines because I'm<br>holding out for 3.3.1, but maybe I'll have to give up waiting so I can use<br>soft cores.</blockquote><div><br>
You _could_ just change that bit of the source code as described in
that post and recompile. I was using a CVS version of 3.2.1 from Nov.
2004 which had that bugfix until very recently and the softcore stuff
was working fine for me.<br>
<br>
It's worth noting, though, that there are number of other bugfixes
implemented in 3.3 which affect free energy calculations; it's not
clear to me whether the version I was using of 3.2.1, or the version
you are using, had all of these correctly implemented. <br>
<br>
And just by way of advertisement, I'm currently using 3.3 and have been
able to compute free energies for turning off charges and vdW for a
number of small molecules in water. I've done so with the AMBER force
field (from the Pande group) and my results match what I was getting
with the AMBER package. This gives me a reasonable amount of confidence
that things are working properly. I've also computed the charging free
energy of spc water in water, the answer of which is known, and I get
the correct answer. <br>
<br>
I have, however, included this bugfix
<a href="http://www.gromacs.org/pipermail/gmx-users/2005-October/017495.html">http://www.gromacs.org/pipermail/gmx-users/2005-October/017495.html</a> in
my 3.3 installation. It's not clear to me that without the bugfix I
necessarily would have had problems, but...<br>
<br>
</div><br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">David M., I guess this is the paper you were talking about, right?<br>T. C. Beutler, A. E. Mark, R. C. van Schaik, P. R. Gerber and W. F. van
<br>Gunsteren. "Avoiding singularities and numerical instabilities in free<br>energy calculations based on molecular simulations." Chem. Phys. Lett., 222,<br>529-539 (1994).<br> <a href="http://dx.doi.org/10.1016/0009-2614%2894%2900397-1" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">
http://dx.doi.org/10.1016/0009-2614(94)00397-1</a></blockquote><div><br>
Yep, that's the one.<br>
<br>
David<br>
<br>
</div><br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Thanks,<br>Jonathan<br>____________________________<br>Jonathan Moore, Ph.D.<br>
Research and Engineering Sciences - New Products<br>Core R&D<br>The Dow Chemical Company
<br>1702 Building, Office 4E<br>Midland, MI 48674 USA<br>Phone: (989) 636-9765<br>Fax: (989) 636-4019<br>E Mail: <a href="mailto:jmoore2@dow.com" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">jmoore2@dow.com
</a><br><br>-----Original Message-----<br>From: David Mobley [mailto:
<a href="mailto:dmobley@gmail.com" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">dmobley@gmail.com</a>]<br>Sent: Friday, November 11, 2005 3:27 PM<br>To: Discussion list for GROMACS users<br>Subject: Re: [gmx-users] free energy calculations blowing up
<br><br><br>Jonathan,
<br><br>To clarify slightly, the singularity for lambda->1 is when the lambda=1<br>Hamiltonian corresponds to a completely disappeared particle. If you are<br>changing on vdW radii into another vdW radii things are somewhat more
<br>complicated and it isn't clear to me that you really have the same sort of<br>problem. I wouldn't be surprised if you did, though, but it would be at some<br>intermediate lambda value. I would still try the soft core option,
<br>especially if your time series of dV/dlambda have anything that looks at all<br>like shot noise.<br><br>David<br><br><br><br>On 11/11/05, David Mobley <<a href="mailto:dmobley@gmail.com" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">
dmobley@gmail.com</a>> wrote:
<br>Jonathan,<br><br>Two comments:<br>(1) Although this wasn't your question, you might not want to run the<br>simulations sequentially. I think previous work has shown that you can<br>introduce some hysteresis this way, although the reference eludes me at the
<br>moment.<br>(2) When changing vdW parameters, soft core helps a lot. From what you're<br>describing it sounds like you are doing that. If you don't use soft core,<br>dV/dlambda actually has a singularity for particle insertion as lambda->1
<br>(assuming the exponent of lambda is 1, which it is by default) and is also<br>very susceptible to shot-type noise, and susceptible to crashes when forces<br>get too large (which often happens near lambda=1). You can reduce the
<br>crashes somewhat by using smaller timesteps but this doesn't help with the<br>shot noise or singularity problems. To understand why this happens you can<br>look up van Gunsteren's paper from around 1992 introducing the soft core
<br>potential, or contact me directly for a brief explanation. If you have any<br>trouble finding the reference I can look it up for you, but I don't have it<br>in front of me right this instant.<br><br>David Mobley<br>UCSF
<br><br><br><br><br><br>On 11/11/05, Moore, Jonathan (J) <<a href="mailto:JMoore2@dow.com" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">JMoore2@dow.com</a> > wrote:<br><br>I have a question about free energy calculations.
<br><br>I'm performing several different simulations where the hydrogen of certain
<br>hydroxyls on my molecule is replaced with a methyl group. It's a united<br>atom model, so the change from hydrogen to methyl doesn't change the number<br>of sites. I have a single molecule of interest in SPC water (see more
<br>details below). I'm experiencing the problem that on occasion the system<br>blows up. I've had no problems like this previously when simulating only<br>one state or the other (methylated or unmethylated), so I think the problem
<br>is due to the free energy calculation. Unless I've made an error in setting<br>up the free energy calculation, I assume my problem is that the system isn't<br>stable for the timestep that I'm using (2 fs) and how fast I'm changing
<br>lambda (I'm incrementing lambda by 0.05 followed by 200 ps of simulation<br>before increasing it again).<br><br>Is it common to either decrease the time step or use lambda increments<br>smaller than 0.05 when doing free energy calculations this way? I'm
<br>planning to try both to see if they fix the problem, but I'm curious if<br>anyone else has had a problem like this. I'm not using soft cores, but<br>maybe I have to?<br><br>Also, I'm specifying the atom types like this (
i.e., with the mass being<br>taken from the .atp file):<br><br>[ atoms ]<br>;
nr type
resnr residu atom cgnr charge<br>
6
H
1 H000
HO3 2 0.410<br><br>Therefore, when I specify the B state, I do it like this (i.e. with out<br>specifying the change in mass, assuming the GROMACS will also get the B<br>state mass properly from the .atp file):<br>
<br>[ atoms ]<br>;
nr type
resnr residu atom cgnr charge<br>
6
H
1 H000
HO3 2
0.410 CH3 0.180<br><br>Should that work OK? It wasn't obvious to me where to look in the output to<br>make sure it is. I assume I'll be able to answer that question myself when<br>3.3.1 is released and I can use "gmxcheck -ab"
<br><br>More details:<br>v. 3.2.1, single precision<br>2 fs timestep<br>All bond lengths constrained (LINCS)<br>Nosé-Hoover thermostat with the polymer and solvent coupled separately; 0.4<br>ps time constant<br>Berendsen barostat with
0.5 ps time constant and 4.5x10-5 bar -1<br>compressibility<br>Triple-range cut-off method<br><br>Thanks,<br>Jonathan<br><br>____________________________<br>Jonathan Moore, Ph.D .<br>Research and Engineering Sciences - New Products
<br>Core R&D<br>The Dow Chemical Company<br>1702 Building, Office 4E<br>Midland, MI 48674 USA<br>Phone: (989) 636-9765<br>Fax: (989) 636-4019<br>E Mail: <a href="mailto:jmoore2@dow.com" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">
jmoore2@dow.com</a><br><br>_______________________________________________
<br>gmx-users mailing list<br><a href="mailto:gmx-users@gromacs.org" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">gmx-users@gromacs.org</a><br><a href="http://www.gromacs.org/mailman/listinfo/gmx-users" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">
http://www.gromacs.org/mailman/listinfo/gmx-users</a><br>Please don't post (un)subscribe requests to the list. Use the
<br>www interface or send it to <a href="mailto:gmx-users-request@gromacs.org" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">gmx-users-request@gromacs.org</a> .<br>_______________________________________________
<br>gmx-users mailing list<br><a href="mailto:gmx-users@gromacs.org" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">
gmx-users@gromacs.org</a><br><a href="http://www.gromacs.org/mailman/listinfo/gmx-users" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">http://www.gromacs.org/mailman/listinfo/gmx-users</a><br>Please don't post (un)subscribe requests to the list. Use the
<br>www interface or send it to
<a href="mailto:gmx-users-request@gromacs.org" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">gmx-users-request@gromacs.org</a>.<br></blockquote></div><br>