Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionWithValue.h"
23 : #include "core/ActionAtomistic.h"
24 : #include "core/ActionWithArguments.h"
25 : #include "core/ActionRegister.h"
26 : #include "core/PlumedMain.h"
27 : #include "tools/Matrix.h"
28 : #include "tools/Communicator.h"
29 : #include "core/GenericMolInfo.h"
30 : #include "core/ActionSet.h"
31 : #include "tools/File.h"
32 : #include "tools/Random.h"
33 : #include "tools/OpenMP.h"
34 :
35 : #include <string>
36 : #include <map>
37 : #include <numeric>
38 : #include <ctime>
39 : #include <algorithm>
40 : #include <utility>
41 :
42 : namespace PLMD {
43 : namespace isdb {
44 :
45 : //+PLUMEDOC ISDB_COLVAR EMMI
46 : /*
47 : Calculate the fit of a structure or ensemble of structures with a cryo-EM density map.
48 :
49 : This action implements the multi-scale Bayesian approach to cryo-EM data fitting introduced in the first paper cited below.
50 : This method allows efficient and accurate structural modeling of cryo-electron microscopy density maps at multiple scales, from coarse-grained to atomistic resolution, by addressing the presence of random and systematic errors in the data, sample heterogeneity, data correlation, and noise correlation.
51 :
52 : The experimental density map is fit by a Gaussian Mixture Model (GMM), which is provided as an external file specified by the keyword
53 : GMM_FILE. We are currently working on a web server to perform
54 : this operation. In the meantime, the user can request a stand-alone version of the GMM code at massimiliano.bonomi_AT_gmail.com.
55 :
56 : When run in single-replica mode, this action allows atomistic, flexible refinement of an individual structure into a density map.
57 : Combined with a multi-replica framework (such as the -multi option in GROMACS), the user can model an ensemble of structures using
58 : the Metainference approach that is discussed in the second paper cited below.
59 :
60 : !!! warning ""
61 :
62 : To use [EMMI](EMMI.md), the user should always add a [MOLINFO](MOLINFO.md) line and specify a pdb file of the system.
63 :
64 : !!! note ""
65 :
66 : To enhance sampling in single-structure refinement, one can use a Replica Exchange Method, such as Parallel Tempering.
67 : In this case, the user should add the NO_AVER flag to the input line. To use a replica-based enhanced sampling scheme such as
68 : Parallel-Bias Metadynamics ([PBMETAD](PBMETAD.md)), one should use the REWEIGHT flag and pass the Metadynamics bias using the ARG keyword.
69 :
70 : !!! note ""
71 :
72 : [EMMI](EMMI.md) can be used in combination with periodic and non-periodic systems. In the latter case, one should
73 : add the NOPBC flag to the input line
74 :
75 : ## Examples
76 :
77 : In this example, we perform a single-structure refinement based on an experimental cryo-EM map. The map is fit with a GMM, whose
78 : parameters are listed in the file GMM_fit.dat. This file contains one line per GMM component in the following format:
79 :
80 : ````
81 : #! FIELDS Id Weight Mean_0 Mean_1 Mean_2 Cov_00 Cov_01 Cov_02 Cov_11 Cov_12 Cov_22 Beta
82 : 0 2.9993805e+01 6.54628 10.37820 -0.92988 2.078920e-02 1.216254e-03 5.990827e-04 2.556246e-02 8.411835e-03 2.486254e-02 1
83 : 1 2.3468312e+01 6.56095 10.34790 -0.87808 1.879859e-02 6.636049e-03 3.682865e-04 3.194490e-02 1.750524e-03 3.017100e-02 1
84 : ...
85 : ````
86 :
87 : To accelerate the computation of the Bayesian score, one can:
88 :
89 : - use neighbor lists, specified by the keywords NL_CUTOFF and NL_STRIDE;
90 : - calculate the restraint every other step (or more).
91 :
92 : All the heavy atoms of the system are used to calculate the density map. This list can conveniently be provided
93 : using a GROMACS index file.
94 :
95 : The input file looks as follows:
96 :
97 : ```plumed
98 : # include pdb info
99 : MOLINFO STRUCTURE=prot.pdb
100 :
101 : # all heavy atoms
102 : protein-h: GROUP NDX_FILE=index.ndx NDX_GROUP=Protein-H
103 :
104 : # create EMMI score
105 : gmm: EMMI NOPBC SIGMA_MIN=0.01 TEMP=300.0 NL_STRIDE=100 NL_CUTOFF=0.01 GMM_FILE=GMM_fit.dat ATOMS=protein-h
106 :
107 : # translate into bias - apply every 2 steps
108 : emr: BIASVALUE ARG=gmm.scoreb STRIDE=2
109 :
110 : PRINT ARG=emr.* FILE=COLVAR STRIDE=500 FMT=%20.10f
111 : ```
112 :
113 : */
114 : //+ENDPLUMEDOC
115 :
116 : class EMMI :
117 : public ActionAtomistic,
118 : public ActionWithArguments,
119 : public ActionWithValue {
120 : private:
121 :
122 : // temperature in kbt
123 : double kbt_;
124 : // model GMM - atom types
125 : std::vector<unsigned> GMM_m_type_;
126 : // model GMM - list of atom sigmas - one per atom type
127 : std::vector<double> GMM_m_s_;
128 : // model GMM - list of atom weights - one per atom type
129 : std::vector<double> GMM_m_w_;
130 : // data GMM - means, weights, and covariances + beta option
131 : std::vector<Vector> GMM_d_m_;
132 : std::vector<double> GMM_d_w_;
133 : std::vector< VectorGeneric<6> > GMM_d_cov_;
134 : std::vector<int> GMM_d_beta_;
135 : std::vector < std::vector<int> > GMM_d_grps_;
136 : // overlaps
137 : std::vector<double> ovmd_;
138 : std::vector<double> ovdd_;
139 : std::vector<double> ovmd_ave_;
140 : // and derivatives
141 : std::vector<Vector> ovmd_der_;
142 : std::vector<Vector> atom_der_;
143 : std::vector<double> GMMid_der_;
144 : // constants
145 : double cfact_;
146 : double inv_sqrt2_, sqrt2_pi_;
147 : // metainference
148 : unsigned nrep_;
149 : unsigned replica_;
150 : std::vector<double> sigma_;
151 : std::vector<double> sigma_min_;
152 : std::vector<double> sigma_max_;
153 : std::vector<double> dsigma_;
154 : // list of prefactors for overlap between two Gaussians
155 : // pre_fact = 1.0 / (2pi)**1.5 / std::sqrt(det_md) * Wm * Wd
156 : std::vector<double> pre_fact_;
157 : // inverse of the sum of model and data covariances matrices
158 : std::vector< VectorGeneric<6> > inv_cov_md_;
159 : // neighbor list
160 : double nl_cutoff_;
161 : unsigned nl_stride_;
162 : bool first_time_;
163 : bool no_aver_;
164 : // neighbor list pairs [nl_begin_,nl_end_) are computed by this rank
165 : unsigned nl_begin_;
166 : unsigned nl_end_;
167 : // components
168 : Value* valueScore_;
169 : Value* valueAcc_;
170 : std::vector<Value*> valueSigma_;
171 : std::vector<unsigned> nl_;
172 : // parallel stuff
173 : unsigned size_;
174 : unsigned rank_;
175 : // pbc
176 : bool pbc_;
177 : // Monte Carlo stuff
178 : int MCstride_;
179 : double MCaccept_;
180 : double MCtrials_;
181 : Random random_;
182 : // status stuff
183 : unsigned int statusstride_;
184 : std::string statusfilename_;
185 : OFile statusfile_;
186 : bool first_status_;
187 : // regression
188 : unsigned nregres_;
189 : double scale_;
190 : double scale_min_;
191 : double scale_max_;
192 : double dscale_;
193 : // tabulated exponential
194 : double dpcutoff_;
195 : double dexp_;
196 : unsigned nexp_;
197 : std::vector<double> tab_exp_;
198 : // simulated annealing
199 : unsigned nanneal_;
200 : double kanneal_;
201 : double anneal_;
202 : // prior exponent
203 : double prior_;
204 : // noise type
205 : unsigned noise_;
206 : // total score and virial;
207 : double ene_;
208 : Tensor virial_;
209 : // model overlap file
210 : unsigned int ovstride_;
211 : std::string ovfilename_;
212 :
213 : // Reweighting additions
214 : bool do_reweight_;
215 : bool first_time_w_;
216 : std::vector<double> forces;
217 : std::vector<double> forcesToApply;
218 :
219 : // average weights
220 : double decay_w_;
221 : std::vector<double> average_weights_;
222 :
223 : // write file with model overlap
224 : void write_model_overlap(long long int step);
225 : // get median of std::vector
226 : double get_median(std::vector<double> &v);
227 : // annealing
228 : double get_annealing(long long int step);
229 : // do regression
230 : double scaleEnergy(double s);
231 : double doRegression();
232 : // read and write status
233 : void read_status();
234 : void print_status(long long int step);
235 : // accept or reject
236 : bool doAccept(double oldE, double newE, double kbt);
237 : // do MonteCarlo
238 : void doMonteCarlo();
239 : void doMonteCarloMove();
240 : // read error file
241 : std::vector<double> read_exp_errors(const std::string & errfile);
242 : // read experimental overlaps
243 : std::vector<double> read_exp_overlaps(const std::string & ovfile);
244 : // calculate model GMM weights and covariances
245 : std::vector<double> get_GMM_m(std::vector<AtomNumber> &atoms);
246 : // read data GMM file
247 : void get_GMM_d(const std::string & gmm_file);
248 : // check GMM data
249 : void check_GMM_d(const VectorGeneric<6> &cov, double w);
250 : // auxiliary method
251 : void calculate_useful_stuff(double reso);
252 : // get pref_fact and inv_cov_md
253 : double get_prefactor_inverse (const VectorGeneric<6> &GMM_cov_0, const VectorGeneric<6> &GMM_cov_1,
254 : double GMM_w_0, double GMM_w_1,
255 : VectorGeneric<6> &sum, VectorGeneric<6> &inv_sum);
256 : // calculate self overlaps between data GMM components - ovdd_
257 : double get_self_overlap(unsigned id);
258 : // calculate overlap between two Gaussians
259 : double get_overlap(const Vector &m_m, const Vector &d_m, double pre_fact,
260 : const VectorGeneric<6> &inv_cov_md, Vector &ov_der);
261 : // calculate exponent of overlap for neighbor list update
262 : double get_exp_overlap(const Vector &m_m, const Vector &d_m,
263 : const VectorGeneric<6> &inv_cov_md);
264 : // update the neighbor list
265 : void update_neighbor_list();
266 : // calculate overlap
267 : void calculate_overlap();
268 : // Gaussian noise
269 : void calculate_Gauss();
270 : // Outliers noise
271 : void calculate_Outliers();
272 : // Marginal noise
273 : void calculate_Marginal();
274 :
275 : // See MetainferenceBase
276 : void get_weights(double &weight, double &norm, double &neff);
277 :
278 : public:
279 : static void registerKeywords( Keywords& keys );
280 : explicit EMMI(const ActionOptions&);
281 : // needed for reweighting
282 : void setDerivatives();
283 : void turnOnDerivatives() override;
284 : unsigned getNumberOfDerivatives() override;
285 : void lockRequests() override;
286 : void unlockRequests() override;
287 : void calculateNumericalDerivatives( ActionWithValue* a ) override;
288 : void apply() override;
289 : void setArgDerivatives(Value *v, const double &d);
290 : void setAtomsDerivatives(Value*v, const unsigned i, const Vector&d);
291 : void setBoxDerivatives(Value*v, const Tensor&d);
292 : // active methods:
293 : void prepare() override;
294 : void calculate() override;
295 : };
296 :
297 : inline
298 28 : void EMMI::setDerivatives() {
299 : // Get appropriate number of derivatives
300 : // Derivatives are first for arguments and then for atoms
301 : unsigned nder;
302 28 : if( getNumberOfAtoms()>0 ) {
303 28 : nder = 3*getNumberOfAtoms() + 9 + getNumberOfArguments();
304 : } else {
305 0 : nder = getNumberOfArguments();
306 : }
307 :
308 : // Resize all derivative arrays
309 28 : forces.resize( nder );
310 28 : forcesToApply.resize( nder );
311 164 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
312 136 : getPntrToComponent(i)->resizeDerivatives(nder);
313 : }
314 28 : }
315 :
316 : inline
317 24 : void EMMI::turnOnDerivatives() {
318 24 : ActionWithValue::turnOnDerivatives();
319 24 : }
320 :
321 : inline
322 132 : unsigned EMMI::getNumberOfDerivatives() {
323 132 : if( getNumberOfAtoms()>0 ) {
324 132 : return 3*getNumberOfAtoms() + 9 + getNumberOfArguments();
325 : }
326 0 : return getNumberOfArguments();
327 : }
328 :
329 : inline
330 108 : void EMMI::lockRequests() {
331 : ActionAtomistic::lockRequests();
332 : ActionWithArguments::lockRequests();
333 108 : }
334 :
335 : inline
336 108 : void EMMI::unlockRequests() {
337 : ActionAtomistic::unlockRequests();
338 : ActionWithArguments::unlockRequests();
339 108 : }
340 :
341 : inline
342 0 : void EMMI::calculateNumericalDerivatives( ActionWithValue* a=NULL ) {
343 0 : if( getNumberOfArguments()>0 ) {
344 0 : ActionWithArguments::calculateNumericalDerivatives( a );
345 : }
346 0 : if( getNumberOfAtoms()>0 ) {
347 0 : Matrix<double> save_derivatives( getNumberOfComponents(), getNumberOfArguments() );
348 0 : for(unsigned j=0; j<getNumberOfComponents(); ++j) {
349 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
350 0 : if(getPntrToComponent(j)->hasDerivatives()) {
351 0 : save_derivatives(j,i)=getPntrToComponent(j)->getDerivative(i);
352 : }
353 : }
354 : }
355 0 : calculateAtomicNumericalDerivatives( a, getNumberOfArguments() );
356 0 : for(unsigned j=0; j<getNumberOfComponents(); ++j) {
357 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
358 0 : if(getPntrToComponent(j)->hasDerivatives()) {
359 0 : getPntrToComponent(j)->addDerivative( i, save_derivatives(j,i) );
360 : }
361 : }
362 : }
363 : }
364 0 : }
365 :
366 : inline
367 108 : void EMMI::apply() {
368 : bool wasforced=false;
369 108 : forcesToApply.assign(forcesToApply.size(),0.0);
370 710 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
371 602 : if( getPntrToComponent(i)->applyForce( forces ) ) {
372 : wasforced=true;
373 480 : for(unsigned ii=0; ii<forces.size(); ++ii) {
374 468 : forcesToApply[ii]+=forces[ii];
375 : }
376 : }
377 : }
378 108 : if( wasforced ) {
379 12 : unsigned ind=0;
380 12 : addForcesOnArguments( 0, forcesToApply, ind );
381 12 : if( getNumberOfAtoms()>0 ) {
382 12 : setForcesOnAtoms( forcesToApply, ind );
383 : }
384 : }
385 108 : }
386 :
387 : inline
388 : void EMMI::setArgDerivatives(Value *v, const double &d) {
389 : v->addDerivative(0,d);
390 : }
391 :
392 : inline
393 29496 : void EMMI::setAtomsDerivatives(Value*v, const unsigned i, const Vector&d) {
394 29496 : const unsigned noa=getNumberOfArguments();
395 29496 : v->addDerivative(noa+3*i+0,d[0]);
396 29496 : v->addDerivative(noa+3*i+1,d[1]);
397 29496 : v->addDerivative(noa+3*i+2,d[2]);
398 29496 : }
399 :
400 : inline
401 108 : void EMMI::setBoxDerivatives(Value* v,const Tensor&d) {
402 108 : const unsigned noa=getNumberOfArguments();
403 : const unsigned nat=getNumberOfAtoms();
404 108 : v->addDerivative(noa+3*nat+0,d(0,0));
405 108 : v->addDerivative(noa+3*nat+1,d(0,1));
406 108 : v->addDerivative(noa+3*nat+2,d(0,2));
407 108 : v->addDerivative(noa+3*nat+3,d(1,0));
408 108 : v->addDerivative(noa+3*nat+4,d(1,1));
409 108 : v->addDerivative(noa+3*nat+5,d(1,2));
410 108 : v->addDerivative(noa+3*nat+6,d(2,0));
411 108 : v->addDerivative(noa+3*nat+7,d(2,1));
412 108 : v->addDerivative(noa+3*nat+8,d(2,2));
413 108 : }
414 :
415 : PLUMED_REGISTER_ACTION(EMMI,"EMMI")
416 :
417 32 : void EMMI::registerKeywords( Keywords& keys ) {
418 32 : Action::registerKeywords( keys );
419 32 : ActionAtomistic::registerKeywords( keys );
420 32 : ActionWithValue::registerKeywords( keys );
421 32 : ActionWithArguments::registerKeywords( keys );
422 64 : keys.addInputKeyword("optional","ARG","scalar","the labels of the values from which the function is calculated");
423 32 : keys.add("atoms","ATOMS","atoms for which we calculate the density map, typically all heavy atoms");
424 32 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
425 32 : keys.add("compulsory","GMM_FILE","file with the parameters of the GMM components");
426 32 : keys.add("compulsory","NL_CUTOFF","The cutoff in overlap for the neighbor list");
427 32 : keys.add("compulsory","NL_STRIDE","The frequency with which we are updating the neighbor list");
428 32 : keys.add("compulsory","SIGMA_MIN","minimum uncertainty");
429 32 : keys.add("compulsory","RESOLUTION", "Cryo-EM map resolution");
430 32 : keys.add("compulsory","NOISETYPE","functional form of the noise (GAUSS, OUTLIERS, MARGINAL)");
431 32 : keys.add("optional","SIGMA0","initial value of the uncertainty");
432 32 : keys.add("optional","DSIGMA","MC step for uncertainties");
433 32 : keys.add("optional","MC_STRIDE", "Monte Carlo stride");
434 32 : keys.add("optional","ERR_FILE","file with experimental or GMM fit errors");
435 32 : keys.add("optional","OV_FILE","file with experimental overlaps");
436 32 : keys.add("optional","NORM_DENSITY","integral of the experimental density");
437 32 : keys.add("optional","STATUS_FILE","write a file with all the data useful for restart");
438 32 : keys.add("optional","WRITE_STRIDE","write the status to a file every N steps, this can be used for restart");
439 32 : keys.add("optional","REGRESSION","regression stride");
440 32 : keys.add("optional","REG_SCALE_MIN","regression minimum scale");
441 32 : keys.add("optional","REG_SCALE_MAX","regression maximum scale");
442 32 : keys.add("optional","REG_DSCALE","regression maximum scale MC move");
443 32 : keys.add("optional","SCALE","scale factor");
444 32 : keys.add("optional","ANNEAL", "Length of annealing cycle");
445 32 : keys.add("optional","ANNEAL_FACT", "Annealing temperature factor");
446 32 : keys.add("optional","TEMP","temperature");
447 32 : keys.add("optional","PRIOR", "exponent of uncertainty prior");
448 32 : keys.add("optional","WRITE_OV_STRIDE","write model overlaps every N steps");
449 32 : keys.add("optional","WRITE_OV","write a file with model overlaps");
450 32 : keys.add("optional","AVERAGING", "Averaging window for weights");
451 32 : keys.add("optional","SEED","seed of the random number generator used by the Monte Carlo sampling (by default it is taken from the clock at every run)");
452 32 : keys.addFlag("NO_AVER",false,"don't do ensemble averaging in multi-replica mode");
453 32 : keys.addFlag("REWEIGHT",false,"simple REWEIGHT using the ARG as energy");
454 64 : keys.addOutputComponent("scoreb","default","scalar","Bayesian score");
455 64 : keys.addOutputComponent("acc", "NOISETYPE","scalar","MC acceptance for uncertainty");
456 64 : keys.addOutputComponent("scale", "REGRESSION","scalar","scale factor");
457 64 : keys.addOutputComponent("accscale", "REGRESSION","scalar","MC acceptance for scale regression");
458 64 : keys.addOutputComponent("enescale", "REGRESSION","scalar","MC energy for scale regression");
459 64 : keys.addOutputComponent("anneal","ANNEAL","scalar","annealing factor");
460 64 : keys.addOutputComponent("weight", "REWEIGHT","scalar", "weights of the weighted average");
461 64 : keys.addOutputComponent("biasDer", "REWEIGHT","scalar", "derivatives with respect to the bias");
462 64 : keys.addOutputComponent("sigma", "NOISETYPE","scalar", "uncertainty in the forward models and experiment");
463 64 : keys.addOutputComponent("neff", "default","scalar", "effective number of replicas");
464 32 : keys.addDOI("10.1101/113951");
465 32 : keys.addDOI("10.1126/sciadv.1501177");
466 32 : }
467 :
468 28 : EMMI::EMMI(const ActionOptions&ao):
469 : Action(ao),
470 : ActionAtomistic(ao),
471 : ActionWithArguments(ao),
472 : ActionWithValue(ao),
473 28 : inv_sqrt2_(0.707106781186548),
474 28 : sqrt2_pi_(0.797884560802865),
475 28 : first_time_(true), no_aver_(false), nl_begin_(0), nl_end_(0),
476 28 : valueScore_(NULL), valueAcc_(NULL), pbc_(true),
477 28 : MCstride_(1), MCaccept_(0.), MCtrials_(0.),
478 28 : statusstride_(0), first_status_(true),
479 28 : nregres_(0), scale_(1.),
480 28 : dpcutoff_(15.0), nexp_(1000000), nanneal_(0),
481 28 : kanneal_(0.), anneal_(1.), prior_(1.), ovstride_(0),
482 28 : do_reweight_(false), first_time_w_(true), decay_w_(1.) {
483 : // periodic boundary conditions
484 28 : bool nopbc=!pbc_;
485 28 : parseFlag("NOPBC",nopbc);
486 28 : pbc_=!nopbc;
487 :
488 : // list of atoms
489 28 : std::vector<AtomNumber> atoms;
490 56 : parseAtomList("ATOMS", atoms);
491 :
492 : // file with data GMM
493 : std::string GMM_file;
494 56 : parse("GMM_FILE", GMM_file);
495 :
496 : // type of data noise
497 : std::string noise;
498 28 : parse("NOISETYPE",noise);
499 28 : if (noise=="GAUSS") {
500 20 : noise_ = 0;
501 8 : } else if(noise=="OUTLIERS") {
502 5 : noise_ = 1;
503 3 : } else if(noise=="MARGINAL") {
504 3 : noise_ = 2;
505 : } else {
506 0 : error("Unknown noise type!");
507 : }
508 :
509 : // minimum value for error
510 : double sigma_min;
511 28 : parse("SIGMA_MIN", sigma_min);
512 28 : if(sigma_min<0) {
513 0 : error("SIGMA_MIN should be greater or equal to zero");
514 : }
515 :
516 : // the following parameters must be specified with noise type 0 and 1
517 28 : double sigma_ini=-1.0;
518 28 : double dsigma=-1.0;
519 28 : if(noise_!=2) {
520 : // initial value of the uncertainty
521 25 : parse("SIGMA0", sigma_ini);
522 25 : if(sigma_ini<=0) {
523 0 : error("you must specify a positive SIGMA0");
524 : }
525 : // MC parameters
526 25 : parse("DSIGMA", dsigma);
527 25 : if(dsigma<0) {
528 0 : error("you must specify a positive DSIGMA");
529 : }
530 25 : parse("MC_STRIDE", MCstride_);
531 25 : if(dsigma>0 && MCstride_<=0) {
532 0 : error("you must specify a positive MC_STRIDE");
533 : }
534 : // status file parameters
535 25 : parse("WRITE_STRIDE", statusstride_);
536 25 : if(statusstride_==0) {
537 0 : error("you must specify a positive WRITE_STRIDE");
538 : }
539 25 : parse("STATUS_FILE", statusfilename_);
540 25 : if(statusfilename_=="") {
541 50 : statusfilename_ = "MISTATUS"+getLabel();
542 : } else {
543 0 : statusfilename_ = statusfilename_+getLabel();
544 : }
545 : }
546 :
547 : // error file
548 : std::string errfile;
549 56 : parse("ERR_FILE", errfile);
550 :
551 : // file with experimental overlaps
552 : std::string ovfile;
553 28 : parse("OV_FILE", ovfile);
554 :
555 : // integral of the experimetal density
556 28 : double norm_d = 0.0;
557 28 : parse("NORM_DENSITY", norm_d);
558 :
559 : // temperature
560 28 : kbt_ = getkBT();
561 :
562 : // exponent of uncertainty prior
563 28 : parse("PRIOR",prior_);
564 :
565 : // simulated annealing stuff
566 28 : parse("ANNEAL", nanneal_);
567 28 : parse("ANNEAL_FACT", kanneal_);
568 28 : if(nanneal_>0 && kanneal_<=1.0) {
569 0 : error("with ANNEAL, ANNEAL_FACT must be greater than 1");
570 : }
571 :
572 : // regression stride
573 28 : parse("REGRESSION",nregres_);
574 : // other regression parameters
575 28 : if(nregres_>0) {
576 2 : parse("REG_SCALE_MIN",scale_min_);
577 2 : parse("REG_SCALE_MAX",scale_max_);
578 2 : parse("REG_DSCALE",dscale_);
579 : // checks
580 2 : if(scale_max_<=scale_min_) {
581 0 : error("with REGRESSION, REG_SCALE_MAX must be greater than REG_SCALE_MIN");
582 : }
583 2 : if(dscale_<=0.) {
584 0 : error("with REGRESSION, REG_DSCALE must be positive");
585 : }
586 : }
587 :
588 : // scale factor
589 28 : parse("SCALE", scale_);
590 :
591 : // read map resolution
592 : double reso;
593 28 : parse("RESOLUTION", reso);
594 28 : if(reso<=0.) {
595 0 : error("RESOLUTION should be strictly positive");
596 : }
597 :
598 : // neighbor list stuff
599 28 : parse("NL_CUTOFF",nl_cutoff_);
600 28 : if(nl_cutoff_<=0.0) {
601 0 : error("NL_CUTOFF should be explicitly specified and positive");
602 : }
603 28 : parse("NL_STRIDE",nl_stride_);
604 28 : if(nl_stride_==0) {
605 0 : error("NL_STRIDE should be explicitly specified and positive");
606 : }
607 :
608 : // averaging or not
609 28 : parseFlag("NO_AVER",no_aver_);
610 :
611 : // write overlap file
612 28 : parse("WRITE_OV_STRIDE", ovstride_);
613 28 : parse("WRITE_OV", ovfilename_);
614 28 : if(ovstride_>0 && ovfilename_=="") {
615 0 : error("With WRITE_OV_STRIDE you must specify WRITE_OV");
616 : }
617 :
618 : // set parallel stuff
619 28 : size_=comm.Get_size();
620 28 : rank_=comm.Get_rank();
621 :
622 : // get number of replicas
623 28 : if(no_aver_) {
624 18 : nrep_ = 1;
625 : } else {
626 10 : nrep_ = getNumberOfReplicas();
627 : }
628 28 : replica_ = getReplicaIndex();
629 :
630 : // Reweighting flag
631 28 : parseFlag("REWEIGHT", do_reweight_);
632 28 : if(do_reweight_&&getNumberOfArguments()!=1) {
633 0 : error("To REWEIGHT one must provide one single bias as an argument");
634 : }
635 28 : if(do_reweight_&&no_aver_) {
636 0 : error("REWEIGHT cannot be used with NO_AVER");
637 : }
638 28 : if(do_reweight_&&nrep_<2) {
639 0 : error("REWEIGHT can only be used in parallel with 2 or more replicas");
640 : }
641 28 : if(!getRestart()) {
642 28 : average_weights_.resize(nrep_, 1./static_cast<double>(nrep_));
643 : } else {
644 0 : average_weights_.resize(nrep_, 0.);
645 : }
646 :
647 28 : unsigned averaging=0;
648 28 : parse("AVERAGING", averaging);
649 28 : if(averaging>0) {
650 2 : decay_w_ = 1./static_cast<double> (averaging);
651 : }
652 :
653 28 : int seed=-1;
654 28 : parse("SEED",seed);
655 28 : if(seed<-1) {
656 0 : error("SEED must be a non-negative integer");
657 : }
658 :
659 28 : checkRead();
660 :
661 28 : log.printf(" atoms involved : ");
662 10964 : for(unsigned i=0; i<atoms.size(); ++i) {
663 10936 : log.printf("%d ",atoms[i].serial());
664 : }
665 28 : log.printf("\n");
666 28 : log.printf(" GMM data file : %s\n", GMM_file.c_str());
667 28 : if(no_aver_) {
668 18 : log.printf(" without ensemble averaging\n");
669 : }
670 28 : log.printf(" type of data noise : %s\n", noise.c_str());
671 28 : log.printf(" neighbor list cutoff : %lf\n", nl_cutoff_);
672 28 : log.printf(" neighbor list stride : %u\n", nl_stride_);
673 28 : log.printf(" minimum uncertainty : %f\n",sigma_min);
674 28 : log.printf(" scale factor : %lf\n",scale_);
675 28 : if(nregres_>0) {
676 2 : log.printf(" regression stride : %u\n", nregres_);
677 2 : log.printf(" regression minimum scale : %lf\n", scale_min_);
678 2 : log.printf(" regression maximum scale : %lf\n", scale_max_);
679 2 : log.printf(" regression maximum scale MC move : %lf\n", dscale_);
680 : }
681 28 : if(noise_!=2) {
682 25 : log.printf(" initial value of the uncertainty : %f\n",sigma_ini);
683 25 : log.printf(" max MC move in uncertainty : %f\n",dsigma);
684 25 : log.printf(" MC stride : %d\n", MCstride_);
685 25 : log.printf(" reading/writing to status file : %s\n",statusfilename_.c_str());
686 25 : log.printf(" with stride : %u\n",statusstride_);
687 : }
688 28 : if(errfile.size()>0) {
689 2 : log.printf(" reading experimental errors from file : %s\n", errfile.c_str());
690 : }
691 28 : if(ovfile.size()>0) {
692 2 : log.printf(" reading experimental overlaps from file : %s\n", ovfile.c_str());
693 : }
694 28 : log.printf(" temperature of the system in energy unit : %f\n",kbt_);
695 28 : log.printf(" prior exponent : %f\n",prior_);
696 28 : log.printf(" number of replicas for averaging: %u\n",nrep_);
697 28 : log.printf(" id of the replica : %u\n",replica_);
698 28 : if(nanneal_>0) {
699 4 : log.printf(" length of annealing cycle : %u\n",nanneal_);
700 4 : log.printf(" annealing factor : %f\n",kanneal_);
701 : }
702 28 : if(ovstride_>0) {
703 4 : log.printf(" stride for writing model overlaps : %u\n",ovstride_);
704 4 : log.printf(" file for writing model overlaps : %s\n", ovfilename_.c_str());
705 : }
706 :
707 : // set constant quantity before calculating stuff
708 28 : cfact_ = 1.0/pow( 2.0*pi, 1.5 );
709 :
710 : // calculate model GMM constant parameters
711 28 : std::vector<double> GMM_m_w = get_GMM_m(atoms);
712 :
713 : // read data GMM parameters
714 28 : get_GMM_d(GMM_file);
715 28 : log.printf(" number of GMM components : %u\n", static_cast<unsigned>(GMM_d_m_.size()));
716 :
717 : // normalize atom weight map
718 28 : if(norm_d <= 0.0) {
719 28 : norm_d = accumulate(GMM_d_w_.begin(), GMM_d_w_.end(), 0.0);
720 : }
721 : double norm_m = accumulate(GMM_m_w.begin(), GMM_m_w.end(), 0.0);
722 : // renormalization
723 140 : for(unsigned i=0; i<GMM_m_w_.size(); ++i) {
724 112 : GMM_m_w_[i] *= norm_d / norm_m;
725 : }
726 :
727 : // read experimental errors
728 28 : std::vector<double> exp_err;
729 28 : if(errfile.size()>0) {
730 4 : exp_err = read_exp_errors(errfile);
731 : }
732 :
733 : // get self overlaps between data GMM components
734 28 : if(ovfile.size()>0) {
735 4 : ovdd_ = read_exp_overlaps(ovfile);
736 : } else {
737 7388 : for(unsigned i=0; i<GMM_d_m_.size(); ++i) {
738 7362 : double ov = get_self_overlap(i);
739 7362 : ovdd_.push_back(ov);
740 : }
741 : }
742 :
743 28 : log.printf(" number of GMM groups : %u\n", static_cast<unsigned>(GMM_d_grps_.size()));
744 : // cycle on GMM groups
745 64 : for(unsigned Gid=0; Gid<GMM_d_grps_.size(); ++Gid) {
746 36 : log.printf(" group %d\n", Gid);
747 : // calculate median overlap and experimental error
748 36 : std::vector<double> ovdd;
749 36 : std::vector<double> err;
750 : // cycle on the group members
751 9198 : for(unsigned i=0; i<GMM_d_grps_[Gid].size(); ++i) {
752 : // GMM id
753 9162 : int GMMid = GMM_d_grps_[Gid][i];
754 : // add to experimental error
755 9162 : if(errfile.size()>0) {
756 1800 : err.push_back(exp_err[GMMid]);
757 : } else {
758 7362 : err.push_back(0.);
759 : }
760 : // add to GMM overlap
761 9162 : ovdd.push_back(ovdd_[GMMid]);
762 : }
763 : // calculate median quantities
764 36 : double ovdd_m = get_median(ovdd);
765 36 : double err_m = get_median(err);
766 : // print out statistics
767 36 : log.printf(" # of members : %zu\n", GMM_d_grps_[Gid].size());
768 36 : log.printf(" median overlap : %lf\n", ovdd_m);
769 36 : log.printf(" median error : %lf\n", err_m);
770 : // add minimum value of sigma for this group of GMMs
771 36 : sigma_min_.push_back(std::sqrt(err_m*err_m+sigma_min*ovdd_m*sigma_min*ovdd_m));
772 : // these are only needed with Gaussian and Outliers noise models
773 36 : if(noise_!=2) {
774 : // set dsigma
775 33 : dsigma_.push_back(dsigma * ovdd_m);
776 : // set sigma max
777 33 : sigma_max_.push_back(10.0*ovdd_m + sigma_min_[Gid] + dsigma_[Gid]);
778 : // initialize sigma
779 66 : sigma_.push_back(std::max(sigma_min_[Gid],std::min(sigma_ini*ovdd_m,sigma_max_[Gid])));
780 : }
781 : }
782 :
783 : // read status file if restarting
784 28 : if(getRestart() && noise_!=2) {
785 0 : read_status();
786 : }
787 :
788 : // calculate auxiliary stuff
789 28 : calculate_useful_stuff(reso);
790 :
791 : // prepare data and derivative std::vectors
792 28 : ovmd_.resize(ovdd_.size());
793 28 : ovmd_ave_.resize(ovdd_.size());
794 28 : atom_der_.resize(GMM_m_type_.size());
795 28 : GMMid_der_.resize(ovdd_.size());
796 :
797 : // clear things that are no longer needed
798 28 : GMM_d_cov_.clear();
799 :
800 : // add components
801 56 : addComponentWithDerivatives("scoreb");
802 28 : componentIsNotPeriodic("scoreb");
803 28 : valueScore_=getPntrToComponent("scoreb");
804 :
805 28 : if(noise_!=2) {
806 50 : addComponent("acc");
807 25 : componentIsNotPeriodic("acc");
808 25 : valueAcc_=getPntrToComponent("acc");
809 : }
810 :
811 28 : if(nregres_>0) {
812 4 : addComponent("scale");
813 4 : componentIsNotPeriodic("scale");
814 4 : addComponent("accscale");
815 4 : componentIsNotPeriodic("accscale");
816 4 : addComponent("enescale");
817 4 : componentIsNotPeriodic("enescale");
818 : }
819 :
820 28 : if(nanneal_>0) {
821 8 : addComponent("anneal");
822 8 : componentIsNotPeriodic("anneal");
823 : }
824 :
825 28 : if(do_reweight_) {
826 12 : addComponent("biasDer");
827 12 : componentIsNotPeriodic("biasDer");
828 12 : addComponent("weight");
829 12 : componentIsNotPeriodic("weight");
830 : }
831 :
832 56 : addComponent("neff");
833 28 : componentIsNotPeriodic("neff");
834 :
835 61 : for(unsigned i=0; i<sigma_.size(); ++i) {
836 : std::string num;
837 33 : Tools::convert(i, num);
838 66 : addComponent("sigma-"+num);
839 33 : componentIsNotPeriodic("sigma-"+num);
840 66 : valueSigma_.push_back(getPntrToComponent("sigma-"+num));
841 33 : valueSigma_.back()->set(sigma_[i]);
842 : }
843 :
844 : // initialize random seed: the Monte Carlo moves are done by the first rank of each replica.
845 : // With SEED replica r gets the stream 1+(SEED*nreplicas+r) mod (2^31-1)
846 : unsigned iseed;
847 28 : if(seed>=0) {
848 8 : const unsigned long long nr = getNumberOfReplicas();
849 8 : iseed = static_cast<unsigned>(1 + (static_cast<unsigned long long>(seed)*nr+replica_) % 2147483646ULL);
850 20 : } else if(rank_==0) {
851 12 : iseed = time(NULL)+replica_;
852 : } else {
853 8 : iseed = 0;
854 : }
855 28 : if(seed<0) {
856 20 : comm.Sum(&iseed, 1);
857 : }
858 28 : random_.setSeed(-static_cast<int>(iseed));
859 :
860 : // request the atoms
861 28 : requestAtoms(atoms);
862 28 : setDerivatives();
863 :
864 : // print bibliography
865 28 : log<<" Bibliography "<<plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)");
866 28 : log<<plumed.cite("Hanot, Bonomi, Greenberg, Sali, Nilges, Vendruscolo, Pellarin, bioRxiv doi: 10.1101/113951 (2017)");
867 28 : log<<plumed.cite("Bonomi, Pellarin, Vendruscolo, Biophys. J. 114, 1604 (2018)");
868 28 : if(do_reweight_) {
869 12 : log<<plumed.cite("Bonomi, Camilloni, Vendruscolo, Sci. Rep. 6, 31232 (2016)");
870 : }
871 28 : if(!no_aver_ && nrep_>1) {
872 12 : log<<plumed.cite("Bonomi, Camilloni, Cavalli, Vendruscolo, Sci. Adv. 2, e150117 (2016)");
873 : }
874 28 : log<<"\n";
875 28 : }
876 :
877 6 : void EMMI::write_model_overlap(long long int step) {
878 6 : OFile ovfile;
879 6 : ovfile.link(*this);
880 : std::string num;
881 6 : Tools::convert(step,num);
882 6 : std::string name = ovfilename_+"-"+num;
883 6 : ovfile.open(name);
884 : ovfile.setHeavyFlush();
885 6 : ovfile.fmtField("%10.7e ");
886 : // write overlaps
887 1842 : for(unsigned i=0; i<ovmd_.size(); ++i) {
888 1836 : ovfile.printField("Model", ovmd_[i]);
889 3672 : ovfile.printField("ModelScaled", scale_ * ovmd_[i]);
890 1836 : ovfile.printField("Data", ovdd_[i]);
891 1836 : ovfile.printField();
892 : }
893 6 : ovfile.close();
894 6 : }
895 :
896 72 : double EMMI::get_median(std::vector<double> &v) {
897 : // dimension of std::vector
898 72 : unsigned size = v.size();
899 : // in case of only one entry
900 72 : if (size==1) {
901 0 : return v[0];
902 : } else {
903 : // reorder std::vector
904 72 : sort(v.begin(), v.end());
905 : // odd or even?
906 72 : if (size%2==0) {
907 36 : return (v[size/2-1]+v[size/2])/2.0;
908 : } else {
909 36 : return v[size/2];
910 : }
911 : }
912 : }
913 :
914 0 : void EMMI::read_status() {
915 : double MDtime;
916 : // open file
917 : auto ifile = Tools::make_unique<IFile>();
918 0 : ifile->link(*this);
919 0 : if(ifile->FileExist(statusfilename_)) {
920 0 : ifile->open(statusfilename_);
921 0 : while(ifile->scanField("MD_time", MDtime)) {
922 0 : for(unsigned i=0; i<sigma_.size(); ++i) {
923 : // convert i to std::string
924 : std::string num;
925 0 : Tools::convert(i,num);
926 : // read entries
927 0 : ifile->scanField("s"+num, sigma_[i]);
928 : }
929 : // new line
930 0 : ifile->scanField();
931 : }
932 0 : ifile->close();
933 : } else {
934 0 : error("Cannot find status file "+statusfilename_+"\n");
935 : }
936 0 : }
937 :
938 25 : void EMMI::print_status(long long int step) {
939 : // if first time open the file
940 25 : if(first_status_) {
941 25 : first_status_ = false;
942 25 : statusfile_.link(*this);
943 25 : statusfile_.open(statusfilename_);
944 : statusfile_.setHeavyFlush();
945 50 : statusfile_.fmtField("%6.3e ");
946 : }
947 : // write fields
948 25 : double MDtime = static_cast<double>(step)*getTimeStep();
949 25 : statusfile_.printField("MD_time", MDtime);
950 58 : for(unsigned i=0; i<sigma_.size(); ++i) {
951 : // convert i to std::string
952 : std::string num;
953 33 : Tools::convert(i,num);
954 : // print entry
955 66 : statusfile_.printField("s"+num, sigma_[i]);
956 : }
957 25 : statusfile_.printField();
958 25 : }
959 :
960 200012 : bool EMMI::doAccept(double oldE, double newE, double kbt) {
961 : bool accept = false;
962 : // calculate delta energy
963 200012 : double delta = ( newE - oldE ) / kbt;
964 : // if delta is negative always accept move
965 200012 : if( delta < 0.0 ) {
966 : accept = true;
967 : } else {
968 : // otherwise extract random number
969 171964 : double s = random_.RandU01();
970 171964 : if( s < std::exp(-delta) ) {
971 : accept = true;
972 : }
973 : }
974 200012 : return accept;
975 : }
976 :
977 24 : void EMMI::doMonteCarlo() {
978 : // the move is done by the first rank and then broadcast
979 24 : if(rank_==0) {
980 12 : doMonteCarloMove();
981 : }
982 : // increment number of trials
983 24 : MCtrials_ += 1.0;
984 24 : if(size_>1) {
985 24 : comm.Bcast(sigma_,0);
986 24 : comm.Bcast(MCaccept_,0);
987 : }
988 : // update sigma output
989 72 : for(unsigned i=0; i<sigma_.size(); ++i) {
990 48 : valueSigma_[i]->set(sigma_[i]);
991 : }
992 24 : }
993 :
994 12 : void EMMI::doMonteCarloMove() {
995 : // extract random GMM group
996 12 : unsigned nGMM = static_cast<unsigned>(std::floor(random_.RandU01()*static_cast<double>(GMM_d_grps_.size())));
997 12 : if(nGMM==GMM_d_grps_.size()) {
998 0 : nGMM -= 1;
999 : }
1000 :
1001 : // generate random move
1002 12 : double shift = dsigma_[nGMM] * ( 2.0 * random_.RandU01() - 1.0 );
1003 : // new sigma
1004 12 : double new_s = sigma_[nGMM] + shift;
1005 : // check boundaries
1006 12 : if(new_s > sigma_max_[nGMM]) {
1007 0 : new_s = 2.0 * sigma_max_[nGMM] - new_s;
1008 : }
1009 12 : if(new_s < sigma_min_[nGMM]) {
1010 0 : new_s = 2.0 * sigma_min_[nGMM] - new_s;
1011 : }
1012 : // old s2
1013 12 : double old_inv_s2 = 1.0 / sigma_[nGMM] / sigma_[nGMM];
1014 : // new s2
1015 12 : double new_inv_s2 = 1.0 / new_s / new_s;
1016 :
1017 : // cycle on GMM group and calculate old and new energy
1018 : double old_ene = 0.0;
1019 : double new_ene = 0.0;
1020 12 : double ng = static_cast<double>(GMM_d_grps_[nGMM].size());
1021 :
1022 : // in case of Gaussian noise
1023 12 : if(noise_==0) {
1024 : double chi2 = 0.0;
1025 2706 : for(unsigned i=0; i<GMM_d_grps_[nGMM].size(); ++i) {
1026 : // id GMM component
1027 2700 : int GMMid = GMM_d_grps_[nGMM][i];
1028 : // deviation
1029 2700 : double dev = ( scale_*ovmd_ave_[GMMid]-ovdd_[GMMid] );
1030 : // add to chi2
1031 2700 : chi2 += dev * dev;
1032 : }
1033 : // final energy calculation: add normalization and prior
1034 6 : old_ene = 0.5 * kbt_ * ( chi2 * old_inv_s2 - (ng+prior_) * std::log(old_inv_s2) );
1035 6 : new_ene = 0.5 * kbt_ * ( chi2 * new_inv_s2 - (ng+prior_) * std::log(new_inv_s2) );
1036 : }
1037 :
1038 : // in case of Outliers noise
1039 12 : if(noise_==1) {
1040 2706 : for(unsigned i=0; i<GMM_d_grps_[nGMM].size(); ++i) {
1041 : // id GMM component
1042 2700 : int GMMid = GMM_d_grps_[nGMM][i];
1043 : // calculate deviation
1044 2700 : double dev = ( scale_*ovmd_ave_[GMMid]-ovdd_[GMMid] );
1045 : // add to energies
1046 2700 : old_ene += std::log1p( 0.5 * dev * dev * old_inv_s2);
1047 2700 : new_ene += std::log1p( 0.5 * dev * dev * new_inv_s2);
1048 : }
1049 : // final energy calculation: add normalization and prior
1050 6 : old_ene = kbt_ * ( old_ene + (ng+prior_) * std::log(sigma_[nGMM]) );
1051 6 : new_ene = kbt_ * ( new_ene + (ng+prior_) * std::log(new_s) );
1052 : }
1053 :
1054 : // accept or reject
1055 12 : bool accept = doAccept(old_ene/anneal_, new_ene/anneal_, kbt_);
1056 12 : if(accept) {
1057 5 : sigma_[nGMM] = new_s;
1058 5 : MCaccept_ += 1.0;
1059 : }
1060 12 : }
1061 :
1062 2 : std::vector<double> EMMI::read_exp_errors(const std::string & errfile) {
1063 : int nexp, idcomp;
1064 : double err;
1065 2 : std::vector<double> exp_err;
1066 : // open file
1067 2 : IFile *ifile = new IFile();
1068 2 : if(ifile->FileExist(errfile)) {
1069 2 : ifile->open(errfile);
1070 : // scan for number of experimental errors
1071 2 : ifile->scanField("Nexp", nexp);
1072 : // cycle on GMM components
1073 3604 : while(ifile->scanField("Id",idcomp)) {
1074 : // total experimental error
1075 1800 : double err_tot = 0.0;
1076 : // cycle on number of experimental overlaps
1077 5400 : for(int i=0; i<nexp; ++i) {
1078 : std::string ss;
1079 3600 : Tools::convert(i,ss);
1080 3600 : ifile->scanField("Err"+ss, err);
1081 : // add to total error
1082 3600 : err_tot += err*err;
1083 : }
1084 : // new line
1085 1800 : ifile->scanField();
1086 : // calculate RMSE
1087 1800 : err_tot = std::sqrt(err_tot/static_cast<double>(nexp));
1088 : // add to global
1089 1800 : exp_err.push_back(err_tot);
1090 : }
1091 2 : ifile->close();
1092 : } else {
1093 0 : error("Cannot find ERR_FILE "+errfile+"\n");
1094 : }
1095 2 : return exp_err;
1096 : }
1097 :
1098 2 : std::vector<double> EMMI::read_exp_overlaps(const std::string & ovfile) {
1099 : int idcomp;
1100 : double ov;
1101 2 : std::vector<double> ovdd;
1102 : // open file
1103 2 : IFile *ifile = new IFile();
1104 2 : if(ifile->FileExist(ovfile)) {
1105 2 : ifile->open(ovfile);
1106 : // cycle on GMM components
1107 3604 : while(ifile->scanField("Id",idcomp)) {
1108 : // read experimental overlap
1109 1800 : ifile->scanField("Overlap", ov);
1110 : // add to ovdd
1111 1800 : ovdd.push_back(ov);
1112 : // new line
1113 1800 : ifile->scanField();
1114 : }
1115 2 : ifile->close();
1116 : } else {
1117 0 : error("Cannot find OV_FILE "+ovfile+"\n");
1118 : }
1119 2 : return ovdd;
1120 : }
1121 :
1122 28 : std::vector<double> EMMI::get_GMM_m(std::vector<AtomNumber> &atoms) {
1123 : // list of weights - one per atom
1124 28 : std::vector<double> GMM_m_w;
1125 :
1126 28 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
1127 : // map of atom types to A and B coefficients of scattering factor
1128 : // f(s) = A * exp(-B*s**2)
1129 : // B is in Angstrom squared
1130 : // map between an atom type and an index
1131 : std::map<std::string, unsigned> type_map;
1132 28 : type_map["C"]=0;
1133 28 : type_map["O"]=1;
1134 28 : type_map["N"]=2;
1135 28 : type_map["S"]=3;
1136 : // fill in sigma std::vector
1137 28 : GMM_m_s_.push_back(0.01*15.146); // type 0
1138 28 : GMM_m_s_.push_back(0.01*8.59722); // type 1
1139 28 : GMM_m_s_.push_back(0.01*11.1116); // type 2
1140 28 : GMM_m_s_.push_back(0.01*15.8952); // type 3
1141 : // fill in weight std::vector
1142 28 : GMM_m_w_.push_back(2.49982); // type 0
1143 28 : GMM_m_w_.push_back(1.97692); // type 1
1144 28 : GMM_m_w_.push_back(2.20402); // type 2
1145 28 : GMM_m_w_.push_back(5.14099); // type 3
1146 :
1147 : // check if MOLINFO line is present
1148 28 : if( moldat ) {
1149 28 : log<<" MOLINFO DATA found with label " <<moldat->getLabel()<<", using proper atom names\n";
1150 10964 : for(unsigned i=0; i<atoms.size(); ++i) {
1151 : // get atom name
1152 10936 : std::string name = moldat->getAtomName(atoms[i]);
1153 : char type;
1154 : // get atom type
1155 10936 : char first = name.at(0);
1156 : // GOLDEN RULE: type is first letter, if not a number
1157 10936 : if (!isdigit(first)) {
1158 : type = first;
1159 : // otherwise is the second
1160 : } else {
1161 0 : type = name.at(1);
1162 : }
1163 : // check if key in map
1164 : std::string type_s = std::string(1,type);
1165 10936 : if(type_map.find(type_s) != type_map.end()) {
1166 : // save atom type
1167 10936 : GMM_m_type_.push_back(type_map[type_s]);
1168 : // this will be normalized in the final density
1169 10936 : GMM_m_w.push_back(GMM_m_w_[type_map[type_s]]);
1170 : } else {
1171 0 : error("Wrong atom type "+type_s+" from atom name "+name+"\n");
1172 : }
1173 : }
1174 : } else {
1175 0 : error("MOLINFO DATA not found\n");
1176 : }
1177 28 : return GMM_m_w;
1178 : }
1179 :
1180 9162 : void EMMI::check_GMM_d(const VectorGeneric<6> &cov, double w) {
1181 :
1182 : // check if positive defined, by calculating the 3 leading principal minors
1183 9162 : const double pm1 = cov[0];
1184 9162 : const double pm2 = cov[0]*cov[3]-cov[1]*cov[1];
1185 9162 : const double pm3 = cov[0]*(cov[3]*cov[5]-cov[4]*cov[4])-cov[1]*(cov[1]*cov[5]-cov[4]*cov[2])+cov[2]*(cov[1]*cov[4]-cov[3]*cov[2]);
1186 : // apply Sylvester’s criterion
1187 9162 : if(pm1<=0.0 || pm2<=0.0 || pm3<=0.0) {
1188 0 : error("check data GMM: covariance matrix is not positive defined");
1189 : }
1190 :
1191 : // check if weight is positive
1192 9162 : if(w<=0) {
1193 0 : error("check data GMM: weight must be positive");
1194 : }
1195 9162 : }
1196 :
1197 : // read GMM data file in PLUMED format:
1198 28 : void EMMI::get_GMM_d(const std::string & GMM_file) {
1199 28 : VectorGeneric<6> cov;
1200 :
1201 : // open file
1202 : auto ifile=Tools::make_unique<IFile>();
1203 28 : if(ifile->FileExist(GMM_file)) {
1204 28 : ifile->open(GMM_file);
1205 : int idcomp;
1206 18380 : while(ifile->scanField("Id",idcomp)) {
1207 : int beta;
1208 : double w, m0, m1, m2;
1209 9162 : ifile->scanField("Weight",w);
1210 9162 : ifile->scanField("Mean_0",m0);
1211 9162 : ifile->scanField("Mean_1",m1);
1212 9162 : ifile->scanField("Mean_2",m2);
1213 9162 : ifile->scanField("Cov_00",cov[0]);
1214 9162 : ifile->scanField("Cov_01",cov[1]);
1215 9162 : ifile->scanField("Cov_02",cov[2]);
1216 9162 : ifile->scanField("Cov_11",cov[3]);
1217 9162 : ifile->scanField("Cov_12",cov[4]);
1218 9162 : ifile->scanField("Cov_22",cov[5]);
1219 9162 : ifile->scanField("Beta",beta);
1220 : // check input
1221 9162 : check_GMM_d(cov, w);
1222 : // check beta
1223 9162 : if(beta<0) {
1224 0 : error("Beta must be positive!");
1225 : }
1226 : // center of the Gaussian
1227 9162 : GMM_d_m_.push_back(Vector(m0,m1,m2));
1228 : // covariance matrix
1229 9162 : GMM_d_cov_.push_back(cov);
1230 : // weight
1231 9162 : GMM_d_w_.push_back(w);
1232 : // beta
1233 9162 : GMM_d_beta_.push_back(beta);
1234 : // new line
1235 9162 : ifile->scanField();
1236 : }
1237 : } else {
1238 0 : error("Cannot find GMM_FILE "+GMM_file+"\n");
1239 : }
1240 : // now create a set from beta (unique set of values)
1241 28 : std::set<int> bu(GMM_d_beta_.begin(), GMM_d_beta_.end());
1242 : // now prepare the group std::vector
1243 28 : GMM_d_grps_.resize(bu.size());
1244 : // and fill it in
1245 9190 : for(unsigned i=0; i<GMM_d_beta_.size(); ++i) {
1246 : //I do not know if GMM_d_beta is always >0, if I can assumo that GMM_d_beta is >=0, i can skip the signedness check
1247 9162 : if((GMM_d_beta_[i]>=0) &&
1248 9162 : (static_cast<unsigned>(GMM_d_beta_[i])>=GMM_d_grps_.size())) {
1249 0 : error("Check Beta values");
1250 : }
1251 9162 : GMM_d_grps_[GMM_d_beta_[i]].push_back(i);
1252 : }
1253 28 : }
1254 :
1255 28 : void EMMI::calculate_useful_stuff(double reso) {
1256 : // We use the following definition for resolution:
1257 : // the Fourier transform of the density distribution in real space
1258 : // f(s) falls to 1/e of its maximum value at wavenumber 1/resolution
1259 : // i.e. from f(s) = A * exp(-B*s**2) -> Res = std::sqrt(B).
1260 : // average value of B
1261 : double Bave = 0.0;
1262 10964 : for(unsigned i=0; i<GMM_m_type_.size(); ++i) {
1263 10936 : Bave += GMM_m_s_[GMM_m_type_[i]];
1264 : }
1265 28 : Bave /= static_cast<double>(GMM_m_type_.size());
1266 : // calculate blur factor
1267 : double blur = 0.0;
1268 28 : if(reso*reso>Bave) {
1269 10 : blur = reso*reso-Bave;
1270 : } else {
1271 36 : warning("PLUMED should not be used with maps at resolution better than 0.3 nm");
1272 : }
1273 : // add blur to B
1274 140 : for(unsigned i=0; i<GMM_m_s_.size(); ++i) {
1275 112 : GMM_m_s_[i] += blur;
1276 : }
1277 : // calculate average resolution
1278 : double ave_res = 0.0;
1279 10964 : for(unsigned i=0; i<GMM_m_type_.size(); ++i) {
1280 10936 : ave_res += std::sqrt(GMM_m_s_[GMM_m_type_[i]]);
1281 : }
1282 28 : ave_res = ave_res / static_cast<double>(GMM_m_type_.size());
1283 28 : log.printf(" experimental map resolution : %3.2f\n", reso);
1284 28 : log.printf(" predicted map resolution : %3.2f\n", ave_res);
1285 28 : log.printf(" blur factor : %f\n", blur);
1286 : // now calculate useful stuff
1287 28 : VectorGeneric<6> cov, sum, inv_sum;
1288 : // cycle on all atoms types (4 for the moment)
1289 140 : for(unsigned i=0; i<GMM_m_s_.size(); ++i) {
1290 : // the Gaussian in density (real) space is the FT of scattering factor
1291 : // f(r) = A * (pi/B)**1.5 * exp(-pi**2/B*r**2)
1292 112 : double s = std::sqrt ( 0.5 * GMM_m_s_[i] ) / pi;
1293 : // covariance matrix for spherical Gaussian
1294 112 : cov[0]=s*s;
1295 : cov[1]=0.0;
1296 : cov[2]=0.0;
1297 112 : cov[3]=s*s;
1298 : cov[4]=0.0;
1299 112 : cov[5]=s*s;
1300 : // cycle on all data GMM
1301 36760 : for(unsigned j=0; j<GMM_d_m_.size(); ++j) {
1302 : // we need the sum of the covariance matrices
1303 256536 : for(unsigned k=0; k<6; ++k) {
1304 219888 : sum[k] = cov[k] + GMM_d_cov_[j][k];
1305 : }
1306 : // and to calculate its determinant
1307 36648 : double det = sum[0]*(sum[3]*sum[5]-sum[4]*sum[4]);
1308 36648 : det -= sum[1]*(sum[1]*sum[5]-sum[4]*sum[2]);
1309 36648 : det += sum[2]*(sum[1]*sum[4]-sum[3]*sum[2]);
1310 : // calculate prefactor - model weights are already normalized
1311 36648 : double pre_fact = cfact_ / std::sqrt(det) * GMM_d_w_[j] * GMM_m_w_[i];
1312 : // and its inverse
1313 36648 : inv_sum[0] = (sum[3]*sum[5] - sum[4]*sum[4])/det;
1314 36648 : inv_sum[1] = (sum[2]*sum[4] - sum[1]*sum[5])/det;
1315 36648 : inv_sum[2] = (sum[1]*sum[4] - sum[2]*sum[3])/det;
1316 36648 : inv_sum[3] = (sum[0]*sum[5] - sum[2]*sum[2])/det;
1317 36648 : inv_sum[4] = (sum[2]*sum[1] - sum[0]*sum[4])/det;
1318 36648 : inv_sum[5] = (sum[0]*sum[3] - sum[1]*sum[1])/det;
1319 : // now we store the prefactor
1320 36648 : pre_fact_.push_back(pre_fact);
1321 : // and the inverse of the sum
1322 36648 : inv_cov_md_.push_back(inv_sum);
1323 : }
1324 : }
1325 : // tabulate exponential
1326 28 : dexp_ = dpcutoff_ / static_cast<double> (nexp_-1);
1327 28000028 : for(unsigned i=0; i<nexp_; ++i) {
1328 28000000 : tab_exp_.push_back(std::exp(-static_cast<double>(i) * dexp_));
1329 : }
1330 28 : }
1331 :
1332 : // get prefactors
1333 6481458 : double EMMI::get_prefactor_inverse
1334 : (const VectorGeneric<6> &GMM_cov_0, const VectorGeneric<6> &GMM_cov_1,
1335 : double GMM_w_0, double GMM_w_1,
1336 : VectorGeneric<6> &sum, VectorGeneric<6> &inv_sum) {
1337 : // we need the sum of the covariance matrices
1338 45370206 : for(unsigned k=0; k<6; ++k) {
1339 38888748 : sum[k] = GMM_cov_0[k] + GMM_cov_1[k];
1340 : }
1341 :
1342 : // and to calculate its determinant
1343 6481458 : double det = sum[0]*(sum[3]*sum[5]-sum[4]*sum[4]);
1344 6481458 : det -= sum[1]*(sum[1]*sum[5]-sum[4]*sum[2]);
1345 6481458 : det += sum[2]*(sum[1]*sum[4]-sum[3]*sum[2]);
1346 :
1347 : // the prefactor is
1348 6481458 : double pre_fact = cfact_ / std::sqrt(det) * GMM_w_0 * GMM_w_1;
1349 :
1350 : // and its inverse
1351 6481458 : inv_sum[0] = (sum[3]*sum[5] - sum[4]*sum[4])/det;
1352 6481458 : inv_sum[1] = (sum[2]*sum[4] - sum[1]*sum[5])/det;
1353 6481458 : inv_sum[2] = (sum[1]*sum[4] - sum[2]*sum[3])/det;
1354 6481458 : inv_sum[3] = (sum[0]*sum[5] - sum[2]*sum[2])/det;
1355 6481458 : inv_sum[4] = (sum[2]*sum[1] - sum[0]*sum[4])/det;
1356 6481458 : inv_sum[5] = (sum[0]*sum[3] - sum[1]*sum[1])/det;
1357 :
1358 : // return pre-factor
1359 6481458 : return pre_fact;
1360 : }
1361 :
1362 7362 : double EMMI::get_self_overlap(unsigned id) {
1363 : double ov_tot = 0.0;
1364 7362 : VectorGeneric<6> sum, inv_sum;
1365 7362 : Vector ov_der;
1366 : // start loop
1367 6488820 : for(unsigned i=0; i<GMM_d_m_.size(); ++i) {
1368 : // call auxiliary method
1369 6481458 : double pre_fact = get_prefactor_inverse(GMM_d_cov_[id], GMM_d_cov_[i],
1370 6481458 : GMM_d_w_[id], GMM_d_w_[i], sum, inv_sum);
1371 : // add overlap to ov_tot
1372 6481458 : ov_tot += get_overlap(GMM_d_m_[id], GMM_d_m_[i], pre_fact, inv_sum, ov_der);
1373 : }
1374 : // and return it
1375 7362 : return ov_tot;
1376 : }
1377 :
1378 : // get overlap and derivatives
1379 6634067 : double EMMI::get_overlap(const Vector &m_m, const Vector &d_m, double pre_fact,
1380 : const VectorGeneric<6> &inv_cov_md, Vector &ov_der) {
1381 : Vector md;
1382 : // calculate std::vector difference m_m-d_m with/without pbc
1383 6634067 : if(pbc_) {
1384 0 : md = pbcDistance(d_m, m_m);
1385 : } else {
1386 : md = delta(d_m, m_m);
1387 : }
1388 : // calculate product of transpose of md and inv_cov_md
1389 6634067 : double p_x = md[0]*inv_cov_md[0]+md[1]*inv_cov_md[1]+md[2]*inv_cov_md[2];
1390 6634067 : double p_y = md[0]*inv_cov_md[1]+md[1]*inv_cov_md[3]+md[2]*inv_cov_md[4];
1391 6634067 : double p_z = md[0]*inv_cov_md[2]+md[1]*inv_cov_md[4]+md[2]*inv_cov_md[5];
1392 : // calculate product of prod and md
1393 6634067 : double ov = md[0]*p_x+md[1]*p_y+md[2]*p_z;
1394 : // final calculation
1395 6634067 : ov = pre_fact * std::exp(-0.5*ov);
1396 : // derivatives
1397 6634067 : ov_der = ov * Vector(p_x, p_y, p_z);
1398 6634067 : return ov;
1399 : }
1400 :
1401 : // get the exponent of the overlap
1402 144180 : double EMMI::get_exp_overlap(const Vector &m_m, const Vector &d_m,
1403 : const VectorGeneric<6> &inv_cov_md) {
1404 : Vector md;
1405 : // calculate std::vector difference m_m-d_m with/without pbc
1406 144180 : if(pbc_) {
1407 0 : md = pbcDistance(d_m, m_m);
1408 : } else {
1409 : md = delta(d_m, m_m);
1410 : }
1411 : // calculate product of transpose of md and inv_cov_md
1412 144180 : double p_x = md[0]*inv_cov_md[0]+md[1]*inv_cov_md[1]+md[2]*inv_cov_md[2];
1413 144180 : double p_y = md[0]*inv_cov_md[1]+md[1]*inv_cov_md[3]+md[2]*inv_cov_md[4];
1414 144180 : double p_z = md[0]*inv_cov_md[2]+md[1]*inv_cov_md[4]+md[2]*inv_cov_md[5];
1415 : // calculate product of prod and md
1416 144180 : double ov = md[0]*p_x+md[1]*p_y+md[2]*p_z;
1417 144180 : return ov;
1418 : }
1419 :
1420 36 : void EMMI::update_neighbor_list() {
1421 : // dimension of GMM and atom std::vectors
1422 36 : const unsigned GMM_d_size = GMM_d_m_.size();
1423 36 : const unsigned GMM_m_size = GMM_m_type_.size();
1424 : // each rank does a contiguous block of GMM components, so that
1425 : // the neighbor list is sorted by GMM component
1426 36 : const unsigned id_begin = static_cast<unsigned>((static_cast<unsigned long long>(GMM_d_size)*rank_)/size_);
1427 36 : const unsigned id_end = static_cast<unsigned>((static_cast<unsigned long long>(GMM_d_size)*(rank_+1))/size_);
1428 36 : const unsigned nt = OpenMP::getNumThreads();
1429 36 : std::vector< std::vector<unsigned> > nl_t(nt);
1430 :
1431 36 : #pragma omp parallel num_threads(nt)
1432 : {
1433 : std::vector<unsigned> & nl_l = nl_t[OpenMP::getThreadNum()];
1434 : // overlaps and atoms (a map keyed by the overlap lost the atoms with the same tabulated overlap)
1435 : std::vector< std::pair<double,unsigned> > ov_l;
1436 : #pragma omp for schedule(static)
1437 : for(unsigned id=id_begin; id<id_end; ++id) {
1438 : ov_l.clear();
1439 : // total overlap with id
1440 : double ov_tot = 0.0;
1441 : // cycle on all atoms
1442 : for(unsigned im=0; im<GMM_m_size; ++im) {
1443 : // get index in auxiliary lists
1444 : unsigned kaux = GMM_m_type_[im] * GMM_d_size + id;
1445 : // calculate exponent of overlap
1446 : double expov = get_exp_overlap(GMM_d_m_[id], getPosition(im), inv_cov_md_[kaux]);
1447 : // get index of 0.5*expov in tabulated exponential
1448 : unsigned itab = static_cast<unsigned> (round( 0.5*expov/dexp_ ));
1449 : // check boundaries and skip atom in case
1450 : if(itab >= tab_exp_.size()) {
1451 : continue;
1452 : }
1453 : // in case calculate overlap
1454 : double ov = pre_fact_[kaux] * tab_exp_[itab];
1455 : // add to list
1456 : ov_l.push_back(std::make_pair(ov,im));
1457 : // increase ov_tot
1458 : ov_tot += ov;
1459 : }
1460 : // check if zero size -> ov_tot = 0
1461 : if(ov_l.size()==0) {
1462 : continue;
1463 : }
1464 : // define cutoff
1465 : double ov_cut = ov_tot * nl_cutoff_;
1466 : // sort ov_l in ascending order
1467 : std::sort(ov_l.begin(), ov_l.end());
1468 : // integrate ov_l and skip the smallest overlaps until reaching the cutoff
1469 : double res = 0.0;
1470 : unsigned first = 0;
1471 : for(first=0; first<ov_l.size(); ++first) {
1472 : res += ov_l[first].first;
1473 : if(res >= ov_cut) {
1474 : break;
1475 : }
1476 : }
1477 : // now add atoms to neighborlist
1478 : for(unsigned i=first; i<ov_l.size(); ++i) {
1479 : nl_l.push_back(id*GMM_m_size+ov_l[i].second);
1480 : }
1481 : }
1482 : }
1483 : // local neighbor list, in the order of the GMM components
1484 36 : std::vector<unsigned> nl_l;
1485 108 : for(unsigned t=0; t<nt; ++t) {
1486 72 : nl_l.insert(nl_l.end(), nl_t[t].begin(), nl_t[t].end());
1487 : }
1488 : // find total dimension of neighborlist
1489 36 : std::vector <int> recvcounts(size_, 0);
1490 36 : recvcounts[rank_] = nl_l.size();
1491 36 : comm.Sum(&recvcounts[0], size_);
1492 : int tot_size = accumulate(recvcounts.begin(), recvcounts.end(), 0);
1493 : // resize neighbor stuff
1494 36 : nl_.resize(tot_size);
1495 : // calculate std::vector of displacement
1496 36 : std::vector<int> disp(size_);
1497 36 : disp[0] = 0;
1498 : int rank_size = 0;
1499 68 : for(unsigned i=0; i<size_-1; ++i) {
1500 32 : rank_size += recvcounts[i];
1501 32 : disp[i+1] = rank_size;
1502 : }
1503 : // Allgather neighbor list
1504 36 : if(size_>1) {
1505 32 : comm.Allgatherv(nl_l.data(), recvcounts[rank_], nl_.data(), &recvcounts[0], &disp[0]);
1506 : } else {
1507 4 : nl_ = nl_l;
1508 : }
1509 : // now resize derivatives
1510 36 : ovmd_der_.resize(tot_size);
1511 : // contiguous block of the neighbor list computed by this rank
1512 36 : nl_begin_ = static_cast<unsigned>((static_cast<unsigned long long>(tot_size)*rank_)/size_);
1513 36 : nl_end_ = static_cast<unsigned>((static_cast<unsigned long long>(tot_size)*(rank_+1))/size_);
1514 36 : }
1515 :
1516 108 : void EMMI::prepare() {
1517 108 : if(getExchangeStep()) {
1518 0 : first_time_=true;
1519 : }
1520 108 : }
1521 :
1522 : // overlap calculator
1523 108 : void EMMI::calculate_overlap() {
1524 :
1525 108 : if(first_time_ || getExchangeStep() || getStep()%nl_stride_==0) {
1526 36 : update_neighbor_list();
1527 36 : first_time_=false;
1528 : }
1529 :
1530 : // clean temporary std::vectors
1531 : std::fill(ovmd_.begin(), ovmd_.end(), 0.0);
1532 :
1533 : // we have to cycle over all model and data GMM components in the neighbor list:
1534 : // each rank does the block [nl_begin_,nl_end_), and the derivatives ovmd_der_
1535 : // are used only by the same rank (in calculate), so they are not communicated
1536 108 : const unsigned GMM_d_size = GMM_d_m_.size();
1537 108 : const unsigned GMM_m_size = GMM_m_type_.size();
1538 108 : const unsigned nt = OpenMP::getNumThreads();
1539 152609 : auto overlap=[&](const unsigned i, std::vector<double> & ovmd) {
1540 : // get data (id) and atom (im) indexes
1541 152609 : const unsigned id = nl_[i] / GMM_m_size;
1542 152609 : const unsigned im = nl_[i] % GMM_m_size;
1543 : // get index in auxiliary lists
1544 152609 : const unsigned kaux = GMM_m_type_[im] * GMM_d_size + id;
1545 : // add overlap with im component of model GMM
1546 305218 : ovmd[id] += get_overlap(GMM_d_m_[id], getPosition(im), pre_fact_[kaux],
1547 152609 : inv_cov_md_[kaux], ovmd_der_[i]);
1548 152717 : };
1549 108 : if(nt==1) {
1550 0 : for(unsigned i=nl_begin_; i<nl_end_; ++i) {
1551 0 : overlap(i,ovmd_);
1552 : }
1553 : } else {
1554 : // one array per thread, summed in thread order (reproducible results)
1555 108 : std::vector< std::vector<double> > ovmd_t(nt, std::vector<double>(ovmd_.size(),0.0));
1556 108 : #pragma omp parallel num_threads(nt)
1557 : {
1558 : std::vector<double> & ovmd_l = ovmd_t[OpenMP::getThreadNum()];
1559 : #pragma omp for schedule(static)
1560 : for(unsigned i=nl_begin_; i<nl_end_; ++i) {
1561 : overlap(i,ovmd_l);
1562 : }
1563 : }
1564 324 : for(unsigned t=0; t<nt; ++t) {
1565 109080 : for(unsigned id=0; id<ovmd_.size(); ++id) {
1566 108864 : ovmd_[id] += ovmd_t[t][id];
1567 : }
1568 : }
1569 108 : }
1570 : // communicate stuff
1571 108 : if(size_>1) {
1572 94 : comm.Sum(&ovmd_[0], ovmd_.size());
1573 : }
1574 108 : }
1575 :
1576 200002 : double EMMI::scaleEnergy(double s) {
1577 : double ene = 0.0;
1578 180201802 : for(unsigned i=0; i<ovdd_.size(); ++i) {
1579 180001800 : ene += std::log( std::fabs ( s * ovmd_ave_[i] - ovdd_[i] ) );
1580 : }
1581 200002 : return ene;
1582 : }
1583 :
1584 4 : double EMMI::doRegression() {
1585 : // the simulated annealing is done once, by the first rank of the first replica
1586 : // (all the replicas have the same average overlaps), and the result is broadcast
1587 4 : double res[3]= {0.0,0.0,0.0};
1588 4 : const bool do_mc = (rank_==0 && (no_aver_ || nrep_==1 || replica_==0));
1589 : if(do_mc) {
1590 : // standard MC parameters
1591 : const unsigned MCsteps = 100000;
1592 : const double kbtmin = 1.0;
1593 : const double kbtmax = 10.0;
1594 : const unsigned ncold = 5000;
1595 : const unsigned nhot = 2000;
1596 : double MCacc = 0.0;
1597 :
1598 : // initial value of scale factor and energy
1599 2 : double scale = random_.RandU01() * ( scale_max_ - scale_min_ ) + scale_min_;
1600 2 : double ene = scaleEnergy(scale);
1601 : // set best energy and scale
1602 : double ebest = ene;
1603 : double scale_best = scale;
1604 :
1605 : // MC loop
1606 200002 : for(unsigned istep=0; istep<MCsteps; ++istep) {
1607 : // get temperature
1608 200000 : const double kbt = (istep%(ncold+nhot)<ncold) ? kbtmin : kbtmax;
1609 : // propose move in scale
1610 200000 : double ds = dscale_ * ( 2.0 * random_.RandU01() - 1.0 );
1611 200000 : double new_scale = scale + ds;
1612 : // check boundaries
1613 200000 : if(new_scale > scale_max_) {
1614 0 : new_scale = 2.0 * scale_max_ - new_scale;
1615 : }
1616 200000 : if(new_scale < scale_min_) {
1617 35666 : new_scale = 2.0 * scale_min_ - new_scale;
1618 : }
1619 : // new energy
1620 200000 : double new_ene = scaleEnergy(new_scale);
1621 : // accept or reject
1622 200000 : bool accept = doAccept(ene, new_ene, kbt);
1623 : // in case of acceptance
1624 200000 : if(accept) {
1625 : scale = new_scale;
1626 : ene = new_ene;
1627 56226 : MCacc += 1.0;
1628 : }
1629 : // save best
1630 200000 : if(ene<ebest) {
1631 : ebest = ene;
1632 : scale_best = scale;
1633 : }
1634 : }
1635 2 : res[0] = scale_best;
1636 2 : res[1] = ebest;
1637 2 : res[2] = MCacc / static_cast<double>(MCsteps);
1638 : }
1639 : // global communication
1640 4 : if(!no_aver_ && nrep_>1 && rank_==0) {
1641 0 : multi_sim_comm.Bcast(&res[0],3,0);
1642 : }
1643 : // local communication
1644 4 : if(size_>1) {
1645 4 : comm.Bcast(&res[0],3,0);
1646 : }
1647 : // set scale parameters
1648 4 : getPntrToComponent("accscale")->set(res[2]);
1649 4 : getPntrToComponent("enescale")->set(res[1]);
1650 : // return scale value
1651 4 : return res[0];
1652 : }
1653 :
1654 20 : double EMMI::get_annealing(long long int step) {
1655 : // default no annealing
1656 : double fact = 1.0;
1657 : // position in annealing cycle
1658 20 : unsigned nc = step%(4*nanneal_);
1659 : // useful doubles
1660 20 : double ncd = static_cast<double>(nc);
1661 20 : double nn = static_cast<double>(nanneal_);
1662 : // set fact
1663 20 : if(nc>=nanneal_ && nc<2*nanneal_) {
1664 6 : fact = (kanneal_-1.0) / nn * ( ncd - nn ) + 1.0;
1665 : }
1666 20 : if(nc>=2*nanneal_ && nc<3*nanneal_) {
1667 4 : fact = kanneal_;
1668 : }
1669 20 : if(nc>=3*nanneal_) {
1670 2 : fact = (1.0-kanneal_) / nn * ( ncd - 3.0*nn) + kanneal_;
1671 : }
1672 20 : return fact;
1673 : }
1674 :
1675 108 : void EMMI::get_weights(double &weight, double &norm, double &neff) {
1676 108 : const double dnrep = static_cast<double>(nrep_);
1677 : // calculate the weights either from BIAS
1678 108 : if(do_reweight_) {
1679 36 : std::vector<double> bias(nrep_,0);
1680 36 : if(rank_==0) {
1681 24 : bias[replica_] = getArgument(0);
1682 24 : if(nrep_>1) {
1683 24 : multi_sim_comm.Sum(&bias[0], nrep_);
1684 : }
1685 : }
1686 36 : comm.Sum(&bias[0], nrep_);
1687 :
1688 : // accumulate weights
1689 36 : if(!first_time_w_) {
1690 90 : for(unsigned i=0; i<nrep_; ++i) {
1691 60 : const double delta=bias[i]-average_weights_[i];
1692 : // FIXME: multiplying by fractional decay here causes problems with numerical derivatives,
1693 : // probably because we're making several calls to calculate(), causing accumulation
1694 : // of epsilons. Maybe we can work on a temporary copy of the action instead?
1695 60 : average_weights_[i]+=decay_w_*delta;
1696 : }
1697 : } else {
1698 6 : first_time_w_ = false;
1699 18 : for(unsigned i=0; i<nrep_; ++i) {
1700 12 : average_weights_[i] = bias[i];
1701 : }
1702 : }
1703 :
1704 : // set average back into bias and set norm to one
1705 36 : const double maxbias = *(std::max_element(average_weights_.begin(), average_weights_.end()));
1706 108 : for(unsigned i=0; i<nrep_; ++i) {
1707 72 : bias[i] = std::exp((average_weights_[i]-maxbias)/kbt_);
1708 : }
1709 : // set local weight, norm and weight variance
1710 36 : weight = bias[replica_];
1711 : double w2=0.;
1712 108 : for(unsigned i=0; i<nrep_; ++i) {
1713 72 : w2 += bias[i]*bias[i];
1714 72 : norm += bias[i];
1715 : }
1716 36 : neff = norm*norm/w2;
1717 72 : getPntrToComponent("weight")->set(weight/norm);
1718 : } else {
1719 : // or arithmetic ones
1720 72 : neff = dnrep;
1721 72 : weight = 1.0;
1722 72 : norm = dnrep;
1723 : }
1724 108 : getPntrToComponent("neff")->set(neff);
1725 108 : }
1726 :
1727 108 : void EMMI::calculate() {
1728 :
1729 : // calculate CV
1730 108 : calculate_overlap();
1731 :
1732 : // rescale factor for ensemble average
1733 108 : double weight = 0.;
1734 108 : double neff = 0.;
1735 108 : double norm = 0.;
1736 108 : get_weights(weight, norm, neff);
1737 :
1738 : // in case of ensemble averaging, calculate average overlap
1739 108 : if(!no_aver_ && nrep_>1) {
1740 : // if master node, calculate average across replicas
1741 36 : if(rank_==0) {
1742 21624 : for(unsigned i=0; i<ovmd_.size(); ++i) {
1743 21600 : ovmd_ave_[i] = weight / norm * ovmd_[i];
1744 : }
1745 24 : multi_sim_comm.Sum(&ovmd_ave_[0], ovmd_ave_.size());
1746 : } else {
1747 10812 : for(unsigned i=0; i<ovmd_ave_.size(); ++i) {
1748 10800 : ovmd_ave_[i] = 0.0;
1749 : }
1750 : }
1751 : // local communication
1752 36 : if(size_>1) {
1753 24 : comm.Sum(&ovmd_ave_[0], ovmd_ave_.size());
1754 : }
1755 : } else {
1756 22104 : for(unsigned i=0; i<ovmd_.size(); ++i) {
1757 22032 : ovmd_ave_[i] = ovmd_[i];
1758 : }
1759 : }
1760 :
1761 : // get time step
1762 108 : long long int step = getStep();
1763 :
1764 : // do regression
1765 108 : if(nregres_>0) {
1766 12 : if(step%nregres_==0 && !getExchangeStep()) {
1767 4 : scale_ = doRegression();
1768 : }
1769 : // set scale component
1770 24 : getPntrToComponent("scale")->set(scale_);
1771 : }
1772 :
1773 : // write model overlap to file
1774 108 : if(ovstride_>0 && step%ovstride_==0) {
1775 6 : write_model_overlap(step);
1776 : }
1777 :
1778 : // clear energy and virial
1779 108 : ene_ = 0.0;
1780 : virial_.zero();
1781 :
1782 : // Gaussian noise
1783 108 : if(noise_==0) {
1784 90 : calculate_Gauss();
1785 : }
1786 :
1787 : // Outliers noise
1788 108 : if(noise_==1) {
1789 15 : calculate_Outliers();
1790 : }
1791 :
1792 : // Marginal noise
1793 108 : if(noise_==2) {
1794 3 : calculate_Marginal();
1795 : }
1796 :
1797 : // get annealing rescale factor
1798 108 : if(nanneal_>0) {
1799 20 : anneal_ = get_annealing(step);
1800 40 : getPntrToComponent("anneal")->set(anneal_);
1801 : }
1802 :
1803 : // annealing rescale
1804 108 : ene_ /= anneal_;
1805 :
1806 : // The total energy is the sum of the energies of the replicas, which depend on the
1807 : // replica average of the overlaps: its derivative with respect to the average is
1808 : // the sum over the replicas of GMMid_der_, and the derivative of the average with
1809 : // respect to the overlaps of this replica is weight/norm
1810 108 : std::vector<double> GMMid_der_sum(GMMid_der_.size(), 0.0);
1811 : // in case of ensemble averaging
1812 108 : if(!no_aver_ && nrep_>1) {
1813 : // if master node, sum der_GMMid derivatives and ene
1814 36 : if(rank_==0) {
1815 24 : GMMid_der_sum = GMMid_der_;
1816 24 : multi_sim_comm.Sum(&GMMid_der_sum[0], GMMid_der_sum.size());
1817 24 : multi_sim_comm.Sum(&ene_, 1);
1818 : } else {
1819 : // set energy to zero
1820 12 : ene_ = 0.0;
1821 : }
1822 : // local communication
1823 36 : if(size_>1) {
1824 24 : comm.Sum(&GMMid_der_sum[0], GMMid_der_sum.size());
1825 24 : comm.Sum(&ene_, 1);
1826 : }
1827 : } else {
1828 72 : GMMid_der_sum = GMMid_der_;
1829 : }
1830 : // derivative of the average with respect to the overlaps of this replica
1831 108 : const double dave = weight / norm;
1832 :
1833 : // get derivatives of bias with respect to atoms, in the block of the neighbor list of this rank
1834 108 : const unsigned GMM_m_size = GMM_m_type_.size();
1835 108 : const unsigned nt = OpenMP::getNumThreads();
1836 152609 : auto derivative=[&](const unsigned i, std::vector<Vector> & atom_der, Tensor & virial) {
1837 : // get indexes of data and model component
1838 152609 : const unsigned id = nl_[i] / GMM_m_size;
1839 152609 : const unsigned im = nl_[i] % GMM_m_size;
1840 : // chain rule + replica normalization
1841 152609 : const Vector tot_der = dave * GMMid_der_sum[id] * ovmd_der_[i] * scale_ / anneal_;
1842 152609 : Vector pos;
1843 152609 : if(pbc_) {
1844 0 : pos = pbcDistance(GMM_d_m_[id], getPosition(im)) + GMM_d_m_[id];
1845 : } else {
1846 152609 : pos = getPosition(im);
1847 : }
1848 : // increment derivatives and virial
1849 152609 : atom_der[im] += tot_der;
1850 152609 : virial += Tensor(pos, -tot_der);
1851 108 : };
1852 : std::fill(atom_der_.begin(), atom_der_.end(), Vector(0,0,0));
1853 108 : if(nt==1) {
1854 0 : for(unsigned i=nl_begin_; i<nl_end_; ++i) {
1855 0 : derivative(i, atom_der_, virial_);
1856 : }
1857 : } else {
1858 : // one array per thread, summed in thread order (reproducible results)
1859 108 : std::vector< std::vector<Vector> > atom_der_t(nt, std::vector<Vector>(atom_der_.size()));
1860 108 : std::vector<Tensor> virial_t(nt);
1861 108 : #pragma omp parallel num_threads(nt)
1862 : {
1863 : const unsigned t=OpenMP::getThreadNum();
1864 : #pragma omp for schedule(static)
1865 : for(unsigned i=nl_begin_; i<nl_end_; ++i) {
1866 : derivative(i, atom_der_t[t], virial_t[t]);
1867 : }
1868 : }
1869 324 : for(unsigned t=0; t<nt; ++t) {
1870 59208 : for(unsigned im=0; im<atom_der_.size(); ++im) {
1871 58992 : atom_der_[im] += atom_der_t[t][im];
1872 : }
1873 216 : virial_ += virial_t[t];
1874 : }
1875 108 : }
1876 :
1877 : // communicate local derivatives and virial
1878 108 : if(size_>1) {
1879 94 : comm.Sum(&atom_der_[0][0], 3*atom_der_.size());
1880 94 : comm.Sum(virial_);
1881 : }
1882 :
1883 : // set derivatives, virial, and score
1884 29604 : for(unsigned i=0; i<atom_der_.size(); ++i) {
1885 29496 : setAtomsDerivatives(valueScore_, i, atom_der_[i]);
1886 : }
1887 108 : setBoxDerivatives(valueScore_, virial_);
1888 108 : valueScore_->set(ene_);
1889 :
1890 108 : if (do_reweight_) {
1891 : // derivative of the total energy with respect to the bias of this replica
1892 : double w_tmp = 0.;
1893 32436 : for (unsigned i = 0; i < ovmd_.size(); ++i) {
1894 32400 : w_tmp += (ovmd_[i] - ovmd_ave_[i]) * GMMid_der_sum[i];
1895 : }
1896 36 : w_tmp *= scale_ * dave / kbt_ * decay_w_ / anneal_;
1897 :
1898 36 : setArgDerivatives(valueScore_, w_tmp);
1899 72 : getPntrToComponent("biasDer")->set(w_tmp);
1900 : }
1901 :
1902 : // This part is needed only for Gaussian and Outliers noise models
1903 108 : if(noise_!=2) {
1904 :
1905 : // do Montecarlo
1906 105 : if(dsigma_[0]>0 && step%MCstride_==0 && !getExchangeStep()) {
1907 24 : doMonteCarlo();
1908 : }
1909 :
1910 : // print status
1911 105 : if(step%statusstride_==0) {
1912 25 : print_status(step);
1913 : }
1914 :
1915 : // calculate acceptance ratio
1916 105 : const double acc = (MCtrials_>0.0) ? MCaccept_ / MCtrials_ : 0.0;
1917 :
1918 : // set value
1919 105 : valueAcc_->set(acc);
1920 :
1921 : }
1922 :
1923 108 : }
1924 :
1925 90 : void EMMI::calculate_Gauss() {
1926 : // cycle on all the GMM groups
1927 216 : for(unsigned i=0; i<GMM_d_grps_.size(); ++i) {
1928 : double eneg = 0.0;
1929 : // cycle on all the members of the group
1930 43704 : for(unsigned j=0; j<GMM_d_grps_[i].size(); ++j) {
1931 : // id of the GMM component
1932 43578 : int GMMid = GMM_d_grps_[i][j];
1933 : // calculate deviation
1934 43578 : double dev = ( scale_*ovmd_ave_[GMMid]-ovdd_[GMMid] ) / sigma_[i];
1935 : // add to group energy
1936 43578 : eneg += 0.5 * dev * dev;
1937 : // store derivative for later
1938 43578 : GMMid_der_[GMMid] = kbt_ * dev / sigma_[i];
1939 : }
1940 : // add to total energy along with normalizations and prior
1941 126 : ene_ += kbt_ * ( eneg + (static_cast<double>(GMM_d_grps_[i].size())+prior_) * std::log(sigma_[i]) );
1942 : }
1943 90 : }
1944 :
1945 15 : void EMMI::calculate_Outliers() {
1946 : // cycle on all the GMM groups
1947 42 : for(unsigned i=0; i<GMM_d_grps_.size(); ++i) {
1948 : // cycle on all the members of the group
1949 : double eneg = 0.0;
1950 10854 : for(unsigned j=0; j<GMM_d_grps_[i].size(); ++j) {
1951 : // id of the GMM component
1952 10827 : int GMMid = GMM_d_grps_[i][j];
1953 : // calculate deviation
1954 10827 : double dev = ( scale_*ovmd_ave_[GMMid]-ovdd_[GMMid] ) / sigma_[i];
1955 : // add to group energy
1956 10827 : eneg += std::log1p( 0.5 * dev * dev );
1957 : // store derivative for later
1958 10827 : GMMid_der_[GMMid] = kbt_ / ( 1.0 + 0.5 * dev * dev ) * dev / sigma_[i];
1959 : }
1960 : // add to total energy along with normalizations and prior
1961 27 : ene_ += kbt_ * ( eneg + (static_cast<double>(GMM_d_grps_[i].size())+prior_) * std::log(sigma_[i]) );
1962 : }
1963 15 : }
1964 :
1965 3 : void EMMI::calculate_Marginal() {
1966 : // cycle on all the GMM groups
1967 6 : for(unsigned i=0; i<GMM_d_grps_.size(); ++i) {
1968 : // cycle on all the members of the group
1969 30 : for(unsigned j=0; j<GMM_d_grps_[i].size(); ++j) {
1970 : // id of the GMM component
1971 27 : int GMMid = GMM_d_grps_[i][j];
1972 : // calculate deviation
1973 27 : double dev = ( scale_*ovmd_ave_[GMMid]-ovdd_[GMMid] );
1974 27 : if(dev==0.0) {
1975 : // limit for dev -> 0: erf(dev/sqrt(2)/smin)/dev -> sqrt(2/pi)/smin, and zero derivative
1976 0 : ene_ += -kbt_ * std::log ( 0.5 * sqrt2_pi_ / sigma_min_[i] );
1977 0 : GMMid_der_[GMMid] = 0.0;
1978 0 : continue;
1979 : }
1980 : // calculate errf
1981 27 : double errf = erf ( dev * inv_sqrt2_ / sigma_min_[i] );
1982 : // add to group energy
1983 27 : ene_ += -kbt_ * std::log ( 0.5 / dev * errf ) ;
1984 : // store derivative for later
1985 27 : GMMid_der_[GMMid] = - kbt_/errf*sqrt2_pi_*std::exp(-0.5*dev*dev/sigma_min_[i]/sigma_min_[i])/sigma_min_[i]+kbt_/dev;
1986 : }
1987 : }
1988 3 : }
1989 :
1990 : }
1991 : }
|
