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

Debashish Banerjee deb.aerospace at gmail.com
Thu May 7 17:19:21 CEST 2020


Hello gmx users,

I am stuck at a very basic command line and need your expertise in dealing
through it.

Basically, I have a clay sheet containing 192 surface oxygen. Since I am
doing umbrella sampling for PMF calculations, I want to get the progression
of COM distances of my fixed Cs atom from the center of mass of the clay
sheet (basically these 192 surface oxygen) over time.

I have already created special index file which has relevant groups such as
;
group 18 --> *referenceSurface* ~containing 192 surface oxygen
group 19 --> *restrainedCs* ~containing 1 Cs atom (basically the one that I
have fixed)
group 20 --> another Cs atom (this one is not fixed and can move freely in
all sample space during the entire course of simulation)

Command line:

*Case 1*

*gmx_mpi distance -n index.ndx -f conf0.gro -s pull.tpr -select 'com of
group 19 plus com of group 18' -oxyz dist00.xvg*

Inconsistency in user input:
Invalid index group reference(s)
  Group 'referenceSurface' cannot be used in selections except as a full
value
  of the selection, because atom indices in it are not sorted and/or it
  contains duplicate atoms.

*Case 2*
*gmx_mpi distance -n index.ndx -f conf0.gro -s pull.tpr -select 'com of
group 19 plus com of group 20' -oxyz dist2.xvg*

Reading file npt.tpr, VERSION 2016.4 (single precision)
Reading frames from gro file 'Generated by trjconv : clay + Cs + acetate +
water t=   2.00000', 24221 atoms.
Last frame          0 time    2.000
Analyzed 1 frames, last time 2.000
com of group 19 plus com of group 20:
 Number of samples:  1
 Average distance:   0.52370  nm
 Standard deviation: 0.00000  nm

My understanding is that if I compare the above 2 cases, I would say that
for *case 2*, Cs is measuring it's distance from the only other Cs atom. No
problem here, whereas for *case 1*, Cs atom is measuring it's distance from
all 192 oxygen atoms, somehow not understanding how to take the center of
mass of the whole sheet containing 192 atoms.

Approach taken:
1. Modified my command line from this:
 gmx_mpi distance -s pull.tpr -f conf0.gro -n index.ndx -select 'com of
group "restrainedCs" plus com of group "referenceSurface"' -oall dist0.xvg
2. Used all output options like -oav / -oxyz

I have followed the following threads, but still couldn't manage to rectify
this error
https://www.mail-archive.com/gromacs.org_gmx-users@maillist.sys.kth.se/msg16576.html

https://www.mail-archive.com/gromacs.org_gmx-users@maillist.sys.kth.se/msg30552.html


Please help. I am pretty much sure that I am overseeing a very basic step.

Regards,
Debashish


On Wed, May 6, 2020 at 12:40 AM Debashish Banerjee <deb.aerospace at gmail.com>
wrote:

> Thank you Justin, you are a life saver..!!!!!
> It's working now. You have correctly pointed out the issue.
> Thank you so much...... :)
>
> Regards,
> Debashish
>
> On Wed, May 6, 2020 at 12:11 AM Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> 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
>>
>> ==================================================
>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-request at gromacs.org.
>
>
>
> --
> *Best Regards*
>
> *Debashish Banerjee*
>
> *Ph.D. (Nuclear Materials, Subatech)*
> *MS(Sustainable Nuclear Engineering)*
> *Advance Nuclear Waste Management*
> *Institut Mines-Telecom, **France*
>


-- 
*Best Regards*

*Debashish Banerjee*

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


More information about the gromacs.org_gmx-users mailing list