<commit_before>// $Id$
//
//  Copyright (C) 2004-2006 Rational Discovery LLC
//
//   @@ All Rights Reserved  @@
//
#include "AngleBend.h"
#include "BondStretch.h"
#include "Params.h"
#include <math.h>
#include <ForceField/ForceField.h>
#include <RDGeneral/Invariant.h>

namespace ForceFields {
  namespace UFF {

    namespace Utils {
      double calcAngleForceConstant(double theta0,
				    double bondOrder12,double bondOrder23,
				    const AtomicParams *at1Params,
				    const AtomicParams *at2Params,
				    const AtomicParams *at3Params){
	double cosTheta0=cos(theta0);
	double r12 = calcBondRestLength(bondOrder12,at1Params,at2Params);
	double r23 = calcBondRestLength(bondOrder23,at2Params,at3Params);
	double r13 = sqrt(r12*r12 + r23*r23 - 2.*r12*r23*cosTheta0);
	double beta = 2.*Params::G/(r12*r23);

	double preFactor = beta*at1Params->Z1*at3Params->Z1 / pow(r13,5);
	double rTerm = r12*r23;
	double innerBit = r12*r23*(1-cosTheta0*cosTheta0) - r13*r13*cosTheta0;
            
	double res=preFactor*rTerm*innerBit;
	return res;
      }
    } // end of namespace Utils
  
    AngleBendContrib::AngleBendContrib(ForceField *owner,
				       unsigned int idx1,unsigned int idx2,unsigned int idx3,
				       double bondOrder12,double bondOrder23,
				       const AtomicParams *at1Params,
				       const AtomicParams *at2Params,
				       const AtomicParams *at3Params,
				       unsigned int order){
      PRECONDITION(owner,"bad owner");
      PRECONDITION(at1Params,"bad params pointer");
      PRECONDITION(at2Params,"bad params pointer");
      PRECONDITION(at3Params,"bad params pointer");
      PRECONDITION((idx1!=idx2&&idx2!=idx3&&idx1!=idx3),"degenerate points");
      RANGE_CHECK(0,idx1,owner->positions().size()-1);
      RANGE_CHECK(0,idx2,owner->positions().size()-1);
      RANGE_CHECK(0,idx3,owner->positions().size()-1);
      dp_forceField = owner;
      d_at1Idx = idx1;
      d_at2Idx = idx2;
      d_at3Idx = idx3;
      d_order = order;
      this->d_forceConstant = Utils::calcAngleForceConstant(at2Params->theta0,
							    bondOrder12,bondOrder23,
							    at1Params,at2Params,at3Params);
      if(order==0){
	double sinTheta0=sin(at2Params->theta0);
	double cosTheta0=cos(at2Params->theta0);
	this->d_C2 = 1./(4.*sinTheta0*sinTheta0);
	this->d_C1 = -4.*this->d_C2*cosTheta0;
	this->d_C0 = this->d_C2*(2.*cosTheta0*cosTheta0 + 1.);
      }
    }

    double AngleBendContrib::getEnergy(double *pos) const {
      PRECONDITION(dp_forceField,"no owner");
      PRECONDITION(pos,"bad vector");

      double dist1=this->dp_forceField->distance(this->d_at1Idx,this->d_at2Idx,pos);
      double dist2=this->dp_forceField->distance(this->d_at2Idx,this->d_at3Idx,pos);

      RDGeom::Point3D p1(pos[3*this->d_at1Idx],
			 pos[3*this->d_at1Idx+1],
			 pos[3*this->d_at1Idx+2]);
      RDGeom::Point3D p2(pos[3*this->d_at2Idx],
			 pos[3*this->d_at2Idx+1],
			 pos[3*this->d_at2Idx+2]);
      RDGeom::Point3D p3(pos[3*this->d_at3Idx],
			 pos[3*this->d_at3Idx+1],
			 pos[3*this->d_at3Idx+2]);
      RDGeom::Point3D p12=p1-p2;
      RDGeom::Point3D p32=p3-p2;
      double cosTheta = p12.dotProduct(p32)/(dist1*dist2);
      // we need sin^2(theta) to get cos(2*theta), so compute that:
      double sinThetaSq = 1-cosTheta*cosTheta;
    
      double angleTerm = this->getEnergyTerm(cosTheta,sinThetaSq);
      double res = this->d_forceConstant*angleTerm;

      return res;
    }

