[gmx-users] Constraining two virtual sites

[Mark Abraham]

On 6/12/2010 9:43 PM, Sebastian Fritsch wrote:
> Hi everybody,
>
>
> I have some trouble setting up topology for my, I admit, quite unusual 
> system. I need to constrain the center of masses of two groups of 
> atoms to a fixed 'bond' length. Since I have to do this for every 
> molecule I cannot use the 'pull constraint' of the free energy code.
>
> Thus in my topology I defined a com vsite for each group and tried to 
> apply a constraint between them. The atoms within each group are 
> connected by harmonic springs forming a ring polymer (for those of you 
> who know it: this is supposed to become the path integral 
> representation in the end...)
> The topol.top for a single molecule test system looks like this:
>
> [ defaults ]
> ; nbfunc   comb-rule    gen-pairs    fudgeLJ    fudgeQQ
>  1       1        no        1.0    1.0
> [ atomtypes ]
> ;name  mass        charge    ptype  C6             C12           ; 
> sigma     epsilon
> A      1.00800     0  A      0           0
> B      1.00800     0  A      0           0
> C      1.00800     0  A      0           0
> VS     0.000       0   V      0           0
> [ moleculetype ]
> ; molname      nrexcl
> SOL        0
> [ atoms ]
> ; id at resnr resnm atnm       cgnr      charge
> 1  VS    1  A         V1        1         0
> 2  VS    1  B         V2        2         0
> 3  A     1  A         A          3         0 4  B     1  A         
> B          4         0 5  C     1  A         C          5         0 6  
> A     1  B         A          6         0 7  B     1  B         
> B          7         0 8  C     1  B         C          8         0 [ 
> virtual_sitesN ]
> ;
> 1 2 3 4 5
> 2 2 6 7 8
> [ bonds ]
> ; ai aj
> 3 4 1 0 1687
> 4 5 1 0 1687
> 5 3 1 0 1687
> 6 7 1 0 1687
> 7 8 1 0 1687
> 8 6 1 0 1687
> [ constraints ]
> ; ai   aj funct   length_A
> 1 2 1  0.1633
> [ system ]
> TEST
> [ molecules ]
> SOL     1
>
> The conf.gro is setup such that the constraint is satisfied initially.
>
> When I run the simulation with LINCS the constraint seem to be not 
> applied at all, the distance beetween the two vsites increases just as 
> in free diffusion. With SHAKE the system explodes immediately (without 
> error message, but every value is inf or nan).
> I am running the system with the sd integrator in (vers. 4.5.3) but 
> also tried md and no thermostat, with no change in results.
>
> I appreciate any advice on this!
>
> Thanks,
> Sebastian

I can't see anything that looks wrong. Make sure you inspect the output 
of grompp thoroughly for clues about how it's interpreting the .top 
file. You may like to use gmxcheck to compare the .tpr files produced 
under various permutations to see that things look like they're working 
correctly (remember atom numbers will start from zero in the output). 
Try putting the virtual site atoms after the real atoms in the ordering.

Mark