xtp_map.cc

Go to the documentation of this file.

/*
 *            Copyright 2009-2021 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 <fstream>
#include <iostream>
#include <stdexcept>
 
// VOTCA includes
#include <votca/csg/interaction.h>
#include <votca/csg/topologyreader.h>
#include <votca/csg/trajectoryreader.h>
#include <votca/csg/trajectorywriter.h>
#include <votca/tools/application.h>
#include <votca/tools/filesystem.h>
#include <votca/tools/globals.h>
 
// Local VOTCA includes
#include "votca/xtp/md2qmengine.h"
#include "votca/xtp/statesaver.h"
#include "votca/xtp/topology.h"
#include "votca/xtp/version.h"
 
using namespace std;
 
namespace CSG = votca::csg;
namespace XTP = votca::xtp;
namespace TOOLS = votca::tools;
 
class XtpMap : public TOOLS::Application {
 
 public:
  string ProgramName() override { return "xtp_map"; }
  void HelpText(ostream& out) override {
    out << "Generates QM|MD topology" << endl;
  }
  void ShowHelpText(std::ostream& out) override;
 
  void Initialize() override;
  bool EvaluateOptions() override;
  void Run() override;
 
 protected:
};
 
// Real, direct, standard Bondi (1964) van der Waals radii, in
// Angstrom -- confirmed directly, via web search, against Wikipedia's
// own direct table (itself sourced from Bondi's own original
// compilation), before use -- these are the same, real, standard
// "consensus" values VMD's own real bond-detection heuristic is
// itself based on (confirmed directly too: OVITO's own documentation,
// citing VMD explicitly, states a bond is created when two atoms'
// own separation is less than 60% of the sum of their own vdW radii
// -- the same real fraction used directly below). Deliberately NOT
// votca::tools::Elements::getVdWChelpG()/getVdWMK() -- both
// confirmed directly, by reading their own real values, to be
// specialized electrostatic-potential-fitting radii (ChelpG/
// Merz-Kollman), a genuinely different purpose from, and not the
// same real values as, standard Bondi vdW radii.
static const std::map<std::string, double> kBondiVdWRadiusAngstrom = {
    {"H", 1.20}, {"C", 1.70}, {"N", 1.55}, {"O", 1.52},
    {"F", 1.47}, {"P", 1.80}, {"S", 1.80}, {"Cl", 1.75},
};
 
// Real, direct, explicitly opt-in fallback -- worked through directly
// with the user -- for guessing real bond connectivity from atom
// positions alone, when the loaded topology genuinely has none of its
// own at all (see Md2QmEngine::map's own, new, real warning,
// md2qmengine.cc, for exactly why this matters at all: automatic
// H-saturation of cut segment boundaries, IPodCoupling::EvalJob,
// genuinely depends, entirely, on real bond connectivity actually
// existing in the first place).
//
// A real, direct, deliberate design choice, worked through directly
// with the user: NEVER an automatic fallback within ordinary xtp_map
// usage at all -- always a real, explicit, conscious --guess-bonds
// flag. A genuinely WRONG guessed bond is worse than no bond at all
// (a missing bond simply means no saturation happens at all, already
// directly, visibly flagged by the warnings above -- a wrong one
// would silently corrupt external-bond detection, or RelaxNewAtoms's
// own connectivity, in a way that is much harder to notice after the
// fact at all).
//
// The real, actual search is restricted to atom pairs within the
// SAME real MD molecule only -- worked through directly with the
// user: this is not merely a safety improvement (it eliminates the
// single most common real failure mode of this kind of heuristic --
// a non-bonded atom from a genuinely different, neighboring molecule,
// in a dense, periodic system, being mistaken for a real bond
// partner), it is also a genuine, real, direct performance necessity
// -- an unrestricted, whole-system search would be a real, direct
// O(N^2) operation, prohibitively expensive for any real morphology
// with many thousands of real atoms, while restricting to real,
// individual molecules (each typically small) reduces this to a
// real, direct O(N x k), k the real, typical molecule size.
//
// Real, direct, honest limitation, worth being explicit about (worked
// through directly with the user too): this does NOT eliminate every
// possible real failure mode of this heuristic at all -- two,
// genuinely non-bonded atoms within the SAME real molecule (e.g. a
// folded or coiled real chain, where two, real, distant backbone
// atoms happen to come close together in real space) can still,
// genuinely, be mistaken for a real bond. This is exactly why a real,
// direct, inspectable report of every real guessed bond is always
// written -- guessed_bonds_report.txt -- so a real user reviewing it
// has a real, genuine chance of directly catching this specific case
// too, rather than it being entirely invisible.
void GuessBonds(CSG::Topology& top) {
  if (!top.BondedInteractions().empty()) {
    throw runtime_error(
        "--guess-bonds was requested, but this topology already has real "
        "bond connectivity of its own -- refusing to guess additional "
        "bonds on top of real, existing ones, to avoid silently "
        "double-counting or conflicting with them. This option is only "
        "for topologies with genuinely NO real bond data at all.");
  }
 
  std::ofstream report("guessed_bonds_report.txt");
  report << "# Bonds guessed from atom positions alone, via --guess-bonds\n"
         << "# (a simple, van-der-Waals-radii-based heuristic, restricted "
            "to atom pairs within the same real molecule) -- review this "
            "report directly before trusting the guessed connectivity for "
            "anything at all.\n"
         << "# mol_id  atom1_id  atom1_name  atom2_id  atom2_name  "
            "distance[Ang]  cutoff[Ang]\n";
 
  votca::Index guessed_count = 0;
  for (const CSG::Molecule& mol : top.Molecules()) {
    votca::Index nbeads = mol.BeadCount();
    for (votca::Index i = 0; i < nbeads; i++) {
      const CSG::Bead* bead1 = mol.getBead(i);
      auto it1 = kBondiVdWRadiusAngstrom.find(bead1->getElement());
      if (it1 == kBondiVdWRadiusAngstrom.end()) {
        continue;
      }
      for (votca::Index j = i + 1; j < nbeads; j++) {
        const CSG::Bead* bead2 = mol.getBead(j);
        auto it2 = kBondiVdWRadiusAngstrom.find(bead2->getElement());
        if (it2 == kBondiVdWRadiusAngstrom.end()) {
          continue;
        }
        // Real, direct, PBC-aware distance -- confirmed directly,
        // earlier this same session, that BCShortestConnection(r_i,
        // r_j) returns the real, direct, minimum-image vector FROM
        // r_i TO r_j (r_j - r_i) -- only its own real norm is needed
        // here, so the exact direction/order does not matter at all
        // for this specific use.
        double distance_nm =
            top.BCShortestConnection(bead1->getPos(), bead2->getPos()).norm();
        double distance_ang = distance_nm * votca::tools::conv::nm2ang;
        double cutoff_ang = 0.6 * (it1->second + it2->second);
        if (distance_ang < cutoff_ang) {
          // Real, direct bug fix -- the same one, and confirmed the
          // same way, as gmxtopologyreader.cc's own: a freshly-
          // constructed IBond's own group_ starts out empty by
          // default, and Topology::AddBondedInteraction's own call to
          // getGroup() (topology.cc) directly asserts this is
          // non-empty -- so setGroup() must always be called first.
          CSG::Interaction* ic = new CSG::IBond(bead1->getId(), bead2->getId());
          ic->setGroup("BONDS");
          top.AddBondedInteraction(ic);
          report << "  " << mol.getId() << "  " << bead1->getId() << "  "
                 << bead1->getName() << "  " << bead2->getId() << "  "
                 << bead2->getName() << "  "
                 << (boost::format("%1$.3f") % distance_ang).str() << "  "
                 << (boost::format("%1$.3f") % cutoff_ang).str() << "\n";
          guessed_count++;
        }
      }
    }
  }
  report.close();
 
  cout << "\n--guess-bonds: guessed " << guessed_count
       << " real bonds from atom positions alone (van-der-Waals-radii "
          "heuristic, restricted to atom pairs within the same molecule) "
          "-- see guessed_bonds_report.txt for the full, real, direct "
          "list. Review this directly before trusting it for anything at "
          "all."
       << endl;
}
 
namespace propt = boost::program_options;
 
void XtpMap::Initialize() {
 
  CSG::TrajectoryWriter::RegisterPlugins();
  CSG::TrajectoryReader::RegisterPlugins();
  CSG::TopologyReader::RegisterPlugins();
 
  AddProgramOptions()("topology,t", propt::value<string>(), "  topology");
  AddProgramOptions()("coordinates,c", propt::value<string>(),
                      "  coordinates or trajectory");
  AddProgramOptions()("segments,s", propt::value<string>(),
                      "  definition of segments and fragments");
  AddProgramOptions()("makesegments,m", "  write out a skeleton segments file");
  AddProgramOptions()("file,f", propt::value<string>(), "  state file");
  AddProgramOptions()("first-frame,i",
                      propt::value<votca::Index>()->default_value(0),
                      "  start from this frame");
  AddProgramOptions()("begin,b", propt::value<double>()->default_value(0.0),
                      "  start time in simulation");
  AddProgramOptions()("nframes,n",
                      propt::value<votca::Index>()->default_value(1),
                      "  number of frames to process");
  AddProgramOptions()(
      "guess-bonds",
      "  guess real bond connectivity from atom positions alone (a simple, "
      "van-der-Waals-radii-based heuristic, matching the one VMD itself "
      "uses for visualization -- restricted to atom pairs within the same "
      "real MD molecule only) -- ONLY if the loaded topology genuinely "
      "has no real bond connectivity of its own at all; never used if "
      "real bonds are already present. Writes a real, direct, inspectable "
      "report of every guessed bond to guessed_bonds_report.txt -- always "
      "review this before trusting the guessed connectivity for anything "
      "at all, since this heuristic can be genuinely wrong (see VMD's own "
      "developers' documented caveats about it).");
}
 
bool XtpMap::EvaluateOptions() {
 
  CheckRequired("topology", "Missing topology file");
  CheckRequired("segments", "Missing segment definition file");
  CheckRequired("coordinates", "Missing trajectory input");
  if (!(OptionsMap().count("makesegments"))) {
    CheckRequired("file", "Missing state file");
  }
  return 1;
}
 
void XtpMap::Run() {
 
  std::string name = ProgramName();
  if (VersionString() != "") {
    name = name + ", version " + VersionString();
  }
  XTP::HelpTextHeader(name);
 
  // ++++++++++++++++++++++++++++ //
  // Create MD topology from file //
  // ++++++++++++++++++++++++++++ //
 
  // Create topology reader
  string topfile = OptionsMap()["topology"].as<string>();
  std::unique_ptr<CSG::TopologyReader> topread =
      CSG::TopReaderFactory().Create(topfile);
 
  if (topread == nullptr) {
    throw runtime_error(string("Input format not supported: ") +
                        OptionsMap()["topology"].as<string>());
  }
  CSG::Topology mdtopol;
  topread->ReadTopology(topfile, mdtopol);
  if (votca::Log::verbose()) {
    cout << "Read MD topology from " << topfile << ": Found "
         << mdtopol.BeadCount() << " atoms in " << mdtopol.MoleculeCount()
         << " molecules. " << endl;
  }
 
  // ++++++++++++++++++++++++++++++ //
  // Create MD trajectory from file //
  // ++++++++++++++++++++++++++++++ //
 
  // Create trajectory reader and initialize
  string trjfile = OptionsMap()["coordinates"].as<string>();
  std::unique_ptr<CSG::TrajectoryReader> trjread =
      CSG::TrjReaderFactory().Create(trjfile);
 
  if (trjread == nullptr) {
    throw runtime_error(string("Input format not supported: ") +
                        OptionsMap()["coordinates"].as<string>());
  }
  trjread->Open(trjfile);
  trjread->FirstFrame(mdtopol);
 
  if (OptionsMap().count("guess-bonds")) {
    GuessBonds(mdtopol);
  }
 
  string mapfile = OptionsMap()["segments"].as<string>();
  if (OptionsMap().count("makesegments")) {
    if (TOOLS::filesystem::FileExists(mapfile)) {
      cout << endl
           << "xtp_map : map file '" << mapfile
           << "' already in use. Delete the current mapfile or specify a "
              "different name."
           << endl;
      return;
    }
 
    cout << " Writing template mapfile to " << mapfile << std::endl;
 
    TOOLS::Property mapfile_prop("topology", "", "");
    TOOLS::Property& molecules = mapfile_prop.add("molecules", "");
 
    std::map<std::string, const CSG::Molecule*> firstmolecule;
 
    std::map<std::string, votca::Index> molecule_names;
    for (const CSG::Molecule& mol : mdtopol.Molecules()) {
      if (!molecule_names.count(mol.getName())) {
        firstmolecule[mol.getName()] = &mol;
      }
      molecule_names[mol.getName()]++;
    }
    for (const auto& mol : molecule_names) {
      std::cout << "Found " << mol.second << " with name " << mol.first
                << std::endl;
    }
    for (const auto& mol : molecule_names) {
      TOOLS::Property& molecule = molecules.add("molecule", "");
      molecule.add("mdname", mol.first);
      TOOLS::Property& segments = molecule.add("segments", "");
      TOOLS::Property& segment = segments.add("segment", "");
      segment.add("name", "UPTOYOU_BUTUNIQUE");
      segment.add("qmcoords_n", "XYZFILE_GROUNDSTATE");
      segment.add("multipoles_n", "MPSFILE_GROUNDSTATE");
      segment.add("map2md", "WANTTOMAPTOMDGEOMETRY");
      segment.add("U_xX_nN_h", "REORG1_hole");
      segment.add("U_nX_nN_h", "REORG2_hole");
      segment.add("U_xN_xX_h", "REORG3_hole");
      TOOLS::Property& fragments = segment.add("fragments", "");
      TOOLS::Property& fragment = fragments.add("fragment", "");
      std::string atomnames = "";
      const CSG::Molecule* csgmol = firstmolecule[mol.first];
      std::vector<const CSG::Bead*> sortedbeads;
      sortedbeads.reserve(csgmol->BeadCount());
      for (const CSG::Bead* bead : csgmol->Beads()) {
        sortedbeads.push_back(bead);
      }
      std::sort(sortedbeads.begin(), sortedbeads.end(),
                [&](const CSG::Bead* b1, const CSG::Bead* b2) {
                  return b1->getId() < b2->getId();
                });
 
      for (const CSG::Bead* bead : sortedbeads) {
        atomnames += " " + std::to_string(bead->getResnr()) + ":" +
                     bead->getName() + ":" + std::to_string(bead->getId());
      }
      fragment.add("name", "UPTOYOU_BUTUNIQUE");
      fragment.add("mdatoms", atomnames);
      fragment.add("qmatoms", "IDS of QMATOMS i.e 0:C 1:H 2:C");
      fragment.add("mpoles", "IDS of MPOLES i.e 0:C 1:H 2:C");
      fragment.add("weights",
                   "weights for mapping(often atomic mass) i.e. 12  1 12");
      fragment.add("localframe", "IDs of up to 3 qmatoms or mpoles i.e. 0 1 2");
      std::ofstream template_mapfile(mapfile);
      template_mapfile << mapfile_prop << std::flush;
      template_mapfile.close();
 
      std::cout << "MOLECULETYPE " << csgmol->getName() << std::endl;
      std::cout << "SAMPLECOORDINATES" << std::endl;
      std::cout << "ID NAME COORDINATES[Angstroem] " << std::endl;
      for (const CSG::Bead* bead : sortedbeads) {
        Eigen::Vector3d pos = bead->getPos() * votca::tools::conv::nm2ang;
        std::string output =
            (boost::format("%1$i %2$s %3$+1.4f %4$+1.4f %5$+1.4f\n") %
             bead->getId() % bead->getName() % pos[0] % pos[1] % pos[2])
                .str();
        std::cout << output;
      }
    }
    std::cout << std::flush;
    return;
  }
 
  if (!TOOLS::filesystem::FileExists(mapfile)) {
    cout << endl
         << "xtp_map : map file '" << mapfile << "' could not be found."
         << endl;
    return;
  }
  XTP::Md2QmEngine md2qm(mapfile);
 
  votca::Index firstFrame = OptionsMap()["first-frame"].as<votca::Index>();
  votca::Index nFrames = OptionsMap()["nframes"].as<votca::Index>();
  bool beginAt = false;
  double time = OptionsMap()["begin"].as<double>();
  double startTime = mdtopol.getTime();
  if (time > 0.0) {
    beginAt = true;
    startTime = time;
  }
 
  // Extract first frame specified
  bool hasFrame;
  votca::Index frames_found = 0;
  votca::Index firstframecounter = firstFrame;
  for (hasFrame = true; hasFrame == true;
       hasFrame = trjread->NextFrame(mdtopol)) {
    frames_found++;
    if (((mdtopol.getTime() < startTime) && beginAt) || firstframecounter > 0) {
      firstframecounter--;
      continue;
    }
    break;
  }
  if (!hasFrame) {
    trjread->Close();
 
    throw runtime_error("Time or frame number exceeds trajectory length");
  }
  if (votca::Log::verbose()) {
    cout << "Read MD trajectory from " << trjfile << ": found " << frames_found
         << " frames, starting from frame " << firstFrame << endl;
  }
  // +++++++++++++++++++++++++ //
  // Convert MD to QM Topology //
  // +++++++++++++++++++++++++ //
 
  string statefile = OptionsMap()["file"].as<string>();
  if (TOOLS::filesystem::FileExists(statefile)) {
    cout << endl
         << "xtp_map : state file '" << statefile
         << "' already in use. Delete the current statefile or specify a "
            "different name."
         << endl;
    return;
  }
 
  XTP::StateSaver statsav(statefile);
  votca::Index laststep =
      -1;  // for some formats no step is given out so we check if the step
  for (votca::Index saved = 0; hasFrame && saved < nFrames;
       hasFrame = trjread->NextFrame(mdtopol), saved++) {
    if (mdtopol.getStep() == laststep) {
      mdtopol.setStep(laststep + 1);
    }
    laststep = mdtopol.getStep();
    XTP::Topology qmtopol = md2qm.map(mdtopol);
    statsav.WriteFrame(qmtopol);
  }
}
 
void XtpMap::ShowHelpText(std::ostream& out) {
  string name = ProgramName();
  if (VersionString() != "") {
    name = name + ", version " + VersionString();
  }
  XTP::HelpTextHeader(name);
  HelpText(out);
  out << "\n\n" << VisibleOptions() << endl;
}
 
int main(int argc, char** argv) {
  XtpMap xtpmap;
  return xtpmap.Exec(argc, argv);
}