Add examples for ASE interface (fix pyscf#624)

Author: sunqm

examples/geomopt/01-ase.py

@@ -0,0 +1,16 @@
+#!/usr/bin/env python
+
+'''
+Geometry optimization with ASE
+
+ASE uses the BFGS optimizer by default. This general-purpose algorithm is often
+not as efficient as dedicated geometry optimizers such as Berny or geomeTRIC.
+'''
+
+import pyscf
+
+mol = pyscf.M(atom='N 0 0 0; N 0 0 1.2', basis='ccpvdz')
+mf_opt = mol.RHF().Gradients().optimizer(solver='ase')
+
+mol_eq = mf_opt.kernel()
+print(mol_eq.tostring())

examples/pbc/09-ase_geometry_optimization.py

@@ -0,0 +1,47 @@
+#!/usr/bin/env python
+
+import pyscf
+from pyscf.pbc.tools.pyscf_ase import PySCF, cell_from_ase
+from ase.build import bulk
+from ase.optimize import BFGS
+from ase.filters import UnitCellFilter, StrainFilter
+
+atoms = bulk('Si', 'diamond', a=5.43)
+# Only .atoms and .a are defined in this cell instance. It's necessary to assign
+# basis set and pseudo potential
+cell = cell_from_ase(atoms)
+cell.basis = 'gth-dzvp'
+cell.pseudo = 'gth-pade'
+
+kmesh = [2, 2, 2]
+kpts = cell.make_kpts(kmesh)
+# Create a PySCF method, which will be used as a template in the following ASE
+# calculator. Parameters (such as convergence tol) can be defined in the
+# template.
+mf = cell.KRKS(xc='pbe', kpts=kpts)
+mf.conv_tol = 1e-6
+
+# PySCF(method=...) creates an ASE calculator for the specified PySCF method.
+# All methods implemented in PySCF can be used as templates to create a Calculator.
+# However, some methods may not support gradients or stress tensors. Geometry
+# optimization or lattice optimization that rely on these functionalities will
+# be halted.
+atoms.calc = PySCF(method=mf)
+
+# optimize atom positions
+opt = BFGS(atoms, logfile='atom_pos.log')
+opt.run()
+print(atoms.get_positions())
+
+
+# Optimize both lattice and atom positions
+opt = BFGS(UnitCellFilter(atoms), logfile='lattice_atom.log')
+opt.run()
+print(atoms.get_positions())
+print(atoms.cell)
+
+
+# Optimize lattice only. Atom positions (fractional coordinates) are frozen.
+opt = BFGS(StrainFilter(atoms), logfile='lattice_atom.log')
+opt.run()
+print(atoms.cell)

examples/pbc/09-talk_to_ase.py

@@ -38,7 +38,7 @@
 mf_dict = { 'xc' : 'lda,vwn' }
 
 # Once this is setup, ASE is used for everything from this point on
-ase_atom.set_calculator(pyscf_ase.PySCF(molcell=cell, mf_class=mf_class, mf_dict=mf_dict))
+ase_atom.calc = pyscf_ase.PySCF(method=mf)
 
 print("ASE energy", ase_atom.get_potential_energy())
 print("ASE energy (should avoid re-evaluation)", ase_atom.get_potential_energy())

examples/tools/08-ase_calculator.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python
+
+'''
+Using PySCF as an ASE Calculator
+
+Any PySCF mean-field or post-mean-field method (e.g., DFT, HF, CCSD) can be
+wrapped as an ASE-compatible calculator using the pyscf_ase interface.
+
+More examples can be found in examples/pbc/09-ase_geometry_optimization.py and
+examples/pbc/09-talk_to_ase.py
+'''
+
+import pyscf
+from pyscf.pbc.tools.pyscf_ase import PySCF, ase_atoms_to_pyscf
+from ase import Atoms
+
+# Define molecules using ase.Atoms class, the ase_atoms_to_pyscf function can
+# convert the geometry to pyscf input format
+atoms = Atoms('N2', [(0, 0, 0), (0, 0, 1.2)])
+mol = pyscf.M(atom=ase_atoms_to_pyscf(atoms), basis='ccpvdz')
+
+# Example 1: Use a DFT method as an ASE calculator
+mf = mol.RKS(xc='pbe0')
+atoms.calc = PySCF(method=mf)
+
+# Example 2: Use a post-HF method (CCSD) as an ASE calculator
+atoms.calc = PySCF(method=mol.RHF().CCSD())

pyscf/pbc/tools/pyscf_ase.py

@@ -72,7 +72,7 @@ def cell_from_ase(ase_atoms):
 
 class PySCF(Calculator):
     implemented_properties = ['energy', 'forces', 'stress',
-                              'dipole', 'magmom']
+                              'dipole', 'magmom', 'polarizability']
 
     default_parameters = {}
 
@@ -169,6 +169,20 @@ def calculate(self, atoms=None, properties=['energy'],
             # in Gaussian cgs unit
             self.results['dipole'] = base_method.dip_moment() * Debye
 
+        if 'polarizability' in properties:
+            assert hasattr(base_method, 'istype') and base_method.istype('SCF'), \
+                    'Polarizability can only be computed with mean-field methods'
+            if self.pbc:
+                from pyscf.pbc.prop.polarizability.rhf import Polarizability
+                p = Polarizability(base_method).polarizability()
+            else:
+                from pyscf.prop.polarizability import rhf, uhf
+                if base_method.istype('UHF'):
+                    p = uhf.Polarizability(base_method).polarizability()
+                else:
+                    p = rhf.Polarizability(base_method).polarizability()
+            self.results['polarizability'] = p * (BOHR**3)
+
         if 'magmom' in properties:
             magmom = self.mol.spin
             self.results['magmom'] = magmom