Build HistCustomModes with strain and then displacement.

Test that buildHist -> analyze gives back the amplitudes generated randomly (Tested only at Gamma for 1 1 1 cell)
piti-diablotin · Mar 12, 2021 · 62d5ed4 · 62d5ed4
1 parent 131b8b6
commit 62d5ed4
Show file tree

Hide file tree

Showing 4 changed files with 114 additions and 11 deletions.
diff --git a/include/hist/histcustommodes.hpp b/include/hist/histcustommodes.hpp
@@ -287,6 +287,21 @@ class HistCustomModes : public HistDataDtset {
      */
     void setInstableModes(const InstableModes instableModes);
 
+    /**
+     * @brief Get the amplitudes of the displacements that can be condensed 
+     * at time itime
+     * @param itime the time of the HIST
+     * @return A map of qpt / vector of qMode informations (number,amplitude,energy)
+     */
+    const DispDB::qptTree& getDispAmplitudes(const int itime);
+
+    /**
+     * @brief Get the strain matrix used to distord the lattice at itime 
+     * @param itime the time of the HIST
+     * @return the strain matrix
+     */
+    const geometry::mat3d& getStrainMatrix(const int itime);
+
 };
 
 #endif  // HISTCUSTOMMODES_HPP
diff --git a/src/hist/histcustommodes.cpp b/src/hist/histcustommodes.cpp
@@ -581,6 +581,11 @@ void HistCustomModes::push(const Dtset& dtset)
 }
 
 void HistCustomModes::buildInsert(Supercell& supercell, const unsigned itime) {
+  // First apply strain, then displacement so that
+  // X = (1+eta)(R+tau) + u
+  if ( itime < _strainDist.size() ) {
+    supercell.applyStrain(_strainDist[itime]);
+  }
   if ( itime < _condensedModes.size() ) {
     if (_db.natom()!=_reference.natom())
       throw EXCEPTION("natoms are different in DB and reference structure",ERRDIV);
@@ -590,9 +595,6 @@ void HistCustomModes::buildInsert(Supercell& supercell, const unsigned itime) {
       }
     }
   }
-  if ( itime < _strainDist.size() ) {
-    supercell.applyStrain(_strainDist[itime]);
-  }
   if ( itime >= _ntime )
     throw EXCEPTION("itime larger than _ntime.",ERRWAR);
   this->insert(itime,supercell);
@@ -620,3 +622,15 @@ void HistCustomModes::initRandomEngine() {
   }
   _randomEngine.seed(seed);
 }
+
+const DispDB::qptTree& HistCustomModes::getDispAmplitudes(const int itime) {
+  if ( itime >= _condensedModes.size() ) 
+    throw EXCEPTION("Bad value for itime",ERRDIV);
+  return _condensedModes[itime];
+}
+
+const geometry::mat3d& HistCustomModes::getStrainMatrix(const int itime) {
+  if ( itime >= _strainDist.size() ) 
+    throw EXCEPTION("Bad value for itime",ERRDIV);
+  return _strainDist[itime];
+}
diff --git a/src/io/dtset.cpp b/src/io/dtset.cpp
@@ -1075,4 +1075,5 @@ void Dtset::applyStrain(const geometry::mat3d& strainMatrix) {
   const mat3d buffRprim(_rprim);
   _rprim = strainMatrix*_rprim;
   _rprim = _rprim + buffRprim;
+  geometry::changeBasis(_rprim, _xcart, _xred, false);
 }
diff --git a/test/histcustommodes_phonons.h b/test/histcustommodes_phonons.h
@@ -65,14 +65,14 @@ class HistCMPhonons : public CxxTest::TestSuite
       histC.setRandomType(HistCustomModes::Normal);
       histC.setStatistics(HistCustomModes::Classical);
       histC.setInstableModes(HistCustomModes::Ignore);
-      
+
       histQ.setRandomType(HistCustomModes::Normal);
       histQ.setStatistics(HistCustomModes::Quantum);
       histQ.setInstableModes(HistCustomModes::Ignore);
-     
+
       histC.buildHist(qptGrid, temperature, strainBoundsIso, ntime);
       histQ.buildHist(qptGrid, temperature, strainBoundsIso, ntime);
-     
+
       Supercell first(histC,0);
       first.findReference(ref);
       DispDB::qptTree modes = { { {0,0,0}, { { 7, 1, 1 } } } };  // mode 7 with 1 amplitude 1 energye (we don't care here)
@@ -84,23 +84,96 @@ class HistCMPhonons : public CxxTest::TestSuite
         currentC.setReference(first);
         projC[itime] = currentC.projectOnModes(ref, db,modes,Supercell::NONE,false)[0];
       }
