Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2026 Pengchao Zhang
       3             : 
       4             :    See http://www.plumed.org for more information.
       5             : 
       6             :    This file is part of plumed, version 2.
       7             : 
       8             :    plumed is free software: you can redistribute it and/or modify
       9             :    it under the terms of the GNU Lesser General Public License as published by
      10             :    the Free Software Foundation, either version 3 of the License, or
      11             :    (at your option) any later version.
      12             : 
      13             :    plumed is distributed in the hope that it will be useful,
      14             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      15             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      16             :    GNU Lesser General Public License for more details.
      17             : 
      18             :    You should have received a copy of the GNU Lesser General Public License
      19             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      20             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      21             : #include "function/Function.h"
      22             : #include "core/ActionRegister.h"
      23             : #include <algorithm>
      24             : #include <cmath>
      25             : 
      26             : namespace PLMD {
      27             : namespace pathbias {
      28             : 
      29             : //+PLUMEDOC FUNCTION PROBABILITY_MIX
      30             : /*
      31             : Mix two frozen path potentials with a fixed global normalization.
      32             : 
      33             : ARG contains two nonperiodic energy scalars V_c,V_A. The output is
      34             : -kBT log[(1-eta) exp(-V_c/kBT)+eta exp(-V_A/kBT)/C], where COUPLING is
      35             : eta in [0,1) and LOG_NORMALIZER is log C. With C=Z_A/Z_c this samples the
      36             : normalized mixture (1-eta)q_c+eta q_A. An approximate frozen C still defines
      37             : a conservative potential but changes the actual mixture fraction. C is
      38             : global, not the coordinate-dependent normalizer of CONDITIONAL_PATH.
      39             : Derivatives are (1-chi,chi). All replicas must use identical parameters and
      40             : the same shared full-path inputs. This function neither learns nor biases
      41             : directly; apply the result using BIASVALUE.
      42             : 
      43             : ```plumed
      44             : c: DISTANCE ATOMS=1,2 NOPBC
      45             : vc: CUSTOM ARG=c FUNC=x*x PERIODIC=NO
      46             : va: CUSTOM ARG=c FUNC=2*x*x PERIODIC=NO
      47             : total: PROBABILITY_MIX ARG=vc,va KBT=2.5 COUPLING=0.5 LOG_NORMALIZER=0
      48             : bias: BIASVALUE ARG=total
      49             : ```
      50             : 
      51             : The example illustrates composition only; its C is not a fitted partition
      52             : function ratio. KBT and both arguments must use the same energy units.
      53             : At eta=0 this function returns V_c with derivatives (1,0), but PLUMED still
      54             : evaluates upstream actions: generate a centroid-only graph to avoid evaluating
      55             : an invalid or expensive inactive field. Pure V_A is a separate graph.
      56             : */
      57             : //+ENDPLUMEDOC
      58             : 
      59             : class ProbabilityMix : public function::Function {
      60             :   double kbt=0;
      61             :   double coupling=0;
      62             :   double logNormalizer=0;
      63             : public:
      64             :   explicit ProbabilityMix(const ActionOptions&);
      65             :   static void registerKeywords(Keywords&);
      66             :   void calculate() override;
      67             : };
      68             : 
      69             : PLUMED_REGISTER_ACTION(ProbabilityMix,"PROBABILITY_MIX")
      70             : 
      71          76 : void ProbabilityMix::registerKeywords(Keywords& keys) {
      72          76 :   Function::registerKeywords(keys);
      73          76 :   keys.add("compulsory","KBT","positive thermal energy in the argument energy units");
      74          76 :   keys.add("compulsory","COUPLING","mixture fraction in [0,1)");
      75          76 :   keys.add("compulsory","LOG_NORMALIZER","fixed dimensionless global log partition-function ratio");
      76         152 :   keys.setValueDescription("scalar","the complete-path mixed bias energy");
      77          76 : }
      78             : 
      79          72 : ProbabilityMix::ProbabilityMix(const ActionOptions& ao):Action(ao),Function(ao) {
      80          72 :   parse("KBT",kbt);
      81          72 :   parse("COUPLING",coupling);
      82         144 :   parse("LOG_NORMALIZER",logNormalizer);
      83          72 :   if(getNumberOfArguments()!=2) {
      84           0 :     error("ARG requires the centroid and arithmetic path energies");
      85             :   }
      86         216 :   for(unsigned i=0; i<2; ++i) {
      87         144 :     if(getPntrToArgument(i)->getRank()!=0 || getPntrToArgument(i)->isPeriodic()) {
      88           0 :       error("arguments must be nonperiodic scalars");
      89             :     }
      90             :   }
      91          72 :   if(!std::isfinite(kbt) || kbt<=0 || !std::isfinite(coupling) || coupling<0 || coupling>=1 || !std::isfinite(logNormalizer)) {
      92           0 :     error("KBT must be positive, COUPLING in [0,1), and all parameters finite");
      93             :   }
      94          72 :   checkRead();
      95          72 :   addValueWithDerivatives();
      96          72 :   setNotPeriodic();
      97          72 : }
      98             : 
      99         612 : void ProbabilityMix::calculate() {
     100         612 :   const double vc=getArgument(0);
     101         612 :   if(!std::isfinite(vc)) {
     102           0 :     error("centroid energy must be finite");
     103             :   }
     104         612 :   double value=vc,chi=0;
     105         612 :   if(coupling>0) {
     106         408 :     const double va=getArgument(1);
     107         408 :     const double lhs=std::log1p(-coupling)-vc/kbt;
     108         408 :     const double rhs=std::log(coupling)-va/kbt-logNormalizer;
     109         408 :     if(!std::isfinite(va) || !std::isfinite(lhs) || !std::isfinite(rhs)) {
     110           0 :       error("mixture component log densities must be finite");
     111             :     }
     112         408 :     const double maximum=std::max(lhs,rhs);
     113         408 :     const double sum=std::exp(lhs-maximum)+std::exp(rhs-maximum);
     114         408 :     value=-kbt*(maximum+std::log(sum));
     115         408 :     chi=std::exp(rhs-maximum)/sum;
     116         408 :     if(!std::isfinite(value)) {
     117           0 :       error("mixed energy exceeds floating-point range");
     118             :     }
     119             :   }
     120         612 :   setValue(value);
     121         612 :   setDerivative(0,1-chi);
     122             :   setDerivative(1,chi);
     123         612 : }
     124             : 
     125             : }
     126             : }