Setting spherical restraints

GROMACS version: 5.1.5
GROMACS modification: No
Dear all,

I want to set flat-bottom spherical restrains to TIP3P molecules (oxygen atom) in a radius = 20A from the system’s COM, so all those molecules stay trapped in the sphere. I have tried several values for the force, but the molecules are not restrained and diffuse all over the box.

Am I missing something? In the [ position_restrains ] declaration, the atom value refers to the atom number in the input structure or in the TIP3.itp file.

My TIP3.itp is:

[ moleculetype ]
; name	nrexcl
TIP3	     2

[ atoms ]
; nr	type	resnr	residu	atom	cgnr	charge	mass
     1         OT      1     TIP3    OH2      1    -0.834000    15.9994   ; qtot -0.834
     2         HT      1     TIP3     H1      2     0.417000     1.0080   ; qtot -0.417
     3         HT      1     TIP3     H2      3     0.417000     1.0080   ; qtot  0.000

[ settles ]
; OW	funct	doh	dhh
    1     1  9.572000e-02  1.513900e-01

[ exclusions ]
    1    2    3
    2    1    3
    3    1    2

My topology is:

; Include forcefield parameters
#include "toppar/forcefield.itp"

; Include forcefield parameters for TIP3 water
#include "toppar/TIP3.itp"

#ifdef POSRES_WATER_SPH
[ position_restraints ]
;  	atom	type	g	r	k
	1	2	1	2	1000
#endif

[ system ]
; Name
Spherical restrained water in a box

[ molecules ]
; Compound	#mols
TIP3  	        6841

and my .mdp file is:

define                  = -DPOSRES_WATER_SPH
integrator              = md
dt                      = 0.002
nsteps                  = 500000
nstxtcout               = 5000		; 100 
nstvout                 = 5000		; steps 
nstfout                 = 5000		; writen
nstcalcenergy           = 100
nstenergy               = 1000
nstlog                  = 1000
;
cutoff-scheme           = Verlet
nstlist                 = 20
rlist                   = 1.2
coulombtype             = pme
rcoulomb                = 1.2
vdwtype                 = Cut-off
vdw-modifier            = Force-switch
rvdw_switch             = 1.0
rvdw                    = 1.2
;
tcoupl                  = Nose-Hoover
tc_grps                 = SYSTEM
tau_t                   = 1.0
ref_t                   = 300
;
pcoupl                  = Parrinello-Rahman
pcoupltype              = isotropic
tau_p                   = 5.0
compressibility         = 4.5e-5
ref_p                   = 1.0
;
constraints             = h-bonds
constraint_algorithm    = LINCS
continuation            = no
;
refcoord_scaling        = com

Thanks a lot in advance.

Best.

Are the coordinate passed to grompp -r set up correctly? If you want the origin of the restraint to be the center of the system, then all the restrained atoms have to have their coordinates set to this value. This is the somewhat awkward requirement of the flat-bottom restraints.

Thanks for your reply.

Yes, my grompp command was:

gmx grompp -f md_f-b_restrained.mdp -o md_f-b_restrained.tpr -c 1ns_equilibration.gro -r 1ns_equilibration.gro -p topol.top -maxwarn 1

1ns_equilibration.gro is a 1 ns run to equilibrate the density of the system and its COM is the one for the restraints.

On the other hand, it’s there a way to tweak the Heavy side function? I want the restraints to be applied sharply in the boundary, so the molecules can approach very close to the boundary, but bounce back into the sphere.

Best.

This is not correct. What you’re telling grompp is that you want your water oxygen atoms to be allowed to move 2 nm from the positions they currently occupy in 1ns_equilibration.gro. If you want the reference to be the COM of the system, all water oxygen atoms have to have their (x,y,z) coordinates set to the COM. It requires manual (or scripted, better) manipulation of the coordinates passed to grompp -r.

That sounds more like a reflective wall, which is not what the restraint code does, so no, I do not think something like this is supported. You could massively increased the force constant to cause a very hard reversal of diffusion, but that’s about as close as you can come, I think.

So basically I have to set the COM of the system as the origin of the coordinates, right?

As the origin of the restraint potential, yes. That is accomplished by changing the coordinates of the atoms, which then serve as the reference state for which Vrestraint = 0.

So, to be sure if I got the whole idea.

I have to change the coordinates of the whole system, setting its COM as the origin. From this transformed system, I save the coordinates of the water molecules with oxygens up to radius r, and pass them to grompp with the -r option. This will be the restrained subsystem, where the

[ position_restraints ]
; atom   type   g   r   k
  1      2      1   1   1000

directive will apply. Correct?

No, this is not correct.

What you have laid out is that you want a sphere of water in which all the molecules remain within a 2-nm sphere. The restraint potential is applied to the oxygen atom. Therefore, the coordinates passed to grompp -r (1) contain all the coordinates of all the atoms in the system and (2) must be manually edited so that each O atom has (x,y,z) coordinates set to the COM of the system. There is no GROMACS utility that will prepare the coordinate file passed to grompp -r. You have to create/edit it yourself.

Setting the coordinates to the COM of the system is not a problem, I can prepare a script for that. What I was trying to understand is how to apply the restraint only the waters (oxygens) in a sphere around the COM, because in the directive atom = 1 refers to all oxygens.

When I run:

gmx grompp -f md_f-b_restrained.mdp -o md_f-b_restrained.tpr -c 1ns_equilibration.gro -r 1ns_equilibration.gro -p topol.top -maxwarn 1

with

[ position_restraints ]
; atom   type   g   r   k
  1      2      1   1   1000

that restraints each oxygen atom of the system inside its “own” sphere of radius = 10A, got it.

Then, maybe flat-bottom restrain is not the optimal approach for this. So allow me to rephrase my question, how can I encapsulate water molecules inside a sphere around a point in the system (ej. COM)? (Although you said that a hard wall potential is not supported.)

If you only want the restraint to apply to a subset of the system, you have to generate a distinct topology for those water molecules, even though the contents will be redundant.

Thank you very much.

Hello,

Following up with this issue, I set the coordinates of the water oxygens of interest to the COM of the system. However, when I run mdrun get:

Fatal error:
Step 0: The total potential energy is 1.28914e+60, which is extremely high.
The LJ and electrostatic contributions to the energy are 1.28914e+60 and
1.57237e+10, respectively. A very high potential energy can be caused by
overlapping interactions in bonded interactions or very large coordinate
values. Usually this is caused by a badly- or non-equilibrated initial
configuration, incorrect interactions or parameters in the topology.

This is kind of expectable since all those molecules share the same oxygen coordinates. Is there any way to skip/solve that error, without doing a minimization/equilibration?

The manipulated coordinates should be passed to grompp -r, not grompp -c.

Thanks a lot! it worked.

One last question. How would be declaration for the inverted flat-bottom potential (http://manual.gromacs.org/documentation/2020/reference-manual/functions/restraints.html#flat-bottomed-position-restraints, Fig. 30B)?

As noted below the figure, " In addition, it is possible to invert the restrained region with the unrestrained region, leading to a potential that acts to keep the particle outside of the volume defined by Ri, g, and rfb. That feature is switched on by defining a negative rfb in the topology."

Thanks!