[gmx-users] Huge acceleration needed to reproduce results!
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Tue Jul 28 18:27:42 CEST 2009
As for your posted problem, I don't think that I can be of any more
assistance. I picked the low-hanging fruit, but probably you need some
assistance from somebody who has actually done constant acceleration
simulations.
I will take this opportunity to note, though, that you will never get
diffusion information from a constant acceleration simulation (as far
as I can figure). Although I am certain that you require a resolution
to your posted problem, I suggest that you concurrently ensure that
your approach is even warranted.
Chris.
Quoting Jennifer Williams <Jennifer.Williams at ed.ac.uk>:
>
>
> 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