[gmx-users] Huge acceleration needed to reproduce results!

chris.neale at utoronto.ca chris.neale at utoronto.ca
Tue Jul 28 16:30:19 CEST 2009


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.
>
>
>






More information about the gromacs.org_gmx-users mailing list