soecp_eval_reference.cpp File Reference

Include dependency graph for soecp_eval_reference.cpp:

Go to the source code of this file.

Classes
class	cubicBSpline

Functions
double	B (double x, int k, int i, const std::vector< double > &t)
double	bspline (double x, const std::vector< double > &t, const std::vector< double > &c, int k)
void	testJastrow ()
std::complex< double >	Ylm (int l, int m, const std::vector< double > &sph)
std::vector< double >	cart2sph (const std::vector< double > &cart)
std::vector< double >	sph2cart (const std::vector< double > &sph)
std::complex< double >	chiu (double s)
std::complex< double >	chid (double s)
std::complex< double >	orbu (const std::vector< double > &sph)
std::complex< double >	orbd (const std::vector< double > &sph)
std::complex< double >	spinor0 (const std::vector< double > &sph, double s)
std::complex< double >	spinor1 (const std::vector< double > &sph, double s)
std::complex< double >	sMatrixElement (double s1, double s2, int d)
int	kroneckerDelta (int x, int y)
std::complex< double >	lMatrixElement (int l, int m1, int m2, int d)
double	Wso (int l, double r)
std::complex< double >	TWF (std::vector< std::vector< double >> &positions, std::vector< double > &spins, const cubicBSpline &J2)
std::complex< double >	dratioTWF (const int ip, std::vector< std::vector< double >> &positions, std::vector< double > &spins, cubicBSpline &J2)
std::complex< double >	calcAngInt (int iel, double s2, std::vector< std::vector< double >> &positions, std::vector< double > &spins, const cubicBSpline &J2)
std::complex< double >	calcVal (int npts, int iel, std::vector< std::vector< double >> &positions, std::vector< double > &spins, const cubicBSpline &J2)
void	calcSOECP (int npts)
void	evaluateValueAndDerivative ()
int	main ()

Variables
double	PI = 4. * std::atan(1)
std::vector< std::vector< double > >	quad
std::vector< double >	wt
std::complex< double >	I = std::complex<double>(0, 1)

Function Documentation

B()

double B	(	double	x,
		int	k,
		int	i,
		const std::vector< double > &	t
	)

Definition at line 27 of file soecp_eval_reference.cpp.

Referenced by MultiQuinticSpline1D< T >::add_spline(), qmcplusplus::bessel_steed_array_cpu(), bspline(), TWFFastDerivWrapper::buildX(), LatticeAnalyzer< T, 3 >::calcSimulationCellRadius(), CubicFormula(), qmcplusplus::MatrixOperators::diag_product(), RotatedSPOs::evaluateDerivatives(), LRHandlerBase::evaluateGrad(), NonLocalECPotential::evaluateOneBodyOpMatrix(), BareKineticEnergy::evaluateOneBodyOpMatrix(), NonLocalECPComponent::evaluateOneBodyOpMatrixContribution(), qmcplusplus::cuBLAS::geam(), qmcplusplus::compute::BLAS::gemm(), BLAS::gemm(), qmcplusplus::syclBLAS::gemm(), qmcplusplus::ompBLAS::gemm< double >(), qmcplusplus::ompBLAS::gemm< float >(), qmcplusplus::ompBLAS::gemm< std::complex< double > >(), qmcplusplus::ompBLAS::gemm< std::complex< float > >(), qmcplusplus::compute::BLAS::gemm_batched(), qmcplusplus::ompBLAS::gemm_batched< double >(), qmcplusplus::ompBLAS::gemm_batched< float >(), qmcplusplus::ompBLAS::gemm_batched< std::complex< double > >(), qmcplusplus::ompBLAS::gemm_batched< std::complex< float > >(), qmcplusplus::ompBLAS::gemm_batched_impl(), qmcplusplus::ompBLAS::gemm_impl(), LinearMethod::getLowestEigenvector_Gen(), LinearMethod::getLowestEigenvector_Inv(), qmcplusplus::cusolver::getrs(), qmcplusplus::rocsolver::getrs(), qmcplusplus::MatrixOperators::half_outerProduct(), qmcplusplus::MatrixOperators::innerProduct(), TestFillBufferRngReal< T >::operator()(), TestFillVecRngReal< T >::operator()(), TestGetVecRngReal< T >::operator()(), TestFillBufferRngComplex< T >::operator()(), TestFillVecRngComplex< T >::operator()(), TestGetVecRngComplex< T >::operator()(), operator*(), operator+(), operator+=(), operator-(), operator-=(), SymmArray< T >::operator=(), qmcplusplus::MatrixOperators::other_half_outerProduct(), qmcplusplus::MatrixOperators::product(), qmcplusplus::MatrixOperators::product_ABt(), qmcplusplus::Product_ABt(), qmcplusplus::MatrixOperators::product_AtB(), RadialWF::PseudoDerivs(), QuinticSplineSolve(), qmcplusplus::simd::remapCopy(), SimCellRad(), Eigensolver::solveGeneralizedEigenvalues(), Eigensolver::solveGeneralizedEigenvalues_Inv(), SymmArray< T >::SymmArray(), RotatedSPOs::table_method_eval(), qmcplusplus::TEST_CASE(), qmcplusplus::test_gemm(), qmcplusplus::test_gemv(), qmcplusplus::test_one_gemm(), qmcplusplus::test_one_gemv(), TWFFastDerivWrapper::trAB(), qmcplusplus::simd::transpose(), qmcplusplus::MatrixOperators::transpose(), and OrbitalImages::write_orbital_xsf().

