TABLE OF CONTENTS


ABINIT/m_crystal [ Modules ]

[ Top ] [ Modules ]

NAME

 m_crystal

FUNCTION

 Module containing the definition of the crystal_t data type and methods used to handle it.

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (MG, YP, MJV)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 MODULE m_crystal
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32  use m_atomdata
33 
34  use m_numeric_tools,  only : set2unit
35  use m_symtk,          only : mati3inv, sg_multable, symatm, print_symmetries
36  use m_spgdata,        only : spgdata
37  use m_geometry,       only : metric, xred2xcart, remove_inversion, getspinrot
38  use m_io_tools,       only : open_file
39  use m_fstrings,       only : int2char10
40 
41  implicit none
42 
43  private

m_crystal/adata_type [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

  adata_type

FUNCTION

  Return atomic data from the itypat index

PARENTS

SOURCE

830 function adata_type(crystal,itypat) result(atom)
831 
832 
833 !This section has been created automatically by the script Abilint (TD).
834 !Do not modify the following lines by hand.
835 #undef ABI_FUNC
836 #define ABI_FUNC 'adata_type'
837 !End of the abilint section
838 
839  implicit none
840 
841 !Arguments ------------------------------------
842 !scalars
843  integer,intent(in) :: itypat
844  type(crystal_t),intent(in) :: crystal
845  type(atomdata_t) :: atom
846 
847 ! *************************************************************************
848 
849  call atomdata_from_znucl(atom,crystal%znucl(itypat))
850 
851 end function adata_type

m_crystal/crystal_free [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

  crystal_free

FUNCTION

  Destroy the dynamic arrays in a crystal_t data type.

PARENTS

      anaddb,bethe_salpeter,cut3d,dfpt_looppert,eig2tot,eph,fold2Bloch,gstate
      gwls_hamiltonian,m_ddk,m_dvdb,m_effective_potential
      m_effective_potential_file,m_gruneisen,m_ioarr,m_iowf,m_wfd,m_wfk
      mlwfovlp_qp,mover,mrgscr,optic,outscfcv,respfn,screening,sigma,vtorho
      wfk_analyze

CHILDREN

      mati3inv,sg_multable

SOURCE

435 subroutine crystal_free(Cryst)
436 
437 
438 !This section has been created automatically by the script Abilint (TD).
439 !Do not modify the following lines by hand.
440 #undef ABI_FUNC
441 #define ABI_FUNC 'crystal_free'
442 !End of the abilint section
443 
444  implicit none
445 
446 !Arguments ------------------------------------
447  type(crystal_t),intent(inout) :: Cryst
448 
449 ! *********************************************************************
450 
451  DBG_ENTER("COLL")
452 
453  !@crystal_t
454 
455 !integer
456  if (allocated(Cryst%indsym))  then
457    ABI_FREE(Cryst%indsym)
458  end if
459  if (allocated(Cryst%symafm))  then
460    ABI_FREE(Cryst%symafm)
461  end if
462  if (allocated(Cryst%symrec))  then
463    ABI_FREE(Cryst%symrec)
464  end if
465  if (allocated(Cryst%symrel))  then
466    ABI_FREE(Cryst%symrel)
467  end if
468  if (allocated(Cryst%atindx))  then
469    ABI_FREE(Cryst%atindx)
470  end if
471  if (allocated(Cryst%atindx1))  then
472    ABI_FREE(Cryst%atindx1)
473  end if
474  if (allocated(Cryst%typat  ))  then
475    ABI_FREE(Cryst%typat)
476  end if
477  if (allocated(Cryst%nattyp))  then
478    ABI_FREE(Cryst%nattyp)
479  end if
480 
481 !real
482  if (allocated(Cryst%tnons))  then
483    ABI_FREE(Cryst%tnons)
484  end if
485  if (allocated(Cryst%xcart))  then
486    ABI_FREE(Cryst%xcart)
487  end if
488  if (allocated(Cryst%xred))  then
489    ABI_FREE(Cryst%xred)
490  end if
491  if (allocated(Cryst%zion))  then
492    ABI_FREE(Cryst%zion)
493  end if
494  if (allocated(Cryst%znucl))  then
495    ABI_FREE(Cryst%znucl)
496  end if
497  if (allocated(Cryst%amu))  then
498    ABI_FREE(Cryst%amu)
499  end if
500  if (allocated(Cryst%spinrot)) then
501     ABI_FREE(Cryst%spinrot)
502  end if
503 
504 !character
505  if (allocated(Cryst%title))  then
506    ABI_FREE(Cryst%title)
507  end if
508 
509  DBG_EXIT("COLL")
510 
511 end subroutine crystal_free

m_crystal/crystal_init [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

  crystal_init

FUNCTION

  Initialize a crystal_t data type.
  Ideally the routine should work in two different modes:
  Either the symmetries are directly supplied or the space group
  is determined starting from the definition of the unit cell.
  Only the first method is implemented, the second one should be
  a wrapper for the symmetry finder library. To implement the
  second case I have to add additional entries in the object
  and I have also to pass an object describing the (optional) geometry builder.

INPUTS

  natom=number of atom
  ntypat=number of type of atoms
  nsym=number of symmetry operations
  rprimd(3,3)=dimensional lattive vector (real space)
  typat(natom)=type of each atom
  xred(3,natom)=reduced coordinates of each atom
  symrel(3,3,nsym) [optional]=symmetry operations in real space
  space_group=Space group (0 if not available)
  tnons(3,nsym) [optional]=fractional Translations
  symafm(nsym) [optional]=  ferromagnetic symmetries
  remove_inv [optional]= if .TRUE. the inversion is removed from the set of symmetries
  timrev ==2 => take advantage of time-reversal symmetry
         ==1 ==> do not use time-reversal symmetry

OUTPUT

  Cryst<crystal_t>= the object completely initialized.

TODO

  Add additional entries in the class:
  1) Info on space and point group (generators?).
  2) alchemy
  3) masses and nuclear (pseudo&AE) charge
  4) forces stresses, velocities.
  5) constraints for the relaxation
  6) Likely I will need also info on the electric field and berryopt

