wagnergroup / pyqmc Goto Github PK

It'd be nice to have a couple accumulator functions that give information about the data that will be produced. For example:

• keys()
• shapes -> return a dictionary of shapes that will be returned
• information -> any information that should be stored but does not depend on the wave function or electronic positions. For example, ijkl for tbdm.

File "/Users/lkwagner/work/h2_excited_states/functions.py", line 45, in vmc_gs
File "/Users/lkwagner/work/pyqmc/pyqmc/init.py", line 187, in generate_wf
File "/Users/lkwagner/work/pyqmc/pyqmc/init.py", line 187, in
File "/Users/lkwagner/work/pyqmc/pyqmc/init.py", line 90, in default_jastrow
index 0 is out of bounds for axis 0 with size 0

Occurs when na is set to zero.

Here is some code that exhibits the problem:

from pyscf import gto,scf 
import pyqmc

r=1.54/.529177
mol = gto.M(atom='H 0. 0. 0.; H 0. 0. %g'%r, ecp='bfd',basis='bfd_vtz',unit='bohr',verbose=1)
mf=scf.RHF(mol).run()
wf=pyqmc.slater_jastrow(mol,mf)
enacc=pyqmc.EnergyAccumulator(mol)
print(list(wf.parameters.keys()))

nconfig=500
coords=pyqmc.initial_guess(mol,nconfig)
wf,graddata=pyqmc.gradient_descent(wf,coords,warmup=1,
        params=['wf1mo_coeff','wf2acoeff','wf2bcoeff'],
                             accumulators={'pgrad':pyqmc.PGradAccumulator(enacc)},
                             step=0.1,maxiters=10)

When the number of parameters gets large, we start to see a slowdown:

1. Profile orbital optimization for a problem like benzene-ish in size.
2. It may be due to the fact that we construct the full overlap matrix, which is large. There are techniques in the literature to avoid construction of this matrix, look into using those.

This involves:

1. changing the ECP calculation so that it can return all possible moves.

2. Figuring out what to do with the accumulators and how to give DMC access to the ECP calculation.

For the line minimization algorithm, use correlated sampling to compute the minimum. This will improve the efficiency dramatically.

A simple algorithm is:
run VMC for a couple steps to get a new sample, with x0.
Then for each xprime, use an EnergyAccumulator to compute the energy and wave function. Average energy is given by sum( w E_L) )/ sum(w) where w is (psi(xprime)**2)/(psi(xprime))

We should have a new index on 'acoeff' so that the parameters can be different for each atom.

line_minimization and the recipes VMC and DMC can load in equilibrated configs from the optimization output instead of starting from initial_guess. If the number of configs is different on restarting, the saved configs can be sampled from it.

We need some basis sets that smoothly go to zero. Can use the ones that are used in QWalk. There is one with a cusp and one without (CUTOFF_CUSP and POLYPADE).

We need documentation for the helper functions in init

In particular, the parallel implementation does some unnecessary blocking and unblocking in the VMC part. We can have the VMC part just return the averages, similar to the regular VMC.

The tests that use VMC do not compute their uncertainties correctly, and so they fail too often.

Accumulators should be able to return metadata not in the returns.

IndexError: boolean index did not match indexed array along dimension 0; dimension is 1000 but corresponding boolean dimension is 47563810862832

Error message:

distributed.worker - WARNING - Compute Failed
Function: execute_task
args: ((<function apply at 0x2b4222f0be60>, <function vmc at 0x2b4235b7a710>, (<class 'tuple'>, [<pyqmc.slaterpbc.PySCFSlaterPBC object at 0x2b423b930c50>, <pyqmc.coord.PeriodicConfigs object at 0x2b423bc71b50>]), {'nsteps': 50, 'accumulators': {'energy': <pyqmc.accumulators.EnergyAccumulator object at 0x2b423bc74050>}, 'stepoffset': 150}))
kwargs: {}
Exception: IndexError('boolean index did not match indexed array along dimension 0; dimension is 1000 but corresponding boolean dimension is 47563810862832')

