diff --git a/include/phonons/phononmode.hpp b/include/phonons/phononmode.hpp index f52b7ab..a5cd052 100644 --- a/include/phonons/phononmode.hpp +++ b/include/phonons/phononmode.hpp @@ -128,13 +128,12 @@ class PhononMode { /** * Calculate linear response from ddb - * @param _qpt The Q - Point -> always GAMMA * @param Edir direction vector for electric field in cartesian coordinates * @param Eamp Amplitude of electric field * @param ddb complete ddb */ - std::vector lin_res(const geometry::vec3d& _qpt, geometry::vec3d &Edir, double Eamp, const Ddb& ddb); + std::vector lin_res(const geometry::vec3d &Edir, double Eamp, const Ddb& ddb); }; diff --git a/src/phonons/dispdb.cpp b/src/phonons/dispdb.cpp index 6a22d11..5eb8a71 100644 --- a/src/phonons/dispdb.cpp +++ b/src/phonons/dispdb.cpp @@ -324,8 +324,7 @@ void DispDB::loadFromEigParserPhonon(EigParserPhonons& eigparser) { void DispDB::linearResponseE(std::vector &Edir, double A, Ddb &ddb) { PhononMode respE; geometry::vec3d E_dir = {{ Edir[0], Edir[1], Edir[2] }}; - geometry::vec3d gamma = {{ 0,0,0 }}; - auto disp_E = respE.lin_res(gamma, E_dir, A, ddb); + auto disp_E = respE.lin_res(E_dir, A, ddb); _linResE.resize(3*ddb.natom()); for ( unsigned i = 0 ; i < 3*ddb.natom() ; ++i ) { _linResE[i].real(disp_E[i]); diff --git a/src/phonons/phononmode.cpp b/src/phonons/phononmode.cpp index 6fd126e..60eca22 100644 --- a/src/phonons/phononmode.cpp +++ b/src/phonons/phononmode.cpp @@ -133,10 +133,9 @@ void PhononMode::computeForceCst(const geometry::vec3d& qpt, const Ddb& ddb) { /*Calculate Linear Repsonse of Phonons to a static dielectric field from DFPT*/ -std::vector PhononMode::lin_res(const geometry::vec3d& _qpt, geometry::vec3d &E_vec, double E_Amp, const Ddb& ddb) { +std::vector PhononMode::lin_res(const geometry::vec3d &E_vec, double E_Amp, const Ddb& ddb) { #ifndef HAVE_EIGEN throw EXCEPTION("To use this functionnality you need to compile with EIGEN support",ERRABT); - (void) _qpt; (void) E_vec; (void) E_Amp; (void) ddb; @@ -153,16 +152,7 @@ std::vector PhononMode::lin_res(const geometry::vec3d& _qpt, geometry::v for ( unsigned iatom = 0 ; iatom < _natom ; ++iatom ) _zeff[iatom] = ddb.getZeff(iatom); /// get BEC-Tensors - for (unsigned i = 0; i < _zeff.size(); i++) { - if (_zeff[i][geometry::mat3dind( 2,1 )] == 0 ){ - _zeff[i][geometry::mat3dind( 1, 2)] = 0; - } - if (_zeff[i][geometry::mat3dind( 3,1 )] == 0 ){ - _zeff[i][geometry::mat3dind( 1, 3)] = 0; - } - } - - geometry::mat3d rprim = ddb.rprim(); /// get rprim + const geometry::mat3d rprim = ddb.rprim(); /// get rprim _rprim << rprim[0], rprim[1],rprim[2], rprim[3], rprim[4], rprim[5], rprim[6], rprim[7], rprim[8]; @@ -182,12 +172,13 @@ std::vector PhononMode::lin_res(const geometry::vec3d& _qpt, geometry::v if (freq_gamma[m] < -1) throw EXCEPTION("NEGATIVE PHONON FREQUENCY FOUND: The linear response to an Electric-Field calculation makes only sense in stable structures. Fully relax your structure",ERRDIV); } - for ( unsigned iatom = 0 ; iatom < _natom ; ++iatom ) { - _mass[iatom] = MendeTable.mass[ddb.znucl().at(ddb.typat().at(iatom)-1)]; // type starts at 1 - } - for ( unsigned i = 0; i < disp_gamma.size(); ++i ) { - disp_gamma[i] = disp_gamma[i].real()*pow(phys::amu_emass,0.5); - } + _mass.resize(_natom); + for ( unsigned iatom = 0 ; iatom < _natom ; ++iatom ) { + _mass[iatom] = MendeTable.mass[ddb.znucl().at(ddb.typat().at(iatom)-1)]; // type starts at 1 + } + for ( unsigned i = 0; i < disp_gamma.size(); ++i ) { + disp_gamma[i] = disp_gamma[i]*sqrt(phys::amu_emass); + } /*--------- Calculate diplacement under electric filed ---------*/ /* 1. Calculate polarity of each mode and store*/