Forgot column labels on energies. <div><br></div><div> Average RMSD Fluct. Drift Tot-Drift<div>2 fs time step Total Energy 46497.2 36445.6 2630.7 -148.833 -148833</div>
<div>4 fs time step Total Energy 41580.8 37693.4 2856.08 -130.198 -130198</div><div><br></div></div><div><div class="gmail_quote">On Tue, Mar 10, 2009 at 7:17 PM, Ilya Chorny <span dir="ltr"><<a href="mailto:ichorny@gmail.com">ichorny@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">Hello All,<div><br></div><div>I am trying to run some calibration calculations with 2/4 fs time steps. I am trying to reproduce the results in the Gromacs 4.0 paper on a protein/water (not the same protein) system with ~100K atoms. I ran 1 ns simulation in the NVE ensemble. My mdp params are shown below. Topolgy was created using -vsite h. My system has 218441 degrees of freedom according to GMXDUMP. My P.E. energy is ~ -1E5 and has not full relaxed after 1ns but is close. </div>
<div><br></div><div>2 fs time step Total Energy 46497.2 36445.6 2630.7 -148.833 -148833</div><div>4 fs time step Total Energy 41580.8 37693.4 2856.08 -130.198 -130198</div>
<div><br></div><div>Drift (kt/ns) 2fs/4fs = .34/.29</div><div><br></div><div>The drift seems to be much larger for both the 2fs and 4fs time steps then in the paper. Do I have a clear bug in my params? Should I wait till the system fully relaxes before doing the drift calculation?</div>
<div><br></div><div><br></div><div><div>integrator = md</div><div>; Start time and timestep in ps</div><div>tinit = 0</div><div>dt = 0.002/0.004</div><div>nsteps = 500000/250000</div>
<div><div>; For exact run continuation or redoing part of a run</div><div>; Part index is updated automatically on checkpointing (keeps files separate)</div><div>simulation_part = 1</div><div>init_step = 0</div>
<div>; mode for center of mass motion removal</div><div>comm_mode = linear</div><div>; number of steps for center of mass motion removal</div><div>nstcomm = 1</div><div>; group(s) for center of mass motion removal</div>
<div>comm_grps = system</div><div><div>; NEIGHBORSEARCHING PARAMETERS</div><div>; nblist update frequency</div><div>nstlist = 10/5</div><div>; ns algorithm (simple or grid)</div><div>ns_type = grid</div>
<div>; Periodic boundary conditions: xyz, no, xy</div><div>pbc = xyz</div><div>periodic_molecules = no</div><div>; nblist cut-off</div><div>rlist = 1.0</div><div><div>; OPTIONS FOR ELECTROSTATICS AND VDW</div>
<div>; Method for doing electrostatics</div><div>coulombtype = PME</div><div>rcoulomb-switch = 0</div><div>rcoulomb = 1.0</div><div>; Relative dielectric constant for the medium and the reaction field</div>
<div>epsilon-r = 80</div><div>epsilon_rf = 1</div><div>; Method for doing Van der Waals</div><div>vdw-type = Switch</div><div>; cut-off lengths</div><div>rvdw-switch = .9</div>
<div>rvdw = 1.0</div><div>; Apply long range dispersion corrections for Energy and Pressure</div><div>DispCorr = EnerPres</div><div>; Extension of the potential lookup tables beyond the cut-off</div>
<div>table-extension = 1</div><div>; Seperate tables between energy group pairs</div><div>energygrp_table =</div><div>; Spacing for the PME/PPPM FFT grid</div><div>fourierspacing = 0.12</div><div>
; FFT grid size, when a value is 0 fourierspacing will be used</div><div>fourier_nx = 0</div><div>fourier_ny = 0</div><div>fourier_nz = 0</div><div>; EWALD/PME/PPPM parameters</div>
<div>pme_order = 4</div><div>ewald_rtol = 1.e-04</div><div>ewald_geometry = 3d</div><div>epsilon_surface = 0</div><div>optimize_fft = no</div><div><div>; OPTIONS FOR WEAK COUPLING ALGORITHMS</div>
<div>; Temperature coupling</div><div>Tcoupl = no</div><div>; Groups to couple separately</div><div>tc-grps = System</div><div>; Time constant (ps) and reference temperature (K)</div><div>
tau_t = 0.1</div><div>ref_t = 298.0</div><div>; Pressure coupling</div><div>Pcoupl = No</div><div>Pcoupltype = Isotropic</div><div>; Time constant (ps), compressibility (1/bar) and reference P (bar)</div>
<div>tau_p = 1.67</div><div>compressibility = 4.5e-5</div><div>ref_p = 1.01325</div><div>; Scaling of reference coordinates, No, All or COM</div><div>refcoord_scaling = No</div>
<div>; Random seed for Andersen thermostat</div><div>andersen_seed = 815131</div><div><div>; GENERATE VELOCITIES FOR STARTUP RUN</div><div>gen_vel = yes</div><div>gen_temp = 298.0</div>
<div>gen-seed = 173529</div><div><div>; OPTIONS FOR BONDS</div><div>constraints = hbonds/all-bonds</div><div>; Type of constraint algorithm</div><div>constraint-algorithm = lincs</div><div>
; Do not constrain the start configuration</div><div>continuation = no</div><div>; Use successive overrelaxation to reduce the number of shake iterations</div><div>Shake-SOR = no</div><div>; Relative tolerance of shake</div>
<div>shake-tol = 0.0001</div><div>; Highest order in the expansion of the constraint coupling matrix</div><div>lincs-order = 6</div><div>; Number of iterations in the final step of LINCS. 1 is fine for</div>
<div>; normal simulations, but use 2 to conserve energy in NVE runs.</div><div>; For energy minimization with constraints it should be 4 to 8.</div><div>lincs-iter = 2</div><div>; Lincs will write a warning to the stderr if in one step a bond</div>
<div>; rotates over more degrees than</div><div>lincs-warnangle = 30</div><div>; Convert harmonic bonds to morse potentials</div><div>morse = no</div><div><br></div></div></div></div></div></div>
</div><br>-- <br>Ilya Chorny Ph.D.<br><br>
</div>
</blockquote></div><br><br clear="all"><br>-- <br>Ilya Chorny Ph.D.<br><br>
</div>