Line data Source code
1 : !!****m* ABINIT/m_esymm
2 : !! NAME
3 : !! m_esymm
4 : !!
5 : !! FUNCTION
6 : !! Objects and procedures to find the irreducible representations associated to electronic eigenstates.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_esymm
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 :
28 : use m_io_tools, only : file_exists
29 : use m_matrix, only : matr3inv
30 : use m_symtk, only : sg_multable, symrelrot, littlegroup_q
31 : use m_symfind, only : symbrav
32 : use m_fstrings, only : int2char10, itoa, sjoin
33 : use m_numeric_tools, only : print_arr, set2unit, get_trace
34 : use m_hide_lapack, only : xgeev, xginv
35 : use m_crystal, only : crystal_t
36 : use m_defs_ptgroups, only : point_group_t, irrep_t, irrep_free, copy_irrep
37 : use m_ptgroups, only : get_classes, point_group_init, mult_table, sum_irreps
38 :
39 : implicit none
40 :
41 : private
42 : !!***
43 :
44 : ! Error codes.
45 : integer,private,parameter :: ESYM_NOERROR = 0
46 : integer,private,parameter :: ESYM_ACCDEG_ERROR = 10
47 : integer,private,parameter :: ESYM_CLASSIFICATION_ERROR= 11
48 : integer,private,parameter :: ESYM_ORTHO_ERROR = 12
49 : integer,private,parameter :: ESYM_UNITARY_ERROR = 13
50 : integer,private,parameter :: ESYM_PTG_WRONG_MAPPING = 20
51 : integer,private,parameter :: ESYM_HERRING_WRONG_TEST = 30
52 : integer,private,parameter :: ESYM_HEUR_WRONG_NCLASSES = 40
53 : integer,private,parameter :: ESYM_HEUR_WRONG_DIMS = 41
54 :
55 : !----------------------------------------------------------------------
56 :
57 : !!****t* m_esymm/esymm_t
58 : !! NAME
59 : !! esymm_t
60 : !!
61 : !! FUNCTION
62 : !! Dataype gathering data and tables needed to analize the symmetries
63 : !! of electronic states at a given k-point via Group Theory.
64 : !!
65 : !! SOURCE
66 :
67 : type,public :: esymm_t
68 :
69 : integer :: nspinor
70 : ! Number of spinor components.
71 :
72 : integer :: first_ib
73 : ! Index of the first treated band.
74 :
75 : integer :: nbnds
76 : ! Number of bands for this k-point and spin.
77 :
78 : integer :: nclass
79 : ! The number of classes in the group of k.
80 :
81 : integer :: nsym_gk
82 : ! Number of symmetries in the group of k. Namely that the set of symmetries such that Sk = k +G0.
83 :
84 : integer :: nsym_trgk
85 : ! Number of symmetries in the extended group of k. Namely that the set of symmetries such that -Sk = k + G0.
86 :
87 : integer :: err_status = ESYM_NOERROR
88 : ! Flag signaling if the classification algorithm succeed or not.
89 :
90 : real(dp) :: tol_deg
91 : ! Energy tolerance below which two states are considered degenerate.
92 :
93 : logical :: can_use_tr
94 : ! .TRUE. if time-reversal can be used
95 :
96 : logical :: only_trace
97 : ! if .TRUE. only the trace of a single matrix per class is calculated
98 : ! this is the standard way used to analyze bands symmetries. If .FALSE.
99 : ! the full matrices of the irreducible representations are calculated and stored
100 :
101 : logical :: has_spatial_inv
102 : ! .TRUE. if the inversion belongs to the space group
103 :
104 : logical :: nonsymmorphic_at_zoneborder
105 : ! if .TRUE. analysis cannot be performed since kpt is
106 : ! at border zone and non-zero fractional translations are present in the space group
107 :
108 : logical :: has_chtabs
109 : ! True if Ref_irreps and character tables are available (tables are initialized either
110 : ! from point group irreps or from an external database downloaded from the Bilbao server)
111 :
112 : real(dp) :: kpt(3)
113 : ! The crystalline momentum of the wavefunctions in reduced coordinates.
114 :
115 : character(len=500) :: err_msg = "None"
116 : ! Error message:
117 :
118 : integer,allocatable :: g0(:,:)
119 : ! (3,nsym_gk)
120 : ! The umklapp g0 vector associated to each little group operation.
121 :
122 : integer,allocatable :: tr_g0(:,:)
123 : ! (3,nsym_trgk)
124 : ! The umklapp g0 vector associated to each little group operation.
125 :
126 : integer :: ndegs
127 : ! Number of degenerate states.
128 :
129 : integer,allocatable :: nelements(:)
130 : ! (nclass)
131 : ! Number of symmetry operations in each class.
132 :
133 : integer,allocatable :: sgk2symrec(:)
134 : ! (nsym_gk)
135 : ! Mapping between the symmetries of the group of k and the symrec(l) array.
136 : ! The symmetries of the little group are always packed in classes to facilitate
137 : ! the calculation of the character of the irrep. Abinit symmetries are randomly ordered.
138 :
139 : integer,allocatable :: tr_sgk2symrec(:)
140 : ! (nsym_trgk)
141 : ! Mapping between the symmetries of the group of k and the symrec(l) array.
142 : ! The symmetries of the little group are always packed in classes to facilitate
143 : ! the calculation of the character of the irrep. Abinit symmetries are randomly ordered.
144 :
145 : integer,allocatable :: herring_test(:)
146 : ! (nclass)
147 : ! The result of Herring test for each irreducible representantion of the group of k.
148 : ! Possible values are: +1, 0, -1
149 :
150 : integer,allocatable :: b2irrep(:)
151 : ! (nbnds)
152 : ! For each band, the index of the irreducible representation in Ref_irreps.
153 :
154 : type(coeffi1_type),allocatable :: irrep2b(:)
155 : ! irrep2b(0:nclass)%value(:)
156 : ! Ragged arrays with the mapping between the set of irreducible representation and the band indices.
157 : ! irrep2b(irp)%value(:) gives the indices of the states belonging to irrep irp, irp=1,nclass
158 : ! irrep2b(0)%value(:) stores the indices of the states that have not been classified due to
159 : ! the presence of an accidental degeneracy.
160 :
161 : integer,allocatable :: degs_bounds(:,:)
162 : ! degs_bounds(2,ndegs)
163 : ! degs_bounds(1,idg)= first band index of the degenerate set idg=1,ndegs
164 : ! degs_bounds(2,idg)= final band index of the degenerate set idg=1,ndegs
165 :
166 : integer,allocatable :: degs_dim(:)
167 : ! (ndegs)
168 : ! Number of states in each degenerate subspace. Cannot be larger that nclass provided
169 : ! that no accidental degeneracy occurs.
170 :
171 : !% integer,allocatable :: class_ids(:,:)
172 : ! class_ids(2,nclass)
173 : ! (1,icl) = index of the first symmetry of class icl
174 : ! (2,icl) = index of the last symmetry of class icl
175 : ! Note that symmetries in sym are packed in classes.
176 :
177 : type(irrep_t),allocatable :: Calc_irreps(:)
178 : ! (ndegs)
179 : ! The representations of the little group of k calculated from the wavefunctions. <\phi_nk|R_t|\phi_mk>
180 : ! where R_t belong to the little group of k.
181 : ! They represent an unitary irreducible representation provided that no accidental degeneracy occurs.
182 :
183 : type(irrep_t),allocatable :: trCalc_irreps(:)
184 : ! (ndegs)
185 : ! The representations of the little group of k calculated from the wavefunctions. <\phi_nk|R_t|\phi_mk>
186 : ! where R_t belong to the little group of k.
187 : ! They represent an unitary irreducible representation provided that no accidental degeneracy occurs.
188 :
189 : type(irrep_t),allocatable :: Ref_irreps(:)
190 : ! (nclass)
191 : ! Reference irreducible representations of the group of k derived from the point group
192 : ! or from the external database downloaded from the Bilbao web site.
193 :
194 : contains
195 : procedure :: init => esymm_init ! Initialize the object
196 : procedure :: print => esymm_print ! Print info
197 : procedure :: finalize => esymm_finalize ! Finalize the object
198 : procedure :: failed => esymm_failed ! True if symmetry analysis failed.
199 : procedure :: symmetrize_mels => esymm_symmetrize_mels ! Symmetrize given matrix elements
200 : procedure :: free => esymm_free_0D ! Free memory
201 : !procedure :: which_irrep => which_irrep
202 : end type esymm_t
203 :
204 : public :: esymm_free ! Free memory
205 : !!***
206 :
207 : !----------------------------------------------------------------------
208 :
209 : !public :: polish_irreps ! TODO method of Irreps_t, therefore should be moved to m_ptgroups.
210 : ! but first one has to solve the dependency on m_abilasi and scalapack
211 : interface esymm_free
212 : module procedure esymm_free_0D
213 : module procedure esymm_free_2D
214 : end interface esymm_free
215 :
216 : contains
217 : !!***
218 :
219 : !----------------------------------------------------------------------
220 :
221 : !!****f* m_esymm/esymm_init
222 : !! NAME
223 : !! esymm_init
224 : !!
225 : !! FUNCTION
226 : !! Initialize a esymm_t datatype containing data and parameters
227 : !! needed to analyze the irreducible representations at a particular k-point.
228 : !!
229 : !! INPUTS
230 : !! kpt_in(3)=The k-point where the classification of bands is required.
231 : !! Cryst<crystal_t>=Datatype describing the unit cell and its symmetries.
232 : !! nspinor=number of spinorial components
233 : !! nsppol=number of independent polarizations
234 : !! first_ib=Index of the first band.
235 : !! nbnds=Number of bands for this k-point.
236 : !! ene_k(nbnds)=energies for this k-point. ene_k(1) corresponds to band first_ib.
237 : !! EDIFF_TOL=tolerance below which two states are considered to belong to the same irreducible representation
238 : !!
239 : !! NOTES
240 : !! The present implementation does NOT work at zone border if the little group of
241 : !! kpt_in is non-symmorphic namely there is at least a symmetry operation with non-zero tnons.
242 : !!
243 : !! SOURCE
244 :
245 0 : subroutine esymm_init(esymm, kpt_in, Cryst, only_trace, nspinor, first_ib, nbnds, EDIFF_TOL, ene_k, tolsym)
246 :
247 : !Arguments ------------------------------------
248 : !scalars
249 : class(esymm_t),intent(out) :: esymm
250 : integer,intent(in) :: nbnds,nspinor,first_ib
251 : real(dp),intent(in) :: EDIFF_TOL,tolsym
252 : logical,intent(in) :: only_trace
253 : type(crystal_t),intent(in) :: Cryst
254 : !arrays
255 : real(dp),intent(in) :: ene_k(nbnds),kpt_in(3)
256 :
257 : !Local variables-------------------------------
258 : !scalars
259 : integer :: dim_degs,iband,idg,irp,nacc_deg,isym_gk,grp_ierr
260 : integer :: nsym_fm,idx_fm,idx_gk,idx_trgk,isym,jsym,dummy_timrev !,iholohedry
261 : integer :: iel,icls,msym,iord !isym1,!iprod,dim_irrep,icls2, isym2,isym_tr,
262 : integer :: spgroup,chkprim !,ptgroupma
263 : real(dp) :: mkt
264 : !complex(dp) :: phase_k
265 : character(len=5) :: ptgroup,ptgroup_name
266 : character(len=10) :: spgroup_str
267 : character(len=1000) :: msg
268 : character(len=fnlen) :: lgroup_fname
269 : !arrays
270 : integer :: inversion(3,3), bravais(11),sym_axis(3)
271 0 : integer,allocatable :: degs_bounds(:,:),dim_irreps(:)
272 : real(dp) :: pmat1(3,3),pmat2(3,3),pmat3(3,3),pmat4(3,3),pmat5(3,3),pmat6(3,3)
273 : !real(dp) :: genafm(3)
274 : !integer :: rot2(3,3)
275 0 : integer,allocatable :: elements_idx(:,:),tmp_nelements(:), found(:),symrec_fm(:,:,:),fm2symrec(:)
276 0 : integer,allocatable :: ksym_table(:,:,:),sgk(:,:,:),tr_sgk(:,:,:),dum_symafm(:)
277 0 : integer,allocatable :: new_idx(:),new_g0(:,:),tmp_symrec(:,:,:),conv_symrec(:,:,:) !,tr_conv_symrec(:,:,:)
278 0 : integer,allocatable :: dummy_symafm(:) !, mtab(:,:)
279 : real(dp) :: conv_gprimd(3,3),axes(3,3) !,tau2(3)
280 : !complex(dp),allocatable :: her_test(:) !,mat_test(:,:)
281 0 : complex(dp),allocatable :: phase_mkt(:)
282 0 : type(point_group_t) :: Ptg
283 : ! *************************************************************************
284 :
285 0 : esymm%err_status= ESYM_NOERROR
286 : inversion=RESHAPE((/-1,0,0,0,-1,0,0,0,-1/),(/3,3/))
287 :
288 : ! Initialize basic variables
289 0 : esymm%nspinor = nspinor
290 0 : esymm%first_ib = first_ib
291 0 : esymm%nbnds = nbnds
292 0 : esymm%only_trace = only_trace
293 0 : esymm%tol_deg = EDIFF_TOL
294 0 : esymm%has_spatial_inv= (cryst%idx_spatial_inversion() /= 0)
295 0 : esymm%can_use_tr = .TRUE. !TODO this should be input
296 0 : esymm%has_chtabs = .FALSE.
297 0 : esymm%kpt = kpt_in(:)
298 0 : esymm%nonsymmorphic_at_zoneborder=.FALSE.
299 :
300 : ! Locate degenerate_bands ===
301 0 : esymm%ndegs=1
302 0 : ABI_MALLOC(degs_bounds,(2,nbnds))
303 0 : degs_bounds=0; degs_bounds(1,1)=1
304 :
305 0 : do iband=2,nbnds
306 0 : if (ABS(ene_k(iband)-ene_k(iband-1))>EDIFF_TOL) then
307 0 : degs_bounds(2,esymm%ndegs) = iband-1 + (first_ib-1)
308 0 : esymm%ndegs=esymm%ndegs+1
309 0 : degs_bounds(1,esymm%ndegs) = iband + (first_ib-1)
310 : end if
311 : end do
312 0 : degs_bounds(2,esymm%ndegs)=nbnds + (first_ib-1)
313 :
314 0 : ABI_MALLOC(esymm%degs_bounds,(2,esymm%ndegs))
315 0 : esymm%degs_bounds = degs_bounds(:,1:esymm%ndegs)
316 0 : ABI_FREE(degs_bounds)
317 :
318 : ! Each band is initialized as "Unknown".
319 0 : ABI_MALLOC(esymm%b2irrep,(esymm%nbnds))
320 0 : esymm%b2irrep = 0
321 :
322 : ! Find the group of kpt_in.
323 : ! The small point group is the subset of symrec such that $ S q = q + g0 $
324 : ! Symmetries are packed in classes.
325 : ! For the time being, AFM symmetries are not treated.
326 :
327 0 : write(msg,'(a,3(1x,f7.4))')" Finding the little group of k-point: ",esymm%kpt
328 0 : call wrtout(std_out,msg)
329 :
330 : ! Only FM symmetries are used.
331 0 : nsym_fm = COUNT(Cryst%symafm==1)
332 :
333 0 : if (nsym_fm /= Cryst%nsym) then
334 0 : write(msg,'(4a)')ch10,&
335 0 : "Band classification in terms of magnetic space groups not coded! ",ch10,&
336 0 : "Only the ferromagnetic subgroup will be used "
337 0 : ABI_COMMENT(msg)
338 : end if
339 :
340 0 : ABI_MALLOC(symrec_fm, (3,3,nsym_fm))
341 0 : ABI_MALLOC(fm2symrec, (nsym_fm))
342 :
343 0 : idx_fm = 0
344 0 : do isym=1,Cryst%nsym
345 0 : if (Cryst%symafm(isym) == 1) then
346 0 : idx_fm = idx_fm + 1
347 0 : symrec_fm(:,:,idx_fm) = Cryst%symrec(:,:,isym)
348 0 : fm2symrec(idx_fm) = isym
349 : end if
350 : end do
351 :
352 : ! Find symmetries that preserve k.
353 0 : ABI_MALLOC(ksym_table,(4,2,nsym_fm))
354 0 : ABI_MALLOC(dummy_symafm,(nsym_fm))
355 :
356 0 : dummy_symafm = 1
357 0 : call littlegroup_q(nsym_fm, esymm%kpt, ksym_table, symrec_fm, dummy_symafm, dummy_timrev, prtvol=0)
358 :
359 0 : esymm%nsym_gk =COUNT(ksym_table(4,1,:)==1) ! # S such that S k = k +G0
360 :
361 0 : esymm%nsym_trgk=0
362 0 : if (esymm%can_use_tr) esymm%nsym_trgk=COUNT(ksym_table(4,2,:)==1) ! # S such that -S k = k +G0
363 :
364 : ! Allocate workspace arrays.
365 0 : ABI_MALLOC(sgk,(3,3,esymm%nsym_gk))
366 0 : ABI_MALLOC(tr_sgk,(3,3,esymm%nsym_trgk))
367 :
368 : ! Allocate mapping little-group --> symrec and table for umklapps.
369 0 : ABI_MALLOC(esymm%sgk2symrec,(esymm%nsym_gk))
370 0 : ABI_MALLOC(esymm%g0,(3,esymm%nsym_gk))
371 0 : ABI_MALLOC(esymm%tr_sgk2symrec,(esymm%nsym_trgk))
372 0 : ABI_MALLOC(esymm%tr_g0,(3,esymm%nsym_trgk))
373 :
374 : ! Important NOTE:
375 : ! If nonsymmorphic_at_zoneborder symmetry analysis cannot be performed unless
376 : ! an external database retrieved from the bilbao server (REPRES) is found.
377 0 : idx_gk=0; idx_trgk=0
378 0 : esymm%sgk2symrec=-999; esymm%tr_sgk2symrec=-999
379 :
380 0 : do isym=1,nsym_fm
381 0 : if (ksym_table(4,1,isym)==1) then ! S k = k +G0
382 0 : idx_gk=idx_gk+1
383 0 : sgk(:,:,idx_gk)=symrec_fm(:,:,isym)
384 0 : esymm%g0(:,idx_gk)=ksym_table(1:3,1,isym)
385 0 : esymm%sgk2symrec(idx_gk)=fm2symrec(isym)
386 0 : if (ANY(ksym_table(1:3,1,isym)/=0).and.(ANY(ABS(Cryst%tnons(:,fm2symrec(isym)))>tol6))) then
387 0 : esymm%nonsymmorphic_at_zoneborder=.TRUE.
388 : end if
389 : end if
390 :
391 0 : if (esymm%can_use_tr.and.ksym_table(4,2,isym)==1) then ! -S k = k +G0
392 0 : idx_trgk=idx_trgk+1
393 0 : tr_sgk(:,:,idx_trgk)=symrec_fm(:,:,isym)
394 0 : esymm%tr_g0(:,idx_trgk)=ksym_table(1:3,2,isym)
395 0 : esymm%tr_sgk2symrec(idx_trgk)=fm2symrec(isym)
396 : end if
397 : end do
398 :
399 0 : ABI_FREE(ksym_table)
400 0 : ABI_FREE(symrec_fm)
401 0 : ABI_FREE(fm2symrec)
402 :
403 : ! Divide the operations into classes.
404 0 : ABI_MALLOC(dum_symafm, (esymm%nsym_gk))
405 0 : dum_symafm = 1
406 :
407 : ! Check group closure
408 0 : call sg_multable(esymm%nsym_gk, dum_symafm, sgk, grp_ierr)
409 0 : ABI_CHECK_IEQ(grp_ierr, 0, "sg_multable failed")
410 0 : ABI_FREE(dum_symafm)
411 :
412 0 : ABI_MALLOC(tmp_nelements, (esymm%nsym_gk))
413 0 : ABI_MALLOC(elements_idx, (esymm%nsym_gk, esymm%nsym_gk))
414 :
415 0 : call get_classes(esymm%nsym_gk, sgk, esymm%nclass, tmp_nelements, elements_idx)
416 :
417 0 : ABI_MALLOC(esymm%nelements, (esymm%nclass))
418 0 : esymm%nelements = tmp_nelements(1:esymm%nclass)
419 0 : ABI_FREE(tmp_nelements)
420 :
421 : ! From the list of symmetry operations and the lattice vectors, determine the
422 : ! Bravais information including the holohedry, the centering, the coordinate of
423 : ! the primitive vectors in the conventional vectors, as well as the point group,
424 0 : msym=192; if (allocated(Cryst%symrec)) msym=size(Cryst%symrec,3)
425 0 : ABI_MALLOC(tmp_symrec,(3,3,msym))
426 0 : tmp_symrec(:,:,1:esymm%nsym_gk)=sgk
427 :
428 0 : call symbrav(bravais,msym,esymm%nsym_gk,ptgroup,Cryst%gprimd,tmp_symrec,tolsym,axis=sym_axis)
429 :
430 0 : ABI_FREE(tmp_symrec)
431 :
432 0 : write(std_out,'(a)')" symptgroup returned point group: "//TRIM(ptgroup)
433 0 : write(std_out,'(a,i2)')" iholohedry ",bravais(1)
434 0 : write(std_out,'(a,i2)')" center ",bravais(2)
435 0 : write(std_out,'(a,9i3)')" gprimd in the axes of the conventional bravais lattice (*2 if center/=0)",bravais(3:11)
436 0 : write(std_out,'(a,3i3)')" sym_axis ",sym_axis
437 :
438 : ! Branching:
439 : ! 1) If the little group is not symmorphic_at_zoneborder we can
440 : ! classify the states using the irreducible representation of the point group.
441 : !
442 : ! 2) If the little group is symmorphic_at_zoneborder, we have to rely on
443 : ! an external database retrieved from the Bilbao server in order to classify the states.
444 : ! If the file is not available, we only know the number of classes but neither their
445 : ! character nor the dimension of the irreducible representation.
446 : !
447 0 : if (esymm%nonsymmorphic_at_zoneborder) then
448 :
449 0 : spgroup=0
450 0 : chkprim=1 ! Cell must be primitive.
451 : !call symlatt(bravais,std_out,msym,nptsym,ptsymrel,rprimd,tolsym)
452 : !call symspgr(bravais,Cryst%nsym,spgroup,Cryst%symrel,Cryst%tnons,tolsym)
453 :
454 : !call symanal(bravais,chkprim,genafm,msym,nsym,ptgroupma,rprimd,spgroup,symafm,symrel,tnons,tolsym)
455 :
456 0 : call int2char10(spgroup, spgroup_str)
457 0 : lgroup_fname = "lgroup_"//TRIM(spgroup_str)
458 :
459 0 : if (file_exists(lgroup_fname)) then
460 0 : ABI_ERROR("Not coded")
461 : ! Read little groups from the external database.
462 : !% call init_groupk_from_file(Lgrp,spgroup,lgroup_fname,ierr)
463 :
464 : ! Save the irreducible representations in esymm.
465 : ! Reorder symmetries such that they correspond to the Bilbao database.
466 : !% allocate(esymm%Ref_irreps(esymm%nclass))
467 : !% call copy_irrep(Irreps, esymm%Ref_irreps)
468 :
469 : else
470 : write(msg,'(7a)')&
471 0 : "Non-symmorphic small group and zone border. ",ch10,&
472 0 : "External file: ",TRIM(lgroup_fname)," containing Bilbao tables not found ",ch10,&
473 0 : "Character analysis cannot be performed. Accidental degeneracies cannot be detected. "
474 0 : ABI_WARNING(msg)
475 :
476 0 : esymm%has_chtabs = .FALSE.
477 :
478 : ! Reorder indices such that symmetries are packed in classes.
479 0 : ABI_MALLOC(new_idx,(esymm%nsym_gk))
480 0 : ABI_MALLOC(new_g0,(3,esymm%nsym_gk))
481 0 : new_g0=0; iord = 0
482 0 : do icls=1,esymm%nclass
483 0 : do iel=1,esymm%nelements(icls)
484 0 : iord = iord+1
485 0 : jsym = elements_idx(iel,icls)
486 0 : new_idx(iord) = esymm%sgk2symrec(jsym)
487 0 : new_g0(:,iord) = esymm%g0(:,jsym)
488 : end do
489 : end do
490 :
491 0 : esymm%sgk2symrec = new_idx
492 0 : esymm%g0 = new_g0
493 :
494 0 : ABI_FREE(new_idx)
495 0 : ABI_FREE(new_g0)
496 : end if ! file exists
497 :
498 : else
499 : !
500 : ! **** This part is still under development. It might not work for particular ****
501 : ! **** orientations of the unit cell or particular lattices. ****
502 : !
503 : ! The symmetries in the Bilbao database refer to the conventional unit cells.
504 : ! Therefore we have to map the abinit symmetries (in reduced coordinates)
505 : ! onto the Bilbao dataset. Bilbao standard settings are:
506 : !
507 : ! * unique axis b (cell choice 1) for space groups withing the monoclinic system
508 : ! * obverse triple hexagonal unit cell R space groups.
509 : ! * origin choice two - inversion center at (0, 0, 0) - for the centrosymmetric
510 : ! space groups for which there are two origins choices, within the
511 : ! orthorombic, tetragonal and cubic system.
512 :
513 : ! 1) Retrieve the rotation matrices and the irreducible representations (Bilbao setting).
514 0 : call point_group_init(Ptg,ptgroup)
515 :
516 0 : esymm%has_chtabs = .TRUE.
517 0 : ABI_CHECK(esymm%nclass == Ptg%nclass,"esymm%nclass/=Ptg%nclass!")
518 :
519 0 : do icls=1,esymm%nclass ! FIXME this is awful, should be done in a cleaner way.
520 0 : esymm%nelements(icls)=Ptg%class_ids(2,icls) - Ptg%class_ids(1,icls) + 1
521 : end do
522 :
523 : ! 2) Generate the symmetry operations in the conventional vector coordinates.
524 0 : conv_gprimd(:,1)=bravais(3:5)
525 0 : conv_gprimd(:,2)=bravais(6:8)
526 0 : conv_gprimd(:,3)=bravais(9:11)
527 :
528 0 : axes = conv_gprimd
529 0 : call matr3inv(conv_gprimd,axes) !; axes=TRANSPOSE(axes)
530 :
531 0 : conv_gprimd=MATMUL(Cryst%gprimd,TRANSPOSE(axes))
532 : !conv_gprimd=MATMUL(axes,Cryst%gprimd)
533 : !conv_gprimd=MATMUL(TRANSPOSE(axes),Cryst%gprimd)
534 : !write(std_out,*)"conv_gprimd:", conv_gprimd
535 :
536 0 : ptgroup_name = ADJUSTL(ptgroup)
537 :
538 0 : select case (ptgroup_name)
539 :
540 : case ("3m","-3m")
541 0 : call wrtout(std_out," Changing the conventional cell: rhombohedral --> triple hexagonal")
542 : ! Transformation matrices: primitive rhombohedral --> triple hexagonal cell obverse setting. Table 5.1.3.1 ITA page 81.
543 0 : pmat1 = RESHAPE( (/ 1,-1, 0, 0, 1,-1, 1, 1, 1/), (/3,3/) ) ! R1
544 : pmat2 = RESHAPE( (/ 0, 1,-1,-1, 0, 1, 1, 1, 1/), (/3,3/) ) ! R2
545 : pmat3 = RESHAPE( (/-1, 0, 1, 1,-1, 0, 1, 1, 1/), (/3,3/) ) ! R3
546 : pmat4 = RESHAPE( (/-1, 1, 0, 0,-1, 1, 1, 1, 1/), (/3,3/) ) ! R1 reverse setting.
547 : pmat5 = RESHAPE( (/ 0,-1, 1, 1, 0,-1, 1, 1, 1/), (/3,3/) ) ! R2 reverse setting.
548 : pmat6 = RESHAPE( (/ 1, 0,-1,-1, 1, 0, 1, 1, 1/), (/3,3/) ) ! R3 reverse setting.
549 0 : conv_gprimd = MATMUL(conv_gprimd,pmat1)
550 : !conv_gprimd = MATMUL(conv_gprimd,pmat2)
551 : !conv_gprimd = MATMUL(conv_gprimd,pmat3)
552 : !conv_gprimd = MATMUL(conv_gprimd,pmat4)
553 : !conv_gprimd = MATMUL(conv_gprimd,pmat5)
554 : !conv_gprimd = MATMUL(conv_gprimd,pmat6)
555 : !write(std_out,*)" New conv_gprimd:", conv_gprimd
556 :
557 : case ("mm2")
558 0 : call wrtout(std_out," Changing the conventional cell: unconventional orthorhombic setting --> conventional")
559 : ! Transformation matrices: unconvential orthorhombic --> conventional orthorhombic. Table 5.1.3.1 ITA page 81.
560 : pmat1 = RESHAPE( (/ 0, 1, 0, 1, 0, 0, 0, 0,-1/), (/3,3/) ) ! ( b, a,-c) --> (a,b,c)
561 0 : pmat2 = RESHAPE( (/ 0, 1, 0, 0, 0, 1, 1, 0, 0/), (/3,3/) ) ! ( c, a, b) --> (a,b,c)
562 : pmat3 = RESHAPE( (/ 0, 0, 1, 0, 1, 0,-1, 0, 0/), (/3,3/) ) ! (-c, b, a) --> (a,b,c)
563 : pmat4 = RESHAPE( (/ 0, 0, 1, 1, 0, 0, 0, 1, 0/), (/3,3/) ) ! ( b, c, a) --> (a,b,c)
564 : pmat5 = RESHAPE( (/ 1, 0, 0, 0, 0, 1, 0,-1, 0/), (/3,3/) ) ! ( a,-c, b) --> (a,b,c)
565 0 : conv_gprimd = MATMUL(conv_gprimd,pmat2)
566 : !write(std_out,*)" New conv_gprimd:", conv_gprimd
567 : case default
568 0 : continue
569 : end select
570 :
571 0 : ABI_MALLOC(conv_symrec,(3,3,esymm%nsym_gk))
572 0 : conv_symrec = sgk
573 :
574 : !axes=zero; axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one
575 : !call symrelrot(esymm%nsym_gk,conv_gprimd,axes,conv_symrec,tolsym)
576 0 : call symrelrot(esymm%nsym_gk,Cryst%gprimd,conv_gprimd,conv_symrec,tolsym)
577 :
578 : ! 3) Reorder indices such that symmetries are packed in classes.
579 0 : ABI_MALLOC(found,(esymm%nsym_gk))
580 0 : ABI_MALLOC(new_idx,(esymm%nsym_gk))
581 0 : ABI_MALLOC(new_g0,(3,esymm%nsym_gk))
582 0 : new_g0=0; found=0
583 :
584 0 : do isym=1,esymm%nsym_gk
585 0 : do jsym=1,esymm%nsym_gk
586 0 : if (ALL(Ptg%sym(:,:,isym) == conv_symrec(:,:,jsym) )) then
587 0 : found(isym) = found(isym) + 1
588 0 : new_idx(isym) = esymm%sgk2symrec(jsym)
589 0 : new_g0(:,isym) = esymm%g0(:,jsym)
590 : !EXIT
591 : end if
592 : end do
593 : end do
594 : !
595 : ! DEBUGGING SECTION
596 : !do isym=1,esymm%nsym_gk
597 : ! jsym=esymm%sgk2symrec(isym)
598 : ! call print_symmetries(1,Cryst%symrec(:,:,jsym),Cryst%tnons(:,jsym),Cryst%symafm(jsym))
599 : ! write(std_out,*)esymm%g0(:,isym)
600 : !end do
601 :
602 0 : if ( Ptg%nsym/=esymm%nsym_gk .or. ANY(found/=1) ) then
603 : !write(std_out,*)Ptg%nsym, esymm%nsym_gk
604 : !write(std_out,'(a,(i2))')" found = ",found
605 0 : write(std_out,*)" Ptg%sym list, conv_symrec list, found Ptg% "
606 0 : do isym=1,Ptg%nsym
607 0 : write(std_out,'(a,i2,a,9i2,4x,a,9i2)')" found ",found(isym)," Ptg ",Ptg%sym(:,:,isym),"conv_symrec ",conv_symrec(:,:,isym)
608 : end do
609 0 : msg = " sgk and esymm%Ptg are inconsistent. Check tables or source"
610 0 : ABI_WARNING(msg)
611 0 : esymm%err_msg = msg(1:500)
612 0 : esymm%err_status = ESYM_PTG_WRONG_MAPPING
613 0 : esymm%has_chtabs = .FALSE.
614 :
615 : else ! Reorder symmetries.
616 0 : esymm%sgk2symrec = new_idx
617 0 : esymm%g0 = new_g0
618 : end if
619 :
620 0 : ABI_FREE(new_idx)
621 0 : ABI_FREE(new_g0)
622 0 : ABI_FREE(found)
623 0 : ABI_FREE(conv_symrec)
624 :
625 0 : if (esymm%has_chtabs) then
626 : ! Multiply the point group irreps by e^{-ik.\tau} to have the irreps of the little group.
627 : ! Store the results in esymm%Ref_irreps so that one can classify the states afterwards.
628 0 : ABI_MALLOC(esymm%Ref_irreps,(esymm%nclass))
629 0 : ABI_MALLOC(phase_mkt,(esymm%nsym_gk))
630 :
631 0 : do isym_gk=1,esymm%nsym_gk
632 0 : isym = esymm%sgk2symrec(isym_gk)
633 0 : mkt = -two_pi * DOT_PRODUCT(esymm%kpt, Cryst%tnons(:,isym))
634 0 : phase_mkt(isym_gk) = CMPLX(DCOS(mkt), DSIN(mkt))
635 : end do
636 :
637 0 : call copy_irrep(Ptg%Irreps, esymm%Ref_irreps, phase_mkt)
638 0 : ABI_FREE(phase_mkt)
639 : end if
640 :
641 : #if 0
642 : ! Herring test requires the evaluation of the expression:
643 : !
644 : ! sum_{S,\tau} \chi^{k,\alpha} ({S|\tau}^2)
645 : !
646 : ! where Sk = -k + g0, and \chi is the trace of the \alpha-th
647 : ! irreducible representation of the little group of k.
648 : ! \chi^{k,\alpha} = e^{-ik.\tau} \chi(\alpha) provided that
649 : ! we are not at zone border with a non-symmorphic operation.
650 : ! The expression is always real and it can only be equal to \pm Ptg%nsym or zero.
651 : ! FIXME this part has to be rewritten from scratch.
652 : !if (esymm%err_status/=esymm_NOERROR) then
653 : ! write(std_out,*)" Skipping Herring test"
654 : ! goto 110
655 : !end if
656 :
657 : if (esymm%can_use_tr) then
658 : ABI_MALLOC(her_test,(esymm%nclass))
659 :
660 : ABI_MALLOC(tr_conv_symrec,(3,3,esymm%nsym_trgk))
661 : do isym_tr=1,esymm%nsym_trgk
662 : isym = esymm%tr_sgk2symrec(isym_tr)
663 : tr_conv_symrec(:,:,isym_tr)=Cryst%symrec(:,:,isym)
664 : end do
665 :
666 : call symrelrot(esymm%nsym_trgk,Cryst%gprimd,conv_gprimd,tr_conv_symrec_tr,tolsym)
667 :
668 : do isym_tr=1,esymm%nsym_trgk
669 : isym = esymm%tr_sgk2symrec(isym_tr)
670 : !rot2 = MATMUL(tr_sgk(:,:,isym),tr_sgk(:,:,isym))
671 : !tau2 = MATMUL(tr_sgk(:,:,isym),Cryst%tnons(:,isym)) + Cryst%tnons(:,isym)
672 :
673 : rot2 = MATMUL(tr_conv_symrec(:,:,isym_tr),tr_conv_symrec(:,:,isym_tr))
674 : tau2 = MATMUL(tr_conv_symrec(:,:,isym_tr),Cryst%tnons(:,isym)) + Cryst%tnons(:,isym)
675 :
676 : phase_k = EXP(-j_dpc*two_pi*DOT_PRODUCT(kpoint,tau2))
677 : call locate_sym(Ptg,rot2,isym2,icls2)
678 :
679 : do irp=1,esymm%nclass
680 : her_test(irp) = her_test(irp) + phase_k * Ptg%Irreps(irp)%trace(icls2)
681 : end do
682 : end do
683 :
684 : ABI_FREE(tr_conv_symrec)
685 :
686 : ! FIXME
687 : ABI_MALLOC(esymm%herring_test,(esymm%nclass))
688 :
689 : do irp=1,esymm%nclass
690 : if ( ABS(her_test(irp) - Ptg%nsym) < tol6 ) then
691 : esymm%herring_test(irp) = +1
692 : else if ( ABS(her_test(irp)) < tol6 ) then
693 : esymm%herring_test(irp) = 0
694 : else if ( ABS(her_test(irp) + Ptg%nsym) < tol6 ) then
695 : esymm%herring_test(irp) = -1
696 : else
697 : write(msg,'(a,i0,2a,i0,a,i0)')&
698 : "Herring test for the irreducible representation number ",irp,ch10,&
699 : "gave ",esymm%herring_test(irp),", while it should be 0 or +- ",Ptg%nsym
700 : ABI_WARNING(msg)
701 : esymm%err_msg = msg
702 : esymm%err_status = ESYMM_HERRING_WRONG_TEST
703 : end if
704 : end do
705 :
706 : ABI_FREE(her_test)
707 : end if ! can_use_tr
708 : #endif
709 : !
710 : ! Final check
711 : !allocate(mtab(esymm%nsym_gk,esymm%nsym_gk))
712 : !call mult_table(esymm%nsym_gk,Ptg%sym,mtab)
713 :
714 : !do isym=1,esymm%nsym_gk
715 : ! isym1 = esymm%sgk2symrec(isym)
716 : ! do jsym=1,esymm%nsym_gk
717 : ! isym2 = esymm%sgk2symrec(jsym)
718 : ! rot2 = MATMUL(Cryst%symrec(:,:,isym1),Cryst%symrec(:,:,isym2))
719 :
720 : ! iprod = mtab(isym,jsym)
721 :
722 : ! do irp=1,esymm%nclass
723 : ! dim_irrep = Ptg%Irreps(irp)%dim
724 : ! allocate(mat_test(dim_irrep,dim_irrep))
725 : ! mat_test = Ptg%Irreps(irp)%mat(:,:,isym) * Ptg%Irreps(irp)%mat(:,:,jsym)
726 : ! !call locate_sym(Ptg,rot2,isym2,icls2)
727 : ! write(std_out,*)mat_test - Ptg%Irreps(irp)%mat(:,:,iprod)
728 : ! deallocate(mat_test)
729 : ! end do
730 : !
731 : ! end do
732 : !end do
733 : !
734 : !deallocate(mtab)
735 : end if
736 :
737 0 : ABI_FREE(sgk)
738 0 : ABI_FREE(tr_sgk)
739 0 : ABI_FREE(elements_idx)
740 :
741 : !% allocate(esymm%irrep2b(0:esymm%nclass))
742 : !% call nullify_coeff(esymm%irrep2b)
743 : !
744 : ! 1) Allocate space for the irreducible representations.
745 :
746 : ! 2) Try to determine if we are in presence of an accidental degeneracy. Sufficient condition:
747 : ! There exists a set of degenerate states whose dimension is greater than the dimension
748 : ! of the irreducible representations of the point group. The check can be done only
749 : ! if Character tables are available.
750 :
751 0 : if (esymm%has_chtabs) then
752 0 : ABI_MALLOC(dim_irreps,(esymm%nclass))
753 0 : dim_irreps = (/(esymm%Ref_irreps(irp)%dim, irp=1,esymm%nclass)/)
754 : end if
755 :
756 0 : nacc_deg=0
757 0 : ABI_MALLOC(esymm%degs_dim,(esymm%ndegs))
758 0 : ABI_MALLOC(esymm%Calc_irreps,(esymm%ndegs))
759 :
760 0 : if (esymm%can_use_tr) then
761 0 : ABI_MALLOC(esymm%trCalc_irreps,(esymm%ndegs))
762 : end if
763 :
764 0 : do idg=1,esymm%ndegs
765 0 : dim_degs=esymm%degs_bounds(2,idg)-esymm%degs_bounds(1,idg)+1
766 :
767 0 : if (esymm%has_chtabs) then
768 0 : if (ALL(dim_degs /= dim_irreps)) then ! An accidental degeneracy is present.
769 0 : nacc_deg=nacc_deg+1
770 : end if
771 : end if
772 :
773 0 : esymm%degs_dim(idg) = dim_degs
774 :
775 0 : call esymm%Calc_irreps(idg)%init(esymm%nsym_gk, dim_degs)
776 0 : if (esymm%can_use_tr) call esymm%trCalc_irreps(idg)%init(esymm%nsym_trgk, dim_degs)
777 : end do ! idg
778 :
779 0 : if (esymm%has_chtabs) then
780 0 : ABI_FREE(dim_irreps)
781 0 : if (nacc_deg /= 0) then
782 0 : write(msg,'(a,i0,a)')" Detected ",nacc_deg," accidental degeneracies."
783 0 : ABI_WARNING(msg)
784 0 : esymm%err_status = ESYM_ACCDEG_ERROR
785 : ! TODO this should signal to the caller that we have to decompose the calculated representation.
786 0 : esymm%err_msg = msg(1:500)
787 : end if
788 : end if
789 :
790 :
791 0 : ABI_FREE(dummy_symafm)
792 0 : call Ptg%free()
793 :
794 : DBG_EXIT("COLL")
795 :
796 0 : end subroutine esymm_init
797 : !!***
798 :
799 : !----------------------------------------------------------------------
800 :
801 : !!****f* m_esymm/esymm_print
802 : !! NAME
803 : !! esymm_print
804 : !!
805 : !! FUNCTION
806 : !!
807 : !! INPUTS
808 : !!
809 : !! OUTPUT
810 : !! only printing
811 : !!
812 : !! SOURCE
813 :
814 0 : subroutine esymm_print(esymm, units, prtvol)
815 :
816 : !Arguments ------------------------------------
817 : class(esymm_t),intent(in) :: esymm
818 : integer,intent(in) :: units(:), prtvol
819 :
820 : !Local variables-------------------------------
821 : integer :: icl, idg, irr_idx, nstates, nunknown, istart, istop, ii
822 : character(len=1000) :: fmt, msg, msg0
823 : ! *********************************************************************
824 :
825 0 : write(fmt, *)'(2a,3f8.4,3a,i4,2a,i3,2a,i2,2a,i2,a,',esymm%nclass,'i2,a)'
826 0 : write(msg, fmt) ch10,&
827 0 : ' ===== Character of bands at k-point: ',esymm%kpt,' ===== ',ch10,&
828 0 : ' Total number of bands analyzed .................. ',esymm%nbnds,ch10,&
829 0 : ' Number of degenerate sets detected .............. ',esymm%ndegs,ch10,&
830 0 : ' Number of operations in the little group of k ... ',esymm%nsym_gk,ch10,&
831 0 : ' Number of classes (irreps) in the group of k .... ',esymm%nclass,' (',(esymm%nelements(icl),icl=1,esymm%nclass),' )'
832 0 : call wrtout(units, msg)
833 :
834 0 : if (esymm%nonsymmorphic_at_zoneborder) then
835 0 : call wrtout(units," Non-symmorphic small group at zone border. Character analysis not available ")
836 : end if
837 :
838 0 : if (esymm_failed(esymm)) then
839 0 : write(std_out,'(3a)')"Band classification algorithm failed with the error:",ch10,TRIM(esymm%err_msg)
840 0 : write(msg,'(3a)')"Band classification algorithm failed with the error:",ch10,TRIM(esymm%err_msg)
841 0 : call wrtout(units, msg)
842 : end if
843 :
844 : !nunknown=0
845 : !do iband=1,esymm%nbnds
846 : ! irr_idx = esymm%b2irrep(iband)
847 : ! if (irr_idx /= 0) then
848 : ! if ( esymm%has_chtabs) irr_name = esymm%Ref_Irreps(irr_idx)%name
849 : ! if (.not.esymm%has_chtabs) write(irr_name,'(i0)')irr_idx ! use the index instead of the name.
850 : ! else
851 : ! irr_name = "???"
852 : ! nunknown = nunknown +1
853 : ! end if
854 : ! write(msg,'(a,i3,2a)')' Band ',iband,' belongs to irrep ',TRIM(irr_name)
855 : ! call wrtout(units, msg)
856 : !end do
857 :
858 0 : do irr_idx=1,esymm%nclass
859 0 : nstates = size(esymm%irrep2b(irr_idx)%value)
860 0 : if (esymm%has_chtabs) then
861 0 : write(msg0,'(a,i0,3a)')" Found ",nstates," states with character ",TRIM(esymm%Ref_irreps(irr_idx)%name),": "
862 : else
863 0 : write(msg0,'(2(a,i0),a)')" Found ",nstates," states with character index ",irr_idx,": "
864 : end if
865 0 : do istart=1,nstates,20
866 0 : istop=istart+11; if (istop>nstates) istop=nstates
867 0 : write(msg,'(20(1x,i0))')(esymm%irrep2b(irr_idx)%value(ii), ii=istart,istop)
868 0 : if (istart==1) msg = TRIM(msg0)//TRIM(msg)
869 0 : if (istart/=1) msg = " "//TRIM(msg)
870 0 : call wrtout(units, msg)
871 : end do
872 : end do
873 :
874 0 : nunknown = size(esymm%irrep2b(0)%value)
875 0 : if (nunknown > 0) then
876 0 : write(msg0,'(a,i0,a)')" WARNING: ",nunknown," states have not been classified:"
877 0 : do istart=1,nunknown,20
878 0 : istop=istart+11; if (istop>nunknown) istop=nunknown
879 0 : write(msg,'(20(1x,i0))')(esymm%irrep2b(0)%value(ii), ii=istart,istop)
880 0 : if (istart==1) msg = TRIM(msg0)//TRIM(msg)
881 0 : if (istart/=1) msg = " "//TRIM(msg)
882 0 : call wrtout(units, msg)
883 : end do
884 : end if
885 :
886 0 : if (prtvol > 0 .or. nunknown > 0 .or. .not.esymm%has_chtabs) then
887 : ! print the calculated character table.
888 0 : call wrtout(units,ch10//" Calculated character table ")
889 : !write(fmt,*)'(i2,a,i2,1x,',esymm%nclass,'(a,2f6.3),a)'
890 0 : write(fmt,*)'(i2,a,i2,1x,',esymm%nclass,'(a,2f5.2),a)'
891 0 : do idg=1,esymm%ndegs
892 : write(msg, fmt) &
893 0 : esymm%degs_bounds(1,idg),'-',esymm%degs_bounds(2,idg),&
894 0 : ('|',esymm%Calc_irreps(idg)%trace(esymm%nelements(icl)), icl=1,esymm%nclass),'|'
895 0 : call wrtout(units, msg)
896 : end do
897 : end if
898 :
899 0 : end subroutine esymm_print
900 : !!***
901 :
902 : !----------------------------------------------------------------------
903 :
904 : !!****f* m_esymm/esymm_free_0D
905 : !! NAME
906 : !! esymm_free_0D
907 : !!
908 : !! FUNCTION
909 : !! Deallocate the memory allocated in the esymm_t datatype (scalar version)
910 : !!
911 : !! SOURCE
912 :
913 1232 : subroutine esymm_free_0D(esymm)
914 :
915 : !Arguments ------------------------------------
916 : class(esymm_t),intent(inout) :: esymm
917 :
918 : !Local variables ------------------------------
919 : integer :: ii
920 : ! *************************************************************************
921 :
922 1232 : ABI_SFREE(esymm%g0)
923 1232 : ABI_SFREE(esymm%tr_g0)
924 1232 : ABI_SFREE(esymm%nelements)
925 1232 : ABI_SFREE(esymm%sgk2symrec)
926 1232 : ABI_SFREE(esymm%tr_sgk2symrec)
927 1232 : ABI_SFREE(esymm%herring_test)
928 1232 : ABI_SFREE(esymm%b2irrep)
929 1232 : ABI_SFREE(esymm%degs_bounds)
930 1232 : ABI_SFREE(esymm%degs_dim)
931 :
932 1232 : if (allocated(esymm%irrep2b)) then
933 0 : do ii=LBOUND(esymm%irrep2b,DIM=1),UBOUND(esymm%irrep2b,DIM=1)
934 0 : ABI_FREE(esymm%irrep2b(ii)%value)
935 : end do
936 0 : ABI_FREE(esymm%irrep2b)
937 : end if
938 :
939 1232 : if (allocated(esymm%Calc_irreps)) then
940 0 : call irrep_free(esymm%Calc_irreps)
941 0 : ABI_FREE(esymm%Calc_irreps)
942 : end if
943 1232 : if (allocated(esymm%trCalc_irreps)) then
944 0 : call irrep_free(esymm%trCalc_irreps)
945 0 : ABI_FREE(esymm%trCalc_irreps)
946 : end if
947 1232 : if (allocated(esymm%Ref_irreps)) then
948 0 : call irrep_free(esymm%Ref_irreps)
949 0 : ABI_FREE(esymm%Ref_irreps)
950 : end if
951 :
952 1232 : end subroutine esymm_free_0D
953 : !!***
954 :
955 : !----------------------------------------------------------------------
956 :
957 : !!****f* m_esymm/esymm_free_2D
958 : !! NAME
959 : !! esymm_free_2D
960 : !!
961 : !! FUNCTION
962 : !! Deallocate the memory allocated in the esymm_t datatype (2D version)
963 : !!
964 : !! SOURCE
965 :
966 201 : subroutine esymm_free_2D(esymm)
967 :
968 : !Arguments ------------------------------------
969 : class(esymm_t),intent(inout) :: esymm(:,:)
970 :
971 : !Local variables ------------------------------
972 : integer :: id1,id2
973 : ! *************************************************************************
974 :
975 406 : do id2=1,SIZE(esymm,DIM=2)
976 1638 : do id1=1,SIZE(esymm,DIM=1)
977 1437 : call esymm_free_0D(esymm(id1,id2))
978 : end do
979 : end do
980 :
981 201 : end subroutine esymm_free_2D
982 : !!***
983 :
984 : !----------------------------------------------------------------------
985 :
986 : !!****f* m_esymm/esymm_finalize
987 : !! NAME
988 : !! esymm_finalize
989 : !!
990 : !! FUNCTION
991 : !!
992 : !! INPUTS
993 : !!
994 : !! SOURCE
995 :
996 0 : subroutine esymm_finalize(esymm, prtvol)
997 :
998 : !Arguments ------------------------------------
999 : !scalars
1000 : class(esymm_t),target,intent(inout) :: esymm
1001 : integer,intent(in) :: prtvol
1002 :
1003 : !Local variables-------------------------------
1004 : integer :: idg,ib1,ib2,idx,nunknown,dg_dim
1005 : integer :: try,irep,nitems,nseen,isn
1006 : integer :: isym,idg1,idg2,dim_mat,irr_idx2,irr_idx1
1007 : real(dp),parameter :: TOL_TRACE=0.1_dp,TOL_ORTHO=0.1_dp,TOL_UNITARY=0.1_dp ! Large tolerance is needed to avoid problems.
1008 : !real(dp),parameter :: TOL_TRACE=0.01_dp,TOL_ORTHO=0.01_dp,TOL_UNITARY=0.01_dp ! Large tolerance is needed to avoid problems.
1009 : !real(dp),parameter :: TOL_TRACE=tol3,TOL_ORTHO=tol3,TOL_UNITARY=tol3 ! Large tolerance is needed to avoid problems.
1010 : real(dp) :: uerr,max_err
1011 : complex(dp) :: ctest
1012 : logical :: isnew
1013 : character(len=500) :: msg
1014 : !arrays
1015 0 : integer,allocatable :: dims_seen(:)
1016 0 : complex(dp),allocatable :: traces_seen(:,:), cidentity(:,:)
1017 0 : complex(dp),pointer :: trace(:), calc_mat(:,:),trace1(:),trace2(:)
1018 : ! *************************************************************************
1019 :
1020 : ! Each band is initialized as "Unknown".
1021 0 : esymm%b2irrep = 0
1022 :
1023 : ! Force the matrices to be unitary.
1024 0 : call polish_irreps(esymm%Calc_irreps)
1025 :
1026 0 : if (.not. esymm%has_chtabs) then
1027 :
1028 : write(msg,'(5a)')&
1029 0 : "Reference character table not available. ",ch10,&
1030 0 : "Symmetry analysis not available. Using heuristic method to classify the states.",ch10,&
1031 0 : "It might not work, especially if accidental degeneracies are present."
1032 0 : ABI_WARNING(msg)
1033 :
1034 : ! The simplest thing we can do here is using the calculated matrices to get the
1035 : ! character and comparing the results hoping everything is OK.
1036 0 : ABI_MALLOC(traces_seen,(esymm%nsym_gk,esymm%ndegs))
1037 0 : ABI_MALLOC(dims_seen,(esymm%ndegs))
1038 :
1039 0 : traces_seen=czero; nseen=1
1040 0 : traces_seen(:,1) = esymm%Calc_irreps(1)%trace
1041 0 : dims_seen(1) = esymm%Calc_irreps(1)%dim
1042 :
1043 0 : do idg=2,esymm%ndegs
1044 0 : dg_dim = esymm%Calc_irreps(idg)%dim
1045 0 : trace => esymm%Calc_irreps(idg)%trace
1046 0 : isnew=.TRUE.
1047 0 : do isn=1,nseen
1048 0 : if (ALL (ABS(trace - traces_seen(:,isn)) < TOL_TRACE) ) then
1049 : isnew=.FALSE.; EXIT
1050 : end if
1051 : end do
1052 :
1053 0 : if (isnew) then
1054 0 : nseen = nseen+1
1055 0 : traces_seen(:,nseen) = trace
1056 0 : dims_seen(nseen) = dg_dim
1057 : end if
1058 : end do
1059 :
1060 0 : if (nseen > esymm%nclass) then
1061 : write(msg, '(3a)') &
1062 0 : "The number of different calculated traces is found to be greater than nclasses!",ch10,&
1063 0 : "Heuristic method clearly failed. Symmetry analysis cannot be performed."
1064 0 : ABI_WARNING(msg)
1065 0 : esymm%err_status = ESYM_HEUR_WRONG_NCLASSES
1066 0 : esymm%err_msg = msg
1067 :
1068 0 : do isn=1,nseen
1069 0 : write(msg,'(a,i0)')" Representation: ",isn
1070 0 : call wrtout(std_out,msg)
1071 0 : call print_arr([std_out], traces_seen(:,isn),max_r=esymm%nsym_gk)
1072 : end do
1073 :
1074 : else ! It seems that the Heuristic method succeeded.
1075 0 : do idg=1,esymm%ndegs
1076 0 : ib1=esymm%degs_bounds(1,idg)
1077 0 : ib2=esymm%degs_bounds(2,idg)
1078 0 : trace => esymm%Calc_irreps(idg)%trace
1079 0 : do isn=1,nseen
1080 0 : if (ALL (ABS(trace - traces_seen(:,isn)) < TOL_TRACE) ) then
1081 0 : esymm%b2irrep(ib1:ib2)=isn
1082 0 : if (esymm%Calc_irreps(idg)%dim /= dims_seen(isn)) then
1083 : write(msg,'(3a)')&
1084 0 : "Found two set of degenerate states with same character but different dimension!",ch10,&
1085 0 : "heuristic method clearly failed. Symmetry analysis cannot be performed."
1086 0 : ABI_ERROR(msg)
1087 0 : esymm%err_status = ESYM_HEUR_WRONG_DIMS
1088 0 : esymm%err_msg = msg
1089 : end if
1090 : EXIT
1091 : end if
1092 : end do
1093 : end do
1094 : end if
1095 :
1096 0 : ABI_FREE(traces_seen)
1097 0 : ABI_FREE(dims_seen)
1098 :
1099 : else
1100 : ! Search in the lookup table definining the irreducible representation
1101 0 : nunknown = 0
1102 0 : do idg=1,esymm%ndegs
1103 :
1104 0 : ib1=esymm%degs_bounds(1,idg)
1105 0 : ib2=esymm%degs_bounds(2,idg)
1106 0 : trace => esymm%Calc_irreps(idg)%trace
1107 :
1108 0 : try = which_irrep(esymm, trace, tol3)
1109 0 : if (try==0) try = which_irrep(esymm, trace, 0.1_dp) ! try again with increased tolerance.
1110 0 : if (try/=0) then
1111 0 : esymm%b2irrep(ib1:ib2)=try
1112 : else
1113 0 : esymm%b2irrep(ib1:ib2)=0
1114 0 : nunknown = nunknown + (ib2-ib1+1)
1115 : end if
1116 : end do
1117 : end if
1118 :
1119 : ! %irrep2b(0)) gives the indices of the states that have not been classified.
1120 0 : ABI_MALLOC(esymm%irrep2b,(0:esymm%nclass))
1121 : !write(std_out,*)"b2irrep",esymm%b2irrep
1122 :
1123 0 : do irep=0,esymm%nclass
1124 0 : nitems = COUNT(esymm%b2irrep==irep)
1125 0 : ABI_MALLOC(esymm%irrep2b(irep)%value,(nitems))
1126 0 : idx=0
1127 0 : do ib1=1,esymm%nbnds
1128 0 : if (esymm%b2irrep(ib1) == irep) then
1129 0 : idx = idx + 1
1130 0 : esymm%irrep2b(irep)%value(idx) = ib1
1131 : end if
1132 : end do
1133 : end do
1134 :
1135 0 : if (size(esymm%irrep2b(0)%value) /= 0) then
1136 0 : write(msg,'(a,i0,a)')" Band classification algorithm was not able to classify ",size(esymm%irrep2b(0)%value)," states."
1137 0 : ABI_WARNING(msg)
1138 0 : esymm%err_status = ESYM_CLASSIFICATION_ERROR
1139 0 : esymm%err_msg = msg
1140 : end if
1141 :
1142 : ! Test basic properties of irreducible representations.
1143 0 : if (.not. esymm%failed()) then
1144 : ! 1) \sum_R \chi^*_a(R)\chi_b(R)= N_R \delta_{ab}
1145 : !call wrtout(std_out," \sum_R \chi^*_a(R)\chi_b(R) = N_R \delta_{ab} ")
1146 0 : max_err=zero
1147 0 : do idg2=1,esymm%ndegs
1148 0 : trace2 => esymm%Calc_irreps(idg2)%trace(1:esymm%nsym_gk)
1149 0 : ib2 = esymm%degs_bounds(1,idg2)
1150 0 : irr_idx2 = esymm%b2irrep(ib2)
1151 0 : if (irr_idx2 == 0) CYCLE
1152 :
1153 0 : do idg1=1,idg2
1154 0 : trace1 => esymm%Calc_irreps(idg1)%trace(1:esymm%nsym_gk)
1155 0 : ib1 = esymm%degs_bounds(1,idg1)
1156 0 : irr_idx1 = esymm%b2irrep(ib1)
1157 0 : if (irr_idx1 == 0) CYCLE
1158 0 : ctest=DOT_PRODUCT(trace1,trace2)/esymm%nsym_gk
1159 0 : if (irr_idx1==irr_idx2) ctest=ctest-one
1160 0 : max_err = MAX(max_err,ABS(ctest))
1161 0 : if (.FALSE..and.ABS(ctest)>tol3) then
1162 : write(msg,'(a,4i3,2es16.8)')&
1163 : ' WARNING: should be delta_ij: cx1 cx2, irr1, irr2, ctest: ',idg1,idg2,irr_idx1,irr_idx2,ctest
1164 : call wrtout(std_out,msg)
1165 : end if
1166 : end do
1167 : end do
1168 :
1169 0 : if (max_err>TOL_ORTHO) then
1170 0 : write(msg,'(a,es10.2)')" Too large maximum error on \sum_R \chi^*_a(R)\chi_b(R) = N_R \delta_{ab}: ",max_err
1171 0 : ABI_WARNING(msg)
1172 0 : esymm%err_status = ESYM_ORTHO_ERROR
1173 0 : esymm%err_msg = msg
1174 : else
1175 0 : write(msg,'(a,es10.2)')" maximum error on \sum_R \chi^*_a(R)\chi_b(R) = N_R \delta_{ab}: ",max_err
1176 0 : call wrtout(std_out,msg)
1177 : end if
1178 :
1179 0 : if (.not.esymm%only_trace) then
1180 : !call wrtout(std_out," **** Testing the unitary of the calculated irreps ****")
1181 0 : max_err=zero
1182 0 : do idg1=1,esymm%ndegs
1183 0 : ib1 = esymm%degs_bounds(1,idg1)
1184 0 : irr_idx1 = esymm%b2irrep(ib1)
1185 0 : if (irr_idx1 == 0) CYCLE
1186 :
1187 0 : do isym=1,esymm%nsym_gk
1188 0 : calc_mat => esymm%Calc_irreps(idg1)%mat(:,:,isym)
1189 0 : dim_mat = esymm%Calc_irreps(idg1)%dim
1190 0 : ABI_MALLOC(cidentity,(dim_mat,dim_mat))
1191 0 : call set2unit(cidentity)
1192 0 : uerr = MAXVAL( ABS(MATMUL(calc_mat,TRANSPOSE(DCONJG(calc_mat))) - cidentity) )
1193 0 : max_err = MAX(max_err,uerr)
1194 0 : ABI_FREE(cidentity)
1195 0 : if (.FALSE..and.prtvol>=10) then
1196 : write(std_out,'(a,i3,a,i2,a,es16.8,a)')&
1197 : " === idg: ",idg1,", isym: ",isym,", Error on U^* U = 1: ",uerr," ==="
1198 : call print_arr([std_out], calc_mat,dim_mat,dim_mat)
1199 : end if
1200 : end do
1201 : end do
1202 :
1203 0 : if (max_err>TOL_UNITARY) then
1204 0 : write(msg,'(a,es10.2)')" Too large maximum error on the unitary of representions matrices: ",max_err
1205 0 : ABI_WARNING(msg)
1206 0 : esymm%err_msg = msg
1207 0 : esymm%err_status = ESYM_UNITARY_ERROR
1208 : else
1209 0 : write(msg,'(a,es10.2)')" maximum error on the unitary of representions matrices: ",max_err
1210 0 : call wrtout(std_out,msg)
1211 : end if
1212 : end if
1213 : end if
1214 :
1215 0 : end subroutine esymm_finalize
1216 : !!***
1217 :
1218 : !----------------------------------------------------------------------
1219 :
1220 : !!****f* m_esymm/which_irrep
1221 : !! NAME
1222 : !! m_esymm
1223 : !!
1224 : !! FUNCTION
1225 : !! Return the index of the irreducible representation with character charact. 0 if not found.
1226 : !!
1227 : !! INPUTS
1228 : !! esymm<esymm_t>
1229 : !! trace(%nsym_gk)=The trace of the representation to be compared with the internal database (if present).
1230 : !! tolerr=Absolute error on the character.
1231 : !!
1232 : !! SOURCE
1233 :
1234 0 : pure integer function which_irrep(esymm, trace, tolerr)
1235 :
1236 : !Arguments ------------------------------------
1237 : !scalars
1238 : class(esymm_t),intent(in) :: esymm
1239 : real(dp),intent(in) :: tolerr
1240 : !arrays
1241 : complex(dp),intent(in) :: trace(esymm%nsym_gk)
1242 :
1243 : !Local variables-------------------------------
1244 : integer :: irp
1245 : ! *********************************************************************
1246 :
1247 0 : which_irrep = 0
1248 0 : if (esymm%has_chtabs) then ! Symmetry analysis can be performed.
1249 0 : do irp=1,esymm%nclass
1250 0 : if (ALL(ABS(esymm%Ref_irreps(irp)%trace(:) - trace(:)) < tolerr)) then
1251 : which_irrep = irp; exit
1252 : end if
1253 : end do
1254 : end if
1255 :
1256 0 : end function which_irrep
1257 : !!***
1258 :
1259 : !----------------------------------------------------------------------
1260 :
1261 : !!****f* m_esymm/esymm_symmetrize_mels
1262 : !! NAME
1263 : !! esymm_symmetrize_mels
1264 : !!
1265 : !! FUNCTION
1266 : !!
1267 : !! INPUTS
1268 : !!
1269 : !! SOURCE
1270 :
1271 0 : subroutine esymm_symmetrize_mels(esymm, lbnd, ubnd, in_me, out_me)
1272 :
1273 : !Arguments ------------------------------------
1274 : !scalars
1275 : class(esymm_t),target,intent(in) :: esymm
1276 : integer :: lbnd,ubnd
1277 : !arrays
1278 : complex(dp),intent(in) :: in_me(2,lbnd:ubnd,lbnd:ubnd)
1279 : complex(dp),intent(out) :: out_me(lbnd:ubnd,lbnd:ubnd)
1280 :
1281 : !Local variables-------------------------------
1282 : !scalars
1283 : integer :: idg1,b1_start,b1_stop,irp1, idg2,b2_start,b2_stop,irp2
1284 : integer :: ii,jj,ib,jb,kk,kb,lb,ll
1285 : complex(dp) :: tr_ofd,ofd,dsd,tr_dsd
1286 : type(irrep_t),pointer :: Irrep1, Irrep2, tr_Irrep1, tr_Irrep2
1287 : ! *********************************************************************
1288 :
1289 0 : if (esymm_failed(esymm)) then
1290 0 : ABI_ERROR("Symmetrization cannot be performed. You should not be here!")
1291 : end if
1292 :
1293 0 : do idg1=1,esymm%ndegs ! First loop over set of degenerate states.
1294 0 : b1_start = esymm%degs_bounds(1,idg1)
1295 0 : b1_stop = esymm%degs_bounds(2,idg1)
1296 :
1297 : !if (b1_stop<lbnd .or. b2_start >ubnd) then
1298 : ! ABI_ERROR("Wrong band indices, check esymm initialization")
1299 : !end if
1300 :
1301 0 : Irrep1 => esymm%Calc_irreps(idg1)
1302 0 : if (esymm%can_use_tr) tr_Irrep1 => esymm%trCalc_irreps(idg1)
1303 0 : irp1 = esymm%b2irrep(b1_start)
1304 :
1305 0 : do idg2=1,esymm%ndegs ! Second loop over set of degenerate states.
1306 : !write(std_out,*)" ==> Symmetrizing degenerate set ",idg1,idg2
1307 0 : b2_start = esymm%degs_bounds(1,idg2)
1308 0 : b2_stop = esymm%degs_bounds(2,idg2)
1309 0 : irp2 = esymm%b2irrep(b2_start)
1310 :
1311 0 : if (irp1/=irp2 .or. idg1==idg2) CYCLE ! Skip diago elements or elements belonging to different irreps.
1312 :
1313 0 : Irrep2 => esymm%Calc_irreps(idg2)
1314 0 : if (esymm%can_use_tr) tr_Irrep2 => esymm%trCalc_irreps(idg2)
1315 : !
1316 : ! Symmetrize the off-diagonal matrix elements.
1317 : ! summing over kk and ll. ii and jj are the indices of the bands that are symmetrized
1318 0 : do ii=1,b1_stop-b1_start+1
1319 0 : ib= ii+b1_start-1
1320 0 : do jj=1,b2_stop-b2_start+1
1321 0 : jb= jj+b2_start-1
1322 : !write(std_out,*)" ====> Symmetrizing ",ib,jb
1323 :
1324 0 : ofd= czero; tr_ofd=czero
1325 0 : do kk=1,b1_stop-b1_start+1
1326 0 : kb= kk+b1_start-1
1327 0 : do ll=1,b2_stop-b2_start+1
1328 0 : lb= ll+b2_start-1
1329 0 : dsd = sum_irreps(Irrep1,Irrep2,kk,ii,ll,jj)
1330 0 : ofd = ofd + dsd * in_me(1,kb,lb)
1331 0 : if (esymm%can_use_tr) then
1332 0 : tr_dsd = sum_irreps(tr_Irrep1,tr_Irrep2,kk,jj,ll,ii) ! Exchange of band indices.
1333 0 : tr_ofd = tr_ofd + tr_dsd * in_me(2,kb,lb) ! Contribution obtained from TR.
1334 : end if
1335 : end do
1336 : end do
1337 :
1338 0 : out_me(ib,jb)= ofd/esymm%nsym_gk
1339 0 : if (esymm%can_use_tr .and. esymm%nsym_trgk>0) out_me(ib,jb)= out_me(ib,jb) + tr_ofd/esymm%nsym_trgk
1340 : end do
1341 : end do
1342 : end do
1343 : end do
1344 :
1345 0 : end subroutine esymm_symmetrize_mels
1346 : !!***
1347 :
1348 : !----------------------------------------------------------------------
1349 :
1350 : !!****f* m_esymm/esymm_failed
1351 : !! NAME
1352 : !! esymm_failed
1353 : !!
1354 : !! FUNCTION
1355 : !!
1356 : !! SOURCE
1357 :
1358 0 : pure logical function esymm_failed(esymm)
1359 :
1360 : !Arguments ------------------------------------
1361 : class(esymm_t),intent(in) :: esymm
1362 : ! *********************************************************************
1363 :
1364 0 : esymm_failed = (esymm%err_status /= ESYM_NOERROR)
1365 :
1366 0 : end function esymm_failed
1367 : !!***
1368 :
1369 : !----------------------------------------------------------------------
1370 :
1371 : !!****f* m_esymm/polish_irreps
1372 : !! NAME
1373 : !! polish_irreps
1374 : !!
1375 : !! FUNCTION
1376 : !!
1377 : !! INPUTS
1378 : !!
1379 : !! SOURCE
1380 :
1381 0 : subroutine polish_irreps(Irreps)
1382 :
1383 : !Arguments ------------------------------------
1384 : type(irrep_t),intent(inout) :: Irreps(:)
1385 :
1386 : !Local variables-------------------------------
1387 : !scalars
1388 : integer,parameter :: ldvl1=1, ldvr1=1
1389 : integer :: irp,sym,dim,ldvr,ii,ivec,jvec,info
1390 : !character(len=500) :: msg
1391 : !arrays
1392 0 : complex(dp),allocatable :: vl(:,:),vr(:,:),vrm1(:,:),overlap(:,:), cmat(:,:),eigval(:)
1393 : ! *********************************************************************
1394 :
1395 : ! Eigen decomposition: A = V D V^{-1}.
1396 0 : do irp=1,SIZE(Irreps)
1397 0 : dim = Irreps(irp)%dim
1398 0 : ABI_MALLOC(cmat,(dim,dim))
1399 0 : ABI_MALLOC(eigval,(dim))
1400 0 : ldvr=dim
1401 0 : ABI_MALLOC(vl,(ldvl1,dim))
1402 0 : ABI_MALLOC(vr,(ldvr,dim))
1403 0 : ABI_MALLOC(vrm1,(dim,dim))
1404 0 : ABI_MALLOC(overlap,(dim,dim))
1405 0 : do sym=1,Irreps(irp)%nsym
1406 0 : cmat = Irreps(irp)%mat(:,:,sym)
1407 0 : call xgeev("No vectors","Vectors",dim,cmat,dim,eigval,vl,ldvl1,vr,ldvr)
1408 : ! Orthogonalize the eigenvectors using Cholesky orthogonalization.
1409 0 : do jvec=1,dim
1410 0 : do ivec=1,jvec
1411 0 : overlap(ivec,jvec) = DOT_PRODUCT(vr(:,ivec),vr(:,jvec))
1412 : end do
1413 : end do
1414 : ! 2) Cholesky factorization: overlap = U^H U with U upper triangle matrix.
1415 0 : call ZPOTRF('U',dim,overlap,dim,info)
1416 0 : ABI_CHECK(info == 0, sjoin('ZPOTRF returned info=', itoa(info)))
1417 :
1418 : ! 3) Solve X U = Vr, on exit the Vr treated by this node is orthonormalized.
1419 0 : call ZTRSM('R','U','N','N',dim,dim,cone,overlap,dim,vr,dim)
1420 : !write(std_out,*)"After ortho",MATMUL(TRANSPOSE(CONJG(vr)),vr)
1421 :
1422 0 : vrm1 = vr
1423 0 : call xginv(vrm1,dim)
1424 0 : do ii=1,dim
1425 0 : eigval(ii) = eigval(ii)/ABS(eigval(ii)) ! Rescale the eigevalues.
1426 0 : vrm1(ii,:) = eigval(ii) * vrm1(ii,:)
1427 : end do
1428 0 : Irreps(irp)%mat(:,:,sym) = MATMUL(vr,vrm1)
1429 0 : Irreps(irp)%trace(sym) = get_trace(Irreps(irp)%mat(:,:,sym))
1430 : end do
1431 0 : ABI_FREE(cmat)
1432 0 : ABI_FREE(eigval)
1433 0 : ABI_FREE(vl)
1434 0 : ABI_FREE(vr)
1435 0 : ABI_FREE(vrm1)
1436 0 : ABI_FREE(overlap)
1437 : end do
1438 :
1439 0 : end subroutine polish_irreps
1440 : !!***
1441 :
1442 : !----------------------------------------------------------------------
1443 :
1444 0 : end module m_esymm
|
