29 #ifndef CDPL_FORCEFIELD_MMFF94GRADIENTFUNCTIONS_HPP
30 #define CDPL_FORCEFIELD_MMFF94GRADIENTFUNCTIONS_HPP
50 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
53 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
94 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
96 const ValueType& force_const,
const ValueType& ref_length);
99 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
102 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
154 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
156 GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad,
bool linear,
157 const ValueType& force_const,
const ValueType& ref_angle);
160 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
163 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
215 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
218 GradVec& ctr_atom_grad, GradVec& term_atom2_grad,
const ValueType& ijk_force_const,
const ValueType& kji_force_const,
219 const ValueType& ref_angle,
const ValueType& ref_length1,
const ValueType& ref_length2);
222 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
225 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
265 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
267 const CoordsVec& oop_atom_pos, GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad,
268 GradVec& oop_atom_grad,
const ValueType& force_const);
271 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
274 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
316 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
318 const CoordsVec& term_atom2_pos, GradVec& term_atom1_grad, GradVec& ctr_atom1_grad, GradVec& ctr_atom2_grad,
319 GradVec& term_atom2_grad,
const ValueType& tor_param1,
const ValueType& tor_param2,
const ValueType& tor_param3);
322 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
325 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
369 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
371 const ValueType& atom1_chg,
const ValueType& atom2_chg,
const ValueType& scale_fact,
372 const ValueType& de_const,
const ValueType& dist_expo);
375 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
378 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
436 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
438 const ValueType& e_IJ,
const ValueType& r_IJ,
const ValueType& r_IJ_7);
446 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
449 return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
450 static_cast<ValueType (*)(
const MMFF94BondStretchingInteraction&,
const CoordsArray&, GradVector&)
>(
451 &calcMMFF94BondStretchingGradient<ValueType, CoordsArray, GradVector>));
454 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
457 return calcMMFF94BondStretchingGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
458 grad[iaction.getAtom2Index()], iaction.getForceConstant(), iaction.getReferenceLength());
461 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
463 const ValueType& force_const,
const ValueType& ref_length)
465 ValueType dist_atom1_grad[3];
466 ValueType dist_atom2_grad[3];
468 ValueType dr_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad) - ref_length;
469 ValueType dr_ij_2 = dr_ij * dr_ij;
471 ValueType grad_fact = (ValueType(167.92125 * 4) * dr_ij_2 * dr_ij - ValueType(215.89875 * 2) * dr_ij_2 + ValueType(143.9325) * dr_ij) * force_const;
473 Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
474 Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
476 ValueType e_b = ValueType(143.9325 * 0.5) * force_const * dr_ij_2 * (1 - 2 * dr_ij + 28 * dr_ij_2 / 12);
482 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
485 return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
486 static_cast<ValueType (*)(
const MMFF94AngleBendingInteraction&,
const CoordsArray&, GradVector&)
>(
487 &calcMMFF94AngleBendingGradient<ValueType, CoordsArray, GradVector>));
490 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
493 return calcMMFF94AngleBendingGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
494 coords[iaction.getTerminalAtom2Index()], grad[iaction.getTerminalAtom1Index()],
495 grad[iaction.getCenterAtomIndex()], grad[iaction.getTerminalAtom2Index()], iaction.isLinearAngle(),
496 iaction.getForceConstant(), iaction.getReferenceAngle());
499 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
501 GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad,
bool linear,
502 const ValueType& force_const,
const ValueType& ref_angle)
504 ValueType ac_term1_grad[3];
505 ValueType ac_ctr_grad[3];
506 ValueType ac_term2_grad[3];
508 ValueType a_ijk_cos = calcBondAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, ac_term1_grad, ac_ctr_grad, ac_term2_grad);
509 ValueType grad_fact = ValueType(1);
510 ValueType e_a = ValueType(0);
513 grad_fact = ValueType(143.9325) * force_const;
514 e_a = ValueType(143.9325) * force_const * (1 + a_ijk_cos);
517 ValueType a_ijk_cos_2 = a_ijk_cos * a_ijk_cos;
518 ValueType a_ijk = std::acos(a_ijk_cos);
519 ValueType div = std::sqrt(1 - a_ijk_cos_2);
521 if (div < ValueType(0.0000001))
522 div = ValueType(0.0000001);
524 grad_fact = force_const / div *
525 (a_ijk * (ValueType(86.58992538) * a_ijk - ValueType(143.9313616)) -
526 ref_angle * (ValueType(3.022558594) * a_ijk - ValueType(0.02637679965) * ref_angle - ValueType(2.512076157)));
528 ValueType da_ijk = a_ijk * ValueType(180 / M_PI) - ref_angle;
530 e_a = ValueType(0.043844 * 0.5) * force_const * da_ijk * da_ijk * (1 - ValueType(0.007) * da_ijk);
533 Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
534 Detail::scaleAddVector(ac_ctr_grad, grad_fact, ctr_atom_grad);
535 Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
541 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
544 return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
545 static_cast<ValueType (*)(
const MMFF94StretchBendInteraction&,
const CoordsArray&, GradVector&)
>(
546 &calcMMFF94StretchBendGradient<ValueType, CoordsArray, GradVector>));
549 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
552 return calcMMFF94StretchBendGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
553 coords[iaction.getTerminalAtom2Index()], grad[iaction.getTerminalAtom1Index()],
554 grad[iaction.getCenterAtomIndex()], grad[iaction.getTerminalAtom2Index()], iaction.getIJKForceConstant(),
555 iaction.getKJIForceConstant(), iaction.getReferenceAngle(), iaction.getReferenceLength1(),
556 iaction.getReferenceLength2());
559 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
561 GradVec& term_atom1_grad, GradVec& ctr_atom_grad, GradVec& term_atom2_grad,
562 const ValueType& ijk_force_const,
const ValueType& kji_force_const,
const ValueType& ref_angle,
563 const ValueType& ref_length1,
const ValueType& ref_length2)
565 ValueType ac_term1_grad[3];
566 ValueType ac_ctr_grad[3];
567 ValueType ac_term2_grad[3];
569 ValueType dist_term1_grad[3];
570 ValueType dist_ctr_grad1[3];
571 ValueType dist_ctr_grad2[3];
572 ValueType dist_term2_grad[3];
574 ValueType r_ij = calcDistanceDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, dist_term1_grad, dist_ctr_grad1);
575 ValueType r_kj = calcDistanceDerivatives<ValueType>(term_atom2_pos, ctr_atom_pos, dist_term2_grad, dist_ctr_grad2);
576 ValueType a_ijk_cos = calcBondAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, ac_term1_grad, ac_ctr_grad, ac_term2_grad);
577 ValueType a_ijk_cos_2 = a_ijk_cos * a_ijk_cos;
578 ValueType a_ijk = std::acos(a_ijk_cos);
580 ValueType dr_ij = r_ij - ref_length1;
581 ValueType dr_kj = r_kj - ref_length2;
582 ValueType da_ijk = a_ijk * ValueType(180 / M_PI) - ref_angle;
583 ValueType div = std::sqrt(1 - a_ijk_cos_2);
585 if (div < ValueType(0.0000001))
586 div = ValueType(0.0000001);
588 ValueType a_ijk_grad_fact = ValueType(-180 * 2.5121 / M_PI) / div * (dr_ij * ijk_force_const + dr_kj * kji_force_const);
590 ValueType r_ij_grad_fact = ValueType(2.5121) * da_ijk * ijk_force_const;
591 ValueType r_kj_grad_fact = ValueType(2.5121) * da_ijk * kji_force_const;
593 Detail::scaleAddVector(dist_term1_grad, r_ij_grad_fact, term_atom1_grad);
594 Detail::scaleAddVector(ac_term1_grad, a_ijk_grad_fact, term_atom1_grad);
596 Detail::scaleAddVector(dist_term2_grad, r_kj_grad_fact, term_atom2_grad);
597 Detail::scaleAddVector(ac_term2_grad, a_ijk_grad_fact, term_atom2_grad);
599 Detail::scaleAddVector(dist_ctr_grad1, r_ij_grad_fact, ctr_atom_grad);
600 Detail::scaleAddVector(dist_ctr_grad2, r_kj_grad_fact, ctr_atom_grad);
601 Detail::scaleAddVector(ac_ctr_grad, a_ijk_grad_fact, ctr_atom_grad);
603 ValueType e_ab = r_ij_grad_fact * dr_ij + r_kj_grad_fact * dr_kj;
609 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
612 return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
613 static_cast<ValueType (*)(
const MMFF94OutOfPlaneBendingInteraction&,
const CoordsArray&, GradVector&)
>(
614 &calcMMFF94OutOfPlaneBendingGradient<ValueType, CoordsArray, GradVector>));
617 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
620 return calcMMFF94OutOfPlaneBendingGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtomIndex()],
621 coords[iaction.getTerminalAtom2Index()], coords[iaction.getOutOfPlaneAtomIndex()],
622 grad[iaction.getTerminalAtom1Index()], grad[iaction.getCenterAtomIndex()],
623 grad[iaction.getTerminalAtom2Index()], grad[iaction.getOutOfPlaneAtomIndex()],
624 iaction.getForceConstant());
627 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
629 const CoordsVec& oop_atom_pos, GradVec& term_atom1_grad, GradVec& ctr_atom_grad,
630 GradVec& term_atom2_grad, GradVec& oop_atom_grad,
const ValueType& force_const)
632 ValueType ac_term1_grad[3];
633 ValueType ac_term2_grad[3];
634 ValueType ac_ctr_grad[3];
635 ValueType ac_oop_grad[3];
637 ValueType chi_ijkl_cos = calcOutOfPlaneAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom_pos, term_atom2_pos, oop_atom_pos, ac_term1_grad,
638 ac_ctr_grad, ac_term2_grad, ac_oop_grad);
639 ValueType chi_ijkl = ValueType(M_PI * 0.5) - std::acos(chi_ijkl_cos);
640 ValueType div = std::sqrt(1 - chi_ijkl_cos * chi_ijkl_cos);
642 if (div < ValueType(0.0000001))
643 div = ValueType(0.0000001);
645 ValueType grad_fact = ValueType(0.043844 * 180 * 180) / div * ValueType(1 / (M_PI * M_PI)) * force_const * chi_ijkl;
647 Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
648 Detail::scaleAddVector(ac_ctr_grad, grad_fact, ctr_atom_grad);
649 Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
650 Detail::scaleAddVector(ac_oop_grad, grad_fact, oop_atom_grad);
652 chi_ijkl *= ValueType(180 / M_PI);
654 ValueType e_oop = ValueType(0.5 * 0.043844) * force_const * chi_ijkl * chi_ijkl;
660 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
663 return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
664 static_cast<ValueType (*)(
const MMFF94TorsionInteraction&,
const CoordsArray&, GradVector&)
>(
665 &calcMMFF94TorsionGradient<ValueType, CoordsArray, GradVector>));
668 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
671 return calcMMFF94TorsionGradient<ValueType>(coords[iaction.getTerminalAtom1Index()], coords[iaction.getCenterAtom1Index()],
672 coords[iaction.getCenterAtom2Index()], coords[iaction.getTerminalAtom2Index()],
673 grad[iaction.getTerminalAtom1Index()], grad[iaction.getCenterAtom1Index()], grad[iaction.getCenterAtom2Index()],
674 grad[iaction.getTerminalAtom2Index()], iaction.getTorsionParameter1(), iaction.getTorsionParameter2(),
675 iaction.getTorsionParameter3());
678 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
681 const CoordsVec& term_atom2_pos, GradVec& term_atom1_grad, GradVec& ctr_atom1_grad, GradVec& ctr_atom2_grad,
682 GradVec& term_atom2_grad,
const ValueType& tor_param1,
const ValueType& tor_param2,
const ValueType& tor_param3)
684 ValueType ac_term1_grad[3];
685 ValueType ac_ctr1_grad[3];
686 ValueType ac_ctr2_grad[3];
687 ValueType ac_term2_grad[3];
689 ValueType phi_cos = calcDihedralAngleCosDerivatives<ValueType>(term_atom1_pos, ctr_atom1_pos, ctr_atom2_pos, term_atom2_pos, ac_term1_grad, ac_ctr1_grad,
690 ac_ctr2_grad, ac_term2_grad);
691 ValueType phi = std::acos(phi_cos);
692 ValueType phi_cos_2 = phi_cos * phi_cos;
693 ValueType div = std::sqrt(1 - phi_cos_2);
695 if (div < ValueType(0.0000001))
696 div = ValueType(0.0000001);
698 ValueType grad_fact =
699 ValueType(-1) / div * (tor_param2 * std::sin(2 * phi) - ValueType(0.5) * tor_param1 * std::sin(phi) - ValueType(1.5) * tor_param3 * std::sin(3 * phi));
701 Detail::scaleAddVector(ac_term1_grad, grad_fact, term_atom1_grad);
702 Detail::scaleAddVector(ac_ctr1_grad, grad_fact, ctr_atom1_grad);
703 Detail::scaleAddVector(ac_ctr2_grad, grad_fact, ctr_atom2_grad);
704 Detail::scaleAddVector(ac_term2_grad, grad_fact, term_atom2_grad);
706 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)));
712 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
715 return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
716 static_cast<ValueType (*)(
const MMFF94ElectrostaticInteraction&,
const CoordsArray&, GradVector&)
>(
717 &calcMMFF94ElectrostaticGradient<ValueType, CoordsArray, GradVector>));
720 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
723 return calcMMFF94ElectrostaticGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
724 grad[iaction.getAtom2Index()], iaction.getAtom1Charge(), iaction.getAtom2Charge(),
725 iaction.getScalingFactor(), iaction.getDielectricConstant(), iaction.getDistanceExponent());
728 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
730 const ValueType& atom1_chg,
const ValueType& atom2_chg,
const ValueType& scale_fact,
731 const ValueType& de_const,
const ValueType& dist_expo)
733 ValueType dist_atom1_grad[3];
734 ValueType dist_atom2_grad[3];
736 ValueType r_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad);
738 ValueType tmp1 = r_ij + ValueType(0.05);
739 ValueType tmp2 = std::pow(tmp1, dist_expo);
740 ValueType tmp3 = scale_fact * atom1_chg * atom2_chg / (de_const * tmp2);
742 ValueType grad_fact = ValueType(-332.0716) * dist_expo * tmp3 / tmp1;
744 Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
745 Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
747 double e_q = ValueType(332.0716) * tmp3;
753 template <
typename ValueType,
typename Iter,
typename CoordsArray,
typename GradVector>
756 return Detail::calcInteractionGradient<ValueType>(beg, end, coords, grad,
757 static_cast<ValueType (*)(
const MMFF94VanDerWaalsInteraction&,
const CoordsArray&, GradVector&)
>(
758 &calcMMFF94VanDerWaalsGradient<ValueType, CoordsArray, GradVector>));
761 template <
typename ValueType,
typename CoordsArray,
typename GradVector>
764 return calcMMFF94VanDerWaalsGradient<ValueType>(coords[iaction.getAtom1Index()], coords[iaction.getAtom2Index()], grad[iaction.getAtom1Index()],
765 grad[iaction.getAtom2Index()], iaction.getEIJ(), iaction.getRIJ(), iaction.getRIJPow7());
768 template <
typename ValueType,
typename CoordsVec,
typename GradVec>
770 const ValueType& e_IJ,
const ValueType& r_IJ,
const ValueType& r_IJ_7)
772 ValueType dist_atom1_grad[3];
773 ValueType dist_atom2_grad[3];
775 ValueType r_ij = calcDistanceDerivatives<ValueType>(atom1_pos, atom2_pos, dist_atom1_grad, dist_atom2_grad);
776 ValueType r_ij_2 = r_ij * r_ij;
777 ValueType r_ij_6 = r_ij_2 * r_ij_2 * r_ij_2;
778 ValueType r_ij_7 = r_ij_6 * r_ij;
780 ValueType tmp1 = r_ij + ValueType(0.07) * r_IJ;
781 ValueType tmp1_2 = tmp1 * tmp1;
782 ValueType tmp1_4 = tmp1_2 * tmp1_2;
784 ValueType tmp2 = r_ij_7 + ValueType(0.12) * r_IJ_7;
786 ValueType tmp3 = ValueType(1.07) * r_IJ / (r_ij + ValueType(0.07) * r_IJ);
787 ValueType tmp3_2 = tmp3 * tmp3;
788 ValueType tmp3_7 = tmp3_2 * tmp3_2 * tmp3_2 * tmp3;
790 ValueType grad_fact = -r_IJ_7 * e_IJ / (tmp1_4 * tmp1_4 * tmp2 * tmp2) *
791 (ValueType(-22.48094067) * r_ij_7 * r_ij_7 + ValueType(19.78322779) * r_ij_7 * r_IJ_7 +
792 ValueType(0.8812528743) * r_ij_6 * r_IJ_7 * r_IJ + ValueType(1.186993667) * r_IJ_7 * r_IJ_7);
794 Detail::scaleAddVector(dist_atom1_grad, grad_fact, atom1_grad);
795 Detail::scaleAddVector(dist_atom2_grad, grad_fact, atom2_grad);
797 ValueType e_vdw = e_IJ * tmp3_7 * (ValueType(1.12) * r_IJ_7 / tmp2 - 2);
804 #endif // CDPL_FORCEFIELD_MMFF94GRADIENTFUNCTIONS_HPP