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