1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
|
// ----------------------------------------------------
// $Maintainer: Marcel Schumann $
// $Authors: Marcel Schumann $
// ----------------------------------------------------
#include <BALL/SCORING/COMPONENTS/PB.h>
#include <BALL/SCORING/COMMON/scoringFunction.h>
#include <BALL/SYSTEM/timer.h>
#include <BALL/KERNEL/PTE.h>
using namespace BALL;
using namespace std;
PB::PB(ScoringFunction& sf)
: ScoringComponent(sf)
{
pb_solver_ = new FDPB;
// FDPB uses Atom::getRadius() instead of Element::getVanDerWaalsRadius(), for whatever reason.
// However, Atom::radius_ is not set by default, so we have to do it here.
for (AtomIterator it = scoring_function_->getLigand()->beginAtom(); +it; it++)
{
it->setRadius(it->getElement().getVanDerWaalsRadius());
}
for (AtomIterator it = scoring_function_->getReceptor()->beginAtom(); +it; it++)
{
it->setRadius(it->getElement().getVanDerWaalsRadius());
}
}
PB::~PB()
{
delete pb_solver_;
}
void PB::setupLigand()
{
// FDPB uses Atom::getRadius() instead of Element::getVanDerWaalsRadius(), for whatever reason.
// However, Atom::radius_ is not set by default, so we have to do it here.
for (AtomIterator it = scoring_function_->getLigand()->beginAtom(); +it; it++)
{
it->setRadius(it->getElement().getVanDerWaalsRadius());
}
}
void PB::update(const vector<std::pair<Atom*, Atom*> >& pair_vector)
{
// Clear-functions will automatically erase all data that was dynamically allocated by the last call of this function
receptor_atoms_.clear();
ligand_atoms_.clear();
complex_atoms_.clear();
Molecule* receptor_mol = new Molecule;
Molecule* ligand_mol = new Molecule;
Molecule* complex_mol = new Molecule;
receptor_atoms_.insert(*receptor_mol);
ligand_atoms_.insert(*ligand_mol);
complex_atoms_.insert(*complex_mol);
// At the moment we need to copy the relevant atoms and insert them into a dummy-system, because FDPB can only work on Systems (not an atom-list).
set<Atom*> copied_atoms;
for (vector < std::pair < Atom*, Atom* > > ::const_iterator it = pair_vector.begin(); it != pair_vector.end(); it++)
{
if (copied_atoms.find(it->first) == copied_atoms.end())
{
ligand_mol->insert(*new Atom(*it->first));
complex_mol->insert(*new Atom(*it->first));
copied_atoms.insert(it->first);
}
if (copied_atoms.find(it->second) == copied_atoms.end())
{
receptor_mol->insert(*new Atom(*it->second));
complex_mol->insert(*new Atom(*it->second));
copied_atoms.insert(it->second);
}
}
}
double PB::updateScore()
{
Timer timer;
timer.start();
Options options = scoring_function_->getOptions();
// calculate energy for receptor in vacuum
options.set(FDPB::Option::SOLVENT_DC, 1);
pb_solver_->setup(receptor_atoms_, options);
pb_solver_->solve();
double rec_vac = pb_solver_->getEnergy();
// calculate energy for solvated receptor
options.set(FDPB::Option::SOLVENT_DC, 78);
pb_solver_->setup(receptor_atoms_, options);
pb_solver_->solve();
double rec_solv = pb_solver_->getEnergy();
// calculate energy for ligand in vacuum
options.set(FDPB::Option::SOLVENT_DC, 1);
pb_solver_->setup(ligand_atoms_, options);
pb_solver_->solve();
double lig_vac = pb_solver_->getEnergy();
// calculate energy for solvated ligand
options.set(FDPB::Option::SOLVENT_DC, 78);
pb_solver_->setup(ligand_atoms_, options);
pb_solver_->solve();
double lig_solv = pb_solver_->getEnergy();
// calculate energy for complex in vacuum
options.set(FDPB::Option::SOLVENT_DC, 1);
pb_solver_->setup(complex_atoms_, options);
pb_solver_->solve();
double com_vac = pb_solver_->getEnergy();
// calculate energy for solvated complex
options.set(FDPB::Option::SOLVENT_DC, 78);
pb_solver_->setup(complex_atoms_, options);
pb_solver_->solve();
double com_solv = pb_solver_->getEnergy();
score_ = (com_solv-com_vac)-(rec_solv-rec_vac)-(lig_solv-lig_vac);
//scaleScore();
cout<<" rec_vac = "<<rec_vac<<endl;
cout<<" rec_solv = "<<rec_solv<<endl;
cout<<" lig_vac = "<<lig_vac<<endl;
cout<<" lig_solv = "<<lig_solv<<endl;
cout<<" com_vac = "<<com_vac<<endl;
cout<<" com_solv = "<<com_solv<<endl;
cout<<" ----------------------"<<endl;
cout<<" total PB score = "<< getScaledScore() <<endl<<flush;
timer.stop();
Log<<"Time used by PB::updateScore(): "<<timer.getClockTime()<<endl;
//return score_;
return getScaledScore();
}
|