Maybe the easiest way is to modify the python psi4 driver running the scf:
You would need to convert your numpy matrix to a psi4.Matrix (see example here: Interface to NumPy) and the you can add it like in the linked-line above.
At the start I would avoid symmetry.
You can also do self.H().add(MyPotential).
It may be useful to expose an interface that does not require modifying the python code. At least I don’t think we have one yet.
Last week we came finally back to this project. [Student graduated, project mothballed, …] But with Jonathan’s help, we were able to implement a psithon solution (see below).
Suggestions for improvement and hint regarding pitfalls with our approach greatly appreciated.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jul 18 15:58:46 2024
@author: Thomas Sommerfeld
Psi4numpy implementation for adding an arbitrary one-electron potential
to the one-electron Hamiltonian:
H_lambda = H + lambda*U
where U is an artifical - normally stabilizing - potential, and
lambda is a scaling parameter.
Example:
U = V(el-nuc) [because it is already there]
lambda = 0.01
Effectively, this scales the nuclear charges by 1.01.
Purpose:
Compute energies of temporary anions by analytical
continuation; see DOI 10.1088/1402-4896/ab7941
Remark:
In the provided example only V(el-nuc) is scaled,
because V(nuc-nuc) cancels when computing vertical electron binging energies.
Computed energies in the example below:
ESCF_unscaled (lambda = 0)
EMP2_unscaled (lambda = 0)
ESCF_Z_scaled (lambda = 0.01)
EMP2_Z_scaled (lambda = 0.01)
"""
import psi4
import numpy as np
"""
Memory and output file notebook style
"""
Mega=1e6
psi4.set_memory(int(500*Mega))
numpy_memory = 2 # GB
psi4.core.set_output_file('output.dat', False)
"""
Build base-wavefunction
"""
mol = psi4.geometry("""
N
N 1 1.1
symmetry D2h
""")
psi4.set_options({'basis': 'aug-cc-pvdz',
'scf_type': 'pk',
'e_convergence': 1e-8})
base_wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option('BASIS'))
"""
useful for arbitrary U
n_irrep = base_wfn.nirrep()
print(f'n_Irreps = {n_irrep}')
irreps = mol.irrep_labels()
print(irreps)
"""
"""
SCF and MP2 without any artificial potential
save occupation for SCF with changed V(el-nuc)
"""
scf_wfn = psi4.driver.scf_wavefunction_factory('scf', base_wfn, psi4.core.get_option('SCF', 'REFERENCE'))
ESCF_unscaled = scf_wfn.compute_energy()
EMP2_unscaled = psi4.energy('mp2', ref_wfn=scf_wfn)
doc=scf_wfn.doccpi()
soc=scf_wfn.soccpi()
"""
SCF and MP2 with V(el-nuc) *= scaling_factor
* start with fresh ref-wfn
* force occupation
* initialize, so the one-electron Hamiltonian H exists
* get U - in SO rep. - from the integrals helper
* scale U
* add U to H
* do the scf iteration and compute the scf energy
* compute the MP2 energy using the converged scf
"""
scf_wfn = psi4.driver.scf_wavefunction_factory('scf', base_wfn, psi4.core.get_option('SCF', 'REFERENCE'))
scf_wfn.force_occpi(doc, soc)
scf_wfn.initialize()
mints = psi4.core.MintsHelper(scf_wfn.basisset())
U = mints.so_potential()
U.scale(0.01)
scf_wfn.H().add(U)
scf_wfn.iterations()
ESCF_Z_scaled=scf_wfn.finalize_energy()
EMP2_Z_scaled=psi4.energy('mp2', ref_wfn=scf_wfn)
print('lambda E(SCF) E(MP2)')
print(f' 0 {ESCF_unscaled} {EMP2_unscaled}')
print(f' 0.01 {ESCF_Z_scaled} {EMP2_Z_scaled}')