Skip to content

Commit

Permalink
more documentation updates.
Browse files Browse the repository at this point in the history
  • Loading branch information
lkwagner committed Jan 3, 2021
1 parent 91028ab commit b105e4c
Show file tree
Hide file tree
Showing 3 changed files with 85 additions and 0 deletions.
1 change: 1 addition & 0 deletions doc/source/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ Welcome to pyQMC's documentation!

install
tutorial
howto
specific_instructions
developer

Expand Down
49 changes: 49 additions & 0 deletions doc/source/reduced_density_matrix_snippets.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,49 @@
Reduced density matrix snippets
------------------------------------


### 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([])
```
35 changes: 35 additions & 0 deletions doc/source/snippets.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
Snippets
--------------------

### 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()
```

```
import pyqmc.recipes
import numpy as np
def run_si_qmc(chkfile="si_scf.chk"):
# Define periodic supercell in PyQMC
conventional_S = np.ones((3, 3)) - 2 * np.eye(3)
S = 2 * conventional_S
pyqmc.recipes.OPTIMIZE(chkfile, "si_opt.chk", S=S)
pyqmc.recipes.DMC(chkfile, "si_dmc.chk", start_from="si_opt.chk", S=S)
```

0 comments on commit b105e4c

Please sign in to comment.