Virtual sites for repulsive potential

User discussions
1

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
(ε = -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,

2

No - you need to convert Rmin to σ.

3

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…

4

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.

5

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

6

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

7

ε = -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

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…

8

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.

9

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

10

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.

11

Hi jalemkul,

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