[gmx-users] inconsistency in obtained temperature values by traj and g_energy tools

Rasoul Nasiri nasiri1355 at gmail.com
Sat Jul 12 00:10:15 CEST 2014


On Fri, Jul 11, 2014 at 10:27 PM, Justin Lemkul <jalemkul at vt.edu> wrote:

>
>
> On 7/10/14, 9:30 PM, Rasoul Nasiri wrote:
>
>> On Fri, Jul 11, 2014 at 12:58 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>
>>
>>>
>>> On 7/10/14, 7:09 PM, Rasoul Nasiri wrote:
>>>
>>>  Because
>>>>
>>>> 1- Estimation of T for different sub-systems say phases of gas,
>>>> interface
>>>> and liquid.
>>>> Can be obtained with g_energy?
>>>>
>>>>
>>>>  No, not unless they are assigned as separate tc-grps in the .mdp file,
>>>
>>
>>
>> as sub-systems are changing during the simulation in terms of number of
>> atoms in a non-equilibrium situation, I don't think g_energy can give
>> relevant temperature values, even tc-grps are specified in mdp. Isn't?
>>
>>
> No, you can't then.  This also makes any sort of interpretation
> non-trivial; if an atom dissociates from one group to which it was linked,
> the interpretation of its new velocity when bonding to another group is now
> somewhat discontinuous, is it not?
>
>
>
So it could make clrea why I did not select g_energy for determination of T
and replied your question!


>
>>
>>  but then (1) you shouldn't be designing something affecting the physics
>>> purely based on what you hope to analyze
>>>
>>
>>
>>
>>
>>
>>  and (2) you then already know the outcome.
>>>
>>> This is actually what g_traj can do; like I said, you just have to
>>> account
>>> for restraints.
>>>
>>
>>
>> can you clarify how it can be practically done?
>>
>>
> You need to recalculate the number of degrees of freedom based on whatever
> constraints there might be.  You said there are none, so that may not be
> the case.  Please excuse the typo in the earlier message, "restraints"
> should be "constraints."


Which DOFs, do I need consider translational and rotational as well as
vibrational ones? And how this number will affact crude values obtained by
g_traj? Please consider constrained MD simulation done using OPLS.


>
>
>
>>
>>>
>>>   2- My trajectory is not constrain one as it has been recorded using
>>>
>>>> reactive FF.
>>>>
>>>> I need to estimate T with converted trajectory (reaxff---->gro).
>>>>
>>>>
>>>>  Converted trajectory?
>>>
>>
>>
>> yes, since non-reactive FF at high temperatures cannot be reliable due to
>> vibrational effects, all my MD has been done by reactive one.
>>
>>
>>    Your first post suggested the run was done with Gromacs (being that you
>>> had .tpr and .edr files).
>>>
>>
>>
>> I have produced .tpr file but not .edr.
>>
>>
>>    The outcomes you posted agree perfectly with constrained molecules,
>>>
>>
>>
>> The posted data refereed to un-reliability of g_traj for obtaining of T. I
>> just wanted to know the reason via comparison with g_energy. Those
>> obtained
>> by non-reactive FF of OPLS.
>>
>>
> This is where I was getting confused.  Your very first message said you
> had the "same trajectory and input files" but clearly this is not true.
>  You're comparing apples and oranges; an accurate statement of what's
> really going on is very important and saves everyone time.
>
>
That means all posted values in my first message obtained by ONLY OPLS FF
just to understand how much g_traj and g_energy give us different T values
on same trajectory and input file.



