[gmx-users] ligan minimization in vacuo

[Justin Lemkul]

On 7/10/17 5:57 PM, ‪farial tavakoli‬ ‪ wrote:
> Dear gmx-users
> I have a problem in equilibration my protein-ligand complex and encountered to this error after 2 steps of 50000 steps:one or more water molecules can not be settled. check for bad contacts or reduce the time steps.
> so I decided to create a topology for my designed drug (50 atoms) as a ligand to minimize it in the vacuo in the absence of the protein to check if it minimizes  , so i tried to change the .itp file obtained from PRODRG to create a .top file,  in this way:I added  " ; Include forcefield parameters
>      #include "gromos43a1.ff/forcefield.itp " above the [ moleculetype ] directive as a first line , then added [ system ] and [ molecule ] directives under the ligand, respectively :
>    23  24  26  27   1    180.0   33.5 2    180.0   33.5 2 ; dih   CAS  CAR  NAQ  CAN
>    32  27  26  24   1    180.0   33.5 2    180.0   33.5 2 ; dih   CAM  CAN  NAQ  CAR
>    39  37  36  30   1    180.0    5.9 2    180.0    5.9 2 ; dih   SAG  CAH  CAK  CAP
>    40  48  49  50   1    180.0    7.1 2    180.0    7.1 2 ; dih   CAA  CAF  OAJ  HAJ
> 
> ; Ligand position restraints
> #ifdef POSRES
> #include "posre_DRG.itp"
> #endif
> 
> ; Include water topology
> #include "gromos43a1.ff/spc.itp"
> 
> [ system ]
> ; Name
> DRG in water
> 
> [ molecules ]
> ; Compound        #mols
> DRG                      1
> issued this command:gmx grompp -v -f em.mdp -c DRG.gro -p DRG.top -o DRG-EM-vacuum.tprgmx mdrun -v -deffnm DRG-EM-vacuum -c DRG-EM-vacuum.gro
> and simulation results:
> Steepest Descents:
>     Tolerance (Fmax)   =  1.00000e+03
>     Number of steps    =        50000
> Step=    0, Dmax= 1.0e-02 nm, Epot=  1.41636e+03 Fmax= 2.81150e+04, atom= 25
> Step=    1, Dmax= 1.0e-02 nm, Epot=  1.01416e+03 Fmax= 1.21118e+04, atom= 25
> Step=    2, Dmax= 1.2e-02 nm, Epot=  7.87016e+02 Fmax= 4.71639e+03, atom= 25
> Step=    3, Dmax= 1.4e-02 nm, Epot=  6.94638e+02 Fmax= 5.88719e+03, atom= 15
> Step=    5, Dmax= 8.6e-03 nm, Epot=  6.79808e+02 Fmax= 7.13390e+03, atom= 15
> Step=    7, Dmax= 5.2e-03 nm, Epot=  6.29169e+02 Fmax= 1.31714e+03, atom= 25
> Step=    8, Dmax= 6.2e-03 nm, Epot=  6.12279e+02 Fmax= 2.91712e+03, atom= 14
> Step=   10, Dmax= 3.7e-03 nm, Epot=  6.04849e+02 Fmax= 2.77294e+03, atom= 27
> Step=   12, Dmax= 2.2e-03 nm, Epot=  5.86062e+02 Fmax= 9.77798e+02, atom= 24
> 
> writing lowest energy coordinates.
> 
> Back Off! I just backed up DRG-EM-vacuum.gro to ./#DRG-EM-vacuum.gro.1#
> 
> Steepest Descents converged to Fmax < 1000 in 13 steps
> Potential Energy  =  5.8606232e+02
> Maximum force     =  9.7779791e+02 on atom 24
> Norm of force     =  3.8248090e+02
> 
> Simulation ended prematurely, no performance report will be written.
> 
> 
>   here's my em.mdp file:
> 
> integrator      = steep        ; Algorithm (steep = steepest descent minimization)
> emtol            = 1000.0      ; Stop minimization when the maximum force < 10.0 kJ/mol
> emstep         = 0.01         ; Energy step size
> nsteps          = 50000       ; Maximum number of (minimization) steps to perform
> energygrps    = system    ; Which energy group(s) to write to disk
> 
> ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
> nstlist               = 1            ; Frequency to update the neighbor list and long range forces
> cutoff-scheme   = Verlet
> ns_type            = grid        ; Method to determine neighbor list (simple, grid)
> rlist                  = 1.0        ; Cut-off for making neighbor list (short range forces)
> coulombtype     = PME      ; Treatment of long range electrostatic interactions
> rcoulomb          = 1.0        ; long range electrostatic cut-off
> rvdw                = 1.0        ; long range Van der Waals cut-off
> pbc                 = xyz        ; Periodic Boundary Conditions
> 

To minimize in vacuo, cutoffs should all be set to zero, no PBC, no PME, simple 
neighbor searching and group cutoff scheme.

> I just wanted to know if I did correct?In addition after minimization in vacuo , checked my ligand in pymol and noticed it is disintegrated.  Is there anyone help me how come it is disintegrated ,however i edited the .itp file obtained from PRODRG ( charges and charg groups)? How can i fix this problem to equilibrate my complex?

You'll have to define "disintegrated" because bonds can't break or form in 
classical molecular mechanical processes, so perhaps PyMOL is just having a hard 
time guessing the bonded structure.  You have to be judicious in what 
modifications you make to the topology.  PRODRG topologies require complete 
reparametrization of (at least) the charges, but without knowing what you did or 
why you did it, there's little to go on here.

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

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