[gmx-users] Huge acceleration needed to reproduce results!
Jennifer Williams
Jennifer.Williams at ed.ac.uk
Tue Jul 28 18:11:47 CEST 2009
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.neale at 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.Williams at 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.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.
>>
>>
>>
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.
More information about the gromacs.org_gmx-users
mailing list