1 /********************************************************************* 2 forcefieldmmff94.cpp - MMFF94 force field 3 4 Copyright (C) 2006-2008 by Tim Vandermeersch <tim.vandermeersch@gmail.com> 5 6 This file is part of the Open Babel project. 7 For more information, see <http://openbabel.org/> 8 9 This program is free software; you can redistribute it and/or modify 10 it under the terms of the GNU General Public License as published by 11 the Free Software Foundation version 2 of the License. 12 13 This program is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 GNU General Public License for more details. 17 ***********************************************************************/ 18 19 /* 20 * Source code layout: 21 * - Functions to calculate the actual interactions 22 * - Parse parameter files 23 * - Setup Functions 24 * - Validation functions 25 * - Calculate bond type, angle type, stretch-bend type, torsion type 26 * - Various tables & misc. functions 27 * 28 */ 29 30 #include <openbabel/babelconfig.h> 31 #include <openbabel/obconversion.h> 32 #include <openbabel/mol.h> 33 #include <openbabel/locale.h> 34 #include <openbabel/elements.h> 35 #include <openbabel/bond.h> 36 #include <openbabel/obiter.h> 37 #include <openbabel/ring.h> 38 39 #ifndef M_PI 40 #define M_PI 3.14159265358979323846 41 #endif 42 43 #include <iomanip> 44 #include "forcefieldmmff94.h" 45 46 using namespace std; 47 48 namespace OpenBabel 49 { 50 //////////////////////////////////////////////////////////////////////////////// 51 //////////////////////////////////////////////////////////////////////////////// 52 // 53 // Functions to calculate the actual interactions 54 // 55 //////////////////////////////////////////////////////////////////////////////// 56 //////////////////////////////////////////////////////////////////////////////// 57 Energy(bool gradients)58 double OBForceFieldMMFF94::Energy(bool gradients) 59 { 60 double energy; 61 62 IF_OBFF_LOGLVL_MEDIUM 63 OBFFLog("\nE N E R G Y\n\n"); 64 65 if (gradients) { 66 ClearGradients(); 67 energy = E_Bond<true>(); 68 energy += E_Angle<true>(); 69 energy += E_StrBnd<true>(); 70 energy += E_Torsion<true>(); 71 energy += E_OOP<true>(); 72 energy += E_VDW<true>(); 73 energy += E_Electrostatic<true>(); 74 } else { 75 energy = E_Bond<false>(); 76 energy += E_Angle<false>(); 77 energy += E_StrBnd<false>(); 78 energy += E_Torsion<false>(); 79 energy += E_OOP<false>(); 80 energy += E_VDW<false>(); 81 energy += E_Electrostatic<false>(); 82 } 83 84 IF_OBFF_LOGLVL_MEDIUM { 85 snprintf(_logbuf, BUFF_SIZE, "\nTOTAL ENERGY = %8.5f %s\n", energy, GetUnit().c_str()); 86 OBFFLog(_logbuf); 87 } 88 89 return energy; 90 } 91 92 // 93 // MMFF part I - page 494 94 // 95 // kb_ij 7 96 // EB_ij = 143.9325 ------- /\r_ij^2 (1 + cs /\_rij + ---- cs^2 r_ij^2) 97 // 2 12 98 // 99 // kb_ij force constant (md/A) 100 // 101 // /\r_ij r_ij - r0_ij (A) 102 // 103 // cs cubic stretch constant = -2 A^(-1) 104 // 105 template<bool gradients> Compute()106 inline void OBFFBondCalculationMMFF94::Compute() 107 { 108 if (OBForceField::IgnoreCalculation(idx_a, idx_b)) { 109 energy = 0.0; 110 return; 111 } 112 113 double delta2; 114 115 if (gradients) { 116 rab = OBForceField::VectorBondDerivative(pos_a, pos_b, force_a, force_b); 117 delta = rab - r0; 118 delta2 = delta * delta; 119 120 const double dE = 143.9325 * kb * delta * (1.0 - 3.0 * delta + 14.0/3.0 * delta2); 121 122 OBForceField::VectorSelfMultiply(force_a, dE); 123 OBForceField::VectorSelfMultiply(force_b, dE); 124 } else { 125 rab = OBForceField::VectorDistance(pos_a, pos_b); 126 delta = rab - r0; 127 delta2 = delta * delta; 128 } 129 130 energy = kb * delta2 * (1.0 - 2.0 * delta + 7.0/3.0 * delta2); 131 } 132 133 template<bool gradients> E_Bond()134 double OBForceFieldMMFF94::E_Bond() 135 { 136 double energy = 0.0; 137 138 IF_OBFF_LOGLVL_HIGH { 139 OBFFLog("\nB O N D S T R E T C H I N G\n\n"); 140 OBFFLog("ATOM TYPES FF BOND IDEAL FORCE\n"); 141 OBFFLog(" I J CLASS LENGTH LENGTH CONSTANT DELTA ENERGY\n"); 142 OBFFLog("------------------------------------------------------------------------\n"); 143 } 144 145 #ifdef _OPENMP 146 #pragma omp parallel for reduction(+:energy) 147 #endif 148 for (int i = 0; i < _bondcalculations.size(); ++i) { 149 _bondcalculations[i].template Compute<gradients>(); 150 energy += _bondcalculations[i].energy; 151 152 #ifndef _OPENMP 153 if (gradients) { 154 AddGradient(_bondcalculations[i].force_a, _bondcalculations[i].idx_a); 155 AddGradient(_bondcalculations[i].force_b, _bondcalculations[i].idx_b); 156 } 157 #endif 158 159 IF_OBFF_LOGLVL_HIGH { 160 snprintf(_logbuf, BUFF_SIZE, "%2d %2d %d %8.3f %8.3f %8.3f %8.3f %8.3f\n", 161 atoi(_bondcalculations[i].a->GetType()), atoi(_bondcalculations[i].b->GetType()), 162 _bondcalculations[i].bt, _bondcalculations[i].rab, _bondcalculations[i].r0, 163 _bondcalculations[i].kb, _bondcalculations[i].delta, 164 143.9325 * 0.5 * _bondcalculations[i].energy); 165 OBFFLog(_logbuf); 166 } 167 } 168 169 #ifdef _OPENMP 170 for (int i = 0; i < _bondcalculations.size(); ++i) { 171 if (gradients) { 172 AddGradient(_bondcalculations[i].force_a, _bondcalculations[i].idx_a); 173 AddGradient(_bondcalculations[i].force_b, _bondcalculations[i].idx_b); 174 } 175 } 176 #endif 177 178 IF_OBFF_LOGLVL_MEDIUM { 179 snprintf(_logbuf, BUFF_SIZE, " TOTAL BOND STRETCHING ENERGY = %8.5f %s\n", 143.9325 * 0.5 * energy, GetUnit().c_str()); 180 OBFFLog(_logbuf); 181 } 182 183 return (143.9325 * 0.5 * energy); 184 } 185 186 // 187 // MMFF part I - page 495 188 // 189 // ka_ijk 190 // EA_ijk = 0.438449325 -------- /\0_ijk^2 (1 + cs /\0_ijk) 191 // 2 192 // 193 // ka_ijk force constant (md A/rad^2) 194 // 195 // /\0_ijk 0_ijk - 00_ijk (degrees) 196 // 197 // cs cubic bend constant = -0.007 deg^-1 = -0.4 rad^-1 198 // 199 template<bool gradients> Compute()200 inline void OBFFAngleCalculationMMFF94::Compute() 201 { 202 if (OBForceField::IgnoreCalculation(idx_a, idx_b, idx_c)) { 203 energy = 0.0; 204 return; 205 } 206 207 double delta2, dE; 208 209 if (gradients) { 210 theta = OBForceField::VectorAngleDerivative(pos_a, pos_b, pos_c, force_a, force_b, force_c); 211 212 if (!isfinite(theta)) 213 theta = 0.0; // doesn't explain why GetAngle is returning NaN but solves it for us; 214 215 delta = theta - theta0; 216 217 if (linear) { 218 energy = 143.9325 * ka * (1.0 + cos((theta) * DEG_TO_RAD)); 219 dE = -sin((theta) * DEG_TO_RAD) * 143.9325 * ka; 220 } else { 221 delta2 = delta * delta; 222 energy = 0.043844 * 0.5 * ka * delta2 * (1.0 - 0.007 * delta); 223 dE = RAD_TO_DEG * 0.043844 * ka * delta * (1.0 - 1.5 * 0.007 * delta); 224 } 225 226 OBForceField::VectorSelfMultiply(force_a, dE); 227 OBForceField::VectorSelfMultiply(force_b, dE); 228 OBForceField::VectorSelfMultiply(force_c, dE); 229 } else { 230 theta = OBForceField::VectorAngle(pos_a, pos_b, pos_c); 231 232 if (!isfinite(theta)) 233 theta = 0.0; // doesn't explain why GetAngle is returning NaN but solves it for us; 234 235 delta = theta - theta0; 236 237 if (linear) { 238 energy = 143.9325 * ka * (1.0 + cos(theta * DEG_TO_RAD)); 239 } else { 240 delta2 = delta * delta; 241 energy = 0.043844 * 0.5 * ka * delta2 * (1.0 - 0.007 * delta); 242 } 243 244 } 245 246 } 247 248 template<bool gradients> E_Angle()249 double OBForceFieldMMFF94::E_Angle() 250 { 251 double energy = 0.0; 252 253 IF_OBFF_LOGLVL_HIGH { 254 OBFFLog("\nA N G L E B E N D I N G\n\n"); 255 OBFFLog("ATOM TYPES FF VALENCE IDEAL FORCE\n"); 256 OBFFLog(" I J K CLASS ANGLE ANGLE CONSTANT DELTA ENERGY\n"); 257 OBFFLog("-----------------------------------------------------------------------------\n"); 258 } 259 260 #ifdef _OPENMP 261 #pragma omp parallel for reduction(+:energy) 262 #endif 263 for (int i = 0; i < _anglecalculations.size(); ++i) { 264 265 _anglecalculations[i].template Compute<gradients>(); 266 energy += _anglecalculations[i].energy; 267 268 #ifndef _OPENMP 269 if (gradients) { 270 AddGradient(_anglecalculations[i].force_a, _anglecalculations[i].idx_a); 271 AddGradient(_anglecalculations[i].force_b, _anglecalculations[i].idx_b); 272 AddGradient(_anglecalculations[i].force_c, _anglecalculations[i].idx_c); 273 } 274 #endif 275 276 IF_OBFF_LOGLVL_HIGH { 277 snprintf(_logbuf, BUFF_SIZE, "%2d %2d %2d %d %8.3f %8.3f %8.3f %8.3f %8.3f\n", 278 atoi(_anglecalculations[i].a->GetType()), atoi(_anglecalculations[i].b->GetType()), 279 atoi(_anglecalculations[i].c->GetType()), _anglecalculations[i].at, 280 _anglecalculations[i].theta, _anglecalculations[i].theta0, 281 _anglecalculations[i].ka, _anglecalculations[i].delta, 282 _anglecalculations[i].energy); 283 OBFFLog(_logbuf); 284 } 285 } 286 287 #ifdef _OPENMP 288 for (int i = 0; i < _anglecalculations.size(); ++i) { 289 if (gradients) { 290 AddGradient(_anglecalculations[i].force_a, _anglecalculations[i].idx_a); 291 AddGradient(_anglecalculations[i].force_b, _anglecalculations[i].idx_b); 292 AddGradient(_anglecalculations[i].force_c, _anglecalculations[i].idx_c); 293 } 294 } 295 #endif 296 297 IF_OBFF_LOGLVL_MEDIUM { 298 snprintf(_logbuf, BUFF_SIZE, " TOTAL ANGLE BENDING ENERGY = %8.5f %s\n", energy, GetUnit().c_str()); 299 OBFFLog(_logbuf); 300 } 301 302 return energy; 303 } 304 305 // 306 // MMFF part I - page 495 307 // 308 // EBA_ijk = 2.51210 (kba_ijk /\r_ij + kba_kji /\r_kj) /\0_ijk 309 // 310 // kba_ijk force constant (md/rad) 311 // kba_kji force constant (md/rad) 312 // 313 // /\r_xx see above 314 // /\0_ijk see above 315 // 316 template<bool gradients> Compute()317 inline void OBFFStrBndCalculationMMFF94::Compute() 318 { 319 if (OBForceField::IgnoreCalculation(idx_a, idx_b, idx_c)) { 320 energy = 0.0; 321 return; 322 } 323 324 if (gradients) { 325 theta = OBForceField::VectorAngleDerivative(pos_a, pos_b, pos_c, 326 force_abc_a, force_abc_b, force_abc_c); 327 rab = OBForceField::VectorDistanceDerivative(pos_a, pos_b, force_ab_a, force_ab_b); 328 rbc = OBForceField::VectorDistanceDerivative(pos_b, pos_c, force_bc_b, force_bc_c); 329 } else { 330 theta = OBForceField::VectorAngle(pos_a, pos_b, pos_c); 331 rab = OBForceField::VectorDistance(pos_a, pos_b); 332 rbc = OBForceField::VectorDistance(pos_b, pos_c); 333 } 334 335 if (!isfinite(theta)) 336 theta = 0.0; // doesn't explain why GetAngle is returning NaN but solves it for us; 337 338 delta_theta = theta - theta0; 339 delta_rab = rab - rab0; 340 delta_rbc = rbc - rbc0; 341 const double factor = RAD_TO_DEG * (kbaABC * delta_rab + kbaCBA * delta_rbc); 342 343 energy = DEG_TO_RAD * factor * delta_theta; 344 if (gradients) { 345 //grada = 2.51210 * (kbaABC * rab_da * delta_theta + RAD_TO_DEG * theta_da * (kbaABC * delta_rab + kbaCBA * delta_rbc)); 346 OBForceField::VectorSelfMultiply(force_ab_a, (kbaABC*delta_theta)); 347 OBForceField::VectorSelfMultiply(force_abc_a, factor); 348 OBForceField::VectorAdd(force_ab_a, force_abc_a, force_ab_a); 349 OBForceField::VectorMultiply(force_ab_a, 2.51210, force_a); 350 //gradc = 2.51210 * (kbaCBA * rbc_dc * delta_theta + RAD_TO_DEG * theta_dc * (kbaABC * delta_rab + kbaCBA * delta_rbc)); 351 OBForceField::VectorSelfMultiply(force_bc_c, (kbaCBA*delta_theta)); 352 OBForceField::VectorSelfMultiply(force_abc_c, factor); 353 OBForceField::VectorAdd(force_bc_c, force_abc_c, force_bc_c); 354 OBForceField::VectorMultiply(force_bc_c, 2.51210, force_c); 355 //gradb = -grada - gradc; 356 OBForceField::VectorAdd(force_a, force_c, force_b); 357 OBForceField::VectorSelfMultiply(force_b, -1.0); 358 } 359 } 360 361 template<bool gradients> E_StrBnd()362 double OBForceFieldMMFF94::E_StrBnd() 363 { 364 double energy = 0.0; 365 366 IF_OBFF_LOGLVL_HIGH { 367 OBFFLog("\nS T R E T C H B E N D I N G\n\n"); 368 OBFFLog("ATOM TYPES FF VALENCE DELTA FORCE CONSTANT\n"); 369 OBFFLog(" I J K CLASS ANGLE ANGLE I J J K ENERGY\n"); 370 OBFFLog("---------------------------------------------------------------------------\n"); 371 } 372 373 #ifdef _OPENMP 374 #pragma omp parallel for reduction(+:energy) 375 #endif 376 for (int i = 0; i < _strbndcalculations.size(); ++i) { 377 378 _strbndcalculations[i].template Compute<gradients>(); 379 energy += _strbndcalculations[i].energy; 380 381 #ifndef _OPENMP 382 if (gradients) { 383 AddGradient(_strbndcalculations[i].force_a, _strbndcalculations[i].idx_a); 384 AddGradient(_strbndcalculations[i].force_b, _strbndcalculations[i].idx_b); 385 AddGradient(_strbndcalculations[i].force_c, _strbndcalculations[i].idx_c); 386 } 387 #endif 388 389 IF_OBFF_LOGLVL_HIGH { 390 snprintf(_logbuf, BUFF_SIZE, "%2d %2d %2d %2d %8.3f %8.3f %8.3f %8.3f %8.3f\n", 391 atoi(_strbndcalculations[i].a->GetType()), atoi(_strbndcalculations[i].b->GetType()), 392 atoi(_strbndcalculations[i].c->GetType()), _strbndcalculations[i].sbt, 393 _strbndcalculations[i].theta, _strbndcalculations[i].delta_theta, 394 _strbndcalculations[i].kbaABC, _strbndcalculations[i].kbaCBA, 395 2.51210 * _strbndcalculations[i].energy); 396 OBFFLog(_logbuf); 397 } 398 } 399 400 #ifdef _OPENMP 401 for (int i = 0; i < _strbndcalculations.size(); ++i) { 402 if (gradients) { 403 AddGradient(_strbndcalculations[i].force_a, _strbndcalculations[i].idx_a); 404 AddGradient(_strbndcalculations[i].force_b, _strbndcalculations[i].idx_b); 405 AddGradient(_strbndcalculations[i].force_c, _strbndcalculations[i].idx_c); 406 } 407 } 408 #endif 409 410 IF_OBFF_LOGLVL_MEDIUM { 411 snprintf(_logbuf, BUFF_SIZE, " TOTAL STRETCH BENDING ENERGY = %8.5f %s\n", 2.51210 * energy, GetUnit().c_str()); 412 OBFFLog(_logbuf); 413 } 414 415 return (2.51210 * energy); 416 } 417 GetElementRow(OBAtom * atom)418 int OBForceFieldMMFF94::GetElementRow(OBAtom *atom) 419 { 420 int row; 421 422 row = 0; 423 424 if (atom->GetAtomicNum() > 2) 425 row++; 426 if (atom->GetAtomicNum() > 10) 427 row++; 428 if (atom->GetAtomicNum() > 18) 429 row++; 430 if (atom->GetAtomicNum() > 36) 431 row++; 432 if (atom->GetAtomicNum() > 54) 433 row++; 434 if (atom->GetAtomicNum() > 86) 435 row++; 436 437 return row; 438 } 439 440 // 441 // MMFF part I - page 495 442 // 443 // ET_ijkl = 0.5 ( V1 (1 + cos(0_ijkl)) + V2 (1 - cos(2 0_ijkl)) + V3 (1 + cos(3 0_ijkl)) ) 444 // 445 // V1 force constant (md/rad) 446 // V2 force constant (md/rad) 447 // V3 force constant (md/rad) 448 // 449 // 0_ijkl torsion angle (degrees) 450 // 451 template<bool gradients> Compute()452 inline void OBFFTorsionCalculationMMFF94::Compute() 453 { 454 if (OBForceField::IgnoreCalculation(idx_a, idx_b, idx_c, idx_d)) { 455 energy = 0.0; 456 return; 457 } 458 459 double cosine, cosine2, cosine3; 460 double phi1, phi2, phi3; 461 double dE, sine, sine2, sine3; 462 463 if (gradients) { 464 tor = OBForceField::VectorTorsionDerivative(pos_a, pos_b, pos_c, pos_d, 465 force_a, force_b, force_c, force_d); 466 if (!isfinite(tor)) 467 tor = 1.0e-3; 468 469 sine = sin(DEG_TO_RAD * tor); 470 sine2 = sin(2.0 * DEG_TO_RAD * tor); 471 sine3 = sin(3.0 * DEG_TO_RAD * tor); 472 473 dE = 0.5 * (v1 * sine - 2.0 * v2 * sine2 + 3.0 * v3 * sine3); // MMFF 474 475 OBForceField::VectorSelfMultiply(force_a, dE); 476 OBForceField::VectorSelfMultiply(force_b, dE); 477 OBForceField::VectorSelfMultiply(force_c, dE); 478 OBForceField::VectorSelfMultiply(force_d, dE); 479 } else { 480 tor = OBForceField::VectorTorsion(pos_a, pos_b, pos_c, pos_d); 481 if (!isfinite(tor)) 482 tor = 1.0e-3; 483 } 484 485 cosine = cos(DEG_TO_RAD * tor); 486 cosine2 = cos(DEG_TO_RAD * 2 * tor); 487 cosine3 = cos(DEG_TO_RAD * 3 * tor); 488 489 phi1 = 1.0 + cosine; 490 phi2 = 1.0 - cosine2; 491 phi3 = 1.0 + cosine3; 492 493 energy = (v1 * phi1 + v2 * phi2 + v3 * phi3); 494 495 } 496 497 template<bool gradients> E_Torsion()498 double OBForceFieldMMFF94::E_Torsion() 499 { 500 double energy = 0.0; 501 502 IF_OBFF_LOGLVL_HIGH { 503 OBFFLog("\nT O R S I O N A L\n\n"); 504 OBFFLog("ATOM TYPES FF TORSION FORCE CONSTANT\n"); 505 OBFFLog(" I J K L CLASS ANGLE V1 V2 V3 ENERGY\n"); 506 OBFFLog("--------------------------------------------------------------------\n"); 507 } 508 509 #ifdef _OPENMP 510 #pragma omp parallel for reduction(+:energy) 511 #endif 512 for (int i = 0; i < _torsioncalculations.size(); ++i) { 513 514 _torsioncalculations[i].template Compute<gradients>(); 515 energy += _torsioncalculations[i].energy; 516 517 #ifndef _OPENMP 518 if (gradients) { 519 AddGradient(_torsioncalculations[i].force_a, _torsioncalculations[i].idx_a); 520 AddGradient(_torsioncalculations[i].force_b, _torsioncalculations[i].idx_b); 521 AddGradient(_torsioncalculations[i].force_c, _torsioncalculations[i].idx_c); 522 AddGradient(_torsioncalculations[i].force_d, _torsioncalculations[i].idx_d); 523 } 524 #endif 525 526 IF_OBFF_LOGLVL_HIGH { 527 snprintf(_logbuf, BUFF_SIZE, "%2d %2d %2d %2d %d %8.3f %6.3f %6.3f %6.3f %8.3f\n", 528 atoi(_torsioncalculations[i].a->GetType()), atoi(_torsioncalculations[i].b->GetType()), 529 atoi(_torsioncalculations[i].c->GetType()), atoi(_torsioncalculations[i].d->GetType()), 530 _torsioncalculations[i].tt, _torsioncalculations[i].tor, _torsioncalculations[i].v1, 531 _torsioncalculations[i].v2, _torsioncalculations[i].v3, 0.5 * _torsioncalculations[i].energy); 532 OBFFLog(_logbuf); 533 } 534 535 } 536 537 #ifdef _OPENMP 538 for (int i = 0; i < _torsioncalculations.size(); ++i) { 539 if (gradients) { 540 AddGradient(_torsioncalculations[i].force_a, _torsioncalculations[i].idx_a); 541 AddGradient(_torsioncalculations[i].force_b, _torsioncalculations[i].idx_b); 542 AddGradient(_torsioncalculations[i].force_c, _torsioncalculations[i].idx_c); 543 AddGradient(_torsioncalculations[i].force_d, _torsioncalculations[i].idx_d); 544 } 545 } 546 #endif 547 548 IF_OBFF_LOGLVL_MEDIUM { 549 snprintf(_logbuf, BUFF_SIZE, " TOTAL TORSIONAL ENERGY = %8.5f %s\n", 0.5 * energy, GetUnit().c_str()); 550 OBFFLog(_logbuf); 551 } 552 553 return (0.5 * energy); 554 } 555 556 // // 557 // a // 558 // \ // 559 // b---d plane = a-b-c // 560 // / // 561 // c // 562 // // 563 template<bool gradients> Compute()564 void OBFFOOPCalculationMMFF94::Compute() 565 { 566 if (OBForceField::IgnoreCalculation(idx_a, idx_b, idx_c, idx_d)) { 567 energy = 0.0; 568 return; 569 } 570 571 double angle2, dE; 572 573 if (gradients) { 574 angle = OBForceField::VectorOOPDerivative(pos_a, pos_b, pos_c, pos_d, 575 force_a, force_b, force_c, force_d); 576 577 dE = (-1.0 * RAD_TO_DEG * 0.043844 * angle * koop) / cos(angle * DEG_TO_RAD); 578 579 OBForceField::VectorSelfMultiply(force_a, dE); 580 OBForceField::VectorSelfMultiply(force_b, dE); 581 OBForceField::VectorSelfMultiply(force_c, dE); 582 OBForceField::VectorSelfMultiply(force_d, dE); 583 } else { 584 angle = OBForceField::VectorOOP(pos_a, pos_b, pos_c, pos_d); 585 } 586 587 if (!isfinite(angle)) 588 angle = 0.0; // doesn't explain why GetAngle is returning NaN but solves it for us; 589 590 angle2 = angle * angle; 591 energy = koop * angle2; 592 593 } 594 595 template<bool gradients> E_OOP()596 double OBForceFieldMMFF94::E_OOP() 597 { 598 double energy = 0.0; 599 600 IF_OBFF_LOGLVL_HIGH { 601 OBFFLog("\nO U T - O F - P L A N E B E N D I N G\n\n"); 602 OBFFLog("ATOM TYPES FF OOP FORCE\n"); 603 OBFFLog(" I J K L CLASS ANGLE CONSTANT ENERGY\n"); 604 OBFFLog("----------------------------------------------------------\n"); 605 } 606 607 #ifdef _OPENMP 608 #pragma omp parallel for reduction(+:energy) 609 #endif 610 for (int i = 0; i < _oopcalculations.size(); ++i) { 611 612 _oopcalculations[i].template Compute<gradients>(); 613 energy += _oopcalculations[i].energy; 614 615 #ifndef _OPENMP 616 if (gradients) { 617 AddGradient(_oopcalculations[i].force_a, _oopcalculations[i].idx_a); 618 AddGradient(_oopcalculations[i].force_b, _oopcalculations[i].idx_b); 619 AddGradient(_oopcalculations[i].force_c, _oopcalculations[i].idx_c); 620 AddGradient(_oopcalculations[i].force_d, _oopcalculations[i].idx_d); 621 } 622 #endif 623 624 IF_OBFF_LOGLVL_HIGH { 625 snprintf(_logbuf, BUFF_SIZE, "%2d %2d %2d %2d 0 %8.3f %8.3f %8.3f\n", 626 atoi(_oopcalculations[i].a->GetType()), atoi(_oopcalculations[i].b->GetType()), 627 atoi(_oopcalculations[i].c->GetType()), atoi(_oopcalculations[i].d->GetType()), 628 _oopcalculations[i].angle, _oopcalculations[i].koop, 629 0.043844 * 0.5 * _oopcalculations[i].energy); 630 OBFFLog(_logbuf); 631 } 632 } 633 634 #ifdef _OPENMP 635 for (int i = 0; i < _oopcalculations.size(); ++i) { 636 if (gradients) { 637 AddGradient(_oopcalculations[i].force_a, _oopcalculations[i].idx_a); 638 AddGradient(_oopcalculations[i].force_b, _oopcalculations[i].idx_b); 639 AddGradient(_oopcalculations[i].force_c, _oopcalculations[i].idx_c); 640 AddGradient(_oopcalculations[i].force_d, _oopcalculations[i].idx_d); 641 } 642 } 643 #endif 644 645 IF_OBFF_LOGLVL_MEDIUM { 646 snprintf(_logbuf, BUFF_SIZE, " TOTAL OUT-OF-PLANE BENDING ENERGY = %8.5f %s\n", 0.043844 * 0.5 * energy, GetUnit().c_str()); 647 OBFFLog(_logbuf); 648 } 649 650 return (0.043844 * 0.5 * energy); 651 } 652 653 template<bool gradients> Compute()654 inline void OBFFVDWCalculationMMFF94::Compute() 655 { 656 if (OBForceField::IgnoreCalculation(idx_a, idx_b)) { 657 energy = 0.0; 658 return; 659 } 660 661 if (gradients) { 662 rab = OBForceField::VectorDistanceDerivative(pos_a, pos_b, force_a, force_b); 663 } else { 664 rab = OBForceField::VectorDistance(pos_a, pos_b); 665 } 666 667 const double rab7 = rab*rab*rab*rab*rab*rab*rab; 668 669 double erep = (1.07 * R_AB) / (rab + 0.07 * R_AB); //*** 670 double erep7 = erep*erep*erep*erep*erep*erep*erep; 671 672 double eattr = (((1.12 * R_AB7) / (rab7 + 0.12 * R_AB7)) - 2.0); 673 674 energy = epsilon * erep7 * eattr; 675 676 if (gradients) { 677 const double q = rab / R_AB; 678 const double q6 = q*q*q*q*q*q; 679 const double q7 = q6 * q; 680 erep = 1.07 / (q + 0.07); 681 erep7 = erep*erep*erep*erep*erep*erep*erep; 682 const double term = q7 + 0.12; 683 const double term2 = term * term; 684 eattr = (-7.84 * q6) / term2 + ((-7.84 / term) + 14) / (q + 0.07); 685 const double dE = (epsilon / R_AB) * erep7 * eattr; 686 OBForceField::VectorSelfMultiply(force_a, dE); 687 OBForceField::VectorSelfMultiply(force_b, dE); 688 } 689 } 690 691 template<bool gradients> E_VDW()692 double OBForceFieldMMFF94::E_VDW() 693 { 694 double energy = 0.0; 695 696 IF_OBFF_LOGLVL_HIGH { 697 OBFFLog("\nV A N D E R W A A L S\n\n"); 698 OBFFLog("ATOM TYPES\n"); 699 OBFFLog(" I J Rij R*IJ EPSILON ENERGY\n"); 700 OBFFLog("--------------------------------------------------\n"); 701 // XX XX -000.000 -000.000 -000.000 -000.000 702 } 703 704 #ifdef _OPENMP 705 #pragma omp parallel for reduction(+:energy) 706 #endif 707 for (int i = 0; i < _vdwcalculations.size(); ++i) { 708 // Cut-off check 709 if (_cutoff) 710 if (!_vdwpairs.BitIsSet(_vdwcalculations[i].pairIndex)) 711 continue; 712 713 _vdwcalculations[i].template Compute<gradients>(); 714 energy += _vdwcalculations[i].energy; 715 716 #ifndef _OPENMP 717 if (gradients) { 718 AddGradient(_vdwcalculations[i].force_a, _vdwcalculations[i].idx_a); 719 AddGradient(_vdwcalculations[i].force_b, _vdwcalculations[i].idx_b); 720 } 721 #endif 722 723 IF_OBFF_LOGLVL_HIGH { 724 snprintf(_logbuf, BUFF_SIZE, "%2d %2d %8.3f %8.3f %8.3f %8.3f\n", 725 atoi(_vdwcalculations[i].a->GetType()), atoi(_vdwcalculations[i].b->GetType()), 726 _vdwcalculations[i].rab, _vdwcalculations[i].R_AB, _vdwcalculations[i].epsilon, _vdwcalculations[i].energy); 727 OBFFLog(_logbuf); 728 } 729 730 } 731 732 #ifdef _OPENMP 733 for (int i = 0; i < _vdwcalculations.size(); ++i) { 734 // Cut-off check 735 if (_cutoff) 736 if (!_vdwpairs.BitIsSet(i)) 737 continue; 738 739 if (gradients) { 740 AddGradient(_vdwcalculations[i].force_a, _vdwcalculations[i].idx_a); 741 AddGradient(_vdwcalculations[i].force_b, _vdwcalculations[i].idx_b); 742 } 743 } 744 #endif 745 746 IF_OBFF_LOGLVL_MEDIUM { 747 snprintf(_logbuf, BUFF_SIZE, " TOTAL VAN DER WAALS ENERGY = %8.5f %s\n", energy, GetUnit().c_str()); 748 OBFFLog(_logbuf); 749 } 750 751 return energy; 752 } 753 754 template<bool gradients> Compute()755 inline void OBFFElectrostaticCalculationMMFF94::Compute() 756 { 757 if (OBForceField::IgnoreCalculation(idx_a, idx_b)) { 758 energy = 0.0; 759 return; 760 } 761 762 if (gradients) { 763 rab = OBForceField::VectorDistanceDerivative(pos_a, pos_b, force_a, force_b); 764 rab += 0.05; // ?? 765 const double rab2 = rab * rab; 766 const double dE = -qq / rab2; 767 OBForceField::VectorSelfMultiply(force_a, dE); 768 OBForceField::VectorSelfMultiply(force_b, dE); 769 } else { 770 rab = OBForceField::VectorDistance(pos_a, pos_b); 771 rab += 0.05; // ?? 772 } 773 774 energy = qq / rab; 775 } 776 777 template<bool gradients> E_Electrostatic()778 double OBForceFieldMMFF94::E_Electrostatic() 779 { 780 double energy = 0.0; 781 782 IF_OBFF_LOGLVL_HIGH { 783 OBFFLog("\nE L E C T R O S T A T I C I N T E R A C T I O N S\n\n"); 784 OBFFLog("ATOM TYPES\n"); 785 OBFFLog(" I J Rij Qi Qj ENERGY\n"); 786 OBFFLog("-----------------------------------------------------\n"); 787 // XX XX XXXXXXXX XXXXXXXX XXXXXXXX XXXXXXXX 788 } 789 790 #ifdef _OPENMP 791 #pragma omp parallel for reduction(+:energy) 792 #endif 793 for (int i = 0; i < _electrostaticcalculations.size(); ++i) { 794 // Cut-off check 795 if (_cutoff) 796 if (!_elepairs.BitIsSet(_electrostaticcalculations[i].pairIndex)) 797 continue; 798 799 _electrostaticcalculations[i].template Compute<gradients>(); 800 energy += _electrostaticcalculations[i].energy; 801 802 #ifndef _OPENMP 803 if (gradients) { 804 AddGradient(_electrostaticcalculations[i].force_a, _electrostaticcalculations[i].idx_a); 805 AddGradient(_electrostaticcalculations[i].force_b, _electrostaticcalculations[i].idx_b); 806 } 807 #endif 808 809 IF_OBFF_LOGLVL_HIGH { 810 snprintf(_logbuf, BUFF_SIZE, "%2d %2d %8.3f %8.3f %8.3f %8.3f\n", 811 atoi(_electrostaticcalculations[i].a->GetType()), atoi(_electrostaticcalculations[i].b->GetType()), 812 _electrostaticcalculations[i].rab, _electrostaticcalculations[i].a->GetPartialCharge(), 813 _electrostaticcalculations[i].b->GetPartialCharge(), _electrostaticcalculations[i].energy); 814 OBFFLog(_logbuf); 815 } 816 } 817 818 #ifdef _OPENMP 819 for (int i = 0; i < _electrostaticcalculations.size(); ++i) { 820 // Cut-off check 821 if (_cutoff) 822 if (!_elepairs.BitIsSet(i)) 823 continue; 824 825 if (gradients) { 826 AddGradient(_electrostaticcalculations[i].force_a, _electrostaticcalculations[i].idx_a); 827 AddGradient(_electrostaticcalculations[i].force_b, _electrostaticcalculations[i].idx_b); 828 } 829 } 830 #endif 831 832 IF_OBFF_LOGLVL_MEDIUM { 833 snprintf(_logbuf, BUFF_SIZE, " TOTAL ELECTROSTATIC ENERGY = %8.5f %s\n", energy, GetUnit().c_str()); 834 OBFFLog(_logbuf); 835 } 836 837 return energy; 838 } 839 840 // 841 // OBForceFieldMMFF member functions 842 // 843 //*********************************************** 844 //Make a global instance 845 OBForceFieldMMFF94 theForceFieldMMFF94("MMFF94", false); 846 OBForceFieldMMFF94 theForceFieldMMFF94s("MMFF94s", false); 847 //*********************************************** 848 ~OBForceFieldMMFF94()849 OBForceFieldMMFF94::~OBForceFieldMMFF94() 850 { 851 } 852 operator =(OBForceFieldMMFF94 & src)853 OBForceFieldMMFF94 &OBForceFieldMMFF94::operator=(OBForceFieldMMFF94 &src) 854 { 855 _mol = src._mol; 856 _init = src._init; 857 return *this; 858 } 859 860 //////////////////////////////////////////////////////////////////////////////// 861 //////////////////////////////////////////////////////////////////////////////// 862 // 863 // Parse parameter files 864 // 865 //////////////////////////////////////////////////////////////////////////////// 866 //////////////////////////////////////////////////////////////////////////////// 867 ParseParamFile()868 bool OBForceFieldMMFF94::ParseParamFile() 869 { 870 // Set the locale for number parsing to avoid locale issues: PR#1785463 871 obLocale.SetLocale(); 872 873 vector<string> vs; 874 char buffer[80]; 875 876 // open data/_parFile 877 ifstream ifs; 878 if (OpenDatafile(ifs, _parFile).length() == 0) { 879 obErrorLog.ThrowError(__FUNCTION__, "Cannot open parameter file", obError); 880 return false; 881 } 882 883 while (ifs.getline(buffer, 80)) { 884 if (EQn(buffer, "#", 1)) continue; 885 886 tokenize(vs, buffer); 887 if (vs.size() < 2) 888 continue; 889 890 if (vs[0] == "prop") 891 ParseParamProp(vs[1]); 892 if (vs[0] == "def") 893 ParseParamDef(vs[1]); 894 if (vs[0] == "bond") 895 ParseParamBond(vs[1]); 896 if (vs[0] == "ang") 897 ParseParamAngle(vs[1]); 898 if (vs[0] == "bndk") 899 ParseParamBndk(vs[1]); 900 if (vs[0] == "chg") 901 ParseParamCharge(vs[1]); 902 if (vs[0] == "dfsb") 903 ParseParamDfsb(vs[1]); 904 if (vs[0] == "oop") 905 ParseParamOOP(vs[1]); 906 if (vs[0] == "pbci") 907 ParseParamPbci(vs[1]); 908 if (vs[0] == "stbn") 909 ParseParamStrBnd(vs[1]); 910 if (vs[0] == "tor") 911 ParseParamTorsion(vs[1]); 912 if (vs[0] == "vdw") 913 ParseParamVDW(vs[1]); 914 } 915 916 if (ifs) 917 ifs.close(); 918 919 // return the locale to the original one 920 obLocale.RestoreLocale(); 921 return true; 922 } 923 ParseParamBond(std::string & filename)924 bool OBForceFieldMMFF94::ParseParamBond(std::string &filename) 925 { 926 vector<string> vs; 927 char buffer[80]; 928 929 OBFFParameter parameter; 930 931 // open data/mmffbond.par 932 ifstream ifs; 933 if (OpenDatafile(ifs, filename).length() == 0) { 934 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffbond.par", obError); 935 return false; 936 } 937 938 while (ifs.getline(buffer, 80)) { 939 if (EQn(buffer, "*", 1)) continue; 940 if (EQn(buffer, "$", 1)) continue; 941 942 tokenize(vs, buffer); 943 944 parameter.clear(); 945 parameter._ipar.push_back(atoi(vs[0].c_str())); // FF class 946 parameter.a = atoi(vs[1].c_str()); 947 parameter.b = atoi(vs[2].c_str()); 948 parameter._dpar.push_back(atof(vs[3].c_str())); // kb 949 parameter._dpar.push_back(atof(vs[4].c_str())); // r0 950 _ffbondparams.push_back(parameter); 951 } 952 953 if (ifs) 954 ifs.close(); 955 956 return 0; 957 } 958 ParseParamBndk(std::string & filename)959 bool OBForceFieldMMFF94::ParseParamBndk(std::string &filename) 960 { 961 vector<string> vs; 962 char buffer[80]; 963 964 OBFFParameter parameter; 965 966 // open data/mmffbndk.par 967 ifstream ifs; 968 if (OpenDatafile(ifs, filename).length() == 0) { 969 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffbndk.par", obError); 970 return false; 971 } 972 973 while (ifs.getline(buffer, 80)) { 974 if (EQn(buffer, "*", 1)) continue; 975 if (EQn(buffer, "$", 1)) continue; 976 977 tokenize(vs, buffer); 978 979 parameter.clear(); 980 parameter.a = atoi(vs[0].c_str()); 981 parameter.b = atoi(vs[1].c_str()); 982 parameter._dpar.push_back(atof(vs[2].c_str())); // r0-ref 983 parameter._dpar.push_back(atof(vs[3].c_str())); // kb-ref 984 _ffbndkparams.push_back(parameter); 985 } 986 987 if (ifs) 988 ifs.close(); 989 990 return 0; 991 } 992 ParseParamAngle(std::string & filename)993 bool OBForceFieldMMFF94::ParseParamAngle(std::string &filename) 994 { 995 vector<string> vs; 996 char buffer[80]; 997 998 OBFFParameter parameter; 999 1000 // open data/mmffang.par 1001 ifstream ifs; 1002 if (OpenDatafile(ifs, filename).length() == 0) { 1003 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffang.par", obError); 1004 return false; 1005 } 1006 1007 while (ifs.getline(buffer, 80)) { 1008 if (EQn(buffer, "*", 1)) continue; 1009 if (EQn(buffer, "$", 1)) continue; 1010 1011 tokenize(vs, buffer); 1012 1013 parameter.clear(); 1014 parameter._ipar.push_back(atoi(vs[0].c_str())); // FF class 1015 parameter.a = atoi(vs[1].c_str()); 1016 parameter.b = atoi(vs[2].c_str()); 1017 parameter.c = atoi(vs[3].c_str()); 1018 parameter._dpar.push_back(atof(vs[4].c_str())); // ka 1019 parameter._dpar.push_back(atof(vs[5].c_str())); // theta0 1020 _ffangleparams.push_back(parameter); 1021 } 1022 1023 if (ifs) 1024 ifs.close(); 1025 1026 return 0; 1027 } 1028 ParseParamStrBnd(std::string & filename)1029 bool OBForceFieldMMFF94::ParseParamStrBnd(std::string &filename) 1030 { 1031 vector<string> vs; 1032 char buffer[80]; 1033 1034 OBFFParameter parameter; 1035 1036 // open data/mmffstbn.par 1037 ifstream ifs; 1038 if (OpenDatafile(ifs, filename).length() == 0) { 1039 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffstbn.par", obError); 1040 return false; 1041 } 1042 1043 while (ifs.getline(buffer, 80)) { 1044 if (EQn(buffer, "*", 1)) continue; 1045 if (EQn(buffer, "$", 1)) continue; 1046 1047 tokenize(vs, buffer); 1048 1049 parameter.clear(); 1050 parameter._ipar.push_back(atoi(vs[0].c_str())); // FF class 1051 parameter.a = atoi(vs[1].c_str()); 1052 parameter.b = atoi(vs[2].c_str()); 1053 parameter.c = atoi(vs[3].c_str()); 1054 parameter._dpar.push_back(atof(vs[4].c_str())); // kbaIJK 1055 parameter._dpar.push_back(atof(vs[5].c_str())); // kbaKJI 1056 _ffstrbndparams.push_back(parameter); 1057 } 1058 1059 if (ifs) 1060 ifs.close(); 1061 1062 return 0; 1063 } 1064 ParseParamDfsb(std::string & filename)1065 bool OBForceFieldMMFF94::ParseParamDfsb(std::string &filename) 1066 { 1067 vector<string> vs; 1068 char buffer[80]; 1069 1070 OBFFParameter parameter; 1071 1072 // open data/mmffdfsb.par 1073 ifstream ifs; 1074 if (OpenDatafile(ifs, filename).length() == 0) { 1075 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffdfsb.par", obError); 1076 return false; 1077 } 1078 1079 while (ifs.getline(buffer, 80)) { 1080 if (EQn(buffer, "*", 1)) continue; 1081 if (EQn(buffer, "$", 1)) continue; 1082 1083 tokenize(vs, buffer); 1084 1085 parameter.clear(); 1086 parameter.a = atoi(vs[0].c_str()); 1087 parameter.b = atoi(vs[1].c_str()); 1088 parameter.c = atoi(vs[2].c_str()); 1089 parameter._dpar.push_back(atof(vs[3].c_str())); // kbaIJK 1090 parameter._dpar.push_back(atof(vs[4].c_str())); // kbaKJI 1091 _ffdfsbparams.push_back(parameter); 1092 } 1093 1094 if (ifs) 1095 ifs.close(); 1096 1097 return 0; 1098 } 1099 ParseParamOOP(std::string & filename)1100 bool OBForceFieldMMFF94::ParseParamOOP(std::string &filename) 1101 { 1102 vector<string> vs; 1103 char buffer[80]; 1104 1105 OBFFParameter parameter; 1106 1107 // open data/mmffoop.par 1108 ifstream ifs; 1109 if (OpenDatafile(ifs, filename).length() == 0) { 1110 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffoop.par", obError); 1111 return false; 1112 } 1113 1114 while (ifs.getline(buffer, 80)) { 1115 if (EQn(buffer, "*", 1)) continue; 1116 if (EQn(buffer, "$", 1)) continue; 1117 1118 tokenize(vs, buffer); 1119 1120 parameter.clear(); 1121 parameter.a = atoi(vs[0].c_str()); 1122 parameter.b = atoi(vs[1].c_str()); 1123 parameter.c = atoi(vs[2].c_str()); 1124 parameter.d = atoi(vs[3].c_str()); 1125 parameter._dpar.push_back(atof(vs[4].c_str())); // koop 1126 _ffoopparams.push_back(parameter); 1127 } 1128 1129 if (ifs) 1130 ifs.close(); 1131 1132 return 0; 1133 } 1134 ParseParamTorsion(std::string & filename)1135 bool OBForceFieldMMFF94::ParseParamTorsion(std::string &filename) 1136 { 1137 vector<string> vs; 1138 char buffer[80]; 1139 1140 OBFFParameter parameter; 1141 1142 // open data/mmfftor.par 1143 ifstream ifs; 1144 if (OpenDatafile(ifs, filename).length() == 0) { 1145 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmfftor.par", obError); 1146 return false; 1147 } 1148 1149 while (ifs.getline(buffer, 80)) { 1150 if (EQn(buffer, "*", 1)) continue; 1151 if (EQn(buffer, "$", 1)) continue; 1152 1153 tokenize(vs, buffer); 1154 1155 parameter.clear(); 1156 parameter._ipar.push_back(atoi(vs[0].c_str())); // FF class 1157 parameter.a = atoi(vs[1].c_str()); 1158 parameter.b = atoi(vs[2].c_str()); 1159 parameter.c = atoi(vs[3].c_str()); 1160 parameter.d = atoi(vs[4].c_str()); 1161 parameter._dpar.push_back(atof(vs[5].c_str())); // v1 1162 parameter._dpar.push_back(atof(vs[6].c_str())); // v2 1163 parameter._dpar.push_back(atof(vs[7].c_str())); // v3 1164 _fftorsionparams.push_back(parameter); 1165 } 1166 1167 if (ifs) 1168 ifs.close(); 1169 1170 return 0; 1171 } 1172 ParseParamVDW(std::string & filename)1173 bool OBForceFieldMMFF94::ParseParamVDW(std::string &filename) 1174 { 1175 vector<string> vs; 1176 char buffer[80]; 1177 1178 OBFFParameter parameter; 1179 1180 // open data/mmffvdw.par 1181 ifstream ifs; 1182 if (OpenDatafile(ifs, filename).length() == 0) { 1183 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffvdw.par", obError); 1184 return false; 1185 } 1186 1187 while (ifs.getline(buffer, 80)) { 1188 if (EQn(buffer, "*", 1)) continue; 1189 if (EQn(buffer, "$", 1)) continue; 1190 1191 tokenize(vs, buffer); 1192 1193 parameter.clear(); 1194 parameter.a = atoi(vs[0].c_str()); 1195 parameter._dpar.push_back(atof(vs[1].c_str())); // alpha-i 1196 parameter._dpar.push_back(atof(vs[2].c_str())); // N-i 1197 parameter._dpar.push_back(atof(vs[3].c_str())); // A-i 1198 parameter._dpar.push_back(atof(vs[4].c_str())); // G-i 1199 if (EQn(vs[5].c_str(), "-", 1)) 1200 parameter._ipar.push_back(0); 1201 else if (EQn(vs[5].c_str(), "D", 1)) 1202 parameter._ipar.push_back(1); // hydrogen bond donor 1203 else if (EQn(vs[5].c_str(), "A", 1)) 1204 parameter._ipar.push_back(2); // hydrogen bond acceptor 1205 _ffvdwparams.push_back(parameter); 1206 } 1207 1208 if (ifs) 1209 ifs.close(); 1210 1211 return 0; 1212 } 1213 ParseParamCharge(std::string & filename)1214 bool OBForceFieldMMFF94::ParseParamCharge(std::string &filename) 1215 { 1216 vector<string> vs; 1217 char buffer[80]; 1218 1219 OBFFParameter parameter; 1220 1221 // open data/mmffchg.par 1222 ifstream ifs; 1223 if (OpenDatafile(ifs, filename).length() == 0) { 1224 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffchg.par", obError); 1225 return false; 1226 } 1227 1228 while (ifs.getline(buffer, 80)) { 1229 if (EQn(buffer, "*", 1)) continue; 1230 if (EQn(buffer, "$", 1)) continue; 1231 1232 tokenize(vs, buffer); 1233 1234 parameter.clear(); 1235 parameter._ipar.push_back(atoi(vs[0].c_str())); // FF class 1236 parameter.a = atoi(vs[1].c_str()); 1237 parameter.b = atoi(vs[2].c_str()); 1238 parameter._dpar.push_back(atof(vs[3].c_str())); // bci 1239 _ffchgparams.push_back(parameter); 1240 } 1241 1242 if (ifs) 1243 ifs.close(); 1244 1245 return 0; 1246 } 1247 ParseParamPbci(std::string & filename)1248 bool OBForceFieldMMFF94::ParseParamPbci(std::string &filename) 1249 { 1250 vector<string> vs; 1251 char buffer[80]; 1252 1253 OBFFParameter parameter; 1254 1255 // open data/mmffpbci.par 1256 ifstream ifs; 1257 if (OpenDatafile(ifs, filename).length() == 0) { 1258 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffpbci", obError); 1259 return false; 1260 } 1261 1262 while (ifs.getline(buffer, 80)) { 1263 if (EQn(buffer, "*", 1)) continue; 1264 if (EQn(buffer, "$", 1)) continue; 1265 1266 tokenize(vs, buffer); 1267 1268 parameter.clear(); 1269 parameter.a = atoi(vs[1].c_str()); 1270 parameter._dpar.push_back(atof(vs[2].c_str())); // pbci 1271 parameter._dpar.push_back(atof(vs[3].c_str())); // fcadj 1272 _ffpbciparams.push_back(parameter); 1273 } 1274 1275 if (ifs) 1276 ifs.close(); 1277 1278 return 0; 1279 } 1280 ParseParamProp(std::string & filename)1281 bool OBForceFieldMMFF94::ParseParamProp(std::string &filename) 1282 { 1283 vector<string> vs; 1284 char buffer[80]; 1285 1286 OBFFParameter parameter; 1287 1288 // open data/mmffprop.par 1289 ifstream ifs; 1290 if (OpenDatafile(ifs, filename).length() == 0) { 1291 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffprop.par", obError); 1292 return false; 1293 } 1294 1295 while (ifs.getline(buffer, 80)) { 1296 if (EQn(buffer, "*", 1)) continue; 1297 if (EQn(buffer, "$", 1)) continue; 1298 1299 tokenize(vs, buffer); 1300 1301 parameter.clear(); 1302 parameter.a = atoi(vs[0].c_str()); // atom type 1303 parameter._ipar.push_back(atoi(vs[1].c_str())); // at.no 1304 parameter._ipar.push_back(atoi(vs[2].c_str())); // crd 1305 parameter._ipar.push_back(atoi(vs[3].c_str())); // val 1306 parameter._ipar.push_back(atoi(vs[4].c_str())); // pilp 1307 parameter._ipar.push_back(atoi(vs[5].c_str())); // mltb 1308 parameter._ipar.push_back(atoi(vs[6].c_str())); // arom 1309 parameter._ipar.push_back(atoi(vs[7].c_str())); // linh 1310 parameter._ipar.push_back(atoi(vs[8].c_str())); // sbmb 1311 1312 if (parameter._ipar[3]) 1313 _ffpropPilp.SetBitOn(parameter.a); 1314 if (parameter._ipar[5]) 1315 _ffpropArom.SetBitOn(parameter.a); 1316 if (parameter._ipar[6]) 1317 _ffpropLin.SetBitOn(parameter.a); 1318 if (parameter._ipar[7]) 1319 _ffpropSbmb.SetBitOn(parameter.a); 1320 1321 _ffpropparams.push_back(parameter); 1322 } 1323 1324 if (ifs) 1325 ifs.close(); 1326 1327 return 0; 1328 } 1329 ParseParamDef(std::string & filename)1330 bool OBForceFieldMMFF94::ParseParamDef(std::string &filename) 1331 { 1332 vector<string> vs; 1333 char buffer[80]; 1334 1335 OBFFParameter parameter; 1336 1337 // open data/mmffdef.par 1338 ifstream ifs; 1339 if (OpenDatafile(ifs, filename).length() == 0) { 1340 obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffdef.par", obError); 1341 return false; 1342 } 1343 1344 while (ifs.getline(buffer, 80)) { 1345 if (EQn(buffer, "*", 1)) continue; 1346 if (EQn(buffer, "$", 1)) continue; 1347 1348 tokenize(vs, buffer); 1349 1350 parameter.clear(); 1351 parameter._ipar.push_back(atoi(vs[1].c_str())); // level 1 1352 parameter._ipar.push_back(atoi(vs[2].c_str())); // level 2 1353 parameter._ipar.push_back(atoi(vs[3].c_str())); // level 3 1354 parameter._ipar.push_back(atoi(vs[4].c_str())); // level 4 1355 parameter._ipar.push_back(atoi(vs[5].c_str())); // level 5 1356 _ffdefparams.push_back(parameter); 1357 } 1358 1359 if (ifs) 1360 ifs.close(); 1361 1362 return 0; 1363 } 1364 1365 //////////////////////////////////////////////////////////////////////////////// 1366 //////////////////////////////////////////////////////////////////////////////// 1367 // 1368 // Setup Functions 1369 // 1370 //////////////////////////////////////////////////////////////////////////////// 1371 //////////////////////////////////////////////////////////////////////////////// 1372 1373 // The MMFF94 article doesn't seem to include information about how 1374 // aromaticity is perceived. This function was written by studying the 1375 // MMFF_opti.log file, trail-and-error and using the MMFF94 validation 1376 // set to check the results (If all atom types are assigned correctly, 1377 // aromatic rings are probably detected correctly) PerceiveAromatic()1378 bool OBForceFieldMMFF94::PerceiveAromatic() 1379 { 1380 bool done = false; // not done actually.... 1381 OBAtom *ringatom; 1382 OBBond *ringbond; 1383 vector<OBRing*> vr; 1384 vr = _mol.GetSSSR(); 1385 1386 vector<OBRing*>::iterator ri; 1387 vector<int>::iterator rj; 1388 int n, index, ringsize, first_rj, prev_rj, pi_electrons, c60; 1389 for (ri = vr.begin();ri != vr.end();++ri) { // for each ring 1390 ringsize = (*ri)->Size(); 1391 1392 n = 1; 1393 pi_electrons = 0; 1394 c60 = 0; // we have a special case to get c60 right (all atom type 37) 1395 for(rj = (*ri)->_path.begin();rj != (*ri)->_path.end();++rj) { // for each ring atom 1396 index = *rj; 1397 ringatom = _mol.GetAtom(index); 1398 1399 // is the bond to the previous ring atom double? 1400 if (n > 1) { 1401 ringbond = _mol.GetBond(prev_rj, index); 1402 if (!ringbond) { 1403 prev_rj = index; 1404 continue; 1405 } 1406 if (ringbond->GetBondOrder() == 2) { 1407 pi_electrons += 2; 1408 prev_rj = index; 1409 n++; 1410 continue; 1411 } 1412 prev_rj = index; 1413 } else { 1414 prev_rj = index; 1415 first_rj = index; 1416 } 1417 1418 // does the current ring atom have a exocyclic double bond? 1419 FOR_NBORS_OF_ATOM (nbr, ringatom) { 1420 if ((*ri)->IsInRing(nbr->GetIdx())) 1421 continue; 1422 1423 if (!nbr->IsAromatic()) { 1424 if (ringatom->GetAtomicNum() == OBElements::Carbon && ringatom->IsInRingSize(5) 1425 && ringatom->IsInRingSize(6) && nbr->GetAtomicNum() == OBElements::Carbon && nbr->IsInRingSize(5) 1426 && nbr->IsInRingSize(6)) { 1427 c60++; 1428 } else { 1429 continue; 1430 } 1431 } 1432 1433 ringbond = _mol.GetBond(nbr->GetIdx(), index); 1434 if (!ringbond) { 1435 continue; 1436 } 1437 if (ringbond->GetBondOrder() == 2) 1438 pi_electrons++; 1439 } 1440 1441 // is the atom N, O or S in 5 rings 1442 if (ringsize == 5 && 1443 ringatom->GetIdx() == (*ri)->GetRootAtom()) 1444 pi_electrons += 2; 1445 1446 n++; 1447 1448 } // for each ring atom 1449 1450 // is the bond from the first to the last atom double? 1451 ringbond = _mol.GetBond(first_rj, index); 1452 if (ringbond) { 1453 if (ringbond->GetBondOrder() == 2) 1454 pi_electrons += 2; 1455 } 1456 1457 if (((pi_electrons == 6) && ((ringsize == 5) || (ringsize == 6))) 1458 || ((pi_electrons == 5) && (c60 == 5))) { 1459 // mark ring atoms as aromatic 1460 for(rj = (*ri)->_path.begin();rj != (*ri)->_path.end();++rj) { 1461 if (!_mol.GetAtom(*rj)->IsAromatic()) 1462 done = true; 1463 _mol.GetAtom(*rj)->SetAromatic(); 1464 } 1465 // mark all ring bonds as aromatic 1466 FOR_BONDS_OF_MOL (bond, _mol) 1467 if((*ri)->IsMember(&*bond)) 1468 bond->SetAromatic(); 1469 } 1470 } 1471 1472 return done; 1473 } 1474 1475 // Symbolic atom typing is skipped 1476 // 1477 // atom typing is based on: 1478 // MMFF94 I - Table III 1479 // MMFF94 V - Table I 1480 // GetType(OBAtom * atom)1481 int OBForceFieldMMFF94::GetType(OBAtom *atom) 1482 { 1483 OBBond *bond; 1484 int oxygenCount, nitrogenCount, sulphurCount, doubleBondTo; 1485 //////////////////////////////// 1486 // Aromatic Atoms 1487 //////////////////////////////// 1488 if (atom->IsAromatic()) { 1489 if (atom->IsInRingSize(5)) { 1490 bool IsAromatic = false; 1491 vector<OBAtom*> alphaPos, betaPos; 1492 vector<OBAtom*> alphaAtoms, betaAtoms; 1493 1494 if (atom->GetAtomicNum() == OBElements::Sulfur) { 1495 return 44; // Aromatic 5-ring sulfur with pi lone pair (STHI) 1496 } 1497 if (atom->GetAtomicNum() == OBElements::Oxygen) { 1498 return 59; // Aromatic 5-ring oxygen with pi lone pair (OFUR) 1499 } 1500 if (atom->GetAtomicNum() == OBElements::Nitrogen) { 1501 FOR_NBORS_OF_ATOM (nbr, atom) { 1502 if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) { 1503 return 82; // N-oxide nitrogen in 5-ring alpha position, 1504 // N-oxide nitrogen in 5-ring beta position, 1505 // N-oxide nitrogen in other 5-ring position, 1506 // (N5AX, N5BX, N5OX) 1507 } 1508 } 1509 } 1510 FOR_NBORS_OF_ATOM (nbr, atom) { 1511 if (!((_mol.GetBond(atom, &*nbr))->IsAromatic()) || !nbr->IsInRingSize(5)) 1512 continue; 1513 1514 if (IsInSameRing(atom, &*nbr)) { 1515 alphaPos.push_back(&*nbr); 1516 } 1517 1518 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 1519 if (nbrNbr->GetIdx() == atom->GetIdx()) 1520 continue; 1521 if (!((_mol.GetBond(&*nbr, &*nbrNbr))->IsAromatic()) || !nbrNbr->IsInRingSize(5)) 1522 continue; 1523 1524 IsAromatic = true; 1525 1526 if (IsInSameRing(atom, &*nbrNbr)) { 1527 betaPos.push_back(&*nbrNbr); 1528 } 1529 } 1530 } 1531 1532 if (IsAromatic) { 1533 1534 1535 for (unsigned int i = 0; i < alphaPos.size(); i++) { 1536 if (alphaPos[i]->GetAtomicNum() == OBElements::Sulfur) { 1537 alphaAtoms.push_back(alphaPos[i]); 1538 } else if (alphaPos[i]->GetAtomicNum() == OBElements::Oxygen) { 1539 alphaAtoms.push_back(alphaPos[i]); 1540 } else if (alphaPos[i]->GetAtomicNum() == OBElements::Nitrogen && (alphaPos[i]->GetExplicitDegree() == 3)) { 1541 bool IsNOxide = false; 1542 FOR_NBORS_OF_ATOM (nbr, alphaPos[i]) { 1543 if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) { 1544 IsNOxide = true; 1545 } 1546 } 1547 1548 if (!IsNOxide) { 1549 alphaAtoms.push_back(alphaPos[i]); 1550 } 1551 } 1552 } 1553 for (unsigned int i = 0; i < betaPos.size(); i++) { 1554 if (betaPos[i]->GetAtomicNum() == OBElements::Sulfur) { 1555 betaAtoms.push_back(betaPos[i]); 1556 } else if (betaPos[i]->GetAtomicNum() == OBElements::Oxygen) { 1557 betaAtoms.push_back(betaPos[i]); 1558 } else if (betaPos[i]->GetAtomicNum() == OBElements::Nitrogen && (betaPos[i]->GetExplicitDegree() == 3)) { 1559 bool IsNOxide = false; 1560 FOR_NBORS_OF_ATOM (nbr, betaPos[i]) { 1561 if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) { 1562 IsNOxide = true; 1563 } 1564 } 1565 1566 if (!IsNOxide) { 1567 betaAtoms.push_back(betaPos[i]); 1568 } 1569 } 1570 } 1571 if (!betaAtoms.size()) { 1572 nitrogenCount = 0; 1573 FOR_NBORS_OF_ATOM (nbr, atom) { 1574 if (nbr->GetAtomicNum() == OBElements::Nitrogen && (nbr->GetExplicitDegree() == 3)) { 1575 if ((nbr->GetExplicitValence() == 4) && nbr->IsAromatic()) { 1576 nitrogenCount++; 1577 } else if ((nbr->GetExplicitValence() == 3) && !nbr->IsAromatic()) { 1578 nitrogenCount++; 1579 } 1580 } 1581 } 1582 if (nitrogenCount >= 2) { 1583 return 80; // Aromatic carbon between N's in imidazolium (CIM+) 1584 } 1585 } 1586 if (!alphaAtoms.size() && !betaAtoms.size()) { 1587 if (atom->GetAtomicNum() == OBElements::Carbon) { 1588 bool c60 = true; // special case to ensure c60 is typed correctly -- Paolo Tosco 1589 FOR_NBORS_OF_ATOM (nbr, atom) { 1590 if (!(nbr->GetAtomicNum() == OBElements::Carbon && nbr->IsAromatic() && nbr->IsInRingSize(6))) 1591 c60 = false; 1592 } 1593 if (c60) 1594 return 37; // correct atom type for c in c60 (all atoms symmetric) 1595 // there is no S:, O:, or N: 1596 // this is the case for anions with only carbon and nitrogen in the ring 1597 return 78; // General carbon in 5-membered aromatic ring (C5) 1598 } else if (atom->GetAtomicNum() == OBElements::Nitrogen) { 1599 if (atom->GetExplicitDegree() == 3) { 1600 // this is the N: atom 1601 return 39; // Aromatic 5 ring nitrogen with pi lone pair (NPYL) 1602 } else { 1603 // again, no S:, O:, or N: 1604 return 76; // Nitrogen in 5-ring aromatic anion (N5M) 1605 } 1606 } 1607 } 1608 if (alphaAtoms.size() == 2) { 1609 if (atom->GetAtomicNum() == OBElements::Carbon && IsInSameRing(alphaAtoms[0], alphaAtoms[1])) { 1610 if (alphaAtoms[0]->GetAtomicNum() == OBElements::Nitrogen && alphaAtoms[1]->GetAtomicNum() == OBElements::Nitrogen) { 1611 if ((alphaAtoms[0]->GetExplicitDegree() == 3) && (alphaAtoms[1]->GetExplicitDegree() == 3)) { 1612 return 80; // Aromatic carbon between N's in imidazolium (CIM+) 1613 } 1614 } 1615 } 1616 } 1617 if (alphaAtoms.size() && !betaAtoms.size()) { 1618 if (atom->GetAtomicNum() == OBElements::Carbon) { 1619 return 63; // Aromatic 5-ring C, alpha to N:, O:, or S: (C5A) 1620 } else if (atom->GetAtomicNum() == OBElements::Nitrogen) { 1621 if (atom->GetExplicitDegree() == 3) { 1622 return 81; // Posivite nitrogen in 5-ring alpha position (N5A+) 1623 } else { 1624 return 65; // Aromatic 5-ring N, alpha to N:, O:, or S: (N5A) 1625 } 1626 } 1627 } 1628 if (!alphaAtoms.size() && betaAtoms.size()) { 1629 if (atom->GetAtomicNum() == OBElements::Carbon) { 1630 return 64; // Aromatic 5-ring C, beta to N:, O:, or S: (C5B) 1631 } else if (atom->GetAtomicNum() == OBElements::Nitrogen) { 1632 if (atom->GetExplicitDegree() == 3) { 1633 return 81; // Posivite nitrogen in 5-ring beta position (N5B+) 1634 } else { 1635 return 66; // Aromatic 5-ring N, beta to N:, O:, or S: (N5B) 1636 } 1637 } 1638 } 1639 if (alphaAtoms.size() && betaAtoms.size()) { 1640 for (unsigned int i = 0; i < alphaAtoms.size(); i++) { 1641 for (unsigned int j = 0; j < betaAtoms.size(); j++) { 1642 if (!IsInSameRing(alphaAtoms[i], betaAtoms[j])) { 1643 if (atom->GetAtomicNum() == OBElements::Carbon) { 1644 return 78; // General carbon in 5-membered aromatic ring (C5) 1645 } else if (atom->GetAtomicNum() == OBElements::Nitrogen) { 1646 return 79; // General nitrogen in 5-membered aromatic ring (N5) 1647 } 1648 } 1649 } 1650 } 1651 for (unsigned int i = 0; i < alphaAtoms.size(); i++) { 1652 if (alphaAtoms[i]->GetAtomicNum() == OBElements::Sulfur || alphaAtoms[i]->GetAtomicNum() == OBElements::Oxygen) { 1653 if (atom->GetAtomicNum() == OBElements::Carbon) { 1654 return 63; // Aromatic 5-ring C, alpha to N:, O:, or S: (C5A) 1655 } else if (atom->GetAtomicNum() == OBElements::Nitrogen) { 1656 return 65; // Aromatic 5-ring N, alpha to N:, O:, or S: (N5A) 1657 } 1658 } 1659 } 1660 for (unsigned int i = 0; i < betaAtoms.size(); i++) { 1661 if (betaAtoms[i]->GetAtomicNum() == OBElements::Sulfur || betaAtoms[i]->GetAtomicNum() == OBElements::Oxygen) { 1662 if (atom->GetAtomicNum() == OBElements::Carbon) { 1663 return 64; // Aromatic 5-ring C, beta to N:, O:, or S: (C5B) 1664 } else if (atom->GetAtomicNum() == OBElements::Nitrogen) { 1665 return 66; // Aromatic 5-ring N, beta to N:, O:, or S: (N5B) 1666 } 1667 } 1668 } 1669 1670 if (atom->GetAtomicNum() == OBElements::Carbon) { 1671 return 78; // General carbon in 5-membered aromatic ring (C5) 1672 } else if (atom->GetAtomicNum() == OBElements::Nitrogen) { 1673 return 79; // General nitrogen in 5-membered aromatic ring (N5) 1674 } 1675 } 1676 } 1677 } 1678 1679 if (atom->IsInRingSize(6)) { 1680 1681 if (atom->GetAtomicNum() == OBElements::Carbon) { 1682 return 37; // Aromatic carbon, e.g., in benzene (CB) 1683 } else if (atom->GetAtomicNum() == OBElements::Nitrogen) { 1684 FOR_NBORS_OF_ATOM (nbr, atom) { 1685 if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) { 1686 return 69; // Pyridinium N-oxide nitrogen (NPOX) 1687 } 1688 } 1689 1690 if (atom->GetExplicitDegree() == 3) { 1691 return 58; // Aromatic nitrogen in pyridinium (NPD+) 1692 } else { 1693 return 38; // Aromatic nitrogen with sigma lone pair (NPYD) 1694 } 1695 } 1696 } 1697 } 1698 1699 //////////////////////////////// 1700 // Hydrogen 1701 //////////////////////////////// 1702 if (atom->GetAtomicNum() == 1) { 1703 FOR_NBORS_OF_ATOM (nbr, atom) { 1704 if (nbr->GetAtomicNum() == OBElements::Carbon) { 1705 return 5; // Hydrogen attatched to carbon (HC) 1706 } 1707 if (nbr->GetAtomicNum() == 14) { 1708 return 5; // Hydrogen attatched to silicon (HSI) 1709 } 1710 if (nbr->GetAtomicNum() == OBElements::Oxygen) { 1711 if (nbr->GetExplicitValence() == 3) { 1712 if (nbr->GetExplicitDegree() == 3) { 1713 return 50; // Hydrogen on oxonium oxygen (HO+) 1714 } else { 1715 return 52; // Hydrogen on oxenium oxygen (HO=+) 1716 } 1717 } 1718 1719 int hydrogenCount = 0; 1720 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 1721 if (nbrNbr->GetAtomicNum() == OBElements::Hydrogen) { 1722 hydrogenCount++; 1723 continue; 1724 } 1725 if (nbrNbr->GetAtomicNum() == OBElements::Carbon) { 1726 if (nbrNbr->IsAromatic()) { 1727 return 29; // phenol 1728 } 1729 1730 FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) { 1731 if (nbrNbrNbr->GetIdx() == nbr->GetIdx()) 1732 continue; 1733 1734 bond = _mol.GetBond(&*nbrNbr, &*nbrNbrNbr); 1735 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 1736 if (nbrNbrNbr->GetAtomicNum() == OBElements::Oxygen) { 1737 return 24; // Hydroxyl hydrogen in carboxylic acids (HOCO) 1738 } 1739 if (nbrNbrNbr->GetAtomicNum() == OBElements::Carbon || nbrNbrNbr->GetAtomicNum() == OBElements::Nitrogen) { 1740 return 29; // Enolic or phenolic hydroxyl hydrogen, 1741 // Hydroxyl hydrogen in HO-C=N moiety (HOCC, HOCN) 1742 } 1743 } 1744 } 1745 } 1746 if (nbrNbr->GetAtomicNum() == OBElements::Phosphorus) { 1747 return 24; // Hydroxyl hydrogen in H-O-P moiety (HOP) 1748 } 1749 if (nbrNbr->GetAtomicNum() == OBElements::Sulfur) { 1750 return 33; // Hydrogen on oxygen attached to sulfur (HOS) 1751 } 1752 1753 } 1754 if (hydrogenCount == 2) { 1755 return 31; // Hydroxyl hydrogen in water (HOH) 1756 } 1757 1758 return 21; // Hydroxyl hydrogen in alcohols, Generic hydroxyl hydrogen (HOR, HO) 1759 } 1760 if (nbr->GetAtomicNum() == OBElements::Nitrogen) { 1761 switch (GetType(&*nbr)) { 1762 case 81: 1763 return 36; // Hydrogen on imidazolium nitrogen (HIM+) 1764 case 68: 1765 return 23; // Hydrogen on N in N-oxide (HNOX) 1766 case 67: 1767 return 23; // Hydrogen on N in N-oxide (HNOX) 1768 case 62: 1769 return 23; // Generic hydrogen on sp3 nitrogen, e.g., in amines (HNR) 1770 case 56: 1771 return 36; // Hydrogen on guanimdinium nitrogen (HGD+) 1772 case 55: 1773 return 36; // Hydrogen on amidinium nitrogen (HNN+) 1774 case 43: 1775 return 28; // Hydrogen on NSO, NSO2, or NSO3 nitrogen, Hydrogen on N triply bonded to C (HNSO, HNC%) 1776 case 39: 1777 return 23; // Hydrogen on nitrogen in pyrrole (HPYL) 1778 case 8: 1779 return 23; // Generic hydrogen on sp3 nitrogen, e.g., in amines, Hydrogen on nitrogen in ammonia (HNR, H3N) 1780 } 1781 1782 if (nbr->GetExplicitValence() == 4) { 1783 if (nbr->GetExplicitDegree() == 2) { 1784 return 28; // Hydrogen on N triply bonded to C (HNC%) 1785 } else { 1786 return 36; // Hydrogen on pyridinium nitrogen, Hydrogen on protonated imine nitrogen (HPD+, HNC+) 1787 } 1788 } 1789 1790 if (nbr->GetExplicitDegree() == 2) { 1791 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 1792 if (nbrNbr->GetAtomicNum() == OBElements::Hydrogen) 1793 continue; 1794 1795 bond = _mol.GetBond(&*nbr, &*nbrNbr); 1796 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 1797 if (nbrNbr->GetAtomicNum() == OBElements::Carbon || nbrNbr->GetAtomicNum() == OBElements::Nitrogen) { 1798 return 27; // Hydrogen on imine nitrogen, Hydrogen on azo nitrogen (HN=C, HN=N) 1799 } 1800 1801 return 28; // Generic hydrogen on sp2 nitrogen (HSP2) 1802 } 1803 } 1804 } 1805 1806 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 1807 if (nbrNbr->GetAtomicNum() == OBElements::Hydrogen) 1808 continue; 1809 1810 if (nbrNbr->GetAtomicNum() == OBElements::Carbon) { 1811 if (nbrNbr->IsAromatic()) { 1812 return 28; // deloc. lp pair 1813 } 1814 1815 FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) { 1816 if (nbrNbrNbr->GetIdx() == nbr->GetIdx()) 1817 continue; 1818 1819 bond = _mol.GetBond(&*nbrNbr, &*nbrNbrNbr); 1820 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 1821 if (nbrNbrNbr->GetAtomicNum() == OBElements::Carbon || nbrNbrNbr->GetAtomicNum() == OBElements::Nitrogen || nbrNbrNbr->GetAtomicNum() == OBElements::Oxygen || nbrNbrNbr->GetAtomicNum() == OBElements::Sulfur) { 1822 return 28; // Hydrogen on amide nitrogen, Hydrogen on thioamide nitrogen, 1823 // Hydrogen on enamine nitrogen, Hydrogen in H-N-C=N moiety (HNCO, HNCS, HNCC, HNCN) 1824 } 1825 } 1826 } 1827 } 1828 if (nbrNbr->GetAtomicNum() == OBElements::Nitrogen) { 1829 FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) { 1830 if (nbrNbrNbr->GetIdx() == nbr->GetIdx()) 1831 continue; 1832 1833 bond = _mol.GetBond(&*nbrNbr, &*nbrNbrNbr); 1834 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 1835 if (nbrNbrNbr->GetAtomicNum() == OBElements::Carbon || nbrNbrNbr->GetAtomicNum() == OBElements::Nitrogen) { 1836 return 28; // Hydrogen in H-N-N=C moiety, Hydrogen in H-N-N=N moiety (HNNC, HNNN) 1837 } 1838 } 1839 } 1840 } 1841 if (nbrNbr->GetAtomicNum() == OBElements::Sulfur) { 1842 FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) { 1843 if (nbrNbrNbr->GetIdx() == nbr->GetIdx()) 1844 continue; 1845 1846 if (nbrNbrNbr->GetAtomicNum() == OBElements::Oxygen || (nbrNbrNbr->GetExplicitDegree() == 1)) { 1847 return 28; // Hydrogen on NSO, NSO2 or NSO3 nitrogen (HNSO) 1848 } 1849 } 1850 } 1851 } 1852 1853 return 23; // Generic hydrogen on sp3 nitrogen e.g., in amines, 1854 // Hydrogen on nitrogen in pyrrole, Hydrogen in ammonia, 1855 // Hydrogen on N in N-oxide (HNR, HPYL, H3N, HNOX) 1856 } 1857 if (nbr->GetAtomicNum() == OBElements::Sulfur || nbr->GetAtomicNum() == OBElements::Phosphorus) { 1858 return 71; // Hydrogen attached to sulfur, Hydrogen attached to >S= sulfur doubly bonded to N, 1859 // Hydrogen attached to phosphorus (HS, HS=N, HP) 1860 } 1861 } 1862 } 1863 1864 //////////////////////////////// 1865 // Lithium 1866 //////////////////////////////// 1867 if (atom->GetAtomicNum() == 3) { 1868 // 0 neighbours 1869 if (atom->GetExplicitDegree() == 0) { 1870 return 92; // Lithium cation (LI+) 1871 } 1872 } 1873 1874 //////////////////////////////// 1875 // Carbon 1876 //////////////////////////////// 1877 if (atom->GetAtomicNum() == 6) { 1878 // 4 neighbours 1879 if (atom->GetExplicitDegree() == 4) { 1880 if (atom->IsInRingSize(3)) { 1881 return 22; // Aliphatic carbon in 3-membered ring (CR3R) 1882 } 1883 1884 if (atom->IsInRingSize(4)) { 1885 return 20; // Aliphatic carbon in 4-membered ring (CR4R) 1886 } 1887 1888 return 1; // Alkyl carbon (CR) 1889 } 1890 // 3 neighbours 1891 if (atom->GetExplicitDegree() == 3) { 1892 int N2count = 0; 1893 int N3count = 0; 1894 int N3fcharge = 0; 1895 oxygenCount = sulphurCount = doubleBondTo = 0; 1896 1897 FOR_NBORS_OF_ATOM (nbr, atom) { 1898 bond = _mol.GetBond(&*nbr, atom); 1899 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 1900 doubleBondTo = nbr->GetAtomicNum(); 1901 } 1902 1903 if (nbr->GetExplicitDegree() == 1) { 1904 if (nbr->GetAtomicNum() == OBElements::Oxygen) { 1905 oxygenCount++; 1906 } else if (nbr->GetAtomicNum() == OBElements::Sulfur) { 1907 sulphurCount++; 1908 } 1909 } else if (nbr->GetExplicitDegree() == 3) { 1910 if (nbr->GetAtomicNum() == OBElements::Nitrogen) { 1911 N3fcharge += nbr->GetFormalCharge(); 1912 N3count++; 1913 } 1914 } else if ((nbr->GetExplicitDegree() == 2) && !bond->IsAromatic() && bond->GetBondOrder() == 2) { 1915 if (nbr->GetAtomicNum() == OBElements::Nitrogen) { 1916 N2count++; 1917 } 1918 } 1919 } 1920 if ((N3count >= 2) && (doubleBondTo == 7 || (!(doubleBondTo == 6) && atom->GetExplicitDegree() == 3 1921 && N3fcharge == 1)) && !N2count && !oxygenCount && !sulphurCount) { 1922 // N3==C--N3 1923 return 57; // Guanidinium carbon, Carbon in +N=C-N: resonance structures (CGD+, CNN+) 1924 } 1925 if ((oxygenCount == 2) || (sulphurCount == 2)) { 1926 // O1-?-C-?-O1 or S1-?-C-?-S1 1927 return 41; // Carbon in carboxylate anion, Carbon in thiocarboxylate anion (CO2M, CS2M) 1928 } 1929 if (atom->IsInRingSize(4) && (doubleBondTo == 6)) { 1930 return 30; // Olefinic carbon in 4-membered ring (CR4E) 1931 } 1932 if ((doubleBondTo == 7) || (doubleBondTo == 8) || 1933 (doubleBondTo == 15) || (doubleBondTo == 16)) { 1934 // C==N, C==O, C==P, C==S 1935 return 3; // Generic carbonyl carbon, Imine-type carbon, Guanidine carbon, 1936 // Ketone or aldehyde carbonyl carbon, Amide carbonyl carbon, 1937 // Carboxylic acid or ester carbonyl carbon, Carbamate carbonyl carbon, 1938 // Carbonic acid or ester carbonyl carbon, Thioester carbonyl (double 1939 // bonded to O or S), Thioamide carbon (double bonded to S), Carbon 1940 // in >C=SO2, Sulfinyl carbon in >C=S=O, Thiocarboxylic acid or ester 1941 // carbon, Carbon doubly bonded to P (C=O, C=N, CGD, C=OR, C=ON, COO, 1942 // COON, COOO, C=OS, C=S, C=SN, CSO2, CS=O, CSS, C=P) 1943 } 1944 1945 return 2; // Vinylic Carbon, Generic sp2 carbon (C=C, CSP2) 1946 1947 } 1948 // 2 neighbours 1949 if (atom->GetExplicitDegree() == 2) { 1950 return 4; // Acetylenic carbon, Allenic caron (CSP, =C=) 1951 } 1952 // 1 neighbours 1953 if (atom->GetExplicitDegree() == 1) { 1954 return 60; // Isonitrile carbon (C%-) 1955 } 1956 } 1957 1958 //////////////////////////////// 1959 // Nitrogen 1960 //////////////////////////////// 1961 if (atom->GetAtomicNum() == 7) { 1962 // 4 neighbours 1963 if (atom->GetExplicitDegree() == 4) { 1964 FOR_NBORS_OF_ATOM (nbr, atom) { 1965 if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) { 1966 return 68; // sp3-hybridized N-oxide nitrogen (N3OX) 1967 } 1968 } 1969 1970 return 34; // Quaternary nitrogen (NR+) 1971 } 1972 // 3 neighbours 1973 if (atom->GetExplicitDegree() == 3) { 1974 if (atom->GetExplicitValence() >= 4) { // > because we accept -N(=O)=O as a valid nitro group 1975 oxygenCount = nitrogenCount = doubleBondTo = 0; 1976 1977 FOR_NBORS_OF_ATOM (nbr, atom) { 1978 if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) { 1979 oxygenCount++; 1980 } 1981 if (nbr->GetAtomicNum() == OBElements::Nitrogen) { 1982 bond = _mol.GetBond(&*nbr, atom); 1983 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 1984 doubleBondTo = 7; 1985 } 1986 } 1987 if (nbr->GetAtomicNum() == OBElements::Carbon) { 1988 bond = _mol.GetBond(&*nbr, atom); 1989 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 1990 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 1991 if (nbrNbr->GetAtomicNum() == OBElements::Nitrogen && (nbrNbr->GetExplicitDegree() == 3)) { 1992 nitrogenCount++; 1993 } 1994 } 1995 } 1996 } 1997 } 1998 1999 if (oxygenCount == 1) { 2000 return 67; // sp2-hybridized N-oxide nitrogen (N2OX) 2001 } 2002 if (oxygenCount >= 2) { 2003 return 45; // Nitrogen in nitro group, Nitrogen in nitrate group (NO2, NO3) 2004 } 2005 2006 if (nitrogenCount == 1) { 2007 return 54; // Iminium nitrogen (N+=C) 2008 } 2009 if (nitrogenCount == 2) { 2010 return 55; // Either nitrogen in N+=C-N: (NCN+) 2011 } 2012 if (nitrogenCount == 3) { 2013 return 56; // Guanidinium nitrogen (NGD+) 2014 } 2015 2016 if (doubleBondTo == 7) { 2017 return 54; // Positivly charged nitrogen doubly bonded to nitrogen (N+=N) 2018 } 2019 } 2020 2021 if (atom->GetExplicitValence() == 3) { 2022 bool IsAmide = false; 2023 bool IsSulfonAmide = false; 2024 bool IsNNNorNNC = false; 2025 int tripleBondTo = 0; 2026 doubleBondTo = 0; 2027 2028 FOR_NBORS_OF_ATOM (nbr, atom) { 2029 if (nbr->GetAtomicNum() == OBElements::Sulfur || nbr->GetAtomicNum() == OBElements::Phosphorus) { 2030 oxygenCount = 0; 2031 2032 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 2033 if (nbrNbr->GetAtomicNum() == OBElements::Oxygen && (nbrNbr->GetExplicitDegree() == 1)) { 2034 oxygenCount++; 2035 } 2036 } 2037 if (oxygenCount >= 2) { 2038 IsSulfonAmide = true; 2039 //return 43; // Sulfonamide nitrogen (NSO2, NSO3) 2040 } 2041 } 2042 } 2043 2044 FOR_NBORS_OF_ATOM (nbr, atom) { 2045 if (nbr->GetAtomicNum() == OBElements::Carbon) { 2046 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 2047 bond = _mol.GetBond(&*nbr, &*nbrNbr); 2048 if (!bond->IsAromatic() && bond->GetBondOrder() == 2 && (nbrNbr->GetAtomicNum() == OBElements::Oxygen || nbrNbr->GetAtomicNum() == OBElements::Sulfur)) { 2049 IsAmide = true; 2050 //return 10; // Amide nitrogen, Thioamide nitrogen (NC=O, NC=S) 2051 } 2052 } 2053 } 2054 } 2055 2056 FOR_NBORS_OF_ATOM (nbr, atom) { 2057 if (nbr->GetAtomicNum() == OBElements::Carbon) { 2058 int N2count = 0; 2059 int N3count = 0; 2060 oxygenCount = sulphurCount = 0; 2061 2062 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 2063 bond = _mol.GetBond(&*nbr, &*nbrNbr); 2064 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 2065 doubleBondTo = nbrNbr->GetAtomicNum(); 2066 } 2067 if (bond->IsAromatic()) { 2068 if ((nbrNbr->GetAtomicNum() == 7) || (nbrNbr->GetAtomicNum() == 6)) { 2069 doubleBondTo = nbrNbr->GetAtomicNum(); 2070 } 2071 } 2072 if (!bond->IsAromatic() && bond->GetBondOrder() == 3) { 2073 tripleBondTo = nbrNbr->GetAtomicNum(); 2074 } 2075 if (nbrNbr->GetAtomicNum() == OBElements::Nitrogen && (nbrNbr->GetExplicitDegree() == 3)) { 2076 int nbrOxygen = 0; 2077 FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) { 2078 if (nbrNbrNbr->GetAtomicNum() == OBElements::Oxygen) { 2079 nbrOxygen++; 2080 } 2081 } 2082 if (nbrOxygen < 2) { 2083 N3count++; 2084 } 2085 } 2086 if (nbrNbr->GetAtomicNum() == OBElements::Nitrogen && (nbrNbr->GetExplicitDegree() == 2) && (bond->GetBondOrder() == 2 || bond->IsAromatic())) { 2087 N2count++; 2088 } 2089 if (nbrNbr->IsAromatic()) { 2090 if (nbrNbr->GetAtomicNum() == OBElements::Oxygen) { 2091 oxygenCount++; 2092 } 2093 if (nbrNbr->GetAtomicNum() == OBElements::Sulfur) { 2094 sulphurCount++; 2095 } 2096 } 2097 } 2098 if (N3count == 3) { 2099 return 56; // Guanidinium nitrogen (NGD+) 2100 } 2101 2102 if (!IsAmide && !IsSulfonAmide && !oxygenCount && !sulphurCount && nbr->IsAromatic()) { 2103 return 40; 2104 } 2105 2106 if ((N3count == 2) && (doubleBondTo == 7) && !N2count) { 2107 return 55; // Either nitrogen in N+=C-N: (NCN+) 2108 } 2109 } 2110 2111 if (nbr->GetAtomicNum() == OBElements::Nitrogen) { 2112 nitrogenCount = 0; 2113 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 2114 bond = _mol.GetBond(&*nbr, &*nbrNbr); 2115 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 2116 if (nbrNbr->GetAtomicNum() == OBElements::Carbon) { 2117 oxygenCount = sulphurCount = 0; 2118 FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) { 2119 if (nbrNbrNbr->GetAtomicNum() == OBElements::Oxygen) { 2120 oxygenCount++; 2121 } 2122 if (nbrNbrNbr->GetAtomicNum() == OBElements::Sulfur) { 2123 sulphurCount++; 2124 } 2125 if (nbrNbrNbr->GetAtomicNum() == OBElements::Sulfur) { 2126 nitrogenCount++; 2127 } 2128 } 2129 if (!oxygenCount && !sulphurCount && (nitrogenCount == 1)) { 2130 bool bondToAromC = false; 2131 FOR_NBORS_OF_ATOM (nbr2, atom) { 2132 if (nbr2->IsAromatic() && nbr2->GetAtomicNum() == OBElements::Carbon && nbr2->IsInRingSize(6)) { 2133 bondToAromC = true; 2134 } 2135 } 2136 if (!bondToAromC) { 2137 IsNNNorNNC = true; 2138 } 2139 } 2140 } 2141 if (nbrNbr->GetAtomicNum() == OBElements::Nitrogen) { 2142 bool bondToAromC = false; 2143 FOR_NBORS_OF_ATOM (nbr2, atom) { 2144 if (nbr2->IsAromatic() && nbr2->GetAtomicNum() == OBElements::Carbon && nbr2->IsInRingSize(6)) { 2145 bondToAromC = true; 2146 } 2147 } 2148 if (!bondToAromC) { 2149 IsNNNorNNC = true; 2150 } 2151 } 2152 } 2153 } 2154 } 2155 } 2156 2157 if (IsSulfonAmide) { 2158 return 43; // Sulfonamide nitrogen (NSO2, NSO3) 2159 } 2160 if (IsAmide) { 2161 return 10; // Amide nitrogen, Thioamide nitrogen (NC=O, NC=S) 2162 } 2163 2164 if ((doubleBondTo == 6) || (doubleBondTo == 7) ||(doubleBondTo == 15) || (tripleBondTo == 6)) { 2165 return 40; // Enamine or aniline nitrogen (deloc. lp), Nitrogen in N-C=N with deloc. lp, 2166 // Nitrogen in N-C=N with deloc. lp, Nitrogen attached to C-C triple bond 2167 // (NC=C, NC=N, NC=P, NC%C) 2168 } 2169 if (tripleBondTo == 7) { 2170 return 43; // Nitrogen attached to cyano group (NC%N) 2171 } 2172 if (IsNNNorNNC) { 2173 return 10; // Nitrogen in N-N=C moiety with deloc. lp 2174 // Nitrogen in N-N=N moiety with deloc. lp (NN=C, NN=N) 2175 } 2176 2177 return 8; // Amine nitrogen (NR) 2178 } 2179 } 2180 // 2 neighbours 2181 if (atom->GetExplicitDegree() == 2) { 2182 if (atom->GetExplicitValence() >= 4) { 2183 bool isNbrCarbon = false; 2184 bool isBondTriple = false; 2185 FOR_NBORS_OF_ATOM (nbr, atom) { 2186 if (!isNbrCarbon) 2187 isNbrCarbon = nbr->GetAtomicNum() == OBElements::Carbon; 2188 bond = _mol.GetBond(&*nbr, atom); 2189 if (!isBondTriple) 2190 isBondTriple = !bond->IsAromatic() && bond->GetBondOrder() == 3; 2191 } 2192 if (isBondTriple && isNbrCarbon) 2193 return 61; // Isonitrile nitrogen (NR%) 2194 2195 return 53; // Central nitrogen in C=N=N or N=N=N (=N=) 2196 } 2197 2198 if (atom->GetExplicitValence() == 3) { 2199 doubleBondTo = 0; 2200 2201 FOR_NBORS_OF_ATOM (nbr, atom) { 2202 bond = _mol.GetBond(&*nbr, atom); 2203 if (nbr->GetAtomicNum() == OBElements::Oxygen && !bond->IsAromatic() && bond->GetBondOrder() == 2 && (nbr->GetExplicitDegree() == 1)) { 2204 return 46; // Nitrogen in nitroso group (N=O) 2205 } 2206 if ((nbr->GetAtomicNum() == OBElements::Carbon || nbr->GetAtomicNum() == OBElements::Nitrogen) && !bond->IsAromatic() && bond->GetBondOrder() == 2) { 2207 return 9; // Iminie nitrogen, Azo-group nitrogen (N=C, N=N) 2208 } 2209 } 2210 FOR_NBORS_OF_ATOM (nbr, atom) { 2211 if (nbr->GetAtomicNum() == OBElements::Sulfur) { 2212 oxygenCount = 0; 2213 2214 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 2215 if (nbrNbr->GetAtomicNum() == OBElements::Oxygen && (nbrNbr->GetExplicitDegree() == 1)) { 2216 oxygenCount++; 2217 } 2218 } 2219 if (oxygenCount >= 2) { 2220 return 43; // Sulfonamide nitrogen (NSO2, NSO3) 2221 } 2222 } 2223 } 2224 } 2225 2226 if (atom->GetExplicitValence() >= 2) { // Bug reported by Paolo Tosco 2227 oxygenCount = sulphurCount = 0; 2228 2229 FOR_NBORS_OF_ATOM (nbr, atom) { 2230 if (nbr->GetAtomicNum() == OBElements::Sulfur) { 2231 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 2232 if (nbrNbr->GetAtomicNum() == OBElements::Oxygen && (nbrNbr->GetExplicitDegree() == 1)) { 2233 oxygenCount++; 2234 } 2235 } 2236 if (oxygenCount == 1) { 2237 return 48; // Divalent nitrogen replacing monovalent O in SO2 group (NSO) 2238 } 2239 } 2240 } 2241 2242 return 62; // Anionic divalent nitrogen (NM) 2243 } 2244 } 2245 // 1 neighbours 2246 if (atom->GetExplicitDegree() == 1) { 2247 FOR_NBORS_OF_ATOM (nbr, atom) { 2248 bond = _mol.GetBond(&*nbr, atom); 2249 if (!bond->IsAromatic() && bond->GetBondOrder() == 3) { 2250 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 2251 if (atom != &*nbrNbr && !(nbr->GetAtomicNum() == OBElements::Nitrogen 2252 && nbrNbr->GetAtomicNum() == OBElements::Nitrogen && nbrNbr->GetExplicitDegree() == 2)) { 2253 return 42; // Triply bonded nitrogen (NSP) 2254 } 2255 } 2256 } 2257 if (nbr->GetAtomicNum() == OBElements::Nitrogen && (nbr->GetExplicitDegree() == 2)) { 2258 return 47; // Terminal nitrogen in azido group (NAZT) 2259 } 2260 } 2261 } 2262 return 8; // generic amine nitrogen 2263 } 2264 2265 //////////////////////////////// 2266 // Oxygen 2267 //////////////////////////////// 2268 if (atom->GetAtomicNum() == 8) { 2269 // 3 neighbours 2270 if (atom->GetExplicitDegree() == 3) { 2271 return 49; // Oxonium oxygen (O+) 2272 } 2273 // 2 neighbours 2274 if (atom->GetExplicitDegree() == 2) { 2275 int hydrogenCount = 0; 2276 FOR_NBORS_OF_ATOM (nbr, atom) { 2277 if (nbr->GetAtomicNum() == OBElements::Hydrogen) { 2278 hydrogenCount++; 2279 } 2280 } 2281 2282 if (hydrogenCount == 2) { 2283 // H--O--H 2284 return 70; // Oxygen in water (OH2) 2285 } 2286 if (atom->GetExplicitValence() == 3) { 2287 return 51; // Oxenium oxygen (O=+) 2288 } 2289 2290 return 6; // Generic divalent oxygen, Ether oxygen, Carboxylic acid or ester oxygen, 2291 // Enolic or phenolic oxygen, Oxygen in -O-C=N- moiety, Divalent oxygen in 2292 // thioacid or ester, Divalent nitrate "ether" oxygen, Divalent oxygen in 2293 // sulfate group, Divalent oxygen in sulfite group, One of two divalent 2294 // oxygens attached to sulfur, Divalent oxygen in R(RO)S=O, Other divalent 2295 // oxygen attached to sulfur, Divalent oxygen in phosphate group, Divalent 2296 // oxygen in phosphite group, Divalent oxygen (one of two oxygens attached 2297 // to P), Other divalent oxygen (-O-, OR, OC=O, OC=C, OC=N, OC=S, ONO2, 2298 // ON=O, OSO3, OSO2, OSO, OS=O, -OS, OPO3, OPO2, OPO, -OP) 2299 2300 // 59 ar 2301 } 2302 // 1 neighbour 2303 if (atom->GetExplicitDegree() == 1) { 2304 oxygenCount = sulphurCount = 0; 2305 2306 FOR_NBORS_OF_ATOM (nbr, atom) { 2307 bond = _mol.GetBond(&*nbr, atom); 2308 2309 if (nbr->GetAtomicNum() == OBElements::Carbon || nbr->GetAtomicNum() == OBElements::Nitrogen) { 2310 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 2311 if (nbrNbr->GetAtomicNum() == OBElements::Oxygen && (nbrNbr->GetExplicitDegree() == 1)) { 2312 oxygenCount++; 2313 } 2314 if (nbrNbr->GetAtomicNum() == OBElements::Sulfur && (nbrNbr->GetExplicitDegree() == 1)) { 2315 sulphurCount++; 2316 } 2317 } 2318 } 2319 // O---H 2320 if (nbr->GetAtomicNum() == OBElements::Hydrogen) { 2321 return 35; 2322 } 2323 // O-?-C 2324 if (nbr->GetAtomicNum() == OBElements::Carbon) { 2325 if (oxygenCount == 2) { 2326 // O-?-C-?-O 2327 return 32; // Oxygen in carboxylate group (O2CM) 2328 } 2329 if (!bond->IsAromatic() && bond->GetBondOrder() == 1) { 2330 // O--C 2331 return 35; // Oxide oxygen on sp3 carbon, Oxide oxygen on sp2 carbon (OM, OM2) 2332 } else { 2333 // O==C 2334 return 7; // Generic carbonyl oxygen, Carbonyl oxygen in amides, 2335 // Carbonyl oxygen in aldehydes and ketones, Carbonyl 2336 // oxygen in acids or esters (O=C, O=CN, O=CR, O=CO) 2337 } 2338 } 2339 // O-?-N 2340 if (nbr->GetAtomicNum() == OBElements::Nitrogen) { 2341 if (oxygenCount >= 2) { 2342 // O-?-N-?-O 2343 return 32; // Oxygen in nitro group, Nitro-group oxygen in nitrate, 2344 // Nitrate anion oxygen (O2N, O2NO, O3N) 2345 } 2346 if (!bond->IsAromatic() && bond->GetBondOrder() == 1) { 2347 if ((nbr->GetExplicitDegree() == 2) || (nbr->GetExplicitValence() == 3)) 2348 // O(-)--N 2349 return 35; 2350 else 2351 // O--N 2352 return 32; // Oxygen in N-oxides (ONX) 2353 } else { 2354 // sometimes ONX bonds are labelled as double bonds 2355 // (e.g. by MOE, noticed by Paolo Tosco) 2356 int ndab = 0; 2357 FOR_BONDS_OF_ATOM(bond2, &*nbr) { 2358 if (bond2->GetBondOrder() == 2 || bond2->GetBondOrder() == 5) 2359 ndab++; 2360 } 2361 if (ndab + nbr->GetExplicitDegree() == 5) 2362 // O--N 2363 return 32; // Oxygen in N-oxides (ONX) 2364 else 2365 // O==N 2366 return 7; // Nitroso oxygen (O=N) 2367 } 2368 } 2369 // O-?-S 2370 if (nbr->GetAtomicNum() == OBElements::Sulfur) { 2371 if (sulphurCount == 1) { 2372 // O1-?-S-?-S1 2373 return 32; // Terminal oxygen in thiosulfinate anion (OSMS) 2374 } 2375 if (!bond->IsAromatic() && bond->GetBondOrder() == 1) { 2376 // O--S 2377 return 32; // Single terminal oxygen on sulfur, One of 2 terminal O's on sulfur, 2378 // One of 3 terminal O's on sulfur, Terminal O in sulfate anion, 2379 // (O-S, O2S, O3S, O4S) 2380 } else { 2381 // O==S 2382 2383 // are all sulfur nbr atoms carbon? 2384 bool isSulfoxide = true; 2385 int oxygenBoundToSulfur = 0; 2386 FOR_NBORS_OF_ATOM (nbr2, &*nbr) { 2387 if (atom == &*nbr2) 2388 continue; 2389 2390 if (nbr2->GetAtomicNum() == OBElements::Oxygen) 2391 ++oxygenBoundToSulfur; 2392 } 2393 FOR_NBORS_OF_ATOM (nbr2, &*nbr) { 2394 if (atom == &*nbr2) 2395 continue; 2396 OBBond* bond = nbr->GetBond(&*nbr2); 2397 if (!bond->IsAromatic() && bond->GetBondOrder() == 2 2398 && nbr2->GetAtomicNum() == OBElements::Carbon && oxygenBoundToSulfur == 1) 2399 isSulfoxide = false; // O=S on sulfur doubly bonded to, e.g., C (O=S=) 2400 2401 if ((nbr2->GetAtomicNum() == OBElements::Oxygen && nbr2->GetExplicitDegree() == 1) 2402 || (nbr2->GetAtomicNum() == OBElements::Nitrogen && nbr2->GetExplicitDegree() == 2)) 2403 isSulfoxide = false; 2404 } 2405 2406 if (isSulfoxide) 2407 return 7; // Doubly bonded sulfoxide oxygen (O=S) 2408 else 2409 return 32; // (O2S, O2S=C, O3S, O4S) 2410 } 2411 } 2412 2413 return 32; // Oxygen in phosphine oxide, One of 2 terminal O's on sulfur, 2414 // One of 3 terminal O's on sulfur, One of 4 terminal O's on sulfur, 2415 // Oxygen in perchlorate anion (OP, O2P, O3P, O4P, O4Cl) 2416 } 2417 } 2418 } 2419 2420 //////////////////////////////// 2421 // Flourine 2422 //////////////////////////////// 2423 if (atom->GetAtomicNum() == 9) { 2424 // 1 neighbour 2425 if (atom->GetExplicitDegree() == 1) { 2426 return 11; // Fluorine (F) 2427 } 2428 // 0 neighbours 2429 if (atom->GetExplicitDegree() == 0) { 2430 return 89; // Fluoride anion (F-) 2431 } 2432 } 2433 2434 //////////////////////////////// 2435 // Sodium 2436 //////////////////////////////// 2437 if (atom->GetAtomicNum() == 11) { 2438 return 93; // Sodium cation (NA+) 2439 } 2440 2441 //////////////////////////////// 2442 // Magnesium 2443 //////////////////////////////// 2444 if (atom->GetAtomicNum() == 12) { 2445 return 99; // Dipositive magnesium cation (MG+2) 2446 } 2447 2448 //////////////////////////////// 2449 // Silicon 2450 //////////////////////////////// 2451 if (atom->GetAtomicNum() == 14) { 2452 return 19; // Silicon (SI) 2453 } 2454 2455 //////////////////////////////// 2456 // Phosphorus 2457 //////////////////////////////// 2458 if (atom->GetAtomicNum() == 15) { 2459 if (atom->GetExplicitDegree() == 4) { 2460 return 25; // Phosphate group phosphorus, Phosphorus with 3 attached oxygens, 2461 // Phosphorus with 2 attached oxygens, Phosphine oxide phosphorus, 2462 // General tetracoordinate phosphorus (PO4, PO3, PO2, PO, PTET) 2463 } 2464 if (atom->GetExplicitDegree() == 3) { 2465 return 26; // Phosphorus in phosphines (P) 2466 } 2467 if (atom->GetExplicitDegree() == 2) { 2468 return 75; // Phosphorus doubly bonded to C (-P=C) 2469 } 2470 } 2471 2472 //////////////////////////////// 2473 // Sulfur 2474 //////////////////////////////// 2475 if (atom->GetAtomicNum() == 16) { 2476 // 4 neighbours 2477 if (atom->GetExplicitDegree() == 4) { 2478 return 18; // Sulfone sulfur, Sulfonamide sulfur, Sulfonate group sulfur, 2479 // Sulfate group sulfur, Sulfur in nitrogen analog of sulfone 2480 // (SO2, SO2N, SO3, SO4, SNO) 2481 } 2482 // 3 neighbours 2483 if (atom->GetExplicitDegree() == 3) { 2484 oxygenCount = sulphurCount = doubleBondTo = 0; 2485 2486 FOR_NBORS_OF_ATOM (nbr, atom) { 2487 bond = _mol.GetBond(&*nbr, atom); 2488 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 2489 if (nbr->GetAtomicNum() == 6) 2490 doubleBondTo = 6; 2491 } 2492 2493 if (nbr->GetExplicitDegree() == 1) { 2494 if (nbr->GetAtomicNum() == OBElements::Oxygen) { 2495 oxygenCount++; 2496 } else if (nbr->GetAtomicNum() == OBElements::Sulfur) { 2497 sulphurCount++; 2498 } 2499 } 2500 } 2501 2502 if (oxygenCount == 2) { 2503 if (doubleBondTo == 6) { 2504 return 18; // Sulfone sulfur, doubly bonded to carbon (=SO2) 2505 } 2506 return 73; // Sulfur in anionic sulfinate group (SO2M) 2507 } 2508 if (oxygenCount && sulphurCount) 2509 return 73; // Tricoordinate sulfur in anionic thiosulfinate group (SSOM) 2510 2511 //if ((doubleBondTo == 6) || (doubleBondTo == 8)) 2512 return 17; // Sulfur doubly bonded to carbon, Sulfoxide sulfur (S=C, S=O) 2513 } 2514 // 2 neighbours 2515 if (atom->GetExplicitDegree() == 2) { 2516 doubleBondTo = 0; 2517 2518 FOR_NBORS_OF_ATOM (nbr, atom) { 2519 if (nbr->GetAtomicNum() == OBElements::Oxygen) { 2520 bond = _mol.GetBond(&*nbr, atom); 2521 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 2522 doubleBondTo = 8; 2523 } 2524 } 2525 } 2526 2527 if (doubleBondTo == 8) 2528 return 74; // Sulfinyl sulfur, e.g., in C=S=O (=S=O) 2529 2530 return 15; // Thiol, sulfide, or disulfide sulfor (S) 2531 } 2532 // 1 neighbour 2533 if (atom->GetExplicitDegree() == 1) { 2534 sulphurCount = doubleBondTo = 0; 2535 2536 FOR_NBORS_OF_ATOM (nbr, atom) { 2537 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) { 2538 if (nbrNbr->GetAtomicNum() == OBElements::Sulfur && (nbrNbr->GetExplicitDegree() == 1)) { 2539 sulphurCount++; 2540 } 2541 } 2542 bond = _mol.GetBond(&*nbr, atom); 2543 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) { 2544 doubleBondTo = nbr->GetAtomicNum(); 2545 } 2546 } 2547 2548 if ((doubleBondTo == 6) && (sulphurCount != 2)) { 2549 return 16; // Sulfur doubly bonded to carbon (S=C) 2550 } 2551 2552 return 72; // Terminal sulfur bonded to P, Anionic terminal sulfur, 2553 // Terminal sulfur in thiosulfinate group (S-P, SM, SSMO) 2554 } 2555 2556 // 44 ar 2557 } 2558 2559 //////////////////////////////// 2560 // Clorine 2561 //////////////////////////////// 2562 if (atom->GetAtomicNum() == 17) { 2563 // 4 neighbour 2564 if (atom->GetExplicitDegree() == 4) { 2565 oxygenCount = 0; 2566 2567 FOR_NBORS_OF_ATOM (nbr, atom) { 2568 if (nbr->GetAtomicNum() == OBElements::Oxygen) { 2569 oxygenCount++; 2570 } 2571 } 2572 if (oxygenCount == 4) 2573 return 77; // Perchlorate anion chlorine (CLO4) 2574 } 2575 // 1 neighbour 2576 if (atom->GetExplicitDegree() == 1) { 2577 return 12; // Chlorine (CL) 2578 } 2579 // 0 neighbours 2580 if (atom->GetExplicitDegree() == 0) { 2581 return 90; // Chloride anion (CL-) 2582 } 2583 } 2584 2585 //////////////////////////////// 2586 // Potasium 2587 //////////////////////////////// 2588 if (atom->GetAtomicNum() == 19) { 2589 return 94; // Potasium cation (K+) 2590 } 2591 2592 //////////////////////////////// 2593 // Calcium 2594 //////////////////////////////// 2595 if (atom->GetAtomicNum() == 20) { 2596 // 0 neighbours 2597 if (atom->GetExplicitDegree() == 0) { 2598 return 96; // Dipositive calcium cation (CA+2) 2599 } 2600 } 2601 2602 //////////////////////////////// 2603 // Iron 2604 //////////////////////////////// 2605 if (atom->GetAtomicNum() == 26) { 2606 if (atom->GetFormalCharge() == 2) 2607 return 87; // Dipositive iron (FE+2) 2608 else 2609 return 88; // Tripositive iron (FE+3) 2610 } 2611 2612 //////////////////////////////// 2613 // Copper 2614 //////////////////////////////// 2615 if (atom->GetAtomicNum() == 29) { 2616 if (atom->GetFormalCharge() == 1) 2617 return 97; // Monopositive copper cation (CU+1) 2618 else 2619 return 98; // Dipositive copper cation (CU+2) 2620 } 2621 2622 //////////////////////////////// 2623 // Zinc 2624 //////////////////////////////// 2625 if (atom->GetAtomicNum() == 30) { 2626 return 95; // Dipositive zinc cation (ZN+2) 2627 } 2628 2629 //////////////////////////////// 2630 // Bromine 2631 //////////////////////////////// 2632 if (atom->GetAtomicNum() == 35) { 2633 // 1 neighbour 2634 if (atom->GetExplicitDegree() == 1) { 2635 return 13; // Bromine (BR) 2636 } 2637 // 0 neighbours 2638 if (atom->GetExplicitDegree() == 0) { 2639 return 91; // Bromide anion (BR-) 2640 } 2641 } 2642 2643 //////////////////////////////// 2644 // Iodine 2645 //////////////////////////////// 2646 if (atom->GetAtomicNum() == 53) { 2647 // 1 neighbour 2648 if (atom->GetExplicitDegree() == 1) { 2649 return 14; // Iodine (I) 2650 } 2651 } 2652 2653 2654 2655 return 0; 2656 } 2657 SetTypes()2658 bool OBForceFieldMMFF94::SetTypes() 2659 { 2660 char type[4]; 2661 2662 _mol.SetAtomTypesPerceived(); 2663 2664 // mark all atoms and bonds as non-aromatic 2665 _mol.SetAromaticPerceived(); 2666 FOR_BONDS_OF_MOL (bond, _mol) 2667 bond->SetAromatic(false); 2668 FOR_ATOMS_OF_MOL (atom, _mol) 2669 atom->SetAromatic(false); 2670 2671 // It might be needed to run this function more than once... 2672 bool done = true; 2673 while (done) { 2674 done = PerceiveAromatic(); 2675 } 2676 2677 FOR_ATOMS_OF_MOL (atom, _mol) { 2678 snprintf(type, 3, "%d", GetType(&*atom)); 2679 atom->SetType(type); 2680 } 2681 2682 PrintTypes(); 2683 2684 return true; 2685 } 2686 SetupCalculations()2687 bool OBForceFieldMMFF94::SetupCalculations() 2688 { 2689 OBFFParameter *parameter; 2690 OBAtom *a, *b, *c, *d; 2691 int type_a, type_b, type_c, type_d; 2692 bool found; 2693 int order; 2694 2695 IF_OBFF_LOGLVL_LOW 2696 OBFFLog("\nS E T T I N G U P C A L C U L A T I O N S\n\n"); 2697 2698 // 2699 // Bond Calculations 2700 // 2701 // no "step-down" procedure 2702 // MMFF part V - page 625 (empirical rule) 2703 // 2704 IF_OBFF_LOGLVL_LOW 2705 OBFFLog("SETTING UP BOND CALCULATIONS...\n"); 2706 2707 OBFFBondCalculationMMFF94 bondcalc; 2708 int bondtype; 2709 2710 _bondcalculations.clear(); 2711 2712 FOR_BONDS_OF_MOL(bond, _mol) { 2713 a = bond->GetBeginAtom(); 2714 b = bond->GetEndAtom(); 2715 2716 // skip this bond if the atoms are ignored 2717 if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) ) 2718 continue; 2719 2720 // if there are any groups specified, check if the two bond atoms are in a single intraGroup 2721 if (HasGroups()) { 2722 bool validBond = false; 2723 for (unsigned int i=0; i < _intraGroup.size(); ++i) { 2724 if (_intraGroup[i].BitIsSet(a->GetIdx()) && _intraGroup[i].BitIsSet(b->GetIdx())) { 2725 validBond = true; 2726 break; 2727 } 2728 } 2729 if (!validBond) 2730 continue; 2731 } 2732 2733 bondtype = GetBondType(a, b); 2734 2735 parameter = GetTypedParameter2Atom(bondtype, atoi(a->GetType()), atoi(b->GetType()), _ffbondparams); // from mmffbond.par 2736 if (parameter == nullptr) { 2737 parameter = GetParameter2Atom(a->GetAtomicNum(), b->GetAtomicNum(), _ffbndkparams); // from mmffbndk.par - emperical rules 2738 if (parameter == nullptr) { 2739 IF_OBFF_LOGLVL_LOW { 2740 // This should never happen 2741 snprintf(_logbuf, BUFF_SIZE, " COULD NOT FIND PARAMETERS FOR BOND %d-%d (IDX)...\n", a->GetIdx(), b->GetIdx()); 2742 OBFFLog(_logbuf); 2743 } 2744 return false; 2745 } else { 2746 IF_OBFF_LOGLVL_LOW { 2747 snprintf(_logbuf, BUFF_SIZE, " USING EMPIRICAL RULE FOR BOND STRETCHING %d-%d (IDX)...\n", a->GetIdx(), b->GetIdx()); 2748 OBFFLog(_logbuf); 2749 } 2750 2751 double rr, rr2, rr4, rr6; 2752 bondcalc.a = a; 2753 bondcalc.b = b; 2754 bondcalc.r0 = GetRuleBondLength(a, b); 2755 2756 rr = parameter->_dpar[0] / bondcalc.r0; // parameter->_dpar[0] = r0-ref 2757 rr2 = rr * rr; 2758 rr4 = rr2 * rr2; 2759 rr6 = rr4 * rr2; 2760 2761 bondcalc.kb = parameter->_dpar[1] * rr6; // parameter->_dpar[1] = kb-ref 2762 bondcalc.bt = bondtype; 2763 bondcalc.SetupPointers(); 2764 2765 _bondcalculations.push_back(bondcalc); 2766 } 2767 } else { 2768 bondcalc.a = a; 2769 bondcalc.b = b; 2770 bondcalc.kb = parameter->_dpar[0]; 2771 bondcalc.r0 = parameter->_dpar[1]; 2772 bondcalc.bt = bondtype; 2773 bondcalc.SetupPointers(); 2774 2775 _bondcalculations.push_back(bondcalc); 2776 } 2777 } 2778 2779 // 2780 // Angle Calculations 2781 // 2782 // MMFF part I - page 513 ("step-down" prodedure) 2783 // MMFF part I - page 519 (reference 68 is actually a footnote) 2784 // MMFF part V - page 627 (empirical rule) 2785 // 2786 // First try and find an exact match, if this fails, step down using the equivalences from mmffdef.par 2787 // five-stage protocol: 1-1-1, 2-2-2, 3-2-3, 4-2-4, 5-2-5 2788 // If this fails, use empirical rules 2789 // Since 1-1-1 = 2-2-2, we will only try 1-1-1 before going to 3-2-3 2790 // 2791 // Stretch-Bend Calculations 2792 // 2793 IF_OBFF_LOGLVL_LOW 2794 OBFFLog("SETTING UP ANGLE & STRETCH-BEND CALCULATIONS...\n"); 2795 2796 OBFFAngleCalculationMMFF94 anglecalc; 2797 OBFFStrBndCalculationMMFF94 strbndcalc; 2798 int angletype, strbndtype, bondtype1, bondtype2; 2799 2800 _anglecalculations.clear(); 2801 _strbndcalculations.clear(); 2802 2803 FOR_ANGLES_OF_MOL(angle, _mol) { 2804 b = _mol.GetAtom((*angle)[0] + 1); 2805 a = _mol.GetAtom((*angle)[1] + 1); 2806 c = _mol.GetAtom((*angle)[2] + 1); 2807 2808 type_a = atoi(a->GetType()); 2809 type_b = atoi(b->GetType()); 2810 type_c = atoi(c->GetType()); 2811 2812 // skip this angle if the atoms are ignored 2813 if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) || _constraints.IsIgnored(c->GetIdx()) ) 2814 continue; 2815 2816 // if there are any groups specified, check if the three angle atoms are in a single intraGroup 2817 if (HasGroups()) { 2818 bool validAngle = false; 2819 for (unsigned int i=0; i < _intraGroup.size(); ++i) { 2820 if (_intraGroup[i].BitIsSet(a->GetIdx()) && _intraGroup[i].BitIsSet(b->GetIdx()) && 2821 _intraGroup[i].BitIsSet(c->GetIdx())) { 2822 validAngle = true; 2823 break; 2824 } 2825 } 2826 if (!validAngle) 2827 continue; 2828 } 2829 2830 angletype = GetAngleType(a, b, c); 2831 strbndtype = GetStrBndType(a, b, c); 2832 bondtype1 = GetBondType(a, b); 2833 bondtype2 = GetBondType(b, c); 2834 2835 if (HasLinSet(type_b)) { 2836 anglecalc.linear = true; 2837 } else { 2838 anglecalc.linear = false; 2839 } 2840 2841 // try exact match 2842 parameter = GetTypedParameter3Atom(angletype, type_a, type_b, type_c, _ffangleparams); 2843 if (parameter == nullptr) // try 3-2-3 2844 parameter = GetTypedParameter3Atom(angletype, EqLvl3(type_a), type_b, EqLvl3(type_c), _ffangleparams); 2845 if (parameter == nullptr) // try 4-2-4 2846 parameter = GetTypedParameter3Atom(angletype, EqLvl4(type_a), type_b, EqLvl4(type_c), _ffangleparams); 2847 if (parameter == nullptr) // try 5-2-5 2848 parameter = GetTypedParameter3Atom(angletype, EqLvl5(type_a), type_b, EqLvl5(type_c), _ffangleparams); 2849 2850 if (parameter) { 2851 anglecalc.ka = parameter->_dpar[0]; 2852 anglecalc.theta0 = parameter->_dpar[1]; 2853 strbndcalc.theta0 = parameter->_dpar[1]; // ** 2854 } else { 2855 IF_OBFF_LOGLVL_LOW { 2856 snprintf(_logbuf, BUFF_SIZE, " USING DEFAULT ANGLE FOR %d-%d-%d (IDX)...\n", a->GetIdx(), b->GetIdx(), c->GetIdx()); 2857 snprintf(_logbuf, BUFF_SIZE, " USING EMPIRICAL RULE FOR ANGLE BENDING %d-%d-%d (IDX)...\n", a->GetIdx(), b->GetIdx(), c->GetIdx()); 2858 OBFFLog(_logbuf); 2859 } 2860 2861 anglecalc.ka = 0.0; 2862 anglecalc.theta0 = 120.0; 2863 2864 if (GetCrd(type_b) == 4) 2865 anglecalc.theta0 = 109.45; 2866 2867 if ((GetCrd(type_b) == 2) && b->GetAtomicNum() == OBElements::Oxygen) 2868 anglecalc.theta0 = 105.0; 2869 2870 if (b->GetAtomicNum() > 10) 2871 anglecalc.theta0 = 95.0; 2872 2873 if (HasLinSet(type_b)) 2874 anglecalc.theta0 = 180.0; 2875 2876 if ((GetCrd(type_b) == 3) && (GetVal(type_b) == 3) && !GetMltb(type_b)) { 2877 if (b->GetAtomicNum() == OBElements::Nitrogen) { 2878 anglecalc.theta0 = 107.0; 2879 } else { 2880 anglecalc.theta0 = 92.0; 2881 } 2882 } 2883 2884 if (a->IsInRingSize(3) && b->IsInRingSize(3) && c->IsInRingSize(3) && IsInSameRing(a, c)) 2885 anglecalc.theta0 = 60.0; 2886 2887 if (a->IsInRingSize(4) && b->IsInRingSize(4) && c->IsInRingSize(4) && IsInSameRing(a, c)) 2888 anglecalc.theta0 = 90.0; 2889 2890 strbndcalc.theta0 = anglecalc.theta0; // ** 2891 } 2892 2893 // empirical rule for 0-b-0 and standard angles 2894 if (anglecalc.ka == 0.0) { 2895 IF_OBFF_LOGLVL_LOW { 2896 snprintf(_logbuf, BUFF_SIZE, " USING EMPIRICAL RULE FOR ANGLE BENDING FORCE CONSTANT %d-%d-%d (IDX)...\n", a->GetIdx(), b->GetIdx(), c->GetIdx()); 2897 OBFFLog(_logbuf); 2898 } 2899 2900 double beta, Za, Zc, Cb, r0ab, r0bc, theta, theta2, D, rr, rr2; 2901 Za = GetZParam(a); 2902 Cb = GetCParam(b); // Fixed typo -- PR#2741658 2903 Zc = GetZParam(c); 2904 2905 r0ab = GetBondLength(a, b); 2906 r0bc = GetBondLength(b, c); 2907 rr = r0ab + r0bc; 2908 rr2 = rr * rr; 2909 D = (r0ab - r0bc) / rr2; 2910 2911 theta = anglecalc.theta0; 2912 theta2 = theta * theta; 2913 2914 beta = 1.75; 2915 if (a->IsInRingSize(4) && b->IsInRingSize(4) && c->IsInRingSize(4) && IsInSameRing(a, c)) 2916 beta = 0.85 * beta; 2917 if (a->IsInRingSize(3) && b->IsInRingSize(3) && c->IsInRingSize(3) && IsInSameRing(a, c)) 2918 beta = 0.05 * beta; 2919 2920 // Theta2 is in Degrees^2, but parameters are expecting radians 2921 // PR#2741669 2922 anglecalc.ka = (beta * Za * Cb * Zc * exp(-2 * D)) / (rr * theta2 * DEG_TO_RAD * DEG_TO_RAD); 2923 } 2924 2925 anglecalc.a = a; 2926 anglecalc.b = b; 2927 anglecalc.c = c; 2928 anglecalc.at = angletype; 2929 2930 anglecalc.SetupPointers(); 2931 _anglecalculations.push_back(anglecalc); 2932 2933 if (anglecalc.linear) 2934 continue; 2935 2936 parameter = GetTypedParameter3Atom(strbndtype, type_a, type_b, type_c, _ffstrbndparams); 2937 if (parameter == nullptr) { 2938 int rowa, rowb, rowc; 2939 2940 // This is not a real empirical rule... 2941 //IF_OBFF_LOGLVL_LOW { 2942 // snprintf(_logbuf, BUFF_SIZE, " USING EMPIRICAL RULE FOR STRETCH-BENDING FORCE CONSTANT %d-%d-%d (IDX)...\n", a->GetIdx(), b->GetIdx(), c->GetIdx()); 2943 // OBFFLog(_logbuf); 2944 //} 2945 2946 rowa = GetElementRow(a); 2947 rowb = GetElementRow(b); 2948 rowc = GetElementRow(c); 2949 2950 parameter = GetParameter3Atom(rowa, rowb, rowc, _ffdfsbparams); 2951 2952 if (parameter == nullptr) { 2953 // This should never happen 2954 IF_OBFF_LOGLVL_LOW { 2955 snprintf(_logbuf, BUFF_SIZE, " COULD NOT FIND PARAMETERS FOR STRETCH-BEND %d-%d-%d (IDX)...\n", a->GetIdx(), b->GetIdx(), c->GetIdx()); 2956 OBFFLog(_logbuf); 2957 } 2958 return false; 2959 } 2960 2961 if (rowa == parameter->a) { 2962 strbndcalc.kbaABC = parameter->_dpar[0]; 2963 strbndcalc.kbaCBA = parameter->_dpar[1]; 2964 } else { 2965 strbndcalc.kbaABC = parameter->_dpar[1]; 2966 strbndcalc.kbaCBA = parameter->_dpar[0]; 2967 } 2968 } else { 2969 if (type_a == parameter->a) { 2970 strbndcalc.kbaABC = parameter->_dpar[0]; 2971 strbndcalc.kbaCBA = parameter->_dpar[1]; 2972 } else { 2973 strbndcalc.kbaABC = parameter->_dpar[1]; 2974 strbndcalc.kbaCBA = parameter->_dpar[0]; 2975 } 2976 } 2977 2978 strbndcalc.rab0 = GetBondLength(a, b); 2979 strbndcalc.rbc0 = GetBondLength(b ,c); 2980 strbndcalc.a = a; 2981 strbndcalc.b = b; 2982 strbndcalc.c = c; 2983 strbndcalc.sbt = strbndtype; 2984 strbndcalc.SetupPointers(); 2985 // Set the pointers to addresses in the anglecalc, find the matching bondcalcs and do the same. 2986 // This should improve performance by not calculating all this twice. We could do the same 2987 // for torsion and angles since the bond lengths are calculated for bond stretching first. 2988 //bool found_angle = false; 2989 /* 2990 for (unsigned int ai = 0; ai < _anglecalculations.size(); ++ai) { 2991 if ( (_anglecalculations[ai].a->GetIdx() == a->GetIdx()) && 2992 (_anglecalculations[ai].b->GetIdx() == b->GetIdx()) && 2993 (_anglecalculations[ai].c->GetIdx() == c->GetIdx()) ) { 2994 strbndcalc.theta = &(_anglecalculations[ai].theta); 2995 strbndcalc.force_abc_a = _anglecalculations[ai].force_a; 2996 strbndcalc.force_abc_b = _anglecalculations[ai].force_b; 2997 strbndcalc.force_abc_c = _anglecalculations[ai].force_c; 2998 found_angle = true; 2999 break; 3000 } else if ( (_anglecalculations[ai].a->GetIdx() == c->GetIdx()) && 3001 (_anglecalculations[ai].b->GetIdx() == b->GetIdx()) && 3002 (_anglecalculations[ai].c->GetIdx() == a->GetIdx()) ) { 3003 strbndcalc.theta = &(_anglecalculations[ai].theta); 3004 strbndcalc.force_abc_a = _anglecalculations[ai].force_c; 3005 strbndcalc.force_abc_b = _anglecalculations[ai].force_b; 3006 strbndcalc.force_abc_c = _anglecalculations[ai].force_a; 3007 found_angle = true; 3008 break; 3009 } 3010 } 3011 */ 3012 3013 /* 3014 vector<OBFFAngleCalculationMMFF94>::iterator ai; 3015 for (ai = _anglecalculations.begin(); ai != _anglecalculations.end(); ++ai) { 3016 if ( (((*ai).a)->GetIdx() == a->GetIdx()) && (((*ai).b)->GetIdx() == b->GetIdx()) && (((*ai).c)->GetIdx() == c->GetIdx()) ) { 3017 strbndcalc.theta = (*ai).theta; 3018 cout << "theta prt = " << (*ai).theta << endl; 3019 cout << "delta prt = " << &((*ai).delta) << endl; 3020 cout << "GetThetaPointer = " << ai->GetThetaPointer() << endl; 3021 strbndcalc.force_abc_a = (*ai).force_a; 3022 strbndcalc.force_abc_b = (*ai).force_b; 3023 strbndcalc.force_abc_c = (*ai).force_c; 3024 found_angle = true; 3025 break; 3026 } else if ( (((*ai).a)->GetIdx() == c->GetIdx()) && (((*ai).b)->GetIdx() == b->GetIdx()) && (((*ai).c)->GetIdx() == a->GetIdx()) ) { 3027 strbndcalc.theta = (*ai).theta; 3028 strbndcalc.force_abc_a = (*ai).force_c; 3029 strbndcalc.force_abc_b = (*ai).force_b; 3030 strbndcalc.force_abc_c = (*ai).force_a; 3031 found_angle = true; 3032 break; 3033 } 3034 } 3035 if (!found_angle) // didn't find matching angle, shouldn't happen, but continue to be safe 3036 continue; 3037 3038 3039 bool found_rab = false; 3040 bool found_rbc = false; 3041 vector<OBFFBondCalculationMMFF94>::iterator bi; 3042 for (bi = _bondcalculations.begin(); bi != _bondcalculations.end(); ++bi) { 3043 // find rab 3044 if ( (((*bi).a)->GetIdx() == a->GetIdx()) && (((*bi).b)->GetIdx() == b->GetIdx()) ) { 3045 strbndcalc.rab = &((*bi).rab); 3046 strbndcalc.force_ab_a = (*bi).force_a; 3047 strbndcalc.force_ab_b = (*bi).force_b; 3048 found_rab = true; 3049 } else if ( (((*bi).a)->GetIdx() == b->GetIdx()) && (((*bi).b)->GetIdx() == a->GetIdx()) ) { 3050 strbndcalc.rab = &((*bi).rab); 3051 strbndcalc.force_ab_a = (*bi).force_b; 3052 strbndcalc.force_ab_b = (*bi).force_a; 3053 found_rab = true; 3054 } 3055 // find rbc 3056 if ( (((*bi).a)->GetIdx() == b->GetIdx()) && (((*bi).b)->GetIdx() == c->GetIdx()) ) { 3057 strbndcalc.rbc = &(bondcalc.rab); 3058 strbndcalc.force_ab_a = (*bi).force_a; 3059 strbndcalc.force_ab_b = (*bi).force_b; 3060 found_rbc = true; 3061 } else if ( (((*bi).a)->GetIdx() == c->GetIdx()) && (((*bi).b)->GetIdx() == b->GetIdx()) ) { 3062 strbndcalc.rbc = &(bondcalc.rab); 3063 strbndcalc.force_ab_a = (*bi).force_b; 3064 strbndcalc.force_ab_b = (*bi).force_a; 3065 found_rbc = true; 3066 } 3067 3068 if (found_rab && found_rbc) 3069 break; 3070 } 3071 if (!found_rab || !found_rbc) // didn't find matching bond, or atoms overlap 3072 continue; 3073 */ 3074 _strbndcalculations.push_back(strbndcalc); 3075 3076 } 3077 3078 // 3079 // Torsion Calculations 3080 // 3081 // MMFF part I - page 513 ("step-down" prodedure) 3082 // MMFF part I - page 519 (reference 68 is actually a footnote) 3083 // MMFF part IV - page 631 (empirical rule) 3084 // 3085 // First try and find an exact match, if this fails, step down using the equivalences from mmffdef.par 3086 // five-stage protocol: 1-1-1-1, 2-2-2-2, 3-2-2-5, 5-2-2-3, 5-2-2-5 3087 // If this fails, use empirical rules 3088 // Since 1-1-1-1 = 2-2-2-2, we will only try 1-1-1-1 before going to 3-2-2-5 3089 // 3090 IF_OBFF_LOGLVL_LOW 3091 OBFFLog("SETTING UP TORSION CALCULATIONS...\n"); 3092 3093 OBFFTorsionCalculationMMFF94 torsioncalc; 3094 int torsiontype; 3095 3096 _torsioncalculations.clear(); 3097 3098 FOR_TORSIONS_OF_MOL(t, _mol) { 3099 a = _mol.GetAtom((*t)[0] + 1); 3100 b = _mol.GetAtom((*t)[1] + 1); 3101 c = _mol.GetAtom((*t)[2] + 1); 3102 d = _mol.GetAtom((*t)[3] + 1); 3103 3104 type_a = atoi(a->GetType()); 3105 type_b = atoi(b->GetType()); 3106 type_c = atoi(c->GetType()); 3107 type_d = atoi(d->GetType()); 3108 3109 // skip this torsion if the atoms are ignored 3110 if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) || 3111 _constraints.IsIgnored(c->GetIdx()) || _constraints.IsIgnored(d->GetIdx()) ) 3112 continue; 3113 3114 // if there are any groups specified, check if the four torsion atoms are in a single intraGroup 3115 if (HasGroups()) { 3116 bool validTorsion = false; 3117 for (unsigned int i=0; i < _intraGroup.size(); ++i) { 3118 if (_intraGroup[i].BitIsSet(a->GetIdx()) && _intraGroup[i].BitIsSet(b->GetIdx()) && 3119 _intraGroup[i].BitIsSet(c->GetIdx()) && _intraGroup[i].BitIsSet(d->GetIdx())) { 3120 validTorsion = true; 3121 break; 3122 } 3123 } 3124 if (!validTorsion) 3125 continue; 3126 } 3127 3128 torsiontype = GetTorsionType(a, b, c, d); 3129 // CXT = MC*(J*MA**3 + K*MA**2 + I*MA + L) + TTijkl MC = 6, MA = 136 3130 order = (type_c*2515456 + type_b*18496 + type_d*136 + type_a) 3131 - (type_b*2515456 + type_c*18496 + type_a*136 + type_d); 3132 3133 if (order >= 0) { 3134 // try exact match 3135 parameter = GetTypedParameter4Atom(torsiontype, type_a, type_b, type_c, type_d, _fftorsionparams); 3136 if (parameter == nullptr) // try 3-2-2-5 3137 parameter = GetTypedParameter4Atom(torsiontype, EqLvl3(type_a), type_b, type_c, EqLvl5(type_d), _fftorsionparams); 3138 if (parameter == nullptr) // try 5-2-2-3 3139 parameter = GetTypedParameter4Atom(torsiontype, EqLvl5(type_a), type_b, type_c, EqLvl3(type_d), _fftorsionparams); 3140 if (parameter == nullptr) // try 5-2-2-5 3141 parameter = GetTypedParameter4Atom(torsiontype, EqLvl5(type_a), type_b, type_c, EqLvl5(type_d), _fftorsionparams); 3142 } else { 3143 // try exact match 3144 parameter = GetTypedParameter4Atom(torsiontype, type_d, type_c, type_b, type_a, _fftorsionparams); 3145 if (parameter == nullptr) // try 3-2-2-5 3146 parameter = GetTypedParameter4Atom(torsiontype, EqLvl3(type_d), type_c, type_b, EqLvl5(type_a), _fftorsionparams); 3147 if (parameter == nullptr) // try 5-2-2-3 3148 parameter = GetTypedParameter4Atom(torsiontype, EqLvl5(type_d), type_c, type_b, EqLvl3(type_a), _fftorsionparams); 3149 if (parameter == nullptr) // try 5-2-2-5 3150 parameter = GetTypedParameter4Atom(torsiontype, EqLvl5(type_d), type_c, type_b, EqLvl5(type_a), _fftorsionparams); 3151 } 3152 3153 if (parameter) { 3154 torsioncalc.v1 = parameter->_dpar[0]; 3155 torsioncalc.v2 = parameter->_dpar[1]; 3156 torsioncalc.v3 = parameter->_dpar[2]; 3157 } else { 3158 bool found_rule = false; 3159 3160 //IF_OBFF_LOGLVL_LOW { 3161 // snprintf(_logbuf, BUFF_SIZE, " USING EMPIRICAL RULE FOR TORSION FORCE CONSTANT %d-%d-%d-%d (IDX)...\n", 3162 // a->GetIdx(), b->GetIdx(), c->GetIdx(), d->GetIdx()); 3163 // OBFFLog(_logbuf); 3164 //} 3165 3166 // rule (a) page 631 3167 if (HasLinSet(type_b) || HasLinSet(type_c)) 3168 continue; 3169 3170 // rule (b) page 631 3171 if (b->GetBond(c)->IsAromatic()) { 3172 double Ub, Uc, pi_bc, beta; 3173 Ub = GetUParam(b); 3174 Uc = GetUParam(c); 3175 3176 if (!HasPilpSet(type_b) && !HasPilpSet(type_c)) 3177 pi_bc = 0.5; 3178 else 3179 pi_bc = 0.3; 3180 3181 if (((GetVal(type_b) == 3) && (GetVal(type_c) == 4)) || 3182 ((GetVal(type_b) == 4) && (GetVal(type_c) == 3))) 3183 beta = 3.0; 3184 else 3185 beta = 6.0; 3186 3187 torsioncalc.v1 = 0.0; 3188 torsioncalc.v2 = beta * pi_bc * sqrt(Ub * Uc); 3189 torsioncalc.v3 = 0.0; 3190 found_rule = true; 3191 } else { 3192 // rule (c) page 631 3193 double Ub, Uc, pi_bc, beta; 3194 Ub = GetUParam(b); 3195 Uc = GetUParam(c); 3196 OBBond *bond = a->GetBond(b); 3197 if (((GetMltb(type_b) == 2) && (GetMltb(type_c) == 2)) && !bond->IsAromatic() && bond->GetBondOrder() == 2) 3198 pi_bc = 1.0; 3199 else 3200 pi_bc = 0.4; 3201 3202 beta = 6.0; 3203 torsioncalc.v1 = 0.0; 3204 torsioncalc.v2 = beta * pi_bc * sqrt(Ub * Uc); 3205 torsioncalc.v3 = 0.0; 3206 found_rule = true; 3207 } 3208 3209 // rule (d) page 632 3210 if (!found_rule) 3211 if (((GetCrd(type_b) == 4) && (GetCrd(type_c) == 4))) { 3212 double Vb, Vc; 3213 Vb = GetVParam(b); 3214 Vc = GetVParam(c); 3215 3216 torsioncalc.v1 = 0.0; 3217 torsioncalc.v2 = 0.0; 3218 torsioncalc.v3 = sqrt(Vb * Vc) / 9.0; 3219 found_rule = true; 3220 } 3221 3222 // rule (e) page 632 3223 if (!found_rule) 3224 if (((GetCrd(type_b) == 4) && (GetCrd(type_c) != 4))) { 3225 if (GetCrd(type_c) == 3) // case (1) 3226 if ((GetVal(type_c) == 4) || (GetVal(type_c) == 34) || (GetMltb(type_c) != 0)) 3227 continue; 3228 3229 if (GetCrd(type_c) == 2) // case (2) 3230 if ((GetVal(type_c) == 3) || (GetMltb(type_c) != 0)) 3231 continue; 3232 3233 // case (3) saturated bonds -- see rule (h) 3234 } 3235 3236 // rule (f) page 632 3237 if (!found_rule) 3238 if (((GetCrd(type_b) != 4) && (GetCrd(type_c) == 4))) { 3239 if (GetCrd(type_b) == 3) // case (1) 3240 if ((GetVal(type_b) == 4) || (GetVal(type_b) == 34) || (GetMltb(type_b) != 0)) 3241 continue; 3242 3243 if (GetCrd(type_b) == 2) // case (2) 3244 if ((GetVal(type_b) == 3) || (GetMltb(type_b) != 0)) 3245 continue; 3246 3247 // case (3) saturated bonds 3248 } 3249 3250 // rule (g) page 632 3251 if (!found_rule) { 3252 OBBond *bond = b->GetBond(c); 3253 if (!bond->IsAromatic() && bond->GetBondOrder() == 1 && ( 3254 (GetMltb(type_b) && GetMltb(type_c)) || 3255 (GetMltb(type_b) && HasPilpSet(type_c)) || 3256 (GetMltb(type_c) && HasPilpSet(type_b)))) { 3257 if (HasPilpSet(type_b) && HasPilpSet(type_c)) // case (1) 3258 continue; 3259 3260 double Ub, Uc, pi_bc, beta; 3261 Ub = GetUParam(b); 3262 Uc = GetUParam(c); 3263 beta = 6.0; 3264 3265 if (HasPilpSet(type_b) && GetMltb(type_c)) { // case (2) 3266 if (GetMltb(type_c) == 1) 3267 pi_bc = 0.5; 3268 else if ((GetElementRow(b) == 1) && (GetElementRow(c) == 1)) 3269 pi_bc = 0.3; 3270 else 3271 pi_bc = 0.15; 3272 found_rule = true; 3273 } 3274 3275 if (HasPilpSet(type_c) && GetMltb(type_b)) { // case (3) 3276 if (GetMltb(type_b) == 1) 3277 pi_bc = 0.5; 3278 else if ((GetElementRow(b) == 1) && (GetElementRow(c) == 1)) 3279 pi_bc = 0.3; 3280 else 3281 pi_bc = 0.15; 3282 found_rule = true; 3283 } 3284 3285 if (!found_rule) 3286 if (((GetMltb(type_b) == 1) || (GetMltb(type_c) == 1)) && (b->GetAtomicNum() != OBElements::Carbon || c->GetAtomicNum() != OBElements::Carbon)) { 3287 pi_bc = 0.4; 3288 found_rule = true; 3289 } 3290 3291 if (!found_rule) 3292 pi_bc = 0.15; 3293 3294 torsioncalc.v1 = 0.0; 3295 torsioncalc.v2 = beta * pi_bc * sqrt(Ub * Uc); 3296 torsioncalc.v3 = 0.0; 3297 found_rule = true; 3298 } 3299 } 3300 3301 // rule (h) page 632 3302 if (!found_rule) { 3303 if ((b->GetAtomicNum() == OBElements::Oxygen || b->GetAtomicNum() == OBElements::Sulfur) && (c->GetAtomicNum() == OBElements::Oxygen || c->GetAtomicNum() == OBElements::Sulfur)) { 3304 double Wb, Wc; 3305 3306 if (b->GetAtomicNum() == OBElements::Oxygen) { 3307 Wb = 2.0; 3308 } 3309 else { 3310 Wb = 8.0; 3311 } 3312 3313 if (c->GetAtomicNum() == OBElements::Oxygen) { 3314 Wc = 2.0; 3315 } 3316 else { 3317 Wc = 8.0; 3318 } 3319 3320 torsioncalc.v1 = 0.0; 3321 torsioncalc.v2 = -sqrt(Wb * Wc); 3322 torsioncalc.v3 = 0.0; 3323 } else { 3324 double Vb, Vc, Nbc; 3325 Vb = GetVParam(b); 3326 Vc = GetVParam(c); 3327 3328 IF_OBFF_LOGLVL_LOW { 3329 snprintf(_logbuf, BUFF_SIZE, " USING EMPIRICAL RULE FOR TORSION FORCE CONSTANT %d-%d-%d-%d (IDX)...\n", 3330 a->GetIdx(), b->GetIdx(), c->GetIdx(), d->GetIdx()); 3331 OBFFLog(_logbuf); 3332 } 3333 3334 Nbc = GetCrd(type_b) * GetCrd(type_c); 3335 3336 torsioncalc.v1 = 0.0; 3337 torsioncalc.v2 = 0.0; 3338 torsioncalc.v3 = sqrt(Vb * Vc) / Nbc; 3339 } 3340 } 3341 } 3342 3343 torsioncalc.a = a; 3344 torsioncalc.b = b; 3345 torsioncalc.c = c; 3346 torsioncalc.d = d; 3347 torsioncalc.SetupPointers(); 3348 torsioncalc.tt = torsiontype; 3349 3350 _torsioncalculations.push_back(torsioncalc); 3351 } 3352 3353 // 3354 // Out-Of-Plane Calculations 3355 // 3356 IF_OBFF_LOGLVL_LOW 3357 OBFFLog("SETTING UP OOP CALCULATIONS...\n"); 3358 3359 OBFFOOPCalculationMMFF94 oopcalc; 3360 3361 _oopcalculations.clear(); 3362 3363 FOR_ATOMS_OF_MOL(atom, _mol) { 3364 b = (OBAtom*) &*atom; 3365 3366 found = false; 3367 3368 type_b = atoi(b->GetType()); 3369 3370 for (unsigned int idx=0; idx < _ffoopparams.size(); idx++) { 3371 if (type_b == _ffoopparams[idx].b) { 3372 a = nullptr; 3373 c = nullptr; 3374 d = nullptr; 3375 3376 FOR_NBORS_OF_ATOM(nbr, b) { 3377 if (a ==nullptr) 3378 a = (OBAtom*) &*nbr; 3379 else if (c == nullptr) 3380 c = (OBAtom*) &*nbr; 3381 else 3382 d = (OBAtom*) &*nbr; 3383 } 3384 3385 if (a == nullptr || c == nullptr || d == nullptr) 3386 break; 3387 3388 type_a = atoi(a->GetType()); 3389 type_c = atoi(c->GetType()); 3390 type_d = atoi(d->GetType()); 3391 3392 // skip this oop if the atoms are ignored 3393 if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) || 3394 _constraints.IsIgnored(c->GetIdx()) || _constraints.IsIgnored(d->GetIdx()) ) 3395 continue; 3396 3397 // if there are any groups specified, check if the four oop atoms are in a single intraGroup 3398 if (HasGroups()) { 3399 bool validOOP = false; 3400 for (unsigned int i=0; i < _intraGroup.size(); ++i) { 3401 if (_intraGroup[i].BitIsSet(a->GetIdx()) && _intraGroup[i].BitIsSet(b->GetIdx()) && 3402 _intraGroup[i].BitIsSet(c->GetIdx()) && _intraGroup[i].BitIsSet(d->GetIdx())) { 3403 validOOP = true; 3404 break; 3405 } 3406 } 3407 if (!validOOP) 3408 continue; 3409 } 3410 3411 if (((type_a == _ffoopparams[idx].a) && (type_c == _ffoopparams[idx].c) && (type_d == _ffoopparams[idx].d)) || 3412 ((type_c == _ffoopparams[idx].a) && (type_a == _ffoopparams[idx].c) && (type_d == _ffoopparams[idx].d)) || 3413 ((type_c == _ffoopparams[idx].a) && (type_d == _ffoopparams[idx].c) && (type_a == _ffoopparams[idx].d)) || 3414 ((type_d == _ffoopparams[idx].a) && (type_c == _ffoopparams[idx].c) && (type_a == _ffoopparams[idx].d)) || 3415 ((type_a == _ffoopparams[idx].a) && (type_d == _ffoopparams[idx].c) && (type_c == _ffoopparams[idx].d)) || 3416 ((type_d == _ffoopparams[idx].a) && (type_a == _ffoopparams[idx].c) && (type_c == _ffoopparams[idx].d))) 3417 { 3418 found = true; 3419 3420 oopcalc.koop = _ffoopparams[idx]._dpar[0]; 3421 3422 // A-B-CD || C-B-AD PLANE = ABC 3423 oopcalc.a = a; 3424 oopcalc.b = b; 3425 oopcalc.c = c; 3426 oopcalc.d = d; 3427 3428 oopcalc.SetupPointers(); 3429 _oopcalculations.push_back(oopcalc); 3430 3431 // C-B-DA || D-B-CA PLANE BCD 3432 oopcalc.a = d; 3433 oopcalc.d = a; 3434 3435 oopcalc.SetupPointers(); 3436 _oopcalculations.push_back(oopcalc); 3437 3438 // A-B-DC || D-B-AC PLANE ABD 3439 oopcalc.a = a; 3440 oopcalc.c = d; 3441 oopcalc.d = c; 3442 3443 oopcalc.SetupPointers(); 3444 _oopcalculations.push_back(oopcalc); 3445 } 3446 3447 if ((_ffoopparams[idx].a == 0) && (_ffoopparams[idx].c == 0) && (_ffoopparams[idx].d == 0) && !found) // *-XX-*-* 3448 { 3449 oopcalc.koop = _ffoopparams[idx]._dpar[0]; 3450 3451 // A-B-CD || C-B-AD PLANE = ABC 3452 oopcalc.a = a; 3453 oopcalc.b = b; 3454 oopcalc.c = c; 3455 oopcalc.d = d; 3456 3457 oopcalc.SetupPointers(); 3458 _oopcalculations.push_back(oopcalc); 3459 3460 // C-B-DA || D-B-CA PLANE BCD 3461 oopcalc.a = d; 3462 oopcalc.d = a; 3463 3464 oopcalc.SetupPointers(); 3465 _oopcalculations.push_back(oopcalc); 3466 3467 // A-B-DC || D-B-AC PLANE ABD 3468 oopcalc.a = a; 3469 oopcalc.c = d; 3470 oopcalc.d = c; 3471 3472 oopcalc.SetupPointers(); 3473 _oopcalculations.push_back(oopcalc); 3474 } 3475 } 3476 } 3477 } 3478 3479 // 3480 // VDW Calculations 3481 // 3482 IF_OBFF_LOGLVL_LOW 3483 OBFFLog("SETTING UP VAN DER WAALS CALCULATIONS...\n"); 3484 3485 OBFFVDWCalculationMMFF94 vdwcalc; 3486 3487 _vdwcalculations.clear(); 3488 3489 int pairIndex = -1; 3490 FOR_PAIRS_OF_MOL(p, _mol) { 3491 ++pairIndex; 3492 a = _mol.GetAtom((*p)[0]); 3493 b = _mol.GetAtom((*p)[1]); 3494 3495 // skip this vdw if the atoms are ignored 3496 if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) ) 3497 continue; 3498 3499 // if there are any groups specified, check if the two atoms are in a single _interGroup or if 3500 // two two atoms are in one of the _interGroups pairs. 3501 if (HasGroups()) { 3502 bool validVDW = false; 3503 for (unsigned int i=0; i < _interGroup.size(); ++i) { 3504 if (_interGroup[i].BitIsSet(a->GetIdx()) && _interGroup[i].BitIsSet(b->GetIdx())) { 3505 validVDW = true; 3506 break; 3507 } 3508 } 3509 if (!validVDW) { 3510 for (unsigned int i=0; i < _interGroups.size(); ++i) { 3511 if (_interGroups[i].first.BitIsSet(a->GetIdx()) && _interGroups[i].second.BitIsSet(b->GetIdx())) { 3512 validVDW = true; 3513 break; 3514 } 3515 if (_interGroups[i].first.BitIsSet(b->GetIdx()) && _interGroups[i].second.BitIsSet(a->GetIdx())) { 3516 validVDW = true; 3517 break; 3518 } 3519 } 3520 } 3521 3522 if (!validVDW) 3523 continue; 3524 } 3525 3526 OBFFParameter *parameter_a, *parameter_b; 3527 parameter_a = GetParameter1Atom(atoi(a->GetType()), _ffvdwparams); 3528 parameter_b = GetParameter1Atom(atoi(b->GetType()), _ffvdwparams); 3529 if (parameter_a == nullptr || parameter_b == nullptr) { 3530 IF_OBFF_LOGLVL_LOW { 3531 snprintf(_logbuf, BUFF_SIZE, " COULD NOT FIND VAN DER WAALS PARAMETERS FOR %d-%d (IDX)...\n", a->GetIdx(), b->GetIdx()); 3532 OBFFLog(_logbuf); 3533 } 3534 3535 return false; 3536 } 3537 3538 vdwcalc.a = a; 3539 vdwcalc.alpha_a = parameter_a->_dpar[0]; 3540 vdwcalc.Na = parameter_a->_dpar[1]; 3541 vdwcalc.Aa = parameter_a->_dpar[2]; 3542 vdwcalc.Ga = parameter_a->_dpar[3]; 3543 vdwcalc.aDA = parameter_a->_ipar[0]; 3544 3545 vdwcalc.b = b; 3546 vdwcalc.alpha_b = parameter_b->_dpar[0]; 3547 vdwcalc.Nb = parameter_b->_dpar[1]; 3548 vdwcalc.Ab = parameter_b->_dpar[2]; 3549 vdwcalc.Gb = parameter_b->_dpar[3]; 3550 vdwcalc.bDA = parameter_b->_ipar[0]; 3551 3552 //these calculations only need to be done once for each pair, 3553 //we do them now and save them for later use 3554 double R_AA, R_BB, R_AB6, g_AB, g_AB2; 3555 double R_AB2, R_AB4, /*R_AB7,*/ sqrt_a, sqrt_b; 3556 3557 R_AA = vdwcalc.Aa * pow(vdwcalc.alpha_a, 0.25); 3558 R_BB = vdwcalc.Ab * pow(vdwcalc.alpha_b, 0.25); 3559 sqrt_a = sqrt(vdwcalc.alpha_a / vdwcalc.Na); 3560 sqrt_b = sqrt(vdwcalc.alpha_b / vdwcalc.Nb); 3561 3562 if (vdwcalc.aDA == 1) { // hydrogen bond donor 3563 vdwcalc.R_AB = 0.5 * (R_AA + R_BB); 3564 R_AB2 = vdwcalc.R_AB * vdwcalc.R_AB; 3565 R_AB4 = R_AB2 * R_AB2; 3566 R_AB6 = R_AB4 * R_AB2; 3567 3568 if (vdwcalc.bDA == 2) { // hydrogen bond acceptor 3569 vdwcalc.epsilon = 0.5 * (181.16 * vdwcalc.Ga * vdwcalc.Gb * vdwcalc.alpha_a * vdwcalc.alpha_b) / (sqrt_a + sqrt_b) * (1.0 / R_AB6); 3570 // R_AB is scaled to 0.8 for D-A interactions. The value used in the calculation of epsilon is not scaled. 3571 vdwcalc.R_AB = 0.8 * vdwcalc.R_AB; 3572 } else 3573 vdwcalc.epsilon = (181.16 * vdwcalc.Ga * vdwcalc.Gb * vdwcalc.alpha_a * vdwcalc.alpha_b) / (sqrt_a + sqrt_b) * (1.0 / R_AB6); 3574 3575 R_AB2 = vdwcalc.R_AB * vdwcalc.R_AB; 3576 R_AB4 = R_AB2 * R_AB2; 3577 R_AB6 = R_AB4 * R_AB2; 3578 vdwcalc.R_AB7 = R_AB6 * vdwcalc.R_AB; 3579 } else if (vdwcalc.bDA == 1) { // hydrogen bond donor 3580 vdwcalc.R_AB = 0.5 * (R_AA + R_BB); 3581 R_AB2 = vdwcalc.R_AB * vdwcalc.R_AB; 3582 R_AB4 = R_AB2 * R_AB2; 3583 R_AB6 = R_AB4 * R_AB2; 3584 3585 if (vdwcalc.aDA == 2) { // hydrogen bond acceptor 3586 vdwcalc.epsilon = 0.5 * (181.16 * vdwcalc.Ga * vdwcalc.Gb * vdwcalc.alpha_a * vdwcalc.alpha_b) / (sqrt_a + sqrt_b) * (1.0 / R_AB6); 3587 // R_AB is scaled to 0.8 for D-A interactions. The value used in the calculation of epsilon is not scaled. 3588 vdwcalc.R_AB = 0.8 * vdwcalc.R_AB; 3589 } else 3590 vdwcalc.epsilon = (181.16 * vdwcalc.Ga * vdwcalc.Gb * vdwcalc.alpha_a * vdwcalc.alpha_b) / (sqrt_a + sqrt_b) * (1.0 / R_AB6); 3591 3592 R_AB2 = vdwcalc.R_AB * vdwcalc.R_AB; 3593 R_AB4 = R_AB2 * R_AB2; 3594 R_AB6 = R_AB4 * R_AB2; 3595 vdwcalc.R_AB7 = R_AB6 * vdwcalc.R_AB; 3596 } else { 3597 g_AB = (R_AA - R_BB) / ( R_AA + R_BB); 3598 g_AB2 = g_AB * g_AB; 3599 vdwcalc.R_AB = 0.5 * (R_AA + R_BB) * (1.0 + 0.2 * (1.0 - exp(-12.0 * g_AB2))); 3600 R_AB2 = vdwcalc.R_AB * vdwcalc.R_AB; 3601 R_AB4 = R_AB2 * R_AB2; 3602 R_AB6 = R_AB4 * R_AB2; 3603 vdwcalc.R_AB7 = R_AB6 * vdwcalc.R_AB; 3604 vdwcalc.epsilon = (181.16 * vdwcalc.Ga * vdwcalc.Gb * vdwcalc.alpha_a * vdwcalc.alpha_b) / (sqrt_a + sqrt_b) * (1.0 / R_AB6); 3605 } 3606 3607 vdwcalc.pairIndex = pairIndex; 3608 vdwcalc.SetupPointers(); 3609 _vdwcalculations.push_back(vdwcalc); 3610 } 3611 3612 // 3613 // Electrostatic Calculations 3614 // 3615 IF_OBFF_LOGLVL_LOW 3616 OBFFLog("SETTING UP ELECTROSTATIC CALCULATIONS...\n"); 3617 3618 OBFFElectrostaticCalculationMMFF94 elecalc; 3619 3620 _electrostaticcalculations.clear(); 3621 3622 pairIndex = -1; 3623 FOR_PAIRS_OF_MOL(p, _mol) { 3624 ++pairIndex; 3625 a = _mol.GetAtom((*p)[0]); 3626 b = _mol.GetAtom((*p)[1]); 3627 3628 // skip this ele if the atoms are ignored 3629 if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) ) 3630 continue; 3631 3632 // if there are any groups specified, check if the two atoms are in a single _interGroup or if 3633 // two two atoms are in one of the _interGroups pairs. 3634 if (HasGroups()) { 3635 bool validEle = false; 3636 for (unsigned int i=0; i < _interGroup.size(); ++i) { 3637 if (_interGroup[i].BitIsSet(a->GetIdx()) && _interGroup[i].BitIsSet(b->GetIdx())) { 3638 validEle = true; 3639 break; 3640 } 3641 } 3642 if (!validEle) { 3643 for (unsigned int i=0; i < _interGroups.size(); ++i) { 3644 if (_interGroups[i].first.BitIsSet(a->GetIdx()) && _interGroups[i].second.BitIsSet(b->GetIdx())) { 3645 validEle = true; 3646 break; 3647 } 3648 if (_interGroups[i].first.BitIsSet(b->GetIdx()) && _interGroups[i].second.BitIsSet(a->GetIdx())) { 3649 validEle = true; 3650 break; 3651 } 3652 } 3653 } 3654 3655 if (!validEle) 3656 continue; 3657 } 3658 3659 elecalc.qq = 332.0716 * a->GetPartialCharge() * b->GetPartialCharge() / _epsilon; 3660 3661 if (elecalc.qq) { 3662 elecalc.a = &*a; 3663 elecalc.b = &*b; 3664 3665 // 1-4 scaling 3666 if (a->IsOneFour(b)) 3667 elecalc.qq *= 0.75; 3668 3669 elecalc.pairIndex = pairIndex; 3670 elecalc.SetupPointers(); 3671 _electrostaticcalculations.push_back(elecalc); 3672 } 3673 } 3674 3675 return true; 3676 } 3677 SetupPointers()3678 bool OBForceFieldMMFF94::SetupPointers() 3679 { 3680 for (unsigned int i = 0; i < _bondcalculations.size(); ++i) 3681 _bondcalculations[i].SetupPointers(); 3682 for (unsigned int i = 0; i < _anglecalculations.size(); ++i) 3683 _anglecalculations[i].SetupPointers(); 3684 for (unsigned int i = 0; i < _strbndcalculations.size(); ++i) 3685 _strbndcalculations[i].SetupPointers(); 3686 for (unsigned int i = 0; i < _torsioncalculations.size(); ++i) 3687 _torsioncalculations[i].SetupPointers(); 3688 for (unsigned int i = 0; i < _oopcalculations.size(); ++i) 3689 _oopcalculations[i].SetupPointers(); 3690 for (unsigned int i = 0; i < _vdwcalculations.size(); ++i) 3691 _vdwcalculations[i].SetupPointers(); 3692 for (unsigned int i = 0; i < _electrostaticcalculations.size(); ++i) 3693 _electrostaticcalculations[i].SetupPointers(); 3694 3695 return true; 3696 } 3697 3698 3699 // we set the the formal charge with SetPartialCharge because formal charges 3700 // in MMFF94 are not always and integer SetFormalCharges()3701 bool OBForceFieldMMFF94::SetFormalCharges() 3702 { 3703 _mol.SetAutomaticPartialCharge(false); 3704 3705 FOR_ATOMS_OF_MOL (atom, _mol) { 3706 int type = atoi(atom->GetType()); 3707 atom->SetPartialCharge(0.0); 3708 3709 bool done = false; 3710 switch (type) { 3711 case 34: 3712 case 49: 3713 case 51: 3714 case 54: 3715 case 58: 3716 case 92: 3717 case 93: 3718 case 94: 3719 case 97: 3720 atom->SetPartialCharge(1.0); 3721 done = true; 3722 break; 3723 case 35: 3724 case 62: 3725 case 89: 3726 case 90: 3727 case 91: 3728 atom->SetPartialCharge(-1.0); 3729 done = true; 3730 break; 3731 case 55: 3732 atom->SetPartialCharge(0.5); 3733 done = true; 3734 break; 3735 case 87: 3736 case 95: 3737 case 96: 3738 case 98: 3739 case 99: 3740 atom->SetPartialCharge(2.0); 3741 done = true; 3742 break; 3743 case 88: 3744 atom->SetPartialCharge(3.0); 3745 done = true; 3746 break; 3747 //case 98: 3748 // atom->SetPartialCharge(3.0); 3749 default: 3750 break; 3751 } 3752 3753 if (done) 3754 continue; 3755 3756 if (type == 56) { 3757 int n_count = 0; 3758 int temp_type; 3759 FOR_ATOMS_OF_MOL (atom2, _mol) { 3760 temp_type = atoi(atom2->GetType()); 3761 if (temp_type == 56 || temp_type == 81) 3762 ++n_count; 3763 } 3764 atom->SetPartialCharge((double)((n_count + 1) / 3) / (double)n_count); 3765 } else if (type == 32) { 3766 int o_count = 0; 3767 bool sulfonamide = false; 3768 bool sulfone_s_c = false; 3769 int s_count = 0; 3770 3771 FOR_NBORS_OF_ATOM(nbr, &*atom) { 3772 FOR_NBORS_OF_ATOM(nbr2, &*nbr) { 3773 if (nbr2->GetAtomicNum() == OBElements::Oxygen && (nbr2->GetExplicitDegree() == 1)) 3774 o_count++; 3775 if (nbr2->GetAtomicNum() == OBElements::Sulfur && (nbr2->GetExplicitDegree() == 1)) 3776 s_count++; 3777 if (nbr2->GetAtomicNum() == OBElements::Nitrogen && !nbr2->IsAromatic()) 3778 sulfonamide = true; 3779 OBBond *bond = nbr->GetBond(&*nbr2); 3780 if (nbr2->GetAtomicNum() == OBElements::Carbon && !bond->IsAromatic() && bond->GetBondOrder() == 2) 3781 sulfone_s_c = true; 3782 } 3783 3784 if (nbr->GetAtomicNum() == OBElements::Carbon) 3785 atom->SetPartialCharge(-0.5); // O2CM 3786 3787 if (nbr->GetAtomicNum() == OBElements::Nitrogen && (o_count == 3)) 3788 atom->SetPartialCharge(-1.0 / o_count); // O3N 3789 3790 if (nbr->GetAtomicNum() == OBElements::Sulfur && !sulfonamide) { 3791 if (((o_count + s_count) == 2) && (nbr->GetExplicitDegree() == 3) 3792 && (nbr->GetExplicitValence() >= 3) && !sulfone_s_c) { 3793 atom->SetPartialCharge(-0.5); // O2S (sulfinate) 3794 } 3795 else if ((o_count + s_count) == 3) { 3796 atom->SetPartialCharge(-1.0 / 3.0); // O3S 3797 } 3798 else if ((o_count + s_count) == 4) { 3799 atom->SetPartialCharge(-0.5); // O4S 3800 } 3801 } 3802 3803 if (nbr->GetAtomicNum() == OBElements::Phosphorus) { 3804 if ((o_count + s_count) == 2) { 3805 atom->SetPartialCharge(-0.5); // O2P 3806 } 3807 else if ((o_count + s_count) == 3) { 3808 atom->SetPartialCharge(-2.0 / 3.0); // O3P 3809 } 3810 else if ((o_count + s_count) == 4) { 3811 atom->SetPartialCharge(-0.25); // O4P 3812 } 3813 } 3814 3815 if (atoi(nbr->GetType()) == 77) 3816 atom->SetPartialCharge(-0.25); // O4CL 3817 } 3818 } else if (type == 61) { 3819 FOR_BONDS_OF_ATOM(bond, &*atom) { 3820 OBAtom *nbr = bond->GetNbrAtom(&*atom); 3821 if (!bond->IsAromatic() && bond->GetBondOrder() == 3 && nbr->GetAtomicNum() == OBElements::Nitrogen) 3822 atom->SetPartialCharge(1.0); 3823 } 3824 } else if (type == 72) { 3825 int s_count = 0; 3826 3827 FOR_NBORS_OF_ATOM(nbr, &*atom) { 3828 if (nbr->GetAtomicNum() == OBElements::Sulfur) 3829 s_count++; 3830 3831 if (nbr->GetAtomicNum() == OBElements::Phosphorus || nbr->GetAtomicNum() == OBElements::Sulfur) { 3832 FOR_NBORS_OF_ATOM(nbr2, &*nbr) 3833 if ((nbr2->GetAtomicNum() == OBElements::Sulfur || nbr2->GetAtomicNum() == OBElements::Oxygen) && (nbr2->GetExplicitDegree() == 1) && (atom->GetIdx() != nbr2->GetIdx())) 3834 atom->SetPartialCharge(-0.5); 3835 } else 3836 atom->SetPartialCharge(-1.0); 3837 3838 if (nbr->GetAtomicNum() == OBElements::Carbon) 3839 FOR_NBORS_OF_ATOM(nbr2, &*nbr) 3840 if (nbr2->GetAtomicNum() == OBElements::Sulfur && (nbr2->GetExplicitDegree() == 1) && (atom->GetIdx() != nbr2->GetIdx())) 3841 atom->SetPartialCharge(-0.5); // SSMO 3842 3843 if (s_count >= 2) 3844 atom->SetPartialCharge(-0.5); // SSMO 3845 } 3846 } else if (type == 76) { 3847 vector<OBRing*> vr; 3848 vr = _mol.GetSSSR(); 3849 vector<OBRing*>::iterator ri; 3850 vector<int>::iterator rj; 3851 int n_count; 3852 3853 for (ri = vr.begin();ri != vr.end();++ri) { // for each ring 3854 n_count = 0; 3855 3856 if ((*ri)->IsAromatic() && (*ri)->IsMember(&*atom) && ((*ri)->Size() == 5)) { 3857 for(rj = (*ri)->_path.begin();rj != (*ri)->_path.end();++rj) // for each ring atom 3858 if (_mol.GetAtom(*rj)->GetAtomicNum() == OBElements::Nitrogen) 3859 n_count++; 3860 3861 if (n_count > 1) 3862 atom->SetPartialCharge(-1.0 / n_count); 3863 } 3864 } 3865 } else if (type == 81) { 3866 atom->SetPartialCharge(1.0); 3867 3868 vector<OBRing*> vr; 3869 vr = _mol.GetSSSR(); 3870 vector<OBRing*>::iterator ri; 3871 vector<int>::iterator rj; 3872 for (ri = vr.begin();ri != vr.end();++ri) // for each ring 3873 if ((*ri)->IsAromatic() && (*ri)->IsMember(&*atom) && ((*ri)->Size() == 5)) { 3874 int n_count = 0; 3875 for(rj = (*ri)->_path.begin();rj != (*ri)->_path.end();++rj) // for each ring atom 3876 if (_mol.GetAtom(*rj)->GetAtomicNum() == OBElements::Nitrogen && (_mol.GetAtom(*rj)->GetExplicitDegree() == 3)) 3877 n_count++; 3878 3879 if (n_count) // coverity defensive testing 3880 atom->SetPartialCharge(1.0 / n_count); // NIM+ 3881 3882 FOR_NBORS_OF_ATOM(nbr, &*atom) 3883 FOR_NBORS_OF_ATOM(nbr2, &*nbr) 3884 if (atoi(nbr2->GetType()) == 56) 3885 atom->SetPartialCharge(1.0 / 3.0); 3886 3887 FOR_NBORS_OF_ATOM(nbr, &*atom) 3888 FOR_NBORS_OF_ATOM(nbr2, &*nbr) 3889 if (atoi(nbr2->GetType()) == 55) 3890 atom->SetPartialCharge(1.0 / (1.0 + n_count)); 3891 } 3892 } 3893 3894 } 3895 3896 PrintFormalCharges(); 3897 3898 return true; 3899 } 3900 SetPartialCharges()3901 bool OBForceFieldMMFF94::SetPartialCharges() 3902 { 3903 vector<double> charges(_mol.NumAtoms()+1, 0); 3904 double M, Wab, factor, q0a, q0b, Pa, Pb; 3905 3906 FOR_ATOMS_OF_MOL (atom, _mol) { 3907 int type = atoi(atom->GetType()); 3908 3909 switch (type) { 3910 case 32: 3911 case 35: 3912 case 72: 3913 factor = 0.5; 3914 break; 3915 case 62: 3916 case 76: 3917 factor = 0.25; 3918 break; 3919 default: 3920 factor = 0.0; 3921 break; 3922 } 3923 3924 M = GetCrd(type); 3925 q0a = atom->GetPartialCharge(); 3926 3927 // charge sharing 3928 if (!factor) 3929 FOR_NBORS_OF_ATOM (nbr, &*atom) 3930 if (nbr->GetPartialCharge() < 0.0) 3931 q0a += nbr->GetPartialCharge() / (2.0 * (double)(nbr->GetExplicitDegree())); 3932 3933 // needed for SEYWUO, positive charge sharing? 3934 if (type == 62) 3935 FOR_NBORS_OF_ATOM (nbr, &*atom) 3936 if (nbr->GetPartialCharge() > 0.0) 3937 q0a -= nbr->GetPartialCharge() / 2.0; 3938 3939 q0b = 0.0; 3940 Wab = 0.0; 3941 Pa = Pb = 0.0; 3942 FOR_NBORS_OF_ATOM (nbr, &*atom) { 3943 int nbr_type = atoi(nbr->GetType()); 3944 3945 q0b += nbr->GetPartialCharge(); 3946 3947 bool bci_found = false; 3948 for (unsigned int idx=0; idx < _ffchgparams.size(); idx++) 3949 if (GetBondType(&*atom, &*nbr) == _ffchgparams[idx]._ipar[0]) { 3950 if ((type == _ffchgparams[idx].a) && (nbr_type == _ffchgparams[idx].b)) { 3951 Wab += -_ffchgparams[idx]._dpar[0]; 3952 bci_found = true; 3953 } else if ((type == _ffchgparams[idx].b) && (nbr_type == _ffchgparams[idx].a)) { 3954 Wab += _ffchgparams[idx]._dpar[0]; 3955 bci_found = true; 3956 } 3957 } 3958 3959 if (!bci_found) { 3960 for (unsigned int idx=0; idx < _ffpbciparams.size(); idx++) { 3961 if (type == _ffpbciparams[idx].a) 3962 Pa = _ffpbciparams[idx]._dpar[0]; 3963 if (nbr_type == _ffpbciparams[idx].a) 3964 Pb = _ffpbciparams[idx]._dpar[0]; 3965 } 3966 Wab += Pa - Pb; 3967 } 3968 } 3969 if (factor) 3970 charges[atom->GetIdx()] = (1.0 - M * factor) * q0a + factor * q0b + Wab; 3971 else 3972 charges[atom->GetIdx()] = q0a + Wab; 3973 } 3974 3975 FOR_ATOMS_OF_MOL (atom, _mol) 3976 atom->SetPartialCharge(charges[atom->GetIdx()]); 3977 3978 PrintPartialCharges(); 3979 3980 return true; 3981 } 3982 3983 //////////////////////////////////////////////////////////////////////////////// 3984 //////////////////////////////////////////////////////////////////////////////// 3985 // 3986 // Validation functions 3987 // 3988 //////////////////////////////////////////////////////////////////////////////// 3989 //////////////////////////////////////////////////////////////////////////////// 3990 3991 // used to validate the implementation Validate()3992 bool OBForceFieldMMFF94::Validate () 3993 { 3994 OBConversion conv; 3995 OBFormat *format_in = conv.FindFormat("mol2"); 3996 vector<string> vs; 3997 vector<int> types; 3998 vector<double> fcharges, pcharges; 3999 vector<double> bond_lengths; 4000 char buffer[150]; 4001 bool molfound, atomfound, bondfound, fchgfound, pchgfound; 4002 double etot, ebond, eangle, eoop, estbn, etor, evdw, eeq; 4003 double termcount; //1=bond, 2=angle, 3=strbnd, 4=torsion, 5=oop 4004 int n = 0; 4005 4006 if (!format_in || !conv.SetInFormat(format_in)) { 4007 obErrorLog.ThrowError(__FUNCTION__, "Could not set mol2 input format", obError); 4008 return false; 4009 } 4010 4011 ifstream ifs, ifs2; 4012 ofstream ofs; 4013 4014 ifs.open("MMFF94_dative.mol2"); 4015 if (!ifs) { 4016 obErrorLog.ThrowError(__FUNCTION__, "Could not open ./MMFF94_dative.mol2", obError); 4017 return false; 4018 } 4019 4020 ifs2.open("MMFF94_opti.log"); 4021 if (!ifs2) { 4022 obErrorLog.ThrowError(__FUNCTION__, "Coulg not open ./MMFF_opti.log", obError); 4023 return false; 4024 } 4025 4026 ofs.open("MMFF94_openbabel.log"); 4027 if (!ofs) { 4028 obErrorLog.ThrowError(__FUNCTION__, "Coulg not open ./MMFF_openbabel.log", obError); 4029 return false; 4030 } 4031 4032 if (!_init) { 4033 ParseParamFile(); 4034 _init = true; 4035 } 4036 4037 4038 SetLogFile(&ofs); 4039 SetLogLevel(OBFF_LOGLVL_HIGH); 4040 4041 for (unsigned int c=1;; c++) { 4042 _mol.Clear(); 4043 types.clear(); 4044 fcharges.clear(); 4045 pcharges.clear(); 4046 bond_lengths.clear(); 4047 4048 if (!conv.Read(&_mol, &ifs)) 4049 break; 4050 if (_mol.Empty()) 4051 break; 4052 4053 _ncoords = _mol.NumAtoms() * 3; 4054 _gradientPtr = new double[_ncoords]; 4055 4056 SetTypes(); 4057 4058 if ((c == 98) || (c == 692)) // CUDPAS & VUWXUG 4059 continue; 4060 4061 termcount = 0; 4062 molfound = false; 4063 atomfound = false; 4064 bondfound = false; 4065 fchgfound = false; 4066 pchgfound = false; 4067 4068 while (ifs2.getline(buffer, 150)) { 4069 tokenize(vs, buffer); 4070 if (vs.size() == 0) { 4071 bondfound = false; 4072 continue; 4073 } 4074 4075 string str(buffer); 4076 if (string::npos != str.find(_mol.GetTitle(),0)) 4077 molfound = true; 4078 4079 if (atomfound) { 4080 if (n) { 4081 types.push_back(atoi(vs[2].c_str())); 4082 types.push_back(atoi(vs[5].c_str())); 4083 types.push_back(atoi(vs[8].c_str())); 4084 types.push_back(atoi(vs[11].c_str())); 4085 } else { 4086 if (vs.size() > 2) 4087 types.push_back(atoi(vs[2].c_str())); 4088 if (vs.size() > 5) 4089 types.push_back(atoi(vs[5].c_str())); 4090 if (vs.size() > 8) 4091 types.push_back(atoi(vs[8].c_str())); 4092 4093 atomfound = false; 4094 } 4095 n--; 4096 } 4097 4098 if (fchgfound) { 4099 if (n) { 4100 fcharges.push_back(atof(vs[2].c_str())); 4101 fcharges.push_back(atof(vs[5].c_str())); 4102 fcharges.push_back(atof(vs[8].c_str())); 4103 fcharges.push_back(atof(vs[11].c_str())); 4104 } else { 4105 if (vs.size() > 2) 4106 fcharges.push_back(atof(vs[2].c_str())); 4107 if (vs.size() > 5) 4108 fcharges.push_back(atof(vs[5].c_str())); 4109 if (vs.size() > 8) 4110 fcharges.push_back(atof(vs[8].c_str())); 4111 4112 fchgfound = false; 4113 } 4114 n--; 4115 } 4116 4117 if (pchgfound) { 4118 if (n) { 4119 pcharges.push_back(atof(vs[2].c_str())); 4120 pcharges.push_back(atof(vs[5].c_str())); 4121 pcharges.push_back(atof(vs[8].c_str())); 4122 pcharges.push_back(atof(vs[11].c_str())); 4123 } else { 4124 if (vs.size() > 2) 4125 pcharges.push_back(atof(vs[2].c_str())); 4126 if (vs.size() > 5) 4127 pcharges.push_back(atof(vs[5].c_str())); 4128 if (vs.size() > 8) 4129 pcharges.push_back(atof(vs[8].c_str())); 4130 4131 pchgfound = false; 4132 } 4133 n--; 4134 } 4135 4136 if (molfound && EQn(buffer, " ATOM NAME TYPE", 16)) { 4137 atomfound = true; 4138 n = _mol.NumAtoms() / 4; 4139 } 4140 if (molfound && EQn(buffer, " ATOM FCHARGE", 17)) { 4141 fchgfound = true; 4142 n = _mol.NumAtoms() / 4; 4143 } 4144 if (molfound && EQn(buffer, " ATOM CHARGE", 17)) { 4145 pchgfound = true; 4146 n = _mol.NumAtoms() / 4; 4147 } 4148 4149 if (bondfound) 4150 bond_lengths.push_back(atof(vs[7].c_str())); 4151 4152 if (molfound) { 4153 if (EQn(buffer, " Total ENERGY", 13)) 4154 etot = atof(vs[3].c_str()); 4155 if (EQn(buffer, " Bond Stretching", 16)) 4156 ebond = atof(vs[2].c_str()); 4157 if (EQn(buffer, " Angle Bending", 14)) 4158 eangle = atof(vs[2].c_str()); 4159 if (EQn(buffer, " Out-of-Plane Bending", 21)) 4160 eoop = atof(vs[2].c_str()); 4161 if (EQn(buffer, " Stretch-Bend", 13)) 4162 estbn = atof(vs[1].c_str()); 4163 if (EQn(buffer, " Total Torsion", 18)) 4164 etor = atof(vs[2].c_str()); 4165 if (EQn(buffer, " Net vdW", 12)) 4166 evdw = atof(vs[2].c_str()); 4167 if (EQn(buffer, " Electrostatic", 14)) 4168 eeq = atof(vs[1].c_str()); 4169 if (EQn(buffer, " ---------------------", 22) && (termcount == 0)) { 4170 termcount++; 4171 bondfound = true; 4172 } 4173 if (EQn(buffer, " OPTIMOL> # read next", 22)) 4174 break; 4175 } 4176 4177 4178 4179 } // while (getline) 4180 4181 vector<int>::iterator i; 4182 vector<double>::iterator di; 4183 unsigned int ni; 4184 bool failed; 4185 4186 cout << "--------------------------------------------------------------------------------" << endl; 4187 cout << " " << endl; 4188 cout << " VALIDATE MOLECULE " << c << ": " << _mol.GetTitle() << endl; 4189 cout << " " << endl; 4190 cout << "IDX HYB AROM OB_TYPE LOG_TYPE RESULT " << endl; 4191 cout << "---------------------------------------------- " << endl; 4192 4193 ni = 1; 4194 failed = false; 4195 for (i = types.begin(); i != types.end();++i) { 4196 if (ni > _mol.NumAtoms()) 4197 continue; 4198 4199 if ( (atoi(_mol.GetAtom(ni)->GetType()) == 87) || 4200 (atoi(_mol.GetAtom(ni)->GetType()) == 97) 4201 ) continue; 4202 4203 if (atoi(_mol.GetAtom(ni)->GetType()) == (*i)) 4204 snprintf(_logbuf, BUFF_SIZE, "%2d %3d %4d %3d %3d PASSED", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetHyb(), 4205 _mol.GetAtom(ni)->IsAromatic(), atoi(_mol.GetAtom(ni)->GetType()), *i); 4206 else { 4207 snprintf(_logbuf, BUFF_SIZE, "%2d %3d %4d %3d %3d XXX FAILED XXX", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetHyb(), 4208 _mol.GetAtom(ni)->IsAromatic(), atoi(_mol.GetAtom(ni)->GetType()), *i); 4209 failed = true; 4210 } 4211 4212 cout << _logbuf << endl; 4213 4214 ni++; 4215 } 4216 4217 if (failed) { 4218 cout << "Could not successfully assign atom types" << endl; 4219 return false; 4220 //continue; 4221 } 4222 4223 SetFormalCharges(); 4224 cout << endl; 4225 cout << "IDX OB_FCARGE LOG_FCHARGE RESULT" << endl; 4226 cout << "----------------------------------------" << endl; 4227 4228 ni = 1; 4229 for (di = fcharges.begin(); di != fcharges.end(); ++di) { 4230 if (ni > _mol.NumAtoms()) 4231 continue; 4232 4233 if ( (atoi(_mol.GetAtom(ni)->GetType()) == 87) || 4234 (atoi(_mol.GetAtom(ni)->GetType()) == 97) 4235 ) continue; 4236 4237 if (fabs((*di) - _mol.GetAtom(ni)->GetPartialCharge()) <= 0.001) 4238 snprintf(_logbuf, BUFF_SIZE, "%2d %7.4f %7.4f PASSED", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetPartialCharge(), *di); 4239 else { 4240 snprintf(_logbuf, BUFF_SIZE, "%2d %7.4f %7.4f XXX FAILED XXX", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetPartialCharge(), *di); 4241 failed = true; 4242 } 4243 4244 cout << _logbuf << endl; 4245 4246 ni++; 4247 } 4248 4249 if (failed) { 4250 cout << "Could not successfully assign formal charges" << endl; 4251 //return false; 4252 continue; 4253 } 4254 4255 SetPartialCharges(); 4256 cout << endl; 4257 cout << "IDX OB_PCARGE LOG_PCHARGE RESULT" << endl; 4258 cout << "----------------------------------------" << endl; 4259 4260 ni = 1; 4261 for (di = pcharges.begin(); di != pcharges.end(); ++di) { 4262 if (ni > _mol.NumAtoms()) 4263 continue; 4264 4265 if ( (atoi(_mol.GetAtom(ni)->GetType()) == 87) || 4266 (atoi(_mol.GetAtom(ni)->GetType()) == 97) 4267 ) continue; 4268 4269 if (fabs((*di) - _mol.GetAtom(ni)->GetPartialCharge()) <= 0.001) 4270 snprintf(_logbuf, BUFF_SIZE, "%2d %7.4f %7.4f PASSED", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetPartialCharge(), *di); 4271 else { 4272 snprintf(_logbuf, BUFF_SIZE, "%2d %7.4f %7.4f XXX FAILED XXX", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetPartialCharge(), *di); 4273 failed = true; 4274 } 4275 4276 cout << _logbuf << endl; 4277 4278 ni++; 4279 } 4280 4281 if (failed) { 4282 cout << "Could not successfully assign partial charges" << endl; 4283 //return false; 4284 continue; 4285 } 4286 4287 4288 4289 if (!SetupCalculations()) { 4290 cout << "Could not setup calculations (missing parameters...)" << endl; 4291 return false; 4292 //continue; 4293 } 4294 4295 double delta; 4296 cout << endl; 4297 cout << "TERM OB ENERGY LOG ENERGY DELTA" << endl; 4298 cout << "---------------------------------------------------------------" << endl; 4299 4300 delta = (E_Bond() - ebond); 4301 snprintf(_logbuf, BUFF_SIZE, "Bond Stretching %11.5f %11.5f %11.5f", E_Bond(), ebond, delta); 4302 cout << _logbuf << endl; 4303 4304 delta = (E_Angle() - eangle); 4305 snprintf(_logbuf, BUFF_SIZE, "Angle Bending %11.5f %11.5f %11.5f", E_Angle(), eangle, delta); 4306 cout << _logbuf << endl; 4307 4308 delta = (E_StrBnd() - estbn); 4309 snprintf(_logbuf, BUFF_SIZE, "Stretch-Bending %11.5f %11.5f %11.5f", E_StrBnd(), estbn, delta); 4310 cout << _logbuf << endl; 4311 4312 delta = (E_OOP() - eoop); 4313 snprintf(_logbuf, BUFF_SIZE, "Out-Of-Plane Bending %11.5f %11.5f %11.5f", E_OOP(), eoop, delta); 4314 cout << _logbuf << endl; 4315 4316 delta = (E_Torsion() - etor); 4317 snprintf(_logbuf, BUFF_SIZE, "Torsional %11.5f %11.5f %11.5f", E_Torsion(), etor, delta); 4318 cout << _logbuf << endl; 4319 4320 delta = (E_VDW() - evdw); 4321 snprintf(_logbuf, BUFF_SIZE, "Van der Waals %11.5f %11.5f %11.5f", E_VDW(), evdw, delta); 4322 cout << _logbuf << endl; 4323 4324 delta = (E_Electrostatic() - eeq); 4325 snprintf(_logbuf, BUFF_SIZE, "Electrostatic %11.5f %11.5f %11.5f", E_Electrostatic(), eeq, delta); 4326 cout << _logbuf << endl; 4327 4328 cout << endl; 4329 delta = (Energy() - etot); 4330 snprintf(_logbuf, BUFF_SIZE, "Total ENERGY %11.5f %11.5f %11.5f", Energy(), etot, delta); 4331 cout << _logbuf << endl; 4332 4333 } // for (unsigned int c;; c++ ) 4334 4335 if (ifs) 4336 ifs.close(); 4337 if (ifs2) 4338 ifs2.close(); 4339 4340 return true; 4341 } 4342 ValidateGradients()4343 bool OBForceFieldMMFF94::ValidateGradients () 4344 { 4345 vector3 numgrad, anagrad, err; 4346 int coordIdx; 4347 4348 bool passed = true; // set to false if any component fails 4349 4350 cout << "----------------------------------------------------------------------------------------" << endl; 4351 cout << " " << endl; 4352 cout << " VALIDATE GRADIENTS : " << _mol.GetTitle() << endl; 4353 cout << " " << endl; 4354 cout << " " << endl; 4355 cout << "ATOM IDX NUMERICAL GRADIENT ANALYTICAL GRADIENT REL. ERROR (%) " << endl; 4356 cout << "----------------------------------------------------------------------------------------" << endl; 4357 // "XX (000.000, 000.000, 000.000) (000.000, 000.000, 000.000) (00.00, 00.00, 00.00)" 4358 4359 FOR_ATOMS_OF_MOL (a, _mol) { 4360 coordIdx = (a->GetIdx() - 1) * 3; 4361 4362 // OBFF_ENERGY 4363 numgrad = NumericalDerivative(&*a, OBFF_ENERGY); 4364 Energy(); // compute 4365 anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]); 4366 err = ValidateGradientError(numgrad, anagrad); 4367 4368 snprintf(_logbuf, BUFF_SIZE, "%2d (%7.3f, %7.3f, %7.3f) (%7.3f, %7.3f, %7.3f) (%5.2f, %5.2f, %5.2f)\n", a->GetIdx(), numgrad.x(), numgrad.y(), numgrad.z(), 4369 anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z()); 4370 OBFFLog(_logbuf); 4371 4372 // OBFF_EBOND 4373 numgrad = NumericalDerivative(&*a, OBFF_EBOND); 4374 ClearGradients(); 4375 E_Bond(); // compute 4376 anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]); 4377 err = ValidateGradientError(numgrad, anagrad); 4378 4379 snprintf(_logbuf, BUFF_SIZE, " bond (%7.3f, %7.3f, %7.3f) (%7.3f, %7.3f, %7.3f) (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(), 4380 anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z()); 4381 OBFFLog(_logbuf); 4382 if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0) 4383 passed = false; 4384 4385 // OBFF_EANGLE 4386 numgrad = NumericalDerivative(&*a, OBFF_EANGLE); 4387 ClearGradients(); 4388 E_Angle(); // compute 4389 anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]); 4390 err = ValidateGradientError(numgrad, anagrad); 4391 4392 snprintf(_logbuf, BUFF_SIZE, " angle (%7.3f, %7.3f, %7.3f) (%7.3f, %7.3f, %7.3f) (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(), 4393 anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z()); 4394 OBFFLog(_logbuf); 4395 if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0) 4396 passed = false; 4397 4398 // OBFF_ESTRBND 4399 numgrad = NumericalDerivative(&*a, OBFF_ESTRBND); 4400 ClearGradients(); 4401 E_StrBnd(); // compute 4402 anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]); 4403 err = ValidateGradientError(numgrad, anagrad); 4404 4405 snprintf(_logbuf, BUFF_SIZE, " strbnd (%7.3f, %7.3f, %7.3f) (%7.3f, %7.3f, %7.3f) (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(), 4406 anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z()); 4407 OBFFLog(_logbuf); 4408 if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0) 4409 passed = false; 4410 4411 // OBFF_ETORSION 4412 numgrad = NumericalDerivative(&*a, OBFF_ETORSION); 4413 ClearGradients(); 4414 E_Torsion(); // compute 4415 anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]); 4416 err = ValidateGradientError(numgrad, anagrad); 4417 4418 snprintf(_logbuf, BUFF_SIZE, " torsion (%7.3f, %7.3f, %7.3f) (%7.3f, %7.3f, %7.3f) (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(), 4419 anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z()); 4420 OBFFLog(_logbuf); 4421 if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0) 4422 passed = false; 4423 4424 // OBFF_EOOP 4425 numgrad = NumericalDerivative(&*a, OBFF_EOOP); 4426 ClearGradients(); 4427 E_OOP(); // compute 4428 anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]); 4429 err = ValidateGradientError(numgrad, anagrad); 4430 4431 snprintf(_logbuf, BUFF_SIZE, " oop (%7.3f, %7.3f, %7.3f) (%7.3f, %7.3f, %7.3f) (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(), 4432 anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z()); 4433 OBFFLog(_logbuf); 4434 // disable OOP gradient validation for now -- some small errors, but nothing major 4435 // if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0) 4436 // passed = false; 4437 4438 // OBFF_EVDW 4439 numgrad = NumericalDerivative(&*a, OBFF_EVDW); 4440 ClearGradients(); 4441 E_VDW(); // compute 4442 anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]); 4443 err = ValidateGradientError(numgrad, anagrad); 4444 4445 snprintf(_logbuf, BUFF_SIZE, " vdw (%7.3f, %7.3f, %7.3f) (%7.3f, %7.3f, %7.3f) (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(), 4446 anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z()); 4447 OBFFLog(_logbuf); 4448 if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0) 4449 passed = false; 4450 4451 // OBFF_EELECTROSTATIC 4452 numgrad = NumericalDerivative(&*a, OBFF_EELECTROSTATIC); 4453 ClearGradients(); 4454 E_Electrostatic(); // compute 4455 anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]); 4456 err = ValidateGradientError(numgrad, anagrad); 4457 4458 snprintf(_logbuf, BUFF_SIZE, " electro (%7.3f, %7.3f, %7.3f) (%7.3f, %7.3f, %7.3f) (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(), 4459 anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z()); 4460 OBFFLog(_logbuf); 4461 if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0) 4462 passed = false; 4463 } 4464 4465 return passed; // did we pass every single component? 4466 } 4467 4468 //////////////////////////////////////////////////////////////////////////////// 4469 //////////////////////////////////////////////////////////////////////////////// 4470 // 4471 // Calculate bond type, angle type, stretch-bend type, torsion type 4472 // 4473 //////////////////////////////////////////////////////////////////////////////// 4474 //////////////////////////////////////////////////////////////////////////////// 4475 4476 // 4477 // MMFF part V - page 620 4478 // 4479 // BTij is 1 when: 4480 // a) single bond between atoms i and j, both i and j are not aromatic and both types have sbmb set in mmffprop.par, or 4481 // b) bewtween two aromatic atoms, but the bond is not aromatic (e.g. connecting bond in biphenyl) 4482 // GetBondType(OBAtom * a,OBAtom * b)4483 int OBForceFieldMMFF94::GetBondType(OBAtom* a, OBAtom* b) 4484 { 4485 OBBond *bond = _mol.GetBond(a, b); 4486 if (bond->GetBondOrder() != 1 || bond->IsAromatic()) 4487 return 0; 4488 4489 if (HasAromSet(atoi(a->GetType())) && HasAromSet(atoi(b->GetType()))) 4490 return 1; 4491 4492 if (HasSbmbSet(atoi(a->GetType())) && HasSbmbSet(atoi(b->GetType()))) 4493 return 1; 4494 4495 return 0; 4496 } 4497 GetAngleType(OBAtom * a,OBAtom * b,OBAtom * c)4498 int OBForceFieldMMFF94::GetAngleType(OBAtom* a, OBAtom* b, OBAtom *c) 4499 { 4500 int sumbondtypes; 4501 4502 sumbondtypes = GetBondType(a,b) + GetBondType(b, c); 4503 4504 if (a->IsInRingSize(3) && b->IsInRingSize(3) && c->IsInRingSize(3) && IsInSameRing(a, c)) 4505 switch (sumbondtypes) { 4506 case 0: 4507 return 3; 4508 case 1: 4509 return 5; 4510 case 2: 4511 return 6; 4512 } 4513 4514 if (a->IsInRingSize(4) && b->IsInRingSize(4) && c->IsInRingSize(4) && IsInSameRing(a, c)) 4515 switch (sumbondtypes) { 4516 case 0: 4517 return 4; 4518 case 1: 4519 return 7; 4520 case 2: 4521 return 8; 4522 } 4523 4524 return sumbondtypes; 4525 } 4526 GetStrBndType(OBAtom * a,OBAtom * b,OBAtom * c)4527 int OBForceFieldMMFF94::GetStrBndType(OBAtom* a, OBAtom* b, OBAtom *c) 4528 { 4529 int btab, btbc, atabc; 4530 bool inverse; 4531 4532 btab = GetBondType(a, b); 4533 btbc = GetBondType(b, c); 4534 atabc = GetAngleType(a, b, c); 4535 4536 if (atoi(a->GetType()) <= atoi(c->GetType())) 4537 inverse = false; 4538 else 4539 inverse = true; 4540 4541 switch (atabc) { 4542 case 0: 4543 return 0; 4544 4545 case 1: 4546 if (btab) { 4547 if (!inverse) { 4548 return 1; 4549 } else { 4550 return 2; 4551 } 4552 } 4553 if (btbc) { 4554 if (!inverse) { 4555 return 2; 4556 } else { 4557 return 1; 4558 } 4559 } 4560 4561 case 2: 4562 return 3; 4563 4564 case 3: 4565 return 5; 4566 4567 case 4: 4568 return 4; 4569 4570 case 5: 4571 if (btab) { 4572 if (!inverse) { 4573 return 6; 4574 } else { 4575 return 7; 4576 } 4577 } 4578 if (btbc) { 4579 if (!inverse) { 4580 return 7; 4581 } else { 4582 return 6; 4583 } 4584 } 4585 4586 case 6: 4587 return 8; 4588 4589 case 7: 4590 if (btab) { 4591 if (!inverse) { 4592 return 9; 4593 } else { 4594 return 10; 4595 } 4596 } 4597 if (btbc) { 4598 if (!inverse) { 4599 return 10; 4600 } else { 4601 return 9; 4602 } 4603 } 4604 4605 case 8: 4606 return 11; 4607 } 4608 4609 return 0; 4610 } 4611 4612 // 4613 // MMFF part IV - page 609 4614 // 4615 // TTijkl = 1 when BTjk = 1 4616 // TTijkl = 2 when BTjk = 0 but BTij and/or BTkl = 1 4617 // TTijkl = 4 when i, j, k and l are all members of the same four-membered ring 4618 // TTijkl = 5 when i, j, k and l are members of a five-membered ring and at least one is a sp3-hybridized carbon (MMFF atom type 1) 4619 // GetTorsionType(OBAtom * a,OBAtom * b,OBAtom * c,OBAtom * d)4620 int OBForceFieldMMFF94::GetTorsionType(OBAtom* a, OBAtom* b, OBAtom *c, OBAtom *d) 4621 { 4622 int btab, btbc, btcd; 4623 4624 btab = GetBondType(a, b); 4625 btbc = GetBondType(b, c); 4626 btcd = GetBondType(c, d); 4627 4628 if (btbc == 1) 4629 return 1; 4630 4631 if (a->IsInRingSize(4) && b->IsInRingSize(4) && c->IsInRingSize(4) && d->IsInRingSize(4)) 4632 if (IsInSameRing(a,b) && IsInSameRing(b,c) && IsInSameRing(c,d)) 4633 return 4; 4634 4635 OBBond *bond = _mol.GetBond(b, c); 4636 if (bond->GetBondOrder() == 1 && !bond->IsAromatic()) { 4637 if (btab || btcd) 4638 return 2; 4639 /* 4640 unsigned int order1 = GetCXT(0, atoi(d->GetType()), atoi(c->GetType()), atoi(b->GetType()), atoi(a->GetType())); 4641 unsigned int order2 = GetCXT(0, atoi(a->GetType()), atoi(b->GetType()), atoi(c->GetType()), atoi(d->GetType())); 4642 4643 cout << "GetTorsionType(" << a->GetType() << ", " << b->GetType() << ", " << c->GetType() << ", " << d->GetType() << ")" << endl; 4644 cout << " order1 = " << order1 << endl; 4645 cout << " order2 = " << order2 << endl; 4646 cout << " btab = " << btab << endl; 4647 cout << " btbc = " << btbc << endl; 4648 cout << " btcd = " << btcd << endl; 4649 */ 4650 } 4651 4652 if (a->IsInRingSize(5) && b->IsInRingSize(5) && c->IsInRingSize(5) && d->IsInRingSize(5)) { 4653 vector<OBRing*> vr; 4654 vr = _mol.GetSSSR(); 4655 4656 if( !((atoi(a->GetType()) == 1) || (atoi(b->GetType()) == 1) || (atoi(c->GetType()) == 1) || (atoi(d->GetType()) == 1)) ) 4657 return 0; 4658 4659 vector<OBRing*>::iterator ri; 4660 vector<int>::iterator rj; 4661 for (ri = vr.begin();ri != vr.end();++ri) { // for each ring 4662 if ((*ri)->IsAromatic()) 4663 continue; 4664 4665 if ((*ri)->Size() != 5) 4666 continue; 4667 4668 if (!(*ri)->IsMember(a) || !(*ri)->IsMember(b) || !(*ri)->IsMember(c) || !(*ri)->IsMember(d)) 4669 continue; 4670 4671 return 5; 4672 } 4673 } 4674 4675 4676 return 0; 4677 } 4678 4679 // CXB = MC * (I * MA + J) + BTij GetCXB(int type,int a,int b)4680 unsigned int OBForceFieldMMFF94::GetCXB(int type, int a, int b) 4681 { 4682 unsigned int cxb; 4683 cxb = 2 * (a * 136 + b) + type; 4684 return cxb; 4685 } 4686 4687 // CXA = MC * (J * MA^2 + I * MA + K) + ATijk GetCXA(int type,int a,int b,int c)4688 unsigned int OBForceFieldMMFF94::GetCXA(int type, int a, int b, int c) 4689 { 4690 unsigned int cxa; 4691 cxa = 9 * (b * 18496 + a * 136 + c) + type; 4692 return cxa; 4693 } 4694 4695 // CXS = MC * (J * MA^2 + I * MA + K) + STijk GetCXS(int type,int a,int b,int c)4696 unsigned int OBForceFieldMMFF94::GetCXS(int type, int a, int b, int c) 4697 { 4698 unsigned int cxs; 4699 cxs = 12 * (b * 18496 + a * 136 + c) + type; 4700 return cxs; 4701 } 4702 4703 // CXO = J * MA^3 + I * MA^2 + K * MA + L GetCXO(int a,int b,int c,int d)4704 unsigned int OBForceFieldMMFF94::GetCXO(int a, int b, int c, int d) 4705 { 4706 unsigned int cxo; 4707 cxo = b * 2515456 + a * 18496 + c * 136 + d; 4708 return cxo; 4709 } 4710 4711 // CXT = MC * (J * MA^3 + K * MA^2 + I * MA + L) + TTijkl GetCXT(int type,int a,int b,int c,int d)4712 unsigned int OBForceFieldMMFF94::GetCXT(int type, int a, int b, int c, int d) 4713 { 4714 unsigned int cxt; 4715 cxt = 6 * (b * 2515456 + c * 18496 + a * 136 + d) + type; 4716 return cxt; 4717 } 4718 4719 // CXQ = MC * (I * MA + J) + BTij GetCXQ(int type,int a,int b)4720 unsigned int OBForceFieldMMFF94::GetCXQ(int type, int a, int b) 4721 { 4722 unsigned int cxq; 4723 cxq = 2 * (a * 136 + b) + type; 4724 return cxq; 4725 } 4726 4727 //////////////////////////////////////////////////////////////////////////////// 4728 //////////////////////////////////////////////////////////////////////////////// 4729 // 4730 // Various tables & misc. functions 4731 // 4732 //////////////////////////////////////////////////////////////////////////////// 4733 //////////////////////////////////////////////////////////////////////////////// 4734 4735 // MMFF part V - TABLE I HasLinSet(int atomtype)4736 bool OBForceFieldMMFF94::HasLinSet(int atomtype) 4737 { 4738 return _ffpropLin.BitIsSet(atomtype); 4739 } 4740 4741 // MMFF part V - TABLE I HasPilpSet(int atomtype)4742 bool OBForceFieldMMFF94::HasPilpSet(int atomtype) 4743 { 4744 return _ffpropPilp.BitIsSet(atomtype); 4745 } 4746 4747 // MMFF part V - TABLE I HasAromSet(int atomtype)4748 bool OBForceFieldMMFF94::HasAromSet(int atomtype) 4749 { 4750 return _ffpropArom.BitIsSet(atomtype); 4751 } 4752 4753 // MMFF part V - TABLE I HasSbmbSet(int atomtype)4754 bool OBForceFieldMMFF94::HasSbmbSet(int atomtype) 4755 { 4756 return _ffpropSbmb.BitIsSet(atomtype); 4757 } 4758 4759 // MMFF part V - TABLE I GetCrd(int atomtype)4760 int OBForceFieldMMFF94::GetCrd(int atomtype) 4761 { 4762 OBFFParameter *par; 4763 4764 par = GetParameter1Atom(atomtype, _ffpropparams); // from mmffprop.par 4765 if (par) 4766 return par->_ipar[1]; 4767 4768 return 0; 4769 } 4770 4771 // MMFF part V - TABLE I GetVal(int atomtype)4772 int OBForceFieldMMFF94::GetVal(int atomtype) 4773 { 4774 OBFFParameter *par; 4775 4776 par = GetParameter1Atom(atomtype, _ffpropparams); // from mmffprop.par 4777 if (par) 4778 return par->_ipar[2]; 4779 4780 return 0; 4781 } 4782 4783 // MMFF part V - TABLE I GetMltb(int atomtype)4784 int OBForceFieldMMFF94::GetMltb(int atomtype) 4785 { 4786 OBFFParameter *par; 4787 4788 par = GetParameter1Atom(atomtype, _ffpropparams); // from mmffprop.par 4789 if (par) 4790 return par->_ipar[4]; 4791 4792 return 0; 4793 } 4794 4795 // MMFF part I - TABLE IV EqLvl2(int type)4796 int OBForceFieldMMFF94::EqLvl2(int type) 4797 { 4798 for (unsigned int idx=0; idx < _ffdefparams.size(); idx++) 4799 if (_ffdefparams[idx]._ipar[0] == type) 4800 return _ffdefparams[idx]._ipar[1]; 4801 4802 return type; 4803 } 4804 4805 // MMFF part I - TABLE IV EqLvl3(int type)4806 int OBForceFieldMMFF94::EqLvl3(int type) 4807 { 4808 for (unsigned int idx=0; idx < _ffdefparams.size(); idx++) 4809 if (_ffdefparams[idx]._ipar[0] == type) 4810 return _ffdefparams[idx]._ipar[2]; 4811 4812 return type; 4813 } 4814 4815 // MMFF part I - TABLE IV EqLvl4(int type)4816 int OBForceFieldMMFF94::EqLvl4(int type) 4817 { 4818 for (unsigned int idx=0; idx < _ffdefparams.size(); idx++) 4819 if (_ffdefparams[idx]._ipar[0] == type) 4820 return _ffdefparams[idx]._ipar[3]; 4821 4822 return type; 4823 } 4824 4825 // MMFF part I - TABLE IV EqLvl5(int type)4826 int OBForceFieldMMFF94::EqLvl5(int type) 4827 { 4828 for (unsigned int idx=0; idx < _ffdefparams.size(); idx++) 4829 if (_ffdefparams[idx]._ipar[0] == type) 4830 return _ffdefparams[idx]._ipar[4]; 4831 4832 return type; 4833 } 4834 4835 // MMFF part V - TABLE VI GetZParam(OBAtom * atom)4836 double OBForceFieldMMFF94::GetZParam(OBAtom* atom) 4837 { 4838 switch (atom->GetAtomicNum()) { 4839 case OBElements::Hydrogen: 4840 return 1.395; 4841 case OBElements::Carbon: 4842 return 2.494; 4843 case OBElements::Nitrogen: 4844 return 2.711; 4845 case OBElements::Oxygen: 4846 return 3.045; 4847 case OBElements::Fluorine: 4848 return 2.847; 4849 case OBElements::Silicon: 4850 return 2.350; 4851 case OBElements::Phosphorus: 4852 return 2.350; 4853 case OBElements::Sulfur: 4854 return 2.980; 4855 case OBElements::Chlorine: 4856 return 2.909; 4857 case OBElements::Bromine: 4858 return 3.017; 4859 case OBElements::Iodine: 4860 return 3.086; 4861 } 4862 4863 return 0.0; 4864 } 4865 4866 // MMFF part V - TABLE VI GetCParam(OBAtom * atom)4867 double OBForceFieldMMFF94::GetCParam(OBAtom* atom) 4868 { 4869 switch (atom->GetAtomicNum()) { 4870 case OBElements::Boron: 4871 return 0.704; 4872 case OBElements::Carbon: 4873 return 1.016; 4874 case OBElements::Nitrogen: 4875 return 1.113; 4876 case OBElements::Oxygen: 4877 return 1.337; 4878 case OBElements::Silicon: 4879 return 0.811; 4880 case OBElements::Phosphorus: 4881 return 1.068; 4882 case OBElements::Sulfur: 4883 return 1.249; 4884 case OBElements::Chlorine: 4885 return 1.078; 4886 case OBElements::Arsenic: 4887 return 0.825; 4888 } 4889 4890 return 0.0; 4891 } 4892 4893 // MMFF part V - TABLE X GetUParam(OBAtom * atom)4894 double OBForceFieldMMFF94::GetUParam(OBAtom* atom) 4895 { 4896 switch (atom->GetAtomicNum()) { 4897 case OBElements::Carbon: 4898 return 2.0; 4899 case OBElements::Nitrogen: 4900 return 2.0; 4901 case OBElements::Oxygen: 4902 return 2.0; 4903 case OBElements::Silicon: 4904 return 1.25; 4905 case OBElements::Phosphorus: 4906 return 1.25; 4907 case OBElements::Sulfur: 4908 return 1.25; 4909 } 4910 4911 return 0.0; 4912 } 4913 4914 // MMFF part V - TABLE X GetVParam(OBAtom * atom)4915 double OBForceFieldMMFF94::GetVParam(OBAtom* atom) 4916 { 4917 switch (atom->GetAtomicNum()) { 4918 case OBElements::Carbon: 4919 return 2.12; 4920 case OBElements::Nitrogen: 4921 return 1.5; 4922 case OBElements::Oxygen: 4923 return 0.2; 4924 case OBElements::Silicon: 4925 return 1.22; 4926 case OBElements::Phosphorus: 4927 return 2.4; 4928 case OBElements::Sulfur: 4929 return 0.49; 4930 } 4931 4932 return 0.0; 4933 } 4934 4935 // R Blom and A Haaland, J. Mol. Struct., 128, 21-27 (1985) GetCovalentRadius(OBAtom * a)4936 double OBForceFieldMMFF94::GetCovalentRadius(OBAtom* a) { 4937 4938 switch (a->GetAtomicNum()) { 4939 case 1: 4940 return 0.33; // corrected value from MMFF part V 4941 case 5: 4942 return 0.81; 4943 case 6: 4944 return 0.77; // corrected value from MMFF part V 4945 case 7: 4946 return 0.73; 4947 case 8: 4948 return 0.72; 4949 case 9: 4950 return 0.74; 4951 case 13: 4952 return 1.22; 4953 case 14: 4954 return 1.15; 4955 case 15: 4956 return 1.09; 4957 case 16: 4958 return 1.03; 4959 case 17: 4960 return 1.01; 4961 case 31: 4962 return 1.19; 4963 case 32: 4964 return 1.20; 4965 case 33: 4966 return 1.20; 4967 case 34: 4968 return 1.16; 4969 case 35: 4970 return 1.15; 4971 case 44: 4972 return 1.46; 4973 case 50: 4974 return 1.40; 4975 case 51: 4976 return 1.41; 4977 case 52: 4978 return 1.35; 4979 case 53: 4980 return 1.33; 4981 case 81: 4982 return 1.51; 4983 case 82: 4984 return 1.53; 4985 case 83: 4986 return 1.55; 4987 default: 4988 return OBElements::GetCovalentRad(a->GetAtomicNum()); 4989 } 4990 } 4991 GetBondLength(OBAtom * a,OBAtom * b)4992 double OBForceFieldMMFF94::GetBondLength(OBAtom* a, OBAtom* b) 4993 { 4994 OBFFParameter *parameter; 4995 double rab; 4996 4997 parameter = GetTypedParameter2Atom(GetBondType(a, b), atoi(a->GetType()), atoi(b->GetType()), _ffbondparams); 4998 if (parameter == nullptr) 4999 rab = GetRuleBondLength(a, b); 5000 else 5001 rab = parameter->_dpar[1]; 5002 5003 return rab; 5004 } 5005 5006 // MMFF part V - page 625 GetRuleBondLength(OBAtom * a,OBAtom * b)5007 double OBForceFieldMMFF94::GetRuleBondLength(OBAtom* a, OBAtom* b) 5008 { 5009 double r0ab, r0a, r0b, c, Xa, Xb; 5010 int Ha, Hb, BOab; 5011 r0a = GetCovalentRadius(a); 5012 r0b = GetCovalentRadius(b); 5013 Xa = OBElements::GetAllredRochowElectroNeg(a->GetAtomicNum()); 5014 Xb = OBElements::GetAllredRochowElectroNeg(b->GetAtomicNum()); 5015 5016 5017 if (a->GetAtomicNum() == OBElements::Hydrogen) 5018 r0a = 0.33; 5019 if (b->GetAtomicNum() == OBElements::Hydrogen) 5020 r0b = 0.33; 5021 5022 if (a->GetAtomicNum() == OBElements::Hydrogen || b->GetAtomicNum() == OBElements::Hydrogen) 5023 c = 0.050; 5024 else 5025 c = 0.085; 5026 5027 if (GetMltb(atoi(a->GetType()) == 3)) 5028 Ha = 1; 5029 else if ((GetMltb(atoi(a->GetType())) == 1) || (GetMltb(atoi(a->GetType())) == 2)) 5030 Ha = 2; 5031 else 5032 Ha = 3; 5033 5034 if (GetMltb(atoi(b->GetType()) == 3)) 5035 Hb = 1; 5036 else if ((GetMltb(atoi(b->GetType())) == 1) || (GetMltb(atoi(b->GetType())) == 2)) 5037 Hb = 2; 5038 else 5039 Hb = 3; 5040 5041 BOab = a->GetBond(b)->GetBondOrder(); 5042 if ((GetMltb(atoi(a->GetType())) == 1) && (GetMltb(atoi(b->GetType())) == 1)) 5043 BOab = 4; 5044 if ((GetMltb(atoi(a->GetType())) == 1) && (GetMltb(atoi(b->GetType())) == 2)) 5045 BOab = 5; 5046 if ((GetMltb(atoi(a->GetType())) == 2) && (GetMltb(atoi(b->GetType())) == 1)) 5047 BOab = 5; 5048 if (a->GetBond(b)->IsAromatic()) { 5049 if (!HasPilpSet(atoi(a->GetType())) && !HasPilpSet(atoi(b->GetType()))) { 5050 BOab = 4; 5051 } else { 5052 BOab = 5; 5053 } 5054 } 5055 5056 switch (BOab) { 5057 case 5: 5058 r0a -= 0.04; 5059 r0b -= 0.04; 5060 break; 5061 case 4: 5062 r0a -= 0.075; 5063 r0b -= 0.075; 5064 break; 5065 case 3: 5066 r0a -= 0.17; 5067 r0b -= 0.17; 5068 break; 5069 case 2: 5070 r0a -= 0.10; 5071 r0b -= 0.10; 5072 break; 5073 case 1: 5074 if (Ha == 1) 5075 r0a -= 0.08; 5076 if (Ha == 2) 5077 r0a -= 0.03; 5078 if (Hb == 1) 5079 r0b -= 0.08; 5080 if (Hb == 2) 5081 r0b -= 0.03; 5082 } 5083 5084 /* 5085 cout << "Ha=" << Ha << " Hb=" << Hb << " BOab=" << BOab << endl; 5086 cout << "r0a=" << r0a << " Xa=" << Xa << endl; 5087 cout << "r0b=" << r0b << " Xb=" << Xb << endl; 5088 cout << "r0a + r0b=" << r0a +r0b << endl; 5089 cout << "c=" << c << " |Xa-Xb|=" << fabs(Xa-Xb) << " |Xa-Xb|^1.4=" << pow(fabs(Xa-Xb), 1.4) << endl; 5090 */ 5091 r0ab = r0a + r0b - c * pow(fabs(Xa - Xb), 1.4) - 0.008; 5092 5093 return r0ab; 5094 } 5095 GetParameter1Atom(int a,std::vector<OBFFParameter> & parameter)5096 OBFFParameter* OBForceFieldMMFF94::GetParameter1Atom(int a, std::vector<OBFFParameter> ¶meter) 5097 { 5098 OBFFParameter *par; 5099 5100 for (unsigned int idx=0; idx < parameter.size(); idx++) 5101 if (a == parameter[idx].a) { 5102 par = ¶meter[idx]; 5103 return par; 5104 } 5105 5106 return nullptr; 5107 } 5108 GetParameter2Atom(int a,int b,std::vector<OBFFParameter> & parameter)5109 OBFFParameter* OBForceFieldMMFF94::GetParameter2Atom(int a, int b, std::vector<OBFFParameter> ¶meter) 5110 { 5111 OBFFParameter *par; 5112 5113 for (unsigned int idx=0; idx < parameter.size(); idx++) 5114 if (((a == parameter[idx].a) && (b == parameter[idx].b)) || 5115 ((a == parameter[idx].b) && (b == parameter[idx].a))) 5116 { 5117 par = ¶meter[idx]; 5118 return par; 5119 } 5120 5121 return nullptr; 5122 } 5123 GetParameter3Atom(int a,int b,int c,std::vector<OBFFParameter> & parameter)5124 OBFFParameter* OBForceFieldMMFF94::GetParameter3Atom(int a, int b, int c, std::vector<OBFFParameter> ¶meter) 5125 { 5126 OBFFParameter *par; 5127 5128 for (unsigned int idx=0; idx < parameter.size(); idx++) 5129 if (((a == parameter[idx].a) && (b == parameter[idx].b) && (c == parameter[idx].c)) || 5130 ((a == parameter[idx].c) && (b == parameter[idx].b) && (c == parameter[idx].a))) 5131 { 5132 par = ¶meter[idx]; 5133 return par; 5134 } 5135 5136 return nullptr; 5137 } 5138 GetTypedParameter2Atom(int ffclass,int a,int b,std::vector<OBFFParameter> & parameter)5139 OBFFParameter* OBForceFieldMMFF94::GetTypedParameter2Atom(int ffclass, int a, int b, std::vector<OBFFParameter> ¶meter) 5140 { 5141 OBFFParameter *par; 5142 5143 for (unsigned int idx=0; idx < parameter.size(); idx++) 5144 if (((a == parameter[idx].a) && (b == parameter[idx].b) && (ffclass == parameter[idx]._ipar[0])) || 5145 ((a == parameter[idx].b) && (b == parameter[idx].a) && (ffclass == parameter[idx]._ipar[0]))) 5146 { 5147 par = ¶meter[idx]; 5148 return par; 5149 } 5150 5151 return nullptr; 5152 } 5153 GetTypedParameter3Atom(int ffclass,int a,int b,int c,std::vector<OBFFParameter> & parameter)5154 OBFFParameter* OBForceFieldMMFF94::GetTypedParameter3Atom(int ffclass, int a, int b, int c, std::vector<OBFFParameter> ¶meter) 5155 { 5156 OBFFParameter *par; 5157 5158 for (unsigned int idx=0; idx < parameter.size(); idx++) 5159 if (((a == parameter[idx].a) && (b == parameter[idx].b) && (c == parameter[idx].c) && (ffclass == parameter[idx]._ipar[0])) || 5160 ((a == parameter[idx].c) && (b == parameter[idx].b) && (c == parameter[idx].a) && (ffclass == parameter[idx]._ipar[0])) ) 5161 { 5162 par = ¶meter[idx]; 5163 return par; 5164 } 5165 5166 return nullptr; 5167 } 5168 GetTypedParameter4Atom(int ffclass,int a,int b,int c,int d,std::vector<OBFFParameter> & parameter)5169 OBFFParameter* OBForceFieldMMFF94::GetTypedParameter4Atom(int ffclass, int a, int b, int c, int d, std::vector<OBFFParameter> ¶meter) 5170 { 5171 OBFFParameter *par; 5172 5173 for (unsigned int idx=0; idx < parameter.size(); idx++) 5174 if (((a == parameter[idx].a) && (b == parameter[idx].b) && (c == parameter[idx].c) && 5175 (d == parameter[idx].d) && (ffclass == parameter[idx]._ipar[0])) 5176 /* || ((a == parameter[idx].d) && (b == parameter[idx].c) && (c == parameter[idx].b) && 5177 (d == parameter[idx].a) && (ffclass == parameter[idx]._ipar[0])) */ ) 5178 { 5179 par = ¶meter[idx]; 5180 return par; 5181 } 5182 5183 return nullptr; 5184 } 5185 5186 } // end namespace OpenBabel 5187 5188 //! \file forcefieldmmff94.cpp 5189 //! \brief MMFF94 force field 5190