Setting up a virtual site to represent a dummy(=massless) LJ sphere

User discussions
1

GROMACS version: 2021.3
GROMACS modification: No

Sorry, for this long post, just trying to summarise and link all the concepts, as I am struggling to wrap my head around things.

I need to create a repulsive area of a given size, it should be able to move around the simulation box as the system relaxes.

I thought I could make use of a dummy/virtual LJ spere with some radius, but w/o mass or charge. Obviously, I need a mass for the MD (=forces), so, honestly, I am confused here.
→ Am I correct to be using dummies for massless repulsive LJ sphere?

My approach is as follows, consulting manual obviously:

  1. add a particle into ffnonbonded.itp:
[ atomtypes ]
;  type atnum   mass        charge  ptype sigma    epsilon
    LJ     0     0.0         0.000     V  0.5             0.001      ; LJ test
 CG331     6    12.011000    0.000     A  0.365268474438  0.3263520  ; aliphatic C for CH3
 CG321     6    12.011000    0.000     A  0.358141284692  0.2343040  ; aliphatic C for CH2
  HGA3     1     1.008000    0.000     A  0.238760856462  0.1004160  ; alphatic proton, CH3
  HGA2     1     1.008000    0.000     A  0.238760856462  0.1464400  ; alphatic proton, CH2

I declare ptype as V for the virtual site, as per Table 10:
→ Am I correct that I need particle to be a virtual site, to be able to have no mass, unlike A for an atom?
→ Do I need to make an entry in atomtypes.atp? or is it only for ptype A? I am reading Atom types entry in the manual and I guess no?

  1. make an LJ.itp for it:
[ moleculetype ]
; Name            nrexcl
LJ                 0

[ atoms ]
; id	at type	res nr 	residu name	at name	cg nr	charge
1          LJ        1       LJ      LJ     1       0.0

Here, I am relying on the mixing rules of the forcefield, as declared under forcefield.itp, for instance, these are CHARMM:

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
       1           2         yes       1.0     1.0

#include "ffnonbonded.itp"
#include "ffbonded.itp"

If I do not want to rely on whatever I get from the mixing rules, I need to declare the values myself. As per Table 13, I add [ nonbond_params ] into the ffnonbonded.itp, for example:

[ nonbond_params ]
  ; i    j func       sigma        epsilon
    LJ    CG331    1   0.5         0.001      ; LJ test
    LJ    CG321    1   0.5         0.001      ; LJ test
    LJ    HGA3     1   0.5         0.001      ; LJ test
    LJ    HGA2     1   0.5         0.001      ; LJ test

Self note: I do not use [ pairs ] & [ pairtypes ] as those are to modify interactions within the molecule itself, see Table 13.

  1. I build a system, for example, for 5 LJ and 20 alkanes in a box my topology.top looks like this:
#include "../charmm36.ff/forcefield.itp"

#include "./c5h.itp"
#include "./LJ.itp"

[ system ]
alkane + LJ

[ molecules ]
; Compound        #mols
C5H 20
LJ 5

Now, I am struggling to understand what is [ virtual_sites1 ] entry is for. In my mind, it was something to deal with the massless dummy? In the manual Virtual sites it is linked to type 1 interaction, but no such interaction is shown on the Fig 34 in the Virtual interaction sites. It also seems to involve 2 atoms and falls under intramolecular interactions, Table 14

I guess, all I want to understand is:
→ Do I need to use [ virtual_sites1 ] when introducing a massless dummy LJ?
→ Is there anything I am missing in my setup?
→ More general Q - If I have a massless atom floating around the box (such as in this example), how are the dynamics of it computed? Would MD only be feasible when there are other real (=with mass) atoms in the box?

Thank you for your advice and sorry for such long post!
V

@mholmboe Michael, I am just tagging you here, since I saw your earlier post and you may have some insight. My Q is a little different, hence, I thought I’d start a new thread rather than piggyback on yours.

2

Im answering since I was tagged, not necessarily meaning I can help you…

→ Am I correct to be using dummies for massless repulsive LJ sphere?

I’m unsure what exactly what you want to achieve. What would be the benefit of using a massless particle? Seems better to use something with a mass, with or without using position restraints and such, and then use a [ nonbond_params ] section. Alternatively, have you considered using (inverted) flat-bottom distance restraints (Fig 30.B)? while using refcoord_scaling=COM
https://manual.gromacs.org/documentation/2021/reference-manual/functions/restraints.html#flat-bottomed-position-restraints

→ Am I correct that I need particle to be a virtual site, to be able to have no mass, unlike A for an atom?

I would think any virtual site has to have one or more parent/s (constructing atoms). Not sure if virtual sites can be single particles, do not think so. Why not use an actual dummy atom with some mass, and/or position restraints?

→ Do I need to make an entry in atomtypes.atp ? or is it only for ptype A ?

I’m unsure. Try running gmx grompp and then gmx dump with the generated .tpr file, to inspect what mass is used for your dummy particle.

→ Do I need to use [ virtual_sites1 ] when introducing a massless dummy LJ?

a virtual_site1 site will only sit on top of its one constructing/parent atom site, as I understand, which means you need an actual atom anyway.

→ Is there anything I am missing in my setup?

I do not think so, but I would not recommend using a massless site or a virtual site.

→ More general Q - If I have a massless atom floating around the box (such as in this example), how are the dynamics of it computed? Would MD only be feasible when there are other real (=with mass) atoms in the box?

I have not tried, but would assume the actual physics would be unreal since the massless particle might achieve insane velocities? I assume Gromacs would not be happy.

All the best,

Michael Holmboe