1 // ENG1_SF.H : calculations for simplified protein models (molecular mechanics). 2 3 // Copyright (C) 1998 Tommi Hassinen. 4 5 // This package is free software; you can redistribute it and/or modify 6 // it under the terms of the GNU General Public License as published by 7 // the Free Software Foundation; either version 2 of the License, or 8 // (at your option) any later version. 9 10 // This package is distributed in the hope that it will be useful, 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU General Public License for more details. 14 15 // You should have received a copy of the GNU General Public License 16 // along with this package; if not, write to the Free Software 17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 18 19 /*################################################################################################*/ 20 21 #ifndef ENG1_SF_H 22 #define ENG1_SF_H 23 24 struct eng1_sf_param; 25 26 class sf_chn; // chain 27 class sf_res; // residue 28 29 class sf_dsb; // disulphide bridge 30 31 class sf_helix4info; 32 class sf_strandinfo; 33 class sf_strandpair; 34 35 class setup1_sf; 36 37 struct sf_bt1_data; // saved distance results. 38 struct sf_bt2_data; // saved angle results. 39 struct sf_nonbonded_lookup; 40 41 struct sf_bt1; // virtual-bond stretching. 42 struct sf_bt2; // virtual-angle bending. 43 struct sf_bt3; // main-chain torsion/dipole terms. 44 struct sf_bt4; // 1st side-chain virtual-atom terms. 45 46 struct sf_nbt1; // nonbonded terms. 47 48 struct sf_nbt3_nd; 49 struct sf_nbt3_ipd; 50 51 struct sf_nbt3_nl; 52 struct sf_nbt3_coi; 53 struct sf_nbt3_ips; 54 struct sf_nbt3_arc; 55 56 #define LAYERS 3 // NOT YET PROPERLY TESTED!!! 2 layers for all atoms a better compromise??? 57 #define LAYER_LIMIT 0.10 // THEREFORE DISABLED IN CODE!!! but still makes the arrays -> memory!!! 58 59 static const int size_nl[LAYERS] = { 100, 200, 400 }; 60 #define MAX_SIZE_NL 400 // size_nl[2] = 400 61 62 #define SIZE_NT 100 63 #define SIZE_COI 100 64 #define SIZE_IPD 50 65 #define SIZE_IPS 100 66 #define SIZE_ARC 100 67 68 #define INDEX_FLAG 0x8000000 // index of the point 69 #define ORDER_FLAG 0x4000000 // 0 = starting point, 1 = ending point 70 #define FLAG_MASK ~(INDEX_FLAG | ORDER_FLAG) 71 72 #define BETA 52.0 // beta-angle of the bt4-term... 73 74 class eng1_sf; 75 76 /*################################################################################################*/ 77 78 #include "engine.h" 79 #include "model.h" 80 81 #include <vector> 82 using namespace std; 83 84 /*################################################################################################*/ 85 86 struct eng1_sf_param 87 { 88 f64 vdwrad; 89 f64 lenjon; 90 f64 wescc; 91 f64 wescd; 92 f64 dipole1; 93 f64 dipole2; 94 f64 epsilon1; 95 f64 epsilon2; 96 f64 epsilon3; 97 f64 epsilon4; 98 f64 epsilon5; 99 f64 epsilon9; 100 f64 exp_solv_1n; 101 f64 exp_solv_1p; 102 f64 exp_solv_2; 103 f64 imp_solv_1n; 104 f64 imp_solv_1p; 105 f64 imp_solv_2; 106 f64 solvrad; 107 f64 wang; 108 f64 wtor1; 109 f64 wtor2; 110 f64 wrep; 111 112 f64 charge_wes[4]; 113 f64 charge_sasa1[4]; 114 f64 charge_sasa2[4]; 115 f64 charge_pKa[9]; 116 bool charge_acid[9]; 117 118 f64 rms; // for random search only... ERROR??? 119 120 f64 pH; ///< This is the ONLY value that can be modified by user!!! 121 122 bool operator<(const eng1_sf_param & p1) const { return (rms < p1.rms); } 123 }; 124 125 void eng1_sf_param_SetDefaultValues(eng1_sf_param &); 126 void eng1_sf_param_MakeCopy(eng1_sf_param &, eng1_sf_param &); 127 128 /*################################################################################################*/ 129 130 class sf_chn // chain 131 { 132 // protected: // strict... 133 public: // loose ; this is for ghemical... 134 135 vector<sf_res> res_vector; 136 137 friend void CopyCRD(model *, engine *, i32u); 138 friend void CopyCRD(engine *, model *, i32u); 139 140 friend class setup1_sf; 141 friend class eng1_sf; 142 143 public: 144 145 sf_chn(void); 146 sf_chn(const sf_chn &); 147 ~sf_chn(void); 148 }; 149 150 /*################################################################################################*/ 151 152 #define SF_STATE_HELIX4 0 153 #define SF_STATE_STRAND 1 154 #define SF_STATE_LOOP 2 155 156 class sf_res // residue 157 { 158 protected: 159 160 char symbol; 161 atom * peptide[5]; // pointers to atoms of the peptide unit. 162 163 i32s natm; 164 atom * atmr[3]; // pointers to "simplified" atoms. 165 i32s loc_varind[3]; // SF-atom indices of "simplified" atoms. 166 167 i32u state; 168 169 friend void CopyCRD(model *, engine *, i32u); 170 friend void CopyCRD(engine *, model *, i32u); 171 172 friend class setup1_sf; 173 friend class eng1_sf; 174 175 public: 176 177 sf_res(void); 178 sf_res(char, atom *, atom *, atom *, atom *, atom *, i32s, atom *, atom *, atom *, i32s, i32s, i32s); 179 sf_res(const sf_res &); 180 ~sf_res(void); 181 GetNumA(void)182 i32s GetNumA(void) { return natm; } GetRefA(i32u atmi)183 atom * GetRefA(i32u atmi) { return (atmi < 3 ? atmr[atmi] : NULL); } 184 GetState(void)185 i32u GetState(void) { return state; } 186 GetSymbol(void)187 char GetSymbol(void) { return symbol; } 188 }; 189 190 /*################################################################################################*/ 191 192 class sf_dsb // disulphide bridge 193 { 194 protected: 195 196 i32s chn[2]; 197 i32s res[2]; 198 199 friend class setup1_sf; 200 friend class eng1_sf; 201 202 public: 203 sf_dsb(void)204 sf_dsb(void) { } ~sf_dsb(void)205 ~sf_dsb(void) { } 206 GetChn(i32s i)207 i32s GetChn(i32s i) 208 { 209 if (i < 0 || i > 1) 210 { 211 cout << "sf_dsb::GetChn() : index overflow!" << endl; 212 exit(EXIT_FAILURE); 213 } 214 215 return chn[i]; 216 } 217 GetRes(i32s i)218 i32s GetRes(i32s i) 219 { 220 if (i < 0 || i > 1) 221 { 222 cout << "sf_dsb::GetRes() : index overflow!" << endl; 223 exit(EXIT_FAILURE); 224 } 225 226 return res[i]; 227 } 228 }; 229 230 /*################################################################################################*/ 231 232 class sf_helix4info 233 { 234 protected: 235 236 i32u chn; 237 i32u res[2]; 238 239 vector<atom *> mc_H_don; 240 vector<atom *> mc_H_acc; 241 242 vector<atom *> ca_H_don; 243 vector<atom *> ca_H_acc; 244 245 friend class setup1_sf; 246 247 friend class project; // this is for ghemical... 248 249 public: 250 sf_helix4info(i32u c,i32u rb,i32u re)251 sf_helix4info(i32u c, i32u rb, i32u re) 252 { 253 chn = c; 254 res[0] = rb; 255 res[1] = re; 256 } 257 ~sf_helix4info(void)258 ~sf_helix4info(void) 259 { 260 } 261 GetChn(void)262 i32u GetChn(void) { return chn; } GetResBgn(void)263 i32u GetResBgn(void) { return res[0]; } GetResEnd(void)264 i32u GetResEnd(void) { return res[1]; } 265 }; 266 267 class sf_strandinfo 268 { 269 protected: 270 271 i32u chn; 272 i32u res[2]; 273 274 public: 275 sf_strandinfo(i32u c,i32u rb,i32u re)276 sf_strandinfo(i32u c, i32u rb, i32u re) 277 { 278 chn = c; 279 res[0] = rb; 280 res[1] = re; 281 282 // cout << "DEBUG sf_strandinfo ; created record "; 283 // cout << "c= " << chn << " rb= " << res[0] << " re= " << res[1]; 284 // cout << endl; 285 } 286 ~sf_strandinfo(void)287 ~sf_strandinfo(void) 288 { 289 } 290 GetChn(void)291 i32u GetChn(void) { return chn; } GetResBgn(void)292 i32u GetResBgn(void) { return res[0]; } GetResEnd(void)293 i32u GetResEnd(void) { return res[1]; } 294 }; 295 296 class sf_strandpair 297 { 298 protected: 299 300 i32u chn[2]; // the first residues in chain direction... 301 i32u res[2]; // the first residues in chain direction... 302 303 i32u length; 304 bool parallel; 305 306 friend class setup1_sf; 307 308 friend class project; // this is for ghemical... 309 310 public: 311 312 vector<atom *> mc_S_don; 313 vector<atom *> mc_S_acc; 314 315 vector<atom *> ca_S_2x; 316 vector<atom *> cx_S_2x; 317 318 public: 319 sf_strandpair(i32u c1,i32u r1,i32u c2,i32u r2,i32u l,bool p)320 sf_strandpair(i32u c1, i32u r1, i32u c2, i32u r2, i32u l, bool p) 321 { 322 chn[0] = c1; res[0] = r1; 323 chn[1] = c2; res[1] = r2; 324 325 length = l; 326 parallel = p; 327 } 328 ~sf_strandpair(void)329 ~sf_strandpair(void) 330 { 331 } 332 ContainsPair(i32u c1,i32u r1,i32u c2,i32u r2)333 bool ContainsPair(i32u c1, i32u r1, i32u c2, i32u r2) 334 { 335 for (i32s n1 = 0;n1 < (i32s) length;n1++) 336 { 337 for (i32s dir = 0;dir < 2;dir++) 338 { 339 if (c1 != chn[dir] || c2 != chn[!dir]) continue; 340 341 i32s tmpr[2]; tmpr[dir] = ((i32s) res[0]) + n1; 342 tmpr[!dir] = ((i32s) res[1]) + (parallel ? +n1 : -n1); 343 344 if (r1 != (i32u) tmpr[0] || r2 != (i32u) tmpr[1]) continue; 345 346 return true; 347 } 348 } 349 350 return false; 351 } 352 IsParallel(void)353 bool IsParallel(void) { return parallel; } IsAntiParallel(void)354 bool IsAntiParallel(void) { return !parallel; } 355 }; 356 357 /*################################################################################################*/ 358 359 /// A setup class for SF submodels. 360 361 class setup1_sf : virtual public setup 362 { 363 // protected: // strict... 364 public: // loose ; this is for ghemical... 365 366 vector<sf_chn> chn_vector; 367 vector<sf_dsb> dsb_vector; 368 369 enum SFmode { modeUA = 0, modeP5 = 1, modeP3 = 2 } mode; 370 371 vector<sf_helix4info> hi_vector; 372 vector<sf_strandinfo> si_vector; 373 vector<sf_strandpair> sp_vector; 374 375 vector<atom *> tb_vis1_vector; 376 vector<float> tb_vis2_vector; 377 378 protected: 379 380 eng1_sf_param prm; 381 382 friend void CopyCRD(model *, engine *, i32u); 383 friend void CopyCRD(engine *, model *, i32u); 384 385 friend class eng1_sf; 386 387 public: 388 389 setup1_sf(model *, SFmode, bool = true); 390 ~setup1_sf(void); 391 392 void UpdateAtomFlags(void); // virtual 393 394 void GetReducedCRD(iter_al *, vector<i32s> &, fGL *, i32u); 395 396 void StorePStatesToModel(eng1_sf *); 397 398 static i32u static_GetEngineCount(void); 399 static i32u static_GetEngineIDNumber(i32u); 400 static const char * static_GetEngineName(i32u); 401 static const char * static_GetClassName(void); 402 403 i32u GetEngineCount(void); // virtual 404 i32u GetEngineIDNumber(i32u); // virtual 405 const char * GetEngineName(i32u); // virtual 406 const char * GetClassName_lg(void); // virtual 407 408 engine * CreateEngineByIndex(i32u); // virtual 409 }; 410 411 /*################################################################################################*/ 412 413 struct sf_bt1_data // saved distance results. 414 { 415 f64 data1; // len???? 416 f64 data2[2][3]; // dlen???? 417 }; 418 419 struct sf_bt2_data // saved angle results. 420 { 421 f64 data1; // ang???? 422 f64 data2[3][3]; // dang???? 423 }; 424 425 struct sf_nonbonded_lookup 426 { 427 char s1; i32u a1; 428 char s2; i32u a2; 429 430 f64 opte; 431 }; 432 433 #define TTYPE_LOOP 0x00 434 #define TTYPE_HELIX 0x01 435 #define TTYPE_STRAND 0x02 436 437 #define TTYPE_SIDE 0x10 438 #define TTYPE_BRIDGE 0x11 439 #define TTYPE_TERM 0x12 440 441 struct sf_bt1 // virtual-bond stretching. 442 { 443 i32s atmi[2]; 444 445 f64 opt; f64 fc; 446 GetIndexsf_bt1447 i32s GetIndex(i32s p1, bool p2) 448 { 449 return atmi[p2 ? p1 : !p1]; 450 } 451 }; 452 453 struct sf_bt2 // virtual-angle bending. 454 { 455 i32s atmi[3]; 456 457 i32s index1[2]; 458 bool dir1[2]; 459 460 i32s ttype; 461 462 f64 opt; f64 fc[2]; 463 }; 464 465 struct sf_bt3 // main-chain torsion/dipole terms. 466 { 467 i32s atmi[4]; 468 469 i32s index2[2]; 470 i32s index1[4]; 471 bool dir1[4]; 472 473 i32s tor_ttype; 474 475 f64 torc[3]; 476 f64 tors[3]; 477 478 // dipole term starts here... dipole term starts here... dipole term starts here... 479 // dipole term starts here... dipole term starts here... dipole term starts here... 480 // dipole term starts here... dipole term starts here... dipole term starts here... 481 482 f64 pbdd; 483 484 f64 bv[2][3]; 485 f64 dbv[2][3][3][3]; 486 487 f64 dv[3]; 488 f64 ddv[4][3][3]; 489 490 i32s dip_ttype; 491 492 bool skip; 493 494 f64 dipc[3]; 495 f64 dips[3]; 496 f64 dipk[2]; 497 }; 498 499 struct sf_bt4 // 1st side-chain virtual-atom terms. 500 { 501 i32s index1; 502 i32s index2; 503 504 f64 opt; f64 fc; 505 506 f64 fscos[3]; 507 f64 fssin[3]; 508 }; 509 510 struct sf_nbt1 // nonbonded terms. 511 { 512 i32s atmi[2]; 513 514 f64 data[2]; 515 516 bool operator==(const sf_nbt1 & p1) const 517 { 518 if ((atmi[0] == p1.atmi[0]) && (atmi[1] == p1.atmi[1])) return true; 519 if ((atmi[0] == p1.atmi[1]) && (atmi[1] == p1.atmi[0])) return true; 520 521 return false; 522 } 523 }; 524 525 struct sf_nbt3_nd ///< SASA neighbor data. 526 { 527 i32s index; 528 f64 distance; 529 530 // these are sorted in reverse order, from large to small... 531 // these are sorted in reverse order, from large to small... 532 // these are sorted in reverse order, from large to small... 533 534 bool operator<(const sf_nbt3_nd & p1) const 535 { 536 return (distance > p1.distance); 537 } 538 }; 539 540 struct sf_nbt3_ipd ///< SASA intersection point data. 541 { 542 f64 angle; 543 i32u ipdata; 544 545 bool operator<(const sf_nbt3_ipd & p1) const 546 { 547 return (angle < p1.angle); 548 } 549 }; 550 551 struct sf_nbt3_nl ///< SASA neighbor list. 552 { 553 i32s index_count; 554 i32s * index; 555 }; 556 557 struct sf_nbt3_coi ///< SASA circle of intersection. 558 { 559 i32s index; 560 561 i32s ipd_count; 562 sf_nbt3_ipd ipdt[SIZE_IPD]; 563 564 f64 refv[3]; 565 566 f64 dist; 567 f64 dv[3]; f64 ddv[3][3]; 568 569 f64 g; f64 dg[3]; 570 f64 ct; f64 dct[3]; 571 572 bool flag; 573 AddIPDsf_nbt3_coi574 void AddIPD(f64 * p1, i32u p2) 575 { 576 ipdt[ipd_count].ipdata = p2; 577 578 if (!ipd_count) 579 { 580 f64 t1a[3]; 581 t1a[0] = dv[0] * p1[0]; 582 t1a[1] = dv[1] * p1[1]; 583 t1a[2] = dv[2] * p1[2]; 584 585 f64 t1b = t1a[0] + t1a[1] + t1a[2]; 586 587 refv[0] = p1[0] - dv[0] * t1b; 588 refv[1] = p1[1] - dv[1] * t1b; 589 refv[2] = p1[2] - dv[2] * t1b; 590 591 f64 t1c = sqrt(refv[0] * refv[0] + refv[1] * refv[1] + refv[2] * refv[2]); 592 refv[0] /= t1c; refv[1] /= t1c; refv[2] /= t1c; 593 594 ipdt[ipd_count].angle = 0.0; 595 } 596 else 597 { 598 f64 t1a[3]; 599 t1a[0] = dv[0] * p1[0]; 600 t1a[1] = dv[1] * p1[1]; 601 t1a[2] = dv[2] * p1[2]; 602 603 f64 t1b = t1a[0] + t1a[1] + t1a[2]; 604 605 f64 t2a[3]; 606 t2a[0] = p1[0] - dv[0] * t1b; 607 t2a[1] = p1[1] - dv[1] * t1b; 608 t2a[2] = p1[2] - dv[2] * t1b; 609 610 f64 t1c = sqrt(t2a[0] * t2a[0] + t2a[1] * t2a[1] + t2a[2] * t2a[2]); 611 t2a[0] /= t1c; t2a[1] /= t1c; t2a[2] /= t1c; 612 613 f64 t1d = refv[0] * t2a[0] + refv[1] * t2a[1] + refv[2] * t2a[2]; 614 if (t1d < -1.0) t1d = -1.0; // domain check... 615 if (t1d > +1.0) t1d = +1.0; // domain check... 616 617 f64 t9a = acos(t1d); 618 619 f64 t3a[3]; 620 t3a[0] = dv[1] * t2a[2] - dv[2] * t2a[1]; 621 t3a[1] = dv[2] * t2a[0] - dv[0] * t2a[2]; 622 t3a[2] = dv[0] * t2a[1] - dv[1] * t2a[0]; 623 624 f64 t9b = refv[0] * t3a[0] + refv[1] * t3a[1] + refv[2] * t3a[2]; 625 626 if (t9b < 0.0) ipdt[ipd_count].angle = -t9a; 627 else ipdt[ipd_count].angle = +t9a; 628 } 629 630 ipd_count++; 631 if (ipd_count >= SIZE_IPD) { cout << "BUG: IPD overflow!!!" << endl; exit(EXIT_FAILURE); } 632 } 633 }; 634 635 struct sf_nbt3_ips ///< SASA intersection points. 636 { 637 i32s coi[2]; 638 639 f64 ipv[2][3]; 640 f64 dipv[2][2][3][3]; 641 }; 642 643 struct sf_nbt3_arc ///< SASA positively oriented arc. 644 { 645 i32s coi; 646 i32s index[2][2]; 647 648 i32u ipdata[2]; 649 650 f64 tv[2][3]; 651 f64 dtv[2][2][3][3]; 652 653 bool flag; 654 }; 655 656 /*################################################################################################*/ 657 658 class eng1_sf : virtual public engine_bp 659 { 660 protected: 661 662 i32u * l2g_sf; ///< The local-to-global lookup table. 663 i32s * index_chn; ///< This maps the (local) atoms to setup1_sf chains/residues. 664 i32s * index_res; ///< This maps the (local) atoms to setup1_sf chains/residues. 665 i32s num_solvent; 666 667 eng1_sf_param * myprm; 668 bool use_explicit_solvent; 669 bool use_implicit_solvent; 670 671 vector<sf_bt1> bt1_vector; 672 vector<sf_bt2> bt2_vector; 673 vector<sf_bt3> bt3_vector; 674 vector<sf_bt4> bt4_vector; 675 676 sf_bt1_data * bt1data; 677 sf_bt2_data * bt2data; 678 679 vector<sf_nbt1> nbt1_vector; 680 681 f64 * mass; // MD??? 682 f64 * vdwr; 683 f64 * charge1; // actual charges. 684 f64 * charge2; // neutralized charges. 685 686 i32s * dist1; f64 * dist2; 687 sf_nbt3_nl * nbt3_nl[LAYERS]; 688 689 f64 * vdwr1[LAYERS]; 690 f64 * vdwr2[LAYERS]; 691 f64 * sasaE[LAYERS]; 692 693 fGL * solv_exp[LAYERS]; 694 695 public: 696 697 f64 energy_bt1; 698 f64 energy_bt2; 699 f64 energy_bt3a; 700 f64 energy_bt3b; 701 f64 energy_bt4a; 702 f64 energy_bt4b; 703 f64 energy_nbt1a; 704 f64 energy_nbt1b; 705 f64 energy_nbt1c; 706 f64 energy_nbt2a; 707 f64 energy_nbt2b; 708 f64 energy_nbt2c; 709 f64 energy_nbt3a; 710 f64 energy_nbt3b; 711 f64 energy_pnlt; 712 713 i32s * tmp_vartab; 714 f64 * tmp_parames; 715 f64 * tmp_paramsa1; 716 f64 * tmp_paramsa2; 717 f64 * tmp_newpKa; 718 719 f64 constraints; 720 721 friend void CopyCRD(model *, engine *, i32u); 722 friend void CopyCRD(engine *, model *, i32u); 723 724 friend class moldyn_langevin; 725 726 public: 727 728 eng1_sf(setup *, i32u, bool, bool); 729 virtual ~eng1_sf(void); 730 731 bool SetTorsionConstraint(atom *, atom *, atom *, atom *, f64, f64, bool); // virtual 732 bool RemoveTorsionConstraint(atom *, atom *, atom *, atom *); // virtual 733 734 void SetupCharges(void); 735 void GetChgGrpVar(i32s, i32s &, i32s &); 736 737 void Compute(i32u, bool = false); // virtual 738 GetOrbitalCount(void)739 i32s GetOrbitalCount(void) { return 0; } // virtual GetOrbitalEnergy(i32s)740 f64 GetOrbitalEnergy(i32s) { return 0.0; } // virtual 741 GetElectronCount(void)742 i32s GetElectronCount(void) { return 0; } // virtual 743 SetupPlotting(void)744 void SetupPlotting(void) { } // virtual 745 746 fGL GetVDWSurf(fGL *, fGL *); // virtual 747 748 fGL GetESP(fGL *, fGL *); // virtual 749 GetElDens(fGL *,fGL *)750 fGL GetElDens(fGL *, fGL *) { return 0.0; } // virtual 751 GetOrbital(fGL *,fGL *)752 fGL GetOrbital(fGL *, fGL *) { return 0.0; } // virtual GetOrbDens(fGL *,fGL *)753 fGL GetOrbDens(fGL *, fGL *) { return 0.0; } // virtual 754 755 protected: 756 757 bool InitNBT1(sf_nbt1 *, vector<sf_nonbonded_lookup> &); 758 void InitLenJon(sf_nbt1 *, f64, f64); 759 760 void ComputeBT1(i32u); 761 void ComputeBT2(i32u); 762 void ComputeBT3(i32u); 763 void ComputeBT4(i32u); 764 765 void ComputeNBT1(i32u); 766 void ComputeNBT2(i32u); 767 void ComputeNBT3(i32u); 768 }; 769 770 /*################################################################################################*/ 771 772 #endif // ENG1_SF_H 773 774 // eof 775