Merge pull request #133 from mir-group/phonopy_fix

Phonopy fix

src/delta_function.cpp

@@ -53,13 +53,14 @@ double GaussianDeltaFunction::getSmearing(const double &energy,
 }
 
 AdaptiveGaussianDeltaFunction::AdaptiveGaussianDeltaFunction(
-    BaseBandStructure &bandStructure) {
+    BaseBandStructure &bandStructure, double broadeningCutoff_) {
   auto tup = bandStructure.getPoints().getMesh();
   auto mesh = std::get<0>(tup);
   qTensor = bandStructure.getPoints().getCrystal().getReciprocalUnitCell();
   qTensor.row(0) /= mesh(0);
   qTensor.row(1) /= mesh(1);
   qTensor.row(2) /= mesh(2);
+  broadeningCutoff = broadeningCutoff_;
 }
 
 double
@@ -77,8 +78,14 @@ AdaptiveGaussianDeltaFunction::getSmearing(const double &energy,
   }
   sigma = prefactor * sqrt(sigma / 6.);
 
-  if (sigma == 0.)
+  if (sigma == 0.) {
     return 0.;
+  }
+
+  // if the smearing is smaller than 0.1 meV, we renormalize it
+  if (sigma < broadeningCutoff ) {
+    sigma = broadeningCutoff;
+  }
 
   if (abs(energy) > 2. * sigma)
     return 0.;

src/delta_function.h