Traceback (most recent call last):
File "hfvmc.py", line 67, in
wf, configs, nsteps=4000, nsteps_per=50, accumulators={"energy": enacc}, hdf_file=hdf_file, verbose=True, client=client,
File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/dasktools.py", line 68, in distvmc
res = r.result()
File "/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/distributed/client.py", line 221, in result
six.reraise(*result)
File "/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/six.py", line 692, in reraise
raise value.with_traceback(tb)
File "/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/dask/compatibility.py", line 107, in apply
return func(*args, **kwargs)
File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/mc.py", line 170, in vmc
ratio = np.multiply(wf.testvalue(e, newcoorde) ** 2, t_prob)
File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/slaterpbc.py", line 239, in testvalue
aos = self.evaluate_orbitals(epos, mask)
File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/slaterpbc.py", line 142, in evaluate_orbitals
mycoords = mycoords[mask]
IndexError: boolean index did not match indexed array along dimension 0; dimension is 1000 but corresponding boolean dimension is 47563810862832

An accumulator object that accumulates 2-RDMs. Some required features:

• ability to only evaluate elements from a list of 2-RDM matrix elements

Implement a parallel DMC algorithm similar to the VMC implementation in parsltools.py

slaterpbc.py imports cubetools, which are not distributed anywhere, and there is no reason to distribute. It's not clear to me why we need generate_cube and/or plot_orbitals in slaterpbc. These are things that should be separate, or reproduce features already in pyscf.

Also, the testing code should be in tests.

ADAM was part of some tests, and is no longer necessary.

Depends on #73. Some constraints should be automatically set:

• Jastrow factor cusp conditions
• 'pivot' for orbital optimization
• 'pivot' for determinant weights

There should be a convenience function that sets up the determinant, Jastrow, etc with standard settings for an optimization, depending on whether we want to optimize various combinations.

This is a good first issue -- essentially implement the multiply object but with a sum with arbitary coefficients. Would be useful for creating samples for DMD.

The only trick is that you don't want to ever evaluate the total wave function. So for example, suppose that you want to evaluate the sum of two wave functions. They return s1=sign(psi), alpha1=ln(|psi|), s2, alpha2. We'd like to evaluate ln(s1e^alpha1+s2e^alpha2). But evaluating e^alpha1 or e^alpha2 may cause an overflow/underflow. Instead, we can evaluate it as ln(s1e^(alpha1-alpha2)+s2) + alpha2. If we choose things such that alpha2 is the larger of the two (or multiple), then the errors due to floating point, will be minimized, I think.

Similar formulae should be developed for the derivatives.

lm_sampler returns all quantities for all walkers, which is more than we need. Possibilities:

• Return (and compute) only energies (easy fix)
• Return averages only

This will help with optimizing lots of parameters when distributed parallelization is used.

testvalue(), etc need to be tested with masks.

The implementation of masks in the J3 wave function looks suspect. Sometimes there are divide by zero errors and something looks off..

This should be an object that:
init: uses a Cell object and initializes partial sums

ewald(configs): computes the ewald summed energy for a set of configurations.

QWalk has an implementation in Periodic_system which does the computation in O(N) time.

Right now there are only rules for 6 atomic integration points, should have at least 12

Implement code similar to the linear optimizer to make sure that the configurations are equilibriated before running optimization.

If beyond a cutoff, the Jastrow basis functions need not be computed. This should improve the speed quite a bit.

will need np.abs for the acceptance probability.

Matt Turk suggested we look into using cupy for GPU acceleration.

One way to test this would be to implement the Slater or Jastrow computation. Arrays get allocated on the GPU in evaluate() and then the input/output are the electron positions/wave function values and derivatives.

VMC and DMC are computing the gradient of the wave function twice per step. This can be reduced to once per step.

This is relatively simple; right now the code selects an electron at random; in the new implementation, select an up and a down electron at random, and average in separate matrices.

The name don't make sense.

In distances.py make a distance calculator that computes the minimal image.

