Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 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 :
23 : #ifdef __PLUMED_HAS_LIBTORCH
24 : #include "colvar/Colvar.h"
25 : #include "core/ActionRegister.h"
26 : #include "core/PlumedMain.h"
27 : #include "tools/Communicator.h"
28 : #include "tools/Matrix.h"
29 : #include "core/GenericMolInfo.h"
30 : #include "core/ActionSet.h"
31 : #include "tools/File.h"
32 : #include "tools/OpenMP.h"
33 : #include <string>
34 : #include <cmath>
35 : #include <map>
36 : #include <numeric>
37 : #include <ctime>
38 : #include <algorithm>
39 : #include <utility>
40 : #include "tools/Random.h"
41 :
42 : //this avoids nasty warnings on which we have no control
43 : //Here I am using the pragma for gcc/clang, for visual studio you'll need something different
44 : //This should be the correct direction,
45 : //but the warning that it is blocking everithing is:
46 : // error: ISO C++11 requires at least one argument for the "..." in a variadic macro [-Werror]
47 : //that is not named like [-Werror=shadow]
48 : //so , at time of writing this I did not found any solution but deactivating warning in the makefile
49 :
50 : #pragma GCC diagnostic push
51 : #pragma GCC diagnostic ignored "-Wshadow"
52 :
53 : #include <torch/torch.h>
54 : #include <torch/script.h>
55 :
56 : #pragma GCC diagnostic push
57 :
58 : #ifndef M_PI
59 : #define M_PI 3.14159265358979323846
60 : #endif
61 :
62 : namespace PLMD {
63 : namespace isdb {
64 :
65 : //+PLUMEDOC ISDB_COLVAR EMMIVOX
66 : /*
67 : Bayesian single-structure and ensemble refinement with cryo-EM maps.
68 :
69 : This action implements the Bayesian approach for single-structure and ensemble refinement from cryo-EM maps introduced <a href="https://doi.org/10.1371/journal.pcbi.1012180">here</a>.
70 : EMMIVox does not require fitting the cryo-EM map with a Gaussian Mixture Model, as done in [EMMI](EMMI.md), but uses directly the voxels in the deposited map.
71 :
72 : When run in single-replica mode, this action allows atomistic, flexible refinement (and B-factors inference) of an individual structure into a density map.
73 : A coarse-grained forward model can also be used in combination with the MARTINI force field.
74 : Combined with a multi-replica framework (such as the -multi option in GROMACS), the user can model an ensemble of structures using
75 : the Metainference approach that is introduced in the paper cited below. The approach can be used to model continous dynamics of flexible regions as well as semi-ordered waters, lipids, and ions.
76 :
77 : !!! note "installing libtorch"
78 :
79 : To use EMMIVOX, PLUMED must be linked against the LibTorch library using the instructions [on the module page](module_isdb.md).
80 :
81 : ## Examples
82 :
83 : Complete tutorials for single-structure and ensemble refinement can be found <a href="https://github.com/COSBlab/EMMIVox">here</a>.
84 :
85 : */
86 : //+ENDPLUMEDOC
87 :
88 : class EMMIVOX : public Colvar {
89 :
90 : private:
91 :
92 : // temperature in kbt
93 : double kbt_;
94 : // model - atom types
95 : std::vector<unsigned> Model_type_;
96 : // model - list of atom sigmas - one per atom type
97 : std::vector<Vector5d> Model_s_;
98 : // model - list of atom weights - one per atom type
99 : std::vector<Vector5d> Model_w_;
100 : // model - map between residue/chain IDs and list of atoms
101 : std::map< std::pair<int,std::string>, std::vector<unsigned> > Model_resmap_;
102 : // model - list of residue/chain IDs per atom
103 : std::vector< std::pair<int,std::string> > Model_res_;
104 : // model - list of neighboring voxels per atom
105 : std::vector< std::vector<unsigned> > Model_nb_;
106 : // model - map between residue/chain ID and bfactor
107 : std::map< std::pair<int, std::string>, double> Model_b_;
108 : // model - global list of residue/chain IDs
109 : std::vector< std::pair<int,std::string> > Model_rlist_;
110 : // model density
111 : std::vector<double> ovmd_;
112 :
113 : // data map - voxel position
114 : std::vector<Vector> Map_m_;
115 : // data map - density
116 : std::vector<double> ovdd_;
117 : // data map - error
118 : std::vector<double> exp_err_;
119 :
120 : // derivatives
121 : std::vector<Vector> ovmd_der_;
122 : std::vector<Vector> atom_der_;
123 : std::vector<double> score_der_;
124 : // constants
125 : double inv_sqrt2_, sqrt2_pi_, inv_pi2_;
126 : std::vector<Vector5d> pref_;
127 : std::vector<Vector5d> invs2_;
128 : std::vector<Vector5d> cfact_;
129 : std::vector<double> cut_;
130 : // metainference
131 : unsigned nrep_;
132 : unsigned replica_;
133 : std::vector<double> ismin_;
134 : // neighbor list
135 : double nl_dist_cutoff_;
136 : double nl_gauss_cutoff_;
137 : unsigned nl_stride_;
138 : bool first_time_;
139 : std::vector< std::pair<unsigned,unsigned> > nl_;
140 : std::vector< std::pair<unsigned,unsigned> > ns_;
141 : std::vector<Vector> refpos_;
142 : // averaging
143 : bool no_aver_;
144 : // correlation;
145 : bool do_corr_;
146 : // Monte Carlo stuff
147 : Random random_;
148 : // Scale and Offset
149 : double scale_;
150 : double offset_;
151 : // Bfact Monte Carlo
152 : double dbfact_;
153 : double bfactmin_;
154 : double bfactmax_;
155 : double bfactsig_;
156 : bool bfactnoc_;
157 : bool bfactread_;
158 : int MCBstride_;
159 : double MCBaccept_;
160 : double MCBtrials_;
161 : // residue neighbor list
162 : std::vector< std::vector<unsigned> > nl_res_;
163 : bool bfactemin_;
164 : // Martini scattering factors
165 : bool martini_;
166 : // status stuff
167 : unsigned int statusstride_;
168 : std::string statusfilename_;
169 : OFile statusfile_;
170 : bool first_status_;
171 : // total energy and virial
172 : double ene_;
173 : Tensor virial_;
174 : double eps_;
175 : // model density file
176 : unsigned int mapstride_;
177 : std::string mapfilename_;
178 : // Libtorch stuff
179 : bool gpu_;
180 : torch::Tensor ovmd_gpu_;
181 : torch::Tensor ovmd_der_gpu_;
182 : torch::Tensor ismin_gpu_;
183 : torch::Tensor ovdd_gpu_;
184 : torch::Tensor Map_m_gpu_;
185 : torch::Tensor pref_gpu_;
186 : torch::Tensor invs2_gpu_;
187 : torch::Tensor nl_id_gpu_;
188 : torch::Tensor nl_im_gpu_;
189 : torch::Tensor pref_nl_gpu_;
190 : torch::Tensor invs2_nl_gpu_;
191 : torch::Tensor Map_m_nl_gpu_;
192 : torch::DeviceType device_t_;
193 : //
194 : // write file with model density
195 : void write_model_density(long int step);
196 : // get median of vector
197 : double get_median(std::vector<double> v);
198 : // read and write status
199 : void read_status();
200 : void print_status(long int step);
201 : // accept or reject
202 : bool doAccept(double oldE, double newE, double kbt);
203 : // vector of close residues
204 : void get_close_residues();
205 : // do MonteCarlo for Bfactor
206 : void doMonteCarloBfact();
207 : // calculate model parameters
208 : std::vector<double> get_Model_param(std::vector<AtomNumber> &atoms);
209 : // read data file
210 : void get_exp_data(const std::string &datafile);
211 : // auxiliary methods
212 : void prepare_gpu();
213 : void initialize_Bfactor(double reso);
214 : void get_auxiliary_vectors();
215 : void push_auxiliary_gpu();
216 : // calculate overlap between two Gaussians
217 : double get_overlap(const Vector &d_m, const Vector &m_m,
218 : const Vector5d &cfact, const Vector5d &m_s, double bfact);
219 : // update the neighbor list
220 : void update_neighbor_list();
221 : // update data on device
222 : void update_gpu();
223 : // update the neighbor sphere
224 : void update_neighbor_sphere();
225 : bool do_neighbor_sphere();
226 : // calculate forward model and score on device
227 : void calculate_fmod();
228 : void calculate_score();
229 : // calculate correlation
230 : void calculate_corr();
231 :
232 : public:
233 : static void registerKeywords( Keywords& keys );
234 : explicit EMMIVOX(const ActionOptions&);
235 : // active methods:
236 : void prepare() override;
237 : void calculate() override;
238 : };
239 :
240 : PLUMED_REGISTER_ACTION(EMMIVOX,"EMMIVOX")
241 :
242 6 : void EMMIVOX::registerKeywords( Keywords& keys ) {
243 6 : Colvar::registerKeywords( keys );
244 6 : keys.add("atoms","ATOMS","atoms used in the calculation of the density map, typically all heavy atoms");
245 6 : keys.add("compulsory","DATA_FILE","file with cryo-EM map");
246 6 : keys.add("compulsory","RESOLUTION", "cryo-EM map resolution");
247 6 : keys.add("compulsory","NORM_DENSITY","integral of experimental density");
248 6 : keys.add("compulsory","WRITE_STRIDE","stride for writing status file");
249 6 : keys.add("optional","NL_DIST_CUTOFF","neighbor list distance cutoff");
250 6 : keys.add("optional","NL_GAUSS_CUTOFF","neighbor list Gaussian sigma cutoff");
251 6 : keys.add("optional","NL_STRIDE","neighbor list update frequency");
252 6 : keys.add("optional","SIGMA_MIN","minimum density error");
253 6 : keys.add("optional","DBFACT","Bfactor MC step");
254 6 : keys.add("optional","BFACT_MAX","Bfactor maximum value");
255 6 : keys.add("optional","MCBFACT_STRIDE", "Bfactor MC stride");
256 6 : keys.add("optional","BFACT_SIGMA","Bfactor sigma prior");
257 6 : keys.add("optional","STATUS_FILE","write a file with all the data useful for restart");
258 6 : keys.add("optional","SCALE","scale factor");
259 6 : keys.add("optional","OFFSET","offset");
260 6 : keys.add("optional","TEMP","temperature");
261 6 : keys.add("optional","WRITE_MAP","file with model density");
262 6 : keys.add("optional","WRITE_MAP_STRIDE","stride for writing model density to file");
263 6 : keys.addFlag("NO_AVER",false,"no ensemble averaging in multi-replica mode");
264 6 : keys.addFlag("CORRELATION",false,"calculate correlation coefficient");
265 6 : keys.addFlag("GPU",false,"calculate EMMIVOX on GPU with Libtorch");
266 6 : keys.addFlag("BFACT_NOCHAIN",false,"Do not use chain ID for Bfactor MC");
267 6 : keys.addFlag("BFACT_READ",false,"Read Bfactor on RESTART (automatic with DBFACT>0)");
268 6 : keys.addFlag("BFACT_MINIMIZE",false,"Accept only moves that decrease energy");
269 6 : keys.addFlag("MARTINI",false,"Use Martini scattering factors");
270 6 : keys.add("optional","SEED","seed of the random number generator used by the Monte Carlo sampling of the Bfactors (by default it is taken from the clock at every run)");
271 12 : keys.addOutputComponent("scoreb","default","scalar","Bayesian score");
272 12 : keys.addOutputComponent("scale", "default","scalar","scale factor");
273 12 : keys.addOutputComponent("offset","default","scalar","offset");
274 12 : keys.addOutputComponent("accB", "default","scalar", "Bfactor MC acceptance");
275 12 : keys.addOutputComponent("kbt", "default","scalar", "temperature in energy unit");
276 12 : keys.addOutputComponent("corr", "CORRELATION","scalar", "correlation coefficient");
277 6 : keys.addDOI("10.1126/sciadv.1501177");
278 6 : }
279 :
280 2 : EMMIVOX::EMMIVOX(const ActionOptions&ao):
281 : PLUMED_COLVAR_INIT(ao),
282 2 : nl_dist_cutoff_(1.0), nl_gauss_cutoff_(3.0), nl_stride_(50),
283 2 : first_time_(true), no_aver_(false), do_corr_(false),
284 2 : scale_(1.), offset_(0.),
285 2 : dbfact_(0.0), bfactmin_(0.05), bfactmax_(5.0),
286 2 : bfactsig_(0.1), bfactnoc_(false), bfactread_(false),
287 2 : MCBstride_(1), MCBaccept_(0.), MCBtrials_(0.), bfactemin_(false),
288 2 : martini_(false), statusstride_(0), first_status_(true),
289 2 : eps_(0.0001), mapstride_(0), gpu_(false) {
290 : // set constants
291 2 : inv_sqrt2_ = 1.0/sqrt(2.0);
292 2 : sqrt2_pi_ = sqrt(2.0 / M_PI);
293 2 : inv_pi2_ = 0.5 / M_PI / M_PI;
294 :
295 : // list of atoms
296 2 : std::vector<AtomNumber> atoms;
297 4 : parseAtomList("ATOMS", atoms);
298 :
299 : // file with experimental cryo-EM map
300 : std::string datafile;
301 2 : parse("DATA_FILE", datafile);
302 :
303 : // neighbor list cutoffs
304 2 : parse("NL_DIST_CUTOFF",nl_dist_cutoff_);
305 2 : parse("NL_GAUSS_CUTOFF",nl_gauss_cutoff_);
306 : // checks
307 2 : if(nl_dist_cutoff_<=0. && nl_gauss_cutoff_<=0.) {
308 0 : error("You must specify either NL_DIST_CUTOFF or NL_GAUSS_CUTOFF or both");
309 : }
310 2 : if(nl_gauss_cutoff_<=0.) {
311 0 : nl_gauss_cutoff_ = 1.0e+10;
312 : }
313 2 : if(nl_dist_cutoff_<=0.) {
314 0 : nl_dist_cutoff_ = 1.0e+10;
315 : }
316 : // neighbor list update stride
317 2 : parse("NL_STRIDE",nl_stride_);
318 2 : if(nl_stride_<=0) {
319 0 : error("NL_STRIDE must be explicitly specified and positive");
320 : }
321 :
322 : // minimum value for error
323 2 : double sigma_min = 0.2;
324 2 : parse("SIGMA_MIN", sigma_min);
325 2 : if(sigma_min<0.) {
326 0 : error("SIGMA_MIN must be greater or equal to zero");
327 : }
328 :
329 : // status file parameters
330 2 : parse("WRITE_STRIDE", statusstride_);
331 2 : if(statusstride_<=0) {
332 0 : error("you must specify a positive WRITE_STRIDE");
333 : }
334 2 : parse("STATUS_FILE", statusfilename_);
335 2 : if(statusfilename_=="") {
336 0 : statusfilename_ = "EMMIStatus"+getLabel();
337 : }
338 :
339 : // integral of the experimetal density
340 : double norm_d;
341 2 : parse("NORM_DENSITY", norm_d);
342 :
343 : // temperature
344 2 : kbt_ = getkBT();
345 :
346 : // scale and offset
347 2 : parse("SCALE", scale_);
348 2 : parse("OFFSET",offset_);
349 :
350 : // B-factors MC
351 2 : parse("DBFACT",dbfact_);
352 : // read Bfactors
353 2 : parseFlag("BFACT_READ",bfactread_);
354 : // do not use chains
355 2 : parseFlag("BFACT_NOCHAIN",bfactnoc_);
356 : // other parameters
357 2 : if(dbfact_>0.) {
358 1 : parse("MCBFACT_STRIDE",MCBstride_);
359 1 : parse("BFACT_MAX",bfactmax_);
360 1 : parse("BFACT_SIGMA",bfactsig_);
361 1 : parseFlag("BFACT_MINIMIZE",bfactemin_);
362 : // checks
363 1 : if(MCBstride_<=0) {
364 0 : error("you must specify a positive MCBFACT_STRIDE");
365 : }
366 1 : if(bfactmax_<=bfactmin_) {
367 0 : error("you must specify a positive BFACT_MAX");
368 : }
369 1 : if(MCBstride_%nl_stride_!=0) {
370 0 : error("MCBFACT_STRIDE must be multiple of NL_STRIDE");
371 : }
372 1 : if(bfactsig_<=0.) {
373 0 : error("you must specify a positive BFACT_SIGMA");
374 : }
375 : }
376 :
377 : // read map resolution
378 : double reso;
379 2 : parse("RESOLUTION", reso);
380 2 : if(reso<=0.) {
381 0 : error("RESOLUTION must be strictly positive");
382 : }
383 :
384 : // averaging or not
385 2 : parseFlag("NO_AVER",no_aver_);
386 :
387 : // calculate correlation coefficient
388 2 : parseFlag("CORRELATION",do_corr_);
389 :
390 : // write density file
391 2 : parse("WRITE_MAP_STRIDE", mapstride_);
392 2 : parse("WRITE_MAP", mapfilename_);
393 2 : if(mapstride_>0 && mapfilename_=="") {
394 0 : error("With WRITE_MAP_STRIDE you must specify WRITE_MAP");
395 : }
396 :
397 : // use GPU?
398 2 : parseFlag("GPU",gpu_);
399 : // set device
400 2 : if (gpu_ && torch::cuda::is_available()) {
401 0 : device_t_ = torch::kCUDA;
402 : } else {
403 2 : device_t_ = torch::kCPU;
404 2 : gpu_ = false;
405 : }
406 :
407 : // Martini model
408 2 : parseFlag("MARTINI",martini_);
409 :
410 : // seed of the random number generator
411 2 : int seed=-1;
412 2 : parse("SEED",seed);
413 2 : if(seed<-1) {
414 0 : error("SEED must be a non-negative integer");
415 : }
416 :
417 : // check read
418 2 : checkRead();
419 :
420 : // set parallel stuff
421 2 : unsigned mpisize=comm.Get_size();
422 2 : if(mpisize>1) {
423 0 : error("EMMIVOX supports only OpenMP parallelization");
424 : }
425 :
426 : // get number of replicas
427 2 : if(no_aver_) {
428 0 : nrep_ = 1;
429 : } else {
430 2 : nrep_ = multi_sim_comm.Get_size();
431 : }
432 2 : replica_ = multi_sim_comm.Get_rank();
433 :
434 2 : if(nrep_>1 && dbfact_>0) {
435 0 : error("Bfactor sampling not supported with ensemble averaging");
436 : }
437 :
438 2 : log.printf(" number of atoms involved : %u\n", atoms.size());
439 2 : log.printf(" experimental density map : %s\n", datafile.c_str());
440 2 : if(no_aver_) {
441 0 : log.printf(" without ensemble averaging\n");
442 : }
443 2 : if(gpu_) {
444 0 : log.printf(" running on GPU \n");
445 : } else {
446 2 : log.printf(" running on CPU \n");
447 : }
448 2 : if(nl_dist_cutoff_ <1.0e+10) {
449 2 : log.printf(" neighbor list distance cutoff : %lf\n", nl_dist_cutoff_);
450 : }
451 2 : if(nl_gauss_cutoff_<1.0e+10) {
452 2 : log.printf(" neighbor list Gaussian sigma cutoff : %lf\n", nl_gauss_cutoff_);
453 : }
454 2 : log.printf(" neighbor list update stride : %u\n", nl_stride_);
455 2 : log.printf(" minimum density error : %f\n", sigma_min);
456 2 : log.printf(" scale factor : %lf\n", scale_);
457 2 : log.printf(" offset : %lf\n", offset_);
458 2 : log.printf(" reading/writing to status file : %s\n", statusfilename_.c_str());
459 2 : log.printf(" with stride : %u\n", statusstride_);
460 2 : if(dbfact_>0) {
461 1 : log.printf(" maximum Bfactor MC move : %f\n", dbfact_);
462 1 : log.printf(" stride MC move : %u\n", MCBstride_);
463 1 : log.printf(" using prior with sigma : %f\n", bfactsig_);
464 : }
465 2 : if(bfactread_) {
466 0 : log.printf(" reading Bfactors from file : %s\n", statusfilename_.c_str());
467 : }
468 2 : log.printf(" temperature of the system in energy unit : %f\n", kbt_);
469 2 : if(nrep_>1) {
470 0 : log.printf(" number of replicas for averaging: %u\n", nrep_);
471 0 : log.printf(" replica ID : %u\n", replica_);
472 : }
473 2 : if(mapstride_>0) {
474 1 : log.printf(" writing model density to file : %s\n", mapfilename_.c_str());
475 1 : log.printf(" with stride : %u\n", mapstride_);
476 : }
477 2 : if(martini_) {
478 0 : log.printf(" using Martini scattering factors\n");
479 : }
480 :
481 : // calculate model constant parameters
482 2 : std::vector<double> Model_w = get_Model_param(atoms);
483 :
484 : // read experimental map and errors
485 2 : get_exp_data(datafile);
486 2 : log.printf(" number of voxels : %u\n", static_cast<unsigned>(ovdd_.size()));
487 :
488 : // normalize atom weight map
489 : double norm_m = accumulate(Model_w.begin(), Model_w.end(), 0.0);
490 : // renormalization and constant factor on atom types
491 130 : for(unsigned i=0; i<Model_w_.size(); ++i) {
492 128 : Vector5d cf;
493 768 : for(unsigned j=0; j<5; ++j) {
494 640 : Model_w_[i][j] *= norm_d / norm_m;
495 640 : cf[j] = Model_w_[i][j]/pow( 2.0*pi, 1.5 );
496 : }
497 128 : cfact_.push_back(cf);
498 : }
499 :
500 : // median density
501 2 : double ovdd_m = get_median(ovdd_);
502 : // median experimental error
503 2 : double err_m = get_median(exp_err_);
504 : // minimum error
505 2 : double minerr = sigma_min*ovdd_m;
506 : // print out statistics
507 2 : log.printf(" median density : %lf\n", ovdd_m);
508 2 : log.printf(" minimum error : %lf\n", minerr);
509 2 : log.printf(" median error : %lf\n", err_m);
510 : // populate ismin: cycle on all voxels
511 14668 : for(unsigned id=0; id<ovdd_.size(); ++id) {
512 : // define smin
513 14666 : double smin = std::max(minerr, exp_err_[id]);
514 : // and to ismin_
515 14666 : ismin_.push_back(1.0/smin);
516 : }
517 :
518 : // prepare gpu stuff: map centers, data, and error
519 2 : prepare_gpu();
520 :
521 : // initialize Bfactors
522 2 : initialize_Bfactor(reso);
523 :
524 : // read status file if restarting
525 2 : if(getRestart() || bfactread_) {
526 0 : read_status();
527 : }
528 :
529 : // prepare auxiliary vectors
530 2 : get_auxiliary_vectors();
531 :
532 : // prepare other vectors: data and derivatives
533 2 : ovmd_.resize(ovdd_.size());
534 2 : atom_der_.resize(Model_type_.size());
535 2 : score_der_.resize(ovdd_.size());
536 :
537 : // add components
538 4 : addComponentWithDerivatives("scoreb");
539 4 : componentIsNotPeriodic("scoreb");
540 4 : addComponent("scale");
541 4 : componentIsNotPeriodic("scale");
542 4 : addComponent("offset");
543 4 : componentIsNotPeriodic("offset");
544 4 : addComponent("kbt");
545 2 : componentIsNotPeriodic("kbt");
546 2 : if(dbfact_>0) {
547 2 : addComponent("accB");
548 2 : componentIsNotPeriodic("accB");
549 : }
550 2 : if(do_corr_) {
551 2 : addComponent("corr");
552 2 : componentIsNotPeriodic("corr");
553 : }
554 :
555 : // initialize random seed
556 : // with SEED replica r gets the stream 1+(SEED*nreplicas+r) mod (2^31-1)
557 : unsigned iseed;
558 2 : if(seed>=0) {
559 1 : const unsigned long long nr = multi_sim_comm.Get_size();
560 1 : iseed = static_cast<unsigned>(1 + (static_cast<unsigned long long>(seed)*nr+replica_) % 2147483646ULL);
561 : } else {
562 1 : iseed = time(NULL)+replica_;
563 : }
564 2 : random_.setSeed(-static_cast<int>(iseed));
565 :
566 : // request atoms
567 2 : requestAtoms(atoms);
568 :
569 : // print bibliography
570 2 : log<<" Bibliography "<<plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)");
571 2 : log<<plumed.cite("Hoff, Thomasen, Lindorff-Larsen, Bonomi, PLoS Comput. Biol. 20 (2024) e1012180");
572 2 : if(!no_aver_ && nrep_>1) {
573 0 : log<<plumed.cite("Bonomi, Camilloni, Cavalli, Vendruscolo, Sci. Adv. 2, e150117 (2016)");
574 : }
575 2 : log<<"\n";
576 2 : }
577 :
578 2 : void EMMIVOX::prepare_gpu() {
579 : // number of data points
580 2 : int nd = ovdd_.size();
581 : // 1) put ismin_ on device_t_
582 8 : ismin_gpu_ = torch::from_blob(ismin_.data(), {nd}, torch::kFloat64).to(torch::kFloat32).to(device_t_);
583 : // 2) put ovdd_ on device_t_
584 6 : ovdd_gpu_ = torch::from_blob(ovdd_.data(), {nd}, torch::kFloat64).to(torch::kFloat32).to(device_t_);
585 : // 3) put Map_m_ on device_t_
586 2 : std::vector<double> Map_m_gpu(3*nd);
587 2 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
588 : for(int i=0; i<nd; ++i) {
589 : Map_m_gpu[i] = Map_m_[i][0];
590 : Map_m_gpu[i+nd] = Map_m_[i][1];
591 : Map_m_gpu[i+2*nd] = Map_m_[i][2];
592 : }
593 : // libtorch tensor
594 8 : Map_m_gpu_ = torch::from_blob(Map_m_gpu.data(), {3,nd}, torch::kFloat64).clone().to(torch::kFloat32).to(device_t_);
595 2 : }
596 :
597 2 : void EMMIVOX::write_model_density(long int step) {
598 2 : OFile ovfile;
599 2 : ovfile.link(*this);
600 : std::string num;
601 2 : Tools::convert(step,num);
602 2 : std::string name = mapfilename_+"-"+num;
603 2 : ovfile.open(name);
604 : ovfile.setHeavyFlush();
605 2 : ovfile.fmtField("%10.7e ");
606 : // write density
607 14668 : for(unsigned i=0; i<ovmd_.size(); ++i) {
608 14666 : ovfile.printField("Model", ovmd_[i]);
609 29332 : ovfile.printField("ModelScaled", scale_ * ovmd_[i] + offset_);
610 14666 : ovfile.printField("Data", ovdd_[i]);
611 14666 : ovfile.printField();
612 : }
613 2 : ovfile.close();
614 2 : }
615 :
616 4 : double EMMIVOX::get_median(std::vector<double> v) {
617 : // dimension of vector
618 4 : unsigned size = v.size();
619 : // in case of only one entry
620 4 : if (size==1) {
621 0 : return v[0];
622 : } else {
623 : // reorder vector
624 4 : sort(v.begin(), v.end());
625 : // odd or even?
626 4 : if (size%2==0) {
627 0 : return (v[size/2-1]+v[size/2])/2.0;
628 : } else {
629 4 : return v[size/2];
630 : }
631 : }
632 : }
633 :
634 0 : void EMMIVOX::read_status() {
635 : double MDtime;
636 : // open file
637 0 : IFile *ifile = new IFile();
638 0 : ifile->link(*this);
639 0 : if(ifile->FileExist(statusfilename_)) {
640 0 : ifile->open(statusfilename_);
641 0 : while(ifile->scanField("MD_time", MDtime)) {
642 : // read scale and offset
643 0 : ifile->scanField("scale", scale_);
644 0 : ifile->scanField("offset", offset_);
645 : // read bfactors if doing fitting of reading it at restart
646 0 : if(dbfact_>0 || bfactread_) {
647 : // cycle on residues
648 0 : for(unsigned ir=0; ir<Model_rlist_.size(); ++ir) {
649 : // key: pair of residue/chain IDs
650 : std::pair<int,std::string> key = Model_rlist_[ir];
651 : // convert ires to std::string
652 : std::string num;
653 0 : Tools::convert(key.first,num);
654 : // read entry
655 0 : std::string ch = key.second;
656 0 : if(ch==" ") {
657 : ch="";
658 : }
659 0 : ifile->scanField("bf-"+num+":"+ch, Model_b_[key]);
660 : }
661 : }
662 : // new line
663 0 : ifile->scanField();
664 : }
665 0 : ifile->close();
666 : } else {
667 0 : error("Cannot find status file "+statusfilename_+"\n");
668 : }
669 0 : delete ifile;
670 0 : }
671 :
672 6 : void EMMIVOX::print_status(long int step) {
673 : // if first time open the file
674 6 : if(first_status_) {
675 2 : first_status_ = false;
676 2 : statusfile_.link(*this);
677 2 : statusfile_.open(statusfilename_);
678 : statusfile_.setHeavyFlush();
679 4 : statusfile_.fmtField("%10.7e ");
680 : }
681 : // write fields
682 6 : double MDtime = static_cast<double>(step)*getTimeStep();
683 6 : statusfile_.printField("MD_time", MDtime);
684 : // write scale and offset
685 6 : statusfile_.printField("scale", scale_);
686 6 : statusfile_.printField("offset", offset_);
687 : // write bfactors only if doing fitting or reading bfactors
688 6 : if(dbfact_>0 || bfactread_) {
689 : // cycle on residues
690 93 : for(unsigned ir=0; ir<Model_rlist_.size(); ++ir) {
691 : // key: pair of residue/chain IDs
692 : std::pair<int,std::string> key = Model_rlist_[ir];
693 : // bfactor from map
694 90 : double bf = Model_b_[key];
695 : // convert ires to std::string
696 : std::string num;
697 90 : Tools::convert(key.first,num);
698 : // print entry
699 270 : statusfile_.printField("bf-"+num+":"+key.second, bf);
700 : }
701 : }
702 6 : statusfile_.printField();
703 6 : }
704 :
705 60 : bool EMMIVOX::doAccept(double oldE, double newE, double kbt) {
706 : bool accept = false;
707 : // calculate delta energy
708 60 : double delta = ( newE - oldE ) / kbt;
709 : // if delta is negative always accept move
710 60 : if( delta < 0.0 ) {
711 : accept = true;
712 : } else {
713 : // otherwise extract random number
714 34 : double s = random_.RandU01();
715 34 : if( s < exp(-delta) ) {
716 : accept = true;
717 : }
718 : }
719 60 : return accept;
720 : }
721 :
722 2 : std::vector<double> EMMIVOX::get_Model_param(std::vector<AtomNumber> &atoms) {
723 : // check if MOLINFO line is present
724 2 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
725 2 : if(!moldat) {
726 0 : error("MOLINFO DATA not found\n");
727 : }
728 2 : log<<" MOLINFO DATA found with label " <<moldat->getLabel()<<", using proper atom names\n";
729 :
730 : // list of weights - one per atom
731 2 : std::vector<double> Model_w;
732 : // 5-Gaussians parameters
733 : // map of atom types to A and B coefficients of scattering factor
734 : // f(s) = A * exp(-B*s**2)
735 : // B is in Angstrom squared
736 : // Elastic atomic scattering factors of electrons for neutral atoms
737 : // and s up to 6.0 A^-1: as implemented in PLUMED
738 : // map between an atom type and an index
739 : std::map<std::string, unsigned> type_map;
740 : // atomistic types
741 2 : type_map["C"]=0;
742 2 : type_map["O"]=1;
743 2 : type_map["N"]=2;
744 2 : type_map["S"]=3;
745 2 : type_map["P"]=4;
746 2 : type_map["F"]=5;
747 2 : type_map["NA"]=6;
748 2 : type_map["MG"]=7;
749 2 : type_map["CL"]=8;
750 2 : type_map["CA"]=9;
751 2 : type_map["K"]=10;
752 2 : type_map["ZN"]=11;
753 : // Martini types
754 2 : type_map["ALA_BB"]=12;
755 2 : type_map["ALA_SC1"]=13;
756 2 : type_map["CYS_BB"]=14;
757 2 : type_map["CYS_SC1"]=15;
758 2 : type_map["ASP_BB"]=16;
759 2 : type_map["ASP_SC1"]=17;
760 2 : type_map["GLU_BB"]=18;
761 2 : type_map["GLU_SC1"]=19;
762 2 : type_map["PHE_BB"]=20;
763 2 : type_map["PHE_SC1"]=21;
764 2 : type_map["PHE_SC2"]=22;
765 2 : type_map["PHE_SC3"]=23;
766 2 : type_map["GLY_BB"]=24;
767 2 : type_map["HIS_BB"]=25;
768 2 : type_map["HIS_SC1"]=26;
769 2 : type_map["HIS_SC2"]=27;
770 2 : type_map["HIS_SC3"]=28;
771 2 : type_map["ILE_BB"]=29;
772 2 : type_map["ILE_SC1"]=30;
773 2 : type_map["LYS_BB"]=31;
774 2 : type_map["LYS_SC1"]=32;
775 2 : type_map["LYS_SC2"]=33;
776 2 : type_map["LEU_BB"]=34;
777 2 : type_map["LEU_SC1"]=35;
778 2 : type_map["MET_BB"]=36;
779 2 : type_map["MET_SC1"]=37;
780 2 : type_map["ASN_BB"]=38;
781 2 : type_map["ASN_SC1"]=39;
782 2 : type_map["PRO_BB"]=40;
783 2 : type_map["PRO_SC1"]=41;
784 2 : type_map["GLN_BB"]=42;
785 2 : type_map["GLN_SC1"]=43;
786 2 : type_map["ARG_BB"]=44;
787 2 : type_map["ARG_SC1"]=45;
788 2 : type_map["ARG_SC2"]=46;
789 2 : type_map["SER_BB"]=47;
790 2 : type_map["SER_SC1"]=48;
791 2 : type_map["THR_BB"]=49;
792 2 : type_map["THR_SC1"]=50;
793 2 : type_map["VAL_BB"]=51;
794 2 : type_map["VAL_SC1"]=52;
795 2 : type_map["TRP_BB"]=53;
796 2 : type_map["TRP_SC1"]=54;
797 2 : type_map["TRP_SC2"]=55;
798 2 : type_map["TRP_SC3"]=56;
799 2 : type_map["TRP_SC4"]=57;
800 2 : type_map["TRP_SC5"]=58;
801 2 : type_map["TYR_BB"]=59;
802 2 : type_map["TYR_SC1"]=60;
803 2 : type_map["TYR_SC2"]=61;
804 2 : type_map["TYR_SC3"]=62;
805 2 : type_map["TYR_SC4"]=63;
806 : // fill in sigma vector for atoms
807 4 : Model_s_.push_back(0.01*Vector5d(0.114,1.0825,5.4281,17.8811,51.1341)); // C
808 4 : Model_s_.push_back(0.01*Vector5d(0.0652,0.6184,2.9449,9.6298,28.2194)); // O
809 4 : Model_s_.push_back(0.01*Vector5d(0.0541,0.5165,2.8207,10.6297,34.3764)); // N
810 4 : Model_s_.push_back(0.01*Vector5d(0.0838,0.7788,4.3462,15.5846,44.63655)); // S
811 4 : Model_s_.push_back(0.01*Vector5d(0.0977,0.9084,4.9654,18.5471,54.3648)); // P
812 4 : Model_s_.push_back(0.01*Vector5d(0.0613,0.5753,2.6858,8.8214,25.6668)); // F
813 4 : Model_s_.push_back(0.01*Vector5d(0.1684,1.7150,8.8386,50.8265,147.2073)); // NA
814 4 : Model_s_.push_back(0.01*Vector5d(0.1356,1.3579,6.9255,32.3165,92.1138)); // MG
815 4 : Model_s_.push_back(0.01*Vector5d(0.0694,0.6443,3.5351,12.5058,35.8633)); // CL
816 4 : Model_s_.push_back(0.01*Vector5d(0.1742,1.8329,8.8407,47.4583,134.9613)); // CA
817 4 : Model_s_.push_back(0.01*Vector5d(0.1660,1.6906,8.7447,46.7825,165.6923)); // K
818 4 : Model_s_.push_back(0.01*Vector5d(0.0876,0.8650,3.8612,18.8726,64.7016)); // ZN
819 : // fill in sigma vector for Martini beads
820 4 : Model_s_.push_back(0.01*Vector5d(22.000000,22.000000,22.000000,22.000000,22.000000)); // ALA_BB
821 4 : Model_s_.push_back(0.01*Vector5d(0.500000,0.500000,0.500000,0.500000,0.500000)); // ALA_SC1
822 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // CYS_BB
823 4 : Model_s_.push_back(0.01*Vector5d(8.500000,8.500000,8.500000,8.500000,8.500000)); // CYS_SC1
824 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // ASP_BB
825 4 : Model_s_.push_back(0.01*Vector5d(17.000000,17.000000,17.000000,17.000000,17.000000)); // ASP_SC1
826 4 : Model_s_.push_back(0.01*Vector5d(22.000000,22.000000,22.000000,22.000000,22.000000)); // GLU_BB
827 4 : Model_s_.push_back(0.01*Vector5d(24.000000,24.000000,24.000000,24.000000,24.000000)); // GLU_SC1
828 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // PHE_BB
829 4 : Model_s_.push_back(0.01*Vector5d(17.500000,17.500000,17.500000,17.500000,17.500000)); // PHE_SC1
830 4 : Model_s_.push_back(0.01*Vector5d(11.000000,11.000000,11.000000,11.000000,11.000000)); // PHE_SC2
831 4 : Model_s_.push_back(0.01*Vector5d(11.000000,11.000000,11.000000,11.000000,11.000000)); // PHE_SC3
832 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // GLY_BB
833 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // HIS_BB
834 4 : Model_s_.push_back(0.01*Vector5d(11.500000,11.500000,11.500000,11.500000,11.500000)); // HIS_SC1
835 4 : Model_s_.push_back(0.01*Vector5d(9.000000,9.000000,9.000000,9.000000,9.000000)); // HIS_SC2
836 4 : Model_s_.push_back(0.01*Vector5d(8.500000,8.500000,8.500000,8.500000,8.500000)); // HIS_SC3
837 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // ILE_BB
838 4 : Model_s_.push_back(0.01*Vector5d(25.500000,25.500000,25.500000,25.500000,25.500000)); // ILE_SC1
839 4 : Model_s_.push_back(0.01*Vector5d(22.000000,22.000000,22.000000,22.000000,22.000000)); // LYS_BB
840 4 : Model_s_.push_back(0.01*Vector5d(18.000000,18.000000,18.000000,18.000000,18.000000)); // LYS_SC1
841 4 : Model_s_.push_back(0.01*Vector5d(11.000000,11.000000,11.000000,11.000000,11.000000)); // LYS_SC2
842 4 : Model_s_.push_back(0.01*Vector5d(22.000000,22.000000,22.000000,22.000000,22.000000)); // LEU_BB
843 4 : Model_s_.push_back(0.01*Vector5d(21.500000,21.500000,21.500000,21.500000,21.500000)); // LEU_SC1
844 4 : Model_s_.push_back(0.01*Vector5d(22.000000,22.000000,22.000000,22.000000,22.000000)); // MET_BB
845 4 : Model_s_.push_back(0.01*Vector5d(22.500000,22.500000,22.500000,22.500000,22.500000)); // MET_SC1
846 4 : Model_s_.push_back(0.01*Vector5d(22.000000,22.000000,22.000000,22.000000,22.000000)); // ASN_BB
847 4 : Model_s_.push_back(0.01*Vector5d(18.500000,18.500000,18.500000,18.500000,18.500000)); // ASN_SC1
848 4 : Model_s_.push_back(0.01*Vector5d(23.500000,23.500000,23.500000,23.500000,23.500000)); // PRO_BB
849 4 : Model_s_.push_back(0.01*Vector5d(17.500000,17.500000,17.500000,17.500000,17.500000)); // PRO_SC1
850 4 : Model_s_.push_back(0.01*Vector5d(22.000000,22.000000,22.000000,22.000000,22.000000)); // GLN_BB
851 4 : Model_s_.push_back(0.01*Vector5d(24.500000,24.500000,24.500000,24.500000,24.500000)); // GLN_SC1
852 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // ARG_BB
853 4 : Model_s_.push_back(0.01*Vector5d(18.000000,18.000000,18.000000,18.000000,18.000000)); // ARG_SC1
854 4 : Model_s_.push_back(0.01*Vector5d(18.000000,18.000000,18.000000,18.000000,18.000000)); // ARG_SC2
855 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // SER_BB
856 4 : Model_s_.push_back(0.01*Vector5d(9.000000,9.000000,9.000000,9.000000,9.000000)); // SER_SC1
857 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // THR_BB
858 4 : Model_s_.push_back(0.01*Vector5d(17.000000,17.000000,17.000000,17.000000,17.000000)); // THR_SC1
859 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // VAL_BB
860 4 : Model_s_.push_back(0.01*Vector5d(18.000000,18.000000,18.000000,18.000000,18.000000)); // VAL_SC1
861 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // TRP_BB
862 4 : Model_s_.push_back(0.01*Vector5d(11.500000,11.500000,11.500000,11.500000,11.500000)); // TRP_SC1
863 4 : Model_s_.push_back(0.01*Vector5d(9.000000,9.000000,9.000000,9.000000,9.000000)); // TRP_SC2
864 4 : Model_s_.push_back(0.01*Vector5d(11.000000,11.000000,11.000000,11.000000,11.000000)); // TRP_SC3
865 4 : Model_s_.push_back(0.01*Vector5d(11.000000,11.000000,11.000000,11.000000,11.000000)); // TRP_SC4
866 4 : Model_s_.push_back(0.01*Vector5d(9.500000,9.500000,9.500000,9.500000,9.500000)); // TRP_SC5
867 4 : Model_s_.push_back(0.01*Vector5d(23.000000,23.000000,23.000000,23.000000,23.000000)); // TYR_BB
868 4 : Model_s_.push_back(0.01*Vector5d(12.000000,12.000000,12.000000,12.000000,12.000000)); // TYR_SC1
869 4 : Model_s_.push_back(0.01*Vector5d(11.000000,11.000000,11.000000,11.000000,11.000000)); // TYR_SC2
870 4 : Model_s_.push_back(0.01*Vector5d(11.000000,11.000000,11.000000,11.000000,11.000000)); // TYR_SC3
871 4 : Model_s_.push_back(0.01*Vector5d(8.500000,8.500000,8.500000,8.500000,8.500000)); // TYR_SC4
872 : // fill in weight vector for atoms
873 2 : Model_w_.push_back(Vector5d(0.0489,0.2091,0.7537,1.1420,0.3555)); // C
874 2 : Model_w_.push_back(Vector5d(0.0365,0.1729,0.5805,0.8814,0.3121)); // O
875 2 : Model_w_.push_back(Vector5d(0.0267,0.1328,0.5301,1.1020,0.4215)); // N
876 2 : Model_w_.push_back(Vector5d(0.0915,0.4312,1.0847,2.4671,1.0852)); // S
877 2 : Model_w_.push_back(Vector5d(0.1005,0.4615,1.0663,2.5854,1.2725)); // P
878 2 : Model_w_.push_back(Vector5d(0.0382,0.1822,0.5972,0.7707,0.2130)); // F
879 2 : Model_w_.push_back(Vector5d(0.1260,0.6442,0.8893,1.8197,1.2988)); // NA
880 2 : Model_w_.push_back(Vector5d(0.1130,0.5575,0.9046,2.1580,1.4735)); // MG
881 2 : Model_w_.push_back(Vector5d(0.0799,0.3891,1.0037,2.3332,1.0507)); // CL
882 2 : Model_w_.push_back(Vector5d(0.2355,0.9916,2.3959,3.7252,2.5647)); // CA
883 2 : Model_w_.push_back(Vector5d(0.2149,0.8703,2.4999,2.3591,3.0318)); // K
884 2 : Model_w_.push_back(Vector5d(0.1780,0.8096,1.6744,1.9499,1.4495)); // ZN
885 : // fill in weight vector for Martini beads
886 2 : Model_w_.push_back(Vector5d(1.800000,1.800000,1.800000,1.800000,1.800000)); // ALA_BB
887 2 : Model_w_.push_back(Vector5d(0.100000,0.100000,0.100000,0.100000,0.100000)); // ALA_SC1
888 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // CYS_BB
889 2 : Model_w_.push_back(Vector5d(1.100000,1.100000,1.100000,1.100000,1.100000)); // CYS_SC1
890 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // ASP_BB
891 2 : Model_w_.push_back(Vector5d(1.700000,1.700000,1.700000,1.700000,1.700000)); // ASP_SC1
892 2 : Model_w_.push_back(Vector5d(1.800000,1.800000,1.800000,1.800000,1.800000)); // GLU_BB
893 2 : Model_w_.push_back(Vector5d(2.300000,2.300000,2.300000,2.300000,2.300000)); // GLU_SC1
894 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // PHE_BB
895 2 : Model_w_.push_back(Vector5d(1.400000,1.400000,1.400000,1.400000,1.400000)); // PHE_SC1
896 2 : Model_w_.push_back(Vector5d(0.900000,0.900000,0.900000,0.900000,0.900000)); // PHE_SC2
897 2 : Model_w_.push_back(Vector5d(0.900000,0.900000,0.900000,0.900000,0.900000)); // PHE_SC3
898 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // GLY_BB
899 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // HIS_BB
900 2 : Model_w_.push_back(Vector5d(0.900000,0.900000,0.900000,0.900000,0.900000)); // HIS_SC1
901 2 : Model_w_.push_back(Vector5d(0.800000,0.800000,0.800000,0.800000,0.800000)); // HIS_SC2
902 2 : Model_w_.push_back(Vector5d(0.800000,0.800000,0.800000,0.800000,0.800000)); // HIS_SC3
903 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // ILE_BB
904 2 : Model_w_.push_back(Vector5d(2.000000,2.000000,2.000000,2.000000,2.000000)); // ILE_SC1
905 2 : Model_w_.push_back(Vector5d(1.800000,1.800000,1.800000,1.800000,1.800000)); // LYS_BB
906 2 : Model_w_.push_back(Vector5d(1.400000,1.400000,1.400000,1.400000,1.400000)); // LYS_SC1
907 2 : Model_w_.push_back(Vector5d(0.900000,0.900000,0.900000,0.900000,0.900000)); // LYS_SC2
908 2 : Model_w_.push_back(Vector5d(1.800000,1.800000,1.800000,1.800000,1.800000)); // LEU_BB
909 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // LEU_SC1
910 2 : Model_w_.push_back(Vector5d(1.800000,1.800000,1.800000,1.800000,1.800000)); // MET_BB
911 2 : Model_w_.push_back(Vector5d(2.300000,2.300000,2.300000,2.300000,2.300000)); // MET_SC1
912 2 : Model_w_.push_back(Vector5d(1.800000,1.800000,1.800000,1.800000,1.800000)); // ASN_BB
913 2 : Model_w_.push_back(Vector5d(1.800000,1.800000,1.800000,1.800000,1.800000)); // ASN_SC1
914 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // PRO_BB
915 2 : Model_w_.push_back(Vector5d(1.400000,1.400000,1.400000,1.400000,1.400000)); // PRO_SC1
916 2 : Model_w_.push_back(Vector5d(1.800000,1.800000,1.800000,1.800000,1.800000)); // GLN_BB
917 2 : Model_w_.push_back(Vector5d(2.300000,2.300000,2.300000,2.300000,2.300000)); // GLN_SC1
918 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // ARG_BB
919 2 : Model_w_.push_back(Vector5d(1.400000,1.400000,1.400000,1.400000,1.400000)); // ARG_SC1
920 2 : Model_w_.push_back(Vector5d(1.800000,1.800000,1.800000,1.800000,1.800000)); // ARG_SC2
921 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // SER_BB
922 2 : Model_w_.push_back(Vector5d(0.800000,0.800000,0.800000,0.800000,0.800000)); // SER_SC1
923 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // THR_BB
924 2 : Model_w_.push_back(Vector5d(1.400000,1.400000,1.400000,1.400000,1.400000)); // THR_SC1
925 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // VAL_BB
926 2 : Model_w_.push_back(Vector5d(1.400000,1.400000,1.400000,1.400000,1.400000)); // VAL_SC1
927 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // TRP_BB
928 2 : Model_w_.push_back(Vector5d(0.900000,0.900000,0.900000,0.900000,0.900000)); // TRP_SC1
929 2 : Model_w_.push_back(Vector5d(0.800000,0.800000,0.800000,0.800000,0.800000)); // TRP_SC2
930 2 : Model_w_.push_back(Vector5d(0.900000,0.900000,0.900000,0.900000,0.900000)); // TRP_SC3
931 2 : Model_w_.push_back(Vector5d(0.900000,0.900000,0.900000,0.900000,0.900000)); // TRP_SC4
932 2 : Model_w_.push_back(Vector5d(0.800000,0.800000,0.800000,0.800000,0.800000)); // TRP_SC5
933 2 : Model_w_.push_back(Vector5d(1.900000,1.900000,1.900000,1.900000,1.900000)); // TYR_BB
934 2 : Model_w_.push_back(Vector5d(0.900000,0.900000,0.900000,0.900000,0.900000)); // TYR_SC1
935 2 : Model_w_.push_back(Vector5d(0.900000,0.900000,0.900000,0.900000,0.900000)); // TYR_SC2
936 2 : Model_w_.push_back(Vector5d(0.900000,0.900000,0.900000,0.900000,0.900000)); // TYR_SC3
937 2 : Model_w_.push_back(Vector5d(0.800000,0.800000,0.800000,0.800000,0.800000)); // TYR_SC4
938 : // cycle on atoms
939 470 : for(unsigned i=0; i<atoms.size(); ++i) {
940 : // get atom name
941 468 : std::string name = moldat->getAtomName(atoms[i]);
942 : // get residue name
943 468 : std::string resname = moldat->getResidueName(atoms[i]);
944 : // type of atoms/bead
945 : std::string type_s;
946 : // Martini model
947 468 : if(martini_) {
948 0 : type_s = resname+"_"+name;
949 : // Atomistic model
950 : } else {
951 : char type;
952 : // get atom type
953 468 : char first = name.at(0);
954 : // GOLDEN RULE: type is first letter, if not a number
955 468 : if (!isdigit(first)) {
956 : type = first;
957 : // otherwise is the second
958 : } else {
959 0 : type = name.at(1);
960 : }
961 : // convert to std::string
962 468 : type_s = std::string(1,type);
963 : // special cases
964 468 : if(name=="SOD" || name=="NA" || name =="Na") {
965 : type_s = "NA";
966 : }
967 468 : if(name=="MG" || name=="Mg") {
968 : type_s = "MG";
969 : }
970 468 : if(name=="CLA" || name=="CL" || name =="Cl") {
971 : type_s = "CL";
972 : }
973 468 : if((resname=="CAL" || resname=="CA") && (name=="CAL" || name=="CA" || name =="C0")) {
974 : type_s = "CA";
975 : }
976 468 : if(name=="POT" || name=="K") {
977 : type_s = "K";
978 : }
979 468 : if(name=="ZN" || name=="Zn") {
980 : type_s = "ZN";
981 : }
982 : }
983 : // check if key in map
984 468 : if(type_map.find(type_s) != type_map.end()) {
985 : // save atom type
986 468 : Model_type_.push_back(type_map[type_s]);
987 : // this will be normalized in the final density
988 468 : Vector5d w = Model_w_[type_map[type_s]];
989 468 : Model_w.push_back(w[0]+w[1]+w[2]+w[3]+w[4]);
990 : // get residue id
991 468 : const int ires = moldat->getResidueNumber(atoms[i]);
992 : // and chain
993 468 : std::string c ("*");
994 468 : if(!bfactnoc_) {
995 936 : c = moldat->getChainID(atoms[i]);
996 : }
997 : // define pair residue/chain IDs
998 : std::pair<int,std::string> key = std::make_pair(ires,c);
999 : // add to map between residue/chain and list of atoms
1000 468 : Model_resmap_[key].push_back(i);
1001 : // and global list of residue/chain per atom
1002 468 : Model_res_.push_back(key);
1003 : // initialize Bfactor map
1004 468 : Model_b_[key] = 0.0;
1005 : } else {
1006 0 : error("Wrong atom type "+type_s+" from atom name "+name+"\n");
1007 : }
1008 : }
1009 : // create ordered vector of residue-chain IDs
1010 470 : for(unsigned i=0; i<Model_res_.size(); ++i) {
1011 : std::pair<int,std::string> key = Model_res_[i];
1012 : // search in Model_rlist_
1013 468 : if(find(Model_rlist_.begin(), Model_rlist_.end(), key) == Model_rlist_.end()) {
1014 60 : Model_rlist_.push_back(key);
1015 : }
1016 : }
1017 : // return weights
1018 2 : return Model_w;
1019 0 : }
1020 :
1021 : // read experimental data file in PLUMED format:
1022 2 : void EMMIVOX::get_exp_data(const std::string &datafile) {
1023 2 : Vector pos;
1024 : double dens, err;
1025 : int idcomp;
1026 :
1027 : // open file
1028 2 : IFile *ifile = new IFile();
1029 2 : if(ifile->FileExist(datafile)) {
1030 2 : ifile->open(datafile);
1031 29336 : while(ifile->scanField("Id",idcomp)) {
1032 14666 : ifile->scanField("Pos_0",pos[0]);
1033 14666 : ifile->scanField("Pos_1",pos[1]);
1034 14666 : ifile->scanField("Pos_2",pos[2]);
1035 14666 : ifile->scanField("Density",dens);
1036 14666 : ifile->scanField("Error",err);
1037 : // voxel center
1038 14666 : Map_m_.push_back(pos);
1039 : // experimental density
1040 14666 : ovdd_.push_back(dens);
1041 : // error
1042 14666 : exp_err_.push_back(err);
1043 : // new line
1044 14666 : ifile->scanField();
1045 : }
1046 2 : ifile->close();
1047 : } else {
1048 0 : error("Cannot find DATA_FILE "+datafile+"\n");
1049 : }
1050 2 : delete ifile;
1051 2 : }
1052 :
1053 2 : void EMMIVOX::initialize_Bfactor(double reso) {
1054 2 : double bfactini = 0.0;
1055 : // if doing Bfactor Monte Carlo
1056 2 : if(dbfact_>0) {
1057 : // initialize B factor based on empirical relation between resolution and average bfactor
1058 : // calculated on ~8000 cryo-EM data with resolution < 5 Ang
1059 : // Bfact = A*reso**2+B; with A=6.95408 B=-2.45697/100.0 nm^2
1060 1 : bfactini = 6.95408*reso*reso - 0.01*2.45697;
1061 : // check for min and max
1062 3 : bfactini = std::min(bfactmax_, std::max(bfactmin_, bfactini));
1063 : }
1064 : // set initial Bfactor
1065 62 : for(std::map< std::pair<int,std::string>, double>::iterator it=Model_b_.begin(); it!=Model_b_.end(); ++it) {
1066 60 : it->second = bfactini;
1067 : }
1068 2 : log.printf(" experimental map resolution : %3.2f\n", reso);
1069 : // if doing Bfactor Monte Carlo
1070 2 : if(dbfact_>0) {
1071 1 : log.printf(" minimum Bfactor value : %3.2f\n", bfactmin_);
1072 1 : log.printf(" maximum Bfactor value : %3.2f\n", bfactmax_);
1073 1 : log.printf(" initial Bfactor value : %3.2f\n", bfactini);
1074 : }
1075 2 : }
1076 :
1077 : // prepare auxiliary vectors
1078 4 : void EMMIVOX::get_auxiliary_vectors() {
1079 : // number of atoms
1080 4 : unsigned natoms = Model_res_.size();
1081 : // clear lists
1082 4 : pref_.clear();
1083 4 : invs2_.clear();
1084 4 : cut_.clear();
1085 : // resize
1086 4 : pref_.resize(natoms);
1087 4 : invs2_.resize(natoms);
1088 4 : cut_.resize(natoms);
1089 : // cycle on all atoms
1090 4 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1091 : for(unsigned im=0; im<natoms; ++im) {
1092 : // get atom type
1093 : unsigned atype = Model_type_[im];
1094 : // get residue/chain IDs
1095 : std::pair<int,std::string> key = Model_res_[im];
1096 : // get bfactor
1097 : double bfact = Model_b_[key];
1098 : // sigma for 5 gaussians
1099 : Vector5d m_s = Model_s_[atype];
1100 : // calculate constant quantities
1101 : Vector5d pref, invs2;
1102 : // calculate cutoff
1103 : double n = 0.0;
1104 : double d = 0.0;
1105 : for(unsigned j=0; j<5; ++j) {
1106 : // total value of b
1107 : double m_b = m_s[j] + bfact/4.0;
1108 : // calculate invs2
1109 : invs2[j] = 1.0/(inv_pi2_*m_b);
1110 : // prefactor
1111 : pref[j] = cfact_[atype][j] * pow(invs2[j],1.5);
1112 : // cutoff
1113 : n += pref[j] / invs2[j];
1114 : d += pref[j];
1115 : }
1116 : // put into global lists
1117 : pref_[im] = pref;
1118 : invs2_[im] = invs2;
1119 : cut_[im] = std::min(nl_dist_cutoff_, sqrt(n/d)*nl_gauss_cutoff_);
1120 : }
1121 : // push to GPU
1122 4 : push_auxiliary_gpu();
1123 4 : }
1124 :
1125 4 : void EMMIVOX::push_auxiliary_gpu() {
1126 : // 1) create vector of pref_ and invs2_
1127 4 : int natoms = Model_type_.size();
1128 4 : std::vector<double> pref(5*natoms), invs2(5*natoms);
1129 4 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1130 : for(int i=0; i<natoms; ++i) {
1131 : for(int j=0; j<5; ++j) {
1132 : pref[i+j*natoms] = pref_[i][j];
1133 : invs2[i+j*natoms] = invs2_[i][j];
1134 : }
1135 : }
1136 : // 2) initialize gpu tensors
1137 16 : pref_gpu_ = torch::from_blob(pref.data(), {5,natoms}, torch::kFloat64).clone().to(torch::kFloat32).to(device_t_);
1138 16 : invs2_gpu_ = torch::from_blob(invs2.data(), {5,natoms}, torch::kFloat64).clone().to(torch::kFloat32).to(device_t_);
1139 4 : }
1140 :
1141 2 : void EMMIVOX::get_close_residues() {
1142 : // two residues are neighbors if two of their atoms are closer than this
1143 : const double dcut = 0.5;
1144 2 : const unsigned nres = Model_rlist_.size();
1145 : // clear neighbor list
1146 2 : nl_res_.clear();
1147 2 : nl_res_.resize(nres);
1148 :
1149 : // atoms, center and radius of each residue: two residues whose spheres are
1150 : // farther than dcut cannot have atoms closer than dcut
1151 2 : std::vector< const std::vector<unsigned>* > ratoms(nres);
1152 2 : std::vector<Vector> rcenter(nres);
1153 2 : std::vector<double> rradius(nres,0.0);
1154 62 : for(unsigned i=0; i<nres; ++i) {
1155 60 : ratoms[i] = &Model_resmap_[Model_rlist_[i]];
1156 : Vector c;
1157 528 : for(const auto im : *ratoms[i]) {
1158 468 : c += getPosition(im);
1159 : }
1160 60 : if(ratoms[i]->size()>0) {
1161 60 : c /= static_cast<double>(ratoms[i]->size());
1162 : }
1163 60 : rcenter[i] = c;
1164 528 : for(const auto im : *ratoms[i]) {
1165 634 : rradius[i] = std::max(rradius[i], delta(c,getPosition(im)).modulo());
1166 : }
1167 : }
1168 :
1169 : // loop in parallel
1170 2 : const unsigned nt = OpenMP::getNumThreads();
1171 2 : std::vector< std::vector< std::vector<unsigned> > > nl_res_t(nt);
1172 2 : #pragma omp parallel num_threads(nt)
1173 : {
1174 : // private variable
1175 : std::vector< std::vector<unsigned> > & nl_res_l = nl_res_t[OpenMP::getThreadNum()];
1176 : nl_res_l.resize(nres);
1177 : // cycle on residues/chains #1
1178 : #pragma omp for schedule(dynamic)
1179 : for(unsigned i=0; i<nres; ++i) {
1180 : // cycle over residues/chains #2
1181 : for(unsigned j=i+1; j<nres; ++j) {
1182 : if(delta(rcenter[i],rcenter[j]).modulo() > rradius[i]+rradius[j]+dcut) {
1183 : continue;
1184 : }
1185 : // set flag neighbor
1186 : bool neigh = false;
1187 : // cycle over all the atoms belonging to residue i
1188 : for(const auto im1 : *ratoms[i]) {
1189 : // get atom position #1
1190 : const Vector pos1 = getPosition(im1);
1191 : // cycle over all the atoms belonging to residue j
1192 : for(const auto im2 : *ratoms[j]) {
1193 : // if closer than 0.5 nm, then residues i and j are neighbors
1194 : if(delta(pos1,getPosition(im2)).modulo()<dcut) {
1195 : neigh = true;
1196 : break;
1197 : }
1198 : }
1199 : // check if neighbor already found
1200 : if(neigh) {
1201 : break;
1202 : }
1203 : }
1204 : // if neighbors, add to local list
1205 : if(neigh) {
1206 : nl_res_l[i].push_back(j);
1207 : nl_res_l[j].push_back(i);
1208 : }
1209 : }
1210 : }
1211 : }
1212 : // add to global list, sorted (reproducible results); the lists of the threads
1213 : // that did not run (e.g. without OpenMP) are empty and are skipped
1214 6 : for(unsigned t=0; t<nt; ++t) {
1215 4 : if(nl_res_t[t].empty()) {
1216 0 : continue;
1217 : }
1218 124 : for(unsigned i=0; i<nres; ++i) {
1219 120 : nl_res_[i].insert(nl_res_[i].end(), nl_res_t[t][i].begin(), nl_res_t[t][i].end());
1220 : }
1221 : }
1222 62 : for(unsigned i=0; i<nres; ++i) {
1223 60 : std::sort(nl_res_[i].begin(), nl_res_[i].end());
1224 : }
1225 6 : }
1226 :
1227 2 : void EMMIVOX::doMonteCarloBfact() {
1228 : // update residue neighbor list
1229 2 : get_close_residues();
1230 :
1231 : // cycle over residues/chains
1232 62 : for(unsigned ir=0; ir<Model_rlist_.size(); ++ir) {
1233 :
1234 : // key: pair of residue/chain IDs
1235 : std::pair<int,std::string> key = Model_rlist_[ir];
1236 : // old bfactor
1237 60 : double bfactold = Model_b_[key];
1238 :
1239 : // propose move in bfactor
1240 60 : double bfactnew = bfactold + dbfact_ * ( 2.0 * random_.RandU01() - 1.0 );
1241 : // check boundaries
1242 60 : if(bfactnew > bfactmax_) {
1243 0 : bfactnew = 2.0*bfactmax_ - bfactnew;
1244 : }
1245 60 : if(bfactnew < bfactmin_) {
1246 0 : bfactnew = 2.0*bfactmin_ - bfactnew;
1247 : }
1248 :
1249 : // useful quantities
1250 : std::map<unsigned, double> deltaov;
1251 60 : const unsigned nt = OpenMP::getNumThreads();
1252 60 : std::vector< std::map<unsigned, double> > deltaov_t(nt);
1253 :
1254 60 : #pragma omp parallel num_threads(nt)
1255 : {
1256 : // private variables
1257 : std::map<unsigned, double> & deltaov_l = deltaov_t[OpenMP::getThreadNum()];
1258 : #pragma omp for schedule(static)
1259 : // cycle over all the atoms belonging to key (residue/chain)
1260 : for(unsigned ia=0; ia<Model_resmap_[key].size(); ++ia) {
1261 :
1262 : // get atom id
1263 : unsigned im = Model_resmap_[key][ia];
1264 : // get atom type
1265 : unsigned atype = Model_type_[im];
1266 : // sigma for 5 Gaussians
1267 : Vector5d m_s = Model_s_[atype];
1268 : // prefactors
1269 : Vector5d cfact = cfact_[atype];
1270 : // and position
1271 : Vector pos = getPosition(im);
1272 :
1273 : // cycle on all the neighboring voxels affected by a change in Bfactor
1274 : for(unsigned i=0; i<Model_nb_[im].size(); ++i) {
1275 : // voxel id
1276 : unsigned id = Model_nb_[im][i];
1277 : // get contribution to density in id before change
1278 : double dold = get_overlap(Map_m_[id], pos, cfact, m_s, bfactold);
1279 : // get contribution after change
1280 : double dnew = get_overlap(Map_m_[id], pos, cfact, m_s, bfactnew);
1281 : // update delta density
1282 : deltaov_l[id] += dnew-dold;
1283 : }
1284 : }
1285 : }
1286 : // add to global list, in thread order (reproducible results)
1287 180 : for(unsigned t=0; t<nt; ++t) {
1288 18445 : for(const auto & itov : deltaov_t[t]) {
1289 18325 : deltaov[itov.first] += itov.second;
1290 : }
1291 : }
1292 :
1293 : // now calculate new and old score
1294 : double old_ene = 0.0;
1295 : double new_ene = 0.0;
1296 :
1297 : // cycle on all affected voxels
1298 13661 : for(std::map<unsigned,double>::iterator itov=deltaov.begin(); itov!=deltaov.end(); ++itov) {
1299 : // id of the component
1300 13601 : unsigned id = itov->first;
1301 : // new value
1302 13601 : double ovmdnew = ovmd_[id]+itov->second;
1303 : // deviations
1304 13601 : double devold = scale_ * ovmd_[id] + offset_ - ovdd_[id];
1305 13601 : double devnew = scale_ * ovmdnew + offset_ - ovdd_[id];
1306 : // inverse of sigma_min
1307 13601 : double ismin = ismin_[id];
1308 : // scores
1309 13601 : if(devold==0.0) {
1310 0 : old_ene += -kbt_ * std::log( 0.5 * sqrt2_pi_ * ismin );
1311 : } else {
1312 13601 : old_ene += -kbt_ * std::log( 0.5 / devold * erf ( devold * inv_sqrt2_ * ismin ));
1313 : }
1314 13601 : if(devnew==0.0) {
1315 0 : new_ene += -kbt_ * std::log( 0.5 * sqrt2_pi_ * ismin );
1316 : } else {
1317 13601 : new_ene += -kbt_ * std::log( 0.5 / devnew * erf ( devnew * inv_sqrt2_ * ismin ));
1318 : }
1319 : }
1320 :
1321 : // list of neighboring residues
1322 60 : std::vector<unsigned> close = nl_res_[ir];
1323 : // add restraint to keep Bfactors of close residues close
1324 484 : for(unsigned i=0; i<close.size(); ++i) {
1325 : // residue/chain IDs of neighbor
1326 424 : std::pair<int,std::string> keyn = Model_rlist_[close[i]];
1327 : // deviations
1328 424 : double devold = bfactold - Model_b_[keyn];
1329 424 : double devnew = bfactnew - Model_b_[keyn];
1330 : // inverse of sigma_min
1331 424 : double ismin = 1.0 / bfactsig_;
1332 : // scores
1333 424 : if(devold==0.0) {
1334 204 : old_ene += -kbt_ * std::log( 0.5 * sqrt2_pi_ * ismin );
1335 : } else {
1336 220 : old_ene += -kbt_ * std::log( 0.5 / devold * erf ( devold * inv_sqrt2_ * ismin ));
1337 : }
1338 424 : if(devnew==0.0) {
1339 0 : new_ene += -kbt_ * std::log( 0.5 * sqrt2_pi_ * ismin );
1340 : } else {
1341 424 : new_ene += -kbt_ * std::log( 0.5 / devnew * erf ( devnew * inv_sqrt2_ * ismin ));
1342 : }
1343 : }
1344 :
1345 : // increment number of trials
1346 60 : MCBtrials_ += 1.0;
1347 :
1348 : // accept or reject
1349 : bool accept = false;
1350 60 : if(bfactemin_) {
1351 0 : if(new_ene < old_ene) {
1352 : accept = true;
1353 : }
1354 : } else {
1355 60 : accept = doAccept(old_ene, new_ene, kbt_);
1356 : }
1357 :
1358 : // in case of acceptance
1359 60 : if(accept) {
1360 : // update acceptance rate
1361 28 : MCBaccept_ += 1.0;
1362 : // update bfactor
1363 28 : Model_b_[key] = bfactnew;
1364 : // change all the ovmd_ affected
1365 6432 : for(std::map<unsigned,double>::iterator itov=deltaov.begin(); itov!=deltaov.end(); ++itov) {
1366 6404 : ovmd_[itov->first] += itov->second;
1367 : }
1368 : }
1369 :
1370 60 : } // end cycle on bfactors
1371 :
1372 : // update auxiliary lists (to update pref_gpu_, invs2_gpu_, and cut_ on CPU/GPU)
1373 2 : get_auxiliary_vectors();
1374 : // update neighbor list (new cut_ + update pref_nl_gpu_ and invs2_nl_gpu_ on GPU)
1375 2 : update_neighbor_list();
1376 : // recalculate fmod (to update derivatives)
1377 2 : calculate_fmod();
1378 2 : }
1379 :
1380 : // get overlap
1381 59636 : double EMMIVOX::get_overlap(const Vector &d_m, const Vector &m_m,
1382 : const Vector5d &cfact, const Vector5d &m_s, double bfact) {
1383 : // calculate vector difference
1384 : Vector md = delta(m_m, d_m);
1385 : // norm squared
1386 59636 : double md2 = md[0]*md[0]+md[1]*md[1]+md[2]*md[2];
1387 : // cycle on 5 Gaussians
1388 : double ov_tot = 0.0;
1389 357816 : for(unsigned j=0; j<5; ++j) {
1390 : // total value of b
1391 298180 : double m_b = m_s[j]+bfact/4.0;
1392 : // calculate invs2
1393 298180 : double invs2 = 1.0/(inv_pi2_*m_b);
1394 : // final calculation
1395 298180 : ov_tot += cfact[j] * pow(invs2, 1.5) * std::exp(-0.5 * md2 * invs2);
1396 : }
1397 59636 : return ov_tot;
1398 : }
1399 :
1400 4 : void EMMIVOX::update_neighbor_sphere() {
1401 : // number of atoms
1402 4 : const unsigned natoms = Model_type_.size();
1403 : // clear neighbor sphere
1404 4 : ns_.clear();
1405 : // store reference positions
1406 4 : refpos_ = getPositions();
1407 4 : if(natoms==0) {
1408 0 : return;
1409 : }
1410 :
1411 : // the atoms are binned in cells larger than the largest sphere radius (2*cut_),
1412 : // so that only the 27 cells around a voxel have to be searched
1413 4 : const double rc = 2.0 * (*std::max_element(cut_.begin(), cut_.end()));
1414 4 : Vector lo = refpos_[0], hi = refpos_[0];
1415 936 : for(unsigned im=1; im<natoms; ++im) {
1416 3728 : for(unsigned k=0; k<3; ++k) {
1417 3026 : lo[k] = std::min(lo[k], refpos_[im][k]);
1418 2796 : hi[k] = std::max(hi[k], refpos_[im][k]);
1419 : }
1420 : }
1421 : int ncell[3];
1422 : double csize[3];
1423 16 : for(unsigned k=0; k<3; ++k) {
1424 : // at most 200 cells per direction
1425 12 : const double len = hi[k]-lo[k];
1426 24 : ncell[k] = std::max(1, std::min(200, static_cast<int>(std::floor(len/rc))));
1427 12 : csize[k] = std::max(len/static_cast<double>(ncell[k]), rc);
1428 : }
1429 90804 : auto cellIndex=[&](const Vector & x, const unsigned k) {
1430 90804 : return static_cast<int>(std::floor((x[k]-lo[k])/csize[k]));
1431 4 : };
1432 4 : const unsigned ncells = ncell[0]*ncell[1]*ncell[2];
1433 : // atoms of each cell, in ascending order
1434 4 : std::vector<unsigned> cell_start(ncells+1,0), cell_atoms(natoms), atom_cell(natoms);
1435 940 : for(unsigned im=0; im<natoms; ++im) {
1436 : int c[3];
1437 3744 : for(unsigned k=0; k<3; ++k) {
1438 8113 : c[k] = std::min(ncell[k]-1, std::max(0, cellIndex(refpos_[im],k)));
1439 : }
1440 936 : atom_cell[im] = (c[0]*ncell[1]+c[1])*ncell[2]+c[2];
1441 936 : cell_start[atom_cell[im]+1]++;
1442 : }
1443 54244 : for(unsigned ic=0; ic<ncells; ++ic) {
1444 54240 : cell_start[ic+1] += cell_start[ic];
1445 : }
1446 : {
1447 4 : std::vector<unsigned> fill(cell_start.begin(), cell_start.end()-1);
1448 940 : for(unsigned im=0; im<natoms; ++im) {
1449 936 : cell_atoms[fill[atom_cell[im]]++] = im;
1450 : }
1451 4 : }
1452 :
1453 : // cycle on voxels - in parallel; the pairs are sorted by voxel and then by atom
1454 4 : const unsigned nt = OpenMP::getNumThreads();
1455 4 : std::vector< std::vector< std::pair<unsigned,unsigned> > > ns_t(nt);
1456 4 : #pragma omp parallel num_threads(nt)
1457 : {
1458 : // private variables
1459 : std::vector< std::pair<unsigned,unsigned> > & ns_l = ns_t[OpenMP::getThreadNum()];
1460 : std::vector<unsigned> close;
1461 : #pragma omp for schedule(static)
1462 : for(unsigned id=0; id<ovdd_.size(); ++id) {
1463 : // grid point
1464 : const Vector d_m = Map_m_[id];
1465 : int c[3];
1466 : bool outside = false;
1467 : for(unsigned k=0; k<3; ++k) {
1468 : c[k] = cellIndex(d_m,k);
1469 : // more than one cell away from the atoms
1470 : if(c[k]<-1 || c[k]>ncell[k]) {
1471 : outside = true;
1472 : }
1473 : }
1474 : if(outside) {
1475 : continue;
1476 : }
1477 : close.clear();
1478 : for(int i0=std::max(0,c[0]-1); i0<=std::min(ncell[0]-1,c[0]+1); ++i0) {
1479 : for(int i1=std::max(0,c[1]-1); i1<=std::min(ncell[1]-1,c[1]+1); ++i1) {
1480 : for(int i2=std::max(0,c[2]-1); i2<=std::min(ncell[2]-1,c[2]+1); ++i2) {
1481 : const unsigned ic = (i0*ncell[1]+i1)*ncell[2]+i2;
1482 : for(unsigned j=cell_start[ic]; j<cell_start[ic+1]; ++j) {
1483 : const unsigned im = cell_atoms[j];
1484 : // calculate distance
1485 : const double dist = delta(getPosition(im), d_m).modulo();
1486 : // add to local list
1487 : if(dist<=2.0*cut_[im]) {
1488 : close.push_back(im);
1489 : }
1490 : }
1491 : }
1492 : }
1493 : }
1494 : std::sort(close.begin(), close.end());
1495 : for(const auto im : close) {
1496 : ns_l.push_back(std::make_pair(id,im));
1497 : }
1498 : }
1499 : }
1500 : // add to global list, in thread order
1501 12 : for(unsigned t=0; t<nt; ++t) {
1502 8 : ns_.insert(ns_.end(), ns_t[t].begin(), ns_t[t].end());
1503 : }
1504 4 : }
1505 :
1506 7 : bool EMMIVOX::do_neighbor_sphere() {
1507 7 : std::vector<double> dist(getPositions().size());
1508 : bool update = false;
1509 :
1510 : // calculate displacement
1511 7 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1512 : for(unsigned im=0; im<dist.size(); ++im) {
1513 : dist[im] = delta(getPosition(im),refpos_[im]).modulo()/cut_[im];
1514 : }
1515 :
1516 : // check if update or not
1517 7 : double maxdist = *max_element(dist.begin(), dist.end());
1518 7 : if(maxdist>=1.0) {
1519 : update=true;
1520 : }
1521 :
1522 : // return if update or not
1523 7 : return update;
1524 : }
1525 :
1526 11 : void EMMIVOX::update_neighbor_list() {
1527 : // number of atoms
1528 : unsigned natoms = Model_type_.size();
1529 : // clear neighbor list
1530 11 : nl_.clear();
1531 :
1532 : // cycle on neighbour sphere - in parallel, the order of the pairs is preserved
1533 11 : const unsigned nt = OpenMP::getNumThreads();
1534 11 : std::vector< std::vector< std::pair<unsigned,unsigned> > > nl_t(nt);
1535 11 : #pragma omp parallel num_threads(nt)
1536 : {
1537 : // private variables
1538 : std::vector< std::pair<unsigned,unsigned> > & nl_l = nl_t[OpenMP::getThreadNum()];
1539 : #pragma omp for schedule(static)
1540 : for(unsigned long long i=0; i<ns_.size(); ++i) {
1541 : // calculate distance
1542 : double dist = delta(Map_m_[ns_[i].first], getPosition(ns_[i].second)).modulo();
1543 : // add to local neighbour list
1544 : if(dist<=cut_[ns_[i].second]) {
1545 : nl_l.push_back(ns_[i]);
1546 : }
1547 : }
1548 : }
1549 : // add to global list, in thread order
1550 33 : for(unsigned t=0; t<nt; ++t) {
1551 22 : nl_.insert(nl_.end(), nl_t[t].begin(), nl_t[t].end());
1552 : }
1553 :
1554 : // new dimension of neighbor list
1555 : unsigned long long nl_size = nl_.size();
1556 : // now resize derivatives
1557 11 : ovmd_der_.resize(nl_size);
1558 :
1559 : // in case of B-factors sampling - at the right step
1560 11 : if(dbfact_>0 && getStep()%MCBstride_==0) {
1561 : // list of voxels per atom, in the order of the neighbor list
1562 5 : Model_nb_.clear();
1563 5 : Model_nb_.resize(natoms);
1564 75426 : for(unsigned long long i=0; i<nl_size; ++i) {
1565 75421 : Model_nb_[nl_[i].second].push_back(nl_[i].first);
1566 : }
1567 : }
1568 :
1569 : // transfer data to gpu
1570 11 : update_gpu();
1571 11 : }
1572 :
1573 11 : void EMMIVOX::update_gpu() {
1574 : // dimension of neighbor list
1575 11 : unsigned nl_size = nl_.size();
1576 : // create useful vectors
1577 11 : std::vector<int> nl_id(nl_size), nl_im(nl_size);
1578 11 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1579 : for(unsigned long long i=0; i<nl_size; ++i) {
1580 : nl_id[i] = static_cast<int>(nl_[i].first);
1581 : nl_im[i] = static_cast<int>(nl_[i].second);
1582 : }
1583 : // create tensors on device
1584 44 : nl_id_gpu_ = torch::from_blob(nl_id.data(), {nl_size}, torch::kInt32).clone().to(device_t_);
1585 33 : nl_im_gpu_ = torch::from_blob(nl_im.data(), {nl_size}, torch::kInt32).clone().to(device_t_);
1586 : // now we need to create pref_nl_gpu_ [5,nl_size]
1587 11 : pref_nl_gpu_ = torch::index_select(pref_gpu_,1,nl_im_gpu_);
1588 : // and invs2_nl_gpu_ [5,nl_size]
1589 11 : invs2_nl_gpu_ = torch::index_select(invs2_gpu_,1,nl_im_gpu_);
1590 : // and Map_m_nl_gpu_ [3,nl_size]
1591 11 : Map_m_nl_gpu_ = torch::index_select(Map_m_gpu_,1,nl_id_gpu_);
1592 11 : }
1593 :
1594 12 : void EMMIVOX::prepare() {
1595 12 : if(getExchangeStep()) {
1596 0 : first_time_=true;
1597 : }
1598 12 : }
1599 :
1600 : // calculate forward model on gpu
1601 14 : void EMMIVOX::calculate_fmod() {
1602 : // number of atoms
1603 14 : int natoms = Model_type_.size();
1604 : // number of data points
1605 14 : int nd = ovdd_.size();
1606 :
1607 : // fill positions in parallel
1608 14 : std::vector<double> posg(3*natoms);
1609 14 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1610 : for (int i=0; i<natoms; ++i) {
1611 : // fill vectors
1612 : posg[i] = getPosition(i)[0];
1613 : posg[i+natoms] = getPosition(i)[1];
1614 : posg[i+2*natoms] = getPosition(i)[2];
1615 : }
1616 : // transfer positions to pos_gpu [3,natoms]
1617 28 : torch::Tensor pos_gpu = torch::from_blob(posg.data(), {3,natoms}, torch::kFloat64).to(torch::kFloat32).to(device_t_);
1618 : // create pos_nl_gpu_ [3,nl_size]
1619 14 : torch::Tensor pos_nl_gpu = torch::index_select(pos_gpu,1,nl_im_gpu_);
1620 : // calculate vector difference [3,nl_size]
1621 14 : torch::Tensor md = Map_m_nl_gpu_ - pos_nl_gpu;
1622 : // calculate norm squared by column [1,nl_size]
1623 28 : torch::Tensor md2 = torch::sum(md*md,0);
1624 : // calculate density [5,nl_size]
1625 42 : torch::Tensor ov = pref_nl_gpu_ * torch::exp(-0.5 * md2 * invs2_nl_gpu_);
1626 : // and derivatives [5,nl_size]
1627 14 : ovmd_der_gpu_ = invs2_nl_gpu_ * ov;
1628 : // sum density over 5 columns [1,nl_size]
1629 28 : ov = torch::sum(ov,0);
1630 : // sum contributions from the same atom
1631 14 : auto options = torch::TensorOptions().device(device_t_).dtype(torch::kFloat32);
1632 14 : ovmd_gpu_ = torch::zeros({nd}, options);
1633 14 : ovmd_gpu_.index_add_(0, nl_id_gpu_, ov);
1634 : // sum derivatives over 5 rows [1,nl_size] and multiply by md [3,nl_size]
1635 14 : ovmd_der_gpu_ = md * torch::sum(ovmd_der_gpu_,0);
1636 :
1637 : // in case of metainference: average them across replicas
1638 14 : if(!no_aver_ && nrep_>1) {
1639 : // communicate ovmd_gpu_ to CPU [1, nd]
1640 0 : torch::Tensor ovmd_cpu = ovmd_gpu_.detach().to(torch::kCPU).to(torch::kFloat64);
1641 : // and put them in ovmd_
1642 0 : ovmd_ = std::vector<double>(ovmd_cpu.data_ptr<double>(), ovmd_cpu.data_ptr<double>() + ovmd_cpu.numel());
1643 : // sum across replicas
1644 0 : multi_sim_comm.Sum(&ovmd_[0], nd);
1645 : // and divide by number of replicas
1646 0 : double escale = 1.0 / static_cast<double>(nrep_);
1647 0 : for(int i=0; i<nd; ++i) {
1648 0 : ovmd_[i] *= escale;
1649 : }
1650 : // put back on device
1651 0 : ovmd_gpu_ = torch::from_blob(ovmd_.data(), {nd}, torch::kFloat64).to(torch::kFloat32).to(device_t_);
1652 : }
1653 :
1654 : // communicate back model density
1655 : // this is needed only in certain situations
1656 14 : long int step = getStep();
1657 : bool do_comm = false;
1658 14 : if(mapstride_>0 && step%mapstride_==0) {
1659 : do_comm = true;
1660 : }
1661 14 : if(dbfact_>0 && step%MCBstride_==0) {
1662 : do_comm = true;
1663 : }
1664 14 : if(do_corr_) {
1665 : do_comm = true;
1666 : }
1667 : // in case of metainference: already communicated
1668 14 : if(!no_aver_ && nrep_>1) {
1669 : do_comm = false;
1670 : }
1671 14 : if(do_comm) {
1672 : // communicate ovmd_gpu_ to CPU [1, nd]
1673 33 : torch::Tensor ovmd_cpu = ovmd_gpu_.detach().to(torch::kCPU).to(torch::kFloat64);
1674 : // and put them in ovmd_
1675 33 : ovmd_ = std::vector<double>(ovmd_cpu.data_ptr<double>(), ovmd_cpu.data_ptr<double>() + ovmd_cpu.numel());
1676 : }
1677 14 : }
1678 :
1679 : // calculate score
1680 12 : void EMMIVOX::calculate_score() {
1681 : // number of atoms
1682 12 : int natoms = Model_type_.size();
1683 :
1684 : // calculate deviation model/data [1, nd]
1685 36 : torch::Tensor dev = scale_ * ovmd_gpu_ + offset_ - ovdd_gpu_;
1686 : // error function [1, nd]
1687 24 : torch::Tensor errf = torch::erf( dev * inv_sqrt2_ * ismin_gpu_ );
1688 : // take care of dev = zero
1689 12 : torch::Tensor zeros_d = torch::ne(dev, 0.0);
1690 : // redefine dev
1691 60 : dev = dev * zeros_d + eps_ * torch::logical_not(zeros_d);
1692 : // take care of errf = zero
1693 12 : torch::Tensor zeros_e = torch::ne(errf, 0.0);
1694 : // redefine errf
1695 60 : errf = errf * zeros_e + eps_ * torch::logical_not(zeros_e);
1696 : // logical AND: both dev and errf different from zero
1697 : torch::Tensor zeros = torch::logical_and(zeros_d, zeros_e);
1698 : // energy - with limit dev going to zero
1699 60 : torch::Tensor ene = 0.5 * ( errf / dev * zeros + torch::logical_not(zeros) * sqrt2_pi_ * ismin_gpu_);
1700 : // logarithm and sum
1701 : // the sum over the voxels is done in double precision
1702 36 : ene = -kbt_ * torch::sum(torch::log(ene).to(torch::kFloat64));
1703 : // and derivatives [1, nd]
1704 72 : torch::Tensor d_der = -kbt_ * zeros * ( sqrt2_pi_ * torch::exp( -0.5 * dev * dev * ismin_gpu_ * ismin_gpu_ ) * ismin_gpu_ / errf - 1.0 / dev );
1705 : // tensor for derivatives wrt atoms [1, nl_size]
1706 12 : torch::Tensor der_gpu = torch::index_select(d_der,0,nl_id_gpu_);
1707 : // multiply by ovmd_der_gpu_ and scale [3, nl_size]
1708 24 : der_gpu = ovmd_der_gpu_ * scale_ * der_gpu;
1709 : // sum contributions for each atom
1710 12 : auto options = torch::TensorOptions().device(device_t_).dtype(torch::kFloat32);
1711 12 : torch::Tensor atoms_der_gpu = torch::zeros({3,natoms}, options);
1712 24 : atoms_der_gpu.index_add_(1, nl_im_gpu_, der_gpu);
1713 :
1714 : // FINAL STUFF
1715 : //
1716 : // 1) communicate total energy to CPU
1717 36 : torch::Tensor ene_cpu = ene.detach().to(torch::kCPU).to(torch::kFloat64);
1718 12 : ene_ = *ene_cpu.data_ptr<double>();
1719 : // with marginal, simply multiply by number of replicas!
1720 12 : if(!no_aver_ && nrep_>1) {
1721 0 : ene_ *= static_cast<double>(nrep_);
1722 : }
1723 : //
1724 : // 2) communicate derivatives to CPU
1725 36 : torch::Tensor atom_der_cpu = atoms_der_gpu.detach().to(torch::kCPU).to(torch::kFloat64);
1726 : // convert to std::vector<double>
1727 24 : std::vector<double> atom_der = std::vector<double>(atom_der_cpu.data_ptr<double>(), atom_der_cpu.data_ptr<double>() + atom_der_cpu.numel());
1728 : // and put in atom_der_
1729 12 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1730 : for(int i=0; i<natoms; ++i) {
1731 : atom_der_[i] = Vector(atom_der[i],atom_der[i+natoms],atom_der[i+2*natoms]);
1732 : }
1733 : //
1734 : // 3) calculate virial on CPU
1735 12 : Tensor virial;
1736 : // declare omp reduction for Tensors
1737 : #pragma omp declare reduction( sumTensor : Tensor : omp_out += omp_in )
1738 :
1739 12 : #pragma omp parallel for num_threads(OpenMP::getNumThreads()) reduction (sumTensor : virial)
1740 : for(int i=0; i<natoms; ++i) {
1741 : virial += Tensor(getPosition(i), -atom_der_[i]);
1742 : }
1743 : // store virial
1744 12 : virial_ = virial;
1745 12 : }
1746 :
1747 12 : void EMMIVOX::calculate() {
1748 : // get time step
1749 12 : long int step = getStep();
1750 :
1751 : // set temperature value
1752 12 : getPntrToComponent("kbt")->set(kbt_);
1753 :
1754 : // neighbor list update
1755 12 : if(first_time_ || getExchangeStep() || step%nl_stride_==0) {
1756 : // check if time to update neighbor sphere
1757 : bool update = false;
1758 9 : if(first_time_ || getExchangeStep()) {
1759 : update = true;
1760 : } else {
1761 7 : update = do_neighbor_sphere();
1762 : }
1763 : // update neighbor sphere
1764 7 : if(update) {
1765 4 : update_neighbor_sphere();
1766 : }
1767 : // update neighbor list
1768 9 : update_neighbor_list();
1769 : // set flag
1770 9 : first_time_=false;
1771 : }
1772 :
1773 : // calculate forward model
1774 12 : calculate_fmod();
1775 :
1776 : // Monte Carlo on bfactors
1777 12 : if(dbfact_>0) {
1778 : double acc = 0.0;
1779 : // do Monte Carlo
1780 6 : if(step%MCBstride_==0 && !getExchangeStep() && step>0) {
1781 2 : doMonteCarloBfact();
1782 : }
1783 : // calculate acceptance ratio
1784 6 : if(MCBtrials_>0) {
1785 4 : acc = MCBaccept_ / MCBtrials_;
1786 : }
1787 : // set value
1788 12 : getPntrToComponent("accB")->set(acc);
1789 : }
1790 :
1791 : // calculate score
1792 12 : calculate_score();
1793 :
1794 : // set score, virial, and derivatives
1795 12 : Value* score = getPntrToComponent("scoreb");
1796 12 : score->set(ene_);
1797 12 : setBoxDerivatives(score, virial_);
1798 12 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1799 : for(unsigned i=0; i<atom_der_.size(); ++i) {
1800 : setAtomsDerivatives(score, i, atom_der_[i]);
1801 : }
1802 : // set scale and offset value
1803 12 : getPntrToComponent("scale")->set(scale_);
1804 12 : getPntrToComponent("offset")->set(offset_);
1805 : // calculate correlation coefficient
1806 12 : if(do_corr_) {
1807 6 : calculate_corr();
1808 : }
1809 : // PRINT other quantities to files
1810 : // - status file
1811 12 : if(step%statusstride_==0) {
1812 6 : print_status(step);
1813 : }
1814 : // - density file
1815 12 : if(mapstride_>0 && step%mapstride_==0) {
1816 2 : write_model_density(step);
1817 : }
1818 12 : }
1819 :
1820 6 : void EMMIVOX::calculate_corr() {
1821 : // number of data points
1822 6 : double nd = static_cast<double>(ovdd_.size());
1823 : // average ovmd_ and ovdd_
1824 6 : double ave_md = std::accumulate(ovmd_.begin(), ovmd_.end(), 0.) / nd;
1825 6 : double ave_dd = std::accumulate(ovdd_.begin(), ovdd_.end(), 0.) / nd;
1826 : // calculate correlation
1827 : double num = 0.;
1828 : double den1 = 0.;
1829 : double den2 = 0.;
1830 6 : #pragma omp parallel for num_threads(OpenMP::getNumThreads()) reduction( + : num, den1, den2)
1831 : for(unsigned i=0; i<ovdd_.size(); ++i) {
1832 : double md = ovmd_[i]-ave_md;
1833 : double dd = ovdd_[i]-ave_dd;
1834 : num += md*dd;
1835 : den1 += md*md;
1836 : den2 += dd*dd;
1837 : }
1838 : // correlation coefficient
1839 6 : double cc = num / sqrt(den1*den2);
1840 : // set plumed
1841 6 : getPntrToComponent("corr")->set(cc);
1842 6 : }
1843 :
1844 :
1845 : }
1846 : }
1847 :
1848 : #endif
|
