Add differentiable parametric backbone minimization

source/src/core.5.src.settings

@@ -26,6 +26,8 @@ sources = {
 		"MinimizerOptions",
 		"NelderMeadSimplex",
 		"NumericalDerivCheckResult",
+		"ParametricAtomTreeMultifunc",
+		"parametric_minimize_util",
 		"ParticleSwarmMinimizer",
 		"SingleResidueMultifunc",
 	],

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() ) );
 	}
 }
 

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 );

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  <seqpos>
@@ -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;

source/src/core/optimization/ParametricAtomTreeMultifunc.cc

@@ -0,0 +1,234 @@
+// -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*-
+// vi: set ts=2 noet:
+//
+// (c) Copyright Rosetta Commons Member Institutions.
+// (c) This file is part of the Rosetta software suite and is made available under license.
+// (c) The Rosetta software is developed by the contributing members of the Rosetta Commons.
+// (c) For more information, see http://www.rosettacommons.org. Questions about this can be
+// (c) addressed to University of Washington CoMotion, email: license@uw.edu.
+
+/// @file   core/optimization/ParametricAtomTreeMultifunc.cc
+/// @brief  Multifunc for co-minimizing parametric and atom-tree DOFs.
+/// @author Andy Watkins
+
+#include <core/optimization/ParametricAtomTreeMultifunc.hh>
+#include <core/optimization/MinimizerMap.hh>
+#include <core/optimization/atom_tree_minimize.hh>
+#include <core/optimization/parametric_minimize_util.hh>
+
+#include <core/pose/Pose.hh>
+#include <core/scoring/ScoreFunction.hh>
+#include <core/scoring/DerivVectorPair.hh>
+
+#include <core/conformation/Conformation.hh>
+#include <core/conformation/Residue.hh>
+#include <core/conformation/parametric/Parameters.hh>
+#include <core/conformation/parametric/ParametersSet.hh>
+#include <core/conformation/parametric/RealValuedParameter.hh>
+#include <core/conformation/parametric/RealVectorValuedParameter.hh>
+#include <core/conformation/parametric/SizeValuedParameter.hh>
+#include <core/conformation/parametric/SizeVectorValuedParameter.hh>
+
+#include <numeric/crick_equations/BundleParams.hh>
+#include <numeric/crick_equations/BundleParams_derivatives.hh>
+
+#include <basic/Tracer.hh>
+#include <basic/prof.hh>
+
+namespace core {
+namespace optimization {
+
+static basic::Tracer TR( "core.optimization.ParametricAtomTreeMultifunc" );
+
+ParametricAtomTreeMultifunc::ParametricAtomTreeMultifunc(
+	pose::Pose & pose_in,
+	MinimizerMap & min_map_in,
+	scoring::ScoreFunction const & scorefxn_in,
+	utility::vector1< ParametricDOFInfo > const & parametric_dofs_in,
+	RebuildCallback rebuild_callback,
+	bool const deriv_check_in,
+	bool const deriv_check_verbose_in
+) :
+	pose_( pose_in ),
+	min_map_( min_map_in ),
+	score_function_( scorefxn_in ),
+	parametric_dofs_( parametric_dofs_in ),
+	n_standard_dofs_( min_map_in.nangles() ),
+	rebuild_callback_( std::move( rebuild_callback ) ),
+	deriv_check_( deriv_check_in ),
+	deriv_check_verbose_( deriv_check_verbose_in )
+{}
+
+ParametricAtomTreeMultifunc::~ParametricAtomTreeMultifunc() = default;
+
+void
+ParametricAtomTreeMultifunc::apply_parametric_dofs( Multivec const & vars ) const {
+	for ( Size p = 1; p <= parametric_dofs_.size(); ++p ) {
+		Real const val = vars[ n_standard_dofs_ + p ];
+		set_parametric_dof_value( pose_, parametric_dofs_[p], val );
+	}
+	// Rebuild the full backbone (and sidechain frames) from current parameters.
+	// The callback is provided by protocols-level code and calls the appropriate
+	// ParametrizationCalculator::build_helix/build_strand, which correctly
+	// places all atoms — mainchain and sidechain — via the Crick equations
+	// and atom tree operations.
+	rebuild_callback_( pose_ );
+}
+
+Real
+ParametricAtomTreeMultifunc::operator()( Multivec const & vars ) const {
+	PROF_START( basic::FUNC );
+
+	apply_parametric_dofs( vars );
+	min_map_.copy_dofs_to_pose( pose_, vars );
+	Real const score = score_function_( pose_ );
+
+	++eval_count_;
+	if ( trajectory_stride_ > 0 && ( eval_count_ % trajectory_stride_ == 0 || eval_count_ == 1 ) ) {
+		std::string const fname = trajectory_prefix_ + "_" + std::to_string( eval_count_ ) + ".pdb";
+		pose_.dump_pdb( fname );
+		TR << "Trajectory frame " << eval_count_ << " score=" << score << " -> " << fname << std::endl;
+	}
+
+	PROF_STOP( basic::FUNC );
+	return score;
+}
+
+void
+ParametricAtomTreeMultifunc::dfunc( Multivec const & vars, Multivec & dE_dvars ) const {
+	PROF_START( basic::DFUNC );
+
+	apply_parametric_dofs( vars );
+
+	// Compute standard DOF derivatives via the normal atom_tree_dfunc pipeline.
+	// This fills dE_dvars[1..n_standard_dofs_] and populates
+	// min_map_.atom_derivatives_ with per-atom DerivVectorPairs.
+	atom_tree_dfunc( pose_, min_map_, score_function_, vars, dE_dvars );
+
+	// Extend dE_dvars to include parametric DOF slots
+	Size const ntotal = n_standard_dofs_ + parametric_dofs_.size();
+	dE_dvars.resize( ntotal, 0.0 );
+
+	// Compute parametric DOF derivatives via the Crick Jacobian chain rule:
+	// dE/dParam = Σ_atoms (dE/dXYZ · dXYZ/dParam)
+	for ( Size p = 1; p <= parametric_dofs_.size(); ++p ) {
+		ParametricDOFInfo const & info = parametric_dofs_[p];
+		Real dE_dparam = 0.0;
+
+		// Get current Crick parameter values for Jacobian computation
+		conformation::parametric::ParametersSetCOP params_set =
+			pose_.conformation().parameters_set( info.params_set_index );
+		conformation::parametric::ParametersCOP params =
+			params_set->parameters( info.params_index );
+
+		// Extract the Crick parameters needed for derivative computation.
+		// We need: r0, omega0, delta_omega0, omega1, z1, delta_omega1_all, plus
+		// per-atom r1, delta_omega1, delta_z1.
+		// These are stored as parameters in the Parameters object.
+		// The exact enum indices depend on the calculator type (BPC or BBPC).
+		// For generality, we look up by name.
+		Real r0_val = 0, omega0_val = 0, delta_omega0_val = 0;
+		Real omega1_val = 0, z1_val = 0, delta_omega1_all_val = 0;
+		utility::vector1< Real > r1_vals, delta_omega1_vals, delta_z1_vals;
+		//Size residues_per_repeat = 1;
+
+		// Read parameter values by iterating and checking names
+		for ( Size i = 1; i <= params->num_parameters(); ++i ) {
+			conformation::parametric::ParameterCOP param = params->parameter_cop( i );
+			std::string const & name = param->parameter_name();
+			if ( name == "r0" ) {
+				r0_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value();
+			} else if ( name == "omega0" ) {
+				omega0_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value();
+			} else if ( name == "delta_omega0" ) {
+				delta_omega0_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value();
+			} else if ( name == "omega1" ) {
+				omega1_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value();
+			} else if ( name == "z1" ) {
+				z1_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value();
+			} else if ( name == "delta_omega1" ) {
+				delta_omega1_all_val = utility::pointer::static_pointer_cast< conformation::parametric::RealValuedParameter const >( param )->value();
+			} else if ( name == "r1_peratom" ) {
+				r1_vals = utility::pointer::static_pointer_cast< conformation::parametric::RealVectorValuedParameter const >( param )->value();
+			} else if ( name == "delta_omega1_peratom" ) {
+				delta_omega1_vals = utility::pointer::static_pointer_cast< conformation::parametric::RealVectorValuedParameter const >( param )->value();
+			} else if ( name == "delta_z1_peratom" ) {
+				delta_z1_vals = utility::pointer::static_pointer_cast< conformation::parametric::RealVectorValuedParameter const >( param )->value();
+			} /*else if ( name == "residues_per_repeat" ) {
+				residues_per_repeat = utility::pointer::static_pointer_cast< conformation::parametric::SizeValuedParameter const >( param )->value();
+			}*/
+		}
+
+		if ( r1_vals.empty() || delta_omega1_vals.empty() || delta_z1_vals.empty() ) continue;
+
+		// Compute omega1 relative to omega0 (as done in generate_atom_positions)
+		Real const omega1_relative = omega1_val - omega0_val;
+
+		// Compute t values and iterate over mainchain atoms
+		Size const helix_length = info.helix_end_resid - info.helix_start_resid + 1;
+		Real t = -static_cast<Real>( 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