>
>
>>  but without more details on what you're actually doing, I'm not going to
>>> hazard a guess.  Recall that even if you set "constraints = none," SETTLE
>>> is still applied to water molecules unless you use -DFLEXIBLE (which you
>>> shouldn't for most normal purposes).
>>>
>>>
>>>
>>>  My molecules are not water and g_traj gave following T values for same
>> system using two different FFs
>>
>> 1- ReaxFF
>>   0      399.321
>>   1      400.25
>>   2      401.08
>>   3      401.271
>>   4      400.67
>>   5      400.264
>>   6      400.752
>>   7      401.178
>>   8      401.391
>>   9      399.065
>> 10      401.78
>>
>> As i said already I just converted reaxff trajectory to gro one and also
>> used tpr file (constrain=none in mdp) for estimation of T with g_traj.
>>
>> 2- OPLS
>>   0      269.581
>>   1      271.064
>>   2      271.309
>>   3      271.205
>>   4      271.499
>>   5      270.868
>>   6      272.247
>>   7      271.59
>>   8      271.614
>>   9      270.598
>>   10     271.662
>>
>>
>> How constrained MD results (OPLS) can be corrected?
>>
>>
> Now you're losing me again.  You're talking about converting some ReaxFF
> trajectory without constraints, but were there constraints in the OPLS
> trajectory?


here there are two different values of T obtained by ReaxFF and OPLS. As
you told there are no constrains in ReaxFF trajectory BUT there should be
constrain in OPLS trajectory.


>  Can you give an actual listing of commands and .mdp files (if applicable)
> for what you're doing?  This is all very confusing.
>
>
>

Please just tell me how values of obtained T by g_traj can be corrected
when I use trajectory and tpr files of constrain MD simulations (OPLS) as
an input?


command:
g_traj_d -f md.trr  -s topol  -n n.ndx -ot T.xvg


mdp file:

; RUN CONTROL PARAMETERS
integrator               = md
; Start time and timestep in ps
tinit                    = 0
dt                       = 0.002
nsteps                   = 10000000 ;20ns
; For exact run continuation or redoing part of a run
init_step                = 0
; Part index is updated automatically on checkpointing (keeps files
separate)
simulation_part          = 1
; mode for center of mass motion removal
comm-mode                = Linear
; number of steps for center of mass motion removal
nstcomm                  = 10
; 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

; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol                    = 10
emstep                   = 0.01
; Max number of iterations in relax_shells
niter                    = 20
; Step size (ps^2) for minimization of flexible constraints
fcstep                   = 0
; Frequency of steepest descents steps when doing CG
nstcgsteep               = 1000
nbfgscorr                = 10

; TEST PARTICLE INSERTION OPTIONS
rtpi                     = 0.05

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout                  = 100000 ; 200ps
nstvout                  = 100000
nstfout                  = 0
; Output frequency for energies to log file and energy file
nstlog                   = 5000   ;10ps
nstcalcenergy            = 10
nstenergy                = 5000   ;10ps
; Output frequency and precision for xtc file
nstxtcout                = 5000   ;10ps
xtc-precision            = 5000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc_grps                 = System
; 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.5
; long-range cut-off for switched potentials
rlistlong                = 1.5

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype              = PME
rcoulomb-switch          = 0
rcoulomb                 = 1.5
; Relative dielectric constant for the medium and the reaction field
epsilon_r                = 1
epsilon_rf               = 1
; Method for doing Van der Waals
vdw-type                 = switch
; cut-off lengths
rvdw-switch              = 1.1
rvdw                     = 1.3
; 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                = 4
ewald_rtol               = 1e-05
ewald_geometry           = 3d
epsilon_surface          = 0
optimize_fft             = no


; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl                   = v-rescale
nsttcouple               = -1
nh-chain-length          = 10
; Groups to couple separately
tc-grps                  = System
; Time constant (ps) and reference temperature (K)
tau_t                    = 0.1
ref_t                    = 400
; Pressure coupling
Pcoupl                   = parrinello-rahman
Pcoupltype               = isotropic
nstpcouple               = -1
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p                    = 4.0
compressibility          = 13.31e-5 ; cyclopentane
ref_p                    = 1.0
; Scaling of reference coordinates, No, All or COM
refcoord_scaling         = No
; Random seed for Andersen thermostat
andersen_seed            = 815131

; OPTIONS FOR BONDS
constraints              = all-bonds
; 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              = 6
; 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
; Convert harmonic bonds to morse potentials
morse                    = no


---------------------------------------------------------------------------


Rasoul





> -Justin
>
> --
> ==================================================
>
> Justin A. Lemkul, Ph.D.
> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 601
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
>
> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul
>
> ==================================================
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>


More information about the gromacs.org_gmx-users mailing list