{
  double c1, c2;
  if (k == 0)
    return (x >= t[i] && x < t[i + 1]) ? 1.0 : 0.0;
  if (t[i + k] == t[i])
    c1 = 0.0;
  else
    c1 = (x - t[i]) / (t[i + k] - t[i]) * B(x, k - 1, i, t);
  if (t[i + k + 1] == t[i + 1])
    c2 = 0.0;
  else
    c2 = (t[i + k + 1] - x) / (t[i + k + 1] - t[i + 1]) * B(x, k - 1, i + 1, t);
  return c1 + c2;
}

bspline()

double bspline	(	double	x,
		const std::vector< double > &	t,
		const std::vector< double > &	c,
		int	k
	)

Definition at line 43 of file soecp_eval_reference.cpp.

References B(), and qmcplusplus::n.

Referenced by BsplineReader::check_twists(), HybridRepSetReader< SA >::create_atomic_centers_Gspace(), HybridRepSetReader< SA >::create_spline_set(), SplineSetReader< typename SA::SplineBase >::create_spline_set(), SplineSetReader< typename SA::SplineBase >::createSplineDataSpaceLookforDumpFile(), cubicBSpline::getVal(), HybridRepSetReader< SA >::initialize_hybrid_pio_gather(), HybridRepSetReader< SA >::initialize_hybridrep_atomic_centers(), and SplineSetReader< typename SA::SplineBase >::initialize_spline_pio_gather().

{
  int n = t.size() - k - 1;
  assert(n >= k + 1);
  assert(n <= c.size());
  double val = 0.0;
  for (int i = 0; i < n; i++)
    val += c[i] * B(x, k, i, t);
  return val;
}

calcAngInt()

std::complex<double> calcAngInt	(	int	iel,
		double	s2,
		std::vector< std::vector< double >> &	positions,
		std::vector< double > &	spins,
		const cubicBSpline &	J2
	)

Definition at line 405 of file soecp_eval_reference.cpp.

References cart2sph(), qmcplusplus::conj(), lMatrixElement(), PI, quad, sMatrixElement(), TWF(), Wso(), wt, and Ylm().

Referenced by calcVal().

{
  std::complex<double> angint(0.0, 0.0);
  for (int i = 0; i < quad.size(); i++)
  {
    std::vector<double> sph2 = cart2sph(quad[i]);
    sph2[0] *= positions[iel][0]; //now scaled to appropriate distance

    std::vector<std::vector<double>> newpositions = positions;
    std::vector<double> newspins                  = spins;
    newpositions[iel]                             = sph2;
    newspins[iel]                                 = s2;

    std::complex<double> integrand(0.0, 0.0);
    for (int l = 1; l <= 3; l++)
    {
      std::complex<double> msum(0.0, 0.0);
      for (int m1 = -l; m1 <= l; m1++)
      {
        for (int m2 = -l; m2 <= l; m2++)
        {
          std::complex<double> ldots(0.0, 0.0);
          for (int d = 0; d < 3; d++)
            ldots += lMatrixElement(l, m1, m2, d) * sMatrixElement(spins[iel], newspins[iel], d);
          msum += Ylm(l, m1, positions[iel]) * std::conj(Ylm(l, m2, newpositions[iel])) * ldots;
        }
      }
      integrand += Wso(l, positions[iel][0]) * msum;
    }

    std::complex<double> num   = TWF(newpositions, newspins, J2);
    std::complex<double> denom = TWF(positions, spins, J2);
    integrand *= num / denom;
    angint += integrand * wt[i] * 4.0 * PI;
  }
  return angint;
}

calcSOECP()

void calcSOECP

(

int

npts

)

Definition at line 476 of file soecp_eval_reference.cpp.

References calcVal(), cart2sph(), and qmcplusplus::imag().

