Re: Problems with temperature control.

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

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 >= $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 <
jobamo1_at_etsid.upv.es> wrote:

> Hi Brian, sorry for my stupidity and thank you very much for being so
> helpful. This is how I've implemented the rescalevels, the same error of
> unestable system appears. I don't understand why...
>
> Is it correct? Or is it wrong?
>
> Thank you very much.
>
>
> binvelocities testD2.vel
>
> run 100
>
> for { set TEMP 600 } { $TEMP >= 300 } { incr TEMP -1 } {
>
> reassignTemp $TEMP
>
> rescalevels sqrt($TEMP+300)
>
> run 10
> }
>
> reassignFreq 40
>
>
> # Constant Temperature Control
>
> langevin off ;# do langevin dynamics
>
>
>
>
>
> Quoting Brian Radak <brian.radak_at_gmail.com>:
>
> You should only need to adjust langevinTemp in order to accomplish
>> temperature control. I might also recommend rescaling the velocities
>> (rescalevels <scale factor>) at each interval by the sqrt ratio of the old
>> and new temperatures as is done for simulated tempering, but this is
>> probably overkill for a simple annealing run.
>>
>> Again, pressure control is not your friend when the temperature changes
>> quickly. The exclusion errors you are seeing almost always come up when
>> the
>> simulation cell changes rapidly for any reason.
>>
>> On Tue, Jul 31, 2018, 7:37 AM Josep Ivan Balaguer Molins <
>> jobamo1_at_etsid.upv.es> wrote:
>>
>> Hell again Brian, I think the problem is using the reassign
>>> temperature control. I would like to try doing the same with Langvine
>>> if posible:
>>>
>>> for { set TEMP 600 } { $TEMP >= 300 } { incr TEMP -10 } {
>>>
>>> langevin on ;# do langevin dynamics
>>> langevinDamping 2 ;# damping coefficient (gamma) of 1/ps
>>> langevinTemp TEMP
>>> langevinHydrogen no ;# don't couple langevin bath to hydrogens
>>> langevinPiston on
>>> langevinPistonTarget 1.01325 ;# in bar -> 1 atm
>>> langevinPistonPeriod 200
>>> langevinPistonDecay 150
>>> langevinPistonTemp TEMP
>>> }
>>>
>>> This doesn't work, how should I code it? I don't see any example on
>>> internet nor in the manual.
>>>
>>> Best regards,
>>>
>>> Ivan
>>> Quoting Josep Ivan Balaguer Molins <jobamo1_at_etsid.upv.es>:
>>>
>>> > Thank you very much Brian. I've turned off the pressure control. And
>>> > I've changed the temperature increment from -10 to -0.1 and even to
>>> > -0.001:
>>> >
>>> > run 1000
>>> >
>>> > for { set TEMP 600 } { $TEMP >= 300 } { incr TEMP -0.1 } {
>>> > reassignTemp $TEMP
>>> >
>>> > run 1000
>>> > }
>>> >
>>> > reassignFreq 100
>>> >
>>> > In any case I get the warning "low global exclusion count" as you
>>> > expected before. And the simulation stops:
>>> >
>>> > ERROR: Constraint failure in RATTLE algorithm for atom 1347!
>>> > ERROR: Constraint failure; simulation has become unstable.
>>> > ERROR: Constraint failure in RATTLE algorithm for atom 9897!
>>> > ERROR: Constraint failure; simulation has become unstable.
>>> > ERROR: Constraint failure in RATTLE algorithm for atom 12147!
>>> > ERROR: Constraint failure; simulation has become unstable.
>>> > ERROR: Constraint failure in RATTLE algorithm for atom 8997!
>>> > ERROR: Constraint failure; simulation has become unstable.
>>> >
>>> > I've read this:
>>> > http://www.ks.uiuc.edu/Research/namd/wiki/index.cgi?NamdTrou
>>> bleshooting
>>> >
>>> > I've increased the margin from 4 to even 20 without succes. The same
>>> > atoms become unestable. I've decreased the reassingFreq but it
>>> > doesn't solve the problem either.
>>> >
>>> > Any idea?
>>> >
>>> > Thank you very much.
>>> >
>>> > Ivan
>>> >
>>> >
>>> >
>>> >
>>> > Quoting Brian Radak <brian.radak_at_gmail.com>:
>>> >
>>> >> Sorry, I meant reassignFreq. In any event, that is not what is
>>> apparently
>>> >> causing the error.
>>> >>
>>> >> I'm actually not sure why you would get a low global bond count (I
>>> might
>>> >> have expected a low global exclusion count). My guess is the system is
>>> >> cooling too quickly and thus pressure control is causing issues. In
>>> >> retrospect, the pressure control is probably not a good idea in this
>>> case,
>>> >> since your temperature change is very large.
>>> >>
>>> >> On Mon, Jul 30, 2018 at 10:41 AM, Josep Ivan Balaguer Molins <
>>> >> jobamo1_at_etsid.upv.es> wrote:
>>> >>
>>> >>> Sorry, Brian but I don't understand... Here is the code were I get
>>> the
>>> >>> error.
>>> >>>
>>> >>>
>>> >>> #############################################################
>>> >>> ## ADJUSTABLE PARAMETERS ##
>>> >>> #############################################################
>>> >>> gromacs on
>>> >>> grotopfile 32x50LA_crystal.top
>>> >>> grocoorfile 32x50LA_crystal.gro
>>> >>> bincoordinates testD2.coor
>>> >>>
>>> >>> paraTypeCharmm on
>>> >>> exclude scaled1-4
>>> >>> 1-4scaling 0.5
>>> >>>
>>> >>> switching on
>>> >>> switchdist 8 # at 8Å we start to smooth electrostatic to 0
>>> >>> cutoff 12.0
>>> >>> pairlistdist 26.0
>>> >>> margin 4
>>> >>>
>>> >>> # Integrator Parameters
>>> >>>
>>> >>>
>>> >>>
>>> >>> rigidBonds all ;# needed for 2fs steps
>>> >>> nonbondedFreq 1
>>> >>> vdwGeometricSigma yes
>>> >>> fullElectFrequency 2
>>> >>> stepspercycle 20
>>> >>> pairlistsperCycle 2
>>> >>>
>>> >>> wrapAll on
>>> >>>
>>> >>> timestep 1.0 ;# 2fs/step 2 en simulación, 1 minimizacion
>>> >>>
>>> >>> outputEnergies 100
>>> >>> outputtiming 100
>>> >>> outputPressure 100
>>> >>> binaryoutput yes
>>> >>> outputname testD3
>>> >>> dcdfreq 200
>>> >>> restartfreq 1000
>>> >>> restartname rest_TestD3
>>> >>>
>>> >>> binvelocities testD2.vel
>>> >>>
>>> >>> reassignFreq 100
>>> >>> run 80
>>> >>> for { set TEMP 550 } { $TEMP >= 300 } { incr TEMP -10 } {
>>> >>> reassignTemp $TEMP
>>> >>> run 100
>>> >>> }
>>> >>>
>>> >>>
>>> >>>
>>> >>>
>>> >>> # Constant Temperature Control
>>> >>>
>>> >>> langevin off ;# do langevin dynamics
>>> >>> langevinDamping 2 ;# damping coefficient (gamma) of 1/ps
>>> >>> langevinTemp 300
>>> >>> langevinHydrogen no ;# don't couple langevin bath to
>>> hydrogens
>>> >>> langevinPiston on
>>> >>> langevinPistonTarget 1.01325 ;# in bar -> 1 atm
>>> >>> langevinPistonPeriod 10000
>>> >>> langevinPistonDecay 15000
>>> >>> langevinPistonTemp 300
>>> >>>
>>> >>> useGroupPressure yes ;# needed for rigidBonds
>>> >>> useFlexibleCell no
>>> >>> useConstantArea no
>>> >>>
>>> >>> extendedSystem testD2.xsc
>>> >>>
>>> >>> run 300000
>>> >>>
>>> >>>
>>> >>> I don't have a reassignIncr here... or do you mean I have to make it
>>> >>> smaller on the previous minimization and simulation?
>>> >>> If I turn langevin off, tha same error appears.
>>> >>>
>>> >>> Thank you very much for your answer.
>>> >>>
>>> >>> Best regards
>>> >>>
>>> >>> Ivan
>>> >>>
>>> >>>
>>> >>> Quoting Brian Radak <brian.radak_at_gmail.com>:
>>> >>>
>>> >>> Note that the Langevin and reassign thermostats are completely
>>> decoupled
>>> >>>> code paths. My guess is that the reassignIncr is too large compared
>>> to the
>>> >>>> relaxation rate dictated by Langevin damping and the two are
>>> redundant -
>>> >>>> personally I would only use one or the other (you can also script
>>> changes
>>> >>>> to langevinTemp). You could also try using a smaller damping
>>> coefficient,
>>> >>>> but this really isn't any different from not using Langevin at all--000000000000a0842505724f5056--

This archive was generated by hypermail 2.1.6 : Wed Sep 18 2019 - 23:20:12 CDT