Re: Problems with temperature control.

From: Josep Ivan Balaguer Molins (jobamo1_at_etsid.upv.es)
Date: Wed Aug 01 2018 - 12:34:30 CDT

Thank you very much for your help, I really appreciate it. This is how
I've copied pasted the code:

binvelocities testD2.vel

set rampInterval 20000
set T1 600
set TN 300
set dT -0.1

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
}

useGroupPressure yes # needed for rigidBonds
useFlexibleCell no
useConstantArea no

extendedSystem testD2.xsc

run 300000

As you can see I didin't used 'langevin on'. I don't mind cooling
slowly. The error is the same but with other atoms:

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.
ERROR: Exiting prematurely; see error messages above

I'm doing this to get the glass transition of a polymeric chain. For
this, I've run a minimization from 200 to 600K without problems, and
after that, I've done an isothermal simulation of 400000 integration
steps at 600K using 'lagenvin on' without any problem and using the
outputfiles from the minimization. Now, using the outputfiles of the
600K isothermal simulation, I have to cool it down to 300K but I need
to control the speed rate of this cooling to see the glass transition.
If I set langevin to 300K using the 600K simulation output files, the
temperature goes down in less than 2000 integration steps, and with
rescale I get a cool down by steps, the temperature output un vmd
looks like a stair. But when I try using the tcl scripting to make a
y=-mx+n type ramp I get the error.

I don't know what to do next...

Best regards,

Ivan

Quoting Brian Radak <brian.radak_at_gmail.com>:

> 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

This archive was generated by hypermail 2.1.6 : Mon Nov 18 2019 - 23:20:02 CST