Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2023 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 "multicolvar/AtomValuePack.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/SwitchingFunction.h"
26 : #include "tools/HistogramBead.h"
27 : #include "tools/Angle.h"
28 : #include "tools/Matrix.h"
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 : As discussed in the section of the manual on \ref contactmatrix a useful tool for developing complex collective variables is the notion of the
35 : so called adjacency matrix. An adjacency matrix is an \f$N \times N\f$ matrix in which the \f$i\f$th, \f$j\f$th element tells you whether
36 : or not the \f$i\f$th and \f$j\f$th atoms/molecules from a set of \f$N\f$ atoms/molecules are adjacent or not. These matrices can then be further
37 : analyzed using a number of other algorithms as is detailed in \cite tribello-clustering.
38 :
39 : For this action the elements of the adjacency matrix are calculated using:
40 :
41 : \f[
42 : a_{ij} = \sigma_{oo}( |\mathbf{r}_{ij}| ) \sum_{k=1}^N \sigma_{oh}( |\mathbf{r}_{ik}| ) \sigma_{\theta}( \theta_{kij} )
43 : \f]
44 :
45 : This expression was derived by thinking about how to detect if there is a hydrogen bond between atoms \f$i\f$ and \f$j\f$. The notion is that
46 : if the hydrogen bond is present atoms \f$i\f$ and \f$j\f$ should be within a certain cutoff distance. In addition, there should be a hydrogen
47 : within a certain cutoff distance of atom \f$i\f$ and this hydrogen should lie on or close to the vector connecting atoms \f$i\f$ and \f$j\f$.
48 : As such \f$\sigma_{oo}( |\mathbf{r}_{ij}| )\f$ is a \ref switchingfunction that acts on the modulus of the vector connecting atom \f$i\f$ to atom
49 : \f$j\f$. The sum over \f$k\f$ then runs over all the hydrogen atoms that are specified using using HYDROGEN keyword. \f$\sigma_{oh}(|\mathbf{r}_{ik}|)\f$
50 : is a \ref switchingfunction that acts on the modulus of the vector connecting atom \f$i\f$ to atom \f$k\f$ and \f$\sigma_{\theta}(\theta_{kij})\f$
51 : is a \ref switchingfunction that acts on the angle between the vector connecting atoms \f$i\f$ and \f$j\f$ and the vector connecting atoms \f$i\f$ and
52 : \f$k\f$.
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 \ref HBOND_MATRIX are not
56 : symmetric like those calculated by \ref CONTACT_MATRIX. One consequence of this fact is that the quantities found by performing \ref ROWSUMS and
57 : \ref COLUMNSUMS on a square \ref HBOND_MATRIX are not the same as they would be if you performed \ref ROWSUMS and
58 : \ref COLUMNSUMS on a square \ref CONTACT_MATRIX.
59 :
60 : \par Examples
61 :
62 : 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
63 : 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
64 : 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
65 : 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
66 : 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
67 : the number of accepted/donated hydrogen bonds.
68 :
69 : \plumedfile
70 : 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} SUM
71 : rsums: ROWSUMS MATRIX=mat MEAN
72 : csums: COLUMNSUMS MATRIX=mat MEAN
73 : DUMPMULTICOLVAR DATA=rsums FILE=donors.xyz
74 : DUMPMULTICOLVAR DATA=csums FILE=acceptors.xyz
75 : \endplumedfile
76 :
77 : */
78 : //+ENDPLUMEDOC
79 :
80 :
81 : namespace PLMD {
82 : namespace adjmat {
83 :
84 : class HBondMatrix : public AdjacencyMatrixBase {
85 : private:
86 : unsigned ndonor_types;
87 : /// switching function
88 : Matrix<SwitchingFunction> distanceOOSwitch;
89 : Matrix<SwitchingFunction> distanceOHSwitch;
90 : Matrix<SwitchingFunction> angleSwitch;
91 : public:
92 : /// Create manual
93 : static void registerKeywords( Keywords& keys );
94 : /// Constructor
95 : explicit HBondMatrix(const ActionOptions&);
96 : /// Create the ith, ith switching function
97 : void setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) override;
98 : /// This actually calculates the value of the contact function
99 : double calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const override;
100 : /// This does nothing
101 : double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const override;
102 : ///
103 : double calculateForThree( const unsigned& iat, const unsigned& ano, const unsigned& dno, const Vector& ood,
104 : const double& ood_df, const double& ood_sw, multicolvar::AtomValuePack& myatoms ) const;
105 : };
106 :
107 13789 : PLUMED_REGISTER_ACTION(HBondMatrix,"HBOND_MATRIX")
108 :
109 6 : void HBondMatrix::registerKeywords( Keywords& keys ) {
110 6 : AdjacencyMatrixBase::registerKeywords( keys );
111 12 : keys.add("atoms","ATOMS","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 "
112 : "hydrogen bond is assumed to be the same as the set of atoms that can form hydrogen bonds. The atoms involved must be specified "
113 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
114 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
115 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
116 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
117 12 : 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, "
118 : "an index range or by using a \\ref GROUP. A list of hydrogen atoms is always required even if you specify the other atoms using "
119 : "DONORS and ACCEPTORS as described below.");
120 12 : keys.add("atoms-2","DONORS","The list of atoms which can donate a hydrogen bond. The atoms involved must be specified "
121 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
122 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
123 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
124 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
125 12 : keys.add("atoms-2","ACCEPTORS","The list of atoms which can accept a hydrogen bond. The atoms involved must be specified "
126 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
127 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
128 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
129 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
130 12 : keys.add("numbered","SWITCH","The \\ref switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them");
131 12 : keys.add("numbered","HSWITCH","The \\ref switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be "
132 : "considered a hydrogen bond");
133 12 : keys.add("numbered","ASWITCH","A \\ref switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and "
134 : "the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond");
135 6 : keys.use("SUM");
136 6 : }
137 :
138 2 : HBondMatrix::HBondMatrix( const ActionOptions& ao ):
139 : Action(ao),
140 2 : AdjacencyMatrixBase(ao) {
141 4 : readMaxThreeSpeciesMatrix( "ATOMS", "DONORS", "ACCEPTORS", "HYDROGENS", false );
142 : unsigned nrows, ncols;
143 2 : retrieveTypeDimensions( nrows, ncols, ndonor_types );
144 2 : distanceOOSwitch.resize( nrows, ncols );
145 2 : distanceOHSwitch.resize( nrows, ncols );
146 2 : angleSwitch.resize( nrows, ncols );
147 2 : parseConnectionDescriptions("SWITCH",false,ndonor_types);
148 2 : parseConnectionDescriptions("HSWITCH",false,ndonor_types);
149 4 : parseConnectionDescriptions("ASWITCH",false,ndonor_types);
150 :
151 : // Find the largest sf cutoff
152 2 : double sfmax=distanceOOSwitch(0,0).get_dmax();
153 4 : for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
154 4 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
155 2 : double tsf=distanceOOSwitch(i,j).get_dmax();
156 2 : if( tsf>sfmax ) {
157 0 : sfmax=tsf;
158 : }
159 : }
160 : }
161 : // Set the link cell cutoff
162 2 : setLinkCellCutoff( sfmax );
163 2 : }
164 :
165 6 : void HBondMatrix::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
166 6 : plumed_assert( id<3 && desc.size()==1 );
167 6 : if( id==0 ) {
168 : std::string errors;
169 2 : distanceOOSwitch(j,i).set(desc[0],errors);
170 2 : if( errors.length()!=0 ) {
171 0 : error("problem reading switching function description " + errors);
172 : }
173 2 : if( j!=i) {
174 0 : distanceOOSwitch(i,j).set(desc[0],errors);
175 : }
176 4 : log.printf(" atoms of type %u and %u must be within %s\n",i+1,j+1,(distanceOOSwitch(i,j).description()).c_str() );
177 4 : } else if( id==1 ) {
178 : std::string errors;
179 2 : distanceOHSwitch(j,i).set(desc[0],errors);
180 2 : if( errors.length()!=0 ) {
181 0 : error("problem reading switching function description " + errors);
182 : }
183 2 : if( j!=i) {
184 0 : distanceOHSwitch(i,j).set(desc[0],errors);
185 : }
186 4 : log.printf(" for atoms of type %u and %u the OH distance must be less than %s \n",i+1,j+1,(distanceOHSwitch(i,j).description()).c_str() );
187 2 : } else if( id==2 ) {
188 : std::string errors;
189 2 : angleSwitch(j,i).set(desc[0],errors);
190 2 : if( errors.length()!=0 ) {
191 0 : error("problem reading switching function description " + errors);
192 : }
193 2 : if( j!=i) {
194 0 : angleSwitch(i,j).set(desc[0],errors);
195 : }
196 4 : log.printf(" for atoms of type %u and %u the OOH angle must be less than %s \n",i+1,j+1,(angleSwitch(i,j).description()).c_str() );
197 : }
198 6 : }
199 :
200 4038 : double HBondMatrix::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
201 : // Ensure we skip diagonal elements of square matrix
202 4038 : if( myatoms.getIndex(0)==myatoms.getIndex(1) ) {
203 : return 0.0;
204 : }
205 :
206 4038 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
207 8076 : if( distance.modulo2()<distanceOOSwitch( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax2() ) {
208 4038 : return 1.0;
209 : }
210 : return 0.0;
211 : }
212 :
213 4038 : double HBondMatrix::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
214 4038 : Vector ood = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
215 4038 : double ood_l = ood.modulo(); // acceptor - donor
216 : double ood_df, ood_sw=distanceOOSwitch( getBaseColvarNumber( myatoms.getIndex(0) ),
217 4038 : getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( ood_l, ood_df );
218 :
219 : // Get the base colvar numbers
220 4038 : unsigned ano, dno = getBaseColvarNumber( myatoms.getIndex(0) );
221 4038 : if( ndonor_types==0 ) {
222 4038 : ano = getBaseColvarNumber( myatoms.getIndex(1) );
223 : } else {
224 0 : ano = getBaseColvarNumber( myatoms.getIndex(1) ) - ndonor_types;
225 : }
226 :
227 : double value=0;
228 4038 : if( myatoms.getNumberOfAtoms()>3 ) {
229 520170 : for(unsigned i=2; i<myatoms.getNumberOfAtoms(); ++i) {
230 516132 : value+=calculateForThree( i, ano, dno, ood, ood_df, ood_sw, myatoms );
231 : }
232 : } else {
233 : plumed_dbg_assert( myatoms.getNumberOfAtoms()==3 );
234 0 : value=calculateForThree( 2, ano, dno, ood, ood_df, ood_sw, myatoms );
235 : }
236 4038 : return value;
237 : }
238 :
239 516132 : double HBondMatrix::calculateForThree( const unsigned& iat, const unsigned& ano, const unsigned& dno, const Vector& ood,
240 : const double& ood_df, const double& ood_sw, multicolvar::AtomValuePack& myatoms ) const {
241 516132 : Vector ohd=getSeparation( myatoms.getPosition(0), myatoms.getPosition(iat) );
242 516132 : double ohd_l=ohd.modulo();
243 : double ohd_df, ohd_sw=distanceOHSwitch( getBaseColvarNumber( myatoms.getIndex(0) ),
244 516132 : getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( ohd_l, ohd_df );
245 :
246 : Angle a;
247 516132 : Vector ood_adf, ohd_adf;
248 516132 : double angle=a.compute( ood, ohd, ood_adf, ohd_adf );
249 : double angle_df, angle_sw=angleSwitch( getBaseColvarNumber( myatoms.getIndex(0) ),
250 516132 : getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( angle, angle_df );
251 :
252 516132 : if( !doNotCalculateDerivatives() ) {
253 36 : addAtomDerivatives( 1, 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), myatoms );
254 36 : addAtomDerivatives( 1, 1, angle_sw*ohd_sw*(+ood_df)*ood + ood_sw*ohd_sw*angle_df*angle*ood_adf, myatoms );
255 36 : addAtomDerivatives( 1, iat, angle_sw*ood_sw*(+ohd_df)*ohd + ood_sw*ohd_sw*angle_df*angle*ohd_adf, myatoms );
256 72 : myatoms.addBoxDerivatives( 1, angle_sw*ohd_sw*(-ood_df)*Tensor(ood,ood) + angle_sw*ood_sw*(-ohd_df)*Tensor(ohd,ohd)
257 108 : -ood_sw*ohd_sw*angle_df*angle*(Tensor(ood,ood_adf)+Tensor(ohd,ohd_adf)) );
258 : }
259 516132 : return ood_sw*ohd_sw*angle_sw;
260 : }
261 :
262 : }
263 : }
264 :
|
