Chemical Data Processing Library C++ API - Version 1.2.1
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 
46  CDPL_FORCEFIELD_API void filterInteractions(const MMFF94InteractionData& ia_data, MMFF94InteractionData& filtered_ia_data, const Util::BitSet& inc_atom_mask);
47 
63  template <typename ValueType, typename CoordsVec>
64  ValueType calcSquaredDistance(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos);
65 
81  template <typename ValueType, typename CoordsVec>
82  ValueType calcDistance(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos);
83 
107  template <typename ValueType, typename CoordsVec>
108  ValueType calcBondLengthsAndAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
109  ValueType& bond_length1, ValueType& bond_length2);
110 
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);
137 
156  template <typename ValueType, typename CoordsVec>
157  ValueType calcBondAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos);
158 
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);
182 
201  template <typename ValueType, typename CoordsVec>
202  ValueType calcBondAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos);
203 
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);
227 
250  template <typename ValueType, typename CoordsVec>
251  ValueType calcOutOfPlaneAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
252  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos);
253 
277  template <typename ValueType, typename CoordsVec>
278  ValueType calcOutOfPlaneAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
279  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos,
280  const ValueType& r_jl);
281 
305  template <typename ValueType, typename CoordsVec>
306  ValueType calcDihedralAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos,
307  const CoordsVec& ctr_atom2_pos, const CoordsVec& term_atom2_pos);
308 
331  template <typename ValueType, typename CoordsVec, typename GradVec>
332  ValueType calcDistanceDerivatives(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos,
333  GradVec& atom1_deriv, GradVec& atom2_deriv);
334 
365  template <typename ValueType, typename CoordsVec, typename GradVec>
366  ValueType calcBondAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
367  GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv, GradVec& term_atom2_deriv);
368 
410  template <typename ValueType, typename CoordsVec, typename GradVec>
411  ValueType calcDihedralAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos,
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);
415 
475  template <typename ValueType, typename CoordsVec, typename GradVec>
476  ValueType calcOutOfPlaneAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
477  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos,
478  GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv,
479  GradVec& term_atom2_deriv, GradVec& oop_atom_deriv);
480  } // namespace ForceField
481 } // namespace CDPL
482 
483 
484 // Implementation
485 // \cond DOC_IMPL_DETAILS
486 
487 namespace CDPL
488 {
489 
490  namespace ForceField
491  {
492 
493  namespace Detail
494  {
495 
496  template <typename VecType1, typename VecType2, typename VecType3>
497  void addVectors(const VecType1& vec1, const VecType2& vec2, VecType3& res)
498  {
499  res[0] = vec2[0] + vec1[0];
500  res[1] = vec2[1] + vec1[1];
501  res[2] = vec2[2] + vec1[2];
502  }
503 
504  template <typename VecType1, typename VecType2>
505  void subVectors(const VecType1& vec1, const VecType1& vec2, VecType2& res)
506  {
507  res[0] = vec2[0] - vec1[0];
508  res[1] = vec2[1] - vec1[1];
509  res[2] = vec2[2] - vec1[2];
510  }
511 
512  template <typename ValueType, typename VecType>
513  ValueType calcDotProduct(const VecType& vec1, const VecType& vec2)
514  {
515  return (vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]);
516  }
517 
518  template <typename VecType, typename ResVecType>
519  void calcCrossProduct(const VecType& vec1, const VecType& vec2, ResVecType& cross_prod)
520  {
521  cross_prod[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
522  cross_prod[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
523  cross_prod[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
524  }
525 
526  template <typename VecType1, typename VecType2>
527  void copyVector(const VecType1& vec1, VecType2& vec2)
528  {
529  vec2[0] = vec1[0];
530  vec2[1] = vec1[1];
531  vec2[2] = vec1[2];
532  }
533 
534  template <typename VecType>
535  void negateVector(VecType& vec)
536  {
537  vec[0] = -vec[0];
538  vec[1] = -vec[1];
539  vec[2] = -vec[2];
540  }
541 
542  template <typename VecType1, typename VecType2>
543  void negateCopyVector(const VecType1& vec1, VecType2& vec2)
544  {
545  vec2[0] = -vec1[0];
546  vec2[1] = -vec1[1];
547  vec2[2] = -vec1[2];
548  }
549 
550  template <typename VecType, typename T>
551  void scaleVector(VecType& vec, const T& factor)
552  {
553  vec[0] *= factor;
554  vec[1] *= factor;
555  vec[2] *= factor;
556  }
557 
558  template <typename VecType, typename T>
559  void invScaleVector(VecType& vec, const T& factor)
560  {
561  vec[0] /= factor;
562  vec[1] /= factor;
563  vec[2] /= factor;
564  }
565 
566  template <typename VecType1, typename VecType2, typename T>
567  void scaleAddVector(const VecType1& vec1, const T& factor, VecType2& vec2)
568  {
569  vec2[0] += vec1[0] * factor;
570  vec2[1] += vec1[1] * factor;
571  vec2[2] += vec1[2] * factor;
572  }
573 
574  template <typename VecType1, typename VecType2, typename T>
575  void scaleCopyVector(const VecType1& vec1, const T& factor, VecType2& vec2)
576  {
577  vec2[0] = vec1[0] * factor;
578  vec2[1] = vec1[1] * factor;
579  vec2[2] = vec1[2] * factor;
580  }
581 
582  template <typename VecType1, typename VecType2, typename T>
583  void invScaleCopyVector(const VecType1& vec1, const T& factor, VecType2& vec2)
584  {
585  vec2[0] = vec1[0] / factor;
586  vec2[1] = vec1[1] / factor;
587  vec2[2] = vec1[2] / factor;
588  }
589 
590  template <typename ValueType>
591  ValueType clampCosine(const ValueType& v)
592  {
593  if (v > ValueType(1))
594  return ValueType(1);
595 
596  if (v < ValueType(-1))
597  return ValueType(-1);
598 
599  return v;
600  }
601 
602  template <typename ValueType, typename Iter, typename CoordsArray, typename FuncType>
603  ValueType accumInteractionEnergies(Iter& beg, const Iter& end, const CoordsArray& coords, const FuncType& func)
604  {
605  ValueType e = ValueType();
606 
607  for (; beg != end; ++beg)
608  e += func(*beg, coords);
609 
610  return e;
611  }
612 
613  template <typename ValueType, typename Iter, typename CoordsArray, typename GradVector, typename FuncType>
614  ValueType calcInteractionGradient(Iter& beg, const Iter& end, const CoordsArray& coords, GradVector& grad, const FuncType& func)
615  {
616  ValueType e = ValueType();
617 
618  for (; beg != end; ++beg)
619  e += func(*beg, coords, grad);
620 
621  return e;
622  }
623  } // namespace Detail
624  } // namespace ForceField
625 } // namespace CDPL
626 
627 
628 template <typename ValueType, typename CoordsVec>
629 ValueType CDPL::ForceField::calcSquaredDistance(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos)
630 {
631  ValueType pos_diff[3];
632 
633  Detail::subVectors(atom1_pos, atom2_pos, pos_diff);
634 
635  return Detail::calcDotProduct<ValueType>(pos_diff, pos_diff);
636 }
637 
638 template <typename ValueType, typename CoordsVec>
639 ValueType CDPL::ForceField::calcDistance(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos)
640 {
641  return std::sqrt(calcSquaredDistance<ValueType>(atom1_pos, atom2_pos));
642 }
643 
644 template <typename ValueType, typename CoordsVec>
645 ValueType CDPL::ForceField::calcBondLengthsAndAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
646  ValueType& bond_length1, ValueType& bond_length2)
647 {
648  ValueType bond_vec1[3];
649  ValueType bond_vec2[3];
650 
651  Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
652  Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
653 
654  bond_length1 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec1));
655  bond_length2 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec2, bond_vec2));
656 
657  return Detail::clampCosine(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2) / (bond_length1 * bond_length2));
658 }
659 
660 template <typename ValueType, typename CoordsVec>
661 ValueType CDPL::ForceField::calcBondLengthsAndAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
662  ValueType& bond_length1, ValueType& bond_length2)
663 {
664  return std::acos(calcBondLengthsAndAngleCos(term_atom1_pos, ctr_atom_pos, term_atom2_pos, bond_length1, bond_length2));
665 }
666 
667 template <typename ValueType, typename CoordsVec>
668 ValueType CDPL::ForceField::calcBondAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos)
669 {
670  ValueType bl1, bl2;
671 
672  return calcBondLengthsAndAngleCos(term_atom1_pos, ctr_atom_pos, term_atom2_pos, bl1, bl2);
673 }
674 
675 template <typename ValueType, typename CoordsVec>
676 ValueType CDPL::ForceField::calcBondAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
677  const ValueType& r_ij, const ValueType& r_jk)
678 {
679  ValueType bond_vec1[3];
680  ValueType bond_vec2[3];
681 
682  Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
683  Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
684 
685  return Detail::clampCosine(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2) / (r_ij * r_jk));
686 }
687 
688 template <typename ValueType, typename CoordsVec>
689 ValueType CDPL::ForceField::calcBondAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos)
690 {
691  return std::acos(calcBondAngleCos<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos));
692 }
693 
694 template <typename ValueType, typename CoordsVec>
695 ValueType CDPL::ForceField::calcBondAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
696  const ValueType& r_ij, const ValueType& r_jk)
697 {
698  return std::acos(calcBondAngleCos<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, r_ij, r_jk));
699 }
700 
701 template <typename ValueType, typename CoordsVec>
702 ValueType CDPL::ForceField::calcOutOfPlaneAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
703  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos)
704 {
705  ValueType term_bond1_vec[3];
706  ValueType term_bond2_vec[3];
707  ValueType oop_bond_vec[3];
708  ValueType plane_normal[3];
709 
710  Detail::subVectors(ctr_atom_pos, term_atom1_pos, term_bond1_vec);
711  Detail::subVectors(ctr_atom_pos, term_atom2_pos, term_bond2_vec);
712  Detail::subVectors(ctr_atom_pos, oop_atom_pos, oop_bond_vec);
713  Detail::calcCrossProduct(term_bond1_vec, term_bond2_vec, plane_normal);
714 
715  ValueType pn_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal, plane_normal));
716  ValueType oop_bnd_len = std::sqrt(Detail::calcDotProduct<ValueType>(oop_bond_vec, oop_bond_vec));
717  ValueType ang_cos = Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal, oop_bond_vec) / (pn_len * oop_bnd_len));
718 
719  return (ValueType(M_PI * 0.5) - std::acos(ang_cos));
720 }
721 
722 template <typename ValueType, typename CoordsVec>
723 ValueType CDPL::ForceField::calcOutOfPlaneAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
724  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos,
725  const ValueType& r_jl)
726 {
727  ValueType term_bond1_vec[3];
728  ValueType term_bond2_vec[3];
729  ValueType oop_bond_vec[3];
730  ValueType plane_normal[3];
731 
732  Detail::subVectors(ctr_atom_pos, term_atom1_pos, term_bond1_vec);
733  Detail::subVectors(ctr_atom_pos, term_atom2_pos, term_bond2_vec);
734  Detail::subVectors(ctr_atom_pos, oop_atom_pos, oop_bond_vec);
735  Detail::calcCrossProduct(term_bond1_vec, term_bond2_vec, plane_normal);
736 
737  ValueType pn_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal, plane_normal));
738  ValueType ang_cos = Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal, oop_bond_vec) / (pn_len * r_jl));
739 
740  return (ValueType(M_PI * 0.5) - std::acos(ang_cos));
741 }
742 
743 template <typename ValueType, typename CoordsVec>
744 ValueType CDPL::ForceField::calcDihedralAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos,
745  const CoordsVec& ctr_atom2_pos, const CoordsVec& term_atom2_pos)
746 {
747  ValueType term_bond1_vec[3];
748  ValueType ctr_bond_vec[3];
749  ValueType term_bond2_vec[3];
750  ValueType plane_normal1[3];
751  ValueType plane_normal2[3];
752 
753  Detail::subVectors(ctr_atom1_pos, term_atom1_pos, term_bond1_vec);
754  Detail::subVectors(ctr_atom1_pos, ctr_atom2_pos, ctr_bond_vec);
755  Detail::subVectors(term_atom2_pos, ctr_atom2_pos, term_bond2_vec);
756 
757  Detail::calcCrossProduct(term_bond1_vec, ctr_bond_vec, plane_normal1);
758  Detail::calcCrossProduct(ctr_bond_vec, term_bond2_vec, plane_normal2);
759 
760  ValueType pn1_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal1));
761  ValueType pn2_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal2, plane_normal2));
762 
763  return Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal2) / (pn1_len * pn2_len));
764 }
765 
766 template <typename ValueType, typename CoordsVec, typename GradVec>
767 ValueType CDPL::ForceField::calcDistanceDerivatives(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos,
768  GradVec& atom1_deriv, GradVec& atom2_deriv)
769 {
770  Detail::subVectors(atom1_pos, atom2_pos, atom1_deriv);
771 
772  ValueType dist = std::sqrt(Detail::calcDotProduct<ValueType>(atom1_deriv, atom1_deriv));
773 
774  Detail::invScaleVector(atom1_deriv, -dist);
775  Detail::negateCopyVector(atom1_deriv, atom2_deriv);
776 
777  return dist;
778 }
779 
780 template <typename ValueType, typename CoordsVec, typename GradVec>
781 ValueType CDPL::ForceField::calcBondAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
782  GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv, GradVec& term_atom2_deriv)
783 {
784  ValueType bond_vec1[3];
785  ValueType bond_vec2[3];
786 
787  Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
788  Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
789 
790  ValueType bond_length1 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec1));
791  ValueType bond_length2 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec2, bond_vec2));
792 
793  ValueType dot_prod = Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2);
794  ValueType bl_prod = bond_length1 * bond_length2;
795 
796  Detail::invScaleCopyVector(bond_vec2, bl_prod, term_atom1_deriv);
797  Detail::scaleCopyVector(bond_vec1, dot_prod / (bond_length1 * bond_length1 * bl_prod), ctr_atom_deriv);
798  Detail::subVectors(ctr_atom_deriv, term_atom1_deriv, term_atom1_deriv);
799 
800  Detail::invScaleCopyVector(bond_vec1, bl_prod, term_atom2_deriv);
801  Detail::scaleCopyVector(bond_vec2, dot_prod / (bond_length2 * bond_length2 * bl_prod), ctr_atom_deriv);
802  Detail::subVectors(ctr_atom_deriv, term_atom2_deriv, term_atom2_deriv);
803 
804  Detail::negateCopyVector(term_atom1_deriv, ctr_atom_deriv);
805  Detail::subVectors(term_atom2_deriv, ctr_atom_deriv, ctr_atom_deriv);
806 
807  return Detail::clampCosine(dot_prod / bl_prod);
808 }
809 
810 template <typename ValueType, typename CoordsVec, typename GradVec>
811 ValueType CDPL::ForceField::calcDihedralAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos,
812  const CoordsVec& ctr_atom2_pos, const CoordsVec& term_atom2_pos,
813  GradVec& term_atom1_deriv, GradVec& ctr_atom1_deriv,
814  GradVec& ctr_atom2_deriv, GradVec& term_atom2_deriv)
815 {
816  ValueType term_bond1_vec[3];
817  ValueType ctr_bond_vec[3];
818  ValueType term_bond2_vec[3];
819  ValueType term1_ctr2_vec[3];
820  ValueType plane_normal1[3];
821  ValueType plane_normal2[3];
822 
823  Detail::subVectors(ctr_atom1_pos, term_atom1_pos, term_bond1_vec);
824  Detail::subVectors(ctr_atom1_pos, ctr_atom2_pos, ctr_bond_vec);
825  Detail::subVectors(term_atom2_pos, ctr_atom2_pos, term_bond2_vec);
826  Detail::subVectors(ctr_atom2_pos, term_atom1_pos, term1_ctr2_vec);
827 
828  Detail::calcCrossProduct(term_bond1_vec, ctr_bond_vec, plane_normal1);
829  Detail::calcCrossProduct(ctr_bond_vec, term_bond2_vec, plane_normal2);
830 
831  ValueType pn1_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal1));
832  ValueType pn2_len = std::sqrt(Detail::calcDotProduct<ValueType>(plane_normal2, plane_normal2));
833 
834  Detail::invScaleVector(plane_normal1, pn1_len);
835  Detail::invScaleVector(plane_normal2, pn2_len);
836 
837  ValueType ang_cos = Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal2));
838 
839  ValueType a[3];
840  ValueType b[3];
841 
842  Detail::scaleCopyVector(plane_normal1, -ang_cos, a);
843  Detail::addVectors(a, plane_normal2, a);
844  Detail::invScaleVector(a, pn1_len);
845 
846  Detail::scaleCopyVector(plane_normal2, -ang_cos, b);
847  Detail::addVectors(b, plane_normal1, b);
848  Detail::invScaleVector(b, pn2_len);
849 
850  Detail::calcCrossProduct(ctr_bond_vec, a, term_atom1_deriv);
851  Detail::calcCrossProduct(ctr_bond_vec, b, term_atom2_deriv);
852 
853  Detail::calcCrossProduct(term1_ctr2_vec, a, ctr_atom1_deriv);
854  Detail::calcCrossProduct(term_bond2_vec, b, ctr_atom2_deriv);
855 
856  Detail::subVectors(ctr_atom2_deriv, ctr_atom1_deriv, ctr_atom1_deriv);
857 
858  Detail::addVectors(term_atom1_deriv, ctr_atom1_deriv, ctr_atom2_deriv);
859  Detail::addVectors(ctr_atom2_deriv, term_atom2_deriv, ctr_atom2_deriv);
860  Detail::negateVector(ctr_atom2_deriv);
861 
862  return ang_cos;
863 }
864 
865 template <typename ValueType, typename CoordsVec, typename GradVec>
866 ValueType CDPL::ForceField::calcOutOfPlaneAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
867  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos,
868  GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv,
869  GradVec& term_atom2_deriv, GradVec& oop_atom_deriv)
870 {
871  ValueType term_bond1_vec[3];
872  ValueType term_bond2_vec[3];
873  ValueType oop_bond_vec[3];
874  ValueType term1_oop_vec[3];
875  ValueType term2_oop_vec[3];
876  ValueType ijk_pn[3];
877 
878  Detail::subVectors(ctr_atom_pos, term_atom1_pos, term_bond1_vec);
879  Detail::subVectors(ctr_atom_pos, term_atom2_pos, term_bond2_vec);
880  Detail::subVectors(ctr_atom_pos, oop_atom_pos, oop_bond_vec);
881  Detail::subVectors(term_atom2_pos, oop_atom_pos, term2_oop_vec);
882  Detail::subVectors(term_atom1_pos, oop_atom_pos, term1_oop_vec);
883 
884  Detail::calcCrossProduct(term_bond1_vec, term_bond2_vec, ijk_pn);
885 
886  ValueType ijk_pn_len_2 = Detail::calcDotProduct<ValueType>(ijk_pn, ijk_pn);
887  ValueType ijk_pn_len = std::sqrt(ijk_pn_len_2);
888  ValueType oop_bnd_len = std::sqrt(Detail::calcDotProduct<ValueType>(oop_bond_vec, oop_bond_vec));
889  ValueType oop_ijk_pn_len_prod = ijk_pn_len * oop_bnd_len;
890  ValueType oop_ijk_pn_dot_prod = Detail::calcDotProduct<ValueType>(ijk_pn, oop_bond_vec);
891  ValueType ang_cos = Detail::clampCosine(oop_ijk_pn_dot_prod / oop_ijk_pn_len_prod);
892 
893  ValueType kjl_pn[3];
894  ValueType lji_pn[3];
895 
896  Detail::calcCrossProduct(term_bond2_vec, oop_bond_vec, kjl_pn);
897  Detail::calcCrossProduct(oop_bond_vec, term_bond1_vec, lji_pn);
898 
899  ctr_atom_deriv[0] = ijk_pn[1] * -term_bond2_vec[2] + ijk_pn[2] * term_bond2_vec[1];
900  ctr_atom_deriv[1] = ijk_pn[0] * term_bond2_vec[2] + ijk_pn[2] * -term_bond2_vec[0];
901  ctr_atom_deriv[2] = ijk_pn[0] * -term_bond2_vec[1] + ijk_pn[1] * term_bond2_vec[0];
902 
903  Detail::scaleVector(ctr_atom_deriv, ang_cos / ijk_pn_len_2);
904 
905  Detail::invScaleCopyVector(kjl_pn, oop_ijk_pn_len_prod, term_atom1_deriv);
906  Detail::subVectors(ctr_atom_deriv, term_atom1_deriv, term_atom1_deriv);
907 
908  ctr_atom_deriv[0] = ijk_pn[1] * term_bond1_vec[2] + ijk_pn[2] * -term_bond1_vec[1];
909  ctr_atom_deriv[1] = ijk_pn[0] * -term_bond1_vec[2] + ijk_pn[2] * term_bond1_vec[0];
910  ctr_atom_deriv[2] = ijk_pn[0] * term_bond1_vec[1] + ijk_pn[1] * -term_bond1_vec[0];
911 
912  Detail::scaleVector(ctr_atom_deriv, ang_cos / ijk_pn_len_2);
913 
914  Detail::invScaleCopyVector(lji_pn, oop_ijk_pn_len_prod, term_atom2_deriv);
915  Detail::subVectors(ctr_atom_deriv, term_atom2_deriv, term_atom2_deriv);
916 
917  Detail::calcCrossProduct(term2_oop_vec, term1_oop_vec, ctr_atom_deriv);
918 
919  Detail::scaleCopyVector(oop_bond_vec, oop_ijk_pn_dot_prod / (oop_bnd_len * oop_bnd_len), oop_atom_deriv);
920  Detail::addVectors(oop_atom_deriv, lji_pn, oop_atom_deriv);
921  Detail::addVectors(oop_atom_deriv, kjl_pn, oop_atom_deriv);
922  Detail::addVectors(oop_atom_deriv, ctr_atom_deriv, oop_atom_deriv);
923  Detail::invScaleVector(oop_atom_deriv, -oop_ijk_pn_len_prod);
924 
925  Detail::copyVector(term_atom1_deriv, ctr_atom_deriv);
926  Detail::addVectors(term_atom2_deriv, ctr_atom_deriv, ctr_atom_deriv);
927  Detail::addVectors(oop_atom_deriv, ctr_atom_deriv, ctr_atom_deriv);
928  Detail::negateVector(ctr_atom_deriv);
929 
930  return ang_cos;
931 }
932 
933 // \endcond
934 
935 #endif // CDPL_FORCEFIELD_UTILITYFUNCTIONS_HPP
BitSet.hpp
Definition of the 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
Definition: MMFF94InteractionData.hpp:51
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)
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)
Definition: QuaternionAdapter.hpp:237
CDPL::Util::BitSet
boost::dynamic_bitset BitSet
A dynamic bitset class.
Definition: BitSet.hpp:46
CDPL
The namespace of the Chemical Data Processing Library.
Generated by doxygen 1.9.1