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

"""
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Reply via email to