Merge branch 'develop' into test_estimator_py3

CMake/python.cmake

@@ -5,7 +5,7 @@
 #   MODULE_PRESENT - output - True/False based on success of the import
 FUNCTION (TEST_PYTHON_MODULE MODULE_NAME MODULE_PRESENT)
   EXECUTE_PROCESS(
-    COMMAND python ${qmcpack_SOURCE_DIR}/tests/scripts/test_import.py ${MODULE_NAME}
+    COMMAND ${qmcpack_SOURCE_DIR}/tests/scripts/test_import.py ${MODULE_NAME}
     OUTPUT_VARIABLE TMP_OUTPUT_VAR
     OUTPUT_STRIP_TRAILING_WHITESPACE
   )

CMakeLists.txt

@@ -54,6 +54,18 @@ IF (COMPILING_ON_BGQ)
   cmake_policy(SET CMP0060 OLD)
 ENDIF()
 
+######################################################################
+# Verify Python3 available
+######################################################################
+INCLUDE(CMake/python.cmake)
+IF(${CMAKE_VERSION} VERSION_LESS "3.12.0")
+  TEST_PYTHON_MODULE( "os"  Python3_FOUND )
+ELSE()
+  find_package(Python3)
+ENDIF()
+IF ( NOT Python3_FOUND )
+  MESSAGE( FATAL_ERROR "Could not find required python3" )
+ENDIF ( NOT Python3_FOUND )
 ######################################################################
 # Verify QE executables present if QE_BIN specified
 ######################################################################
@@ -66,7 +78,6 @@ ENDIF()
 # Verify PYSCF package is present
 ######################################################################
 IF ( DEFINED PYSCF_BIN )
-  INCLUDE(CMake/python.cmake)
   TEST_PYTHON_MODULE(pyscf HAVE_PYSCF)
   IF (NOT HAVE_PYSCF)
     MESSAGE( FATAL_ERROR "PYSCF_BIN was specified but could not import pyscf." )

README.md

@@ -13,6 +13,7 @@
  * BOOST, peer-reviewed portable C++ source libraries, http://www.boost.org
  * FFTW, FFT library, http://www.fftw.org/
  * MPI, parallel library. Optional, but a near requirement for production calculations.
+ * Python3. Older versions are not supported as of January 2020.
 
 We aim to support open source compilers and libraries released within two years of each QMCPACK release. Use of software versions over
 two years old may work but is discouraged and untested. Proprietary compilers (Intel, PGI) are generally supported over the same

examples/afqmc/01-neon_atom/scf.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 from pyscf import gto, scf, cc
 from pyscf.cc import ccsd_t

examples/afqmc/02-neon_frozen_core/scf.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 from pyscf import gto, scf, cc
 from pyscf.cc import ccsd_t

examples/afqmc/03-carbon_triplet_uhf/scf.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 # Triplet UHF ground state of carbon atom.
 
 from pyscf import gto, scf

examples/afqmc/04-N2_nomsd/scf.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 import sys
 import h5py

examples/afqmc/05-N2_phmsd/scf.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 import sys
 import h5py

examples/afqmc/06-methane_converge_back_prop/scf.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 from math import cos, sin, pi, acos
 import numpy

examples/solids/pyscf-inputs/diamondC_1x1x1_pp_Bspline.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 '''
 Gamma point post-HF calculation needs only real integrals.

examples/solids/pyscf-inputs/diamondC_1x1x1_pp_LCAO.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 '''
 Gamma point post-HF calculation needs only real integrals.

