Repeat energy() with changed core integrals

I would like to

  1. Do a standard SCF.
  2. Get the core Hamiltonian from the wfn object and modify it.
  3. Do a SCF with the modified core Hamiltonian.

Just as an example, I was hoping the following would return a different energy, but it doesn’t. Still, you see where I’m going:

h2 = psi4.geometry("""
0 1
H
H 1 0.78
symmetry c1
“”")
e, wf = psi4.energy(‘scf’, return_wfn=True)
m=np.asarray(wf.H())
m[:,:]*=0.5
psi4.energy(‘scf’)

It seems psi4.energy() computes its integrals anew. Is there any way around that?

Hcore gets re-computed everytime.

Is this like an external potential you’d like to add or a more fundamental change of H?

Yes, an external potential.

Seems external potential solutions exist.

  • Holger has a solution that was discussed before.
  • there is a build-in EXTERN option that seems to be able to read a potential from a grid (ONEPOT_GRID_READ)

Yet, the documentation on ONEPOT_GRID_READ is sparse- it even has a question mark- and I wasn’t able to find an example.

I’d greatly appreciate if anyone could share an example.

  • Establishing the grid (not dense at the nuclei, but away from the nuclei- it’s an external potential).
  • Formatting the “.dx” file.

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.

Very late thank you.

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}')

This topic was automatically closed 60 days after the last reply. New replies are no longer allowed.