From fe876bc0a1d2db5a673bd0094ba88f04df84b1fe Mon Sep 17 00:00:00 2001 From: meister Date: Wed, 27 Nov 2024 15:02:50 -0500 Subject: [PATCH] Assorted fixes Changed how energy-scale and orientations work --- include/cando/chem/energyAnchorRestraint.h | 1 + include/cando/chem/energyAngle.h | 1 + include/cando/chem/energyChiralRestraint.h | 1 + include/cando/chem/energyComponent.h | 28 ++ include/cando/chem/energyDihedral.h | 1 + include/cando/chem/energyDihedralRestraint.h | 23 +- include/cando/chem/energyFixedNonbond.h | 23 +- include/cando/chem/energyFunction.h | 61 +---- include/cando/chem/energyNonbond.h | 57 ++-- include/cando/chem/energyOutOfZPlane.h | 23 +- .../energyPeriodicBoundaryConditionsNonbond.h | 2 + .../cando/chem/energyPointToLineRestraint.h | 1 + include/cando/chem/energyRigidBodyComponent.h | 23 +- include/cando/chem/energyRigidBodyNonbond.h | 1 + include/cando/chem/energyRigidBodyStaple.h | 1 + include/cando/chem/energySketchNonbond.h | 2 + include/cando/chem/energySketchStretch.h | 1 + include/cando/chem/energyStretch.h | 1 + include/cando/chem/minimizer.h | 35 ++- include/cando/chem/rigidBodyEnergyFunction.h | 3 +- include/cando/chem/scoringFunction.h | 66 ++++- include/cando/chem/sketchFunction.h | 11 +- include/cando/kinematics/bondedJoint.h | 2 +- include/cando/kinematics/jumpJoint.h | 2 + src/chem/cluster.cc | 1 - src/chem/energyAnchorRestraint.cc | 23 +- src/chem/energyAngle.cc | 1 + src/chem/energyChiralRestraint.cc | 23 +- src/chem/energyComponent.cc | 16 +- src/chem/energyDihedral.cc | 23 +- src/chem/energyDihedralRestraint.cc | 23 +- src/chem/energyFixedNonbond.cc | 9 +- src/chem/energyFunction.cc | 100 +++---- src/chem/energyNonbond.cc | 140 +++++----- src/chem/energyOutOfZPlane.cc | 23 +- ...energyPeriodicBoundaryConditionsNonbond.cc | 6 +- src/chem/energyPointToLineRestraint.cc | 23 +- src/chem/energyRigidBodyNonbond.cc | 25 +- src/chem/energyRigidBodyStaple.cc | 23 +- src/chem/energySketchNonbond.cc | 46 ++-- src/chem/energySketchStretch.cc | 1 + src/chem/energyStretch.cc | 23 +- src/chem/minimizer.cc | 141 ++++++---- src/chem/rigidBodyEnergyFunction.cc | 9 +- src/chem/scoringFunction.cc | 87 ++++-- src/chem/sketchFunction.cc | 48 ++-- src/kinematics/bondedJoint.cc | 12 + src/kinematics/xyzJoint.cc | 2 +- src/lisp/cando/dynamics.lisp | 32 +-- src/lisp/cando/molecules.lisp | 4 +- src/lisp/chem-extras/chem.lisp | 2 +- src/lisp/regression-tests/energy.lisp | 6 +- src/lisp/sketch2d/workbench.lisp | 2 +- src/lisp/topology/assembler.lisp | 252 +++++++++++++----- src/lisp/topology/internals.lisp | 50 ++-- src/lisp/topology/joint-templates.lisp | 2 +- src/lisp/topology/jupyter.lisp | 2 +- src/lisp/topology/linearize-internals.lisp | 7 +- src/lisp/topology/oligomer.lisp | 23 +- src/lisp/topology/packages.lisp | 11 +- src/lisp/topology/shape.lisp | 89 +++---- src/lisp/topology/topology-classes.lisp | 3 +- 62 files changed, 1007 insertions(+), 676 deletions(-) diff --git a/include/cando/chem/energyAnchorRestraint.h b/include/cando/chem/energyAnchorRestraint.h index fa5ba8e8..57a0b037 100644 --- a/include/cando/chem/energyAnchorRestraint.h +++ b/include/cando/chem/energyAnchorRestraint.h @@ -161,6 +161,7 @@ class EnergyAnchorRestraint_O : public EnergyComponent_O core::T_sp activeAtomMask ); virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/energyAngle.h b/include/cando/chem/energyAngle.h index 361979df..6e039239 100644 --- a/include/cando/chem/energyAngle.h +++ b/include/cando/chem/energyAngle.h @@ -204,6 +204,7 @@ class EnergyAngle_O : public EnergyComponent_O core::T_sp activeAtomMask ); virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/energyChiralRestraint.h b/include/cando/chem/energyChiralRestraint.h index 0a2709a3..8d1289e7 100644 --- a/include/cando/chem/energyChiralRestraint.h +++ b/include/cando/chem/energyChiralRestraint.h @@ -174,6 +174,7 @@ class EnergyChiralRestraint_O : public EnergyComponent_O core::T_sp activeAtomMask ); virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/energyComponent.h b/include/cando/chem/energyComponent.h index c216c6ec..c914bcb2 100644 --- a/include/cando/chem/energyComponent.h +++ b/include/cando/chem/energyComponent.h @@ -71,6 +71,30 @@ This is an open source license for the CANDO software from Temple University, bu #endif #include // energyComponent.h wants QDomNode needs quickDom.fwd.h +namespace chem { +class KahanSummation { +private: + double sum; + double c; // A running compensation for lost low-order bits. + +public: + // Constructor initializes sum and compensation to zero. + KahanSummation() : sum(0.0), c(0.0) {} + + // Adds a new value to the running total using Kahan summation algorithm. + void add(double input) { + double y = input - c; + double t = sum + y; + c = (t - sum) - y; + sum = t; + } + + // Returns the current value of the sum. + double getSum() const { + return sum; + } +}; +}; namespace chem { @@ -280,6 +304,9 @@ class EnergyComponent_O : public core::CxxObject_O //protected: // Define these in subclasses // vector _Terms; // vector _BeyondThresholdTerms; +public: + bool fieldsp() const { return true; }; + void fields(core::Record_sp node); public: CL_DEFMETHOD virtual size_t numberOfTerms() {_OF(); SUBCLASS_MUST_IMPLEMENT();}; void setScale(double s) {this->_Scale = s; }; @@ -322,6 +349,7 @@ class EnergyComponent_O : public core::CxxObject_O CL_DEFMETHOD virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/energyDihedral.h b/include/cando/chem/energyDihedral.h index e55ad842..ef069ea7 100644 --- a/include/cando/chem/energyDihedral.h +++ b/include/cando/chem/energyDihedral.h @@ -233,6 +233,7 @@ class EnergyDihedral_O : public EnergyComponent_O core::T_sp activeAtomMask ); virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/energyDihedralRestraint.h b/include/cando/chem/energyDihedralRestraint.h index 0229b085..b48a25a4 100644 --- a/include/cando/chem/energyDihedralRestraint.h +++ b/include/cando/chem/energyDihedralRestraint.h @@ -173,17 +173,18 @@ class EnergyDihedralRestraint_O : public EnergyComponent_O AbstractLargeSquareMatrix_sp m, core::T_sp activeAtomMask ); virtual double evaluateAllComponent( ScoringFunction_sp scorer, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ); + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ); virtual void compareAnalyticalAndNumericalForceAndHessianTermByTerm( NVector_sp pos ); diff --git a/include/cando/chem/energyFixedNonbond.h b/include/cando/chem/energyFixedNonbond.h index bfdbd377..091fe984 100644 --- a/include/cando/chem/energyFixedNonbond.h +++ b/include/cando/chem/energyFixedNonbond.h @@ -198,17 +198,18 @@ class EnergyFixedNonbondRestraint_O : public EnergyComponent_O //virtual void setupHessianPreconditioner(NVector_sp nvPosition, // AbstractLargeSquareMatrix_sp m ); virtual double evaluateAllComponent( ScoringFunction_sp scorer, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) override; + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) override; virtual void compareAnalyticalAndNumericalForceAndHessianTermByTerm( NVector_sp pos ); diff --git a/include/cando/chem/energyFunction.h b/include/cando/chem/energyFunction.h index cca47293..036d230e 100644 --- a/include/cando/chem/energyFunction.h +++ b/include/cando/chem/energyFunction.h @@ -69,7 +69,6 @@ namespace chem FORWARD(FFPtorDb); FORWARD(FFItorDb); FORWARD(FFNonbondDb); - FORWARD(EnergyScale); class EnergyAtom; @@ -96,10 +95,6 @@ namespace chem FORWARD(BoundingBox); -#define DefaultChiralRestraintOffset 0.2 -#define DefaultChiralRestraintWeight 100000.0 -#define DefaultAnchorRestraintWeight 10.0 - core::T_sp specializeKeepInteractionFactory( core::T_sp keepInteractionFactory, core::T_sp aclass ); bool skipInteraction( core::T_sp keepInteractionFunction, @@ -118,47 +113,6 @@ bool skipInteraction( core::T_sp keepInteractionFunction, FORWARD(EnergyFunction); }; -namespace chem { - FORWARD(EnergyScale); -class EnergyScale_O : public core::CxxObject_O - { - LISP_CLASS(chem,ChemPkg,EnergyScale_O,"EnergyScale",core::CxxObject_O); - public: - double _DielectricConstant; - double _NonbondCutoff; - double _ScaleVdw; - double _ScaleElectrostatic; - double _ChiralRestraintWeight; - double _ChiralRestraintOffset; - double _AnchorRestraintWeight; - double _FixedNonbondRestraintWeight; - EnergyScale_O() : _DielectricConstant(1.0), - _NonbondCutoff(16.0), - _ScaleVdw(1.0), - _ScaleElectrostatic(1.0), - _ChiralRestraintWeight(1.0), - _ChiralRestraintOffset(1.0), - _AnchorRestraintWeight(1.0), - _FixedNonbondRestraintWeight(1.0) - { - this->_ChiralRestraintWeight = DefaultChiralRestraintWeight; - this->_ChiralRestraintOffset = DefaultChiralRestraintOffset; - this->_AnchorRestraintWeight = DefaultAnchorRestraintWeight; - }; - static EnergyScale_sp make(); - CL_DEFMETHOD void setVdwScale(double d) { this->_ScaleVdw = d; }; - CL_DEFMETHOD double getVdwScale() {return this->_ScaleVdw; }; - CL_DEFMETHOD void setElectrostaticScale(double d) { this->_ScaleElectrostatic = d; }; - CL_DEFMETHOD double getElectrostaticScale() {return this->_ScaleElectrostatic; }; - CL_DEFMETHOD void setDielectricConstant(double d) { this->_DielectricConstant = d; }; - CL_DEFMETHOD double getDielectricConstant() { return this->_DielectricConstant; }; - CL_DEFMETHOD void setNonbondCutoff(double d) { this->_NonbondCutoff = d; }; - CL_DEFMETHOD double getNonbondCutoff() { return this->_NonbondCutoff; }; - -}; - -}; - template <> struct gctools::GCInfo { @@ -176,7 +130,6 @@ namespace chem { static EnergyFunction_sp make(core::T_sp matter, core::T_sp disableComponents=nil(), core::List_sp enableComponents=nil(), - core::T_sp energyScale=unbound(), bool useExcludedAtoms=false, core::T_sp keepInteractionFactory=nil(), bool assign_types=false ); @@ -187,7 +140,6 @@ namespace chem { void fields(core::Record_sp node); public: Matter_sp _Matter; // Aggregate or Molecule - core::T_sp _EnergyScale; /*! Stores cross terms for evaluating nonbond interactions */ FFNonbondCrossTermTable_sp _NonbondCrossTermTable; @@ -227,8 +179,6 @@ namespace chem { core::List_sp allComponents() const; - EnergyScale_sp energyScale(); - void setEnergyScale(core::T_sp energyScale); string energyTermsEnabled() ; void loadCoordinatesIntoVector(NVector_sp pos); void saveCoordinatesFromVector(NVector_sp pos); @@ -244,7 +194,7 @@ namespace chem { void setBoundingBox(BoundingBox_sp bounding_box); void makUnboundBoundingBox(); - ForceMatchReport_sp checkIfAnalyticalForceMatchesNumericalForce( NVector_sp pos, NVector_sp force, core::T_sp activeAtomMask ); + ForceMatchReport_sp checkIfAnalyticalForceMatchesNumericalForce( NVector_sp pos, core::T_sp energyScale, NVector_sp force, core::T_sp activeAtomMask ); CL_LISPIFY_NAME("getMatter"); CL_DEFMETHOD Matter_sp getMatter() { return this->_Matter;}; @@ -326,9 +276,10 @@ namespace chem { void addTermsForListOfRestraints( ForceField_sp forceField, core::List_sp restraintList, core::T_sp keepInteractionFactory, core::HashTable_sp atomTypes ); - double calculateNumericalDerivative(NVector_sp pos, double delta, uint i, core::T_sp activeAtomMask ); - double calculateNumericalSecondDerivative(NVector_sp pos, double delta, uint i, uint j, core::T_sp activeAtomMask ); + double calculateNumericalDerivative(NVector_sp pos, core::T_sp energyScale, double delta, uint i, core::T_sp activeAtomMask ); + double calculateNumericalSecondDerivative(NVector_sp pos, core::T_sp energyScale, double delta, uint i, uint j, core::T_sp activeAtomMask ); double evaluateAll(NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, @@ -348,8 +299,8 @@ namespace chem { // adapt::QDomNode_sp accumulateTermsBeyondThresholdAsXml(); uint countTermsBeyondThreshold(); - void evaluateNumericalForce(NVector_sp pos, NVector_sp numForce, double delta, core::T_sp activeAtomMask ); - void evaluateNumericalHessian(NVector_sp pos, AbstractLargeSquareMatrix_sp numHessian, bool calcOffDiagonalElements, double delta, core::T_sp activeAtomMask); + void evaluateNumericalForce(NVector_sp pos, core::T_sp energyScale, NVector_sp numForce, double delta, core::T_sp activeAtomMask ); + void evaluateNumericalHessian(NVector_sp pos, core::T_sp energyScale, AbstractLargeSquareMatrix_sp numHessian, bool calcOffDiagonalElements, double delta, core::T_sp activeAtomMask); core::List_sp checkForBeyondThresholdInteractions(double threshold); diff --git a/include/cando/chem/energyNonbond.h b/include/cando/chem/energyNonbond.h index 79e55cf0..3b21d909 100644 --- a/include/cando/chem/energyNonbond.h +++ b/include/cando/chem/energyNonbond.h @@ -142,7 +142,7 @@ struct to_object namespace chem { -double _evaluateEnergyOnly_Nonbond(ScoringFunction_sp score, +double _evaluateEnergyOnly_Nonbond(ScoringFunction_sp score, core::T_sp energyScale, int I1, int I2, core::T_sp activeAtomMask, num_real x1, num_real y1, num_real z1, num_real x2, num_real y2, num_real z2, @@ -219,35 +219,38 @@ class EnergyNonbond_O : public EnergyComponent_O void verifyExcludedAtoms(Matter_sp matter, ScoringFunction_sp score); virtual double evaluateAllComponent( ScoringFunction_sp scorer, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ); + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ); double debugAllComponent( ScoringFunction_sp scorer, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ); + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ); virtual void compareAnalyticalAndNumericalForceAndHessianTermByTerm( ScoringFunction_sp score, NVector_sp pos, + core::T_sp energyScale, core::T_sp activeAtomMask ); - void expandExcludedAtomsToTerms(ScoringFunction_sp score); + void expandExcludedAtomsToTerms(ScoringFunction_sp score, core::T_sp energyScale ); virtual string beyondThresholdInteractionsAsString(); @@ -280,7 +283,11 @@ class EnergyNonbond_O : public EnergyComponent_O {}; }; -EnergyFunction_sp energyFunctionNonbondParameters(ScoringFunction_sp score, double& dielectricConstant, double& dQ1Q2Scale, double& cutoff ); +EnergyFunction_sp energyFunctionNonbondParameters(ScoringFunction_sp score, + core::T_sp energyScale, + double& dielectricConstant, + double& dQ1Q2Scale, + double& cutoff ); }; diff --git a/include/cando/chem/energyOutOfZPlane.h b/include/cando/chem/energyOutOfZPlane.h index 4d6b4b26..fff419a2 100644 --- a/include/cando/chem/energyOutOfZPlane.h +++ b/include/cando/chem/energyOutOfZPlane.h @@ -156,17 +156,18 @@ class EnergyOutOfZPlane_O : public EnergyComponent_O AbstractLargeSquareMatrix_sp m, core::T_sp activeAtomMask ); virtual double evaluateAllComponent( ScoringFunction_sp scorer, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ); + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ); virtual void compareAnalyticalAndNumericalForceAndHessianTermByTerm( NVector_sp pos ); diff --git a/include/cando/chem/energyPeriodicBoundaryConditionsNonbond.h b/include/cando/chem/energyPeriodicBoundaryConditionsNonbond.h index 94239797..3dfd7e02 100644 --- a/include/cando/chem/energyPeriodicBoundaryConditionsNonbond.h +++ b/include/cando/chem/energyPeriodicBoundaryConditionsNonbond.h @@ -66,6 +66,7 @@ class EnergyPeriodicBoundaryConditionsNonbond_O : public EnergyNonbond_O virtual void evaluateUsingExcludedAtoms( ScoringFunction_sp score, NVector_sp pos, + core::T_sp energyScale, bool calcForce, gc::Nilable force, bool calcDiagonalHessian, @@ -76,6 +77,7 @@ class EnergyPeriodicBoundaryConditionsNonbond_O : public EnergyNonbond_O core::T_sp activeAtomMask ); virtual void evaluateTerms( ScoringFunction_sp score, NVector_sp pos, + core::T_sp energyScale, bool calcForce, gc::Nilable force, bool calcDiagonalHessian, diff --git a/include/cando/chem/energyPointToLineRestraint.h b/include/cando/chem/energyPointToLineRestraint.h index 2c534a04..006b227a 100644 --- a/include/cando/chem/energyPointToLineRestraint.h +++ b/include/cando/chem/energyPointToLineRestraint.h @@ -77,6 +77,7 @@ class EnergyPointToLineRestraint_O : public EnergyComponent_O public: virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/energyRigidBodyComponent.h b/include/cando/chem/energyRigidBodyComponent.h index cc187c8d..9543c08e 100644 --- a/include/cando/chem/energyRigidBodyComponent.h +++ b/include/cando/chem/energyRigidBodyComponent.h @@ -56,17 +56,18 @@ class EnergyRigidBodyComponent_O : public EnergyComponent_O DEFAULT_CTOR_DTOR(EnergyRigidBodyComponent_O); virtual double evaluateAllComponent( ScoringFunction_sp scorer, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) = 0; + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) = 0; }; diff --git a/include/cando/chem/energyRigidBodyNonbond.h b/include/cando/chem/energyRigidBodyNonbond.h index b5311de9..22c8ed47 100644 --- a/include/cando/chem/energyRigidBodyNonbond.h +++ b/include/cando/chem/energyRigidBodyNonbond.h @@ -183,6 +183,7 @@ class EnergyRigidBodyNonbond_O : public EnergyRigidBodyComponent_O virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/energyRigidBodyStaple.h b/include/cando/chem/energyRigidBodyStaple.h index e1c6bb53..b057aa34 100644 --- a/include/cando/chem/energyRigidBodyStaple.h +++ b/include/cando/chem/energyRigidBodyStaple.h @@ -140,6 +140,7 @@ class EnergyRigidBodyStaple_O : public EnergyRigidBodyComponent_O core::T_sp activeAtomMask ); virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/energySketchNonbond.h b/include/cando/chem/energySketchNonbond.h index e515a11e..40786586 100644 --- a/include/cando/chem/energySketchNonbond.h +++ b/include/cando/chem/energySketchNonbond.h @@ -129,6 +129,7 @@ class EnergySketchNonbond_O : public EnergyComponent_O virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, @@ -141,6 +142,7 @@ class EnergySketchNonbond_O : public EnergyComponent_O core::T_sp debugInteractions ); double evaluateTerms( NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/energySketchStretch.h b/include/cando/chem/energySketchStretch.h index 72c543dd..c804c1e9 100644 --- a/include/cando/chem/energySketchStretch.h +++ b/include/cando/chem/energySketchStretch.h @@ -166,6 +166,7 @@ class EnergySketchStretch_O : public EnergyComponent_O core::T_sp activeAtomMask ); virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/energyStretch.h b/include/cando/chem/energyStretch.h index e711e04a..e5839eee 100644 --- a/include/cando/chem/energyStretch.h +++ b/include/cando/chem/energyStretch.h @@ -193,6 +193,7 @@ class EnergyStretch_O : public EnergyComponent_O core::T_sp activeAtomMask); virtual double evaluateAllComponent( ScoringFunction_sp scorer, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/minimizer.h b/include/cando/chem/minimizer.h index a8b7ce35..b76ce361 100644 --- a/include/cando/chem/minimizer.h +++ b/include/cando/chem/minimizer.h @@ -171,7 +171,9 @@ struct RestartMinimizer {}; void set_initial_line_search_step(double step); void lineSearchInitialReport( StepReport_sp report, - NVector_sp nvPos, NVector_sp nvDir, + NVector_sp nvPos, + core::T_sp energyScale, + NVector_sp nvDir, NVector_sp nvForce, double xa, double xx, double xb, double fa, double fx, double fb, @@ -181,11 +183,12 @@ struct RestartMinimizer {}; void stepReport( StepReport_sp report, double energy,NVector_sp force,core::T_sp activeAtomMask ); void getPosition(NVector_sp nvResult, NVector_sp nvOrigin, NVector_sp nvDirection, double x, core::T_sp activeAtomMask); - double dTotalEnergy(NVector_sp pos,core::T_sp activeAtomMask); - double dTotalEnergyForce( NVector_sp nvPos, NVector_sp nvForce, core::T_sp activeAtomMask); - double d1DTotalEnergy(double x, core::T_sp activeAtomMask ); - double d1DTotalEnergyForce( double x, double* fx, double* dfx, core::T_sp activeAtomMask ); + double dTotalEnergy(NVector_sp pos,core::T_sp energyScale,core::T_sp activeAtomMask); + double dTotalEnergyForce( NVector_sp nvPos, core::T_sp energyScale, NVector_sp nvForce, core::T_sp activeAtomMask); + double d1DTotalEnergy(double x, core::T_sp energyScale, core::T_sp activeAtomMask ); + double d1DTotalEnergyForce( double x, core::T_sp energyScale, double* fx, double* dfx, core::T_sp activeAtomMask ); void minBracket( NVector_sp nvOrigin, + core::T_sp energyScale, NVector_sp nvDir, double *dPxa, double *dPxb, @@ -195,6 +198,7 @@ struct RestartMinimizer {}; double *dPfc, core::T_sp activeAtomMask ); double dbrent( double ax, double bx, double cx, + core::T_sp energyScale, double tol, double& step, int& energyEvals, int& forceEvals, @@ -203,6 +207,7 @@ struct RestartMinimizer {}; void lineSearch( double *dPxnew, double *dPfnew, NVector_sp nvOrigin, + core::T_sp energyScale, NVector_sp nvDirection, NVector_sp nvForce, NVector_sp nvTemp1, @@ -220,25 +225,32 @@ struct RestartMinimizer {}; /*! Return true on success */ void _steepestDescent( int numSteps, NVector_sp p, + core::T_sp energyScale, double rmsGradientTol, core::T_sp activeAtomMask, core::T_sp callback ); /*! Return true on success */ void _conjugateGradient( int numSteps, NVector_sp p, + core::T_sp energyScale, double rmsGradientTol, core::T_sp activeAtomMask, core::T_sp callback ); void _truncatedNewton( int numSteps, NVector_sp p, + core::T_sp energyScale, double rmsGradientTol, core::T_sp activeAtomMask, core::T_sp callback ); - void _evaluateEnergyAndForceManyTimes( int numSteps, NVector_sp nvPos, + void _evaluateEnergyAndForceManyTimes( int numSteps, + NVector_sp nvPos, + core::T_sp energyScale, core::T_sp activeAtomMask ); - void validateForce(NVector_sp pos, NVector_sp force, + void validateForce(NVector_sp pos, + core::T_sp energyScale, + NVector_sp force, core::T_sp activeAtomMask ); void debugBeginIteration(int i); @@ -260,10 +272,11 @@ struct RestartMinimizer {}; void _truncatedNewtonInnerLoop( int k, NVector_sp xj, + core::T_sp energyScale, SparseLargeSquareMatrix_sp mprecon, SparseLargeSquareMatrix_sp ldlt, NVector_sp force, - double rmsForceMag, + double rmsForceMag, NVector_sp pj, NVector_sp pjNext, NVector_sp rj, @@ -332,11 +345,11 @@ takes a single argument, the NVECTOR position of the atoms.)dx"); void minimizeSteepestDescent(); void minimizeConjugateGradient(); - void resetAndMinimize(core::T_sp activeAtomMask, core::T_sp callback ); - core::T_mv minimize(core::T_sp activeAtomMask, core::T_sp callback); + void resetAndMinimize(core::T_sp energyScale, core::T_sp activeAtomMask, core::T_sp callback ); + core::T_mv minimize(core::T_sp energyScale, core::T_sp activeAtomMask, core::T_sp callback); // If the minimization is aborted the intermediate results can be recovered void writeIntermediateResultsToEnergyFunction(); - void evaluateEnergyAndForceManyTimes(int times,core::T_sp activeAtomMask ); + void evaluateEnergyAndForceManyTimes(core::T_sp energyScale, int times,core::T_sp activeAtomMask ); adapt::QDomNode_sp asXml(); diff --git a/include/cando/chem/rigidBodyEnergyFunction.h b/include/cando/chem/rigidBodyEnergyFunction.h index f991d57e..d09abeeb 100644 --- a/include/cando/chem/rigidBodyEnergyFunction.h +++ b/include/cando/chem/rigidBodyEnergyFunction.h @@ -142,7 +142,7 @@ class RigidBodyEnergyFunction_O : public ScoringFunction_O adapt::QDomNode_sp identifyTermsBeyondThreshold(); // uint countBadVdwInteractions(double scaleSumOfVdwRadii, geom::DisplayList_sp displayIn); - ForceMatchReport_sp checkIfAnalyticalForceMatchesNumericalForce( NVector_sp pos, NVector_sp force, core::T_sp activeAtomMask ); + ForceMatchReport_sp checkIfAnalyticalForceMatchesNumericalForce( NVector_sp pos, core::T_sp energyScale, NVector_sp force, core::T_sp activeAtomMask ); void useDefaultSettings(); @@ -154,6 +154,7 @@ class RigidBodyEnergyFunction_O : public ScoringFunction_O void setOptions( core::List_sp options ); double evaluateAll(NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/include/cando/chem/scoringFunction.h b/include/cando/chem/scoringFunction.h index 7c2580d0..6c553a67 100644 --- a/include/cando/chem/scoringFunction.h +++ b/include/cando/chem/scoringFunction.h @@ -65,6 +65,48 @@ namespace chem { FORWARD(ScoringFunction); }; +#define DefaultChiralRestraintOffset 0.2 +#define DefaultChiralRestraintWeight 100000.0 +#define DefaultAnchorRestraintWeight 10.0 + +namespace chem { + FORWARD(EnergyScale); +class EnergyScale_O : public core::CxxObject_O +{ + LISP_CLASS(chem,ChemPkg,EnergyScale_O,"EnergyScale",core::CxxObject_O); +public: + double _DielectricConstant; + double _NonbondCutoff; + double _ScaleVdw; + double _ScaleElectrostatic; + double _FixedNonbondRestraintWeight; + EnergyScale_O() : _DielectricConstant(1.0), + _NonbondCutoff(16.0), + _ScaleVdw(1.0), + _ScaleElectrostatic(1.0), + _FixedNonbondRestraintWeight(1.0) + { + }; + + static EnergyScale_sp make(); + CL_DEFMETHOD void setVdwScale(double d) { this->_ScaleVdw = d; }; + CL_DEFMETHOD double getVdwScale() {return this->_ScaleVdw; }; + CL_DEFMETHOD void setElectrostaticScale(double d) { this->_ScaleElectrostatic = d; }; + CL_DEFMETHOD double getElectrostaticScale() {return this->_ScaleElectrostatic; }; + CL_DEFMETHOD void setDielectricConstant(double d) { this->_DielectricConstant = d; }; + CL_DEFMETHOD double getDielectricConstant() { return this->_DielectricConstant; }; + CL_DEFMETHOD void setNonbondCutoff(double d) { this->_NonbondCutoff = d; }; + CL_DEFMETHOD double getNonbondCutoff() { return this->_NonbondCutoff; }; + +}; + +extern double energyScaleDielectricConstant(core::T_sp energyScale); +extern double energyScaleVdwScale(core::T_sp energyScale); +extern double energyScaleElectrostaticScale(core::T_sp energyScale); +extern double energyScaleNonbondCutoff(core::T_sp energyScale); + +}; + template <> struct gctools::GCInfo { static bool constexpr NeedsInitialization = false; @@ -115,7 +157,7 @@ class ScoringFunction_O : public core::CxxObject_O NVector_sp makeCoordinates() const; // virtual double evaluate( NVector_sp pos, NVector_sp force, bool calculateForce ) = 0; - virtual ForceMatchReport_sp checkIfAnalyticalForceMatchesNumericalForce( NVector_sp pos, NVector_sp force, core::T_sp activeAtomMask ) = 0; + virtual ForceMatchReport_sp checkIfAnalyticalForceMatchesNumericalForce( NVector_sp pos, core::T_sp energyScale, NVector_sp force, core::T_sp activeAtomMask ) = 0; CL_LISPIFY_NAME("useDefaultSettings"); CL_DEFMETHOD virtual void useDefaultSettings() = 0; @@ -176,14 +218,15 @@ class ScoringFunction_O : public core::CxxObject_O CL_DEFMETHOD double calculateNumericalSecondDerivative(NVector_sp pos, double delta, uint i, uint j, core::T_sp activeAtomMask ); #endif +#if 0 double calculateEnergy(core::T_sp activeAtomMask, core::T_sp debugInteractions ); - - core::T_mv calculateEnergyAndForce(core::T_sp activeAtomMask, core::T_sp debugInteractions ); - + core::T_mv calculateEnergyAndForce(core::T_sp energyScale, core::T_sp activeAtomMask, core::T_sp debugInteractions ); +#endif CL_LISPIFY_NAME("evaluateAll"); - CL_LAMBDA((scoring-function chem:scoring-function) pos &optional component-energy calc-force force calc-diagonal-hessian calc-off-diagonal-hessian hessian hdvec dvec active-atom-mask debug-interactions); + CL_LAMBDA((scoring-function chem:scoring-function) pos &key energy-scale component-energy calc-force force calc-diagonal-hessian calc-off-diagonal-hessian hessian hdvec dvec active-atom-mask debug-interactions); CL_DEFMETHOD virtual double evaluateAll( NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, @@ -194,22 +237,28 @@ class ScoringFunction_O : public core::CxxObject_O gc::Nilable dvec, core::T_sp activeAtomMask, core::T_sp debugInteractions = nil()) = 0; - virtual double evaluateEnergy( NVector_sp pos, core::T_sp activeAtomMask, core::T_sp debugInteractions=nil() ); + virtual double evaluateEnergy( NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + core::T_sp activeAtomMask, + core::T_sp debugInteractions=nil() ); virtual double evaluateEnergyForce( NVector_sp pos, + core::T_sp energyScale, bool calcForce, NVector_sp force, core::T_sp activeAtomMask, core::T_sp debugInteractions = nil() ); virtual double evaluateEnergyForceFullHessian(NVector_sp pos, + core::T_sp energyScale, bool calcForce, NVector_sp force, bool calcDiagonalHessian, bool calcOffDiagonalHessian, AbstractLargeSquareMatrix_sp hessian, core::T_sp activeAtomMask, core::T_sp debugInteractions=nil()); - virtual double evaluateEnergyForceFullHessianForDebugging(core::T_sp activeAtomMask, core::T_sp debugInteractions ); + virtual double evaluateEnergyForceFullHessianForDebugging(core::T_sp energyScale, core::T_sp activeAtomMask, core::T_sp debugInteractions ); - void evaluateFiniteDifferenceForce(NVector_sp pos, NVector_sp force, double delta=0.00001, core::T_sp activeAtomMask=nil() ); + void evaluateFiniteDifferenceForce(NVector_sp pos, core::T_sp energyScale, NVector_sp force, double delta=0.00001, core::T_sp activeAtomMask=nil() ); #if 0 string summarizeBeyondThresholdInteractionsAsString(); @@ -258,6 +307,7 @@ class EnergyComponents_O : public core::CxxObject_O void maybeSetEnergy( core::T_sp componentEnergy, core::T_sp energyComponentName, double energy ); + }; #endif diff --git a/include/cando/chem/sketchFunction.h b/include/cando/chem/sketchFunction.h index a19ecf9e..dbe3bfdf 100644 --- a/include/cando/chem/sketchFunction.h +++ b/include/cando/chem/sketchFunction.h @@ -147,7 +147,7 @@ class SketchFunction_O : public ScoringFunction_O adapt::QDomNode_sp identifyTermsBeyondThreshold(); // uint countBadVdwInteractions(double scaleSumOfVdwRadii, geom::DisplayList_sp displayIn); - ForceMatchReport_sp checkIfAnalyticalForceMatchesNumericalForce( NVector_sp pos, NVector_sp force, core::T_sp activeAtomMask ); + ForceMatchReport_sp checkIfAnalyticalForceMatchesNumericalForce( NVector_sp pos, core::T_sp energyScale,NVector_sp force, core::T_sp activeAtomMask ); CL_LISPIFY_NAME("getGraph"); CL_DEFMETHOD core::T_sp getGraph() { return this->_Graph;}; @@ -188,9 +188,10 @@ class SketchFunction_O : public ScoringFunction_O void dumpTerms(core::HashTable_sp atomTypes); CL_DEFMETHOD core::T_sp getMessage() { return this->_Message;}; - double calculateNumericalDerivative(NVector_sp pos, double delta, uint i, core::T_sp activeAtomMask ); - double calculateNumericalSecondDerivative(NVector_sp pos, double delta, uint i, uint j, core::T_sp activeAtomMask ); + double calculateNumericalDerivative(NVector_sp pos, core::T_sp energyScale, double delta, uint i, core::T_sp activeAtomMask ); + double calculateNumericalSecondDerivative(NVector_sp pos, core::T_sp energyScale, double delta, uint i, uint j, core::T_sp activeAtomMask ); double evaluateAll(NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, @@ -210,8 +211,8 @@ class SketchFunction_O : public ScoringFunction_O // adapt::QDomNode_sp accumulateTermsBeyondThresholdAsXml(); uint countTermsBeyondThreshold(); - void evaluateNumericalForce(NVector_sp pos, NVector_sp numForce, double delta, core::T_sp activeAtomMask ); - void evaluateNumericalHessian(NVector_sp pos, AbstractLargeSquareMatrix_sp numHessian, bool calcOffDiagonalElements, double delta, core::T_sp activeAtomMask ); + void evaluateNumericalForce(NVector_sp pos, core::T_sp energyScale, NVector_sp numForce, double delta, core::T_sp activeAtomMask ); + void evaluateNumericalHessian(NVector_sp pos, core::T_sp energyScale, AbstractLargeSquareMatrix_sp numHessian, bool calcOffDiagonalElements, double delta, core::T_sp activeAtomMask ); string debugLogAsString(); diff --git a/include/cando/kinematics/bondedJoint.h b/include/cando/kinematics/bondedJoint.h index b488997c..4c3d1ae6 100644 --- a/include/cando/kinematics/bondedJoint.h +++ b/include/cando/kinematics/bondedJoint.h @@ -112,7 +112,7 @@ FORWARD(BondedJoint); /*! Get the stub and update the XYZ coordinate */ void updateXyzCoord(chem::NVector_sp coords); - + virtual Stub getInputStub(chem::NVector_sp coords) const; virtual bool definedp() const; diff --git a/include/cando/kinematics/jumpJoint.h b/include/cando/kinematics/jumpJoint.h index 7d2ff6c4..82c19c0a 100644 --- a/include/cando/kinematics/jumpJoint.h +++ b/include/cando/kinematics/jumpJoint.h @@ -75,6 +75,8 @@ namespace kinematics virtual core::Symbol_sp typeSymbol() const; Stub getInputStub(chem::NVector_sp coords) const; + + virtual bool definedp() const; virtual void _updateInternalCoord(chem::NVector_sp coords); diff --git a/src/chem/cluster.cc b/src/chem/cluster.cc index b2d11b71..20357990 100644 --- a/src/chem/cluster.cc +++ b/src/chem/cluster.cc @@ -155,7 +155,6 @@ void Kmeans_O::Cluster(core::SimpleVector_sp centers, Clusters clusters) /* Calculate the center of each cluster. * Return the number of empty clusters. */ -__attribute__((optnone)) CL_DEFMETHOD int Kmeans_O::Center(core::SimpleVector_sp centers, Clusters clusters) { int emptyClusters = 0; diff --git a/src/chem/energyAnchorRestraint.cc b/src/chem/energyAnchorRestraint.cc index e38cb59b..c74e170b 100644 --- a/src/chem/energyAnchorRestraint.cc +++ b/src/chem/energyAnchorRestraint.cc @@ -252,17 +252,18 @@ void EnergyAnchorRestraint_O::setupHessianPreconditioner(NVector_sp nvPosition, double EnergyAnchorRestraint_O::evaluateAllComponent( ScoringFunction_sp score, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { MAYBE_SETUP_ACTIVE_ATOM_MASK(); MAYBE_SETUP_DEBUG_INTERACTIONS(debugInteractions.notnilp()); diff --git a/src/chem/energyAngle.cc b/src/chem/energyAngle.cc index b7af8a39..b5370526 100644 --- a/src/chem/energyAngle.cc +++ b/src/chem/energyAngle.cc @@ -439,6 +439,7 @@ bool calcOffDiagonalHessian = true; double EnergyAngle_O::evaluateAllComponent( ScoringFunction_sp score, chem::NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/src/chem/energyChiralRestraint.cc b/src/chem/energyChiralRestraint.cc index 8de1a965..3c4e28ba 100644 --- a/src/chem/energyChiralRestraint.cc +++ b/src/chem/energyChiralRestraint.cc @@ -246,17 +246,18 @@ bool calcOffDiagonalHessian = true; double EnergyChiralRestraint_O::evaluateAllComponent( ScoringFunction_sp score, - chem::NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) + chem::NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { MAYBE_SETUP_ACTIVE_ATOM_MASK(); MAYBE_SETUP_DEBUG_INTERACTIONS(debugInteractions.notnilp()); diff --git a/src/chem/energyComponent.cc b/src/chem/energyComponent.cc index 7fc9bcb6..e95be1fc 100644 --- a/src/chem/energyComponent.cc +++ b/src/chem/energyComponent.cc @@ -61,6 +61,14 @@ CL_DEFMETHOD string EnergyComponent_O::debugLogAsString() #endif } +void EnergyComponent_O::fields(core::Record_sp node) +{ + node->field( INTERN_(kw,enabled), this->_Enabled ); + node->field( INTERN_(kw,scale), this->_Scale); + this->Base::fields(node); +} + + string EnergyComponent_O::enabledAsString() { stringstream ss; @@ -83,18 +91,20 @@ string EnergyComponent_O::enabledAsString() CL_DOCSTRING(R"dx(Evaluate the energy of a component)dx"); -CL_LAMBDA(energy-function component pos &optional active-atom-mask debug-interactions); +CL_LAMBDA(energy-function component pos &key energy-scale active-atom-mask debug-interactions); DOCGROUP(cando); CL_DEFUN num_real chem__energy_component_evaluate_energy(EnergyFunction_sp energy_function, EnergyComponent_sp component, chem::NVector_sp pos, + core::T_sp energyScale, core::T_sp activeAtomMask, core::T_sp debugInteractions ) { num_real val = component->evaluateAllComponent(energy_function, pos, + energyScale, nil(), false,nil(), false,false, @@ -107,18 +117,20 @@ num_real chem__energy_component_evaluate_energy(EnergyFunction_sp energy_functio }; CL_DOCSTRING(R"dx(Evaluate the energy and force of a component)dx"); -CL_LAMBDA(energy-function component pos force &optional active-atom-mask debug-interactions); +CL_LAMBDA(energy-function component pos &key energy-scale force active-atom-mask debug-interactions); DOCGROUP(cando); CL_DEFUN num_real chem__energy_component_evaluate_energy_force(EnergyFunction_sp energy_function, EnergyComponent_sp component, NVector_sp pos, + core::T_sp energyScale, NVector_sp force, core::T_sp activeAtomMask, core::T_sp debugInteractions ) { num_real val = component->evaluateAllComponent(energy_function, pos, + energyScale, nil(), true,force, false,false, diff --git a/src/chem/energyDihedral.cc b/src/chem/energyDihedral.cc index f5ff074a..550bae8a 100644 --- a/src/chem/energyDihedral.cc +++ b/src/chem/energyDihedral.cc @@ -1732,17 +1732,18 @@ double EnergyDihedral_O::evaluateAllComponentSimd2( #endif // !_TARGET_OS_DARWIN double EnergyDihedral_O::evaluateAllComponent(ScoringFunction_sp score, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { num_real energy = 0.0; this->_Evaluations++; diff --git a/src/chem/energyDihedralRestraint.cc b/src/chem/energyDihedralRestraint.cc index 0d68ee82..3911c22c 100644 --- a/src/chem/energyDihedralRestraint.cc +++ b/src/chem/energyDihedralRestraint.cc @@ -251,17 +251,18 @@ bool calcOffDiagonalHessian = true; double EnergyDihedralRestraint_O::evaluateAllComponent( ScoringFunction_sp score, - chem::NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) + chem::NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { MAYBE_SETUP_ACTIVE_ATOM_MASK(); MAYBE_SETUP_DEBUG_INTERACTIONS(debugInteractions.notnilp()); diff --git a/src/chem/energyFixedNonbond.cc b/src/chem/energyFixedNonbond.cc index a8493846..865ff09d 100644 --- a/src/chem/energyFixedNonbond.cc +++ b/src/chem/energyFixedNonbond.cc @@ -211,14 +211,17 @@ SYMBOL_EXPORT_SC_(ChemPkg, energyElectrostatic); SYMBOL_EXPORT_SC_(ChemPkg, energyVdw14); SYMBOL_EXPORT_SC_(ChemPkg, energyElectrostatic14); -double EnergyFixedNonbondRestraint_O::evaluateAllComponent(ScoringFunction_sp score, NVector_sp pos, core::T_sp componentEnergy, +double EnergyFixedNonbondRestraint_O::evaluateAllComponent(ScoringFunction_sp score, + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, bool calcForce, gc::Nilable force, bool calcDiagonalHessian, bool calcOffDiagonalHessian, gc::Nilable hessian, gc::Nilable hdvec, gc::Nilable dvec, core::T_sp activeAtomMask, core::T_sp debugInteractions) { - EnergyFunction_sp energyFunction = gc::As(score); - double dielectricConstant = energyFunction->energyScale()->getDielectricConstant(); + + double dielectricConstant = energyScaleDielectricConstant(energyScale); double amber_charge_conversion_18dot2223 = core::Number_O::as_double_float(gc::As(_sym_STARamber_charge_conversion_18_DOT_2223STAR->symbolValue())); double dQ1Q2Scale = amber_charge_conversion_18dot2223 * amber_charge_conversion_18dot2223; diff --git a/src/chem/energyFunction.cc b/src/chem/energyFunction.cc index 2c9db9cb..947b768f 100644 --- a/src/chem/energyFunction.cc +++ b/src/chem/energyFunction.cc @@ -172,17 +172,12 @@ evaluates nothing) and (lambda (component-class)) returns a (lambda (&rest atoms that will be called for each interaction and the function returns T or NIL if each interaction should be added to the energy function. : assign-types - T (default) assign atom types as part of generating the energy function. [I don't know what will happen if assing-types is NIL.)doc"); -CL_LAMBDA(&key matter disable-components enable-components energy-scale (use-excluded-atoms t) (keep-interaction-factory t) (assign-types t)); +CL_LAMBDA(&key matter disable-components enable-components (use-excluded-atoms t) (keep-interaction-factory t) (assign-types t)); CL_LISPIFY_NAME(make_energy_function); CL_DEF_CLASS_METHOD EnergyFunction_sp EnergyFunction_O::make(core::T_sp matter, core::T_sp disableComponents, core::List_sp enableComponents, - core::T_sp energyScale, bool useExcludedAtoms, core::T_sp keepInteractionFactory, bool assign_types ) + bool useExcludedAtoms, core::T_sp keepInteractionFactory, bool assign_types ) { - if (energyScale.unboundp() || energyScale.nilp()) { - energyScale = gctools::GC::allocate_with_default_constructor(); - }; auto me = gctools::GC::allocate_with_default_constructor(); - me->setEnergyScale(energyScale); - if ( matter.notnilp() ) me->defineForMatter(gc::As(matter),useExcludedAtoms,keepInteractionFactory,assign_types); // // Disable and then enable components // @@ -238,22 +233,11 @@ CL_DEF_CLASS_METHOD EnergyFunction_sp EnergyFunction_O::make(core::T_sp matter, } } } + if ( matter.notnilp() ) me->defineForMatter(gc::As(matter),useExcludedAtoms,keepInteractionFactory,assign_types); return me; }; -CL_DEFMETHOD EnergyScale_sp EnergyFunction_O::energyScale() { - return gc::As(this->_EnergyScale); -} - -void EnergyFunction_O::setEnergyScale(core::T_sp energyScale) { - if (!gc::IsA(energyScale)) { - TYPE_ERROR(energyScale,_sym_EnergyScale_O); - } - this->_EnergyScale = energyScale; -} - - CL_DOCSTRING("Return (values enabled-component-class-names disabled-component-class-names)"); CL_DEFMETHOD core::T_mv EnergyFunction_O::enabledDisabled() const { @@ -352,7 +336,6 @@ void EnergyFunction_O::useDefaultSettings() void EnergyFunction_O::fields(core::Record_sp node) { - node->field_if_not_unbound(INTERN_(kw,EnergyScale),this->_EnergyScale); node->field_if_not_unbound(INTERN_(kw,AtomTable),this->_AtomTable); node->field_if_not_unbound(INTERN_(kw,Stretch),this->_Stretch); #if USE_ALL_ENERGY_COMPONENTS @@ -578,6 +561,7 @@ uint EnergyFunction_O::countTermsBeyondThreshold() // double EnergyFunction_O::evaluateAll( NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, @@ -613,6 +597,7 @@ double EnergyFunction_O::evaluateAll( NVector_sp pos, if (this->_Stretch->isEnabled()) { totalEnergy += this->_Stretch->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, @@ -623,6 +608,7 @@ double EnergyFunction_O::evaluateAll( NVector_sp pos, if (this->_Angle->isEnabled()) { totalEnergy += this->_Angle->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, @@ -632,7 +618,9 @@ double EnergyFunction_O::evaluateAll( NVector_sp pos, } if(this->_Dihedral->isEnabled()) { totalEnergy += this->_Dihedral->evaluateAllComponent( this->asSmartPtr(), - pos, componentEnergy, + pos, + energyScale, + componentEnergy, calcForce, force, calcDiagonalHessian, calcOffDiagonalHessian, @@ -642,6 +630,7 @@ double EnergyFunction_O::evaluateAll( NVector_sp pos, if(this->_Nonbond->isEnabled()) { totalEnergy += this->_Nonbond->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, calcOffDiagonalHessian, hessian, hdvec, dvec, activeAtomMask, @@ -650,6 +639,7 @@ double EnergyFunction_O::evaluateAll( NVector_sp pos, if(this->_DihedralRestraint.boundp() && this->_DihedralRestraint->isEnabled()) { totalEnergy += this->_DihedralRestraint->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, @@ -663,6 +653,7 @@ double EnergyFunction_O::evaluateAll( NVector_sp pos, if(this->_ChiralRestraint->isEnabled()) { totalEnergy += this->_ChiralRestraint->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, @@ -677,6 +668,7 @@ double EnergyFunction_O::evaluateAll( NVector_sp pos, if(this->_AnchorRestraint->isEnabled()) { totalEnergy += this->_AnchorRestraint->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, @@ -691,6 +683,7 @@ double EnergyFunction_O::evaluateAll( NVector_sp pos, if(this->_FixedNonbondRestraint->isEnabled()) { totalEnergy += this->_FixedNonbondRestraint->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, @@ -707,6 +700,7 @@ double EnergyFunction_O::evaluateAll( NVector_sp pos, if (component->isEnabled()) { totalEnergy+= component->evaluateAllComponent(this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, @@ -846,33 +840,33 @@ string EnergyFunction_O::energyTermsEnabled() #define DELTA 0.00000001 -double EnergyFunction_O::calculateNumericalDerivative(NVector_sp pos, double delta, uint i, core::T_sp activeAtomMask ) +double EnergyFunction_O::calculateNumericalDerivative(NVector_sp pos, core::T_sp energyScale, double delta, uint i, core::T_sp activeAtomMask ) { double x, ylow, yhigh, fval; double deltaDiv2 = delta/2.0; x = pos->element(i); pos->setElement(i,x-deltaDiv2); - ylow = this->evaluateEnergy(pos,activeAtomMask); + ylow = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x+deltaDiv2); - yhigh = this->evaluateEnergy(pos,activeAtomMask); + yhigh = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x); fval = (yhigh-ylow)/delta; return fval; } -double EnergyFunction_O::calculateNumericalSecondDerivative(NVector_sp pos, double delta, uint i, uint j, core::T_sp activeAtomMask ) +double EnergyFunction_O::calculateNumericalSecondDerivative(NVector_sp pos, core::T_sp energyScale, double delta, uint i, uint j, core::T_sp activeAtomMask ) { double x, fxmh, fx, fxph, f2; double y, fpipj, fpimj, fmipj, fmimj, fp, fm; if ( i==j ) { x = pos->element(i); pos->setElement(i,x-delta); - fxmh = this->evaluateEnergy(pos,activeAtomMask); + fxmh = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x+delta); - fxph = this->evaluateEnergy(pos,activeAtomMask); + fxph = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x); - fx = this->evaluateEnergy(pos,activeAtomMask); + fx = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); f2 = (fxph+fxmh-2.0*(fx))/(delta*delta); } else { double deltaDiv2 = delta/2.0; @@ -880,16 +874,16 @@ double EnergyFunction_O::calculateNumericalSecondDerivative(NVector_sp pos, doub y = pos->element(j); pos->setElement(i,x+deltaDiv2); pos->setElement(j,y+deltaDiv2); - fpipj = this->evaluateEnergy(pos,activeAtomMask); + fpipj = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x+deltaDiv2); pos->setElement(j,y-deltaDiv2); - fpimj = this->evaluateEnergy(pos,activeAtomMask); + fpimj = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x-deltaDiv2); pos->setElement(j,y+deltaDiv2); - fmipj = this->evaluateEnergy(pos,activeAtomMask); + fmipj = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x-deltaDiv2); pos->setElement(j,y-deltaDiv2); - fmimj = this->evaluateEnergy(pos,activeAtomMask); + fmimj = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x); pos->setElement(j,y); LOG("fpipj = {}" , fpipj ); @@ -911,13 +905,13 @@ double EnergyFunction_O::calculateNumericalSecondDerivative(NVector_sp pos, doub /*! Calculate the force numerically */ -void EnergyFunction_O::evaluateNumericalForce(NVector_sp pos, NVector_sp numForce, double delta, core::T_sp activeAtomMask ) +void EnergyFunction_O::evaluateNumericalForce(NVector_sp pos, core::T_sp energyScale, NVector_sp numForce, double delta, core::T_sp activeAtomMask ) { double fval; uint i; for (i=0; isize(); i++ ) { - fval = -this->calculateNumericalDerivative(pos,delta,i,activeAtomMask); + fval = -this->calculateNumericalDerivative(pos,energyScale,delta,i,activeAtomMask); numForce->setElement(i,fval); } } @@ -925,7 +919,7 @@ void EnergyFunction_O::evaluateNumericalForce(NVector_sp pos, NVector_sp numForc /*! Calculate the hessian numerically */ -void EnergyFunction_O::evaluateNumericalHessian(NVector_sp pos, AbstractLargeSquareMatrix_sp hessian, bool calcOffDiagonal, double delta, core::T_sp activeAtomMask ) +void EnergyFunction_O::evaluateNumericalHessian(NVector_sp pos, core::T_sp energyScale, AbstractLargeSquareMatrix_sp hessian, bool calcOffDiagonal, double delta, core::T_sp activeAtomMask ) { double fval; uint c, r; @@ -935,14 +929,14 @@ void EnergyFunction_O::evaluateNumericalHessian(NVector_sp pos, AbstractLargeSqu } hessian->zero(); for ( c=0; csize(); c++ ) { - fval = this->calculateNumericalSecondDerivative(pos,delta,c,c,activeAtomMask); + fval = this->calculateNumericalSecondDerivative(pos,energyScale,delta,c,c,activeAtomMask); hessian->setElement(c,c,fval); } if ( !calcOffDiagonal ) return; for ( c=0; csize(); c++ ) { for ( r=0; rsize(); r++ ) { if ( c!=r) { - fval = this->calculateNumericalSecondDerivative(pos,delta,c,r,activeAtomMask); + fval = this->calculateNumericalSecondDerivative(pos,energyScale,delta,c,r,activeAtomMask); hessian->setElement(c,r,fval); } } @@ -957,7 +951,7 @@ void EnergyFunction_O::evaluateNumericalHessian(NVector_sp pos, AbstractLargeSqu * If there is a mis-match then dump the EnergyFunction into the result. * */ -ForceMatchReport_sp EnergyFunction_O::checkIfAnalyticalForceMatchesNumericalForce(NVector_sp pos, NVector_sp analyticalForce, core::T_sp activeAtomMask ) +ForceMatchReport_sp EnergyFunction_O::checkIfAnalyticalForceMatchesNumericalForce(NVector_sp pos, core::T_sp energyScale, NVector_sp analyticalForce, core::T_sp activeAtomMask ) { ForceMatchReport_sp report; NVector_sp numForce, tempForce; @@ -968,13 +962,13 @@ ForceMatchReport_sp EnergyFunction_O::checkIfAnalyticalForceMatchesNumericalForc report = ForceMatchReport_O::create(); numForce = NVector_O::create(pos->size()); - this->evaluateNumericalForce(pos,numForce,DELTA,activeAtomMask); + this->evaluateNumericalForce(pos,energyScale,numForce,DELTA,activeAtomMask); dot = dotProductWithActiveAtomMask(numForce,analyticalForce,nil()); numericalMag = magnitudeWithActiveAtomMask(numForce,nil()); analyticalMag = magnitudeWithActiveAtomMask(analyticalForce,nil()); tempForce = NVector_O::create(pos->size()); // Evaluate the force at pos again - this->evaluateEnergyForce(pos,true,tempForce,activeAtomMask); + this->evaluateEnergyForce(pos,energyScale,true,tempForce,activeAtomMask); avg = (analyticalMag+numericalMag)/2.0; if ( analyticalMag < VERYSMALL && numericalMag < VERYSMALL ) { result.str(""); @@ -1140,7 +1134,7 @@ int EnergyFunction_O::_applyRestraints(core::T_sp nonbondDb, core::Iterator_sp r iterm.term.xa = anchorPos.getX(); iterm.term.ya = anchorPos.getY(); iterm.term.za = anchorPos.getZ(); - iterm.term.ka = this->energyScale()->_AnchorRestraintWeight; + iterm.term.ka = DefaultAnchorRestraintWeight; iterm.term.I1 = ea1->coordinateIndexTimes3(); this->_AnchorRestraint->addTerm(iterm); ++terms; @@ -1274,7 +1268,13 @@ CL_DEFMETHOD void EnergyFunction_O::defineForMatter(Matter_sp matter, bool useEx // if (chem__verbose(0)) core::clasp_write_string("Searching for rings.\n"); - core::T_sp rings = RingFinder_O::identifyRings(matter); + // If *current-rings* is bound then reuse it, otherwise, calculate rings and bind those + core::T_sp rings = unbound(); + if (_sym_STARcurrent_ringsSTAR->boundP()) { + rings = _sym_STARcurrent_ringsSTAR->symbolValue(); + } else { + rings = RingFinder_O::identifyRings(matter); + } core::DynamicScopeManager ring_scope(_sym_STARcurrent_ringsSTAR, rings ); // @@ -1893,8 +1893,8 @@ CL_DEFMETHOD void EnergyFunction_O::generateRestraintEnergyFunctionTables(Matter ichiral.term.I2 = ea2->coordinateIndexTimes3(); ichiral.term.I3 = eaCenter->coordinateIndexTimes3(); ichiral.term.I4 = ea3->coordinateIndexTimes3(); - ichiral.term.K = this->energyScale()->_ChiralRestraintWeight * side; - ichiral.term.CO = this->energyScale()->_ChiralRestraintOffset; + ichiral.term.K = DefaultChiralRestraintWeight * side; + ichiral.term.CO = DefaultChiralRestraintOffset; this->_ChiralRestraint->addTerm(ichiral); // Now apply it to the other atom // on the chiral center, just flip the sign @@ -1905,8 +1905,8 @@ CL_DEFMETHOD void EnergyFunction_O::generateRestraintEnergyFunctionTables(Matter if ( !skipInteraction( keepInteraction, ichiral._Atom1, ichiral._Atom2, ichiral._Atom3, ichiral._Atom4 ) ) { ichiral.term.I4 = ea4->coordinateIndexTimes3(); // flip the sign of the chiral restraint - ichiral.term.K = this->energyScale()->_ChiralRestraintWeight * side * -1.0; - ichiral.term.CO = this->energyScale()->_ChiralRestraintOffset; + ichiral.term.K = DefaultChiralRestraintWeight * side * -1.0; + ichiral.term.CO = DefaultChiralRestraintOffset; this->_ChiralRestraint->addTerm(ichiral); } // To try and increase the number of molecules that @@ -1925,8 +1925,8 @@ CL_DEFMETHOD void EnergyFunction_O::generateRestraintEnergyFunctionTables(Matter ichiral.term.I2 = ea4->coordinateIndexTimes3(); ichiral.term.I3 = eaCenter->coordinateIndexTimes3(); ichiral.term.I4 = ea3->coordinateIndexTimes3(); - ichiral.term.K = this->energyScale()->_ChiralRestraintWeight * side; - ichiral.term.CO = this->energyScale()->_ChiralRestraintOffset; + ichiral.term.K = DefaultChiralRestraintWeight * side; + ichiral.term.CO = DefaultChiralRestraintOffset; this->_ChiralRestraint->addTerm(ichiral); } // Now apply it to the other atom @@ -1936,8 +1936,8 @@ CL_DEFMETHOD void EnergyFunction_O::generateRestraintEnergyFunctionTables(Matter if ( !skipInteraction( keepInteraction, ichiral._Atom1, ichiral._Atom2, ichiral._Atom3, ichiral._Atom4 ) ) { ichiral.term.I4 = ea1->coordinateIndexTimes3(); // flip the sign of the chiral restraint - ichiral.term.K = this->energyScale()->_ChiralRestraintWeight * side * -1.0; - ichiral.term.CO = this->energyScale()->_ChiralRestraintOffset; + ichiral.term.K = DefaultChiralRestraintWeight * side * -1.0; + ichiral.term.CO = DefaultChiralRestraintOffset; this->_ChiralRestraint->addTerm(ichiral); } } else { diff --git a/src/chem/energyNonbond.cc b/src/chem/energyNonbond.cc index 06564e2a..53d1de69 100644 --- a/src/chem/energyNonbond.cc +++ b/src/chem/energyNonbond.cc @@ -99,11 +99,14 @@ namespace chem { #define BAIL_OUT_IF_CUTOFF(deltaSquared) // #define BAIL_OUT_IF_CUTOFF(deltaSquared) if (deltaSquared>CUTOFF_SQUARED) goto SKIP_term; -EnergyFunction_sp energyFunctionNonbondParameters(ScoringFunction_sp score, double &dielectricConstant, double &dQ1Q2Scale, +EnergyFunction_sp energyFunctionNonbondParameters(ScoringFunction_sp score, + core::T_sp energyScale, + double &dielectricConstant, + double &dQ1Q2Scale, double &cutoff) { - EnergyFunction_sp energyFunction = gc::As(score); - cutoff = energyFunction->energyScale()->getNonbondCutoff(); - dielectricConstant = energyFunction->energyScale()->getDielectricConstant(); + auto energyFunction = gc::As(score); + cutoff = energyScaleNonbondCutoff(energyScale); + dielectricConstant = energyScaleDielectricConstant(energyScale); double amber_charge_conversion_18dot2223 = core::Number_O::as_double_float(gc::As(_sym_STARamber_charge_conversion_18_DOT_2223STAR->symbolValue())); dQ1Q2Scale = amber_charge_conversion_18dot2223 * amber_charge_conversion_18dot2223; @@ -156,7 +159,7 @@ core::List_sp EnergyNonbond::encode() const { void EnergyNonbond::decode(core::List_sp alist) { SIMPLE_ERROR("Implement decode of EnergyNonbond"); } -#define ENERGY_FUNCTION score, I1, I2, activeAtomMask +#define ENERGY_FUNCTION score, energyScale, I1, I2, activeAtomMask #define NONBOND_DEBUG_INTERACTIONS(I1, I2) \ if (doDebugInteractions) { \ @@ -166,45 +169,45 @@ void EnergyNonbond::decode(core::List_sp alist) { SIMPLE_ERROR("Implement decode core::make_fixnum(I1), core::make_fixnum(I2)); \ } -double _evaluateEnergyOnly_Nonbond(ScoringFunction_sp score, int I1, int I2, core::T_sp activeAtomMask, num_real x1, num_real y1, +double _evaluateEnergyOnly_Nonbond(ScoringFunction_sp score, core::T_sp energyScale, int I1, int I2, core::T_sp activeAtomMask, num_real x1, num_real y1, num_real z1, num_real x2, num_real y2, num_real z2, num_real dA, num_real dC, num_real dQ1Q2); -void calculateFiniteDifferenceForces(ScoringFunction_sp score, core::T_sp activeAtomMask, int I1, int I2, double x1, double y1, +void calculateFiniteDifferenceForces(ScoringFunction_sp score, core::T_sp energyScale, core::T_sp activeAtomMask, int I1, int I2, double x1, double y1, double z1, double x2, double y2, double z2, double dA, double dC, double dQ1Q2, double &fdfx1, double &fdfy1, double &fdfz1, double &fdfx2, double &fdfy2, double &fdfz2) { double delta = 0.00001; double deltaTimes2 = delta * 2.0; double ehigh; double elow; - ehigh = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1 + delta, y1, z1, x2, y2, z2, dA, dC, dQ1Q2); - elow = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1 - delta, y1, z1, x2, y2, z2, dA, dC, dQ1Q2); + ehigh = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1 + delta, y1, z1, x2, y2, z2, dA, dC, dQ1Q2); + elow = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1 - delta, y1, z1, x2, y2, z2, dA, dC, dQ1Q2); fdfx1 = -(ehigh - elow) / deltaTimes2; - ehigh = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1, y1 + delta, z1, x2, y2, z2, dA, dC, dQ1Q2); - elow = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1, y1 - delta, z1, x2, y2, z2, dA, dC, dQ1Q2); + ehigh = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1, y1 + delta, z1, x2, y2, z2, dA, dC, dQ1Q2); + elow = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1, y1 - delta, z1, x2, y2, z2, dA, dC, dQ1Q2); fdfy1 = -(ehigh - elow) / deltaTimes2; - ehigh = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1, y1, z1 + delta, x2, y2, z2, dA, dC, dQ1Q2); - elow = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1, y1, z1 - delta, x2, y2, z2, dA, dC, dQ1Q2); + ehigh = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1, y1, z1 + delta, x2, y2, z2, dA, dC, dQ1Q2); + elow = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1, y1, z1 - delta, x2, y2, z2, dA, dC, dQ1Q2); fdfz1 = -(ehigh - elow) / deltaTimes2; - ehigh = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1, y1, z1, x2 + delta, y2, z2, dA, dC, dQ1Q2); - elow = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1, y1, z1, x2 - delta, y2, z2, dA, dC, dQ1Q2); + ehigh = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1, y1, z1, x2 + delta, y2, z2, dA, dC, dQ1Q2); + elow = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1, y1, z1, x2 - delta, y2, z2, dA, dC, dQ1Q2); fdfx2 = -(ehigh - elow) / deltaTimes2; - ehigh = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1, y1, z1, x2, y2 + delta, z2, dA, dC, dQ1Q2); - elow = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1, y1, z1, x2, y2 - delta, z2, dA, dC, dQ1Q2); + ehigh = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1, y1, z1, x2, y2 + delta, z2, dA, dC, dQ1Q2); + elow = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1, y1, z1, x2, y2 - delta, z2, dA, dC, dQ1Q2); fdfy2 = -(ehigh - elow) / deltaTimes2; - ehigh = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1, y1, z1, x2, y2, z2 + delta, dA, dC, dQ1Q2); - elow = _evaluateEnergyOnly_Nonbond(score, I1, I2, activeAtomMask, x1, y1, z1, x2, y2, z2 - delta, dA, dC, dQ1Q2); + ehigh = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1, y1, z1, x2, y2, z2 + delta, dA, dC, dQ1Q2); + elow = _evaluateEnergyOnly_Nonbond(score, energyScale, I1, I2, activeAtomMask, x1, y1, z1, x2, y2, z2 - delta, dA, dC, dQ1Q2); fdfz2 = -(ehigh - elow) / deltaTimes2; } struct NoFiniteDifference { - static void maybeTestFiniteDifference(ScoringFunction_sp score, int I1, int I2, core::T_sp activeAtomMask, double x1, double y1, + static void maybeTestFiniteDifference(ScoringFunction_sp score, core::T_sp energyScale, int I1, int I2, core::T_sp activeAtomMask, double x1, double y1, double z1, double x2, double y2, double z2, double dA, double dC, double dQ1Q2, double fx1, double fy1, double fz1, double fx2, double fy2, double fz2, int index, size_t &fails, bool debugForce) {} }; struct DebugFiniteDifference { - static void maybeTestFiniteDifference(ScoringFunction_sp score, int I1, int I2, core::T_sp activeAtomMask, double x1, double y1, + static void maybeTestFiniteDifference(ScoringFunction_sp score, core::T_sp energyScale, int I1, int I2, core::T_sp activeAtomMask, double x1, double y1, double z1, double x2, double y2, double z2, double dA, double dC, double dQ1Q2, double fx1, double fy1, double fz1, double fx2, double fy2, double fz2, int index, size_t &fails, bool debugForce) { @@ -212,7 +215,7 @@ struct DebugFiniteDifference { #define TEST_DIAGONAL_HESSIAN(x, y, o1, o2, o3, d, i) #undef TEST_OFF_DIAGONAL_HESSIAN #define TEST_OFF_DIAGONAL_HESSIAN(x, y, o1, o2, o3, o4, d, i) -#define ENERGY_FUNCTION score, I1, I2, activeAtomMask +#define ENERGY_FUNCTION score, energyScale, I1, I2, activeAtomMask #include #undef ENERGY_FUNCTION if (debugForce) { @@ -222,7 +225,7 @@ struct DebugFiniteDifference { double fdfx2; double fdfy2; double fdfz2; - calculateFiniteDifferenceForces(score, activeAtomMask, I1, I2, x1, y1, z1, x2, y2, z2, dA, dC, dQ1Q2, fdfx1, fdfy1, fdfz1, + calculateFiniteDifferenceForces(score, energyScale, activeAtomMask, I1, I2, x1, y1, z1, x2, y2, z2, dA, dC, dQ1Q2, fdfx1, fdfy1, fdfz1, fdfx2, fdfy2, fdfz2); core::lisp_write(fmt::format("({0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11} {12} {13} {14} {15} {16} {17} {18} {19} " "{20} {21} {22}) ; {23}:{24}:{25} \n", @@ -245,8 +248,8 @@ inline double calculate_dQ1Q2(double electrostaticScale, double electrostaticMod It is a template function so that template arguments can inline or elide testing code. */ template -//__attribute__((optnone)) double template_evaluateUsingExcludedAtoms(EnergyNonbond_O *mthis, ScoringFunction_sp score, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, bool calcDiagonalHessian, bool calcOffDiagonalHessian, gc::Nilable hessian, gc::Nilable hdvec, @@ -255,7 +258,7 @@ double template_evaluateUsingExcludedAtoms(EnergyNonbond_O *mthis, ScoringFuncti double dielectricConstant; double dQ1Q2Scale; double cutoff; - EnergyFunction_sp energyFunction = energyFunctionNonbondParameters(score, dielectricConstant, dQ1Q2Scale, cutoff); + EnergyFunction_sp energyFunction = energyFunctionNonbondParameters(score, energyScale, dielectricConstant, dQ1Q2Scale, cutoff); #define CUTOFF_SQUARED (cutoff * cutoff) MAYBE_SETUP_ACTIVE_ATOM_MASK(); MAYBE_SETUP_DEBUG_INTERACTIONS(debugInteractions.notnilp()); @@ -272,7 +275,7 @@ double template_evaluateUsingExcludedAtoms(EnergyNonbond_O *mthis, ScoringFuncti bool hasHessian = hessian.notnilp(); bool hasHdAndD = (hdvec.notnilp()) && (dvec.notnilp()); double energyElectrostatic = 0.0; - double energyVdw = 0.0; + KahanSummation energyVdw; #define NONBOND_CALC_FORCE #define NONBOND_CALC_DIAGONAL_HESSIAN #define NONBOND_CALC_OFF_DIAGONAL_HESSIAN @@ -287,7 +290,7 @@ double template_evaluateUsingExcludedAtoms(EnergyNonbond_O *mthis, ScoringFuncti { energyElectrostatic += (e); } #undef NONBOND_EVDW_ENERGY_ACCUMULATE #define NONBOND_EVDW_ENERGY_ACCUMULATE(e) \ - { energyVdw += (e); } + { energyVdw.add(e); } #undef NONBOND_ENERGY_ACCUMULATE #define NONBOND_ENERGY_ACCUMULATE(e) {}; #undef NONBOND_FORCE_ACCUMULATE @@ -306,9 +309,9 @@ double template_evaluateUsingExcludedAtoms(EnergyNonbond_O *mthis, ScoringFuncti #pragma clang diagnostic pop // printf("%s:%d:%s Entering\n", __FILE__, __LINE__, __FUNCTION__ ); num_real x1, y1, z1, x2, y2, z2, dA, dC, dQ1Q2; - double vdwScale = energyFunction->energyScale()->getVdwScale(); - double eelScale = energyFunction->energyScale()->getElectrostaticScale(); - double DIELECTRIC = energyFunction->energyScale()->getDielectricConstant(); + double vdwScale = energyScaleVdwScale(energyScale); + double eelScale = energyScaleElectrostaticScale(energyScale); + double DIELECTRIC = energyScaleDielectricConstant(energyScale); int I1, I2; int i = 0; int endIndex = pos->length() / 3; @@ -387,7 +390,7 @@ double template_evaluateUsingExcludedAtoms(EnergyNonbond_O *mthis, ScoringFuncti #undef CUTOFF_SQUARED #undef DIELECTRIC - MaybeFiniteDiff::maybeTestFiniteDifference(score, I1, I2, activeAtomMask, x1, y1, z1, x2, y2, z2, dA, dC, dQ1Q2, fx1, fy1, + MaybeFiniteDiff::maybeTestFiniteDifference(score, energyScale, I1, I2, activeAtomMask, x1, y1, z1, x2, y2, z2, dA, dC, dQ1Q2, fx1, fy1, fz1, fx2, fy2, fz2, index, fails, debugForce); index++; @@ -449,16 +452,16 @@ double template_evaluateUsingExcludedAtoms(EnergyNonbond_O *mthis, ScoringFuncti LOG("Nonbond energy vdw({}) electrostatic({})\n", (double)mthis->_EnergyVdw, mthis->_EnergyElectrostatic); LOG("Nonbond energy }\n"); maybeSetEnergy(componentEnergy, _sym_energyElectrostaticExcludedAtoms, energyElectrostatic); - maybeSetEnergy(componentEnergy, _sym_energyVdwExcludedAtoms, energyVdw); - return energyElectrostatic + energyVdw; + maybeSetEnergy(componentEnergy, _sym_energyVdwExcludedAtoms, energyVdw.getSum()); + return energyElectrostatic + energyVdw.getSum(); } /*! The core nonbond code using pairwise terms. It as a template function so that template arguments can inline or elide testing code. */ template -//__attribute__((optnone)) -double template_evaluateUsingTerms(EnergyNonbond_O *mthis, ScoringFunction_sp score, NVector_sp pos, core::T_sp componentEnergy, +double template_evaluateUsingTerms(EnergyNonbond_O *mthis, ScoringFunction_sp score, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, bool calcDiagonalHessian, bool calcOffDiagonalHessian, gc::Nilable hessian, gc::Nilable hdvec, gc::Nilable dvec, core::T_sp activeAtomMask, @@ -466,7 +469,7 @@ double template_evaluateUsingTerms(EnergyNonbond_O *mthis, ScoringFunction_sp sc double dielectricConstant; double dQ1Q2Scale; double cutoff; - EnergyFunction_sp energyFunction = energyFunctionNonbondParameters(score, dielectricConstant, dQ1Q2Scale, cutoff); + EnergyFunction_sp energyFunction = energyFunctionNonbondParameters(score, energyScale, dielectricConstant, dQ1Q2Scale, cutoff); double nonbondCutoffSquared = cutoff * cutoff; #define CUTOFF_SQUARED nonbondCutoffSquared MAYBE_SETUP_ACTIVE_ATOM_MASK(); @@ -476,12 +479,12 @@ double template_evaluateUsingTerms(EnergyNonbond_O *mthis, ScoringFunction_sp sc ANN(hdvec); ANN(dvec); double energyElectrostatic = 0.0; - double energyVdw = 0.0; + KahanSummation energyVdw; double energyVdw14 = 0.0; double energyElectrostatic14 = 0.0; - double vdwScale = energyFunction->energyScale()->getVdwScale(); - double eelScale = energyFunction->energyScale()->getElectrostaticScale(); - double DIELECTRIC = energyFunction->energyScale()->getDielectricConstant(); + double vdwScale = energyScaleVdwScale(energyScale); + double eelScale = energyScaleElectrostaticScale(energyScale); + double DIELECTRIC = energyScaleDielectricConstant(energyScale); bool hasForce = force.notnilp(); bool hasHessian = hessian.notnilp(); bool hasHdAndD = (hdvec.notnilp()) && (dvec.notnilp()); @@ -509,7 +512,7 @@ double template_evaluateUsingTerms(EnergyNonbond_O *mthis, ScoringFunction_sp sc if (InteractionIs14) { \ energyVdw14 += (e); \ } else { \ - energyVdw += (e); \ + energyVdw.add(e); \ } \ } @@ -551,7 +554,7 @@ double template_evaluateUsingTerms(EnergyNonbond_O *mthis, ScoringFunction_sp sc #include #undef CUTOFF_SQUARED #undef DIELECTRIC - MaybeFiniteDiff::maybeTestFiniteDifference(score, I1, I2, activeAtomMask, x1, y1, z1, x2, y2, z2, dA, dC, dQ1Q2, fx1, fy1, + MaybeFiniteDiff::maybeTestFiniteDifference(score, energyScale, I1, I2, activeAtomMask, x1, y1, z1, x2, y2, z2, dA, dC, dQ1Q2, fx1, fy1, fz1, fx2, fy2, fz2, index, fails, debugForce); index++; @@ -612,16 +615,19 @@ double template_evaluateUsingTerms(EnergyNonbond_O *mthis, ScoringFunction_sp sc LOG("Nonbond energy }\n"); maybeSetEnergy(componentEnergy, _sym_energyElectrostatic, energyElectrostatic); - maybeSetEnergy(componentEnergy, _sym_energyVdw, energyVdw); + maybeSetEnergy(componentEnergy, _sym_energyVdw, energyVdw.getSum()); maybeSetEnergy(componentEnergy, _sym_energyElectrostatic14, energyElectrostatic14); maybeSetEnergy(componentEnergy, _sym_energyVdw14, energyVdw14); - return energyElectrostatic + energyVdw + energyElectrostatic14 + energyVdw14; + return energyElectrostatic + energyVdw.getSum() + energyElectrostatic14 + energyVdw14; } SYMBOL_EXPORT_SC_(ChemPkg, find_type); -//__attribute__((optnone)) -double EnergyNonbond_O::evaluateAllComponent(ScoringFunction_sp score, NVector_sp pos, core::T_sp componentEnergy, bool calcForce, +double EnergyNonbond_O::evaluateAllComponent(ScoringFunction_sp score, + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, gc::Nilable force, bool calcDiagonalHessian, bool calcOffDiagonalHessian, gc::Nilable hessian, gc::Nilable hdvec, gc::Nilable dvec, core::T_sp activeAtomMask, @@ -633,16 +639,16 @@ double EnergyNonbond_O::evaluateAllComponent(ScoringFunction_sp score, NVector_s size_t index = 0; if (this->_UsesExcludedAtoms) { // Evaluate the nonbonds using the excluded atom list - energy = template_evaluateUsingExcludedAtoms(this, score, pos, componentEnergy, calcForce, force, + energy = template_evaluateUsingExcludedAtoms(this, score, pos, energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, calcOffDiagonalHessian, hessian, hdvec, dvec, activeAtomMask, debugInteractions, fails, index); // Evaluate the 1-4 terms - energy += template_evaluateUsingTerms(this, score, pos, componentEnergy, calcForce, force, + energy += template_evaluateUsingTerms(this, score, pos, energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, calcOffDiagonalHessian, hessian, hdvec, dvec, activeAtomMask, debugInteractions, fails, index); } else { // Evaluate everything using terms - energy = template_evaluateUsingTerms(this, score, pos, componentEnergy, calcForce, force, + energy = template_evaluateUsingTerms(this, score, pos, energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, calcOffDiagonalHessian, hessian, hdvec, dvec, activeAtomMask, debugInteractions, fails, index); } @@ -650,7 +656,11 @@ double EnergyNonbond_O::evaluateAllComponent(ScoringFunction_sp score, NVector_s } CL_DEFMETHOD -double EnergyNonbond_O::debugAllComponent(ScoringFunction_sp score, NVector_sp pos, core::T_sp componentEnergy, bool calcForce, +double EnergyNonbond_O::debugAllComponent(ScoringFunction_sp score, + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, gc::Nilable force, bool calcDiagonalHessian, bool calcOffDiagonalHessian, gc::Nilable hessian, gc::Nilable hdvec, gc::Nilable dvec, core::T_sp activeAtomMask, core::T_sp debugInteractions) { @@ -661,30 +671,30 @@ double EnergyNonbond_O::debugAllComponent(ScoringFunction_sp score, NVector_sp p if (this->_UsesExcludedAtoms) { // Evaluate the nonbonds using the excluded atom list energy = template_evaluateUsingExcludedAtoms( - this, score, pos, componentEnergy, calcForce, force, calcDiagonalHessian, calcOffDiagonalHessian, hessian, hdvec, dvec, + this, score, pos, energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, calcOffDiagonalHessian, hessian, hdvec, dvec, activeAtomMask, debugInteractions, fails, index, true); // Evaluate the 1-4 terms - energy += template_evaluateUsingTerms(this, score, pos, componentEnergy, calcForce, force, + energy += template_evaluateUsingTerms(this, score, pos, energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, calcOffDiagonalHessian, hessian, hdvec, dvec, activeAtomMask, debugInteractions, fails, index, true); } else { // Evaluate everything using terms - energy = template_evaluateUsingTerms(this, score, pos, componentEnergy, calcForce, force, + energy = template_evaluateUsingTerms(this, score, pos, energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, calcOffDiagonalHessian, hessian, hdvec, dvec, activeAtomMask, debugInteractions, fails, index, true); } return energy; } -double _evaluateEnergyOnly_Nonbond(ScoringFunction_sp score, int I1, int I2, core::T_sp activeAtomMask, num_real x1, num_real y1, +double _evaluateEnergyOnly_Nonbond(ScoringFunction_sp score, core::T_sp energyScale, int I1, int I2, core::T_sp activeAtomMask, num_real x1, num_real y1, num_real z1, num_real x2, num_real y2, num_real z2, num_real dA, num_real dC, num_real dQ1Q2) { double dielectricConstant; double dQ1Q2Scale; double cutoff; - auto energyFunction = energyFunctionNonbondParameters(score, dielectricConstant, dQ1Q2Scale, cutoff); - double vdwScale = energyFunction->energyScale()->getVdwScale(); - double eelScale = energyFunction->energyScale()->getElectrostaticScale(); - double DIELECTRIC = energyFunction->energyScale()->getDielectricConstant(); + auto energyFunction = energyFunctionNonbondParameters(score, energyScale, dielectricConstant, dQ1Q2Scale, cutoff); + double vdwScale = energyScaleVdwScale(energyScale); + double eelScale = energyScaleElectrostaticScale(energyScale); + double DIELECTRIC = energyScaleDielectricConstant(energyScale); #define CUTOFF_SQUARED (cutoff * cutoff) MAYBE_SETUP_ACTIVE_ATOM_MASK(); #undef NONBOND_SET_PARAMETER @@ -727,7 +737,9 @@ double _evaluateEnergyOnly_Nonbond(ScoringFunction_sp score, int I1, int I2, cor } CL_LAMBDA((energy-nonbond chem:energy-nonbond) score pos &optional active-atom-mask); -CL_DEFMETHOD void EnergyNonbond_O::compareAnalyticalAndNumericalForceAndHessianTermByTerm(ScoringFunction_sp score, NVector_sp pos, +CL_DEFMETHOD void EnergyNonbond_O::compareAnalyticalAndNumericalForceAndHessianTermByTerm(ScoringFunction_sp score, + NVector_sp pos, + core::T_sp energyScale, core::T_sp activeAtomMask) { double energy = 0.0; size_t fails = 0; @@ -735,16 +747,16 @@ CL_DEFMETHOD void EnergyNonbond_O::compareAnalyticalAndNumericalForceAndHessianT auto force = NVector_O::make(pos->size(), 0.0, true); if (this->_UsesExcludedAtoms) { // Evaluate the nonbonds using the excluded atom list - energy = template_evaluateUsingExcludedAtoms(this, score, pos, nil(), true, force, false, + energy = template_evaluateUsingExcludedAtoms(this, score, pos, energyScale, nil(), true, force, false, false, nil(), nil(), nil(), nil(), nil(), fails, index); // Evaluate the 1-4 terms - energy += template_evaluateUsingTerms(this, score, pos, nil(), true, force, false, false, + energy += template_evaluateUsingTerms(this, score, pos, energyScale, nil(), true, force, false, false, nil(), nil(), nil(), nil(), nil(), fails, index); } else { // Evaluate everything using terms - energy = template_evaluateUsingTerms(this, score, pos, nil(), true, force, false, false, + energy = template_evaluateUsingTerms(this, score, pos, energyScale, nil(), true, force, false, false, nil(), nil(), nil(), nil(), nil(), fails, index); } @@ -863,11 +875,11 @@ void EnergyNonbond_O::setupHessianPreconditioner(NVector_sp nvPosition, Abstract "components of energy are used for the precondition to keep it sparse"); } -CL_DEFMETHOD void EnergyNonbond_O::expandExcludedAtomsToTerms(ScoringFunction_sp score) { +CL_DEFMETHOD void EnergyNonbond_O::expandExcludedAtomsToTerms(ScoringFunction_sp score, core::T_sp energyScale ) { double dielectricConstant; double dQ1Q2Scale; double cutoff; - EnergyFunction_sp energyFunction = energyFunctionNonbondParameters(score, dielectricConstant, dQ1Q2Scale, cutoff); + EnergyFunction_sp energyFunction = energyFunctionNonbondParameters(score, energyScale, dielectricConstant, dQ1Q2Scale, cutoff); // printf("%s:%d In evaluateUsingExcludedAtoms starting this->_DebugEnergy -> %d\n", __FILE__, __LINE__, this->_DebugEnergy ); if (!this->_iac_vec) { SIMPLE_ERROR("The nonbonded excluded atoms parameters have not been set up"); diff --git a/src/chem/energyOutOfZPlane.cc b/src/chem/energyOutOfZPlane.cc index a8041d91..e93100ba 100644 --- a/src/chem/energyOutOfZPlane.cc +++ b/src/chem/energyOutOfZPlane.cc @@ -206,17 +206,18 @@ bool calcOffDiagonalHessian = true; double EnergyOutOfZPlane_O::evaluateAllComponent( ScoringFunction_sp score, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { MAYBE_SETUP_ACTIVE_ATOM_MASK(); MAYBE_SETUP_DEBUG_INTERACTIONS(debugInteractions.notnilp()); diff --git a/src/chem/energyPeriodicBoundaryConditionsNonbond.cc b/src/chem/energyPeriodicBoundaryConditionsNonbond.cc index b5f2b44a..71dd99e4 100644 --- a/src/chem/energyPeriodicBoundaryConditionsNonbond.cc +++ b/src/chem/energyPeriodicBoundaryConditionsNonbond.cc @@ -165,6 +165,7 @@ void EnergyPeriodicBoundaryConditionsNonbond_O::setupHessianPreconditioner(NVect void EnergyPeriodicBoundaryConditionsNonbond_O::evaluateTerms(ScoringFunction_sp score, NVector_sp pos, + core::T_sp energyScale, bool calcForce, gc::Nilable force, bool calcDiagonalHessian, @@ -179,7 +180,7 @@ void EnergyPeriodicBoundaryConditionsNonbond_O::evaluateTerms(ScoringFunction_sp double Evdw = 0.0; double Eeel = 0.0; double cutoff; - energyFunctionNonbondParameters(score,dielectricConstant,dQ1Q2Scale,cutoff); + energyFunctionNonbondParameters(score,energyScale,dielectricConstant,dQ1Q2Scale,cutoff); #define CUTOFF_SQUARED (cutoff*cutoff) double nonbondCutoffSquared = cutoff*cutoff; double nonbondCutoffReciprocal6 = 1.0/(cutoff*cutoff*cutoff*cutoff*cutoff*cutoff); @@ -330,6 +331,7 @@ void EnergyPeriodicBoundaryConditionsNonbond_O::evaluateTerms(ScoringFunction_sp void EnergyPeriodicBoundaryConditionsNonbond_O::evaluateUsingExcludedAtoms(ScoringFunction_sp score, NVector_sp pos, + core::T_sp energyScale, bool calcForce, gc::Nilable force, bool calcDiagonalHessian, @@ -342,7 +344,7 @@ void EnergyPeriodicBoundaryConditionsNonbond_O::evaluateUsingExcludedAtoms(Scori double dielectricConstant; double dQ1Q2Scale; double cutoff; - EnergyFunction_sp energyFunction = energyFunctionNonbondParameters(score,dielectricConstant,dQ1Q2Scale,cutoff); + EnergyFunction_sp energyFunction = energyFunctionNonbondParameters(score,energyScale,dielectricConstant,dQ1Q2Scale,cutoff); double nonbondCutoffSquared = cutoff*cutoff; #define CUTOFF_SQUARED nonbondCutoffSquared double nonbondCutoffReciprocal6 = 1.0/(cutoff*cutoff*cutoff*cutoff*cutoff*cutoff); diff --git a/src/chem/energyPointToLineRestraint.cc b/src/chem/energyPointToLineRestraint.cc index 1ddb837c..868438ab 100644 --- a/src/chem/energyPointToLineRestraint.cc +++ b/src/chem/energyPointToLineRestraint.cc @@ -52,17 +52,18 @@ EnergyPointToLineRestraint_sp EnergyPointToLineRestraint_O::create(EnergySketchS double EnergyPointToLineRestraint_O::evaluateAllComponent( ScoringFunction_sp score, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { MAYBE_SETUP_ACTIVE_ATOM_MASK(); MAYBE_SETUP_DEBUG_INTERACTIONS(debugInteractions.notnilp()); diff --git a/src/chem/energyRigidBodyNonbond.cc b/src/chem/energyRigidBodyNonbond.cc index 4d8c6c49..589b24f3 100644 --- a/src/chem/energyRigidBodyNonbond.cc +++ b/src/chem/energyRigidBodyNonbond.cc @@ -432,22 +432,23 @@ inline num_real periodic_boundary_adjust(const num_real& delta, const num_real& } double EnergyRigidBodyNonbond_O::evaluateAllComponent( ScoringFunction_sp score, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { double dielectricConstant; double dQ1Q2Scale; double cutoff; - energyFunctionNonbondParameters(score,dielectricConstant,dQ1Q2Scale,cutoff); + energyFunctionNonbondParameters(score,energyScale,dielectricConstant,dQ1Q2Scale,cutoff); MAYBE_SETUP_ACTIVE_ATOM_MASK(); MAYBE_SETUP_DEBUG_INTERACTIONS(debugInteractions.notnilp()); // SIMPLE_WARN("FIXactiveAtomMask How do I deal with activeAtomMask"); diff --git a/src/chem/energyRigidBodyStaple.cc b/src/chem/energyRigidBodyStaple.cc index d02e65f8..e0f8a86a 100644 --- a/src/chem/energyRigidBodyStaple.cc +++ b/src/chem/energyRigidBodyStaple.cc @@ -182,17 +182,18 @@ void EnergyRigidBodyStaple_O::setupHessianPreconditioner( double EnergyRigidBodyStaple_O::evaluateAllComponent( ScoringFunction_sp score, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { MAYBE_SETUP_ACTIVE_ATOM_MASK(); MAYBE_SETUP_DEBUG_INTERACTIONS(debugInteractions.notnilp()); diff --git a/src/chem/energySketchNonbond.cc b/src/chem/energySketchNonbond.cc index 3382b2bd..4b07f091 100644 --- a/src/chem/energySketchNonbond.cc +++ b/src/chem/energySketchNonbond.cc @@ -99,20 +99,23 @@ void EnergySketchNonbond_O::setupHessianPreconditioner( } double EnergySketchNonbond_O::evaluateAllComponent( ScoringFunction_sp score, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask, - core::T_sp debugInteractions ) + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { // Evaluate everything using terms - double totalEnergy = this->evaluateTerms(pos,componentEnergy, + double totalEnergy = this->evaluateTerms(pos, + energyScale, + componentEnergy, calcForce,force, calcDiagonalHessian, calcOffDiagonalHessian, hessian,hdvec,dvec,activeAtomMask); return totalEnergy; @@ -121,15 +124,16 @@ double EnergySketchNonbond_O::evaluateAllComponent( ScoringFunction_sp score, double EnergySketchNonbond_O::evaluateTerms(NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask) + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask) { MAYBE_SETUP_ACTIVE_ATOM_MASK(); core::T_sp debugInteractions = nil(); diff --git a/src/chem/energySketchStretch.cc b/src/chem/energySketchStretch.cc index 910b42f1..6798544c 100644 --- a/src/chem/energySketchStretch.cc +++ b/src/chem/energySketchStretch.cc @@ -179,6 +179,7 @@ void EnergySketchStretch_O::setupHessianPreconditioner( double EnergySketchStretch_O::evaluateAllComponent( ScoringFunction_sp score, NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, diff --git a/src/chem/energyStretch.cc b/src/chem/energyStretch.cc index 1af4b66c..c269706d 100644 --- a/src/chem/energyStretch.cc +++ b/src/chem/energyStretch.cc @@ -358,17 +358,18 @@ CL_DEFMETHOD core::T_sp EnergyStretch_O::stretchTermBetweenAtoms(Atom_sp x, Atom } double EnergyStretch_O::evaluateAllComponent( ScoringFunction_sp score, - NVector_sp pos, - core::T_sp componentEnergy, - bool calcForce, - gc::Nilable force, - bool calcDiagonalHessian, - bool calcOffDiagonalHessian, - gc::Nilable hessian, - gc::Nilable hdvec, - gc::Nilable dvec, - core::T_sp activeAtomMask , - core::T_sp debugInteractions ) + NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + bool calcForce, + gc::Nilable force, + bool calcDiagonalHessian, + bool calcOffDiagonalHessian, + gc::Nilable hessian, + gc::Nilable hdvec, + gc::Nilable dvec, + core::T_sp activeAtomMask , + core::T_sp debugInteractions ) { MAYBE_SETUP_ACTIVE_ATOM_MASK(); MAYBE_SETUP_DEBUG_INTERACTIONS(debugInteractions.notnilp()); diff --git a/src/chem/minimizer.cc b/src/chem/minimizer.cc index 77a39237..6bca9045 100644 --- a/src/chem/minimizer.cc +++ b/src/chem/minimizer.cc @@ -308,11 +308,16 @@ void Minimizer_O::getPosition( NVector_sp nvResult, * in nvForce (which must be an initialized vector, * but does not need to be filled with zeros). */ -double Minimizer_O::dTotalEnergy( NVector_sp nvPos, core::T_sp activeAtomMask ) +double Minimizer_O::dTotalEnergy( NVector_sp pos, + core::T_sp energyScale, + core::T_sp activeAtomMask ) { double dEnergy; - dEnergy = this->_ScoringFunction->evaluateEnergy( nvPos, activeAtomMask ); + dEnergy = this->_ScoringFunction->evaluateEnergy( pos, + energyScale, + nil(), + activeAtomMask ); return(dEnergy); } @@ -327,9 +332,13 @@ double Minimizer_O::dTotalEnergy( NVector_sp nvPos, core::T_sp activeAtomMask ) * in nvForce (which must be an initialized vector, * but does not need to be filled with zeros). */ -double Minimizer_O::dTotalEnergyForce( NVector_sp nvPos, NVector_sp nvForce, core::T_sp activeAtomMask ) +double Minimizer_O::dTotalEnergyForce( NVector_sp pos, + core::T_sp energyScale, + NVector_sp nvForce, core::T_sp activeAtomMask ) { - return this->_ScoringFunction->evaluateEnergyForce(nvPos,true,nvForce,activeAtomMask); + return this->_ScoringFunction->evaluateEnergyForce(pos, + energyScale, + true,nvForce,activeAtomMask); } @@ -338,7 +347,7 @@ double Minimizer_O::dTotalEnergyForce( NVector_sp nvPos, NVector_sp nvForce, cor * * Calculate the energy along a 1D coordinate */ -double Minimizer_O::d1DTotalEnergy( double x, core::T_sp activeAtomMask ) +double Minimizer_O::d1DTotalEnergy( double x, core::T_sp energyScale, core::T_sp activeAtomMask ) { #ifdef DEBUG_ON // this->nvP1DSearchOrigin->debugDump("origin"); @@ -346,7 +355,7 @@ double Minimizer_O::d1DTotalEnergy( double x, core::T_sp activeAtomMask ) // this->nvP1DSearchDirection->debugDump("direction"); #endif this->getPosition(this->nvP1DSearchTemp1, this->nvP1DSearchOrigin, this->nvP1DSearchDirection,x, activeAtomMask ); - return this->_ScoringFunction->evaluateEnergy(this->nvP1DSearchTemp1,activeAtomMask); + return this->_ScoringFunction->evaluateEnergy(this->nvP1DSearchTemp1,energyScale,nil(),activeAtomMask); } @@ -355,7 +364,7 @@ double Minimizer_O::d1DTotalEnergy( double x, core::T_sp activeAtomMask ) * * Calculate the energy/derivative along a 1D coordinate */ -double Minimizer_O::d1DTotalEnergyForce( double x, double* fx, double * dfx, core::T_sp activeAtomMask ) +double Minimizer_O::d1DTotalEnergyForce( double x, core::T_sp energyScale, double* fx, double * dfx, core::T_sp activeAtomMask ) { #ifdef DEBUG_ON // this->nvP1DSearchOrigin->debugDump("origin"); @@ -366,6 +375,7 @@ double Minimizer_O::d1DTotalEnergyForce( double x, double* fx, double * dfx, cor this->nvP1DSearchOrigin, this->nvP1DSearchDirection, x, activeAtomMask ); *fx = this->_ScoringFunction->evaluateEnergyForce( this->nvP1DSearchTemp1, + energyScale, true, this->nvP1DSearchTemp2, activeAtomMask ); *dfx = -dotProductWithActiveAtomMask(this->nvP1DSearchTemp2,this->nvP1DSearchDirection,activeAtomMask); @@ -387,6 +397,7 @@ double Minimizer_O::d1DTotalEnergyForce( double x, double* fx, double * dfx, cor */ void Minimizer_O::minBracket( NVector_sp nvOrigin, + core::T_sp energyScale, NVector_sp nvDir, double *dPxa, double *dPxb, @@ -402,15 +413,15 @@ void Minimizer_O::minBracket( this->_MinBracketSteps = 0; xa = *dPxa; xb = *dPxb; - fa = this->d1DTotalEnergy(xa,activeAtomMask); - fb = this->d1DTotalEnergy(xb,activeAtomMask); + fa = this->d1DTotalEnergy(xa,energyScale,activeAtomMask); + fb = this->d1DTotalEnergy(xb,energyScale,activeAtomMask); // Make sure that we are going downhill a->b if ( fb > fa ) { SWAP_VALUES( xa, xb, temp ); SWAP_VALUES( fa, fb, temp ); } xc = xb+GOLD*(xb-xa); - fc = this->d1DTotalEnergy(xc,activeAtomMask); + fc = this->d1DTotalEnergy(xc,energyScale,activeAtomMask); LOG("Start: xa({}) xb({}) xc({}) | fa({}) fb({}) fc({})" , xa , xb , xc , fa , fb , fc ); while ( fb > fc ) { numSteps++; @@ -426,7 +437,7 @@ void Minimizer_O::minBracket( (2.0*SIGN(MAX(fabs(q-r),TINY),q-r)); ulim = (xb)+GLIMIT*(xc-xb); if (( xb-u)*(u-xc)>0.0) { - fu = this->d1DTotalEnergy(u,activeAtomMask); + fu = this->d1DTotalEnergy(u,energyScale,activeAtomMask); if ( fud1DTotalEnergy(u,activeAtomMask); + fu = this->d1DTotalEnergy(u,energyScale,activeAtomMask); } else if ((xc-u)*(u-ulim) > 0.0) { - fu = this->d1DTotalEnergy(u,activeAtomMask); + fu = this->d1DTotalEnergy(u,energyScale,activeAtomMask); if ( fu < fc ) { LEFT_SHIFT_VALUES( xb, xc, u, xc+GOLD*(xc-xb) ); - LEFT_SHIFT_VALUES( fb, fc, fu, this->d1DTotalEnergy(u,activeAtomMask)); + LEFT_SHIFT_VALUES( fb, fc, fu, this->d1DTotalEnergy(u,energyScale,activeAtomMask)); } // Limit parabolic u to maximum allowed value } else if ((u-ulim)*(ulim-xc)>=0.0 ){ u = ulim; - fu = this->d1DTotalEnergy(u,activeAtomMask); + fu = this->d1DTotalEnergy(u,energyScale,activeAtomMask); } else { u = xc+GOLD*(xc-xb); - fu = this->d1DTotalEnergy(u,activeAtomMask); + fu = this->d1DTotalEnergy(u,energyScale,activeAtomMask); } LEFT_SHIFT_VALUES( xa, xb, xc, u ); LEFT_SHIFT_VALUES( fa, fb, fc, fu ); @@ -488,6 +499,7 @@ void Minimizer_O::minBracket( * new minimum in (xmin) */ double Minimizer_O::dbrent( double ax, double bx, double cx, + core::T_sp energyScale, double tol, double& lineStep, int& energyEvals, @@ -516,7 +528,7 @@ double Minimizer_O::dbrent( double ax, double bx, double cx, // // Calculate the derivative along the search direction // - ft = d1DTotalEnergyForce( x, &fx, &dx, activeAtomMask ); + ft = d1DTotalEnergyForce( x, energyScale, &fx, &dx, activeAtomMask ); forceEvals++; LOG("dbrent: derivative x,fx,dx = {} {} {}" , x , fx , dx ); fw=fv=fx; @@ -597,11 +609,11 @@ double Minimizer_O::dbrent( double ax, double bx, double cx, } if (fabs(_d) >= tol1) { u=x+_d; - fu=d1DTotalEnergy(u,activeAtomMask); + fu=d1DTotalEnergy(u,energyScale,activeAtomMask); energyEvals++; } else { u = x+SIGN(tol1,_d); - fu=d1DTotalEnergy(u,activeAtomMask); + fu=d1DTotalEnergy(u,energyScale,activeAtomMask); energyEvals++; if ( fu>fx) { // If the minimum lineStep in the downhill direction takes us // uphill, then we are done @@ -617,7 +629,7 @@ double Minimizer_O::dbrent( double ax, double bx, double cx, // // Calculate the force along the search direction // - ft = d1DTotalEnergyForce( u, &ft, &du, activeAtomMask ); // Now housekeeping + ft = d1DTotalEnergyForce( u, energyScale, &ft, &du, activeAtomMask ); // Now housekeeping forceEvals++; // if (fu<=fx) { @@ -649,7 +661,9 @@ double Minimizer_O::dbrent( double ax, double bx, double cx, void Minimizer_O::lineSearchInitialReport( StepReport_sp report, - NVector_sp nvPos, NVector_sp nvDir, NVector_sp nvForce, + NVector_sp nvPos, + core::T_sp energyScale, + NVector_sp nvDir, NVector_sp nvForce, double xa, double xb, double xc, double fa, double fb, double fc, core::T_sp activeAtomMask ) @@ -674,7 +688,7 @@ void Minimizer_O::lineSearchInitialReport( StepReport_sp report, report->_Fc = fc; report->_MinBracketSteps = this->_MinBracketSteps; report->_EnergyTermsEnabled = this->_ScoringFunction->energyTermsEnabled(); - report->_TotalEnergy = this->d1DTotalEnergy(0.0,activeAtomMask); + report->_TotalEnergy = this->d1DTotalEnergy(0.0,energyScale,activeAtomMask); report->_DirectionMagnitude = magnitudeWithActiveAtomMask(nvDir,activeAtomMask); report->_ForceMagnitude = magnitudeWithActiveAtomMask(nvForce,activeAtomMask); report->_MinimizerStatus = this->statusAsString(); @@ -706,7 +720,7 @@ void Minimizer_O::lineSearchInitialReport( StepReport_sp report, report->_FixedNonbondRestraintEnergyFn = NumericalFunction_O::create(core::SimpleBaseString_O::make("Alpha"),core::SimpleBaseString_O::make("FixedNonbondRestraint"),xmin,xinc); for ( zx=dxa;zx<=dxc;zx+=(dxc-dxa)/100.0 ) { - zy = this->d1DTotalEnergy(zx,activeAtomMask); + zy = this->d1DTotalEnergy(zx,energyScale,activeAtomMask); report->_TotalEnergyFn->appendValue(zy); #if 0 // skipping components - it's not general @@ -768,6 +782,7 @@ void Minimizer_O::stepReport( StepReport_sp report, double energy, NVector_sp fo void Minimizer_O::lineSearch( double *dPstep, double *dPfnew, NVector_sp nvOrigin, + core::T_sp energyScale, NVector_sp nvDirection, NVector_sp nvForce, NVector_sp nvTemp1, @@ -804,7 +819,7 @@ void Minimizer_O::lineSearch( double *dPstep, core::clasp_write_string(fmt::format("Trying to bracket min this->_InitialLineSearchStep,directionMag,xb -> {},{},{}\n" , this->_InitialLineSearchStep , directionMag , xb)); } - this->minBracket( nvOrigin, nvDirection, + this->minBracket( nvOrigin, energyScale, nvDirection, &xa, &xb, &xc, &fa, &fb, &fc, activeAtomMask ); LOG("Bracketed minimum" ); LOG("xa,xb,xc = {} {} {} " , xa , xb , xc ); @@ -816,7 +831,11 @@ void Minimizer_O::lineSearch( double *dPstep, if ( this->_DebugOn ) { - this->lineSearchInitialReport(report,nvOrigin,nvDirection,nvForce, + this->lineSearchInitialReport(report, + nvOrigin, + energyScale, + nvDirection, + nvForce, xa,xb,xc,fa,fb,fc,activeAtomMask); } @@ -828,7 +847,7 @@ void Minimizer_O::lineSearch( double *dPstep, int energyEvals = 0; int forceEvals = 0; int dbrentSteps = 0; - fb = dbrent( xa, xb, xc, TOL, step, energyEvals, forceEvals, dbrentSteps, activeAtomMask ); + fb = dbrent( xa, xb, xc, energyScale, TOL, step, energyEvals, forceEvals, dbrentSteps, activeAtomMask ); LOG("brent bracketed step = {}" , step ); if ( this->_DebugOn ) { @@ -937,6 +956,7 @@ bool Minimizer_O::_displayIntermediateMessage(NVector_sp pos, */ void Minimizer_O::_steepestDescent( int numSteps, NVector_sp pos, + core::T_sp energyScale, double forceTolerance, core::T_sp activeAtomMask, core::T_sp callback ) @@ -988,7 +1008,7 @@ void Minimizer_O::_steepestDescent( int numSteps, // Done innerSteps = MIN(iRestartSteps,ITMAX); LOG("step" ); - double fp = this->dTotalEnergyForce( pos, forceVec, activeAtomMask ); + double fp = this->dTotalEnergyForce( pos, energyScale, forceVec, activeAtomMask ); LOG("step" ); // r->inPlaceTimesScalar(-1.0); // no preconditioning @@ -1032,7 +1052,7 @@ void Minimizer_O::_steepestDescent( int numSteps, // Lets save the current conformation // before throwing this higher // - fp = dTotalEnergyForce( pos, forceVec, activeAtomMask ); + fp = dTotalEnergyForce( pos, energyScale, forceVec, activeAtomMask ); this->_ScoringFunction->saveCoordinatesAndForcesFromVectors(pos,forceVec); ERROR(_sym_MinimizerExceededSD_MaxSteps, (ql::list() << kw::_sym_minimizer << this->asSmartPtr() @@ -1081,7 +1101,7 @@ void Minimizer_O::_steepestDescent( int numSteps, steepestDescent = true; cosAngle = 0.0; } - this->lineSearch( &step, &fnew, pos, dirVec, forceVec, tv1, tv2, localSteps, stepReport, activeAtomMask ); + this->lineSearch( &step, &fnew, pos, energyScale, dirVec, forceVec, tv1, tv2, localSteps, stepReport, activeAtomMask ); prevStep = step; printedLatestMessage = false; if ( this->_PrintIntermediateResults ) { @@ -1104,7 +1124,7 @@ void Minimizer_O::_steepestDescent( int numSteps, inPlaceAddTimesScalarWithActiveAtomMask( pos, dirVec, step, activeAtomMask ); // r = -f'(x) r == force!!!! - fp = dTotalEnergyForce( pos, forceVec, activeAtomMask ); + fp = dTotalEnergyForce( pos, energyScale, forceVec, activeAtomMask ); if (callback.notnilp()) core::eval::funcall( callback, _sym_steepest_descent, this->_ScoringFunction, pos, activeAtomMask ); if ( this->_DebugOn ) { this->stepReport(stepReport,fp,forceVec,activeAtomMask); @@ -1130,7 +1150,7 @@ void Minimizer_O::_steepestDescent( int numSteps, if ( this->_PrintIntermediateResults && !printedLatestMessage ) { this->_displayIntermediateMessage(pos,step,fnew,forceRmsMag,cosAngle,steepestDescent,false,activeAtomMask); } - fp = dTotalEnergyForce( pos, forceVec, activeAtomMask ); + fp = dTotalEnergyForce( pos, energyScale, forceVec, activeAtomMask ); this->_ScoringFunction->saveCoordinatesAndForcesFromVectors(pos,forceVec); LOG("Wrote coordinates and force to atoms" ); if ( this->_DebugOn ) @@ -1154,6 +1174,7 @@ void Minimizer_O::_steepestDescent( int numSteps, void Minimizer_O::_conjugateGradient(int numSteps, NVector_sp x, + core::T_sp energyScale, double forceTolerance, core::T_sp activeAtomMask, core::T_sp callback ) @@ -1199,7 +1220,7 @@ void Minimizer_O::_conjugateGradient(int numSteps, tv2 = NVector_O::create(iRestartSteps); // Done innerSteps = MIN(iRestartSteps,ITMAX); - double fp = dTotalEnergyForce( x, force, activeAtomMask ); + double fp = dTotalEnergyForce( x, energyScale, force, activeAtomMask ); // r->inPlaceTimesScalar(-1.0); // TODO calculate preconditioner here // s = M^(-1)r rather than just copying it from r @@ -1274,7 +1295,7 @@ void Minimizer_O::_conjugateGradient(int numSteps, // Lets save the current conformation // before throwing this higher // - fp = dTotalEnergyForce( x, force, activeAtomMask ); + fp = dTotalEnergyForce( x, energyScale, force, activeAtomMask ); this->_ScoringFunction->saveCoordinatesAndForcesFromVectors(x,force); ERROR(_sym_MinimizerExceededCG_MaxSteps, (ql::list() << kw::_sym_minimizer << this->asSmartPtr() @@ -1294,7 +1315,7 @@ void Minimizer_O::_conjugateGradient(int numSteps, // Lets save the current conformation // before throwing this higher // - fp = dTotalEnergyForce( x, force, activeAtomMask ); + fp = dTotalEnergyForce( x, energyScale, force, activeAtomMask ); this->_ScoringFunction->saveCoordinatesAndForcesFromVectors(x,force); MINIMIZER_STUCK_ERROR("Stuck in conjugate gradients"); } @@ -1355,7 +1376,7 @@ void Minimizer_O::_conjugateGradient(int numSteps, steepestDescent = true; } - this->lineSearch( &step, &fnew, x, d, force, + this->lineSearch( &step, &fnew, x, energyScale, d, force, tv1, tv2, localSteps, stepReport, activeAtomMask ); @@ -1367,7 +1388,7 @@ void Minimizer_O::_conjugateGradient(int numSteps, // x = x + (step)d // r = -f'(x) r == force!!!! inPlaceAddTimesScalarWithActiveAtomMask(x, d, step, activeAtomMask ); - fp = dTotalEnergyForce( x, force, activeAtomMask ); + fp = dTotalEnergyForce( x, energyScale, force, activeAtomMask ); if (callback.notnilp()) core::eval::funcall( callback, _sym_conjugate_gradient, this->_ScoringFunction, x, activeAtomMask ); if ( this->_DebugOn ) @@ -1436,7 +1457,7 @@ void Minimizer_O::_conjugateGradient(int numSteps, if ( this->_PrintIntermediateResults && !printedLatestMessage ) { this->_displayIntermediateMessage(x,step,fnew,forceRmsMag,cosAngle,steepestDescent,false,activeAtomMask); } - fp = dTotalEnergyForce( x, force, activeAtomMask ); + fp = dTotalEnergyForce( x, energyScale, force, activeAtomMask ); this->_ScoringFunction->saveCoordinatesAndForcesFromVectors(x,force); if ( this->_DebugOn ) { @@ -1450,6 +1471,7 @@ void Minimizer_O::_conjugateGradient(int numSteps, void Minimizer_O::_truncatedNewtonInnerLoop(int kk, NVector_sp xk, + core::T_sp energyScale, SparseLargeSquareMatrix_sp mprecon, SparseLargeSquareMatrix_sp ldlt, NVector_sp force, @@ -1525,6 +1547,7 @@ void Minimizer_O::_truncatedNewtonInnerLoop(int kk, // this->_ScoringFunction->evaluateAll( xk, + energyScale, nil(), true, nvDummy, true, true, nmDummy, @@ -1634,9 +1657,9 @@ void Minimizer_O::_truncatedNewtonInnerLoop(int kk, #define SQRT_EPSILONF 1.0e-5 #define CUBERT_EPSILONF 4.6416e-4 -__attribute__((optnone)) void Minimizer_O::_truncatedNewton(int numSteps, NVector_sp xK, + core::T_sp energyScale, double forceTolerance, core::T_sp activeAtomMask, core::T_sp callback ) @@ -1693,7 +1716,7 @@ void Minimizer_O::_truncatedNewton(int numSteps, // Evaluate initial energy and force // LOG("Evaluating initial energy and force" ); - energyXkNext = dTotalEnergyForce( xK, forceK, activeAtomMask ); + energyXkNext = dTotalEnergyForce( xK, energyScale, forceK, activeAtomMask ); rmsForceMag = rmsMagnitudeWithActiveAtomMask(forceK,activeAtomMask); if ( this->_PrintIntermediateResults ) { @@ -1739,7 +1762,7 @@ void Minimizer_O::_truncatedNewton(int numSteps, // // Inner loop // - _truncatedNewtonInnerLoop( kk, xK, opt_mprecon, opt_ldlt, + _truncatedNewtonInnerLoop( kk, xK, energyScale, opt_mprecon, opt_ldlt, forceK, rmsForceMag, dirVec, dirVecNext, rj, dj, zj, qj, innerLoopDeltaJ1, @@ -1775,7 +1798,7 @@ void Minimizer_O::_truncatedNewton(int numSteps, } energyXk = energyXkNext; - this->lineSearch( &alphaK, &energyXkNext, xK, dirVec, forceK, zj, qj, kk, stepReport, activeAtomMask ); + this->lineSearch( &alphaK, &energyXkNext, xK, energyScale, dirVec, forceK, zj, qj, kk, stepReport, activeAtomMask ); if (this->_StepCallback.notnilp()) { core::DoubleFloat_sp dstep = core::DoubleFloat_O::create(alphaK); core::eval::funcall(this->_StepCallback, _sym_truncated_newton, xK, forceK, dstep, dirVec, opt_mprecon, opt_ldlt ); @@ -1785,7 +1808,7 @@ void Minimizer_O::_truncatedNewton(int numSteps, // // Evaluate the force at the new position // - fp = dTotalEnergyForce( xKNext, forceK, activeAtomMask ); + fp = dTotalEnergyForce( xKNext, energyScale, forceK, activeAtomMask ); if ( this->_DebugOn ) { this->stepReport(stepReport,fp,forceK,activeAtomMask); @@ -1869,7 +1892,7 @@ void Minimizer_O::_truncatedNewton(int numSteps, // Lets save the current conformation // before throwing this higher // - dTotalEnergyForce( xK, forceK, activeAtomMask ); + dTotalEnergyForce( xK, energyScale, forceK, activeAtomMask ); this->_ScoringFunction->saveCoordinatesAndForcesFromVectors(xK,forceK); ERROR(_sym_MinimizerExceededTN_MaxSteps, (ql::list() << kw::_sym_minimizer << this->asSmartPtr() @@ -1885,7 +1908,7 @@ void Minimizer_O::_truncatedNewton(int numSteps, } } copyVector(xK,xKNext); - dTotalEnergyForce( xK, forceK, activeAtomMask ); + dTotalEnergyForce( xK, energyScale, forceK, activeAtomMask ); this->_ScoringFunction->saveCoordinatesAndForcesFromVectors(xK,forceK); if ( this->_DebugOn ) { @@ -1905,7 +1928,7 @@ void Minimizer_O::_truncatedNewton(int numSteps, */ -void Minimizer_O::_evaluateEnergyAndForceManyTimes( int numSteps, NVector_sp nvPos, core::T_sp activeAtomMask ) +void Minimizer_O::_evaluateEnergyAndForceManyTimes( int numSteps, NVector_sp nvPos, core::T_sp energyScale, core::T_sp activeAtomMask ) { #define MINSLOPE 0.000001 #define MINCHANGE 0.01 @@ -1928,7 +1951,7 @@ void Minimizer_O::_evaluateEnergyAndForceManyTimes( int numSteps, NVector_sp nv iCount = 0; this->_Iteration = 1; do { - dEnergy = this->dTotalEnergyForce( nvPos, nvNewForce, activeAtomMask ); + dEnergy = this->dTotalEnergyForce( nvPos, energyScale, nvNewForce, activeAtomMask ); if ( iCount % 10000 == 0 ) { core::clasp_writeln_string(fmt::format("Evaluating energy step#{}" , iCount )); } @@ -1940,11 +1963,11 @@ void Minimizer_O::_evaluateEnergyAndForceManyTimes( int numSteps, NVector_sp nv } -void Minimizer_O::validateForce(NVector_sp pos, NVector_sp force, core::T_sp activeAtomMask ) +void Minimizer_O::validateForce(NVector_sp pos, core::T_sp energyScale, NVector_sp force, core::T_sp activeAtomMask ) { ForceMatchReport_sp report; if ( this->_DebugOn ) { - report = this->_ScoringFunction->checkIfAnalyticalForceMatchesNumericalForce(pos,force,activeAtomMask); + report = this->_ScoringFunction->checkIfAnalyticalForceMatchesNumericalForce(pos,energyScale,force,activeAtomMask); this->_Log->addReport(report); } } @@ -2054,30 +2077,30 @@ CL_DEFMETHOD void Minimizer_O::setEnergyFunction(ScoringFunction_sp f) } CL_LISPIFY_NAME("evaluateEnergyAndForceManyTimes"); -CL_DEFMETHOD void Minimizer_O::evaluateEnergyAndForceManyTimes(int numSteps, core::T_sp activeAtomMask) +CL_DEFMETHOD void Minimizer_O::evaluateEnergyAndForceManyTimes(core::T_sp energyScale, int numSteps, core::T_sp activeAtomMask) { NVector_sp pos; ASSERT(this->_ScoringFunction); this->_Iteration = 1; pos = NVector_O::create(this->_ScoringFunction->getNVectorSize()); this->_ScoringFunction->loadCoordinatesIntoVector(pos); - this->_evaluateEnergyAndForceManyTimes(numSteps,pos,activeAtomMask); + this->_evaluateEnergyAndForceManyTimes(numSteps,pos,energyScale,activeAtomMask); } CL_LISPIFY_NAME("resetAndMinimize"); -CL_LAMBDA((minimizer chem:minimizer) &optional active-atom-mask callback ); -CL_DEFMETHOD void Minimizer_O::resetAndMinimize(core::T_sp activeAtomMask, core::T_sp callback ) +CL_LAMBDA((minimizer chem:minimizer) &key energy-scale active-atom-mask callback ); +CL_DEFMETHOD void Minimizer_O::resetAndMinimize(core::T_sp energyScale, core::T_sp activeAtomMask, core::T_sp callback ) { this->_Status = minimizerIdle; - this->minimize(activeAtomMask, callback ); + this->minimize(energyScale, activeAtomMask, callback ); } CL_LISPIFY_NAME("minimize"); -CL_LAMBDA((minimizer chem:minimizer) &optional active-atom-mask callback); -CL_DEFMETHOD core::T_mv Minimizer_O::minimize(core::T_sp activeAtomMask, core::T_sp callback) +CL_LAMBDA((minimizer chem:minimizer) &key energy-scale active-atom-mask callback); +CL_DEFMETHOD core::T_mv Minimizer_O::minimize(core::T_sp energyScale, core::T_sp activeAtomMask, core::T_sp callback) { NVector_sp pos; int retries; @@ -2097,7 +2120,7 @@ CL_DEFMETHOD core::T_mv Minimizer_O::minimize(core::T_sp activeAtomMask, core::T } if ( this->_NumberOfSteepestDescentSteps > 0 ) { this->_steepestDescent( this->_NumberOfSteepestDescentSteps, - pos, this->_SteepestDescentTolerance, activeAtomMask, callback ); + pos, energyScale, this->_SteepestDescentTolerance, activeAtomMask, callback ); } else { if ( this->_PrintIntermediateResults ) { core::clasp_writeln_string("======= Skipping Steepest Descent #steps = 0"); @@ -2105,7 +2128,7 @@ CL_DEFMETHOD core::T_mv Minimizer_O::minimize(core::T_sp activeAtomMask, core::T } if ( this->_NumberOfConjugateGradientSteps > 0 ) { this->_conjugateGradient( this->_NumberOfConjugateGradientSteps, - pos, this->_ConjugateGradientTolerance, activeAtomMask, callback ); + pos, energyScale, this->_ConjugateGradientTolerance, activeAtomMask, callback ); } else { if ( this->_PrintIntermediateResults ) { core::clasp_writeln_string("======= Skipping Conjugate Gradients #steps = 0"); @@ -2113,7 +2136,7 @@ CL_DEFMETHOD core::T_mv Minimizer_O::minimize(core::T_sp activeAtomMask, core::T } if ( this->_NumberOfTruncatedNewtonSteps > 0 ) { this->_truncatedNewton( this->_NumberOfTruncatedNewtonSteps, - pos, this->_TruncatedNewtonTolerance, activeAtomMask, callback ); + pos, energyScale, this->_TruncatedNewtonTolerance, activeAtomMask, callback ); } else { if ( this->_PrintIntermediateResults ) { core::clasp_writeln_string("======= Skipping Truncated Newton #steps = 0"); @@ -2138,7 +2161,7 @@ CL_DEFMETHOD core::T_mv Minimizer_O::minimize(core::T_sp activeAtomMask, core::T << kw::_sym_minimizer << this->asSmartPtr() << kw::_sym_coordinates << pos).result()); DONE: - return Values(pos,mk_double_float(this->dTotalEnergy(pos,activeAtomMask))); + return Values(pos,mk_double_float(this->dTotalEnergy(pos,energyScale,activeAtomMask))); } diff --git a/src/chem/rigidBodyEnergyFunction.cc b/src/chem/rigidBodyEnergyFunction.cc index 3f7daa8a..28020061 100644 --- a/src/chem/rigidBodyEnergyFunction.cc +++ b/src/chem/rigidBodyEnergyFunction.cc @@ -184,7 +184,7 @@ double RigidBodyEnergyFunction_O::evaluateRaw( NVector_sp pos, NVector_sp force adapt::QDomNode_sp identifyTermsBeyondThreshold(); // uint countBadVdwInteractions(double scaleSumOfVdwRadii, geom::DisplayList_sp displayIn); -ForceMatchReport_sp RigidBodyEnergyFunction_O::checkIfAnalyticalForceMatchesNumericalForce( NVector_sp pos, NVector_sp force, core::T_sp activeAtomMask ) { +ForceMatchReport_sp RigidBodyEnergyFunction_O::checkIfAnalyticalForceMatchesNumericalForce( NVector_sp pos, core::T_sp energyScale, NVector_sp force, core::T_sp activeAtomMask ) { IMPLEMENT_ME(); } @@ -213,6 +213,7 @@ void RigidBodyEnergyFunction_O::setOptions( core::List_sp options ) double RigidBodyEnergyFunction_O::evaluateAll( NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, @@ -263,6 +264,7 @@ double RigidBodyEnergyFunction_O::evaluateAll( NVector_sp pos, if (term->isEnabled()) { totalEnergy += term->evaluateAllComponent(this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, @@ -368,6 +370,7 @@ CL_LISPIFY_NAME("rigid-body-velocity-verlet-step-limit-displacement"); DOCGROUP(cando); CL_DEFUN size_t chem__rigid_body_velocity_verlet_step_limit_displacement(ScoringFunction_sp scoringFunc, NVector_sp position, + core::T_sp energyScale, NVector_sp velocity, NVector_sp force, NVector_sp force_dt, @@ -431,7 +434,9 @@ CL_DEFUN size_t chem__rigid_body_velocity_verlet_step_limit_displacement(Scoring } body_idx++; } - scoringFunc->evaluateEnergyForce(position,true,force_dt,tunfrozen); + scoringFunc->evaluateEnergyForce(position, + energyScale, + true,force_dt,tunfrozen); body_idx = 0; for ( size_t idx = 0; idxsize(); idx+=7 ) { if (!unfrozen || unfrozen->testBit(body_idx)==1) { diff --git a/src/chem/scoringFunction.cc b/src/chem/scoringFunction.cc index ec2b5864..91c47579 100644 --- a/src/chem/scoringFunction.cc +++ b/src/chem/scoringFunction.cc @@ -62,6 +62,38 @@ __END_DOC #include +namespace chem { + +double energyScaleDielectricConstant(core::T_sp energyScale) { + if (gc::IsA(energyScale)) { + return gc::As_unsafe(energyScale)->getDielectricConstant(); + } + return 1.0; +} + +double energyScaleVdwScale(core::T_sp energyScale) { + if (gc::IsA(energyScale)) { + return gc::As_unsafe(energyScale)->getVdwScale(); + } + return 1.0; +} + +double energyScaleElectrostaticScale(core::T_sp energyScale) { + if (gc::IsA(energyScale)) { + return gc::As_unsafe(energyScale)->getElectrostaticScale(); + } + return 1.0; +} + +double energyScaleNonbondCutoff(core::T_sp energyScale) { + if (gc::IsA(energyScale)) { + return gc::As_unsafe(energyScale)->getNonbondCutoff(); + } + return 1.0; +} + +}; + namespace chem { @@ -93,12 +125,17 @@ CL_DEFMETHOD void ScoringFunction_O::setAtomTypes(core::HashTable_sp atomTypes) } CL_LISPIFY_NAME("evaluateEnergy"); -CL_LAMBDA((scoring-function chem:scoring-function) positions &optional active-atom-mask debug-interactions); -CL_DEFMETHOD double ScoringFunction_O::evaluateEnergy( NVector_sp pos, core::T_sp activeAtomMask, core::T_sp debugInteractions ) +CL_LAMBDA((scoring-function chem:scoring-function) positions &key energy-scale component-energy active-atom-mask debug-interactions); +CL_DEFMETHOD double ScoringFunction_O::evaluateEnergy( NVector_sp pos, + core::T_sp energyScale, + core::T_sp componentEnergy, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { double energy; energy = this->evaluateAll(pos, - nil(), + energyScale, + componentEnergy, false, nil(), false, false, @@ -113,8 +150,13 @@ CL_DEFMETHOD double ScoringFunction_O::evaluateEnergy( NVector_sp pos, core::T_s CL_LISPIFY_NAME("evaluateEnergyForce"); -CL_LAMBDA((scoring-function chem:scoring-function) positions calc-force force &optional active-atom-mask debug-interactions); -CL_DEFMETHOD double ScoringFunction_O::evaluateEnergyForce( NVector_sp pos, bool calcForce, NVector_sp force, core::T_sp activeAtomMask, core::T_sp debugInteractions ) +CL_LAMBDA((scoring-function chem:scoring-function) positions &key energy-scale calc-force force active-atom-mask debug-interactions); +CL_DEFMETHOD double ScoringFunction_O::evaluateEnergyForce( NVector_sp pos, + core::T_sp energyScale, + bool calcForce, + NVector_sp force, + core::T_sp activeAtomMask, + core::T_sp debugInteractions ) { double energy; gc::Nilable rawGrad; @@ -124,6 +166,7 @@ CL_DEFMETHOD double ScoringFunction_O::evaluateEnergyForce( NVector_sp pos, bool rawGrad = nil(); } energy = this->evaluateAll(pos, + energyScale, nil(), calcForce, rawGrad, @@ -145,9 +188,10 @@ CL_DEFMETHOD NVector_sp ScoringFunction_O::makeCoordinates() const CL_LISPIFY_NAME("evaluateEnergyForceFullHessian"); -CL_LAMBDA((scoring-function chem:scoring-function) positions calc-force force calc-diagonal-hessian calc-off-diagonal-hessian hessian &optional active-atom-mask debug-interactions); +CL_LAMBDA((scoring-function chem:scoring-function) positions &key energy-scale calc-force force calc-diagonal-hessian calc-off-diagonal-hessian hessian active-atom-mask debug-interactions); CL_DEFMETHOD double ScoringFunction_O::evaluateEnergyForceFullHessian( NVector_sp pos, + core::T_sp energyScale, bool calcForce, NVector_sp force, bool calcDiagonalHessian, bool calcOffDiagonalHessian, @@ -160,6 +204,7 @@ CL_DEFMETHOD double ScoringFunction_O::evaluateEnergyForceFullHessian( if ( calcForce ) rawGrad = force; else rawGrad = nil(); energy = this->evaluateAll( pos, + energyScale, nil(), calcForce, rawGrad, calcDiagonalHessian, @@ -173,8 +218,8 @@ CL_DEFMETHOD double ScoringFunction_O::evaluateEnergyForceFullHessian( } CL_LISPIFY_NAME("evaluateEnergyForceFullHessianForDebugging"); -CL_LAMBDA((scoring-function chem:scoring-function) &optional active-atom-mask debug-interactions); -CL_DEFMETHOD double ScoringFunction_O::evaluateEnergyForceFullHessianForDebugging(core::T_sp activeAtomMask, core::T_sp debugInteractions ) +CL_LAMBDA((scoring-function chem:scoring-function) &key energy-scale active-atom-mask debug-interactions); +CL_DEFMETHOD double ScoringFunction_O::evaluateEnergyForceFullHessianForDebugging(core::T_sp energyScale, core::T_sp activeAtomMask, core::T_sp debugInteractions ) { NVector_sp pos, force; AbstractLargeSquareMatrix_sp hessian; @@ -184,6 +229,7 @@ CL_DEFMETHOD double ScoringFunction_O::evaluateEnergyForceFullHessianForDebuggin hessian = FullLargeSquareMatrix_O::create(this->getNVectorSize(),SymmetricDiagonalLower); this->loadCoordinatesIntoVector(pos); energy = this->evaluateAll(pos, + energyScale, nil(), true, force, @@ -197,6 +243,7 @@ CL_DEFMETHOD double ScoringFunction_O::evaluateEnergyForceFullHessianForDebuggin return energy; } +#if 0 CL_LISPIFY_NAME("calculateEnergy"); CL_LAMBDA((scoring-function chem:scoring-function) &optional active-atom-mask debug-interactions); CL_DEFMETHOD double ScoringFunction_O::calculateEnergy(core::T_sp activeAtomMask, core::T_sp debugInteractions) @@ -209,8 +256,8 @@ CL_DEFMETHOD double ScoringFunction_O::calculateEnergy(core::T_sp activeAtomMask CL_LISPIFY_NAME("calculateEnergyAndForce"); -CL_LAMBDA((scoring-function chem:scoring-function) &optional active-atom-mask debug-interactions); -CL_DEFMETHOD core::T_mv ScoringFunction_O::calculateEnergyAndForce( core::T_sp activeAtomMask, core::T_sp debugInteractions ) +CL_LAMBDA((scoring-function chem:scoring-function) &key energy-scale active-atom-mask debug-interactions); +CL_DEFMETHOD core::T_mv ScoringFunction_O::calculateEnergyAndForce( core::T_sp energyScale, core::T_sp activeAtomMask, core::T_sp debugInteractions ) { NVector_sp pos; NVector_sp force; @@ -218,16 +265,16 @@ CL_DEFMETHOD core::T_mv ScoringFunction_O::calculateEnergyAndForce( core::T_sp a pos = NVector_O::create(this->getNVectorSize()); force = NVector_O::create(this->getNVectorSize()); this->loadCoordinatesIntoVector(pos); - energy = this->evaluateEnergyForce(pos,true,force,activeAtomMask,debugInteractions); + energy = this->evaluateEnergyForce(pos,energyScale,true,force,activeAtomMask,debugInteractions); // To calculate the force magnitude use force->magnitude(); // To calculate the force rmsMagnitude use force->rmsMagnitude(); // this->writeForceToAtoms(force); return Values(core::clasp_make_double_float(energy),force,pos); } +#endif - -CL_LAMBDA((scoring-function chem:scoring-function) coords force &optional (delta 0.00001) active-atom-mask); -CL_DEFMETHOD void ScoringFunction_O::evaluateFiniteDifferenceForce(NVector_sp coords, NVector_sp force, double delta, core::T_sp activeAtomMask ) { +CL_LAMBDA((scoring-function chem:scoring-function) coords &key energy-scale force (delta 0.00001) active-atom-mask); +CL_DEFMETHOD void ScoringFunction_O::evaluateFiniteDifferenceForce(NVector_sp coords, core::T_sp energyScale, NVector_sp force, double delta, core::T_sp activeAtomMask ) { { double x, ylow, yhigh, grad; double deltaDiv2 = delta/2.0; @@ -236,9 +283,9 @@ CL_DEFMETHOD void ScoringFunction_O::evaluateFiniteDifferenceForce(NVector_sp co NVector_sp tpos = NVector_O::make(coords->size(),0.0,false,coords->size(),(Vector_real*)&(*coords)[0]); x = tpos->element(ii); tpos->setElement(ii,x-deltaDiv2); - ylow = this->evaluateEnergy(tpos,activeAtomMask); + ylow = this->evaluateEnergy(tpos,energyScale,nil(),activeAtomMask); tpos->setElement(ii,x+deltaDiv2); - yhigh = this->evaluateEnergy(tpos,activeAtomMask); + yhigh = this->evaluateEnergy(tpos,energyScale,nil(),activeAtomMask); //tpos->setElement(ii,x); // restore grad = (yhigh-ylow)/delta; force->setElement(ii, -grad); @@ -263,6 +310,7 @@ CL_LISPIFY_NAME("velocity-verlet-step"); DOCGROUP(cando); CL_DEFUN void chem__velocity_verlet_step(ScoringFunction_sp scoringFunc, NVector_sp position, + core::T_sp energyScale, NVector_sp velocity, NVector_sp force, NVector_sp force_dt, @@ -295,7 +343,7 @@ CL_DEFUN void chem__velocity_verlet_step(ScoringFunction_sp scoringFunc, atom_idx++; } } - scoringFunc->evaluateEnergyForce(position,true,force_dt,tunfrozen); + scoringFunc->evaluateEnergyForce(position,energyScale,true,force_dt,tunfrozen); atom_idx = 0; for ( size_t idx = 0; idxsize(); idx+=3 ) { if (unfrozen->testBit(atom_idx)==1) { @@ -318,7 +366,7 @@ CL_DEFUN void chem__velocity_verlet_step(ScoringFunction_sp scoringFunc, (*position)[idx+2] = (*position)[idx+2] + offsetz; atom_idx++; } - scoringFunc->evaluateEnergyForce(position,true,force_dt,tunfrozen); + scoringFunc->evaluateEnergyForce(position,energyScale,true,force_dt,tunfrozen); atom_idx = 0; for ( size_t idx = 0; idxsize(); idx+=3 ) { (*velocity)[idx+0] = ((*velocity)[idx+0] + (*delta_t_over_mass)[atom_idx]*0.5*((*force)[idx+0]+(*force_dt)[idx+0]))*scoringFunc->_VelocityScale.getX(); @@ -340,6 +388,7 @@ CL_LISPIFY_NAME("velocity-verlet-step-limit-displacement"); DOCGROUP(cando); CL_DEFUN size_t chem__velocity_verlet_step_limit_displacement(ScoringFunction_sp scoringFunc, NVector_sp position, + core::T_sp energyScale, NVector_sp velocity, NVector_sp force, NVector_sp force_dt, @@ -385,7 +434,7 @@ CL_DEFUN size_t chem__velocity_verlet_step_limit_displacement(ScoringFunction_sp } atom_idx++; } - scoringFunc->evaluateEnergyForce(position,true,force_dt,tunfrozen); + scoringFunc->evaluateEnergyForce(position,energyScale,true,force_dt,tunfrozen); atom_idx = 0; for ( size_t idx = 0; idxsize(); idx+=3 ) { if (!unfrozen || unfrozen->testBit(atom_idx)==1) { diff --git a/src/chem/sketchFunction.cc b/src/chem/sketchFunction.cc index f1097312..860d4f28 100644 --- a/src/chem/sketchFunction.cc +++ b/src/chem/sketchFunction.cc @@ -307,6 +307,7 @@ void maybe_dump_force(const string& msg, NVector_sp force) // double SketchFunction_O::evaluateAll( NVector_sp pos, + core::T_sp energyScale, core::T_sp componentEnergy, bool calcForce, gc::Nilable force, @@ -403,6 +404,7 @@ double SketchFunction_O::evaluateAll( NVector_sp pos, if (this->_Stretch->isEnabled()) totalEnergy += this->_Stretch->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, @@ -418,6 +420,7 @@ double SketchFunction_O::evaluateAll( NVector_sp pos, if (this->_Nonbond->isEnabled()) totalEnergy += this->_Nonbond->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, @@ -436,6 +439,7 @@ double SketchFunction_O::evaluateAll( NVector_sp pos, if (this->_PointToLineRestraint->isEnabled()) totalEnergy += this->_PointToLineRestraint->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, @@ -451,6 +455,7 @@ double SketchFunction_O::evaluateAll( NVector_sp pos, if (this->_OutOfZPlane->isEnabled()) totalEnergy += this->_OutOfZPlane->evaluateAllComponent( this->asSmartPtr(), pos, + energyScale, componentEnergy, calcForce, force, calcDiagonalHessian, @@ -513,33 +518,33 @@ string SketchFunction_O::energyTermsEnabled() #define DELTA 0.00000001 -double SketchFunction_O::calculateNumericalDerivative(NVector_sp pos, double delta, uint i, core::T_sp activeAtomMask ) +double SketchFunction_O::calculateNumericalDerivative(NVector_sp pos, core::T_sp energyScale, double delta, uint i, core::T_sp activeAtomMask ) { double x, ylow, yhigh, fval; double deltaDiv2 = delta/2.0; x = pos->element(i); pos->setElement(i,x-deltaDiv2); - ylow = this->evaluateEnergy(pos,activeAtomMask); + ylow = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x+deltaDiv2); - yhigh = this->evaluateEnergy(pos,activeAtomMask); + yhigh = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x); fval = (yhigh-ylow)/delta; return fval; } -double SketchFunction_O::calculateNumericalSecondDerivative(NVector_sp pos, double delta, uint i, uint j, core::T_sp activeAtomMask ) +double SketchFunction_O::calculateNumericalSecondDerivative(NVector_sp pos, core::T_sp energyScale, double delta, uint i, uint j, core::T_sp activeAtomMask ) { double x, fxmh, fx, fxph, f2; double y, fpipj, fpimj, fmipj, fmimj, fp, fm; if ( i==j ) { x = pos->element(i); pos->setElement(i,x-delta); - fxmh = this->evaluateEnergy(pos,activeAtomMask); + fxmh = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x+delta); - fxph = this->evaluateEnergy(pos,activeAtomMask); + fxph = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x); - fx = this->evaluateEnergy(pos,activeAtomMask); + fx = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); f2 = (fxph+fxmh-2.0*(fx))/(delta*delta); } else { double deltaDiv2 = delta/2.0; @@ -547,16 +552,16 @@ double SketchFunction_O::calculateNumericalSecondDerivative(NVector_sp pos, doub y = pos->element(j); pos->setElement(i,x+deltaDiv2); pos->setElement(j,y+deltaDiv2); - fpipj = this->evaluateEnergy(pos,activeAtomMask); + fpipj = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x+deltaDiv2); pos->setElement(j,y-deltaDiv2); - fpimj = this->evaluateEnergy(pos,activeAtomMask); + fpimj = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x-deltaDiv2); pos->setElement(j,y+deltaDiv2); - fmipj = this->evaluateEnergy(pos,activeAtomMask); + fmipj = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x-deltaDiv2); pos->setElement(j,y-deltaDiv2); - fmimj = this->evaluateEnergy(pos,activeAtomMask); + fmimj = this->evaluateEnergy(pos,energyScale,nil(),activeAtomMask); pos->setElement(i,x); pos->setElement(j,y); LOG("fpipj = {}" , fpipj ); @@ -578,12 +583,12 @@ double SketchFunction_O::calculateNumericalSecondDerivative(NVector_sp pos, doub /*! Calculate the force numerically */ -void SketchFunction_O::evaluateNumericalForce(NVector_sp pos, NVector_sp numForce, double delta, core::T_sp activeAtomMask ) +void SketchFunction_O::evaluateNumericalForce(NVector_sp pos, core::T_sp energyScale, NVector_sp numForce, double delta, core::T_sp activeAtomMask ) { double fval; uint i; for (i=0; isize(); i++ ) { - fval = -this->calculateNumericalDerivative(pos,delta,i,activeAtomMask); + fval = -this->calculateNumericalDerivative(pos,energyScale,delta,i,activeAtomMask); numForce->setElement(i,fval); } } @@ -591,7 +596,7 @@ void SketchFunction_O::evaluateNumericalForce(NVector_sp pos, NVector_sp numForc /*! Calculate the hessian numerically */ -void SketchFunction_O::evaluateNumericalHessian(NVector_sp pos, AbstractLargeSquareMatrix_sp hessian, bool calcOffDiagonal, double delta, core::T_sp activeAtomMask ) +void SketchFunction_O::evaluateNumericalHessian(NVector_sp pos, core::T_sp energyScale, AbstractLargeSquareMatrix_sp hessian, bool calcOffDiagonal, double delta, core::T_sp activeAtomMask ) { double fval; uint c, r; @@ -601,14 +606,14 @@ void SketchFunction_O::evaluateNumericalHessian(NVector_sp pos, AbstractLargeSqu } hessian->zero(); for ( c=0; csize(); c++ ) { - fval = this->calculateNumericalSecondDerivative(pos,delta,c,c,activeAtomMask); + fval = this->calculateNumericalSecondDerivative(pos,energyScale,delta,c,c,activeAtomMask); hessian->setElement(c,c,fval); } if ( !calcOffDiagonal ) return; for ( c=0; csize(); c++ ) { for ( r=0; rsize(); r++ ) { if ( c!=r) { - fval = this->calculateNumericalSecondDerivative(pos,delta,c,r,activeAtomMask); + fval = this->calculateNumericalSecondDerivative(pos,energyScale,delta,c,r,activeAtomMask); hessian->setElement(c,r,fval); } } @@ -623,7 +628,7 @@ void SketchFunction_O::evaluateNumericalHessian(NVector_sp pos, AbstractLargeSqu * If there is a mis-match then dump the SketchFunction into the result. * */ -ForceMatchReport_sp SketchFunction_O::checkIfAnalyticalForceMatchesNumericalForce(NVector_sp pos, NVector_sp analyticalForce, core::T_sp activeAtomMask ) +ForceMatchReport_sp SketchFunction_O::checkIfAnalyticalForceMatchesNumericalForce(NVector_sp pos, core::T_sp energyScale, NVector_sp analyticalForce, core::T_sp activeAtomMask ) { ForceMatchReport_sp report; NVector_sp numForce, tempForce; @@ -634,13 +639,13 @@ ForceMatchReport_sp SketchFunction_O::checkIfAnalyticalForceMatchesNumericalForc report = ForceMatchReport_O::create(); numForce = NVector_O::create(pos->size()); - this->evaluateNumericalForce(pos,numForce,DELTA,activeAtomMask); + this->evaluateNumericalForce(pos,energyScale,numForce,DELTA,activeAtomMask); dot = dotProductWithActiveAtomMask(numForce,analyticalForce,activeAtomMask); numericalMag = magnitudeWithActiveAtomMask(numForce,activeAtomMask); analyticalMag = magnitudeWithActiveAtomMask(analyticalForce,activeAtomMask); tempForce = NVector_O::create(pos->size()); // Evaluate the force at pos again - this->evaluateEnergyForce(pos,true,tempForce,activeAtomMask); + this->evaluateEnergyForce(pos,energyScale,true,tempForce,activeAtomMask); avg = (analyticalMag+numericalMag)/2.0; if ( analyticalMag < VERYSMALL && numericalMag < VERYSMALL ) { result.str(""); @@ -798,11 +803,12 @@ void SketchFunction_O::dealWithProblem(core::Symbol_sp error_symbol, core::T_sp } -CL_LAMBDA(scoring-function position velocity force force-dt delta-t-over-mass delta-t &optional unfrozen) +CL_LAMBDA(scoring-function position &key energy-scale velocity force force-dt delta-t-over-mass delta-t unfrozen) CL_LISPIFY_NAME("sketch-function-velocity-verlet-step"); DOCGROUP(cando); CL_DEFUN void chem__SketchFunction_velocity_verlet_step(SketchFunction_sp sketchFunc, NVector_sp position, + core::T_sp energyScale, NVector_sp velocity, NVector_sp force, NVector_sp force_dt, @@ -838,7 +844,7 @@ CL_DEFUN void chem__SketchFunction_velocity_verlet_step(SketchFunction_sp sketch } atom_idx++; } - sketchFunc->evaluateEnergyForce(position,true,force_dt,tunfrozen); + sketchFunc->evaluateEnergyForce(position,energyScale,true,force_dt,tunfrozen); atom_idx = 0; for ( size_t idx = 0; idxsize(); idx+=3 ) { if (!unfrozen || unfrozen->testBit(atom_idx)==1) { diff --git a/src/kinematics/bondedJoint.cc b/src/kinematics/bondedJoint.cc index 6ea7dc99..32ada9d7 100644 --- a/src/kinematics/bondedJoint.cc +++ b/src/kinematics/bondedJoint.cc @@ -95,6 +95,18 @@ CL_DEFMETHOD Joint_sp BondedJoint_O::inputStubJoint1() const { return this->pare CL_DEFMETHOD Joint_sp BondedJoint_O::inputStubJoint2() const { return this->parent()->parent()->parent(); }; +CL_DEFUN void kin__fillInternalsFromSimpleVectorDouble(BondedJoint_sp joint, core::SimpleVector_double_sp values, size_t index3) { + joint->setDistance((*values)[index3]); + joint->setTheta((*values)[index3+1]); + joint->setPhi((*values)[index3+2]); +} + +CL_DEFUN void kin__fillInternalsFromSimpleVectorSingle(BondedJoint_sp joint, core::SimpleVector_float_sp values, size_t index3) { + joint->setDistance((*values)[index3]); + joint->setTheta((*values)[index3+1]); + joint->setPhi((*values)[index3+2]); +} + CL_LAMBDA(atom-id name atom-table); CL_LISPIFY_NAME("make_BondedJoint"); CL_DEF_CLASS_METHOD diff --git a/src/kinematics/xyzJoint.cc b/src/kinematics/xyzJoint.cc index ba55562d..bd2309e5 100644 --- a/src/kinematics/xyzJoint.cc +++ b/src/kinematics/xyzJoint.cc @@ -73,12 +73,12 @@ CL_DEFMETHOD Joint_sp XyzJoint_O::inputStubJoint1() const { return this->parent( /*! Return the stubJoint3 */ CL_DEFMETHOD Joint_sp XyzJoint_O::inputStubJoint2() const { return this->parent()->parent()->parent(); }; - CL_LAMBDA(atom-id name atom-table &optional atom) CL_LISPIFY_NAME("make_XyzJoint"); CL_DEF_CLASS_METHOD XyzJoint_sp XyzJoint_O::make(const chem::AtomId& atomId, core::T_sp name, chem::AtomTable_sp atomTable, core::T_sp atom) { if (atom.nilp()) { + SIMPLE_ERROR("You must define the atom for XyzJoint creation"); auto obj = gctools::GC::allocate(atomId,name,atomTable); return obj; } else { diff --git a/src/lisp/cando/dynamics.lisp b/src/lisp/cando/dynamics.lisp index c9341c2a..9f963a43 100644 --- a/src/lisp/cando/dynamics.lisp +++ b/src/lisp/cando/dynamics.lisp @@ -81,7 +81,7 @@ Methods are specialized on this class in cando-nglview.lisp.")) for dt-over-m = (/ delta-t mass) do (setf (aref delta-t-over-mass index) (geom:vecreal dt-over-m))) (chem:load-coordinates-into-vector scoring-function coordinates) - (chem:evaluate-energy-force scoring-function coordinates t forces) + (chem:evaluate-energy-force scoring-function coordinates :calc-force t :force forces) (let ((sim (make-instance 'simulation :scoring-function scoring-function :coordinates coordinates @@ -111,7 +111,7 @@ Methods are specialized on this class in cando-nglview.lisp.")) for dt-over-m = delta-t do (setf (aref delta-t-over-mass index) dt-over-m)) (chem:load-coordinates-into-vector scoring-function coordinates) - (chem:evaluate-energy-force scoring-function coordinates t forces) + (chem:evaluate-energy-force scoring-function coordinates :calc-force t :force forces) (let ((sim (make-instance 'simulation :scoring-function scoring-function :coordinates coordinates @@ -146,12 +146,13 @@ Methods are specialized on this class in cando-nglview.lisp.")) (funcall velocity-verlet-function (scoring-function dynamics) (coordinates dynamics) - (velocity dynamics) - (forces dynamics) - (temp-forces dynamics) - (delta-t-over-mass dynamics) - (delta-t dynamics) - unfrozen) + :energy-scale nil + :velocity (velocity dynamics) + :force (forces dynamics) + :force-dt (temp-forces dynamics) + :delta-t-over-mass (delta-t-over-mass dynamics) + :delta-t (delta-t dynamics) + :unfrozen unfrozen) (when (accumulate-coordinates dynamics) (do-accumulate-coordinates dynamics (coordinates dynamics))) (incf (current-time dynamics) (delta-t dynamics))) @@ -164,13 +165,14 @@ Methods are specialized on this class in cando-nglview.lisp.")) (funcall velocity-verlet-function (scoring-function dynamics) (coordinates dynamics) - (velocity dynamics) - (forces dynamics) - (temp-forces dynamics) - (delta-t-over-mass dynamics) - (delta-t dynamics) - unfrozen - limit-displacement) + :energy-scale nil + :velocity (velocity dynamics) + :force (forces dynamics) + :force-dt (temp-forces dynamics) + :delta-t-over-mass (delta-t-over-mass dynamics) + :delta-t (delta-t dynamics) + :tunfrozen unfrozen + :limit-displacement limit-displacement) (when (accumulate-coordinates dynamics) (do-accumulate-coordinates dynamics (coordinates dynamics))) (incf (current-time dynamics) (delta-t dynamics))) diff --git a/src/lisp/cando/molecules.lisp b/src/lisp/cando/molecules.lisp index 24775d13..b0daed94 100644 --- a/src/lisp/cando/molecules.lisp +++ b/src/lisp/cando/molecules.lisp @@ -167,7 +167,7 @@ Example: (set-stereoisomer-mapping *agg* '((:C1 :R) (:C2 :S))" (invoke-restart 'cando:skip-rest-of-minimization err)))) (with-handle-linear-angles-dihedrals (:max-times 3 :verbose verbose) (ext:with-float-traps-masked (:underflow :overflow :invalid :inexact :divide-by-zero) - (chem:minimize minimizer active-atoms-mask)))) + (chem:minimize minimizer :active-atom-mask active-atoms-mask)))) ;; skip-rest-of-minimization can also be triggered by the user from the debugger (skip-rest-of-minimization (err) :report "Skip the rest of the current minimization - continue processing" @@ -288,7 +288,6 @@ Disabling happens before enabling - so you can disable all with T and then selec (chem:enable comp)))) (validate-energy-components disable-components valid-components) (validate-energy-components enable-components valid-components) - (chem:set-electrostatic-scale energy-function electrostatic-scale) (configure-minimizer min :sd-tolerance sd-tolerance :cg-tolerance cg-tolerance @@ -346,7 +345,6 @@ Disabling happens before enabling - so you can disable all with T and then selec (chem:enable comp)))) (validate-energy-components disable-components valid-components) (validate-energy-components enable-components valid-components) - (chem:set-electrostatic-scale energy-function electrostatic-scale) (configure-minimizer min :sd-tolerance sd-tolerance :cg-tolerance cg-tolerance diff --git a/src/lisp/chem-extras/chem.lisp b/src/lisp/chem-extras/chem.lisp index 1550844a..7c36cfb2 100644 --- a/src/lisp/chem-extras/chem.lisp +++ b/src/lisp/chem-extras/chem.lisp @@ -343,7 +343,7 @@ evaluates the body in that dynamic environment." (setf *save-pos* pos) (with-output-to-string (sout) (let* ((energy-components (chem:make-energy-components)) - (total-energy (chem:evaluate-all energy-function pos energy-components))) + (total-energy (chem:evaluate-all energy-function pos :component-energy energy-components))) (format sout "Total E. ~10,5f mask ~a~%" total-energy (not (null active-atom-mask))) (format sout "~s~%" (chem:energy-components/component-energies energy-components))))) #|| diff --git a/src/lisp/regression-tests/energy.lisp b/src/lisp/regression-tests/energy.lisp index 980b008d..9f28ea69 100644 --- a/src/lisp/regression-tests/energy.lisp +++ b/src/lisp/regression-tests/energy.lisp @@ -18,9 +18,9 @@ (format t "pos = ~s~%" pos) #+(or)(progn ;; (gctools:wait-for-user-signal "signal") - (time (dotimes (i 1) (defparameter energy (chem:evaluate-energy-force ef pos t force))))) + (time (dotimes (i 1) (defparameter energy (chem:evaluate-energy-force ef pos :calc-force t :force force))))) (format t "About to calculate energy~%") -(defparameter energy (chem:evaluate-energy-force ef pos t force)) +(defparameter energy (chem:evaluate-energy-force ef pos :calc-force t :force force)) (defparameter expected-energy 20937032.6233387) (format t " energy = ~f~%" energy) (format t "expected-energy = ~f~%" expected-energy) @@ -133,7 +133,7 @@ (defun energy-components (name en pos mask) (let ((components (chem:all-components en)) (ht (make-hash-table))) - (format t "~s energy total = ~f~%" name (chem:evaluate-energy en pos mask)) + (format t "~s energy total = ~f~%" name (chem:evaluate-energy en pos :active-atom-mask mask)) (loop for comp in components for energy = (chem:energy-component-evaluate-energy en comp pos mask) for class-name = (class-name (class-of comp)) diff --git a/src/lisp/sketch2d/workbench.lisp b/src/lisp/sketch2d/workbench.lisp index b300686a..a232621f 100644 --- a/src/lisp/sketch2d/workbench.lisp +++ b/src/lisp/sketch2d/workbench.lisp @@ -27,7 +27,7 @@ (loop for x from 0.5 below 100.0 by 0.5 for y = (progn (setf (elt *v* 0) x) - (chem:evaluate-energy-force *sk* *v* t *f* )) + (chem:evaluate-energy-force *sk* *v* :calc-force t :force *f* )) do (format t "~a ~a~%" x y))) ) diff --git a/src/lisp/topology/assembler.lisp b/src/lisp/topology/assembler.lisp index 10f1f81d..028f4158 100644 --- a/src/lisp/topology/assembler.lisp +++ b/src/lisp/topology/assembler.lisp @@ -2,12 +2,13 @@ (defclass orientation () - ((to-origin :initarg :to-origin :accessor to-origin + ((to-origin :initarg :to-origin :reader to-origin :documentation "This is the transform that takes the object to the origin") - (adjust :initarg :adjust :accessor adjust + (adjust :initarg :adjust :reader adjust :documentation "Apply this to the coordinates after the to-origin transform was applied") - (from-origin :initarg :from-origin :accessor from-origin - :documentation "This is the second transform that takes the object after to-origin is applied")) + (from-origin :initarg :from-origin :reader from-origin + :documentation "This is the second transform that takes the object after to-origin is applied") + (transform-cache :initarg :transform-cache :reader transform-cache)) (:documentation "Represent the orientation of a molecule in a design calculation. There are three transforms that define the orientation. The `to-origin` transform moves something on the molecule onto the origin. @@ -26,10 +27,12 @@ default to the identity matrix." (error "make-orientation adjust must be non-nil")) (unless to-origin (error "make-orientation to-origin must be non-nil")) - (make-instance 'orientation - :from-origin from-origin - :adjust adjust - :to-origin to-origin)) + (let ((transform (calculate-orientation-transform from-origin adjust to-origin))) + (make-instance 'orientation + :from-origin from-origin + :adjust adjust + :to-origin to-origin + :transform-cache transform))) (defun copy-orientation (orientation) "Copy an ORIENTATION." @@ -46,23 +49,19 @@ default to the identity matrix." (defgeneric kin:orientation-transform (orientation)) -(defmethod kin:orientation-transform ((orientation orientation)) - "Combine the components of the ORIENTATION into a 4x4 homogeneous matrix." - (let* ((from-origin (topology:from-origin orientation)) - (adjust (adjust orientation)) - (to-origin (to-origin orientation)) - (transform0 (geom:m*m adjust to-origin)) +(defun calculate-orientation-transform (from-origin adjust to-origin) + (let* ((transform0 (geom:m*m adjust to-origin)) (transform (geom:m*m from-origin transform0))) transform)) +(defmethod kin:orientation-transform ((orientation orientation)) + "Combine the components of the ORIENTATION into a 4x4 homogeneous matrix." + (transform-cache orientation)) + (defclass orientations () ((orientations :initform (make-hash-table) :initarg :orientations :accessor orientations)) (:documentation "A hash-table that maps OLIGOMER-SHAPEs to ORIENTATIONs")) -(defun orientation-for-oligomer-shape (assembler oligomer-shape) - "Return the orientation for the OLIGOMER-SHAPE in the ASSEMBLER" - (gethash oligomer-shape (orientations (orientations assembler)))) - (defun ensure-complete-orientations (orientations oligomer-shapes) (loop for oligomer-shape in oligomer-shapes unless (gethash oligomer-shape (orientations orientations)) @@ -81,6 +80,13 @@ default to the identity matrix." do (setf (gethash oligomer-shape ht) orientation)) (make-instance 'orientations :orientations ht))) +(defun oligomer-space-orientations (orientations) + "Convert the orientations hash-table to be keyed on oligomer-space, rather than oligomer or oligomer-shape" + (let ((ht (make-hash-table))) + (maphash (lambda (key value) + (setf (gethash (oligomer-space key) ht) value)) + (orientations orientations)) + (make-instance 'orientations :orientations ht))) (defclass monomer-position () ((molecule-index :initarg :molecule-index :reader molecule-index) @@ -150,9 +156,18 @@ The most important functions are UPDATE-INTERNALS and UPDATE-EXTERNALS." (receptor-oligomer-shape-orientation assembler) (values ligand-oligomer-shape ligand-orientation receptor-oligomer-shape receptor-orientation)))) -(defun lookup-orientation (assembler oligomer-shape) - "Return the orientation for the OLIGOMER-SHAPE in the ASSEMBLER." - (gethash oligomer-shape (orientations (orientations assembler)))) +(defgeneric lookup-orientation (assembler-or-orientations oligomer-thing) + (:documentation "Return the orientation for the OLIGOMER-THING in the ASSEMBLER-OR-ORIENTATION")) + +(defmethod lookup-orientation ((assembler assembler) oligomer-thing) + "Return the orientation for the OLIGOMER-THING (oligomer-space, oligomer, or oligomer-shape) in the ASSEMBLER." + (or (gethash oligomer-thing (orientations (orientations assembler))) + (error "Could not find ~s as a key in ~s" oligomer-thing assembler))) + + +(defmethod lookup-orientation ((orientations orientations) oligomer-thing) + (or (gethash oligomer-thing (orientations orientations)) + (error "Could not find ~s as a key in ~s" oligomer-thing orientations))) (defclass subset-assembler (assembler) () @@ -209,6 +224,18 @@ The most important functions are UPDATE-INTERNALS and UPDATE-EXTERNALS." (error "Could not find ~s in ~s" monomer assembler) nil)) +(defgeneric orientation-for-oligomer-shape (thing oligomer-shape)) + +(defmethod orientation-for-oligomer-shape ((assembler assembler) oligomer-shape) + "Return the orientation for the OLIGOMER-SHAPE in the ASSEMBLER" + (or (gethash oligomer-shape (orientations (orientations assembler))) + (error "Could not find the orientation for ~s in ~s" oligomer-shape assembler))) + +(defmethod orientation-for-oligomer-shape ((orientations orientations) oligomer-shape) + "Return the orientation for the OLIGOMER-SHAPE in the ORIENTATIONS" + (or (gethash oligomer-shape (orientations orientations)) + (error "Could not find the orientation for ~s in ~s" oligomer-shape orientations))) + (defun make-coordinates-for-number-of-atoms (number-of-atoms) (make-array (* 3 number-of-atoms) :element-type (geom:vecreal-type) :initial-element (geom:vecreal 0.0))) @@ -263,13 +290,13 @@ Specialize the foldamer argument to provide methods")) :coordinates coordinates :monomers ht))) -(defun in-monomer-subset (monomer-subset &rest monomers) +(defun in-monomer-subset (monomer-subset monomer) (cond ((null monomer-subset) ;; If monomer-subset is NULL then everything is in the subset t) ((typep monomer-subset 'monomer-subset) - (every (lambda (mon) (gethash mon (monomers monomer-subset))) monomers)) + (gethash monomer (monomers monomer-subset))) (t (error "Illegal value for monomer-subset ~s - must be NIL or a hash-table")))) (defun maybe-update-backbone-dihedral-cache (assembler) @@ -308,15 +335,13 @@ Specialize the foldamer argument to provide methods")) (setf (backbone-dihedral-cache-deg monomer-shape) dihedral-cache) )))))) -(defun make-assembler (oligomer-shapes &key orientations monomer-subset energy-function-factory) +(defun make-assembler (oligomer-shapes &key (orientations nil orientations-p) monomer-subset energy-function-factory (monomer-contexts nil monomer-contexts-p)) "Build a assembler for the OLIGOMER-SHAPES. - OLIGOMER-SHAPES - A list of OLIGOMER-SHAPEs that the ASSEMBLER will build. +MONOMER-CONTEXTS - A map of monomers to monomer-contexts copied from another assembler (avoids recalculating them). ORIENTATIONS - An ORIENTATIONS object that maps OLIGOMER-SHAPES to ORIENTATIONs. USE-EXCLUDED-ATOMS - A parameter passed to make-energy-function. -ENERGY-FUNCTION-FACTORY - If defined, call this with the aggregate to make the energy-function. -/*KEEP-INTERACTION-FACTORY-FACTORY - T or a lambda that takes a list of (cons oligomer-shape molecule) and returns a keep-interaction-factory function for the energy-function.*/ -/*TUNE-ENERGY-FUNCTION - A function that takes the energy-function and an assembler and modifies the energy-function.*/" +ENERGY-FUNCTION-FACTORY - If defined, call this with the aggregate to make the energy-function." (cond ((not (or (= (length oligomer-shapes) 1) orientations)) (error "You must provide orientations when there is more than one oligomer-shape")) @@ -348,10 +373,12 @@ ENERGY-FUNCTION-FACTORY - If defined, call this with the aggregate to make the e :monomer-positions-accumulator monomer-positions) do (chem:setf-force-field-name molecule (oligomer-force-field-name foldamer)) collect (cons oligomer-shape molecule)))) - (let* ((monomer-contexts (make-hash-table)) - (ataggregate (let ((atagg (make-instance 'ataggregate :aggregate aggregate))) + (let* ((ataggregate (let ((atagg (make-instance 'ataggregate :aggregate aggregate))) (resize-atmolecules atagg (length oligomer-shapes)) atagg)) + (new-monomer-contexts (if monomer-contexts-p + monomer-contexts + (make-hash-table))) (joint-tree (make-joint-tree)) (adjustments (make-instance 'adjustments)) (energy-function (if energy-function-factory @@ -365,14 +392,18 @@ ENERGY-FUNCTION-FACTORY - If defined, call this with the aggregate to make the e for foldamer = (foldamer oligomer-space) for molecule-index from 0 for monomer-to-monomer-shape-map = (monomer-shape-map oligomer-shape) - do (if monomer-subset - (loop for monomer in (alexandria:hash-table-keys (topology:monomers monomer-subset)) - when (topology:monomer-position monomer (topology:oligomer-space oligomer)) - do (let ((monomer-context (foldamer-monomer-context monomer oligomer foldamer))) - (setf (gethash monomer monomer-contexts) monomer-context))) - (loop for monomer across (monomers oligomer-space) - for monomer-context = (foldamer-monomer-context monomer oligomer foldamer) - do (setf (gethash monomer monomer-contexts) monomer-context))) + do (cond + (monomer-contexts-p nil) + (monomer-subset + (let ((monomer-contexts (make-hash-table))) + (loop for monomer in (alexandria:hash-table-keys (topology:monomers monomer-subset)) + when (topology:monomer-position monomer (topology:oligomer-space oligomer)) + do (let ((monomer-context (foldamer-monomer-context monomer oligomer foldamer))) + (setf (gethash monomer new-monomer-contexts) monomer-context))))) + (t + (loop for monomer across (monomers oligomer-space) + for monomer-context = (foldamer-monomer-context monomer oligomer foldamer) + do (setf (gethash monomer new-monomer-contexts) monomer-context)))) ;. This is where I would invoke Conformation_O::buildMoleculeUsingOligomer ;; Use the monomers-to-topologys do (let* ((atmolecule (build-atmolecule-using-oligomer oligomer @@ -389,7 +420,7 @@ ENERGY-FUNCTION-FACTORY - If defined, call this with the aggregate to make the e 'subset-assembler 'assembler) :monomer-positions monomer-positions - :monomer-contexts monomer-contexts + :monomer-contexts new-monomer-contexts :oligomer-shapes oligomer-shapes :orientations orientations :aggregate aggregate @@ -582,25 +613,6 @@ ENERGY-FUNCTION-FACTORY - If defined, call this with the aggregate to make the e (declare (ignore atom-id)) (kin:update-internal-coord joint coordinates))))) -(defun build-external-coordinates (assembler &key oligomer-shape - (orientation :identity orientationp) - (coords (topology:make-coordinates-for-assembler assembler))) - (if oligomer-shape - (progn - (when (and oligomer-shape (not orientationp)) - (error "You must provide orientation when you provide oligomer-shape")) - (build-atom-tree-external-coordinates* assembler coords oligomer-shape orientation) - (adjust-atom-tree-external-coordinates assembler coords oligomer-shape)) - (loop for oligomer-shape in (oligomer-shapes assembler) - do (build-external-coordinates assembler :oligomer-shape oligomer-shape - :orientation oligomer-shape - :coords coords)))) - -(defun update-externals (assembler &key oligomer-shape - (orientation :identity orientationp) - (coords (topology:make-coordinates-for-assembler assembler))) - (build-external-coordinates assembler :oligomer-shape oligomer-shape :orientation orientation :coords coords)) - (defun build-atom-tree-external-coordinates* (assembler coords oligomer-shape maybe-orientation) (let* ((orientation (orientation maybe-orientation assembler)) @@ -612,12 +624,16 @@ ENERGY-FUNCTION-FACTORY - If defined, call this with the aggregate to make the e (loop for joint in joints do (kin:update-xyz-coords joint coords))))) -(defun build-all-external-coordinates (assembler &key (coords (topology:make-coordinates-for-assembler assembler))) - (loop for oligomer-shape in (oligomer-shapes assembler) - do (build-external-coordinates assembler - :coords coords - :oligomer-shape oligomer-shape - :orientation oligomer-shape))) + +(defun build-atom-tree-for-monomer-shape-external-coordinates* (assembler coords oligomer-shape monomer-shape maybe-orientation) + (let* ((orientation (orientation maybe-orientation assembler)) + (one-oligomer (oligomer oligomer-shape)) + (joints (gethash one-oligomer (root-map (joint-tree assembler))))) + (when (null joints) + (error "Could not find oligomer ~s in root-map ~s" one-oligomer (root-map (joint-tree assembler)))) + (with-orientation orientation + (loop for joint in joints + do (kin:update-xyz-coords joint coords))))) (defun adjust-atom-tree-external-coordinates (assembler coords oligomer-shape) @@ -669,21 +685,24 @@ ENERGY-FUNCTION-FACTORY - If defined, call this with the aggregate to make the e (loop for oligomer-shape in (oligomer-shapes assembler) do (adjust-internals assembler oligomer-shape))) -(defgeneric apply-monomer-shape-to-atresidue-internals (assembler oligomer-shape monomer-shape monomer-context atresidue coordinates &key verbose) +(defgeneric apply-monomer-shape-to-atresidue-internals (assembler oligomer-shape monomer-shape monomer-context atresidue &key verbose) (:documentation "Specialize this on different monomer-shape classes. Fill in internal coordinates into the atresidue for the monomer-shape. Some specialized methods will need coordinates for the assembler")) - -(defun fill-internals-from-oligomer-shape (assembler oligomer-shape &optional verbose) - "Fill internal coordinates from the fragments" +#+(or) +(defun fill-internals-from-oligomer-shape (assembler oligomer-shape + &key (root-monomer (topology:root-monomer (topology:oligomer oligomer-shape)) root-monomer-p) + verbose) + "Fill internal coordinates from the monomer-shapes in OLIGOMER-SHAPE of the ASSEMBLER. +If ROOT-MONOMER is provided, then start from that." (when verbose (let ((*print-pretty* nil)) (format t "fill-internals-from-oligomer-shape ~s~%" oligomer-shape))) (let ((coordinates (topology:make-coordinates-for-assembler assembler))) (loop for ass-oligomer-shape in (oligomer-shapes assembler) when (eq ass-oligomer-shape oligomer-shape) do (let* ((oligomer (oligomer ass-oligomer-shape)) - (ordered-monomers (ordered-monomers oligomer))) + (ordered-monomers (ordered-monomers oligomer :root-monomer root-monomer))) (when verbose (let ((*print-pretty* nil)) (format t "fill-internals-from-oligomer-shape ~s ordered-monomers: ~s~%" ass-oligomer-shape ordered-monomers))) @@ -705,6 +724,62 @@ Some specialized methods will need coordinates for the assembler")) (when verbose (format t "applying internals for monomer: ~s~%" monomer)) (apply-monomer-shape-to-atresidue-internals assembler oligomer-shape monomer-shape monomer-context atres coordinates :verbose verbose))))))) +(defun fill-internals-from-oligomer-shape (assembler oligomer-shape &key verbose) + "Fill internal coordinates from the monomer-shapes in OLIGOMER-SHAPE of the ASSEMBLER." + (when verbose (let ((*print-pretty* nil)) (format t "fill-internals-from-oligomer-shape ~s~%" oligomer-shape))) + (loop for ass-oligomer-shape in (oligomer-shapes assembler) + when (eq ass-oligomer-shape oligomer-shape) + do (let* ((oligomer (oligomer ass-oligomer-shape)) + (monomers (monomers (oligomer-space oligomer)))) + (loop with atagg = (ataggregate assembler) + ;; IGNORING THE FOLLOWING + ;; It's really important that we use the ordered-monomers so that the monomer-shapes + ;; will install internal coordinates in the order from the root outwards. + ;; so that any preceeding monomer-shapes are built before any following ones. + for monomer across monomers + when (in-monomer-subset (monomer-subset assembler) monomer) + do (let* ((monomer-context (gethash monomer (monomer-contexts assembler))) + (monomer-position (gethash monomer (monomer-positions assembler))) + (molecule-index (molecule-index monomer-position)) + (residue-index (residue-index monomer-position)) + (atmol (elt (atmolecules atagg) molecule-index)) + (atres (elt (atresidues atmol) residue-index)) + (monomer-shape (let ((ms (gethash monomer (monomer-shape-map oligomer-shape)))) + (unless ms (error "Could not get monomer-shape for monomer ~a" monomer)) + ms))) + (when verbose (format t "applying internals for monomer: ~s~%" monomer)) + (apply-monomer-shape-to-atresidue-internals assembler oligomer-shape monomer-shape monomer-context atres :verbose verbose)))))) + +(defun update-internals (assembler oligomer-shape &key verbose) + "Update the internal coordinates of the assembler for the oligomer-shape +or whatever type the second argument is. +The MONOMER-SHAPEs of the OLIGOMER-SHAPE are used to update the internal coordinates in +the atom-tree. + +ASSEMBLER - the assembler to update. +OLIGOMER-SHAPE - An oligomer-shape (or permissible-rotamers - I think this is wrong) in the assembler." + (fill-internals-from-oligomer-shape assembler oligomer-shape :verbose verbose) + (topology:with-orientation (topology:lookup-orientation assembler oligomer-shape) + (adjust-internals assembler oligomer-shape))) + +(defun update-internals-for-monomer-shape (assembler oligomer-shape monomer-shape &key verbose) + "Fill internal coordinates from the MONOMER-SHAPE in OLIGOMER-SHAPE of the ASSEMBLER." + (let* ((atagg (ataggregate assembler)) + (monomer-shape-pos (or (gethash monomer-shape (monomer-shape-to-index oligomer-shape)) + (error "Could not find monomer-shape ~s in oligomer-shape ~s" monomer-shape oligomer-shape))) + (monomer-info (aref (monomer-shape-info-vector oligomer-shape) monomer-shape-pos)) + (monomer (monomer monomer-info)) + (monomer-context (gethash monomer (monomer-contexts assembler))) + (monomer-position (gethash monomer (monomer-positions assembler))) + (molecule-index (molecule-index monomer-position)) + (residue-index (residue-index monomer-position)) + (atmol (elt (atmolecules atagg) molecule-index)) + (atres (elt (atresidues atmol) residue-index))) + (when verbose (format t "applying internals for monomer: ~s~%" monomer)) + (apply-monomer-shape-to-atresidue-internals assembler oligomer-shape monomer-shape monomer-context atres :verbose verbose)) + (topology:with-orientation (topology:lookup-orientation assembler oligomer-shape) + (adjust-internals assembler oligomer-shape))) + #| ;;;Idea to make monomer-shape subclasses control how internal coordinates get generated. @@ -897,3 +972,42 @@ Some specialized methods will need coordinates for the assembler")) do (incf centerz z) do (incf num)) (geom:vec (/ centerx num) (/ centery num) (/ centerz num)))) + +(defun update-externals (assembler &key oligomer-shape + (orientation :identity orientationp) + (coords (topology:make-coordinates-for-assembler assembler))) + "Update the external coordinates in COORDS using the ASSEMBLER and ORIENTATION. +IF OLIGOMER-SHAPE is provided then just build externals for that OLIGOMER-SHAPE. +If OLIGOMER-SHAPE is not provided then build them all and use the OLIGOMER-SHAPE as the ORIENTATION key. +Return the COORDS." + (if oligomer-shape + (progn + (when (and oligomer-shape (not orientationp)) + (error "You must provide orientation when you provide oligomer-shape")) + (build-atom-tree-external-coordinates* assembler coords oligomer-shape orientation) + (adjust-atom-tree-external-coordinates assembler coords oligomer-shape)) + (loop for oligomer-shape in (oligomer-shapes assembler) + do (update-externals assembler :oligomer-shape oligomer-shape + :orientation oligomer-shape + :coords coords))) + coords) + + +(defun update-externals-for-monomer-shape (assembler oligomer-shape monomer-shape coords orientation) + "Update the external coordinates for one MONOMER-SHAPE in the OLIGOMER-SHAPE in the ASSEMBLER using the ORIENTATION." + (let* ((atagg (ataggregate assembler)) + (monomer-shape-pos (or (gethash monomer-shape (monomer-shape-to-index oligomer-shape)) + (error "Could not find monomer-shape ~s in oligomer-shape ~s" monomer-shape oligomer-shape))) + (monomer-info (aref (monomer-shape-info-vector oligomer-shape) monomer-shape-pos)) + (monomer (monomer monomer-info)) + (monomer-context (gethash monomer (monomer-contexts assembler))) + (monomer-position (gethash monomer (monomer-positions assembler))) + (molecule-index (molecule-index monomer-position)) + (residue-index (residue-index monomer-position)) + (atmol (elt (atmolecules atagg) molecule-index)) + (atres (elt (atresidues atmol) residue-index)) + (joints (joints atres)) + (joint0 (elt joints 0))) + (with-orientation orientation + (kin:update-xyz-coords joint0 coords)) + (adjust-atom-tree-external-coordinates assembler coords oligomer-shape))) diff --git a/src/lisp/topology/internals.lisp b/src/lisp/topology/internals.lisp index f298c360..52f78891 100644 --- a/src/lisp/topology/internals.lisp +++ b/src/lisp/topology/internals.lisp @@ -370,15 +370,11 @@ changing the BACKBONE-DIHEDRAL-CACHE." (defgeneric apply-fragment-internals-to-atresidue (fragment-internals rotamer-index atresidue)) -(defmethod apply-fragment-internals-to-atresidue ((obj rotamer) rotamer-index atresidue) +(defmethod apply-fragment-internals-to-atresidue ((rotamer rotamer) rotamer-index atresidue) (setf (rotamer-index atresidue) rotamer-index) (loop for joint across (topology:joints atresidue) - for index from 0 - do (multiple-value-bind (bond angle-rad dihedral-rad) - (topology:extract-bond-angle-rad-dihedral-rad obj index) - (when (typep joint 'kin:bonded-joint) - (topology:fill-joint-internals joint bond angle-rad dihedral-rad))))) - + for index3 from 0 by 3 + do (fill-joint-internals-from-vector joint (internals-values rotamer) index3))) (defmethod monomer-context-to-context-rotamers ((obj rotamers-database)) (context-to-rotamers obj)) @@ -813,7 +809,6 @@ No checking is done to make sure that the list of clusterable-context-rotamers a (format finternals "end-conformation~%") )) - (defgeneric fill-joint-internals (joint bond angle-rad dihedral-rad)) (defmethod fill-joint-internals ((joint kin:jump-joint) bond angle-rad dihedral-rad) @@ -823,9 +818,33 @@ No checking is done to make sure that the list of clusterable-context-rotamers a "Do nothing with xyz-joints" ) +(defgeneric fill-joint-internals-from-vector (joint vector index3)) + +(defmethod fill-joint-internals-from-vector ((joint kin:jump-joint) vector index3) #|do nothing|#) +(defmethod fill-joint-internals-from-vector ((joint kin:xyz-joint) vector index3) #|do nothing|#) + +(defmethod fill-joint-internals-from-vector ((joint kin:bonded-joint) vector index3) + (multiple-value-bind (bond angle-rad dihedral-rad) + (values (aref vector index3) + (aref vector (+ 1 index3)) + (aref vector (+ 2 index3))) + (kin:set-distance joint bond) + (kin:set-theta joint angle-rad) + (kin:bonded-joint/set-phi joint dihedral-rad))) + +(defmethod fill-joint-internals-from-vector ((joint kin:bonded-joint) (vector sys:simple-vector-double) index3) + (kin:fill-internals-from-simple-vector-double joint vector index3)) + +(defmethod fill-joint-internals-from-vector ((joint kin:bonded-joint) (vector sys:simple-vector-float) index3) + (kin:fill-internals-from-simple-vector-single joint vector index3)) + + (defun fill-joint-phi (joint phi) + "Write a dihedral angle in radians in PHI into the JOINT. +I have this to trap writing specific phi values for debugging." (kin:bonded-joint/set-phi joint phi)) +#+(or) (defmethod fill-joint-internals ((joint kin:bonded-joint) bond angle-rad dihedral-rad) (kin:set-distance joint bond) (kin:set-theta joint angle-rad) @@ -833,6 +852,8 @@ No checking is done to make sure that the list of clusterable-context-rotamers a ) (defun extract-bond-angle-rad-dihedral-rad (rotamer index) +"Extract the bond, angle-rad, and dihedral-rad from the ROTAMER at INDEX. +INDEX counts by 1 and _not_ by 3." (let ((index3 (* 3 index))) (values (aref (internals-values rotamer) index3) (aref (internals-values rotamer) (+ 1 index3)) @@ -842,10 +863,8 @@ No checking is done to make sure that the list of clusterable-context-rotamers a (defmethod apply-fragment-internals-to-atresidue ((fragment-internals fragment-internals) rotamer-index atresidue) (setf (rotamer-index atresidue) rotamer-index) (loop for joint across (joints atresidue) - for index from 0 - do (multiple-value-bind (bond angle-rad dihedral-rad) - (extract-bond-angle-rad-dihedral-rad fragment-internals index) - (fill-joint-internals joint bond angle-rad dihedral-rad)))) + for index3 from 0 by 3 + do (fill-joint-internals-from-vector joint (internals-values fragment-internals) index3))) (defun analyze-atresidue (atresidue) (let ((internals-count 0) @@ -881,13 +900,12 @@ No checking is done to make sure that the list of clusterable-context-rotamers a for index from 0 for index3 from 0 by 3 if (or (null updated-internals-mask) - (= (aref updated-internals-mask index) 1)) + (= (aref updated-internals-mask index) 1)) do (typecase joint (kin:bonded-joint (when verbose (format t "Updating bonded-joint ~s~%" joint)) - (fill-joint-internals joint (aref internals index3) - (aref internals (+ 1 index3)) - (aref internals (+ 2 index3)))) + (error "An error is about to occur - you need to pass the internals-values of internals and not internals") + (fill-joint-internals-from-vector joint internals index3)) (kin:xyz-joint (if (kin:xyz-joint/definedp joint) (warn "Skipping defining xyz-joint") diff --git a/src/lisp/topology/joint-templates.lisp b/src/lisp/topology/joint-templates.lisp index 5a57ea55..be4cdd61 100644 --- a/src/lisp/topology/joint-templates.lisp +++ b/src/lisp/topology/joint-templates.lisp @@ -493,7 +493,7 @@ (let* ((constitution-atoms-index (constitution-atoms-index joint-template)) (atom-name (atom-name joint-template)) (atomid (list atmolecule-index atresidue-index constitution-atoms-index)) - (new-joint (kin:make-xyz-joint atomid atom-name atom-table))) + (new-joint (kin:make-xyz-joint atomid atom-name atom-table (chem:make-atom atom-name :DU)))) (put-joint atresidue new-joint constitution-atoms-index) (when parent-joint (kin:joint/add-child parent-joint new-joint)) new-joint))) diff --git a/src/lisp/topology/jupyter.lisp b/src/lisp/topology/jupyter.lisp index 7d4bc210..38fc680a 100644 --- a/src/lisp/topology/jupyter.lisp +++ b/src/lisp/topology/jupyter.lisp @@ -158,7 +158,7 @@ (let* ((assembler (topology:make-assembler (list object))) (coords (topology:make-coordinates-for-assembler assembler))) (topology:update-internals assembler object) - (topology:build-external-coordinates assembler :coords coords + (topology:update-externals assembler :coords coords :oligomer-shape object :orientation :identity) (cando-widgets::show-on-pane pane-instance (cando:mol (topology:aggregate* assembler coords) 0)))) diff --git a/src/lisp/topology/linearize-internals.lisp b/src/lisp/topology/linearize-internals.lisp index 295380d4..f2fc9400 100644 --- a/src/lisp/topology/linearize-internals.lisp +++ b/src/lisp/topology/linearize-internals.lisp @@ -614,8 +614,5 @@ monomer-context-index indexes into ... (defmethod apply-fragment-internals-to-atresidue ((fragment-internals linearized-fragment-internals) rotamer-index atresidue) (setf (rotamer-index atresidue) rotamer-index) (loop for joint across (joints atresidue) - for index from 0 by 3 - for bond = (elt (internals-values fragment-internals) index) - for angle-rad = (elt (internals-values fragment-internals) (+ 1 index)) - for dihedral-rad = (elt (internals-values fragment-internals) (+ 2 index)) - do (fill-joint-internals joint bond angle-rad dihedral-rad))) + for index3 from 0 by 3 + do (fill-internals-from-vector joint internals-values index3))) diff --git a/src/lisp/topology/oligomer.lisp b/src/lisp/topology/oligomer.lisp index 98eb25c1..3630d774 100644 --- a/src/lisp/topology/oligomer.lisp +++ b/src/lisp/topology/oligomer.lisp @@ -177,18 +177,17 @@ unique-ring-couplings))))))) result)) -(defun ordered-monomers (oligomer) - "Return a canonical sequence of monomers" - (let ((root-monomer (root-monomer oligomer))) - (let ((monomer-out-couplings (make-hash-table)) - (unique-ring-couplings nil)) - (loop for index below (length (couplings oligomer)) - for coupling = (elt (couplings oligomer) index) - if (typep coupling 'directional-coupling) - do (push coupling (gethash (source-monomer coupling) monomer-out-couplings)) - else - do (pushnew coupling unique-ring-couplings)) - (ordered-sequence-monomers oligomer nil root-monomer monomer-out-couplings unique-ring-couplings)))) +(defun ordered-monomers (oligomer &key (root-monomer (root-monomer oligomer))) + "Return a canonical sequence of monomers in the OLIGOMER. If ROOT-MONOMER is provided then start from that." + (let ((monomer-out-couplings (make-hash-table)) + (unique-ring-couplings nil)) + (loop for index below (length (couplings oligomer)) + for coupling = (elt (couplings oligomer) index) + if (typep coupling 'directional-coupling) + do (push coupling (gethash (source-monomer coupling) monomer-out-couplings)) + else + do (pushnew coupling unique-ring-couplings)) + (ordered-sequence-monomers oligomer nil root-monomer monomer-out-couplings unique-ring-couplings))) (defun build-residue (oligomer diff --git a/src/lisp/topology/packages.lisp b/src/lisp/topology/packages.lisp index ee515290..19cd8ebc 100644 --- a/src/lisp/topology/packages.lisp +++ b/src/lisp/topology/packages.lisp @@ -373,12 +373,9 @@ #:assembler-base #:oligomer-shapes #:walk-atmolecule-joints - #:build-all-oligomer-shapes-in-coordinates - #:build-external-coordinates #:aggregate* #:make-monomer-subset - #:update-internals-for-atresidue #:update-internals #:update-externals #:build-oligomer-shape-externals @@ -486,7 +483,13 @@ #:validate-indexes-for-shape-keys #:lookup-orientation #:find-atom-for-joint - #:monomer-position)) + #:monomer-position + #:monomer-contexts + #:oligomer-space-contains-monomer + #:orientation-for-oligomer-shape + #:oligomer-space-orientations + #:update-internals-for-monomer-shape + #:update-externals-for-monomer-shape)) (defpackage #:topology.dag (:use #:common-lisp) diff --git a/src/lisp/topology/shape.lisp b/src/lisp/topology/shape.lisp index 43b3c841..f39f1af3 100644 --- a/src/lisp/topology/shape.lisp +++ b/src/lisp/topology/shape.lisp @@ -56,7 +56,7 @@ If ORIGINAL-ROTAMER-SHAPE is defined then it must be a ROTAMER-SHAPE and we copy (context-rotamers (topology:lookup-rotamers-for-context rotamers-database monomer-context))) (make-instance 'rotamer-shape :rotamer-shape-context-rotamers context-rotamers))) -(defmethod apply-monomer-shape-to-atresidue-internals (assembler oligomer-shape (rotamer-shape rotamer-shape) monomer-context atresidue coordinates &key verbose) +(defmethod apply-monomer-shape-to-atresidue-internals (assembler oligomer-shape (rotamer-shape rotamer-shape) monomer-context atresidue &key verbose) (when verbose (let ((*print-pretty* nil)) (format t "apply-monomer-shape-to-atresidue-internals ~% oligomer-shape ~s~% rotamer-shape ~s~% monomer-context ~s~%" @@ -92,7 +92,7 @@ If ORIGINAL-ROTAMER-SHAPE is defined then it must be a ROTAMER-SHAPE and we copy (format stream "~s" (chem:get-name (residue obj))))) -(defmethod apply-monomer-shape-to-atresidue-internals (assembler oligomer-shape (residue-shape residue-shape) monomer-context atresidue coordinates &key verbose) +(defmethod apply-monomer-shape-to-atresidue-internals (assembler oligomer-shape (residue-shape residue-shape) monomer-context atresidue &key verbose) (when verbose (let ((*print-pretty* nil)) (format t "apply-monomer-shape-to-atresidue-internals ~% oligomer-shape ~s~% residue-shape ~s~% monomer-context ~s~%" oligomer-shape residue-shape monomer-context))) @@ -117,18 +117,20 @@ If ORIGINAL-ROTAMER-SHAPE is defined then it must be a ROTAMER-SHAPE and we copy (defclass oligomer-shape () - ((name :initarg :name :accessor name) + ((name :initarg :name :reader name) (rotamers-state :initarg :rotamers-state :accessor rotamers-state - :type (member :undefined :incomplete-no-rotamers :incomplete-backbone-rotamers :complete-sidechain-and-backbone-rotamers)) - (oligomer :initarg :oligomer :accessor oligomer) - (rotamers-database :initarg :rotamers-database :accessor rotamers-database) + :type (member :undefined :incomplete-no-rotamers :incomplete-backbone-rotamers :complete-sidechain-and-backbone-rotamers) + :documentation "Mutable slot") + (oligomer :initarg :oligomer :reader oligomer) + (rotamers-database :initarg :rotamers-database :reader rotamers-database) (monomer-shape-info-vector :initarg :monomer-shape-info-vector :reader monomer-shape-info-vector) - (the-root-monomer :initarg :the-root-monomer :accessor the-root-monomer) - (in-monomers :initarg :in-monomers :accessor in-monomers) - (out-monomers :initarg :out-monomers :accessor out-monomers) + (the-root-monomer :initarg :the-root-monomer :reader the-root-monomer) + (in-monomers :initarg :in-monomers :reader in-monomers) + (out-monomers :initarg :out-monomers :reader out-monomers) (monomer-shape-vector :initarg :monomer-shape-vector :accessor monomer-shape-vector :documentation "This is mutable") - (monomer-shape-map :initarg :monomer-shape-map :accessor monomer-shape-map) + (monomer-shape-to-index :initarg :monomer-shape-to-index :reader monomer-shape-to-index) + (monomer-shape-map :initarg :monomer-shape-map :reader monomer-shape-map) )) (defmacro do-oligomer-shape ((monomer-shape monomer monomer-context monomer-shape-kind) oligomer-shape &body body) @@ -161,14 +163,16 @@ to the components and evaluate the BODY in that scope." (defmethod copy-oligomer-shape ((oligomer-shape oligomer-shape)) (let ((monomer-shape-vector (make-array (length (monomer-shape-vector oligomer-shape)))) - (monomer-shape-map (make-hash-table))) + (monomer-shape-map (make-hash-table)) + (monomer-shape-to-index (make-hash-table))) (loop for monomer-shape across (monomer-shape-vector oligomer-shape) for index from 0 for monomer-shape-info across (monomer-shape-info-vector oligomer-shape) for monomer = (monomer monomer-shape-info) for new-monomer-shape = (copy-monomer-shape monomer-shape) do (setf (aref monomer-shape-vector index) new-monomer-shape - (gethash monomer monomer-shape-map) new-monomer-shape)) + (gethash monomer monomer-shape-map) new-monomer-shape + (gethash monomer-shape monomer-shape-to-index) index)) (make-instance 'oligomer-shape :name (name oligomer-shape) :rotamers-state (rotamers-state oligomer-shape) @@ -176,13 +180,16 @@ to the components and evaluate the BODY in that scope." :rotamers-database (rotamers-database oligomer-shape) :monomer-shape-info-vector (monomer-shape-info-vector oligomer-shape) :the-root-monomer (the-root-monomer oligomer-shape) + :monomer-shape-to-index monomer-shape-to-index :in-monomers (in-monomers oligomer-shape) :out-monomers (out-monomers oligomer-shape) :monomer-shape-vector monomer-shape-vector :monomer-shape-map monomer-shape-map))) +(defmethod oligomer-space ((oligomer-shape oligomer-shape)) + "Return the oligomer-space for the oligomer of the OLIGOMER-SHAPE" + (oligomer-space (oligomer oligomer-shape))) - (defmethod print-object ((obj oligomer-shape) stream) (print-unreadable-object (obj stream :type t) (format stream "~s" (name obj)))) @@ -316,9 +323,11 @@ Use the CALLBACK-BACKBONE-ROTAMER-INDEXES and CALLBACK-SIDECHAIN-ROTAMER-INDEXES (shape-info (shape-info foldamer)) #+(or)(kind-order (loop for kind-keys in (kind-keys shape-info) collect (kind kind-keys)))) - (multiple-value-bind (monomer-shape-vector monomer-shape-info-vector the-root-monomer in-monomers out-monomers monomer-shape-map) + (multiple-value-bind (monomer-shape-vector monomer-shape-info-vector monomer-shape-to-index + the-root-monomer in-monomers out-monomers monomer-shape-map) (loop with monomer-shape-vector = (make-array (length (monomers oligomer))) with monomer-shape-info-vector = (make-array (length (monomers oligomer))) + with monomer-shape-to-index = (make-hash-table) with in-monomers = (make-hash-table) with out-monomers = (make-hash-table) with the-root-monomer = nil @@ -364,15 +373,17 @@ Use the CALLBACK-BACKBONE-ROTAMER-INDEXES and CALLBACK-SIDECHAIN-ROTAMER-INDEXES do (unless in-monomer (setf the-root-monomer monomer)) do (setf (gethash monomer out-monomers) out-mons) do (setf (aref monomer-shape-vector index) monomer-shape) + do (setf (gethash monomer-shape monomer-shape-to-index) index) do (setf (aref monomer-shape-info-vector index) monomer-shape-info) ;; do (format t "monomer-context ~a~%" monomer-context) - finally (return (values monomer-shape-vector monomer-shape-info-vector the-root-monomer in-monomers out-monomers monomer-shape-map))) + finally (return (values monomer-shape-vector monomer-shape-info-vector monomer-shape-to-index the-root-monomer in-monomers out-monomers monomer-shape-map))) (let* ((os (apply 'make-instance oligomer-shape-class-name :oligomer oligomer :rotamers-database rotamers-database :monomer-shape-info-vector monomer-shape-info-vector :monomer-shape-vector monomer-shape-vector + :monomer-shape-to-index monomer-shape-to-index :monomer-shape-map monomer-shape-map :the-root-monomer the-root-monomer :in-monomers in-monomers @@ -511,6 +522,7 @@ If UNINITIALIZED then leave them unbound." (random-rotamer-index-impl root-monomer-shape oligomer-shape)))) +#+(or) (defun build-shapes (oligomer-shapes assembler &key monomer-order) (let ((coordinates (chem:make-coordinates (topology:energy-function assembler)))) (chem:energy-function/load-coordinates-into-vector (topology:energy-function assembler) coordinates) @@ -584,33 +596,29 @@ If UNINITIALIZED then leave them unbound." (defun lookup-dihedral-cache (oligomer-shape monomer-shape dihedral-name &key ignore-degrees) (lookup-dihedral-cache-impl oligomer-shape monomer-shape dihedral-name :ignore-degrees ignore-degrees)) - +#+(or) (defun build-externals-from-internals (assembler &key oligomer-shape into-coords) (let ((coords (make-coordinates-for-assembler assembler))) (build-atom-tree-external-coordinates-and-adjust assembler coords oligomer-shape) (copy-joint-positions-into-atoms assembler coords oligomer-shape))) +#+(or) (defun build-oligomer-shape-in-aggregate (assembler oligomer-shape) (let ((coords (make-coordinates-for-assembler assembler))) (update-internals assembler oligomer-shape) (build-atom-tree-external-coordinates-and-adjust assembler coords oligomer-shape) (copy-joint-positions-into-atoms assembler coords oligomer-shape))) -(defun update-internals-for-atresidue (assembler rotamer rotamer-index atresidue) - "This is for writing/adjusting the internals for a single atresidue" - (apply-fragment-internals-to-atresidue rotamer rotamer-index atresidue) - (let ((adjustments (gethash atresidue (internal-adjustments (adjustments assembler))))) - (loop for adjustment in adjustments - do (internal-adjust adjustment assembler)))) - - +#+(or) (defgeneric build-initial-oligomer-shape-externals (assembler oligomer-shape &key orientation coords) (:documentation "Build the oligomer-shape in the coordinates by generating the initial externals")) +#+(or) (defmethod build-initial-oligomer-shape-externals (assembler oligomer-shape &key (coords (make-coordinates-for-assembler assembler)) orientation) (build-atom-tree-external-coordinates-and-adjust assembler coords oligomer-shape orientation) coords) +#+(or) (defun build-all-externals (assembler &key (coords (make-coordinates-for-assembler assembler))) (loop for oligomer-shape in (oligomer-shapes assembler) do (build-atom-tree-external-coordinates-and-adjust assembler coords oligomer-shape oligomer-shape)) @@ -619,16 +627,6 @@ If UNINITIALIZED then leave them unbound." (defun build-oligomer-shape-keep-internals-in-coordinates (assembler oligomer-shape coords) (build-atom-tree-external-coordinates-and-adjust assembler coords oligomer-shape)) -(defun build-all-oligomer-shapes-in-coordinates (assembler &optional (coords (make-coordinates-for-assembler assembler))) - (loop for oligomer-shape in (topology:oligomer-shapes assembler) - do (update-internals assembler oligomer-shape) - do (build-atom-tree-external-coordinates-and-adjust assembler coords oligomer-shape)) - coords) - -(defun build-all-oligomer-shapes-from-internals-in-coordinates (assembler coords) - (loop for oligomer-shape in (topology:oligomer-shapes assembler) - do (build-atom-tree-external-coordinates-and-adjust assembler coords oligomer-shape oligomer-shape))) - (defmethod aggregate (oligomer-shape) "Generate a random oligomer-shape and return the aggregate" (let* ((bs (make-permissible-backbone-rotamers oligomer-shape))) @@ -639,11 +637,14 @@ If UNINITIALIZED then leave them unbound." (coords (make-coordinates-for-assembler ass)) ) (update-internals ass oligomer-shape) - (build-all-oligomer-shapes-from-internals-in-coordinates ass coords) + (update-externals ass :oligomer-shape oligomer-shape + :orientation oligomer-shape + :coords coords) (copy-all-joint-positions-into-atoms ass coords) (aggregate ass))))) +#+(or) (defun assembler-aggregate (assembler oligomer-shapes) (let ((coords (make-coordinates-for-assembler assembler))) (loop for oligomer-shape in oligomer-shapes @@ -667,6 +668,7 @@ If UNINITIALIZED then leave them unbound." (let* ((oligomer-shape (make-oligomer-shape-of-rotamer-shapes oligomer))) (aggregate oligomer-shape))) +#+(or) (defun analyze-oligomer-shape (oligomer-shape) "Print an analysis of the oligomer-shape. Print the number of backbone and sidechain monomer-shapes that have defined rotamer-index's " (let ((backbone-count 0) @@ -725,25 +727,14 @@ If UNINITIALIZED then leave them unbound." (coords (make-coordinates-for-assembler ass)) ) (update-internals ass oligomer-shape) - (build-all-oligomer-shapes-from-internals-in-coordinates ass coords) + (update-externals ass :oligomer-shape oligomer-shape + :orientation oligomer-shape + :coords coords) (copy-all-joint-positions-into-atoms ass coords) (aggregate ass))))) -(defun update-internals (assembler oligomer-shape &key verbose) - "Update the internal coordinates of the assembler for the oligomer-shape -or whatever type the second argument is. -The MONOMER-SHAPEs of the OLIGOMER-SHAPE are used to update the internal coordinates in -the atom-tree. - -ASSEMBLER - the assembler to update. -OLIGOMER-SHAPE - An oligomer-shape (or permissible-rotamers - I think this is wrong) in the assembler." - (fill-internals-from-oligomer-shape assembler oligomer-shape verbose) - (topology:with-orientation (topology:lookup-orientation assembler oligomer-shape) - (adjust-internals assembler oligomer-shape))) - - (defun build-rotamer-shape-map (oligomer) "Build a monomer-to-monomer-shape-map of only rotamer-shapes. The ROTAMER-INDEX slots of the ROTAMER-SHAPEs will be unbound on return." diff --git a/src/lisp/topology/topology-classes.lisp b/src/lisp/topology/topology-classes.lisp index 533ad26e..2b758146 100644 --- a/src/lisp/topology/topology-classes.lisp +++ b/src/lisp/topology/topology-classes.lisp @@ -453,7 +453,6 @@ that is not avoid-out-coupling-plug-name. Otherwise signal an error" (lambda (obj stream) (print-unreadable-object (obj stream :type t)))) - (defun topologys-in-oligomer-space (oligomer-space) (loop for monomer across (monomers oligomer-space) for monomer-names = (monomers monomer) @@ -510,6 +509,8 @@ Examples: (topology:verify-oligomer-space oligomer-space) oligomer-space)) +(defun oligomer-space-contains-monomer (oligomer-space monomer) + (position monomer (monomers oligomer-space))) (defun monomer-position (monomer oligomer-space) "Return the position of the MONOMER in the OLIGOMER-SPACE or NIL"