Repeat energy() with changed core integrals

tsommerfeld · December 17, 2021, 7:17pm

I would like to

Do a standard SCF.
Get the core Hamiltonian from the wfn object and modify it.
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?

hokru · December 17, 2021, 7:54pm

Hcore gets re-computed everytime.

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

tsommerfeld · December 18, 2021, 3:23pm

Yes, an external potential.

tsommerfeld · December 19, 2021, 6:14pm

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.

hokru · December 20, 2021, 10:32am

Maybe the easiest way is to modify the python psi4 driver running the scf:

github.com

psi4/psi4/blob/master/psi4/driver/procrouting/scf_proc/scf_iterator.py#L305

    
      
          
          
upcm = 0.0
          if core.get_option('SCF', 'PCM'):
              calc_type = core.PCM.CalcType.Total
              if core.get_option("PCM", "PCM_SCF_TYPE") == "SEPARATE":
                  calc_type = core.PCM.CalcType.NucAndEle
              Dt = self.Da().clone()
              Dt.add(self.Db())
              upcm, Vpcm = self.get_PCM().compute_PCM_terms(Dt, calc_type)
              SCFE += upcm
              self.push_back_external_potential(Vpcm)
          self.set_variable("PCM POLARIZATION ENERGY", upcm)  # P::e PCM
          self.set_energies("PCM Polarization", upcm)
          
          
upe = 0.0
          if core.get_option('SCF', 'PE'):
              Dt = self.Da().clone()
              Dt.add(self.Db())
              upe, Vpe = self.pe_state.get_pe_contribution(
                  Dt, elec_only=False
              )

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.