Repository: WagnerGroup/pyqmc (language: python)

Issue:
Slater determinant refactor
Will and I have discussed about how to clean up the current situation. A potential plan: 

1) Start with mulislater.py
2) factor out the orbital construction evaluation from slater.py. This is the only part that cares about PBCs.
3) Insert the orbital stuff into multislater.py.
4) Fix all the tests to make sure they now use multislater.py
5) Change default_slater so that it uses multislater.

This should likely be done before #239 is fixed.

fix twist bugs in orbital evaluator



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index 88f07c3b..956719d3 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -36,4 +36,4 @@ jobs:
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest
       run: |
-        PYTHONPATH="." pytest --durations=100 --reruns 3
+        PYTHONPATH="." pytest --durations=100 --reruns 3 -k "not slow"
diff --git a/benchmarks/c_solid_benchmark.py b/benchmarks/c_solid_benchmark.py
new file mode 100644
index 00000000..936765a2
--- /dev/null
+++ b/benchmarks/c_solid_benchmark.py
@@ -0,0 +1,53 @@
+# Write the benchmarking functions here.
+# See "Writing benchmarks" in the asv docs for more information.
+import pyscf
+import pyqmc
+import pyscf.pbc
+import numpy as np
+import pyqmc.eval_ecp
+np.random.seed(1234)
+
+class SiSuite:
+    """
+    
+    """
+    def setup(self):
+
+        self.cell = pyscf.pbc.gto.Cell()
+        self.cell.atom = '''C     0.      0.      0.    
+              C     0.8917  0.8917  0.8917
+              C     1.7834  1.7834  0.    
+              C     2.6751  2.6751  0.8917
+              C     1.7834  0.      1.7834
+              C     2.6751  0.8917  2.6751
+              C     0.      1.7834  1.7834
+              C     0.8917  2.6751  2.6751'''
+        self.cell.basis = f'ccecpccpvdz'
+        self.cell.pseudo = 'ccecp'
+        self.cell.a = np.eye(3)*3.5668
+        self.cell.exp_to_discard=0.2
+        self.cell.build()
+
+        self.configs = pyqmc.initial_guess(self.cell, 500)
+        self.jastrow, self.jastrow_to_opt =  pyqmc.default_jastrow(self.cell)
+        self.jastrow.recompute(self.configs)
+
+    def time_ecp_jastrow(self):
+        pyqmc.eval_ecp.ecp(self.cell,self.configs,self.jastrow, 10)
+        
+    def time_jastrow_recompute(self):
+        self.jastrow.recompute(self.configs)
+
+    def time_jastrow_laplacian(self):
+        self.jastrow.gradient_laplacian(0, self.configs.electron(0))
+
+    def time_jastrow_derivatives(self):
+        self.jastrow.pgradient()
+
+
+
+
+if __name__=="__main__":
+    h2o = SiSuite()
+    h2o.setup()
+    h2o.time_jastrow_recompute()
\ No newline at end of file
diff --git a/pyqmc/__init__.py b/pyqmc/__init__.py
index a43a196a..100f303a 100644
--- a/pyqmc/__init__.py
+++ b/pyqmc/__init__.py
@@ -1,9 +1,7 @@
 name = "pyqmc"
 from pyqmc.recipes import OPTIMIZE, VMC, DMC
 from pyqmc.mc import vmc, initial_guess
-from pyqmc.slater import PySCFSlater
-from pyqmc.multislater import MultiSlater
-
+from pyqmc.slater import Slater
 from pyqmc.multiplywf import MultiplyWF
 from pyqmc.jastrowspin import JastrowSpin
 from pyqmc.manybody_jastrow import J3
