Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2023 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.60 // 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 : #include <map>
38 : #include <algorithm>
39 :
40 : #include "MetainferenceBase.h"
41 : #include "core/ActionRegister.h"
42 : #include "tools/Pbc.h"
43 : #include "tools/PDB.h"
44 : #include "tools/Torsion.h"
45 : #include "tools/Communicator.h"
46 :
47 : namespace PLMD {
48 : namespace isdb {
49 :
50 : //+PLUMEDOC ISDB_COLVAR CS2BACKBONE
51 : /*
52 : Calculates the backbone chemical shifts for a protein.
53 :
54 : The functional form is that of CamShift, which is discussed in the first paper cited below. The chemical shift
55 : of the selected nuclei can be saved as components. Alternatively one can calculate either
56 : the CAMSHIFT score (useful as a collective variable as you can see in the second paper cited below or as a scoring
57 : function as discussed in the third paper cited below) or a [METAINFERENCE](METAINFERENCE.md) score (using DOSCORE).
58 : For these two latter cases experimental chemical shifts must be provided.
59 :
60 : CS2BACKBONE calculation can be relatively heavy because it often uses a large number of atoms, it can
61 : be run in parallel using MPI and Openmp (see [here](actions.md) for more details).
62 :
63 : As a general rule, when using [CS2BACKBONE](CS2BACKBONE.md) or other experimental restraints it may be better to
64 : increase the accuracy of the constraint algorithm due to the increased strain on the bonded structure.
65 : In the case of GROMACS it is safer to use lincs-iter=2 and lincs-order=6.
66 :
67 : In general the system for which chemical shifts are calculated must be completely included in
68 : ATOMS and a TEMPLATE pdb file for the same atoms should be provided as well in the folder DATADIR.
69 : The system is made automatically whole unless NOPBC is used, in particular if the system is made
70 : by multiple chains it is usually better to use NOPBC and make the molecule whole [WHOLEMOLECULES](WHOLEMOLECULES.md)
71 : selecting an appropriate order of the atoms. The pdb file is needed to the generate a simple topology of the protein.
72 : For histidine residues in protonation states different from D the HIE/HSE HIP/HSP name should be used.
73 : GLH and ASH can be used for the alternative protonation of GLU and ASP. Non-standard amino acids and other
74 : molecules are not yet supported, but in principle they can be named UNK. If multiple chains are present
75 : the chain identifier must be in the standard PDB format, together with the TER keyword at the end of each chain.
76 : Termini groups like ACE or NME should be removed from the TEMPLATE pdb because they are not recognized by
77 : CS2BACKBONE.
78 :
79 : Atoms indices in the TEMPLATE file should be numbered from 1 to N where N is the number of atoms used in ATOMS.
80 : This is not a problem for simple cases where atoms goes from 1 to N but is instead something to be carefull in case
81 : that a terminal group is removed from the PDB file.
82 :
83 : In addition to a pdb file one needs to provide a list of chemical shifts to be calculated using one
84 : file per nucleus type (CAshifts.dat, CBshifts.dat, Cshifts.dat, Hshifts.dat, HAshifts.dat, Nshifts.dat),
85 : add only the files for the nuclei you need, but each file should include all protein residues.
86 : A chemical shift for a nucleus is calculated if a value greater than 0 is provided.
87 : For practical purposes the value can correspond to the experimental value.
88 : Residues numbers should match that used in the pdb file, but must be positive, so double check the pdb.
89 : The first and last residue of each chain should be preceded by a # character.
90 :
91 : ````
92 : CAshifts.dat:
93 : #1 0.0
94 : 2 55.5
95 : 3 58.4
96 : .
97 : .
98 : #last 0.0
99 : #first of second chain
100 : .
101 : #last of second chain
102 : ````
103 :
104 : The default behavior is to store the values for the active nuclei in components (ca-#, cb-#,
105 : co-#, ha-#, hn-#, nh-# and expca-#, expcb-#, expco-#, expha-#, exphn-#, exp-nh#) with NOEXP it is possible
106 : to only store the back-calculated values, where # includes a chain and residue number.
107 :
108 : One additional file is always needed in the folder DATADIR: camshift.db. This file includes all the parameters needed to
109 : calculate the chemical shifts and can be found in regtest/isdb/rt-cs2backbone/data/ .
110 :
111 : Additional material and examples can be also found in the tutorials as well as in the cs2backbone regtests
112 : in the isdb folder.
113 :
114 : ## Examples
115 :
116 : In this first example the chemical shifts are used to calculate a collective variable to be used
117 : in NMR driven Metadynamics that is similar to what was done in the second paper cited below:
118 :
119 : ```plumed
120 : #SETTINGS
121 : whole: GROUP ATOMS=2612-2514:-1,961-1:-1,2466-962:-1,2513-2467:-1
122 : WHOLEMOLECULES ENTITY0=whole
123 : cs: CS2BACKBONE ...
124 : ATOMS=1-2612
125 : DATADIR=regtest/isdb/rt-cs2backbone/data/
126 : TEMPLATE=template.pdb CAMSHIFT NOPBC
127 : ...
128 :
129 : metad: METAD ARG=cs HEIGHT=0.5 SIGMA=0.1 PACE=200 BIASFACTOR=10
130 : PRINT ARG=cs,metad.bias FILE=COLVAR STRIDE=100
131 : ```
132 :
133 : In this second example the chemical shifts are used as replica-averaged restrained as was done in the fourth and fifth paper cited below.
134 :
135 : ```plumed
136 : #SETTINGS NREPLICAS=2
137 : cs: CS2BACKBONE ATOMS=1-174 DATADIR=regtest/isdb/rt-cs2backbone/data/
138 : encs: ENSEMBLE ARG=(cs\.hn-.*),(cs\.nh-.*)
139 : stcs: STATS ARG=encs.* SQDEVSUM PARARG=(cs\.exphn-.*),(cs\.expnh-.*)
140 : RESTRAINT ARG=stcs.sqdevsum AT=0 KAPPA=0 SLOPE=24
141 :
142 : PRINT ARG=(cs\.hn-.*),(cs\.nh-.*) FILE=RESTRAINT STRIDE=100
143 : ```
144 :
145 : This third example show how to use chemical shifts to calculate a [METAINFERENCE](METAINFERENCE.md) score .
146 :
147 : ```plumed
148 : cs: CS2BACKBONE ...
149 : ATOMS=1-174
150 : DATADIR=regtest/isdb/rt-cs2backbone/data/
151 : SIGMA_MEAN0=1.0 DOSCORE
152 : ...
153 : csbias: BIASVALUE ARG=cs.score
154 :
155 : PRINT ARG=(cs\.hn-.*),(cs\.nh-.*) FILE=CS.dat STRIDE=1000
156 : PRINT ARG=cs.score FILE=BIAS STRIDE=100
157 : ```
158 :
159 : */
160 : //+ENDPLUMEDOC
161 :
162 : class CS2BackboneDB {
163 : enum { STD, GLY, PRO};
164 : enum { HA_ATOM, H_ATOM, N_ATOM, CA_ATOM, CB_ATOM, C_ATOM };
165 : static const unsigned aa_kind = 3;
166 : static const unsigned atm_kind = 6;
167 : static const unsigned numXtraDists = 27;
168 :
169 : // ALA, ARG, ASN, ASP, CYS, GLU, GLN, GLY, HIS, ILE, LEU, LYS, MET, PHE, PRO, SER, THR, TRP, TYR, VAL
170 : double c_aa[aa_kind][atm_kind][20];
171 : double c_aa_prev[aa_kind][atm_kind][20];
172 : double c_aa_succ[aa_kind][atm_kind][20];
173 : double co_bb[aa_kind][atm_kind][16];
174 : double co_sc_[aa_kind][atm_kind][20][20];
175 : double co_xd[aa_kind][atm_kind][numXtraDists];
176 : double co_sphere[aa_kind][atm_kind][2][8];
177 : // for ring current effects
178 : // Phe, Tyr, Trp_1, Trp_2, His
179 : double co_ring[aa_kind][atm_kind][5];
180 : // for dihedral angles
181 : // co * (a * cos(3 * omega + c) + b * cos(omega + d))
182 : double co_da[aa_kind][atm_kind][3];
183 : double pars_da[aa_kind][atm_kind][3][5];
184 :
185 : public:
186 :
187 11406 : inline unsigned kind(const std::string &s) {
188 11406 : if(s=="GLY") {
189 : return GLY;
190 9708 : } else if(s=="PRO") {
191 264 : return PRO;
192 : }
193 : return STD;
194 : }
195 :
196 11406 : inline unsigned atom_kind(const std::string &s) {
197 11406 : if(s=="HA") {
198 : return HA_ATOM;
199 11328 : } else if(s=="H") {
200 : return H_ATOM;
201 8346 : } else if(s=="N") {
202 : return N_ATOM;
203 5364 : } else if(s=="CA") {
204 : return CA_ATOM;
205 2288 : } else if(s=="CB") {
206 : return CB_ATOM;
207 78 : } else if(s=="C") {
208 78 : return C_ATOM;
209 : }
210 : return -1;
211 : }
212 :
213 : unsigned get_numXtraDists() {
214 : return numXtraDists;
215 : }
216 :
217 : //PARAMETERS
218 : inline double * CONSTAACURR(const unsigned a_kind, const unsigned at_kind) {
219 : return c_aa[a_kind][at_kind];
220 : }
221 : inline double * CONSTAANEXT(const unsigned a_kind, const unsigned at_kind) {
222 : return c_aa_succ[a_kind][at_kind];
223 : }
224 : inline double * CONSTAAPREV(const unsigned a_kind, const unsigned at_kind) {
225 : return c_aa_prev[a_kind][at_kind];
226 : }
227 : inline double * CONST_BB2(const unsigned a_kind, const unsigned at_kind) {
228 : return co_bb[a_kind][at_kind];
229 : }
230 : inline double * CONST_SC2(const unsigned a_kind, const unsigned at_kind, unsigned res_type) {
231 : return co_sc_[a_kind][at_kind][res_type];
232 : }
233 : inline double * CONST_XD(const unsigned a_kind, const unsigned at_kind) {
234 : return co_xd[a_kind][at_kind];
235 : }
236 : inline double * CO_SPHERE(const unsigned a_kind, const unsigned at_kind, unsigned exp_type) {
237 : return co_sphere[a_kind][at_kind][exp_type];
238 : }
239 : inline double * CO_RING(const unsigned a_kind, const unsigned at_kind) {
240 : return co_ring[a_kind][at_kind];
241 : }
242 : inline double * CO_DA(const unsigned a_kind, const unsigned at_kind) {
243 : return co_da[a_kind][at_kind];
244 : }
245 : inline double * PARS_DA(const unsigned a_kind, const unsigned at_kind, const unsigned ang_kind) {
246 : return pars_da[a_kind][at_kind][ang_kind];
247 : }
248 :
249 26 : void parse(const std::string &file, const double dscale) {
250 26 : std::ifstream in;
251 26 : in.open(file.c_str());
252 26 : if(!in) {
253 0 : plumed_merror("Unable to open DB file: " + file);
254 : }
255 :
256 : unsigned c_kind = 0;
257 : unsigned c_atom = 0;
258 :
259 104 : for(unsigned i=0; i<3; i++)
260 546 : for(unsigned j=0; j<6; j++) {
261 9828 : for(unsigned k=0; k<20; k++) {
262 9360 : c_aa[i][j][k]=0.;
263 9360 : c_aa_prev[i][j][k]=0.;
264 9360 : c_aa_succ[i][j][k]=0.;
265 196560 : for(unsigned m=0; m<20; m++) {
266 187200 : co_sc_[i][j][k][m]=0.;
267 : }
268 : }
269 7956 : for(unsigned k=0; k<16; k++) {
270 7488 : co_bb[i][j][k]=0.;
271 : }
272 4212 : for(unsigned k=0; k<8; k++) {
273 3744 : co_sphere[i][j][0][k]=0.;
274 3744 : co_sphere[i][j][1][k]=0.;
275 : }
276 1872 : for(unsigned k=0; k<3; k++) {
277 1404 : co_da[i][j][k]=0.;
278 8424 : for(unsigned l=0; l<5; l++) {
279 7020 : pars_da[i][j][k][l]=0.;
280 : }
281 : }
282 2808 : for(unsigned k=0; k<5; k++) {
283 2340 : co_ring[i][j][k]=0.;
284 : }
285 13104 : for(unsigned k=0; k<numXtraDists; k++) {
286 12636 : co_xd[i][j][k]=0.;
287 : }
288 : }
289 :
290 54470 : while(!in.eof()) {
291 : std::string line;
292 54444 : getline(in,line);
293 54444 : if(line.compare(0,1,"#")==0) {
294 24934 : continue;
295 : }
296 29510 : std::vector<std::string> tok;
297 29510 : std::vector<std::string> tmp;
298 59020 : tok = split(line,' ');
299 378196 : for(unsigned q=0; q<tok.size(); q++)
300 348686 : if(tok[q].size()) {
301 314496 : tmp.push_back(tok[q]);
302 : }
303 29510 : tok = tmp;
304 29510 : if(tok.size()==0) {
305 26 : continue;
306 : }
307 29484 : if(tok[0]=="PAR") {
308 468 : c_kind = kind(tok[2]);
309 468 : c_atom = atom_kind(tok[1]);
310 468 : continue;
311 29016 : } else if(tok[0]=="WEIGHT") {
312 468 : continue;
313 28548 : } else if(tok[0]=="FLATBTM") {
314 468 : continue;
315 28080 : } else if (tok[0] == "SCALEHARM") {
316 468 : continue;
317 27612 : } else if (tok[0] == "TANHAMPLI") {
318 468 : continue;
319 27144 : } else if (tok[0] == "ENDHARMON") {
320 468 : continue;
321 26676 : } else if (tok[0] == "MAXRCDEVI") {
322 468 : continue;
323 26208 : } else if (tok[0] == "RANDCOIL") {
324 468 : continue;
325 25740 : } else if (tok[0] == "CONST") {
326 468 : continue;
327 25272 : } else if (tok[0] == "CONSTAA") {
328 468 : assign(c_aa[c_kind][c_atom],tok,1);
329 468 : continue;
330 24804 : } else if (tok[0] == "CONSTAA-1") {
331 468 : assign(c_aa_prev[c_kind][c_atom],tok,1);
332 468 : continue;
333 24336 : } else if (tok[0] == "CONSTAA+1") {
334 468 : assign(c_aa_succ[c_kind][c_atom],tok,1);
335 468 : continue;
336 23868 : } else if (tok[0] == "COBB1") {
337 468 : continue;
338 23400 : } else if (tok[0] == "COBB2") {
339 : //angstrom to nm
340 468 : assign(co_bb[c_kind][c_atom],tok,dscale);
341 468 : continue;
342 22932 : } else if (tok[0] == "SPHERE1") {
343 : // angstrom^-3 to nm^-3
344 468 : assign(co_sphere[c_kind][c_atom][0],tok,1./(dscale*dscale*dscale));
345 468 : continue;
346 22464 : } else if (tok[0] == "SPHERE2") {
347 : //angstrom to nm
348 468 : assign(co_sphere[c_kind][c_atom][1],tok,dscale);
349 468 : continue;
350 21996 : } else if (tok[0] == "DIHEDRALS") {
351 468 : assign(co_da[c_kind][c_atom],tok,1);
352 468 : continue;
353 21528 : } else if (tok[0] == "RINGS") {
354 : // angstrom^-3 to nm^-3
355 468 : assign(co_ring[c_kind][c_atom],tok,1./(dscale*dscale*dscale));
356 2808 : for(unsigned i=1; i<tok.size(); i++) {
357 2340 : co_ring[c_kind][c_atom][i-1] *= 1000;
358 : }
359 468 : continue;
360 21528 : } else if (tok[0] == "HBONDS") {
361 468 : continue;
362 20592 : } else if (tok[0] == "XTRADISTS") {
363 : //angstrom to nm
364 468 : assign(co_xd[c_kind][c_atom],tok,dscale);
365 468 : continue;
366 20124 : } else if(tok[0]=="DIHEDPHI") {
367 468 : assign(pars_da[c_kind][c_atom][0],tok,1);
368 468 : continue;
369 19656 : } else if(tok[0]=="DIHEDPSI") {
370 468 : assign(pars_da[c_kind][c_atom][1],tok,1);
371 468 : continue;
372 19188 : } else if(tok[0]=="DIHEDCHI1") {
373 468 : assign(pars_da[c_kind][c_atom][2],tok,1);
374 468 : continue;
375 : }
376 :
377 : bool ok = false;
378 : const std::string scIdent1 [] = {"COSCALA1", "COSCARG1", "COSCASN1", "COSCASP1", "COSCCYS1", "COSCGLN1", "COSCGLU1",
379 : "COSCGLY1", "COSCHIS1", "COSCILE1", "COSCLEU1", "COSCLYS1", "COSCMET1", "COSCPHE1",
380 : "COSCPRO1", "COSCSER1", "COSCTHR1", "COSCTRP1", "COSCTYR1", "COSCVAL1"
381 393120 : };
382 :
383 294840 : for(unsigned scC = 0; scC < 20; scC++) {
384 285480 : if(tok[0]==scIdent1[scC]) {
385 : ok = true;
386 : break;
387 : }
388 : }
389 18720 : if(ok) {
390 402480 : continue;
391 : }
392 :
393 : const std::string scIdent2 [] = {"COSCALA2", "COSCARG2", "COSCASN2", "COSCASP2", "COSCCYS2", "COSCGLN2", "COSCGLU2",
394 : "COSCGLY2", "COSCHIS2", "COSCILE2", "COSCLEU2", "COSCLYS2", "COSCMET2", "COSCPHE2",
395 : "COSCPRO2", "COSCSER2", "COSCTHR2", "COSCTRP2", "COSCTYR2", "COSCVAL2"
396 196560 : };
397 :
398 98280 : for(unsigned scC = 0; scC < 20; scC++) {
399 98280 : if(tok[0]==scIdent2[scC]) {
400 : //angstrom to nm
401 9360 : assign(co_sc_[c_kind][c_atom][scC],tok,dscale);
402 : ok = true;
403 : break;
404 : }
405 : }
406 9360 : if(ok) {
407 196560 : continue;
408 : }
409 :
410 0 : if(tok.size()) {
411 0 : std::string str_err = "DB WARNING: unrecognized token: " + tok[0];
412 0 : plumed_merror(str_err);
413 : }
414 619190 : }
415 26 : in.close();
416 26 : }
417 :
418 : private:
419 :
420 29510 : std::vector<std::string> &split(const std::string &s, char delim, std::vector<std::string> &elems) {
421 29510 : std::stringstream ss(s);
422 : std::string item;
423 378196 : while (getline(ss, item, delim)) {
424 348686 : elems.push_back(item);
425 : }
426 29510 : return elems;
427 29510 : }
428 :
429 29510 : std::vector<std::string> split(const std::string &s, char delim) {
430 29510 : std::vector<std::string> elems;
431 29510 : split(s, delim, elems);
432 29510 : return elems;
433 0 : }
434 :
435 14976 : void assign(double * f, const std::vector<std::string> & v, const double scale) {
436 177372 : for(unsigned i=1; i<v.size(); i++) {
437 162396 : f[i-1] = scale*(atof(v[i].c_str()));
438 162396 : if(std::abs(f[i-1])<0.000001) {
439 91468 : f[i-1]=0.;
440 : }
441 : }
442 14976 : }
443 : };
444 :
445 : class CS2Backbone : public MetainferenceBase {
446 : struct ChemicalShift {
447 : double exp_cs; // a reference chemical shifts
448 : Value *comp; // a pointer to the component
449 : unsigned res_kind; // residue type (STD/GLY/PRO)
450 : unsigned atm_kind; // nuclues (HA/CA/CB/CO/NH/HN)
451 : unsigned res_type_prev; // previous residue (ALA/VAL/..)
452 : unsigned res_type_curr; // current residue (ALA/VAL/..)
453 : unsigned res_type_next; // next residue (ALA/VAL/..)
454 : std::string res_name; // residue name
455 : std::string nucleus; // chemical shift
456 : bool has_chi1; // does we have a chi1
457 : unsigned csatoms; // fixed number of atoms used
458 : unsigned totcsatoms; // number of atoms used
459 : int res_num; // residue number
460 : unsigned chain; // chain number
461 : unsigned ipos; // index of the atom for which we are calculating the chemical shifts
462 : std::vector<unsigned> bb; // atoms for the previous, current and next backbone
463 : std::vector<unsigned> side_chain;// atoms for the current sidechain
464 : std::vector<int> xd1; // additional couple of atoms
465 : std::vector<int> xd2; // additional couple of atoms
466 : std::vector<unsigned> box_nb; // non-bonded atoms
467 :
468 10938 : ChemicalShift():
469 10938 : exp_cs(0.),
470 10938 : comp(NULL),
471 10938 : res_kind(0),
472 10938 : atm_kind(0),
473 10938 : res_type_prev(0),
474 10938 : res_type_curr(0),
475 10938 : res_type_next(0),
476 10938 : res_name(""),
477 10938 : nucleus(""),
478 10938 : has_chi1(true),
479 10938 : csatoms(0),
480 10938 : totcsatoms(0),
481 10938 : res_num(0),
482 10938 : chain(0),
483 10938 : ipos(0) {
484 10938 : xd1.reserve(26);
485 10938 : xd2.reserve(26);
486 10938 : box_nb.reserve(150);
487 10938 : }
488 : };
489 :
490 : struct RingInfo {
491 : enum {R_PHE, R_TYR, R_TRP1, R_TRP2, R_HIS};
492 : unsigned rtype; // one out of five different types
493 : unsigned atom[6]; // up to six member per ring
494 : unsigned numAtoms; // number of ring members (5 or 6)
495 : Vector position; // center of ring coordinates
496 : Vector normVect; // ring plane normal std::vector
497 : Vector g[6]; // std::vector of the std::vectors used for normVect
498 : double lengthN2; // square of length of normVect
499 : double lengthNV; // length of normVect
500 376 : RingInfo():
501 376 : rtype(0),
502 376 : numAtoms(0),
503 376 : lengthN2(NAN),
504 376 : lengthNV(NAN) {
505 2632 : for(unsigned i=0; i<6; i++) {
506 2256 : atom[i]=0;
507 : }
508 376 : }
509 : };
510 :
511 : enum aa_t {ALA, ARG, ASN, ASP, CYS, GLN, GLU, GLY, HIS, ILE, LEU, LYS, MET, PHE, PRO, SER, THR, TRP, TYR, VAL, UNK};
512 : enum sequence_t {Np, CAp, HAp, Cp, Op, Nc, Hc, CAc, HAc, Cc, Oc, Nn, Hn, CAn, HAn, Cn, CBc, CGc};
513 :
514 : CS2BackboneDB db;
515 : std::vector<ChemicalShift> chemicalshifts;
516 :
517 : std::vector<RingInfo> ringInfo;
518 : std::vector<unsigned> type;
519 : std::vector<int> res_num;
520 : // chain of each atom (an index), used together with res_num
521 : std::vector<unsigned> res_chain;
522 : // position in the ATOMS list of each atom of the template (by AtomNumber::index())
523 : std::vector<int> atom_pos_;
524 : // the derivatives of chemical shift cs are stored in [cs_offset_[cs],cs_offset_[cs+1])
525 : std::vector<unsigned> cs_offset_;
526 : // the atom of each derivative (it changes only when the neighbour list is updated)
527 : std::vector<unsigned> cs_atoms_;
528 : // the chemical shifts computed by this rank are [cs_begin_,cs_end_)
529 : unsigned cs_begin_;
530 : unsigned cs_end_;
531 : unsigned box_nupdate;
532 : unsigned box_count;
533 : bool camshift;
534 : bool pbc;
535 : bool serial;
536 :
537 : void init_cs(const std::string &file, const std::string &k, const PDB &pdb);
538 : void update_neighb();
539 : void set_layout(const unsigned cs);
540 : void get_block(const unsigned r, const unsigned size, unsigned &b, unsigned &e) const;
541 : unsigned atomPos(const AtomNumber &a) const;
542 : void compute_ring_parameters();
543 : void init_types(const PDB &pdb);
544 : void init_rings(const PDB &pdb);
545 : static aa_t frag2enum(const std::string &aa);
546 : static std::vector<std::string> side_chain_atoms(const std::string &s);
547 : static bool isSP2(const std::string & resType, const std::string & atomName);
548 : static bool is_chi1_cx(const std::string & frg, const std::string & atm);
549 : static void xdist_name_map(std::string & name);
550 :
551 : public:
552 :
553 : explicit CS2Backbone(const ActionOptions&);
554 : static void registerKeywords( Keywords& keys );
555 : void calculate() override;
556 : };
557 :
558 : PLUMED_REGISTER_ACTION(CS2Backbone,"CS2BACKBONE")
559 :
560 30 : void CS2Backbone::registerKeywords( Keywords& keys ) {
561 30 : MetainferenceBase::registerKeywords( keys );
562 30 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
563 30 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
564 30 : keys.add("atoms","ATOMS","The atoms to be included in the calculation, e.g. the whole protein.");
565 30 : keys.add("compulsory","DATADIR","data/","The folder with the experimental chemical shifts.");
566 30 : keys.add("compulsory","TEMPLATE","template.pdb","A PDB file of the protein system.");
567 30 : keys.add("compulsory","NEIGH_FREQ","20","Period in step for neighbor list update.");
568 30 : keys.addFlag("CAMSHIFT",false,"Set to TRUE if you to calculate a single CamShift score.");
569 30 : keys.addFlag("NOEXP",false,"Set to TRUE if you don't want to have fixed components with the experimental values.");
570 60 : keys.addOutputComponent("ha","default","scalar","the calculated Ha hydrogen chemical shifts");
571 60 : keys.addOutputComponent("hn","default","scalar","the calculated H hydrogen chemical shifts");
572 60 : keys.addOutputComponent("nh","default","scalar","the calculated N nitrogen chemical shifts");
573 60 : keys.addOutputComponent("ca","default","scalar","the calculated Ca carbon chemical shifts");
574 60 : keys.addOutputComponent("cb","default","scalar","the calculated Cb carbon chemical shifts");
575 60 : keys.addOutputComponent("co","default","scalar","the calculated C' carbon chemical shifts");
576 60 : keys.addOutputComponent("expha","default","scalar","the experimental Ha hydrogen chemical shifts");
577 60 : keys.addOutputComponent("exphn","default","scalar","the experimental H hydrogen chemical shifts");
578 60 : keys.addOutputComponent("expnh","default","scalar","the experimental N nitrogen chemical shifts");
579 60 : keys.addOutputComponent("expca","default","scalar","the experimental Ca carbon chemical shifts");
580 60 : keys.addOutputComponent("expcb","default","scalar","the experimental Cb carbon chemical shifts");
581 60 : keys.addOutputComponent("expco","default","scalar","the experimental C' carbon chemical shifts");
582 60 : keys.setValueDescription("scalar","the backbone chemical shifts");
583 30 : keys.addDOI("10.1021/ja903772t");
584 30 : keys.addDOI("10.1073/pnas.1218350110");
585 30 : keys.addDOI("10.1016/j.str.2010.04.016");
586 30 : keys.addDOI("10.1021/ja210951z");
587 30 : keys.addDOI("10.1021/jp3106666");
588 30 : }
589 :
590 26 : CS2Backbone::CS2Backbone(const ActionOptions&ao):
591 : PLUMED_METAINF_INIT(ao),
592 26 : cs_begin_(0),
593 26 : cs_end_(0),
594 26 : camshift(false),
595 26 : pbc(true),
596 26 : serial(false) {
597 26 : std::vector<AtomNumber> used_atoms;
598 26 : parseAtomList("ATOMS",used_atoms);
599 :
600 26 : parseFlag("CAMSHIFT",camshift);
601 26 : if(camshift&&getDoScore()) {
602 0 : plumed_merror("It is not possible to use CAMSHIFT and DOSCORE at the same time");
603 : }
604 :
605 26 : bool nopbc=!pbc;
606 26 : parseFlag("NOPBC",nopbc);
607 26 : pbc=!nopbc;
608 :
609 26 : parseFlag("SERIAL",serial);
610 :
611 26 : bool noexp=false;
612 52 : parseFlag("NOEXP",noexp);
613 :
614 : std::string stringa_data;
615 52 : parse("DATADIR",stringa_data);
616 :
617 : std::string stringa_template;
618 26 : parse("TEMPLATE",stringa_template);
619 :
620 26 : box_count=0;
621 26 : box_nupdate=20;
622 26 : parse("NEIGH_FREQ", box_nupdate);
623 :
624 26 : std::string stringadb = stringa_data + std::string("/camshift.db");
625 26 : std::string stringapdb = stringa_data + std::string("/") + stringa_template;
626 :
627 : /* Length conversion (parameters are tuned for angstrom) */
628 : double scale=1.;
629 26 : if(!usingNaturalUnits()) {
630 26 : scale = 10.*getUnits().getLength();
631 : }
632 :
633 26 : log.printf(" Initialization of the predictor ...\n");
634 26 : db.parse(stringadb,scale);
635 :
636 26 : PDB pdb;
637 26 : if( !pdb.read(stringapdb,usingNaturalUnits(),1./scale) ) {
638 0 : plumed_merror("missing input file " + stringapdb);
639 : }
640 :
641 : // the template and ATOMS must contain the same atoms (in any order):
642 : // the template atoms are then addressed by their position in ATOMS
643 : {
644 26 : const std::vector<AtomNumber> & tatoms = pdb.getAtomNumbers();
645 26 : if(tatoms.size()!=used_atoms.size()) {
646 0 : plumed_merror("ATOMS and the TEMPLATE " + stringapdb + " must contain the same atoms");
647 : }
648 26 : unsigned maxindex=0;
649 49378 : for(const auto & a : used_atoms) {
650 98678 : maxindex=std::max(maxindex,a.index());
651 : }
652 49378 : for(const auto & a : tatoms) {
653 49352 : maxindex=std::max(maxindex,a.index());
654 : }
655 26 : atom_pos_.assign(maxindex+1,-1);
656 49378 : for(unsigned i=0; i<used_atoms.size(); i++) {
657 49352 : if(atom_pos_[used_atoms[i].index()]!=-1) {
658 0 : plumed_merror("an atom is listed twice in ATOMS");
659 : }
660 49352 : atom_pos_[used_atoms[i].index()]=i;
661 : }
662 49378 : for(const auto & a : tatoms) {
663 49352 : if(atom_pos_[a.index()]==-1) {
664 : std::string num;
665 0 : Tools::convert(a.serial(),num);
666 0 : plumed_merror("atom " + num + " of the TEMPLATE " + stringapdb + " is not in ATOMS");
667 : }
668 : }
669 : }
670 :
671 : // first of all we build the list of chemical shifts we want to predict
672 26 : log.printf(" Reading experimental data ...\n");
673 26 : log.flush();
674 26 : stringadb = stringa_data + std::string("/CAshifts.dat");
675 26 : log.printf(" Initializing CA shifts %s\n", stringadb.c_str());
676 26 : init_cs(stringadb, "CA", pdb);
677 26 : stringadb = stringa_data + std::string("/CBshifts.dat");
678 26 : log.printf(" Initializing CB shifts %s\n", stringadb.c_str());
679 26 : init_cs(stringadb, "CB", pdb);
680 26 : stringadb = stringa_data + std::string("/Cshifts.dat");
681 26 : log.printf(" Initializing C' shifts %s\n", stringadb.c_str());
682 26 : init_cs(stringadb, "C", pdb);
683 26 : stringadb = stringa_data + std::string("/HAshifts.dat");
684 26 : log.printf(" Initializing HA shifts %s\n", stringadb.c_str());
685 26 : init_cs(stringadb, "HA", pdb);
686 26 : stringadb = stringa_data + std::string("/Hshifts.dat");
687 26 : log.printf(" Initializing H shifts %s\n", stringadb.c_str());
688 26 : init_cs(stringadb, "H", pdb);
689 26 : stringadb = stringa_data + std::string("/Nshifts.dat");
690 26 : log.printf(" Initializing N shifts %s\n", stringadb.c_str());
691 52 : init_cs(stringadb, "N", pdb);
692 :
693 26 : if(chemicalshifts.size()==0) {
694 0 : plumed_merror("There are no chemical shifts to calculate, there must be at least a not empty file (CA|CB|C|HA|H|N|shifts.dat)");
695 : }
696 :
697 26 : init_types(pdb);
698 26 : init_rings(pdb);
699 :
700 26 : log<<" Bibliography "
701 26 : <<plumed.cite("Kohlhoff K, Robustelli P, Cavalli A, Salvatella A, Vendruscolo M, J. Am. Chem. Soc. 131, 13894 (2009)");
702 26 : if(camshift) {
703 18 : log<<plumed.cite("Granata D, Camilloni C, Vendruscolo M, Laio A, Proc. Natl. Acad. Sci. USA 110, 6817 (2013)");
704 : } else {
705 34 : log<<plumed.cite("Camilloni C, Robustelli P, De Simone A, Cavalli A, Vendruscolo M, J. Am. Chem. Soc. 134, 3968 (2012)");
706 : }
707 26 : log<<plumed.cite("Bonomi M, Camilloni C, Bioinformatics, 33, 3999 (2017)");
708 26 : log<<"\n";
709 :
710 26 : if(camshift) {
711 9 : noexp = true;
712 9 : addValueWithDerivatives();
713 9 : setNotPeriodic();
714 : } else {
715 7842 : for(unsigned cs=0; cs<chemicalshifts.size(); cs++) {
716 : std::string num;
717 7825 : Tools::convert(chemicalshifts[cs].res_num,num);
718 : std::string chain_num;
719 7825 : Tools::convert(chemicalshifts[cs].chain,chain_num);
720 7825 : if(getDoScore()) {
721 7068 : addComponent(chemicalshifts[cs].nucleus+chain_num+"-"+num);
722 7068 : componentIsNotPeriodic(chemicalshifts[cs].nucleus+chain_num+"-"+num);
723 7068 : chemicalshifts[cs].comp = getPntrToComponent(chemicalshifts[cs].nucleus+chain_num+"-"+num);
724 4712 : setParameter(chemicalshifts[cs].exp_cs);
725 : } else {
726 16407 : addComponentWithDerivatives(chemicalshifts[cs].nucleus+chain_num+"-"+num);
727 16407 : componentIsNotPeriodic(chemicalshifts[cs].nucleus+chain_num+"-"+num);
728 16407 : chemicalshifts[cs].comp = getPntrToComponent(chemicalshifts[cs].nucleus+chain_num+"-"+num);
729 : }
730 : }
731 17 : if(getDoScore()) {
732 4 : Initialise(chemicalshifts.size());
733 : }
734 : }
735 :
736 26 : if(!noexp) {
737 7842 : for(unsigned cs=0; cs<chemicalshifts.size(); cs++) {
738 : std::string num;
739 7825 : Tools::convert(chemicalshifts[cs].res_num,num);
740 : std::string chain_num;
741 7825 : Tools::convert(chemicalshifts[cs].chain,chain_num);
742 23475 : addComponent("exp"+chemicalshifts[cs].nucleus+chain_num+"-"+num);
743 23475 : componentIsNotPeriodic("exp"+chemicalshifts[cs].nucleus+chain_num+"-"+num);
744 23475 : Value* comp=getPntrToComponent("exp"+chemicalshifts[cs].nucleus+chain_num+"-"+num);
745 7825 : comp->set(chemicalshifts[cs].exp_cs);
746 : }
747 : }
748 :
749 26 : requestAtoms(used_atoms, false);
750 26 : setDerivatives();
751 26 : checkRead();
752 52 : }
753 :
754 156 : void CS2Backbone::init_cs(const std::string &file, const std::string &nucl, const PDB &pdb) {
755 : // number of chains
756 156 : std::vector<std::string> chains;
757 156 : pdb.getChainNames( chains );
758 : unsigned ichain=0;
759 :
760 156 : std::ifstream in;
761 156 : in.open(file.c_str());
762 156 : if(!in) {
763 : return;
764 : }
765 156 : std::istream_iterator<std::string> iter(in), end;
766 : unsigned begin=0;
767 :
768 156 : while(iter!=end) {
769 20064 : std::string tok = *iter;
770 : ++iter;
771 20064 : if(tok[0]=='#') {
772 : ++iter;
773 624 : if(begin==1) {
774 : begin=0;
775 312 : ichain++;
776 : } else {
777 : begin=1;
778 : }
779 624 : continue;
780 : }
781 : const int resnum = atoi(tok.c_str());
782 : tok = *iter;
783 : ++iter;
784 : double cs = atof(tok.c_str());
785 19440 : if(cs==0) {
786 8034 : continue;
787 : }
788 :
789 : int fres, lres;
790 : std::string errmsg;
791 11406 : pdb.getResidueRange(chains[ichain], fres, lres, errmsg);
792 11406 : if(resnum==fres||resnum==lres) {
793 0 : plumed_merror("First and Last residue of each chain should be annotated as # in " + file + " Remember that residue numbers should match");
794 : }
795 :
796 : // check in the PDB for the chain/residue/atom and enable the chemical shift
797 11406 : std::string RES = pdb.getResidueName(resnum, chains[ichain]);
798 11406 : if(RES=="HIE"||RES=="HIP"||RES=="HIS"||RES=="HSP"||RES=="HSE"||RES=="CYS"||RES=="GLH"||RES=="ASH"||RES=="UNK") {
799 288 : continue;
800 : }
801 11118 : if(RES=="GLN"&&nucl=="CB") {
802 0 : continue;
803 : }
804 11118 : if(RES=="ILE"&&nucl=="CB") {
805 0 : continue;
806 : }
807 11118 : if(RES=="PRO"&&nucl=="N") {
808 0 : continue;
809 : }
810 11118 : if(RES=="PRO"&&nucl=="H") {
811 0 : continue;
812 : }
813 11118 : if(RES=="PRO"&&nucl=="CB") {
814 108 : continue;
815 : }
816 11010 : if(RES=="GLY"&&nucl=="HA") {
817 0 : continue;
818 : }
819 11010 : if(RES=="GLY"&&nucl=="CB") {
820 72 : continue;
821 : }
822 :
823 10938 : ChemicalShift tmp_cs;
824 :
825 10938 : tmp_cs.exp_cs = cs;
826 10938 : if(nucl=="CA") {
827 : tmp_cs.nucleus = "ca-";
828 7940 : } else if(nucl=="CB") {
829 : tmp_cs.nucleus = "cb-";
830 5808 : } else if(nucl=="C") {
831 : tmp_cs.nucleus = "co-";
832 5808 : } else if(nucl=="HA") {
833 : tmp_cs.nucleus = "ha-";
834 5808 : } else if(nucl=="H") {
835 : tmp_cs.nucleus = "hn-";
836 2904 : } else if(nucl=="N") {
837 : tmp_cs.nucleus = "nh-";
838 : }
839 10938 : tmp_cs.chain = ichain;
840 10938 : tmp_cs.res_num = resnum;
841 10938 : tmp_cs.res_type_curr = frag2enum(RES);
842 10938 : tmp_cs.res_type_prev = frag2enum(pdb.getResidueName(resnum-1, chains[ichain]));
843 10938 : tmp_cs.res_type_next = frag2enum(pdb.getResidueName(resnum+1, chains[ichain]));
844 : tmp_cs.res_name = RES;
845 10938 : tmp_cs.res_kind = db.kind(RES);
846 10938 : tmp_cs.atm_kind = db.atom_kind(nucl);
847 10938 : if(RES!="ALA"&&RES!="GLY") {
848 8684 : tmp_cs.bb.resize(18);
849 8684 : tmp_cs.has_chi1=true;
850 : } else {
851 2254 : tmp_cs.bb.resize(16);
852 2254 : tmp_cs.has_chi1=false;
853 : }
854 :
855 10938 : std::vector<AtomNumber> res_atoms = pdb.getAtomsInResidue(resnum, chains[ichain]);
856 : // find the position of the nucleus and of the other backbone atoms as well as for phi/psi/chi
857 178164 : for(unsigned a=0; a<res_atoms.size(); a++) {
858 167226 : std::string AN = pdb.getAtomName(res_atoms[a]);
859 167226 : if(nucl=="HA"&&(AN=="HA"||AN=="HA1"||AN=="HA3")) {
860 0 : tmp_cs.ipos = atomPos(res_atoms[a]);
861 167226 : } else if(nucl=="H"&&(AN=="H"||AN=="HN")) {
862 2904 : tmp_cs.ipos = atomPos(res_atoms[a]);
863 164322 : } else if(nucl=="N"&&AN=="N") {
864 2904 : tmp_cs.ipos = atomPos(res_atoms[a]);
865 161418 : } else if(nucl=="CA"&&AN=="CA") {
866 2998 : tmp_cs.ipos = atomPos(res_atoms[a]);
867 158420 : } else if(nucl=="CB"&&AN=="CB") {
868 2132 : tmp_cs.ipos = atomPos(res_atoms[a]);
869 156288 : } else if(nucl=="C"&&AN=="C" ) {
870 0 : tmp_cs.ipos = atomPos(res_atoms[a]);
871 : }
872 : }
873 :
874 10938 : std::vector<AtomNumber> prev_res_atoms = pdb.getAtomsInResidue(resnum-1, chains[ichain]);
875 : // find the position of the previous residues backbone atoms
876 173650 : for(unsigned a=0; a<prev_res_atoms.size(); a++) {
877 162712 : std::string AN = pdb.getAtomName(prev_res_atoms[a]);
878 162712 : if(AN=="N") {
879 10938 : tmp_cs.bb[Np] = atomPos(prev_res_atoms[a]);
880 151774 : } else if(AN=="CA") {
881 10938 : tmp_cs.bb[CAp] = atomPos(prev_res_atoms[a]);
882 140836 : } else if(AN=="HA"||AN=="HA1"||AN=="HA3") {
883 10938 : tmp_cs.bb[HAp] = atomPos(prev_res_atoms[a]);
884 129898 : } else if(AN=="C" ) {
885 10938 : tmp_cs.bb[Cp] = atomPos(prev_res_atoms[a]);
886 118960 : } else if(AN=="O" ) {
887 10938 : tmp_cs.bb[Op] = atomPos(prev_res_atoms[a]);
888 : }
889 : }
890 :
891 178164 : for(unsigned a=0; a<res_atoms.size(); a++) {
892 167226 : std::string AN = pdb.getAtomName(res_atoms[a]);
893 167226 : if(AN=="N") {
894 10938 : tmp_cs.bb[Nc] = atomPos(res_atoms[a]);
895 156288 : } else if(AN=="H" ||AN=="HN"||(AN=="CD"&&RES=="PRO")) {
896 10938 : tmp_cs.bb[Hc] = atomPos(res_atoms[a]);
897 145350 : } else if(AN=="CA") {
898 10938 : tmp_cs.bb[CAc] = atomPos(res_atoms[a]);
899 134412 : } else if(AN=="HA"||AN=="HA1"||AN=="HA3") {
900 10938 : tmp_cs.bb[HAc] = atomPos(res_atoms[a]);
901 123474 : } else if(AN=="C" ) {
902 10938 : tmp_cs.bb[Cc] = atomPos(res_atoms[a]);
903 112536 : } else if(AN=="O" ) {
904 10938 : tmp_cs.bb[Oc] = atomPos(res_atoms[a]);
905 : }
906 :
907 167226 : if(RES!="ALA"&&RES!="GLY") {
908 149312 : if(AN=="CB") {
909 8684 : tmp_cs.bb[CBc] = atomPos(res_atoms[a]);
910 : }
911 149312 : if(is_chi1_cx(RES,AN)) {
912 8684 : tmp_cs.bb[CGc] = atomPos(res_atoms[a]);
913 : }
914 : }
915 : }
916 :
917 10938 : std::vector<AtomNumber> next_res_atoms = pdb.getAtomsInResidue(resnum+1, chains[ichain]);
918 10938 : std::string NRES = pdb.getResidueName(resnum+1, chains[ichain]);
919 : // find the position of the previous residues backbone atoms
920 174516 : for(unsigned a=0; a<next_res_atoms.size(); a++) {
921 163578 : std::string AN = pdb.getAtomName(next_res_atoms[a]);
922 163578 : if(AN=="N") {
923 10938 : tmp_cs.bb[Nn] = atomPos(next_res_atoms[a]);
924 152640 : } else if(AN=="H" ||AN=="HN"||(AN=="CD"&&NRES=="PRO")) {
925 10938 : tmp_cs.bb[Hn] = atomPos(next_res_atoms[a]);
926 141702 : } else if(AN=="CA") {
927 10938 : tmp_cs.bb[CAn] = atomPos(next_res_atoms[a]);
928 130764 : } else if(AN=="HA"||AN=="HA1"||AN=="HA3") {
929 10938 : tmp_cs.bb[HAn] = atomPos(next_res_atoms[a]);
930 119826 : } else if(AN=="C" ) {
931 10938 : tmp_cs.bb[Cn] = atomPos(next_res_atoms[a]);
932 : }
933 : }
934 :
935 : // set sidechain atoms
936 10938 : std::vector<std::string> sc_atm = side_chain_atoms(RES);
937 :
938 144892 : for(unsigned sc=0; sc<sc_atm.size(); sc++) {
939 2566264 : for(unsigned aa=0; aa<res_atoms.size(); aa++) {
940 2432310 : if(pdb.getAtomName(res_atoms[aa])==sc_atm[sc]) {
941 101706 : tmp_cs.side_chain.push_back(atomPos(res_atoms[aa]));
942 : }
943 : }
944 : }
945 :
946 : // find atoms for extra distances
947 295326 : const std::string atomsP1[] = {"H", "H", "H", "C", "C", "C", "O", "O", "O", "N", "N", "N", "O", "O", "O", "N", "N", "N", "CG", "CG", "CG", "CG", "CG", "CG", "CG", "CA"};
948 10938 : const int resOffsetP1[] = { 0, 0, 0, -1, -1, -1, 0, 0, 0, 1, 1, 1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, -1};
949 :
950 295326 : const std::string atomsP2[] = {"HA", "C", "CB", "HA", "C", "CB", "HA", "N", "CB", "HA", "N", "CB", "HA", "N", "CB", "HA", "N", "CB", "HA", "N", "C", "C", "N", "CA", "CA", "CA"};
951 10938 : 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};
952 :
953 295326 : for(unsigned q=0; q<db.get_numXtraDists()-1; q++) {
954 284388 : std::vector<AtomNumber> at1;
955 284388 : if(resOffsetP1[q]== 0) {
956 153132 : at1 = res_atoms;
957 : }
958 284388 : if(resOffsetP1[q]==-1) {
959 87504 : at1 = prev_res_atoms;
960 : }
961 284388 : if(resOffsetP1[q]==+1) {
962 43752 : at1 = next_res_atoms;
963 : }
964 :
965 284388 : std::vector<AtomNumber> at2;
966 284388 : if(resOffsetP2[q]== 0) {
967 185946 : at2 = res_atoms;
968 : }
969 284388 : if(resOffsetP2[q]==-1) {
970 76566 : at2 = prev_res_atoms;
971 : }
972 284388 : if(resOffsetP2[q]==+1) {
973 21876 : at2 = next_res_atoms;
974 : }
975 :
976 284388 : int tmp1 = -1;
977 2264686 : for(unsigned a=0; a<at1.size(); a++) {
978 2248044 : std::string name = pdb.getAtomName(at1[a]);
979 2248044 : xdist_name_map(name);
980 :
981 2248044 : if(name==atomsP1[q]) {
982 267746 : tmp1 = atomPos(at1[a]);
983 : break;
984 : }
985 : }
986 :
987 284388 : int tmp2 = -1;
988 1471896 : for(unsigned a=0; a<at2.size(); a++) {
989 1461964 : std::string name = pdb.getAtomName(at2[a]);
990 1461964 : xdist_name_map(name);
991 :
992 1461964 : if(name==atomsP2[q]) {
993 274456 : tmp2 = atomPos(at2[a]);
994 : break;
995 : }
996 : }
997 :
998 284388 : tmp_cs.xd1.push_back(tmp1);
999 284388 : tmp_cs.xd2.push_back(tmp2);
1000 284388 : }
1001 :
1002 : // ready to add a new chemical shifts
1003 10938 : tmp_cs.csatoms = 1 + 16 + tmp_cs.side_chain.size() + 2*tmp_cs.xd1.size();
1004 10938 : if(tmp_cs.res_name!="ALA"&&tmp_cs.res_name!="GLY") {
1005 8684 : tmp_cs.csatoms += 2;
1006 : }
1007 10938 : chemicalshifts.push_back(tmp_cs);
1008 612528 : }
1009 :
1010 156 : in.close();
1011 156 : }
1012 :
1013 : // this assigns an atom-type to each atom of the pdb
1014 26 : void CS2Backbone::init_types(const PDB &pdb) {
1015 : enum atom_t {D_C, D_H, D_N, D_O, D_S, D_C2, D_N2, D_O2};
1016 26 : std::vector<AtomNumber> aa = pdb.getAtomNumbers();
1017 : // these are indexed by the position of the atom in ATOMS
1018 26 : res_num.assign(aa.size(),0);
1019 26 : res_chain.assign(aa.size(),0);
1020 26 : type.assign(aa.size(),0);
1021 : std::map<std::string,unsigned> chain_index;
1022 49378 : for(unsigned i=0; i<aa.size(); i++) {
1023 49352 : const unsigned ipos = atomPos(aa[i]);
1024 49352 : const std::string chain = pdb.getChainID(aa[i]);
1025 49352 : if(chain_index.find(chain)==chain_index.end()) {
1026 52 : const unsigned nchains=chain_index.size();
1027 52 : chain_index[chain]=nchains;
1028 : }
1029 49352 : res_chain[ipos] = chain_index[chain];
1030 49352 : const int frag = pdb.getResidueNumber(aa[i]);
1031 49352 : std::string fragName = pdb.getResidueName(aa[i]);
1032 49352 : std::string atom_name = pdb.getAtomName(aa[i]);
1033 49352 : char atom_type = atom_name[0];
1034 49352 : if(isdigit(atom_name[0])) {
1035 4426 : atom_type = atom_name[1];
1036 : }
1037 49352 : res_num[ipos] = frag;
1038 : unsigned t = 0;
1039 49352 : if (!isSP2(fragName, atom_name)) {
1040 : if (atom_type == 'C') {
1041 : t = D_C;
1042 : } else if (atom_type == 'O') {
1043 : t = D_O;
1044 : } else if (atom_type == 'H') {
1045 : t = D_H;
1046 : } else if (atom_type == 'N') {
1047 : t = D_N;
1048 : } else if (atom_type == 'S') {
1049 : t = D_S;
1050 : } else {
1051 0 : plumed_merror("Unknown atom type: " + atom_name);
1052 : }
1053 : } else {
1054 : if (atom_type == 'C') {
1055 : t = D_C2;
1056 : } else if (atom_type == 'O') {
1057 : t = D_O2;
1058 : } else if (atom_type == 'N') {
1059 : t = D_N2;
1060 : } else {
1061 0 : plumed_merror("Unknown atom type: " + atom_name);
1062 : }
1063 : }
1064 49352 : type[ipos] = t;
1065 : }
1066 26 : }
1067 :
1068 26 : void CS2Backbone::init_rings(const PDB &pdb) {
1069 182 : const std::string pheTyr_n[] = {"CG","CD1","CE1","CZ","CE2","CD2"};
1070 182 : const std::string trp1_n[] = {"CD2","CE2","CZ2","CH2","CZ3","CE3"};
1071 156 : const std::string trp2_n[] = {"CG","CD1","NE1","CE2","CD2"};
1072 156 : const std::string his_n[] = {"CG","ND1","CD2","CE1","NE2"};
1073 :
1074 : // number of chains
1075 26 : std::vector<std::string> chains;
1076 26 : pdb.getChainNames( chains );
1077 : unsigned total_rings_atoms = 0;
1078 :
1079 : // cycle over chains
1080 78 : for(unsigned i=0; i<chains.size(); i++) {
1081 : // cycle over the residues present in the chain, excluding the last one
1082 : // (this works also with negative residue numbers and with gaps in the numbering)
1083 52 : const std::vector<int> chain_res = pdb.getResiduesInChain( chains[i] );
1084 3344 : for(std::size_t ires=0; ires+1<chain_res.size(); ires++) {
1085 3292 : const int res = chain_res[ires];
1086 3292 : std::string frg = pdb.getResidueName(res, chains[i]);
1087 6226 : if(!((frg=="PHE")||(frg=="TYR")||(frg=="TRP")||
1088 2934 : (frg=="HIS")||(frg=="HIP")||(frg=="HID")||
1089 2934 : (frg=="HIE")||(frg=="HSD")||(frg=="HSE")||
1090 2934 : (frg=="HSP"))) {
1091 : continue;
1092 : }
1093 :
1094 358 : std::vector<AtomNumber> frg_atoms = pdb.getAtomsInResidue(res,chains[i]);
1095 :
1096 358 : if(frg=="PHE"||frg=="TYR") {
1097 340 : RingInfo ri;
1098 7176 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1099 6836 : unsigned atm = atomPos(frg_atoms[a]);
1100 40712 : for(unsigned aa=0; aa<6; aa++) {
1101 35916 : if(pdb.getAtomName(frg_atoms[a])==pheTyr_n[aa]) {
1102 2040 : ri.atom[aa] = atm;
1103 2040 : break;
1104 : }
1105 : }
1106 : }
1107 340 : ri.numAtoms = 6;
1108 340 : total_rings_atoms += 6;
1109 340 : if(frg=="PHE") {
1110 304 : ri.rtype = RingInfo::R_PHE;
1111 : }
1112 340 : if(frg=="TYR") {
1113 36 : ri.rtype = RingInfo::R_TYR;
1114 : }
1115 340 : ringInfo.push_back(ri);
1116 :
1117 18 : } else if(frg=="TRP") {
1118 : //First ring
1119 18 : RingInfo ri;
1120 450 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1121 432 : unsigned atm = atomPos(frg_atoms[a]);
1122 2646 : for(unsigned aa=0; aa<6; aa++) {
1123 2322 : if(pdb.getAtomName(frg_atoms[a])==trp1_n[aa]) {
1124 108 : ri.atom[aa] = atm;
1125 108 : break;
1126 : }
1127 : }
1128 : }
1129 18 : ri.numAtoms = 6;
1130 : total_rings_atoms += 6;
1131 18 : ri.rtype = RingInfo::R_TRP1;
1132 18 : ringInfo.push_back(ri);
1133 : //Second Ring
1134 18 : RingInfo ri2;
1135 450 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1136 432 : unsigned atm = atomPos(frg_atoms[a]);
1137 2322 : for(unsigned aa=0; aa<5; aa++) {
1138 1980 : if(pdb.getAtomName(frg_atoms[a])==trp2_n[aa]) {
1139 90 : ri2.atom[aa] = atm;
1140 90 : break;
1141 : }
1142 : }
1143 : }
1144 18 : ri2.numAtoms = 5;
1145 18 : total_rings_atoms += 3;
1146 18 : ri2.rtype = RingInfo::R_TRP2;
1147 18 : ringInfo.push_back(ri2);
1148 :
1149 0 : } else if((frg=="HIS")||(frg=="HIP")||(frg=="HID")||
1150 0 : (frg=="HIE")||(frg=="HSD")||(frg=="HSE")||
1151 0 : (frg=="HSP")) {//HIS case
1152 0 : RingInfo ri;
1153 0 : for(unsigned a=0; a<frg_atoms.size(); a++) {
1154 0 : unsigned atm = atomPos(frg_atoms[a]);
1155 0 : for(unsigned aa=0; aa<5; aa++) {
1156 0 : if(pdb.getAtomName(frg_atoms[a])==his_n[aa]) {
1157 0 : ri.atom[aa] = atm;
1158 0 : break;
1159 : }
1160 : }
1161 : }
1162 0 : ri.numAtoms = 5;
1163 0 : total_rings_atoms += 3;
1164 0 : ri.rtype = RingInfo::R_HIS;
1165 0 : ringInfo.push_back(ri);
1166 : } else {
1167 0 : plumed_merror("Unknown Ring Fragment: " + frg);
1168 : }
1169 358 : }
1170 52 : }
1171 :
1172 10964 : for(unsigned cs=0; cs<chemicalshifts.size(); cs++) {
1173 10938 : chemicalshifts[cs].csatoms += total_rings_atoms;
1174 : }
1175 702 : }
1176 :
1177 26 : void CS2Backbone::calculate() {
1178 26 : if(pbc) {
1179 0 : makeWhole();
1180 : }
1181 26 : if(getExchangeStep()) {
1182 0 : box_count=0;
1183 : }
1184 26 : if(box_count==0) {
1185 26 : update_neighb();
1186 : }
1187 26 : compute_ring_parameters();
1188 :
1189 26 : std::vector<double> camshift_sigma2(6);
1190 26 : camshift_sigma2[0] = 0.08; // HA
1191 26 : camshift_sigma2[1] = 0.30; // HN
1192 26 : camshift_sigma2[2] = 9.00; // NH
1193 26 : camshift_sigma2[3] = 1.30; // CA
1194 26 : camshift_sigma2[4] = 1.56; // CB
1195 26 : camshift_sigma2[5] = 1.70; // CO
1196 :
1197 : // derivatives of each chemical shift, with the layout set in update_neighb
1198 26 : std::vector<Vector> cs_derivs(cs_offset_.back(),Vector(0,0,0));
1199 26 : std::vector<Vector> aa_derivs;
1200 26 : std::vector<double> all_shifts(chemicalshifts.size(),0);
1201 26 : if(camshift||getDoScore()) {
1202 13 : aa_derivs.resize(getNumberOfAtoms(),Vector(0,0,0));
1203 : }
1204 :
1205 26 : const unsigned stride = serial ? 1 : comm.Get_size();
1206 26 : const unsigned rank = serial ? 0 : comm.Get_rank();
1207 :
1208 26 : unsigned nt=OpenMP::getNumThreads();
1209 26 : if(static_cast<std::size_t>(nt)*stride*2>chemicalshifts.size()) {
1210 : nt=1;
1211 : }
1212 :
1213 : // a single loop over the chemical shifts of this rank
1214 26 : #pragma omp parallel num_threads(nt)
1215 : {
1216 : #pragma omp for schedule(dynamic)
1217 : for(unsigned cs=cs_begin_; cs<cs_end_; cs++) {
1218 : const unsigned kdx=cs_offset_[cs];
1219 : const ChemicalShift *myfrag = &chemicalshifts[cs];
1220 : const unsigned aa_kind = myfrag->res_kind;
1221 : const unsigned at_kind = myfrag->atm_kind;
1222 :
1223 : double shift = db.CONSTAAPREV(aa_kind,at_kind)[myfrag->res_type_prev] +
1224 : db.CONSTAACURR(aa_kind,at_kind)[myfrag->res_type_curr] +
1225 : db.CONSTAANEXT(aa_kind,at_kind)[myfrag->res_type_next];
1226 :
1227 : const unsigned ipos = myfrag->ipos;
1228 : unsigned atom_counter = 1;
1229 :
1230 : //BACKBONE (PREV CURR NEXT)
1231 : const double * CONST_BB2 = db.CONST_BB2(aa_kind,at_kind);
1232 : const unsigned bbsize = 16;
1233 : for(unsigned q=0; q<bbsize; q++) {
1234 : const double cb2q = CONST_BB2[q];
1235 : if(cb2q==0.) {
1236 : continue;
1237 : }
1238 : const unsigned jpos = myfrag->bb[q];
1239 : if(ipos==jpos) {
1240 : continue;
1241 : }
1242 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
1243 : const double d = distance.modulo();
1244 : const double fact = cb2q/d;
1245 :
1246 : shift += cb2q*d;
1247 : const Vector der = fact*distance;
1248 :
1249 : cs_derivs[kdx+0] += der;
1250 : cs_derivs[kdx+q+atom_counter] = -der;
1251 : }
1252 :
1253 : atom_counter += bbsize;
1254 :
1255 : //DIHEDRAL ANGLES
1256 : const double *CO_DA = db.CO_DA(aa_kind,at_kind);
1257 : //Phi
1258 : {
1259 : const Vector d0 = delta(getPosition(myfrag->bb[Nc]), getPosition(myfrag->bb[Cp]));
1260 : const Vector d1 = delta(getPosition(myfrag->bb[CAc]), getPosition(myfrag->bb[Nc]));
1261 : const Vector d2 = delta(getPosition(myfrag->bb[Cc]), getPosition(myfrag->bb[CAc]));
1262 : Torsion t;
1263 : Vector dd0, dd1, dd2;
1264 : const double t_phi = t.compute(d0,d1,d2,dd0,dd1,dd2);
1265 : const double *PARS_DA = db.PARS_DA(aa_kind,at_kind,0);
1266 : const double val1 = 3.*t_phi+PARS_DA[3];
1267 : const double val2 = t_phi+PARS_DA[4];
1268 : shift += CO_DA[0]*(PARS_DA[0]*std::cos(val1)+PARS_DA[1]*std::cos(val2)+PARS_DA[2]);
1269 : const double fact = -CO_DA[0]*(+3.*PARS_DA[0]*std::sin(val1)+PARS_DA[1]*std::sin(val2));
1270 :
1271 : cs_derivs[kdx+Cp+1] += fact*dd0;
1272 : cs_derivs[kdx+Nc+1] += fact*(dd1-dd0);
1273 : cs_derivs[kdx+CAc+1]+= fact*(dd2-dd1);
1274 : cs_derivs[kdx+Cc+1] += -fact*dd2;
1275 : }
1276 :
1277 : //Psi
1278 : {
1279 : const Vector d0 = delta(getPosition(myfrag->bb[CAc]), getPosition(myfrag->bb[Nc]));
1280 : const Vector d1 = delta(getPosition(myfrag->bb[Cc]), getPosition(myfrag->bb[CAc]));
1281 : const Vector d2 = delta(getPosition(myfrag->bb[Nn]), getPosition(myfrag->bb[Cc]));
1282 : Torsion t;
1283 : Vector dd0, dd1, dd2;
1284 : const double t_psi = t.compute(d0,d1,d2,dd0,dd1,dd2);
1285 : const double *PARS_DA = db.PARS_DA(aa_kind,at_kind,1);
1286 : const double val1 = 3.*t_psi+PARS_DA[3];
1287 : const double val2 = t_psi+PARS_DA[4];
1288 : shift += CO_DA[1]*(PARS_DA[0]*std::cos(val1)+PARS_DA[1]*std::cos(val2)+PARS_DA[2]);
1289 : const double fact = -CO_DA[1]*(+3.*PARS_DA[0]*std::sin(val1)+PARS_DA[1]*std::sin(val2));
1290 :
1291 : cs_derivs[kdx+Nc+1] += fact*dd0;
1292 : cs_derivs[kdx+CAc+1] += fact*(dd1-dd0);
1293 : cs_derivs[kdx+Cc+1] += fact*(dd2-dd1);
1294 : cs_derivs[kdx+Nn+1] += -fact*dd2;
1295 : }
1296 :
1297 : //Chi
1298 : if(myfrag->has_chi1) {
1299 : const Vector d0 = delta(getPosition(myfrag->bb[CAc]), getPosition(myfrag->bb[Nc]));
1300 : const Vector d1 = delta(getPosition(myfrag->bb[CBc]), getPosition(myfrag->bb[CAc]));
1301 : const Vector d2 = delta(getPosition(myfrag->bb[CGc]), getPosition(myfrag->bb[CBc]));
1302 : Torsion t;
1303 : Vector dd0, dd1, dd2;
1304 : const double t_chi1 = t.compute(d0,d1,d2,dd0,dd1,dd2);
1305 : const double *PARS_DA = db.PARS_DA(aa_kind,at_kind,2);
1306 : const double val1 = 3.*t_chi1+PARS_DA[3];
1307 : const double val2 = t_chi1+PARS_DA[4];
1308 : shift += CO_DA[2]*(PARS_DA[0]*std::cos(val1)+PARS_DA[1]*std::cos(val2)+PARS_DA[2]);
1309 : const double fact = -CO_DA[2]*(+3.*PARS_DA[0]*std::sin(val1)+PARS_DA[1]*std::sin(val2));
1310 :
1311 : cs_derivs[kdx+Nc+1] += fact*dd0;
1312 : cs_derivs[kdx+CAc+1] += fact*(dd1-dd0);
1313 : cs_derivs[kdx+CBc+1] += fact*(dd2-dd1);
1314 : cs_derivs[kdx+CGc+1] += -fact*dd2;
1315 :
1316 : atom_counter += 2;
1317 : }
1318 : //END OF DIHE
1319 :
1320 : //SIDE CHAIN
1321 : const double * CONST_SC2 = db.CONST_SC2(aa_kind,at_kind,myfrag->res_type_curr);
1322 : const unsigned sidsize = myfrag->side_chain.size();
1323 : for(unsigned q=0; q<sidsize; q++) {
1324 : const double cs2q = CONST_SC2[q];
1325 : if(cs2q==0.) {
1326 : continue;
1327 : }
1328 : const unsigned jpos = myfrag->side_chain[q];
1329 : if(ipos==jpos) {
1330 : continue;
1331 : }
1332 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
1333 : const double d = distance.modulo();
1334 : const double fact = cs2q/d;
1335 :
1336 : shift += cs2q*d;
1337 : const Vector der = fact*distance;
1338 : cs_derivs[kdx+0] += der;
1339 : cs_derivs[kdx+q+atom_counter] = -der;
1340 : }
1341 :
1342 : atom_counter += sidsize;
1343 :
1344 : //EXTRA DIST
1345 : const double * CONST_XD = db.CONST_XD(aa_kind,at_kind);
1346 : const unsigned xdsize=myfrag->xd1.size();
1347 : for(unsigned q=0; q<xdsize; q++) {
1348 : const double cxdq = CONST_XD[q];
1349 : if(cxdq==0.) {
1350 : continue;
1351 : }
1352 : if(myfrag->xd1[q]==-1||myfrag->xd2[q]==-1) {
1353 : continue;
1354 : }
1355 : const Vector distance = delta(getPosition(myfrag->xd1[q]),getPosition(myfrag->xd2[q]));
1356 : const double d = distance.modulo();
1357 : const double fact = cxdq/d;
1358 :
1359 : shift += cxdq*d;
1360 : const Vector der = fact*distance;
1361 : cs_derivs[kdx+2*q+atom_counter ] = der;
1362 : cs_derivs[kdx+2*q+atom_counter+1] = -der;
1363 : }
1364 :
1365 : atom_counter += 2*xdsize;
1366 :
1367 : //RINGS
1368 : const double *rc = db.CO_RING(aa_kind,at_kind);
1369 : const unsigned rsize = ringInfo.size();
1370 : // cycle over the list of rings
1371 : for(unsigned q=0; q<rsize; q++) {
1372 : // compute angle from ring middle point to current atom position
1373 : // get distance std::vector from query atom to ring center and normal std::vector to ring plane
1374 : const Vector n = ringInfo[q].normVect;
1375 : const double nL = ringInfo[q].lengthNV;
1376 : const double inL2 = ringInfo[q].lengthN2;
1377 :
1378 : const Vector d = delta(ringInfo[q].position, getPosition(ipos));
1379 : const double dL2 = d.modulo2();
1380 : double dL = std::sqrt(dL2);
1381 : const double idL3 = 1./(dL2*dL);
1382 :
1383 : const double dn = dotProduct(d,n);
1384 : const double dn2 = dn*dn;
1385 : const double dLnL = dL*nL;
1386 : const double dL_nL = dL/nL;
1387 :
1388 : const double ang2 = dn2*inL2/dL2;
1389 : const double u = 1.-3.*ang2;
1390 : const double cc = rc[ringInfo[q].rtype];
1391 :
1392 : shift += cc*u*idL3;
1393 :
1394 : const double fUU = -6.*dn*inL2;
1395 : const double fUQ = fUU/dL;
1396 : const Vector gradUQ = fUQ*(dL2*n - dn*d);
1397 : const Vector gradVQ = (3.*dL*u)*d;
1398 :
1399 : const double fact = cc*idL3*idL3;
1400 : cs_derivs[kdx+0] += fact*(gradUQ - gradVQ);
1401 :
1402 : const double fU = fUU/nL;
1403 : double OneOverN = 1./6.;
1404 : if(ringInfo[q].numAtoms==5) {
1405 : OneOverN=1./3.;
1406 : }
1407 : const Vector factor2 = OneOverN*n;
1408 : const Vector factor4 = (OneOverN/dL_nL)*d;
1409 :
1410 : const Vector gradV = -OneOverN*gradVQ;
1411 :
1412 : if(ringInfo[q].numAtoms==6) {
1413 : // update forces on ring atoms
1414 : for(unsigned at=0; at<6; at++) {
1415 : const Vector ab = crossProduct(d,ringInfo[q].g[at]);
1416 : const Vector c = crossProduct(n,ringInfo[q].g[at]);
1417 : const Vector factor3 = 0.5*dL_nL*c;
1418 : const Vector factor1 = 0.5*ab;
1419 : const Vector gradU = fU*( dLnL*(factor1 - factor2) -dn*(factor3 - factor4) );
1420 : cs_derivs[kdx+at+atom_counter] = fact*(gradU - gradV);
1421 : }
1422 : atom_counter += 6;
1423 : } else {
1424 : for(unsigned at=0; at<3; at++) {
1425 : const Vector ab = crossProduct(d,ringInfo[q].g[at]);
1426 : const Vector c = crossProduct(n,ringInfo[q].g[at]);
1427 : const Vector factor3 = dL_nL*c;
1428 : const Vector factor1 = ab;
1429 : const Vector gradU = fU*( dLnL*(factor1 - factor2) -dn*(factor3 - factor4) );
1430 : cs_derivs[kdx+at+atom_counter] = fact*(gradU - gradV);
1431 : }
1432 : atom_counter += 3;
1433 : }
1434 : }
1435 : //END OF RINGS
1436 :
1437 : //NON BOND
1438 : const double * CONST_CO_SPHERE3 = db.CO_SPHERE(aa_kind,at_kind,0);
1439 : const double * CONST_CO_SPHERE = db.CO_SPHERE(aa_kind,at_kind,1);
1440 : const unsigned boxsize = myfrag->box_nb.size();
1441 : for(unsigned q=0; q<boxsize; q++) {
1442 : const unsigned jpos = myfrag->box_nb[q];
1443 : const Vector distance = delta(getPosition(jpos),getPosition(ipos));
1444 : const double d2 = distance.modulo2();
1445 :
1446 : if(d2<cutOffDist2) {
1447 : double factor1 = std::sqrt(d2);
1448 : double dfactor1 = 1./factor1;
1449 : double factor3 = dfactor1*dfactor1*dfactor1;
1450 : double dfactor3 = -3.*factor3*dfactor1*dfactor1;
1451 :
1452 : if(d2>cutOnDist2) {
1453 : const double af = cutOffDist2 - d2;
1454 : const double bf = cutOffDist2 - 3.*cutOnDist2 + 2.*d2;
1455 : const double cf = invswitch*af;
1456 : const double df = cf*af*bf;
1457 : factor1 *= df;
1458 : factor3 *= df;
1459 :
1460 : const double d4 = d2*d2;
1461 : const double af1 = 15.*cutOnDist2*d2;
1462 : const double bf1 = -14.*d4;
1463 : const double cf1 = -3.*cutOffDist2*cutOnDist2 + cutOffDist2*d2;
1464 : const double df1 = af1+bf1+cf1;
1465 : dfactor1 *= cf*(cutOffDist4+df1);
1466 :
1467 : const double af3 = +2.*cutOffDist2*cutOnDist2;
1468 : const double bf3 = d2*(cutOffDist2+cutOnDist2);
1469 : const double cf3 = -2.*d4;
1470 : const double df3 = (af3+bf3+cf3)*d2;
1471 : dfactor3 *= invswitch*(cutMixed+df3);
1472 : }
1473 :
1474 : const unsigned t = type[jpos];
1475 : shift += factor1*CONST_CO_SPHERE[t] + factor3*CONST_CO_SPHERE3[t] ;
1476 : const double fact = dfactor1*CONST_CO_SPHERE[t]+dfactor3*CONST_CO_SPHERE3[t];
1477 : const Vector der = fact*distance;
1478 :
1479 : cs_derivs[kdx+0] += der;
1480 : cs_derivs[kdx+q+atom_counter] = -der;
1481 : }
1482 : }
1483 : //END NON BOND
1484 :
1485 : atom_counter += boxsize;
1486 : all_shifts[cs] = shift;
1487 : }
1488 : }
1489 :
1490 26 : ++box_count;
1491 26 : if(box_count == box_nupdate) {
1492 26 : box_count = 0;
1493 : }
1494 :
1495 : // every rank owns a contiguous block of chemical shifts: gather the blocks
1496 : // (the atoms of the derivatives are known on every rank)
1497 26 : if(!serial && stride>1) {
1498 24 : const unsigned size=comm.Get_size();
1499 24 : std::vector<int> counts(size), displ(size);
1500 72 : for(unsigned r=0; r<size; r++) {
1501 : unsigned b, e;
1502 48 : get_block(r,size,b,e);
1503 48 : counts[r]=e-b;
1504 48 : displ[r]=b;
1505 : }
1506 24 : const std::vector<double> my_shifts(all_shifts.begin()+cs_begin_, all_shifts.begin()+cs_end_);
1507 24 : comm.Allgatherv(my_shifts.data(), cs_end_-cs_begin_, all_shifts.data(), counts.data(), displ.data());
1508 24 : if(!camshift && !getDoScore()) {
1509 36 : for(unsigned r=0; r<size; r++) {
1510 : unsigned b, e;
1511 24 : get_block(r,size,b,e);
1512 24 : counts[r]=3*(cs_offset_[e]-cs_offset_[b]);
1513 24 : displ[r]=3*cs_offset_[b];
1514 : }
1515 12 : const std::vector<Vector> my_derivs(cs_derivs.begin()+cs_offset_[cs_begin_], cs_derivs.begin()+cs_offset_[cs_end_]);
1516 12 : comm.Allgatherv(&my_derivs[0][0], 3*(cs_offset_[cs_end_]-cs_offset_[cs_begin_]), &cs_derivs[0][0], counts.data(), displ.data());
1517 12 : }
1518 24 : }
1519 :
1520 26 : if(!camshift) {
1521 7842 : for(unsigned cs=0; cs<chemicalshifts.size(); cs++) {
1522 7825 : Value *comp = chemicalshifts[cs].comp;
1523 7825 : comp->set(all_shifts[cs]);
1524 7825 : if(getDoScore()) {
1525 2356 : setCalcData(cs, all_shifts[cs]);
1526 : } else {
1527 5469 : Tensor csvirial;
1528 1292323 : for(unsigned i=cs_offset_[cs]; i<cs_offset_[cs+1]; i++) {
1529 1286854 : setAtomsDerivatives(comp,cs_atoms_[i],cs_derivs[i]);
1530 2573708 : csvirial-=Tensor(getPosition(cs_atoms_[i]),cs_derivs[i]);
1531 : }
1532 5469 : setBoxDerivatives(comp,csvirial);
1533 : }
1534 : }
1535 17 : if(!getDoScore()) {
1536 : return;
1537 : }
1538 : }
1539 :
1540 13 : double score = 0.;
1541 :
1542 : /* Metainference */
1543 13 : if(getDoScore()) {
1544 4 : score = getScore();
1545 1182 : for(unsigned cs=cs_begin_; cs<cs_end_; cs++) {
1546 : const double fact = getMetaDer(cs);
1547 282215 : for(unsigned i=cs_offset_[cs]; i<cs_offset_[cs+1]; i++) {
1548 281037 : aa_derivs[cs_atoms_[i]] += cs_derivs[i]*fact;
1549 : }
1550 : }
1551 : }
1552 :
1553 : /* camshift */
1554 13 : if(camshift) {
1555 1860 : for(unsigned cs=cs_begin_; cs<cs_end_; cs++) {
1556 1851 : score += (all_shifts[cs] - chemicalshifts[cs].exp_cs)*(all_shifts[cs] - chemicalshifts[cs].exp_cs)/camshift_sigma2[chemicalshifts[cs].atm_kind];
1557 1851 : double fact = 2.0*(all_shifts[cs] - chemicalshifts[cs].exp_cs)/camshift_sigma2[chemicalshifts[cs].atm_kind];
1558 434580 : for(unsigned i=cs_offset_[cs]; i<cs_offset_[cs+1]; i++) {
1559 432729 : aa_derivs[cs_atoms_[i]] += cs_derivs[i]*fact;
1560 : }
1561 : }
1562 : }
1563 :
1564 13 : if(!serial) {
1565 13 : comm.Sum(&aa_derivs[0][0], 3*aa_derivs.size());
1566 13 : if(camshift) {
1567 9 : comm.Sum(&score, 1);
1568 : }
1569 : }
1570 :
1571 13 : Tensor virial;
1572 13657 : for(unsigned i=rank; i<getNumberOfAtoms(); i+=stride) {
1573 27288 : virial += Tensor(getPosition(i), aa_derivs[i]);
1574 : }
1575 :
1576 13 : if(!serial) {
1577 13 : comm.Sum(&virial[0][0], 9);
1578 : }
1579 :
1580 : /* calculate final derivatives */
1581 : Value* val;
1582 13 : if(getDoScore()) {
1583 4 : val=getPntrToComponent("score");
1584 4 : setScore(score);
1585 : } else {
1586 : val=getPntrToValue();
1587 9 : setValue(score);
1588 : }
1589 :
1590 : /* at this point we cycle over all atoms */
1591 24689 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
1592 24676 : setAtomsDerivatives(val, i, aa_derivs[i]);
1593 : }
1594 13 : setBoxDerivatives(val,-virial);
1595 26 : }
1596 :
1597 904274 : unsigned CS2Backbone::atomPos(const AtomNumber &a) const {
1598 904274 : plumed_assert(a.index()<atom_pos_.size() && atom_pos_[a.index()]>=0);
1599 904274 : return atom_pos_[a.index()];
1600 : }
1601 :
1602 98 : void CS2Backbone::get_block(const unsigned r, const unsigned size, unsigned &b, unsigned &e) const {
1603 : // contiguous blocks of chemical shifts with about the same number of derivatives
1604 98 : const unsigned ncs = chemicalshifts.size();
1605 98 : const unsigned long long total = cs_offset_.back();
1606 196 : auto bound=[&](const unsigned k) {
1607 196 : if(k==0) {
1608 : return 0U;
1609 : }
1610 146 : if(k>=size) {
1611 50 : return ncs;
1612 : }
1613 96 : const unsigned long long target = (total*k)/size;
1614 96 : const unsigned pos = std::lower_bound(cs_offset_.begin(), cs_offset_.end(), target) - cs_offset_.begin();
1615 96 : return std::min(pos, ncs);
1616 98 : };
1617 98 : b = bound(r);
1618 98 : e = bound(r+1);
1619 98 : }
1620 :
1621 10938 : void CS2Backbone::set_layout(const unsigned cs) {
1622 : // the atom of each derivative, in the same order used in calculate()
1623 10938 : const ChemicalShift & myfrag = chemicalshifts[cs];
1624 10938 : unsigned * atoms = &cs_atoms_[cs_offset_[cs]];
1625 10938 : atoms[0] = myfrag.ipos;
1626 : // backbone distances and dihedrals (slots 1-16)
1627 185946 : for(unsigned q=0; q<16; q++) {
1628 175008 : atoms[q+1] = myfrag.bb[q];
1629 : }
1630 : unsigned atom_counter = 17;
1631 10938 : if(myfrag.has_chi1) {
1632 8684 : atoms[CBc+1] = myfrag.bb[CBc];
1633 8684 : atoms[CGc+1] = myfrag.bb[CGc];
1634 : atom_counter += 2;
1635 : }
1636 112644 : for(unsigned q=0; q<myfrag.side_chain.size(); q++) {
1637 101706 : atoms[q+atom_counter] = myfrag.side_chain[q];
1638 : }
1639 10938 : atom_counter += myfrag.side_chain.size();
1640 295326 : for(unsigned q=0; q<myfrag.xd1.size(); q++) {
1641 284388 : if(myfrag.xd1[q]==-1||myfrag.xd2[q]==-1) {
1642 26574 : continue;
1643 : }
1644 257814 : atoms[2*q+atom_counter] = myfrag.xd2[q];
1645 257814 : atoms[2*q+atom_counter+1] = myfrag.xd1[q];
1646 : }
1647 10938 : atom_counter += 2*myfrag.xd1.size();
1648 223650 : for(unsigned q=0; q<ringInfo.size(); q++) {
1649 212712 : if(ringInfo[q].numAtoms==6) {
1650 1414770 : for(unsigned at=0; at<6; at++) {
1651 1212660 : atoms[at+atom_counter] = ringInfo[q].atom[at];
1652 : }
1653 202110 : atom_counter += 6;
1654 : } else {
1655 10602 : atoms[atom_counter] = ringInfo[q].atom[0];
1656 10602 : atoms[atom_counter+1] = ringInfo[q].atom[2];
1657 10602 : atoms[atom_counter+2] = ringInfo[q].atom[3];
1658 10602 : atom_counter += 3;
1659 : }
1660 : }
1661 466384 : for(unsigned q=0; q<myfrag.box_nb.size(); q++) {
1662 455446 : atoms[q+atom_counter] = myfrag.box_nb[q];
1663 : }
1664 10938 : atom_counter += myfrag.box_nb.size();
1665 10938 : plumed_assert(atom_counter==myfrag.totcsatoms);
1666 10938 : }
1667 :
1668 26 : void CS2Backbone::update_neighb() {
1669 : // cycle over chemical shifts
1670 : #ifdef _OPENMP
1671 : //nt is unused if openmp is not declared
1672 26 : unsigned nt=OpenMP::getNumThreads();
1673 : #endif //_OPENMP
1674 26 : #pragma omp parallel for num_threads(nt)
1675 : for(unsigned cs=0; cs<chemicalshifts.size(); cs++) {
1676 : const unsigned boxsize = getNumberOfAtoms();
1677 : chemicalshifts[cs].box_nb.clear();
1678 : chemicalshifts[cs].box_nb.reserve(150);
1679 : const unsigned ipos = chemicalshifts[cs].ipos;
1680 : const int res_curr = res_num[ipos];
1681 : for(unsigned bat=0; bat<boxsize; bat++) {
1682 : // atoms of the same and of the adjacent residues of the same chain are excluded
1683 : if(res_chain[bat]==res_chain[ipos] && std::abs(res_curr-res_num[bat])<2) {
1684 : continue;
1685 : }
1686 : const Vector distance = delta(getPosition(bat),getPosition(ipos));
1687 : const double d2=distance.modulo2();
1688 : if(d2<cutOffNB2) {
1689 : chemicalshifts[cs].box_nb.push_back(bat);
1690 : }
1691 : }
1692 : chemicalshifts[cs].totcsatoms = chemicalshifts[cs].csatoms + chemicalshifts[cs].box_nb.size();
1693 : }
1694 : // layout of the derivatives
1695 26 : cs_offset_.assign(chemicalshifts.size()+1,0);
1696 10964 : for(unsigned cs=0; cs<chemicalshifts.size(); cs++) {
1697 10938 : cs_offset_[cs+1] = cs_offset_[cs] + chemicalshifts[cs].totcsatoms;
1698 : }
1699 26 : cs_atoms_.assign(cs_offset_.back(),0);
1700 26 : #pragma omp parallel for num_threads(nt)
1701 : for(unsigned cs=0; cs<chemicalshifts.size(); cs++) {
1702 : set_layout(cs);
1703 : }
1704 : // chemical shifts computed by this rank
1705 26 : if(serial) {
1706 0 : cs_begin_ = 0;
1707 0 : cs_end_ = chemicalshifts.size();
1708 : } else {
1709 26 : get_block(comm.Get_rank(),comm.Get_size(),cs_begin_,cs_end_);
1710 : }
1711 26 : }
1712 :
1713 26 : void CS2Backbone::compute_ring_parameters() {
1714 402 : for(unsigned i=0; i<ringInfo.size(); i++) {
1715 376 : const unsigned size = ringInfo[i].numAtoms;
1716 376 : if(size==6) {
1717 358 : ringInfo[i].g[0] = delta(getPosition(ringInfo[i].atom[4]),getPosition(ringInfo[i].atom[2]));
1718 358 : ringInfo[i].g[1] = delta(getPosition(ringInfo[i].atom[5]),getPosition(ringInfo[i].atom[3]));
1719 358 : ringInfo[i].g[2] = delta(getPosition(ringInfo[i].atom[0]),getPosition(ringInfo[i].atom[4]));
1720 358 : ringInfo[i].g[3] = delta(getPosition(ringInfo[i].atom[1]),getPosition(ringInfo[i].atom[5]));
1721 358 : ringInfo[i].g[4] = delta(getPosition(ringInfo[i].atom[2]),getPosition(ringInfo[i].atom[0]));
1722 358 : ringInfo[i].g[5] = delta(getPosition(ringInfo[i].atom[3]),getPosition(ringInfo[i].atom[1]));
1723 : // ring center
1724 358 : Vector midP = getPosition(ringInfo[i].atom[0]);
1725 2148 : for(unsigned j=1; j<size; j++) {
1726 1790 : midP += getPosition(ringInfo[i].atom[j]);
1727 : }
1728 358 : ringInfo[i].position = midP/6.;
1729 : // compute normal std::vector to plane
1730 358 : Vector n1 = crossProduct(ringInfo[i].g[2], -ringInfo[i].g[4]);
1731 358 : Vector n2 = crossProduct(ringInfo[i].g[5], -ringInfo[i].g[1]);
1732 358 : ringInfo[i].normVect = 0.5*(n1 + n2);
1733 : } else {
1734 18 : ringInfo[i].g[0] = delta(getPosition(ringInfo[i].atom[3]),getPosition(ringInfo[i].atom[2]));
1735 18 : ringInfo[i].g[1] = delta(getPosition(ringInfo[i].atom[0]),getPosition(ringInfo[i].atom[3]));
1736 18 : ringInfo[i].g[2] = delta(getPosition(ringInfo[i].atom[2]),getPosition(ringInfo[i].atom[0]));
1737 : // ring center
1738 18 : ringInfo[i].position = (getPosition(ringInfo[i].atom[0])+getPosition(ringInfo[i].atom[2])+getPosition(ringInfo[i].atom[3]))/3.;
1739 : // ring plane normal std::vector
1740 18 : ringInfo[i].normVect = crossProduct(ringInfo[i].g[1],-ringInfo[i].g[2]);
1741 :
1742 : }
1743 : // calculate squared length and length of normal std::vector
1744 376 : ringInfo[i].lengthN2 = 1./ringInfo[i].normVect.modulo2();
1745 376 : ringInfo[i].lengthNV = 1./std::sqrt(ringInfo[i].lengthN2);
1746 : }
1747 26 : }
1748 :
1749 32814 : CS2Backbone::aa_t CS2Backbone::frag2enum(const std::string &aa) {
1750 : aa_t type = ALA;
1751 32814 : if (aa == "ALA") {
1752 : type = ALA;
1753 30798 : } else if (aa == "ARG") {
1754 : type = ARG;
1755 29252 : } else if (aa == "ASN") {
1756 : type = ASN;
1757 27614 : } else if (aa == "ASP") {
1758 : type = ASP;
1759 25924 : } else if (aa == "ASH") {
1760 : type = ASP;
1761 25924 : } else if (aa == "CYS") {
1762 : type = CYS;
1763 25402 : } else if (aa == "CYM") {
1764 : type = CYS;
1765 25402 : } else if (aa == "GLN") {
1766 : type = GLN;
1767 24934 : } else if (aa == "GLU") {
1768 : type = GLU;
1769 22612 : } else if (aa == "GLH") {
1770 : type = GLU;
1771 22612 : } else if (aa == "GLY") {
1772 : type = GLY;
1773 17326 : } else if (aa == "HIS") {
1774 : type = HIS;
1775 17326 : } else if (aa == "HSE") {
1776 : type = HIS;
1777 17326 : } else if (aa == "HIE") {
1778 : type = HIS;
1779 17326 : } else if (aa == "HSP") {
1780 : type = HIS;
1781 17326 : } else if (aa == "HIP") {
1782 : type = HIS;
1783 17326 : } else if (aa == "HSD") {
1784 : type = HIS;
1785 17326 : } else if (aa == "HID") {
1786 : type = HIS;
1787 17326 : } else if (aa == "ILE") {
1788 : type = ILE;
1789 15526 : } else if (aa == "LEU") {
1790 : type = LEU;
1791 13744 : } else if (aa == "LYS") {
1792 : type = LYS;
1793 10954 : } else if (aa == "MET") {
1794 : type = MET;
1795 10144 : } else if (aa == "PHE") {
1796 : type = PHE;
1797 7004 : } else if (aa == "PRO") {
1798 : type = PRO;
1799 6020 : } else if (aa == "SER") {
1800 : type = SER;
1801 4302 : } else if (aa == "THR") {
1802 : type = THR;
1803 2196 : } else if (aa == "TRP") {
1804 : type = TRP;
1805 1980 : } else if (aa == "TYR") {
1806 : type = TYR;
1807 1602 : } else if (aa == "VAL") {
1808 : type = VAL;
1809 0 : } else if (aa == "UNK") {
1810 : type = UNK;
1811 : } else {
1812 0 : plumed_merror("Error converting std::string " + aa + " into amino acid index: not a valid 3-letter code");
1813 : }
1814 32814 : return type;
1815 : }
1816 :
1817 10938 : std::vector<std::string> CS2Backbone::side_chain_atoms(const std::string &s) {
1818 10938 : std::vector<std::string> sc;
1819 :
1820 10938 : if(s=="ALA") {
1821 712 : sc.push_back( "CB" );
1822 712 : sc.push_back( "HB1" );
1823 712 : sc.push_back( "HB2" );
1824 712 : sc.push_back( "HB3" );
1825 712 : return sc;
1826 10226 : } else if(s=="ARG") {
1827 500 : sc.push_back( "CB" );
1828 500 : sc.push_back( "CG" );
1829 500 : sc.push_back( "CD" );
1830 500 : sc.push_back( "NE" );
1831 500 : sc.push_back( "CZ" );
1832 500 : sc.push_back( "NH1" );
1833 500 : sc.push_back( "NH2" );
1834 500 : sc.push_back( "NH3" );
1835 500 : sc.push_back( "HB1" );
1836 500 : sc.push_back( "HB2" );
1837 500 : sc.push_back( "HB3" );
1838 500 : sc.push_back( "HG1" );
1839 500 : sc.push_back( "HG2" );
1840 500 : sc.push_back( "HG3" );
1841 500 : sc.push_back( "HD1" );
1842 500 : sc.push_back( "HD2" );
1843 500 : sc.push_back( "HD3" );
1844 500 : sc.push_back( "HE" );
1845 500 : sc.push_back( "HH11" );
1846 500 : sc.push_back( "HH12" );
1847 500 : sc.push_back( "HH21" );
1848 500 : sc.push_back( "HH22" );
1849 500 : sc.push_back( "1HH1" );
1850 500 : sc.push_back( "2HH1" );
1851 500 : sc.push_back( "1HH2" );
1852 500 : sc.push_back( "2HH2" );
1853 500 : return sc;
1854 9726 : } else if(s=="ASN") {
1855 594 : sc.push_back( "CB" );
1856 594 : sc.push_back( "CG" );
1857 594 : sc.push_back( "OD1" );
1858 594 : sc.push_back( "ND2" );
1859 594 : sc.push_back( "HB1" );
1860 594 : sc.push_back( "HB2" );
1861 594 : sc.push_back( "HB3" );
1862 594 : sc.push_back( "HD21" );
1863 594 : sc.push_back( "HD22" );
1864 594 : sc.push_back( "1HD2" );
1865 594 : sc.push_back( "2HD2" );
1866 594 : return sc;
1867 9132 : } else if(s=="ASP"||s=="ASH") {
1868 606 : sc.push_back( "CB" );
1869 606 : sc.push_back( "CG" );
1870 606 : sc.push_back( "OD1" );
1871 606 : sc.push_back( "OD2" );
1872 606 : sc.push_back( "HB1" );
1873 606 : sc.push_back( "HB2" );
1874 606 : sc.push_back( "HB3" );
1875 606 : return sc;
1876 8526 : } else if(s=="CYS"||s=="CYM") {
1877 0 : sc.push_back( "CB" );
1878 0 : sc.push_back( "SG" );
1879 0 : sc.push_back( "HB1" );
1880 0 : sc.push_back( "HB2" );
1881 0 : sc.push_back( "HB3" );
1882 0 : sc.push_back( "HG1" );
1883 0 : sc.push_back( "HG" );
1884 0 : return sc;
1885 8526 : } else if(s=="GLN") {
1886 162 : sc.push_back( "CB" );
1887 162 : sc.push_back( "CG" );
1888 162 : sc.push_back( "CD" );
1889 162 : sc.push_back( "OE1" );
1890 162 : sc.push_back( "NE2" );
1891 162 : sc.push_back( "HB1" );
1892 162 : sc.push_back( "HB2" );
1893 162 : sc.push_back( "HB3" );
1894 162 : sc.push_back( "HG1" );
1895 162 : sc.push_back( "HG2" );
1896 162 : sc.push_back( "HG3" );
1897 162 : sc.push_back( "HE21" );
1898 162 : sc.push_back( "HE22" );
1899 162 : sc.push_back( "1HE2" );
1900 162 : sc.push_back( "2HE2" );
1901 162 : return sc;
1902 8364 : } else if(s=="GLU"||s=="GLH") {
1903 792 : sc.push_back( "CB" );
1904 792 : sc.push_back( "CG" );
1905 792 : sc.push_back( "CD" );
1906 792 : sc.push_back( "OE1" );
1907 792 : sc.push_back( "OE2" );
1908 792 : sc.push_back( "HB1" );
1909 792 : sc.push_back( "HB2" );
1910 792 : sc.push_back( "HB3" );
1911 792 : sc.push_back( "HG1" );
1912 792 : sc.push_back( "HG2" );
1913 792 : sc.push_back( "HG3" );
1914 792 : return sc;
1915 7572 : } else if(s=="GLY") {
1916 1542 : sc.push_back( "HA2" );
1917 1542 : return sc;
1918 6030 : } else if(s=="HIS"||s=="HSE"||s=="HIE"||s=="HSD"||s=="HID"||s=="HIP"||s=="HSP") {
1919 0 : sc.push_back( "CB" );
1920 0 : sc.push_back( "CG" );
1921 0 : sc.push_back( "ND1" );
1922 0 : sc.push_back( "CD2" );
1923 0 : sc.push_back( "CE1" );
1924 0 : sc.push_back( "NE2" );
1925 0 : sc.push_back( "HB1" );
1926 0 : sc.push_back( "HB2" );
1927 0 : sc.push_back( "HB3" );
1928 0 : sc.push_back( "HD1" );
1929 0 : sc.push_back( "HD2" );
1930 0 : sc.push_back( "HE1" );
1931 0 : sc.push_back( "HE2" );
1932 0 : return sc;
1933 6030 : } else if(s=="ILE") {
1934 540 : sc.push_back( "CB" );
1935 540 : sc.push_back( "CG1" );
1936 540 : sc.push_back( "CG2" );
1937 540 : sc.push_back( "CD" );
1938 540 : sc.push_back( "HB" );
1939 540 : sc.push_back( "HG11" );
1940 540 : sc.push_back( "HG12" );
1941 540 : sc.push_back( "HG21" );
1942 540 : sc.push_back( "HG22" );
1943 540 : sc.push_back( "HG23" );
1944 540 : sc.push_back( "1HG1" );
1945 540 : sc.push_back( "2HG1" );
1946 540 : sc.push_back( "1HG2" );
1947 540 : sc.push_back( "2HG2" );
1948 540 : sc.push_back( "3HG2" );
1949 540 : sc.push_back( "HD1" );
1950 540 : sc.push_back( "HD2" );
1951 540 : sc.push_back( "HD3" );
1952 540 : return sc;
1953 5490 : } else if(s=="LEU") {
1954 712 : sc.push_back( "CB" );
1955 712 : sc.push_back( "CG" );
1956 712 : sc.push_back( "CD1" );
1957 712 : sc.push_back( "CD2" );
1958 712 : sc.push_back( "HB1" );
1959 712 : sc.push_back( "HB2" );
1960 712 : sc.push_back( "HB3" );
1961 712 : sc.push_back( "HG" );
1962 712 : sc.push_back( "HD11" );
1963 712 : sc.push_back( "HD12" );
1964 712 : sc.push_back( "HD13" );
1965 712 : sc.push_back( "HD21" );
1966 712 : sc.push_back( "HD22" );
1967 712 : sc.push_back( "HD23" );
1968 712 : sc.push_back( "1HD1" );
1969 712 : sc.push_back( "2HD1" );
1970 712 : sc.push_back( "3HD1" );
1971 712 : sc.push_back( "1HD2" );
1972 712 : sc.push_back( "2HD2" );
1973 712 : sc.push_back( "3HD2" );
1974 712 : return sc;
1975 4778 : } else if(s=="LYS") {
1976 1008 : sc.push_back( "CB" );
1977 1008 : sc.push_back( "CG" );
1978 1008 : sc.push_back( "CD" );
1979 1008 : sc.push_back( "CE" );
1980 1008 : sc.push_back( "NZ" );
1981 1008 : sc.push_back( "HB1" );
1982 1008 : sc.push_back( "HB2" );
1983 1008 : sc.push_back( "HB3" );
1984 1008 : sc.push_back( "HG1" );
1985 1008 : sc.push_back( "HG2" );
1986 1008 : sc.push_back( "HG3" );
1987 1008 : sc.push_back( "HD1" );
1988 1008 : sc.push_back( "HD2" );
1989 1008 : sc.push_back( "HD3" );
1990 1008 : sc.push_back( "HE1" );
1991 1008 : sc.push_back( "HE2" );
1992 1008 : sc.push_back( "HE3" );
1993 1008 : sc.push_back( "HZ1" );
1994 1008 : sc.push_back( "HZ2" );
1995 1008 : sc.push_back( "HZ3" );
1996 1008 : return sc;
1997 3770 : } else if(s=="MET") {
1998 288 : sc.push_back( "CB" );
1999 288 : sc.push_back( "CG" );
2000 288 : sc.push_back( "SD" );
2001 288 : sc.push_back( "CE" );
2002 288 : sc.push_back( "HB1" );
2003 288 : sc.push_back( "HB2" );
2004 288 : sc.push_back( "HB3" );
2005 288 : sc.push_back( "HG1" );
2006 288 : sc.push_back( "HG2" );
2007 288 : sc.push_back( "HG3" );
2008 288 : sc.push_back( "HE1" );
2009 288 : sc.push_back( "HE2" );
2010 288 : sc.push_back( "HE3" );
2011 288 : return sc;
2012 3482 : } else if(s=="PHE") {
2013 1114 : sc.push_back( "CB" );
2014 1114 : sc.push_back( "CG" );
2015 1114 : sc.push_back( "CD1" );
2016 1114 : sc.push_back( "CD2" );
2017 1114 : sc.push_back( "CE1" );
2018 1114 : sc.push_back( "CE2" );
2019 1114 : sc.push_back( "CZ" );
2020 1114 : sc.push_back( "HB1" );
2021 1114 : sc.push_back( "HB2" );
2022 1114 : sc.push_back( "HB3" );
2023 1114 : sc.push_back( "HD1" );
2024 1114 : sc.push_back( "HD2" );
2025 1114 : sc.push_back( "HD3" );
2026 1114 : sc.push_back( "HE1" );
2027 1114 : sc.push_back( "HE2" );
2028 1114 : sc.push_back( "HE3" );
2029 1114 : sc.push_back( "HZ" );
2030 1114 : return sc;
2031 2368 : } else if(s=="PRO") {
2032 108 : sc.push_back( "CB" );
2033 108 : sc.push_back( "CG" );
2034 108 : sc.push_back( "CD" );
2035 108 : sc.push_back( "HB1" );
2036 108 : sc.push_back( "HB2" );
2037 108 : sc.push_back( "HB3" );
2038 108 : sc.push_back( "HG1" );
2039 108 : sc.push_back( "HG2" );
2040 108 : sc.push_back( "HG3" );
2041 108 : sc.push_back( "HD1" );
2042 108 : sc.push_back( "HD2" );
2043 108 : sc.push_back( "HD3" );
2044 108 : return sc;
2045 2260 : } else if(s=="SER") {
2046 694 : sc.push_back( "CB" );
2047 694 : sc.push_back( "OG" );
2048 694 : sc.push_back( "HB1" );
2049 694 : sc.push_back( "HB2" );
2050 694 : sc.push_back( "HB3" );
2051 694 : sc.push_back( "HG1" );
2052 694 : sc.push_back( "HG" );
2053 694 : return sc;
2054 1566 : } else if(s=="THR") {
2055 774 : sc.push_back( "CB" );
2056 774 : sc.push_back( "OG1" );
2057 774 : sc.push_back( "CG2" );
2058 774 : sc.push_back( "HB" );
2059 774 : sc.push_back( "HG1" );
2060 774 : sc.push_back( "HG21" );
2061 774 : sc.push_back( "HG22" );
2062 774 : sc.push_back( "HG23" );
2063 774 : sc.push_back( "1HG2" );
2064 774 : sc.push_back( "2HG2" );
2065 774 : sc.push_back( "3HG2" );
2066 774 : return sc;
2067 792 : } else if(s=="TRP") {
2068 72 : sc.push_back( "CB" );
2069 72 : sc.push_back( "CG" );
2070 72 : sc.push_back( "CD1" );
2071 72 : sc.push_back( "CD2" );
2072 72 : sc.push_back( "NE1" );
2073 72 : sc.push_back( "CE2" );
2074 72 : sc.push_back( "CE3" );
2075 72 : sc.push_back( "CZ2" );
2076 72 : sc.push_back( "CZ3" );
2077 72 : sc.push_back( "CH2" );
2078 72 : sc.push_back( "HB1" );
2079 72 : sc.push_back( "HB2" );
2080 72 : sc.push_back( "HB3" );
2081 72 : sc.push_back( "HD1" );
2082 72 : sc.push_back( "HE1" );
2083 72 : sc.push_back( "HE3" );
2084 72 : sc.push_back( "HZ2" );
2085 72 : sc.push_back( "HZ3" );
2086 72 : sc.push_back( "HH2" );
2087 72 : return sc;
2088 720 : } else if(s=="TYR") {
2089 144 : sc.push_back( "CB" );
2090 144 : sc.push_back( "CG" );
2091 144 : sc.push_back( "CD1" );
2092 144 : sc.push_back( "CD2" );
2093 144 : sc.push_back( "CE1" );
2094 144 : sc.push_back( "CE2" );
2095 144 : sc.push_back( "CZ" );
2096 144 : sc.push_back( "OH" );
2097 144 : sc.push_back( "HB1" );
2098 144 : sc.push_back( "HB2" );
2099 144 : sc.push_back( "HB3" );
2100 144 : sc.push_back( "HD1" );
2101 144 : sc.push_back( "HD2" );
2102 144 : sc.push_back( "HD3" );
2103 144 : sc.push_back( "HE1" );
2104 144 : sc.push_back( "HE2" );
2105 144 : sc.push_back( "HE3" );
2106 144 : sc.push_back( "HH" );
2107 144 : return sc;
2108 576 : } else if(s=="VAL") {
2109 576 : sc.push_back( "CB" );
2110 576 : sc.push_back( "CG1" );
2111 576 : sc.push_back( "CG2" );
2112 576 : sc.push_back( "HB" );
2113 576 : sc.push_back( "HG11" );
2114 576 : sc.push_back( "HG12" );
2115 576 : sc.push_back( "HG13" );
2116 576 : sc.push_back( "HG21" );
2117 576 : sc.push_back( "HG22" );
2118 576 : sc.push_back( "HG23" );
2119 576 : sc.push_back( "1HG1" );
2120 576 : sc.push_back( "2HG1" );
2121 576 : sc.push_back( "3HG1" );
2122 576 : sc.push_back( "1HG2" );
2123 576 : sc.push_back( "2HG2" );
2124 576 : sc.push_back( "3HG2" );
2125 576 : return sc;
2126 : } else {
2127 0 : plumed_merror("Sidechain atoms unknown: " + s);
2128 : }
2129 0 : }
2130 :
2131 49352 : bool CS2Backbone::isSP2(const std::string & resType, const std::string & atomName) {
2132 : bool sp2 = false;
2133 49352 : if (atomName == "C") {
2134 : return true;
2135 : }
2136 46008 : if (atomName == "O") {
2137 : return true;
2138 : }
2139 :
2140 42716 : if(resType == "TRP") {
2141 396 : if (atomName == "CG") {
2142 : sp2 = true;
2143 378 : } else if (atomName == "CD1") {
2144 : sp2 = true;
2145 360 : } else if (atomName == "CD2") {
2146 : sp2 = true;
2147 342 : } else if (atomName == "CE2") {
2148 : sp2 = true;
2149 324 : } else if (atomName == "CE3") {
2150 : sp2 = true;
2151 306 : } else if (atomName == "CZ2") {
2152 : sp2 = true;
2153 288 : } else if (atomName == "CZ3") {
2154 : sp2 = true;
2155 270 : } else if (atomName == "CH2") {
2156 : sp2 = true;
2157 : }
2158 42320 : } else if (resType == "ASP") {
2159 2008 : if (atomName == "CG") {
2160 : sp2 = true;
2161 1814 : } else if (atomName == "OD1") {
2162 : sp2 = true;
2163 1620 : } else if (atomName == "OD2") {
2164 : sp2 = true;
2165 : }
2166 40312 : } else if (resType == "GLU") {
2167 2844 : if (atomName == "CD") {
2168 : sp2 = true;
2169 2628 : } else if (atomName == "OE1") {
2170 : sp2 = true;
2171 2412 : } else if (atomName == "OE2") {
2172 : sp2 = true;
2173 : }
2174 37468 : } else if (resType == "ARG") {
2175 3124 : if (atomName == "CZ") {
2176 : sp2 = true;
2177 : }
2178 34344 : } else if (resType == "HIS") {
2179 0 : if (atomName == "CG") {
2180 : sp2 = true;
2181 0 : } else if (atomName == "ND1") {
2182 : sp2 = true;
2183 0 : } else if (atomName == "CD2") {
2184 : sp2 = true;
2185 0 : } else if (atomName == "CE1") {
2186 : sp2 = true;
2187 0 : } else if (atomName == "NE2") {
2188 : sp2 = true;
2189 : }
2190 34344 : } else if (resType == "PHE") {
2191 5472 : if (atomName == "CG") {
2192 : sp2 = true;
2193 5168 : } else if (atomName == "CD1") {
2194 : sp2 = true;
2195 4864 : } else if (atomName == "CD2") {
2196 : sp2 = true;
2197 4560 : } else if (atomName == "CE1") {
2198 : sp2 = true;
2199 4256 : } else if (atomName == "CE2") {
2200 : sp2 = true;
2201 3952 : } else if (atomName == "CZ") {
2202 : sp2 = true;
2203 : }
2204 28872 : } else if (resType == "TYR") {
2205 684 : if (atomName == "CG") {
2206 : sp2 = true;
2207 648 : } else if (atomName == "CD1") {
2208 : sp2 = true;
2209 612 : } else if (atomName == "CD2") {
2210 : sp2 = true;
2211 576 : } else if (atomName == "CE1") {
2212 : sp2 = true;
2213 540 : } else if (atomName == "CE2") {
2214 : sp2 = true;
2215 504 : } else if (atomName == "CZ") {
2216 : sp2 = true;
2217 : }
2218 28188 : } else if (resType == "ASN") {
2219 1944 : if (atomName == "CG") {
2220 : sp2 = true;
2221 1782 : } else if (atomName == "OD1") {
2222 : sp2 = true;
2223 : }
2224 26244 : } else if (resType == "GLN") {
2225 810 : if (atomName == "CD") {
2226 : sp2 = true;
2227 756 : } else if (atomName == "OE1") {
2228 : sp2 = true;
2229 : }
2230 : }
2231 :
2232 : return sp2;
2233 : }
2234 :
2235 149312 : bool CS2Backbone::is_chi1_cx(const std::string & frg, const std::string & atm) {
2236 149312 : if(atm=="CG") {
2237 : return true;
2238 : }
2239 143212 : if((frg == "CYS")&&(atm =="SG")) {
2240 : return true;
2241 : }
2242 143212 : if(((frg == "ILE")||(frg == "VAL"))&&(atm == "CG1")) {
2243 : return true;
2244 : }
2245 142096 : if((frg == "SER")&&(atm == "OG")) {
2246 : return true;
2247 : }
2248 141402 : if((frg == "THR")&&(atm == "OG1")) {
2249 : return true;
2250 : }
2251 :
2252 : return false;
2253 : }
2254 :
2255 3710008 : void CS2Backbone::xdist_name_map(std::string & name) {
2256 3710008 : if((name == "OT1")||(name == "OC1")) {
2257 : name = "O";
2258 3710008 : } else if ((name == "HN") || (name == "HT1") || (name == "H1")) {
2259 : name = "H";
2260 3684678 : } else if ((name == "CG1")|| (name == "OG")||
2261 7379810 : (name == "SG") || (name == "OG1")) {
2262 : name = "CG";
2263 3655334 : } else if ((name == "HA1")|| (name == "HA3")) {
2264 : name = "HA";
2265 : }
2266 3710008 : }
2267 :
2268 : }
2269 : }
|
