Plane-Wave Mode

examples/pwscf/README.md

@@ -0,0 +1,35 @@
+# Plane-Wave Mode
+
+The `pyscf.pbc.pwscf` module provides experimental support for Hartree-Fock (HF),
+density functional theory (DFT), second-order Møller-Plesset perturbation theory (MP2),
+and coupled cluster singles doubles (CCSD) in a plane-wave basis.
+The CCECP and GTH pseudopotentials are supported for these methods,
+and SG15 pseudopotentials are supported for HF and DFT calculations.
+Occupation smearing and symmetry reduction of k-point meshes are implemented for HF and DFT.
+
+## Feature Overview
+
+The following self-consistent field (SCF) calculations are supported:
+* Hartree-fock (Restricted and Unrestricted)
+* Kohn-Sham DFT (Restricted and Unrestricted), with LDA, GGA, MGGA, and global hybrid functionals
+
+Currently, the Davidson algorithm is implemented for the effective Hamiltonian diagonalization.
+There are two mixing schemes for the effective potential, "Simple" and "Anderson" (DIIS).
+Symmetry reduction of k-points is supported for SCF calculations, along with occupation
+smearing. The default plane-wave basis set and integration grid are determined by `cell.mesh`,
+but these can be customized using the energy cutoffs `ecut_wf` and `ecut_rho` or by setting
+meshes directly using `PWKSCF.set_meshes()`.
+
+The following post-SCF calculations are supported:
+* MP2 (Restricted and Unrestricted)
+* CCSD (Restricted only)
+
+K-point symmetry and occupation smearing are currently not supported for post-SCF
+methods. The `PWKSCF.get_cpw_virtual()` method can be used to create virtual
+molecular orbitals in a GTO basis for use in post-SCF calculations.
+
+Plane-wave calculations can be performed with GTH or CCECP pseudopotentials
+(or all-electron for very small atoms). There is also basic support for SG15
+norm-conserving pseudopotentials. The post-SCF methods have been tested with GTH
+and CCECP but not SG15, while the SCF methods have been tested with GTH, CCECP, and SG15.
+

examples/pwscf/al.py

@@ -0,0 +1,38 @@
+from pyscf.pbc import gto
+from pyscf.pbc.pwscf.smearing import smearing_
+from pyscf.pbc.pwscf import kpt_symm
+import numpy as np
+
+"""
+Simple examples of running DFT for FCC Al. Uses
+k-point symmetrt and smearing
+"""
+
+cell = gto.Cell()
+cell.a = 4.0 * np.eye(3)
+x = 2.0
+cell.atom = f"""
+Al 0 0 0
+Al 0 {x} {x}
+Al {x} 0 {x}
+Al {x} {x} 0
+"""
+cell.pseudo = 'gth-pade'
+cell.basis = 'gth-szv'
+cell.verbose = 4
+cell.space_group_symmetry = True
+cell.symmorphic = True
+cell.build()
+
+kpts = cell.make_kpts(
+    [4, 4, 4],
+    time_reversal_symmetry=True,
+    space_group_symmetry=True
+)
+kmf = kpt_symm.KsymAdaptedPWKRKS(cell, kpts, ecut_wf=40)
+kmf = smearing_(kmf, sigma=0.01, method='gauss')
+kmf.xc = "PBE"
+kmf.nvir = 2
+kmf.conv_tol = 1e-7
+kmf.kernel()
+kmf.dump_scf_summary()

examples/pwscf/kccsd.py

