Chemical Data Processing Library C++ API - Version 1.4.0
ForceField/UtilityFunctions.hpp
Go to the documentation of this file.
1 /*
2  * UtilityFunctions.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_UTILITYFUNCTIONS_HPP
30 #define CDPL_FORCEFIELD_UTILITYFUNCTIONS_HPP
31 
32 #include <cmath>
33 
34 #include "CDPL/ForceField/APIPrefix.hpp"
35 #include "CDPL/Util/BitSet.hpp"
36 
37 
38 namespace CDPL
39 {
40 
41  namespace ForceField
42  {
43 
44  class MMFF94InteractionData;
45 
52  CDPL_FORCEFIELD_API void filterInteractions(const MMFF94InteractionData& ia_data, MMFF94InteractionData& filtered_ia_data,
53  const Util::BitSet& inc_atom_mask);
54 
70  template <typename ValueType, typename CoordsVec>
71  ValueType calcSquaredDistance(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos);
72 
88  template <typename ValueType, typename CoordsVec>
89  ValueType calcDistance(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos);
90 
114  template <typename ValueType, typename CoordsVec>
115  ValueType calcBondLengthsAndAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
116  ValueType& bond_length1, ValueType& bond_length2);
117 
141  template <typename ValueType, typename CoordsVec>
142  ValueType calcBondLengthsAndAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
143  ValueType& bond_length1, ValueType& bond_length2);
144 
163  template <typename ValueType, typename CoordsVec>
164  ValueType calcBondAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos);
165 
186  template <typename ValueType, typename CoordsVec>
187  ValueType calcBondAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
188  const ValueType& r_ij, const ValueType& r_jk);
189 
208  template <typename ValueType, typename CoordsVec>
209  ValueType calcBondAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos);
210 
231  template <typename ValueType, typename CoordsVec>
232  ValueType calcBondAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
233  const ValueType& r_ij, const ValueType& r_jk);
234 
257  template <typename ValueType, typename CoordsVec>
258  ValueType calcOutOfPlaneAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
259  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos);
260 
284  template <typename ValueType, typename CoordsVec>
285  ValueType calcOutOfPlaneAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
286  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos,
287  const ValueType& r_jl);
288 
312  template <typename ValueType, typename CoordsVec>
313  ValueType calcDihedralAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos,
314  const CoordsVec& ctr_atom2_pos, const CoordsVec& term_atom2_pos);
315 
338  template <typename ValueType, typename CoordsVec, typename GradVec>
339  ValueType calcDistanceDerivatives(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos,
340  GradVec& atom1_deriv, GradVec& atom2_deriv);
341 
372  template <typename ValueType, typename CoordsVec, typename GradVec>
373  ValueType calcBondAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
374  GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv, GradVec& term_atom2_deriv);
375 
417  template <typename ValueType, typename CoordsVec, typename GradVec>
418  ValueType calcDihedralAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos,
419  const CoordsVec& ctr_atom2_pos, const CoordsVec& term_atom2_pos,
420  GradVec& term_atom1_deriv, GradVec& ctr_atom1_deriv,
421  GradVec& ctr_atom2_deriv, GradVec& term_atom2_deriv);
422 
482  template <typename ValueType, typename CoordsVec, typename GradVec>
483  ValueType calcOutOfPlaneAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
484  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos,
485  GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv,
486  GradVec& term_atom2_deriv, GradVec& oop_atom_deriv);
487  } // namespace ForceField
488 } // namespace CDPL
489 
490 
491 // Implementation
492 // \cond DOC_IMPL_DETAILS
493 
494 namespace CDPL
495 {
496 
497  namespace ForceField
498  {
499 
500  namespace Detail
501  {
502 
503  template <typename VecType1, typename VecType2, typename VecType3>
504  void addVectors(const VecType1& vec1, const VecType2& vec2, VecType3& res)
505  {
506  res[0] = vec2[0] + vec1[0];
507  res[1] = vec2[1] + vec1[1];
508  res[2] = vec2[2] + vec1[2];
509  }
510 
511  template <typename VecType1, typename VecType2>
512  void subVectors(const VecType1& vec1, const VecType1& vec2, VecType2& res)
513  {
514  res[0] = vec2[0] - vec1[0];
515  res[1] = vec2[1] - vec1[1];
516  res[2] = vec2[2] - vec1[2];
517  }
518 
519  template <typename ValueType, typename VecType>
520  ValueType calcDotProduct(const VecType& vec1, const VecType& vec2)
521  {
522  return (vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]);
523  }
524 
525  template <typename VecType, typename ResVecType>
526  void calcCrossProduct(const VecType& vec1, const VecType& vec2, ResVecType& cross_prod)
527  {
528  cross_prod[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
529  cross_prod[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
530  cross_prod[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
531  }
532 
533  template <typename VecType1, typename VecType2>
534  void copyVector(const VecType1& vec1, VecType2& vec2)
535  {
536  vec2[0] = vec1[0];
537  vec2[1] = vec1[1];
538  vec2[2] = vec1[2];
539  }
540 
541  template <typename VecType>
542  void negateVector(VecType& vec)
543  {
544  vec[0] = -vec[0];
545  vec[1] = -vec[1];
546  vec[2] = -vec[2];
547  }
548 
549  template <typename VecType1, typename VecType2>
550  void negateCopyVector(const VecType1& vec1, VecType2& vec2)
551  {
552  vec2[0] = -vec1[0];
553  vec2[1] = -vec1[1];
554  vec2[2] = -vec1[2];
555  }
556 
557  template <typename VecType, typename T>
558  void scaleVector(VecType& vec, const T& factor)
559  {
560  vec[0] *= factor;
561  vec[1] *= factor;
562  vec[2] *= factor;
563  }
564 
565  template <typename VecType, typename T>
566  void invScaleVector(VecType& vec, const T& factor)
567  {
568  vec[0] /= factor;
569  vec[1] /= factor;
570  vec[2] /= factor;
571  }
572 
573  template <typename VecType1, typename VecType2, typename T>
574  void scaleAddVector(const VecType1& vec1, const T& factor, VecType2& vec2)
575  {
576  vec2[0] += vec1[0] * factor;
577  vec2[1] += vec1[1] * factor;
578  vec2[2] += vec1[2] * factor;
579  }
580 
581  template <typename VecType1, typename VecType2, typename T>
582  void scaleCopyVector(const VecType1& vec1, const T& factor, VecType2& vec2)
583  {
584  vec2[0] = vec1[0] * factor;
585  vec2[1] = vec1[1] * factor;
586  vec2[2] = vec1[2] * factor;
587  }
588 
589  template <typename VecType1, typename VecType2, typename T>
590  void invScaleCopyVector(const VecType1& vec1, const T& factor, VecType2& vec2)
591  {
592  vec2[0] = vec1[0] / factor;
593  vec2[1] = vec1[1] / factor;
594  vec2[2] = vec1[2] / factor;
595  }
596 
597  template <typename ValueType>
598  ValueType clampCosine(const ValueType& v)
599  {
600  if (v > ValueType(1))
601  return ValueType(1);
602 
603  if (v < ValueType(-1))
604  return ValueType(-1);
605 
606  return v;
607  }
608 
609  template <typename ValueType, typename Iter, typename CoordsArray, typename FuncType>
610  ValueType accumInteractionEnergies(Iter& beg, const Iter& end, const CoordsArray& coords, const FuncType& func)
611  {
612  ValueType e = ValueType();
613 
614  for (; beg != end; ++beg)
615  e += func(*beg, coords);
616 
617  return e;
618  }
619 
620  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector, typename FuncType>
621  ValueType calcInteractionGradient(Iter& beg, const Iter& end, const CoordsArray& coords, GradVector& grad, const FuncType& func)
622  {
623  ValueType e = ValueType();
624 
625  for (; beg != end; ++beg)
626  e += func(*beg, coords, grad);
627 
628  return e;
629  }
630  } // namespace Detail
631  } // namespace ForceField
632 } // namespace CDPL
633 
634 
635 template <typename ValueType, typename CoordsVec>
636 ValueType CDPL::ForceField::calcSquaredDistance(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos)
637 {
638  ValueType pos_diff[3];
639 
640  Detail::subVectors(atom1_pos, atom2_pos, pos_diff);
641 
642  return Detail::calcDotProduct<ValueType>(pos_diff, pos_diff);
643 }
644 
645 template <typename ValueType, typename CoordsVec>
646 ValueType CDPL::ForceField::calcDistance(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos)
647 {
648  return std::sqrt(calcSquaredDistance<ValueType>(atom1_pos, atom2_pos));
649 }
650 
651 template <typename ValueType, typename CoordsVec>
652 ValueType CDPL::ForceField::calcBondLengthsAndAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
653  ValueType& bond_length1, ValueType& bond_length2)
654 {
655  ValueType bond_vec1[3];
656  ValueType bond_vec2[3];
657 
658  Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
659  Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
660 
661  bond_length1 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec1));
662  bond_length2 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec2, bond_vec2));
663 
664  return Detail::clampCosine(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2) / (bond_length1 * bond_length2));
665 }
666 
667 template <typename ValueType, typename CoordsVec>
668 ValueType CDPL::ForceField::calcBondLengthsAndAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
669  ValueType& bond_length1, ValueType& bond_length2)
670 {
671  return std::acos(calcBondLengthsAndAngleCos(term_atom1_pos, ctr_atom_pos, term_atom2_pos, bond_length1, bond_length2));
672 }
673 
674 template <typename ValueType, typename CoordsVec>
675 ValueType CDPL::ForceField::calcBondAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos)
676 {
677  ValueType bl1, bl2;
678 
679  return calcBondLengthsAndAngleCos(term_atom1_pos, ctr_atom_pos, term_atom2_pos, bl1, bl2);
680 }
681 
682 template <typename ValueType, typename CoordsVec>
683 ValueType CDPL::ForceField::calcBondAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
684  const ValueType& r_ij, const ValueType& r_jk)
685 {
686  ValueType bond_vec1[3];
687  ValueType bond_vec2[3];
688 
689  Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
690  Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
691 
692  return Detail::clampCosine(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2) / (r_ij * r_jk));
693 }
694 
695 template <typename ValueType, typename CoordsVec>
696 ValueType CDPL::ForceField::calcBondAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos)
697 {
698  return std::acos(calcBondAngleCos<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos));
699 }
700 
701 template <typename ValueType, typename CoordsVec>
702 ValueType CDPL::ForceField::calcBondAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
703  const ValueType& r_ij, const ValueType& r_jk)
704 {
705  return std::acos(calcBondAngleCos<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, r_ij, r_jk));
706 }
707 
708 template <typename ValueType, typename CoordsVec>
709 ValueType CDPL::ForceField::calcOutOfPlaneAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
710  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos)
711 {
712  ValueType term_bond1_vec[3];
713  ValueType term_bond2_vec[3];
714  ValueType oop_bond_vec[3];
715  ValueType plane_normal[3];
716 
717  Detail::subVectors(ctr_atom_pos, term_atom1_pos, term_bond1_vec);
718  Detail::subVectors(ctr_atom_pos, term_atom2_pos, term_bond2_vec);
719  Detail::subVectors(ctr_atom_pos, oop_atom_pos, oop_bond_vec);
720  Detail::calcCrossProduct(term_bond1_vec, term_bond2_vec, plane_normal);
721 
722  ValueType pn_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal, plane_normal));
723  ValueType oop_bnd_len = std::sqrt(Detail::calcDotProduct<ValueType>(oop_bond_vec, oop_bond_vec));
724  ValueType ang_cos = Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal, oop_bond_vec) / (pn_len * oop_bnd_len));
725 
726  return (ValueType(M_PI * 0.5) - std::acos(ang_cos));
727 }
728 
729 template <typename ValueType, typename CoordsVec>
730 ValueType CDPL::ForceField::calcOutOfPlaneAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
731  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos,
732  const ValueType& r_jl)
733 {
734  ValueType term_bond1_vec[3];
735  ValueType term_bond2_vec[3];
736  ValueType oop_bond_vec[3];
737  ValueType plane_normal[3];
738 
739  Detail::subVectors(ctr_atom_pos, term_atom1_pos, term_bond1_vec);
740  Detail::subVectors(ctr_atom_pos, term_atom2_pos, term_bond2_vec);
741  Detail::subVectors(ctr_atom_pos, oop_atom_pos, oop_bond_vec);
742  Detail::calcCrossProduct(term_bond1_vec, term_bond2_vec, plane_normal);
743 
744  ValueType pn_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal, plane_normal));
745  ValueType ang_cos = Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal, oop_bond_vec) / (pn_len * r_jl));
746 
747  return (ValueType(M_PI * 0.5) - std::acos(ang_cos));
748 }
749 
750 template <typename ValueType, typename CoordsVec>
751 ValueType CDPL::ForceField::calcDihedralAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos,
752  const CoordsVec& ctr_atom2_pos, const CoordsVec& term_atom2_pos)
753 {
754  ValueType term_bond1_vec[3];
755  ValueType ctr_bond_vec[3];
756  ValueType term_bond2_vec[3];
757  ValueType plane_normal1[3];
758  ValueType plane_normal2[3];
759 
760  Detail::subVectors(ctr_atom1_pos, term_atom1_pos, term_bond1_vec);
761  Detail::subVectors(ctr_atom1_pos, ctr_atom2_pos, ctr_bond_vec);
762  Detail::subVectors(term_atom2_pos, ctr_atom2_pos, term_bond2_vec);
763 
764  Detail::calcCrossProduct(term_bond1_vec, ctr_bond_vec, plane_normal1);
765  Detail::calcCrossProduct(ctr_bond_vec, term_bond2_vec, plane_normal2);
766 
767  ValueType pn1_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal1));
768  ValueType pn2_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal2, plane_normal2));
769 
770  return Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal2) / (pn1_len * pn2_len));
771 }
772 
773 template <typename ValueType, typename CoordsVec, typename GradVec>
774 ValueType CDPL::ForceField::calcDistanceDerivatives(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos,
775  GradVec& atom1_deriv, GradVec& atom2_deriv)
776 {
777  Detail::subVectors(atom1_pos, atom2_pos, atom1_deriv);
778 
779  ValueType dist = std::sqrt(Detail::calcDotProduct<ValueType>(atom1_deriv, atom1_deriv));
780 
781  Detail::invScaleVector(atom1_deriv, -dist);
782  Detail::negateCopyVector(atom1_deriv, atom2_deriv);
783 
784  return dist;
785 }
786 
787 template <typename ValueType, typename CoordsVec, typename GradVec>
788 ValueType CDPL::ForceField::calcBondAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
789  GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv, GradVec& term_atom2_deriv)
790 {
791  ValueType bond_vec1[3];
792  ValueType bond_vec2[3];
793 
794  Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
795  Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
796 
797  ValueType bond_length1 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec1));
798  ValueType bond_length2 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec2, bond_vec2));
799 
800  ValueType dot_prod = Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2);
801  ValueType bl_prod = bond_length1 * bond_length2;
802 
803  Detail::invScaleCopyVector(bond_vec2, bl_prod, term_atom1_deriv);
804  Detail::scaleCopyVector(bond_vec1, dot_prod / (bond_length1 * bond_length1 * bl_prod), ctr_atom_deriv);
805  Detail::subVectors(ctr_atom_deriv, term_atom1_deriv, term_atom1_deriv);
806 
807  Detail::invScaleCopyVector(bond_vec1, bl_prod, term_atom2_deriv);
808  Detail::scaleCopyVector(bond_vec2, dot_prod / (bond_length2 * bond_length2 * bl_prod), ctr_atom_deriv);
809  Detail::subVectors(ctr_atom_deriv, term_atom2_deriv, term_atom2_deriv);
810 
811  Detail::negateCopyVector(term_atom1_deriv, ctr_atom_deriv);
812  Detail::subVectors(term_atom2_deriv, ctr_atom_deriv, ctr_atom_deriv);
813 
814  return Detail::clampCosine(dot_prod / bl_prod);
815 }
816 
817 template <typename ValueType, typename CoordsVec, typename GradVec>
818 ValueType CDPL::ForceField::calcDihedralAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos,
819  const CoordsVec& ctr_atom2_pos, const CoordsVec& term_atom2_pos,
820  GradVec& term_atom1_deriv, GradVec& ctr_atom1_deriv,
821  GradVec& ctr_atom2_deriv, GradVec& term_atom2_deriv)
822 {
823  ValueType term_bond1_vec[3];
824  ValueType ctr_bond_vec[3];
825  ValueType term_bond2_vec[3];
826  ValueType term1_ctr2_vec[3];
827  ValueType plane_normal1[3];
828  ValueType plane_normal2[3];
829 
830  Detail::subVectors(ctr_atom1_pos, term_atom1_pos, term_bond1_vec);
831  Detail::subVectors(ctr_atom1_pos, ctr_atom2_pos, ctr_bond_vec);
832  Detail::subVectors(term_atom2_pos, ctr_atom2_pos, term_bond2_vec);
833  Detail::subVectors(ctr_atom2_pos, term_atom1_pos, term1_ctr2_vec);
834 
835  Detail::calcCrossProduct(term_bond1_vec, ctr_bond_vec, plane_normal1);
836  Detail::calcCrossProduct(ctr_bond_vec, term_bond2_vec, plane_normal2);
837 
838  ValueType pn1_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal1));
839  ValueType pn2_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal2, plane_normal2));
840 
841  Detail::invScaleVector(plane_normal1, pn1_len);
842  Detail::invScaleVector(plane_normal2, pn2_len);
843 
844  ValueType ang_cos = Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal2));
845 
846  ValueType a[3];
847  ValueType b[3];
848 
849  Detail::scaleCopyVector(plane_normal1, -ang_cos, a);
850  Detail::addVectors(a, plane_normal2, a);
851  Detail::invScaleVector(a, pn1_len);
852 
853  Detail::scaleCopyVector(plane_normal2, -ang_cos, b);
854  Detail::addVectors(b, plane_normal1, b);
855  Detail::invScaleVector(b, pn2_len);
856 
857  Detail::calcCrossProduct(ctr_bond_vec, a, term_atom1_deriv);
858  Detail::calcCrossProduct(ctr_bond_vec, b, term_atom2_deriv);
859 
860  Detail::calcCrossProduct(term1_ctr2_vec, a, ctr_atom1_deriv);
861  Detail::calcCrossProduct(term_bond2_vec, b, ctr_atom2_deriv);
862 
863  Detail::subVectors(ctr_atom2_deriv, ctr_atom1_deriv, ctr_atom1_deriv);
864 
865  Detail::addVectors(term_atom1_deriv, ctr_atom1_deriv, ctr_atom2_deriv);
866  Detail::addVectors(ctr_atom2_deriv, term_atom2_deriv, ctr_atom2_deriv);
867  Detail::negateVector(ctr_atom2_deriv);
868 
869  return ang_cos;
870 }
871 
872 template <typename ValueType, typename CoordsVec, typename GradVec>
873 ValueType CDPL::ForceField::calcOutOfPlaneAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
874  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos,
875  GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv,
876  GradVec& term_atom2_deriv, GradVec& oop_atom_deriv)
877 {
878  ValueType term_bond1_vec[3];
879  ValueType term_bond2_vec[3];
880  ValueType oop_bond_vec[3];
881  ValueType term1_oop_vec[3];
882  ValueType term2_oop_vec[3];
883  ValueType ijk_pn[3];
884 
885  Detail::subVectors(ctr_atom_pos, term_atom1_pos, term_bond1_vec);
886  Detail::subVectors(ctr_atom_pos, term_atom2_pos, term_bond2_vec);
887  Detail::subVectors(ctr_atom_pos, oop_atom_pos, oop_bond_vec);
888  Detail::subVectors(term_atom2_pos, oop_atom_pos, term2_oop_vec);
889  Detail::subVectors(term_atom1_pos, oop_atom_pos, term1_oop_vec);
890 
891  Detail::calcCrossProduct(term_bond1_vec, term_bond2_vec, ijk_pn);
892 
893  ValueType ijk_pn_len_2 = Detail::calcDotProduct<ValueType>(ijk_pn, ijk_pn);
894  ValueType ijk_pn_len = std::sqrt(ijk_pn_len_2);
895  ValueType oop_bnd_len = std::sqrt(Detail::calcDotProduct<ValueType>(oop_bond_vec, oop_bond_vec));
896  ValueType oop_ijk_pn_len_prod = ijk_pn_len * oop_bnd_len;
897  ValueType oop_ijk_pn_dot_prod = Detail::calcDotProduct<ValueType>(ijk_pn, oop_bond_vec);
898  ValueType ang_cos = Detail::clampCosine(oop_ijk_pn_dot_prod / oop_ijk_pn_len_prod);
899 
900  ValueType kjl_pn[3];
901  ValueType lji_pn[3];
902 
903  Detail::calcCrossProduct(term_bond2_vec, oop_bond_vec, kjl_pn);
904  Detail::calcCrossProduct(oop_bond_vec, term_bond1_vec, lji_pn);
905 
906  ctr_atom_deriv[0] = ijk_pn[1] * -term_bond2_vec[2] + ijk_pn[2] * term_bond2_vec[1];
907  ctr_atom_deriv[1] = ijk_pn[0] * term_bond2_vec[2] + ijk_pn[2] * -term_bond2_vec[0];
908  ctr_atom_deriv[2] = ijk_pn[0] * -term_bond2_vec[1] + ijk_pn[1] * term_bond2_vec[0];
909 
910  Detail::scaleVector(ctr_atom_deriv, ang_cos / ijk_pn_len_2);
911 
912  Detail::invScaleCopyVector(kjl_pn, oop_ijk_pn_len_prod, term_atom1_deriv);
913  Detail::subVectors(ctr_atom_deriv, term_atom1_deriv, term_atom1_deriv);
914 
915  ctr_atom_deriv[0] = ijk_pn[1] * term_bond1_vec[2] + ijk_pn[2] * -term_bond1_vec[1];
916  ctr_atom_deriv[1] = ijk_pn[0] * -term_bond1_vec[2] + ijk_pn[2] * term_bond1_vec[0];
917  ctr_atom_deriv[2] = ijk_pn[0] * term_bond1_vec[1] + ijk_pn[1] * -term_bond1_vec[0];
918 
919  Detail::scaleVector(ctr_atom_deriv, ang_cos / ijk_pn_len_2);
920 
921  Detail::invScaleCopyVector(lji_pn, oop_ijk_pn_len_prod, term_atom2_deriv);
922  Detail::subVectors(ctr_atom_deriv, term_atom2_deriv, term_atom2_deriv);
923 
924  Detail::calcCrossProduct(term2_oop_vec, term1_oop_vec, ctr_atom_deriv);
925 
926  Detail::scaleCopyVector(oop_bond_vec, oop_ijk_pn_dot_prod / (oop_bnd_len * oop_bnd_len), oop_atom_deriv);
927  Detail::addVectors(oop_atom_deriv, lji_pn, oop_atom_deriv);
928  Detail::addVectors(oop_atom_deriv, kjl_pn, oop_atom_deriv);
929  Detail::addVectors(oop_atom_deriv, ctr_atom_deriv, oop_atom_deriv);
930  Detail::invScaleVector(oop_atom_deriv, -oop_ijk_pn_len_prod);
931 
932  Detail::copyVector(term_atom1_deriv, ctr_atom_deriv);
933  Detail::addVectors(term_atom2_deriv, ctr_atom_deriv, ctr_atom_deriv);
934  Detail::addVectors(oop_atom_deriv, ctr_atom_deriv, ctr_atom_deriv);
935  Detail::negateVector(ctr_atom_deriv);
936 
937  return ang_cos;
938 }
939 
940 // \endcond
941 
942 #endif // CDPL_FORCEFIELD_UTILITYFUNCTIONS_HPP
BitSet.hpp
Declaration of type CDPL::Util::BitSet.
APIPrefix.hpp
Definition of the preprocessor macro CDPL_FORCEFIELD_API.
CDPL_FORCEFIELD_API
#define CDPL_FORCEFIELD_API
Tells the compiler/linker which classes, functions and variables are part of the library API.
CDPL::ForceField::MMFF94InteractionData
Container holding the full set of MMFF94 interaction parameter records for a molecular graph.
Definition: MMFF94InteractionData.hpp:62
CDPL::Chem::AtomType::T
constexpr unsigned int T
Specifies Hydrogen (Tritium).
Definition: AtomType.hpp:67
CDPL::ForceField::calcOutOfPlaneAngle
ValueType calcOutOfPlaneAngle(const CoordsVec &term_atom1_pos, const CoordsVec &ctr_atom_pos, const CoordsVec &term_atom2_pos, const CoordsVec &oop_atom_pos)
Calculates the out-of-plane angle between the bond j-l and the plane defined by the atoms i-j-k.
CDPL::ForceField::calcDihedralAngleCosDerivatives
ValueType calcDihedralAngleCosDerivatives(const CoordsVec &term_atom1_pos, const CoordsVec &ctr_atom1_pos, const CoordsVec &ctr_atom2_pos, const CoordsVec &term_atom2_pos, GradVec &term_atom1_deriv, GradVec &ctr_atom1_deriv, GradVec &ctr_atom2_deriv, GradVec &term_atom2_deriv)
Calculates the partial derivatives of the cosine of the angle between the planes defined by the ato...
CDPL::ForceField::calcSquaredDistance
ValueType calcSquaredDistance(const CoordsVec &atom1_pos, const CoordsVec &atom2_pos)
Calculates the squared distance between two atoms i and j.
CDPL::ForceField::calcBondAngle
ValueType calcBondAngle(const CoordsVec &term_atom1_pos, const CoordsVec &ctr_atom_pos, const CoordsVec &term_atom2_pos)
Calculates the bond angle between the two bonds i-j and j-k.
CDPL::ForceField::calcBondLengthsAndAngleCos
ValueType calcBondLengthsAndAngleCos(const CoordsVec &term_atom1_pos, const CoordsVec &ctr_atom_pos, const CoordsVec &term_atom2_pos, ValueType &bond_length1, ValueType &bond_length2)
Calculates bond lengths and and the cosine of the bond angle between the two bonds i-j and j-k.
CDPL::ForceField::calcOutOfPlaneAngleCosDerivatives
ValueType calcOutOfPlaneAngleCosDerivatives(const CoordsVec &term_atom1_pos, const CoordsVec &ctr_atom_pos, const CoordsVec &term_atom2_pos, const CoordsVec &oop_atom_pos, GradVec &term_atom1_deriv, GradVec &ctr_atom_deriv, GradVec &term_atom2_deriv, GradVec &oop_atom_deriv)
Calculates the partial derivatives of the cosine of the angle between the bond j-l and the normal o...
CDPL::ForceField::calcDihedralAngleCos
ValueType calcDihedralAngleCos(const CoordsVec &term_atom1_pos, const CoordsVec &ctr_atom1_pos, const CoordsVec &ctr_atom2_pos, const CoordsVec &term_atom2_pos)
Calculates the cosine of the dihedral angle between the planes defined by the atom triplets i-j-k an...
CDPL::ForceField::calcBondLengthsAndAngle
ValueType calcBondLengthsAndAngle(const CoordsVec &term_atom1_pos, const CoordsVec &ctr_atom_pos, const CoordsVec &term_atom2_pos, ValueType &bond_length1, ValueType &bond_length2)
Calculates bond lengths and and the bond angle between the two bonds i-j and j-k.
CDPL::ForceField::calcBondAngleCosDerivatives
ValueType calcBondAngleCosDerivatives(const CoordsVec &term_atom1_pos, const CoordsVec &ctr_atom_pos, const CoordsVec &term_atom2_pos, GradVec &term_atom1_deriv, GradVec &ctr_atom_deriv, GradVec &term_atom2_deriv)
Calculates the partial derivatives of the of the cosine of the angle between the bonds i-j and j-k.
CDPL::ForceField::calcDistance
ValueType calcDistance(const CoordsVec &atom1_pos, const CoordsVec &atom2_pos)
Calculates the distance between two atoms i and j.
CDPL::ForceField::filterInteractions
CDPL_FORCEFIELD_API void filterInteractions(const MMFF94InteractionData &ia_data, MMFF94InteractionData &filtered_ia_data, const Util::BitSet &inc_atom_mask)
Filters an MMFF94 interaction data set, retaining only those interactions that exclusively reference ...
CDPL::ForceField::calcBondAngleCos
ValueType calcBondAngleCos(const CoordsVec &term_atom1_pos, const CoordsVec &ctr_atom_pos, const CoordsVec &term_atom2_pos)
Calculates the cosine of the bond angle between the two bonds i-j and j-k.
CDPL::ForceField::calcDistanceDerivatives
ValueType calcDistanceDerivatives(const CoordsVec &atom1_pos, const CoordsVec &atom2_pos, GradVec &atom1_deriv, GradVec &atom2_deriv)
Calculates the partial derivatives of the distance between two atoms i and j.
CDPL::Math::vec
QuaternionVectorAdapter< E > vec(QuaternionExpression< E > &e)
Creates a mutable Math::QuaternionVectorAdapter view of the quaternion expression e.
Definition: QuaternionAdapter.hpp:404
CDPL::Util::BitSet
boost::dynamic_bitset BitSet
Dynamic bitset class.
Definition: BitSet.hpp:46
CDPL
The namespace of the Chemical Data Processing Library.
Generated by doxygen 1.9.1