Symmetric Molecular Dynamics.

We derive a formulation of molecular dynamics that generates only symmetric configurations. We implement it for all 2D planar and 3D space groups. An atlas of 2D Lennard-Jones crystals under all planar groups is created with symmetric molecular dynamics.

Molecular dynamics has long been proposed as a method for predicting or understanding crystal structures.[1] However, any practitioner will confess it is near impossible to observe point group symmetries in molecular dynamics. Here, we derive a constraint formulation of molecular dynamics where the symmetry group is an input. There is a finite number of symmetry groups. We simply simulate under all symmetry groups to generate symmetric structures.

There are two key ideas to our formulation that correspond to the two components of a space group: the point group symmetry and the Bravais lattice. The point group symmetry is treated as a holonomic constraint. The constraint equation is a function of positions that is zero when the positions are symmetric. Holonomic constraints are a relatively solved problem, and we follow previous approaches.[2−7] The Bravais lattice is a constraint on the simulation lattice vectors that ensures the point group will tile space. Namely, the Bravais lattice specifies the relative lattice vector magnitudes and directions. We ensure our simulations are consistent by working in an unconstrained lattice vector space that is mapped to the correct Bravais lattice via a precomputed tensor. This frees us to use any NPT method in the unconstrained lattice vector space while still matching the Bravais lattice.

The concept of directly simulating under a symmetry group is unknown to us. The closest examples are methods like symmetry restraints.[8] These harmonic restraints generally keep the system close to symmetric, but unlike the method we propose here, no single configuration is actually symmetric. Symmetry has certainly been considered as a measure of molecular configurations. For example, Zabrodsky et al.[9] proposed a continuous symmetry measure, which is used to quantify the symmetry of atoms. This has been used to directly optimize Lennard-Jones clusters with symmetry.[10] Of course, the direct use of symmetry for crystal structure prediction with Monte Carlo is common,[11,12] and generative models with explicitly included symmetry are common.[13] There are no molecular dynamics methods though that can directly sample space groups, which would be useful for crystal structure prediction and modeling biological assemblies.[14] Symmetric molecular dynamics may also be viewed as a special case of objective molecular dynamics, which is a general method that encompasses any infinite or finite periodic tiling of a simulation.[15,16] Similarly, others have explored generalizing periodic boundary conditions to other tilings.[17−20]

Below we derive our equations of motion and discuss implementation details. To assess the method, we show that it conserves energy and is capable of working in arbitrary space groups. Then we demonstrate its use to enumerate crystal structures of the Lennard-Jones potential under all planar groups with NPT simulations.

Theory

Equation of Motion

Consider the dynamics of N indistinguishable particles in D dimensions under a Hamiltonian H(p(t), q(t)). We wish to constrain H so that q(t) is symmetric at all times. Symmetry is a property of q(t) and a specific symmetry group of position transformations G, like mirrors along the x axis. q(t) is point group symmetric if applying any element of the group results in no change to the positions (ignoring ordering of particles)where g· means applying the group element to each particle individually, ∼ means row equivalence, and G is a finite group. Group elements are represented as affine matrices in space and planar groups.

Eq may hold trivially. For example, all particles are at the origin. Such special positions that are invariant to group elements are known as special Wyckoff positions.[21] We remove this assumption in Section IIC, but for now additionally assumewhere I is the identity transformation.

Assuming eqs and 2 hold at t = 0, the particles can be partitioned into N/|G| = n group orbits. A group orbit is the set generated by applying all elements of group G to positions q(t)

One member of all orbits will be q(t) itself, because G contains the identity element. We can label the particles as q(t) where i indicates the orbit and j indicates the group element. In crystallography, the q particles are called the asymmetric unit. We can satisfy eq at all t by specifying the following holonomic constraint:

There are |G| – 1 of these constraints per group orbit, and each removes D degrees of freedom. This means the degrees of freedom of the dynamics is D × (n – 1). We can simulate dynamics under the holonomic constraints by simply only modeling the asymmetric unit—they are the generalized coordinates.[22]

Thus, our algorithm is to only integrate the asymmetric unit and explicitly consider the remaining (N – n) particles only when computing forces. This is similar to Dayal et al.[15] Practically this is done by setting these constrained particles’ positions just before computing forces. Similar to work on periodic boundary conditions, these equations of motion may lead to linear momentum conservation problems.[23,24]

