[gmx-users] Crash in minimization that involves dummy atoms
Davide Bonanni
davide.bonanni at unito.it
Thu Jul 13 15:57:37 CEST 2017
Hi gromacs users,
I'm doing a relative free energy calculation of a ligand-protein complex.
I'm transforming a hydrogen atom into methyl. I have generated input
topology with FESetup. When I try to run minimization I get that error:
"""
Energy minimization has stopped, but the forces have not converged to the
Requested precision Fmax It stopped because the algorithm tried to make a
new step whose size was too
Small, or there was no change in the energy since last step. Either way, we
Regard the minimization as converged to within the available machine
Precision, given your starting configuration and EM parameters.
Double precision normally gives you more accuracy, but this is often not
Needed for preparing to run molecular dynamics.
You may need to increase your constraint accuracy, or turn
Off constraints altogether (set constraints = none in mdp file)
Steepest Descents converged to machine precision in 50 steps,
But did not reach the requested Fmax <100.
Potential Energy = 1.1698845e + 07
Maximum force = 7.4549694e + 05 is atom 5745
Norm of force = 3.5688174e + 03
"""
As you can see the Potential Energy is too high and if a look at the system
I got is completely messed up. Exactly the same complex with the same
parameters but a different perturbation (H -> Cl) gave me no problems. So
the minimization crash should be related to the presence of dummy atoms.
I read some discussion about the same topic but I was not able to find a
solution. Everything looks good to me, I do not see anything wrong in the
file .gro or .top.
I have tried also to run the same minimization on T4 lysozyme tutorial
input (https://ccpforge.cse.rl.ac.uk/gf/project/ccpbiosim/wiki/?
pagename=T4+lysozyme), and when there are dummy atom involved I obtain the
same problem.
Probably I'm missing something in the .mdp or there is some error I can not
see.
I tried to modify constraint but nothing changed.
Any kind of suggestion is really appreciated, thank you in advance.
Cheers,
Davide Bonanni
This is the first part of ligand's topology:
"""
[ moleculetype ]
LIG 3
[ atoms ]
; nr type resno resnm atom cgnr charge mass typeB
chargeB massB
1 c 1 LIG C1 1 0.667751 12.0100 c
0.667700 12.0100
2 cc 1 LIG C2 2 -0.469349 12.0100 cc
-0.471400 12.0100
3 ca 1 LIG C3 3 0.034751 12.0100 ca
0.039700 12.0100
4 ca 1 LIG C4 4 -0.083249 12.0100 ca
-0.032600 12.0100
5 c 1 LIG C5 5 0.742351 12.0100 c
0.742300 12.0100
6 ca 1 LIG C6 6 0.121651 12.0100 ca
0.122600 12.0100
7 ca 1 LIG C7 7 0.029951 12.0100 ca
0.029400 12.0100
8 ca 1 LIG C8 8 0.029951 12.0100 ca
0.029400 12.0100
9 ca 1 LIG C9 9 0.134451 12.0100 ca
0.134400 12.0100
10 ca 1 LIG C10 10 0.134451 12.0100 ca
0.134400 12.0100
11 cp 1 LIG C11 11 -0.137949 12.0100 cp
-0.139000 12.0100
12 cp 1 LIG C12 12 -0.012949 12.0100 cp
-0.013000 12.0100
13 ca 1 LIG C13 13 -0.108449 12.0100 ca
-0.109000 12.0100
14 ca 1 LIG C14 14 -0.108449 12.0100 ca
-0.109000 12.0100
15 ca 1 LIG C15 15 -0.136949 12.0100 ca
-0.137000 12.0100
16 ca 1 LIG C16 16 -0.136949 12.0100 ca
-0.137000 12.0100
17 ca 1 LIG C17 17 -0.134949 12.0100 ca
-0.135000 12.0100
18 ca 1 LIG C18 18 -0.155949 12.0100 ca
-0.158000 12.0100
19 ca 1 LIG C19 19 -0.122949 12.0100 ca
-0.122000 12.0100
20 ca 1 LIG C20 20 -0.226949 12.0100 ca
-0.223000 12.0100
21 f 1 LIG F21 21 -0.109349 19.0000 f
-0.109400 19.0000
22 f 1 LIG F22 22 -0.130849 19.0000 f
-0.131400 19.0000
23 f 1 LIG F23 23 -0.130849 19.0000 f
-0.131400 19.0000
24 f 1 LIG F24 24 -0.109349 19.0000 f
-0.109400 19.0000
25 nc 1 LIG N25 25 -0.632449 14.0100 nc
-0.632500 14.0100
26 na 1 LIG N26 26 0.281751 14.0100 na
0.281700 14.0100
27 n 1 LIG N27 27 -0.465049 14.0100 n
-0.464100 14.0100
28 o 1 LIG O28 28 -0.704049 16.0000 o
-0.705100 16.0000
29 o 1 LIG O29 29 -0.651049 16.0000 o
-0.652100 16.0000
30 h4 1 LIG H30 30 0.158051 1.0080 c3
-0.035800 12.0100
31 ha 1 LIG H31 31 0.141051 1.0080 ha
0.141000 1.0080
32 ha 1 LIG H32 32 0.141051 1.0080 ha
0.141000 1.0080
33 ha 1 LIG H33 33 0.125051 1.0080 ha
0.125000 1.0080
34 ha 1 LIG H34 34 0.125051 1.0080 ha
0.125000 1.0080
35 ha 1 LIG H35 35 0.121051 1.0080 ha
0.121000 1.0080
36 ha 1 LIG H36 36 0.167051 1.0080 ha
0.168000 1.0080
37 ha 1 LIG H37 37 0.117051 1.0080 ha
0.117000 1.0080
38 ha 1 LIG H38 38 0.124051 1.0080 ha
0.124000 1.0080
39 hn 1 LIG H39 39 0.371551 1.0080 hn
0.371500 1.0080
40 du 1 LIG DU40 40 0.000000 1.0000 hc
0.047367 1.0080
41 du 1 LIG DU41 41 0.000000 1.0000 hc
0.047367 1.0080
42 du 1 LIG DU42 42 0.000000 1.0000 hc
0.047367 1.0080
...
...
...
"""
This is the em.mdp:
"""
; Run control
integrator = steep
nsteps = 5000
; EM criteria and other stuff
emtol = 100
emstep = 0.01
niter = 20
nbfgscorr = 10
; Output control
nstlog = 1
nstenergy = 1
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 1
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-5
epsilon_surface = 0
; Temperature and pressure coupling are off during EM
tcoupl = no
pcoupl = no
; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for
each simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
40
fep-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.10 0.20
0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
mass-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
nstdhdl = 100
calc-lambda-neighbors = -1
sc-alpha = 0.5
sc-coul = yes
sc-power = 1
sc-r-power = 6
sc-sigma = 0.3
dhdl-derivatives = yes
dhdl-print-energy = no
separate-dhdl-file = yes
dh_hist_size = 0
dh_hist_spacing = 0.1
; No velocities during EM
gen_vel = no
; options for bonds
constraints = all-bonds
; Type of constraint algorithm
constraint-algorithm = lincs
; Do not constrain the starting configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
lincs-warnangle = 30 ; maximum angle that a bond can rotate
energygrps = Protein Non-Protein
"""
More information about the gromacs.org_gmx-users
mailing list