I think most would agree that our current implementation with dasktools being separate from the Monte Carlo functions is awkward.

It's also a bit inefficient, since the wave functions get communicated to the workers every cycle.

This is a note to remind us to think about this design, and also to look into the load balancing that dask provides.

from pyscf import gto,scf   
import pyqmc
mol = gto.M(atom='O 0 0 0; H 0 -2.757 2.587; H 0 2.757 2.587',basis='bfd_vtz',ecp='bfd')     
mf = scf.RHF(mol).run()

wf=pyqmc.slater_jastrow(mol,mf)
acc=pyqmc.gradient_generator(mol,wf,['wf2acoeff','wf2bcoeff'])
nconf=1000
configs = pyqmc.initial_guess(mol, nconf)
wf,df=pyqmc.gradient_descent(wf,configs,acc,vmcoptions={'nsteps':30} )
dfdmc,configs,weights=pyqmc.dmc(wf,configs,accumulators={'energy':pyqmc.EnergyAccumulator(mol)},verbose=True)

Eventually starts to have troubles:

branch step 37
step 185 acceptance 0.988625 weight 1.002265253751489 weightstd 0.01820438950134978 eref -16.65434726384129
eref_mean -16.65434726384129 logwavg 0.002262691932260555 wmax 1.0758757453898113 wmin 0.8229786537756235
step 186 acceptance 0.989875 weight 1.0032058211691777 weightstd 0.02916024713709346 eref -16.65660995577355
eref_mean -16.65434726384129 logwavg 0.0032006934805391073 wmax 1.1302671678521425 wmin 0.7285751671920164
step 187 acceptance 0.9901249999999999 weight 1.0045723964471802 weightstd 0.03823234826351207 eref -16.65754795732183
eref_mean -16.65434726384129 logwavg 0.004561974798409794 wmax 1.166106805268598 wmin 0.6732892705719967
step 188 acceptance 0.990125 weight 1.0057451155420054 weightstd 0.046748909578166915 eref -16.6589092386397
eref_mean -16.65434726384129 logwavg 0.005728675303037307 wmax 1.1823796638716413 wmin 0.6573241150278046
step 189 acceptance 0.988625 weight 1.0069558278408637 weightstd 0.05513319068074363 eref -16.660075939144328
eref_mean -16.65434726384129 logwavg 0.0069317476709110275 wmax 1.2266356485676178 wmin 0.6500999395092298
branch step 38
step 190 acceptance 0.991 weight 1.0069842948273953 weightstd 0.012836025754139401 eref -16.708865614132332
eref_mean -16.708865614132332 logwavg 0.006960017614205829 wmax 1.0534934348009695 wmin 0.9031003104291069
step 191 acceptance 0.9898750000000001 weight 1.0070757972340736 weightstd 0.023825066022017866 eref -16.715825631746537
eref_mean -16.708865614132332 logwavg 0.007050881245437033 wmax 1.1194207577325734 wmin 0.8678979852350928
step 192 acceptance 0.98875 weight 1.0070303895909605 weightstd 0.034105183054996546 eref -16.71591649537777
eref_mean -16.708865614132332 logwavg 0.007005791623630759 wmax 1.201341353463013 wmin 0.8708230278973654
step 193 acceptance 0.9884999999999999 weight 1.0069453120083678 weightstd 0.04365548791801301 eref -16.715871405755962
eref_mean -16.708865614132332 logwavg 0.006921304424926493 wmax 1.2672532609547271 wmin 0.7943114765685841
step 194 acceptance 0.99 weight 1.0069492432723728 weightstd 0.05185078986959848 eref -16.71578691855726
eref_mean -16.708865614132332 logwavg 0.006925208565781031 wmax 1.3002273341306485 wmin 0.7962907775735976
branch step 39
step 195 acceptance 0.9886250000000001 weight 1.0079522474066884 weightstd 0.03386150757369647 eref -16.718853761207153
eref_mean -16.718853761207153 logwavg 0.007920794922536467 wmax 2.0138984865447456 wmin 0.9385844294912731
step 196 acceptance 0.9878750000000001 weight 1.010088703699885 weightstd 0.09814526131165632 eref -16.72677455612969
eref_mean -16.718853761207153 logwavg 0.010038152442170042 wmax 4.027796973089491 wmin 0.88815044664178
step 197 acceptance 0.99 weight 1.014528166180343 weightstd 0.2251286532797101 eref -16.728891913649324
eref_mean -16.718853761207153 logwavg 0.014423643506357081 wmax 8.055593946178982 wmin 0.8334937175602243
step 198 acceptance 0.9881249999999999 weight 1.0227932919623288 weightstd 0.47917092905759806 eref -16.73327740471351
eref_mean -16.718853761207153 logwavg 0.02253740590960254 wmax 16.111187892357965 wmin 0.8467052920379534
step 199 acceptance 0.989375 weight 1.0389576502131912 weightstd 0.9878564682334 eref -16.741391167116756
eref_mean -16.718853761207153 logwavg 0.038217951145075035 wmax 32.22237578471593 wmin 0.8242939303422971
branch step 40
step 200 acceptance 0.983375 weight 0.9584401203235264 weightstd 0.19029612390084147 eref -22.502547194938348
eref_mean -22.502547194938348 logwavg -0.0424481907238362 wmax 2.0779153004263824 wmin 0.8798714040620701
step 201 acceptance 0.983875 weight 0.9189299162832892 weightstd 0.5480598006378088 eref -22.460099004214513
eref_mean -22.502547194938348 logwavg -0.08454542037962559 wmax 4.155830600852765 wmin 0.7440853111114489
step 202 acceptance 0.98325 weight 0.9293413715121385 weightstd 1.153291476971568 eref -22.41800177455872
eref_mean -22.502547194938348 logwavg -0.07327914640393272 wmax 8.31166120170553 wmin 0.5907851588260982
step 203 acceptance 0.9835 weight 0.966349941387286 weightstd 1.9803071220355255 eref -22.429268048534414
eref_mean -22.502547194938348 logwavg -0.03422925220437351 wmax 16.62332240341106 wmin 0.4776376846535633
step 204 acceptance 0.9866250000000001 weight 1.0522604650215572 weightstd 3.3581503894845515 eref -22.468317942733975
eref_mean -22.502547194938348 logwavg 0.050940673994766095 wmax 33.24664480682212 wmin 0.35687993370174975
branch step 41
step 205 acceptance 0.932375 weight 0.915895005626507 weightstd 0.531035171957418 eref -42.182751568672224
eref_mean -42.182751568672224 logwavg -0.08785354355617378 wmax 2.1045209300431145 wmin 0.5653645796645618
step 206 acceptance 0.944375 weight 0.9700715797996006 weightstd 1.2554002479285529 eref -42.09489802511605
eref_mean -42.182751568672224 logwavg -0.030385416599487246 wmax 4.209041860086229 wmin 0.31978301574296153
step 207 acceptance 0.952875 weight 1.210007604066212 weightstd 2.389486928182594 eref -42.15236615207274
eref_mean -42.182751568672224 logwavg 0.1906266439411446 wmax 8.418083720172458 wmin 0.15842245571385768
step 208 acceptance 0.960125 weight 1.606131474789253 weightstd 4.153749320436266 eref -42.37337821261337
eref_mean -42.182751568672224 logwavg 0.47382847692430297 wmax 16.836167440344916 wmin 0.09501552555794791
step 209 acceptance 0.9670000000000001 weight 2.1704981170674973 weightstd 6.926497595539483 eref -42.65658004559653
eref_mean -42.182751568672224 logwavg 0.7749566882461915 wmax 33.67233488068983 wmin 0.05614032733754875
branch step 42
step 210 acceptance 0.836625 weight 1.9053564394715203 weightstd 1.417599642239037 eref -81.46036777849413
eref_mean -81.46036777849413 logwavg 0.6446690984104017 wmax 4.340996234134995 wmin 0.5479912974833153
step 211 acceptance 0.86575 weight 2.1642143837051573 weightstd 3.017803368006532 eref -82.10503687690453
eref_mean -81.46036777849413 logwavg 0.7720574243329175 wmax 8.68199246826999 wmin 0.1354797986220252
step 212 acceptance 0.897125 weight 2.7481142227548023 weightstd 5.319825118910757 eref -82.23242520282705
eref_mean -81.46036777849413 logwavg 1.0109149392731267 wmax 17.36398493653998 wmin 0.03469770105180635
step 213 acceptance 0.91825 weight 3.61073910441904 weightstd 9.047217728833163 eref -82.47128271776725
eref_mean -81.46036777849413 logwavg 1.283912489455028 wmax 34.72796987307996 wmin 0.008887489567370847
step 214 acceptance 0.934875 weight 4.768701978786448 weightstd 15.191135511762694 eref -82.74428026794916
eref_mean -81.46036777849413 logwavg 1.5620741460206768 wmax 69.45593974615991 wmin 0.0020757623685565703
branch step 43

