Shake algorithm not converging

I am trying to use the shake algorithm with the all-angle constraint option. My system is crashing and I receive a message:

“Shake did not converge in 1000 steps”

I am aware that this error has been reported online, but have not been able to find a solution in the discussions that helps me.

For the simplest system where the error occurs, I am simulating a mixture of water (have tried many different models) and a one glycerol molecule (.itp and .gro files included below)

The glycerol molecule is taken from the PDB and I (and others) have used it successfully in many contexts. So I do not think it is a poor configuration.

I have tried to reduce time stepping and have performed extensive energy minimisation.

Varying the shake-tol does have an effect. The system always crashes for default shake-tol = 0.0001. For values of shake-tol = 0.005 then the systems runs. However for my specific use, (which involves additional molecules) I need the tolerance to be closer to the default value.

Should I Increase the number of shake iterations? I would think that this is not the solution. Actually, I am not sure how to change this variable, I could not find this information either.

Interesting is the following from the .log file for the crashing run:

Constraining the starting coordinates (step 0)
Shake did not converge in 1000 steps
i mi j mj before after should be
1 1.01 2 16.00 0.10000 0.10000 0.10000
1 1.01 3 12.01 0.19929 0.19929 0.19928
2 16.00 3 12.01 0.14345 0.14345 0.14350
2 16.00 4 1.01 0.20815 0.20815 0.20804
2 16.00 5 1.01 0.20485 0.20485 0.20474
2 16.00 6 12.01 0.25013 0.25013 0.25012
3 12.01 4 1.01 0.10980 0.10980 0.11000
3 12.01 5 1.01 0.10882 0.10882 0.10900
3 12.01 6 12.01 0.15302 0.15302 0.15300
3 12.01 7 1.01 0.21420 0.21420 0.21427
3 12.01 8 16.00 0.24219 0.24219 0.24220
3 12.01 10 12.01 0.25217 0.25217 0.25218
4 1.01 5 1.01 0.17708 0.17708 0.17718
4 1.01 6 12.01 0.21622 0.21622 0.21621
5 1.01 6 12.01 0.21093 0.21093 0.21091
6 12.01 7 1.01 0.11012 0.11012 0.11000
6 12.01 8 16.00 0.14351 0.14351 0.14350
6 12.01 9 1.01 0.19928 0.19928 0.19928
6 12.01 10 12.01 0.15289 0.15289 0.15300
6 12.01 11 1.01 0.21451 0.21451 0.21427
6 12.01 12 1.01 0.21649 0.21649 0.21621
6 12.01 13 16.00 0.24443 0.24443 0.24441
7 1.01 8 16.00 0.20794 0.20794 0.20792
7 1.01 10 12.01 0.21615 0.21615 0.21621
8 16.00 9 1.01 0.10000 0.10000 0.10000
8 16.00 10 12.01 0.24213 0.24213 0.24220
10 12.01 11 1.01 0.10950 0.10950 0.11000
10 12.01 12 1.01 0.10951 0.10951 0.11000
10 12.01 13 16.00 0.14347 0.14347 0.14350
10 12.01 14 1.01 0.19928 0.19928 0.19928
11 1.01 12 1.01 0.17773 0.17773 0.17798
11 1.01 13 16.00 0.20551 0.20551 0.20552
12 1.01 13 16.00 0.20891 0.20891 0.20891
13 16.00 14 1.01 0.09999 0.09999 0.10000
Constraint error in algorithm Shake at step 0
Wrote pdb files with previous and current coordinates

From my interpretation of the “before” and “after” columns, it does not appear that any of the constraints are being update (or evolved) during the iteration.

My basic run file looks like this (with an addition energy minimisation step):

gmx solvate -maxsol 1 -box 3.3 3.3 3.3 -cs glycerol.gro -o glycerol_box.gro -p topol.top
gmx solvate -maxsol 800 -cp glycerol_box.gro -cs spc216.gro -o glycerol_solvate.gro -p topol.top

gmx grompp -f em.mdp -c glycerol_solvate.gro -p topol.top -o en_min.tpr
gmx mdrun -deffnm en_min