PARENTS

      dfpt_looppert,eig2tot,gwls_hamiltonian,m_crystal_io,m_ddb
      m_effective_potential,m_effective_potential_file,m_tdep_abitypes,mover
      optic,outscfcv,respfn,vtorho

CHILDREN

      mati3inv,sg_multable

SOURCE

250 subroutine crystal_init(amu,Cryst,space_group,natom,npsp,ntypat,nsym,rprimd,typat,xred,&
251 & zion,znucl,timrev,use_antiferro,remove_inv,title,&
252 & symrel,tnons,symafm) ! Optional
253 
254 
255 !This section has been created automatically by the script Abilint (TD).
256 !Do not modify the following lines by hand.
257 #undef ABI_FUNC
258 #define ABI_FUNC 'crystal_init'
259 !End of the abilint section
260 
261  implicit none
262 
263 !Arguments ------------------------------------
264 !scalars
265  integer,intent(in) :: natom,ntypat,nsym,timrev,space_group,npsp
266  type(crystal_t),intent(inout) :: Cryst
267  logical,intent(in) :: remove_inv,use_antiferro
268 !arrays
269  integer,intent(in) :: typat(natom)
270  integer,optional,intent(in) :: symrel(3,3,nsym),symafm(nsym)
271  real(dp),intent(in) :: amu(ntypat),xred(3,natom),rprimd(3,3),zion(ntypat),znucl(npsp)
272  real(dp),optional,intent(in) :: tnons(3,nsym)
273  character(len=*),intent(in) :: title(ntypat)
274 
275 !Local variables-------------------------------
276 !scalars
277  integer :: iat,indx,itypat,pinv,isym,nsym_noI
278  real(dp) :: tolsym8,ucvol
279  character(len=500) :: msg
280 !arrays
281  integer :: symrec(3,3)
282  real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3)
283  integer,pointer :: symrel_noI(:,:,:)
284  integer,allocatable :: indsym(:,:,:)
285  real(dp),pointer :: tnons_noI(:,:)
286 ! *************************************************************************
287 
288  !@crystal_t
289  Cryst%natom  = natom
290  Cryst%ntypat = ntypat
291  Cryst%npsp   = npsp
292 
293  Cryst%space_group = space_group
294 
295  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
296 
297  Cryst%ucvol  = ucvol
298  Cryst%rprimd = rprimd
299  Cryst%rmet   = rmet
300  Cryst%gmet   = gmet
301  Cryst%gprimd = gprimd
302 
303  Cryst%angdeg(1)=ACOS(Cryst%rmet(2,3)/SQRT(Cryst%rmet(2,2)*Cryst%rmet(3,3)))/two_pi*360.0d0
304  Cryst%angdeg(2)=ACOS(Cryst%rmet(1,3)/SQRT(Cryst%rmet(1,1)*Cryst%rmet(3,3)))/two_pi*360.0d0
305  Cryst%angdeg(3)=ACOS(Cryst%rmet(1,2)/SQRT(Cryst%rmet(1,1)*Cryst%rmet(2,2)))/two_pi*360.0d0
306 
307  ABI_MALLOC(Cryst%typat,(natom))
308  ABI_MALLOC(Cryst%xred,(3,natom))
309  ABI_MALLOC(Cryst%xcart,(3,natom))
310  ABI_MALLOC(Cryst%zion,(ntypat))
311  ABI_MALLOC(Cryst%znucl,(npsp))
312  ABI_MALLOC(Cryst%amu, (ntypat))
313 
314  Cryst%amu   = amu
315  Cryst%typat = typat
316  Cryst%xred  = xred
317  Cryst%zion  = zion
318  Cryst%znucl = znucl
319 
320  call xred2xcart(natom,rprimd,Cryst%xcart,Cryst%xred)
321 
322  ABI_MALLOC(Cryst%title,(ntypat))
323  Cryst%title = title
324  !
325  ! === Generate index table of atoms, in order for them to be used type after type ===
326  ABI_MALLOC(Cryst%atindx,(natom))
327  ABI_MALLOC(Cryst%atindx1,(natom))
328  ABI_MALLOC(Cryst%nattyp,(ntypat))
329 
330  indx=1
331  do itypat=1,ntypat
332    Cryst%nattyp(itypat)=0
333    do iat=1,natom
334      if (Cryst%typat(iat)==itypat) then
335        Cryst%atindx (iat )=indx
336        Cryst%atindx1(indx)=iat
337        indx=indx+1
338        Cryst%nattyp(itypat)=Cryst%nattyp(itypat)+1
339      end if
340    end do
341  end do
342 
343  Cryst%timrev = timrev
344 
345  if (PRESENT(symrel).and.PRESENT(tnons).and.PRESENT(symafm)) then
346    if (.not.remove_inv) then
347      ! * Just a copy
348      Cryst%nsym= nsym
349      ABI_MALLOC(Cryst%symrel,(3,3,nsym))
350      ABI_MALLOC(Cryst%symrec,(3,3,nsym))
351      ABI_MALLOC(Cryst%tnons,(3,nsym))
352      ABI_MALLOC(Cryst%symafm,(nsym))
353      Cryst%symrel=symrel
354      Cryst%tnons=tnons
355      Cryst%symafm=symafm
356      Cryst%use_antiferro = use_antiferro
357      do isym=1,nsym
358        call mati3inv(symrel(:,:,isym),symrec)
359        Cryst%symrec(:,:,isym)=symrec
360      end do
361    else
362      ! * Remove inversion, just to be compatible with old GW implementation
363      ! TODO should be removed!
364      call remove_inversion(nsym,symrel,tnons,nsym_noI,symrel_noI,tnons_noI,pinv)
365      Cryst%nsym=nsym_noI
366      ABI_MALLOC(Cryst%symrel,(3,3,nsym_noI))
367      ABI_MALLOC(Cryst%symrec,(3,3,nsym_noI))
368      ABI_MALLOC(Cryst%tnons,(3,nsym_noI))
369      ABI_MALLOC(Cryst%symafm,(nsym_noI))
370      Cryst%symrel=symrel_noI
371      Cryst%tnons=tnons_noI
372      if (ANY(symafm==-1)) then
373        msg = 'Solve the problem with inversion before adding ferromagnetic symmetries '
374        MSG_BUG(msg)
375      end if
376      Cryst%symafm=1
377      Cryst%use_antiferro=use_antiferro
378      do isym=1,nsym_noI
379        call mati3inv(symrel_noI(:,:,isym),symrec)
380        Cryst%symrec(:,:,isym)=symrec
381      end do
382      ABI_FREE(symrel_noI)
383      ABI_FREE(tnons_noI)
384    end if
385 
386  else
387    ! * Find symmetries symrec,symrel,tnons,symafm
388    ! TODO This should be a wrapper around the abinit library whose usage is not so straightforward
389    MSG_BUG('NotImplememented: symrel, symrec and tnons should be specied')
390  end if
391 
392  ! === Obtain a list of rotated atoms ===
393  ! $ R^{-1} (xred(:,iat)-\tau) = xred(:,iat_sym) + R_0 $
394  ! * indsym(4,  isym,iat) gives iat_sym in the original unit cell.
395  ! * indsym(1:3,isym,iat) gives the lattice vector $R_0$.
396  !
397  ABI_MALLOC(indsym,(4,Cryst%nsym,natom)); indsym(:,:,:)=zero
398  tolsym8=tol8
399  call symatm(indsym,natom,Cryst%nsym,Cryst%symrec,Cryst%tnons,tolsym8,Cryst%typat,Cryst%xred)
400 
401  ABI_MALLOC(Cryst%indsym,(4,Cryst%nsym,natom))
402  Cryst%indsym=indsym
403  ABI_FREE(indsym)
404 
405  ! === Rotation in spinor space ===
406  ABI_MALLOC(Cryst%spinrot,(4,Cryst%nsym))
407  do isym=1,Cryst%nsym
408    call getspinrot(Cryst%rprimd,Cryst%spinrot(:,isym),Cryst%symrel(:,:,isym))
409  end do
410 
411 end subroutine crystal_init

