Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2024 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 :
23 : #include "colvar/Colvar.h"
24 : #include "core/ActionRegister.h"
25 : #include "core/PlumedMain.h"
26 : #include "tools/Communicator.h"
27 : #include "tools/Matrix.h"
28 : #include "core/GenericMolInfo.h"
29 : #include "core/ActionSet.h"
30 : #include "tools/File.h"
31 : #include "tools/OpenMP.h"
32 : #include <string>
33 : #include <cmath>
34 : #include <map>
35 : #include <ctime>
36 : #include "tools/Random.h"
37 :
38 : #ifndef M_PI
39 : //why do not use PLMD::pi?
40 : #define M_PI 3.14159265358979323846
41 : #endif
42 :
43 : namespace PLMD {
44 : namespace isdb {
45 :
46 : //+PLUMEDOC ISDB_COLVAR BAIES
47 : /*
48 : Bayesian refinement of AF models.
49 :
50 : This action implements the Bayesian approach to refine AF models introduced here and here.
51 : It can be used to generate conformational ensembles of IDPs or refine AF models prior to small-molecule virtual screening.
52 :
53 : ## Examples
54 :
55 : Complete tutorials can be found <a href="https://github.com/COSBlab/bAIes-IDP">here</a> and here.
56 :
57 : */
58 : //+ENDPLUMEDOC
59 :
60 : class BAIES : public Colvar {
61 :
62 : private:
63 :
64 : // bool pbc
65 : bool pbc_;
66 : // temperature in kbt
67 : double kbt_;
68 : // number of threads
69 : unsigned nt_;
70 : // positions
71 : std::vector<Vector> pos_;
72 : // derivatives
73 : std::vector<Vector> atom_der_;
74 : // AF2 parameters variables
75 : std::string mtype_;
76 : std::vector < std::pair<unsigned, unsigned> > atom_pairs_;
77 : std::vector<double> mus_;
78 : std::vector<double> sigmas_;
79 : // constants
80 : double inv_sqrt2_, sqrt2_pi_, inv_pi2_;
81 : // prior type
82 : unsigned prior_;
83 : enum { NONE, JEFFREYS, CAUCHY };
84 : // private functions
85 : void read_data_file(std::string datafile, std::vector<AtomNumber> atoms, double smin);
86 : // energie
87 : double getGauss();
88 : double getGaussJeffreys();
89 : double getGaussCauchy();
90 : double getLognorm();
91 : double getLognormJeffreys();
92 : double getLognormCauchy();
93 :
94 : public:
95 : static void registerKeywords( Keywords& keys );
96 : explicit BAIES(const ActionOptions&);
97 : // active methods:
98 : void calculate() override;
99 : };
100 :
101 : PLUMED_REGISTER_ACTION(BAIES,"BAIES")
102 :
103 10 : void BAIES::registerKeywords( Keywords& keys ) {
104 10 : Colvar::registerKeywords( keys );
105 10 : keys.add("atoms","ATOMS","atoms used in the calculation of bAIes energy");
106 10 : keys.add("compulsory","DATA_FILE","file with AF2 fit parameters");
107 10 : keys.add("compulsory","PRIOR", "type of prior to use (NONE, JEFFREYS, CAUCHY");
108 10 : keys.add("optional","TEMP", "temperature of the system, by default the one passed by the MD engine");
109 10 : keys.add("optional","SIGMA_MIN", "minimum value of sigma");
110 20 : keys.addOutputComponent("ene","default","scalar","Bayesian bAIes energy");
111 10 : }
112 :
113 6 : BAIES::BAIES(const ActionOptions&ao):
114 : PLUMED_COLVAR_INIT(ao),
115 6 : pbc_(true) {
116 : // set constants
117 6 : inv_sqrt2_ = 1.0/sqrt(2.0);
118 6 : sqrt2_pi_ = sqrt(2.0 / M_PI);
119 6 : inv_pi2_ = 0.5 / M_PI / M_PI;
120 :
121 : // list of atoms
122 6 : std::vector<AtomNumber> atoms;
123 12 : parseAtomList("ATOMS", atoms);
124 :
125 : // file with AF2 fit parameters
126 : std::string datafile;
127 12 : parse("DATA_FILE", datafile);
128 :
129 : // prior type
130 : std::string prior;
131 6 : parse("PRIOR", prior);
132 : // check priors allowed
133 6 : if(prior=="NONE") {
134 2 : prior_ = NONE;
135 4 : } else if(prior=="JEFFREYS") {
136 2 : prior_ = JEFFREYS;
137 2 : } else if(prior=="CAUCHY") {
138 2 : prior_ = CAUCHY;
139 : } else {
140 0 : error("Unknown PRIOR type!");
141 : }
142 :
143 6 : bool nopbc=!pbc_;
144 6 : parseFlag("NOPBC",nopbc);
145 6 : pbc_=!nopbc;
146 :
147 : // set kbt
148 : // getkBT reads TEMP (in kelvin) and converts it to energy units
149 6 : kbt_ = getkBT();
150 6 : if(kbt_<=0.0) {
151 0 : error("Unless the MD engine passes the temperature to plumed, you must specify it using TEMP");
152 : }
153 :
154 : // set sigma_min
155 6 : double sigma_min = -1.0;
156 6 : parse("SIGMA_MIN",sigma_min);
157 :
158 : // check read
159 6 : checkRead();
160 :
161 : // set parallel stuff
162 6 : unsigned mpisize=comm.Get_size();
163 6 : if(mpisize>1) {
164 0 : error("BAIES supports only OpenMP parallelization");
165 : }
166 : // set number of OpenMP threads
167 6 : nt_ = OpenMP::getNumThreads();
168 :
169 : // read the data file
170 12 : read_data_file(datafile, atoms, sigma_min);
171 :
172 : // print stuff to log
173 6 : log.printf(" number of atoms involved : %zu\n", atoms.size());
174 6 : log.printf(" AF2 fit parameters file : %s\n", datafile.c_str());
175 6 : log.printf(" prior type : %s\n", prior.c_str());
176 : // print stuff from data file
177 6 : log.printf(" noise model type : %s\n", mtype_.c_str());
178 6 : log.printf(" number of pairs involved : %zu\n", atom_pairs_.size());
179 : // other info
180 6 : log.printf(" temperature of the system in energy unit : %f\n", kbt_);
181 6 : log.printf(" minimum value of sigma : %f\n", sigma_min);
182 6 : if(pbc_) {
183 6 : log.printf(" using periodic boundary conditions\n");
184 : } else {
185 0 : log.printf(" without periodic boundary conditions\n");
186 : }
187 :
188 : // prepare other vectors: data and derivatives
189 6 : atom_der_.resize(atoms.size());
190 :
191 : // add components
192 12 : addComponentWithDerivatives("ene");
193 6 : componentIsNotPeriodic("ene");
194 :
195 : // request atoms
196 6 : requestAtoms(atoms);
197 :
198 : // print bibliography
199 6 : log<<" Bibliography "<<plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)");
200 6 : log<<plumed.cite("Schnapka, Morozova, Sen, Bonomi, bioRxiv (2025) doi: XXX");
201 6 : log<<"\n";
202 6 : }
203 :
204 6 : void BAIES::read_data_file(const std::string datafile, const std::vector<AtomNumber> atoms, double smin) {
205 : unsigned id, ai, aj;
206 : double mu, sigma;
207 : // map serials to index in position array
208 : std::map< unsigned, unsigned > index;
209 204 : for(unsigned i=0; i<atoms.size(); ++i) {
210 198 : index[atoms[i].serial()]=i;
211 : }
212 : // open file
213 : auto ifile=Tools::make_unique<IFile>();
214 6 : if(ifile->FileExist(datafile)) {
215 6 : ifile->open(datafile);
216 : // read constant fields
217 12 : if( ifile->FieldExist("model") ) {
218 6 : ifile->scanField("model",mtype_);
219 6 : if( mtype_!="gaussian" && mtype_!="lognormal" ) {
220 0 : error("Unknown noise model type");
221 : }
222 : } else {
223 0 : error("Missing noise model type in DATA_FILE");
224 : }
225 : // position of an atom (given by its serial) in the ATOMS list
226 2112 : auto getIndex=[&](const unsigned serial) {
227 2112 : const auto it=index.find(serial);
228 2112 : if(it==index.end()) {
229 : std::string s;
230 0 : Tools::convert(serial,s);
231 0 : error("atom "+s+" is in DATA_FILE but not in ATOMS");
232 : }
233 2112 : return it->second;
234 6 : };
235 : // read line-by-line
236 2124 : while(ifile->scanField("Id",id)) {
237 1056 : ifile->scanField("atom_i",ai);
238 1056 : ifile->scanField("atom_j",aj);
239 1056 : ifile->scanField("mu",mu);
240 1056 : ifile->scanField("sigma",sigma);
241 1056 : if(ai==aj) {
242 0 : error("atom_i and atom_j must be different in DATA_FILE");
243 : }
244 : // list of pairs of atoms
245 1056 : atom_pairs_.push_back(std::make_pair(getIndex(ai),getIndex(aj)));
246 : // list of mus
247 1056 : mus_.push_back(mu);
248 : // list of sigmas
249 2035 : sigmas_.push_back(std::max(smin,sigma));
250 1056 : if(!(sigmas_.back()>0.)) {
251 0 : error("sigma must be positive in DATA_FILE (or use SIGMA_MIN)");
252 : }
253 : // new line
254 1056 : ifile->scanField();
255 : }
256 6 : ifile->close();
257 : } else {
258 0 : error("Cannot find DATA_FILE "+datafile+"\n");
259 : }
260 12 : }
261 :
262 : // Energy models
263 6 : double BAIES::getGauss() {
264 : double ene = 0.0;
265 : // loop over atoms pairs
266 6 : #pragma omp parallel num_threads(nt_) shared(ene)
267 : {
268 : std::vector<Vector> omp_der(atom_der_.size());
269 : #pragma omp for reduction( + : ene) nowait
270 : for(unsigned i=0; i<atom_pairs_.size(); ++i) {
271 : // get indexes
272 : unsigned ii = atom_pairs_[i].first;
273 : unsigned jj = atom_pairs_[i].second;
274 : // get distance
275 : Vector distance=delta(pos_[ii], pos_[jj]);
276 : const double value=distance.modulo();
277 : const double invvalue=1.0/value;
278 : // get energy
279 : double arg = ( value - mus_[i] ) / sigmas_[i];
280 : ene += arg * arg;
281 : // calculate derivatives
282 : double arg_der = kbt_ * arg / sigmas_[i];
283 : omp_der[ii] += - arg_der * invvalue * distance;
284 : omp_der[jj] += arg_der * invvalue * distance;
285 : }
286 : // gather from tasks
287 : #pragma omp critical
288 : for(unsigned i=0; i<atom_der_.size(); ++i) {
289 : atom_der_[i]+=omp_der[i];
290 : }
291 : }
292 : // missing Gaussian normalization constant
293 6 : return 0.5 * kbt_ * ene;
294 : }
295 :
296 6 : double BAIES::getGaussJeffreys() {
297 : double ene = 0.0;
298 : // loop over all atom pairs
299 6 : #pragma omp parallel num_threads(nt_) shared(ene)
300 : {
301 : std::vector<Vector> omp_der(atom_der_.size());
302 : #pragma omp for reduction( + : ene) nowait
303 : for(unsigned i=0; i<atom_pairs_.size(); ++i) {
304 : // get indexes
305 : unsigned ii = atom_pairs_[i].first;
306 : unsigned jj = atom_pairs_[i].second;
307 : // get distance
308 : Vector distance=delta(pos_[ii], pos_[jj]);
309 : const double value=distance.modulo();
310 : const double invvalue=1.0/value;
311 : // get energy
312 : double val_mu = value - mus_[i];
313 : double inv_val_mu = 1.0 / val_mu;
314 : double inv_sigma = 1.0 / sigmas_[i];
315 : double erfcontent = val_mu * inv_sigma * inv_sqrt2_;
316 : double erfval = std::erf(erfcontent);
317 : // increase energy
318 : ene += -std::log( erfval * inv_val_mu );
319 : // calculate derivatives
320 : double arg_der = - kbt_ * val_mu * ( (sqrt2_pi_ * std::exp(-erfcontent * erfcontent) * inv_val_mu * inv_sigma ) - ( erfval * inv_val_mu * inv_val_mu ) ) / erfval ;
321 : // increase derivatives
322 : omp_der[ii] += - arg_der * invvalue * distance;
323 : omp_der[jj] += arg_der * invvalue * distance;
324 : }
325 : // gather from tasks
326 : #pragma omp critical
327 : for(unsigned i=0; i<atom_der_.size(); ++i) {
328 : atom_der_[i]+=omp_der[i];
329 : }
330 : }
331 : // missing normalization constants
332 6 : return kbt_ * ene;
333 : }
334 :
335 6 : double BAIES::getGaussCauchy() {
336 : double ene = 0.0;
337 : // loop over all atom pairs
338 6 : #pragma omp parallel num_threads(nt_) shared(ene)
339 : {
340 : std::vector<Vector> omp_der(atom_der_.size());
341 : #pragma omp for reduction( + : ene) nowait
342 : for(unsigned i=0; i<atom_pairs_.size(); ++i) {
343 : // get indexes
344 : unsigned ii = atom_pairs_[i].first;
345 : unsigned jj = atom_pairs_[i].second;
346 : // get distance
347 : Vector distance=delta(pos_[ii], pos_[jj]);
348 : const double value=distance.modulo();
349 : const double invvalue=1.0/value;
350 : // get energy
351 : double val_mu = value - mus_[i];
352 : double lncontent = 0.5*val_mu*val_mu + sigmas_[i]*sigmas_[i];
353 : ene += std::log( lncontent );
354 : // calculate derivatives
355 : double arg_der = kbt_ * val_mu / lncontent;
356 : omp_der[ii] += - arg_der * invvalue * distance;
357 : omp_der[jj] += arg_der * invvalue * distance;
358 : }
359 : // gather from tasks
360 : #pragma omp critical
361 : for(unsigned i=0; i<atom_der_.size(); ++i) {
362 : atom_der_[i]+=omp_der[i];
363 : }
364 : }
365 6 : return kbt_ * ene;
366 : }
367 :
368 6 : double BAIES::getLognorm() {
369 : double ene = 0.0;
370 : // loop over atoms pairs
371 6 : #pragma omp parallel num_threads(nt_) shared(ene)
372 : {
373 : std::vector<Vector> omp_der(atom_der_.size());
374 : #pragma omp for reduction( + : ene) nowait
375 : for(unsigned i=0; i<atom_pairs_.size(); ++i) {
376 : // get indexes
377 : unsigned ii = atom_pairs_[i].first;
378 : unsigned jj = atom_pairs_[i].second;
379 : // get distance
380 : Vector distance=delta(pos_[ii], pos_[jj]);
381 : const double value=distance.modulo();
382 : const double invvalue=1.0/value;
383 : // get energy
384 : double arg = ( std::log(value) - mus_[i] ) / sigmas_[i];
385 : ene += arg * arg;
386 : // calculate derivatives
387 : double arg_der = kbt_ * arg / (value * sigmas_[i]);
388 : omp_der[ii] += - arg_der * invvalue * distance;
389 : omp_der[jj] += arg_der * invvalue * distance;
390 : }
391 : // gather from tasks
392 : #pragma omp critical
393 : for(unsigned i=0; i<atom_der_.size(); ++i) {
394 : atom_der_[i]+=omp_der[i];
395 : }
396 : }
397 6 : return 0.5 * kbt_ * ene;
398 : }
399 :
400 6 : double BAIES::getLognormJeffreys() {
401 : double ene = 0.0;
402 : // loop over all atom pairs
403 6 : #pragma omp parallel num_threads(nt_) shared(ene)
404 : {
405 : std::vector<Vector> omp_der(atom_der_.size());
406 : #pragma omp for reduction( + : ene) nowait
407 : for(unsigned i=0; i<atom_pairs_.size(); ++i) {
408 : // get indexes
409 : unsigned ii = atom_pairs_[i].first;
410 : unsigned jj = atom_pairs_[i].second;
411 : // get distance
412 : Vector distance=delta(pos_[ii], pos_[jj]);
413 : const double value=distance.modulo();
414 : const double invvalue=1.0/value;
415 : // get energy
416 : double val_mu = std::log(value) - mus_[i];
417 : double inv_val_mu = 1 / val_mu;
418 : double inv_sigma = 1 / sigmas_[i];
419 : double erfcontent = val_mu * inv_sigma * inv_sqrt2_;
420 : double erfval = std::erf(erfcontent);
421 : // increase energy
422 : ene += -std::log( erfval * inv_val_mu );
423 : // calculate derivatives
424 : double arg_der = - kbt_ * val_mu * ( (sqrt2_pi_ * std::exp(-erfcontent * erfcontent) * inv_val_mu * inv_sigma ) - ( erfval * inv_val_mu * inv_val_mu ) ) / erfval ;
425 : // chain rule
426 : arg_der *= invvalue;
427 : // increase derivatives
428 : omp_der[ii] += - arg_der * invvalue * distance;
429 : omp_der[jj] += arg_der * invvalue * distance;
430 : }
431 : // gather from tasks
432 : #pragma omp critical
433 : for(unsigned i=0; i<atom_der_.size(); ++i) {
434 : atom_der_[i]+=omp_der[i];
435 : }
436 : }
437 : // missing constants
438 6 : return kbt_ * ene;
439 : }
440 :
441 6 : double BAIES::getLognormCauchy() {
442 : double ene = 0.0;
443 : // loop over all atom pairs
444 6 : #pragma omp parallel num_threads(nt_) shared(ene)
445 : {
446 : std::vector<Vector> omp_der(atom_der_.size());
447 : #pragma omp for reduction( + : ene) nowait
448 : for(unsigned i=0; i<atom_pairs_.size(); ++i) {
449 : // get indexes
450 : unsigned ii = atom_pairs_[i].first;
451 : unsigned jj = atom_pairs_[i].second;
452 : // get distance
453 : Vector distance=delta(pos_[ii], pos_[jj]);
454 : const double value=distance.modulo();
455 : const double invvalue=1.0/value;
456 : // get energy
457 : double val_mu = std::log(value) - mus_[i];
458 : double lncontent = 0.5*val_mu*val_mu + sigmas_[i]*sigmas_[i];
459 : ene += std::log( lncontent );
460 : // calculate derivatives
461 : double arg_der = kbt_ * val_mu * invvalue / lncontent;
462 : omp_der[ii] += - arg_der * invvalue * distance;
463 : omp_der[jj] += arg_der * invvalue * distance;
464 : }
465 : // gather from tasks
466 : #pragma omp critical
467 : for(unsigned i=0; i<atom_der_.size(); ++i) {
468 : atom_der_[i]+=omp_der[i];
469 : }
470 : }
471 : // missing constants
472 6 : return kbt_ * ene;
473 : }
474 :
475 36 : void BAIES::calculate() {
476 : // fix PBCs
477 36 : if(pbc_) {
478 36 : makeWhole();
479 : }
480 :
481 : // get the positions at this step
482 36 : pos_ = getPositions();
483 :
484 : // reset derivatives
485 36 : #pragma omp parallel for num_threads(nt_)
486 : for(unsigned i=0; i<atom_der_.size(); ++i) {
487 : atom_der_[i] = Vector(0,0,0);
488 : }
489 :
490 : // calculate bAIes energy and derivatives
491 : double ene = 0.0;
492 : // Gaussian noise model
493 36 : if(mtype_=="gaussian") {
494 18 : switch(prior_) {
495 6 : case NONE:
496 6 : ene = getGauss();
497 : break;
498 6 : case JEFFREYS:
499 6 : ene = getGaussJeffreys();
500 : break;
501 6 : case CAUCHY:
502 6 : ene = getGaussCauchy();
503 : break;
504 : }
505 : } else {
506 18 : switch(prior_) {
507 6 : case NONE:
508 6 : ene = getLognorm();
509 : break;
510 6 : case JEFFREYS:
511 6 : ene = getLognormJeffreys();
512 : break;
513 6 : case CAUCHY:
514 6 : ene = getLognormCauchy();
515 : break;
516 : }
517 : }
518 : // set score
519 72 : Value* score = getPntrToComponent("ene");
520 : score->set(ene);
521 : // calculate virial
522 36 : Tensor virial;
523 : // declare omp reduction for Tensors
524 : #pragma omp declare reduction( sumTensor : Tensor : omp_out += omp_in )
525 :
526 36 : #pragma omp parallel for num_threads(nt_) reduction (sumTensor : virial)
527 : for(unsigned i=0; i<pos_.size(); ++i) {
528 : virial += Tensor(pos_[i], -atom_der_[i]);
529 : }
530 : // set virial
531 36 : setBoxDerivatives(score, virial);
532 : // set derivatives
533 36 : #pragma omp parallel for num_threads(nt_)
534 : for(unsigned i=0; i<atom_der_.size(); ++i) {
535 : setAtomsDerivatives(score, i, atom_der_[i]);
536 : }
537 36 : }
538 :
539 : }
540 : }
|
