[gmx-users] Tabulated potential between atoms and virtual sites

[Thomas Tarenzi]

Dear Gromacs users,

I would like to use a tabulated potential between the C-alphas of a protein and the virtual sites corresponding to the center of mass of the water molecules. I am also using a tabulated potential to define the interactions between these virtual sites of the solvent. In the mdp file I added the following lines:

energygrps = B CA atoms (B is the water virtual site and atoms includes all the rest of the system)
energygrp_table = B B B CA

However, I see in the log file that only the file table_B_B.xvg is read, but not the file table_B_CA.xvg (and in fact, I see from the simulation that I don’t have interactions). On the opposite, if I apply the same table to the interactions between CA and atoms of the water molecules (not virtual sites), it works. Does it mean that I cannot apply tabulated potentials between atoms and virtual sites, or that I should change something in the definition of the virtual site?

They are defined in an itp file in this way:

[ moleculetype ]
; molname      nrexcl
SOL             2

[ atoms ]
; id    at type res nr  residu name     at name
1 ow    1  SOL       ow        1        -0.8476
2 hw    1  SOL       hw1       1         0.4238
3 hw    1  SOL       hw2       1         0.4238
4 B     1  SOL       B         2         0.0000

[ virtual_sites3 ]
; Site from funct a d
4 1 2 3 1 0.05595E+00 0.05595E+00

#ifdef FLEXIBLE
[ bonds ]
; i     j       funct   length  force.c.
  1     2       1       0.1     345000
  1     3       1       0.1     345000
  
[ angles ]
; i      j      k       funct   angle   force.c.
  2      1      3       1       109.47  383

#else
[ settles ]
; i     j       funct   length
1       1       0.1     0.1633
#endif

[ exclusions ]
1       2       3       4
2       1       3       4
3       1       2       4

In the topology I added these lines:

[ atomtypes ]
;name  mass        charge    ptype  C6             C12           
 ow    15.99940   -0.8476    A      2.617065E-03   2.633456E-06  
 hw    1.00800     0.4238    A      0.000000E+00   0.000000E+00
 B    0.00000     0.0000    V      0.000000E+00   0.000000E+00

[ nonbond_params ]
; i   j   func c6             c12
 B B  1    0.000000E+00   1.000000E+00