Chemical Data Processing Library C++ API - Version 1.2.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 
604  template <typename ValueType, typename Iter, typename CoordsArray, typename FuncType>
605  ValueType accumInteractionEnergies(Iter& beg, const Iter& end, const CoordsArray& coords, const FuncType& func)
606  {
607  ValueType e = ValueType();
608 
609  for (; beg != end; ++beg)
610  e += func(*beg, coords);
611 
612  return e;
613  }
614 
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)
617  {
618  ValueType e = ValueType();
619 
620  for (; beg != end; ++beg)
621  e += func(*beg, coords, grad);
622 
623  return e;
624  }
625  } // namespace Detail
626  } // namespace ForceField
627 } // namespace CDPL
628 
629 
630 template <typename ValueType, typename CoordsVec>
631 ValueType CDPL::ForceField::calcSquaredDistance(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos)
632 {
633  ValueType pos_diff[3];
634 
635  Detail::subVectors(atom1_pos, atom2_pos, pos_diff);
636 
637  return Detail::calcDotProduct<ValueType>(pos_diff, pos_diff);
638 }
639 
640 template <typename ValueType, typename CoordsVec>
641 ValueType CDPL::ForceField::calcDistance(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos)
642 {
643  return std::sqrt(calcSquaredDistance<ValueType>(atom1_pos, atom2_pos));
644 }
645 
646 template <typename ValueType, typename CoordsVec>
647 ValueType CDPL::ForceField::calcBondLengthsAndAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
648  ValueType& bond_length1, ValueType& bond_length2)
649 {
650  ValueType bond_vec1[3];
651  ValueType bond_vec2[3];
652 
653  Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
654  Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
655 
656  bond_length1 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec1));
657  bond_length2 = std::sqrt(Detail::calcDotProduct<ValueType>(bond_vec2, bond_vec2));
658 
659  return Detail::clampCosine(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2) / (bond_length1 * bond_length2));
660 }
661 
662 template <typename ValueType, typename CoordsVec>
663 ValueType CDPL::ForceField::calcBondLengthsAndAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
664  ValueType& bond_length1, ValueType& bond_length2)
665 {
666  return std::acos(calcBondLengthsAndAngleCos(term_atom1_pos, ctr_atom_pos, term_atom2_pos, bond_length1, bond_length2));
667 }
668 
669 template <typename ValueType, typename CoordsVec>
670 ValueType CDPL::ForceField::calcBondAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos)
671 {
672  ValueType bl1, bl2;
673 
674  return calcBondLengthsAndAngleCos(term_atom1_pos, ctr_atom_pos, term_atom2_pos, bl1, bl2);
675 }
676 
677 template <typename ValueType, typename CoordsVec>
678 ValueType CDPL::ForceField::calcBondAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
679  const ValueType& r_ij, const ValueType& r_jk)
680 {
681  ValueType bond_vec1[3];
682  ValueType bond_vec2[3];
683 
684  Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
685  Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
686 
687  return Detail::clampCosine(Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2) / (r_ij * r_jk));
688 }
689 
690 template <typename ValueType, typename CoordsVec>
691 ValueType CDPL::ForceField::calcBondAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos)
692 {
693  return std::acos(calcBondAngleCos<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos));
694 }
695 
696 template <typename ValueType, typename CoordsVec>
697 ValueType CDPL::ForceField::calcBondAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
698  const ValueType& r_ij, const ValueType& r_jk)
699 {
700  return std::acos(calcBondAngleCos<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, r_ij, r_jk));
701 }
702 
703 template <typename ValueType, typename CoordsVec>
704 ValueType CDPL::ForceField::calcOutOfPlaneAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
705  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos)
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 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));
720 
721  return (ValueType(M_PI * 0.5) - std::acos(ang_cos));
722 }
723 
724 template <typename ValueType, typename CoordsVec>
725 ValueType CDPL::ForceField::calcOutOfPlaneAngle(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
726  const CoordsVec& term_atom2_pos, const CoordsVec& oop_atom_pos,
727  const ValueType& r_jl)
728 {
729  ValueType term_bond1_vec[3];
730  ValueType term_bond2_vec[3];
731  ValueType oop_bond_vec[3];
732  ValueType plane_normal[3];
733 
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);
738 
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));
741 
742  return (ValueType(M_PI * 0.5) - std::acos(ang_cos));
743 }
744 
745 template <typename ValueType, typename CoordsVec>
746 ValueType CDPL::ForceField::calcDihedralAngleCos(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos,
747  const CoordsVec& ctr_atom2_pos, const CoordsVec& term_atom2_pos)
748 {
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];
754 
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);
758 
759  Detail::calcCrossProduct(term_bond1_vec, ctr_bond_vec, plane_normal1);
760  Detail::calcCrossProduct(ctr_bond_vec, term_bond2_vec, plane_normal2);
761 
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));
764 
765  return Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal2) / (pn1_len * pn2_len));
766 }
767 
768 template <typename ValueType, typename CoordsVec, typename GradVec>
769 ValueType CDPL::ForceField::calcDistanceDerivatives(const CoordsVec& atom1_pos, const CoordsVec& atom2_pos,
770  GradVec& atom1_deriv, GradVec& atom2_deriv)
771 {
772  Detail::subVectors(atom1_pos, atom2_pos, atom1_deriv);
773 
774  ValueType dist = std::sqrt(Detail::calcDotProduct<ValueType>(atom1_deriv, atom1_deriv));
775 
776  Detail::invScaleVector(atom1_deriv, -dist);
777  Detail::negateCopyVector(atom1_deriv, atom2_deriv);
778 
779  return dist;
780 }
781 
782 template <typename ValueType, typename CoordsVec, typename GradVec>
783 ValueType CDPL::ForceField::calcBondAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos, const CoordsVec& term_atom2_pos,
784  GradVec& term_atom1_deriv, GradVec& ctr_atom_deriv, GradVec& term_atom2_deriv)
785 {
786  ValueType bond_vec1[3];
787  ValueType bond_vec2[3];
788 
789  Detail::subVectors(ctr_atom_pos, term_atom1_pos, bond_vec1);
790  Detail::subVectors(ctr_atom_pos, term_atom2_pos, bond_vec2);
791 
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));
794 
795  ValueType dot_prod = Detail::calcDotProduct<ValueType>(bond_vec1, bond_vec2);
796  ValueType bl_prod = bond_length1 * bond_length2;
797 
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);
801 
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);
805 
806  Detail::negateCopyVector(term_atom1_deriv, ctr_atom_deriv);
807  Detail::subVectors(term_atom2_deriv, ctr_atom_deriv, ctr_atom_deriv);
808 
809  return Detail::clampCosine(dot_prod / bl_prod);
810 }
811 
812 template <typename ValueType, typename CoordsVec, typename GradVec>
813 ValueType CDPL::ForceField::calcDihedralAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom1_pos,
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)
817 {
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];
824 
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);
829 
830  Detail::calcCrossProduct(term_bond1_vec, ctr_bond_vec, plane_normal1);
831  Detail::calcCrossProduct(ctr_bond_vec, term_bond2_vec, plane_normal2);
832 
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));
835 
836  Detail::invScaleVector(plane_normal1, pn1_len);
837  Detail::invScaleVector(plane_normal2, pn2_len);
838 
839  ValueType ang_cos = Detail::clampCosine(Detail::calcDotProduct<ValueType>(plane_normal1, plane_normal2));
840 
841  ValueType a[3];
842  ValueType b[3];
843 
844  Detail::scaleCopyVector(plane_normal1, -ang_cos, a);
845  Detail::addVectors(a, plane_normal2, a);
846  Detail::invScaleVector(a, pn1_len);
847 
848  Detail::scaleCopyVector(plane_normal2, -ang_cos, b);
849  Detail::addVectors(b, plane_normal1, b);
850  Detail::invScaleVector(b, pn2_len);
851 
852  Detail::calcCrossProduct(ctr_bond_vec, a, term_atom1_deriv);
853  Detail::calcCrossProduct(ctr_bond_vec, b, term_atom2_deriv);
854 
855  Detail::calcCrossProduct(term1_ctr2_vec, a, ctr_atom1_deriv);
856  Detail::calcCrossProduct(term_bond2_vec, b, ctr_atom2_deriv);
857 
858  Detail::subVectors(ctr_atom2_deriv, ctr_atom1_deriv, ctr_atom1_deriv);
859 
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);
863 
864  return ang_cos;
865 }
866 
867 template <typename ValueType, typename CoordsVec, typename GradVec>
868 ValueType CDPL::ForceField::calcOutOfPlaneAngleCosDerivatives(const CoordsVec& term_atom1_pos, const CoordsVec& ctr_atom_pos,
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)
872 {
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];
878  ValueType ijk_pn[3];
879 
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);
885 
886  Detail::calcCrossProduct(term_bond1_vec, term_bond2_vec, ijk_pn);
887 
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);
894 
895  ValueType kjl_pn[3];
896  ValueType lji_pn[3];
897 
898  Detail::calcCrossProduct(term_bond2_vec, oop_bond_vec, kjl_pn);
899  Detail::calcCrossProduct(oop_bond_vec, term_bond1_vec, lji_pn);
900 
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];
904 
905  Detail::scaleVector(ctr_atom_deriv, ang_cos / ijk_pn_len_2);
906 
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);
909 
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];
913 
914  Detail::scaleVector(ctr_atom_deriv, ang_cos / ijk_pn_len_2);
915 
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);
918 
919  Detail::calcCrossProduct(term2_oop_vec, term1_oop_vec, ctr_atom_deriv);
920 
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);
926 
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);
931 
932  return ang_cos;
933 }
934 
935 // \endcond
936 
937 #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