@@ -0,0 +1,34 @@
+import numpy as np
+from pyscf.pbc import cc
+from pyscf.pbc.pwscf.pw_helper import gtomf2pwmf
+from pyscf.pbc.pwscf.kccsd_rhf import PWKRCCSD
+
+"""
+Simple CCSD calculation
+"""
+
+a0 = 1.78339987
+atom = "C 0 0 0; C %.10f %.10f %.10f" % (a0*0.5, a0*0.5, a0*0.5)
+a = np.asarray([
+        [0., a0, a0],
+        [a0, 0., a0],
+        [a0, a0, 0.]])
+
+from pyscf.pbc import gto, scf, pwscf
+cell = gto.Cell(atom=atom, a=a, basis="gth-szv", pseudo="gth-pade",
+                ke_cutoff=50)
+cell.build()
+cell.verbose = 5
+
+kpts = cell.make_kpts([2,1,1])
+
+mf = scf.KRHF(cell, kpts)
+mf.kernel()
+
+mcc = cc.kccsd_rhf.RCCSD(mf)
+mcc.kernel()
+
+pwmf = gtomf2pwmf(mf)
+pwmcc = PWKRCCSD(pwmf).kernel()
+
+assert(np.abs(mcc.e_corr - pwmcc.e_corr) < 1e-5)

examples/pwscf/kmp2.py

@@ -0,0 +1,40 @@
+from pyscf.pbc import gto, pwscf
+from pyscf.pbc.pwscf.kmp2 import PWKRMP2
+import numpy as np
+
+"""
+Simple MP2 calculation
+"""
+
+atom = "H 0 0 0; H 0.9 0 0"
+a = np.eye(3) * 3
+basis = "gth-szv"
+pseudo = "gth-pade"
+
+ke_cutoff = 50
+
+cell = gto.Cell(atom=atom, a=a, basis=basis, pseudo=pseudo,
+                ke_cutoff=ke_cutoff)
+cell.build()
+cell.verbose = 6
+
+nk = 2
+kmesh = [nk] * 3
+kpts = cell.make_kpts(kmesh)
+nkpts = len(kpts)
+
+pwmf = pwscf.PWKRHF(cell, kpts)
+pwmf.nvir = 20
+pwmf.kernel()
+
+es = {"5": -0.01363871, "10": -0.01873622, "20": -0.02461560}
+
+pwmp = PWKRMP2(pwmf)
+pwmp.kernel(nvir_lst=[5,10,20])
+pwmp.dump_mp2_summary()
+nvir_lst = pwmp.mp2_summary["nvir_lst"]
+ecorr_lst = pwmp.mp2_summary["e_corr_lst"]
+for nvir,ecorr in zip(nvir_lst,ecorr_lst):
+    err = abs(ecorr - es["%d"%nvir])
+    print(err)
+    assert(err < 1e-5)

examples/pwscf/kpt_symm.py

@@ -0,0 +1,56 @@
+from pyscf.pbc import gto
+from pyscf.pbc.pwscf.khf import PWKRHF
+from pyscf.pbc.pwscf.kpt_symm import KsymAdaptedPWKRHF
+import numpy as np
+import time
+
+"""
+Demonstrate the speedup of symmetry-adapted HF over
+non-symmetry-adapted HF.
+"""
+
+cell = gto.Cell(
+    atom = "C 0 0 0; C 0.89169994 0.89169994 0.89169994",
+    a = np.asarray([
+            [0.       , 1.78339987, 1.78339987],
+            [1.78339987, 0.        , 1.78339987],
+            [1.78339987, 1.78339987, 0.        ]]),
+    basis="gth-szv",
+    ke_cutoff=50,
+    pseudo="gth-pade",
+    verbose=4,
+    space_group_symmetry=True,
+    symmorphic=True,
+)
+cell.build()
+
+kmesh = [2, 2, 2]
+center = [0, 0, 0]
+kpts = cell.make_kpts(kmesh)
+skpts = cell.make_kpts(
+    kmesh,
+    scaled_center=center,
+    space_group_symmetry=True,
+    time_reversal_symmetry=True,
+)
+
+mf = PWKRHF(cell, kpts, ecut_wf=40)
+mf.nvir = 4
+t0 = time.monotonic()
+mf.kernel()
+t1 = time.monotonic()
+
+mf2 = KsymAdaptedPWKRHF(cell, skpts, ecut_wf=20)
+mf2.damp_type = "simple"
+mf2.damp_factor = 0.7
+mf2.nvir = 4
+t2 = time.monotonic()
+mf2.kernel()
+t3 = time.monotonic()
+
+print(mf.e_tot, mf2.e_tot)
+mf.dump_scf_summary()
+mf2.dump_scf_summary()
+print("nkpts in BZ and IBZ", skpts.nkpts, skpts.nkpts_ibz)
+print("Runtime without symmmetry", t1 - t0)
+print("Runtime with symmetry", t3 - t2)

