feat: support arbitrary lattice vectors

siqad · Feb 12, 2024 · 74a3536 · 74a3536
1 parent 89c4091
commit 74a3536
Show file tree

Hide file tree

Showing 5 changed files with 43 additions and 17 deletions.
diff --git a/libs/siqadconn b/libs/siqadconn
diff --git a/src/interface.cc b/src/interface.cc
@@ -70,11 +70,33 @@ SimParams SimAnnealInterface::loadSimParams()
 
   SimParams sp;
 
+  // find the lattice layer in the layer list and get the lattice vector
+  bool valid_lat_vec_found = false;
+  for (auto layer : sqconn->getLayers()) {
+    log.debug() << "Layer name: " << layer.name << ", layer type: " << layer.type << std::endl;
+    if (layer.type.compare("Lattice") == 0) {
+      if (!layer.lat_vec.isValid()) {
+        continue;
+      }
+      valid_lat_vec_found = true;
+      sp.lat_vec.a1 = layer.lat_vec.a1;
+      sp.lat_vec.a2 = layer.lat_vec.a2;
+      sp.lat_vec.atoms = layer.lat_vec.atoms;
+      log.debug() << "Lattice vector found: " << layer.lat_vec.name << std::endl;
+      log.debug() << "Lattice vector: a1=(" << sp.lat_vec.a1.first << "," << sp.lat_vec.a1.second << ")"
+          << ", a2=(" << sp.lat_vec.a2.first << "," << sp.lat_vec.a2.second << ")" << std::endl;
+      break;
+    }
+  }
+  if (!valid_lat_vec_found || !sp.lat_vec.isValid()) {
+    throw std::runtime_error("No valid lattice vector found.");
+  }
+
   // grab all physical locations
   log.debug() << "Grab all physical locations..." << std::endl;
   std::vector<EuclCoord> db_locs;
   for(auto db : *(sqconn->dbCollection())) {
-    db_locs.push_back(SimParams::latToEuclCoord(db->n, db->m, db->l));
+    db_locs.push_back(SimParams::latToEuclCoord(db->n, db->m, db->l, sp.lat_vec));
     log.debug() << "DB loc: x=" << db_locs.back().first
         << ", y=" << db_locs.back().second << std::endl;
   }
