[gmx-users] Constraining two virtual sites

[Sebastian Fritsch]

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