scnc
/
qmcpack
mirror of https://github.com/QMCPACK/qmcpack.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3713 from rcclay/twf_prototype_stage2 

[Fast Derivatives] Stage 2:  Introduce TrialWavefunction Wrapper
This commit is contained in:
Ye Luo 2022-01-21 15:44:40 -06:00 committed by GitHub
parent ebd2da1746 06919bff5d
commit 60d050c70b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
15 changed files with 596 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
src/QMCHamiltonians/tests/test_ion_derivs.cpp
View File

@ -20,6 +20,7 @@
#include "QMCHamiltonians/tests/MinimalHamiltonianPool.h"
#include "ParticleIO/XMLParticleIO.h"
#include "Utilities/RandomGenerator.h"
#include "QMCWaveFunctions/TWFFastDerivWrapper.h"
namespace qmcplusplus
{
@ -766,9 +767,9 @@ TEST_CASE("Eloc_Derivatives:proto_sd_noj", "[hamiltonian]")
REQUIRE(psi != nullptr);
//end incantation
// TWFPrototype twf;
TWFFastDerivWrapper twf;
// psi->initialize_TWF_Prototype(elec, twf);
psi->initializeTWFFastDerivWrapper(elec, twf);
SPOSet::ValueVector_t values;
SPOSet::GradVector_t dpsi;
SPOSet::ValueVector_t d2psi;
@ -833,7 +834,7 @@ TEST_CASE("Eloc_Derivatives:proto_sd_noj", "[hamiltonian]")
X.push_back(upmat);
X.push_back(dnmat);
// twf.get_M(elec, matlist);
twf.getM(elec, matlist);
OperatorBase* kinop = ham.getHamiltonian(KINETIC);
@ -851,17 +852,17 @@ TEST_CASE("Eloc_Derivatives:proto_sd_noj", "[hamiltonian]")
minv.push_back(dnmat);
// twf.get_M(elec, matlist);
// twf.getM(elec, matlist);
std::vector<std::vector<ValueMatrix_t>> dB_gs;
std::vector<std::vector<ValueMatrix_t>> dM_gs;
std::vector<ValueMatrix_t> tmp_gs;
// twf.get_gs_matrix(B, B_gs);
// twf.get_gs_matrix(matlist, M_gs);
// twf.invert_M(M_gs, minv);
// twf.build_X(minv, B_gs, X);
twf.getGSMatrices(B, B_gs);
twf.getGSMatrices(matlist, M_gs);
twf.invertMatrices(M_gs, minv);
twf.buildX(minv, B_gs, X);
for (int id = 0; id < matlist.size(); id++)
{
// int ptclnum = twf.num_particles(id);
// int ptclnum = twf.numParticles(id);
int ptclnum = (id==0 ? Nup : Ndn); //hard coded until twf interface comes online.
ValueMatrix_t gs_m;
gs_m.resize(ptclnum, ptclnum);
@ -887,18 +888,18 @@ TEST_CASE("Eloc_Derivatives:proto_sd_noj", "[hamiltonian]")
{
for (int idim = 0; idim < OHMMS_DIM; idim++)
{
// twf.wipe_matrix(dB[idim]);
// twf.wipe_matrix(dM[idim]);
twf.wipeMatrices(dB[idim]);
twf.wipeMatrices(dM[idim]);
}
// twf.get_igrad_M(elec, ions, ionid, dM);
twf.getIonGradM(elec, ions, ionid, dM);
// kinop->evaluateOneBodyOpMatrixForceDeriv(elec, ions, twf, ionid, dB);
for (int idim = 0; idim < OHMMS_DIM; idim++)
{
// twf.get_gs_matrix(dB[idim], dB_gs[idim]);
// twf.get_gs_matrix(dM[idim], dM_gs[idim]);
// fkin_complex[ionid][idim] = twf.compute_gs_derivative(minv, X, dM_gs[idim], dB_gs[idim]);
twf.getGSMatrices(dB[idim], dB_gs[idim]);
twf.getGSMatrices(dM[idim], dM_gs[idim]);
fkin_complex[ionid][idim] = twf.computeGSDerivative(minv, X, dM_gs[idim], dB_gs[idim]);
}
convertToReal(fkin_complex[ionid], fkin[ionid]);
}
@ -906,7 +907,7 @@ TEST_CASE("Eloc_Derivatives:proto_sd_noj", "[hamiltonian]")
ValueType keval = 0.0;
RealType keobs = 0.0;
// keval = twf.trAB(minv, B_gs);
keval = twf.trAB(minv, B_gs);
convertToReal(keval, keobs);
// CHECK(keobs == Approx(9.1821937928e+00));
#if defined(MIXED_PRECISION)
@ -933,12 +934,12 @@ TEST_CASE("Eloc_Derivatives:proto_sd_noj", "[hamiltonian]")
app_log() << " Evaluated. Calling evaluteOneBodyOpMatrix\n";
// twf.wipe_matrix(B);
// twf.wipe_matrix(B_gs);
// twf.wipe_matrix(X);
// twf.wipeMatrices(B);
// twf.wipeMatrices(B_gs);
// twf.wipeMatrices(X);
// nlppop->evaluateOneBodyOpMatrix(elec, twf, B);
// twf.get_gs_matrix(B, B_gs);
// twf.build_X(minv, B_gs, X);
// twf.getGSMatrices(B, B_gs);
// twf.buildX(minv, B_gs, X);
ValueType nlpp = 0.0;
RealType nlpp_obs = 0.0;
@ -955,18 +956,18 @@ TEST_CASE("Eloc_Derivatives:proto_sd_noj", "[hamiltonian]")
{
for (int idim = 0; idim < OHMMS_DIM; idim++)
{
// twf.wipe_matrix(dB[idim]);
// twf.wipe_matrix(dM[idim]);
// twf.wipeMatrices(dB[idim]);
// twf.wipeMatrices(dM[idim]);
}
// twf.get_igrad_M(elec, ions, ionid, dM);
// twf.getIonGradM(elec, ions, ionid, dM);
// nlppop->evaluateOneBodyOpMatrixForceDeriv(elec, ions, twf, ionid, dB);
for (int idim = 0; idim < OHMMS_DIM; idim++)
{
// twf.get_gs_matrix(dB[idim], dB_gs[idim]);
// twf.get_gs_matrix(dM[idim], dM_gs[idim]);
// fnlpp_complex[ionid][idim] = twf.compute_gs_derivative(minv, X, dM_gs[idim], dB_gs[idim]);
// twf.getGSMatrices(dB[idim], dB_gs[idim]);
// twf.getGSMatrices(dM[idim], dM_gs[idim]);
// fnlpp_complex[ionid][idim] = twf.computeGSDerivative(minv, X, dM_gs[idim], dB_gs[idim]);
}
convertToReal(fnlpp_complex[ionid], fnlpp[ionid]);
}

1
src/QMCWaveFunctions/CMakeLists.txt
View File

@ -169,6 +169,7 @@ set(FERMION_SRCS
SPOSetBuilderFactory.cpp
TrialWaveFunction.cpp
TWFdispatcher.cpp
TWFFastDerivWrapper.cpp
WaveFunctionFactory.cpp)
if(ENABLE_CUDA)

7
src/QMCWaveFunctions/Fermion/DiracDeterminant.cpp
View File

@ -21,6 +21,7 @@
#include "CPU/SIMD/simd.hpp"
#include "Numerics/DeterminantOperators.h"
#include "Numerics/MatrixOperators.h"
#include "QMCWaveFunctions/TWFFastDerivWrapper.h"
namespace qmcplusplus
{
@ -344,6 +345,12 @@ void DiracDeterminant<DU_TYPE>::copyFromBuffer(ParticleSet& P, WFBufferType& buf
updateEng.initializeInv(psiM);
}
template<typename DU_TYPE>
void DiracDeterminant<DU_TYPE>::registerTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const
{
twf.addGroup(P, P.getGroupID(FirstIndex), Phi.get());
}
/** return the ratio only for the iat-th partcle move
* @param P current configuration
* @param iat the particle thas is being moved

4
src/QMCWaveFunctions/Fermion/DiracDeterminant.h
View File

@ -76,6 +76,10 @@ public:
void copyFromBuffer(ParticleSet& P, WFBufferType& buf) override;
/** Finds the SPOSet associated with this determinant, and registers it with WFN wrapper
*/
void registerTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const final;
/** return the ratio only for the iat-th partcle move
* @param P current configuration
* @param iat the particle thas is being moved

4
src/QMCWaveFunctions/Fermion/DiracDeterminantBase.h
View File

@ -85,6 +85,10 @@ public:
inline void reportStatus(std::ostream& os) final {}
virtual void registerTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const override
{
APP_ABORT("DiracDeterminantBase::registerTWFFastDerivWrapper must be overridden\n");
}
// expose CPU interfaces
using WaveFunctionComponent::evaluateDerivatives;
using WaveFunctionComponent::evaluateGL;

7
src/QMCWaveFunctions/Fermion/DiracDeterminantBatched.cpp
View File

@ -16,6 +16,7 @@
#include "CPU/BLAS.hpp"
#include "OhmmsPETE/OhmmsMatrix.h"
#include "Numerics/MatrixOperators.h"
#include "QMCWaveFunctions/TWFFastDerivWrapper.h"
#include "CPU/SIMD/simd.hpp"
#include <cassert>
@ -956,6 +957,12 @@ void DiracDeterminantBatched<DET_ENGINE>::evaluateDerivatives(ParticleSet& P,
Phi->evaluateDerivatives(P, active, dlogpsi, dhpsioverpsi, FirstIndex, LastIndex);
}
template<typename DET_ENGINE>
void DiracDeterminantBatched<DET_ENGINE>::registerTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const
{
twf.addGroup(P, P.getGroupID(FirstIndex), Phi.get());
}
template<typename DET_ENGINE>
std::unique_ptr<DiracDeterminantBase> DiracDeterminantBatched<DET_ENGINE>::makeCopy(std::shared_ptr<SPOSet>&& spo) const
{

4
src/QMCWaveFunctions/Fermion/DiracDeterminantBatched.h
View File

@ -29,6 +29,9 @@
namespace qmcplusplus
{
//forward declaration
class TWFFastDerivWrapper;
template<typename DET_ENGINE = MatrixUpdateOMPTarget<QMCTraits::ValueType, QMCTraits::QTFull::ValueType>>
class DiracDeterminantBatched : public DiracDeterminantBase
{
@ -209,6 +212,7 @@ public:
void releaseResource(ResourceCollection& collection,
const RefVectorWithLeader<WaveFunctionComponent>& wfc_list) const override;
void registerTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const override;
/** cloning function
* @param tqp target particleset
* @param spo spo set

8
src/QMCWaveFunctions/Fermion/SlaterDet.cpp
View File

@ -253,4 +253,12 @@ std::unique_ptr<WaveFunctionComponent> SlaterDet::makeClone(ParticleSet& tqp) co
return myclone;
}
void SlaterDet::registerTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const
{
for (int i = 0; i < Dets.size(); ++i)
{
Dets[i]->registerTWFFastDerivWrapper(P, twf);
}
}
} // namespace qmcplusplus

4
src/QMCWaveFunctions/Fermion/SlaterDet.h
View File

@ -18,6 +18,7 @@
#ifndef QMCPLUSPLUS_SLATERDETERMINANT_WITHBASE_H
#define QMCPLUSPLUS_SLATERDETERMINANT_WITHBASE_H
#include "QMCWaveFunctions/Fermion/DiracDeterminantBase.h"
#include <map>
namespace qmcplusplus
@ -28,6 +29,7 @@ namespace qmcplusplus
// then change SlaterDet to SlaterDet<false>
// and SlaterDeterminantWithBackflow to SlaterDet<true>
// and remove all virtuals and inline them
class TWFFastDerivWrapper;
class SlaterDet : public WaveFunctionComponent
{
@ -53,6 +55,8 @@ public:
void reportStatus(std::ostream& os) override;
void registerTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const override;
LogValueType evaluateLog(const ParticleSet& P,
ParticleSet::ParticleGradient_t& G,
ParticleSet::ParticleLaplacian_t& L) override;

266
src/QMCWaveFunctions/TWFFastDerivWrapper.cpp Normal file
View File

@ -0,0 +1,266 @@
//////////////////////////////////////////////////////////////////////////////////////
// This file is distributed under the University of Illinois/NCSA Open Source License.
// See LICENSE file in top directory for details.
//
// Copyright (c) 2021 QMCPACK developers.
//
// File developed by: Raymond Clay III, rclay@sandia.gov, Sandia National Laboratories
//
// File created by: Raymond Clay III, rclay@sandia.gov, Sandia National Laboratories
//////////////////////////////////////////////////////////////////////////////////////
#include "QMCWaveFunctions/TWFFastDerivWrapper.h"
#include "Numerics/DeterminantOperators.h"
#include <iostream>
namespace qmcplusplus
{
TWFFastDerivWrapper::IndexType TWFFastDerivWrapper::getTWFGroupIndex(const IndexType gid)
{
IndexType return_group_index(-1);
for (IndexType i = 0; i < groups_.size(); i++)
if (gid == groups_[i])
return_group_index=i;
assert(return_group_index != -1);
return return_group_index;
}
void TWFFastDerivWrapper::addGroup(const ParticleSet& P, const IndexType gid, SPOSet* spo)
{
if (std::find(groups_.begin(), groups_.end(), gid) == groups_.end())
{
groups_.push_back(gid);
spos_.push_back(spo);
}
}
void TWFFastDerivWrapper::getM(const ParticleSet& P, std::vector<ValueMatrix_t>& mvec)
{
IndexType ngroups = spos_.size();
for (IndexType i = 0; i < ngroups; i++)
{
const IndexType gid = groups_[i];
const IndexType first = P.first(i);
const IndexType last = P.last(i);
const IndexType nptcls = last - first;
const IndexType norbs = spos_[i]->getOrbitalSetSize();
GradMatrix_t tmpgmat;
ValueMatrix_t tmplmat;
tmpgmat.resize(nptcls, norbs);
tmplmat.resize(nptcls, norbs);
spos_[i]->evaluate_notranspose(P, first, last, mvec[i], tmpgmat, tmplmat);
}
}
void TWFFastDerivWrapper::getEGradELaplM(const ParticleSet& P,
std::vector<ValueMatrix_t>& mvec,
std::vector<GradMatrix_t>& gmat,
std::vector<ValueMatrix_t>& lmat)
{
IndexType ngroups = mvec.size();
for (IndexType i = 0; i < ngroups; i++)
{
const IndexType gid = groups_[i];
const IndexType first = P.first(i);
const IndexType last = P.last(i);
const IndexType nptcls = last - first;
const IndexType norbs = spos_[i]->getOrbitalSetSize();
spos_[i]->evaluate_notranspose(P, first, last, mvec[i], gmat[i], lmat[i]);
}
}
void TWFFastDerivWrapper::getIonGradM(const ParticleSet& P,
const ParticleSet& source,
const int iat,
std::vector<std::vector<ValueMatrix_t>>& dmvec)
{
IndexType ngroups = dmvec[0].size();
for (IndexType i = 0; i < ngroups; i++)
{
const IndexType gid = groups_[i];
const IndexType first = P.first(i);
const IndexType last = P.last(i);
const IndexType nptcls = last - first;
const IndexType norbs = spos_[i]->getOrbitalSetSize();
GradMatrix_t grad_phi;
grad_phi.resize(nptcls, norbs);
spos_[i]->evaluateGradSource(P, first, last, source, iat, grad_phi);
for (IndexType idim = 0; idim < OHMMS_DIM; idim++)
for (IndexType iptcl = 0; iptcl < nptcls; iptcl++)
for (IndexType iorb = 0; iorb < norbs; iorb++)
{
dmvec[idim][i][iptcl][iorb] += grad_phi[iptcl][iorb][idim];
}
}
}
void TWFFastDerivWrapper::getIonGradIonGradELaplM(const ParticleSet& P,
const ParticleSet& source,
int iat,
std::vector<std::vector<ValueMatrix_t>>& dmvec,
std::vector<std::vector<ValueMatrix_t>>& dlmat)
{
IndexType ngroups = dmvec[0].size();
for (IndexType i = 0; i < ngroups; i++)
{
const IndexType gid = groups_[i];
const IndexType first = P.first(i);
const IndexType last = P.last(i);
const IndexType nptcls = last - first;
const IndexType norbs = spos_[i]->getOrbitalSetSize();
GradMatrix_t grad_phi;
HessMatrix_t grad_grad_phi;
GradMatrix_t grad_lapl_phi;
grad_phi.resize(nptcls, norbs);
grad_grad_phi.resize(nptcls, norbs);
grad_lapl_phi.resize(nptcls, norbs);
spos_[i]->evaluateGradSource(P, first, last, source, iat, grad_phi, grad_grad_phi, grad_lapl_phi);
for (IndexType idim = 0; idim < OHMMS_DIM; idim++)
for (IndexType iptcl = 0; iptcl < nptcls; iptcl++)
for (IndexType iorb = 0; iorb < norbs; iorb++)
{
dmvec[idim][i][iptcl][iorb] += grad_phi[iptcl][iorb][idim];
dlmat[idim][i][iptcl][iorb] += grad_lapl_phi[iptcl][iorb][idim];
}
}
}
TWFFastDerivWrapper::ValueType TWFFastDerivWrapper::computeGSDerivative(const std::vector<ValueMatrix_t>& Minv,
const std::vector<ValueMatrix_t>& X,
const std::vector<ValueMatrix_t>& dM,
const std::vector<ValueMatrix_t>& dB)
{
IndexType nspecies = Minv.size();
ValueType dval = 0.0;
for (int id = 0; id < nspecies; id++)
{
int ptclnum = Minv[id].rows();
ValueType dval_id = 0.0;
for (int i = 0; i < ptclnum; i++)
for (int j = 0; j < ptclnum; j++)
{
//Tr[M^{-1} dB - X * dM ]
dval_id += Minv[id][i][j] * dB[id][j][i] - X[id][i][j] * dM[id][j][i];
}
dval += dval_id;
}
return dval;
}
void TWFFastDerivWrapper::invertMatrices(const std::vector<ValueMatrix_t>& M, std::vector<ValueMatrix_t>& Minv)
{
IndexType nspecies = M.size();
for (IndexType id = 0; id < nspecies; id++)
{
assert(M[id].cols() == M[id].rows());
Minv[id] = M[id];
invert_matrix(Minv[id]);
}
}
void TWFFastDerivWrapper::buildX(const std::vector<ValueMatrix_t>& Minv,
const std::vector<ValueMatrix_t>& B,
std::vector<ValueMatrix_t>& X)
{
IndexType nspecies = Minv.size();
for (IndexType id = 0; id < nspecies; id++)
{
int ptclnum = Minv[id].rows();
assert(Minv[id].rows()==Minv[id].cols());
ValueMatrix_t tmpmat;
X[id].resize(ptclnum, ptclnum);
tmpmat.resize(ptclnum, ptclnum);
//(B*A^-1)
for (int i = 0; i < ptclnum; i++)
for (int j = 0; j < ptclnum; j++)
for (int k = 0; k < ptclnum; k++)
{
tmpmat[i][j] += B[id][i][k] * Minv[id][k][j];
}
//A^{-1}*B*A^{-1}
for (int i = 0; i < ptclnum; i++)
for (int j = 0; j < ptclnum; j++)
for (int k = 0; k < ptclnum; k++)
{
X[id][i][j] += Minv[id][i][k] * tmpmat[k][j];
}
}
}
void TWFFastDerivWrapper::wipeMatrices(std::vector<ValueMatrix_t>& A)
{
for (IndexType id = 0; id < A.size(); id++)
{
A[id] = 0.0;
}
}
TWFFastDerivWrapper::ValueType TWFFastDerivWrapper::trAB(const std::vector<ValueMatrix_t>& A, const std::vector<ValueMatrix_t>& B)
{
IndexType nspecies = A.size();
assert(A.size() == B.size());
ValueType val = 0.0;
//Now to compute the kinetic energy
for (IndexType id = 0; id < nspecies; id++)
{
int ptclnum = A[id].rows();
ValueType val_id = 0.0;
assert(A[id].cols() == B[id].rows() && A[id].rows() == B[id].cols());
for (int i = 0; i < A[id].rows(); i++)
for (int j = 0; j < A[id].cols(); j++)
{
val_id += A[id][i][j] * B[id][j][i];
}
val += val_id;
}
return val;
}
void TWFFastDerivWrapper::getGSMatrices(const std::vector<ValueMatrix_t>& A, std::vector<ValueMatrix_t>& Aslice)
{
IndexType nspecies = A.size();
Aslice.resize(nspecies);
for (IndexType id = 0; id < nspecies; id++)
{
IndexType ptclnum = A[id].rows();
Aslice[id].resize(ptclnum, ptclnum);
for (IndexType i = 0; i < ptclnum; i++)
for (IndexType j = 0; j < ptclnum; j++)
Aslice[id][i][j] = A[id][i][j];
}
}
TWFFastDerivWrapper::IndexType TWFFastDerivWrapper::getRowM(const ParticleSet& P, const IndexType iel, ValueVector_t& val)
{
IndexType gid = P.getGroupID(iel);
IndexType sid = getTWFGroupIndex(gid);
GradVector_t tempg;
ValueVector_t templ;
IndexType norbs = spos_[sid]->getOrbitalSetSize();
tempg.resize(norbs);
templ.resize(norbs);
spos_[sid]->evaluateVGL(P, iel, val, tempg, templ);
return sid;
}
} // namespace qmcplusplus

229
src/QMCWaveFunctions/TWFFastDerivWrapper.h Normal file
View File

@ -0,0 +1,229 @@
//////////////////////////////////////////////////////////////////////////////////////
// This file is distributed under the University of Illinois/NCSA Open Source License.
// See LICENSE file in top directory for details.
//
// Copyright (c) 2021 QMCPACK developers.
//
// File developed by: Raymond Clay III, rclay@sandia.gov, Sandia National Laboratories
//
// File created by: Raymond Clay III, rclay@sandia.gov, Sandia National Laboratories
//////////////////////////////////////////////////////////////////////////////////////
#ifndef QMCPLUSPLUS_TWFFASTDERIVWRAPPER_H
#define QMCPLUSPLUS_TWFFASTDERIVWRAPPER_H
#include "QMCWaveFunctions/WaveFunctionComponent.h"
#include "QMCWaveFunctions/SPOSet.h"
#include "Configuration.h"
#include "Particle/ParticleSet.h"
namespace qmcplusplus
{
/**
* TWFFastDerivWrapper is a wrapper class for TrialWavefunction that provides separate and low level access to the Jastrow and
* SPOSet objects. This is so that observables can be recast in matrix form and their derivatives taken efficiently.
* Currently this is hard coded for ground state slater jastrow wave functions, but generalization to
* arbitrary occupations and multideterminants are straightforward and will come online as the code is tested and validated.
*
* Please see : J. Chem. Phys. 144, 194105 (2016) https://doi.org/10.1063/1.4948778 for implementation details and formalism.
*/
class TWFFastDerivWrapper
{
public:
using ValueMatrix_t = SPOSet::ValueMatrix_t;
using GradMatrix_t = SPOSet::GradMatrix_t;
using HessMatrix_t = SPOSet::HessMatrix_t;
using IndexType = QMCTraits::IndexType;
using RealType = QMCTraits::RealType;
using ValueType = QMCTraits::ValueType;
using ValueVector_t = SPOSet::ValueVector_t;
using GradVector_t = SPOSet::GradVector_t;
/** @brief Add a particle group.
*
* Here, a "group" corresponds to a subset of particles which are antisymmetric with
* respect to eachother. TWFFastDerivWrapper ensures that there is a binding between the groupid
* in ParticleSet and the SPOSet associated with that particle group. This function stores
* the ParticleSet groupid and SPOSet in a vector for lookup and communication with QMCPACK conventions,
* but is agnostic to the order of group registration or evaluation.
*
* @param[in] P. ParticleSet
* @param[in] groupid. ParticleSet groupid to which SPOSet corresponds.
* @param[in] spo. Pointer to SPOSet.
* @return void.
*/
void addGroup(const ParticleSet& P, const IndexType groupid, SPOSet* spo);
inline void addJastrow(WaveFunctionComponent* j) { jastrow_list_.push_back(j); };
/** @brief Takes particle set groupID and returns the TWF internal index for it.
*
* ParticleSet groups can be registered in whichever order. However, the internal indexing
* of TWFFastDerivWrapper keeps track on a first-come, first serve basis. That is, if I register
* particle groups 3, 1, and 0 in that order, then the internal indexing goes like
* 0->3, 1->1, 2->0. Thus, this function looks up where in the internal indexing scheme
* ParticleSet gid is located. This is necessary, since the matrix list data structures follow
* the internal TWF indexing.
*
* @param[in] gid. ParticleSet group ID to look up.
* @return The corresponding internal groupID.
*/
IndexType getTWFGroupIndex(const IndexType gid);
/** @ingroup Query functions
* @{
* @brief These are query functions to the internal state of various groups and SPOSet info.
* All group indices will refer to the internal group indices. Convert from ParticleSet to
* TWF group first!
*
* Source of truth for orbital sizes will be the individual SPOSets. Particle group sizes
* will be ParticleSet in conjunction with groupID maps.
*/
inline IndexType numGroups() { return spos_.size(); };
SPOSet* getSPOSet(const IndexType sid) { return spos_[sid]; };
inline IndexType numOrbitals(const IndexType sid) { return spos_[sid]->size(); };
/** @} */
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//These are convenience functions/wrappers to SPOSet calls. Idea being that observables just need //
//to make calls to this object to build the auxiliary matrices required for fast derivative computation.//
//On the other hand, it wouldn't be unreasonable to make the observables do the SPOSet calls themselves.//
//////////////////////////////////////////////////////////////////////////////////////////////////////////
/** @brief Returns the non-rectangular slater matrices M_ij=phi_j(r_i) (N_particle x N_orbs) for each species group.
*
* @param[in] P particle set.
* @param[in,out] mmat Output vector of slater matrices. Each vector entry corresponds to a different particle group.
* @return Void
*/
void getM(const ParticleSet& P, std::vector<ValueMatrix_t>& mmat);
/** @brief Returns value of all orbitals (relevant to given species/group) at a particular particle coordinate.
*
* @param[in] P particle set.
* @param[in] iel particle ID.
* @param[in,out] val Vector of phi_i(r_iel) for all i=0,Norbs.
* @return Void
*/
IndexType getRowM(const ParticleSet& P, const IndexType iel, ValueVector_t& val);
/** @brief Returns value, gradient, and laplacian matrices for all orbitals and all particles, species by species.
*
* @param[in] P particle set.
* @param[in,out] mvec Slater matrix M_ij=phi_j(r_i) for each species group.
* @param[in,out] gmat electron gradient of slater matrix [G_ij]_a = d/dr_a,i phi_j(r_i). a=x,y,z
* @param[in,out] lmat electron laplacian of slater matrix [L_ij] = \nabla^2_i phi_j(r_i).
* @return Void
*/
void getEGradELaplM(const ParticleSet& P,
std::vector<ValueMatrix_t>& mvec,
std::vector<GradMatrix_t>& gmat,
std::vector<ValueMatrix_t>& lmat);
/** @brief Returns x,y,z components of ion gradient of slater matrices.
*
* @param[in] P particle set.
* @param[in] source ion particle set.
* @param[in]*** iat ion ID w.r.t. which to take derivative.
* @param[in,out] dmvec Slater matrix d/dR_{iat,a} M_ij=d/dR_{iat,a} phi_j(r_i) for each species group. First index is a=x,y,z.
* @return Void
*/
void getIonGradM(const ParticleSet& P,
const ParticleSet& source,
const int iat,
std::vector<std::vector<ValueMatrix_t>>& dmvec);
/** @brief Returns x,y,z components of ion gradient of slater matrices and their laplacians..
*
* @param[in] P particle set.
* @param[in] source ion particle set.
* @param[in] iat ion ID w.r.t. which to take derivative.
* @param[in,out] dmvec Slater matrix d/dR_{iat,a} M_ij=d/dR_{iat,a} phi_j(r_i) for each species group.
* First index is a=x,y,z.
* @param[in,out] dlmat Slater matrix d/dR_{iat,a} L_ij=d/dR_{iat,a} \nabla^2_i phi_j(r_i) for each species group.
* First index is a=x,y,z.
* @return Void
*/
void getIonGradIonGradELaplM(const ParticleSet& P,
const ParticleSet& source,
int iat,
std::vector<std::vector<ValueMatrix_t>>& dmvec,
std::vector<std::vector<ValueMatrix_t>>& dlmat);
/** @brief Takes sub matrices of full SPOSet quantities (evaluated on all particles and all orbitals), consistent with ground
* state occupations.
*
* @param[in] A non-rectangular matrices of SPOSet derived quantities, evaluated on all orbitals and all particles.
* @param[in,out] Aslice square matrices consistent with a ground state occupation.
* @return Void
*/
void getGSMatrices(const std::vector<ValueMatrix_t>& A, std::vector<ValueMatrix_t>& Aslice);
/** @brief Calculates derivative of observable via Tr[M^{-1} dB - X * dM ]. Consistent with ground state occupation.
*
* @param[in] Minv. inverse of slater matrices for ground state occupations.
* @param[in] X. X=M^-1 B M^-1. B observable matrix, and is consistent with ground state occupation.
* @param[in] dM. Target derivative of M, and is consistent with ground state occupation.
* @param[in] dB. Target derivative of B, and is consistent with ground state occupation.
* @return Derivative of O psi/psi = Tr[M^{-1} dB - X * dM ]
*/
ValueType computeGSDerivative(const std::vector<ValueMatrix_t>& Minv,
const std::vector<ValueMatrix_t>& X,
const std::vector<ValueMatrix_t>& dM,
const std::vector<ValueMatrix_t>& dB);
//////////////////////////////////////////////////////////////////////////////////////////////
//And now we just have some helper functions for doing useful math on our lists of matrices.//
//////////////////////////////////////////////////////////////////////////////////////////////
/** @brief Helper function that inverts all slater matrices in our species list.
*
* @param[in] M. List of slater matrices for each species. These are square and consistent with some occupation.
* @param[in,out] Minv. The species by species list of inverted matrices from M.
* @return Void.
*/
void invertMatrices(const std::vector<ValueMatrix_t>& M, std::vector<ValueMatrix_t>& Minv);
/** @brief Build the auxiliary X=M^-1 B M^-1 matrix.
*
* @param[in] Minv. List of slater matrix inverses M^-1 for a given occupation.
* @param[in] B. Observable auxiliary matrix for a given occupation.
* @param[in,out] X. M^-1*B*M^-1 is stored in this list of matrices.
* @return Void.
*/
void buildX(const std::vector<ValueMatrix_t>& Minv,
const std::vector<ValueMatrix_t>& B,
std::vector<ValueMatrix_t>& X);
/** @brief Goes through a list of matrices and zeros them out. Does this in place.
*
* @param[in,out] A. The list of matrices to be zeroed out. After call, A is all zeros.
* @return Void.
*/
void wipeMatrices(std::vector<ValueMatrix_t>& A);
/** @brief Returns Tr(A*B). Actually, we assume a block diagonal structure, so this is
* really Sum_i Tr(A_i*B_i), where i is the species index.
*
* @param[in] A. Vector of matrices A.
* @param[in] B. Vector of matrices B.
* @return Value of Sum_i Tr(A_i*B_i).
*/
ValueType trAB(const std::vector<ValueMatrix_t>& A, const std::vector<ValueMatrix_t>& B);
private:
std::vector<SPOSet*> spos_;
std::vector<IndexType> groups_;
std::vector<ValueMatrix_t> psi_M_;
std::vector<ValueMatrix_t> psi_M_inv_;
std::vector<WaveFunctionComponent*> jastrow_list_;
};
/**@}*/
} // namespace qmcplusplus
#endif

16
src/QMCWaveFunctions/TrialWaveFunction.cpp
View File

@ -1340,5 +1340,21 @@ RefVector<ParticleSet::ParticleLaplacian_t> TrialWaveFunction::extractLRefList(
return l_list;
}
void TrialWaveFunction::initializeTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const
{
for (int i = 0; i < Z.size(); ++i)
{
if (Z[i]->is_fermionic)
{
//OK, so this is a hack only for SlaterDeterminant objects.
//Needs a bit of logic and protection before this reaches production.
//SlaterDet* det = dynamic_cast<SlaterDet*>(Z[i].get());
//det->registerTWFFastDerivWrapper(P, twf);
Z[i]->registerTWFFastDerivWrapper(P,twf);
}
else //Is Jastrow, so do nothing right now.
{}
}
}
} // namespace qmcplusplus

6
src/QMCWaveFunctions/TrialWaveFunction.h
View File

@ -29,6 +29,7 @@
#include "Utilities/TimerManager.h"
#include "type_traits/template_types.hpp"
#include "Containers/MinimalContainers/RecordArray.hpp"
#include "QMCWaveFunctions/TWFFastDerivWrapper.h"
#ifdef QMC_CUDA
#include "type_traits/CUDATypes.h"
#endif
@ -168,6 +169,9 @@ public:
*/
void reportStatus(std::ostream& os);
/** Initialize a TWF wrapper for fast derivative evaluation
*/
void initializeTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const;
/** evalaute the log (internally gradients and laplacian) of the trial wavefunction. gold reference */
RealType evaluateLog(ParticleSet& P);
@ -539,6 +543,8 @@ private:
///a list of WaveFunctionComponents constituting many-body wave functions
std::vector<std::unique_ptr<WaveFunctionComponent>> Z;
/// For now, TrialWaveFunction will own the wrapper.
TWFFastDerivWrapper twf_prototype;
/// timers at TrialWaveFunction function call level
TimerList_t TWF_timers_;
/// timers at WaveFunctionComponent function call level

7
src/QMCWaveFunctions/WaveFunctionComponent.cpp
View File

@ -226,4 +226,11 @@ void WaveFunctionComponent::evaluateDerivRatios(VirtualParticleSet& VP,
evaluateRatios(VP, ratios);
}
void WaveFunctionComponent::registerTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const
{
std::ostringstream o;
o << "WaveFunctionComponent::registerTWFFastDerivWrapper is not implemented by " << ClassName;
APP_ABORT(o.str());
}
} // namespace qmcplusplus

6
src/QMCWaveFunctions/WaveFunctionComponent.h
View File

@ -50,7 +50,7 @@ struct NLjob
class WaveFunctionComponent;
struct DiffWaveFunctionComponent;
class ResourceCollection;
class TWFFastDerivWrapper;
/**@defgroup WaveFunctionComponent group
* @brief Classes which constitute a many-body trial wave function
*
@ -170,6 +170,10 @@ public:
/** print the state, e.g., optimizables */
virtual void reportStatus(std::ostream& os) = 0;
/** Register the component with the TWFFastDerivWrapper wrapper.
*/
virtual void registerTWFFastDerivWrapper(const ParticleSet& P, TWFFastDerivWrapper& twf) const;
/** evaluate the value of the WaveFunctionComponent from scratch
* \param[in] P active ParticleSet
* \param[out] G Gradients, \f$\nabla\ln\Psi\f$