[gmx-users] Re: Implementing Radial Distance Constraints

[Mark Abraham]

Bob Johnson wrote:
>>please be more clear.
>>
>>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 
"angres" will help you find all the places you need to add the new 
function and the supporting structure to work out when it gets called.

Mark