[gmx-users] How to apply position restrain on selected atoms in .gro files in Gromacs?

[Debashish Banerjee]

Dear gmx users,

I have been stuck in a problem for a month now and have scrutinized all the
past discussions and forms related to it dating as back as 9 years. Now I
really need your help.

My system has Clays (layered philosillacates), Cesium ions (Cs), acetate
molecules as well as water.  In total 24241 atoms in the whole system. I
want to perform free energy calculations *(PMF)* via *umbrella sampling*.

The whole idea is that I want to fix one Cs atom (atom number/index number~
5790) near to the basal surface (top layer) of clay and put in use of the*
pull* code as described in umbrella sampling and thereafter perform *WHAM*
calculations  to see the progression of COM distance between (Cs) and the
surface oxygen's present on the top layer of clay sheet over time.

When, I see it in my .gro files, I can find all 16 Cs atoms ranging from
index number (5785-5780). I just want to fix a particular single Cs atom
(index number 5790) as mentioned above. For this I create a special index
file which now has an extra group containing just 1 Cs atom. After this, I
take in the use of gmx genrestr command :
*gmx genrestr -f npt.gro -n index.ndx -o posre_Cs.itp -fc  0  0  1000*
This commands redirects to the newly created group from the index file
which I created above and put constraints in it.

 *Step 1   *
The *posre_Cs.itp* file looks like this:
[ position_restraints ]
;  i            funct       fcx        fcy        fcz
   5790     1             0            0         1000

* Step 2 *
*My topology should include an if statement, so I also define it like this:*

#include "../charmm36.ff/forcefield.itp"

; include params for ClayFF
#include "../ClayFF.ff/ffnonbonded.itp"

; Include topology for Clay Montmorillonite (MMT)
#include "../ClayFF.ff/MMT_UC/UC2.itp"
#include "../ClayFF.ff/MMT_UC/UC1.itp"
#include "../ClayFF.ff/MMT_UC/UC3B.itp"
#include "../ClayFF.ff/MMT_UC/UC3T.itp"

; Added the part of defining ions.itp which contains Cs atom and defined
position restrain parameters




*#include "../ClayFF.ff/ions.itp"#ifdef POSRES_Cs#include
"posres_Cs.itp"#endif*
#include "../charmm36.ff/spc.itp"

#include "../charmm36.ff/organics/act.itp"

[ system ]
; Name
clay + Cs + acetate + water

[ molecules ]
; Compound        #mols
UC2        1
UC1        1
UC3B      1
UC2        1
UC2        1
UC1        1
UC1        1
UC2        1
UC2        1
UC3T      1
UC1        1
UC2        1
UC2        1
UC3B      1
UC3T      1
UC2        1
UC2        1
UC1        1
UC3B      1
UC2        1
UC2        1
UC1        1
UC1        1
UC2        1
UC2        1
UC3T      1
UC1        1
UC2        1
UC2        1
UC3B     1
UC3T     1
UC2       1
Na          24
SOL       312
UC2       1
UC1       1
UC3B     1
UC2       1
UC2       1
UC1       1
UC1       1
UC2       1
UC2       1
UC3T     1
UC1       1
UC2       1
UC2       1
UC3B     1
UC3T     1
UC2       1
UC2       1
UC1       1
UC3B     1
UC2       1
UC2       1
UC1       1
UC1       1
UC2       1
UC2       1
UC3T     1
UC1       1
UC2       1
UC2       1
UC3B     1
UC3T     1
UC2       1
Na          24
SOL        312
UC2       1
UC1       1
UC3B     1
UC2       1
UC2       1
UC1       1
UC1       1
UC2       1
UC2       1
UC3T     1
UC1       1
UC2       1
UC2       1
UC3B     1
UC3T     1
UC2       1
UC2       1
UC1       1
UC3B     1
UC2       1
UC2       1
UC1       1
UC1       1
UC2       1
UC2       1
UC3T     1
UC1       1
UC2       1
UC2       1
UC3B     1
UC3T     1
UC2       1
Na          24
Cs          16
act          16
SOL        6103

*Step 3*
In my mdp file, I defined it as follows:
*define = -DPOSRES_Cs*


*The problem : *
gmx_mpi   grompp  -f npt.mdp  -c  em.gro   -p topol.top   -o npt.tpr

Setting the LD random seed to 1967282441
Generated 83825 of the 83845 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 1
Generated 57459 of the 83845 1-4 parameter combinations

ERROR 1 [file posres_Cs.itp, line 5]:
Atom index (5790) in position_restraints out of bounds (1-1).
This probably means that you have inserted topology section
 "position_restraints"  in a part belonging to a different molecule than
you intended to.
  In that case move the "position_restraints" section to the right molecule.

There were 2 notes
-------------------------------------------------------
Program:     gmx grompp, version 2016.4
Source file: src/gromacs/gmxpreprocess/toppush.cpp (line 1880)

Fatal error:
There was 1 error in input file(s)


Kindly help me fix it please.

With lots of hope of finding an answer....
-- 
*Best Regards*

*Debashish Banerjee*

*Ph.D.  3rd year (Nuclear Materials, Subatech)*
*MS(Sustainable Nuclear Engineering)*
*Advance Nuclear Waste Management*
*Institut Mines-Telecom, **France*