diff --git a/include/hist/histcustommodes.hpp b/include/hist/histcustommodes.hpp index 630b122..ba07642 100644 --- 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 index 12d28bc..ba8c76d 100644 --- 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 index a47c168..791d190 100644 --- 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 index ba3b0e0..c389bb3 100644 --- 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 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 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) + } + } };