# [gmx-users] Re: Implementing Radial Distance Constraints

Mark Abraham Mark.Abraham at anu.edu.au
Tue Apr 18 03:39:43 CEST 2006

```Bob Johnson wrote:
>>
>>if we assume your cnt is centered along the z-axis, given roughly by
>>x^2+y^2 = d0. Am I correct in assuming that you want all other particles
>>to have x^2+y^2 <= d1, where d1 < d0 ?
>>
>>--
>>David.
>
>
> The molecules in my system adsorb to the tube's outer wall. They are not inside
> the tube. Typically, the adsorbates reside 3.5A above the nanotubes surface (a
> typical distance for molecules interacting via VdW with an aromatic molecule).
> I want to constrain molecules to always adopt a monolayer about the nanotube.
> Therefore, I want to put a constraint on certain atoms to reside 3.5A above the
> nanotube. Let's say my nanotube has a radius of 0.5A. I would then want to
> implement a potential given by V = 1/2 k (R-0.5A)^2 where R = sqrt(x^2+y^2).
> This way, the atoms are able to slide along the z and cirumferential
> directions.

That potential will attempt to drag the adsorbates down to the radius of
the nanotube - vdW forces will oppose this and equilibrium will occur
somewhere, but probably inside 3.5A. I think there are less perturbative
ways of applying these constraints. A potential like half of Figure
4.1.3 in the gromacs manual would make more sense, only switching on
when R > 3.5+0.5+x for some small x that acknowledges that 3.5 is an
average distance. That is

V = 0 for R <= R0
V = 0.5 k ( sqrt(x^2+y^2) - R0 ) ^2 for R > R0

Since you expect R-R0 to be small, you might consider implementing the
approximate potential V = 0.5 k ( R^2 - R0^2 ) for computational efficiency.

As far as finding out what to change goes, you'll need to change grompp
to recognize a new type of cylindrical constraint, and then provide a
function that evaluates the potential. There are a bunch of directives
that affect the function of distance and orientation restraints for MD
which won't want to concern you, so I suggest using the code for angle
restraints as a template to adapt to your purposes. In the gromacs src
directory, "grep angle_restraints */*.[ch]" and also with "*/*/*.[ch]"
to find the headers and code sections where grompp interprets the .top
file. That will give a clue for the names of relevant data structures
and what needs to be added where. You will need to modify
src/gmxlib/bondfree.c to add a new potential type - grepping for