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 and
24 : Carlo Camilloni
25 : */
26 :
27 : #include "MetainferenceBase.h"
28 : #include "core/ActionRegister.h"
29 : #include "core/ActionSet.h"
30 : #include "core/SetupMolInfo.h"
31 : #include "tools/Communicator.h"
32 : #include "tools/Pbc.h"
33 :
34 : #include <string>
35 : #include <cmath>
36 : #include <map>
37 :
38 : #ifndef M_PI
39 : #define M_PI 3.14159265358979323846
40 : #endif
41 :
42 : using namespace std;
43 :
44 : namespace PLMD {
45 : namespace isdb {
46 :
47 : //+PLUMEDOC ISDB_COLVAR SAXS
48 : /*
49 : Calculates SAXS scattered intensity using the Debye equation.
50 :
51 : Intensities are calculated for a set of scattering lenght set using QVALUES numbered keywords, QVALUE cannot be 0.
52 : Structure factors can be either assigned using a polynomial expansion to any order using the PARAMETERS keywords;
53 : automatically assigned to atoms using the ATOMISTIC flag reading a PDB file, a correction for the water density is automatically added;
54 : automatically assigned to Martini pseudoatoms usign the MARTINI flag.
55 : The calculated intensities can be scaled using the SCEXP keywords. This is applied by rescaling the structure factors.
56 : Experimental reference intensities can be added using the ADDEXP and EXPINT flag and keywords.
57 : \ref METAINFERENCE can be activated using DOSCORE and the other relevant keywords.
58 :
59 : \par Examples
60 : in the following example the saxs intensities for a martini model are calculated. structure factors
61 : are obtained from the pdb file indicated in the MOLINFO.
62 :
63 : \plumedfile
64 : MOLINFO STRUCTURE=template.pdb
65 :
66 : SAXS ...
67 : LABEL=saxs
68 : ATOMS=1-355
69 : ADDEXP
70 : SCEXP=3920000
71 : MARTINI
72 : QVALUE1=0.02 EXPINT1=1.0902
73 : QVALUE2=0.05 EXPINT2=0.790632
74 : QVALUE3=0.08 EXPINT3=0.453808
75 : QVALUE4=0.11 EXPINT4=0.254737
76 : QVALUE5=0.14 EXPINT5=0.154928
77 : QVALUE6=0.17 EXPINT6=0.0921503
78 : QVALUE7=0.2 EXPINT7=0.052633
79 : QVALUE8=0.23 EXPINT8=0.0276557
80 : QVALUE9=0.26 EXPINT9=0.0122775
81 : QVALUE10=0.29 EXPINT10=0.00880634
82 : QVALUE11=0.32 EXPINT11=0.0137301
83 : QVALUE12=0.35 EXPINT12=0.0180036
84 : QVALUE13=0.38 EXPINT13=0.0193374
85 : QVALUE14=0.41 EXPINT14=0.0210131
86 : QVALUE15=0.44 EXPINT15=0.0220506
87 : ... SAXS
88 :
89 : PRINT ARG=(saxs\.q_.*),(saxs\.exp_.*) FILE=colvar STRIDE=1
90 :
91 : \endplumedfile
92 :
93 : */
94 : //+ENDPLUMEDOC
95 :
96 24 : class SAXS :
97 : public MetainferenceBase
98 : {
99 : private:
100 : bool pbc;
101 : bool serial;
102 : vector<double> q_list;
103 : vector<double> FF_rank;
104 : vector<vector<double> > FF_value;
105 :
106 : void getMartiniSFparam(const vector<AtomNumber> &atoms, vector<vector<long double> > ¶meter);
107 : void calculateASF(const vector<AtomNumber> &atoms, vector<vector<long double> > &FF_tmp, const double rho);
108 :
109 : public:
110 : static void registerKeywords( Keywords& keys );
111 : explicit SAXS(const ActionOptions&);
112 : virtual void calculate();
113 : void update();
114 : };
115 :
116 6460 : PLUMED_REGISTER_ACTION(SAXS,"SAXS")
117 :
118 9 : void SAXS::registerKeywords(Keywords& keys) {
119 9 : componentsAreNotOptional(keys);
120 9 : useCustomisableComponents(keys);
121 9 : MetainferenceBase::registerKeywords(keys);
122 27 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
123 27 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
124 27 : keys.addFlag("ATOMISTIC",false,"calculate SAXS for an atomistic model");
125 27 : keys.addFlag("MARTINI",false,"calculate SAXS for a Martini model");
126 36 : keys.add("atoms","ATOMS","The atoms to be included in the calculation, e.g. the whole protein.");
127 36 : keys.add("numbered","QVALUE","Selected scattering lenghts in Angstrom are given as QVALUE1, QVALUE2, ... .");
128 36 : keys.add("numbered","PARAMETERS","Used parameter Keywords like PARAMETERS1, PARAMETERS2. These are used to calculate the structure factor for the i-th atom/bead.");
129 45 : keys.add("compulsory","WATERDENS","0.334","Density of the water to be used for the correction of atomistic structure factors.");
130 27 : keys.addFlag("ADDEXP",false,"Set to TRUE if you want to have fixed components with the experimental values.");
131 36 : keys.add("numbered","EXPINT","Add an experimental value for each q value.");
132 45 : keys.add("compulsory","SCEXP","1.0","SCALING value of the experimental data. Usefull to simplify the comparison.");
133 36 : keys.addOutputComponent("q","default","the # SAXS of q");
134 36 : keys.addOutputComponent("exp","ADDEXP","the # experimental intensity");
135 9 : }
136 :
137 8 : SAXS::SAXS(const ActionOptions&ao):
138 : PLUMED_METAINF_INIT(ao),
139 : pbc(true),
140 24 : serial(false)
141 : {
142 : vector<AtomNumber> atoms;
143 16 : parseAtomList("ATOMS",atoms);
144 8 : const unsigned size = atoms.size();
145 :
146 16 : parseFlag("SERIAL",serial);
147 :
148 8 : bool nopbc=!pbc;
149 16 : parseFlag("NOPBC",nopbc);
150 8 : pbc=!nopbc;
151 :
152 8 : double scexp = 0;
153 16 : parse("SCEXP",scexp);
154 8 : if(scexp==0) scexp=1.0;
155 :
156 : unsigned ntarget=0;
157 : for(unsigned i=0;; ++i) {
158 : double t_list;
159 256 : if( !parseNumbered( "QVALUE", i+1, t_list) ) break;
160 120 : q_list.push_back(t_list);
161 : ntarget++;
162 120 : }
163 : const unsigned numq = ntarget;
164 :
165 :
166 8 : bool atomistic=false;
167 16 : parseFlag("ATOMISTIC",atomistic);
168 8 : bool martini=false;
169 16 : parseFlag("MARTINI",martini);
170 :
171 8 : if(martini&&atomistic) error("You cannot use martini and atomistic at the same time");
172 :
173 8 : double rho = 0.334;
174 16 : parse("WATERDENS", rho);
175 :
176 8 : vector<vector<long double> > FF_tmp;
177 16 : FF_tmp.resize(numq,vector<long double>(size));
178 8 : if(!atomistic&&!martini) {
179 : //read in parameter vector
180 0 : vector<vector<long double> > parameter;
181 0 : parameter.resize(size);
182 : ntarget=0;
183 0 : for(unsigned i=0; i<size; ++i) {
184 0 : if( !parseNumberedVector( "PARAMETERS", i+1, parameter[i]) ) break;
185 : ntarget++;
186 : }
187 0 : if( ntarget!=size ) error("found wrong number of parameter vectors");
188 0 : for(unsigned i=0; i<size; ++i) {
189 0 : for(unsigned k=0; k<numq; ++k) {
190 0 : for(unsigned j=0; j<parameter[i].size(); ++j) {
191 0 : FF_tmp[k][i]+= parameter[i][j]*powl(static_cast<long double>(q_list[k]),j);
192 : }
193 : }
194 0 : }
195 8 : } else if(martini) {
196 : //read in parameter vector
197 8 : vector<vector<long double> > parameter;
198 8 : parameter.resize(size);
199 8 : getMartiniSFparam(atoms, parameter);
200 5688 : for(unsigned i=0; i<size; ++i) {
201 88040 : for(unsigned k=0; k<numq; ++k) {
202 979800 : for(unsigned j=0; j<parameter[i].size(); ++j) {
203 894600 : FF_tmp[k][i]+= parameter[i][j]*powl(static_cast<long double>(q_list[k]),j);
204 : }
205 : }
206 : }
207 0 : } else if(atomistic) {
208 0 : calculateASF(atoms, FF_tmp, rho);
209 : }
210 :
211 : // Calculate Rank of FF_matrix
212 8 : FF_rank.resize(numq);
213 16 : FF_value.resize(numq,vector<double>(size));
214 248 : for(unsigned k=0; k<numq; ++k) {
215 85320 : for(unsigned i=0; i<size; i++) {
216 127800 : FF_value[k][i] = static_cast<double>(FF_tmp[k][i])/sqrt(scexp);
217 85200 : FF_rank[k]+=FF_value[k][i]*FF_value[k][i];
218 : }
219 : }
220 :
221 8 : bool exp=false;
222 16 : parseFlag("ADDEXP",exp);
223 8 : if(getDoScore()) exp=true;
224 :
225 : vector<double> expint;
226 8 : expint.resize( numq );
227 : ntarget=0;
228 128 : for(unsigned i=0; i<numq; ++i) {
229 360 : if( !parseNumbered( "EXPINT", i+1, expint[i] ) ) break;
230 : ntarget++;
231 : }
232 8 : if( ntarget!=numq && exp==true) error("found wrong number of EXPINT values");
233 :
234 8 : if(pbc) log.printf(" using periodic boundary conditions\n");
235 0 : else log.printf(" without periodic boundary conditions\n");
236 248 : for(unsigned i=0; i<numq; i++) {
237 240 : if(q_list[i]==0.) error("it is not possible to set q=0\n");
238 240 : log.printf(" my q: %lf \n",q_list[i]);
239 : }
240 :
241 8 : if(!getDoScore()) {
242 124 : for(unsigned i=0; i<numq; i++) {
243 60 : std::string num; Tools::convert(i,num);
244 120 : addComponentWithDerivatives("q_"+num);
245 120 : componentIsNotPeriodic("q_"+num);
246 : }
247 4 : if(exp) {
248 124 : for(unsigned i=0; i<numq; i++) {
249 60 : std::string num; Tools::convert(i,num);
250 120 : addComponent("exp_"+num);
251 120 : componentIsNotPeriodic("exp_"+num);
252 120 : Value* comp=getPntrToComponent("exp_"+num);
253 120 : comp->set(expint[i]);
254 : }
255 : }
256 : } else {
257 124 : for(unsigned i=0; i<numq; i++) {
258 60 : std::string num; Tools::convert(i,num);
259 120 : addComponent("q_"+num);
260 120 : componentIsNotPeriodic("q_"+num);
261 : }
262 124 : for(unsigned i=0; i<numq; i++) {
263 60 : std::string num; Tools::convert(i,num);
264 120 : addComponent("exp_"+num);
265 120 : componentIsNotPeriodic("exp_"+num);
266 120 : Value* comp=getPntrToComponent("exp_"+num);
267 120 : comp->set(expint[i]);
268 : }
269 : }
270 :
271 : // convert units to nm^-1
272 248 : for(unsigned i=0; i<numq; ++i) {
273 240 : q_list[i]=q_list[i]*10.0; //factor 10 to convert from A^-1 to nm^-1
274 : }
275 8 : log<<" Bibliography ";
276 24 : log<<plumed.cite("Jussupow, et al. (in preparation)");
277 24 : if(martini) log<<plumed.cite("Niebling, Björling, Westenhoff, J Appl Crystallogr 47, 1190–1198 (2014).");
278 8 : if(atomistic) {
279 0 : log<<plumed.cite("Fraser, MacRae, Suzuki, J. Appl. Crystallogr., 11, 693–694 (1978).");
280 0 : log<<plumed.cite("Brown, Fox, Maslen, O'Keefe, Willis, International Tables for Crystallography C, 554–595 (International Union of Crystallography, 2006).");
281 : }
282 24 : log<< plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)");
283 8 : log<<"\n";
284 :
285 8 : requestAtoms(atoms);
286 8 : if(getDoScore()) {
287 4 : setParameters(expint);
288 4 : Initialise(numq);
289 : }
290 8 : setDerivatives();
291 8 : checkRead();
292 8 : }
293 :
294 88 : void SAXS::calculate() {
295 88 : if(pbc) makeWhole();
296 :
297 : const unsigned size = getNumberOfAtoms();
298 88 : const unsigned numq = q_list.size();
299 :
300 88 : unsigned stride = comm.Get_size();
301 88 : unsigned rank = comm.Get_rank();
302 88 : if(serial) {
303 : stride = 1;
304 : rank = 0;
305 : }
306 :
307 88 : vector<Vector> deriv(numq*size);
308 88 : vector<double> sum(numq,0);
309 :
310 264 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
311 : for (unsigned k=0; k<numq; k++) {
312 1320 : const unsigned kdx=k*size;
313 457980 : for (unsigned i=rank; i<size-1; i+=stride) {
314 456660 : const double FF=2.*FF_value[k][i];
315 456660 : const Vector posi=getPosition(i);
316 228330 : Vector dsum;
317 40756905 : for (unsigned j=i+1; j<size ; j++) {
318 81057150 : const Vector c_distances = delta(posi,getPosition(j));
319 40528575 : const double m_distances = c_distances.modulo();
320 40528575 : const double qdist = q_list[k]*m_distances;
321 81057150 : const double FFF = FF*FF_value[k][j];
322 40528575 : const double tsq = FFF*sin(qdist)/qdist;
323 40528575 : const double tcq = FFF*cos(qdist);
324 40528575 : const double tmp = (tcq-tsq)/(m_distances*m_distances);
325 40528575 : const Vector dd = c_distances*tmp;
326 40528575 : dsum += dd;
327 81057150 : deriv[kdx+j] += dd;
328 81057150 : sum[k] += tsq;
329 : }
330 456660 : deriv[kdx+i] -= dsum;
331 : }
332 : }
333 :
334 88 : if(!serial) {
335 176 : comm.Sum(&deriv[0][0], 3*deriv.size());
336 176 : comm.Sum(&sum[0], numq);
337 : }
338 :
339 2728 : for (unsigned k=0; k<numq; k++) {
340 3960 : sum[k]+=FF_rank[k];
341 1320 : string num; Tools::convert(k,num);
342 2640 : Value* val=getPntrToComponent("q_"+num);
343 1320 : val->set(sum[k]);
344 2580 : if(getDoScore()) setCalcData(k, sum[k]);
345 : }
346 :
347 88 : if(getDoScore()) {
348 : /* Metainference */
349 84 : double score = getScore();
350 : setScore(score);
351 : }
352 :
353 2728 : for (unsigned k=0; k<numq; k++) {
354 1320 : const unsigned kdx=k*size;
355 1320 : Tensor deriv_box;
356 : Value* val;
357 1320 : if(!getDoScore()) {
358 60 : string num; Tools::convert(k,num);
359 120 : val=getPntrToComponent("q_"+num);
360 42660 : for(unsigned i=0; i<size; i++) {
361 42600 : setAtomsDerivatives(val, i, deriv[kdx+i]);
362 42600 : deriv_box += Tensor(getPosition(i),deriv[kdx+i]);
363 : }
364 : } else {
365 2520 : val=getPntrToComponent("score");
366 895860 : for(unsigned i=0; i<size; i++) {
367 1341900 : setAtomsDerivatives(val, i, deriv[kdx+i]*getMetaDer(k));
368 894600 : deriv_box += Tensor(getPosition(i),deriv[kdx+i]*getMetaDer(k));
369 : }
370 : }
371 1320 : setBoxDerivatives(val, -deriv_box);
372 : }
373 88 : }
374 :
375 88 : void SAXS::update() {
376 : // write status file
377 176 : if(getWstride()>0&& (getStep()%getWstride()==0 || getCPT()) ) writeStatus();
378 88 : }
379 :
380 8 : void SAXS::getMartiniSFparam(const vector<AtomNumber> &atoms, vector<vector<long double> > ¶meter)
381 : {
382 16 : vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
383 8 : if( moldat.size()==1 ) {
384 8 : log<<" MOLINFO DATA found, using proper atom names\n";
385 8536 : for(unsigned i=0; i<atoms.size(); ++i) {
386 5680 : string Aname = moldat[0]->getAtomName(atoms[i]);
387 5680 : string Rname = moldat[0]->getResidueName(atoms[i]);
388 2840 : if(Rname=="ALA") {
389 48 : if(Aname=="BB") {
390 96 : parameter[i].push_back(9.045);
391 96 : parameter[i].push_back(-0.098114);
392 96 : parameter[i].push_back(7.54281);
393 96 : parameter[i].push_back(-1.97438);
394 96 : parameter[i].push_back(-8.32689);
395 96 : parameter[i].push_back(6.09318);
396 96 : parameter[i].push_back(-1.18913);
397 0 : } else error("Atom name not known");
398 2792 : } else if(Rname=="ARG") {
399 192 : if(Aname=="BB") {
400 128 : parameter[i].push_back(10.729);
401 128 : parameter[i].push_back(-0.0392574);
402 128 : parameter[i].push_back(1.15382);
403 128 : parameter[i].push_back(-0.155999);
404 128 : parameter[i].push_back(-2.43619);
405 128 : parameter[i].push_back(1.72922);
406 128 : parameter[i].push_back(-0.33799);
407 128 : } else if(Aname=="SC1") {
408 128 : parameter[i].push_back(-2.796);
409 128 : parameter[i].push_back(0.472403);
410 128 : parameter[i].push_back(8.07424);
411 128 : parameter[i].push_back(4.37299);
412 128 : parameter[i].push_back(-10.7398);
413 128 : parameter[i].push_back(4.95677);
414 128 : parameter[i].push_back(-0.725797);
415 64 : } else if(Aname=="SC2") {
416 128 : parameter[i].push_back(15.396);
417 128 : parameter[i].push_back(0.0636736);
418 128 : parameter[i].push_back(-1.258);
419 128 : parameter[i].push_back(1.93135);
420 128 : parameter[i].push_back(-4.45031);
421 128 : parameter[i].push_back(2.49356);
422 128 : parameter[i].push_back(-0.410721);
423 0 : } else error("Atom name not known");
424 2600 : } else if(Rname=="ASN") {
425 64 : if(Aname=="BB") {
426 64 : parameter[i].push_back(10.738);
427 64 : parameter[i].push_back(-0.0402162);
428 64 : parameter[i].push_back(1.03007);
429 64 : parameter[i].push_back(-0.254174);
430 64 : parameter[i].push_back(-2.12015);
431 64 : parameter[i].push_back(1.55535);
432 64 : parameter[i].push_back(-0.30963);
433 32 : } else if(Aname=="SC1") {
434 64 : parameter[i].push_back(9.249);
435 64 : parameter[i].push_back(-0.0148678);
436 64 : parameter[i].push_back(5.52169);
437 64 : parameter[i].push_back(0.00853212);
438 64 : parameter[i].push_back(-6.71992);
439 64 : parameter[i].push_back(3.93622);
440 64 : parameter[i].push_back(-0.64973);
441 0 : } else error("Atom name not known");
442 2536 : } else if(Rname=="ASP") {
443 160 : if(Aname=="BB") {
444 160 : parameter[i].push_back(10.695);
445 160 : parameter[i].push_back(-0.0410247);
446 160 : parameter[i].push_back(1.03656);
447 160 : parameter[i].push_back(-0.298558);
448 160 : parameter[i].push_back(-2.06064);
449 160 : parameter[i].push_back(1.53495);
450 160 : parameter[i].push_back(-0.308365);
451 80 : } else if(Aname=="SC1") {
452 160 : parameter[i].push_back(9.476);
453 160 : parameter[i].push_back(-0.0254664);
454 160 : parameter[i].push_back(5.57899);
455 160 : parameter[i].push_back(-0.395027);
456 160 : parameter[i].push_back(-5.9407);
457 160 : parameter[i].push_back(3.48836);
458 160 : parameter[i].push_back(-0.569402);
459 0 : } else error("Atom name not known");
460 2376 : } else if(Rname=="CYS") {
461 0 : if(Aname=="BB") {
462 0 : parameter[i].push_back(10.698);
463 0 : parameter[i].push_back(-0.0233493);
464 0 : parameter[i].push_back(1.18257);
465 0 : parameter[i].push_back(0.0684464);
466 0 : parameter[i].push_back(-2.792);
467 0 : parameter[i].push_back(1.88995);
468 0 : parameter[i].push_back(-0.360229);
469 0 : } else if(Aname=="SC1") {
470 0 : parameter[i].push_back(8.199);
471 0 : parameter[i].push_back(-0.0261569);
472 0 : parameter[i].push_back(6.79677);
473 0 : parameter[i].push_back(-0.343845);
474 0 : parameter[i].push_back(-5.03578);
475 0 : parameter[i].push_back(2.7076);
476 0 : parameter[i].push_back(-0.420714);
477 0 : } else error("Atom name not known");
478 2376 : } else if(Rname=="GLN") {
479 192 : if(Aname=="BB") {
480 192 : parameter[i].push_back(10.728);
481 192 : parameter[i].push_back(-0.0391984);
482 192 : parameter[i].push_back(1.09264);
483 192 : parameter[i].push_back(-0.261555);
484 192 : parameter[i].push_back(-2.21245);
485 192 : parameter[i].push_back(1.62071);
486 192 : parameter[i].push_back(-0.322325);
487 96 : } else if(Aname=="SC1") {
488 192 : parameter[i].push_back(8.317);
489 192 : parameter[i].push_back(-0.229045);
490 192 : parameter[i].push_back(12.6338);
491 192 : parameter[i].push_back(-7.6719);
492 192 : parameter[i].push_back(-5.8376);
493 192 : parameter[i].push_back(5.53784);
494 192 : parameter[i].push_back(-1.12604);
495 0 : } else error("Atom name not known");
496 2184 : } else if(Rname=="GLU") {
497 192 : if(Aname=="BB") {
498 192 : parameter[i].push_back(10.694);
499 192 : parameter[i].push_back(-0.0521961);
500 192 : parameter[i].push_back(1.11153);
501 192 : parameter[i].push_back(-0.491995);
502 192 : parameter[i].push_back(-1.86236);
503 192 : parameter[i].push_back(1.45332);
504 192 : parameter[i].push_back(-0.29708);
505 96 : } else if(Aname=="SC1") {
506 192 : parameter[i].push_back(8.544);
507 192 : parameter[i].push_back(-0.249555);
508 192 : parameter[i].push_back(12.8031);
509 192 : parameter[i].push_back(-8.42696);
510 192 : parameter[i].push_back(-4.66486);
511 192 : parameter[i].push_back(4.90004);
512 192 : parameter[i].push_back(-1.01204);
513 0 : } else error("Atom name not known");
514 1992 : } else if(Rname=="GLY") {
515 104 : if(Aname=="BB") {
516 208 : parameter[i].push_back(9.977);
517 208 : parameter[i].push_back(-0.0285799);
518 208 : parameter[i].push_back(1.84236);
519 208 : parameter[i].push_back(-0.0315192);
520 208 : parameter[i].push_back(-2.88326);
521 208 : parameter[i].push_back(1.87323);
522 208 : parameter[i].push_back(-0.345773);
523 0 : } else error("Atom name not known");
524 1888 : } else if(Rname=="HIS") {
525 256 : if(Aname=="BB") {
526 128 : parameter[i].push_back(10.721);
527 128 : parameter[i].push_back(-0.0379337);
528 128 : parameter[i].push_back(1.06028);
529 128 : parameter[i].push_back(-0.236143);
530 128 : parameter[i].push_back(-2.17819);
531 128 : parameter[i].push_back(1.58357);
532 128 : parameter[i].push_back(-0.31345);
533 192 : } else if(Aname=="SC1") {
534 128 : parameter[i].push_back(-0.424);
535 128 : parameter[i].push_back(0.665176);
536 128 : parameter[i].push_back(3.4369);
537 128 : parameter[i].push_back(2.93795);
538 128 : parameter[i].push_back(-5.18288);
539 128 : parameter[i].push_back(2.12381);
540 128 : parameter[i].push_back(-0.284224);
541 128 : } else if(Aname=="SC2") {
542 128 : parameter[i].push_back(5.363);
543 128 : parameter[i].push_back(-0.0176945);
544 128 : parameter[i].push_back(2.9506);
545 128 : parameter[i].push_back(-0.387018);
546 128 : parameter[i].push_back(-1.83951);
547 128 : parameter[i].push_back(0.9703);
548 128 : parameter[i].push_back(-0.1458);
549 64 : } else if(Aname=="SC3") {
550 128 : parameter[i].push_back(5.784);
551 128 : parameter[i].push_back(-0.0293129);
552 128 : parameter[i].push_back(2.74167);
553 128 : parameter[i].push_back(-0.520875);
554 128 : parameter[i].push_back(-1.62949);
555 128 : parameter[i].push_back(0.902379);
556 128 : parameter[i].push_back(-0.139957);
557 0 : } else error("Atom name not known");
558 1632 : } else if(Rname=="ILE") {
559 224 : if(Aname=="BB") {
560 224 : parameter[i].push_back(10.699);
561 224 : parameter[i].push_back(-0.0188962);
562 224 : parameter[i].push_back(1.217);
563 224 : parameter[i].push_back(0.242481);
564 224 : parameter[i].push_back(-3.13898);
565 224 : parameter[i].push_back(2.07916);
566 224 : parameter[i].push_back(-0.392574);
567 112 : } else if(Aname=="SC1") {
568 224 : parameter[i].push_back(-4.448);
569 224 : parameter[i].push_back(1.20996);
570 224 : parameter[i].push_back(11.5141);
571 224 : parameter[i].push_back(6.98895);
572 224 : parameter[i].push_back(-19.1948);
573 224 : parameter[i].push_back(9.89207);
574 224 : parameter[i].push_back(-1.60877);
575 0 : } else error("Atom name not known");
576 1408 : } else if(Rname=="LEU") {
577 288 : if(Aname=="BB") {
578 288 : parameter[i].push_back(10.692);
579 288 : parameter[i].push_back(-0.0414917);
580 288 : parameter[i].push_back(1.1077);
581 288 : parameter[i].push_back(-0.288062);
582 288 : parameter[i].push_back(-2.17187);
583 288 : parameter[i].push_back(1.59879);
584 288 : parameter[i].push_back(-0.318545);
585 144 : } else if(Aname=="SC1") {
586 288 : parameter[i].push_back(-4.448);
587 288 : parameter[i].push_back(2.1063);
588 288 : parameter[i].push_back(6.72381);
589 288 : parameter[i].push_back(14.6954);
590 288 : parameter[i].push_back(-23.7197);
591 288 : parameter[i].push_back(10.7247);
592 288 : parameter[i].push_back(-1.59146);
593 0 : } else error("Atom name not known");
594 1120 : } else if(Rname=="LYS") {
595 336 : if(Aname=="BB") {
596 224 : parameter[i].push_back(10.706);
597 224 : parameter[i].push_back(-0.0468629);
598 224 : parameter[i].push_back(1.09477);
599 224 : parameter[i].push_back(-0.432751);
600 224 : parameter[i].push_back(-1.94335);
601 224 : parameter[i].push_back(1.49109);
602 224 : parameter[i].push_back(-0.302589);
603 224 : } else if(Aname=="SC1") {
604 224 : parameter[i].push_back(-2.796);
605 224 : parameter[i].push_back(0.508044);
606 224 : parameter[i].push_back(7.91436);
607 224 : parameter[i].push_back(4.54097);
608 224 : parameter[i].push_back(-10.8051);
609 224 : parameter[i].push_back(4.96204);
610 224 : parameter[i].push_back(-0.724414);
611 112 : } else if(Aname=="SC2") {
612 224 : parameter[i].push_back(3.070);
613 224 : parameter[i].push_back(-0.0101448);
614 224 : parameter[i].push_back(4.67994);
615 224 : parameter[i].push_back(-0.792529);
616 224 : parameter[i].push_back(-2.09142);
617 224 : parameter[i].push_back(1.02933);
618 224 : parameter[i].push_back(-0.137787);
619 0 : } else error("Atom name not known");
620 784 : } else if(Rname=="MET") {
621 32 : if(Aname=="BB") {
622 32 : parameter[i].push_back(10.671);
623 32 : parameter[i].push_back(-0.0433724);
624 32 : parameter[i].push_back(1.13784);
625 32 : parameter[i].push_back(-0.40768);
626 32 : parameter[i].push_back(-2.00555);
627 32 : parameter[i].push_back(1.51673);
628 32 : parameter[i].push_back(-0.305547);
629 16 : } else if(Aname=="SC1") {
630 32 : parameter[i].push_back(5.85);
631 32 : parameter[i].push_back(-0.0485798);
632 32 : parameter[i].push_back(17.0391);
633 32 : parameter[i].push_back(-3.65327);
634 32 : parameter[i].push_back(-13.174);
635 32 : parameter[i].push_back(8.68286);
636 32 : parameter[i].push_back(-1.56095);
637 0 : } else error("Atom name not known");
638 752 : } else if(Rname=="PHE") {
639 128 : if(Aname=="BB") {
640 64 : parameter[i].push_back(10.741);
641 64 : parameter[i].push_back(-0.0317275);
642 64 : parameter[i].push_back(1.15599);
643 64 : parameter[i].push_back(0.0276187);
644 64 : parameter[i].push_back(-2.74757);
645 64 : parameter[i].push_back(1.88783);
646 64 : parameter[i].push_back(-0.363525);
647 96 : } else if(Aname=="SC1") {
648 64 : parameter[i].push_back(-0.636);
649 64 : parameter[i].push_back(0.527882);
650 64 : parameter[i].push_back(6.77612);
651 64 : parameter[i].push_back(3.18508);
652 64 : parameter[i].push_back(-8.92826);
653 64 : parameter[i].push_back(4.29752);
654 64 : parameter[i].push_back(-0.65187);
655 64 : } else if(Aname=="SC2") {
656 64 : parameter[i].push_back(-0.424);
657 64 : parameter[i].push_back(0.389174);
658 64 : parameter[i].push_back(4.11761);
659 64 : parameter[i].push_back(2.29527);
660 64 : parameter[i].push_back(-4.7652);
661 64 : parameter[i].push_back(1.97023);
662 64 : parameter[i].push_back(-0.262318);
663 32 : } else if(Aname=="SC3") {
664 64 : parameter[i].push_back(-0.424);
665 64 : parameter[i].push_back(0.38927);
666 64 : parameter[i].push_back(4.11708);
667 64 : parameter[i].push_back(2.29623);
668 64 : parameter[i].push_back(-4.76592);
669 64 : parameter[i].push_back(1.97055);
670 64 : parameter[i].push_back(-0.262381);
671 0 : } else error("Atom name not known");
672 624 : } else if(Rname=="PRO") {
673 96 : if(Aname=="BB") {
674 96 : parameter[i].push_back(11.434);
675 96 : parameter[i].push_back(-0.033323);
676 96 : parameter[i].push_back(0.472014);
677 96 : parameter[i].push_back(-0.290854);
678 96 : parameter[i].push_back(-1.81409);
679 96 : parameter[i].push_back(1.39751);
680 96 : parameter[i].push_back(-0.280407);
681 48 : } else if(Aname=="SC1") {
682 96 : parameter[i].push_back(-2.796);
683 96 : parameter[i].push_back(0.95668);
684 96 : parameter[i].push_back(6.84197);
685 96 : parameter[i].push_back(6.43774);
686 96 : parameter[i].push_back(-12.5068);
687 96 : parameter[i].push_back(5.64597);
688 96 : parameter[i].push_back(-0.825206);
689 0 : } else error("Atom name not known");
690 528 : } else if(Rname=="SER") {
691 112 : if(Aname=="BB") {
692 112 : parameter[i].push_back(10.699);
693 112 : parameter[i].push_back(-0.0325828);
694 112 : parameter[i].push_back(1.20329);
695 112 : parameter[i].push_back(-0.0674351);
696 112 : parameter[i].push_back(-2.60749);
697 112 : parameter[i].push_back(1.80318);
698 112 : parameter[i].push_back(-0.346803);
699 56 : } else if(Aname=="SC1") {
700 112 : parameter[i].push_back(3.298);
701 112 : parameter[i].push_back(-0.0366801);
702 112 : parameter[i].push_back(5.11077);
703 112 : parameter[i].push_back(-1.46774);
704 112 : parameter[i].push_back(-1.48421);
705 112 : parameter[i].push_back(0.800326);
706 112 : parameter[i].push_back(-0.108314);
707 0 : } else error("Atom name not known");
708 416 : } else if(Rname=="THR") {
709 224 : if(Aname=="BB") {
710 224 : parameter[i].push_back(10.697);
711 224 : parameter[i].push_back(-0.0242955);
712 224 : parameter[i].push_back(1.24671);
713 224 : parameter[i].push_back(0.146423);
714 224 : parameter[i].push_back(-2.97429);
715 224 : parameter[i].push_back(1.97513);
716 224 : parameter[i].push_back(-0.371479);
717 112 : } else if(Aname=="SC1") {
718 224 : parameter[i].push_back(2.366);
719 224 : parameter[i].push_back(0.0297604);
720 224 : parameter[i].push_back(11.9216);
721 224 : parameter[i].push_back(-9.32503);
722 224 : parameter[i].push_back(1.9396);
723 224 : parameter[i].push_back(0.0804861);
724 224 : parameter[i].push_back(-0.0302721);
725 0 : } else error("Atom name not known");
726 192 : } else if(Rname=="TRP") {
727 0 : if(Aname=="BB") {
728 0 : parameter[i].push_back(10.689);
729 0 : parameter[i].push_back(-0.0265879);
730 0 : parameter[i].push_back(1.17819);
731 0 : parameter[i].push_back(0.0386457);
732 0 : parameter[i].push_back(-2.75634);
733 0 : parameter[i].push_back(1.88065);
734 0 : parameter[i].push_back(-0.360217);
735 0 : } else if(Aname=="SC1") {
736 0 : parameter[i].push_back(0.084);
737 0 : parameter[i].push_back(0.752407);
738 0 : parameter[i].push_back(5.3802);
739 0 : parameter[i].push_back(4.09281);
740 0 : parameter[i].push_back(-9.28029);
741 0 : parameter[i].push_back(4.45923);
742 0 : parameter[i].push_back(-0.689008);
743 0 : } else if(Aname=="SC2") {
744 0 : parameter[i].push_back(5.739);
745 0 : parameter[i].push_back(0.0298492);
746 0 : parameter[i].push_back(4.60446);
747 0 : parameter[i].push_back(1.34463);
748 0 : parameter[i].push_back(-5.69968);
749 0 : parameter[i].push_back(2.84924);
750 0 : parameter[i].push_back(-0.433781);
751 0 : } else if(Aname=="SC3") {
752 0 : parameter[i].push_back(-0.424);
753 0 : parameter[i].push_back(0.388576);
754 0 : parameter[i].push_back(4.11859);
755 0 : parameter[i].push_back(2.29485);
756 0 : parameter[i].push_back(-4.76255);
757 0 : parameter[i].push_back(1.96849);
758 0 : parameter[i].push_back(-0.262015);
759 0 : } else if(Aname=="SC4") {
760 0 : parameter[i].push_back(-0.424);
761 0 : parameter[i].push_back(0.387685);
762 0 : parameter[i].push_back(4.12153);
763 0 : parameter[i].push_back(2.29144);
764 0 : parameter[i].push_back(-4.7589);
765 0 : parameter[i].push_back(1.96686);
766 0 : parameter[i].push_back(-0.261786);
767 0 : } else error("Atom name not known");
768 192 : } else if(Rname=="TYR") {
769 64 : if(Aname=="BB") {
770 32 : parameter[i].push_back(10.689);
771 32 : parameter[i].push_back(-0.0193526);
772 32 : parameter[i].push_back(1.18241);
773 32 : parameter[i].push_back(0.207318);
774 32 : parameter[i].push_back(-3.0041);
775 32 : parameter[i].push_back(1.99335);
776 32 : parameter[i].push_back(-0.376482);
777 48 : } else if(Aname=="SC1") {
778 32 : parameter[i].push_back(-0.636);
779 32 : parameter[i].push_back(0.528902);
780 32 : parameter[i].push_back(6.78168);
781 32 : parameter[i].push_back(3.17769);
782 32 : parameter[i].push_back(-8.93667);
783 32 : parameter[i].push_back(4.30692);
784 32 : parameter[i].push_back(-0.653993);
785 32 : } else if(Aname=="SC2") {
786 32 : parameter[i].push_back(-0.424);
787 32 : parameter[i].push_back(0.388811);
788 32 : parameter[i].push_back(4.11851);
789 32 : parameter[i].push_back(2.29545);
790 32 : parameter[i].push_back(-4.7668);
791 32 : parameter[i].push_back(1.97131);
792 32 : parameter[i].push_back(-0.262534);
793 16 : } else if(Aname=="SC3") {
794 32 : parameter[i].push_back(4.526);
795 32 : parameter[i].push_back(-0.00381305);
796 32 : parameter[i].push_back(5.8567);
797 32 : parameter[i].push_back(-0.214086);
798 16 : parameter[i].push_back(-4.63649);
799 32 : parameter[i].push_back(2.52869);
800 32 : parameter[i].push_back(-0.39894);
801 0 : } else error("Atom name not known");
802 128 : } else if(Rname=="VAL") {
803 128 : if(Aname=="BB") {
804 128 : parameter[i].push_back(10.691);
805 128 : parameter[i].push_back(-0.0162929);
806 128 : parameter[i].push_back(1.24446);
807 128 : parameter[i].push_back(0.307914);
808 128 : parameter[i].push_back(-3.27446);
809 128 : parameter[i].push_back(2.14788);
810 128 : parameter[i].push_back(-0.403259);
811 64 : } else if(Aname=="SC1") {
812 128 : parameter[i].push_back(-3.516);
813 128 : parameter[i].push_back(1.62307);
814 128 : parameter[i].push_back(5.43064);
815 128 : parameter[i].push_back(9.28809);
816 128 : parameter[i].push_back(-14.9927);
817 128 : parameter[i].push_back(6.6133);
818 128 : parameter[i].push_back(-0.964977);
819 0 : } else error("Atom name not known");
820 0 : } else error("Residue not known");
821 : }
822 : } else {
823 0 : error("MOLINFO DATA not found\n");
824 : }
825 8 : }
826 :
827 0 : void SAXS::calculateASF(const vector<AtomNumber> &atoms, vector<vector<long double> > &FF_tmp, const double rho)
828 : {
829 : enum { H, C, N, O, P, S, NTT };
830 : map<string, unsigned> AA_map;
831 0 : AA_map["H"] = H;
832 0 : AA_map["C"] = C;
833 0 : AA_map["N"] = N;
834 0 : AA_map["O"] = O;
835 0 : AA_map["P"] = P;
836 0 : AA_map["S"] = S;
837 :
838 0 : vector<vector<double> > param_a;
839 0 : vector<vector<double> > param_b;
840 : vector<double> param_c;
841 : vector<double> param_v;
842 :
843 0 : param_a.resize(NTT, vector<double>(5));
844 0 : param_b.resize(NTT, vector<double>(5));
845 0 : param_c.resize(NTT);
846 0 : param_v.resize(NTT);
847 :
848 0 : param_a[H][0] = 0.493002; param_b[H][0] = 10.5109; param_c[H] = 0.003038;
849 0 : param_a[H][1] = 0.322912; param_b[H][1] = 26.1257; param_v[H] = 5.15;
850 0 : param_a[H][2] = 0.140191; param_b[H][2] = 3.14236;
851 0 : param_a[H][3] = 0.040810; param_b[H][3] = 57.7997;
852 0 : param_a[H][4] = 0.0; param_b[H][4] = 1.0;
853 :
854 0 : param_a[C][0] = 2.31000; param_b[C][0] = 20.8439; param_c[C] = 0.215600;
855 0 : param_a[C][1] = 1.02000; param_b[C][1] = 10.2075; param_v[C] = 16.44;
856 0 : param_a[C][2] = 1.58860; param_b[C][2] = 0.56870;
857 0 : param_a[C][3] = 0.86500; param_b[C][3] = 51.6512;
858 0 : param_a[C][4] = 0.0; param_b[C][4] = 1.0;
859 :
860 0 : param_a[N][0] = 12.2126; param_b[N][0] = 0.00570; param_c[N] = -11.529;
861 0 : param_a[N][1] = 3.13220; param_b[N][1] = 9.89330; param_v[N] = 2.49;
862 0 : param_a[N][2] = 2.01250; param_b[N][2] = 28.9975;
863 0 : param_a[N][3] = 1.16630; param_b[N][3] = 0.58260;
864 0 : param_a[N][4] = 0.0; param_b[N][4] = 1.0;
865 :
866 0 : param_a[O][0] = 3.04850; param_b[O][0] = 13.2771; param_c[O] = 0.250800 ;
867 0 : param_a[O][1] = 2.28680; param_b[O][1] = 5.70110; param_v[O] = 9.13;
868 0 : param_a[O][2] = 1.54630; param_b[O][2] = 0.32390;
869 0 : param_a[O][3] = 0.86700; param_b[O][3] = 32.9089;
870 0 : param_a[O][4] = 0.0; param_b[O][4] = 1.0;
871 :
872 0 : param_a[P][0] = 6.43450; param_b[P][0] = 1.90670; param_c[P] = 1.11490;
873 0 : param_a[P][1] = 4.17910; param_b[P][1] = 27.1570; param_v[P] = 5.73;
874 0 : param_a[P][2] = 1.78000; param_b[P][2] = 0.52600;
875 0 : param_a[P][3] = 1.49080; param_b[P][3] = 68.1645;
876 0 : param_a[P][4] = 0.0; param_b[P][4] = 1.0;
877 :
878 0 : param_a[S][0] = 6.90530; param_b[S][0] = 1.46790; param_c[S] = 0.866900;
879 0 : param_a[S][1] = 5.20340; param_b[S][1] = 22.2151; param_v[S] = 19.86;
880 0 : param_a[S][2] = 1.43790; param_b[S][2] = 0.25360;
881 0 : param_a[S][3] = 1.58630; param_b[S][3] = 56.1720;
882 0 : param_a[S][4] = 0.0; param_b[S][4] = 1.0;
883 :
884 0 : vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
885 :
886 0 : if( moldat.size()==1 ) {
887 0 : log<<" MOLINFO DATA found, using proper atom names\n";
888 0 : for(unsigned i=0; i<atoms.size(); ++i) {
889 : // get atom name
890 0 : string name = moldat[0]->getAtomName(atoms[i]);
891 : char type;
892 : // get atom type
893 0 : char first = name.at(0);
894 : // GOLDEN RULE: type is first letter, if not a number
895 0 : if (!isdigit(first)) {
896 : type = first;
897 : // otherwise is the second
898 : } else {
899 0 : type = name.at(1);
900 : }
901 0 : std::string type_s = std::string(1,type);
902 0 : if(AA_map.find(type_s) != AA_map.end()) {
903 0 : const unsigned index=AA_map[type_s];
904 0 : const double volr = pow(param_v[index], (2.0/3.0)) /(4. * M_PI);
905 0 : for(unsigned k=0; k<q_list.size(); ++k) {
906 0 : const double q = q_list[k];
907 0 : const double s = q / (4. * M_PI);
908 0 : FF_tmp[k][i] = param_c[index];
909 : // SUM [a_i * EXP( - b_i * (q/4pi)^2 )] Waasmaier and Kirfel (1995)
910 0 : for(unsigned j=0; j<4; j++) {
911 0 : FF_tmp[k][i] += param_a[index][j]*exp(-param_b[index][j]*s*s);
912 : }
913 : // subtract solvation: rho * v_i * EXP( (- v_i^(2/3) / (4pi)) * q^2 ) // since D in Fraser 1978 is 2*s
914 0 : FF_tmp[k][i] -= rho*param_v[index]*exp(-volr*q*q);
915 : }
916 : } else {
917 0 : error("Wrong atom type "+type_s+" from atom name "+name+"\n");
918 : }
919 : }
920 : } else {
921 0 : error("MOLINFO DATA not found\n");
922 : }
923 0 : }
924 :
925 : }
926 4839 : }
|
