# Virtual sites for repulsive potential

**URL:** <https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092>\
**Category:** User discussions\
**Created:** [November 6, 2021, 1:41pm UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092 "2021-11-06T13:41:00Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![rainbow](https://avatars.discourse-cdn.com/v4/letter/r/7993a0/32.png) [@rainbow](https://gromacs.bioexcel.eu/u/rainbow)\
**Post date:** [November 6, 2021, 1:41pm UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/1 "2021-11-06T13:41:00Z")

</div>

GROMACS version: 2021.2  
GROMACS modification: No

Dear gmx users,  
I’m doing simulations using hundreds of benzene and propane.  
To prevent aggregation of these molecules, I trying to apply virtual sites.

===== Benzene.gro =====  
13  
1ben C1 1 -0.001 0.001 -0.006  
1ben C2 2 0.003 -0.067 -0.128  
1ben C3 3 0.007 -0.206 -0.130  
1ben C4 4 0.006 -0.277 -0.011  
1ben C5 5 0.002 -0.210 0.110  
1ben C6 6 -0.001 -0.071 0.113  
1ben H1 7 -0.003 0.109 -0.004  
1ben H2 8 0.004 -0.011 -0.221  
1ben H3 9 0.010 -0.258 -0.225  
1ben H4 10 0.009 -0.385 -0.013  
1ben H5 11 0.002 -0.266 0.203  
1ben H6 12 -0.004 -0.019 0.208  
1ben M1 13 0.003 -0.138 -0.009  
0.27600 9.88200 8.65000

===== Benzene.itp =====  
[atomtypes]  
;name bond\_type mass charge ptype sigma epsilon  
Mm Mm 0.00000 0.00000 D 0.00000e+00 0.00000e+00

* * *

[virtual\_sitesn]  
; Site funct from  
13 2 1 4

or

[virtual\_sites2]  
; Site from funct a  
13 1 4 1 0.5

I wrote as the above to create a virtual site for center of mass.  
According to the manual, there were two ways. Is both ways are possible? if not, how can I modify it?

* * *

[nonbond\_params]  
; i j func V(c6) W(c12)  
Mm Mm 1 1.2000e+00 -0.4184e-01

I used V, W values from the paper, [https://doi.org/10.1021/ci100462t](https://doi.org/10.1021/ci100462t)  
(ε = -0.01 kcal/mol; Rmin = 12.0 Å).  
I am using amber99sb-ildn forcefield. It uses combination rule 2.  
Is it right?

===== Propane.gro =====  
12  
1pan C1 1 0.000 0.000 0.000  
1pan C2 2 0.150 0.000 0.000 \<—  
1pan C3 3 0.200 -0.141 0.000  
1pan H1 4 -0.036 0.103 0.000  
1pan H2 5 -0.036 -0.051 -0.089  
1pan H3 6 -0.036 -0.051 0.089  
1pan H4 7 0.186 0.051 -0.089  
1pan H5 8 0.186 0.051 0.089  
1pan H6 9 0.309 -0.141 0.000  
1pan H7 10 0.164 -0.193 0.089  
1pan H8 11 0.164 -0.193 -0.089  
1pan M2 12 0.150 0.000 0.000 \<—  
6.90600 5.91200 3.56000

===== Propane.itp =====  
[atoms]  
; nr type resi res atom cgnr charge mass ; qtot bond\_type  
12 Mm 1 pan M2 1 0.000000 0.00000

* * *

[virtual\_sitesn]  
; Site funct from  
12 2 2

or

;[virtual\_sites1]  
; Site from funct  
12 2 1

This dummy atoms is for middle carbon of propane.  
Is this right?

* * *

using the above parameters,  
The following errors obtained; after energy minimization and production run, respectively.

_Energy minimization has stopped because the force on at least one atom is not_  
_finite. This usually means atoms are overlapping. Modify the input_  
_coordinates to remove atom overlap or use soft-core potentials with the free_  
_energy code to avoid infinite forces._

_Step 0, time 0 (ps) LINCS WARNING_  
_relative constraint deviation after LINCS:_  
_rms 0.125466, max 19.224945 (between atoms 9928 and 9934)_  
_bonds that rotated more than 30 degrees:_

- atom 1 atom 2 angle previous, current, constraint length\*
- 5733 5739 90.0 0.1086 0.2882 0.1086\*

How can i solve this problem?  
Will it works if i applied constraints for virtual sites?  
There is two atoms for this case in the manual, but in my case only one dummy atom.

[constraints]  
; ai aj funct length\_A length\_B

No matter how much I search, I could’t find a solution…  
Please guide me,  
Kind regards,

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [November 6, 2021, 4:52pm UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/2 "2021-11-06T16:52:54Z")

</div>

> [@rainbow](#):
>
> [nonbond\_params]  
> ; i j func V(c6) W(c12)  
> Mm Mm 1 1.2000e+00 -0.4184e-01
> 
> I used V, W values from the paper, [https://doi.org/10.1021/ci100462t](https://doi.org/10.1021/ci100462t)  
> (ε = -0.01 kcal/mol; Rmin = 12.0 Å).  
> I am using amber99sb-ildn forcefield. It uses combination rule 2.  
> Is it right?

No - you need to convert Rmin to σ.

---

<div class="post-metadata">

**Author:** ![rainbow](https://avatars.discourse-cdn.com/v4/letter/r/7993a0/32.png) [@rainbow](https://gromacs.bioexcel.eu/u/rainbow)\
**Post date:** [November 7, 2021, 11:23am UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/3 "2021-11-07T11:23:45Z")

</div>

Where can I find this information? I’ve read in the manual but I didn’t find a solution.  
Could you give me a little more detail?

And, are the rest things correct?  
Then are the above errors caused by not converting the Rmin value? or should I consider constraints or something else?

I’d appreciate it if you could tell me in more detail.  
Thanks…

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [November 8, 2021, 2:37pm UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/4 "2021-11-08T14:37:47Z")

</div>

Chapter 5 of the manual is your best friend if you’re doing anything complicated related to topologies. Convert to a proper σ value and see if you get a sensible energy. If not, you need to simplify to a toy system whose energy you could compute by hand if needed, to check the stability of the topology.

---

<div class="post-metadata">

**Author:** ![rainbow](https://avatars.discourse-cdn.com/v4/letter/r/7993a0/32.png) [@rainbow](https://gromacs.bioexcel.eu/u/rainbow)\
**Post date:** [November 12, 2021, 2:32am UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/5 "2021-11-12T02:32:57Z")

</div>

Thanks but, no matter how much I read the manual, it is still difficult for me.  
I tried convert the value but still the aggregation keeps happening.  
It’s been stuck for several weeks without a solution.  
I would really appreciate it if you could teach me…  
Thanks…

[Benzene]

13  
1ben C1 1 -0.001 0.001 -0.006  
1ben C2 2 0.003 -0.067 -0.128  
1ben C3 3 0.007 -0.206 -0.130  
1ben C4 4 0.006 -0.277 -0.011  
1ben C5 5 0.002 -0.210 0.110  
1ben C6 6 -0.001 -0.071 0.113  
1ben H1 7 -0.003 0.109 -0.004  
1ben H2 8 0.004 -0.011 -0.221  
1ben H3 9 0.010 -0.258 -0.225  
1ben H4 10 0.009 -0.385 -0.013  
1ben H5 11 0.002 -0.266 0.203  
1ben H6 12 -0.004 -0.019 0.208  
1ben M1 13 0.003 -0.138 -0.009  
0.27600 9.88200 8.65000

[atomtypes]  
;name bond\_type mass charge ptype sigma epsilon Amb  
ca ca 0.00000 0.00000 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860  
c3 c3 0.00000 0.00000 A 3.39967e-01 4.57730e-01 ; 1.91 0.1094  
ha ha 0.00000 0.00000 A 2.59964e-01 6.27600e-02 ; 1.46 0.0150  
hc hc 0.00000 0.00000 A 2.64953e-01 6.56888e-02 ; 1.49 0.0157  
Vb Vb 0.00000 0.00000 V 0.00000e+00 0.00000e+00

[nonbond\_params]  
; i j func V(c6) W(c12)  
Vb Vb 1 1.0691e-01 4.1840e-02

[moleculetype]  
;name nrexcl  
ben 3

[virtual\_sites2]  
; Site from funct a  
13 1 4 1 0.5

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [November 12, 2021, 2:48am UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/6 "2021-11-12T02:48:45Z")

</div>

You simply calculated σ wrong. It’s 10x too small; it should be 1.0691, not 0.10691.

---

<div class="post-metadata">

**Author:** ![rainbow](https://avatars.discourse-cdn.com/v4/letter/r/7993a0/32.png) [@rainbow](https://gromacs.bioexcel.eu/u/rainbow)\
**Post date:** [November 12, 2021, 8:16am UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/7 "2021-11-12T08:16:31Z")

</div>

ε = -0.01 kcal/mol; Rmin = 12.0 Å

For benzene system,  
I already used 1.0691 (Rmin = 1.2 nm). The value is calculated from the  
following formula.  
 ![image](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/b/ba0909ca26ab492976018345801a89e1bc6d8268.png)

But is gives LINCS warning;

Step 8, time 0.016 (ps) LINCS WARNING  
relative constraint deviation after LINCS:  
rms 21543786.000000, max 3833009152.000000 (between atoms 8662 and 8668)  
bonds that rotated more than 30 degrees:  
atom 1 atom 2 angle previous, current, constraint length  
5304 5310 90.0 0.1489 1.4612 0.1086  
5308 5314 90.0 0.1086 0.5957 0.1086  
5603 5609 90.0 0.1086 1.0625 0.1086  
5606 5612 90.0 0.1086 0.9689 0.1086  
6308 6314 90.0 2.2018 2.2630 0.1086  
7049 7055 35.7 0.1086 0.1086 0.1086  
7683 7689 90.0 0.1386 0.8829 0.1086

Then I set it to 0.8691, (by reducing 0.2nm), it successfully running (sigma 0.8691e+00, epsilon 4.1840e-02)  
\*\* If it is set to -4.1840e-02, a LINCS warning occurs, so I changed to a positive number.

However, when benzene and propane are present together, a LINCS warning occurs if the same dummy atom type is set for benzene (center of mass) and propane (center carbon) using sigma for 0.8691e+00, epsilon for 4.1840e-02.

Even if the dummy atom type is set differently for each molecule, LINCS warning also occurs (Vb for benzene, Vp for propane).  
Vb Vb 1 0.8691e+00 4.1840e-02  
Vb Vp 1 0.8691e+00 4.1840e-02  
Vp Vp 1 0.8691e+00 4.1840e-02

I don’t know what the problem is. I would be grateful if you could help me get through this step.  
Thanks…

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [November 12, 2021, 1:37pm UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/8 "2021-11-12T13:37:09Z")

</div>

The value of σ I quoted comes directly from the SilcsBio distribution, which is from the paper you referred to earlier. ε should be positive by GROMACS convention (it is the magnitude of the well depth, not the actual energy itself).

If that’s still crashing, you have a different problem and you need to simplify the system and verify that molecules without virtual sites work, that you can successfully energy-minimize a single molecule (to ensure that the virtual site is being constructed properly) and a dimer system in which the virtual sites interact such that you have a predictable energy.

---

<div class="post-metadata">

**Author:** ![Audreyl](https://avatars.discourse-cdn.com/v4/letter/a/e19b73/32.png) [@Audreyl](https://gromacs.bioexcel.eu/u/Audreyl)\
**Post date:** [December 17, 2021, 1:07pm UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/9 "2021-12-17T13:07:48Z")

</div>

Hi rainbow,  
Did you use a script to include the virtual site in your Benzene.gro? I am not a programming person and will be glad if you could help me with this. I am curently struggling with these virtual sites for months now.

Thanks

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [December 17, 2021, 4:42pm UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/10 "2021-12-17T16:42:55Z")

</div>

A virtual site placed in the center of a ring is as simple as assigning the average carbon coordinates to the virtual site. Or you can actually just give it some relatively random coordinates (sensible but not an exact value) and perform energy minimization, the virtual site will be reconstructed at the center of the ring.

---

<div class="post-metadata">

**Author:** ![Audreyl](https://avatars.discourse-cdn.com/v4/letter/a/e19b73/32.png) [@Audreyl](https://gromacs.bioexcel.eu/u/Audreyl)\
**Post date:** [December 18, 2021, 2:36pm UTC](https://gromacs.bioexcel.eu/t/virtual-sites-for-repulsive-potential/3092/11 "2021-12-18T14:36:29Z")

</div>

Hi jalemkul,

Thank you so much for the explanation it really helped me.