Constructing potential energy surface of a periodic system, how to deal with atoms moving out of the supercell?

Question

I'm using NVE ensemble for doing molecular dynamics (MD) simulations. The periodic system I'm working with contains $50$ atoms and the "calculator" for calculating potential energies and forces is a trained neural network. The platform I'm using for MD simulation is ase (atomic simulation environment).

After I ran the MD simulation for about $3$ps under $380$K, the atoms are moving out of the supercell. The ase module for MD simulation just faithfully updates the positions of atoms in the supercell (the supercell contains 50 atoms) based on the velocity. From its source code, it seems the module does not take periodicity of the system into account.

EDIT

@LonelyProf 's helpful comment reminds me periodic boundary conditions (PBC) should be implemented in the neural network calculator. The neural network calculator consists of two parts:

1. Generate Behler-Parrinello fingerprint given the Cartesian coordinates of atoms in the structure model. The fingerprint is actually a vector of symmetry functions that describes two-body and three-body interactions. Two-body symmetry functions are functions of atomic spacing, three-body symmetry functions are functions of both atomic spacing and angle.

   In my original implementation (decorate the neural network potential to be a readable calculator by ase), there is not a module that adjusts the Cartesian coordinates of atoms "run wild" so as to move them back to the supercell.

2. Calculate the energy of a structure given its fingerprints. The calculated fingerprint of a structure in step one is treated as the neural network input to calculate the energy of the structure.

When constructing the potential energy surface of a periodic system whose structure model has a fixed supercell, should I always keep the atoms within the supercell?

Answer 1

I think that it is best to consider the calculation of the potential energy (respecting the periodic boundary conditions) as a separate matter from keeping the atoms within the supercell. For simplicity I'm going to assume a cuboidal box with sides $(L_x,L_y,L_z)$, with periodicity in all three Cartesian directions.

Some people like to put the atoms back into the primary cell whenever they move out. This is done by replacing the atomic position vector $(x,y,z)$ by $(x,y,z)+(n_xL_x,n_yL_y,n_zL_z)$, where $(n_x,n_y,n_z)$ are integers, chosen to make the new coordinate values lie in the desired range, often $(\pm\frac{1}{2}L_x,\pm\frac{1}{2}L_y,\pm\frac{1}{2}L_z)$. If you do this every MD step, the integers will usually just need to take the values $0,\pm1$. Other people prefer not to apply this correction, and just keep track of the copy of an atom that diffuses into the periodic-image supercells.

Even if you do this correction to the atom positions, you still need to apply a similar correction every time you calculate a vector between atoms, for the purposes of computing two-body, three-body, potentials and forces. So, a properly written calculator for these quantities will not rely on the atoms lying within the same supercell. If the atomic coordinates are allowed to evolve freely, without periodic boundary corrections, then the integers $(n_x,n_y,n_z)$ will take whatever values are needed to bring the corrected interatomic vector into "minimum image" form, with components in the range $(\pm\frac{1}{2}L_x,\pm\frac{1}{2}L_y,\pm\frac{1}{2}L_z)$. There are standard prescriptions for doing this.

Often, the energy/force calculation is speeded up by a neighbour list calculation, which involves assigning each atom to a subcell of the original supercell. On the face of it, this seems to require that we put the atoms back into the supercell, before doing the assignment. But actually, a properly written function to identify nearest neighbours will have this mapping built in to it, and may just work with copies of the original coordinates. So, again, there should be no need to put the atoms back into the original supercell "by hand", if the neighbour identification is done this way, and if the vectors between neighbours are corrected as just described.

My impression, from a cursory look at ase, is that this mapping is indeed built in to the functions used to identify neighbours "within range". There is a simple Lennard-Jones calculator, for instance, which seems to rely on the neighbour list function for this purpose. The relevant lines are

        neighbors, offsets = self.nl.get_neighbors(a1)
        cells = np.dot(offsets, cell)
        d = positions[neighbors] + cells - positions[a1]

where I think that offsets contains the integers $(n_x,n_y,n_z)$ and I think that cell contains the cell lengths $(L_x,L_y,L_z)$ in the form of a $3\times3$ matrix (just the diagonal elements would be nonzero in this case).

If you are not using such a function to locate neighbours "within range" and pre-calculate the necessary minimum image corrections to the interatomic vectors, you will need to put it into your calculation of pair and three-body energies and forces, somewhere, yourself. As I mentioned, this is not the same as putting the atoms back into the supercell, when they move out.

Caveat! I am not an expert in ase, and if you want specific advice on that package, you should look at the package documentation, and at simple examples of calculators provided within it. Generally speaking, I believe that advice on debugging programs, and on specific packages, is considered off-topic for this site, so I think that's as far as I should go in any case.