Position restraints are not applied

User discussions
1

GROMACS version:gmx, 2022.1
GROMACS modification: No

I am running an equilibration simulation (step6.3) using the Dry Martini force field in GROMACS, where I need to apply position restraints on the headgroups of DOPC and cholesterol in my vesicle system.

I have done the following to enforce position restraints:

  1. Created posre_dopc.itp and posre_chol.itp with force constants (1000 1000 1000).
  2. Included these files in system.top using #ifdef POSRES.
  3. Defined -DPOSRES in step6.3_equilibration.mdp.
  4. Checked that my index.ndx contains the correct groups (DOPC_HEAD, CHOL_HEAD).
  5. Ran grompp to generate step6.3_equilibration.tpr.

However, when I check the .tpr file using:
gmx dump -s step6.3_equilibration.tpr | grep posres

I see that position restraints are not applied (posres = 0.00000e+00 instead of 1000 1000 1000).

Despite setting everything correctly, position restraints are missing in the .tpr file.
Could you suggest how to debug and fix this issue?

Am I missing any step in applying position restraints in Dry Martini systems?

Thank you for your guidance!

2

Hi @ehsank

How are the constraints defined in the .itp files?

Also as far as I know the restraints are independent from whatever you define in the index, as they are read from the topology if the keyword (POSRES in your case) is defined.

3

Dear obZehn

I have been working on applying position restraints on the headgroups of DOPC and cholesterol in my vesicle system using the Dry Martini force field in GROMACS. Below is a summary of my approach and the encountered issue:

Steps Taken:

  1. Created Position Restraint Files:
  • I generated posre_dopc.itp and posre_chol.itp, where I applied force constants of 1000 kJ/mol/nm² on selected atoms extracted from dopc_posres_atoms.txt and chol_posres_atoms.txt.
  • Example entries:
[ position_restraints ]
; i funct       fcx        fcy        fcz
2 1 1000 1000 1000
14 1 1000 1000 1000
26 1 1000 1000 1000
...
  1. Modified system.top to Include Restraints:
  • Used #ifdef POSRES to include the .itp files only when -DPOSRES is defined in the .mdp file.
#ifdef POSRES
#include "posre_dopc.itp"
#include "posre_chol.itp"
#endif
  1. Ensured Correct Groups in index.ndx:
  • Used gmx make_ndx to generate groups DOPC_HEAD and CHOL_HEAD and verified them.
  1. Checked the .tpr File for Restraints:
  • After running:
gmx dump -s step6.3_equilibration.tpr | grep posres

the output showed posres = 0.00000e+00, indicating that the restraints were not applied.

Issue and Diagnosis:

  • Error Message:
ERROR 1 [file posre_dopc_fixed.itp, line 5]:
Atom index (14) in position_restraints out of bounds (1-12).
  • This suggests that the atom indices in posre_dopc.itp are incorrect because Dry Martini represents DOPC molecules using only 12 beads and CHOL molecules using 8 beads.
  • The indices in my .itp file were extracted from the full system, not per molecule, causing the index out-of-bounds error.
4

How generate itp file for CHOL and DOPC in dry martini when used corse grain?

5

The error you are getting is related to the fact that you have to define the restraints based on the molecule’s atoms within the molecule’s definition and not by referring to the whole box.

For example, cholesterol will be defined in its topology as a molecule with 8 beads, so it will have 8 atoms [1-8], each one with a type, charge, etc. Now you have to define the position restraints with respect to THESE numbers (so the individual molecules), and not the whole box. As such, if say atom 2 of cholesterol and atom 3 of DOPC are their (beadified) heads, then your position restraints will look simply like this

; restraints for cholesterol head
[ position_restraints ]
; i funct       fcx        fcy        fcz
 2       1       1000       1000       1000

and

; restraints for DOPC head
[ position_restraints ]
; i funct       fcx        fcy        fcz
 3       1       1000       1000       1000
6

thank you so much for your considerataion for clarification, it is my system.top

; This file was generated by TS Back Mapping

; Include the Martini force field and necessary lipid files
#include “dry_martini_v2.1.itp” ; Martini 2.1 force field
#include “dry_martini_v2.1_cholesterol.itp” ; Cholesterol specific parameters
#include “dry_martini_v2.1_lipids.itp” ; DOPC and other lipids

[ system ]
Expect a large membrane

[ molecules ]
; domain 0
; in the upper monolayer
DOPC 1389
#ifdef POSRES
#include “posre_dopc.itp”
#endif

CHOL 926
#ifdef POSRES
#include “posre_chol.itp”
#endif

; domain 0
; in the lower monolayer
DOPC 533
#ifdef POSRES
#include “posre_dopc.itp”
#endif

CHOL 355
#ifdef POSRES
#include “posre_chol.itp”
#endif

I got this error again :
Command line:
gmx grompp -f step6.3_equilibration.mdp -o step6.3_equilibration.tpr -c step6.2_equilibration.gro -p system.top -n index.ndx

NOTE 1 [file step6.3_equilibration.mdp]:
You have set rlist larger than the interaction cut-off, but you also have
verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.

NOTE 2 [file step6.3_equilibration.mdp]:
Setting tcoupl from ‘V-rescale’ to ‘no’. sd handles temperature coupling
implicitly. See the documentation for more information on which
parameters affect temperature for sd.

