Simulated annealing LINCS Errors

GROMACS version: 2020.3
GROMACS modification: No

I’m having the similar issues with a simulation from one of previous posts. I’ve heated up 2 polystyrene 100-mer chains to 1200 K with no problem, but when I try to run simulated annealing to cool it down to 200 K, there’s always a LINCS error to indicate my system’s unstable. I tried running an energy minimization with an emtol = 10 and even tried starting at a lower temperature like 900 K. Going up in temperature is no problem though. It’s just cooling down that has issues. Could anyone help?

Thanks!

Log file:

> GROMACS:      gmx mdrun, version 2020.3
> Executable:   /usr/local/gromacs/bin/gmx
> Data prefix:  /usr/local/gromacs
> Working dir:  /Users/profile1/Desktop/gromacs-related_work/NPs_BBB/PS/3nm-PS
> Process ID:   1754
> Command line:
>   gmx mdrun -deffnm testcool -v
> 
> GROMACS version:    2020.3
> Verified release checksum is c0599e547549c2d0ef4fc678dc5a26ad0000eab045e938fed756f9ff5b99a197
> Precision:          single
> Memory model:       64 bit
> MPI library:        thread_mpi
> OpenMP support:     enabled (GMX_OPENMP_MAX_THREADS = 64)
> GPU support:        OpenCL
> SIMD instructions:  AVX2_256
> FFT library:        fftw-3.3.8-sse2
> RDTSCP usage:       enabled
> TNG support:        enabled
> Hwloc support:      disabled
> Tracing support:    disabled
> C compiler:         /usr/local/bin/icc Intel 19.1.2.20200623
> C compiler flags:   -march=core-avx2 -std=gnu99 -ip -funroll-all-loops -alias-const -ansi-alias -no-prec-div -fimf-domain-exclusion=14 -qoverride-limits -O3 -DNDEBUG
> C++ compiler:       /usr/local/bin/icpc Intel 19.1.2.20200623
> C++ compiler flags: -march=core-avx2 -ip -funroll-all-loops -alias-const -ansi-alias -no-prec-div -fimf-domain-exclusion=14 -qoverride-limits -qopenmp -O3 -DNDEBUG
> OpenCL include dir: /Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.15.sdk/System/Library/Frameworks/OpenCL.framework
> OpenCL library:     /Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.15.sdk/System/Library/Frameworks/OpenCL.framework
> OpenCL version:     1.2
> 
> 
> Running on 1 node with total 4 cores, 4 logical cores, 0 compatible GPUs
> Hardware detected:
>   CPU info:
>     Vendor: Intel
>     Brand:  Intel(R) Core(TM) i5-4570R CPU @ 2.70GHz
>     Family: 6   Model: 70   Stepping: 1
>     Features: aes apic avx avx2 clfsh cmov cx8 cx16 f16c fma htt intel lahf mmx msr nonstop_tsc pcid pclmuldq pdcm pdpe1gb popcnt pse rdrnd rdtscp sse2 sse3 sse4.1 sse4.2 ssse3 tdt x2apic
>   Hardware topology: Only logical processor count
>   GPU info:
>     Number of GPUs detected: 1
>     #0: name: Iris Pro, vendor: Intel, device version: OpenCL 1.2 , stat: incompatible (please recompile with GMX_OPENCL_NB_CLUSTER_SIZE=4)
> 
> 
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, E.
> Lindahl
> GROMACS: High performance molecular simulations through multi-level
> parallelism from laptops to supercomputers
> SoftwareX 1 (2015) pp. 19-25
> -------- -------- --- Thank You --- -------- --------
> 
> 
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> S. Páll, M. J. Abraham, C. Kutzner, B. Hess, E. Lindahl
> Tackling Exascale Software Challenges in Molecular Dynamics Simulations with
> GROMACS
> In S. Markidis & E. Laure (Eds.), Solving Software Challenges for Exascale 8759 (2015) pp. 3-27
> -------- -------- --- Thank You --- -------- --------
> 
> 
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R.
> Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess, and E. Lindahl
> GROMACS 4.5: a high-throughput and highly parallel open source molecular
> simulation toolkit
> Bioinformatics 29 (2013) pp. 845-54
> -------- -------- --- Thank You --- -------- --------
> 
> 
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> B. Hess and C. Kutzner and D. van der Spoel and E. Lindahl
> GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable
> molecular simulation
> J. Chem. Theory Comput. 4 (2008) pp. 435-447
> -------- -------- --- Thank You --- -------- --------
> 
> 
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C.
> Berendsen
> GROMACS: Fast, Flexible and Free
> J. Comp. Chem. 26 (2005) pp. 1701-1719
> -------- -------- --- Thank You --- -------- --------
> 
> 
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> E. Lindahl and B. Hess and D. van der Spoel
> GROMACS 3.0: A package for molecular simulation and trajectory analysis
> J. Mol. Mod. 7 (2001) pp. 306-317
> -------- -------- --- Thank You --- -------- --------
> 
> 
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> H. J. C. Berendsen, D. van der Spoel and R. van Drunen
> GROMACS: A message-passing parallel molecular dynamics implementation
> Comp. Phys. Comm. 91 (1995) pp. 43-56
> -------- -------- --- Thank You --- -------- --------
> 
> 
> ++++ PLEASE CITE THE DOI FOR THIS VERSION OF GROMACS ++++
> https://doi.org/10.5281/zenodo.3923645
> -------- -------- --- Thank You --- -------- --------
> 
> Input Parameters:
>    integrator                     = md
>    tinit                          = 0
>    dt                             = 0.01
>    nsteps                         = 500000
>    init-step                      = 0
>    simulation-part                = 1
>    comm-mode                      = Angular
>    nstcomm                        = 100
>    bd-fric                        = 0
>    ld-seed                        = 1528453215
>    emtol                          = 10
>    emstep                         = 0.01
>    niter                          = 20
>    fcstep                         = 0
>    nstcgsteep                     = 1000
>    nbfgscorr                      = 10
>    rtpi                           = 0.05
>    nstxout                        = 0
>    nstvout                        = 0
>    nstfout                        = 0
>    nstlog                         = 1000
>    nstcalcenergy                  = 100
>    nstenergy                      = 100
>    nstxout-compressed             = 1000
>    compressed-x-precision         = 100
>    cutoff-scheme                  = Verlet
>    nstlist                        = 20
>    pbc                            = xyz
>    periodic-molecules             = false
>    verlet-buffer-tolerance        = 0.005
>    rlist                          = 1.1
>    coulombtype                    = Reaction-Field
>    coulomb-modifier               = Potential-shift
>    rcoulomb-switch                = 0
>    rcoulomb                       = 1.1
>    epsilon-r                      = 15
>    epsilon-rf                     = inf
>    vdw-type                       = Cut-off
>    vdw-modifier                   = Potential-shift
>    rvdw-switch                    = 0
>    rvdw                           = 1.1
>    DispCorr                       = No
>    table-extension                = 1
>    fourierspacing                 = 0.12
>    fourier-nx                     = 0
>    fourier-ny                     = 0
>    fourier-nz                     = 0
>    pme-order                      = 4
>    ewald-rtol                     = 1e-05
>    ewald-rtol-lj                  = 0.001
>    lj-pme-comb-rule               = Geometric
>    ewald-geometry                 = 0
>    epsilon-surface                = 0
>    tcoupl                         = Nose-Hoover
>    nsttcouple                     = 20
>    nh-chain-length                = 1
>    print-nose-hoover-chain-variables = false
>    pcoupl                         = No
>    pcoupltype                     = Isotropic
>    nstpcouple                     = -1
>    tau-p                          = 1
>    compressibility (3x3):
>       compressibility[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>       compressibility[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>       compressibility[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>    ref-p (3x3):
>       ref-p[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>       ref-p[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>       ref-p[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>    refcoord-scaling               = No
>    posres-com (3):
>       posres-com[0]= 0.00000e+00
>       posres-com[1]= 0.00000e+00
>       posres-com[2]= 0.00000e+00
>    posres-comB (3):
>       posres-comB[0]= 0.00000e+00
>       posres-comB[1]= 0.00000e+00
>       posres-comB[2]= 0.00000e+00
>    QMMM                           = false
>    QMconstraints                  = 0
>    QMMMscheme                     = 0
>    MMChargeScaleFactor            = 1
> qm-opts:
>    ngQM                           = 0
>    constraint-algorithm           = Lincs
>    continuation                   = false
>    Shake-SOR                      = false
>    shake-tol                      = 0.0001
>    lincs-order                    = 4
>    lincs-iter                     = 1
>    lincs-warnangle                = 30
>    nwall                          = 0
>    wall-type                      = 9-3
>    wall-r-linpot                  = -1
>    wall-atomtype[0]               = -1
>    wall-atomtype[1]               = -1
>    wall-density[0]                = 0
>    wall-density[1]                = 0
>    wall-ewald-zfac                = 3
>    pull                           = false
>    awh                            = false
>    rotation                       = false
>    interactiveMD                  = false
>    disre                          = No
>    disre-weighting                = Conservative
>    disre-mixed                    = false
>    dr-fc                          = 1000
>    dr-tau                         = 0
>    nstdisreout                    = 100
>    orire-fc                       = 0
>    orire-tau                      = 0
>    nstorireout                    = 100
>    free-energy                    = no
>    cos-acceleration               = 0
>    deform (3x3):
>       deform[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>       deform[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>       deform[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>    simulated-tempering            = false
>    swapcoords                     = no
>    userint1                       = 0
>    userint2                       = 0
>    userint3                       = 0
>    userint4                       = 0
>    userreal1                      = 0
>    userreal2                      = 0
>    userreal3                      = 0
>    userreal4                      = 0
>    applied-forces:
>      electric-field:
>        x:
>          E0                       = 0
>          omega                    = 0
>          t0                       = 0
>          sigma                    = 0
>        y:
>          E0                       = 0
>          omega                    = 0
>          t0                       = 0
>          sigma                    = 0
>        z:
>          E0                       = 0
>          omega                    = 0
>          t0                       = 0
>          sigma                    = 0
>      density-guided-simulation:
>        active                     = false
>        group                      = protein
>        similarity-measure         = inner-product
>        atom-spreading-weight      = unity
>        force-constant             = 1e+09
>        gaussian-transform-spreading-width = 0.2
>        gaussian-transform-spreading-range-in-multiples-of-width = 4
>        reference-density-filename = reference.mrc
>        nst                        = 1
>        normalize-densities        = true
>        adaptive-force-scaling     = false
>        adaptive-force-scaling-time-constant = 4
> grpopts:
>    nrdf:        1794
>    ref-t:         200
>    tau-t:           4
> annealing:      Single
> annealing-npoints:           2
> annealing-time [0]:	         0.0      5000.0
> annealing-temp [0]:	      1200.0       200.0
>    acc:	           0           0           0
>    nfreeze:           N           N           N
>    energygrp-flags[  0]: 0
> 
> Changing nstlist from 20 to 100, rlist from 1.1 to 1.174
> 
> Using 1 MPI thread
> Using 4 OpenMP threads 
> 
> System total charge: 0.000
> Reaction-Field:
> epsRF = 0, rc = 1.1, krf = 0.375657, crf = 1.36364, epsfac = 9.26236
> The electrostatics potential has its minimum at r = 1.1
> Potential shift: LJ r^-12: -3.186e-01 r^-6: -5.645e-01
> 
> Using SIMD 4x8 nonbonded short-range kernels
> 
> Using a dual 4x8 pair-list setup updated with dynamic pruning:
>   outer list: updated every 100 steps, buffer 0.074 nm, rlist 1.174 nm
>   inner list: updated every  67 steps, buffer 0.001 nm, rlist 1.101 nm
> At tolerance 0.005 kJ/mol/ps per atom, equivalent classical 1x1 list would be:
>   outer list: updated every 100 steps, buffer 0.247 nm, rlist 1.347 nm
>   inner list: updated every  67 steps, buffer 0.118 nm, rlist 1.218 nm
> 
> Using full Lennard-Jones parameter combination matrix
> 
> Removing pbc first time
> 
> Initializing LINear Constraint Solver
> 
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> B. Hess
> P-LINCS: A Parallel Linear Constraint Solver for molecular simulation
> J. Chem. Theory Comput. 4 (2008) pp. 116-122
> -------- -------- --- Thank You --- -------- --------
> 
> The number of constraints is 600
> 600 constraints are involved in constraint triangles,
> will apply an additional matrix expansion of order 4 for couplings
> between constraints inside triangles
> There are: 800 Atoms
> 
> Constraining the starting coordinates (step 0)
> 
> Constraining the coordinates at t0-dt (step 0)
> Center of mass motion removal mode is Angular
> We have the following groups for center of mass motion removal:
>   0:  System
> Group System with mass  3.60000e+04, Ekrot  9.43645e+00 Det(I) =  5.62837e+13
>   COM:      5.00001       5.00038       4.99989
>   P:        0.00011      -0.00009       0.00003
>   V:        0.00000      -0.00000       0.00000
>   J:     -233.79065    1014.41174     592.03381
>   w:       -0.00054       0.01368       0.00823
> Inertia tensor (3x3):
>    Inertia tensor[    0]={ 3.06927e+04,  6.06788e+03,  1.22821e+04}
>    Inertia tensor[    1]={ 6.06788e+03,  4.65128e+04,  3.63906e+03}
>    Inertia tensor[    2]={ 1.22821e+04,  3.63906e+03,  4.54163e+04}
> RMS relative constraint deviation after constraining: 5.71e-06
> Initial temperature: 1193.78 K
> 
> Started mdrun on rank 0 Wed Sep 16 12:33:26 2020
> 
>            Step           Time
>               0        0.00000
> 
> Current ref_t for group System:   1200.0
> 
>    Energies (kJ/mol)
>            Bond          Angle        LJ (SR)   Coulomb (SR)      Potential
>     2.00229e+03    4.86993e+03   -6.47130e+03    0.00000e+00    4.00924e+02
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     9.00372e+03    9.40464e+03    9.40464e+03    1.20724e+03    1.98238e+00
>    Constr. rmsd
>     5.72779e-06
> 
>            Step           Time
>            1000       10.00000
> 
> Current ref_t for group System:   1198.0
> 
>    Energies (kJ/mol)
>            Bond          Angle        LJ (SR)   Coulomb (SR)      Potential
>     2.02806e+03    4.82889e+03   -5.59476e+03    0.00000e+00    1.26219e+03
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     9.14918e+03    1.04114e+04    9.18513e+03    1.22675e+03    4.56846e+01
>    Constr. rmsd
>     5.27939e-06
> 
>            Step           Time
>            2000       20.00000
> 
> Current ref_t for group System:   1196.0
> 
>    Energies (kJ/mol)
>            Bond          Angle        LJ (SR)   Coulomb (SR)      Potential
>     2.19078e+03    5.03808e+03   -5.36282e+03    0.00000e+00    1.86604e+03
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     8.97587e+03    1.08419e+04    9.14384e+03    1.20351e+03   -2.86821e+00
>    Constr. rmsd
>     4.18652e-06
> 
>            Step           Time
>            3000       30.00000
> 
> Current ref_t for group System:   1194.0
> 
>    Energies (kJ/mol)
>            Bond          Angle        LJ (SR)   Coulomb (SR)      Potential
>     2.10494e+03    4.69843e+03   -5.79579e+03    0.00000e+00    1.00759e+03
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     9.31693e+03    1.03245e+04    8.84918e+03    1.24924e+03   -6.22543e+00
>    Constr. rmsd
>     5.25796e-06
> 
>            Step           Time
>            4000       40.00000
> 
> Current ref_t for group System:   1192.0
> 
>    Energies (kJ/mol)
>            Bond          Angle        LJ (SR)   Coulomb (SR)      Potential
>     2.05014e+03    5.10215e+03   -5.73250e+03    0.00000e+00    1.41979e+03
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     9.04671e+03    1.04665e+04    8.60389e+03    1.21301e+03    6.16145e+00
>    Constr. rmsd
>     6.82307e-06
> 
>            Step           Time
>            5000       50.00000
> 
> Current ref_t for group System:   1190.0
> 
>    Energies (kJ/mol)
>            Bond          Angle        LJ (SR)   Coulomb (SR)      Potential
>     1.81029e+03    4.76037e+03   -5.85207e+03    0.00000e+00    7.18587e+02
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     9.06558e+03    9.78416e+03    8.33397e+03    1.21554e+03    2.72711e+01
>    Constr. rmsd
>     6.04991e-06
> 
>            Step           Time
>            6000       60.00000
> 
> Current ref_t for group System:   1188.0
> 
>    Energies (kJ/mol)
>            Bond          Angle        LJ (SR)   Coulomb (SR)      Potential
>     2.09104e+03    5.04229e+03   -5.67435e+03    0.00000e+00    1.45898e+03
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     8.92920e+03    1.03882e+04    8.08663e+03    1.19725e+03   -7.92521e+00
>    Constr. rmsd
>     4.78571e-06
> 
> Constraint error in algorithm Lincs at step 6900
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6902
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6903
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6904
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6905
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6906
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6907
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6908
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6909
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6910
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6911
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6912
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6913
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6914
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6915
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6916
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6917
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6918
> Wrote pdb files with previous and current coordinates
> Constraint error in algorithm Lincs at step 6919
> Wrote pdb files with previous and current coordinates
> 
> -------------------------------------------------------
> Program:     gmx mdrun, version 2020.3
> Source file: src/gromacs/mdlib/constr.cpp (line 224)
> 
> Fatal error:
> Too many LINCS warnings (1114)
> If you know what you are doing you can adjust the lincs warning threshold in
> your mdp file
> or set the environment variable GMX_MAXCONSTRWARN to -1,
> but normally it is better to fix the problem
> 
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors
> -------------------------------------------------------

Hi,

there could be different reason for the lincs error. It is always good to visualize the step where it occurs, it can help to have an idea why occurs.
What do you constrain? Are you running atomistic or CG model?

Best regards
Alessandra

Hi there,

I believe the only thing being constrained is the bond angles. I’m using the Martini CG model.

Trying to keep angles rigid is very fragile and prone to failure. Warming the system may work because the system gradually acclimates to the high temperature, but the opposite direction is pretty likely to fail due to simple physical instability at extreme temperature, for which the force field was likely not designed nor validated.

I see. I performed the same thing but with several simulations of 4 PS-100-mers instead of several simulations of 2 PS-100-mers and there were no issues. It’s odd how it only failed with 2 PS-100mers. Maybe it’s due to the smaller system size?

Not likely size directly but with fewer molecules, you’re less likely to have a distorted geometry, though these could pop up randomly.