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 "multicolvar/MultiColvarBase.h"
23 : #include "HBPammObject.h"
24 : #include "tools/NeighborList.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/SwitchingFunction.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : using namespace std;
32 :
33 : namespace PLMD {
34 : namespace pamm {
35 :
36 : //+PLUMEDOC MCOLVAR HBPAMM_SH
37 : /*
38 : Number of HBPAMM hydrogen bonds formed by each hydrogen atom in the system
39 :
40 : \par Examples
41 :
42 : */
43 : //+ENDPLUMEDOC
44 :
45 :
46 : class HBPammHydrogens : public multicolvar::MultiColvarBase {
47 : private:
48 : double rcut2;
49 : unsigned block1upper,block2lower;
50 : Matrix<HBPammObject> hbpamm_obj;
51 : public:
52 : static void registerKeywords( Keywords& keys );
53 : explicit HBPammHydrogens(const ActionOptions&);
54 : // active methods:
55 : double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const override;
56 : /// Returns the number of coordinates of the field
57 0 : bool isPeriodic() override {
58 0 : return false;
59 : }
60 : };
61 :
62 13789 : PLUMED_REGISTER_ACTION(HBPammHydrogens,"HBPAMM_SH")
63 :
64 6 : void HBPammHydrogens::registerKeywords( Keywords& keys ) {
65 6 : multicolvar::MultiColvarBase::registerKeywords( keys );
66 12 : keys.add("atoms-1","HYDROGENS","The list of hydrogen atoms that can form part of a hydrogen bond. The atoms must be specified using a comma separated list.");
67 12 : keys.add("atoms-1","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 "
68 : "hydrogen bond is assumed to be the same as the set of atoms that can form hydrogen bonds. The atoms involved must be specified "
69 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
70 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
71 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
72 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
73 12 : keys.add("atoms-2","DONORS","The list of atoms which can donate a hydrogen bond. The atoms involved must be specified "
74 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
75 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
76 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
77 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
78 12 : keys.add("atoms-2","ACCEPTORS","The list of atoms which can accept a hydrogen bond. The atoms involved must be specified "
79 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
80 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
81 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
82 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
83 12 : keys.add("numbered","CLUSTERS","the name of the file that contains the definitions of all the kernels for PAMM");
84 12 : keys.add("compulsory","REGULARISE","0.001","don't allow the denominator to be smaller then this value");
85 12 : keys.reset_style("CLUSTERS","compulsory");
86 : // Use actionWithDistributionKeywords
87 6 : keys.use("MEAN");
88 6 : keys.use("MORE_THAN");
89 6 : keys.use("LESS_THAN");
90 6 : keys.use("MAX");
91 6 : keys.use("MIN");
92 6 : keys.use("BETWEEN");
93 6 : keys.use("HISTOGRAM");
94 6 : keys.use("MOMENTS");
95 6 : keys.use("ALT_MIN");
96 6 : keys.use("LOWEST");
97 6 : keys.use("HIGHEST");
98 6 : keys.use("SUM");
99 6 : }
100 :
101 2 : HBPammHydrogens::HBPammHydrogens(const ActionOptions&ao):
102 : Action(ao),
103 2 : MultiColvarBase(ao) {
104 : // Read in the atoms
105 2 : usespecies=true;
106 2 : weightHasDerivatives=false;
107 : // Read in hydrogen atom indicees
108 : std::vector<AtomNumber> all_atoms;
109 4 : parseMultiColvarAtomList("HYDROGENS",-1,all_atoms);
110 2 : if( atom_lab.size()==0 ) {
111 0 : error("no hydrogens specified in input file");
112 : }
113 : // Now create a task list - one task per hydrogen
114 2 : unsigned nH = atom_lab.size();
115 136 : for(unsigned i=0; i<nH; ++i) {
116 134 : addTaskToList(i);
117 : }
118 : // Read in other atoms in hydrogen bond
119 2 : ablocks.resize(1);
120 4 : parseMultiColvarAtomList("SITES",-1,all_atoms);
121 2 : if( atom_lab.size()>nH ) {
122 2 : block1upper=atom_lab.size() - nH + 1;
123 2 : block2lower=0;
124 2 : ablocks[0].resize( atom_lab.size() - nH );
125 69 : for(unsigned i=nH; i<atom_lab.size(); ++i) {
126 67 : ablocks[0][i-nH]=i;
127 : }
128 : } else {
129 0 : parseMultiColvarAtomList("DONORS",-1,all_atoms);
130 0 : block1upper=block2lower=atom_lab.size() - nH + 1;
131 0 : for(unsigned i=nH; i<atom_lab.size(); ++i) {
132 0 : ablocks[0].push_back(i);
133 : }
134 0 : parseMultiColvarAtomList("ACCEPTORS",-1,all_atoms);
135 0 : if( atom_lab.size()>(block1upper+nH-1) || (block1upper-1)>0 ) {
136 0 : error("no acceptor donor pairs in input specified therefore no hydrogen bonds");
137 : }
138 0 : for(unsigned i=nH+block2lower-1; i<atom_lab.size(); ++i) {
139 0 : ablocks[0].push_back(i);
140 : }
141 : }
142 2 : setupMultiColvarBase( all_atoms );
143 :
144 : double reg;
145 4 : parse("REGULARISE",reg);
146 2 : unsigned nnode_t=mybasemulticolvars.size();
147 2 : if( nnode_t==0 ) {
148 : nnode_t=1;
149 : }
150 0 : if( nnode_t==1 ) {
151 : std::string errormsg, desc;
152 4 : parse("CLUSTERS",desc);
153 : hbpamm_obj.resize(1,1);
154 2 : hbpamm_obj(0,0).setup(desc, reg, this, errormsg );
155 2 : if( errormsg.length()>0 ) {
156 0 : error( errormsg );
157 : }
158 : } else {
159 : unsigned nr=nnode_t, nc=nnode_t;
160 : hbpamm_obj.resize( nr, nc );
161 0 : for(unsigned i=0; i<nr; ++i) {
162 : // Retrieve the base number
163 : unsigned ibase;
164 0 : if( nc<10 ) {
165 0 : ibase=(i+1)*10;
166 0 : } else if ( nc<100 ) {
167 0 : ibase=(i+1)*100;
168 : } else {
169 0 : error("wow this is an error I never would have expected");
170 : }
171 :
172 0 : for(unsigned j=i; j<nc; ++j) {
173 : std::string errormsg, desc;
174 0 : parseNumbered("CLUSTERS",ibase+j+1,desc);
175 0 : if( i==j ) {
176 0 : hbpamm_obj(i,j).setup( desc, reg, this, errormsg );
177 0 : if( errormsg.length()>0 ) {
178 0 : error( errormsg );
179 : }
180 : } else {
181 0 : hbpamm_obj(i,j).setup( desc, reg, this, errormsg );
182 0 : hbpamm_obj(j,i).setup( desc, reg, this, errormsg );
183 0 : if( errormsg.length()>0 ) {
184 0 : error( errormsg );
185 : }
186 : }
187 : }
188 : }
189 : }
190 :
191 : // Set the link cell cutoff
192 2 : double sfmax=0;
193 4 : for(unsigned i=0; i<hbpamm_obj.ncols(); ++i) {
194 4 : for(unsigned j=i; j<hbpamm_obj.nrows(); ++j) {
195 2 : double rcut=hbpamm_obj(i,j).get_cutoff();
196 2 : if( rcut>sfmax ) {
197 2 : sfmax=rcut;
198 : }
199 : }
200 : }
201 2 : setLinkCellCutoff( sfmax );
202 2 : rcut2 = sfmax*sfmax;
203 :
204 : // Setup the multicolvar base
205 2 : setupMultiColvarBase( all_atoms );
206 : // And setup the ActionWithVessel
207 2 : checkRead();
208 2 : }
209 :
210 134 : double HBPammHydrogens::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
211 : double value=0, md_da ;
212 :
213 : // Calculate the coordination number
214 8344 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
215 8210 : if( i>block1upper ) {
216 0 : continue;
217 : }
218 532552 : for(unsigned j=1; j<myatoms.getNumberOfAtoms(); ++j) {
219 524342 : if( i==j || j<block2lower ) {
220 8210 : continue ;
221 : }
222 : // Get the base colvar numbers
223 516132 : unsigned dno = atom_lab[myatoms.getIndex(i)].first;
224 516132 : unsigned ano = atom_lab[myatoms.getIndex(j)].first;
225 516132 : Vector d_da=getSeparation( myatoms.getPosition(i), myatoms.getPosition(j) );
226 1032264 : if ( (md_da=d_da[0]*d_da[0])<rcut2 &&
227 516132 : (md_da+=d_da[1]*d_da[1])<rcut2 &&
228 412196 : (md_da+=d_da[2]*d_da[2])<rcut2) {
229 270884 : value += hbpamm_obj(dno,ano).evaluate( i, j, 0, d_da, sqrt(md_da), myatoms );
230 : }
231 : }
232 : }
233 :
234 134 : return value;
235 : }
236 :
237 : }
238 : }
239 :
|