Setting the LD random seed to -1470997761

Generated 0 of the 703 non-bonded parameter combinations

Excluding 1 bonded neighbours molecule type ‘DOPC’

WARNING 1 [file system.top, line 19]:
Too few parameters on line (source file
/admin/build/admin/rpms/frontera/BUILD/gromacs-2022.1/src/gromacs/gmxpreprocess/toppush.cpp, line 1900)

WARNING 2 [file system.top, line 26]:
Too few parameters on line (source file
/admin/build/admin/rpms/frontera/BUILD/gromacs-2022.1/src/gromacs/gmxpreprocess/toppush.cpp, line 1900)

WARNING 3 [file system.top, line 31]:
Too few parameters on line (source file
/admin/build/admin/rpms/frontera/BUILD/gromacs-2022.1/src/gromacs/gmxpreprocess/toppush.cpp, line 1900)

Program: gmx grompp, version 2022.1
Source file: src/gromacs/gmxpreprocess/grompp.cpp (line 703)

Fatal error:
number of coordinates in coordinate file (step6.2_equilibration.gro, 33312)
does not match topology (system.top, 16668)

For more information and tips for troubleshooting, please check the GROMACS
website at Common Errors — GROMACS webpage https://www.gromacs.org documentation

login2.frontera(1017)$

my step6.3_equilibration.mdp is:
; Integrator and time settings
define = -DPOSRES ; This enables the #ifdef POSRES directive in system.top
integrator = sd ; Stochastic dynamics (SD) integrator for Dry Martini
tinit = 0.0
dt = 0.03 ; Time step of 30 fs, typical for Martini simulations
nsteps = 3000000 ; Run for 2 million steps (total time = 60 ns)

; Output control
nstlog = 1000
nstenergy = 1000
nstxout-compressed = 1000
compressed-x-precision = 100

; Neighbor searching settings
cutoff-scheme = Verlet
nstlist = 10 ; Update neighbor list every 10 steps
rlist = 1.4 ; Neighbor list cutoff at 1.4 nm
pbc = xyz
verlet-buffer-tolerance = 0.005

; Electrostatics and van der Waals settings
epsilon_r = 15 ; Relative permittivity for implicit solvent
coulombtype = reaction-field
rcoulomb = 1.2 ; Coulomb interactions cutoff at 1.2 nm
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet ; Smoothly shift van der Waals potential to zero
rvdw = 1.2 ; Van der Waals cutoff at 1.2 nm

; No pressure coupling (NVT ensemble)
Pcoupl = no ; No pressure coupling for implicit solvent

; Minimal T-coupling to satisfy GROMACS requirements
tcoupl = v-rescale
tc-grps = System
tau_t = 1.0 ; Minimal coupling constant
ref_t = 310 ; Target temperature for the dummy group

; Generate initial velocities if needed (only for equilibration)
gen_vel = yes
gen_temp = 310
gen_seed = 0395104856

refcoord_scaling = all

I want impose position retrain for my head drop of CHOL and DOPC in my vesicle because when I run equilibration my vesicle has 2 bumper on it , I try by this way remove them . I used dry martini altho

Best Regards
Ehsan

7

Where did you get this .mdp file? I am not familiar with martini dry, so I can help you only partially. Please also refer to the manual.

Firstly, you get two notes. These are not fundamental problems, but are still telling you that something is not fully okay and GROMACS is handling the problem internally.

Note 1 tells you that you are both setting rlist and verlet-buffer-tolerance, while generally one excludes the other. GROMACS is telling you that it is ignoring the rlist and using the verlet-buffer-tolerance value to calculate rlist. If you know that you need an rlist of 1.4nm then you can force it by setting verlet-buffer-tolerance=-1. Since I don’t think this is the case, then you can comment off rlist=1.4 as it is not used.

Note 2 tells you that using the v-rescale thermostat is incorrect as you are requiring integrator=sd, which handles internally the temperature. If you need the sd integrator then you can set tcoupl = no or remove altogether the line, as this is the default. Otherwise, you integrate with integrator = md and you keep the v-rescale thermostat.

Then you get three errors. It is hard to go through them as these are due to a incorrectly defined topology. You first start by defining the force field details and the molecules of interest with the respective restraints (as these are part of the individual molecular topology) and then you put the specifics of the [ system ] and the [ molecules ]. Try with something like this

; general force field definitions
#include “dry_martini_v2.1.itp”
; definition of dopc
#include “dry_martini_v2.1_lipids.itp"
#ifdef POSRES
#include “posre_dopc.itp”
#endif
; definition of cholesterol
#include “dry_martini_v2.1_cholesterol.itp”
#ifdef POSRES
#include “posre_chol.itp”
#endif

[ system ]
Expect a large membrane

[ molecules ]
DOPC 1389
CHOL 926
DOPC 533
CHOL 355

I think this still wont work as the dry_martini_v2.1_lipids.itp will include several lipid definitions and not only that of DOPC, so grompp will complain again. If that’s the case then you will have to modify the file and add these lines

#ifdef POSRES
#include “posre_dopcl.itp”
#endif

to the very end of the DOPC molecule definition. It is very likely that these official lipid parameters will already have a restraint definition, though.

8

Dear @obZehn
Thank you so much for your consideration and your time, I will going to do that
Best Regards
Ehsan