{
  //2el system, atom at origin, 2body jastrow only
  std::vector<double> pos_e1                 = {0.138, -0.24, 0.216};
  std::vector<double> pos_e2                 = {-0.216, 0.24, -0.138};
  std::vector<std::vector<double>> positions = {cart2sph(pos_e1), cart2sph(pos_e2)};
  double spin_e1                             = 0.2;
  double spin_e2                             = 0.51;
  std::vector<double> spins                  = {spin_e1, spin_e2};

  double rcut                = 5.0;
  double cusp                = -0.25;
  std::vector<double> coeffs = {0.02904699284, -0.1004179, -0.1752703883, -0.2232576505, -0.2728029201};
  cubicBSpline J2(rcut, cusp, coeffs);

  std::complex<double> val;
  for (int iel = 0; iel < 2; iel++)
    val += calcVal(npts, iel, positions, spins, J2);
  std::cout << npts << " " << std::setprecision(10) << std::real(val) << " " << std::imag(val) << std::endl;
}

calcVal()

std::complex<double> calcVal	(	int	npts,
		int	iel,
		std::vector< std::vector< double >> &	positions,
		std::vector< double > &	spins,
		const cubicBSpline &	J2
	)

Definition at line 447 of file soecp_eval_reference.cpp.

References calcAngInt(), and PI.

Referenced by calcSOECP().

{
  std::complex<double> sint(0.0, 0.0);
  double smin = 0.0;
  double smax = 2 * PI;
  double h    = (smax - smin) / npts;
  for (int k = 1; k <= npts - 1; k += 2)
  {
    double s2                   = smin + k * h;
    std::complex<double> angint = calcAngInt(iel, s2, positions, spins, J2);
    sint += 4 * h / 3. * angint;
  }
  for (int k = 2; k <= npts - 2; k += 2)
  {
    double s2                   = smin + k * h;
    std::complex<double> angint = calcAngInt(iel, s2, positions, spins, J2);
    sint += 2 * h / 3. * angint;
  }
  sint += h / 3. * calcAngInt(iel, smin, positions, spins, J2);
  sint += h / 3. * calcAngInt(iel, smax, positions, spins, J2);
  sint /= (2.0 * PI);

  return sint;
}

cart2sph()

std::vector<double> cart2sph

(

const std::vector< double > &

cart

)

Definition at line 276 of file soecp_eval_reference.cpp.

References qmcplusplus::acos(), qmcplusplus::atan2(), and qmcplusplus::sqrt().

Referenced by calcAngInt(), calcSOECP(), and evaluateValueAndDerivative().

{
  //Physics convention
  std::vector<double> sph(3);
  sph[0] = std::sqrt(cart[0] * cart[0] + cart[1] * cart[1] + cart[2] * cart[2]);
  sph[1] = std::acos(cart[2] / sph[0]);
  sph[2] = std::atan2(cart[1], cart[0]);
  return sph;
}

chid()

std::complex<double> chid

(

double

s

)

Definition at line 298 of file soecp_eval_reference.cpp.

References qmcplusplus::exp(), I, and qmcplusplus::Units::time::s.

Referenced by spinor0(), and spinor1().

{ return std::exp(-I * s); }

chiu()

std::complex<double> chiu

(

double

s

)

Definition at line 297 of file soecp_eval_reference.cpp.

References qmcplusplus::exp(), I, and qmcplusplus::Units::time::s.

Referenced by spinor0(), and spinor1().

{ return std::exp(I * s); }

dratioTWF()

std::complex<double> dratioTWF	(	const int	ip,
		std::vector< std::vector< double >> &	positions,
		std::vector< double > &	spins,
		cubicBSpline &	J2
	)

Definition at line 390 of file soecp_eval_reference.cpp.

References cubicBSpline::parmDeriv(), sph2cart(), and qmcplusplus::sqrt().

Referenced by evaluateValueAndDerivative().

{
  std::vector<double> cart0 = sph2cart(positions[0]);
  std::vector<double> cart1 = sph2cart(positions[1]);
  double dx                 = (cart1[0] - cart0[0]);
  double dy                 = (cart1[1] - cart0[1]);
  double dz                 = (cart1[2] - cart0[2]);
  double r12                = std::sqrt(dx * dx + dy * dy + dz * dz);
  //Since the det doesn't have any parameters to optimize, the dratioTWF is only -dJ
  return -J2.parmDeriv(ip, r12);
}

evaluateValueAndDerivative()

void evaluateValueAndDerivative

(

)

Definition at line 497 of file soecp_eval_reference.cpp.

References cart2sph(), qmcplusplus::conj(), dratioTWF(), lMatrixElement(), norm(), PI, quad, sMatrixElement(), TWF(), Wso(), wt, and Ylm().

Referenced by main().

