diff --git a/nexus/documentation/nexus_user_guide.pdf b/nexus/documentation/nexus_user_guide.pdf deleted file mode 100644 index 124a15d75..000000000 Binary files a/nexus/documentation/nexus_user_guide.pdf and /dev/null differ diff --git a/nexus/documentation/user_guide_source/README.md b/nexus/documentation/user_guide_source/README.md new file mode 100644 index 000000000..761f7e01f --- /dev/null +++ b/nexus/documentation/user_guide_source/README.md @@ -0,0 +1,9 @@ +# NEXUS User Guide + +This directory contains the LaTeX source for the NEXUS user guide. The +PDF guide is not currently created by the build process and must be +created manually. + +- A script, build_nexus_user_guide.sh provides the current "best" version of the user guide. +- TeX Live 2017 or equivalent, or more recent versions are recommended +- Python and the python module pygmentize are required for the source listings. diff --git a/nexus/documentation/user_guide_source/build_nexus_user_guide.sh b/nexus/documentation/user_guide_source/build_nexus_user_guide.sh new file mode 100755 index 000000000..efb6d83bd --- /dev/null +++ b/nexus/documentation/user_guide_source/build_nexus_user_guide.sh @@ -0,0 +1,24 @@ +#!/bin/sh +if [ -z "$(command -v pdflatex)" ]; then + echo "Error: pdflatex is required for build_nexus_user_guide.sh" + exit 1 +fi + +# From QMCPACK build_manual.sh, for eventual versioning +#if [ $# -eq 1 ]; then +# QMCPACK_VER="$1" +# cp -p version.tex version.save.tex +# sed -i "s/Development Version/$QMCPACK_VER/" version.tex +#fi + + +pdflatex -shell-escape nexus_user_guide.tex +pdflatex -shell-escape nexus_user_guide.tex +pdflatex -shell-escape nexus_user_guide.tex + +# From QMCPACK build_manual.sh, for eventual versioning +#if [ ! -z "$QMCPACK_VER" ]; then +# mv version.save.tex version.tex +# echo Added QMCPACK version $QMCPACK_VER +#fi +exit 0 diff --git a/nexus/documentation/user_guide_source/nexus_user_guide.idx b/nexus/documentation/user_guide_source/nexus_user_guide.idx deleted file mode 100644 index e69de29bb..000000000 diff --git a/nexus/documentation/user_guide_source/nexus_user_guide.ilg b/nexus/documentation/user_guide_source/nexus_user_guide.ilg deleted file mode 100644 index e366a40cd..000000000 --- a/nexus/documentation/user_guide_source/nexus_user_guide.ilg +++ /dev/null @@ -1,6 +0,0 @@ -This is makeindex, version 2.15 [TeX Live 2009] (kpathsea + Thai support). -Scanning input file project_suite_user_guide.idx....done (7 entries accepted, 0 rejected). -Sorting entries....done (22 comparisons). -Generating output file project_suite_user_guide.ind....done (23 lines written, 0 warnings). -Output written in project_suite_user_guide.ind. -Transcript written in project_suite_user_guide.ilg. diff --git a/nexus/documentation/user_guide_source/nexus_user_guide.ind b/nexus/documentation/user_guide_source/nexus_user_guide.ind deleted file mode 100644 index 4e68b252e..000000000 --- a/nexus/documentation/user_guide_source/nexus_user_guide.ind +++ /dev/null @@ -1,23 +0,0 @@ -\begin{theindex} - - \item basic\_qmc, \hyperpage{8} - - \indexspace - - \item generate\_physical\_system, \hyperpage{8} - - \indexspace - - \item load\_analyzer\_image, \hyperpage{8} - - \indexspace - - \item run\_project, \hyperpage{8} - - \indexspace - - \item settings, \hyperpage{8} - \item standard\_qmc, \hyperpage{8} - \item Structure, \hyperpage{8} - -\end{theindex} diff --git a/nexus/documentation/user_guide_source/nexus_user_guide.pdf b/nexus/documentation/user_guide_source/nexus_user_guide.pdf deleted file mode 100644 index 124a15d75..000000000 Binary files a/nexus/documentation/user_guide_source/nexus_user_guide.pdf and /dev/null differ diff --git a/nexus/documentation/user_guide_source/nexus_user_guide.toc b/nexus/documentation/user_guide_source/nexus_user_guide.toc deleted file mode 100644 index 9baa38c7e..000000000 --- a/nexus/documentation/user_guide_source/nexus_user_guide.toc +++ /dev/null @@ -1,54 +0,0 @@ -\contentsline {chapter}{Contents}{ii}{section*.1} -\contentsline {chapter}{\chapternumberline {1}Using this document}{1}{chapter.1} -\contentsline {chapter}{\chapternumberline {2}Overview of Nexus}{2}{chapter.2} -\contentsline {section}{\numberline {2.1}What Nexus is}{2}{section.2.1} -\contentsline {section}{\numberline {2.2}What Nexus can do}{2}{section.2.2} -\contentsline {section}{\numberline {2.3}How Nexus is used}{2}{section.2.3} -\contentsline {chapter}{\chapternumberline {3}Nexus Installation}{3}{chapter.3} -\contentsline {section}{\numberline {3.1}Setting environment variables}{3}{section.3.1} -\contentsline {section}{\numberline {3.2}Installing Python dependencies}{3}{section.3.2} -\contentsline {chapter}{\chapternumberline {4}Nexus User Scripts}{4}{chapter.4} -\contentsline {section}{\numberline {4.1}Nexus imports}{4}{section.4.1} -\contentsline {section}{\numberline {4.2}Nexus settings: global state and user-specific information}{5}{section.4.2} -\contentsline {section}{\numberline {4.3}Physical system specificaton}{6}{section.4.3} -\contentsline {section}{\numberline {4.4}Workflow specification}{7}{section.4.4} -\contentsline {subsection}{Generating simulation objects}{7}{section*.2} -\contentsline {subsection}{Composing workflows from simulation objects}{11}{section*.3} -\contentsline {section}{\numberline {4.5}Workflow execution}{13}{section.4.5} -\contentsline {section}{\numberline {4.6}Data analysis}{14}{section.4.6} -\contentsline {chapter}{\chapternumberline {5}Complete Examples}{15}{chapter.5} -\contentsline {section}{\numberline {5.1}Example 1: Bulk Diamond VMC}{17}{section.5.1} -\contentsline {section}{\numberline {5.2}Example 2: Graphene Sheet DMC}{23}{section.5.2} -\contentsline {section}{\numberline {5.3}Example 3: C 20 Molecule DMC}{33}{section.5.3} -\contentsline {section}{\numberline {5.4}Example 4: Automated oxygen dimer binding curve}{41}{section.5.4} -\contentsline {chapter}{\chapternumberline {6}Recommended Reading}{48}{chapter.6} -\contentsline {section}{\numberline {6.1}Helpful Links for Installing Python Modules}{48}{section.6.1} -\contentsline {section}{\numberline {6.2}Helpful Links for Installing Electronic Structure Codes}{48}{section.6.2} -\contentsline {subsection}{PWSCF: pw.x, pw2qmcpack.x, pw2casino.x}{48}{section*.4} -\contentsline {subsection}{QMCPACK: qmcapp, qmcapp\_complex, convert4qmc, ppconvert, sqd}{48}{section*.5} -\contentsline {subsection}{Wfconvert: wfconvert}{48}{section*.6} -\contentsline {subsection}{VASP}{49}{section*.7} -\contentsline {subsection}{GAMESS}{49}{section*.8} -\contentsline {section}{\numberline {6.3}Brushing Up On Python}{49}{section.6.3} -\contentsline {subsection}{Python}{49}{section*.9} -\contentsline {subsection}{NumPy}{49}{section*.10} -\contentsline {subsection}{Matplotlib}{49}{section*.11} -\contentsline {subsection}{Scipy and H5Py}{50}{section*.12} -\contentsline {section}{\numberline {6.4}Quantum Monte Carlo: Theory and Practice}{50}{section.6.4} -\contentsline {appendix}{\chapternumberline {A}Basic Python constructs}{51}{appendix.A} -\contentsline {section}{\numberline {A.1}Intrinsic types: \texttt {int, float, str}}{51}{section.A.1} -\contentsline {section}{\numberline {A.2}Container types: \texttt {tuple, list, array, dict, obj}}{51}{section.A.2} -\contentsline {section}{\numberline {A.3}Conditional Statements: \texttt {if/elif/else}}{53}{section.A.3} -\contentsline {section}{\numberline {A.4}Iteration: \texttt {for}}{54}{section.A.4} -\contentsline {section}{\numberline {A.5}Functions: \texttt {def}, argument syntax}{55}{section.A.5} -\contentsline {appendix}{\chapternumberline {B}QMC Practice in a Nutshell}{57}{appendix.B} -\contentsline {section}{\numberline {B.1}VMC and DMC in the abstract}{57}{section.B.1} -\contentsline {section}{\numberline {B.2}From expectation values to random walks}{58}{section.B.2} -\contentsline {section}{\numberline {B.3}Quality orbitals: planewaves, cutoffs, splines, and meshes}{58}{section.B.3} -\contentsline {section}{\numberline {B.4}Quality Jastrows: less variance = more efficient}{59}{section.B.4} -\contentsline {section}{\numberline {B.5}Finite size effects: k-points, supercells, and corrections}{60}{section.B.5} -\contentsline {section}{\numberline {B.6}Imaginary time discretization: the DMC timestep}{60}{section.B.6} -\contentsline {section}{\numberline {B.7}Population control bias: safety in numbers}{61}{section.B.7} -\contentsline {section}{\numberline {B.8}The fixed node/phase approximation: varying the nodes/phase}{62}{section.B.8} -\contentsline {section}{\numberline {B.9}Pseudopotentials: theoretical dissonance, the locality approximation, and T-moves}{62}{section.B.9} -\contentsline {section}{\numberline {B.10}Other approximations: what else is missing?}{63}{section.B.10} diff --git a/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.h b/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.h index b865c076f..98487ae57 100644 --- a/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.h +++ b/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.h @@ -66,6 +66,36 @@ public: nElms=nIntPnts; } + void setPhiAndEta(SPOSetBasePtr Phi, SPOSetBasePtr Eta) + { + Psi1 = Phi; + Psi2 = Eta; + int norb = Psi1->OrbitalSetSize; + + val1.resize(norb); + grad1.resize(norb); + lapl1.resize(norb); + + val2.resize(norb); + grad2.resize(norb); + lapl2.resize(norb); + } + + void allocateELspace() + { + ELideal.resize(nElms); + ELorig.resize(nElms); + ELcurr.resize(nElms); + pos.resize(nElms); + } + + void computeValAtZero() + { + TinyVector ddr=0; + evaluate(Psi1,ddr,val1,grad1,lapl1); + valAtZero = val1[0]; + } + ~CuspCorr() {} /* @@ -285,6 +315,7 @@ private: SPOSetBasePtr Psi1,Psi2; +public: // cutoff RealType beta0,DX,eta0, ELorigAtRc; RealType Rc_init,Rc,C,sg,Z,valAtZero,valAtRc,Rc_max; diff --git a/src/QMCWaveFunctions/tests/CMakeLists.txt b/src/QMCWaveFunctions/tests/CMakeLists.txt index de68ee93f..01dee87da 100644 --- a/src/QMCWaveFunctions/tests/CMakeLists.txt +++ b/src/QMCWaveFunctions/tests/CMakeLists.txt @@ -36,8 +36,8 @@ FOREACH(fname ${FILES_TO_COPY}) EXECUTE_PROCESS(COMMAND ${CMAKE_COMMAND} -E copy "${CMAKE_CURRENT_SOURCE_DIR}/${fname}" ${UTEST_DIR}) ENDFOREACH() -IF ((NOT QMC_COMPLEX) AND (NOT QMC_MIXED_PRECISION)) - SET(MO_SRCS test_MO.cpp) +IF ((NOT QMC_COMPLEX) AND (NOT QMC_MIXED_PRECISION) AND (NOT QMC_CUDA)) + SET(MO_SRCS test_MO.cpp test_cusp_corr.cpp) ENDIF() ADD_EXECUTABLE(${UTEST_EXE} test_wf.cpp test_bspline_jastrow.cpp test_einset.cpp test_pw.cpp diff --git a/src/QMCWaveFunctions/tests/gen_cusp_corr.py b/src/QMCWaveFunctions/tests/gen_cusp_corr.py new file mode 100644 index 000000000..cf9a3b7da --- /dev/null +++ b/src/QMCWaveFunctions/tests/gen_cusp_corr.py @@ -0,0 +1,231 @@ +from __future__ import print_function + +# Cusp corrections for gaussian orbitals + +# From "Scheme for adding electron-nucleus cusps to Gaussian orbitals" A. Ma, D. Towler, N. D. Drummond, R. J. Needs, Journal of Chemical Physics 122, 224322(2005) https://doi.org/10.1063/1.1940588 + +# Also qmc_algorithms/Wavefunctions/CuspCorrection.ipynb + + +from sympy import * +import gaussian_orbitals +import read_qmcpack + + + +alpha = IndexedBase('alpha') +rc = Symbol('r_c') +X1,X2,X3,X4,X5 = symbols('X_1 X_2 X_3 X_4 X_5') +r = Symbol('r') +Zeff = Symbol('Z_eff') + + + + +def solve_for_alpha(): + p = alpha[0] + alpha[1]*r + alpha[2]*r**2 + alpha[3]*r**3 + alpha[4]*r**4 + + + # Constraint equations + # Value matches at r_c + eq1 = Eq(p.subs(r,rc), X1) + # 1st derivative matches at r_c + eq2 = Eq(diff(p,r).subs(r,rc), X2) + # 2nd derivative matches at r_c + eq3 = Eq((diff(p,r,2)+diff(p,r)**2).subs(r,rc),X3) + # Cusp condition + eq4 = Eq(diff(p,r).subs(r,0),X4) + # Value of phi tilde at r=0 + eq5 = Eq(p.subs(r,0),X5) + + sln = solve([eq1, eq2, eq3, eq4, eq5],[alpha[0], alpha[1], alpha[2], alpha[3], alpha[4]]) + + print('Symbolic solution for alpha:') + for i,s in enumerate(sln[0]): + print(alpha[i],' = ',s) + print() + + return sln[0] + +alpha_sln = solve_for_alpha() + + +def solve_for_alpha(Xvals, rc_val): + svals = {rc: rc_val, X1:Xvals[1], X2:Xvals[2], X3:Xvals[3], X4:Xvals[4], X5:Xvals[5]} + alpha_vals = [] + for i,s in enumerate(alpha_sln): + alpha_vals.append(s.subs(svals)) + return alpha_vals + +def simple_X_vals(): + rc_val = 0.1 + Xvals = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0] # first entry is unused + + print('For X: ',Xvals[1:]) + alpha_vals = solve_for_alpha(Xvals, rc_val) + for i,a in enumerate(alpha_vals): + print(alpha[i], ' = ', a) + #svals = {rc: rc_val, X1:Xvals[1], X2:Xvals[2], X3:Xvals[3], X4:Xvals[4], X5:Xvals[5]} + #for i,s in enumerate(alpha_sln): + # print(alpha[i], ' = ', s.subs(svals)) + +def evalX(gto, Z_val, rc_val): + Xvals = [0.0]*6 + val_rc, grad_rc, lap_rc = [v[0] for v in gto.eval_vgl(rc_val, 0.0, 0.0)] + val_zero, grad_zero, lap_zero = [v[0] for v in gto.eval_vgl(0.0, 0.0, 0.0)] + Xvals[1] = log(val_rc) + Xvals[2] = grad_rc[0]/val_rc + Xvals[3] = (lap_rc - 2.0*grad_rc[0]/rc_val)/val_rc + Xvals[4] = -Z_val + Xvals[5] = log(abs(val_zero)) # initially use phi at 0 + return Xvals + + + + return he_gto, Xvals, alpha_vals + +def output_required_xvals_and_alphas(Xvals, alpha_vals): + print('Xvals = ',Xvals) + print(' // From gen_cusp_corr.py') + for i in range(5): + print(' REQUIRE(X[%d] == Approx(%.15f));'%(i,Xvals[i+1])) + + print() + print(' // From gen_cusp_corr.py') + for i in range(5): + print(' REQUIRE(cusp.alpha[%d] == Approx(%.15f));'%(i,alpha_vals[i])) + +def output_array(v, name): + print(' // From gen_cusp_corr.py') + for i,a in enumerate(v): + print(' REQUIRE(%s[%d] == Approx(%.15f));'%(name, i, a)) + print() + + +def del_spherical(e, r): + """Compute Laplacian for expression e with respect to symbol r. + Currently works only with radial dependence""" + t1 = r*r*diff(e, r) + t2 = diff(t1, r)/(r*r) + return simplify(t2) + +def get_symbolic_effective_local_energy(): + p = alpha[0] + alpha[1]*r + alpha[2]*r**2 + alpha[3]*r**3 + alpha[4]*r**4 + R_sym = exp(p) + phi_tilde = R_sym + effEl_sym = -S.Half * del_spherical(phi_tilde, r)/phi_tilde - Zeff/r + return effEl_sym + + +def eval_El(El_sym, r_val, Zeff_val, alpha_vals): + slist = {alpha[0]:alpha_vals[0], alpha[1]:alpha_vals[1], alpha[2]: alpha_vals[2], alpha[3]:alpha_vals[3], + alpha[4]:alpha_vals[4], Zeff:Zeff_val, r:r_val} + val = El_sym.subs(slist).evalf() + return val + + +def get_grid(): + pos = [] + for i in range(10): + rval = .012*(i+1) + pos.append(rval) + return pos + +def get_original_local_energy(pos, gto, El_sym, alpha_vals, rc_val, Zeff_val): + vals = [] + for rval in pos: + val, grad, lap = gto.eval_vgl(rval, 0.0, 0.0) + real_el = -.5*lap[0]/val[0] - Zeff_val/rval + vals.append(real_el) + + return vals + + +def get_current_local_energy(pos, gto, El_sym, alpha_vals, rc_val, dE, Zeff_val): + vals = [] + for rval in pos: + el = eval_El(El_sym, rval, 2.0, alpha_vals) + if rval < rc_val: + vals.append(el + dE) + else: + val, grad, lap = gto.eval_vgl(rval, 0.0, 0.0) + real_el = -.5*lap[0]/val[0] - Zeff_val/rval + vals.append(real_el + dE) + + return vals + + + +def get_ideal_local_energy(pos, rc_val, beta0_val=0.0): + ideal_EL = [] + beta0 = Symbol('beta_0') + beta_vals = [beta0, 3.25819, -15.0126, 33.7308, -42.8705, 31.2276, -12.1316, 1.94692] + El_terms = [beta_vals[n]*r**(n+1) for n in range(1,8)] + Z = Symbol('Z') + Z_val = 2.0 + EL_ideal_sym = Z*Z*(beta0 + sum(El_terms)) + + slist = {beta0: 0.0, Z:Z_val, r: rc_val} + + idealEL_at_rc = EL_ideal_sym.subs(slist).evalf() + #print('idealEL at rc = ',idealEL_at_rc) + beta0_val = (-idealEL_at_rc)/Z_val/Z_val + #print('beta0_val = ',beta0_val) + + for rval in pos: + slist = {beta0: beta0_val, Z:Z_val, r: rval} + v = EL_ideal_sym.subs(slist).evalf() + ideal_EL.append(v) + + return ideal_EL + + +def values_for_He(): + rc_val = 0.1 + Z_val = 2 + + basis_set,MO = read_qmcpack.parse_qmc_wf('he_sto3g.wfj.xml') + he_gto = gaussian_orbitals.GTO(basis_set['He']) + + Xvals = evalX(he_gto, Z_val, rc_val) + + alpha_vals = solve_for_alpha(Xvals, rc_val) + output_required_xvals_and_alphas(Xvals, alpha_vals) + + El_sym = get_symbolic_effective_local_energy() + print("El_sym = ",El_sym) + + Zeff_val = 2.0 + el_at_rc = -eval_El(El_sym, rc_val, Zeff_val, alpha_vals) + dE = el_at_rc + print('el at rc_val = ',el_at_rc) + + + pos = get_grid() + + current_EL = get_current_local_energy(pos, he_gto, El_sym, alpha_vals, rc_val, dE, Zeff_val) + original_EL = get_original_local_energy(pos, he_gto, El_sym, alpha_vals, rc_val, Zeff_val) + #print('Current effective local energy') + #for (p,v) in zip(pos,current_EL): + # print(p,v) + print(" // Grid for local energy evaluations") + output_array(pos, "cusp.pos"); + + print(" // Original local energy") + output_array(original_EL, "cusp.ELorig") + + print(" // Current local energy") + output_array(current_EL, "cusp.ELcurr") + + print(" // Ideal local energy") + ideal_EL = get_ideal_local_energy(pos,rc_val) + output_array(ideal_EL, "cusp.ELideal") + + chi2 = 0.0 + for rval, ideal, curr in zip(pos, ideal_EL, current_EL): + chi2 += (ideal - curr)**2 + print(' REQUIRE(chi2 == Approx(%.10f)); '%chi2) + +if __name__ == '__main__': + #simple_X_vals() + values_for_He() diff --git a/src/QMCWaveFunctions/tests/test_cusp_corr.cpp b/src/QMCWaveFunctions/tests/test_cusp_corr.cpp new file mode 100644 index 000000000..d3eec8892 --- /dev/null +++ b/src/QMCWaveFunctions/tests/test_cusp_corr.cpp @@ -0,0 +1,235 @@ + +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2018 Jeongnim Kim and QMCPACK developers. +// +// File developed by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory +// +// File created by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + + +#include "catch.hpp" + +#include "Configuration.h" + +#include "Message/Communicate.h" +#include "Numerics/OneDimGridBase.h" +#include "Particle/DistanceTableData.h" +#include "ParticleIO/XMLParticleIO.h" +#include "Numerics/GaussianBasisSet.h" + +#include "QMCWaveFunctions/MolecularOrbitals/LocalizedBasisSet.h" +#include "QMCWaveFunctions/MolecularOrbitals/SphericalBasisSet.h" +#include "QMCWaveFunctions/MolecularOrbitals/NGOBuilder.h" +#include "QMCWaveFunctions/MolecularOrbitals/CuspCorr.h" + +#include "QMCWaveFunctions/SPOSetBuilderFactory.h" + +namespace qmcplusplus +{ + +TEST_CASE("CuspCorrection He", "[wavefunction]") +{ + OHMMS::Controller->initialize(0, NULL); + Communicate *c = OHMMS::Controller; + + ParticleSet elec; + std::vector agroup(2); + agroup[0] = 1; + agroup[1] = 1; + elec.setName("e"); + elec.create(agroup); + elec.R[0] = 0.0; + + SpeciesSet &tspecies = elec.getSpeciesSet(); + int upIdx = tspecies.addSpecies("u"); + int downIdx = tspecies.addSpecies("d"); + int massIdx = tspecies.addAttribute("mass"); + tspecies(massIdx, upIdx) = 1.0; + tspecies(massIdx, downIdx) = 1.0; + + ParticleSet ions; + ions.setName("ion0"); + ions.create(1); + ions.R[0] = 0.0; + SpeciesSet &ispecies = ions.getSpeciesSet(); + int heIdx = ispecies.addSpecies("He"); + ions.update(); + + + #ifdef ENABLE_SOA + elec.addTable(ions,DT_SOA); + #else + elec.addTable(ions,DT_AOS); + #endif + elec.update(); + + + Libxml2Document doc; + bool okay = doc.parse("he_sto3g.wfj.xml"); + REQUIRE(okay); + xmlNodePtr root = doc.getRoot(); + + TrialWaveFunction psi(c); + + OrbitalBuilderBase::PtclPoolType particle_set_map; + particle_set_map["e"] = &elec; + particle_set_map["ion0"] = &ions; + + + SPOSetBuilderFactory bf(elec, psi, particle_set_map); + + OhmmsXPathObject MO_base("//determinantset", doc.getXPathContext()); + REQUIRE(MO_base.size() == 1); + + SPOSetBuilder *bb = bf.createSPOSetBuilder(MO_base[0]); + REQUIRE(bb != NULL); + + OhmmsXPathObject slater_base("//determinant", doc.getXPathContext()); + bb->loadBasisSetFromXML(MO_base[0]); + SPOSetBase *sposet = bb->createSPOSet(slater_base[0]); + + + typedef OneDimGridBase GridType; + typedef LCOrbitalSet>, false> OrbType; + OrbType *lcob = dynamic_cast(sposet); + REQUIRE(lcob != NULL); + + + typedef CuspCorr>> CuspCorrType; + typedef OrbitalSetTraits::ValueVector_t ValueVector_t; + + double rc = 0.1; + int npts = 10; + CuspCorrType cusp(rc, npts, &elec, &ions); + + typedef NGOBuilder::CenteredOrbitalType COT; + OrbType bs_phi(lcob->myBasisSet); + bs_phi.setOrbitalSetSize(lcob->OrbitalSetSize); + bs_phi.BasisSetSize = lcob->BasisSetSize; + bs_phi.setIdentity(false); + + *(bs_phi.C) = *(lcob->C); + + OrbType bs_eta(lcob->myBasisSet); + bs_eta.setOrbitalSetSize(lcob->OrbitalSetSize); + bs_eta.BasisSetSize = lcob->BasisSetSize; + bs_eta.setIdentity(false); + *(bs_eta.C) = *(lcob->C); + // For He in a minimal basis, there are only s-type orbitals + (*bs_eta.C)(0,0) = 0.0; + + cusp.setPhiAndEta(&bs_phi, &bs_eta); + + cusp.curCenter = 0; + cusp.curOrb = 0; + cusp.Z = 2.0; + cusp.computeValAtZero(); + + ValueVector_t X(5); + cusp.evalX(X); + + //std::cout << "X = " << X << std::endl; + + // From gen_cusp_corr.py + REQUIRE(X[0] == Approx(-0.033436891110336)); + REQUIRE(X[1] == Approx(-0.653568722769692)); + REQUIRE(X[2] == Approx(-5.819488164002633)); + REQUIRE(X[3] == Approx(-2.000000000000000)); + REQUIRE(X[4] == Approx(-0.000396345019839)); + + cusp.X2alpha(X); + + // From gen_cusp_corr.py + REQUIRE(cusp.alpha[0] == Approx(-0.000396345019839)); + REQUIRE(cusp.alpha[1] == Approx(-2.000000000000000)); + REQUIRE(cusp.alpha[2] == Approx(56.659413909100188)); + REQUIRE(cusp.alpha[3] == Approx(-599.993590267020409)); + REQUIRE(cusp.alpha[4] == Approx(2003.589050855219512)); + + + cusp.allocateELspace(); + + // compute original EL + cusp.getELorig(); + + // compute current EL + cusp.getELcurr(); + + // compute ideal EL + cusp.getELideal(cusp.ELorigAtRc); + + // Grid for local energy evaluations + // From gen_cusp_corr.py + REQUIRE(cusp.pos[0] == Approx(0.012000000000000)); + REQUIRE(cusp.pos[1] == Approx(0.024000000000000)); + REQUIRE(cusp.pos[2] == Approx(0.036000000000000)); + REQUIRE(cusp.pos[3] == Approx(0.048000000000000)); + REQUIRE(cusp.pos[4] == Approx(0.060000000000000)); + REQUIRE(cusp.pos[5] == Approx(0.072000000000000)); + REQUIRE(cusp.pos[6] == Approx(0.084000000000000)); + REQUIRE(cusp.pos[7] == Approx(0.096000000000000)); + REQUIRE(cusp.pos[8] == Approx(0.108000000000000)); + REQUIRE(cusp.pos[9] == Approx(0.120000000000000)); + + // Original local energy + // From gen_cusp_corr.py + REQUIRE(cusp.ELorig[0] == Approx(-156.654088753559449)); + REQUIRE(cusp.ELorig[1] == Approx(-73.346068180623860)); + REQUIRE(cusp.ELorig[2] == Approx(-45.610385939854496)); + REQUIRE(cusp.ELorig[3] == Approx(-31.780236703094037)); + REQUIRE(cusp.ELorig[4] == Approx(-23.522092887496903)); + REQUIRE(cusp.ELorig[5] == Approx(-18.057926774366479)); + REQUIRE(cusp.ELorig[6] == Approx(-14.196956436578184)); + REQUIRE(cusp.ELorig[7] == Approx(-11.343582162638119)); + REQUIRE(cusp.ELorig[8] == Approx(-9.166698588588746)); + REQUIRE(cusp.ELorig[9] == Approx(-7.467419605354912)); + + + // Current local energy + // From gen_cusp_corr.py + REQUIRE(cusp.ELcurr[0] == Approx(-119.501377810849448)); + REQUIRE(cusp.ELcurr[1] == Approx(-84.586558430353278)); + REQUIRE(cusp.ELcurr[2] == Approx(-55.798846293994899)); + REQUIRE(cusp.ELcurr[3] == Approx(-32.804136849327627)); + REQUIRE(cusp.ELcurr[4] == Approx(-15.556451408317010)); + REQUIRE(cusp.ELcurr[5] == Approx(-4.108843171781601)); + REQUIRE(cusp.ELcurr[6] == Approx(1.506652537070609)); + REQUIRE(cusp.ELcurr[7] == Approx(1.330023830008599)); + REQUIRE(cusp.ELcurr[8] == Approx(1.387870101712975)); + REQUIRE(cusp.ELcurr[9] == Approx(3.087149084946809)); + + // Ideal local energy + // From gen_cusp_corr.py + REQUIRE(cusp.ELideal[0] == Approx(-0.080399129819489)); + REQUIRE(cusp.ELideal[1] == Approx(-0.075454680124444)); + REQUIRE(cusp.ELideal[2] == Approx(-0.067869571712510)); + REQUIRE(cusp.ELideal[3] == Approx(-0.058114416794904)); + REQUIRE(cusp.ELideal[4] == Approx(-0.046606832483210)); + REQUIRE(cusp.ELideal[5] == Approx(-0.033715675742608)); + REQUIRE(cusp.ELideal[6] == Approx(-0.019765043739301)); + REQUIRE(cusp.ELideal[7] == Approx(-0.005038047738043)); + REQUIRE(cusp.ELideal[8] == Approx(0.010219631429280)); + REQUIRE(cusp.ELideal[9] == Approx(0.025796398438142)); + + double chi2 = cusp.getchi2(); + REQUIRE(chi2 == Approx(25854.2846426019)); + + + double data[9]; + Vector xgrid; + Vector rad_orb; + xgrid.resize(1); + xgrid[0] = 0.012; + rad_orb.resize(1); + + cusp.execute(0, 0, 2.0, &bs_phi, &bs_eta, xgrid, rad_orb, "none", rc, data); + + SPOSetBuilderFactory::clear(); + +} + +}