ZMP

The Zhao, Morrison, and Parr method works by updating self-consistently a potential generated by the difference between the current density and the target density.

\[v_c^i(r) = \lambda \int \frac{n_i(r') - n_{0}(r')}{|r-r'|} dr'\]

The lambda parameter is a Lagrange multiplier that when it goes to infinity, it allows the previous potential to become the exact Kohn-Sham Potential. In practice we won’t be able to increase lambda arbitrarily but values like 25, 50 or 75 are commonly used. An iterative approach also lets you perform an extrapolation technique iteratively. In order to trigger this calculation, specify more than one lambda. A reliable option to choose is to begin with lambda=10 and make enough iterations to achieve desired convergence in density.

The full method’s description can be found: From electron densities to Kohn-Sham kinetic energies, orbital energies, exchange-correlation potentials, and exchange-correlation energies Qingsheng Zhao, Robert C. Morrison, and Robert G. Parr Phys. Rev. A 50, 2138 – Published 1 September 1994

In n2v, we would can request the method in the following way. Let’s assume we want the exchange-correlation potential of a Neon Atom.

The algorithm will need an initial potential components to guide the inversion procedure. It is common to use the “Fermi Amaldi” (FA) Potential. Let us now use all of the previous information to perform the desired inversions.

Ne = psi4.geometry(
"""
0 1
Ne 0.0 0.0 0.0
noreorient
nocom
units bohr
symmetry c1
""" )


#n2v is driven by psi4's reference option. Make sure you set it accordingly.
psi4.set_options({"reference" : "rhf"})

#Perform a calculation for a target density.
#Remember that for post scf calculations, Psi4 does not update the density.
#Thus make sure you obtain something like a dipole in order to do so.
e, wfn_ne = psi4.properties("ccsd/cc-pcvqz", return_wfn=True, properties=["dipole"], molecule=Ne)

#Define inverter objects for each molcule. Simply use the wnf object from psi4 as an argument.
ine = n2v.Inverter(wfn_ne)

# Invert the density
  ine.invert("zmp", opt_max_iter=400, opt_tol=1e-7, zmp_mixing=1,
             lambda_list=[1,10,20,40], guide_potential_components=["fermi_amaldi"])

If we want to visualize the potential, we need to compute the potential according to:

\[v_{xc}(r) = v_c(r) - (1-1/N) \cdot v_{FA}(r)\]