Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2018 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 : #define cutOffNB 0.70 // buffer distance for neighbour-lists
24 : #define cutOffDist 0.50 // cut off distance for non-bonded pairwise forces
25 : #define cutOnDist 0.32 // cut off distance for non-bonded pairwise forces
26 : #define cutOffNB2 cutOffNB*cutOffNB // squared buffer distance for neighbour-lists
27 : #define cutOffDist2 cutOffDist*cutOffDist
28 : #define cutOnDist2 cutOnDist*cutOnDist
29 : #define invswitch 1.0/((cutOffDist2-cutOnDist2)*(cutOffDist2-cutOnDist2)*(cutOffDist2-cutOnDist2))
30 : #define cutOffDist4 cutOffDist2*cutOffDist2
31 : #define cutMixed cutOffDist2*cutOffDist2*cutOffDist2 -3.*cutOffDist2*cutOffDist2*cutOnDist2
32 :
33 : #include <string>
34 : #include <fstream>
35 : #include <iterator>
36 : #include <sstream>
37 :
38 : #include "Colvar.h"
39 : #include "ActionRegister.h"
40 : #include "core/PlumedMain.h"
41 : #include "tools/OpenMP.h"
42 : #include "tools/Pbc.h"
43 : #include "tools/PDB.h"
44 : #include "tools/Torsion.h"
45 :
46 : using namespace std;
47 :
48 : namespace PLMD {
49 : namespace colvar {
50 :
51 : //+PLUMEDOC COLVAR CS2BACKBONE
52 : /*
53 : This collective variable calculates the backbone chemical shifts for a protein.
54 :
55 : The functional form is that of CamShift \cite Kohlhoff:2009us. The chemical shifts
56 : of the selected nuclei/residues are saved as components. Reference experimental values
57 : can also be stored as components. The two components can then be used to calculate
58 : either a scoring function as in \cite Robustelli:2010dn \cite Granata:2013dk, using
59 : the keyword CAMSHIFT or to calculate ensemble averages chemical shift per chemical
60 : shift as in \cite Camilloni:2012je \cite Camilloni:2013hs (see \ref STATS and
61 : \ref ENSEMBLE).
62 :
63 : CamShift calculation is relatively heavy because it often uses a large number of atoms, in order
64 : to make it faster it is currently parallelised with \ref Openmp.
65 :
66 : As a general rule, when using \ref CS2BACKBONE or other experimental restraints it is better to
67 : increase the accuracy of the constraint algorithm due to the increased strain on the bonded structure.
68 : In the case of GROMACS it is safer to use lincs-iter=2 and lincs-order=6.
69 :
70 : In general the system for which chemical shifts are calculated must be completly included in
71 : ATOMS and a TEMPLATE pdb file for the same atoms should be provided as well in the folder DATA.
72 : The atoms are made automatically whole unless NOPBC is used, in particular if the system is made of
73 : by multiple chains it is usually better to use NOPBC and make the molecule whole \ref WHOLEMOLECULES
74 : selecting an appropriate order.
75 :
76 : In addition to a pdb file one needs to provide a list of chemical shifts to be calculated using one
77 : file per nucleus type (CAshifts.dat, CBshifts.dat, Cshifts.dat, Hshifts.dat, HAshifts.dat, Nshifts.dat),
78 : all the six files should always be present. A chemical shift for a nucleus is calculated if a value
79 : greater than 0 is provided. For practical purposes the value can correspond to the experimental value.
80 : Residues numbers should go from 1 to N irrespectively of the numbers used in the pdb file. The first and
81 : last residue of each chain should be preceeded by a # character. Termini groups like ACE or NME should
82 : be removed from the PDB.
83 :
84 : \verbatim
85 : CAshifts.dat:
86 : #1 0.0
87 : 2 55.5
88 : 3 58.4
89 : .
90 : .
91 : #last 0.0
92 : #last+1 (first) of second chain
93 : .
94 : #last of second chain
95 : \endverbatim
96 :
97 : The default behaviour is to store the values for the active nuclei in components (ca_#, cb_#,
98 : co_#, ha_#, hn_#, nh_# and expca_#, expcb_#, expco_#, expha_#, exphn_#, exp_nh#) with NOEXP it is possible
99 : to only store the backcalculated values.
100 :
101 : A pdb file is needed to the generate a simple topology of the protein. For histidines in protonation
102 : states different from D the HIE/HSE HIP/HSP name should be used. GLH and ASH can be used for the alternative
103 : protonation of GLU and ASP. Non-standard amino acids and other molecules are not yet supported, but in principle
104 : they can be named UNK. If multiple chains are present the chain identifier must be in the standard PDB format,
105 : together with the TER keyword at the end of each chain.
106 :
107 : One more standard file is also needed in the folder DATA: camshift.db. This file includes all the CamShift parameters
108 : and can be found in regtest/basic/rt45/data/ .
109 :
110 : All the above files must be in a single folder that must be specified with the keyword DATA.
111 :
112 : Additional material and examples can be also found in the tutorial \ref belfast-9
113 :
114 : \par Examples
115 :
116 : In this first example the chemical shifts are used to calculate a scoring function to be used
117 : in NMR driven Metadynamics \cite Granata:2013dk :
118 :
119 : \verbatim
120 : whole: GROUP ATOMS=2612-2514:-1,961-1:-1,2466-962:-1,2513-2467:-1
121 : WHOLEMOLECULES ENTITY0=whole
122 : cs: CS2BACKBONE ATOMS=1-2612 NRES=176 DATA=../data/ TEMPLATE=template.pdb CAMSHIFT NOPBC
123 : metad: METAD ARG=cs HEIGHT=0.5 SIGMA=0.1 PACE=200 BIASFACTOR=10
124 : PRINT ARG=cs,metad.bias FILE=COLVAR STRIDE=100
125 : \endverbatim
126 :
127 : In this second example the chemical shifts are used as replica-averaged restrained as in \cite Camilloni:2012je \cite Camilloni:2013hs.
128 :
129 : \verbatim
130 : cs: CS2BACKBONE ATOMS=1-174 DATA=data/ NRES=13
131 : encs: ENSEMBLE ARG=(cs\.hn_.*),(cs\.nh_.*)
132 : stcs: STATS ARG=encs.* SQDEVSUM PARARG=(cs\.exphn_.*),(cs\.expnh_.*)
133 : RESTRAINT ARG=stcs.sqdevsum AT=0 KAPPA=0 SLOPE=24
134 :
135 : PRINT ARG=(cs\.hn_.*),(cs\.nh_.*) FILE=RESTRAINT STRIDE=100
136 :
137 : \endverbatim
138 :
139 : (See also \ref WHOLEMOLECULES, \ref STATS, \ref METAD, \ref RESTRAINT and \ref PRINT)
140 :
141 : */
142 : //+ENDPLUMEDOC
143 :
144 : class CS2BackboneDB {
145 : enum { STD, GLY, PRO};
146 : enum { HA_ATOM, H_ATOM, N_ATOM, CA_ATOM, CB_ATOM, C_ATOM };
147 : static const unsigned aa_kind = 3;
148 : static const unsigned atm_kind = 6;
149 : static const unsigned numXtraDists = 27;
150 :
151 : // ALA, ARG, ASN, ASP, CYS, GLU, GLN, GLY, HIS, ILE, LEU, LYS, MET, PHE, PRO, SER, THR, TRP, TYR, VAL
152 : double c_aa[aa_kind][atm_kind][20];
153 : double c_aa_prev[aa_kind][atm_kind][20];
154 : double c_aa_succ[aa_kind][atm_kind][20];
155 : double co_bb[aa_kind][atm_kind][16];
156 : double co_sc_[aa_kind][atm_kind][20][20];
157 : double co_xd[aa_kind][atm_kind][numXtraDists];
158 : double co_sphere[aa_kind][atm_kind][2][8];
159 : // for ring current effects
160 : // Phe, Tyr, Trp_1, Trp_2, His
161 : double co_ring[aa_kind][atm_kind][5];
162 : // for dihedral angles
163 : // co * (a * cos(3 * omega + c) + b * cos(omega + d))
164 : double co_da[aa_kind][atm_kind][3];
165 : double pars_da[aa_kind][atm_kind][3][5];
166 :
167 : public:
168 :
169 1164 : inline unsigned kind(const string &s) {
170 1164 : if(s=="GLY") return GLY;
171 948 : else if(s=="PRO") return PRO;
172 870 : return STD;
173 : }
174 :
175 108 : inline unsigned atom_kind(const string &s) {
176 108 : if(s=="HA")return HA_ATOM;
177 90 : else if(s=="H") return H_ATOM;
178 72 : else if(s=="N") return N_ATOM;
179 54 : else if(s=="CA")return CA_ATOM;
180 36 : else if(s=="CB")return CB_ATOM;
181 18 : else if(s=="C") return C_ATOM;
182 0 : return -1;
183 : }
184 :
185 27864 : unsigned get_numXtraDists() {return numXtraDists;}
186 :
187 : //PARAMETERS
188 3533 : inline double * CONSTAACURR(const unsigned a_kind, const unsigned at_kind) {return c_aa[a_kind][at_kind];}
189 3534 : inline double * CONSTAANEXT(const unsigned a_kind, const unsigned at_kind) {return c_aa_succ[a_kind][at_kind];}
190 3534 : inline double * CONSTAAPREV(const unsigned a_kind, const unsigned at_kind) {return c_aa_prev[a_kind][at_kind];}
191 3534 : inline double * CONST_BB2_PREV(const unsigned a_kind, const unsigned at_kind) {return co_bb[a_kind][at_kind];}
192 3533 : inline double * CONST_BB2_CURR(const unsigned a_kind, const unsigned at_kind) {return co_bb[a_kind][at_kind]+5;}
193 3533 : inline double * CONST_BB2_NEXT(const unsigned a_kind, const unsigned at_kind) {return co_bb[a_kind][at_kind]+11;}
194 3534 : inline double * CONST_SC2(const unsigned a_kind, const unsigned at_kind, unsigned res_type) { return co_sc_[a_kind][at_kind][res_type];}
195 3534 : inline double * CONST_XD(const unsigned a_kind, const unsigned at_kind) { return co_xd[a_kind][at_kind];}
196 7068 : inline double * CO_SPHERE(const unsigned a_kind, const unsigned at_kind, unsigned exp_type) { return co_sphere[a_kind][at_kind][exp_type];}
197 3534 : inline double * CO_RING(const unsigned a_kind, const unsigned at_kind) { return co_ring[a_kind][at_kind];}
198 3534 : inline double * CO_DA(const unsigned a_kind, const unsigned at_kind) { return co_da[a_kind][at_kind];}
199 9887 : inline double * PARS_DA(const unsigned a_kind, const unsigned at_kind, const unsigned ang_kind) { return pars_da[a_kind][at_kind][ang_kind];}
200 :
201 6 : void parse(const string &file, const double dscale) {
202 6 : ifstream in;
203 6 : in.open(file.c_str());
204 6 : if(!in) plumed_merror("Unable to open CS2Backbone DB file " +file);
205 :
206 6 : unsigned c_kind = 0;
207 6 : unsigned c_atom = 0;
208 6 : unsigned nline = 0;
209 :
210 114 : for(unsigned i=0; i<3; i++) for(unsigned j=0; j<6; j++) {
211 2268 : for(unsigned k=0; k<20; k++) {
212 2160 : c_aa[i][j][k]=0.;
213 2160 : c_aa_prev[i][j][k]=0.;
214 2160 : c_aa_succ[i][j][k]=0.;
215 2160 : for(unsigned m=0; m<20; m++) co_sc_[i][j][k][m]=0.;
216 : }
217 108 : for(unsigned k=0; k<16; k++) {co_bb[i][j][k]=0.; }
218 108 : for(unsigned k=0; k<8; k++) { co_sphere[i][j][0][k]=0.; co_sphere[i][j][1][k]=0.; }
219 432 : for(unsigned k=0; k<3; k++) {
220 324 : co_da[i][j][k]=0.;
221 324 : for(unsigned l=0; l<5; l++) pars_da[i][j][k][l]=0.;
222 : }
223 108 : for(unsigned k=0; k<5; k++) co_ring[i][j][k]=0.;
224 108 : for(unsigned k=0; k<numXtraDists; k++) co_xd[i][j][k]=0.;
225 : }
226 :
227 12576 : while(!in.eof()) {
228 12564 : string line;
229 12564 : getline(in,line);
230 12564 : ++nline;
231 12564 : if(line.compare(0,1,"#")==0) continue;
232 6810 : vector<string> tok;
233 6810 : vector<string> tmp;
234 6810 : tok = split(line,' ');
235 87276 : for(unsigned q=0; q<tok.size(); q++)
236 80466 : if(tok[q].size()) tmp.push_back(tok[q]);
237 6810 : tok = tmp;
238 6810 : if(tok.size()==0) continue;
239 6804 : if(tok[0]=="PAR") {
240 108 : c_kind = kind(tok[2]);
241 108 : c_atom = atom_kind(tok[1]);
242 108 : continue;
243 : }
244 6696 : else if(tok[0]=="WEIGHT") {
245 108 : continue;
246 : }
247 6588 : else if(tok[0]=="FLATBTM") {
248 108 : continue;
249 : }
250 6480 : else if (tok[0] == "SCALEHARM") {
251 108 : continue;
252 : }
253 6372 : else if (tok[0] == "TANHAMPLI") {
254 108 : continue;
255 : }
256 6264 : else if (tok[0] == "ENDHARMON") {
257 108 : continue;
258 : }
259 6156 : else if (tok[0] == "MAXRCDEVI") {
260 108 : continue;
261 : }
262 6048 : else if (tok[0] == "RANDCOIL") {
263 108 : continue;
264 : }
265 5940 : else if (tok[0] == "CONST") {
266 108 : continue;
267 : }
268 5832 : else if (tok[0] == "CONSTAA") {
269 108 : assign(c_aa[c_kind][c_atom],tok,1);
270 108 : continue;
271 : }
272 5724 : else if (tok[0] == "CONSTAA-1") {
273 108 : assign(c_aa_prev[c_kind][c_atom],tok,1);
274 108 : continue;
275 : }
276 5616 : else if (tok[0] == "CONSTAA+1") {
277 108 : assign(c_aa_succ[c_kind][c_atom],tok,1);
278 108 : continue;
279 : }
280 5508 : else if (tok[0] == "COBB1") {
281 108 : continue;
282 : }
283 5400 : else if (tok[0] == "COBB2") {
284 : //angstrom to nm
285 108 : assign(co_bb[c_kind][c_atom],tok,dscale);
286 108 : continue;
287 : }
288 5292 : else if (tok[0] == "SPHERE1") {
289 : // angstrom^-3 to nm^-3
290 108 : assign(co_sphere[c_kind][c_atom][0],tok,1./(dscale*dscale*dscale));
291 108 : continue;
292 : }
293 5184 : else if (tok[0] == "SPHERE2") {
294 : //angstrom to nm
295 108 : assign(co_sphere[c_kind][c_atom][1],tok,dscale);
296 108 : continue;
297 : }
298 5076 : else if (tok[0] == "DIHEDRALS") {
299 108 : assign(co_da[c_kind][c_atom],tok,1);
300 108 : continue;
301 : }
302 4968 : else if (tok[0] == "RINGS") {
303 : // angstrom^-3 to nm^-3
304 108 : assign(co_ring[c_kind][c_atom],tok,1./(dscale*dscale*dscale));
305 648 : for(unsigned i=1; i<tok.size(); i++)
306 540 : co_ring[c_kind][c_atom][i-1] *= 1000;
307 108 : continue;
308 : }
309 4860 : else if (tok[0] == "HBONDS") {
310 108 : continue;
311 : }
312 4752 : else if (tok[0] == "XTRADISTS") {
313 : //angstrom to nm
314 108 : assign(co_xd[c_kind][c_atom],tok,dscale);
315 108 : continue;
316 : }
317 4644 : else if(tok[0]=="DIHEDPHI") {
318 108 : assign(pars_da[c_kind][c_atom][0],tok,1);
319 108 : continue;
320 : }
321 4536 : else if(tok[0]=="DIHEDPSI") {
322 108 : assign(pars_da[c_kind][c_atom][1],tok,1);
323 108 : continue;
324 : }
325 4428 : else if(tok[0]=="DIHEDCHI1") {
326 108 : assign(pars_da[c_kind][c_atom][2],tok,1);
327 108 : continue;
328 : }
329 :
330 4320 : bool ok = false;
331 : string scIdent1 [] = {"COSCALA1", "COSCARG1", "COSCASN1", "COSCASP1", "COSCCYS1", "COSCGLN1", "COSCGLU1",
332 : "COSCGLY1", "COSCHIS1", "COSCILE1", "COSCLEU1", "COSCLYS1", "COSCMET1", "COSCPHE1",
333 : "COSCPRO1", "COSCSER1", "COSCTHR1", "COSCTRP1", "COSCTYR1", "COSCVAL1"
334 90720 : };
335 :
336 68040 : for(unsigned scC = 0; scC < 20; scC++) {
337 65880 : if(tok[0]==scIdent1[scC]) {
338 2160 : ok = true;
339 2160 : break;
340 : }
341 : }
342 4320 : if(ok) continue;
343 :
344 : string scIdent2 [] = {"COSCALA2", "COSCARG2", "COSCASN2", "COSCASP2", "COSCCYS2", "COSCGLN2", "COSCGLU2",
345 : "COSCGLY2", "COSCHIS2", "COSCILE2", "COSCLEU2", "COSCLYS2", "COSCMET2", "COSCPHE2",
346 : "COSCPRO2", "COSCSER2", "COSCTHR2", "COSCTRP2", "COSCTYR2", "COSCVAL2"
347 45360 : };
348 :
349 22680 : for(unsigned scC = 0; scC < 20; scC++) {
350 22680 : if(tok[0]==scIdent2[scC]) {
351 : //angstrom to nm
352 2160 : assign(co_sc_[c_kind][c_atom][scC],tok,dscale);
353 2160 : ok = true; break;
354 : }
355 : }
356 2160 : if(ok) continue;
357 :
358 0 : if(tok.size()) {
359 0 : string str_err = "CS2Backbone DB WARNING: unrecognized token: " + tok[0];
360 0 : plumed_merror(str_err);
361 : }
362 0 : }
363 6 : in.close();
364 6 : }
365 :
366 : private:
367 :
368 6810 : vector<string> &split(const string &s, char delim, vector<string> &elems) {
369 6810 : stringstream ss(s);
370 13620 : string item;
371 94086 : while (getline(ss, item, delim)) {
372 80466 : elems.push_back(item);
373 : }
374 13620 : return elems;
375 : }
376 :
377 6810 : vector<string> split(const string &s, char delim) {
378 6810 : vector<string> elems;
379 6810 : split(s, delim, elems);
380 6810 : return elems;
381 : }
382 :
383 3456 : void assign(double * f, const vector<string> & v, const double scale) {
384 40932 : for(unsigned i=1; i<v.size(); i++)
385 37476 : f[i-1] = scale*(atof(v[i].c_str()));
386 3456 : }
387 : };
388 :
389 12 : class CS2Backbone : public Colvar {
390 10500 : struct Fragment {
391 : vector<Value*> comp;
392 : vector<double> exp_cs;
393 : unsigned res_type_prev;
394 : unsigned res_type_curr;
395 : unsigned res_type_next;
396 : unsigned res_kind;
397 : unsigned fd;
398 : string res_name;
399 : vector<int> pos;
400 : vector<int> prev;
401 : vector<int> curr;
402 : vector<int> next;
403 : vector<int> side_chain;
404 : vector<int> xd1;
405 : vector<int> xd2;
406 : vector<unsigned> box_nb;
407 : vector<int> phi;
408 : vector<int> psi;
409 : vector<int> chi1;
410 : double t_phi;
411 : double t_psi;
412 : double t_chi1;
413 : vector<Vector> dd0, dd10, dd21, dd2;
414 :
415 1056 : Fragment() {
416 1056 : comp.resize(6);
417 1056 : exp_cs.resize(6,0);
418 1056 : res_type_prev = res_type_curr = res_type_next = 0;
419 1056 : res_kind = 0;
420 1056 : fd = 0;
421 1056 : res_name = "";
422 1056 : pos.resize(6,-1);
423 1056 : prev.reserve(5);
424 1056 : curr.reserve(6);
425 1056 : next.reserve(5);
426 1056 : side_chain.reserve(20);
427 1056 : xd1.reserve(27);
428 1056 : xd2.reserve(27);
429 1056 : box_nb.reserve(250);
430 1056 : phi.reserve(4);
431 1056 : psi.reserve(4);
432 1056 : chi1.reserve(4);
433 1056 : t_phi = t_psi = t_chi1 = 0;
434 1056 : dd0.resize(3);
435 1056 : dd10.resize(3);
436 1056 : dd21.resize(3);
437 1056 : dd2.resize(3);
438 1056 : }
439 :
440 : };
441 :
442 : struct RingInfo {
443 : enum {R_PHE, R_TYR, R_TRP1, R_TRP2, R_HIS};
444 : unsigned rtype; // one out of five different types
445 : unsigned atom[6]; // up to six member per ring
446 : unsigned numAtoms; // number of ring members (5 or 6)
447 : Vector position; // center of ring coordinates
448 : Vector normVect; // ring plane normal vector
449 : Vector n1, n2; // two atom plane normal vectors used to compute ring plane normal
450 : Vector g[6]; // vector of the vectors used to calculate n1,n2
451 : double lengthN2; // square of length of normVect
452 : double lengthNV; // length of normVect
453 120 : RingInfo():
454 : rtype(0),numAtoms(0),
455 120 : lengthN2(NAN),lengthNV(NAN)
456 120 : {for(unsigned i=0; i<6; i++) atom[i]=0;}
457 : };
458 :
459 : enum aa_t {ALA, ARG, ASN, ASP, CYS, GLN, GLU, GLY, HIS, ILE, LEU, LYS, MET, PHE, PRO, SER, THR, TRP, TYR, VAL, UNK};
460 : enum atom_t {D_C, D_H, D_N, D_O, D_S, D_C2, D_N2, D_O2};
461 :
462 : CS2BackboneDB db;
463 : vector<vector<Fragment> > atom;
464 : vector<RingInfo> ringInfo;
465 : vector<unsigned> seg_last;
466 : vector<unsigned> type;
467 : vector<unsigned> res_num;
468 : unsigned box_nupdate;
469 : unsigned box_count;
470 : bool camshift;
471 : bool pbc;
472 :
473 : void remove_problematic(const string &res, const string &nucl);
474 : void read_cs(const string &file, const string &k);
475 : void update_neighb();
476 : void compute_ring_parameters();
477 : void compute_dihedrals();
478 : void init_backbone(const PDB &pdb);
479 : void init_sidechain(const PDB &pdb);
480 : void init_xdist(const PDB &pdb);
481 : void init_types(const PDB &pdb);
482 : void init_rings(const PDB &pdb);
483 : aa_t frag2enum(const string &aa);
484 : vector<string> side_chain_atoms(const string &s);
485 : bool isSP2(const string & resType, const string & atomName);
486 : bool is_chi1_cx(const string & frg, const string & atm);
487 : unsigned frag_segment(const unsigned p);
488 : unsigned frag_relitive_index(const unsigned p, const unsigned s);
489 : void debug_report();
490 : void xdist_name_map(string & name);
491 :
492 : public:
493 :
494 : explicit CS2Backbone(const ActionOptions&);
495 : static void registerKeywords( Keywords& keys );
496 : virtual void calculate();
497 : };
498 :
499 2529 : PLUMED_REGISTER_ACTION(CS2Backbone,"CS2BACKBONE")
500 :
501 7 : void CS2Backbone::registerKeywords( Keywords& keys ) {
502 7 : Colvar::registerKeywords( keys );
503 7 : componentsAreNotOptional(keys);
504 7 : useCustomisableComponents(keys);
505 7 : keys.add("atoms","ATOMS","The atoms to be included in the calculation, e.g. the whole protein.");
506 7 : keys.add("compulsory","DATA","data/","The folder with the experimental chemical shifts.");
507 7 : keys.add("compulsory","TEMPLATE","template.pdb","A PDB file of the protein system to initialise ALMOST.");
508 7 : keys.add("compulsory","NEIGH_FREQ","20","Period in step for neighbour list update.");
509 7 : keys.add("compulsory","NRES","Number of residues, corresponding to the number of chemical shifts.");
510 7 : keys.addFlag("CAMSHIFT",false,"Set to TRUE if you to calculate a single CamShift score.");
511 7 : keys.addFlag("NOEXP",false,"Set to TRUE if you don't want to have fixed components with the experimetnal values.");
512 7 : keys.addOutputComponent("ha","default","the calculated Ha hydrogen chemical shifts");
513 7 : keys.addOutputComponent("hn","default","the calculated H hydrogen chemical shifts");
514 7 : keys.addOutputComponent("nh","default","the calculated N nitrogen chemical shifts");
515 7 : keys.addOutputComponent("ca","default","the calculated Ca carbon chemical shifts");
516 7 : keys.addOutputComponent("cb","default","the calculated Cb carbon chemical shifts");
517 7 : keys.addOutputComponent("co","default","the calculated C' carbon chemical shifts");
518 7 : keys.addOutputComponent("expha","default","the experimental Ha hydrogen chemical shifts");
519 7 : keys.addOutputComponent("exphn","default","the experimental H hydrogen chemical shifts");
520 7 : keys.addOutputComponent("expnh","default","the experimental N nitrogen chemical shifts");
521 7 : keys.addOutputComponent("expca","default","the experimental Ca carbon chemical shifts");
522 7 : keys.addOutputComponent("expcb","default","the experimental Cb carbon chemical shifts");
523 7 : keys.addOutputComponent("expco","default","the experimental C' carbon chemical shifts");
524 7 : }
525 :
526 6 : CS2Backbone::CS2Backbone(const ActionOptions&ao):
527 : PLUMED_COLVAR_INIT(ao),
528 : camshift(false),
529 6 : pbc(true)
530 : {
531 6 : string stringadb;
532 12 : string stringapdb;
533 :
534 6 : parseFlag("CAMSHIFT",camshift);
535 :
536 6 : bool nopbc=!pbc;
537 6 : parseFlag("NOPBC",nopbc);
538 6 : pbc=!nopbc;
539 :
540 6 : bool noexp=false;
541 6 : parseFlag("NOEXP",noexp);
542 :
543 12 : string stringa_data;
544 6 : parse("DATA",stringa_data);
545 :
546 12 : string stringa_template;
547 6 : parse("TEMPLATE",stringa_template);
548 :
549 6 : box_count=0;
550 6 : box_nupdate=20;
551 6 : parse("NEIGH_FREQ", box_nupdate);
552 :
553 : unsigned numResidues;
554 6 : parse("NRES", numResidues);
555 :
556 6 : stringadb = stringa_data + string("/camshift.db");
557 6 : stringapdb = stringa_data + string("/") + stringa_template;
558 :
559 : /* Lenght conversion (parameters are tuned for angstrom) */
560 6 : double scale=1.;
561 6 : if(!plumed.getAtoms().usingNaturalUnits()) {
562 6 : scale = 10.*atoms.getUnits().getLength();
563 : }
564 :
565 6 : log.printf(" Initialization of the predictor ...\n"); log.flush();
566 6 : db.parse(stringadb,scale);
567 12 : PDB pdb;
568 6 : if( !pdb.read(stringapdb,plumed.getAtoms().usingNaturalUnits(),1./scale) ) error("missing input file " + stringapdb);
569 6 : init_backbone(pdb);
570 6 : init_sidechain(pdb);
571 6 : init_xdist(pdb);
572 6 : init_types(pdb);
573 6 : init_rings(pdb);
574 : #ifndef NDEBUG
575 : debug_report();
576 : #endif
577 6 : log.printf(" Reading experimental data ...\n"); log.flush();
578 6 : stringadb = stringa_data + string("/CAshifts.dat");
579 6 : log.printf(" Initializing CA shifts %s\n", stringadb.c_str()); log.flush();
580 6 : read_cs(stringadb, "CA");
581 6 : stringadb = stringa_data + string("/CBshifts.dat");
582 6 : log.printf(" Initializing CB shifts %s\n", stringadb.c_str()); log.flush();
583 6 : read_cs(stringadb, "CB");
584 6 : stringadb = stringa_data + string("/Cshifts.dat");
585 6 : log.printf(" Initializing C' shifts %s\n", stringadb.c_str()); log.flush();
586 6 : read_cs(stringadb, "C");
587 6 : stringadb = stringa_data + string("/HAshifts.dat");
588 6 : log.printf(" Initializing HA shifts %s\n", stringadb.c_str()); log.flush();
589 6 : read_cs(stringadb, "HA");
590 6 : stringadb = stringa_data + string("/Hshifts.dat");
591 6 : log.printf(" Initializing H shifts %s\n", stringadb.c_str()); log.flush();
592 6 : read_cs(stringadb, "H");
593 6 : stringadb = stringa_data + string("/Nshifts.dat");
594 6 : log.printf(" Initializing N shifts %s\n", stringadb.c_str()); log.flush();
595 6 : read_cs(stringadb, "N");
596 :
597 : /* this is a workaround for those chemical shifts that can result in too large forces */
598 6 : remove_problematic("GLN", "CB");
599 6 : remove_problematic("ILE", "CB");
600 6 : remove_problematic("PRO", "N");
601 6 : remove_problematic("PRO", "H");
602 6 : remove_problematic("PRO", "CB");
603 6 : remove_problematic("GLY", "HA");
604 6 : remove_problematic("GLY", "CB");
605 : /* this is a workaround for those chemical shifts that are not parameterized */
606 6 : remove_problematic("HIE", "HA"); remove_problematic("HIP", "HA"); remove_problematic("HSP", "HA");
607 6 : remove_problematic("HIE", "H"); remove_problematic("HIP", "H"); remove_problematic("HSP", "H");
608 6 : remove_problematic("HIE", "N"); remove_problematic("HIP", "N"); remove_problematic("HSP", "N");
609 6 : remove_problematic("HIE", "CA"); remove_problematic("HIP", "CA"); remove_problematic("HSP", "CA");
610 6 : remove_problematic("HIE", "CB"); remove_problematic("HIP", "CB"); remove_problematic("HSP", "CB");
611 6 : remove_problematic("HIE", "C"); remove_problematic("HIP", "C"); remove_problematic("HSP", "C");
612 6 : remove_problematic("GLH", "HA"); remove_problematic("ASH", "HA"); remove_problematic("HSE", "HA");
613 6 : remove_problematic("GLH", "H"); remove_problematic("ASH", "H"); remove_problematic("HSE", "H");
614 6 : remove_problematic("GLH", "N"); remove_problematic("ASH", "N"); remove_problematic("HSE", "N");
615 6 : remove_problematic("GLH", "CA"); remove_problematic("ASH", "CA"); remove_problematic("HSE", "CA");
616 6 : remove_problematic("GLH", "CB"); remove_problematic("ASH", "CB"); remove_problematic("HSE", "CB");
617 6 : remove_problematic("GLH", "C"); remove_problematic("ASH", "C"); remove_problematic("HSE", "C");
618 6 : remove_problematic("UNK", "HA");
619 6 : remove_problematic("UNK", "H");
620 6 : remove_problematic("UNK", "N");
621 6 : remove_problematic("UNK", "CA");
622 6 : remove_problematic("UNK", "CB");
623 6 : remove_problematic("UNK", "C");
624 6 : remove_problematic("CYS", "HA");
625 6 : remove_problematic("CYS", "H");
626 6 : remove_problematic("CYS", "N");
627 6 : remove_problematic("CYS", "CA");
628 6 : remove_problematic("CYS", "CB");
629 6 : remove_problematic("CYS", "C");
630 6 : remove_problematic("HIS", "HA");
631 6 : remove_problematic("HIS", "H");
632 6 : remove_problematic("HIS", "N");
633 6 : remove_problematic("HIS", "CA");
634 6 : remove_problematic("HIS", "CB");
635 6 : remove_problematic("HIS", "C");
636 : /* done */
637 :
638 12 : vector<AtomNumber> atoms;
639 6 : parseAtomList("ATOMS",atoms);
640 6 : checkRead();
641 :
642 6 : log<<" Bibliography "
643 18 : <<plumed.cite("Kohlhoff K, Robustelli P, Cavalli A, Salvatella A, Vendruscolo M, J. Am. Chem. Soc. 131, 13894 (2009)")
644 18 : <<plumed.cite("Camilloni C, Robustelli P, De Simone A, Cavalli A, Vendruscolo M, J. Am. Chem. Soc. 134, 3968 (2012)")
645 18 : <<plumed.cite("Granata D, Camilloni C, Vendruscolo M, Laio A, Proc. Natl. Acad. Sci. USA 110, 6817 (2013)")<<"\n";
646 :
647 12 : const string str_cs[] = {"ha_","hn_","nh_","ca_","cb_","co_"};
648 6 : unsigned index=0;
649 6 : if(camshift) noexp = true;
650 6 : if(!camshift) {
651 15 : for(unsigned i=0; i<atom.size(); i++) {
652 890 : for(unsigned a=0; a<atom[i].size(); a++) {
653 880 : unsigned res=index+a;
654 880 : std::string num; Tools::convert(res,num);
655 6160 : for(unsigned at_kind=0; at_kind<6; at_kind++) {
656 5280 : if(atom[i][a].exp_cs[at_kind]!=0) {
657 2945 : addComponentWithDerivatives(str_cs[at_kind]+num);
658 2945 : componentIsNotPeriodic(str_cs[at_kind]+num);
659 2945 : atom[i][a].comp[at_kind] = getPntrToComponent(str_cs[at_kind]+num);
660 : }
661 : }
662 880 : }
663 10 : index += atom[i].size();
664 : }
665 : } else {
666 1 : addValueWithDerivatives();
667 1 : setNotPeriodic();
668 3 : for(unsigned i=0; i<atom.size(); i++) {
669 2 : index += atom[i].size();
670 : }
671 : }
672 :
673 6 : if(!noexp) {
674 5 : index = 0;
675 15 : for(unsigned i=0; i<atom.size(); i++) {
676 890 : for(unsigned a=0; a<atom[i].size(); a++) {
677 880 : unsigned res=index+a;
678 880 : std::string num; Tools::convert(res,num);
679 6160 : for(unsigned at_kind=0; at_kind<6; at_kind++) {
680 5280 : if(atom[i][a].exp_cs[at_kind]!=0) {
681 2945 : addComponent("exp"+str_cs[at_kind]+num);
682 2945 : componentIsNotPeriodic("exp"+str_cs[at_kind]+num);
683 2945 : Value* comp=getPntrToComponent("exp"+str_cs[at_kind]+num);
684 2945 : comp->set(atom[i][a].exp_cs[at_kind]);
685 : }
686 : }
687 880 : }
688 10 : index += atom[i].size();
689 : }
690 : }
691 :
692 : /* temporary check, the idea is that I can remove NRES completely */
693 6 : if(index!=numResidues) error("NRES and the number of residues in the PDB do not match!");
694 :
695 12 : requestAtoms(atoms);
696 6 : }
697 :
698 366 : void CS2Backbone::remove_problematic(const string &res, const string &nucl) {
699 : unsigned n;
700 366 : if(nucl=="HA") n=0;
701 306 : else if(nucl=="H") n=1;
702 246 : else if(nucl=="N") n=2;
703 186 : else if(nucl=="CA")n=3;
704 132 : else if(nucl=="CB")n=4;
705 54 : else if(nucl=="C") n=5;
706 366 : else return;
707 :
708 1098 : for(unsigned i=0; i<atom.size(); i++) {
709 65148 : for(unsigned a=0; a<atom[i].size(); a++) {
710 64416 : if(atom[i][a].res_name.c_str()==res) {
711 708 : atom[i][a].exp_cs[n] = 0;
712 : }
713 : }
714 : }
715 : }
716 :
717 36 : void CS2Backbone::read_cs(const string &file, const string &nucl) {
718 36 : ifstream in;
719 36 : in.open(file.c_str());
720 36 : if(!in) error("CS2Backbone: Unable to open " + file);
721 72 : istream_iterator<string> iter(in), end;
722 :
723 : unsigned n;
724 36 : if(nucl=="HA") n=0;
725 30 : else if(nucl=="H") n=1;
726 24 : else if(nucl=="N") n=2;
727 18 : else if(nucl=="CA")n=3;
728 12 : else if(nucl=="CB")n=4;
729 6 : else if(nucl=="C") n=5;
730 36 : else return;
731 :
732 36 : int oldseg = -1;
733 36 : int oldp = -1;
734 6408 : while(iter!=end) {
735 6336 : string tok;
736 6336 : tok = *iter; ++iter;
737 6336 : if(tok[0]=='#') { ++iter; continue;}
738 6192 : unsigned p = atoi(tok.c_str());
739 6192 : p = p - 1;
740 6192 : const unsigned seg = frag_segment(p);
741 6192 : p = frag_relitive_index(p,seg);
742 6192 : if(oldp==-1) oldp=p;
743 6192 : if(oldseg==-1) oldseg=seg;
744 6192 : if(p<oldp&&seg==oldseg) {
745 0 : string errmsg = "Error while reading " + file + "! The same residue number has been used twice!";
746 0 : error(errmsg);
747 : }
748 6192 : tok = *iter; ++iter;
749 6192 : double cs = atof(tok.c_str());
750 6192 : if(atom[seg][p].pos[n]<=0) cs=0;
751 6018 : else atom[seg][p].exp_cs[n] = cs;
752 6192 : oldseg = seg;
753 6192 : oldp = p;
754 6192 : }
755 72 : in.close();
756 : }
757 :
758 6 : void CS2Backbone::calculate()
759 : {
760 6 : if(pbc) makeWhole();
761 6 : if(getExchangeStep()) box_count=0;
762 6 : if(box_count==0) update_neighb();
763 :
764 6 : compute_ring_parameters();
765 6 : compute_dihedrals();
766 :
767 6 : double score = 0.;
768 :
769 6 : vector<double> camshift_sigma2(6);
770 6 : camshift_sigma2[0] = 0.08; // HA
771 6 : camshift_sigma2[1] = 0.30; // HN
772 6 : camshift_sigma2[2] = 9.00; // NH
773 6 : camshift_sigma2[3] = 1.30; // CA
774 6 : camshift_sigma2[4] = 1.56; // CB
775 6 : camshift_sigma2[5] = 1.70; // CO
776 :
777 6 : const unsigned chainsize = atom.size();
778 6 : const unsigned atleastned = 72+ringInfo.size()*6;
779 :
780 : // CYCLE OVER MULTIPLE CHAINS
781 17 : #pragma omp parallel num_threads(OpenMP::getNumThreads())
782 35 : for(unsigned s=0; s<chainsize; s++) {
783 23 : const unsigned psize = atom[s].size();
784 23 : vector<Vector> omp_deriv;
785 27 : if(camshift) omp_deriv.resize(getNumberOfAtoms(), Vector(0,0,0));
786 50 : #pragma omp for reduction(+:score)
787 : // SKIP FIRST AND LAST RESIDUE OF EACH CHAIN
788 25 : for(unsigned a=1; a<psize-1; a++) {
789 :
790 1033 : const Fragment *myfrag = &atom[s][a];
791 1031 : const unsigned aa_kind = myfrag->res_kind;
792 1031 : const unsigned res_type_curr = myfrag->res_type_curr;
793 1031 : const unsigned res_type_prev = myfrag->res_type_prev;
794 1031 : const unsigned res_type_next = myfrag->res_type_next;
795 1031 : const unsigned needed_atoms = atleastned+myfrag->box_nb.size();
796 :
797 : /* Extra Distances are the same for each residue */
798 1032 : const unsigned xdsize=myfrag->xd1.size();
799 1032 : vector<Vector> ext_distances(xdsize);
800 2064 : vector<double> ext_d(xdsize);
801 27862 : for(unsigned q=0; q<xdsize; q++) {
802 29638 : if(myfrag->xd1[q]==-1||myfrag->xd2[q]==-1) continue;
803 24024 : const Vector distance = delta(getPosition(myfrag->xd1[q]),getPosition(myfrag->xd2[q]));
804 24024 : ext_d[q] = distance.modulo();
805 24024 : ext_distances[q] = distance/ext_d[q];
806 : }
807 :
808 : // CYCLE OVER THE SIX BACKBONE CHEMICAL SHIFTS
809 7224 : for(unsigned at_kind=0; at_kind<6; at_kind++) {
810 6192 : if(atom[s][a].exp_cs[at_kind]!=0) {
811 : // Common constant and AATYPE
812 3534 : const double * CONSTAACURR = db.CONSTAACURR(aa_kind,at_kind);
813 3534 : const double * CONSTAANEXT = db.CONSTAANEXT(aa_kind,at_kind);
814 3534 : const double * CONSTAAPREV = db.CONSTAAPREV(aa_kind,at_kind);
815 :
816 7068 : double cs = CONSTAACURR[res_type_curr] +
817 7068 : CONSTAANEXT[res_type_next] +
818 7068 : CONSTAAPREV[res_type_prev];
819 : // this is the atom for which we are calculating the chemical shift
820 3534 : const unsigned ipos = myfrag->pos[at_kind];
821 :
822 3534 : vector<unsigned> list;
823 3533 : list.reserve(needed_atoms);
824 3534 : list.push_back(ipos);
825 7031 : vector<Vector> ff;
826 3534 : ff.reserve(needed_atoms);
827 3534 : ff.push_back(Vector(0,0,0));
828 :
829 : //PREV
830 3534 : const double * CONST_BB2_PREV = db.CONST_BB2_PREV(aa_kind,at_kind);
831 3534 : const unsigned presize = myfrag->prev.size();
832 21202 : for(unsigned q=0; q<presize; q++) {
833 17669 : const double cb2pq = CONST_BB2_PREV[q];
834 20477 : if(cb2pq==0.) continue;
835 14861 : const unsigned jpos = myfrag->prev[q];
836 14861 : list.push_back(jpos);
837 14862 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
838 14862 : const double d = distance.modulo();
839 14862 : const double fact = cb2pq/d;
840 :
841 14862 : cs += cb2pq*d;
842 14862 : const Vector der = fact*distance;
843 14862 : ff[0] += der;
844 14862 : ff.push_back(-der);
845 : }
846 :
847 : //CURR
848 3533 : const double * CONST_BB2_CURR = db.CONST_BB2_CURR(aa_kind,at_kind);
849 3533 : const unsigned cursize = myfrag->curr.size();
850 24735 : for(unsigned q=0; q<cursize; q++) {
851 21201 : const double cb2cq = CONST_BB2_CURR[q];
852 34035 : if(cb2cq==0.) continue;
853 14784 : const unsigned jpos = myfrag->curr[q];
854 14784 : if(ipos==jpos) continue;
855 14784 : list.push_back(jpos);
856 14782 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
857 14783 : const double d = distance.modulo();
858 14782 : const double fact = cb2cq/d;
859 :
860 14782 : cs += cb2cq*d;
861 14782 : const Vector der = fact*distance;
862 14784 : ff[0] += der;
863 14784 : ff.push_back(-der);
864 : }
865 :
866 : //NEXT
867 3534 : const double * CONST_BB2_NEXT = db.CONST_BB2_NEXT(aa_kind,at_kind);
868 3533 : const unsigned nexsize = myfrag->next.size();
869 21202 : for(unsigned q=0; q<nexsize; q++) {
870 17668 : const double cb2nq = CONST_BB2_NEXT[q];
871 18604 : if(cb2nq==0.) continue;
872 16732 : const unsigned jpos = myfrag->next[q];
873 16733 : list.push_back(jpos);
874 16734 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
875 16733 : const double d = distance.modulo();
876 16734 : const double fact = cb2nq/d;
877 :
878 16734 : cs += cb2nq*d;
879 16734 : const Vector der = fact*distance;
880 16734 : ff[0] += der;
881 16733 : ff.push_back(-der);
882 : }
883 :
884 : //SIDE CHAIN
885 3534 : const double * CONST_SC2 = db.CONST_SC2(aa_kind,at_kind,res_type_curr);
886 3534 : const unsigned sidsize = myfrag->side_chain.size();
887 36587 : for(unsigned q=0; q<sidsize; q++) {
888 33053 : const double cs2q = CONST_SC2[q];
889 70815 : if(cs2q==0.) continue;
890 14172 : const unsigned jpos = myfrag->side_chain[q];
891 14172 : if(ipos==jpos) continue;
892 14172 : list.push_back(jpos);
893 14172 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
894 14172 : const double d = distance.modulo();
895 14172 : const double fact = cs2q/d;
896 :
897 14172 : cs += cs2q*d;
898 14172 : const Vector der = fact*distance;
899 14172 : ff[0] += der;
900 14172 : ff.push_back(-der);
901 : }
902 :
903 : //EXTRA DIST
904 3534 : const double * CONST_XD = db.CONST_XD(aa_kind,at_kind);
905 95411 : for(unsigned q=0; q<xdsize; q++) {
906 91877 : const double cxdq = CONST_XD[q];
907 102467 : if(cxdq==0.) continue;
908 96167 : if(myfrag->xd1[q]==-1||myfrag->xd2[q]==-1) continue;
909 83393 : list.push_back(myfrag->xd2[q]);
910 83387 : list.push_back(myfrag->xd1[q]);
911 83370 : cs += cxdq*ext_d[q];
912 83381 : const Vector der = cxdq*ext_distances[q];
913 83391 : ff.push_back( der);
914 83392 : ff.push_back(-der);
915 : }
916 :
917 : //NON BOND
918 : {
919 3534 : const double * CONST_CO_SPHERE3 = db.CO_SPHERE(aa_kind,at_kind,0);
920 3534 : const double * CONST_CO_SPHERE = db.CO_SPHERE(aa_kind,at_kind,1);
921 3534 : const unsigned boxsize = myfrag->box_nb.size();
922 378308 : for(unsigned bat=0; bat<boxsize; bat++) {
923 374774 : const unsigned jpos = myfrag->box_nb[bat];
924 374801 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
925 374885 : const double d2 = distance.modulo2();
926 :
927 374747 : if(d2<cutOffDist2) {
928 65121 : double factor1 = sqrt(d2);
929 65121 : double dfactor1 = 1./factor1;
930 65121 : double factor3 = dfactor1*dfactor1*dfactor1;
931 65121 : double dfactor3 = -3.*factor3*dfactor1*dfactor1;
932 :
933 65121 : if(d2>cutOnDist2) {
934 61026 : const double af = cutOffDist2 - d2;
935 61026 : const double bf = cutOffDist2 - 3.*cutOnDist2 + 2.*d2;
936 61026 : const double cf = invswitch*af;
937 61026 : const double df = cf*af*bf;
938 61026 : factor1 *= df;
939 61026 : factor3 *= df;
940 :
941 61026 : const double d4 = d2*d2;
942 61026 : const double af1 = 15.*cutOnDist2*d2;
943 61026 : const double bf1 = -14.*d4;
944 61026 : const double cf1 = -3.*cutOffDist2*cutOnDist2 + cutOffDist2*d2;
945 61026 : const double df1 = af1+bf1+cf1;
946 61026 : dfactor1 *= cf*(cutOffDist4+df1);
947 :
948 61026 : const double af3 = +2.*cutOffDist2*cutOnDist2;
949 61026 : const double bf3 = d2*(cutOffDist2+cutOnDist2);
950 61026 : const double cf3 = -2.*d4;
951 61026 : const double df3 = (af3+bf3+cf3)*d2;
952 61026 : dfactor3 *= invswitch*(cutMixed+df3);
953 : }
954 :
955 65121 : const unsigned t = type[jpos];
956 65122 : cs += factor1*CONST_CO_SPHERE[t] + factor3*CONST_CO_SPHERE3[t] ;
957 65122 : const double fact = dfactor1*CONST_CO_SPHERE[t]+dfactor3*CONST_CO_SPHERE3[t];
958 65122 : const Vector der = fact*distance;
959 :
960 65122 : list.push_back(jpos);
961 65119 : ff[0] += der;
962 65121 : ff.push_back(-der);
963 : }
964 : }
965 : }
966 : //END NON BOND
967 :
968 : //RINGS
969 : {
970 3534 : const double *rc = db.CO_RING(aa_kind,at_kind);
971 3534 : const unsigned rsize = ringInfo.size();
972 74210 : for(unsigned i=0; i<rsize; i++) {
973 : // compute angle from ring middle point to current atom position
974 : // get distance vector from query atom to ring center and normal vector to ring plane
975 70676 : const Vector n = ringInfo[i].normVect;
976 70677 : const double nL = ringInfo[i].lengthNV;
977 70679 : const double inL2 = ringInfo[i].lengthN2;
978 :
979 70677 : const Vector d = delta(ringInfo[i].position, getPosition(ipos));
980 70676 : const double dL2 = d.modulo2();
981 70673 : const double dL = sqrt(dL2);
982 70673 : const double idL3 = 1./(dL2*dL);
983 :
984 70673 : const double dn = dotProduct(d,n);
985 70680 : const double dn2 = dn*dn;
986 70680 : const double dLnL = dL*nL;
987 70680 : const double dL_nL = dL/nL;
988 :
989 70680 : const double ang2 = dn2*inL2/dL2;
990 70680 : const double u = 1.-3.*ang2;
991 70680 : const double cc = rc[ringInfo[i].rtype];
992 :
993 70680 : cs += cc*u*idL3;
994 :
995 70680 : const double fUU = -6*dn*inL2;
996 70680 : const double fUQ = fUU/dL;
997 70680 : const Vector gradUQ = fUQ*(dL2*n - dn*d);
998 70680 : const Vector gradVQ = (3*dL*u)*d;
999 :
1000 70680 : const double fact = cc*idL3*idL3;
1001 70680 : ff[0] += fact*(gradUQ - gradVQ);
1002 :
1003 70679 : const double fU = fUU/nL;
1004 70679 : double OneOverN = 1./6.;
1005 74213 : if(ringInfo[i].numAtoms==5) OneOverN=1./3.;
1006 70680 : const Vector factor2 = OneOverN*n;
1007 70679 : const Vector factor4 = (OneOverN/dL_nL)*d;
1008 :
1009 70680 : const Vector gradV = -OneOverN*gradVQ;
1010 :
1011 70680 : if(ringInfo[i].numAtoms==6) {
1012 : // update forces on ring atoms
1013 469857 : for(unsigned at=0; at<6; at++) {
1014 402718 : const Vector ab = crossProduct(d,ringInfo[i].g[at]);
1015 402821 : const Vector c = crossProduct(n,ringInfo[i].g[at]);
1016 402813 : const Vector factor3 = 0.5*dL_nL*c;
1017 402852 : const Vector factor1 = 0.5*ab;
1018 402857 : const Vector gradU = fU*( dLnL*(factor1 - factor2) -dn*(factor3 - factor4) );
1019 402835 : ff.push_back(fact*(gradU - gradV));
1020 402759 : list.push_back(ringInfo[i].atom[at]);
1021 : }
1022 : } else {
1023 14136 : for(unsigned at=0; at<3; at++) {
1024 10602 : const Vector ab = crossProduct(d,ringInfo[i].g[at]);
1025 10602 : const Vector c = crossProduct(n,ringInfo[i].g[at]);
1026 10602 : const Vector factor3 = dL_nL*c;
1027 10602 : const Vector factor1 = ab;
1028 10602 : const Vector gradU = fU*( dLnL*(factor1 - factor2) -dn*(factor3 - factor4) );
1029 10602 : ff.push_back(fact*(gradU - gradV));
1030 : }
1031 3534 : list.push_back(ringInfo[i].atom[0]);
1032 3534 : list.push_back(ringInfo[i].atom[2]);
1033 3534 : list.push_back(ringInfo[i].atom[3]);
1034 : }
1035 : }
1036 : }
1037 : //END OF RINGS
1038 :
1039 : //DIHEDRAL ANGLES
1040 : {
1041 3534 : const double *CO_DA = db.CO_DA(aa_kind,at_kind);
1042 3534 : if(myfrag->phi.size()==4) {
1043 3534 : const double *PARS_DA = db.PARS_DA(aa_kind,at_kind,0);
1044 3534 : const double val1 = 3.*myfrag->t_phi+PARS_DA[3];
1045 3534 : const double val2 = myfrag->t_phi+PARS_DA[4];
1046 3534 : cs += CO_DA[0]*(PARS_DA[0]*cos(val1)+PARS_DA[1]*cos(val2)+PARS_DA[2]);
1047 3534 : const double fact = -CO_DA[0]*(+3.*PARS_DA[0]*sin(val1)+PARS_DA[1]*sin(val2));
1048 :
1049 3534 : ff.push_back(fact*myfrag->dd0[0]);
1050 3534 : ff.push_back(fact*myfrag->dd10[0]);
1051 3534 : ff.push_back(fact*myfrag->dd21[0]);
1052 3534 : ff.push_back(-fact*myfrag->dd2[0]);
1053 3534 : list.push_back(myfrag->phi[0]);
1054 3534 : list.push_back(myfrag->phi[1]);
1055 3534 : list.push_back(myfrag->phi[2]);
1056 3534 : list.push_back(myfrag->phi[3]);
1057 : }
1058 :
1059 3534 : if(myfrag->psi.size()==4) {
1060 3534 : const double *PARS_DA = db.PARS_DA(aa_kind,at_kind,1);
1061 3534 : const double val1 = 3.*myfrag->t_psi+PARS_DA[3];
1062 3534 : const double val2 = myfrag->t_psi+PARS_DA[4];
1063 3534 : cs += CO_DA[1]*(PARS_DA[0]*cos(val1)+PARS_DA[1]*cos(val2)+PARS_DA[2]);
1064 3534 : const double fact = -CO_DA[1]*(+3.*PARS_DA[0]*sin(val1)+PARS_DA[1]*sin(val2));
1065 :
1066 3534 : ff.push_back(fact*myfrag->dd0[1]);
1067 3534 : ff.push_back(fact*myfrag->dd10[1]);
1068 3534 : ff.push_back(fact*myfrag->dd21[1]);
1069 3533 : ff.push_back(-fact*myfrag->dd2[1]);
1070 3534 : list.push_back(myfrag->psi[0]);
1071 3533 : list.push_back(myfrag->psi[1]);
1072 3534 : list.push_back(myfrag->psi[2]);
1073 3534 : list.push_back(myfrag->psi[3]);
1074 : }
1075 :
1076 : //Chi
1077 3534 : if(myfrag->chi1.size()==4) {
1078 2819 : const double *PARS_DA = db.PARS_DA(aa_kind,at_kind,2);
1079 2819 : const double val1 = 3.*myfrag->t_chi1+PARS_DA[3];
1080 2819 : const double val2 = myfrag->t_chi1+PARS_DA[4];
1081 2819 : cs += CO_DA[2]*(PARS_DA[0]*cos(val1)+PARS_DA[1]*cos(val2)+PARS_DA[2]);
1082 2819 : const double fact = -CO_DA[2]*(+3.*PARS_DA[0]*sin(val1)+PARS_DA[1]*sin(val2));
1083 :
1084 2819 : ff.push_back(fact*myfrag->dd0[2]);
1085 2820 : ff.push_back(fact*myfrag->dd10[2]);
1086 2820 : ff.push_back(fact*myfrag->dd21[2]);
1087 2820 : ff.push_back(-fact*myfrag->dd2[2]);
1088 2820 : list.push_back(myfrag->chi1[0]);
1089 2820 : list.push_back(myfrag->chi1[1]);
1090 2820 : list.push_back(myfrag->chi1[2]);
1091 2820 : list.push_back(myfrag->chi1[3]);
1092 : }
1093 : }
1094 : //END OF DIHE
1095 :
1096 : Value * comp;
1097 3533 : double fact = 1.0;
1098 3533 : if(!camshift) {
1099 2944 : comp = atom[s][a].comp[at_kind];
1100 2944 : comp->set(cs);
1101 2944 : Tensor virial;
1102 627063 : for(unsigned i=0; i<list.size(); i++) {
1103 624094 : setAtomsDerivatives(comp,list[i],fact*ff[i]);
1104 623823 : virial-=Tensor(getPosition(list[i]),fact*ff[i]);
1105 : }
1106 2945 : setBoxDerivatives(comp,virial);
1107 : } else {
1108 : // but I would also divide for the weights derived with metainference
1109 589 : comp = getPntrToValue();
1110 588 : score += (cs - atom[s][a].exp_cs[at_kind])*(cs - atom[s][a].exp_cs[at_kind])/camshift_sigma2[at_kind];
1111 589 : fact = 2.0*(cs - atom[s][a].exp_cs[at_kind])/camshift_sigma2[at_kind];
1112 589 : for(unsigned i=0; i<list.size(); i++) omp_deriv[list[i]] += fact*ff[i];
1113 3534 : }
1114 : }
1115 1032 : }
1116 : }
1117 48 : #pragma omp critical
1118 28 : if(camshift) for(int i=0; i<getPositions().size(); i++) setAtomsDerivatives(i,omp_deriv[i]);
1119 24 : }
1120 :
1121 : // in the case of camshift we calculate the virial at the end
1122 6 : if(camshift) {
1123 1 : setBoxDerivativesNoPbc();
1124 1 : setValue(score);
1125 : }
1126 :
1127 6 : ++box_count;
1128 6 : if(box_count == box_nupdate) box_count = 0;
1129 6 : }
1130 :
1131 6 : void CS2Backbone::update_neighb() {
1132 6 : const unsigned chainsize = atom.size();
1133 18 : for(unsigned s=0; s<chainsize; s++) {
1134 12 : const unsigned psize = atom[s].size();
1135 36 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1136 24 : for(unsigned a=1; a<psize-1; a++) {
1137 1031 : const unsigned boxsize = getNumberOfAtoms();
1138 1031 : atom[s][a].box_nb.clear();
1139 1031 : atom[s][a].box_nb.reserve(300);
1140 1032 : const unsigned res_curr = res_num[atom[s][a].pos[0]];
1141 2662960 : for(unsigned bat=0; bat<boxsize; bat++) {
1142 2661928 : const unsigned res_dist = abs(static_cast<int>(res_curr-res_num[bat]));
1143 2694210 : if(res_dist<2) continue;
1144 17643363 : for(unsigned at_kind=0; at_kind<6; at_kind++) {
1145 21530664 : if(atom[s][a].exp_cs[at_kind]==0.) continue;
1146 8570568 : const unsigned ipos = atom[s][a].pos[at_kind];
1147 8559910 : const Vector distance = delta(getPosition(bat),getPosition(ipos));
1148 8547849 : const double d2=distance.modulo2();
1149 8560869 : if(d2<cutOffNB2) {
1150 103367 : atom[s][a].box_nb.push_back(bat);
1151 103358 : break;
1152 : }
1153 : }
1154 : }
1155 : }
1156 : }
1157 6 : }
1158 :
1159 6 : void CS2Backbone::compute_ring_parameters() {
1160 126 : for(unsigned i=0; i<ringInfo.size(); i++) {
1161 120 : const unsigned size = ringInfo[i].numAtoms;
1162 120 : if(size==6) {
1163 114 : ringInfo[i].g[0] = delta(getPosition(ringInfo[i].atom[4]),getPosition(ringInfo[i].atom[2]));
1164 114 : ringInfo[i].g[1] = delta(getPosition(ringInfo[i].atom[5]),getPosition(ringInfo[i].atom[3]));
1165 114 : ringInfo[i].g[2] = delta(getPosition(ringInfo[i].atom[0]),getPosition(ringInfo[i].atom[4]));
1166 114 : ringInfo[i].g[3] = delta(getPosition(ringInfo[i].atom[1]),getPosition(ringInfo[i].atom[5]));
1167 114 : ringInfo[i].g[4] = delta(getPosition(ringInfo[i].atom[2]),getPosition(ringInfo[i].atom[0]));
1168 114 : ringInfo[i].g[5] = delta(getPosition(ringInfo[i].atom[3]),getPosition(ringInfo[i].atom[1]));
1169 114 : vector<Vector> a(size);
1170 114 : a[0] = getPosition(ringInfo[i].atom[0]);
1171 : // ring center
1172 114 : Vector midP = a[0];
1173 684 : for(unsigned j=1; j<size; j++) {
1174 570 : a[j] = getPosition(ringInfo[i].atom[j]);
1175 570 : midP += a[j];
1176 : }
1177 114 : midP /= (double) size;
1178 114 : ringInfo[i].position = midP;
1179 : // compute normal vector to plane containing first three atoms in array
1180 114 : ringInfo[i].n1 = crossProduct(delta(a[0],a[4]), delta(a[0],a[2]));
1181 : // compute normal vector to plane containing last three atoms in array
1182 : // NB: third atom of five-membered ring used for both computations above
1183 114 : ringInfo[i].n2 = crossProduct(delta(a[3],a[1]), delta(a[3],a[5]));
1184 : // ring plane normal vector is average of n1 and n2
1185 114 : ringInfo[i].normVect = 0.5*(ringInfo[i].n1 + ringInfo[i].n2);
1186 : } else {
1187 6 : ringInfo[i].g[0] = delta(getPosition(ringInfo[i].atom[3]),getPosition(ringInfo[i].atom[2]));
1188 6 : ringInfo[i].g[1] = delta(getPosition(ringInfo[i].atom[0]),getPosition(ringInfo[i].atom[3]));
1189 6 : ringInfo[i].g[2] = delta(getPosition(ringInfo[i].atom[2]),getPosition(ringInfo[i].atom[0]));
1190 6 : vector<Vector> a(size);
1191 36 : for(unsigned j=0; j<size; j++) {
1192 30 : a[j] = getPosition(ringInfo[i].atom[j]);
1193 : }
1194 : // ring center
1195 6 : Vector midP = (a[0]+a[2]+a[3])/3.;
1196 6 : ringInfo[i].position = midP;
1197 : // compute normal vector to plane containing first three atoms in array
1198 6 : ringInfo[i].n1 = crossProduct(delta(a[0],a[3]), delta(a[0],a[2]));
1199 : // ring plane normal vector is average of n1 and n2
1200 6 : ringInfo[i].normVect = ringInfo[i].n1;
1201 :
1202 : }
1203 : // calculate squared length and length of normal vector
1204 120 : ringInfo[i].lengthN2 = 1./ringInfo[i].normVect.modulo2();
1205 120 : ringInfo[i].lengthNV = 1./sqrt(ringInfo[i].lengthN2);
1206 : }
1207 6 : }
1208 :
1209 6 : void CS2Backbone::compute_dihedrals() {
1210 6 : const unsigned chainsize = atom.size();
1211 18 : for(unsigned s=0; s<chainsize; s++) {
1212 12 : const unsigned psize = atom[s].size();
1213 36 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1214 24 : for(unsigned a=1; a<psize-1; a++) {
1215 1030 : const Fragment *myfrag = &atom[s][a];
1216 1031 : if(myfrag->phi.size()==4) {
1217 1030 : const Vector d0 = delta(getPosition(myfrag->phi[1]), getPosition(myfrag->phi[0]));
1218 1032 : const Vector d1 = delta(getPosition(myfrag->phi[2]), getPosition(myfrag->phi[1]));
1219 1032 : const Vector d2 = delta(getPosition(myfrag->phi[3]), getPosition(myfrag->phi[2]));
1220 : Torsion t;
1221 1032 : Vector dd0, dd1, dd2;
1222 1032 : atom[s][a].t_phi = t.compute(d0,d1,d2,dd0,dd1,dd2);
1223 1032 : atom[s][a].dd0[0] = dd0;
1224 1032 : atom[s][a].dd10[0] = dd1-dd0;
1225 1029 : atom[s][a].dd21[0] = dd2-dd1;
1226 1029 : atom[s][a].dd2[0] = dd2;
1227 : }
1228 1026 : if(myfrag->psi.size()==4) {
1229 1026 : const Vector d0 = delta(getPosition(myfrag->psi[1]), getPosition(myfrag->psi[0]));
1230 1032 : const Vector d1 = delta(getPosition(myfrag->psi[2]), getPosition(myfrag->psi[1]));
1231 1032 : const Vector d2 = delta(getPosition(myfrag->psi[3]), getPosition(myfrag->psi[2]));
1232 : Torsion t;
1233 1032 : Vector dd0, dd1, dd2;
1234 1032 : atom[s][a].t_psi = t.compute(d0,d1,d2,dd0,dd1,dd2);
1235 1032 : atom[s][a].dd0[1] = dd0;
1236 1032 : atom[s][a].dd10[1] = dd1-dd0;
1237 1032 : atom[s][a].dd21[1] = dd2-dd1;
1238 1032 : atom[s][a].dd2[1] = dd2;
1239 : }
1240 1032 : if(myfrag->chi1.size()==4) {
1241 796 : const Vector d0 = delta(getPosition(myfrag->chi1[1]), getPosition(myfrag->chi1[0]));
1242 798 : const Vector d1 = delta(getPosition(myfrag->chi1[2]), getPosition(myfrag->chi1[1]));
1243 798 : const Vector d2 = delta(getPosition(myfrag->chi1[3]), getPosition(myfrag->chi1[2]));
1244 : Torsion t;
1245 798 : Vector dd0, dd1, dd2;
1246 798 : atom[s][a].t_chi1 = t.compute(d0,d1,d2,dd0,dd1,dd2);
1247 798 : atom[s][a].dd0[2] = dd0;
1248 798 : atom[s][a].dd10[2] = dd1-dd0;
1249 798 : atom[s][a].dd21[2] = dd2-dd1;
1250 798 : atom[s][a].dd2[2] = dd2;
1251 : }
1252 : }
1253 : }
1254 6 : }
1255 :
1256 6 : void CS2Backbone::init_backbone(const PDB &pdb) {
1257 : // number of chains
1258 6 : vector<string> chains;
1259 6 : pdb.getChainNames( chains );
1260 6 : seg_last.resize(chains.size());
1261 6 : unsigned old_size=0;
1262 :
1263 18 : for(unsigned i=0; i<chains.size(); i++) {
1264 : unsigned start, end;
1265 12 : string errmsg;
1266 12 : pdb.getResidueRange( chains[i], start, end, errmsg );
1267 :
1268 12 : unsigned res_offset = start;
1269 12 : unsigned resrange = end-start+1;
1270 :
1271 12 : if(i==0) seg_last[i] = resrange;
1272 6 : else seg_last[i] = seg_last[i-1]+resrange;
1273 :
1274 24 : vector<int> N_;
1275 24 : vector<int> H_;
1276 24 : vector<int> CA_;
1277 24 : vector<int> CB_;
1278 24 : vector<int> HA_;
1279 24 : vector<int> C_;
1280 24 : vector<int> O_;
1281 24 : vector<int> CX_;
1282 12 : N_.resize (resrange,-1);
1283 12 : H_.resize (resrange,-1);
1284 12 : CA_.resize(resrange,-1);
1285 12 : CB_.resize(resrange,-1);
1286 12 : HA_.resize(resrange,-1);
1287 12 : C_.resize (resrange,-1);
1288 12 : O_.resize (resrange,-1);
1289 12 : CX_.resize(resrange,-1);
1290 :
1291 24 : vector<AtomNumber> allatoms = pdb.getAtomsInChain(chains[i]);
1292 : // cycle over all the atoms in the chain
1293 15684 : for(unsigned a=0; a<allatoms.size(); a++) {
1294 15672 : unsigned atm_index=a+old_size;
1295 15672 : unsigned f = pdb.getResidueNumber(allatoms[a]);
1296 15672 : unsigned f_idx = f-res_offset;
1297 15672 : string AN = pdb.getAtomName(allatoms[a]);
1298 31344 : string RES = pdb.getResidueName(allatoms[a]);
1299 15672 : if(AN=="N") N_ [f_idx] = atm_index;
1300 14616 : else if(AN=="H" ||AN=="HN" ) H_ [f_idx] = atm_index;
1301 13614 : else if(AN=="HA"||AN=="HA1") HA_[f_idx] = atm_index;
1302 12558 : else if(AN=="CA" ) CA_[f_idx] = atm_index;
1303 11502 : else if(AN=="CB" ) CB_[f_idx] = atm_index;
1304 10626 : else if(AN=="C" ) C_ [f_idx] = atm_index;
1305 9570 : else if(AN=="O" ) O_ [f_idx] = atm_index;
1306 8526 : else if(AN=="CD"&&RES=="PRO") H_[f_idx] = atm_index;
1307 15672 : if(is_chi1_cx(RES,AN)) CX_[f_idx] = atm_index;
1308 15672 : }
1309 12 : old_size+=allatoms.size();
1310 :
1311 : // vector of residues for a given chain
1312 24 : vector<Fragment> atm_;
1313 : // cycle over all residues in the chain
1314 1068 : for(unsigned a=start; a<=end; a++) {
1315 1056 : unsigned f_idx = a - res_offset;
1316 1056 : Fragment at;
1317 1056 : at.pos[0] = HA_[f_idx];
1318 1056 : at.pos[1] = H_[f_idx];
1319 1056 : at.pos[2] = N_[f_idx];
1320 1056 : at.pos[3] = CA_[f_idx];
1321 1056 : at.pos[4] = CB_[f_idx];
1322 1056 : at.pos[5] = C_[f_idx];
1323 1056 : at.res_type_prev = at.res_type_curr = at.res_type_next = 0;
1324 1056 : at.res_name = pdb.getResidueName(a, chains[i]);
1325 1056 : at.res_kind = db.kind(at.res_name);
1326 1056 : at.fd = a;
1327 : //REGISTER PREV CURR NEXT
1328 : {
1329 1056 : if(a>start) {
1330 1044 : at.prev.push_back( N_[f_idx-1]);
1331 1044 : at.prev.push_back(CA_[f_idx-1]);
1332 1044 : at.prev.push_back(HA_[f_idx-1]);
1333 1044 : at.prev.push_back( C_[f_idx-1]);
1334 1044 : at.prev.push_back( O_[f_idx-1]);
1335 1044 : at.res_type_prev = frag2enum(pdb.getResidueName(a-1, chains[i]));
1336 : }
1337 :
1338 1056 : at.curr.push_back( N_[f_idx]);
1339 1056 : at.curr.push_back( H_[f_idx]);
1340 1056 : at.curr.push_back(CA_[f_idx]);
1341 1056 : at.curr.push_back(HA_[f_idx]);
1342 1056 : at.curr.push_back( C_[f_idx]);
1343 1056 : at.curr.push_back( O_[f_idx]);
1344 1056 : at.res_type_curr = frag2enum(pdb.getResidueName(a, chains[i]));
1345 :
1346 1056 : if(a<end) {
1347 1044 : at.next.push_back (N_[f_idx+1]);
1348 1044 : at.next.push_back (H_[f_idx+1]);
1349 1044 : at.next.push_back(CA_[f_idx+1]);
1350 1044 : at.next.push_back(HA_[f_idx+1]);
1351 1044 : at.next.push_back( C_[f_idx+1]);
1352 1044 : at.res_type_next = frag2enum(pdb.getResidueName(a+1, chains[i]));
1353 : }
1354 :
1355 : //PHI | PSI | CH1
1356 1056 : if(a>start) {
1357 1044 : at.phi.push_back( C_[f_idx-1]);
1358 1044 : at.phi.push_back( N_[f_idx]);
1359 1044 : at.phi.push_back(CA_[f_idx]);
1360 1044 : at.phi.push_back( C_[f_idx]);
1361 : }
1362 :
1363 1056 : if(a<end) {
1364 1044 : at.psi.push_back( N_[f_idx]);
1365 1044 : at.psi.push_back(CA_[f_idx]);
1366 1044 : at.psi.push_back( C_[f_idx]);
1367 1044 : at.psi.push_back( N_[f_idx+1]);
1368 : }
1369 :
1370 1056 : if(CX_[f_idx]!=-1&&CB_[f_idx]!=-1) {
1371 816 : at.chi1.push_back( N_[f_idx]);
1372 816 : at.chi1.push_back(CA_[f_idx]);
1373 816 : at.chi1.push_back(CB_[f_idx]);
1374 816 : at.chi1.push_back(CX_[f_idx]);
1375 : }
1376 : }
1377 1056 : atm_.push_back(at);
1378 1056 : }
1379 12 : atom.push_back(atm_);
1380 18 : }
1381 6 : }
1382 :
1383 6 : void CS2Backbone::init_sidechain(const PDB &pdb) {
1384 6 : vector<string> chains;
1385 6 : pdb.getChainNames( chains );
1386 6 : unsigned old_size=0;
1387 : // cycle over chains
1388 18 : for(unsigned s=0; s<atom.size(); s++) {
1389 12 : AtomNumber astart, aend;
1390 12 : string errmsg;
1391 12 : pdb.getAtomRange( chains[s], astart, aend, errmsg );
1392 12 : unsigned atom_offset = astart.index();
1393 : // cycle over residues
1394 1068 : for(unsigned a=0; a<atom[s].size(); a++) {
1395 1056 : if(atom[s][a].res_name=="UNK") continue;
1396 1056 : vector<AtomNumber> atm = pdb.getAtomsInResidue(atom[s][a].fd, chains[s]);
1397 2112 : vector<string> sc_atm = side_chain_atoms(atom[s][a].res_name);
1398 :
1399 13338 : for(unsigned sc=0; sc<sc_atm.size(); sc++) {
1400 231804 : for(unsigned aa=0; aa<atm.size(); aa++) {
1401 219522 : if(pdb.getAtomName(atm[aa])==sc_atm[sc]) {
1402 9342 : atom[s][a].side_chain.push_back(atm[aa].index()-atom_offset+old_size);
1403 : }
1404 : }
1405 : }
1406 :
1407 1056 : }
1408 12 : old_size = aend.index()+1;
1409 18 : }
1410 6 : }
1411 :
1412 6 : void CS2Backbone::init_xdist(const PDB &pdb) {
1413 : const string atomsP1[] = {"H", "H", "H", "C", "C", "C",
1414 : "O", "O", "O", "N", "N", "N",
1415 : "O", "O", "O", "N", "N", "N",
1416 : "CG", "CG", "CG", "CG", "CG",
1417 : "CG", "CG", "CA"
1418 162 : };
1419 :
1420 : const int resOffsetP1 [] = {0, 0, 0, -1, -1, -1,
1421 : 0, 0, 0, 1, 1, 1,
1422 : -1, -1, -1, 0, 0, 0,
1423 : 0, 0, 0, 0, 0, -1, 1, -1
1424 6 : };
1425 :
1426 : const string atomsP2[] = {"HA", "C", "CB", "HA", "C", "CB",
1427 : "HA", "N", "CB", "HA", "N", "CB",
1428 : "HA", "N", "CB", "HA", "N", "CB",
1429 : "HA", "N", "C", "C", "N", "CA", "CA", "CA"
1430 12 : };
1431 :
1432 6 : const int resOffsetP2 [] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -1, 0, 0, 0, -1, 1, 0, 0, 1};
1433 :
1434 12 : vector<string> chains;
1435 6 : pdb.getChainNames( chains );
1436 6 : unsigned old_size=0;
1437 18 : for(unsigned s=0; s<atom.size(); s++) {
1438 12 : AtomNumber astart, aend;
1439 12 : string errmsg;
1440 12 : pdb.getAtomRange( chains[s], astart, aend, errmsg );
1441 12 : unsigned atom_offset = astart.index();
1442 :
1443 1044 : for(unsigned a=1; a<atom[s].size()-1; a++) {
1444 1032 : vector<AtomNumber> atm_curr = pdb.getAtomsInResidue(atom[s][a].fd,chains[s]);
1445 2064 : vector<AtomNumber> atm_prev = pdb.getAtomsInResidue(atom[s][a].fd-1,chains[s]);
1446 2064 : vector<AtomNumber> atm_next = pdb.getAtomsInResidue(atom[s][a].fd+1,chains[s]);
1447 :
1448 27864 : for(unsigned q=0; q<db.get_numXtraDists()-1; q++) {
1449 26832 : vector<AtomNumber>::iterator at1, at1_end;
1450 26832 : vector<AtomNumber>::iterator at2, at2_end;
1451 :
1452 26832 : bool init_p1=false;
1453 26832 : AtomNumber p1;
1454 26832 : bool init_p2=false;
1455 26832 : AtomNumber p2;
1456 :
1457 26832 : if(resOffsetP1[q]== 0) { at1 = atm_curr.begin(); at1_end = atm_curr.end();}
1458 26832 : if(resOffsetP1[q]==-1) { at1 = atm_prev.begin(); at1_end = atm_prev.end();}
1459 26832 : if(resOffsetP1[q]==+1) { at1 = atm_next.begin(); at1_end = atm_next.end();}
1460 239766 : while(at1!=at1_end) {
1461 211170 : AtomNumber aa = *at1;
1462 211170 : ++at1;
1463 211170 : string name = pdb.getAtomName(aa);
1464 :
1465 211170 : xdist_name_map(name);
1466 :
1467 211170 : if(name==atomsP1[q]) {
1468 25068 : p1 = aa;
1469 25068 : init_p1=true;
1470 25068 : break;
1471 : }
1472 186102 : }
1473 :
1474 26832 : if(resOffsetP2[q]== 0) { at2 = atm_curr.begin(); at2_end = atm_curr.end();}
1475 26832 : if(resOffsetP2[q]==-1) { at2 = atm_prev.begin(); at2_end = atm_prev.end();}
1476 26832 : if(resOffsetP2[q]==+1) { at2 = atm_next.begin(); at2_end = atm_next.end();}
1477 166716 : while(at2!=at2_end) {
1478 138840 : AtomNumber aa = *at2;
1479 138840 : ++at2;
1480 138840 : string name = pdb.getAtomName(aa);
1481 :
1482 138840 : xdist_name_map(name);
1483 :
1484 138840 : if(name==atomsP2[q]) {
1485 25788 : p2 = aa;
1486 25788 : init_p2=true;
1487 25788 : break;
1488 : }
1489 113052 : }
1490 26832 : int add1 = -1;
1491 26832 : int add2 = -1;
1492 26832 : if(init_p1) add1=p1.index()-atom_offset+old_size;
1493 26832 : if(init_p2) add2=p2.index()-atom_offset+old_size;
1494 26832 : atom[s][a].xd1.push_back(add1);
1495 26832 : atom[s][a].xd2.push_back(add2);
1496 : }
1497 1032 : }
1498 12 : old_size = aend.index()+1;
1499 174 : }
1500 6 : }
1501 :
1502 6 : void CS2Backbone::init_types(const PDB &pdb) {
1503 6 : vector<AtomNumber> aa = pdb.getAtomNumbers();
1504 15678 : for(unsigned i=0; i<aa.size(); i++) {
1505 15672 : unsigned frag = pdb.getResidueNumber(aa[i]);
1506 15672 : string fragName = pdb.getResidueName(aa[i]);
1507 31344 : string atom_name = pdb.getAtomName(aa[i]);
1508 15672 : char atom_type = atom_name[0];
1509 15672 : if(isdigit(atom_name[0])) atom_type = atom_name[1];
1510 15672 : res_num.push_back(frag);
1511 15672 : unsigned t = 0;
1512 15672 : if (!isSP2(fragName, atom_name)) {
1513 12312 : if (atom_type == 'C') t = D_C;
1514 9306 : else if (atom_type == 'O') t = D_O;
1515 9150 : else if (atom_type == 'H') t = D_H;
1516 1398 : else if (atom_type == 'N') t = D_N;
1517 54 : else if (atom_type == 'S') t = D_S;
1518 0 : else plumed_merror("CS2Backbone:init_type: unknown atom type!\n");
1519 : } else {
1520 3360 : if (atom_type == 'C') t = D_C2;
1521 1368 : else if (atom_type == 'O') t = D_O2;
1522 0 : else if (atom_type == 'N') t = D_N2;
1523 0 : else plumed_merror("CS2Backbone:init_type: unknown atom type!\n");
1524 : }
1525 15672 : type.push_back(t);
1526 15678 : }
1527 6 : }
1528 :
1529 6 : void CS2Backbone::init_rings(const PDB &pdb) {
1530 :
1531 42 : const string pheTyr_n[] = {"CG","CD1","CE1","CZ","CE2","CD2"};
1532 12 : const string trp1_n[] = {"CD2","CE2","CZ2","CH2","CZ3","CE3"};
1533 12 : const string trp2_n[] = {"CG","CD1","NE1","CE2","CD2"};
1534 12 : const string his_n[] = {"CG","ND1","CD2","CE1","NE2"};
1535 :
1536 12 : vector<string> chains;
1537 6 : pdb.getChainNames( chains );
1538 12 : vector<AtomNumber> allatoms = pdb.getAtomNumbers();
1539 6 : unsigned old_size=0;
1540 :
1541 18 : for(unsigned s=0; s<atom.size(); s++) {
1542 12 : AtomNumber astart, aend;
1543 12 : string errmsg;
1544 12 : pdb.getAtomRange( chains[s], astart, aend, errmsg );
1545 12 : unsigned atom_offset = astart.index();
1546 1068 : for(unsigned r=0; r<atom[s].size(); r++) {
1547 1056 : string frg = pdb.getResidueName(atom[s][r].fd);
1548 4014 : if(!((frg=="PHE")||(frg=="TYR")||(frg=="TRP")||
1549 2826 : (frg=="HIS")||(frg=="HIP")||(frg=="HID")||
1550 1884 : (frg=="HIE")||(frg=="HSD")||(frg=="HSE")||
1551 2940 : (frg=="HSP"))) continue;
1552 :
1553 228 : vector<AtomNumber> frg_atoms = pdb.getAtomsInResidue(atom[s][r].fd,chains[s]);
1554 :
1555 114 : if(frg=="PHE"||frg=="TYR") {
1556 108 : RingInfo ri;
1557 2280 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1558 2172 : unsigned atm = frg_atoms[a].index()-atom_offset+old_size;
1559 12936 : for(unsigned aa=0; aa<6; aa++) {
1560 11412 : if(pdb.getAtomName(frg_atoms[a])==pheTyr_n[aa]) {
1561 648 : ri.atom[aa] = atm;
1562 648 : break;
1563 : }
1564 : }
1565 : }
1566 108 : ri.numAtoms = 6;
1567 108 : if(frg=="PHE") ri.rtype = RingInfo::R_PHE;
1568 108 : if(frg=="TYR") ri.rtype = RingInfo::R_TYR;
1569 108 : ringInfo.push_back(ri);
1570 :
1571 6 : } else if(frg=="TRP") {
1572 : //First ring
1573 6 : RingInfo ri;
1574 150 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1575 144 : unsigned atm = frg_atoms[a].index()-atom_offset+old_size;
1576 882 : for(unsigned aa=0; aa<6; aa++) {
1577 774 : if(pdb.getAtomName(frg_atoms[a])==trp1_n[aa]) {
1578 36 : ri.atom[aa] = atm;
1579 36 : break;
1580 : }
1581 : }
1582 : }
1583 6 : ri.numAtoms = 6;
1584 6 : ri.rtype = RingInfo::R_TRP1;
1585 6 : ringInfo.push_back(ri);
1586 : //Second Ring
1587 6 : RingInfo ri2;
1588 150 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1589 144 : unsigned atm = frg_atoms[a].index()-atom_offset+old_size;
1590 774 : for(unsigned aa=0; aa<5; aa++) {
1591 660 : if(pdb.getAtomName(frg_atoms[a])==trp2_n[aa]) {
1592 30 : ri2.atom[aa] = atm;
1593 30 : break;
1594 : }
1595 : }
1596 : }
1597 6 : ri2.numAtoms = 5;
1598 6 : ri2.rtype = RingInfo::R_TRP2;
1599 6 : ringInfo.push_back(ri2);
1600 :
1601 0 : } else if((frg=="HIS")||(frg=="HIP")||(frg=="HID")||
1602 0 : (frg=="HIE")||(frg=="HSD")||(frg=="HSE")||
1603 0 : (frg=="HSP")) {//HIS case
1604 0 : RingInfo ri;
1605 0 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1606 0 : unsigned atm = frg_atoms[a].index()-atom_offset+old_size;
1607 0 : for(unsigned aa=0; aa<5; aa++) {
1608 0 : if(pdb.getAtomName(frg_atoms[a])==his_n[aa]) {
1609 0 : ri.atom[aa] = atm;
1610 0 : break;
1611 : }
1612 : }
1613 : }
1614 0 : ri.numAtoms = 5;
1615 0 : ri.rtype = RingInfo::R_HIS;
1616 0 : ringInfo.push_back(ri);
1617 : } else {
1618 0 : plumed_merror("Unkwown Ring Fragment");
1619 : }
1620 114 : }
1621 12 : old_size = aend.index()+1;
1622 54 : }
1623 6 : }
1624 :
1625 3144 : CS2Backbone::aa_t CS2Backbone::frag2enum(const string &aa) {
1626 :
1627 3144 : aa_t type = ALA;
1628 3144 : if (aa == "ALA") type = ALA;
1629 2964 : else if (aa == "ARG") type = ARG;
1630 2838 : else if (aa == "ASN") type = ASN;
1631 2676 : else if (aa == "ASP") type = ASP;
1632 2520 : else if (aa == "ASH") type = ASP;
1633 2520 : else if (aa == "CYS") type = CYS;
1634 2448 : else if (aa == "CYM") type = CYS;
1635 2448 : else if (aa == "GLN") type = GLN;
1636 2394 : else if (aa == "GLU") type = GLU;
1637 2184 : else if (aa == "GLH") type = GLU;
1638 2184 : else if (aa == "GLY") type = GLY;
1639 1650 : else if (aa == "HIS") type = HIS;
1640 1650 : else if (aa == "HSE") type = HIS;
1641 1650 : else if (aa == "HIE") type = HIS;
1642 1650 : else if (aa == "HSP") type = HIS;
1643 1650 : else if (aa == "HIP") type = HIS;
1644 1650 : else if (aa == "HSD") type = HIS;
1645 1650 : else if (aa == "HID") type = HIS;
1646 1650 : else if (aa == "ILE") type = ILE;
1647 1470 : else if (aa == "LEU") type = LEU;
1648 1308 : else if (aa == "LYS") type = LYS;
1649 1056 : else if (aa == "MET") type = MET;
1650 972 : else if (aa == "PHE") type = PHE;
1651 684 : else if (aa == "PRO") type = PRO;
1652 558 : else if (aa == "SER") type = SER;
1653 396 : else if (aa == "THR") type = THR;
1654 198 : else if (aa == "TRP") type = TRP;
1655 180 : else if (aa == "TYR") type = TYR;
1656 144 : else if (aa == "VAL") type = VAL;
1657 0 : else if (aa == "UNK") type = UNK;
1658 0 : else plumed_merror("CS2Backbone: Error converting string " + aa + " into amino acid index: not a valid 3-letter code");
1659 3144 : return type;
1660 : }
1661 :
1662 1056 : vector<string> CS2Backbone::side_chain_atoms(const string &s) {
1663 1056 : vector<string> sc;
1664 :
1665 1056 : if(s=="ALA") {
1666 60 : sc.push_back( "CB" );
1667 60 : sc.push_back( "HB1" );
1668 60 : sc.push_back( "HB2" );
1669 60 : sc.push_back( "HB3" );
1670 60 : return sc;
1671 996 : } else if(s=="ARG") {
1672 42 : sc.push_back( "CB" );
1673 42 : sc.push_back( "CG" );
1674 42 : sc.push_back( "CD" );
1675 42 : sc.push_back( "NE" );
1676 42 : sc.push_back( "CZ" );
1677 42 : sc.push_back( "NH1" );
1678 42 : sc.push_back( "NH2" );
1679 42 : sc.push_back( "NH3" );
1680 42 : sc.push_back( "HB1" );
1681 42 : sc.push_back( "HB2" );
1682 42 : sc.push_back( "HB3" );
1683 42 : sc.push_back( "HG1" );
1684 42 : sc.push_back( "HG2" );
1685 42 : sc.push_back( "HG3" );
1686 42 : sc.push_back( "HD1" );
1687 42 : sc.push_back( "HD2" );
1688 42 : sc.push_back( "HD3" );
1689 42 : sc.push_back( "HE" );
1690 42 : sc.push_back( "HH11" );
1691 42 : sc.push_back( "HH12" );
1692 42 : sc.push_back( "HH21" );
1693 42 : sc.push_back( "HH22" );
1694 42 : sc.push_back( "1HH1" );
1695 42 : sc.push_back( "2HH1" );
1696 42 : sc.push_back( "1HH2" );
1697 42 : sc.push_back( "2HH2" );
1698 42 : return sc;
1699 954 : } else if(s=="ASN") {
1700 54 : sc.push_back( "CB" );
1701 54 : sc.push_back( "CG" );
1702 54 : sc.push_back( "OD1" );
1703 54 : sc.push_back( "ND2" );
1704 54 : sc.push_back( "HB1" );
1705 54 : sc.push_back( "HB2" );
1706 54 : sc.push_back( "HB3" );
1707 54 : sc.push_back( "HD21" );
1708 54 : sc.push_back( "HD22" );
1709 54 : sc.push_back( "1HD2" );
1710 54 : sc.push_back( "2HD2" );
1711 54 : return sc;
1712 900 : } else if(s=="ASP"||s=="ASH") {
1713 54 : sc.push_back( "CB" );
1714 54 : sc.push_back( "CG" );
1715 54 : sc.push_back( "OD1" );
1716 54 : sc.push_back( "OD2" );
1717 54 : sc.push_back( "HB1" );
1718 54 : sc.push_back( "HB2" );
1719 54 : sc.push_back( "HB3" );
1720 54 : return sc;
1721 846 : } else if(s=="CYS"||s=="CYM") {
1722 24 : sc.push_back( "CB" );
1723 24 : sc.push_back( "SG" );
1724 24 : sc.push_back( "HB1" );
1725 24 : sc.push_back( "HB2" );
1726 24 : sc.push_back( "HB3" );
1727 24 : sc.push_back( "HG1" );
1728 24 : sc.push_back( "HG" );
1729 24 : return sc;
1730 822 : } else if(s=="GLN") {
1731 18 : sc.push_back( "CB" );
1732 18 : sc.push_back( "CG" );
1733 18 : sc.push_back( "CD" );
1734 18 : sc.push_back( "OE1" );
1735 18 : sc.push_back( "NE2" );
1736 18 : sc.push_back( "HB1" );
1737 18 : sc.push_back( "HB2" );
1738 18 : sc.push_back( "HB3" );
1739 18 : sc.push_back( "HG1" );
1740 18 : sc.push_back( "HG2" );
1741 18 : sc.push_back( "HG3" );
1742 18 : sc.push_back( "HE21" );
1743 18 : sc.push_back( "HE22" );
1744 18 : sc.push_back( "1HE2" );
1745 18 : sc.push_back( "2HE2" );
1746 18 : return sc;
1747 804 : } else if(s=="GLU"||s=="GLH") {
1748 72 : sc.push_back( "CB" );
1749 72 : sc.push_back( "CG" );
1750 72 : sc.push_back( "CD" );
1751 72 : sc.push_back( "OE1" );
1752 72 : sc.push_back( "OE2" );
1753 72 : sc.push_back( "HB1" );
1754 72 : sc.push_back( "HB2" );
1755 72 : sc.push_back( "HB3" );
1756 72 : sc.push_back( "HG1" );
1757 72 : sc.push_back( "HG2" );
1758 72 : sc.push_back( "HG3" );
1759 72 : return sc;
1760 732 : } else if(s=="GLY") {
1761 180 : sc.push_back( "HA2" );
1762 180 : return sc;
1763 552 : } else if(s=="HIS"||s=="HSE"||s=="HIE"||s=="HSD"||s=="HID"||s=="HIP"||s=="HSP") {
1764 0 : sc.push_back( "CB" );
1765 0 : sc.push_back( "CG" );
1766 0 : sc.push_back( "ND1" );
1767 0 : sc.push_back( "CD2" );
1768 0 : sc.push_back( "CE1" );
1769 0 : sc.push_back( "NE2" );
1770 0 : sc.push_back( "HB1" );
1771 0 : sc.push_back( "HB2" );
1772 0 : sc.push_back( "HB3" );
1773 0 : sc.push_back( "HD1" );
1774 0 : sc.push_back( "HD2" );
1775 0 : sc.push_back( "HE1" );
1776 0 : sc.push_back( "HE2" );
1777 0 : return sc;
1778 552 : } else if(s=="ILE") {
1779 60 : sc.push_back( "CB" );
1780 60 : sc.push_back( "CG1" );
1781 60 : sc.push_back( "CG2" );
1782 60 : sc.push_back( "CD" );
1783 60 : sc.push_back( "HB" );
1784 60 : sc.push_back( "HG11" );
1785 60 : sc.push_back( "HG12" );
1786 60 : sc.push_back( "HG21" );
1787 60 : sc.push_back( "HG22" );
1788 60 : sc.push_back( "HG23" );
1789 60 : sc.push_back( "1HG1" );
1790 60 : sc.push_back( "2HG1" );
1791 60 : sc.push_back( "1HG2" );
1792 60 : sc.push_back( "2HG2" );
1793 60 : sc.push_back( "3HG2" );
1794 60 : sc.push_back( "HD1" );
1795 60 : sc.push_back( "HD2" );
1796 60 : sc.push_back( "HD3" );
1797 60 : return sc;
1798 492 : } else if(s=="LEU") {
1799 54 : sc.push_back( "CB" );
1800 54 : sc.push_back( "CG" );
1801 54 : sc.push_back( "CD1" );
1802 54 : sc.push_back( "CD2" );
1803 54 : sc.push_back( "HB1" );
1804 54 : sc.push_back( "HB2" );
1805 54 : sc.push_back( "HB3" );
1806 54 : sc.push_back( "HG" );
1807 54 : sc.push_back( "HD11" );
1808 54 : sc.push_back( "HD12" );
1809 54 : sc.push_back( "HD13" );
1810 54 : sc.push_back( "HD21" );
1811 54 : sc.push_back( "HD22" );
1812 54 : sc.push_back( "HD23" );
1813 54 : sc.push_back( "1HD1" );
1814 54 : sc.push_back( "2HD1" );
1815 54 : sc.push_back( "3HD1" );
1816 54 : sc.push_back( "1HD2" );
1817 54 : sc.push_back( "2HD2" );
1818 54 : sc.push_back( "3HD2" );
1819 54 : return sc;
1820 438 : } else if(s=="LYS") {
1821 84 : sc.push_back( "CB" );
1822 84 : sc.push_back( "CG" );
1823 84 : sc.push_back( "CD" );
1824 84 : sc.push_back( "CE" );
1825 84 : sc.push_back( "NZ" );
1826 84 : sc.push_back( "HB1" );
1827 84 : sc.push_back( "HB2" );
1828 84 : sc.push_back( "HB3" );
1829 84 : sc.push_back( "HG1" );
1830 84 : sc.push_back( "HG2" );
1831 84 : sc.push_back( "HG3" );
1832 84 : sc.push_back( "HD1" );
1833 84 : sc.push_back( "HD2" );
1834 84 : sc.push_back( "HD3" );
1835 84 : sc.push_back( "HE1" );
1836 84 : sc.push_back( "HE2" );
1837 84 : sc.push_back( "HE3" );
1838 84 : sc.push_back( "HZ1" );
1839 84 : sc.push_back( "HZ2" );
1840 84 : sc.push_back( "HZ3" );
1841 84 : return sc;
1842 354 : } else if(s=="MET") {
1843 30 : sc.push_back( "CB" );
1844 30 : sc.push_back( "CG" );
1845 30 : sc.push_back( "SD" );
1846 30 : sc.push_back( "CE" );
1847 30 : sc.push_back( "HB1" );
1848 30 : sc.push_back( "HB2" );
1849 30 : sc.push_back( "HB3" );
1850 30 : sc.push_back( "HG1" );
1851 30 : sc.push_back( "HG2" );
1852 30 : sc.push_back( "HG3" );
1853 30 : sc.push_back( "HE1" );
1854 30 : sc.push_back( "HE2" );
1855 30 : sc.push_back( "HE3" );
1856 30 : return sc;
1857 324 : } else if(s=="PHE") {
1858 96 : sc.push_back( "CB" );
1859 96 : sc.push_back( "CG" );
1860 96 : sc.push_back( "CD1" );
1861 96 : sc.push_back( "CD2" );
1862 96 : sc.push_back( "CE1" );
1863 96 : sc.push_back( "CE2" );
1864 96 : sc.push_back( "CZ" );
1865 96 : sc.push_back( "HB1" );
1866 96 : sc.push_back( "HB2" );
1867 96 : sc.push_back( "HB3" );
1868 96 : sc.push_back( "HD1" );
1869 96 : sc.push_back( "HD2" );
1870 96 : sc.push_back( "HD3" );
1871 96 : sc.push_back( "HE1" );
1872 96 : sc.push_back( "HE2" );
1873 96 : sc.push_back( "HE3" );
1874 96 : sc.push_back( "HZ" );
1875 96 : return sc;
1876 228 : } else if(s=="PRO") {
1877 42 : sc.push_back( "CB" );
1878 42 : sc.push_back( "CG" );
1879 42 : sc.push_back( "CD" );
1880 42 : sc.push_back( "HB1" );
1881 42 : sc.push_back( "HB2" );
1882 42 : sc.push_back( "HB3" );
1883 42 : sc.push_back( "HG1" );
1884 42 : sc.push_back( "HG2" );
1885 42 : sc.push_back( "HG3" );
1886 42 : sc.push_back( "HD1" );
1887 42 : sc.push_back( "HD2" );
1888 42 : sc.push_back( "HD3" );
1889 42 : return sc;
1890 186 : } else if(s=="SER") {
1891 54 : sc.push_back( "CB" );
1892 54 : sc.push_back( "OG" );
1893 54 : sc.push_back( "HB1" );
1894 54 : sc.push_back( "HB2" );
1895 54 : sc.push_back( "HB3" );
1896 54 : sc.push_back( "HG1" );
1897 54 : sc.push_back( "HG" );
1898 54 : return sc;
1899 132 : } else if(s=="THR") {
1900 66 : sc.push_back( "CB" );
1901 66 : sc.push_back( "OG1" );
1902 66 : sc.push_back( "CG2" );
1903 66 : sc.push_back( "HB" );
1904 66 : sc.push_back( "HG1" );
1905 66 : sc.push_back( "HG21" );
1906 66 : sc.push_back( "HG22" );
1907 66 : sc.push_back( "HG23" );
1908 66 : sc.push_back( "1HG2" );
1909 66 : sc.push_back( "2HG2" );
1910 66 : sc.push_back( "3HG2" );
1911 66 : return sc;
1912 66 : } else if(s=="TRP") {
1913 6 : sc.push_back( "CB" );
1914 6 : sc.push_back( "CG" );
1915 6 : sc.push_back( "CD1" );
1916 6 : sc.push_back( "CD2" );
1917 6 : sc.push_back( "NE1" );
1918 6 : sc.push_back( "CE2" );
1919 6 : sc.push_back( "CE3" );
1920 6 : sc.push_back( "CZ2" );
1921 6 : sc.push_back( "CZ3" );
1922 6 : sc.push_back( "CH2" );
1923 6 : sc.push_back( "HB1" );
1924 6 : sc.push_back( "HB2" );
1925 6 : sc.push_back( "HB3" );
1926 6 : sc.push_back( "HD1" );
1927 6 : sc.push_back( "HE1" );
1928 6 : sc.push_back( "HE3" );
1929 6 : sc.push_back( "HZ2" );
1930 6 : sc.push_back( "HZ3" );
1931 6 : sc.push_back( "HH2" );
1932 6 : return sc;
1933 60 : } else if(s=="TYR") {
1934 12 : sc.push_back( "CB" );
1935 12 : sc.push_back( "CG" );
1936 12 : sc.push_back( "CD1" );
1937 12 : sc.push_back( "CD2" );
1938 12 : sc.push_back( "CE1" );
1939 12 : sc.push_back( "CE2" );
1940 12 : sc.push_back( "CZ" );
1941 12 : sc.push_back( "OH" );
1942 12 : sc.push_back( "HB1" );
1943 12 : sc.push_back( "HB2" );
1944 12 : sc.push_back( "HB3" );
1945 12 : sc.push_back( "HD1" );
1946 12 : sc.push_back( "HD2" );
1947 12 : sc.push_back( "HD3" );
1948 12 : sc.push_back( "HE1" );
1949 12 : sc.push_back( "HE2" );
1950 12 : sc.push_back( "HE3" );
1951 12 : sc.push_back( "HH" );
1952 12 : return sc;
1953 48 : } else if(s=="VAL") {
1954 48 : sc.push_back( "CB" );
1955 48 : sc.push_back( "CG1" );
1956 48 : sc.push_back( "CG2" );
1957 48 : sc.push_back( "HB" );
1958 48 : sc.push_back( "HG11" );
1959 48 : sc.push_back( "HG12" );
1960 48 : sc.push_back( "HG13" );
1961 48 : sc.push_back( "HG21" );
1962 48 : sc.push_back( "HG22" );
1963 48 : sc.push_back( "HG23" );
1964 48 : sc.push_back( "1HG1" );
1965 48 : sc.push_back( "2HG1" );
1966 48 : sc.push_back( "3HG1" );
1967 48 : sc.push_back( "1HG2" );
1968 48 : sc.push_back( "2HG2" );
1969 48 : sc.push_back( "3HG2" );
1970 48 : return sc;
1971 0 : } else plumed_merror("CS2Backbone: side_chain_atoms unknown");
1972 : }
1973 :
1974 15672 : bool CS2Backbone::isSP2(const string & resType, const string & atomName) {
1975 15672 : bool sp2 = false;
1976 15672 : if (atomName == "C") return true;
1977 14616 : if (atomName == "O") return true;
1978 :
1979 13572 : if(resType == "TRP") {
1980 :
1981 132 : if (atomName == "CG") sp2 = true;
1982 126 : else if (atomName == "CD1") sp2 = true;
1983 120 : else if (atomName == "CD2") sp2 = true;
1984 114 : else if (atomName == "CE2") sp2 = true;
1985 108 : else if (atomName == "CE3") sp2 = true;
1986 102 : else if (atomName == "CZ2") sp2 = true;
1987 96 : else if (atomName == "CZ3") sp2 = true;
1988 90 : else if (atomName == "CH2") sp2 = true;
1989 :
1990 13440 : } else if (resType == "ASP") {
1991 :
1992 552 : if (atomName == "CG") sp2 = true;
1993 498 : else if (atomName == "OD1") sp2 = true;
1994 444 : else if (atomName == "OD2") sp2 = true;
1995 :
1996 12888 : } else if (resType == "GLU") {
1997 :
1998 948 : if (atomName == "CD") sp2 = true;
1999 876 : else if (atomName == "OE1") sp2 = true;
2000 804 : else if (atomName == "OE2") sp2 = true;
2001 :
2002 11940 : } else if (resType == "ARG") {
2003 :
2004 924 : if (atomName == "CZ") sp2 = true;
2005 :
2006 11016 : } else if (resType == "HIS") {
2007 :
2008 0 : if (atomName == "CG") sp2 = true;
2009 0 : else if (atomName == "ND1") sp2 = true;
2010 0 : else if (atomName == "CD2") sp2 = true;
2011 0 : else if (atomName == "CE1") sp2 = true;
2012 0 : else if (atomName == "NE2") sp2 = true;
2013 :
2014 11016 : } else if (resType == "PHE") {
2015 :
2016 1728 : if (atomName == "CG") sp2 = true;
2017 1632 : else if (atomName == "CD1") sp2 = true;
2018 1536 : else if (atomName == "CD2") sp2 = true;
2019 1440 : else if (atomName == "CE1") sp2 = true;
2020 1344 : else if (atomName == "CE2") sp2 = true;
2021 1248 : else if (atomName == "CZ") sp2 = true;
2022 :
2023 9288 : } else if (resType == "TYR") {
2024 :
2025 228 : if (atomName == "CG") sp2 = true;
2026 216 : else if (atomName == "CD1") sp2 = true;
2027 204 : else if (atomName == "CD2") sp2 = true;
2028 192 : else if (atomName == "CE1") sp2 = true;
2029 180 : else if (atomName == "CE2") sp2 = true;
2030 168 : else if (atomName == "CZ") sp2 = true;
2031 :
2032 9060 : } else if (resType == "ASN") {
2033 :
2034 648 : if (atomName == "CG") sp2 = true;
2035 594 : else if (atomName == "OD1") sp2 = true;
2036 :
2037 8412 : } else if (resType == "GLN") {
2038 :
2039 270 : if (atomName == "CD") sp2 = true;
2040 252 : else if (atomName == "OE1") sp2 = true;
2041 :
2042 : }
2043 :
2044 13572 : return sp2;
2045 : }
2046 :
2047 15672 : bool CS2Backbone::is_chi1_cx(const string & frg, const string & atm) {
2048 15672 : if(atm=="CG") return true;
2049 15108 : if((frg == "CYS")&&(atm =="SG")) return true;
2050 15084 : if(((frg == "ILE")||(frg == "VAL"))&&(atm == "CG1")) return true;
2051 14976 : if((frg == "SER")&&(atm == "OG")) return true;
2052 14922 : if((frg == "THR")&&(atm == "OG1")) return true;
2053 :
2054 14856 : return false;
2055 : }
2056 :
2057 6192 : unsigned CS2Backbone::frag_segment(const unsigned p) {
2058 6192 : unsigned s = 0;
2059 6516 : for(unsigned i=0; i<seg_last.size()-1; i++) {
2060 6192 : if(p>seg_last[i]) s = i+1;
2061 5868 : else break;
2062 : }
2063 6192 : return s;
2064 : }
2065 :
2066 6192 : unsigned CS2Backbone::frag_relitive_index(const unsigned p, const unsigned s) {
2067 6192 : if(s==0) return p;
2068 324 : return p-seg_last[s-1];
2069 : }
2070 :
2071 0 : void CS2Backbone::debug_report() {
2072 0 : printf("\t CS2Backbone Initialization report: \n");
2073 0 : printf("\t -------------------------------\n");
2074 0 : printf("\t Number of segments: %u\n", static_cast<unsigned>(atom.size()));
2075 0 : printf("\t Segments size: ");
2076 0 : for(unsigned i=0; i<atom.size(); i++) printf("%u ", static_cast<unsigned>(atom[i].size())); printf("\n");
2077 : printf("\t%8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s \n",
2078 0 : "Seg","N","AA","Prev","Curr","Next","SC","XD1","XD2","Phi","Psi","Chi1");
2079 0 : for(unsigned i=0; i<atom.size(); i++) {
2080 0 : for(unsigned j=0; j<atom[i].size(); j++) {
2081 : printf("\t%8u %8u %8s %8u %8u %8u %8u %8u %8u %8u %8u %8u \n",
2082 : i+1, j+1,
2083 0 : atom[i][j].res_name.c_str(),
2084 0 : (unsigned)atom[i][j].prev.size(),
2085 0 : (unsigned)atom[i][j].curr.size(),
2086 0 : (unsigned)atom[i][j].next.size(),
2087 0 : (unsigned)atom[i][j].side_chain.size(),
2088 0 : (unsigned)atom[i][j].xd1.size(),
2089 0 : (unsigned)atom[i][j].xd2.size(),
2090 0 : (unsigned)atom[i][j].phi.size(),
2091 0 : (unsigned)atom[i][j].psi.size(),
2092 0 : (unsigned)atom[i][j].chi1.size());
2093 :
2094 0 : for(unsigned k=0; k<atom[i][j].prev.size(); k++) printf("%8i ", atom[i][j].prev[k]); printf("\n");
2095 0 : for(unsigned k=0; k<atom[i][j].curr.size(); k++) printf("%8i ", atom[i][j].curr[k]); printf("\n");
2096 0 : for(unsigned k=0; k<atom[i][j].next.size(); k++) printf("%8i ", atom[i][j].next[k]); printf("\n");
2097 0 : for(unsigned k=0; k<atom[i][j].side_chain.size(); k++) printf("%8i ", atom[i][j].side_chain[k]); printf("\n");
2098 0 : for(unsigned k=0; k<atom[i][j].xd1.size(); k++) printf("%8i ", atom[i][j].xd1[k]); printf("\n");
2099 0 : for(unsigned k=0; k<atom[i][j].xd2.size(); k++) printf("%8i ", atom[i][j].xd2[k]); printf("\n");
2100 0 : for(unsigned k=0; k<atom[i][j].phi.size(); k++) printf("%8i ", atom[i][j].phi[k]); printf("\n");
2101 0 : for(unsigned k=0; k<atom[i][j].psi.size(); k++) printf("%8i ", atom[i][j].psi[k]); printf("\n");
2102 0 : for(unsigned k=0; k<atom[i][j].chi1.size(); k++) printf("%8i ", atom[i][j].chi1[k]); printf("\n");
2103 :
2104 : }
2105 : }
2106 :
2107 0 : printf("\t Rings: \n");
2108 0 : printf("\t ------ \n");
2109 0 : printf("\t Number of rings: %u\n", static_cast<unsigned>(ringInfo.size()));
2110 0 : printf("\t%8s %8s %8s %8s\n", "Num","Type","RType","N.atoms");
2111 0 : for(unsigned i=0; i<ringInfo.size(); i++) {
2112 0 : printf("\t%8u %8u %8u \n",i+1,ringInfo[i].rtype,ringInfo[i].numAtoms);
2113 0 : for(unsigned j=0; j<ringInfo[i].numAtoms; j++) printf("%8u ", ringInfo[i].atom[j]); printf("\n");
2114 : }
2115 0 : }
2116 :
2117 350010 : void CS2Backbone::xdist_name_map(string & name) {
2118 350010 : if((name == "OT1")||(name == "OC1")) name = "O";
2119 350010 : else if ((name == "HN") || (name == "HT1") || (name == "H1")) name = "H";
2120 1394034 : else if ((name == "CG1")|| (name == "OG")||
2121 1047654 : (name == "SG") || (name == "OG1")) name = "CG";
2122 344826 : else if ((name == "HA1")) name = "HA";
2123 350010 : }
2124 :
2125 : }
2126 2523 : }
|
