Greetings Grimaces community,
I am trying to execute the following “The first step is 1000 steps of SD minimization with
strong positional restraints applied to the heavy (i.e., non- hydrogen) atoms of the large molecules using a force constant of 5.0 kcal/mol Å and the initial coordinates as a reference. No other constraints (e.g., SHAKE) should be applied during this step”
I am using a typical em.mdp file. H0ow do I add restriction for the heavy atoms?
Here is my em.mdp file
; LINES STARTING WITH ‘;’ ARE COMMENTS
title = Minimization ; Step#1 J. Chem. Phys. 153, 054123 (2020)
; Parameters describing what to do, when to stop and what to save
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 = 1000; Maximum number of (minimization)
; 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.2 ; Cut-off for making neighbor list (short range forces)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.2 ; long range electrostatic cut-off
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; long range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions
DispCorr = no
There is no general rule for how to apply position restraints. They are (usually) defined in the force field parameters for your system. That said, if you have followed usual system setup procedures (gmx pdb2gmx etc) you should be able be able to turn on restraints using define = -DPOSRES in the MDP file. The exact strengths of the restraints, and which atoms are restrained, are usually found in one or more posre_[…].itp files where you generated your topology.
Gromacs Community,
Thanks for your response. Is this procedure OK
(1) Add define = -DPOSRES to my em.mdp file
(2) Create position restraint file for the protein with entires as
[ position_restraints ]
; atom type fx fy fz
1 1 5 5 5
4 1 5 5 5
(3) Update the tool.top file to include only the proteins that will be position restraint
; Include Position restraint file #ifdef POSRES #include “posre_005.itp” #endif
I think you’re almost there. I think it’s only the units of the restraints you need to consider. GROMACS uses kJ/mol nm^2. You wrote above that you wanted a restraint of 5 kcal/mol Å (I presume kcal/mol Å^2), which would translate to 2092 kJ/mol nm^2.
The publication that I am following “A protocol for preparing explicitly solvated systems for stable molecular dynamics simulations, J. Chem. Phys. 153, 054123 (2020)” - uses
kcal/mol Å throughout as the force constant unit. I hope this is an oversight!
As they also state that they are using harmonic positional restraints, I presume it’s a mistake on their part, yes. But you could contact the authors if you want to be sure.