gmx grompp -f npt.mdp -c en_min.gro -p topol.top -o npt.tpr -maxwarn 1
gmx mdrun -deffnm npt

gmx grompp -f em.mdp -c npt.gro -p topol.top -o en_min2.tpr -maxwarn 1
gmx mdrun -deffnm en_min2

gmx grompp -f nvt.mdp -c en_min2.gro -p topol.top -o nvt.tpr -maxwarn 1
gmx mdrun -deffnm nvt

The nvt.mdp is:

; Running Parameters ===============================================
integrator = md ; leap-frog integrator
nsteps = 25000 ; 1000000 = 1 microsecond
dt = 0.0002 ; fs

; Output control ===================================================
nstxout-compressed = 1
xtc-precision = 1000000
nstlog = 250000

; Constraints ======================================================
constraint-algorithm = SHAKE
; shake-tol = 0.005
constraints = all-angles

; Neighborsearching ===============================================
cutoff-scheme = Verlet
ns_type = grid
rlist = 1.041
nstlist = 100 ; neighbour list calculate

; VdW and electrostatics ===========================================
rvdw = 1.0
rvdw-switch = 1.0
rcoulomb = 1.0
rcoulomb-switch = 0
periodic_molecules = no
coulombtype = pme
vdw-modifier = Potential-shift

; velocity generation
gen_vel = yes

; Dispersion correction
DispCorr = EnerPres

; Temperature Coupling =============================================
tcoupl = v-rescale ; modified Berendsen thermostat
ld_seed = -1
tau_t = 1.0
ref_t = 300.0
tc_grps = System

; Pressure Coupling ================================================
Pcoupl = no
;compressibility = 4.5E-5
;tau_p = 1.0
;ref_p = 1.01325
;Pcoupltype = isotropic

; FFT grid size, when a value is 0 fourierspacing will be used
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0

; EWALD/PME/PPPM parameters =====================================
pme_order = 4
ewald_rtol = 1e-04
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no

The glycerol.itp is:

[ moleculetype ]
; Name nrexcl
6KSG 3
[ atoms ]
; nr type resnr resid atom cgnr charge mass total_charge
1 H 1 6KSG H3 1 0.394 1.0080
2 OA 1 6KSG O1 1 -0.625 15.9994
3 C 1 6KSG C1 1 0.164 12.0110
4 HC 1 6KSG H1 1 0.038 1.0080
5 HC 1 6KSG H2 1 0.038 1.0080 ; 0.009
6 C 1 6KSG C2 2 0.092 12.0110
7 HC 1 6KSG H4 2 0.064 1.0080
8 OA 1 6KSG O2 2 -0.597 15.9994
9 H 1 6KSG H5 2 0.423 1.0080 ; -0.018
10 C 1 6KSG C3 3 0.164 12.0110
11 HC 1 6KSG H6 3 0.038 1.0080
12 HC 1 6KSG H7 3 0.038 1.0080
13 OA 1 6KSG O3 3 -0.625 15.9994
14 H 1 6KSG H8 3 0.394 1.0080 ; 0.009
; total charge of the molecule: 0.000
[ bonds ]
; ai aj funct c0 c1
1 2 2 0.1000 1.5700e+07
2 3 2 0.1435 6.1000e+06
3 4 2 0.1100 1.2100e+07
3 5 2 0.1090 1.2300e+07
3 6 2 0.1530 7.1500e+06
6 7 2 0.1100 1.2100e+07
6 8 2 0.1435 6.1000e+06
6 10 2 0.1530 7.1500e+06
8 9 2 0.1000 1.5700e+07
10 11 2 0.1100 1.2100e+07
10 12 2 0.1100 1.2100e+07
10 13 2 0.1435 6.1000e+06
13 14 2 0.1000 1.5700e+07
[ pairs ]
; ai aj funct ; all 1-4 pairs but the ones excluded in GROMOS itp
1 4 1
1 5 1
1 6 1
2 7 1
2 8 1
2 10 1
3 9 1
3 11 1
3 12 1
3 13 1
4 7 1
4 8 1
4 10 1
5 7 1
5 8 1
5 10 1
6 14 1
7 9 1
7 11 1
7 12 1
7 13 1
8 11 1
8 12 1
8 13 1
9 10 1
11 14 1
12 14 1
[ angles ]
; ai aj ak funct angle fc
1 2 3 2 108.53 443.00
2 3 4 2 109.60 450.00
2 3 5 2 107.57 484.00
2 3 6 2 115.00 610.00
4 3 5 2 108.00 465.00
4 3 6 2 109.50 285.00
5 3 6 2 106.00 848.00
3 6 7 2 108.00 740.00
3 6 8 2 109.50 520.00
3 6 10 2 111.00 530.00
7 6 8 2 109.50 618.00
7 6 10 2 109.50 618.00
8 6 10 2 109.50 520.00
6 8 9 2 108.53 443.00
6 10 11 2 108.00 740.00
6 10 12 2 109.50 285.00
6 10 13 2 111.00 530.00
11 10 12 2 108.00 465.00
11 10 13 2 107.60 507.00
12 10 13 2 110.30 524.00
10 13 14 2 108.53 443.00
[ dihedrals ]
; GROMOS improper dihedrals
; ai aj ak al funct angle fc
[ dihedrals ]
; ai aj ak al funct ph0 cp mult
1 2 3 6 1 0.00 1.26 3
2 3 6 8 1 0.00 5.92 3
3 6 8 9 1 180.00 1.00 3
3 6 10 13 1 0.00 5.92 3
6 10 13 14 1 0.00 1.26 3
[ exclusions ]
; ai aj funct ; GROMOS 1-4 exclusions

