Core: fixed Geometry resize for multi-atom cells.

Changed the Geometry constructor to not change the atom positions. Instead, the Configparser and Geometry API transform by the bravais_vectors as appropriate.

core/docs/c-api/Geometry.md

@@ -52,6 +52,7 @@ void Geometry_Set_Cell_Atoms(State *state, int n_atoms, float ** atoms)
 ```
 
 Set the number and positions of atoms in a basis cell.
+Positions are in units of the bravais vectors.
 
 
 

core/include/Spirit/Geometry.h

@@ -57,6 +57,7 @@ PREFIX void Geometry_Set_N_Cells(State * state, int n_cells[3]) SUFFIX;
 
 /*
 Set the number and positions of atoms in a basis cell.
+Positions are in units of the bravais vectors.
 */
 PREFIX void Geometry_Set_Cell_Atoms(State *state, int n_atoms, float ** atoms) SUFFIX;
 

core/src/Spirit/Geometry.cpp

@@ -207,7 +207,11 @@ try
     {
         for (int i=0; i<n_atoms; ++i)
         {
-            cell_atoms.push_back({atoms[i][0], atoms[i][1], atoms[i][2]});
+            Vector3 cell_atom =
+                  old_geometry.bravais_vectors[0] * atoms[i][0]
+                + old_geometry.bravais_vectors[1] * atoms[i][1]
+                + old_geometry.bravais_vectors[2] * atoms[i][2];
+            cell_atoms.push_back(cell_atom);
             iatom.push_back(i);
             if( i < old_geometry.n_cell_atoms )
             {

core/src/data/Geometry.cpp

@@ -10,6 +10,8 @@
 #include "QhullFacetList.h"
 #include "QhullVertexSet.h"
 
+#include <fmt/ostream.h>
+
 #include <random>
 #include <array>
 
@@ -25,15 +27,6 @@ namespace Data
         n_cells_total(n_cells[0] * n_cells[1] * n_cells[2]),
         pinning(pinning), defects(defects)
     {
-        // Get x,y,z of component of atom positions in unit of length (instead of in units of a,b,c)
-        for (int iatom = 0; iatom < this->n_cell_atoms; ++iatom)
-        {
-            Vector3 build_array = this->bravais_vectors[0] * this->cell_atoms[iatom][0]
-                                + this->bravais_vectors[1] * this->cell_atoms[iatom][1]
-                                + this->bravais_vectors[2] * this->cell_atoms[iatom][2];
-            this->cell_atoms[iatom] = this->lattice_constant * build_array;
-        }
-
         // Generate positions and atom types
         this->positions = vectorfield(this->nos);
         this->generatePositions();
@@ -117,8 +110,10 @@ namespace Data
                                 std::abs(diff[1]) < epsilon &&
                                 std::abs(diff[2]) < epsilon )
                             {
-                                spirit_throw(Utility::Exception_Classifier::System_not_Initialized, Utility::Log_Level::Severe, fmt::format(
-                                    "Unable to initialize Spin-System, since 2 spins occupy the same space within a margin of {}.\nPlease check the config file!", epsilon ));
+                                std::string message = fmt::format(
+                                    "Unable to initialize Spin-System, since 2 spins occupy the same space within a margin of {} at position ({}).\n"
+                                    "Please check the config file!", epsilon, cell_atoms[i].transpose());
+                                spirit_throw(Utility::Exception_Classifier::System_not_Initialized, Utility::Log_Level::Severe, message);
                             }
                         }
                     }
@@ -502,16 +497,20 @@ namespace Data
         const scalar epsilon = 1e-6;
 
         // ----- Find dimensionality of the basis -----
-        if     ( n_cell_atoms == 1 ) dims_basis = 0;
-        else if( n_cell_atoms == 2 ) dims_basis = 1;
-        else if( n_cell_atoms == 3 ) dims_basis = 2;
+        if     ( n_cell_atoms == 1 )
+            dims_basis = 0;
+        else if( n_cell_atoms == 2 )
+        {
+            dims_basis = 1;
+            test_vec_basis = cell_atoms[0] - cell_atoms[1];
+        }
         else
         {
             // Get basis atoms relative to the first atom
             Vector3 v0 = cell_atoms[0];
             std::vector<Vector3> b_vectors(n_cell_atoms-1);
             for( int i = 1; i < n_cell_atoms; ++i )
-                b_vectors[i-1] = cell_atoms[i] - v0;
+                b_vectors[i-1] = (cell_atoms[i] - v0).normalized();
 
             // Calculate basis dimensionality
             // test vec is along line

core/src/io/Configparser.cpp

@@ -443,13 +443,26 @@ namespace IO
             Log(Log_Level::Parameter, Log_Sender::IO, fmt::format("        a = {}", bravais_vectors[0].transpose()));
             Log(Log_Level::Parameter, Log_Sender::IO, fmt::format("        b = {}", bravais_vectors[1].transpose()));
             Log(Log_Level::Parameter, Log_Sender::IO, fmt::format("        c = {}", bravais_vectors[2].transpose()));
-            Log(Log_Level::Parameter, Log_Sender::IO, fmt::format("Basis: {}  atom(s) at the following positions:", n_cell_atoms));
-            if( !cell_composition.disordered )
+            Log(Log_Level::Parameter, Log_Sender::IO, fmt::format("Basis cell: {}  atom(s)", n_cell_atoms));
+            Log(Log_Level::Parameter, Log_Sender::IO, "Relative positions (first 10):");
+            for( int iatom = 0; iatom < n_cell_atoms && iatom < 10; ++iatom )
+                Log(Log_Level::Parameter, Log_Sender::IO, fmt::format("        atom {} at ({}), mu_s={}", iatom, cell_atoms[iatom].transpose(), cell_composition.mu_s[iatom]));
+
+            // Get x,y,z of component of atom positions in unit of length (instead of in units of a,b,c)
+            for( int iatom = 0; iatom < n_cell_atoms; ++iatom )
             {
-                for (int iatom = 0; iatom < n_cell_atoms; ++iatom)
-                    Log(Log_Level::Parameter, Log_Sender::IO, fmt::format("        atom {} at ({}), mu_s={}", iatom, cell_atoms[iatom].transpose(), cell_composition.mu_s[iatom]));
+                Vector3 cell_atom =
+                      bravais_vectors[0] * cell_atoms[iatom][0]
+                    + bravais_vectors[1] * cell_atoms[iatom][1]
+                    + bravais_vectors[2] * cell_atoms[iatom][2];
+                cell_atoms[iatom] = lattice_constant * cell_atom;
             }
-            else
+
+            Log(Log_Level::Parameter, Log_Sender::IO, "Absolute atom positions (first 10):", n_cell_atoms);
+            for( int iatom = 0; iatom < n_cell_atoms && iatom < 10; ++iatom )
+                Log(Log_Level::Parameter, Log_Sender::IO, fmt::format("        atom {} at ({}), mu_s={}", iatom, cell_atoms[iatom].transpose(), cell_composition.mu_s[iatom]));
+
+            if( cell_composition.disordered )
                 Log(Log_Level::Parameter, Log_Sender::IO, "Note: the lattice has some disorder!");
 
             // Defects