1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2    Copyright (c) 2017-2021 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  This class was originally written by Alexander Jussupow
24  Extension for the middleman algorithm by Max Muehlbauer
25  Martini beads strucutre factor for Nucleic Acids by Cristina Paissoni
26 */
27 
28 #include "MetainferenceBase.h"
29 #include "core/ActionRegister.h"
30 #include "core/ActionSet.h"
31 #include "core/GenericMolInfo.h"
32 #include "tools/Communicator.h"
33 #include "tools/Pbc.h"
34 
35 #include <string>
36 #include <cmath>
37 #include <map>
38 
39 #ifdef __PLUMED_HAS_ARRAYFIRE
40 #include <arrayfire.h>
41 #include <af/util.h>
42 #endif
43 
44 #ifndef M_PI
45 #define M_PI           3.14159265358979323846
46 #endif
47 
48 using namespace std;
49 
50 namespace PLMD {
51 namespace isdb {
52 
53 //+PLUMEDOC ISDB_COLVAR SAXS
54 /*
55 Calculates SAXS scattered intensity using either the Debye equation.
56 
57 Intensities are calculated for a set of scattering length set using QVALUE keywords that are numbered starting from 0.
58 Structure factors can be either assigned using a polynomial expansion to any order using the PARAMETERS keywords;
59 automatically assigned to atoms using the ATOMISTIC flag reading a PDB file, a correction for the water density is
60 automatically added, with water density that by default is 0.334 but that can be set otherwise using WATERDENS;
61 automatically assigned to Martini pseudo atoms using the MARTINI flag.
62 The calculated intensities can be scaled using the SCALEINT keywords. This is applied by rescaling the structure factors.
63 Experimental reference intensities can be added using the EXPINT keywords.
64 By default SAXS is calculated using Debye on CPU, by adding the GPU flag it is possible to solve the equation on a GPU
65 if the ARRAYFIRE libraries are installed and correctly linked.
66 \ref METAINFERENCE can be activated using DOSCORE and the other relevant keywords.
67 
68 \par Examples
69 in the following example the saxs intensities for a martini model are calculated. structure factors
70 are obtained from the pdb file indicated in the MOLINFO.
71 
72 \plumedfile
73 MOLINFO STRUCTURE=template.pdb
74 
75 SAXS ...
76 LABEL=saxs
77 ATOMS=1-355
78 SCALEINT=3920000
79 MARTINI
80 QVALUE1=0.02 EXPINT1=1.0902
81 QVALUE2=0.05 EXPINT2=0.790632
82 QVALUE3=0.08 EXPINT3=0.453808
83 QVALUE4=0.11 EXPINT4=0.254737
84 QVALUE5=0.14 EXPINT5=0.154928
85 QVALUE6=0.17 EXPINT6=0.0921503
86 QVALUE7=0.2 EXPINT7=0.052633
87 QVALUE8=0.23 EXPINT8=0.0276557
88 QVALUE9=0.26 EXPINT9=0.0122775
89 QVALUE10=0.29 EXPINT10=0.00880634
90 QVALUE11=0.32 EXPINT11=0.0137301
91 QVALUE12=0.35 EXPINT12=0.0180036
92 QVALUE13=0.38 EXPINT13=0.0193374
93 QVALUE14=0.41 EXPINT14=0.0210131
94 QVALUE15=0.44 EXPINT15=0.0220506
95 ... SAXS
96 
97 PRINT ARG=(saxs\.q-.*),(saxs\.exp-.*) FILE=colvar STRIDE=1
98 
99 \endplumedfile
100 
101 */
102 //+ENDPLUMEDOC
103 
104 class SAXS :
105   public MetainferenceBase
106 {
107 private:
108   enum { H, C, N, O, P, S, NTT };
109   enum { ALA_BB, ARG_BB, ARG_SC1, ARG_SC2, ASN_BB, ASN_SC1, ASP_BB, ASP_SC1, CYS_BB, CYS_SC1,
110          GLN_BB, GLN_SC1, GLU_BB, GLU_SC1, GLY_BB, HIS_BB, HIS_SC1, HIS_SC2, HIS_SC3, ILE_BB,
111          ILE_SC1, LEU_BB, LEU_SC1, LYS_BB, LYS_SC1, LYS_SC2, MET_BB, MET_SC1, PHE_BB, PHE_SC1,
112          PHE_SC2, PHE_SC3, PRO_BB, PRO_SC1, SER_BB, SER_SC1, THR_BB, THR_SC1, TRP_BB, TRP_SC1,
113          TRP_SC2, TRP_SC3, TRP_SC4, TYR_BB, TYR_SC1, TYR_SC2, TYR_SC3, VAL_BB, VAL_SC1, A_BB1,
114          A_BB2, A_BB3, A_SC1, A_SC2, A_SC3, A_SC4, A_3TE, A_5TE, A_TE3, A_TE5, C_BB1, C_BB2,
115          C_BB3, C_SC1, C_SC2, C_SC3, C_3TE, C_5TE, C_TE3, C_TE5, G_BB1, G_BB2, G_BB3, G_SC1,
116          G_SC2, G_SC3, G_SC4, G_3TE, G_5TE, G_TE3, G_TE5, U_BB1, U_BB2, U_BB3, U_SC1, U_SC2,
117          U_SC3, U_3TE, U_5TE, U_TE3, U_TE5, DA_BB1, DA_BB2, DA_BB3, DA_SC1, DA_SC2, DA_SC3,
118          DA_SC4, DA_3TE, DA_5TE, DA_TE3, DA_TE5, DC_BB1, DC_BB2, DC_BB3, DC_SC1, DC_SC2, DC_SC3,
119          DC_3TE, DC_5TE, DC_TE3, DC_TE5, DG_BB1, DG_BB2, DG_BB3, DG_SC1, DG_SC2, DG_SC3, DG_SC4,
120          DG_3TE, DG_5TE, DG_TE3, DG_TE5, DT_BB1, DT_BB2, DT_BB3, DT_SC1, DT_SC2, DT_SC3, DT_3TE,
121          DT_5TE, DT_TE3, DT_TE5, NMARTINI
122        };
123   bool                       pbc;
124   bool                       serial;
125   bool                       gpu;
126   int                        deviceid;
127   vector<unsigned>           atoi;
128   vector<double>             q_list;
129   vector<double>             FF_rank;
130   vector<vector<double> >    FF_value;
131   vector<vector<float> >     FFf_value;
132 
133   void calculate_gpu(vector<Vector> &deriv);
134   void calculate_cpu(vector<Vector> &deriv);
135   void getMartiniSFparam(const vector<AtomNumber> &atoms, vector<vector<long double> > &parameter);
136   double calculateASF(const vector<AtomNumber> &atoms, vector<vector<long double> > &FF_tmp, const double rho);
137 
138 public:
139   static void registerKeywords( Keywords& keys );
140   explicit SAXS(const ActionOptions&);
141   void calculate() override;
142   void update() override;
143 };
144 
145 PLUMED_REGISTER_ACTION(SAXS,"SAXS")
146 
registerKeywords(Keywords & keys)147 void SAXS::registerKeywords(Keywords& keys) {
148   componentsAreNotOptional(keys);
149   MetainferenceBase::registerKeywords(keys);
150   keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
151   keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
152   keys.add("compulsory","DEVICEID","0","Identifier of the GPU to be used");
153   keys.addFlag("GPU",false,"calculate SAXS using ARRAYFIRE on an accelerator device");
154   keys.addFlag("ATOMISTIC",false,"calculate SAXS for an atomistic model");
155   keys.addFlag("MARTINI",false,"calculate SAXS for a Martini model");
156   keys.add("atoms","ATOMS","The atoms to be included in the calculation, e.g. the whole protein.");
157   keys.add("numbered","QVALUE","Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ... .");
158   keys.add("numbered","PARAMETERS","Used parameter Keywords like PARAMETERS1, PARAMETERS2. These are used to calculate the structure factor for the \\f$i\\f$th atom/bead.");
159   keys.add("compulsory","WATERDENS","0.334","Density of the water to be used for the correction of atomistic structure factors.");
160   keys.add("numbered","EXPINT","Add an experimental value for each q value.");
161   keys.add("compulsory","SCALEINT","1.0","SCALING value of the calculated data. Useful to simplify the comparison.");
162   keys.addOutputComponent("q","default","the # SAXS of q");
163   keys.addOutputComponent("exp","EXPINT","the # experimental intensity");
164 }
165 
SAXS(const ActionOptions & ao)166 SAXS::SAXS(const ActionOptions&ao):
167   PLUMED_METAINF_INIT(ao),
168   pbc(true),
169   serial(false),
170   gpu(false),
171   deviceid(0)
172 {
173   vector<AtomNumber> atoms;
174   parseAtomList("ATOMS",atoms);
175   const unsigned size = atoms.size();
176 
177   parseFlag("SERIAL",serial);
178 
179   bool nopbc=!pbc;
180   parseFlag("NOPBC",nopbc);
181   pbc=!nopbc;
182   if(pbc)      log.printf("  using periodic boundary conditions\n");
183   else         log.printf("  without periodic boundary conditions\n");
184 
185   parseFlag("GPU",gpu);
186 #ifndef  __PLUMED_HAS_ARRAYFIRE
187   if(gpu) error("To use the GPU mode PLUMED must be compiled with ARRAYFIRE");
188 #endif
189 
190   parse("DEVICEID",deviceid);
191 #ifdef  __PLUMED_HAS_ARRAYFIRE
192   if(gpu) {
193     af::setDevice(deviceid);
194     af::info();
195   }
196 #endif
197 
198   unsigned ntarget=0;
199   for(unsigned i=0;; ++i) {
200     double t_list;
201     if( !parseNumbered( "QVALUE", i+1, t_list) ) break;
202     if(t_list<=0.) error("QVALUE cannot be less or equal to zero!\n");
203     q_list.push_back(t_list);
204     ntarget++;
205   }
206   const unsigned numq = ntarget;
207 
208   for(unsigned i=0; i<numq; i++) {
209     if(q_list[i]==0.) error("it is not possible to set q=0\n");
210     if(i>0&&q_list[i]<q_list[i-1]) error("QVALUE must be in ascending order");
211     log.printf("  my q: %lf \n",q_list[i]);
212   }
213 
214   bool atomistic=false;
215   parseFlag("ATOMISTIC",atomistic);
216   bool martini=false;
217   parseFlag("MARTINI",martini);
218 
219   if(martini&&atomistic) error("You cannot use MARTINI and ATOMISTIC at the same time");
220 
221   double rho = 0.334;
222   parse("WATERDENS", rho);
223 
224   double Iq0=0;
225   vector<vector<long double> > FF_tmp;
226   atoi.resize(size);
227   if(!atomistic&&!martini) {
228     //read in parameter vector
229     vector<vector<long double> > parameter;
230     parameter.resize(size);
231     ntarget=0;
232     for(unsigned i=0; i<size; ++i) {
233       if( !parseNumberedVector( "PARAMETERS", i+1, parameter[i]) ) break;
234       ntarget++;
235     }
236     if( ntarget!=size ) error("found wrong number of parameter vectors");
237     FF_tmp.resize(numq,vector<long double>(size));
238     for(unsigned i=0; i<size; ++i) {
239       atoi[i]=i;
240       for(unsigned k=0; k<numq; ++k) {
241         for(unsigned j=0; j<parameter[i].size(); ++j) {
242           FF_tmp[k][i]+= parameter[i][j]*powl(static_cast<long double>(q_list[k]),j);
243         }
244       }
245     }
246     for(unsigned i=0; i<size; ++i) Iq0+=parameter[i][0];
247   } else if(martini) {
248     //read in parameter vector
249     FF_tmp.resize(numq,vector<long double>(NMARTINI));
250     vector<vector<long double> > parameter;
251     parameter.resize(NMARTINI);
252     getMartiniSFparam(atoms, parameter);
253     for(unsigned i=0; i<NMARTINI; ++i) {
254       for(unsigned k=0; k<numq; ++k) {
255         for(unsigned j=0; j<parameter[i].size(); ++j) {
256           FF_tmp[k][i]+= parameter[i][j]*powl(static_cast<long double>(q_list[k]),j);
257         }
258       }
259     }
260     for(unsigned i=0; i<size; ++i) Iq0+=parameter[atoi[i]][0];
261   } else if(atomistic) {
262     FF_tmp.resize(numq,vector<long double>(NTT));
263     Iq0=calculateASF(atoms, FF_tmp, rho);
264   }
265   double scale_int = Iq0*Iq0;
266 
267   vector<double> expint;
268   expint.resize( numq );
269   ntarget=0;
270   for(unsigned i=0; i<numq; ++i) {
271     if( !parseNumbered( "EXPINT", i+1, expint[i] ) ) break;
272     ntarget++;
273   }
274   bool exp=false;
275   if(ntarget!=numq && ntarget!=0) error("found wrong number of EXPINT values");
276   if(ntarget==numq) exp=true;
277   if(getDoScore()&&!exp) error("with DOSCORE you need to set the EXPINT values");
278 
279   double tmp_scale_int=1.;
280   parse("SCALEINT",tmp_scale_int);
281 
282   if(tmp_scale_int!=1) scale_int /= tmp_scale_int;
283   else {
284     if(exp) scale_int /= expint[0];
285   }
286 
287   if(!gpu) {
288     FF_rank.resize(numq);
289     unsigned n_atom_types=size;
290     if(atomistic) n_atom_types=NTT;
291     else if(martini) n_atom_types=NMARTINI;
292     FF_value.resize(n_atom_types,vector<double>(numq));
293     for(unsigned k=0; k<numq; ++k) {
294       for(unsigned i=0; i<n_atom_types; i++) FF_value[i][k] = static_cast<double>(FF_tmp[k][i])/sqrt(scale_int);
295       for(unsigned i=0; i<size; i++) FF_rank[k] += FF_value[atoi[i]][k]*FF_value[atoi[i]][k];
296     }
297   } else {
298     FFf_value.resize(numq,vector<float>(size));
299     for(unsigned k=0; k<numq; ++k) {
300       for(unsigned i=0; i<size; i++) {
301         FFf_value[k][i] = static_cast<float>(FF_tmp[k][atoi[i]])/sqrt(scale_int);
302       }
303     }
304   }
305 
306   if(!getDoScore()) {
307     for(unsigned i=0; i<numq; i++) {
308       std::string num; Tools::convert(i,num);
309       addComponentWithDerivatives("q-"+num);
310       componentIsNotPeriodic("q-"+num);
311     }
312     if(exp) {
313       for(unsigned i=0; i<numq; i++) {
314         std::string num; Tools::convert(i,num);
315         addComponent("exp-"+num);
316         componentIsNotPeriodic("exp-"+num);
317         Value* comp=getPntrToComponent("exp-"+num);
318         comp->set(expint[i]);
319       }
320     }
321   } else {
322     for(unsigned i=0; i<numq; i++) {
323       std::string num; Tools::convert(i,num);
324       addComponent("q-"+num);
325       componentIsNotPeriodic("q-"+num);
326     }
327     for(unsigned i=0; i<numq; i++) {
328       std::string num; Tools::convert(i,num);
329       addComponent("exp-"+num);
330       componentIsNotPeriodic("exp-"+num);
331       Value* comp=getPntrToComponent("exp-"+num);
332       comp->set(expint[i]);
333     }
334   }
335 
336   // convert units to nm^-1
337   for(unsigned i=0; i<numq; ++i) {
338     q_list[i]=q_list[i]*10.0;    //factor 10 to convert from A^-1 to nm^-1
339   }
340   log<<"  Bibliography ";
341   if(martini) {
342     log<<plumed.cite("Niebling, Björling, Westenhoff, J Appl Crystallogr 47, 1190–1198 (2014).");
343     log<<plumed.cite("Paissoni, Jussupow, Camilloni, J Appl Crystallogr 52, 394-402 (2019).");
344   }
345   if(atomistic) {
346     log<<plumed.cite("Fraser, MacRae, Suzuki, J. Appl. Crystallogr., 11, 693–694 (1978).");
347     log<<plumed.cite("Brown, Fox, Maslen, O'Keefe, Willis, International Tables for Crystallography C, 554–595 (International Union of Crystallography, 2006).");
348   }
349   log<< plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)");
350   log<<"\n";
351 
352   requestAtoms(atoms, false);
353   if(getDoScore()) {
354     setParameters(expint);
355     Initialise(numq);
356   }
357   setDerivatives();
358   checkRead();
359 }
360 
calculate_gpu(vector<Vector> & deriv)361 void SAXS::calculate_gpu(vector<Vector> &deriv)
362 {
363 #ifdef __PLUMED_HAS_ARRAYFIRE
364   const unsigned size = getNumberOfAtoms();
365   const unsigned numq = q_list.size();
366 
367   std::vector<float> sum;
368   sum.resize(numq);
369 
370   std::vector<float> dd;
371   dd.resize(size*3*numq);
372 
373   // on gpu only the master rank run the calculation
374   if(comm.Get_rank()==0) {
375     vector<float> posi;
376     posi.resize(3*size);
377     #pragma omp parallel for num_threads(OpenMP::getNumThreads())
378     for (unsigned i=0; i<size; i++) {
379       const Vector tmp = getPosition(i);
380       posi[3*i]   = static_cast<float>(tmp[0]);
381       posi[3*i+1] = static_cast<float>(tmp[1]);
382       posi[3*i+2] = static_cast<float>(tmp[2]);
383     }
384 
385     // create array a and b containing atomic coordinates
386     af::setDevice(deviceid);
387     // 3,size,1,1
388     af::array pos_a = af::array(3, size, &posi.front());
389     // size,3,1,1
390     pos_a = af::moddims(pos_a.T(), size, 1, 3);
391     // size,3,1,1
392     af::array pos_b = pos_a(af::span, af::span);
393     // size,1,3,1
394     pos_a = af::moddims(pos_a, size, 1, 3);
395     // 1,size,3,1
396     pos_b = af::moddims(pos_b, 1, size, 3);
397 
398     // size,size,3,1
399     af::array xyz_dist = af::tile(pos_a, 1, size, 1) - af::tile(pos_b, size, 1, 1);
400     // size,size,1,1
401     af::array square = af::sum(xyz_dist*xyz_dist,2);
402     // size,size,1,1
403     af::array dist_sqrt = af::sqrt(square);
404     // replace the zero of square with one to avoid nan in the derivatives (the number does not matter because this are multiplied by zero)
405     af::replace(square,!(af::iszero(square)),1.);
406     // size,size,3,1
407     xyz_dist = xyz_dist / af::tile(square, 1, 1, 3);
408     // numq,1,1,1
409     af::array sum_device   = af::constant(0, numq, f32);
410     // numq,size,3,1
411     af::array deriv_device = af::constant(0, numq, size, 3, f32);
412 
413     for (unsigned k=0; k<numq; k++) {
414       // calculate FF matrix
415       // size,1,1,1
416       af::array AFF_value(size, &FFf_value[k].front());
417       // size,size,1,1
418       af::array FFdist_mod = af::tile(AFF_value(af::span), 1, size)*af::transpose(af::tile(AFF_value(af::span), 1, size));
419 
420       // get q
421       const float qvalue = static_cast<float>(q_list[k]);
422       // size,size,1,1
423       af::array dist_q = qvalue*dist_sqrt;
424       // size,size,1
425       af::array dist_sin = af::sin(dist_q)/dist_q;
426       af::replace(dist_sin,!(af::isNaN(dist_sin)),1.);
427       // 1,1,1,1
428       sum_device(k) = af::sum(af::flat(dist_sin)*af::flat(FFdist_mod));
429 
430       // size,size,1,1
431       af::array tmp = FFdist_mod*(dist_sin - af::cos(dist_q));
432       // size,size,3,1
433       af::array dd_all = af::tile(tmp, 1, 1, 3)*xyz_dist;
434       // it should become 1,size,3
435       deriv_device(k, af::span, af::span) = af::sum(dd_all,0);
436     }
437 
438     // read out results
439     sum_device.host(&sum.front());
440 
441     deriv_device = af::reorder(deriv_device, 2, 1, 0);
442     deriv_device = af::flat(deriv_device);
443     deriv_device.host(&dd.front());
444   }
445 
446   comm.Bcast(dd, 0);
447   comm.Bcast(sum, 0);
448 
449   for(unsigned k=0; k<numq; k++) {
450     string num; Tools::convert(k,num);
451     Value* val=getPntrToComponent("q-"+num);
452     val->set(sum[k]);
453     if(getDoScore()) setCalcData(k, sum[k]);
454     for(unsigned i=0; i<size; i++) {
455       const unsigned di = k*size*3+i*3;
456       deriv[k*size+i] = Vector(2.*dd[di+0],2.*dd[di+1],2.*dd[di+2]);
457     }
458   }
459 #endif
460 }
461 
calculate_cpu(vector<Vector> & deriv)462 void SAXS::calculate_cpu(vector<Vector> &deriv)
463 {
464   const unsigned size = getNumberOfAtoms();
465   const unsigned numq = q_list.size();
466 
467   unsigned stride = comm.Get_size();
468   unsigned rank   = comm.Get_rank();
469   if(serial) {
470     stride = 1;
471     rank   = 0;
472   }
473 
474   vector<double> sum(numq,0);
475   unsigned nt=OpenMP::getNumThreads();
476   #pragma omp parallel num_threads(nt)
477   {
478     std::vector<Vector> omp_deriv(deriv.size());
479     std::vector<double> omp_sum(numq,0);
480     #pragma omp for nowait
481     for (unsigned i=rank; i<size-1; i+=stride) {
482       Vector posi = getPosition(i);
483       for (unsigned j=i+1; j<size ; j++) {
484         Vector c_distances = delta(posi,getPosition(j));
485         double m_distances = c_distances.modulo();
486         c_distances = c_distances/m_distances/m_distances;
487         for (unsigned k=0; k<numq; k++) {
488           unsigned kdx=k*size;
489           double qdist = q_list[k]*m_distances;
490           double FFF = 2.*FF_value[atoi[i]][k]*FF_value[atoi[j]][k];
491           double tsq = sin(qdist)/qdist;
492           double tcq = cos(qdist);
493           double tmp = FFF*(tcq-tsq);
494           Vector dd  = c_distances*tmp;
495           if(nt>1) {
496             omp_deriv[kdx+i] -=dd;
497             omp_deriv[kdx+j] +=dd;
498             omp_sum[k] += FFF*tsq;
499           } else {
500             deriv[kdx+i] -= dd;
501             deriv[kdx+j] += dd;
502             sum[k]       += FFF*tsq;
503           }
504         }
505       }
506     }
507     #pragma omp critical
508     if(nt>1) {
509       for(unsigned i=0; i<deriv.size(); i++) deriv[i]+=omp_deriv[i];
510       for(unsigned k=0; k<numq; k++) sum[k]+=omp_sum[k];
511     }
512   }
513 
514   if(!serial) {
515     comm.Sum(&deriv[0][0], 3*deriv.size());
516     comm.Sum(&sum[0], numq);
517   }
518 
519   for (unsigned k=0; k<numq; k++) {
520     sum[k]+=FF_rank[k];
521     string num; Tools::convert(k,num);
522     Value* val=getPntrToComponent("q-"+num);
523     val->set(sum[k]);
524     if(getDoScore()) setCalcData(k, sum[k]);
525   }
526 }
527 
calculate()528 void SAXS::calculate()
529 {
530   if(pbc) makeWhole();
531 
532   const unsigned size = getNumberOfAtoms();
533   const unsigned numq = q_list.size();
534 
535   vector<Vector> deriv(numq*size);
536   if(gpu) calculate_gpu(deriv);
537   else calculate_cpu(deriv);
538 
539   if(getDoScore()) {
540     /* Metainference */
541     double score = getScore();
542     setScore(score);
543   }
544 
545   for (unsigned k=0; k<numq; k++) {
546     const unsigned kdx=k*size;
547     Tensor deriv_box;
548     Value* val;
549     if(!getDoScore()) {
550       string num; Tools::convert(k,num);
551       val=getPntrToComponent("q-"+num);
552       for(unsigned i=0; i<size; i++) {
553         setAtomsDerivatives(val, i, deriv[kdx+i]);
554         deriv_box += Tensor(getPosition(i),deriv[kdx+i]);
555       }
556     } else {
557       val=getPntrToComponent("score");
558       for(unsigned i=0; i<size; i++) {
559         setAtomsDerivatives(val, i, deriv[kdx+i]*getMetaDer(k));
560         deriv_box += Tensor(getPosition(i),deriv[kdx+i]*getMetaDer(k));
561       }
562     }
563     setBoxDerivatives(val, -deriv_box);
564   }
565 }
566 
update()567 void SAXS::update() {
568   // write status file
569   if(getWstride()>0&& (getStep()%getWstride()==0 || getCPT()) ) writeStatus();
570 }
571 
getMartiniSFparam(const vector<AtomNumber> & atoms,vector<vector<long double>> & parameter)572 void SAXS::getMartiniSFparam(const vector<AtomNumber> &atoms, vector<vector<long double> > &parameter)
573 {
574   parameter[ALA_BB].push_back(9.045);
575   parameter[ALA_BB].push_back(-0.098114);
576   parameter[ALA_BB].push_back(7.54281);
577   parameter[ALA_BB].push_back(-1.97438);
578   parameter[ALA_BB].push_back(-8.32689);
579   parameter[ALA_BB].push_back(6.09318);
580   parameter[ALA_BB].push_back(-1.18913);
581 
582   parameter[ARG_BB].push_back(10.729);
583   parameter[ARG_BB].push_back(-0.0392574);
584   parameter[ARG_BB].push_back(1.15382);
585   parameter[ARG_BB].push_back(-0.155999);
586   parameter[ARG_BB].push_back(-2.43619);
587   parameter[ARG_BB].push_back(1.72922);
588   parameter[ARG_BB].push_back(-0.33799);
589 
590   parameter[ARG_SC1].push_back(-2.796);
591   parameter[ARG_SC1].push_back(0.472403);
592   parameter[ARG_SC1].push_back(8.07424);
593   parameter[ARG_SC1].push_back(4.37299);
594   parameter[ARG_SC1].push_back(-10.7398);
595   parameter[ARG_SC1].push_back(4.95677);
596   parameter[ARG_SC1].push_back(-0.725797);
597 
598   parameter[ARG_SC2].push_back(15.396);
599   parameter[ARG_SC2].push_back(0.0636736);
600   parameter[ARG_SC2].push_back(-1.258);
601   parameter[ARG_SC2].push_back(1.93135);
602   parameter[ARG_SC2].push_back(-4.45031);
603   parameter[ARG_SC2].push_back(2.49356);
604   parameter[ARG_SC2].push_back(-0.410721);
605 
606   parameter[ASN_BB].push_back(10.738);
607   parameter[ASN_BB].push_back(-0.0402162);
608   parameter[ASN_BB].push_back(1.03007);
609   parameter[ASN_BB].push_back(-0.254174);
610   parameter[ASN_BB].push_back(-2.12015);
611   parameter[ASN_BB].push_back(1.55535);
612   parameter[ASN_BB].push_back(-0.30963);
613 
614   parameter[ASN_SC1].push_back(9.249);
615   parameter[ASN_SC1].push_back(-0.0148678);
616   parameter[ASN_SC1].push_back(5.52169);
617   parameter[ASN_SC1].push_back(0.00853212);
618   parameter[ASN_SC1].push_back(-6.71992);
619   parameter[ASN_SC1].push_back(3.93622);
620   parameter[ASN_SC1].push_back(-0.64973);
621 
622   parameter[ASP_BB].push_back(10.695);
623   parameter[ASP_BB].push_back(-0.0410247);
624   parameter[ASP_BB].push_back(1.03656);
625   parameter[ASP_BB].push_back(-0.298558);
626   parameter[ASP_BB].push_back(-2.06064);
627   parameter[ASP_BB].push_back(1.53495);
628   parameter[ASP_BB].push_back(-0.308365);
629 
630   parameter[ASP_SC1].push_back(9.476);
631   parameter[ASP_SC1].push_back(-0.0254664);
632   parameter[ASP_SC1].push_back(5.57899);
633   parameter[ASP_SC1].push_back(-0.395027);
634   parameter[ASP_SC1].push_back(-5.9407);
635   parameter[ASP_SC1].push_back(3.48836);
636   parameter[ASP_SC1].push_back(-0.569402);
637 
638   parameter[CYS_BB].push_back(10.698);
639   parameter[CYS_BB].push_back(-0.0233493);
640   parameter[CYS_BB].push_back(1.18257);
641   parameter[CYS_BB].push_back(0.0684464);
642   parameter[CYS_BB].push_back(-2.792);
643   parameter[CYS_BB].push_back(1.88995);
644   parameter[CYS_BB].push_back(-0.360229);
645 
646   parameter[CYS_SC1].push_back(8.199);
647   parameter[CYS_SC1].push_back(-0.0261569);
648   parameter[CYS_SC1].push_back(6.79677);
649   parameter[CYS_SC1].push_back(-0.343845);
650   parameter[CYS_SC1].push_back(-5.03578);
651   parameter[CYS_SC1].push_back(2.7076);
652   parameter[CYS_SC1].push_back(-0.420714);
653 
654   parameter[GLN_BB].push_back(10.728);
655   parameter[GLN_BB].push_back(-0.0391984);
656   parameter[GLN_BB].push_back(1.09264);
657   parameter[GLN_BB].push_back(-0.261555);
658   parameter[GLN_BB].push_back(-2.21245);
659   parameter[GLN_BB].push_back(1.62071);
660   parameter[GLN_BB].push_back(-0.322325);
661 
662   parameter[GLN_SC1].push_back(8.317);
663   parameter[GLN_SC1].push_back(-0.229045);
664   parameter[GLN_SC1].push_back(12.6338);
665   parameter[GLN_SC1].push_back(-7.6719);
666   parameter[GLN_SC1].push_back(-5.8376);
667   parameter[GLN_SC1].push_back(5.53784);
668   parameter[GLN_SC1].push_back(-1.12604);
669 
670   parameter[GLU_BB].push_back(10.694);
671   parameter[GLU_BB].push_back(-0.0521961);
672   parameter[GLU_BB].push_back(1.11153);
673   parameter[GLU_BB].push_back(-0.491995);
674   parameter[GLU_BB].push_back(-1.86236);
675   parameter[GLU_BB].push_back(1.45332);
676   parameter[GLU_BB].push_back(-0.29708);
677 
678   parameter[GLU_SC1].push_back(8.544);
679   parameter[GLU_SC1].push_back(-0.249555);
680   parameter[GLU_SC1].push_back(12.8031);
681   parameter[GLU_SC1].push_back(-8.42696);
682   parameter[GLU_SC1].push_back(-4.66486);
683   parameter[GLU_SC1].push_back(4.90004);
684   parameter[GLU_SC1].push_back(-1.01204);
685 
686   parameter[GLY_BB].push_back(9.977);
687   parameter[GLY_BB].push_back(-0.0285799);
688   parameter[GLY_BB].push_back(1.84236);
689   parameter[GLY_BB].push_back(-0.0315192);
690   parameter[GLY_BB].push_back(-2.88326);
691   parameter[GLY_BB].push_back(1.87323);
692   parameter[GLY_BB].push_back(-0.345773);
693 
694   parameter[HIS_BB].push_back(10.721);
695   parameter[HIS_BB].push_back(-0.0379337);
696   parameter[HIS_BB].push_back(1.06028);
697   parameter[HIS_BB].push_back(-0.236143);
698   parameter[HIS_BB].push_back(-2.17819);
699   parameter[HIS_BB].push_back(1.58357);
700   parameter[HIS_BB].push_back(-0.31345);
701 
702   parameter[HIS_SC1].push_back(-0.424);
703   parameter[HIS_SC1].push_back(0.665176);
704   parameter[HIS_SC1].push_back(3.4369);
705   parameter[HIS_SC1].push_back(2.93795);
706   parameter[HIS_SC1].push_back(-5.18288);
707   parameter[HIS_SC1].push_back(2.12381);
708   parameter[HIS_SC1].push_back(-0.284224);
709 
710   parameter[HIS_SC2].push_back(5.363);
711   parameter[HIS_SC2].push_back(-0.0176945);
712   parameter[HIS_SC2].push_back(2.9506);
713   parameter[HIS_SC2].push_back(-0.387018);
714   parameter[HIS_SC2].push_back(-1.83951);
715   parameter[HIS_SC2].push_back(0.9703);
716   parameter[HIS_SC2].push_back(-0.1458);
717 
718   parameter[HIS_SC3].push_back(5.784);
719   parameter[HIS_SC3].push_back(-0.0293129);
720   parameter[HIS_SC3].push_back(2.74167);
721   parameter[HIS_SC3].push_back(-0.520875);
722   parameter[HIS_SC3].push_back(-1.62949);
723   parameter[HIS_SC3].push_back(0.902379);
724   parameter[HIS_SC3].push_back(-0.139957);
725 
726   parameter[ILE_BB].push_back(10.699);
727   parameter[ILE_BB].push_back(-0.0188962);
728   parameter[ILE_BB].push_back(1.217);
729   parameter[ILE_BB].push_back(0.242481);
730   parameter[ILE_BB].push_back(-3.13898);
731   parameter[ILE_BB].push_back(2.07916);
732   parameter[ILE_BB].push_back(-0.392574);
733 
734   parameter[ILE_SC1].push_back(-4.448);
735   parameter[ILE_SC1].push_back(1.20996);
736   parameter[ILE_SC1].push_back(11.5141);
737   parameter[ILE_SC1].push_back(6.98895);
738   parameter[ILE_SC1].push_back(-19.1948);
739   parameter[ILE_SC1].push_back(9.89207);
740   parameter[ILE_SC1].push_back(-1.60877);
741 
742   parameter[LEU_BB].push_back(10.692);
743   parameter[LEU_BB].push_back(-0.0414917);
744   parameter[LEU_BB].push_back(1.1077);
745   parameter[LEU_BB].push_back(-0.288062);
746   parameter[LEU_BB].push_back(-2.17187);
747   parameter[LEU_BB].push_back(1.59879);
748   parameter[LEU_BB].push_back(-0.318545);
749 
750   parameter[LEU_SC1].push_back(-4.448);
751   parameter[LEU_SC1].push_back(2.1063);
752   parameter[LEU_SC1].push_back(6.72381);
753   parameter[LEU_SC1].push_back(14.6954);
754   parameter[LEU_SC1].push_back(-23.7197);
755   parameter[LEU_SC1].push_back(10.7247);
756   parameter[LEU_SC1].push_back(-1.59146);
757 
758   parameter[LYS_BB].push_back(10.706);
759   parameter[LYS_BB].push_back(-0.0468629);
760   parameter[LYS_BB].push_back(1.09477);
761   parameter[LYS_BB].push_back(-0.432751);
762   parameter[LYS_BB].push_back(-1.94335);
763   parameter[LYS_BB].push_back(1.49109);
764   parameter[LYS_BB].push_back(-0.302589);
765 
766   parameter[LYS_SC1].push_back(-2.796);
767   parameter[LYS_SC1].push_back(0.508044);
768   parameter[LYS_SC1].push_back(7.91436);
769   parameter[LYS_SC1].push_back(4.54097);
770   parameter[LYS_SC1].push_back(-10.8051);
771   parameter[LYS_SC1].push_back(4.96204);
772   parameter[LYS_SC1].push_back(-0.724414);
773 
774   parameter[LYS_SC2].push_back(3.070);
775   parameter[LYS_SC2].push_back(-0.0101448);
776   parameter[LYS_SC2].push_back(4.67994);
777   parameter[LYS_SC2].push_back(-0.792529);
778   parameter[LYS_SC2].push_back(-2.09142);
779   parameter[LYS_SC2].push_back(1.02933);
780   parameter[LYS_SC2].push_back(-0.137787);
781 
782   parameter[MET_BB].push_back(10.671);
783   parameter[MET_BB].push_back(-0.0433724);
784   parameter[MET_BB].push_back(1.13784);
785   parameter[MET_BB].push_back(-0.40768);
786   parameter[MET_BB].push_back(-2.00555);
787   parameter[MET_BB].push_back(1.51673);
788   parameter[MET_BB].push_back(-0.305547);
789 
790   parameter[MET_SC1].push_back(5.85);
791   parameter[MET_SC1].push_back(-0.0485798);
792   parameter[MET_SC1].push_back(17.0391);
793   parameter[MET_SC1].push_back(-3.65327);
794   parameter[MET_SC1].push_back(-13.174);
795   parameter[MET_SC1].push_back(8.68286);
796   parameter[MET_SC1].push_back(-1.56095);
797 
798   parameter[PHE_BB].push_back(10.741);
799   parameter[PHE_BB].push_back(-0.0317275);
800   parameter[PHE_BB].push_back(1.15599);
801   parameter[PHE_BB].push_back(0.0276187);
802   parameter[PHE_BB].push_back(-2.74757);
803   parameter[PHE_BB].push_back(1.88783);
804   parameter[PHE_BB].push_back(-0.363525);
805 
806   parameter[PHE_SC1].push_back(-0.636);
807   parameter[PHE_SC1].push_back(0.527882);
808   parameter[PHE_SC1].push_back(6.77612);
809   parameter[PHE_SC1].push_back(3.18508);
810   parameter[PHE_SC1].push_back(-8.92826);
811   parameter[PHE_SC1].push_back(4.29752);
812   parameter[PHE_SC1].push_back(-0.65187);
813 
814   parameter[PHE_SC2].push_back(-0.424);
815   parameter[PHE_SC2].push_back(0.389174);
816   parameter[PHE_SC2].push_back(4.11761);
817   parameter[PHE_SC2].push_back(2.29527);
818   parameter[PHE_SC2].push_back(-4.7652);
819   parameter[PHE_SC2].push_back(1.97023);
820   parameter[PHE_SC2].push_back(-0.262318);
821 
822   parameter[PHE_SC3].push_back(-0.424);
823   parameter[PHE_SC3].push_back(0.38927);
824   parameter[PHE_SC3].push_back(4.11708);
825   parameter[PHE_SC3].push_back(2.29623);
826   parameter[PHE_SC3].push_back(-4.76592);
827   parameter[PHE_SC3].push_back(1.97055);
828   parameter[PHE_SC3].push_back(-0.262381);
829 
830   parameter[PRO_BB].push_back(11.434);
831   parameter[PRO_BB].push_back(-0.033323);
832   parameter[PRO_BB].push_back(0.472014);
833   parameter[PRO_BB].push_back(-0.290854);
834   parameter[PRO_BB].push_back(-1.81409);
835   parameter[PRO_BB].push_back(1.39751);
836   parameter[PRO_BB].push_back(-0.280407);
837 
838   parameter[PRO_SC1].push_back(-2.796);
839   parameter[PRO_SC1].push_back(0.95668);
840   parameter[PRO_SC1].push_back(6.84197);
841   parameter[PRO_SC1].push_back(6.43774);
842   parameter[PRO_SC1].push_back(-12.5068);
843   parameter[PRO_SC1].push_back(5.64597);
844   parameter[PRO_SC1].push_back(-0.825206);
845 
846   parameter[SER_BB].push_back(10.699);
847   parameter[SER_BB].push_back(-0.0325828);
848   parameter[SER_BB].push_back(1.20329);
849   parameter[SER_BB].push_back(-0.0674351);
850   parameter[SER_BB].push_back(-2.60749);
851   parameter[SER_BB].push_back(1.80318);
852   parameter[SER_BB].push_back(-0.346803);
853 
854   parameter[SER_SC1].push_back(3.298);
855   parameter[SER_SC1].push_back(-0.0366801);
856   parameter[SER_SC1].push_back(5.11077);
857   parameter[SER_SC1].push_back(-1.46774);
858   parameter[SER_SC1].push_back(-1.48421);
859   parameter[SER_SC1].push_back(0.800326);
860   parameter[SER_SC1].push_back(-0.108314);
861 
862   parameter[THR_BB].push_back(10.697);
863   parameter[THR_BB].push_back(-0.0242955);
864   parameter[THR_BB].push_back(1.24671);
865   parameter[THR_BB].push_back(0.146423);
866   parameter[THR_BB].push_back(-2.97429);
867   parameter[THR_BB].push_back(1.97513);
868   parameter[THR_BB].push_back(-0.371479);
869 
870   parameter[THR_SC1].push_back(2.366);
871   parameter[THR_SC1].push_back(0.0297604);
872   parameter[THR_SC1].push_back(11.9216);
873   parameter[THR_SC1].push_back(-9.32503);
874   parameter[THR_SC1].push_back(1.9396);
875   parameter[THR_SC1].push_back(0.0804861);
876   parameter[THR_SC1].push_back(-0.0302721);
877 
878   parameter[TRP_BB].push_back(10.689);
879   parameter[TRP_BB].push_back(-0.0265879);
880   parameter[TRP_BB].push_back(1.17819);
881   parameter[TRP_BB].push_back(0.0386457);
882   parameter[TRP_BB].push_back(-2.75634);
883   parameter[TRP_BB].push_back(1.88065);
884   parameter[TRP_BB].push_back(-0.360217);
885 
886   parameter[TRP_SC1].push_back(0.084);
887   parameter[TRP_SC1].push_back(0.752407);
888   parameter[TRP_SC1].push_back(5.3802);
889   parameter[TRP_SC1].push_back(4.09281);
890   parameter[TRP_SC1].push_back(-9.28029);
891   parameter[TRP_SC1].push_back(4.45923);
892   parameter[TRP_SC1].push_back(-0.689008);
893 
894   parameter[TRP_SC2].push_back(5.739);
895   parameter[TRP_SC2].push_back(0.0298492);
896   parameter[TRP_SC2].push_back(4.60446);
897   parameter[TRP_SC2].push_back(1.34463);
898   parameter[TRP_SC2].push_back(-5.69968);
899   parameter[TRP_SC2].push_back(2.84924);
900   parameter[TRP_SC2].push_back(-0.433781);
901 
902   parameter[TRP_SC3].push_back(-0.424);
903   parameter[TRP_SC3].push_back(0.388576);
904   parameter[TRP_SC3].push_back(4.11859);
905   parameter[TRP_SC3].push_back(2.29485);
906   parameter[TRP_SC3].push_back(-4.76255);
907   parameter[TRP_SC3].push_back(1.96849);
908   parameter[TRP_SC3].push_back(-0.262015);
909 
910   parameter[TRP_SC4].push_back(-0.424);
911   parameter[TRP_SC4].push_back(0.387685);
912   parameter[TRP_SC4].push_back(4.12153);
913   parameter[TRP_SC4].push_back(2.29144);
914   parameter[TRP_SC4].push_back(-4.7589);
915   parameter[TRP_SC4].push_back(1.96686);
916   parameter[TRP_SC4].push_back(-0.261786);
917 
918   parameter[TYR_BB].push_back(10.689);
919   parameter[TYR_BB].push_back(-0.0193526);
920   parameter[TYR_BB].push_back(1.18241);
921   parameter[TYR_BB].push_back(0.207318);
922   parameter[TYR_BB].push_back(-3.0041);
923   parameter[TYR_BB].push_back(1.99335);
924   parameter[TYR_BB].push_back(-0.376482);
925 
926   parameter[TYR_SC1].push_back(-0.636);
927   parameter[TYR_SC1].push_back(0.528902);
928   parameter[TYR_SC1].push_back(6.78168);
929   parameter[TYR_SC1].push_back(3.17769);
930   parameter[TYR_SC1].push_back(-8.93667);
931   parameter[TYR_SC1].push_back(4.30692);
932   parameter[TYR_SC1].push_back(-0.653993);
933 
934   parameter[TYR_SC2].push_back(-0.424);
935   parameter[TYR_SC2].push_back(0.388811);
936   parameter[TYR_SC2].push_back(4.11851);
937   parameter[TYR_SC2].push_back(2.29545);
938   parameter[TYR_SC2].push_back(-4.7668);
939   parameter[TYR_SC2].push_back(1.97131);
940   parameter[TYR_SC2].push_back(-0.262534);
941 
942   parameter[TYR_SC3].push_back(4.526);
943   parameter[TYR_SC3].push_back(-0.00381305);
944   parameter[TYR_SC3].push_back(5.8567);
945   parameter[TYR_SC3].push_back(-0.214086);
946   parameter[TYR_SC3].push_back(-4.63649);
947   parameter[TYR_SC3].push_back(2.52869);
948   parameter[TYR_SC3].push_back(-0.39894);
949 
950   parameter[VAL_BB].push_back(10.691);
951   parameter[VAL_BB].push_back(-0.0162929);
952   parameter[VAL_BB].push_back(1.24446);
953   parameter[VAL_BB].push_back(0.307914);
954   parameter[VAL_BB].push_back(-3.27446);
955   parameter[VAL_BB].push_back(2.14788);
956   parameter[VAL_BB].push_back(-0.403259);
957 
958   parameter[VAL_SC1].push_back(-3.516);
959   parameter[VAL_SC1].push_back(1.62307);
960   parameter[VAL_SC1].push_back(5.43064);
961   parameter[VAL_SC1].push_back(9.28809);
962   parameter[VAL_SC1].push_back(-14.9927);
963   parameter[VAL_SC1].push_back(6.6133);
964   parameter[VAL_SC1].push_back(-0.964977);
965 
966   parameter[A_BB1].push_back(32.88500000);
967   parameter[A_BB1].push_back(0.08339900);
968   parameter[A_BB1].push_back(-7.36054400);
969   parameter[A_BB1].push_back(2.19220300);
970   parameter[A_BB1].push_back(-3.56523400);
971   parameter[A_BB1].push_back(2.33326900);
972   parameter[A_BB1].push_back(-0.39785500);
973 
974   parameter[A_BB2].push_back(3.80600000);
975   parameter[A_BB2].push_back(-0.10727600);
976   parameter[A_BB2].push_back(9.58854100);
977   parameter[A_BB2].push_back(-6.23740500);
978   parameter[A_BB2].push_back(-0.48267300);
979   parameter[A_BB2].push_back(1.14119500);
980   parameter[A_BB2].push_back(-0.21385600);
981 
982   parameter[A_BB3].push_back(3.59400000);
983   parameter[A_BB3].push_back(0.04537300);
984   parameter[A_BB3].push_back(9.59178900);
985   parameter[A_BB3].push_back(-1.29202200);
986   parameter[A_BB3].push_back(-7.10851000);
987   parameter[A_BB3].push_back(4.05571200);
988   parameter[A_BB3].push_back(-0.63372500);
989 
990   parameter[A_SC1].push_back(6.67100000);
991   parameter[A_SC1].push_back(-0.00855300);
992   parameter[A_SC1].push_back(1.63222400);
993   parameter[A_SC1].push_back(-0.06466200);
994   parameter[A_SC1].push_back(-1.48694200);
995   parameter[A_SC1].push_back(0.78544600);
996   parameter[A_SC1].push_back(-0.12083500);
997 
998   parameter[A_SC2].push_back(5.95100000);
999   parameter[A_SC2].push_back(-0.02606600);
1000   parameter[A_SC2].push_back(2.54399900);
1001   parameter[A_SC2].push_back(-0.48436900);
1002   parameter[A_SC2].push_back(-1.55357400);
1003   parameter[A_SC2].push_back(0.86466900);
1004   parameter[A_SC2].push_back(-0.13509000);
1005 
1006   parameter[A_SC3].push_back(11.39400000);
1007   parameter[A_SC3].push_back(0.00871300);
1008   parameter[A_SC3].push_back(-0.23891300);
1009   parameter[A_SC3].push_back(0.48919400);
1010   parameter[A_SC3].push_back(-1.75289400);
1011   parameter[A_SC3].push_back(0.99267500);
1012   parameter[A_SC3].push_back(-0.16291300);
1013 
1014   parameter[A_SC4].push_back(6.45900000);
1015   parameter[A_SC4].push_back(0.01990600);
1016   parameter[A_SC4].push_back(4.17970400);
1017   parameter[A_SC4].push_back(0.97629900);
1018   parameter[A_SC4].push_back(-5.03297800);
1019   parameter[A_SC4].push_back(2.55576700);
1020   parameter[A_SC4].push_back(-0.39150500);
1021 
1022   parameter[A_3TE].push_back(4.23000000);
1023   parameter[A_3TE].push_back(0.00064800);
1024   parameter[A_3TE].push_back(0.92124600);
1025   parameter[A_3TE].push_back(0.08064300);
1026   parameter[A_3TE].push_back(-0.39054400);
1027   parameter[A_3TE].push_back(0.12429100);
1028   parameter[A_3TE].push_back(-0.01122700);
1029 
1030   parameter[A_5TE].push_back(4.23000000);
1031   parameter[A_5TE].push_back(0.00039300);
1032   parameter[A_5TE].push_back(0.92305100);
1033   parameter[A_5TE].push_back(0.07747500);
1034   parameter[A_5TE].push_back(-0.38792100);
1035   parameter[A_5TE].push_back(0.12323800);
1036   parameter[A_5TE].push_back(-0.01106600);
1037 
1038   parameter[A_TE3].push_back(7.82400000);
1039   parameter[A_TE3].push_back(-0.04881000);
1040   parameter[A_TE3].push_back(8.21557900);
1041   parameter[A_TE3].push_back(-0.89491400);
1042   parameter[A_TE3].push_back(-9.54293700);
1043   parameter[A_TE3].push_back(6.33122200);
1044   parameter[A_TE3].push_back(-1.16672900);
1045 
1046   parameter[A_TE5].push_back(8.03600000);
1047   parameter[A_TE5].push_back(0.01641200);
1048   parameter[A_TE5].push_back(5.14902200);
1049   parameter[A_TE5].push_back(0.83419700);
1050   parameter[A_TE5].push_back(-7.59068300);
1051   parameter[A_TE5].push_back(4.52063200);
1052   parameter[A_TE5].push_back(-0.78260800);
1053 
1054   parameter[C_BB1].push_back(32.88500000);
1055   parameter[C_BB1].push_back(0.08311100);
1056   parameter[C_BB1].push_back(-7.35432100);
1057   parameter[C_BB1].push_back(2.18610000);
1058   parameter[C_BB1].push_back(-3.55788300);
1059   parameter[C_BB1].push_back(2.32918700);
1060   parameter[C_BB1].push_back(-0.39720000);
1061 
1062   parameter[C_BB2].push_back(3.80600000);
1063   parameter[C_BB2].push_back(-0.10808100);
1064   parameter[C_BB2].push_back(9.61612600);
1065   parameter[C_BB2].push_back(-6.28595400);
1066   parameter[C_BB2].push_back(-0.45187000);
1067   parameter[C_BB2].push_back(1.13326000);
1068   parameter[C_BB2].push_back(-0.21320300);
1069 
1070   parameter[C_BB3].push_back(3.59400000);
1071   parameter[C_BB3].push_back(0.04484200);
1072   parameter[C_BB3].push_back(9.61919800);
1073   parameter[C_BB3].push_back(-1.33582800);
1074   parameter[C_BB3].push_back(-7.07200400);
1075   parameter[C_BB3].push_back(4.03952900);
1076   parameter[C_BB3].push_back(-0.63098200);
1077 
1078   parameter[C_SC1].push_back(5.95100000);
1079   parameter[C_SC1].push_back(-0.02911300);
1080   parameter[C_SC1].push_back(2.59700400);
1081   parameter[C_SC1].push_back(-0.55507700);
1082   parameter[C_SC1].push_back(-1.56344600);
1083   parameter[C_SC1].push_back(0.88956200);
1084   parameter[C_SC1].push_back(-0.14061300);
1085 
1086   parameter[C_SC2].push_back(11.62100000);
1087   parameter[C_SC2].push_back(0.01366100);
1088   parameter[C_SC2].push_back(-0.25959200);
1089   parameter[C_SC2].push_back(0.48918300);
1090   parameter[C_SC2].push_back(-1.52550500);
1091   parameter[C_SC2].push_back(0.83644100);
1092   parameter[C_SC2].push_back(-0.13407300);
1093 
1094   parameter[C_SC3].push_back(5.01900000);
1095   parameter[C_SC3].push_back(-0.03276100);
1096   parameter[C_SC3].push_back(5.53776900);
1097   parameter[C_SC3].push_back(-0.95105000);
1098   parameter[C_SC3].push_back(-3.71130800);
1099   parameter[C_SC3].push_back(2.16146000);
1100   parameter[C_SC3].push_back(-0.34918600);
1101 
1102   parameter[C_3TE].push_back(4.23000000);
1103   parameter[C_3TE].push_back(0.00057300);
1104   parameter[C_3TE].push_back(0.92174800);
1105   parameter[C_3TE].push_back(0.07964500);
1106   parameter[C_3TE].push_back(-0.38965700);
1107   parameter[C_3TE].push_back(0.12392500);
1108   parameter[C_3TE].push_back(-0.01117000);
1109 
1110   parameter[C_5TE].push_back(4.23000000);
1111   parameter[C_5TE].push_back(0.00071000);
1112   parameter[C_5TE].push_back(0.92082800);
1113   parameter[C_5TE].push_back(0.08150600);
1114   parameter[C_5TE].push_back(-0.39127000);
1115   parameter[C_5TE].push_back(0.12455900);
1116   parameter[C_5TE].push_back(-0.01126300);
1117 
1118   parameter[C_TE3].push_back(7.82400000);
1119   parameter[C_TE3].push_back(-0.05848300);
1120   parameter[C_TE3].push_back(8.29319900);
1121   parameter[C_TE3].push_back(-1.12563800);
1122   parameter[C_TE3].push_back(-9.42197600);
1123   parameter[C_TE3].push_back(6.35441700);
1124   parameter[C_TE3].push_back(-1.18356900);
1125 
1126   parameter[C_TE5].push_back(8.03600000);
1127   parameter[C_TE5].push_back(0.00493500);
1128   parameter[C_TE5].push_back(4.92622000);
1129   parameter[C_TE5].push_back(0.64810700);
1130   parameter[C_TE5].push_back(-7.05100000);
1131   parameter[C_TE5].push_back(4.26064400);
1132   parameter[C_TE5].push_back(-0.74819100);
1133 
1134   parameter[G_BB1].push_back(32.88500000);
1135   parameter[G_BB1].push_back(0.08325400);
1136   parameter[G_BB1].push_back(-7.35736000);
1137   parameter[G_BB1].push_back(2.18914800);
1138   parameter[G_BB1].push_back(-3.56154800);
1139   parameter[G_BB1].push_back(2.33120600);
1140   parameter[G_BB1].push_back(-0.39752300);
1141 
1142   parameter[G_BB2].push_back(3.80600000);
1143   parameter[G_BB2].push_back(-0.10788300);
1144   parameter[G_BB2].push_back(9.60930800);
1145   parameter[G_BB2].push_back(-6.27402500);
1146   parameter[G_BB2].push_back(-0.46192700);
1147   parameter[G_BB2].push_back(1.13737000);
1148   parameter[G_BB2].push_back(-0.21383100);
1149 
1150   parameter[G_BB3].push_back(3.59400000);
1151   parameter[G_BB3].push_back(0.04514500);
1152   parameter[G_BB3].push_back(9.61234700);
1153   parameter[G_BB3].push_back(-1.31542100);
1154   parameter[G_BB3].push_back(-7.09150500);
1155   parameter[G_BB3].push_back(4.04706200);
1156   parameter[G_BB3].push_back(-0.63201000);
1157 
1158   parameter[G_SC1].push_back(6.67100000);
1159   parameter[G_SC1].push_back(-0.00863200);
1160   parameter[G_SC1].push_back(1.63252300);
1161   parameter[G_SC1].push_back(-0.06567200);
1162   parameter[G_SC1].push_back(-1.48680500);
1163   parameter[G_SC1].push_back(0.78565600);
1164   parameter[G_SC1].push_back(-0.12088900);
1165 
1166   parameter[G_SC2].push_back(11.39400000);
1167   parameter[G_SC2].push_back(0.00912200);
1168   parameter[G_SC2].push_back(-0.22869000);
1169   parameter[G_SC2].push_back(0.49616400);
1170   parameter[G_SC2].push_back(-1.75039000);
1171   parameter[G_SC2].push_back(0.98649200);
1172   parameter[G_SC2].push_back(-0.16141600);
1173 
1174   parameter[G_SC3].push_back(10.90100000);
1175   parameter[G_SC3].push_back(0.02208700);
1176   parameter[G_SC3].push_back(0.17032800);
1177   parameter[G_SC3].push_back(0.73280800);
1178   parameter[G_SC3].push_back(-1.95292000);
1179   parameter[G_SC3].push_back(0.98357600);
1180   parameter[G_SC3].push_back(-0.14790900);
1181 
1182   parameter[G_SC4].push_back(6.45900000);
1183   parameter[G_SC4].push_back(0.02023700);
1184   parameter[G_SC4].push_back(4.17655400);
1185   parameter[G_SC4].push_back(0.98731800);
1186   parameter[G_SC4].push_back(-5.04352800);
1187   parameter[G_SC4].push_back(2.56059400);
1188   parameter[G_SC4].push_back(-0.39234300);
1189 
1190   parameter[G_3TE].push_back(4.23000000);
1191   parameter[G_3TE].push_back(0.00066300);
1192   parameter[G_3TE].push_back(0.92118800);
1193   parameter[G_3TE].push_back(0.08062700);
1194   parameter[G_3TE].push_back(-0.39041600);
1195   parameter[G_3TE].push_back(0.12419400);
1196   parameter[G_3TE].push_back(-0.01120500);
1197 
1198   parameter[G_5TE].push_back(4.23000000);
1199   parameter[G_5TE].push_back(0.00062800);
1200   parameter[G_5TE].push_back(0.92133500);
1201   parameter[G_5TE].push_back(0.08029900);
1202   parameter[G_5TE].push_back(-0.39015300);
1203   parameter[G_5TE].push_back(0.12411600);
1204   parameter[G_5TE].push_back(-0.01119900);
1205 
1206   parameter[G_TE3].push_back(7.82400000);
1207   parameter[G_TE3].push_back(-0.05177400);
1208   parameter[G_TE3].push_back(8.34606700);
1209   parameter[G_TE3].push_back(-1.02936300);
1210   parameter[G_TE3].push_back(-9.55211900);
1211   parameter[G_TE3].push_back(6.37776600);
1212   parameter[G_TE3].push_back(-1.17898000);
1213 
1214   parameter[G_TE5].push_back(8.03600000);
1215   parameter[G_TE5].push_back(0.00525100);
1216   parameter[G_TE5].push_back(4.71070600);
1217   parameter[G_TE5].push_back(0.66746900);
1218   parameter[G_TE5].push_back(-6.72538700);
1219   parameter[G_TE5].push_back(4.03644100);
1220   parameter[G_TE5].push_back(-0.70605700);
1221 
1222   parameter[U_BB1].push_back(32.88500000);
1223   parameter[U_BB1].push_back(0.08321400);
1224   parameter[U_BB1].push_back(-7.35634900);
1225   parameter[U_BB1].push_back(2.18826800);
1226   parameter[U_BB1].push_back(-3.56047400);
1227   parameter[U_BB1].push_back(2.33064700);
1228   parameter[U_BB1].push_back(-0.39744000);
1229 
1230   parameter[U_BB2].push_back(3.80600000);
1231   parameter[U_BB2].push_back(-0.10773100);
1232   parameter[U_BB2].push_back(9.60099900);
1233   parameter[U_BB2].push_back(-6.26131900);
1234   parameter[U_BB2].push_back(-0.46668300);
1235   parameter[U_BB2].push_back(1.13698100);
1236   parameter[U_BB2].push_back(-0.21351600);
1237 
1238   parameter[U_BB3].push_back(3.59400000);
1239   parameter[U_BB3].push_back(0.04544300);
1240   parameter[U_BB3].push_back(9.59625900);
1241   parameter[U_BB3].push_back(-1.29222200);
1242   parameter[U_BB3].push_back(-7.11143200);
1243   parameter[U_BB3].push_back(4.05687700);
1244   parameter[U_BB3].push_back(-0.63382800);
1245 
1246   parameter[U_SC1].push_back(5.95100000);
1247   parameter[U_SC1].push_back(-0.02924500);
1248   parameter[U_SC1].push_back(2.59668700);
1249   parameter[U_SC1].push_back(-0.56118700);
1250   parameter[U_SC1].push_back(-1.56477100);
1251   parameter[U_SC1].push_back(0.89265100);
1252   parameter[U_SC1].push_back(-0.14130800);
1253 
1254   parameter[U_SC2].push_back(10.90100000);
1255   parameter[U_SC2].push_back(0.02178900);
1256   parameter[U_SC2].push_back(0.18839000);
1257   parameter[U_SC2].push_back(0.72223100);
1258   parameter[U_SC2].push_back(-1.92581600);
1259   parameter[U_SC2].push_back(0.96654300);
1260   parameter[U_SC2].push_back(-0.14501300);
1261 
1262   parameter[U_SC3].push_back(5.24600000);
1263   parameter[U_SC3].push_back(-0.04586500);
1264   parameter[U_SC3].push_back(5.89978100);
1265   parameter[U_SC3].push_back(-1.50664700);
1266   parameter[U_SC3].push_back(-3.17054400);
1267   parameter[U_SC3].push_back(1.93717100);
1268   parameter[U_SC3].push_back(-0.31701000);
1269 
1270   parameter[U_3TE].push_back(4.23000000);
1271   parameter[U_3TE].push_back(0.00067500);
1272   parameter[U_3TE].push_back(0.92102300);
1273   parameter[U_3TE].push_back(0.08100800);
1274   parameter[U_3TE].push_back(-0.39084300);
1275   parameter[U_3TE].push_back(0.12441900);
1276   parameter[U_3TE].push_back(-0.01124900);
1277 
1278   parameter[U_5TE].push_back(4.23000000);
1279   parameter[U_5TE].push_back(0.00059000);
1280   parameter[U_5TE].push_back(0.92154600);
1281   parameter[U_5TE].push_back(0.07968200);
1282   parameter[U_5TE].push_back(-0.38950100);
1283   parameter[U_5TE].push_back(0.12382500);
1284   parameter[U_5TE].push_back(-0.01115100);
1285 
1286   parameter[U_TE3].push_back(7.82400000);
1287   parameter[U_TE3].push_back(-0.02968100);
1288   parameter[U_TE3].push_back(7.93783200);
1289   parameter[U_TE3].push_back(-0.33078100);
1290   parameter[U_TE3].push_back(-10.14120200);
1291   parameter[U_TE3].push_back(6.63334700);
1292   parameter[U_TE3].push_back(-1.22111200);
1293 
1294   parameter[U_TE5].push_back(8.03600000);
1295   parameter[U_TE5].push_back(-0.00909700);
1296   parameter[U_TE5].push_back(4.33193500);
1297   parameter[U_TE5].push_back(0.43416500);
1298   parameter[U_TE5].push_back(-5.80831400);
1299   parameter[U_TE5].push_back(3.52438800);
1300   parameter[U_TE5].push_back(-0.62382400);
1301 
1302   parameter[DA_BB1].push_back(32.88500000);
1303   parameter[DA_BB1].push_back(0.08179900);
1304   parameter[DA_BB1].push_back(-7.31735900);
1305   parameter[DA_BB1].push_back(2.15614500);
1306   parameter[DA_BB1].push_back(-3.52263200);
1307   parameter[DA_BB1].push_back(2.30604700);
1308   parameter[DA_BB1].push_back(-0.39270100);
1309 
1310   parameter[DA_BB2].push_back(3.80600000);
1311   parameter[DA_BB2].push_back(-0.10597700);
1312   parameter[DA_BB2].push_back(9.52537500);
1313   parameter[DA_BB2].push_back(-6.12991000);
1314   parameter[DA_BB2].push_back(-0.54092600);
1315   parameter[DA_BB2].push_back(1.15429100);
1316   parameter[DA_BB2].push_back(-0.21503500);
1317 
1318   parameter[DA_BB3].push_back(-1.35600000);
1319   parameter[DA_BB3].push_back(0.58928300);
1320   parameter[DA_BB3].push_back(6.71894100);
1321   parameter[DA_BB3].push_back(4.14050900);
1322   parameter[DA_BB3].push_back(-9.65859900);
1323   parameter[DA_BB3].push_back(4.43185000);
1324   parameter[DA_BB3].push_back(-0.64657300);
1325 
1326   parameter[DA_SC1].push_back(6.67100000);
1327   parameter[DA_SC1].push_back(-0.00871400);
1328   parameter[DA_SC1].push_back(1.63289100);
1329   parameter[DA_SC1].push_back(-0.06637700);
1330   parameter[DA_SC1].push_back(-1.48632900);
1331   parameter[DA_SC1].push_back(0.78551800);
1332   parameter[DA_SC1].push_back(-0.12087300);
1333 
1334   parameter[DA_SC2].push_back(5.95100000);
1335   parameter[DA_SC2].push_back(-0.02634300);
1336   parameter[DA_SC2].push_back(2.54864300);
1337   parameter[DA_SC2].push_back(-0.49015800);
1338   parameter[DA_SC2].push_back(-1.55386900);
1339   parameter[DA_SC2].push_back(0.86630200);
1340   parameter[DA_SC2].push_back(-0.13546200);
1341 
1342   parameter[DA_SC3].push_back(11.39400000);
1343   parameter[DA_SC3].push_back(0.00859500);
1344   parameter[DA_SC3].push_back(-0.25471400);
1345   parameter[DA_SC3].push_back(0.48718800);
1346   parameter[DA_SC3].push_back(-1.74520000);
1347   parameter[DA_SC3].push_back(0.99246200);
1348   parameter[DA_SC3].push_back(-0.16351900);
1349 
1350   parameter[DA_SC4].push_back(6.45900000);
1351   parameter[DA_SC4].push_back(0.01991800);
1352   parameter[DA_SC4].push_back(4.17962300);
1353   parameter[DA_SC4].push_back(0.97469100);
1354   parameter[DA_SC4].push_back(-5.02950400);
1355   parameter[DA_SC4].push_back(2.55371800);
1356   parameter[DA_SC4].push_back(-0.39113400);
1357 
1358   parameter[DA_3TE].push_back(4.23000000);
1359   parameter[DA_3TE].push_back(0.00062600);
1360   parameter[DA_3TE].push_back(0.92142000);
1361   parameter[DA_3TE].push_back(0.08016400);
1362   parameter[DA_3TE].push_back(-0.39000300);
1363   parameter[DA_3TE].push_back(0.12402500);
1364   parameter[DA_3TE].push_back(-0.01117900);
1365 
1366   parameter[DA_5TE].push_back(4.23000000);
1367   parameter[DA_5TE].push_back(0.00055500);
1368   parameter[DA_5TE].push_back(0.92183900);
1369   parameter[DA_5TE].push_back(0.07907600);
1370   parameter[DA_5TE].push_back(-0.38895100);
1371   parameter[DA_5TE].push_back(0.12359600);
1372   parameter[DA_5TE].push_back(-0.01111600);
1373 
1374   parameter[DA_TE3].push_back(2.87400000);
1375   parameter[DA_TE3].push_back(0.00112900);
1376   parameter[DA_TE3].push_back(12.51167200);
1377   parameter[DA_TE3].push_back(-7.67548000);
1378   parameter[DA_TE3].push_back(-2.02234000);
1379   parameter[DA_TE3].push_back(2.50837100);
1380   parameter[DA_TE3].push_back(-0.49458500);
1381 
1382   parameter[DA_TE5].push_back(8.03600000);
1383   parameter[DA_TE5].push_back(0.00473100);
1384   parameter[DA_TE5].push_back(4.65554400);
1385   parameter[DA_TE5].push_back(0.66424100);
1386   parameter[DA_TE5].push_back(-6.62131300);
1387   parameter[DA_TE5].push_back(3.96107400);
1388   parameter[DA_TE5].push_back(-0.69075800);
1389 
1390   parameter[DC_BB1].push_back(32.88500000);
1391   parameter[DC_BB1].push_back(0.08189900);
1392   parameter[DC_BB1].push_back(-7.32493500);
1393   parameter[DC_BB1].push_back(2.15976900);
1394   parameter[DC_BB1].push_back(-3.52612100);
1395   parameter[DC_BB1].push_back(2.31058600);
1396   parameter[DC_BB1].push_back(-0.39402700);
1397 
1398   parameter[DC_BB2].push_back(3.80600000);
1399   parameter[DC_BB2].push_back(-0.10559800);
1400   parameter[DC_BB2].push_back(9.52527700);
1401   parameter[DC_BB2].push_back(-6.12131700);
1402   parameter[DC_BB2].push_back(-0.54899400);
1403   parameter[DC_BB2].push_back(1.15592900);
1404   parameter[DC_BB2].push_back(-0.21494500);
1405 
1406   parameter[DC_BB3].push_back(-1.35600000);
1407   parameter[DC_BB3].push_back(0.55525700);
1408   parameter[DC_BB3].push_back(6.80305500);
1409   parameter[DC_BB3].push_back(4.05924700);
1410   parameter[DC_BB3].push_back(-9.61034700);
1411   parameter[DC_BB3].push_back(4.41253800);
1412   parameter[DC_BB3].push_back(-0.64315100);
1413 
1414   parameter[DC_SC1].push_back(5.95100000);
1415   parameter[DC_SC1].push_back(-0.02899900);
1416   parameter[DC_SC1].push_back(2.59587800);
1417   parameter[DC_SC1].push_back(-0.55388300);
1418   parameter[DC_SC1].push_back(-1.56395100);
1419   parameter[DC_SC1].push_back(0.88967400);
1420   parameter[DC_SC1].push_back(-0.14062500);
1421 
1422   parameter[DC_SC2].push_back(11.62100000);
1423   parameter[DC_SC2].push_back(0.01358100);
1424   parameter[DC_SC2].push_back(-0.24913000);
1425   parameter[DC_SC2].push_back(0.48787200);
1426   parameter[DC_SC2].push_back(-1.52867300);
1427   parameter[DC_SC2].push_back(0.83694900);
1428   parameter[DC_SC2].push_back(-0.13395300);
1429 
1430   parameter[DC_SC3].push_back(5.01900000);
1431   parameter[DC_SC3].push_back(-0.03298400);
1432   parameter[DC_SC3].push_back(5.54242800);
1433   parameter[DC_SC3].push_back(-0.96081500);
1434   parameter[DC_SC3].push_back(-3.71051600);
1435   parameter[DC_SC3].push_back(2.16500200);
1436   parameter[DC_SC3].push_back(-0.35023400);
1437 
1438   parameter[DC_3TE].push_back(4.23000000);
1439   parameter[DC_3TE].push_back(0.00055700);
1440   parameter[DC_3TE].push_back(0.92181400);
1441   parameter[DC_3TE].push_back(0.07924000);
1442   parameter[DC_3TE].push_back(-0.38916400);
1443   parameter[DC_3TE].push_back(0.12369900);
1444   parameter[DC_3TE].push_back(-0.01113300);
1445 
1446   parameter[DC_5TE].push_back(4.23000000);
1447   parameter[DC_5TE].push_back(0.00066500);
1448   parameter[DC_5TE].push_back(0.92103900);
1449   parameter[DC_5TE].push_back(0.08064600);
1450   parameter[DC_5TE].push_back(-0.39034900);
1451   parameter[DC_5TE].push_back(0.12417600);
1452   parameter[DC_5TE].push_back(-0.01120600);
1453 
1454   parameter[DC_TE3].push_back(2.87400000);
1455   parameter[DC_TE3].push_back(-0.05235500);
1456   parameter[DC_TE3].push_back(13.09201200);
1457   parameter[DC_TE3].push_back(-9.48128200);
1458   parameter[DC_TE3].push_back(-0.14958600);
1459   parameter[DC_TE3].push_back(1.75537200);
1460   parameter[DC_TE3].push_back(-0.39347500);
1461 
1462   parameter[DC_TE5].push_back(8.03600000);
1463   parameter[DC_TE5].push_back(-0.00513600);
1464   parameter[DC_TE5].push_back(4.67705700);
1465   parameter[DC_TE5].push_back(0.48333300);
1466   parameter[DC_TE5].push_back(-6.34511000);
1467   parameter[DC_TE5].push_back(3.83388500);
1468   parameter[DC_TE5].push_back(-0.67367800);
1469 
1470   parameter[DG_BB1].push_back(32.88500000);
1471   parameter[DG_BB1].push_back(0.08182900);
1472   parameter[DG_BB1].push_back(-7.32133900);
1473   parameter[DG_BB1].push_back(2.15767900);
1474   parameter[DG_BB1].push_back(-3.52369700);
1475   parameter[DG_BB1].push_back(2.30839600);
1476   parameter[DG_BB1].push_back(-0.39348300);
1477 
1478   parameter[DG_BB2].push_back(3.80600000);
1479   parameter[DG_BB2].push_back(-0.10618100);
1480   parameter[DG_BB2].push_back(9.54169000);
1481   parameter[DG_BB2].push_back(-6.15177600);
1482   parameter[DG_BB2].push_back(-0.53462400);
1483   parameter[DG_BB2].push_back(1.15581300);
1484   parameter[DG_BB2].push_back(-0.21567000);
1485 
1486   parameter[DG_BB3].push_back(-1.35600000);
1487   parameter[DG_BB3].push_back(0.57489100);
1488   parameter[DG_BB3].push_back(6.75164700);
1489   parameter[DG_BB3].push_back(4.11300900);
1490   parameter[DG_BB3].push_back(-9.63394600);
1491   parameter[DG_BB3].push_back(4.41675400);
1492   parameter[DG_BB3].push_back(-0.64339900);
1493 
1494   parameter[DG_SC1].push_back(6.67100000);
1495   parameter[DG_SC1].push_back(-0.00886600);
1496   parameter[DG_SC1].push_back(1.63333000);
1497   parameter[DG_SC1].push_back(-0.06892100);
1498   parameter[DG_SC1].push_back(-1.48683500);
1499   parameter[DG_SC1].push_back(0.78670800);
1500   parameter[DG_SC1].push_back(-0.12113900);
1501 
1502   parameter[DG_SC2].push_back(11.39400000);
1503   parameter[DG_SC2].push_back(0.00907900);
1504   parameter[DG_SC2].push_back(-0.22475500);
1505   parameter[DG_SC2].push_back(0.49535100);
1506   parameter[DG_SC2].push_back(-1.75324900);
1507   parameter[DG_SC2].push_back(0.98767400);
1508   parameter[DG_SC2].push_back(-0.16150800);
1509 
1510   parameter[DG_SC3].push_back(10.90100000);
1511   parameter[DG_SC3].push_back(0.02207600);
1512   parameter[DG_SC3].push_back(0.17932200);
1513   parameter[DG_SC3].push_back(0.73253200);
1514   parameter[DG_SC3].push_back(-1.95554900);
1515   parameter[DG_SC3].push_back(0.98339900);
1516   parameter[DG_SC3].push_back(-0.14763600);
1517 
1518   parameter[DG_SC4].push_back(6.45900000);
1519   parameter[DG_SC4].push_back(0.02018400);
1520   parameter[DG_SC4].push_back(4.17705400);
1521   parameter[DG_SC4].push_back(0.98531700);
1522   parameter[DG_SC4].push_back(-5.04354900);
1523   parameter[DG_SC4].push_back(2.56123700);
1524   parameter[DG_SC4].push_back(-0.39249300);
1525 
1526   parameter[DG_3TE].push_back(4.23000000);
1527   parameter[DG_3TE].push_back(0.00061700);
1528   parameter[DG_3TE].push_back(0.92140100);
1529   parameter[DG_3TE].push_back(0.08016400);
1530   parameter[DG_3TE].push_back(-0.39003500);
1531   parameter[DG_3TE].push_back(0.12406900);
1532   parameter[DG_3TE].push_back(-0.01119200);
1533 
1534   parameter[DG_5TE].push_back(4.23000000);
1535   parameter[DG_5TE].push_back(0.00064900);
1536   parameter[DG_5TE].push_back(0.92110500);
1537   parameter[DG_5TE].push_back(0.08031500);
1538   parameter[DG_5TE].push_back(-0.38997000);
1539   parameter[DG_5TE].push_back(0.12401200);
1540   parameter[DG_5TE].push_back(-0.01118100);
1541 
1542   parameter[DG_TE3].push_back(2.87400000);
1543   parameter[DG_TE3].push_back(0.00182000);
1544   parameter[DG_TE3].push_back(12.41507000);
1545   parameter[DG_TE3].push_back(-7.47384800);
1546   parameter[DG_TE3].push_back(-2.11864700);
1547   parameter[DG_TE3].push_back(2.50112600);
1548   parameter[DG_TE3].push_back(-0.48652200);
1549 
1550   parameter[DG_TE5].push_back(8.03600000);
1551   parameter[DG_TE5].push_back(0.00676400);
1552   parameter[DG_TE5].push_back(4.65989200);
1553   parameter[DG_TE5].push_back(0.78482500);
1554   parameter[DG_TE5].push_back(-6.86460600);
1555   parameter[DG_TE5].push_back(4.11675400);
1556   parameter[DG_TE5].push_back(-0.72249100);
1557 
1558   parameter[DT_BB1].push_back(32.88500000);
1559   parameter[DT_BB1].push_back(0.08220100);
1560   parameter[DT_BB1].push_back(-7.33006800);
1561   parameter[DT_BB1].push_back(2.16636500);
1562   parameter[DT_BB1].push_back(-3.53465700);
1563   parameter[DT_BB1].push_back(2.31447600);
1564   parameter[DT_BB1].push_back(-0.39445400);
1565 
1566   parameter[DT_BB2].push_back(3.80600000);
1567   parameter[DT_BB2].push_back(-0.10723000);
1568   parameter[DT_BB2].push_back(9.56675000);
1569   parameter[DT_BB2].push_back(-6.20236100);
1570   parameter[DT_BB2].push_back(-0.49550400);
1571   parameter[DT_BB2].push_back(1.14300600);
1572   parameter[DT_BB2].push_back(-0.21420000);
1573 
1574   parameter[DT_BB3].push_back(-1.35600000);
1575   parameter[DT_BB3].push_back(0.56737900);
1576   parameter[DT_BB3].push_back(6.76595400);
1577   parameter[DT_BB3].push_back(4.08976100);
1578   parameter[DT_BB3].push_back(-9.61512500);
1579   parameter[DT_BB3].push_back(4.40975100);
1580   parameter[DT_BB3].push_back(-0.64239800);
1581 
1582   parameter[DT_SC1].push_back(5.95100000);
1583   parameter[DT_SC1].push_back(-0.02926500);
1584   parameter[DT_SC1].push_back(2.59630300);
1585   parameter[DT_SC1].push_back(-0.56152200);
1586   parameter[DT_SC1].push_back(-1.56532600);
1587   parameter[DT_SC1].push_back(0.89322800);
1588   parameter[DT_SC1].push_back(-0.14142900);
1589 
1590   parameter[DT_SC2].push_back(10.90100000);
1591   parameter[DT_SC2].push_back(0.02183400);
1592   parameter[DT_SC2].push_back(0.19463000);
1593   parameter[DT_SC2].push_back(0.72393000);
1594   parameter[DT_SC2].push_back(-1.93199500);
1595   parameter[DT_SC2].push_back(0.96856300);
1596   parameter[DT_SC2].push_back(-0.14512600);
1597 
1598   parameter[DT_SC3].push_back(4.31400000);
1599   parameter[DT_SC3].push_back(-0.07745600);
1600   parameter[DT_SC3].push_back(12.49820300);
1601   parameter[DT_SC3].push_back(-7.64994200);
1602   parameter[DT_SC3].push_back(-3.00359600);
1603   parameter[DT_SC3].push_back(3.26263300);
1604   parameter[DT_SC3].push_back(-0.64498600);
1605 
1606   parameter[DT_3TE].push_back(4.23000000);
1607   parameter[DT_3TE].push_back(0.00062000);
1608   parameter[DT_3TE].push_back(0.92141100);
1609   parameter[DT_3TE].push_back(0.08030900);
1610   parameter[DT_3TE].push_back(-0.39021500);
1611   parameter[DT_3TE].push_back(0.12414000);
1612   parameter[DT_3TE].push_back(-0.01120100);
1613 
1614   parameter[DT_5TE].push_back(4.23000000);
1615   parameter[DT_5TE].push_back(0.00063700);
1616   parameter[DT_5TE].push_back(0.92130800);
1617   parameter[DT_5TE].push_back(0.08026900);
1618   parameter[DT_5TE].push_back(-0.39007500);
1619   parameter[DT_5TE].push_back(0.12406600);
1620   parameter[DT_5TE].push_back(-0.01118800);
1621 
1622   parameter[DT_TE3].push_back(2.87400000);
1623   parameter[DT_TE3].push_back(-0.00251200);
1624   parameter[DT_TE3].push_back(12.43576400);
1625   parameter[DT_TE3].push_back(-7.55343800);
1626   parameter[DT_TE3].push_back(-2.07363500);
1627   parameter[DT_TE3].push_back(2.51279300);
1628   parameter[DT_TE3].push_back(-0.49437100);
1629 
1630   parameter[DT_TE5].push_back(8.03600000);
1631   parameter[DT_TE5].push_back(0.00119900);
1632   parameter[DT_TE5].push_back(4.91762300);
1633   parameter[DT_TE5].push_back(0.65637000);
1634   parameter[DT_TE5].push_back(-7.23392500);
1635   parameter[DT_TE5].push_back(4.44636600);
1636   parameter[DT_TE5].push_back(-0.79467800);
1637 
1638   auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
1639   if( moldat ) {
1640     log<<"  MOLINFO DATA found with label " <<moldat->getLabel()<<", using proper atom names\n";
1641     for(unsigned i=0; i<atoms.size(); ++i) {
1642       string Aname = moldat->getAtomName(atoms[i]);
1643       string Rname = moldat->getResidueName(atoms[i]);
1644       if(Rname=="ALA") {
1645         if(Aname=="BB") {
1646           atoi[i]=ALA_BB;
1647         } else error("Atom name not known: "+Aname);
1648       } else if(Rname=="ARG") {
1649         if(Aname=="BB") {
1650           atoi[i]=ARG_BB;
1651         } else if(Aname=="SC1") {
1652           atoi[i]=ARG_SC1;
1653         } else if(Aname=="SC2") {
1654           atoi[i]=ARG_SC2;
1655         } else error("Atom name not known: "+Aname);
1656       } else if(Rname=="ASN") {
1657         if(Aname=="BB") {
1658           atoi[i]=ASN_BB;
1659         } else if(Aname=="SC1") {
1660           atoi[i]=ASN_SC1;
1661         } else error("Atom name not known: "+Aname);
1662       } else if(Rname=="ASP") {
1663         if(Aname=="BB") {
1664           atoi[i]=ASP_BB;
1665         } else if(Aname=="SC1") {
1666           atoi[i]=ASP_SC1;
1667         } else error("Atom name not known: "+Aname);
1668       } else if(Rname=="CYS") {
1669         if(Aname=="BB") {
1670           atoi[i]=CYS_BB;
1671         } else if(Aname=="SC1") {
1672           atoi[i]=CYS_SC1;
1673         } else error("Atom name not known: "+Aname);
1674       } else if(Rname=="GLN") {
1675         if(Aname=="BB") {
1676           atoi[i]=GLN_BB;
1677         } else if(Aname=="SC1") {
1678           atoi[i]=GLN_SC1;
1679         } else error("Atom name not known: "+Aname);
1680       } else if(Rname=="GLU") {
1681         if(Aname=="BB") {
1682           atoi[i]=GLU_BB;
1683         } else if(Aname=="SC1") {
1684           atoi[i]=GLU_SC1;
1685         } else error("Atom name not known: "+Aname);
1686       } else if(Rname=="GLY") {
1687         if(Aname=="BB") {
1688           atoi[i]=GLY_BB;
1689         } else error("Atom name not known: "+Aname);
1690       } else if(Rname=="HIS") {
1691         if(Aname=="BB") {
1692           atoi[i]=HIS_BB;
1693         } else if(Aname=="SC1") {
1694           atoi[i]=HIS_SC1;
1695         } else if(Aname=="SC2") {
1696           atoi[i]=HIS_SC2;
1697         } else if(Aname=="SC3") {
1698           atoi[i]=HIS_SC3;
1699         } else error("Atom name not known: "+Aname);
1700       } else if(Rname=="ILE") {
1701         if(Aname=="BB") {
1702           atoi[i]=ILE_BB;
1703         } else if(Aname=="SC1") {
1704           atoi[i]=ILE_SC1;
1705         } else error("Atom name not known: "+Aname);
1706       } else if(Rname=="LEU") {
1707         if(Aname=="BB") {
1708           atoi[i]=LEU_BB;
1709         } else if(Aname=="SC1") {
1710           atoi[i]=LEU_SC1;
1711         } else error("Atom name not known: "+Aname);
1712       } else if(Rname=="LYS") {
1713         if(Aname=="BB") {
1714           atoi[i]=LYS_BB;
1715         } else if(Aname=="SC1") {
1716           atoi[i]=LYS_SC1;
1717         } else if(Aname=="SC2") {
1718           atoi[i]=LYS_SC2;
1719         } else error("Atom name not known: "+Aname);
1720       } else if(Rname=="MET") {
1721         if(Aname=="BB") {
1722           atoi[i]=MET_BB;
1723         } else if(Aname=="SC1") {
1724           atoi[i]=MET_SC1;
1725         } else error("Atom name not known: "+Aname);
1726       } else if(Rname=="PHE") {
1727         if(Aname=="BB") {
1728           atoi[i]=PHE_BB;
1729         } else if(Aname=="SC1") {
1730           atoi[i]=PHE_SC1;
1731         } else if(Aname=="SC2") {
1732           atoi[i]=PHE_SC2;
1733         } else if(Aname=="SC3") {
1734           atoi[i]=PHE_SC3;
1735         } else error("Atom name not known: "+Aname);
1736       } else if(Rname=="PRO") {
1737         if(Aname=="BB") {
1738           atoi[i]=PRO_BB;
1739         } else if(Aname=="SC1") {
1740           atoi[i]=PRO_SC1;
1741         } else error("Atom name not known: "+Aname);
1742       } else if(Rname=="SER") {
1743         if(Aname=="BB") {
1744           atoi[i]=SER_BB;
1745         } else if(Aname=="SC1") {
1746           atoi[i]=SER_SC1;
1747         } else error("Atom name not known: "+Aname);
1748       } else if(Rname=="THR") {
1749         if(Aname=="BB") {
1750           atoi[i]=THR_BB;
1751         } else if(Aname=="SC1") {
1752           atoi[i]=THR_SC1;
1753         } else error("Atom name not known: "+Aname);
1754       } else if(Rname=="TRP") {
1755         if(Aname=="BB") {
1756           atoi[i]=TRP_BB;
1757         } else if(Aname=="SC1") {
1758           atoi[i]=TRP_SC1;
1759         } else if(Aname=="SC2") {
1760           atoi[i]=TRP_SC2;
1761         } else if(Aname=="SC3") {
1762           atoi[i]=TRP_SC3;
1763         } else if(Aname=="SC4") {
1764           atoi[i]=TRP_SC4;
1765         } else error("Atom name not known: "+Aname);
1766       } else if(Rname=="TYR") {
1767         if(Aname=="BB") {
1768           atoi[i]=TYR_BB;
1769         } else if(Aname=="SC1") {
1770           atoi[i]=TYR_SC1;
1771         } else if(Aname=="SC2") {
1772           atoi[i]=TYR_SC2;
1773         } else if(Aname=="SC3") {
1774           atoi[i]=TYR_SC3;
1775         } else error("Atom name not known: "+Aname);
1776       } else if(Rname=="VAL") {
1777         if(Aname=="BB") {
1778           atoi[i]=VAL_BB;
1779         } else if(Aname=="SC1") {
1780           atoi[i]=VAL_SC1;
1781         } else error("Atom name not known: "+Aname);
1782       } else if(Rname=="  A") {
1783         if(Aname=="BB1") {
1784           atoi[i]=A_BB1;
1785         } else if(Aname=="BB2") {
1786           atoi[i]=A_BB2;
1787         } else if(Aname=="BB3") {
1788           atoi[i]=A_BB3;
1789         } else if(Aname=="SC1") {
1790           atoi[i]=A_SC1;
1791         } else if(Aname=="SC2") {
1792           atoi[i]=A_SC2;
1793         } else if(Aname=="SC3") {
1794           atoi[i]=A_SC3;
1795         } else if(Aname=="SC4") {
1796           atoi[i]=A_SC4;
1797         } else if(Aname=="3TE") {
1798           atoi[i]=A_3TE;
1799         } else if(Aname=="5TE") {
1800           atoi[i]=A_5TE;
1801         } else if(Aname=="TE3") {
1802           atoi[i]=A_TE3;
1803         } else if(Aname=="TE5") {
1804           atoi[i]=A_TE5;
1805         } else error("Atom name not known: "+Aname);
1806       } else if(Rname=="  C") {
1807         if(Aname=="BB1") {
1808           atoi[i]=C_BB1;
1809         } else if(Aname=="BB2") {
1810           atoi[i]=C_BB2;
1811         } else if(Aname=="BB3") {
1812           atoi[i]=C_BB3;
1813         } else if(Aname=="SC1") {
1814           atoi[i]=C_SC1;
1815         } else if(Aname=="SC2") {
1816           atoi[i]=C_SC2;
1817         } else if(Aname=="SC3") {
1818           atoi[i]=C_SC3;
1819         } else if(Aname=="3TE") {
1820           atoi[i]=C_3TE;
1821         } else if(Aname=="5TE") {
1822           atoi[i]=C_5TE;
1823         } else if(Aname=="TE3") {
1824           atoi[i]=C_TE3;
1825         } else if(Aname=="TE5") {
1826           atoi[i]=C_TE5;
1827         } else error("Atom name not known: "+Aname);
1828       } else if(Rname=="  G") {
1829         if(Aname=="BB1") {
1830           atoi[i]=G_BB1;
1831         } else if(Aname=="BB2") {
1832           atoi[i]=G_BB2;
1833         } else if(Aname=="BB3") {
1834           atoi[i]=G_BB3;
1835         } else if(Aname=="SC1") {
1836           atoi[i]=G_SC1;
1837         } else if(Aname=="SC2") {
1838           atoi[i]=G_SC2;
1839         } else if(Aname=="SC3") {
1840           atoi[i]=G_SC3;
1841         } else if(Aname=="SC4") {
1842           atoi[i]=G_SC4;
1843         } else if(Aname=="3TE") {
1844           atoi[i]=G_3TE;
1845         } else if(Aname=="5TE") {
1846           atoi[i]=G_5TE;
1847         } else if(Aname=="TE3") {
1848           atoi[i]=G_TE3;
1849         } else if(Aname=="TE5") {
1850           atoi[i]=G_TE5;
1851         } else error("Atom name not known: "+Aname);
1852       } else if(Rname=="  U") {
1853         if(Aname=="BB1") {
1854           atoi[i]=U_BB1;
1855         } else if(Aname=="BB2") {
1856           atoi[i]=U_BB2;
1857         } else if(Aname=="BB3") {
1858           atoi[i]=U_BB3;
1859         } else if(Aname=="SC1") {
1860           atoi[i]=U_SC1;
1861         } else if(Aname=="SC2") {
1862           atoi[i]=U_SC2;
1863         } else if(Aname=="SC3") {
1864           atoi[i]=U_SC3;
1865         } else if(Aname=="3TE") {
1866           atoi[i]=U_3TE;
1867         } else if(Aname=="5TE") {
1868           atoi[i]=U_5TE;
1869         } else if(Aname=="TE3") {
1870           atoi[i]=U_TE3;
1871         } else if(Aname=="TE5") {
1872           atoi[i]=U_TE5;
1873         } else error("Atom name not known: "+Aname);
1874       } else if(Rname==" DA") {
1875         if(Aname=="BB1") {
1876           atoi[i]=DA_BB1;
1877         } else if(Aname=="BB2") {
1878           atoi[i]=DA_BB2;
1879         } else if(Aname=="BB3") {
1880           atoi[i]=DA_BB3;
1881         } else if(Aname=="SC1") {
1882           atoi[i]=DA_SC1;
1883         } else if(Aname=="SC2") {
1884           atoi[i]=DA_SC2;
1885         } else if(Aname=="SC3") {
1886           atoi[i]=DA_SC3;
1887         } else if(Aname=="SC4") {
1888           atoi[i]=DA_SC4;
1889         } else if(Aname=="3TE") {
1890           atoi[i]=DA_3TE;
1891         } else if(Aname=="5TE") {
1892           atoi[i]=DA_5TE;
1893         } else if(Aname=="TE3") {
1894           atoi[i]=DA_TE3;
1895         } else if(Aname=="TE5") {
1896           atoi[i]=DA_TE5;
1897         } else error("Atom name not known: "+Aname);
1898       } else if(Rname==" DC") {
1899         if(Aname=="BB1") {
1900           atoi[i]=DC_BB1;
1901         } else if(Aname=="BB2") {
1902           atoi[i]=DC_BB2;
1903         } else if(Aname=="BB3") {
1904           atoi[i]=DC_BB3;
1905         } else if(Aname=="SC1") {
1906           atoi[i]=DC_SC1;
1907         } else if(Aname=="SC2") {
1908           atoi[i]=DC_SC2;
1909         } else if(Aname=="SC3") {
1910           atoi[i]=DC_SC3;
1911         } else if(Aname=="3TE") {
1912           atoi[i]=DC_3TE;
1913         } else if(Aname=="5TE") {
1914           atoi[i]=DC_5TE;
1915         } else if(Aname=="TE3") {
1916           atoi[i]=DC_TE3;
1917         } else if(Aname=="TE5") {
1918           atoi[i]=DC_TE5;
1919         } else error("Atom name not known: "+Aname);
1920       } else if(Rname==" DG") {
1921         if(Aname=="BB1") {
1922           atoi[i]=DG_BB1;
1923         } else if(Aname=="BB2") {
1924           atoi[i]=DG_BB2;
1925         } else if(Aname=="BB3") {
1926           atoi[i]=DG_BB3;
1927         } else if(Aname=="SC1") {
1928           atoi[i]=DG_SC1;
1929         } else if(Aname=="SC2") {
1930           atoi[i]=DG_SC2;
1931         } else if(Aname=="SC3") {
1932           atoi[i]=DG_SC3;
1933         } else if(Aname=="SC4") {
1934           atoi[i]=DG_SC4;
1935         } else if(Aname=="3TE") {
1936           atoi[i]=DG_3TE;
1937         } else if(Aname=="5TE") {
1938           atoi[i]=DG_5TE;
1939         } else if(Aname=="TE3") {
1940           atoi[i]=DG_TE3;
1941         } else if(Aname=="TE5") {
1942           atoi[i]=DG_TE5;
1943         } else error("Atom name not known: "+Aname);
1944       } else if(Rname==" DT") {
1945         if(Aname=="BB1") {
1946           atoi[i]=DT_BB1;
1947         } else if(Aname=="BB2") {
1948           atoi[i]=DT_BB2;
1949         } else if(Aname=="BB3") {
1950           atoi[i]=DT_BB3;
1951         } else if(Aname=="SC1") {
1952           atoi[i]=DT_SC1;
1953         } else if(Aname=="SC2") {
1954           atoi[i]=DT_SC2;
1955         } else if(Aname=="SC3") {
1956           atoi[i]=DT_SC3;
1957         } else if(Aname=="3TE") {
1958           atoi[i]=DT_3TE;
1959         } else if(Aname=="5TE") {
1960           atoi[i]=DT_5TE;
1961         } else if(Aname=="TE3") {
1962           atoi[i]=DT_TE3;
1963         } else if(Aname=="TE5") {
1964           atoi[i]=DT_TE5;
1965         } else error("Atom name not known: "+Aname);
1966       } else error("Residue not known: "+Rname);
1967     }
1968   } else {
1969     error("MOLINFO DATA not found\n");
1970   }
1971 }
1972 
calculateASF(const vector<AtomNumber> & atoms,vector<vector<long double>> & FF_tmp,const double rho)1973 double SAXS::calculateASF(const vector<AtomNumber> &atoms, vector<vector<long double> > &FF_tmp, const double rho)
1974 {
1975   map<string, unsigned> AA_map;
1976   AA_map["H"] = H;
1977   AA_map["C"] = C;
1978   AA_map["N"] = N;
1979   AA_map["O"] = O;
1980   AA_map["P"] = P;
1981   AA_map["S"] = S;
1982 
1983   vector<vector<double> > param_a;
1984   vector<vector<double> > param_b;
1985   vector<double> param_c;
1986   vector<double> param_v;
1987 
1988   param_a.resize(NTT, vector<double>(5));
1989   param_b.resize(NTT, vector<double>(5));
1990   param_c.resize(NTT);
1991   param_v.resize(NTT);
1992 
1993   param_a[H][0] = 0.493002; param_b[H][0] = 10.5109; param_c[H] = 0.003038;
1994   param_a[H][1] = 0.322912; param_b[H][1] = 26.1257; param_v[H] = 5.15;
1995   param_a[H][2] = 0.140191; param_b[H][2] = 3.14236;
1996   param_a[H][3] = 0.040810; param_b[H][3] = 57.7997;
1997   param_a[H][4] = 0.0;      param_b[H][4] = 1.0;
1998 
1999   param_a[C][0] = 2.31000; param_b[C][0] = 20.8439; param_c[C] = 0.215600;
2000   param_a[C][1] = 1.02000; param_b[C][1] = 10.2075; param_v[C] = 16.44;
2001   param_a[C][2] = 1.58860; param_b[C][2] = 0.56870;
2002   param_a[C][3] = 0.86500; param_b[C][3] = 51.6512;
2003   param_a[C][4] = 0.0;     param_b[C][4] = 1.0;
2004 
2005   param_a[N][0] = 12.2126; param_b[N][0] = 0.00570; param_c[N] = -11.529;
2006   param_a[N][1] = 3.13220; param_b[N][1] = 9.89330; param_v[N] = 2.49;
2007   param_a[N][2] = 2.01250; param_b[N][2] = 28.9975;
2008   param_a[N][3] = 1.16630; param_b[N][3] = 0.58260;
2009   param_a[N][4] = 0.0;     param_b[N][4] = 1.0;
2010 
2011   param_a[O][0] = 3.04850; param_b[O][0] = 13.2771; param_c[O] = 0.250800 ;
2012   param_a[O][1] = 2.28680; param_b[O][1] = 5.70110; param_v[O] = 9.13;
2013   param_a[O][2] = 1.54630; param_b[O][2] = 0.32390;
2014   param_a[O][3] = 0.86700; param_b[O][3] = 32.9089;
2015   param_a[O][4] = 0.0;     param_b[O][4] = 1.0;
2016 
2017   param_a[P][0] = 6.43450; param_b[P][0] = 1.90670; param_c[P] = 1.11490;
2018   param_a[P][1] = 4.17910; param_b[P][1] = 27.1570; param_v[P] = 5.73;
2019   param_a[P][2] = 1.78000; param_b[P][2] = 0.52600;
2020   param_a[P][3] = 1.49080; param_b[P][3] = 68.1645;
2021   param_a[P][4] = 0.0;     param_b[P][4] = 1.0;
2022 
2023   param_a[S][0] = 6.90530; param_b[S][0] = 1.46790; param_c[S] = 0.866900;
2024   param_a[S][1] = 5.20340; param_b[S][1] = 22.2151; param_v[S] = 19.86;
2025   param_a[S][2] = 1.43790; param_b[S][2] = 0.25360;
2026   param_a[S][3] = 1.58630; param_b[S][3] = 56.1720;
2027   param_a[S][4] = 0.0;     param_b[S][4] = 1.0;
2028 
2029   auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
2030 
2031   double Iq0=0.;
2032   if( moldat ) {
2033     log<<"  MOLINFO DATA found with label " <<moldat->getLabel()<<", using proper atom names\n";
2034     // cycle over the atom types
2035     for(unsigned i=0; i<NTT; i++) {
2036       const double volr = pow(param_v[i], (2.0/3.0)) /(4. * M_PI);
2037       // cycle over q
2038       for(unsigned k=0; k<q_list.size(); ++k) {
2039         const double q = q_list[k];
2040         const double s = q / (4. * M_PI);
2041         FF_tmp[k][i] = param_c[i];
2042         // SUM [a_i * EXP( - b_i * (q/4pi)^2 )] Waasmaier and Kirfel (1995)
2043         for(unsigned j=0; j<4; j++) {
2044           FF_tmp[k][i] += param_a[i][j]*exp(-param_b[i][j]*s*s);
2045         }
2046         // subtract solvation: rho * v_i * EXP( (- v_i^(2/3) / (4pi)) * q^2  ) // since  D in Fraser 1978 is 2*s
2047         FF_tmp[k][i] -= rho*param_v[i]*exp(-volr*q*q);
2048       }
2049     }
2050     // cycle over the atoms to assign the atom type and calculate I0
2051     for(unsigned i=0; i<atoms.size(); ++i) {
2052       // get atom name
2053       string name = moldat->getAtomName(atoms[i]);
2054       char type;
2055       // get atom type
2056       char first = name.at(0);
2057       // GOLDEN RULE: type is first letter, if not a number
2058       if (!isdigit(first)) {
2059         type = first;
2060         // otherwise is the second
2061       } else {
2062         type = name.at(1);
2063       }
2064       std::string type_s = std::string(1,type);
2065       if(AA_map.find(type_s) != AA_map.end()) {
2066         const unsigned index=AA_map[type_s];
2067         atoi[i] = AA_map[type_s];
2068         for(unsigned j=0; j<4; j++) Iq0 += param_a[index][j];
2069         Iq0 = Iq0 -rho*param_v[index] + param_c[index];
2070       } else {
2071         error("Wrong atom type "+type_s+" from atom name "+name+"\n");
2072       }
2073     }
2074   } else {
2075     error("MOLINFO DATA not found\n");
2076   }
2077 
2078   return Iq0;
2079 }
2080 
2081 }
2082 }
2083