[gmx-users] How to add a repulsive harmonic potential

[atsutoshi.okabe at takeda.com]

>Depending on how structurally complex the ligands are, it may be very easy to 
>define a dummy atom within the structure using a virtual site, then specifying a 
>nonbonded interaction only between these virtual sites using [nonbond_params]. 
>That way, these sites will not affect ligand-water, ligand-protein, etc but will 
>repel ligands from one another.

Thanks for your reply.
I tried to use [virtual_sitesn] and [nonbond_params] below.
VS1 means virtual site( dummy atoms) of center of mass of one ligand_1.
VS2 means virtual site( dummy atoms) of center of mass of one ligand_2.
Complex.top file

[atomtypes]
;name  bond_type mass charge  ptype  sigma        epsilon
.............
VS1     VS1    0.0000 0.0000  V   0.00000e+00   0.00000e+00
VS2     VS2    0.0000 0.0000  V   0.00000e+00   0.00000e+00

[ nonbond_params ]
;i   j   func      V(c6)    W(c12)
VS1  VS2   1   0.00000e+00  5.8954e+03

[atoms]
;nr    type   resnr residue atom cgnr  charge  mass
.............
100      ha     2   LIG    H7    100    0.21132    1.00000
101      VS1    3   VIR    VS1   101    0.00000    0.00000   
102      VS2    4   VIR    VS2   102    0.00000    0.00000

[ virtual_sitesn ]
; Site funct from
101    2      25 26 27 28
102    2      75 76 77 78

However, I got fatal error below.
Fatal error:
Number of coordinates in coordinate file(NVT_1.gro, 77111)
            Does not match topology (complex.top, 77113)

Does it mean that I have to add coordinate information of dummy atoms(VS1, VS2) in .gro file ?
Or is how to define virtual site incorrect?

Best regards,
Atsutoshi Okabe


> Thank you for your help.
> Actually, I want to predict the binding mode of ligand for protein/ligand complex during long time simulation, so I want to add only repulsive potential between ligands to avoid ligand aggregation, not fix these ligands.
> The reason why I add some ligand molecules(not one ligand) in simulation is to enhance ligand sampling.
>
> I tried using [distance_restaints] below, not use [bonds] because of domain decomposition error. Then, the MD simulation could run normally without domain decomposition error!
> [ distance_restraints ]
> ;ai aj type index type' low up1 up2 fac
> 13  63  1    1     1    20   100 100 10000
>
> However, the distance between 13 atom and 63 atom seems to be less than 20 angstrom in every steps of MD simulation.
> So I tried to change the fac value from 10000 to -10000, but the distance was less than 20 angstrom during MD simulation
> Why the ligands were restrained less than 20 angstrom although I think I add repulsive force between these atoms?
>

Depending on how structurally complex the ligands are, it may be very easy to 
define a dummy atom within the structure using a virtual site, then specifying a 
nonbonded interaction only between these virtual sites using [nonbond_params]. 
That way, these sites will not affect ligand-water, ligand-protein, etc but will 
repel ligands from one another.

-Justin