When running the examples, cvmc_h2.py and parsl_h20.py an error is raised:

Traceback (most recent call last):
  File "cvmc_h2.py", line 111, in <module>
    sys["wf"], configs, accumulators={}, nsteps=40, nsteps_per=40, npartitions=ncore
  File "/mnt/home/bbusemeyer/programs/pyqmc/pyqmc/parsltools.py", line 64, in distvmc
    coord = np.split(coords, npartitions)
  File "/cm/shared/sw/pkg/devel/python3/3.7.3/lib/python3.7/site-packages/numpy/lib/shape_base.py", line 846, in split
    N = ary.shape[axis]
AttributeError: 'OpenConfigs' object has no attribute 'shape'

Around line 238, the code below computes the function for all r's, then zeros it out. This is slow and also can lead to floating point exceptions which can also be slow. We should fix it so it's like the other functions of this type.

        mask = r >= self.parameters["rcut"]
        r = r[..., np.newaxis]
        rvec = rvec
        z = r / self.parameters["rcut"]
        beta = self.parameters["beta"]

        p = z * z * (6 - 8 * z + 3 * z * z)
        dpdz = 12 * z * (z * z - 2 * z + 1)
        dbdp = -(1 + beta) / (1 + beta * p) ** 2
        dzdx = rvec / (r * self.parameters["rcut"])
        grad = dbdp * dpdz * dzdx
        d2pdz2_over_dpdz = (3 * z - 1) / (z * (z - 1))
        d2bdp2_over_dbdp = -2 * beta / (1 + beta * p)
        d2zdx2 = (1 - (rvec / r) ** 2) / (r * self.parameters["rcut"])
        lap = grad * (d2bdp2_over_dbdp * dpdz * dzdx + d2pdz2_over_dpdz * dzdx)
        lap += dbdp * dpdz * d2zdx2
        grad[mask] = 0
        lap[mask] = 0
        return grad, lap

