[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