Jonathan,<br>
<br>
To clarify slightly, the singularity for lambda->1 is when the
lambda=1 Hamiltonian corresponds to a completely disappeared particle.
If you are changing on vdW radii into another vdW radii things are
somewhat more complicated and it isn't clear to me that you really have
the same sort of problem. I wouldn't be surprised if you did, though,
but it would be at some intermediate lambda value. I would still try
the soft core option, especially if your time series of dV/dlambda have
anything that looks at all like shot noise.<br>
<br>
David<br>
<br>
<br><div><span class="gmail_quote">On 11/11/05, <b class="gmail_sendername">David Mobley</b> <<a href="mailto:dmobley@gmail.com">dmobley@gmail.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;">
Jonathan,<br>
<br>
Two comments:<br>
(1) Although this wasn't your question, you might not want to run the
simulations sequentially. I think previous work has shown that you can
introduce some hysteresis this way, although the reference eludes me at
the moment. <br>
(2) When changing vdW parameters, soft core helps a lot. From what
you're describing it sounds like you are doing that. If you don't use
soft core, dV/dlambda actually has a singularity for particle insertion
as lambda->1 (assuming the exponent of lambda is 1, which it is by
default) and is also very susceptible to shot-type noise, and
susceptible to crashes when forces get too large (which often happens
near lambda=1). You can reduce the crashes somewhat by using smaller
timesteps but this doesn't help with the shot noise or singularity
problems. To understand why this happens you can look up van
Gunsteren's paper from around 1992 introducing the soft core potential,
or contact me directly for a brief explanation. If you have any trouble
finding the reference I can look it up for you, but I don't have it in
front of me right this instant. <br><span class="sg">
<br>
David Mobley<br>
UCSF</span><div><span class="e" id="q_10780ff3f6c488fa_2"><br>
<br>
<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;">
<br>I have a question about free energy calculations.<br><br>I'm
performing several different simulations where the hydrogen of certain
hydroxyls on my molecule is replaced with a methyl
group. It's a united atom model, so the change from hydrogen
to methyl doesn't change the number of sites. I have a
single molecule of interest in SPC water (see more details
below). I'm experiencing the problem that on occasion the
system blows up. I've had no problems like this previously
when simulating only one state or the other (methylated or
unmethylated), so I think the problem is due to the free energy
calculation. Unless I've made an error in setting up the
free energy calculation, I assume my problem is that the system isn't
stable for the timestep that I'm using (2 fs) and how fast I'm changing
lambda (I'm incrementing lambda by 0.05 followed by 200 ps of
simulation before increasing it again).<br><br>Is it common to either
decrease the time step or use lambda increments smaller than 0.05 when
doing free energy calculations this way? I'm planning to try
both to see if they fix the problem, but I'm curious if anyone else has
had a problem like this. I'm not using soft cores, but maybe
I have to?<br><br>Also, I'm specifying the atom types like this (i.e., with the mass being 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 specifying
the change in mass, assuming the GROMACS will also get the B 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 make sure it is. I assume I'll be able to answer that
question myself when 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 ps time constant<br>Berendsen barostat with 0.5 ps time constant and 4.5x10-5 bar -1 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></blockquote></div><br>
</span></div></blockquote></div><br>