@@ -49,7 +47,7 @@ def default_slater(
     Returns:
       slater, to_opt
     """
-    wf = PySCFSlater(mol, mf, twist=twist)
+    wf = Slater(mol, mf, twist=twist)
     to_opt = {}
     if optimize_orbitals:
         for k in ["mo_coeff_alpha", "mo_coeff_beta"]:
@@ -65,7 +63,7 @@ def default_multislater(
 ):
     import numpy as np
 
-    wf = MultiSlater(mol, mf, mc, tol)
+    wf = Slater(mol, mf, mc, tol)
     to_opt = ["det_coeff"]
     to_opt = {"det_coeff": np.ones(wf.parameters["det_coeff"].shape).astype(bool)}
     to_opt["det_coeff"][0] = False  # Determinant coefficient pivot
diff --git a/pyqmc/coord.py b/pyqmc/coord.py
index 96940cdb..ff5469ba 100644
--- a/pyqmc/coord.py
+++ b/pyqmc/coord.py
@@ -4,6 +4,12 @@
 import copy
 
 
+class OpenElectron:
+    def __init__(self, epos, dist):
+        self.configs=epos
+        self.dist=dist
+
+
 class OpenConfigs:
     def __init__(self, configs):
         self.configs = configs
@@ -15,14 +21,14 @@ def electron(self, e):
     def mask(self, mask):
         return OpenConfigs(self.configs[mask])
 
-    def make_irreducible(self, e, vec):
+    def make_irreducible(self, e, vec, mask=True):
         """ 
           Input: 
             e: unused electron index
             vec: a (nconfig, 3) vector 
           Output: OpenConfigs object with just one electron
         """
-        return OpenConfigs(vec)
+        return OpenElectron(vec, self.dist)
 
     def move(self, e, new, accept):
         """
@@ -92,6 +98,22 @@ def load_hdf(self, hdf):
         self.configs = np.array(hdf["configs"])
 
 
+class PeriodicElectron:
+    """
+    Represents the coordinates of a test electron position, for many walkers and 
+    potentially several different points.
+
+    configs is a 2D or 3D vector with elements [config, point, dimension]
+    wrap is same shape as configs
+    lvec and dist will most likely be references to the parent object
+    """
+    def __init__(self, epos, lattice_vectors, dist, wrap=None):
+        self.configs=epos
+        self.lvec=lattice_vectors
+        self.wrap = wrap if wrap is not None else np.zeros_like(epos)
+        self.dist=dist
+    
+
 class PeriodicConfigs:
     def __init__(self, configs, lattice_vectors, wrap=None):
         configs, wrap_ = enforce_pbc(lattice_vectors, configs)
@@ -108,16 +130,21 @@ def electron(self, e):
     def mask(self, mask):
         return PeriodicConfigs(self.configs[mask], self.lvecs, wrap=self.wrap[mask])
 
-    def make_irreducible(self, e, vec):
+    def make_irreducible(self, e, vec, mask=None):
         """ 
-         Input: a (nconfig, 3) vector 
-         Output: a tuple with the wrapped vector and the number of wraps
-        """
-        epos, wrap = enforce_pbc(self.lvecs, vec)
-        currentwrap = self.wrap if len(self.wrap.shape) == 2 else self.wrap[:, e]
-        if len(vec.shape) == 3:
-            currentwrap = currentwrap[:, np.newaxis]
-        return PeriodicConfigs(epos, self.lvecs, wrap=wrap + currentwrap)
+         Input: a (nconfig, 3) vector or a (nconfig, N, 3) vector
+         Output: A Periodic Electron
+        """
+        if mask is None:
+            mask = np.ones(vec.shape[0:-1],dtype=bool)
+        epos_, wrap_ = enforce_pbc(self.lvecs, vec[mask])
+        epos = vec.copy()
+        epos[mask] = epos_
+        wrap = self.wrap[:,e,:].copy()
+        if len(vec.shape)==3:
+            wrap=np.repeat(self.wrap[:,e][:,np.newaxis],vec.shape[1], axis=1)
+        wrap[mask] += wrap_
+        return PeriodicElectron(epos, self.lvecs, wrap=wrap, dist=self.dist)
 
     def move(self, e, new, accept):
         """
diff --git a/pyqmc/cvmc.py b/pyqmc/cvmc.py
deleted file mode 100644
index a03b9ccf..00000000
--- a/pyqmc/cvmc.py
+++ /dev/null
@@ -1,372 +0,0 @@
-import pyqmc
-import numpy as np
-import h5py
-
-
-class DescriptorFromOBDM:
-    """
-    The reason that this has to be an object here is that parsl doesn't support 
-    functions and objects that are defined in __main__.
-    """
-
-    def __init__(self, mapping, norm=1.0):
-        """mapping should be a dictionary such that each descriptor has 
-        nret lists of weights and indices to add together
-        For example, 
-        {'t': [ 
-                 [ (1.0, (0,1)), 
-                   (1.0, (1,0)) 
-                   ], 
-                 [ (1.0, (0,1)), 
-                    (1.0,(1,0)) ] 
-               ] 
-        }
-        """
-        self.norm = norm
-        self.mapping = mapping
-        pass
-
-    def __call__(self, rets):
-        """
-        Return a dictionary of descriptors
-        """
-        avgvals = {}
-        for k, mapping in self.mapping.items():
-            n = rets[0]["value"].shape[0]
-            totsum = np.zeros(n)
-            for ret, ellist in zip(rets, mapping):
-                for w, ind in ellist:
-                    totsum += self.norm * w * ret["value"][:, ind[0], ind[1]]
-            avgvals[k] = totsum
-
-        return avgvals
-
-
-class DescriptorFromTBDM:
-    """
-    The reason that this has to be an object here is that parsl doesn't support 
-    functions and objects that are defined in __main__. 
-    """
-
-    def __init__(self, mapping, norm=1.0):
-        """mapping should be a dictionary such that each descriptor has 
-        nret lists of weights and indices to add together
-        For example, 
-        {'U': [
-                [], #up - up
-                [(0.5, (0))], #up - dn
-                [(0.5, (0))], #dn - up
-                [], #up - up
-        ]
-        """
-        self.norm = norm
-        self.mapping = mapping
-        pass
-
-    def __call__(self, rets):
-        """
-        Return a dictionary of descriptors
-        """
-        avgvals = {}
-        for k, mapping in self.mapping.items():
-            n = rets[0]["value"].shape[0]
-            totsum = np.zeros(n)
-            for ret, ellist in zip(rets, mapping):
-                for w, ind in ellist:
-                    totsum += self.norm * w * ret["value"][:, ind]
-            avgvals[k] = totsum
-
-        return avgvals
-
-
-class PGradDescriptor:
-    """   """
-
-    def __init__(self, enacc, transform, dm_evaluators, descriptors, nodal_cutoff=1e-3):
-        """ 
-        dm_evaluators: dictionary of density matrix accumulator objects   
-        descriptors : dictionary of function-like objects that translate dm_evaluators[key] to 
-                      a dictionary of descriptors
-        """
-        self.enacc = enacc
-        self.transform = transform
-        self.nodal_cutoff = nodal_cutoff
-        self.dm_evaluators = dm_evaluators
-        self.descriptors = descriptors
-
-    def _node_regr(self, configs, wf):
-        """ 
-        Return true if a given configuration is within nodal_cutoff 
-        of the node 
-        Also return the regularization polynomial if true, 
-        f = a * r ** 2 + b * r ** 4 + c * r ** 3
-        """
-        ne = configs.configs.shape[1]
-        d2 = 0.0
-        for e in range(ne):
-            d2 += np.sum(wf.gradient(e, configs.electron(e)) ** 2, axis=0)
-        r = 1.0 / d2
-        mask = r < self.nodal_cutoff ** 2
-
-        c = 7.0 / (self.nodal_cutoff ** 6)
-        b = -15.0 / (self.nodal_cutoff ** 4)
-        a = 9.0 / (self.nodal_cutoff ** 2)
-
-        f = a * r + b * r ** 2 + c * r ** 3
-        f[np.logical_not(mask)] = 1.0
-
-        return mask, f
-
-    def __call__(self, configs, wf):
-        pgrad = wf.pgradient()
-        d = self.enacc(configs, wf)
-        energy = d["total"]
-        dp = self.transform.serialize_gradients(pgrad)
-
-        node_cut, f = self._node_regr(configs, wf)
-        dp_regularized = dp * f[:, np.newaxis]
-
-        d["dpH"] = np.einsum("i,ij->ij", energy, dp_regularized)
-        d["dppsi"] = dp_regularized
-        d["dpidpj"] = np.einsum("ij,ik->ijk", dp, dp_regularized)
-        raise NotImplementedError("define __call__ for PGradOBDMTransform")
-        return d
-
-    def avg(self, configs, wf):
-        nconf = configs.configs.shape[0]
-        pgrad = wf.pgradient()
-        den = self.enacc(configs, wf)
-        energy = den["total"]
-        dp = self.transform.serialize_gradients(pgrad)
-
-        descript = {}
-        for dm_type, evaluators in self.dm_evaluators.items():
-            dms = [evaluate(configs, wf) for evaluate in evaluators]
-            descript.update(self.descriptors[dm_type](dms))
-
-        node_cut, f = self._node_regr(configs, wf)
-        dp_regularized = dp * f[:, np.newaxis]
-
-        d = {}
-        for k, it in den.items():
-            d[k] = np.mean(it, axis=0)
-        d["dpH"] = np.einsum("i,ij->j", energy, dp_regularized) / nconf
-        d["dppsi"] = np.mean(dp_regularized, axis=0)
-        d["dpidpj"] = np.einsum("ij,ik->jk", dp, dp_regularized) / nconf
-
-        for di, desc in descript.items():
-            d["dp" + di] = np.einsum("ij,i->j", dp_regularized, desc) / nconf
-            d["avg" + di] = np.sum(desc) / nconf
-        d["nodal_cutoff"] = np.sum(node_cut)
-
-        return d
-
-
-def cvmc_optimize(
-    wf,
-    configs,
-    acc,
-    objective,
-    forcing,
-    iters=10,
-    tstep=0.5,
-    npts=10,
-    datafile=None,
-    vmc=None,
-    vmcoptions=None,
-    lm=None,
-    lmoptions=None,
-    update=pyqmc.linemin.sr_update,
-    update_kws=None,
-    hdf_file=None,
-):
-    """
-    Args:
-
-       wf : a wave function object
-
-       configs : starting configurations for VMC. Ideally equilibrated with wf.
-
-       acc : A PGradDescriptor object which generates descriptors
-
-       objective : A dictionary which has one value for every descriptor returned by acc
-
-       forcing : A dictionary which has one value for every descriptor returned by acc
-    """
-    import pandas as pd
-
-    if vmc is None:
-        vmc = pyqmc.vmc
-    if vmcoptions is None:
-        vmcoptions = {}
-    if lm is None:
-        lm = lm_cvmc
-    if lmoptions is None:
-        lmoptions = {}
-    if update_kws is None:
-        update_kws = {}
-
-    # Restart
-    if hdf_file is not None:
-        with h5py.File(hdf_file, "a") as hdf:
-            if "wf" in hdf.keys():
-                grp = hdf["wf"]
-                for k in grp.keys():
-                    wf.parameters[k] = np.array(grp[k])
-
-    # Attributes for cvmc solve
-    attr = dict(iters=iters, npts=npts, tstep=tstep)
-    for k, it in objective.items():
-        attr["objective_" + k] = it
-    for k, it in forcing.items():
-        attr["forcing_" + k] = it
-
-    def get_obj_deriv(x):
-        nonlocal configs
-        for k, p in acc.transform.deserialize(x).items():
-            wf.parameters[k] = p
-        df, configs = vmc(wf, configs, accumulators={"grad": acc}, **vmcoptions)
-
-        df = pd.DataFrame(df)
-        dpavg = np.mean(df["graddppsi"], axis=0)
-        Sij = np.mean(df["graddpidpj"], axis=0) - np.einsum("i,j->ij", dpavg, dpavg)
-
-        havg = np.mean(df["gradtotal"], axis=0)
-        dpH = np.mean(df["graddpH"], axis=0)
-        dEdp = dpH - dpavg * havg
-
-        qavg = {}
-        qdp = {}
-        distfromobj = 0.0
-        objderiv = dEdp
-        objfunc = havg
-
-        for k, force in forcing.items():
-            qavg[k] = np.mean(df["gradavg" + k])
-            qdp[k] = np.mean(df["graddp" + k]) - dpavg * qavg[k]
-            distobj = qavg[k] - objective[k]
-            objderiv += 2 * force * distobj * qdp[k]
-            distfromobj += distobj
-            objfunc += force * distobj ** 2
-
-        dret = {
-            "objderiv": objderiv,
-            "energy": havg,
-            "objfunc": objfunc,
-            "dist": distfromobj,
-            "dEdp": dEdp,
-        }
-        for k, avg in qavg.items():
-            dret["avg" + k] = avg
-        for k, avg in qdp.items():
-            dret["dp" + k] = avg
-        return dret, Sij
-
-    x0 = acc.transform.serialize_parameters(wf.parameters)
-
-    df = []
-    for it in range(iters):
-        grad, Sij = get_obj_deriv(x0)
-        grad["iteration"] = it
-        grad["parameters"] = x0.copy()
-        for k, force in forcing.items():
-            print(k, grad["avg" + k], grad["dp" + k], flush=True)
-        xfit = []
-        yfit = []
-
-        taus = np.linspace(0, tstep, npts + 1)
-        taus[0] = -tstep / (npts - 1)
-        params = [x0 + update(grad["objderiv"], Sij, tau, **update_kws) for tau in taus]
-        stepsdata = lm(wf, configs, params, acc, **lmoptions)
-
-        for data, p, tau in zip(stepsdata, params, taus):
-            en = np.mean(data["total"] * data["weight"]) / np.mean(data["weight"])
-
-            qavg = {}
-            distfromobj = 0.0
-            objfunc = en
-            for k, force in forcing.items():
-                qavg[k] = np.mean(data[k] * data["weight"]) / np.mean(data["weight"])
-                distobj = qavg[k] - objective[k]
-                distfromobj += distobj
-                objfunc += force * distobj ** 2
-
-            xfit.append(tau)
-            yfit.append(objfunc)
-
-        est_min = pyqmc.linemin.stable_fit(xfit, yfit)
-        x0 = x0 + update(grad["objderiv"], Sij, est_min, **update_kws)
-
-        grad["yfit"] = yfit
-        grad["taus"] = xfit
-        pyqmc.linemin.opt_hdf(
-            hdf_file, grad, attr, configs, acc.transform.deserialize(x0)
-        )
-
-        df.append(grad)
-        if datafile is not None:
-            pd.DataFrame(df).to_json(datafile)
-
-    for k, p in acc.transform.deserialize(x0).items():
-        wf.parameters[k] = p
-
-    return wf, df
-
-
-def lm_cvmc(wf, configs, params, acc):
-    """ 
-    Evaluates accumulator on the same set of configs for correlated sampling of different wave function parameters
-
-    Args:
-        wf: wave function object
-        configs: (nconf, nelec, 3) array
-        params: (nsteps, nparams) array 
-            list of arrays of parameters (serialized) at each step
-        acc: PGradDescriptor 
-
-    Returns:
-        data: list of dicts, one dict for each sample
-            each dict contains arrays returned from PGradDescriptor, weighted by psi**2/psi0**2
-    """
-
-    import copy
-    import numpy as np
-
-    data = []
-    psi0 = wf.recompute(configs)[1]  # recompute gives logdet
-
-    # Aggregate evaluator configurations
-    extra_configs = []
-    auxassignments = []
-    for dm_type, evaluators in acc.dm_evaluators.items():
-        for evaluate in evaluators:
-            res = evaluate.get_extra_configs(configs)
-            extra_configs.append(res[0])
-            auxassignments.append(res[1])
-
-    # Run the correlated evaluation
-    for p in params:
-        newparms = acc.transform.deserialize(p)
-        for k in newparms:
-            wf.parameters[k] = newparms[k]
-        psi = wf.recompute(configs)[1]  # recompute gives logdet
-        rawweights = np.exp(2 * (psi - psi0))  # convert from log(|psi|) to |psi|**2
-
-        df = acc.enacc(configs, wf)
-
-        descript = {}
-        i = 0
-        for dm_type, evaluators in acc.dm_evaluators.items():
-            dms = []
-            for evaluate in evaluators:
-                dms += [evaluate(configs, wf, extra_configs[i], auxassignments[i])]
-                i += 1
-            descript.update(acc.descriptors[dm_type](dms))
-
-        for di, desc in descript.items():
-            df[di] = desc
-        df["weight"] = rawweights
-
-        data.append(df)
-
-    return data
diff --git a/pyqmc/determinant_tools.py b/pyqmc/determinant_tools.py
new file mode 100644
index 00000000..8dfc8770
--- /dev/null
+++ b/pyqmc/determinant_tools.py
@@ -0,0 +1,111 @@
+import numpy as np
+from pyscf import fci
+
+
+def binary_to_occ(S, ncore):
+    """
+  Converts the binary cistring for a given determinant
+  to occupation values for molecular orbitals within
+  the determinant.
+  """
+    occup = [int(i) for i in range(ncore)]
+    occup += [int(i + ncore) for i, c in enumerate(reversed(S)) if c == "1"]
+    max_orb = max(occup)
+    return (occup, max_orb)
+
+def determinants_from_mean_field(mf):
+    """
+    mf can be a few different things:
+    RHF on molecule or single k-point:
+        mf.mo_coeff is [nao, nmo]
+    ROHF on molecule or single k-point:
+        mf.mo_coeff is [nao, nmo]
+    UHF on molecule or single k-point
+        mf.mo_coeff is [spin][nao,nmo]
+    KRHF:
+        mf.mo_coeff is [k][nao,nmo]
+    KUHF:
+        mf.mo_coeff is [spin][k][nao,nmo]
+    """
+    detwt=np.array([1.0])
+    print(np.nonzero(mf.mo_occ[0] > 0.9))
+    occup = [[list(np.nonzero(mf.mo_occ[0] > 0.9)[0])],[list(np.nonzero(mf.mo_occ[1] > 0.9)[0])]]
+    map_dets = np.array([[0],[0]])
+    return detwt, occup, map_dets
+
+
+def interpret_ci(mc, tol):
+    """       
+    Copies over determinant coefficients and MO occupations
+    for a multi-configuration calculation mc.
+
+    This implementation separates the up and down determinants, so that we only have to compute 
+
+
+    returns:
+    detwt: array of weights for each determinant
+    occup: which orbitals go in which determinants
+    map_dets: given a determinant in detwt, which determinant in occup it corresponds to
+    """
+    from pyscf import fci
+
+    ncore = mc.ncore if hasattr(mc, "ncore") else 0
+
+    # find multi slater determinant occupation
+    if hasattr(mc, "_strs"):
+        # if this is a HCI object, it will have _strs
+        bigcis = np.abs(mc.ci) > tol
+        nstrs = int(mc._strs.shape[1] / 2)
+        # old code for single strings.
+        # deters = [(c,bin(s[0]), bin(s[1])) for c, s in zip(mc.ci[bigcis],mc._strs[bigcis,:])]
+        deters = []
+        # In pyscf, the first n/2 strings represent the up determinant and the second
+        # represent the down determinant.
+        for c, s in zip(mc.ci[bigcis], mc._strs[bigcis, :]):
+            s1 = "".join(str(bin(p)).replace("0b", "") for p in s[0:nstrs])
+            s2 = "".join(str(bin(p)).replace("0b", "") for p in s[nstrs:])
+            deters.append((c, s1, s2))
+    else:
+        deters = fci.addons.large_ci(mc.ci, mc.ncas, mc.nelecas, tol=-1)
+
+    # Create map and occupation objects
+    detwt = []
+    map_dets = [[], []]
+    occup = [[], []]
+    for x in deters:
+        if np.abs(x[0]) > tol:
+            detwt.append(x[0])
+            alpha_occ, __ = binary_to_occ(x[1], ncore)
+            beta_occ, __ = binary_to_occ(x[2], ncore)
+            if alpha_occ not in occup[0]:
+                map_dets[0].append(len(occup[0]))
+                occup[0].append(alpha_occ)
+            else:
+                map_dets[0].append(occup[0].index(alpha_occ))
+
+            if beta_occ not in occup[1]:
+                map_dets[1].append(len(occup[1]))
+                occup[1].append(beta_occ)
+            else:
+                map_dets[1].append(occup[1].index(beta_occ))
+    return np.array(detwt), occup, np.array(map_dets)
+
+
+def compute_value(updets, dndets, det_coeffs):
+    """
+    Given the up and down determinant values, safely compute the total log wave function.
+    """
+    upref = np.amax(updets[1]).real
+    dnref = np.amax(dndets[1]).real
+    phases = updets[0] * dndets[0]
+    logvals = updets[1] - upref + dndets[1] - dnref
+
+    phases = updets[0] * dndets[0]
+    wf_val = np.einsum(
+        "d,id->i", det_coeffs, phases * np.exp(logvals)
+    )
+
+    wf_sign = wf_val/np.abs(wf_val)
+    wf_logval = np.log(np.abs(wf_val)) + upref + dnref
+    return wf_sign, wf_logval
+
diff --git a/pyqmc/eval_ecp.py b/pyqmc/eval_ecp.py
index 1f07249c..78120e2d 100644
--- a/pyqmc/eval_ecp.py
+++ b/pyqmc/eval_ecp.py
@@ -34,6 +34,8 @@ def ecp_ea(mol, configs, wf, e, atom, threshold):
     masked_v_l = v_l[mask]
     masked_v_l[:, :-1] /= prob[mask, np.newaxis]
 
+    #print(np.sum(mask))
+
     # Use masked objects internally
     r_ea = r_ea[mask]
     r_ea_vec = r_ea_vec[mask]
@@ -41,13 +43,11 @@ def ecp_ea(mol, configs, wf, e, atom, threshold):
 
     # Note: epos_rot is not just apos+r_ea_i because of the boundary;
     # positions of the samples are relative to the electron, not atom.
-    epos_rot = (configs.configs[mask, e, :] - r_ea_vec)[:, np.newaxis] + r_ea_i
+    epos_rot = np.repeat(configs.configs[:,e,:][:,np.newaxis,:],P_l.shape[1], axis=1)
+    epos_rot[mask] = (configs.configs[mask, e, :] - r_ea_vec)[:, np.newaxis] + r_ea_i
 
-    # Expand externally
-    expanded_epos_rot = np.zeros((nconf, P_l.shape[1], 3))
-    expanded_epos_rot[mask] = epos_rot
-    epos = configs.make_irreducible(e, expanded_epos_rot)
-    ratio = wf.testvalue(e, epos, mask)
+    epos = configs.make_irreducible(e, epos_rot, mask)
+    ratio = wf.testvalue(e,epos, mask)
 
     # Compute local and non-local parts
     ecp_val[mask] = np.einsum("ij,ik,ijk->i", ratio, masked_v_l, P_l)
diff --git a/pyqmc/jastrowspin.py b/pyqmc/jastrowspin.py
index e13b8e3e..d41b5462 100644
--- a/pyqmc/jastrowspin.py
+++ b/pyqmc/jastrowspin.py
@@ -151,6 +151,8 @@ def _b_update(self, e, epos, mask):
 
     def _b_update_many(self, e, epos, mask, spin):
         r"""
+        Compute the update to b for each electron moving to epos. 
+
           Calculate b (e-e) partial sums for electron e
         _b_partial_e is the array :math:`B^p_{ils} = \sum_s b_l(r^i_{es}`, with e fixed; :math:`s` indexes over :math:`\uparrow` (:math:`\alpha`) and :math:`\downarrow` (:math:`\beta`) sums, not including electron e. 
           :math:`i` is the configuration index.
@@ -159,16 +161,21 @@ def _b_update_many(self, e, epos, mask, spin):
               epos: configs object for electron e
               mask: mask over configs axis, only return values for configs where mask==True. b_partial_e might have a smaller configs axis than epos, _configscurrent, and _b_partial because of the mask.
         """
+        #print(type(epos), epos.configs.shape)
         nup = self._mol.nelec[0]
-        d = epos.dist.dist_i(self._configscurrent.configs[mask], epos.configs[mask])
+        d = epos.dist.dist_i(self._configscurrent.configs[mask], epos.configs[mask,:])
+        #print("d shape", d.shape)
         r = np.linalg.norm(d, axis=-1)
+        #print("electrons", e)
         b_partial_e = np.zeros((e.shape[0], *r.shape[:-1], *self._b_partial.shape[2:]))
-
+        #print('partial', b_partial_e.shape)
         for l, b in enumerate(self.b_basis):
             bval = b.value(d, r)
+            #print("bval", bval.shape, d.shape, r.shape)
             b_partial_e[..., l, 0] = bval[..., :nup].sum(axis=-1)
             b_partial_e[..., l, 1] = bval[..., nup:].sum(axis=-1)
-            b_partial_e[..., l, spin] -= bval[..., e].T
+            #b_partial_e[..., l, spin] -= bval[..., e].T
+            b_partial_e[..., l, spin] -= np.moveaxis(bval[..., e],-1,0)
 
         return b_partial_e
 
@@ -307,7 +314,8 @@ def testvalue(self, e, epos, mask=None):
 
     def testvalue_many(self, e, epos, mask=None):
         r"""
-        Compute the ratio :math:`\Psi_{\rm new}/\Psi_{\rm old}` for moving electron e to epos.
+        Compute the ratio :math:`\Psi_{\rm new}/\Psi_{\rm old}` for moving electrons in e to epos.
+
         _avalues is the array for current configurations :math:`A_{Iks} = \sum_s a_{k}(r_{Is})` where :math:`s` indexes over :math:`\uparrow` (:math:`\alpha`) and :math:`\downarrow` (:math:`\beta`) sums.
         _bvalues is the array for current configurations :math:`B_{ls} = \sum_s b_{l}(r_{s})` where :math:`s` indexes over :math:`\uparrow\uparrow` (:math:`\alpha_1 < \alpha_2`), :math:`\uparrow\downarrow` (:math:`\alpha, \beta`), and :math:`\downarrow\downarrow` (:math:`\beta_1 < \beta_2`)  sums.
         The update for _avalues and _b_values from moving one electron only requires computing the new sum for that electron. The sums for the electron in the current configuration are stored in _a_partial and _b_partial.
diff --git a/pyqmc/multislater.py b/pyqmc/multislater.py
deleted file mode 100644
index 5ebb5e5d..00000000
--- a/pyqmc/multislater.py
+++ /dev/null
@@ -1,371 +0,0 @@
-import numpy as np
-from pyqmc.slater import sherman_morrison_row, get_complex_phase
-
-
-def sherman_morrison_ms(e, inv, vec):
-    ratio = np.einsum("idj,idj->id", vec, inv[:, :, :, e])
-    tmp = np.einsum("edk,edkj->edj", vec, inv)
-    invnew = (
-        inv
-        - np.einsum("kdi,kdj->kdij", inv[:, :, :, e], tmp)
-        / ratio[:, :, np.newaxis, np.newaxis]
-    )
-    invnew[:, :, :, e] = inv[:, :, :, e] / ratio[:, :, np.newaxis]
-    return ratio, invnew
-
-
-def binary_to_occ(S, ncore):
-    """
-  Converts the binary cistring for a given determinant
-  to occupation values for molecular orbitals within
-  the determinant.
-  """
-    occup = [int(i) for i in range(ncore)]
-    occup += [int(i + ncore) for i, c in enumerate(reversed(S)) if c == "1"]
-    max_orb = max(occup)
-    return (occup, max_orb)
-
-
-class MultiSlater:
-    """
-    A multi-determinant wave function object initialized
-    via an SCF calculation. Methods and structure are very similar
-    to the PySCFSlaterUHF class.
-
-    How to use with hci
-
-    .. code-block:: python
-
-        cisolver = pyscf.hci.SCI(mol)
-        cisolver.select_cutoff=0.1
-        nmo = mf.mo_coeff.shape[1]
-        nelec = mol.nelec
-        h1 = mf.mo_coeff.T.dot(mf.get_hcore()).dot(mf.mo_coeff)
-        h2 = pyscf.ao2mo.full(mol, mf.mo_coeff)
-        e, civec = cisolver.kernel(h1, h2, nmo, nelec, verbose=4)
-        cisolver.ci = civec[0]
-        wf = pyqmc.multislater.MultiSlater(mol, mf, cisolver, tol=0.1)
-
-
-    """
-
-    def __init__(self, mol, mf, mc, tol=None, freeze_orb=None):
-        self.tol = -1 if tol is None else tol
-        self.parameters = {}
-        self._mol = mol
-        if hasattr(mc, "nelecas"):
-            # In case nelecas overrode the information from the molecule object.
-            self._nelec = (mc.nelecas[0] + mc.ncore, mc.nelecas[1] + mc.ncore)
-        else:
-            self._nelec = mol.nelec
-        self._copy_ci(mc)
-        mo_coeff = mc.mo_coeff if hasattr(mc, "mo_coeff") else mf.mo_coeff
-        mo_cutoff_alpha = np.max(self._det_occup[0]) + 1
-        mo_cutoff_beta = np.max(self._det_occup[1]) + 1
-
-        if len(mo_coeff.shape) == 3:
-            self.parameters["mo_coeff_alpha"] = mo_coeff[0][:, :mo_cutoff_alpha]
-            self.parameters["mo_coeff_beta"] = mo_coeff[1][:, :mo_cutoff_beta]
-        else:
-            self.parameters["mo_coeff_alpha"] = mo_coeff[:, :mo_cutoff_alpha]
-            self.parameters["mo_coeff_beta"] = mo_coeff[:, :mo_cutoff_beta]
-        self._coefflookup = ("mo_coeff_alpha", "mo_coeff_beta")
-        self.pbc_str = "PBC" if hasattr(mol, "a") else ""
-        self.iscomplex = bool(sum(map(np.iscomplexobj, self.parameters.values())))
-        self.get_phase = get_complex_phase if self.iscomplex else np.sign
-        self.freeze_orb = [[], []] if freeze_orb is None else freeze_orb
-
-    def _copy_ci(self, mc):
-        """       
-        Copies over determinant coefficients and MO occupations
-        for a multi-configuration calculation mc.
-        """
-        from pyscf import fci
-
-        ncore = mc.ncore if hasattr(mc, "ncore") else 0
-
-        # find multi slater determinant occupation
-        if hasattr(mc, "_strs"):
-            # if this is a HCI object, it will have _strs
-            bigcis = np.abs(mc.ci) > self.tol
-            nstrs = int(mc._strs.shape[1] / 2)
-            # old code for single strings.
-            # deters = [(c,bin(s[0]), bin(s[1])) for c, s in zip(mc.ci[bigcis],mc._strs[bigcis,:])]
-            deters = []
-            # In pyscf, the first n/2 strings represent the up determinant and the second
-            # represent the down determinant.
-            for c, s in zip(mc.ci[bigcis], mc._strs[bigcis, :]):
-                s1 = "".join(str(bin(p)).replace("0b", "") for p in s[0:nstrs])
-                s2 = "".join(str(bin(p)).replace("0b", "") for p in s[nstrs:])
-                deters.append((c, s1, s2))
-        else:
-            deters = fci.addons.large_ci(mc.ci, mc.ncas, mc.nelecas, tol=-1)
-
-        # Create map and occupation objects
-        detwt = []
-        map_dets = [[], []]
-        occup = [[], []]
-        for x in deters:
-            if np.abs(x[0]) > self.tol:
-                detwt.append(x[0])
-                alpha_occ, __ = binary_to_occ(x[1], ncore)
-                beta_occ, __ = binary_to_occ(x[2], ncore)
-                if alpha_occ not in occup[0]:
-                    map_dets[0].append(len(occup[0]))
-                    occup[0].append(alpha_occ)
-                else:
-                    map_dets[0].append(occup[0].index(alpha_occ))
-
-                if beta_occ not in occup[1]:
-                    map_dets[1].append(len(occup[1]))
-                    occup[1].append(beta_occ)
-                else:
-                    map_dets[1].append(occup[1].index(beta_occ))
-
-        self.parameters["det_coeff"] = np.array(detwt)
-        self._det_occup = occup  # Spin, [Ndet_up_unique, Ndet_dn_unique]
-        self._det_map = np.array(map_dets)  # Spin, N_det
-
-    def recompute(self, configs):
-        """This computes the value from scratch. Returns the logarithm of the wave function as
-        (phase,logdet). If the wf is real, phase will be +/- 1."""
-
-        nconf, nelec, ndim = configs.configs.shape
-        mycoords = configs.configs.reshape((nconf * nelec, ndim))
-        ao = np.real_if_close(
-            self._mol.eval_gto(self.pbc_str + "GTOval_sph", mycoords).reshape(
-                (nconf, nelec, -1)
-            ),
-            tol=1e4,
-        )
-
-        self._aovals = ao
-        self._dets = []
-        self._inverse = []
-        for s in [0, 1]:
-            mo = ao[:, self._nelec[0] * s : self._nelec[0] + self._nelec[1] * s, :].dot(
-                self.parameters[self._coefflookup[s]]
-            )
-            mo_vals = np.swapaxes(mo[:, :, self._det_occup[s]], 1, 2)
-            self._dets.append(
-                np.array(np.linalg.slogdet(mo_vals))
-            )  # Spin, (sign, val), nconf, [ndet_up, ndet_dn]
-            self._inverse.append(
-                np.linalg.inv(mo_vals)
-            )  # spin, Nconf, [ndet_up, ndet_dn], nelec, nelec
-        return self.value()
-
-    def updateinternals(self, e, epos, mask=None):
-        """Update any internals given that electron e moved to epos. mask is a Boolean array 
-        which allows us to update only certain walkers"""
-        # MAY want to vectorize later if it really hangs here, shouldn't!
-
-        s = int(e >= self._nelec[0])
-        if mask is None:
-            mask = [True] * epos.configs.shape[0]
-        eeff = e - s * self._nelec[0]
-        ao = np.real_if_close(
-            self._mol.eval_gto(self.pbc_str + "GTOval_sph", epos.configs[mask]), tol=1e4
-        )
-        self._aovals[mask, e, :] = ao
-        mo = ao.dot(self.parameters[self._coefflookup[s]])
-
-        mo_vals = mo[:, self._det_occup[s]]
-        det_ratio, self._inverse[s][mask, :, :, :] = sherman_morrison_ms(
-            eeff, self._inverse[s][mask, :, :, :], mo_vals
-        )
-
-        self._updateval(det_ratio, s, mask)
-
-    def value(self):
-        """Return logarithm of the wave function as noted in recompute()"""
-        updets = self._dets[0][:, :, self._det_map[0]]
-        dndets = self._dets[1][:, :, self._det_map[1]]
-        upref = np.amax(self._dets[0][1])
-        dnref = np.amax(self._dets[1][1])
-        phases = updets[0] * dndets[0]
-        logvals = updets[1] - upref + dndets[1] - dnref
-
-        wf_val = np.einsum(
-            "d,id->i", self.parameters["det_coeff"], phases * np.exp(logvals)
-        )
-
-        wf_sign = self.get_phase(wf_val)
-        wf_logval = np.log(np.abs(wf_val)) + upref + dnref
-        return wf_sign, wf_logval
-
-    def _updateval(self, ratio, s, mask):
-        self._dets[s][0, mask, :] *= self.get_phase(ratio)
-        self._dets[s][1, mask, :] += np.log(np.abs(ratio))
-
-    def _testrow(self, e, vec, mask=None, spin=None):
-        """vec is a nconfig,nmo vector which replaces row e"""
-        s = int(e >= self._nelec[0]) if spin is None else spin
-        if mask is None:
-            mask = [True] * vec.shape[0]
-
-        ratios = np.einsum(
-            "i...dj,idj...->i...d",
-            vec,
-            self._inverse[s][mask][..., e - s * self._nelec[0]],
-        )
-        det_array = (
-            self._dets[0][0, :, self._det_map[0]][:, mask]
-            * self._dets[1][0, :, self._det_map[1]][:, mask]
-            * np.exp(
-                self._dets[0][1, :, self._det_map[0]][:, mask]
-                + self._dets[1][1, :, self._det_map[1]][:, mask]
-            )
-        )
-        numer = np.einsum(
-            "i...d,d,di->i...",
-            ratios[..., self._det_map[s]],
-            self.parameters["det_coeff"],
-            det_array,
-        )
-
-        curr_val = self.value()
-        denom = curr_val[0][mask] * np.exp(curr_val[1][mask])
-        if len(numer.shape) == 2:
-            denom = denom[:, np.newaxis]
-        return numer / denom
-
-    def _testcol(self, det, i, s, vec):
-        """vec is a nconfig,nmo vector which replaces column i 
-        of spin s in determinant det"""
-
-        return np.einsum(
-            "ij...,ij->i...", vec, self._inverse[s][:, det, i, :], optimize="greedy"
-        )
-
-    def gradient(self, e, epos):
-        """ Compute the gradient of the log wave function 
-        Note that this can be called even if the internals have not been updated for electron e,
-        if epos differs from the current position of electron e."""
-        s = int(e >= self._nelec[0])
-        aograd = np.real_if_close(
-            self._mol.eval_gto("GTOval_sph_deriv1", epos.configs), tol=1e4
-        )
-        mograd = aograd.dot(self.parameters[self._coefflookup[s]])
-        mograd_vals = mograd[:, :, self._det_occup[s]]
-
-        ratios = np.asarray([self._testrow(e, x) for x in mograd_vals])
-        return ratios[1:] / ratios[:1]
-
-    def laplacian(self, e, epos):
-        """ Compute the laplacian Psi/ Psi. """
-        s = int(e >= self._nelec[0])
-        ao = np.real_if_close(
-            self._mol.eval_gto(self.pbc_str + "GTOval_sph_deriv2", epos.configs)[
-                [0, 4, 7, 9]
-            ],
-            tol=1e4,
-        )
-        molap = np.dot(
-            [ao[0], ao[1:].sum(axis=0)], self.parameters[self._coefflookup[s]]
-        )
-        molap_vals = self._testrow(e, molap[1][:, self._det_occup[s]])
-        testvalue = self._testrow(e, molap[0][:, self._det_occup[s]])
-
-        return molap_vals / testvalue
-
-    def gradient_laplacian(self, e, epos):
-        s = int(e >= self._nelec[0])
-        ao = np.real_if_close(
-            self._mol.eval_gto(self.pbc_str + "GTOval_sph_deriv2", epos.configs)[
-                [0, 1, 2, 3, 4, 7, 9]
-            ],
-            tol=1e4,
-        )
-        ao = np.concatenate([ao[0:4], ao[4:].sum(axis=0, keepdims=True)])
-        mo = np.dot(ao, self.parameters[self._coefflookup[s]])
-        mo_vals = mo[:, :, self._det_occup[s]]
-        ratios = np.asarray([self._testrow(e, x) for x in mo_vals])
-        return ratios[1:-1] / ratios[:1], ratios[-1] / ratios[0]
-
-    def testvalue(self, e, epos, mask=None):
-        """ return the ratio between the current wave function and the wave function if 
-        electron e's position is replaced by epos"""
-        s = int(e >= self._nelec[0])
-        if mask is None:
-            mask = [True] * epos.configs.shape[0]
-        eposmask = epos.configs[mask]
-        if len(eposmask) == 0:
-            return np.zeros(eposmask.shape[:2])
-        ao = self._mol.eval_gto(
-            self.pbc_str + "GTOval_sph", eposmask.reshape((-1, 3))
-        ).reshape((*eposmask.shape[:-1], -1))
-        mo = ao.dot(self.parameters[self._coefflookup[s]])
-        mo_vals = mo[..., self._det_occup[s]]
-        return self._testrow(e, mo_vals, mask)
-
-    def testvalue_many(self, e, epos, mask=None):
-        """ return the ratio between the current wave function and the wave function if 
-        electron e's position is replaced by epos for each electron"""
-        s = (e >= self._nelec[0]).astype(int)
-        if mask is None:
-            mask = [True] * epos.configs.shape[0]
-        eposmask = epos.configs[mask]
-        if len(eposmask) == 0:
-            return np.zeros(eposmask.shape[:2])
-        ao = self._mol.eval_gto(
-            self.pbc_str + "GTOval_sph", eposmask.reshape((-1, 3))
-        ).reshape((*eposmask.shape[:-1], -1))
-
-        ratios = np.zeros((epos.configs.shape[0], e.shape[0]))
-        for spin in [0, 1]:
-            ind = s == spin
-            mo = ao.dot(self.parameters[self._coefflookup[spin]])
-            mo_vals = mo[..., self._det_occup[spin]]
-            ratios[:, ind] = self._testrow(e[ind], mo_vals, mask, spin=spin)
-
-        return ratios
-
-    def pgradient(self):
-        r"""Compute the parameter gradient of Psi. 
-        Returns $$d_p \Psi/\Psi$$ as a dictionary of numpy arrays,
-        which correspond to the parameter dictionary."""
-        d = {}
-
-        # Det coeff
-        det_coeff_grad = (
-            self._dets[0][0, :, self._det_map[0]]
-            * self._dets[1][0, :, self._det_map[1]]
-            * np.exp(
-                self._dets[0][1, :, self._det_map[0]]
-                + self._dets[1][1, :, self._det_map[1]]
-            )
-        )
-
-        curr_val = self.value()
-        d["det_coeff"] = (
-            det_coeff_grad.T / (curr_val[0] * np.exp(curr_val[1]))[:, np.newaxis]
-        )
-
-        # Mo_coeff, adapted from SlaterUHF
-        for parm in ["mo_coeff_alpha", "mo_coeff_beta"]:
-            s = 0
-            if "beta" in parm:
-                s = 1
-
-            ao = self._aovals[
-                :, s * self._nelec[0] : self._nelec[s] + s * self._nelec[0], :
-            ]
-            pgrad_shape = (ao.shape[0],) + self.parameters[parm].shape
-            pgrad = np.zeros(pgrad_shape)
-
-            largest_mo = np.max(np.ravel(self._det_occup[s]))
-            for i in range(largest_mo + 1):  # MO loop
-                if i not in self.freeze_orb[s]:
-                    for det in range(self.parameters["det_coeff"].shape[0]):  # Det loop
-                        if (
-                            i in self._det_occup[s][self._det_map[s][det]]
-                        ):  # Check if MO in det
-                            col = self._det_occup[s][self._det_map[s][det]].index(i)
-                            pgrad[:, :, i] += (
-                                self.parameters["det_coeff"][det]
-                                * d["det_coeff"][:, det, np.newaxis]
-                                * self._testcol(self._det_map[s][det], col, s, ao)
-                            )
-            d[parm] = np.array(pgrad)
-        return d
diff --git a/pyqmc/multislaterpbc.py b/pyqmc/multislaterpbc.py
deleted file mode 100644
index f27e4aaa..00000000
--- a/pyqmc/multislaterpbc.py
+++ /dev/null
@@ -1,325 +0,0 @@
-import numpy as np
-from pyqmc.multislater import sherman_morrison_ms
-from pyqmc.slater import get_k_indices
-from pyqmc.supercell import get_supercell_kpts
-from pyqmc import pbc, slater
-
-
-class MultiSlaterPBC:
-    """
-    A multi-determinant wave function object initialized
-    via an SCF calculation. Methods and structure are very similar
-    to the PySCFSlaterUHF class.
-    """
-
-    def __init__(
-        self, supercell, mf, tol=-1, detwt=(1,), occup=None, map_dets=None, twist=None
-    ):
-        """
-        detwt: list of determinant weights
-        occup: list (spin, det, dict{kind: occupation list})
-        map_dets: list (spin, ndet) to identify which determinant of each spin to use (e.g. may use the same up-determinant in multiple products)
-        """
-        for attribute in ["original_cell", "S"]:
-            if not hasattr(supercell, attribute):
-                print('Warning: supercell is missing attribute "%s"' % attribute)
-                print("setting original_cell=supercell and S=np.eye(3)")
-                supercell.original_cell = supercell
-                supercell.S = np.eye(3)
-
-        assert occup is not None
-        assert len(map_dets[0]) == len(detwt) and len(map_dets[1]) == len(detwt)
-        self.parameters = {"det_coeff": np.array(detwt)}
-
-        self.tol = tol
-        self.real_tol = 1e4
-        self._mol = supercell.original_cell
-        self._nelec = tuple(int(sum(len(v) for v in o[0].values())) for o in occup)
-
-        self.supercell = supercell
-        if twist is None:
-            twist = np.zeros(3)
-        else:
-            twist = np.dot(np.linalg.inv(supercell.a), np.mod(twist, 1.0)) * 2 * np.pi
-        self.kinds = get_k_indices(self._mol, mf, get_supercell_kpts(supercell) + twist)
-        self._kpts = mf.kpts[self.kinds]
-        assert len(self.kinds) == len(self._kpts), (self._kpts, mf.kpts)
-        self.nk = len(self._kpts)
-        print("nk", self.nk, self.kinds)
-
-        maxorb = [
-            {
-                k: np.amax([np.amax(list(d[k]) + [-1]) for d in o]) + 1
-                for k in self.kinds
-            }
-            for o in occup
-        ]
-
-        self.param_split = {}
-        self._coefflookup = ("mo_coeff_alpha", "mo_coeff_beta")
-        for s, lookup in enumerate(self._coefflookup):
-            mclist = []
-            for kind in self.kinds:
-                if len(mf.mo_coeff[0][0].shape) == 2:
-                    mca = mf.mo_coeff[s][kind][:, : maxorb[s][kind]]
-                else:
-                    mca = mf.mo_coeff[kind][:, : maxorb[s][kind]]
-                mca = np.real_if_close(mca, tol=self.real_tol)
-                mclist.append(mca / np.sqrt(self.nk))
-            self.param_split[lookup] = np.cumsum([m.shape[1] for m in mclist])
-            self.parameters[lookup] = np.concatenate(mclist, axis=-1)
-
-        # _det_occup: Spin, [(Ndet_up_unique, nup), (Ndet_dn_unique, ndn)]
-        self._det_occup = [
-            np.zeros((len(o), n), dtype=int) for o, n in zip(occup, self._nelec)
-        ]
-        for s, o in enumerate(occup):
-            orb_inds = np.cumsum([0] + [maxorb[s][k] for k in self.kinds[:-1]])
-            for d, det in enumerate(o):
-                self._det_occup[s][d] = np.concatenate(
-                    [np.array(det[k]) + ind for k, ind in zip(self.kinds, orb_inds)]
-                )
-        self._det_map = np.asarray(map_dets)  # Spin, N_det
-        print("det_coeff", self.parameters["det_coeff"])
-        print("_det_occup\n", self._det_occup)
-        print("_det_map\n", self._det_map)
-
-        self.iscomplex = bool(sum(map(np.iscomplexobj, self.parameters.values())))
-        self.iscomplex = self.iscomplex or np.linalg.norm(self._kpts) > 1e-12
-        self.dtype = complex if self.iscomplex else float
-        if self.iscomplex:
-            self.get_phase = slater.get_complex_phase
-            self.get_wrapphase = slater.get_wrapphase_complex
-        else:
-            self.get_phase = np.sign
-            self.get_wrapphase = slater.get_wrapphase_real
-
-    def evaluate_orbitals(self, configs, mask=None, eval_str="PBCGTOval_sph"):
-        mycoords = configs.configs
-        configswrap = configs.wrap
-        if mask is not None:
-            mycoords = mycoords[mask]
-            configswrap = configswrap[mask]
-        mycoords = mycoords.reshape((-1, mycoords.shape[-1]))
-        # wrap supercell positions into primitive cell
-        prim_coords, prim_wrap = pbc.enforce_pbc(self._mol.lattice_vectors(), mycoords)
-        configswrap = configswrap.reshape(prim_wrap.shape)
-        wrap = prim_wrap + np.dot(configswrap, self.supercell.S)
-        kdotR = np.linalg.multi_dot((self._kpts, self._mol.lattice_vectors().T, wrap.T))
-        wrap_phase = self.get_wrapphase(kdotR)
-        # evaluate AOs for all electron positions
-        ao = self._mol.eval_gto(eval_str, prim_coords, kpts=self._kpts)
-        ao = [ao[k] * wrap_phase[k][:, np.newaxis] for k in range(self.nk)]
-        return ao
-
-    def evaluate_mos(self, aos, s):
-        l = self._coefflookup[s]
-        p = np.split(self.parameters[l], self.param_split[l], axis=-1)
-        mo = [np.dot(ao, p[k]) for k, ao in enumerate(aos)]
-        mo = np.concatenate(mo, axis=-1)
-        return mo[..., self._det_occup[s]]
-
-    def recompute(self, configs):
-        """This computes the value from scratch. Returns the logarithm of the wave function as
-        (phase,logdet). If the wf is real, phase will be +/- 1."""
-        nconf, nelec, ndim = configs.configs.shape
-        aos = self.evaluate_orbitals(configs)
-        aos = np.reshape(aos, (self.nk, nconf, nelec, -1))
-        self._aovals = aos
-        self._dets = []
-        self._inverse = []
-        for s in [0, 1]:
-            i0, i1 = s * self._nelec[0], self._nelec[0] + s * self._nelec[1]
-            mo_vals = self.evaluate_mos(aos[:, :, i0:i1], s).swapaxes(1, 2)
-            self._dets.append(np.array(np.linalg.slogdet(mo_vals)))
-            self._inverse.append(np.linalg.inv(mo_vals))
-        return self.value()
-
-    def updateinternals(self, e, epos, mask=None):
-        """Update any internals given that electron e moved to epos. mask is a Boolean array 
-        which allows us to update only certain walkers"""
-        s = int(e >= self._nelec[0])
-        if mask is None:
-            mask = [True] * epos.configs.shape[0]
-        eeff = e - s * self._nelec[0]
-        aos = self.evaluate_orbitals(epos, mask=mask)
-        self._aovals[:, mask, e, :] = np.asarray(aos)
-        mo = self.evaluate_mos(aos, s)
-        ratio, self._inverse[s][mask] = sherman_morrison_ms(
-            eeff, self._inverse[s][mask], mo
-        )
-        self._updateval(ratio, s, mask)
-
-    # identical to multislater.py
-    def value(self):
-        """Return logarithm of the wave function as noted in recompute()"""
-        wf_val = 0
-        wf_sign = 0
-
-        wf_val = np.einsum(
-            "d,di->i",
-            self.parameters["det_coeff"],
-            self._dets[0][0, :, self._det_map[0]]
-            * self._dets[1][0, :, self._det_map[1]]
-            * np.exp(
-                self._dets[0][1, :, self._det_map[0]]
-                + self._dets[1][1, :, self._det_map[1]]
-            ),
-        )
-
-        wf_sign = self.get_phase(wf_val)
-        wf_val = np.log(np.abs(wf_val))
-        return wf_sign, wf_val
-
-    # identical to slateruhf
-    def _updateval(self, ratio, s, mask):
-        self._dets[s][0, mask, :] *= self.get_phase(ratio)
-        self._dets[s][1, mask, :] += np.log(np.abs(ratio))
-
-    # identical to multislater.py
-    def _testrow(self, e, vec, mask=None, spin=None):
-        """vec is a nconfig,nmo vector which replaces row e"""
-        s = int(e >= self._nelec[0]) if spin is None else spin
-        if mask is None:
-            mask = [True] * vec.shape[0]
-
-        ratios = np.einsum(
-            "i...dj,idj...->i...d",
-            vec,
-            self._inverse[s][mask][..., e - s * self._nelec[0]],
-        )
-        det_array = (
-            self._dets[0][0, :, self._det_map[0]][:, mask]
-            * self._dets[1][0, :, self._det_map[1]][:, mask]
-            * np.exp(
-                self._dets[0][1, :, self._det_map[0]][:, mask]
-                + self._dets[1][1, :, self._det_map[1]][:, mask]
-            )
-        )
-        numer = np.einsum(
-            "i...d,d,di->i...",
-            ratios[..., self._det_map[s]],
-            self.parameters["det_coeff"],
-            det_array,
-        )
-
-        curr_val = self.value()
-        denom = curr_val[0][mask] * np.exp(curr_val[1][mask])
-        if len(numer.shape) == 2:
-            denom = denom[:, np.newaxis]
-        return numer / denom
-
-    # identical to multislater.py
-    def _testcol(self, det, i, s, vec):
-        """vec is a nconfig,nmo vector which replaces column i 
-        of spin s in determinant det"""
-
-        return np.einsum("ij...,ij->i...", vec, self._inverse[s][:, det, i, :])
-
-    # identical to slaterpbc
-    def testvalue(self, e, epos, mask=None):
-        """ return the ratio between the current wave function and the wave function if 
-        electron e's position is replaced by epos"""
-        s = int(e >= self._nelec[0])
-        if mask is None:
-            mask = [True] * epos.configs.shape[0]
-        nmask = np.sum(mask)
-        if nmask == 0:
-            return np.zeros((0, epos.configs.shape[1]))
-        ao = self.evaluate_orbitals(epos, mask=mask)
-        mo = self.evaluate_mos(ao, s)
-        mo = mo.reshape((nmask, *epos.configs.shape[1:-1], -1, self._nelec[s]))
-        return self._testrow(e, mo, mask)
-
-    # identical to slaterpbc
-    def testvalue_many(self, e, epos, mask=None):
-        """ return the ratio between the current wave function and the wave function if 
-        electron e's position is replaced by epos for each electron"""
-        s = (e >= self._nelec[0]).astype(int)
-        if mask is None:
-            mask = [True] * epos.configs.shape[0]
-        nmask = np.sum(mask)
-        if nmask == 0:
-            return np.zeros((0, epos.configs.shape[1]))
-
-        ao = self.evaluate_orbitals(epos, mask=mask)
-        ratios = np.zeros((epos.configs.shape[0], e.shape[0]), dtype=self.dtype)
-        for spin in [0, 1]:
-            ind = s == spin
-            mo = self.evaluate_mos(ao, spin)
-            mo = mo.reshape(nmask, *epos.configs.shape[1:-1], -1, self._nelec[spin])
-            ratios[:, ind] = self._testrow(e[ind], mo, mask=mask, spin=spin)
-        return ratios
-
-    # identical to slaterpbc
-    def gradient(self, e, epos):
-        """ Compute the gradient of the log wave function 
-        Note that this can be called even if the internals have not been updated for electron e,
-        if epos differs from the current position of electron e."""
-        s = int(e >= self._nelec[0])
-        aograd = self.evaluate_orbitals(epos, eval_str="PBCGTOval_sph_deriv1")
-        mograd = self.evaluate_mos(aograd, s)
-        ratios = np.asarray([self._testrow(e, x) for x in mograd])
-        return ratios[1:] / ratios[:1]
-
-    # identical to slaterpbc
-    def laplacian(self, e, epos):
-        """ Compute the laplacian Psi/ Psi. """
-        s = int(e >= self._nelec[0])
-        ao = self.evaluate_orbitals(epos, eval_str="PBCGTOval_sph_deriv2")
-        aostack = [np.stack([ak[0], ak[[4, 7, 9]].sum(axis=0)], axis=0) for ak in ao]
-        molap = self.evaluate_mos(aostack, s)
-        ratios = np.asarray([self._testrow(e, x) for x in molap])
-        return ratios[1] / ratios[0]
-
-    # identical to slaterpbc
-    def gradient_laplacian(self, e, epos):
-        s = int(e >= self._nelec[0])
-        ao = self.evaluate_orbitals(epos, eval_str="PBCGTOval_sph_deriv2")
-        aostack = [
-            np.concatenate([ak[0:4], ak[[4, 7, 9]].sum(axis=0, keepdims=True)], axis=0)
-            for ak in ao
-        ]
-        mo_vals = self.evaluate_mos(aostack, s)
-        ratios = np.asarray([self._testrow(e, x) for x in mo_vals])
-        return ratios[1:-1] / ratios[:1], ratios[-1] / ratios[0]
-
-    def pgradient(self):
-        r"""Compute the parameter gradient of Psi. 
-        Returns :math:`\frac{\partial_p \Psi}{\Psi}` as a dictionary of numpy arrays,
-        which correspond to the parameter dictionary."""
-        d = {}
-        # Det coeff
-        det_coeff_grad = (
-            self._dets[0][0, :, self._det_map[0]]
-            * self._dets[1][0, :, self._det_map[1]]
-            * np.exp(
-                self._dets[0][1, :, self._det_map[0]]
-                + self._dets[1][1, :, self._det_map[1]]
-            )
-        )
-
-        curr_val = self.value()
-        d["det_coeff"] = (
-            det_coeff_grad.T / (curr_val[0] * np.exp(curr_val[1]))[:, np.newaxis]
-        )
-
-        for s, parm in enumerate(["mo_coeff_alpha", "mo_coeff_beta"]):
-            i0, i1 = s * self._nelec[0], self._nelec[0] + s * self._nelec[1]
-            ao = self._aovals[:, :, i0:i1, :]  # nk, nconf, nelec, nao
-            pgrad_shape = (ao.shape[1],) + self.parameters[parm].shape
-            pgrad = np.zeros(pgrad_shape, dtype=self.dtype)
-            split_sizes = np.diff([0] + list(self.param_split[parm]))
-            k = np.repeat(np.arange(self.nk), split_sizes)
-
-            largest_mo = np.max(np.ravel(self._det_occup[s]))
-            for det, detwt in enumerate(self.parameters["det_coeff"]):  # Det loop
-                mapdet = self._det_map[s][det]
-                for col, i in enumerate(self._det_occup[s][mapdet]):
-                    pgrad[:, :, i] += (
-                        detwt
-                        * d["det_coeff"][:, det, np.newaxis]
-                        * self._testcol(mapdet, col, s, ao[k[i]])
-                    )
-            d[parm] = np.array(pgrad)
-        return d
diff --git a/pyqmc/obdm.py b/pyqmc/obdm.py
index b6c92ad8..e0fec3fc 100644
--- a/pyqmc/obdm.py
+++ b/pyqmc/obdm.py
@@ -1,8 +1,10 @@
 """ Evaluate the OBDM for a wave function object. """
+from pyqmc.orbitals import MoleculeOrbitalEvaluator, PBCOrbitalEvaluatorKpoints
 import numpy as np
-from copy import deepcopy
 from pyqmc.mc import initial_guess
-from pyqmc import pbc, slater
+
+import pyqmc.supercell as supercell
+
 
 
 class OBDMAccumulator:
@@ -70,47 +72,28 @@ def __init__(
             self._electrons = np.arange(0, np.sum(mol.nelec))
 
         self.iscomplex = bool(sum(map(np.iscomplexobj, orb_coeff)))
-        if hasattr(mol, "lattice_vectors"):
-            assert kpts is not None
-            assert len(orb_coeff) == len(kpts)
-            self.iscomplex = self.iscomplex or np.linalg.norm(kpts) > 1e-12
-            self._kpts = kpts
-            self.evaluate_orbitals = self._evaluate_orbitals_pbc
-            if self.iscomplex:
-                self.get_wrapphase = slater.get_wrapphase_complex
-            else:
-                self.get_wrapphase = slater.get_wrapphase_real
-            for attribute in ["original_cell", "S"]:
-                if not hasattr(mol, attribute):
-                    from pyqmc.supercell import get_supercell
-
-                    mol = get_supercell(mol, np.eye(3))
-            self.supercell = mol
-            self._mol = mol.original_cell
-        else:
-            self._mol = mol
-            self.evaluate_orbitals = self._evaluate_orbitals
-            assert (
-                len(orb_coeff.shape) == 2
-            ), "orb_coeff should be a list of orbital coefficients."
-
-        self._orb_coeff = orb_coeff
-        if hasattr(self._orb_coeff, "shape"):
-            self.norb = self._orb_coeff.shape[1]
-        else:
-            self.norb = sum(c.shape[1] for c in self._orb_coeff)
+
+        if kpts is None:
+            self.orbitals = MoleculeOrbitalEvaluator(mol, [orb_coeff, orb_coeff])
+        else: 
+            if not hasattr(mol, "original_cell"):
+                mol = supercell.get_supercell(mol, np.eye(3))
+            self.orbitals = PBCOrbitalEvaluatorKpoints(mol, [orb_coeff, orb_coeff], kpts)
+
+
         self._tstep = tstep
         self.nelec = len(self._electrons)
         self._nsweeps = nsweeps
         self._nstep = nsweeps * self.nelec
+        self.norb = self.orbitals.parameters['mo_coeff_alpha'].shape[-1]
 
         self._extra_config = initial_guess(mol, int(naux / self.nelec) + 1)
-        self._extra_config.reshape((-1, 3))
+        self._extra_config.reshape((-1,1,3))
 
-        accept, extra_configs = self.sample_onebody(self._extra_config, warmup)
+        accept, extra_configs, _ = sample_onebody(self._extra_config, self.orbitals, spin=0, nsamples=warmup, tstep=self._tstep)
         self._extra_config = extra_configs[-1]
 
-    def __call__(self, configs, wf, extra_configs=None, auxassignments=None):
+    def __call__(self, configs, wf):
         """ 
         """
 
@@ -121,48 +104,40 @@ def __call__(self, configs, wf, extra_configs=None, auxassignments=None):
             "norm": np.zeros((nconf, self.norb)),
             "acceptance": np.zeros(nconf),
         }
-        acceptance = 0
         naux = self._extra_config.configs.shape[0]
 
-        if extra_configs is None:
-            auxassignments = np.random.randint(0, naux, size=(self._nsweeps, nconf))
-            accept, extra_configs = self.sample_onebody(
-                self._extra_config, self._nsweeps
-            )
-            self._extra_config = extra_configs[-1]
-            results["acceptance"] += np.sum(accept) / naux
-        else:
-            assert auxassignments is not None
+        auxassignments = np.random.randint(0, naux, size=(self._nsweeps, nconf))
+        accept, extra_configs, borb_aux = sample_onebody(
+            self._extra_config, self.orbitals, spin=0, nsamples=self._nsweeps, tstep=self._tstep
+        )
+        self._extra_config = extra_configs[-1]
+
+        for conf,assign in zip(extra_configs,auxassignments):
+            conf.resample(assign)
+        borb_aux = [orb[assign,...] for orb, assign in zip(borb_aux,auxassignments)]
+
+        results["acceptance"] += np.sum(accept) / naux
 
         borb_configs = self.evaluate_orbitals(configs.electron(self._electrons))
         borb_configs = borb_configs.reshape(nconf, self.nelec, -1)
-        # Orbital evaluations at extra coordinate.
-        all_extra_configs = extra_configs[0].mask(
-            np.zeros(naux * self._nsweeps, dtype=int)
-        )
-        all_extra_configs.join(extra_configs)
-        borb_aux = self.evaluate_orbitals(all_extra_configs)
-        borb_aux = borb_aux.reshape(self._nsweeps, naux, -1)
+
         bauxsquared = np.abs(borb_aux) ** 2
         fsum = np.sum(bauxsquared, axis=-1, keepdims=True)
         norm = bauxsquared / fsum
         baux_f = borb_aux / fsum
 
         for sweep, aux in enumerate(auxassignments):
-            epos = extra_configs[sweep].configs[aux]
-            newconfigs = configs.make_irreducible(0, epos)
-            wfratio = wf.testvalue_many(self._electrons, newconfigs)
-
+            wfratio = wf.testvalue_many(self._electrons, extra_configs[sweep].electron(0))
             ratio = np.einsum(
                 "ie,ij,iek->ijk",
                 wfratio.conj(),
-                baux_f[sweep, aux],
+                baux_f[sweep, :],
                 borb_configs.conj(),
                 optimize=True,
             )
 
             results["value"] += ratio
-            results["norm"] += norm[sweep, aux]
+            results["norm"] += norm[sweep, :]
 
         results["value"] /= self._nstep
         results["norm"] = results["norm"] / self._nstep
@@ -173,61 +148,9 @@ def __call__(self, configs, wf, extra_configs=None, auxassignments=None):
     def avg(self, configs, wf):
         return {k: np.mean(it, axis=0) for k, it in self(configs, wf).items()}
 
-    def get_extra_configs(self, configs):
-        """ Returns an nstep length array of configurations
-            starting from self._extra_config """
-        nconf = configs.configs.shape[0]
-        naux = self._extra_config.configs.shape[0]
-        auxassignments = np.random.randint(0, naux, size=(self._nsweeps, nconf))
-        accept, extra_configs = self.sample_onebody(self._extra_config, self._nsweeps)
-        self._extra_config = extra_configs[-1]
-        return extra_configs, auxassignments
-
-    def sample_onebody(self, configs, nsamples=1):
-        r""" For a set of orbitals defined by orb_coeff, return samples from :math:`f(r) = \sum_i \phi_i(r)^2`. """
-        n = configs.configs.shape[0]
-        borb = self.evaluate_orbitals(configs)
-        fsum = (np.abs(borb) ** 2).sum(axis=1)
-
-        allaccept = np.zeros((nsamples, n))
-        allconfigs = []
-        for s in range(nsamples):
-            shift = np.sqrt(self._tstep) * np.random.randn(*configs.configs.shape)
-            newconfigs = configs.make_irreducible(None, configs.configs + shift)
-            borbnew = self.evaluate_orbitals(newconfigs)
-            fsumnew = (np.abs(borbnew) ** 2).sum(axis=1)
-            accept = fsumnew / fsum > np.random.rand(n)
-            configs.move_all(newconfigs, accept)
-            borb[accept] = borbnew[accept]
-            fsum[accept] = fsumnew[accept]
-            allconfigs.append(configs.copy())
-            allaccept[s] = accept
-
-        # config_pack = configs.mask(np.tile(np.arange(n), 2)) # two copies of configs
-        # config_pack.join([configs, newconfigs])
-
-        # borb = self.evaluate_orbitals(config_pack)
-
-        return allaccept, allconfigs
-
-    def _evaluate_orbitals(self, configs):
-        ao = self._mol.eval_gto("GTOval_sph", configs.configs.reshape(-1, 3))
-        return ao.dot(self._orb_coeff)
-
-    def _evaluate_orbitals_pbc(self, configs):
-        # wrap supercell positions into primitive cell
-        lvecs = self._mol.lattice_vectors()
-        prim_coords, prim_wrap = pbc.enforce_pbc(lvecs, configs.configs)
-        wrap = prim_wrap + np.dot(configs.wrap, self.supercell.S)
-        wrap = wrap.reshape(-1, 3)
-        prim_coords = prim_coords.reshape(-1, 3)
-        kdotR = np.linalg.multi_dot((self._kpts, lvecs.T, wrap.T))
-        wrap_phase = self.get_wrapphase(kdotR)
-        # evaluate AOs for all electron positions
-        ao = self._mol.eval_gto("PBCGTOval_sph", prim_coords, kpts=self._kpts)
-        ao = [aok * wrap_phase[k][:, np.newaxis] for k, aok in enumerate(ao)]
-        borb = [aok.dot(ock) for aok, ock in zip(ao, self._orb_coeff)]
-        return np.concatenate(borb, axis=1)
+    def evaluate_orbitals(self, configs):
+        ao = self.orbitals.aos("GTOval_sph",configs)
+        return self.orbitals.mos(ao, spin=0)
 
     def keys(self):
         return set(["value", "norm", "acceptance"])
@@ -237,19 +160,33 @@ def shapes(self):
         return {"value": (norb, norb), "norm": (norb,), "acceptance": ()}
 
 
-def sample_onebody(mol, orb_coeff, configs, tstep=2.0):
-    r""" For a set of orbitals defined by orb_coeff, return samples from :math:`f(r) = \sum_i \phi_i(r)^2`. """
-    n = configs.shape[0]
-    config_pack = np.concatenate([configs, configs], axis=0)
-    config_pack[n:] += np.sqrt(tstep) * np.random.randn(*configs.shape)
 
-    ao = mol.eval_gto("GTOval_sph", config_pack)
-    borb = ao.dot(orb_coeff)
+
+def sample_onebody(configs, orbitals, spin, nsamples=1, tstep=0.5):
+    r""" For a set of orbitals defined by orb_coeff, return samples from :math:`f(r) = \sum_i \phi_i(r)^2`. """
+    n = configs.configs.shape[0]
+    ao = orbitals.aos("GTOval_sph",configs)
+    borb= orbitals.mos(ao, spin=spin)
     fsum = (np.abs(borb) ** 2).sum(axis=1)
 
-    accept = fsum[n:] / fsum[0:n] > np.random.rand(n)
-    configs[accept] = config_pack[n:][accept]
-    return accept, configs
+    allaccept = np.zeros((nsamples, n))
+    allconfigs = []
+    allorbs = []
+    for s in range(nsamples):
+        shift = np.sqrt(tstep) * np.random.randn(*configs.configs.shape)
+        newconfigs = configs.make_irreducible(0, configs.configs + shift)
+        ao = orbitals.aos("GTOval_sph", newconfigs)
+        borbnew = orbitals.mos(ao, spin=spin)
+        fsumnew = (np.abs(borbnew) ** 2).sum(axis=1)
+        accept = fsumnew / fsum > np.random.rand(n)
+        configs.move_all(newconfigs, accept)
+        borb[accept] = borbnew[accept]
+        fsum[accept] = fsumnew[accept]
+        allconfigs.append(configs.copy())
+        allaccept[s] = accept
+        allorbs.append(borb.copy())
+
+    return allaccept, allconfigs, allorbs
 
 
 def normalize_obdm(obdm, norm):
diff --git a/pyqmc/orbitals.py b/pyqmc/orbitals.py
new file mode 100644
index 00000000..a68e89eb
--- /dev/null
+++ b/pyqmc/orbitals.py
@@ -0,0 +1,227 @@
+import numpy as np
+import pyqmc.pbc as pbc
+from pyqmc.supercell import get_supercell_kpts, get_supercell
+import pyqmc.determinant_tools
+
+"""
+The evaluators have the concept of a 'set' of atomic orbitals, that may apply to 
+different sets of molecular orbitals
+
+For example, for the PBC evaluator, each k-point is a set, since each molecular 
+orbital is only a sum over the k-point of its type.
+
+In the future, this could apply to orbitals of a given point group symmetry, for example.
+"""
+
+
+def get_wrapphase_real(x):
+    return (-1) ** np.round(x / np.pi)
+
+def get_wrapphase_complex(x):
+    return np.exp(1j * x)
+
+def get_complex_phase(x):
+    return x / np.abs(x)
+
+def choose_evaluator_from_pyscf(mol, mf, mc=None, twist=None):
+    """
+    Returns:
+    an orbital evaluator chosen based on the inputs. 
+    """
+
+    if hasattr(mol, "a"): 
+        if mc is not None:
+            raise NotImplementedError("Do not support multiple determinants for k-points orbital evaluator")
+        return PBCOrbitalEvaluatorKpoints.from_mean_field(mol, mf, twist)
+    if mc is None:
+        return MoleculeOrbitalEvaluator.from_pyscf(mol, mf)
+    return MoleculeOrbitalEvaluator.from_pyscf(mol, mf, mc)
+
+
+
+
+class MoleculeOrbitalEvaluator:
+    def __init__(self, mol, mo_coeff):
+        self.iscomplex=False
+        self.parameters={'mo_coeff_alpha':mo_coeff[0],
+                        'mo_coeff_beta':mo_coeff[1]}
+        self.parm_names=['_alpha','_beta']
+
+        self._mol = mol
+
+    @classmethod
+    def from_pyscf(self, mol, mf, mc=None, tol=-1):
+        """
+        mol: A Mole object
+        mf: An object with mo_coeff and mo_occ. 
+        mc: (optional) a CI object from pyscf
+
+        """
+        obj = mc if hasattr(mc, 'mo_coeff') else mf
+        if mc is not None:
+            detcoeff, occup, det_map = pyqmc.determinant_tools.interpret_ci(mc, tol)
+        else: 
+            detcoeff = np.array([1.0])
+            det_map = np.array([[0],[0]])
+            #occup
+            if len(mf.mo_occ.shape) == 2:
+                occup = [[list(np.argwhere(mf.mo_occ[spin] > 0.5)[:,0])] for spin in [0,1]]
+            else:
+                occup = [[list(np.argwhere(mf.mo_occ > 1.5-spin)[:,0])] for spin in [0,1]]
+
+        max_orb = [int(np.max(occup[s], initial=0)+1) for s in [0,1]]
+        if len(obj.mo_coeff[0].shape) == 2:
+            mo_coeff = [obj.mo_coeff[spin][:,0:max_orb[spin]] for spin in [0,1]]
+        else:
+            mo_coeff = [obj.mo_coeff[:,0:max_orb[spin]] for spin in [0,1]]
+
+        return detcoeff, occup, det_map, MoleculeOrbitalEvaluator(mol, mo_coeff)
+
+    def nmo(self):
+        return [self.parameters['mo_coeff_alpha'].shape[-1], self.parameters['mo_coeff_beta'].shape[-1]]
+
+    def aos(self, eval_str, configs, mask=None):
+        """
+        
+        """
+        mycoords = configs.configs if mask is None else configs.configs[mask]
+        mycoords = mycoords.reshape((-1, mycoords.shape[-1]))
+        aos= np.asarray([self._mol.eval_gto(eval_str, mycoords)])
+        if len(aos.shape)==4: # if derivatives are included
+            return aos.reshape((1,aos.shape[1], *mycoords.shape[:-1], aos.shape[-1]))
+        else:
+            return aos.reshape((1,*mycoords.shape[:-1], aos.shape[-1]))
+
+    def mos(self, ao, spin):
+        return ao[0].dot(self.parameters[f'mo_coeff{self.parm_names[spin]}'])
+
+    def pgradient(self, ao, spin):
+        return np.array([self.parameters[f'mo_coeff{self.parm_names[spin]}'].shape[1]]),ao
+    
+
+    
+def get_k_indices(cell, mf, kpts, tol=1e-6):
+    """Given a list of kpts, return inds such that mf.kpts[inds] is a list of kpts equivalent to the input list"""
+    kdiffs = mf.kpts[np.newaxis] - kpts[:, np.newaxis]
+    frac_kdiffs = np.dot(kdiffs, cell.lattice_vectors().T) / (2 * np.pi)
+    kdiffs = np.mod(frac_kdiffs + 0.5, 1) - 0.5
+    return np.nonzero(np.linalg.norm(kdiffs, axis=-1) < tol)[1]
+
+class PBCOrbitalEvaluatorKpoints:
+    """
+    Evaluate orbitals from a PBC object. 
+    cell is expected to be one made with make_supercell().
+    mo_coeff should be in [spin][k][ao,mo] order
+    kpts should be a list of the k-points corresponding to mo_coeff  
+
+    """
+    def __init__(self, cell, mo_coeff, kpts=None):
+        self.iscomplex=True
+        self._cell = cell.original_cell
+        self.S = cell.S
+
+        self._kpts = [0,0,0] if kpts is None else kpts 
+        self.param_split = [np.cumsum([m.shape[1] for m in mo_coeff[spin]]) for spin in [0,1]]
+        self.parm_names=['_alpha','_beta']
+        self.parameters={'mo_coeff_alpha': np.concatenate(mo_coeff[0], axis=1),
+                         'mo_coeff_beta': np.concatenate(mo_coeff[1], axis=1)}
+
+    @classmethod
+    def from_mean_field(self, cell, mf, twist=None):
+        """
+        mf is expected to be a KUHF, KRHF, or equivalent DFT objects. 
+        Selects occupied orbitals from a given twist 
+        If cell is a supercell, will automatically choose the folded k-points that correspond to that twist.
+        """
+
+        cell = cell if hasattr(cell, 'original_cell') else get_supercell(cell,np.asarray([[1,0,0],[0,1,0],[0,0,1]]))
+
+        if twist is None:
+            twist = np.zeros(3)
+        else:
+            twist = np.dot(np.linalg.inv(cell.a), np.mod(twist, 1.0)) * 2 * np.pi
+        kinds = list(set(get_k_indices(cell, mf, get_supercell_kpts(cell) + twist)))
+        if len(kinds) != cell.scale:
+            raise ValueError("Did not find the right number of k-points for this supercell")
+
+        detcoeff = np.array([1.0])
+        det_map = np.array([[0],[0]])
+
+        if len(mf.mo_coeff[0][0].shape) == 2:
+            occup_k = [[[list(np.argwhere(mf.mo_occ[spin][k] > 0.5)[:,0])] for k in kinds ] for spin in [0,1]]
+        elif len(mf.mo_coeff[0][0].shape) == 1:
+            occup_k = [[[list(np.argwhere(mf.mo_occ[k] > 1.5-spin)[:,0])] for k in kinds ] for spin in [0,1]]
+
+        occup = [[],[]]
+        for spin in [0,1]:
+            count = 0
+            for occ_k in occup_k[spin]:
+                occup[spin] += [o+count for o in occ_k[0]]
+                count+=len(occ_k[0])
+
+        kpts = mf.kpts[kinds]
+        if len(mf.mo_coeff[0][0].shape) == 2:
+            mo_coeff = [[mf.mo_coeff[spin][k][:,occup_k[spin][kinds.index(k)][0]] for k in kinds] for spin in [0,1]]
+        elif len(mf.mo_coeff[0][0].shape) == 1:
+            mo_coeff = [[mf.mo_coeff[k][:, occup_k[spin][kinds.index(k)][0]] for k in kinds] for spin in [0,1]]
+        else:
+            raise ValueError("Did not expect an scf object of type", type(mf))
+        for s in [0,1]:
+            occup[s] = [occup[s]]
+        return detcoeff, occup, det_map, PBCOrbitalEvaluatorKpoints(cell, mo_coeff, kpts)
+        
+    def nmo(self):
+        return [self.parameters['mo_coeff_alpha'].shape[-1], self.parameters['mo_coeff_beta'].shape[-1]]
+
+    def aos(self,eval_str,configs, mask=None):
+        """
+        Returns an ndarray in order [k,..., orbital] of the ao's if value is requested
+
+        if a derivative is requested, will instead return [k,d,...,orbital]
+        """
+        mycoords = configs.configs if mask is None else configs.configs[mask]
+        mycoords = mycoords.reshape((-1, mycoords.shape[-1]))
+
+        #coordinate, dimension
+        wrap = configs.wrap if mask is None else configs.wrap[mask]
+        wrap = np.dot(wrap, self.S)
+        wrap = wrap.reshape((-1,wrap.shape[-1]))
+        kdotR = np.linalg.multi_dot(
+            (self._kpts, self._cell.lattice_vectors().T, wrap.T)
+        )
+        # k, coordinate
+        wrap_phase = get_wrapphase_complex(kdotR)
+        # k,coordinate, orbital
+        ao = np.asarray(self._cell.eval_gto("PBC"+eval_str, mycoords, kpts=self._kpts))
+        ao= np.einsum("k...,k...a->k...a",wrap_phase, ao)
+        if len(ao.shape)==4: # if derivatives are included
+            return ao.reshape((ao.shape[0],ao.shape[1], *mycoords.shape[:-1], ao.shape[-1]))
+        else:
+            return ao.reshape((ao.shape[0],*mycoords.shape[:-1], ao.shape[-1]))
+
+
+        
+    def mos(self, ao, spin):
+        """ao should be [k,[d],...,ao].
+        Returns a concatenated list of all molecular orbitals in form [..., mo]
+
+        In the derivative case, returns [d,..., mo]
+        """
+        p = np.split(self.parameters[f'mo_coeff{self.parm_names[spin]}'], self.param_split[spin], axis=-1)
+        return np.concatenate([ak.dot(mok) for ak,mok in zip(ao,p[0:-1])], axis=-1)
+
+    def pgradient(self,ao, spin):
+        """
+        returns:
+        N sets of atomic orbitals
+        split: which molecular orbitals correspond to which set
+
+        You can construct the determinant by doing, for example:
+        split, aos = pgradient(self.aos)
+        mos = np.split(range(nmo),split)
+        for ao, mo in zip(aos,mos):
+            for i in mo:
+                pgrad[:,:,i] = self._testcol(i,spin,ao)
+                
+        """
+        return self.param_split[spin], ao
diff --git a/pyqmc/pbc.py b/pyqmc/pbc.py
index 44a0a7ff..f0483e9b 100644
--- a/pyqmc/pbc.py
+++ b/pyqmc/pbc.py
@@ -18,10 +18,18 @@ def enforce_pbc(lattvecs, epos):
     # Writes epos in terms of (lattice vecs) fractional coordinates
     recpvecs = np.linalg.inv(lattvecs)
     epos_lvecs_coord = np.dot(epos, recpvecs)
+    to_wrap = np.any((epos_lvecs_coord < 0) | (epos_lvecs_coord > 1), axis=-1)
+    #print(to_wrap.shape)
+    wrapped = np.divmod(epos_lvecs_coord[to_wrap,:],1)
 
+    wraparound=np.zeros(epos.shape)
+    wraparound[to_wrap,:] = wrapped[0]
+
+    final_epos = epos.copy()
+    final_epos[to_wrap,:] = np.dot(wrapped[1],lattvecs)
     # Finds position inside box and wraparound vectors (in lattice vector coordinates)
-    tmp = np.divmod(epos_lvecs_coord, 1)
-    wraparound = tmp[0]
-    final_epos = np.dot(tmp[1], lattvecs)
+    #tmp = np.divmod(epos_lvecs_coord, 1)
+    #wraparound = tmp[0]
+    #final_epos = np.dot(tmp[1], lattvecs)
 
     return final_epos, wraparound
diff --git a/pyqmc/recipes.py b/pyqmc/recipes.py
index c24a500a..d690a36e 100644
--- a/pyqmc/recipes.py
+++ b/pyqmc/recipes.py
@@ -74,7 +74,6 @@ def OPTIMIZE(
 def generate_accumulators(mol, mf, energy=True, rdm1=False, extra_accumulators=None):
     acc = {} if extra_accumulators is None else extra_accumulators
 
-
     if energy:
         if "energy" in acc.keys():
             raise Exception("Found energy in extra_accumulators and energy is True")
diff --git a/pyqmc/slater.py b/pyqmc/slater.py
index 50be2b3c..17eec47c 100644
--- a/pyqmc/slater.py
+++ b/pyqmc/slater.py
@@ -1,331 +1,344 @@
-import numpy as np
-from pyqmc import pbc
-
-
-def sherman_morrison_row(e, inv, vec):
-    tmp = np.einsum("ek,ekj->ej", vec, inv)
-    ratio = tmp[:, e]
-    inv_ratio = inv[:, :, e] / ratio[:, np.newaxis]
-    invnew = inv - np.einsum("ki,kj->kij", inv_ratio, tmp)
-    invnew[:, :, e] = inv_ratio
-    return ratio, invnew
-
-
-_gldict = {"laplacian": np.s_[:1], "gradient_laplacian": np.s_[0:4]}
-
-
-def _aostack_mol(ao, gl):
-    return np.concatenate(
-        [ao[_gldict[gl]], ao[[4, 7, 9]].sum(axis=0, keepdims=True)], axis=0
-    )
-
-
-def _aostack_pbc(ao, gl):
-    return [_aostack_mol(ak, gl) for ak in ao]
-
-
-def get_wrapphase_real(x):
-    return (-1) ** np.round(x / np.pi)
-
-
-def get_wrapphase_complex(x):
-    return np.exp(1j * x)
-
-
-def get_complex_phase(x):
-    return x / np.abs(x)
-
-
-def get_k_indices(cell, mf, kpts, tol=1e-6):
-    """Given a list of kpts, return inds such that mf.kpts[inds] is a list of kpts equivalent to the input list"""
-    kdiffs = mf.kpts[np.newaxis] - kpts[:, np.newaxis]
-    frac_kdiffs = np.dot(kdiffs, cell.lattice_vectors().T) / (2 * np.pi)
-    kdiffs = np.mod(frac_kdiffs + 0.5, 1) - 0.5
-    return np.nonzero(np.linalg.norm(kdiffs, axis=-1) < tol)[1]
-
-
-class PySCFSlater:
-    """A wave function object has a state defined by a reference configuration of electrons.
-    The functions recompute() and updateinternals() change the state of the object, and 
-    the rest compute and return values from that state. """
-
-    def __init__(self, mol, mf, twist=None):
-        """
-        Inputs:
-          supercell: object returned by get_supercell(cell, S)
-          mf: scf object of primitive cell calculation. scf calculation must include k points that fold onto the gamma point of the supercell
-          twist: (3,) array, twisted boundary condition in fractional coordinates, i.e. as coefficients of the reciprocal lattice vectors of the supercell. Integer values are equivalent to zero.
-        """
-        self.parameters = {"det_coeff": np.array([1.0])}
-        self.real_tol = 1e4
-        self._coefflookup = ("mo_coeff_alpha", "mo_coeff_beta")
-
-        if hasattr(mol, "a"):
-            self._init_pbc(mol, mf, twist)
-        else:
-            self._init_mol(mol, mf)
-        self.pbc_str = "PBC" if hasattr(mol, "a") else ""
-        self._aostack = _aostack_pbc if hasattr(mol, "a") else _aostack_mol
-
-        self.dtype = complex if self.iscomplex else float
-        if self.iscomplex:
-            self.get_phase = get_complex_phase
-            self.get_wrapphase = get_wrapphase_complex
-        else:
-            self.get_phase = np.sign
-            self.get_wrapphase = get_wrapphase_real
-
-    def _init_mol(self, mol, mf):
-        from pyscf import scf
-
-        for s, lookup in enumerate(self._coefflookup):
-            if len(mf.mo_occ.shape) == 2:
-                self.parameters[lookup] = mf.mo_coeff[s][
-                    :, np.asarray(mf.mo_occ[s] > 0.9)
-                ]
-            else:
-                minocc = (0.9, 1.1)[s]
-                self.parameters[lookup] = mf.mo_coeff[:, np.asarray(mf.mo_occ > minocc)]
-        self._nelec = tuple(mol.nelec)
-        self._mol = mol
-        self.iscomplex = bool(sum(map(np.iscomplexobj, self.parameters.values())))
-        self.evaluate_orbitals = self._evaluate_orbitals_mol
-        self.evaluate_mos = self._evaluate_mos_mol
-
-    def _init_pbc(self, cell, mf, twist):
-        from pyscf.pbc import scf
-        from pyqmc.supercell import get_supercell_kpts
-
-        # Make sure supercell has attributes S and original_cell
-        for attribute in ["original_cell", "S", "scale"]:
-            if not hasattr(cell, attribute):
-                print('Warning: supercell is missing attribute "%s"' % attribute)
-                print("setting original_cell=supercell and S=np.eye(3)")
-                cell.original_cell = cell
-                cell.S = np.eye(3)
-                cell.scale = 1
-        self.supercell = cell
-        self._cell = cell.original_cell
-
-        # Define kpts
-        if twist is None:
-            twist = np.zeros(3)
-        else:
-            twist = np.dot(np.linalg.inv(cell.a), np.mod(twist, 1.0)) * 2 * np.pi
-        self.kinds = get_k_indices(self._cell, mf, get_supercell_kpts(cell) + twist)
-        self._kpts = mf.kpts[self.kinds]
-        assert len(self.kinds) == len(self._kpts), (self._kpts, mf.kpts)
-        if len(self.kinds)!= cell.scale:
-            raise Exception(f"Expected to find {cell.scale} k-points, but only found {len(self.kinds)}.")
-        self.nk = len(self._kpts)
-
-        # Define parameters
-        self.param_split = {}
-        for s, lookup in enumerate(self._coefflookup):
-            mclist = []
-            for kind in self.kinds:
-                if len(mf.mo_coeff[0][0].shape) == 2:
-                    mca = mf.mo_coeff[s][kind][:, np.asarray(mf.mo_occ[s][kind] > 0.9)]
-                else:
-                    minocc = (0.9, 1.1)[s]
-                    mca = mf.mo_coeff[kind][:, np.asarray(mf.mo_occ[kind] > minocc)]
-                mca = np.real_if_close(mca, tol=self.real_tol)
-                mclist.append(mca / np.sqrt(self.nk))
-            self.param_split[lookup] = np.cumsum([m.shape[1] for m in mclist])
-            self.parameters[lookup] = np.concatenate(mclist, axis=-1)
-
-        self.iscomplex = bool(sum(map(np.iscomplexobj, self.parameters.values())))
-        self.iscomplex = self.iscomplex or np.linalg.norm(self._kpts) > 1e-12
-
-        # Define nelec
-        if len(mf.mo_coeff[0][0].shape) == 2:
-            # Then indices are (spin, kpt, basis, mo)
-            self._nelec = tuple(int(np.sum([o[k] for k in self.kinds])) for o in mf.mo_occ)
-        elif len(mf.mo_coeff[0][0].shape) == 1:
-            # Then indices are (kpt, basis, mo)
-            self._nelec = tuple(
-                int(np.sum([mf.mo_occ[k] > t for k in self.kinds])) for t in (0.9, 1.1)
-            )
-        else:
-            print("Warning: PySCFSlater not expecting scf object of type", type(mf))
-            scale = self.supercell.scale
-            self._nelec = tuple(int(np.round(n * scale)) for n in self._cell.nelec)
-
-        self.evaluate_orbitals = self._evaluate_orbitals_pbc
-        self.evaluate_mos = self._evaluate_mos_pbc
-
-    def _evaluate_orbitals_mol(self, configs, mask=None, eval_str="GTOval_sph"):
-        mycoords = configs.configs if mask is None else configs.configs[mask]
-        mycoords = mycoords.reshape((-1, mycoords.shape[-1]))
-        return self._mol.eval_gto(eval_str, mycoords)
-
-    def _evaluate_mos_mol(self, ao, s):
-        return ao.dot(self.parameters[self._coefflookup[s]])
-
-    def _evaluate_orbitals_pbc(self, configs, mask=None, eval_str="GTOval_sph"):
-        mycoords = configs.configs
-        configswrap = configs.wrap
-        if mask is not None:
-            mycoords = mycoords[mask]
-            configswrap = configswrap[mask]
-        mycoords = mycoords.reshape((-1, mycoords.shape[-1]))
-        # wrap supercell positions into primitive cell
-        prim_coords, prim_wrap = pbc.enforce_pbc(self._cell.lattice_vectors(), mycoords)
-        configswrap = configswrap.reshape(prim_wrap.shape)
-        wrap = prim_wrap + np.dot(configswrap, self.supercell.S)
-        kdotR = np.linalg.multi_dot(
-            (self._kpts, self._cell.lattice_vectors().T, wrap.T)
-        )
-        wrap_phase = self.get_wrapphase(kdotR)
-        # evaluate AOs for all electron positions
-        ao = self._cell.eval_gto("PBC" + eval_str, prim_coords, kpts=self._kpts)
-        ao = [ao[k] * wrap_phase[k][:, np.newaxis] for k in range(self.nk)]
-        return ao
-
-    def _evaluate_mos_pbc(self, aos, s):
-        """
-        Evaluate MOs for spin s given aos
-        """
-        c = self._coefflookup[s]
-        p = np.split(self.parameters[c], self.param_split[c], axis=-1)
-        mo = [ao.dot(p[k]) for k, ao in enumerate(aos)]
-        return np.concatenate(mo, axis=-1)
-
-    def recompute(self, configs):
-        """This computes the value from scratch. Returns the logarithm of the wave function as
-        (phase,logdet). If the wf is real, phase will be +/- 1."""
-        nconf, nelec, ndim = configs.configs.shape
-        aos = self.evaluate_orbitals(configs)
-        if hasattr(self, "nk"):
-            aos_shape = (self.nk, nconf, nelec, -1)
-        else:
-            aos_shape = (1, nconf, nelec, -1)
-        aos = np.reshape(aos, aos_shape)
-        self._aovals = aos
-        self._dets = []
-        self._inverse = []
-        for s in [0, 1]:
-            i0, i1 = s * self._nelec[0], self._nelec[0] + s * self._nelec[1]
-            ne = self._nelec[s]
-            mo = self.evaluate_mos(aos[:, :, i0:i1], s).reshape(nconf, ne, ne)
-            phase, mag = np.linalg.slogdet(mo)
-            self._dets.append((phase, mag))
-            self._inverse.append(np.linalg.inv(mo))
-
-        return self.value()
-
-    def updateinternals(self, e, epos, mask=None):
-        s = int(e >= self._nelec[0])
-        if mask is None:
-            mask = [True] * epos.configs.shape[0]
-        eeff = e - s * self._nelec[0]
-        aos = self.evaluate_orbitals(epos, mask=mask)
-        self._aovals[:, mask, e, :] = np.asarray(aos)  # (kpt, config, ao)
-        mo = self.evaluate_mos(aos, s)
-        ratio, self._inverse[s][mask, :, :] = sherman_morrison_row(
-            eeff, self._inverse[s][mask, :, :], mo
-        )
-        self._updateval(ratio, s, mask)
-
-    def _updateval(self, ratio, s, mask):
-        self._dets[s][0][mask] *= self.get_phase(ratio)
-        self._dets[s][1][mask] += np.log(np.abs(ratio))
-
-    ### not state-changing functions
-
-    def value(self):
-        """Return logarithm of the wave function as noted in recompute()"""
-        return (
-            self._dets[0][0] * self._dets[1][0],
-            self._dets[0][1]
-            + self._dets[1][1]
-            + np.log(np.abs(self.parameters["det_coeff"][0])),
-        )
-
-    def _testrow(self, e, vec, mask=None, spin=None):
-        """vec is a nconfig,nmo vector which replaces row e"""
-        s = int(e >= self._nelec[0]) if spin is None else spin
-        elec = e - s * self._nelec[0]
-        if mask is None:
-            return np.einsum("i...j,ij...->i...", vec, self._inverse[s][:, :, elec])
-
-        return np.einsum("i...j,ij...->i...", vec, self._inverse[s][mask, :, elec])
-
-    def _testcol(self, i, s, vec):
-        """vec is a nconfig,nmo vector which replaces column i"""
-        return np.einsum("ij...,ij->i...", vec, self._inverse[s][:, i, :])
-
-    def testvalue(self, e, epos, mask=None):
-        """ return the ratio between the current wave function and the wave function if 
-        electron e's position is replaced by epos"""
-        s = int(e >= self._nelec[0])
-        nmask = epos.configs.shape[0] if mask is None else np.sum(mask)
-        if nmask == 0:
-            return np.zeros((0, epos.configs.shape[1]))
-        aos = self.evaluate_orbitals(epos, mask)
-        mo = self.evaluate_mos(aos, s)
-        mo = mo.reshape(nmask, *epos.configs.shape[1:-1], self._nelec[s])
-        return self._testrow(e, mo, mask)
-
-    def testvalue_many(self, e, epos, mask=None):
-        """ return the ratio between the current wave function and the wave function if 
-        an electron's position is replaced by epos for each electron"""
-        s = (e >= self._nelec[0]).astype(int)
-        nmask = epos.configs.shape[0] if mask is None else np.sum(mask)
-        if nmask == 0:
-            return np.zeros((0, epos.configs.shape[1]))
-
-        aos = self.evaluate_orbitals(epos, mask)
-        ratios = np.zeros((epos.configs.shape[0], e.shape[0]), dtype=self.dtype)
-        for spin in [0, 1]:
-            ind = s == spin
-            mo = self.evaluate_mos(aos, spin)
-            mo = mo.reshape(nmask, *epos.configs.shape[1:-1], self._nelec[spin])
-            ratios[:, ind] = self._testrow(e[ind], mo, mask=mask, spin=spin)
-        return ratios
-
-    def gradient(self, e, epos):
-        """ Compute the gradient of the log wave function 
-        Note that this can be called even if the internals have not been updated for electron e,
-        if epos differs from the current position of electron e."""
-        s = int(e >= self._nelec[0])
-        aograd = self.evaluate_orbitals(epos, eval_str="GTOval_sph_deriv1")
-        mograd = self.evaluate_mos(aograd, s)
-        ratios = np.asarray([self._testrow(e, x) for x in mograd])
-        return ratios[1:] / ratios[:1]
-
-    def laplacian(self, e, epos):
-        s = int(e >= self._nelec[0])
-        ao = self.evaluate_orbitals(epos, eval_str="GTOval_sph_deriv2")
-        mo = self.evaluate_mos(self._aostack(ao, "laplacian"), s)
-        ratios = np.asarray([self._testrow(e, x) for x in mo])
-        return ratios[1] / ratios[0]
-
-    def gradient_laplacian(self, e, epos):
-        s = int(e >= self._nelec[0])
-        ao = self.evaluate_orbitals(epos, eval_str="GTOval_sph_deriv2")
-        mo = self.evaluate_mos(self._aostack(ao, "gradient_laplacian"), s)
-        ratios = np.asarray([self._testrow(e, x) for x in mo])
-        return ratios[1:-1] / ratios[:1], ratios[-1] / ratios[0]
-
-    def pgradient(self):
-        d = {"det_coeff": np.zeros(self._aovals.shape[-3])}
-        for parm in ["mo_coeff_alpha", "mo_coeff_beta"]:
-            s = int("beta" in parm)
-            # Get AOs for our spin channel only
-            i0, i1 = s * self._nelec[0], self._nelec[0] + s * self._nelec[1]
-            ao = self._aovals[:, :, i0:i1, :]  # (kpt, config, electron, ao)
-            pgrad_shape = (ao.shape[-3],) + self.parameters[parm].shape
-            pgrad = np.zeros(pgrad_shape, dtype=self.dtype)  # (nconf, coeff)
-            # Compute derivatives w.r.t. MO coefficients
-            if ao.shape[0] > 1:  # multiple kpts
-                split_sizes = np.diff([0] + list(self.param_split[parm]))
-                k = np.repeat(np.arange(self.nk), split_sizes)
-                for i in range(self._nelec[s]):  # MO loop
-                    pgrad[:, :, i] = self._testcol(i, s, ao[k[i]])
-            else:
-                ao = ao[0]
-                for i in range(self._nelec[s]):  # MO loop
-                    pgrad[:, :, i] = self._testcol(i, s, ao)
-            d[parm] = np.asarray(pgrad)
-        return d
+import numpy as np
+#from pyqmc.slater import sherman_morrison_row, get_complex_phase
+import pyqmc.determinant_tools as determinant_tools
+import pyqmc.orbitals
+
+
+def sherman_morrison_row(e, inv, vec):
+    tmp = np.einsum("ek,ekj->ej", vec, inv)
+    ratio = tmp[:, e]
+    inv_ratio = inv[:, :, e] / ratio[:, np.newaxis]
+    invnew = inv - np.einsum("ki,kj->kij", inv_ratio, tmp)
+    invnew[:, :, e] = inv_ratio
+    return ratio, invnew
+
+
+def get_complex_phase(x):
+    return x / np.abs(x)
+
+class JoinParameters:
+    """
+    This class provides a dict-like interface that actually references 
+    other dictionaries in the background.
+    If keys collide, then the first dictionary that matches the key will be returned.
+    However, some bad things may happen if you have colliding keys.
+    """
+    def __init__(self, dicts):
+        self.data = {}
+        self.data = dicts
+
+
+    def find_i(self,idx):
+        for i,d in enumerate(self.data):
+            if idx in d:
+                return i
+
+
+    def __setitem__(self, idx, value):
+        i = self.find_i(idx)
+        self.data[i][idx] = value
+
+    def __getitem__(self, idx):
+        i = self.find_i(idx)
+        return self.data[i][idx]
+
+    def __delitem__(self, idx):
+        i = self.find_i(idx)
+        del self.data[i][idx]
+
+    def __iter__(self):
+        for d in self.data:
+            yield from d.keys()
+
+    def __len__(self):
+        return sum(len(i) for i in self.data)
+
+    def items(self):
+        for d in self.data:
+            yield from d.items()
+
+    def __repr__(self):
+        return self.data.__repr__()
+
+    def keys(self):
+        for d in self.data:
+            yield from d.keys()
+    
+    def values(self):
+        for d in self.data:
+            yield from d.values()
+
+
+
+def sherman_morrison_ms(e, inv, vec):
+    ratio = np.einsum("idj,idj->id", vec, inv[:, :, :, e])
+    tmp = np.einsum("edk,edkj->edj", vec, inv)
+    invnew = (
+        inv
+        - np.einsum("kdi,kdj->kdij", inv[:, :, :, e], tmp)
+        / ratio[:, :, np.newaxis, np.newaxis]
+    )
+    invnew[:, :, :, e] = inv[:, :, :, e] / ratio[:, :, np.newaxis]
+    return ratio, invnew
+
+
+
+
+class Slater:
+    """
+    A multi-determinant wave function object initialized
+    via an SCF calculation.
+
+    How to use with hci
+
+    .. code-block:: python
+
+        cisolver = pyscf.hci.SCI(mol)
+        cisolver.select_cutoff=0.1
+        nmo = mf.mo_coeff.shape[1]
+        nelec = mol.nelec
+        h1 = mf.mo_coeff.T.dot(mf.get_hcore()).dot(mf.mo_coeff)
+        h2 = pyscf.ao2mo.full(mol, mf.mo_coeff)
+        e, civec = cisolver.kernel(h1, h2, nmo, nelec, verbose=4)
+        cisolver.ci = civec[0]
+        wf = pyqmc.multislater.MultiSlater(mol, mf, cisolver, tol=0.1)
+
+
+    """
+
+    def __init__(self, mol, mf, mc=None, tol=None, twist=None):
+        self.tol = -1 if tol is None else tol
+        self._mol = mol
+        if hasattr(mc, "nelecas"):
+            # In case nelecas overrode the information from the molecule object.
+            self._nelec = (mc.nelecas[0] + mc.ncore, mc.nelecas[1] + mc.ncore)
+        else:
+            self._nelec = mol.nelec
+
+        self.myparameters={}
+        self.myparameters["det_coeff"], self._det_occup, self._det_map,\
+        self.orbitals = pyqmc.orbitals.choose_evaluator_from_pyscf(mol, mf, mc, twist=twist)
+        self.parameters=JoinParameters([self.myparameters,self.orbitals.parameters])
+
+        self.iscomplex = bool(sum(map(np.iscomplexobj, self.parameters.values())))
+        self.dtype = complex if self.iscomplex else float
+        self.get_phase = get_complex_phase if self.iscomplex else np.sign
+
+
+    def recompute(self, configs):
+        """This computes the value from scratch. Returns the logarithm of the wave function as
+        (phase,logdet). If the wf is real, phase will be +/- 1."""
+
+        nconf, nelec, ndim = configs.configs.shape
+
+        aos = self.orbitals.aos('GTOval_sph', configs)
+        self._aovals = aos.reshape(-1,nconf,nelec, aos.shape[-1])
+        self._dets = []
+        self._inverse = []
+        for s in [0, 1]:
+            begin = self._nelec[0]*s
+            end = self._nelec[0] + self._nelec[1] * s
+            mo = self.orbitals.mos(self._aovals[:,:,begin:end, :],s)
+            mo_vals = np.swapaxes(mo[:, :, self._det_occup[s]], 1, 2)
+            self._dets.append(
+                np.array(np.linalg.slogdet(mo_vals))
+            )  # Spin, (sign, val), nconf, [ndet_up, ndet_dn]
+            self._inverse.append(
+                np.linalg.inv(mo_vals)
+            )  # spin, Nconf, [ndet_up, ndet_dn], nelec, nelec
+        return self.value()
+
+    def updateinternals(self, e, epos, mask=None):
+        """Update any internals given that electron e moved to epos. mask is a Boolean array 
+        which allows us to update only certain walkers"""
+
+        s = int(e >= self._nelec[0])
+        if mask is None:
+            mask = [True] * epos.configs.shape[0]
+        eeff = e - s * self._nelec[0]
+        ao = self.orbitals.aos("GTOval_sph",epos,mask)
+        self._aovals[:,mask, e, :] = ao
+        mo = self.orbitals.mos(ao,s)
+
+        mo_vals = mo[:, self._det_occup[s]]
+        det_ratio, self._inverse[s][mask, :, :, :] = sherman_morrison_ms(
+            eeff, self._inverse[s][mask, :, :, :], mo_vals
+        )
+
+        self._updateval(det_ratio, s, mask)
+
+    def value(self):
+        """Return logarithm of the wave function as noted in recompute()"""
+        updets = self._dets[0][:, :, self._det_map[0]]
+        dndets = self._dets[1][:, :, self._det_map[1]]
+        return determinant_tools.compute_value(updets,dndets, self.parameters["det_coeff"])
+
+    def _updateval(self, ratio, s, mask):
+        self._dets[s][0, mask, :] *= self.get_phase(ratio)
+        self._dets[s][1, mask, :] += np.log(np.abs(ratio))
+
+    def _testrow(self, e, vec, mask=None, spin=None):
+        """vec is a nconfig,nmo vector which replaces row e"""
+        s = int(e >= self._nelec[0]) if spin is None else spin
+        if mask is None:
+            mask = [True] * vec.shape[0]
+
+        ratios = np.einsum(
+            "i...dj,idj...->i...d",
+            vec,
+            self._inverse[s][mask][..., e - s * self._nelec[0]],
+        )
+
+        upref = np.amax(self._dets[0][1]).real
+        dnref = np.amax(self._dets[1][1]).real
+        
+        det_array = (
+            self._dets[0][0, :, self._det_map[0]][:, mask]
+            * self._dets[1][0, :, self._det_map[1]][:, mask]
+            * np.exp(
+                self._dets[0][1, :, self._det_map[0]][:, mask]
+                + self._dets[1][1, :, self._det_map[1]][:, mask]-upref-dnref
+            )
+        )
+        numer = np.einsum(
+            "i...d,d,di->i...",
+            ratios[..., self._det_map[s]],
+            self.parameters["det_coeff"],
+            det_array,
+        )
+        denom = np.einsum(
+            "d,di->i...",
+            self.parameters["det_coeff"],
+            det_array,
+        )
+        #curr_val = self.value()
+        
+        if len(numer.shape) == 2:
+            denom = denom[:, np.newaxis]
+        return numer / denom
+
+    def _testcol(self, det, i, s, vec):
+        """vec is a nconfig,nmo vector which replaces column i 
+        of spin s in determinant det"""
+
+        return np.einsum(
+            "ij...,ij->i...", vec, self._inverse[s][:, det, i, :], optimize="greedy"
+        )
+
+    def gradient(self, e, epos):
+        """ Compute the gradient of the log wave function 
+        Note that this can be called even if the internals have not been updated for electron e,
+        if epos differs from the current position of electron e."""
+        s = int(e >= self._nelec[0])
+        aograd = self.orbitals.aos('GTOval_sph_deriv1', epos)
+        mograd = self.orbitals.mos(aograd, s)
+
+        mograd_vals = mograd[:, :, self._det_occup[s]]
+
+        ratios = np.asarray([self._testrow(e, x) for x in mograd_vals])
+        return ratios[1:] / ratios[0]
+
+    def laplacian(self, e, epos):
+        """ Compute the laplacian Psi/ Psi. """
+        s = int(e >= self._nelec[0])
+        ao = self.orbitals.aos("GTOval_sph_deriv2", epos)
+        ao_val = ao[:,0,:,:]
+        ao_lap = np.sum(ao[:,[4,7,9],:,:],axis=1)
+        mos = [self.orbitals.mos(x,s)[...,self._det_occup[s]] for x in [ao_val, ao_lap]]
+        ratios = [self._testrow(e,mo) for mo in mos]
+        return ratios[1]/ratios[0]
+
+    def gradient_laplacian(self, e, epos):
+        s = int(e >= self._nelec[0])
+        ao = self.orbitals.aos("GTOval_sph_deriv2", epos)
+        ao = np.concatenate([ao[:,0:4,...], ao[:,[4,7,9],...].sum(axis=1,keepdims=True)],axis=1)
+        mo = self.orbitals.mos(ao,s)
+        mo_vals = mo[:, :, self._det_occup[s]]
+        ratios = np.asarray([self._testrow(e, x) for x in mo_vals])
+        return ratios[1:-1] / ratios[:1], ratios[-1] / ratios[0]
+
+    def testvalue(self, e, epos, mask=None):
+        """ return the ratio between the current wave function and the wave function if 
+        electron e's position is replaced by epos"""
+        s = int(e >= self._nelec[0])
+        ao = self.orbitals.aos('GTOval_sph',epos, mask)
+        mo = self.orbitals.mos(ao, s)
+        mo_vals = mo[..., self._det_occup[s]]
+        if len(epos.configs.shape) > 2:
+            mo_vals = mo_vals.reshape(-1, epos.configs.shape[1], mo_vals.shape[1], mo_vals.shape[2])
+        return self._testrow(e, mo_vals, mask)
+
+    def testvalue_many(self, e, epos, mask=None):
+        """ return the ratio between the current wave function and the wave function if 
+        electron e's position is replaced by epos for each electron"""
+        s = (e >= self._nelec[0]).astype(int)
+        ao = self.orbitals.aos('GTOval_sph', epos, mask)
+        ratios = np.zeros((epos.configs.shape[0], e.shape[0]),dtype=self.dtype)
+        for spin in [0, 1]:
+            ind = s == spin
+            mo = self.orbitals.mos(ao, spin)
+            mo = mo.reshape(-1, *epos.configs.shape[1:-1], self._nelec[spin])
+            mo_vals = mo[..., self._det_occup[spin]]
+            ratios[:, ind] = self._testrow(e[ind], mo_vals, mask, spin=spin)
+
+        return ratios
+
+    def pgradient(self):
+        r"""Compute the parameter gradient of Psi. 
+        Returns $$d_p \Psi/\Psi$$ as a dictionary of numpy arrays,
+        which correspond to the parameter dictionary.
+        
+        The wave function is given by ci Di, with an implicit sum
+
+        We have two sets of parameters:
+
+        Determinant coefficients: 
+        di psi/psi = Dui Ddi/psi
+
+        Orbital coefficients:
+        dj psi/psi = ci dj (Dui Ddi)/psi
+
+        Let's suppose that j corresponds to an up orbital coefficient. Then 
+        dj (Dui Ddi) = (dj Dui)/Dui Dui Ddi/psi = (dj Dui)/Dui di psi/psi
+        where di psi/psi is the derivative defined above.
+        """
+        d = {}
+
+        # Det coeff
+        curr_val = self.value()
+        d["det_coeff"] = (
+            self._dets[0][0, :, self._det_map[0]]
+            * self._dets[1][0, :, self._det_map[1]]
+            * np.exp(
+                self._dets[0][1, :, self._det_map[0]]
+                + self._dets[1][1, :, self._det_map[1]]
+                - curr_val[1]
+            )/curr_val[0]
+        ).T
+
+        for s,parm in zip([0,1],["mo_coeff_alpha", "mo_coeff_beta"]):
+            ao = self._aovals[
+                :, :, s * self._nelec[0] : self._nelec[s] + s * self._nelec[0], :
+            ]
+
+            split, aos = self.orbitals.pgradient(ao, s)
+            mos = np.split(range(split[-1]), split)
+            # Compute dj Diu/Diu
+            nao = aos[0].shape[-1]
+            nconf = aos[0].shape[0]
+            nmo = split[-1]
+            deriv = np.zeros((len(self._det_occup[s]), nconf, nao, nmo),dtype=curr_val[0].dtype)
+            for det, occ in enumerate(self._det_occup[s]):
+                for ao, mo in zip(aos, mos):
+                    for i in mo:
+                        if i in occ:
+                            col = occ.index(i)
+                            deriv[det, :, :, i]= self._testcol(det, col, s, ao) 
+
+            # now we reduce over determinants
+            d[parm] = np.zeros(deriv.shape[1:], dtype=curr_val[0].dtype)
+            for di,coeff in enumerate(self.parameters['det_coeff']):
+                whichdet = self._det_map[s][di]
+                d[parm] += deriv[whichdet]*coeff*d["det_coeff"][:,di, np.newaxis, np.newaxis]
+
+        return d
diff --git a/pyqmc/tbdm.py b/pyqmc/tbdm.py
index cf9f4ed4..e46be544 100644
--- a/pyqmc/tbdm.py
+++ b/pyqmc/tbdm.py
@@ -4,12 +4,18 @@
 from pyqmc.mc import initial_guess
 from pyqmc.obdm import sample_onebody
 from sys import stdout
-
+from pyqmc.orbitals import MoleculeOrbitalEvaluator, PBCOrbitalEvaluatorKpoints
+import pyqmc.supercell as supercell
 
 class TBDMAccumulator:
     """ Returns one spin sector of the tbdm[s1,s2] as an array (norb_s1,norb_s1,norb_s2,norb_s2) with indices (using pySCF's 
     convention): tbdm[s1,s2][i,j,k,l] = < c^+_{s1,i} c^+_{s2,k} c_{s2,l} c_{s1,j} > = \phi*_{s1,j} \phi*_{s2,l} \phi_{s2,k} \phi_{s1,i}.
 
+
+    We use pySCF's index convention (while Eq. 10 in DOI:10.1063/1.4793531 uses QWalk's)
+    QWalk -> tbdm[s1,s2,i,j,k,l] = < c^+_{s1,i} c^+_{s2,j} c_{s2,l} c_{s1,k} > = \phi*_{s1,k} \phi*_{s2,l} \phi_{s2,j} \phi_{s1,i}
+    pySCF -> tbdm[s1,s2,i,j,k,l] = < c^+_{s1,i} c^+_{s2,k} c_{s2,l} c_{s1,j} > = \phi*_{s1,j} \phi*_{s2,l} \phi_{s2,k} \phi_{s1,i}
+
     Args:
 
       mol (Mole): PySCF Mole object.
@@ -42,47 +48,38 @@ def __init__(
         warmup=200,
         naux=500,
         ijkl=None,
+        kpts=None
     ):
         assert (
             len(orb_coeff.shape) == 3
         ), "orb_coeff should be a list of orbital coefficients with size (2,num_mobasis,num_orb)."
 
         self._mol = mol
-        self._orb_coeff = orb_coeff
         self._tstep = tstep
         self._nsweeps = nsweeps
         self._spin = spin
 
+        if kpts is None:
+            self.orbitals = MoleculeOrbitalEvaluator(mol, orb_coeff)
+        else: 
+            if not hasattr(mol, "original_cell"):
+                mol = supercell.get_supercell(mol, np.eye(3))
+            self.orbitals = PBCOrbitalEvaluatorKpoints(mol, orb_coeff, kpts)
+
         self._spin_sector = spin
-        self._electrons_a = np.arange(
-            spin[0] * mol.nelec[0], mol.nelec[0] + spin[0] * mol.nelec[1]
-        )
-        self._electrons_b = np.arange(
-            spin[1] * mol.nelec[0], mol.nelec[0] + spin[1] * mol.nelec[1]
-        )
-        self._pairs = np.array(
-            np.meshgrid(self._electrons_a, self._electrons_b)
-        ).T.reshape(-1, 2)
-        self._pairs = self._pairs[
-            self._pairs[:, 0] != self._pairs[:, 1]
-        ]  # Removes repeated electron pairs
-
-        # Initialization and warmup of aux_configs_a
-        self._aux_configs_a = initial_guess(
-            mol, int(naux / sum(self._mol.nelec))
-        ).configs.reshape(-1, 3)
-        for i in range(warmup):
-            accept_a, self._aux_configs_a = sample_onebody(
-                mol, orb_coeff[self._spin_sector[0]], self._aux_configs_a, tstep
-            )
-        # Initialization and warmup of aux_configs_b
-        self._aux_configs_b = initial_guess(
-            mol, int(naux / sum(self._mol.nelec))
-        ).configs.reshape(-1, 3)
-        for i in range(warmup):
-            accept_b, self._aux_configs_b = sample_onebody(
-                mol, orb_coeff[self._spin_sector[1]], self._aux_configs_b, tstep
-            )
+        self._electrons=[np.arange(spin[s] * mol.nelec[0], mol.nelec[0] + spin[s] * mol.nelec[1]) for s in [0,1]]
+
+        # Initialization and warmup of configurations
+        nwalkers = int(naux / sum(self._mol.nelec))
+        self._aux_configs=[]
+        for spin in [0,1]:
+            self._aux_configs.append(initial_guess(mol, nwalkers))
+            self._aux_configs[spin].reshape((-1,1, 3))
+            _, self._aux_configs[spin], _ = \
+                sample_onebody(self._aux_configs[spin], self.orbitals, 0, nsamples=warmup)
+            self._aux_configs[spin] = self._aux_configs[spin][-1]
+
+
         # Default to full 2rdm if ijkl not specified
         if ijkl is None:
             norb_up = orb_coeff[0].shape[1]
@@ -96,129 +93,103 @@ def __init__(
             ]
         self._ijkl = np.array(ijkl).T
 
-    def __call__(self, configs, wf, extra_configs=None, auxassignments=None):
-        """Gathers quantities from equation (10) of DOI:10.1063/1.4793531."""
+    def get_configurations(self, nconf):
+        """
+        Obtain a sequence of auxilliary configurations. This function returns one auxilliary configuration 
+        for each nconf.
+        Changes internal state: self._aux_configs is updated to the last sampled location. 
+
+        This will resample the auxilliary configurations to match the number of walkers.
+
+        returns a dictionary with the following elements, separated by spin:
+            assignments: [nsweeps, nconf]: assignment of configurations for each sweep to an auxilliary walker.
+            orbs: [nsweeps, conf, norb]: orbital values
+            configs: [nsweeps] Configuration object with nconf configurations of 1 electron
+            acceptance: [nsweeps, naux] acceptance probability for each auxilliary walker
+
+        TODO: Should we just resize the configurations to nconf instead of taking naux as an input?
+        """
+        configs=[]
+        assignments = []
+        orbs = []
+        acceptance = []
+
+        for spin in [0,1]:
+            naux = self._aux_configs[spin].configs.shape[0]
+            accept, tmp_config, tmp_orbs = sample_onebody(self._aux_configs[spin], self.orbitals, spin, self._nsweeps, tstep=self._tstep)
+            assignments.append(np.random.randint(0, naux, size=(self._nsweeps, nconf)))
+            self._aux_configs[spin] = tmp_config[-1].copy()
+            acceptance.append(accept)
+            for conf,assign in zip(tmp_config,assignments[-1]):
+                conf.resample(assign)
+            configs.append(tmp_config)
+            orbs.append([orb[assign,...] for orb, assign in zip(tmp_orbs,assignments[-1])])
+
+        return {'acceptance':acceptance, 'orbs':orbs, 'configs':configs, 'assignments':assignments}
+
+
+    def __call__(self, configs, wf):
+        """Gathers quantities from equation (10) of DOI:10.1063/1.4793531.
+        
+        assignments maps from the auxilliary walkers onto the main walkers. 
+        It should be of length [nsweeps,nconf], and contain integers between 0 and naux.
+        """
 
-        # Constructs results dictionary
         nconf = configs.configs.shape[0]
-        results = {}
-        orb_a_size = self._orb_coeff[self._spin_sector[0]].shape[1]
-        orb_b_size = self._orb_coeff[self._spin_sector[1]].shape[1]
-        results["value"] = np.zeros((nconf, self._ijkl.shape[1]))
-        for i, e in enumerate(["a", "b"]):
-            results["norm_%s" % e] = np.zeros(
-                (nconf, self._orb_coeff[self._spin_sector[i]].shape[1])
-            )
-            results["acceptance_%s" % e] = np.zeros(nconf)
-
-        # Returns empty arrays if no electron pairs
-        if len(self._pairs) == 0:
-            return results
-
-        if extra_configs is None:
-            # Generates aux_configs_a and aux_configs_b
-            aux_configs_a = []
-            aux_configs_b = []
-            for step in range(self._nsweeps):
-                aux_configs_a.append(np.copy(self._aux_configs_a))
-                accept_a, self._aux_configs_a = sample_onebody(
-                    self._mol,
-                    self._orb_coeff[self._spin_sector[0]],
-                    self._aux_configs_a,
-                    tstep=self._tstep,
-                )
-                aux_configs_b.append(np.copy(self._aux_configs_b))
-                accept_b, self._aux_configs_b = sample_onebody(
-                    self._mol,
-                    self._orb_coeff[self._spin_sector[1]],
-                    self._aux_configs_b,
-                    tstep=self._tstep,
-                )
-                results["acceptance_a"] += np.mean(accept_a)
-                results["acceptance_b"] += np.mean(accept_b)
-            results["acceptance_a"] /= self._nsweeps
-            results["acceptance_b"] /= self._nsweeps
-            aux_configs_a = np.array(aux_configs_a)
-            aux_configs_b = np.array(aux_configs_b)
-            # Generates random choice of aux_config_a and aux_config_b for moving electron_a and electron_b
-            naux_a = self._aux_configs_a.shape[0]
-            naux_b = self._aux_configs_b.shape[0]
-            auxassignments_a = np.random.randint(0, naux_a, size=(self._nsweeps, nconf))
-            auxassignments_b = np.random.randint(0, naux_b, size=(self._nsweeps, nconf))
-        else:
-            assert auxassignments is not None
-            aux_configs_a = extra_configs[0]
-            aux_configs_b = extra_configs[1]
-            naux_a = self._aux_configs_a.shape[0]
-            naux_b = self._aux_configs_b.shape[0]
-            auxassignments_a = auxassignments[0]
-            auxassignments_b = auxassignments[1]
-
-        # Evaluate VMC configurations
-        coords = configs.configs.reshape(
-            (configs.configs.shape[0] * configs.configs.shape[1], -1)
-        )
-        ao_configs = self._mol.eval_gto("GTOval_sph", coords)
-        orb_a_configs = ao_configs.dot(self._orb_coeff[self._spin_sector[0]]).reshape(
-            (configs.configs.shape[0], configs.configs.shape[1], -1)
-        )
-        orb_b_configs = ao_configs.dot(self._orb_coeff[self._spin_sector[1]]).reshape(
-            (configs.configs.shape[0], configs.configs.shape[1], -1)
-        )
-        orb_a_configs = orb_a_configs[:, self._pairs[:, 0], :]
-        orb_b_configs = orb_b_configs[:, self._pairs[:, 1], :]
-
-        # Sweeps over electron pairs
+        aux = self.get_configurations(nconf)
+
+        # Evaluate orbital values for the primary samples
+        ao_configs = self.orbitals.aos("GTOval_sph", configs)
+        ao_configs = ao_configs.reshape((ao_configs.shape[0], nconf, -1, ao_configs.shape[-1]))
+        orb_configs = [self.orbitals.mos(ao_configs[...,self._electrons[spin], :], spin) for spin in [0,1]]
+        results = {"value": np.zeros((nconf, self._ijkl.shape[1])),
+            'norm_a': np.zeros((nconf, orb_configs[0].shape[-1])),
+            'norm_b': np.zeros((nconf, orb_configs[1].shape[-1])),
+            'acceptance_a': np.mean(aux['acceptance'][0], axis=0),
+            'acceptance_b': np.mean(aux['acceptance'][0],axis=0) }
+        orb_configs = [orb_configs[s][:,:, self._ijkl[2*s]] for s in [0,1]]
+
+        down_start = [np.min(self._electrons[s]) for s in [0,1]]
         for sweep in range(self._nsweeps):
-            ao_a_aux = self._mol.eval_gto("GTOval_sph", aux_configs_a[sweep])
-            ao_b_aux = self._mol.eval_gto("GTOval_sph", aux_configs_b[sweep])
-            orb_a_aux = ao_a_aux.dot(self._orb_coeff[self._spin_sector[0]])
-            orb_b_aux = ao_b_aux.dot(self._orb_coeff[self._spin_sector[1]])
-            fsum_a = np.sum(orb_a_aux * orb_a_aux, axis=1)
-            fsum_b = np.sum(orb_b_aux * orb_b_aux, axis=1)
-            norm_a = orb_a_aux * orb_a_aux / fsum_a[:, np.newaxis]
-            norm_b = orb_b_aux * orb_b_aux / fsum_b[:, np.newaxis]
-
-            # We use pySCF's index convention (while Eq. 10 in DOI:10.1063/1.4793531 uses QWalk's)
-            # QWalk -> tbdm[s1,s2,i,j,k,l] = < c^+_{s1,i} c^+_{s2,j} c_{s2,l} c_{s1,k} > = \phi*_{s1,k} \phi*_{s2,l} \phi_{s2,j} \phi_{s1,i}
-            # pySCF -> tbdm[s1,s2,i,j,k,l] = < c^+_{s1,i} c^+_{s2,k} c_{s2,l} c_{s1,j} > = \phi*_{s1,j} \phi*_{s2,l} \phi_{s2,k} \phi_{s1,i}
-            orbratio = (
-                (
-                    orb_a_aux[auxassignments_a[sweep]][:, self._ijkl[1]]
-                    / fsum_a[auxassignments_a[sweep], np.newaxis]
-                )[:, np.newaxis, :]
-                * (
-                    orb_b_aux[auxassignments_b[sweep]][:, self._ijkl[3]]
-                    / fsum_b[auxassignments_b[sweep], np.newaxis]
-                )[:, np.newaxis, :]
-                * orb_a_configs[..., self._ijkl[0]]
-                * orb_b_configs[..., self._ijkl[2]]
-            )
-
-            # Calculation of wf ratio (no McMillan trick yet)
-            epos_a = configs.make_irreducible(
-                -1, aux_configs_a[sweep][auxassignments_a[sweep]]
-            )
-            epos_b = configs.make_irreducible(
-                -1, aux_configs_b[sweep][auxassignments_b[sweep]]
-            )
+            fsum = [np.sum(np.abs(aux['orbs'][spin][sweep])**2, axis=1) for spin in [0,1]]
+            norm = [np.abs(aux['orbs'][spin][sweep])**2/fsum[spin][:,np.newaxis] for spin in [0,1]]
 
             wfratio = []
-            for ea in self._electrons_a:
-                electrons_b = self._electrons_b[self._electrons_b != ea]
+            electrons_a_ind = []
+            electrons_b_ind = []
+            for ea in self._electrons[0]:
+                # Don't move the same electron twice
+                electrons_b = self._electrons[1][self._electrons[1] != ea]
+                epos_a = aux['configs'][0][sweep].electron(0)
+                epos_b = aux['configs'][1][sweep].electron(0)
                 wfratio_a = wf.testvalue(ea, epos_a)
                 wf.updateinternals(ea, epos_a)
                 wfratio_b = wf.testvalue_many(electrons_b, epos_b)
                 wf.updateinternals(ea, configs.electron(ea))
                 wfratio.append(wfratio_a[:, np.newaxis] * wfratio_b)
+                electrons_a_ind.extend([ea-down_start[0]]*len(electrons_b))
+                electrons_b_ind.extend(electrons_b-down_start[1])
             wfratio = np.concatenate(wfratio, axis=1)
 
-            # Adding to results
+            """
+            orbratio collects 
+            phi_i(r1) phi_j(r1') phi_k(r2) phi_l(r2')/rho(r1') rho(r2')
+            """
+            phi_j_r1p = aux['orbs'][0][sweep][...,self._ijkl[1]]
+            phi_l_r2p = aux['orbs'][1][sweep][...,self._ijkl[3]]
+            rho1rho2 = 1.0/(fsum[0]*fsum[1])
+            #n is the walker number, i is the electron pair index, o is the orbital
+            orbratio = np.einsum("nio,nio,no,no,n ->nio", 
+                                 orb_configs[0][:,electrons_a_ind,:], #phi_i(r1)
+                                 orb_configs[1][:,electrons_b_ind,:], #phi_k(r2)
+                                 phi_j_r1p, #phi_j
+                                 phi_l_r2p, #phi_l
+                                 rho1rho2)
+
             results["value"] += np.einsum("in,inj->ij", wfratio, orbratio)
-            results["norm_a"] += norm_a[auxassignments_a[sweep]]
-            results["norm_b"] += norm_b[auxassignments_b[sweep]]
+            results["norm_a"] += norm[0]
+            results["norm_b"] += norm[1]
 
-        # Average over sweeps and pairs
         results["value"] /= self._nsweeps
         for e in ["a", "b"]:
             results["norm_%s" % e] /= self._nsweeps
@@ -232,8 +203,9 @@ def keys(self):
 
     def shapes(self):
         d = {"value": (self._ijkl.shape[1],), "ijkl": self._ijkl.T.shape}
+        nmo = self.orbitals.nmo()
         for e, s in zip(["a", "b"], self._spin_sector):
-            d["norm_%s" % e] = (self._orb_coeff[s].shape[1],)
+            d["norm_%s" % e] = (nmo[s],)
             d["acceptance_%s" % e] = ()
         return d
 
@@ -242,38 +214,6 @@ def avg(self, configs, wf):
         d["ijkl"] = self._ijkl.T
         return d
 
-    def get_extra_configs(self, configs):
-        """ Returns an nstep length array of configurations
-        starting from self._extra_config """
-        nconf = configs.configs.shape[0]
-
-        aux_configs_a = []
-        aux_configs_b = []
-        for step in range(self._nsweeps):
-            aux_configs_a.append(np.copy(self._aux_configs_a))
-            accept_a, self._aux_configs_a = sample_onebody(
-                self._mol,
-                self._orb_coeff[self._spin_sector[0]],
-                self._aux_configs_a,
-                tstep=self._tstep,
-            )
-            aux_configs_b.append(np.copy(self._aux_configs_b))
-            accept_b, self._aux_configs_b = sample_onebody(
-                self._mol,
-                self._orb_coeff[self._spin_sector[1]],
-                self._aux_configs_b,
-                tstep=self._tstep,
-            )
-        aux_configs_a = np.array(aux_configs_a)
-        aux_configs_b = np.array(aux_configs_b)
-
-        # Generates random choice of aux_config_a and aux_config_b for moving electron_a and electron_b
-        naux_a = self._aux_configs_a.shape[0]
-        naux_b = self._aux_configs_b.shape[0]
-        auxassignments_a = np.random.randint(0, naux_a, size=(self._nsweeps, nconf))
-        auxassignments_b = np.random.randint(0, naux_b, size=(self._nsweeps, nconf))
-        return [aux_configs_a, aux_configs_b], [auxassignments_a, auxassignments_b]
-
 
 def normalize_tbdm(tbdm, norm_a, norm_b):
     """Returns tbdm by taking the ratio of the averages in Eq. (10) of DOI:10.1063/1.4793531."""


Test patch:
diff --git a/pyqmc/testwf.py b/pyqmc/testwf.py
index 331800fe..e7ca6c66 100644
--- a/pyqmc/testwf.py
+++ b/pyqmc/testwf.py
@@ -197,12 +197,12 @@ def test_wf_gradient_laplacian(wf, configs):
         tt1 = time.time()
         tsep += ts1 - ts0
         ttog += tt1 - tt0
-        rel_grad = np.abs((andgrad - grad) / grad)
-        rel_lap = np.abs((andlap - lap) / lap)
-        rmax_grad = np.max(rel_grad)
-        rmax_lap = np.max(rel_lap)
+    rel_grad = np.abs((andgrad - grad) / grad)
+    rel_lap = np.abs((andlap - lap) / lap)
+    rmax_grad = np.max(rel_grad)
+    rmax_lap = np.max(rel_lap)
 
-    print("separate", tsep)
-    print("together", ttog)
+    print("time evaluated separately", tsep)
+    print("time evaluated together", ttog)
 
     return {"grad": rmax_grad, "lap": rmax_lap}
diff --git a/pytest.ini b/pytest.ini
new file mode 100644
index 00000000..d415aa40
--- /dev/null
+++ b/pytest.ini
@@ -0,0 +1,5 @@
+[pytest]
+markers =
+    slow: marks tests as slow (deselect with '-m "not slow"')
+    serial
+
diff --git a/tests/integration/test_complex_linemin.py b/tests/integration/test_complex_linemin.py
index a7f0f138..d8486f32 100644
--- a/tests/integration/test_complex_linemin.py
+++ b/tests/integration/test_complex_linemin.py
@@ -1,7 +1,7 @@
 import numpy as np
 import os
 import pyqmc
-
+import pytest
 
 def run_scf(chkfile):
     from pyscf import gto, scf
@@ -15,7 +15,7 @@ def run_scf(chkfile):
     mf = remove_linear_dep_(mf)
     energy = mf.kernel()
 
-
+@pytest.mark.slow
 def test():
     chkfile = "h2.hdf5"
     optfile = "linemin.hdf5"
diff --git a/tests/integration/test_dmc.py b/tests/integration/test_dmc.py
index fda9dbf4..8cc86189 100644
--- a/tests/integration/test_dmc.py
+++ b/tests/integration/test_dmc.py
@@ -5,17 +5,14 @@
 os.environ["MKL_NUM_THREADS"] = "1"
 os.environ["NUMEXPR_NUM_THREADS"] = "1"
 os.environ["OMP_NUM_THREADS"] = "1"
-import sys
 import numpy as np
-import pyqmc.testwf as testwf
-import pytest
 from pyqmc import reblock
+import pytest 
 
-
+@pytest.mark.slow
 def test():
     """ Ensure that DMC obtains the exact result for a hydrogen atom """
-    from pyscf import lib, gto, scf
-    from pyqmc.slater import PySCFSlater
+    from pyscf import gto, scf
     from pyqmc.jastrowspin import JastrowSpin
     from pyqmc.dmc import limdrift, rundmc
     from pyqmc.mc import vmc
@@ -23,13 +20,14 @@ def test():
     from pyqmc.func3d import CutoffCuspFunction
     from pyqmc.multiplywf import MultiplyWF
     from pyqmc.coord import OpenConfigs
+    from pyqmc import Slater
     import pandas as pd
 
     mol = gto.M(atom="H 0. 0. 0.", basis="sto-3g", unit="bohr", spin=1)
     mf = scf.UHF(mol).run()
     nconf = 1000
     configs = OpenConfigs(np.random.randn(nconf, 1, 3))
-    wf1 = PySCFSlater(mol, mf)
+    wf1 = Slater(mol, mf)
     wf = wf1
     wf2 = JastrowSpin(mol, a_basis=[CutoffCuspFunction(5, 0.2)], b_basis=[])
     wf2.parameters["acoeff"] = np.asarray([[[1.0, 0]]])
diff --git a/tests/integration/test_ecp.py b/tests/integration/test_ecp.py
index efef753b..a2160040 100644
--- a/tests/integration/test_ecp.py
+++ b/tests/integration/test_ecp.py
@@ -1,7 +1,7 @@
 import pandas as pd
 from pyqmc.mc import vmc, initial_guess
 from pyscf import gto, scf
-from pyqmc.slater import PySCFSlater
+from pyqmc import Slater
 from pyqmc.accumulators import EnergyAccumulator
 import numpy as np
 
@@ -11,7 +11,7 @@ def test_ecp():
     mol = gto.M(atom="C 0. 0. 0.", ecp="bfd", basis="bfd_vtz")
     mf = scf.RHF(mol).run()
     nconf = 5000
-    wf = PySCFSlater(mol, mf)
+    wf = Slater(mol, mf)
     coords = initial_guess(mol, nconf)
     df, coords = vmc(
         wf, coords, nsteps=100, accumulators={"energy": EnergyAccumulator(mol)}
diff --git a/tests/integration/test_multislater.py b/tests/integration/test_multislater.py
index 6fa3b30e..d37f2364 100644
--- a/tests/integration/test_multislater.py
+++ b/tests/integration/test_multislater.py
@@ -10,7 +10,7 @@
 import pyqmc.testwf as testwf
 from pyqmc.mc import vmc, initial_guess
 from pyqmc.accumulators import EnergyAccumulator
-from pyqmc.multislater import MultiSlater
+from pyqmc import Slater
 from pyqmc.coord import OpenConfigs
 from pyqmc.reblock import reblock
 import pyqmc
@@ -27,11 +27,11 @@ def test():
     delta = 1e-5
     nsteps = 200
     warmup = 10
-    for mf in [scf.RHF(mol).run(), scf.ROHF(mol).run(), scf.UHF(mol).run()]:
+    for mf in [scf.UHF(mol).run()]:#[scf.RHF(mol).run(), scf.ROHF(mol).run(), scf.UHF(mol).run()]:
         # Test same number of elecs
         mc = mcscf.CASCI(mf, ncas=4, nelecas=(1, 1))
         mc.kernel()
-        wf = MultiSlater(mol, mf, mc)
+        wf = Slater(mol, mf, mc)
 
         nconf = 10
 
@@ -61,7 +61,7 @@ def test():
         # Test different number of elecs
         mc = mcscf.CASCI(mf, ncas=4, nelecas=(2, 0))
         mc.kernel()
-        wf = MultiSlater(mol, mf, mc)
+        wf = Slater(mol, mf, mc)
 
         nelec = np.sum(mol.nelec)
         epos = initial_guess(mol, nconf)
diff --git a/tests/integration/test_obdm.py b/tests/integration/test_obdm.py
index 8b13722b..0f3b851b 100644
--- a/tests/integration/test_obdm.py
+++ b/tests/integration/test_obdm.py
@@ -5,13 +5,14 @@
 os.environ["OMP_NUM_THREADS"] = "1"
 
 import numpy as np
+np.random.seed(12534234)
 from pyscf import gto, scf, lo
 from numpy.linalg import solve
-from pyqmc import PySCFSlater
+from pyqmc import Slater
 from pyqmc.mc import initial_guess, vmc
 from pandas import DataFrame
 from pyqmc.obdm import OBDMAccumulator, normalize_obdm
-
+import pytest
 
 def test():
 
@@ -33,7 +34,7 @@ def test():
     nconf = 500
     nsteps = 400
     warmup = 15
-    wf = PySCFSlater(mol, mf)
+    wf = Slater(mol, mf)
     configs = initial_guess(mol, nconf)
     obdm_dict = dict(mol=mol, orb_coeff=lowdin, nsweeps=5, warmup=15)
     obdm = OBDMAccumulator(**obdm_dict)
@@ -60,7 +61,7 @@ def test():
     print(obdm_est["obdm_down"].diagonal().round(3))
     assert np.mean(np.abs(obdm_est["obdm_up"] + obdm_est["obdm_down"] - mfobdm)) < 0.05
 
-
+@pytest.mark.slow
 def test_pbc():
     from pyscf.pbc import gto, scf
     from pyqmc import supercell
@@ -99,7 +100,7 @@ def test_pbc():
     nconf = 800
     nsteps = 50
     warmup = 6
-    wf = PySCFSlater(mol, mf)
+    wf = Slater(mol, mf)
     configs = initial_guess(mol, nconf)
     obdm_dict = dict(mol=mol, orb_coeff=lowdin, kpts=kpts, nsweeps=4, warmup=10)
     obdm = OBDMAccumulator(**obdm_dict)
diff --git a/tests/integration/test_periodic.py b/tests/integration/test_periodic.py
index f79c115f..e483e478 100644
--- a/tests/integration/test_periodic.py
+++ b/tests/integration/test_periodic.py
@@ -114,7 +114,7 @@ def runtest(mol, mf, kind=0, do_mc=False):
             dm = np.sum(dm, axis=0)
     else:
         kpt = mf.kpts[kind]
-        wf = pyqmc.PySCFSlater(mol, mf)
+        wf = pyqmc.Slater(mol, mf)
         dm = mf.make_rdm1()
         print("original dm shape", dm.shape)
         if len(dm.shape) == 4:
diff --git a/tests/integration/test_shci_wf.py b/tests/integration/test_shci_wf.py
index 10084d78..df667128 100644
--- a/tests/integration/test_shci_wf.py
+++ b/tests/integration/test_shci_wf.py
@@ -5,7 +5,7 @@
 os.environ["OMP_NUM_THREADS"] = "1"
 import numpy as np
 import pyscf
-import pyqmc.multislater
+import pyqmc
 import pyscf.hci
 
 
@@ -16,7 +16,7 @@ def avg(vec):
     return avg, std / np.sqrt(nblock)
 
 
-def run_test():
+def test_shci_wf():
     mol = pyscf.gto.M(
         atom="O 0. 0. 0.; H 0. 0. 2.0",
         basis="ccecpccpvtz",
@@ -39,7 +39,7 @@ def run_test():
 
     tol = 0.0
     configs = pyqmc.initial_guess(mol, 1000)
-    wf = pyqmc.multislater.MultiSlater(mol, mf, cisolver, tol=tol)
+    wf = pyqmc.Slater(mol, mf, cisolver, tol=tol)
     data, configs = pyqmc.vmc(
         wf,
         configs,
@@ -55,4 +55,4 @@ def run_test():
 
 
 if __name__ == "__main__":
-    run_test()
+    test_shci_wf()
diff --git a/tests/integration/test_tbdm.py b/tests/integration/test_tbdm.py
index 61ce3ead..30d77a73 100644
--- a/tests/integration/test_tbdm.py
+++ b/tests/integration/test_tbdm.py
@@ -1,6 +1,6 @@
 import numpy as np
 from pyscf import gto, scf, lo
-from pyqmc import PySCFSlater
+from pyqmc import Slater
 from pyqmc.mc import initial_guess, vmc
 from pyqmc.accumulators import EnergyAccumulator
 from pyqmc.tbdm import TBDMAccumulator, normalize_tbdm
@@ -164,7 +164,7 @@ def test(atom="He", total_spin=0, total_charge=0, scf_basis="sto-3g"):
     tbdm_sweeps = 4
     tbdm_tstep = 0.5
 
-    wf = PySCFSlater(mol, mf)  # Single-Slater (no jastrow) wf
+    wf = Slater(mol, mf)  # Single-Slater (no jastrow) wf
     configs = initial_guess(mol, nconf)
     energy = EnergyAccumulator(mol)
     obdm_up = OBDMAccumulator(mol=mol, orb_coeff=iaos[0], nsweeps=tbdm_sweeps, spin=0)
@@ -216,17 +216,16 @@ def test(atom="He", total_spin=0, total_charge=0, scf_basis="sto-3g"):
     # Compares tbdm from QMC and MF
     avg_norm = {}
     avg_tbdm = {}
-    tbdm_est = {}
     for t in ["tbdm_upup", "tbdm_updown", "tbdm_downup", "tbdm_downdown"]:
         for k in df.keys():
             if k.startswith(t + "norm_"):
                 avg_norm[k.split("_")[-1]] = np.mean(df[k][vmc_warmup:], axis=0)
             if k.startswith(t + "value"):
                 avg_tbdm[k.split("_")[-1]] = np.mean(df[k][vmc_warmup:], axis=0)
-    for k in avg_tbdm:
-        tbdm_est[k] = normalize_tbdm(
+
+    tbdm_est = {k: normalize_tbdm(
             avg_tbdm[k].reshape(2, 2, 2, 2), avg_norm["a"], avg_norm["b"]
-        )
+        ) for k in avg_tbdm}
     qmctbdm = np.array(
         [
             [tbdm_est["upupvalue"], tbdm_est["updownvalue"]],
@@ -235,8 +234,9 @@ def test(atom="He", total_spin=0, total_charge=0, scf_basis="sto-3g"):
     )
     print("\nComparing QMC and MF tbdm:")
     for sa, sb in [[0, 0], [0, 1], [1, 0], [1, 1]]:
-        # print('QMC tbdm[%d,%d]:\n'%(sa,sb),qmctbdm[sa,sb])
-        # print('MF tbdm[%d,%d]:\n'%(sa,sb),mftbdm[sa,sb])
+        print("spins", sa, sb)
+        print('QMC tbdm[%d,%d]:\n'%(sa,sb),qmctbdm[sa,sb].round(3))
+        print('MF tbdm[%d,%d]:\n'%(sa,sb),mftbdm[sa,sb].round(3))
         diff = qmctbdm[sa, sb] - mftbdm[sa, sb]
         print("diff[%d,%d]:\n" % (sa, sb), diff)
         assert np.max(np.abs(diff)) < 0.05
diff --git a/tests/integration/test_twist.py b/tests/integration/test_twist.py
index 054311af..19eb4568 100644
--- a/tests/integration/test_twist.py
+++ b/tests/integration/test_twist.py
@@ -7,7 +7,7 @@
 from pyscf.pbc import gto, scf
 from pyscf.pbc.dft.multigrid import multigrid
 from pyscf.scf.addons import remove_linear_dep_
-
+import pyqmc
 
 def test_cubic_with_ecp(kind=1, nk=(2, 2, 2)):
     from pyscf.pbc.dft.multigrid import multigrid
@@ -78,14 +78,14 @@ def runtest(mol, mf, kind=0):
     twist = np.dot(kpt, mol.lattice_vectors().T / (2 * np.pi))
     print("kpt", kpt)
     print("twist", twist)
-    wf0 = pyqmc.PySCFSlater(mol, mf)
-    wft = pyqmc.PySCFSlater(mol, mf, twist=twist)
+    wf0 = pyqmc.Slater(mol, mf)
+    wft = pyqmc.Slater(mol, mf, twist=twist)
 
     #####################################
     ## compare values across boundary
     ## psi, KE, ecp,
     #####################################
-    nconfig = 100
+    nconfig = 50
     coords = pyqmc.initial_guess(mol, nconfig, 1)
     nelec = coords.configs.shape[1]
     epos, wrap = enforce_pbc(coords.lvecs, coords.configs)
@@ -113,6 +113,7 @@ def runtest(mol, mf, kind=0):
     assert np.linalg.norm(rat0 - 1) < 1e-10, rat0 - 1
     ratt = wft.testvalue(e, newcoords.electron(e))
     rattdiff = ratt - phase[:, e]
+    print("phase", phase[:,e])
     assert np.linalg.norm(rattdiff) < 1e-9, [
         np.round(rattdiff, 10),
         np.amax(np.abs(rattdiff)),
diff --git a/tests/integration/test_vmc.py b/tests/integration/test_vmc.py
index 20257e32..35a51a34 100644
--- a/tests/integration/test_vmc.py
+++ b/tests/integration/test_vmc.py
@@ -6,12 +6,13 @@
 import numpy as np
 import pandas as pd
 from pyqmc.mc import vmc, initial_guess
-from pyscf import lib, gto, scf
+from pyscf import gto, scf
 from pyqmc.reblock import reblock
-from pyqmc.slater import PySCFSlater
+from pyqmc import Slater
 from pyqmc.accumulators import EnergyAccumulator
+import pytest
 
-
+@pytest.mark.slow
 def test_vmc():
     """
     Test that a VMC calculation of a Slater determinant matches Hartree-Fock within error bars.
@@ -27,8 +28,8 @@ def test_vmc():
     warmup = 30
 
     for wf, mf in [
-        (PySCFSlater(mol, mf_rhf), mf_rhf),
-        (PySCFSlater(mol, mf_uhf), mf_uhf),
+        (Slater(mol, mf_rhf), mf_rhf),
+        (Slater(mol, mf_uhf), mf_uhf),
     ]:
 
         # Without blocks
@@ -70,19 +71,12 @@ def test_vmc():
 def test_accumulator():
     """ Tests that the accumulator gets inserted into the data output correctly.
     """
-    import pandas as pd
-    from pyqmc.mc import vmc, initial_guess
-    from pyscf import gto, scf
-    from pyqmc.energy import energy
-    from pyqmc.slater import PySCFSlater
-    from pyqmc.accumulators import EnergyAccumulator
-
     mol = gto.M(
         atom="Li 0. 0. 0.; Li 0. 0. 1.5", basis="cc-pvtz", unit="bohr", verbose=5
     )
     mf = scf.RHF(mol).run()
     nconf = 5000
-    wf = PySCFSlater(mol, mf)
+    wf = Slater(mol, mf)
     coords = initial_guess(mol, nconf)
 
     df, coords = vmc(
diff --git a/tests/integration/test_vmc_parallel.py b/tests/integration/test_vmc_parallel.py
index 7fe14edd..471901f0 100644
--- a/tests/integration/test_vmc_parallel.py
+++ b/tests/integration/test_vmc_parallel.py
@@ -1,6 +1,7 @@
 import pyqmc.recipes
 import concurrent.futures
 import os
+import pytest
 
 ncore = 2
 nconfig = ncore * 400
@@ -16,7 +17,7 @@ def run_scf(chkfile):
     mf.chkfile = chkfile
     mf.kernel()
 
-
+@pytest.mark.slow
 def test_parallel():
     run_scf("h2.hdf5")
     with concurrent.futures.ProcessPoolExecutor(max_workers=ncore) as client:
diff --git a/tests/unit/test_accumulators.py b/tests/unit/test_accumulators.py
index 0ab80530..0167da46 100644
--- a/tests/unit/test_accumulators.py
+++ b/tests/unit/test_accumulators.py
@@ -40,7 +40,7 @@ def test_info_functions_mol():
     from pyqmc.tbdm import TBDMAccumulator
 
     mol = gto.Mole()
-    mol.atom = """He 0.00 0.00 0.00 """
+    mol.atom = """Li 0.00 0.00 0.00; H 0. 0. 1.5 """
     mol.basis = "ccpvdz"
     mol.build()
 
diff --git a/tests/unit/test_derivatives.py b/tests/unit/test_derivatives.py
index 8c2d9a6a..1e12977b 100644
--- a/tests/unit/test_derivatives.py
+++ b/tests/unit/test_derivatives.py
@@ -16,7 +16,7 @@ def test_wfs():
     """
 
     from pyscf import lib, gto, scf
-    from pyqmc.slater import PySCFSlater
+    from pyqmc import Slater
     from pyqmc.multiplywf import MultiplyWF
     from pyqmc.manybody_jastrow import J3
     from pyqmc import default_jastrow
@@ -32,11 +32,11 @@ def test_wfs():
     for wf in [
         default_jastrow(mol)[0],
         J3(mol),
-        MultiplyWF(PySCFSlater(mol, mf), default_jastrow(mol)[0]),
-        MultiplyWF(PySCFSlater(mol, mf), default_jastrow(mol)[0], J3(mol)),
-        PySCFSlater(mol, mf_uhf),
-        PySCFSlater(mol, mf),
-        PySCFSlater(mol, mf_rohf),
+        MultiplyWF(Slater(mol, mf), default_jastrow(mol)[0]),
+        MultiplyWF(Slater(mol, mf), default_jastrow(mol)[0], J3(mol)),
+        Slater(mol, mf_uhf),
+        Slater(mol, mf),
+        Slater(mol, mf_rohf),
     ]:
         for k in wf.parameters:
             if k != "mo_coeff":
@@ -71,8 +71,7 @@ def test_pbc_wfs():
 
     from pyscf.pbc import lib, gto, scf
     from pyqmc.supercell import get_supercell
-    from pyqmc.slater import PySCFSlater
-    from pyqmc.multislaterpbc import MultiSlaterPBC
+    from pyqmc.slater import Slater
     from pyqmc.multiplywf import MultiplyWF
     from pyqmc import default_jastrow
     import pyqmc
@@ -91,19 +90,16 @@ def test_pbc_wfs():
     supercell = get_supercell(mol, S=(np.ones((3, 3)) - 2 * np.eye(3)))
     epos = pyqmc.initial_guess(supercell, nconf)
     # For multislaterpbc
-    kinds = 0, 3, 5, 6  # G, X, Y, Z
-    d1 = {kind: [0] for kind in kinds}
-    d2 = d1.copy()
-    d2.update({0: [], 3: [0, 1]})
-    detwt = [2 ** 0.5, 2 ** 0.5]
-    occup = [[d1, d2], [d1]]
-    map_dets = [[0, 1], [0, 0]]
+    #kinds = 0, 3, 5, 6  # G, X, Y, Z
+    #d1 = {kind: [0] for kind in kinds}
+    #d2 = d1.copy()
+    #d2.update({0: [], 3: [0, 1]})
+    #detwt = [2 ** 0.5, 2 ** 0.5]
+    #occup = [[d1, d2], [d1]]
+    #map_dets = [[0, 1], [0, 0]]
     for wf in [
-        MultiplyWF(PySCFSlater(supercell, mf), default_jastrow(supercell)[0]),
-        PySCFSlater(supercell, mf),
-        MultiSlaterPBC(supercell, mf, detwt=detwt, occup=occup, map_dets=map_dets),
-        # PySCFSlaterPBC(supercell, mf_uhf),
-        # PySCFSlaterPBC(supercell, mf_rohf),
+        MultiplyWF(Slater(supercell, mf), default_jastrow(supercell)[0]),
+        Slater(supercell, mf),
     ]:
         for k in wf.parameters:
             if "mo_coeff" not in k and k != "det_coeff":
diff --git a/tests/unit/test_ecp_cutoff.py b/tests/unit/test_ecp_cutoff.py
index aa1c9dfe..d5c4b9b9 100644
--- a/tests/unit/test_ecp_cutoff.py
+++ b/tests/unit/test_ecp_cutoff.py
@@ -1,10 +1,10 @@
 import pandas as pd
 from pyqmc.mc import vmc, initial_guess
 from pyscf import gto, scf, mcscf
-from pyqmc.slater import PySCFSlater
+from pyqmc import Slater
 from pyqmc.accumulators import EnergyAccumulator
 from pyqmc.multiplywf import MultiplyWF
-from pyqmc.multislater import MultiSlater
+from pyqmc import Slater
 from pyqmc import default_jastrow
 import numpy as np
 import time
@@ -25,9 +25,9 @@ def test_ecp():
     label = ["S", "J", "SJ"]
     ind = 0
     for wf in [
-        PySCFSlater(mol, mf),
+        Slater(mol, mf),
         default_jastrow(mol)[0],
-        MultiplyWF(PySCFSlater(mol, mf), default_jastrow(mol)[0]),
+        MultiplyWF(Slater(mol, mf), default_jastrow(mol)[0]),
     ]:
         wf.recompute(coords)
         print(label[ind])
@@ -43,7 +43,7 @@ def test_ecp():
 
     label = ["MS"]
     ind = 0
-    for wf in [MultiSlater(mol, mf, mc)]:
+    for wf in [Slater(mol, mf, mc)]:
         wf.recompute(coords)
         print(label[ind])
         ind += 1
diff --git a/tests/unit/test_line_minimization.py b/tests/unit/test_line_minimization.py
index 03cab8da..123f5726 100644
--- a/tests/unit/test_line_minimization.py
+++ b/tests/unit/test_line_minimization.py
@@ -8,9 +8,9 @@
 import pandas as pd
 from pyscf import lib, gto, scf
 from pyqmc import default_sj, line_minimization, initial_guess, gradient_generator
-import h5py
-
+import pytest
 
+@pytest.mark.slow
 def test():
     """ Optimize a Helium atom's wave function and check that it's 
     better than Hartree-Fock"""
diff --git a/tests/unit/test_overlap_derivative.py b/tests/unit/test_overlap_derivative.py
index 775de165..a303e7dc 100644
--- a/tests/unit/test_overlap_derivative.py
+++ b/tests/unit/test_overlap_derivative.py
@@ -1,4 +1,3 @@
-import copy
 import numpy as np
 from pyscf import gto, scf
 import pyqmc
@@ -67,7 +66,7 @@ def test():
     mf.stdout = None
 
     print("setting up wfs", flush=True)
-    wf0 = pyqmc.PySCFSlater(mol, mf)
+    wf0 = pyqmc.Slater(mol, mf)
     mf.mo_coeff[0][:, 0] = np.mean(mf.mo_coeff[0][:, :2], axis=1)
     wf1, to_opt = pyqmc.default_slater(mol, mf, optimize_orbitals=True)