{
  const int npts = 8; //number of spin integral pts
  //2el system, atom at origin, 2body jastrow only
  std::vector<double> pos_e1                 = {0.138, -0.24, 0.216};
  std::vector<double> pos_e2                 = {-0.216, 0.24, -0.138};
  std::vector<std::vector<double>> positions = {cart2sph(pos_e1), cart2sph(pos_e2)};
  double spin_e1                             = 0.2;
  double spin_e2                             = 0.51;
  std::vector<double> spins                  = {spin_e1, spin_e2};

  double rcut                = 5.0;
  double cusp                = -0.25;
  std::vector<double> coeffs = {0.02904699284, -0.1004179, -0.1752703883, -0.2232576505, -0.2728029201};
  cubicBSpline spl(rcut, cusp, coeffs);

  //these are the reference values from evaluateDerivatives, which gets dhpsioverpsi for the
  //kinetic energy component only. Going to print the accumulated value for the SOECP component
  //and the kinetic energy and print to add to unit test in test_ecp
  std::vector<double> dhpsioverpsi = {-3.807451601, 0.1047251267, 3.702726474, 0, 0};

  std::complex<double> soecp_value;
  //loop over each particle
  for (int iel = 0; iel < 2; iel++)
  {
    //set positions, spins, and weights for all integral points for particle iel
    std::vector<std::vector<std::vector<double>>> positions_vp;
    std::vector<std::vector<double>> spins_vp;
    std::vector<double> spin_quad_wts;

    auto add_quad_points = [&](const double spinval, const double spin_wt) {
      for (int i = 0; i < quad.size(); i++)
      {
        std::vector<double> sph2 = cart2sph(quad[i]);
        sph2[0] *= positions[iel][0]; //now scaled to appropriate distance
        positions_vp.push_back(positions);
        spins_vp.push_back(spins);
        positions_vp.back()[iel] = sph2;
        spins_vp.back()[iel]     = spinval;
        double spin_norm         = 1.0 / (2.0 * PI);
        double spatial_norm      = 4.0 * PI;
        double norm              = spin_norm * spatial_norm;
        spin_quad_wts.push_back(spin_wt * wt[i] * norm);
      }
    };
    double smin = 0.0;
    double smax = 2 * PI;
    double h    = (smax - smin) / npts;
    for (int k = 1; k <= npts - 1; k += 2)
    {
      double newspin = smin + k * h;
      double spinwt  = 4. * h / 3.;
      add_quad_points(newspin, spinwt);
    }
    for (int k = 2; k <= npts - 2; k += 2)
    {
      double newspin = smin + k * h;
      double spinwt  = 2. * h / 3.;
      add_quad_points(newspin, spinwt);
    }
    add_quad_points(smin, h / 3.);
    add_quad_points(smax, h / 3.);
    assert(positions_vp.size() == spins_vp.size());
    assert(spins_vp.size() == spin_quad_wts.size());
    assert(spin_quad_wts.size() == (npts + 1) * quad.size());

    std::vector<std::complex<double>> dlogpsi(coeffs.size());
    for (int ip = 0; ip < dlogpsi.size(); ip++)
      dlogpsi[ip] = dratioTWF(ip, positions, spins, spl);
    if (iel == 0) //only print for first electron
    {
      std::cout << "std::vector<RealType> dlogpsi_refs = { ";
      for (int ip = 0; ip < dlogpsi.size(); ip++)
      {
        double real = std::real(dlogpsi[ip]);
        (ip != dlogpsi.size() - 1) ? std::cout << std::setprecision(10) << real << ", "
                                   : std::cout << std::setprecision(10) << real << "};" << std::endl;
      }
    }

    std::vector<std::vector<std::complex<double>>> dratio;
    for (int iq = 0; iq < spin_quad_wts.size(); iq++)
    {
      std::vector<std::complex<double>> dratio_row(coeffs.size());
      for (int ip = 0; ip < dratio_row.size(); ip++)
        dratio_row[ip] = dratioTWF(ip, positions_vp[iq], spins_vp[iq], spl) - dlogpsi[ip];
      dratio.push_back(dratio_row);
    }

    //Now we can evaluate the contribution to the potential from each quadrature/spin point
    auto SOECP_spin_quad = [&](std::vector<std::vector<double>>& newpos, std::vector<std::vector<double>>& oldpos,
                               std::vector<double>& newspin, std::vector<double>& oldspin, const cubicBSpline& J2) {
      std::complex<double> integrand;
      for (int l = 1; l <= 3; l++)
      {
        std::complex<double> msum;
        for (int m1 = -l; m1 <= l; m1++)
        {
          for (int m2 = -l; m2 <= l; m2++)
          {
            std::complex<double> ldots;
            for (int id = 0; id < 3; id++)
              ldots += lMatrixElement(l, m1, m2, id) * sMatrixElement(oldspin[iel], newspin[iel], id);
            msum += Ylm(l, m1, oldpos[iel]) * std::conj(Ylm(l, m2, newpos[iel])) * ldots;
          }
        }
        integrand += Wso(l, oldpos[iel][0]) * msum;
      }
      std::complex<double> psinew = TWF(newpos, newspin, J2);
      std::complex<double> psiold = TWF(oldpos, oldspin, J2);
      integrand *= psinew / psiold;
      return integrand;
    };

    std::vector<std::complex<double>> wvec(spin_quad_wts.size());
    for (int iq = 0; iq < spin_quad_wts.size(); iq++)
    {
      wvec[iq] = spin_quad_wts[iq] * SOECP_spin_quad(positions_vp[iq], positions, spins_vp[iq], spins, spl);
      soecp_value += wvec[iq];
    }

    //Now do final evaluation of dhpsioverpsi, using a gemv
    for (int ip = 0; ip < dhpsioverpsi.size(); ip++)
    {
      std::complex<double> sum;
      for (int iq = 0; iq < wvec.size(); iq++)
        sum += wvec[iq] * dratio[iq][ip];
      dhpsioverpsi[ip] += std::real(sum);
    }
  } // loop over each particle
  std::cout << std::setprecision(10) << "SOECP Value: " << soecp_value << std::endl;
  std::cout << "std::vector<RealType> dhpsioverpsi_refs = { ";
  for (int ip = 0; ip < dhpsioverpsi.size(); ip++)
    (ip != dhpsioverpsi.size() - 1) ? std::cout << std::setprecision(10) << dhpsioverpsi[ip] << ", "
                                    : std::cout << std::setprecision(10) << dhpsioverpsi[ip] << "};" << std::endl;
}

