Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2021-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 "CoordinationBase.h"
23 : #include "tools/SwitchingFunction.h"
24 : #include "core/ActionRegister.h"
25 : #include "core/PlumedMain.h"
26 : #include "tools/IFile.h"
27 :
28 : #include <string>
29 :
30 : namespace PLMD {
31 : namespace colvar {
32 :
33 : //+PLUMEDOC COLVAR GHBFIX
34 : /*
35 : Calculate the GHBFIX interaction energy between GROUPA and GROUPB
36 :
37 : The GHBFIX interaction energy between GROUPA and GROUPB is computed
38 : using a potential that is defined in the paper cited in the bibliography below. The potential introduced in that paper is
39 : essentially a switching function that is -1 for small distances and 0 for large distances with a smooth interpolation in the middle.
40 : This function can be scaled as desired by specifying interaction scaling parameters and energy units.
41 :
42 : This collective variable can be used to analyze hydrogen bond interactions, or to generate bias potentials.
43 : Notice that the value of the GHBFIX is returned in plumed units (see [UNITS](UNITS.md)), unless the ENERGY_UNITS keyword is employed.
44 :
45 : This example prints the GHBFIX interaction in kcal/mol between two groups of atoms using D_0, D_MAX and C
46 : It is applied in the functional form introduced in the pioneering paper.
47 : The types of atoms 1-6 should be defined in typesTable_examples.dat while their interaction parameters should be defined in scalingParameters_examples.dat in kBT units.
48 :
49 : ```plumed
50 : #SETTINGS INPUTFILES=regtest/basic/rt-ghbfix/typesTable_examples.dat,regtest/basic/rt-ghbfix/scalingParameters_examples.dat
51 : gh: GHBFIX ...
52 : PAIR GROUPA=1,2,3 GROUP=4,5,6 D_0=0.2 D_MAX=0.3 C=0.8
53 : TYPES=regtest/basic/rt-ghbfix/typesTable_examples.dat
54 : PARAMS=regtest/basic/rt-ghbfix/scalingParameters_examples.dat
55 : ENERGY_UNITS=kcal/mol
56 : ...
57 : PRINT FILE=output ARG=gh
58 : ```
59 :
60 : */
61 : //+ENDPLUMEDOC
62 :
63 : class GHBFIX : public CoordinationBase {
64 :
65 : double dmax;
66 : double dmax_squared;
67 : double d0;
68 : double c;
69 :
70 : std::vector<unsigned> typesTable;
71 :
72 : std::vector<double> etas;
73 :
74 : unsigned n;
75 :
76 : double dmax2;
77 : double A;
78 : double B;
79 : double C;
80 : double D;
81 :
82 : public:
83 : explicit GHBFIX(const ActionOptions&);
84 : // active methods:
85 : static void registerKeywords( Keywords& keys );
86 : double pairing(double distance,double&dfunc,unsigned i,unsigned j)const override;
87 : };
88 :
89 : PLUMED_REGISTER_ACTION(GHBFIX,"GHBFIX")
90 :
91 9 : void GHBFIX::registerKeywords( Keywords& keys ) {
92 9 : CoordinationBase::registerKeywords(keys);
93 :
94 9 : keys.add("optional","ENERGY_UNITS","the value of ENERGY_UNITS in the switching function");
95 9 : keys.add("compulsory","TYPES","the value of TYPES in the switching function");
96 9 : keys.add("compulsory","PARAMS","the value of PARAMS in the switching function");
97 9 : keys.add("compulsory","D_MAX","the value of D_MAX in the switching function");
98 9 : keys.add("compulsory","D_0","the value of D_0 in the switching function");
99 9 : keys.add("compulsory","C","the value of C in the switching function");
100 18 : keys.setValueDescription("scalar","the GHBFIX interaction energy between the atoms in GROUPA and GROUPB");
101 9 : keys.addDOI("10.1021/acs.jctc.8b00955");
102 9 : }
103 :
104 5 : GHBFIX::GHBFIX(const ActionOptions&ao):
105 : Action(ao),
106 5 : CoordinationBase(ao) {
107 : std::string types;
108 : std::string params;
109 7 : std::string energy_units ="plumed" ;
110 :
111 5 : parse("D_MAX",dmax);
112 5 : dmax_squared = dmax*dmax;
113 5 : parse("D_0",d0);
114 5 : parse("C",c);
115 5 : parse("TYPES",types);
116 5 : parse("PARAMS",params);
117 5 : parse("ENERGY_UNITS",energy_units);
118 5 : checkRead();
119 :
120 5 : if( d0<0.0 || dmax<=d0 ) {
121 1 : error("D_MAX should be larger than D_0 and D_0 should not be negative");
122 : }
123 4 : if( !(c>0.0 && c<1.0) ) {
124 2 : error("C should be larger than 0 and smaller than 1");
125 : }
126 4 : dmax2 = dmax-d0;
127 :
128 4 : A = (-c*dmax2*dmax2)/((1-c)*dmax2*dmax2);
129 4 : B = (2*dmax2)/((1-c)*dmax2*dmax2);
130 4 : C = -1/((1-c)*dmax2*dmax2);
131 4 : D = 1/(c*dmax2*dmax2);
132 :
133 : std::map<std::string,unsigned> MapTypesTable;
134 :
135 : //typesTable
136 4 : IFile typesfile;
137 4 : typesfile.link(*this);
138 4 : typesfile.open(types);
139 : std::string itype;
140 24 : while(typesfile.scanField("itype",itype).scanField()) {
141 8 : plumed_assert(itype.empty()==false)<<"itype is empty";
142 :
143 8 : if (MapTypesTable.empty()) {
144 4 : MapTypesTable.insert({itype, 0});
145 : } else if (MapTypesTable.count(itype) == 0) {
146 : unsigned currentMax = 0;
147 8 : for(auto it = MapTypesTable.cbegin(); it != MapTypesTable.cend(); ++it ) {
148 4 : if (it ->second > currentMax) {
149 : currentMax = it->second;
150 : }
151 : }
152 8 : MapTypesTable.insert({itype, currentMax+1});
153 : }
154 :
155 8 : typesTable.push_back(MapTypesTable[itype]);
156 : }
157 :
158 4 : if( typesTable.empty() ) {
159 0 : error("no atom types were found in the file " + types);
160 : }
161 : // the atoms of the two groups should all have a type
162 11 : for(const auto & index : getAbsoluteIndexes() ) {
163 8 : if( index.index()>=typesTable.size() ) {
164 4 : error("your types table only covers " + std::to_string(typesTable.size()) + " atoms, but atom " + std::to_string(index.serial()) + " is used");
165 : }
166 : }
167 3 : n = (int)*std::max_element(std::begin(typesTable), std::end(typesTable));
168 3 : n+=1;
169 :
170 : //scalingParameters
171 4 : etas.resize(static_cast<std::size_t>(n)*n,0.0);
172 3 : IFile etafile;
173 3 : etafile.open(params);
174 : std::string it,jt;
175 : double eta;
176 : unsigned nskipped=0;
177 32 : while(etafile.scanField("itype",it).scanField("jtype",jt).scanField("eta",eta).scanField()) {
178 13 : plumed_assert(it.empty()==false)<<"itype is empty";
179 13 : plumed_assert(jt.empty()==false)<<"jtype is empty";
180 : // The file of the parameters can contain types that are not used by the atoms
181 : // in the types table: these cannot be used and must be skipped.
182 : // Using operator[] here would add them to the map as type 0 and the parameters
183 : // of the (wrong) type 0 would be overwritten.
184 : const auto iti = MapTypesTable.find(it);
185 : const auto itj = MapTypesTable.find(jt);
186 13 : if( iti==MapTypesTable.end() || itj==MapTypesTable.end() ) {
187 10 : nskipped++;
188 : continue;
189 : }
190 3 : etas[n*iti->second+itj->second]=eta;
191 : }
192 3 : if( nskipped>0 ) {
193 2 : log.printf(" skipped %u lines of %s that contain types that are not in %s\n",nskipped,params.c_str(),types.c_str());
194 : }
195 :
196 3 : if(energy_units!="plumed") {
197 3 : Units units;
198 3 : units.setEnergy(energy_units);
199 15 : for(unsigned i=0; i<etas.size(); i++) {
200 12 : etas[i]*=units.getEnergy()/getUnits().getEnergy();
201 : }
202 3 : }
203 :
204 18 : }
205 :
206 :
207 450 : double GHBFIX::pairing(double distance2,double&dfunc,unsigned i,unsigned j)const {
208 : // Most of the pairs are beyond the cutoff: check it before looking up the types
209 450 : if(distance2>dmax_squared) {
210 3 : dfunc=0.0;
211 3 : return 0.;
212 : }
213 : // that the atoms are in the table of the types is checked in the constructor
214 447 : const auto i1=getAbsoluteIndex(i).index();
215 447 : const auto i2=getAbsoluteIndex(j).index();
216 : plumed_dbg_assert(i1<typesTable.size() && i2<typesTable.size());
217 447 : const double scale=etas[n*typesTable[i1]+typesTable[i2]];
218 447 : double distance=std::sqrt(distance2);
219 447 : const double rdist = (distance-d0);
220 : double result=-1.;
221 447 : if(rdist<=0.) {
222 300 : dfunc=0.0;
223 : } else {
224 147 : dfunc=0.0;
225 147 : if (rdist > c*dmax2) {
226 27 : result+=(A + B*rdist + C*rdist*rdist);
227 27 : dfunc+=B+2*C*rdist;
228 : } else {
229 120 : result+=D*(rdist*rdist);
230 120 : dfunc+=2*D*rdist;
231 : }
232 147 : dfunc/=distance;
233 : }
234 447 : result*=scale;
235 447 : dfunc*=scale;
236 447 : return result;
237 : }
238 :
239 : }
240 : }
|
