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