HOWTOs -------------------- ### Periodic systems ``` def run_si_scf(chkfile="si_scf.chk", a=5.43): cell = gto.Cell( atom="Si 0. 0. 0.; Si {0} {0} {0}".format(a / 4), unit="angstrom", basis="ccecp-ccpvtz", ecp="ccecp", a=(np.ones((3, 3)) - np.eye(3)) * a / 2, ) cell.exp_to_discard = 0.1 cell.build() kpts = cell.make_kpts([8, 8, 8]) mf = scf.KRKS(cell, kpts=kpts) mf.chkfile = chkfile mf.run() ``` Run QMC on `n_conventional` conventional unit cells of silicon. ``` import pyqmc.api as pyq import numpy as np def run_si_qmc(chkfile="si_scf.chk", n_conventional=2): # Define periodic supercell in PyQMC conventional_S = np.ones((3, 3)) - 2 * np.eye(3) S = n_conventional * conventional_S pyq.OPTIMIZE(chkfile, "si_opt.chk", S=S) pyq.DMC(chkfile, "si_dmc.chk", load_parameters="si_opt.chk", S=S, slater_kws={'twist':0} ) ``` To get the available twists: ``` import pyqmc.api as pyq cell, mf = pyq.recover_pyqmc("si_scf.chk") conventional_S = np.ones((3, 3)) - 2 * np.eye(3) S = n_conventional * conventional_S cell = pyq.get_supercell(cell, S) twists=pyq.create_supercell_twists(cell,mf) print(twists) ``` ### Orbital optimization ``` pyq.OPTIMIZE(chkfile, "opt.chk",slater_kws={'optimize_orbitals':True,'optimize_zeros':False}) ``` ### Selected CI wave function Most of the effort is in setting up and saving the CI coefficients correctly, which is done in `run_hci()` here. You can copy `run_hci()` and use it on any system. If you are using pyscf 2.0 or above, you will need to run ``` pip install git+git://github.com/pyscf/naive-hci ``` to get the simple selected CI method. ``` import pyqmc.api as pyq import pyscf import pyscf.hci import numpy as np from pyscf import gto, scf def run_mf(chkfile): mol = gto.M( atom=""" O 0.0000000, 0.000000, 0.00000000 H 0.761561 , 0.478993, 0.00000000 H -0.761561, 0.478993, 0.00000000""", basis="ccecp-ccpvdz", ecp={"O": "ccecp"}, ) mf = scf.RHF(mol) mf.chkfile = chkfile mf.run() def run_hci(hf_chkfile, chkfile, select_cutoff=0.1, nroots=4): mol, mf = pyqmc.recover_pyscf(hf_chkfile, cancel_outputs=False) cisolver = pyscf.hci.SCI(mol) cisolver.select_cutoff=select_cutoff cisolver.nroots=nroots 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) pyscf.lib.chkfile.save(chkfile,'ci', {'ci':np.array(civec), 'nmo':nmo, 'nelec':nelec, '_strs':cisolver._strs, 'select_cutoff':select_cutoff, 'energy':e+mol.energy_nuc(), }) if __name__=="__main__": run_mf("mf.chk") run_hci("mf.chk","hci.chk") pyq.OPTIMIZE("mf.chk", "opt.chk", ci_checkfile="hci.chk", nconfig=1000) ``` ### CASCI/CASSCF wave function Note that when using CASCI/CASSCF, the `mc` object needs to have these four attributes: - `nelecas` - list of 2 ints - `ci` - matrix of ci coefficients - `ncas` - int - `mo_coeff` - array of orbital coefficients ``` def run_casci(scf_checkfile, ci_checkfile): mol, mf = pyq.recover_pyscf(scf_checkfile, cancel_outputs=False) mc = mcscf.CASCI(mf, 2, 2) mc.kernel() with h5py.File(ci_checkfile, "a") as f: f.create_group("ci") f["ci/ncas"] = mc.ncas f["ci/nelecas"] = list(mc.nelecas) f["ci/ci"] = mc.ci f["ci/mo_coeff"] = mc.mo_coeff return mc def run_casscf(scf_checkfile, ci_checkfile): mol, mf = pyq.recover_pyscf(scf_checkfile, cancel_outputs=False) mc = mcscf.CASSCF(mf, 2, 2) mc.chkfile = ci_checkfile mc.kernel() with h5py.File(mc.chkfile, "a") as f: f["mcscf/nelecas"] = list(mc.nelecas) f["mcscf/ci"] = mc.ci return mc if __name__ == "__main__": run_mf("mf.chk") run_casscf("mf.chk", "cas.chk") pyq.OPTIMIZE("mf.chk", "opt.chk", ci_checkfile="cas.chk", nconfig=1000) ``` You can also generate the wave function directly from the object without saving the results. ``` if __name__ == "__main__": mol = gto.M( atom="""H 0.0, 0.0, 0.0; H 0.7, 0.0, 0.0 """, basis="ccecp-ccpvdz", ) mf = scf.RHF(mol).run() mc = mcscf.CASSCF(mf, 2, 2) mc.kernel() wf, _ = pyq.generate_wf(mol, mf, mc=mc) ``` ### Create a wavefunction with 3 body Jastrow factor. ``` import pyscf.gto as gto import pyscf.scf as scf import pyqmc.api as pyq import pyqmc.wftools as wftools import pyqmc.mc as mc def linemin(mol,mf): #create wavefunction object. note jastrow_kws have to be passed even if empty wf, to_opt = wftools.generate_wf(mol, mf,jastrow=[wftools.generate_jastrow,wftools.generate_jastrow3],jastrow_kws=[{},{}] ) nconf = 100 configs = mc.initial_guess(mol, nconf) wf, dfgrad = pyq.line_minimization( wf, configs, pyq.gradient_generator(mol, wf, to_opt),max_iterations=50,verbose=True,hdf_file='3_jastrow.chk' ) def H2_ccecp_uhf(): r = 2 mol = gto.M( atom="H 0. 0. 0.; H 0. 0. %g" % r, ecp="ccecp", basis="ccpvdz", unit="bohr", verbose=1, ) mf = scf.UHF(mol).run() return mol, mf if __name__ == "__main__": mol,mf = H2_ccecp_uhf() linemin(mol,mf) ``` ### MPI parallelization You must put the execution in an `if __name__=="__main__"` block; otherwise `mpi4py` will not work. ``` from mpi4py.futures import MPIPoolExecutor import mpi4py.MPI import pyqmc.recipes if __name__=="__main__": comm = mpi4py.MPI.COMM_WORLD npartitions= comm.Get_size()-1 with MPIPoolExecutor(max_workers=npartitions) as client: pyqmc.recipes.OPTIMIZE("h2o.hdf5", "h2o_opt_mpi.hdf5", client=client, npartitions=npartitions) ``` One of the MPI ranks gets used as the main thread. So if you want to run it on 2 processors, you can do: ``` mpiexec -n 3 python -m mpi4py.futures test_mpi.py ``` ### Single node parallelization (without MPI) ``` from concurrent.futures import ProcessPoolExecutor if __name__=="__main__": npartitions = 2 with ProcessPoolExecutor(max_workers=npartitions) as client: pyqmc.recipes.OPTIMIZE("h2o.hdf5", "h2o_opt_mpi.hdf5", client=client, npartitions=npartitions) ``` ### Define a new accumulator This simple accumulator computes the average x,y, and z position of electrons. This is done in the `__call__` function; `configs.configs` is a numpy array of [walker, electron, xyz]. Note that if you run in parallel, you must define the class at the global scope. ``` import numpy as np class DipoleAccumulator: def __init__(self): pass def __call__(self, configs, wf): return {'electric_dipole':np.sum(configs.configs,axis=1) } def shapes(self): return {"electric_dipole": (3,)} def avg(self, configs, wf): d = {} for k, it in self(configs, wf).items(): d[k] = np.mean(it, axis=0) return d def keys(self): return self.shapes().keys() import pyqmc.recipes pyqmc.recipes.VMC("h2o.hdf5", "dipole.hdf5", load_parameters="h2o_sj_800.hdf5", accumulators={'extra_accumulators':{'dipole':DipoleAccumulator()}}) ``` ### Read in the 1-RDM You should have computed the 1-RDM by setting `accumulators={'rdm1':True}` in pyqmc.recipes.VMC or DMC. ``` import pyqmc.obdm def avg(vec): nblock = vec.shape[0] avg = np.mean(vec,axis=0) std = np.std(vec,axis=0) return avg, std/np.sqrt(nblock) with h5py.File("h2o_sj_vmc.hdf5") as f: warmup=2 en, en_err = avg(f['energytotal'][warmup:,...]) rdm1, rdm1_err=avg(f['rdm1value'][warmup:,...]) rdm1norm, rdm1norm_err = avg(f['rdm1norm'][warmup:,...]) rdm1=pyqmc.obdm.normalize_obdm(rdm1,rdm1norm) rdm1_err=pyqmc.obdm.normalize_obdm(rdm1_err,rdm1norm) ``` ### Compute the density from the 1-RDM ``` mol, mf = pyqmc.recover_pyscf("h2o.hdf5") import pyscf.tools ao_rdm1 = np.einsum('pi,ij,qj->pq', mf.mo_coeff, rdm1, mf.mo_coeff.conj()) resolution=0.05 dens=pyscf.tools.cubegen.density(mol, "h2o_sj_density.cube",ao_rdm1,resolution=resolution) ``` ### Plot the density from a cube file ``` import ase.io data = ase.io.cube.read_cube(open("h2o_sj_density.cube")) print(data.keys()) print(data['data'].shape) yzplane = data['data'][int(data['data'].shape[0]/2),:,:] fig = plt.figure(figsize=(8,8)) vmax=np.max(yzplane) plt.contourf(yzplane,levels=[vmax*x for x in np.logspace(-6,-1,80)], cmap='magma') plt.xticks([]) plt.yticks([]) ```