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