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 : \par Examples
35 :
36 : */
37 : //+ENDPLUMEDOC
38 :
39 : //+PLUMEDOC COLVAR HBPAMM_SA
40 : /*
41 : Calculate the number of hydrogen bonds each acceptor participates in using the HBPamm method
42 :
43 : \par Examples
44 :
45 : */
46 : //+ENDPLUMEDOC
47 :
48 : //+PLUMEDOC COLVAR HBPAMM_SD
49 : /*
50 : Calculate the number of hydrogen bonds each donor participates in using the HBPamm method
51 :
52 : \par Examples
53 :
54 : */
55 : //+ENDPLUMEDOC
56 :
57 : //+PLUMEDOC COLVAR HBPAMM_SH
58 : /*
59 : Calculate the number of hydrogen bonds each hydrogen participates in using the HBPamm method
60 :
61 : \par Examples
62 :
63 : */
64 : //+ENDPLUMEDOC
65 :
66 : namespace PLMD {
67 : namespace pamm {
68 :
69 : class HBPammMatrix : public adjmat::AdjacencyMatrixBase {
70 : private:
71 : double regulariser;
72 : Tensor incoord_to_hbcoord;
73 : std::vector<double> weight;
74 : std::vector<Vector> centers;
75 : std::vector<Tensor> kmat;
76 : public:
77 : /// Create manual
78 : static void registerKeywords( Keywords& keys );
79 : /// Constructor
80 : explicit HBPammMatrix(const ActionOptions&);
81 : ///
82 : double calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const ;
83 : };
84 :
85 : PLUMED_REGISTER_ACTION(HBPammMatrix,"HBPAMM_MATRIX")
86 :
87 31 : void HBPammMatrix::registerKeywords( Keywords& keys ) {
88 31 : adjmat::AdjacencyMatrixBase::registerKeywords( keys );
89 31 : keys.use("GROUPC");
90 62 : keys.add("compulsory","ORDER","dah","the order in which the groups are specified in the input. Can be dah (donor/acceptor/hydrogens), "
91 : "adh (acceptor/donor/hydrogens) or hda (hydrogens/donor/hydrogens");
92 62 : keys.add("compulsory","CLUSTERS","the name of the file that contains the definitions of all the kernels for PAMM");
93 62 : keys.add("compulsory","REGULARISE","0.001","don't allow the denominator to be smaller then this value");
94 62 : 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))");
95 31 : keys.needsAction("PAMM");
96 31 : keys.needsAction("ONES");
97 31 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
98 31 : }
99 :
100 6 : HBPammMatrix::HBPammMatrix(const ActionOptions& ao):
101 : Action(ao),
102 6 : AdjacencyMatrixBase(ao) {
103 : double DP2CUTOFF;
104 12 : parse("GAUSS_CUTOFF",DP2CUTOFF);
105 : std::string sorder;
106 12 : parse("ORDER",sorder);
107 6 : if( sorder=="dah" ) {
108 2 : incoord_to_hbcoord(0,0)=1;
109 2 : incoord_to_hbcoord(0,1)=-1;
110 2 : incoord_to_hbcoord(0,2)=0;
111 2 : incoord_to_hbcoord(1,0)=1;
112 2 : incoord_to_hbcoord(1,1)=1;
113 2 : incoord_to_hbcoord(1,2)=0;
114 2 : incoord_to_hbcoord(2,0)=0;
115 2 : incoord_to_hbcoord(2,1)=0;
116 2 : incoord_to_hbcoord(2,2)=1;
117 2 : log.printf(" GROUPA is list of donor atoms \n");
118 4 : } else if( sorder=="adh" ) {
119 2 : incoord_to_hbcoord(0,0)=-1;
120 2 : incoord_to_hbcoord(0,1)=1;
121 2 : incoord_to_hbcoord(0,2)=0;
122 2 : incoord_to_hbcoord(1,0)=1;
123 2 : incoord_to_hbcoord(1,1)=1;
124 2 : incoord_to_hbcoord(1,2)=0;
125 2 : incoord_to_hbcoord(2,0)=0;
126 2 : incoord_to_hbcoord(2,1)=0;
127 2 : incoord_to_hbcoord(2,2)=1;
128 2 : log.printf(" GROUPA is list of acceptor atoms \n");
129 2 : } else if( sorder=="hda" ) {
130 2 : incoord_to_hbcoord(0,0)=-1;
131 2 : incoord_to_hbcoord(0,1)=0;
132 2 : incoord_to_hbcoord(0,2)=1;
133 2 : incoord_to_hbcoord(1,0)=1;
134 2 : incoord_to_hbcoord(1,1)=0;
135 2 : incoord_to_hbcoord(1,2)=1;
136 2 : incoord_to_hbcoord(2,0)=0;
137 2 : incoord_to_hbcoord(2,1)=1;
138 2 : incoord_to_hbcoord(2,2)=0;
139 2 : log.printf(" GROUPA is list of hydrogen atoms \n");
140 : } else {
141 0 : plumed_error();
142 : }
143 : // Read in the regularisation parameter
144 12 : parse("REGULARISE",regulariser);
145 :
146 : // Read in the kernels
147 : double sqr2pi = sqrt(2*pi);
148 : double sqrt2pi3 = sqr2pi*sqr2pi*sqr2pi;
149 : std::string fname;
150 6 : parse("CLUSTERS", fname);
151 6 : double sfmax=0, ww;
152 6 : Vector cent;
153 6 : Tensor covar;
154 6 : IFile ifile;
155 6 : ifile.open(fname);
156 6 : ifile.allowIgnoredFields();
157 : for(unsigned k=0;; ++k) {
158 144 : if( !ifile.scanField("height",ww) ) {
159 : break;
160 : }
161 66 : ifile.scanField("ptc",cent[0]);
162 66 : ifile.scanField("ssc",cent[1]);
163 66 : ifile.scanField("adc",cent[2]);
164 66 : ifile.scanField("sigma_ptc_ptc",covar[0][0]);
165 66 : ifile.scanField("sigma_ptc_ssc",covar[0][1]);
166 66 : ifile.scanField("sigma_ptc_adc",covar[0][2]);
167 66 : covar[1][0] = covar[0][1];
168 66 : ifile.scanField("sigma_ssc_ssc",covar[1][1]);
169 66 : ifile.scanField("sigma_ssc_adc",covar[1][2]);
170 66 : covar[2][0] = covar[0][2];
171 66 : covar[2][1] = covar[1][2];
172 66 : ifile.scanField("sigma_adc_adc",covar[2][2]);
173 66 : weight.push_back( ww / ( sqrt2pi3 * sqrt(covar.determinant()) ) );
174 66 : centers.push_back( cent );
175 66 : kmat.push_back( covar.inverse() );
176 :
177 66 : Vector eigval;
178 66 : Tensor eigvec;
179 66 : diagMatSym( covar, eigval, eigvec );
180 : unsigned ind_maxeval=0;
181 66 : double max_eval=eigval[0];
182 198 : for(unsigned i=1; i<3; ++i) {
183 132 : if( eigval[i]>max_eval ) {
184 : max_eval=eigval[i];
185 : ind_maxeval=i;
186 : }
187 : }
188 66 : double rcut = cent[2] + sqrt(2.0*DP2CUTOFF)*fabs(sqrt(max_eval)*eigvec(2,ind_maxeval));
189 66 : if( rcut > sfmax ) {
190 24 : sfmax = rcut;
191 : }
192 66 : ifile.scanField();
193 66 : }
194 6 : ifile.close();
195 6 : setLinkCellCutoff( false, sfmax );
196 12 : }
197 :
198 16286 : double HBPammMatrix::calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const {
199 16286 : Vector ddij, ddik, ddin, in_dists, hb_pamm_dists, hb_pamm_ders, real_ders;
200 16286 : ddin = pbcDistance( pos1, pos2 );
201 16286 : in_dists[2] = ddin.modulo();
202 16286 : if( in_dists[2]<epsilon ) {
203 : return 0;
204 : }
205 :
206 : double tot=0;
207 16286 : Vector disp, der, tmp_der;
208 1572892 : for(unsigned i=0; i<natoms; ++i) {
209 1556606 : ddij = getPosition(i,myvals);
210 1556606 : in_dists[0] = ddij.modulo();
211 1556606 : ddik = pbcDistance( pos2, getPosition(i,myvals) );
212 1556606 : in_dists[1] = ddik.modulo();
213 1556606 : if( in_dists[1]<epsilon ) {
214 8210 : continue;
215 : }
216 :
217 1548396 : hb_pamm_dists = matmul( incoord_to_hbcoord, in_dists );
218 1548396 : disp = hb_pamm_dists - centers[0];
219 1548396 : der = matmul( kmat[0], disp );
220 1548396 : double vv = weight[0]*exp( -dotProduct( disp, der ) / 2. );
221 1548396 : der *= -vv;
222 :
223 1548396 : double denom = regulariser + vv;
224 6193584 : for(unsigned j=0; j<3; ++j) {
225 4645188 : hb_pamm_ders[j] = der[j];
226 : }
227 17032356 : for(unsigned k=1; k<weight.size(); ++k) {
228 15483960 : disp = hb_pamm_dists - centers[k];
229 15483960 : tmp_der = matmul( kmat[k], disp );
230 15483960 : double tval = weight[k]*exp( -dotProduct( disp, tmp_der ) / 2. );
231 15483960 : denom += tval;
232 15483960 : hb_pamm_ders += -tmp_der*tval;
233 : }
234 1548396 : double vf = vv / denom;
235 1548396 : tot += vf;
236 1548396 : if( fabs(vf)<epsilon ) {
237 1546155 : continue;
238 : }
239 : // Now get derivatives
240 2241 : real_ders = matmul( der / denom - vf*hb_pamm_ders/denom, incoord_to_hbcoord );
241 :
242 : // And add the derivatives to the underlying atoms
243 2241 : addAtomDerivatives( 0, -(real_ders[0]/in_dists[0])*ddij - (real_ders[2]/in_dists[2])*ddin, myvals );
244 2241 : addAtomDerivatives( 1, -(real_ders[1]/in_dists[1])*ddik + (real_ders[2]/in_dists[2])*ddin, myvals );
245 2241 : addThirdAtomDerivatives( i, (real_ders[0]/in_dists[0])*ddij + (real_ders[1]/in_dists[1])*ddik, myvals );
246 4482 : addBoxDerivatives( -(real_ders[0]/in_dists[0])*Tensor( ddij, ddij )
247 4482 : -(real_ders[1]/in_dists[1])*Tensor( ddik, ddik )
248 6723 : -(real_ders[2]/in_dists[2])*Tensor( ddin, ddin ), myvals );
249 : }
250 16286 : return tot;
251 : }
252 :
253 : class HBPammShortcut : public ActionShortcut {
254 : public:
255 : static void registerKeywords( Keywords& keys );
256 : HBPammShortcut(const ActionOptions&);
257 : };
258 :
259 : PLUMED_REGISTER_ACTION(HBPammShortcut,"HBPAMM_SD")
260 : PLUMED_REGISTER_ACTION(HBPammShortcut,"HBPAMM_SA")
261 : PLUMED_REGISTER_ACTION(HBPammShortcut,"HBPAMM_SH")
262 :
263 17 : void HBPammShortcut::registerKeywords( Keywords& keys ) {
264 17 : HBPammMatrix::registerKeywords( keys );
265 17 : keys.remove("GROUP");
266 17 : keys.remove("GROUPA");
267 17 : keys.remove("GROUPB");
268 17 : keys.remove("COMPONENTS");
269 34 : keys.add("optional","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 "
270 : "hydrogen bond is assumed to be the same as the set of atoms that can form hydrogen bonds. The atoms involved must be specified"
271 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
272 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
273 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
274 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
275 34 : keys.add("optional","DONORS","The list of atoms which can donate a hydrogen bond. The atoms involved must be specified "
276 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
277 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
278 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
279 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
280 34 : keys.add("optional","ACCEPTORS","The list of atoms which can accept a hydrogen bond. The atoms involved must be specified "
281 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
282 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
283 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
284 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
285 34 : keys.add("compulsory","HYDROGENS","The list of hydrogen atoms that can form part of a hydrogen bond. The atoms must be specified using a comma separated list, "
286 : "an index range or by using a \\ref GROUP");
287 17 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
288 17 : keys.needsAction("HBPAMM_MATRIX");
289 17 : }
290 :
291 6 : HBPammShortcut::HBPammShortcut(const ActionOptions&ao):
292 : Action(ao),
293 6 : ActionShortcut(ao) {
294 6 : std::string mwords = getShortcutLabel() + "_mat: HBPAMM_MATRIX";
295 6 : if( getName()=="HBPAMM_SD" ) {
296 : std::string site_str;
297 4 : parse("SITES",site_str);
298 2 : if( site_str.length()>0 ) {
299 4 : mwords += " GROUP=" + site_str;
300 : } else {
301 : std::string d_str;
302 0 : parse("DONORS",d_str);
303 0 : mwords += " GROUPA=" + d_str;
304 : std::string a_str;
305 0 : parse("ACCEPTORS",a_str);
306 0 : mwords += " GROUPB=" + a_str;
307 : }
308 : std::string h_str;
309 2 : parse("HYDROGENS",h_str);
310 4 : mwords += " GROUPC=" + h_str + " ORDER=dah";
311 4 : } else if( getName()=="HBPAMM_SA" ) {
312 : std::string site_str;
313 4 : parse("SITES",site_str);
314 2 : if( site_str.length()>0 ) {
315 4 : mwords += " GROUP=" + site_str;
316 : } else {
317 : std::string a_str;
318 0 : parse("ACCEPTORS",a_str);
319 0 : mwords += " GROUPA=" + a_str;
320 : std::string d_str;
321 0 : parse("DONORS",d_str);
322 0 : mwords += " GROUPB=" + d_str;
323 : }
324 : std::string h_str;
325 2 : parse("HYDROGENS",h_str);
326 4 : mwords += " GROUPC=" + h_str + " ORDER=adh";
327 2 : } else if( getName()=="HBPAMM_SH" ) {
328 : std::string h_str;
329 2 : parse("HYDROGENS",h_str);
330 4 : mwords += " GROUPA=" + h_str + " ORDER=hda";
331 : std::string site_str;
332 4 : parse("SITES",site_str);
333 2 : if( site_str.length()>0 ) {
334 2 : mwords += " GROUPB=" + site_str;
335 4 : mwords += " GROUPC=" + site_str;
336 : } else {
337 : std::string d_str;
338 0 : parse("DONORS",d_str);
339 0 : mwords += " GROUPB=" + d_str;
340 : std::string a_str;
341 0 : parse("ACCEPTORS",a_str);
342 0 : mwords += " GROUPC=" + a_str;
343 : }
344 : }
345 : std::map<std::string,std::string> keymap;
346 6 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
347 12 : readInputLine( mwords + " " + convertInputLineToString() );
348 6 : ActionWithValue* mb=plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
349 6 : plumed_assert( mb );
350 : std::string nsize;
351 6 : Tools::convert( (mb->copyOutput(getShortcutLabel() + "_mat"))->getShape()[1], nsize );
352 12 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + nsize );
353 12 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat," + getShortcutLabel() + "_ones");
354 12 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
355 6 : }
356 :
357 : }
358 : }
|
