Hi Chris,
Thanks for the input-its good to get another person's perspective.
Let me explain a bit more about my system...
MCM stands for MCM-41, a mesoporous silica which consits of 1071 atoms
so even though it is only 1 molecule, it is the more massive group.
I am purposely accelerating the methane to try and calculate their
diffusion thourgh the pore of the MCM.
As I understand it, the acceleration I specify should be applied to
each and every molecule of CH4. I.e if I specify an acceleration of
0.1 nm ps-2, each of my 340 molecules should experience an
acceleration of 0.1. I.e the acceleration is not in any way divided up
between the total number of accelerated molecules?
Jenny
Quoting chris.ne...@utoronto.ca:
Title indicates you think that you have CH4 in MCM:
[ system ]
; Name
CH4 in MCM
Actual system is MCM in CH4:
[ molecules ]
; Compound #mols
MCM 1
CH4 340
.mdp accelerates the CH4:
acc-grps = CH4
Now I'm not sure if this is the source of any problem, equal and
opposite being true, but you are certainly trying to accelerate the
more massive group and it seems like a strange thing to do. Could this
be it?
Chris.
-- original message --
Quoting Jennifer Williams <jennifer.willi...@ed.ac.uk>:
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.ne...@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.
Dr. Jennifer Williams
Institute for Materials and Processes
School of Engineering
University of Edinburgh
Sanderson Building
The King's Buildings
Mayfield Road
Edinburgh, EH9 3JL, United Kingdom
Phone: ++44 (0)131 650 4 861
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
_______________________________________________
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
