diff --git a/source/src/core.5.src.settings b/source/src/core.5.src.settings index 7df47feb3c7..3ddd3a0a611 100644 --- a/source/src/core.5.src.settings +++ b/source/src/core.5.src.settings @@ -26,6 +26,8 @@ sources = { "MinimizerOptions", "NelderMeadSimplex", "NumericalDerivCheckResult", + "ParametricAtomTreeMultifunc", + "parametric_minimize_util", "ParticleSwarmMinimizer", "SingleResidueMultifunc", ], diff --git a/source/src/core/conformation/parametric/SizeValuedParameter.cc b/source/src/core/conformation/parametric/SizeValuedParameter.cc index d1d05b5552c..3943456aa47 100644 --- a/source/src/core/conformation/parametric/SizeValuedParameter.cc +++ b/source/src/core/conformation/parametric/SizeValuedParameter.cc @@ -93,6 +93,9 @@ SizeValuedParameter::set_parameter_type( ) { runtime_assert( type_in == PT_generic_natural_number || type_in == PT_generic_whole_number ); Parameter::set_parameter_type( type_in ); + if ( type_in == PT_generic_natural_number && default_value_ == 0) { + default_value_ = 1; // can't be zero anymore. + } correct_range(); } @@ -160,7 +163,11 @@ SizeValuedParameter::provide_xsd_information( using namespace utility::tag; if ( provide_setting ) { - xsd + XMLSchemaAttribute::attribute_w_default( parameter_name(), parameter_type() == PT_generic_natural_number ? xsct_positive_integer : xsct_non_negative_integer, parameter_description(), std::to_string( default_value() ) ); + xsd + XMLSchemaAttribute::attribute_w_default( + parameter_name(), + parameter_type() == PT_generic_natural_number ? xsct_positive_integer : xsct_non_negative_integer, + parameter_description(), + std::to_string( default_value() ) ); } } diff --git a/source/src/core/kinematics/MoveMap.cc b/source/src/core/kinematics/MoveMap.cc index fa3ee1e8a2f..8d24d47b8d0 100644 --- a/source/src/core/kinematics/MoveMap.cc +++ b/source/src/core/kinematics/MoveMap.cc @@ -86,6 +86,9 @@ MoveMap::clear() dof_id_map_.clear(); jump_id_map_.clear(); + + parametric_ = false; + parametric_set_ = false; } /// set/get for JumpIDs --- fold-tree independent definition of jumps @@ -756,6 +759,8 @@ core::kinematics::MoveMap::save( Archive & arc ) const { arc( CEREAL_NVP( dof_id_map_ ) ); // DOF_ID_Map arc( CEREAL_NVP( jump_id_map_ ) ); // JumpID_Map arc( CEREAL_NVP( atom_id_map_) ); //AtomID_Map + arc( CEREAL_NVP( parametric_ ) ); + arc( CEREAL_NVP( parametric_set_ ) ); } /// @brief Automatically generated deserialization method @@ -769,6 +774,8 @@ core::kinematics::MoveMap::load( Archive & arc ) { arc( dof_id_map_ ); // DOF_ID_Map arc( jump_id_map_ ); // JumpID_Map arc( atom_id_map_ ); //AtomID_Map + arc( parametric_ ); + arc( parametric_set_ ); } SAVE_AND_LOAD_SERIALIZABLE( core::kinematics::MoveMap ); diff --git a/source/src/core/kinematics/MoveMap.hh b/source/src/core/kinematics/MoveMap.hh index 9e9d97a5ad5..45c67cbfa57 100644 --- a/source/src/core/kinematics/MoveMap.hh +++ b/source/src/core/kinematics/MoveMap.hh @@ -540,6 +540,18 @@ public: void set( DOF_ID const & id, bool const setting ); + /// @brief Set whether parametric DOFs (Crick parameters for helical bundles, barrels, etc.) + /// should be minimized. When enabled, residues under parametric control will have their + /// backbone torsional DOFs automatically excluded (to prevent redundant DOF control), + /// while their chi angles remain free if set_chi is also true. + inline void set_parametric( bool const setting ) { parametric_ = setting; parametric_set_ = true; } + + /// @brief Get whether parametric DOFs are enabled for minimization. + inline bool get_parametric() const { return parametric_set_ ? parametric_ : false; } + + /// @brief Was the parametric setting explicitly set? + inline bool parametric_is_set() const { return parametric_set_; } + public: // accessors /// @brief Returns if BB torsions are movable or not for residue @@ -854,6 +866,10 @@ private: /// Used for fine control of cartesian minimization/protocols std::map< AtomID, bool > atom_id_map_; + /// @brief Whether parametric DOFs (Crick parameters) should be included in minimization. + bool parametric_ = false; + bool parametric_set_ = false; + #ifdef SERIALIZATION public: template< class Archive > void save( Archive & arc ) const; diff --git a/source/src/core/optimization/ParametricAtomTreeMultifunc.cc b/source/src/core/optimization/ParametricAtomTreeMultifunc.cc new file mode 100644 index 00000000000..8b57e38c2e9 --- /dev/null +++ b/source/src/core/optimization/ParametricAtomTreeMultifunc.cc @@ -0,0 +1,234 @@ +// -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*- +// vi: set ts=2 noet: +// +// (c) Copyright Rosetta Commons Member Institutions. +// (c) This file is part of the Rosetta software suite and is made available under license. +// (c) The Rosetta software is developed by the contributing members of the Rosetta Commons. +// (c) For more information, see http://www.rosettacommons.org. Questions about this can be +// (c) addressed to University of Washington CoMotion, email: license@uw.edu. + +/// @file core/optimization/ParametricAtomTreeMultifunc.cc +/// @brief Multifunc for co-minimizing parametric and atom-tree DOFs. +/// @author Andy Watkins + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace core { +namespace optimization { + +static basic::Tracer TR( "core.optimization.ParametricAtomTreeMultifunc" ); + +ParametricAtomTreeMultifunc::ParametricAtomTreeMultifunc( + pose::Pose & pose_in, + MinimizerMap & min_map_in, + scoring::ScoreFunction const & scorefxn_in, + utility::vector1< ParametricDOFInfo > const & parametric_dofs_in, + RebuildCallback rebuild_callback, + bool const deriv_check_in, + bool const deriv_check_verbose_in +) : + pose_( pose_in ), + min_map_( min_map_in ), + score_function_( scorefxn_in ), + parametric_dofs_( parametric_dofs_in ), + n_standard_dofs_( min_map_in.nangles() ), + rebuild_callback_( std::move( rebuild_callback ) ), + deriv_check_( deriv_check_in ), + deriv_check_verbose_( deriv_check_verbose_in ) +{} + +ParametricAtomTreeMultifunc::~ParametricAtomTreeMultifunc() = default; + +void +ParametricAtomTreeMultifunc::apply_parametric_dofs( Multivec const & vars ) const { + for ( Size p = 1; p <= parametric_dofs_.size(); ++p ) { + Real const val = vars[ n_standard_dofs_ + p ]; + set_parametric_dof_value( pose_, parametric_dofs_[p], val ); + } + // Rebuild the full backbone (and sidechain frames) from current parameters. + // The callback is provided by protocols-level code and calls the appropriate + // ParametrizationCalculator::build_helix/build_strand, which correctly + // places all atoms — mainchain and sidechain — via the Crick equations + // and atom tree operations. + rebuild_callback_( pose_ ); +} + +Real +ParametricAtomTreeMultifunc::operator()( Multivec const & vars ) const { + PROF_START( basic::FUNC ); + + apply_parametric_dofs( vars ); + min_map_.copy_dofs_to_pose( pose_, vars ); + Real const score = score_function_( pose_ ); + + ++eval_count_; + if ( trajectory_stride_ > 0 && ( eval_count_ % trajectory_stride_ == 0 || eval_count_ == 1 ) ) { + std::string const fname = trajectory_prefix_ + "_" + std::to_string( eval_count_ ) + ".pdb"; + pose_.dump_pdb( fname ); + TR << "Trajectory frame " << eval_count_ << " score=" << score << " -> " << fname << std::endl; + } + + PROF_STOP( basic::FUNC ); + return score; +} + +void +ParametricAtomTreeMultifunc::dfunc( Multivec const & vars, Multivec & dE_dvars ) const { + PROF_START( basic::DFUNC ); + + apply_parametric_dofs( vars ); + + // Compute standard DOF derivatives via the normal atom_tree_dfunc pipeline. + // This fills dE_dvars[1..n_standard_dofs_] and populates + // min_map_.atom_derivatives_ with per-atom DerivVectorPairs. + atom_tree_dfunc( pose_, min_map_, score_function_, vars, dE_dvars ); + + // Extend dE_dvars to include parametric DOF slots + Size const ntotal = n_standard_dofs_ + parametric_dofs_.size(); + dE_dvars.resize( ntotal, 0.0 ); + + // Compute parametric DOF derivatives via the Crick Jacobian chain rule: + // dE/dParam = Σ_atoms (dE/dXYZ · dXYZ/dParam) + for ( Size p = 1; p <= parametric_dofs_.size(); ++p ) { + ParametricDOFInfo const & info = parametric_dofs_[p]; + Real dE_dparam = 0.0; + + // Get current Crick parameter values for Jacobian computation + conformation::parametric::ParametersSetCOP params_set = + pose_.conformation().parameters_set( info.params_set_index ); + conformation::parametric::ParametersCOP params = + params_set->parameters( info.params_index ); + + // Extract the Crick parameters needed for derivative computation. + // We need: r0, omega0, delta_omega0, omega1, z1, delta_omega1_all, plus + // per-atom r1, delta_omega1, delta_z1. + // These are stored as parameters in the Parameters object. + // The exact enum indices depend on the calculator type (BPC or BBPC). + // For generality, we look up by name. + Real r0_val = 0, omega0_val = 0, delta_omega0_val = 0; + Real omega1_val = 0, z1_val = 0, delta_omega1_all_val = 0; + utility::vector1< Real > r1_vals, delta_omega1_vals, delta_z1_vals; + //Size residues_per_repeat = 1; + + // Read parameter values by iterating and checking names + for ( Size i = 1; i <= params->num_parameters(); ++i ) { + conformation::parametric::ParameterCOP param = params->parameter_cop( i ); + std::string const & name = param->parameter_name(); + if ( name == "r0" ) { + r0_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value(); + } else if ( name == "omega0" ) { + omega0_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value(); + } else if ( name == "delta_omega0" ) { + delta_omega0_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value(); + } else if ( name == "omega1" ) { + omega1_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value(); + } else if ( name == "z1" ) { + z1_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value(); + } else if ( name == "delta_omega1" ) { + delta_omega1_all_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value(); + } else if ( name == "r1_peratom" ) { + r1_vals = utility::pointer::static_pointer_cast< conformation::parametric::RealVectorValuedParameter const >( param )->value(); + } else if ( name == "delta_omega1_peratom" ) { + delta_omega1_vals = utility::pointer::static_pointer_cast< conformation::parametric::RealVectorValuedParameter const >( param )->value(); + } else if ( name == "delta_z1_peratom" ) { + delta_z1_vals = utility::pointer::static_pointer_cast< conformation::parametric::RealVectorValuedParameter const >( param )->value(); + } /*else if ( name == "residues_per_repeat" ) { + residues_per_repeat = utility::pointer::static_pointer_cast< conformation::parametric::SizeValuedParameter const >( param )->value(); + }*/ + } + + if ( r1_vals.empty() || delta_omega1_vals.empty() || delta_z1_vals.empty() ) continue; + + // Compute omega1 relative to omega0 (as done in generate_atom_positions) + Real const omega1_relative = omega1_val - omega0_val; + + // Compute t values and iterate over mainchain atoms + Size const helix_length = info.helix_end_resid - info.helix_start_resid + 1; + Real t = -static_cast( helix_length + 2 ) / 2.0; + + Size atom_counter = 0; + for ( Size resid = info.helix_start_resid; resid <= info.helix_end_resid; ++resid ) { + Size const n_mainchain = pose_.residue( resid ).n_mainchain_atoms(); + for ( Size iatom = 1; iatom <= n_mainchain; ++iatom ) { + ++atom_counter; + Size const atom_index_in_repeat = ((atom_counter - 1) % r1_vals.size()) + 1; + + Real const r1 = r1_vals[ atom_index_in_repeat ]; + Real const dw1 = delta_omega1_vals[ atom_index_in_repeat ] + delta_omega1_all_val; + Real const dz1 = delta_z1_vals[ atom_index_in_repeat ]; + + // Compute Crick derivatives for this atom + bool deriv_failed = false; + numeric::crick_equations::CrickDerivatives crick_derivs = + numeric::crick_equations::compute_crick_derivatives( + t, r0_val, omega0_val, delta_omega0_val, + r1, omega1_relative, z1_val, + dw1, dz1, deriv_failed ); + + if ( deriv_failed ) continue; + + // Get dE/dXYZ for this atom from the MinimizerMap's atom_derivatives + Size const real_atomno = pose_.residue_type( resid ).mainchain_atom( iatom ); + scoring::DerivVectorPair const & dvp = min_map_.atom_derivatives( resid )[ real_atomno ]; + numeric::xyzVector< Real > const & dE_dXYZ = dvp.f2(); + + // Select the correct derivative component based on which parametric DOF this is + numeric::xyzVector< Real > dXYZ_dParam( 0.0, 0.0, 0.0 ); + std::string const & dof_name = info.param_name; + if ( dof_name == "r0" ) { + dXYZ_dParam = crick_derivs.dXYZ_dr0; + } else if ( dof_name == "omega0" ) { + dXYZ_dParam = crick_derivs.dXYZ_domega0; + } else if ( dof_name == "delta_omega0" ) { + dXYZ_dParam = crick_derivs.dXYZ_ddelta_omega0; + } else if ( dof_name == "delta_omega1" ) { + dXYZ_dParam = crick_derivs.dXYZ_ddelta_omega1; + } else if ( dof_name == "delta_t" || dof_name == "delta_z0" ) { + dXYZ_dParam = crick_derivs.dXYZ_ddelta_t; + } + + dE_dparam += dE_dXYZ.dot( dXYZ_dParam ); + } + // Advance t by 1 per residue (matching generate_atom_positions logic) + t += 1.0; + } + + dE_dvars[ n_standard_dofs_ + p ] = dE_dparam; + } + + if ( deriv_check_ ) { + numerical_derivative_check( min_map_, *this, vars, dE_dvars, nullptr, deriv_check_verbose_ ); + } + + PROF_STOP( basic::DFUNC ); +} + +void +ParametricAtomTreeMultifunc::dump( Multivec const & /*vars*/, Multivec const & /*vars2*/ ) const { + // Minimal implementation — extend if debugging is needed +} + +} // namespace optimization +} // namespace core diff --git a/source/src/core/optimization/ParametricAtomTreeMultifunc.hh b/source/src/core/optimization/ParametricAtomTreeMultifunc.hh new file mode 100644 index 00000000000..8cb4bc258d8 --- /dev/null +++ b/source/src/core/optimization/ParametricAtomTreeMultifunc.hh @@ -0,0 +1,99 @@ +// -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*- +// vi: set ts=2 noet: +// +// (c) Copyright Rosetta Commons Member Institutions. +// (c) This file is part of the Rosetta software suite and is made available under license. +// (c) The Rosetta software is developed by the contributing members of the Rosetta Commons. +// (c) For more information, see http://www.rosettacommons.org. Questions about this can be +// (c) addressed to University of Washington CoMotion, email: license@uw.edu. + +/// @file core/optimization/ParametricAtomTreeMultifunc.hh +/// @brief Multifunc that co-minimizes parametric DOFs (Crick parameters) alongside standard +/// atom-tree DOFs (chi angles, backbone torsions, jumps). +/// @details Parametric DOFs are appended to the end of the standard DOF vector. The operator() +/// rebuilds parametric backbone coordinates via the Crick equations before applying standard DOFs +/// and scoring. The dfunc() computes standard derivatives via atom_tree_dfunc, then adds +/// parametric derivatives via the analytical Crick Jacobian chain rule. +/// @author Andy Watkins + +#ifndef INCLUDED_core_optimization_ParametricAtomTreeMultifunc_hh +#define INCLUDED_core_optimization_ParametricAtomTreeMultifunc_hh + +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace core { +namespace optimization { + +class ParametricAtomTreeMultifunc : public Multifunc { + +public: + + /// @brief The rebuild callback receives a Pose and must rebuild all parametric + /// backbone coordinates (and sidechain frames) from the current parameter values. + /// This is provided by protocols-level code that knows about the specific + /// ParametrizationCalculator (e.g., BundleParametrizationCalculator::build_helix). + using RebuildCallback = std::function< void( pose::Pose & ) >; + + ParametricAtomTreeMultifunc( + pose::Pose & pose_in, + MinimizerMap & min_map_in, + scoring::ScoreFunction const & scorefxn_in, + utility::vector1< ParametricDOFInfo > const & parametric_dofs_in, + RebuildCallback rebuild_callback, + bool deriv_check_in = false, + bool deriv_check_verbose_in = false + ); + + ~ParametricAtomTreeMultifunc() override; + + Real operator()( Multivec const & vars ) const override; + + void dfunc( Multivec const & vars, Multivec & dE_dvars ) const override; + + void dump( Multivec const & vars, Multivec const & vars2 ) const override; + + Size n_standard_dofs() const { return n_standard_dofs_; } + Size n_parametric_dofs() const { return parametric_dofs_.size(); } + Size total_dofs() const { return n_standard_dofs_ + parametric_dofs_.size(); } + + /// @brief Enable trajectory dumping: write a PDB every stride evaluations. + void set_trajectory_dump( std::string const & prefix, Size stride = 1 ) { + trajectory_prefix_ = prefix; + trajectory_stride_ = stride; + } + +private: + + void apply_parametric_dofs( Multivec const & vars ) const; + + pose::Pose & pose_; + MinimizerMap & min_map_; + scoring::ScoreFunction const & score_function_; + + utility::vector1< ParametricDOFInfo > parametric_dofs_; + Size n_standard_dofs_; + + RebuildCallback rebuild_callback_; + + bool deriv_check_; + bool deriv_check_verbose_; + + mutable Size eval_count_ = 0; + std::string trajectory_prefix_; + Size trajectory_stride_ = 0; + +}; // ParametricAtomTreeMultifunc + +} // namespace optimization +} // namespace core + +#endif diff --git a/source/src/core/optimization/parametric_minimize_util.cc b/source/src/core/optimization/parametric_minimize_util.cc new file mode 100644 index 00000000000..506b3274446 --- /dev/null +++ b/source/src/core/optimization/parametric_minimize_util.cc @@ -0,0 +1,229 @@ +// -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*- +// vi: set ts=2 noet: +// +// (c) Copyright Rosetta Commons Member Institutions. +// (c) This file is part of the Rosetta software suite and is made available under license. +// (c) The Rosetta software is developed by the contributing members of the Rosetta Commons. +// (c) For more information, see http://www.rosettacommons.org. Questions about this can be +// (c) addressed to University of Washington CoMotion, email: license@uw.edu. + +/// @file core/optimization/parametric_minimize_util.cc +/// @brief Utility functions for minimizing over parametric DOFs (e.g. Crick parameters for helical bundles). +/// @author Andy Watkins (andy.watkins2@gmail.com) + +// Unit headers +#include + +// Package headers +#include +#include +#include +#include +#include +#include +#include + +// Basic headers +#include + +// Utility headers +#include + +// C++ headers +#include + +static basic::Tracer TR( "core.optimization.parametric_minimize_util" ); + +namespace core { +namespace optimization { + +/// @brief Enumerate all parametric DOFs from a pose that has ParametersSets. +/// @details For each ParametersSet, for each Parameters, for each perturbable +/// RealValuedParameter (those with can_be_perturbed()==true), creates a ParametricDOFInfo. +/// For parameters marked as global_for_parameters_set, only one DOF is created +/// (associated with the first Parameters object). +void +enumerate_parametric_dofs( + pose::Pose const & pose, + utility::vector1< ParametricDOFInfo > & dof_infos +) { + using namespace core::conformation::parametric; + + dof_infos.clear(); + + core::Size const n_param_sets( pose.conformation().n_parameters_sets() ); + if ( n_param_sets == 0 ) { + TR.Warning << "No ParametersSets found in the pose conformation. No parametric DOFs to enumerate." << std::endl; + return; + } + + for ( core::Size ps_index = 1; ps_index <= n_param_sets; ++ps_index ) { + + // Const access to the ParametersSet. + ParametersSetCOP params_set( pose.conformation().parameters_set( ps_index ) ); + runtime_assert_string_msg( params_set != nullptr, + "Error in core::optimization::enumerate_parametric_dofs(): Null ParametersSet encountered." ); + + core::Size const n_params( params_set->n_parameters() ); + + // Keep track of which global parameters we have already added (by parameter index within a Parameters object). + // We only want to add each global parameter once per ParametersSet. + std::set< core::Size > global_params_added; + + for ( core::Size p_index = 1; p_index <= n_params; ++p_index ) { + + ParametersCOP params( params_set->parameters( p_index ) ); + runtime_assert_string_msg( params != nullptr, + "Error in core::optimization::enumerate_parametric_dofs(): Null Parameters encountered." ); + + core::Size const n_individual_params( params->num_parameters() ); + + for ( core::Size ip = 1; ip <= n_individual_params; ++ip ) { + + ParameterCOP param( params->parameter_cop( ip ) ); + if ( param == nullptr ) continue; + + // Only consider RealValuedParameters. + RealValuedParameterCOP real_param( + utility::pointer::dynamic_pointer_cast< RealValuedParameter const >( param ) ); + if ( real_param == nullptr ) continue; + + // Only consider parameters that can be perturbed. + if ( !real_param->can_be_perturbed() ) continue; + + // For global parameters, only add once (associated with the first Parameters object). + if ( real_param->global_for_parameters_set() ) { + if ( global_params_added.count( ip ) > 0 ) continue; + global_params_added.insert( ip ); + } + + ParametricDOFInfo info; + info.params_set_index = ps_index; + info.params_index = p_index; + info.param_enum = ip; + info.helix_start_resid = params->first_residue_index(); + info.helix_end_resid = params->last_residue_index(); + info.is_global = real_param->global_for_parameters_set(); + info.scale_factor = 1.0; + info.param_name = real_param->parameter_name(); + + dof_infos.push_back( info ); + + TR.Debug << "Enumerated parametric DOF: params_set=" << ps_index + << " params=" << p_index + << " param=" << ip + << " name=" << real_param->parameter_name() + << " global=" << ( info.is_global ? "true" : "false" ) + << " residues=" << info.helix_start_resid << "-" << info.helix_end_resid + << std::endl; + } + } + } + + TR << "Enumerated " << dof_infos.size() << " parametric DOFs from " << n_param_sets << " ParametersSet(s)." << std::endl; +} + +/// @brief Get the set of residue indices that are under parametric control. +/// @details Reads the residue lists from all Parameters objects in all ParametersSets. +std::set< Size > +get_parametric_residues( + pose::Pose const & pose +) { + using namespace core::conformation::parametric; + + std::set< Size > parametric_resids; + + core::Size const n_param_sets( pose.conformation().n_parameters_sets() ); + for ( core::Size ps_index = 1; ps_index <= n_param_sets; ++ps_index ) { + + ParametersSetCOP params_set( pose.conformation().parameters_set( ps_index ) ); + if ( params_set == nullptr ) continue; + + core::Size const n_params( params_set->n_parameters() ); + for ( core::Size p_index = 1; p_index <= n_params; ++p_index ) { + + ParametersCOP params( params_set->parameters( p_index ) ); + if ( params == nullptr ) continue; + + core::Size const first_res( params->first_residue_index() ); + core::Size const last_res( params->last_residue_index() ); + for ( core::Size res = first_res; res <= last_res; ++res ) { + parametric_resids.insert( res ); + } + } + } + + return parametric_resids; +} + +/// @brief Get the current value of a parametric DOF from the pose. +Real +get_parametric_dof_value( + pose::Pose const & pose, + ParametricDOFInfo const & info +) { + using namespace core::conformation::parametric; + + runtime_assert_string_msg( info.params_set_index >= 1 && info.params_set_index <= pose.conformation().n_parameters_sets(), + "Error in core::optimization::get_parametric_dof_value(): params_set_index out of range." ); + + ParametersSetCOP params_set( pose.conformation().parameters_set( info.params_set_index ) ); + runtime_assert_string_msg( params_set != nullptr, + "Error in core::optimization::get_parametric_dof_value(): Null ParametersSet." ); + + runtime_assert_string_msg( info.params_index >= 1 && info.params_index <= params_set->n_parameters(), + "Error in core::optimization::get_parametric_dof_value(): params_index out of range." ); + + ParametersCOP params( params_set->parameters( info.params_index ) ); + runtime_assert_string_msg( params != nullptr, + "Error in core::optimization::get_parametric_dof_value(): Null Parameters." ); + + runtime_assert_string_msg( info.param_enum >= 1 && info.param_enum <= params->num_parameters(), + "Error in core::optimization::get_parametric_dof_value(): param_enum out of range." ); + + ParameterCOP param( params->parameter_cop( info.param_enum ) ); + RealValuedParameterCOP real_param( + utility::pointer::dynamic_pointer_cast< RealValuedParameter const >( param ) ); + runtime_assert_string_msg( real_param != nullptr, + "Error in core::optimization::get_parametric_dof_value(): Parameter is not a RealValuedParameter." ); + + return real_param->value(); +} + +/// @brief Set a parametric DOF value in the pose's ParametersSet. +void +set_parametric_dof_value( + pose::Pose & pose, + ParametricDOFInfo const & info, + Real value +) { + using namespace core::conformation::parametric; + + runtime_assert_string_msg( info.params_set_index >= 1 && info.params_set_index <= pose.conformation().n_parameters_sets(), + "Error in core::optimization::set_parametric_dof_value(): params_set_index out of range." ); + + ParametersSetOP params_set( pose.conformation().parameters_set( info.params_set_index ) ); + runtime_assert_string_msg( params_set != nullptr, + "Error in core::optimization::set_parametric_dof_value(): Null ParametersSet." ); + + runtime_assert_string_msg( info.params_index >= 1 && info.params_index <= params_set->n_parameters(), + "Error in core::optimization::set_parametric_dof_value(): params_index out of range." ); + + ParametersOP params( params_set->parameters( info.params_index ) ); + runtime_assert_string_msg( params != nullptr, + "Error in core::optimization::set_parametric_dof_value(): Null Parameters." ); + + runtime_assert_string_msg( info.param_enum >= 1 && info.param_enum <= params->num_parameters(), + "Error in core::optimization::set_parametric_dof_value(): param_enum out of range." ); + + ParameterOP param( params->parameter_op( info.param_enum ) ); + RealValuedParameterOP real_param( + utility::pointer::dynamic_pointer_cast< RealValuedParameter >( param ) ); + runtime_assert_string_msg( real_param != nullptr, + "Error in core::optimization::set_parametric_dof_value(): Parameter is not a RealValuedParameter." ); + + real_param->set_value( value, true /*ignore_use_degrees*/ ); +} + +} // namespace optimization +} // namespace core diff --git a/source/src/core/optimization/parametric_minimize_util.hh b/source/src/core/optimization/parametric_minimize_util.hh new file mode 100644 index 00000000000..3911d2e998b --- /dev/null +++ b/source/src/core/optimization/parametric_minimize_util.hh @@ -0,0 +1,86 @@ +// -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*- +// vi: set ts=2 noet: +// +// (c) Copyright Rosetta Commons Member Institutions. +// (c) This file is part of the Rosetta software suite and is made available under license. +// (c) The Rosetta software is developed by the contributing members of the Rosetta Commons. +// (c) For more information, see http://www.rosettacommons.org. Questions about this can be +// (c) addressed to University of Washington CoMotion, email: license@uw.edu. + +/// @file core/optimization/parametric_minimize_util.hh +/// @brief Utility functions for minimizing over parametric DOFs (e.g. Crick parameters for helical bundles). +/// @author Andy Watkins (andy.watkins2@gmail.com) + +#ifndef INCLUDED_core_optimization_parametric_minimize_util_hh +#define INCLUDED_core_optimization_parametric_minimize_util_hh + +// Project headers +#include +#include + +// Utility headers +#include + +// C++ headers +#include +#include + +namespace core { +namespace optimization { + +/// @brief Information about a single parametric DOF to be minimized. +struct ParametricDOFInfo { + Size params_set_index; // Index into pose.conformation().parameters_set() + Size params_index; // Index into ParametersSet::parameters() + Size param_enum; // Parameter enum value (e.g., BPC_r0, BBPC_r0) + Size helix_start_resid; // First residue of this helix/strand in the pose + Size helix_end_resid; // Last residue of this helix/strand in the pose + bool is_global; // True if this param is global for the ParametersSet + Real scale_factor; // Scale factor for DOF vector (default 1.0) + std::string param_name; // Parameter name (e.g., "r0", "omega0") for Jacobian dispatch + + ParametricDOFInfo() : + params_set_index( 0 ), + params_index( 0 ), + param_enum( 0 ), + helix_start_resid( 0 ), + helix_end_resid( 0 ), + is_global( false ), + scale_factor( 1.0 ), + param_name( "" ) + {} +}; + +/// @brief Enumerate all parametric DOFs from a pose that has ParametersSets. +/// @details For each ParametersSet, for each Parameters, for each perturbable +/// RealValuedParameter (those with can_be_perturbed()==true), creates a ParametricDOFInfo. +/// For parameters marked as global_for_parameters_set, only one DOF is created +/// (associated with the first Parameters object). +void enumerate_parametric_dofs( + pose::Pose const & pose, + utility::vector1< ParametricDOFInfo > & dof_infos +); + +/// @brief Get the set of residue indices that are under parametric control. +/// @details Reads the residue lists from all Parameters objects in all ParametersSets. +std::set< Size > get_parametric_residues( + pose::Pose const & pose +); + +/// @brief Get the current value of a parametric DOF from the pose. +Real get_parametric_dof_value( + pose::Pose const & pose, + ParametricDOFInfo const & info +); + +/// @brief Set a parametric DOF value in the pose's ParametersSet. +void set_parametric_dof_value( + pose::Pose & pose, + ParametricDOFInfo const & info, + Real value +); + +} // namespace optimization +} // namespace core + +#endif // INCLUDED_core_optimization_parametric_minimize_util_hh diff --git a/source/src/core/scoring/constraints/SequenceProfileConstraint.cc b/source/src/core/scoring/constraints/SequenceProfileConstraint.cc index f0e137d09a9..10c82fc8084 100644 --- a/source/src/core/scoring/constraints/SequenceProfileConstraint.cc +++ b/source/src/core/scoring/constraints/SequenceProfileConstraint.cc @@ -154,7 +154,8 @@ SequenceProfileConstraint::read_def( seqpos_ = residue_index; const_cast(sequence_profile_.get())->prof_row( aa_scores, residue_index ); } else { - auto residue_index(boost::lexical_cast(version)); + debug_assert(version >= 0); + auto residue_index(static_cast(version)); std::string profile_filename; is >> profile_filename; diff --git a/source/src/core/select/movemap/MoveMapFactory.cc b/source/src/core/select/movemap/MoveMapFactory.cc index 42eeae22e35..708523960e6 100644 --- a/source/src/core/select/movemap/MoveMapFactory.cc +++ b/source/src/core/select/movemap/MoveMapFactory.cc @@ -150,6 +150,8 @@ void MoveMapFactory::add_bondangles_action( void MoveMapFactory::all_bondlengths( bool setting ) { use_all_bondlengths_ = true; all_bondlengths_setting_ = setting; } +void MoveMapFactory::all_parametric( bool setting ) { use_all_parametric_ = true; all_parametric_setting_ = setting; } + void MoveMapFactory::add_bondlengths_action( move_map_action action, residue_selector::ResidueSelectorCOP selector @@ -203,6 +205,7 @@ MoveMapFactory::edit_movemap_given_pose( if ( use_all_jumps_ ) mm.set_jump( all_jumps_setting_ ); if ( use_all_bondangles_ ) mm.set( core::id::THETA, all_bondangles_setting_ ); if ( use_all_bondlengths_ ) mm.set( core::id::THETA, all_bondlengths_setting_ ); + if ( use_all_parametric_ ) mm.set_parametric( all_parametric_setting_ ); ///////// Backbone ///////// for ( auto const & bb_act : bb_actions_ ) { @@ -387,6 +390,9 @@ void MoveMapFactory::parse_my_tag( } set_cartesian(tag->getOption< bool >("cartesian", cartesian_specific_)); + if ( tag->hasOption( "parametric" ) ) { + all_parametric( tag->getOption< bool >( "parametric" ) ); + } for ( auto const & subtag : tag->getTags() ) { move_map_action enable = move_map_action( subtag->getOption< bool >( "enable", true ) ); @@ -485,7 +491,8 @@ void MoveMapFactory::provide_xml_schema( utility::tag::XMLSchemaDefinition & xsd + XMLSchemaAttribute( "nu", xsct_rosetta_bool, "Enable or disable movement for all nu torsions." ) + XMLSchemaAttribute( "branches", xsct_rosetta_bool, "Enable or disable movement for all branch torsions." ) + XMLSchemaAttribute( "jumps", xsct_rosetta_bool, "Enable or disable movement for all jump DOFs." ) - + XMLSchemaAttribute::attribute_w_default( "cartesian", xsct_rosetta_bool, "Set the MMF for specific cartesian overrides. Currently is only used for glycans in order to maintain IUPAC nomenclature during moves", "false" ); + + XMLSchemaAttribute::attribute_w_default( "cartesian", xsct_rosetta_bool, "Set the MMF for specific cartesian overrides. Currently is only used for glycans in order to maintain IUPAC nomenclature during moves", "false" ) + + XMLSchemaAttribute( "parametric", xsct_rosetta_bool, "Enable or disable minimization over parametric DOFs (Crick parameters for helical bundles, barrels, etc.)." ); XMLSchemaComplexTypeGenerator xsct; xsct.element_name( element_name() ) diff --git a/source/src/core/select/movemap/MoveMapFactory.hh b/source/src/core/select/movemap/MoveMapFactory.hh index 5614cd100e9..a01a4a3bfcc 100644 --- a/source/src/core/select/movemap/MoveMapFactory.hh +++ b/source/src/core/select/movemap/MoveMapFactory.hh @@ -88,6 +88,8 @@ public: void add_bondlengths_action( move_map_action action, residue_selector::ResidueSelectorCOP selector ); + void all_parametric( bool setting ); + ///@brief Option to use specific cartesian settings for the movemap. /// Used for glycans as dihedral vs cart MM settings are completely different to move the same /// underlying coordinates in IUPAC definitions @@ -182,9 +184,10 @@ private: bool all_bondlengths_setting_ = false; std::list< MMResAction > bondlengths_actions_; - bool cartesian_specific_ = false; //Option to use specific cartesian settings for the movemap. Used for glycans as dihedral vs cart MM settings are completely different to move the same underlying coordinates in IUPAC definitions - // TO DO: - // Specific torsion selectors! + bool cartesian_specific_ = false; + + bool use_all_parametric_ = false; + bool all_parametric_setting_ = false; }; diff --git a/source/src/numeric.src.settings b/source/src/numeric.src.settings index 958d3e447cc..dcda1956234 100644 --- a/source/src/numeric.src.settings +++ b/source/src/numeric.src.settings @@ -35,6 +35,7 @@ sources = { ], "numeric/crick_equations":[ "BundleParams", + "BundleParams_derivatives", "HelixParams", ], "numeric/expression_parser": [ diff --git a/source/src/numeric/crick_equations/BundleParams_derivatives.cc b/source/src/numeric/crick_equations/BundleParams_derivatives.cc new file mode 100644 index 00000000000..1cd98a3cef7 --- /dev/null +++ b/source/src/numeric/crick_equations/BundleParams_derivatives.cc @@ -0,0 +1,192 @@ +// -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*- +// vi: set ts=2 noet: +// +// (c) Copyright Rosetta Commons Member Institutions. +// (c) This file is part of the Rosetta software suite and is made available under license. +// (c) The Rosetta software is developed by the contributing members of the Rosetta Commons. +// (c) For more information, see http://www.rosettacommons.org. Questions about this can be +// (c) addressed to University of Washington CoMotion, email: license@uw.edu. + +/// @file numeric/crick_equations/BundleParams_derivatives.cc +/// @brief Analytical derivatives of the Crick equations with respect to bundle parameters. +/// @details For the epsilon=1 (circular cross-section) case, the Crick equations are: +/// +/// Pω0 = 2π · √(z1² - (r0·ω0)²) +/// α = atan2(2π·r0·ω0, Pω0) +/// gn = z1 (simplifies when epsilon=1) +/// +/// X = r0·c0 + r1·c1·c0 - r1·cos(α)·s0·s1 - (r0·ω0·s0/gn)·δz1 +/// Y = r0·s0 + r1·c1·s0 + r1·cos(α)·c0·s1 + (r0·ω0·c0/gn)·δz1 +/// Z = (Pω0·t/2π) - r1·sin(α)·s1 + (Pω0/2π/gn)·δz1 +/// +/// where c0=cos(ω0·t+δω0), s0=sin(ω0·t+δω0), c1=cos(ω1·t+δω1), s1=sin(ω1·t+δω1). +/// +/// @author Andy Watkins + +#include +#include + +#include + +namespace numeric { +namespace crick_equations { + +CrickDerivatives compute_crick_derivatives( + Real const t, + Real const r0, + Real const omega0, + Real const delta_omega0, + Real const r1, + Real const omega1, + Real const z1, + Real const delta_omega1, + Real const delta_z1, + bool & failed +) { + CrickDerivatives derivs; + failed = false; + + Real const twopi = numeric::constants::d::pi_2; + + // Intermediate quantities + Real const r0w0 = r0 * omega0; + Real const r0w0_sq = r0w0 * r0w0; + Real const z1_sq = z1 * z1; + Real const val = z1_sq - r0w0_sq; + if ( val < 0 ) { + failed = true; + return derivs; + } + Real const sqrt_val = std::sqrt( val ); + Real const Pw0 = twopi * sqrt_val; // P*omega0 + Real const Pw0_over_twopi = sqrt_val; // Pw0 / 2pi = sqrt(z1^2 - (r0*w0)^2) + Real const gn = z1; // gradnorm = z1 when epsilon=1 + + // alpha and its trig + Real const alpha = std::atan2( twopi * r0w0, Pw0 ); // = atan2(r0*w0, sqrt_val) + Real const cos_alpha = std::cos( alpha ); + Real const sin_alpha = std::sin( alpha ); + + // Trig basis + Real const c0 = std::cos( omega0 * t + delta_omega0 ); + Real const s0 = std::sin( omega0 * t + delta_omega0 ); + Real const c1 = std::cos( omega1 * t + delta_omega1 ); + Real const s1 = std::sin( omega1 * t + delta_omega1 ); + + // Ratio used in delta_z1 terms + Real const w0_over_gn = omega0 / gn; // omega0 / z1 + Real const r0w0_over_gn = r0w0 / gn; // r0*omega0 / z1 + Real const Pw0_over_twopi_gn = Pw0_over_twopi / gn; // sqrt(z1^2-(r0*w0)^2) / z1 + + // ======================== + // dXYZ / d(delta_omega0) + // ======================== + // delta_omega0 only appears in c0, s0: dc0/ddw0 = -s0, ds0/ddw0 = c0 + { + Real const dX = -r0 * s0 - r1 * c1 * s0 - r1 * cos_alpha * c0 * s1 - r0w0_over_gn * c0 * delta_z1; + Real const dY = r0 * c0 + r1 * c1 * c0 - r1 * cos_alpha * s0 * s1 - r0w0_over_gn * s0 * delta_z1; + Real const dZ = 0.0; + derivs.dXYZ_ddelta_omega0 = xyzVector( dX, dY, dZ ); + } + + // ======================== + // dXYZ / d(delta_omega1) + // ======================== + // delta_omega1 only appears in c1, s1: dc1/ddw1 = -s1, ds1/ddw1 = c1 + { + Real const dX = -r1 * s1 * c0 - r1 * cos_alpha * s0 * c1; + Real const dY = -r1 * s1 * s0 + r1 * cos_alpha * c0 * c1; + Real const dZ = -r1 * sin_alpha * c1; + derivs.dXYZ_ddelta_omega1 = xyzVector( dX, dY, dZ ); + } + + // ======================== + // dXYZ / d(delta_t) + // ======================== + // delta_t shifts t → t + delta_t. So dXYZ/d(delta_t) = dXYZ/dt. + // This is the tangent vector to the helix at position t. + { + // dc0/dt = -omega0*s0, ds0/dt = omega0*c0 + // dc1/dt = -omega1*s1, ds1/dt = omega1*c1 + Real const dX = -r0 * omega0 * s0 + + r1 * (-omega1 * s1 * c0 - omega0 * s0 * c1) + - r1 * cos_alpha * (omega0 * c0 * s1 + omega1 * s0 * c1) + - r0w0_over_gn * omega0 * c0 * delta_z1; + Real const dY = r0 * omega0 * c0 + + r1 * (-omega1 * s1 * s0 + omega0 * c0 * c1) + + r1 * cos_alpha * (-omega0 * s0 * s1 + omega1 * c0 * c1) + - r0w0_over_gn * omega0 * s0 * delta_z1; + Real const dZ = Pw0_over_twopi - r1 * sin_alpha * omega1 * c1; + derivs.dXYZ_ddelta_t = xyzVector( dX, dY, dZ ); + } + + // ======================== + // dXYZ / d(r0) + // ======================== + // r0 appears directly AND through Pw0 and alpha. + // dPw0/dr0 = 2pi * d/dr0[sqrt(z1^2-(r0*w0)^2)] = 2pi * (-r0*w0^2) / sqrt_val = -twopi*r0*w0^2/sqrt_val + // d(alpha)/dr0: alpha = atan2(twopi*r0*w0, Pw0) + // Let u = twopi*r0*w0, v = Pw0 = twopi*sqrt_val + // d(atan2(u,v))/dr0 = (v*du/dr0 - u*dv/dr0) / (u^2+v^2) + // du/dr0 = twopi*w0 + // dv/dr0 = twopi*(-r0*w0^2)/sqrt_val = dPw0/dr0 + // u^2+v^2 = (twopi*r0*w0)^2 + (twopi*sqrt_val)^2 = twopi^2 * z1^2 + // So d(alpha)/dr0 = [Pw0*twopi*w0 - twopi*r0*w0*(-twopi*r0*w0^2/sqrt_val)] / (twopi^2*z1^2) + // = twopi*w0*[Pw0 + twopi*r0^2*w0^2/sqrt_val] / (twopi^2*z1^2) + // = w0*[sqrt_val + r0^2*w0^2/sqrt_val] / (twopi*z1^2) + // = w0*z1^2/(sqrt_val*twopi*z1^2) + // = w0/(twopi*sqrt_val) + { + Real const dPw0_dr0 = -twopi * r0 * omega0 * omega0 / sqrt_val; + Real const dalpha_dr0 = omega0 / ( twopi * sqrt_val ); + Real const dcos_alpha_dr0 = -sin_alpha * dalpha_dr0; + Real const dsin_alpha_dr0 = cos_alpha * dalpha_dr0; + + // d(r0*c0)/dr0 = c0 + // d(r0*w0*s0/gn)/dr0 = w0*s0/gn (gn=z1 independent of r0) + Real const dX = c0 + r1 * 0.0 - r1 * dcos_alpha_dr0 * s0 * s1 - w0_over_gn * s0 * delta_z1; + Real const dY = s0 + r1 * 0.0 + r1 * dcos_alpha_dr0 * c0 * s1 + w0_over_gn * c0 * delta_z1; + Real const dZ = (dPw0_dr0 * t / twopi) - r1 * dsin_alpha_dr0 * s1 + (dPw0_dr0 / twopi / gn) * delta_z1; + derivs.dXYZ_dr0 = xyzVector( dX, dY, dZ ); + } + + // ======================== + // dXYZ / d(omega0) + // ======================== + // omega0 appears through c0, s0 (via t*omega0+delta_omega0), Pw0, alpha, and the delta_z1 terms. + // dc0/domega0 = -t*s0, ds0/domega0 = t*c0 + // dPw0/domega0 = 2pi * d/domega0[sqrt(z1^2-r0^2*w0^2)] = 2pi*(-r0^2*w0)/sqrt_val + // d(alpha)/domega0 = r0/(twopi*sqrt_val) [analogous derivation to dr0 case] + { + Real const dPw0_dw0 = -twopi * r0 * r0 * omega0 / sqrt_val; + Real const dalpha_dw0 = r0 / ( twopi * sqrt_val ); + Real const dcos_alpha_dw0 = -sin_alpha * dalpha_dw0; + Real const dsin_alpha_dw0 = cos_alpha * dalpha_dw0; + + // Direct trig derivatives from omega0 appearing in the argument of c0, s0: + Real const dc0 = -t * s0; + Real const ds0 = t * c0; + + // X = r0*c0 + r1*c1*c0 - r1*cos(alpha)*s0*s1 - (r0*w0*s0/gn)*dz1 + Real const dX = r0 * dc0 + + r1 * c1 * dc0 + - r1 * (dcos_alpha_dw0 * s0 + cos_alpha * ds0) * s1 + - (r0 * s0 / gn + r0w0_over_gn * ds0) * delta_z1; + + Real const dY = r0 * ds0 + + r1 * c1 * ds0 + + r1 * (dcos_alpha_dw0 * c0 + cos_alpha * dc0) * s1 + + (r0 * c0 / gn + r0w0_over_gn * dc0) * delta_z1; + + Real const dZ = (dPw0_dw0 * t / twopi) + - r1 * dsin_alpha_dw0 * s1 + + (dPw0_dw0 / twopi / gn) * delta_z1; + + derivs.dXYZ_domega0 = xyzVector( dX, dY, dZ ); + } + + return derivs; +} + +} // namespace crick_equations +} // namespace numeric diff --git a/source/src/numeric/crick_equations/BundleParams_derivatives.hh b/source/src/numeric/crick_equations/BundleParams_derivatives.hh new file mode 100644 index 00000000000..0d1de023685 --- /dev/null +++ b/source/src/numeric/crick_equations/BundleParams_derivatives.hh @@ -0,0 +1,63 @@ +// -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*- +// vi: set ts=2 noet: +// +// (c) Copyright Rosetta Commons Member Institutions. +// (c) This file is part of the Rosetta software suite and is made available under license. +// (c) The Rosetta software is developed by the contributing members of the Rosetta Commons. +// (c) For more information, see http://www.rosettacommons.org. Questions about this can be +// (c) addressed to University of Washington CoMotion, email: license@uw.edu. + +/// @file numeric/crick_equations/BundleParams_derivatives.hh +/// @brief Analytical derivatives of the Crick equations with respect to bundle parameters. +/// @details Computes the Jacobian dXYZ/dParam for each perturbable Crick parameter, +/// enabling gradient-based minimization in parametric space. +/// @author Andy Watkins + +#ifndef INCLUDED_numeric_crick_equations_BundleParams_derivatives_hh +#define INCLUDED_numeric_crick_equations_BundleParams_derivatives_hh + +#include +#include + +namespace numeric { +namespace crick_equations { + +/// @brief Derivatives of Crick XYZ coordinates with respect to each perturbable parameter. +struct CrickDerivatives { + xyzVector dXYZ_dr0; + xyzVector dXYZ_domega0; + xyzVector dXYZ_ddelta_omega0; + xyzVector dXYZ_ddelta_omega1; + xyzVector dXYZ_ddelta_t; +}; + +/// @brief Compute analytical derivatives of the Crick bundle equations. +/// @details Currently supports the epsilon=1 (circular cross-section) case only. +/// @param[in] t Position along the major helix (residue index, can be fractional) +/// @param[in] r0 Major helix radius (Angstroms) +/// @param[in] omega0 Major helix angular frequency (radians per residue) +/// @param[in] delta_omega0 Major helix phase offset (radians) +/// @param[in] r1 Minor helix radius for this atom (Angstroms) +/// @param[in] omega1 Minor helix angular frequency (radians per residue) +/// @param[in] z1 Minor helix rise per residue (Angstroms) +/// @param[in] delta_omega1 Minor helix phase offset for this atom (radians) +/// @param[in] delta_z1 Axial offset for this atom (Angstroms) +/// @param[out] failed Set to true if the derivative computation fails +/// @return CrickDerivatives struct with dXYZ/dParam for each parameter +CrickDerivatives compute_crick_derivatives( + Real t, + Real r0, + Real omega0, + Real delta_omega0, + Real r1, + Real omega1, + Real z1, + Real delta_omega1, + Real delta_z1, + bool & failed +); + +} // namespace crick_equations +} // namespace numeric + +#endif diff --git a/source/src/protocols/frag_picker/scores/AtomBasedConstraintsScore.cc b/source/src/protocols/frag_picker/scores/AtomBasedConstraintsScore.cc index 80b7d3f0a7d..3009d0ed65f 100644 --- a/source/src/protocols/frag_picker/scores/AtomBasedConstraintsScore.cc +++ b/source/src/protocols/frag_picker/scores/AtomBasedConstraintsScore.cc @@ -74,16 +74,6 @@ AtomBasedConstraintsScore::AtomBasedConstraintsScore(core::Size priority, constrainable_atoms_.insert(std::pair("2HB", 10)); } -std::string AtomBasedConstraintsScore::get_constrained_atom_name( core::Size atom_id ) { - for ( auto const & elem : constrainable_atoms_ ) { - if ( elem.second == atom_id ) { - return elem.first; - } - } - - return nullptr; // TODO this is a string return. why not just use std::find? -} - void AtomBasedConstraintsScore::do_caching(VallChunkOP chunk) { trAtomBasedConstraintsScore.Debug << "caching backbone atoms for " diff --git a/source/src/protocols/frag_picker/scores/AtomBasedConstraintsScore.hh b/source/src/protocols/frag_picker/scores/AtomBasedConstraintsScore.hh index 458bdb619ff..513a071de58 100644 --- a/source/src/protocols/frag_picker/scores/AtomBasedConstraintsScore.hh +++ b/source/src/protocols/frag_picker/scores/AtomBasedConstraintsScore.hh @@ -82,14 +82,10 @@ public: /// @brief returns an internal ID assigned to a given atom name /// @details this ID remains the same for all residues - inline core::Size get_constrained_atom_id(std::string atom_name) { + inline core::Size get_constrained_atom_id(std::string const & atom_name) { return constrainable_atoms_.find(atom_name)->second; } - - /// @brief returns a name of a constrained atom when its internal ID is known - /// @details this is the oposite to get_constrained_atom_id(std::string) - std::string get_constrained_atom_name(core::Size atom_id); - + /// @brief provides an access to the size of the length of a query sequence inline core::Size get_query_size() { return query_size_; diff --git a/source/src/protocols/helical_bundle/BundleParametrizationCalculator.cc b/source/src/protocols/helical_bundle/BundleParametrizationCalculator.cc index 0798ac0bcfc..15d3ac58cd0 100644 --- a/source/src/protocols/helical_bundle/BundleParametrizationCalculator.cc +++ b/source/src/protocols/helical_bundle/BundleParametrizationCalculator.cc @@ -478,7 +478,7 @@ BundleParametrizationCalculator::parameter_properties_from_enum( ParameterizationCalculatorProperties( true, true, true, true, false ), //z0_offset ParameterizationCalculatorProperties( true, true, true, true, false ), //z1_offset ParameterizationCalculatorProperties( true, true, true, true, false ), //epsilon - ParameterizationCalculatorProperties( false, false, false, false, false ), //residues_per_repeat + ParameterizationCalculatorProperties( true, false, false, false, false ), //residues_per_repeat ParameterizationCalculatorProperties( true, false, false, false, false ), //repeating_unit_offset ParameterizationCalculatorProperties( false, false, false, false, false ), //atoms_per_residue ParameterizationCalculatorProperties( true, false, false, false, false ), //r1_peratom diff --git a/source/src/protocols/minimization_packing/MinMover.cc b/source/src/protocols/minimization_packing/MinMover.cc index 97739185a3b..bfdec0e0793 100644 --- a/source/src/protocols/minimization_packing/MinMover.cc +++ b/source/src/protocols/minimization_packing/MinMover.cc @@ -36,6 +36,17 @@ #include #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include @@ -309,7 +320,80 @@ MinMover::minimize(pose::Pose & pose, core::kinematics::MoveMap const & active_m void MinMover::inner_run_minimizer( core::pose::Pose & pose, core::kinematics::MoveMap const & active_movemap ) { - if ( !cartesian( ) ) { + using namespace core::optimization; + + bool const use_parametric = active_movemap.get_parametric() && + pose.conformation().n_parameters_sets() > 0; + + if ( use_parametric && !cartesian() ) { + // Parametric minimization: co-optimize parametric DOFs alongside standard DOFs. + // This path is handled here (protocols layer) because the rebuild callback + // needs access to build_helix() which lives in protocols/. + + utility::vector1< ParametricDOFInfo > parametric_dofs; + enumerate_parametric_dofs( pose, parametric_dofs ); + std::set< core::Size > param_residues = get_parametric_residues( pose ); + + // Exclude backbone torsions of parametric residues to prevent redundant DOF control + core::kinematics::MoveMap effective_movemap( active_movemap ); + for ( core::Size resid : param_residues ) { + effective_movemap.set_bb( resid, false ); + } + + // Create the rebuild callback that properly rebuilds ALL atoms + // from the current parameter values by calling build_helix() + auto rebuild_callback = []( core::pose::Pose & p ) { + using namespace core::conformation::parametric; + for ( core::Size ps = 1; ps <= p.conformation().n_parameters_sets(); ++ps ) { + ParametersSetCOP params_set = p.conformation().parameters_set( ps ); + if ( params_set == nullptr ) continue; + for ( core::Size pidx = 1; pidx <= params_set->n_parameters(); ++pidx ) { + ParametersCOP params = params_set->parameters( pidx ); + if ( params == nullptr || params->n_residue() == 0 ) continue; + core::Size const start = params->first_residue_index(); + core::Size const end = params->last_residue_index(); + protocols::helical_bundle::BundleParametrizationCalculator calc( + false, utility::pointer::dynamic_pointer_cast< + protocols::helical_bundle::parameters::BundleParameters const >( params ) ); + calc.build_helix( p, start, end ); + } + } + }; + + score_before_minimization_ = (*scorefxn_)( pose ); + + MinimizerMap min_map; + min_map.setup( pose, effective_movemap ); + + bool const use_nblist = min_options_->use_nblist(); + if ( use_nblist ) { + pose.energies().set_use_nblist( pose, min_map.domain_map(), min_options_->nblist_auto_update() ); + } + scorefxn_->setup_for_minimizing( pose, min_map ); + + ParametricAtomTreeMultifunc f( pose, min_map, *scorefxn_, parametric_dofs, + rebuild_callback, min_options_->deriv_check(), min_options_->deriv_check_verbose() ); + + if ( TR.Debug.visible() ) { + f.set_trajectory_dump( "parametric_traj", 1 ); + } + + Multivec dofs( f.total_dofs() ); + min_map.copy_dofs_from_pose( pose, dofs ); + for ( core::Size p = 1; p <= parametric_dofs.size(); ++p ) { + dofs[ min_map.nangles() + p ] = get_parametric_dof_value( pose, parametric_dofs[p] ); + } + + Minimizer minimizer( f, *min_options_ ); + minimizer.run( dofs ); + f( dofs ); // apply final state + + if ( use_nblist ) pose.energies().reset_nblist(); + min_map.reset_jump_rb_deltas( pose, dofs ); + scorefxn_->finalize_after_minimizing( pose ); + score_after_minimization_ = (*scorefxn_)( pose ); + + } else if ( !cartesian() ) { AtomTreeMinimizerOP minimizer; if ( core::pose::symmetry::is_symmetric( pose ) ) { minimizer = utility::pointer::make_shared< core::optimization::symmetry::SymAtomTreeMinimizer >(); @@ -449,6 +533,9 @@ void MinMover::parse_movemap_factory( TagCOP const tag, basic::datacache::DataMa bool const value( tag->getOption("bondlength") ); mmf->all_bondlengths( value ); } + if ( tag->hasOption("parametric") ) { + mmf->all_parametric( tag->getOption("parametric") ); + } movemap_factory( protocols::rosetta_scripts::parse_movemap_factory_legacy( tag, data, false, mmf ) ); } @@ -559,7 +646,8 @@ MinMover::complex_type_generator_for_min_mover( utility::tag::XMLSchemaDefinitio attributes + XMLSchemaAttribute( "chi", xsct_rosetta_bool , "Minimize chi angles?" ) + XMLSchemaAttribute( "bb", xsct_rosetta_bool , "Minimize backbone torsion angles?" ) - + XMLSchemaAttribute::attribute_w_default( "omega", xsct_rosetta_bool , "Minimize omega torsions?", "true" ); + + XMLSchemaAttribute::attribute_w_default( "omega", xsct_rosetta_bool , "Minimize omega torsions?", "true" ) + + XMLSchemaAttribute( "parametric", xsct_rosetta_bool , "Minimize over parametric DOFs (Crick parameters for helical bundles, barrels, etc.)? Backbone torsions of parametric residues are automatically excluded." ); //All of these are lists of task operations, but none use parse_task_operations attributes + XMLSchemaAttribute( "bb_task_operations", xsct_task_operation_comma_separated_list , "Task operations specifying residues for backbone minimization" ) diff --git a/source/src/protocols/rosetta_scripts/util.cc b/source/src/protocols/rosetta_scripts/util.cc index 4605a3363c0..4f8b3ecb518 100644 --- a/source/src/protocols/rosetta_scripts/util.cc +++ b/source/src/protocols/rosetta_scripts/util.cc @@ -335,6 +335,7 @@ parse_movemap_tag( if ( in_tag->hasOption( "bb" ) ) mmf->all_bb( in_tag->getOption< bool >( "bb" ) ); if ( in_tag->hasOption( "chi" ) ) mmf->all_chi( in_tag->getOption< bool >( "chi" ) ); if ( in_tag->hasOption( "jump" ) ) mmf->all_jumps( in_tag->getOption< bool >( "jump" ) ); + if ( in_tag->hasOption( "parametric" ) ) mmf->all_parametric( in_tag->getOption< bool >( "parametric" ) ); } @@ -475,7 +476,8 @@ common_movemap_complex_type_def( utility::tag::XMLSchemaComplexTypeGenerator & c movemap_tag_attributes + XMLSchemaAttribute( "bb", xsct_rosetta_bool , bb_desc ) + XMLSchemaAttribute( "chi", xsct_rosetta_bool , chi_desc ) - + XMLSchemaAttribute( "jump", xsct_rosetta_bool , "move all jumps?" ); + + XMLSchemaAttribute( "jump", xsct_rosetta_bool , "move all jumps?" ) + + XMLSchemaAttribute( "parametric", xsct_rosetta_bool , "Minimize over parametric DOFs (Crick parameters for helical bundles, barrels, etc.)?" ); ct_gen.element_name( "MoveMap" ) .description( "MoveMaps dicate what degrees of freedom are mobile to other bits of code, especially minimization - they do NOT work with packing" ) diff --git a/source/src/protocols/sparta/PDB.cc b/source/src/protocols/sparta/PDB.cc index 35e62276ba2..09bc856d99d 100644 --- a/source/src/protocols/sparta/PDB.cc +++ b/source/src/protocols/sparta/PDB.cc @@ -74,7 +74,7 @@ string PDB::getField(const string &str, int index) case 11 : return str.substr(60,6); // B-factor } - return nullptr; + return "You've provided an illegal argument and this function is only permitted to return a string."; } diff --git a/source/test/numeric.test.settings b/source/test/numeric.test.settings index a13fc1b68a8..4f10ac1a474 100644 --- a/source/test/numeric.test.settings +++ b/source/test/numeric.test.settings @@ -42,6 +42,9 @@ sources = { "xyzTriple", "xyzVector", ], + "crick_equations" : [ + "BundleParams_derivatives", + ], "fourier" : [ "fft", ], diff --git a/source/test/numeric/crick_equations/BundleParams_derivatives.cxxtest.hh b/source/test/numeric/crick_equations/BundleParams_derivatives.cxxtest.hh new file mode 100644 index 00000000000..cae27a61b36 --- /dev/null +++ b/source/test/numeric/crick_equations/BundleParams_derivatives.cxxtest.hh @@ -0,0 +1,149 @@ +// -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*- +// vi: set ts=2 noet: +// +// (c) Copyright Rosetta Commons Member Institutions. +// (c) This file is part of the Rosetta software suite and is made available under license. +// (c) The Rosetta software is developed by the contributing members of the Rosetta Commons. +// (c) For more information, see http://www.rosettacommons.org. Questions about this can be +// (c) addressed to University of Washington CoMotion, email: license@uw.edu. + +/// @file numeric/crick_equations/BundleParams_derivatives.cxxtest.hh +/// @brief Unit tests for analytical Crick equation derivatives. +/// @author Andy Watkins + +#include +#include + +#include +#include + +#include + +static basic::Tracer TR("CrickDerivativesTests"); + +class CrickDerivativesTests : public CxxTest::TestSuite { + +public: + + void setUp() { + core_init(); + } + + void tearDown() {} + + numeric::xyzVector numerical_XYZ( + numeric::Real t, numeric::Real r0, numeric::Real omega0, numeric::Real delta_omega0, + numeric::Real r1, numeric::Real omega1, numeric::Real z1, + numeric::Real delta_omega1, numeric::Real delta_z1 + ) { + bool failed = false; + return numeric::crick_equations::XYZ_BUNDLE( + t, r0, omega0, delta_omega0, r1, omega1, z1, + delta_omega1, delta_z1, 1.0 /*epsilon*/, failed ); + } + + void check_derivative( + std::string const & param_name, + numeric::Real t, numeric::Real r0, numeric::Real omega0, numeric::Real delta_omega0, + numeric::Real r1, numeric::Real omega1, numeric::Real z1, + numeric::Real delta_omega1, numeric::Real delta_z1 + ) { + using numeric::Real; + + bool failed = false; + numeric::crick_equations::CrickDerivatives derivs = + numeric::crick_equations::compute_crick_derivatives( + t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1, failed ); + TS_ASSERT( !failed ); + if ( failed ) return; + + Real const eps = 1e-6; + numeric::xyzVector analytical_deriv; + numeric::xyzVector fwd, bwd; + + if ( param_name == "r0" ) { + analytical_deriv = derivs.dXYZ_dr0; + fwd = numerical_XYZ(t, r0+eps, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + bwd = numerical_XYZ(t, r0-eps, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + } else if ( param_name == "omega0" ) { + analytical_deriv = derivs.dXYZ_domega0; + fwd = numerical_XYZ(t, r0, omega0+eps, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + bwd = numerical_XYZ(t, r0, omega0-eps, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + } else if ( param_name == "delta_omega0" ) { + analytical_deriv = derivs.dXYZ_ddelta_omega0; + fwd = numerical_XYZ(t, r0, omega0, delta_omega0+eps, r1, omega1, z1, delta_omega1, delta_z1); + bwd = numerical_XYZ(t, r0, omega0, delta_omega0-eps, r1, omega1, z1, delta_omega1, delta_z1); + } else if ( param_name == "delta_omega1" ) { + analytical_deriv = derivs.dXYZ_ddelta_omega1; + fwd = numerical_XYZ(t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1+eps, delta_z1); + bwd = numerical_XYZ(t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1-eps, delta_z1); + } else if ( param_name == "delta_t" ) { + analytical_deriv = derivs.dXYZ_ddelta_t; + fwd = numerical_XYZ(t+eps, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + bwd = numerical_XYZ(t-eps, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + } else { + TS_FAIL("Unknown parameter name: " + param_name); + return; + } + + numeric::xyzVector numerical_deriv = (fwd - bwd) / (2.0 * eps); + + Real const tol = 1e-4; + TS_ASSERT_DELTA( analytical_deriv.x(), numerical_deriv.x(), tol ); + TS_ASSERT_DELTA( analytical_deriv.y(), numerical_deriv.y(), tol ); + TS_ASSERT_DELTA( analytical_deriv.z(), numerical_deriv.z(), tol ); + + if ( std::abs(analytical_deriv.x() - numerical_deriv.x()) > tol || + std::abs(analytical_deriv.y() - numerical_deriv.y()) > tol || + std::abs(analytical_deriv.z() - numerical_deriv.z()) > tol ) { + TR << "MISMATCH for " << param_name << " at t=" << t + << " r0=" << r0 << " omega0=" << omega0 << std::endl; + TR << " analytical: " << analytical_deriv.x() << " " << analytical_deriv.y() << " " << analytical_deriv.z() << std::endl; + TR << " numerical: " << numerical_deriv.x() << " " << numerical_deriv.y() << " " << numerical_deriv.z() << std::endl; + } + } + + void test_derivatives_alpha_helix_params() { + TR << "Testing Crick derivatives against numerical for alpha helix parameters." << std::endl; + + // Alpha helix-like parameters + numeric::Real const r0 = 5.0; + numeric::Real const omega0 = -0.045; // radians per residue + numeric::Real const delta_omega0 = 0.5; + numeric::Real const r1 = 2.26; + numeric::Real const omega1 = 1.72; + numeric::Real const z1 = 1.51; + numeric::Real const delta_omega1 = 0.0; + numeric::Real const delta_z1 = 0.0; + + for ( numeric::Real t = -5.0; t <= 5.0; t += 2.5 ) { + check_derivative("r0", t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + check_derivative("omega0", t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + check_derivative("delta_omega0", t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + check_derivative("delta_omega1", t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + check_derivative("delta_t", t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + } + } + + void test_derivatives_with_nonzero_delta_z1() { + TR << "Testing Crick derivatives with nonzero delta_z1." << std::endl; + + numeric::Real const r0 = 6.5; + numeric::Real const omega0 = -0.045; + numeric::Real const delta_omega0 = 1.2; + numeric::Real const r1 = 1.5; + numeric::Real const omega1 = 1.72; + numeric::Real const z1 = 1.51; + numeric::Real const delta_omega1 = 0.48; + numeric::Real const delta_z1 = 1.05; + + for ( numeric::Real t = -3.0; t <= 3.0; t += 1.5 ) { + check_derivative("r0", t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + check_derivative("omega0", t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + check_derivative("delta_omega0", t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + check_derivative("delta_omega1", t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + check_derivative("delta_t", t, r0, omega0, delta_omega0, r1, omega1, z1, delta_omega1, delta_z1); + } + } + +}; diff --git a/source/tools/build/basic.settings b/source/tools/build/basic.settings index 29f1e3d6d01..9e93d9068f6 100644 --- a/source/tools/build/basic.settings +++ b/source/tools/build/basic.settings @@ -2209,6 +2209,23 @@ settings = { }, }, + "clang, |cxx_ver:>=21|" : { + "appends" : { + "flags" : { + "warn" : [ + "Wno-error=deprecated-declarations", + ], + }, + }, + "removes" : { + "flags" : { + "warn" : [ + "Wno-error=enum-constexpr-conversion", + ], + }, + }, + }, + # OSs ##################################################################### "clang, linux" : { diff --git a/tests/integration/tests/parametric_minimize/README.txt b/tests/integration/tests/parametric_minimize/README.txt new file mode 100644 index 00000000000..9bc3e0563a3 --- /dev/null +++ b/tests/integration/tests/parametric_minimize/README.txt @@ -0,0 +1,16 @@ +Parametric Minimization Integration Test + +Tests simultaneous minimization of Crick parametric DOFs (r0, omega0, +delta_omega0, etc.) alongside standard atom-tree DOFs (chi angles). + +Creates a 3-helix bundle using MakeBundle with a realistic ELLKAIA heptad +repeat sequence (21 residues per helix, 63 total), then minimizes using +MinMover with the new parametric DOF support. + +Two runs are performed: + 1. Parametric DOFs only (chi angles frozen) + 2. Parametric DOFs + chi angles simultaneously + +Both runs should converge to lower energy than the starting structure. +The parametric+chi run should achieve at least as low energy as the +parametric-only run. diff --git a/tests/integration/tests/parametric_minimize/command b/tests/integration/tests/parametric_minimize/command new file mode 100644 index 00000000000..799d34b3376 --- /dev/null +++ b/tests/integration/tests/parametric_minimize/command @@ -0,0 +1,34 @@ +# +# Integration test for parametric minimization. +# Builds a 3-helix bundle with ELLKAIA heptad repeat and minimizes: +# Run 1: Parametric DOFs only (no chi angles) +# Run 2: Parametric DOFs + chi angles simultaneously +# +# Enable trajectory PDB dumps via Debug tracer level. +# + +cd %(workdir)s + +[ -x %(bin)s/rosetta_scripts.%(binext)s ] || exit 1 + +# Run 1: parametric-only minimization (with trajectory dump) +%(bin)s/rosetta_scripts.%(binext)s %(additional_flags)s @flags \ + -parser:protocol inputs/param_only.xml \ + -out:prefix param_only_ \ + -out:levels core.optimization.AtomTreeMinimizer:debug \ + -database %(database)s -testing:INTEGRATION_TEST 2>&1 \ + | egrep -vf ../../ignore_list \ + > log_param_only + +test "${PIPESTATUS[0]}" != '0' && exit 1 || true + +# Run 2: parametric + chi minimization (with trajectory dump) +%(bin)s/rosetta_scripts.%(binext)s %(additional_flags)s @flags \ + -parser:protocol inputs/param_plus_chi.xml \ + -out:prefix param_plus_chi_ \ + -out:levels core.optimization.AtomTreeMinimizer:debug \ + -database %(database)s -testing:INTEGRATION_TEST 2>&1 \ + | egrep -vf ../../ignore_list \ + > log_param_plus_chi + +test "${PIPESTATUS[0]}" != '0' && exit 1 || true diff --git a/tests/integration/tests/parametric_minimize/flags b/tests/integration/tests/parametric_minimize/flags new file mode 100644 index 00000000000..1485df056f7 --- /dev/null +++ b/tests/integration/tests/parametric_minimize/flags @@ -0,0 +1,4 @@ +-nstruct 1 +-overwrite +-out:file:scorefile score.sc +-input_empty_pose true diff --git a/tests/integration/tests/parametric_minimize/inputs/param_only.xml b/tests/integration/tests/parametric_minimize/inputs/param_only.xml new file mode 100644 index 00000000000..0ba81054bf4 --- /dev/null +++ b/tests/integration/tests/parametric_minimize/inputs/param_only.xml @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml b/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml new file mode 100644 index 00000000000..d4e33dbe370 --- /dev/null +++ b/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +