Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2021 Omar Valsson
3 :
4 : This file is part of S2 contact model module
5 :
6 : The S2 contact model module is free software: you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation, either version 3 of the License, or
9 : (at your option) any later version.
10 :
11 : The S2 contact model module is distributed in the hope that it will be useful,
12 : but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : GNU Lesser General Public License for more details.
15 :
16 : You should have received a copy of the GNU Lesser General Public License
17 : along with the S2 contact model module. If not, see <http://www.gnu.org/licenses/>.
18 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
19 :
20 : #include "colvar/Colvar.h"
21 : #include "tools/NeighborList.h"
22 : #include "tools/Communicator.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 :
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : namespace PLMD {
31 : namespace s2cm {
32 :
33 :
34 : //
35 :
36 : //+PLUMEDOC COLVAR S2CM
37 : /*
38 : S2 contact model CV.
39 :
40 : This CV was used in the first paper cited below. It is based on NH order parameter from the second paper cited below and the methyl order parameter from the third. Input parameters can be found in the relevant papers.
41 :
42 :
43 : */
44 : //+ENDPLUMEDOC
45 :
46 :
47 :
48 :
49 :
50 :
51 :
52 : class S2ContactModel : public Colvar {
53 : bool pbc_;
54 : bool serial_;
55 : std::unique_ptr<NeighborList> nl;
56 : bool invalidateList;
57 : bool firsttime;
58 : //
59 : double r_eff_;
60 : double inv_r_eff_;
61 : double prefactor_a_;
62 : double exp_b_;
63 : double offset_c_;
64 : double n_i_;
65 : double total_prefactor_;
66 : double r_globalshift_;
67 :
68 : enum ModelType {methyl,nh} modeltype_;
69 :
70 : //
71 : public:
72 : explicit S2ContactModel(const ActionOptions&);
73 : ~S2ContactModel();
74 : virtual void calculate();
75 : virtual void prepare();
76 : static void registerKeywords( Keywords& keys );
77 : };
78 :
79 : PLUMED_REGISTER_ACTION(S2ContactModel,"S2CM")
80 :
81 28 : void S2ContactModel::registerKeywords( Keywords& keys ) {
82 28 : Colvar::registerKeywords(keys);
83 28 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
84 28 : keys.addFlag("NLIST",false,"Use a neighbour list to speed up the calculation");
85 28 : keys.add("optional","NL_CUTOFF","The cutoff for the neighbour list");
86 28 : keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbour list");
87 28 : keys.add("atoms","METHYL_ATOM","the methyl carbon atom of the residue (i)");
88 28 : keys.add("atoms","NH_ATOMS","the hydrogen atom of the NH group of the residue (i) and carbonyl oxygen of the preceding residue (i-1)");
89 28 : keys.add("atoms","HEAVY_ATOMS","the heavy atoms to be included in the calculation");
90 : //
91 28 : keys.add("compulsory","R_EFF","the effective distance, r_eff in the equation, given in nm.");
92 28 : keys.add("compulsory","PREFACTOR_A","the prefactor, a in the equation");
93 28 : keys.add("compulsory","EXPONENT_B","the exponent, b in the equation");
94 28 : keys.add("compulsory","OFFSET_C","the offset, c in the equation");
95 28 : keys.add("compulsory","N_I"," n_i in the equation");
96 28 : keys.add("optional","R_SHIFT","shift all distances by given amount");
97 56 : keys.setValueDescription("scalar","the value of the CV");
98 28 : keys.addDOI("10.1021/acs.jpclett.7b01770");
99 28 : keys.addDOI("10.1021/ja027847a");
100 28 : keys.addDOI("10.1023/B:JNMR.0000032612.70767.35");
101 28 : }
102 :
103 24 : S2ContactModel::S2ContactModel(const ActionOptions&ao):
104 : PLUMED_COLVAR_INIT(ao),
105 24 : pbc_(true),
106 24 : serial_(false),
107 24 : invalidateList(true),
108 24 : firsttime(true),
109 24 : r_eff_(0.0),
110 24 : inv_r_eff_(0.0),
111 24 : prefactor_a_(0.0),
112 24 : exp_b_(0.0),
113 24 : offset_c_(0.0),
114 24 : n_i_(0.0),
115 24 : total_prefactor_(0.0),
116 24 : r_globalshift_(0.0),
117 24 : modeltype_(methyl) {
118 :
119 24 : parseFlag("SERIAL",serial_);
120 :
121 24 : std::vector<AtomNumber> methyl_atom;
122 24 : parseAtomList("METHYL_ATOM",methyl_atom);
123 24 : std::vector<AtomNumber> nh_atoms;
124 48 : parseAtomList("NH_ATOMS",nh_atoms);
125 :
126 24 : if(methyl_atom.size()==0 && nh_atoms.size()==0) {
127 0 : error("you have to give either METHYL_ATOM or NH_ATOMS");
128 : }
129 24 : if(methyl_atom.size()>0 && nh_atoms.size()>0) {
130 0 : error("you cannot give both METHYL_ATOM or NH_ATOMS");
131 : }
132 24 : if(methyl_atom.size()>0 && methyl_atom.size()!=1) {
133 0 : error("you should give one atom in METHYL_ATOM, the methyl carbon atom of the residue");
134 : }
135 24 : if(nh_atoms.size()>0 && nh_atoms.size()!=2) {
136 0 : error("you should give two atoms in NH_ATOMS, the hydrogen atom of the NH group of the residue (i) and carbonyl oxygen of the preceding residue (i-1)");
137 : }
138 :
139 24 : std::vector<AtomNumber> heavy_atoms;
140 48 : parseAtomList("HEAVY_ATOMS",heavy_atoms);
141 24 : if(heavy_atoms.size()==0) {
142 0 : error("HEAVY_ATOMS is not given");
143 : }
144 :
145 24 : std::vector<AtomNumber> main_atoms;
146 24 : if(methyl_atom.size()>0) {
147 0 : modeltype_= methyl;
148 0 : main_atoms = methyl_atom;
149 24 : } else if(nh_atoms.size()>0) {
150 24 : modeltype_= nh;
151 24 : main_atoms = nh_atoms;
152 : }
153 :
154 24 : bool nopbc=!pbc_;
155 24 : parseFlag("NOPBC",nopbc);
156 24 : pbc_=!nopbc;
157 :
158 : // neighbor list stuff
159 24 : bool doneigh=false;
160 24 : double nl_cut=0.0;
161 24 : int nl_st=0;
162 24 : parseFlag("NLIST",doneigh);
163 24 : if(doneigh) {
164 12 : parse("NL_CUTOFF",nl_cut);
165 12 : if(nl_cut<=0.0) {
166 0 : error("NL_CUTOFF should be explicitly specified and positive");
167 : }
168 12 : parse("NL_STRIDE",nl_st);
169 12 : if(nl_st<=0) {
170 0 : error("NL_STRIDE should be explicitly specified and positive");
171 : }
172 : }
173 :
174 24 : parse("R_EFF",r_eff_);
175 24 : inv_r_eff_ = 1.0/r_eff_;
176 24 : parse("PREFACTOR_A",prefactor_a_);
177 24 : parse("EXPONENT_B",exp_b_);
178 24 : parse("OFFSET_C",offset_c_);
179 : unsigned int n_i_int;
180 24 : parse("N_I",n_i_int);
181 24 : n_i_ = static_cast<double>(n_i_int);
182 24 : total_prefactor_ = prefactor_a_/pow(n_i_,exp_b_);
183 : //
184 24 : parse("R_SHIFT",r_globalshift_);
185 :
186 24 : checkRead();
187 :
188 24 : addValueWithDerivatives();
189 24 : setNotPeriodic();
190 :
191 24 : bool dopair=false;
192 24 : if(doneigh) {
193 24 : nl=Tools::make_unique<NeighborList>(main_atoms,heavy_atoms,serial_,dopair,pbc_,getPbc(),comm,nl_cut,nl_st);
194 : } else {
195 24 : nl=Tools::make_unique<NeighborList>(main_atoms,heavy_atoms,serial_,dopair,pbc_,getPbc(),comm);
196 : }
197 :
198 24 : requestAtoms(nl->getFullAtomList());
199 :
200 :
201 24 : log.printf(" NMR S2 contact model CV, please read and cite ");
202 24 : log << plumed.cite("Palazzesi, Valsson, and Parrinello, J. Phys. Chem. Lett. 8, 4752 (2017) - DOI:10.1021/acs.jpclett.7b01770");
203 :
204 24 : if(modeltype_==methyl) {
205 0 : log << plumed.cite("Ming and Bruschweiler, J. Biomol. NMR, 29, 363 (2004) - DOI:10.1023/B:JNMR.0000032612.70767.35");
206 0 : log.printf("\n");
207 0 : log.printf(" calculation of methyl order parameter using atom %d \n",methyl_atom[0].serial());
208 24 : } else if(modeltype_==nh) {
209 24 : log << plumed.cite("Zhang and Bruschweiler, J. Am. Chem. Soc. 124, 12654 (2002) - DOI:10.1021/ja027847a");
210 24 : log.printf("\n");
211 24 : log.printf(" calculation of NH order parameter using atoms %d and %d\n",nh_atoms[0].serial(),nh_atoms[1].serial());
212 : }
213 24 : log.printf(" heavy atoms used in the calculation (%u in total):\n",static_cast<unsigned int>(heavy_atoms.size()));
214 1872 : for(unsigned int i=0; i<heavy_atoms.size(); ++i) {
215 1848 : if( (i+1) % 25 == 0 ) {
216 72 : log.printf(" \n");
217 : }
218 1848 : log.printf(" %d", heavy_atoms[i].serial());
219 : }
220 24 : log.printf(" \n");
221 24 : log.printf(" total number of distances: %u\n",nl->size());
222 : //
223 24 : log.printf(" using parameters");
224 24 : log.printf(" r_eff=%f,",r_eff_);
225 24 : log.printf(" a=%f,",prefactor_a_);
226 24 : log.printf(" b=%f,",exp_b_);
227 24 : log.printf(" c=%f,",offset_c_);
228 24 : log.printf(" n_i=%u",n_i_int);
229 24 : if(r_globalshift_!=0.0) {
230 4 : log.printf(", r_shift=%f",r_globalshift_);
231 : }
232 24 : log.printf("\n");
233 24 : if(pbc_) {
234 12 : log.printf(" using periodic boundary conditions\n");
235 : } else {
236 12 : log.printf(" without periodic boundary conditions\n");
237 : }
238 24 : if(doneigh) {
239 12 : log.printf(" using neighbor lists with\n");
240 12 : log.printf(" update every %d steps and cutoff %f\n",nl_st,nl_cut);
241 : }
242 24 : if(serial_) {
243 0 : log.printf(" calculation done in serial\n");
244 : }
245 24 : }
246 :
247 48 : S2ContactModel::~S2ContactModel() {
248 48 : }
249 :
250 48 : void S2ContactModel::prepare() {
251 48 : std::tie(firsttime,invalidateList) = nl->prepare(this,firsttime,invalidateList).get();
252 48 : }
253 :
254 : // calculator
255 5952 : void S2ContactModel::calculate() {
256 :
257 5952 : Tensor virial;
258 5952 : std::vector<Vector> deriv(getNumberOfAtoms());
259 :
260 5952 : if(nl->getStride()>0 && invalidateList) {
261 2976 : nl->update(getPositions());
262 : }
263 :
264 5952 : unsigned int stride=comm.Get_size();
265 5952 : unsigned int rank=comm.Get_rank();
266 5952 : if(serial_) {
267 : stride=1;
268 : rank=0;
269 : }
270 :
271 5952 : double contact_sum = 0.0;
272 :
273 5952 : const unsigned int nn=nl->size();
274 :
275 321408 : for(unsigned int i=rank; i<nn; i+=stride) {
276 315456 : Vector distance;
277 315456 : unsigned int i0=nl->getClosePair(i).first;
278 315456 : unsigned int i1=nl->getClosePair(i).second;
279 315456 : if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) {
280 0 : continue;
281 : }
282 :
283 315456 : if(pbc_) {
284 86304 : distance=pbcDistance(getPosition(i0),getPosition(i1));
285 : } else {
286 229152 : distance=delta(getPosition(i0),getPosition(i1));
287 : }
288 :
289 315456 : double exp_arg = exp(-(distance.modulo()-r_globalshift_)*inv_r_eff_);
290 315456 : contact_sum += exp_arg;
291 :
292 315456 : exp_arg /= distance.modulo();
293 : Vector dd(exp_arg*distance);
294 315456 : Tensor vv(dd,distance);
295 315456 : deriv[i0]-=dd;
296 315456 : deriv[i1]+=dd;
297 : virial-=vv;
298 : }
299 :
300 5952 : if(!serial_) {
301 5952 : comm.Sum(contact_sum);
302 5952 : if(!deriv.empty()) {
303 5952 : comm.Sum(&deriv[0][0],3*deriv.size());
304 : }
305 5952 : comm.Sum(virial);
306 : }
307 :
308 5952 : double value = tanh(total_prefactor_*contact_sum);
309 : // using that d/dx[tanh(x)]= 1-[tanh(x)]^2
310 5952 : double deriv_f = -inv_r_eff_*total_prefactor_*(1.0-value*value);
311 5952 : value -= offset_c_;
312 :
313 476160 : for(unsigned int i=0; i<deriv.size(); ++i) {
314 470208 : setAtomsDerivatives(i,deriv_f*deriv[i]);
315 : }
316 5952 : setValue(value);
317 11904 : setBoxDerivatives(deriv_f*virial);
318 :
319 5952 : }
320 : }
321 : }
|