In particular, test_pbc() is quite expensive, several factors more than the normal tests.

If it's not possible to speed it up, we could make an 'expensive' set of tests that are not run in the CI.

Make an adaptor object that wraps parameters derivatives.

init should be a representation of what parameters should be set equal or constant.

There should be three functions:

linearize_parameters:
input: dictionary of wave function parameters
output: one-dimensional numpy array of parameters

unlinearize_parameters:
input: 1-D array
output: dictionary of parameters

linearize_gradient:
input: dictionary of wave function parameter derivatives
output: one-dimensional array of reduced parameter derivatives.

A function that vmc, eval_ecp, obdm, and dmc can use that has
Arguments: lattice vectors, epos (nconfig,3) and translation (nconfig,3)
Returns: displaced epos and wraparound

the wraparound return gives which lattice vectors were used to wrap the electron back into the simulation cell. This can be combined with the twisted boundary conditions (k-point) to compute the overall phase of the wave function.

Brian reports that the OBDM evaluator doesn't work correctly. Check to see whether it does give correct results for a Slater determinant, in comparison to pyscf 1-RDM.

I'm running line_minimization on a large system (64 Si atoms) using the following setup

#!/usr/bin/env python

import numpy as np
import pyqmc
from pyqmc.dasktools import distvmc as vmc, line_minimization
from dask.distributed import Client, LocalCluster
import generate


