Coverage for gpaw/zero_field_splitting.py: 76%
136 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1"""Zero-field splitting.
3See::
5 Spin decontamination for magnetic dipolar coupling calculations:
6 Application to high-spin molecules and solid-state spin qubits
8 Timur Biktagirov, Wolf Gero Schmidt, and Uwe Gerstmann
10 Phys. Rev. Research 2, 022024(R) – Published 30 April 2020
12"""
13from math import pi
14from typing import List, Tuple, Dict
16import numpy as np
17from ase.units import Bohr, Ha, _c, _e, _hplanck
19from gpaw.calculator import GPAW
20from gpaw.grid_descriptor import GridDescriptor
21from gpaw.typing import Array1D, Array2D, Array4D
22from gpaw.hyperfine import alpha # fine-structure constant: ~ 1 / 137
23from gpaw.setup import Setup
24from gpaw.pw.lfc import PWLFC
25from gpaw.pw.descriptor import PWDescriptor
26from gpaw.mpi import serial_comm
29def zfs(calc: GPAW,
30 method: int = 1) -> Array2D:
31 """Zero-field splitting.
33 Calculate magnetic dipole coupling tensor in eV.
34 """
35 (kpt1, kpt2), = calc.wfs.kpt_qs # spin-polarized and gamma only
37 nocc1 = (kpt1.f_n > 0.5).sum()
38 nocc2 = (kpt2.f_n > 0.5).sum()
40 assert nocc1 == nocc2 + 2, (nocc1, nocc2)
42 if method == 1:
43 wf1 = WaveFunctions.from_calc(calc, 0, nocc1 - 2, nocc1)
44 wf12 = [wf1]
45 else:
46 wf1 = WaveFunctions.from_calc(calc, 0, 0, nocc1)
47 wf2 = WaveFunctions.from_calc(calc, 1, 0, nocc2)
48 wf12 = [wf1, wf2]
50 D_vv = np.zeros((3, 3))
52 if calc.world.rank == 0:
53 compensation_charge = create_compensation_charge(wf1.setups,
54 wf1.pd,
55 calc.spos_ac)
56 for wfa in wf12:
57 for wfb in wf12:
58 d_vv = zfs1(wfa, wfb, compensation_charge)
59 D_vv += d_vv
61 calc.world.broadcast(D_vv, 0)
63 return D_vv
66class WaveFunctions:
67 def __init__(self,
68 psit_nR: Array4D,
69 P_ani: Dict[int, Array2D],
70 spin: int,
71 setups: List[Setup],
72 gd: GridDescriptor = None,
73 pd: PWDescriptor = None):
74 """Container for wave function in real-space and projections."""
76 self.pd = pd or PWDescriptor(ecut=None, gd=gd)
77 self.psit_nR = psit_nR
78 self.P_ani = P_ani
79 self.spin = spin
80 self.setups = setups
82 @staticmethod
83 def from_calc(calc: GPAW, spin: int, n1: int, n2: int) -> 'WaveFunctions':
84 """Create WaveFunctions object GPAW calculation."""
85 kpt = calc.wfs.kpt_qs[0][spin]
86 gd = calc.wfs.gd.new_descriptor(pbc_c=np.ones(3, bool),
87 comm=serial_comm)
88 psit_nR = gd.empty(n2 - n1)
89 for band, psit_R in enumerate(psit_nR):
90 psit_R[:] = calc.get_pseudo_wave_function(
91 band + n1,
92 spin=spin) * Bohr**1.5
94 return WaveFunctions(psit_nR,
95 kpt.projections.as_dict_on_master(n1, n2),
96 spin,
97 calc.setups,
98 gd=gd)
100 def __len__(self) -> int:
101 return len(self.psit_nR)
104def create_compensation_charge(setups: List[Setup],
105 pd: PWDescriptor,
106 spos_ac: Array2D) -> PWLFC:
107 compensation_charge = PWLFC([data.ghat_l for data in setups], pd)
108 compensation_charge.set_positions(spos_ac)
109 return compensation_charge
112def zfs1(wf1: WaveFunctions,
113 wf2: WaveFunctions,
114 compensation_charge: PWLFC) -> Array2D:
115 """Compute dipole coupling."""
116 pd = wf1.pd
117 setups = wf1.setups
118 N2 = len(wf2)
120 G_G = pd.G2_qG[0]**0.5
121 G_G[0] = 1.0
122 G_Gv = pd.get_reciprocal_vectors(add_q=False) / G_G[:, np.newaxis]
124 n_sG = pd.zeros(2)
125 for n_G, wf in zip(n_sG, [wf1, wf2]):
126 D_aii = {}
127 Q_aL = {}
128 for a, P_ni in wf.P_ani.items():
129 D_ii = np.einsum('ni, nj -> ij', P_ni, P_ni)
130 D_aii[a] = D_ii
131 Q_aL[a] = np.einsum('ij, ijL -> L', D_ii, setups[a].Delta_iiL)
133 for psit_R in wf.psit_nR:
134 n_G += pd.fft(psit_R**2)
136 compensation_charge.add(n_G, Q_aL)
138 nn_G = (n_sG[0] * n_sG[1].conj()).real
139 D_vv = zfs2(pd, G_Gv, nn_G)
141 n_nG = pd.empty(N2)
142 for n1, psit1_R in enumerate(wf1.psit_nR):
143 D_anii = {}
144 Q_anL = {}
145 for a, P1_ni in wf1.P_ani.items():
146 D_nii = np.einsum('i, nj -> nij', P1_ni[n1], wf2.P_ani[a])
147 D_anii[a] = D_nii
148 Q_anL[a] = np.einsum('nij, ijL -> nL',
149 D_nii, setups[a].Delta_iiL)
151 for n_G, psit2_R in zip(n_nG, wf2.psit_nR):
152 n_G[:] = pd.fft(psit1_R * psit2_R)
154 compensation_charge.add(n_nG, Q_anL)
156 nn_G = (n_nG * n_nG.conj()).sum(axis=0).real
157 D_vv -= zfs2(pd, G_Gv, nn_G)
159 D_vv -= np.trace(D_vv) / 3 * np.eye(3) # should be traceless
161 sign = 1.0 if wf1.spin == wf2.spin else -1.0
163 return sign * alpha**2 * pi * D_vv * Ha
166def zfs2(pd: PWDescriptor,
167 G_Gv: Array2D,
168 nn_G: Array1D) -> Array2D:
169 """Integral."""
170 D_vv = np.einsum('gv, gw, g -> vw', G_Gv, G_Gv, nn_G)
171 D_vv *= 2 * pd.gd.dv / pd.gd.N_c.prod()
172 return D_vv
175def convert_tensor(D_vv: Array2D,
176 unit: str = 'eV') -> Tuple[float, float, Array1D, Array2D]:
177 """Convert 3x3 tensor to D, E and easy axis.
179 Input tensor must be in eV and the result can be returned in
180 eV, μeV, MHz or 1/cm according to the value of *unit*
181 (must be one of "eV", "ueV", "MHz", "1/cm").
183 >>> D_vv = np.diag([1, 2, 3])
184 >>> D, E, axis, _ = convert_tensor(D_vv)
185 >>> D
186 4.5
187 >>> E
188 0.5
189 >>> axis
190 array([0., 0., 1.])
191 """
192 if unit == 'ueV':
193 scale = 1e6
194 elif unit == 'MHz':
195 scale = _e / _hplanck * 1e-6
196 elif unit == '1/cm':
197 scale = _e / _hplanck / _c / 100
198 elif unit == 'eV':
199 scale = 1.0
200 else:
201 raise ValueError(f'Unknown unit: {unit}')
203 (e1, e2, e3), U = np.linalg.eigh(D_vv * scale)
205 if abs(e1) > abs(e3):
206 D = 1.5 * e1
207 E = 0.5 * (e2 - e3)
208 axis = U[:, 0]
209 else:
210 D = 1.5 * e3
211 E = 0.5 * (e2 - e1)
212 axis = U[:, 2]
214 return float(D), float(E), axis, D_vv * scale
217def main(argv: List[str] = None) -> Array2D:
218 """CLI interface."""
219 import argparse
221 parser = argparse.ArgumentParser(
222 prog='python3 -m gpaw.zero_field_splitting',
223 description='...')
224 add = parser.add_argument
225 add('file', metavar='input-file',
226 help='GPW-file with wave functions.')
227 add('-u', '--unit', default='ueV', choices=['ueV', 'MHz', '1/cm'],
228 help='Unit. Must be "ueV" (micro-eV, default), "MHz" or "1/cm".')
229 add('-m', '--method', type=int, default=1)
231 args = parser.parse_intermixed_args(argv)
233 calc = GPAW(args.file)
235 D_vv = zfs(calc, args.method)
236 D, E, axis, D_vv = convert_tensor(D_vv, args.unit)
238 unit = args.unit
239 if unit == 'ueV':
240 unit = 'μeV'
242 print('D_ij = (' +
243 ',\n '.join('(' + ', '.join(f'{d:10.3f}' for d in D_v) + ')'
244 for D_v in D_vv) +
245 ') ', unit)
246 print('i, j = x, y, z')
247 print()
248 print(f'D = {D:.3f} {unit}')
249 print(f'E = {E:.3f} {unit}')
250 x, y, z = axis
251 print(f'axis = ({x:.3f}, {y:.3f}, {z:.3f})')
253 return D_vv
256if __name__ == '__main__':
257 main()