<br><br>----- Original Message -----<br>From: teklebrh@ualberta.ca<br>Date: Monday, June 28, 2010 6:01<br>Subject: [gmx-users] Error in MD simulation<br>To: jalemkul@vt.edu, Discussion list for GROMACS users &lt;gmx-users@gromacs.org&gt;<br><br>&gt; Dear Gromacs users,<br>&gt; <br>&gt; First I have to mention few points. I read all the LINCS error <br>&gt; and tried a number of options.<br>&gt; <br>&gt; 1, Minimized well the system. and the Epot and Max Forces are <br>&gt; very resonable.<br>&gt; <br>&gt; Steepest Descents converged to Fmax &lt; 100 in 379 steps<br>&gt; Potential Energy&nbsp; = -1.0573750e+05<br>&gt; Maximum force&nbsp;&nbsp;&nbsp;&nbsp; =&nbsp; 8.0876366e+01 on <br>&gt; atom 1117<br>&gt; Norm of force&nbsp;&nbsp;&nbsp;&nbsp; =&nbsp; 8.7360611e+00<br>&gt; <br>&gt; 2, Equlibrate well the system. It works for 100ps<br>&gt; <br>&gt; Once I run NPT simulation using the following .mdp parameter set <br>&gt; it shows the following error.<br>&gt; <br>&gt; ; The etire lines starting with ';' ARE COMMENTS<br>&gt; <br>&gt; title&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;        = NPT <br>&gt; equilibration using the berendsen thermostat and barostat         <br>&gt; Title of run<br>&gt; cpp&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = /usr/bin/cpp ; location of cpp on linux<br>&gt; <br>&gt; ; The following lines tell the program the standard locations <br>&gt; where to find certain files<br>&gt; <br>&gt; ; VARIOUS PREPROCESSING OPTIONS<br>&gt; ; Preprocessor information: use cpp syntax.<br>&gt; ; e.g.: -I/home/joe/doe -I/home/mary/hoe<br>&gt; include&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; =<br>&gt; <br>&gt; ; e.g.: -DI_Want_Cookies -DMe_Too<br>&gt; <br>&gt; ; RUN CONTROL PARAMETERS<br>&gt; <br>&gt; integrator&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = md ; leap-frog integrator<br>&gt; ; Start time and timestep in ns<br>&gt; tinit&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 0<br>&gt; dt&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 0.002&nbsp;&nbsp;&nbsp;&nbsp; ; 2fs<br>&gt; nsteps&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 500000&nbsp;&nbsp; ; 2*500 = 1ns<br>&gt; <br>&gt; ; THOSE OPTIONS REMOVE MOTION OF THE ASPHALTENE MODEL COMPUNDS <br>&gt; RELATIVE TO THE SOLVENT/IONS<br>&gt; <br>&gt; ; mode for center of mass motion removal<br>&gt; comm-<br>&gt; mode&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = Linear<br>&gt; ; number of steps for center of mass motion removal<br>&gt; nstcomm&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1<br>&gt; ; group(s) for center of mass motion removal<br>&gt; comm-<br>&gt; grps&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = System<br>&gt; <br>&gt; ; LANGEVIN DYNAMICS OPTIONS<br>&gt; ; Friction coefficient (amu/ps) and random seed<br>&gt; bd-<br>&gt; fric&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 0<br>&gt; ld-<br>&gt; seed&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1993<br>&gt; <br>&gt; ; ENERGY MINIMIZATION OPTIONS<br>&gt; ; Force tolerance and initial step-size<br>&gt; emtol&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 200<br>&gt; emstep&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 0.01<br>&gt; <br>&gt; ; Frequency of steepest descents steps when doing CG<br>&gt; nstcgsteep&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1000<br>&gt; nbfgscorr&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 10<br>&gt; <br>&gt; ; OUTPUT CONTROL OPTIONS<br>&gt; <br>&gt; ; Output frequency for coords (x), velocities (v) and forces (f)<br>&gt; nstxout&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 2000<br>&gt; nstvout&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 2000<br>&gt; nstfout&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 0<br>&gt; ; Output frequency for energies to log file and energy file<br>&gt; nstlog&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1000<br>&gt; nstenergy&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1000<br>&gt; ; Output frequency and precision for xtc file<br>&gt; nstxtcout&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1000<br>&gt; xtc-<br>&gt; precision&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1000<br>&gt; <br>&gt; <br>&gt; ; This selects the subset of atoms for the xtc file. You can<br>&gt; ; select multiple groups. By default all atoms will be written.<br>&gt; xtc-<br>&gt; grps&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = BIS&nbsp; HEP<br>&gt; ; Selection of energy groups<br>&gt; energygrps&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = BIS&nbsp; HEP<br>&gt; <br>&gt; ; NEIGHBORSEARCHING PARAMETERS<br>&gt; <br>&gt; ; nblist update frequency<br>&gt; nstlist&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 5<br>&gt; ; ns algorithm (simple or grid)<br>&gt; ns-<br>&gt; type&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = grid<br>&gt; ; Periodic boundary conditions: xyz, no, xy<br>&gt; pbc&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = xyz<br>&gt; ; nblist cut-off<br>&gt; rlist&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1.2<br>&gt; <br>&gt; ; OPTIONS FOR ELECTROSTATICS AND VDW<br>&gt; <br>&gt; ; Method for doing electrostatics<br>&gt; coulombtype&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = PME<br>&gt; rcoulomb&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1.2<br>&gt; <br>&gt; <br>&gt; ; Method for doing Van der Waals<br>&gt; vdw-<br>&gt; type&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = Cut-off<br>&gt; ; cut-off lengths<br>&gt; rvdw&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1.2<br>&gt; <br>&gt; ; Apply long range dispersion corrections for Energy and Pressure<br>&gt; DispCorr&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = EnerPres ; account for the cut-off vdw scheme<br>&gt; <br>&gt; ; Spacing for the PME/PPPM FFT grid<br>&gt; fourierspacing&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 0.16<br>&gt; <br>&gt; ; EWALD/PME/PPPM parameters<br>&gt; pme_order&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 4<br>&gt; ewald_rtol&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1e-05<br>&gt; ewald_geometry&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 3d<br>&gt; <br>&gt; ; IMPLICIT SOLVENT ALGORITHM<br>&gt; implicit_solvent&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <br>&gt; = No<br>&gt; <br>&gt; ; OPTIONS FOR WEAK COUPLING ALGORITHMS<br>&gt; ; Temperature coupling is on<br>&gt; tcoupl&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = V-rescale; modified berendsen thermostat<br>&gt; ; Groups to couple separately<br>&gt; tc-<br>&gt; grps&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = BIS&nbsp; HEP ; three coupling gropus<br>&gt; ; Time constant (ps) and reference temperature (K)<br>&gt; tau-<br>&gt; t&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 0.1&nbsp; 0.1&nbsp;&nbsp; ; time constant<br>&gt; ref-<br>&gt; t&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 298&nbsp; 298&nbsp;&nbsp; ; reference temprature<br>&gt; ; Pressure coupling is on<br>&gt; Pcoupl&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = berendsen; Pressure coupling is on<br>&gt; Pcoupltype&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = Isotropic<br>&gt; ; Time constant (ps), compressibility (1/bar) and reference P (bar)<br>&gt; tau-<br>&gt; p&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 3&nbsp; 3<br>&gt; compressibility&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1.47e-4&nbsp; 1.47e-4&nbsp; ; isothermal compressibility, bar^-1<br>&gt; ref-<br>&gt; p&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 1&nbsp; 1<br>&gt; ; Random seed for Andersen thermostat<br>&gt; andersen_seed&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 815131<br>&gt; <br>&gt; ; GENERATE VELOCITIES FOR STARTUP RUN<br>&gt; gen-<br>&gt; vel&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = yes ; assign velocities from Maxwell distribution (MXD)<br>&gt; gen-<br>&gt; temp&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 298 ; temperature for MXD<br>&gt; gen-<br>&gt; seed&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 173529&nbsp;&nbsp;&nbsp; ; used to initialize random generator for random velocities<br>&gt; <br>&gt; ; OPTIONS FOR BONDS<br>&gt; constraints&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = all-bonds<br>&gt; ; Type of constraint algorithm<br>&gt; constraint-algorithm&nbsp;&nbsp;&nbsp;&nbsp; = Lincs<br>&gt; ; Do not constrain the start configuration<br>&gt; ; Highest order in the expansion of the constraint coupling matrix<br>&gt; lincs-<br>&gt; order&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 4<br>&gt; ; Number of iterations in the final step of LINCS. 1 is fine for<br>&gt; ; normal simulations, but use 2 to conserve energy in NVE runs.<br>&gt; ; For energy minimization with constraints it should be 4 to 8.<br>&gt; lincs-<br>&gt; iter&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 2<br>&gt; ; Lincs will write a warning to the stderr if in one step a bond<br>&gt; ; rotates over more degrees than<br>&gt; lincs-<br>&gt; warnangle&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = 30<br>&gt; ; Convert harmonic bonds to morse potentials<br>&gt; morse&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = no<br>&gt; <br>&gt; ==========<br>&gt; <br>&gt; It run the simulation for quite sometime but it gives me an <br>&gt; error of LINCS some how in the middle.<br>&gt; <br>&gt; I try changing some parameter sets like step size,tau-p , tau-t, <br>&gt; lincs-warnangle but did not overcame the problem.<br><br>You are trying to equilibrate straight into NPT, which can be problematic if (for example) the density isn't right. Consider the advice in point 8 of http://www.gromacs.org/Documentation/How-tos/Steps_to_Perform_a_Simulation<br><br>Mark<br><br>&gt; Step 197685, LINCS WARNING<br>&gt; &gt;relative constraint deviation after LINCS:<br>&gt; &gt;rms 0.000000, max 0.000001 (between atoms 1856 and 15982)<br>&gt; &gt;bonds that rotated more than 30 degrees:<br>&gt; &gt; atom 1 atom 2&nbsp; angle&nbsp; previous, current, constraint length<br>&gt; &gt;&nbsp;&nbsp;&nbsp; 97&nbsp;&nbsp;&nbsp;&nbsp; 98&nbsp;&nbsp; <br>&gt; 35.3&nbsp;&nbsp;&nbsp; 0.1000&nbsp;&nbsp; <br>&gt; 0.1000&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.1000<br>&gt; <br>&gt; can some one give me an advice or help?<br>&gt; <br>&gt; For my other 10 MD set up it works fine but with diffrent solvents.<br>&gt; <br>&gt; rob<br>&gt; -- <br>&gt; gmx-users mailing list&nbsp;&nbsp;&nbsp; gmx-users@gromacs.org<br>&gt; http://lists.gromacs.org/mailman/listinfo/gmx-users<br>&gt; Please search the archive at http://www.gromacs.org/search <br>&gt; before posting!<br>&gt; Please don't post (un)subscribe requests to the list. Use thewww <br>&gt; interface or send it to gmx-users-request@gromacs.org.<br>&gt; Can't post? Read http://www.gromacs.org/mailing_lists/users.php