Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 : #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 "MetainferenceBase.h"
39 : #include "core/ActionRegister.h"
40 : #include "tools/Pbc.h"
41 : #include "tools/PDB.h"
42 : #include "tools/Torsion.h"
43 :
44 : using namespace std;
45 :
46 : namespace PLMD {
47 : namespace isdb {
48 :
49 : //+PLUMEDOC ISDB_COLVAR CS2BACKBONE
50 : /*
51 : Calculates the backbone chemical shifts for a protein.
52 :
53 : The functional form is that of CamShift \cite Kohlhoff:2009us. The chemical shifts
54 : of the selected nuclei/residues are saved as components. Reference experimental values
55 : can also be stored as components. The two sets of components can then be used to calculate
56 : either a scoring function as in \cite Robustelli:2010dn \cite Granata:2013dk, using
57 : the keyword CAMSHIFT or to calculate ensemble averaged chemical shift as in \cite Camilloni:2012je
58 : \cite Camilloni:2013hs (see \ref ENSEMBLE, \ref STATS and \ref RESTRAINT). Finally they can
59 : also be used as input for \ref METAINFERENCE, \cite Bonomi:2016ip . In the current implementation there is
60 : no need to pass the data to \ref METAINFERENCE because \ref CS2BACKBONE can internally enable Metainference
61 : using the keywork DOSCORE.
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 DATADIR.
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 DATADIR: camshift.db. This file includes all the CamShift parameters
108 : and can be found in regtest/isdb/rt-cs2backbone/data/ .
109 :
110 : All the above files must be in a single folder that must be specified with the keyword DATADIR.
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 : \plumedfile
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 DATADIR=../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 : \endplumedfile
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 : \plumedfile
130 : cs: CS2BACKBONE ATOMS=1-174 DATADIR=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 : \endplumedfile
138 :
139 : This third example show how to use chemical shifts to calculate a \ref METAINFERENCE score .
140 :
141 : \plumedfile
142 : cs: CS2BACKBONE ATOMS=1-174 DATADIR=data/ NRES=13 DOSCORE NDATA=24
143 : csbias: BIASVALUE ARG=cs.score
144 :
145 : PRINT ARG=(cs\.hn_.*),(cs\.nh_.*) FILE=CS.dat STRIDE=1000
146 : PRINT ARG=cs.score FILE=BIAS STRIDE=100
147 : \endplumedfile
148 :
149 : */
150 : //+ENDPLUMEDOC
151 :
152 : class CS2BackboneDB {
153 : enum { STD, GLY, PRO};
154 : enum { HA_ATOM, H_ATOM, N_ATOM, CA_ATOM, CB_ATOM, C_ATOM };
155 : static const unsigned aa_kind = 3;
156 : static const unsigned atm_kind = 6;
157 : static const unsigned numXtraDists = 27;
158 :
159 : // ALA, ARG, ASN, ASP, CYS, GLU, GLN, GLY, HIS, ILE, LEU, LYS, MET, PHE, PRO, SER, THR, TRP, TYR, VAL
160 : double c_aa[aa_kind][atm_kind][20];
161 : double c_aa_prev[aa_kind][atm_kind][20];
162 : double c_aa_succ[aa_kind][atm_kind][20];
163 : double co_bb[aa_kind][atm_kind][16];
164 : double co_sc_[aa_kind][atm_kind][20][20];
165 : double co_xd[aa_kind][atm_kind][numXtraDists];
166 : double co_sphere[aa_kind][atm_kind][2][8];
167 : // for ring current effects
168 : // Phe, Tyr, Trp_1, Trp_2, His
169 : double co_ring[aa_kind][atm_kind][5];
170 : // for dihedral angles
171 : // co * (a * cos(3 * omega + c) + b * cos(omega + d))
172 : double co_da[aa_kind][atm_kind][3];
173 : double pars_da[aa_kind][atm_kind][3][5];
174 :
175 : public:
176 :
177 2716 : inline unsigned kind(const string &s) {
178 2716 : if(s=="GLY") return GLY;
179 2212 : else if(s=="PRO") return PRO;
180 2030 : return STD;
181 : }
182 :
183 252 : inline unsigned atom_kind(const string &s) {
184 252 : if(s=="HA")return HA_ATOM;
185 210 : else if(s=="H") return H_ATOM;
186 168 : else if(s=="N") return N_ATOM;
187 126 : else if(s=="CA")return CA_ATOM;
188 84 : else if(s=="CB")return CB_ATOM;
189 42 : else if(s=="C") return C_ATOM;
190 0 : return -1;
191 : }
192 :
193 : unsigned get_numXtraDists() {return numXtraDists;}
194 :
195 : //PARAMETERS
196 8246 : inline double * CONSTAACURR(const unsigned a_kind, const unsigned at_kind) {return c_aa[a_kind][at_kind];}
197 8246 : inline double * CONSTAANEXT(const unsigned a_kind, const unsigned at_kind) {return c_aa_succ[a_kind][at_kind];}
198 8246 : inline double * CONSTAAPREV(const unsigned a_kind, const unsigned at_kind) {return c_aa_prev[a_kind][at_kind];}
199 8246 : inline double * CONST_BB2_PREV(const unsigned a_kind, const unsigned at_kind) {return co_bb[a_kind][at_kind];}
200 : inline double * CONST_BB2_CURR(const unsigned a_kind, const unsigned at_kind) {return co_bb[a_kind][at_kind]+5;}
201 : inline double * CONST_BB2_NEXT(const unsigned a_kind, const unsigned at_kind) {return co_bb[a_kind][at_kind]+11;}
202 8246 : inline double * CONST_SC2(const unsigned a_kind, const unsigned at_kind, unsigned res_type) { return co_sc_[a_kind][at_kind][res_type];}
203 8246 : inline double * CONST_XD(const unsigned a_kind, const unsigned at_kind) { return co_xd[a_kind][at_kind];}
204 8246 : inline double * CO_SPHERE(const unsigned a_kind, const unsigned at_kind, unsigned exp_type) { return co_sphere[a_kind][at_kind][exp_type];}
205 8246 : inline double * CO_RING(const unsigned a_kind, const unsigned at_kind) { return co_ring[a_kind][at_kind];}
206 8246 : inline double * CO_DA(const unsigned a_kind, const unsigned at_kind) { return co_da[a_kind][at_kind];}
207 23072 : 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];}
208 :
209 14 : void parse(const string &file, const double dscale) {
210 28 : ifstream in;
211 14 : in.open(file.c_str());
212 14 : if(!in) plumed_merror("Unable to open CS2Backbone DB file " +file);
213 :
214 : unsigned c_kind = 0;
215 : unsigned c_atom = 0;
216 : unsigned nline = 0;
217 :
218 308 : for(unsigned i=0; i<3; i++) for(unsigned j=0; j<6; j++) {
219 10332 : for(unsigned k=0; k<20; k++) {
220 5040 : c_aa[i][j][k]=0.;
221 5040 : c_aa_prev[i][j][k]=0.;
222 5040 : c_aa_succ[i][j][k]=0.;
223 105840 : for(unsigned m=0; m<20; m++) co_sc_[i][j][k][m]=0.;
224 : }
225 4284 : for(unsigned k=0; k<16; k++) {co_bb[i][j][k]=0.; }
226 2268 : for(unsigned k=0; k<8; k++) { co_sphere[i][j][0][k]=0.; co_sphere[i][j][1][k]=0.; }
227 1764 : for(unsigned k=0; k<3; k++) {
228 756 : co_da[i][j][k]=0.;
229 4536 : for(unsigned l=0; l<5; l++) pars_da[i][j][k][l]=0.;
230 : }
231 1512 : for(unsigned k=0; k<5; k++) co_ring[i][j][k]=0.;
232 7056 : for(unsigned k=0; k<numXtraDists; k++) co_xd[i][j][k]=0.;
233 : }
234 :
235 29330 : while(!in.eof()) {
236 : string line;
237 29316 : getline(in,line);
238 : ++nline;
239 29316 : if(line.compare(0,1,"#")==0) continue;
240 0 : vector<string> tok;
241 0 : vector<string> tmp;
242 31780 : tok = split(line,' ');
243 595042 : for(unsigned q=0; q<tok.size(); q++)
244 187754 : if(tok[q].size()) tmp.push_back(tok[q]);
245 15890 : tok = tmp;
246 31780 : if(tok.size()==0) continue;
247 16128 : if(tok[0]=="PAR") {
248 252 : c_kind = kind(tok[2]);
249 252 : c_atom = atom_kind(tok[1]);
250 252 : continue;
251 : }
252 15624 : else if(tok[0]=="WEIGHT") {
253 252 : continue;
254 : }
255 15372 : else if(tok[0]=="FLATBTM") {
256 252 : continue;
257 : }
258 15120 : else if (tok[0] == "SCALEHARM") {
259 252 : continue;
260 : }
261 14868 : else if (tok[0] == "TANHAMPLI") {
262 252 : continue;
263 : }
264 14616 : else if (tok[0] == "ENDHARMON") {
265 252 : continue;
266 : }
267 14364 : else if (tok[0] == "MAXRCDEVI") {
268 252 : continue;
269 : }
270 14112 : else if (tok[0] == "RANDCOIL") {
271 252 : continue;
272 : }
273 13860 : else if (tok[0] == "CONST") {
274 252 : continue;
275 : }
276 13860 : else if (tok[0] == "CONSTAA") {
277 252 : assign(c_aa[c_kind][c_atom],tok,1);
278 252 : continue;
279 : }
280 13608 : else if (tok[0] == "CONSTAA-1") {
281 252 : assign(c_aa_prev[c_kind][c_atom],tok,1);
282 252 : continue;
283 : }
284 13356 : else if (tok[0] == "CONSTAA+1") {
285 252 : assign(c_aa_succ[c_kind][c_atom],tok,1);
286 252 : continue;
287 : }
288 12852 : else if (tok[0] == "COBB1") {
289 252 : continue;
290 : }
291 12852 : else if (tok[0] == "COBB2") {
292 : //angstrom to nm
293 252 : assign(co_bb[c_kind][c_atom],tok,dscale);
294 252 : continue;
295 : }
296 12600 : else if (tok[0] == "SPHERE1") {
297 : // angstrom^-3 to nm^-3
298 252 : assign(co_sphere[c_kind][c_atom][0],tok,1./(dscale*dscale*dscale));
299 252 : continue;
300 : }
301 12348 : else if (tok[0] == "SPHERE2") {
302 : //angstrom to nm
303 252 : assign(co_sphere[c_kind][c_atom][1],tok,dscale);
304 252 : continue;
305 : }
306 12096 : else if (tok[0] == "DIHEDRALS") {
307 252 : assign(co_da[c_kind][c_atom],tok,1);
308 252 : continue;
309 : }
310 11592 : else if (tok[0] == "RINGS") {
311 : // angstrom^-3 to nm^-3
312 252 : assign(co_ring[c_kind][c_atom],tok,1./(dscale*dscale*dscale));
313 4284 : for(unsigned i=1; i<tok.size(); i++)
314 1260 : co_ring[c_kind][c_atom][i-1] *= 1000;
315 252 : continue;
316 : }
317 11340 : else if (tok[0] == "HBONDS") {
318 252 : continue;
319 : }
320 11340 : else if (tok[0] == "XTRADISTS") {
321 : //angstrom to nm
322 252 : assign(co_xd[c_kind][c_atom],tok,dscale);
323 252 : continue;
324 : }
325 11088 : else if(tok[0]=="DIHEDPHI") {
326 252 : assign(pars_da[c_kind][c_atom][0],tok,1);
327 252 : continue;
328 : }
329 10836 : else if(tok[0]=="DIHEDPSI") {
330 252 : assign(pars_da[c_kind][c_atom][1],tok,1);
331 252 : continue;
332 : }
333 10584 : else if(tok[0]=="DIHEDCHI1") {
334 252 : assign(pars_da[c_kind][c_atom][2],tok,1);
335 252 : continue;
336 : }
337 :
338 : bool ok = false;
339 : string scIdent1 [] = {"COSCALA1", "COSCARG1", "COSCASN1", "COSCASP1", "COSCCYS1", "COSCGLN1", "COSCGLU1",
340 : "COSCGLY1", "COSCHIS1", "COSCILE1", "COSCLEU1", "COSCLYS1", "COSCMET1", "COSCPHE1",
341 : "COSCPRO1", "COSCSER1", "COSCTHR1", "COSCTRP1", "COSCTYR1", "COSCVAL1"
342 10080 : };
343 :
344 307440 : for(unsigned scC = 0; scC < 20; scC++) {
345 307440 : if(tok[0]==scIdent1[scC]) {
346 : ok = true;
347 : break;
348 : }
349 : }
350 221760 : if(ok) continue;
351 :
352 : string scIdent2 [] = {"COSCALA2", "COSCARG2", "COSCASN2", "COSCASP2", "COSCCYS2", "COSCGLN2", "COSCGLU2",
353 : "COSCGLY2", "COSCHIS2", "COSCILE2", "COSCLEU2", "COSCLYS2", "COSCMET2", "COSCPHE2",
354 : "COSCPRO2", "COSCSER2", "COSCTHR2", "COSCTRP2", "COSCTYR2", "COSCVAL2"
355 5040 : };
356 :
357 100800 : for(unsigned scC = 0; scC < 20; scC++) {
358 105840 : if(tok[0]==scIdent2[scC]) {
359 : //angstrom to nm
360 5040 : assign(co_sc_[c_kind][c_atom][scC],tok,dscale);
361 : ok = true; break;
362 : }
363 : }
364 110880 : if(ok) continue;
365 :
366 0 : if(tok.size()) {
367 0 : string str_err = "CS2Backbone DB WARNING: unrecognized token: " + tok[0];
368 0 : plumed_merror(str_err);
369 : }
370 : }
371 14 : in.close();
372 14 : }
373 :
374 : private:
375 :
376 15890 : vector<string> &split(const string &s, char delim, vector<string> &elems) {
377 31780 : stringstream ss(s);
378 : string item;
379 407288 : while (getline(ss, item, delim)) {
380 187754 : elems.push_back(item);
381 : }
382 15890 : return elems;
383 : }
384 :
385 15890 : vector<string> split(const string &s, char delim) {
386 : vector<string> elems;
387 15890 : split(s, delim, elems);
388 15890 : return elems;
389 : }
390 :
391 8064 : void assign(double * f, const vector<string> & v, const double scale) {
392 278460 : for(unsigned i=1; i<v.size(); i++) {
393 87444 : f[i-1] = scale*(atof(v[i].c_str()));
394 87444 : if(fabs(f[i-1])<0.000001) f[i-1]=0.;
395 : }
396 8064 : }
397 : };
398 :
399 42 : class CS2Backbone : public MetainferenceBase {
400 47320 : struct Fragment {
401 : vector<Value*> comp;
402 : vector<double> exp_cs;
403 : unsigned res_type_prev;
404 : unsigned res_type_curr;
405 : unsigned res_type_next;
406 : unsigned res_kind;
407 : unsigned fd;
408 : string res_name;
409 : vector<int> pos;
410 : vector<int> prev;
411 : vector<int> curr;
412 : vector<int> next;
413 : vector<int> side_chain;
414 : vector<int> xd1;
415 : vector<int> xd2;
416 : vector<unsigned> box_nb;
417 : vector<int> phi;
418 : vector<int> psi;
419 : vector<int> chi1;
420 : double t_phi;
421 : double t_psi;
422 : double t_chi1;
423 : vector<Vector> dd0, dd10, dd21, dd2;
424 :
425 44352 : Fragment() {
426 2464 : comp.resize(6);
427 2464 : exp_cs.resize(6,0);
428 2464 : res_type_prev = res_type_curr = res_type_next = 0;
429 2464 : res_kind = 0;
430 2464 : fd = 0;
431 : res_name = "";
432 2464 : pos.resize(6,-1);
433 2464 : prev.reserve(5);
434 2464 : curr.reserve(6);
435 2464 : next.reserve(5);
436 2464 : side_chain.reserve(20);
437 2464 : xd1.reserve(27);
438 2464 : xd2.reserve(27);
439 2464 : box_nb.reserve(250);
440 2464 : phi.reserve(4);
441 2464 : psi.reserve(4);
442 2464 : chi1.reserve(4);
443 2464 : t_phi = t_psi = t_chi1 = 0;
444 2464 : dd0.resize(3);
445 2464 : dd10.resize(3);
446 2464 : dd21.resize(3);
447 2464 : dd2.resize(3);
448 2464 : }
449 :
450 : };
451 :
452 : struct RingInfo {
453 : enum {R_PHE, R_TYR, R_TRP1, R_TRP2, R_HIS};
454 : unsigned rtype; // one out of five different types
455 : unsigned atom[6]; // up to six member per ring
456 : unsigned numAtoms; // number of ring members (5 or 6)
457 : Vector position; // center of ring coordinates
458 : Vector normVect; // ring plane normal vector
459 : Vector n1, n2; // two atom plane normal vectors used to compute ring plane normal
460 : Vector g[6]; // vector of the vectors used to calculate n1,n2
461 : double lengthN2; // square of length of normVect
462 : double lengthNV; // length of normVect
463 280 : RingInfo():
464 : rtype(0),numAtoms(0),
465 280 : lengthN2(NAN),lengthNV(NAN)
466 1960 : {for(unsigned i=0; i<6; i++) atom[i]=0;}
467 : };
468 :
469 : enum aa_t {ALA, ARG, ASN, ASP, CYS, GLN, GLU, GLY, HIS, ILE, LEU, LYS, MET, PHE, PRO, SER, THR, TRP, TYR, VAL, UNK};
470 : enum atom_t {D_C, D_H, D_N, D_O, D_S, D_C2, D_N2, D_O2};
471 :
472 : CS2BackboneDB db;
473 : vector<vector<Fragment> > atom;
474 : vector<vector<vector<unsigned> > > index_cs;
475 : vector<RingInfo> ringInfo;
476 : vector<unsigned> seg_last;
477 : vector<unsigned> type;
478 : vector<unsigned> res_num;
479 : unsigned box_nupdate;
480 : unsigned box_count;
481 : bool camshift;
482 : bool pbc;
483 :
484 : void remove_problematic(const string &res, const string &nucl);
485 : void read_cs(const string &file, const string &k);
486 : void update_neighb();
487 : void compute_ring_parameters();
488 : void compute_dihedrals();
489 : void init_backbone(const PDB &pdb);
490 : void init_sidechain(const PDB &pdb);
491 : void init_xdist(const PDB &pdb);
492 : void init_types(const PDB &pdb);
493 : void init_rings(const PDB &pdb);
494 : aa_t frag2enum(const string &aa);
495 : vector<string> side_chain_atoms(const string &s);
496 : bool isSP2(const string & resType, const string & atomName);
497 : bool is_chi1_cx(const string & frg, const string & atm);
498 : unsigned frag_segment(const unsigned p);
499 : unsigned frag_relitive_index(const unsigned p, const unsigned s);
500 : void debug_report();
501 : void xdist_name_map(string & name);
502 :
503 : public:
504 :
505 : explicit CS2Backbone(const ActionOptions&);
506 : static void registerKeywords( Keywords& keys );
507 : void calculate();
508 : void update();
509 : };
510 :
511 6466 : PLUMED_REGISTER_ACTION(CS2Backbone,"CS2BACKBONE")
512 :
513 15 : void CS2Backbone::registerKeywords( Keywords& keys ) {
514 15 : componentsAreNotOptional(keys);
515 15 : useCustomisableComponents(keys);
516 15 : MetainferenceBase::registerKeywords( keys );
517 45 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
518 60 : keys.add("atoms","ATOMS","The atoms to be included in the calculation, e.g. the whole protein.");
519 75 : keys.add("compulsory","DATADIR","data/","The folder with the experimental chemical shifts.");
520 75 : keys.add("compulsory","TEMPLATE","template.pdb","A PDB file of the protein system to initialise ALMOST.");
521 75 : keys.add("compulsory","NEIGH_FREQ","20","Period in step for neighbour list update.");
522 60 : keys.add("compulsory","NRES","Number of residues, corresponding to the number of chemical shifts.");
523 45 : keys.addFlag("CAMSHIFT",false,"Set to TRUE if you to calculate a single CamShift score.");
524 45 : keys.addFlag("NOEXP",false,"Set to TRUE if you don't want to have fixed components with the experimetnal values.");
525 60 : keys.addOutputComponent("ha","default","the calculated Ha hydrogen chemical shifts");
526 60 : keys.addOutputComponent("hn","default","the calculated H hydrogen chemical shifts");
527 60 : keys.addOutputComponent("nh","default","the calculated N nitrogen chemical shifts");
528 60 : keys.addOutputComponent("ca","default","the calculated Ca carbon chemical shifts");
529 60 : keys.addOutputComponent("cb","default","the calculated Cb carbon chemical shifts");
530 60 : keys.addOutputComponent("co","default","the calculated C' carbon chemical shifts");
531 60 : keys.addOutputComponent("expha","default","the experimental Ha hydrogen chemical shifts");
532 60 : keys.addOutputComponent("exphn","default","the experimental H hydrogen chemical shifts");
533 60 : keys.addOutputComponent("expnh","default","the experimental N nitrogen chemical shifts");
534 60 : keys.addOutputComponent("expca","default","the experimental Ca carbon chemical shifts");
535 60 : keys.addOutputComponent("expcb","default","the experimental Cb carbon chemical shifts");
536 60 : keys.addOutputComponent("expco","default","the experimental C' carbon chemical shifts");
537 15 : }
538 :
539 14 : CS2Backbone::CS2Backbone(const ActionOptions&ao):
540 : PLUMED_METAINF_INIT(ao),
541 : camshift(false),
542 42 : pbc(true)
543 : {
544 : string stringadb;
545 : string stringapdb;
546 :
547 28 : parseFlag("CAMSHIFT",camshift);
548 15 : if(camshift&&getDoScore()) error("It is not possible to use CAMSHIFT together with DOSCORE");
549 :
550 14 : bool nopbc=!pbc;
551 28 : parseFlag("NOPBC",nopbc);
552 14 : pbc=!nopbc;
553 :
554 14 : bool noexp=false;
555 28 : parseFlag("NOEXP",noexp);
556 :
557 : string stringa_data;
558 28 : parse("DATADIR",stringa_data);
559 :
560 : string stringa_template;
561 28 : parse("TEMPLATE",stringa_template);
562 :
563 14 : box_count=0;
564 14 : box_nupdate=20;
565 28 : parse("NEIGH_FREQ", box_nupdate);
566 :
567 : unsigned numResidues;
568 28 : parse("NRES", numResidues);
569 :
570 42 : stringadb = stringa_data + string("/camshift.db");
571 56 : stringapdb = stringa_data + string("/") + stringa_template;
572 :
573 : /* Lenght conversion (parameters are tuned for angstrom) */
574 : double scale=1.;
575 14 : if(!plumed.getAtoms().usingNaturalUnits()) {
576 14 : scale = 10.*atoms.getUnits().getLength();
577 : }
578 :
579 14 : log.printf(" Initialization of the predictor ...\n"); log.flush();
580 14 : db.parse(stringadb,scale);
581 28 : PDB pdb;
582 28 : if( !pdb.read(stringapdb,plumed.getAtoms().usingNaturalUnits(),1./scale) ) error("missing input file " + stringapdb);
583 14 : init_backbone(pdb);
584 14 : init_sidechain(pdb);
585 14 : init_xdist(pdb);
586 14 : init_types(pdb);
587 14 : init_rings(pdb);
588 : #ifndef NDEBUG
589 : debug_report();
590 : #endif
591 14 : log.printf(" Reading experimental data ...\n"); log.flush();
592 42 : stringadb = stringa_data + string("/CAshifts.dat");
593 28 : log.printf(" Initializing CA shifts %s\n", stringadb.c_str()); log.flush();
594 28 : read_cs(stringadb, "CA");
595 42 : stringadb = stringa_data + string("/CBshifts.dat");
596 28 : log.printf(" Initializing CB shifts %s\n", stringadb.c_str()); log.flush();
597 28 : read_cs(stringadb, "CB");
598 42 : stringadb = stringa_data + string("/Cshifts.dat");
599 28 : log.printf(" Initializing C' shifts %s\n", stringadb.c_str()); log.flush();
600 28 : read_cs(stringadb, "C");
601 42 : stringadb = stringa_data + string("/HAshifts.dat");
602 28 : log.printf(" Initializing HA shifts %s\n", stringadb.c_str()); log.flush();
603 28 : read_cs(stringadb, "HA");
604 42 : stringadb = stringa_data + string("/Hshifts.dat");
605 28 : log.printf(" Initializing H shifts %s\n", stringadb.c_str()); log.flush();
606 28 : read_cs(stringadb, "H");
607 42 : stringadb = stringa_data + string("/Nshifts.dat");
608 28 : log.printf(" Initializing N shifts %s\n", stringadb.c_str()); log.flush();
609 28 : read_cs(stringadb, "N");
610 :
611 : /* this is a workaround for those chemical shifts that can result in too large forces */
612 42 : remove_problematic("GLN", "CB");
613 42 : remove_problematic("ILE", "CB");
614 42 : remove_problematic("PRO", "N");
615 42 : remove_problematic("PRO", "H");
616 42 : remove_problematic("PRO", "CB");
617 42 : remove_problematic("GLY", "HA");
618 42 : remove_problematic("GLY", "CB");
619 : /* this is a workaround for those chemical shifts that are not parameterized */
620 98 : remove_problematic("HIE", "HA"); remove_problematic("HIP", "HA"); remove_problematic("HSP", "HA");
621 98 : remove_problematic("HIE", "H"); remove_problematic("HIP", "H"); remove_problematic("HSP", "H");
622 98 : remove_problematic("HIE", "N"); remove_problematic("HIP", "N"); remove_problematic("HSP", "N");
623 98 : remove_problematic("HIE", "CA"); remove_problematic("HIP", "CA"); remove_problematic("HSP", "CA");
624 98 : remove_problematic("HIE", "CB"); remove_problematic("HIP", "CB"); remove_problematic("HSP", "CB");
625 98 : remove_problematic("HIE", "C"); remove_problematic("HIP", "C"); remove_problematic("HSP", "C");
626 98 : remove_problematic("GLH", "HA"); remove_problematic("ASH", "HA"); remove_problematic("HSE", "HA");
627 98 : remove_problematic("GLH", "H"); remove_problematic("ASH", "H"); remove_problematic("HSE", "H");
628 98 : remove_problematic("GLH", "N"); remove_problematic("ASH", "N"); remove_problematic("HSE", "N");
629 98 : remove_problematic("GLH", "CA"); remove_problematic("ASH", "CA"); remove_problematic("HSE", "CA");
630 98 : remove_problematic("GLH", "CB"); remove_problematic("ASH", "CB"); remove_problematic("HSE", "CB");
631 98 : remove_problematic("GLH", "C"); remove_problematic("ASH", "C"); remove_problematic("HSE", "C");
632 42 : remove_problematic("UNK", "HA");
633 42 : remove_problematic("UNK", "H");
634 42 : remove_problematic("UNK", "N");
635 42 : remove_problematic("UNK", "CA");
636 42 : remove_problematic("UNK", "CB");
637 42 : remove_problematic("UNK", "C");
638 42 : remove_problematic("CYS", "HA");
639 42 : remove_problematic("CYS", "H");
640 42 : remove_problematic("CYS", "N");
641 42 : remove_problematic("CYS", "CA");
642 42 : remove_problematic("CYS", "CB");
643 42 : remove_problematic("CYS", "C");
644 42 : remove_problematic("HIS", "HA");
645 42 : remove_problematic("HIS", "H");
646 42 : remove_problematic("HIS", "N");
647 42 : remove_problematic("HIS", "CA");
648 42 : remove_problematic("HIS", "CB");
649 42 : remove_problematic("HIS", "C");
650 : /* done */
651 :
652 : vector<AtomNumber> atoms;
653 28 : parseAtomList("ATOMS",atoms);
654 :
655 14 : log<<" Bibliography "
656 42 : <<plumed.cite("Kohlhoff K, Robustelli P, Cavalli A, Salvatella A, Vendruscolo M, J. Am. Chem. Soc. 131, 13894 (2009)");
657 16 : if(camshift) log<<plumed.cite("Granata D, Camilloni C, Vendruscolo M, Laio A, Proc. Natl. Acad. Sci. USA 110, 6817 (2013)");
658 39 : else log<<plumed.cite("Camilloni C, Robustelli P, De Simone A, Cavalli A, Vendruscolo M, J. Am. Chem. Soc. 134, 3968 (2012)");
659 42 : log<<plumed.cite("Bonomi M, Camilloni C, Bioinformatics, 33, 3999 (2017)");
660 14 : log<<"\n";
661 :
662 112 : const string str_cs[] = {"ha_","hn_","nh_","ca_","cb_","co_"};
663 : unsigned index=0;
664 14 : if(camshift) {
665 1 : noexp = true;
666 1 : addValueWithDerivatives();
667 1 : setNotPeriodic();
668 : } else {
669 13 : if(getDoScore()) {
670 4 : index_cs.resize(atom.size(), vector<vector<unsigned> >());
671 20 : for(unsigned i=0; i<atom.size(); i++) {
672 24 : index_cs[i].resize(atom[i].size(), vector<unsigned>(6));
673 : }
674 : }
675 : unsigned l=0;
676 104 : for(unsigned i=0; i<atom.size(); i++) {
677 6916 : for(unsigned a=0; a<atom[i].size(); a++) {
678 2288 : unsigned res=index+a;
679 2288 : std::string num; Tools::convert(res,num);
680 29744 : for(unsigned at_kind=0; at_kind<6; at_kind++) {
681 27456 : if(atom[i][a].exp_cs[at_kind]!=0) {
682 7657 : if(getDoScore()) {
683 4712 : addComponent(str_cs[at_kind]+num);
684 4712 : componentIsNotPeriodic(str_cs[at_kind]+num);
685 7068 : atom[i][a].comp[at_kind] = getPntrToComponent(str_cs[at_kind]+num);
686 2356 : index_cs[i][a][at_kind]=l;
687 2356 : setParameter(atom[i][a].exp_cs[at_kind]);
688 2356 : l++;
689 : } else {
690 10602 : addComponentWithDerivatives(str_cs[at_kind]+num);
691 10602 : componentIsNotPeriodic(str_cs[at_kind]+num);
692 15903 : atom[i][a].comp[at_kind] = getPntrToComponent(str_cs[at_kind]+num);
693 : }
694 : }
695 : }
696 : }
697 26 : index += atom[i].size();
698 : }
699 13 : if(getDoScore()) Initialise(l);
700 : }
701 :
702 14 : if(!noexp) {
703 : index = 0;
704 104 : for(unsigned i=0; i<atom.size(); i++) {
705 6916 : for(unsigned a=0; a<atom[i].size(); a++) {
706 2288 : unsigned res=index+a;
707 2288 : std::string num; Tools::convert(res,num);
708 29744 : for(unsigned at_kind=0; at_kind<6; at_kind++) {
709 27456 : if(atom[i][a].exp_cs[at_kind]!=0) {
710 22971 : addComponent("exp"+str_cs[at_kind]+num);
711 22971 : componentIsNotPeriodic("exp"+str_cs[at_kind]+num);
712 22971 : Value* comp=getPntrToComponent("exp"+str_cs[at_kind]+num);
713 7657 : comp->set(atom[i][a].exp_cs[at_kind]);
714 : }
715 : }
716 : }
717 26 : index += atom[i].size();
718 : }
719 : }
720 :
721 : /* temporary check, the idea is that I can remove NRES completely */
722 : index=0;
723 112 : for(unsigned i=0; i<atom.size(); i++) {
724 28 : index += atom[i].size();
725 : }
726 14 : if(index!=numResidues) error("NRES and the number of residues in the PDB do not match!");
727 :
728 14 : requestAtoms(atoms);
729 14 : setDerivatives();
730 14 : checkRead();
731 14 : }
732 :
733 854 : void CS2Backbone::remove_problematic(const string &res, const string &nucl) {
734 : unsigned n;
735 854 : if(nucl=="HA") n=0;
736 714 : else if(nucl=="H") n=1;
737 574 : else if(nucl=="N") n=2;
738 434 : else if(nucl=="CA")n=3;
739 308 : else if(nucl=="CB")n=4;
740 126 : else if(nucl=="C") n=5;
741 : else return;
742 :
743 6832 : for(unsigned i=0; i<atom.size(); i++) {
744 454328 : for(unsigned a=0; a<atom[i].size(); a++) {
745 150304 : if(atom[i][a].res_name.c_str()==res) {
746 3304 : atom[i][a].exp_cs[n] = 0;
747 : }
748 : }
749 : }
750 : }
751 :
752 84 : void CS2Backbone::read_cs(const string &file, const string &nucl) {
753 168 : ifstream in;
754 84 : in.open(file.c_str());
755 84 : if(!in) error("CS2Backbone: Unable to open " + file);
756 84 : istream_iterator<string> iter(in), end;
757 :
758 : unsigned n;
759 84 : if(nucl=="HA") n=0;
760 70 : else if(nucl=="H") n=1;
761 56 : else if(nucl=="N") n=2;
762 42 : else if(nucl=="CA")n=3;
763 28 : else if(nucl=="CB")n=4;
764 14 : else if(nucl=="C") n=5;
765 : else return;
766 :
767 : int oldseg = -1;
768 : int oldp = -1;
769 : while(iter!=end) {
770 : string tok;
771 : tok = *iter; ++iter;
772 15120 : if(tok[0]=='#') { ++iter; continue;}
773 : int p = atoi(tok.c_str());
774 14448 : p = p - 1;
775 14448 : const int seg = frag_segment(p);
776 14448 : p = frag_relitive_index(p,seg);
777 14448 : if(oldp==-1) oldp=p;
778 14448 : if(oldseg==-1) oldseg=seg;
779 14448 : if(p<oldp&&seg==oldseg) {
780 0 : string errmsg = "Error while reading " + file + "! The same residue number has been used twice!";
781 0 : error(errmsg);
782 : }
783 : tok = *iter; ++iter;
784 : double cs = atof(tok.c_str());
785 28896 : if(atom[seg][p].pos[n]<=0) cs=0;
786 14042 : else atom[seg][p].exp_cs[n] = cs;
787 : oldseg = seg;
788 : oldp = p;
789 : }
790 84 : in.close();
791 : }
792 :
793 14 : void CS2Backbone::calculate()
794 : {
795 14 : if(pbc) makeWhole();
796 14 : if(getExchangeStep()) box_count=0;
797 14 : if(box_count==0) update_neighb();
798 :
799 14 : compute_ring_parameters();
800 14 : compute_dihedrals();
801 :
802 14 : double score = 0.;
803 :
804 14 : vector<double> camshift_sigma2(6);
805 14 : camshift_sigma2[0] = 0.08; // HA
806 14 : camshift_sigma2[1] = 0.30; // HN
807 14 : camshift_sigma2[2] = 9.00; // NH
808 14 : camshift_sigma2[3] = 1.30; // CA
809 14 : camshift_sigma2[4] = 1.56; // CB
810 14 : camshift_sigma2[5] = 1.70; // CO
811 :
812 14 : const unsigned chainsize = atom.size();
813 14 : const unsigned atleastned = 72+ringInfo.size()*6;
814 :
815 14 : vector<vector<unsigned> > all_list;
816 14 : vector<vector<Vector> > all_ff;
817 14 : if(getDoScore()) {
818 8 : all_list.resize(getNarg(), vector<unsigned>());
819 8 : all_ff.resize(getNarg(), vector<Vector>());
820 : }
821 :
822 : // CYCLE OVER MULTIPLE CHAINS
823 42 : #pragma omp parallel num_threads(OpenMP::getNumThreads())
824 140 : for(unsigned s=0; s<chainsize; s++) {
825 112 : const unsigned psize = atom[s].size();
826 : vector<Vector> omp_deriv;
827 60 : if(camshift) omp_deriv.resize(getNumberOfAtoms(), Vector(0,0,0));
828 56 : #pragma omp for reduction(+:score)
829 : // SKIP FIRST AND LAST RESIDUE OF EACH CHAIN
830 56 : for(unsigned a=1; a<psize-1; a++) {
831 :
832 2408 : const Fragment *myfrag = &atom[s][a];
833 2408 : const unsigned aa_kind = myfrag->res_kind;
834 2408 : const unsigned needed_atoms = atleastned+myfrag->box_nb.size();
835 :
836 : /* Extra Distances are the same for each residue */
837 2408 : const unsigned xdsize=myfrag->xd1.size();
838 2408 : vector<Vector> ext_distances(xdsize);
839 2408 : vector<double> ext_d(xdsize);
840 127624 : for(unsigned q=0; q<xdsize; q++) {
841 190260 : if(myfrag->xd1[q]==-1||myfrag->xd2[q]==-1) continue;
842 56056 : const Vector distance = delta(getPosition(myfrag->xd1[q]),getPosition(myfrag->xd2[q]));
843 56056 : ext_d[q] = distance.modulo();
844 112112 : ext_distances[q] = distance/ext_d[q];
845 : }
846 :
847 : // CYCLE OVER THE SIX BACKBONE CHEMICAL SHIFTS
848 31304 : for(unsigned at_kind=0; at_kind<6; at_kind++) {
849 28896 : if(atom[s][a].exp_cs[at_kind]!=0) {
850 : // Common constant and AATYPE
851 : const double * CONSTAACURR = db.CONSTAACURR(aa_kind,at_kind);
852 : const double * CONSTAANEXT = db.CONSTAANEXT(aa_kind,at_kind);
853 : const double * CONSTAAPREV = db.CONSTAAPREV(aa_kind,at_kind);
854 :
855 24738 : double cs = CONSTAACURR[myfrag->res_type_curr] +
856 8246 : CONSTAANEXT[myfrag->res_type_next] +
857 8246 : CONSTAAPREV[myfrag->res_type_prev];
858 : // this is the atom for which we are calculating the chemical shift
859 8246 : const unsigned ipos = myfrag->pos[at_kind];
860 :
861 : vector<unsigned> list;
862 8246 : list.reserve(needed_atoms);
863 8246 : list.push_back(ipos);
864 : vector<Vector> ff;
865 8246 : ff.reserve(needed_atoms);
866 8246 : ff.push_back(Vector(0,0,0));
867 :
868 :
869 : //PREV
870 : const double * CONST_BB2_PREV = db.CONST_BB2_PREV(aa_kind,at_kind);
871 8246 : const unsigned presize = myfrag->prev.size();
872 90706 : for(unsigned q=0; q<presize; q++) {
873 41230 : const double cb2pq = CONST_BB2_PREV[q];
874 47782 : if(cb2pq==0.) continue;
875 34678 : const unsigned jpos = myfrag->prev[q];
876 34678 : list.push_back(jpos);
877 104034 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
878 34678 : const double d = distance.modulo();
879 34678 : const double fact = cb2pq/d;
880 :
881 34678 : cs += cb2pq*d;
882 34678 : const Vector der = fact*distance;
883 34678 : ff[0] += der;
884 69356 : ff.push_back(-der);
885 : }
886 :
887 : //CURR
888 : const double * CONST_BB2_CURR = db.CONST_BB2_CURR(aa_kind,at_kind);
889 8246 : const unsigned cursize = myfrag->curr.size();
890 107198 : for(unsigned q=0; q<cursize; q++) {
891 49476 : const double cb2cq = CONST_BB2_CURR[q];
892 79436 : if(cb2cq==0.) continue;
893 34496 : const unsigned jpos = myfrag->curr[q];
894 34496 : if(ipos==jpos) continue;
895 34496 : list.push_back(jpos);
896 103488 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
897 34496 : const double d = distance.modulo();
898 34496 : const double fact = cb2cq/d;
899 :
900 34496 : cs += cb2cq*d;
901 34496 : const Vector der = fact*distance;
902 34496 : ff[0] += der;
903 68992 : ff.push_back(-der);
904 : }
905 :
906 : //NEXT
907 : const double * CONST_BB2_NEXT = db.CONST_BB2_NEXT(aa_kind,at_kind);
908 8246 : const unsigned nexsize = myfrag->next.size();
909 90706 : for(unsigned q=0; q<nexsize; q++) {
910 41230 : const double cb2nq = CONST_BB2_NEXT[q];
911 43414 : if(cb2nq==0.) continue;
912 39046 : const unsigned jpos = myfrag->next[q];
913 39046 : list.push_back(jpos);
914 117138 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
915 39046 : const double d = distance.modulo();
916 39046 : const double fact = cb2nq/d;
917 :
918 39046 : cs += cb2nq*d;
919 39046 : const Vector der = fact*distance;
920 39046 : ff[0] += der;
921 78092 : ff.push_back(-der);
922 : }
923 :
924 : //SIDE CHAIN
925 8246 : const double * CONST_SC2 = db.CONST_SC2(aa_kind,at_kind,myfrag->res_type_curr);
926 8246 : const unsigned sidsize = myfrag->side_chain.size();
927 162498 : for(unsigned q=0; q<sidsize; q++) {
928 77126 : const double cs2q = CONST_SC2[q];
929 165242 : if(cs2q==0.) continue;
930 33068 : const unsigned jpos = myfrag->side_chain[q];
931 33068 : if(ipos==jpos) continue;
932 33068 : list.push_back(jpos);
933 99204 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
934 33068 : const double d = distance.modulo();
935 33068 : const double fact = cs2q/d;
936 :
937 33068 : cs += cs2q*d;
938 33068 : const Vector der = fact*distance;
939 33068 : ff[0] += der;
940 66136 : ff.push_back(-der);
941 : }
942 :
943 : //EXTRA DIST
944 : const double * CONST_XD = db.CONST_XD(aa_kind,at_kind);
945 437038 : for(unsigned q=0; q<xdsize; q++) {
946 214396 : const double cxdq = CONST_XD[q];
947 239106 : if(cxdq==0.) continue;
948 421820 : if(myfrag->xd1[q]==-1||myfrag->xd2[q]==-1) continue;
949 389172 : list.push_back(myfrag->xd2[q]);
950 389172 : list.push_back(myfrag->xd1[q]);
951 194586 : cs += cxdq*ext_d[q];
952 194586 : const Vector der = cxdq*ext_distances[q];
953 194586 : ff.push_back( der);
954 389172 : ff.push_back(-der);
955 : }
956 :
957 : //NON BOND
958 : {
959 : const double * CONST_CO_SPHERE3 = db.CO_SPHERE(aa_kind,at_kind,0);
960 : const double * CONST_CO_SPHERE = db.CO_SPHERE(aa_kind,at_kind,1);
961 8246 : const unsigned boxsize = myfrag->box_nb.size();
962 1755970 : for(unsigned bat=0; bat<boxsize; bat++) {
963 1747724 : const unsigned jpos = myfrag->box_nb[bat];
964 2621586 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
965 873862 : const double d2 = distance.modulo2();
966 :
967 873862 : if(d2<cutOffDist2) {
968 151762 : double factor1 = sqrt(d2);
969 151762 : double dfactor1 = 1./factor1;
970 151762 : double factor3 = dfactor1*dfactor1*dfactor1;
971 151762 : double dfactor3 = -3.*factor3*dfactor1*dfactor1;
972 :
973 151762 : if(d2>cutOnDist2) {
974 142214 : const double af = cutOffDist2 - d2;
975 142214 : const double bf = cutOffDist2 - 3.*cutOnDist2 + 2.*d2;
976 142214 : const double cf = invswitch*af;
977 142214 : const double df = cf*af*bf;
978 142214 : factor1 *= df;
979 142214 : factor3 *= df;
980 :
981 142214 : const double d4 = d2*d2;
982 142214 : const double af1 = 15.*cutOnDist2*d2;
983 142214 : const double bf1 = -14.*d4;
984 142214 : const double cf1 = -3.*cutOffDist2*cutOnDist2 + cutOffDist2*d2;
985 142214 : const double df1 = af1+bf1+cf1;
986 142214 : dfactor1 *= cf*(cutOffDist4+df1);
987 :
988 : const double af3 = +2.*cutOffDist2*cutOnDist2;
989 142214 : const double bf3 = d2*(cutOffDist2+cutOnDist2);
990 142214 : const double cf3 = -2.*d4;
991 142214 : const double df3 = (af3+bf3+cf3)*d2;
992 142214 : dfactor3 *= invswitch*(cutMixed+df3);
993 : }
994 :
995 303524 : const unsigned t = type[jpos];
996 151762 : cs += factor1*CONST_CO_SPHERE[t] + factor3*CONST_CO_SPHERE3[t] ;
997 151762 : const double fact = dfactor1*CONST_CO_SPHERE[t]+dfactor3*CONST_CO_SPHERE3[t];
998 151762 : const Vector der = fact*distance;
999 :
1000 151762 : list.push_back(jpos);
1001 151762 : ff[0] += der;
1002 303524 : ff.push_back(-der);
1003 : }
1004 : }
1005 : }
1006 : //END NON BOND
1007 :
1008 : //RINGS
1009 : {
1010 : const double *rc = db.CO_RING(aa_kind,at_kind);
1011 8246 : const unsigned rsize = ringInfo.size();
1012 338086 : for(unsigned i=0; i<rsize; i++) {
1013 : // compute angle from ring middle point to current atom position
1014 : // get distance vector from query atom to ring center and normal vector to ring plane
1015 329840 : const Vector n = ringInfo[i].normVect;
1016 164920 : const double nL = ringInfo[i].lengthNV;
1017 164920 : const double inL2 = ringInfo[i].lengthN2;
1018 :
1019 329840 : const Vector d = delta(ringInfo[i].position, getPosition(ipos));
1020 164920 : const double dL2 = d.modulo2();
1021 164920 : double dL = sqrt(dL2);
1022 164920 : const double idL3 = 1./(dL2*dL);
1023 :
1024 164920 : const double dn = dotProduct(d,n);
1025 164920 : const double dn2 = dn*dn;
1026 164920 : const double dLnL = dL*nL;
1027 164920 : const double dL_nL = dL/nL;
1028 :
1029 164920 : const double ang2 = dn2*inL2/dL2;
1030 164920 : const double u = 1.-3.*ang2;
1031 164920 : const double cc = rc[ringInfo[i].rtype];
1032 :
1033 164920 : cs += cc*u*idL3;
1034 :
1035 164920 : const double fUU = -6*dn*inL2;
1036 164920 : const double fUQ = fUU/dL;
1037 164920 : const Vector gradUQ = fUQ*(dL2*n - dn*d);
1038 164920 : const Vector gradVQ = (3*dL*u)*d;
1039 :
1040 164920 : const double fact = cc*idL3*idL3;
1041 329840 : ff[0] += fact*(gradUQ - gradVQ);
1042 :
1043 164920 : const double fU = fUU/nL;
1044 : double OneOverN = 1./6.;
1045 164920 : if(ringInfo[i].numAtoms==5) OneOverN=1./3.;
1046 164920 : const Vector factor2 = OneOverN*n;
1047 164920 : const Vector factor4 = (OneOverN/dL_nL)*d;
1048 :
1049 164920 : const Vector gradV = -OneOverN*gradVQ;
1050 :
1051 164920 : if(ringInfo[i].numAtoms==6) {
1052 : // update forces on ring atoms
1053 2036762 : for(unsigned at=0; at<6; at++) {
1054 940044 : const Vector ab = crossProduct(d,ringInfo[i].g[at]);
1055 940044 : const Vector c = crossProduct(n,ringInfo[i].g[at]);
1056 940044 : const Vector factor3 = 0.5*dL_nL*c;
1057 940044 : const Vector factor1 = 0.5*ab;
1058 940044 : const Vector gradU = fU*( dLnL*(factor1 - factor2) -dn*(factor3 - factor4) );
1059 1880088 : ff.push_back(fact*(gradU - gradV));
1060 940044 : list.push_back(ringInfo[i].atom[at]);
1061 : }
1062 : } else {
1063 57722 : for(unsigned at=0; at<3; at++) {
1064 24738 : const Vector ab = crossProduct(d,ringInfo[i].g[at]);
1065 24738 : const Vector c = crossProduct(n,ringInfo[i].g[at]);
1066 24738 : const Vector factor3 = dL_nL*c;
1067 24738 : const Vector factor1 = ab;
1068 24738 : const Vector gradU = fU*( dLnL*(factor1 - factor2) -dn*(factor3 - factor4) );
1069 49476 : ff.push_back(fact*(gradU - gradV));
1070 : }
1071 8246 : list.push_back(ringInfo[i].atom[0]);
1072 8246 : list.push_back(ringInfo[i].atom[2]);
1073 8246 : list.push_back(ringInfo[i].atom[3]);
1074 : }
1075 : }
1076 : }
1077 : //END OF RINGS
1078 :
1079 : //DIHEDRAL ANGLES
1080 : {
1081 : const double *CO_DA = db.CO_DA(aa_kind,at_kind);
1082 8246 : if(myfrag->phi.size()==4) {
1083 : const double *PARS_DA = db.PARS_DA(aa_kind,at_kind,0);
1084 8246 : const double val1 = 3.*myfrag->t_phi+PARS_DA[3];
1085 8246 : const double val2 = myfrag->t_phi+PARS_DA[4];
1086 8246 : cs += CO_DA[0]*(PARS_DA[0]*cos(val1)+PARS_DA[1]*cos(val2)+PARS_DA[2]);
1087 8246 : const double fact = -CO_DA[0]*(+3.*PARS_DA[0]*sin(val1)+PARS_DA[1]*sin(val2));
1088 :
1089 16492 : ff.push_back(fact*myfrag->dd0[0]);
1090 16492 : ff.push_back(fact*myfrag->dd10[0]);
1091 16492 : ff.push_back(fact*myfrag->dd21[0]);
1092 16492 : ff.push_back(-fact*myfrag->dd2[0]);
1093 16492 : list.push_back(myfrag->phi[0]);
1094 16492 : list.push_back(myfrag->phi[1]);
1095 16492 : list.push_back(myfrag->phi[2]);
1096 16492 : list.push_back(myfrag->phi[3]);
1097 : }
1098 :
1099 8246 : if(myfrag->psi.size()==4) {
1100 : const double *PARS_DA = db.PARS_DA(aa_kind,at_kind,1);
1101 8246 : const double val1 = 3.*myfrag->t_psi+PARS_DA[3];
1102 8246 : const double val2 = myfrag->t_psi+PARS_DA[4];
1103 8246 : cs += CO_DA[1]*(PARS_DA[0]*cos(val1)+PARS_DA[1]*cos(val2)+PARS_DA[2]);
1104 8246 : const double fact = -CO_DA[1]*(+3.*PARS_DA[0]*sin(val1)+PARS_DA[1]*sin(val2));
1105 :
1106 16492 : ff.push_back(fact*myfrag->dd0[1]);
1107 16492 : ff.push_back(fact*myfrag->dd10[1]);
1108 16492 : ff.push_back(fact*myfrag->dd21[1]);
1109 16492 : ff.push_back(-fact*myfrag->dd2[1]);
1110 16492 : list.push_back(myfrag->psi[0]);
1111 16492 : list.push_back(myfrag->psi[1]);
1112 16492 : list.push_back(myfrag->psi[2]);
1113 16492 : list.push_back(myfrag->psi[3]);
1114 : }
1115 :
1116 : //Chi
1117 8246 : if(myfrag->chi1.size()==4) {
1118 : const double *PARS_DA = db.PARS_DA(aa_kind,at_kind,2);
1119 6580 : const double val1 = 3.*myfrag->t_chi1+PARS_DA[3];
1120 6580 : const double val2 = myfrag->t_chi1+PARS_DA[4];
1121 6580 : cs += CO_DA[2]*(PARS_DA[0]*cos(val1)+PARS_DA[1]*cos(val2)+PARS_DA[2]);
1122 6580 : const double fact = -CO_DA[2]*(+3.*PARS_DA[0]*sin(val1)+PARS_DA[1]*sin(val2));
1123 :
1124 13160 : ff.push_back(fact*myfrag->dd0[2]);
1125 13160 : ff.push_back(fact*myfrag->dd10[2]);
1126 13160 : ff.push_back(fact*myfrag->dd21[2]);
1127 13160 : ff.push_back(-fact*myfrag->dd2[2]);
1128 13160 : list.push_back(myfrag->chi1[0]);
1129 13160 : list.push_back(myfrag->chi1[1]);
1130 13160 : list.push_back(myfrag->chi1[2]);
1131 13160 : list.push_back(myfrag->chi1[3]);
1132 : }
1133 : }
1134 : //END OF DIHE
1135 :
1136 8246 : if(getDoScore()) {
1137 2356 : setCalcData(index_cs[s][a][at_kind], cs);
1138 4712 : all_list[index_cs[s][a][at_kind]] = list;
1139 4712 : all_ff[index_cs[s][a][at_kind]] = ff;
1140 2356 : Value *comp = atom[s][a].comp[at_kind];
1141 : comp->set(cs);
1142 : } else {
1143 5890 : if(camshift) {
1144 1178 : score += (cs - atom[s][a].exp_cs[at_kind])*(cs - atom[s][a].exp_cs[at_kind])/camshift_sigma2[at_kind];
1145 589 : double fact = 2.0*(cs - atom[s][a].exp_cs[at_kind])/camshift_sigma2[at_kind];
1146 625603 : for(unsigned i=0; i<list.size(); i++) omp_deriv[list[i]] += fact*ff[i];
1147 : } else {
1148 5301 : Value *comp = atom[s][a].comp[at_kind];
1149 : comp->set(cs);
1150 5301 : Tensor virial;
1151 3380793 : for(unsigned i=0; i<list.size(); i++) {
1152 2246794 : setAtomsDerivatives(comp,list[i],ff[i]);
1153 2246794 : virial-=Tensor(getPosition(list[i]),ff[i]);
1154 : }
1155 5301 : setBoxDerivatives(comp,virial);
1156 : }
1157 : }
1158 : }
1159 : }
1160 : }
1161 112 : #pragma omp critical
1162 31408 : if(camshift) for(unsigned i=0; i<getPositions().size(); i++) setAtomsDerivatives(getPntrToValue(),i,omp_deriv[i]);
1163 : }
1164 :
1165 14 : if(getDoScore()) {
1166 : /* Metainference */
1167 4 : double score = getScore();
1168 : setScore(score);
1169 :
1170 : /* calculate final derivatives */
1171 8 : Value* val=getPntrToComponent("score");
1172 :
1173 4 : Tensor virial;
1174 :
1175 7076 : for(unsigned i=0; i<all_list.size(); i++) {
1176 : const double fact = getMetaDer(i);
1177 1502480 : for(unsigned j=0; j<all_list[i].size(); j++) {
1178 998512 : setAtomsDerivatives(val, all_list[i][j], all_ff[i][j]*fact);
1179 1497768 : virial-=Tensor(getPosition(all_list[i][j]), all_ff[i][j]*fact);
1180 : }
1181 : }
1182 4 : setBoxDerivatives(val,virial);
1183 : }
1184 :
1185 : // in the case of camshift we calculate the virial at the end
1186 14 : if(camshift) {
1187 1 : Tensor virial;
1188 : unsigned nat=getNumberOfAtoms();
1189 1 : Value* val=getPntrToValue();
1190 10449 : for(unsigned i=0; i<nat; i++) virial-=Tensor(getPosition(i),
1191 7836 : Vector(val->getDerivative(3*i+0),
1192 : val->getDerivative(3*i+1),
1193 2612 : val->getDerivative(3*i+2)));
1194 1 : setBoxDerivatives(val,virial);
1195 1 : setValue(score);
1196 : }
1197 :
1198 14 : ++box_count;
1199 14 : if(box_count == box_nupdate) box_count = 0;
1200 14 : }
1201 :
1202 14 : void CS2Backbone::update_neighb() {
1203 14 : const unsigned chainsize = atom.size();
1204 70 : for(unsigned s=0; s<chainsize; s++) {
1205 56 : const unsigned psize = atom[s].size();
1206 84 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1207 56 : for(unsigned a=1; a<psize-1; a++) {
1208 : const unsigned boxsize = getNumberOfAtoms();
1209 2408 : atom[s][a].box_nb.clear();
1210 2408 : atom[s][a].box_nb.reserve(300);
1211 4816 : const unsigned res_curr = res_num[atom[s][a].pos[0]];
1212 6292104 : for(unsigned bat=0; bat<boxsize; bat++) {
1213 12579392 : const unsigned res_dist = abs(static_cast<int>(res_curr-res_num[bat]));
1214 6397062 : if(res_dist<2) continue;
1215 78235586 : for(unsigned at_kind=0; at_kind<6; at_kind++) {
1216 88150694 : if(atom[s][a].exp_cs[at_kind]==0.) continue;
1217 20652670 : const unsigned ipos = atom[s][a].pos[at_kind];
1218 41305340 : const Vector distance = delta(getPosition(bat),getPosition(ipos));
1219 20652670 : const double d2=distance.modulo2();
1220 20652670 : if(d2<cutOffNB2) {
1221 241160 : atom[s][a].box_nb.push_back(bat);
1222 241160 : break;
1223 : }
1224 : }
1225 : }
1226 : }
1227 : }
1228 14 : }
1229 :
1230 14 : void CS2Backbone::compute_ring_parameters() {
1231 868 : for(unsigned i=0; i<ringInfo.size(); i++) {
1232 280 : const unsigned size = ringInfo[i].numAtoms;
1233 280 : if(size==6) {
1234 798 : ringInfo[i].g[0] = delta(getPosition(ringInfo[i].atom[4]),getPosition(ringInfo[i].atom[2]));
1235 798 : ringInfo[i].g[1] = delta(getPosition(ringInfo[i].atom[5]),getPosition(ringInfo[i].atom[3]));
1236 798 : ringInfo[i].g[2] = delta(getPosition(ringInfo[i].atom[0]),getPosition(ringInfo[i].atom[4]));
1237 798 : ringInfo[i].g[3] = delta(getPosition(ringInfo[i].atom[1]),getPosition(ringInfo[i].atom[5]));
1238 798 : ringInfo[i].g[4] = delta(getPosition(ringInfo[i].atom[2]),getPosition(ringInfo[i].atom[0]));
1239 798 : ringInfo[i].g[5] = delta(getPosition(ringInfo[i].atom[3]),getPosition(ringInfo[i].atom[1]));
1240 266 : vector<Vector> a(size);
1241 532 : a[0] = getPosition(ringInfo[i].atom[0]);
1242 : // ring center
1243 266 : Vector midP = a[0];
1244 2926 : for(unsigned j=1; j<size; j++) {
1245 3990 : a[j] = getPosition(ringInfo[i].atom[j]);
1246 1330 : midP += a[j];
1247 : }
1248 266 : midP /= (double) size;
1249 266 : ringInfo[i].position = midP;
1250 : // compute normal vector to plane containing first three atoms in array
1251 798 : ringInfo[i].n1 = crossProduct(delta(a[0],a[4]), delta(a[0],a[2]));
1252 : // compute normal vector to plane containing last three atoms in array
1253 : // NB: third atom of five-membered ring used for both computations above
1254 798 : ringInfo[i].n2 = crossProduct(delta(a[3],a[1]), delta(a[3],a[5]));
1255 : // ring plane normal vector is average of n1 and n2
1256 532 : ringInfo[i].normVect = 0.5*(ringInfo[i].n1 + ringInfo[i].n2);
1257 : } else {
1258 42 : ringInfo[i].g[0] = delta(getPosition(ringInfo[i].atom[3]),getPosition(ringInfo[i].atom[2]));
1259 42 : ringInfo[i].g[1] = delta(getPosition(ringInfo[i].atom[0]),getPosition(ringInfo[i].atom[3]));
1260 42 : ringInfo[i].g[2] = delta(getPosition(ringInfo[i].atom[2]),getPosition(ringInfo[i].atom[0]));
1261 14 : vector<Vector> a(size);
1262 154 : for(unsigned j=0; j<size; j++) {
1263 210 : a[j] = getPosition(ringInfo[i].atom[j]);
1264 : }
1265 : // ring center
1266 14 : Vector midP = (a[0]+a[2]+a[3])/3.;
1267 14 : ringInfo[i].position = midP;
1268 : // compute normal vector to plane containing first three atoms in array
1269 42 : ringInfo[i].n1 = crossProduct(delta(a[0],a[3]), delta(a[0],a[2]));
1270 : // ring plane normal vector is average of n1 and n2
1271 14 : ringInfo[i].normVect = ringInfo[i].n1;
1272 :
1273 : }
1274 : // calculate squared length and length of normal vector
1275 280 : ringInfo[i].lengthN2 = 1./ringInfo[i].normVect.modulo2();
1276 560 : ringInfo[i].lengthNV = 1./sqrt(ringInfo[i].lengthN2);
1277 : }
1278 14 : }
1279 :
1280 14 : void CS2Backbone::compute_dihedrals() {
1281 14 : const unsigned chainsize = atom.size();
1282 70 : for(unsigned s=0; s<chainsize; s++) {
1283 56 : const unsigned psize = atom[s].size();
1284 84 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
1285 56 : for(unsigned a=1; a<psize-1; a++) {
1286 2408 : const Fragment *myfrag = &atom[s][a];
1287 2408 : if(myfrag->phi.size()==4) {
1288 7224 : const Vector d0 = delta(getPosition(myfrag->phi[1]), getPosition(myfrag->phi[0]));
1289 7224 : const Vector d1 = delta(getPosition(myfrag->phi[2]), getPosition(myfrag->phi[1]));
1290 7224 : const Vector d2 = delta(getPosition(myfrag->phi[3]), getPosition(myfrag->phi[2]));
1291 : Torsion t;
1292 2408 : Vector dd0, dd1, dd2;
1293 2408 : atom[s][a].t_phi = t.compute(d0,d1,d2,dd0,dd1,dd2);
1294 2408 : atom[s][a].dd0[0] = dd0;
1295 2408 : atom[s][a].dd10[0] = dd1-dd0;
1296 2408 : atom[s][a].dd21[0] = dd2-dd1;
1297 2408 : atom[s][a].dd2[0] = dd2;
1298 : }
1299 2408 : if(myfrag->psi.size()==4) {
1300 7224 : const Vector d0 = delta(getPosition(myfrag->psi[1]), getPosition(myfrag->psi[0]));
1301 7224 : const Vector d1 = delta(getPosition(myfrag->psi[2]), getPosition(myfrag->psi[1]));
1302 7224 : const Vector d2 = delta(getPosition(myfrag->psi[3]), getPosition(myfrag->psi[2]));
1303 : Torsion t;
1304 2408 : Vector dd0, dd1, dd2;
1305 2408 : atom[s][a].t_psi = t.compute(d0,d1,d2,dd0,dd1,dd2);
1306 2408 : atom[s][a].dd0[1] = dd0;
1307 2408 : atom[s][a].dd10[1] = dd1-dd0;
1308 2408 : atom[s][a].dd21[1] = dd2-dd1;
1309 2408 : atom[s][a].dd2[1] = dd2;
1310 : }
1311 2408 : if(myfrag->chi1.size()==4) {
1312 5586 : const Vector d0 = delta(getPosition(myfrag->chi1[1]), getPosition(myfrag->chi1[0]));
1313 5586 : const Vector d1 = delta(getPosition(myfrag->chi1[2]), getPosition(myfrag->chi1[1]));
1314 5586 : const Vector d2 = delta(getPosition(myfrag->chi1[3]), getPosition(myfrag->chi1[2]));
1315 : Torsion t;
1316 1862 : Vector dd0, dd1, dd2;
1317 1862 : atom[s][a].t_chi1 = t.compute(d0,d1,d2,dd0,dd1,dd2);
1318 1862 : atom[s][a].dd0[2] = dd0;
1319 1862 : atom[s][a].dd10[2] = dd1-dd0;
1320 1862 : atom[s][a].dd21[2] = dd2-dd1;
1321 1862 : atom[s][a].dd2[2] = dd2;
1322 : }
1323 : }
1324 : }
1325 14 : }
1326 :
1327 14 : void CS2Backbone::init_backbone(const PDB &pdb) {
1328 : // number of chains
1329 14 : vector<string> chains;
1330 14 : pdb.getChainNames( chains );
1331 28 : seg_last.resize(chains.size());
1332 : unsigned old_size=0;
1333 :
1334 112 : for(unsigned i=0; i<chains.size(); i++) {
1335 : unsigned start, end;
1336 : string errmsg;
1337 28 : pdb.getResidueRange( chains[i], start, end, errmsg );
1338 :
1339 28 : unsigned res_offset = start;
1340 28 : unsigned resrange = end-start+1;
1341 :
1342 42 : if(i==0) seg_last[i] = resrange;
1343 42 : else seg_last[i] = seg_last[i-1]+resrange;
1344 :
1345 : vector<int> N_;
1346 : vector<int> H_;
1347 : vector<int> CA_;
1348 : vector<int> CB_;
1349 : vector<int> HA_;
1350 : vector<int> C_;
1351 : vector<int> O_;
1352 : vector<int> CX_;
1353 28 : N_.resize (resrange,-1);
1354 28 : H_.resize (resrange,-1);
1355 28 : CA_.resize(resrange,-1);
1356 28 : CB_.resize(resrange,-1);
1357 28 : HA_.resize(resrange,-1);
1358 28 : C_.resize (resrange,-1);
1359 28 : O_.resize (resrange,-1);
1360 28 : CX_.resize(resrange,-1);
1361 :
1362 28 : vector<AtomNumber> allatoms = pdb.getAtomsInChain(chains[i]);
1363 : // cycle over all the atoms in the chain
1364 109760 : for(unsigned a=0; a<allatoms.size(); a++) {
1365 36568 : unsigned atm_index=a+old_size;
1366 36568 : unsigned f = pdb.getResidueNumber(allatoms[a]);
1367 36568 : unsigned f_idx = f-res_offset;
1368 36568 : string AN = pdb.getAtomName(allatoms[a]);
1369 36568 : string RES = pdb.getResidueName(allatoms[a]);
1370 39032 : if(AN=="N") N_[f_idx] = atm_index;
1371 68208 : else if(AN=="H" ||AN=="HN" ) H_[f_idx] = atm_index;
1372 95718 : else if(AN=="HA"||AN=="HA1"||AN=="HA3") HA_[f_idx] = atm_index;
1373 34230 : else if(AN=="CA") CA_[f_idx] = atm_index;
1374 28882 : else if(AN=="CB") CB_[f_idx] = atm_index;
1375 27258 : else if(AN=="C" ) C_ [f_idx] = atm_index;
1376 24766 : else if(AN=="O" ) O_ [f_idx] = atm_index;
1377 20734 : else if(AN=="CD"&&RES=="PRO") H_[f_idx] = atm_index;
1378 38472 : if(is_chi1_cx(RES,AN)) CX_[f_idx] = atm_index;
1379 : }
1380 28 : old_size+=allatoms.size();
1381 :
1382 : // vector of residues for a given chain
1383 28 : vector<Fragment> atm_;
1384 : // cycle over all residues in the chain
1385 2492 : for(unsigned a=start; a<=end; a++) {
1386 2464 : unsigned f_idx = a - res_offset;
1387 4928 : Fragment at;
1388 4928 : at.pos[0] = HA_[f_idx];
1389 2464 : at.pos[1] = H_[f_idx];
1390 2464 : at.pos[2] = N_[f_idx];
1391 2464 : at.pos[3] = CA_[f_idx];
1392 2464 : at.pos[4] = CB_[f_idx];
1393 2464 : at.pos[5] = C_[f_idx];
1394 2464 : at.res_type_prev = at.res_type_curr = at.res_type_next = 0;
1395 4928 : at.res_name = pdb.getResidueName(a, chains[i]);
1396 2464 : at.res_kind = db.kind(at.res_name);
1397 2464 : at.fd = a;
1398 : //REGISTER PREV CURR NEXT
1399 : {
1400 2464 : if(a>start) {
1401 4872 : at.prev.push_back( N_[f_idx-1]);
1402 2436 : at.prev.push_back(CA_[f_idx-1]);
1403 2436 : at.prev.push_back(HA_[f_idx-1]);
1404 2436 : at.prev.push_back( C_[f_idx-1]);
1405 2436 : at.prev.push_back( O_[f_idx-1]);
1406 4872 : at.res_type_prev = frag2enum(pdb.getResidueName(a-1, chains[i]));
1407 : }
1408 :
1409 2464 : at.curr.push_back( N_[f_idx]);
1410 2464 : at.curr.push_back( H_[f_idx]);
1411 2464 : at.curr.push_back(CA_[f_idx]);
1412 2464 : at.curr.push_back(HA_[f_idx]);
1413 2464 : at.curr.push_back( C_[f_idx]);
1414 2464 : at.curr.push_back( O_[f_idx]);
1415 4928 : at.res_type_curr = frag2enum(pdb.getResidueName(a, chains[i]));
1416 :
1417 2464 : if(a<end) {
1418 4872 : at.next.push_back (N_[f_idx+1]);
1419 2436 : at.next.push_back (H_[f_idx+1]);
1420 2436 : at.next.push_back(CA_[f_idx+1]);
1421 2436 : at.next.push_back(HA_[f_idx+1]);
1422 2436 : at.next.push_back( C_[f_idx+1]);
1423 4872 : at.res_type_next = frag2enum(pdb.getResidueName(a+1, chains[i]));
1424 : }
1425 :
1426 : //PHI | PSI | CH1
1427 2464 : if(a>start) {
1428 4872 : at.phi.push_back( C_[f_idx-1]);
1429 2436 : at.phi.push_back( N_[f_idx]);
1430 2436 : at.phi.push_back(CA_[f_idx]);
1431 2436 : at.phi.push_back( C_[f_idx]);
1432 : }
1433 :
1434 2464 : if(a<end) {
1435 2436 : at.psi.push_back( N_[f_idx]);
1436 2436 : at.psi.push_back(CA_[f_idx]);
1437 2436 : at.psi.push_back( C_[f_idx]);
1438 4872 : at.psi.push_back( N_[f_idx+1]);
1439 : }
1440 :
1441 4368 : if(CX_[f_idx]!=-1&&CB_[f_idx]!=-1) {
1442 1904 : at.chi1.push_back( N_[f_idx]);
1443 1904 : at.chi1.push_back(CA_[f_idx]);
1444 1904 : at.chi1.push_back(CB_[f_idx]);
1445 1904 : at.chi1.push_back(CX_[f_idx]);
1446 : }
1447 : }
1448 2464 : atm_.push_back(at);
1449 : }
1450 28 : atom.push_back(atm_);
1451 : }
1452 14 : }
1453 :
1454 14 : void CS2Backbone::init_sidechain(const PDB &pdb) {
1455 14 : vector<string> chains;
1456 14 : pdb.getChainNames( chains );
1457 : unsigned old_size=0;
1458 : // cycle over chains
1459 112 : for(unsigned s=0; s<atom.size(); s++) {
1460 : AtomNumber astart, aend;
1461 : string errmsg;
1462 28 : pdb.getAtomRange( chains[s], astart, aend, errmsg );
1463 : unsigned atom_offset = astart.index();
1464 : // cycle over residues
1465 7448 : for(unsigned a=0; a<atom[s].size(); a++) {
1466 4928 : if(atom[s][a].res_name=="UNK") continue;
1467 2464 : vector<AtomNumber> atm = pdb.getAtomsInResidue(atom[s][a].fd, chains[s]);
1468 4928 : vector<string> sc_atm = side_chain_atoms(atom[s][a].res_name);
1469 :
1470 90902 : for(unsigned sc=0; sc<sc_atm.size(); sc++) {
1471 1593970 : for(unsigned aa=0; aa<atm.size(); aa++) {
1472 1024436 : if(pdb.getAtomName(atm[aa])==sc_atm[sc]) {
1473 65394 : atom[s][a].side_chain.push_back(atm[aa].index()-atom_offset+old_size);
1474 : }
1475 : }
1476 : }
1477 :
1478 : }
1479 28 : old_size = aend.index()+1;
1480 : }
1481 14 : }
1482 :
1483 14 : void CS2Backbone::init_xdist(const PDB &pdb) {
1484 : const string atomsP1[] = {"H", "H", "H", "C", "C", "C",
1485 : "O", "O", "O", "N", "N", "N",
1486 : "O", "O", "O", "N", "N", "N",
1487 : "CG", "CG", "CG", "CG", "CG",
1488 : "CG", "CG", "CA"
1489 392 : };
1490 :
1491 14 : const int resOffsetP1 [] = {0, 0, 0, -1, -1, -1,
1492 : 0, 0, 0, 1, 1, 1,
1493 : -1, -1, -1, 0, 0, 0,
1494 : 0, 0, 0, 0, 0, -1, 1, -1
1495 : };
1496 :
1497 : const string atomsP2[] = {"HA", "C", "CB", "HA", "C", "CB",
1498 : "HA", "N", "CB", "HA", "N", "CB",
1499 : "HA", "N", "CB", "HA", "N", "CB",
1500 : "HA", "N", "C", "C", "N", "CA", "CA", "CA"
1501 392 : };
1502 :
1503 14 : 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};
1504 :
1505 14 : vector<string> chains;
1506 14 : pdb.getChainNames( chains );
1507 : unsigned old_size=0;
1508 112 : for(unsigned s=0; s<atom.size(); s++) {
1509 : AtomNumber astart, aend;
1510 : string errmsg;
1511 28 : pdb.getAtomRange( chains[s], astart, aend, errmsg );
1512 : unsigned atom_offset = astart.index();
1513 :
1514 7280 : for(unsigned a=1; a<atom[s].size()-1; a++) {
1515 2408 : vector<AtomNumber> atm_curr = pdb.getAtomsInResidue(atom[s][a].fd,chains[s]);
1516 2408 : vector<AtomNumber> atm_prev = pdb.getAtomsInResidue(atom[s][a].fd-1,chains[s]);
1517 2408 : vector<AtomNumber> atm_next = pdb.getAtomsInResidue(atom[s][a].fd+1,chains[s]);
1518 :
1519 127624 : for(unsigned q=0; q<db.get_numXtraDists()-1; q++) {
1520 : vector<AtomNumber>::iterator at1, at1_end;
1521 : vector<AtomNumber>::iterator at2, at2_end;
1522 :
1523 : bool init_p1=false;
1524 : AtomNumber p1;
1525 : bool init_p2=false;
1526 : AtomNumber p2;
1527 :
1528 62608 : if(resOffsetP1[q]== 0) { at1 = atm_curr.begin(); at1_end = atm_curr.end();}
1529 62608 : if(resOffsetP1[q]==-1) { at1 = atm_prev.begin(); at1_end = atm_prev.end();}
1530 62608 : if(resOffsetP1[q]==+1) { at1 = atm_next.begin(); at1_end = atm_next.end();}
1531 496846 : while(at1!=at1_end) {
1532 492730 : AtomNumber aa = *at1;
1533 : ++at1;
1534 492730 : string name = pdb.getAtomName(aa);
1535 :
1536 492730 : xdist_name_map(name);
1537 :
1538 492730 : if(name==atomsP1[q]) {
1539 58492 : p1 = aa;
1540 : init_p1=true;
1541 : break;
1542 : }
1543 : }
1544 :
1545 62608 : if(resOffsetP2[q]== 0) { at2 = atm_curr.begin(); at2_end = atm_curr.end();}
1546 62608 : if(resOffsetP2[q]==-1) { at2 = atm_prev.begin(); at2_end = atm_prev.end();}
1547 62608 : if(resOffsetP2[q]==+1) { at2 = atm_next.begin(); at2_end = atm_next.end();}
1548 326396 : while(at2!=at2_end) {
1549 323960 : AtomNumber aa = *at2;
1550 : ++at2;
1551 323960 : string name = pdb.getAtomName(aa);
1552 :
1553 323960 : xdist_name_map(name);
1554 :
1555 323960 : if(name==atomsP2[q]) {
1556 60172 : p2 = aa;
1557 : init_p2=true;
1558 : break;
1559 : }
1560 : }
1561 62608 : int add1 = -1;
1562 62608 : int add2 = -1;
1563 62608 : if(init_p1) add1=p1.index()-atom_offset+old_size;
1564 62608 : if(init_p2) add2=p2.index()-atom_offset+old_size;
1565 62608 : atom[s][a].xd1.push_back(add1);
1566 62608 : atom[s][a].xd2.push_back(add2);
1567 : }
1568 : }
1569 28 : old_size = aend.index()+1;
1570 : }
1571 14 : }
1572 :
1573 14 : void CS2Backbone::init_types(const PDB &pdb) {
1574 14 : vector<AtomNumber> aa = pdb.getAtomNumbers();
1575 109732 : for(unsigned i=0; i<aa.size(); i++) {
1576 36568 : unsigned frag = pdb.getResidueNumber(aa[i]);
1577 36568 : string fragName = pdb.getResidueName(aa[i]);
1578 36568 : string atom_name = pdb.getAtomName(aa[i]);
1579 36568 : char atom_type = atom_name[0];
1580 36568 : if(isdigit(atom_name[0])) atom_type = atom_name[1];
1581 36568 : res_num.push_back(frag);
1582 36568 : unsigned t = 0;
1583 36568 : if (!isSP2(fragName, atom_name)) {
1584 28728 : if (atom_type == 'C') t = D_C;
1585 21714 : else if (atom_type == 'O') t = D_O;
1586 21350 : else if (atom_type == 'H') t = D_H;
1587 3262 : else if (atom_type == 'N') t = D_N;
1588 126 : else if (atom_type == 'S') t = D_S;
1589 0 : else plumed_merror("CS2Backbone:init_type: unknown atom type!\n");
1590 : } else {
1591 7840 : if (atom_type == 'C') t = D_C2;
1592 3192 : else if (atom_type == 'O') t = D_O2;
1593 0 : else if (atom_type == 'N') t = D_N2;
1594 0 : else plumed_merror("CS2Backbone:init_type: unknown atom type!\n");
1595 : }
1596 36568 : type.push_back(t);
1597 : }
1598 14 : }
1599 :
1600 14 : void CS2Backbone::init_rings(const PDB &pdb) {
1601 :
1602 112 : const string pheTyr_n[] = {"CG","CD1","CE1","CZ","CE2","CD2"};
1603 112 : const string trp1_n[] = {"CD2","CE2","CZ2","CH2","CZ3","CE3"};
1604 98 : const string trp2_n[] = {"CG","CD1","NE1","CE2","CD2"};
1605 98 : const string his_n[] = {"CG","ND1","CD2","CE1","NE2"};
1606 :
1607 14 : vector<string> chains;
1608 14 : pdb.getChainNames( chains );
1609 14 : vector<AtomNumber> allatoms = pdb.getAtomNumbers();
1610 : unsigned old_size=0;
1611 :
1612 112 : for(unsigned s=0; s<atom.size(); s++) {
1613 : AtomNumber astart, aend;
1614 : string errmsg;
1615 28 : pdb.getAtomRange( chains[s], astart, aend, errmsg );
1616 : unsigned atom_offset = astart.index();
1617 7448 : for(unsigned r=0; r<atom[s].size(); r++) {
1618 2464 : string frg = pdb.getResidueName(atom[s][r].fd);
1619 11312 : if(!((frg=="PHE")||(frg=="TYR")||(frg=="TRP")||
1620 6594 : (frg=="HIS")||(frg=="HIP")||(frg=="HID")||
1621 4396 : (frg=="HIE")||(frg=="HSD")||(frg=="HSE")||
1622 : (frg=="HSP"))) continue;
1623 :
1624 266 : vector<AtomNumber> frg_atoms = pdb.getAtomsInResidue(atom[s][r].fd,chains[s]);
1625 :
1626 308 : if(frg=="PHE"||frg=="TYR") {
1627 252 : RingInfo ri;
1628 15708 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1629 5068 : unsigned atm = frg_atoms[a].index()-atom_offset+old_size;
1630 55300 : for(unsigned aa=0; aa<6; aa++) {
1631 79884 : if(pdb.getAtomName(frg_atoms[a])==pheTyr_n[aa]) {
1632 1512 : ri.atom[aa] = atm;
1633 1512 : break;
1634 : }
1635 : }
1636 : }
1637 252 : ri.numAtoms = 6;
1638 252 : if(frg=="PHE") ri.rtype = RingInfo::R_PHE;
1639 252 : if(frg=="TYR") ri.rtype = RingInfo::R_TYR;
1640 252 : ringInfo.push_back(ri);
1641 :
1642 14 : } else if(frg=="TRP") {
1643 : //First ring
1644 14 : RingInfo ri;
1645 1036 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1646 336 : unsigned atm = frg_atoms[a].index()-atom_offset+old_size;
1647 3780 : for(unsigned aa=0; aa<6; aa++) {
1648 5418 : if(pdb.getAtomName(frg_atoms[a])==trp1_n[aa]) {
1649 84 : ri.atom[aa] = atm;
1650 84 : break;
1651 : }
1652 : }
1653 : }
1654 14 : ri.numAtoms = 6;
1655 14 : ri.rtype = RingInfo::R_TRP1;
1656 14 : ringInfo.push_back(ri);
1657 : //Second Ring
1658 14 : RingInfo ri2;
1659 1036 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1660 336 : unsigned atm = frg_atoms[a].index()-atom_offset+old_size;
1661 3276 : for(unsigned aa=0; aa<5; aa++) {
1662 4620 : if(pdb.getAtomName(frg_atoms[a])==trp2_n[aa]) {
1663 70 : ri2.atom[aa] = atm;
1664 70 : break;
1665 : }
1666 : }
1667 : }
1668 14 : ri2.numAtoms = 5;
1669 14 : ri2.rtype = RingInfo::R_TRP2;
1670 14 : ringInfo.push_back(ri2);
1671 :
1672 0 : } else if((frg=="HIS")||(frg=="HIP")||(frg=="HID")||
1673 0 : (frg=="HIE")||(frg=="HSD")||(frg=="HSE")||
1674 : (frg=="HSP")) {//HIS case
1675 0 : RingInfo ri;
1676 0 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1677 0 : unsigned atm = frg_atoms[a].index()-atom_offset+old_size;
1678 0 : for(unsigned aa=0; aa<5; aa++) {
1679 0 : if(pdb.getAtomName(frg_atoms[a])==his_n[aa]) {
1680 0 : ri.atom[aa] = atm;
1681 0 : break;
1682 : }
1683 : }
1684 : }
1685 0 : ri.numAtoms = 5;
1686 0 : ri.rtype = RingInfo::R_HIS;
1687 0 : ringInfo.push_back(ri);
1688 : } else {
1689 0 : plumed_merror("Unkwown Ring Fragment");
1690 : }
1691 : }
1692 28 : old_size = aend.index()+1;
1693 : }
1694 14 : }
1695 :
1696 7336 : CS2Backbone::aa_t CS2Backbone::frag2enum(const string &aa) {
1697 :
1698 : aa_t type = ALA;
1699 7336 : if (aa == "ALA") type = ALA;
1700 6916 : else if (aa == "ARG") type = ARG;
1701 6622 : else if (aa == "ASN") type = ASN;
1702 6244 : else if (aa == "ASP") type = ASP;
1703 5880 : else if (aa == "ASH") type = ASP;
1704 5880 : else if (aa == "CYS") type = CYS;
1705 5712 : else if (aa == "CYM") type = CYS;
1706 5712 : else if (aa == "GLN") type = GLN;
1707 5586 : else if (aa == "GLU") type = GLU;
1708 5096 : else if (aa == "GLH") type = GLU;
1709 5096 : else if (aa == "GLY") type = GLY;
1710 3850 : else if (aa == "HIS") type = HIS;
1711 3850 : else if (aa == "HSE") type = HIS;
1712 3850 : else if (aa == "HIE") type = HIS;
1713 3850 : else if (aa == "HSP") type = HIS;
1714 3850 : else if (aa == "HIP") type = HIS;
1715 3850 : else if (aa == "HSD") type = HIS;
1716 3850 : else if (aa == "HID") type = HIS;
1717 3850 : else if (aa == "ILE") type = ILE;
1718 3430 : else if (aa == "LEU") type = LEU;
1719 3052 : else if (aa == "LYS") type = LYS;
1720 2464 : else if (aa == "MET") type = MET;
1721 2268 : else if (aa == "PHE") type = PHE;
1722 1596 : else if (aa == "PRO") type = PRO;
1723 1302 : else if (aa == "SER") type = SER;
1724 924 : else if (aa == "THR") type = THR;
1725 462 : else if (aa == "TRP") type = TRP;
1726 420 : else if (aa == "TYR") type = TYR;
1727 336 : else if (aa == "VAL") type = VAL;
1728 0 : else if (aa == "UNK") type = UNK;
1729 0 : else plumed_merror("CS2Backbone: Error converting string " + aa + " into amino acid index: not a valid 3-letter code");
1730 7336 : return type;
1731 : }
1732 :
1733 2464 : vector<string> CS2Backbone::side_chain_atoms(const string &s) {
1734 : vector<string> sc;
1735 :
1736 2464 : if(s=="ALA") {
1737 280 : sc.push_back( "CB" );
1738 280 : sc.push_back( "HB1" );
1739 280 : sc.push_back( "HB2" );
1740 280 : sc.push_back( "HB3" );
1741 140 : return sc;
1742 2324 : } else if(s=="ARG") {
1743 196 : sc.push_back( "CB" );
1744 196 : sc.push_back( "CG" );
1745 196 : sc.push_back( "CD" );
1746 196 : sc.push_back( "NE" );
1747 196 : sc.push_back( "CZ" );
1748 196 : sc.push_back( "NH1" );
1749 196 : sc.push_back( "NH2" );
1750 196 : sc.push_back( "NH3" );
1751 196 : sc.push_back( "HB1" );
1752 196 : sc.push_back( "HB2" );
1753 196 : sc.push_back( "HB3" );
1754 196 : sc.push_back( "HG1" );
1755 196 : sc.push_back( "HG2" );
1756 196 : sc.push_back( "HG3" );
1757 196 : sc.push_back( "HD1" );
1758 196 : sc.push_back( "HD2" );
1759 196 : sc.push_back( "HD3" );
1760 196 : sc.push_back( "HE" );
1761 196 : sc.push_back( "HH11" );
1762 196 : sc.push_back( "HH12" );
1763 196 : sc.push_back( "HH21" );
1764 196 : sc.push_back( "HH22" );
1765 196 : sc.push_back( "1HH1" );
1766 196 : sc.push_back( "2HH1" );
1767 196 : sc.push_back( "1HH2" );
1768 196 : sc.push_back( "2HH2" );
1769 98 : return sc;
1770 2226 : } else if(s=="ASN") {
1771 252 : sc.push_back( "CB" );
1772 252 : sc.push_back( "CG" );
1773 252 : sc.push_back( "OD1" );
1774 252 : sc.push_back( "ND2" );
1775 252 : sc.push_back( "HB1" );
1776 252 : sc.push_back( "HB2" );
1777 252 : sc.push_back( "HB3" );
1778 252 : sc.push_back( "HD21" );
1779 252 : sc.push_back( "HD22" );
1780 252 : sc.push_back( "1HD2" );
1781 252 : sc.push_back( "2HD2" );
1782 126 : return sc;
1783 4074 : } else if(s=="ASP"||s=="ASH") {
1784 252 : sc.push_back( "CB" );
1785 252 : sc.push_back( "CG" );
1786 252 : sc.push_back( "OD1" );
1787 252 : sc.push_back( "OD2" );
1788 252 : sc.push_back( "HB1" );
1789 252 : sc.push_back( "HB2" );
1790 252 : sc.push_back( "HB3" );
1791 126 : return sc;
1792 3892 : } else if(s=="CYS"||s=="CYM") {
1793 112 : sc.push_back( "CB" );
1794 112 : sc.push_back( "SG" );
1795 112 : sc.push_back( "HB1" );
1796 112 : sc.push_back( "HB2" );
1797 112 : sc.push_back( "HB3" );
1798 112 : sc.push_back( "HG1" );
1799 112 : sc.push_back( "HG" );
1800 56 : return sc;
1801 1918 : } else if(s=="GLN") {
1802 84 : sc.push_back( "CB" );
1803 84 : sc.push_back( "CG" );
1804 84 : sc.push_back( "CD" );
1805 84 : sc.push_back( "OE1" );
1806 84 : sc.push_back( "NE2" );
1807 84 : sc.push_back( "HB1" );
1808 84 : sc.push_back( "HB2" );
1809 84 : sc.push_back( "HB3" );
1810 84 : sc.push_back( "HG1" );
1811 84 : sc.push_back( "HG2" );
1812 84 : sc.push_back( "HG3" );
1813 84 : sc.push_back( "HE21" );
1814 84 : sc.push_back( "HE22" );
1815 84 : sc.push_back( "1HE2" );
1816 84 : sc.push_back( "2HE2" );
1817 42 : return sc;
1818 3584 : } else if(s=="GLU"||s=="GLH") {
1819 336 : sc.push_back( "CB" );
1820 336 : sc.push_back( "CG" );
1821 336 : sc.push_back( "CD" );
1822 336 : sc.push_back( "OE1" );
1823 336 : sc.push_back( "OE2" );
1824 336 : sc.push_back( "HB1" );
1825 336 : sc.push_back( "HB2" );
1826 336 : sc.push_back( "HB3" );
1827 336 : sc.push_back( "HG1" );
1828 336 : sc.push_back( "HG2" );
1829 336 : sc.push_back( "HG3" );
1830 168 : return sc;
1831 1708 : } else if(s=="GLY") {
1832 840 : sc.push_back( "HA2" );
1833 420 : return sc;
1834 9016 : } else if(s=="HIS"||s=="HSE"||s=="HIE"||s=="HSD"||s=="HID"||s=="HIP"||s=="HSP") {
1835 0 : sc.push_back( "CB" );
1836 0 : sc.push_back( "CG" );
1837 0 : sc.push_back( "ND1" );
1838 0 : sc.push_back( "CD2" );
1839 0 : sc.push_back( "CE1" );
1840 0 : sc.push_back( "NE2" );
1841 0 : sc.push_back( "HB1" );
1842 0 : sc.push_back( "HB2" );
1843 0 : sc.push_back( "HB3" );
1844 0 : sc.push_back( "HD1" );
1845 0 : sc.push_back( "HD2" );
1846 0 : sc.push_back( "HE1" );
1847 0 : sc.push_back( "HE2" );
1848 0 : return sc;
1849 1288 : } else if(s=="ILE") {
1850 280 : sc.push_back( "CB" );
1851 280 : sc.push_back( "CG1" );
1852 280 : sc.push_back( "CG2" );
1853 280 : sc.push_back( "CD" );
1854 280 : sc.push_back( "HB" );
1855 280 : sc.push_back( "HG11" );
1856 280 : sc.push_back( "HG12" );
1857 280 : sc.push_back( "HG21" );
1858 280 : sc.push_back( "HG22" );
1859 280 : sc.push_back( "HG23" );
1860 280 : sc.push_back( "1HG1" );
1861 280 : sc.push_back( "2HG1" );
1862 280 : sc.push_back( "1HG2" );
1863 280 : sc.push_back( "2HG2" );
1864 280 : sc.push_back( "3HG2" );
1865 280 : sc.push_back( "HD1" );
1866 280 : sc.push_back( "HD2" );
1867 280 : sc.push_back( "HD3" );
1868 140 : return sc;
1869 1148 : } else if(s=="LEU") {
1870 252 : sc.push_back( "CB" );
1871 252 : sc.push_back( "CG" );
1872 252 : sc.push_back( "CD1" );
1873 252 : sc.push_back( "CD2" );
1874 252 : sc.push_back( "HB1" );
1875 252 : sc.push_back( "HB2" );
1876 252 : sc.push_back( "HB3" );
1877 252 : sc.push_back( "HG" );
1878 252 : sc.push_back( "HD11" );
1879 252 : sc.push_back( "HD12" );
1880 252 : sc.push_back( "HD13" );
1881 252 : sc.push_back( "HD21" );
1882 252 : sc.push_back( "HD22" );
1883 252 : sc.push_back( "HD23" );
1884 252 : sc.push_back( "1HD1" );
1885 252 : sc.push_back( "2HD1" );
1886 252 : sc.push_back( "3HD1" );
1887 252 : sc.push_back( "1HD2" );
1888 252 : sc.push_back( "2HD2" );
1889 252 : sc.push_back( "3HD2" );
1890 126 : return sc;
1891 1022 : } else if(s=="LYS") {
1892 392 : sc.push_back( "CB" );
1893 392 : sc.push_back( "CG" );
1894 392 : sc.push_back( "CD" );
1895 392 : sc.push_back( "CE" );
1896 392 : sc.push_back( "NZ" );
1897 392 : sc.push_back( "HB1" );
1898 392 : sc.push_back( "HB2" );
1899 392 : sc.push_back( "HB3" );
1900 392 : sc.push_back( "HG1" );
1901 392 : sc.push_back( "HG2" );
1902 392 : sc.push_back( "HG3" );
1903 392 : sc.push_back( "HD1" );
1904 392 : sc.push_back( "HD2" );
1905 392 : sc.push_back( "HD3" );
1906 392 : sc.push_back( "HE1" );
1907 392 : sc.push_back( "HE2" );
1908 392 : sc.push_back( "HE3" );
1909 392 : sc.push_back( "HZ1" );
1910 392 : sc.push_back( "HZ2" );
1911 392 : sc.push_back( "HZ3" );
1912 196 : return sc;
1913 826 : } else if(s=="MET") {
1914 140 : sc.push_back( "CB" );
1915 140 : sc.push_back( "CG" );
1916 140 : sc.push_back( "SD" );
1917 140 : sc.push_back( "CE" );
1918 140 : sc.push_back( "HB1" );
1919 140 : sc.push_back( "HB2" );
1920 140 : sc.push_back( "HB3" );
1921 140 : sc.push_back( "HG1" );
1922 140 : sc.push_back( "HG2" );
1923 140 : sc.push_back( "HG3" );
1924 140 : sc.push_back( "HE1" );
1925 140 : sc.push_back( "HE2" );
1926 140 : sc.push_back( "HE3" );
1927 70 : return sc;
1928 756 : } else if(s=="PHE") {
1929 448 : sc.push_back( "CB" );
1930 448 : sc.push_back( "CG" );
1931 448 : sc.push_back( "CD1" );
1932 448 : sc.push_back( "CD2" );
1933 448 : sc.push_back( "CE1" );
1934 448 : sc.push_back( "CE2" );
1935 448 : sc.push_back( "CZ" );
1936 448 : sc.push_back( "HB1" );
1937 448 : sc.push_back( "HB2" );
1938 448 : sc.push_back( "HB3" );
1939 448 : sc.push_back( "HD1" );
1940 448 : sc.push_back( "HD2" );
1941 448 : sc.push_back( "HD3" );
1942 448 : sc.push_back( "HE1" );
1943 448 : sc.push_back( "HE2" );
1944 448 : sc.push_back( "HE3" );
1945 448 : sc.push_back( "HZ" );
1946 224 : return sc;
1947 532 : } else if(s=="PRO") {
1948 196 : sc.push_back( "CB" );
1949 196 : sc.push_back( "CG" );
1950 196 : sc.push_back( "CD" );
1951 196 : sc.push_back( "HB1" );
1952 196 : sc.push_back( "HB2" );
1953 196 : sc.push_back( "HB3" );
1954 196 : sc.push_back( "HG1" );
1955 196 : sc.push_back( "HG2" );
1956 196 : sc.push_back( "HG3" );
1957 196 : sc.push_back( "HD1" );
1958 196 : sc.push_back( "HD2" );
1959 196 : sc.push_back( "HD3" );
1960 98 : return sc;
1961 434 : } else if(s=="SER") {
1962 252 : sc.push_back( "CB" );
1963 252 : sc.push_back( "OG" );
1964 252 : sc.push_back( "HB1" );
1965 252 : sc.push_back( "HB2" );
1966 252 : sc.push_back( "HB3" );
1967 252 : sc.push_back( "HG1" );
1968 252 : sc.push_back( "HG" );
1969 126 : return sc;
1970 308 : } else if(s=="THR") {
1971 308 : sc.push_back( "CB" );
1972 308 : sc.push_back( "OG1" );
1973 308 : sc.push_back( "CG2" );
1974 308 : sc.push_back( "HB" );
1975 308 : sc.push_back( "HG1" );
1976 308 : sc.push_back( "HG21" );
1977 308 : sc.push_back( "HG22" );
1978 308 : sc.push_back( "HG23" );
1979 308 : sc.push_back( "1HG2" );
1980 308 : sc.push_back( "2HG2" );
1981 308 : sc.push_back( "3HG2" );
1982 154 : return sc;
1983 154 : } else if(s=="TRP") {
1984 28 : sc.push_back( "CB" );
1985 28 : sc.push_back( "CG" );
1986 28 : sc.push_back( "CD1" );
1987 28 : sc.push_back( "CD2" );
1988 28 : sc.push_back( "NE1" );
1989 28 : sc.push_back( "CE2" );
1990 28 : sc.push_back( "CE3" );
1991 28 : sc.push_back( "CZ2" );
1992 28 : sc.push_back( "CZ3" );
1993 28 : sc.push_back( "CH2" );
1994 28 : sc.push_back( "HB1" );
1995 28 : sc.push_back( "HB2" );
1996 28 : sc.push_back( "HB3" );
1997 28 : sc.push_back( "HD1" );
1998 28 : sc.push_back( "HE1" );
1999 28 : sc.push_back( "HE3" );
2000 28 : sc.push_back( "HZ2" );
2001 28 : sc.push_back( "HZ3" );
2002 28 : sc.push_back( "HH2" );
2003 14 : return sc;
2004 140 : } else if(s=="TYR") {
2005 56 : sc.push_back( "CB" );
2006 56 : sc.push_back( "CG" );
2007 56 : sc.push_back( "CD1" );
2008 56 : sc.push_back( "CD2" );
2009 56 : sc.push_back( "CE1" );
2010 56 : sc.push_back( "CE2" );
2011 56 : sc.push_back( "CZ" );
2012 56 : sc.push_back( "OH" );
2013 56 : sc.push_back( "HB1" );
2014 56 : sc.push_back( "HB2" );
2015 56 : sc.push_back( "HB3" );
2016 56 : sc.push_back( "HD1" );
2017 56 : sc.push_back( "HD2" );
2018 56 : sc.push_back( "HD3" );
2019 56 : sc.push_back( "HE1" );
2020 56 : sc.push_back( "HE2" );
2021 56 : sc.push_back( "HE3" );
2022 56 : sc.push_back( "HH" );
2023 28 : return sc;
2024 112 : } else if(s=="VAL") {
2025 224 : sc.push_back( "CB" );
2026 224 : sc.push_back( "CG1" );
2027 224 : sc.push_back( "CG2" );
2028 224 : sc.push_back( "HB" );
2029 224 : sc.push_back( "HG11" );
2030 224 : sc.push_back( "HG12" );
2031 224 : sc.push_back( "HG13" );
2032 224 : sc.push_back( "HG21" );
2033 224 : sc.push_back( "HG22" );
2034 224 : sc.push_back( "HG23" );
2035 224 : sc.push_back( "1HG1" );
2036 224 : sc.push_back( "2HG1" );
2037 224 : sc.push_back( "3HG1" );
2038 224 : sc.push_back( "1HG2" );
2039 224 : sc.push_back( "2HG2" );
2040 224 : sc.push_back( "3HG2" );
2041 112 : return sc;
2042 0 : } else plumed_merror("CS2Backbone: side_chain_atoms unknown");
2043 : }
2044 :
2045 36568 : bool CS2Backbone::isSP2(const string & resType, const string & atomName) {
2046 : bool sp2 = false;
2047 36568 : if (atomName == "C") return true;
2048 34104 : if (atomName == "O") return true;
2049 :
2050 31668 : if(resType == "TRP") {
2051 :
2052 308 : if (atomName == "CG") sp2 = true;
2053 294 : else if (atomName == "CD1") sp2 = true;
2054 280 : else if (atomName == "CD2") sp2 = true;
2055 266 : else if (atomName == "CE2") sp2 = true;
2056 252 : else if (atomName == "CE3") sp2 = true;
2057 238 : else if (atomName == "CZ2") sp2 = true;
2058 224 : else if (atomName == "CZ3") sp2 = true;
2059 210 : else if (atomName == "CH2") sp2 = true;
2060 :
2061 31360 : } else if (resType == "ASP") {
2062 :
2063 1288 : if (atomName == "CG") sp2 = true;
2064 1162 : else if (atomName == "OD1") sp2 = true;
2065 1036 : else if (atomName == "OD2") sp2 = true;
2066 :
2067 30072 : } else if (resType == "GLU") {
2068 :
2069 2212 : if (atomName == "CD") sp2 = true;
2070 2044 : else if (atomName == "OE1") sp2 = true;
2071 1876 : else if (atomName == "OE2") sp2 = true;
2072 :
2073 27860 : } else if (resType == "ARG") {
2074 :
2075 2156 : if (atomName == "CZ") sp2 = true;
2076 :
2077 25704 : } else if (resType == "HIS") {
2078 :
2079 0 : if (atomName == "CG") sp2 = true;
2080 0 : else if (atomName == "ND1") sp2 = true;
2081 0 : else if (atomName == "CD2") sp2 = true;
2082 0 : else if (atomName == "CE1") sp2 = true;
2083 0 : else if (atomName == "NE2") sp2 = true;
2084 :
2085 25704 : } else if (resType == "PHE") {
2086 :
2087 4032 : if (atomName == "CG") sp2 = true;
2088 3808 : else if (atomName == "CD1") sp2 = true;
2089 3584 : else if (atomName == "CD2") sp2 = true;
2090 3360 : else if (atomName == "CE1") sp2 = true;
2091 3136 : else if (atomName == "CE2") sp2 = true;
2092 2912 : else if (atomName == "CZ") sp2 = true;
2093 :
2094 21672 : } else if (resType == "TYR") {
2095 :
2096 532 : if (atomName == "CG") sp2 = true;
2097 504 : else if (atomName == "CD1") sp2 = true;
2098 476 : else if (atomName == "CD2") sp2 = true;
2099 448 : else if (atomName == "CE1") sp2 = true;
2100 420 : else if (atomName == "CE2") sp2 = true;
2101 392 : else if (atomName == "CZ") sp2 = true;
2102 :
2103 21140 : } else if (resType == "ASN") {
2104 :
2105 1512 : if (atomName == "CG") sp2 = true;
2106 1386 : else if (atomName == "OD1") sp2 = true;
2107 :
2108 19628 : } else if (resType == "GLN") {
2109 :
2110 630 : if (atomName == "CD") sp2 = true;
2111 588 : else if (atomName == "OE1") sp2 = true;
2112 :
2113 : }
2114 :
2115 : return sp2;
2116 : }
2117 :
2118 36568 : bool CS2Backbone::is_chi1_cx(const string & frg, const string & atm) {
2119 36568 : if(atm=="CG") return true;
2120 35868 : if((frg == "CYS")&&(atm =="SG")) return true;
2121 72184 : if(((frg == "ILE")||(frg == "VAL"))&&(atm == "CG1")) return true;
2122 36330 : if((frg == "SER")&&(atm == "OG")) return true;
2123 36974 : if((frg == "THR")&&(atm == "OG1")) return true;
2124 :
2125 : return false;
2126 : }
2127 :
2128 14448 : unsigned CS2Backbone::frag_segment(const unsigned p) {
2129 : unsigned s = 0;
2130 31164 : for(unsigned i=0; i<seg_last.size()-1; i++) {
2131 14448 : if(p>seg_last[i]) s = i+1;
2132 : else break;
2133 : }
2134 14448 : return s;
2135 : }
2136 :
2137 14448 : unsigned CS2Backbone::frag_relitive_index(const unsigned p, const unsigned s) {
2138 14448 : if(s==0) return p;
2139 1512 : return p-seg_last[s-1];
2140 : }
2141 :
2142 0 : void CS2Backbone::debug_report() {
2143 : printf("\t CS2Backbone Initialization report: \n");
2144 : printf("\t -------------------------------\n");
2145 0 : printf("\t Number of segments: %u\n", static_cast<unsigned>(atom.size()));
2146 : printf("\t Segments size: ");
2147 0 : for(unsigned i=0; i<atom.size(); i++) {printf("%u ", static_cast<unsigned>(atom[i].size()));} printf("\n");
2148 : printf("\t%8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s \n",
2149 : "Seg","N","AA","Prev","Curr","Next","SC","XD1","XD2","Phi","Psi","Chi1");
2150 0 : for(unsigned i=0; i<atom.size(); i++) {
2151 0 : for(unsigned j=0; j<atom[i].size(); j++) {
2152 0 : printf("\t%8u %8u %8s %8u %8u %8u %8u %8u %8u %8u %8u %8u \n",
2153 : i+1, j+1,
2154 : atom[i][j].res_name.c_str(),
2155 : (unsigned)atom[i][j].prev.size(),
2156 : (unsigned)atom[i][j].curr.size(),
2157 : (unsigned)atom[i][j].next.size(),
2158 : (unsigned)atom[i][j].side_chain.size(),
2159 : (unsigned)atom[i][j].xd1.size(),
2160 : (unsigned)atom[i][j].xd2.size(),
2161 : (unsigned)atom[i][j].phi.size(),
2162 : (unsigned)atom[i][j].psi.size(),
2163 : (unsigned)atom[i][j].chi1.size());
2164 :
2165 0 : for(unsigned k=0; k<atom[i][j].prev.size(); k++) { printf("%8i ", atom[i][j].prev[k]);} printf("\n");
2166 0 : for(unsigned k=0; k<atom[i][j].curr.size(); k++) { printf("%8i ", atom[i][j].curr[k]);} printf("\n");
2167 0 : for(unsigned k=0; k<atom[i][j].next.size(); k++) { printf("%8i ", atom[i][j].next[k]);} printf("\n");
2168 0 : for(unsigned k=0; k<atom[i][j].side_chain.size(); k++) { printf("%8i ", atom[i][j].side_chain[k]);} printf("\n");
2169 0 : for(unsigned k=0; k<atom[i][j].xd1.size(); k++) { printf("%8i ", atom[i][j].xd1[k]);} printf("\n");
2170 0 : for(unsigned k=0; k<atom[i][j].xd2.size(); k++) { printf("%8i ", atom[i][j].xd2[k]);} printf("\n");
2171 0 : for(unsigned k=0; k<atom[i][j].phi.size(); k++) { printf("%8i ", atom[i][j].phi[k]);} printf("\n");
2172 0 : for(unsigned k=0; k<atom[i][j].psi.size(); k++) { printf("%8i ", atom[i][j].psi[k]);} printf("\n");
2173 0 : for(unsigned k=0; k<atom[i][j].chi1.size(); k++) { printf("%8i ", atom[i][j].chi1[k]);} printf("\n");
2174 :
2175 : }
2176 : }
2177 :
2178 : printf("\t Rings: \n");
2179 : printf("\t ------ \n");
2180 0 : printf("\t Number of rings: %u\n", static_cast<unsigned>(ringInfo.size()));
2181 : printf("\t%8s %8s %8s %8s\n", "Num","Type","RType","N.atoms");
2182 0 : for(unsigned i=0; i<ringInfo.size(); i++) {
2183 0 : printf("\t%8u %8u %8u \n",i+1,ringInfo[i].rtype,ringInfo[i].numAtoms);
2184 0 : for(unsigned j=0; j<ringInfo[i].numAtoms; j++) {printf("%8u ", ringInfo[i].atom[j]);} printf("\n");
2185 : }
2186 0 : }
2187 :
2188 816690 : void CS2Backbone::xdist_name_map(string & name) {
2189 1633380 : if((name == "OT1")||(name == "OC1")) name = "O";
2190 2450070 : else if ((name == "HN") || (name == "HT1") || (name == "H1")) name = "H";
2191 1620094 : else if ((name == "CG1")|| (name == "OG")||
2192 1624000 : (name == "SG") || (name == "OG1")) name = "CG";
2193 1595790 : else if ((name == "HA1")|| (name == "HA3")) name = "HA";
2194 816690 : }
2195 :
2196 14 : void CS2Backbone::update() {
2197 : // write status file
2198 28 : if(getWstride()>0&& (getStep()%getWstride()==0 || getCPT()) ) writeStatus();
2199 14 : }
2200 :
2201 : }
2202 4839 : }
|
