From f84a8d234b4b9854783474b432ca98d0cba1aa82 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Fri, 1 May 2026 16:10:25 +0000 Subject: [PATCH 01/13] Add differentiable parametric backbone minimization Enables MinMover to co-optimize parametric DOFs (r0, omega0, delta_omega0, etc.) alongside standard atom-tree DOFs (chi angles, backbone torsions for non-parametric residues, jumps). Components: - Analytical Crick equation derivatives (dXYZ/dParam) in numeric/crick_equations/BundleParams_derivatives for r0, omega0, delta_omega0, delta_omega1, and delta_t - MoveMap::set_parametric(true) flag to enable parametric DOFs - ParametricAtomTreeMultifunc: extends the DOF vector with parametric variables, computes gradients via Jacobian chain rule (dE/dParam = sum_atoms dE/dXYZ . dXYZ/dParam) - AtomTreeMinimizer dispatch: detects parametric poses and automatically excludes backbone torsions of parametric residues to prevent redundant DOF control - Integration test: 3-helix ELLKAIA bundle minimized with and without chi angle co-optimization --- source/src/core.5.src.settings | 2 + source/src/core/kinematics/MoveMap.cc | 7 + source/src/core/kinematics/MoveMap.hh | 16 + .../core/optimization/AtomTreeMinimizer.cc | 76 +++- .../ParametricAtomTreeMultifunc.cc | 220 ++++++++++++ .../ParametricAtomTreeMultifunc.hh | 79 +++++ .../optimization/parametric_minimize_util.cc | 325 ++++++++++++++++++ .../optimization/parametric_minimize_util.hh | 94 +++++ source/src/numeric.src.settings | 1 + .../BundleParams_derivatives.cc | 192 +++++++++++ .../BundleParams_derivatives.hh | 63 ++++ source/test/numeric.test.settings | 3 + .../BundleParams_derivatives.cxxtest.hh | 149 ++++++++ .../tests/parametric_minimize/README.txt | 16 + .../tests/parametric_minimize/command | 30 ++ .../tests/parametric_minimize/flags | 3 + .../parametric_minimize/inputs/param_only.xml | 31 ++ .../inputs/param_plus_chi.xml | 31 ++ 18 files changed, 1321 insertions(+), 17 deletions(-) create mode 100644 source/src/core/optimization/ParametricAtomTreeMultifunc.cc create mode 100644 source/src/core/optimization/ParametricAtomTreeMultifunc.hh create mode 100644 source/src/core/optimization/parametric_minimize_util.cc create mode 100644 source/src/core/optimization/parametric_minimize_util.hh create mode 100644 source/src/numeric/crick_equations/BundleParams_derivatives.cc create mode 100644 source/src/numeric/crick_equations/BundleParams_derivatives.hh create mode 100644 source/test/numeric/crick_equations/BundleParams_derivatives.cxxtest.hh create mode 100644 tests/integration/tests/parametric_minimize/README.txt create mode 100644 tests/integration/tests/parametric_minimize/command create mode 100644 tests/integration/tests/parametric_minimize/flags create mode 100644 tests/integration/tests/parametric_minimize/inputs/param_only.xml create mode 100644 tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml 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/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/AtomTreeMinimizer.cc b/source/src/core/optimization/AtomTreeMinimizer.cc index 245904ff794..62995a06993 100644 --- a/source/src/core/optimization/AtomTreeMinimizer.cc +++ b/source/src/core/optimization/AtomTreeMinimizer.cc @@ -22,9 +22,12 @@ #include #include #include +#include +#include // Project headers #include +#include #include #include #include @@ -79,9 +82,31 @@ AtomTreeMinimizer::run( pose.energies().show( TR.Trace ); } - // setup the map of the degrees of freedom + // Check whether parametric DOFs should be included + bool const use_parametric = move_map.get_parametric() && pose.conformation().n_parameters_sets() > 0; + utility::vector1< ParametricDOFInfo > parametric_dofs; + std::set< Size > parametric_residues; + if ( use_parametric ) { + enumerate_parametric_dofs( pose, parametric_dofs ); + parametric_residues = get_parametric_residues( pose ); + if ( TR.Debug.visible() ) { + TR.Debug << "Parametric minimization: " << parametric_dofs.size() << " parametric DOFs, " + << parametric_residues.size() << " residues under parametric control." << std::endl; + } + } + + // Setup the map of the degrees of freedom. + // If parametric minimization is active, we need a MoveMap that excludes backbone + // torsions of parametric residues (to prevent redundant DOF control). + kinematics::MoveMap effective_move_map( move_map ); + if ( use_parametric ) { + for ( Size resid : parametric_residues ) { + effective_move_map.set_bb( resid, false ); + } + } + MinimizerMap min_map; - min_map.setup( pose, move_map ); + min_map.setup( pose, effective_move_map ); // if we are using the nblist, set it up if ( use_nblist ) { @@ -91,27 +116,44 @@ AtomTreeMinimizer::run( scorefxn.setup_for_minimizing( pose, min_map ); - // setup the function that we will pass to the low-level minimizer - AtomTreeMultifunc f( pose, min_map, scorefxn, - options.deriv_check(), options.deriv_check_verbose() ); + // Setup the multifunc — use parametric version if parametric DOFs are present + Real start_func, end_func; + if ( use_parametric && !parametric_dofs.empty() ) { + ParametricAtomTreeMultifunc f( pose, min_map, scorefxn, parametric_dofs, + options.deriv_check(), options.deriv_check_verbose() ); - if ( deriv_check_result_ ) f.set_deriv_check_result( deriv_check_result_ ); + Multivec dofs( f.total_dofs() ); + min_map.copy_dofs_from_pose( pose, dofs ); + for ( Size p = 1; p <= parametric_dofs.size(); ++p ) { + dofs[ min_map.nangles() + p ] = get_parametric_dof_value( pose, parametric_dofs[p] ); + } - // starting position -- "dofs" = Degrees Of Freedom - Multivec dofs( min_map.nangles() ); - min_map.copy_dofs_from_pose( pose, dofs ); + start_func = f( dofs ); + if ( TR.Trace.visible() && !use_nblist ) { + pose.energies().show( TR.Trace ); + } - Real const start_func( f( dofs ) ); + Minimizer minimizer( f, options ); + minimizer.run( dofs ); + end_func = f( dofs ); + } else { + AtomTreeMultifunc f( pose, min_map, scorefxn, + options.deriv_check(), options.deriv_check_verbose() ); - if ( TR.Trace.visible() && ! use_nblist ) { - pose.energies().show( TR.Trace ); - } + if ( deriv_check_result_ ) f.set_deriv_check_result( deriv_check_result_ ); - // now do the optimization with the low-level minimizer function - Minimizer minimizer( f, options ); - minimizer.run( dofs ); + Multivec dofs( min_map.nangles() ); + min_map.copy_dofs_from_pose( pose, dofs ); - Real const end_func( f( dofs ) ); + start_func = f( dofs ); + if ( TR.Trace.visible() && !use_nblist ) { + pose.energies().show( TR.Trace ); + } + + Minimizer minimizer( f, options ); + minimizer.run( dofs ); + end_func = f( dofs ); + } if ( TR.Debug.visible() && ! use_nblist ) { TR.Debug << "end_func: " << end_func << std::endl; diff --git a/source/src/core/optimization/ParametricAtomTreeMultifunc.cc b/source/src/core/optimization/ParametricAtomTreeMultifunc.cc new file mode 100644 index 00000000000..923c7884e28 --- /dev/null +++ b/source/src/core/optimization/ParametricAtomTreeMultifunc.cc @@ -0,0 +1,220 @@ +// -*- 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, + 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() ), + 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_parametric_backbone( 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_ ); + + 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..a8f2b4b0917 --- /dev/null +++ b/source/src/core/optimization/ParametricAtomTreeMultifunc.hh @@ -0,0 +1,79 @@ +// -*- 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 + +namespace core { +namespace optimization { + +class ParametricAtomTreeMultifunc : public Multifunc { + +public: + + ParametricAtomTreeMultifunc( + pose::Pose & pose_in, + MinimizerMap & min_map_in, + scoring::ScoreFunction const & scorefxn_in, + utility::vector1< ParametricDOFInfo > const & parametric_dofs_in, + 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(); } + +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_; + + bool deriv_check_; + bool deriv_check_verbose_; + +}; // 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..efa64a1b8aa --- /dev/null +++ b/source/src/core/optimization/parametric_minimize_util.cc @@ -0,0 +1,325 @@ +// -*- 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 +#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*/ ); +} + +void +rebuild_parametric_backbone( + pose::Pose & pose +) { + using namespace core::conformation::parametric; + + for ( core::Size ps = 1; ps <= pose.conformation().n_parameters_sets(); ++ps ) { + ParametersSetCOP params_set = pose.conformation().parameters_set( ps ); + if ( params_set == nullptr ) continue; + for ( core::Size p = 1; p <= params_set->n_parameters(); ++p ) { + ParametersCOP params = params_set->parameters( p ); + if ( params == nullptr ) continue; + if ( params->n_residue() == 0 ) continue; + + core::Size const start = params->first_residue_index(); + core::Size const end = params->last_residue_index(); + + // Rebuild this parametric element's backbone coordinates. + // We use the helical_bundle utility functions directly: generate new atom + // positions from current parameter values, then place them. + // For now, this is a simplified rebuild that re-applies the Crick equations. + // A more complete version would use the ParametrizationCalculator's build_helix/build_strand. + + // Extract parameters by name for the Crick equation call + Real r0_val = 0, omega0_val = 0, delta_omega0_val = 0; + Real omega1_val = 0, z1_val = 0, delta_omega1_all_val = 0; + Real delta_t_val = 0, epsilon_val = 1.0; + bool invert = false; + utility::vector1< Real > r1_vals, delta_omega1_vals, delta_z1_vals; + core::Size residues_per_repeat = 1, repeating_unit_offset = 0; + utility::vector1< core::Size > atoms_per_residue_vals; + + for ( core::Size i = 1; i <= params->num_parameters(); ++i ) { + ParameterCOP param = params->parameter_cop( i ); + std::string const & name = param->parameter_name(); + if ( name == "r0" ) { + r0_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); + } else if ( name == "omega0" ) { + omega0_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); + } else if ( name == "delta_omega0" ) { + delta_omega0_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); + } else if ( name == "omega1" ) { + omega1_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); + } else if ( name == "z1" ) { + z1_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); + } else if ( name == "delta_omega1" ) { + delta_omega1_all_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); + } else if ( name == "delta_t" || name == "delta_z0" ) { + delta_t_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); + } else if ( name == "epsilon" ) { + epsilon_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); + } else if ( name == "invert" ) { + invert = utility::pointer::static_pointer_cast< BooleanValuedParameter const >( param )->value(); + } else if ( name == "r1_peratom" ) { + r1_vals = utility::pointer::static_pointer_cast< RealVectorValuedParameter const >( param )->value(); + } else if ( name == "delta_omega1_peratom" ) { + delta_omega1_vals = utility::pointer::static_pointer_cast< RealVectorValuedParameter const >( param )->value(); + } else if ( name == "delta_z1_peratom" ) { + delta_z1_vals = utility::pointer::static_pointer_cast< RealVectorValuedParameter const >( param )->value(); + } else if ( name == "residues_per_repeat" ) { + residues_per_repeat = utility::pointer::static_pointer_cast< SizeValuedParameter const >( param )->value(); + } else if ( name == "atoms_per_residue" ) { + atoms_per_residue_vals = utility::pointer::static_pointer_cast< SizeVectorValuedParameter const >( param )->value(); + } + } + + if ( r1_vals.empty() ) continue; + + // Use helical_bundle utilities for the actual rebuild + utility::vector1< utility::vector1< numeric::xyzVector< Real > > > atom_positions; + bool failed = false; + protocols::helical_bundle::generate_atom_positions( + atom_positions, pose, start, end, + r0_val, omega0_val, delta_omega0_val, delta_t_val, + 0.0 /*z1_offset*/, 0.0 /*z0_offset*/, epsilon_val, + invert, r1_vals, omega1_val, z1_val, + delta_omega1_vals, delta_omega1_all_val, + delta_z1_vals, residues_per_repeat, + atoms_per_residue_vals, repeating_unit_offset, failed ); + + if ( failed ) { + TR.Warning << "rebuild_parametric_backbone: Failed to rebuild parametric element " + << p << " in ParametersSet " << ps << std::endl; + continue; + } + + protocols::helical_bundle::place_atom_positions( pose, atom_positions, start, end ); + } + } +} + +} // 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..539018cbe2a --- /dev/null +++ b/source/src/core/optimization/parametric_minimize_util.hh @@ -0,0 +1,94 @@ +// -*- 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 +); + +/// @brief Rebuild backbone coordinates for all parametric elements in the pose. +/// @details Iterates over all ParametersSets, reconstructing backbone atom positions +/// from current parameter values via the ParametrizationCalculator. +/// This is a potentially expensive operation (full Crick equation evaluation). +void rebuild_parametric_backbone( + pose::Pose & pose +); + +} // namespace optimization +} // namespace core + +#endif // INCLUDED_core_optimization_parametric_minimize_util_hh 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/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/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..d82b3f60b7c --- /dev/null +++ b/tests/integration/tests/parametric_minimize/command @@ -0,0 +1,30 @@ +# +# 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 +# + +cd %(workdir)s + +[ -x %(bin)s/rosetta_scripts.%(binext)s ] || exit 1 + +# Run 1: parametric-only minimization +%(bin)s/rosetta_scripts.%(binext)s %(additional_flags)s @flags \ + -parser:protocol inputs/param_only.xml \ + -out:prefix param_only_ \ + -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 +%(bin)s/rosetta_scripts.%(binext)s %(additional_flags)s @flags \ + -parser:protocol inputs/param_plus_chi.xml \ + -out:prefix param_plus_chi_ \ + -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..a59650654ed --- /dev/null +++ b/tests/integration/tests/parametric_minimize/flags @@ -0,0 +1,3 @@ +-nstruct 1 +-overwrite +-out:file:scorefile score.sc 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..07b4750635e --- /dev/null +++ b/tests/integration/tests/parametric_minimize/inputs/param_only.xml @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + 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..946e0716596 --- /dev/null +++ b/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + From 5e24380a22cc2c51890cde93d2f86b56705c701e Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Fri, 1 May 2026 12:52:57 -0700 Subject: [PATCH 02/13] warning gone --- source/tools/build/basic.settings | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) 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" : { From 9a5d9fed9f9321d4e223cd9cc93934de656bccf8 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Fri, 1 May 2026 14:03:26 -0700 Subject: [PATCH 03/13] get building --- .../core/optimization/AtomTreeMinimizer.cc | 45 +++++++++++++------ .../ParametricAtomTreeMultifunc.cc | 6 +-- .../constraints/SequenceProfileConstraint.cc | 3 +- .../scores/AtomBasedConstraintsScore.cc | 10 ----- .../scores/AtomBasedConstraintsScore.hh | 8 +--- source/src/protocols/sparta/PDB.cc | 2 +- 6 files changed, 39 insertions(+), 35 deletions(-) diff --git a/source/src/core/optimization/AtomTreeMinimizer.cc b/source/src/core/optimization/AtomTreeMinimizer.cc index 62995a06993..a0cd1142034 100644 --- a/source/src/core/optimization/AtomTreeMinimizer.cc +++ b/source/src/core/optimization/AtomTreeMinimizer.cc @@ -121,8 +121,8 @@ AtomTreeMinimizer::run( if ( use_parametric && !parametric_dofs.empty() ) { ParametricAtomTreeMultifunc f( pose, min_map, scorefxn, parametric_dofs, options.deriv_check(), options.deriv_check_verbose() ); - Multivec dofs( f.total_dofs() ); + min_map.copy_dofs_from_pose( pose, dofs ); for ( Size p = 1; p <= parametric_dofs.size(); ++p ) { dofs[ min_map.nangles() + p ] = get_parametric_dof_value( pose, parametric_dofs[p] ); @@ -136,6 +136,22 @@ AtomTreeMinimizer::run( Minimizer minimizer( f, options ); minimizer.run( dofs ); end_func = f( dofs ); + + if ( TR.Debug.visible() && ! use_nblist ) { + TR.Debug << "end_func: " << end_func << std::endl; + pose.energies().show( TR.Trace ); + } + + // turn off nblist + if ( use_nblist ) pose.energies().reset_nblist(); + + // if we were doing rigid-body minimization, fold the rotation and + // translation offsets into the jump transforms + // + // also sets rb dofs to 0.0, so in principle func value should be the same + // + min_map.reset_jump_rb_deltas( pose, dofs ); + } else { AtomTreeMultifunc f( pose, min_map, scorefxn, options.deriv_check(), options.deriv_check_verbose() ); @@ -153,22 +169,23 @@ AtomTreeMinimizer::run( Minimizer minimizer( f, options ); minimizer.run( dofs ); end_func = f( dofs ); - } + + if ( TR.Debug.visible() && ! use_nblist ) { + TR.Debug << "end_func: " << end_func << std::endl; + pose.energies().show( TR.Trace ); + } - if ( TR.Debug.visible() && ! use_nblist ) { - TR.Debug << "end_func: " << end_func << std::endl; - pose.energies().show( TR.Trace ); - } + // turn off nblist + if ( use_nblist ) pose.energies().reset_nblist(); - // turn off nblist - if ( use_nblist ) pose.energies().reset_nblist(); + // if we were doing rigid-body minimization, fold the rotation and + // translation offsets into the jump transforms + // + // also sets rb dofs to 0.0, so in principle func value should be the same + // + min_map.reset_jump_rb_deltas( pose, dofs ); - // if we were doing rigid-body minimization, fold the rotation and - // translation offsets into the jump transforms - // - // also sets rb dofs to 0.0, so in principle func value should be the same - // - min_map.reset_jump_rb_deltas( pose, dofs ); + } // rescore Real const end_score( scorefxn( pose ) ); diff --git a/source/src/core/optimization/ParametricAtomTreeMultifunc.cc b/source/src/core/optimization/ParametricAtomTreeMultifunc.cc index 923c7884e28..363ca7a0d94 100644 --- a/source/src/core/optimization/ParametricAtomTreeMultifunc.cc +++ b/source/src/core/optimization/ParametricAtomTreeMultifunc.cc @@ -116,7 +116,7 @@ ParametricAtomTreeMultifunc::dfunc( Multivec const & vars, Multivec & dE_dvars ) 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; + //Size residues_per_repeat = 1; // Read parameter values by iterating and checking names for ( Size i = 1; i <= params->num_parameters(); ++i ) { @@ -140,9 +140,9 @@ ParametricAtomTreeMultifunc::dfunc( Multivec const & vars, Multivec & dE_dvars ) 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" ) { + } /*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; 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/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/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."; } From ae639ae8cc56a638e2783f4f1ca61d9c4f2d0f05 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Fri, 1 May 2026 21:40:39 +0000 Subject: [PATCH 04/13] Fix library layer violation and compilation error in parametric minimization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Remove protocols/helical_bundle/util.hh include from core-level parametric_minimize_util.cc — core cannot depend on protocols - Rewrite rebuild_parametric_backbone() to call numeric::crick_equations::XYZ_BUNDLE() directly and set atom coordinates via pose.set_xyz(), using only core/numeric dependencies - Fix dofs scoping in AtomTreeMinimizer::run() (variable declared inside if/else branches but used after) - All modified files pass g++ -fsyntax-only with no errors --- .../core/optimization/AtomTreeMinimizer.cc | 5 +- .../optimization/parametric_minimize_util.cc | 63 +++++++++++++------ 2 files changed, 47 insertions(+), 21 deletions(-) diff --git a/source/src/core/optimization/AtomTreeMinimizer.cc b/source/src/core/optimization/AtomTreeMinimizer.cc index 62995a06993..03a01392440 100644 --- a/source/src/core/optimization/AtomTreeMinimizer.cc +++ b/source/src/core/optimization/AtomTreeMinimizer.cc @@ -117,12 +117,13 @@ AtomTreeMinimizer::run( scorefxn.setup_for_minimizing( pose, min_map ); // Setup the multifunc — use parametric version if parametric DOFs are present + Multivec dofs; Real start_func, end_func; if ( use_parametric && !parametric_dofs.empty() ) { ParametricAtomTreeMultifunc f( pose, min_map, scorefxn, parametric_dofs, options.deriv_check(), options.deriv_check_verbose() ); - Multivec dofs( f.total_dofs() ); + dofs.resize( f.total_dofs() ); min_map.copy_dofs_from_pose( pose, dofs ); for ( Size p = 1; p <= parametric_dofs.size(); ++p ) { dofs[ min_map.nangles() + p ] = get_parametric_dof_value( pose, parametric_dofs[p] ); @@ -142,7 +143,7 @@ AtomTreeMinimizer::run( if ( deriv_check_result_ ) f.set_deriv_check_result( deriv_check_result_ ); - Multivec dofs( min_map.nangles() ); + dofs.resize( min_map.nangles() ); min_map.copy_dofs_from_pose( pose, dofs ); start_func = f( dofs ); diff --git a/source/src/core/optimization/parametric_minimize_util.cc b/source/src/core/optimization/parametric_minimize_util.cc index efa64a1b8aa..55523763862 100644 --- a/source/src/core/optimization/parametric_minimize_util.cc +++ b/source/src/core/optimization/parametric_minimize_util.cc @@ -26,7 +26,10 @@ #include #include -#include +#include +#include + +#include // Basic headers #include @@ -298,25 +301,47 @@ rebuild_parametric_backbone( if ( r1_vals.empty() ) continue; - // Use helical_bundle utilities for the actual rebuild - utility::vector1< utility::vector1< numeric::xyzVector< Real > > > atom_positions; - bool failed = false; - protocols::helical_bundle::generate_atom_positions( - atom_positions, pose, start, end, - r0_val, omega0_val, delta_omega0_val, delta_t_val, - 0.0 /*z1_offset*/, 0.0 /*z0_offset*/, epsilon_val, - invert, r1_vals, omega1_val, z1_val, - delta_omega1_vals, delta_omega1_all_val, - delta_z1_vals, residues_per_repeat, - atoms_per_residue_vals, repeating_unit_offset, failed ); - - if ( failed ) { - TR.Warning << "rebuild_parametric_backbone: Failed to rebuild parametric element " - << p << " in ParametersSet " << ps << std::endl; - continue; + // Compute omega1 relative to omega0 (matching generate_atom_positions convention) + Real const omega1_relative = omega1_val - omega0_val; + + // Compute atom positions directly via Crick equations (numeric layer, no protocols dependency). + core::Size const helix_length = end - start + 1; + Real t = -static_cast( helix_length + 2 ) / 2.0 + delta_t_val; + + bool rebuild_failed = false; + core::Size atom_counter = 0; + for ( core::Size resid = start; resid <= end && !rebuild_failed; ++resid ) { + core::Size const n_mc = pose.residue( resid ).n_mainchain_atoms(); + for ( core::Size iatom = 1; iatom <= n_mc && !rebuild_failed; ++iatom ) { + ++atom_counter; + core::Size const idx = ((atom_counter - 1) % r1_vals.size()) + 1; + + Real const r1 = r1_vals[ idx ]; + Real const dw1 = delta_omega1_vals[ idx ] + delta_omega1_all_val; + Real const dz1 = delta_z1_vals[ idx ]; + + bool failed = false; + numeric::xyzVector< Real > xyz = numeric::crick_equations::XYZ_BUNDLE( + t, r0_val, omega0_val, delta_omega0_val, + r1, omega1_relative, z1_val, dw1, dz1, epsilon_val, failed ); + + if ( failed ) { + TR.Warning << "rebuild_parametric_backbone: Crick equation failed for element " + << p << " residue " << resid << " atom " << iatom << std::endl; + rebuild_failed = true; + break; + } + + if ( invert ) { + xyz.x( -xyz.x() ); + xyz.z( -xyz.z() ); + } + + core::Size const real_atomno = pose.residue_type( resid ).mainchain_atom( iatom ); + pose.set_xyz( id::AtomID( real_atomno, resid ), xyz ); + } + t += 1.0; } - - protocols::helical_bundle::place_atom_positions( pose, atom_positions, start, end ); } } } From f971f5794aa3de2152dc17bc8f81f45c7bbd2fc5 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Fri, 1 May 2026 22:58:14 +0000 Subject: [PATCH 05/13] Add parametric DOF support to MoveMap/MoveMapFactory XML schema The sub-element was not recognized by the MoveMap XSD. Instead, add parametric as a top-level attribute (like bb, chi, jump) on both MoveMapFactory and the legacy MoveMap: or in MoveMapFactory style: Changes: - MoveMapFactory: add all_parametric(bool), data members, parse in parse_my_tag(), apply in edit_movemap_given_pose() - MoveMapFactory XSD: add "parametric" attribute - Legacy MoveMap schema (rosetta_scripts/util.cc): add "parametric" attribute and parse in parse_movemap_tag() - Integration test XMLs: use attribute form --- source/src/core/select/movemap/MoveMapFactory.cc | 9 ++++++++- source/src/core/select/movemap/MoveMapFactory.hh | 9 ++++++--- source/src/protocols/rosetta_scripts/util.cc | 4 +++- .../tests/parametric_minimize/inputs/param_only.xml | 8 ++------ .../tests/parametric_minimize/inputs/param_plus_chi.xml | 8 ++------ 5 files changed, 21 insertions(+), 17 deletions(-) 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/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/tests/integration/tests/parametric_minimize/inputs/param_only.xml b/tests/integration/tests/parametric_minimize/inputs/param_only.xml index 07b4750635e..2ba188daffa 100644 --- a/tests/integration/tests/parametric_minimize/inputs/param_only.xml +++ b/tests/integration/tests/parametric_minimize/inputs/param_only.xml @@ -14,12 +14,8 @@ - - - - - + type="lbfgs_armijo_nonmonotone" tolerance="0.01" max_iter="200" + bb="false" chi="false" parametric="true"/> diff --git a/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml b/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml index 946e0716596..6f4d61c4b17 100644 --- a/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml +++ b/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml @@ -14,12 +14,8 @@ - - - - - + type="lbfgs_armijo_nonmonotone" tolerance="0.01" max_iter="200" + bb="false" chi="true" parametric="true"/> From eabeff84c4f83dd0005927a6ca8dc1d3d86cf721 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Sat, 2 May 2026 00:03:21 +0000 Subject: [PATCH 06/13] Add parametric attribute to MinMover XSD schema and parsing --- source/src/protocols/minimization_packing/MinMover.cc | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/source/src/protocols/minimization_packing/MinMover.cc b/source/src/protocols/minimization_packing/MinMover.cc index 97739185a3b..159c5650fb9 100644 --- a/source/src/protocols/minimization_packing/MinMover.cc +++ b/source/src/protocols/minimization_packing/MinMover.cc @@ -449,6 +449,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 +562,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" ) From 79740c2d078ef55c8e1189857095cbb0eab0dec5 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Fri, 1 May 2026 17:04:56 -0700 Subject: [PATCH 07/13] empty pose --- tests/integration/tests/parametric_minimize/flags | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/integration/tests/parametric_minimize/flags b/tests/integration/tests/parametric_minimize/flags index a59650654ed..1485df056f7 100644 --- a/tests/integration/tests/parametric_minimize/flags +++ b/tests/integration/tests/parametric_minimize/flags @@ -1,3 +1,4 @@ -nstruct 1 -overwrite -out:file:scorefile score.sc +-input_empty_pose true From 557a0832b9b7639c3176a88193f2494b76c3ba3f Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Fri, 1 May 2026 17:50:47 -0700 Subject: [PATCH 08/13] parameter fix --- .../core/conformation/parametric/SizeValuedParameter.cc | 9 ++++++++- .../helical_bundle/BundleParametrizationCalculator.cc | 2 +- 2 files changed, 9 insertions(+), 2 deletions(-) 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/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 From 68c1865739e5bd7aa812b9368a0b4678d16ffb60 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Sat, 2 May 2026 00:54:41 +0000 Subject: [PATCH 09/13] Fix integration test: thread realistic sequence onto poly-ALA bundle MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit MakeBundle builds poly-ALA; use SimpleThreadingMover to install ELLKAIA heptad repeat (Glu/Leu/Lys/Ala/Ile — gives chi angles and inter-helix packing). Add DumpPdb before minimization for comparison. --- .../parametric_minimize/inputs/param_only.xml | 15 ++++++++++++++- .../parametric_minimize/inputs/param_plus_chi.xml | 15 ++++++++++++++- 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/tests/integration/tests/parametric_minimize/inputs/param_only.xml b/tests/integration/tests/parametric_minimize/inputs/param_only.xml index 2ba188daffa..0ba81054bf4 100644 --- a/tests/integration/tests/parametric_minimize/inputs/param_only.xml +++ b/tests/integration/tests/parametric_minimize/inputs/param_only.xml @@ -6,13 +6,22 @@ + + + + + + + @@ -20,6 +29,10 @@ + + + + diff --git a/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml b/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml index 6f4d61c4b17..d4e33dbe370 100644 --- a/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml +++ b/tests/integration/tests/parametric_minimize/inputs/param_plus_chi.xml @@ -6,13 +6,22 @@ + + + + + + + @@ -20,6 +29,10 @@ + + + + From 66ea3c3c82c9b9a8645bd6966946f26478975b02 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Sat, 2 May 2026 01:02:06 +0000 Subject: [PATCH 10/13] Fix sidechain drift during parametric minimization, add trajectory dump Sidechain fix: after setting backbone XYZ from Crick equations, sync the atom tree's internal coordinates by reading phi/psi/omega back from the new Cartesian positions and writing them as torsions. Without this, the atom tree retains the old backbone frame and sidechains don't follow the backbone when it moves. Trajectory visualization: ParametricAtomTreeMultifunc can now dump a PDB at each function evaluation via set_trajectory_dump(). Enabled automatically when the AtomTreeMinimizer tracer is at Debug level (-out:levels core.optimization.AtomTreeMinimizer:debug). --- .../core/optimization/AtomTreeMinimizer.cc | 6 ++++++ .../ParametricAtomTreeMultifunc.cc | 7 +++++++ .../ParametricAtomTreeMultifunc.hh | 10 +++++++++ .../optimization/parametric_minimize_util.cc | 21 +++++++++++++++++++ .../tests/parametric_minimize/command | 8 +++++-- 5 files changed, 50 insertions(+), 2 deletions(-) diff --git a/source/src/core/optimization/AtomTreeMinimizer.cc b/source/src/core/optimization/AtomTreeMinimizer.cc index 6697e89f68d..94261d1a949 100644 --- a/source/src/core/optimization/AtomTreeMinimizer.cc +++ b/source/src/core/optimization/AtomTreeMinimizer.cc @@ -122,6 +122,12 @@ AtomTreeMinimizer::run( if ( use_parametric && !parametric_dofs.empty() ) { ParametricAtomTreeMultifunc f( pose, min_map, scorefxn, parametric_dofs, options.deriv_check(), options.deriv_check_verbose() ); + + // Enable trajectory dumping at Debug tracer level + if ( TR.Debug.visible() ) { + f.set_trajectory_dump( "parametric_traj", 1 ); + } + dofs.resize( f.total_dofs() ); min_map.copy_dofs_from_pose( pose, dofs ); for ( Size p = 1; p <= parametric_dofs.size(); ++p ) { diff --git a/source/src/core/optimization/ParametricAtomTreeMultifunc.cc b/source/src/core/optimization/ParametricAtomTreeMultifunc.cc index 363ca7a0d94..3e26bbde436 100644 --- a/source/src/core/optimization/ParametricAtomTreeMultifunc.cc +++ b/source/src/core/optimization/ParametricAtomTreeMultifunc.cc @@ -76,6 +76,13 @@ ParametricAtomTreeMultifunc::operator()( Multivec const & vars ) const { 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; } diff --git a/source/src/core/optimization/ParametricAtomTreeMultifunc.hh b/source/src/core/optimization/ParametricAtomTreeMultifunc.hh index a8f2b4b0917..ade36471280 100644 --- a/source/src/core/optimization/ParametricAtomTreeMultifunc.hh +++ b/source/src/core/optimization/ParametricAtomTreeMultifunc.hh @@ -57,6 +57,12 @@ public: 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; @@ -71,6 +77,10 @@ private: bool deriv_check_; bool deriv_check_verbose_; + mutable Size eval_count_ = 0; + std::string trajectory_prefix_; + Size trajectory_stride_ = 0; + }; // ParametricAtomTreeMultifunc } // namespace optimization diff --git a/source/src/core/optimization/parametric_minimize_util.cc b/source/src/core/optimization/parametric_minimize_util.cc index eca6798c02c..89e34a64225 100644 --- a/source/src/core/optimization/parametric_minimize_util.cc +++ b/source/src/core/optimization/parametric_minimize_util.cc @@ -28,6 +28,7 @@ #include #include +#include #include @@ -342,6 +343,26 @@ rebuild_parametric_backbone( } t += 1.0; } + + if ( rebuild_failed ) continue; + + // Sync the atom tree's internal coordinates with the new backbone XYZ. + // Without this, sidechain atoms stay at their old positions because + // the atom tree still has the pre-move backbone torsions and rebuilds + // sidechains from those when chi angles are applied. + // Reading phi/psi/omega back from the new XYZ and writing them as + // torsions forces the atom tree to adopt the new backbone frame. + for ( core::Size resid = start; resid <= end; ++resid ) { + if ( resid > start ) { + pose.set_phi( resid, pose.phi( resid ) ); + } + if ( resid < end ) { + pose.set_psi( resid, pose.psi( resid ) ); + } + if ( resid < end ) { + pose.set_omega( resid, pose.omega( resid ) ); + } + } } } } diff --git a/tests/integration/tests/parametric_minimize/command b/tests/integration/tests/parametric_minimize/command index d82b3f60b7c..799d34b3376 100644 --- a/tests/integration/tests/parametric_minimize/command +++ b/tests/integration/tests/parametric_minimize/command @@ -4,25 +4,29 @@ # 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 +# 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 +# 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 From 21d90cb52d1abe277b363e833609842144201b21 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Sat, 2 May 2026 01:17:02 +0000 Subject: [PATCH 11/13] Fix sidechain rebuild: use set_phi/psi/omega to sync atom tree The previous approach (calling private/protected AtomTree methods) didn't compile. Use public Pose::set_phi/set_psi/set_omega instead: reading each backbone torsion from the new XYZ and writing it back via set_torsion updates the atom tree's internal coordinates, which triggers a coordinate rebuild that repositions all downstream atoms including sidechains. --- .../optimization/parametric_minimize_util.cc | 36 ++++++++++--------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/source/src/core/optimization/parametric_minimize_util.cc b/source/src/core/optimization/parametric_minimize_util.cc index 89e34a64225..6bf1fe6daf0 100644 --- a/source/src/core/optimization/parametric_minimize_util.cc +++ b/source/src/core/optimization/parametric_minimize_util.cc @@ -28,7 +28,7 @@ #include #include -#include +#include #include @@ -344,24 +344,26 @@ rebuild_parametric_backbone( t += 1.0; } - if ( rebuild_failed ) continue; + } + } - // Sync the atom tree's internal coordinates with the new backbone XYZ. - // Without this, sidechain atoms stay at their old positions because - // the atom tree still has the pre-move backbone torsions and rebuilds - // sidechains from those when chi angles are applied. - // Reading phi/psi/omega back from the new XYZ and writing them as - // torsions forces the atom tree to adopt the new backbone frame. + // After all parametric backbone atoms have been placed via set_xyz(), + // sync the atom tree by reading backbone torsions from the new XYZ and + // writing them back via set_torsion. This updates the internal coordinate + // representation and triggers a coordinate rebuild that repositions + // sidechain atoms to match the new backbone frame. + for ( core::Size ps = 1; ps <= pose.conformation().n_parameters_sets(); ++ps ) { + ParametersSetCOP ps_obj = pose.conformation().parameters_set( ps ); + if ( ps_obj == nullptr ) continue; + for ( core::Size p = 1; p <= ps_obj->n_parameters(); ++p ) { + ParametersCOP params = ps_obj->parameters( p ); + if ( params == nullptr || params->n_residue() == 0 ) continue; + core::Size const start = params->first_residue_index(); + core::Size const end = params->last_residue_index(); for ( core::Size resid = start; resid <= end; ++resid ) { - if ( resid > start ) { - pose.set_phi( resid, pose.phi( resid ) ); - } - if ( resid < end ) { - pose.set_psi( resid, pose.psi( resid ) ); - } - if ( resid < end ) { - pose.set_omega( resid, pose.omega( resid ) ); - } + pose.set_phi( resid, pose.phi( resid ) ); + pose.set_psi( resid, pose.psi( resid ) ); + pose.set_omega( resid, pose.omega( resid ) ); } } } From b980b67e983c155a2521411cd6bb336ba257baa2 Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Sat, 2 May 2026 01:26:38 +0000 Subject: [PATCH 12/13] Use build_helix() for backbone rebuild during parametric minimization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The previous approach (setting backbone XYZ then trying to sync the atom tree) failed to correctly rebuild sidechain atoms. The fix: - ParametricAtomTreeMultifunc now takes a RebuildCallback (std::function) that is called to rebuild the full pose from current parameter values - MinMover::inner_run_minimizer() provides this callback — it creates a temporary BundleParametrizationCalculator and calls build_helix(), which correctly places ALL atoms (mainchain + sidechain) via the Crick equations + atom tree operations - Parametric dispatch moved from AtomTreeMinimizer (core/) to MinMover (protocols/), respecting library layering - AtomTreeMinimizer reverted to original non-parametric code - Removed broken rebuild_parametric_backbone() from parametric_minimize_util --- .../core/optimization/AtomTreeMinimizer.cc | 123 ++++----------- .../ParametricAtomTreeMultifunc.cc | 9 +- .../ParametricAtomTreeMultifunc.hh | 10 ++ .../optimization/parametric_minimize_util.cc | 144 ------------------ .../optimization/parametric_minimize_util.hh | 8 - .../minimization_packing/MinMover.cc | 86 ++++++++++- 6 files changed, 132 insertions(+), 248 deletions(-) diff --git a/source/src/core/optimization/AtomTreeMinimizer.cc b/source/src/core/optimization/AtomTreeMinimizer.cc index 94261d1a949..245904ff794 100644 --- a/source/src/core/optimization/AtomTreeMinimizer.cc +++ b/source/src/core/optimization/AtomTreeMinimizer.cc @@ -22,12 +22,9 @@ #include #include #include -#include -#include // Project headers #include -#include #include #include #include @@ -82,31 +79,9 @@ AtomTreeMinimizer::run( pose.energies().show( TR.Trace ); } - // Check whether parametric DOFs should be included - bool const use_parametric = move_map.get_parametric() && pose.conformation().n_parameters_sets() > 0; - utility::vector1< ParametricDOFInfo > parametric_dofs; - std::set< Size > parametric_residues; - if ( use_parametric ) { - enumerate_parametric_dofs( pose, parametric_dofs ); - parametric_residues = get_parametric_residues( pose ); - if ( TR.Debug.visible() ) { - TR.Debug << "Parametric minimization: " << parametric_dofs.size() << " parametric DOFs, " - << parametric_residues.size() << " residues under parametric control." << std::endl; - } - } - - // Setup the map of the degrees of freedom. - // If parametric minimization is active, we need a MoveMap that excludes backbone - // torsions of parametric residues (to prevent redundant DOF control). - kinematics::MoveMap effective_move_map( move_map ); - if ( use_parametric ) { - for ( Size resid : parametric_residues ) { - effective_move_map.set_bb( resid, false ); - } - } - + // setup the map of the degrees of freedom MinimizerMap min_map; - min_map.setup( pose, effective_move_map ); + min_map.setup( pose, move_map ); // if we are using the nblist, set it up if ( use_nblist ) { @@ -116,82 +91,42 @@ AtomTreeMinimizer::run( scorefxn.setup_for_minimizing( pose, min_map ); - // Setup the multifunc — use parametric version if parametric DOFs are present - Multivec dofs; - Real start_func, end_func; - if ( use_parametric && !parametric_dofs.empty() ) { - ParametricAtomTreeMultifunc f( pose, min_map, scorefxn, parametric_dofs, - options.deriv_check(), options.deriv_check_verbose() ); - - // Enable trajectory dumping at Debug tracer level - if ( TR.Debug.visible() ) { - f.set_trajectory_dump( "parametric_traj", 1 ); - } - - dofs.resize( f.total_dofs() ); - min_map.copy_dofs_from_pose( pose, dofs ); - for ( Size p = 1; p <= parametric_dofs.size(); ++p ) { - dofs[ min_map.nangles() + p ] = get_parametric_dof_value( pose, parametric_dofs[p] ); - } - - start_func = f( dofs ); - if ( TR.Trace.visible() && !use_nblist ) { - pose.energies().show( TR.Trace ); - } + // setup the function that we will pass to the low-level minimizer + AtomTreeMultifunc f( pose, min_map, scorefxn, + options.deriv_check(), options.deriv_check_verbose() ); - Minimizer minimizer( f, options ); - minimizer.run( dofs ); - end_func = f( dofs ); - - if ( TR.Debug.visible() && ! use_nblist ) { - TR.Debug << "end_func: " << end_func << std::endl; - pose.energies().show( TR.Trace ); - } + if ( deriv_check_result_ ) f.set_deriv_check_result( deriv_check_result_ ); - // turn off nblist - if ( use_nblist ) pose.energies().reset_nblist(); - - // if we were doing rigid-body minimization, fold the rotation and - // translation offsets into the jump transforms - // - // also sets rb dofs to 0.0, so in principle func value should be the same - // - min_map.reset_jump_rb_deltas( pose, dofs ); + // starting position -- "dofs" = Degrees Of Freedom + Multivec dofs( min_map.nangles() ); + min_map.copy_dofs_from_pose( pose, dofs ); - } else { - AtomTreeMultifunc f( pose, min_map, scorefxn, - options.deriv_check(), options.deriv_check_verbose() ); + Real const start_func( f( dofs ) ); - if ( deriv_check_result_ ) f.set_deriv_check_result( deriv_check_result_ ); + if ( TR.Trace.visible() && ! use_nblist ) { + pose.energies().show( TR.Trace ); + } - dofs.resize( min_map.nangles() ); - min_map.copy_dofs_from_pose( pose, dofs ); + // now do the optimization with the low-level minimizer function + Minimizer minimizer( f, options ); + minimizer.run( dofs ); - start_func = f( dofs ); - if ( TR.Trace.visible() && !use_nblist ) { - pose.energies().show( TR.Trace ); - } + Real const end_func( f( dofs ) ); - Minimizer minimizer( f, options ); - minimizer.run( dofs ); - end_func = f( dofs ); - - if ( TR.Debug.visible() && ! use_nblist ) { - TR.Debug << "end_func: " << end_func << std::endl; - pose.energies().show( TR.Trace ); - } - - // turn off nblist - if ( use_nblist ) pose.energies().reset_nblist(); + if ( TR.Debug.visible() && ! use_nblist ) { + TR.Debug << "end_func: " << end_func << std::endl; + pose.energies().show( TR.Trace ); + } - // if we were doing rigid-body minimization, fold the rotation and - // translation offsets into the jump transforms - // - // also sets rb dofs to 0.0, so in principle func value should be the same - // - min_map.reset_jump_rb_deltas( pose, dofs ); + // turn off nblist + if ( use_nblist ) pose.energies().reset_nblist(); - } + // if we were doing rigid-body minimization, fold the rotation and + // translation offsets into the jump transforms + // + // also sets rb dofs to 0.0, so in principle func value should be the same + // + min_map.reset_jump_rb_deltas( pose, dofs ); // rescore Real const end_score( scorefxn( pose ) ); diff --git a/source/src/core/optimization/ParametricAtomTreeMultifunc.cc b/source/src/core/optimization/ParametricAtomTreeMultifunc.cc index 3e26bbde436..8b57e38c2e9 100644 --- a/source/src/core/optimization/ParametricAtomTreeMultifunc.cc +++ b/source/src/core/optimization/ParametricAtomTreeMultifunc.cc @@ -45,6 +45,7 @@ ParametricAtomTreeMultifunc::ParametricAtomTreeMultifunc( 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 ) : @@ -53,6 +54,7 @@ ParametricAtomTreeMultifunc::ParametricAtomTreeMultifunc( 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 ) {} @@ -65,7 +67,12 @@ ParametricAtomTreeMultifunc::apply_parametric_dofs( Multivec const & vars ) cons Real const val = vars[ n_standard_dofs_ + p ]; set_parametric_dof_value( pose_, parametric_dofs_[p], val ); } - rebuild_parametric_backbone( pose_ ); + // 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 diff --git a/source/src/core/optimization/ParametricAtomTreeMultifunc.hh b/source/src/core/optimization/ParametricAtomTreeMultifunc.hh index ade36471280..8cb4bc258d8 100644 --- a/source/src/core/optimization/ParametricAtomTreeMultifunc.hh +++ b/source/src/core/optimization/ParametricAtomTreeMultifunc.hh @@ -28,6 +28,7 @@ #include #include +#include namespace core { namespace optimization { @@ -36,11 +37,18 @@ 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 ); @@ -74,6 +82,8 @@ private: utility::vector1< ParametricDOFInfo > parametric_dofs_; Size n_standard_dofs_; + RebuildCallback rebuild_callback_; + bool deriv_check_; bool deriv_check_verbose_; diff --git a/source/src/core/optimization/parametric_minimize_util.cc b/source/src/core/optimization/parametric_minimize_util.cc index 6bf1fe6daf0..506b3274446 100644 --- a/source/src/core/optimization/parametric_minimize_util.cc +++ b/source/src/core/optimization/parametric_minimize_util.cc @@ -21,16 +21,7 @@ #include #include #include -#include -#include -#include -#include - #include -#include -#include - -#include // Basic headers #include @@ -234,140 +225,5 @@ set_parametric_dof_value( real_param->set_value( value, true /*ignore_use_degrees*/ ); } -void -rebuild_parametric_backbone( - pose::Pose & pose -) { - using namespace core::conformation::parametric; - - for ( core::Size ps = 1; ps <= pose.conformation().n_parameters_sets(); ++ps ) { - ParametersSetCOP params_set = pose.conformation().parameters_set( ps ); - if ( params_set == nullptr ) continue; - for ( core::Size p = 1; p <= params_set->n_parameters(); ++p ) { - ParametersCOP params = params_set->parameters( p ); - if ( params == nullptr ) continue; - if ( params->n_residue() == 0 ) continue; - - core::Size const start = params->first_residue_index(); - core::Size const end = params->last_residue_index(); - - // Rebuild this parametric element's backbone coordinates. - // We use the helical_bundle utility functions directly: generate new atom - // positions from current parameter values, then place them. - // For now, this is a simplified rebuild that re-applies the Crick equations. - // A more complete version would use the ParametrizationCalculator's build_helix/build_strand. - - // Extract parameters by name for the Crick equation call - Real r0_val = 0, omega0_val = 0, delta_omega0_val = 0; - Real omega1_val = 0, z1_val = 0, delta_omega1_all_val = 0; - Real delta_t_val = 0, epsilon_val = 1.0; - bool invert = false; - utility::vector1< Real > r1_vals, delta_omega1_vals, delta_z1_vals; - core::Size /*residues_per_repeat = 1,*/ repeating_unit_offset = 0; - utility::vector1< core::Size > atoms_per_residue_vals; - - for ( core::Size i = 1; i <= params->num_parameters(); ++i ) { - ParameterCOP param = params->parameter_cop( i ); - std::string const & name = param->parameter_name(); - if ( name == "r0" ) { - r0_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); - } else if ( name == "omega0" ) { - omega0_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); - } else if ( name == "delta_omega0" ) { - delta_omega0_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); - } else if ( name == "omega1" ) { - omega1_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); - } else if ( name == "z1" ) { - z1_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); - } else if ( name == "delta_omega1" ) { - delta_omega1_all_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); - } else if ( name == "delta_t" || name == "delta_z0" ) { - delta_t_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); - } else if ( name == "epsilon" ) { - epsilon_val = utility::pointer::static_pointer_cast< RealValuedParameter const >( param )->value(); - } else if ( name == "invert" ) { - invert = utility::pointer::static_pointer_cast< BooleanValuedParameter const >( param )->value(); - } else if ( name == "r1_peratom" ) { - r1_vals = utility::pointer::static_pointer_cast< RealVectorValuedParameter const >( param )->value(); - } else if ( name == "delta_omega1_peratom" ) { - delta_omega1_vals = utility::pointer::static_pointer_cast< RealVectorValuedParameter const >( param )->value(); - } else if ( name == "delta_z1_peratom" ) { - delta_z1_vals = utility::pointer::static_pointer_cast< RealVectorValuedParameter const >( param )->value(); - } /*else if ( name == "residues_per_repeat" ) { - residues_per_repeat = utility::pointer::static_pointer_cast< SizeValuedParameter const >( param )->value(); - }*/ else if ( name == "atoms_per_residue" ) { - atoms_per_residue_vals = utility::pointer::static_pointer_cast< SizeVectorValuedParameter const >( param )->value(); - } - } - - if ( r1_vals.empty() ) continue; - - // Compute omega1 relative to omega0 (matching generate_atom_positions convention) - Real const omega1_relative = omega1_val - omega0_val; - - // Compute atom positions directly via Crick equations (numeric layer, no protocols dependency). - core::Size const helix_length = end - start + 1; - Real t = -static_cast( helix_length + 2 ) / 2.0 + delta_t_val; - - bool rebuild_failed = false; - core::Size atom_counter = 0; - for ( core::Size resid = start; resid <= end && !rebuild_failed; ++resid ) { - core::Size const n_mc = pose.residue( resid ).n_mainchain_atoms(); - for ( core::Size iatom = 1; iatom <= n_mc && !rebuild_failed; ++iatom ) { - ++atom_counter; - core::Size const idx = ((atom_counter - 1) % r1_vals.size()) + 1; - - Real const r1 = r1_vals[ idx ]; - Real const dw1 = delta_omega1_vals[ idx ] + delta_omega1_all_val; - Real const dz1 = delta_z1_vals[ idx ]; - - bool failed = false; - numeric::xyzVector< Real > xyz = numeric::crick_equations::XYZ_BUNDLE( - t, r0_val, omega0_val, delta_omega0_val, - r1, omega1_relative, z1_val, dw1, dz1, epsilon_val, failed ); - - if ( failed ) { - TR.Warning << "rebuild_parametric_backbone: Crick equation failed for element " - << p << " residue " << resid << " atom " << iatom << std::endl; - rebuild_failed = true; - break; - } - - if ( invert ) { - xyz.x( -xyz.x() ); - xyz.z( -xyz.z() ); - } - - core::Size const real_atomno = pose.residue_type( resid ).mainchain_atom( iatom ); - pose.set_xyz( id::AtomID( real_atomno, resid ), xyz ); - } - t += 1.0; - } - - } - } - - // After all parametric backbone atoms have been placed via set_xyz(), - // sync the atom tree by reading backbone torsions from the new XYZ and - // writing them back via set_torsion. This updates the internal coordinate - // representation and triggers a coordinate rebuild that repositions - // sidechain atoms to match the new backbone frame. - for ( core::Size ps = 1; ps <= pose.conformation().n_parameters_sets(); ++ps ) { - ParametersSetCOP ps_obj = pose.conformation().parameters_set( ps ); - if ( ps_obj == nullptr ) continue; - for ( core::Size p = 1; p <= ps_obj->n_parameters(); ++p ) { - ParametersCOP params = ps_obj->parameters( p ); - if ( params == nullptr || params->n_residue() == 0 ) continue; - core::Size const start = params->first_residue_index(); - core::Size const end = params->last_residue_index(); - for ( core::Size resid = start; resid <= end; ++resid ) { - pose.set_phi( resid, pose.phi( resid ) ); - pose.set_psi( resid, pose.psi( resid ) ); - pose.set_omega( resid, pose.omega( resid ) ); - } - } - } -} - } // namespace optimization } // namespace core diff --git a/source/src/core/optimization/parametric_minimize_util.hh b/source/src/core/optimization/parametric_minimize_util.hh index 539018cbe2a..3911d2e998b 100644 --- a/source/src/core/optimization/parametric_minimize_util.hh +++ b/source/src/core/optimization/parametric_minimize_util.hh @@ -80,14 +80,6 @@ void set_parametric_dof_value( Real value ); -/// @brief Rebuild backbone coordinates for all parametric elements in the pose. -/// @details Iterates over all ParametersSets, reconstructing backbone atom positions -/// from current parameter values via the ParametrizationCalculator. -/// This is a potentially expensive operation (full Crick equation evaluation). -void rebuild_parametric_backbone( - pose::Pose & pose -); - } // namespace optimization } // namespace core diff --git a/source/src/protocols/minimization_packing/MinMover.cc b/source/src/protocols/minimization_packing/MinMover.cc index 159c5650fb9..8c0f6a12576 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 = [&pose]( 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 >(); From 3d65a0eaed7318f77428e83fad3ef96b61dbb68f Mon Sep 17 00:00:00 2001 From: Andy Watkins Date: Sat, 2 May 2026 01:33:29 +0000 Subject: [PATCH 13/13] Remove unused lambda capture --- source/src/protocols/minimization_packing/MinMover.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/src/protocols/minimization_packing/MinMover.cc b/source/src/protocols/minimization_packing/MinMover.cc index 8c0f6a12576..bfdec0e0793 100644 --- a/source/src/protocols/minimization_packing/MinMover.cc +++ b/source/src/protocols/minimization_packing/MinMover.cc @@ -342,7 +342,7 @@ MinMover::inner_run_minimizer( core::pose::Pose & pose, core::kinematics::MoveMa // Create the rebuild callback that properly rebuilds ALL atoms // from the current parameter values by calling build_helix() - auto rebuild_callback = [&pose]( core::pose::Pose & p ) { + 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 );