Coverage for gpaw/xc/qna.py: 97%
159 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-19 00:19 +0000
1from gpaw.xc.gga import GGA
2import numpy as np
3from gpaw.lfc import LFC
4from gpaw.spline import Spline
5from gpaw.xc.gga import gga_x, gga_c
8class QNAKernel:
9 def __init__(self, qna):
10 self.qna = qna
11 self.type = 'GGA'
12 self.name = 'QNA'
13 self.kappa = 0.804
15 def calculate(self, e_g, n_sg, dedn_sg,
16 sigma_xg, dedsigma_xg,
17 tau_sg=None, dedtau_sg=None, mu_g=None, beta_g=None,
18 dedmu_g=None, dedbeta_g=None):
20 e_g[:] = 0.
21 dedsigma_xg[:] = 0.
23 if self.qna.override_atoms is not None:
24 atoms = self.qna.override_atoms
25 self.qna.Pa.set_positions(atoms.get_scaled_positions() % 1.0)
26 else:
27 atoms = self.qna.atoms
29 if len(n_sg.shape) > 2:
30 # 3D xc calculation
31 mu_g, beta_g = self.qna.calculate_spatial_parameters(atoms)
32 dedmu_g = self.qna.dedmu_g
33 dedbeta_g = self.qna.dedbeta_g
34 else:
35 # Atomic xc calculation: use always atomwise mu and beta parameters
36 mu, beta = self.qna.parameters[atoms[self.qna.current_atom].symbol]
37 mu_g = np.zeros_like(n_sg[0])
38 beta_g = np.zeros_like(n_sg[0])
39 mu_g[:] = mu
40 beta_g[:] = beta
41 dedmu_g = None
42 dedbeta_g = None
44 # Enable to use PBE always
45 # mu_g[:] = 0.2195149727645171
46 # beta_g[:] = 0.06672455060314922
48 # Write mu and beta fields
49 if 0:
50 from ase.io import write
51 write('mu_g.cube', atoms, data=mu_g)
52 write('beta_g.cube', atoms, data=beta_g)
53 raise SystemExit
55 # spin-paired: XXX Copy-paste from gga.py to prevent
56 # distruptions to pyGGA
57 if len(n_sg) == 1:
58 n = n_sg[0]
59 n[n < 1e-20] = 1e-40
61 # exchange
62 res = gga_x(self.name, 0, n, sigma_xg[0], self.kappa, mu_g,
63 dedmu_g=dedmu_g)
64 ex, rs, dexdrs, dexda2 = res
66 if dedmu_g is not None:
67 dedmu_g[:] = n * dedmu_g
69 # correlation
70 res = gga_c(self.name, 0, n, sigma_xg[0], 0, beta_g,
71 decdbeta_g=dedbeta_g)
72 ec, rs_, decdrs, decda2, decdzeta = res
73 e_g[:] += n * (ex + ec)
74 dedn_sg[:] += ex + ec - rs * (dexdrs + decdrs) / 3.
75 dedsigma_xg[:] += n * (dexda2 + decda2)
76 # spin-polarized:
77 else:
78 na = 2. * n_sg[0]
79 na[na < 1e-20] = 1e-40
81 nb = 2. * n_sg[1]
82 nb[nb < 1e-20] = 1e-40
84 n = 0.5 * (na + nb)
85 zeta = 0.5 * (na - nb) / n
87 if dedmu_g is not None:
88 dedmua_g = dedmu_g.copy()
89 dedmub_g = dedmu_g.copy()
90 else:
91 dedmua_g = None
92 dedmub_g = None
94 # exchange
95 exa, rsa, dexadrs, dexada2 = gga_x(
96 self.name, 1, na, 4.0 * sigma_xg[0], self.kappa, mu_g,
97 dedmu_g=dedmua_g)
98 exb, rsb, dexbdrs, dexbda2 = gga_x(
99 self.name, 1, nb, 4.0 * sigma_xg[2], self.kappa, mu_g,
100 dedmu_g=dedmub_g)
101 a2 = sigma_xg[0] + 2.0 * sigma_xg[1] + sigma_xg[2]
102 if dedmu_g is not None:
103 dedmu_g[:] = 0.5 * (na * dedmua_g + nb * dedmub_g)
105 # correlation
106 ec, rs, decdrs, decda2, decdzeta = gga_c(
107 self.name, 1, n, a2, zeta, beta_g, decdbeta_g=dedbeta_g)
108 e_g[:] += 0.5 * (na * exa + nb * exb) + n * ec
109 dedn_sg[0][:] += (exa + ec - (rsa * dexadrs + rs * decdrs) / 3.0 -
110 (zeta - 1.0) * decdzeta)
111 dedn_sg[1][:] += (exb + ec - (rsb * dexbdrs + rs * decdrs) / 3.0 -
112 (zeta + 1.0) * decdzeta)
113 dedsigma_xg[0][:] += 2.0 * na * dexada2 + n * decda2
114 dedsigma_xg[1][:] += 2.0 * n * decda2
115 dedsigma_xg[2][:] += 2.0 * nb * dexbda2 + n * decda2
117 if dedbeta_g is not None:
118 dedbeta_g[:] = dedbeta_g * n
121class QNA(GGA):
122 def __init__(self, atoms, parameters, qna_setup_name='PBE', alpha=2.0,
123 override_atoms=None, stencil=2):
124 # override_atoms is only used to test the partial derivatives
125 # of xc-functional
126 kernel = QNAKernel(self)
127 GGA.__init__(self, kernel, stencil=stencil)
128 self.atoms = atoms
129 self.parameters = parameters
130 self.qna_setup_name = qna_setup_name
131 self.alpha = alpha
132 self.override_atoms = override_atoms
133 self.orbital_dependent = False
135 def todict(self):
136 dct = dict(type='qna-gga',
137 name='QNA',
138 setup_name=self.qna_setup_name,
139 parameters=self.parameters,
140 alpha=self.alpha,
141 orbital_dependent=False)
142 return dct
144 def set_grid_descriptor(self, gd):
145 GGA.set_grid_descriptor(self, gd)
146 self.dedmu_g = gd.zeros()
147 self.dedbeta_g = gd.zeros()
148 # Create gaussian LFC
149 l_lim = 1.0e-30
150 rcut = 12
151 points = 200
152 r_i = np.linspace(0, rcut, points + 1)
153 rcgauss = 1.2
154 g_g = (2 / rcgauss**3 / np.pi *
155 np.exp(-((r_i / rcgauss)**2)**self.alpha))
157 # Values too close to zero can cause numerical problems especially with
158 # forces (some parts of the mu and beta field can become negative)
159 g_g[np.where(g_g < l_lim)] = l_lim
160 spline = Spline.from_data(l=0, rmax=rcut, f_g=g_g)
161 spline_j = [[spline]] * len(self.atoms)
162 self.Pa = LFC(gd, spline_j)
164 def set_positions(self, spos_ac, atom_partition=None):
165 self.Pa.set_positions(spos_ac)
167 def calculate_spatial_parameters(self, atoms):
168 mu_g = self.gd.zeros()
169 beta_g = self.gd.zeros()
170 denominator = self.gd.zeros()
171 mu_a = {}
172 beta_a = {}
173 eye_a = {}
174 for atom in atoms:
175 mu, beta = self.parameters[atom.symbol]
176 mu_a[atom.index] = np.array([mu])
177 beta_a[atom.index] = np.array([beta])
178 eye_a[atom.index] = np.array(1.0)
179 self.Pa.add(mu_g, mu_a)
180 self.Pa.add(beta_g, beta_a)
181 self.Pa.add(denominator, eye_a)
182 mu_g /= denominator
183 beta_g /= denominator
184 return mu_g, beta_g
186 def calculate_paw_correction(self, setup, D_sp, dEdD_sp=None,
187 addcoredensity=True, a=None):
188 self.current_atom = a
189 return GGA.calculate_paw_correction(self, setup, D_sp, dEdD_sp,
190 addcoredensity, a)
192 def get_setup_name(self):
193 return self.qna_setup_name
195 def get_description(self):
196 return 'QNA Parameters: ' + str(self.parameters)
198 def add_forces(self, F_av):
199 assert self.gd.comm.size == 1, 'Domain decomposition is not supported'
200 mu_g = self.gd.zeros()
201 beta_g = self.gd.zeros()
202 denominator = self.gd.zeros()
203 mu_a = {}
204 beta_a = {}
205 eye_a = {}
206 for atom in self.atoms:
207 mu, beta = self.parameters[atom.symbol]
208 mu_a[atom.index] = np.array([mu])
209 beta_a[atom.index] = np.array([beta])
210 eye_a[atom.index] = np.array(1.0)
211 self.Pa.add(mu_g, mu_a)
212 self.Pa.add(beta_g, beta_a)
213 self.Pa.add(denominator, eye_a)
214 mu_g /= denominator
215 beta_g /= denominator
217 # mu
218 part1 = -self.dedmu_g / denominator
219 part2 = -part1 * mu_g
220 c_axiv = self.Pa.dict(derivative=True)
221 self.Pa.derivative(part1, c_axiv)
223 for atom in self.atoms:
224 F_av[atom.index] -= c_axiv[atom.index][0][:] * mu_a[atom.index][0]
225 c_axiv = self.Pa.dict(derivative=True)
226 self.Pa.derivative(part2, c_axiv)
227 for atom in self.atoms:
228 F_av[atom.index] -= c_axiv[atom.index][0][:]
230 # beta
231 part1 = -self.dedbeta_g / denominator
232 part2 = -part1 * beta_g
233 c_axiv = self.Pa.dict(derivative=True)
234 self.Pa.derivative(part1, c_axiv)
235 for atom in self.atoms:
236 F_av[atom.index] -= c_axiv[atom.index][0] * beta_a[atom.index][0]
237 c_axiv = self.Pa.dict(derivative=True)
238 self.Pa.derivative(part2, c_axiv)
239 for atom in self.atoms:
240 F_av[atom.index] -= c_axiv[atom.index][0][:]