fragmentsaturator.h

Go to the documentation of this file.

/*
 *            Copyright 2009-2026 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.
 *
 */
 
#pragma once
#ifndef VOTCA_XTP_FRAGMENTSATURATOR_H
#define VOTCA_XTP_FRAGMENTSATURATOR_H
 
// Local VOTCA includes
#include "qmmolecule.h"
 
namespace votca {
namespace xtp {

class FragmentSaturator {
 public:
  // A generic, "typical" C-H bond length (Angstrom) -- not specific to
  // any one hybridization/element (this class only ever adds H atoms,
  // never any other element, matching the "H-saturation" scope this
  // whole feature was designed around from the start). Deliberately
  // just a reasonable, generic default, not a carefully-chosen,
  // atom-specific value -- the caller can override this directly via
  // SaturateExternalBonds's own bond_length_angstrom parameter, and,
  // in any case, this initial placement is only ever a starting point
  // for a later, separate relaxation step, not the final, trusted
  // position.
  static constexpr double kDefaultCHBondLengthAngstrom = 1.09;
 
  // Direct answer to a real, genuine gap surfaced while designing
  // IPodCoupling::EvalJob, worked through directly with the user
  // before implementing this: PODCoupling's own constructor needs to
  // know which fragment (A or B) every atom belongs to, including
  // every new, saturating H atom -- but which fragment a given H
  // belongs to is NOT recoverable from its own final position or
  // index alone (it must inherit whichever ORIGINAL atom it is
  // saturating). new_atom_parent_ids is sized to mol's own,
  // resulting, saturated size (result.mol.size()), and directly
  // indexed by atom id: new_atom_parent_ids[i] is -1 for every atom
  // that was already in the original, input molecule (i <
  // original mol.size()), and the ORIGINAL atom's own id (i.e. the
  // atom whose own external bond this new H is saturating) for every
  // new H atom appended (i >= original mol.size()) -- deliberately
  // NOT "which fragment", since SaturateExternalBonds itself has no
  // concept of fragments at all (a generic, "saturate whatever
  // molecule you are given" method); the caller (which does know
  // which atom ids belong to which fragment) can translate this
  // directly, itself.
  struct SaturationResult {
    QMMolecule mol;
    std::vector<Index> new_atom_parent_ids;
  };
 
  // Returns a NEW QMMolecule (the input mol is never modified) with an
  // additional H atom appended for every atom in mol that has
  // hasExternalBond() == true, placed at bond_length_angstrom along
  // that atom's own, already-recorded, already-transformed external-
  // bond direction. Atoms without an external bond are copied
  // unchanged, at their own, original indices -- the new H atoms are
  // appended strictly after all of mol's own, original atoms, so
  // every original atom's own index is preserved exactly (matching
  // the same, direct requirement already confirmed, earlier this
  // session, to be necessary for fragment_A/fragment_B atom-index
  // definitions computed before this saturation step to remain valid
  // afterward). See SaturationResult's own comment above for the new
  // atoms' own parent-atom tracking.
  static SaturationResult SaturateExternalBonds(
      const QMMolecule& mol,
      double bond_length_angstrom = kDefaultCHBondLengthAngstrom);
  // Returns a NEW QMMolecule (the input mol is never modified) with a
  // constrained MMFF94 (or, if MMFF94 typing fails for some atom,
  // UFF -- see this method's own .cc implementation for exactly when
  // that fallback is used) relaxation applied: every atom with index
  // < n_original_atoms is held fixed at its own, current position;
  // every atom with index >= n_original_atoms (i.e. the new H atom(s)
  // SaturateExternalBonds appended, per its own, documented
  // convention of appending strictly after all original atoms) is
  // free to move. Genuinely needs mol's own, real, full internal bond
  // connectivity (via QMAtom::getBondedPartnerIds(), populated by the
  // SegmentMapper/Md2QmEngine pipeline already built earlier this
  // session) to build a real OpenBabel OBMol and correctly set up its
  // own force field at all -- confirmed directly, earlier this
  // session, that OpenBabel's own OBMol::PerceiveBondOrders() (used
  // here to derive bond order/hybridization from this known
  // connectivity plus real 3D geometry) needs that connectivity as
  // real input, it does not itself guess which atoms are bonded at
  // all -- and confirmed, separately, that geometry-only bond
  // perception (no known connectivity at all) is genuinely unreliable
  // for exactly this kind of situation (a newly-added atom), via a
  // real, earlier, direct RDKit test this same session that gave a
  // chemically wrong result for an analogous case.
  //
  // n_steps is the maximum number of conjugate-gradient steps
  // (OBForceField::ConjugateGradients's own argument, passed through
  // directly, unmodified).
  static QMMolecule RelaxNewAtoms(const QMMolecule& mol, Index n_original_atoms,
                                  Index n_steps = 500);
};
 
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_FRAGMENTSATURATOR_H