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 1 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 assume

where I is the identity transformation.

Assuming eqs 1 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 i (t)

One member of all orbits will be q i (t) itself, because G contains the identity element. We can label the particles as q ij (t) where i indicates the orbit and j indicates the group element. In crystallography, the q i0 particles are called the asymmetric unit. We can satisfy eq 1 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 (Nn) 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-equivariant

because 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 3 and 5 to rewrite the potential as

Uu G uGu qq q qq qq () ( ) 2 () ( ) ij ij i n j G ii j i n ki n j G ik j 01 0 01 0 0 ∑ ∑∑ ∑ ∑ ∑ =-= || -+ | | -< == || == += || (6)

where 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. B. Bravais Lattice. A space group consists of both a point group and a Bravais lattice. The Bravais lattice is specified with DD -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 aD 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 via

where 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 3 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, L 2011 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′.

C. Wyckoff Positions. It is possible to have particles that violate eq 2 while still satisfying eq 1 if q i0 is in a special position called a Wyckoff positionlike the origin. 21 To perform constrained molecular dynamics of particle q i0 (t) occupying a Wyckoff position, we define a subgroup G′ that contains the elements of G which do not leave q i0 (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 i0 (t) must stay in a Wyckoff position at all times to satisfy eq 1. This can be accomplished via traditional constrained molecular dynamics of Lagrange multipliers. 28 Omitting the indices on q i0 (t), our holonomic constraint is

and the force from the constraint will be

where 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)] = 0

where Δt is the time step, m is the mass of the particle, and q′(t + Δt)isq(t) integrated without the constraint force by Δt.

All terms are constant except σ[q′(t + Δt)], which simplifies computation.

We use the BAOAB Langevin dynamics integrator described in refs 29 and 30. Eq 4 is applied during position updates, and eq 3 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 3 2D 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][33][34] We define our results in reduced units, as defined in ref 35. Specifically, energy (ϵ) is fundamental, and τ is a derived unit of the form

where L, m, and ϵ are the fundamental units of length, mass, and energy, respectively.

We first consider if our implementation conserves energy. Figure 1 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 2 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 2 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.

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.

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 10.

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.