Hi,
I’ve searched to understand how this works but it’s not clear to me:
let us assume I have a simulation box that is smaller than the force cutoff for the reaxFF force field implemented via reax/c. If I use periodic boundary conditions and half of a molecule crosses the boundary and gets wrapped around, is the direction and magnitude of the force preserved correctly, crossing the boundary, or does the wrap around create unwanted forces acting across the box?
Thank you in advance,
Mattia Siviero
Hi,
I’ve searched to understand how this works but it’s not clear to me:
let us assume I have a simulation box that is smaller than the force cutoff for the reaxFF force field implemented via reax/c. If I use periodic boundary conditions and half of a molecule crosses the boundary and gets wrapped around, is the direction and magnitude of the force preserved correctly, crossing the boundary, or does the wrap around create unwanted forces acting across the box?
well, first of all. This is a trivial enough question that you could simply set up a test computation output the relevant data and see for yourself.
LAMMPS does not employ minimum image conventions. rather it does domain decomposition where atoms from neighboring subdomains (and in the case of a single processor, those are periodic copies of the system) are included as so-called ghost atoms up until the communication cutoff.
now if the box is smaller than the force cutoff, the central atom will “see” multiple periodic copies of the same atom.
take this simple illustrative input with just two atoms for example:
atom_modify map array
region box block -1 1 0 1 0 1
create_box 1 box
create_atoms 1 single 0.2 0.0 0.0
create_atoms 1 single 0.5 0.0 0.0
mass 1 1.0
pair_style lj/cut 0.4
pair_coeff 1 1 1.0 1.0
run 0
if you look at the neighbor list statistics at the end of the output you’ll see:
Nlocal: 2.00000 ave 2 max 2 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 10.0000 ave 10 max 10 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 1.00000 ave 1 max 1 min
Histogram: 1 0 0 0 0 0 0 0 0 0
there are 2 local atoms (i.e atoms in the primary subdomain), but 10 ghost atoms (the communication/ghost atom cutoff is 0.7 sigma)
but each atom has exactly one neighbor. now if you increase the cutoff to 1.0 sigma:
Nlocal: 2.00000 ave 2 max 2 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 78.0000 ave 78 max 78 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 9.00000 ave 9 max 9 min
Histogram: 1 0 0 0 0 0 0 0 0 0
there are still 2 local atoms, but now you have 78 ghost atoms and each atom has 9 neighbors, i.e. the other atom and 8 periodic copies.
thus for your system where part of the molecule is outside the box, LAMMPS records the position in the principal box (plus image flags for counting how often an atom has passed through PBC) but will always consider all periodic copies up to the communication cutoff when building the neighbor lists and take those distances when computing forces.
now back to the original input and let’s look at the forces. you need to add at the bottom:
print “atom 1: (fx[1]) atom 2: (fx[2])”
this generates the output:
atom 1: -300958488.33642864227 atom 2: 300958488.33642864227
now I am creating the atoms so they are across periodic boundaries but have the same relative position:
create_atoms 1 single 0.8 0.0 0.0
create_atoms 1 single -0.9 0.0 0.0
and you’ll see that the forces are the same.
hope that explains things.
Axel.