LCOV - code coverage report
Current view: top level - src/43_ptgroups - m_ptgroups.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 0.0 % 186 0
Test Date: 2026-09-21 19:39:32 Functions: 0.0 % 7 0

            Line data    Source code
       1              : !!****m* ABINIT/m_ptgroups
       2              : !! NAME
       3              : !! m_ptgroups
       4              : !!
       5              : !! FUNCTION
       6              : !!  This module contains the irreducible representations and the
       7              : !!  character tables of the 32 point groups.
       8              : !!
       9              : !! COPYRIGHT
      10              : !! Copyright (C) 2010-2026 ABINIT group (MG)
      11              : !! This file is distributed under the terms of the
      12              : !! GNU General Public License, see ~abinit/COPYING
      13              : !! or http://www.gnu.org/copyleft/gpl.txt .
      14              : !!
      15              : 
      16              : #if defined HAVE_CONFIG_H
      17              : #include "config.h"
      18              : #endif
      19              : 
      20              : #include "abi_common.h"
      21              : 
      22              : module m_ptgroups
      23              : 
      24              :  use defs_basis
      25              :  use m_errors
      26              :  use m_abicore
      27              :  use m_defs_ptgroups
      28              : 
      29              :  use m_io_tools,       only : open_file
      30              :  use m_fstrings,       only : sjoin
      31              :  use m_numeric_tools,  only : get_trace, cmplx_sphcart
      32              :  use m_matrix,         only : mati3inv
      33              :  use m_symtk,          only : sg_multable
      34              : 
      35              :  ! Import group tables
      36              :  use m_ptg_C1
      37              :  use m_ptg_Ci
      38              :  use m_ptg_C2
      39              :  use m_ptg_Cs
      40              :  use m_ptg_C2h
      41              :  use m_ptg_D2
      42              :  use m_ptg_C2v
      43              :  use m_ptg_D2h
      44              :  use m_ptg_C4
      45              :  use m_ptg_S4
      46              :  use m_ptg_C4h
      47              :  use m_ptg_D4
      48              :  use m_ptg_C4v
      49              :  use m_ptg_D2d
      50              :  use m_ptg_D4h
      51              :  use m_ptg_C3
      52              :  use m_ptg_C3i
      53              :  use m_ptg_D3
      54              :  use m_ptg_C3v
      55              :  use m_ptg_D3d
      56              :  use m_ptg_C6
      57              :  use m_ptg_C3h
      58              :  use m_ptg_C6h
      59              :  use m_ptg_D6
      60              :  use m_ptg_C6v
      61              :  use m_ptg_D3h
      62              :  use m_ptg_D6h
      63              :  use m_ptg_T
      64              :  use m_ptg_Th
      65              :  use m_ptg_O
      66              :  use m_ptg_Td
      67              :  use m_ptg_Oh
      68              : 
      69              :  implicit none
      70              : 
      71              :  private
      72              : 
      73              :  public :: get_point_group
      74              :  public :: get_classes
      75              :  public :: show_character_tables
      76              :  public :: point_group_init
      77              :  public :: sum_irreps
      78              :  public :: mult_table
      79              :  !public :: polish_irreps
      80              : 
      81              : 
      82              : contains
      83              : !!***
      84              : 
      85              : !----------------------------------------------------------------------
      86              : 
      87              : !!****f* m_ptgroups/get_point_group
      88              : !! NAME
      89              : !!  get_point_group
      90              : !!
      91              : !! FUNCTION
      92              : !!
      93              : !! INPUTS
      94              : !! ptg_name=point group name as returned by symptgroup.
      95              : !!
      96              : !! OUTPUT
      97              : !! nsym=Number of symmetries in the point group.
      98              : !! nclass=Number of classes.
      99              : !! sym(3,3,nsym)=Elements of the point group ordered by classe.
     100              : !! class_ids(2,nclass)=Initial and final index in sym, for each
     101              : !! Irreps(nclass)=Datatype gathering data on the different irreducible representations.
     102              : !!
     103              : !! SOURCE
     104              : 
     105            0 : subroutine get_point_group(ptg_name, nsym, nclass, sym, class_ids, class_names, Irreps)
     106              : 
     107              : !Arguments ------------------------------------
     108              : !scalars
     109              :  character(len=*),intent(in) :: ptg_name
     110              :  integer,intent(out) :: nsym, nclass
     111              : !arrays
     112              :  integer,allocatable,intent(out) :: sym(:,:,:),class_ids(:,:)
     113              :  character(len=5),allocatable,intent(out) :: class_names(:)
     114              :  type(irrep_t),allocatable,intent(out) :: Irreps(:)
     115              : 
     116              : !Local variables-------------------------------
     117              :  integer :: irp,isym
     118              : ! *************************************************************************
     119              : 
     120            0 :  SELECT CASE (TRIM(ADJUSTL(ptg_name)))
     121              :  CASE ('1')
     122            0 :    call ptg_C1  (nsym,nclass,sym,class_ids,class_names,Irreps)
     123              :  CASE ('-1')
     124            0 :    call ptg_Ci  (nsym,nclass,sym,class_ids,class_names,Irreps)
     125              :  CASE ('2')
     126            0 :    call ptg_C2  (nsym,nclass,sym,class_ids,class_names,Irreps)
     127              :  CASE ('m',"-2") ! Abinit uses "-2"
     128            0 :    call ptg_Cs  (nsym,nclass,sym,class_ids,class_names,Irreps)
     129              :  CASE ('2/m')
     130            0 :    call ptg_C2h (nsym,nclass,sym,class_ids,class_names,Irreps)
     131              :  CASE ('222')
     132            0 :    call ptg_D2  (nsym,nclass,sym,class_ids,class_names,Irreps)
     133              :  CASE ('mm2')
     134            0 :    call ptg_C2v (nsym,nclass,sym,class_ids,class_names,Irreps)
     135              :  CASE ('mmm')
     136            0 :    call ptg_D2h (nsym,nclass,sym,class_ids,class_names,Irreps)
     137              :  CASE ('4')
     138            0 :    call ptg_C4  (nsym,nclass,sym,class_ids,class_names,Irreps)
     139              :  CASE ('-4')
     140            0 :    call ptg_S4  (nsym,nclass,sym,class_ids,class_names,Irreps)
     141              :  CASE ('4/m')
     142            0 :    call ptg_C4h (nsym,nclass,sym,class_ids,class_names,Irreps)
     143              :  CASE ('422')
     144            0 :    call ptg_D4  (nsym,nclass,sym,class_ids,class_names,Irreps)
     145              :  CASE ('4mm')
     146            0 :    call ptg_C4v (nsym,nclass,sym,class_ids,class_names,Irreps)
     147              :  CASE ('-42m')
     148            0 :    call ptg_D2d (nsym,nclass,sym,class_ids,class_names,Irreps)
     149              :  CASE ('4/mmm')
     150            0 :    call ptg_D4h (nsym,nclass,sym,class_ids,class_names,Irreps)
     151              :  CASE ('3')
     152            0 :    call ptg_C3  (nsym,nclass,sym,class_ids,class_names,Irreps)
     153              :  CASE ('-3')
     154            0 :    call ptg_C3i (nsym,nclass,sym,class_ids,class_names,Irreps)
     155              :  CASE ('32')
     156            0 :    call ptg_D3  (nsym,nclass,sym,class_ids,class_names,Irreps)
     157              :  CASE ('3m')
     158            0 :    call ptg_C3v (nsym,nclass,sym,class_ids,class_names,Irreps)
     159              :  CASE ('-3m')
     160            0 :    call ptg_D3d (nsym,nclass,sym,class_ids,class_names,Irreps)
     161              :  CASE ('6')
     162            0 :    call ptg_C6  (nsym,nclass,sym,class_ids,class_names,Irreps)
     163              :  CASE ('-6')
     164            0 :    call ptg_C3h (nsym,nclass,sym,class_ids,class_names,Irreps)
     165              :  CASE ('6/m')
     166            0 :    call ptg_C6h (nsym,nclass,sym,class_ids,class_names,Irreps)
     167              :  CASE ('622')
     168            0 :    call ptg_D6  (nsym,nclass,sym,class_ids,class_names,Irreps)
     169              :  CASE ('6mm')
     170            0 :    call ptg_C6v (nsym,nclass,sym,class_ids,class_names,Irreps)
     171              :  CASE ('-62m')
     172            0 :    call ptg_D3h (nsym,nclass,sym,class_ids,class_names,Irreps)
     173              :  CASE ('6/mmm')
     174            0 :    call ptg_D6h (nsym,nclass,sym,class_ids,class_names,Irreps)
     175              :  CASE ('23')
     176            0 :    call ptg_T   (nsym,nclass,sym,class_ids,class_names,Irreps)
     177              :  CASE ('m-3')
     178            0 :    call ptg_Th  (nsym,nclass,sym,class_ids,class_names,Irreps)
     179              :  CASE ('432')
     180            0 :    call ptg_O   (nsym,nclass,sym,class_ids,class_names,Irreps)
     181              :  CASE ('-43m')
     182            0 :    call ptg_Td  (nsym,nclass,sym,class_ids,class_names,Irreps)
     183              :  CASE ('m-3m')
     184            0 :    call ptg_Oh  (nsym,nclass,sym,class_ids,class_names,Irreps)
     185              :  CASE DEFAULT
     186            0 :    ABI_BUG(sjoin("Unknown value for ptg_name:", ptg_name))
     187              :  END SELECT
     188              : 
     189              :  ! Calculate the trace of each irreducible representation in order to have the character at hand.
     190            0 :  do irp=1,SIZE(Irreps)
     191            0 :    ABI_MALLOC(Irreps(irp)%trace, (nsym))
     192            0 :    do isym=1,nsym
     193            0 :      Irreps(irp)%trace(isym) = get_trace(Irreps(irp)%mat(:,:,isym))
     194              :    end do
     195              :  end do
     196              : 
     197            0 : end subroutine get_point_group
     198              : !!***
     199              : 
     200              : !----------------------------------------------------------------------
     201              : 
     202              : !!****f* m_ptgroups/get_classes
     203              : !! NAME
     204              : !! get_classes
     205              : !!
     206              : !! FUNCTION
     207              : !!  Given a set of nsym 3x3 operations in reciprocal or real space,
     208              : !!  which are supposed to form a group, this routine divides the group into classes.
     209              : !!
     210              : !! INPUTS
     211              : !! nsym=number of symmetry operation.
     212              : !! sym(3,3,nsym)=The symmetry operations.
     213              : !!
     214              : !! OUTPUT
     215              : !! nclass=The number of classes
     216              : !! nelements(1:nclass)=For each class, the number of elements
     217              : !! elements_idx(ii,1:nclass)=For each class, this table gives the index of its elements (ii=1,..,nelements(iclass))
     218              : !!
     219              : !! NOTES
     220              : !!  * A class is defined as the set of distinct elements obtained by
     221              : !!    considering for each element, S, of the group all its conjugate
     222              : !!    elements X^-1 S X where X ranges over all the elements of the group.
     223              : !!
     224              : !!  * It does not work in case of non-collinear magnetism.
     225              : !!
     226              : !!  * The routine assumes that anti-ferromagnetic symmetries (if any) have been removed by the caller.
     227              : !!
     228              : !! SOURCE
     229              : 
     230            0 : subroutine get_classes(nsym, sym, nclass, nelements, elements_idx)
     231              : 
     232              : !Arguments ------------------------------------
     233              : !scalars
     234              :  integer,intent(in) :: nsym
     235              :  integer,intent(out) :: nclass
     236              : !arrays
     237              :  integer,intent(in) :: sym(3,3,nsym)
     238              :  integer,intent(out) :: nelements(nsym),elements_idx(nsym,nsym)
     239              : 
     240              : !Local variables-------------------------------
     241              : !scalars
     242              :  integer :: isym,jsym,ksym,identity_idx,ierr
     243              :  character(len=500) :: msg
     244              : !arrays
     245            0 :  integer :: cjg(3,3),ss(3,3),xx(3,3),xxm1(3,3),test(3,3), identity(3,3), dummy_symafm(nsym)
     246            0 :  logical :: found(nsym),found_identity
     247              : !************************************************************************
     248              : 
     249              :  ! Check if identity is present in the first position
     250            0 :  identity=RESHAPE((/1,0,0,0,1,0,0,0,1/),(/3,3/)); found_identity=.FALSE.
     251              : 
     252            0 :  do isym=1,nsym
     253            0 :    if (all(sym(:,:,isym) == identity)) then
     254            0 :      found_identity=.TRUE.; identity_idx=isym; exit
     255              :    end if
     256              :  end do
     257              : 
     258            0 :  if (.not. found_identity .or. identity_idx /= 1) then
     259              :    write(msg, '(3a)')&
     260            0 :     'Either identity is not present or it is not the first operation ',ch10,'Check set of symmetry operations '
     261            0 :    ABI_ERROR(msg)
     262              :  end if
     263              : 
     264            0 :  dummy_symafm = 1
     265            0 :  call sg_multable(nsym, dummy_symafm, sym, ierr)
     266            0 :  ABI_CHECK(ierr == 0, "Error in group closure")
     267              : 
     268            0 :  nclass=0; nelements(:)=0; elements_idx(:,:)=0; found(:)=.FALSE.
     269            0 :  do isym=1,nsym
     270            0 :    if (.not.found(isym)) then
     271            0 :      nclass=nclass+1
     272            0 :      ss(:,:)=sym(:,:,isym)
     273              : 
     274            0 :      do jsym=1,nsym ! Form conjugate.
     275            0 :        xx(:,:)=sym(:,:,jsym)
     276            0 :        call mati3inv(xx,xxm1) ; xxm1=TRANSPOSE(xxm1)
     277            0 :        cjg(:,:)=MATMUL(xxm1,MATMUL(ss,xx))
     278            0 :        do ksym=1,nsym ! Is it already found?
     279            0 :          test(:,:)=sym(:,:,ksym)
     280            0 :          if (.not.found(ksym).and.(ALL((test-cjg)==0))) then
     281            0 :            found(ksym)=.TRUE.
     282            0 :            nelements(nclass)=nelements(nclass)+1
     283            0 :            elements_idx(nelements(nclass),nclass)=ksym
     284              :          end if
     285              :        end do
     286              :      end do
     287              : 
     288              :    end if
     289              :  end do
     290              : 
     291            0 : end subroutine get_classes
     292              : !!***
     293              : 
     294              : !----------------------------------------------------------------------
     295              : 
     296              : !!****f* m_ptgroups/show_character_tables
     297              : !! NAME
     298              : !!  show_character_tables
     299              : !!
     300              : !! FUNCTION
     301              : !!   Printout of the caracter tables of the 32 point groups.
     302              : !!
     303              : !! INPUTS
     304              : !!  [unit]=Unit number of output file. Defaults to std_out
     305              : !!
     306              : !! OUTPUT
     307              : !!  Only writing.
     308              : !!
     309              : !! SOURCE
     310              : 
     311            0 : subroutine show_character_tables(unit)
     312              : 
     313              : !Arguments ------------------------------------
     314              : !scalars
     315              :  integer,optional,intent(in) :: unit
     316              : 
     317              : !Local variables-------------------------------
     318              :  integer :: igrp,my_unt
     319              :  character(len=5) :: ptg_name
     320            0 :  type(point_group_t) :: Ptg
     321              : !arrays
     322              :  !integer,allocatable :: elements_idx(:,:),nelements(:)
     323              : ! *********************************************************************
     324              : 
     325            0 :  my_unt = std_out; if (PRESENT(unit)) my_unt=unit
     326              : 
     327            0 :  do igrp=1,SIZE(ptgroup_names)
     328            0 :    ptg_name = ptgroup_names(igrp)
     329            0 :    call point_group_init(Ptg, ptg_name)
     330            0 :    call Ptg%print(unit=my_unt)
     331              :    !allocate(nelements(Ptg%nsym),elements_idx(Ptg%nsym,Ptg%nsym))
     332              :    !call get_classes(Ptg%nsym,Ptg%sym,nclass,nelements,elements_idx)
     333              :    !deallocate(nelements,elements_idx)
     334            0 :    call Ptg%free()
     335              :  end do
     336              : 
     337            0 : end subroutine show_character_tables
     338              : !!***
     339              : 
     340              : !----------------------------------------------------------------------
     341              : 
     342              : !!****f* m_ptgroups/point_group_init
     343              : !! NAME
     344              : !! point_group_init
     345              : !!
     346              : !! FUNCTION
     347              : !!  Creation method for the point_group_t datatype.
     348              : !!
     349              : !! INPUTS
     350              : !!  ptg_name=The name of the point group (International conventions).
     351              : !!
     352              : !! OUTPUT
     353              : !!  The datatype completely initialized.
     354              : !!
     355              : !! SOURCE
     356              : 
     357            0 : subroutine point_group_init(Ptg, ptg_name)
     358              : 
     359              : !Arguments ------------------------------------
     360              :  class(point_group_t),intent(inout) :: Ptg
     361              :  character(len=5),intent(in) :: ptg_name
     362              : ! *********************************************************************
     363              : 
     364            0 :  Ptg%gname = ptg_name
     365            0 :  call get_point_group(Ptg%gname,Ptg%nsym,Ptg%nclass,Ptg%sym,Ptg%class_ids,Ptg%class_names,Ptg%Irreps)
     366              : 
     367            0 : end subroutine point_group_init
     368              : !!***
     369              : 
     370              : !----------------------------------------------------------------------
     371              : 
     372              : !!****f* m_ptgroups/mult_table
     373              : !! NAME
     374              : !! mult_table
     375              : !!
     376              : !! FUNCTION
     377              : !!  Given a set of nsym 3x3 operations which are supposed to form a group,
     378              : !!  this routine constructs the multiplication table of the group.
     379              : !!
     380              : !! INPUTS
     381              : !! nsym=number of symmetry operation
     382              : !! sym(3,3,nsym)=the operations
     383              : !!
     384              : !! OUTPUT
     385              : !!  mtab(nsym,nsym)=The index of the product S_i * S_j in the input set sym.
     386              : !!
     387              : !! SOURCE
     388              : 
     389            0 : subroutine mult_table(nsym, sym, mtab)
     390              : 
     391              : !Arguments ------------------------------------
     392              : !scalars
     393              :  integer,intent(in) :: nsym
     394              : !arrays
     395              :  integer,intent(in) :: sym(3,3,nsym)
     396              :  integer,intent(out) :: mtab(nsym,nsym)
     397              : 
     398              : !Local variables-------------------------------
     399              : !scalars
     400              :  integer :: isym,jsym,ksym
     401              :  !character(len=500) :: msg
     402              : !arrays
     403            0 :  integer :: prod_ij(3,3),found(nsym)
     404              : !************************************************************************
     405              : 
     406            0 :  do jsym=1,nsym
     407            0 :    found(:)=0 ! Each symmetry should compare only once in a given (row|col).
     408              : 
     409            0 :    do isym=1,nsym
     410            0 :      prod_ij = MATMUL(sym(:,:,isym),sym(:,:,jsym))
     411            0 :      do ksym=1,nsym
     412            0 :        if ( ALL(prod_ij == sym(:,:,ksym)) ) then
     413            0 :          found(ksym)=found(ksym)+1
     414            0 :          mtab(isym,jsym) = ksym
     415              :        end if
     416              :      end do
     417              :    end do ! jsym
     418              : 
     419            0 :    if (any(found /= 1)) then
     420            0 :      write(std_out,*)"found = ",found
     421            0 :      ABI_ERROR("Input elements do not form a group")
     422              :    end if
     423              :  end do ! isym
     424              : 
     425            0 : end subroutine mult_table
     426              : !!***
     427              : 
     428              : !----------------------------------------------------------------------
     429              : 
     430              : !!****f* m_ptgroups/groupk_from_file
     431              : !! NAME
     432              : !!  groupk_from_file
     433              : !!
     434              : !! FUNCTION
     435              : !!  Initialize the group_k_t datatype from an external database retrieved from
     436              : !!  the Bilbao server via the ptg.py script.
     437              : !!
     438              : !! INPUTS
     439              : !!  fname(len=*)=file name
     440              : !!
     441              : !! OUTPUT
     442              : !!  ierr=Status error
     443              : !!  Lgrps<group_k_t>=The structure completely initialized.
     444              : !!
     445              : !! TODO
     446              : !!   This is a stub. I still have to complete the fileformat for the Bilbao database.
     447              : !!
     448              : !! SOURCE
     449              : 
     450            0 : subroutine groupk_from_file(Lgrps, spgroup, fname, nkpt, klist, ierr)
     451              : 
     452              : !Arguments ------------------------------------
     453              : !scalars
     454              :  integer,intent(in) :: spgroup
     455              :  integer,intent(out) :: ierr,nkpt
     456              :  character(len=fnlen),intent(in) :: fname
     457              : !arrays
     458              :  type(group_k_t),target,allocatable :: Lgrps(:)
     459              :  real(dp),pointer :: klist(:,:)
     460              : 
     461              : !Local variables-------------------------------
     462              : !scalars
     463              :  integer,parameter :: last_file_version=1
     464              :  integer :: unt,fvers,ik,nsym_ltgk,isym,nirreps_k,irp,icls
     465              :  integer :: irrep_idx,irrep_dim,sym_idx,ita_spgroup,nels,prev,now
     466              :  character(len=IRREPNAME_LEN) :: irrep_name
     467              :  character(len=1) :: basis
     468              :  character(len=500) :: msg
     469              :  type(group_k_t),pointer :: Gk
     470              : !arrays
     471            0 :  integer,allocatable :: nelements(:),elements_idx(:,:)
     472              :  real(dp) :: kpt(3)
     473            0 :  character(len=10),allocatable :: kname(:)
     474              :  type(irrep_t),pointer :: OneIrr
     475              : ! *************************************************************************
     476              : 
     477            0 :  ierr=0
     478            0 :  if (open_file(fname,msg,newunit=unt,form="formatted") /=0) then
     479            0 :    ABI_ERROR(msg)
     480              :  end if
     481              : 
     482            0 :  read(unt,*,ERR=10)          ! Skip the header.
     483            0 :  read(unt,*,ERR=10) fvers    ! File version.
     484            0 :  if (fvers > last_file_version) then
     485            0 :    write(msg,"(2(a,i0))")" Found file format= ",fvers," but the latest supported version is: ",last_file_version
     486            0 :    ABI_ERROR(msg)
     487              :  end if
     488            0 :  read(unt,*,ERR=10) ita_spgroup
     489            0 :  read(unt,*,ERR=10) basis
     490              : 
     491            0 :  if (spgroup/=ita_spgroup) then
     492            0 :    write(msg,'(a,2i0)')" Input space group does not match with the value reported on file: ",spgroup,ita_spgroup
     493            0 :    ABI_ERROR(msg)
     494              :  end if
     495              : 
     496            0 :  if (basis /= "b") then
     497            0 :    ABI_ERROR(" Wrong value for basis: "//TRIM(basis))
     498              :  end if
     499              : 
     500              :  ! Read the list of the k-points.
     501            0 :  read(unt,*,ERR=10) nkpt
     502              : 
     503            0 :  ABI_MALLOC(Lgrps,(nkpt))
     504            0 :  ABI_MALLOC(klist,(3,nkpt))
     505            0 :  ABI_MALLOC(kname,(nkpt))
     506            0 :  do ik=1,nkpt
     507            0 :    read(unt,*,ERR=10) klist(:,ik), kname(ik)
     508              :  end do
     509              : 
     510              :  ! Read tables for each k-point
     511            0 :  do ik=1,nkpt
     512              : 
     513            0 :    read(unt,*,ERR=10) kpt
     514            0 :    read(unt,*,ERR=10) nsym_ltgk
     515            0 :    Gk  => Lgrps(ik)
     516              : 
     517            0 :    Gk%spgroup = ita_spgroup
     518            0 :    Gk%nsym    = nsym_ltgk
     519            0 :    Gk%point   = kpt
     520            0 :    ABI_MALLOC(Gk%sym,(3,3,nsym_ltgk))
     521            0 :    ABI_MALLOC(Gk%tnons,(3,nsym_ltgk))
     522              : 
     523            0 :    do isym=1,nsym_ltgk ! Read symmetries of the little group.
     524            0 :      read(unt,*,ERR=10) Gk%sym(:,:,isym)
     525            0 :      read(unt,*,ERR=10) Gk%tnons(:,isym)
     526              :    end do
     527              : 
     528            0 :    ABI_MALLOC(nelements,(nsym_ltgk))
     529            0 :    ABI_MALLOC(elements_idx,(nsym_ltgk,nsym_ltgk))
     530              : 
     531            0 :    call get_classes(nsym_ltgk,Gk%sym,Gk%nclass,nelements,elements_idx)
     532              : 
     533              :    ! The operations reported on the file are supposed to be packed in classes
     534              :    ! otherwise one should perform a rearrangement of the indices.
     535            0 :    prev = 0
     536            0 :    do icls=1,Gk%nclass
     537            0 :      do isym=1,nelements(icls)
     538            0 :        now = elements_idx(isym,icls)
     539            0 :        if ( (now-prev) /= 1 ) then
     540            0 :          write(msg,"(2(a,i0))")" Symmetries on file are not ordered in classes. icls= ",icls,", isym= ",isym
     541            0 :          ABI_ERROR(msg)
     542              :        else
     543            0 :          prev = now
     544              :        end if
     545              :      end do
     546              :    end do
     547              : 
     548            0 :    ABI_MALLOC(Gk%class_ids,(2,Gk%nclass))
     549            0 :    do icls=1,Gk%nclass
     550            0 :      nels = nelements(icls)
     551            0 :      Gk%class_ids(1,icls) = elements_idx(1,   icls)
     552            0 :      Gk%class_ids(2,icls) = elements_idx(nels,icls)
     553              :    end do
     554              : 
     555            0 :    ABI_FREE(nelements)
     556            0 :    ABI_FREE(elements_idx)
     557              : 
     558              :    ! Read the irreducible representations.
     559            0 :    read(unt,*,ERR=10) nirreps_k
     560            0 :    ABI_CHECK(Gk%nclass == nirreps_k,"Gk%nclass /= nirreps_k")
     561              : 
     562              :    !$$ allocate(Gk%class_names(Gk%nclass))
     563            0 :    ABI_MALLOC(Gk%Irreps, (nirreps_k))
     564              : 
     565            0 :    do irp=1,nirreps_k
     566            0 :      OneIrr =>  Gk%Irreps(irp)
     567            0 :      read(unt,*,ERR=10) irrep_idx, irrep_dim, irrep_name
     568            0 :      call OneIrr%init(nsym_ltgk, irrep_dim, irrep_name)
     569            0 :      do isym=1,nsym_ltgk
     570            0 :        read(unt,*,ERR=10) sym_idx, OneIrr%mat(:,:,isym)
     571            0 :        ABI_CHECK(sym_idx==irp,"sym_idx/=irp!")
     572              :        ! Matrix elements on file are in the form (rho, theta) with theta given in degrees.
     573            0 :        call cmplx_sphcart(OneIrr%mat(:,:,isym),from="Sphere",units="Degrees")
     574            0 :        OneIrr%trace(isym) = get_trace(OneIrr%mat(:,:,isym))
     575              :      end do
     576              :    end do
     577              : 
     578              :  end do
     579              : 
     580            0 :  close(unt)
     581            0 :  RETURN
     582              : 
     583              :  ! Handle IO-error.
     584            0 : 10 ierr=1
     585            0 :  close(unt)
     586            0 :  RETURN
     587              : 
     588            0 : end subroutine groupk_from_file
     589              : !!***
     590              : 
     591              : !----------------------------------------------------------------------
     592              : 
     593              : !!****f* m_ptgroups/sum_irreps
     594              : !! NAME
     595              : !!  sum_irreps
     596              : !!
     597              : !! FUNCTION
     598              : !!
     599              : !! INPUTS
     600              : !!
     601              : !! OUTPUT
     602              : !!
     603              : !! SOURCE
     604              : 
     605            0 : function sum_irreps(Irrep1,Irrep2,ii,jj,kk,ll) result(res)
     606              : 
     607              : !Arguments ------------------------------------
     608              : !scalars
     609              :  class(irrep_t),intent(in) :: Irrep1,Irrep2
     610              :  integer,intent(in) :: ii,jj,kk,ll
     611              : !arrays
     612              :  complex(dp) :: res
     613              : 
     614              : !Local variables-------------------------------
     615              :  integer :: isym,nsym,ierr
     616              :  !character(len=500) :: msg
     617              : ! *********************************************************************
     618              : 
     619            0 :  ierr=0; res = czero
     620              : 
     621            0 :  nsym = Irrep1%nsym
     622            0 :  if (nsym /= Irrep2%nsym) then
     623            0 :    ABI_WARNING("Irreps have different nsym")
     624            0 :    ierr=ierr+1
     625              :  end if
     626              : 
     627            0 :  if (Irrep1%dim /= Irrep2%dim) then
     628            0 :    ABI_WARNING("Irreps have different dimensions")
     629            0 :    write(std_out,*)Irrep1%dim,Irrep2%dim
     630            0 :    ierr=ierr+1
     631              :  end if
     632              : 
     633            0 :  if (ii > Irrep2%dim .or. jj > Irrep2%dim .or. kk > Irrep1%dim .or. ll > Irrep1%dim) then
     634            0 :    ABI_WARNING("Wrong indices")
     635            0 :    write(std_out,*)ii,Irrep2%dim,jj,Irrep2%dim,kk>Irrep1%dim,ll,Irrep1%dim
     636            0 :    ierr=ierr+1
     637              :  end if
     638              : 
     639            0 :  if (ierr/=0) RETURN
     640              : 
     641            0 :  do isym=1,nsym
     642            0 :    res = res + DCONJG(Irrep1%mat(ii,jj,isym)) * Irrep2%mat(kk,ll,isym)
     643              :  end do
     644              : 
     645            0 : end function sum_irreps
     646              : !!***
     647              : 
     648              : end module m_ptgroups
     649              : !!***
        

Generated by: LCOV version 2.3-1