Branch data 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 : : !!***
|
