29 #ifndef CDPL_FORCEFIELD_UTILITYFUNCTIONS_HPP
30 #define CDPL_FORCEFIELD_UTILITYFUNCTIONS_HPP
44 class MMFF94InteractionData;
63 template <
typename ValueType,
typename CoordsVec>
81 template <
typename ValueType,
typename CoordsVec>
82 ValueType
calcDistance(
const CoordsVec& atom1_pos,
const CoordsVec& atom2_pos);
107 template <
typename ValueType,
typename CoordsVec>
109 ValueType& bond_length1, ValueType& bond_length2);
134 template <
typename ValueType,
typename CoordsVec>
135 ValueType
calcBondLengthsAndAngle(
const CoordsVec& term_atom1_pos,
const CoordsVec& ctr_atom_pos,
const CoordsVec& term_atom2_pos,
136 ValueType& bond_length1, ValueType& bond_length2);
156 template <
typename ValueType,
typename CoordsVec>
157 ValueType
calcBondAngleCos(
const CoordsVec& term_atom1_pos,
const CoordsVec& ctr_atom_pos,
const CoordsVec& term_atom2_pos);
179 template <
typename ValueType,
typename CoordsVec>
180 ValueType
calcBondAngleCos(
const CoordsVec& term_atom1_pos,
const CoordsVec& ctr_atom_pos,
const CoordsVec& term_atom2_pos,
181 const ValueType& r_ij,
const ValueType& r_jk);
201 template <
typename ValueType,
typename CoordsVec>
202 ValueType
calcBondAngle(
const CoordsVec& term_atom1_pos,
const CoordsVec& ctr_atom_pos,
const CoordsVec& term_atom2_pos);
224 template <
typename ValueType,
typename CoordsVec>
225 ValueType
calcBondAngle(
const CoordsVec& term_atom1_pos,
const CoordsVec& ctr_atom_pos,
const CoordsVec& term_atom2_pos,
226 const ValueType& r_ij,
const ValueType& r_jk);
250 template <
typename ValueType,
typename CoordsVec>
252 const CoordsVec& term_atom2_pos,
const CoordsVec& oop_atom_pos);
277 template <
typename ValueType,
typename CoordsVec>
279 const CoordsVec& term_atom2_pos,
const CoordsVec& oop_atom_pos,
280 const ValueType& r_jl);
305 template <
typename ValueType,
typename CoordsVec>
307 const CoordsVec& ctr_atom2_pos,
const CoordsVec& term_atom2_pos);
331 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
333 GradVec& atom1_deriv, GradVec& atom2_deriv);
365 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
367 GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv, GradVec& term_atom2_deriv);
410 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
412 const CoordsVec& ctr_atom2_pos,
const CoordsVec& term_atom2_pos,
413 GradVec& term_atom1_deriv, GradVec& ctr_atom1_deriv,
414 GradVec& ctr_atom2_deriv, GradVec& term_atom2_deriv);
477 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
479 const CoordsVec& term_atom2_pos,
const CoordsVec& oop_atom_pos,
480 GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv,
481 GradVec& term_atom2_deriv, GradVec& oop_atom_deriv);
498 template <
typename VecType1,
typename VecType2,
typename VecType3>
499 void addVectors(
const VecType1& vec1,
const VecType2& vec2, VecType3& res)
501 res[0] = vec2[0] + vec1[0];
502 res[1] = vec2[1] + vec1[1];
503 res[2] = vec2[2] + vec1[2];
506 template <
typename VecType1,
typename VecType2>
507 void subVectors(
const VecType1& vec1,
const VecType1& vec2, VecType2& res)
509 res[0] = vec2[0] - vec1[0];
510 res[1] = vec2[1] - vec1[1];
511 res[2] = vec2[2] - vec1[2];
514 template <
typename ValueType,
typename VecType>
515 ValueType calcDotProduct(
const VecType& vec1,
const VecType& vec2)
517 return (vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]);
520 template <
typename VecType,
typename ResVecType>
521 void calcCrossProduct(
const VecType& vec1,
const VecType& vec2, ResVecType& cross_prod)
523 cross_prod[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
524 cross_prod[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
525 cross_prod[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
528 template <
typename VecType1,
typename VecType2>
529 void copyVector(
const VecType1& vec1, VecType2& vec2)
536 template <
typename VecType>
537 void negateVector(VecType&
vec)
544 template <
typename VecType1,
typename VecType2>
545 void negateCopyVector(
const VecType1& vec1, VecType2& vec2)
552 template <
typename VecType,
typename T>
553 void scaleVector(VecType&
vec,
const T& factor)
560 template <
typename VecType,
typename T>
561 void invScaleVector(VecType&
vec,
const T& factor)
568 template <
typename VecType1,
typename VecType2,
typename T>
569 void scaleAddVector(
const VecType1& vec1,
const T& factor, VecType2& vec2)
571 vec2[0] += vec1[0] * factor;
572 vec2[1] += vec1[1] * factor;
573 vec2[2] += vec1[2] * factor;
576 template <
typename VecType1,
typename VecType2,
typename T>
577 void scaleCopyVector(
const VecType1& vec1,
const T& factor, VecType2& vec2)
579 vec2[0] = vec1[0] * factor;
580 vec2[1] = vec1[1] * factor;
581 vec2[2] = vec1[2] * factor;
584 template <
typename VecType1,
typename VecType2,
typename T>
585 void invScaleCopyVector(
const VecType1& vec1,
const T& factor, VecType2& vec2)
587 vec2[0] = vec1[0] / factor;
588 vec2[1] = vec1[1] / factor;
589 vec2[2] = vec1[2] / factor;
592 template <
typename ValueType>
593 ValueType clampCosine(
const ValueType& v)
595 if (v > ValueType(1))
598 if (v < ValueType(-1))
599 return ValueType(-1);
604 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename FuncType>
605 ValueType accumInteractionEnergies(Iter& beg,
const Iter& end,
const CoordsArray& coords,
const FuncType& func)
607 ValueType e = ValueType();
609 for (; beg != end; ++beg)
610 e += func(*beg, coords);
615 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector,
typename FuncType>
616 ValueType calcInteractionGradient(Iter& beg,
const Iter& end,
const CoordsArray& coords, GradVector& grad,
const FuncType& func)
618 ValueType e = ValueType();
620 for (; beg != end; ++beg)
621 e += func(*beg, coords, grad);
630 template <
typename ValueType,
typename CoordsVec>
633 ValueType pos_diff[3];
635 Detail::subVectors(atom1_pos, atom2_pos, pos_diff);
637 return Detail::calcDotProduct<ValueType>(pos_diff, pos_diff);
640 template <
typename ValueType,
typename CoordsVec>
643 return std::sqrt(calcSquaredDistance<ValueType>(atom1_pos, atom2_pos));
646 template <
typename ValueType,
typename CoordsVec>
648 ValueType& bond_length1, ValueType& bond_length2)
650 ValueType bond_vec1[3];
651 ValueType bond_vec2[3];
653 Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
654 Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
656 bond_length1 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec1));
657 bond_length2 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec2, bond_vec2));
659 return Detail::clampCosine(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2) / (bond_length1 * bond_length2));
662 template <
typename ValueType,
typename CoordsVec>
664 ValueType& bond_length1, ValueType& bond_length2)
669 template <
typename ValueType,
typename CoordsVec>
677 template <
typename ValueType,
typename CoordsVec>
679 const ValueType& r_ij,
const ValueType& r_jk)
681 ValueType bond_vec1[3];
682 ValueType bond_vec2[3];
684 Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
685 Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
687 return Detail::clampCosine(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2) / (r_ij * r_jk));
690 template <
typename ValueType,
typename CoordsVec>
693 return std::acos(calcBondAngleCos<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos));
696 template <
typename ValueType,
typename CoordsVec>
698 const ValueType& r_ij,
const ValueType& r_jk)
700 return std::acos(calcBondAngleCos<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, r_ij, r_jk));
703 template <
typename ValueType,
typename CoordsVec>
705 const CoordsVec& term_atom2_pos,
const CoordsVec& oop_atom_pos)
707 ValueType term_bond1_vec[3];
708 ValueType term_bond2_vec[3];
709 ValueType oop_bond_vec[3];
710 ValueType plane_normal[3];
712 Detail::subVectors(ctr_atom_pos, term_atom1_pos, term_bond1_vec);
713 Detail::subVectors(ctr_atom_pos, term_atom2_pos, term_bond2_vec);
714 Detail::subVectors(ctr_atom_pos, oop_atom_pos, oop_bond_vec);
715 Detail::calcCrossProduct(term_bond1_vec, term_bond2_vec, plane_normal);
717 ValueType pn_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal, plane_normal));
718 ValueType oop_bnd_len = std::sqrt(Detail::calcDotProduct<ValueType>(oop_bond_vec, oop_bond_vec));
719 ValueType ang_cos = Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal, oop_bond_vec) / (pn_len * oop_bnd_len));
721 return (ValueType(M_PI * 0.5) - std::acos(ang_cos));
724 template <
typename ValueType,
typename CoordsVec>
726 const CoordsVec& term_atom2_pos,
const CoordsVec& oop_atom_pos,
727 const ValueType& r_jl)
729 ValueType term_bond1_vec[3];
730 ValueType term_bond2_vec[3];
731 ValueType oop_bond_vec[3];
732 ValueType plane_normal[3];
734 Detail::subVectors(ctr_atom_pos, term_atom1_pos, term_bond1_vec);
735 Detail::subVectors(ctr_atom_pos, term_atom2_pos, term_bond2_vec);
736 Detail::subVectors(ctr_atom_pos, oop_atom_pos, oop_bond_vec);
737 Detail::calcCrossProduct(term_bond1_vec, term_bond2_vec, plane_normal);
739 ValueType pn_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal, plane_normal));
740 ValueType ang_cos = Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal, oop_bond_vec) / (pn_len * r_jl));
742 return (ValueType(M_PI * 0.5) - std::acos(ang_cos));
745 template <
typename ValueType,
typename CoordsVec>
747 const CoordsVec& ctr_atom2_pos,
const CoordsVec& term_atom2_pos)
749 ValueType term_bond1_vec[3];
750 ValueType ctr_bond_vec[3];
751 ValueType term_bond2_vec[3];
752 ValueType plane_normal1[3];
753 ValueType plane_normal2[3];
755 Detail::subVectors(ctr_atom1_pos, term_atom1_pos, term_bond1_vec);
756 Detail::subVectors(ctr_atom1_pos, ctr_atom2_pos, ctr_bond_vec);
757 Detail::subVectors(term_atom2_pos, ctr_atom2_pos, term_bond2_vec);
759 Detail::calcCrossProduct(term_bond1_vec, ctr_bond_vec, plane_normal1);
760 Detail::calcCrossProduct(ctr_bond_vec, term_bond2_vec, plane_normal2);
762 ValueType pn1_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal1));
763 ValueType pn2_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal2, plane_normal2));
765 return Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal2) / (pn1_len * pn2_len));
768 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
770 GradVec& atom1_deriv, GradVec& atom2_deriv)
772 Detail::subVectors(atom1_pos, atom2_pos, atom1_deriv);
774 ValueType dist = std::sqrt(Detail::calcDotProduct<ValueType>(atom1_deriv, atom1_deriv));
776 Detail::invScaleVector(atom1_deriv, -dist);
777 Detail::negateCopyVector(atom1_deriv, atom2_deriv);
782 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
784 GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv, GradVec& term_atom2_deriv)
786 ValueType bond_vec1[3];
787 ValueType bond_vec2[3];
789 Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
790 Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
792 ValueType bond_length1 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec1));
793 ValueType bond_length2 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec2, bond_vec2));
795 ValueType dot_prod = Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2);
796 ValueType bl_prod = bond_length1 * bond_length2;
798 Detail::invScaleCopyVector(bond_vec2, bl_prod, term_atom1_deriv);
799 Detail::scaleCopyVector(bond_vec1, dot_prod / (bond_length1 * bond_length1 * bl_prod), ctr_atom_deriv);
800 Detail::subVectors(ctr_atom_deriv, term_atom1_deriv, term_atom1_deriv);
802 Detail::invScaleCopyVector(bond_vec1, bl_prod, term_atom2_deriv);
803 Detail::scaleCopyVector(bond_vec2, dot_prod / (bond_length2 * bond_length2 * bl_prod), ctr_atom_deriv);
804 Detail::subVectors(ctr_atom_deriv, term_atom2_deriv, term_atom2_deriv);
806 Detail::negateCopyVector(term_atom1_deriv, ctr_atom_deriv);
807 Detail::subVectors(term_atom2_deriv, ctr_atom_deriv, ctr_atom_deriv);
809 return Detail::clampCosine(dot_prod / bl_prod);
812 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
814 const CoordsVec& ctr_atom2_pos,
const CoordsVec& term_atom2_pos,
815 GradVec& term_atom1_deriv, GradVec& ctr_atom1_deriv,
816 GradVec& ctr_atom2_deriv, GradVec& term_atom2_deriv)
818 ValueType term_bond1_vec[3];
819 ValueType ctr_bond_vec[3];
820 ValueType term_bond2_vec[3];
821 ValueType term1_ctr2_vec[3];
822 ValueType plane_normal1[3];
823 ValueType plane_normal2[3];
825 Detail::subVectors(ctr_atom1_pos, term_atom1_pos, term_bond1_vec);
826 Detail::subVectors(ctr_atom1_pos, ctr_atom2_pos, ctr_bond_vec);
827 Detail::subVectors(term_atom2_pos, ctr_atom2_pos, term_bond2_vec);
828 Detail::subVectors(ctr_atom2_pos, term_atom1_pos, term1_ctr2_vec);
830 Detail::calcCrossProduct(term_bond1_vec, ctr_bond_vec, plane_normal1);
831 Detail::calcCrossProduct(ctr_bond_vec, term_bond2_vec, plane_normal2);
833 ValueType pn1_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal1));
834 ValueType pn2_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal2, plane_normal2));
836 Detail::invScaleVector(plane_normal1, pn1_len);
837 Detail::invScaleVector(plane_normal2, pn2_len);
839 ValueType ang_cos = Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal2));
844 Detail::scaleCopyVector(plane_normal1, -ang_cos, a);
845 Detail::addVectors(a, plane_normal2, a);
846 Detail::invScaleVector(a, pn1_len);
848 Detail::scaleCopyVector(plane_normal2, -ang_cos, b);
849 Detail::addVectors(b, plane_normal1, b);
850 Detail::invScaleVector(b, pn2_len);
852 Detail::calcCrossProduct(ctr_bond_vec, a, term_atom1_deriv);
853 Detail::calcCrossProduct(ctr_bond_vec, b, term_atom2_deriv);
855 Detail::calcCrossProduct(term1_ctr2_vec, a, ctr_atom1_deriv);
856 Detail::calcCrossProduct(term_bond2_vec, b, ctr_atom2_deriv);
858 Detail::subVectors(ctr_atom2_deriv, ctr_atom1_deriv, ctr_atom1_deriv);
860 Detail::addVectors(term_atom1_deriv, ctr_atom1_deriv, ctr_atom2_deriv);
861 Detail::addVectors(ctr_atom2_deriv, term_atom2_deriv, ctr_atom2_deriv);
862 Detail::negateVector(ctr_atom2_deriv);
867 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
869 const CoordsVec& term_atom2_pos,
const CoordsVec& oop_atom_pos,
870 GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv,
871 GradVec& term_atom2_deriv, GradVec& oop_atom_deriv)
873 ValueType term_bond1_vec[3];
874 ValueType term_bond2_vec[3];
875 ValueType oop_bond_vec[3];
876 ValueType term1_oop_vec[3];
877 ValueType term2_oop_vec[3];
880 Detail::subVectors(ctr_atom_pos, term_atom1_pos, term_bond1_vec);
881 Detail::subVectors(ctr_atom_pos, term_atom2_pos, term_bond2_vec);
882 Detail::subVectors(ctr_atom_pos, oop_atom_pos, oop_bond_vec);
883 Detail::subVectors(term_atom2_pos, oop_atom_pos, term2_oop_vec);
884 Detail::subVectors(term_atom1_pos, oop_atom_pos, term1_oop_vec);
886 Detail::calcCrossProduct(term_bond1_vec, term_bond2_vec, ijk_pn);
888 ValueType ijk_pn_len_2 = Detail::calcDotProduct<ValueType>(ijk_pn, ijk_pn);
889 ValueType ijk_pn_len = std::sqrt(ijk_pn_len_2);
890 ValueType oop_bnd_len = std::sqrt(Detail::calcDotProduct<ValueType>(oop_bond_vec, oop_bond_vec));
891 ValueType oop_ijk_pn_len_prod = ijk_pn_len * oop_bnd_len;
892 ValueType oop_ijk_pn_dot_prod = Detail::calcDotProduct<ValueType>(ijk_pn, oop_bond_vec);
893 ValueType ang_cos = Detail::clampCosine(oop_ijk_pn_dot_prod / oop_ijk_pn_len_prod);
898 Detail::calcCrossProduct(term_bond2_vec, oop_bond_vec, kjl_pn);
899 Detail::calcCrossProduct(oop_bond_vec, term_bond1_vec, lji_pn);
901 ctr_atom_deriv[0] = ijk_pn[1] * -term_bond2_vec[2] + ijk_pn[2] * term_bond2_vec[1];
902 ctr_atom_deriv[1] = ijk_pn[0] * term_bond2_vec[2] + ijk_pn[2] * -term_bond2_vec[0];
903 ctr_atom_deriv[2] = ijk_pn[0] * -term_bond2_vec[1] + ijk_pn[1] * term_bond2_vec[0];
905 Detail::scaleVector(ctr_atom_deriv, ang_cos / ijk_pn_len_2);
907 Detail::invScaleCopyVector(kjl_pn, oop_ijk_pn_len_prod, term_atom1_deriv);
908 Detail::subVectors(ctr_atom_deriv, term_atom1_deriv, term_atom1_deriv);
910 ctr_atom_deriv[0] = ijk_pn[1] * term_bond1_vec[2] + ijk_pn[2] * -term_bond1_vec[1];
911 ctr_atom_deriv[1] = ijk_pn[0] * -term_bond1_vec[2] + ijk_pn[2] * term_bond1_vec[0];
912 ctr_atom_deriv[2] = ijk_pn[0] * term_bond1_vec[1] + ijk_pn[1] * -term_bond1_vec[0];
914 Detail::scaleVector(ctr_atom_deriv, ang_cos / ijk_pn_len_2);
916 Detail::invScaleCopyVector(lji_pn, oop_ijk_pn_len_prod, term_atom2_deriv);
917 Detail::subVectors(ctr_atom_deriv, term_atom2_deriv, term_atom2_deriv);
919 Detail::calcCrossProduct(term2_oop_vec, term1_oop_vec, ctr_atom_deriv);
921 Detail::scaleCopyVector(oop_bond_vec, oop_ijk_pn_dot_prod / (oop_bnd_len * oop_bnd_len), oop_atom_deriv);
922 Detail::addVectors(oop_atom_deriv, lji_pn, oop_atom_deriv);
923 Detail::addVectors(oop_atom_deriv, kjl_pn, oop_atom_deriv);
924 Detail::addVectors(oop_atom_deriv, ctr_atom_deriv, oop_atom_deriv);
925 Detail::invScaleVector(oop_atom_deriv, -oop_ijk_pn_len_prod);
927 Detail::copyVector(term_atom1_deriv, ctr_atom_deriv);
928 Detail::addVectors(term_atom2_deriv, ctr_atom_deriv, ctr_atom_deriv);
929 Detail::addVectors(oop_atom_deriv, ctr_atom_deriv, ctr_atom_deriv);
930 Detail::negateVector(ctr_atom_deriv);
937 #endif // CDPL_FORCEFIELD_UTILITYFUNCTIONS_HPP