m_crystal/crystal_point_group [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

  crystal_point_group

FUNCTION

  Return the symmetries of the point group of the crystal.

INPUTS

  [include_timrev]=If True, time-reversal symmetry is included in the point group unless
    the system has spatial inversion. Default: False

OUTPUT

  ptg_nsym=Number of symmetries in the point group
  ptg_symrel(3,3,ptg_nsym)=Rotations in real space
  ptg_symrec(3,3,ptg_nsym)=Rotations in reciprocal space
  has_inversion=True if spatial inversion is present in the point group.

PARENTS

      m_skw

CHILDREN

      mati3inv,sg_multable

SOURCE

 923 subroutine crystal_point_group(cryst, ptg_nsym, ptg_symrel, ptg_symrec, has_inversion, include_timrev)
 924 
 925 
 926 !This section has been created automatically by the script Abilint (TD).
 927 !Do not modify the following lines by hand.
 928 #undef ABI_FUNC
 929 #define ABI_FUNC 'crystal_point_group'
 930 !End of the abilint section
 931 
 932  implicit none
 933 
 934 !Arguments ------------------------------------
 935 !scalars
 936  type(crystal_t),intent(in) :: cryst
 937  integer,intent(out) :: ptg_nsym
 938  logical,optional,intent(in) :: include_timrev
 939  logical,intent(out) :: has_inversion
 940 !arrays
 941  integer,allocatable,intent(out) :: ptg_symrel(:,:,:),ptg_symrec(:,:,:)
 942 
 943 !Local variables-------------------------------
 944 !scalars
 945  logical :: debug
 946  integer :: isym,search,tmp_nsym,ierr
 947  logical :: found,my_include_timrev
 948 !arrays
 949  integer :: work_symrel(3,3,cryst%nsym)
 950  integer,allocatable :: symafm(:)
 951  real(dp),allocatable :: tnons(:,:)
 952 
 953 ! *************************************************************************
 954 
 955  my_include_timrev = .False.; if (present(include_timrev)) my_include_timrev = include_timrev
 956 
 957  tmp_nsym = 1; work_symrel(:,:,1) = cryst%symrel(:,:,1)
 958  do isym=2,cryst%nsym
 959    if (cryst%symafm(isym) == -1) cycle
 960    do search=1,tmp_nsym
 961      found = all(work_symrel(:,:,search) == cryst%symrel(:,:,isym))
 962      if (found) exit
 963    end do
 964    if (.not. found) then
 965      tmp_nsym = tmp_nsym + 1
 966      work_symrel(:,:,tmp_nsym) = cryst%symrel(:,:,isym)
 967    end if
 968  end do
 969 
 970  has_inversion = .False.
 971  do isym=1,tmp_nsym
 972    if (all(work_symrel(:,:,isym) == inversion_3d) ) then
 973      has_inversion = .True.; exit
 974    end if
 975  end do
 976 
 977  ! Now we know the symmetries of the point group.
 978  ptg_nsym = tmp_nsym; if (.not. has_inversion .and. my_include_timrev) ptg_nsym = 2 * tmp_nsym
 979  ABI_MALLOC(ptg_symrel, (3,3,ptg_nsym))
 980  ABI_MALLOC(ptg_symrec, (3,3,ptg_nsym))
 981 
 982  ptg_symrel(:,:,1:tmp_nsym) = work_symrel(:,:,1:tmp_nsym)
 983  do isym=1,tmp_nsym
 984    call mati3inv(ptg_symrel(:,:,isym), ptg_symrec(:,:,isym))
 985  end do
 986 
 987  if (.not. has_inversion .and. my_include_timrev) then
 988    ptg_symrel(:,:,tmp_nsym+1:) = -work_symrel(:,:,1:tmp_nsym)
 989    do isym=tmp_nsym+1,ptg_nsym
 990      call mati3inv(ptg_symrel(:,:,isym), ptg_symrec(:,:,isym))
 991    end do
 992  end if
 993 
 994  debug = .False.
 995  if (debug) then
 996    ABI_CALLOC(tnons, (3, ptg_nsym))
 997    ABI_MALLOC(symafm, (ptg_nsym))
 998    symafm = 1
 999    call sg_multable(ptg_nsym, symafm, ptg_symrel, tnons, tol12, ierr)
1000    ABI_CHECK(ierr == 0, "point group is not a group! See messages above")
1001    ABI_FREE(tnons)
1002    ABI_FREE(symafm)
1003  end if
1004 
1005 end subroutine crystal_point_group

m_crystal/crystal_print [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

  crystal_print

FUNCTION

  Print the content of crystal_t data type

INPUTS

  Cryst<crystal_t>=The structure.
  [unit]=Unit number for output. Defaults to std_out
  [prtvol]=Verbosity level. If prtvol== -1, only lattice parameters are printed. Defaults to 0
  [mode_paral]=Either "COLL" or "PERS"
  [header]=String to be printed as header for additional info.

OUTPUT

  Only printing

PARENTS

      eph,gwls_hamiltonian,m_dvdb,m_gruneisen,setup_bse,setup_screening
      setup_sigma,wfk_analyze

CHILDREN

      mati3inv,sg_multable

SOURCE

542 subroutine crystal_print(Cryst,header,unit,mode_paral,prtvol)
543 
544 
545 !This section has been created automatically by the script Abilint (TD).
546 !Do not modify the following lines by hand.
547 #undef ABI_FUNC
548 #define ABI_FUNC 'crystal_print'
549 !End of the abilint section
550 
551  implicit none
552 
553 !Arguments ------------------------------------
554 !scalars
555  integer,optional,intent(in) :: unit,prtvol
556  character(len=*),optional,intent(in) :: mode_paral
557  character(len=*),optional,intent(in) :: header
558  type(crystal_t),intent(in) :: Cryst
559 
560 !Local variables-------------------------------
561  integer :: my_unt,my_prtvol,nu,iatom
562  character(len=4) :: my_mode
563  character(len=500) :: msg
564 ! *********************************************************************
565 
566  my_unt   =std_out; if (PRESENT(unit      )) my_unt   =unit
567  my_prtvol=0      ; if (PRESENT(prtvol    )) my_prtvol=prtvol
568  my_mode  ='COLL' ; if (PRESENT(mode_paral)) my_mode  =mode_paral
569 
570  msg=' ==== Info on the Cryst% object ==== '
571  if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
572  call wrtout(my_unt,msg,my_mode)
573 
574  write(msg,'(a)')' Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):'
575  call wrtout(my_unt,msg,my_mode)
576  do nu=1,3
577    write(msg,'(1x,a,i1,a,3f11.7,2x,a,i1,a,3f11.7)')&
578     'R(',nu,')=',Cryst%rprimd(:,nu)+tol10,&
579     'G(',nu,')=',Cryst%gprimd(:,nu)+tol10 !tol10 is used to be consistent with metric.F90
580    call wrtout(my_unt,msg,my_mode)
581  end do
582 
583  write(msg,'(a,1p,e15.7,a)')' Unit cell volume ucvol=',Cryst%ucvol+tol10,' bohr^3'
584  call wrtout(my_unt,msg,my_mode)
585 
586  write(msg,'(a,3es16.8,a)')' Angles (23,13,12)=',Cryst%angdeg(1:3),' degrees'
587  call wrtout(my_unt,msg,my_mode)
588 
589  if (Cryst%timrev==1) then
590    msg = ' Time-reversal symmetry is not present '
591  else if (Cryst%timrev==2) then
592    msg = ' Time-reversal symmetry is present '
593  else
594    MSG_BUG('Wrong value for timrev')
595  end if
596  call wrtout(my_unt,msg,my_mode)
597  if (my_prtvol == -1) return
598 
599  call print_symmetries(Cryst%nsym,Cryst%symrel,Cryst%tnons,Cryst%symafm,unit=my_unt,mode_paral=my_mode)
600 
601  if (Cryst%use_antiferro) call wrtout(my_unt,' System has magnetic symmetries ',my_mode)
602 
603  call wrtout(my_unt,"Reduced atomic positions [iatom, xred, symbol]:",my_mode)
604  do iatom=1,cryst%natom
605    write(msg,"(i5,a,2x,3f11.7,2x,a)")iatom,")",cryst%xred(:,iatom),symbol_type(cryst,cryst%typat(iatom))
606    call wrtout(my_unt,msg,my_mode)
607  end do
608 
609 end subroutine crystal_print

m_crystal/crystal_t [ Types ]

[ Top ] [ m_crystal ] [ Types ]

NAME

 crystal_t

FUNCTION

 Structure defining the unit cell (geometry, atomic positions and symmetry operations in real and reciprocal space)

SOURCE

 57  type,public :: crystal_t
 58 
 59 !scalars
 60   !integer :: point_group                    ! Point group
 61   !integer :: bravais,crystsys               ! Bravais lattice, Crystal system
 62   !integer :: nptsym                         ! No of point symmetries of the Bravais lattice
 63   !integer :: bravais(11)                    ! bravais(1)=iholohedry, bravais(2)=center
 64                                              ! bravais(3:11)=coordinates of rprim in the axes of the conventional
 65                                              ! bravais lattice (*2 if center/=0)
 66 
 67   !integer :: vacuum(3)
 68   !integer,pointer ptsymrel(:,:,:)
 69   !ptsymrel(3,3,nptsym)
 70   ! nptsym point-symmetry operations of the Bravais lattice in real space in terms of primitive translations.
 71 
 72   integer :: natom
 73   ! Number of atoms
 74 
 75   integer :: nsym
 76   ! Number of symmetry operations
 77 
 78   integer :: ntypat
 79   ! Number of type of atoms
 80 
 81   !$integer :: ntypalch,ntyppure
 82 
 83   integer :: npsp
 84   ! No. of pseudopotentials
 85 
 86   integer :: space_group
 87   ! Space group
 88 
 89   integer :: timrev
 90   ! TODO BE CAREFUL here, as the convention used in abinit is different.
 91   ! 1 => do not use time-reversal symmetry.
 92   ! 2 => take advantage of time-reversal symmetry.
 93 
 94   real(dp) :: ucvol
 95   ! Real space unit cell volume.
 96 
 97   logical :: use_antiferro
 98   ! .TRUE. if AFM symmetries are present and used.
 99 
100 !arrays
101   real(dp) :: angdeg(3)
102   ! Angles among rprim (degree).
103 
104   real(dp) :: gmet(3,3)
105   ! Reciprocal space metric ($\textrm{bohr}^{-2}$).
106 
107   real(dp) :: gprimd(3,3)
108   ! Dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
109 
110   real(dp) :: rmet(3,3)
111   ! Metric in real space.
112 
113   real(dp) :: rprimd(3,3)
114   ! Direct lattice vectors, Bohr units.
115 
116   integer,allocatable :: indsym(:,:,:)
117   ! indsym(4,nsym,natom)
118   ! indirect indexing array for atoms, see symatm.F90.
119 
120   integer,allocatable :: symafm(:)
121   ! symafm(nsym)
122   ! (Anti)Ferromagnetic symmetries. +1/-1
123 
124   integer,allocatable :: symrec(:,:,:)
125   ! symrec(3,3,nsym)
126   ! Symmetry operation in reciprocal space (reduced coordinates)
127 
128   integer,allocatable :: symrel(:,:,:)
129   ! symrel(3,3,nsym)
130   ! Symmetry operations in direct space (reduced coordinates).
131 
132   integer,allocatable :: atindx(:)
133   integer,allocatable :: atindx1(:)
134   ! atindx(natom), atindx1(natom)
135   ! Index tables for atoms useful to treat atoms type after type.
136 
137   integer,allocatable :: typat(:)
138   integer,allocatable :: nattyp(:)
139   ! typat(natom), nattyp(ntypat)
140   ! Type of each natom and number of atoms of each type.
141 
142   real(dp),allocatable :: tnons(:,:)
143   ! tnons(3,nsym)
144   ! Fractional translations (reduced coordinates)
145 
146   real(dp),allocatable :: xcart(:,:)
147   ! xcart(3,natom)
148   ! Cartesian coordinates.
149 
150   real(dp),allocatable :: xred(:,:)
151   ! xred(3,natom)
152   ! Reduced coordinates.
153 
154   real(dp),allocatable :: spinrot(:,:)
155   ! spinrot(4,nsym)
156   ! spinor rotation matrices.
157 
158   ! Useful quantities that might be added in the future
159   real(dp),allocatable :: amu(:)
160   !  amu(ntypat)
161   !  mass of the atoms (atomic mass unit)
162 
163   real(dp),allocatable :: zion(:)
164   ! zion(ntypat)
165   ! Charge of the pseudo-ion
166   ! (No of valence electrons needed to screen exactly the pseudopotential).
167 
168   !real(dp),allocatable :: znucltypat(:)
169    ! znucltypat(ntypat)
170    ! The atomic number of each type of atom (might be alchemy wrt psps)
171 
172   real(dp),allocatable :: znucl(:)
173   ! znucl(npsp)
174   ! Nuclear charge for each type of pseudopotential
175 
176   character(len=132),allocatable :: title(:)
177    ! title(ntypat)
178    ! The content of first line read from the psp file
179 
180  end type crystal_t
181 
182  public :: crystal_init            ! Main Creation method.
183  public :: crystal_free            ! Free memory.
184  public :: crystal_print           ! Print dimensions and basic info stored in the object
185  public :: idx_spatial_inversion   ! Return the index of the spatial inversion, 0 if not present.
186  public :: isymmorphic             ! True if space group is symmorphic.
187  public :: isalchemical            ! True if we are using alchemical pseudopotentials.
188  public :: adata_type              ! Return atomic data from the itypat index.
189  public :: symbol_type             ! Return the atomic symbol from the itypat index.
190  public :: symbols_crystal         ! Return an array with the atomic symbol:["Sr","Ru","O1","O2","O3"]
191  public :: crystal_point_group     ! Return the symmetries of the point group of the crystal.
192  public :: prt_cif                 ! Print CIF file.
193  public :: prtposcar               ! output VASP style POSCAR and FORCES files.

m_crystal/idx_spatial_inversion [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

  idx_spatial_inversion

FUNCTION

  Return the index of the spatial inversion, 0 if not present

PARENTS

CHILDREN

SOURCE

709 pure function idx_spatial_inversion(Cryst) result(inv_idx)
710 
711 
712 !This section has been created automatically by the script Abilint (TD).
713 !Do not modify the following lines by hand.
714 #undef ABI_FUNC
715 #define ABI_FUNC 'idx_spatial_inversion'
716 !End of the abilint section
717 
718  implicit none
719 
720 !Arguments ------------------------------------
721 !scalars
722  integer :: inv_idx
723  type(crystal_t),intent(in) :: Cryst
724 
725 !Local variables-------------------------------
726 !scalars
727  integer :: isym
728 
729 ! *************************************************************************
730 
731  inv_idx=0
732  do isym=1,cryst%nsym
733    if (all(cryst%symrel(:,:,isym) == inversion_3d)) then
734     inv_idx=isym; return
735    end if
736  end do
737 
738 end function idx_spatial_inversion

m_crystal/isalchemical [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

  isalchemical

FUNCTION

  Returns .TRUE. if we are using alchemical pseudopotentials

PARENTS

CHILDREN

SOURCE

794 pure function isalchemical(Cryst) result(ans)
795 
796 
797 !This section has been created automatically by the script Abilint (TD).
798 !Do not modify the following lines by hand.
799 #undef ABI_FUNC
800 #define ABI_FUNC 'isalchemical'
801 !End of the abilint section
802 
803  implicit none
804 
805 !Arguments ------------------------------------
806 !scalars
807  logical :: ans
808  type(crystal_t),intent(in) :: Cryst
809 
810 ! *************************************************************************
811 
812  ans = (Cryst%npsp /= Cryst%ntypat)
813 
814 end function isalchemical

m_crystal/isymmorphic [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

  isymmorphic

FUNCTION

  Returns .TRUE. if space group is symmorphic, i.e. all fractional translations are zero.

PARENTS

CHILDREN

SOURCE

756 pure function isymmorphic(Cryst) result(ans)
757 
758 
759 !This section has been created automatically by the script Abilint (TD).
760 !Do not modify the following lines by hand.
761 #undef ABI_FUNC
762 #define ABI_FUNC 'isymmorphic'
763 !End of the abilint section
764 
765  implicit none
766 
767 !Arguments ------------------------------------
768 !scalars
769  logical :: ans
770  type(crystal_t),intent(in) :: Cryst
771 
772 ! *************************************************************************
773 
774  ans = ALL (ABS(Cryst%tnons)<tol6)
775 
776 end function isymmorphic

m_crystal/prt_cif [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

 prt_cif

FUNCTION

   print out CIF format file

INPUTS

OUTPUT

NOTES

PARENTS

      outscfcv

CHILDREN

SOURCE

1030 subroutine prt_cif(brvltt, ciffname, natom, nsym, ntypat, rprimd, &
1031 &   spgaxor, spgroup, spgorig, symrel, tnon, typat, xred, znucl)
1032 
1033 
1034 !This section has been created automatically by the script Abilint (TD).
1035 !Do not modify the following lines by hand.
1036 #undef ABI_FUNC
1037 #define ABI_FUNC 'prt_cif'
1038 !End of the abilint section
1039 
1040  implicit none
1041 
1042 !Arguments ------------------------------------
1043 !scalars
1044  integer,intent(in) :: natom, ntypat, nsym
1045  integer, intent(in) :: brvltt, spgaxor, spgroup, spgorig
1046 !arrays
1047  integer, intent(in) :: typat(natom)
1048  integer, intent(in) :: symrel(3,3,nsym)
1049  character(len=*), intent(in) :: ciffname
1050  real(dp), intent(in) :: tnon(3,nsym)
1051  real(dp), intent(in) :: rprimd(3,3)
1052  real(dp), intent(in) :: xred(3,natom)
1053  real(dp), intent(in) :: znucl(ntypat)
1054 
1055 !Local variables -------------------------------
1056 !scalars
1057  integer :: unitcif
1058  integer :: iatom, isym
1059  integer :: sporder
1060  integer :: itypat, nat_this_type
1061  real(dp) :: ucvol
1062  type(atomdata_t) :: atom
1063 
1064 !arrays
1065  character(len=80) :: tmpstring
1066 
1067  character(len=1) :: brvsb
1068  character(len=15) :: intsb,ptintsb,ptschsb,schsb
1069  character(len=35) :: intsbl
1070 
1071  character(len=10) :: str_nat_type
1072  character(len=100) :: chemformula
1073  character(len=500) :: msg
1074 
1075  real(dp) :: angle(3)
1076  real(dp) :: gprimd(3,3)
1077  real(dp) :: rmet(3,3), gmet(3,3)
1078 
1079 !*************************************************************************
1080 
1081 !open file in append mode xlf and other compilers refuse append mode
1082  if (open_file(ciffname,msg,newunit=unitcif) /=0) then
1083    MSG_WARNING(msg)
1084    return
1085  end if
1086 
1087 !print title for dataset
1088  write (unitcif,'(a)') 'data_set'
1089 
1090 !print cell parameters a,b,c, angles, volume
1091  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1092  angle(1)=acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))/two_pi*360.0_dp
1093  angle(2)=acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))/two_pi*360.0_dp
1094  angle(3)=acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))/two_pi*360.0_dp
1095 
1096  write (unitcif,'(a,E20.10)') '_cell_length_a                     ', sqrt(rmet(1,1))*Bohr_Ang
1097  write (unitcif,'(a,E20.10)') '_cell_length_b                     ', sqrt(rmet(2,2))*Bohr_Ang
1098  write (unitcif,'(a,E20.10)') '_cell_length_c                     ', sqrt(rmet(3,3))*Bohr_Ang
1099  write (unitcif,'(a,E20.10)') '_cell_angle_alpha                  ', angle(1)
1100  write (unitcif,'(a,E20.10)') '_cell_angle_beta                   ', angle(2)
1101  write (unitcif,'(a,E20.10)') '_cell_angle_gamma                  ', angle(3)
1102  write (unitcif,'(a,E20.10)') '_cell_volume                       ', ucvol*(Bohr_Ang)**3
1103 
1104 !print reduced positions
1105  write (unitcif,'(a)') 'loop_'
1106  write (unitcif,'(a,E20.10)') '  _atom_site_label                   '
1107  write (unitcif,'(a,E20.10)') '  _atom_site_fract_x                 '
1108  write (unitcif,'(a,E20.10)') '  _atom_site_fract_y                 '
1109  write (unitcif,'(a,E20.10)') '  _atom_site_fract_z                 '
1110  do iatom = 1, natom
1111    call atomdata_from_znucl(atom,znucl(typat(iatom)))
1112    write (unitcif,'(2a,3E20.10)') '  ', atom%symbol, xred(:,iatom)
1113  end do
1114 
1115 !
1116 !other specs in CIF dictionary which may be useful:
1117 !GEOM_BOND GEOM_ANGLE GEOM_TORSION
1118 !
1119 
1120 !print chemical composition in simplest form
1121  chemformula = "'"
1122  do itypat = 1, ntypat
1123    nat_this_type = 0
1124    do iatom = 1, natom
1125      if (typat(iatom) == itypat) nat_this_type = nat_this_type+1
1126    end do
1127    call atomdata_from_znucl(atom,znucl(itypat))
1128    call int2char10(nat_this_type, str_nat_type)
1129    chemformula = trim(chemformula) // atom%symbol // trim(str_nat_type) // "  "
1130  end do
1131  chemformula = trim(chemformula) // "'"
1132  write (unitcif,'(2a)') '_chemical_formula_analytical              ', chemformula
1133 
1134 !FIXME: check that brvltt is correctly used here - is it equal to bravais(1) in the invars routines?
1135  if     (brvltt==1)then
1136    write (unitcif,'(a)') '_symmetry_cell_setting             triclinic'
1137  else if(brvltt==2)then
1138    write (unitcif,'(a)') '_symmetry_cell_setting             monoclinic'
1139  else if(brvltt==3)then
1140    write (unitcif,'(a)') '_symmetry_cell_setting             orthorhombic'
1141  else if(brvltt==4)then
1142    write (unitcif,'(a)') '_symmetry_cell_setting             tetragonal'
1143  else if(brvltt==5)then
1144    write (unitcif,'(a)') '_symmetry_cell_setting             rhombohedral'
1145  else if(brvltt==6)then
1146    write (unitcif,'(a)') '_symmetry_cell_setting             hexagonal'
1147  else if(brvltt==7)then
1148    write (unitcif,'(a)') '_symmetry_cell_setting             cubic'
1149  end if
1150 
1151  call spgdata(brvsb,intsb,intsbl,ptintsb,ptschsb,schsb,spgaxor,spgroup,sporder,spgorig)
1152 
1153 !print symmetry operations
1154  write (unitcif,'(a,I6)') "_symmetry_Int_Tables_number          ", spgroup
1155  write (unitcif,'(5a)') "_symmetry_space_group_name_H-M        '", brvsb, " ", trim(intsb), "'"
1156  write (unitcif,'(a)') ''
1157  write (unitcif,'(a)') 'loop_'
1158  write (unitcif,'(a)') '  _symmetry_equiv_pos_as_xyz           '
1159  do isym = 1, nsym
1160    call  symrel2string(symrel(:,:,isym), tnon(:,isym), tmpstring)
1161    write (unitcif,'(2a)') '  ', trim(tmpstring)
1162  end do
1163 
1164  close (unitcif)
1165 
1166 end subroutine prt_cif

