Patch error with qpoints.yaml files from phonopy

TODO: Fix equivalent atoms in BORN file

Author: Jordan Bieder

src/canvas/canvasphonons.cpp

@@ -161,9 +161,9 @@ bool CanvasPhonons::readDdb(const std::string& filename) {
      * we can assume this is true. The user should know what he does.
      * We cannot do in an other way since qpoints.yaml does not always provide
      * the structure (1.10 at least does not)
-     */
     if ( !(*_ddb==_reference) ) 
       throw EXCEPTION("It seems reference and DDB are different structures",ERRDIV);
+     */
     _ddb->buildFrom(_reference); 
     DispDB disp;
     try {

src/io/ddb.cpp

@@ -199,10 +199,8 @@ geometry::mat3d Ddb::getZeff(const unsigned iatom) const {
     throw EXCEPTION("Atom "+utils::to_string(iatom)+" is not in DDB", ERRDIV);
 
   auto data = this->getDdb(qpt);
-  mat3d zeff;
-  mat3d count;
-  for ( auto &e : count ) e = 0e0;
-  for ( auto &e : zeff ) e = 0e0;
+  mat3d zeff = {0};
+  mat3d count = {0};
   const double twopi = 2*phys::pi;
 
 	/* Read values from ddb into _zeff*/
@@ -383,6 +381,20 @@ void Ddb::setZeff(const unsigned iatom, const geometry::mat3d &zeff) {
   mat3d rprimTranspose = transpose(_rprim);
   geometry::mat3d d2red = (rprimTranspose * (zeff  * _gprim))*twopi;
 
+  for ( int i = 0 ; i < 2 ; ++i ) {
+    for ( unsigned idir1 = 0 ; idir1 < 3 ; ++idir1 ) {
+      for ( unsigned idir2 = 0 ; idir2 < 3 ; ++idir2 ) {
+        auto alreadyIn = std::find_if(block.begin(),block.end(),[=](
+              const std::pair<std::array<unsigned,4>,std::complex<double>>& entry) { 
+            auto& array = entry.first;
+            return ( array[0] == idir1 && array[1] ==(_natom+1) && array[2] == idir2 && array[3]==iatom) || 
+            ( array[0] == idir2 && array[1] ==iatom && array[2] == idir1 && array[3]==(_natom+1)); 
+            });
+        if ( alreadyIn != block.end() ) block.erase(alreadyIn);
+      }
+    }
+  }
+
   for ( unsigned idir1 = 0 ; idir1 < 3 ; ++idir1 )
     for ( unsigned idir2 = 0 ; idir2 < 3 ; ++idir2 ) {
       block.push_back(

src/io/ddbphonopy.cpp

@@ -118,7 +118,7 @@ void DdbPhonopy::readFromFile(const std::string& filename) {
     std::string bornfile = utils::dirname(filename)+"/BORN";
     std::ifstream born(bornfile);
     if ( born ) {
-      _zion.resize(_typat.size()); 
+      _zion.resize(_typat.size(),0); 
       unsigned iline = 0;
       std::string line;
       utils::getline(born,line,iline); // First line is an conversion factor eventually.
@@ -161,15 +161,16 @@ void DdbPhonopy::buildFrom(const Dtset& dtset) {
       break;
     }
   }
-  if ( currentZ.size() != _natom ) {
-    throw EXCEPTION("Number of Zeff is smaller than natom is not supported yet",ERRDIV);
+  if ( currentZ.size() > 0 && currentZ.size() != _natom ) {
+    auto e = EXCEPTION("Number of Zeff is smaller than natom is not supported yet",ERRWAR);
+    std::clog << e.fullWhat() << std::endl;
 #ifdef HAVE_SPGLIB
     SpglibDataset* spgDtset = this->getSpgDtset(0.001*phys::A2b);
     if ( spgDtset != nullptr ) {
-      unsigned prevAtom = 0;
+      const geometry::mat3d zero = {0};
       for ( unsigned iatom = 0, izeff=0; iatom < _natom; ++iatom ) {
         unsigned eqatom = static_cast<unsigned>(spgDtset->equivalent_atoms[iatom]);
-        this->setZeff(iatom,eqatom == prevAtom ? currentZ[izeff] : currentZ[++izeff]);
+        this->setZeff(iatom,(iatom == eqatom ? currentZ[izeff++]:zero));
       }
     }
     spg_free_dataset(spgDtset);

src/io/procar.cpp

@@ -70,7 +70,7 @@ void Procar::readFromFile(const std::string& filename) {
     std::getline(procar,tmp,':');
     procar >> natom;
 
-    if ( procar.eof() ) break;
+    if ( procar && procar.eof() ) break;
     if ( !procar )
       throw EXCEPTION("Cannot read properly the parameters",ERRDIV);