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 ); keys.use("GROUPC");
89 62 : keys.add("compulsory","ORDER","dah","the order in which the groups are specified in the input. Can be dah (donor/acceptor/hydrogens), "
90 : "adh (acceptor/donor/hydrogens) or hda (hydrogens/donor/hydrogens");
91 62 : keys.add("compulsory","CLUSTERS","the name of the file that contains the definitions of all the kernels for PAMM");
92 62 : keys.add("compulsory","REGULARISE","0.001","don't allow the denominator to be smaller then this value");
93 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))");
94 93 : keys.needsAction("PAMM"); keys.needsAction("ONES"); keys.needsAction("MATRIX_VECTOR_PRODUCT");
95 31 : }
96 :
97 6 : HBPammMatrix::HBPammMatrix(const ActionOptions& ao):
98 : Action(ao),
99 6 : AdjacencyMatrixBase(ao)
100 : {
101 18 : double DP2CUTOFF; parse("GAUSS_CUTOFF",DP2CUTOFF); std::string sorder; parse("ORDER",sorder);
102 6 : if( sorder=="dah" ) {
103 2 : incoord_to_hbcoord(0,0)=1; incoord_to_hbcoord(0,1)=-1; incoord_to_hbcoord(0,2)=0;
104 2 : incoord_to_hbcoord(1,0)=1; incoord_to_hbcoord(1,1)=1; incoord_to_hbcoord(1,2)=0;
105 2 : incoord_to_hbcoord(2,0)=0; incoord_to_hbcoord(2,1)=0; incoord_to_hbcoord(2,2)=1;
106 2 : log.printf(" GROUPA is list of donor atoms \n");
107 4 : } else if( sorder=="adh" ) {
108 2 : incoord_to_hbcoord(0,0)=-1; incoord_to_hbcoord(0,1)=1; incoord_to_hbcoord(0,2)=0;
109 2 : incoord_to_hbcoord(1,0)=1; incoord_to_hbcoord(1,1)=1; incoord_to_hbcoord(1,2)=0;
110 2 : incoord_to_hbcoord(2,0)=0; incoord_to_hbcoord(2,1)=0; incoord_to_hbcoord(2,2)=1;
111 2 : log.printf(" GROUPA is list of acceptor atoms \n");
112 2 : } else if( sorder=="hda" ) {
113 2 : incoord_to_hbcoord(0,0)=-1; incoord_to_hbcoord(0,1)=0; incoord_to_hbcoord(0,2)=1;
114 2 : incoord_to_hbcoord(1,0)=1; incoord_to_hbcoord(1,1)=0; incoord_to_hbcoord(1,2)=1;
115 2 : incoord_to_hbcoord(2,0)=0; incoord_to_hbcoord(2,1)=1; incoord_to_hbcoord(2,2)=0;
116 2 : log.printf(" GROUPA is list of hydrogen atoms \n");
117 0 : } else plumed_error();
118 : // Read in the regularisation parameter
119 12 : parse("REGULARISE",regulariser);
120 :
121 : // Read in the kernels
122 : double sqr2pi = sqrt(2*pi); double sqrt2pi3 = sqr2pi*sqr2pi*sqr2pi;
123 12 : std::string fname; parse("CLUSTERS", fname); double sfmax=0, ww; Vector cent; Tensor covar;
124 6 : IFile ifile; ifile.open(fname); ifile.allowIgnoredFields();
125 : for(unsigned k=0;; ++k) {
126 144 : if( !ifile.scanField("height",ww) ) break;
127 198 : ifile.scanField("ptc",cent[0]); ifile.scanField("ssc",cent[1]); ifile.scanField("adc",cent[2]);
128 198 : ifile.scanField("sigma_ptc_ptc",covar[0][0]); ifile.scanField("sigma_ptc_ssc",covar[0][1]); ifile.scanField("sigma_ptc_adc",covar[0][2]);
129 132 : covar[1][0] = covar[0][1]; ifile.scanField("sigma_ssc_ssc",covar[1][1]); ifile.scanField("sigma_ssc_adc",covar[1][2]);
130 66 : covar[2][0] = covar[0][2]; covar[2][1] = covar[1][2]; ifile.scanField("sigma_adc_adc",covar[2][2]);
131 66 : weight.push_back( ww / ( sqrt2pi3 * sqrt(covar.determinant()) ) );
132 66 : centers.push_back( cent ); kmat.push_back( covar.inverse() );
133 :
134 66 : Vector eigval; Tensor eigvec; diagMatSym( covar, eigval, eigvec );
135 66 : unsigned ind_maxeval=0; double max_eval=eigval[0];
136 198 : for(unsigned i=1; i<3; ++i) {
137 132 : if( eigval[i]>max_eval ) { max_eval=eigval[i]; ind_maxeval=i; }
138 : }
139 66 : double rcut = cent[2] + sqrt(2.0*DP2CUTOFF)*fabs(sqrt(max_eval)*eigvec(2,ind_maxeval));
140 66 : if( rcut > sfmax ) sfmax = rcut;
141 66 : ifile.scanField();
142 66 : }
143 6 : ifile.close(); setLinkCellCutoff( false, sfmax );
144 12 : }
145 :
146 16286 : double HBPammMatrix::calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const {
147 16286 : Vector ddij, ddik, ddin, in_dists, hb_pamm_dists, hb_pamm_ders, real_ders;
148 16286 : ddin = pbcDistance( pos1, pos2 ); in_dists[2] = ddin.modulo();
149 16286 : if( in_dists[2]<epsilon ) return 0;
150 :
151 16286 : double tot=0; Vector disp, der, tmp_der;
152 1572892 : for(unsigned i=0; i<natoms; ++i) {
153 1556606 : ddij = getPosition(i,myvals); in_dists[0] = ddij.modulo();
154 1556606 : ddik = pbcDistance( pos2, getPosition(i,myvals) ); in_dists[1] = ddik.modulo();
155 1556606 : if( in_dists[1]<epsilon ) continue;
156 :
157 1548396 : hb_pamm_dists = matmul( incoord_to_hbcoord, in_dists );
158 1548396 : disp = hb_pamm_dists - centers[0]; der = matmul( kmat[0], disp );
159 1548396 : double vv = weight[0]*exp( -dotProduct( disp, der ) / 2. ); der *= -vv;
160 :
161 6193584 : double denom = regulariser + vv; for(unsigned j=0; j<3; ++j) hb_pamm_ders[j] = der[j];
162 17032356 : for(unsigned k=1; k<weight.size(); ++k) {
163 15483960 : disp = hb_pamm_dists - centers[k]; tmp_der = matmul( kmat[k], disp );
164 15483960 : double tval = weight[k]*exp( -dotProduct( disp, tmp_der ) / 2. );
165 15483960 : denom += tval; hb_pamm_ders += -tmp_der*tval;
166 : }
167 1548396 : double vf = vv / denom; tot += vf;
168 1548396 : if( fabs(vf)<epsilon ) continue;
169 : // Now get derivatives
170 2241 : real_ders = matmul( der / denom - vf*hb_pamm_ders/denom, incoord_to_hbcoord );
171 :
172 : // And add the derivatives to the underlying atoms
173 2241 : addAtomDerivatives( 0, -(real_ders[0]/in_dists[0])*ddij - (real_ders[2]/in_dists[2])*ddin, myvals );
174 2241 : addAtomDerivatives( 1, -(real_ders[1]/in_dists[1])*ddik + (real_ders[2]/in_dists[2])*ddin, myvals );
175 2241 : addThirdAtomDerivatives( i, (real_ders[0]/in_dists[0])*ddij + (real_ders[1]/in_dists[1])*ddik, myvals );
176 4482 : addBoxDerivatives( -(real_ders[0]/in_dists[0])*Tensor( ddij, ddij )
177 4482 : -(real_ders[1]/in_dists[1])*Tensor( ddik, ddik )
178 6723 : -(real_ders[2]/in_dists[2])*Tensor( ddin, ddin ), myvals );
179 : }
180 16286 : return tot;
181 : }
182 :
183 : class HBPammShortcut : public ActionShortcut {
184 : public:
185 : static void registerKeywords( Keywords& keys );
186 : HBPammShortcut(const ActionOptions&);
187 : };
188 :
189 : PLUMED_REGISTER_ACTION(HBPammShortcut,"HBPAMM_SD")
190 : PLUMED_REGISTER_ACTION(HBPammShortcut,"HBPAMM_SA")
191 : PLUMED_REGISTER_ACTION(HBPammShortcut,"HBPAMM_SH")
192 :
193 17 : void HBPammShortcut::registerKeywords( Keywords& keys ) {
194 68 : HBPammMatrix::registerKeywords( keys ); keys.remove("GROUP"); keys.remove("GROUPA"); keys.remove("GROUPB"); keys.remove("COMPONENTS");
195 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 "
196 : "hydrogen bond is assumed to be the same as the set of atoms that can form hydrogen bonds. The atoms involved must be specified"
197 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
198 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
199 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
200 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
201 34 : keys.add("optional","DONORS","The list of atoms which can donate a hydrogen bond. The atoms involved must be specified "
202 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
203 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
204 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
205 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
206 34 : keys.add("optional","ACCEPTORS","The list of atoms which can accept a hydrogen bond. The atoms involved must be specified "
207 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
208 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
209 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
210 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
211 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, "
212 : "an index range or by using a \\ref GROUP");
213 17 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys ); keys.needsAction("HBPAMM_MATRIX");
214 17 : }
215 :
216 6 : HBPammShortcut::HBPammShortcut(const ActionOptions&ao):
217 : Action(ao),
218 6 : ActionShortcut(ao)
219 : {
220 6 : std::string mwords = getShortcutLabel() + "_mat: HBPAMM_MATRIX";
221 6 : if( getName()=="HBPAMM_SD" ) {
222 4 : std::string site_str; parse("SITES",site_str);
223 4 : if( site_str.length()>0 ) mwords += " GROUP=" + site_str;
224 : else {
225 0 : std::string d_str; parse("DONORS",d_str); mwords += " GROUPA=" + d_str;
226 0 : std::string a_str; parse("ACCEPTORS",a_str); mwords += " GROUPB=" + a_str;
227 : }
228 6 : std::string h_str; parse("HYDROGENS",h_str); mwords += " GROUPC=" + h_str + " ORDER=dah";
229 4 : } else if( getName()=="HBPAMM_SA" ) {
230 4 : std::string site_str; parse("SITES",site_str);
231 4 : if( site_str.length()>0 ) mwords += " GROUP=" + site_str;
232 : else {
233 0 : std::string a_str; parse("ACCEPTORS",a_str); mwords += " GROUPA=" + a_str;
234 0 : std::string d_str; parse("DONORS",d_str); mwords += " GROUPB=" + d_str;
235 : }
236 6 : std::string h_str; parse("HYDROGENS",h_str); mwords += " GROUPC=" + h_str + " ORDER=adh";
237 2 : } else if( getName()=="HBPAMM_SH" ) {
238 6 : std::string h_str; parse("HYDROGENS",h_str); mwords += " GROUPA=" + h_str + " ORDER=hda";
239 4 : std::string site_str; parse("SITES",site_str);
240 2 : if( site_str.length()>0 ) {
241 2 : mwords += " GROUPB=" + site_str;
242 4 : mwords += " GROUPC=" + site_str;
243 : } else {
244 0 : std::string d_str; parse("DONORS",d_str); mwords += " GROUPB=" + d_str;
245 0 : std::string a_str; parse("ACCEPTORS",a_str); mwords += " GROUPC=" + a_str;
246 : }
247 : }
248 6 : std::map<std::string,std::string> keymap; multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
249 12 : readInputLine( mwords + " " + convertInputLineToString() );
250 6 : ActionWithValue* mb=plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
251 6 : plumed_assert( mb ); std::string nsize; Tools::convert( (mb->copyOutput(getShortcutLabel() + "_mat"))->getShape()[1], nsize );
252 12 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + nsize );
253 12 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat," + getShortcutLabel() + "_ones");
254 12 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
255 6 : }
256 :
257 : }
258 : }
|