if __name__=="__main__":
    slater = generate.get_slater(S=2 * (np.ones((3,3)) - 2 * np.eye(3)), excitations=[4])
    wf, pgrad = generate.make_sj_pgrad(slater, pyqmc.JastrowSpin, rcut=7.5)
    generate.load_wf_from_linemin("../linemin_conv_spin.hdf", wf, n=20, nprim=1)

    ncore = 20
    nconfig = ncore * 30
    configs = pyqmc.initial_guess(slater.supercell, nconfig, r=1.0)

    cluster = LocalCluster(n_workers=ncore, threads_per_worker=1)
    client = Client(cluster)

    # Warm up VMC
    df, configs = vmc(wf, configs, nsteps=3, verbose=True, client=client)
    print("Finished warmup VMC", flush=True)

    # Optimize jastrow
    hdf_file="exc4_linemin.hdf"
    lm_dict = dict(hdf_file=hdf_file, verbose=True, client=client, )
    lm_dict.update(dict(vmcoptions={"nsteps":40}), maxiters=50)
    wf, lm_df = line_minimization(wf, configs, pgrad, **lm_dict)

I get the following warnings and eventual error

/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/contextlib.py:119: UserWarning: Creating scratch directories is taking a surprisingly long time. This is often due to running workers on a network file system. Consider specifying a local-directory to point workers to write scratch data to a local disk.
  next(self.gen)
distributed.worker - WARNING - Memory use is high but worker has no data to store to disk.  Perhaps some other process is leaking memory?  Process memory: 2.48 GB -- Worker memory limit: 3.36 GB
/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/contextlib.py:119: UserWarning: Creating scratch directories is taking a surprisingly long time. This is often due to running workers on a network file system. Consider specifying a local-directory to point workers to write scratch data to a local disk.
  next(self.gen)
distributed.worker - WARNING - Memory use is high but worker has no data to store to disk.  Perhaps some other process is leaking memory?  Process memory: 2.57 GB -- Worker memory limit: 3.36 GB
... (many more lines like this) ...
distributed.worker - WARNING - Memory use is high but worker has no data to store to disk.  Perhaps some other process is leaking memory?  Process memory: 3.11 GB -- Worker memory limit: 3.36 GB
distributed.worker - WARNING - Memory use is high but worker has no data to store to disk.  Perhaps some other process is leaking memory?  Process memory: 3.29 GB -- Worker memory limit: 3.36 GB
distributed.nanny - WARNING - Worker exceeded 95% memory budget. Restarting
distributed.nanny - WARNING - Restarting worker
distributed.nanny - WARNING - Worker exceeded 95% memory budget. Restarting
distributed.nanny - WARNING - Restarting worker
distributed.worker - WARNING - Memory use is high but worker has no data to store to disk.  Perhaps some other process is leaking memory?  Process memory: 2.41 GB -- Worker memory limit: 3.36 GB
distributed.worker - WARNING - Memory use is high but worker has no data to store to disk.  Perhaps some other process is leaking memory?  Process memory: 2.65 GB -- Worker memory limit: 3.36 GB
distributed.worker - WARNING - Worker is at 85% memory usage. Pausing worker.  Process memory: 2.89 GB -- Worker memory limit: 3.36 GB
distributed.worker - WARNING - Memory use is high but worker has no data to store to disk.  Perhaps some other process is leaking memory?  Process memory: 2.89 GB -- Worker memory limit: 3.36 GB
distributed.worker - WARNING - Memory use is high but worker has no data to store to disk.  Perhaps some other process is leaking memory?  Process memory: 3.13 GB -- Worker memory limit: 3.36 GB
distributed.nanny - WARNING - Worker exceeded 95% memory budget. Restarting
Traceback (most recent call last):
  File "exc4_opt.py", line 29, in <module>
    wf, lm_df = line_minimization(wf, configs, pgrad, **lm_dict)
  File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/dasktools.py", line 171, in line_minimization
    return pyqmc.line_minimization(*args, vmc=distvmc, lm=dist_lm_sampler, **kwargs)
  File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/linemin.py", line 198, in line_minimization
    coords, last_en, pgrad, Sij, en, en_err = gradient_energy_function(x0, coords)
  File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/linemin.py", line 177, in gradient_energy_function
    data, coords = vmc(wf, coords, accumulators={"pgrad": pgrad_acc}, **vmcoptions)
  File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/dasktools.py", line 96, in distvmc
    allresults = list(zip(*[r.result() for r in runs]))
  File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/dasktools.py", line 96, in <listcomp>
    allresults = list(zip(*[r.result() for r in runs]))
  File "/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/distributed/client.py", line 221, in result
    six.reraise(*result)
  File "/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/six.py", line 703, in reraise
    raise value
