[gmx-developers] Change restraints in run-time?

[Gerrit Groenhof]

> 
> Hi all,
> 
> Every n md steps i like to make changes to (distance,angular) restraints, that 
> is: add, remove or change restraints. Since the data of these restraints is 
> located in the t_topology structure, I wondered if it is possible to do this 
> by just adding information in this structure. I expect problems here since 
> this structure is kept staticly(?) during simulation. 
> 
> I suspect that I have to hack a workaround to get these restraints somewhat 
> dynamic. Honestly, i would not like to touch the internal force calculations, 
> but simply stick to a more higher level. Am I overlooking essential 
> information here (is pull code relevant?), or is something just not possible?
> 
> Gerwin de Haan
> gerwin at mindbus.nl
>


Hi Gerwin,

If I understand you correctly, you want to modify restraints during
a simulaiton. In the case of restraints I think it is not too
difficult.

In the idef struct there is a so-called ilist struct for each
interaction type in gromacs (F_NRE). Restraints is one of those (see
the enum at the top of the idef.h file). 

In the idef (idef.h) there are two important entries you have to
consider for what you want: the array of iparams structs and the array
of ilist structs. there is also an functype array, but that you don't
have to modify i think.

Now the idef "works" as follows. I might be instructive to keep an eye
on http://md.chem.rug.nl/~groenhof/idef.html while reading this.

ilist is a struct which has an array called iatoms[] taht contains
grouped blocks of atoms that have a certain interaction of type
F_whatever
(see the enum at the top of idef.h) with each other (bond, angle, etc)
AND also a number (or index) that points to a location in the iparams
and functype arrays of the idef. These arrays contain resp. the parameters
(t_iparams in idef.h) and the functiontpe (just an int). The latter is
needed by gromacs to select the proper funciton while computing the
interaction and in that computation it uses the parameters of the
iparams array. So with the content of the iatoms array of a certain
interaction type (idef.il[F_whatever].iatoms) gromacs can compute the
actual interaction. See the idef sheet. For instance, the function
that does the disres computation is ta_disres (in disres.c) and from
its arguments you can see that indeed the iparam (ip[])and iatom
(forceatoms[]) structs are passed to this function. Moreover in the
loop in this function you can clearly see how first the type is
exctracted (type = forceatoms[fa]) and then how this "type" is used to
acces the parameters (up1 = ip[type].disres.up1 and so on)

If you look at the iparams struct in the idef.h you will see it
contains a struct called disres. This struct contains the restraint
parameters. So I believe that you shoul somewhere in the md loop
(mdrunner.c if it hasn't changes in the mean time....) add some
statements to change the content of the elements of the
top->idef.iparams[n].disres struct (low,up1,up2,kfac).  The index n
stands for which of the disres interactions you want to change. you
can obtain n from the top->idef.il[F_DISRES].iatoms[] array (like in
the ta_disres). If you want to change the disres parameters of only
certain groups of atoms, you can go through the iatoms[] array and scan
(with step of three (one for the index that pointS to the iparams
array and two for the two atoms involved)) for the ones you need. For
the latter I suggest you read my email to Hank the Haan (a far
relative maybe?) about defining and accessing group in gromacs.

I hope this will work. And if not(actually I am not 100% sure) i hope that  at least I helped you somewhat in the right direction.

Good luck,

Gerrit