@@ -104,7 +126,7 @@ SimParams SimAnnealInterface::loadSimParams()
     } else if (defect->has_lat_coord) {
       log.debug() << "**** HAS LATCOORD ****" << std::endl;
       // TODO find center point and calculate that x y z
-      EuclCoord anchor_tl = SimParams::latToEuclCoord(defect->n, defect->m, defect->l);
+      EuclCoord anchor_tl = SimParams::latToEuclCoord(defect->n, defect->m, defect->l, sp.lat_vec);
       log.debug() << "tl -- n=" << defect->n << ", m=" << defect->m << ", l=" << defect->l << std::endl;
       log.debug() << "tl -- x=" << anchor_tl.first << ", y=" << anchor_tl.second << std::endl;
       if (defect->w == 1 && defect->h == 1) {
@@ -114,7 +136,7 @@ SimParams SimAnnealInterface::loadSimParams()
         int m_br = defect->m + (int) floor((defect->h - defect->l)/ 2);
         int l_br = (defect->h - defect->l) % 2;
         log.debug() << "br -- n=" << n_br << ", m=" << m_br << ", l=" << l_br << std::endl;
-        EuclCoord anchor_br = SimParams::latToEuclCoord(n_br, m_br, l_br);
+        EuclCoord anchor_br = SimParams::latToEuclCoord(n_br, m_br, l_br, sp.lat_vec);
         log.debug() << "br -- x=" << anchor_br.first << ", y=" << anchor_br.second << std::endl;
         FPType x_mean = (anchor_br.first + anchor_tl.first) / 2;
         FPType y_mean = (anchor_br.second + anchor_tl.second)  / 2;

diff --git a/src/simanneal.cc b/src/simanneal.cc
@@ -51,20 +51,20 @@ void SimParams::setDBLocs(const std::vector<EuclCoord> &t_db_locs)
   v_fc.resize(n_dbs);
 }
 
-void SimParams::setDBLocs(const std::vector<LatCoord> &t_db_locs)
+void SimParams::setDBLocs(const std::vector<LatCoord> &t_db_locs, const LatticeVector &lat_vec)
 {
   std::vector<EuclCoord> db_locs;
   for (LatCoord lat_coord : t_db_locs) {
     assert(lat_coord.size() == 3);
-    db_locs.push_back(latToEuclCoord(lat_coord[0], lat_coord[1], lat_coord[2]));
+    db_locs.push_back(latToEuclCoord(lat_coord[0], lat_coord[1], lat_coord[2], lat_vec));
   }
   setDBLocs(db_locs);
 }
 
-EuclCoord SimParams::latToEuclCoord(const int &n, const int &m, const int &l)
+EuclCoord SimParams::latToEuclCoord(const int &n, const int &m, const int &l, const LatticeVector &lat_vec)
 {
-  FPType x = n * constants::lat_a;
-  FPType y = m * constants::lat_b + l * constants::lat_c;
+  FPType x = n * (lat_vec.a1.first + lat_vec.a2.first) + lat_vec.atoms[l].first;
+  FPType y = m * (lat_vec.a1.second + lat_vec.a2.second) + lat_vec.atoms[l].second;
   return std::make_pair(x, y);
 }
 

diff --git a/src/simanneal.h b/src/simanneal.h
@@ -17,6 +17,7 @@
 #include <cmath>
 #include <mutex>
 #include <thread>
+#include "libs/siqadconn/src/siqadconn.h"
 
 //#include <boost/thread.hpp>
 #include <boost/random.hpp>
@@ -26,11 +27,6 @@
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 
 namespace constants{
-  // lattice
-  const FPType lat_a = 3.84;  // lattice vector in x, angstroms (intra dimer row)
-  const FPType lat_b = 7.68;  // lattice vector in y, angstroms (inter dimer row)
-  const FPType lat_c = 2.25;  // dimer pair separation, angstroms
-
   // energy band diagram
   const FPType eta = 0.59;    // TODO enter true value; energy difference between (0/-) and (+/0) levels
 
@@ -75,9 +71,9 @@ namespace phys {
 
     void setDBLocs(const std::vector<EuclCoord> &);
 
-    void setDBLocs(const std::vector<LatCoord> &);
+    void setDBLocs(const std::vector<LatCoord> &, const LatticeVector &);
 
-    static EuclCoord latToEuclCoord(const int &n, const int &m, const int &l);
+    static EuclCoord latToEuclCoord(const int &n, const int &m, const int &l, const LatticeVector &lat_unit_cell);
 
     // Clear old fixed charges and set new ones
     void setFixedCharges(const std::vector<EuclCoord3d> &t_fc_locs, 
@@ -128,6 +124,7 @@ namespace phys {
     FPType eps_r=5.6;        // Relative premittivity on the surface
     FPType debye_length=5.0;     // Debye Length (nm)
     int n_dbs;                   // Number of DBs in the simulation
+    LatticeVector lat_vec;
     std::vector<std::pair<FPType,FPType>> db_locs;  // Location of DBs
     ublas::matrix<FPType> db_r;  // Matrix of distances between all DBs
     ublas::matrix<FPType> v_ij;  // Matrix of coulombic repulsion between occupied DBs

diff --git a/tests/simanneal.test.cc b/tests/simanneal.test.cc
@@ -5,6 +5,13 @@ namespace ublas = boost::numeric::ublas;
 
 TEST_CASE( "OR_mod 00 test" ) {
     auto sim_params = phys::SimParams();
+    // 100 lat vec
+    sim_params.lat_vec.a1 = std::pair<FPType, FPType>(3.84, 0);
+    sim_params.lat_vec.a2 = std::pair<FPType, FPType>(0, 7.68);
+    sim_params.lat_vec.atoms = std::vector<std::pair<float, float>>({
+        {0, 0},
+        {0, 2.25}
+    });
     sim_params.setDBLocs({
         {-4,-2, 0},
         {-2,-1, 0},
@@ -13,7 +20,7 @@ TEST_CASE( "OR_mod 00 test" ) {
         { 0, 1, 0},
         { 0, 2, 1},
         { 0, 4, 1}
-    });
+    }, sim_params.lat_vec);
     auto annealer = phys::SimAnneal(sim_params);
     annealer.invokeSimAnneal();
     auto results = annealer.suggestedConfigResults(true);