Bonds crossing periodic boundaries

Hi all,

I know that the following has been discussed in the past but it still confuses me. I am modelling a zeolite with a force field that has harmonic springs for bonds and angles (e.g. a CHARMM type of a FF). Inevitably due to periodicity I will have bonds crossing periodic boundaries and, if I understood the previous debate this is not a concern in LAMMPS. I also get the warning of inconsistent image flags and bond extending over more than a half of box dimension.
When I minimize the structure and run the dynamics the energy seems fine but the system shrinks. I was wondering if this could be the issue due to the bonds.
My approach for building the topology was:

  1. start from the crystal structure
  2. find all bonds between atoms that should be connected by taking periodicity into account
  3. i have explicitly assigned image flags to 0 to all atoms
    I am aware that such issue could be a matter of e.g. incorrect unit conversion but I am quite sure it is not. Also I am trying to replicate a paper where this force field is claimed to work. I also tried to play around with force constants but found nothing obvious.
    Did anyone observe a similar thing or does anyone see a flaw in my logic?

Thanks a lot.

That one is expected for bonds crossing periodic boundaries.

That means that your simulation cell is too small (and thus the choice of periodic image ambiguous) or the bond assignments are not proper. Here is a minimal example showing an infinite 1-d bond chain through periodic boundaries.

atom_style bond
region box block 0.0 4.0 0.0 4.0 0.0 4.0
create_box 1 box bond/types 1 extra/bond/per/atom 2 extra/special/per/atom 8

mass * 1.0
pair_style zero 2.0
pair_coeff * *

bond_style harmonic
bond_coeff 1 1.0 1.0

create_atoms 1 single  0.5 1.0 1.0
create_atoms 1 single  1.5 1.0 1.0
create_atoms 1 single  2.5 1.0 1.0
create_atoms 1 single  3.5 1.0 1.0

create_bonds single/bond 1 1 2
create_bonds single/bond 1 2 3
create_bonds single/bond 1 3 4
create_bonds single/bond 1 4 1

write_data chain.data
minimize 0.0 0.0 100 1000

As can be seen from the output, only the first warning shows up and there is no energy/force since all bonds are at their equilibrium position.

minimize 0.0 0.0 100 1000
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
WARNING: Inconsistent image flags (src/domain.cpp:815)
Per MPI rank memory allocation (min/avg/max) = 6.313 | 6.313 | 6.313 Mbytes
   Step          Temp          E_pair         E_mol          TotEng         Press     
         0   0              0              0              0              0            
         1   0              0              0              0              0            
Loop time of 1.3817e-05 on 1 procs for 1 steps with 4 atoms

Hi @AnzeH,

As Axel has stated, one of the message is expected and the other one problematic. It is always good to have a look at your crystal structure with a visualisation software to check the topology output from LAMMPS.

Your procedure is also very vague. On paper your steps sound nice but they barely describe what you actually do.

  • “start from the crystal structure” what structure? A CIF file? An XYZ? A single crystal cell? Several?
  • “find all bonds between atoms” How? Using Lammps? An homebrew script?
  • “i have explicitly assigned image flags to 0 to all atoms” This seems unnecessary to me or you are using weird box settings with regards to the coordinates of your atoms.

A work flow that worked fine for me until now in these situations:

  • Start with minimal unit cell, either custom of from a reference file (no bonds defined)
  • Replicate the structure the desired number of cells in each directions
  • Only then create the needed bond structure. The create_bonds many command is a good way to do it in a controlled way. If any angles or dihedral terms are needed, then this is another story.

Take a good look at Axel’s example, do not forget to use the relevant options when creating your box or reading data file (extra/bond/per/atom and extra/special/per/atom etc…).

1 Like

An important issue to note here is that you must not use the replicate command after you have created the bonds across periodic boundaries. LAMMPS does not know how to handle this case and would generate invalid topologies.

1 Like

A good exercise to understand the principles would be to take the example and extend it to add infinite chains that go through atom 1 in y- and z- direction. There should be just the warning about the image flags and no forces, if done correctly.

1 Like

I added a note stating this in my WF proposition, you are absolutely right emphasising it.

Thanks to both of you.
Yes, I started from the .cif file, replicated it using Atomsk and then used my script to find bonds and angles. It is much more clear to me now.

Best,
Anze

Hi again,

I checked Axel’s minimal example. I tried changing the region command to:

region box block 0.0 5.0 0.0 4.0 0.0 4.0

which again gives the warning about the extent of bonds. Of course in this case the energy/forces will not be zero. But, the 1-4 bond length is 2.0. As I have only bonds in this system the maximum extent should be 2.2, which is less than 5.0/2 = 2.5. Unless of course LAMMPS took the distance between 1-4 by not crossing the boundary. In that case the distance is 3.0 which would indeed make sense to print out the warning. But, by visualisation this does not seem to be the case (i dumped only atoms 1 and 4).


I also tried varying the cutoff but to no effect. Any additional comment would be appreciated.

EDIT: I just realised that going to

region box block 0.0 5.0 0.0 5.0 0.0 5.0

indeed removes the said warning.
Best,
Anze

Yes, because the check in question determines the maximum extent by taking the maximum existing bond distance, multiplies it by 1.1 and then compares that against half of the box length in x-, y-, and z-direction if they are all periodic. Thus even if half the box is 2.5 in x after your initial modification, it is still 2.0 in y and z. Only after the second modification are all box lengths larger.