H2 ENS PDFT InversionΒΆ
[7]:
import numpy as np
import matplotlib.pyplot as plt
from CADMium import Pssolver, Psgrid, Partition, Inverter
import CADMium
a = 1.446/2
Za, Zb = 1,1
pol = 2
#Set up grid
NP = 7
NM = [4,4]
L = np.arccosh(10/a)
loc = np.array(range(-4,5)) #Stencil outline
grid = Psgrid(NP, NM, a, L, loc)
grid.initialize()
# ALPHA FRAGMENT
Nmo_a = [[1,0]]; Nmo_A = [[0,1]]
N_a = [[1,0]]; N_A = [[0,1]]
nu_a = 0.5
#Fragment b electrons
Nmo_b = [[1,0]]; Nmo_B = [[0,1]]
N_b = [[1,0]]; N_B = [[0,1]]
nu_b = 0.5
#Molecular elctron configuration
Nmo_m = [[1,1]]
N_m = [[1,1]]
part = Partition(grid, Za, Zb, pol, [Nmo_a, Nmo_A], [N_a, N_A], nu_a, [Nmo_b, Nmo_B], [N_b, N_B], nu_b, { "AB_SYM" : True,
# "ENS_SPIN_SYM" : False,
"interaction_type" : "dft",
"kinetic_part_type" : "libxcke",
"hxc_part_type" : "exact",
"k_family" : "gga",
"ke_func_id" : 500,
})
print("be careful! Treating ensembles with nu_x=1.0 will break things")
#Setup inverter object
mol_solver = Pssolver(grid, Nmo_m, N_m)
part.inverter = Inverter(grid, mol_solver, { "AB_SYM" : True,
# "ENS_SPIN_SYM" : False,
"use_iterative" : False,
"invert_type" : "wuyang",
"DISP" : False,
})
part.optPartition.isolated = True
part.scf({"disp" : True,
"alpha" : [0.6],
"e_tol" : 1e-8})
part.optPartition.isolated = False
part.scf({"disp" : True,
"alpha" : [0.6],
"max_iter" : 200,
"e_tol" : 2e-8,
"iterative" : False,
"continuing" : True})
----> Active Ensemble:
Fragment A electrons bewteen: [[1, 0]] and [[0, 1]]
Fragment B electrons between: [[1, 0]] and [[0, 1]]
be careful! Treating ensembles with nu_x=1.0 will break things
----> Begin SCF calculation for *Isolated* Fragments
Total Energy (a.u.)
__________________
Iteration A B res
_______________________________________________________
1 -0.49358 -0.49358 1.000e+00
2 -0.48337 -0.48337 7.421e-02
3 -0.48014 -0.48014 2.405e-02
4 -0.47914 -0.47914 7.514e-03
5 -0.47883 -0.47883 2.293e-03
6 -0.47874 -0.47874 6.811e-04
7 -0.47872 -0.47872 1.950e-04
8 -0.47871 -0.47871 5.238e-05
9 -0.47871 -0.47871 1.242e-05
10 -0.47871 -0.47871 4.630e-06
11 -0.47871 -0.47871 1.939e-06
12 -0.47871 -0.47871 8.214e-07
13 -0.47871 -0.47871 3.511e-07
14 -0.47871 -0.47871 1.529e-07
15 -0.47871 -0.47871 7.679e-08
16 -0.47871 -0.47871 3.658e-08
17 -0.47871 -0.47871 1.688e-08
18 -0.47871 -0.47871 7.631e-09
----> Begin SCF calculation for *Interacting* Fragments
Total Energy (a.u.)
__________________
Iteration A B res
_______________________________________________________
1 -0.37067 -0.37067 1.000e+00
2 -0.43876 -0.43876 2.064e-01
3 -0.45279 -0.45279 4.174e-02
4 -0.45536 -0.45536 7.563e-03
5 -0.45584 -0.45584 1.352e-03
6 -0.45593 -0.45593 2.493e-04
7 -0.45596 -0.45596 8.718e-05
8 -0.45596 -0.45596 3.578e-05
9 -0.45597 -0.45597 1.473e-05
10 -0.45597 -0.45597 6.109e-06
11 -0.45597 -0.45597 2.536e-06
12 -0.45597 -0.45597 1.054e-06
13 -0.45597 -0.45597 4.388e-07
14 -0.45597 -0.45597 1.833e-07
15 -0.45597 -0.45597 7.671e-08
16 -0.45597 -0.45597 3.212e-08
17 -0.45597 -0.45597 1.347e-08
[8]:
print("Separation Distance:", 2*a)
print("Fragment Energy:", part.E.Ef)
print("Partition Energy:", part.E.Ep)
print("Vnn Energy", part.E.Vnn)
print("Total Energy:", part.E.E)
Separation Distance: 1.446
Fragment Energy: -0.9119352602603984
Partition Energy: -0.9173176478006648
Vnn Energy 0.6915629322268326
Total Energy: -1.1376899758342305
[9]:
vars(part.E)
[9]:
{'Ea': -0.4559676301301992,
'Eb': -0.4559676301301992,
'Ef': -0.9119352602603984,
'Tsf': 1.2370581655198534,
'Eksf': array([[-0.37308829, -0.37308829]]),
'Enucf': -2.2015612712154513,
'Exf': -0.597456290786234,
'Ecf': -0.04649334414591891,
'Ehf': 0.6965174803673524,
'Vhxcf': 0.5432458236153052,
'Ep': -0.9173176478006648,
'Ep_pot': -1.3476199523713455,
'Ep_kin': -0.15206910662943796,
'Ep_hxc': 0.5823714112001186,
'Et': -1.829252908061063,
'Vnn': 0.6915629322268326,
'E': -1.1376899758342305,
'evals_a': array([-3.73088291e-01, -4.50359963e+15, -4.50359963e+15, -3.73088291e-01]),
'evals_b': array([-3.73088291e-01, -4.50359963e+15, -4.50359963e+15, -3.73088291e-01]),
'Ep_h': 0.5842341814267367,
'Ep_x': 0.04513385334754405,
'Ep_c': -0.046996623574162164}