kroneckerDelta()

int kroneckerDelta	(	int	x,
		int	y
	)

Definition at line 342 of file soecp_eval_reference.cpp.

Referenced by lMatrixElement().

{ return (x == y) ? 1 : 0; }

lMatrixElement()

std::complex<double> lMatrixElement	(	int	l,
		int	m1,
		int	m2,
		int	d
	)

Definition at line 345 of file soecp_eval_reference.cpp.

References kroneckerDelta(), and qmcplusplus::sqrt().

Referenced by calcAngInt(), and evaluateValueAndDerivative().

{
  double pref1, pref2, val;
  switch (d)
  {
  case 0:
    pref1 = std::sqrt(l * (l + 1) - m2 * (m2 + 1));
    pref2 = std::sqrt(l * (l + 1) - m2 * (m2 - 1));
    val   = 0.5 * (pref1 * kroneckerDelta(m1, m2 + 1) + pref2 * kroneckerDelta(m1, m2 - 1));
    return std::complex<double>(val, 0.0);
    break;
  case 1:
    pref1 = std::sqrt(l * (l + 1) - m2 * (m2 - 1));
    pref2 = std::sqrt(l * (l + 1) - m2 * (m2 + 1));
    val   = 0.5 * (pref1 * kroneckerDelta(m1, m2 - 1) - pref2 * kroneckerDelta(m1, m2 + 1));
    return std::complex<double>(0.0, val);
    break;
  case 2:
    return std::complex<double>(m2 * kroneckerDelta(m1, m2), 0.0);
    break;
  default:
    exit(1);
    break;
  }
}

main()

int main

(

)

Definition at line 634 of file soecp_eval_reference.cpp.

References evaluateValueAndDerivative().

{
  //used to validate bspline jastrow implementation
  //testJastrow();
  //for (int n = 2; n <= 100; n += 2)
  //  calcSOECP(n);
  //looks converged by n=8 if you uncomment above,, used in evaluateValueAndDerivative
  evaluateValueAndDerivative();
}

orbd()

std::complex<double> orbd

(

const std::vector< double > &

sph

)

Definition at line 306 of file soecp_eval_reference.cpp.

References qmcplusplus::exp(), I, and sph2cart().

Referenced by spinor0(), and spinor1().

{
  std::vector<double> cart = sph2cart(sph);
  return std::exp(2. * I * (cart[0] + cart[1] + cart[2]));
}

orbu()

std::complex<double> orbu

(

const std::vector< double > &

sph

)

Definition at line 301 of file soecp_eval_reference.cpp.

References qmcplusplus::exp(), I, and sph2cart().

Referenced by spinor0(), and spinor1().

{
  std::vector<double> cart = sph2cart(sph);
  return std::exp(I * (cart[0] + cart[1] + cart[2]));
}

sMatrixElement()

std::complex<double> sMatrixElement	(	double	s1,
		double	s2,
		int	d
	)

