[gmx-users] Re: Implementing Radial Distance Constraints
Mark.Abraham at anu.edu.au
Tue Apr 18 03:39:43 CEST 2006
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 ?
> 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
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.
More information about the gromacs.org_gmx-users