examples/pwscf/li.py

@@ -0,0 +1,61 @@
+from pyscf.pbc import gto
+from pyscf.gto.basis import parse_cp2k_pp
+from pyscf.pbc.pwscf import kpt_symm, khf, smearing
+import numpy as np
+
+"""
+Simple examples of running HF and DFT for BCC Li.
+Both calculations run with smearing.
+The DFT calculation handles a larger k-mesh with
+k-point symmetry.
+"""
+
+cell = gto.Cell()
+a = 3.4393124531669552
+cell.a = a * np.eye(3)
+x = 0.5 * a
+cell.atom = f"""
+Li 0   0   0
+Li {x} {x} {x}
+"""
+cell.pseudo = {'Li': parse_cp2k_pp.parse("""
+#PSEUDOPOTENTIAL
+Li GTH2-HF-q1
+    1    0    0    0
+    0.75910286326041       2   -1.83343584669401    0.32295157976066
+       2
+    0.66792517034256       1    1.83367870276199
+    1.13098354939590       1   -0.00004141168540
+""")}
+cell.basis = 'gth-szv'
+cell.verbose = 4
+cell.space_group_symmetry = True
+cell.symmorphic = True
+cell.mesh = [10, 10, 10]
+cell.build()
+
+# Center at the Baldereshi point
+kpts = cell.make_kpts(
+    [2, 2, 2],
+    scaled_center=[1.0/6, 1.0/6, 0.5]
+)
+kmf = khf.PWKRHF(cell, kpts, ecut_wf=40)
+kmf = smearing.smearing_(kmf, sigma=0.02, method='gauss')
+kmf.xc = "PBE"
+kmf.conv_tol = 1e-7
+kmf.conv_tol_grad = 2e-3
+ehf = kmf.kernel()
+
+kpts = cell.make_kpts(
+    [4, 4, 4],
+    scaled_center=[1.0/6, 1.0/6, 0.5],
+    time_reversal_symmetry=True,
+    space_group_symmetry=True
+)
+kmf = kpt_symm.KsymAdaptedPWKRKS(cell, kpts, ecut_wf=40)
+kmf = smearing.smearing_(kmf, sigma=0.01, method='gauss')
+kmf.xc = "PBE"
+kmf.conv_tol = 1e-7
+kmf.conv_tol_grad = 2e-3
+ehf = kmf.kernel()
+

examples/pwscf/set_meshes.py

@@ -0,0 +1,32 @@
+from pyscf.pbc import gto
+from pyscf.pbc.pwscf.krks import PWKRKS
+import numpy as np
+
+"""
+This example demonstrates converging the energy with respect to
+the mesh size. Note that 
+"""
+
+cell = gto.Cell(
+    atom = "C 0 0 0; C 0.89169994 0.89169994 0.89169994",
+    a = np.asarray([
+            [0.       , 1.78339987, 1.78339987],
+            [1.78339987, 0.        , 1.78339987],
+            [1.78339987, 1.78339987, 0.        ]]),
+    basis="gth-szv",
+    ke_cutoff=50,
+    pseudo="gth-pade",
+)
+
+kmesh = [2, 1, 1]
+kpts = cell.make_kpts(kmesh)
+
+mf = PWKRKS(cell, kpts)
+
+# defaults
+print(cell.mesh, mf.wf_mesh, mf.xc_mesh)
+mf.kernel()
+
+mf.set_meshes(wf_mesh=[15, 15, 15], xc_mesh=[33, 33, 33])
+print(mf.wf_mesh, mf.xc_mesh)
+mf.kernel()

