Re: Problems with temperature control.

From: Brian Radak (brian.radak_at_gmail.com)
Date: Tue Jul 31 2018 - 14:57:09 CDT

The thermostat response is probably a bit faster than most people expect,
especially with damping coefficients on the order of ps^-1. Try allocating
initial velocities at 300 K but then using langevinTemp 400, you'll
probably see few, if any, frames where the instantaneous temperature is
close to 300. The difference would be even more dramatic if you start at
0-100 K (which is just silly IMO).

However, just because the velocities are from the "correct" distribution
that doesn't mean the coordinates make any sense at all. This is why
rapidly heating a system (over the course of 10s of ps, say) doesn't really
accomplish much in the way of equilibration. Even at high temperatures you
need quite a bit of time for the system to respond to a temperature change.

On Tue, Jul 31, 2018 at 3:18 PM, Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov>
wrote:

> Just an interested observer here... If you change langevinTemp alone, why
> would there be a kinetic energy jump? Your random forces would start
> getting larger/smaller, and the damping would bring the system up to the
> new temperature.
>
> I think the history of where these ramping profiles came from is
> informative. CHARMM and AMBER used to recommend velocity quenching
> minimization when their minimizers weren't all that great. The temperature
> ramps were just a way of letting the system settle itself down before
> production, as jumping directly to large temperature could blow up your
> system. There isn't much point to temperature ramping nowadays, except
> perhaps to seed temperature replica exchange simulations.
>
> -Josh
>
>
>
> On 2018-07-31 12:02:12-06:00 owner-namd-l_at_ks.uiuc.edu wrote:
>
> This isn't stupidity, just a really tricky usage of the code. I think
> you're making this harder than it needs to be - aren't you just trying to
> enhance relaxation to a thermodynamically reasonable state?
>
> The "reassign" keywords are very old and kind of unnecessary IMO. Were I to
> implement a heating protocol I might do it as follows (this is completely
> untested and I can't guarantee that it will work, let alone apply to your
> given system/problem):
>
> set rampInterval 5000
> set T1 600
> set TN 300
> set dT -10
>
> reinitvels $T1
> set oldTemp $T1
> for {set temp $T1} {$temp &gt;= $TN} {incr temp $dT} {
> langevinTemp $temp
> rescalevels [expr {sqrt($temp/$oldTemp)}] ;# this effectively "pre-heats"
> the system so that you don't get sudden changes in kinetic energy
> run $rampInterval
> set oldTemp $temp
> }
>
> Note that "rampInterval" and "dT" are fairly critical values in order for
> this to work. The former needs to be suitably long so that the system can
> actually respond to the thermostat, while the latter needs to be suitably
> small so that heating is not too strong, but also large enough that you
> don't waste all of your time at intermediate temperatures.
>
> HTH,
> BKR
>
>
>
> On Tue, Jul 31, 2018 at 11:54 AM, Josep Ivan Balaguer Molins &lt;jobamo1_at_etsid.upv.es&gt; wrote:
>
> &gt; Hi Brian, sorry for my stupidity and thank you very much for being so
> &gt; helpful. This is how I've implemented the rescalevels, the same error of
> &gt; unestable system appears. I don't understand why...
> &gt;
> &gt; Is it correct? Or is it wrong?
> &gt;
> &gt; Thank you very much.
> &gt;
> &gt;
> &gt; binvelocities testD2.vel
> &gt;
> &gt; run 100
> &gt;
> &gt; for { set TEMP 600 } { $TEMP &gt;= 300 } { incr TEMP -1 } {
> &gt;
> &gt; reassignTemp $TEMP
> &gt;
> &gt; rescalevels sqrt($TEMP+300)
> &gt;
> &gt; run 10
> &gt; }
> &gt;
> &gt; reassignFreq 40
> &gt;
> &gt;
> &gt; # Constant Temperature Control
> &gt;
> &gt; langevin off ;# do langevin dynamics
> &gt;
> &gt;
> &gt;
> &gt;
> &gt;
> &gt; Quoting Brian Radak <brian.radak_at_gmail.com>:
> &gt;
> &gt; You should only need to adjust langevinTemp in order to accomplish
> &gt;&gt; temperature control. I might also recommend rescaling the velocities
> &gt;&gt; (rescalevels <scale factor="">) at each interval by the sqrt ratio of the old
> &gt;&gt; and new temperatures as is done for simulated tempering, but this is
> &gt;&gt; probably overkill for a simple annealing run.
> &gt;&gt;
> &gt;&gt; Again, pressure control is not your friend when the temperature changes
> &gt;&gt; quickly. The exclusion errors you are seeing almost always come up when
> &gt;&gt; the
> &gt;&gt; simulation cell changes rapidly for any reason.
> &gt;&gt;
> &gt;&gt; On Tue, Jul 31, 2018, 7:37 AM Josep Ivan Balaguer Molins &lt;
> &gt;&gt; jobamo1_at_etsid.upv.es&gt; wrote:
> &gt;&gt;
> &gt;&gt; Hell again Brian, I think the problem is using the reassign
> &gt;&gt;&gt; temperature control. I would like to try doing the same with Langvine
> &gt;&gt;&gt; if posible:
> &gt;&gt;&gt;
> &gt;&gt;&gt; for { set TEMP 600 } { $TEMP &gt;= 300 } { incr TEMP -10 } {
> &gt;&gt;&gt;
> &gt;&gt;&gt; langevin on ;# do langevin dynamics
> &gt;&gt;&gt; langevinDamping 2 ;# damping coefficient (gamma) of 1/ps
> &gt;&gt;&gt; langevinTemp TEMP
> &gt;&gt;&gt; langevinHydrogen no ;# don't couple langevin bath to hydrogens
> &gt;&gt;&gt; langevinPiston on
> &gt;&gt;&gt; langevinPistonTarget 1.01325 ;# in bar -&gt; 1 atm
> &gt;&gt;&gt; langevinPistonPeriod 200
> &gt;&gt;&gt; langevinPistonDecay 150
> &gt;&gt;&gt; langevinPistonTemp TEMP
> &gt;&gt;&gt; }
> &gt;&gt;&gt;
> &gt;&gt;&gt; This doesn't work, how should I code it? I don't see any example on
> &gt;&gt;&gt; internet nor in the manual.
> &gt;&gt;&gt;
> &gt;&gt;&gt; Best regards,
> &gt;&gt;&gt;
> &gt;&gt;&gt; Ivan
> &gt;&gt;&gt; Quoting Josep Ivan Balaguer Molins <jobamo1_at_etsid.upv.es>:
> &gt;&gt;&gt;
> &gt;&gt;&gt; &gt; Thank you very much Brian. I've turned off the pressure control. And
> &gt;&gt;&gt; &gt; I've changed the temperature increment from -10 to -0.1 and even to
> &gt;&gt;&gt; &gt; -0.001:
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; run 1000
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; for { set TEMP 600 } { $TEMP &gt;= 300 } { incr TEMP -0.1 } {
> &gt;&gt;&gt; &gt; reassignTemp $TEMP
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; run 1000
> &gt;&gt;&gt; &gt; }
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; reassignFreq 100
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; In any case I get the warning "low global exclusion count" as you
> &gt;&gt;&gt; &gt; expected before. And the simulation stops:
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; ERROR: Constraint failure in RATTLE algorithm for atom 1347!
> &gt;&gt;&gt; &gt; ERROR: Constraint failure; simulation has become unstable.
> &gt;&gt;&gt; &gt; ERROR: Constraint failure in RATTLE algorithm for atom 9897!
> &gt;&gt;&gt; &gt; ERROR: Constraint failure; simulation has become unstable.
> &gt;&gt;&gt; &gt; ERROR: Constraint failure in RATTLE algorithm for atom 12147!
> &gt;&gt;&gt; &gt; ERROR: Constraint failure; simulation has become unstable.
> &gt;&gt;&gt; &gt; ERROR: Constraint failure in RATTLE algorithm for atom 8997!
> &gt;&gt;&gt; &gt; ERROR: Constraint failure; simulation has become unstable.
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; I've read this:
> &gt;&gt;&gt; &gt; https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.ks.uiuc.edu%2FResearch%2Fnamd%2Fwiki%2Findex.cgi%3FNamdTrou&amp;data=02%7C01%7CJoshua.Vermaas%40nrel.gov%7Cd6ed5e5a2cac4ad70bba08d5f70fba3c%7Ca0f29d7e28cd4f5484427885aee7c080%7C0%7C0%7C636686569307280344&amp;sdata=VyV0T8wz2jeqPpWkcX4QbOYlQtaDvCZQAF0NndHKDc4%3D&amp;reserved=0
> &gt;&gt;&gt <https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.ks.uiuc.edu%2FResearch%2Fnamd%2Fwiki%2Findex.cgi%3FNamdTrou&amp;data=02%7C01%7CJoshua.Vermaas%40nrel.gov%7Cd6ed5e5a2cac4ad70bba08d5f70fba3c%7Ca0f29d7e28cd4f5484427885aee7c080%7C0%7C0%7C636686569307280344&amp;sdata=VyV0T8wz2jeqPpWkcX4QbOYlQtaDvCZQAF0NndHKDc4%3D&amp;reserved=0&gt;&gt;&gt>; bleshooting
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; I've increased the margin from 4 to even 20 without succes. The same
> &gt;&gt;&gt; &gt; atoms become unestable. I've decreased the reassingFreq but it
> &gt;&gt;&gt; &gt; doesn't solve the problem either.
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; Any idea?
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; Thank you very much.
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; Ivan
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt; Quoting Brian Radak <brian.radak_at_gmail.com>:
> &gt;&gt;&gt; &gt;
> &gt;&gt;&gt; &gt;&gt; Sorry, I meant reassignFreq. In any event, that is not what is
> &gt;&gt;&gt; apparently
> &gt;&gt;&gt; &gt;&gt; causing the error.
> &gt;&gt;&gt; &gt;&gt;
> &gt;&gt;&gt; &gt;&gt; I'm actually not sure why you would get a low global bond count (I
> &gt;&gt;&gt; might
> &gt;&gt;&gt; &gt;&gt; have expected a low global exclusion count). My guess is the system is
> &gt;&gt;&gt; &gt;&gt; cooling too quickly and thus pressure control is causing issues. In
> &gt;&gt;&gt; &gt;&gt; retrospect, the pressure control is probably not a good idea in this
> &gt;&gt;&gt; case,
> &gt;&gt;&gt; &gt;&gt; since your temperature change is very large.
> &gt;&gt;&gt; &gt;&gt;
> &gt;&gt;&gt; &gt;&gt; On Mon, Jul 30, 2018 at 10:41 AM, Josep Ivan Balaguer Molins &lt;
> &gt;&gt;&gt; &gt;&gt; jobamo1_at_etsid.upv.es&gt; wrote:
> &gt;&gt;&gt; &gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; Sorry, Brian but I don't understand... Here is the code were I get
> &gt;&gt;&gt; the
> &gt;&gt;&gt; &gt;&gt;&gt; error.
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; #############################################################
> &gt;&gt;&gt; &gt;&gt;&gt; ## ADJUSTABLE PARAMETERS ##
> &gt;&gt;&gt; &gt;&gt;&gt; #############################################################
> &gt;&gt;&gt; &gt;&gt;&gt; gromacs on
> &gt;&gt;&gt; &gt;&gt;&gt; grotopfile 32x50LA_crystal.top
> &gt;&gt;&gt; &gt;&gt;&gt; grocoorfile 32x50LA_crystal.gro
> &gt;&gt;&gt; &gt;&gt;&gt; bincoordinates testD2.coor
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; paraTypeCharmm on
> &gt;&gt;&gt; &gt;&gt;&gt; exclude scaled1-4
> &gt;&gt;&gt; &gt;&gt;&gt; 1-4scaling 0.5
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; switching on
> &gt;&gt;&gt; &gt;&gt;&gt; switchdist 8 # at 8Å we start to smooth electrostatic to 0
> &gt;&gt;&gt; &gt;&gt;&gt; cutoff 12.0
> &gt;&gt;&gt; &gt;&gt;&gt; pairlistdist 26.0
> &gt;&gt;&gt; &gt;&gt;&gt; margin 4
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; # Integrator Parameters
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; rigidBonds all ;# needed for 2fs steps
> &gt;&gt;&gt; &gt;&gt;&gt; nonbondedFreq 1
> &gt;&gt;&gt; &gt;&gt;&gt; vdwGeometricSigma yes
> &gt;&gt;&gt; &gt;&gt;&gt; fullElectFrequency 2
> &gt;&gt;&gt; &gt;&gt;&gt; stepspercycle 20
> &gt;&gt;&gt; &gt;&gt;&gt; pairlistsperCycle 2
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; wrapAll on
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; timestep 1.0 ;# 2fs/step 2 en simulación, 1 minimizacion
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; outputEnergies 100
> &gt;&gt;&gt; &gt;&gt;&gt; outputtiming 100
> &gt;&gt;&gt; &gt;&gt;&gt; outputPressure 100
> &gt;&gt;&gt; &gt;&gt;&gt; binaryoutput yes
> &gt;&gt;&gt; &gt;&gt;&gt; outputname testD3
> &gt;&gt;&gt; &gt;&gt;&gt; dcdfreq 200
> &gt;&gt;&gt; &gt;&gt;&gt; restartfreq 1000
> &gt;&gt;&gt; &gt;&gt;&gt; restartname rest_TestD3
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; binvelocities testD2.vel
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; reassignFreq 100
> &gt;&gt;&gt; &gt;&gt;&gt; run 80
> &gt;&gt;&gt; &gt;&gt;&gt; for { set TEMP 550 } { $TEMP &gt;= 300 } { incr TEMP -10 } {
> &gt;&gt;&gt; &gt;&gt;&gt; reassignTemp $TEMP
> &gt;&gt;&gt; &gt;&gt;&gt; run 100
> &gt;&gt;&gt; &gt;&gt;&gt; }
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; # Constant Temperature Control
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; langevin off ;# do langevin dynamics
> &gt;&gt;&gt; &gt;&gt;&gt; langevinDamping 2 ;# damping coefficient (gamma) of 1/ps
> &gt;&gt;&gt; &gt;&gt;&gt; langevinTemp 300
> &gt;&gt;&gt; &gt;&gt;&gt; langevinHydrogen no ;# don't couple langevin bath to
> &gt;&gt;&gt; hydrogens
> &gt;&gt;&gt; &gt;&gt;&gt; langevinPiston on
> &gt;&gt;&gt; &gt;&gt;&gt; langevinPistonTarget 1.01325 ;# in bar -&gt; 1 atm
> &gt;&gt;&gt; &gt;&gt;&gt; langevinPistonPeriod 10000
> &gt;&gt;&gt; &gt;&gt;&gt; langevinPistonDecay 15000
> &gt;&gt;&gt; &gt;&gt;&gt; langevinPistonTemp 300
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; useGroupPressure yes ;# needed for rigidBonds
> &gt;&gt;&gt; &gt;&gt;&gt; useFlexibleCell no
> &gt;&gt;&gt; &gt;&gt;&gt; useConstantArea no
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; extendedSystem testD2.xsc
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; run 300000
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; I don't have a reassignIncr here... or do you mean I have to make it
> &gt;&gt;&gt; &gt;&gt;&gt; smaller on the previous minimization and simulation?
> &gt;&gt;&gt; &gt;&gt;&gt; If I turn langevin off, tha same error appears.
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; Thank you very much for your answer.
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; Best regards
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; Ivan
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; Quoting Brian Radak <brian.radak_at_gmail.com>:
> &gt;&gt;&gt; &gt;&gt;&gt;
> &gt;&gt;&gt; &gt;&gt;&gt; Note that the Langevin and reassign thermostats are completely
> &gt;&gt;&gt; decoupled
> &gt;&gt;&gt; &gt;&gt;&gt;&gt; code paths. My guess is that the reassignIncr is too large compared
> &gt;&gt;&gt; to the
> &gt;&gt;&gt; &gt;&gt;&gt;&gt; relaxation rate dictated by Langevin damping and the two are
> &gt;&gt;&gt; redundant -
> &gt;&gt;&gt; &gt;&gt;&gt;&gt; personally I would only use one or the other (you can also script
> &gt;&gt;&gt; changes
> &gt;&gt;&gt; &gt;&gt;&gt;&gt; to langevinTemp). You could also try using a smaller damping
> &gt;&gt;&gt; coefficient,
> &gt;&gt;&gt; &gt;&gt;&gt;&gt; but this really isn't any different from not using Langevin at all</brian.radak_at_gmail.com></brian.radak_at_gmail.com></jobamo1_at_etsid.upv.es></scale></brian.radak_at_gmail.com>
>
>

This archive was generated by hypermail 2.1.6 : Tue Sep 17 2019 - 23:20:07 CDT