[gmx-users] Switching from v4.0.7 to v 4.5.3 - being able to get the correct terms from a edr file when continuing a simulation

[Anna Duncan]

Hi,

I'm doing a coarse-grained simulation, using the MARTINI forcefield,  
of a protein in a lipid bilayer.  I carried out the equilibration  
stages using Gromacs 4.0.7.  Equilibration was done in several stages  
but the last stage was with a Nosé-Hoover thermostat for temperature  
coupling and Parrinello-Rahman thermostat for pressure coupling.  The  
next stage I carried out was a simulation of system with nothing  
except the backbone of the protein restrained, to relax the  
sidechains.  For this stage I used a newer, faster server which has  
Gromacs 4.5.3 installed on it.  The mdp file is below.
I used the command:
g_grompp -f ../md_T296.mdp -p ../prot_memb_system2.top -c npt_bPR.gro - 
e npt_bPR.edr -t npt_bPR.trr -n prot_bilayer.ndx -o md_schain

I got the error message:
-------------------------------------------------------
Program g_grompp, VERSION 4.5.3
Source code file: /builddir/build/BUILD/gromacs-4.5.3/src/gmxlib/ 
enxio.c, line: 1056

Fatal error:
Could not find energy term named 'Xi-0-Protein'
For more information and tips for troubleshooting, please check the  
GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------


I assumed that this was something to do with my switching between  
versions, maybe that the .edr file that's written by version 4.0.7  
doesn't have the terms that are required by version 4.5.3, so I  
removed the -e option from the command line and ran the simulation  
using the commands:
g_grompp -f ../md_T296.mdp -p ../prot_memb_system2.top -c npt_bPR.gro - 
t npt_bPR.trr -n prot_bilayer.ndx -o md_schain
g_mdrun -v -nt 1 -deffnm md_schain

For the next stage, I want to simulate the whole system, for which I  
use the mdp file below except with no -DPOSREBB defined (line 10) and  
with nsteps = 10000000, and version 4.5.3 again.
When I use the command:
g_grompp -f ../md_T296.mdp -p ../prot_memb_system2.top -c  
md_schain.gro -e md_schain.edr -t md_schain.trr -n prot_bilayer.ndx -o  
md1_t296

I again get the error message:
-------------------------------------------------------
Program g_grompp, VERSION 4.5.3
Source code file: /builddir/build/BUILD/gromacs-4.5.3/src/gmxlib/ 
enxio.c, line: 1056

Fatal error:
Could not find energy term named 'Xi-0-Protein'
For more information and tips for troubleshooting, please check the  
GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------


Can you help me with exactly what's going wrong?  I'm not sure if I  
should be carrying out the simulation without being able to pass on  
the information from the .edr file of the sidechain-relaxing  
simulation?  If not, what can I do to pass on the information  
necessary from the .edr file?

Many thanks in advance,

Anna

***********************************************************
mdp file
***********************************************************
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.0
;
; for use with GROMACS 3.3
;

; VARIOUS PREPROCESSING OPTIONS =
title                    = Martini
cpp                      = /usr/bin/cpp
define			 = -DPOSREBB

; RUN CONTROL PARAMETERS =
; MARTINI - Most simulations are stable with dt=40 fs,
; some (especially rings) require 20-30 fs.
; The range of time steps used for parametrization
; is 20-40 fs, using smaller time steps is therefore not recommended.

integrator               = md
; start time and timestep in ps
tinit                    = 0.0
dt                       = 0.030
nsteps                   = 1000000
; number of steps for center of mass motion removal =
nstcomm                  = 1
comm-mode		 = Linear
comm-grps		 = Protein_Lipids W

; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout                  = 5000
nstvout                  = 5000
nstfout                  = 0
; Output frequency for energies to log file and energy file =
nstlog                   = 1000
nstenergy                = 1000
; Output frequency and precision for xtc file =
nstxtcout                = 1000
xtc_precision            = 100
; This selects the subset of atoms for the xtc file. You can =
; select multiple groups. By default all atoms will be written. =
xtc-grps                 =
; Selection of energy groups =
energygrps               =

; NEIGHBORSEARCHING PARAMETERS =
; MARTINI - no need for more frequent updates
; or larger neighborlist cut-off due
; to the use of shifted potential energy functions.

; nblist update frequency =
nstlist                  = 10
; ns algorithm (simple or grid) =
ns_type                  = grid
; Periodic boundary conditions: xyz or none =
pbc                      = xyz
; nblist cut-off         =
rlist                    = 1.2

; OPTIONS FOR ELECTROSTATICS AND VDW =
; MARTINI - vdw and electrostatic interactions are used
; in their shifted forms. Changing to other types of
; electrostatics will affect the general performance of
; the model.

; Method for doing electrostatics =
coulombtype              = Shift
rcoulomb_switch          = 0.0
rcoulomb                 = 1.2
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r                = 15
; Method for doing Van der Waals =
vdw_type                 = Shift
; cut-off lengths        =
rvdw_switch              = 0.9
rvdw                     = 1.2
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr                 = No

; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; MARTINI - normal temperature and pressure coupling schemes
; can be used. It is recommended to couple individual groups
; in your system seperately.

; Temperature coupling   =
tcoupl                   = Nose-Hoover
nsttcouple		 = 1
; Groups to couple separately =
tc-grps                  = Protein W Lipids
; Time constant (ps) and reference temperature (K) =
tau_t                    = 1.2 1.2 1.2
ref_t                    = 296 296 296
; Pressure coupling      =
Pcoupl                   = Parrinello-Rahman
Pcoupltype               = semiisotropic
nstpcouple		 = 1
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau_p                    = 12.0
compressibility          = 3e-5 3e-5
ref_p                    = 1.0 1.0

; GENERATE VELOCITIES FOR STARTUP RUN =
gen_vel                  = no
gen_temp                 = 296
gen_seed                 = 666

; Distance restraints type: No, Simple or Ensemble
disre                    = Simple
; Force weighting of pairs in one distance restraint: Conservative or  
Equal
disre-weighting          = Conservative
; Use sqrt of the time averaged times the instantaneous violation
disre-mixed              = no
disre-fc                 = 500
disre-tau                = 0
; Output frequency for pair distances to energy file
nstdisreout              = 100

; OPTIONS FOR BONDS     =
; MARTINI - for ring systems constraints are defined
; which are best handled using Lincs.


constraints              = none
; Type of constraint algorithm =
constraint_algorithm     = Lincs
; Do not constrain the start configuration =
unconstrained_start      = yes
; Highest order in the expansion of the constraint coupling matrix =
lincs_order              = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs_warnangle          = 30

; FREEZE GROUPS
freezegrps		 =
freezedim		 =

; CONTINUATION
; this is specified by using -t -e flags in grompp and  
unconstrained_start = yes