Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2020 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "adjmat/AdjacencyMatrixBase.h"
23 : #include "multicolvar/MultiColvarShortcuts.h"
24 : #include "core/ActionRegister.h"
25 : #include "core/ActionShortcut.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/ActionSet.h"
28 : #include "tools/IFile.h"
29 :
30 : //+PLUMEDOC MATRIX HBPAMM_MATRIX
31 : /*
32 : Adjacency matrix in which two electronegative atoms are adjacent if they are hydrogen bonded
33 :
34 : This method allows you to calculate an adjacency matrix and use all the methods discussed in the documentation for
35 : the [CONTACT_MATRIX](CONTACT_MATRIX.md) upon it to define CVs. The $i,j$ element of the matrix that is calculated
36 : by this action is one if there is a hydrogen bond connecting atom $i$ to atom $j$. Furthermore, we determine whether
37 : there is a hydrogen atom between these two atoms by using the PAMM technique that is discussed in the articles from the
38 : bibliography below and in the documentation for the [PAMM](PAMM.md) action.
39 :
40 : The example shown below illustrates how the method is used in practice
41 :
42 : ```plumed
43 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
44 : m: HBPAMM_MATRIX ...
45 : GROUP=1-192:3 GROUPC=2-192:3,3-192:3
46 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
47 : ...
48 : PRINT ARG=m FILE=colvar
49 : ```
50 :
51 : The input above is outputting the full hbpamm matrix. However, this action is perhaps more usefully used to investigate the connectivity between
52 : a collection of water molecules in liquid water. Importantly, however,
53 : the output matrix here is __not__ symmetric. We thus calculate the number of hydrogen bonds that these atoms are donating using the
54 : following input:
55 :
56 : ```plumed
57 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
58 : m: HBPAMM_MATRIX ...
59 : GROUP=1-192:3 GROUPC=2-192:3,3-192:3
60 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
61 : ...
62 : ones: ONES SIZE=64
63 : rsums: MATRIX_VECTOR_PRODUCT ARG=m,ones
64 : DUMPATOMS ATOMS=1-192:3 ARG=rsums FILE=donors.xyz
65 : ```
66 :
67 : To calculate the number of hydrogen bonds these atoms accept we would use the following input:
68 :
69 : ```plumed
70 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
71 : m: HBPAMM_MATRIX ...
72 : GROUP=1-192:3 GROUPC=2-192:3,3-192:3
73 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
74 : ...
75 : mT: TRANSPOSE ARG=m
76 : ones: ONES SIZE=64
77 : rsums: MATRIX_VECTOR_PRODUCT ARG=mT,ones
78 : DUMPATOMS ATOMS=1-192:3 ARG=rsums FILE=acceptors.xyz
79 : ```
80 :
81 : To calculate the total number of hydorgen bonds these atoms participate in you would use an input like this:
82 :
83 : ```plumed
84 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
85 : m: HBPAMM_MATRIX ...
86 : GROUP=1-192:3 GROUPC=2-192:3,3-192:3
87 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
88 : ...
89 : mT: TRANSPOSE ARG=m
90 : hbmat: CUSTOM ARG=m,mT FUNC=x+y PERIODIC=NO
91 : ones: ONES SIZE=64
92 : rsums: MATRIX_VECTOR_PRODUCT ARG=hbmat,ones
93 : DUMPATOMS ATOMS=1-192:3 ARG=rsums FILE=hbonds.xyz
94 : ```
95 :
96 : If you want to investigate whether there are hydrogen bonds between two groups of molecules you can use an input like this:
97 :
98 : ```plumed
99 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
100 : m: HBPAMM_MATRIX ...
101 : GROUPA=1 GROUPB=2-192:3 GROUPC=2-192:3,3-192:3
102 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
103 : ...
104 : PRINT ARG=m FILE=colvar
105 : ```
106 :
107 : This input outputs a $1\times 63$ matrix in which the $1,i$th element tells you whether or not atom 1 donates a hydrogen bond
108 : to the $i$th element in the group of 63 atoms that was specified using the ACCEPTORS keyword. The $i,1$th element of the
109 : transpose of this matrix tells you if the $i$th center donates a hydrogen bond to atom 1.
110 :
111 : In general, it is better to use this action through the [HBPAMM_SA](HBPAMM_SA.md), [HBPAMM_SD](HBPAMM_SD.md) and [HBPAMM_SH](HBPAMM_SH.md)
112 : keywords, which can be used to calculate the number of hydrogen bonds each donor, acceptor or hydrogen atom in your system participates in.
113 :
114 : ## Periodic boundary conditions
115 :
116 : Notice that in all the inputs above the distances values that enter the pamm expressions are calculated in a way that takes the
117 : periodic boundary conditions into account. If you want to ignore the periodic boundary conditions you can use the NOPBC flag as shown below.
118 :
119 : ```plumed
120 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
121 : m: HBPAMM_MATRIX ...
122 : GROUP=1-192:3 GROUPC=2-192:3,3-192:3
123 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
124 : NOPBC
125 : ...
126 : ```
127 :
128 : ## COMPONENTS flag
129 :
130 : If you add the flag COMPONENTS to the input as shown below:
131 :
132 : ```plumed
133 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
134 : c4: HBPAMM_MATRIX ...
135 : GROUP=1-192:3 GROUPC=2-192:3,3-192:3
136 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
137 : COMPONENTS
138 : ...
139 : ```
140 :
141 : then four matrices with the labels `c4.w`, `c4.x`, `c4.y` and `c4.z` are output by the action. The matrix with the label `c4.w` is the adjacency matrix
142 : that would be output if you had not added the COMPONENTS flag. The $i,j$ component of the matrices `c4.x`, `c4.y` and `c4.z` contain the $x$, $y$ and $z$
143 : components of the vector connecting atoms $j$ and $k$. Importantly, however, the components of these vectors are only stored in `c4.x`, `c4.y` and `c4.z`
144 : if the elements of `c4.w` are non-zero. Using the COMPONENTS flag in this way ensures that you can use HBPAMM_MATRIX in tandem with many of the functionalities
145 : that are part of the [symfunc module](module_symfunc.md). Remember, however, that the $i,j$ element of the HBPAMM_MATRIX is only non-zero if atom $i$ donates
146 : a hydrogen bond to atom $j$. __You cannot use HBPAMM_MATRIX to identify the set of atoms that each atom is hydrogen bonded to.__
147 :
148 : ## The MASK keyword
149 :
150 : You use the MASK keyword with HBPAMM_MATRIX in the same way that is used in [CONTACT_MATRIX](CONTACT_MATRIX.md). This keyword thus expects a vector in input,
151 : which tells HBOND_MATRIX that it is safe to not calculate certain rows of the output matrix. An example where this keyword is used is shown below:
152 :
153 : ```plumed
154 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
155 : # The atoms that are of interest
156 : ow: GROUP ATOMS=1-1650
157 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
158 : center: FIXEDATOM AT=2.5,2.5,2.5
159 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
160 : sphere: INSPHERE ATOMS=ow CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
161 : # Calculates cooordination numbers
162 : cmap: HBPAMM_MATRIX ...
163 : GROUP=ow GROUPC=1650-3000
164 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
165 : MASK=sphere
166 : ...
167 : ones: ONES SIZE=1650
168 : cc: MATRIX_VECTOR_PRODUCT ARG=cmap,ones
169 : # Multiply coordination numbers by sphere vector
170 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
171 : # Sum of coordination numbers for atoms that are in the sphere of interest
172 : numer: SUM ARG=prod PERIODIC=NO
173 : # Number of atoms that are in sphere of interest
174 : denom: SUM ARG=sphere PERIODIC=NO
175 : # Average coordination number for atoms in sphere of interest
176 : av: CUSTOM ARG=prod,denom FUNC=x/y PERIODIC=NO
177 : # And print out final CV to a file
178 : PRINT ARG=av FILE=colvar STRIDE=1
179 : ```
180 :
181 : This input calculates the average number of hydrogen bonds each of the atoms that are within a spherical region that is centered on the point
182 : $(2.5,2.5,2.5)$ donate.
183 :
184 : */
185 : //+ENDPLUMEDOC
186 :
187 : //+PLUMEDOC COLVAR HBPAMM_SA
188 : /*
189 : Calculate the number of hydrogen bonds each acceptor participates in using the HBPamm method
190 :
191 : This shortcut action allows you to calculate the number of hydrogen bonds each of the atoms specified using the SITES
192 : keyword accepts from its neighbours. The number of hydrogen bonds that a particular site accepts is determined by using the
193 : PAMM tehcnique that is discussed in the articles from the bibliography below and in the documentation for the [PAMM](PAMM.md) action
194 :
195 : The following example shows how you can use this action to calculate how many hydrogen bonds each of the water moecules in
196 : a box of water accepts from its neighbours.
197 :
198 : ```plumed
199 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
200 : acceptors: HBPAMM_SA ...
201 : SITES=1-192:3 HYDROGENS=2-192:3,3-192:3
202 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
203 : REGULARISE=0.001 GAUSS_CUTOFF=6.25
204 : ...
205 : DUMPATOMS ARG=acceptors ATOMS=1-192:3 FILE=acceptors.xyz
206 : ```
207 :
208 : The output here is an xyz with five columns. As explained in the documentation for [DUMPATOMS](DUMPATOMS.md), the first four
209 : columns are the usual columns that you would expect in an xyz file. The fifth column then contains the number of hydrogen bonds
210 : that have been accepted by each of the atoms.
211 :
212 : In the example input above all the atoms specified using the SITE keyword can both accept and donate hydrogen bonds. If the group
213 : of atoms that can accept hydrogen bonds is different from the group of atoms that can donate hydrogen bonds then you use the ACCEPTORS and
214 : DONORS keyword as shown below:
215 :
216 : ```plumed
217 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
218 : acceptors: HBPAMM_SA ...
219 : ACCEPTORS=1-9:3 DONORS=1-192:3
220 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
221 : HYDROGENS=2-192:3,3-192:3
222 : ...
223 : DUMPATOMS ARG=acceptors ATOMS=1-9:3 FILE=acceptors.xyz
224 : ```
225 :
226 : This input will still output an xyz file with five columns. The fifth column in this file will give the number of hydrogen bonds that each
227 : of the three atoms that were specified using the ACCEPTORS keyword accept from the atoms that were specified by the DONORS keyword. Notice also that the
228 : as the REGULARISE and GAUSS_CUTOFF keywords were not used in this input these quantities are set to the default values that are given in the table below.
229 :
230 : ## The NOPBC flag
231 :
232 : When you use the [HBPAMM_MATRIX](HBPAMM_MATRIX.md) you need to calculate the distances between various pairs of atoms. In all the inputs that have been provided
233 : above we calculate these distances between pairs of atoms in a way that takes the periodic boundary conditions into account. If you want to ignore the periodic
234 : boundary conditions you can use the NOPBC flag as shown below.
235 :
236 : ```plumed
237 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
238 : acceptors: HBPAMM_SA ...
239 : SITES=1-192:3 HYDROGENS=2-192:3,3-192:3
240 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
241 : NOPBC
242 : ...
243 : DUMPATOMS ARG=acceptors ATOMS=1-192:3 FILE=acceptors.xyz
244 : ```
245 :
246 : ## The MASK keyword
247 :
248 : You can use the MASK keyword with this action in the same way that it is used with [COORDINATIONNUMBER](COORDINATIONNUMBER.md). This keyword thus expects a vector in
249 : input, which tells PLUMED the atoms for which you do not need to calculate the number of hydrogen bonds. As illustrated below, this is useful if you are using functionality
250 : from the [volumes module](module_volumes.md) to calculate the average number of hydrogen bonds those atoms that lie in a certain part of the simulation box accept from their neighbors.
251 :
252 : ```plumed
253 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
254 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
255 : center: FIXEDATOM AT=2.5,2.5,2.5
256 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
257 : sphere: INSPHERE ATOMS=1-192:3 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
258 : # Calculate the number of hydrogen bonds that have been accepted
259 : acceptors: HBPAMM_SA ...
260 : SITES=1-192:3 HYDROGENS=2-192:3,3-192:3
261 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
262 : MASK=sphere
263 : ...
264 : # Multiply number of hydrogen bonds by sphere vector
265 : prod: CUSTOM ARG=acceptors,sphere FUNC=x*y PERIODIC=NO
266 : # Total number of hydrogen bonds accepted for atoms that are in the sphere of interest
267 : numer: SUM ARG=prod PERIODIC=NO
268 : # Number of atoms that are in sphere of interest
269 : denom: SUM ARG=sphere PERIODIC=NO
270 : # Average number of accepted hydrogen bonds for atoms in sphere of interest
271 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
272 : # And print out final CV to a file
273 : PRINT ARG=av FILE=colvar STRIDE=1
274 : ```
275 :
276 : This input calculate the number of hydrogen bonds that are accepted by those atoms that are within a spherical region that is centered on the point
277 : $(2.5,2.5,2.5)$.
278 :
279 : ## Optimisation details
280 :
281 : If you expand the inputs above you will see that they all the [HBPAMM_MATRIX](HBPAMM_MATRIX.md) action. The matrix that this action computes is sparse
282 : as each atom only forms hydroen bonds with a relatively small number of neighbors. The vast majority of the elements in the [HBPAMM_MATRIX](HBPAMM_MATRIX.md)
283 : are thsu zero. To reduce the amount of memory that PLUMED requires PLUMED uses sparse matrix storage. Consequently, whenever you calculate and store a contact
284 : matrix only the elements of the matrix that are non-zero are stored.
285 :
286 : We also use the sparsity of the [HBPAMM_MATRIX](HBPAMM_MATRIX.md) to make the time required to compute the matrix scale linearly rather than quadratically with
287 : the number of atoms. We know that element $i,j$ is only non-zero if the donor and acceptor atoms are within:
288 :
289 : $$
290 : r_c = \textrm{max}_i \left( \mu_i^{(ad)} + \sqrt(2 g_{cut}) | \sqrt{\lambda_{i,max}} v_\textrm{i,max}^{(ad)} | \right)
291 : $$
292 :
293 : where $i$ runs over all the PAMM kernels. $g_{cut}$ is the gaussian cutoff that is specified using the REGULARIZE
294 : keyword, $\mu_i^{(ad)}$ is the mean acceptor-donor distance for the $i$th kernel, $\lambda_{i,max}$ is the largest
295 : eigenvalue of the covariance for the $i$th kernel and $v_\textrm{i,max}^{(ad)}$ is the acceptor-donor distance component
296 : of its corresponding eigenvector. We can determine that many pairs of atoms are further apart than $r_c$ without computing the
297 : distance between these atoms by using divide and conquer strategies such as linked lists and neighbour lists. Furthermore, we do not even
298 : need to set a `D_MAX` parameter as we would have to if we were using [CONTACT_MATRIX](CONTACT_MATRIX.md) as the above cutoff can be determined
299 : from the parameters of the Gaussian kernels directly.
300 :
301 : */
302 : //+ENDPLUMEDOC
303 :
304 : //+PLUMEDOC COLVAR HBPAMM_SD
305 : /*
306 : Calculate the number of hydrogen bonds each donor participates in using the HBPamm method
307 :
308 : This shortcut action allows you to calculate the number of hydrogen bonds each of the atoms specified using the SITES
309 : keyword donates to its neighbours. The number of hydrogen bonds that a particular site donates is determined by using the
310 : PAMM tehcnique that is discussed in the articles from the bibliography below and in the documentation for the [PAMM](PAMM.md) action
311 :
312 : The following example shows how you can use this action to calculate how many hydrogen bonds each of the water moecules in
313 : a box of water donates to its neighbours.
314 :
315 : ```plumed
316 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
317 : donors: HBPAMM_SD ...
318 : SITES=1-192:3 HYDROGENS=2-192:3,3-192:3
319 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
320 : REGULARISE=0.001 GAUSS_CUTOFF=6.25
321 : ...
322 : DUMPATOMS ARG=donors ATOMS=1-192:3 FILE=donors.xyz
323 : ```
324 :
325 : The output here is an xyz with five columns. As explained in the documentation for [DUMPATOMS](DUMPATOMS.md), the first four
326 : columns are the usual columns that you would expect in an xyz file. The fifth column then contains the number of hydrogen bonds
327 : that have been donated by each of the atoms.
328 :
329 : In the example input above all the atoms specified using the SITE keyword can both accept and donate hydrogen bonds. If the group
330 : of atoms that can accept hydrogen bonds is different from the group of atoms that can donate hydrogen bonds then you use the ACCEPTORS and
331 : DONORS keyword as shown below:
332 :
333 : ```plumed
334 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
335 : donors: HBPAMM_SD ...
336 : DONORS=1-9:3 ACCEPTORS=1-192:3
337 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
338 : HYDROGENS=2-192:3,3-192:3
339 : ...
340 : DUMPATOMS ARG=donors ATOMS=1-9:3 FILE=donors.xyz
341 : ```
342 :
343 : This input will still output an xyz file with five columns. The fifth column in this file will give the number of hydrogen bonds that each
344 : of the three atoms that were specified using the DONORS keyword donate to atoms that were specified by the ACCEPTORS keyword. Notice also that the
345 : as the REGULARISE and GAUSS_CUTOFF keywords were not used in this input these quantities are set to the default values that are given in the table below.
346 :
347 : ## The NOPBC flag
348 :
349 : When you use the [HBPAMM_MATRIX](HBPAMM_MATRIX.md) you need to calculate the distances between various pairs of atoms. In all the inputs that have been provided
350 : above we calculate these distances between pairs of atoms in a way that takes the periodic boundary conditions into account. If you want to ignore the periodic
351 : boundary conditions you can use the NOPBC flag as shown below.
352 :
353 : ```plumed
354 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
355 : donors: HBPAMM_SD ...
356 : SITES=1-192:3 HYDROGENS=2-192:3,3-192:3
357 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
358 : NOPBC
359 : ...
360 : DUMPATOMS ARG=donors ATOMS=1-192:3 FILE=donors.xyz
361 : ```
362 :
363 : ## The MASK keyword
364 :
365 : You can use the MASK keyword with this action in the same way that it is used with [COORDINATIONNUMBER](COORDINATIONNUMBER.md). This keyword thus expects a vector in
366 : input, which tells PLUMED the atoms for which you do not need to calculate the number of hydrogen bonds. As illustrated below, this is useful if you are using functionality
367 : from the [volumes module](module_volumes.md) to calculate the average number of hydrogen bonds those atoms that lie in a certain part of the simulation box donoate to their neighbors.
368 :
369 : ```plumed
370 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
371 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
372 : center: FIXEDATOM AT=2.5,2.5,2.5
373 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
374 : sphere: INSPHERE ATOMS=1-192:3 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
375 : # Calculate the number of hydrogen bonds that have been donated
376 : donors: HBPAMM_SD ...
377 : SITES=1-192:3 HYDROGENS=2-192:3,3-192:3
378 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
379 : MASK=sphere
380 : ...
381 : # Multiply number of hydrogen bonds by sphere vector
382 : prod: CUSTOM ARG=donors,sphere FUNC=x*y PERIODIC=NO
383 : # Total number of hydrogen bonds donated for atoms that are in the sphere of interest
384 : numer: SUM ARG=prod PERIODIC=NO
385 : # Number of atoms that are in sphere of interest
386 : denom: SUM ARG=sphere PERIODIC=NO
387 : # Average number of donated hydrogen bonds for atoms in sphere of interest
388 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
389 : # And print out final CV to a file
390 : PRINT ARG=av FILE=colvar STRIDE=1
391 : ```
392 :
393 : This input calculate the number of hydrogen bonds that have been donated by those atoms that are within a spherical region that is centered on the point
394 : $(2.5,2.5,2.5)$.
395 :
396 : ## Optimisation details
397 :
398 : If you expand the inputs above you will see that they all the [HBPAMM_MATRIX](HBPAMM_MATRIX.md) action. The matrix that this action computes is sparse
399 : as each atom only forms hydroen bonds with a relatively small number of neighbors. The vast majority of the elements in the [HBPAMM_MATRIX](HBPAMM_MATRIX.md)
400 : are thsu zero. To reduce the amount of memory that PLUMED requires PLUMED uses sparse matrix storage. Consequently, whenever you calculate and store a contact
401 : matrix only the elements of the matrix that are non-zero are stored.
402 :
403 : We also use the sparsity of the [HBPAMM_MATRIX](HBPAMM_MATRIX.md) to make the time required to compute the matrix scale linearly rather than quadratically with
404 : the number of atoms. We know that element $i,j$ is only non-zero if the donor and acceptor atoms are within:
405 :
406 : $$
407 : r_c = \textrm{max}_i \left( \mu_i^{(ad)} + \sqrt(2 g_{cut}) | \sqrt{\lambda_{i,max}} v_\textrm{i,max}^{(ad)} | \right)
408 : $$
409 :
410 : where $i$ runs over all the PAMM kernels. $g_{cut}$ is the gaussian cutoff that is specified using the REGULARIZE
411 : keyword, $\mu_i^{(ad)}$ is the mean acceptor-donor distance for the $i$th kernel, $\lambda_{i,max}$ is the largest
412 : eigenvalue of the covariance for the $i$th kernel and $v_\textrm{i,max}^{(ad)}$ is the acceptor-donor distance component
413 : of its corresponding eigenvector. We can determine that many pairs of atoms are further apart than $r_c$ without computing the
414 : distance between these atoms by using divide and conquer strategies such as linked lists and neighbour lists. Furthermore, we do not even
415 : need to set a `D_MAX` parameter as we would have to if we were using [CONTACT_MATRIX](CONTACT_MATRIX.md) as the above cutoff can be determined
416 : from the parameters of the Gaussian kernels directly.
417 :
418 : */
419 : //+ENDPLUMEDOC
420 :
421 : //+PLUMEDOC COLVAR HBPAMM_SH
422 : /*
423 : Calculate the number of hydrogen bonds each hydrogen participates in using the HBPamm method
424 :
425 : This shortcut action allows you to calculate the number of hydrogen bonds each of the atoms specified using the SITES
426 : keyword donates to its neighbours. The number of hydrogen bonds that a particular site donates is determined by using the
427 : PAMM tehcnique that is discussed in the articles from the bibliography below and in the documentation for the [PAMM](PAMM.md) action
428 :
429 : The following example shows how you can use this action to calculate how many hydrogen bonds each of the water moecules in
430 : a box of water donates to its neighbours.
431 :
432 : ```plumed
433 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
434 : hyd: HBPAMM_SH ...
435 : SITES=1-192:3 HYDROGENS=2-192:3,3-192:3
436 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
437 : REGULARISE=0.001 GAUSS_CUTOFF=6.25
438 : ...
439 : DUMPATOMS ARG=hyd ATOMS=2-192:3,3-192:3 FILE=hydrogens.xyz
440 : ```
441 :
442 : The output here is an xyz with five columns. As explained in the documentation for [DUMPATOMS](DUMPATOMS.md), the first four
443 : columns are the usual columns that you would expect in an xyz file. The fifth column then contains the number of hydrogen bonds
444 : that each hydrogen atom participates in.
445 :
446 : In the example input above all the atoms specified using the SITE keyword can both accept and donate hydrogen bonds. If the group
447 : of atoms that can accept hydrogen bonds is different from the group of atoms that can donate hydrogen bonds then you use the ACCEPTORS and
448 : DONORS keyword as shown below:
449 :
450 : ```plumed
451 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
452 : hyd: HBPAMM_SH ...
453 : ACCEPTORS=1-9:3 DONORS=1-192:3
454 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
455 : HYDROGENS=2-192:3,3-192:3
456 : ...
457 : DUMPATOMS ARG=hyd ATOMS=2-192:3,3-192:3 FILE=hydrogens.xyz
458 : ```
459 :
460 : This input will still output an xyz file with five columns. The fifth column in this file will give the number of hydrogen bonds that each of the hydrogen
461 : atoms that were specified in this input participate within. For this atom for a hydrogen atom to particpiate in a hydrogen bond the second closest atom to
462 : it must be either atom 1, 4 or 7 as these are the only atoms that can accept hydrogen bonds according to the above input. Notice also that the
463 : as the REGULARISE and GAUSS_CUTOFF keywords were not used in this input these quantities are set to the default values that are given in the table below.
464 :
465 : ## The NOPBC flag
466 :
467 : When you use the [HBPAMM_MATRIX](HBPAMM_MATRIX.md) you need to calculate the distances between various pairs of atoms. In all the inputs that have been provided
468 : above we calculate these distances between pairs of atoms in a way that takes the periodic boundary conditions into account. If you want to ignore the periodic
469 : boundary conditions you can use the NOPBC flag as shown below.
470 :
471 : ```plumed
472 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
473 : hyd: HBPAMM_SH ...
474 : SITES=1-192:3 HYDROGENS=2-192:3,3-192:3
475 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
476 : NOPBC
477 : ...
478 : DUMPATOMS ARG=hyd ATOMS=2-192:3,3-192:3 FILE=hydrogens.xyz
479 : ```
480 :
481 : ## The MASK keyword
482 :
483 : You can use the MASK keyword with this action in the same way that it is used with [COORDINATIONNUMBER](COORDINATIONNUMBER.md). This keyword thus expects a vector in
484 : input, which tells PLUMED the atoms for which you do not need to calculate the number of hydrogen bonds. As illustrated below, this is useful if you are using functionality
485 : from the [volumes module](module_volumes.md) to calculate the average number of hydrogen bonds those hydrogen atoms that lie in a certain part of the simulation box partipipate in.
486 :
487 : ```plumed
488 : #SETTINGS INPUTFILES=regtest/pamm/rt-hbpamm/b3lyp.pamm
489 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
490 : center: FIXEDATOM AT=2.5,2.5,2.5
491 : # Vector in which element i is one if hydrogen i is in sphere of interest and zero otherwise
492 : sphere: INSPHERE ATOMS=2-192:3,3-192:3 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
493 : # Calculate the number of bonds each hydrogen participates in
494 : hyd: HBPAMM_SH ...
495 : SITES=1-192:3 HYDROGENS=2-192:3,3-192:3
496 : CLUSTERS=regtest/pamm/rt-hbpamm/b3lyp.pamm
497 : MASK=sphere
498 : ...
499 : # Multiply number of hydrogen bonds by sphere vector
500 : prod: CUSTOM ARG=hyd,sphere FUNC=x*y PERIODIC=NO
501 : # Total number of hydrogn bonds that the atoms in the sphere of interest participate in
502 : numer: SUM ARG=prod PERIODIC=NO
503 : # Number of hydrogen atoms that are in sphere of interest
504 : denom: SUM ARG=sphere PERIODIC=NO
505 : # Average number of hydrogen bonds each hydrogen atoms in sphere of interest participates in
506 : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
507 : # And print out final CV to a file
508 : PRINT ARG=av FILE=colvar STRIDE=1
509 : ```
510 :
511 : This input calculate the average number of hydrogen bonds each of the hydrogen atoms that are within a spherical region that is centered on the point
512 : $(2.5,2.5,2.5)$ participates in.
513 :
514 : ## Optimisation details
515 :
516 : If you expand the inputs above you will see that they all the [HBPAMM_MATRIX](HBPAMM_MATRIX.md) action. The matrix that this action computes is sparse
517 : as each atom only forms hydroen bonds with a relatively small number of neighbors. The vast majority of the elements in the [HBPAMM_MATRIX](HBPAMM_MATRIX.md)
518 : are thsu zero. To reduce the amount of memory that PLUMED requires PLUMED uses sparse matrix storage. Consequently, whenever you calculate and store a contact
519 : matrix only the elements of the matrix that are non-zero are stored.
520 :
521 : We also use the sparsity of the [HBPAMM_MATRIX](HBPAMM_MATRIX.md) to make the time required to compute the matrix scale linearly rather than quadratically with
522 : the number of atoms. We know that element $i,j$ is only non-zero if the donor and acceptor atoms are within:
523 :
524 : $$
525 : r_c = \textrm{max}_i \left( \mu_i^{(ad)} + \sqrt(2 g_{cut}) | \sqrt{\lambda_{i,max}} v_\textrm{i,max}^{(ad)} | \right)
526 : $$
527 :
528 : where $i$ runs over all the PAMM kernels. $g_{cut}$ is the gaussian cutoff that is specified using the REGULARIZE
529 : keyword, $\mu_i^{(ad)}$ is the mean acceptor-donor distance for the $i$th kernel, $\lambda_{i,max}$ is the largest
530 : eigenvalue of the covariance for the $i$th kernel and $v_\textrm{i,max}^{(ad)}$ is the acceptor-donor distance component
531 : of its corresponding eigenvector. We can determine that many pairs of atoms are further apart than $r_c$ without computing the
532 : distance between these atoms by using divide and conquer strategies such as linked lists and neighbour lists. Furthermore, we do not even
533 : need to set a `D_MAX` parameter as we would have to if we were using [CONTACT_MATRIX](CONTACT_MATRIX.md) as the above cutoff can be determined
534 : from the parameters of the Gaussian kernels directly.
535 :
536 : */
537 : //+ENDPLUMEDOC
538 :
539 : namespace PLMD {
540 : namespace pamm {
541 :
542 : class HBPammMatrix {
543 : public:
544 : double regulariser=0.0;
545 : Tensor incoord_to_hbcoord;
546 : std::vector<double> weight;
547 : std::vector<Vector> centers;
548 : std::vector<Tensor> kmat;
549 : static void registerKeywords( Keywords& keys );
550 : void parseInput( adjmat::AdjacencyMatrixBase<HBPammMatrix>* action );
551 : // the assignment is customized (see below): copy construction is the default one
552 6 : HBPammMatrix() = default;
553 : HBPammMatrix(const HBPammMatrix&) = default;
554 6 : HBPammMatrix& operator=( const HBPammMatrix& m ) {
555 6 : regulariser = m.regulariser;
556 6 : incoord_to_hbcoord = m.incoord_to_hbcoord;
557 6 : weight = m.weight;
558 6 : centers = m.centers;
559 6 : kmat = m.kmat;
560 6 : return *this;
561 : }
562 : static void calculateWeight( const HBPammMatrix& data, const adjmat::AdjacencyMatrixInput& input, adjmat::MatrixOutput output );
563 : };
564 :
565 : typedef adjmat::AdjacencyMatrixBase<HBPammMatrix> hbpmap;
566 : PLUMED_REGISTER_ACTION(hbpmap,"HBPAMM_MATRIX")
567 :
568 36 : void HBPammMatrix::registerKeywords( Keywords& keys ) {
569 36 : keys.use("GROUPC");
570 36 : keys.add("compulsory","ORDER","dah","the order in which the groups are specified in the input. Can be dah (donor/acceptor/hydrogens), "
571 : "adh (acceptor/donor/hydrogens) or hda (hydrogens/donor/acceptors");
572 36 : keys.add("compulsory","CLUSTERS","the name of the file that contains the definitions of all the kernels for PAMM");
573 36 : keys.add("compulsory","REGULARISE","0.001","don't allow the denominator to be smaller then this value");
574 36 : keys.add("compulsory","GAUSS_CUTOFF","6.25","the cutoff at which to stop evaluating the kernel function is set equal to sqrt(2*x)*(max(adc)+cov(adc))");
575 36 : keys.needsAction("PAMM");
576 36 : keys.needsAction("ONES");
577 36 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
578 36 : keys.addDOI("10.1063/1.4900655");
579 36 : keys.addDOI("10.1021/acs.jctc.7b00993");
580 36 : }
581 :
582 6 : void HBPammMatrix::parseInput( adjmat::AdjacencyMatrixBase<HBPammMatrix>* action ) {
583 : double DP2CUTOFF;
584 12 : action->parse("GAUSS_CUTOFF",DP2CUTOFF);
585 : std::string sorder;
586 6 : action->parse("ORDER",sorder);
587 6 : if( sorder=="dah" ) {
588 2 : incoord_to_hbcoord(0,0)=1;
589 2 : incoord_to_hbcoord(0,1)=-1;
590 2 : incoord_to_hbcoord(0,2)=0;
591 2 : incoord_to_hbcoord(1,0)=1;
592 2 : incoord_to_hbcoord(1,1)=1;
593 2 : incoord_to_hbcoord(1,2)=0;
594 2 : incoord_to_hbcoord(2,0)=0;
595 2 : incoord_to_hbcoord(2,1)=0;
596 2 : incoord_to_hbcoord(2,2)=1;
597 2 : action->log.printf(" GROUPA is list of donor atoms \n");
598 4 : } else if( sorder=="adh" ) {
599 2 : incoord_to_hbcoord(0,0)=-1;
600 2 : incoord_to_hbcoord(0,1)=1;
601 2 : incoord_to_hbcoord(0,2)=0;
602 2 : incoord_to_hbcoord(1,0)=1;
603 2 : incoord_to_hbcoord(1,1)=1;
604 2 : incoord_to_hbcoord(1,2)=0;
605 2 : incoord_to_hbcoord(2,0)=0;
606 2 : incoord_to_hbcoord(2,1)=0;
607 2 : incoord_to_hbcoord(2,2)=1;
608 2 : action->log.printf(" GROUPA is list of acceptor atoms \n");
609 2 : } else if( sorder=="hda" ) {
610 2 : incoord_to_hbcoord(0,0)=-1;
611 2 : incoord_to_hbcoord(0,1)=0;
612 2 : incoord_to_hbcoord(0,2)=1;
613 2 : incoord_to_hbcoord(1,0)=1;
614 2 : incoord_to_hbcoord(1,1)=0;
615 2 : incoord_to_hbcoord(1,2)=1;
616 2 : incoord_to_hbcoord(2,0)=0;
617 2 : incoord_to_hbcoord(2,1)=1;
618 2 : incoord_to_hbcoord(2,2)=0;
619 2 : action->log.printf(" GROUPA is list of hydrogen atoms \n");
620 : } else {
621 0 : plumed_error();
622 : }
623 : // Read in the regularisation parameter
624 12 : action->parse("REGULARISE",regulariser);
625 :
626 : // Read in the kernels
627 : double sqr2pi = sqrt(2*pi);
628 : double sqrt2pi3 = sqr2pi*sqr2pi*sqr2pi;
629 : std::string fname;
630 6 : action->parse("CLUSTERS", fname);
631 6 : double sfmax=0, ww;
632 6 : Vector cent;
633 6 : Tensor covar;
634 6 : IFile ifile;
635 6 : ifile.open(fname);
636 6 : ifile.allowIgnoredFields();
637 : while(true) {
638 144 : if( !ifile.scanField("height",ww) ) {
639 : break;
640 : }
641 66 : ifile.scanField("ptc",cent[0]);
642 66 : ifile.scanField("ssc",cent[1]);
643 66 : ifile.scanField("adc",cent[2]);
644 66 : ifile.scanField("sigma_ptc_ptc",covar[0][0]);
645 66 : ifile.scanField("sigma_ptc_ssc",covar[0][1]);
646 66 : ifile.scanField("sigma_ptc_adc",covar[0][2]);
647 66 : covar[1][0] = covar[0][1];
648 66 : ifile.scanField("sigma_ssc_ssc",covar[1][1]);
649 66 : ifile.scanField("sigma_ssc_adc",covar[1][2]);
650 66 : covar[2][0] = covar[0][2];
651 66 : covar[2][1] = covar[1][2];
652 66 : ifile.scanField("sigma_adc_adc",covar[2][2]);
653 66 : weight.push_back( ww / ( sqrt2pi3 * sqrt(covar.determinant()) ) );
654 66 : centers.push_back( cent );
655 66 : kmat.push_back( covar.inverse() );
656 :
657 66 : Vector eigval;
658 66 : Tensor eigvec;
659 66 : diagMatSym( covar, eigval, eigvec );
660 : unsigned ind_maxeval=0;
661 66 : double max_eval=eigval[0];
662 198 : for(unsigned i=1; i<3; ++i) {
663 132 : if( eigval[i]>max_eval ) {
664 : max_eval=eigval[i];
665 : ind_maxeval=i;
666 : }
667 : }
668 66 : double rcut = cent[2] + sqrt(2.0*DP2CUTOFF)*fabs(sqrt(max_eval)*eigvec(2,ind_maxeval));
669 66 : if( rcut > sfmax ) {
670 24 : sfmax = rcut;
671 : }
672 66 : ifile.scanField();
673 66 : }
674 6 : ifile.close();
675 6 : action->setLinkCellCutoff( false, sfmax );
676 12 : }
677 :
678 16316 : void HBPammMatrix::calculateWeight( const HBPammMatrix& data, const adjmat::AdjacencyMatrixInput& input, adjmat::MatrixOutput output ) {
679 16316 : Vector ddik, ddin, in_dists, hb_pamm_dists, hb_pamm_ders, real_ders;
680 16316 : ddin = input.pos;
681 16316 : in_dists[2] = ddin.modulo();
682 16316 : if( in_dists[2]<epsilon ) {
683 0 : return;
684 : }
685 :
686 16316 : output.val[0]=0;
687 16316 : Vector disp, der, tmp_der;
688 1573048 : for(unsigned i=0; i<input.natoms; ++i) {
689 1556732 : Vector ddij( input.extra_positions[i][0], input.extra_positions[i][1], input.extra_positions[i][2] );
690 1556732 : in_dists[0] = ddij.modulo();
691 1556732 : ddik = input.pbc->distance( input.pos, ddij );
692 1556732 : in_dists[1] = ddik.modulo();
693 1556732 : if( in_dists[1]<epsilon ) {
694 1554476 : continue;
695 : }
696 :
697 1548504 : hb_pamm_dists = matmul( data.incoord_to_hbcoord, in_dists );
698 1548504 : disp = hb_pamm_dists - data.centers[0];
699 1548504 : der = matmul( data.kmat[0], disp );
700 1548504 : double vv = data.weight[0]*exp( -dotProduct( disp, der ) / 2. );
701 1548504 : der *= -vv;
702 :
703 1548504 : double denom = data.regulariser + vv;
704 6194016 : for(unsigned j=0; j<3; ++j) {
705 4645512 : hb_pamm_ders[j] = der[j];
706 : }
707 17033544 : for(unsigned k=1; k<data.weight.size(); ++k) {
708 15485040 : disp = hb_pamm_dists - data.centers[k];
709 15485040 : tmp_der = matmul( data.kmat[k], disp );
710 15485040 : double tval = data.weight[k]*exp( -dotProduct( disp, tmp_der ) / 2. );
711 15485040 : denom += tval;
712 : hb_pamm_ders += -tmp_der*tval;
713 : }
714 1548504 : double vf = vv / denom;
715 1548504 : output.val[0] += vf;
716 1548504 : if( fabs(vf)<epsilon ) {
717 1546248 : continue;
718 : }
719 : // Now get derivatives
720 2256 : real_ders = matmul( der / denom - vf*hb_pamm_ders/denom, data.incoord_to_hbcoord );
721 :
722 : // And add the derivatives to the underlying atoms
723 2256 : Vector d1 = -(real_ders[0]/in_dists[0])*ddij - (real_ders[2]/in_dists[2])*ddin;
724 2256 : output.deriv[0] += d1[0];
725 2256 : output.deriv[1] += d1[1];
726 2256 : output.deriv[2] += d1[2];
727 2256 : Vector d2 = -(real_ders[1]/in_dists[1])*ddik + (real_ders[2]/in_dists[2])*ddin;
728 2256 : output.deriv[3] += d2[0];
729 2256 : output.deriv[4] += d2[1];
730 2256 : output.deriv[5] += d2[2];
731 : Vector d3 = (real_ders[0]/in_dists[0])*ddij + (real_ders[1]/in_dists[1])*ddik;
732 2256 : output.deriv[6+i*3+0] = d3[0];
733 2256 : output.deriv[6+i*3+1] = d3[1];
734 2256 : output.deriv[6+i*3+2] = d3[2];
735 2256 : Tensor vir = -(real_ders[0]/in_dists[0])*Tensor( ddij, ddij )
736 2256 : -(real_ders[1]/in_dists[1])*Tensor( ddik, ddik )
737 2256 : -(real_ders[2]/in_dists[2])*Tensor( ddin, ddin );
738 2256 : output.deriv[6 + 3*input.natoms + 0] += vir[0][0];
739 2256 : output.deriv[6 + 3*input.natoms + 1] += vir[0][1];
740 2256 : output.deriv[6 + 3*input.natoms + 2] += vir[0][2];
741 2256 : output.deriv[6 + 3*input.natoms + 3] += vir[1][0];
742 2256 : output.deriv[6 + 3*input.natoms + 4] += vir[1][1];
743 2256 : output.deriv[6 + 3*input.natoms + 5] += vir[1][2];
744 2256 : output.deriv[6 + 3*input.natoms + 6] += vir[2][0];
745 2256 : output.deriv[6 + 3*input.natoms + 7] += vir[2][1];
746 2256 : output.deriv[6 + 3*input.natoms + 8] += vir[2][2];
747 : }
748 : }
749 :
750 : class HBPammShortcut : public ActionShortcut {
751 : public:
752 : static void registerKeywords( Keywords& keys );
753 : HBPammShortcut(const ActionOptions&);
754 : };
755 :
756 : PLUMED_REGISTER_ACTION(HBPammShortcut,"HBPAMM_SD")
757 : PLUMED_REGISTER_ACTION(HBPammShortcut,"HBPAMM_SA")
758 : PLUMED_REGISTER_ACTION(HBPammShortcut,"HBPAMM_SH")
759 :
760 20 : void HBPammShortcut::registerKeywords( Keywords& keys ) {
761 20 : adjmat::AdjacencyMatrixBase<HBPammMatrix>::registerKeywords( keys );
762 20 : keys.remove("GROUP");
763 20 : keys.remove("GROUPA");
764 20 : keys.remove("GROUPB");
765 20 : keys.remove("GROUPC");
766 20 : keys.remove("ORDER");
767 20 : keys.remove("COMPONENTS");
768 40 : keys.reset_style("NL_CUTOFF","hidden");
769 40 : keys.reset_style("NL_STRIDE","hidden");
770 20 : keys.add("atoms","SITES","The list of atoms which can be part of a hydrogen bond. When this command is used the set of atoms that can donate a "
771 : "hydrogen bond is assumed to be the same as the set of atoms that can form hydrogen bonds. The atoms involved must be specified "
772 : "as a list of labels of multicolvar actions. If you would just like to use "
773 : "the atomic positions you can use a DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
774 : "other multicolvar commands is useful, however, as you can then exploit a much wider "
775 : "variety of functions of the contact matrix as described in the documentation for CONTACT_MATRIX");
776 20 : keys.add("atoms","DONORS","The list of atoms which can donate a hydrogen bond. The atoms involved must be specified "
777 : "as a list of labels of multicolvar actions. If you would just like to use "
778 : "the atomic positions you can use a DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
779 : "other multicolvar commands is useful, however, as you can then exploit a much wider "
780 : "variety of functions of the contact matrix as described in the documentation for CONTACT_MATRIX");
781 20 : keys.add("atoms","ACCEPTORS","The list of atoms which can accept a hydrogen bond. The atoms involved must be specified "
782 : "as a list of labels of multicolvar actions. If you would just like to use "
783 : "the atomic positions you can use a DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
784 : "other multicolvar commands is useful, however, as you can then exploit a much wider "
785 : "variety of functions of the contact matrix as described in the documentation for CONTACT_MATRIX");
786 20 : keys.add("atoms","HYDROGENS","The list of hydrogen atoms that can form part of a hydrogen bond. The atoms must be specified using a comma separated list, "
787 : "an index range or by using a GROUP action");
788 20 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
789 20 : keys.needsAction("HBPAMM_MATRIX");
790 40 : keys.setValueDescription("vector","a vector specifiying the number of hydrogen bonds each of the specified atoms participates within");
791 20 : }
792 :
793 6 : HBPammShortcut::HBPammShortcut(const ActionOptions&ao):
794 : Action(ao),
795 6 : ActionShortcut(ao) {
796 6 : std::string mwords = getShortcutLabel() + "_mat: HBPAMM_MATRIX";
797 6 : if( getName()=="HBPAMM_SD" ) {
798 : std::string site_str;
799 4 : parse("SITES",site_str);
800 2 : if( site_str.length()>0 ) {
801 2 : mwords += " GROUP=" + site_str;
802 : } else {
803 : std::string d_str;
804 0 : parse("DONORS",d_str);
805 0 : mwords += " GROUPA=" + d_str;
806 : std::string a_str;
807 0 : parse("ACCEPTORS",a_str);
808 0 : mwords += " GROUPB=" + a_str;
809 : }
810 : std::string h_str;
811 2 : parse("HYDROGENS",h_str);
812 2 : mwords += " GROUPC=" + h_str + " ORDER=dah";
813 4 : } else if( getName()=="HBPAMM_SA" ) {
814 : std::string site_str;
815 4 : parse("SITES",site_str);
816 2 : if( site_str.length()>0 ) {
817 2 : mwords += " GROUP=" + site_str;
818 : } else {
819 : std::string a_str;
820 0 : parse("ACCEPTORS",a_str);
821 0 : mwords += " GROUPA=" + a_str;
822 : std::string d_str;
823 0 : parse("DONORS",d_str);
824 0 : mwords += " GROUPB=" + d_str;
825 : }
826 : std::string h_str;
827 2 : parse("HYDROGENS",h_str);
828 2 : mwords += " GROUPC=" + h_str + " ORDER=adh";
829 2 : } else if( getName()=="HBPAMM_SH" ) {
830 : std::string h_str;
831 2 : parse("HYDROGENS",h_str);
832 2 : mwords += " GROUPA=" + h_str + " ORDER=hda";
833 : std::string site_str;
834 4 : parse("SITES",site_str);
835 2 : if( site_str.length()>0 ) {
836 2 : mwords += " GROUPB=" + site_str;
837 2 : mwords += " GROUPC=" + site_str;
838 : } else {
839 : std::string d_str;
840 0 : parse("DONORS",d_str);
841 0 : mwords += " GROUPB=" + d_str;
842 : std::string a_str;
843 0 : parse("ACCEPTORS",a_str);
844 0 : mwords += " GROUPC=" + a_str;
845 : }
846 : }
847 : std::map<std::string,std::string> keymap;
848 6 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
849 12 : readInputLine( mwords + " " + convertInputLineToString() );
850 6 : ActionWithValue* mb=plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
851 6 : plumed_assert( mb );
852 : std::string nsize;
853 6 : Tools::convert( (mb->copyOutput(getShortcutLabel() + "_mat"))->getShape()[1], nsize );
854 6 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + nsize );
855 18 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat," + getShortcutLabel() + "_ones");
856 12 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
857 20 : }
858 :
859 : }
860 : }
|