One feature of nearly all potentials used in molecular dynamics is that they are G-invariant, where G is any planar, space, or permutation group: U(g·q) = U(q). That makes the forces, F(q), G-equivariantbecause the potentials are composed of angles and distances, which are invariant to rotations, mirrors, and permutations.[25,26] For a pairwise potential, we can use eqs and 5 to rewrite the potential aswhere the |G| factor accounts for intragroup orbit interactions that are not explicitly computed. This translates an algorithm of an outer loop over the asymmetric unit and an inner loop over all particles.

Bravais Lattice

A space group consists of both a point group and a Bravais lattice. The Bravais lattice is specified with D D-dimensional unit cell vectors. Particles always remain in one cell among the lattice cells, which are called images. For example, we could simulate the “root” cell and its 26 neighboring cells in 3 dimensions. We follow the approach above and treat each image of the system with virtual particles while only integrating the root cell. This means all images of the system are explicit, and we can violate the minimum image convention. We were not signatories of the minimum image convention anyway. This approach allows the cell vectors to shrink well below the distance cutoff of the potential, provided we have enough virtual particles to populate past the cutoff of the asymmetric unit of the origin cell. You can simulate 3 images to allow the cells to shrink to at least 1/a the cutoff distance.

We need to convert between the fractional coordinates, which are used to tile the particles and apply the point group symmetry, to the Cartesian coordinates, which are used for integration and computed potentials. Given the box vectors in row-form B, we can transform between the representations viawhere s(t) is the fraction of each lattice vector (i.e., fractional coordinates). Wrapping is trivial with fractional coordinates: s(t) fmod 1.0 will wrap the coordinates. All point group transformations are applied in s(t); however, a B–1 term should be added to eq so that it operates on fractional coordinates.

Bravais lattices include more than just the usual cubic and triclinic lattices commonly seen in molecular dynamics barostats. To ensure the cell vectors are consistent with the Bravais lattice while changing box size, we define a tensor L of shape D × D × D × D that maps from a triclinic box vector to the proper Bravais lattice box vectors of the space group. For example, L2011 is the contribution to Bravais lattice vector 2’s x component from triclinic box vector 1’s y component. There are many choices that could be made for L. For example, to make a cubic Bravais lattice from a triclinic box vector, we require a single parameter a to define the three lattice vectors (a, 0, 0),(0, a, 0),(0, 0, a). We could set a by averaging all the vector lengths, averaging all vector components, or selecting a to be the first element of the first vector. Each of these choices gives a different L, and some have large null spaces. NPT is then accomplished via scaling Monte Carlo moves in the triclinic box vectors (B′) following Frenkel and Eppenga,[27] and the proper Bravais lattice is computed via B = LB′.

Wyckoff Positions

It is possible to have particles that violate eq while still satisfying eq if q is in a special position called a Wyckoff position—like the origin.[21] To perform constrained molecular dynamics of particle q(t) occupying a Wyckoff position, we define a subgroup G′ that contains the elements of G which do not leave q(t) invariant plus an identity group element. The identity of this subgroup is not the identity transform but instead a transform that projects from a general position into the Wyckoff position. For example, the Wyckoff position may be the vertical line x = 0, and the identity group element would be the transform x′ = 0, y′ = y. We will denote this group element as P to hint it is a projection.

The group orbit is similarly defined on the subgroup, and the other procedures above apply. However, q(t) must stay in a Wyckoff position at all times to satisfy eq . This can be accomplished via traditional constrained molecular dynamics of Lagrange multipliers.[28] Omitting the indices on q(t), our holonomic constraint isand the force from the constraint will bewhere J[σ] is the Jacobian of σ with respect to constraint dimension and element of q(t). We can solve for λ by knowing that σ[q(t + Δt)] = 0where Δt is the time step, m is the mass of the particle, and q′(t + Δt) is q(t) integrated without the constraint force by Δt. All terms are constant except σ[q′(t + Δt)], which simplifies computation.

Methods

