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 : The following input can be used to analyze the number of hydrogen bonds each of the oxygen atoms in a box of water participates in. Each
59 : 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
60 : it can accept a bond between the hydrogen of a neighboring water molecule and its own oxygen. The input below allows you to output information
61 : on the number of hydrogen bonds each of the water molecules donates and accepts. This information is output in two xyz files which each contain
62 : 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 accepted/donated hydrogen bonds.
64 :
65 : ```plumed
66 : mat: HBOND_MATRIX ATOMS=1-192:3 HYDROGENS=2-192:3,3-192:3 SWITCH={RATIONAL R_0=3.20} HSWITCH={RATIONAL R_0=2.30} ASWITCH={RATIONAL R_0=0.167pi}
67 : ones: ONES SIZE=64
68 : rsums: MATRIX_VECTOR_PRODUCT ARG=mat,ones
69 : DUMPATOMS ATOMS=1-192:3 ARG=rsums FILE=donors.xyz
70 : ```
71 :
72 : */
73 : //+ENDPLUMEDOC
74 :
75 : namespace PLMD {
76 : namespace adjmat {
77 :
78 : class HbondMatrix : public AdjacencyMatrixBase {
79 : private:
80 : SwitchingFunction distanceOOSwitch;
81 : SwitchingFunction distanceOHSwitch;
82 : SwitchingFunction angleSwitch;
83 : public:
84 : static void registerKeywords( Keywords& keys );
85 : explicit HbondMatrix(const ActionOptions&);
86 : // active methods:
87 : double calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const override ;
88 : };
89 :
90 : PLUMED_REGISTER_ACTION(HbondMatrix,"HBOND_MATRIX")
91 :
92 4 : void HbondMatrix::registerKeywords( Keywords& keys ) {
93 4 : AdjacencyMatrixBase::registerKeywords( keys );
94 4 : keys.add("atoms-2","DONORS","The list of atoms which can donate a hydrogen bond");
95 4 : keys.add("atoms-2","ACCEPTORS","The list of atoms which can accept a hydrogen bond");
96 4 : keys.add("atoms","HYDROGENS","The list of atoms that can form the bridge between the two interesting parts "
97 : "of the structure.");
98 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");
99 8 : keys.linkActionInDocs("SWITCH","LESS_THAN");
100 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 "
101 : "considered a hydrogen bond");
102 8 : keys.linkActionInDocs("HSWITCH","LESS_THAN");
103 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 "
104 : "the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond");
105 8 : keys.linkActionInDocs("ASWITCH","LESS_THAN");
106 4 : }
107 :
108 2 : HbondMatrix::HbondMatrix(const ActionOptions&ao):
109 : Action(ao),
110 2 : AdjacencyMatrixBase(ao) {
111 : std::string sfinput,errors;
112 4 : parse("SWITCH",sfinput);
113 2 : if( sfinput.length()==0 ) {
114 0 : error("could not find SWITCH keyword");
115 : }
116 2 : distanceOOSwitch.set(sfinput,errors);
117 2 : if( errors.length()!=0 ) {
118 0 : error("problem reading SWITCH keyword : " + errors );
119 : }
120 :
121 : std::string hsfinput;
122 4 : parse("HSWITCH",hsfinput);
123 2 : if( hsfinput.length()==0 ) {
124 0 : error("could not find HSWITCH keyword");
125 : }
126 2 : distanceOHSwitch.set(hsfinput,errors);
127 2 : if( errors.length()!=0 ) {
128 0 : error("problem reading HSWITCH keyword : " + errors );
129 : }
130 :
131 : std::string asfinput;
132 4 : parse("ASWITCH",asfinput);
133 2 : if( asfinput.length()==0 ) {
134 0 : error("could not find SWITCH keyword");
135 : }
136 2 : angleSwitch.set(asfinput,errors);
137 2 : if( errors.length()!=0 ) {
138 0 : error("problem reading SWITCH keyword : " + errors );
139 : }
140 :
141 : // Setup link cells
142 2 : setLinkCellCutoff( false, distanceOOSwitch.get_dmax() );
143 :
144 : // And check everything has been read in correctly
145 2 : checkRead();
146 2 : }
147 :
148 4038 : double HbondMatrix::calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const {
149 4038 : Vector ood = pos2;
150 4038 : double ood_l = ood.modulo2(); // acceptor - donor
151 4038 : if( ood_l<epsilon) {
152 : return 0;
153 : }
154 4038 : double ood_df, ood_sw=distanceOOSwitch.calculateSqr( ood_l, ood_df );
155 :
156 : double value=0;
157 520170 : for(unsigned i=0; i<natoms; ++i) {
158 516132 : Vector ohd=getPosition(i,myvals);
159 516132 : double ohd_l=ohd.modulo2();
160 516132 : double ohd_df, ohd_sw=distanceOHSwitch.calculateSqr( ohd_l, ohd_df );
161 :
162 : Angle a;
163 516132 : Vector ood_adf, ohd_adf;
164 516132 : double angle=a.compute( ood, ohd, ood_adf, ohd_adf );
165 516132 : double angle_df, angle_sw=angleSwitch.calculate( angle, angle_df );
166 516132 : value += ood_sw*ohd_sw*angle_sw;
167 :
168 516132 : if( !doNotCalculateDerivatives() ) {
169 36 : addAtomDerivatives( 0, angle_sw*ohd_sw*(-ood_df)*ood + angle_sw*ood_sw*(-ohd_df)*ohd + ood_sw*ohd_sw*angle_df*angle*(-ood_adf-ohd_adf), myvals );
170 36 : addAtomDerivatives( 1, angle_sw*ohd_sw*(+ood_df)*ood + ood_sw*ohd_sw*angle_df*angle*ood_adf, myvals );
171 36 : addThirdAtomDerivatives( i, angle_sw*ood_sw*(+ohd_df)*ohd + ood_sw*ohd_sw*angle_df*angle*ohd_adf, myvals );
172 72 : addBoxDerivatives( angle_sw*ohd_sw*(-ood_df)*Tensor(ood,ood) + angle_sw*ood_sw*(-ohd_df)*Tensor(ohd,ohd) -
173 108 : ood_sw*ohd_sw*angle_df*angle*(Tensor(ood,ood_adf)+Tensor(ohd,ohd_adf)), myvals );
174 : }
175 : }
176 4038 : return value;
177 : }
178 :
179 : }
180 : }
|