Definition at line 323 of file soecp_eval_reference.cpp.

References qmcplusplus::cos(), and qmcplusplus::sin().

Referenced by calcAngInt(), and evaluateValueAndDerivative().

{
  switch (d)
  {
  case 0:
    return std::complex<double>(std::cos(s1 + s2), 0.0);
    break;
  case 1:
    return std::complex<double>(std::sin(s1 + s2), 0.0);
    break;
  case 2:
    return std::complex<double>(0.0, std::sin(s1 - s2));
    break;
  default:
    exit(1);
    break;
  }
}

sph2cart()

std::vector<double> sph2cart

(

const std::vector< double > &

sph

)

Definition at line 286 of file soecp_eval_reference.cpp.

References qmcplusplus::cos(), and qmcplusplus::sin().

Referenced by dratioTWF(), orbd(), orbu(), and TWF().

{
  //Physics convention
  std::vector<double> cart(3);
  cart[0] = sph[0] * std::sin(sph[1]) * std::cos(sph[2]);
  cart[1] = sph[0] * std::sin(sph[1]) * std::sin(sph[2]);
  cart[2] = sph[0] * std::cos(sph[1]);
  return cart;
}

spinor0()

std::complex<double> spinor0	(	const std::vector< double > &	sph,
		double	s
	)

Definition at line 313 of file soecp_eval_reference.cpp.

References chid(), chiu(), orbd(), orbu(), and qmcplusplus::Units::time::s.

Referenced by TWF().

{
  return orbu(sph) * chiu(s) + orbd(sph) * chid(s);
}

spinor1()

std::complex<double> spinor1	(	const std::vector< double > &	sph,
		double	s
	)

Definition at line 317 of file soecp_eval_reference.cpp.

References chid(), chiu(), orbd(), orbu(), and qmcplusplus::Units::time::s.

Referenced by TWF().

{
  return orbd(sph) * chiu(s) + orbu(sph) * chid(s);
}

testJastrow()

void testJastrow

(

)

Definition at line 113 of file soecp_eval_reference.cpp.

References qmcplusplus::abs(), and cubicBSpline::getVal().

{
  //testing 1b cusp from gen_bspline_jastrow.py
  {
    const double rcut = 10;
    const double cusp = 2.0;

    std::vector<double> coeffs = {-0.2032153051,  -0.1625595974,  -0.143124599,   -0.1216434956,
                                  -0.09919771951, -0.07111729038, -0.04445345869, -0.02135082917};

    cubicBSpline spl(rcut, cusp, coeffs);

    //from gen_bspline_jastrow.py
    std::vector<double> refVals = {-0.9304041433,   -0.252599792,   -0.1637586749,  -0.1506226948,  -0.1394848415,
                                   -0.128023472,    -0.1161729491,  -0.1036884223,  -0.08992443283, -0.07519614609,
                                   -0.06054074137,  -0.04654631918, -0.03347994129, -0.0211986378,  -0.01004416026,
                                   -0.002594125744, -0.000166024047};

    std::vector<double> rvals(refVals.size());
    for (int i = 0; i < refVals.size(); i++)
      rvals[i] = i * 0.60;

    bool good = true;
    for (int i = 0; i < refVals.size(); i++)
      if (std::abs(refVals[i] - spl.getVal(rvals[i])) > 1E-6)
      {
        std::cout << "error: " << rvals[i] << " " << refVals[i] << " != " << spl.getVal(rvals[i]) << std::endl;
        good = false;
      }
    std::cout << "1b jastrow w/ cusp correct: " << good << std::endl;
  }

  //test J2 from gen_bspline_jastrow.py
  {
    const double rcut = 10;
    const double cusp = -0.5;

    std::vector<double> coeffs = {0.02904699284, -0.1004179,    -0.1752703883, -0.2232576505, -0.2728029201,
                                  -0.3253286875, -0.3624525145, -0.3958223107, -0.4268582166, -0.4394531176};

    cubicBSpline spl(rcut, cusp, coeffs);

    //from gen_bspline_jastrow.py
    std::vector<double> refVals = {0.1374071801,   -0.04952403966, -0.121361995,  -0.1695590431, -0.2058414025,
                                   -0.2382237097,  -0.2712606182,  -0.3047843679, -0.3347515004, -0.3597048574,
                                   -0.3823503292,  -0.4036800017,  -0.4219818468, -0.4192355508, -0.3019238309,
                                   -0.09726352421, -0.006239062395};


    double r   = 0.6935647049;
    double val = spl.getVal(r);
    std::cout << r << " " << val << std::endl;
    std::vector<double> rvals(refVals.size());
    for (int i = 0; i < refVals.size(); i++)
      rvals[i] = i * 0.60;

    bool good = true;
    for (int i = 0; i < refVals.size(); i++)
      if (std::abs(refVals[i] - spl.getVal(rvals[i])) > 1E-6)
      {
        std::cout << "error: " << rvals[i] << " " << refVals[i] << " != " << spl.getVal(rvals[i]) << std::endl;
        good = false;
      }
    std::cout << "2b jastrow correct: " << good << std::endl;
  }
}

