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