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