@@ -103,7 +103,8 @@ class AdaptiveGaussianDeltaFunction : public DeltaFunction {
    * @param bandStructure: mainly to see what mesh of points and crystal is
    * being used, and to prepare a suitable scaling of velocities.
    */
-  explicit AdaptiveGaussianDeltaFunction(BaseBandStructure &bandStructure);
+  explicit AdaptiveGaussianDeltaFunction(BaseBandStructure &bandStructure,
+                                         double broadeningCutoff_=0.0001 / energyRyToEv);
 
   /** Method to obtain the value of smearing.
    * @param energy: the energy difference.
@@ -126,6 +127,8 @@ class AdaptiveGaussianDeltaFunction : public DeltaFunction {
 protected:
   int id = DeltaFunction::adaptiveGaussian;
 
+  double broadeningCutoff;
+
   const double prefactor = 1.; // could be exposed to the user
   Eigen::Matrix3d qTensor;     // to normalize by the number of wavevectors.
 

src/parser/ifc3_parser.cpp

@@ -389,6 +389,39 @@ Interaction3Ph IFC3Parser::parseFromPhono3py(Context &context,
   // set up so that the first nAtoms*dimSup of the superCell are unit cell
   // atom #1, and the second nAtoms*dimSup are atom #2, and so on.
 
+
+
+  // First, we open the phonopyDispFileName to check the distance units.
+  // Phono3py uses Bohr for QE calculations and Angstrom for VASP
+  // so, here we decide which is the correct unit,
+  // by parsing the phono3py_disp.yaml file
+  double distanceConversion = 1. / distanceBohrToAng;
+  {
+    std::string fileName = context.getPhonopyDispFileName();
+    std::ifstream infile;
+    if (fileName.empty()) {
+      Error("Phono3py parser required file phonopyDispFileName "
+            "(phonon3py_disp.yaml) not specified in input file.");
+    }
+    infile.open(fileName);
+    if (not infile.is_open()) {
+      Error("Phono3py parser couldn't open the file (phono3py_disp.yaml)");
+    }
+    if (mpi->mpiHead())
+      std::cout << "Reading in " + fileName << std::endl;
+
+    std::string line;
+    while (std::getline(infile, line)) {
+      if (line.find("length") != std::string::npos) {
+        if (line.find("au") != std::string::npos) {
+          distanceConversion = 1.; // distances already in Bohr
+          break;
+        }
+      }
+    }
+    infile.close();
+  }
+
   // First, read in the information form disp_fc3.yaml
   int numAtoms = crystal.getNumAtoms();
 
@@ -437,12 +470,12 @@ Interaction3Ph IFC3Parser::parseFromPhono3py(Context &context,
       std::string temp = line.substr(5, 62);// just the elements
       int idx1 = temp.find(',');
       lattice(iLattice, 0) =
-          std::stod(temp.substr(0, idx1)) / distanceBohrToAng;
+          std::stod(temp.substr(0, idx1)) * distanceConversion;
       int idx2 = temp.find(',', idx1 + 1);
       lattice(iLattice, 1) =
-          std::stod(temp.substr(idx1 + 1, idx2)) / distanceBohrToAng;
+          std::stod(temp.substr(idx1 + 1, idx2)) * distanceConversion;
       lattice(iLattice, 2) =
-          std::stod(temp.substr(idx2 + 1)) / distanceBohrToAng;
+          std::stod(temp.substr(idx2 + 1)) * distanceConversion;
       iLattice++;
     }
     if (line.find("lattice:") != std::string::npos) {

src/parser/phonopy_input_parser.cpp

@@ -44,6 +44,8 @@ std::tuple<Crystal, PhononH0> PhonopyParser::parsePhHarmonic(Context &context) {
   // If both superCells were the same, both are written to disp_fc3.yaml,
   // and disp_fc2.yaml will not have been created.
 
+  double distanceConversion = 1. / distanceBohrToAng;
+
   // open input file
   auto fileName = context.getDispFC2FileName();
   std::ifstream infile(fileName);
@@ -104,12 +106,12 @@ std::tuple<Crystal, PhononH0> PhonopyParser::parsePhHarmonic(Context &context) {
       std::string temp = line.substr(5, 62); // just the elements
       int idx1 = temp.find(',');
       supLattice(ilatt, 0) =
-          std::stod(temp.substr(0, idx1)) / distanceBohrToAng;
+          std::stod(temp.substr(0, idx1));
       int idx2 = temp.find(',', idx1 + 1);
       supLattice(ilatt, 1) =
-          std::stod(temp.substr(idx1 + 1, idx2)) / distanceBohrToAng;
+          std::stod(temp.substr(idx1 + 1, idx2));
       supLattice(ilatt, 2) =
-          std::stod(temp.substr(idx2 + 1)) / distanceBohrToAng;
+          std::stod(temp.substr(idx2 + 1));
       ilatt++;
     }
     if (line.find("lattice:") != std::string::npos) {
@@ -145,8 +147,7 @@ std::tuple<Crystal, PhononH0> PhonopyParser::parsePhHarmonic(Context &context) {
   // we have to do this first, because we need to use this info
   // to allocate the below data storage.
   Eigen::Vector3i qCoarseGrid;
-  while (infile) {
-    getline(infile, line);
+  while (std::getline(infile, line)) {
 
     // In the case where force constants where generated with different
     // superCells for fc2 and fc3, the label we need is dim_fc2.
@@ -231,12 +232,12 @@ std::tuple<Crystal, PhononH0> PhonopyParser::parsePhHarmonic(Context &context) {
       std::string temp = line.substr(9, 67); // just the elements
       int idx1 = temp.find(',');
       directUnitCell(ilatt, 0) =
-          std::stod(temp.substr(0, idx1)) / distanceBohrToAng;
+          std::stod(temp.substr(0, idx1));
       int idx2 = temp.find(',', idx1 + 1);
       directUnitCell(ilatt, 1) =
-          std::stod(temp.substr(idx1 + 1, idx2)) / distanceBohrToAng;
+          std::stod(temp.substr(idx1 + 1, idx2));
       directUnitCell(ilatt, 2) =
-          std::stod(temp.substr(idx2 + 1, temp.find(']'))) / distanceBohrToAng;
+          std::stod(temp.substr(idx2 + 1, temp.find(']')));
       ilatt++;
     }
     if (line.find("lattice:") != std::string::npos) {
@@ -246,9 +247,19 @@ std::tuple<Crystal, PhononH0> PhonopyParser::parsePhHarmonic(Context &context) {
     if (line.find("reciprocal_lattice:") != std::string::npos) {
       break;
     }
+
+    if (line.find("length") != std::string::npos) {
+      if (line.find("au") != std::string::npos) {
+        distanceConversion = 1.; // distances already in Bohr
+      }
+    }
   }
   infile.close();
 
+  // convert distances to Bohr
+  supLattice *= distanceConversion;
+  directUnitCell *= distanceConversion;
+
   // Process the information that has been read in
   // =================================================
   // calculate species mass values (in Ry) for later use

src/statistics_sweep.cpp

@@ -377,9 +377,11 @@ void StatisticsSweep::printInfo() {
   if (particle.isElectron()) {
     std::cout << "Fermi level: " << fermiLevel * energyRyToEv << " (eV)"
               << std::endl;
+    std::cout << "Index, temperature, chemical potential, doping concentration\n";
+  } else {
+    std::cout << "Index, temperature\n";
   }
 
-  std::cout << "Index, temperature, chemical potential, doping concentration\n";
 
   for (int iCalc = 0; iCalc < numCalculations; iCalc++) {
     double temp = infoCalculations(iCalc, 0);