[gmx-users] Lincs error and explosion of the system

[Justin A. Lemkul]

Why don't you do energy minimization?  It's pretty clear that you've generated
some atomic overlap in your procedure, from what you've described.

For general reference on such problems, check the hundreds of related posts
(LINCS warnings) in the list archive or the following link to the wiki site:

http://wiki.gromacs.org/index.php/blowing_up

-Justin

Quoting Egidijus Kuprusevicius <ekuprusevicius at yahoo.com>:

> Dear users,
> I have a problem on running MD of 146 LC molecules system. I do not do the
> energy minimisation. I use NVT ensemble to equilibrate my system. My topology
> was created using PRODRG beta server. After the creation of my systems box
> (on top of each other)
>
> 1. editconf -f lc -bt d -box 1.3 0.02 0.02 -princ -c -d 0.5 -rvdw 0.12 -o lcc
> 2. genconf -f lcc -nbox 8 8 8 -dist 0 0 0 -nmolat 27 -o lc512
>
> I put inside another (a little bigger) molecule
>
> 3. editconf -f sl -bt c -box 6 -princ -c -d 2.5 -rvdw 0.12 -o sll
> 4. genbox -cp sll -cs lc512 -o slw -p sl
>
> then create tpr
>
> 5. grompp -v -np 32 -shuffle -sort -f eq.mdp -c slw -p sl -o topol.tpr
> and submit the job
>
> I always get my system exploding
>
> Initializing LINear Constraint Solver
> number of constraints is 122
> average number of constraints coupled to one constraint is 3.0
> Rel. Constraint Deviation: Max between atoms RMS
> Before LINCS 0.086035 40 39 0.026421
> After LINCS 0.000704 103 104 0.000292
> Energies (kJ/mol)
> G96Angle Proper Dih. Improper Dih. LJ-14 Coulomb-14
> 2.76783e+03 8.11377e+01 8.21662e+01 3.98602e+04 5.57164e+03
> LJ (SR) LJ (LR) Coulomb (SR) Coul. recip. Potential
> 2.71048e+03 -1.88575e+02 -3.38992e+03 -3.62073e+04 1.12877e+04
> Kinetic En. Total Energy Temperature Pressure (bar)
> 6.65517e+04 7.78394e+04 2.06669e+03 1.75872e+05
>
> Step 86, time 0.172 (ps) LINCS WARNING
> relative constraint deviation after LINCS:
> max 96.909290 (between atoms 85 and 86) rms 10.933811
> bonds that rotated more than 30 degrees:
> atom 1 atom 2 angle previous, current, constraint length
> 77 78 113.7 0.1390 0.5917 0.1390
> 77 87 116.3 0.1390 0.5570 0.1390
> 78 79 109.8 0.1090 0.5917 0.1090
> 78 80 104.4 0.1390 1.2899 0.1390
> 80 81 93.6 0.1090 5.0458 0.1090
> 82 80 99.8 0.1390 3.9224 0.1390
> 82 83 100.5 0.1430 3.8921 0.1430
> 82 85 99.9 0.1390 3.9638 0.1390
> 83 84 103.4 0.1140 3.1119 0.1140
> 85 86 91.7 0.1090 10.6721 0.1090
> 85 87 100.2 0.1390 1.7797 0.1390
> 87 88 110.3 0.1090 0.5717 0.1090
> Constraint error in algorithm Lincs at step 86
> Wrote pdb files with previous and current coordinates
> Large VCM(group rest):
>
-949890684120872786061526585043023133196136319944012857455888493019588547485855908650717825011817737590092000256377384500897124503083675335640562545921274578670325230339034153228858928503740516247758241500986457604229169152.00000,
>
-7700177673405372508341068198495354437944151718908865834595177256517607987036524874348348597888871194101005699859793685064099787214560958505758284878533509661264882479706501126070295847653568962516830194295490734095350104064.00000,
>
-17482030675841594076114066100363479435124387445824631353703195421406478322250794440904806833364056039859394683424411714469447896344546474426178109668279235955517800609185251692588873229116161214420413237478774727134943379456.00000,
> ekin-cm: inf
>
> MDP file
>
> title = molecular dynamics
> cpp = /usr/bin/cpp
> include = -I../top
> ; RUN CONTROL PARAMETERS
> integrator = md
> dt = 0.002 ;ps
> nsteps = 5000000 ;10 ns
> ;comm_mode = Angular ;remove center of mass translation and rotation around
> the center of mass
> nstcomm = 1 ;number of steps for center of mass motion removal
> ; OUTPUT CONTROL OPTIONS
> nstxout = 100 ;every 0.2 ps
> nstvout = 100
> nstfout = 0
> nstlog = 100 ;Output frequency to write energies to log file
> nstenergy = 100 ;Output frequency to write energies to energy file
> energygrps = MMM DRG ;Selection of energy groups
> ; NEIGHBOR SEARCHING PARAMETERS
> nstlist = 1 ;update frequency
> ns_type = grid ;algorithm (simple or grid)
> pbc = xyz
> rlist = 1.2 ;cut-off distance nm
> ; OPTIONS FOR ELECTROSTATICS
> coulombtype = PME ;Method for doing electrostatics
> rcoulomb = 1.2 ;Coulomb cut-off distance nm
> epsilon-r = 1.0 ;relative dielectric constant
> ; VDW
> vdw-type = Cut-off ;Method for doing VDW
> rvdw = 1.5 ;cut-off lengths nm
> ;DispCorr = EnerPres ;Apply long range dispersion corrections for Energy and
> Pressure
> fourierspacing = 0.1 ;Spacing for the PME/PPPM FFT grid
> ; EWALD/PME/PPPM parameters
> pme_order = 4 ;cubic interpolation
> ewald_rtol = 1e-05 ;relative Evald-shifted potential at cut-off
> optimize_fft = yes ;calculate at startup
> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
> Tcoupl = Berendsen ;Temperature coupling
> Tc-grps = System ;Groups to couple separately
> tau_T = 5.0 ;Time constant (ps)
> ref_T = 299 ;reference temperature (K)
> Pcoupl = no ;Pressure coupling
> Pcoupltype = Isotropic
> tau_P = 5.0 ;Time constant (ps)
> compressibility = 4.5e-5 ;compressibility (1/bar)
> ref_P = 1.0 ;reference P (bar)
> ; OPTIONS FOR BONDS
> constraints = all-bonds
> constraint-algorithm = lincs
> ;shake_tol = 0.0001
> unconstrained-start = yes ;don't constrain the start configuration
> lincs-order = 4 ;Highest order in the expansion of the constraint coupling
> matrix
> lincs-warnangle = 30
> lincs-iter = 2
> morse = no ;Convert harmonic bonds to morse potentials
> ;GENERATE VELOCITES
> gen_vel = no ;generates velocities
> gen_temp = 299
> gen_seed = 173529
>
> Could anyone help me?
>
> Thanks
>
> Egidijus Kuprusevicius
>
>
>
>
>       Be a better friend, newshound, and
> know-it-all with Yahoo! Mobile.  Try it now.
>



========================================

Justin A. Lemkul
Graduate Research Assistant
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul at vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/

========================================