distributed.scheduler.KilledWorker: ('vmc-8198f95b6e97c8f539177b02d94100a9', <Worker 'tcp://127.0.0.1:35107', memory: 0, processing: 15>)

It's clear that the problem is in running distvmc, and it seems to be in the PGradAccumulator. For now, I'm just documenting the problem.

The ECP calculation spends a bit too much time computing some things. In particular get_P_l is fairly slow and involves explicit loops in Python. This could be cleaned up for some speed improvements.

Timer unit: 1e-06 s

Total time: 0.013537 s
File: /Users/lkwagner/scr/pyqmc/pyqmc/eval_ecp.py
Function: ecp_ea at line 177

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   177                                           def ecp_ea(mol, configs, wf, e, at, threshold):
   178                                               """ 
   179                                               Returns the ECP value between electron e and atom at, local+nonlocal.
   180                                               """
   181         1          4.0      4.0      0.0      nconf = configs.configs.shape[0]
   182         1         12.0     12.0      0.1      ecp_val = np.zeros(nconf)
   183                                           
   184         1        952.0    952.0      7.0      l_list, v_l = get_v_l(mol, configs, e, at)
   185         1        183.0    183.0      1.4      mask, prob = ecp_mask(v_l, threshold)
   186         1         13.0     13.0      0.1      masked_v_l = v_l[mask]
   187         1         41.0     41.0      0.3      masked_v_l[:, :-1] /= prob[mask, np.newaxis]
   188         1         45.0     45.0      0.3      masked_configs = configs.mask(mask)
   189                                           
   190         1          1.0      1.0      0.0      naip = 6
   191         1          2.0      2.0      0.0      if len(l_list) > 2:
   192                                                   naip = 12
   193                                           
   194                                               # Use masked objects internally
   195         1       1257.0   1257.0      9.3      weights, epos_rot = get_rot(mol, masked_configs, e, at, naip)
   196         1       2798.0   2798.0     20.7      P_l = get_P_l(mol, masked_configs, weights, epos_rot, l_list, e, at)
   197                                           
   198                                               # Expand externally
   199         1         26.0     26.0      0.2      expanded_epos_rot = np.zeros((nconf, naip, 3))
   200         1         54.0     54.0      0.4      expanded_epos_rot[mask] = epos_rot
   201         1       8123.0   8123.0     60.0      ratio = get_wf_ratio(wf, configs, expanded_epos_rot, e, mask)
   202                                           
   203                                               # Compute local and non-local parts
   204         1         17.0     17.0      0.1      ecp_val[mask] = np.einsum("ij,ik,ijk->i", ratio, masked_v_l, P_l)
   205         1          1.0      1.0      0.0      local_l = -1
   206         1          7.0      7.0      0.1      ecp_val += v_l[:, local_l]
   207         1          1.0      1.0      0.0      return ecp_val

Related to closed #44, read that as well.

