[gmx-users] Problem with temperature and stochastic dynamics in FEP
Justin A. Lemkul
jalemkul at vt.edu
Wed Aug 3 03:23:33 CEST 2011
Sanku M wrote:
> Hi Justin,
> I first performed two minimization using steep and l-bfgs as the
> method on the solvated system for each lambda. Then I ran a 500ps NVT
> simulation on minimized system ( for each lambda) where essentially same
> .mdp file was used except the following change:
>
> gen-vel = yes
> gen-temp = 280
> gen-seed = -1
>
> So, essentially, I heated the minimized system up from 280 K with ref_t
> 300K using sd integrator. But, here also I found the average temperature
> goes to 303 K( instead of 300 K) .
Unless you did simulated annealing in between, you didn't heat the system, you
suddenly jolted the temperature, which is probably not stable. Maybe the
integrator/Langevin thermostat does not react well to these conditions.
Equilibrate under the conditions you wish to collect data. Your data collection
uses NPT, so running NVT at a different temperature and suddenly changing the
temperature and introducing pressure coupling is not an appropriate procedure.
If nothing else, your data will be unreliable, even if your temperature was correct.
> One thing I don't understand that why using md as integrator gives the
> desired average temperature( when using Nose-Hoover thermostat) while
> the sd integrator does not.
I can only hazard a guess that for some reason Nose-Hoover can arrive at the
correct temperature more quickly, but that seems uncharacteristic for it given
the nature of its fluctuations. Perhaps someone with more knowledge of the
algorithms and code can comment. But until you've established a properly
equilibrated system under appropriate conditions, it's not worth investigating
these differences much further.
-Justin
> Any further suggestion will be appreciated.
>
> Jagannath
> ------------------------------------------------------------------------
> *From:* Justin A. Lemkul <jalemkul at vt.edu>
> *To:* Discussion list for GROMACS users <gmx-users at gromacs.org>
> *Sent:* Tue, August 2, 2011 8:08:20 PM
> *Subject:* Re: [gmx-users] Problem with temperature and stochastic
> dynamics in FEP
>
>
>
> Sanku M wrote:
> > Hi,
> > I am trying to calculate the solvation free energy using
> thermodynamic integration(TI) method . I am using gromacs 4.0.7 . But, I
> am having a problem in getting accurate average temperature .
> > The following is what I did:
> > I found that when doing TI, grompp recommends using 'sd' (
> stochastic dynamics ) as integrator ( provides an warning that for
> decoupled system, it is better to use sd in stead of md ).
> > Also I went through Justin Lemkul's tutorial which also recommends
> using 'sd' for TI method. So, I tried to use sd integrator.
>
> The tutorial uses BAR for calculation free energy, not TI. Though the
> protocols are similar, there is a difference, FYI.
>
> > I wanted to run NVT simulation at 300 K.
> > I also found from manual and also from Justin's website on FEP, sd
> integrator implicitly controls temperature and so there is no need to
> specify a thermostat. So, I did not specify any thermostat ( I wrote
> tcoupl = No ). But, unfortunately, I found that after a long simulation
> , the average temperature actually goes to 303 K in stead of 300 K( the
> desired temperature). This happened for all Lambda values where the
> average temperature turned out to be 303 K. The .mdp file is pasted at
> the end of the email and I felt I am using reasonable cutoffs and PME as
> electrostatics.
> >
> > Here is the output from g_energy command to get the temperature.
> Statistics over 4000001 steps [ 0.0000 thru 8000.0005 ps ], 1 data sets
> > All averages are exact over 4000001 steps
> >
> > Energy Average RMSD Fluct. Drift
> Tot-Drift
> >
> -------------------------------------------------------------------------------
> > Temperature 303.333 4.78084 4.78082 5.46937e-06
> 0.043755
> > Heat Capacity Cv: 12.4764 J/mol K (factor = 0.00024841)
> >
> >
> > I tried three further tests:
> > a) I removed tcoupl = No option . But, it is giving same 303K as
> average temperature when using sd.
> >
> > a) I tried to specify the Nose-Hoover thermostat along with sd. But
> still, I found that when using sd, the average temperature is going to
> 303 K in stead of 300 K.
> >
>
> From the manual description of the sd integrator: "The parameter tcoupl
> is ignored." This explains (a) and (b).
>
> > c) Finally, I tried to overlook grompp warning and went ahead and
> used 'md' along with Nose-Hoover thermostat. This time I found the right
> average temperature 300 K is being achieved.B
> > But, I understand that for a decoupled system like here, I need to
> use a method like stochastic dynamics.
> > But , I was wondering why it is reaching 303 K in stead of 300 K
> when using sd but 300K is achieved when using md. I looked at the
> mailing list where people had issues with sd and temperature control.
> But, I could not find a good solution. So, any help on the right
> protocol will be really appreciated.
> >
>
> What type of equilibration did you do prior to the data collection? If
> your system isn't sampling the desired ensemble, then you shouldn't proceed.
>
> -Justin
>
> >
> > Here is my .mdp file.
> >
> > ; RUN CONTROL PARAMETERS
> > integrator = sd
> > ; Start time and timestep in ps
> > tinit = 0.0
> > dt = 0.002
> > nsteps = 4000000
> > ; For exact run continuation or redoing part of a run
> > ; Part index is updated automatically on checkpointing (keeps files
> separate)
> > simulation_part = 1
> > init_step = 0
> > ; mode for center of mass motion removal
> > comm-mode = Linear
> > ; number of steps for center of mass motion removal
> > nstcomm = 1
> > ; group(s) for center of mass motion removal
> > comm-grps =
> > ; LANGEVIN DYNAMICS OPTIONS
> > ; Friction coefficient (amu/ps) and random seed
> > bd-fric = 0
> > ld-seed = 1993
> >
> >
> >
> > ; OUTPUT CONTROL OPTIONS
> > ; OUTPUT CONTROL OPTIONS
> > ; Output frequency for coords (x), velocities (v) and forces (f)
> > nstxout = 10000
> > nstvout = 10000
> > nstfout = 0
> > ; Output frequency for energies to log file and energy file
> > nstlog = 1000
> > nstenergy = 10
> > ; Output frequency and precision for xtc file
> > nstxtcout = 1000
> > xtc-precision = 1000
> > ; This selects the subset of atoms for the xtc file. You can
> > ; select multiple groups. By default all atoms will be written.
> > xtc-grps =
> > ; Selection of energy groups
> > energygrps =
> >
> > ; NEIGHBORSEARCHING PARAMETERS
> > ; nblist update frequency
> > nstlist = 10
> > ; ns algorithm (simple or grid)
> > ns_type = grid
> > ; Periodic boundary conditions: xyz, no, xy
> > pbc = xyz
> > periodic_molecules = no
> > ; nblist cut-off
> > rlist = 1.4
> >
> > ; OPTIONS FOR ELECTROSTATICS AND VDW
> > ; Method for doing electrostatics
> > coulombtype = pme
> > rcoulomb-switch = 0
> > rcoulomb = 1.4
> > ; Relative dielectric constant for the medium and the reaction field
> > epsilon_r = 1
> > epsilon_rf = 1
> > ; Method for doing Van der Waals
> > epsilon_rf = 1
> > ; Method for doing Van der Waals
> > vdw-type = shift
> > ; cut-off lengths
> > rvdw-switch = 1.0
> > rvdw = 1.4
> > ; Apply long range dispersion corrections for Energy and Pressure
> > DispCorr = EnerPres
> > ; Extension of the potential lookup tables beyond the cut-off
> > table-extension = 1
> > ; Seperate tables between energy group pairs
> > energygrp_table =
> > ; Spacing for the PME/PPPM FFT grid
> > fourierspacing = 0.12
> > ; FFT grid size, when a value is 0 fourierspacing will be used
> > fourier_nx = 0
> > fourier_ny = 0
> > fourier_nz = 0
> > ; EWALD/PME/PPPM parameters
> > pme_order = 6
> > ewald_rtol = 1e-6
> > ewald_geometry = 3d
> > epsilon_surface = 0
> > ; OPTIONS FOR WEAK COUPLING ALGORITHMS
> > ; Temperature coupling
> > tcoupl = No
> > ; Groups to couple separately
> > tc_grps = System
> > ; Time constant (ps) and reference temperature (K)
> > tau_t = 1.0
> > ref_t = 300
> > ; Pressure coupling
> > Pcoupl = no
> > Pcoupltype = isotropic
> > ; Time constant (ps), compressibility (1/bar) and reference P (bar)
> > tau-p = 1.0
> > compressibility = 4.5e-5
> > ref-p = 1
> > ; Scaling of reference coordinates, No, All or COM
> > refcoord_scaling = No
> > ; Random seed for Andersen thermostat
> > andersen_seed = 124821
> >
> > gen-vel = no
> > gen-temp = 300
> > gen-seed = -1
> >
> > ; OPTIONS FOR BONDS
> > constraints = hbonds
> > ; Type of constraint algorithm
> > constraint-algorithm = Lincs
> > ; Do not constrain the start configuration
> > continuation = yes
> > ; Use successive overrelaxation to reduce the number of shake iterations
> > Shake-SOR = no
> > ; Relative tolerance of shake
> > shake-tol = 1e-04
> > ; Highest order in the expansion of the constraint coupling matrix
> > lincs-order = 12
> > ; Number of iterations in the final step of LINCS. 1 is fine for
> > ; normal simulations, but use 2 to conserve energy in NVE runs.
> > ; For energy minimization with constraints it should be 4 to 8.
> > lincs-iter = 1
> > ; Lincs will write a warning to the stderr if in one step a bond
> > ; rotates over more degrees than
> > lincs-warnangle = 30
> >
> > ; Free energy control stuff
> > free-energy = yes
> > init-lambda = .05
> > delta-lambda = 0
> > sc-alpha = 0.0
> > sc-power = 0
> > sc-sigma = 0.0
> > couple-moltype = Protein
> > couple-lambda0 = vdw
> > couple-lambda1 = vdw-q
> > couple-intramol = no
> >
>
> -- ========================================
>
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> ========================================
> -- gmx-users mailing list gmx-users at gromacs.org
> <mailto:gmx-users at gromacs.org>
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-request at gromacs.org
> <mailto:gmx-users-request at gromacs.org>.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list