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 "AdjacencyMatrixBase.h"
23 : #include "tools/SwitchingFunction.h"
24 : #include "tools/Angle.h"
25 : #include "core/ActionRegister.h"
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : //+PLUMEDOC MATRIX HBOND_MATRIX
31 : /*
32 : Adjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them.
33 :
34 : A useful tool for developing complex collective variables is the notion of the
35 : so called adjacency matrix. An adjacency matrix is an $N \times N$ matrix in which the $i$th, $j$th element tells you whether
36 : or not the $i$th and $j$th atoms/molecules from a set of $N$ atoms/molecules are adjacent or not. As detailed in the documentation
37 : for [CONTACT_MATRIX](CONTACT_MATRIX.md) there are then a range of further analyses that you can perform on these matrices.
38 :
39 : For this action the elements of the adjacency matrix are calculated using:
40 :
41 : $$
42 : a_{ij} = \sigma_{oo}( |\mathbf{r}_{ij}| ) \sum_{k=1}^N \sigma_{oh}( |\mathbf{r}_{ik}| ) \sigma_{\theta}( \theta_{kij} )
43 : $$
44 :
45 : This expression was derived by thinking about how to detect if there is a hydrogen bond between atoms $i$ and $j$. The notion is that
46 : if the hydrogen bond is present atoms $i$ and $j$ should be within a certain cutoff distance. In addition, there should be a hydrogen
47 : within a certain cutoff distance of atom $i$ and this hydrogen should lie on or close to the vector connecting atoms $i$ and $j$.
48 : As such $\sigma_{oo}(r_{ij})$ is a switching function that acts on the modulus of the vector connecting atom $i$ to atom
49 : $j$. The sum over $k$ then runs over all the hydrogen atoms that are specified using using HYDROGEN keyword. $\sigma_{oh}(r_{ik})$
50 : is a switching function that acts on the modulus of the vector connecting atom $i$ to atom $k$ and $\sigma_{\theta}(\theta_{kij})$
51 : is a switching function that acts on the angle between the vector connecting atoms $i$ and $j$ and the vector connecting atoms $i$ and
52 : $k$.
53 :
54 : It is important to note that hydrogen bonds, unlike regular bonds, are asymmetric. In other words, the hydrogen atom does not sit at the
55 : mid point between the two other atoms in this three-center bond. As a result of this adjacency matrices calculated using HBOND_MATRIX are not
56 : symmetric like those calculated by [CONTACT_MATRIX](CONTACT_MATRIX.md).
57 :
58 : Each water molecule can participate in a hydrogen bond in one of two ways. It can either donate one of its hydrogen atom to the neighboring oxygen or
59 : it can accept a bond between the hydrogen of a neighboring water molecule and its own oxygen.
60 :
61 : The following input can be used to analyze the number of hydrogen bonds each of the oxygen atoms in a box of water donates. This information is output in an
62 : xyz files which contains five columns of data. The first four of these columns are a label for the atom and the x, y and z position of the oxygen. The last column is then
63 : the number of hydrogen bond that water molecule donates.
64 :
65 : ```plumed
66 : mat: HBOND_MATRIX ...
67 : DONORS=1-192:3 ACCEPTORS=1-192:3 HYDROGENS=2-192:3,3-192:3
68 : SWITCH={RATIONAL R_0=3.20} HSWITCH={RATIONAL R_0=2.30}
69 : ASWITCH={RATIONAL R_0=0.167pi}
70 : ...
71 : ones: ONES SIZE=64
72 : rsums: MATRIX_VECTOR_PRODUCT ARG=mat,ones
73 : DUMPATOMS ATOMS=1-192:3 ARG=rsums FILE=donors.xyz
74 : ```
75 :
76 : If you want to calculate the number of hydorgen bonds each of the oxygen atoms accepts you would use an input like the one below:
77 :
78 : ```plumed
79 : mat: HBOND_MATRIX ...
80 : DONORS=1-192:3 ACCEPTORS=1-192:3 HYDROGENS=2-192:3,3-192:3
81 : SWITCH={RATIONAL R_0=3.20} HSWITCH={RATIONAL R_0=2.30}
82 : ASWITCH={RATIONAL R_0=0.167pi}
83 : ...
84 : matT: TRANSPOSE ARG=mat
85 : ones: ONES SIZE=64
86 : rsums: MATRIX_VECTOR_PRODUCT ARG=matT,ones
87 : DUMPATOMS ATOMS=1-192:3 ARG=rsums FILE=acceptors.xyz
88 : ```
89 :
90 : Consequently, if you want the total number of hydrogen bonds each oxygen atom participates in you need to use the following input:
91 :
92 : ```plumed
93 : mat: HBOND_MATRIX ...
94 : DONORS=1-192:3 ACCEPTORS=1-192:3 HYDROGENS=2-192:3,3-192:3
95 : SWITCH={RATIONAL R_0=3.20} HSWITCH={RATIONAL R_0=2.30}
96 : ASWITCH={RATIONAL R_0=0.167pi}
97 : ...
98 : matT: TRANSPOSE ARG=mat
99 : hbmat: CUSTOM ARG=mat,matT FUNC=x+y PERIODIC=NO
100 : ones: ONES SIZE=64
101 : rsums: MATRIX_VECTOR_PRODUCT ARG=hbmat,ones
102 : DUMPATOMS ATOMS=1-192:3 ARG=rsums FILE=hbonds.xyz
103 : ```
104 :
105 : Notice that in all the inputs above the $r_{ij}$ and $r_{ik}$ values that enter the formula above are calculated in a way that takes the
106 : periodic boundary conditions into account. If you want to ignore the periodic boundary conditions you can use the NOPBC flag as shown below.
107 :
108 : ```plumed
109 : mat: HBOND_MATRIX ...
110 : DONORS=1-192:3 ACCEPTORS=1-192:3 HYDROGENS=2-192:3,3-192:3
111 : SWITCH={RATIONAL R_0=3.20} HSWITCH={RATIONAL R_0=2.30}
112 : ASWITCH={RATIONAL R_0=0.167pi}
113 : NOPBC
114 : ...
115 : ```
116 :
117 : ## COMPONENTS flag
118 :
119 : If you add the flag COMPONENTS to the input as shown below:
120 :
121 : ```plumed
122 : c4: HBOND_MATRIX ...
123 : DONORS=1-192:3 ACCEPTORS=1-192:3 HYDROGENS=2-192:3,3-192:3
124 : SWITCH={RATIONAL R_0=3.20} HSWITCH={RATIONAL R_0=2.30}
125 : ASWITCH={RATIONAL R_0=0.167pi}
126 : COMPONENTS
127 : ...
128 : ```
129 :
130 : 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
131 : 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$
132 : 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`
133 : if the elements of `c4.w` are non-zero. Using the COMPONENTS flag in this way ensures that you can use HBOND_MATRIX in tandem with many of the functionalities
134 : that are part of the [symfunc module](module_symfunc.md). Remember, however, that the $i,j$ element of the HBOND_MATRIX is only non-zero if atom $i$ donates
135 : a hydrogen bond to atom $j$. __You cannot use HBOND_MATRIX to identify the set of atoms that each atom is hydrogen bonded to.__
136 :
137 : ## The MASK keyword
138 :
139 : You use the MASK keyword with HBOND_MATRIX in the same way that is used in [CONTACT_MATRIX](CONTACT_MATRIX.md). This keyword thus expects a vector in input,
140 : 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:
141 :
142 : ```plumed
143 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
144 : center: FIXEDATOM AT=2.5,2.5,2.5
145 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
146 : sphere: INSPHERE ATOMS=1-192:3 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
147 : # Calculates cooordination numbers
148 : cmap: HBOND_MATRIX ...
149 : DONORS=1-192:3 ACCEPTORS=1-192:3 HYDROGENS=2-192:3,3-192:3
150 : SWITCH={RATIONAL R_0=3.20} HSWITCH={RATIONAL R_0=2.30}
151 : ASWITCH={RATIONAL R_0=0.167pi} MASK=sphere
152 : ...
153 : ones: ONES SIZE=64
154 : cc: MATRIX_VECTOR_PRODUCT ARG=cmap,ones
155 : # Multiply coordination numbers by sphere vector
156 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
157 : # Sum of coordination numbers for atoms that are in the sphere of interest
158 : numer: SUM ARG=prod PERIODIC=NO
159 : # Number of atoms that are in sphere of interest
160 : denom: SUM ARG=sphere PERIODIC=NO
161 : # Average coordination number for atoms in sphere of interest
162 : av: CUSTOM ARG=prod,sphere FUNC=x/y PERIODIC=NO
163 : # And print out final CV to a file
164 : PRINT ARG=av FILE=colvar STRIDE=1
165 : ```
166 :
167 : This input calculates the average number of hydrogen bonds each of the aatoms that are within a spherical region that is centered on the point
168 : $(2.5,2.5,2.5)$ donate.
169 :
170 : ## Optimisation details
171 :
172 : Adjacency matrices are sparse. Each atom is only be connected to a small number of neighbours and the vast majority of the elements of the contact matrix are thus zero. To reduce
173 : the amount of memory that PLUMED requires PLUMED uses sparse matrix storage. Consequently, whenever you calculate and store a contact matrix only the elements of the matrix that are
174 : non-zero are stored. The same thing holds for the additional values that are created when you use the COMPONENTS flag. The components of the vectors connecting atoms are only stored
175 : when the elements of `c4.w` are non-zero.
176 :
177 : We can also use the sparsity of the adjacency matrix to make the time required to compute a contact matrix scale linearly rather than quadratically with the number of atoms. Element
178 : $i,j$ of the contact matrix is only non-zero if two atoms are within a cutoff, $r_c$. We can determine that many pairs of atoms are further appart than $r_c$ without computing the
179 : distance between these atoms by using divide and conquer strategies such as linked lists and neighbour lists. __To turn on these features you need to set the `D_MAX` parameter in the
180 : switching functions.__ The value you pass to the `D_MAX` keyword is used as the cutoff in the link cell algorithm.
181 :
182 : In theory we could further optimize the implementation of the HBOND_MATRIX action by exploiting neighbor lists. If we were to do this we would likely add two further keywords as shown
183 : below:
184 :
185 : ```plumed
186 : cmap: HBOND_MATRIX ...
187 : DONORS=1-192:3 ACCEPTORS=1-192:3 HYDROGENS=2-192:3,3-192:3
188 : SWITCH={RATIONAL R_0=3.20 D_MAX=4.0} HSWITCH={RATIONAL R_0=2.30 D_MAX=3.5}
189 : ASWITCH={RATIONAL R_0=0.167pi}
190 : NL_CUTOFF=5.0 NL_STRIDE=5
191 : ...
192 : ```
193 :
194 : The `NL_CUTOFF` keyword would be used to specify the cutoff (in nm) to use when constructing neighbor lists. This value would need to be slightly larger than the D_MAX parameter for the switching function that
195 : acts on the acceptor donor distance.
196 : The `NL_STRIDE` keyword would then be used to specify how frequently the neighbour list should be updated. Thus far we have not found it necessary to implement this algorithm. We have been happy with the
197 : performance even if we use the linked list algorithm to update the neighbors on every step. If you feel that you need this CV to perform better please get in touch as adding a neighbor list for this action
198 : should be relatively straightforward.
199 :
200 :
201 : */
202 : //+ENDPLUMEDOC
203 :
204 : namespace PLMD {
205 : namespace adjmat {
206 :
207 : class HbondMatrix {
208 : public:
209 : SwitchingFunction distanceOOSwitch;
210 : SwitchingFunction distanceOHSwitch;
211 : SwitchingFunction angleSwitch;
212 : static void registerKeywords( Keywords& keys );
213 : void parseInput( AdjacencyMatrixBase<HbondMatrix>* action );
214 : static void calculateWeight( const HbondMatrix& data,
215 : const AdjacencyMatrixInput& input,
216 : MatrixOutput output );
217 : };
218 :
219 : typedef AdjacencyMatrixBase<HbondMatrix> hmap;
220 : PLUMED_REGISTER_ACTION(hmap,"HBOND_MATRIX")
221 :
222 4 : void HbondMatrix::registerKeywords( Keywords& keys ) {
223 8 : keys.reset_style("GROUP","deprecated");
224 4 : keys.remove("ATOMS");
225 4 : keys.remove("GROUPA");
226 8 : keys.remove("GROUPB");
227 4 : keys.add("atoms-2","DONORS","The list of atoms which can donate a hydrogen bond");
228 4 : keys.add("atoms-2","ACCEPTORS","The list of atoms which can accept a hydrogen bond");
229 4 : keys.add("atoms","HYDROGENS","The list of atoms that can form the bridge between the two interesting parts "
230 : "of the structure.");
231 4 : keys.add("numbered","SWITCH","The switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them");
232 8 : keys.linkActionInDocs("SWITCH","LESS_THAN");
233 4 : keys.add("numbered","HSWITCH","The switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be "
234 : "considered a hydrogen bond");
235 8 : keys.linkActionInDocs("HSWITCH","LESS_THAN");
236 4 : keys.add("numbered","ASWITCH","A switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and "
237 : "the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond");
238 8 : keys.linkActionInDocs("ASWITCH","LESS_THAN");
239 4 : }
240 :
241 2 : void HbondMatrix::parseInput( AdjacencyMatrixBase<HbondMatrix>* action ) {
242 : std::string errors;
243 : std::string OOinput;
244 4 : action->parse("SWITCH",OOinput);
245 2 : if( OOinput.length()==0 ) {
246 0 : action->error("could not find SWITCH keyword");
247 : }
248 2 : distanceOOSwitch.set(OOinput,errors);
249 2 : if( errors.length()!=0 ) {
250 0 : action->error("problem reading SWITCH keyword : " + errors );
251 : }
252 :
253 : std::string OHinput;
254 4 : action->parse("HSWITCH",OHinput);
255 2 : if( OHinput.length()==0 ) {
256 0 : action->error("could not find HSWITCH keyword");
257 : }
258 2 : distanceOHSwitch.set(OHinput,errors);
259 2 : if( errors.length()!=0 ) {
260 0 : action->error("problem reading HSWITCH keyword : " + errors );
261 : }
262 :
263 : std::string anginput;
264 4 : action->parse("ASWITCH",anginput);
265 2 : if( anginput.length()==0 ) {
266 0 : action->error("could not find SWITCH keyword");
267 : }
268 2 : angleSwitch.set(anginput,errors);
269 2 : if( errors.length()!=0 ) {
270 0 : action->error("problem reading SWITCH keyword : " + errors );
271 : }
272 :
273 : // Setup link cells
274 2 : action->setLinkCellCutoff( false, distanceOOSwitch.get_dmax() );
275 2 : }
276 :
277 4108 : void HbondMatrix::calculateWeight( const HbondMatrix& data,
278 : const AdjacencyMatrixInput& input,
279 : MatrixOutput output ) {
280 4108 : Vector ood = input.pos;
281 : double ood_l = ood.modulo2(); // acceptor - donor
282 4108 : if( ood_l<epsilon) {
283 64 : return;
284 : }
285 : double ood_df;
286 4044 : double ood_sw=data.distanceOOSwitch.calculateSqr( ood_l, ood_df );
287 :
288 520212 : for(unsigned i=0; i<input.natoms; ++i) {
289 : Vector ohd(input.extra_positions[i][0],
290 : input.extra_positions[i][1],
291 516168 : input.extra_positions[i][2]);
292 : double ohd_l=ohd.modulo2();
293 : double ohd_df;
294 516168 : double ohd_sw=data.distanceOHSwitch.calculateSqr( ohd_l, ohd_df );
295 :
296 : Angle a;
297 : Vector ood_adf;
298 : Vector ohd_adf;
299 516168 : double angle=a.compute( ood, ohd, ood_adf, ohd_adf );
300 : double angle_df;
301 516168 : double angle_sw=data.angleSwitch.calculate( angle, angle_df );
302 516168 : output.val[0] += ood_sw*ohd_sw*angle_sw;
303 :
304 516168 : if( !input.noderiv ) {
305 36 : Vector d1 = angle_sw*ohd_sw*(-ood_df)*ood
306 36 : + angle_sw*ood_sw*(-ohd_df)*ohd
307 36 : + ood_sw*ohd_sw*angle_df*angle*(-ood_adf-ohd_adf);
308 36 : output.deriv[0] += d1[0];
309 36 : output.deriv[1] += d1[1];
310 36 : output.deriv[2] += d1[2];
311 : Vector d2 = angle_sw*ohd_sw*(+ood_df)*ood
312 : + ood_sw*ohd_sw*angle_df*angle*ood_adf;
313 36 : output.deriv[3] += d2[0];
314 36 : output.deriv[4] += d2[1];
315 36 : output.deriv[5] += d2[2];
316 : Vector d3 = angle_sw*ood_sw*(+ohd_df)*ohd
317 : + ood_sw*ohd_sw*angle_df*angle*ohd_adf;
318 36 : output.deriv[6+i*3+0] = d3[0];
319 36 : output.deriv[6+i*3+1] = d3[1];
320 36 : output.deriv[6+i*3+2] = d3[2];
321 36 : Tensor vir = angle_sw*ohd_sw*(-ood_df)*Tensor(ood,ood)
322 36 : + angle_sw*ood_sw*(-ohd_df)*Tensor(ohd,ohd)
323 36 : - ood_sw*ohd_sw*angle_df*angle*(Tensor(ood,ood_adf)
324 36 : + Tensor(ohd,ohd_adf));
325 36 : output.deriv[6 + 3*input.natoms + 0] += vir[0][0];
326 36 : output.deriv[6 + 3*input.natoms + 1] += vir[0][1];
327 36 : output.deriv[6 + 3*input.natoms + 2] += vir[0][2];
328 36 : output.deriv[6 + 3*input.natoms + 3] += vir[1][0];
329 36 : output.deriv[6 + 3*input.natoms + 4] += vir[1][1];
330 36 : output.deriv[6 + 3*input.natoms + 5] += vir[1][2];
331 36 : output.deriv[6 + 3*input.natoms + 6] += vir[2][0];
332 36 : output.deriv[6 + 3*input.natoms + 7] += vir[2][1];
333 36 : output.deriv[6 + 3*input.natoms + 8] += vir[2][2];
334 : }
335 : }
336 : }
337 :
338 : }
339 : }
|
