Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2026 Rangsiman Ketkaew
3 :
4 : The xsprint module is free software: you can redistribute it and/or modify
5 : it under the terms of the GNU Lesser General Public License as published by
6 : the Free Software Foundation, either version 3 of the License, or
7 : (at your option) any later version.
8 :
9 : The xsprint module is distributed in the hope that it will be useful,
10 : but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 : GNU Lesser General Public License for more details.
13 :
14 : You should have received a copy of the GNU Lesser General Public License
15 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
16 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
17 : #include "core/ActionShortcut.h"
18 : #include "core/ActionRegister.h"
19 : #include "core/ActionWithValue.h"
20 : #include "core/PlumedMain.h"
21 : #include "core/ActionSet.h"
22 :
23 : //+PLUMEDOC MATRIXF XSPRINT
24 : /*
25 : Calculate eXtended SPRINT (xSPRINT) topological variables from adjacency matrices over multiple shells.
26 :
27 : The xSPRINT coordinates, introduced in \cite ketkaew2022, extend the standard [SPRINT](SPRINT.md) coordinates
28 : by incorporating structural information from multiple concentric spherical shells around each atom.
29 : For a system divided into $N$ shells, the xSPRINT coordinate of atom $i$ is defined as:
30 :
31 : $$
32 : S_i^x = \begin{cases}
33 : S_i, & \text{if } r_{ij} < r_0 \\
34 : \frac{1}{\omega} \sum_{n=1}^{N} r_n\, S_i^{(n)}, & \text{otherwise}
35 : \end{cases}
36 : $$
37 :
38 : where $S_i^{(n)}$ is the SPRINT coordinate computed from the adjacency matrix of the $n$-th shell,
39 : $r_n = n\,r_0$ is the shell radius, $r_0$ is the user-defined cutoff for the first shell, and
40 : $\omega = \left(\sum_{n=1}^{N} r_n\right)^2$ is a normalization weight. Adjacency matrix elements
41 : from inner shells are zeroed out before computing each shell's SPRINT to avoid redundancy.
42 :
43 : When `NSHELLS=1` the xSPRINT coordinate reduces to the standard SPRINT coordinate.
44 :
45 : The following example computes the xSPRINT coordinates for a 7-atom cluster using 3 shells
46 : with a shell radius of 0.1:
47 :
48 : ```plumed
49 : xs: XSPRINT GROUP=1-7 SWITCH={RATIONAL R_0=0.1} R0=0.1 NSHELLS=3
50 : PRINT ARG=xs.* FILE=colvar
51 : ```
52 :
53 : Note that `R0` defines the xSPRINT shell radius ($r_0$ in the equation) and is independent
54 : of `R_0` inside `SWITCH`, which controls the shape of the switching function. For shell $n$,
55 : the switching function `R_0` is automatically set to $n \times$ `R0`.
56 :
57 : This example computes xSPRINT with 2 shells for two groups of atoms:
58 :
59 : ```plumed
60 : xs: XSPRINT ...
61 : GROUP1=1-7 GROUP2=8-14
62 : SWITCH11={RATIONAL R_0=2.6 NN=6 MM=12}
63 : SWITCH12={RATIONAL R_0=2.2 NN=6 MM=12}
64 : SWITCH22={RATIONAL R_0=2.2 NN=6 MM=12}
65 : R0=2.0
66 : NSHELLS=2
67 : ...
68 :
69 : PRINT ARG=xs.* FILE=colvar
70 : ```
71 :
72 : You can also provide a pre-built matrix (for a single-shell case this is equivalent to SPRINT):
73 :
74 : ```plumed
75 : b: BRIDGE_MATRIX GROUP=1-7 BRIDGING_ATOMS=8-100 SWITCH={RATIONAL R_0=0.2}
76 : xs: XSPRINT MATRIX=b NSHELLS=1
77 : PRINT ARG=xs.* FILE=colvar
78 : ```
79 :
80 : */
81 : //+ENDPLUMEDOC
82 :
83 : namespace PLMD {
84 : namespace xsprint {
85 :
86 : class Xsprint : public ActionShortcut {
87 : public:
88 : static void registerKeywords(Keywords& keys);
89 : explicit Xsprint(const ActionOptions&);
90 : };
91 :
92 : // Function to replace all R_0=<value> entries in a string with R_0=<newval>
93 3 : static std::string replaceR0(const std::string& input, double newval) {
94 3 : std::string result = input;
95 3 : std::string needle = "R_0=";
96 : std::string newvalstr;
97 3 : Tools::convert(newval, newvalstr);
98 : size_t pos = 0;
99 6 : while((pos = result.find(needle, pos)) != std::string::npos) {
100 3 : size_t valstart = pos + needle.length();
101 3 : size_t valend = result.find_first_of(" \t},", valstart);
102 3 : if(valend == std::string::npos) {
103 : valend = result.length();
104 : }
105 3 : result.replace(valstart, valend - valstart, newvalstr);
106 3 : pos = valstart + newvalstr.length();
107 : }
108 3 : return result;
109 : }
110 :
111 : PLUMED_REGISTER_ACTION(Xsprint,"XSPRINT")
112 :
113 8 : void Xsprint::registerKeywords(Keywords& keys) {
114 8 : ActionShortcut::registerKeywords( keys );
115 8 : keys.add("optional","MATRIX","the matrix that you would like to perform xSPRINT on (used only when NSHELLS=1)");
116 8 : keys.add("numbered","GROUP","specifies the list of atoms that should be assumed indistinguishable");
117 8 : keys.add("numbered","SWITCH","specify the switching function to use between two sets of indistinguishable atoms");
118 8 : keys.add("compulsory","NSHELLS","1","the number of concentric shells to use for the xSPRINT calculation");
119 8 : keys.add("optional","R0","the radius of the first shell for xSPRINT (shell n has radius n*R0). Required when NSHELLS>1. This is independent of the R_0 in SWITCH which controls the switching function shape.");
120 8 : keys.needsAction("CONTACT_MATRIX");
121 8 : keys.needsAction("SPRINT");
122 8 : keys.needsAction("CUSTOM");
123 8 : keys.needsAction("SELECT_COMPONENTS");
124 8 : keys.needsAction("SORT");
125 8 : keys.needsAction("COMBINE");
126 16 : keys.addOutputComponent("coord","default","scalar","the xsprint coordinates");
127 8 : keys.addDOI("10.1021/acs.jpclett.1c04004");
128 8 : }
129 :
130 4 : Xsprint::Xsprint(const ActionOptions& ao):
131 : Action(ao),
132 4 : ActionShortcut(ao) {
133 4 : const std::string label = getShortcutLabel();
134 :
135 4 : unsigned nshells = 1;
136 4 : parse("NSHELLS", nshells);
137 4 : if( nshells<1 ) {
138 2 : error("NSHELLS should be at least 1");
139 : }
140 : std::string matinp;
141 6 : parse("MATRIX", matinp);
142 3 : if( !matinp.empty() && nshells > 1 ) {
143 0 : error("MATRIX keyword can only be used with NSHELLS=1");
144 : }
145 3 : double r0 = 0.0;
146 3 : parse("R0", r0);
147 3 : if( nshells > 1 && r0 <= 0.0 ) {
148 0 : plumed_merror("R0 keyword is required and must be positive when NSHELLS > 1");
149 : }
150 :
151 : // NSHELLS == 1 : use existing SPRINT
152 3 : if( nshells == 1 ) {
153 2 : if( !matinp.empty() ) {
154 2 : readInputLine( label + "_sprint: SPRINT MATRIX=" + matinp );
155 : } else {
156 2 : readInputLine( label + "_sprint: SPRINT " + convertInputLineToString() );
157 : }
158 : // Forward SPRINT's coord-N outputs as our own coord-N outputs
159 2 : for(unsigned i = 0;; ++i) {
160 : std::string inum;
161 30 : Tools::convert( i, inum );
162 60 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>(
163 30 : label + "_sprint_coord-" + inum );
164 30 : if( !av ) {
165 : break;
166 : }
167 84 : readInputLine( label + "_coord-" + inum + ": COMBINE ARG=" +
168 56 : label + "_sprint_coord-" + inum + " PERIODIC=NO" );
169 28 : }
170 : return;
171 : }
172 :
173 : // NSHELLS > 1 : compute xSPRINT using SPRINT per shell
174 1 : std::string baseargs = convertInputLineToString();
175 :
176 : // Shell 1: create CONTACT_MATRIX with R_0=r0 and compute SPRINT
177 1 : std::string shell1args = replaceR0(baseargs, r0);
178 1 : readInputLine( label + "_shell1_mat: CONTACT_MATRIX " + shell1args );
179 2 : readInputLine( label + "_shell1_sprint: SPRINT MATRIX=" + label + "_shell1_mat" );
180 :
181 : // Determine group sizes from shell-1 matrix (this is needed for final sorting)
182 1 : std::string mat1 = label + "_shell1_mat";
183 1 : std::vector<unsigned> nin_group;
184 1 : for(unsigned i = 1;; ++i) {
185 : std::string inum;
186 1 : Tools::convert( i, inum );
187 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( mat1 + inum + inum );
188 1 : if( !av ) {
189 : break;
190 : }
191 0 : nin_group.push_back( (av->copyOutput(0))->getShape()[0] );
192 0 : }
193 1 : if( nin_group.empty() ) {
194 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( mat1 );
195 1 : nin_group.push_back( (av->copyOutput(0))->getShape()[0] );
196 : }
197 :
198 : // For shells 2..N: create CONTACT_MATRIX with R_0=s*r0,
199 : // subtract inner shell, compute SPRINT, respectively
200 3 : for(unsigned s = 2; s <= nshells; ++s) {
201 : std::string snum, sprev;
202 2 : Tools::convert( s, snum );
203 2 : Tools::convert( s - 1, sprev );
204 :
205 2 : double shell_r0 = static_cast<double>(s) * r0;
206 2 : std::string shellargs = replaceR0(baseargs, shell_r0);
207 4 : readInputLine( label + "_shell" + snum + "_mat: CONTACT_MATRIX " + shellargs );
208 :
209 : // Subtract inner-shell matrix to isolate contacts in this shell only
210 6 : readInputLine( label + "_shell" + snum + "_diff: CUSTOM ARG=" +
211 4 : label + "_shell" + snum + "_mat," +
212 2 : label + "_shell" + sprev + "_mat FUNC=x-y PERIODIC=NO" );
213 :
214 6 : readInputLine( label + "_shell" + snum + "_sprint: SPRINT MATRIX=" +
215 4 : label + "_shell" + snum + "_diff" );
216 : }
217 :
218 : // Combine per-shell SPRINT vectors with weights r_n / omega
219 : // where r_n = n*r0, omega = (sum_{n=1}^{N} r_n)^2
220 1 : double sum_rn = r0 * static_cast<double>(nshells) * static_cast<double>(nshells + 1) / 2.0;
221 1 : double omega = sum_rn * sum_rn;
222 :
223 : std::string combine_args, coefficients;
224 4 : for(unsigned s = 1; s <= nshells; ++s) {
225 : std::string snum;
226 3 : Tools::convert( s, snum );
227 3 : if( s > 1 ) {
228 : combine_args += ",";
229 : coefficients += ",";
230 : }
231 6 : combine_args += label + "_shell" + snum + "_sprint_sp";
232 : std::string wstr;
233 3 : Tools::convert( static_cast<double>(s) * r0 / omega, wstr );
234 : coefficients += wstr;
235 : }
236 :
237 2 : readInputLine( label + "_xsp: COMBINE ARG=" + combine_args +
238 1 : " COEFFICIENTS=" + coefficients + " PERIODIC=NO" );
239 :
240 : // Sort xSPRINT coordinates per group
241 1 : unsigned k = 0, kk = 0;
242 2 : for(unsigned j = 0; j < nin_group.size(); ++j) {
243 : std::string jnum, knum;
244 1 : Tools::convert( j+1, jnum );
245 1 : Tools::convert( k+1, knum );
246 : k++;
247 3 : std::string sort_act = label + "_selection" + jnum + ": SELECT_COMPONENTS ARG=" + label + "_xsp COMPONENTS=" + knum;
248 14 : for(unsigned n = 1; n < nin_group[j]; ++n) {
249 13 : Tools::convert( k+1, knum );
250 13 : sort_act += "," + knum;
251 : k++;
252 : }
253 1 : readInputLine( sort_act );
254 3 : readInputLine( label + jnum + ": SORT ARG=" + label + "_selection" + jnum );
255 15 : for(unsigned n = 0; n < nin_group[j]; ++n) {
256 : std::string nnum;
257 14 : Tools::convert( kk, knum );
258 14 : Tools::convert( n+1, nnum );
259 14 : kk++;
260 70 : readInputLine( label + "_coord-" + knum + ": COMBINE ARG=" + label + jnum + "." + nnum + " PERIODIC=NO" );
261 : }
262 : }
263 16 : }
264 :
265 : }
266 : }
|
