Chemical Data Processing Library C++ API - Version 1.0.0
MMFF94GradientFunctions.hpp
Go to the documentation of this file.
1 /*
2  * MMFF94GradientFunctions.hpp
3  *
4  * This file is part of the Chemical Data Processing Toolkit
5  *
6  * Copyright (C) 2003 Thomas Seidel <thomas.seidel@univie.ac.at>
7  *
8  * This library is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU Lesser General Public
10  * License as published by the Free Software Foundation; either
11  * version 2 of the License, or (at your option) any later version.
12  *
13  * This library 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 GNU
16  * Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public License
19  * along with this library; see the file COPYING. If not, write to
20  * the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
21  * Boston, MA 02111-1307, USA.
22  */
23 
29 #ifndef CDPL_FORCEFIELD_MMFF94GRADIENTFUNCTIONS_HPP
30 #define CDPL_FORCEFIELD_MMFF94GRADIENTFUNCTIONS_HPP
31 
32 #include <algorithm>
33 
34 #include "CDPL/ForceField/MMFF94BondStretchingInteraction.hpp"
35 #include "CDPL/ForceField/MMFF94AngleBendingInteraction.hpp"
36 #include "CDPL/ForceField/MMFF94StretchBendInteraction.hpp"
37 #include "CDPL/ForceField/MMFF94OutOfPlaneBendingInteraction.hpp"
38 #include "CDPL/ForceField/MMFF94TorsionInteraction.hpp"
39 #include "CDPL/ForceField/MMFF94VanDerWaalsInteraction.hpp"
40 #include "CDPL/ForceField/MMFF94ElectrostaticInteraction.hpp"
41 #include "CDPL/ForceField/UtilityFunctions.hpp"
42 
43 
44 namespace CDPL
45 {
46 
47  namespace ForceField
48  {
49 
50  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
51  ValueType calcMMFF94BondStretchingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
52 
53  template <typename ValueType, typename CoordsArray, typename GradVector>
54  ValueType calcMMFF94BondStretchingGradient(const MMFF94BondStretchingInteraction& iaction, const CoordsArray& coords, GradVector& grad);
55 
94  template <typename ValueType, typename CoordsVec, typename GradVec>
95  ValueType calcMMFF94BondStretchingGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
96  const ValueType& force_const, const ValueType& ref_length);
97 
98 
99  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
100  ValueType calcMMFF94AngleBendingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
101 
102  template <typename ValueType, typename CoordsArray, typename GradVector>
103  ValueType calcMMFF94AngleBendingGradient(const MMFF94AngleBendingInteraction& iaction, const CoordsArray& coords, GradVector& grad);
104 
154  template <typename ValueType, typename CoordsVec, typename GradVec>
155  ValueType calcMMFF94AngleBendingGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
156  GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad, bool linear,
157  const ValueType& force_const, const ValueType& ref_angle);
158 
159 
160  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
161  ValueType calcMMFF94StretchBendGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
162 
163  template <typename ValueType, typename CoordsArray, typename GradVector>
164  ValueType calcMMFF94StretchBendGradient(const MMFF94StretchBendInteraction& iaction, const CoordsArray& coords, GradVector& grad);
165 
215  template <typename ValueType, typename CoordsVec, typename GradVec>
216  ValueType
217  calcMMFF94StretchBendGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos, GradVec& term_atom1_grad,
218  GradVec& ctr_atom_grad, GradVec& term_atom2_grad, const ValueType& ijk_force_const, const ValueType& kji_force_const,
219  const ValueType& ref_angle, const ValueType& ref_length1, const ValueType& ref_length2);
220 
221 
222  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
223  ValueType calcMMFF94OutOfPlaneBendingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
224 
225  template <typename ValueType, typename CoordsArray, typename GradVector>
226  ValueType calcMMFF94OutOfPlaneBendingGradient(const MMFF94OutOfPlaneBendingInteraction& iaction, const CoordsArray& coords, GradVector& grad);
227 
265  template <typename ValueType, typename CoordsVec, typename GradVec>
266  ValueType calcMMFF94OutOfPlaneBendingGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
267  const CoordsVec& oop_atom_pos, GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad,
268  GradVec& oop_atom_grad, const ValueType& force_const);
269 
270 
271  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
272  ValueType calcMMFF94TorsionGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
273 
274  template <typename ValueType, typename CoordsArray, typename GradVector>
275  ValueType calcMMFF94TorsionGradient(const MMFF94TorsionInteraction& iaction, const CoordsArray& coords, GradVector& grad);
276 
316  template <typename ValueType, typename CoordsVec, typename GradVec>
317  ValueType calcMMFF94TorsionGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos, const CoordsVec& ctr_atom2_pos,
318  const CoordsVec& term_atom2_pos, GradVec& term_atom1_grad, GradVec& ctr_atom1_grad, GradVec& ctr_atom2_grad,
319  GradVec& term_atom2_grad, const ValueType& tor_param1, const ValueType& tor_param2, const ValueType& tor_param3);
320 
321 
322  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
323  ValueType calcMMFF94ElectrostaticGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
324 
325  template <typename ValueType, typename CoordsArray, typename GradVector>
326  ValueType calcMMFF94ElectrostaticGradient(const MMFF94ElectrostaticInteraction& iaction, const CoordsArray& coords, GradVector& grad);
327 
369  template <typename ValueType, typename CoordsVec, typename GradVec>
370  ValueType calcMMFF94ElectrostaticGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
371  const ValueType& atom1_chg, const ValueType& atom2_chg, const ValueType& scale_fact,
372  const ValueType& de_const, const ValueType& dist_expo);
373 
374 
375  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
376  ValueType calcMMFF94VanDerWaalsGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
377 
378  template <typename ValueType, typename CoordsArray, typename GradVector>
379  ValueType calcMMFF94VanDerWaalsGradient(const MMFF94VanDerWaalsInteraction& iaction, const CoordsArray& coords, GradVector& grad);
380  /*
381  * dEvdwij/dVi = dEvdwij/dRij * dRij/dVi
382  * dEvdwij/dRij = -R*IJ^7 * eIJ / (Rij + 0.07 * R*IJ)^8 / (Rij^7 + 0.12 * R*IJ^7)^2 *
383  * (-22.48094067 * Rij^14 + 19.78322779 * Rij^7 * R*IJ^7 + 0.8812528743 * Rij^6 * R*IJ^8 + 1.186993667 * R*IJ^14)
384  */
385 
436  template <typename ValueType, typename CoordsVec, typename GradVec>
437  ValueType calcMMFF94VanDerWaalsGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
438  const ValueType& e_IJ, const ValueType& r_IJ, const ValueType& r_IJ_7);
439  } // namespace ForceField
440 } // namespace CDPL
441 
442 
443 // Implementation
444 // \cond DOC_IMPL_DETAILS
445 
446 namespace CDPL
447 {
448 
449  namespace ForceField
450  {
451 
452  namespace Detail
453  {
454 
455  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector, typename FuncType>
456  ValueType calcMMFF94InteractionGradient(Iter& beg, const Iter& end, const CoordsArray& coords, GradVector& grad, const FuncType& func)
457  {
458  ValueType e = ValueType();
459 
460  for (; beg != end; ++beg)
461  e += func(*beg, coords, grad);
462 
463  return e;
464  }
465  } // namespace Detail
466  } // namespace ForceField
467 } // namespace CDPL
468 
469 
470 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
471 ValueType CDPL::ForceField::calcMMFF94BondStretchingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
472 {
473  return Detail::calcMMFF94InteractionGradient<ValueType>(beg, end, coords, grad,
474  static_cast<ValueType (*)(const MMFF94BondStretchingInteraction&, const CoordsArray&, GradVector&)>(
475  &calcMMFF94BondStretchingGradient<ValueType, CoordsArray, GradVector>));
476 }
477 
478 template <typename ValueType, typename CoordsArray, typename GradVector>
479 ValueType CDPL::ForceField::calcMMFF94BondStretchingGradient(const MMFF94BondStretchingInteraction& iaction, const CoordsArray& coords, GradVector& grad)
480 {
481  return calcMMFF94BondStretchingGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
482  grad[iaction.getAtom2Index()], iaction.getForceConstant(), iaction.getReferenceLength());
483 }
484 
485 template <typename ValueType, typename CoordsVec, typename GradVec>
486 ValueType CDPL::ForceField::calcMMFF94BondStretchingGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
487  const ValueType& force_const, const ValueType& ref_length)
488 {
489  ValueType dist_atom1_grad[3];
490  ValueType dist_atom2_grad[3];
491 
492  ValueType dr_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad) - ref_length;
493  ValueType dr_ij_2 = dr_ij * dr_ij;
494 
495  ValueType grad_fact = (ValueType(167.92125 * 4) * dr_ij_2 * dr_ij - ValueType(215.89875 * 2) * dr_ij_2 + ValueType(143.9325) * dr_ij) * force_const;
496 
497  Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
498  Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
499 
500  ValueType e_b = ValueType(143.9325 * 0.5) * force_const * dr_ij_2 * (1 - 2 * dr_ij + 28 * dr_ij_2 / 12);
501 
502  return e_b;
503 }
504 
505 
506 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
507 ValueType CDPL::ForceField::calcMMFF94AngleBendingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
508 {
509  return Detail::calcMMFF94InteractionGradient<ValueType>(beg, end, coords, grad,
510  static_cast<ValueType (*)(const MMFF94AngleBendingInteraction&, const CoordsArray&, GradVector&)>(
511  &calcMMFF94AngleBendingGradient<ValueType, CoordsArray, GradVector>));
512 }
513 
514 template <typename ValueType, typename CoordsArray, typename GradVector>
515 ValueType CDPL::ForceField::calcMMFF94AngleBendingGradient(const MMFF94AngleBendingInteraction& iaction, const CoordsArray& coords, GradVector& grad)
516 {
517  return calcMMFF94AngleBendingGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
518  coords[iaction.getTerminalAtom2Index()], grad[iaction.getTerminalAtom1Index()],
519  grad[iaction.getCenterAtomIndex()], grad[iaction.getTerminalAtom2Index()], iaction.isLinearAngle(),
520  iaction.getForceConstant(), iaction.getReferenceAngle());
521 }
522 
523 template <typename ValueType, typename CoordsVec, typename GradVec>
524 ValueType CDPL::ForceField::calcMMFF94AngleBendingGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
525  GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad, bool linear,
526  const ValueType& force_const, const ValueType& ref_angle)
527 {
528  ValueType ac_term1_grad[3];
529  ValueType ac_ctr_grad[3];
530  ValueType ac_term2_grad[3];
531 
532  ValueType a_ijk_cos = calcBondAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, ac_term1_grad, ac_ctr_grad, ac_term2_grad);
533  ValueType grad_fact = ValueType(1);
534  ValueType e_a = ValueType(0);
535 
536  if (linear) {
537  grad_fact = ValueType(143.9325) * force_const;
538  e_a = ValueType(143.9325) * force_const * (1 + a_ijk_cos);
539 
540  } else {
541  ValueType a_ijk_cos_2 = a_ijk_cos * a_ijk_cos;
542  ValueType a_ijk = std::acos(a_ijk_cos);
543  ValueType div = std::sqrt(1 - a_ijk_cos_2);
544 
545  if (div < ValueType(0.0000001))
546  div = ValueType(0.0000001);
547 
548  grad_fact = force_const / div *
549  (a_ijk * (ValueType(86.58992538) * a_ijk - ValueType(143.9313616)) -
550  ref_angle * (ValueType(3.022558594) * a_ijk - ValueType(0.02637679965) * ref_angle - ValueType(2.512076157)));
551 
552  ValueType da_ijk = a_ijk * ValueType(180 / M_PI) - ref_angle;
553 
554  e_a = ValueType(0.043844 * 0.5) * force_const * da_ijk * da_ijk * (1 - ValueType(0.007) * da_ijk);
555  }
556 
557  Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
558  Detail::scaleAddVector(ac_ctr_grad, grad_fact, ctr_atom_grad);
559  Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
560 
561  return e_a;
562 }
563 
564 
565 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
566 ValueType CDPL::ForceField::calcMMFF94StretchBendGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
567 {
568  return Detail::calcMMFF94InteractionGradient<ValueType>(beg, end, coords, grad,
569  static_cast<ValueType (*)(const MMFF94StretchBendInteraction&, const CoordsArray&, GradVector&)>(
570  &calcMMFF94StretchBendGradient<ValueType, CoordsArray, GradVector>));
571 }
572 
573 template <typename ValueType, typename CoordsArray, typename GradVector>
574 ValueType CDPL::ForceField::calcMMFF94StretchBendGradient(const MMFF94StretchBendInteraction& iaction, const CoordsArray& coords, GradVector& grad)
575 {
576  return calcMMFF94StretchBendGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
577  coords[iaction.getTerminalAtom2Index()], grad[iaction.getTerminalAtom1Index()],
578  grad[iaction.getCenterAtomIndex()], grad[iaction.getTerminalAtom2Index()], iaction.getIJKForceConstant(),
579  iaction.getKJIForceConstant(), iaction.getReferenceAngle(), iaction.getReferenceLength1(),
580  iaction.getReferenceLength2());
581 }
582 
583 template <typename ValueType, typename CoordsVec, typename GradVec>
584 ValueType CDPL::ForceField::calcMMFF94StretchBendGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
585  GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad,
586  const ValueType& ijk_force_const, const ValueType& kji_force_const, const ValueType& ref_angle,
587  const ValueType& ref_length1, const ValueType& ref_length2)
588 {
589  ValueType ac_term1_grad[3];
590  ValueType ac_ctr_grad[3];
591  ValueType ac_term2_grad[3];
592 
593  ValueType dist_term1_grad[3];
594  ValueType dist_ctr_grad1[3];
595  ValueType dist_ctr_grad2[3];
596  ValueType dist_term2_grad[3];
597 
598  ValueType r_ij = calcDistanceDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, dist_term1_grad, dist_ctr_grad1);
599  ValueType r_kj = calcDistanceDerivatives<ValueType>(term_atom2_pos, ctr_atom_pos, dist_term2_grad, dist_ctr_grad2);
600  ValueType a_ijk_cos = calcBondAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, ac_term1_grad, ac_ctr_grad, ac_term2_grad);
601  ValueType a_ijk_cos_2 = a_ijk_cos * a_ijk_cos;
602  ValueType a_ijk = std::acos(a_ijk_cos);
603 
604  ValueType dr_ij = r_ij - ref_length1;
605  ValueType dr_kj = r_kj - ref_length2;
606  ValueType da_ijk = a_ijk * ValueType(180 / M_PI) - ref_angle;
607  ValueType div = std::sqrt(1 - a_ijk_cos_2);
608 
609  if (div < ValueType(0.0000001))
610  div = ValueType(0.0000001);
611 
612  ValueType a_ijk_grad_fact = ValueType(-180 * 2.5121 / M_PI) / div * (dr_ij * ijk_force_const + dr_kj * kji_force_const);
613 
614  ValueType r_ij_grad_fact = ValueType(2.5121) * da_ijk * ijk_force_const;
615  ValueType r_kj_grad_fact = ValueType(2.5121) * da_ijk * kji_force_const;
616 
617  Detail::scaleAddVector(dist_term1_grad, r_ij_grad_fact, term_atom1_grad);
618  Detail::scaleAddVector(ac_term1_grad, a_ijk_grad_fact, term_atom1_grad);
619 
620  Detail::scaleAddVector(dist_term2_grad, r_kj_grad_fact, term_atom2_grad);
621  Detail::scaleAddVector(ac_term2_grad, a_ijk_grad_fact, term_atom2_grad);
622 
623  Detail::scaleAddVector(dist_ctr_grad1, r_ij_grad_fact, ctr_atom_grad);
624  Detail::scaleAddVector(dist_ctr_grad2, r_kj_grad_fact, ctr_atom_grad);
625  Detail::scaleAddVector(ac_ctr_grad, a_ijk_grad_fact, ctr_atom_grad);
626 
627  ValueType e_ab = r_ij_grad_fact * dr_ij + r_kj_grad_fact * dr_kj;
628 
629  return e_ab;
630 }
631 
632 
633 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
634 ValueType CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
635 {
636  return Detail::calcMMFF94InteractionGradient<ValueType>(beg, end, coords, grad,
637  static_cast<ValueType (*)(const MMFF94OutOfPlaneBendingInteraction&, const CoordsArray&,
638  GradVector&)>(
639  &calcMMFF94OutOfPlaneBendingGradient<ValueType, CoordsArray, GradVector>));
640 }
641 
642 template <typename ValueType, typename CoordsArray, typename GradVector>
643 ValueType CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient(const MMFF94OutOfPlaneBendingInteraction& iaction, const CoordsArray& coords, GradVector& grad)
644 {
645  return calcMMFF94OutOfPlaneBendingGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
646  coords[iaction.getTerminalAtom2Index()], coords[iaction.getOutOfPlaneAtomIndex()],
647  grad[iaction.getTerminalAtom1Index()], grad[iaction.getCenterAtomIndex()],
648  grad[iaction.getTerminalAtom2Index()], grad[iaction.getOutOfPlaneAtomIndex()],
649  iaction.getForceConstant());
650 }
651 
652 template <typename ValueType, typename CoordsVec, typename GradVec>
653 ValueType CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
654  const CoordsVec& oop_atom_pos, GradVec& term_atom1_grad, GradVec& ctr_atom_grad,
655  GradVec& term_atom2_grad, GradVec& oop_atom_grad, const ValueType& force_const)
656 {
657  ValueType ac_term1_grad[3];
658  ValueType ac_term2_grad[3];
659  ValueType ac_ctr_grad[3];
660  ValueType ac_oop_grad[3];
661 
662  ValueType chi_ijkl_cos = calcOutOfPlaneAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, oop_atom_pos, ac_term1_grad,
663  ac_ctr_grad, ac_term2_grad, ac_oop_grad);
664  ValueType chi_ijkl = ValueType(M_PI * 0.5) - std::acos(chi_ijkl_cos);
665  ValueType div = std::sqrt(1 - chi_ijkl_cos * chi_ijkl_cos);
666 
667  if (div < ValueType(0.0000001))
668  div = ValueType(0.0000001);
669 
670  ValueType grad_fact = ValueType(0.043844 * 180 * 180) / div * ValueType(1 / (M_PI * M_PI)) * force_const * chi_ijkl;
671 
672  Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
673  Detail::scaleAddVector(ac_ctr_grad, grad_fact, ctr_atom_grad);
674  Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
675  Detail::scaleAddVector(ac_oop_grad, grad_fact, oop_atom_grad);
676 
677  chi_ijkl *= ValueType(180 / M_PI);
678 
679  ValueType e_oop = ValueType(0.5 * 0.043844) * force_const * chi_ijkl * chi_ijkl;
680 
681  return e_oop;
682 }
683 
684 
685 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
686 ValueType CDPL::ForceField::calcMMFF94TorsionGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
687 {
688  return Detail::calcMMFF94InteractionGradient<ValueType>(beg, end, coords, grad,
689  static_cast<ValueType (*)(const MMFF94TorsionInteraction&, const CoordsArray&, GradVector&)>(
690  &calcMMFF94TorsionGradient<ValueType, CoordsArray, GradVector>));
691 }
692 
693 template <typename ValueType, typename CoordsArray, typename GradVector>
694 ValueType CDPL::ForceField::calcMMFF94TorsionGradient(const MMFF94TorsionInteraction& iaction, const CoordsArray& coords, GradVector& grad)
695 {
696  return calcMMFF94TorsionGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtom1Index()],
697  coords[iaction.getCenterAtom2Index()], coords[iaction.getTerminalAtom2Index()],
698  grad[iaction.getTerminalAtom1Index()], grad[iaction.getCenterAtom1Index()], grad[iaction.getCenterAtom2Index()],
699  grad[iaction.getTerminalAtom2Index()], iaction.getTorsionParameter1(), iaction.getTorsionParameter2(),
700  iaction.getTorsionParameter3());
701 }
702 
703 template <typename ValueType, typename CoordsVec, typename GradVec>
704 ValueType
705 CDPL::ForceField::calcMMFF94TorsionGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos, const CoordsVec& ctr_atom2_pos,
706  const CoordsVec& term_atom2_pos, GradVec& term_atom1_grad, GradVec& ctr_atom1_grad, GradVec& ctr_atom2_grad,
707  GradVec& term_atom2_grad, const ValueType& tor_param1, const ValueType& tor_param2, const ValueType& tor_param3)
708 {
709  ValueType ac_term1_grad[3];
710  ValueType ac_ctr1_grad[3];
711  ValueType ac_ctr2_grad[3];
712  ValueType ac_term2_grad[3];
713 
714  ValueType phi_cos = calcDihedralAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom1_pos, ctr_atom2_pos, term_atom2_pos, ac_term1_grad, ac_ctr1_grad,
715  ac_ctr2_grad, ac_term2_grad);
716  ValueType phi = std::acos(phi_cos);
717  ValueType phi_cos_2 = phi_cos * phi_cos;
718  ValueType div = std::sqrt(1 - phi_cos_2);
719 
720  if (div < ValueType(0.0000001))
721  div = ValueType(0.0000001);
722 
723  ValueType grad_fact =
724  ValueType(-1) / div * (tor_param2 * std::sin(2 * phi) - ValueType(0.5) * tor_param1 * std::sin(phi) - ValueType(1.5) * tor_param3 * std::sin(3 * phi));
725 
726  Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
727  Detail::scaleAddVector(ac_ctr1_grad, grad_fact, ctr_atom1_grad);
728  Detail::scaleAddVector(ac_ctr2_grad, grad_fact, ctr_atom2_grad);
729  Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
730 
731  ValueType e_t = ValueType(0.5) * (tor_param1 * (1 + phi_cos) + tor_param2 * (1 - std::cos(2 * phi)) + tor_param3 * (1 + std::cos(3 * phi)));
732 
733  return e_t;
734 }
735 
736 
737 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
738 ValueType CDPL::ForceField::calcMMFF94ElectrostaticGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
739 {
740  return Detail::calcMMFF94InteractionGradient<ValueType>(beg, end, coords, grad,
741  static_cast<ValueType (*)(const MMFF94ElectrostaticInteraction&, const CoordsArray&, GradVector&)>(
742  &calcMMFF94ElectrostaticGradient<ValueType, CoordsArray, GradVector>));
743 }
744 
745 template <typename ValueType, typename CoordsArray, typename GradVector>
746 ValueType CDPL::ForceField::calcMMFF94ElectrostaticGradient(const MMFF94ElectrostaticInteraction& iaction, const CoordsArray& coords, GradVector& grad)
747 {
748  return calcMMFF94ElectrostaticGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
749  grad[iaction.getAtom2Index()], iaction.getAtom1Charge(), iaction.getAtom2Charge(),
750  iaction.getScalingFactor(), iaction.getDielectricConstant(), iaction.getDistanceExponent());
751 }
752 
753 template <typename ValueType, typename CoordsVec, typename GradVec>
754 ValueType CDPL::ForceField::calcMMFF94ElectrostaticGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
755  const ValueType& atom1_chg, const ValueType& atom2_chg, const ValueType& scale_fact,
756  const ValueType& de_const, const ValueType& dist_expo)
757 {
758  ValueType dist_atom1_grad[3];
759  ValueType dist_atom2_grad[3];
760 
761  ValueType r_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad);
762 
763  ValueType tmp1 = r_ij + ValueType(0.05);
764  ValueType tmp2 = std::pow(tmp1, dist_expo);
765  ValueType tmp3 = scale_fact * atom1_chg * atom2_chg / (de_const * tmp2);
766 
767  ValueType grad_fact = ValueType(-332.0716) * dist_expo * tmp3 / tmp1;
768 
769  Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
770  Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
771 
772  double e_q = ValueType(332.0716) * tmp3;
773 
774  return e_q;
775 }
776 
777 
778 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
779 ValueType CDPL::ForceField::calcMMFF94VanDerWaalsGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
780 {
781  return Detail::calcMMFF94InteractionGradient<ValueType>(beg, end, coords, grad,
782  static_cast<ValueType (*)(const MMFF94VanDerWaalsInteraction&, const CoordsArray&, GradVector&)>(
783  &calcMMFF94VanDerWaalsGradient<ValueType, CoordsArray, GradVector>));
784 }
785 
786 template <typename ValueType, typename CoordsArray, typename GradVector>
787 ValueType CDPL::ForceField::calcMMFF94VanDerWaalsGradient(const MMFF94VanDerWaalsInteraction& iaction, const CoordsArray& coords, GradVector& grad)
788 {
789  return calcMMFF94VanDerWaalsGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
790  grad[iaction.getAtom2Index()], iaction.getEIJ(), iaction.getRIJ(), iaction.getRIJPow7());
791 }
792 
793 template <typename ValueType, typename CoordsVec, typename GradVec>
794 ValueType CDPL::ForceField::calcMMFF94VanDerWaalsGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
795  const ValueType& e_IJ, const ValueType& r_IJ, const ValueType& r_IJ_7)
796 {
797  ValueType dist_atom1_grad[3];
798  ValueType dist_atom2_grad[3];
799 
800  ValueType r_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad);
801  ValueType r_ij_2 = r_ij * r_ij;
802  ValueType r_ij_6 = r_ij_2 * r_ij_2 * r_ij_2;
803  ValueType r_ij_7 = r_ij_6 * r_ij;
804 
805  ValueType tmp1 = r_ij + ValueType(0.07) * r_IJ;
806  ValueType tmp1_2 = tmp1 * tmp1;
807  ValueType tmp1_4 = tmp1_2 * tmp1_2;
808 
809  ValueType tmp2 = r_ij_7 + ValueType(0.12) * r_IJ_7;
810 
811  ValueType tmp3 = ValueType(1.07) * r_IJ / (r_ij + ValueType(0.07) * r_IJ);
812  ValueType tmp3_2 = tmp3 * tmp3;
813  ValueType tmp3_7 = tmp3_2 * tmp3_2 * tmp3_2 * tmp3;
814 
815  ValueType grad_fact = -r_IJ_7 * e_IJ / (tmp1_4 * tmp1_4 * tmp2 * tmp2) *
816  (ValueType(-22.48094067) * r_ij_7 * r_ij_7 + ValueType(19.78322779) * r_ij_7 * r_IJ_7 +
817  ValueType(0.8812528743) * r_ij_6 * r_IJ_7 * r_IJ + ValueType(1.186993667) * r_IJ_7 * r_IJ_7);
818 
819  Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
820  Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
821 
822  ValueType e_vdw = e_IJ * tmp3_7 * (ValueType(1.12) * r_IJ_7 / tmp2 - 2);
823 
824  return e_vdw;
825 }
826 
827 // \endcond
828 
829 #endif // CDPL_FORCEFIELD_MMFF94GRADIENTFUNCTIONS_HPP
MMFF94OutOfPlaneBendingInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94OutOfPlaneBendingInteraction.
MMFF94ElectrostaticInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94ElectrostaticInteraction.
CDPL::ForceField::calcMMFF94StretchBendGradient
ValueType calcMMFF94StretchBendGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
CDPL::ForceField::MMFF94BondStretchingInteraction
Definition: MMFF94BondStretchingInteraction.hpp:44
UtilityFunctions.hpp
Utility functions used in the calculation of force field energies and gradients.
MMFF94VanDerWaalsInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94VanDerWaalsInteraction.
MMFF94BondStretchingInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94BondStretchingInteraction.
CDPL::ForceField::MMFF94StretchBendInteraction
Definition: MMFF94StretchBendInteraction.hpp:44
MMFF94StretchBendInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94StretchBendInteraction.
CDPL::ForceField::calcMMFF94AngleBendingGradient
ValueType calcMMFF94AngleBendingGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
MMFF94TorsionInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94TorsionInteraction.
CDPL::ForceField::MMFF94OutOfPlaneBendingInteraction
Definition: MMFF94OutOfPlaneBendingInteraction.hpp:44
CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient
ValueType calcMMFF94OutOfPlaneBendingGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
CDPL::ForceField::MMFF94VanDerWaalsInteraction
Definition: MMFF94VanDerWaalsInteraction.hpp:45
CDPL::ForceField::calcMMFF94VanDerWaalsGradient
ValueType calcMMFF94VanDerWaalsGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
CDPL::ForceField::calcMMFF94BondStretchingGradient
ValueType calcMMFF94BondStretchingGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
CDPL::ForceField::MMFF94TorsionInteraction
Definition: MMFF94TorsionInteraction.hpp:44
CDPL
The namespace of the Chemical Data Processing Library.
CDPL::ForceField::calcMMFF94TorsionGradient
ValueType calcMMFF94TorsionGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
CDPL::ForceField::calcMMFF94ElectrostaticGradient
ValueType calcMMFF94ElectrostaticGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
CDPL::ForceField::MMFF94AngleBendingInteraction
Definition: MMFF94AngleBendingInteraction.hpp:44
MMFF94AngleBendingInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94AngleBendingInteraction.
CDPL::ForceField::MMFF94ElectrostaticInteraction
Definition: MMFF94ElectrostaticInteraction.hpp:44
Generated by doxygen 1.8.20