We use the BAOAB Langevin dynamics integrator described in refs (29 and 30). Eq is applied during position updates, and eq is applied before velocity updates (force computation). Degrees of freedom is computed from number of asymmetric unit particles and deducted degrees of freedom from Wyckoff position restraints. All simulations are Lennard-Jones potentials with cutoff 3.5 and in reduced units. NVE simulations are conducted with the velocity-verlet integrator. A time step of 0.005 and a Langvin γ of 0.1 were used for all simulations. Since images are explicit in our implementation, it is necessary to specify the number. We use an image radius of 2–meaning 32 images are simulated where D is the dimension. To generate starting configurations, points were randomly generated and filtered to fit into the space group asymmetric unit as specified by Aroyo.[31] Point group generators and Wyckoff sites were taken from the Bilbao crystallography server.[32−34]

We define our results in reduced units, as defined in ref (35). Specifically, energy (ϵ) is fundamental, and τ is a derived unit of the formwhere L, m, and ϵ are the fundamental units of length, mass, and energy, respectively.

Results

We first consider if our implementation conserves energy. Figure shows the total energy of NVE simulations under a subset of space groups with 5 particles in the asymmetric unit. These were done at number densities of 0.2, with a starting temperature 0.5, and for 30k timesteps. The bottom trace (P1) has no symmetry constraints and shows good conservation. There are more fluctuations at other symmetry groups because there are more particles in their unit cells and thus higher energy fluctuations. For example, space group 127 (P4/mbm) has 80 particles in a unit cell when there are 5 in the asymmetric unit, meaning the interaction potential felt has more particles contributing to it.

Figure 1

Total energy from NVE simulations under different symmetry groups in 3D. Groups are indicated with Hall numbers. Four particles are in the asymmetric unit, and the simulations are at a number density of 0.2 and a starting temperature of 0.5. The increase in fluctuations is because the unit cell (total particles) increases with the size of the symmetry group.

Figure shows an enumeration of crystal structures under different symmetry groups for a 2D Lennard-Jones fluid. The structures are generated in 2 steps. First, we simulate under a symmetry group constraint in NPT (P = 0.25, T = 0.1) for 1 M steps. Next, we do a constrained equilibration under NVT for 100k steps at T = 0.05. This structure is then the proposed crystal structure for the given symmetry group. Figure shows the root mean square deviation (RMSD) if the resultant structure is simulated under no symmetry constraint in NVE for 5k steps. The assumption is that if the structure does not collapse (RMSD rise), it is metastable. We indeed find that this protocol under no symmetry constraints (p1) gives the correct hexagonal packing.

Figure 2

An atlas of 258 Lennard-Jones crystal structures. Subplots are broken down by planar group, and each line color/style indicates Wyckoff occupancy and unit cell particle count. Individual plots show RMSD from the starting configuration, which was constrained to match the planar group in subplot titles. RMSD was calculated during the NVE simulation of 5k steps with no symmetry constraints. A rise indicates change of configuration–meaning the starting configuration was not stable. Missing traces indicate either simulation diverged or number density did not exceed 0.5.

To enumerate all planar groups in 2D, we simulate under each group, with 4 choices of particle number (in unit cell) and varying occupancy of Wyckoff sites. As expected, the planar groups with hexagonal Bravais lattices (or permit them) have stable structures: p1, p2, pg, p3, p6, and cmm.[36] Some unusual stable structures are seen without hexagonal close-packing like p4m and p3m1. Metastable structures like these would be nearly impossible to generate without symmetry constraints in molecular dynamics. We find some symmetries have no metastable structures: p6m and most square close-packing (cm, cmm, p4g). Interestingly, voids seem to be the way to stabilize these close-packing structures like in pmg.

Discussion

Here is our advice on implementing symmetric molecular dynamics in a modern molecular dynamics engine. The asymmetric unit should be integrated as usual. Make the nonasymmetric unit particles (images) be ghost particles; ghost particles are nonintegrated particles used in force-field calculations. The ghost particles’ positions should be set using affine matrices defining the group transformations in fractional coordinates. These matrices can be obtained from our library or crystallography tables.

Pressure computed for the asymmetric unit is not meaningful, and NPT should be done using the algorithm described above that does not require computing pressure from a virial. The lattice vectors may be stored separately than the usual lattice vectors and are only used to set the ghost particle positions. The tensor transforms can be loaded from our library. Periodic boundary conditions should be disabled entirely if doing NPT.

The constraints for Wyckoff sites are implemented as Lagrange multiplier constraints. The terms can be computed analytically at each step using eq .

Conclusions

We have formulated a symmetric molecular dynamics algorithm and implemented it. Results show that it can do NPT to enumerate metastable crystal structures. A reference implementation is available at https://github.com/whitead/symd.