examples/pwscf/sg15.py

@@ -0,0 +1,82 @@
+from pyscf.pbc import gto
+from pyscf.pbc.pwscf.krks import PWKRKS
+# from pyscf.pbc.pwscf.kpt_symm import KsymAdaptedPWKRKS as PWKRKS
+from pyscf.pbc.pwscf.ncpp_cell import NCPPCell
+import numpy as np
+
+"""
+This example demonstrates the convergence of the total
+energy with respect to plane-wave energy cutoff for
+GTH and SG15 pseudopotentials. The SG15 converges
+faster, especially up to a 1000 eV cutoff (36.76 Ha),
+because these potentials were designed to converge more
+quickly.
+
+NOTE: Before using this example, you must set
+pbc_pwscf_ncpp_cell_sg15_path in your pyscf config file.
+"""
+
+kwargs = dict(
+    atom = "C 0 0 0; C 0.89169994 0.89169994 0.89169994",
+    a = np.asarray([
+            [0.       , 1.78339987, 1.78339987],
+            [1.78339987, 0.        , 1.78339987],
+            [1.78339987, 1.78339987, 0.        ]]),
+    basis="gth-szv",
+    ke_cutoff=50,
+    pseudo="gth-pade",
+    verbose=0,
+)
+
+cell = gto.Cell(**kwargs)
+cell.build()
+
+kwargs.pop("pseudo")
+nccell = NCPPCell(**kwargs)
+nccell.build()
+
+kmesh = [2, 2, 2]
+kpts = cell.make_kpts(kmesh)
+
+ens1 = []
+ens2 = []
+# A larger set of ecuts below is provided in case it's useful.
+# ecuts = [18.38235294, 22.05882353, 25.73529412, 29.41176471, 33.08823529,
+#          36.76470588, 44.11764706, 55.14705882, 73.52941176, 91.91176471]
+ecuts = [18.38235294, 25.73529412, 33.08823529, 36.76470588, 55.14705882]
+for ecut in ecuts:
+    print("ECUT", ecut)
+    # Run the GTH calculations
+    mf = PWKRKS(cell, kpts, xc="PBE", ecut_wf=ecut)
+    mf.damp_type = "simple"
+    mf.damp_factor = 0.7
+    mf.nvir = 4 # converge first 4 virtual bands
+    mf.kernel()
+    ens1.append(mf.e_tot)
+
+    # Run the SG15 calculations
+    mf2 = PWKRKS(nccell, kpts, xc="PBE", ecut_wf=ecut)
+    mf2.damp_type = "simple"
+    mf2.damp_factor = 0.7
+    mf2.nvir = 4 # converge first 4 virtual bands
+    mf2.init_pp()
+    mf2.init_jk()
+    mf2.kernel()
+    ens2.append(mf2.e_tot)
+    print(mf.e_tot, mf2.e_tot)
+    print()
+
+print()
+print("GTH Total Energies (Ha)")
+print(ens1)
+print("Energy cutoffs (Ha)")
+print(ecuts[:-1])
+print("Differences vs Max Cutoff (Ha)")
+print(np.array(ens1[:-1]) - ens1[-1])
+print()
+print("SG15 Total Energies (Ha)")
+print(ens2)
+print("Energy cutoffs (Ha)")
+print(ecuts[:-1])
+print("Differences vs Max Cutoff (Ha)")
+print(np.array(ens2[:-1]) - ens2[-1])

pyscf/lib/CMakeLists.txt

@@ -199,4 +199,5 @@ set_target_properties (clib_csf PROPERTIES
     OUTPUT_NAME "csf")
 
 add_subdirectory(sfnoci)
+add_subdirectory(pwscf)