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