labs/lab1_qmc_statistics/average/average.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 average = ( 
 (-0.45298911858) + 

labs/lab1_qmc_statistics/average/stddev.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 stddev = (1./(5.*(5.-1.)) * ( 
 (-0.45298911858-(-0.464736835668))**2 + 

labs/lab1_qmc_statistics/average/stddev2.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 stddev = ( (1./(5.-1.)) * (
 (0.60984565298-(-0.45298911858)**2) + 

labs/lab2_qmc_basics/oxygen_atom/ip_conv.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 import os
 import sys

labs/lab2_qmc_basics/oxygen_atom/reference/ip_conv.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 import os
 import sys

labs/lab2_qmc_basics/oxygen_dimer/O_dimer.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,get_machine,run_project

labs/lab2_qmc_basics/oxygen_dimer/dimer_fit.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 import sys
 import os

labs/lab2_qmc_basics/oxygen_dimer/reference/O_dimer.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,get_machine,run_project

labs/lab2_qmc_basics/oxygen_dimer/reference/dimer_fit.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 import sys
 import os

labs/lab2_qmc_basics/your_system/example.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,run_project,get_machine

labs/lab2_qmc_basics/your_system/reference/example.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import project suite functions
 from project import settings,Job,run_project,get_machine

labs/lab4_condensed_matter/Be-16at-qmc.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,run_project,get_machine,obj

labs/lab4_condensed_matter/Be-2at-qmc.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,run_project,get_machine,obj

labs/lab4_condensed_matter/Be-2at-setup.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,run_project,get_machine,obj

labs/lab4_condensed_matter/graphene-final.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,run_project,get_machine

labs/lab4_condensed_matter/graphene-loop-mesh.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,run_project,get_machine

labs/lab4_condensed_matter/graphene-setup.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,run_project,get_machine

labs/lab5_excited_states/band.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 from nexus import settings,job,run_project
 from nexus import generate_physical_system

labs/lab5_excited_states/optical.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 from nexus import settings,job,run_project
 from nexus import generate_physical_system

labs/lab5_excited_states/quasiparticle.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 from nexus import settings,job,run_project
 from nexus import generate_physical_system

labs/misc/jobrun_vesta

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 import os
 import sys

manual/convert4qmc.tex

@@ -511,7 +511,7 @@ For the bash shell, this can be done as follows:\\
 Copy and paste the following code in a file named LiH.py.
 
 \begin{lstlisting}[style=Python]
-#!/usr/bin/env python
+#! /usr/bin/env python3
 from pyscf import gto, scf, df
 import numpy
 

manual/gaussian_orbitals_solids.tex

@@ -115,7 +115,7 @@ can be generated according to the following example for a $2 \times 1 \times 1$
 %requiring complex coefficients will be supported in future releases.
 
 \begin{lstlisting}[style=Python,caption=Example PySCF input for single k-point calculation for a $2 \times 1 \times 1$ carbon supercell.]
-#!/usr/bin/env python
+#! /usr/bin/env python3
 import numpy
 import h5py
 from pyscf.pbc import gto, scf, dft, df
@@ -192,7 +192,7 @@ export PYTHONPATH=QMCPACK_PATH/src/QMCTools:$PYTHONPATH
 %The following example corresponds to the same carbon system ($2 \times 1 \times 1$); however, in this case, we use a primitive simulation cell and a $2 \times 1 \times 1$ k-point mesh.   
 
 %\begin{lstlisting}[style=Python,caption=Example PySCF input for single k-point calculation for a $2 \times 1 \times 1$ carbon supercell.]
-%#!/usr/bin/env python
+%#! /usr/bin/env python3
 
 %import numpy
 %from pyscf.pbc import gto, scf, dft,df

manual/lab_qmc_basics.tex

@@ -760,7 +760,7 @@ Nexus is driven by simple user-defined scripts that resemble keyword-driven inpu
 
 \ifws
 \begin{lstlisting}[style=Python]
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,get_machine,run_project 
@@ -821,7 +821,7 @@ run_project(qmc)                       # write input file and submit job
 \else
 
 \begin{lstlisting}[style=Python]
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # import Nexus functions
 from nexus import settings,job,get_machine,run_project 

nexus/docs/code_sections/settings_examples.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 from nexus import settings
 

nexus/docs/nexus_user_guide.tex

@@ -1114,7 +1114,7 @@ The script is similar to the one discussed in section \ref{sec:user_imports}, di
 %\begin{lstlisting}[caption={Nexus input script for the usage example (\texttt{diamond.py}).  DFT is performed in diamond primitive cell followed by VMC in a 16 atom supercell.},label={ex:d16}]
 \HRule
 \begin{minted}{python}
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 # nexus imports
 from nexus import settings,Job,run_project
@@ -1357,7 +1357,7 @@ Take a moment to study the ``input file'' script
 
 \HRule
 \begin{minted}{python}
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 from nexus import settings,Job,run_project
 from nexus import generate_physical_system
@@ -1811,7 +1811,7 @@ rather than being generated), the boundary conditions (open BC's, see
 
 \HRule
 \begin{minted}{python}
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 from nexus import settings,Job,run_project
 from nexus import Structure,PhysicalSystem
@@ -2439,7 +2439,7 @@ In the "qmc" object given below, pay attention to keyword "excitation   = ['up',
 %\begin{lstlisting}[caption={Nexus input script for the usage example (\texttt{diamond.py}).  DFT is performed in diamond primitive cell followed by VMC in a 16 atom supercell.},label={ex:d16}]
 \HRule
 \begin{minted}{python}
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 from nexus import settings,job,run_project
 from nexus import generate_physical_system

nexus/docs/reference_python/reference_python.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 
 from project import settings

nexus/examples/pyscf/01_h2o_hf/README.rst

@@ -18,7 +18,7 @@ calculation is used for this example.  The contents of this script
 
 .. code-block:: python
 
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     from pyscf import scf
     
@@ -37,7 +37,7 @@ The Nexus script, ``h2o_ae_hf.py``, is shown below:
 
 .. code-block:: python
 
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     from nexus import settings,job,run_project,obj
     from nexus import generate_physical_system
@@ -118,7 +118,7 @@ Next, let's look at the PySCF script produced by Nexus (see
 
 .. code-block:: python
 
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     from pyscf import scf
     

nexus/examples/pyscf/02_diamond_hf/README.rst

@@ -12,7 +12,7 @@ shown below:
 
 .. code-block:: python
 
-    #!/usr/bin/env python
+    #! /usr/bin/env python3
     
     from pyscf.pbc import df, scf
     
@@ -33,7 +33,7 @@ The Nexus script, ``diamond_pp_hf_gamma.py``, is shown below:
 
 .. code-block:: python
 
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     from nexus import settings,job,run_project,obj
     from nexus import generate_physical_system
@@ -117,7 +117,7 @@ Next, let's look at the PySCF script produced by Nexus (see
 
 .. code-block:: python
 
-    #!/usr/bin/env python
+    #! /usr/bin/env python3
     
     from pyscf.pbc import df, scf
     

nexus/examples/pyscf/02_diamond_hf/scf_template.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#! /usr/bin/env python3
 
 from pyscf.pbc import df, scf
 

nexus/examples/pyscf/diamond/dft_template.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#! /usr/bin/env python3
 
 from pyscf.pbc import df, scf
 

nexus/examples/qmcpack/afqmc_pyscf/01_neon_atom/scf_template.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#! /usr/bin/env python3
 
 from pyscf import gto, scf, cc
 from pyscf.cc import ccsd_t

nexus/examples/qmcpack/afqmc_pyscf/02_neon_frozen_core/scf_template.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#! /usr/bin/env python3
 
 from pyscf import gto, scf, cc
 from pyscf.cc import ccsd_t

nexus/examples/qmcpack/afqmc_pyscf/03_carbon_triplet_uhf/scf_template.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#! /usr/bin/env python3
 # Triplet UHF ground state of carbon atom.
 
 from pyscf import gto, scf

nexus/examples/qmcpack/afqmc_pyscf/06_methane_back_prop/scf_template.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#! /usr/bin/env python3
 
 from math import cos, sin, pi, acos
 import numpy

nexus/examples/qmcpack/rsqmc_misc/diamond_pyscf/dft_template.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#! /usr/bin/env python3
 
 from pyscf.pbc import df, scf
 

nexus/examples/qmcpack/rsqmc_pyscf/01_h2o_hf_qmc/README.rst

@@ -18,7 +18,7 @@ Important parts of the Nexus script, ``h2o_ae_hf_qmc.py``, that differ from
 Example 1 are shown below:
 
 .. code-block:: python
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     ...
     from nexus import generate_convert4qmc

nexus/examples/qmcpack/rsqmc_pyscf/02_diamond_hf_qmc/README.rst

@@ -12,7 +12,7 @@ The most important changes relative to Example 2 are shown below:
  
 .. code-block:: python
 
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     ...
     from nexus import generate_convert4qmc

nexus/examples/qmcpack/rsqmc_pyscf/02_diamond_hf_qmc/scf_template.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#! /usr/bin/env python3
 
 from pyscf.pbc import df, scf
 

nexus/examples/qmcpack/rsqmc_quantum_espresso/01_diamond_dft_vmc/README.rst

@@ -37,7 +37,7 @@ Important differences from the prior example are bolded.
 
 .. parsed-literal:: 
 
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     from nexus import settings,job,run_project
     from nexus import generate_physical_system

nexus/examples/qmcpack/rsqmc_quantum_package/01_h2o_hf_qmc/README.rst

@@ -18,7 +18,7 @@ Important parts of the Nexus script, ``h2o_ae_hf_qmc.py``, that differ from
 Example 1 are shown below:
 
 .. code-block:: python
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     ...
     from nexus import generate_convert4qmc

nexus/examples/qmcpack/rsqmc_quantum_package/02_o2_selci_qmc/README.rst

@@ -20,7 +20,7 @@ Important parts of the Nexus script, ``o2_selci_vmc_dmc.py``, that differ
 from Example 2 are shown below:
 
 .. code-block:: python
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     ...
     from nexus import generate_convert4qmc

nexus/examples/quantum_package/01_h2o_hf/README.rst

@@ -20,7 +20,7 @@ the script are shown below:
 
 .. code-block:: python
 
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     from nexus import settings,job,run_project
     from nexus import generate_physical_system

nexus/examples/quantum_package/02_o2_selci/README.rst

@@ -16,7 +16,7 @@ the script are shown below:
 
 .. code-block:: python
 
-    #! /usr/bin/env python
+    #! /usr/bin/env python3
     
     from nexus import settings,job,run_project
     from nexus import generate_physical_system

nexus/install

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 import os
 import sys

nexus/tests/unit/test_physical_system.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 import numpy as np
 import testing

nexus/tests/unit/test_pyscf_input.py

@@ -11,7 +11,7 @@ H  0.000000  0.757160 -0.586260
 '''
 
 
-scf_template = '''#! /usr/bin/env python
+scf_template = '''#! /usr/bin/env python3
 
 from pyscf import scf
 
@@ -223,7 +223,7 @@ def test_write():
     text = open(write_path,'r').read()
 
     ref_text = '''
-        #! /usr/bin/env python
+        #! /usr/bin/env python3
         
         from pyscf import scf
         
@@ -297,7 +297,7 @@ def test_write():
     text = open(write_path,'r').read()
 
     ref_text = '''
-        #! /usr/bin/env python
+        #! /usr/bin/env python3
         
         from pyscf import scf
         

src/Particle/VirtualParticleSet.cpp

@@ -38,14 +38,20 @@ VirtualParticleSet::VirtualParticleSet(const ParticleSet& p, int nptcl) : refPS(
 }
 
 /// move virtual particles to new postions and update distance tables
-void VirtualParticleSet::makeMoves(int jel, const ParticlePos_t& vitualPos, bool sphere, int iat)
+void VirtualParticleSet::makeMoves(int jel,
+                                   const PosType& ref_pos,
+                                   const std::vector<PosType>& deltaV,
+                                   bool sphere,
+                                   int iat)
 {
   if (sphere && iat < 0)
     APP_ABORT("VirtualParticleSet::makeMoves is invoked incorrectly, the flag sphere=true requires iat specified!");
   onSphere      = sphere;
   refPtcl       = jel;
   refSourcePtcl = iat;
-  R             = vitualPos;
+  assert(R.size() == deltaV.size());
+  for (size_t ivp = 0; ivp < R.size(); ivp++)
+    R[ivp] = ref_pos + deltaV[ivp];
   update();
 }
 

src/Particle/VirtualParticleSet.h

@@ -49,11 +49,16 @@ public:
 
   /** move virtual particles to new postions and update distance tables
      * @param jel reference particle that all the VP moves from
-     * @param vitualPos new positions
+     * @param ref_pos reference particle position
+     * @param deltaV Position delta for virtual moves.
      * @param sphere set true if VP are on a sphere around the reference source particle
      * @param iat reference source particle
      */
-  void makeMoves(int jel, const ParticlePos_t& vitualPos, bool sphere = false, int iat = -1);
+  void makeMoves(int jel,
+                 const PosType& ref_pos,
+                 const std::vector<PosType>& deltaV,
+                 bool sphere = false,
+                 int iat     = -1);
 };
 } // namespace qmcplusplus
 #endif

src/QMCDrivers/CMakeLists.txt

@@ -59,6 +59,7 @@ SET(QMCDRIVERS
   DMC/DMCDriverInput.cpp
   DMC/WalkerControlFactory.cpp
   DMC/WalkerReconfiguration.cpp
+  DMC/SODMCUpdatePbyPFast.cpp
   RMC/RMC.cpp
   RMC/RMCUpdatePbyP.cpp
   RMC/RMCUpdateAll.cpp

src/QMCDrivers/DMC/DMC.cpp

@@ -19,6 +19,7 @@
 
 #include "QMCDrivers/DMC/DMC.h"
 #include "QMCDrivers/DMC/DMCUpdatePbyP.h"
+#include "QMCDrivers/DMC/SODMCUpdatePbyP.h"
 #include "QMCDrivers/DMC/DMCUpdateAll.h"
 #include "QMCApp/HamiltonianPool.h"
 #include "Message/Communicate.h"
@@ -48,7 +49,9 @@ DMC::DMC(MCWalkerConfiguration& w,
       KillNodeCrossing(0),
       BranchInterval(-1),
       Reconfiguration("no"),
-      mover_MaxAge(-1)
+      mover_MaxAge(-1),
+      SpinMoves("no"),
+      SpinMass(1.0)
 {
   RootName = "dmc";
   QMCType  = "DMC";
@@ -59,6 +62,8 @@ DMC::DMC(MCWalkerConfiguration& w,
   m_param.add(NonLocalMove, "nonlocalmove", "string");
   m_param.add(NonLocalMove, "nonlocalmoves", "string");
   m_param.add(mover_MaxAge, "MaxAge", "double");
+  m_param.add(SpinMoves,"SpinMoves","string");
+  m_param.add(SpinMass,"SpinMass","double");
 }
 
 void DMC::resetUpdateEngines()
@@ -83,6 +88,13 @@ void DMC::resetUpdateEngines()
     estimatorClones.resize(NumThreads, 0);
     traceClones.resize(NumThreads, 0);
     FairDivideLow(W.getActiveWalkers(), NumThreads, wPerNode);
+
+    tolower(SpinMoves);
+    if (SpinMoves != "yes" && SpinMoves != "no")
+    {
+      APP_ABORT("SpinMoves must be yes/no\n");
+    }
+
     {
       //log file
       std::ostringstream o;
@@ -98,6 +110,8 @@ void DMC::resetUpdateEngines()
         o << "\n  Walkers are killed when a node crossing is detected";
       else
         o << "\n  DMC moves are rejected when a node crossing is detected";
+      if (SpinMoves=="yes")
+        o << "\n  Spins treated as dynamic variable with SpinMass: " << SpinMass;
       app_log() << o.str() << std::endl;
     }
 #pragma omp parallel for
@@ -114,22 +128,40 @@ void DMC::resetUpdateEngines()
       Rng[ip] = new RandomGenerator_t(*RandomNumberControl::Children[ip]);
       hClones[ip]->setRandomGenerator(Rng[ip]);
 #endif
-      if (qmc_driver_mode[QMC_UPDATE_MODE])
+      if (SpinMoves == "yes")
       {
-        Movers[ip] = new DMCUpdatePbyPWithRejectionFast(*wClones[ip], *psiClones[ip], *hClones[ip], *Rng[ip]);
-        Movers[ip]->put(qmcNode);
-        Movers[ip]->resetRun(branchEngine, estimatorClones[ip], traceClones[ip], DriftModifier);
-        Movers[ip]->initWalkersForPbyP(W.begin() + wPerNode[ip], W.begin() + wPerNode[ip + 1]);
+        if (qmc_driver_mode[QMC_UPDATE_MODE])
+        {
+          Movers[ip] = new SODMCUpdatePbyPWithRejectionFast(*wClones[ip], *psiClones[ip], *hClones[ip], *Rng[ip]);
+          Movers[ip]->setSpinMass(SpinMass);
+          Movers[ip]->put(qmcNode);
+          Movers[ip]->resetRun(branchEngine, estimatorClones[ip], traceClones[ip], DriftModifier);
+          Movers[ip]->initWalkersForPbyP(W.begin() + wPerNode[ip], W.begin() + wPerNode[ip+1]);
+        }
+        else 
+        {
+          APP_ABORT("SODMC Driver Mode must be PbyP\n");
+        }
       }
-      else
+      else 
       {
-        if (KillNodeCrossing)
-          Movers[ip] = new DMCUpdateAllWithKill(*wClones[ip], *psiClones[ip], *hClones[ip], *Rng[ip]);
+        if (qmc_driver_mode[QMC_UPDATE_MODE])
+        {
+          Movers[ip] = new DMCUpdatePbyPWithRejectionFast(*wClones[ip], *psiClones[ip], *hClones[ip], *Rng[ip]);
+          Movers[ip]->put(qmcNode);
+          Movers[ip]->resetRun(branchEngine, estimatorClones[ip], traceClones[ip], DriftModifier);
+          Movers[ip]->initWalkersForPbyP(W.begin() + wPerNode[ip], W.begin() + wPerNode[ip + 1]);
+        }
         else
-          Movers[ip] = new DMCUpdateAllWithRejection(*wClones[ip], *psiClones[ip], *hClones[ip], *Rng[ip]);
-        Movers[ip]->put(qmcNode);
-        Movers[ip]->resetRun(branchEngine, estimatorClones[ip], traceClones[ip], DriftModifier);
-        Movers[ip]->initWalkers(W.begin() + wPerNode[ip], W.begin() + wPerNode[ip + 1]);
+        {
+          if (KillNodeCrossing)
+            Movers[ip] = new DMCUpdateAllWithKill(*wClones[ip], *psiClones[ip], *hClones[ip], *Rng[ip]);
+          else
+            Movers[ip] = new DMCUpdateAllWithRejection(*wClones[ip], *psiClones[ip], *hClones[ip], *Rng[ip]);
+          Movers[ip]->put(qmcNode);
+          Movers[ip]->resetRun(branchEngine, estimatorClones[ip], traceClones[ip], DriftModifier);
+          Movers[ip]->initWalkers(W.begin() + wPerNode[ip], W.begin() + wPerNode[ip + 1]);
+        }
       }
     }
   }

src/QMCDrivers/DMC/DMC.h

@@ -58,6 +58,9 @@ private:
   std::string UseFastGrad;
   ///input to control maximum age allowed for walkers.
   IndexType mover_MaxAge;
+  ///turn on spin moves
+  std::string SpinMoves;
+  RealType SpinMass;
 
   void resetUpdateEngines();
   /// Copy Constructor (disabled)

src/QMCDrivers/DMC/DMCBatched.cpp

@@ -272,7 +272,7 @@ void DMCBatched::advanceWalkers(const StateForThread& sft,
   }
 
   std::for_each(crowd.get_walker_twfs().begin(), crowd.get_walker_twfs().end(),
-                [](auto& twf) { twf.get().completeUpdates(); });
+                [](TrialWaveFunction& twf) { twf.completeUpdates(); });
 
   ParticleSet::flex_donePbyP(crowd.get_walker_elecs());
   //dmc_timers.dmc_movePbyP.stop();

src/QMCDrivers/DMC/SODMCUpdatePbyP.h

@@ -0,0 +1,54 @@
+//////////////////////////////////////////////////////////////////////////////////////
+// This file is distributed under the University of Illinois/NCSA Open Source License.
+// See LICENSE file in top directory for details.
+//
+// Copyright (c) 2020 Jeongnim Kim and QMCPACK developers.
+//
+// File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
+//                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
+//                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
+//                    Jeongnim Kim, jeongnim.kim@intel.com, Intel Corp.
+//
+// File created by: Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories
+//////////////////////////////////////////////////////////////////////////////////////
+
+
+#ifndef QMCPLUSPLUS_SODMC_UPDATE_PARTICLEBYPARTCLE_H
+#define QMCPLUSPLUS_SODMC_UPDATE_PARTICLEBYPARTCLE_H
+#include "QMCDrivers/QMCUpdateBase.h"
+#include "Utilities/NewTimer.h"
+namespace qmcplusplus
+{
+class SODMCUpdatePbyPWithRejectionFast : public QMCUpdateBase
+{
+public:
+  /// Constructor.
+  SODMCUpdatePbyPWithRejectionFast(MCWalkerConfiguration& w,
+                                   TrialWaveFunction& psi,
+                                   QMCHamiltonian& h,
+                                   RandomGenerator_t& rg);
+  ///destructor
+  ~SODMCUpdatePbyPWithRejectionFast();
+
+  void advanceWalker(Walker_t& thisWalker, bool recompute);
+
+private:
+  TimerList_t myTimers;
+};
+
+
+enum SODMCTimers
+{
+  SODMC_buffer,
+  SODMC_movePbyP,
+  SODMC_hamiltonian,
+  SODMC_collectables,
+  SODMC_tmoves
+};
+
+extern TimerNameList_t<SODMCTimers> SODMCTimerNames;
+
+
+} // namespace qmcplusplus
+
+#endif

src/QMCDrivers/DMC/SODMCUpdatePbyPFast.cpp

@@ -0,0 +1,190 @@
+//////////////////////////////////////////////////////////////////////////////////////
+// This file is distributed under the University of Illinois/NCSA Open Source License.
+// See LICENSE file in top directory for details.
+//
+// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
+//
+// File developed by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
+//                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
+//                    Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
+//                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
+//
+// File created by: Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories
+//////////////////////////////////////////////////////////////////////////////////////
+
+
+#include "QMCDrivers/DMC/SODMCUpdatePbyP.h"
+#include "Particle/MCWalkerConfiguration.h"
+#include "Particle/HDFWalkerIO.h"
+#include "ParticleBase/ParticleUtility.h"
+#include "ParticleBase/RandomSeqGenerator.h"
+#include "QMCDrivers/DriftOperators.h"
+#if !defined(REMOVE_TRACEMANAGER)
+#include "Estimators/TraceManager.h"
+#else
+typedef int TraceManager;
+#endif
+//#define TEST_INNERBRANCH
+
+
+namespace qmcplusplus
+{
+TimerNameList_t<SODMCTimers> SODMCTimerNames = {{SODMC_buffer, "SODMCUpdatePbyP::Buffer"},
+                                                {SODMC_movePbyP, "SODMCUpdatePbyP::movePbyP"},
+                                                {SODMC_hamiltonian, "SODMCUpdatePbyP::Hamiltonian"},
+                                                {SODMC_collectables, "SODMCUpdatePbyP::Collectables"},
+                                                {SODMC_tmoves, "SODMCUpdatePbyP::Tmoves"}};
+
+
+/// Constructor.
+SODMCUpdatePbyPWithRejectionFast::SODMCUpdatePbyPWithRejectionFast(MCWalkerConfiguration& w,
+                                                                   TrialWaveFunction& psi,
+                                                                   QMCHamiltonian& h,
+                                                                   RandomGenerator_t& rg)
+    : QMCUpdateBase(w, psi, h, rg)
+{
+  setup_timers(myTimers, SODMCTimerNames, timer_level_medium);
+}
+
+/// destructor
+SODMCUpdatePbyPWithRejectionFast::~SODMCUpdatePbyPWithRejectionFast() {}
+
+void SODMCUpdatePbyPWithRejectionFast::advanceWalker(Walker_t& thisWalker, bool recompute)
+{
+  myTimers[SODMC_buffer]->start();
+  Walker_t::WFBuffer_t& w_buffer(thisWalker.DataSet);
+  W.loadWalker(thisWalker, true);
+  Psi.copyFromBuffer(W, w_buffer);
+  myTimers[SODMC_buffer]->stop();
+  //create a 3N-Dimensional Gaussian with variance=1
+  makeGaussRandomWithEngine(deltaR, RandomGen);
+  makeGaussRandomWithEngine(deltaS, RandomGen);
+  int nAcceptTemp(0);
+  int nRejectTemp(0);
+  //copy the old energy and scale factor of drift
+  FullPrecRealType eold(thisWalker.Properties(LOCALENERGY));
+  FullPrecRealType enew(eold);
+  RealType rr_proposed = 0.0;
+  RealType rr_accepted = 0.0;
+  RealType gf_acc      = 1.0;
+  myTimers[SODMC_movePbyP]->start();
+  for (int ig = 0; ig < W.groups(); ++ig) //loop over species
+  {
+    RealType tauovermass = Tau * MassInvS[ig];
+    RealType oneover2tau = 0.5 / (tauovermass);
+    RealType sqrttau     = std::sqrt(tauovermass);
+    for (int iat = W.first(ig); iat < W.last(ig); ++iat)
+    {
+      //get the displacement
+      TrialWaveFunction::LogValueType spingrad_iat;
+      GradType grad_iat = Psi.evalGradWithSpin(W, iat, spingrad_iat);
+      PosType dr;
+      ParticleSet::Scalar_t ds;
+      DriftModifier->getDrift(tauovermass, grad_iat, dr);
+      ds = tauovermass/spinMass*std::real(spingrad_iat); //using raw spin grad, no UNR modifier
+      dr += sqrttau * deltaR[iat];
+      ds += std::sqrt(tauovermass/spinMass)*deltaS[iat];
+      RealType rr = tauovermass * dot(deltaR[iat], deltaR[iat]);
+      rr_proposed += rr;
+      if (rr > m_r2max)
+      {
+        ++nRejectTemp;
+        continue;
+      }
+      if (!W.makeMoveAndCheckWithSpin(iat, dr, ds))
+        continue;
+      ValueType ratio = Psi.calcRatioGradWithSpin(W, iat, grad_iat, spingrad_iat);
+      //node is crossed reject the move
+      if (branchEngine->phaseChanged(Psi.getPhaseDiff()))
+      {
+        ++nRejectTemp;
+        ++nNodeCrossing;
+        W.rejectMove(iat);
+        Psi.rejectMove(iat);
+      }
+      else
+      {
+        FullPrecRealType logGf = -0.5 * dot(deltaR[iat], deltaR[iat]);
+        logGf += -0.5*deltaS[iat]*deltaS[iat];
+        //Use the force of the particle iat
+        DriftModifier->getDrift(tauovermass, grad_iat, dr);
+        ds = tauovermass/spinMass*std::real(spingrad_iat);
+        dr                     = W.R[iat] - W.activePos - dr;
+        ds                     = W.spins[iat] - W.activeSpinVal - ds;
+        FullPrecRealType logGb = -oneover2tau * dot(dr, dr);
+        logGb                 += -spinMass*oneover2tau*ds*ds;
+        RealType prob          = std::norm(ratio) * std::exp(logGb - logGf);
+        if (RandomGen() < prob)
+        {
+          ++nAcceptTemp;
+          Psi.acceptMove(W, iat, true);
+          W.acceptMove(iat, true);
+          rr_accepted += rr;
+          gf_acc *= prob; //accumulate the ratio
+        }
+        else
+        {
+          ++nRejectTemp;
+          W.rejectMove(iat);
+          Psi.rejectMove(iat);
+        }
+      }
+    }
+  }
+  Psi.completeUpdates();
+  W.donePbyP();
+  myTimers[SODMC_movePbyP]->stop();
+
+  if (nAcceptTemp > 0)
+  {
+    //need to overwrite the walker properties
+    myTimers[SODMC_buffer]->start();
+    thisWalker.Age  = 0;
+    RealType logpsi = Psi.updateBuffer(W, w_buffer, recompute);
+    W.saveWalker(thisWalker);
+    myTimers[SODMC_buffer]->stop();
+    myTimers[SODMC_hamiltonian]->start();
+    enew = H.evaluateWithToperator(W);
+    myTimers[SODMC_hamiltonian]->stop();
+    thisWalker.resetProperty(logpsi, Psi.getPhase(), enew, rr_accepted, rr_proposed, 1.0);
+    thisWalker.Weight *= branchEngine->branchWeight(enew, eold);
+    myTimers[SODMC_collectables]->start();
+    H.auxHevaluate(W, thisWalker);
+    H.saveProperty(thisWalker.getPropertyBase());
+    myTimers[SODMC_collectables]->stop();
+  }
+  else
+  {
+    //all moves are rejected: does not happen normally with reasonable wavefunctions
+    thisWalker.Age++;
+    thisWalker.Properties(R2ACCEPTED) = 0.0;
+    //weight is set to 0 for traces
+    // consistent w/ no evaluate/auxHevaluate
+    RealType wtmp     = thisWalker.Weight;
+    thisWalker.Weight = 0.0;
+    H.rejectedMove(W, thisWalker);
+    thisWalker.Weight = wtmp;
+    ++nAllRejected;
+    enew   = eold; //copy back old energy
+    gf_acc = 1.0;
+    thisWalker.Weight *= branchEngine->branchWeight(enew, eold);
+  }
+#if !defined(REMOVE_TRACEMANAGER)
+  Traces->buffer_sample(W.current_step);
+#endif
+  myTimers[SODMC_tmoves]->start();
+  const int NonLocalMoveAcceptedTemp = H.makeNonLocalMoves(W);
+  if (NonLocalMoveAcceptedTemp > 0)
+  {
+    RealType logpsi = Psi.updateBuffer(W, w_buffer, false);
+    W.saveWalker(thisWalker);
+    NonLocalMoveAccepted += NonLocalMoveAcceptedTemp;
+  }
+  myTimers[SODMC_tmoves]->stop();
+  nAccept += nAcceptTemp;
+  nReject += nRejectTemp;
+
+  setMultiplicity(thisWalker);
+}
+
+} // namespace qmcplusplus

src/QMCDrivers/VMC/SOVMCUpdatePbyP.h

@@ -2,7 +2,7 @@
 // This file is distributed under the University of Illinois/NCSA Open Source License.
 // See LICENSE file in top directory for details.
 //
-// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
+// Copyright (c) 2020 Jeongnim Kim and QMCPACK developers.
 //
 // File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
 //                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign

src/QMCDrivers/VMC/VMC.cpp

@@ -182,7 +182,12 @@ void VMC::resetRun()
       Rng[ip] = new RandomGenerator_t(*(RandomNumberControl::Children[ip]));
 #endif
       hClones[ip]->setRandomGenerator(Rng[ip]);
-
+      
+      tolower(SpinMoves);
+      if (SpinMoves != "yes" && SpinMoves != "no") 
+      {
+        APP_ABORT("SpinMoves must be yes/no\n");
+      }
       if (SpinMoves == "yes") 
       {
         if (qmc_driver_mode[QMC_UPDATE_MODE]) 

src/QMCDrivers/tests/test_dmc_driver.cpp

@@ -141,4 +141,104 @@ TEST_CASE("DMC", "[drivers][dmc]")
   delete doc;
   delete p_bke;
 }
+
+TEST_CASE("SODMC", "[drivers][dmc]")
+{
+  Communicate* c;
+  OHMMS::Controller->initialize(0, NULL);
+  c = OHMMS::Controller;
+
+  ParticleSet ions;
+  MCWalkerConfiguration elec;
+
+  ions.setName("ion");
+  ions.create(1);
+  ions.R[0][0] = 0.0;
+  ions.R[0][1] = 0.0;
+  ions.R[0][2] = 0.0;
+
+  elec.setName("elec");
+  std::vector<int> agroup(1);
+  agroup[0] = 1;
+  elec.create(agroup);
+  elec.R[0][0] = 1.0;
+  elec.R[0][1] = 0.0;
+  elec.R[0][2] = 0.0;
+  elec.spins[0] = 0.0;
+  elec.createWalkers(1);
+
+  SpeciesSet& tspecies         = elec.getSpeciesSet();
+  int upIdx                    = tspecies.addSpecies("u");
+  int chargeIdx                = tspecies.addAttribute("charge");
+  int massIdx                  = tspecies.addAttribute("mass");
+  tspecies(chargeIdx, upIdx)   = -1;
+  tspecies(massIdx, upIdx)     = 1.0;
+
+#ifdef ENABLE_SOA
+  elec.addTable(ions, DT_SOA);
+#else
+  elec.addTable(ions, DT_AOS);
+#endif
+  elec.update();
+
+  CloneManager::clear_for_unit_tests();
+
+  TrialWaveFunction psi(c);
+  ConstantOrbital* orb  = new ConstantOrbital;
+  psi.addComponent(orb, "Constant");
+  psi.registerData(elec, elec.WalkerList[0]->DataSet);
+  elec.WalkerList[0]->DataSet.allocate();
+
+  FakeRandom rg;
+
+  QMCHamiltonian h;
+  BareKineticEnergy<double>* p_bke = new BareKineticEnergy<double>(elec);
+  h.addOperator(p_bke, "Kinetic");
+  h.addObservables(elec); // get double free error on 'h.Observables' w/o this
+
+  elec.resetWalkerProperty(); // get memory corruption w/o this
+
+  HamiltonianPool hpool(c);
+
+  WaveFunctionPool wpool(c);
+
+  //EstimatorManagerBase emb(c);
+
+
+  DMC dmc_omp(elec, psi, h, wpool, c);
+
+  const char* dmc_input = "<qmc method=\"dmc\"> \
+   <parameter name=\"steps\">1</parameter> \
+   <parameter name=\"blocks\">1</parameter> \
+   <parameter name=\"timestep\">0.1</parameter> \
+   <parameter name=\"SpinMoves\">yes</parameter> \
+   <parameter name=\"SpinMass\">0.25</parameter> \
+  </qmc> \
+  ";
+  Libxml2Document* doc  = new Libxml2Document;
+  bool okay             = doc->parseFromString(dmc_input);
+  REQUIRE(okay);
+  xmlNodePtr root = doc->getRoot();
+
+  dmc_omp.process(root); // need to call 'process' for QMCDriver, which in turn calls 'put'
+
+  dmc_omp.run();
+
+  // With the constant wavefunction, no moves should be rejected
+  double ar = dmc_omp.acceptRatio();
+  REQUIRE(ar == Approx(1.0));
+
+  // Each electron moved sqrt(tau)*gaussian_rng()
+  //  See Particle>Base/tests/test_random_seq.cpp for the gaussian random numbers
+  //  Values from diffuse.py for moving one step
+
+  REQUIRE(elec[0]->R[0][0] == Approx(0.627670258894097));
+  REQUIRE(elec.R[0][1] == Approx(0.0));
+  REQUIRE(elec.R[0][2] == Approx(-0.372329741105903));
+
+  REQUIRE(elec.spins[0] == Approx(-0.74465948215809097));
+
+  delete doc;
+  delete p_bke;
+}
 } // namespace qmcplusplus

src/QMCHamiltonians/ECPotentialBuilder.cpp

@@ -79,8 +79,6 @@ bool ECPotentialBuilder::put(xmlNodePtr cur)
   bool doForces = (forces == "yes") || (forces == "true");
   if (use_DLA == "yes")
     app_log() << "    Using determinant localization approximation (DLA)" << std::endl;
-  if (NLPP_algo == "batched" && use_DLA == "yes")
-    myComm->barrier_and_abort("Batched pseudopotential evaluation has not been made available to DLA!");
   if (ecpFormat == "xml")
   {
     useXmlFormat(cur);
@@ -131,7 +129,7 @@ bool ECPotentialBuilder::put(xmlNodePtr cur)
     NonLocalECPotential_CUDA* apot = new NonLocalECPotential_CUDA(IonConfig, targetPtcl, targetPsi, usePBC, doForces);
 #else
     NonLocalECPotential* apot =
-        new NonLocalECPotential(IonConfig, targetPtcl, targetPsi, doForces);
+        new NonLocalECPotential(IonConfig, targetPtcl, targetPsi, doForces, use_DLA == "yes");
 #endif
     int nknot_max = 0;
     for (int i = 0; i < nonLocalPot.size(); i++)
@@ -141,8 +139,6 @@ bool ECPotentialBuilder::put(xmlNodePtr cur)
         nknot_max = std::max(nknot_max, nonLocalPot[i]->getNknot());
         if (NLPP_algo == "batched")
           nonLocalPot[i]->initVirtualParticle(targetPtcl);
-        if (use_DLA == "yes")
-          nonLocalPot[i]->enableDLA();
         apot->addComponent(i, nonLocalPot[i]);
       }
     }

src/QMCHamiltonians/NonLocalECPComponent.cpp

@@ -19,13 +19,7 @@
 
 namespace qmcplusplus
 {
-NonLocalECPComponent::NonLocalECPComponent()
-    : lmax(0), nchannel(0), nknot(0), Rmax(-1), VP(nullptr), use_DLA(false)
-{
-#if !defined(REMOVE_TRACEMANAGER)
-  streaming_particles = false;
-#endif
-}
+NonLocalECPComponent::NonLocalECPComponent() : lmax(0), nchannel(0), nknot(0), Rmax(-1), VP(nullptr) {}
 
 NonLocalECPComponent::~NonLocalECPComponent()
 {
@@ -67,6 +61,7 @@ void NonLocalECPComponent::resize_warrays(int n, int m, int l)
   deltaV.resize(n);
   cosgrad.resize(n);
   wfngrad.resize(n);
+  knot_pots.resize(n);
   vrad.resize(m);
   dvrad.resize(m);
   vgrad.resize(m);
@@ -114,48 +109,51 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOne(ParticleSet& W,
                                                                  int iel,
                                                                  RealType r,
                                                                  const PosType& dr,
-                                                                 bool Tmove,
-                                                                 std::vector<NonLocalData>& Txy)
+                                                                 bool use_DLA)
 {
-  constexpr RealType czero(0);
-  constexpr RealType cone(1);
+  buildQuadraturePointDeltaPositions(r, dr, deltaV);
 
   if (VP)
   {
     // Compute ratios with VP
-    ParticleSet::ParticlePos_t VPos(nknot);
-    for (int j = 0; j < nknot; j++)
-    {
-      deltaV[j] = r * rrotsgrid_m[j] - dr;
-      VPos[j]   = deltaV[j] + W.R[iel];
-    }
-    VP->makeMoves(iel, VPos, true, iat);
-    psi.evaluateRatios(*VP, psiratio);
-    for (int j = 0; j < nknot; j++)
-      psiratio[j] *= sgridweight_m[j];
+    VP->makeMoves(iel, W.R[iel], deltaV, true, iat);
+    if (use_DLA)
+      psi.evaluateRatios(*VP, psiratio, TrialWaveFunction::ComputeType::FERMIONIC);
+    else
+      psi.evaluateRatios(*VP, psiratio);
   }
   else
   {
     // Compute ratio of wave functions
     for (int j = 0; j < nknot; j++)
     {
-      deltaV[j] = r * rrotsgrid_m[j] - dr;
       W.makeMove(iel, deltaV[j], false);
-      if(use_DLA)
-        psiratio[j] = psi.calcRatio(W, iel, TrialWaveFunction::ComputeType::FERMIONIC) * sgridweight_m[j];
+      if (use_DLA)
+        psiratio[j] = psi.calcRatio(W, iel, TrialWaveFunction::ComputeType::FERMIONIC);
       else
-        psiratio[j] = psi.calcRatio(W, iel) * sgridweight_m[j];
+        psiratio[j] = psi.calcRatio(W, iel);
       W.rejectMove(iel);
       psi.resetPhaseDiff();
     }
   }
 
+  return calculateProjector(r, dr);
+}
+
+NonLocalECPComponent::RealType NonLocalECPComponent::calculateProjector(RealType r, const PosType& dr)
+{
+  for (int j = 0; j < nknot; j++)
+    psiratio[j] *= sgridweight_m[j];
+
   // Compute radial potential, multiplied by (2l+1) factor.
   for (int ip = 0; ip < nchannel; ip++)
     vrad[ip] = nlpp_m[ip]->splint(r) * wgt_angpp_m[ip];
 
+  constexpr RealType czero(0);
+  constexpr RealType cone(1);
+
   const RealType rinv = cone / r;
-  RealType pairpot    = 0;
+  RealType pairpot    = czero;
   // Compute spherical harmonics on grid
   for (int j = 0; j < nknot; j++)
   {
@@ -176,35 +174,119 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOne(ParticleSet& W,
     ValueType lsum = 0.0;
     for (int l = 0; l < nchannel; l++)
       lsum += vrad[l] * lpol[angpp_m[l]];
-    lsum *= psiratio[j];
-    if (Tmove)
-      Txy.push_back(NonLocalData(iel, std::real(lsum), deltaV[j]));
-    pairpot += std::real(lsum);
+    knot_pots[j] = std::real(lsum * psiratio[j]);
+    pairpot += knot_pots[j];
   }
 
-#if !defined(REMOVE_TRACEMANAGER)
-  if (streaming_particles)
-  {
-    (*Vi_sample)(iat) += .5 * pairpot;
-    (*Ve_sample)(iel) += .5 * pairpot;
-  }
-#endif
   return pairpot;
 }
 
+void NonLocalECPComponent::flex_evaluateOne(const RefVector<NonLocalECPComponent>& ecp_component_list,
+                                            const RefVector<ParticleSet>& p_list,
+                                            const std::vector<int>& iat_list,
+                                            const RefVector<TrialWaveFunction>& psi_list,
+                                            const std::vector<int>& iel_list,
+                                            const std::vector<RealType>& r_list,
+                                            const std::vector<PosType>& dr_list,
+                                            std::vector<RealType>& pairpots,
+                                            bool use_DLA)
+{
+  if (ecp_component_list.size() > 1)
+  {
+    if (ecp_component_list[0].get().VP)
+    {
+      // Compute ratios with VP
+#pragma omp parallel for
+      for (size_t i = 0; i < ecp_component_list.size(); i++)
+      {
+        NonLocalECPComponent& component(ecp_component_list[i]);
+        ParticleSet& W(p_list[i]);
+        int iat(iat_list[i]);
+        int iel(iel_list[i]);
+        auto r   = r_list[i];
+        auto& dr = dr_list[i];
+
+        component.buildQuadraturePointDeltaPositions(r, dr, component.deltaV);
+
+        component.VP->makeMoves(iel, W.R[iel], component.deltaV, true, iat);
+      }
+
+      RefVector<const VirtualParticleSet> vp_list;
+      RefVector<std::vector<ValueType>> psiratios_list;
+      vp_list.reserve(ecp_component_list.size());
+      psiratios_list.reserve(ecp_component_list.size());
+
+      for (size_t i = 0; i < ecp_component_list.size(); i++)
+      {
+        NonLocalECPComponent& component(ecp_component_list[i]);
+        vp_list.push_back(*component.VP);
+        psiratios_list.push_back(component.psiratio);
+      }
+
+      if (use_DLA)
+        TrialWaveFunction::flex_evaluateRatios(psi_list, vp_list, psiratios_list,
+                                               TrialWaveFunction::ComputeType::FERMIONIC);
+      else
+        TrialWaveFunction::flex_evaluateRatios(psi_list, vp_list, psiratios_list);
+    }
+    else
+    {
+      // Compute ratios without VP. This is working but very slow code path.
+#pragma omp parallel for
+      for (size_t i = 0; i < p_list.size(); i++)
+      {
+        NonLocalECPComponent& component(ecp_component_list[i]);
+        auto* VP = component.VP;
+        ParticleSet& W(p_list[i]);
+        int iat(iat_list[i]);
+        TrialWaveFunction& psi(psi_list[i]);
+        int iel(iel_list[i]);
+        auto r   = r_list[i];
+        auto& dr = dr_list[i];
+
+        component.buildQuadraturePointDeltaPositions(r, dr, component.deltaV);
+
+        // Compute ratio of wave functions
+        for (int j = 0; j < component.getNknot(); j++)
+        {
+          W.makeMove(iel, component.deltaV[j], false);
+          if (use_DLA)
+            component.psiratio[j] = psi.calcRatio(W, iel, TrialWaveFunction::ComputeType::FERMIONIC);
+          else
+            component.psiratio[j] = psi.calcRatio(W, iel);
+          W.rejectMove(iel);
+          psi.resetPhaseDiff();
+        }
+      }
+    }
+
+    for (size_t i = 0; i < p_list.size(); i++)
+    {
+      NonLocalECPComponent& component(ecp_component_list[i]);
+      auto r      = r_list[i];
+      auto& dr    = dr_list[i];
+      pairpots[i] = component.calculateProjector(r, dr);
+    }
+  }
+  else if (ecp_component_list.size() == 1)
+    pairpots[0] = ecp_component_list[0].get().evaluateOne(p_list[0], iat_list[0], psi_list[0], iel_list[0], r_list[0],
+                                                          dr_list[0], use_DLA);
+}
+
 NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOneWithForces(ParticleSet& W,
                                                                            int iat,
                                                                            TrialWaveFunction& psi,
                                                                            int iel,
                                                                            RealType r,
                                                                            const PosType& dr,
-                                                                           PosType& force_iat,
-                                                                           bool Tmove,
-                                                                           std::vector<NonLocalData>& Txy)
+                                                                           PosType& force_iat)
 {
   constexpr RealType czero(0);
   constexpr RealType cone(1);
 
+  for (int j = 0; j < nknot; j++)
+    deltaV[j] = r * rrotsgrid_m[j] - dr;
+
   GradType gradtmp_(0);
   PosType realgradtmp_(0);
 
@@ -220,30 +302,20 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOneWithForces(Parti
   {
     APP_ABORT("NonLocalECPComponent::evaluateOneWithForces(...): Forces not implemented with virtual particle moves\n");
     // Compute ratios with VP
-    ParticleSet::ParticlePos_t VPos(nknot);
-    for (int j = 0; j < nknot; j++)
-    {
-      deltaV[j] = r * rrotsgrid_m[j] - dr;
-      VPos[j]   = deltaV[j] + W.R[iel];
-    }
-    VP->makeMoves(iel, VPos, true, iat);
+    VP->makeMoves(iel, W.R[iel], deltaV, true, iat);
     psi.evaluateRatios(*VP, psiratio);
-    for (int j = 0; j < nknot; j++)
-      psiratio[j] *= sgridweight_m[j];
   }
   else
   {
     // Compute ratio of wave functions
     for (int j = 0; j < nknot; j++)
     {
-      deltaV[j] = r * rrotsgrid_m[j] - dr;
       W.makeMove(iel, deltaV[j], false);
-      ValueType ratio = psi.calcRatioGrad(W, iel, gradtmp_);
-      psiratio[j] = ratio * sgridweight_m[j];
+      psiratio[j] = psi.calcRatioGrad(W, iel, gradtmp_);
       //QMCPACK spits out $\nabla\Psi(q)/\Psi(q)$.
       //Multiply times $\Psi(q)/\Psi(r)$ to get
       // $\nabla\Psi(q)/\Psi(r)
-      gradtmp_ *= ratio;
+      gradtmp_ *= psiratio[j];
 #if defined(QMC_COMPLEX)
       //And now we take the real part and save it.
       convert(gradtmp_, gradpsiratio[j]);
@@ -258,6 +330,9 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOneWithForces(Parti
     }
   }
 
+  for (int j = 0; j < nknot; j++)
+    psiratio[j] *= sgridweight_m[j];
+
   // This is just a temporary variable to dump d2/dr2 into for spline evaluation.
   RealType secondderiv(0);
 
@@ -314,25 +389,16 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOneWithForces(Parti
 
     for (int l = 0; l < nchannel; l++)
     {
-      lsum           += std::real(vrad[l]) * lpol[angpp_m[l]] * std::real(psiratio[j]);
-      gradpotterm_   += vgrad[l] * lpol[angpp_m[l]] * std::real(psiratio[j]);
+      lsum += std::real(vrad[l]) * lpol[angpp_m[l]];
+      gradpotterm_ += vgrad[l] * lpol[angpp_m[l]] * std::real(psiratio[j]);
       gradlpolyterm_ += std::real(vrad[l]) * dlpol[angpp_m[l]] * cosgrad[j] * std::real(psiratio[j]);
-      gradwfnterm_   += std::real(vrad[l]) * lpol[angpp_m[l]] * wfngrad[j];
+      gradwfnterm_ += std::real(vrad[l]) * lpol[angpp_m[l]] * wfngrad[j];
     }
-
-    if (Tmove)
-      Txy.push_back(NonLocalData(iel, lsum, deltaV[j]));
-    pairpot += lsum;
+    knot_pots[j] = std::real(lsum * psiratio[j]);
+    pairpot += knot_pots[j];
     force_iat += gradpotterm_ + gradlpolyterm_ - gradwfnterm_;
   }
 
-#if !defined(REMOVE_TRACEMANAGER)
-  if (streaming_particles)
-  {
-    (*Vi_sample)(iat) += .5 * pairpot;
-    (*Ve_sample)(iel) += .5 * pairpot;
-  }
-#endif
   return pairpot;
 }
 
@@ -344,13 +410,14 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOneWithForces(Parti
                                                                            RealType r,
                                                                            const PosType& dr,
                                                                            PosType& force_iat,
-                                                                           ParticleSet::ParticlePos_t& pulay_terms,
-                                                                           bool Tmove,
-                                                                           std::vector<NonLocalData>& Txy)
+                                                                           ParticleSet::ParticlePos_t& pulay_terms)
 {
   constexpr RealType czero(0);
   constexpr RealType cone(1);
 
+  for (int j = 0; j < nknot; j++)
+    deltaV[j] = r * rrotsgrid_m[j] - dr;
+
   GradType gradtmp_(0);
   PosType realgradtmp_(0);
 
@@ -363,8 +430,6 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOneWithForces(Parti
   PosType gradwfnterm_(0);
 
   //Now for the Pulay specific stuff...
-  //Full (potentiall complex) $\Psi(...q...)/\Psi(...r...)$ for all quadrature points q.
-  std::vector<ValueType> psiratiofull_(nknot);
   // $\nabla_I \Psi(...r...)/\Psi(...r...)$
   ParticleSet::ParticlePos_t pulay_ref;
   ParticleSet::ParticlePos_t pulaytmp_;
@@ -376,38 +441,27 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOneWithForces(Parti
   //resize everything.
   pulay_ref.resize(ions.getTotalNum());
   pulaytmp_.resize(ions.getTotalNum());
-  for (unsigned int j = 0; j < nknot; j++)
+  for (size_t j = 0; j < nknot; j++)
     pulay_quad[j].resize(ions.getTotalNum());
 
   if (VP)
   {
     APP_ABORT("NonLocalECPComponent::evaluateOneWithForces(...): Forces not implemented with virtual particle moves\n");
     // Compute ratios with VP
-    ParticleSet::ParticlePos_t VPos(nknot);
-    for (int j = 0; j < nknot; j++)
-    {
-      deltaV[j] = r * rrotsgrid_m[j] - dr;
-      VPos[j]   = deltaV[j] + W.R[iel];
-    }
-    VP->makeMoves(iel, VPos, true, iat);
+    VP->makeMoves(iel, W.R[iel], deltaV, true, iat);
     psi.evaluateRatios(*VP, psiratio);
-    for (int j = 0; j < nknot; j++)
-      psiratio[j] *= sgridweight_m[j];
   }
   else
   {
     // Compute ratio of wave functions
     for (int j = 0; j < nknot; j++)
     {
-      deltaV[j] = r * rrotsgrid_m[j] - dr;
       W.makeMove(iel, deltaV[j], false);
-      ValueType ratio = psi.calcRatioGrad(W, iel, gradtmp_);
-      psiratiofull_[j] = ratio;
-      psiratio[j] = ratio * sgridweight_m[j];
+      psiratio[j] = psi.calcRatioGrad(W, iel, gradtmp_);
       //QMCPACK spits out $\nabla\Psi(q)/\Psi(q)$.
       //Multiply times $\Psi(q)/\Psi(r)$ to get
       // $\nabla\Psi(q)/\Psi(r)
-      gradtmp_ *= ratio;
+      gradtmp_ *= psiratio[j];
 #if defined(QMC_COMPLEX)
       //And now we take the real part and save it.
       convert(gradtmp_, gradpsiratio[j]);
@@ -422,6 +476,9 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOneWithForces(Parti
     }
   }
 
+  for (int j = 0; j < nknot; j++)
+    psiratio[j] *= sgridweight_m[j];
+
   // This is just a temporary variable to dump d2/dr2 into for spline evaluation.
   RealType secondderiv(0);
 
@@ -437,24 +494,24 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOneWithForces(Parti
 
   //Now to construct the 3N dimensional ionic wfn derivatives for pulay terms.
   //This is going to be slow an painful for now.
-  for (unsigned int jat = 0; jat < ions.getTotalNum(); jat++)
+  for (size_t jat = 0; jat < ions.getTotalNum(); jat++)
   {
     pulay_ref[jat] = psi.evalGradSource(W, ions, jat);
     gradpotterm_   = 0;
-    for (unsigned int j = 0; j < nknot; j++)
+    for (size_t j = 0; j < nknot; j++)
     {
       deltaV[j] = r * rrotsgrid_m[j] - dr;
-      //This sequence is necessary to update the distance tables and make the 
-      //inverse matrix available for force computation.  Move the particle to 
+      //This sequence is necessary to update the distance tables and make the
+      //inverse matrix available for force computation.  Move the particle to
       //quadrature point...
       W.makeMove(iel, deltaV[j]);
       psi.calcRatio(W, iel);
       psi.acceptMove(W, iel);
       W.acceptMove(iel); // it only updates the jel-th row of e-e table
-      //Done with the move.  Ready for force computation.  
-    
+      //Done with the move.  Ready for force computation.
+
       iongradtmp_ = psi.evalGradSource(W, ions, jat);
-      iongradtmp_ *= psiratiofull_[j] * sgridweight_m[j];
+      iongradtmp_ *= psiratio[j];
 #ifdef QMC_COMPLEX
       convert(iongradtmp_, pulay_quad[j][jat]);
 #endif
@@ -516,29 +573,22 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateOneWithForces(Parti
     {
       //Note.  Because we are computing "forces", there's a -1 difference between this and
       //direct finite difference calculations.
-      lsum           += std::real(vrad[l]) * lpol[angpp_m[l]] * std::real(psiratio[j]);
-      gradpotterm_   += vgrad[l] * lpol[angpp_m[l]] * std::real(psiratio[j]);
+      lsum += std::real(vrad[l]) * lpol[angpp_m[l]];
+      gradpotterm_ += vgrad[l] * lpol[angpp_m[l]] * std::real(psiratio[j]);
       gradlpolyterm_ += std::real(vrad[l]) * dlpol[angpp_m[l]] * cosgrad[j] * std::real(psiratio[j]);
-      gradwfnterm_   += std::real(vrad[l]) * lpol[angpp_m[l]] * wfngrad[j];
+      gradwfnterm_ += std::real(vrad[l]) * lpol[angpp_m[l]] * wfngrad[j];
       pulaytmp_ -= std::real(vrad[l]) * lpol[angpp_m[l]] * pulay_quad[j];
     }
-    pulaytmp_ += lsum * pulay_ref;
-    if (Tmove)
-      Txy.push_back(NonLocalData(iel, lsum, deltaV[j]));
-    pairpot += lsum;
+    knot_pots[j] = std::real(lsum * psiratio[j]);
+    pulaytmp_ += knot_pots[j] * pulay_ref;
+    pairpot += knot_pots[j];
     force_iat += gradpotterm_ + gradlpolyterm_ - gradwfnterm_;
     pulay_terms += pulaytmp_;
   }
 
-#if !defined(REMOVE_TRACEMANAGER)
-  if (streaming_particles)
-  {
-    (*Vi_sample)(iat) += .5 * pairpot;
-    (*Ve_sample)(iel) += .5 * pairpot;
-  }
-#endif
   return pairpot;
 }
+
 ///Randomly rotate sgrid_m
 void NonLocalECPComponent::randomize_grid(RandomGenerator_t& myRNG)
 {
@@ -574,6 +624,20 @@ void NonLocalECPComponent::randomize_grid(std::vector<T>& sphere, RandomGenerato
       sphere[OHMMS_DIM * i + j] = rrotsgrid_m[i][j];
 }
 
+void NonLocalECPComponent::buildQuadraturePointDeltaPositions(RealType r,
+                                                              const PosType& dr,
+                                                              std::vector<PosType>& deltaV) const
+{
+  for (int j = 0; j < nknot; j++)
+    deltaV[j] = r * rrotsgrid_m[j] - dr;
+}
+
+void NonLocalECPComponent::contributeTxy(int iel, std::vector<NonLocalData>& Txy) const
+{
+  for (int j = 0; j < nknot; j++)
+    Txy.push_back(NonLocalData(iel, knot_pots[j], deltaV[j]));
+}
+
 /// \relates NonLocalEcpComponent
 template void NonLocalECPComponent::randomize_grid(std::vector<float>& sphere, RandomGenerator_t& myRNG);
 template void NonLocalECPComponent::randomize_grid(std::vector<double>& sphere, RandomGenerator_t& myRNG);

src/QMCHamiltonians/NonLocalECPComponent.h

@@ -83,6 +83,8 @@ private:
   std::vector<PosType> cosgrad;
   //This stores grad psi/psi - dot(u,grad psi)
   std::vector<PosType> wfngrad;
+  //This stores potential contribution per knot:
+  std::vector<RealType> knot_pots;
 
   /// scratch spaces used by evaluateValueAndDerivatives
   Matrix<ValueType> dratio;
@@ -99,17 +101,15 @@ private:
 
   ///virtual particle set: delayed initialization
   VirtualParticleSet* VP;
-  ///true, determinant localization approximation(DLA) is enabled
-  bool use_DLA;
+
+  /// build QP position deltas from the reference electron using internally stored random grid points
+  void buildQuadraturePointDeltaPositions(RealType r, const PosType& dr, std::vector<PosType>& deltaV) const;
+
+  /** finalize the calculation of $\frac{V\Psi_T}{\Psi_T}$
+   */
+  RealType calculateProjector(RealType r, const PosType& dr);
 
 public:
-#if !defined(REMOVE_TRACEMANAGER)
-  ///pointers to trace data of containing NonLocalECPotential object
-  Array<TraceReal, 1>* Ve_sample;
-  Array<TraceReal, 1>* Vi_sample;
-  bool streaming_particles;
-#endif
-
   NonLocalECPComponent();
 
   ///destructor
@@ -127,14 +127,18 @@ public:
     sgridweight_m.push_back(weight);
   }
 
-  inline void enableDLA() { use_DLA = true; }
-
   void resize_warrays(int n, int m, int l);
 
   void randomize_grid(RandomGenerator_t& myRNG);
   template<typename T>
   void randomize_grid(std::vector<T>& sphere, RandomGenerator_t& myRNG);
 
+  /** contribute local non-local move data
+   * @param iel reference electron id.
+   * @param Txy nonlocal move data.
+   */
+  void contributeTxy(int iel, std::vector<NonLocalData>& Txy) const;
+
   /** @brief Evaluate the nonlocal pp contribution via randomized quadrature grid
    * to total energy from ion "iat" and electron "iel".
    *
@@ -144,8 +148,7 @@ public:
    * @param iel index of electron
    * @param r the distance between ion iat and electron iel.
    * @param dr displacement from ion iat to electron iel.
-   * @param Tmove flag to compute tmove contributions.
-   * @param Txy nonlocal move data.
+   * @param use_DLA if ture, use determinant localization approximation (DLA).
    *
    * @return RealType Contribution to $\frac{V\Psi_T}{\Psi_T}$ from ion iat and electron iel.
    */
@@ -155,8 +158,33 @@ public:
                        int iel,
                        RealType r,
                        const PosType& dr,
-                       bool Tmove,
-                       std::vector<NonLocalData>& Txy);
+                       bool use_DLA);
+
+  /** @brief Evaluate the nonlocal pp contribution via randomized quadrature grid
+   * to total energy from ion "iat" and electron "iel" for a batch of walkers.
+   *
+   * @param ecp_component_list a list of ECP components
+   * @param p_list a list of electron particle set.
+   * @param iat_list a list of ion indices.
+   * @param psi_list a list of trial wave function object
+   * @param iel_list a list of electron indices
+   * @param r_list a list of the distances between ion iat and electron iel.
+   * @param dr_list a list of displacements from ion iat to electron iel.
+   * @param pairpots a list of contribution to $\frac{V\Psi_T}{\Psi_T}$ from ion iat and electron iel.
+   * @param use_DLA if ture, use determinant localization approximation (DLA).
+   *
+   * Note: ecp_component_list allows including different NLPP component for different walkers.
+   * electrons in iel_list must be of the same group (spin)
+   */
+  static void flex_evaluateOne(const RefVector<NonLocalECPComponent>& ecp_component_list,
+                               const RefVector<ParticleSet>& p_list,
+                               const std::vector<int>& iat_list,
+                               const RefVector<TrialWaveFunction>& psi_list,
+                               const std::vector<int>& iel_list,
+                               const std::vector<RealType>& r_list,
+                               const std::vector<PosType>& dr_list,
+                               std::vector<RealType>& pairpots,
+                               bool use_DLA);
 
   /** @brief Evaluate the nonlocal pp contribution via randomized quadrature grid
    * to total energy from ion "iat" and electron "iel".
@@ -167,9 +195,7 @@ public:
    * @param iel index of electron
    * @param r the distance between ion iat and electron iel.
    * @param dr displacement from ion iat to electron iel.
-   * @param Tmove flag to compute tmove contributions.
    * @param force_iat 3d vector for Hellman-Feynman contribution.  This gets modified.
-   * @param Txy nonlocal move data.
    *
    * @return RealType Contribution to $\frac{V\Psi_T}{\Psi_T}$ from ion iat and electron iel.
    */
@@ -179,9 +205,7 @@ public:
                                  int iel,
                                  RealType r,
                                  const PosType& dr,
-                                 PosType& force_iat,
-                                 bool Tmove,
-                                 std::vector<NonLocalData>& Txy);
+                                 PosType& force_iat);
 
   /** @brief Evaluate the nonlocal pp energy, Hellman-Feynman force, and "Pulay" force contribution
    * via randomized quadrature grid from ion "iat" and electron "iel".
@@ -195,8 +219,6 @@ public:
    * @param dr displacement from ion iat to electron iel.
    * @param force_iat 3d vector for Hellman-Feynman contribution.  This gets modified.
    * @param pulay_terms Nion x 3 object, holding a contribution for each ionic gradient from \Psi_T.
-   * @param Tmove flag to compute tmove contributions.
-   * @param Txy nonlocal move data.
    *
    * @return RealType Contribution to $\frac{V\Psi_T}{\Psi_T}$ from ion iat and electron iel.
    */
@@ -208,9 +230,7 @@ public:
                                  RealType r,
                                  const PosType& dr,
                                  PosType& force_iat,
-                                 ParticleSet::ParticlePos_t& pulay_terms,
-                                 bool Tmove,
-                                 std::vector<NonLocalData>& Txy);
+                                 ParticleSet::ParticlePos_t& pulay_terms);
 
   // This function needs to be updated to SoA. myTableIndex is introduced temporarily.
   RealType evaluateValueAndDerivatives(ParticleSet& P,

src/QMCHamiltonians/NonLocalECPotential.cpp

@@ -33,7 +33,7 @@ NonLocalECPotential::NonLocalECPotential(ParticleSet& ions,
                                          ParticleSet& els,
                                          TrialWaveFunction& psi,
                                          bool computeForces,
-                                         bool useVP)
+                                         bool enable_DLA)
     : ForceBase(ions, els),
       myRNG(nullptr),
       IonConfig(ions),
@@ -43,7 +43,8 @@ NonLocalECPotential::NonLocalECPotential(ParticleSet& ions,
       IonNeighborElecs(ions),
       UseTMove(TMOVE_OFF),
       nonLocalOps(els.getTotalNum()),
-      ComputeForces(computeForces)
+      ComputeForces(computeForces),
+      use_DLA(enable_DLA)
 {
   set_energy_domain(potential);
   two_body_quantum_domain(ions, els);
@@ -55,8 +56,12 @@ NonLocalECPotential::NonLocalECPotential(ParticleSet& ions,
   PPset.resize(IonConfig.getSpeciesSet().getTotalNum(), 0);
   PulayTerm.resize(NumIons);
   UpdateMode.set(NONLOCAL, 1);
-  Ve_samp_tmp.resize(els.getTotalNum());
-  Vi_samp_tmp.resize(ions.getTotalNum());
+  nlpp_jobs.resize(els.groups());
+  for (size_t ig = 0; ig < els.groups(); ig++)
+  {
+    // this should be enough in most calculations assuming that every electron cannot be in more than two pseudo regions.
+    nlpp_jobs[ig].reserve(2 * els.groupsize(ig));
+  }
 }
 
 ///destructor
@@ -80,15 +85,6 @@ void NonLocalECPotential::checkout_particle_quantities(TraceManager& tm)
   {
     Ve_sample = tm.checkout_real<1>(myName, Peln);
     Vi_sample = tm.checkout_real<1>(myName, IonConfig);
-    for (int iat = 0; iat < NumIons; iat++)
-    {
-      if (PP[iat])
-      {
-        PP[iat]->streaming_particles = streaming_particles;
-        //PP[iat]->Ve_sample           = Ve_sample;
-        //PP[iat]->Vi_sample           = Vi_sample;
-      }
-    }
   }
 }
 
@@ -96,15 +92,6 @@ void NonLocalECPotential::delete_particle_quantities()
 {
   if (streaming_particles)
   {
-    for (int iat = 0; iat < NumIons; iat++)
-    {
-      if (PP[iat])
-      {
-        PP[iat]->streaming_particles = false;
-        //PP[iat]->Ve_sample           = NULL;
-        //PP[iat]->Vi_sample           = NULL;
-      }
-    }
     delete Ve_sample;
     delete Vi_sample;
   }
@@ -113,10 +100,295 @@ void NonLocalECPotential::delete_particle_quantities()
 
 NonLocalECPotential::Return_t NonLocalECPotential::evaluate(ParticleSet& P)
 {
-  evaluate(P, false);
+  evaluateImpl(P, false);
   return Value;
 }
 
+void NonLocalECPotential::mw_evaluate(const RefVector<OperatorBase>& O_list, const RefVector<ParticleSet>& P_list)
+{
+  mw_evaluateImpl(O_list, P_list, false);
+}
+
+NonLocalECPotential::Return_t NonLocalECPotential::evaluateWithToperator(ParticleSet& P)
+{
+  if (UseTMove == TMOVE_V0 || UseTMove == TMOVE_V3)
+    evaluateImpl(P, true);
+  else
+    evaluateImpl(P, false);
+  return Value;
+}
+
+void NonLocalECPotential::mw_evaluateWithToperator(const RefVector<OperatorBase>& O_list,
+                                                   const RefVector<ParticleSet>& P_list)
+{
+  if (UseTMove == TMOVE_V0 || UseTMove == TMOVE_V3)
+    mw_evaluateImpl(O_list, P_list, true);
+  else
+    mw_evaluateImpl(O_list, P_list, false);
+}
+
+void NonLocalECPotential::evaluateImpl(ParticleSet& P, bool Tmove)
+{
+  if (Tmove)
+    nonLocalOps.reset();
+  std::vector<NonLocalData>& Txy(nonLocalOps.Txy);
+  Value = 0.0;
+#if !defined(REMOVE_TRACEMANAGER)
+  auto& Ve_samp = *Ve_sample;
+  auto& Vi_samp = *Vi_sample;
+  if (streaming_particles)
+  {
+    Ve_samp = 0.0;
+    Vi_samp = 0.0;
+  }
+#endif
+  for (int ipp = 0; ipp < PPset.size(); ipp++)
+    if (PPset[ipp])
+      PPset[ipp]->randomize_grid(*myRNG);
+  //loop over all the ions
+  const auto& myTable = P.getDistTable(myTableIndex);
+  // clear all the electron and ion neighbor lists
+  for (int iat = 0; iat < NumIons; iat++)
+    IonNeighborElecs.getNeighborList(iat).clear();
+  for (int jel = 0; jel < P.getTotalNum(); jel++)
+    ElecNeighborIons.getNeighborList(jel).clear();
+
+  if (ComputeForces)
+  {
+    forces = 0;
+    if (myTable.DTType == DT_SOA)
+    {
+      for (int jel = 0; jel < P.getTotalNum(); jel++)
+      {
+        const auto& dist               = myTable.getDistRow(jel);
+        const auto& displ              = myTable.getDisplRow(jel);
+        std::vector<int>& NeighborIons = ElecNeighborIons.getNeighborList(jel);
+        for (int iat = 0; iat < NumIons; iat++)
+          if (PP[iat] != nullptr && dist[iat] < PP[iat]->getRmax())
+          {
+            RealType pairpot = PP[iat]->evaluateOneWithForces(P, iat, Psi, jel, dist[iat], -displ[iat], forces[iat]);
+            if (Tmove)
+              PP[iat]->contributeTxy(jel, Txy);
+            Value += pairpot;
+            NeighborIons.push_back(iat);
+            IonNeighborElecs.getNeighborList(iat).push_back(jel);
+          }
+      }
+    }
+    else
+    {
+      APP_ABORT("NonLocalECPotential::evaluate():  Forces not imlpemented for AoS build\n");
+    }
+  }
+  else
+  {
+    if (myTable.DTType == DT_SOA)
+    {
+      for (int jel = 0; jel < P.getTotalNum(); jel++)
+      {
+        const auto& dist               = myTable.getDistRow(jel);
+        const auto& displ              = myTable.getDisplRow(jel);
+        std::vector<int>& NeighborIons = ElecNeighborIons.getNeighborList(jel);
+        for (int iat = 0; iat < NumIons; iat++)
+          if (PP[iat] != nullptr && dist[iat] < PP[iat]->getRmax())
+          {
+            RealType pairpot = PP[iat]->evaluateOne(P, iat, Psi, jel, dist[iat], -displ[iat], use_DLA);
+            if (Tmove)
+              PP[iat]->contributeTxy(jel, Txy);
+            Value += pairpot;
+            NeighborIons.push_back(iat);
+            IonNeighborElecs.getNeighborList(iat).push_back(jel);
+            if (streaming_particles)
+            {
+              Ve_samp(jel) = 0.5 * pairpot;
+              Vi_samp(iat) = 0.5 * pairpot;
+            }
+          }
+      }
+    }
+    else
+    {
+#ifndef ENABLE_SOA
+      for (int iat = 0; iat < NumIons; iat++)
+      {
+        if (PP[iat] == nullptr)
+          continue;
+        std::vector<int>& NeighborElecs = IonNeighborElecs.getNeighborList(iat);
+        for (int nn = myTable.M[iat], iel = 0; nn < myTable.M[iat + 1]; nn++, iel++)
+        {
+          const RealType r(myTable.r(nn));
+          if (r > PP[iat]->getRmax())
+            continue;
+          RealType pairpot = PP[iat]->evaluateOne(P, iat, Psi, iel, r, myTable.dr(nn), use_DLA);
+          if (Tmove)
+            PP[iat]->contributeTxy(iel, Txy);
+          Value += pairpot;
+          NeighborElecs.push_back(iel);
+          ElecNeighborIons.getNeighborList(iel).push_back(iat);
+          if (streaming_particles)
+          {
+            Ve_samp(iel) = 0.5 * pairpot;
+            Vi_samp(iat) = 0.5 * pairpot;
+          }
+        }
+      }
+#endif
+    }
+  }
+
+#if defined(TRACE_CHECK)
+  if (streaming_particles)
+  {
+    Return_t Vnow  = Value;
+    RealType Visum = Vi_sample->sum();
+    RealType Vesum = Ve_sample->sum();
+    RealType Vsum  = Vesum + Visum;
+    if (std::abs(Vsum - Vnow) > TraceManager::trace_tol)
+    {
+      app_log() << "accumtest: NonLocalECPotential::evaluate()" << std::endl;
+      app_log() << "accumtest:   tot:" << Vnow << std::endl;
+      app_log() << "accumtest:   sum:" << Vsum << std::endl;
+      APP_ABORT("Trace check failed");
+    }
+    if (std::abs(Vesum - Visum) > TraceManager::trace_tol)
+    {
+      app_log() << "sharetest: NonLocalECPotential::evaluate()" << std::endl;
+      app_log() << "sharetest:   e share:" << Vesum << std::endl;
+      app_log() << "sharetest:   i share:" << Visum << std::endl;
+      APP_ABORT("Trace check failed");
+    }
+  }
+#endif
+}
+
+void NonLocalECPotential::mw_evaluateImpl(const RefVector<OperatorBase>& O_list,
+                                          const RefVector<ParticleSet>& P_list,
+                                          bool Tmove)
+{
+  const size_t ngroups = P_list[0].get().groups();
+  const size_t nw      = O_list.size();
+  /// maximal number of jobs per spin
+  std::vector<size_t> max_num_jobs(ngroups, 0);
+
+#pragma omp parallel for
+  for (size_t iw = 0; iw < nw; iw++)
+  {
+    NonLocalECPotential& O(static_cast<NonLocalECPotential&>(O_list[iw].get()));
+    ParticleSet& P(P_list[iw]);
+
+    if (Tmove)
+      O.nonLocalOps.reset();
+
+    for (int ipp = 0; ipp < O.PPset.size(); ipp++)
+      if (O.PPset[ipp])
+        O.PPset[ipp]->randomize_grid(*O.myRNG);
+
+    //loop over all the ions
+    const auto& myTable = P.getDistTable(myTableIndex);
+    // clear all the electron and ion neighbor lists
+    for (int iat = 0; iat < NumIons; iat++)
+      O.IonNeighborElecs.getNeighborList(iat).clear();
+    for (int jel = 0; jel < P.getTotalNum(); jel++)
+      O.ElecNeighborIons.getNeighborList(jel).clear();
+
+    if (ComputeForces)
+      APP_ABORT("NonLocalECPotential::mw_evaluateImpl(): Forces not imlpemented\n");
+    if (myTable.DTType != DT_SOA)
+      APP_ABORT("NonLocalECPotential::mw_evaluateImpl(): not imlpemented for AoS builds\n");
+
+    for (int ig = 0; ig < P.groups(); ++ig) //loop over species
+    {
+      auto& joblist = O.nlpp_jobs[ig];
+      joblist.clear();
+
+      for (int jel = P.first(ig); jel < P.last(ig); ++jel)
+      {
+        const auto& dist               = myTable.getDistRow(jel);
+        const auto& displ              = myTable.getDisplRow(jel);
+        std::vector<int>& NeighborIons = O.ElecNeighborIons.getNeighborList(jel);
+        for (int iat = 0; iat < O.NumIons; iat++)
+          if (O.PP[iat] != nullptr && dist[iat] < O.PP[iat]->getRmax())
+          {
+            NeighborIons.push_back(iat);
+            O.IonNeighborElecs.getNeighborList(iat).push_back(jel);
+            joblist.push_back(NLPPJob(iat, jel, P.R[jel], dist[iat], -displ[iat]));
+          }
+      }
+      // find the max number of jobs of all the walkers
+      max_num_jobs[ig] = std::max(max_num_jobs[ig], joblist.size());
+    }
+
+    O.Value = 0.0;
+  }
+
+  RefVector<NonLocalECPotential> ecp_potential_list;
+  RefVector<NonLocalECPComponent> ecp_component_list;
+  RefVector<ParticleSet> p_list;
+  std::vector<int> iat_list;
+  RefVector<TrialWaveFunction> psi_list;
+  std::vector<int> jel_list;
+  std::vector<RealType> r_list;
+  std::vector<PosType> dr_list;
+  std::vector<RealType> pairpots(nw);
+
+  ecp_potential_list.reserve(nw);
+  ecp_component_list.reserve(nw);
+  p_list.reserve(nw);
+  iat_list.reserve(nw);
+  psi_list.reserve(nw);
+  jel_list.reserve(nw);
+  r_list.reserve(nw);
+  dr_list.reserve(nw);
+
+  for (int ig = 0; ig < ngroups; ++ig) //loop over species
+  for (size_t jobid = 0; jobid < max_num_jobs[ig]; jobid++)
+  {
+    ecp_potential_list.clear();
+    ecp_component_list.clear();
+    p_list.clear();
+    iat_list.clear();
+    psi_list.clear();
+    jel_list.clear();
+    r_list.clear();
+    dr_list.clear();
+    for (size_t iw = 0; iw < nw; iw++)
+    {
+      NonLocalECPotential& O(static_cast<NonLocalECPotential&>(O_list[iw].get()));
+      ParticleSet& P(P_list[iw]);
+      if (jobid < O.nlpp_jobs[ig].size())
+      {
+        const auto& job = O.nlpp_jobs[ig][jobid];
+        ecp_potential_list.push_back(O);
+        ecp_component_list.push_back(*O.PP[std::get<0>(job)]);
+        p_list.push_back(P);
+        iat_list.push_back(std::get<0>(job));
+        psi_list.push_back(O.Psi);
+        jel_list.push_back(std::get<1>(job));
+        r_list.push_back(std::get<3>(job));
+        dr_list.push_back(std::get<4>(job));
+      }
+    }
+
+    NonLocalECPComponent::flex_evaluateOne(ecp_component_list, p_list, iat_list, psi_list, jel_list, r_list, dr_list,
+                                           pairpots, use_DLA);
+
+    for (size_t j = 0; j < ecp_potential_list.size(); j++)
+    {
+      if (false)
+      { // code usefully for debugging
+        RealType check_value = ecp_component_list[j].get().evaluateOne(p_list[j], iat_list[j], psi_list[j], jel_list[j],
+                                                                       r_list[j], dr_list[j], use_DLA);
+        if (std::abs(check_value - pairpots[j]) > 1e-5)
+          std::cout << "check " << check_value << " wrong " << pairpots[j] << " diff "
+                    << std::abs(check_value - pairpots[j]) << std::endl;
+      }
+      ecp_potential_list[j].get().Value += pairpots[j];
+      if (Tmove)
+        ecp_component_list[j].get().contributeTxy(jel_list[j], ecp_potential_list[j].get().nonLocalOps.Txy);
+    }
+  }
+}
+
 NonLocalECPotential::Return_t NonLocalECPotential::evaluateWithIonDerivs(ParticleSet& P,
                                                                          ParticleSet& ions,
                                                                          TrialWaveFunction& psi,
@@ -155,8 +427,10 @@ NonLocalECPotential::Return_t NonLocalECPotential::evaluateWithIonDerivs(Particl
       for (int iat = 0; iat < NumIons; iat++)
         if (PP[iat] != nullptr && dist[iat] < PP[iat]->getRmax())
         {
-          Value += PP[iat]->evaluateOneWithForces(P, ions, iat, Psi, jel, dist[iat], -displ[iat],
-                                                  forces[iat], PulayTerm, Tmove, Txy);
+          Value +=
+              PP[iat]->evaluateOneWithForces(P, ions, iat, Psi, jel, dist[iat], -displ[iat], forces[iat], PulayTerm);
+          if (Tmove)
+            PP[iat]->contributeTxy(jel, Txy);
           NeighborIons.push_back(iat);
           IonNeighborElecs.getNeighborList(iat).push_back(jel);
         }
@@ -166,146 +440,11 @@ NonLocalECPotential::Return_t NonLocalECPotential::evaluateWithIonDerivs(Particl
   {
     APP_ABORT("NonLocalECPotential::evaluate():  Forces not imlpemented for AoS build\n");
   }
-  hf_terms    -= forces;
+  hf_terms -= forces;
   pulay_terms -= PulayTerm;
   return Value;
 }
 
-NonLocalECPotential::Return_t NonLocalECPotential::evaluateWithToperator(ParticleSet& P)
-{
-  if (UseTMove == TMOVE_V0 || UseTMove == TMOVE_V3)
-  {
-    nonLocalOps.reset();
-    evaluate(P, true);
-  }
-  else
-    evaluate(P, false);
-  return Value;
-}
-
-void NonLocalECPotential::evaluate(ParticleSet& P, bool Tmove)
-{
-  std::vector<NonLocalData>& Txy(nonLocalOps.Txy);
-  Value = 0.0;
-#if !defined(REMOVE_TRACEMANAGER)
-  if (streaming_particles)
-  {
-    (*Ve_sample) = 0.0;
-    (*Vi_sample) = 0.0;
-  }
-  auto& Ve_samp = Ve_samp_tmp;
-  auto& Vi_samp = Vi_samp_tmp;
-  Ve_samp = 0.0;
-  Vi_samp = 0.0;
-#endif
-  for (int ipp = 0; ipp < PPset.size(); ipp++)
-    if (PPset[ipp])
-      PPset[ipp]->randomize_grid(*myRNG);
-  //loop over all the ions
-  const auto& myTable = P.getDistTable(myTableIndex);
-  // clear all the electron and ion neighbor lists
-  for (int iat = 0; iat < NumIons; iat++)
-    IonNeighborElecs.getNeighborList(iat).clear();
-  for (int jel = 0; jel < P.getTotalNum(); jel++)
-    ElecNeighborIons.getNeighborList(jel).clear();
-
-  if (ComputeForces)
-  {
-    forces = 0;
-    if (myTable.DTType == DT_SOA)
-    {
-      for (int jel = 0; jel < P.getTotalNum(); jel++)
-      {
-        const auto& dist               = myTable.getDistRow(jel);
-        const auto& displ              = myTable.getDisplRow(jel);
-        std::vector<int>& NeighborIons = ElecNeighborIons.getNeighborList(jel);
-        for (int iat = 0; iat < NumIons; iat++)
-          if (PP[iat] != nullptr && dist[iat] < PP[iat]->getRmax())
-          {
-            RealType pairpot = PP[iat]->evaluateOneWithForces(P, iat, Psi, jel, dist[iat], -displ[iat], forces[iat],
-                                                    Tmove, Txy);
-            Value += pairpot;
-            NeighborIons.push_back(iat);
-            IonNeighborElecs.getNeighborList(iat).push_back(jel);
-          }
-      }
-    }
-    else
-    {
-      APP_ABORT("NonLocalECPotential::evaluate():  Forces not imlpemented for AoS build\n");
-    }
-  }
-  else
-  {
-    if (myTable.DTType == DT_SOA)
-    {
-      for (int jel = 0; jel < P.getTotalNum(); jel++)
-      {
-        const auto& dist               = myTable.getDistRow(jel);
-        const auto& displ              = myTable.getDisplRow(jel);
-        std::vector<int>& NeighborIons = ElecNeighborIons.getNeighborList(jel);
-        for (int iat = 0; iat < NumIons; iat++)
-          if (PP[iat] != nullptr && dist[iat] < PP[iat]->getRmax())
-          {
-            RealType pairpot = PP[iat]->evaluateOne(P, iat, Psi, jel, dist[iat], -displ[iat], Tmove, Txy);
-            Value += pairpot;
-            NeighborIons.push_back(iat);
-            IonNeighborElecs.getNeighborList(iat).push_back(jel);
-            Ve_samp(jel) = 0.5*pairpot;
-            Vi_samp(iat) = 0.5*pairpot;
-          }
-      }
-    }
-    else
-    {
-#ifndef ENABLE_SOA
-      for (int iat = 0; iat < NumIons; iat++)
-      {
-        if (PP[iat] == nullptr)
-          continue;
-        std::vector<int>& NeighborElecs = IonNeighborElecs.getNeighborList(iat);
-        for (int nn = myTable.M[iat], iel = 0; nn < myTable.M[iat + 1]; nn++, iel++)
-        {
-          const RealType r(myTable.r(nn));
-          if (r > PP[iat]->getRmax())
-            continue;
-          RealType pairpot = PP[iat]->evaluateOne(P, iat, Psi, iel, r, myTable.dr(nn), Tmove, Txy);
-          Value += pairpot;
-          NeighborElecs.push_back(iel);
-          ElecNeighborIons.getNeighborList(iel).push_back(iat);
-          Ve_samp(iel) = 0.5*pairpot;
-          Vi_samp(iat) = 0.5*pairpot;
-        }
-      }
-#endif
-    }
-  }
-
-#if defined(TRACE_CHECK)
-  if (streaming_particles)
-  {
-    Return_t Vnow  = Value;
-    RealType Visum = Vi_sample->sum();
-    RealType Vesum = Ve_sample->sum();
-    RealType Vsum  = Vesum + Visum;
-    if (std::abs(Vsum - Vnow) > TraceManager::trace_tol)
-    {
-      app_log() << "accumtest: NonLocalECPotential::evaluate()" << std::endl;
-      app_log() << "accumtest:   tot:" << Vnow << std::endl;
-      app_log() << "accumtest:   sum:" << Vsum << std::endl;
-      APP_ABORT("Trace check failed");
-    }
-    if (std::abs(Vesum - Visum) > TraceManager::trace_tol)
-    {
-      app_log() << "sharetest: NonLocalECPotential::evaluate()" << std::endl;
-      app_log() << "sharetest:   e share:" << Vesum << std::endl;
-      app_log() << "sharetest:   i share:" << Visum << std::endl;
-      APP_ABORT("Trace check failed");
-    }
-  }
-#endif
-}
-
 void NonLocalECPotential::computeOneElectronTxy(ParticleSet& P, const int ref_elec)
 {
   nonLocalOps.reset();
@@ -320,7 +459,8 @@ void NonLocalECPotential::computeOneElectronTxy(ParticleSet& P, const int ref_el
     for (int atom_index = 0; atom_index < NeighborIons.size(); atom_index++)
     {
       const int iat = NeighborIons[atom_index];
-      PP[iat]->evaluateOne(P, iat, Psi, ref_elec, dist[iat], -displ[iat], true, Txy);
+      PP[iat]->evaluateOne(P, iat, Psi, ref_elec, dist[iat], -displ[iat], use_DLA);
+      PP[iat]->contributeTxy(ref_elec, Txy);
     }
   }
   else
@@ -330,7 +470,8 @@ void NonLocalECPotential::computeOneElectronTxy(ParticleSet& P, const int ref_el
     {
       const int iat = NeighborIons[atom_index];
       int nn        = myTable.M[iat] + ref_elec;
-      PP[iat]->evaluateOne(P, iat, Psi, ref_elec, myTable.r(nn), myTable.dr(nn), true, Txy);
+      PP[iat]->evaluateOne(P, iat, Psi, ref_elec, myTable.r(nn), myTable.dr(nn), use_DLA);
+      PP[iat]->contributeTxy(ref_elec, Txy);
     }
 #endif
   }
@@ -474,17 +615,13 @@ void NonLocalECPotential::addComponent(int groupID, NonLocalECPComponent* ppot)
 {
   for (int iat = 0; iat < PP.size(); iat++)
     if (IonConfig.GroupID[iat] == groupID)
-    {
       PP[iat] = ppot;
-      ppot->Ve_sample = &Ve_samp_tmp;
-      ppot->Vi_sample = &Vi_samp_tmp;
-    }
   PPset[groupID] = ppot;
 }
 
 OperatorBase* NonLocalECPotential::makeClone(ParticleSet& qp, TrialWaveFunction& psi)
 {
-  NonLocalECPotential* myclone = new NonLocalECPotential(IonConfig, qp, psi, ComputeForces);
+  NonLocalECPotential* myclone = new NonLocalECPotential(IonConfig, qp, psi, ComputeForces, use_DLA);
   for (int ig = 0; ig < PPset.size(); ++ig)
   {
     if (PPset[ig])

src/QMCHamiltonians/NonLocalECPotential.h

@@ -22,39 +22,44 @@
 
 namespace qmcplusplus
 {
-
 class NonLocalECPComponent;
+
 /** @ingroup hamiltonian
  * \brief Evaluate the semi local potentials
  */
 class NonLocalECPotential : public OperatorBase, public ForceBase
 {
 public:
-  NonLocalECPotential(ParticleSet& ions,
-                      ParticleSet& els,
-                      TrialWaveFunction& psi,
-                      bool computeForces = false,
-                      bool useVP         = false);
+  /** tuple type for NLPP calculation of a pair of ion and electron
+   * ion id, electron id, electron position, ion electron distance, ion to electron displacement
+   */
+  using NLPPJob = std::tuple<int, int, PosType, RealType, PosType>;
+
+  NonLocalECPotential(ParticleSet& ions, ParticleSet& els, TrialWaveFunction& psi, bool computeForces, bool enable_DLA);
 
   ~NonLocalECPotential();
 
-  void resetTargetParticleSet(ParticleSet& P);
+  void resetTargetParticleSet(ParticleSet& P) override;
 
 #if !defined(REMOVE_TRACEMANAGER)
-  virtual void contribute_particle_quantities();
-  virtual void checkout_particle_quantities(TraceManager& tm);
-  virtual void delete_particle_quantities();
+  virtual void contribute_particle_quantities() override;
+  virtual void checkout_particle_quantities(TraceManager& tm) override;
+  virtual void delete_particle_quantities() override;
 #endif
 
-  Return_t evaluate(ParticleSet& P);
+  Return_t evaluate(ParticleSet& P) override;
+
+  void mw_evaluate(const RefVector<OperatorBase>& O_list, const RefVector<ParticleSet>& P_list) override;
+
+  Return_t evaluateWithToperator(ParticleSet& P) override;
+
+  void mw_evaluateWithToperator(const RefVector<OperatorBase>& O_list, const RefVector<ParticleSet>& P_list) override;
 
   Return_t evaluateWithIonDerivs(ParticleSet& P,
                                  ParticleSet& ions,
                                  TrialWaveFunction& psi,
                                  ParticleSet::ParticlePos_t& hf_terms,
-                                 ParticleSet::ParticlePos_t& pulay_terms);
-
-  Return_t evaluateWithToperator(ParticleSet& P);
+                                 ParticleSet::ParticlePos_t& pulay_terms) override;
 
   /** set non local moves options
    * @param cur the xml input
@@ -62,14 +67,11 @@ public:
   void setNonLocalMoves(xmlNodePtr cur) { UseTMove = nonLocalOps.put(cur); }
 
   void setNonLocalMoves(const std::string& non_local_move_option,
-                                        const double tau,
-                                        const double alpha,
-                                        const double gamma)
+                        const double tau,
+                        const double alpha,
+                        const double gamma)
   {
-    UseTMove = nonLocalOps.thingsThatShouldBeInMyConstructor(non_local_move_option,
-                                        tau,
-                                        alpha,
-                                        gamma);
+    UseTMove = nonLocalOps.thingsThatShouldBeInMyConstructor(non_local_move_option, tau, alpha, gamma);
   }
   /** make non local moves with particle-by-particle moves
    * @param P particle set
@@ -80,33 +82,33 @@ public:
   Return_t evaluateValueAndDerivatives(ParticleSet& P,
                                        const opt_variables_type& optvars,
                                        const std::vector<ValueType>& dlogpsi,
-                                       std::vector<ValueType>& dhpsioverpsi);
+                                       std::vector<ValueType>& dhpsioverpsi) override;
 
   /** Do nothing */
-  bool put(xmlNodePtr cur) { return true; }
+  bool put(xmlNodePtr cur) override { return true; }
 
-  bool get(std::ostream& os) const
+  bool get(std::ostream& os) const override
   {
     os << "NonLocalECPotential: " << IonConfig.getName();
     return true;
   }
 
-  OperatorBase* makeClone(ParticleSet& qp, TrialWaveFunction& psi);
+  OperatorBase* makeClone(ParticleSet& qp, TrialWaveFunction& psi) override;
 
   void addComponent(int groupID, NonLocalECPComponent* pp);
 
   /** set the internal RNG pointer as the given pointer
    * @param rng input RNG pointer
    */
-  void setRandomGenerator(RandomGenerator_t* rng) { myRNG = rng; }
+  void setRandomGenerator(RandomGenerator_t* rng) override { myRNG = rng; }
 
-  void addObservables(PropertySetType& plist, BufferType& collectables);
+  void addObservables(PropertySetType& plist, BufferType& collectables) override;
 
-  void setObservables(PropertySetType& plist);
+  void setObservables(PropertySetType& plist) override;
 
-  void setParticlePropertyList(PropertySetType& plist, int offset);
+  void setParticlePropertyList(PropertySetType& plist, int offset) override;
 
-  void registerObservables(std::vector<observable_helper*>& h5list, hid_t gid) const;
+  void registerObservables(std::vector<observable_helper*>& h5list, hid_t gid) const override;
 
 protected:
   ///random number generator
@@ -139,22 +141,30 @@ private:
   NonLocalTOperator nonLocalOps;
   ///true if we should compute forces
   bool ComputeForces;
+  ///true, determinant localization approximation(DLA) is enabled
+  bool use_DLA;
   ///Pulay force vector
   ParticleSet::ParticlePos_t PulayTerm;
 #if !defined(REMOVE_TRACEMANAGER)
   ///single particle trace samples
-  Array<TraceReal, 1> Ve_samp_tmp;
-  Array<TraceReal, 1> Vi_samp_tmp;
   Array<TraceReal, 1>* Ve_sample;
   Array<TraceReal, 1>* Vi_sample;
 #endif
+  ///NLPP job list of ion-electron pairs by spin group
+  std::vector<std::vector<NLPPJob>> nlpp_jobs;
 
   /** the actual implementation, used by evaluate and evaluateWithToperator
    * @param P particle set
    * @param Tmove whether Txy for Tmove is updated
    */
-  void evaluate(ParticleSet& P, bool Tmove);
+  void evaluateImpl(ParticleSet& P, bool Tmove);
 
+  /** the actual implementation for batched walkers, used by mw_evaluate and mw_evaluateWithToperator
+   * @param O_list the list of NonLocalECPotential in a walker batch
+   * @param P_list the list of ParticleSet in a walker batch
+   * @param Tmove whether Txy for Tmove is updated
+   */
+  void mw_evaluateImpl(const RefVector<OperatorBase>& O_list, const RefVector<ParticleSet>& P_list, bool Tmove);
 
   /** compute the T move transition probability for a given electron
    * member variable nonLocalOps.Txy is updated
@@ -163,9 +173,10 @@ private:
    */
   void computeOneElectronTxy(ParticleSet& P, const int ref_elec);
 
-  /** mark all the electrons affected by Tmoves
+  /** mark all the electrons affected by Tmoves and update ElecNeighborIons and IonNeighborElecs
    * @param myTable electron ion distance table
    * @param iel reference electron
+   * Note this funtion should be called before acceptMove for a Tmove
    */
   void markAffectedElecs(const DistanceTableData& myTable, int iel);
 };

src/QMCHamiltonians/NonLocalECPotential_CUDA.cpp

@@ -26,7 +26,7 @@ NonLocalECPotential_CUDA::NonLocalECPotential_CUDA(ParticleSet& ions,
                                                    TrialWaveFunction& psi,
                                                    bool usePBC,
                                                    bool doForces)
-    : NonLocalECPotential(ions, els, psi, doForces),
+    : NonLocalECPotential(ions, els, psi, doForces, false),
       UsePBC(usePBC),
       CurrentNumWalkers(0),
       Ions_GPU("NonLocalECPotential_CUDA::Ions_GPU"),

src/QMCHamiltonians/NonLocalTOperator.cpp

@@ -78,12 +78,12 @@ int NonLocalTOperator::thingsThatShouldBeInMyConstructor(const std::string& non_
                                                          const double alpha,
                                                          const double gamma)
 {
-  Tau = tau;
-  Alpha = alpha;
-  Gamma = gamma;
-  plusFactor   = Tau * Gamma;
-  minusFactor  = -Tau * (1.0 - Alpha * (1.0 + Gamma));
-  int v_tmove  = TMOVE_OFF;
+  Tau         = tau;
+  Alpha       = alpha;
+  Gamma       = gamma;
+  plusFactor  = Tau * Gamma;
+  minusFactor = -Tau * (1.0 - Alpha * (1.0 + Gamma));
+  int v_tmove = TMOVE_OFF;
   std::ostringstream o;
 
   if (non_local_move_option == "no")
@@ -111,7 +111,7 @@ int NonLocalTOperator::thingsThatShouldBeInMyConstructor(const std::string& non_
     APP_ABORT("NonLocalTOperator::put unknown nonlocalmove option " + non_local_move_option);
   }
   app_log() << o.str() << std::endl;
-  return v_tmove;  
+  return v_tmove;
 }
 void NonLocalTOperator::reset() { Txy.clear(); }
 

src/QMCHamiltonians/OperatorBase.h

@@ -258,10 +258,7 @@ struct OperatorBase : public QMCTraits
   /** Evaluate the contribution of this component of multiple walkers */
   virtual void mw_evaluateWithToperator(const RefVector<OperatorBase>& O_list, const RefVector<ParticleSet>& P_list)
   {
-    for(int iw = 0; iw < O_list.size(); ++iw)
-    {
-      O_list[iw].get().evaluateWithToperator(P_list[iw]);
-    }
+    mw_evaluate(O_list, P_list);
   }
 
   /** evaluate value and derivatives wrt the optimizables

src/QMCHamiltonians/QMCHamiltonian.cpp

@@ -19,7 +19,8 @@
 #include "QMCHamiltonians/QMCHamiltonian.h"
 #include "Particle/WalkerSetRef.h"
 #include "Particle/DistanceTableData.h"
-#include "QMCHamiltonians/NonLocalECPComponent.h"
+#include "QMCWaveFunctions/TrialWaveFunction.h"
+#include "QMCHamiltonians/NonLocalECPotential.h"
 #include "Utilities/NewTimer.h"
 #ifdef QMC_CUDA
 #include "Particle/MCWalkerConfiguration.h"
@@ -782,6 +783,33 @@ void QMCHamiltonian::setRandomGenerator(RandomGenerator_t* rng)
     nlpp_ptr->setRandomGenerator(rng);
 }
 
+  void QMCHamiltonian::setNonLocalMoves(xmlNodePtr cur)
+  {
+    if (nlpp_ptr != nullptr)
+      nlpp_ptr->setNonLocalMoves(cur);
+  }
+
+  void QMCHamiltonian::setNonLocalMoves(const std::string& non_local_move_option,
+                                        const double tau,
+                                        const double alpha,
+                                        const double gamma)
+  {
+    if (nlpp_ptr != nullptr)
+      nlpp_ptr->setNonLocalMoves(non_local_move_option,
+                                        tau,
+                                        alpha,
+                                        gamma);
+  }
+
+  int QMCHamiltonian::makeNonLocalMoves(ParticleSet& P)
+  {
+    if (nlpp_ptr == nullptr)
+      return 0;
+    else
+      return nlpp_ptr->makeNonLocalMovesPbyP(P);
+  }
+
+
 std::vector<int> QMCHamiltonian::flex_makeNonLocalMoves(RefVector<QMCHamiltonian>& h_list,
                                                         RefVector<ParticleSet>& p_list)
 {

src/QMCHamiltonians/QMCHamiltonian.h

@@ -21,7 +21,6 @@
 #define QMCPLUSPLUS_HAMILTONIAN_H
 #include "Configuration.h"
 #include <QMCHamiltonians/OperatorBase.h>
-#include "QMCHamiltonians/NonLocalECPotential.h"
 #if !defined(REMOVE_TRACEMANAGER)
 #include <Estimators/TraceManager.h>
 #endif
@@ -31,6 +30,7 @@ namespace qmcplusplus
 class MCWalkerConfiguration;
 class NewTimer;
 class HamiltonianFactory;
+class NonLocalECPotential;
 
 /**  Collection of Local Energy Operators
  *
@@ -270,34 +270,18 @@ public:
   /** set non local moves options
    * @param cur the xml input
    */
-  void setNonLocalMoves(xmlNodePtr cur)
-  {
-    if (nlpp_ptr != nullptr)
-      nlpp_ptr->setNonLocalMoves(cur);
-  }
+  void setNonLocalMoves(xmlNodePtr cur);
 
   void setNonLocalMoves(const std::string& non_local_move_option,
                                         const double tau,
                                         const double alpha,
-                                        const double gamma)
-  {
-    if (nlpp_ptr != nullptr)
-      nlpp_ptr->setNonLocalMoves(non_local_move_option,
-                                        tau,
-                                        alpha,
-                                        gamma);
-  }
+                                        const double gamma);
+
   /** make non local moves
    * @param P particle set
    * @return the number of accepted moves
    */
-  int makeNonLocalMoves(ParticleSet& P)
-  {
-    if (nlpp_ptr == nullptr)
-      return 0;
-    else
-      return nlpp_ptr->makeNonLocalMovesPbyP(P);
-  }
+  int makeNonLocalMoves(ParticleSet& P);
 
   static std::vector<int> flex_makeNonLocalMoves(RefVector<QMCHamiltonian>& h_list, RefVector<ParticleSet>& p_list);
   /** evaluate energy 

src/QMCHamiltonians/tests/test_ecp.cpp

@@ -291,11 +291,6 @@ TEST_CASE("Evaluate_ecp", "[hamiltonian]")
   //quadrature Lattice instead...
   copyGridUnrotatedForTest(*nlpp);
 
-
-  //Not testing nonlocal moves here, but the PP functions take this as an argument
-  NonLocalTOperator nonLocalOps(elec.getTotalNum());
-  std::vector<NonLocalData>& Txy(nonLocalOps.Txy);
-
   const int myTableIndex = elec.addTable(ions, DT_SOA_PREFERRED);
 
   const auto& myTable = elec.getDistTable(myTableIndex);
@@ -317,7 +312,7 @@ TEST_CASE("Evaluate_ecp", "[hamiltonian]")
     const auto& displ = myTable.getDisplRow(jel);
     for (int iat = 0; iat < ions.getTotalNum(); iat++)
       if (nlpp != nullptr && dist[iat] < nlpp->getRmax())
-        Value1 += nlpp->evaluateOne(elec, iat, psi, jel, dist[iat], -displ[iat], 0, Txy);
+        Value1 += nlpp->evaluateOne(elec, iat, psi, jel, dist[iat], -displ[iat], false);
   }
 #else
   for (int iat = 0; iat < ions.getTotalNum(); iat++)
@@ -329,7 +324,7 @@ TEST_CASE("Evaluate_ecp", "[hamiltonian]")
       const RealType r(myTable.r(nn));
       if (r > nlpp->getRmax())
         continue;
-      Value1 += nlpp->evaluateOne(elec, iat, psi, iel, r, myTable.dr(nn), 0, Txy);
+      Value1 += nlpp->evaluateOne(elec, iat, psi, iel, r, myTable.dr(nn), false);
     }
   }
 #endif
@@ -413,9 +408,9 @@ TEST_CASE("Evaluate_ecp", "[hamiltonian]")
     for (int iat = 0; iat < ions.getTotalNum(); iat++)
       if (nlpp != nullptr && dist[iat] < nlpp->getRmax())
       {
-        Value2 += nlpp->evaluateOneWithForces(elec, iat, psi, jel, dist[iat], -displ[iat], HFTerm[iat], 0, Txy);
+        Value2 += nlpp->evaluateOneWithForces(elec, iat, psi, jel, dist[iat], -displ[iat], HFTerm[iat]);
         Value3 += nlpp->evaluateOneWithForces(elec, ions, iat, psi, jel, dist[iat], -displ[iat], HFTerm2[iat],
-                                              PulayTerm, 0, Txy);
+                                              PulayTerm);
       }
   }
   //These values are validated against print statements.

src/QMCWaveFunctions/Fermion/DiracDeterminant.cpp

@@ -375,7 +375,7 @@ typename DiracDeterminant<DU_TYPE>::PsiValueType DiracDeterminant<DU_TYPE>::rati
 }
 
 template<typename DU_TYPE>
-void DiracDeterminant<DU_TYPE>::evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios)
+void DiracDeterminant<DU_TYPE>::evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios)
 {
   SPOVTimer.start();
   const int WorkingIndex = VP.refPtcl - FirstIndex;

src/QMCWaveFunctions/Fermion/DiracDeterminant.h

@@ -96,7 +96,7 @@ public:
 
   /** compute multiple ratios for a particle move
    */
-  void evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios) override;
+  void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios) override;
 
   PsiValueType ratioGrad(ParticleSet& P, int iat, GradType& grad_iat) override;
 

src/QMCWaveFunctions/Fermion/SlaterDet.h

@@ -90,7 +90,7 @@ public:
 
   virtual void copyFromBuffer(ParticleSet& P, WFBufferType& buf) override;
 
-  virtual inline void evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios) override
+  virtual inline void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios) override
   {
     return Dets[getDetID(VP.refPtcl)]->evaluateRatios(VP, ratios);
   }

src/QMCWaveFunctions/Jastrow/J1OrbitalSoA.h

@@ -160,7 +160,7 @@ struct J1OrbitalSoA : public WaveFunctionComponent
     return std::exp(static_cast<PsiValueType>(Vat[iat] - curAt));
   }
 
-  inline void evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios)
+  inline void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios)
   {
     for (int k = 0; k < ratios.size(); ++k)
       ratios[k] = std::exp(Vat[VP.refPtcl] - computeU(VP.getDistTable(myTableID).getDistRow(k)));

src/QMCWaveFunctions/Jastrow/J2OrbitalSoA.h

@@ -250,7 +250,7 @@ public:
   void recompute(ParticleSet& P);
 
   PsiValueType ratio(ParticleSet& P, int iat);
-  void evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios)
+  void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios)
   {
     for (int k = 0; k < ratios.size(); ++k)
       ratios[k] =

src/QMCWaveFunctions/Jastrow/JeeIOrbitalSoA.h

@@ -425,7 +425,7 @@ public:
     return std::exp(static_cast<PsiValueType>(DiffVal));
   }
 
-  void evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios)
+  void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios)
   {
     for (int k = 0; k < ratios.size(); ++k)
       ratios[k] = std::exp(Uat[VP.refPtcl] -

src/QMCWaveFunctions/Jastrow/OneBodyJastrowOrbital.h

@@ -262,7 +262,7 @@ public:
     return std::exp(static_cast<PsiValueType>(U[iat] - curVal));
   }
 
-  inline void evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios)
+  inline void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios)
   {
     //vector<RealType> myr(ratios.size(),0.0);
     //const DistanceTableData* d_table=VP.getDistTable(myTableIndex);

src/QMCWaveFunctions/Jastrow/OneBodySpinJastrowOrbital.h

@@ -259,7 +259,7 @@ public:
     return std::exp(static_cast<PsiValueType>(U[iat] - curVal));
   }
 
-  inline void evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios)
+  inline void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios)
   {
     std::vector<RealType> myr(ratios.size(), U[VP.refPtcl]);
     const auto& d_table = VP.getDistTable(myTableIndex);

src/QMCWaveFunctions/Jastrow/TwoBodyJastrowOrbital.h

@@ -310,7 +310,7 @@ public:
     return std::exp(static_cast<PsiValueType>(DiffVal));
   }
 
-  inline void evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios)
+  inline void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios)
   {
     const int iat = VP.refPtcl;
 

src/QMCWaveFunctions/Jastrow/eeI_JastrowOrbital.h

@@ -705,7 +705,7 @@ public:
     return G;
   }
 
-  inline void evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios)
+  inline void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios)
   {
     const int iat = VP.refPtcl;
     const int nk  = ratios.size();

src/QMCWaveFunctions/TrialWaveFunction.cpp

@@ -332,13 +332,12 @@ TrialWaveFunction::ValueType TrialWaveFunction::calcRatio(ParticleSet& P, int ia
 {
   PsiValueType r(1.0);
   for (int i = 0, ii = V_TIMER; i < Z.size(); i++, ii += TIMER_SKIP)
-  {
-    myTimers[ii]->start();
     if (ct == ComputeType::ALL || (Z[i]->is_fermionic && ct == ComputeType::FERMIONIC) ||
         (!Z[i]->is_fermionic && ct == ComputeType::NONFERMIONIC))
+    {
+      ScopedTimer local_timer(myTimers[ii]);
       r *= Z[i]->ratio(P, iat);
-    myTimers[ii]->stop();
-  }
+    }
   return static_cast<ValueType>(r);
 }
 
@@ -640,19 +639,19 @@ void TrialWaveFunction::completeUpdates()
   }
 }
 
-void TrialWaveFunction::flex_completeUpdates(const std::vector<TrialWaveFunction*>& WF_list) const
+void TrialWaveFunction::flex_completeUpdates(const std::vector<TrialWaveFunction*>& wf_list) const
 {
-  if (WF_list.size() > 1)
+  if (wf_list.size() > 1)
   {
     for (int i = 0, ii = ACCEPT_TIMER; i < Z.size(); i++, ii += TIMER_SKIP)
     {
       ScopedTimer local_timer(myTimers[ii]);
-      std::vector<WaveFunctionComponent*> WFC_list(extractWFCPtrList(WF_list, i));
-      Z[i]->mw_completeUpdates(WFC_list);
+      std::vector<WaveFunctionComponent*> wfc_list(extractWFCPtrList(wf_list, i));
+      Z[i]->mw_completeUpdates(wfc_list);
     }
   }
-  else if (WF_list.size() == 1)
-    WF_list[0]->completeUpdates();
+  else if (wf_list.size() == 1)
+    wf_list[0]->completeUpdates();
 }
 
 void TrialWaveFunction::checkInVariables(opt_variables_type& active)
@@ -867,19 +866,54 @@ void TrialWaveFunction::flex_copyFromBuffer(const RefVector<TrialWaveFunction>&
     bufGetTwfLog(buf_list[iw], wf_list[iw]);
 }
 
-void TrialWaveFunction::evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios)
+void TrialWaveFunction::evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios, ComputeType ct)
 {
   assert(VP.getTotalNum() == ratios.size());
   std::vector<ValueType> t(ratios.size());
   std::fill(ratios.begin(), ratios.end(), 1.0);
   for (int i = 0, ii = NL_TIMER; i < Z.size(); ++i, ii += TIMER_SKIP)
+    if (ct == ComputeType::ALL || (Z[i]->is_fermionic && ct == ComputeType::FERMIONIC) ||
+        (!Z[i]->is_fermionic && ct == ComputeType::NONFERMIONIC))
+    {
+      ScopedTimer local_timer(myTimers[ii]);
+      Z[i]->evaluateRatios(VP, t);
+      for (int j = 0; j < ratios.size(); ++j)
+        ratios[j] *= t[j];
+    }
+}
+
+void TrialWaveFunction::flex_evaluateRatios(const RefVector<TrialWaveFunction>& wf_list, const RefVector<const VirtualParticleSet>& vp_list, const RefVector<std::vector<ValueType>>& ratios_list, ComputeType ct)
+{
+  if (wf_list.size() > 1)
   {
-    myTimers[ii]->start();
-    Z[i]->evaluateRatios(VP, t);
-    for (int j = 0; j < ratios.size(); ++j)
-      ratios[j] *= t[j];
-    myTimers[ii]->stop();
+    auto& wavefunction_components = wf_list[0].get().Z;
+    std::vector<std::vector<ValueType>> t(ratios_list.size());
+    for (int iw = 0; iw < wf_list.size(); iw++)
+    {
+      std::vector<ValueType>& ratios = ratios_list[iw];
+      assert(vp_list[iw].get().getTotalNum() == ratios.size());
+      std::fill(ratios.begin(), ratios.end(), 1.0);
+      t[iw].resize(ratios.size());
+    }
+
+    for (int i = 0, ii = NL_TIMER; i < wavefunction_components.size(); i++, ii += TIMER_SKIP)
+      if (ct == ComputeType::ALL || (wavefunction_components[i]->is_fermionic && ct == ComputeType::FERMIONIC) ||
+        (!wavefunction_components[i]->is_fermionic && ct == ComputeType::NONFERMIONIC))
+      {
+        ScopedTimer local_timer(wf_list[0].get().get_timers()[ii]);
+        const auto wfc_list(extractWFCRefList(wf_list, i));
+        wavefunction_components[i]->mw_evaluateRatios(wfc_list, vp_list, t);
+        for (int iw = 0; iw < wf_list.size(); iw++)
+        {
+          std::vector<ValueType>& ratios = ratios_list[iw];
+          for (int j = 0; j < ratios.size(); ++j)
+            ratios[j] *= t[iw][j];
+        }
+      }
   }
+  else if (wf_list.size() == 1)
+    wf_list[0].get().evaluateRatios(vp_list[0], ratios_list[0], ct);
+
 }
 
 void TrialWaveFunction::evaluateDerivRatios(VirtualParticleSet& VP,

src/QMCWaveFunctions/TrialWaveFunction.h

@@ -183,7 +183,9 @@ public:
 
   /** compulte multiple ratios to handle non-local moves and other virtual moves
    */
-  void evaluateRatios(VirtualParticleSet& P, std::vector<ValueType>& ratios);
+  void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios, ComputeType ct = ComputeType::ALL);
+  static void flex_evaluateRatios(const RefVector<TrialWaveFunction>& WF_list, const RefVector<const VirtualParticleSet>& VP_list, const RefVector<std::vector<ValueType>>& ratios_list, ComputeType ct = ComputeType::ALL);
+
   /** compute both ratios and deriatives of ratio with respect to the optimizables*/
   void evaluateDerivRatios(VirtualParticleSet& P,
                            const opt_variables_type& optvars,

src/QMCWaveFunctions/WaveFunctionComponent.cpp

@@ -83,7 +83,7 @@ void WaveFunctionComponent::evaluateRatiosAlltoOne(ParticleSet& P, std::vector<V
     ratios[i] = ratio(P, i);
 }
 
-void WaveFunctionComponent::evaluateRatios(VirtualParticleSet& P, std::vector<ValueType>& ratios)
+void WaveFunctionComponent::evaluateRatios(const VirtualParticleSet& P, std::vector<ValueType>& ratios)
 {
   std::ostringstream o;
   o << "WaveFunctionComponent::evaluateRatios is not implemented by " << ClassName;

src/QMCWaveFunctions/WaveFunctionComponent.h

@@ -615,7 +615,21 @@ struct WaveFunctionComponent : public QMCTraits
    * @param VP VirtualParticleSet
    * @param ratios ratios with new positions VP.R[k] the VP.refPtcl
    */
-  virtual void evaluateRatios(VirtualParticleSet& VP, std::vector<ValueType>& ratios);
+  virtual void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios);
+
+  /** evaluate ratios to evaluate the non-local PP multiple walkers
+   * @param wfc_list the list of WaveFunctionComponent references of the same component in a walker batch
+   * @param vp_list the list of VirtualParticleSet references in a walker batch
+   * @param ratios of all the virtual moves of all the walkers
+   */
+  virtual void mw_evaluateRatios(const RefVector<WaveFunctionComponent>& wfc_list,
+                                 const RefVector<const VirtualParticleSet>& vp_list,
+                                 std::vector<std::vector<ValueType>>& ratios)
+  {
+#pragma omp parallel for
+    for (int iw = 0; iw < wfc_list.size(); iw++)
+      wfc_list[iw].get().evaluateRatios(vp_list[iw], ratios[iw]);
+  }
 
   /** evaluate ratios to evaluate the non-local PP
    * @param VP VirtualParticleSet

src/QMCWaveFunctions/tests/Carbon1x1x1-tw1_gen_mos.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 import numpy
 import h5py
 from pyscf.pbc import gto, scf, dft, df

src/QMCWaveFunctions/tests/test_bspline_jastrow.cpp

@@ -244,10 +244,11 @@ TEST_CASE("BSpline builder Jastrow J2", "[wavefunction]")
   REQUIRE(std::real(ratio_0) == Approx(0.9522052017));
 
   VirtualParticleSet VP(elec_, 2);
-  ParticleSet::ParticlePos_t newpos2(2);
-  newpos2[0] = newpos;
-  newpos2[1] = PosType(0.2, 0.5, 0.3);
-  VP.makeMoves(1, newpos2);
+  std::vector<PosType> newpos2(2);
+  std::vector<ValueType> ratios2(2);
+  newpos2[0] = newpos - elec_.R[1];
+  newpos2[1] = PosType(0.2, 0.5, 0.3) - elec_.R[1];
+  VP.makeMoves(1, elec_.R[1], newpos2);
   j2->evaluateRatios(VP, ratios);
 
   REQUIRE(std::real(ratios[0]) == Approx(0.9871985577));

src/QMCWaveFunctions/tests/test_polynomial_eeI_jastrow.cpp

@@ -190,11 +190,11 @@ TEST_CASE("PolynomialFunctor3D Jastrow", "[wavefunction]")
   REQUIRE(std::real(dhpsioverpsi[43]) == Approx(-2.3246270644e+05));
 
   VirtualParticleSet VP(elec_, 2);
-  ParticleSet::ParticlePos_t newpos2(2);
+  std::vector<PosType> newpos2(2);
   std::vector<ValueType> ratios2(2);
-  newpos2[0] = newpos;
-  newpos2[1] = PosType(0.2, 0.5, 0.3);
-  VP.makeMoves(1, newpos2);
+  newpos2[0] = newpos - elec_.R[1];
+  newpos2[1] = PosType(0.2, 0.5, 0.3) - elec_.R[1];
+  VP.makeMoves(1, elec_.R[1], newpos2);
   j3->evaluateRatios(VP, ratios2);
 
   REQUIRE(std::real(ratios2[0]) == Approx(1.0357541137));

tests/afqmc/C_1x1x1_dzvp/pyscf/scf.rhf.out

@@ -1,5 +1,5 @@
 #INFO: **** input file is /p/lscratchh/mmorale/PYSCF/debug/C/afqmc_rhf/pyscf/scf.py ****
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 import numpy
 from functools import reduce

tests/afqmc/C_1x1x1_dzvp/pyscf/scf.rhf.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 import numpy
 from functools import reduce

tests/afqmc/C_1x1x1_dzvp/pyscf/scf.uhf.out

@@ -1,5 +1,5 @@
 #INFO: **** input file is /g/g90/malone14/projects/qmcpack/tests/afqmc/C_uhf/pyscf/scf.py ****
-#!/usr/bin/env python
+#! /usr/bin/env python3
 
 import numpy
 from functools import reduce