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

[Justin Lemkul]

On 5/5/20 5:54 PM, Debashish Banerjee wrote:
> 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.

http://manual.gromacs.org/current/user-guide/run-time-errors.html#atom-index-n-in-position-restraints-out-of-bounds

Presumably your Cs+ topology defines the [moleculetype] as having one 
ion in it. Therefore, the only valid atom number in 
[position_restraints] is 1. The global atom numbering is meaningless in 
the context of position restraints (which are assigned at the topology 
level) as genrestr warns you (see the help description; if you're trying 
to restrain anything that isn't the first molecule, you can't use it).

The solution is to define a new [moleculetype] for the Cs+ ion you want 
to restrain, and use

[position_restraints]
1 1  0 0 1000

in its topology.

-Justin

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

Justin A. Lemkul, Ph.D.
Assistant Professor
Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

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