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 CONDITIONAL_PATH
30 : /*
31 : Evaluate a fixed conditional path log mixture, with its full chain rule.
32 :
33 : For positive path score a and positive conditional normalizer m(c), this
34 : function returns log R = log[(1-lambda)+lambda*a/m(c)]. ARG must contain
35 : exactly three nonperiodic scalars: c, log(a), log(m(c)), in that order.
36 : COUPLING is lambda in [0,1). LOWER and UPPER bound the admitted domain of c;
37 : evaluation outside this closed interval fails instead of extrapolating.
38 :
39 : The derivatives with respect to the three arguments are (0,chi,-chi),
40 : where chi=lambda*a/(m*R). The explicit c argument only checks the domain.
41 : The normalizer argument MUST retain its dependence on c: freezing a fitted
42 : function does not mean dropping its derivative. All functions must be fixed
43 : during production. This action does not learn m or preserve external file
44 : identity across restarts. Verify frozen input/model hashes before restarting.
45 :
46 : A bias V=B_c-kBT*log R is built using CUSTOM and BIASVALUE. kBT must be in
47 : the same energy units as B_c. With an exact m=E_0[a|c], the selected c marginal
48 : of the fixed B_c ensemble is preserved. Approximate m still defines a fixed
49 : potential with full-path weight exp(beta*V), but does not imply that marginal
50 : identity or any efficiency improvement. Every bead in one path shares the
51 : same path weight; beads are not independent samples.
52 :
53 : The following two-replica example uses a mean Cartesian distance component,
54 : then squares it, rather than averaging the squared components. It assumes
55 : consistent nonperiodic coordinate images across the replicas. This scalar
56 : example does not implement arbitrary molecular Cartesian-centroid CVs.
57 :
58 : ```plumed
59 : #SETTINGS NREPLICAS=2
60 : d: DISTANCE ATOMS=1,2 COMPONENTS NOPBC
61 : mean: ENSEMBLE ARG=d.x
62 : c: CUSTOM ARG=mean.d.x FUNC=x*x PERIODIC=NO
63 : h: CUSTOM ARG=d.x FUNC=exp(-x*x/2) PERIODIC=NO
64 : fraction: ENSEMBLE ARG=h
65 : loga: CUSTOM ARG=fraction.h FUNC=log(0.1+x) PERIODIC=NO
66 : logm: CUSTOM ARG=c FUNC=-1+0.2*x PERIODIC=NO
67 : ratio: CONDITIONAL_PATH ARG=c,loga,logm COUPLING=0.5 LOWER=0 UPPER=4
68 : correction: CUSTOM ARG=ratio FUNC=-2.5*x PERIODIC=NO
69 : bias: BIASVALUE ARG=correction
70 : PRINT ARG=c,loga,logm,ratio,bias.* FILE=COLVAR
71 : ```
72 :
73 : The illustrative logm above is not a fitted conditional reference. Supply a
74 : qualified frozen normalizer and a separate fixed centroid bias for production.
75 : The molecular dynamics adapter must provide consistent dynamical force
76 : scaling and single-path energy accounting. This action does not apply
77 : additional factors of the bead count.
78 : */
79 : //+ENDPLUMEDOC
80 :
81 : class ConditionalPath : public function::Function {
82 : double coupling;
83 : double lower;
84 : double upper;
85 : public:
86 : explicit ConditionalPath(const ActionOptions&);
87 : static void registerKeywords(Keywords& keys);
88 : void calculate() override;
89 : };
90 :
91 : PLUMED_REGISTER_ACTION(ConditionalPath,"CONDITIONAL_PATH")
92 :
93 92 : void ConditionalPath::registerKeywords(Keywords& keys) {
94 92 : Function::registerKeywords(keys);
95 92 : keys.add("compulsory","COUPLING","mixture fraction in the half-open interval [0,1)");
96 92 : keys.add("compulsory","LOWER","inclusive lower limit of the conditioning argument");
97 92 : keys.add("compulsory","UPPER","inclusive upper limit of the conditioning argument");
98 184 : keys.setValueDescription("scalar","the dimensionless conditional path log mixture");
99 92 : }
100 :
101 88 : ConditionalPath::ConditionalPath(const ActionOptions& ao):
102 : Action(ao),
103 : Function(ao),
104 88 : coupling(0),
105 88 : lower(0),
106 88 : upper(0) {
107 88 : parse("COUPLING",coupling);
108 88 : parse("LOWER",lower);
109 176 : parse("UPPER",upper);
110 : if(getNumberOfArguments()!=3) {
111 0 : error("ARG must contain conditioning coordinate, log score, and log normalizer");
112 : }
113 352 : for(unsigned i=0; i<3; ++i) {
114 264 : if(getPntrToArgument(i)->getRank()!=0 || getPntrToArgument(i)->isPeriodic()) {
115 0 : error("arguments must be nonperiodic scalars");
116 : }
117 : }
118 88 : if(!std::isfinite(coupling) || coupling<0 || coupling>=1) {
119 16 : error("COUPLING must be finite and in [0,1)");
120 : }
121 72 : if(!std::isfinite(lower) || !std::isfinite(upper) || lower>=upper) {
122 0 : error("LOWER and UPPER must be finite and strictly increasing");
123 : }
124 72 : checkRead();
125 72 : addValueWithDerivatives();
126 72 : setNotPeriodic();
127 72 : log.printf(" fixed conditional path mixture with coupling %g in [%g,%g]\n",
128 : coupling,lower,upper);
129 104 : }
130 :
131 240 : void ConditionalPath::calculate() {
132 240 : const double c=getArgument(0), loga=getArgument(1), logm=getArgument(2);
133 240 : if(!std::isfinite(c) || !std::isfinite(loga) || !std::isfinite(logm)) {
134 0 : error("conditional path arguments must be finite");
135 : }
136 240 : if(c<lower || c>upper) {
137 8 : error("conditioning coordinate is outside the frozen model domain");
138 : }
139 232 : double logratio=0, chi=0;
140 232 : if(coupling>0) {
141 160 : const double lhs=std::log1p(-coupling);
142 160 : const double rhs=std::log(coupling)+loga-logm;
143 160 : if(!std::isfinite(rhs)) {
144 0 : error("conditional log ratio exceeds floating-point range");
145 : }
146 160 : logratio=std::max(lhs,rhs)+std::log1p(std::exp(-std::abs(lhs-rhs)));
147 160 : chi=std::exp(rhs-logratio);
148 : }
149 232 : setValue(logratio);
150 : setDerivative(0,0);
151 : setDerivative(1,chi);
152 : setDerivative(2,-chi);
153 232 : }
154 :
155 : }
156 : }
|