We would like to be able to constrain certain parameters to not change. This involves adding an option to LinearTransform, which allows selected parameters in a dictionary to be fixed. This will allow us to enforce cusp conditions (which we really want to do)

A simple method would be to allow a dictionary to be passed with a mask that is true for fixed and false for free parameters. This then gets converted into serialize_parameters and so on.

I ran VMC with dask (distvmc) using a HF slater determinant (PySCFSlaterPBC) with BFD ecps for 8 He atoms. An error was raised trying to invert a singular matrix, which followed division by zero warnings in the gradient and laplacian functions. Probably the wf was zero for the initial configuration, so determinant updates didn't work.

Running the attached file hfvmc_bug.py on the campus cluster with 20 processors resulted in the following error (full output file attached as well).

/projects/wagner/wawheel2/code/pyqmc/pyqmc/slaterpbc.py:269: RuntimeWarning: divide by zero encountered in true_divide
return ratios / testvalue
/projects/wagner/wawheel2/code/pyqmc/pyqmc/slaterpbc.py:269: RuntimeWarning: invalid value encountered in true_divide
return ratios / testvalue
/projects/wagner/wawheel2/code/pyqmc/pyqmc/slaterpbc.py:256: RuntimeWarning: divide by zero encountered in true_divide
return ratios[1:] / ratios[:1]
/projects/wagner/wawheel2/code/pyqmc/pyqmc/slaterpbc.py:256: RuntimeWarning: invalid value encountered in true_divide
return ratios[1:] / ratios[:1]
/projects/wagner/wawheel2/code/pyqmc/pyqmc/mc.py:85: RuntimeWarning: invalid value encountered in greater
mask = tot > cutoff
/projects/wagner/wawheel2/code/pyqmc/pyqmc/mc.py:86: RuntimeWarning: invalid value encountered in true_divide
g[mask, :] = cutoff * g[mask, :] / tot[mask, np.newaxis]
/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/numpy/core/_methods.py:75: RuntimeWarning: invalid value encountered in reduce
ret = umr_sum(arr, axis, dtype, out, keepdims)
distributed.worker - WARNING - Compute Failed
Function: execute_task
args: ((<function apply at 0x2ba25ae07e60>, <function vmc at 0x2ba26d6bc440>, (<class 'tuple'>, [<pyqmc.slaterpbc.PySCFSlaterPBC object at 0x2ba293a44b50>, <pyqmc.coord.PeriodicConfigs object at 0x2ba29451fa10>]), {'nsteps': 50, 'accumulators': {'energy': <pyqmc.accumulators.EnergyAccumulator object at 0x2ba294522090>}, 'stepoffset': 50}))
kwargs: {}
Exception: LinAlgError('Singular matrix')

Traceback (most recent call last):
File "hfvmc_bug.py", line 56, in
client=client, nsteps_per=50,
File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/dasktools.py", line 68, in distvmc
res = r.result()
File "/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/distributed/client.py", line 221, in result
six.reraise(*result)
File "/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/six.py", line 692, in reraise
raise value.with_traceback(tb)
File "/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/dask/compatibility.py", line 107, in apply
return func(*args, **kwargs)
File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/mc.py", line 151, in vmc
wf.recompute(configs)
File "/projects/wagner/wawheel2/code/pyqmc/pyqmc/slaterpbc.py", line 177, in recompute
self._inverse.append(np.linalg.inv(mo))
File "/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/numpy/linalg/linalg.py", line 532, in inv
ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
File "/projects/wagner/anaconda3/envs/pyscf/lib/python3.7/site-packages/numpy/linalg/linalg.py", line 89, in _raise_linalgerror_singular
raise LinAlgError("Singular matrix")
numpy.linalg.linalg.LinAlgError: Singular matrix

hfvmc_bug.py.txt
output.txt

This was Lucas's idea: the functionality we are using from dask is using futures to keep track of distributed jobs. We don't need the whole dask infrastructure just for this; we can implement our own simpler "client" with futures, so that pyqmc can also be used in parallel without dask.