    void AngleBendContrib::getGrad(double *pos,double *grad) const {
      PRECONDITION(dp_forceField,"no owner");
      PRECONDITION(pos,"bad vector");
      PRECONDITION(grad,"bad vector");

      double dist1=this->dp_forceField->distance(this->d_at1Idx,this->d_at2Idx,pos);
      double dist2=this->dp_forceField->distance(this->d_at2Idx,this->d_at3Idx,pos);

      //std::cout << "\tAngle("<<this->d_at1Idx<<","<<this->d_at2Idx<<","<<this->d_at3Idx<<") " << dist1 << " " << dist2 << std::endl;
      
      RDGeom::Point3D p1(pos[3*this->d_at1Idx],
			 pos[3*this->d_at1Idx+1],
			 pos[3*this->d_at1Idx+2]);
      RDGeom::Point3D p2(pos[3*this->d_at2Idx],
			 pos[3*this->d_at2Idx+1],
			 pos[3*this->d_at2Idx+2]);
      RDGeom::Point3D p3(pos[3*this->d_at3Idx],
			 pos[3*this->d_at3Idx+1],
			 pos[3*this->d_at3Idx+2]);
      double *g1=&(grad[3*this->d_at1Idx]);
      double *g2=&(grad[3*this->d_at2Idx]);
      double *g3=&(grad[3*this->d_at3Idx]);

      RDGeom::Point3D p12=p1-p2;
      RDGeom::Point3D p32=p3-p2;
      double cosTheta = p12.dotProduct(p32)/(dist1*dist2);
      double sinTheta = sqrt(1.0-cosTheta*cosTheta);


      // ------------------------------------
      // FIX: handle cases where sinTheta is zero
      // ------------------------------------


      //std::cerr << "GRAD: " << cosTheta << " (" << acos(cosTheta)<< "), ";
      //std::cerr << sinTheta << " (" << asin(sinTheta)<< ")" << std::endl;
    
      // use the chain rule:
      // dE/dx = dE/dTheta * dTheta/dx

      // dE/dTheta is independent of cartesians:
      double dE_dTheta=getThetaDeriv(cosTheta,sinTheta);
    <commit_msg>catch a possible divide by zero<commit_after>// $Id$
//
//  Copyright (C) 2004-2006 Rational Discovery LLC
//
//   @@ All Rights Reserved  @@
//
#include "AngleBend.h"
#include "BondStretch.h"
#include "Params.h"
#include <math.h>
#include <ForceField/ForceField.h>
#include <RDGeneral/Invariant.h>

namespace ForceFields {
  namespace UFF {

    namespace Utils {
      double calcAngleForceConstant(double theta0,
				    double bondOrder12,double bondOrder23,
				    const AtomicParams *at1Params,
				    const AtomicParams *at2Params,
				    const AtomicParams *at3Params){
	double cosTheta0=cos(theta0);
	double r12 = calcBondRestLength(bondOrder12,at1Params,at2Params);
	double r23 = calcBondRestLength(bondOrder23,at2Params,at3Params);
	double r13 = sqrt(r12*r12 + r23*r23 - 2.*r12*r23*cosTheta0);
	double beta = 2.*Params::G/(r12*r23);

	double preFactor = beta*at1Params->Z1*at3Params->Z1 / pow(r13,5);
	double rTerm = r12*r23;
	double innerBit = r12*r23*(1-cosTheta0*cosTheta0) - r13*r13*cosTheta0;
            
	double res=preFactor*rTerm*innerBit;
	return res;
      }
    } // end of namespace Utils
  
    AngleBendContrib::AngleBendContrib(ForceField *owner,
				       unsigned int idx1,unsigned int idx2,unsigned int idx3,
				       double bondOrder12,double bondOrder23,
				       const AtomicParams *at1Params,
				       const AtomicParams *at2Params,
				       const AtomicParams *at3Params,
				       unsigned int order){
      PRECONDITION(owner,"bad owner");
      PRECONDITION(at1Params,"bad params pointer");
      PRECONDITION(at2Params,"bad params pointer");
      PRECONDITION(at3Params,"bad params pointer");
      PRECONDITION((idx1!=idx2&&idx2!=idx3&&idx1!=idx3),"degenerate points");
      RANGE_CHECK(0,idx1,owner->positions().size()-1);
      RANGE_CHECK(0,idx2,owner->positions().size()-1);
      RANGE_CHECK(0,idx3,owner->positions().size()-1);
      dp_forceField = owner;
      d_at1Idx = idx1;
      d_at2Idx = idx2;
      d_at3Idx = idx3;
      d_order = order;
      this->d_forceConstant = Utils::calcAngleForceConstant(at2Params->theta0,
							    bondOrder12,bondOrder23,
							    at1Params,at2Params,at3Params);
      if(order==0){
	double sinTheta0=sin(at2Params->theta0);
	double cosTheta0=cos(at2Params->theta0);
	this->d_C2 = 1./(4.*std::max(sinTheta0*sinTheta0,1e-8));
	this->d_C1 = -4.*this->d_C2*cosTheta0;
	this->d_C0 = this->d_C2*(2.*cosTheta0*cosTheta0 + 1.);
      }
    }

    double AngleBendContrib::getEnergy(double *pos) const {
      PRECONDITION(dp_forceField,"no owner");
      PRECONDITION(pos,"bad vector");

      double dist1=this->dp_forceField->distance(this->d_at1Idx,this->d_at2Idx,pos);
      double dist2=this->dp_forceField->distance(this->d_at2Idx,this->d_at3Idx,pos);

      RDGeom::Point3D p1(pos[3*this->d_at1Idx],
			 pos[3*this->d_at1Idx+1],
			 pos[3*this->d_at1Idx+2]);
      RDGeom::Point3D p2(pos[3*this->d_at2Idx],
			 pos[3*this->d_at2Idx+1],
			 pos[3*this->d_at2Idx+2]);
      RDGeom::Point3D p3(pos[3*this->d_at3Idx],
			 pos[3*this->d_at3Idx+1],
			 pos[3*this->d_at3Idx+2]);
      RDGeom::Point3D p12=p1-p2;
      RDGeom::Point3D p32=p3-p2;
      double cosTheta = p12.dotProduct(p32)/(dist1*dist2);
      // we need sin^2(theta) to get cos(2*theta), so compute that:
      double sinThetaSq = 1-cosTheta*cosTheta;
    
      double angleTerm = this->getEnergyTerm(cosTheta,sinThetaSq);
      double res = this->d_forceConstant*angleTerm;

      return res;
    }

    void AngleBendContrib::getGrad(double *pos,double *grad) const {
      PRECONDITION(dp_forceField,"no owner");
      PRECONDITION(pos,"bad vector");
      PRECONDITION(grad,"bad vector");

      double dist1=this->dp_forceField->distance(this->d_at1Idx,this->d_at2Idx,pos);
      double dist2=this->dp_forceField->distance(this->d_at2Idx,this->d_at3Idx,pos);

      //std::cout << "\tAngle("<<this->d_at1Idx<<","<<this->d_at2Idx<<","<<this->d_at3Idx<<") " << dist1 << " " << dist2 << std::endl;
      
      RDGeom::Point3D p1(pos[3*this->d_at1Idx],
			 pos[3*this->d_at1Idx+1],
			 pos[3*this->d_at1Idx+2]);
      RDGeom::Point3D p2(pos[3*this->d_at2Idx],
			 pos[3*this->d_at2Idx+1],
			 pos[3*this->d_at2Idx+2]);
      RDGeom::Point3D p3(pos[3*this->d_at3Idx],
			 pos[3*this->d_at3Idx+1],
			 pos[3*this->d_at3Idx+2]);
      double *g1=&(grad[3*this->d_at1Idx]);
      double *g2=&(grad[3*this->d_at2Idx]);
      double *g3=&(grad[3*this->d_at3Idx]);

      RDGeom::Point3D p12=p1-p2;
      RDGeom::Point3D p32=p3-p2;
      double cosTheta = p12.dotProduct(p32)/(dist1*dist2);
      double sinTheta = std::max(sqrt(1.0-cosTheta*cosTheta),1e-8);

      //std::cerr << "GRAD: " << cosTheta << " (" << acos(cosTheta)<< "), ";
      //std::cerr << sinTheta << " (" << asin(sinTheta)<< ")" << std::endl;
    
      // use the chain rule:
      // dE/dx = dE/dTheta * dTheta/dx

      // dE/dTheta is independent of cartesians:
      double dE_dTheta=getThetaDeriv(cosTheta,sinTheta);
    
{'fraction_non_alphanumeric': 0.16605092361457813, 'fraction_numerical': 0.0527209186220669, 'mean_word_length': 3.4456280514869064, 'pattern_counts': {'":': 0, '<': 93, '<?xml version=': 0, '>': 121, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1058374', 'n_tokens_mistral': 4152, 'n_tokens_neox': 3749, 'n_words': 544}