m_crystal/prtposcar [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

  prtposcar

FUNCTION

  output VASP style POSCAR and FORCES files for use with frozen phonon codes, like
  PHON from Dario Alfe' or frophon
  IMPORTANT: the order of atoms is fixed such that typat is re-grouped. First typat=1
  then typat=2, etc...

INPUTS

  fcart = forces on atoms in cartesian coordinates
  natom = number of atoms
  ntypat = number of types of atoms
  rprimd = lattice vectors for the primitive cell
  typat = type for each of the natom atoms
  ucvol = unit cell volume
  xred = reduced positions of the atoms
  znucl = nuclear charge of each atomic type

 OUTPUTS
   Only files written

PARENTS

      afterscfloop

CHILDREN

      atomdata_from_znucl

SOURCE

1276 subroutine prtposcar(fcart, fnameradix, natom, ntypat, rprimd, typat, ucvol, xred, znucl)
1277 
1278 
1279 !This section has been created automatically by the script Abilint (TD).
1280 !Do not modify the following lines by hand.
1281 #undef ABI_FUNC
1282 #define ABI_FUNC 'prtposcar'
1283 !End of the abilint section
1284 
1285  implicit none
1286 
1287 !Arguments ------------------------------------
1288 !scalars
1289  integer, intent(in) :: natom, ntypat
1290  real(dp), intent(in) :: ucvol
1291 !arrays
1292  integer, intent(in) :: typat(natom)
1293  real(dp), intent(in) :: fcart(3,natom)
1294  real(dp), intent(in) :: rprimd(3,3)
1295  real(dp), intent(in) :: xred(3,natom)
1296  real(dp), intent(in) :: znucl(ntypat)
1297  character(len=fnlen), intent(in) :: fnameradix
1298 
1299 !Local variables-------------------------------
1300 !scalars
1301  integer :: iatom, itypat, iout
1302  type(atomdata_t) :: atom
1303 ! arrays
1304  integer :: natoms_this_type(ntypat)
1305  character(len=2) :: symbol
1306  character(len=7) :: natoms_this_type_str
1307  character(len=100) :: chem_formula, natoms_all_types
1308  character(len=500) :: msg
1309  character(len=fnlen) :: fname
1310 
1311 !************************************************************************
1312 
1313 !output POSCAR file for positions, atom types etc
1314  fname = trim(fnameradix)//"_POSCAR"
1315  if (open_file(fname,msg,newunit=iout) /= 0) then
1316    MSG_ERROR(msg)
1317  end if
1318 
1319  natoms_this_type = 0
1320  do itypat=1, ntypat
1321    do iatom=1,natom
1322      if (typat(iatom) == itypat) natoms_this_type(itypat) = natoms_this_type(itypat)+1
1323    end do
1324  end do
1325 
1326  chem_formula = ""
1327  do itypat=1, ntypat
1328    call atomdata_from_znucl(atom,znucl(itypat))
1329    symbol = atom%symbol
1330    if (natoms_this_type(itypat) < 10) then
1331      write (natoms_this_type_str, '(I1)') natoms_this_type(itypat)
1332    else if (natoms_this_type(itypat) < 100) then
1333      write (natoms_this_type_str, '(I2)') natoms_this_type(itypat)
1334    else if (natoms_this_type(itypat) < 1000) then
1335      write (natoms_this_type_str, '(I3)') natoms_this_type(itypat)
1336    end if
1337    chem_formula = trim(chem_formula) // symbol // trim(natoms_this_type_str)
1338  end do
1339  write (iout,'(2a)') "ABINIT generated POSCAR file. Title string - should be chemical formula... ",&
1340 & trim(chem_formula)
1341 
1342  write (iout,'(E24.14)') -ucvol*Bohr_Ang*Bohr_Ang*Bohr_Ang
1343  write (iout,'(3E24.14,1x)') Bohr_Ang*rprimd(:,1) ! (angstr? bohr?)
1344  write (iout,'(3E24.14,1x)') Bohr_Ang*rprimd(:,2)
1345  write (iout,'(3E24.14,1x)') Bohr_Ang*rprimd(:,3)
1346  natoms_all_types = "   "
1347  do itypat=1, ntypat
1348    write (natoms_this_type_str, '(I7)') natoms_this_type(itypat)
1349    natoms_all_types = trim(natoms_all_types) // "   " // trim(natoms_this_type_str)
1350  end do
1351  write (iout,'(a)') trim(natoms_all_types)
1352  write (iout,'(a)') "Direct"
1353  do itypat=1, ntypat
1354    do iatom=1,natom
1355      if (typat(iatom) /= itypat) cycle
1356      write (iout,'(3(E24.14,1x))') xred(:,iatom)
1357    end do
1358  end do
1359  close (iout)
1360 
1361 
1362 !output FORCES file for forces in same order as positions above
1363  fname = trim(fnameradix)//"_FORCES"
1364  if (open_file(fname,msg,newunit=iout) /= 0 ) then
1365    MSG_ERROR(msg)
1366  end if
1367 !ndisplacements
1368 !iatom_displaced displacement_red_coord(3)
1369 !forces_cart_ev_Angstr(3)
1370 !...
1371 !<repeat for other displaced atoms>
1372  write (iout,'(I7)') 1
1373  write (iout,'(a)') '1 0 0 0        ! TO BE FILLED IN '
1374  do itypat=1, ntypat
1375    do iatom=1,natom
1376      if (typat(iatom) /= itypat) cycle
1377      write (iout,'(3(E24.14,1x))') Ha_eV/Bohr_Ang*fcart(:,iatom)
1378    end do
1379  end do
1380  close (iout)
1381 
1382 end subroutine prtposcar

m_crystal/symbol_type [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

  symbol_type

FUNCTION

  Return the atomic symbol from the itypat index

PARENTS

SOURCE

867 function symbol_type(crystal,itypat) result(symbol)
868 
869 
870 !This section has been created automatically by the script Abilint (TD).
871 !Do not modify the following lines by hand.
872 #undef ABI_FUNC
873 #define ABI_FUNC 'symbol_type'
874 !End of the abilint section
875 
876  implicit none
877 
878 !Arguments ------------------------------------
879 !scalars
880  integer,intent(in) :: itypat
881  character(len=2) :: symbol
882  type(crystal_t),intent(in) :: crystal
883 
884 !Local variables-------------------------------
885 !scalars
886  type(atomdata_t) :: atom
887 
888 ! *************************************************************************
889 
890  atom = adata_type(crystal, itypat)
891  symbol = atom%symbol
892 
893 end function symbol_type

m_crystal/symbols_crystal [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

 symbols_crystal

FUNCTION

 Return a array with the symbol of each atoms
 with indexation
 ["Sr","Ru","O1","O2","O3"] for example

INPUTS

 natom = number of atoms
 ntypat = number of typat
 npsp =  number of pseudopotentials
 znucl = Nuclear charge for each type of pseudopotential

OUTPUT

 symbols = array with the symbol of each atoms

PARENTS

      m_effective_potential_file,m_fit_polynomial_coeff,m_polynomial_coeff

CHILDREN

      mati3inv,sg_multable

SOURCE

641 subroutine symbols_crystal(natom,ntypat,npsp,symbols,typat,znucl)
642 
643 
644 !This section has been created automatically by the script Abilint (TD).
645 !Do not modify the following lines by hand.
646 #undef ABI_FUNC
647 #define ABI_FUNC 'symbols_crystal'
648 !End of the abilint section
649 
650  implicit none
651 
652 !Arguments ------------------------------------
653 !scalars
654  integer,intent(in) :: natom,ntypat,npsp
655 !arrays
656  real(dp),intent(in):: znucl(npsp)
657  integer,intent(in) :: typat(natom)
658  character(len=5),intent(out) :: symbols(natom)
659  character(len=3) :: powerchar
660 !Local variables-------------------------------
661 !scalar
662  integer :: ia,ii,itypat,jj
663 !arrays
664 ! *************************************************************************
665 
666 !  Fill the symbols array
667    do ia=1,natom
668      symbols(ia) = adjustl(znucl2symbol(znucl(typat(ia))))
669    end do
670    itypat = zero
671    do itypat =1,ntypat
672      ii = zero
673      do ia=1,natom
674        if(typat(ia)==itypat) then
675          ii = ii + 1
676        end if
677      end do
678      if(ii>1)then
679        jj=1
680        do ia=1,natom
681          if(typat(ia)==itypat) then
682            write(powerchar,'(I0)') jj
683            symbols(ia) = trim(symbols(ia))//trim(powerchar)
684            jj=jj+1
685          end if
686        end do
687      end if
688    end do
689 
690 end subroutine symbols_crystal

m_crystal/symrel2string [ Functions ]

[ Top ] [ m_crystal ] [ Functions ]

NAME

 symrel2string

FUNCTION

INPUTS

OUTPUT

NOTES

PARENTS

      prt_cif

CHILDREN

SOURCE

1188 subroutine symrel2string(symrel1, tnon, string)
1189 
1190 
1191 !This section has been created automatically by the script Abilint (TD).
1192 !Do not modify the following lines by hand.
1193 #undef ABI_FUNC
1194 #define ABI_FUNC 'symrel2string'
1195 !End of the abilint section
1196 
1197  implicit none
1198 
1199 !Arguments ------------------------------------
1200 !scalars
1201 !arrays
1202  integer, intent(in) :: symrel1(3,3)
1203  real(dp), intent(in) :: tnon(3)
1204  character(len=80), intent(out) :: string
1205 
1206 !Local variables -------------------------
1207 !scalars
1208  integer :: i1,i2
1209  character(len=1) :: xyz(3)
1210 
1211 ! *********************************************************************
1212 
1213  xyz(1) = 'x'
1214  xyz(2) = 'y'
1215  xyz(3) = 'z'
1216 
1217  string = ''
1218  do i1=1,3
1219    if (abs(tnon(i1)) > tol10) then
1220 !    find fraction 1/n for tnon, otherwise do not know what to print
1221      if (abs(one-two*tnon(i1)) < tol10) string = trim(string)//'1/2'
1222      if (abs(one+two*tnon(i1)) < tol10) string = trim(string)//'-1/2'
1223 
1224      if (abs(one-three*tnon(i1)) < tol10) string = trim(string)//'1/3'
1225      if (abs(one+three*tnon(i1)) < tol10) string = trim(string)//'-1/3'
1226      if (abs(two-three*tnon(i1)) < tol10) string = trim(string)//'2/3'
1227      if (abs(two+three*tnon(i1)) < tol10) string = trim(string)//'-2/3'
1228 
1229      if (abs(one-six*tnon(i1)) < tol10) string = trim(string)//'1/6'
1230      if (abs(one+six*tnon(i1)) < tol10) string = trim(string)//'-1/6'
1231      if (abs(five-six*tnon(i1)) < tol10) string = trim(string)//'5/6'
1232      if (abs(five+six*tnon(i1)) < tol10) string = trim(string)//'-5/6'
1233    end if
1234    do i2=1,3
1235 !    FIXME: check if this is correct ordering for symrel(i1,i2) looks ok
1236      if (symrel1(i1,i2) == 1)  string = trim(string)//'+'//xyz(i2)
1237      if (symrel1(i1,i2) == -1) string = trim(string)//'-'//xyz(i2)
1238    end do
1239    if (i1 /= 3) string = trim(string)//','
1240  end do
1241 
1242 end subroutine symrel2string