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