TWF()

std::complex<double> TWF	(	std::vector< std::vector< double >> &	positions,
		std::vector< double > &	spins,
		const cubicBSpline &	J2
	)

Definition at line 375 of file soecp_eval_reference.cpp.

References qmcplusplus::det(), qmcplusplus::exp(), cubicBSpline::getVal(), sph2cart(), spinor0(), spinor1(), and qmcplusplus::sqrt().

Referenced by calcAngInt(), and evaluateValueAndDerivative().

{
  std::vector<double> cart0 = sph2cart(positions[0]);
  std::vector<double> cart1 = sph2cart(positions[1]);
  double dx                 = (cart1[0] - cart0[0]);
  double dy                 = (cart1[1] - cart0[1]);
  double dz                 = (cart1[2] - cart0[2]);
  double r12                = std::sqrt(dx * dx + dy * dy + dz * dz);
  std::complex<double> det  = spinor0(positions[0], spins[0]) * spinor1(positions[1], spins[1]) -
      spinor1(positions[0], spins[0]) * spinor0(positions[1], spins[1]);
  return det * exp(-J2.getVal(r12));
}

Wso()

double Wso	(	int	l,
		double	r
	)

Definition at line 373 of file soecp_eval_reference.cpp.

References qmcplusplus::exp().

Referenced by calcAngInt(), and evaluateValueAndDerivative().

{ return exp(-l * r * r); }

Ylm()

std::complex<double> Ylm	(	int	l,
		int	m,
		const std::vector< double > &	sph
	)

Definition at line 182 of file soecp_eval_reference.cpp.

References qmcplusplus::cos(), qmcplusplus::exp(), I, qmcplusplus::Units::distance::m, PI, qmcplusplus::pow(), qmcplusplus::sin(), and qmcplusplus::sqrt().

Referenced by SCTFunctor< SCT, L >::apply(), SCTFunctor< SCT, 1 >::apply(), SCTFunctor< SCT, 2 >::apply(), SCTFunctor< SCT, 3 >::apply(), calcAngInt(), SphericalTensor< T, Point_t, Tensor_t, GGG_t >::evaluate(), SoaSphericalTensor< ST >::evaluate_bare(), SphericalTensor< T, Point_t, Tensor_t, GGG_t >::evaluateAll(), SphericalTensor< T, Point_t, Tensor_t, GGG_t >::evaluateTest(), evaluateValueAndDerivative(), SoaSphericalTensor< ST >::evaluateVGL_impl(), and SphericalTensor< T, Point_t, Tensor_t, GGG_t >::evaluateWithHessian().

