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 { return false; }
58 : };
59 :
60 10423 : PLUMED_REGISTER_ACTION(HBPammHydrogens,"HBPAMM_SH")
61 :
62 3 : void HBPammHydrogens::registerKeywords( Keywords& keys ) {
63 3 : multicolvar::MultiColvarBase::registerKeywords( keys );
64 6 : 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.");
65 6 : 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 "
66 : "hydrogen bond is assumed to be the same as the set of atoms that can form hydrogen bonds. The atoms involved must be specified "
67 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
68 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
69 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
70 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
71 6 : keys.add("atoms-2","DONORS","The list of atoms which can donate a hydrogen bond. The atoms involved must be specified "
72 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
73 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
74 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
75 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
76 6 : keys.add("atoms-2","ACCEPTORS","The list of atoms which can accept a hydrogen bond. The atoms involved must be specified "
77 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
78 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
79 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
80 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
81 6 : keys.add("numbered","CLUSTERS","the name of the file that contains the definitions of all the kernels for PAMM");
82 6 : keys.add("compulsory","REGULARISE","0.001","don't allow the denominator to be smaller then this value");
83 6 : keys.reset_style("CLUSTERS","compulsory");
84 : // Use actionWithDistributionKeywords
85 12 : keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN"); keys.use("MAX");
86 12 : keys.use("MIN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
87 12 : keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST"); keys.use("SUM");
88 3 : }
89 :
90 2 : HBPammHydrogens::HBPammHydrogens(const ActionOptions&ao):
91 : Action(ao),
92 2 : MultiColvarBase(ao)
93 : {
94 : // Read in the atoms
95 2 : usespecies=true; weightHasDerivatives=false;
96 : // Read in hydrogen atom indicees
97 4 : std::vector<AtomNumber> all_atoms; parseMultiColvarAtomList("HYDROGENS",-1,all_atoms);
98 2 : if( atom_lab.size()==0 ) error("no hydrogens specified in input file");
99 : // Now create a task list - one task per hydrogen
100 136 : unsigned nH = atom_lab.size(); for(unsigned i=0; i<nH; ++i) addTaskToList(i);
101 : // Read in other atoms in hydrogen bond
102 4 : ablocks.resize(1); parseMultiColvarAtomList("SITES",-1,all_atoms);
103 2 : if( atom_lab.size()>nH ) {
104 2 : block1upper=atom_lab.size() - nH + 1; block2lower=0; ablocks[0].resize( atom_lab.size() - nH );
105 69 : for(unsigned i=nH; i<atom_lab.size(); ++i) ablocks[0][i-nH]=i;
106 : } else {
107 0 : parseMultiColvarAtomList("DONORS",-1,all_atoms);
108 0 : block1upper=block2lower=atom_lab.size() - nH + 1;
109 0 : for(unsigned i=nH; i<atom_lab.size(); ++i) ablocks[0].push_back(i);
110 0 : parseMultiColvarAtomList("ACCEPTORS",-1,all_atoms);
111 0 : if( atom_lab.size()>(block1upper+nH-1) || (block1upper-1)>0 ) error("no acceptor donor pairs in input specified therefore no hydrogen bonds");
112 0 : for(unsigned i=nH+block2lower-1; i<atom_lab.size(); ++i) ablocks[0].push_back(i);
113 : }
114 2 : setupMultiColvarBase( all_atoms );
115 :
116 4 : double reg; parse("REGULARISE",reg);
117 2 : unsigned nnode_t=mybasemulticolvars.size(); if( nnode_t==0 ) nnode_t=1;
118 0 : if( nnode_t==1 ) {
119 4 : std::string errormsg, desc; parse("CLUSTERS",desc);
120 : hbpamm_obj.resize(1,1);
121 2 : hbpamm_obj(0,0).setup(desc, reg, this, errormsg );
122 2 : if( errormsg.length()>0 ) error( errormsg );
123 : } else {
124 : unsigned nr=nnode_t, nc=nnode_t;
125 : hbpamm_obj.resize( nr, nc );
126 0 : for(unsigned i=0; i<nr; ++i) {
127 : // Retrieve the base number
128 : unsigned ibase;
129 0 : if( nc<10 ) {
130 0 : ibase=(i+1)*10;
131 0 : } else if ( nc<100 ) {
132 0 : ibase=(i+1)*100;
133 : } else {
134 0 : error("wow this is an error I never would have expected");
135 : }
136 :
137 0 : for(unsigned j=i; j<nc; ++j) {
138 0 : std::string errormsg, desc; parseNumbered("CLUSTERS",ibase+j+1,desc);
139 0 : if( i==j ) {
140 0 : hbpamm_obj(i,j).setup( desc, reg, this, errormsg );
141 0 : if( errormsg.length()>0 ) error( errormsg );
142 : } else {
143 0 : hbpamm_obj(i,j).setup( desc, reg, this, errormsg );
144 0 : hbpamm_obj(j,i).setup( desc, reg, this, errormsg );
145 0 : if( errormsg.length()>0 ) error( errormsg );
146 : }
147 : }
148 : }
149 : }
150 :
151 : // Set the link cell cutoff
152 2 : double sfmax=0;
153 4 : for(unsigned i=0; i<hbpamm_obj.ncols(); ++i) {
154 4 : for(unsigned j=i; j<hbpamm_obj.nrows(); ++j) {
155 2 : double rcut=hbpamm_obj(i,j).get_cutoff();
156 2 : if( rcut>sfmax ) { sfmax=rcut; }
157 : }
158 : }
159 2 : setLinkCellCutoff( sfmax );
160 2 : rcut2 = sfmax*sfmax;
161 :
162 : // Setup the multicolvar base
163 2 : setupMultiColvarBase( all_atoms );
164 : // And setup the ActionWithVessel
165 2 : checkRead();
166 2 : }
167 :
168 134 : double HBPammHydrogens::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
169 : double value=0, md_da ;
170 :
171 : // Calculate the coordination number
172 8344 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
173 8210 : if( i>block1upper ) continue;
174 532552 : for(unsigned j=1; j<myatoms.getNumberOfAtoms(); ++j) {
175 524342 : if( i==j || j<block2lower ) continue ;
176 : // Get the base colvar numbers
177 516132 : unsigned dno = atom_lab[myatoms.getIndex(i)].first;
178 516132 : unsigned ano = atom_lab[myatoms.getIndex(j)].first;
179 516132 : Vector d_da=getSeparation( myatoms.getPosition(i), myatoms.getPosition(j) );
180 1032264 : if ( (md_da=d_da[0]*d_da[0])<rcut2 &&
181 516132 : (md_da+=d_da[1]*d_da[1])<rcut2 &&
182 412196 : (md_da+=d_da[2]*d_da[2])<rcut2) value += hbpamm_obj(dno,ano).evaluate( i, j, 0, d_da, sqrt(md_da), myatoms );
183 : }
184 : }
185 :
186 134 : return value;
187 : }
188 :
189 : }
190 : }
191 :
|
