Re: SMD and CPT? (fwd)

From: Marcos Sotomayor (
Date: Fri Feb 09 2007 - 15:23:38 CST

Hi Sterling,

Let me be more specific. First, when the external force applied to your
system is time-dependent, the Hamiltonian of your system is not conserved anymore
and does not represent the total energy. Moreover, the expression F
\times d, as far as I can tell, is no longer the work done on the system
(usually one assumes a time independent hamiltonian/force to derive that
expression). So my sentence should have read like: " The quantity F
\times d (previously work done) should only be similar to the energy
change". Please correct me if I am wrong on that.

I think you are right in that the total amount of work done by
the external force does match the energy change, but you cannot use Fd to
measure it. Therefore, the first law of thermodynamics is not violated.
There is a nice article on a similar topic that you may want to check:

"Invariants for time-dependent Hamiltonian systems"
Jurgen Struckmeier and Claus Riedel
Physical Review E 64 026503


On Fri, 9 Feb 2007, Sterling Paramore wrote:

> Marcos,
> You've mentioned that the work done on the system should only be similar to
> the energy change. If the simulations are done at NVE, then they should be
> identical according to the First Law (to within the precision of the MD
> program and your integration method). Under NVE, there's no other place for
> the work to go except for into a change in the total energy, regardless of
> whether the force is time-dependent or not.
> -Sterling
> Gianluca Interlandi wrote:
>> I'm sending Marcos' answer to my last e-mail to the whole list as many
>> people might be interested in our discussion.
>> Marcos, Thank you very much for your help.
>> Best,
>> Gianluca
>> ---------- Forwarded message ----------
>> Date: Fri, 9 Feb 2007 13:43:32 -0600 (CST)
>> From: Marcos Sotomayor <>
>> To: Gianluca Interlandi <>
>> Subject: Re: namd-l: SMD and CPT?
>> Dear Gianluca,
>> Two comments on your plots:
>> - I see you are using a 2fs time step. Are you using rigid bonds and/or
>> multiple
>> time stepping? (I guess you are using rigid bonds and not multiple time
>> stepping
>> since the temperature seems to be quite under control, good!)
>> - Check the energy of your system (TOTAL3 in namd output) and compute the
>> difference between the starting value and the final value (due to intrinsic
>> fluctuations of the energy you will be off by ~ +/- 10 kcal/mol). Then plot
>> Force vs distance (I assume the plot you sent is F vs t). The area under
>> Force
>> vs distance (you can integrate F vs d numerically using xmgrace, be sure to
>> use
>> the appropriate units) should be similar (not identical since you are
>> applying a
>> time dependent force) to the energy difference you computed before. That's
>> a
>> good way to check that the simulation is OK in terms of the physics behind
>> it.
>> The temperature rises at the breaking point because that's when the
>> external
>> force starts to do work (F \times d). If you plot the energy of your
>> system,
>> likely you will see that increases at the breaking point as well.
>> Although I feel more comfortable with NVE simulations, do not completely
>> disregard your CPT simulation, stiffer springs will give you a lot of large
>> fluctuations! Likely that's the reason why you see forces on the order of
>> 1200pN
>> and probably you will see something similar in your NVE simulation. Plot F
>> vs t
>> and a running average (just like you did for the temperature) to see at
>> what
>> level is the force.
>> Finally, the choice of spring constant is a tough one. AFM people use
>> really
>> soft cantilevers, while if you want to get a potential of mean force (PMF)
>> you
>> need to use a stiff spring (have a look at Park and Schulten, Journal of
>> Chemical Physics, 120:5946-5961, 2004.)
>> Hope this helps and don't hesitate to contact me again if you have further
>> questions.
>> Regards,
>> Marcos
>> PS: perhaps you could resend this e-mail (without the attachment) to
>> namd-l,
>> since I think many people would benefit from our discussion even without
>> the
>> plot.
>> On Fri, 9 Feb 2007, Gianluca Interlandi wrote:
>>> Dear Marcos,
>>> Thanks a lot for your answer and for offering me your help.
>>> I have already run a bunch of pulling simulations in the NVE ensemble. I
>>> attach a plot of a constant velocity simulation where I'm pulling two
>>> proteins apart. At the point of breaking the temperature raises of a few
>>> degrees. The force peak is at around 670 pN. Recently, I have started a
>>> new simulation of the same system but with a stiffer spring, i.e., 140
>>> pN/A instead of 14 pN/A, and in the CPT (Langevin as thermostat and
>>> Langevin-Nose Hoover as barostat). After 5 ns the force has reached 1200
>>> pN but the two proteins are still bound and the total RMSD is still around
>>> 2 A. I started wondering whether this is due to the stiffer spring only or
>>> whether CPT has introduced some artifacts. To investigate this I have
>>> started a new simulation with exactly the same parameters but in the NVE.
>>> The latter simulation has just started so I'm waiting for the results.
>>> In any case, now I know that I have to stick with NVE.
>>> Many thanks,
>>> Gianluca
>>> On Fri, 9 Feb 2007, Marcos Sotomayor wrote:
>>>> Actually, I have three reasons to disagree with Sterling.
>>>> 1 - Performing simulations in the NVE ensemble permits you to monitor the
>>>> change in energy of your system, and check that it is at least similar to
>>>> the
>>>> work done by the external applied force (and exactly the same when using
>>>> constant forces). I've monitored temperature changes of 0.5 K in
>>>> simulations
>>>> that do not use multiple time stepping, which matches well the work done
>>>> on
>>>> the system by external forces.
>>>> 2 - Algorithms used to control temperature usually are not designed to
>>>> handle
>>>> systems in which external forces induce motions in a preferred direction.
>>>> The
>>>> Langevin thermostat for instance, will apply a net (average) zero force
>>>> to
>>>> atoms that do not move on a preferred direction, but likely will
>>>> introduce
>>>> an
>>>> artificial viscous drag to atoms that are being pulled. Thus, your force
>>>> peaks
>>>> will be larger.
>>>> 3 - Some temperature control methods induced center of mass motion, which
>>>> is
>>>> not good when you have fixed reference points like in SMD.
>>>> The best solution is to perform simulations in both ensembles and
>>>> compare,
>>>> focusing on those results that are "ensemble independent".
>>>> However, with a water box big enough I suggest to perform SMD simulations
>>>> in
>>>> the NVE ensemble.
>>>> Ginaluca, that's the reason why we performed simulations in the NVE
>>>> ensemble
>>>> in our Structure paper. Let me know if you have more questions, I'll be
>>>> glad
>>>> to discuss them with you.
>>>> Regards,
>>>> Marcos
>>>> On Fri, 9 Feb 2007, Sterling Paramore wrote:
>>>>> I would definitely suggest you do SMD with a thermostat. Otherwise, the
>>>>> viscous heating that occurs when you pull on your system will increase
>>>>> the
>>>>> temperature. In any real system, the heat produced would be transferred
>>>>> to
>>>>> the surrounding thermal reservoir. The artificial thermostat is just a
>>>>> way
>>>>> to model this effect.
>>>>> -Sterling
>>>>> On Feb 8, 2007, at 11:49 PM, Gianluca Interlandi wrote:
>>>>>> I have a question concerning steered molecular dynamics simulations
>>>>>> (constant force and constant velocity). Is it appropriate to use a
>>>>>> thermostat and barostat (CPT) while performing a constant force or
>>>>>> constant velocity pulling simulation? I have seen that many people
>>>>>> prefer
>>>>>> NVE, i.e., no thermostat.
>>>>>> I expect my protein to undergo large conformational changes during
>>>>>> pulling. Does a thermostat slow down the sequence of events, since part
>>>>>> of
>>>>>> the applied force is converted into heat?
>>>>>> Many thanks,
>>>>>> Gianluca
>>>>>> -----------------------------------------------------
>>>>>> Dr. Gianluca Interlandi
>>>>>> +1 (206) 685 4435
>>>>>> +1 (206) 714 4303
>>>>>> Postdoc at the Department of Bioengineering
>>>>>> at the University of Washington, Seattle WA U.S.A.
>>>>>> -----------------------------------------------------
>>> -----------------------------------------------------
>>> Dr. Gianluca Interlandi
>>> +1 (206) 685 4435
>>> +1 (206) 714 4303
>>> Postdoc at the Department of Bioengineering
>>> at the University of Washington, Seattle WA U.S.A.
>>> -----------------------------------------------------

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:44:23 CST