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