glycerol.gro:

single glycerol
14
16KSG H3 1 0.320 0.225 0.096 0.3437 -0.7676 0.1009
16KSG O1 2 0.255 0.177 0.038 -0.3373 -0.2037 0.3774
16KSG C1 3 0.128 0.205 0.087 0.2463 0.1211 -0.1218
16KSG H1 4 0.055 0.171 0.012 0.4618 -0.5324 -0.0408
16KSG H2 5 0.110 0.311 0.103 -2.5220 -0.1869 -0.7309
16KSG C2 6 0.095 0.137 0.224 -0.0869 -0.0916 -0.2780
16KSG H4 7 -0.006 0.173 0.247 -0.4511 -0.1175 -1.6964
16KSG O2 8 0.095 -0.010 0.224 0.0493 -0.5914 0.5611
16KSG H5 9 0.101 -0.030 0.126 0.1925 -0.6527 0.5820
16KSG C3 10 0.196 0.174 0.333 -0.2795 -0.5942 -0.2071
16KSG H6 11 0.296 0.137 0.308 0.2189 -0.6429 1.7642
16KSG H7 12 0.167 0.130 0.429 -1.9416 0.7505 -0.0402
16KSG O3 13 0.218 0.318 0.350 0.6131 -0.2762 -0.1693
16KSG H8 14 0.255 0.313 0.443 -0.7638 1.2757 0.5177
0.32600 0.34800 0.43100

Any advice would be great. Thanks in advance for any help. Apologies, but I have no idea how to move forward.

Hi Alevilla,

Thanks for the reply. Actually, the manual page to which you link says that:

“SHAKE is slightly slower and less stable than LINCS, but does work with angle constraints”.

Plus, I have run multiple numerical tests on other molecules (butane, for example) for which I do not have the problems described above. That is, for which the SHAKE all-angle constraint works as expected for the input parameters that I have used.

Regards,
Ben

Sorry, you are right. Post remove. Thank you very much for reporting it.
Alessandra

Hi Alessandra,

Not a problem.

Do you have any other thoughts about what the problem could be? Sorry that the presentation of the input files is so unpleasant but as a new user I was not given permission to upload files. I would be great if you could provide any other insights. I could try send you the input files, if that helps.

Cheers,
Ben

Dear Dalton

I am a new user of the commander I write a diploma associated with relaxation in Glycerin. I could not find the PDB file.
Can you tell you from these file from?
Can you help in the preparation of topology files?