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

Merge branch 'develop' into rmg-conv-fix

This commit is contained in:
Ye Luo 2022-08-23 21:33:24 -05:00 committed by GitHub
parent c47bea5be3 4d9e4b481d
commit e88319f564
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
275 changed files with 79575 additions and 22111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
.github/workflows/ci-github-actions.yaml vendored
View File

@ -30,6 +30,8 @@ jobs:
GCC9-MPI-Gcov-Complex,
GCC11-NoMPI-Werror-Real,
GCC11-NoMPI-Werror-Complex,
GCC11-NoMPI-Werror-Real-Mixed,
GCC11-NoMPI-Werror-Complex-Mixed,
Clang10-NoMPI-ASan-Real,
Clang10-NoMPI-ASan-Complex,
Clang10-NoMPI-UBSan-Real,
@ -76,6 +78,16 @@ jobs:
image: williamfgc/qmcpack-ci:ubuntu2110-serial
options: -u 1001
- jobname: GCC11-NoMPI-Werror-Real-Mixed
container:
image: williamfgc/qmcpack-ci:ubuntu2110-serial
options: -u 1001
- jobname: GCC11-NoMPI-Werror-Complex-Mixed
container:
image: williamfgc/qmcpack-ci:ubuntu2110-serial
options: -u 1001
- jobname: Clang10-NoMPI-ASan-Real
container:
image: williamfgc/qmcpack-ci:ubuntu20-openmpi

4
CMakeLists.txt
View File

@ -18,6 +18,10 @@ project(
VERSION 3.14.9
LANGUAGES C CXX)
# add the automatically determined parts of the RPATH
# which point to directories outside the build tree to the install RPATH
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
#--------------------------------------------------------------------
# Directory where customize cmake files reside
#--------------------------------------------------------------------

95
config/build_alcf_polaris_Clang.sh Executable file
View File

