Coverage for gpaw/utilities/kspot.py: 96%
112 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-19 00:19 +0000
1from math import pi, sqrt
3import numpy as np
5from ase.units import Bohr
7from gpaw.atom.all_electron import AllElectron
8from gpaw.sphere.lebedev import weight_n, R_nv
11# XXX why on earth is this necessary?
12def get_scaled_positions(atoms, positions):
13 """COPY PASTE FROM ASE! Get positions relative to unit cell.
15 Atoms outside the unit cell will be wrapped into the cell in
16 those directions with periodic boundary conditions so that the
17 scaled coordinates are between zero and one."""
19 scaled = np.linalg.solve(atoms.cell.T, positions.T).T
20 for i in range(3):
21 if atoms.pbc[i]:
22 scaled[i] %= 1.0
23 return scaled
26class AllElectronPotential:
27 def __init__(self, paw):
28 self.paw = paw
30 def write_spherical_ks_potentials(self, txt):
31 f = open(txt, 'w')
32 for a in self.paw.density.D_asp:
33 r_g, vKS_g = self.get_spherical_ks_potential(a)
34 setup = self.paw.density.setups[a]
35 # Calculate also atomic LDA for reference
36 g = AllElectron(setup.symbol,
37 xcname='LDA',
38 nofiles=True,
39 scalarrel=True,
40 txt=None)
41 g.run()
42 print(r_g[0], vKS_g[0], g.vr[0], 0.0, file=f)
43 for r, vKS, vr in zip(r_g[1:], vKS_g[1:], g.vr[1:]):
44 print(r, vKS, vr, (vKS - vr) / r, file=f)
46 f.close()
48 def grid_to_radial(self, a, gd, f_g):
50 # Coordinates of an atom
51 atom_c = self.paw.atoms.get_positions()[a]
53 # Get xccorr for atom a
54 setup = self.paw.density.setups[a]
55 xccorr = setup.xc_correction
57 radf_g = np.zeros(xccorr.ng)
59 for w, p in zip(weight_n, R_nv):
60 scaled_nc = []
61 # Very inefficient loop
62 for i, r in enumerate(xccorr.rgd.r_g):
63 # Obtain the position of this integration quadrature point in
64 # specified grid
65 pos_c = atom_c + (r * Bohr) * p
66 # And in scaled coordinates
67 scaled_c = get_scaled_positions(self.paw.atoms, pos_c)
68 scaled_nc.append(scaled_c)
70 scaled_nc = np.array(scaled_nc)
72 target_g = gd.interpolate_grid_points(scaled_nc, f_g)
73 radf_g += w * target_g
75 return radf_g
77 def get_spherical_ks_potential(self, a):
78 # self.paw.restore_state()
80 print("XC:", self.paw.hamiltonian.xc.name)
81 assert self.paw.hamiltonian.xc.type == 'LDA'
83 # If the calculation is just loaded, density needs to be interpolated
84 if self.paw.density.nt_sg is None:
85 print("Interpolating density")
86 self.paw.density.interpolate_pseudo_density()
88 # Get xccorr for atom a
89 setup = self.paw.density.setups[a]
90 xccorr = setup.xc_correction
92 # Get D_sp for atom a
93 D_sp = self.paw.density.D_asp[a]
95 # density a function of L and partial wave radial pair density
96 # coefficient
97 D_sLq = np.inner(D_sp, xccorr.B_pqL.T)
99 # The 'spherical' spherical harmonic
100 Y0 = 1.0 / sqrt(4 * pi)
102 # Generate cartesian fine grid xc-potential
103 print("Generate cartesian fine grid xc-potential")
104 gd = self.paw.density.finegd
105 vxct_sg = gd.zeros(1)
106 xc = self.paw.hamiltonian.xc
107 xc.calculate(gd, self.paw.density.nt_sg, vxct_sg)
109 # ---------------------------------------------
110 # The Electrostatic potential
111 # ---------------------------------------------
112 # V_ES(r) = Vt_ES(r) - Vt^a(r) + V^a(r), where
113 # Vt^a = P[ nt^a(r) + \sum Q_L g_L(r) ]
114 # V^a = P[ -Z\delta(r) + n^a(r) ]
115 # P = Poisson solution
116 # ---------------------------------------------
118 print("Evaluating ES Potential...")
119 # Make sure that the calculation has ES potential
120 # TODO
121 if self.paw.hamiltonian.vHt_g is None:
122 raise "Converge tha Hartree potential first."
124 # Interpolate the smooth electrostatic potential from fine grid to
125 # radial grid
126 radHt_g = self.grid_to_radial(a, gd, self.paw.hamiltonian.vHt_g)
127 radHt_g *= xccorr.rgd.r_g
129 print("D_sp", D_sp)
131 # Calculate the difference in density and pseudo density
132 dn_g = (Y0 * np.dot(D_sLq[0, 0], (xccorr.n_qg - xccorr.nt_qg)) +
133 xccorr.nc_g - xccorr.nct_g)
135 # Add the compensation charge contribution
136 dn_g -= Y0 * self.paw.density.Q_aL[a][0] * setup.g_lg[0]
138 # Calculate the Hartree potential for this
139 vHr = xccorr.rgd.poisson(dn_g)
141 # Add the core potential contribution
142 vHr -= setup.Z
144 radHt_g += vHr
146 # --------------------------------------------
147 # The XC potential
148 # --------------------------------------------
149 # V_xc = Vt_xc(r) - Vt_xc^a(r) + V_xc^a(r)
150 # --------------------------------------------
152 print("Evaluating xc potential")
153 # Interpolate the smooth xc potential from fine grid to radial grid
154 radvxct_g = self.grid_to_radial(a, gd, vxct_sg[0])
156 # Arrays for evaluating radial xc potential slice
157 vxc_sg = np.zeros((len(D_sp), xccorr.ng))
159 for n, Y_L in enumerate(xccorr.Y_nL):
160 n_sLg = np.dot(D_sLq, xccorr.n_qg)
161 n_sLg[:, 0] += sqrt(4 * pi) * xccorr.nc_g
162 vxc_sg[:] = xc.calculate_radial(xccorr.rgd, n_sLg, Y_L)[1]
163 radvxct_g += weight_n[n] * vxc_sg[0]
164 nt_sLg = np.dot(D_sLq, xccorr.nt_qg)
165 nt_sLg[:, 0] += sqrt(4 * pi) * xccorr.nct_g
166 vxc_sg[:] = xc.calculate_radial(xccorr.rgd, nt_sLg, Y_L)[1]
167 radvxct_g -= weight_n[n] * vxc_sg[0]
169 radvks_g = radvxct_g * xccorr.rgd.r_g + radHt_g
170 return (xccorr.rgd.r_g, radvks_g)
173class CoreEigenvalues(AllElectronPotential):
174 def get_core_eigenvalues(self, a, scalarrel=True):
175 """Return the core eigenvalues by solving the radial schrodinger
176 equation.
178 Using AllElectron potential class, the spherically averaged Kohn--Sham
179 potential is obtained around the spesified atom. The eigenstates for
180 this potential are solved, the the resulting core states returned.
181 Still experimental.
183 """
185 r, v_g = self.get_spherical_ks_potential(a)
187 # Get xccorr for atom a
188 setup = self.paw.density.setups[a]
189 symbol = setup.symbol
191 # Create AllElectron object for eigensolver
192 atom = AllElectron(symbol, txt=None, scalarrel=scalarrel)
193 # Calculate initial guess
194 atom.run()
196 # Set the potential
197 atom.vr[:len(v_g)] = v_g
198 # After the setups cutoff, arbitrary barrier is used
199 atom.vr[len(v_g):] = 10.0
201 # Solve the eigenstates
202 atom.solve()
204 # The display format is just copy paste from AllElectron class
205 # TODO: Make it a method in AllElectron class, thus it can be called
206 # directly
207 def t(a):
208 print(a)
210 t('Calculated core eigenvalues of atom ' + str(a) + ':' + symbol)
211 t('state eigenvalue ekin rmax')
212 t('-----------------------------------------------')
213 for m, l, f, e, u in zip(atom.n_j, atom.l_j, atom.f_j, atom.e_j,
214 atom.u_j):
215 # Find kinetic energy:
216 k = e - np.sum((
217 np.where(abs(u) < 1e-160, 0, u)**2 * # XXXNumeric!
218 atom.vr * atom.dr)[1:] / atom.r[1:])
220 # Find outermost maximum:
221 g = atom.N - 4
222 while u[g - 1] >= u[g]:
223 g -= 1
224 x = atom.r[g - 1:g + 2]
225 y = u[g - 1:g + 2]
226 A = np.transpose(np.array([x**i for i in range(3)]))
227 c, b, a = np.linalg.solve(A, y)
228 assert a < 0.0
229 rmax = -0.5 * b / a
231 s = 'spdf' [l]
232 t('%d%s^%-4.1f: %12.6f %12.6f %12.3f' % (m, s, f, e, k, rmax))
233 t('-----------------------------------------------')
234 t('(units: Bohr and Hartree)')
235 return atom.e_j