29 #ifndef CDPL_FORCEFIELD_MMFF94ENERGYFUNCTIONS_HPP
30 #define CDPL_FORCEFIELD_MMFF94ENERGYFUNCTIONS_HPP
48 template <
typename ValueType,
typename Iter,
typename CoordsArray>
51 template <
typename ValueType,
typename CoordsArray>
77 template <
typename ValueType,
typename CoordsVec>
79 const ValueType& force_const,
const ValueType& ref_length);
103 template <
typename ValueType>
107 template <
typename ValueType,
typename Iter,
typename CoordsArray>
110 template <
typename ValueType,
typename CoordsArray>
141 template <
typename ValueType,
typename CoordsVec>
143 bool linear,
const ValueType& force_const,
const ValueType& ref_angle);
175 template <
typename ValueType,
typename CoordsVec>
177 const ValueType& r_ij,
const ValueType& r_jk,
178 bool linear,
const ValueType& force_const,
const ValueType& ref_angle);
181 template <
typename ValueType,
typename Iter,
typename CoordsArray>
184 template <
typename ValueType,
typename CoordsArray>
217 template <
typename ValueType,
typename CoordsVec>
219 const ValueType& ijk_force_const,
const ValueType& kji_force_const,
const ValueType& ref_angle,
220 const ValueType& ref_length1,
const ValueType& ref_length2);
254 template <
typename ValueType,
typename CoordsVec>
256 const ValueType& r_ij,
const ValueType& r_jk,
257 const ValueType& ijk_force_const,
const ValueType& kji_force_const,
const ValueType& ref_angle,
258 const ValueType& ref_length1,
const ValueType& ref_length2);
261 template <
typename ValueType,
typename Iter,
typename CoordsArray>
264 template <
typename ValueType,
typename CoordsArray>
284 template <
typename ValueType,
typename CoordsVec>
286 const CoordsVec& oop_atom_pos,
const ValueType& force_const);
306 template <
typename ValueType,
typename CoordsVec>
308 const CoordsVec& oop_atom_pos,
const ValueType& r_jl,
const ValueType& force_const);
311 template <
typename ValueType,
typename Iter,
typename CoordsArray>
314 template <
typename ValueType,
typename CoordsArray>
338 template <
typename ValueType,
typename CoordsVec>
339 ValueType
calcMMFF94TorsionEnergy(
const CoordsVec& term_atom1_pos,
const CoordsVec& ctr_atom1_pos,
const CoordsVec& ctr_atom2_pos,
340 const CoordsVec& term_atom2_pos,
const ValueType& tor_param1,
const ValueType& tor_param2,
341 const ValueType& tor_param3);
344 template <
typename ValueType,
typename Iter,
typename CoordsArray>
347 template <
typename ValueType,
typename CoordsArray>
376 template <
typename ValueType,
typename CoordsVec>
378 const ValueType& atom2_chg,
const ValueType& scale_fact,
const ValueType& de_const,
379 const ValueType& dist_expo);
406 template <
typename ValueType>
408 const ValueType& scale_fact,
const ValueType& de_const,
const ValueType& dist_expo);
411 template <
typename ValueType,
typename Iter,
typename CoordsArray>
414 template <
typename ValueType,
typename CoordsArray>
452 template <
typename ValueType,
typename CoordsVec>
454 const ValueType& r_IJ,
const ValueType& r_IJ_7);
490 template <
typename ValueType>
492 const ValueType& r_IJ,
const ValueType& r_IJ_7);
500 template <
typename ValueType,
typename Iter,
typename CoordsArray>
503 return Detail::accumInteractionEnergies<ValueType>(beg, end, coords,
504 static_cast<ValueType (*)(
const MMFF94BondStretchingInteraction&,
const CoordsArray&)
>(
505 &calcMMFF94BondStretchingEnergy<ValueType, CoordsArray>));
508 template <
typename ValueType,
typename CoordsArray>
511 return calcMMFF94BondStretchingEnergy<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()],
512 iaction.getForceConstant(), iaction.getReferenceLength());
515 template <
typename ValueType,
typename CoordsVec>
517 const ValueType& force_const,
const ValueType& ref_length)
519 ValueType dr_ij = calcDistance<ValueType>(atom1_pos, atom2_pos) - ref_length;
520 ValueType dr_ij_2 = dr_ij * dr_ij;
521 ValueType e_b = ValueType(143.9325 * 0.5) * force_const * dr_ij_2 * (1 - 2 * dr_ij + 28 * dr_ij_2 / 12);
526 template <
typename ValueType>
529 ValueType dr_ij = r_ij - ref_length;
530 ValueType dr_ij_2 = dr_ij * dr_ij;
531 ValueType e_b = ValueType(143.9325 * 0.5) * force_const * dr_ij_2 * (1 - 2 * dr_ij + 28 * dr_ij_2 / 12);
537 template <
typename ValueType,
typename Iter,
typename CoordsArray>
540 return Detail::accumInteractionEnergies<ValueType>(beg, end, coords,
541 static_cast<ValueType (*)(
const MMFF94AngleBendingInteraction&,
const CoordsArray&)
>(
542 &calcMMFF94AngleBendingEnergy<ValueType, CoordsArray>));
545 template <
typename ValueType,
typename CoordsArray>
548 return calcMMFF94AngleBendingEnergy<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
549 coords[iaction.getTerminalAtom2Index()], iaction.isLinearAngle(), iaction.getForceConstant(),
550 iaction.getReferenceAngle());
553 template <
typename ValueType,
typename CoordsVec>
555 bool linear,
const ValueType& force_const,
const ValueType& ref_angle)
558 return (ValueType(143.9325) * force_const * (1 + calcBondAngleCos<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos)));
560 ValueType a_ijk = calcBondAngle<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos) * ValueType(180 / M_PI);
561 ValueType da_ijk = a_ijk - ref_angle;
562 ValueType e_a = ValueType(0.043844 * 0.5) * force_const * da_ijk * da_ijk * (1 - ValueType(0.007) * da_ijk);
567 template <
typename ValueType,
typename CoordsVec>
569 const ValueType& r_ij,
const ValueType& r_jk,
570 bool linear,
const ValueType& force_const,
const ValueType& ref_angle)
573 return (ValueType(143.9325) * force_const * (1 + calcBondAngleCos<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, r_ij, r_jk)));
575 ValueType a_ijk = calcBondAngle<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, r_ij, r_jk) * ValueType(180 / M_PI);
576 ValueType da_ijk = a_ijk - ref_angle;
577 ValueType e_a = ValueType(0.043844 * 0.5) * force_const * da_ijk * da_ijk * (1 - ValueType(0.007) * da_ijk);
583 template <
typename ValueType,
typename Iter,
typename CoordsArray>
586 return Detail::accumInteractionEnergies<ValueType>(beg, end, coords,
587 static_cast<ValueType (*)(
const MMFF94StretchBendInteraction&,
const CoordsArray&)
>(
588 &calcMMFF94StretchBendEnergy<ValueType, CoordsArray>));
591 template <
typename ValueType,
typename CoordsArray>
594 return calcMMFF94StretchBendEnergy<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
595 coords[iaction.getTerminalAtom2Index()], iaction.getIJKForceConstant(),
596 iaction.getKJIForceConstant(), iaction.getReferenceAngle(), iaction.getReferenceLength1(),
597 iaction.getReferenceLength2());
600 template <
typename ValueType,
typename CoordsVec>
602 const ValueType& ijk_force_const,
const ValueType& kji_force_const,
const ValueType& ref_angle,
603 const ValueType& ref_length1,
const ValueType& ref_length2)
605 ValueType r_ij = ValueType();
606 ValueType r_kj = ValueType();
607 ValueType a_ijk =
calcBondLengthsAndAngle(term_atom1_pos, ctr_atom_pos, term_atom2_pos, r_ij, r_kj) * ValueType(180 / M_PI);
609 ValueType dr_ij = r_ij - ref_length1;
610 ValueType dr_kj = r_kj - ref_length2;
611 ValueType da_ijk = a_ijk - ref_angle;
613 ValueType e_ab = ValueType(2.51210) * (ijk_force_const * dr_ij + kji_force_const * dr_kj) * da_ijk;
618 template <
typename ValueType,
typename CoordsVec>
620 const ValueType& r_ij,
const ValueType& r_jk,
621 const ValueType& ijk_force_const,
const ValueType& kji_force_const,
const ValueType& ref_angle,
622 const ValueType& ref_length1,
const ValueType& ref_length2)
624 ValueType a_ijk =
calcBondAngle(term_atom1_pos, ctr_atom_pos, term_atom2_pos, r_ij, r_jk) * ValueType(180 / M_PI);
626 ValueType dr_ij = r_ij - ref_length1;
627 ValueType dr_kj = r_jk - ref_length2;
628 ValueType da_ijk = a_ijk - ref_angle;
630 ValueType e_ab = ValueType(2.51210) * (ijk_force_const * dr_ij + kji_force_const * dr_kj) * da_ijk;
636 template <
typename ValueType,
typename Iter,
typename CoordsArray>
639 return Detail::accumInteractionEnergies<ValueType>(beg, end, coords,
640 static_cast<ValueType (*)(
const MMFF94OutOfPlaneBendingInteraction&,
const CoordsArray&)
>(
641 &calcMMFF94OutOfPlaneBendingEnergy<ValueType, CoordsArray>));
644 template <
typename ValueType,
typename CoordsArray>
647 return calcMMFF94OutOfPlaneBendingEnergy<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
648 coords[iaction.getTerminalAtom2Index()], coords[iaction.getOutOfPlaneAtomIndex()],
649 iaction.getForceConstant());
652 template <
typename ValueType,
typename CoordsVec>
654 const CoordsVec& oop_atom_pos,
const ValueType& force_const)
656 ValueType chi_ijkl = calcOutOfPlaneAngle<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, oop_atom_pos) * ValueType(180 / M_PI);
657 ValueType e_oop = ValueType(0.5 * 0.043844) * force_const * chi_ijkl * chi_ijkl;
662 template <
typename ValueType,
typename CoordsVec>
664 const CoordsVec& oop_atom_pos,
const ValueType& r_jl,
const ValueType& force_const)
666 ValueType chi_ijkl = calcOutOfPlaneAngle<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, oop_atom_pos, r_jl) * ValueType(180 / M_PI);
667 ValueType e_oop = ValueType(0.5 * 0.043844) * force_const * chi_ijkl * chi_ijkl;
673 template <
typename ValueType,
typename Iter,
typename CoordsArray>
676 return Detail::accumInteractionEnergies<ValueType>(beg, end, coords,
677 static_cast<ValueType (*)(
const MMFF94TorsionInteraction&,
const CoordsArray&)
>(
678 &calcMMFF94TorsionEnergy<ValueType, CoordsArray>));
681 template <
typename ValueType,
typename CoordsArray>
684 return calcMMFF94TorsionEnergy<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtom1Index()],
685 coords[iaction.getCenterAtom2Index()], coords[iaction.getTerminalAtom2Index()],
686 iaction.getTorsionParameter1(), iaction.getTorsionParameter2(), iaction.getTorsionParameter3());
689 template <
typename ValueType,
typename CoordsVec>
691 const CoordsVec& term_atom2_pos,
const ValueType& tor_param1,
const ValueType& tor_param2,
692 const ValueType& tor_param3)
694 ValueType phi_cos = calcDihedralAngleCos<ValueType>(term_atom1_pos, ctr_atom1_pos, ctr_atom2_pos, term_atom2_pos);
695 ValueType phi = std::acos(phi_cos);
696 ValueType e_t = ValueType(0.5) * (tor_param1 * (1 + phi_cos) + tor_param2 * (1 - std::cos(2 * phi)) + tor_param3 * (1 + std::cos(3 * phi)));
702 template <
typename ValueType,
typename Iter,
typename CoordsArray>
705 return Detail::accumInteractionEnergies<ValueType>(beg, end, coords,
706 static_cast<ValueType (*)(
const MMFF94ElectrostaticInteraction&,
const CoordsArray&)
>(
707 &calcMMFF94ElectrostaticEnergy<ValueType, CoordsArray>));
710 template <
typename ValueType,
typename CoordsArray>
713 return calcMMFF94ElectrostaticEnergy<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()],
714 iaction.getAtom1Charge(), iaction.getAtom2Charge(), iaction.getScalingFactor(),
715 iaction.getDielectricConstant(), iaction.getDistanceExponent());
718 template <
typename ValueType,
typename CoordsVec>
720 const ValueType& atom2_chg,
const ValueType& scale_fact,
const ValueType& de_const,
721 const ValueType& dist_expo)
723 ValueType tmp = std::pow(calcDistance<ValueType>(atom1_pos, atom2_pos) + ValueType(0.05), dist_expo);
724 ValueType e_q = scale_fact * ValueType(332.0716) * atom1_chg * atom2_chg / (de_const * tmp);
729 template <
typename ValueType>
731 const ValueType& scale_fact,
const ValueType& de_const,
const ValueType& dist_expo)
733 ValueType tmp = std::pow(r_ij + ValueType(0.05), dist_expo);
734 ValueType e_q = scale_fact * ValueType(332.0716) * atom1_chg * atom2_chg / (de_const * tmp);
740 template <
typename ValueType,
typename Iter,
typename CoordsArray>
743 return Detail::accumInteractionEnergies<ValueType>(beg, end, coords,
744 static_cast<ValueType (*)(
const MMFF94VanDerWaalsInteraction&,
const CoordsArray&)
>(
745 &calcMMFF94VanDerWaalsEnergy<ValueType, CoordsArray>));
748 template <
typename ValueType,
typename CoordsArray>
751 return calcMMFF94VanDerWaalsEnergy<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()],
752 iaction.getEIJ(), iaction.getRIJ(), iaction.getRIJPow7());
755 template <
typename ValueType,
typename CoordsVec>
757 const ValueType& r_IJ,
const ValueType& r_IJ_7)
759 ValueType r_ij_2 = calcSquaredDistance<ValueType>(atom1_pos, atom2_pos);
760 ValueType r_ij = std::sqrt(r_ij_2);
761 ValueType r_ij_7 = r_ij_2 * r_ij_2 * r_ij_2 * r_ij;
763 ValueType tmp = ValueType(1.07) * r_IJ / (r_ij + ValueType(0.07) * r_IJ);
764 ValueType tmp_2 = tmp * tmp;
765 ValueType tmp_7 = tmp_2 * tmp_2 * tmp_2 * tmp;
767 ValueType e_vdw = e_IJ * tmp_7 * (ValueType(1.12) * r_IJ_7 / (r_ij_7 + ValueType(0.12) * r_IJ_7) - 2);
772 template <
typename ValueType>
774 const ValueType& r_IJ,
const ValueType& r_IJ_7)
776 ValueType r_ij_2 = r_ij * r_ij;
777 ValueType r_ij_7 = r_ij_2 * r_ij_2 * r_ij_2 * r_ij;
779 ValueType tmp = ValueType(1.07) * r_IJ / (r_ij + ValueType(0.07) * r_IJ);
780 ValueType tmp_2 = tmp * tmp;
781 ValueType tmp_7 = tmp_2 * tmp_2 * tmp_2 * tmp;
783 ValueType e_vdw = e_IJ * tmp_7 * (ValueType(1.12) * r_IJ_7 / (r_ij_7 + ValueType(0.12) * r_IJ_7) - 2);
790 #endif // CDPL_FORCEFIELD_MMFF94ENERGYFUNCTIONS_HPP