Chemical Data Processing Library C++ API - Version 1.2.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 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
447 ValueType CDPL::ForceField::calcMMFF94BondStretchingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
448 {
449  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
450  static_cast<ValueType (*)(const MMFF94BondStretchingInteraction&, const CoordsArray&, GradVector&)>(
451  &calcMMFF94BondStretchingGradient<ValueType, CoordsArray, GradVector>));
452 }
453 
454 template <typename ValueType, typename CoordsArray, typename GradVector>
455 ValueType CDPL::ForceField::calcMMFF94BondStretchingGradient(const MMFF94BondStretchingInteraction& iaction, const CoordsArray& coords, GradVector& grad)
456 {
457  return calcMMFF94BondStretchingGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
458  grad[iaction.getAtom2Index()], iaction.getForceConstant(), iaction.getReferenceLength());
459 }
460 
461 template <typename ValueType, typename CoordsVec, typename GradVec>
462 ValueType CDPL::ForceField::calcMMFF94BondStretchingGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
463  const ValueType& force_const, const ValueType& ref_length)
464 {
465  ValueType dist_atom1_grad[3];
466  ValueType dist_atom2_grad[3];
467 
468  ValueType dr_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad) - ref_length;
469  ValueType dr_ij_2 = dr_ij * dr_ij;
470 
471  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;
472 
473  Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
474  Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
475 
476  ValueType e_b = ValueType(143.9325 * 0.5) * force_const * dr_ij_2 * (1 - 2 * dr_ij + 28 * dr_ij_2 / 12);
477 
478  return e_b;
479 }
480 
481 
482 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
483 ValueType CDPL::ForceField::calcMMFF94AngleBendingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
484 {
485  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
486  static_cast<ValueType (*)(const MMFF94AngleBendingInteraction&, const CoordsArray&, GradVector&)>(
487  &calcMMFF94AngleBendingGradient<ValueType, CoordsArray, GradVector>));
488 }
489 
490 template <typename ValueType, typename CoordsArray, typename GradVector>
491 ValueType CDPL::ForceField::calcMMFF94AngleBendingGradient(const MMFF94AngleBendingInteraction& iaction, const CoordsArray& coords, GradVector& grad)
492 {
493  return calcMMFF94AngleBendingGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
494  coords[iaction.getTerminalAtom2Index()], grad[iaction.getTerminalAtom1Index()],
495  grad[iaction.getCenterAtomIndex()], grad[iaction.getTerminalAtom2Index()], iaction.isLinearAngle(),
496  iaction.getForceConstant(), iaction.getReferenceAngle());
497 }
498 
499 template <typename ValueType, typename CoordsVec, typename GradVec>
500 ValueType CDPL::ForceField::calcMMFF94AngleBendingGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
501  GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad, bool linear,
502  const ValueType& force_const, const ValueType& ref_angle)
503 {
504  ValueType ac_term1_grad[3];
505  ValueType ac_ctr_grad[3];
506  ValueType ac_term2_grad[3];
507 
508  ValueType a_ijk_cos = calcBondAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, ac_term1_grad, ac_ctr_grad, ac_term2_grad);
509  ValueType grad_fact = ValueType(1);
510  ValueType e_a = ValueType(0);
511 
512  if (linear) {
513  grad_fact = ValueType(143.9325) * force_const;
514  e_a = ValueType(143.9325) * force_const * (1 + a_ijk_cos);
515 
516  } else {
517  ValueType a_ijk_cos_2 = a_ijk_cos * a_ijk_cos;
518  ValueType a_ijk = std::acos(a_ijk_cos);
519  ValueType div = std::sqrt(1 - a_ijk_cos_2);
520 
521  if (div < ValueType(0.0000001))
522  div = ValueType(0.0000001);
523 
524  grad_fact = force_const / div *
525  (a_ijk * (ValueType(86.58992538) * a_ijk - ValueType(143.9313616)) -
526  ref_angle * (ValueType(3.022558594) * a_ijk - ValueType(0.02637679965) * ref_angle - ValueType(2.512076157)));
527 
528  ValueType da_ijk = a_ijk * ValueType(180 / M_PI) - ref_angle;
529 
530  e_a = ValueType(0.043844 * 0.5) * force_const * da_ijk * da_ijk * (1 - ValueType(0.007) * da_ijk);
531  }
532 
533  Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
534  Detail::scaleAddVector(ac_ctr_grad, grad_fact, ctr_atom_grad);
535  Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
536 
537  return e_a;
538 }
539 
540 
541 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
542 ValueType CDPL::ForceField::calcMMFF94StretchBendGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
543 {
544  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
545  static_cast<ValueType (*)(const MMFF94StretchBendInteraction&, const CoordsArray&, GradVector&)>(
546  &calcMMFF94StretchBendGradient<ValueType, CoordsArray, GradVector>));
547 }
548 
549 template <typename ValueType, typename CoordsArray, typename GradVector>
550 ValueType CDPL::ForceField::calcMMFF94StretchBendGradient(const MMFF94StretchBendInteraction& iaction, const CoordsArray& coords, GradVector& grad)
551 {
552  return calcMMFF94StretchBendGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
553  coords[iaction.getTerminalAtom2Index()], grad[iaction.getTerminalAtom1Index()],
554  grad[iaction.getCenterAtomIndex()], grad[iaction.getTerminalAtom2Index()], iaction.getIJKForceConstant(),
555  iaction.getKJIForceConstant(), iaction.getReferenceAngle(), iaction.getReferenceLength1(),
556  iaction.getReferenceLength2());
557 }
558 
559 template <typename ValueType, typename CoordsVec, typename GradVec>
560 ValueType CDPL::ForceField::calcMMFF94StretchBendGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
561  GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad,
562  const ValueType& ijk_force_const, const ValueType& kji_force_const, const ValueType& ref_angle,
563  const ValueType& ref_length1, const ValueType& ref_length2)
564 {
565  ValueType ac_term1_grad[3];
566  ValueType ac_ctr_grad[3];
567  ValueType ac_term2_grad[3];
568 
569  ValueType dist_term1_grad[3];
570  ValueType dist_ctr_grad1[3];
571  ValueType dist_ctr_grad2[3];
572  ValueType dist_term2_grad[3];
573 
574  ValueType r_ij = calcDistanceDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, dist_term1_grad, dist_ctr_grad1);
575  ValueType r_kj = calcDistanceDerivatives<ValueType>(term_atom2_pos, ctr_atom_pos, dist_term2_grad, dist_ctr_grad2);
576  ValueType a_ijk_cos = calcBondAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, ac_term1_grad, ac_ctr_grad, ac_term2_grad);
577  ValueType a_ijk_cos_2 = a_ijk_cos * a_ijk_cos;
578  ValueType a_ijk = std::acos(a_ijk_cos);
579 
580  ValueType dr_ij = r_ij - ref_length1;
581  ValueType dr_kj = r_kj - ref_length2;
582  ValueType da_ijk = a_ijk * ValueType(180 / M_PI) - ref_angle;
583  ValueType div = std::sqrt(1 - a_ijk_cos_2);
584 
585  if (div < ValueType(0.0000001))
586  div = ValueType(0.0000001);
587 
588  ValueType a_ijk_grad_fact = ValueType(-180 * 2.5121 / M_PI) / div * (dr_ij * ijk_force_const + dr_kj * kji_force_const);
589 
590  ValueType r_ij_grad_fact = ValueType(2.5121) * da_ijk * ijk_force_const;
591  ValueType r_kj_grad_fact = ValueType(2.5121) * da_ijk * kji_force_const;
592 
593  Detail::scaleAddVector(dist_term1_grad, r_ij_grad_fact, term_atom1_grad);
594  Detail::scaleAddVector(ac_term1_grad, a_ijk_grad_fact, term_atom1_grad);
595 
596  Detail::scaleAddVector(dist_term2_grad, r_kj_grad_fact, term_atom2_grad);
597  Detail::scaleAddVector(ac_term2_grad, a_ijk_grad_fact, term_atom2_grad);
598 
599  Detail::scaleAddVector(dist_ctr_grad1, r_ij_grad_fact, ctr_atom_grad);
600  Detail::scaleAddVector(dist_ctr_grad2, r_kj_grad_fact, ctr_atom_grad);
601  Detail::scaleAddVector(ac_ctr_grad, a_ijk_grad_fact, ctr_atom_grad);
602 
603  ValueType e_ab = r_ij_grad_fact * dr_ij + r_kj_grad_fact * dr_kj;
604 
605  return e_ab;
606 }
607 
608 
609 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
610 ValueType CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
611 {
612  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
613  static_cast<ValueType (*)(const MMFF94OutOfPlaneBendingInteraction&, const CoordsArray&, GradVector&)>(
614  &calcMMFF94OutOfPlaneBendingGradient<ValueType, CoordsArray, GradVector>));
615 }
616 
617 template <typename ValueType, typename CoordsArray, typename GradVector>
618 ValueType CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient(const MMFF94OutOfPlaneBendingInteraction& iaction, const CoordsArray& coords, GradVector& grad)
619 {
620  return calcMMFF94OutOfPlaneBendingGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
621  coords[iaction.getTerminalAtom2Index()], coords[iaction.getOutOfPlaneAtomIndex()],
622  grad[iaction.getTerminalAtom1Index()], grad[iaction.getCenterAtomIndex()],
623  grad[iaction.getTerminalAtom2Index()], grad[iaction.getOutOfPlaneAtomIndex()],
624  iaction.getForceConstant());
625 }
626 
627 template <typename ValueType, typename CoordsVec, typename GradVec>
628 ValueType CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
629  const CoordsVec& oop_atom_pos, GradVec& term_atom1_grad, GradVec& ctr_atom_grad,
630  GradVec& term_atom2_grad, GradVec& oop_atom_grad, const ValueType& force_const)
631 {
632  ValueType ac_term1_grad[3];
633  ValueType ac_term2_grad[3];
634  ValueType ac_ctr_grad[3];
635  ValueType ac_oop_grad[3];
636 
637  ValueType chi_ijkl_cos = calcOutOfPlaneAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, oop_atom_pos, ac_term1_grad,
638  ac_ctr_grad, ac_term2_grad, ac_oop_grad);
639  ValueType chi_ijkl = ValueType(M_PI * 0.5) - std::acos(chi_ijkl_cos);
640  ValueType div = std::sqrt(1 - chi_ijkl_cos * chi_ijkl_cos);
641 
642  if (div < ValueType(0.0000001))
643  div = ValueType(0.0000001);
644 
645  ValueType grad_fact = ValueType(0.043844 * 180 * 180) / div * ValueType(1 / (M_PI * M_PI)) * force_const * chi_ijkl;
646 
647  Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
648  Detail::scaleAddVector(ac_ctr_grad, grad_fact, ctr_atom_grad);
649  Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
650  Detail::scaleAddVector(ac_oop_grad, grad_fact, oop_atom_grad);
651 
652  chi_ijkl *= ValueType(180 / M_PI);
653 
654  ValueType e_oop = ValueType(0.5 * 0.043844) * force_const * chi_ijkl * chi_ijkl;
655 
656  return e_oop;
657 }
658 
659 
660 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
661 ValueType CDPL::ForceField::calcMMFF94TorsionGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
662 {
663  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
664  static_cast<ValueType (*)(const MMFF94TorsionInteraction&, const CoordsArray&, GradVector&)>(
665  &calcMMFF94TorsionGradient<ValueType, CoordsArray, GradVector>));
666 }
667 
668 template <typename ValueType, typename CoordsArray, typename GradVector>
669 ValueType CDPL::ForceField::calcMMFF94TorsionGradient(const MMFF94TorsionInteraction& iaction, const CoordsArray& coords, GradVector& grad)
670 {
671  return calcMMFF94TorsionGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtom1Index()],
672  coords[iaction.getCenterAtom2Index()], coords[iaction.getTerminalAtom2Index()],
673  grad[iaction.getTerminalAtom1Index()], grad[iaction.getCenterAtom1Index()], grad[iaction.getCenterAtom2Index()],
674  grad[iaction.getTerminalAtom2Index()], iaction.getTorsionParameter1(), iaction.getTorsionParameter2(),
675  iaction.getTorsionParameter3());
676 }
677 
678 template <typename ValueType, typename CoordsVec, typename GradVec>
679 ValueType
680 CDPL::ForceField::calcMMFF94TorsionGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos, const CoordsVec& ctr_atom2_pos,
681  const CoordsVec& term_atom2_pos, GradVec& term_atom1_grad, GradVec& ctr_atom1_grad, GradVec& ctr_atom2_grad,
682  GradVec& term_atom2_grad, const ValueType& tor_param1, const ValueType& tor_param2, const ValueType& tor_param3)
683 {
684  ValueType ac_term1_grad[3];
685  ValueType ac_ctr1_grad[3];
686  ValueType ac_ctr2_grad[3];
687  ValueType ac_term2_grad[3];
688 
689  ValueType phi_cos = calcDihedralAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom1_pos, ctr_atom2_pos, term_atom2_pos, ac_term1_grad, ac_ctr1_grad,
690  ac_ctr2_grad, ac_term2_grad);
691  ValueType phi = std::acos(phi_cos);
692  ValueType phi_cos_2 = phi_cos * phi_cos;
693  ValueType div = std::sqrt(1 - phi_cos_2);
694 
695  if (div < ValueType(0.0000001))
696  div = ValueType(0.0000001);
697 
698  ValueType grad_fact =
699  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));
700 
701  Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
702  Detail::scaleAddVector(ac_ctr1_grad, grad_fact, ctr_atom1_grad);
703  Detail::scaleAddVector(ac_ctr2_grad, grad_fact, ctr_atom2_grad);
704  Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
705 
706  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)));
707 
708  return e_t;
709 }
710 
711 
712 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
713 ValueType CDPL::ForceField::calcMMFF94ElectrostaticGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
714 {
715  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
716  static_cast<ValueType (*)(const MMFF94ElectrostaticInteraction&, const CoordsArray&, GradVector&)>(
717  &calcMMFF94ElectrostaticGradient<ValueType, CoordsArray, GradVector>));
718 }
719 
720 template <typename ValueType, typename CoordsArray, typename GradVector>
721 ValueType CDPL::ForceField::calcMMFF94ElectrostaticGradient(const MMFF94ElectrostaticInteraction& iaction, const CoordsArray& coords, GradVector& grad)
722 {
723  return calcMMFF94ElectrostaticGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
724  grad[iaction.getAtom2Index()], iaction.getAtom1Charge(), iaction.getAtom2Charge(),
725  iaction.getScalingFactor(), iaction.getDielectricConstant(), iaction.getDistanceExponent());
726 }
727 
728 template <typename ValueType, typename CoordsVec, typename GradVec>
729 ValueType CDPL::ForceField::calcMMFF94ElectrostaticGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
730  const ValueType& atom1_chg, const ValueType& atom2_chg, const ValueType& scale_fact,
731  const ValueType& de_const, const ValueType& dist_expo)
732 {
733  ValueType dist_atom1_grad[3];
734  ValueType dist_atom2_grad[3];
735 
736  ValueType r_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad);
737 
738  ValueType tmp1 = r_ij + ValueType(0.05);
739  ValueType tmp2 = std::pow(tmp1, dist_expo);
740  ValueType tmp3 = scale_fact * atom1_chg * atom2_chg / (de_const * tmp2);
741 
742  ValueType grad_fact = ValueType(-332.0716) * dist_expo * tmp3 / tmp1;
743 
744  Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
745  Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
746 
747  double e_q = ValueType(332.0716) * tmp3;
748 
749  return e_q;
750 }
751 
752 
753 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
754 ValueType CDPL::ForceField::calcMMFF94VanDerWaalsGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
755 {
756  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
757  static_cast<ValueType (*)(const MMFF94VanDerWaalsInteraction&, const CoordsArray&, GradVector&)>(
758  &calcMMFF94VanDerWaalsGradient<ValueType, CoordsArray, GradVector>));
759 }
760 
761 template <typename ValueType, typename CoordsArray, typename GradVector>
762 ValueType CDPL::ForceField::calcMMFF94VanDerWaalsGradient(const MMFF94VanDerWaalsInteraction& iaction, const CoordsArray& coords, GradVector& grad)
763 {
764  return calcMMFF94VanDerWaalsGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
765  grad[iaction.getAtom2Index()], iaction.getEIJ(), iaction.getRIJ(), iaction.getRIJPow7());
766 }
767 
768 template <typename ValueType, typename CoordsVec, typename GradVec>
769 ValueType CDPL::ForceField::calcMMFF94VanDerWaalsGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
770  const ValueType& e_IJ, const ValueType& r_IJ, const ValueType& r_IJ_7)
771 {
772  ValueType dist_atom1_grad[3];
773  ValueType dist_atom2_grad[3];
774 
775  ValueType r_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad);
776  ValueType r_ij_2 = r_ij * r_ij;
777  ValueType r_ij_6 = r_ij_2 * r_ij_2 * r_ij_2;
778  ValueType r_ij_7 = r_ij_6 * r_ij;
779 
780  ValueType tmp1 = r_ij + ValueType(0.07) * r_IJ;
781  ValueType tmp1_2 = tmp1 * tmp1;
782  ValueType tmp1_4 = tmp1_2 * tmp1_2;
783 
784  ValueType tmp2 = r_ij_7 + ValueType(0.12) * r_IJ_7;
785 
786  ValueType tmp3 = ValueType(1.07) * r_IJ / (r_ij + ValueType(0.07) * r_IJ);
787  ValueType tmp3_2 = tmp3 * tmp3;
788  ValueType tmp3_7 = tmp3_2 * tmp3_2 * tmp3_2 * tmp3;
789 
790  ValueType grad_fact = -r_IJ_7 * e_IJ / (tmp1_4 * tmp1_4 * tmp2 * tmp2) *
791  (ValueType(-22.48094067) * r_ij_7 * r_ij_7 + ValueType(19.78322779) * r_ij_7 * r_IJ_7 +
792  ValueType(0.8812528743) * r_ij_6 * r_IJ_7 * r_IJ + ValueType(1.186993667) * r_IJ_7 * r_IJ_7);
793 
794  Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
795  Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
796 
797  ValueType e_vdw = e_IJ * tmp3_7 * (ValueType(1.12) * r_IJ_7 / tmp2 - 2);
798 
799  return e_vdw;
800 }
801 
802 // \endcond
803 
804 #endif // CDPL_FORCEFIELD_MMFF94GRADIENTFUNCTIONS_HPP
UtilityFunctions.hpp
Utility functions used in the calculation of force field energies and gradients.
MMFF94AngleBendingInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94AngleBendingInteraction.
MMFF94BondStretchingInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94BondStretchingInteraction.
MMFF94ElectrostaticInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94ElectrostaticInteraction.
MMFF94OutOfPlaneBendingInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94OutOfPlaneBendingInteraction.
MMFF94StretchBendInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94StretchBendInteraction.
MMFF94TorsionInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94TorsionInteraction.
MMFF94VanDerWaalsInteraction.hpp
Definition of the class CDPL::ForceField::MMFF94VanDerWaalsInteraction.
CDPL::ForceField::MMFF94AngleBendingInteraction
Definition: MMFF94AngleBendingInteraction.hpp:42
CDPL::ForceField::MMFF94BondStretchingInteraction
Definition: MMFF94BondStretchingInteraction.hpp:42
CDPL::ForceField::MMFF94ElectrostaticInteraction
Definition: MMFF94ElectrostaticInteraction.hpp:42
CDPL::ForceField::MMFF94OutOfPlaneBendingInteraction
Definition: MMFF94OutOfPlaneBendingInteraction.hpp:42
CDPL::ForceField::MMFF94StretchBendInteraction
Definition: MMFF94StretchBendInteraction.hpp:42
CDPL::ForceField::MMFF94TorsionInteraction
Definition: MMFF94TorsionInteraction.hpp:42
CDPL::ForceField::MMFF94VanDerWaalsInteraction
Definition: MMFF94VanDerWaalsInteraction.hpp:43
CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient
ValueType calcMMFF94OutOfPlaneBendingGradient(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::calcMMFF94VanDerWaalsGradient
ValueType calcMMFF94VanDerWaalsGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
CDPL::ForceField::calcMMFF94AngleBendingGradient
ValueType calcMMFF94AngleBendingGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
CDPL::ForceField::calcMMFF94StretchBendGradient
ValueType calcMMFF94StretchBendGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
CDPL::ForceField::calcMMFF94TorsionGradient
ValueType calcMMFF94TorsionGradient(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
The namespace of the Chemical Data Processing Library.
Generated by doxygen 1.9.1