Fix interaction symmetry operations (#200)

When we check local symmetry operations in the dipole methods or when we generate symmetric interactions in the exchange interaction we were using the rotation matricies of the symmetry operations of the whole crystal. This is the crystal point group. In many cases the point group of the Wyckoff positions we are using is the same as the point group of the crystal. But this is not always the case.

This meant that sometimes our symmetry checks or use of symops to generate exchange interactions, would fail. Generally the fall back then has been to turn off the symmetry operations and avoid checks and write interactions by hand.

The correct solution is actually to calculate the local point group of each site in the unit cell motif and then use these to check/generate interactions. The local point group can be found by considering the rotational symmetry operations of the crystal which leave that point invariant.

This is now implemented and the local point group operations can be returned with atom_motif_local_point_group_symops().

src/jams/core/interactions.cc

@@ -25,7 +25,7 @@ namespace { //anon
 
       for (auto const &J : interactions) {
         auto new_J = J;
-        auto symmetric_points = globals::lattice->generate_symmetric_points(new_J.r_ij, jams::defaults::lattice_tolerance);
+        auto symmetric_points = globals::lattice->generate_symmetric_points(J.unit_cell_pos_i, new_J.r_ij, jams::defaults::lattice_tolerance);
         for (const auto p : symmetric_points) {
           new_J.r_ij = p;
           symops_interaction_data.push_back(new_J);

src/jams/core/lattice.cc

@@ -1004,15 +1004,15 @@ double Lattice::max_interaction_radius() const {
 
 // generate a vector of points which are symmetric to r_cart under the crystal symmetry
 // the tolerance is used to determine if two points are equivalent
-std::vector<Vec3> Lattice::generate_symmetric_points(const Vec3 &r_cart, const double &tolerance = jams::defaults::lattice_tolerance) const {
+std::vector<Vec3> Lattice::generate_symmetric_points(const int motif_position, const Vec3 &r_cart, const double &tolerance = jams::defaults::lattice_tolerance) {
 
   const auto r_frac = cartesian_to_fractional(r_cart);
   std::vector<Vec3> symmetric_points;
 
   // store the original point
   symmetric_points.push_back(r_cart);
   // loop through all of the symmmetry operations
-  for (const auto rotation_matrix : rotations_) {
+  for (const auto rotation_matrix : atom_motif_local_point_group_symops(motif_position)) {
     // apply a symmetry operation
     const auto r_sym = fractional_to_cartesian(rotation_matrix * r_frac);
 
@@ -1026,13 +1026,13 @@ std::vector<Vec3> Lattice::generate_symmetric_points(const Vec3 &r_cart, const d
   return symmetric_points;
 }
 
-bool Lattice::is_a_symmetry_complete_set(const std::vector<Vec3> &points, const double &tolerance = jams::defaults::lattice_tolerance) const {
+bool Lattice::is_a_symmetry_complete_set(const int motif_position, const std::vector<Vec3> &points, const double &tolerance = jams::defaults::lattice_tolerance) {
   // loop over the collection of points
   for (const auto r : points) {
     // for each point generate the symmetric points according to the the crystal symmetry
-    for (const auto r_sym : generate_symmetric_points(r, tolerance)) {
+    for (const auto r_sym : generate_symmetric_points(motif_position, r, tolerance)) {
       // if a symmetry generated point is not in our original collection of points then our original collection was not a complete set
-      // and we return fals
+      // and we return false
       if (!vec_exists_in_container(points, r_sym, tolerance)) {
         return false;
       }
@@ -1116,6 +1116,23 @@ const std::vector<Vec3> &Lattice::atom_cartesian_positions() const {
   return cartesian_positions_;
 }
 
+const std::vector<Mat3> &Lattice::atom_motif_local_point_group_symops(const int &i) {
+    // Pre-calculate the symops the first time the function is called
+    if (motif_local_pg_symops_.empty()) {
+        motif_local_pg_symops_.resize(num_motif_atoms());
+        for (auto m = 0; m < num_motif_atoms(); ++m) {
+            auto motif_position = motif_atom(m).position;
+            for (auto symop : rotations_) {
+                auto new_position = symop * motif_position;
+                if  (approximately_equal(motif_position, new_position, jams::defaults::lattice_tolerance)) {
+                    motif_local_pg_symops_[m].push_back(symop);
+                }
+            }
+        }
+    }
+    return motif_local_pg_symops_[i];
+}
+
 double jams::maximum_interaction_length(const Vec3 &a, const Vec3 &b, const Vec3 &c, const Vec3b& periodic_boundaries) {
   // 3D periodic
   // -----------

src/jams/core/lattice.h

@@ -77,6 +77,7 @@ class Lattice : public Base {
     const std::vector<Vec3>& atom_cartesian_positions() const;
 
     unsigned atom_motif_position(const int &i) const;   // integer index within the motif
+    const std::vector<Mat3>& atom_motif_local_point_group_symops(const int &i); // return a vector of the point group operations local to the site
 
     int atom_unitcell(const int &i) const;
 
@@ -112,13 +113,16 @@ class Lattice : public Base {
 
     Vec3 generate_image_position(const Vec3 &unit_cell_cart_pos, const Vec3i &image_vector) const;
 
-    std::vector<Vec3> generate_symmetric_points(const Vec3 &r, const double &tolerance) const;
+    // Generates a list of points symmetric to r_cart based on the local point group symmetry of motif_position
+    std::vector<Vec3> generate_symmetric_points(const int motif_position, const Vec3 &r_cart, const double &tolerance);
 
     Vec3 cartesian_to_fractional(const Vec3 &r_cart) const;
 
     Vec3 fractional_to_cartesian(const Vec3 &r_frac) const;
 
-    bool is_a_symmetry_complete_set(const std::vector<Vec3> &points, const double &tolerance) const;
+    // Returns true if the points are a symmetry complete set (the symmetry operations on each point, generate
+    // a point in the set), based on the local point group operations of the given motif position.
+    bool is_a_symmetry_complete_set(int motif_position, const std::vector<Vec3> &points, const double &tolerance);
 
     // lookup the site index but unit cell integer coordinates and motif offset
     int site_index_by_unit_cell(const int &i, const int &j, const int &k, const int &m) const;
@@ -184,6 +188,8 @@ class Lattice : public Base {
 
     SpglibDataset *spglib_dataset_ = nullptr;
     std::vector<Mat3> rotations_;
+    std::vector<std::vector<Mat3>> motif_local_pg_symops_;
+
 
 };
 

src/jams/hamiltonian/cuda_dipole_fft.cu

@@ -142,7 +142,7 @@ CudaDipoleFFTHamiltonian::CudaDipoleFFTHamiltonian(const libconfig::Setting &set
 
     }
       if (check_symmetry_ && (globals::lattice->is_periodic(0) && globals::lattice->is_periodic(1) && globals::lattice->is_periodic(2))) {
-        if (!globals::lattice->is_a_symmetry_complete_set(generated_positions, distance_tolerance_)) {
+        if (!globals::lattice->is_a_symmetry_complete_set(pos_i, generated_positions, distance_tolerance_)) {
           throw std::runtime_error("The points included in the dipole tensor do not form set of all symmetric points.\n"
                                    "This can happen if the r_cutoff just misses a point because of floating point arithmetic"
                                    "Check that the lattice vectors are specified to enough precision or increase r_cutoff by a very small amount.");

src/jams/hamiltonian/dipole_fft.cc

@@ -82,7 +82,7 @@ DipoleFFTHamiltonian::DipoleFFTHamiltonian(const libconfig::Setting &settings, c
         kspace_tensors_[pos_i].push_back(generate_kspace_dipole_tensor(pos_i, pos_j, generated_positions));
       }
       if (check_symmetry_ && (globals::lattice->is_periodic(0) && globals::lattice->is_periodic(1) && globals::lattice->is_periodic(2))) {
-        if (!globals::lattice->is_a_symmetry_complete_set(generated_positions, r_distance_tolerance_)) {
+        if (!globals::lattice->is_a_symmetry_complete_set(pos_i, generated_positions, r_distance_tolerance_)) {
           throw std::runtime_error("The points included in the dipole tensor do not form set of all symmetric points.\n"
                                    "This can happen if the r_cutoff just misses a point because of floating point arithmetic"
                                    "Check that the lattice vectors are specified to enough precision or increase r_cutoff by a very small amount.");