<> Navigation: * [[AJMPublic/camcasp|CamCASP]] = IP and AC-shift calculations with Psi4 = The following Psi4 code can be used to easily calculate vertical ionization potentials (IPs) and the shifts needed for the asymptotic correction (AC-Shift). Since we do this sort of calculation often and it can be error-prone, I wrote the following code. It was also a small challenge to me to do everything in one calculation, and use some of the Python constructs in Psi4. This is an example for a molecule taken from the IL16 database. To enter your molecule just alter the contents of the '''molecule''' block. Make sure you set the correct charge and multiplicity. The code should deal with charged systems correctly. You probably want to choose a better basis too. The code defaults to the PBE0 functional. This is a good choice for IPs and is recommended. Though some of the new meta-GGAs (the "//ab initio//" ones) may prove better. [[attachment:IL150-IP-calculation.psi4]] {{{ memory 4000 mb molecule monomer { 1 1 # charge and multiplicity. In this case the molecule has charge +1 units angstrom C1 1.25456 0.04517 -0.22319 N2 -0.03032 0.04506 0.31669 C3 0.07512 0.04497 1.66566 N4 1.36986 0.04500 1.99001 C5 2.12628 0.04514 0.83592 C6 -1.29617 0.04507 -0.43477 H7 1.68755 0.04500 2.95634 H8 -0.72290 0.04489 2.39723 H9 3.20705 0.04520 0.85849 H10 1.42695 0.04526 -1.29040 H11 -1.07074 0.04515 -1.50429 H12 -1.87711 -0.85035 -0.18843 H13 -1.87716 0.94042 -0.18831 # Unfortunately, Psi4 1.x does not write out the orbital energies # correctly if symmetry is used. So we disable it here by forcing C1. # This will have the consequence that the calculation will proceed # more slowly, and, for symmetric systems, the calculation may even # fail to converge. If symmetry is enabled then the e_HOMO and AC-shift # values may be incorrect and these should be checked and corrected manually. Symmetry C1 } # Define the XC functional here dft_func = 'pbe0' set basis aug-cc-pVDZ set freeze_core TRUE # get the charge and multiplicity of the monomer charge = monomer.molecular_charge() mult = monomer.multiplicity() set reference rks EN, wfn = energy(dft_func, return_wfn=True) # We need to find the HOMO eigenvalue of this system nelec_a = wfn.nalpha() epsilon_a = wfn.epsilon_a().np # The .np is need to convert the Psi4 Vector object to an np.array e_homo = epsilon_a[nelec_a-1] # now remove and electron and increase the multiplicity. monomer.set_molecular_charge(charge+1) monomer.set_multiplicity(mult+1) set reference uks EN1 = energy(dft_func) IP = EN1 - EN AC_shift = IP + e_homo psi4.print_out("\n") psi4.print_out("%s energies and IP \n\n" % (dft_func)) psi4.print_out(" Num Electrons Charge E_int [a.u.] \n") psi4.print_out("-----------------------------------------------------\n") psi4.print_out(" N %3d %10.6f \n" % (charge,EN)) psi4.print_out(" N-1 %3d %10.6f \n" % (charge+1,EN1)) psi4.print_out("-----------------------------------------------------\n") psi4.print_out(" IP %10.6f a.u. \n" % (IP)) psi4.print_out(" e_HOMO %10.6f a.u. \n" % (e_homo)) psi4.print_out(" AC-shift %10.6f a.u. \n" % (AC_shift)) psi4.print_out("-----------------------------------------------------\n") }}} The important part of the output will be at the bottom and should look like: [[attachment:IL150 snippet]] {{{ PBE0 energies and IP Num Electrons Charge E_int [a.u.] ----------------------------------------------------- N 1 -265.631128 N-1 2 -265.094762 ----------------------------------------------------- IP 0.536366 a.u. e_HOMO -0.457412 a.u. AC-shift 0.078954 a.u. ----------------------------------------------------- }}} == Known issues == * For systems with symmetry (atoms) this code is known to pick the wrong energy of the HOMO level. I believe it is a bug in Psi4 1.x. as the orbital energy array misses certain energies associated with symmetry-related degenerate levels. * With newer versions of Psi4 the wfn.epsilon_a() lists are not available directly. The fix is to use '''wfn.epsilon_a().np''' to obtain the np.array of up-spin orbital energies. * The IP could be the computed one (listed at the bottom of the output), or the experimental one. Make sure everything is in atomic units. See the example output after this input file. [[attachment:Ar_IP-calculation_Psi4_1.3_version.psi4]] {{{ memory 4 gb molecule monomer { 0 1 # charge and multiplicity. units angstrom Ar 0.0 0.0 0.0 # Unfortunately, Psi4 1.x does not write out the orbital energies # correctly if symmetry is used. So we disable it here by forcing C1. # This will have the consequence that the calculation will proceed # more slowly, and, for symmetric systems, the calculation may even # fail to converge. If symmetry is enabled then the e_HOMO and AC-shift # values may be incorrect and these should be checked and corrected manually. Symmetry C1 } # Define the XC functional here dft_func = 'pbe0' set basis aug-cc-pVTZ set freeze_core TRUE # get the charge and multiplicity of the monomer charge = monomer.molecular_charge() mult = monomer.multiplicity() set reference rks EN, wfn = energy(dft_func, return_wfn=True) # We need to find the HOMO eigenvalue of this system nelec_a = wfn.nalpha() epsilon_a = wfn.epsilon_a().np e_homo = epsilon_a[nelec_a-1] # now remove and electron and increase the multiplicity. monomer.set_molecular_charge(charge+1) monomer.set_multiplicity(mult+1) set reference uks EN1 = energy(dft_func) IP = EN1 - EN AC_shift = IP + e_homo psi4.print_out("\n") psi4.print_out("%s energies and IP \n\n" % (dft_func)) psi4.print_out(" Num Electrons Charge E_int [a.u.] \n") psi4.print_out("-----------------------------------------------------\n") psi4.print_out(" N %3d %10.6f \n" % (charge,EN)) psi4.print_out(" N-1 %3d %10.6f \n" % (charge+1,EN1)) psi4.print_out("-----------------------------------------------------\n") psi4.print_out(" IP %10.6f a.u. \n" % (IP)) psi4.print_out(" e_HOMO %10.6f a.u. \n" % (e_homo)) psi4.print_out(" AC-shift %10.6f a.u. \n" % (AC_shift)) psi4.print_out("-----------------------------------------------------\n") }}} Here are the relevant parts of the output file. The first part is the energy and orbital energies for the species with $N$ electrons (in this case, the neutral Ar atom). We need to use this and not the data for the species with $N-1$ electrons. This is very important! The last part is the summary table produced by the above code. It will be at the end of the output file. {{{ ... ... ==> Geometry <== Molecular point group: c1 Geometry (in Angstrom), charge = 0, multiplicity = 1: Center X Y Z Mass ------------ ----------------- ----------------- ----------------- ----------------- AR 0.000000000000 0.000000000000 0.000000000000 39.962383123700 Running in c1 symmetry. Rotational constants: A = '''''''''''''''''' B = '''''''''''''''''' C = '''''''''''''''''' [cm^-1] Rotational constants: A = '''''''''''''''''' B = '''''''''''''''''' C = '''''''''''''''''' [MHz] Nuclear repulsion = 0.000000000000000 Charge = 0 Multiplicity = 1 Electrons = 18 Nalpha = 9 Nbeta = 9 ... ... SCF Guess: Core (One-Electron) Hamiltonian. ==> Iterations <== Total Energy Delta E RMS |[F,P]| @DF-RKS iter 1: -511.52636091434454 -5.11526e+02 1.20735e-01 DIIS @DF-RKS iter 2: -516.43497634350854 -4.90862e+00 8.70805e-02 DIIS @DF-RKS iter 3: -526.80716133498390 -1.03722e+01 2.38723e-02 DIIS @DF-RKS iter 4: -527.38114680046090 -5.73985e-01 1.15418e-03 DIIS @DF-RKS iter 5: -527.38145757657207 -3.10776e-04 4.30197e-04 DIIS @DF-RKS iter 6: -527.38168212177584 -2.24545e-04 3.76718e-05 DIIS @DF-RKS iter 7: -527.38168337388311 -1.25211e-06 2.39707e-06 DIIS @DF-RKS iter 8: -527.38168337947411 -5.59100e-09 5.69820e-08 DIIS Energy and wave function converged. ==> Post-Iterations <== Orbital Energies [Eh] --------------------- Doubly Occupied: 1A -115.306790 2A -11.219611 3A -8.742421 4A -8.742421 5A -8.742421 6A -0.992776 7A -0.440597 8A -0.440597 9A -0.440597 Virtual: 10A 0.063800 11A 0.084352 12A 0.084352 13A 0.084352 14A 0.279236 15A 0.279236 ... ... ... ... pbe0 energies and IP Num Electrons Charge E_int [a.u.] ----------------------------------------------------- N 0 -527.381683 N-1 1 -526.803729 ----------------------------------------------------- IP 0.577954 a.u. e_HOMO -0.440597 a.u. AC-shift 0.137357 a.u. ----------------------------------------------------- Psi4 stopped on: Monday, 04 January 2021 11:52AM Psi4 wall time for execution: 0:00:01.05 '''* Psi4 exiting successfully. Buy a developer a beer! }}} From here we get the following: * $\epsilon_{\rm HOMO} = -0.440597$ a.u. This is the energy of the highest-occupied molecular orbital of the Ar atom (species with $N$ electrons). Note that the output file will also contain the relevant data for the Ar$^+$ cation (species with $N-1$ electrons): we do not need any information from this. * $\text{IP} = 0.577954$ a.u. is the computed ionization potential defined as $\text{IP} = E[N-1] - E[N]$, where $E[N]$ is the energy of the species with $N$ electrons, and $E[N-1]$ is the energy of the ionised species. * In the table above, the e_HOMO and AC-shift values are incorrect. The correct e_HOMO if the one we have extracted above. The AC-shift is computed as follows:$\text{AC-Shift} = \text{IP} + \epsilon_{\rm HOMO}$ = 0.577954 - 0.440597 = 0.137357 = 0.1374 ~ $\text{a.u. to 4-sig figs.}$ * You can also use the experimental IP if you have it. You can find these at the NIST Chemistry !WebBook[[https://webbook.nist.gov/chemistry/]]