Precompute_Coulomb_Onsite.h

#ifndef QNANO_NEW_PRECOMPUTE_COULOMB_ONSITE
#define QNANO_NEW_PRECOMPUTE_COULOMB_ONSITE

#include "tools/Static_Const_Member.h"
#include "tools/Integrator.h"
#include "tools/QNANO_Constants.h"
#include "GSL/Integrator_VEGAS.h"
#include "orbitals/Slater_Orbital_List.h"

class Precompute_Coulomb_Onsite: public Integrable_Function{
public:

  Slater_Orbital_List solist;
  std::vector<int> cur_indices;



  static double f(double *x, size_t dim, void *params){
    const Precompute_Coulomb_Onsite *cl=(const Precompute_Coulomb_Onsite *)params;
    double *x2=&x[3];
    double DR=sqrt( (x[0]-x2[0])*(x[0]-x2[0]) + (x[1]-x2[1])*(x[1]-x2[1]) + (x[2]-x2[2])*(x[2]-x2[2]) );

    return QNANO_Constants::COULOMB_PREFAC* \
       cl->solist[cl->cur_indices[0]].f(x,3,(void*)&cl->solist[cl->cur_indices[0]])* \
       cl->solist[cl->cur_indices[1]].f(x2,3,(void*)&cl->solist[cl->cur_indices[1]])* \
       cl->solist[cl->cur_indices[2]].f(x2,3,(void*)&cl->solist[cl->cur_indices[2]])* \
       cl->solist[cl->cur_indices[3]].f(x,3,(void*)&cl->solist[cl->cur_indices[3]])/DR;
  }


  void precompute_and_print(const std::string &outfile, double lim, double Ndiscr, double epsilon=1.);


  Precompute_Coulomb_Onsite(const std::string &orbfile) 
   : Integrable_Function(f, 6, this), solist(orbfile), cur_indices(4,0) {
  }

};

#endif