Chemical Data Processing Library C++ API - Version 1.4.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 
62  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
63  ValueType calcMMFF94BondStretchingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
64 
75  template <typename ValueType, typename CoordsArray, typename GradVector>
76  ValueType calcMMFF94BondStretchingGradient(const MMFF94BondStretchingInteraction& iaction, const CoordsArray& coords, GradVector& grad);
77 
117  template <typename ValueType, typename CoordsVec, typename GradVec>
118  ValueType calcMMFF94BondStretchingGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
119  const ValueType& force_const, const ValueType& ref_length);
120 
133  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
134  ValueType calcMMFF94AngleBendingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
135 
146  template <typename ValueType, typename CoordsArray, typename GradVector>
147  ValueType calcMMFF94AngleBendingGradient(const MMFF94AngleBendingInteraction& iaction, const CoordsArray& coords, GradVector& grad);
148 
199  template <typename ValueType, typename CoordsVec, typename GradVec>
200  ValueType calcMMFF94AngleBendingGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
201  GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad, bool linear,
202  const ValueType& force_const, const ValueType& ref_angle);
203 
216  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
217  ValueType calcMMFF94StretchBendGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
218 
229  template <typename ValueType, typename CoordsArray, typename GradVector>
230  ValueType calcMMFF94StretchBendGradient(const MMFF94StretchBendInteraction& iaction, const CoordsArray& coords, GradVector& grad);
231 
282  template <typename ValueType, typename CoordsVec, typename GradVec>
283  ValueType
284  calcMMFF94StretchBendGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos, GradVec& term_atom1_grad,
285  GradVec& ctr_atom_grad, GradVec& term_atom2_grad, const ValueType& ijk_force_const, const ValueType& kji_force_const,
286  const ValueType& ref_angle, const ValueType& ref_length1, const ValueType& ref_length2);
287 
300  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
301  ValueType calcMMFF94OutOfPlaneBendingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
302 
313  template <typename ValueType, typename CoordsArray, typename GradVector>
314  ValueType calcMMFF94OutOfPlaneBendingGradient(const MMFF94OutOfPlaneBendingInteraction& iaction, const CoordsArray& coords, GradVector& grad);
315 
354  template <typename ValueType, typename CoordsVec, typename GradVec>
355  ValueType calcMMFF94OutOfPlaneBendingGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
356  const CoordsVec& oop_atom_pos, GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad,
357  GradVec& oop_atom_grad, const ValueType& force_const);
358 
371  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
372  ValueType calcMMFF94TorsionGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
373 
384  template <typename ValueType, typename CoordsArray, typename GradVector>
385  ValueType calcMMFF94TorsionGradient(const MMFF94TorsionInteraction& iaction, const CoordsArray& coords, GradVector& grad);
386 
427  template <typename ValueType, typename CoordsVec, typename GradVec>
428  ValueType calcMMFF94TorsionGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos, const CoordsVec& ctr_atom2_pos,
429  const CoordsVec& term_atom2_pos, GradVec& term_atom1_grad, GradVec& ctr_atom1_grad, GradVec& ctr_atom2_grad,
430  GradVec& term_atom2_grad, const ValueType& tor_param1, const ValueType& tor_param2, const ValueType& tor_param3);
431 
444  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
445  ValueType calcMMFF94ElectrostaticGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
446 
457  template <typename ValueType, typename CoordsArray, typename GradVector>
458  ValueType calcMMFF94ElectrostaticGradient(const MMFF94ElectrostaticInteraction& iaction, const CoordsArray& coords, GradVector& grad);
459 
502  template <typename ValueType, typename CoordsVec, typename GradVec>
503  ValueType calcMMFF94ElectrostaticGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
504  const ValueType& atom1_chg, const ValueType& atom2_chg, const ValueType& scale_fact,
505  const ValueType& de_const, const ValueType& dist_expo);
506 
519  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
520  ValueType calcMMFF94VanDerWaalsGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad);
521 
532  template <typename ValueType, typename CoordsArray, typename GradVector>
533  ValueType calcMMFF94VanDerWaalsGradient(const MMFF94VanDerWaalsInteraction& iaction, const CoordsArray& coords, GradVector& grad);
534 
586  template <typename ValueType, typename CoordsVec, typename GradVec>
587  ValueType calcMMFF94VanDerWaalsGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
588  const ValueType& e_IJ, const ValueType& r_IJ, const ValueType& r_IJ_7);
589  } // namespace ForceField
590 } // namespace CDPL
591 
592 
593 // Implementation
594 // \cond DOC_IMPL_DETAILS
595 
596 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
597 ValueType CDPL::ForceField::calcMMFF94BondStretchingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
598 {
599  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
600  static_cast<ValueType (*)(const MMFF94BondStretchingInteraction&, const CoordsArray&, GradVector&)>(
601  &calcMMFF94BondStretchingGradient<ValueType, CoordsArray, GradVector>));
602 }
603 
604 template <typename ValueType, typename CoordsArray, typename GradVector>
605 ValueType CDPL::ForceField::calcMMFF94BondStretchingGradient(const MMFF94BondStretchingInteraction& iaction, const CoordsArray& coords, GradVector& grad)
606 {
607  return calcMMFF94BondStretchingGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
608  grad[iaction.getAtom2Index()], iaction.getForceConstant(), iaction.getReferenceLength());
609 }
610 
611 template <typename ValueType, typename CoordsVec, typename GradVec>
612 ValueType CDPL::ForceField::calcMMFF94BondStretchingGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
613  const ValueType& force_const, const ValueType& ref_length)
614 {
615  ValueType dist_atom1_grad[3];
616  ValueType dist_atom2_grad[3];
617 
618  ValueType dr_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad) - ref_length;
619  ValueType dr_ij_2 = dr_ij * dr_ij;
620 
621  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;
622 
623  Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
624  Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
625 
626  ValueType e_b = ValueType(143.9325 * 0.5) * force_const * dr_ij_2 * (1 - 2 * dr_ij + 28 * dr_ij_2 / 12);
627 
628  return e_b;
629 }
630 
631 
632 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
633 ValueType CDPL::ForceField::calcMMFF94AngleBendingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
634 {
635  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
636  static_cast<ValueType (*)(const MMFF94AngleBendingInteraction&, const CoordsArray&, GradVector&)>(
637  &calcMMFF94AngleBendingGradient<ValueType, CoordsArray, GradVector>));
638 }
639 
640 template <typename ValueType, typename CoordsArray, typename GradVector>
641 ValueType CDPL::ForceField::calcMMFF94AngleBendingGradient(const MMFF94AngleBendingInteraction& iaction, const CoordsArray& coords, GradVector& grad)
642 {
643  return calcMMFF94AngleBendingGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
644  coords[iaction.getTerminalAtom2Index()], grad[iaction.getTerminalAtom1Index()],
645  grad[iaction.getCenterAtomIndex()], grad[iaction.getTerminalAtom2Index()], iaction.isLinearAngle(),
646  iaction.getForceConstant(), iaction.getReferenceAngle());
647 }
648 
649 template <typename ValueType, typename CoordsVec, typename GradVec>
650 ValueType CDPL::ForceField::calcMMFF94AngleBendingGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
651  GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad, bool linear,
652  const ValueType& force_const, const ValueType& ref_angle)
653 {
654  ValueType ac_term1_grad[3];
655  ValueType ac_ctr_grad[3];
656  ValueType ac_term2_grad[3];
657 
658  ValueType a_ijk_cos = calcBondAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, ac_term1_grad, ac_ctr_grad, ac_term2_grad);
659  ValueType grad_fact = ValueType(1);
660  ValueType e_a = ValueType(0);
661 
662  if (linear) {
663  grad_fact = ValueType(143.9325) * force_const;
664  e_a = ValueType(143.9325) * force_const * (1 + a_ijk_cos);
665 
666  } else {
667  ValueType a_ijk_cos_2 = a_ijk_cos * a_ijk_cos;
668  ValueType a_ijk = std::acos(a_ijk_cos);
669  ValueType div = std::sqrt(1 - a_ijk_cos_2);
670 
671  if (div < ValueType(0.0000001))
672  div = ValueType(0.0000001);
673 
674  grad_fact = force_const / div *
675  (a_ijk * (ValueType(86.58992538) * a_ijk - ValueType(143.9313616)) -
676  ref_angle * (ValueType(3.022558594) * a_ijk - ValueType(0.02637679965) * ref_angle - ValueType(2.512076157)));
677 
678  ValueType da_ijk = a_ijk * ValueType(180 / M_PI) - ref_angle;
679 
680  e_a = ValueType(0.043844 * 0.5) * force_const * da_ijk * da_ijk * (1 - ValueType(0.007) * da_ijk);
681  }
682 
683  Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
684  Detail::scaleAddVector(ac_ctr_grad, grad_fact, ctr_atom_grad);
685  Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
686 
687  return e_a;
688 }
689 
690 
691 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
692 ValueType CDPL::ForceField::calcMMFF94StretchBendGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
693 {
694  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
695  static_cast<ValueType (*)(const MMFF94StretchBendInteraction&, const CoordsArray&, GradVector&)>(
696  &calcMMFF94StretchBendGradient<ValueType, CoordsArray, GradVector>));
697 }
698 
699 template <typename ValueType, typename CoordsArray, typename GradVector>
700 ValueType CDPL::ForceField::calcMMFF94StretchBendGradient(const MMFF94StretchBendInteraction& iaction, const CoordsArray& coords, GradVector& grad)
701 {
702  return calcMMFF94StretchBendGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
703  coords[iaction.getTerminalAtom2Index()], grad[iaction.getTerminalAtom1Index()],
704  grad[iaction.getCenterAtomIndex()], grad[iaction.getTerminalAtom2Index()], iaction.getIJKForceConstant(),
705  iaction.getKJIForceConstant(), iaction.getReferenceAngle(), iaction.getReferenceLength1(),
706  iaction.getReferenceLength2());
707 }
708 
709 template <typename ValueType, typename CoordsVec, typename GradVec>
710 ValueType CDPL::ForceField::calcMMFF94StretchBendGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
711  GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad,
712  const ValueType& ijk_force_const, const ValueType& kji_force_const, const ValueType& ref_angle,
713  const ValueType& ref_length1, const ValueType& ref_length2)
714 {
715  ValueType ac_term1_grad[3];
716  ValueType ac_ctr_grad[3];
717  ValueType ac_term2_grad[3];
718 
719  ValueType dist_term1_grad[3];
720  ValueType dist_ctr_grad1[3];
721  ValueType dist_ctr_grad2[3];
722  ValueType dist_term2_grad[3];
723 
724  ValueType r_ij = calcDistanceDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, dist_term1_grad, dist_ctr_grad1);
725  ValueType r_kj = calcDistanceDerivatives<ValueType>(term_atom2_pos, ctr_atom_pos, dist_term2_grad, dist_ctr_grad2);
726  ValueType a_ijk_cos = calcBondAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, ac_term1_grad, ac_ctr_grad, ac_term2_grad);
727  ValueType a_ijk_cos_2 = a_ijk_cos * a_ijk_cos;
728  ValueType a_ijk = std::acos(a_ijk_cos);
729 
730  ValueType dr_ij = r_ij - ref_length1;
731  ValueType dr_kj = r_kj - ref_length2;
732  ValueType da_ijk = a_ijk * ValueType(180 / M_PI) - ref_angle;
733  ValueType div = std::sqrt(1 - a_ijk_cos_2);
734 
735  if (div < ValueType(0.0000001))
736  div = ValueType(0.0000001);
737 
738  ValueType a_ijk_grad_fact = ValueType(-180 * 2.5121 / M_PI) / div * (dr_ij * ijk_force_const + dr_kj * kji_force_const);
739 
740  ValueType r_ij_grad_fact = ValueType(2.5121) * da_ijk * ijk_force_const;
741  ValueType r_kj_grad_fact = ValueType(2.5121) * da_ijk * kji_force_const;
742 
743  Detail::scaleAddVector(dist_term1_grad, r_ij_grad_fact, term_atom1_grad);
744  Detail::scaleAddVector(ac_term1_grad, a_ijk_grad_fact, term_atom1_grad);
745 
746  Detail::scaleAddVector(dist_term2_grad, r_kj_grad_fact, term_atom2_grad);
747  Detail::scaleAddVector(ac_term2_grad, a_ijk_grad_fact, term_atom2_grad);
748 
749  Detail::scaleAddVector(dist_ctr_grad1, r_ij_grad_fact, ctr_atom_grad);
750  Detail::scaleAddVector(dist_ctr_grad2, r_kj_grad_fact, ctr_atom_grad);
751  Detail::scaleAddVector(ac_ctr_grad, a_ijk_grad_fact, ctr_atom_grad);
752 
753  ValueType e_ab = r_ij_grad_fact * dr_ij + r_kj_grad_fact * dr_kj;
754 
755  return e_ab;
756 }
757 
758 
759 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
760 ValueType CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
761 {
762  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
763  static_cast<ValueType (*)(const MMFF94OutOfPlaneBendingInteraction&, const CoordsArray&, GradVector&)>(
764  &calcMMFF94OutOfPlaneBendingGradient<ValueType, CoordsArray, GradVector>));
765 }
766 
767 template <typename ValueType, typename CoordsArray, typename GradVector>
768 ValueType CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient(const MMFF94OutOfPlaneBendingInteraction& iaction, const CoordsArray& coords, GradVector& grad)
769 {
770  return calcMMFF94OutOfPlaneBendingGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
771  coords[iaction.getTerminalAtom2Index()], coords[iaction.getOutOfPlaneAtomIndex()],
772  grad[iaction.getTerminalAtom1Index()], grad[iaction.getCenterAtomIndex()],
773  grad[iaction.getTerminalAtom2Index()], grad[iaction.getOutOfPlaneAtomIndex()],
774  iaction.getForceConstant());
775 }
776 
777 template <typename ValueType, typename CoordsVec, typename GradVec>
778 ValueType CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
779  const CoordsVec& oop_atom_pos, GradVec& term_atom1_grad, GradVec& ctr_atom_grad,
780  GradVec& term_atom2_grad, GradVec& oop_atom_grad, const ValueType& force_const)
781 {
782  ValueType ac_term1_grad[3];
783  ValueType ac_term2_grad[3];
784  ValueType ac_ctr_grad[3];
785  ValueType ac_oop_grad[3];
786 
787  ValueType chi_ijkl_cos = calcOutOfPlaneAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, oop_atom_pos, ac_term1_grad,
788  ac_ctr_grad, ac_term2_grad, ac_oop_grad);
789  ValueType chi_ijkl = ValueType(M_PI * 0.5) - std::acos(chi_ijkl_cos);
790  ValueType div = std::sqrt(1 - chi_ijkl_cos * chi_ijkl_cos);
791 
792  if (div < ValueType(0.0000001))
793  div = ValueType(0.0000001);
794 
795  ValueType grad_fact = ValueType(0.043844 * 180 * 180) / div * ValueType(1 / (M_PI * M_PI)) * force_const * chi_ijkl;
796 
797  Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
798  Detail::scaleAddVector(ac_ctr_grad, grad_fact, ctr_atom_grad);
799  Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
800  Detail::scaleAddVector(ac_oop_grad, grad_fact, oop_atom_grad);
801 
802  chi_ijkl *= ValueType(180 / M_PI);
803 
804  ValueType e_oop = ValueType(0.5 * 0.043844) * force_const * chi_ijkl * chi_ijkl;
805 
806  return e_oop;
807 }
808 
809 
810 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
811 ValueType CDPL::ForceField::calcMMFF94TorsionGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
812 {
813  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
814  static_cast<ValueType (*)(const MMFF94TorsionInteraction&, const CoordsArray&, GradVector&)>(
815  &calcMMFF94TorsionGradient<ValueType, CoordsArray, GradVector>));
816 }
817 
818 template <typename ValueType, typename CoordsArray, typename GradVector>
819 ValueType CDPL::ForceField::calcMMFF94TorsionGradient(const MMFF94TorsionInteraction& iaction, const CoordsArray& coords, GradVector& grad)
820 {
821  return calcMMFF94TorsionGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtom1Index()],
822  coords[iaction.getCenterAtom2Index()], coords[iaction.getTerminalAtom2Index()],
823  grad[iaction.getTerminalAtom1Index()], grad[iaction.getCenterAtom1Index()], grad[iaction.getCenterAtom2Index()],
824  grad[iaction.getTerminalAtom2Index()], iaction.getTorsionParameter1(), iaction.getTorsionParameter2(),
825  iaction.getTorsionParameter3());
826 }
827 
828 template <typename ValueType, typename CoordsVec, typename GradVec>
829 ValueType
830 CDPL::ForceField::calcMMFF94TorsionGradient(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos, const CoordsVec& ctr_atom2_pos,
831  const CoordsVec& term_atom2_pos, GradVec& term_atom1_grad, GradVec& ctr_atom1_grad, GradVec& ctr_atom2_grad,
832  GradVec& term_atom2_grad, const ValueType& tor_param1, const ValueType& tor_param2, const ValueType& tor_param3)
833 {
834  ValueType ac_term1_grad[3];
835  ValueType ac_ctr1_grad[3];
836  ValueType ac_ctr2_grad[3];
837  ValueType ac_term2_grad[3];
838 
839  ValueType phi_cos = calcDihedralAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom1_pos, ctr_atom2_pos, term_atom2_pos, ac_term1_grad, ac_ctr1_grad,
840  ac_ctr2_grad, ac_term2_grad);
841  ValueType phi = std::acos(phi_cos);
842  ValueType phi_cos_2 = phi_cos * phi_cos;
843  ValueType div = std::sqrt(1 - phi_cos_2);
844 
845  if (div < ValueType(0.0000001))
846  div = ValueType(0.0000001);
847 
848  ValueType grad_fact =
849  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));
850 
851  Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
852  Detail::scaleAddVector(ac_ctr1_grad, grad_fact, ctr_atom1_grad);
853  Detail::scaleAddVector(ac_ctr2_grad, grad_fact, ctr_atom2_grad);
854  Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
855 
856  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)));
857 
858  return e_t;
859 }
860 
861 
862 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
863 ValueType CDPL::ForceField::calcMMFF94ElectrostaticGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
864 {
865  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
866  static_cast<ValueType (*)(const MMFF94ElectrostaticInteraction&, const CoordsArray&, GradVector&)>(
867  &calcMMFF94ElectrostaticGradient<ValueType, CoordsArray, GradVector>));
868 }
869 
870 template <typename ValueType, typename CoordsArray, typename GradVector>
871 ValueType CDPL::ForceField::calcMMFF94ElectrostaticGradient(const MMFF94ElectrostaticInteraction& iaction, const CoordsArray& coords, GradVector& grad)
872 {
873  return calcMMFF94ElectrostaticGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
874  grad[iaction.getAtom2Index()], iaction.getAtom1Charge(), iaction.getAtom2Charge(),
875  iaction.getScalingFactor(), iaction.getDielectricConstant(), iaction.getDistanceExponent());
876 }
877 
878 template <typename ValueType, typename CoordsVec, typename GradVec>
879 ValueType CDPL::ForceField::calcMMFF94ElectrostaticGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
880  const ValueType& atom1_chg, const ValueType& atom2_chg, const ValueType& scale_fact,
881  const ValueType& de_const, const ValueType& dist_expo)
882 {
883  ValueType dist_atom1_grad[3];
884  ValueType dist_atom2_grad[3];
885 
886  ValueType r_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad);
887 
888  ValueType tmp1 = r_ij + ValueType(0.05);
889  ValueType tmp2 = std::pow(tmp1, dist_expo);
890  ValueType tmp3 = scale_fact * atom1_chg * atom2_chg / (de_const * tmp2);
891 
892  ValueType grad_fact = ValueType(-332.0716) * dist_expo * tmp3 / tmp1;
893 
894  Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
895  Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
896 
897  double e_q = ValueType(332.0716) * tmp3;
898 
899  return e_q;
900 }
901 
902 
903 template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector>
904 ValueType CDPL::ForceField::calcMMFF94VanDerWaalsGradient(Iter beg, const Iter& end, const CoordsArray& coords, GradVector& grad)
905 {
906  return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
907  static_cast<ValueType (*)(const MMFF94VanDerWaalsInteraction&, const CoordsArray&, GradVector&)>(
908  &calcMMFF94VanDerWaalsGradient<ValueType, CoordsArray, GradVector>));
909 }
910 
911 template <typename ValueType, typename CoordsArray, typename GradVector>
912 ValueType CDPL::ForceField::calcMMFF94VanDerWaalsGradient(const MMFF94VanDerWaalsInteraction& iaction, const CoordsArray& coords, GradVector& grad)
913 {
914  return calcMMFF94VanDerWaalsGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
915  grad[iaction.getAtom2Index()], iaction.getEIJ(), iaction.getRIJ(), iaction.getRIJPow7());
916 }
917 
918 template <typename ValueType, typename CoordsVec, typename GradVec>
919 ValueType CDPL::ForceField::calcMMFF94VanDerWaalsGradient(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos, GradVec& atom1_grad, GradVec& atom2_grad,
920  const ValueType& e_IJ, const ValueType& r_IJ, const ValueType& r_IJ_7)
921 {
922  ValueType dist_atom1_grad[3];
923  ValueType dist_atom2_grad[3];
924 
925  ValueType r_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad);
926  ValueType r_ij_2 = r_ij * r_ij;
927  ValueType r_ij_6 = r_ij_2 * r_ij_2 * r_ij_2;
928  ValueType r_ij_7 = r_ij_6 * r_ij;
929 
930  ValueType tmp1 = r_ij + ValueType(0.07) * r_IJ;
931  ValueType tmp1_2 = tmp1 * tmp1;
932  ValueType tmp1_4 = tmp1_2 * tmp1_2;
933 
934  ValueType tmp2 = r_ij_7 + ValueType(0.12) * r_IJ_7;
935 
936  ValueType tmp3 = ValueType(1.07) * r_IJ / (r_ij + ValueType(0.07) * r_IJ);
937  ValueType tmp3_2 = tmp3 * tmp3;
938  ValueType tmp3_7 = tmp3_2 * tmp3_2 * tmp3_2 * tmp3;
939 
940  ValueType grad_fact = -r_IJ_7 * e_IJ / (tmp1_4 * tmp1_4 * tmp2 * tmp2) *
941  (ValueType(-22.48094067) * r_ij_7 * r_ij_7 + ValueType(19.78322779) * r_ij_7 * r_IJ_7 +
942  ValueType(0.8812528743) * r_ij_6 * r_IJ_7 * r_IJ + ValueType(1.186993667) * r_IJ_7 * r_IJ_7);
943 
944  Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
945  Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
946 
947  ValueType e_vdw = e_IJ * tmp3_7 * (ValueType(1.12) * r_IJ_7 / tmp2 - 2);
948 
949  return e_vdw;
950 }
951 
952 // \endcond
953 
954 #endif // CDPL_FORCEFIELD_MMFF94GRADIENTFUNCTIONS_HPP
UtilityFunctions.hpp
Utility functions used in the calculation of force field energies and gradients.
MMFF94AngleBendingInteraction.hpp
Definition of class CDPL::ForceField::MMFF94AngleBendingInteraction.
MMFF94BondStretchingInteraction.hpp
Definition of class CDPL::ForceField::MMFF94BondStretchingInteraction.
MMFF94ElectrostaticInteraction.hpp
Definition of class CDPL::ForceField::MMFF94ElectrostaticInteraction.
MMFF94OutOfPlaneBendingInteraction.hpp
Definition of class CDPL::ForceField::MMFF94OutOfPlaneBendingInteraction.
MMFF94StretchBendInteraction.hpp
Definition of class CDPL::ForceField::MMFF94StretchBendInteraction.
MMFF94TorsionInteraction.hpp
Definition of class CDPL::ForceField::MMFF94TorsionInteraction.
MMFF94VanDerWaalsInteraction.hpp
Definition of class CDPL::ForceField::MMFF94VanDerWaalsInteraction.
CDPL::ForceField::MMFF94AngleBendingInteraction
Stores parameters for a single MMFF94 angle-bending interaction defined over an atom triplet.
Definition: MMFF94AngleBendingInteraction.hpp:46
CDPL::ForceField::MMFF94BondStretchingInteraction
Stores parameters for a single MMFF94 bond-stretching interaction between two bonded atoms.
Definition: MMFF94BondStretchingInteraction.hpp:46
CDPL::ForceField::MMFF94ElectrostaticInteraction
Stores parameters for a single MMFF94 electrostatic interaction between two non-bonded atoms.
Definition: MMFF94ElectrostaticInteraction.hpp:46
CDPL::ForceField::MMFF94OutOfPlaneBendingInteraction
Stores parameters for a single MMFF94 out-of-plane bending interaction at a trigonal center.
Definition: MMFF94OutOfPlaneBendingInteraction.hpp:46
CDPL::ForceField::MMFF94StretchBendInteraction
Stores paramters for a single MMFF94 stretch-bend coupling interaction.
Definition: MMFF94StretchBendInteraction.hpp:46
CDPL::ForceField::MMFF94TorsionInteraction
Stores parameters for a single MMFF94 torsion interaction over an atom quadruplet i-j-k-l.
Definition: MMFF94TorsionInteraction.hpp:46
CDPL::ForceField::MMFF94VanDerWaalsInteraction
Stores parameters for a single MMFF94 Van der Waals interaction between two non-bonded atoms.
Definition: MMFF94VanDerWaalsInteraction.hpp:47
CDPL::ForceField::calcMMFF94OutOfPlaneBendingGradient
ValueType calcMMFF94OutOfPlaneBendingGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
Accumulates the out-of-plane bending interaction energies and the corresponding per-atom gradient con...
CDPL::ForceField::calcMMFF94ElectrostaticGradient
ValueType calcMMFF94ElectrostaticGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
Accumulates the electrostatic interaction energies and the corresponding per-atom gradient contributi...
CDPL::ForceField::calcMMFF94VanDerWaalsGradient
ValueType calcMMFF94VanDerWaalsGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
Accumulates the Van der Waals interaction energies and the corresponding per-atom gradient contributi...
CDPL::ForceField::calcMMFF94AngleBendingGradient
ValueType calcMMFF94AngleBendingGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
Accumulates the angle-bending interaction energies and the corresponding per-atom gradient contributi...
CDPL::ForceField::calcMMFF94StretchBendGradient
ValueType calcMMFF94StretchBendGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
Accumulates the stretch-bend coupling interaction energies and the corresponding per-atom gradient co...
CDPL::ForceField::calcMMFF94TorsionGradient
ValueType calcMMFF94TorsionGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
Accumulates the torsion interaction energies and the corresponding per-atom gradient contributions ca...
CDPL::ForceField::calcMMFF94BondStretchingGradient
ValueType calcMMFF94BondStretchingGradient(Iter beg, const Iter &end, const CoordsArray &coords, GradVector &grad)
Accumulates the bond-stretching interaction energies and the corresponding per-atom gradient contribu...
CDPL
The namespace of the Chemical Data Processing Library.
Generated by doxygen 1.9.1