{
  //From wiki
  double pref;
  switch (l)
  {
  case 0:
    return std::complex<double>(0.5 * std::sqrt(1. / PI), 0.0);
    break;
  case 1:
    switch (m)
    {
    case -1:
      pref = 0.5 * std::sqrt(3. / (2 * PI));
      return pref * std::exp(-I * sph[2]) * std::sin(sph[1]);
      break;
    case 0:
      pref = 0.5 * std::sqrt(3. / (PI));
      return pref * std::cos(sph[1]);
      break;
    case 1:
      pref = -0.5 * std::sqrt(3. / (2 * PI));
      return pref * std::exp(I * sph[2]) * std::sin(sph[1]);
      break;
    default:
      exit(1);
      break;
    }
    break;
  case 2:
    switch (m)
    {
    case -2:
      pref = 0.25 * std::sqrt(15. / (2 * PI));
      return pref * std::exp(-2.0 * I * sph[2]) * std::sin(sph[1]) * std::sin(sph[1]);
      break;
    case -1:
      pref = 0.5 * std::sqrt(15. / (2 * PI));
      return pref * std::exp(-I * sph[2]) * std::sin(sph[1]) * std::cos(sph[1]);
      break;
    case 0:
      pref = 0.25 * std::sqrt(5. / (PI));
      return pref * (3. * std::cos(sph[1]) * std::cos(sph[1]) - 1.);
      break;

    case 1:
      pref = -0.5 * std::sqrt(15. / (2 * PI));
      return pref * std::exp(I * sph[2]) * std::sin(sph[1]) * std::cos(sph[1]);
      break;
    case 2:
      pref = 0.25 * std::sqrt(15. / (2 * PI));
      return pref * std::exp(2.0 * I * sph[2]) * std::sin(sph[1]) * std::sin(sph[1]);
      break;
    }
  case 3:
    switch (m)
    {
    case -3:
      pref = 0.125 * std::sqrt(35. / PI);
      return pref * std::exp(-3. * I * sph[2]) * std::pow(std::sin(sph[1]), 3);
      break;
    case -2:
      pref = 0.25 * std::sqrt(105. / (2 * PI));
      return pref * std::exp(-2. * I * sph[2]) * std::pow(std::sin(sph[1]), 2) * std::cos(sph[1]);
      break;
    case -1:
      pref = 0.125 * std::sqrt(21. / PI);
      return pref * std::exp(-I * sph[2]) * std::sin(sph[1]) * (5 * std::pow(std::cos(sph[1]), 2) - 1);
      break;
    case 0:
      pref = 0.25 * std::sqrt(7. / PI);
      return pref * (5 * std::pow(std::cos(sph[1]), 3) - 3. * std::cos(sph[1]));
      break;
    case 1:
      pref = -0.125 * std::sqrt(21. / PI);
      return pref * std::exp(I * sph[2]) * std::sin(sph[1]) * (5 * std::pow(std::cos(sph[1]), 2) - 1);
      break;
    case 2:
      pref = 0.25 * std::sqrt(105. / (2 * PI));
      return pref * std::exp(2. * I * sph[2]) * std::pow(std::sin(sph[1]), 2) * std::cos(sph[1]);
      break;
    case 3:
      pref = -0.125 * std::sqrt(35. / PI);
      return pref * std::exp(3. * I * sph[2]) * std::pow(std::sin(sph[1]), 3);
      break;
    }
    break;
  default:
    exit(1);
    break;
  }
}

Variable Documentation

I

std::complex<double> I = std::complex<double>(0, 1)

Definition at line 180 of file soecp_eval_reference.cpp.

Referenced by CountingGaussianRegion::acceptMove(), CountingJastrow< RegionType >::acceptMove(), SymmetryBuilder::buildI(), MultiDiracDeterminant::buildTableMatrix_calculateRatios_impl(), chid(), chiu(), CountingJastrowBuilder::createCJ(), CountingGaussianRegion::dLval_saved(), CountingGaussianRegion::evaluate(), CountingGaussianRegion::evaluateDerivatives(), CountingJastrow< RegionType >::evaluateDerivatives(), CountingJastrow< RegionType >::evaluateExponents(), CountingGaussianRegion::evaluateTemp(), CountingGaussianRegion::evaluateTempDerivatives(), CountingJastrow< RegionType >::evaluateTempExponents(), CountingGaussianRegion::getVars(), CountingGaussianRegion::initialize(), CountingJastrow< RegionType >::initialize(), MultiDiracDeterminant::mw_buildTableMatrix_calculateRatios_impl(), orbd(), orbu(), CountingGaussianRegion::reportStatus(), CountingGaussianRegion::restore(), CountingGaussianRegion::total_num_derivs(), and Ylm().

PI

double PI = 4. * std::atan(1)

Definition at line 8 of file soecp_eval_reference.cpp.

Referenced by calcAngInt(), calcVal(), evaluateValueAndDerivative(), and Ylm().

quad

std::vector<std::vector<double> > quad

Initial value:

= {{1, 0, 0},
                                         {-1, 1.224646853e-16, 0},
                                         {0.4472135901, 0.8944271803, 0},
                                         {-0.4472135901, 0.7236068249, 0.5257310867},
                                         {0.4472135901, 0.2763932049, 0.8506507874},
                                         {-0.4472135901, -0.2763932049, 0.8506507874},
                                         {0.4472135901, -0.7236068249, 0.5257310867},
                                         {-0.4472135901, -0.8944271803, 1.095357398e-16},
                                         {0.4472135901, -0.7236068249, -0.5257310867},
                                         {-0.4472135901, -0.2763932049, -0.8506507874},
                                         {0.4472135901, 0.2763932049, -0.8506507874},
                                         {-0.4472135901, 0.7236068249, -0.5257310867}}

Definition at line 11 of file soecp_eval_reference.cpp.

Referenced by calcAngInt(), and evaluateValueAndDerivative().

wt

std::vector<double> wt

Initial value:

= {0.08333333582, 0.08333333582, 0.08333333582, 0.08333333582, 0.08333333582, 0.08333333582,
                          0.08333333582, 0.08333333582, 0.08333333582, 0.08333333582, 0.08333333582, 0.08333333582}

Definition at line 23 of file soecp_eval_reference.cpp.

Referenced by calcAngInt(), and evaluateValueAndDerivative().