-      
+
       double meanC = utils::mean(projC.begin(),projC.end());
       double devC = utils::deviation(projC.begin(),projC.end(),meanC);
 
 
-     for ( unsigned itime = 0 ; itime < ntime ; ++itime ) {
+      for ( unsigned itime = 0 ; itime < ntime ; ++itime ) {
         Supercell currentQ(histQ,itime);
         currentQ.setReference(first);
         projQ[itime] = currentQ.projectOnModes(ref, db,modes,Supercell::NONE,false)[0];
       }
-      
+
       double meanQ = utils::mean(projQ.begin(),projQ.end());
       double devQ = utils::deviation(projQ.begin(),projQ.end(),meanQ);
 
       TS_ASSERT_DELTA(devQ,devC,1e-2)
 
-
-     }   
+
+    }   
+
+    void testCondensationAnalysisNoStrain() {
+#include "PTO_DDB.hxx"
+      Dtset ref;
+      ref.readFromFile("ref_PTO_DDB");
+      Ddb* ddb = Ddb::getDdb("ref_PTO_DDB");
+      DispDB db; 
+      db.computeFromDDB(*ddb);
+      delete ddb;
+      HistCustomModes histQ(ref,db);
+      const geometry::vec3d& qptGrid={1, 1, 1}; 
+      std::map<HistCustomModes::StrainDistBound,double> strainBoundsIso;
+
+      unsigned ntime = 100;
+      const double temperature = 100;
+
+      histQ.setRandomType(HistCustomModes::Normal);
+      histQ.setStatistics(HistCustomModes::Quantum);
+      histQ.setInstableModes(HistCustomModes::Ignore);
+
+      histQ.buildHist(qptGrid, temperature, strainBoundsIso, ntime);
+
+      Supercell first(histQ,0);
+      first.findReference(ref);
+      DispDB::qptTree modes = { { {0,0,0}, { { 7, 1, 1 } } } };  // mode 7 with 1 amplitude 1 energye (we don't care here)
+
+      for ( unsigned itime = 0 ; itime < ntime ; ++itime ) {
+        Supercell currentQ(histQ,itime);
+        currentQ.setReference(first);
+        double projQ = currentQ.projectOnModes(ref, db,modes,Supercell::NONE,false)[0];
+        double dispAmp = histQ.getDispAmplitudes(itime).begin()->second[1].amplitude; // Mode 7 is at position 1 (0 1 2 3 4 5 instables + ASR)
+        TS_ASSERT_DELTA(projQ,dispAmp,1e-5)
+      }
+    }   
+
+    void testCondensationAnalysisStrain() {
+#include "PTO_DDB.hxx"
+      Dtset ref;
+      ref.readFromFile("ref_PTO_DDB");
+      Ddb* ddb = Ddb::getDdb("ref_PTO_DDB");
+      DispDB db; 
+      db.computeFromDDB(*ddb);
+      delete ddb;
+      HistCustomModes histQ(ref,db);
+      const geometry::vec3d& qptGrid={1, 1, 1}; 
+      std::map<HistCustomModes::StrainDistBound,double> strainBoundsIso {
+        {HistCustomModes::IsoMin, 0.0001 }, {HistCustomModes::IsoMax, 0.1},
+        {HistCustomModes::TetraMin, 0.0001 }, {HistCustomModes::TetraMax, 0.1},
+        {HistCustomModes::ShearMin, 0.0001 }, {HistCustomModes::ShearMax, 0.1}
+      };
+
+      unsigned ntime = 100;
+      const double temperature = 100;
+
+      histQ.setRandomType(HistCustomModes::Normal);
+      histQ.setStatistics(HistCustomModes::Quantum);
+      histQ.setInstableModes(HistCustomModes::Ignore);
+
+      histQ.buildHist(qptGrid, temperature, strainBoundsIso, ntime);
+
+      Supercell first(histQ,0);
+      first.findReference(ref);
+      DispDB::qptTree modes = { { {0,0,0}, { { 7, 1, 1 } } } };  // mode 7 with 1 amplitude 1 energye (we don't care here)
+
+      for ( unsigned itime = 0 ; itime < ntime ; ++itime ) {
+        Supercell currentQ(histQ,itime);
+        currentQ.setReference(first);
+        auto disp = currentQ.getDisplacement(ref,false);
+        double projQ = currentQ.projectOnModes(ref, db,modes,Supercell::NONE,false)[0];
+        double dispAmp = histQ.getDispAmplitudes(itime).begin()->second[1].amplitude; // Mode 7 is at position 1 (0 1 2 3 4 5 instables + ASR)
+        TS_ASSERT_DELTA(projQ,dispAmp,1e-5)
+      }
+    }   
 
 };