qmmm.cc

Go to the documentation of this file.

/*
 *            Copyright 2009-2020 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
 
// Standard includes
#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <filesystem>
#include <sstream>
 
// Third party includes
#include <numeric>
#include <stdexcept>
 
// Local VOTCA includes
#include "votca/tools/property.h"
#include "votca/tools/tokenizer.h"
#include "votca/xtp/IndexParser.h"
#include "votca/xtp/jobtopology.h"
#include "votca/xtp/qmregion.h"
 
// Local private VOTCA includes
#include "qmmm.h"
#include "votca/xtp/qmstate.h"
 
namespace votca {
namespace xtp {
 
void QMMM::ParseSpecificOptions(const tools::Property& options) {
 
  print_regions_pdb_ = options.get(".print_regions_pdb").as<bool>();
  max_iterations_ = options.get(".max_iterations").as<Index>();
  regions_def_.second = options.get(".regions");
  regions_def_.first = mapfile_;
  use_gs_for_ex_ = options.get("io_jobfile.use_gs_for_ex").as<bool>();
 
  states_ = options.get("io_jobfile.states").as<std::vector<QMState>>();
  which_segments_ = options.get("io_jobfile.segments").as<std::string>();
 
  bool groundstate_found = std::any_of(
      states_.begin(), states_.end(),
      [](const QMState& state) { return state.Type() == QMStateType::Gstate; });
  // The n job is the reference -j read subtracts. A job with
  // environment_screening converges that ground state itself and reports
  // its site energies relative to it, and refuses state n anyway.
  if (!groundstate_found && !hasEnvironmentScreening()) {
    states_.push_back(QMState("n"));
  }
}
 
Job::JobResult QMMM::EvalJob(const Topology& top, Job& job, QMThread& Thread) {
  std::chrono::time_point<std::chrono::system_clock> start =
      std::chrono::system_clock::now();
 
  std::string qmmm_work_dir = "QMMM";
  if (!this->hasQMRegion()) {
    qmmm_work_dir = "MMMM";
  }
  std::string frame_dir =
      "frame_" + boost::lexical_cast<std::string>(top.getStep());
  std::string job_dir =
      "job_" + std::to_string(job.getId()) + "_" + job.getTag();
  std::filesystem::path arg_path;
  std::string workdir =
      (arg_path / qmmm_work_dir / frame_dir / job_dir).generic_string();
  std::filesystem::create_directories(workdir);
  Job::JobResult jres = Job::JobResult();
  Logger& pLog = Thread.getLogger();
 
  JobTopology jobtop = JobTopology(job, pLog, workdir);
  if (job.getInput().exists("restart")) {
    std::string checkptfile = job.getInput().get("restart").as<std::string>();
    XTP_LOG(Log::error, pLog)
        << TimeStamp() << " Restart job from " << checkptfile << std::flush;
    jobtop.ReadFromHdf5(checkptfile);
  } else {
    jobtop.BuildRegions(top, regions_def_);
  }
 
  if (print_regions_pdb_) {
    std::string pdb_filename = "regions.pdb";
    XTP_LOG(Log::error, pLog) << TimeStamp() << " Writing jobtopology to "
                              << (workdir + "/" + pdb_filename) << std::flush;
    jobtop.WriteToPdb(workdir + "/" + pdb_filename);
  }
 
  // Embedded GW-BSE: its reaction field depends only on geometry and the
  // auxiliary basis, so a bad environment is caught before the loop.
  if (QMRegion* qmregion = dynamic_cast<QMRegion*>(jobtop.Regions()[0].get());
      qmregion != nullptr && qmregion->EnvironmentScreeningEnabled()) {
    try {
      qmregion->CheckScreeningEnvironment(jobtop.Regions());
    } catch (std::exception& e) {
      XTP_LOG(Log::error, pLog)
          << TimeStamp() << " Screening environment check failed: " << e.what()
          << std::flush;
      jres.setStatus(Job::JobStatus::FAILED);
      jres.setError(std::string("Screening environment check failed: ") +
                    e.what());
      return jres;
    }
  }
 
  Index no_static_regions = std::accumulate(
      jobtop.begin(), jobtop.end(), 0, [](Index count, const auto& region) {
        return count += Index(region->Converged());
      });
 
  bool no_top_scf = false;
  if (jobtop.size() - no_static_regions < 2) {
    XTP_LOG(Log::error, pLog)
        << TimeStamp() << " Only " << jobtop.size() - no_static_regions
        << " scf region is used. The remaining regions are static. So no "
           "inter regions scf is required. "
        << std::flush;
    no_top_scf = true;
    max_iterations_ = 1;
  }
  Index iteration = 0;
  for (; iteration < max_iterations_; iteration++) {
 
    XTP_LOG(Log::error, pLog)
        << TimeStamp() << " --Inter Region SCF Iteration " << iteration + 1
        << " of " << max_iterations_ << std::flush;
 
    for (std::unique_ptr<Region>& region : jobtop) {
      XTP_LOG(Log::error, pLog)
          << TimeStamp() << " Evaluating " << region->identify() << " "
          << region->getId() << std::flush;
      region->Reset();
      region->Evaluate(jobtop.Regions());
      if (!region->Successful()) {
        jres.setStatus(Job::JobStatus::FAILED);
        jres.setError(region->ErrorMsg());
        return jres;
      }
    }
 
    std::string checkpointfilename =
        "checkpoint_iter_" + std::to_string(iteration + 1) + ".hdf5";
    XTP_LOG(Log::error, pLog) << TimeStamp() << " Writing checkpoint to "
                              << checkpointfilename << std::flush;
    try {
      jobtop.WriteToHdf5(workdir + "/" + checkpointfilename);
    } catch (std::exception& e) {
      XTP_LOG(Log::error, pLog)
          << TimeStamp()
          << " Could not write checkpoint file due to exception: " << e.what()
          << std::flush;
    }
 
    if (!no_top_scf) {
      std::vector<bool> converged_regions;
      for (std::unique_ptr<Region>& region : jobtop) {
        converged_regions.push_back(region->Converged());
      }
 
      double etot = std::accumulate(
          jobtop.begin(), jobtop.end(), 0.0,
          [](double e, const auto& region) { return region->Etotal() + e; });
 
      XTP_LOG(Log::error, pLog) << TimeStamp() << " --Total Energy all regions "
                                << etot << std::flush;
 
      bool all_regions_converged =
          std::all_of(converged_regions.begin(), converged_regions.end(),
                      [](bool i) { return i; });
 
      if (all_regions_converged) {
        XTP_LOG(Log::error, pLog)
            << TimeStamp() << " Job converged after " << iteration + 1
            << " iterations." << std::flush;
        jres.setStatus(Job::JobStatus::COMPLETE);
        break;
      }
      if (iteration == max_iterations_ - 1) {
        XTP_LOG(Log::error, pLog)
            << TimeStamp() << " Job did not converge after " << iteration + 1
            << " iterations.\n Writing results to jobfile." << std::flush;
        jres.setStatus(Job::JobStatus::FAILED);
        jres.setError("Inter Region SCF did not converge in " +
                      std::to_string(max_iterations_) + " iterations.");
      }
    } else {
      jres.setStatus(Job::JobStatus::COMPLETE);
    }
  }
 
  // Embedded GW-BSE: the loop above converged the ground state; the
  // excitation is screened by the polar regions in one GW-BSE run.
  QMRegion* screened_qmregion = nullptr;
  double e_ground = 0.0;
  if (jres.getStatus() == Job::JobStatus::COMPLETE) {
    QMRegion* qmregion = dynamic_cast<QMRegion*>(jobtop.Regions()[0].get());
    if (qmregion != nullptr && qmregion->EnvironmentScreeningEnabled()) {
      // The converged ground state's total energy, before the screened
      // run replaces the QM region's energy with that of the state.
      for (const std::unique_ptr<Region>& reg : jobtop) {
        e_ground += reg->Etotal();
      }
      screened_qmregion = qmregion;
      try {
        qmregion->EvaluateScreenedGWBSE(jobtop.Regions());
      } catch (std::exception& e) {
        XTP_LOG(Log::error, pLog)
            << TimeStamp() << " Screened GW-BSE failed: " << e.what()
            << std::flush;
        jres.setStatus(Job::JobStatus::FAILED);
        jres.setError(std::string("Screened GW-BSE failed: ") + e.what());
        return jres;
      }
      try {
        jobtop.WriteToHdf5(workdir + "/checkpoint_screened.hdf5");
      } catch (std::exception& e) {
        XTP_LOG(Log::error, pLog)
            << TimeStamp()
            << " Could not write checkpoint file due to exception: " << e.what()
            << std::flush;
      }
    }
  }
 
  tools::Property results;
  tools::Property& jobresult = results.add("output", "");
  tools::Property& regionsresults = jobresult.add("regions", "");
  double etot = 0.0;
  double charge = 0.0;
  for (const std::unique_ptr<Region>& reg : jobtop) {
    reg->AddResults(regionsresults);
    etot += reg->Etotal();
    charge += reg->charge();
  }
  std::chrono::time_point<std::chrono::system_clock> end =
      std::chrono::system_clock::now();
  std::chrono::duration<double> elapsed_time = end - start;
  jobresult.add("E_tot", std::to_string(etot * tools::conv::hrt2ev));
  jobresult.add("Compute_Time", std::to_string(Index(elapsed_time.count())));
  jobresult.add("Total_Charge", std::to_string(charge));
  if (!no_top_scf) {
    jobresult.add("Iterations", std::to_string(iteration + 1));
  }
  if (screened_qmregion != nullptr) {
    jobresult.add("E_ground", std::to_string(e_ground * tools::conv::hrt2ev));
    screened_qmregion->AppendScreenedSiteEnergies(jobresult);
  }
  jres.setOutput(results);
  return jres;
}
 
bool QMMM::hasEnvironmentScreening() const {
  for (const tools::Property& reg : regions_def_.second) {
    if (reg.name() == "qmregion" && reg.exists("environment_screening")) {
      const tools::Property& env = reg.get("environment_screening");
      return env.ifExistsReturnElseReturnDefault<bool>("enabled", true);
    }
  }
  return false;
}
 
bool QMMM::hasQMRegion() const {
  Logger log;
  QMRegion QMdummy(0, log, "");
  return std::any_of(regions_def_.second.begin(), regions_def_.second.end(),
                     [&](const tools::Property& reg) {
                       return reg.name() == QMdummy.identify();
                     });
}
 
std::string QMMM::getFirstRegionName() const {
  for (const auto& reg : regions_def_.second) {
    if (reg.get("id").as<Index>() == 0) {
      return reg.name();
    }
  }
  throw std::runtime_error("region ids do not start at 0");
  return "";
}
 
void QMMM::WriteJobFile(const Topology& top) {
 
  std::cout << std::endl
            << "... ... Writing job file " << jobfile_ << std::flush;
 
  std::ofstream ofs;
  ofs.open(jobfile_, std::ofstream::out);
  if (!ofs.is_open()) {
    throw std::runtime_error("\nERROR: bad file handle: " + jobfile_);
  }
 
  std::vector<Index> segments_to_write;
  if (which_segments_ == "all") {
    for (Index i = 0; i < Index(top.Segments().size()); ++i) {
      segments_to_write.push_back(i);
    }
  } else {
    segments_to_write = IndexParser().CreateIndexVector(which_segments_);
  }
 
  ofs << "<jobs>" << std::endl;
  Index jobid = 0;
  for (Index segID : segments_to_write) {
    const Segment& seg = top.Segments()[segID];
    for (const QMState& state : states_) {
      Job job = createJob(seg, state, jobid);
      job.ToStream(ofs);
      jobid++;
    }
  }
 
  ofs << "</jobs>" << std::endl;
  ofs.close();
  std::cout << std::endl
            << "... ... In total " << jobid << " jobs" << std::flush;
  return;
}
 
Job QMMM::createJob(const Segment& seg, const QMState& state,
                    Index jobid) const {
  std::string marker = std::to_string(seg.getId()) + ":" + state.ToString();
  std::string tag = seg.getType() + "_" + marker;
 
  tools::Property Input;
  tools::Property& pInput = Input.add("input", "");
  pInput.add("site_energies", marker);
  tools::Property& regions = pInput.add("regions", "");
  tools::Property& region = regions.add(getFirstRegionName(), "");
  region.add("id", "0");
  if (hasQMRegion()) {
    region.add("state", state.ToString());
  }
  // Quasiparticle levels are levels of the neutral molecule: its geometry,
  // there being no qmcoords_pqp to look up.
  if ((use_gs_for_ex_ && (state.Type() == QMStateType::Singlet ||
                          state.Type() == QMStateType::Triplet)) ||
      state.Type().isGWState()) {
    region.add("segments", std::to_string(seg.getId()) + ":n");
  } else {
    region.add("segments", marker);
  }
  Job job(jobid, tag, Input, Job::AVAILABLE);
  return job;
}
 
void QMMM::ReadJobFile(Topology& top) {
 
  Index incomplete_jobs = 0;
 
  Eigen::Matrix<double, Eigen::Dynamic, 5> energies =
      Eigen::Matrix<double, Eigen::Dynamic, 5>::Zero(top.Segments().size(), 5);
  Eigen::Matrix<bool, Eigen::Dynamic, 5> found =
      Eigen::Matrix<bool, Eigen::Dynamic, 5>::Zero(top.Segments().size(), 5);
  // Site energies reported directly by screened jobs (relative to their
  // own ground state, Hartree), for e, h, s, t.
  Eigen::Matrix<double, Eigen::Dynamic, 4> screened =
      Eigen::Matrix<double, Eigen::Dynamic, 4>::Zero(top.Segments().size(), 4);
  Eigen::Matrix<bool, Eigen::Dynamic, 4> screened_found =
      Eigen::Matrix<bool, Eigen::Dynamic, 4>::Zero(top.Segments().size(), 4);
 
  tools::Property xml;
  xml.LoadFromXML(jobfile_);
  for (tools::Property* job : xml.Select("jobs.job")) {
 
    Index jobid = job->get("id").as<Index>();
    if (!job->exists("status")) {
      throw std::runtime_error(
          "Jobfile is malformed. <status> tag missing for job " +
          std::to_string(jobid));
    }
    if (job->get("status").as<std::string>() != "COMPLETE" ||
        !job->exists("output")) {
      incomplete_jobs++;
      continue;
    }
 
    std::vector<std::string> split =
        tools::Tokenizer(job->get("input.site_energies").as<std::string>(), ":")
            .ToVector();
 
    Index segid = std::stoi(split[0]);
    if (segid < 0 || segid >= Index(top.Segments().size())) {
      throw std::runtime_error("JobSegment id" + std::to_string(segid) +
                               " is not in topology for job " +
                               std::to_string(jobid));
    }
    QMState state;
    try {
      state.FromString(split[1]);
    } catch (std::runtime_error& e) {
      std::stringstream message;
      message << e.what() << " for job " << jobid;
      throw std::runtime_error(message.str());
    }
    if (job->exists("output.screened_site_energies")) {
      const tools::Property& se = job->get("output.screened_site_energies");
      const std::array<std::pair<const char*, QMStateType::statetype>, 4>
          types = {{{"e", QMStateType::Electron},
                    {"h", QMStateType::Hole},
                    {"s", QMStateType::Singlet},
                    {"t", QMStateType::Triplet}}};
      for (const auto& [tag, type] : types) {
        if (!se.exists(tag)) {
          continue;
        }
        const double value = se.get(tag).as<double>() * tools::conv::ev2hrt;
        if (screened_found(segid, type)) {
          // e and h come with every screened job of a segment (an s1 and a
          // pqp job, say), from the same ground state: they must agree.
          const double diff =
              std::abs(value - screened(segid, type)) * tools::conv::hrt2ev;
          if (diff > 1e-3) {
            throw std::runtime_error(
                "Screened jobs for segment " + std::to_string(segid) +
                " disagree on its " + std::string(tag) + " site energy by " +
                std::to_string(diff) +
                " eV; were they run with different settings?");
          }
          continue;
        }
        screened(segid, type) = value;
        screened_found(segid, type) = true;
      }
      // Its E_tot is not an n/e/h/s/t total energy to be differenced:
      // the site energies above are already relative to the ground state.
      continue;
    }
 
    const Index col = state.Type().Type();
    if (col > QMStateType::Gstate) {
      throw std::runtime_error(
          "Job " + std::to_string(jobid) + " has state " + state.ToString() +
          ", which -j read can only use from a job with environment_screening "
          "(it reports its own site energies).");
    }
    double energy = job->get("output.E_tot").as<double>() * tools::conv::ev2hrt;
    if (found(segid, col) != 0) {
      throw std::runtime_error("There are two entries in jobfile for segment " +
                               std::to_string(segid) +
                               " state:" + state.ToString());
    }
 
    energies(segid, col) = energy;
    found(segid, col) = true;
  }
 
  Eigen::Matrix<Index, 1, 5> found_states = found.colwise().count();
  std::cout << std::endl;
  for (Index i = 0; i < 5; i++) {
    if (found_states(i) > 0) {
      QMStateType type(static_cast<QMStateType::statetype>(i));
      std::cout << "Found " << found_states(i) << " states of type "
                << type.ToString() << std::endl;
    }
  }
  if (incomplete_jobs > 0) {
    std::cout << incomplete_jobs << " incomplete jobs found." << std::endl;
  }
 
  Eigen::Matrix<Index, 1, 4> screened_states = screened_found.colwise().count();
  for (Index i = 0; i < 4; i++) {
    if (screened_states(i) > 0) {
      QMStateType type(static_cast<QMStateType::statetype>(i));
      std::cout << "Found " << screened_states(i) << " screened site energies "
                << "of type " << type.ToString() << std::endl;
    }
  }
 
  Index overridden = 0;
  for (Segment& seg : top.Segments()) {
    Index segid = seg.getId();
    for (Index i = 0; i < 4; i++) {
      QMStateType type(static_cast<QMStateType::statetype>(i));
      if (screened_found(segid, i)) {
        // A screened job's own value is the one to use: it is relative to
        // the ground state of the very same run.
        if (found(segid, i) && found(segid, 4)) {
          ++overridden;
        }
        seg.setEMpoles(type, screened(segid, i));
      } else if (found(segid, i) && found(segid, 4)) {
        double energy = energies(segid, i) - energies(segid, 4);
        seg.setEMpoles(type, energy);
      }
    }
  }
  if (overridden > 0) {
    std::cout << overridden
              << " site energies were available both from n-referenced and "
                 "from screened jobs; the screened values were used."
              << std::endl;
  }
}
 
}  // namespace xtp
}  // namespace votca