Branch data 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 : 1247 : 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 [ - + ]: 1247 : ABI_SFREE(esymm%g0)
923 [ - + ]: 1247 : ABI_SFREE(esymm%tr_g0)
924 [ - + ]: 1247 : ABI_SFREE(esymm%nelements)
925 [ - + ]: 1247 : ABI_SFREE(esymm%sgk2symrec)
926 [ - + ]: 1247 : ABI_SFREE(esymm%tr_sgk2symrec)
927 [ - + ]: 1247 : ABI_SFREE(esymm%herring_test)
928 [ - + ]: 1247 : ABI_SFREE(esymm%b2irrep)
929 [ - + ]: 1247 : ABI_SFREE(esymm%degs_bounds)
930 [ - + ]: 1247 : ABI_SFREE(esymm%degs_dim)
931 : :
932 [ - + ]: 1247 : 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 [ - + ]: 1247 : 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 [ - + ]: 1247 : 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 [ - + ]: 1247 : 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 : 1247 : 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 : 206 : 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 [ + + ]: 416 : do id2=1,SIZE(esymm,DIM=2)
976 [ + + ]: 1663 : do id1=1,SIZE(esymm,DIM=1)
977 : 1457 : call esymm_free_0D(esymm(id1,id2))
978 : : end do
979 : : end do
980 : :
981 : 206 : 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
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
|
