H2+ ENS PDFT InversionΒΆ
[8]:
import numpy as np
import matplotlib.pyplot as plt
from CADMium import Pssolver, Psgrid, Partition, Inverter
import CADMium
a = 2/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 = [[1,0]]
N_a = [[0,0]]; N_A = [[1,0]]
nu_a = 0.5
#Fragment b electrons
Nmo_b = [[1,0]]; Nmo_B = [[1,0]]
N_b = [[0,0]]; N_B = [[1,0]]
nu_b = 0.5
#Molecular elctron configuration
Nmo_m = [[1,0]]
N_m = [[1,0]]
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,
})
#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: [[0, 0]] and [[1, 0]]
Fragment B electrons between: [[0, 0]] and [[1, 0]]
----> Begin SCF calculation for *Isolated* Fragments
Total Energy (a.u.)
__________________
Iteration A B res
_______________________________________________________
1 -0.24679 -0.24679 1.000e+00
2 -0.24169 -0.24169 6.139e-01
3 -0.24007 -0.24007 1.627e-01
4 -0.23957 -0.23957 4.803e-02
5 -0.23942 -0.23942 1.441e-02
6 -0.23937 -0.23937 4.263e-03
7 -0.23936 -0.23936 1.218e-03
8 -0.23935 -0.23935 3.271e-04
9 -0.23935 -0.23935 7.760e-05
10 -0.23935 -0.23935 1.314e-05
11 -0.23935 -0.23935 8.375e-07
12 -0.23935 -0.23935 2.442e-06
13 -0.23935 -0.23935 1.718e-06
14 -0.23935 -0.23935 9.537e-07
15 -0.23935 -0.23935 4.791e-07
16 -0.23935 -0.23935 2.283e-07
17 -0.23935 -0.23935 1.053e-07
18 -0.23935 -0.23935 4.763e-08
19 -0.23935 -0.23935 2.124e-08
20 -0.23935 -0.23935 9.382e-09
----> Begin SCF calculation for *Interacting* Fragments
Total Energy (a.u.)
__________________
Iteration A B res
_______________________________________________________
1 -0.17855 -0.17855 1.000e+00
2 -0.19798 -0.19798 1.815e-01
3 -0.20663 -0.20663 7.202e-02
4 -0.21012 -0.21012 2.775e-02
5 -0.21153 -0.21153 1.075e-02
6 -0.21209 -0.21209 4.201e-03
7 -0.21233 -0.21233 1.651e-03
8 -0.21242 -0.21242 6.522e-04
9 -0.21246 -0.21246 2.588e-04
10 -0.21248 -0.21248 1.032e-04
11 -0.21249 -0.21249 4.133e-05
12 -0.21249 -0.21249 1.665e-05
13 -0.21249 -0.21249 6.746e-06
14 -0.21249 -0.21249 2.752e-06
15 -0.21249 -0.21249 1.130e-06
16 -0.21249 -0.21249 4.678e-07
17 -0.21249 -0.21249 1.951e-07
18 -0.21249 -0.21249 8.207e-08
19 -0.21249 -0.21249 3.481e-08
20 -0.21249 -0.21249 1.489e-08
[7]:
print("Separation Distance:", 2*a)
print("Total Energy:", part.E.E)
print("Fragment Energy:", part.E.Ef)
print("Partition Energy:", part.E.Ep)
print("Vnn Energy", part.E.Vnn)
Separation Distance: 2.0
Total Energy: -0.5838603443185497
Fragment Energy: -0.42498501268318534
Partition Energy: -0.6588753316353645
Vnn Energy 0.5
[ ]: