[gmx-users] Huge acceleration needed to reproduce results!
Jennifer Williams
Jennifer.Williams at ed.ac.uk
Tue Jul 28 12:34:25 CEST 2009
Hi,
Yes I realised that gromacs works in ps. I converted my force in kj
mol-1 A-1 to acceleration in nm/ps2. I also took into account that the
msd.xvg is plotted in nm and ps-2 and the calculated gradient printed
at the top of the msd.xvg file is in cm2/s.
One strange thing that I do get is the message ?There were 228
inconsistent shifts. Check your topology? when I carry out the g_msd
with the ?mol option but not when I don?t use -mol. Why is this?
I also came across a forum post
(http://www.mail-archive.com/gmx-users@gromacs.org/msg11115.html) that
said ?If the distance between two atoms is close to half the box, the
force may arbitrarily change sign. This is an ill-defined situation
for which there is no obvious solution.? Could this somehow be
affecting my simulations?
Below are the relevant parts of my .mdp file and other files. If you
see something suspicious please let me know because I?m stuck,
Thanks
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.001
nsteps = 1000000
; 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 =
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 1000
nstvout = 1000
nstfout = 0
; Output frequency for energies to log file and energy file
nstlog = 1000
nstenergy = 1000
; 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 =
; ns algorithm (simple or grid)
ns_type = grid
; Periodic boundary conditions: xyz, no, xy
pbc = xyz
periodic_molecules = yes
; nblist cut-off
rlist = 0.9
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = PME
rcoulomb-switch = 0
rcoulomb = 0.9
; Relative dielectric constant for the medium and the reaction field
epsilon_r =
epsilon_rf =
; Method for doing Van der Waals
vdw-type = Cut-off
; cut-off lengths
rvdw-switch = 0
rvdw = 0.9
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = No
; Extension of the potential lookup tables beyond the cut-off
table-extension =
; 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 =
ewald_rtol = 1e-05
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = yes
; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl = nose-hoover
; Groups to couple separately
tc-grps = System
; Time constant (ps) and reference temperature (K)
tau_t = 0.1
ref_t = 150
; Pressure coupling
Pcoupl = No
Pcoupltype =
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p =
compressibility =
ref-p =
; Scaling of reference coordinates, No, All or COM
refcoord_scaling = no
; Random seed for Andersen thermostat
andersen_seed =
; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel = yes
gen_temp = 150
gen_seed = 173529
; OPTIONS FOR BONDS
constraints = none
; Type of constraint algorithm
constraint-algorithm = Lincs
; Do not constrain the start configuration
continuation = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 0.0001
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 4
; 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
; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp_excl =
; Non-equilibrium MD stuff
acc-grps = CH4
accelerate = 0.0 0.0 -200.0
freezegrps = SI_O
freezedim = Y Y Y
cos-acceleration = 0
deform =
and a shortened version of my top file:
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 no 1.0 1.0
;
;
[ atomtypes ]
; type mass charge ptype c6 c12
SI 28.08 1.28 A 0.000 0.00
O 15.999 -0.64 A 0.2708 1.538176
OH 15.999 -0.53 A 0.30 1.538176
H 1.008 0.21 A 0.000 0.000
;
; Include forcefield parameters
#include "CH4.itp"
;
;
[ moleculetype ]
; Name nrexcl
MCM 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 SI 1 MCM SI 1 1.2804993 28.086 2 SI 1 MCM SI 2 1.2804993 28.086 ......
287 SI 1 MCM SI 287 1.2804993 28.086
288 SI 1 MCM SI 288 1.2804993 28.086
289 O 1 MCM O 289 -0.64024965 15.9994
290 O 1 MCM O 290 -0.64024965 15.9994
.....
1070 OH 1 MCM OH 1070 -0.52612471 15.9994
1071 H 1 MCM H 1071 0.20599988 1.00797
[ bonds ]
; ai aj funct c0 c1 c2 c3
286 289 6 0.16 260510
8 1058 6 0.16 260510
......
[ constraints ]
; ai aj funct c0 c1 c2 c3
796 797 1 0.0945
798 799 1 0.0945
[ angles ]
; ai aj ak funct c0 c1
365 1 414 1 109.04 416.8167
365 1 631 1 109.04 416.8167
?..
537 288 728 1 109.04 416.8167
592 288 728 1 109.04 416.8167
225 289 286 1 72.364 180
66 290 117 1 72.364 180
....
9 794 239 1 72.364 180
175 795 228 1 72.364 180
218 796 797 1 108.5 460.954
2 798 799 1 108.5 460.954
[ dihedrals ]
; ai aj ak al funct c0 c1
709 8 1058 1059 5 0 0 0 0
403 8 1058 1059 5 0 0 0 0
603 8 1058 1059 5 0 0 0 0
962 4 1044 1045 5 0 0 0 0
366 4 1044 1045 5 0 0 0 0
524 4 1044 1045 5 0 0 0 0
1012 259 295 150 5 0.0000 0.0000 0.0000
[ system ]
; Name
CH4 in MCM
[ molecules ]
; Compound #mols
MCM 1
CH4 340
And my .itp file
[ atomtypes ]
; type mass charge ptype c6 c12
CH4 16.043 0.00 A 0.3732 1.24650457
;
[ moleculetype ]
; name nrexcl
CH4 2
[ atoms ]
; nr type resnr residu atom cgnr charge mass
1 CH4 1 CH4 CH4 1 0.00 16.043
Quoting Chris Neale <chris.neale at utoronto.ca>:
> Gromacs uses nm as the unit of distance. Did you account for that? If
> so, please add some .mdp file snippits and any other relevant files so
> that we can see directly what you are doing.
>
> Chris.
>
> --- original message ---
>
> Hello users,
>
> I am trying to reproduce a calculation that I carried out in DL_POLY.
> It is to calculate the transport diffusion coefficient for CH4 in a
> frozen mesoporous silica.
>
> In DL_POLY I used an external force of 0.1 kJ mol-1 A-1. (0.1 KJ per
> mole per angstrom). This equates to 10 dl_poly internal units which I
> add in this way at each timestep;
>
> Fsum = Fsum + Fex
>
> In Gromacs, I want to apply the same force as I used in DL_POLY so I
> calculated the required acceleration using F=ma. Where I took the mass
> to be the mass of one molecule of methane (16 a.m.u).
>
> The final value for acceleration that I came up with (which
> corresponds to a force of 0.1kj mol-1 A-1 on each molecule) was 0.0625
> nm ps-2.
>
> The first hiccup was when I used this value, the MSD was negative
> (though linear in the negative region of the graph). I assumed that
> this had something to do with the orientation of the unit cell and
> tried applying 0.0 0.0 -0.0625. The plot then looked much better.
>
> The problem is when I calculate the Mean displacement of the CH4
> molecules. (I do this using a slightly altered version of the g_msd
> code). The Mean displacement from gromacs is very different to that
> which I calculate using DL_POLY,
>
> Gromacs gives 95.0, dl_poly 21347.0.
>
> The MSD however (where I don?t add an acceleration are similar) so the
> problem lies with the force I am adding.
>
> To test that it wasn?t some bug in my code to calculate the Mean
> displacement, I also looked at how the acceleration/force altered the
> MSD in DL_POLY and gromacs.
>
> In DL_POLY, adding an external force of 0.1kj mol-1 A-1 would change
> the MSD of methane by 3 orders of magnitude compared to a run with no
> force added.
>
> My equivalent acceleration of -0.0625 in gromacs, in comparison,
> barely changes the MSD from that of a run with no acceleration added.
> In fact it takes an acceleration of -200 in the z direction to cause
> such a difference in the Ds coefficients between runs with and without
> acceleration.
>
> Does anyone have any idea what is going on here? An acceleration of
> 200 ps nm-2 surely is not reasonable is it?. It seems very large
> compared to the example in the manual of 0.1. This would then imply
> that my back of the envelope calculation for relating force and
> acceleration is wrong. Am I missing something? I?m quite sure that the
> force I am adding in DL_POLY is equivalent to 0.1 kJ mol-1 A-1 so why
> are my methane molecules moving so much less in gromacs in response to
> the equivalent acceleration?
>
> Also I noticed that although in the .mdp file I specify:
>
> ; Non-equilibrium MD stuff
> acc-grps = CH4
> accelerate = 0.0 0.0 -200.00
>
> In the md.log file I get the following output
>
> acc: 0 0 -154.549 0 0 45.4514
>
> Can someone clarify what this means?
>
> Any advice/comments appreciated
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the gromacs.org_gmx-users
mailing list