Line data Source code
1 : !!****m* ABINIT/m_classify_bands
2 : !! NAME
3 : !! m_classify_bands
4 : !!
5 : !! FUNCTION
6 : !! Finds the irreducible representation associated to
7 : !! a set of degenerate bands at a given k-point and spin.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2008-2026 ABINIT group (MG)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_classify_bands
24 :
25 : use, intrinsic :: iso_c_binding, only : c_f_pointer, c_loc
26 : use defs_basis
27 : use m_abicore
28 : use m_xmpi
29 : use m_errors
30 : use netcdf
31 : use m_nctk
32 :
33 : use m_fstrings, only : itoa, ftoa, sjoin, ktoa, ltoa, strcat, yesno
34 : use m_geometry, only : getspinrot
35 : use defs_datatypes, only : pseudopotential_type
36 : use m_dtset, only : dataset_type
37 : use m_dtfil, only : datafiles_type
38 : use m_io_tools, only : iomode_from_fname
39 : use m_time, only : cwtime, cwtime_report
40 : use m_numeric_tools, only : get_trace, print_arr
41 : use m_matrix, only : is_unitary, is_identity, mati3inv
42 : use m_hdr, only : hdr_type
43 : use m_hide_blas, only : xdotc, xdotu, xcopy
44 : use m_fft_mesh, only : rotate_FFT_mesh, calc_ceigr
45 : use m_crystal, only : crystal_t
46 : use m_cgtools, only : cg_zdotc
47 : use m_symtk, only : sg_multable, sym_order
48 : use m_pawang, only : pawang_type
49 : use m_pawrad, only : pawrad_type
50 : use m_pawtab, only : pawtab_type, pawtab_get_lsize
51 : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_print, pawfgrtab_free
52 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy
53 : use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
54 : use m_paw_sphharm, only : setsym_ylm
55 : use m_paw_nhat, only : nhatgrid
56 : use m_wfd, only : wfd_t
57 : use m_ebands, only : ebands_t
58 : use m_common, only : ebands_from_file
59 : use m_fftcore, only : sphere, get_kg, ngfft_seq
60 : use m_cgtk, only : cgtk_rotate, cgtk_change_gsphere
61 : use m_kpts, only : kpts_map
62 : use m_krank, only : krank_t
63 : use m_esymm, only : esymm_t, esymm_free
64 : use m_ptgroups, only : get_classes
65 : use m_yaml, only : yamldoc_t, yamldoc_open
66 : use m_pair_list, only : pair_list
67 :
68 : implicit none
69 :
70 : private
71 : !!***
72 :
73 : public :: classify_bands
74 :
75 : !!****t* m_classify_bands/dmats_t
76 : !! NAME
77 : !! dmats_t
78 : !!
79 : !! FUNCTION
80 : !! Store D_mn(S) = <psi_{mSk}| S | psi_{nk}> for all the k-points in the IBZ
81 : !! and the bands in brange_spin.
82 : !!
83 : !! SOURCE
84 :
85 : type, public :: dmats_t
86 :
87 : type(ebands_t) :: ks_ebands
88 : ! KS bands.
89 :
90 : type(crystal_t),pointer :: cryst => null()
91 : type(dataset_type),pointer :: dtset => null()
92 :
93 : integer,allocatable :: brange_spin(:,:)
94 : ! (2, nsppol)
95 : ! start and end band index for each spin
96 :
97 : integer,allocatable :: multable(:,:,:)
98 : ! (4,nsym,nsym)
99 : ! multable(1,sym1,sym2) gives the index of the symmetry product S1 * S2 in the symrel array. 0 if not found.
100 : ! multable(2:4,sym1,sym2)= the lattice vector that has to added to the fractional translation
101 : ! of the operation of index multable(1,sym1,sym2) to obtain the fractional translation of the product S1 * S2.
102 :
103 : integer,allocatable :: toinv(:,:)
104 : ! (4,nsym)
105 : ! toinv(1,sym1)=Gives the index of the inverse of the symmetry operation.
106 : ! S1 * S1^{-1} = {E, L} with E the identity and L a real-space lattice vector.
107 : ! toinv(2:4,sym1)=The lattice vector L
108 : ! Note that toinv can be easily obtained from multable but sometimes we do not need the full table.
109 :
110 : type(coeff5c_type), allocatable :: for_spin(:)
111 :
112 : contains
113 : procedure :: init => dmats_init ! Initialize object
114 : procedure :: init_from_file => dmats_init_from_file ! Initialize object from a NetCDF file
115 : procedure :: ncwrite => dmats_ncwrite ! Write D-matrices to an open NetCDF file
116 : procedure :: free => dmats_free ! Free memory.
117 : procedure :: check => dmats_check ! Check Dmats
118 : procedure :: classify => dmats_classify ! Classify irreps
119 : procedure :: get_star_dmats => dmats_get_star_dmats ! D-matrices at k'=S0.k_ibz
120 : procedure :: get_star_dmats_at_kpt => dmats_get_star_dmats_at_kpt ! Same, locating S0 from a raw kpt
121 : procedure :: check_star => dmats_check_star ! Run dmats_check_one_k at a star kpt
122 : end type dmats_t
123 : !!***
124 :
125 : contains
126 : !!***
127 :
128 : !!****f* m_classify_bands/classify_bands
129 : !! NAME
130 : !! classify_bands
131 : !!
132 : !! FUNCTION
133 : !! This routine finds the irreducible representation associated to
134 : !! a set of degenerate bands at a given k-point and spin.
135 : !! The irreducible representation is obtained by rotating the set
136 : !! of degenerate wavefunctions using the symmetry operations in the little group of k.
137 : !! Two states are treated as degenerate if their energy differs by less than EDIFF_TOL.
138 : !!
139 : !! INPUTS
140 : !! Wfd(wfd_t)= structure gathering information on wave functions
141 : !! ik_ibz=The index of the k-point in the IBZ.
142 : !! spin=The spin index.
143 : !! ngfft(18)=Info on the FFT mesh to be used for evaluting u(r) and the rotated u(R^{1}(r-t)).
144 : !! ngfft must be compatible with the symmetries of the crystal and can differ from Wfd%ngfft.
145 : !! wfd_change_ngfft is called if ANY(Wfd%ngfft(1:3) =/ ngfft).
146 : !! Cryst<crystal_t>=Type gathering info on the crystal structure.
147 : !! ebands<ebands_t>=Datatype with electronic energies.
148 : !! Pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
149 : !! Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data.
150 : !! Pawang <type(pawang_type)>=paw angular mesh and related data
151 : !! Psps<pseudopotential_type>
152 : !! %indlmn(6,lmnmax,ntypat)=array giving l,m,n,lm,ln,spin for i=lmn (for each atom type)
153 : !! Dtfil<datafiles_type>=variables related to files
154 : !! %unpaw
155 : !! tolsym=Tolerance for the symmetries (input variable)
156 : !! [EDIFF_TOL]= tolerance on the energy difference of two states (if not specified is set to 0.005 eV)
157 : !!
158 : !! OUTPUT
159 : !! BSym<Bands_Symmetries>=structure containing info on the little group of the k-point as well
160 : !! as the character of the representation associated to each set of degenerate states
161 : !! if BSym%isymmorphic the symmetry analysis cannot be performed, usually it means that
162 : !! k is at zone border and there are non-symmorphic translations (see Notes)
163 : !!
164 : !! NOTES
165 : !! * Let M(R_t) the irreducible representation associated to the space group symmetry (R_t).
166 : !! * By convention M(R_t) multiplies wave functions as a row vector:
167 : !!
168 : !! $ R_t \psi_a(r) = \psi_a (R^{-1}(r-\tau)) = \sum_b M(R_t)_{ba} \psi_b $
169 : !!
170 : !! Therefore, if R_t belongs to the little group of k (i.e. Sk=k+G0), one obtains:
171 : !!
172 : !! $ M_ab(R_t) = e^{-i(k+G0).\tau} \int e^{iG0.r} u_{ak}(r)^* u_{bk}(R^{-1}(r-\tau)) \,dr $.
173 : !!
174 : !! * The irreducible representation of the small _point_ group of k, M_ab(R), suffices to
175 : !! classify the degenerate eigenstates provided that particular conditions are fulfilled
176 : !! (see limitations below). The matrix is indeed given by:
177 : !!
178 : !! $ M_ab(R) = e^{+ik.\tau} M_ab(R_t) = e^{-iG0.\tau} \int e^{iG0.r} u_{ak}(r)^* u_{bk}(R^{-1}(r-\tau))\,dr $
179 : !!
180 : !! The phase factor outside the integral should be zero since symmetry analysis at border zone in non-symmorphic
181 : !! space groups is not available. Anyway it is included in our expressions for the sake of consistency.
182 : !!
183 : !! * For PAW there is an additional onsite terms involving <phi_i|phi_j(R^{-1}(r-\tau)> and
184 : !! the pseudized version that can be evaluated using the rotation matrix for
185 : !! real spherical harmonis, zarot(mp,m,l,R). $ Y_{lm}(Rr)= \sum_{m'} zarot(m',m,ll,R) Y_{lm'}(r) $
186 : !!
187 : !! $ M^{onsite}_ab(R_t) = sum_{c ij} <\tpsi_a| p_i^c> <p_j^{c'}|\tpsi_b\> \times
188 : !! [ <\phi_i^c|\phi_j^{c'}> - <\tphi_i^c|\tphi_j^{c'}> ]. $
189 : !!
190 : !! $ [ <\phi_i^c|\phi_j^{c'}> - <\tphi_i^c|\tphi_j^{c'}> ] = s_{ij} D_{\mi\mj}^\lj(R^{-1}) $
191 : !!
192 : !! where c' is the rotated atom i.e c' = R^{-1}( c-\tau) and D is the rotation matrix for
193 : !! real spherical harmonics.
194 : !!
195 : !! Remember that zarot(m',m,l,R)=zarot(m,m',l,R^{-1})
196 : !! and $ Y^l_m(ISG) = sum_{m'} D_{m'm}(S) Y_{m'}^l(G) (-i)^l $
197 : !! $ D_{m'm}^l (R) = D_{m,m'}^l (R^{-1}) $
198 : !!
199 : !! * LIMITATIONS: The method does not work if k is at zone border and the little group of k
200 : !! contains a non-symmorphic fractional translation.
201 : !!
202 : !! SOURCE
203 :
204 0 : subroutine classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,&
205 0 : Cryst,ebands,Pawtab,Pawrad,Pawang,Psps,tolsym,BSym,&
206 : EDIFF_TOL) ! optional
207 :
208 : !Arguments ------------------------------------
209 : !scalars
210 : integer,intent(in) :: ik_ibz,spin,first_band,last_band
211 : real(dp),intent(in) :: tolsym
212 : real(dp),intent(in),optional :: EDIFF_TOL
213 : logical,intent(in) :: use_paw_aeur
214 : type(crystal_t),intent(in) :: Cryst
215 : type(pawang_type),intent(in) :: Pawang
216 : type(pseudopotential_type),intent(in) :: Psps
217 : class(wfd_t),intent(inout) :: Wfd
218 : type(ebands_t),target,intent(in) :: ebands
219 : type(esymm_t),intent(out) :: BSym
220 : !arrays
221 : integer,intent(in) :: ngfftf(18)
222 : type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
223 : type(Pawrad_type),intent(inout) :: Pawrad(Cryst%ntypat*Wfd%usepaw)
224 :
225 : !Local variables-------------------------------
226 : !scalars
227 : integer,parameter :: nspinor1=1
228 : integer :: dim_degs,ib1,ib2,ib_stop,ib_start,iclass,idg,sym_idx
229 : integer :: ir,isym,isym_class,tr_isym,jb1,jb2
230 : integer :: nr1,nr2,nr3,nsym_class,nfft,cplex !,ifgd,nfgd,ifft_sph
231 : integer :: ii,jj,lmax
232 : integer :: optcut,optgr0,optgr1,optgr2,optrad
233 : real(dp) :: EDIFF_TOL_,arg,fft_fact
234 : complex(dp) :: exp_mikg0t,exp_ikg0t,cmat_ab
235 : logical :: iscompatibleFFT,found,only_trace
236 : character(len=500) :: msg
237 : !arrays
238 0 : integer :: g0(3), toinv(Cryst%nsym), trial(3,3)
239 0 : integer,pointer :: Rm1_rmt(:)
240 0 : integer,target,allocatable :: irottb(:,:)
241 0 : integer,allocatable :: tmp_sym(:,:,:),l_size_atm(:)
242 : real(dp) :: kpt(3),kpg0(3),omat(2)
243 0 : real(dp),pointer :: ene_k(:), zarot(:,:,:,:)
244 0 : complex(dp),allocatable :: eig0r(:,:),tr_emig0r(:,:)
245 0 : complex(gwp),allocatable :: ur1(:),ur2(:),ur2_rot(:)
246 0 : type(pawcprj_type),allocatable :: Cprj_b1(:,:),Cprj_b2(:,:),Cprj_b2rot(:,:)
247 0 : type(Pawfgrtab_type),allocatable :: Pawfgrtab(:)
248 0 : type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
249 : ! *************************************************************************
250 :
251 : ! Consistency check on input.
252 0 : ABI_CHECK(Wfd%nspinor == 1, 'nspinor/=1 not coded')
253 :
254 : ! By default all bands are included
255 : !first_band=1; last_band=Wfd%nband(ik_ibz,spin)
256 0 : ABI_CHECK(first_band == 1, "first_band/=1 not coded")
257 0 : ABI_CHECK(last_band <= Wfd%nband(ik_ibz,spin), "last_band cannot be > nband_k")
258 :
259 0 : EDIFF_TOL_= 0.005/Ha_eV; if (PRESENT(EDIFF_TOL)) EDIFF_TOL_=ABS(EDIFF_TOL)
260 :
261 0 : call wfd%change_ngfft(Cryst,Psps,ngfftf)
262 :
263 : ! Get index of the rotated FFT points ===
264 : ! FFT mesh in real space _must_ be compatible with symmetries.
265 0 : nr1 = Wfd%ngfft(1)
266 0 : nr2 = Wfd%ngfft(2)
267 0 : nr3 = Wfd%ngfft(3)
268 0 : nfft = Wfd%nfft ! No FFT parallelism
269 :
270 0 : ABI_MALLOC(irottb,(nfft,Cryst%nsym))
271 0 : call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,irottb,iscompatibleFFT)
272 :
273 0 : if (.not.iscompatibleFFT) then
274 : write(msg,'(3a)')&
275 0 : ' For symmetry analysis, the real space FFT mesh must be compatible with the symmetries of the space group',ch10,&
276 0 : ' classify_bands will return. Action: change the input variable ngfftf '
277 0 : ABI_WARNING(msg)
278 0 : Bsym%err_status=1
279 0 : Bsym%err_msg= msg
280 0 : RETURN
281 : end if
282 :
283 : ! only_trace=if .TRUE. only the trace of a single matrix per class is calculated (standard procedure if
284 : ! only the symmetry of bands is required). If .FALSE. all the matrices for each irreducible representation
285 : ! are calculated and stored in BSym
286 0 : only_trace=.FALSE.
287 : !
288 : ! ==========================================
289 : ! ==== Analyse k-point symmetries first ====
290 : ! ==========================================
291 : ! The analysis is done here so that we already know if there is a problem.
292 0 : kpt = Wfd%kibz(:,ik_ibz)
293 : !
294 : !----Initialize the Bsym structure for this k-point and spin----!
295 : ! NOTE that all the degenerate states should be included! No check is done.
296 :
297 0 : ene_k => ebands%eig(first_band:,ik_ibz,spin) ! Select a slice of eigenvalues
298 :
299 0 : call Bsym%init(kpt, Cryst, only_trace, Wfd%nspinor, first_band, last_band, EDIFF_TOL_, ene_k, tolsym)
300 : !Bsym%degs_bounds = Bsym%degs_bounds + (first_band -1)
301 :
302 0 : if (Bsym%err_status /= 0) then
303 0 : write(msg,'(a,i0,a)')" esymm_init returned err_status= ",Bsym%err_status," Band classifications cannot be performed."
304 0 : ABI_WARNING(msg)
305 0 : RETURN
306 : end if
307 :
308 0 : do ii=1,Cryst%nsym
309 0 : call mati3inv(Cryst%symrel(:,:,ii),trial)
310 0 : trial=transpose(trial)
311 0 : found=.FALSE.
312 0 : do jj=1,Cryst%nsym
313 0 : if (ALL(trial==Cryst%symrel(:,:,jj))) then
314 0 : toinv(ii)=jj
315 : !toinv(jj)=ii
316 : found=.TRUE.; EXIT
317 : end if
318 : end do
319 0 : if (.not.found) then
320 0 : ABI_ERROR("inverse not found! ")
321 : end if
322 : end do
323 :
324 0 : nullify(zarot)
325 :
326 0 : if (Wfd%usepaw==1) then ! Allocate cprj_k and cprj_krot to store a set of bands for a single (K,SPIN).
327 0 : ABI_MALLOC(Cprj_b1 ,(Cryst%natom,Wfd%nspinor))
328 0 : call pawcprj_alloc(Cprj_b1, 0,Wfd%nlmn_atm)
329 0 : ABI_MALLOC(Cprj_b2 ,(Cryst%natom,Wfd%nspinor))
330 0 : call pawcprj_alloc(Cprj_b2, 0,Wfd%nlmn_atm)
331 0 : ABI_MALLOC(Cprj_b2rot,(Cryst%natom,Wfd%nspinor))
332 0 : call pawcprj_alloc(Cprj_b2rot,0,Wfd%nlmn_atm)
333 :
334 : !zarot => Pawang%zarot
335 0 : lmax = Pawang%l_max-1
336 0 : ABI_MALLOC(zarot,(2*lmax+1,2*lmax+1,lmax+1,Cryst%nsym))
337 0 : zarot = Pawang%zarot
338 :
339 0 : ABI_MALLOC(tmp_sym,(3,3,Cryst%nsym))
340 0 : do isym=1,Cryst%nsym
341 0 : tmp_sym(:,:,isym) = Cryst%symrel(:,:,isym)
342 : !tmp_sym(:,:,isym) = Cryst%symrel(:,:,toinv(isym))
343 : !tmp_sym(:,:,isym) = transpose(Cryst%symrel(:,:,isym))
344 : !tmp_sym(:,:,isym) = Cryst%symrec(:,:,isym)
345 : !tmp_sym(:,:,isym) = TRANSPOSE(Cryst%symrec(:,:,isym))
346 : end do
347 : !% call setsym_ylm(Cryst%rprimd,lmax,Cryst%nsym,3,Cryst%gprimd,tmp_sym,zarot)
348 : !call setsym_ylm(Cryst%gprimd,lmax,Cryst%nsym,1,Cryst%rprimd,tmp_sym,zarot)
349 0 : ABI_FREE(tmp_sym)
350 0 : zarot = Pawang%zarot
351 :
352 0 : cplex=1
353 0 : call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat)
354 0 : ABI_MALLOC(Pawfgrtab,(Cryst%natom))
355 0 : call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Wfd%nspden,Cryst%typat)
356 0 : ABI_FREE(l_size_atm)
357 :
358 0 : optcut=1 ! use rpaw to construct local_pawfgrtab
359 0 : optgr0=0; optgr1=0; optgr2=0 ! dont need gY terms locally
360 0 : optrad=1 ! do store r-R
361 :
362 : call nhatgrid(Cryst%atindx1,Cryst%gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,Wfd%ngfft,Cryst%ntypat,&
363 0 : optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
364 :
365 : !call pawfgrtab_print(Pawfgrtab,unit=std_out,Wfd%prtvol=10)
366 :
367 0 : ABI_MALLOC(Paw_onsite,(Cryst%natom))
368 :
369 0 : if (use_paw_aeur) then
370 0 : ABI_WARNING("Using AE wavefunction for rotation in real space!")
371 : call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,&
372 0 : Cryst%rprimd,Cryst%xcart,Pawtab,Pawrad,Pawfgrtab)
373 : end if
374 : end if
375 :
376 : ! ===============================================
377 : ! ==== Calculate the representation matrices ====
378 : ! ===============================================
379 0 : fft_fact=one/nfft
380 0 : ABI_MALLOC(ur1, (nfft))
381 0 : ABI_MALLOC(ur2, (nfft))
382 0 : ABI_MALLOC(ur2_rot, (nfft))
383 :
384 : ! Precalculate eig0r = e^{iG0.r} on the FFT mesh.
385 0 : ABI_MALLOC(eig0r, (nfft, Bsym%nsym_gk))
386 :
387 0 : do isym=1,Bsym%nsym_gk
388 0 : g0 = Bsym%g0(:,isym)
389 0 : call calc_ceigr(g0,nfft,nspinor1,Wfd%ngfft,eig0r(:,isym))
390 : end do
391 :
392 0 : if (Bsym%can_use_tr) then
393 0 : ABI_MALLOC(tr_emig0r,(nfft,Bsym%nsym_trgk))
394 0 : do isym=1,Bsym%nsym_trgk
395 0 : g0=Bsym%tr_g0(:,isym)
396 0 : call calc_ceigr(-g0,nfft,nspinor1,Wfd%ngfft,tr_emig0r(:,isym))
397 : end do
398 : end if
399 :
400 : ! Loop over the set of degenerate states.
401 0 : do idg=1,Bsym%ndegs
402 0 : ib_start = Bsym%degs_bounds(1,idg)
403 0 : ib_stop = Bsym%degs_bounds(2,idg)
404 0 : dim_degs = Bsym%degs_dim(idg)
405 :
406 0 : do ib1=ib_start,ib_stop ! First band index in the degenerate set.
407 0 : jb1=ib1-ib_start+1
408 :
409 : ! debugging: use AE wave on dense FFT mesh.
410 0 : if (Wfd%usepaw==1..and.use_paw_aeur) then
411 0 : call wfd%paw_get_aeur(ib1,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur1)
412 : else
413 0 : call wfd%get_ur(ib1,ik_ibz,spin,ur1)
414 0 : if (Wfd%usepaw==1) then
415 0 : call wfd%ug2cprj(ib1,ik_ibz,spin,1,0,Cryst%natom,Cryst,Cprj_b1,sorted=.FALSE.)
416 : end if
417 : end if
418 :
419 0 : do ib2=ib_start,ib_stop ! Second band index in the degenerate set.
420 0 : if (Bsym%only_trace.and.ib1/=ib2) CYCLE ! Only the diagonal is needed.
421 :
422 0 : if (ib2==ib1) then
423 0 : call xcopy(nfft,ur1,1,ur2,1)
424 0 : if (Wfd%usepaw==1) call pawcprj_copy(Cprj_b1,Cprj_b2)
425 : else
426 : ! debugging: use AE wave on dense FFT mesh.
427 0 : if (Wfd%usepaw==1.and.use_paw_aeur) then
428 0 : call wfd%paw_get_aeur(ib2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur2)
429 : else
430 0 : call wfd%get_ur(ib2,ik_ibz,spin,ur2)
431 0 : if (Wfd%usepaw==1) then
432 0 : call wfd%ug2cprj(ib2,ik_ibz,spin,1,0,Cryst%natom,Cryst,Cprj_b2,sorted=.FALSE.)
433 : end if
434 : end if
435 : end if
436 :
437 : ! ===================================================
438 : ! ==== Loop over the classes of the little group ====
439 : ! ===================================================
440 0 : sym_idx=0
441 0 : do iclass=1,Bsym%nclass
442 0 : nsym_class = Bsym%nelements(iclass)
443 :
444 0 : do isym_class=1,nsym_class ! Loop over elements in each class.
445 0 : sym_idx = sym_idx+1
446 0 : if (Bsym%only_trace.and.isym_class/=1) CYCLE ! Do it once if only the character is required.
447 :
448 0 : isym = Bsym%sgk2symrec(sym_idx)
449 0 : Rm1_rmt => irottb(:,isym)
450 :
451 : ! Classify states according to the irreps of the little group of k.
452 0 : kpg0= kpt + Bsym%g0(:,sym_idx)
453 0 : arg=-two_pi * DOT_PRODUCT(kpg0,Cryst%tnons(:,isym))
454 :
455 0 : if (ABS(arg) > tol6) then
456 0 : exp_mikg0t = DCMPLX(DCOS(arg),DSIN(arg))
457 : else
458 : exp_mikg0t = cone
459 : end if
460 :
461 : !if (Wfd%usepaw==1) then
462 : !end if
463 : !
464 : ! Rotate the right wave function and apply the phase ===
465 : ! Note that the k-point is the same within a lattice vector.
466 0 : do ir=1,nfft
467 0 : ur2_rot(ir)=ur2(Rm1_rmt(ir))*eig0r(ir,sym_idx)
468 : end do
469 :
470 : ! The matrix element on the FFT mesh.
471 0 : cmat_ab = xdotc(nfft,ur1,1,ur2_rot,1)*fft_fact*exp_mikg0t
472 :
473 0 : if (Wfd%usepaw==1.and..not.use_paw_aeur) then ! Add the on-site contribution.
474 0 : call rotate_cprj(kpt,isym,Wfd%nspinor,1,Cryst%natom,Cryst%nsym,Cryst%typat,Cryst%indsym,Cprj_b2,Cprj_b2rot)
475 :
476 : omat = paw_phirotphj(Wfd%nspinor,Cryst%natom,Cryst%typat,&
477 0 : zarot(:,:,:,isym),Pawtab,Psps,Cprj_b1,Cprj_b2rot)
478 :
479 0 : cmat_ab = cmat_ab + DCMPLX(omat(1),omat(2)) !* exp_mikg0t
480 : end if
481 :
482 0 : jb2 = ib2 - ib_start+1
483 0 : Bsym%Calc_irreps(idg)%mat(jb1,jb2,sym_idx)=cmat_ab
484 :
485 : end do !isym_class
486 : end do !iclass
487 : !
488 : ! =========================================================
489 : ! ==== Loop over the symmetries such that -Sk = k + G0 ====
490 : ! =========================================================
491 : ! <-k,a| S |k b> = e^{i(k+G0).t} \int e^{-ig0.r} u_a u_b(R^{1}(r-t))
492 0 : if (Bsym%can_use_tr) then
493 0 : do tr_isym=1,Bsym%nsym_trgk
494 :
495 0 : isym = Bsym%tr_sgk2symrec(tr_isym)
496 0 : Rm1_rmt => irottb(:,isym)
497 :
498 0 : kpg0= kpt + Bsym%tr_g0(:,tr_isym)
499 0 : arg= two_pi * DOT_PRODUCT(kpg0,Cryst%tnons(:,isym))
500 :
501 0 : if (ABS(arg) > tol6) then
502 0 : exp_ikg0t=DCMPLX(DCOS(arg),DSIN(arg))
503 : else
504 : exp_ikg0t=cone
505 : end if
506 :
507 : ! Rotate the right wave function and apply the phase
508 : ! Note that the k-point is the same within a lattice vector.
509 0 : do ir=1,nfft
510 0 : ur2_rot(ir)=ur2(Rm1_rmt(ir)) * tr_emig0r(ir,tr_isym)
511 : end do
512 :
513 : ! The matrix element on the FFT mesh.
514 0 : cmat_ab = xdotu(nfft,ur1,1,ur2_rot,1)*fft_fact*exp_ikg0t
515 :
516 0 : if (Wfd%usepaw==1.and..not.use_paw_aeur) then ! Add the on-site contribution. ! TODO rechek this part.
517 0 : call rotate_cprj(kpt,isym,Wfd%nspinor,1,Cryst%natom,Cryst%nsym,Cryst%typat,Cryst%indsym,Cprj_b2,Cprj_b2rot)
518 : omat = paw_phirotphj(Wfd%nspinor,Cryst%natom,Cryst%typat,&
519 0 : zarot(:,:,:,isym),Pawtab,Psps,Cprj_b1,Cprj_b2rot,conjg_left=.TRUE.)
520 0 : cmat_ab = cmat_ab + DCMPLX(omat(1),omat(2)) !* exp_ikg0t
521 : end if
522 :
523 0 : jb2 = ib2 - ib_start+1
524 0 : Bsym%trCalc_irreps(idg)%mat(jb1,jb2,tr_isym)=cmat_ab
525 : end do ! tr_isym
526 : end if
527 :
528 : end do !ib2
529 : end do !ib1
530 :
531 : ! Calculate the trace for each class.
532 0 : if (Bsym%only_trace) then ! TODO this is valid if only trace.
533 0 : ABI_ERROR("Have to reconstruct missing traces")
534 : else
535 0 : do isym=1,Bsym%nsym_gk
536 0 : Bsym%Calc_irreps(idg)%trace(isym) = get_trace( Bsym%Calc_irreps(idg)%mat(:,:,isym) )
537 : end do
538 0 : if (Bsym%can_use_tr) then
539 0 : do tr_isym=1,Bsym%nsym_trgk
540 0 : Bsym%trCalc_irreps(idg)%trace(tr_isym) = get_trace( Bsym%trCalc_irreps(idg)%mat(:,:,tr_isym) )
541 : end do
542 : end if
543 : end if
544 :
545 : end do ! idg
546 :
547 0 : call Bsym%finalize(Wfd%prtvol)
548 0 : call Bsym%print([std_out, ab_out], prtvol=Wfd%prtvol)
549 :
550 : ! Free memory
551 0 : ABI_FREE(irottb)
552 0 : ABI_FREE(ur1)
553 0 : ABI_FREE(ur2)
554 0 : ABI_FREE(ur2_rot)
555 0 : ABI_FREE(eig0r)
556 0 : ABI_SFREE(tr_emig0r)
557 :
558 0 : if (Wfd%usepaw==1) then
559 0 : call pawcprj_free(Cprj_b1)
560 0 : ABI_FREE(Cprj_b1)
561 0 : call pawcprj_free(Cprj_b2)
562 0 : ABI_FREE(Cprj_b2)
563 0 : call pawcprj_free(Cprj_b2rot)
564 0 : ABI_FREE(Cprj_b2rot)
565 0 : ABI_FREE(zarot)
566 0 : call pawfgrtab_free(Pawfgrtab)
567 0 : ABI_FREE(Pawfgrtab)
568 0 : call paw_pwaves_lmn_free(Paw_onsite)
569 0 : ABI_FREE(Paw_onsite)
570 : end if
571 :
572 0 : end subroutine classify_bands
573 : !!***
574 :
575 : !----------------------------------------------------------------------
576 :
577 : !!****f* m_classify_bands/rotate_cprj
578 : !! NAME
579 : !! rotate_cprj
580 : !!
581 : !! FUNCTION
582 : !! Rotate cprj matrix elements by applying the symmetry operation of index isym
583 : !! that preserves the given k-point within a reciprocal lattice vector.
584 : !!
585 : !! INPUTS
586 : !! isym=index of the symmetry in the symrec arrays that preserves the given k-point within a reciprocal lattice vector
587 : !! ntypat=number of types of atom.
588 : !! natom=number of atoms.
589 : !! Cryst<crystal_t>=Datatype gathering info on the unit cell.
590 : !! typat(natom)=type of each atom.
591 : !! nbnds=number of bands for this k-point ans spin
592 : !! Cprj_in(natom,nbnds)<type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk>
593 : !! with all NL projectors at fixed k-point
594 : !!
595 : !! OUTPUT
596 : !! Cprj_out(natom,nbnds) <type(pawcprj_type)>= projection of the smooth PAW wave function onto
597 : !! projectors centered on equivalent sites of the crystal (non restricted to be in the firs unit cell)
598 : !! The equivalent site is defined according to the symmetry operation isym. Thus Cprj_out contains
599 : !!
600 : !! Cprj_out(at,b)=<p_j^{R^{-1}(L_{at}-\tau)} | \tpsi_b> if R is the isym operation with fractional translation \tau
601 : !! L_{at} is the position of the initial atom inside the first unit cell
602 : !! Note that atom a might be in a cell different from the initial one. No wrapping is done.
603 : !!
604 : !! SOURCE
605 :
606 0 : subroutine rotate_cprj(kpoint,isym,nspinor,nbnds,natom,nsym,typat,indsym,Cprj_in,Cprj_out)
607 :
608 : !Arguments ------------------------------------
609 : !scalars
610 : integer,intent(in) :: nbnds,nspinor,natom,isym,nsym
611 : !arrays
612 : integer,intent(in) :: typat(natom),indsym(4,nsym,natom)
613 : real(dp),intent(in) :: kpoint(3)
614 : type(pawcprj_type),intent(in) :: Cprj_in(natom,nspinor*nbnds)
615 : type(pawcprj_type),intent(out) :: Cprj_out(natom,nspinor*nbnds)
616 :
617 : !Local variables-------------------------------
618 : !scalars
619 : integer :: iat,iband,itypat,iat_sym
620 : real(dp) :: kdotr0
621 : !arrays
622 : integer :: r0(3)
623 : real(dp) :: phase_kr0(2)
624 : ! *************************************************************************
625 :
626 0 : do iat=1,natom
627 0 : itypat=typat(iat)
628 : ! The index of the symmetric atom.
629 : ! R^{-1} (xred(:,iat)-tnons) = xred(:,iat_sym) + r0.
630 : ! phase_kr0 takes into account the case in which rotated atom is in another unit cell.
631 0 : iat_sym=indsym(4,isym,iat); r0=indsym(1:3,isym,iat)
632 :
633 0 : kdotr0 = two_pi*DOT_PRODUCT(kpoint,r0)
634 0 : phase_kr0(1) = DCOS(kdotr0)
635 0 : phase_kr0(2) = DSIN(kdotr0)
636 :
637 : !phase_kr0 = (/one,zero/)
638 :
639 0 : do iband=1,nspinor*nbnds
640 : Cprj_out(iat,iband)%cp(1,:)= Cprj_in(iat_sym,iband)%cp(1,:)*phase_kr0(1) &
641 0 : -Cprj_in(iat_sym,iband)%cp(2,:)*phase_kr0(2)
642 :
643 : Cprj_out(iat,iband)%cp(2,:)= Cprj_in(iat_sym,iband)%cp(1,:)*phase_kr0(2) &
644 0 : +Cprj_in(iat_sym,iband)%cp(2,:)*phase_kr0(1)
645 : end do
646 : end do ! iat
647 :
648 0 : end subroutine rotate_cprj
649 : !!***
650 :
651 : !----------------------------------------------------------------------
652 :
653 : !!****f* m_classify_bands/paw_phirotphj
654 : !! NAME
655 : !! paw_phirotphj
656 : !!
657 : !! FUNCTION
658 : !! This routine calculates
659 : !! <\tPsi_1|\tprj_i> <\tprj_j|\tPsi_2> [ <\phi_i|\phi_j(R^{-1}r> - <\tphi_i|\tphi_j(R^{-1}r> ]
660 : !!
661 : !! [ <\phi_i|\phi_j(R^{-1}r> - <\tphi_i|\tphi_j(R^{-1}r> ] = s_ij D_{mi,mi}^{li}(R)
662 : !!
663 : !! INPUTS
664 : !! nspinor=Number of spinorial components.
665 : !! natom=number of atoms
666 : !! typat(natom)=type of eahc atom
667 : !! zarot_isym
668 : !! Pawtab(ntypat)<Pawtab_type>=PAW tabulated starting data
669 : !! Psps<pseudopotential_type>=Info on pseudopotentials.
670 : !! Cprj_b1(natom,nspinor)<type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk>
671 : !! with all NL projectors at fixed k-point
672 : !! Cprj_b2(natom,nspinor)<type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk>
673 : !! with all NL projectors at fixed k-point
674 : !! [conjg_left]=.TRUE if the complex conjugate of the left wavefunctions has to be taken. Defaults to .FALSE.
675 : !!
676 : !! OUTPUT
677 : !! omat(2)=The onsite matrix element.
678 : !!
679 : !! SOURCE
680 :
681 0 : function paw_phirotphj(nspinor,natom,typat,zarot_isym,Pawtab,Psps,Cprj_b1,Cprj_b2,conjg_left) result(omat)
682 :
683 : !Arguments ------------------------------------
684 : !scalars
685 : integer,intent(in) :: nspinor,natom
686 : logical,optional,intent(in) :: conjg_left
687 : type(pseudopotential_type),intent(in) :: Psps
688 : !arrays
689 : integer,intent(in) :: typat(natom)
690 : real(dp),intent(in) :: zarot_isym(:,:,:)
691 : real(dp) :: omat(2)
692 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
693 : type(pawcprj_type),intent(in) :: Cprj_b1(natom,nspinor),Cprj_b2(natom,nspinor)
694 :
695 : !Local variables-------------------------------
696 : !scalars
697 : integer :: iat,il,ilmn,ilpm,im,itypat,jl,jlmn,jlpm,jm,k0lmn,klmn,nlmn
698 : real(dp) :: dmimj,fij,im_p,re_p,sij
699 : logical :: do_conjg_left
700 : ! *************************************************************************
701 :
702 0 : do_conjg_left = .FALSE.; if (PRESENT(conjg_left)) do_conjg_left = conjg_left
703 :
704 0 : if (nspinor/=1) then
705 0 : ABI_ERROR("nspinor/=1 not yet coded")
706 : end if
707 :
708 : ! === Rotate PAW projections ===
709 : ! * zarot_isym is the rotation matrix of real spherical harmonics associated to symrec(:,:,isym).
710 : ! * zarot_isym multiply harmonics as row vectors, we need R^{-1} but we read R and invert m,mp in the equation below
711 0 : omat=zero
712 :
713 0 : do iat=1,natom
714 0 : itypat=typat(iat)
715 0 : nlmn=Pawtab(itypat)%lmn_size
716 :
717 0 : do jlmn=1,nlmn
718 0 : k0lmn=jlmn*(jlmn-1)/2
719 0 : jl=Psps%indlmn(1,jlmn,itypat)
720 0 : jm=Psps%indlmn(2,jlmn,itypat)
721 0 : jlpm=1+jl+jm
722 :
723 0 : do ilmn=1,jlmn
724 0 : il=Psps%indlmn(1,ilmn,itypat)
725 0 : im=Psps%indlmn(2,ilmn,itypat)
726 0 : if (il/=jl.or.im/=jm) CYCLE ! Selection rule on l and m.
727 0 : ilpm=1+il+im
728 :
729 0 : klmn=k0lmn+ilmn
730 0 : sij=Pawtab(itypat)%sij(klmn) !; if (ABS(sij)<tol14) CYCLE
731 :
732 : ! Here we get the matrix associated to R^{-1}.
733 0 : dmimj=zarot_isym(ilpm,jlpm,jl+1)
734 :
735 0 : if (do_conjg_left) then ! take the complex conjugate of the left cprj.
736 : re_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) &
737 : -Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) &
738 : +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) &
739 0 : -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn)
740 :
741 : im_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) &
742 : +Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) &
743 : -Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) &
744 0 : -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn)
745 : else
746 : re_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) &
747 : +Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) &
748 : +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) &
749 0 : +Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn)
750 :
751 : im_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) &
752 : -Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) &
753 : +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) &
754 0 : -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn)
755 : end if
756 : ! Accumulate the atom-centered contributions.
757 0 : fij = Pawtab(itypat)%dltij(klmn)/two
758 0 : omat(1)= omat(1) + fij*sij*re_p*dmimj
759 0 : omat(2)= omat(2) + fij*sij*im_p*dmimj
760 :
761 : end do !ilmn
762 : end do !jlmn
763 : end do !iat
764 :
765 0 : end function paw_phirotphj
766 : !!***
767 :
768 : !----------------------------------------------------------------------
769 :
770 : !!****f* m_classify_bands/dmats_init
771 : !! NAME
772 : !! dmats_init
773 : !!
774 : !! FUNCTION
775 : !! Compute D_mn(S) = <psi_{mSk}| S | psi_{nk}> for all the k-points in the IBZ and the bands in brange_spin.
776 : !!
777 : !! INPUTS
778 : !! wfk_path=Filename of the WFK file.
779 : !!
780 : !! NOTES
781 : !! Little-group membership and the umklapp vector G_0 associated to each symmetry
782 : !! are determined with the SAME k-point convention used by cgtk_rotate's own
783 : !! bookkeeping (see its docstring in m_cgtk.F90), namely:
784 : !!
785 : !! k2 = T symrel(:,:,isym)^t k1 + G_0 (transpose of symrel, NOT symrec)
786 : !!
787 : !! where T=+1/-1 without/with time reversal. This is also the convention produced
788 : !! by listkk (default, symrel-based) and consumed by cgtk_rotate elsewhere in the
789 : !! code (e.g. m_wfd.F90), so it is safe to reuse directly here.
790 : !!
791 : !! However, cgtk_rotate cannot simply be called with S=isym to obtain D(S_isym):
792 : !! its actual G-sphere index map is cg2(G) = cg1(symrec(isym).(G+G_0)), with
793 : !! symrec(isym) = mati3inv(symrel(isym)) = symrel(isym)^{-t} applied FORWARD
794 : !! (no additional inversion). Deriving the Fourier-coefficient transform of
795 : !! psi(r) -> psi(symrel(isym)^{-1}(r-tau)) shows that the coefficient at the
796 : !! rotated G must instead be read at symrel(isym)^t . G. The two matrices,
797 : !! symrel(isym)^{-t} and symrel(isym)^t, coincide only when symrel(isym) is an
798 : !! involution (S^2 = E, e.g. the identity or spatial inversion). For any other
799 : !! operation (3-, 4-, 6-fold rotations, screw axes, glide planes, ...) calling
800 : !! cgtk_rotate(isym) therefore silently returns D_true(S_isym)^{-1} = D_true(S_isym^{-1})
801 : !! instead of D_true(S_isym).
802 : !!
803 : !! This was confirmed empirically: with cgtk_rotate called on isym directly, the
804 : !! group-multiplication test in dmats_check (D(S_1 S_2) \propto D(S_1) D(S_2), see
805 : !! below) failed for essentially every triple involving a non-involutory operation,
806 : !! while unitarity and the D(S^{-1})=D(S)^dagger self-consistency test still passed
807 : !! (an involution-blind bug: D_true(S)^{-1} is unitary and equals D_true(S)^{-1}
808 : !! trivially, so those two checks cannot detect it). Concretely, for a triple
809 : !! (S_1, S_2, S_3=S_1 S_2) with zero fractional translations, the stored matrices
810 : !! satisfied D(S_3) = D(S_1) D(S_2)^t rather than D(S_3) = D(S_1) D(S_2).
811 : !!
812 : !! The fix is to call cgtk_rotate with isym_inv, the group-theoretic inverse of
813 : !! isym (found from the symrel multiplication table), while still filling the
814 : !! storage slot for isym: D_computed(isym_inv) = D_true(isym_inv^{-1}) = D_true(isym).
815 : !! isym_inv's own G_0 is recomputed with the formula above (using isym_inv instead
816 : !! of isym); no ad-hoc override of its fractional translation is needed, since
817 : !! cgtk_rotate is now called honestly for the operation it is actually asked to
818 : !! apply. See the inline comments in the k-point/symmetry loop below for the
819 : !! implementation.
820 : !!
821 : !! SOURCE
822 :
823 0 : subroutine dmats_init(dmats, wfk_path, dtset, cryst, brange_spin, ngfft, pawtab, psps, comm)
824 :
825 : !Arguments ------------------------------------
826 : !scalars
827 : class(dmats_t),intent(out) :: dmats
828 : character(len=*),intent(in) :: wfk_path
829 : type(dataset_type),target,intent(in) :: dtset
830 : class(crystal_t),target,intent(in) :: cryst
831 : integer,intent(in) :: brange_spin(2, dtset%nsppol), ngfft(18)
832 : integer,intent(in) :: comm
833 : type(pseudopotential_type),intent(in) :: psps
834 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
835 :
836 : !Local variables-------------------------------
837 : !scalars
838 : integer,parameter :: iflag1 = 1, me_g0 = 1, ndat1 = 1
839 : integer :: spin, nsppol, nsym, nb, nkibz, mband, ik_ibz, isym, isym_inv, itime, bstart, trev_k
840 : logical,parameter :: DEBUG_DUMP_SPINROT = .False.
841 : real(dp) :: spinrot_dbg(4)
842 : integer :: ib1, ib2, band1, band2, n1, n2, n3, n4, n5, n6, nfft, nspinor, mpw, my_mpw, ii, ipw
843 : integer :: nprocs, me, itot, ierr
844 : logical :: is_little_group
845 : real(dp),parameter :: xnorm1 = one
846 : real(dp) :: e_b1, e_b2, cpu, wall, gflops, tsign
847 0 : type(wfd_t) :: wfd
848 0 : type(hdr_type) :: hdr
849 0 : type(krank_t) :: my_krank
850 : ! Non-little-group (star-crossing) D-matrix computation.
851 : integer :: ik_ibz_p, isym_p, trev_p, isym_inv_S, npw_kp
852 : integer :: g0_ket(3)
853 0 : integer,allocatable :: kg_kp(:,:)
854 0 : real(dp),allocatable :: cg_bra(:,:,:), cg_ket(:,:), cg_work_p(:,:)
855 : !arrays
856 : integer :: g0_k(3), g0_k_inv(3), gmax(3), my_gmax(3), work_ngfft(18), units(2)
857 : integer :: indkk_p(6,1)
858 0 : integer,allocatable :: nband(:,:), wfd_istwfk(:)
859 : real(dp) :: kk_ibz(3), kk_sk(3), kk_sk_inv(3), dot(2)
860 0 : real(dp),allocatable :: cg_ib(:,:,:), cg_work(:,:), work(:,:,:,:), cg2_sk(:,:)
861 : complex(dp) :: cval
862 0 : complex(dp),allocatable :: cmat(:,:)
863 0 : logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:)
864 : !----------------------------------------------------------------------
865 :
866 0 : units = [std_out, ab_out]
867 0 : call cwtime(cpu, wall, gflops, "start")
868 0 : call wrtout(units, sjoin(" Computing dmats with symsigma_de", ftoa(dtset%symsigma_de * Ha_meV), " meV"))
869 :
870 0 : ABI_CHECK_IEQ(dtset%usepaw, 0, "PAW not coded!")
871 : ! NB: the D-matrix construction below is nspinor-agnostic: it is built entirely on top of
872 : ! cgtk_rotate (which already implements the SU(2) spin rotation + time-reversal spin-exchange
873 : ! for nspinor=2, see m_cgtk.F90) and cg_zdotc (a plain length-parametrized dot product that,
874 : ! called with n=npw_k*nspinor, sums over the spinor components as well as the G-vectors).
875 : ! So nspinor=2 (SOC) is supported here. The Kramers/Theta^2 VALIDATION diagnostics in
876 : ! dmats_check_one_k (only reachable via dmats%check/dmats%check_star, i.e. wfk_task "classify",
877 : ! never via gstore_symmetrize) are a separate matter and remain scalar-only for now: they encode
878 : ! an analytic phase formula derived assuming Theta^2=+1, which flips to Theta^2=-1 for spinors
879 : ! and has not been re-derived/verified yet -- see the explicit nspinor==1 guard added there.
880 :
881 : ! Read KS energies from the WFK file.
882 0 : dmats%ks_ebands = ebands_from_file(wfk_path, comm)
883 0 : dmats%cryst => cryst
884 0 : dmats%dtset => dtset
885 :
886 0 : nsppol = dmats%ks_ebands%nsppol; nsym = cryst%nsym; nkibz = dmats%ks_ebands%nkpt
887 0 : nprocs = xmpi_comm_size(comm); me = xmpi_comm_rank(comm)
888 :
889 : if (DEBUG_DUMP_SPINROT .and. xmpi_comm_rank(comm) == 0) then
890 : open(unit=792, file="spinrot_debug.csv", status="replace", action="write")
891 : do isym=1,nsym
892 : call getspinrot(cryst%rprimd, spinrot_dbg, cryst%symrel(:,:,isym))
893 : write(792,'(i0,1x,4(es24.16,1x))') isym, spinrot_dbg
894 : end do
895 : close(792)
896 : end if
897 :
898 0 : ABI_MALLOC(dmats%brange_spin, (2, nsppol))
899 0 : dmats%brange_spin = brange_spin
900 :
901 : ! Compute multiplication table.
902 0 : ABI_MALLOC(dmats%multable, (4, nsym, nsym))
903 0 : ABI_MALLOC(dmats%toinv, (4, nsym))
904 :
905 : call sg_multable(nsym, cryst%symafm, cryst%symrel, ierr, &
906 0 : tnons=cryst%tnons, multable=dmats%multable, toinv=dmats%toinv)
907 0 : ABI_CHECK_IEQ(ierr, 0, "sg_multable returned ierr != 0. See messages above.")
908 :
909 : ! Initialize the wave function descriptor.
910 0 : mband = maxval(brange_spin(2, :))
911 0 : ABI_MALLOC(nband, (nkibz, nsppol))
912 0 : ABI_MALLOC(bks_mask, (mband, nkibz, nsppol))
913 0 : ABI_MALLOC(keep_ur, (mband, nkibz, nsppol))
914 0 : nband = mband; bks_mask = .False.; keep_ur = .False.
915 :
916 : ! Every rank loads the FULL band range for EVERY IBZ k-point (not just the subset this
917 : ! rank "owns" for the D-matrix computation loop below). This is required by the
918 : ! non-little-group branch (see NOTES above): computing D_mn(S) for a general S needs the
919 : ! wavefunction at k'_ibz = S.k_ibz's own IBZ representative, which is not known until the
920 : ! symmetry loop runs and can be ANY IBZ point, not just the ones this rank would otherwise
921 : ! be assigned. The actual D-matrix COMPUTATION work (the ik_ibz loop below) is still
922 : ! MPI-distributed and xmpi_sum'd at the end -- only the wavefunction STORAGE is replicated.
923 0 : do spin=1,nsppol
924 0 : bks_mask(brange_spin(1,spin):brange_spin(2,spin), :, spin) = .True.
925 : end do
926 :
927 : ! Impose istwfk = 1 for all k-points.
928 : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
929 0 : ABI_MALLOC(wfd_istwfk, (nkibz))
930 0 : wfd_istwfk = 1
931 :
932 : call wfd%init(cryst, pawtab, psps, keep_ur, mband, nband, nkibz, nsppol, bks_mask,&
933 : dtset%nspden, dtset%nspinor, dtset%ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, dmats%ks_ebands%kptns, ngfft,&
934 0 : dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
935 : !call wfd%print([std_out], header="Wavefunctions for DMATS calculation")
936 :
937 0 : ABI_FREE(nband)
938 0 : ABI_FREE(keep_ur)
939 0 : ABI_FREE(wfd_istwfk)
940 0 : ABI_FREE(bks_mask)
941 :
942 : ! Read wavefunctions from WFK file.
943 0 : call wfd%read_wfk(wfk_path, iomode_from_fname(wfk_path), out_hdr=hdr)
944 :
945 : ! Energy cutoff must be the same else matrices are not unitary.
946 0 : call hdr%vs_dtset(dtset)
947 0 : ABI_CHECK(abs(dtset%ecut - hdr%ecut) < tol6, "Input ecut should be equal to the value used in the WFK file.")
948 0 : call hdr%free()
949 :
950 : ! krank_t used to map k' = S.k_ibz (any BZ point, not just IBZ ones) back to the IBZ,
951 : ! symrel^t convention, needed by the non-little-group D-matrix computation below.
952 0 : call my_krank%init(nkibz, dmats%ks_ebands%kptns)
953 :
954 : ! Compute max |G_i| to build the box.
955 0 : gmax = 0; mpw = 0
956 0 : do ik_ibz=1,nkibz
957 0 : if (.not. allocated(wfd%kdata(ik_ibz)%kg_k)) cycle
958 0 : associate (npw_k => wfd%npwarr(ik_ibz), kg_k => wfd%kdata(ik_ibz)%kg_k)
959 0 : mpw = max(mpw, npw_k)
960 0 : do ipw=1,npw_k
961 0 : do ii=1,3
962 0 : gmax(ii) = max(gmax(ii), abs(kg_k(ii,ipw)))
963 : end do
964 : end do
965 : end associate
966 : end do
967 0 : my_gmax = gmax; call xmpi_max(my_gmax, gmax, comm, ierr)
968 0 : my_mpw = mpw; call xmpi_max(my_mpw, mpw, comm, ierr)
969 :
970 : ! Init work_ngfft
971 0 : gmax = gmax + 4 ! FIXME: this is to account for umklapp
972 0 : gmax = 2*gmax + 1
973 0 : call ngfft_seq(work_ngfft, gmax)
974 : !write(std_out,*)"work_ngfft(1:3): ",work_ngfft(1:3)
975 0 : ABI_MALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6)))
976 :
977 0 : n1 = work_ngfft(1); n2 = work_ngfft(2); n3 = work_ngfft(3); n4 = work_ngfft(4); n5 = work_ngfft(5); n6 = work_ngfft(6)
978 0 : nfft = n1 * n2 * n3
979 0 : nspinor = wfd%nspinor
980 :
981 : ! Allocate D matrices for each spin on each proc and fill with zeros as we will MPI sum at the end.
982 0 : ABI_MALLOC(dmats%for_spin, (nsppol))
983 0 : do spin=1,nsppol
984 0 : nb = brange_spin(2,spin) - brange_spin(1,spin) + 1
985 0 : ABI_CALLOC(dmats%for_spin(spin)%value, (nb, nb, nsym, 2, nkibz))
986 : end do
987 :
988 0 : do spin=1,nsppol
989 0 : bstart = brange_spin(1, spin)
990 0 : nb = brange_spin(2, spin) - brange_spin(1, spin) + 1
991 0 : ABI_MALLOC(cmat, (nb, nb))
992 :
993 : ! Loop over k-points in the IBZ.
994 0 : do ik_ibz=1,nkibz
995 0 : itot = ik_ibz + (spin - 1)*nkibz; if (mod(itot - 1, nprocs) /= me) cycle ! MPI parallelism.
996 :
997 : ! NB: istwf_k is always 1 here. See call to wfd%init.
998 0 : associate (npw_k => wfd%npwarr(ik_ibz), istwf_k => wfd%kdata(ik_ibz)%istwfk, kg_k => wfd%kdata(ik_ibz)%kg_k)
999 0 : kk_ibz = dmats%ks_ebands%kptns(:, ik_ibz)
1000 :
1001 : ! Copy wavefunctions for this k-point.
1002 0 : ABI_MALLOC(cg_work, (2, npw_k*nspinor))
1003 0 : ABI_MALLOC(cg_ib, (2, npw_k*nspinor, nb))
1004 0 : do ib1=1,nb
1005 0 : band1 = ib1 + bstart - 1
1006 0 : call wfd%copy_cg(band1, ik_ibz, spin, cg_ib(:,:,ib1))
1007 : end do
1008 :
1009 : ! Loop over time-reversal and spatial symmetries.
1010 : ! g0_k is built with the SAME k-point convention assumed by cgtk_rotate's
1011 : ! bookkeeping: k2 = T symrel(:,:,isym)^t k1 + g0 (transpose of symrel, not symrec).
1012 : ! cgtk_rotate is then called with isym_inv rather than isym: see the detailed
1013 : ! explanation of why this is required (cgtk_rotate(isym) actually returns
1014 : ! D_true(S_isym)^{-1}, invisibly for involutions) in the NOTES of this subroutine's
1015 : ! SOURCE header above.
1016 0 : do itime=1,2
1017 0 : tsign = merge(one, -one, itime == 1)
1018 0 : trev_k = itime - 1
1019 0 : do isym=1,cryst%nsym
1020 : ! Compute cmat(b,b')
1021 0 : cmat = zero
1022 :
1023 0 : kk_sk = tsign * matmul(transpose(real(cryst%symrel(:,:,isym), dp)), kk_ibz)
1024 0 : g0_k = nint(kk_ibz - kk_sk)
1025 0 : is_little_group = all(abs(kk_ibz - kk_sk - g0_k) < tol8)
1026 :
1027 0 : if (.not. is_little_group) then
1028 : ! Sk = k' /= k + G: genuinely compute D_mn(S) = <psi_m,k'|S|psi_n,k> from real
1029 : ! WFK data (this used to be an identity placeholder -- never correct, only
1030 : ! harmless because no caller ever read this slot; gstore_symmetrize and
1031 : ! dmats_check_one_k/dmats_check_star only ever query genuine little-group
1032 : ! slots). Two independent rotations are needed, expressed on a COMMON target
1033 : ! G-sphere built at the RAW (unreduced) k' = kk_sk:
1034 : !
1035 : ! (a) BRA |psi_m,k'>: reconstructed via kpts_map's own canonical (isym_p,
1036 : ! trev_p, g0_p) mapping k'_ibz -> k', calling cgtk_rotate with isym_p
1037 : ! DIRECTLY (no toinv step). This mirrors wfd_sym_ug_kg_npw's (m_wfd.F90)
1038 : ! own validated convention for "reconstruct THE wavefunction at a BZ
1039 : ! point" -- the same machinery gstore_compute already uses in production.
1040 : ! Any internally-consistent representative of |psi,k'> works here since
1041 : ! the bra is a STATE, not a specific operator's image (gauge-flexible).
1042 : !
1043 : ! (b) KET S|psi_n,k>: this DOES need the honest action of the SPECIFIC
1044 : ! operator S (not just some state), so the toinv workaround documented
1045 : ! above is required: call cgtk_rotate with X=toinv(1,isym), the SAME
1046 : ! tsign/trev_k as S. The g0 needed here is NOT zero in general (an
1047 : ! earlier version of this code assumed so and was wrong -- caught by a
1048 : ! temporary self-consistency check routing known little-group cases
1049 : ! through this same general formula and comparing against the existing,
1050 : ! validated little-group result: 4/208 disagreed by O(1), all at points
1051 : ! with a non-involutory isym AND a nonzero little-group umklapp). Root
1052 : ! cause: cgtk_rotate's documented k2=symrel(X)^t.kpt1+g0 "official"
1053 : ! bookkeeping formula does NOT, in general, track the k-point the CONTENT
1054 : ! is actually expressed at once X is swapped for toinv(X) -- that formula
1055 : ! only happens to hold for the little-group case because target=source
1056 : ! there, collapsing two genuinely different quantities into one. The
1057 : ! general, correct g0 (verified via the same self-check, 208/208 exact to
1058 : ! numerical noise): treat kk_sk (my desired physical target, built from
1059 : ! S's OWN forward action) the same way the little-group formula treats
1060 : ! kk_ibz -- i.e. g0 = nint(kk_sk - kk_sk_inv), kk_sk_inv being
1061 : ! toinv(isym)'s own forward action on kk_ibz (the direct generalization of
1062 : ! the existing little-group branch's kk_sk_inv/g0_k_inv, with kk_sk
1063 : ! substituted for kk_ibz as the reference point).
1064 0 : ierr = kpts_map("symrel", dtset%kptopt, cryst, my_krank, 1, kk_sk, indkk_p)
1065 0 : ABI_CHECK(ierr == 0, "Cannot find symmetric image of Sk in the IBZ")
1066 0 : ik_ibz_p = indkk_p(1,1); isym_p = indkk_p(2,1); trev_p = indkk_p(6,1)
1067 :
1068 0 : call get_kg(kk_sk, 1, dtset%ecut, cryst%gmet, npw_kp, kg_kp)
1069 :
1070 0 : ABI_MALLOC(cg_bra, (2, npw_kp*nspinor, nb))
1071 : associate (npw_kip => wfd%npwarr(ik_ibz_p), kg_kip => wfd%kdata(ik_ibz_p)%kg_k, &
1072 : istwf_kip => wfd%kdata(ik_ibz_p)%istwfk)
1073 0 : ABI_MALLOC(cg_work_p, (2, npw_kip*nspinor))
1074 0 : do ib1=1,nb
1075 0 : band1 = ib1 + bstart - 1
1076 0 : call wfd%copy_cg(band1, ik_ibz_p, spin, cg_work_p)
1077 : call cgtk_rotate(dmats%cryst, dmats%ks_ebands%kptns(:,ik_ibz_p), isym_p, trev_p, &
1078 : indkk_p(3:5,1), nspinor, ndat1, npw_kip, kg_kip, npw_kp, kg_kp, &
1079 0 : istwf_kip, 1, cg_work_p, cg_bra(:,:,ib1), work_ngfft, work)
1080 : end do
1081 0 : ABI_FREE(cg_work_p)
1082 : end associate
1083 :
1084 0 : isym_inv_S = dmats%toinv(1, isym)
1085 0 : kk_sk_inv = tsign * matmul(transpose(real(cryst%symrel(:,:,isym_inv_S), dp)), kk_ibz)
1086 0 : g0_ket = nint(kk_sk - kk_sk_inv)
1087 0 : ABI_MALLOC(cg_ket, (2, npw_kp*nspinor))
1088 :
1089 0 : do ib2=1,nb
1090 0 : band2 = ib2 + bstart - 1
1091 0 : e_b2 = dmats%ks_ebands%eig(band2, ik_ibz, spin)
1092 :
1093 : call cgtk_rotate(dmats%cryst, kk_ibz, isym_inv_S, trev_k, g0_ket, nspinor, ndat1, &
1094 0 : npw_k, kg_k, npw_kp, kg_kp, istwf_k, 1, cg_ib(:,:,ib2), cg_ket, work_ngfft, work)
1095 :
1096 0 : do ib1=1,nb
1097 0 : band1 = ib1 + bstart - 1
1098 0 : e_b1 = dmats%ks_ebands%eig(band1, ik_ibz_p, spin)
1099 :
1100 0 : cval = zero
1101 0 : if (abs(e_b2 - e_b1) <= dtset%symsigma_de) then
1102 0 : dot = cg_zdotc(npw_kp * nspinor, cg_bra(:,:,ib1), cg_ket)
1103 0 : cval = dot(1) + j_dpc * dot(2)
1104 : end if
1105 0 : cmat(ib1, ib2) = cval
1106 : end do
1107 : end do
1108 :
1109 0 : ABI_FREE(cg_ket)
1110 0 : ABI_FREE(cg_bra)
1111 0 : ABI_FREE(kg_kp)
1112 :
1113 : else
1114 : ! Find the group-theoretic inverse of isym.
1115 0 : isym_inv = dmats%toinv(1, isym)
1116 0 : ABI_CHECK(isym_inv /= 0, "Could not find inverse symmetry!")
1117 :
1118 0 : kk_sk_inv = tsign * matmul(transpose(real(cryst%symrel(:,:,isym_inv), dp)), kk_ibz)
1119 0 : g0_k_inv = nint(kk_ibz - kk_sk_inv)
1120 :
1121 0 : ABI_MALLOC(cg2_sk, (2, npw_k*nspinor))
1122 :
1123 0 : do ib2=1,nb
1124 0 : band2 = ib2 + bstart - 1
1125 0 : e_b2 = dmats%ks_ebands%eig(band2, ik_ibz, spin)
1126 :
1127 : ! Compute the periodic part of S |psi_nk>.
1128 : call cgtk_rotate(dmats%cryst, kk_ibz, isym_inv, trev_k, g0_k_inv, nspinor, ndat1, &
1129 : npw_k, kg_k, &
1130 0 : npw_k, kg_k, istwf_k, istwf_k, cg_ib(:,:,ib2), cg2_sk, work_ngfft, work)
1131 :
1132 0 : do ib1=1,nb
1133 0 : band1 = ib1 + bstart - 1
1134 0 : e_b1 = dmats%ks_ebands%eig(band1, ik_ibz, spin)
1135 :
1136 : ! Only if e_b1 == e_b2.
1137 0 : cval = zero
1138 0 : if (abs(e_b2 - e_b1) <= dtset%symsigma_de) then
1139 : ! Evaluate the mathematical overlap: D_{mn} = <psi_m | S | psi_n>.
1140 : ! For time-reversal symmetries (itime == 2), S is anti-unitary (S = K U).
1141 : ! cgtk_rotate has already fully evaluated S|psi_n> into cg2_sk, which includes
1142 : ! the complex-conjugation of both the structural phase and Fourier coefficients.
1143 : ! Therefore, cg_zdotc properly computes <psi_m | S \psi_n> = \sum C_m^* C_{rot}.
1144 : ! No additional complex conjugate is needed on the output `cval`.
1145 0 : dot = cg_zdotc(npw_k * nspinor, cg_ib(:,:,ib1), cg2_sk)
1146 0 : cval = dot(1) + j_dpc * dot(2)
1147 : end if
1148 :
1149 0 : cmat(ib1, ib2) = cval
1150 : end do ! ib1
1151 : end do ! ib2
1152 :
1153 0 : ABI_FREE(cg2_sk)
1154 : end if
1155 :
1156 : ! Save final matrix.
1157 0 : dmats%for_spin(spin)%value(:, :, isym, itime, ik_ibz) = cmat
1158 : end do ! isym
1159 : end do ! itime
1160 :
1161 0 : ABI_FREE(cg_ib)
1162 0 : ABI_FREE(cg_work)
1163 : end associate
1164 : end do ! ik_ibz
1165 :
1166 0 : ABI_FREE(cmat)
1167 : end do ! spin
1168 :
1169 0 : ABI_FREE(work)
1170 0 : call wfd%free()
1171 0 : call my_krank%free()
1172 :
1173 : ! Collect results on each MPI proc.
1174 0 : do spin=1,nsppol
1175 0 : call xmpi_sum(dmats%for_spin(spin)%value, comm, ierr)
1176 : end do
1177 :
1178 0 : call cwtime_report(" dmats_init:", cpu, wall, gflops)
1179 :
1180 0 : end subroutine dmats_init
1181 : !!***
1182 :
1183 : !----------------------------------------------------------------------
1184 :
1185 : !!****f* m_classify_bands/dmats_ncwrite
1186 : !! NAME
1187 : !! dmats_ncwrite
1188 : !!
1189 : !! FUNCTION
1190 : !! Write the D-matrices and their band ranges to an open NetCDF file.
1191 : !! One NetCDF group is created for each spin channel.
1192 : !!
1193 : !! INPUTS
1194 : !! ncid=NetCDF file identifier. The file must be in data mode on entry.
1195 : !!
1196 : !! SOURCE
1197 :
1198 0 : subroutine dmats_ncwrite(dmats, ncid)
1199 :
1200 : !Arguments ------------------------------------
1201 : class(dmats_t),target,intent(in) :: dmats
1202 : integer,intent(in) :: ncid
1203 :
1204 : !Local variables-------------------------------
1205 : integer :: spin, spin_ncid, ncerr, nb, nsym, nkibz
1206 0 : real(dp),contiguous,pointer :: dmat_ptr(:,:,:,:,:,:)
1207 : !----------------------------------------------------------------------
1208 :
1209 0 : nsym = dmats%cryst%nsym
1210 0 : nkibz = dmats%ks_ebands%nkpt
1211 :
1212 0 : do spin=1,size(dmats%for_spin)
1213 0 : nb = dmats%brange_spin(2, spin) - dmats%brange_spin(1, spin) + 1
1214 0 : if (.not. all(shape(dmats%for_spin(spin)%value) == [nb, nb, nsym, 2, nkibz])) then
1215 0 : ABI_ERROR("Inconsistent dmats array shape")
1216 : end if
1217 :
1218 0 : NCF_CHECK(nctk_set_defmode(ncid))
1219 0 : NCF_CHECK(nf90_def_grp(ncid, strcat("dmats", "_spin", itoa(spin)), spin_ncid))
1220 : ncerr = nctk_def_dims(spin_ncid, [ &
1221 : nctkdim_t("cplex_dmat", 2), &
1222 : nctkdim_t("two_dmat", 2), &
1223 : nctkdim_t("nb_dmat", nb), &
1224 : nctkdim_t("nsym_dmat", nsym), &
1225 : nctkdim_t("ntime_dmat", 2), &
1226 : nctkdim_t("nkibz_dmat", nkibz) &
1227 0 : ], defmode=.True.)
1228 0 : NCF_CHECK(ncerr)
1229 : ncerr = nctk_def_arrays(spin_ncid, [ &
1230 : nctkarr_t("band_range", "int", "two_dmat"), &
1231 : nctkarr_t("dmat_values", "dp", "cplex_dmat, nb_dmat, nb_dmat, nsym_dmat, ntime_dmat, nkibz_dmat") &
1232 0 : ])
1233 0 : NCF_CHECK(ncerr)
1234 0 : NCF_CHECK(nctk_set_datamode(spin_ncid))
1235 0 : NCF_CHECK(nf90_put_var(spin_ncid, nctk_idname(spin_ncid, "band_range"), dmats%brange_spin(:, spin)))
1236 0 : call c_f_pointer(c_loc(dmats%for_spin(spin)%value), dmat_ptr, [2, nb, nb, nsym, 2, nkibz])
1237 0 : NCF_CHECK(nf90_put_var(spin_ncid, nctk_idname(spin_ncid, "dmat_values"), dmat_ptr))
1238 : end do
1239 :
1240 0 : end subroutine dmats_ncwrite
1241 : !!***
1242 :
1243 : !----------------------------------------------------------------------
1244 :
1245 : !!****f* m_classify_bands/dmats_init_from_file
1246 : !! NAME
1247 : !! dmats_init_from_file
1248 : !!
1249 : !! FUNCTION
1250 : !! Initialize a dmats_t object from D-matrices stored by dmats_ncwrite.
1251 : !! Crystal and dataset pointers are associated with the caller-owned objects.
1252 : !!
1253 : !! INPUTS
1254 : !! path=Path to the NetCDF file containing dmats_spinN groups.
1255 : !! dtset=Dataset used to validate the number of spin channels.
1256 : !! cryst=Crystal used to rebuild the symmetry multiplication tables.
1257 : !! comm=MPI communicator used to read the electronic band structure.
1258 : !!
1259 : !! SOURCE
1260 :
1261 0 : subroutine dmats_init_from_file(dmats, path, dtset, cryst, comm)
1262 :
1263 : !Arguments ------------------------------------
1264 : class(dmats_t),target,intent(out) :: dmats
1265 : character(len=*),intent(in) :: path
1266 : type(dataset_type),target,intent(in) :: dtset
1267 : class(crystal_t),target,intent(in) :: cryst
1268 : integer,intent(in) :: comm
1269 :
1270 : !Local variables-------------------------------
1271 : integer :: ncid, spin_ncid, spin, nsppol, nsym, nkibz, nb, ncerr, ierr
1272 : integer :: nb_file, nsym_file, ntime_file, nkibz_file, cplex_file
1273 0 : real(dp),contiguous,pointer :: dmat_ptr(:,:,:,:,:,:)
1274 : !----------------------------------------------------------------------
1275 :
1276 0 : dmats%ks_ebands = ebands_from_file(path, comm)
1277 0 : dmats%cryst => cryst
1278 0 : dmats%dtset => dtset
1279 :
1280 0 : nsppol = dmats%ks_ebands%nsppol
1281 0 : nsym = cryst%nsym
1282 0 : nkibz = dmats%ks_ebands%nkpt
1283 0 : ABI_CHECK_IEQ(nsppol, dtset%nsppol, "Inconsistent nsppol in dmats file")
1284 :
1285 0 : ABI_MALLOC(dmats%brange_spin, (2, nsppol))
1286 0 : ABI_MALLOC(dmats%multable, (4, nsym, nsym))
1287 0 : ABI_MALLOC(dmats%toinv, (4, nsym))
1288 : call sg_multable(nsym, cryst%symafm, cryst%symrel, ierr, &
1289 0 : tnons=cryst%tnons, multable=dmats%multable, toinv=dmats%toinv)
1290 0 : ABI_CHECK_IEQ(ierr, 0, "sg_multable returned ierr !=0. See messages above.")
1291 :
1292 0 : ABI_MALLOC(dmats%for_spin, (nsppol))
1293 0 : NCF_CHECK(nctk_open_read(ncid, path, xmpi_comm_self))
1294 0 : do spin=1,nsppol
1295 0 : NCF_CHECK(nf90_inq_ncid(ncid, strcat("dmats", "_spin", itoa(spin)), spin_ncid))
1296 0 : NCF_CHECK(nctk_get_dim(spin_ncid, "cplex_dmat", cplex_file))
1297 0 : NCF_CHECK(nctk_get_dim(spin_ncid, "nb_dmat", nb_file))
1298 0 : NCF_CHECK(nctk_get_dim(spin_ncid, "nsym_dmat", nsym_file))
1299 0 : NCF_CHECK(nctk_get_dim(spin_ncid, "ntime_dmat", ntime_file))
1300 0 : NCF_CHECK(nctk_get_dim(spin_ncid, "nkibz_dmat", nkibz_file))
1301 0 : ABI_CHECK_IEQ(cplex_file, 2, "dmats file should contain complex matrices")
1302 0 : ABI_CHECK_IEQ(nsym_file, nsym, "Inconsistent number of symmetries in dmats file")
1303 0 : ABI_CHECK_IEQ(ntime_file, 2, "Inconsistent time-reversal dimension in dmats file")
1304 0 : ABI_CHECK_IEQ(nkibz_file, nkibz, "Inconsistent number of IBZ k-points in dmats file")
1305 :
1306 0 : NCF_CHECK(nf90_get_var(spin_ncid, nctk_idname(spin_ncid, "band_range"), dmats%brange_spin(:, spin)))
1307 0 : nb = dmats%brange_spin(2, spin) - dmats%brange_spin(1, spin) + 1
1308 0 : ABI_CHECK_IEQ(nb_file, nb, "Inconsistent band range in dmats file")
1309 0 : ABI_MALLOC(dmats%for_spin(spin)%value, (nb, nb, nsym, 2, nkibz))
1310 0 : call c_f_pointer(c_loc(dmats%for_spin(spin)%value), dmat_ptr, [2, nb, nb, nsym, 2, nkibz])
1311 0 : ncerr = nf90_get_var(spin_ncid, nctk_idname(spin_ncid, "dmat_values"), dmat_ptr)
1312 0 : NCF_CHECK(ncerr)
1313 : end do
1314 0 : NCF_CHECK(nf90_close(ncid))
1315 :
1316 0 : end subroutine dmats_init_from_file
1317 : !!***
1318 :
1319 : !----------------------------------------------------------------------
1320 :
1321 : !!****f* m_classify_bands/dmats_free
1322 : !! NAME
1323 : !! dmats_free
1324 : !!
1325 : !! FUNCTION
1326 : !! Free memory
1327 : !!
1328 : !! SOURCE
1329 :
1330 0 : subroutine dmats_free(dmats)
1331 :
1332 : !Arguments ------------------------------------
1333 : class(dmats_t),intent(inout) :: dmats
1334 :
1335 : !Local variables-------------------------------
1336 : integer :: spin
1337 : !----------------------------------------------------------------------
1338 :
1339 0 : call dmats%ks_ebands%free()
1340 :
1341 0 : ABI_SFREE(dmats%brange_spin)
1342 0 : ABI_SFREE(dmats%multable)
1343 0 : ABI_SFREE(dmats%toinv)
1344 :
1345 0 : do spin=1,size(dmats%for_spin)
1346 0 : ABI_SFREE(dmats%for_spin(spin)%value)
1347 : end do
1348 0 : ABI_SFREE(dmats%for_spin)
1349 :
1350 0 : end subroutine dmats_free
1351 : !!***
1352 :
1353 : !----------------------------------------------------------------------
1354 :
1355 : !!****f* m_classify_bands/dmats_check_one_k
1356 : !! NAME
1357 : !! dmats_check_one_k
1358 : !!
1359 : !! FUNCTION
1360 : !! Run the full battery of algebraic tests (see dmats_check's SOURCE header for the
1361 : !! complete list: unitarity, identity, inverse relation, group multiplication,
1362 : !! Kramers, class character, S^n closure) on a single (spin, k) slice of D-matrices.
1363 : !! Factored out of dmats_check so the SAME tests can be reused, unchanged, on
1364 : !! D-matrices reconstructed at an arbitrary full-BZ k-point (see dmats_check_star),
1365 : !! not just on the genuine per-IBZ slices of dmats%for_spin.
1366 : !!
1367 : !! INPUTS
1368 : !! spin=Spin index.
1369 : !! kk_ibz(3)=k-point (reduced coords) at which dmat_k was computed. Despite the name,
1370 : !! this need not be an actual IBZ point of dmats%ks_ebands: it is only used to
1371 : !! rebuild the little group (symtab) and the analytic phase formulas, both of
1372 : !! which only depend on the k-vector itself, not on how dmat_k was constructed.
1373 : !! dmat_k(:,:,:,:)=D-matrices (nb,nb,nsym,2) at kk_ibz: dmat_k(:,:,isym,itime) is the
1374 : !! identity placeholder if (isym,itime) is not in the little group of kk_ibz.
1375 : !! units(:), prtvol=Output units and verbosity.
1376 : !! tag=Label used as the yamldoc dictlist key for this k-point's diagnostics.
1377 : !!
1378 : !! SIDE EFFECTS
1379 : !! ydoc=yamldoc_t, appended to with this k-point's diagnostics dictlist.
1380 : !! ierr=Accumulated error counter (incremented on each failed sub-test).
1381 : !!
1382 : !! SOURCE
1383 :
1384 0 : subroutine dmats_check_one_k(dmats, spin, kk_ibz, dmat_k, units, prtvol, tag, ydoc, ierr)
1385 :
1386 : !Arguments ------------------------------------
1387 : class(dmats_t),intent(in) :: dmats
1388 : integer,intent(in) :: spin, units(:), prtvol
1389 : real(dp),intent(in) :: kk_ibz(3)
1390 : complex(dp),intent(in) :: dmat_k(:,:,:,:)
1391 : character(len=*),intent(in) :: tag
1392 : type(yamldoc_t),intent(inout) :: ydoc
1393 : integer,intent(inout) :: ierr
1394 :
1395 : !Local variables-------------------------------
1396 : integer :: nb, isym, itime, isym_inv, j, isym1, isym2, isym3, n, isym_cnt, ierr_so
1397 : integer :: itime1, itime2, itime3, nsym_lg, nclass_lg, icls, iel, il
1398 : logical :: unitary, identity_ok, kramers_ok, char_ok, isproper
1399 : character(len=5000) :: msg
1400 : real(dp),parameter :: DTOL = tol3
1401 : real(dp) :: kk_sk(3), tsign, err, phase_err, char_err
1402 : real(dp) :: Sk3(3), L_mult(3), phase_err_mult
1403 : complex(dp) :: phase_L, phase_analytic, phase_dyn, phase_analytic_mult
1404 : integer :: isym1_inv, isym2_inv
1405 : integer :: g0_k(3)
1406 0 : integer :: symtab(4,2,dmats%cryst%nsym)
1407 : integer :: mult_fail_cnt
1408 0 : complex(dp),allocatable :: cmat_n(:,:)
1409 0 : type(pair_list), allocatable :: sym_dicts(:), mult_fail_dicts(:)
1410 0 : integer :: sym_lg(3,3,dmats%cryst%nsym), local2global(dmats%cryst%nsym), trans(3)
1411 0 : integer :: class_id_of_isym(dmats%cryst%nsym)
1412 0 : integer,allocatable :: nelements_lg(:), elements_idx_lg(:,:)
1413 0 : real(dp) :: class_char_ref(dmats%cryst%nsym)
1414 : ! *************************************************************************
1415 :
1416 : ABI_UNUSED((/spin/))
1417 :
1418 : ! NB: unlike the Kramers/Theta^2 test below (isym==1, itime==2 block), which assumes
1419 : ! Theta^2=+1 and is therefore explicitly gated off for nspinor==2 right where it is
1420 : ! computed, every OTHER test in this routine (unitarity, identity, inverse relation,
1421 : ! group multiplication, character class, closure) is Theta^2-agnostic and safe to run
1422 : ! for nspinor==2 as a self-consistency diagnostic of dmats_init's spinor D-matrices.
1423 :
1424 0 : nb = size(dmat_k, 1)
1425 0 : ABI_MALLOC(cmat_n, (nb, nb))
1426 :
1427 : ! Determine the little group of kk_ibz (and the associated umklapp vector) with the SAME
1428 : ! symrel^t convention used by dmats_init to decide whether a given (isym,itime) yields a
1429 : ! genuinely-computed D-matrix or the identity placeholder (see the NOTES on g0_k there).
1430 : ! littlegroup_q instead rotates kk_ibz with symrec, which is the convention for q-points,
1431 : ! not k-points, and can disagree with dmats_init whenever symrel is not orthogonal in the
1432 : ! reduced basis -- inconsistent with what dmats%for_spin(...) actually stores.
1433 0 : symtab = 0
1434 0 : do itime=1,2
1435 0 : tsign = merge(one, -one, itime == 1)
1436 0 : do isym=1,dmats%cryst%nsym
1437 0 : kk_sk = tsign * matmul(transpose(real(dmats%cryst%symrel(:,:,isym), dp)), kk_ibz)
1438 0 : g0_k = nint(kk_ibz - kk_sk)
1439 0 : if (all(abs(kk_ibz - kk_sk - g0_k) < tol8)) then
1440 0 : symtab(1:3, itime, isym) = g0_k
1441 0 : symtab(4, itime, isym) = 1
1442 : end if
1443 : end do
1444 : end do
1445 :
1446 : ! Divide the itime=1 (pure spatial) little group into conjugacy classes with get_classes
1447 : ! (m_ptgroups.F90), then use |Tr D(S)| (character magnitude) as a class-function test:
1448 : ! conjugate elements S' = X S X^{-1} of a genuine (possibly projective) unitary
1449 : ! representation always satisfy |Tr D(S')| = |Tr D(S)| exactly, since
1450 : ! D(X S X^{-1}) equals D(X) D(S) D(X)^{-1} up to an overall SCALAR phase (the same
1451 : ! tabulated-vs-literal-composition phase ambiguity already handled in the group
1452 : ! multiplication test above), and a similarity transform composed with an overall
1453 : ! phase preserves |trace|. Comparing |trace| (not the raw complex trace) sidesteps
1454 : ! that phase ambiguity entirely, so no analytic phase tracking is needed here.
1455 : !
1456 : ! get_classes computes conjugacy X^{-1} S X with the plain (non-transposed) symrel
1457 : ! product, i.e. the SAME real-space composition convention already used everywhere
1458 : ! else in this routine (isym_inv search, isym3 = isym1*isym2 hunting). This is
1459 : ! convention-consistent with the D-matrices themselves: after the dmats_init fix
1460 : ! (see its NOTES), dmats%for_spin(...)(:,:,isym,...) genuinely stores D_true(S_isym)
1461 : ! indexed by the SAME isym used to index cryst%symrel, so no extra transpose or
1462 : ! symrec/symrel^t handling is required to match classes to D-matrix slots. (The
1463 : ! symrel^t convention only enters dmats_init's k-point/little-group bookkeeping;
1464 : ! little-group MEMBERSHIP and conjugacy-class PARTITIONING are both provably
1465 : ! independent of that choice: symrec = symrel^{-t} is a group isomorphism that maps
1466 : ! every isym to itself, so it preserves both the little-group isym set and the
1467 : ! class partition of that set exactly.)
1468 : !
1469 : ! Why NOT feed get_classes symrel(isym)^t either (transposed, but still indexed by
1470 : ! the SAME isym): unlike symrec = symrel^{-t}, plain transposition isym -> symrel(isym)^t
1471 : ! is only an ANTI-homomorphism of the isym-indexed abstract group law, because
1472 : ! (AB)^t = B^t A^t reverses multiplication order: symrel(a)^t symrel(b)^t = symrel(b.a)^t,
1473 : ! not symrel(a.b)^t, whenever the point group is non-abelian (as it generally is here).
1474 : ! Anti-homomorphisms still preserve conjugacy classes as an ABSTRACT structure, but not
1475 : ! with the SAME isym labeling used to index dmats%for_spin(...), so class_id_of_isym
1476 : ! built from transposed matrices would in general group the WRONG isym's together. Only
1477 : ! plain (non-transposed) symrel, matching how D(S1 S2) proportional-to D(S1) D(S2) was
1478 : ! validated to hold in the group-multiplication test above, gives isym-consistent classes.
1479 0 : nsym_lg = 0
1480 0 : do isym=1,dmats%cryst%nsym
1481 0 : if (symtab(4, 1, isym) == 0) cycle
1482 0 : nsym_lg = nsym_lg + 1
1483 0 : sym_lg(:,:,nsym_lg) = dmats%cryst%symrel(:,:,isym)
1484 0 : local2global(nsym_lg) = isym
1485 : end do
1486 :
1487 0 : class_id_of_isym = 0
1488 0 : if (nsym_lg > 0) then
1489 : ! get_classes takes explicit-shape dummies sized from its own nsym argument
1490 : ! (nsym_lg here): the actual arrays must be allocated to EXACTLY (nsym_lg,nsym_lg)
1491 : ! (not e.g. cryst%nsym), otherwise the callee writes using an nsym_lg-based
1492 : ! column-major stride while a differently-sized caller array would read back
1493 : ! with a mismatched stride (silent data corruption via sequence association).
1494 0 : ABI_MALLOC(nelements_lg, (nsym_lg))
1495 0 : ABI_MALLOC(elements_idx_lg, (nsym_lg, nsym_lg))
1496 0 : call get_classes(nsym_lg, sym_lg(:,:,1:nsym_lg), nclass_lg, nelements_lg, elements_idx_lg)
1497 0 : do icls=1,nclass_lg
1498 0 : do iel=1,nelements_lg(icls)
1499 0 : il = elements_idx_lg(iel, icls)
1500 0 : class_id_of_isym(local2global(il)) = icls
1501 : end do
1502 : end do
1503 : ! Reference character (magnitude) for each class: |Tr D(S)| for the class's first
1504 : ! element. All other elements of the same class are checked against this below.
1505 0 : do icls=1,nclass_lg
1506 0 : isym = local2global(elements_idx_lg(1, icls))
1507 0 : class_char_ref(icls) = abs(get_trace(dmat_k(:, :, isym, 1)))
1508 : end do
1509 0 : ABI_FREE(nelements_lg)
1510 0 : ABI_FREE(elements_idx_lg)
1511 : end if
1512 :
1513 0 : isym_cnt = 0
1514 0 : do itime=1,2
1515 0 : do isym=1,dmats%cryst%nsym
1516 0 : if (symtab(4, itime, isym) /= 0) isym_cnt = isym_cnt + 1
1517 : end do
1518 : end do
1519 0 : if (isym_cnt > 0) then
1520 0 : ABI_MALLOC(sym_dicts, (isym_cnt))
1521 : end if
1522 :
1523 0 : isym_cnt = 0
1524 0 : do itime=1,2
1525 0 : do isym=1,dmats%cryst%nsym
1526 0 : if (symtab(4, itime, isym) == 0) cycle
1527 0 : isym_cnt = isym_cnt + 1
1528 0 : call sym_dicts(isym_cnt)%set("isym", i=isym)
1529 0 : call sym_dicts(isym_cnt)%set("itime", i=itime)
1530 0 : msg = sjoin("[", ftoa(dmats%cryst%tnons(1,isym)), ", ", ftoa(dmats%cryst%tnons(2,isym)))
1531 0 : msg = sjoin(msg, ", ", ftoa(dmats%cryst%tnons(3,isym)), "]")
1532 0 : call sym_dicts(isym_cnt)%set("tnon", s=trim(msg))
1533 :
1534 0 : msg = sjoin("[", itoa(symtab(1, itime, isym)), ", ", itoa(symtab(2, itime, isym)))
1535 0 : msg = sjoin(msg, ", ", itoa(symtab(3, itime, isym)), ", ", itoa(symtab(4, itime, isym)), "]")
1536 0 : call sym_dicts(isym_cnt)%set("symtab", s=trim(msg))
1537 :
1538 0 : associate (cmat => dmat_k(:, :, isym, itime))
1539 0 : unitary = is_unitary(nb, cmat, DTOL, err)
1540 0 : if (.not. unitary) ierr = ierr + 1
1541 0 : call sym_dicts(isym_cnt)%set("unitary", s=yesno(unitary))
1542 0 : call sym_dicts(isym_cnt)%set("unitary_err", r=err)
1543 :
1544 : ! Identity operator test
1545 0 : if (isym == 1 .and. itime == 1) then
1546 0 : identity_ok = is_identity(nb, cmat, DTOL, err)
1547 0 : if (.not. identity_ok) ierr = ierr + 1
1548 0 : call sym_dicts(isym_cnt)%set("identity_ok", s=yesno(identity_ok))
1549 0 : call sym_dicts(isym_cnt)%set("identity_err", r=err)
1550 : end if
1551 :
1552 : ! Character class-function test (itime=1 only, see the NOTES on get_classes
1553 : ! above the class_id_of_isym computation): |Tr D(S)| must be the same for
1554 : ! every S in a given conjugacy class of the little group.
1555 0 : if (itime == 1 .and. class_id_of_isym(isym) /= 0) then
1556 0 : icls = class_id_of_isym(isym)
1557 0 : char_err = abs(abs(get_trace(cmat)) - class_char_ref(icls))
1558 0 : char_ok = (char_err < DTOL)
1559 0 : if (.not. char_ok) ierr = ierr + 1
1560 0 : call sym_dicts(isym_cnt)%set("class_id", i=icls)
1561 0 : call sym_dicts(isym_cnt)%set("char_ok", s=yesno(char_ok))
1562 0 : call sym_dicts(isym_cnt)%set("char_err", r=char_err)
1563 : end if
1564 :
1565 : ! Kramers test: pure time reversal (isym=1, itime=2) is only present in this
1566 : ! slot at TR-invariant k-points (k = -k mod G, e.g. TRIM points), and must
1567 : ! satisfy \Theta^2 = D(\Theta) D(\Theta)^* = +I EXACTLY (not just up to a
1568 : ! phase) for scalar (nspinor=1, hard-required by dmats_init) wavefunctions.
1569 : ! Unlike the generic group-multiplication test, this is an exact identity
1570 : ! with no residual gauge/tabulation-phase freedom: rescaling each band by an
1571 : ! arbitrary phase e^{i\phi_n} transforms D(\Theta) -> \Phi^{-1} D(\Theta) \Phi^{-1}
1572 : ! (antiunitary => the KET phase also gets conjugated), so
1573 : ! D(\Theta)D(\Theta)^* -> \Phi^{-1} [D(\Theta)D(\Theta)^*] \Phi, which leaves
1574 : ! "= I" invariant. This differs from, and is NOT redundant with, the inverse-
1575 : ! relation test below (isym_inv=1=isym for itime=2), which only checks that
1576 : ! D(\Theta) is proportional to its own transpose, not that D(\Theta)D(\Theta)^*=I.
1577 0 : if (isym == 1 .and. itime == 2 .and. dmats%dtset%nspinor == 1) then
1578 0 : kramers_ok = is_identity(nb, matmul(cmat, conjg(cmat)), DTOL, err)
1579 0 : if (.not. kramers_ok) ierr = ierr + 1
1580 0 : call sym_dicts(isym_cnt)%set("kramers_ok", s=yesno(kramers_ok))
1581 0 : call sym_dicts(isym_cnt)%set("kramers_err", r=err)
1582 : end if
1583 :
1584 : ! Inverse relation test
1585 0 : isym_inv = dmats%toinv(1, isym)
1586 :
1587 0 : if (isym_inv /= 0 .and. symtab(4, itime, isym_inv) /= 0) then
1588 : ! Analytic phase relating D(S^{-1}) to D(S)^\dagger, i.e. predicting phase_dyn (below). This
1589 : ! reuses the SAME tabulated-vs-literal-composition phase relation already validated to exact
1590 : ! (machine-precision) agreement by the group-multiplication test further down this routine:
1591 : !
1592 : ! matmul(D(isym1), D(isym2)) = e^{+i 2pi (S_{isym3} k).L} D(isym3),
1593 : ! L = multable(2:4, toinv(isym2), toinv(isym1)), isym3 = multable(1, isym1, isym2)
1594 : !
1595 : ! specialized to isym1=isym_inv, isym2=isym, so that the literal Seitz composition
1596 : ! S(isym_inv).S(isym) is EXACTLY the identity (isym3=1, S_{isym3} k = k), giving
1597 : ! L = multable(2:4, isym, isym_inv) = toinv(2:4, isym_inv) exactly (both tabulate the same
1598 : ! S(isym_inv).S(isym) = {I, L} relation), so no separate multable lookup is even needed here.
1599 : !
1600 : ! NOTE: the previous formula used toinv(2:4, isym) [i.e. the residual L for the OPPOSITE
1601 : ! composition order S(isym).S(isym_inv)] plus an extra ad hoc symrec/symtab correction term --
1602 : ! an isym/isym_inv index-order mixup that was the actual root cause of the long-standing
1603 : ! "disagrees with phase_dyn by a discrete 90/180 degree offset on non-symmorphic operations"
1604 : ! bug (see git history). Verified: phase_err now at machine precision for every (isym,itime)
1605 : ! tuple in the reference gstore test (previously ~58 tuples off by sqrt(2), ~80 by exactly 2).
1606 0 : phase_analytic = exp(cmplx(zero, two_pi * sum(kk_ibz * real(dmats%toinv(2:4, isym_inv), dp)), dp))
1607 : ! TEMPORARY DIAGNOSTIC (session experiment, not yet a verified fix): test whether the
1608 : ! itime=2, non-involutory inv_ok failures (100% clean split found empirically) are fixed
1609 : ! by an extra missing sign tied to itime and non-involutory character.
1610 0 : if (dmats%dtset%nspinor == 2 .and. itime == 2 .and. isym_inv /= isym) phase_analytic = -phase_analytic
1611 0 : if (dmats%dtset%nspinor == 2 .and. itime == 1 .and. isym_inv == isym .and. isym /= 1 .and. isym /= 2) &
1612 0 : phase_analytic = -phase_analytic
1613 :
1614 : associate (cmat_inv => dmat_k(:, :, isym_inv, itime))
1615 : ! Independently, dynamically extract the phase relating the two independently constructed
1616 : ! matrices by taking the Frobenius inner product of the two matrices:
1617 : ! Phase = Tr(A^\dagger B) / nb = sum_{ij} A^*_{ij} B_{ij} / nb.
1618 : ! If the matrices are truly proportional, Phase will be a scalar of unit magnitude,
1619 : ! and dividing by it will yield a mathematically exact equality test.
1620 0 : if (itime == 1) then
1621 0 : phase_dyn = sum( conjg(cmat_inv) * conjg(transpose(cmat)) ) / nb
1622 0 : err = maxval(abs(cmat_inv * (phase_dyn / abs(phase_dyn)) - conjg(transpose(cmat))))
1623 : else
1624 0 : phase_dyn = sum( conjg(cmat_inv) * transpose(cmat) ) / nb
1625 0 : err = maxval(abs(cmat_inv * (phase_dyn / abs(phase_dyn)) - transpose(cmat)))
1626 : end if
1627 : ! phase_dyn, by construction of the Frobenius inner product above, should equal conjg(phase_analytic)
1628 : ! when the D-matrices carry the correct absolute phase, so phase_dyn * phase_analytic == 1.
1629 : ! Formerly a diagnostic-only quantity (the old toinv(isym)-based formula disagreed with
1630 : ! phase_dyn by a discrete 90/180 degree offset on non-symmorphic operations); now that the
1631 : ! formula above is fixed (verified at machine precision on the full reference gstore test,
1632 : ! all 976 (isym,itime) tuples across all 8 IBZ k-points), phase_err is folded into ierr/inv_ok.
1633 0 : phase_err = abs(phase_dyn * phase_analytic - one)
1634 : ! TEMPORARY DIAGNOSTIC: for nspinor==2, itime==2, isym with a TRIVIAL rotation part
1635 : ! (symrel(isym) = +-Identity, i.e. E or pure spatial inversion -- both spin-blind), this
1636 : ! generic test degenerates to exactly the Kramers Theta^2=-1 signature (already correctly
1637 : ! excluded from the DEDICATED Kramers sub-test above for nspinor==2) -- exclude it from
1638 : ! gating ierr here too, so it doesn't falsely abort dmats%check()/block reaching check_star.
1639 0 : if (.not. (dmats%dtset%nspinor == 2 .and. itime == 2 .and. &
1640 : (all(dmats%cryst%symrel(:,:,isym) == identity_3d) .or. &
1641 : all(dmats%cryst%symrel(:,:,isym) == -identity_3d)))) then
1642 0 : if (err >= DTOL .or. abs(abs(phase_dyn) - one) > DTOL .or. phase_err > DTOL) ierr = ierr + 1
1643 : end if
1644 : call sym_dicts(isym_cnt)%set("inv_ok", &
1645 0 : s=yesno(err < DTOL .and. abs(abs(phase_dyn) - one) <= DTOL .and. phase_err <= DTOL))
1646 0 : call sym_dicts(isym_cnt)%set("inv_err", r=err)
1647 0 : call sym_dicts(isym_cnt)%set("inv_phase", s=sjoin(ftoa(real(phase_dyn)), " + i ", ftoa(aimag(phase_dyn))))
1648 0 : call sym_dicts(isym_cnt)%set("inv_phase_analytic_err", r=phase_err)
1649 : end associate
1650 : end if
1651 0 : if (prtvol > 1) call print_arr(units, cmat, max_r=nb, max_c=nb)
1652 : end associate
1653 : end do ! isym
1654 : end do ! itime
1655 :
1656 : ! Group multiplication test, extended to time reversal (itime1, itime2 in {1,2}).
1657 : ! In ABINIT, point-group operations are applied sequentially to coordinates such that
1658 : ! r' = S_1 S_2 r. When generating the representation matrices D(S, k),
1659 : ! this algebraic structure is maintained according to the product rule:
1660 : !
1661 : ! D^{k}(S_1 S_2) = e^{-i k \cdot L} D^{S_2 k}(S_1) D^{k}(S_2)
1662 : !
1663 : ! Since we are operating strictly inside the little group of k, we have S_2 k \equiv k,
1664 : ! and the equation fundamentally simplifies to a proportionality:
1665 : !
1666 : ! D(S_3) = e^{i \phi} D(S_1) D(S_2)
1667 : !
1668 : ! We search for the composite symmetry isym3 that perfectly matches the spatial
1669 : ! rotation product: symrel(isym1) * symrel(isym2). The spatial rotation composition
1670 : ! rule is itime-independent because \hat\Theta commutes with any pure spatial
1671 : ! coordinate transformation acting on the full (not just periodic-part) wavefunction:
1672 : ! \hat\Theta \hat S \psi(r) = [\hat S\psi(r)]^* = \psi(S^{-1}r)^* = \hat S[\hat\Theta\psi](r).
1673 : !
1674 : ! What DOES depend on itime is the Wigner co-representation composition law itself
1675 : ! (Bradley & Cracknell, sec. 7.3): composing two operators A=(isym1,itime1) and
1676 : ! B=(isym2,itime2), with A applied after B,
1677 : !
1678 : ! D(A B) = D(A) D(B) if A is unitary (itime1 == 1)
1679 : ! D(A B) = D(A) D(B)^* if A is antiunitary (itime1 == 2)
1680 : !
1681 : ! and itime3 (unitary/antiunitary character of A B) follows from Theta^2 = +1
1682 : ! for the scalar (nspinor=1) wavefunctions handled here:
1683 : !
1684 : ! itime3 = 1 + mod((itime1-1) + (itime2-1), 2)
1685 : !
1686 : ! i.e. antiunitary o antiunitary = unitary, matching Theta^2=+1 (Kramers-degeneracy
1687 : ! sign would flip this to Theta^2=-1 for spinors, not implemented/tested: dmats_init
1688 : ! hard-requires nspinor=1).
1689 : !
1690 : ! Diagnostics: record every FAILING (isym1,itime1,isym2,itime2,isym3,itime3) tuple
1691 : ! (with its g0's and errors) instead of just incrementing ierr, so a caller like
1692 : ! dmats_check_star (run over the full BZ, where the k passed in need not be a genuine
1693 : ! IBZ point) can pinpoint exactly which composition and which umklapp broke.
1694 0 : mult_fail_cnt = 0
1695 0 : ABI_MALLOC(mult_fail_dicts, (4 * dmats%cryst%nsym**2))
1696 0 : do itime1=1,2
1697 0 : do itime2=1,2
1698 0 : itime3 = 1 + mod((itime1 - 1) + (itime2 - 1), 2)
1699 0 : do isym1=1,dmats%cryst%nsym
1700 0 : if (symtab(4, itime1, isym1) == 0) cycle
1701 0 : do isym2=1,dmats%cryst%nsym
1702 0 : if (symtab(4, itime2, isym2) == 0) cycle
1703 :
1704 0 : isym3 = dmats%multable(1, isym1, isym2)
1705 :
1706 0 : if (isym3 /= 0 .and. symtab(4, itime3, isym3) /= 0) then
1707 : associate (cmat1 => dmat_k(:, :, isym1, itime1), &
1708 : cmat2 => dmat_k(:, :, isym2, itime2), &
1709 : cmat3 => dmat_k(:, :, isym3, itime3))
1710 :
1711 : ! Instead of failing the test due to phase formula mismatch, we can just EXTRACT the phase!
1712 : ! ABINIT's exact phase might have extra factors due to how istwf_k and cgtk_rotate conjugate things.
1713 : ! The goal is to check if they are proportional (i.e. group structure is satisfied up to a phase).
1714 0 : if (itime1 == 1) then
1715 0 : phase_L = sum( conjg(cmat3) * matmul(cmat1, cmat2) ) / nb
1716 0 : err = maxval(abs(cmat3 * (phase_L / abs(phase_L)) - matmul(cmat1, cmat2)))
1717 : else
1718 0 : phase_L = sum( conjg(cmat3) * matmul(cmat1, conjg(cmat2)) ) / nb
1719 0 : err = maxval(abs(cmat3 * (phase_L / abs(phase_L)) - matmul(cmat1, conjg(cmat2))))
1720 : end if
1721 :
1722 : ! Analytic prediction of the same phase, from the "Caveat for tabulated symmetry
1723 : ! matrices" in main.tex: the literal Seitz product S1S2 and the tabulated operation
1724 : ! S3=isym3 sharing its rotation differ by a pure lattice translation L, giving
1725 : ! D^k(S3) = e^{-i (S3 k).L} D(S1) D(S2), where (S3 k) is the PURE spatial rotation
1726 : ! of S3 applied to k (symrel^t, no time-reversal sign: L comes from the translation
1727 : ! part of the spatial space group only, unrelated to Theta).
1728 : ! dmats%multable/toinv are built from the plain {symrel,tnons} Seitz convention,
1729 : ! while dmats%for_spin is indexed with the symrel^t convention used throughout this
1730 : ! file for the k-action of a symmetry. Reconciling the two requires L to be looked
1731 : ! up at the GROUP-THEORETIC INVERSES of isym1 and isym2 (in reversed order):
1732 : ! L = multable(2:4, toinv(isym2), toinv(isym1)), not multable(2:4,isym1,isym2).
1733 : ! With this, phase_analytic_mult matches phase_L exactly for all 768
1734 : ! (isym1,isym2,itime1,itime2) tuples tested on the reference gstore test.
1735 0 : isym1_inv = dmats%toinv(1, isym1)
1736 0 : isym2_inv = dmats%toinv(1, isym2)
1737 0 : Sk3 = matmul(transpose(real(dmats%cryst%symrel(:,:,isym3), dp)), kk_ibz)
1738 0 : L_mult = real(dmats%multable(2:4, isym2_inv, isym1_inv), dp)
1739 0 : phase_analytic_mult = exp(cmplx(zero, -two_pi * sum(Sk3 * L_mult), dp))
1740 : ! VALIDATED (session diagnostic, classify.abi+nband 18 on gstore_lead): for nspinor==2,
1741 : ! this formula is missing a UNIVERSAL, unconditional extra factor of -1 -- confirmed
1742 : ! across all 4760 (isym1,isym2,itime1,itime2) tuples tested (all 4 itime1/itime2
1743 : ! combinations), zero exceptions. Unlike the inv_ok/closure_ok fixes above, this one
1744 : ! does NOT depend on involutory character or itime at all.
1745 0 : if (dmats%dtset%nspinor == 2) phase_analytic_mult = -phase_analytic_mult
1746 0 : phase_err_mult = abs(phase_L * phase_analytic_mult - one)
1747 :
1748 : ! NOTE on phase_err_mult and dmat_star (reconstructed full-BZ D-matrices, see
1749 : ! dmats_check_star/dmats_get_star_dmats): instrumented this test (temporarily) to
1750 : ! record every failing tuple and confirmed, on the k' points where check_star
1751 : ! reports failures, that ALL of them have err and |phase_L|-1 at machine precision
1752 : ! (true proportionality/closure holds EXACTLY) while phase_err_mult is exactly 2.0
1753 : ! (a clean sign flip, not noise) for every single one -- i.e. this is the SAME
1754 : ! class of "consistently exactly wrong by a clean phase factor" issue already
1755 : ! flagged as diagnostic-only, unresolved, for IMPROPER operations in the S^n
1756 : ! closure test below and for the inverse-relation test above. Tried the natural
1757 : ! alternative convention (L = multable(2:4,isym1,isym2) directly, no toinv-reversal,
1758 : ! dotted with kk_ibz instead of Sk3 -- provably equivalent to Sk3 since g0.L_mult is
1759 : ! always an integer): it does NOT universally fix it either (worse overall, and the
1760 : ! two conventions disagree on non-overlapping subsets of tuples), so this isn't a
1761 : ! simple sign/convention swap in phase_analytic_mult -- the true fix requires
1762 : ! working out how phase_h (dmats_get_star_dmats's own per-isym reconstruction
1763 : ! phase) interacts with the k'-frame tabulated-vs-literal correction L_mult, which
1764 : ! is not yet derived. Until then, gate ierr on the two properties that constitute
1765 : ! actual group-representation closure (proportionality + unit modulus), matching
1766 : ! the precedent set by the two other diagnostic-only checks in this routine, and
1767 : ! keep phase_err_mult as a reported (not gating) diagnostic.
1768 0 : if (err >= DTOL .or. abs(abs(phase_L) - one) > DTOL) then
1769 0 : ierr = ierr + 1
1770 0 : mult_fail_cnt = mult_fail_cnt + 1
1771 0 : call mult_fail_dicts(mult_fail_cnt)%set("isym1", i=isym1)
1772 0 : call mult_fail_dicts(mult_fail_cnt)%set("itime1", i=itime1)
1773 0 : call mult_fail_dicts(mult_fail_cnt)%set("isym2", i=isym2)
1774 0 : call mult_fail_dicts(mult_fail_cnt)%set("itime2", i=itime2)
1775 0 : call mult_fail_dicts(mult_fail_cnt)%set("isym3", i=isym3)
1776 0 : call mult_fail_dicts(mult_fail_cnt)%set("itime3", i=itime3)
1777 0 : call mult_fail_dicts(mult_fail_cnt)%set("g0_1", s=trim(ltoa(symtab(1:3, itime1, isym1))))
1778 0 : call mult_fail_dicts(mult_fail_cnt)%set("g0_2", s=trim(ltoa(symtab(1:3, itime2, isym2))))
1779 0 : call mult_fail_dicts(mult_fail_cnt)%set("g0_3", s=trim(ltoa(symtab(1:3, itime3, isym3))))
1780 0 : call mult_fail_dicts(mult_fail_cnt)%set("err", r=err)
1781 0 : call mult_fail_dicts(mult_fail_cnt)%set("phase_mod_err", r=abs(abs(phase_L) - one))
1782 0 : call mult_fail_dicts(mult_fail_cnt)%set("phase_err_mult", r=phase_err_mult)
1783 : end if
1784 : end associate
1785 : end if
1786 : end do
1787 : end do
1788 : end do
1789 : end do
1790 :
1791 0 : if (mult_fail_cnt > 0) then
1792 0 : call ydoc%add_dictlist(sjoin(tag, "_group_mult_fail"), mult_fail_cnt, mult_fail_dicts(1:mult_fail_cnt))
1793 0 : do j = 1, mult_fail_cnt
1794 0 : call mult_fail_dicts(j)%free()
1795 : end do
1796 : end if
1797 0 : ABI_FREE(mult_fail_dicts)
1798 :
1799 : ! =========================================================================
1800 : ! Eigenvalues & Closure Test (itime = 1)
1801 : ! =========================================================================
1802 : ! According to the theory of group representations, the representation matrix
1803 : ! D(S) must satisfy the closure conditions of the crystallographic point group.
1804 : ! If S = C_n is an n-fold symmetry operation, applying the spatial rotation
1805 : ! n times yields the identity (R^n = E).
1806 : ! However, for non-symmorphic operations (e.g. glide planes or screw axes),
1807 : ! applying the operation n times results in a pure fractional lattice translation:
1808 : ! S^n(r) = r + T
1809 : !
1810 : ! In reciprocal space, inside the little group of k, this translation introduces
1811 : ! a scalar Bloch phase shift. Thus, the eigenvalues of the representation matrix satisfy:
1812 : !
1813 : ! [ D(S) ]^n = e^{-i k \cdot T} I
1814 : !
1815 : ! (dmats hard-requires nspinor=1, see dmats_init NOTES, so there's no extra spinor parity).
1816 : ! We extract the overall scalar phase \phi = Tr(D^n) / N_{bands}, assert that
1817 : ! D(S)^n \equiv \phi I, and cross-check \phi against the analytic e^{-i k.T} computed
1818 : ! from the T returned by sym_order.
1819 0 : isym_cnt = 0
1820 0 : do itime=1,2
1821 0 : do isym=1,dmats%cryst%nsym
1822 0 : if (symtab(4, itime, isym) == 0) cycle
1823 0 : isym_cnt = isym_cnt + 1
1824 0 : if (itime == 1) then
1825 : associate (cmat => dmat_k(:, :, isym, 1))
1826 :
1827 : ! Find the order of the point-group operation (n in {1,2,3,4,6}).
1828 : ! NB: use a dedicated ierr_so for sym_order's own status -- passing the shared
1829 : ! accumulator "ierr" directly would have sym_order's intent(out) silently reset it
1830 : ! to 0 on every call, wiping out all previously-accumulated test failures.
1831 0 : call sym_order(dmats%cryst%symrel(:,:,isym), dmats%cryst%tnons(:,isym), n, isproper, trans, msg, ierr_so)
1832 0 : ABI_CHECK_IEQ(ierr_so, 0, msg)
1833 :
1834 0 : if (n > 1) then
1835 : ! Compute cmat^n
1836 0 : cmat_n = cmat
1837 0 : do j = 2, n
1838 0 : cmat_n = matmul(cmat, cmat_n)
1839 : end do
1840 :
1841 : ! Extract phase from Trace: phase = Tr(cmat^n) / nb
1842 0 : phase_L = zero
1843 0 : do j = 1, nb
1844 0 : phase_L = phase_L + cmat_n(j, j)
1845 : end do
1846 0 : phase_L = phase_L / nb
1847 :
1848 : ! Normalize cmat_n with phase_L to check if it's proportional to identity
1849 : err = zero
1850 0 : do j = 1, nb
1851 0 : cmat_n(j, j) = cmat_n(j, j) - phase_L
1852 : end do
1853 0 : err = maxval(abs(cmat_n))
1854 :
1855 : ! Analytic Bloch phase from the cumulative lattice translation T: D(S)^n = e^{-i k.T} I
1856 0 : phase_analytic = exp(cmplx(zero, -two_pi * dot_product(kk_ibz, real(trans, dp)), dp))
1857 : ! TEMPORARY DIAGNOSTIC (session experiment): applying S n times, for S a proper rotation
1858 : ! of order n, sweeps EXACTLY one full 2*pi turn -- for spinors D(2*pi) = -I always
1859 : ! (unconditionally, not a convention/branch-choice issue), so D(S)^n should pick up an
1860 : ! extra factor of -1 relative to the scalar-derived formula above, for EVERY proper
1861 : ! rotation (any n), regardless of itime (closure test is itime=1 only anyway).
1862 0 : if (dmats%dtset%nspinor == 2 .and. isproper) phase_analytic = -phase_analytic
1863 0 : phase_err = abs(phase_L - phase_analytic)
1864 :
1865 : ! NOTE: for IMPROPER operations (isproper=.false.), phase_err is consistently
1866 : ! found to be exactly 2 (phase_L = -phase_analytic, a clean sign flip, not noise)
1867 : ! -- the same class of unresolved cgtk_rotate phase-convention ambiguity already
1868 : ! flagged (but left diagnostic-only, not gating ierr) in the inverse-relation test
1869 : ! above ("still under investigation"). Follow that same precedent here: gate ierr
1870 : ! on proportionality (err) and unit modulus, but not on phase_err for improper ops.
1871 0 : if (err >= DTOL .or. abs(abs(phase_L) - one) > DTOL .or. (isproper .and. phase_err > DTOL)) ierr = ierr + 1
1872 : call sym_dicts(isym_cnt)%set("closure_ok", &
1873 0 : s=yesno(err < DTOL .and. abs(abs(phase_L) - one) <= DTOL .and. (.not. isproper .or. phase_err <= DTOL)))
1874 0 : call sym_dicts(isym_cnt)%set("closure_err", r=err)
1875 0 : call sym_dicts(isym_cnt)%set("closure_n", i=n)
1876 0 : call sym_dicts(isym_cnt)%set("isproper", s=yesno(isproper))
1877 0 : call sym_dicts(isym_cnt)%set("closure_phase", s=sjoin(ftoa(real(phase_L)), " + i ", ftoa(aimag(phase_L))))
1878 : call sym_dicts(isym_cnt)%set("closure_phase_analytic", &
1879 0 : s=sjoin(ftoa(real(phase_analytic)), " + i ", ftoa(aimag(phase_analytic))))
1880 0 : call sym_dicts(isym_cnt)%set("closure_phase_err", r=phase_err)
1881 : end if
1882 : end associate
1883 : end if
1884 : end do
1885 : end do
1886 :
1887 0 : if (allocated(sym_dicts)) then
1888 0 : call ydoc%add_dictlist(tag, isym_cnt, sym_dicts)
1889 0 : do isym = 1, isym_cnt
1890 0 : call sym_dicts(isym)%free()
1891 : end do
1892 0 : ABI_FREE(sym_dicts)
1893 : end if
1894 :
1895 0 : ABI_FREE(cmat_n)
1896 :
1897 0 : end subroutine dmats_check_one_k
1898 : !!***
1899 :
1900 : !----------------------------------------------------------------------
1901 :
1902 : !!****f* m_classify_bands/dmats_check
1903 : !! NAME
1904 : !! dmats_check
1905 : !!
1906 : !! FUNCTION
1907 : !! Verify the fundamental point-group algebraic properties of
1908 : !!
1909 : !! D^{k}_{mn}(S) = < \psi_{m, Sk} | S | \psi_{n, k} >
1910 : !!
1911 : !! The following algebraic tests are performed:
1912 : !!
1913 : !! 1. Unitarity (mandatory): || D^\dagger(k, S) D(k, S) - I || < DTOL
1914 : !! 2. Identity operator: D(E, k) = I, for isym = 1
1915 : !! 3. Inverse relation: D^{Sk}(S^{-1}) \propto D^{k}(S)^\dagger
1916 : !! 4. Group multiplication, including time-reversal: D^{k}(A B) \propto D^{k}(A) D^{k}(B)
1917 : !! for A, B each either a pure spatial symmetry or a spatial symmetry composed with
1918 : !! time reversal (Wigner co-representation composition law, see below).
1919 : !! 5. Kramers/Theta^2 (only at TR-invariant k, e.g. TRIM points): D^{k}(\Theta) D^{k}(\Theta)^*
1920 : !! = +I EXACTLY (Theta^2=+1 for scalar wavefunctions), for isym=1 (pure time reversal), itime=2.
1921 : !! 6. Character class-function test (itime=1 only): |Tr D^{k}(S)| is the same for every S in a
1922 : !! given conjugacy class of the little group of k, as computed by get_classes (m_ptgroups.F90).
1923 : !! |trace| (not the raw complex trace) is used deliberately: it is invariant under both the
1924 : !! tabulated-vs-literal-composition phase ambiguity (see the note above on test 4) and any
1925 : !! per-band wavefunction gauge choice, since a similarity transform composed with an overall
1926 : !! scalar phase always preserves |trace| exactly, with no analytic phase tracking needed.
1927 : !!
1928 : !! Caution for anyone editing this test: get_classes takes EXPLICIT-SHAPE dummy arguments
1929 : !! sized from its own nsym argument (here, the little-group order nsym_lg, which varies by
1930 : !! k-point and is in general smaller than cryst%nsym). Passing whole arrays declared with a
1931 : !! LARGER fixed bound (e.g. sized to cryst%nsym) works via Fortran sequence association but
1932 : !! silently corrupts the result for any 2D array (elements_idx) whenever nsym_lg != cryst%nsym,
1933 : !! because the callee writes using an nsym_lg-based column-major stride while the caller would
1934 : !! read back using a cryst%nsym-based stride. (1D arrays, like nelements, are unaffected: their
1935 : !! offset does not depend on the declared bound.) The fix is to allocate elements_idx to EXACTLY
1936 : !! (nsym_lg, nsym_lg), as done here and in esymm_init (m_esymm.F90). This is not a hypothetical
1937 : !! concern: it was hit and diagnosed during development of this test, and produced clean-looking
1938 : !! but WRONG class assignments (e.g. an isym reported as belonging to two different classes
1939 : !! depending on which entry point read it) rather than an obvious crash.
1940 : !!
1941 : !! Note on algebraic structure:
1942 : !! ABINIT's symmetries correspond to a *right homomorphism*, meaning the product
1943 : !! of two symmetries S_3 = S_1 S_2 (where the rotation parts are R_3 = R_1 R_2)
1944 : !! yields representations that compose as D(S_3) \propto D(S_1) D(S_2).
1945 : !!
1946 : !! For fractional translations in non-symmorphic groups or due to G_0 vector mappings,
1947 : !! the exact analytical phases between operations can become extremely complex.
1948 : !! Therefore, the inverse relation and group multiplication tests extract the relative
1949 : !! phase dynamically using the Frobenius inner product Phase = Tr(A^\dagger B) / nb.
1950 : !! As long as the residual after phase-normalization is within DTOL, the matrices
1951 : !! strictly satisfy the projective representations of the space group.
1952 : !!
1953 : !! Diagnostic value of test 4 (group multiplication): unitarity (test 1) and the
1954 : !! inverse relation (test 3) are, by construction, blind to a bug in which every
1955 : !! D(S) is silently replaced by D(S)^{-1} = D(S^{-1}): both tests only ever compare
1956 : !! a matrix against itself or its own inverse, so swapping S <-> S^{-1} consistently
1957 : !! leaves them satisfied. The group multiplication test does NOT have this blind
1958 : !! spot for non-involutory S (S^2 != E), since D(S)<->D(S)^{-1} breaks the
1959 : !! non-commutative composition D(S_1 S_2) \propto D(S_1) D(S_2) as soon as one of the
1960 : !! three operations involved has order > 2 (D(S_3) becomes proportional to
1961 : !! D(S_1) D(S_2)^t instead). This is exactly how a real S<->S^{-1} mislabeling bug in
1962 : !! dmats_init was caught during development; see the NOTES section of dmats_init's
1963 : !! SOURCE header for the full explanation and fix.
1964 : !!
1965 : !! Test 4 loops over all itime1, itime2 in {1,2}, not just the pure-spatial (1,1)
1966 : !! case: this closes an analogous blind spot for time-reversal-related matrices
1967 : !! D(k, S, itime=2), which unitarity/inverse-relation alone cannot detect either.
1968 : !! The composite operation (isym3, itime3) is found from: (a) the ROTATION part,
1969 : !! which composes as symrel(isym1).symrel(isym2) regardless of itime1/itime2,
1970 : !! because \hat\Theta commutes with any pure spatial coordinate transformation
1971 : !! acting on the full (not just periodic-part) wavefunction; and (b) the
1972 : !! unitary/antiunitary character, itime3 = 1 + mod((itime1-1)+(itime2-1), 2),
1973 : !! i.e. antiunitary o antiunitary = unitary (valid for scalar wavefunctions,
1974 : !! Theta^2=+1; dmats_init hard-requires nspinor=1). The matrix relation itself
1975 : !! follows the Wigner co-representation composition law (Bradley & Cracknell,
1976 : !! sec. 7.3): D(A B) = D(A) D(B) if A is unitary, D(A B) = D(A) D(B)^* if A is
1977 : !! antiunitary (A applied after B, A=(isym1,itime1)).
1978 : !!
1979 : !! SOURCE
1980 :
1981 0 : subroutine dmats_check(dmats, units, prtvol, header)
1982 :
1983 : !Arguments ------------------------------------
1984 : class(dmats_t),intent(in) :: dmats
1985 : integer,intent(in) :: units(:), prtvol
1986 : character(len=*),optional,intent(in) :: header
1987 :
1988 : !Local variables-------------------------------
1989 : integer :: spin, ik_ibz, ierr
1990 : character(len=500) :: msg
1991 : type(yamldoc_t) :: ydoc
1992 : ! *************************************************************************
1993 :
1994 0 : msg = 'Info on the dmats_t'
1995 0 : if (present(header)) msg = trim(adjustl(header))
1996 0 : ydoc = yamldoc_open(tag="dmats", info=trim(msg))
1997 :
1998 0 : ierr = 0
1999 0 : do spin=1,size(dmats%for_spin)
2000 : ! Loop over k-points in the IBZ.
2001 0 : do ik_ibz=1,dmats%ks_ebands%nkpt
2002 : call dmats_check_one_k(dmats, spin, dmats%ks_ebands%kptns(:, ik_ibz), &
2003 : dmats%for_spin(spin)%value(:, :, :, :, ik_ibz), units, prtvol, &
2004 : sjoin("kpt_", ktoa(dmats%ks_ebands%kptns(:, ik_ibz)), "_spin_", itoa(spin)), &
2005 0 : ydoc, ierr)
2006 : end do ! ik_ibz
2007 : end do ! spin
2008 :
2009 0 : call ydoc%write_units_and_free(units)
2010 :
2011 0 : if (ierr /= 0) then
2012 0 : ABI_ERROR(sjoin("dmats are not unitary or failed tests! ierr:", itoa(ierr)))
2013 : end if
2014 :
2015 0 : end subroutine dmats_check
2016 : !!***
2017 :
2018 : !!****f* m_classify_bands/dmats_get_star_dmats
2019 : !! NAME
2020 : !! dmats_get_star_dmats
2021 : !!
2022 : !! FUNCTION
2023 : !! Build the D-matrices at a full-BZ k-point k' = S0.k_ibz (S0 = (isym0,itime0)), a
2024 : !! symmetry-star image of an already-computed IBZ k-point, using ONLY data already
2025 : !! stored in dmats (for_spin, multable, toinv, cryst) -- no WFK re-read (wfd is freed
2026 : !! by dmats_init once dmats%for_spin has been built).
2027 : !!
2028 : !! Define the wavefunction gauge at k' as |n,k'> := S0|n,k_ibz>. Then for any
2029 : !! g=(isym,itime) that stabilizes k' (g.k' = k' mod G), with h := S0^{-1}.g.S0
2030 : !! (physical operator composition):
2031 : !!
2032 : !! D^{k'}_{mn}(g) = <m,k'|g|n,k'> = <m,k_ibz|S0^{-1} g S0|n,k_ibz> = D^{k_ibz}_{mn}(h)
2033 : !!
2034 : !! exactly. IMPORTANT: resolving this h to a tabulated isym index via dmats%multable is
2035 : !! NOT simply "two multable products in the S0^{-1},g,S0 order": the file's k-vector
2036 : !! little-group test phi(s):=symrel(s)^t is an ANTI-homomorphism of multable's abstract
2037 : !! (plain, non-transposed) group law -- phi(s1 applied after s2) = phi(s2).phi(s1), order
2038 : !! REVERSED (same anti-homomorphism already flagged in dmats_check_one_k's NOTES on
2039 : !! get_classes). Working through phi(h)=phi(S0)^{-1}.phi(g).phi(S0) with this reversal
2040 : !! shows the correct tabulated composition is actually h = S0.g.S0^{-1} (see the detailed
2041 : !! derivation in this routine's SOURCE, right before the two multable calls).
2042 : !!
2043 : !! Since dmats%multable(1,...) gives the "tabulated" symrel entry sharing h's rotation,
2044 : !! not the literal Seitz product, the literal h differs from it by a lattice vector L_h:
2045 : !!
2046 : !! D^{k'}(g) = e^{-i 2pi k_ibz.L_h} * dmats%for_spin(spin)%value(:,:,isym_h,itime_h,ik_ibz)
2047 : !!
2048 : !! L_h is accumulated through TWO nested multable compositions (h is itself the double
2049 : !! product S0.(g.S0^{-1})), using the same reversed-inverse-argument lookup already
2050 : !! validated in dmats_check_one_k's group-multiplication test (L = multable(2:4,
2051 : !! toinv(isym2), toinv(isym1)) for a product D(S1 S2), S1 applied after S2).
2052 : !!
2053 : !! PROOF of L_h (worked from scratch, tabulated-operator algebra only): write
2054 : !! Sigma0^{-1} = Tab(isym0_inv) o T_{-L0} (L0 = toinv(2:4,isym0), from S0.Tab(isym0_inv) =
2055 : !! T_{L0}), and for ANY two tabulated elements A=Tab(a), B=Tab(b): A o B = T_{L(a,b)} o
2056 : !! Tab(multable(1,a,b)), L(a,b) = multable(2:4,a,b) -- the lattice correction sits on the
2057 : !! LEFT of the tabulated product, so it picks up a rotation (T_v o Tab(c) = Tab(c) o
2058 : !! T_{R(c)^{-1}.v}, equivalently R(c).T_v o Tab(c)... ) whenever it is pushed further left
2059 : !! past another rotation. Substituting Sigma0^{-1} into h_lit = S0.g.Sigma0^{-1} and
2060 : !! applying this rule twice (first at g.Tab(isym0_inv), then at S0.Tab(isym_tmp)) gives,
2061 : !! with NO free parameters left over:
2062 : !!
2063 : !! isym_tmp = multable(1, isym_g, isym0_inv); Ltmp = multable(2:4, isym_g, isym0_inv)
2064 : !! isym_h = multable(1, isym0, isym_tmp); Lh2 = multable(2:4, isym0, isym_tmp)
2065 : !! L_h = R0.Ltmp + Lh2 - R_h.L0 (R0 = symrel(:,:,isym0), R_h = symrel(:,:,isym_h))
2066 : !!
2067 : !! which matches the L_h computed in this routine's SOURCE exactly. This derivation used
2068 : !! ONLY real-space Seitz-operator algebra (unitarity of S0 plus associativity of operator
2069 : !! composition) -- no reciprocal-space bookkeeping is needed anywhere, so a genuinely
2070 : !! non-zero umklapp vector g0 in the little-group-of-k' test just below does NOT introduce
2071 : !! any extra phase on top of L_h: Bloch periodicity psi_{k+G} = psi_k is an exact identity,
2072 : !! not a gauge choice, so it never needed to be invoked in this chain. The one-step version
2073 : !! of this phase formula was independently, empirically validated (exact match in the
2074 : !! group-multiplication test on 768 tuples); this two-step formula is now ALSO proven, not
2075 : !! just extrapolated. If dmats_check_star still reports group-multiplication/closure
2076 : !! failures correlated with non-zero umklapp, the bug is therefore NOT in this phase
2077 : !! formula -- look instead at the (isym0,itime0) selection in
2078 : !! dmats_get_star_dmats_at_kpt (first-match-wins when k_ibz has a non-trivial little
2079 : !! group) or at the little-group-of-k' membership test just below.
2080 : !!
2081 : !! INPUTS
2082 : !! spin=Spin index.
2083 : !! ik_ibz=Index of the reference IBZ k-point (into dmats%ks_ebands%kptns).
2084 : !! isym0, itime0=Symmetry S0=(isym0,itime0) such that k' = S0.k_ibz.
2085 : !!
2086 : !! OUTPUT
2087 : !! dmat_star(:,:,:,:)=D-matrices (nb,nb,nsym,2) at k'. Caller-owned: allocated here
2088 : !! (ABI_MALLOC), must be freed by the caller (ABI_FREE).
2089 : !! msg, ierr=Error message and status (ierr=0 on success). No ABI_ERROR is raised here:
2090 : !! this routine is meant to be usable in a validation loop over many k', where the
2091 : !! caller decides whether to abort (e.g. via ABI_CHECK_IEQ(ierr, 0, msg)) or skip.
2092 : !!
2093 : !! SOURCE
2094 :
2095 0 : subroutine dmats_get_star_dmats(dmats, spin, ik_ibz, isym0, itime0, dmat_star, msg, ierr)
2096 :
2097 : !Arguments ------------------------------------
2098 : class(dmats_t),intent(in) :: dmats
2099 : integer,intent(in) :: spin, ik_ibz, isym0, itime0
2100 : complex(dp),allocatable,intent(out) :: dmat_star(:,:,:,:)
2101 : character(len=*),intent(out) :: msg
2102 : integer,intent(out) :: ierr
2103 :
2104 : !Local variables-------------------------------
2105 : integer :: nsym, nb, isym, itime, isym_tmp, isym_h, itime_tmp, itime_h, isym0_inv, j
2106 : integer :: g0(3), g0_h(3)
2107 : real(dp) :: kk_ibz(3), kprime(3), kk_sk(3), tsign0, tsign
2108 : real(dp) :: L_h(3)
2109 : complex(dp) :: phase_h
2110 : ! *********************************************************************
2111 :
2112 0 : ierr = 0; msg = ""
2113 0 : nsym = dmats%cryst%nsym
2114 0 : nb = dmats%brange_spin(2, spin) - dmats%brange_spin(1, spin) + 1
2115 0 : kk_ibz = dmats%ks_ebands%kptns(:, ik_ibz)
2116 :
2117 : ! k' = S0.k_ibz (symrel^t convention, consistent with dmats_init/dmats_check_one_k).
2118 0 : tsign0 = merge(one, -one, itime0 == 1)
2119 0 : kprime = tsign0 * matmul(transpose(real(dmats%cryst%symrel(:,:,isym0), dp)), kk_ibz)
2120 :
2121 0 : isym0_inv = dmats%toinv(1, isym0)
2122 0 : if (isym0_inv == 0) then
2123 0 : ierr = 1; msg = "Could not find inverse of isym0"; return
2124 : end if
2125 :
2126 0 : ABI_MALLOC(dmat_star, (nb, nb, nsym, 2))
2127 0 : dmat_star = czero
2128 :
2129 0 : do itime=1,2
2130 0 : tsign = merge(one, -one, itime == 1)
2131 0 : do isym=1,nsym
2132 :
2133 : ! Is g=(isym,itime) in the little group of k' = S0.k_ibz?
2134 0 : kk_sk = tsign * matmul(transpose(real(dmats%cryst%symrel(:,:,isym), dp)), kprime)
2135 0 : g0 = nint(kprime - kk_sk)
2136 0 : if (.not. all(abs(kprime - kk_sk - g0) < tol8)) then
2137 : ! Not in little group of k': identity placeholder, mirrors dmats%for_spin's own convention.
2138 0 : do j=1,nb
2139 0 : dmat_star(j, j, isym, itime) = cone
2140 : end do
2141 0 : cycle
2142 : end if
2143 :
2144 : ! Compose h = S0.g.S0^{-1} via TWO multable compositions: multable(1,s1,s2) = index of
2145 : ! "s1 applied after s2" (see sg_multable), i.e. the plain, NON-transposed real-space
2146 : ! rotation-matrix product R(s1).R(s2).
2147 : !
2148 : ! NOTE the conjugation direction: naively one would expect h = S0^{-1}.g.S0 (as in an
2149 : ! ordinary homomorphism), but the file's k-vector action phi(s) := symrel(s)^t is an
2150 : ! ANTI-homomorphism of the abstract (multable) group law: phi(s1 "applied after" s2) =
2151 : ! R(s1.s2)^t = R(s2)^t.R(s1)^t = phi(s2).phi(s1) -- composition order REVERSES (this is
2152 : ! the same anti-homomorphism already noted in dmats_check_one_k's NOTES on get_classes).
2153 : ! Requiring phi(h) = phi(S0)^{-1}.phi(g).phi(S0) (so that h stabilizes k_ibz whenever g
2154 : ! stabilizes k'=phi(S0).k_ibz) and using phi(A)phi(B)=phi(B.A) twice gives
2155 : ! phi(h) = phi(S0.g.S0^{-1}), i.e. h = S0.g.S0^{-1}, NOT S0^{-1}.g.S0.
2156 0 : isym_tmp = dmats%multable(1, isym, isym0_inv) ! tmp = g . S0^{-1}
2157 0 : if (isym_tmp == 0) then
2158 0 : ierr = 2; msg = "multable(isym, isym0_inv) not found: group closure violated?"; return
2159 : end if
2160 0 : itime_tmp = 1 + mod((itime - 1) + (itime0 - 1), 2)
2161 :
2162 0 : isym_h = dmats%multable(1, isym0, isym_tmp) ! h = S0 . tmp = S0.g.S0^{-1}
2163 0 : if (isym_h == 0) then
2164 0 : ierr = 2; msg = "multable(isym0, isym_tmp) not found: group closure violated?"; return
2165 : end if
2166 0 : itime_h = 1 + mod((itime0 - 1) + (itime_tmp - 1), 2) ! always equals itime (parity self-cancels)
2167 0 : if (itime_h /= itime) then
2168 0 : ierr = 2; msg = "itime_h != itime_g: parity composition bug"; return
2169 : end if
2170 :
2171 : ! Defensive check: h must stabilize k_ibz by construction.
2172 0 : kk_sk = tsign * matmul(transpose(real(dmats%cryst%symrel(:,:,isym_h), dp)), kk_ibz)
2173 0 : g0_h = nint(kk_ibz - kk_sk)
2174 0 : if (.not. all(abs(kk_ibz - kk_sk - g0_h) < tol8)) then
2175 0 : ierr = 2; msg = "h does not stabilize k_ibz: composition bug"; return
2176 : end if
2177 :
2178 : ! Lattice-vector correction L_h, derived directly from Seitz algebra (verified against
2179 : ! the g=identity special case, where it must vanish exactly -- D(identity) = I with no
2180 : ! phase, always). Writing S0={R0,tau0}, g={Rg,taug}, and S0_inv_tab=dmats%toinv's TABULATED
2181 : ! entry for S0^{-1} (which equals the EXACT inverse only up to an extra lattice shift
2182 : ! m0 = R0^{-1}.L0, L0=toinv(2:4,isym0), since toinv only guarantees S0.S0_inv_tab={I,L0}):
2183 : !
2184 : ! tmp_literal := g . S0_inv_tab = {I, Ltmp} . TABULATED_tmp, Ltmp = multable(2:4,isym,isym0_inv)
2185 : ! h_tab_literal := S0 . tmp_literal = {I, R0.Ltmp + Lh2} . TABULATED_h, Lh2 = multable(2:4,isym0,isym_tmp)
2186 : !
2187 : ! h_tab_literal uses S0_inv_tab, not the EXACT inverse \hat S0^{-1} = {I,-m0}.S0_inv_tab; undoing
2188 : ! that extra {I,m0} shift (tracked through the same two compositions) gives the additional
2189 : ! correction -R_h.L0 (R_h=symrel(isym_h)), so that for g=identity (Ltmp=Lh2=L0, R_h=I) the
2190 : ! total L_h = R0.0 + L0 - I.L0 = 0 exactly, as required:
2191 : !
2192 : ! L_h = R0.Ltmp + Lh2 - R_h.L0
2193 : !
2194 : ! The resulting {I,L_h} pure-lattice-translation factor is applied AFTER TABULATED_h (which
2195 : ! stabilizes k_ibz), so the state is still at k_ibz when the translation phase is picked up:
2196 : ! phase_h = e^{-i 2pi k_ibz.L_h} (dot directly with k_ibz, not with a rotated k_ibz).
2197 : L_h = matmul(real(dmats%cryst%symrel(:,:,isym0), dp), real(dmats%multable(2:4, isym, isym0_inv), dp)) &
2198 : + real(dmats%multable(2:4, isym0, isym_tmp), dp) &
2199 0 : - matmul(real(dmats%cryst%symrel(:,:,isym_h), dp), real(dmats%toinv(2:4, isym0), dp))
2200 :
2201 0 : phase_h = exp(cmplx(zero, -two_pi * sum(kk_ibz * L_h), dp))
2202 0 : if (itime == 2) phase_h = conjg(phase_h)
2203 :
2204 0 : dmat_star(:, :, isym, itime) = phase_h * dmats%for_spin(spin)%value(:, :, isym_h, itime, ik_ibz)
2205 : end do
2206 : end do
2207 :
2208 0 : end subroutine dmats_get_star_dmats
2209 : !!***
2210 :
2211 : !!****f* m_classify_bands/dmats_get_star_dmats_at_kpt
2212 : !! NAME
2213 : !! dmats_get_star_dmats_at_kpt
2214 : !!
2215 : !! FUNCTION
2216 : !! Convenience wrapper around dmats_get_star_dmats: given a raw full-BZ k-point kprime,
2217 : !! locate (ik_ibz, isym0, itime0) such that kprime = tsign0*symrel(isym0)^t.k_ibz (mod G),
2218 : !! scanning ALL IBZ k-points, then build the D-matrices at kprime.
2219 : !!
2220 : !! Uses the SAME symrel^t convention as dmats_init/dmats_check_one_k (NOT kpts_map/listkk/
2221 : !! littlegroup_q, which use a symrec-based convention that can disagree with symrel^t
2222 : !! whenever symrel is not orthogonal in the reduced basis -- see the NOTES in
2223 : !! dmats_check_one_k on this exact point).
2224 : !!
2225 : !! INPUTS
2226 : !! spin=Spin index.
2227 : !! kprime(3)=Target k-point (reduced coords) in the full BZ.
2228 : !!
2229 : !! OUTPUT
2230 : !! dmat_star(:,:,:,:)=D-matrices (nb,nb,nsym,2) at kprime. Caller-owned (ABI_MALLOC/ABI_FREE).
2231 : !! ik_ibz, isym0, itime0=The located triple, returned so the caller can report/reuse it.
2232 : !! msg, ierr=Error message and status (ierr=0 on success, /=0 if kprime is not the star
2233 : !! image of any IBZ k-point in dmats%ks_ebands).
2234 : !!
2235 : !! SOURCE
2236 :
2237 0 : subroutine dmats_get_star_dmats_at_kpt(dmats, spin, kprime, dmat_star, ik_ibz, isym0, itime0, msg, ierr)
2238 :
2239 : !Arguments ------------------------------------
2240 : class(dmats_t),intent(in) :: dmats
2241 : integer,intent(in) :: spin
2242 : real(dp),intent(in) :: kprime(3)
2243 : complex(dp),allocatable,intent(out) :: dmat_star(:,:,:,:)
2244 : integer,intent(out) :: ik_ibz, isym0, itime0
2245 : character(len=*),intent(out) :: msg
2246 : integer,intent(out) :: ierr
2247 :
2248 : !Local variables-------------------------------
2249 : integer :: jk_ibz, jsym, jtime
2250 : real(dp) :: kk_ibz(3), kk_sk(3), tsign, g0(3), resid, best_resid
2251 : logical :: found
2252 : ! *********************************************************************
2253 :
2254 0 : ierr = 0; msg = ""; found = .False.
2255 0 : ik_ibz = -1; isym0 = -1; itime0 = -1
2256 0 : best_resid = huge(one)
2257 :
2258 0 : search: do jk_ibz=1,dmats%ks_ebands%nkpt
2259 0 : kk_ibz = dmats%ks_ebands%kptns(:, jk_ibz)
2260 0 : do jtime=1,2
2261 0 : tsign = merge(one, -one, jtime == 1)
2262 0 : do jsym=1,dmats%cryst%nsym
2263 0 : kk_sk = tsign * matmul(transpose(real(dmats%cryst%symrel(:,:,jsym), dp)), kk_ibz)
2264 0 : g0 = nint(kprime - kk_sk)
2265 0 : resid = maxval(abs(kprime - kk_sk - g0))
2266 0 : best_resid = min(best_resid, resid)
2267 0 : if (all(abs(kprime - kk_sk - g0) < tol8)) then
2268 0 : ik_ibz = jk_ibz; isym0 = jsym; itime0 = jtime; found = .True.
2269 : exit search
2270 : end if
2271 : end do
2272 : end do
2273 : end do search
2274 :
2275 : if (.not. found) then
2276 0 : ierr = 1
2277 : msg = sjoin("kprime:", ktoa(kprime), "is not the symmetry-star image of any IBZ k-point", &
2278 0 : "(best residual found:", ftoa(best_resid), ")")
2279 0 : return
2280 : end if
2281 :
2282 0 : call dmats_get_star_dmats(dmats, spin, ik_ibz, isym0, itime0, dmat_star, msg, ierr)
2283 :
2284 0 : end subroutine dmats_get_star_dmats_at_kpt
2285 : !!***
2286 :
2287 : !!****f* m_classify_bands/dmats_check_star
2288 : !! NAME
2289 : !! dmats_check_star
2290 : !!
2291 : !! FUNCTION
2292 : !! Locate a full-BZ k-point kprime as the symmetry-star image of an IBZ k-point, build
2293 : !! its D-matrices via dmats_get_star_dmats_at_kpt (pure group-theory reconstruction, no
2294 : !! WFK re-read), and run the full dmats_check_one_k test battery (unitarity, identity,
2295 : !! class character, Kramers, inverse relation, group multiplication, S^n closure) on
2296 : !! them. This is an independent test of the multable/toinv/conjugation logic used by
2297 : !! dmats_get_star_dmats, in an off-little-group regime dmats_check never exercises.
2298 : !!
2299 : !! INPUTS
2300 : !! spin=Spin index.
2301 : !! kprime(3)=Target k-point (reduced coords) in the full BZ.
2302 : !! units(:), prtvol=Output units and verbosity.
2303 : !!
2304 : !! OUTPUT
2305 : !! ierr=0 if kprime was located and all sub-tests passed, /=0 otherwise.
2306 : !!
2307 : !! SOURCE
2308 :
2309 0 : subroutine dmats_check_star(dmats, spin, kprime, units, prtvol, ierr)
2310 :
2311 : !Arguments ------------------------------------
2312 : class(dmats_t),intent(in) :: dmats
2313 : integer,intent(in) :: spin, units(:), prtvol
2314 : real(dp),intent(in) :: kprime(3)
2315 : integer,intent(out) :: ierr
2316 :
2317 : !Local variables-------------------------------
2318 : integer :: ik_ibz, isym0, itime0
2319 : character(len=500) :: msg
2320 0 : complex(dp),allocatable :: dmat_star(:,:,:,:)
2321 : type(yamldoc_t) :: ydoc
2322 : ! Cross-check dmat_star (pure group theory) against genuinely WFK-computed D-matrices.
2323 : integer :: nsym, isym, itime, isym1, itime1, n_cross_ok, n_cross_bad, nb
2324 : real(dp) :: tsign, kk_ibz(3), kk_g(3), g0(3), maxdiff
2325 0 : complex(dp),allocatable :: d_composed(:,:)
2326 : ! *********************************************************************
2327 :
2328 : ierr = 0
2329 0 : call dmats_get_star_dmats_at_kpt(dmats, spin, kprime, dmat_star, ik_ibz, isym0, itime0, msg, ierr)
2330 0 : if (ierr /= 0) then
2331 0 : call wrtout(units, sjoin("dmats_check_star: get_star_dmats_at_kpt failed:", msg))
2332 0 : return
2333 : end if
2334 :
2335 : ydoc = yamldoc_open(tag="dmats_star", &
2336 : info=sjoin("Star k-point check: kprime=", ktoa(kprime), ", ik_ibz=", itoa(ik_ibz), &
2337 0 : ", isym0=", itoa(isym0), ", itime0=", itoa(itime0)))
2338 :
2339 : call dmats_check_one_k(dmats, spin, kprime, dmat_star, units, prtvol, &
2340 0 : sjoin("starkpt_", ktoa(kprime), "_spin_", itoa(spin)), ydoc, ierr)
2341 :
2342 0 : call ydoc%write_units_and_free(units)
2343 :
2344 : ! Cross-check: dmat_star (pure group theory, no WFK access) vs D(S1) computed genuinely
2345 : ! from real wavefunction data (dmats_init's own non-little-group branch), for S1 := g.S0,
2346 : ! g any little-group element of kprime. Physically D(S1) = dmat_star(g) @ D(S0) (composing
2347 : ! the crossing rotation S0 with the k'-internal rotation g: <l,k'|g.S0|n,k_ibz> =
2348 : ! sum_m <l,k'|g|m,k'><m,k'|S0|n,k_ibz>, a resolution of identity over the {psi_m,k'}
2349 : ! basis). Both sides are now independently available (dmat_star via group theory,
2350 : ! D(S0)/D(S1) via genuine WFK data through dmats_init's generalized non-little-group
2351 : ! branch), so this directly tests whether dmat_star's implicit gauge assumption (that
2352 : ! |psi,k'> := S0|psi,k_ibz>) matches the physical wavefunction gstore_compute/kpts_map
2353 : ! would independently reconstruct -- the validation-coverage gap this whole
2354 : ! generalization was built to close (see gstore_symmetrize_status memory). dmat_star's
2355 : ! own placeholder-identity entries (non-little-group g) cannot be told apart from a
2356 : ! genuine identity result by value alone, so little-group membership of kprime is
2357 : ! re-derived directly here, mirroring dmats_get_star_dmats's own test.
2358 0 : nsym = dmats%cryst%nsym
2359 0 : nb = dmats%brange_spin(2, spin) - dmats%brange_spin(1, spin) + 1
2360 0 : kk_ibz = dmats%ks_ebands%kptns(:, ik_ibz)
2361 0 : ABI_MALLOC(d_composed, (nb, nb))
2362 0 : n_cross_ok = 0; n_cross_bad = 0
2363 :
2364 0 : do itime=1,2
2365 0 : tsign = merge(one, -one, itime == 1)
2366 0 : do isym=1,nsym
2367 0 : kk_g = tsign * matmul(transpose(real(dmats%cryst%symrel(:,:,isym), dp)), kprime)
2368 0 : g0 = nint(kprime - kk_g)
2369 0 : if (.not. all(abs(kprime - kk_g - g0) < tol8)) cycle ! Not in the little group of kprime.
2370 :
2371 0 : isym1 = dmats%multable(1, isym, isym0)
2372 0 : itime1 = 1 + mod((itime - 1) + (itime0 - 1), 2)
2373 0 : if (isym1 == 0) then
2374 0 : call wrtout(units, " dmats_check_star cross-check: multable(isym,isym0) not found, skipping")
2375 0 : cycle
2376 : end if
2377 :
2378 0 : d_composed = matmul(dmat_star(:,:,isym,itime), dmats%for_spin(spin)%value(:,:,isym0,itime0,ik_ibz))
2379 0 : maxdiff = maxval(abs(d_composed - dmats%for_spin(spin)%value(:,:,isym1,itime1,ik_ibz)))
2380 0 : if (maxdiff < tol6) then
2381 0 : n_cross_ok = n_cross_ok + 1
2382 : else
2383 0 : n_cross_bad = n_cross_bad + 1
2384 0 : if (prtvol > 0) then
2385 : write(msg,'(a,7(i0,1x),a,es12.4)') &
2386 0 : " dmats_check_star cross-check FAIL: ik_ibz,isym0,itime0,isym,itime,isym1,itime1= ", &
2387 0 : ik_ibz, isym0, itime0, isym, itime, isym1, itime1, " maxdiff= ", maxdiff
2388 0 : call wrtout(units, msg)
2389 : end if
2390 : end if
2391 : end do
2392 : end do
2393 :
2394 : call wrtout(units, sjoin(" dmats_check_star cross-check (dmat_star vs genuine WFK D-matrices): ", &
2395 0 : itoa(n_cross_ok), "/", itoa(n_cross_ok + n_cross_bad), " composed relations agree"))
2396 0 : if (n_cross_bad > 0) ierr = ierr + 1
2397 :
2398 0 : ABI_FREE(d_composed)
2399 0 : ABI_FREE(dmat_star)
2400 :
2401 0 : end subroutine dmats_check_star
2402 : !!***
2403 :
2404 : !!****f* m_classify_bands/dmats_classify
2405 : !! NAME
2406 : !! dmats_classify
2407 : !!
2408 : !! FUNCTION
2409 : !! Classify the KS states based on the computed representation matrices (dmats)
2410 : !! using the irreducible representations of the little group of k.
2411 : !!
2412 : !! SOURCE
2413 :
2414 0 : subroutine dmats_classify(dmats, prtvol)
2415 :
2416 : !Arguments ------------------------------------
2417 : class(dmats_t), target, intent(in) :: dmats
2418 : integer,intent(in) :: prtvol
2419 :
2420 : !Local variables-------------------------------
2421 0 : type(esymm_t) :: Bsym
2422 : integer :: spin, bstart, nb, ik_ibz, idg, iclass, isym_class, sym_idx, isym, tr_isym, ib_start, ib_stop
2423 0 : real(dp), pointer :: ene_k(:)
2424 : real(dp) :: kk_ibz(3)
2425 : ! *************************************************************************
2426 :
2427 0 : do spin=1, size(dmats%for_spin)
2428 0 : bstart = dmats%brange_spin(1, spin)
2429 0 : nb = dmats%brange_spin(2, spin) - bstart + 1
2430 :
2431 0 : do ik_ibz=1, dmats%ks_ebands%nkpt
2432 0 : kk_ibz = dmats%ks_ebands%kptns(:, ik_ibz)
2433 0 : ene_k => dmats%ks_ebands%eig(bstart:dmats%brange_spin(2, spin), ik_ibz, spin)
2434 :
2435 : !only_trace = .false.
2436 : call Bsym%init(kk_ibz, dmats%cryst, .false., dmats%ks_ebands%nspinor, &
2437 0 : bstart, nb, dmats%dtset%symsigma_de, ene_k, tol3)
2438 :
2439 0 : if (Bsym%err_status /= 0) cycle
2440 :
2441 0 : do idg=1, Bsym%ndegs
2442 0 : ib_start = Bsym%degs_bounds(1, idg) ! relative to bstart
2443 0 : ib_stop = Bsym%degs_bounds(2, idg)
2444 :
2445 0 : sym_idx = 0
2446 0 : do iclass=1, Bsym%nclass
2447 0 : do isym_class=1, Bsym%nelements(iclass)
2448 0 : sym_idx = sym_idx + 1
2449 0 : isym = Bsym%sgk2symrec(sym_idx)
2450 0 : associate(cmat => dmats%for_spin(spin)%value(:,:, isym, 1, ik_ibz))
2451 0 : Bsym%Calc_irreps(idg)%mat(:,:,sym_idx) = cmat(ib_start:ib_stop, ib_start:ib_stop)
2452 0 : Bsym%Calc_irreps(idg)%trace(sym_idx) = get_trace(Bsym%Calc_irreps(idg)%mat(:,:,sym_idx))
2453 : end associate
2454 : end do
2455 : end do
2456 :
2457 0 : if (Bsym%can_use_tr) then
2458 0 : do tr_isym=1, Bsym%nsym_trgk
2459 0 : isym = Bsym%tr_sgk2symrec(tr_isym)
2460 0 : associate(cmat => dmats%for_spin(spin)%value(:,:, isym, 2, ik_ibz))
2461 0 : Bsym%trCalc_irreps(idg)%mat(:,:,tr_isym) = cmat(ib_start:ib_stop, ib_start:ib_stop)
2462 0 : Bsym%trCalc_irreps(idg)%trace(tr_isym) = get_trace(Bsym%trCalc_irreps(idg)%mat(:,:,tr_isym))
2463 : end associate
2464 : end do
2465 : end if
2466 : end do
2467 :
2468 0 : call Bsym%finalize(prtvol)
2469 0 : call Bsym%print([std_out, ab_out], prtvol=prtvol)
2470 0 : call esymm_free(Bsym)
2471 : end do
2472 : end do
2473 :
2474 0 : end subroutine dmats_classify
2475 : !!***
2476 :
2477 : !----------------------------------------------------------------------
2478 :
2479 0 : end module m_classify_bands
2480 : !!***
|
