18 #include <basic/options/option.hh>
30 #include <numeric/xyz.functions.hh>
31 #include <numeric/statistics/functions.hh>
40 #include <basic/options/keys/edensity.OptionKeys.gen.hh>
43 #include <basic/Tracer.hh>
46 #include <utility/vector1.hh>
53 #define _USE_MATH_DEFINES
60 namespace electron_density {
80 using namespace core::scoring::methods;
82 static THREAD_LOCAL basic::Tracer
TR(
"core.scoring.electron_density.ElecDensEnergy" );
136 utility_exit_with_message(
"Density scoring function called but no map loaded.");
141 int virt_res_idx = root_edge.
start();
145 if ( root_res.type().name() !=
"VRT" || root_edge.label() < 0 ) {
146 utility_exit_with_message(
"Fold tree is not set properly for density scoring!");
153 bool create_new_lre_container(
false );
155 if ( energies.long_range_container( lr_type ) == 0 ) {
156 create_new_lre_container =
true;
159 OneToAllEnergyContainerOP dec( utility::pointer::static_pointer_cast< core::scoring::OneToAllEnergyContainer > ( lrc ) );
161 if ( dec->size() != pose.
total_residue() || dec->fixed() != virt_res_idx ) {
162 create_new_lre_container =
true;
166 if ( create_new_lre_container ) {
167 TR.Debug <<
"Creating new LRE container (" << pose.
total_residue() <<
")" << std::endl;
169 energies.set_long_range_container( lr_type, new_dec );
186 TR.Debug <<
"ElecDensCenEnergy::setup_for_scoring() returns CC = " <<
structure_score << std::endl;
230 using namespace numeric::statistics;
236 Real z_CC = cc / 0.1;
237 Real p_null = 0.5 * errfc( z_CC/sqrt(2.0) );
238 Real edensScore = log ( p_null );
255 using namespace numeric::statistics;
261 int resid =
id.rsd();
262 int atmid =
id.atomno();
268 numeric::xyzVector<core::Real>
X = pose.
xyz(
id);
269 numeric::xyzVector< core::Real > dCCdx(0,0,0);
270 numeric::xyzMatrix< core::Real >
R = numeric::xyzMatrix<core::Real>::rows(1,0,0, 0,1,0, 0,0,1);
281 core::Size nres_per = symminfo->num_independent_residues();
282 bool remapSymm = basic::options::option[ basic::options::OptionKeys::edensity::score_symm_complex ]();
287 if ( atmid != 2 )
return;
291 int nchildren = edges_i.size();
292 if ( nchildren < 2 )
return;
298 if ( ! symminfo->jump_is_independent( edge_incoming.
label() ) ) {
308 utility::vector1< core::Size > vrtclones;
309 for (
int j=1; j<=nchildren; ++j ) {
310 int basejump = edges_i[j].
label();
311 if ( ! symminfo->jump_is_independent( basejump ) ) {
312 basejump = symminfo->jump_follows( basejump );
314 utility::vector1< core::Size > jumpclones = symminfo->jump_clones( basejump );
315 for (
int k=0; k<=(int)jumpclones.size(); ++k ) {
318 if ( edges_j.size() > 1 && std::find( vrtclones.begin(), vrtclones.end(), upstream ) == vrtclones.end() ) {
319 vrtclones.push_back( upstream );
325 for (
int i=1; i<=(int)vrtclones.size(); ++i ) {
329 for (
int j=1; j<=nchildren; ++j ) {
330 int downstream = edges_i[j].stop();
331 utility::vector1<int> mapping_j;
332 numeric::xyzMatrix< core::Real > R_j;
335 for (
int k=1; k<=(int)nsubunits; ++k ) {
336 if ( mapping_j[k] == 0 )
continue;
339 for (
int l=1, l_end=nres_per; l<=l_end; ++l ) {
341 int source_res = (k-1)*nres_per+l;
342 int target_res = (mapping_j[k]-1)*nres_per+l;
346 numeric::xyzVector<core::Real> X_lm_src = pose.
residue(source_res).
atom(2).
xyz();
347 numeric::xyzVector<core::Real> X_lm_tgt = pose.
residue(target_res).
atom(2).
xyz();
351 Real z_CC = CC / 0.1;
352 Real p_null = 0.5 * errfc( z_CC/sqrt(2.0) );
353 numeric::xyzVector< core::Real > dEdx = 0.5 * ( 1.0 / p_null ) *
355 exp(-
SQ( z_CC/sqrt(2.0) )) *
359 numeric::xyzVector<core::Real> atom_x = X_lm_tgt;
360 numeric::xyzVector<core::Real>
const f2( dEdx );
361 numeric::xyzVector<core::Real> atom_y = -f2 + atom_x;
362 Vector const f1( atom_x.cross( atom_y ) );
371 utility::vector1<int> mapping_i;
372 numeric::xyzMatrix< core::Real > R_i;
374 for (
int k=1; k<=(int)nsubunits; ++k ) {
375 if ( mapping_i[k] == 0 )
continue;
378 for (
int l=1, l_end=nres_per; l<=l_end; ++l ) {
380 int source_res = (k-1)*nres_per+l;
381 int target_res = (mapping_i[k]-1)*nres_per+l;
385 numeric::xyzVector<core::Real> X_lm_src = pose.
residue(source_res).
atom(2).
xyz();
386 numeric::xyzVector<core::Real> X_lm_tgt = pose.
residue(target_res).
atom(2).
xyz();
390 Real z_CC = CC / 0.1;
391 Real p_null = 0.5 * errfc( z_CC/sqrt(2.0) );
392 numeric::xyzVector< core::Real > dEdx = 0.5 * ( 1.0 / p_null ) *
394 exp(-
SQ( z_CC/sqrt(2.0) )) *
398 numeric::xyzVector<core::Real> atom_x = X_lm_tgt;
399 numeric::xyzVector<core::Real>
const f2( dEdx );
400 numeric::xyzVector<core::Real> atom_y = -f2 + atom_x;
401 Vector const f1( atom_x.cross( atom_y ) );
409 if ( ! symminfo->bb_is_independent( resid ) )
return;
412 utility::vector1< Size > myClones = symminfo->bb_clones(resid);
413 for (
int i=0; i<=(int)myClones.size(); ++i ) {
414 numeric::xyzVector<core::Real> X_i = (i==0) ?
X : pose.
xyz(
id::AtomID( 2, myClones[i] ) );
421 Real z_CC = CC / 0.1;
422 Real p_null = 0.5 * errfc( z_CC/sqrt(2.0) );
425 numeric::xyzVector< core::Real > dEdx = 0.5 * ( 1.0 / p_null ) *
427 exp(-
SQ( z_CC/sqrt(2.0) )) *
431 numeric::xyzVector<core::Real> atom_x =
X;
432 numeric::xyzVector<core::Real>
const f2( dEdx );
433 numeric::xyzVector<core::Real> atom_y = -f2 + atom_x;
434 Vector const f1( atom_x.cross( atom_y ) );
441 Real z_CC = CC / 0.1;
442 Real p_null = 0.5 * errfc( z_CC/sqrt(2.0) );
445 numeric::xyzVector< core::Real > dEdx = ( 1.0 / p_null ) *
448 exp(-
SQ( z_CC/sqrt(2.0) )) *
453 numeric::xyzVector<core::Real> atom_x =
X;
454 numeric::xyzVector<core::Real>
const f2( dEdx );
455 numeric::xyzVector<core::Real> atom_y = -f2 + atom_x;
456 Vector const f1( atom_x.cross( atom_y ) );
469 Real z_CC = CC / 0.1;
470 Real p_null = 0.5 * errfc( z_CC/sqrt(2.0) );
475 numeric::xyzVector< core::Real > dEdx = ( 1.0 / p_null ) *
478 exp(-
SQ( z_CC/sqrt(2.0) )) *
483 numeric::xyzVector<core::Real> atom_x =
X;
484 numeric::xyzVector<core::Real>
const f2( dEdx );
485 numeric::xyzVector<core::Real> atom_y = -f2 + atom_x;
486 Vector const f1( atom_x.cross( atom_y ) );
void get_symmMap(int vrtid, utility::vector1< int > &X_map, numeric::xyzMatrix< core::Real > &R)
virtual void setup_for_derivatives(pose::Pose &pose, ScoreFunction const &sf) const
Called immediately before atom- and DOF-derivatives are calculated allowing the derived class a chanc...
virtual void eval_atom_derivative(id::AtomID const &id, pose::Pose const &pose, kinematics::DomainMap const &, ScoreFunction const &sfxn, EnergyMap const &weights, Vector &F1, Vector &F2) const
called during gradient-based minimization inside dfunc
bool is_root(int const seqpos) const
Returns true if the the root.
kinematics::FoldTree const & fold_tree() const
Returns the pose FoldTree.
int start() const
start vertex, return by value
utility functions for handling with symmetric conformations
void set_nres(int nres)
set # of residues
This object defines a ScoreFunction, it contains methods for calculating the various scoring componen...
Edge const & jump_edge(int const jump_number) const
Returns the jump edge with jump number (const)
virtual void residue_pair_energy(conformation::Residue const &rsd1, conformation::Residue const &rsd2, pose::Pose const &pose, ScoreFunction const &sfxn, EnergyMap &emap) const
Evaluate the interaction between a given residue pair accumulating the unweighted energies in an Ener...
A cached energies object.
core::Real SQ(core::Real N)
A molecular system including residues, kinematics, and energies.
void get_R(int subunit, numeric::xyzMatrix< core::Real > &R)
virtual ScoreTypes score_types_for_method() const
Return the set of score types claimed by the EnergyMethod this EnergyMethodCreator creates in its cre...
Conformation const & conformation() const
Returns the pose Conformation (const-access)
utility::vector1< Edge > get_outgoing_edges(int const seqpos) const
Returns all edges that build a residue directly off of
virtual core::Size version() const
Return the version of the energy method.
methods::LongRangeEnergyType long_range_type() const
Size total_residue() const
Returns the number of residues in the pose conformation.
const_iterator begin() const
begin iterator of the edge_list
Declaration for the class that connects ElecDensCenEnergy with the ScoringManager.
PointPosition const & xyz(AtomID const &id) const
Returns the location (xyz) of pose AtomID
ObjexxFCL::FArray1D_int DomainMap
ElectronDensity & getDensityMap(std::string filename, bool force_reload)
The EDM instance.
utility::pointer::shared_ptr< EnergyMethod > EnergyMethodOP
virtual methods::EnergyMethodOP create_energy_method(methods::EnergyMethodOptions const &) const
Instantiate a new ElecDensCenEnergy.
utility::vector1< ScoreType > ScoreTypes
int label() const
label (edge type), return by value
Atom identifier class. Defined by the atom number and the residue number.
static THREAD_LOCAL basic::Tracer TR("core.scoring.methods.ElecDensEnergy")
utility::pointer::shared_ptr< LREnergyContainer > LREnergyContainerOP
numeric::xyzVector< Length > Vector
virtual methods::EnergyMethodOP clone() const
clone
void dCCdx_cen(int resid, numeric::xyzVector< core::Real > const &X, core::pose::Pose const &pose, numeric::xyzVector< core::Real > &gradX)
Return the gradient of CC w.r.t. res X's CA's movement Centroid-mode analogue of dCCdx.
A container interface for storing and scoring long range energies.
utility::pointer::shared_ptr< OneToAllEnergyContainer > OneToAllEnergyContainerOP
virtual void finalize_total_energy(pose::Pose const &pose, ScoreFunction const &, EnergyMap &totals) const
called at the end of energy evaluation
virtual void eval_intrares_energy(conformation::Residue const &rsd, pose::Pose const &pose, ScoreFunction const &sfxn, EnergyMap &emap) const
Evaluate the intra-residue constraint energy for a given residue.
Method declarations and simple accessor definitions for the Residue class.
bool is_symmetric(scoring::Energies const &energies)
utility::pointer::shared_ptr< EnergyMethodCreator > EnergyMethodCreatorOP
A vector for storing energy data, initially all values are 0.
Residue const & residue(Size const seqpos) const
Returns the Residue at position (read access) Note: this method will trigger a refold if eit...
void compute_symm_rotations(core::pose::Pose const &pose, core::conformation::symmetry::SymmetryInfoCOP symmInfo=NULL)
Computes the symmatric rotation matrices.
core::Real matchCentroidPose(core::pose::Pose const &pose, core::conformation::symmetry::SymmetryInfoCOP symmInfo=NULL, bool cacheCCs=false)
Quickly matches a centroid pose into a low-resolution density map by placing a single Gaussian at eac...
Energy graph class declaration.
virtual bool defines_residue_pair_energy(pose::Pose const &pose, Size res1, Size res2) const
single edge of the fold_tree
Edge const & get_residue_edge(int const seqpos) const
Returns the edge that builds the residue Does not work for root atom (will fail) ...
scoring::Energies const & energies() const
Returns the pose Energies (const-access)
virtual void setup_for_scoring(pose::Pose &pose, ScoreFunction const &) const
if an energy method needs to cache something in the pose (e.g. in pose.energies()), before scoring begins, it must do so in this method. All long range energy functions must initialize their LREnergyContainers before scoring begins. The default is to do nothing.