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