#pragma once
/*
* sketcherMinimizer.h
*
* Created by Nicola Zonta on 13/04/2010.
* Copyright Schrodinger, LLC. All rights reserved.
*
*/
#include <map>
#include <vector>
#include <iostream>
#include "CoordgenConfig.hpp"
#include "CoordgenFragmentBuilder.h"
#include "CoordgenMinimizer.h"
#include "sketcherMinimizerFragment.h"
#include "sketcherMinimizerMarchingSquares.h"
#include "sketcherMinimizerMolecule.h"
#include "sketcherMinimizerResidue.h"
#include "sketcherMinimizerResidueInteraction.h"
#include "sketcherMinimizerRing.h"
static const float SKETCHER_STANDARD_PRECISION = 1.f;
static const float SKETCHER_QUICK_PRECISION = 0.2f;
static const float SKETCHER_BEST_PRECISION = 3.f;
class sketcherMinimimizerInteraction;
namespace schrodinger
{
namespace mae
{
class Block;
}
} // namespace schrodinger
typedef struct {
std::vector<sketcherMinimizerPointF> additionVectors;
std::vector<sketcherMinimizerPointF> centers;
std::vector<int> counters;
} proximityData;
/* class to handle templates of common difficult ring structures */
class CoordgenTemplates
{
public:
CoordgenTemplates() = default;
~CoordgenTemplates()
{
for (auto molecule : m_templates) {
for (auto atom : molecule->_atoms) {
delete atom;
}
for (auto bond : molecule->_bonds) {
delete bond;
}
delete molecule;
}
m_templates.clear();
}
std::vector<sketcherMinimizerMolecule*>& getTemplates()
{
return m_templates;
}
void setTemplateDir(std::string&& dir)
{
m_templateDir = std::move(dir);
if (m_templateDir.back() != '/') {
m_templateDir += "/";
}
}
std::string getTemplateDir() { return m_templateDir; }
private:
std::vector<sketcherMinimizerMolecule*> m_templates;
std::string m_templateDir = "";
};
/* main class. Creates 2d coordinates for molecular inputs */
class EXPORT_COORDGEN sketcherMinimizer
{
public:
sketcherMinimizer(float precision = SKETCHER_STANDARD_PRECISION);
~sketcherMinimizer();
/* fields used in rdkit */
std::vector<sketcherMinimizerFragment*> _fragments;
CoordgenMinimizer m_minimizer;
/* run coordinates generation and return true if the pose is considered
* optimal */
bool runGenerateCoordinates();
/* initialize data from given molecule */
void initialize(sketcherMinimizerMolecule* minMol);
/* run a force-field based minimization on the given molecule */
void minimizeMolecule(sketcherMinimizerMolecule* molecule);
// void initializeFromMolecule(ChmMol& mol);
void writeMinimizationData();
/* run a force-field based minimization */
void forceFieldMinimize();
/* if mol contains separate molecules, split them into a vector */
void splitIntoMolecules(sketcherMinimizerMolecule* mol,
std::vector<sketcherMinimizerMolecule*>& mols);
/* split molecules into rigid fragments */
void findFragments();
/* constrain coordinates on all atoms */
void constrainAllAtoms();
/* constrain coordinates on atoms corresponding to true */
void constrainAtoms(const std::vector<bool>& constrained);
/* fix cooordinates (i.e. guarantee they will not change) on atoms marked as
* true */
void fixAtoms(const std::vector<bool>& fixed);
/* set a flag to enable/disable the scoring of interactions with residues */
void setScoreResidueInteractions(bool b);
/* add bond to m_extrabonds */
void addExtraBond(sketcherMinimizerBond* bond);
/* getters */
const std::vector<sketcherMinimizerAtom*>& getAtoms() const { return m_atoms; }
const std::vector<sketcherMinimizerMolecule*>& getMolecules() const { return m_molecules; }
const std::vector<sketcherMinimizerAtom*>& getReferenceAtoms() const { return m_referenceAtoms; }
const std::vector<sketcherMinimizerBond*>& getReferenceBonds() const { return m_referenceBonds; }
const std::vector<sketcherMinimizerFragment*>& getIndependentFragments() const { return m_independentFragments; }
bool getTreatNonterminalBondsToMetalAsZOBs() const {return m_treatNonterminalBondsToMetalAsZOBs;}
const std::vector<sketcherMinimizerBendInteraction*>& getBendInteractions() const;
const std::vector<sketcherMinimizerStretchInteraction*>& getStretchInteractions() const;
const std::vector<sketcherMinimizerInteraction*>& getInteractions() const;
/* setters */
void setFragments(std::vector<sketcherMinimizerFragment*> fragments) { _fragments = fragments; }
void setTreatNonterminalBondsToMetalAsZOBs(bool b) {m_treatNonterminalBondsToMetalAsZOBs = b;}
void setEvenAngles(bool b);
void setSkipMinimization(bool b);
void setForceOpenMacrocycles(bool b);
/// Interactions with CoordgenMinimizer and CoordgenFragmentBuilder
/* find a list of carbons from the backbone C=O of a peptide */
std::set<sketcherMinimizerAtom*>
getChetoCs(const std::vector<sketcherMinimizerAtom*>& allAtoms) const;
/* find a list of nitrogens from the backbon NH of a peptide */
std::set<sketcherMinimizerAtom*>
getAminoNs(const std::vector<sketcherMinimizerAtom*>& allAtoms) const;
/* find a list of alpha carbons of a peptide */
std::set<sketcherMinimizerAtom*>
getAlphaCs(const std::vector<sketcherMinimizerAtom*>& allAtoms,
const std::set<sketcherMinimizerAtom*>& chetoCs,
const std::set<sketcherMinimizerAtom*>& aminoNs) const;
/* clear all the interactions loaded in the minimizer and free memory */
void clearInteractions();
/* setup all constraints of given molecule */
void addInteractionsOfMolecule(sketcherMinimizerMolecule* molecule,
bool intrafragmentClashes = false);
void addExtraInteraction(sketcherMinimizerMolecule* molecule,
sketcherMinimizerInteraction* interaction);
void buildFromFragments(bool b);
bool avoidClashesOfMolecule(
sketcherMinimizerMolecule* molecule,
const std::vector<sketcherMinimizerInteraction*>& extraInteractions =
std::vector<sketcherMinimizerInteraction*>());
/// static members
/* put atoms in a canonical order to reduce dependency from order in the
* input vector */
static void canonicalOrdering(sketcherMinimizerMolecule* minMol);
/* pick one atom out of the vector. Arbitrary criteria such as atomic number
* and connectivity are used */
static sketcherMinimizerAtom*
pickBestAtom(std::vector<sketcherMinimizerAtom*>& atoms);
/* return a score of the alignment between direction and templat.first,
* weight on the angle between the two and templat.second */
static float
testAlignment(const sketcherMinimizerPointF& direction,
const std::pair<sketcherMinimizerPointF, float>& templat);
/* find the best alignment of a fragment to its parent and set invert in
* case the fragment needs to be flipped */
static sketcherMinimizerPointF scoreDirections(
sketcherMinimizerFragment* fragment, float angle,
const std::vector<std::pair<sketcherMinimizerPointF, float>>&
directions,
bool& invert);
/* align the fragment to its parent */
static void
alignWithParentDirection(sketcherMinimizerFragment* f,
const sketcherMinimizerPointF& position,
float angle);
/* align the fragment to its parent in the case of constrained coordinates
*/
static bool
alignWithParentDirectionConstrained(sketcherMinimizerFragment* fragment,
const sketcherMinimizerPointF& position,
float angle);
/* align the fragment to its parent in the case of unconstrained coordinates
*/
static bool
alignWithParentDirectionUnconstrained(sketcherMinimizerFragment* fragment,
float angle);
/* get all bonds to a terminal atom */
static std::vector<sketcherMinimizerBond*>
getAllTerminalBonds(sketcherMinimizerFragment* fragment);
/* return a list of vectors the given fragment can be aligned with and a
* score of the importance of each */
static std::vector<std::pair<sketcherMinimizerPointF, float>>
findDirectionsToAlignWith(sketcherMinimizerFragment* fragment);
/* if the three atoms share a ring, return it */
static sketcherMinimizerRing* sameRing(const sketcherMinimizerAtom* at1,
const sketcherMinimizerAtom* at2,
const sketcherMinimizerAtom* at3);
/* if the two atoms share a ring, return it */
static sketcherMinimizerRing* sameRing(const sketcherMinimizerAtom* at1,
const sketcherMinimizerAtom* at2);
/* if the two atoms share a bond, return it */
static sketcherMinimizerBond* getBond(const sketcherMinimizerAtom* a1,
const sketcherMinimizerAtom* a2);
/* calculate root mean square deviation between templates and points */
static float RMSD(const std::vector<sketcherMinimizerPointF>& templates,
const std::vector<sketcherMinimizerPointF>& points);
/* singular value decomposition for 2x2 matrices.
used for 2D alignment. */
static void svd(float* a, float* U, float* Sig, float* V);
/* set m to a rotation matrix to align ref to points */
static void
alignmentMatrix(const std::vector<sketcherMinimizerPointF>& ref,
const std::vector<sketcherMinimizerPointF>& points,
float* m);
static void
checkIdentity(std::vector<unsigned int> solution, int newSol,
std::vector<bool>& matrix,
std::vector<sketcherMinimizerPointF>& templateCoordinates,
std::vector<std::vector<size_t>>& molBonds,
std::vector<std::vector<size_t>>& templateBonds,
std::vector<std::vector<size_t>>& molCisTransChains,
std::vector<bool>& molIsCis, size_t size, bool& found,
std::vector<unsigned int>& mapping);
/* compare atoms and bonds to template and map which atom is which in case
* of a positive match */
static bool compare(const std::vector<sketcherMinimizerAtom*>& atoms,
const std::vector<sketcherMinimizerBond*>& bonds,
sketcherMinimizerMolecule* templ,
std::vector<unsigned int>& mapping);
/* calculate morgan scores for the given input */
static int morganScores(const std::vector<sketcherMinimizerAtom*>& atoms,
const std::vector<sketcherMinimizerBond*>& bonds,
std::vector<int>& scores);
/* load the templates from the template file */
static void setTemplateFileDir(std::string dir);
static void loadTemplates();
static CoordgenTemplates m_templates;
private:
/*all non-terminal bonds to a metal atom are treated as if they were zero order bonds (this usually results
in a longer bond*/
bool m_treatNonterminalBondsToMetalAsZOBs = true;
CoordgenFragmentBuilder m_fragmentBuilder;
std::vector<sketcherMinimizerAtom*> m_atoms;
std::vector<sketcherMinimizerAtom*> m_referenceAtoms;
std::vector<sketcherMinimizerResidue*> m_residues;
std::vector<sketcherMinimizerResidueInteraction*> m_residueInteractions;
std::vector<sketcherMinimizerFragment*> m_independentFragments;
std::vector<sketcherMinimizerBond*> m_bonds;
std::vector<sketcherMinimizerBond*> m_referenceBonds;
std::vector<sketcherMinimizerBond*> m_proximityRelations;
std::vector<sketcherMinimizerBond*> m_extraBonds;
std::vector<sketcherMinimizerMolecule*> m_molecules;
/* return the position of res, which is part of SSE, given that the first
* residue of SSE is placed at startF and consecutive residues are placed
* increment away from each other. All distances are expressed in floats,
* where 0.f is an arbitrary starting point, 0.5 is the opposite side of the
* curve and 1.0 is again the starting point */
float getResidueDistance(float startF, float increment,
sketcherMinimizerResidue* res,
const std::vector<sketcherMinimizerResidue*>& SSE) const;
/* return the vector index corresponding to floatPosition */
int getShapeIndex(const std::vector<sketcherMinimizerPointF>& shape,
float floatPosition) const;
/* return true if the molecules structure is reasonable (e.g. reasonable
* amount of fused rings) */
bool structurePassSanityCheck() const;
/* clear data and free memory */
void clear();
/* flag atoms that will be drawn with 90° angles (e.g. phosphate P) */
void flagCrossAtoms();
/* assign coordinates to all molecules and residues */
void minimizeAll();
/* find the best angle to rotate each molecule */
void bestRotation();
/* add info to choose the best angle so that, if present, peptide chains are
* horizontal */
void addBestRotationInfoForPeptides(
std::vector<std::pair<float, float>>& angles,
const std::vector<sketcherMinimizerAtom*>& atoms);
/* if a peptide chain is present make sure that the N term is on the left */
void maybeFlipPeptides(const std::vector<sketcherMinimizerAtom*>& atoms,
float& scoreX);
/* consider flipping the molecule horizontally and/or vertically */
void maybeFlip();
/* mark atoms with wedges as above or below the plane to correctly draw
* crossing bonds */
void assignPseudoZ();
/* write wedges and dashed bonds to mark stereochemistry */
void writeStereoChemistry();
/* arrange multiple molecules next to each other */
void arrangeMultipleMolecules();
/* arrange molecules that have parts that interact with each other so that
* they are close */
void placeMoleculesWithProximityRelations(
std::vector<sketcherMinimizerMolecule*> proximityMols);
/* place molecules so that one is in the middle and other are around */
void placeMolResidueLigandStyle(sketcherMinimizerMolecule* mol,
sketcherMinimizerMolecule* parent);
/* if the molecule has more than one interaction and they cross mirror its
coordinates so they don't cross anymore */
void flipIfCrossingInteractions(sketcherMinimizerMolecule* mol);
/* build data vectors to place molecules with proximity relations */
std::vector<proximityData> buildProximityDataVector(
std::vector<sketcherMinimizerMolecule*>& proximityMols,
std::map<sketcherMinimizerMolecule*, sketcherMinimizerAtom*>& molMap);
/* translate molecules with proximity relations */
void translateMoleculesWithProximityRelations(
std::vector<sketcherMinimizerMolecule*>& proximityMols,
std::map<sketcherMinimizerMolecule*, sketcherMinimizerAtom*>& molMap,
std::map<sketcherMinimizerMolecule*, sketcherMinimizerPointF>&
templateCenters,
std::vector<proximityData>& proximityDataVecto);
/* rotate molecules with proximity relations */
void rotateMoleculesWithProximityRelations(
std::vector<sketcherMinimizerMolecule*>& proximityMols,
std::map<sketcherMinimizerMolecule*, sketcherMinimizerAtom*>& molMap,
std::vector<proximityData>& proximityDataVector);
/* place residues in concentric contours around the ligand */
void placeResiduesInCrowns();
/* place residues from the given strands of consecutive residues to fill
the given path */
bool fillShape(std::vector<std::vector<sketcherMinimizerResidue*>>& SSEs,
const std::vector<sketcherMinimizerPointF>& shape,
int shapeN);
/* place a single strand of consecutive residues */
void placeSSE(const std::vector<sketcherMinimizerResidue*>& SSE,
const std::vector<sketcherMinimizerPointF>& shape, int shapeN,
std::vector<bool>& penalties,
std::set<sketcherMinimizerResidue*>& outliers,
bool placeOnlyInteracting = false);
/* score the position of the given strands */
float scoreSSEPosition(const std::vector<sketcherMinimizerResidue*>& SSE,
const std::vector<sketcherMinimizerPointF>& shape,
int shapeN, std::vector<bool>& penalties, float startingPosition,
float increment);
/* score the distance between the two given points of connected residues */
float scoreSSEBondStretch(const sketcherMinimizerPointF& coordinates1,
const sketcherMinimizerPointF& coordinates2);
/* solution represent the placement chosen for residues in SSE. Mark the
* corresponding sections of the crown to prevent other residues to be
* placed
* there */
void markSolution(const std::pair<float, float>& solution,
const std::vector<sketcherMinimizerResidue*>& SSE,
const std::vector<sketcherMinimizerPointF>& shape,
std::vector<bool>& penalties,
std::set<sketcherMinimizerResidue*>& outliers);
/* return a concentric shape around the ligand. CrownN controls how far away
from the ligand the shape is */
std::vector<sketcherMinimizerPointF> shapeAroundLigand(int crownN);
/* group residues in strands of consecutive residues */
std::vector<std::vector<sketcherMinimizerResidue*>>
groupResiduesInSSEs(const std::vector<sketcherMinimizerResidue*>& residues);
/* score the position of given residues */
float
scoreResiduePosition(int index,
const std::vector<sketcherMinimizerPointF>& shape,
int shapeN, std::vector<bool>& penalties,
sketcherMinimizerResidue* residue);
/* assign coordinates to residues */
void placeResidues(const std::vector<sketcherMinimizerAtom*>& atoms =
std::vector<sketcherMinimizerAtom*>(0));
/* assign coordinates to residues in the context of a protein-protein
interaction diagram */
void placeResiduesProteinOnlyMode();
/* assign coordinates to residues in a protein-protein interaction
diagram shaped as a circle */
void placeResiduesProteinOnlyModeCircleStyle(
const std::map<std::string, std::vector<sketcherMinimizerResidue*>>&
chains);
/* assign coordinates to residues in a protein-protein interaction
diagram shaped as a LID */
void placeResiduesProteinOnlyModeLIDStyle(
const std::map<std::string, std::vector<sketcherMinimizerResidue*>>&
chains);
/* order residues for drawing, so that residues interacting together are
* drawn one after the other and residues with more interactions are drawn
* first */
std::vector<sketcherMinimizerResidue*> orderResiduesOfChains(
const std::map<std::string, std::vector<sketcherMinimizerResidue*>>&
chains);
/* find center position for each chain of residues using a meta-molecule
* approach, building a molecule where each atom represents a chain and each
* bond connects two interacting chains */
std::map<std::string, sketcherMinimizerPointF>
computeChainsStartingPositionsMetaMol(
const std::map<std::string, std::vector<sketcherMinimizerResidue*>>&
chains);
/* place interacting residues closer to each other, so they end up at the
* periphery of the chain */
void shortenInteractions(
const std::map<std::string, std::vector<sketcherMinimizerResidue*>>&
chains);
/* explore positions in a grid around the current one to ease clashes */
sketcherMinimizerPointF exploreGridAround(
const sketcherMinimizerPointF& centerOfGrid, unsigned int levels,
float gridD, float dx = 0.f, float dy = 0.f,
float distanceFromAtoms = -1.f, bool watermap = false,
sketcherMinimizerResidue* residueForInteractions = nullptr,
const sketcherMinimizerPointF& direction = sketcherMinimizerPointF(0,
1));
sketcherMinimizerPointF exploreMolPosition(sketcherMinimizerMolecule* mol,
unsigned int levels, float gridD,
float distanceFromAtoms = -1.f);
/* add given angle to a vector of angles, clustering them if a close angle
* is already in the vector */
void addToVector(float weight, float angle,
std::vector<std::pair<float, float>>& angles);
void assignLongestChainFromHere(sketcherMinimizerFragment* f);
void assignNumberOfChildrenAtomsFromHere(sketcherMinimizerFragment* f);
/* for each residue, find the closest atom among atoms, or all atoms
if none are given */
void
findClosestAtomToResidues(const std::vector<sketcherMinimizerAtom*>& atoms =
std::vector<sketcherMinimizerAtom*>(0));
/* initialize data and coordinates for each fragment */
void initializeFragments();
/* for tracking coordinates */
float sin_flip = 0.f;
float cos_flip = 0.f;
float centerX = 0.f;
float centerY = 0.f;
int flipX = 1;
int flipY = 1;
};