@ -0,0 +1,95 @@
#!/bin/bash
# This recipe is intended for ALCF Polaris https://www.alcf.anl.gov/polaris
# It builds all the varaints of QMCPACK in the current directory
# last revision: Aug 17th 2022
#
# How to invoke this script?
# build_alcf_polaris_Clang.sh # build all the variants assuming the current directory is the source directory.
# build_alcf_polaris_Clang.sh <source_dir> # build all the variants with a given source directory <source_dir>
# build_alcf_polaris_Clang.sh <source_dir> <install_dir> # build all the variants with a given source directory <source_dir> and install to <install_dir>
module load mpiwrappers/cray-mpich-llvm llvm/main-20220317
module load cudatoolkit-standalone/11.2.2
module load cray-fftw/3.3.8.13
module load cray-hdf5-parallel/1.12.1.3
module load cmake/3.23.2
export BOOST_ROOT=/soft/applications/qmcpack/boost_1_79_0
export CMAKE_PREFIX_PATH=/soft/libraries/openblas/0.3.20-omp:$CMAKE_PREFIX_PATH
echo "**********************************"
echo '$ clang -v'
clang -v
echo "**********************************"
TYPE=Release
Machine=polaris
Compiler=Clang
if [[ $# -eq 0 ]]; then
source_folder=`pwd`
elif [[ $# -eq 1 ]]; then
source_folder=$1
else
source_folder=$1
install_folder=$2
fi
if [[ -f $source_folder/CMakeLists.txt ]]; then
echo Using QMCPACK source directory $source_folder
else
echo "Source directory $source_folder doesn't contain CMakeLists.txt. Pass QMCPACK source directory as the first argument."
exit
fi
for name in offload_cuda_real_MP offload_cuda_real offload_cuda_cplx_MP offload_cuda_cplx \
cpu_real_MP cpu_real cpu_cplx_MP cpu_cplx
do
CMAKE_FLAGS="-DCMAKE_BUILD_TYPE=$TYPE"
if [[ $name == *"cplx"* ]]; then
CMAKE_FLAGS="$CMAKE_FLAGS -DQMC_COMPLEX=ON"
fi
if [[ $name == *"_MP"* ]]; then
CMAKE_FLAGS="$CMAKE_FLAGS -DQMC_MIXED_PRECISION=ON"
fi
if [[ $name == *"offload"* ]]; then
CMAKE_FLAGS="$CMAKE_FLAGS -DENABLE_OFFLOAD=ON -DUSE_OBJECT_TARGET=ON -DOFFLOAD_ARCH=sm_80"
fi
if [[ $name == *"cuda"* ]]; then
CMAKE_FLAGS="$CMAKE_FLAGS -DENABLE_CUDA=ON -DCMAKE_CUDA_ARCHITECTURES=80"
fi
folder=build_${Machine}_${Compiler}_${name}
if [[ -v install_folder ]]; then
CMAKE_FLAGS="$CMAKE_FLAGS -DCMAKE_INSTALL_PREFIX=$install_folder/$folder"
fi
echo "**********************************"
echo "$folder"
echo "$CMAKE_FLAGS"
echo "**********************************"
mkdir $folder
cd $folder
if [ ! -f CMakeCache.txt ] ; then
cmake $CMAKE_FLAGS -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx $source_folder
fi
if [[ -v install_folder ]]; then
make -j16 install && chmod -R -w $install_folder/$folder
else
make -j16
fi
cd ..
echo
done

2
docs/LCAO.rst
View File

@ -161,7 +161,7 @@ To generate QMCPACK input files, you will need to run ``convert4qmc`` exactly as
::
convert4qmc -pyscf C_diamond-tiled-cplx
convert4qmc -orbitals C_diamond-tiled-cplx
This tool can be used with any option described in convert4qmc. Since
the HDF5 contains all the information needed, there is no need to

6
docs/installation.rst
View File

@ -882,6 +882,12 @@ Each node features a second-generation Intel Xeon Phi 7230 processor and 192 GB
make -j 24
ls -l bin/qmcpack
Installing on ALCF Polaris
~~~~~~~~~~~~~~~~~~~~~~~~~~
Polaris is a HPE Apollo Gen10+ based 44 petaflops system.
Each node features a AMD EPYC 7543P CPU and 4 NVIDIA A100 GPUs.
A build recipe for Polaris can be found at ``<qmcpack_source>/config/build_alcf_polaris_Clang.sh``
Installing on ORNL OLCF Summit
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

3
docs/lab_advanced_molecules.rst
View File

@ -880,7 +880,6 @@ Appendix C: Wavefunction optimization XML block
<parameter name="bigchange">10.0</parameter>
<estimator name="LocalEnergy" hdf5="no"/>
<parameter name="usebuffer"> yes </parameter>
<parameter name="nonlocalpp"> yes </parameter>
<parameter name="MinMethod">quartic</parameter>
<parameter name="exp0">-6</parameter>
<parameter name="alloweddifference"> 1.0e-5 </parameter>
@ -895,8 +894,6 @@ Options:
- usebuffer: (default no) Save useful information during VMC
- nonlocalpp: (default no) Include nonlocal energy on 1-D min
- MinMethod: (default quartic) Method to calculate magnitude of
parameter change quartic: fit quartic polynomial to four values of
the cost function obtained using reweighting along chosen direction

11
docs/lab_qmc_basics.rst
View File

@ -347,8 +347,6 @@ The relevant portion of the input describing the linear optimization process is
<parameter name="warmupSteps" > 50 </parameter>
<parameter name="blocks" > 200 </parameter>
<parameter name="subSteps" > 1 </parameter>
<parameter name="nonlocalpp" > yes </parameter>
<parameter name="useBuffer" > yes </parameter>
...
</qmc>
</loop>
@ -397,13 +395,6 @@ subSteps
measurements can be more efficient. Will be less efficient with many
substeps.
nonlocalpp,useBuffer
If ``nonlocalpp="no,"`` then the nonlocal part of the pseudopotential
is not included when computing the cost function. If
``useBuffer="yes,"`` then temporary data is stored to speed up
nonlocal pseudopotential evaluation at the expense of memory
consumption.
loop max
Number of times to repeat the optimization. Using the resulting
wavefunction from the previous optimization in the next one improves
@ -962,8 +953,6 @@ Objects representing QMCPACK simulations are then constructed with the ``generat
warmupsteps = 50,
blocks = 200,
substeps = 1,
nonlocalpp = True,
usebuffer = True,
walkers = 1,
minwalkers = 0.5,
maxweight = 1e9,

37
docs/methods.rst
View File

@ -519,9 +519,9 @@ The input parameters are listed in the following table.
+--------------------------+--------------+-------------+-------------+--------------------------------------------------+
| **Name** | **Datatype** | **Values** | **Default** | **Description** |
+==========================+==============+=============+=============+==================================================+
| ``nonlocalpp`` | text | yes, no | no | include non-local PP energy in the cost function |
| ``nonlocalpp`` | text | | | No more effective. Will be removed. |
+--------------------------+--------------+-------------+-------------+--------------------------------------------------+
| ``use_nonlocalpp_deriv`` | text | yes, no | yes | Add non-local PP energy derivative contribution |
| ``use_nonlocalpp_deriv`` | text | | | No more effective. Will be removed. |
+--------------------------+--------------+-------------+-------------+--------------------------------------------------+
| ``minwalkers`` | real | 0--1 | 0.3 | Lower bound of the effective weight |
+--------------------------+--------------+-------------+-------------+--------------------------------------------------+
@ -532,13 +532,7 @@ Additional information:
- ``maxWeight`` The default should be good.
- ``nonlocalpp`` The ``nonlocalpp`` contribution to the local energy depends on the
wavefunction. When a new set of parameters is proposed, this
contribution needs to be updated if the cost function consists of local
energy. Fortunately, nonlocal contribution is chosen small when making a
PP for small locality error. We can ignore its change and avoid the
expensive computational cost. An implementation issue with legacy GPU code is
that a large amount of memory is consumed with this option.
- ``nonlocalpp`` and ``use_nonlocalpp_deriv`` are obsolete and will be treated as invalid options (trigger application abort) in future releases. From this point forward, the code behaves as prior versions of qmcpack did when both were set to ``yes``.
- ``minwalkers`` This is a ``critical`` parameter. When the ratio of effective samples to actual number of samples in a reweighting step goes lower than ``minwalkers``,
the proposed set of parameters is invalid.
@ -551,6 +545,31 @@ The cost function consists of three components: energy, unreweighted variance, a
<cost name="unreweightedvariance"> 0.00 </cost>
<cost name="reweightedvariance"> 0.05 </cost>
Varational parameter selection
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The predominant way of selecting varational paramemters is via ``<wavefunction>`` input.
``<coefficients>`` entries support ``optimize="yes"/"no"`` to enable/disable variational parameters in the wavefunction optimization.
The secondary way of selecting varational paramemters is via ``variational_subset`` parameter in the ``<qmc>`` driver input.
It allows controlling optimization granularity at each optimization step.
If ``variational_subset`` is not provided or empty, all the varational paramemters are selected.
If variational paramemters are set as not optimizable in the predominant way, the secondary way won't be able to set them optimizable even they are selected.
The following example shows optimizing subsets of parameters in stages in a single QMCPACK run.
::
<qmc method="linear">
...
<parameter name="variational_subset"> uu ud </parameter>
</qmc>
<qmc method="linear">
...
<parameter name="variational_subset"> uu ud eH </parameter>
</qmc>
<qmc method="linear">
...
<parameter name="variational_subset"> uu ud eH CI </parameter>
</qmc>
Optimizers
~~~~~~~~~~

0
docs/running_docker.rst Executable file → Normal file
View File

1
examples/molecules/He/he_bspline_jastrow.xml
View File

@ -86,7 +86,6 @@
<parameter name="useDrift"> yes </parameter>
<parameter name="bigchange">10.0</parameter>
<estimator name="LocalEnergy" hdf5="no"/>
<parameter name="nonlocalpp"> yes </parameter>
<parameter name="MinMethod">quartic</parameter>
<parameter name="exp0">-6</parameter>
<parameter name="alloweddifference"> 1.0e-5 </parameter>

1
examples/molecules/He/he_example_wf.xml
View File

@ -63,7 +63,6 @@
<parameter name="useDrift"> yes </parameter>
<parameter name="bigchange">10.0</parameter>
<estimator name="LocalEnergy" hdf5="no"/>
<parameter name="nonlocalpp"> yes </parameter>
<parameter name="MinMethod">quartic</parameter>
<parameter name="exp0">-6</parameter>
<parameter name="alloweddifference"> 1.0e-5 </parameter>

1
examples/molecules/He/he_from_gamess.xml
View File

@ -2076,7 +2076,6 @@
<parameter name="useDrift"> yes </parameter>
<parameter name="bigchange">10.0</parameter>
<estimator name="LocalEnergy" hdf5="no"/>
<parameter name="nonlocalpp"> yes </parameter>
<parameter name="MinMethod">quartic</parameter>
<parameter name="exp0">-6</parameter>
<parameter name="alloweddifference"> 1.0e-5 </parameter>

1
examples/molecules/He/he_simple_opt.xml
View File

@ -83,7 +83,6 @@
<parameter name="useDrift"> yes </parameter>
<parameter name="bigchange">10.0</parameter>
<estimator name="LocalEnergy" hdf5="no"/>
<parameter name="nonlocalpp"> yes </parameter>
<parameter name="MinMethod">quartic</parameter>
<parameter name="exp0">-6</parameter>
<parameter name="alloweddifference"> 1.0e-5 </parameter>

32
src/Containers/OhmmsPETE/ompReductionTinyVector.hpp Normal file
View File

@ -0,0 +1,32 @@
//////////////////////////////////////////////////////////////////////////////////////
// This file is distributed under the University of Illinois/NCSA Open Source License.
// See LICENSE file in top directory for details.
//
// Copyright (c) 2022 QMCPACK developers.
//
// File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory
//
// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory
//////////////////////////////////////////////////////////////////////////////////////
#ifndef QMCPLUSPLUS_OMPREDUCTION_TINYVECTOR_H
#define QMCPLUSPLUS_OMPREDUCTION_TINYVECTOR_H
#include <complex>
#include "config.h"
#include "TinyVector.h"
namespace qmcplusplus
{
#if !defined(OPENMP_NO_UDR)
PRAGMA_OFFLOAD("omp declare reduction(+: TinyVector<float, OHMMS_DIM>: omp_out += omp_in)")
PRAGMA_OFFLOAD("omp declare reduction(+: TinyVector<double, OHMMS_DIM>: omp_out += omp_in)")
#endif
#if !defined(OPENMP_NO_COMPLEX) && !defined(OPENMP_NO_UDR)
PRAGMA_OFFLOAD("omp declare reduction(+: TinyVector<std::complex<float>, OHMMS_DIM>: omp_out += omp_in)")
PRAGMA_OFFLOAD("omp declare reduction(+: TinyVector<std::complex<double>, OHMMS_DIM>: omp_out += omp_in)")
#endif
}
#endif // QMCPLUSPLUS_OMPREDUCTION_TINYVECTOR_H

1
src/Estimators/EstimatorManagerInput.cpp
View File

@ -15,6 +15,7 @@
#include "MomentumDistributionInput.h"
#include "OneBodyDensityMatricesInput.h"
#include "SpinDensityInput.h"
#include "ModernStringUtils.hpp"
namespace qmcplusplus
{

5
src/Estimators/InputSection.cpp
View File

@ -11,7 +11,9 @@
#include "InputSection.h"
#include "Message/UniformCommunicateError.h"
#include "ModernStringUtils.hpp"
#include "Utilities/string_utils.h"
namespace qmcplusplus
{
@ -51,8 +53,7 @@ void InputSection::readXML(xmlNodePtr cur)
std::istringstream stream(XMLNodeString{element});
setFromStream(name, stream);
}
else if (isDelegate(ename))
{}
else if (isDelegate(ename)) {}
element = element->next;
}

1
src/Estimators/ScalarEstimatorInputs.cpp
View File

@ -12,6 +12,7 @@
#include "ScalarEstimatorInputs.h"
#include "InputSection.h"
#include "EstimatorInput.h"
#include "ModernStringUtils.hpp"
namespace qmcplusplus
{

1
src/Estimators/tests/test_ScalarEstimatorInputs.cpp
View File

@ -15,6 +15,7 @@
#include "OhmmsData/Libxml2Doc.h"
#include "ScalarEstimatorInputs.h"
#include "ValidScalarEstimatorInput.h"
#include "ModernStringUtils.hpp"
namespace qmcplusplus
{

4
src/Numerics/GaussianTimesRN.h
View File

@ -122,9 +122,9 @@ struct GaussianTimesRN : public OptimizableFunctorBase
void reset() override;
void checkInVariables(opt_variables_type& active) override {}
void checkInVariablesExclusive(opt_variables_type& active) override {}
void checkOutVariables(const opt_variables_type& active) override {}
void resetParameters(const opt_variables_type& active) override
void resetParametersExclusive(const opt_variables_type& active) override
{
///DO NOTHING FOR NOW
}

105
src/Numerics/MatrixOperators.h
View File

@ -26,15 +26,6 @@
namespace qmcplusplus
{
template<typename T>
inline T TESTDOT(const T* restrict f, const T* restrict l, const T* restrict b)
{
T s = 0;
while (f != l)
s += (*f++) * (*b++);
return s;
}
namespace MatrixOperators
{
/** static function to perform C=AB for real matrices
@ -173,45 +164,23 @@ inline void product(const Matrix<double>& A, const Matrix<std::complex<double>>&
/** static function to perform y=Ax for generic matrix and vector
*/
template<typename T>
inline void product(const Matrix<T>& A, const Vector<T>& x, T* restrict yptr)
inline void product(const Matrix<T>& A, const Vector<T>& x, Vector<T>& y)
{
const char transa = 'T';
const T one = 1.0;
const T zero = 0.0;
BLAS::gemv(transa, A.cols(), A.rows(), one, A.data(), A.cols(), x.data(), 1, zero, yptr, 1);
BLAS::gemv(transa, A.cols(), A.rows(), one, A.data(), A.cols(), x.data(), 1, zero, y.data(), 1);
}
/** static function to perform y = A^t x for generic matrix and vector
*/
template<typename T>
inline void product_Atx(const Matrix<T>& A, const Vector<T>& x, T* restrict yptr)
inline void product_Atx(const Matrix<T>& A, const Vector<T>& x, Vector<T>& y)
{
const char transa = 'N';
const T one = 1.0;
const T zero = 0.0;
BLAS::gemv(transa, A.cols(), A.rows(), one, A.data(), A.cols(), x.data(), 1, zero, yptr, 1);
}
/** static function to perform y=Ax for generic matrix and vector
*/
template<typename T>
inline void product(const Matrix<T>& A, const T* restrict xptr, T* restrict yptr)
{
const char transa = 'T';
const T one = 1.0;
const T zero = 0.0;
BLAS::gemv(transa, A.cols(), A.rows(), one, A.data(), A.cols(), xptr, 1, zero, yptr, 1);
}
/** static function to perform y = A^t x for generic matrix and vector
*/
template<typename T>
inline void product_Atx(const Matrix<T>& A, const T* restrict xptr, T* restrict yptr)
{
const char transa = 'N';
const T one = 1.0;
const T zero = 0.0;
BLAS::gemv(transa, A.cols(), A.rows(), one, A.data(), A.cols(), xptr, 1, zero, yptr, 1);
BLAS::gemv(transa, A.cols(), A.rows(), one, A.data(), A.cols(), x.data(), 1, zero, y.data(), 1);
}
/** static function to perform y=Ax for generic matrix and vector
@ -294,73 +263,7 @@ inline void product(const Matrix<std::complex<double>>& A, const Vector<double>&
}
}
inline void product(const Matrix<std::complex<double>>& A,
const std::complex<double>* restrict x,
std::complex<double>* restrict yptr)
{
const char transa = 'T';
const std::complex<double> zone(1.0, 0.0);
const std::complex<double> zero(0.0, 0.0);
zgemv(transa, A.cols(), A.rows(), zone, A.data(), A.cols(), x, 1, zero, yptr, 1);
}
/** static function to perform y=Ax for generic matrix and vector
*/
inline void product(const Matrix<double>& A, const Vector<std::complex<double>>& x, double* restrict yptr)
{
std::cerr << " Undefined C=AB with real A and complex x " << std::endl;
}
template<typename T>
inline void product(const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& C, std::vector<int>& offset)
{
int nr = C.rows();
int nb = offset.size() - 1;
for (int i = 0; i < nr; i++)
{
for (int b = 0; b < nb; b++)
{
int firstK = offset[b];
int lastK = offset[b + 1];
const T* restrict firstY = A[i] + firstK;
const T* restrict lastY = A[i] + lastK;
for (int k = firstK; k < lastK; k++)
{
C(i, k) = TESTDOT(firstY, lastY, B[k] + firstK);
}
}
}
}
// template<typename T>
// inline statis void product(const Matrix<T>& A, const T* restrict x, T* restrict y)
// {
// GEMV<T,0>::apply(A.data(),x,y,A.rows(),A.cols());
// }
} // namespace MatrixOperators
/** API to handle gemv */
namespace simd
{
template<typename T>
inline void gemv(const Matrix<T>& a, const T* restrict v, T* restrict b)
{
MatrixOperators::product(a, v, b);
}
template<typename T, unsigned D>
inline void gemv(const Matrix<T>& a, const TinyVector<T, D>* restrict v, TinyVector<T, D>* restrict b)
{
MatrixOperators::product(a, v, b);
}
template<typename T, unsigned D>
inline void gemv(const Matrix<T>& a, const Tensor<T, D>* restrict v, Tensor<T, D>* restrict b)
{
MatrixOperators::product(a, v, b);
}
} // namespace simd
} // namespace qmcplusplus
#endif

21
src/Particle/ParticleIO/XMLParticleIO.cpp
View File

@ -229,7 +229,8 @@ bool XMLParticleParser::readXML(xmlNodePtr cur)
throw UniformCommunicateError("'ionid' only supports datatype=\"" + stringtype_tag + "\"");
std::vector<std::string> d_in(nat);
putContent(d_in, element);
int storage_index = 0;
bool input_ungrouped = false;
int storage_index = 0;
for (int ig = 0; ig < nat_group.size(); ig++)
{
const auto& group_species_name = tspecies.getSpeciesName(ig);
@ -242,6 +243,8 @@ bool XMLParticleParser::readXML(xmlNodePtr cur)
" doesn't match any species from 'group' XML element nodes.");
if (element_index == ig)
{
if (iat != storage_index)
input_ungrouped = true;
count_group_size++;
map_storage_to_input[storage_index++] = iat;
}
@ -260,6 +263,22 @@ bool XMLParticleParser::readXML(xmlNodePtr cur)
throw UniformCommunicateError(msg.str());
}
}
if (input_ungrouped)
{
app_log() << " Input particle set is not grouped by species. Remapping particle position indices "
"internally."
<< std::endl;
app_debug() << " Species : input particle index -> internal particle index" << std::endl;
for (int new_idx = 0; new_idx < map_storage_to_input.size(); new_idx++)
{
int old_idx = map_storage_to_input[new_idx];
if (new_idx != old_idx)
{
app_debug() << " " << d_in[old_idx] << " : " << old_idx << " -> " << new_idx << std::endl;
}
}
}
}
});

2
src/Platforms/OMPTarget/ompBLAS.cpp
View File

@ -14,7 +14,7 @@
#include <stdexcept>
#include "config.h"
#if !defined(OPENMP_NO_COMPLEX)
#include "ompReduction.hpp"
#include "ompReductionComplex.hpp"
#endif
namespace qmcplusplus

0
src/Platforms/OMPTarget/ompReduction.hpp → src/Platforms/OMPTarget/ompReductionComplex.hpp
View File

2
src/QMCApp/QMCAppBase.cpp
View File

@ -71,7 +71,7 @@ void QMCAppBase::saveXml()
{
if (!XmlDocStack.empty())
{
std::string newxml(myProject.CurrentMainRoot());
std::string newxml(myProject.currentMainRoot());
newxml.append(".cont.xml");
app_log() << "\n========================================================="
<< "\n A new xml input file : " << newxml << std::endl;

4
src/QMCApp/QMCMain.cpp
View File

@ -307,7 +307,7 @@ bool QMCMain::execute()
HDFVersion cur_version;
v_str << cur_version[0] << " " << cur_version[1];
xmlNodePtr newmcptr = xmlNewNode(NULL, (const xmlChar*)"mcwalkerset");
xmlNewProp(newmcptr, (const xmlChar*)"fileroot", (const xmlChar*)myProject.CurrentMainRoot().c_str());
xmlNewProp(newmcptr, (const xmlChar*)"fileroot", (const xmlChar*)myProject.currentMainRoot().c_str());
xmlNewProp(newmcptr, (const xmlChar*)"node", (const xmlChar*)"-1");
xmlNewProp(newmcptr, (const xmlChar*)"nprocs", (const xmlChar*)np_str.str().c_str());
xmlNewProp(newmcptr, (const xmlChar*)"version", (const xmlChar*)v_str.str().c_str());
@ -638,7 +638,7 @@ bool QMCMain::runQMC(xmlNodePtr cur, bool reuse)
if (!FirstQMC && !append_run)
myProject.advance();
qmc_driver->setStatus(myProject.CurrentMainRoot(), "", append_run);
qmc_driver->setStatus(myProject.currentMainRoot(), "", append_run);
// PD:
// Q: How does m_walkerset_in end up being non empty?
// A: Anytime that we aren't doing a restart.

9
src/QMCDrivers/DMC/WalkerControl.cpp
View File

@ -549,14 +549,12 @@ std::vector<WalkerControl::IndexType> WalkerControl::syncFutureWalkersPerRank(Co
bool WalkerControl::put(xmlNodePtr cur)
{
int nw_target = 0, nw_max = 0;
std::string nonblocking;
std::string debug_disable_branching;
ParameterSet params;
params.add(max_copy_, "maxCopy");
params.add(nw_target, "targetwalkers");
params.add(nw_max, "max_walkers");
params.add(nonblocking, "use_nonblocking", {"yes", "no"});
params.add(debug_disable_branching, "debug_disable_branching", {"no", "yes"});
params.add(use_nonblocking_, "use_nonblocking", {true});
params.add(debug_disable_branching_, "debug_disable_branching", {false});
try
{
@ -567,9 +565,6 @@ bool WalkerControl::put(xmlNodePtr cur)
myComm->barrier_and_abort("WalkerControl::put parsing error. " + std::string(re.what()));
}
use_nonblocking_ = nonblocking == "yes";
debug_disable_branching_ = debug_disable_branching == "yes";
setMinMax(nw_target, nw_max);
app_log() << " WalkerControl parameters " << std::endl;

2
src/QMCDrivers/QMCDriverFactory.cpp
View File

@ -102,7 +102,7 @@ QMCDriverFactory::DriverAssemblyState QMCDriverFactory::readSection(xmlNodePtr c
const int nchars = qmc_mode.size();
using DV = ProjectData::DriverVersion;
switch (project_data_.get_driver_version())
switch (project_data_.getDriverVersion())
{
case DV::BATCH:
#if defined(QMC_CUDA)

1
src/QMCDrivers/QMCDriverInput.cpp
View File

@ -17,6 +17,7 @@
#include "QMCDriverInput.h"
#include "EstimatorInputDelegates.h"
#include "Concurrency/Info.hpp"
#include "ModernStringUtils.hpp"
namespace qmcplusplus
{

2
src/QMCDrivers/QMCDriverNew.h
View File

@ -308,7 +308,7 @@ protected:
/// check logpsi and grad and lap against values computed from scratch
static void checkLogAndGL(Crowd& crowd, const std::string_view location);
const std::string& get_root_name() const override { return project_data_.CurrentMainRoot(); }
const std::string& get_root_name() const override { return project_data_.currentMainRoot(); }
/** The timers for the driver.
*

8
src/QMCDrivers/WFOpt/CostFunctionCrowdData.cpp
View File

@ -19,7 +19,6 @@ CostFunctionCrowdData::CostFunctionCrowdData(int crowd_size,
ParticleSet& P,
TrialWaveFunction& Psi,
QMCHamiltonian& H,
std::vector<std::string>& H_KE_node_names,
RandomGenerator& Rng)
: h0_res_("h0 resource"), e0_(0.0), e2_(0.0), wgt_(0.0), wgt2_(0.0)
{
@ -39,9 +38,10 @@ CostFunctionCrowdData::CostFunctionCrowdData(int crowd_size,
// build a temporary H_KE for later calling makeClone
// need makeClone to setup internal my_index_ of a new copy.
QMCHamiltonian H_KE;
for (const std::string& node_name : H_KE_node_names)
H_KE.addOperator(H.getHamiltonian(node_name)->makeClone(P, Psi), node_name);
const auto components = H.getTWFDependentComponents();
QMCHamiltonian H_KE("h_free");
for (OperatorBase& component : components)
H_KE.addOperator(component.makeClone(P, Psi), component.getName());
H_KE.createResource(h0_res_);
for (int ib = 0; ib < crowd_size; ib++)

1
src/QMCDrivers/WFOpt/CostFunctionCrowdData.h
View File

@ -35,7 +35,6 @@ public:
ParticleSet& P,
TrialWaveFunction& Psi,
QMCHamiltonian& H,
std::vector<std::string>& H_KE_node_names,
RandomGenerator& Rng);
/// Set the log_psi_* arrays to zero

55
src/QMCDrivers/WFOpt/HamiltonianRef.cpp
View File

@ -19,21 +19,16 @@ namespace qmcplusplus
{
using FullPrecRealType = HamiltonianRef::FullPrecRealType;
void HamiltonianRef::addOperator(OperatorBase& op) { Hrefs_.emplace_back(op); }
HamiltonianRef::HamiltonianRef(const RefVector<OperatorBase> refs) : Hrefs_(refs) {}
FullPrecRealType HamiltonianRef::evaluateValueAndDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi,
bool compute_deriv)
Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi)
{
FullPrecRealType LocalEnergy = Hrefs_[0].get().evaluate(P);
if (compute_deriv)
for (int i = 1; i < Hrefs_.size(); ++i)
LocalEnergy += Hrefs_[i].get().evaluateValueAndDerivatives(P, optvars, dlogpsi, dhpsioverpsi);
else
for (int i = 1; i < Hrefs_.size(); ++i)
LocalEnergy += Hrefs_[i].get().evaluate(P);
FullPrecRealType LocalEnergy(0);
for (OperatorBase& Href : Hrefs_)
LocalEnergy += Href.evaluateValueAndDerivatives(P, optvars, dlogpsi, dhpsioverpsi);
return LocalEnergy;
}
@ -51,42 +46,4 @@ FullPrecRealType HamiltonianRef::evaluate(ParticleSet& P)
return LocalEnergy;
}
int HamiltonianRef::addObservables(ParticleSet& P)
{
OperatorBase::PropertySetType Observables;
//ParticleSet::mcObservables (large data, e.g. density) are accumulated while evaluations
P.Collectables.clear();
P.Collectables.rewind();
for (int i = 0; i < Hrefs_.size(); ++i)
Hrefs_[i].get().addObservables(Observables, P.Collectables);
const int myIndex = P.PropertyList.add(Observables.Names[0]);
for (int i = 1; i < Observables.size(); ++i)
P.PropertyList.add(Observables.Names[i]);
P.Collectables.size();
app_log() << "\n QMCHamiltonian::add2WalkerProperty added"
<< "\n " << Observables.size() << " to P::PropertyList "
<< "\n " << P.Collectables.size() << " to P::Collectables "
<< "\n starting Index of the observables in P::PropertyList = " << myIndex << std::endl;
return Observables.size();
}
#ifdef QMC_CUDA
void HamiltonianRef::evaluate(MCWalkerConfiguration& W, std::vector<RealType>& energyVector)
{
using WP = WalkerProperties::Indexes;
auto& walkers = W.WalkerList;
const int nw = walkers.size();
std::vector<FullPrecRealType> LocalEnergyVector(nw, 0.0);
for (int i = 0; i < Hrefs_.size(); ++i)
Hrefs_[i].get().addEnergy(W, LocalEnergyVector);
for (int iw = 0; iw < walkers.size(); iw++)
{
walkers[iw]->getPropertyBase()[WP::LOCALENERGY] = LocalEnergyVector[iw];
walkers[iw]->getPropertyBase()[WP::LOCALPOTENTIAL] =
LocalEnergyVector[iw] - walkers[iw]->getPropertyBase()[WP::NUMPROPERTIES];
energyVector[iw] = LocalEnergyVector[iw];
}
}
#endif
} // namespace qmcplusplus

17
src/QMCDrivers/WFOpt/HamiltonianRef.h
View File

@ -30,31 +30,22 @@ public:
using ValueType = OperatorBase::ValueType;
using RealType = OperatorBase::RealType;
/// record operator reference
void addOperator(OperatorBase& op);
HamiltonianRef(const RefVector<OperatorBase>);
/// the same evaluateValueAndDerivatives as QMCHamiltonian
FullPrecRealType evaluateValueAndDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi,
bool compute_deriv);
Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi);
/// the same evaluate as QMCHamiltonian
FullPrecRealType evaluate(ParticleSet& P);
#ifdef QMC_CUDA
/// the same evaluate as QMCHamiltonian
void evaluate(MCWalkerConfiguration& W, std::vector<RealType>& LocalEnergy);
#endif
int addObservables(ParticleSet& P);
int size() const { return Hrefs_.size(); }
private:
/// collected references
RefVector<OperatorBase> Hrefs_;
const RefVector<OperatorBase> Hrefs_;
};
} // namespace qmcplusplus

102
src/QMCDrivers/WFOpt/QMCCostFunction.cpp
View File

@ -179,31 +179,23 @@ void QMCCostFunction::getConfigurations(const std::string& aroot)
HDerivRecords.resize(NumThreads, 0);
}
app_log() << " Using Nonlocal PP in Opt: " << includeNonlocalH << std::endl;
outputManager.pause();
//#pragma omp parallel for
for (int ip = 0; ip < NumThreads; ++ip)
{
if (H_KE_Node[ip] == 0)
if (!H_KE_Node[ip])
{
H_KE_Node[ip] = std::make_unique<HamiltonianRef>();
H_KE_Node[ip]->addOperator(*hClones[ip]->getHamiltonian("Kinetic"));
if (includeNonlocalH != "no")
auto components = hClones[ip]->getTWFDependentComponents();
if (ip == 0)
{
OperatorBase* a(hClones[ip]->getHamiltonian(includeNonlocalH));
if (a)
{
app_log() << " Found non-local Hamiltonian element named " << includeNonlocalH << std::endl;
H_KE_Node[ip]->addOperator(*a);
}
else
app_log() << " Did not find non-local Hamiltonian element named " << includeNonlocalH << std::endl;
app_log() << " Found " << components.size() << " wavefunction dependent components in the Hamiltonian";
if (components.size())
for (const OperatorBase& component : components)
app_log() << " '" << component.getName() << "'";
app_log() << "." << std::endl;
}
H_KE_Node[ip] = std::make_unique<HamiltonianRef>(components);
}
}
//load samples from SampleStack
outputManager.resume();
app_log() << " Number of samples loaded to each thread : ";
wPerRank[0] = 0;
for (int ip = 0; ip < NumThreads; ++ip)
@ -259,8 +251,6 @@ void QMCCostFunction::checkConfigurations(EngineHandle& handle)
HDerivRecords[ip]->resize(wRef.numSamples(), NumOptimizables);
}
}
OperatorBase* nlpp = (includeNonlocalH == "no") ? nullptr : hClones[ip]->getHamiltonian(includeNonlocalH);
bool compute_nlpp = useNLPPDeriv && nlpp;
// synchronize the random number generator with the node
(*MoverRng[ip]) = (*RngSaved[ip]);
hClones[ip]->setRandomGenerator(MoverRng[ip]);
@ -273,20 +263,20 @@ void QMCCostFunction::checkConfigurations(EngineHandle& handle)
wRef.loadSample(wRef, iw);
wRef.update();
Return_rt* restrict saved = (*RecordsOnNode[ip])[iw];
psiClones[ip]->evaluateDeltaLog(wRef, saved[LOGPSI_FIXED], saved[LOGPSI_FREE], *dLogPsi[iwg], *d2LogPsi[iwg]);
psiClones[ip]->evaluateDeltaLogSetup(wRef, saved[LOGPSI_FIXED], saved[LOGPSI_FREE], *dLogPsi[iwg],
*d2LogPsi[iwg]);
saved[REWEIGHT] = 1.0;
Return_rt etmp;
if (needGrads)
{
//allocate vector
std::vector<Return_rt> rDsaved(NumOptimizables, 0.0);
std::vector<Return_rt> rHDsaved(NumOptimizables, 0.0);
Vector<Return_rt> rDsaved(NumOptimizables, 0.0);
Vector<Return_rt> rHDsaved(NumOptimizables, 0.0);
std::vector<Return_t> Dsaved(NumOptimizables, 0.0);
std::vector<Return_t> HDsaved(NumOptimizables, 0.0);
Vector<Return_t> Dsaved(NumOptimizables, 0.0);
Vector<Return_t> HDsaved(NumOptimizables, 0.0);
psiClones[ip]->evaluateDerivatives(wRef, OptVariablesForPsi, Dsaved, HDsaved);
etmp = hClones[ip]->evaluateValueAndDerivatives(wRef, OptVariablesForPsi, Dsaved, HDsaved, compute_nlpp);
etmp = hClones[ip]->evaluateValueAndDerivatives(wRef, OptVariablesForPsi, Dsaved, HDsaved);
//FIXME the ifdef should be removed after the optimizer is made compatible with complex coefficients
@ -295,17 +285,18 @@ void QMCCostFunction::checkConfigurations(EngineHandle& handle)
rDsaved[i] = std::real(Dsaved[i]);
rHDsaved[i] = std::real(HDsaved[i]);
}
copy(rDsaved.begin(), rDsaved.end(), (*DerivRecords[ip])[iw]);
copy(rHDsaved.begin(), rHDsaved.end(), (*HDerivRecords[ip])[iw]);
std::copy(rDsaved.begin(), rDsaved.end(), (*DerivRecords[ip])[iw]);
std::copy(rHDsaved.begin(), rHDsaved.end(), (*HDerivRecords[ip])[iw]);
}
else
etmp = hClones[ip]->evaluate(wRef);
e0 += saved[ENERGY_TOT] = saved[ENERGY_NEW] = etmp;
e2 += etmp * etmp;
saved[ENERGY_FIXED] = hClones[ip]->getLocalPotential();
if (nlpp)
saved[ENERGY_FIXED] -= nlpp->getValue();
saved[ENERGY_FIXED] = saved[ENERGY_TOT];
const auto twf_dependent_components = hClones[ip]->getTWFDependentComponents();
for (const OperatorBase& component : twf_dependent_components)
saved[ENERGY_FIXED] -= component.getValue();
}
//add them all using reduction
et_tot += e0;
@ -382,8 +373,6 @@ void QMCCostFunction::engine_checkConfigurations(cqmc::engine::LMYEngine<Return_
//HDerivRecords[ip]->resize(wRef.numSamples(),NumOptimizables);
}
}
OperatorBase* nlpp = (includeNonlocalH == "no") ? nullptr : hClones[ip]->getHamiltonian(includeNonlocalH.c_str());
bool compute_nlpp = useNLPPDeriv && nlpp;
// synchronize the random number generator with the node
(*MoverRng[ip]) = (*RngSaved[ip]);
hClones[ip]->setRandomGenerator(MoverRng[ip]);
@ -398,17 +387,17 @@ void QMCCostFunction::engine_checkConfigurations(cqmc::engine::LMYEngine<Return_
wRef.loadSample(wRef, iw);
wRef.update();
Return_rt* restrict saved = (*RecordsOnNode[ip])[iw];
psiClones[ip]->evaluateDeltaLog(wRef, saved[LOGPSI_FIXED], saved[LOGPSI_FREE], *dLogPsi[iwg], *d2LogPsi[iwg]);
psiClones[ip]->evaluateDeltaLogSetup(wRef, saved[LOGPSI_FIXED], saved[LOGPSI_FREE], *dLogPsi[iwg],
*d2LogPsi[iwg]);
saved[REWEIGHT] = 1.0;
Return_rt etmp;
if (needGrads)
{
//allocate vector
std::vector<Return_t> Dsaved(NumOptimizables, 0.0);
std::vector<Return_t> HDsaved(NumOptimizables, 0.0);
Vector<Return_t> Dsaved(NumOptimizables, 0.0);
Vector<Return_t> HDsaved(NumOptimizables, 0.0);
psiClones[ip]->evaluateDerivatives(wRef, OptVariablesForPsi, Dsaved, HDsaved);
etmp = hClones[ip]->evaluateValueAndDerivatives(wRef, OptVariablesForPsi, Dsaved, HDsaved, compute_nlpp);
etmp = hClones[ip]->evaluateValueAndDerivatives(wRef, OptVariablesForPsi, Dsaved, HDsaved);
// add non-differentiated derivative vector
std::vector<Return_t> der_rat_samp(NumOptimizables + 1, 0.0);
@ -417,12 +406,12 @@ void QMCCostFunction::engine_checkConfigurations(cqmc::engine::LMYEngine<Return_
// dervative vectors
der_rat_samp.at(0) = 1.0;
for (int i = 0; i < Dsaved.size(); i++)
der_rat_samp.at(i + 1) = Dsaved.at(i);
der_rat_samp[i + 1] = Dsaved[i];
// energy dervivatives
le_der_samp.at(0) = etmp;
for (int i = 0; i < HDsaved.size(); i++)
le_der_samp.at(i + 1) = HDsaved.at(i) + etmp * Dsaved.at(i);
le_der_samp[i + 1] = HDsaved[i] + etmp * Dsaved[i];
#ifdef HAVE_LMY_ENGINE
if (MinMethod == "adaptive")
@ -446,9 +435,11 @@ void QMCCostFunction::engine_checkConfigurations(cqmc::engine::LMYEngine<Return_
e0 += saved[ENERGY_TOT] = etmp;
e2 += etmp * etmp;
saved[ENERGY_FIXED] = hClones[ip]->getLocalPotential();
if (nlpp)
saved[ENERGY_FIXED] -= nlpp->getValue();
saved[ENERGY_FIXED] = saved[ENERGY_TOT];
const auto twf_dependent_components = hClones[ip]->getTWFDependentComponents();
for (const OperatorBase& component : twf_dependent_components)
saved[ENERGY_FIXED] -= component.getValue();
}
//add them all using reduction
@ -509,9 +500,10 @@ void QMCCostFunction::resetPsi(bool final_reset)
//cout << "######### QMCCostFunction::resetPsi " << std::endl;
//OptVariablesForPsi.print(std::cout);
//cout << "-------------------------------------- " << std::endl;
Psi.resetParameters(OptVariablesForPsi);
resetOptimizableObjects(Psi, OptVariablesForPsi);
for (int i = 0; i < psiClones.size(); ++i)
psiClones[i]->resetParameters(OptVariablesForPsi);
resetOptimizableObjects(*psiClones[i], OptVariablesForPsi);
}
QMCCostFunction::EffectiveWeight QMCCostFunction::correlatedSampling(bool needGrad)
@ -523,15 +515,14 @@ QMCCostFunction::EffectiveWeight QMCCostFunction::correlatedSampling(bool needGr
hClones[ip]->setRandomGenerator(MoverRng[ip]);
}
const bool nlpp = (includeNonlocalH != "no");
Return_rt wgt_tot = 0.0;
Return_rt wgt_tot2 = 0.0;
Return_rt inv_n_samples = 1.0 / NumSamples;
#pragma omp parallel reduction(+ : wgt_tot, wgt_tot2)
{
int ip = omp_get_thread_num();
bool compute_nlpp = useNLPPDeriv && (includeNonlocalH != "no");
bool compute_all_from_scratch = (includeNonlocalH != "no"); //true if we have nlpp
const int ip = omp_get_thread_num();
//if we have more than KE depending on TWF, TWF must be fully recomputed.
const bool compute_all_from_scratch = hClones[ip]->getTWFDependentComponents().size() > 1;
MCWalkerConfiguration& wRef(*wClones[ip]);
Return_rt wgt_node = 0.0, wgt_node2 = 0.0;
@ -547,16 +538,14 @@ QMCCostFunction::EffectiveWeight QMCCostFunction::correlatedSampling(bool needGr
Return_rt weight = saved[REWEIGHT] = vmc_or_dmc * (logpsi - saved[LOGPSI_FREE]);
if (needGrad)
{
std::vector<Return_t> Dsaved(NumOptimizables, 0);
std::vector<Return_t> HDsaved(NumOptimizables, 0);
Vector<Return_t> Dsaved(NumOptimizables, 0);
Vector<Return_t> HDsaved(NumOptimizables, 0);
std::vector<Return_rt> rDsaved(NumOptimizables, 0);
std::vector<Return_rt> rHDsaved(NumOptimizables, 0);
psiClones[ip]->evaluateDerivatives(wRef, OptVariablesForPsi, Dsaved, HDsaved);
Vector<Return_rt> rDsaved(NumOptimizables, 0);
Vector<Return_rt> rHDsaved(NumOptimizables, 0);
saved[ENERGY_NEW] =
H_KE_Node[ip]->evaluateValueAndDerivatives(wRef, OptVariablesForPsi, Dsaved, HDsaved, compute_nlpp) +
saved[ENERGY_FIXED];
H_KE_Node[ip]->evaluateValueAndDerivatives(wRef, OptVariablesForPsi, Dsaved, HDsaved) + saved[ENERGY_FIXED];
;
for (int i = 0; i < NumOptimizables; i++)
@ -571,7 +560,6 @@ QMCCostFunction::EffectiveWeight QMCCostFunction::correlatedSampling(bool needGr
(*DerivRecords[ip])(iw, i) = rDsaved[i];
(*HDerivRecords[ip])(iw, i) = rHDsaved[i];
}
//saved[ENERGY_NEW] = H_KE_Node[ip]->evaluate(wRef) + saved[ENERGY_FIXED];
}
else
saved[ENERGY_NEW] = H_KE_Node[ip]->evaluate(wRef) + saved[ENERGY_FIXED];

98
src/QMCDrivers/WFOpt/QMCCostFunctionBase.cpp
View File

@ -21,16 +21,13 @@
#include "OhmmsData/ParameterSet.h"
#include "OhmmsData/XMLParsingString.h"
#include "Message/CommOperators.h"
#include <set>
#include "Message/UniformCommunicateError.h"
//#define QMCCOSTFUNCTION_DEBUG
namespace qmcplusplus
{
QMCCostFunctionBase::QMCCostFunctionBase(ParticleSet& w,
TrialWaveFunction& psi,
QMCHamiltonian& h,
Communicate* comm)
QMCCostFunctionBase::QMCCostFunctionBase(ParticleSet& w, TrialWaveFunction& psi, QMCHamiltonian& h, Communicate* comm)
: MPIObjectBase(comm),
reportH5(false),
CI_Opt(false),
@ -56,7 +53,6 @@ QMCCostFunctionBase::QMCCostFunctionBase(ParticleSet& w,
msg_stream(0),
m_wfPtr(NULL),
m_doc_out(NULL),
includeNonlocalH("no"),
debug_stream(0)
{
GEVType = "mixed";
@ -65,8 +61,7 @@ QMCCostFunctionBase::QMCCostFunctionBase(ParticleSet& w,
//default: don't check fo MinNumWalkers
MinNumWalkers = 0.3;
SumValue.resize(SUM_INDEX_SIZE, 0.0);
IsValid = true;
useNLPPDeriv = false;
IsValid = true;
#if defined(QMCCOSTFUNCTION_DEBUG)
char fname[16];
sprintf(fname, "optdebug.p%d", OHMMS::Controller->mycontext());
@ -131,7 +126,7 @@ QMCCostFunctionBase::Return_rt QMCCostFunctionBase::Cost(bool needGrad)
resetPsi();
//evaluate new local energies
EffectiveWeight effective_weight = correlatedSampling(needGrad);
IsValid = isEffectiveWeightValid(effective_weight);
IsValid = isEffectiveWeightValid(effective_weight);
return computedCost();
}
@ -191,9 +186,9 @@ void QMCCostFunctionBase::Report()
if (msg_stream)
{
msg_stream->precision(8);
*msg_stream << " curCost " << std::setw(5) << ReportCounter << std::setw(16) << CostValue
<< std::setw(16) << curAvg_w << std::setw(16) << curAvg << std::setw(16) << curVar_w
<< std::setw(16) << curVar << std::setw(16) << curVar_abs << std::endl;
*msg_stream << " curCost " << std::setw(5) << ReportCounter << std::setw(16) << CostValue << std::setw(16)
<< curAvg_w << std::setw(16) << curAvg << std::setw(16) << curVar_w << std::setw(16) << curVar
<< std::setw(16) << curVar_abs << std::endl;
*msg_stream << " curVars " << std::setw(5) << ReportCounter;
for (int i = 0; i < OptVariables.size(); i++)
*msg_stream << std::setw(16) << OptVariables[i];
@ -312,36 +307,41 @@ bool QMCCostFunctionBase::checkParameters()
*/
bool QMCCostFunctionBase::put(xmlNodePtr q)
{
std::string includeNonlocalH;
std::string writeXmlPerStep("no");
std::string computeNLPPderiv;
std::string output_override_str("no");
astring variational_subset_str;
ParameterSet m_param;
m_param.add(writeXmlPerStep, "dumpXML");
m_param.add(MinNumWalkers, "minwalkers");
m_param.add(MaxWeight, "maxWeight");
m_param.add(includeNonlocalH, "nonlocalpp");
m_param.add(computeNLPPderiv, "use_nonlocalpp_deriv", {"yes", "no"});
m_param.add(includeNonlocalH, "nonlocalpp", {}, TagStatus::DEPRECATED);
m_param.add(computeNLPPderiv, "use_nonlocalpp_deriv", {}, TagStatus::DEPRECATED);
m_param.add(w_beta, "beta");
m_param.add(GEVType, "GEVMethod");
m_param.add(targetExcitedStr, "targetExcited");
m_param.add(omega_shift, "omega");
m_param.add(output_override_str, "output_vp_override", {"no", "yes"});
m_param.add(variational_subset_str, "variational_subset");
m_param.put(q);
if (!includeNonlocalH.empty())
app_warning() << "'nonlocalpp' no more affects any part of the execution. Please remove it from your input file."
<< std::endl;
if (!computeNLPPderiv.empty())
app_warning()
<< "'use_nonlocalpp_deriv' no more affects any part of the execution. Please remove it from your input file."
<< std::endl;
targetExcitedStr = lowerCase(targetExcitedStr);
targetExcited = (targetExcitedStr == "yes");
if (output_override_str == "yes")
do_override_output = true;
if (includeNonlocalH == "yes")
includeNonlocalH = "NonLocalECP";
variational_subset_names = convertStrToVec<std::string>(variational_subset_str.s);
if (computeNLPPderiv != "no" && includeNonlocalH != "no")
{
app_log() << " Going to include the derivatives of " << includeNonlocalH << std::endl;
useNLPPDeriv = true;
}
// app_log() << " QMCCostFunctionBase::put " << std::endl;
// m_param.get(app_log());
Write2OneXml = (writeXmlPerStep == "no");
@ -399,10 +399,21 @@ bool QMCCostFunctionBase::put(xmlNodePtr q)
}
cur = cur->next;
}
UniqueOptObjRefs opt_obj_refs = extractOptimizableObjects(Psi);
app_log() << " TrialWaveFunction \"" << Psi.getName() << "\" has " << opt_obj_refs.size()
<< " optimizable objects:" << std::endl;
for (OptimizableObject& obj : opt_obj_refs)
app_log() << " '" << obj.getName() << "'" << (obj.isOptimized() ? " optimized" : " fixed") << std::endl;
//build optimizables from the wavefunction
OptVariablesForPsi.clear();
Psi.checkInVariables(OptVariablesForPsi);
for (OptimizableObject& obj : opt_obj_refs)
if (obj.isOptimized())
obj.checkInVariablesExclusive(OptVariablesForPsi);
OptVariablesForPsi.resetIndex();
app_log() << " Variational subset selects " << OptVariablesForPsi.size() << " parameters." << std::endl;
//synchronize OptVariables and OptVariablesForPsi
OptVariables = OptVariablesForPsi;
InitVariables = OptVariablesForPsi;
@ -469,6 +480,9 @@ bool QMCCostFunctionBase::put(xmlNodePtr q)
{
APP_ABORT("QMCCostFunctionBase::put No valid optimizable variables are found.");
}
else
app_log() << " In total " << NumOptimizables << " parameters being optimized after applying constraints."
<< std::endl;
// app_log() << "<active-optimizables> " << std::endl;
// OptVariables.print(app_log());
// app_log() << "</active-optimizables>" << std::endl;
@ -1036,11 +1050,12 @@ void QMCCostFunctionBase::printCJParams(xmlNodePtr cur, std::string& rname)
bool QMCCostFunctionBase::isEffectiveWeightValid(EffectiveWeight effective_weight) const
{
app_log() << "Effective weight of all the samples measured by correlated sampling is "
<< effective_weight << std::endl;
app_log() << "Effective weight of all the samples measured by correlated sampling is " << effective_weight
<< std::endl;
if (effective_weight < MinNumWalkers)
{
WARNMSG(" Smaller than the user specified threshold \"minwalkers\" = " << MinNumWalkers << std::endl
WARNMSG(" Smaller than the user specified threshold \"minwalkers\" = "
<< MinNumWalkers << std::endl
<< " If this message appears frequently. You might have to be cautious. " << std::endl
<< " Find info about parameter \"minwalkers\" in the user manual!");
return false;
@ -1049,6 +1064,39 @@ bool QMCCostFunctionBase::isEffectiveWeightValid(EffectiveWeight effective_weigh
return true;
}
UniqueOptObjRefs QMCCostFunctionBase::extractOptimizableObjects(TrialWaveFunction& psi) const
{
const auto& names(variational_subset_names);
// survey all the optimizable objects
const auto opt_obj_refs = psi.extractOptimizableObjectRefs();
// check if input names are valid
for (auto& name : names)
if (std::find_if(opt_obj_refs.begin(), opt_obj_refs.end(),
[&name](const OptimizableObject& obj) { return name == obj.getName(); }) == opt_obj_refs.end())
{
std::ostringstream msg;
msg << "Variational subset entry '" << name << "' doesn't exist in the trial wavefunction which contains";
for (OptimizableObject& obj : opt_obj_refs)
msg << " '" << obj.getName() << "'";
msg << "." << std::endl;
throw UniformCommunicateError(msg.str());
}
for (OptimizableObject& obj : opt_obj_refs)
obj.setOptimization(names.empty() || std::find_if(names.begin(), names.end(), [&obj](const std::string& name) {
return name == obj.getName();
}) != names.end());
return opt_obj_refs;
}
void QMCCostFunctionBase::resetOptimizableObjects(TrialWaveFunction& psi, const opt_variables_type& opt_variables) const
{
const auto opt_obj_refs = extractOptimizableObjects(psi);
for (OptimizableObject& obj : opt_obj_refs)
if (obj.isOptimized())
obj.resetParametersExclusive(opt_variables);
}
///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief If the LMYEngine is available, returns the cost function calculated by the engine.
/// Otherwise, returns the usual cost function.

10
src/QMCDrivers/WFOpt/QMCCostFunctionBase.h
View File

@ -184,8 +184,6 @@ protected:
///if true, do not write the *.opt.#.xml
bool Write2OneXml;
///if true, use analytic derivatives for the non-local potential component
bool useNLPPDeriv;
/** |E-E_T|^PowerE is used for the cost function
*
* default PowerE=1
@ -276,8 +274,9 @@ protected:
///Random number generators
UPtrVector<RandomGenerator> RngSaved;
std::vector<RandomGenerator*> MoverRng;
std::string includeNonlocalH;
/// optimized parameter names
std::vector<std::string> variational_subset_names;
/** Sum of energies and weights for averages
*
@ -318,6 +317,11 @@ protected:
/// check the validity of the effective weight calculated by correlatedSampling
bool isEffectiveWeightValid(EffectiveWeight effective_weight) const;
/// survey all the optimizable objects
UniqueOptObjRefs extractOptimizableObjects(TrialWaveFunction& psi) const;
void resetOptimizableObjects(TrialWaveFunction& psi, const opt_variables_type& opt_variables) const;
#ifdef HAVE_LMY_ENGINE
virtual Return_rt LMYEngineCost_detail(cqmc::engine::LMYEngine<Return_t>* EngineObj)
{

73
src/QMCDrivers/WFOpt/QMCCostFunctionBatched.cpp
View File

@ -171,27 +171,12 @@ void QMCCostFunctionBatched::GradCost(std::vector<Return_rt>& PGradient,
void QMCCostFunctionBatched::getConfigurations(const std::string& aroot)
{
app_log() << " Using Nonlocal PP in Opt: " << includeNonlocalH << std::endl;
outputManager.pause();
if (H_KE_node_names_.size() == 0)
{
H_KE_node_names_.reserve(2);
H_KE_node_names_.emplace_back("Kinetic");
if (includeNonlocalH != "no")
{
OperatorBase* a(H.getHamiltonian(includeNonlocalH));
outputManager.resume();
if (a)
{
app_log() << " Found non-local Hamiltonian element named " << includeNonlocalH << std::endl;
H_KE_node_names_.emplace_back(includeNonlocalH);
}
else
app_log() << " Did not find non-local Hamiltonian element named " << includeNonlocalH << std::endl;
}
}
outputManager.resume();
auto components = H.getTWFDependentComponents();
app_log() << " Found " << components.size() << " wavefunction dependent components in the Hamiltonian";
if (components.size())
for(const OperatorBase& component: components)
app_log() << " '" << component.getName() << "'";
app_log() << "." << std::endl;
rank_local_num_samples_ = samples_.getNumSamples();
@ -256,8 +241,6 @@ void QMCCostFunctionBatched::checkConfigurations(EngineHandle& handle)
HDerivRecords_.resize(rank_local_num_samples_, NumOptimizables);
}
}
OperatorBase* nlpp = (includeNonlocalH == "no") ? nullptr : H.getHamiltonian(includeNonlocalH);
bool compute_nlpp = useNLPPDeriv && nlpp;
// synchronize the random number generator with the node
(*MoverRng[0]) = (*RngSaved[0]);
H.setRandomGenerator(MoverRng[0]);
@ -269,7 +252,7 @@ void QMCCostFunctionBatched::checkConfigurations(EngineHandle& handle)
opt_eval_.resize(opt_num_crowds);
for (int i = 0; i < opt_num_crowds; i++)
opt_eval_[i] =
std::make_unique<CostFunctionCrowdData>(walkers_per_crowd_[i], W, Psi, H, H_KE_node_names_, *MoverRng[0]);
std::make_unique<CostFunctionCrowdData>(walkers_per_crowd_[i], W, Psi, H, *MoverRng[0]);
outputManager.resume();
@ -286,8 +269,7 @@ void QMCCostFunctionBatched::checkConfigurations(EngineHandle& handle)
std::vector<ParticleGradient*>& gradPsi, std::vector<ParticleLaplacian*>& lapPsi,
Matrix<Return_rt>& RecordsOnNode, Matrix<Return_rt>& DerivRecords,
Matrix<Return_rt>& HDerivRecords, const SampleStack& samples, opt_variables_type& optVars,
bool needGrads, bool compute_nlpp, const std::string& includeNonlocalH,
EngineHandle& handle) {
bool needGrads, EngineHandle& handle) {
CostFunctionCrowdData& opt_data = *opt_crowds[crowd_id];
const int local_samples = samples_per_crowd_offsets[crowd_id + 1] - samples_per_crowd_offsets[crowd_id];
@ -352,12 +334,10 @@ void QMCCostFunctionBatched::checkConfigurations(EngineHandle& handle)
int nparam = optVars.size();
RecordArray<Return_t> dlogpsi_array(nparam, current_batch_size);
RecordArray<Return_t> dhpsioverpsi_array(nparam, current_batch_size);
TrialWaveFunction::mw_evaluateParameterDerivatives(wf_list, p_list, optVars, dlogpsi_array, dhpsioverpsi_array);
auto energy_list =
QMCHamiltonian::mw_evaluateValueAndDerivatives(h_list, wf_list, p_list, optVars, dlogpsi_array,
dhpsioverpsi_array, compute_nlpp);
dhpsioverpsi_array);
for (int ib = 0; ib < current_batch_size; ib++)
{
@ -385,12 +365,9 @@ void QMCCostFunctionBatched::checkConfigurations(EngineHandle& handle)
RecordsOnNode[is][REWEIGHT] = 1.0;
RecordsOnNode[is][ENERGY_FIXED] = h_list[ib].getLocalPotential();
if (includeNonlocalH != "no")
{
OperatorBase* nlpp = h_list[ib].getHamiltonian(includeNonlocalH);
if (nlpp)
RecordsOnNode[is][ENERGY_FIXED] -= nlpp->getValue();
}
const auto twf_dependent_components = h_list[ib].getTWFDependentComponents();
for (const OperatorBase& component: twf_dependent_components)
RecordsOnNode[is][ENERGY_FIXED] -= component.getValue();
}
}
else
@ -417,7 +394,7 @@ void QMCCostFunctionBatched::checkConfigurations(EngineHandle& handle)
ParallelExecutor<> crowd_tasks;
crowd_tasks(opt_num_crowds, evalOptConfig, opt_eval_, samples_per_crowd_offsets, walkers_per_crowd_, dLogPsi,
d2LogPsi, RecordsOnNode_, DerivRecords_, HDerivRecords_, samples_, OptVariablesForPsi, needGrads,
compute_nlpp, includeNonlocalH, handle);
handle);
// Sum energy values over crowds
for (int i = 0; i < opt_eval_.size(); i++)
{
@ -479,14 +456,10 @@ void QMCCostFunctionBatched::resetPsi(bool final_reset)
//cout << "######### QMCCostFunctionBatched::resetPsi " << std::endl;
//OptVariablesForPsi.print(std::cout);
//cout << "-------------------------------------- " << std::endl;
Psi.resetParameters(OptVariablesForPsi);
resetOptimizableObjects(Psi, OptVariablesForPsi);
for (int i = 0; i < opt_eval_.size(); i++)
{
for (int j = 0; j < opt_eval_[i]->get_wf_ptr_list().size(); j++)
{
opt_eval_[i]->get_wf_ptr_list()[j]->resetParameters(OptVariablesForPsi);
}
}
resetOptimizableObjects(*opt_eval_[i]->get_wf_ptr_list()[j], OptVariablesForPsi);
}
QMCCostFunctionBatched::EffectiveWeight QMCCostFunctionBatched::correlatedSampling(bool needGrad)
@ -500,7 +473,6 @@ QMCCostFunctionBatched::EffectiveWeight QMCCostFunctionBatched::correlatedSampli
}
//Return_rt wgt_node = 0.0, wgt_node2 = 0.0;
const bool nlpp = (includeNonlocalH != "no");
Return_rt wgt_tot = 0.0;
Return_rt wgt_tot2 = 0.0;
@ -509,9 +481,6 @@ QMCCostFunctionBatched::EffectiveWeight QMCCostFunctionBatched::correlatedSampli
Return_rt inv_n_samples = 1.0 / samples_.getGlobalNumSamples();
bool compute_nlpp = useNLPPDeriv && (includeNonlocalH != "no");
bool compute_all_from_scratch = (includeNonlocalH != "no"); //true if we have nlpp
const size_t opt_num_crowds = walkers_per_crowd_.size();
// Divide samples among crowds
std::vector<int> samples_per_crowd_offsets(opt_num_crowds + 1);
@ -524,10 +493,9 @@ QMCCostFunctionBatched::EffectiveWeight QMCCostFunctionBatched::correlatedSampli
std::vector<ParticleLaplacian*>& lapPsi, Matrix<Return_rt>& RecordsOnNode,
Matrix<Return_rt>& DerivRecords, Matrix<Return_rt>& HDerivRecords,
const SampleStack& samples, const opt_variables_type& optVars,
bool compute_all_from_scratch, Return_rt vmc_or_dmc, bool needGrad, bool compute_nlpp) {
bool compute_all_from_scratch, Return_rt vmc_or_dmc, bool needGrad) {
CostFunctionCrowdData& opt_data = *opt_crowds[crowd_id];
const int local_samples = samples_per_crowd_offsets[crowd_id + 1] - samples_per_crowd_offsets[crowd_id];
int num_batches;
@ -615,13 +583,10 @@ QMCCostFunctionBatched::EffectiveWeight QMCCostFunctionBatched::correlatedSampli
RecordArray<Return_t> dlogpsi_array(nparam, current_batch_size);
RecordArray<Return_t> dhpsioverpsi_array(nparam, current_batch_size);
TrialWaveFunction::mw_evaluateParameterDerivatives(wf_list, p_list, optVars, dlogpsi_array, dhpsioverpsi_array);
// Energy
auto energy_list =
QMCHamiltonian::mw_evaluateValueAndDerivatives(h0_list, wf_list, p_list, optVars, dlogpsi_array,
dhpsioverpsi_array, compute_nlpp);
dhpsioverpsi_array);
for (int ib = 0; ib < current_batch_size; ib++)
{
@ -652,10 +617,12 @@ QMCCostFunctionBatched::EffectiveWeight QMCCostFunctionBatched::correlatedSampli
}
};
//if we have more than KE depending on TWF, TWF must be fully recomputed.
const bool compute_all_from_scratch = H.getTWFDependentComponents().size() > 1;
ParallelExecutor<> crowd_tasks;
crowd_tasks(opt_num_crowds, evalOptCorrelated, opt_eval_, samples_per_crowd_offsets, walkers_per_crowd_, dLogPsi,
d2LogPsi, RecordsOnNode_, DerivRecords_, HDerivRecords_, samples_, OptVariablesForPsi,
compute_all_from_scratch, vmc_or_dmc, needGrad, compute_nlpp);
compute_all_from_scratch, vmc_or_dmc, needGrad);
// Sum weights over crowds
for (int i = 0; i < opt_eval_.size(); i++)
{

18
src/QMCDrivers/WaveFunctionTester.cpp
View File

@ -781,7 +781,7 @@ bool WaveFunctionTester::checkGradientAtConfiguration(MCWalkerConfiguration::Wal
LogValueType logpsi1 = orb->evaluateLog(W, G, L);
fail_log << "WaveFunctionComponent " << iorb << " " << orb->ClassName << " log psi = " << logpsi1 << std::endl;
fail_log << "WaveFunctionComponent " << iorb << " " << orb->getClassName() << " log psi = " << logpsi1 << std::endl;
FiniteDifference::ValueVector logpsi_vals;
FiniteDifference::PosChangeVector::iterator it;
@ -807,7 +807,7 @@ bool WaveFunctionTester::checkGradientAtConfiguration(MCWalkerConfiguration::Wal
}
fd.computeFiniteDiff(delta, positions, logpsi_vals, G1, L1);
fout << " WaveFunctionComponent " << iorb << " " << orb->ClassName << std::endl;
fout << " WaveFunctionComponent " << iorb << " " << orb->getClassName() << std::endl;
if (!checkGradients(0, nat, G, L, G1, L1, fail_log, 1))
{
@ -1726,8 +1726,8 @@ void WaveFunctionTester::runDerivTest()
wfvar_prime = wfVars;
wfVars.print(fout);
int Nvars = wfVars.size();
std::vector<ValueType> Dsaved(Nvars);
std::vector<ValueType> HDsaved(Nvars);
Vector<ValueType> Dsaved(Nvars);
Vector<ValueType> HDsaved(Nvars);
std::vector<RealType> PGradient(Nvars);
std::vector<RealType> HGradient(Nvars);
Psi.resetParameters(wfVars);
@ -1833,8 +1833,8 @@ void WaveFunctionTester::runDerivNLPPTest()
wfvar_prime = wfVars;
wfVars.print(nlout);
int Nvars = wfVars.size();
std::vector<ValueType> Dsaved(Nvars);
std::vector<ValueType> HDsaved(Nvars);
Vector<ValueType> Dsaved(Nvars);
Vector<ValueType> HDsaved(Nvars);
std::vector<RealType> PGradient(Nvars);
std::vector<RealType> HGradient(Nvars);
Psi.resetParameters(wfVars);
@ -1847,7 +1847,7 @@ void WaveFunctionTester::runDerivNLPPTest()
std::vector<RealType> ene(4), ene_p(4), ene_m(4);
Psi.evaluateDerivatives(W, wfVars, Dsaved, HDsaved);
ene[0] = H.evaluateValueAndDerivatives(W, wfVars, Dsaved, HDsaved, true);
ene[0] = H.evaluateValueAndDerivatives(W, wfVars, Dsaved, HDsaved);
app_log() << "Check the energy " << eloc << " " << H.getLocalEnergy() << " " << ene[0] << std::endl;
RealType FiniteDiff = 1e-6;
@ -1935,8 +1935,8 @@ void WaveFunctionTester::runDerivCloneTest()
wfvar_prime.print(fout);
psi_clone->resetParameters(wfvar_prime);
Psi.resetParameters(wfVars);
std::vector<ValueType> Dsaved(Nvars, 0), og_Dsaved(Nvars, 0);
std::vector<ValueType> HDsaved(Nvars, 0), og_HDsaved(Nvars, 0);
Vector<ValueType> Dsaved(Nvars, 0), og_Dsaved(Nvars, 0);
Vector<ValueType> HDsaved(Nvars, 0), og_HDsaved(Nvars, 0);
std::vector<RealType> PGradient(Nvars, 0), og_PGradient(Nvars, 0);
std::vector<RealType> HGradient(Nvars, 0), og_HGradient(Nvars, 0);
ValueType logpsi2 = psi_clone->evaluateLog(*w_clone);

4
src/QMCDrivers/tests/test_dmc.cpp
View File

@ -82,7 +82,7 @@ TEST_CASE("DMC Particle-by-Particle advanceWalkers ConstantOrbital", "[drivers][
FakeRandom rg;
QMCHamiltonian h;
h.addOperator(std::make_unique<BareKineticEnergy>(elec), "Kinetic");
h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
h.addObservables(elec); // get double free error on 'h.Observables' w/o this
elec.resetWalkerProperty(); // get memory corruption w/o this
@ -174,7 +174,7 @@ TEST_CASE("DMC Particle-by-Particle advanceWalkers LinearOrbital", "[drivers][dm
FakeRandom rg;
QMCHamiltonian h;
h.addOperator(std::make_unique<BareKineticEnergy>(elec), "Kinetic");
h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
h.addObservables(elec); // get double free error on 'h.Observables' w/o this
elec.resetWalkerProperty(); // get memory corruption w/o this

4
src/QMCDrivers/tests/test_dmc_driver.cpp
View File

@ -84,7 +84,7 @@ TEST_CASE("DMC", "[drivers][dmc]")
FakeRandom rg;
QMCHamiltonian h;
std::unique_ptr<BareKineticEnergy> p_bke = std::make_unique<BareKineticEnergy>(elec);
std::unique_ptr<BareKineticEnergy> p_bke = std::make_unique<BareKineticEnergy>(elec, psi);
h.addOperator(std::move(p_bke), "Kinetic");
h.addObservables(elec); // get double free error on 'h.Observables' w/o this
@ -169,7 +169,7 @@ TEST_CASE("SODMC", "[drivers][dmc]")
FakeRandom rg;
QMCHamiltonian h;
std::unique_ptr<BareKineticEnergy> p_bke = std::make_unique<BareKineticEnergy>(elec);
std::unique_ptr<BareKineticEnergy> p_bke = std::make_unique<BareKineticEnergy>(elec, psi);
h.addOperator(std::move(p_bke), "Kinetic");
h.addObservables(elec); // get double free error on 'h.Observables' w/o this

2
src/QMCDrivers/tests/test_vmc.cpp
View File

@ -80,7 +80,7 @@ TEST_CASE("VMC Particle-by-Particle advanceWalkers", "[drivers][vmc]")
FakeRandom rg;
QMCHamiltonian h;
h.addOperator(std::make_unique<BareKineticEnergy>(elec), "Kinetic");
h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
h.addObservables(elec); // get double free error on 'h.Observables' w/o this
elec.resetWalkerProperty(); // get memory corruption w/o this

6
src/QMCDrivers/tests/test_vmc_driver.cpp
View File

@ -77,7 +77,7 @@ TEST_CASE("VMC", "[drivers][vmc]")
FakeRandom rg;
QMCHamiltonian h;
h.addOperator(std::make_unique<BareKineticEnergy>(elec), "Kinetic");
h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
h.addObservables(elec); // get double free error on 'h.Observables' w/o this
elec.resetWalkerProperty(); // get memory corruption w/o this
@ -158,7 +158,7 @@ TEST_CASE("SOVMC", "[drivers][vmc]")
FakeRandom rg;
QMCHamiltonian h;
h.addOperator(std::make_unique<BareKineticEnergy>(elec), "Kinetic");
h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
h.addObservables(elec); // get double free error on 'h.Observables' w/o this
elec.resetWalkerProperty(); // get memory corruption w/o this
@ -244,7 +244,7 @@ TEST_CASE("SOVMC-alle", "[drivers][vmc]")
FakeRandom rg;
QMCHamiltonian h;
h.addOperator(std::make_unique<BareKineticEnergy>(elec), "Kinetic");
h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
h.addObservables(elec); // get double free error on 'h.Observables' w/o this
elec.resetWalkerProperty(); // get memory corruption w/o this

3
src/QMCHamiltonians/ACForce.h
View File

@ -34,6 +34,9 @@ public:
/** Destructor, "final" triggers a clang warning **/
~ACForce() override = default;
bool dependsOnWaveFunction() const override { return true; }
std::string getClassName() const override { return "ACForce"; }
/** I/O Routines */
bool put(xmlNodePtr cur) final;

31
src/QMCHamiltonians/BareKineticEnergy.cpp
View File

@ -37,7 +37,7 @@ using Return_t = BareKineticEnergy::Return_t;
* Store mass per species and use SameMass to choose the methods.
* if SameMass, probably faster and easy to vectorize but no impact on the performance.
*/
BareKineticEnergy::BareKineticEnergy(ParticleSet& p) : Ps(p)
BareKineticEnergy::BareKineticEnergy(ParticleSet& p, TrialWaveFunction& psi) : Ps(p), psi_(psi)
{
setEnergyDomain(KINETIC);
oneBodyQuantumDomain(p);
@ -127,6 +127,33 @@ Return_t BareKineticEnergy::evaluate(ParticleSet& P)
return value_;
}
Return_t BareKineticEnergy::evaluateValueAndDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
const Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi)
{
// const_cast is needed because TWF::evaluateDerivatives calculates dlogpsi.
// KineticEnergy must be the first element in the hamiltonian array.
psi_.evaluateDerivatives(P, optvars, const_cast<Vector<ValueType>&>(dlogpsi), dhpsioverpsi);
return evaluate(P);
}
void BareKineticEnergy::mw_evaluateWithParameterDerivatives(const RefVectorWithLeader<OperatorBase>& o_list,
const RefVectorWithLeader<ParticleSet>& p_list,
const opt_variables_type& optvars,
const RecordArray<ValueType>& dlogpsi,
RecordArray<ValueType>& dhpsioverpsi) const
{
RefVectorWithLeader<TrialWaveFunction> wf_list(o_list.getCastedLeader<BareKineticEnergy>().psi_);
for (int i = 0; i < o_list.size(); i++)
wf_list.push_back(o_list.getCastedElement<BareKineticEnergy>(i).psi_);
mw_evaluate(o_list, wf_list, p_list);
// const_cast is needed because TWF::evaluateDerivatives calculates dlogpsi.
// KineticEnergy must be the first element in the hamiltonian array.
TrialWaveFunction::mw_evaluateParameterDerivatives(wf_list, p_list, optvars,
const_cast<RecordArray<ValueType>&>(dlogpsi), dhpsioverpsi);
}
/**@brief Function to compute the value, direct ionic gradient terms, and pulay terms for the local kinetic energy.
*
* This general function represents the OperatorBase interface for computing. For an operator \hat{O}, this
@ -601,7 +628,7 @@ bool BareKineticEnergy::get(std::ostream& os) const
std::unique_ptr<OperatorBase> BareKineticEnergy::makeClone(ParticleSet& qp, TrialWaveFunction& psi)
{
return std::make_unique<BareKineticEnergy>(*this);
return std::make_unique<BareKineticEnergy>(qp, psi);
}
#ifdef QMC_CUDA

17
src/QMCHamiltonians/BareKineticEnergy.h
View File

@ -73,10 +73,12 @@ public:
* Store mass per species and use SameMass to choose the methods.
* if SameMass, probably faster and easy to vectorize but no impact on the performance.
*/
BareKineticEnergy(ParticleSet& p);
BareKineticEnergy(ParticleSet& p, TrialWaveFunction& psi);
///destructor
~BareKineticEnergy() override;
bool dependsOnWaveFunction() const override { return true; }
std::string getClassName() const override { return "BareKineticEnergy"; }
void resetTargetParticleSet(ParticleSet& P) override {}
#if !defined(REMOVE_TRACEMANAGER)
@ -87,6 +89,17 @@ public:
Return_t evaluate(ParticleSet& P) override;
Return_t evaluateValueAndDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
const Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi) override;
void mw_evaluateWithParameterDerivatives(const RefVectorWithLeader<OperatorBase>& o_list,
const RefVectorWithLeader<ParticleSet>& p_list,
const opt_variables_type& optvars,
const RecordArray<ValueType>& dlogpsi,
RecordArray<ValueType>& dhpsioverpsi) const override;
/** Evaluate the contribution of this component for multiple walkers reporting
* to registered listeners from Estimators.
*/
@ -186,6 +199,8 @@ private:
};
ResourceHandle<MultiWalkerResource> mw_res_;
TrialWaveFunction& psi_;
};
} // namespace qmcplusplus

1
src/QMCHamiltonians/ChiesaCorrection.h
View File

@ -28,6 +28,7 @@ private:
public:
ChiesaCorrection(ParticleSet& ptcl, const TrialWaveFunction& psi) : psi_ref(psi), ptcl_ref(ptcl) {}
std::string getClassName() const override { return "ChiesaCorrection"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;

2
src/QMCHamiltonians/ConservedEnergy.h
View File

@ -76,6 +76,8 @@ struct ConservedEnergy : public OperatorBase
void resetTargetParticleSet(ParticleSet& P) override {}
std::string getClassName() const override { return "ConservedEnergy"; }
Return_t evaluate(ParticleSet& P) override
{
RealType gradsq = Dot(P.G, P.G);

2
src/QMCHamiltonians/CoulombPBCAA.h
View File

@ -92,6 +92,8 @@ struct CoulombPBCAA : public OperatorBase, public ForceBase
~CoulombPBCAA() override;
std::string getClassName() const override { return "CoulombPBCAA"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;

1
src/QMCHamiltonians/CoulombPBCAB.h
View File

@ -134,6 +134,7 @@ public:
void resetTargetParticleSet(ParticleSet& P) override;
std::string getClassName() const override { return "CoulombPBCAB"; }
#if !defined(REMOVE_TRACEMANAGER)
void contributeParticleQuantities() override;

2
src/QMCHamiltonians/CoulombPotential.h
View File

@ -106,6 +106,8 @@ struct CoulombPotential : public OperatorBase, public ForceBase
nCenters = s.getTotalNum();
}
std::string getClassName() const override { return "CoulombPotential"; }
#if !defined(REMOVE_TRACEMANAGER)
void contributeParticleQuantities() override { request_.contribute_array(name_); }

1
src/QMCHamiltonians/DensityEstimator.h
View File

@ -30,6 +30,7 @@ class DensityEstimator : public OperatorBase
public:
DensityEstimator(ParticleSet& elns);
int potentialIndex;
std::string getClassName() const override { return "DensityEstimator"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;

2
src/QMCHamiltonians/DensityMatrices1B.h
View File

@ -159,6 +159,8 @@ public:
DensityMatrices1B(DensityMatrices1B& master, ParticleSet& P, TrialWaveFunction& psi);
~DensityMatrices1B() override;
bool dependsOnWaveFunction() const override { return true; }
std::string getClassName() const override { return "DensityMatrices1B"; }
//standard interface
std::unique_ptr<OperatorBase> makeClone(ParticleSet& P, TrialWaveFunction& psi) final;
bool put(xmlNodePtr cur) override;

1
src/QMCHamiltonians/EnergyDensityEstimator.h
View File

@ -31,6 +31,7 @@ public:
EnergyDensityEstimator(const PSPool& PSP, const std::string& defaultKE);
~EnergyDensityEstimator() override;
std::string getClassName() const override { return "EnergyDensityEstimator"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;
void addObservables(PropertySetType& plist) {}

1
src/QMCHamiltonians/ForceBase.h
View File

@ -88,6 +88,7 @@ private:
public:
BareForce(ParticleSet& ions, ParticleSet& elns);
std::string getClassName() const override { return "BareForce"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;

2
src/QMCHamiltonians/ForceCeperley.cpp
View File

@ -76,7 +76,7 @@ void ForceCeperley::InitMatrix()
// in Numerics/DeterminantOperators.h
invert_matrix(Sinv, false);
// in Numerics/MatrixOperators.h
MatrixOperators::product(Sinv, h.data(), c.data());
MatrixOperators::product(Sinv, h, c);
}
ForceCeperley::Return_t ForceCeperley::evaluate(ParticleSet& P)

2
src/QMCHamiltonians/ForceCeperley.h
View File

@ -41,6 +41,8 @@ public:
ForceCeperley(ParticleSet& ions, ParticleSet& elns);
std::string getClassName() const override { return "ForceCeperley"; }
Return_t evaluate(ParticleSet& P) override;
void InitMatrix();

2
src/QMCHamiltonians/ForceChiesaPBCAA.cpp
View File

@ -68,7 +68,7 @@ void ForceChiesaPBCAA::InitMatrix()
// in Numerics/DeterminantOperators.h
invert_matrix(Sinv, false);
// in Numerics/MatrixOperators.h
MatrixOperators::product(Sinv, h.data(), c.data());
MatrixOperators::product(Sinv, h, c);
}
void ForceChiesaPBCAA::initBreakup(ParticleSet& P)

2
src/QMCHamiltonians/ForceChiesaPBCAA.h
View File

@ -61,6 +61,8 @@ struct ForceChiesaPBCAA : public OperatorBase, public ForceBase
ForceChiesaPBCAA(ParticleSet& ions, ParticleSet& elns, bool firsttime = true);
std::string getClassName() const override { return "ForceChiesaPBCAA"; }
Return_t evaluate(ParticleSet& P) override;
void InitMatrix();

1
src/QMCHamiltonians/ForwardWalking.h
View File

@ -39,6 +39,7 @@ struct ForwardWalking : public OperatorBase
///destructor
~ForwardWalking() override {}
std::string getClassName() const override { return "ForwardWalking"; }
void resetTargetParticleSet(ParticleSet& P) override {}
inline Return_t rejectedMove(ParticleSet& P) override

1
src/QMCHamiltonians/GridExternalPotential.h
View File

@ -42,6 +42,7 @@ struct GridExternalPotential : public OperatorBase
oneBodyQuantumDomain(P);
}
std::string getClassName() const override { return "GridExternalPotential"; }
//unneeded interface functions
void resetTargetParticleSet(ParticleSet& P) override {}

7
src/QMCHamiltonians/HamiltonianFactory.cpp
View File

@ -76,7 +76,6 @@ HamiltonianFactory::HamiltonianFactory(const std::string& hName,
ClassName = "HamiltonianFactory";
myName = hName;
targetPtcl.set_quantum();
targetH->addOperator(std::make_unique<BareKineticEnergy>(targetPtcl), "Kinetic");
}
/** main hamiltonian build function
@ -133,6 +132,9 @@ bool HamiltonianFactory::build(xmlNodePtr cur, bool buildtree)
if (psi_it == psiPool.end())
APP_ABORT("Unknown psi \"" + psiName + "\" for target Psi");
TrialWaveFunction* targetPsi = psi_it->second.get();
// KineticEnergy must be the first element in the hamiltonian array.
if (defaultKE != "no")
targetH->addOperator(std::make_unique<BareKineticEnergy>(targetPtcl, *targetPsi), "Kinetic");
xmlNodePtr cur_saved(cur);
cur = cur->children;
while (cur != NULL)
@ -307,8 +309,7 @@ bool HamiltonianFactory::build(xmlNodePtr cur, bool buildtree)
{
APP_ABORT("Unknown source \"" + source + "\" for DensityMatrices1B");
}
std::unique_ptr<DensityMatrices1B> apot =
std::make_unique<DensityMatrices1B>(targetPtcl, *targetPsi, Pc);
std::unique_ptr<DensityMatrices1B> apot = std::make_unique<DensityMatrices1B>(targetPtcl, *targetPsi, Pc);
apot->put(cur);
targetH->addOperator(std::move(apot), potName, false);
}

1
src/QMCHamiltonians/HarmonicExternalPotential.h
View File

@ -42,6 +42,7 @@ struct HarmonicExternalPotential : public OperatorBase
~HarmonicExternalPotential() override {}
std::string getClassName() const override { return "HarmonicExternalPotential"; }
//unneeded interface functions
void resetTargetParticleSet(ParticleSet& P) override {}

3
src/QMCHamiltonians/L2Potential.h
View File

@ -71,6 +71,9 @@ struct L2Potential : public OperatorBase
L2Potential(const ParticleSet& ions, ParticleSet& els, TrialWaveFunction& psi);
bool dependsOnWaveFunction() const override { return true; }
std::string getClassName() const override { return "L2Potential"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;

1
src/QMCHamiltonians/LatticeDeviationEstimator.h
View File

@ -42,6 +42,7 @@ public:
LatticeDeviationEstimator(ParticleSet& P, ParticleSet& sP, const std::string& tgroup, const std::string& sgroup);
~LatticeDeviationEstimator() override {}
std::string getClassName() const override { return "LatticeDeviationEstimator"; }
bool put(xmlNodePtr cur) override; // read input xml node, required
bool get(std::ostream& os) const override; // class description, required

1
src/QMCHamiltonians/LocalECPotential.h
View File

@ -63,6 +63,7 @@ struct LocalECPotential : public OperatorBase
LocalECPotential(const ParticleSet& ions, ParticleSet& els);
std::string getClassName() const override { return "LocalECPotential"; }
void resetTargetParticleSet(ParticleSet& P) override;
#if !defined(REMOVE_TRACEMANAGER)

1
src/QMCHamiltonians/MPC.h
View File

@ -76,6 +76,7 @@ public:
~MPC() override;
std::string getClassName() const override { return "MPC"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;

2
src/QMCHamiltonians/MomentumEstimator.h
View File

@ -21,6 +21,8 @@ class MomentumEstimator : public OperatorBase
{
public:
MomentumEstimator(ParticleSet& elns, TrialWaveFunction& psi);
bool dependsOnWaveFunction() const override { return true; }
std::string getClassName() const override { return "MomentumEstimator"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;

6
src/QMCHamiltonians/NonLocalECPComponent.h
View File

@ -94,7 +94,7 @@ private:
/// scratch spaces used by evaluateValueAndDerivatives
Matrix<ValueType> dratio;
std::vector<ValueType> dlogpsi_vp;
Vector<ValueType> dlogpsi_vp;
// For Pulay correction to the force
std::vector<RealType> WarpNorm;
@ -241,8 +241,8 @@ public:
RealType r,
const PosType& dr,
const opt_variables_type& optvars,
const std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi);
const Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi);
/**
* @brief Evaluate contribution to B of election iel and ion iat. Filippi scheme for computing fast derivatives.

13
src/QMCHamiltonians/NonLocalECPotential.deriv.cpp
View File

@ -21,13 +21,18 @@ namespace qmcplusplus
{
NonLocalECPotential::Return_t NonLocalECPotential::evaluateValueAndDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
const std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi)
const Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi)
{
value_ = 0.0;
for (int ipp = 0; ipp < PPset.size(); ipp++)
if (PPset[ipp])
PPset[ipp]->randomize_grid(*myRNG);
/* evaluating TWF ratio values requires calling prepareGroup
* In evaluate() we first loop over species and call prepareGroup before looping over all the electrons of a species
* Here it is not necessary because TWF::evaluateLog has been called and precomputed data is up-to-date
*/
const auto& myTable = P.getDistTableAB(myTableIndex);
for (int jel = 0; jel < P.getTotalNum(); jel++)
{
@ -60,8 +65,8 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateValueAndDerivatives
RealType r,
const PosType& dr,
const opt_variables_type& optvars,
const std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi)
const Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi)
{
const size_t num_vars = optvars.num_active_vars;
dratio.resize(nknot, num_vars);

6
src/QMCHamiltonians/NonLocalECPotential.h
View File

@ -37,6 +37,8 @@ public:
NonLocalECPotential(ParticleSet& ions, ParticleSet& els, TrialWaveFunction& psi, bool computeForces, bool enable_DLA);
~NonLocalECPotential() override;
bool dependsOnWaveFunction() const override { return true; }
std::string getClassName() const override { return "NonLocalECPotential"; }
void resetTargetParticleSet(ParticleSet& P) override;
#if !defined(REMOVE_TRACEMANAGER)
@ -98,8 +100,8 @@ public:
Return_t evaluateValueAndDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
const std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi) override;
const Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi) override;
/** Do nothing */
bool put(xmlNodePtr cur) override { return true; }

15
src/QMCHamiltonians/OperatorBase.cpp
View File

@ -126,12 +126,12 @@ void OperatorBase::mw_evaluatePerParticle(const RefVectorWithLeader<OperatorBase
void OperatorBase::mw_evaluateWithParameterDerivatives(const RefVectorWithLeader<OperatorBase>& o_list,
const RefVectorWithLeader<ParticleSet>& p_list,
const opt_variables_type& optvars,
RecordArray<ValueType>& dlogpsi,
const RecordArray<ValueType>& dlogpsi,
RecordArray<ValueType>& dhpsioverpsi) const
{
const int nparam = dlogpsi.nparam();
std::vector<ValueType> tmp_dlogpsi(nparam);
std::vector<ValueType> tmp_dhpsioverpsi(nparam);
Vector<ValueType> tmp_dlogpsi(nparam);
Vector<ValueType> tmp_dhpsioverpsi(nparam);
for (int iw = 0; iw < o_list.size(); iw++)
{
for (int j = 0; j < nparam; j++)
@ -177,9 +177,14 @@ void OperatorBase::mw_evaluatePerParticleWithToperator(
OperatorBase::Return_t OperatorBase::evaluateValueAndDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
const std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi)
const Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi)
{
if (dependsOnWaveFunction())
throw std::logic_error("Bug!! " + getClassName() +
"::evaluateValueAndDerivatives"
"must be overloaded when the OperatorBase depends on a wavefunction.");
return evaluate(P);
}

12
src/QMCHamiltonians/OperatorBase.h
View File

@ -120,6 +120,9 @@ public:
virtual ~OperatorBase() = default;
/// return true if this operator depends on a wavefunction
virtual bool dependsOnWaveFunction() const { return false; }
//////// GETTER AND SETTER FUNCTIONS ////////////////
/**
@ -143,6 +146,9 @@ public:
*/
std::string getName() const noexcept;
/// return class name
virtual std::string getClassName() const = 0;
/**
* @brief Set my_name member, uses small string optimization (pass by value)
*
@ -276,7 +282,7 @@ public:
virtual void mw_evaluateWithParameterDerivatives(const RefVectorWithLeader<OperatorBase>& o_list,
const RefVectorWithLeader<ParticleSet>& p_list,
const opt_variables_type& optvars,
RecordArray<ValueType>& dlogpsi,
const RecordArray<ValueType>& dlogpsi,
RecordArray<ValueType>& dhpsioverpsi) const;
/**
@ -333,8 +339,8 @@ public:
*/
virtual Return_t evaluateValueAndDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
const std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi);
const Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi);
/**
* @brief Evaluate contribution to local energy and derivatives w.r.t ionic coordinates from OperatorBase.

2
src/QMCHamiltonians/OrbitalImages.h
View File

@ -217,6 +217,8 @@ public:
OrbitalImages(ParticleSet& P, const PSPool& PSP, Communicate* mpicomm, const SPOMap& spomap);
OrbitalImages(const OrbitalImages& other);
std::string getClassName() const override { return "OrbitalImages"; }
//standard interface
std::unique_ptr<OperatorBase> makeClone(ParticleSet& P, TrialWaveFunction& psi) final;

1
src/QMCHamiltonians/PairCorrEstimator.h
View File

@ -35,6 +35,7 @@ public:
*/
PairCorrEstimator(ParticleSet& elns, std::string& sources);
std::string getClassName() const override { return "PairCorrEstimator"; }
void resetTargetParticleSet(ParticleSet& P) override;
/* evaluate the pair correlation functions */

2
src/QMCHamiltonians/Pressure.h
View File

@ -48,6 +48,8 @@ struct Pressure : public OperatorBase
///destructor
~Pressure() override {}
bool dependsOnWaveFunction() const override { return true; }
std::string getClassName() const override { return "Pressure"; }
void resetTargetParticleSet(ParticleSet& P) override { pNorm = 1.0 / (P.getLattice().DIM * P.getLattice().Volume); }
inline Return_t evaluate(ParticleSet& P) override

64
src/QMCHamiltonians/QMCHamiltonian.cpp
View File

@ -21,6 +21,7 @@
#include "QMCWaveFunctions/TrialWaveFunction.h"
#include "QMCHamiltonians/NonLocalECPotential.h"
#include "Utilities/TimerManager.h"
#include "BareKineticEnergy.h"
#include "Containers/MinimalContainers/RecordArray.hpp"
#ifdef QMC_CUDA
#include "Particle/MCWalkerConfiguration.h"
@ -29,6 +30,7 @@
namespace qmcplusplus
{
/** constructor
*/
QMCHamiltonian::QMCHamiltonian(const std::string& aname)
@ -37,7 +39,8 @@ QMCHamiltonian::QMCHamiltonian(const std::string& aname)
myName(aname),
nlpp_ptr(nullptr),
l2_ptr(nullptr),
ham_timer_(*timer_manager.createTimer("Hamiltonian:" + aname, timer_level_medium))
ham_timer_(*timer_manager.createTimer("Hamiltonian:" + aname + "::evaluate", timer_level_medium)),
eval_vals_derivs_timer_(*timer_manager.createTimer("Hamiltonian:" + aname + "::ValueParamDerivs", timer_level_medium))
#if !defined(REMOVE_TRACEMANAGER)
,
streaming_position(false),
@ -609,21 +612,29 @@ std::vector<QMCHamiltonian::FullPrecRealType> QMCHamiltonian::mw_evaluate(
QMCHamiltonian::FullPrecRealType QMCHamiltonian::evaluateValueAndDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi,
bool compute_deriv)
Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi)
{
LocalEnergy = KineticEnergy = H[0]->evaluate(P);
if (compute_deriv)
for (int i = 1; i < H.size(); ++i)
LocalEnergy += H[i]->evaluateValueAndDerivatives(P, optvars, dlogpsi, dhpsioverpsi);
else
for (int i = 1; i < H.size(); ++i)
LocalEnergy += H[i]->evaluate(P);
// The first componennt must be BareKineticEnergy for both handling KineticEnergy and dlogpsi computation
// by calling TWF::evaluateDerivatives inside BareKineticEnergy::evaluateValueAndDerivatives
assert(dynamic_cast<BareKineticEnergy*>(H[0].get()) &&
"BUG: The first componennt in Hamiltonian must be BareKineticEnergy.");
ScopedTimer local_timer(eval_vals_derivs_timer_);
{
ScopedTimer h_timer(my_timers_[0]);
LocalEnergy = KineticEnergy = H[0]->evaluateValueAndDerivatives(P, optvars, dlogpsi, dhpsioverpsi);
}
for (int i = 1; i < H.size(); ++i)
{
ScopedTimer h_timer(my_timers_[i]);
LocalEnergy += H[i]->evaluateValueAndDerivatives(P, optvars, dlogpsi, dhpsioverpsi);
}
return LocalEnergy;
}
std::vector<QMCHamiltonian::FullPrecRealType> QMCHamiltonian::mw_evaluateValueAndDerivativesInner(
std::vector<QMCHamiltonian::FullPrecRealType> QMCHamiltonian::mw_evaluateValueAndDerivatives(
const RefVectorWithLeader<QMCHamiltonian>& ham_list,
const RefVectorWithLeader<TrialWaveFunction>& wf_list,
const RefVectorWithLeader<ParticleSet>& p_list,
@ -663,26 +674,6 @@ std::vector<QMCHamiltonian::FullPrecRealType> QMCHamiltonian::mw_evaluateValueAn
return local_energies;
}
std::vector<QMCHamiltonian::FullPrecRealType> QMCHamiltonian::mw_evaluateValueAndDerivatives(
const RefVectorWithLeader<QMCHamiltonian>& ham_list,
const RefVectorWithLeader<TrialWaveFunction>& wf_list,
const RefVectorWithLeader<ParticleSet>& p_list,
const opt_variables_type& optvars,
RecordArray<ValueType>& dlogpsi,
RecordArray<ValueType>& dhpsioverpsi,
bool compute_deriv)
{
std::vector<FullPrecRealType> local_energies(ham_list.size(), 0.0);
if (compute_deriv)
local_energies =
QMCHamiltonian::mw_evaluateValueAndDerivativesInner(ham_list, wf_list, p_list, optvars, dlogpsi, dhpsioverpsi);
else
local_energies = QMCHamiltonian::mw_evaluate(ham_list, wf_list, p_list);
return local_energies;
}
QMCHamiltonian::FullPrecRealType QMCHamiltonian::evaluateVariableEnergy(ParticleSet& P, bool free_nlpp)
{
RealType nlpp = 0.0;
@ -927,6 +918,15 @@ OperatorBase* QMCHamiltonian::getHamiltonian(const std::string& aname)
return nullptr;
}
RefVector<OperatorBase> QMCHamiltonian::getTWFDependentComponents()
{
RefVector<OperatorBase> components;
for (int i = 0; i < H.size(); i++)
if (H[i]->dependsOnWaveFunction())
components.push_back(*H[i]);
return components;
}
void QMCHamiltonian::resetTargetParticleSet(ParticleSet& P)
{
for (int i = 0; i < H.size(); i++)

20
src/QMCHamiltonians/QMCHamiltonian.h
View File

@ -92,6 +92,10 @@ public:
*/
OperatorBase* getHamiltonian(int i) { return H[i].get(); }
/** return components, auxH not included, depending on TWF.
*/
RefVector<OperatorBase> getTWFDependentComponents();
#if !defined(REMOVE_TRACEMANAGER)
///initialize trace data
void initialize_traces(TraceManager& tm, ParticleSet& P);
@ -269,20 +273,10 @@ public:
*/
FullPrecRealType evaluateValueAndDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi,
bool compute_deriv);
Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi);
static std::vector<QMCHamiltonian::FullPrecRealType> mw_evaluateValueAndDerivatives(
const RefVectorWithLeader<QMCHamiltonian>& ham_list,
const RefVectorWithLeader<TrialWaveFunction>& wf_list,
const RefVectorWithLeader<ParticleSet>& p_list,
const opt_variables_type& optvars,
RecordArray<ValueType>& dlogpsi,
RecordArray<ValueType>& dhpsioverpsi,
bool compute_deriv);
static std::vector<QMCHamiltonian::FullPrecRealType> mw_evaluateValueAndDerivativesInner(
const RefVectorWithLeader<QMCHamiltonian>& ham_list,
const RefVectorWithLeader<TrialWaveFunction>& wf_list,
const RefVectorWithLeader<ParticleSet>& p_list,
@ -455,6 +449,8 @@ private:
std::vector<std::unique_ptr<OperatorBase>> auxH;
/// Total timer for H evaluation
NewTimer& ham_timer_;
/// Total timer for H evaluation
NewTimer& eval_vals_derivs_timer_;
/// timers for H components
TimerList_t my_timers_;
///types of component operators

2
src/QMCHamiltonians/SOECPotential.h
View File

@ -22,6 +22,8 @@ class SOECPotential : public OperatorBase
public:
SOECPotential(ParticleSet& ions, ParticleSet& els, TrialWaveFunction& psi);
bool dependsOnWaveFunction() const override { return true; }
std::string getClassName() const override { return "SOECPotential"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;

1
src/QMCHamiltonians/SkAllEstimator.h
View File

@ -31,6 +31,7 @@ class SkAllEstimator : public OperatorBase
public:
SkAllEstimator(ParticleSet& ions, ParticleSet& elns);
std::string getClassName() const override { return "SkAllEstimator"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;

1
src/QMCHamiltonians/SkEstimator.h
View File

@ -30,6 +30,7 @@ class SkEstimator : public OperatorBase
public:
SkEstimator(ParticleSet& elns);
std::string getClassName() const override { return "SkEstimator"; }
void resetTargetParticleSet(ParticleSet& P) override;
Return_t evaluate(ParticleSet& P) override;

1
src/QMCHamiltonians/SkPot.h
View File

@ -29,6 +29,7 @@ class SkPot : public OperatorBase
public:
SkPot(ParticleSet& elns);
std::string getClassName() const override { return "SkPot"; }
void resetTargetParticleSet(ParticleSet& P) override;
[[noreturn]] Return_t evaluate(ParticleSet& P) override;

1
src/QMCHamiltonians/SpeciesKineticEnergy.h
View File

@ -28,6 +28,7 @@ class SpeciesKineticEnergy : public OperatorBase
public:
SpeciesKineticEnergy(ParticleSet& P);
std::string getClassName() const override { return "SpeciesKineticEnergy"; }
bool put(xmlNodePtr cur) override; // read input xml node, required
bool get(std::ostream& os) const override; // class description, required

1
src/QMCHamiltonians/SpinDensity.h
View File

@ -43,6 +43,7 @@ public:
~SpinDensity() override {}
//standard interface
std::string getClassName() const override { return "SpinDensity"; }
std::unique_ptr<OperatorBase> makeClone(ParticleSet& P, TrialWaveFunction& psi) final;
bool put(xmlNodePtr cur) override;
Return_t evaluate(ParticleSet& P) override;

1
src/QMCHamiltonians/StaticStructureFactor.h
View File

@ -37,6 +37,7 @@ public:
~StaticStructureFactor() override {}
//standard interface
std::string getClassName() const override { return "StaticStructureFactor"; }
std::unique_ptr<OperatorBase> makeClone(ParticleSet& P, TrialWaveFunction& psi) final;
bool put(xmlNodePtr cur) override;
Return_t evaluate(ParticleSet& P) override;

2
src/QMCHamiltonians/StressPBC.h
View File

@ -65,6 +65,8 @@ struct StressPBC : public OperatorBase, public ForceBase
bool firstTimeStress;
StressPBC(ParticleSet& ions, ParticleSet& elns, TrialWaveFunction& Psi);
std::string getClassName() const override { return "StressPBC"; }
Return_t evaluate(ParticleSet& P) override;
void initBreakup(ParticleSet& P);

14
src/QMCHamiltonians/tests/test_bare_kinetic.cpp
View File

@ -73,7 +73,8 @@ TEST_CASE("Bare Kinetic Energy", "[hamiltonian]")
xmlNodePtr h1 = xmlFirstElementChild(root);
BareKineticEnergy bare_ke(elec);
TrialWaveFunction psi;
BareKineticEnergy bare_ke(elec, psi);
bare_ke.put(h1);
elec.L[0] = 1.0;
@ -211,7 +212,7 @@ TEST_CASE("Bare KE Pulay PBC", "[hamiltonian]")
xmlNodePtr h1 = xmlFirstElementChild(root);
BareKineticEnergy bare_ke(elec);
BareKineticEnergy bare_ke(elec, psi);
bare_ke.put(h1);
// update all distance tables
@ -298,11 +299,13 @@ void testElecCase(double mass_up,
elec2.R[1][2] = 1.0;
elec2.update();
TrialWaveFunction psi, psi_clone;
RefVector<ParticleSet> ptcls{elec, elec2};
RefVectorWithLeader<ParticleSet> p_list(elec, ptcls);
BareKineticEnergy bare_ke(elec);
BareKineticEnergy bare_ke2(elec);
BareKineticEnergy bare_ke(elec, psi);
BareKineticEnergy bare_ke2(elec, psi_clone);
RefVector<OperatorBase> bare_kes{bare_ke, bare_ke2};
RefVectorWithLeader<OperatorBase> o_list(bare_ke, bare_kes);
@ -311,9 +314,6 @@ void testElecCase(double mass_up,
elec2.L[0] = 1.0;
elec2.L[1] = 0.0;
TrialWaveFunction psi;
TrialWaveFunction psi_clone;
RefVectorWithLeader<TrialWaveFunction> twf_list(psi, {psi, psi_clone});
ResourceCollection bare_ke_res("test_bare_ke_res");

4
src/QMCHamiltonians/tests/test_ecp.cpp
View File

@ -310,8 +310,8 @@ TEST_CASE("Evaluate_ecp", "[hamiltonian]")
}
opt_variables_type optvars;
std::vector<ValueType> dlogpsi;
std::vector<ValueType> dhpsioverpsi;
Vector<ValueType> dlogpsi;
Vector<ValueType> dhpsioverpsi;
psi.checkInVariables(optvars);
optvars.resetIndex();

2
src/QMCTools/PyscfToQmcpack.py
View File

@ -570,7 +570,7 @@ def savetoqmcpack(cell,mf,title="Default",kpts=[],kmesh=[],sp_twist=[],weight=1.
H5_qmcpack.close()
print ('Wavefunction successfully saved to QMCPACK HDF5 Format')
print ('Use: "convert4qmc -pyscf {}.h5" to generate QMCPACK input files'.format(title))
print ('Use: "convert4qmc -orbitals {}.h5" to generate QMCPACK input files'.format(title))
# Close the file before exiting

6
src/QMCTools/QMCGaussianParserBase.cpp
View File

@ -16,8 +16,6 @@
#include "QMCGaussianParserBase.h"
#include "ParticleIO/XMLParticleIO.h"
#include "Numerics/HDFSTLAttrib.h"
#include <iterator>
#include <algorithm>
#include <numeric>
@ -27,7 +25,9 @@
#include <sstream>
#include <bitset>
#include <iomanip>
#include "ParticleIO/XMLParticleIO.h"
#include "Numerics/HDFSTLAttrib.h"
#include "ModernStringUtils.hpp"
//std::vector<std::string> QMCGaussianParserBase::IonName;
const int OhmmsAsciiParser::bufferSize;

4
src/QMCWaveFunctions/AGPDeterminant.cpp
View File

@ -21,9 +21,7 @@ namespace qmcplusplus
{
using std::copy;
AGPDeterminant::AGPDeterminant(BasisSetType* bs)
: WaveFunctionComponent("AGPDeterminant"), GeminalBasis(bs), NumPtcls(0)
{}
AGPDeterminant::AGPDeterminant(BasisSetType* bs) : GeminalBasis(bs), NumPtcls(0) {}
AGPDeterminant::~AGPDeterminant() {}
void AGPDeterminant::resize(int nup, int ndown)

6
src/QMCWaveFunctions/AGPDeterminant.h
View File

@ -46,6 +46,8 @@ public:
///default destructor
~AGPDeterminant() override;
std::string getClassName() const override { return "AGPDeterminant"; }
///reset the size: with the number of particles and number of orbtials
void resize(int nup, int ndown);
@ -93,8 +95,8 @@ public:
void evaluateDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi) override
Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi) override
{}
///Total number of particles

2
src/QMCWaveFunctions/BsplineFactory/SplineC2COMPTarget.cpp
View File

@ -14,7 +14,7 @@
#include "spline2/MultiBsplineEval.hpp"
#include "spline2/MultiBsplineEval_OMPoffload.hpp"
#include "QMCWaveFunctions/BsplineFactory/contraction_helper.hpp"
#include "Platforms/OMPTarget/ompReduction.hpp"
#include "Platforms/OMPTarget/ompReductionComplex.hpp"
#include "OMPTarget/OMPTargetMath.hpp"
namespace qmcplusplus

1
src/QMCWaveFunctions/BsplineFactory/SplineR2R.h
View File

@ -79,6 +79,7 @@ public:
virtual std::string getClassName() const override { return "SplineR2R"; }
virtual std::string getKeyword() const override { return "SplineR2R"; }
bool isComplex() const override { return false; };
bool isRotationSupported() const override { return true; }
std::unique_ptr<SPOSet> makeClone() const override { return std::make_unique<SplineR2R>(*this); }

7
src/QMCWaveFunctions/CompositeSPOSet.cpp
View File

@ -121,13 +121,6 @@ void CompositeSPOSet::evaluate(const ParticleSet& P, PosType& r, ValueVector& ps
}
#endif
//methods to be implemented later
void CompositeSPOSet::resetParameters(const opt_variables_type& optVariables)
{
for (int c = 0; c < components.size(); ++c)
components[c]->resetParameters(optVariables);
}
void CompositeSPOSet::evaluate_notranspose(const ParticleSet& P,
int first,
int last,

2
src/QMCWaveFunctions/CompositeSPOSet.h
View File

@ -64,8 +64,6 @@ public:
APP_ABORT("CompositeSPOSet::" + method + " has not been implemented");
}
void resetParameters(const opt_variables_type& optVariables) override;
//methods to be implemented in the future (possibly)
#ifdef QMC_CUDA
void evaluate(const ParticleSet& P, PosType& r, ValueVector& psi) override;

8
src/QMCWaveFunctions/ConstantOrbital.h
View File

@ -22,7 +22,9 @@ class ConstantOrbital : public WaveFunctionComponent
public:
PsiValueType FakeGradRatio;
ConstantOrbital() : WaveFunctionComponent("ConstantOrbital"), FakeGradRatio(1.0) {}
ConstantOrbital() : FakeGradRatio(1.0) {}
std::string getClassName() const override { return "ConstantOrbital"; }
LogValueType evaluateLog(const ParticleSet& P,
ParticleSet::ParticleGradient& G,
@ -56,8 +58,8 @@ public:
void evaluateDerivatives(ParticleSet& P,
const opt_variables_type& optvars,
std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi) override
Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi) override
{}
};

6
src/QMCWaveFunctions/ExampleHeComponent.cpp
View File

@ -209,7 +209,7 @@ std::unique_ptr<WaveFunctionComponent> ExampleHeComponent::makeClone(ParticleSet
return std::make_unique<ExampleHeComponent>(*this);
}
void ExampleHeComponent::resetParameters(const OptVariablesType& active)
void ExampleHeComponent::resetParametersExclusive(const OptVariablesType& active)
{
if (my_vars_.size())
{
@ -227,8 +227,8 @@ void ExampleHeComponent::resetParameters(const OptVariablesType& active)
void ExampleHeComponent::evaluateDerivatives(ParticleSet& P,
const OptVariablesType& optvars,
std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi)
Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi)
{
using RealGradType = TinyVector<RealType, 3>;

17
src/QMCWaveFunctions/ExampleHeComponent.h
View File

@ -23,11 +23,11 @@
*/
namespace qmcplusplus
{
class ExampleHeComponent : public WaveFunctionComponent
class ExampleHeComponent : public WaveFunctionComponent, OptimizableObject
{
public:
ExampleHeComponent(const ParticleSet& ions, ParticleSet& els)
: WaveFunctionComponent("ExampleHeComponent"),
: OptimizableObject("example"),
ions_(ions),
my_table_ee_idx_(els.addTable(els, DTModes::NEED_TEMP_DATA_ON_HOST | DTModes::NEED_VP_FULL_TABLE_ON_HOST)),
my_table_ei_idx_(els.addTable(ions, DTModes::NEED_VP_FULL_TABLE_ON_HOST)){};
@ -35,13 +35,12 @@ public:
using OptVariablesType = optimize::VariableSet;
using PtclGrpIndexes = QMCTraits::PtclGrpIndexes;
std::string getClassName() const override { return "ExampleHeComponent"; }
void extractOptimizableObjectRefs(UniqueOptObjRefs& opt_obj_refs) override { opt_obj_refs.push_back(*this); }
bool isOptimizable() const override { return true; }
void checkInVariables(OptVariablesType& active) override { active.insertFrom(my_vars_); }
void checkInVariablesExclusive(OptVariablesType& active) override { active.insertFrom(my_vars_); }
void checkOutVariables(const OptVariablesType& active) override { my_vars_.getIndex(active); }
void resetParameters(const OptVariablesType& active) override;
void reportStatus(std::ostream& os) override {}
void resetParametersExclusive(const OptVariablesType& active) override;
LogValueType evaluateLog(const ParticleSet& P,
ParticleSet::ParticleGradient& G,
@ -59,8 +58,8 @@ public:
void evaluateDerivatives(ParticleSet& P,
const OptVariablesType& optvars,
std::vector<ValueType>& dlogpsi,
std::vector<ValueType>& dhpsioverpsi) override;
Vector<ValueType>& dlogpsi,
Vector<ValueType>& dhpsioverpsi) override;
void registerData(ParticleSet& P, WFBufferType& buf) override {}

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

@ -247,7 +247,7 @@ std::unique_ptr<BackflowFunctionBase> BackflowBuilder::addOneBody(xmlNodePtr cur
for (int i = 0; i < funs.size(); i++)
{
// BsplineFunctor<RealType> *bsp = new BsplineFunctor<RealType>(cusps[i]);
auto bsp = std::make_unique<BsplineFunctor<RealType>>();
auto bsp = std::make_unique<BsplineFunctor<RealType>>(extractCoefficientsID(funs[i]));
bsp->cutoff_radius = targetPtcl.getLattice().WignerSeitzRadius;
bsp->put(funs[i]);
if (bsp->cutoff_radius > cutOff)
@ -323,7 +323,7 @@ std::unique_ptr<BackflowFunctionBase> BackflowBuilder::addTwoBody(xmlNodePtr cur
}
app_log() << "Adding radial component for species: " << spA << " " << spB << " " << ia << " " << ib
<< std::endl;
auto bsp = std::make_unique<BsplineFunctor<RealType>>();
auto bsp = std::make_unique<BsplineFunctor<RealType>>(extractCoefficientsID(cur));
bsp->cutoff_radius = targetPtcl.getLattice().WignerSeitzRadius;
bsp->put(cur);
if (bsp->cutoff_radius > cutOff)
@ -625,7 +625,7 @@ void BackflowBuilder::makeShortRange_twoBody(xmlNodePtr cur,
{
APP_ABORT("Unknown correlation type " + type + " in Backflow.");
}
auto bsp = std::make_unique<BsplineFunctor<RealType>>();
auto bsp = std::make_unique<BsplineFunctor<RealType>>(extractCoefficientsID(cur));
if (init == "true" || init == "yes")
{
app_log() << "Initializing backflow radial functions with RPA.";
@ -670,7 +670,7 @@ void BackflowBuilder::makeShortRange_twoBody(xmlNodePtr cur,
// for (int j=0; j<nfitgaussians; j++) y1_c+=gb[j]*(std::exp(-x[1]*x[1]/((j+1)*Rcut*Rcut)));
// for (int j=0; j<nfitgaussians; j++) y2_c+=gb[j]*(std::exp(-x[2]*x[2]/((j+1)*Rcut*Rcut)));
//make a temp functor to ensure right BC's (Necessary?)
auto tmp_bsp = std::make_unique<BsplineFunctor<RealType>>();
auto tmp_bsp = std::make_unique<BsplineFunctor<RealType>>("tmp_bsp");
tmp_bsp->initialize(12, x, y, cusp, Rcut, id, optimize);
// tmp_bsp->print(app_log());
for (int i = 0; i < myGrid->size(); i++)

2
src/QMCWaveFunctions/Fermion/BackflowTransformation.cpp
View File

@ -21,7 +21,7 @@
namespace qmcplusplus
{
BackflowTransformation::BackflowTransformation(ParticleSet& els)
: QP(els), cutOff(0.0), myTableIndex_(els.addTable(els))
: OptimizableObject("bf"), QP(els), cutOff(0.0), myTableIndex_(els.addTable(els))
{
NumTargets = els.getTotalNum();
Bmat.resize(NumTargets);

Some files were not shown because too many files have changed in this diff Show More