TABLE OF CONTENTS


ABINIT/m_symfind [ Modules ]

[ Top ] [ Modules ]

NAME

  m_symfind

FUNCTION

  Symmetry finder high-level API.

COPYRIGHT

  Copyright (C) 2000-2017 ABINIT group (XG, RC)
  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_symfind
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32  use m_symlist
33 
34  use m_symtk,     only : chkprimit, matr3inv, symrelrot, symdet, symcharac, holocell, smallprim
35  use m_geometry,  only : xred2xcart
36  use m_spgdata,   only : getptgroupma, symptgroup, spgdata
37 
38  implicit none
39 
40  private

m_symfind/symanal [ Functions ]

[ Top ] [ m_symfind ] [ Functions ]

NAME

 symanal

FUNCTION

 Find the space group, Bravais lattice, including Shubnikov characteristics
 from the list of symmetries (including magnetic characteristics), and lattice parameters
 Warning: the recognition of the space group might not yet work for the
 Shubnikov group of type IV

INPUTS

 chkprim= if 1 then stop if the cell is not primitive
 msym=default maximal number of symmetries
 nsym=actual number of symmetries
 rprimd(3,3)=dimensional primitive translations for real space (bohr)
 symafm(1:msym)=(anti)ferromagnetic part of symmetry operations
 symrel(3,3,1:msym)=symmetry operations in real space in terms
  of primitive translations
 tnons(3,1:msym)=nonsymmorphic translations for symmetry operations
 tolsym=tolerance for the symmetry operations

OUTPUT

 bravais(11)=characteristics of Bravais lattice (see symlatt.F90)
 genafm(3)=magnetic translation generator (in case of Shubnikov group type IV)
 ptgroupma = magnetic point group number
 spgroup=symmetry space group

PARENTS

      ingeo,m_ab7_symmetry,m_tdep_sym,m_use_ga

CHILDREN

      chkprimit,getptgroupma,symbrav,symlatt,symptgroup,symspgr,wrtout

SOURCE

496 subroutine symanal(bravais,chkprim,genafm,msym,nsym,ptgroupma,rprimd,spgroup,symafm,symrel,tnons,tolsym)
497 
498 
499 !This section has been created automatically by the script Abilint (TD).
500 !Do not modify the following lines by hand.
501 #undef ABI_FUNC
502 #define ABI_FUNC 'symanal'
503 !End of the abilint section
504 
505  implicit none
506 
507 !Arguments ------------------------------------
508 !scalars
509  integer,intent(in) :: chkprim,msym,nsym
510  integer,intent(out) :: ptgroupma,spgroup
511  real(dp),intent(in) :: tolsym
512 !arrays
513  integer,intent(out) :: bravais(11)
514  integer,intent(in) :: symafm(msym),symrel(3,3,msym)
515  real(dp),intent(in) :: rprimd(3,3)
516  real(dp),intent(inout) :: tnons(3,msym)
517  real(dp),intent(out) :: genafm(3)
518 
519 !Local variables-------------------------------
520 !scalars
521  integer :: iholohedry_nomagn,isym,isym_nomagn,multi
522  integer :: nptsym,nsym_nomagn,shubnikov
523  character(len=5) :: ptgroup,ptgroupha
524  character(len=500) :: message
525 !arrays
526  integer :: identity(3,3)
527  integer,allocatable :: ptsymrel(:,:,:),symrel_nomagn(:,:,:)
528  real(dp),allocatable :: tnons_nomagn(:,:)
529 
530 ! *************************************************************************
531 
532 !This routine finds the Bravais characteristics, without actually
533 !looking at the symmetry operations.
534  ABI_ALLOCATE(ptsymrel,(3,3,msym))
535  call symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym)
536  ABI_DEALLOCATE(ptsymrel)
537 
538 !Check whether the cell is primitive or not.
539  call chkprimit(chkprim,multi,nsym,symafm,symrel)
540 
541  spgroup=0 ; ptgroupma=0 ; genafm(:)=zero
542 
543  if(multi>1)then !  Modify bravais if the cell is not primitive ; no determination of the space group
544    bravais(1)=-bravais(1)
545  else
546 
547 !  The cell is primitive, so that the space group can be
548 !  determined. Need to distinguish Fedorov and Shubnikov groups.
549 !  Do not distinguish Shubnikov types I and II.
550 !  Also identify genafm, in case of Shubnikov type IV
551    identity(:,:)=reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/))
552    shubnikov=1
553    do isym=1,nsym
554      if(symafm(isym)==-1)then
555        shubnikov=3
556        if(sum(abs(symrel(:,:,isym)-identity(:,:)))==0)then
557          shubnikov=4
558          genafm(:)=tnons(:,isym)
559 !        DEBUG
560 !        write(std_out,*)' isym=',isym
561 !        write(std_out,*)' symrel(:,:,isym)',symrel(:,:,isym)
562 !        write(std_out,*)' tnons(:,isym)',tnons(:,isym)
563 !        write(std_out,*)' symafm(isym)',symafm(isym)
564 !        ENDDEBUG
565          exit
566        end if
567      end if
568    end do
569 
570    if(shubnikov/=1)then
571      if(shubnikov==3)write(message, '(a)' )' Shubnikov space group type III'
572      if(shubnikov==4)write(message, '(a)' )' Shubnikov space group type IV'
573      call wrtout(std_out,message,'COLL')
574    end if
575 
576    if(shubnikov==1 .or. shubnikov==3)then
577 !    Find the correct Bravais characteristics and point group
578 !    Should also be used for Shubnikov groups of type IV ...
579      call symbrav(bravais,msym,nsym,ptgroup,rprimd,symrel,tolsym)
580 !    Find the space group
581      call symspgr(bravais,nsym,spgroup,symrel,tnons,tolsym)
582    end if
583 
584    if(shubnikov/=1)then
585 
586 !    Determine nonmagnetic symmetry operations
587      nsym_nomagn=nsym/2
588      ABI_ALLOCATE(symrel_nomagn,(3,3,nsym_nomagn))
589      ABI_ALLOCATE(tnons_nomagn,(3,nsym_nomagn))
590      isym_nomagn=0
591      do isym=1,nsym
592        if(symafm(isym)==1)then
593          isym_nomagn=isym_nomagn+1
594          symrel_nomagn(:,:,isym_nomagn)=symrel(:,:,isym)
595          tnons_nomagn(:,isym_nomagn)=tnons(:,isym)
596        end if
597      end do
598 
599      if(shubnikov==3)then
600 
601 !      DEBUG
602 !      write(std_out,*)' symanal : will enter symbrav with halved symmetry set'
603 !      write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
604 !      do isym=1,nsym_nomagn
605 !      write(std_out,'(i3,2x,9i3,3es12.2,i3)')isym,symrel_nomagn(:,:,isym),tnons_nomagn(:,isym)
606 !      end do
607 !      ENDDEBUG
608 
609 !      Find the point group of the halved symmetry set
610        call symptgroup(iholohedry_nomagn,nsym_nomagn,ptgroupha,symrel_nomagn)
611 
612 !      Deduce the magnetic point group (ptgroupma) from ptgroup and ptgroupha
613        call getptgroupma(ptgroup,ptgroupha,ptgroupma)
614 
615      else if(shubnikov==4)then
616 
617 !      Find the Fedorov space group of the halved symmetry set
618        call symspgr(bravais,nsym_nomagn,spgroup,symrel_nomagn,tnons_nomagn,tolsym)
619 
620 !      The magnetic translation generator genafm has already been determined
621 !      write(std_out,*)' genafm =',genafm, ' spgroup=',spgroup
622 
623      end if
624 
625      ABI_DEALLOCATE(symrel_nomagn)
626      ABI_DEALLOCATE(tnons_nomagn)
627    end if ! Shubnikov groups
628 
629  end if
630 
631 end subroutine symanal

m_symfind/symbrav [ Functions ]

[ Top ] [ m_symfind ] [ Functions ]

NAME

 symbrav

FUNCTION

 From the list of symmetry operations, and the lattice vectors,
 determine the Bravais information (including the holohedry, the centering,
 the coordinate of the primitive vectors in the conventional vectors),
 as well as the point group.

INPUTS

 msym=dimension of symrel
 nsym=actual number of symmetries
 rprimd(3,3)=dimensional primitive translations for real space (bohr)
 symrel(3,3,msym)=symmetry operations in real space in terms
                  of primitive translations
 tolsym=tolerance for the symmetries

OUTPUT

 bravais(11): bravais(1)=iholohedry
              bravais(2)=center
              bravais(3:11)=coordinates of rprimd in the axes
              of the conventional bravais lattice (*2 if center/=0)
 ptgroup=symmetry point group
 [axis(3)]=Invariant axis in the conventional vector coordinates
   Set to (/0,0,0/) if the lattice belongs to the same holohedry as the lattice+atoms (+electric field + ...).

PARENTS

      m_esymm,symanal

CHILDREN

      matr3inv,symlatt,symptgroup,symrelrot

SOURCE

669 subroutine symbrav(bravais,msym,nsym,ptgroup,rprimd,symrel,tolsym,axis)
670 
671 
672 !This section has been created automatically by the script Abilint (TD).
673 !Do not modify the following lines by hand.
674 #undef ABI_FUNC
675 #define ABI_FUNC 'symbrav'
676 !End of the abilint section
677 
678  implicit none
679 
680 !Arguments ------------------------------------
681 !scalars
682  integer,intent(in) :: msym,nsym
683  real(dp),intent(in) :: tolsym
684  character(len=5),intent(out) :: ptgroup
685 !arrays
686  integer,intent(in) :: symrel(3,3,msym)
687  integer,optional,intent(out) :: axis(3)
688  integer,intent(out) :: bravais(11)
689  real(dp),intent(in) :: rprimd(3,3)
690 
691 !Local variables-------------------------------
692 !scalars
693  integer :: iaxis,ii,bravais1now,ideform,iholohedry,invariant,isym
694  integer :: jaxis,next_stage,nptsym,problem
695  real(dp) :: norm,scprod
696  character(len=500) :: message
697 !arrays
698  integer :: identity(3,3),axis_trial(3),hexa_axes(3,7),ortho_axes(3,13)
699  integer,allocatable :: ptsymrel(:,:,:),symrelconv(:,:,:)
700  real(dp) :: axes(3,3),axis_cart(3),axis_red(3)
701  real(dp) :: rprimdconv(3,3),rprimdtry(3,3),rprimdnow(3,3)
702  real(dp) :: rprimdconv_invt(3,3)
703 
704 !**************************************************************************
705 
706  identity(:,:)=0
707  identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
708 
709  ortho_axes(:,:)=0
710  ortho_axes(1,1)=1
711  ortho_axes(2,2)=1
712  ortho_axes(3,3)=1
713  ortho_axes(:,4)=(/0,1,1/)
714  ortho_axes(:,5)=(/1,0,1/)
715  ortho_axes(:,6)=(/1,1,0/)
716  ortho_axes(:,7)=(/0,1,-1/)
717  ortho_axes(:,8)=(/-1,0,1/)
718  ortho_axes(:,9)=(/1,-1,0/)
719  ortho_axes(:,10)=(/1,1,1/)
720  ortho_axes(:,11)=(/-1,1,1/)
721  ortho_axes(:,12)=(/1,-1,1/)
722  ortho_axes(:,13)=(/1,1,-1/)
723 
724  hexa_axes(:,:)=0
725  hexa_axes(1,1)=1
726  hexa_axes(2,2)=1
727  hexa_axes(3,3)=1
728  hexa_axes(:,4)=(/1,-1,0/)
729  hexa_axes(:,5)=(/2,1,0/)
730  hexa_axes(:,6)=(/1,1,0/)
731  hexa_axes(:,7)=(/1,2,0/)
732 
733 !Determine the point group from the list of symmetry operations.
734 !Also determine the holohedry, up to one undeterminacy : hR versus hP
735  call symptgroup(iholohedry,nsym,ptgroup,symrel)
736 
737 !Loop over trial deformations
738 !This is needed in case the Bravais lattice determination from the lattice vectors
739 !has a higher holohedry than the real one, in which the symmetry
740 !operations for the atoms (or electric field, etc) are taken into account
741  iaxis=0
742  invariant=0
743  next_stage=0
744  rprimdnow(:,:)=rprimd(:,:)
745  rprimdtry(:,:)=rprimd(:,:)
746  ABI_ALLOCATE(symrelconv,(3,3,nsym))
747 
748 !At most will have to try 65 deformations (13 axes, five stages)
749  do ideform=1,65
750 
751    ABI_ALLOCATE(ptsymrel,(3,3,msym))
752    call symlatt(bravais,msym,nptsym,ptsymrel,rprimdtry,tolsym)
753    ABI_DEALLOCATE(ptsymrel)
754 
755 !  Examine the agreement with bravais(1)
756 !  Warning : might change Bravais lattice hR to hP, if hexagonal axes
757    problem=0
758    select case (bravais(1))
759    case(7)
760      if(iholohedry<6)problem=1
761      if(iholohedry==6)problem=2
762    case(6)
763      if(iholohedry<4)problem=1
764      if(iholohedry==7 .or. iholohedry==4)problem=2
765 !      Here, change hR into hP
766      if(iholohedry==5)iholohedry=6
767    case(5)
768      if(iholohedry<4)problem=1
769      if(iholohedry==7 .or. iholohedry==6 .or. iholohedry==4)problem=2
770    case(4)
771      if(iholohedry<4)problem=1
772      if(iholohedry>4)problem=2
773    case(3)
774      if(iholohedry<3)problem=1
775      if(iholohedry>3)problem=2
776    case(2)
777      if(iholohedry<2)problem=1
778      if(iholohedry>2)problem=2
779    case(1)
780      if(iholohedry>1)problem=2
781    end select
782 
783 !  This is the usual situation, in which the lattice belong to the same holohedry
784 !  as the lattice+atoms (+electric field + ...)
785    if(problem==0)exit
786 
787    if(problem==2)then
788      if(iaxis==0)then
789        write(message, '(3a,i3,3a,i3,7a)' )&
790 &       'The Bravais lattice determined only from the primitive',ch10,&
791 &       'vectors (rprim or angdeg), bravais(1)=',bravais(1),', is not compatible',ch10,&
792 &       'with the real one, iholohedry=',iholohedry,', obtained by taking into',ch10,&
793 &       'account the symmetry operations. This might be due to an insufficient',ch10,&
794 &       'number of digits in the specification of rprim (at least 10),',ch10,&
795 &       'or to an erroneous rprim or angdeg. If this is not the case, then ...'
796        MSG_BUG(message)
797      end if
798      if(iaxis==1)then
799        write(message, '(3a,3i3,2a,i3,2a,i3)' )&
800 &       'Could not succeed to determine the bravais lattice',ch10,&
801 &       'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
802 &       'bravais(1)=',bravais(1),ch10,&
803 &       'iholohedry=',iholohedry
804        MSG_BUG(message)
805      end if
806    end if
807 
808    if(problem==1)then  ! One is left with the problem=1 case, basically iholohedry is lower than bravais(1)
809      if(iaxis==0)then
810        write(message, '(a,a,a,i3,a,a,a,i3,a,a,a)' )&
811 &       'The Bravais lattice determined only from the primitive',ch10,&
812 &       'vectors, bravais(1)=',bravais(1),', is more symmetric',ch10,&
813 &       'than the real one, iholohedry=',iholohedry,', obtained by taking into',ch10,&
814 &       'account the atomic positions. Start deforming the primitive vector set.'
815        MSG_COMMENT(message)
816        next_stage=1
817      else if(iaxis/=0)then
818        if(bravais(1)<bravais1now)then
819          write(message, '(3a,i3,3a,i3,2a)' )&
820 &         'The Bravais lattice determined from modified primitive',ch10,&
821 &         'vectors, bravais(1)=',bravais(1),', has a lower symmetry than before,',ch10,&
822 &         'but is still more symmetric than the real one, iholohedry=',iholohedry,ch10,&
823 &         'obtained by taking into account the atomic positions.'
824          MSG_COMMENT(message)
825          next_stage=1
826        else if(iaxis==1)then
827          write(message, '(3a,3i3,2a,i3,2a,i3)' )&
828 &         'Could not succeed to determine the bravais lattice',ch10,&
829 &         'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
830 &         'bravais(1)=',bravais(1),ch10,&
831 &         'iholohedry=',iholohedry
832          MSG_BUG(message)
833        end if
834      end if
835    end if ! problem==1
836 
837    if(next_stage==1)then
838      bravais1now=bravais(1)
839      rprimdnow(:,:)=rprimdtry(:,:)
840 !    Generate the symmetry operations in the conventional vector coordinates
841      rprimdconv(:,1)=bravais(3:5)
842      rprimdconv(:,2)=bravais(6:8)
843      rprimdconv(:,3)=bravais(9:11)
844      axes(:,:)=zero
845      axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one
846      symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym)
847      call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym)
848      if(bravais(1)/=6)then
849        iaxis=14
850      else
851        iaxis=8
852      end if
853      next_stage=0
854    end if
855 
856    iaxis=iaxis-1
857    do jaxis=iaxis,1,-1
858      if(bravais(1)/=6)then
859        axis_trial(:)=ortho_axes(:,jaxis)
860      else
861        axis_trial(:)=hexa_axes(:,jaxis)
862      end if
863 !    DEBUG
864 !    write(std_out,*)' symbrav : try jaxis=',jaxis
865 !    write(std_out,*)' axis_trial=',axis_trial
866 !    ENDDEBUG
867      invariant=1
868 !    Examine whether all symmetry operations leave the axis invariant (might be reversed, though)
869      do isym=1,nsym
870        if(sum(abs(matmul(symrelconv(:,:,isym),axis_trial)+(-axis_trial(:))))/=0 .and. &
871 &       sum(abs(matmul(symrelconv(:,:,isym),axis_trial)+axis_trial(:)))/=0 )invariant=0
872      end do
873      if(invariant==1)then
874        iaxis=jaxis
875 !      write(message, '(2a,i3)' )ch10,' symbrav : found invariant axis, jaxis=',iaxis
876 !      call wrtout(std_out,message,'COLL')
877        exit
878      end if
879    end do
880 
881    if(invariant==0)then
882 !    Not a single axis was invariant with respect to all operations ?!
883 !    do isym=1,nsym; write(std_out, '(a,10i4)' )' isym,symrelconv=',isym,symrelconv(:,:,isym); enddo
884      write(message, '(3a,3i3,2a,i3,2a,i3)' )&
885 &     'Could not succeed to determine the bravais lattice (not a single invariant)',ch10,&
886 &     'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
887 &     'bravais(1)=',bravais(1),ch10,&
888 &     'iholohedry=',iholohedry
889      MSG_BUG(message)
890    end if
891 
892    call matr3inv(rprimdconv,rprimdconv_invt)
893    axis_red(:)=axis_trial(1)*rprimdconv_invt(1,:)+ &
894 &   axis_trial(2)*rprimdconv_invt(2,:)+ &
895 &   axis_trial(3)*rprimdconv_invt(3,:)
896    axis_cart(:)=axis_red(1)*rprimdnow(:,1)+ &
897 &   axis_red(2)*rprimdnow(:,2)+ &
898 &   axis_red(3)*rprimdnow(:,3)
899    norm=sum(axis_cart(:)**2)
900 !  Expand by a uniform, quite arbitrary, dilatation, along the invariant axis
901 !  Note : make these dilatation different, according to ideform
902 !  XG 20151221  : Still, the interplay between the size of the deformation and the tolsym is not easy to address.
903 !  Indeed the deformation must be sufficiently large to be perceived by symlatt as a real breaking of the
904 !  symmetry of the lattice. In order to deal with all the small values od tolsym, it has been set at a minimum of tol3,
905 !  but it must also be larger than tolsym. Moreover, for some axis choice, the deformation is not aligned with the axis, decreasing
906 !  the effective deformation length. An additional factor of three is thus included, actually increased to six just to be sure...
907    do ii=1,3
908      scprod=axis_cart(1)*rprimdnow(1,ii)+axis_cart(2)*rprimdnow(2,ii)+axis_cart(3)*rprimdnow(3,ii)
909      rprimdtry(:,ii)=rprimdnow(:,ii)+ideform*(max(tol3,six*tolsym)-tol6)*scprod/norm*axis_cart(:)
910    end do
911 
912  end do ! ideform
913 
914  if(bravais(1)/=iholohedry)then
915    write(message, '(3a,3i3,2a,i3,2a,i3)' )&
916 &   'Despite efforts, Could not succeed to determine the bravais lattice :',ch10,&
917 &   'bravais(1)=',bravais(1),ch10,&
918 &   'iholohedry=',iholohedry
919    MSG_BUG(message)
920  end if
921 
922  ABI_DEALLOCATE(symrelconv)
923 
924  if (PRESENT(axis)) then  ! Return symmetry axis.
925    axis=(/0,0,0/)
926    if (iaxis/=0) then
927      if(bravais(1)/=6)then
928        axis=ortho_axes(:,iaxis)
929      else
930        axis=hexa_axes(:,iaxis)
931      end if
932    end if
933  end if
934 
935 end subroutine symbrav

m_symfind/symfind [ Functions ]

[ Top ] [ m_symfind ] [ Functions ]

NAME

 symfind

FUNCTION

 Symmetry finder.
 From the symmetries of the Bravais lattice (ptsymrel),
 select those that leave invariant the system, and generate
 the corresponding tnons vectors.
 The algorithm is explained in T.G. Worlton and J.L. Warren, Comp. Phys. Comm. 3, 88 (1972) [[cite:Worton1972]]

INPUTS

 berryopt    =  4/14, 6/16, 7/17: electric or displacement field
 efield=cartesian coordinates of the electric field
 gprimd(3,3)=dimensional primitive translations for reciprocal space
 msym=default maximal number of symmetries
 natom=number of atoms in cell.
 noncoll=1 if non-collinear magnetism is activated
         else 0
 nptsym=number of point symmetries of the Bravais lattice
 nucdipmom(3,natom) (optional) array of nuclear dipole moments
 nzchempot=if non-zero, means that a z-spatially varying chemical potential is added
 ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
   of the Bravais lattice in real space in terms
   of primitive translations.
 spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
 tolsym=tolerance for the symmetries
 typat(natom)=integer identifying type of atom.
 use_inversion=1 if inversion can be included in set of symmetries
 xred(3,natom)=reduced coordinates of atoms in terms of real space
   primitive translations

OUTPUT

 nsym=actual number of symmetries
 symafm(1:msym)=(anti)ferromagnetic part of nsym symmetry operations
 symrel(3,3,1:msym)= nsym symmetry operations in real space in terms
  of primitive translations
 tnons(3,1:msym)=nonsymmorphic translations for each symmetry (would
  be 0 0 0 each for a symmorphic space group)

PARENTS

      ingeo,inqpt,m_ab7_symmetry,m_effective_potential_file,m_tdep_sym
      m_use_ga,thmeig

CHILDREN

      wrtout

SOURCE

103  subroutine symfind(berryopt,efield,gprimd,jellslab,msym,natom,noncoll,nptsym,nsym,&
104 &  nzchempot,prtvol, ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
105 &  nucdipmom)
106 
107 
108 !This section has been created automatically by the script Abilint (TD).
109 !Do not modify the following lines by hand.
110 #undef ABI_FUNC
111 #define ABI_FUNC 'symfind'
112 !End of the abilint section
113 
114  implicit none
115 
116 !Arguments ------------------------------------
117 !scalars
118  integer,intent(in) :: berryopt,jellslab,msym,natom,noncoll,nptsym,nzchempot,use_inversion
119  integer,intent(in) :: prtvol
120  integer,intent(out) :: nsym
121  real(dp),intent(in) :: tolsym
122 !arrays
123  integer,intent(in) :: ptsymrel(3,3,msym),typat(natom)
124  integer,intent(inout) :: symafm(msym),symrel(3,3,msym) !vz_i
125  real(dp),intent(in) :: efield(3),gprimd(3,3),spinat(3,natom),xred(3,natom)
126  real(dp),optional :: nucdipmom(3,natom)
127  real(dp),intent(inout) :: tnons(3,msym) !vz_i
128 
129 !Local variables-------------------------------
130 ! TRUE if antiferro symmetries are used with non-collinear magnetism
131 !scalars
132  integer :: found3,foundcl,iatom,iatom0,iatom1,iatom2,iatom3,iclass,iclass0,ii
133  integer :: isym,jj,kk,natom0,nclass,ntrial,printed,trialafm,trialok
134  real(dp) :: spinatcl2,spinatcl20,det
135  logical,parameter :: afm_noncoll=.true.
136  logical :: test_sameabscollin,test_sameabsnoncoll,test_samespin
137  character(len=500) :: message
138 !arrays
139  integer,allocatable :: class(:,:),natomcl(:),typecl(:)
140  real(dp) :: diff(3),efieldrot(3),sxred0(3),symnucdipmom2(3)
141  real(dp) :: symspinat2(3),symxred2(3),trialnons(3)
142  real(dp),allocatable :: spinatcl(:,:),spinatred(:,:)
143 
144 !**************************************************************************
145 
146 !DEBUG
147  if (prtvol>1) message="remove me later"
148 ! write(std_out,*)' symfind : enter'
149 ! call flush(6)
150 ! write(std_out,*)' symfind : nzchempot= ',nzchempot
151 ! write(std_out,*)'   ptsymrel matrices are :'
152 ! do isym=1,nptsym
153 ! write(std_out,'(i4,4x,9i4)' )isym,ptsymrel(:,:,isym)
154 ! end do
155 ! write(std_out,*)' symfind : natom=',natom
156 ! do iatom=1,natom
157 ! write(std_out,*)'  atom number',iatom
158 ! write(std_out,*)'   typat   =',typat(iatom)
159 ! write(std_out,*)'   spinat  =',spinat(:,iatom)
160 ! write(std_out,*)'   xred    =',xred(:,iatom)
161 ! end do
162 ! call flush(6)
163 !ENDDEBUG
164 
165 !Find the number of classes of atoms (type and spinat must be identical,
166 !spinat might differ by a sign, if aligned with the z direction)
167 !natomcl(iclass) will contain the number of atoms in the class
168 !typecl(iclass) will contain the type of the atoms in the class
169 !spinatcl(1:3,iclass) will contain the spinat of the atoms in the class
170 !class(1:natomclass(iclass),iclass) will contain the index of the
171 !atoms belonging to the class
172  ABI_ALLOCATE(class,(natom+3,natom))
173  ABI_ALLOCATE(natomcl,(natom))
174  ABI_ALLOCATE(typecl,(natom))
175  ABI_ALLOCATE(spinatcl,(3,natom))
176 
177 !Initialise with the first atom
178  nclass=1
179  natomcl(1)=1
180  typecl(1)=typat(1)
181  spinatcl(:,1)=spinat(:,1)
182  class(1,1)=1
183  if(natom>1)then
184    do iatom=2,natom
185 !    DEBUG
186 !    write(std_out,*)' symfind : examine iatom=',iatom
187 !    ENDDEBUG
188      foundcl=0
189      do iclass=1,nclass
190 !      Compare the typat and spinat of atom iatom with existing ones.
191 !      Admit either identical spinat, or z-aligned spinat with same
192 !      absolute magnitude
193        if( typat(iatom)==typecl(iclass)) then
194          test_samespin=  &
195 &         abs(spinat(1,iatom)-spinatcl(1,iclass))<tolsym .and. &
196 &         abs(spinat(2,iatom)-spinatcl(2,iclass))<tolsym .and. &
197 &         abs(spinat(3,iatom)-spinatcl(3,iclass))<tolsym
198          test_sameabscollin= &
199 &         noncoll==0 .and.&
200 &         abs(spinat(1,iatom))<tolsym .and. abs(spinatcl(1,iclass))<tolsym .and.&
201 &         abs(spinat(2,iatom))<tolsym .and. abs(spinatcl(2,iclass))<tolsym .and.&
202 &         abs(abs(spinat(3,iatom))-abs(spinatcl(3,iclass)))<tolsym
203          test_sameabsnoncoll= &
204 &         noncoll==1 .and. afm_noncoll .and. &
205 &         abs(spinat(1,iatom)+spinatcl(1,iclass))<tolsym .and. &
206 &         abs(spinat(2,iatom)+spinatcl(2,iclass))<tolsym .and. &
207 &         abs(spinat(3,iatom)+spinatcl(3,iclass))<tolsym
208          if( test_samespin .or. test_sameabscollin .or. test_sameabsnoncoll) then
209 !          DEBUG
210 !          write(std_out,*)' symfind : find it belongs to class iclass=',iclass
211 !          write(std_out,*)' symfind : spinat(:,iatom)=',spinat(:,iatom)
212 !          write(std_out,*)' symfind : spinatcl(:,iclass)=',spinatcl(:,iclass)
213 !          write(std_out,*)' symfind : test_samespin,test_sameabscollin,test_sameabsnoncoll=',&
214 !          &      test_samespin,test_sameabscollin,test_sameabsnoncoll
215 !          ENDDEBUG
216            natomcl(iclass)=natomcl(iclass)+1
217            class(natomcl(iclass),iclass)=iatom
218            foundcl=1
219            exit
220          end if
221        end if
222      end do
223 !    If no class with these characteristics exist, create one
224      if(foundcl==0)then
225        nclass=nclass+1
226        natomcl(nclass)=1
227        typecl(nclass)=typat(iatom)
228        spinatcl(:,nclass)=spinat(:,iatom)
229        class(1,nclass)=iatom
230      end if
231    end do
232  end if
233 
234 !DEBUG
235 !write(std_out,*)' symfind : found ',nclass,' nclass of atoms'
236 !do iclass=1,nclass
237 !write(std_out,*)'  class number',iclass
238 !write(std_out,*)'   natomcl =',natomcl(iclass)
239 !write(std_out,*)'   typecl  =',typecl(iclass)
240 !write(std_out,*)'   spinatcl=',spinatcl(:,iclass)
241 !write(std_out,*)'   class   =',(class(iatom,iclass),iatom=1,natomcl(iclass))
242 !end do
243 !ENDDEBUG
244 
245 !Select the class with the least number of atoms, and non-zero spinat if any
246 !It is important to select a magnetic class of atom, if any, otherwise
247 !the determination of the initial (inclusive) set of symmetries takes only
248 !non-magnetic symmetries, and not both magnetic and non-magnetic ones, see later.
249  iclass0=1
250  natom0=natomcl(1)
251  spinatcl20=spinatcl(1,1)**2+spinatcl(2,1)**2+spinatcl(3,1)**2
252  if(nclass>1)then
253    do iclass=2,nclass
254      spinatcl2=spinatcl(1,iclass)**2+spinatcl(2,iclass)**2+spinatcl(3,iclass)**2
255      if( (natomcl(iclass)<natom0 .and. (spinatcl20<tolsym .or. spinatcl2>tolsym))  &
256 &     .or. (spinatcl20<tolsym .and. spinatcl2>tolsym)                         )then
257        iclass0=iclass
258        natom0=natomcl(iclass)
259        spinatcl20=spinatcl2
260      end if
261    end do
262  end if
263 
264  printed=0
265 
266 !DEBUG
267 !write(std_out,*)' symfind : has selected iclass0=',iclass0
268 !write(std_out,*)' #    iatom     xred             spinat '
269 !do iatom0=1,natomcl(iclass0)
270 !iatom=class(iatom0,iclass0)
271 !write(std_out,'(2i4,6f10.4)' )iatom0,iatom,xred(:,iatom),spinat(:,iatom)
272 !end do
273 !ENDDEBUG
274 
275 !If non-collinear spinat have to be used, transfer them in reduced coordinates
276  if (noncoll==1) then
277    ABI_ALLOCATE(spinatred,(3,natom))
278    do iatom=1,natom
279      do ii=1,3
280        spinatred(ii,iatom)=gprimd(1,ii)*spinat(1,iatom) &
281 &       +gprimd(2,ii)*spinat(2,iatom) &
282 &       +gprimd(3,ii)*spinat(3,iatom)
283      end do
284    end do
285  end if
286 
287 !Big loop over each symmetry operation of the Bravais lattice
288  nsym=0
289  do isym=1,nptsym
290 
291 !  ji: Check whether symmetry operation leaves efield invariant
292    if (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. &
293 &   berryopt==14 .or. berryopt==16 .or. berryopt==17) then
294      efieldrot(:) = ptsymrel(:,1,isym)*efield(1) +  &
295 &     ptsymrel(:,2,isym)*efield(2) +  &
296 &     ptsymrel(:,3,isym)*efield(3)
297      diff(:)=efield(:)-efieldrot(:)
298      if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 ) cycle
299    end if
300 
301    if (use_inversion==0) then
302      det=ptsymrel(1,1,isym)*ptsymrel(2,2,isym)*ptsymrel(3,3,isym)+&
303 &     ptsymrel(2,1,isym)*ptsymrel(3,2,isym)*ptsymrel(1,3,isym)+&
304 &     ptsymrel(1,2,isym)*ptsymrel(2,3,isym)*ptsymrel(3,1,isym) - &
305 &     (ptsymrel(3,1,isym)*ptsymrel(2,2,isym)*ptsymrel(1,3,isym)+&
306 &     ptsymrel(2,1,isym)*ptsymrel(1,2,isym)*ptsymrel(3,3,isym)+&
307 &     ptsymrel(3,2,isym)*ptsymrel(2,3,isym)*ptsymrel(1,1,isym))
308      if(det==-1) cycle
309    end if
310 
311 !  jellium slab and spatially varying chemical potential cases:
312 !  (actually, an inversion symmetry/mirror plane perpendicular to z symmetry operation might still be allowed... TO BE DONE !)
313    if (jellslab/=0 .or. nzchempot/=0) then
314 !    check whether symmetry operation produce a rotation only in the xy plane
315      if( ptsymrel(1,3,isym)/=0 .or. ptsymrel(2,3,isym)/=0 .or. &
316 &     ptsymrel(3,1,isym)/=0 .or. ptsymrel(3,2,isym)/=0 ) cycle
317 !    check whether symmetry operation does not change the z
318      if( ptsymrel(3,3,isym)/=1 ) cycle
319    end if
320 
321 !  Select a tentative set of associated translations
322 !  First compute the symmetric of the first atom in the smallest class,
323 !  using the point symmetry
324    iatom0=class(1,iclass0)
325    sxred0(:)=ptsymrel(:,1,isym)*xred(1,iatom0)+ &
326 &   ptsymrel(:,2,isym)*xred(2,iatom0)+ &
327 &   ptsymrel(:,3,isym)*xred(3,iatom0)
328 
329 !  From the set of possible images, deduce tentative translations,
330 !  and magnetic factor then test whether it send each atom on a symmetric one
331    ntrial=0
332    do ii=1,natom0
333      iatom1=class(ii,iclass0)
334 
335 !    The tentative translation is found
336      trialnons(:)=xred(:,iatom1)-sxred0(:)
337      trialafm=1
338      if(abs(spinat(3,iatom1)-spinat(3,iatom0))>tolsym)trialafm=-1
339 
340      if(sum(abs(spinat(:,iatom1)*trialafm-spinat(:,iatom0)))>tolsym)then
341        write(message,'(3a,3i5)')&
342 &       'Problem with matching the spin part within a class.',ch10,&
343 &       'isym,iatom0,iatom1=',isym,iatom0,iatom1
344        MSG_ERROR_CLASS(message, "TolSymError")
345      end if
346 !    jellium slab case: check whether symmetry operation has no translational
347 !    component along z
348      if( jellslab/=0 .and. abs(trialnons(3)) > tolsym ) cycle
349      trialok=1
350 
351 !    DEBUG
352 !    write(std_out,*)' isym,trialnons(:),trialafm =',isym,trialnons(:),trialafm
353 !    ENDDEBUG
354 
355 !    Loop over all classes, then all atoms in the class,
356 !    to find whether they have a symmetric
357      do iclass=1,nclass
358        do jj=1,natomcl(iclass)
359 
360          iatom2=class(jj,iclass)
361 !        Generate the tentative symmetric position of iatom2
362          symxred2(:)=ptsymrel(:,1,isym)*xred(1,iatom2)+ &
363 &         ptsymrel(:,2,isym)*xred(2,iatom2)+ &
364 &         ptsymrel(:,3,isym)*xred(3,iatom2)+ trialnons(:)
365 !        Generate the tentative symmetric spinat of iatom2
366          if (noncoll==0) then
367            symspinat2(:)=trialafm*spinat(:,iatom2)
368          else
369            symspinat2(:)=trialafm*(ptsymrel(:,1,isym)*spinatred(1,iatom2)+ &
370 &           ptsymrel(:,2,isym)*spinatred(2,iatom2)+ &
371 &           ptsymrel(:,3,isym)*spinatred(3,iatom2))
372          end if
373          if(present(nucdipmom)) then
374 !        Generate the tentative symmetric nuclear dipole moment of iatom2
375            symnucdipmom2(:)=(ptsymrel(:,1,isym)*nucdipmom(1,iatom2)+ &
376 &           ptsymrel(:,2,isym)*nucdipmom(2,iatom2)+ &
377 &           ptsymrel(:,3,isym)*nucdipmom(3,iatom2))
378          end if
379 
380 !        DEBUG
381 !        write(std_out,'(a,3f12.4,a,3f12.4)') ' Send atom at xred=',xred(:,iatom2),' to ',symxred2(:)
382 !        ENDDEBUG
383 
384 !        Check whether there exists an atom of the same class at the
385 !        same location, with the correct spinat and nuclear dipole moment
386          do kk=1,natomcl(iclass)
387 
388            found3=1
389            iatom3=class(kk,iclass)
390 !          Check the location
391            diff(:)=xred(:,iatom3)-symxred2(:)
392            diff(:)=diff(:)-nint(diff(:))
393            if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 )found3=0
394 !          Check the spinat
395            if (noncoll==0) then
396              diff(:)=spinat(:,iatom3)-symspinat2(:)
397            else
398              diff(:)=spinatred(:,iatom3)-symspinat2(:)
399            end if
400            if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 )found3=0
401            if(present(nucdipmom)) then
402 !            Check the nuclear dipole moment
403              diff(:) = nucdipmom(:,iatom3) - symnucdipmom2(:)
404              if(any(diff>tolsym))found3=0
405            end if
406 
407            if(found3==1)exit
408          end do ! End loop over iatom3
409 
410          if(found3==0)then
411            trialok=0
412            exit
413          end if
414        end do ! End loop over iatom2
415 
416        if(trialok==0)exit
417      end do ! End loop over all classes
418 
419      if(trialok==1)then
420        nsym=nsym+1
421        if(nsym>msym)then
422          write(message,'(3a,i0,4a)')&
423 &         'The number of symmetries (including non-symmorphic translations)',ch10,&
424 &         'is larger than maxnsym=',msym,ch10,&
425 &         'Action: increase maxnsym in the input, or take a cell that is primitive, ',ch10,&
426 &         'or at least smaller than the present one.'
427          MSG_ERROR(message)
428        end if
429        ntrial=ntrial+1
430        symrel(:,:,nsym)=ptsymrel(:,:,isym)
431        symafm(nsym)=trialafm
432        tnons(:,nsym)=trialnons(:)-nint(trialnons(:)-tolsym)
433      end if
434 
435    end do ! End the loop on tentative translations
436  end do ! End big loop over each symmetry operation of the Bravais lattice
437 
438  ABI_DEALLOCATE(class)
439  ABI_DEALLOCATE(natomcl)
440  ABI_DEALLOCATE(spinatcl)
441  ABI_DEALLOCATE(typecl)
442  if (noncoll==1)   then
443    ABI_DEALLOCATE(spinatred)
444  end if
445 
446 !DEBUG
447 ! write(message,'(a,I0,a)')' symfind : exit, nsym=',nsym,ch10
448 ! write(message,'(2a)') trim(message),'   symrel matrices, symafm and tnons are :'
449 ! call wrtout(std_out,message,'COLL')
450 ! do isym=1,nsym
451 !   write(message,'(i4,4x,3i4,2x,3i4,2x,3i4,4x,i4,4x,3f8.4)' ) isym,symrel(:,:,isym),&
452 !&   symafm(isym),tnons(:,isym)
453 !   call wrtout(std_out,message,'COLL')
454 ! end do
455 !stop
456 !ENDDEBUG
457 
458 end subroutine symfind

m_symfind/symlatt [ Functions ]

[ Top ] [ m_symfind ] [ Functions ]

NAME

 symlatt

FUNCTION

 From the unit cell vectors (rprimd) and the corresponding metric tensor,
 find the Bravais lattice and its symmetry operations (ptsymrel).
 1) Find the shortest possible primitive vectors for the lattice
 2) Determines the holohedral group of the lattice, and the
    axes to be used for the conventional cell
    (this is a delicate part, in which the centering of the
    reduced cell must be taken into account)
    The idea is to determine the basis vectors of the conventional
    cell from the reduced cell basis vectors.
 3) Generate the symmetry operations of the holohedral group

INPUTS

 msym=default maximal number of symmetries
 rprimd(3,3)=dimensional primitive translations for real space (bohr)
 tolsym=tolerance for the symmetries

OUTPUT

  bravais(11): bravais(1)=iholohedry
               bravais(2)=center
               bravais(3:11)=coordinates of rprim in the axes
               of the conventional bravais lattice (*2 if center/=0)
 nptsym=number of point symmetries of the Bravais lattice
 ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
 of the Bravais lattice in real space in terms of primitive translations.

NOTES

 WARNING: bravais(1) might be given a negative value in another
 routine, if the cell is non-primitive.
 The holohedral groups are numbered as follows
 (see international tables for crystallography (1983), p. 13)
 iholohedry=1   triclinic      1bar
 iholohedry=2   monoclinic     2/m
 iholohedry=3   orthorhombic   mmm
 iholohedry=4   tetragonal     4/mmm
 iholohedry=5   trigonal       3bar m
 iholohedry=6   hexagonal      6/mmm
 iholohedry=7   cubic          m3bar m
 Centering
 center=0        no centering
 center=-1       body-centered
 center=-3       face-centered
 center=1        A-face centered
 center=2        B-face centered
 center=3        C-face centered

PARENTS

      ingeo,inqpt,m_ab7_symmetry,m_effective_potential_file,m_tdep_sym
      m_use_ga,symanal,symbrav,thmeig

CHILDREN

      holocell,matr3inv,smallprim,symrelrot,wrtout

SOURCE

1429 subroutine symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym)
1430 
1431 
1432 !This section has been created automatically by the script Abilint (TD).
1433 !Do not modify the following lines by hand.
1434 #undef ABI_FUNC
1435 #define ABI_FUNC 'symlatt'
1436 !End of the abilint section
1437 
1438  implicit none
1439 
1440 !Arguments ------------------------------------
1441 !scalars
1442  integer,intent(in) :: msym
1443  integer,intent(out) :: nptsym
1444  real(dp),intent(in) :: tolsym
1445 !arrays
1446  integer,intent(out) :: bravais(11),ptsymrel(3,3,msym)
1447  real(dp),intent(in) :: rprimd(3,3)
1448 
1449 !Local variables-------------------------------
1450 !scalars
1451  integer,parameter :: mgen=4
1452  integer :: center,fact,found,foundc,ia,ib,icase,igen,iholohedry,ii,index,isym
1453  integer :: itrial,jj,jsym,ngen=0,orthogonal,sign12,sign13,sign23,sumsign
1454  real(dp) :: determinant,norm2a,norm2b,norm2c,norm2trial,reduceda,reducedb,sca
1455  real(dp) :: scalarprod,scb,trace,val
1456  character(len=500) :: message
1457 !arrays
1458  integer,parameter :: list_holo(7)=(/7,6,4,3,5,2,1/)
1459  integer :: ang90(3),equal(3),gen(3,3,mgen),gen2xy(3,3),gen2y(3,3),gen2z(3,3)
1460  integer :: gen3(3,3),gen6(3,3),icoord(3,3),identity(3,3),nvecta(3),nvectb(3)
1461  integer :: order(mgen)
1462  real(dp) :: axes(3,3),axesinvt(3,3),cell_base(3,3),coord(3,3),metmin(3,3)
1463  real(dp) :: minim(3,3),scprods(3,3),vecta(3),vectb(3),vectc(3),vin1(3),vin2(3),vext(3)
1464 
1465 !**************************************************************************
1466 
1467  identity(:,:)=0 ; identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
1468  nvecta(1)=2 ; nvectb(1)=3
1469  nvecta(2)=1 ; nvectb(2)=3
1470  nvecta(3)=1 ; nvectb(3)=2
1471 
1472 !--------------------------------------------------------------------------
1473 !Reduce the input vectors to a set of minimal vectors
1474  call smallprim(metmin,minim,rprimd)
1475 
1476 !DEBUG
1477 !write(std_out,*)' symlatt : minim(:,1)=',minim(:,1)
1478 !write(std_out,*)' symlatt : minim(:,2)=',minim(:,2)
1479 !write(std_out,*)' symlatt : minim(:,3)=',minim(:,3)
1480 !ENDDEBUG
1481 
1482 !--------------------------------------------------------------------------
1483 !Examine the angles and vector lengths
1484  ang90(:)=0
1485  if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1
1486  if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1
1487  if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1
1488  equal(:)=0
1489  if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1
1490  if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1
1491  if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1
1492 
1493 !DEBUG
1494 !write(std_out,*)' ang90=',ang90(:)
1495 !write(std_out,*)' equal=',equal(:)
1496 !ENDDEBUG
1497 
1498 !-----------------------------------------------------------------------
1499 !Identification of the centering
1500 
1501  foundc=0
1502 !Default values
1503  fact=1 ; center=0
1504  cell_base(:,:)=minim(:,:)
1505 
1506 !Examine each holohedral group
1507 !This search is ordered : should not be happy with tetragonal,
1508 !while there is FCC ...
1509  do index=1,6
1510 
1511 !  If the holohedry is already found, exit
1512    if(foundc==1)exit
1513 
1514 !  Initialize the target holohedry
1515    iholohedry=list_holo(index)
1516 
1517 !  DEBUG
1518 !  write(std_out,*)' symlatt : trial holohedry',iholohedry
1519 !  ENDDEBUG
1520 
1521    orthogonal=0
1522    if(iholohedry==7 .or. iholohedry==4 .or. iholohedry==3)orthogonal=1
1523 
1524 !  Now, will examine different working hypothesis.
1525 !  The set of these hypothesis is thought to cover all possible cases ...
1526 
1527 !  Working hypothesis : the basis is orthogonal
1528    if(ang90(1)+ang90(2)+ang90(3)==3 .and. orthogonal==1)then
1529      fact=1 ; center=0
1530      cell_base(:,:)=minim(:,:)
1531 !    Checks that the basis vectors are OK for the target holohedry
1532      call holocell(cell_base,0,foundc,iholohedry,tolsym)
1533    end if
1534 
1535 !  Select one trial direction
1536    do itrial=1,3
1537 
1538 !    If the holohedry is already found, exit
1539      if(foundc==1)exit
1540 
1541      ia=nvecta(itrial) ; ib=nvectb(itrial)
1542 
1543 !    This is in case of hexagonal holohedry
1544      if(foundc==0 .and. iholohedry==6 .and. ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then
1545        reduceda=metmin(ib,ia)/metmin(ia,ia)
1546        fact=1 ; center=0
1547        if(abs(reduceda+0.5d0)<tolsym)then
1548          cell_base(:,1)=minim(:,ia)
1549          cell_base(:,2)=minim(:,ib)
1550          cell_base(:,3)=minim(:,itrial)
1551 !        Checks that the basis vectors are OK for the target holohedry
1552          call holocell(cell_base,0,foundc,iholohedry,tolsym)
1553        else if(abs(reduceda-0.5d0)<tolsym)then
1554          cell_base(:,1)=minim(:,ia)
1555          cell_base(:,2)=minim(:,ib)-minim(:,ia)
1556          cell_base(:,3)=minim(:,itrial)
1557 !        Checks that the basis vectors are OK for the target holohedry
1558          call holocell(cell_base,0,foundc,iholohedry,tolsym)
1559        end if
1560      end if
1561 
1562 !    Working hypothesis : the conventional cell is orthogonal,
1563 !    and the two other vectors are axes of the conventional cell
1564      if(foundc==0 .and. orthogonal==1 .and. ang90(itrial)==1)then
1565 
1566 !      Compute the reduced coordinate of trial vector in the basis
1567 !      of the two other vectors
1568        reduceda=metmin(itrial,ia)/metmin(ia,ia)
1569        reducedb=metmin(itrial,ib)/metmin(ib,ib)
1570        cell_base(:,ia)=minim(:,ia)
1571        cell_base(:,ib)=minim(:,ib)
1572        if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. &
1573 &       ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym)       )then
1574          if(abs(abs(reduceda)-0.5d0)<tolsym)center=ib
1575          if(abs(abs(reducedb)-0.5d0)<tolsym)center=ia
1576          fact=2
1577          cell_base(:,itrial)= &
1578 &         (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0
1579          call holocell(cell_base,0,foundc,iholohedry,tolsym)
1580        else if( abs(abs(reduceda)-0.5d0)<tolsym .and.&
1581 &         abs(abs(reducedb)-0.5d0)<tolsym       ) then
1582          fact=2 ; center=-1
1583          cell_base(:,itrial)= &
1584 &         (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0
1585          call holocell(cell_base,0,foundc,iholohedry,tolsym)
1586        end if
1587      end if
1588 
1589 !    Working hypothesis : the conventional cell is orthogonal, and
1590 !    the trial vector is one of the future axes, and the face perpendicular to it is centered
1591      if(foundc==0 .and. iholohedry==3 .and. &
1592 &     ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then
1593        fact=2 ; center=itrial
1594        cell_base(:,ia)=minim(:,ia)+minim(:,ib)
1595        cell_base(:,ib)=minim(:,ia)-minim(:,ib)
1596        cell_base(:,itrial)=minim(:,itrial)
1597 !      Checks that the basis vectors are OK for the target holohedry
1598        call holocell(cell_base,0,foundc,iholohedry,tolsym)
1599      end if
1600 
1601 !    DEBUG
1602 !    write(std_out,*)' after test_b, foundc=',foundc
1603 !    ENDDEBUG
1604 
1605 !    Working hypothesis : the conventional cell is orthogonal, and
1606 !    the trial vector is one of the future axes
1607      if(foundc==0 .and. orthogonal==1)then
1608 !      Compute the projection of the two other vectors on the trial vector
1609        reduceda=metmin(itrial,ia)/metmin(itrial,itrial)
1610        reducedb=metmin(itrial,ib)/metmin(itrial,itrial)
1611 !      If both projections are half-integer, one might have found an axis
1612        if( abs(abs(reduceda)-0.5d0)<tolsym .and.&
1613 &       abs(abs(reducedb)-0.5d0)<tolsym       ) then
1614          vecta(:)=minim(:,ia)-reduceda*minim(:,itrial)
1615          vectb(:)=minim(:,ib)-reducedb*minim(:,itrial)
1616          norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
1617          norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
1618          scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
1619 !        Note the order of selection : body-centered is prefered
1620 !        over face centered, which is correct for the tetragonal case
1621          if(abs(norm2a-norm2b)<tolsym*half*(norm2a+norm2b))then
1622 !          The lattice is body centered
1623            fact=2 ; center=-1
1624            cell_base(:,ia)=vecta(:)+vectb(:)
1625            cell_base(:,ib)=vecta(:)-vectb(:)
1626            cell_base(:,itrial)=minim(:,itrial)
1627            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1628          else if(abs(scalarprod)<tolsym*half*(norm2a+norm2b))then
1629 !          The lattice is face centered
1630            fact=2 ; center=-3
1631            cell_base(:,ia)=2.0d0*vecta(:)
1632            cell_base(:,ib)=2.0d0*vectb(:)
1633            cell_base(:,itrial)=minim(:,itrial)
1634            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1635          end if
1636        end if
1637      end if
1638 
1639 !    DEBUG
1640 !    write(std_out,*)' after test_c, foundc=',foundc
1641 !    ENDDEBUG
1642 
1643 !    Working hypothesis : the conventional cell is orthogonal,
1644 !    and body centered with no basis vector being an axis,
1645 !    in which case the basis vectors must be equal (even for orthorhombic)
1646      if(foundc==0 .and. orthogonal==1 .and. &
1647 &     equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then
1648 !      Compute the combination of the two other vectors
1649        vecta(:)=minim(:,ia)+minim(:,ib)
1650        vectb(:)=minim(:,ia)-minim(:,ib)
1651        norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
1652        norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
1653 !      Project the trial vector on the first of the two vectors
1654        reduceda=( minim(1,itrial)*vecta(1)+       &
1655 &       minim(2,itrial)*vecta(2)+       &
1656 &       minim(3,itrial)*vecta(3) )/norm2a
1657        reducedb=( minim(1,itrial)*vectb(1)+       &
1658 &       minim(2,itrial)*vectb(2)+       &
1659 &       minim(3,itrial)*vectb(3) )/norm2b
1660        if( abs(abs(reduceda)-0.5d0)<tolsym )then
1661 !        The first vector is an axis
1662          fact=2 ; center=-1
1663          cell_base(:,ia)=vecta(:)
1664          vecta(:)=minim(:,itrial)-reduceda*vecta(:)
1665          vectb(:)=0.5d0*vectb(:)
1666          cell_base(:,ib)=vecta(:)+vectb(:)
1667          cell_base(:,itrial)=vecta(:)-vectb(:)
1668          call holocell(cell_base,0,foundc,iholohedry,tolsym)
1669        else if( abs(abs(reducedb)-0.5d0)<tolsym )then
1670 !        The second vector is an axis
1671          fact=2 ; center=-1
1672          cell_base(:,ib)=vectb(:)
1673          vectb(:)=minim(:,itrial)-reducedb*vectb(:)
1674          vecta(:)=0.5d0*vecta(:)
1675          cell_base(:,ia)=vectb(:)+vecta(:)
1676          cell_base(:,itrial)=vectb(:)-vecta(:)
1677          call holocell(cell_base,0,foundc,iholohedry,tolsym)
1678        end if
1679      end if
1680 
1681 !    Working hypothesis : the conventional cell is orthogonal,
1682 !    and face centered, in the case where two minimal vectors are equal
1683      if(foundc==0 .and. orthogonal==1 .and. equal(itrial)==1 ) then
1684 !      Compute the combination of these two vectors
1685        vecta(:)=minim(:,ia)+minim(:,ib)
1686        vectb(:)=minim(:,ia)-minim(:,ib)
1687        norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
1688        norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
1689 !      Project the trial vector on the two vectors
1690        reduceda=( minim(1,itrial)*vecta(1)+       &
1691 &       minim(2,itrial)*vecta(2)+       &
1692 &       minim(3,itrial)*vecta(3) )/norm2a
1693        reducedb=( minim(1,itrial)*vectb(1)+       &
1694 &       minim(2,itrial)*vectb(2)+       &
1695 &       minim(3,itrial)*vectb(3) )/norm2b
1696        if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. &
1697 &       ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym)       )then
1698          fact=2 ; center=-3
1699          cell_base(:,itrial)= &
1700 &         (minim(:,itrial)-reduceda*vecta(:)-reducedb*vectb(:) )*2.0d0
1701          cell_base(:,ia)=vecta(:)
1702          cell_base(:,ib)=vectb(:)
1703          call holocell(cell_base,0,foundc,iholohedry,tolsym)
1704        end if
1705      end if
1706 
1707 !    Working hypothesis : the conventional cell is orthogonal,
1708 !    face centered, but no two vectors are on the same "square"
1709      if(foundc==0 .and. orthogonal==1)then
1710 !      Compute the combination of these two vectors
1711        vecta(:)=minim(:,ia)+minim(:,ib)
1712        vectb(:)=minim(:,ia)-minim(:,ib)
1713        norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
1714        norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
1715 !      The trial vector length must be equal to one of these lengths
1716        if(abs(metmin(itrial,itrial)-norm2a)<tolsym*norm2a)then
1717          fact=2 ; center=-3
1718          cell_base(:,ia)=vecta(:)+minim(:,itrial)
1719          cell_base(:,ib)=vecta(:)-minim(:,itrial)
1720 !        Project vectb perpendicular to cell_base(:,ia) and cell_base(:,ib)
1721          norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2
1722          norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2
1723          reduceda=( cell_base(1,ia)*vectb(1)+       &
1724 &         cell_base(2,ia)*vectb(2)+       &
1725 &         cell_base(3,ia)*vectb(3) )/norm2a
1726          reducedb=( cell_base(1,ib)*vectb(1)+       &
1727 &         cell_base(2,ib)*vectb(2)+       &
1728 &         cell_base(3,ib)*vectb(3) )/norm2b
1729          if( abs(abs(reduceda)-0.5d0)<tolsym .and.         &
1730 &         abs(abs(reducedb)-0.5d0)<tolsym      )then
1731            cell_base(:,itrial)=vectb(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib)
1732            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1733          end if
1734        else if(abs(metmin(itrial,itrial)-norm2b)<tolsym*norm2b)then
1735          fact=2 ; center=-3
1736          cell_base(:,ia)=vectb(:)+minim(:,itrial)
1737          cell_base(:,ib)=vectb(:)-minim(:,itrial)
1738 !        Project vecta perpendicular to cell_base(:,ia) and cell_base(:,ib)
1739          norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2
1740          norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2
1741          reduceda=( cell_base(1,ia)*vecta(1)+       &
1742 &         cell_base(2,ia)*vecta(2)+       &
1743 &         cell_base(3,ia)*vecta(3) )/norm2a
1744          reducedb=( cell_base(1,ib)*vecta(1)+       &
1745 &         cell_base(2,ib)*vecta(2)+       &
1746 &         cell_base(3,ib)*vecta(3) )/norm2b
1747          if( abs(abs(reduceda)-0.5d0)<tolsym .and.         &
1748 &         abs(abs(reducedb)-0.5d0)<tolsym      )then
1749            cell_base(:,itrial)=vecta(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib)
1750            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1751          end if
1752        end if
1753      end if
1754 
1755 !    Working hypothesis : the cell is rhombohedral, and
1756 !    the three minimal vectors have same length and same absolute scalar product
1757      if(foundc==0 .and. iholohedry==5 .and. &
1758 &     equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then
1759        if( abs(abs(metmin(1,2))-abs(metmin(1,3)))<tolsym*metmin(1,1) .and.     &
1760 &       abs(abs(metmin(1,2))-abs(metmin(2,3)))<tolsym*metmin(2,2)      )then
1761          fact=1 ; center=0
1762          cell_base(:,:)=minim(:,:)
1763 !        One might have to change the sign of one of the vectors
1764          sign12=1 ; sign13=1 ; sign23=1
1765          if(metmin(1,2)<0.0d0)sign12=-1
1766          if(metmin(1,3)<0.0d0)sign13=-1
1767          if(metmin(2,3)<0.0d0)sign23=-1
1768          sumsign=sign12+sign13+sign23
1769          if(sumsign==-1)then
1770            if(sign12==1)cell_base(:,3)=-cell_base(:,3)
1771            if(sign13==1)cell_base(:,2)=-cell_base(:,2)
1772            if(sign23==1)cell_base(:,1)=-cell_base(:,1)
1773          else if(sumsign==1)then
1774            if(sign12==-1)cell_base(:,3)=-cell_base(:,3)
1775            if(sign13==-1)cell_base(:,2)=-cell_base(:,2)
1776            if(sign23==-1)cell_base(:,1)=-cell_base(:,1)
1777          end if
1778          call holocell(cell_base,0,foundc,iholohedry,tolsym)
1779        end if
1780      end if
1781 
1782 !    DEBUG
1783 !    write(std_out,*)' after test_3a, foundc=',foundc
1784 !    write(std_out,*)' after test_3a, itrial=',itrial
1785 !    write(std_out,*)' after test_3a, equal(:)=',equal(:)
1786 !    ENDDEBUG
1787 
1788 !    Working hypothesis : the cell is rhombohedral, one vector
1789 !    is parallel to the trigonal axis
1790      if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 )then
1791        vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib)
1792        norm2trial=minim(1,itrial)**2+minim(2,itrial)**2+minim(3,itrial)**2
1793        reduceda=( minim(1,itrial)*vecta(1)+       &
1794 &       minim(2,itrial)*vecta(2)+       &
1795 &       minim(3,itrial)*vecta(3) )/norm2trial
1796        reducedb=( minim(1,itrial)*vectb(1)+       &
1797 &       minim(2,itrial)*vectb(2)+       &
1798 &       minim(3,itrial)*vectb(3) )/norm2trial
1799 !      DEBUG
1800 !      write(std_out,*)' reduceda,reducedb=',reduceda,reducedb
1801 !      ENDDEBUG
1802        if(abs(abs(reduceda)-1.0d0/3.0d0)<tolsym .and.      &
1803 &       abs(abs(reducedb)-1.0d0/3.0d0)<tolsym      ) then
1804 !        Possibly change of sign to make positive the scalar product with
1805 !        the vector parallel to the trigonal axis
1806          if(reduceda<zero)vecta(:)=-vecta(:)
1807          if(reducedb<zero)vectb(:)=-vectb(:)
1808 !        Projection on the orthogonal plane
1809          vecta(:)=vecta(:)-abs(reduceda)*cell_base(:,itrial)
1810          vectb(:)=vectb(:)-abs(reducedb)*cell_base(:,itrial)
1811 !        These two vectors should have an angle of 120 degrees
1812          norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
1813          scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
1814 !        DEBUG
1815 !        write(std_out,*)' norm2a,scalarprod=',norm2a,scalarprod
1816 !        ENDDEBUG
1817          if(abs(two*scalarprod+norm2a)<tolsym*norm2a)then
1818            fact=1 ; center=0
1819            if(scalarprod>0.0d0)vectb(:)=-vectb(:)
1820 !          Now vecta and vectb have an angle of 120 degrees
1821            cell_base(:,1)=cell_base(:,itrial)/3.0d0+vecta(:)
1822            cell_base(:,2)=cell_base(:,itrial)/3.0d0+vectb(:)
1823            cell_base(:,3)=cell_base(:,itrial)/3.0d0-vecta(:)-vectb(:)
1824 !          DEBUG
1825 !          write(std_out,*)' cell_base(:,1)=',cell_base(:,1)
1826 !          write(std_out,*)' cell_base(:,2)=',cell_base(:,2)
1827 !          write(std_out,*)' cell_base(:,3)=',cell_base(:,3)
1828 !          ENDDEBUG
1829            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1830          end if
1831        end if
1832      end if
1833 
1834 !    Working hypothesis : the cell is rhombohedral, one vector
1835 !    is in the plane perpendicular to the trigonal axis
1836      if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then
1837        vecta(:)=minim(:,ia)+minim(:,ib)
1838        vectb(:)=minim(:,ia)-minim(:,ib)
1839        norm2trial=cell_base(1,itrial)**2+cell_base(2,itrial)**2+cell_base(3,itrial)**2
1840        norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
1841        norm2b=vecta(1)**2+vecta(2)**2+vecta(3)**2
1842        reduceda=( cell_base(1,itrial)*vecta(1)+       &
1843 &       cell_base(2,itrial)*vecta(2)+       &
1844 &       cell_base(3,itrial)*vecta(3) )/norm2trial
1845        reducedb=( cell_base(1,itrial)*vectb(1)+       &
1846 &       cell_base(2,itrial)*vectb(2)+       &
1847 &       cell_base(3,itrial)*vectb(3) )/norm2trial
1848        if(abs(norm2trial-norm2a)<tolsym*norm2a .and. &
1849 &       abs(abs(2*reduceda)-norm2trial)<tolsym*norm2trial    )then
1850          fact=1 ; center=0
1851          cell_base(:,1)=minim(:,ia)
1852          cell_base(:,2)=-minim(:,ib)
1853          cell_base(:,3)=-minim(:,ib)+2*reduceda*minim(:,itrial)
1854          call holocell(cell_base,0,foundc,iholohedry,tolsym)
1855        else if (abs(norm2trial-norm2b)<tolsym*norm2b .and. &
1856 &         abs(abs(2*reducedb)-norm2trial)<tolsym*norm2trial    )then
1857          fact=1 ; center=0
1858          cell_base(:,1)=minim(:,ia)
1859          cell_base(:,2)=minim(:,ib)
1860          cell_base(:,3)=minim(:,ib)+2*reducedb*minim(:,itrial)
1861          call holocell(cell_base,0,foundc,iholohedry,tolsym)
1862        end if
1863      end if
1864 
1865 !    Working hypothesis : the cell is rhombohedral, two vectors
1866 !    are in the plane perpendicular to the trigonal axis
1867      if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then
1868        vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib)
1869        norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
1870        norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
1871        scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
1872        if(abs(abs(2*scalarprod)-norm2a)<tolsym*norm2a)then
1873 !        This is in order to have 120 angle between vecta and vectb
1874          if(scalarprod>0.0d0)vectb(:)=-vectb(:)
1875          reduceda=( cell_base(1,itrial)*vecta(1)+        &
1876 &         cell_base(2,itrial)*vecta(2)+        &
1877 &         cell_base(3,itrial)*vecta(3) )/norm2a
1878          reducedb=( cell_base(1,itrial)*vectb(1)+        &
1879 &         cell_base(2,itrial)*vectb(2)+        &
1880 &         cell_base(3,itrial)*vectb(3) )/norm2b
1881          fact=1 ; center=0
1882          cell_base(:,1)=minim(:,itrial)
1883          if(abs(reduceda-0.5d0)<tolsym .and. abs(reducedb)<tolsym )then
1884            cell_base(:,2)=minim(:,itrial)-vecta(:)
1885            cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:)
1886            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1887          else if(abs(reduceda-0.5d0)<tolsym.and. abs(reducedb+0.5d0)<tolsym )then
1888            cell_base(:,2)=minim(:,itrial)-vecta(:)
1889            cell_base(:,3)=minim(:,itrial)+vectb(:)
1890            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1891          else if(abs(reduceda)<tolsym .and. abs(reducedb+0.5d0)<tolsym )then
1892            cell_base(:,2)=minim(:,itrial)+vectb(:)
1893            cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:)
1894            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1895          else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb)<tolsym)then
1896            cell_base(:,2)=minim(:,itrial)+vecta(:)
1897            cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:)
1898            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1899          else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb-0.5d0)<tolsym)then
1900            cell_base(:,2)=minim(:,itrial)+vecta(:)
1901            cell_base(:,3)=minim(:,itrial)-vectb(:)
1902            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1903          else if(abs(reduceda)<tolsym .and. abs(reducedb-0.5d0)<tolsym )then
1904            cell_base(:,2)=minim(:,itrial)-vectb(:)
1905            cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:)
1906            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1907          end if
1908        end if
1909      end if
1910 
1911 !    Working hypothesis : monoclinic holohedry, primitive. Then, two angles are 90 degrees
1912      if(foundc==0 .and. iholohedry==2 .and. &
1913 &     ang90(ia)==1 .and. ang90(ib)==1 ) then
1914        fact=1 ; center=0
1915        cell_base(:,1)=minim(:,ia)
1916        cell_base(:,2)=minim(:,itrial)
1917        cell_base(:,3)=minim(:,ib)
1918 !      Checks that the basis vectors are OK for the target holohedry
1919        call holocell(cell_base,0,foundc,iholohedry,tolsym)
1920      end if
1921 
1922 !    Monoclinic holohedry, one-face-centered cell
1923 !    Working hypothesis, two vectors have equal length.
1924      do icase=1,5
1925        if(foundc==0 .and. iholohedry==2 .and. equal(itrial)==1 ) then
1926          vecta(:)=cell_base(:,ia)+cell_base(:,ib)
1927          vectb(:)=cell_base(:,ia)-cell_base(:,ib)
1928          norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
1929          norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
1930 !        The minim(:,trial) vector belongs to the
1931 !        plane parallel to the cell_base(:,ia),cell_base(:,ib) plane
1932 !        In that plane, must try minim(:,itrial),
1933 !        as well as the 4 different combinations of
1934 !        minim(:,itrial) with the vectors in the plane
1935          if(icase==1)vectc(:)=minim(:,itrial)
1936          if(icase==2)vectc(:)=minim(:,itrial)+cell_base(:,ia)
1937          if(icase==3)vectc(:)=minim(:,itrial)+cell_base(:,ib)
1938          if(icase==4)vectc(:)=minim(:,itrial)-cell_base(:,ia)
1939          if(icase==5)vectc(:)=minim(:,itrial)-cell_base(:,ib)
1940          norm2c=vectc(1)**2+vectc(2)**2+vectc(3)**2
1941          sca=vectc(1)*vecta(1)+&
1942 &         vectc(2)*vecta(2)+&
1943 &         vectc(3)*vecta(3)
1944          scb=vectc(1)*vectb(1)+&
1945 &         vectc(2)*vectb(2)+&
1946 &         vectc(3)*vectb(3)
1947 !        DEBUG
1948 !        write(std_out,*)' symlatt : test iholohedry=2, sca,scb=',sca,scb
1949 !        ENDDEBUG
1950          if(abs(sca)<tolsym*sqrt(norm2c*norm2a) .or. abs(scb)<tolsym*sqrt(norm2c*norm2b))then
1951            fact=2 ; center=3
1952 !          The itrial direction is centered
1953            cell_base(:,3)=vectc(:)
1954            if(abs(sca)<tolsym*sqrt(norm2c*norm2a))then
1955              cell_base(:,2)=vecta(:)
1956              cell_base(:,1)=vectb(:)
1957              call holocell(cell_base,0,foundc,iholohedry,tolsym)
1958            else if(abs(scb)<tolsym*sqrt(norm2c*norm2b))then
1959              cell_base(:,2)=vectb(:)
1960              cell_base(:,1)=vecta(:)
1961              call holocell(cell_base,0,foundc,iholohedry,tolsym)
1962            end if
1963          end if
1964        end if
1965      end do ! icase=1,5
1966 
1967 !    Monoclinic holohedry, one-face-centered cell, but non equivalent.
1968 !    This case, one pair of vectors is orthogonal
1969      if(foundc==0 .and. iholohedry==2 .and. ang90(itrial)==1) then
1970        vecta(:)=minim(:,ia)
1971        vectb(:)=minim(:,ib)
1972        norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
1973        norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
1974 !      Project the trial vector on the two vectors
1975        reduceda=( minim(1,itrial)*vecta(1)+       &
1976 &       minim(2,itrial)*vecta(2)+       &
1977 &       minim(3,itrial)*vecta(3) )/norm2a
1978        reducedb=( minim(1,itrial)*vectb(1)+       &
1979 &       minim(2,itrial)*vectb(2)+       &
1980 &       minim(3,itrial)*vectb(3) )/norm2b
1981        if(abs(abs(reduceda)-0.5d0)<tolsym .or. abs(abs(reducedb)-0.5d0)<tolsym) then
1982          fact=2 ; center=3
1983          if(abs(abs(reduceda)-0.5d0)<tolsym)then
1984            cell_base(:,2)=vecta(:)
1985            cell_base(:,3)=vectb(:)
1986            cell_base(:,1)=2*(minim(:,itrial)-reduceda*vecta(:))
1987            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1988          else if(abs(abs(reducedb)-0.5d0)<tolsym)then
1989            cell_base(:,2)=vectb(:)
1990            cell_base(:,3)=vecta(:)
1991            cell_base(:,1)=2*(minim(:,itrial)-reducedb*vectb(:))
1992            call holocell(cell_base,0,foundc,iholohedry,tolsym)
1993          end if
1994        end if
1995      end if
1996 
1997 !    Monoclinic holohedry, one-face-centered cell, but non equivalent.
1998 !    This case, no pair of vectors is orthogonal, no pair of vector of equal lentgh
1999      if(foundc==0 .and. iholohedry==2)then
2000 !      Try to find a vector that belongs to the mediator plane, or the binary vector.
2001 !      There must be one such vector, if centered monoclinic and no pair of vectors of equal length,
2002 !      either among the three vectors, or among one of their differences or sums.
2003 !      And there must be, among the two other vectors, one vector whose projection
2004 !      on this vector is half the length of this vector.
2005        vecta(:)=minim(:,ia)
2006        vectb(:)=minim(:,ib)
2007 !      Try the different possibilities for the vector on which the projection will be half ...
2008        do ii=1,5
2009          if(ii==1)vectc(:)=minim(:,itrial)
2010          if(ii==2)vectc(:)=minim(:,itrial)+vecta(:)
2011          if(ii==3)vectc(:)=minim(:,itrial)-vecta(:)
2012          if(ii==4)vectc(:)=minim(:,itrial)+vectb(:)
2013          if(ii==5)vectc(:)=minim(:,itrial)-vectb(:)
2014          norm2trial=vectc(1)**2+vectc(2)**2+vectc(3)**2
2015 !        Project the two vectors on the trial vector
2016          reduceda=( vectc(1)*vecta(1)+vectc(2)*vecta(2)+vectc(3)*vecta(3) )/norm2trial
2017          reducedb=( vectc(1)*vectb(1)+vectc(2)*vectb(2)+vectc(3)*vectb(3) )/norm2trial
2018          found=0
2019          if(abs(abs(reduceda)-0.5d0)<tolsym)then
2020            vin1(:)=vectc(:)
2021            vin2(:)=2.0d0*(vecta(:)-reduceda*vectc(:))
2022            vext(:)=vectb(:)
2023            found=1
2024          else if(abs(abs(reducedb)-0.5d0)<tolsym)then
2025            vin1(:)=vectc(:)
2026            vin2(:)=2.0d0*(vectb(:)-reduceda*vectc(:))
2027            vext(:)=vecta(:)
2028            found=1
2029          end if
2030          if(found==1)exit
2031        end do
2032 !      Now, vin1 and vin2 are perpendicular to each other, and in the plane that contains the binary vector.
2033 !      One of them must be the binary vector if any.
2034 !      On the other hand, vext is out-of-plane. Might belong to the mediator plane or not.
2035 !      If C monoclinc, then the projection of this vext on the binary vector will be either 0 or +1/2 or -1/2.
2036 !      The binary axis must be stored in cell_base(:,2) for conventional C-cell
2037        if(found==1)then
2038          found=0
2039 
2040 !        Test vin1 being the binary axis
2041          norm2trial=vin1(1)**2+vin1(2)**2+vin1(3)**2
2042          reduceda=(vext(1)*vin1(1)+vext(2)*vin1(2)+vext(3)*vin1(3))/norm2trial
2043          if(abs(reduceda)<tolsym)then  ! vin1 is the binary axis and vext is in the mediator plane
2044            found=1
2045            cell_base(:,1)=vin2(:)
2046            cell_base(:,2)=vin1(:)
2047            cell_base(:,3)=vext(:)
2048          else if(abs(abs(reduceda)-0.5d0)<tolsym)then  ! vin1 is the binary axis and vext has +1/2 or -1/2 as projection
2049            found=1
2050            cell_base(:,1)=vin2(:)
2051            cell_base(:,2)=vin1(:)
2052            cell_base(:,3)=vext(:)-reduceda*vin1(:)+vin2(:)*half
2053          else
2054 !          Test vin2 being the binary axis
2055            norm2trial=vin2(1)**2+vin2(2)**2+vin2(3)**2
2056            reduceda=(vext(1)*vin2(1)+vext(2)*vin2(2)+vext(3)*vin2(3))/norm2trial
2057            if(abs(reduceda)<tolsym)then  ! vin2 is the binary axis and vext is in the mediator plane
2058              found=1
2059              cell_base(:,1)=vin1(:)
2060              cell_base(:,2)=vin2(:)
2061              cell_base(:,3)=vext(:)
2062            else if(abs(abs(reduceda)-0.5d0)<tolsym)then  ! vin2 is the binary axis and vext has +1/2 or -1/2 as projection
2063              found=1
2064              cell_base(:,1)=vin1(:)
2065              cell_base(:,2)=vin2(:)
2066              cell_base(:,3)=vext(:)-reduceda*vin2(:)+vin1(:)*half
2067            end if
2068          end if
2069 
2070          if(found==1)then
2071            fact=2 ; center=3
2072            call holocell(cell_base,0,foundc,iholohedry,tolsym)
2073          end if
2074        end if
2075      end if
2076 
2077    end do ! Do-loop on three different directions
2078  end do !  Do-loop on different target holohedries
2079 
2080  if(foundc==0)then
2081    iholohedry=1 ; fact=1 ; center=0
2082    cell_base(:,:)=minim(:,:)
2083  end if
2084 
2085 !DEBUG
2086 !write(std_out,*)' symlatt : done with centering tests, foundc=',foundc
2087 !write(std_out,*)'  center=',center
2088 !write(std_out,*)'  iholohedry=',iholohedry
2089 !ENDDEBUG
2090 
2091 !--------------------------------------------------------------------------
2092 !Final check on the Bravais lattice, using the basis vectors
2093 
2094 !Recompute the metric tensor
2095  if(foundc==1)then
2096    do ii=1,3
2097      metmin(:,ii)=cell_base(1,:)*cell_base(1,ii)+&
2098 &     cell_base(2,:)*cell_base(2,ii)+&
2099 &     cell_base(3,:)*cell_base(3,ii)
2100    end do
2101  end if
2102 
2103 !Examine the angles and vector lengths
2104  ang90(:)=0
2105  if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1
2106  if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1
2107  if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1
2108  equal(:)=0
2109  if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1
2110  if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1
2111  if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1
2112 
2113 !DEBUG
2114 !write(std_out,*)' symlatt : recompute the  metric tensor '
2115 !write(std_out,*)'  ang90=',ang90
2116 !write(std_out,*)'  equal=',equal
2117 !ENDDEBUG
2118 
2119 !The axes will be aligned with the previously determined
2120 !basis vectors, EXCEPT for the tetragonal cell, see later
2121  axes(:,:)=cell_base(:,:)
2122 
2123  found=0
2124 !Check orthogonal conventional cells
2125  if(ang90(1)+ang90(2)+ang90(3)==3)then
2126 
2127 !  Cubic system
2128    if(equal(1)+equal(2)+equal(3)==3)then
2129 !    However, one-face centered is not admitted
2130      if(center==0 .or. center==-1 .or. center==-3)then
2131        iholohedry=7 ; found=1
2132        if(center==0)then
2133          write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is cP (primitive cubic)'
2134        else if(center==-1)then
2135          write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is cI (body-centered cubic)'
2136        else if(center==-3)then
2137          write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is cF (face-centered cubic)'
2138        end if
2139      end if
2140    end if
2141 
2142 !  Tetragonal system
2143    if(found==0 .and. &
2144 &   (equal(1)==1 .or. equal(2)==1 .or. equal(3)==1) )then
2145 !    However, one-face centered or face-centered is not admitted
2146      if(center==0 .or. center==-1)then
2147        iholohedry=4 ; found=1
2148        if(equal(1)==1)then
2149          axes(:,3)=cell_base(:,1) ; axes(:,1)=cell_base(:,2) ; axes(:,2)=cell_base(:,3)
2150        else if(equal(2)==1)then
2151          axes(:,3)=cell_base(:,2) ; axes(:,2)=cell_base(:,1) ; axes(:,1)=cell_base(:,3)
2152        else if(equal(3)==1)then
2153          axes(:,:)=cell_base(:,:)
2154        end if
2155        if(center==0)then
2156          write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is tP (primitive tetragonal)'
2157        else if(center==-1)then
2158          write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is tI (body-centered tetragonal)'
2159        end if
2160      end if
2161    end if
2162 
2163 !  Orthorhombic system
2164    if(found==0)then
2165      iholohedry=3 ; found=1
2166      axes(:,:)=cell_base(:,:)
2167      if(center==0)then
2168        write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oP (primitive orthorhombic)'
2169      else if(center==-1)then
2170        write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oI (body-centered orthorhombic)'
2171      else if(center==1 .or. center==2 .or. center==3)then
2172        write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oC (one-face-centered orthorhombic)'
2173      else if(center==-3)then
2174        write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oF (face-centered orthorhombic)'
2175      end if
2176    end if
2177 
2178  else
2179 
2180 !  Hexagonal system
2181    if(found==0 .and. ang90(1)==1 .and. ang90(2)==1 .and. equal(3)==1 .and. (2*metmin(2,1)+metmin(1,1))<tolsym*metmin(1,1))then
2182      iholohedry=6 ; found=1
2183      write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is hP (primitive hexagonal)'
2184    end if
2185 
2186 !  Rhombohedral system
2187    if(found==0 .and. equal(1)+equal(2)+equal(3)==3 .and.       &
2188 &   abs(metmin(2,1)-metmin(3,2))<tolsym*metmin(2,2)             .and.       &
2189 &   abs(metmin(2,1)-metmin(3,1))<tolsym*metmin(1,1) )then
2190      iholohedry=5 ; found=1
2191      write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is hR (rhombohedral)'
2192    end if
2193 
2194 !  Monoclinic system
2195    if(found==0 .and. ang90(1)+ang90(2)+ang90(3)==2 )then
2196      iholohedry=2 ; found=1
2197      if(center==0)then
2198        write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is mP (primitive monoclinic)'
2199      else if(center==3)then
2200        write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is mC (one-face-centered monoclinic)'
2201      end if
2202    end if
2203 
2204 !  Triclinic system
2205    if(found==0)then
2206      iholohedry=1 ; found=1
2207      write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is aP (primitive triclinic)'
2208    end if
2209 
2210  end if
2211 
2212  call wrtout(std_out,message,'COLL')
2213 
2214 !--------------------------------------------------------------------------
2215 !Make sure that axes form a right-handed coordinate system
2216 !(Note : this should be done in the body of the routine,
2217 !by making changes that leave the sign of the mixed product of the three
2218 !vectors invariant)
2219  determinant=axes(1,1)*axes(2,2)*axes(3,3) &
2220 & +axes(1,2)*axes(2,3)*axes(3,1) &
2221 & +axes(1,3)*axes(3,2)*axes(2,1) &
2222 & -axes(1,1)*axes(3,2)*axes(2,3) &
2223 & -axes(1,3)*axes(2,2)*axes(3,1) &
2224 & -axes(1,2)*axes(2,1)*axes(3,3)
2225  if(determinant<0.0d0)then
2226    axes(:,:)=-axes(:,:)
2227  end if
2228 
2229 !--------------------------------------------------------------------------
2230 !Prefer symmetry axes on the same side as the primitive axes,
2231 !when the changes are allowed
2232  do
2233    do ia=1,3
2234      scprods(ia,:)=axes(1,ia)*rprimd(1,:)+&
2235 &     axes(2,ia)*rprimd(2,:)+&
2236 &     axes(3,ia)*rprimd(3,:)
2237      norm2trial=sum(axes(:,ia)**2)
2238      scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial)
2239    end do
2240    do ia=1,3
2241      norm2trial=sum(rprimd(:,ia)**2)
2242      scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial)
2243    end do
2244 
2245 !  One should now try all the generators of the
2246 !  proper rotations of each Bravais lattice, coupled with change of
2247 !  signs of each vector. This is not done systematically in what follows ...
2248 !  Here, the third axis is left unchanged
2249    if(iholohedry/=5)then
2250      if(scprods(1,1)<-tolsym .and. scprods(2,2)<-tolsym)then
2251        axes(:,1)=-axes(:,1) ; axes(:,2)=-axes(:,2)
2252        cycle
2253      end if
2254    end if
2255 !  The first (or second) axis is left unchanged
2256    if(iholohedry/=5 .and. iholohedry/=6)then
2257      if(scprods(2,2)<-tolsym .and. scprods(3,3)<-tolsym)then
2258        axes(:,2)=-axes(:,2) ; axes(:,3)=-axes(:,3)
2259        cycle
2260      end if
2261      if(scprods(1,1)<-tolsym .and. scprods(3,3)<-tolsym)then
2262        axes(:,1)=-axes(:,1) ; axes(:,3)=-axes(:,3)
2263        cycle
2264      end if
2265    end if
2266 !  Permutation of the three axis
2267    if(iholohedry==5 .or. iholohedry==7)then
2268      trace=scprods(1,1)+scprods(2,2)+scprods(3,3)
2269      if(trace+tolsym< scprods(1,2)+scprods(2,3)+scprods(3,1))then
2270        vecta(:)=axes(:,1) ; axes(:,1)=axes(:,3)
2271        axes(:,3)=axes(:,2); axes(:,2)=vecta(:)
2272        cycle
2273      end if
2274      if(trace+tolsym < scprods(1,3)+scprods(2,1)+scprods(3,2))then
2275        vecta(:)=axes(:,1) ; axes(:,1)=axes(:,2)
2276        axes(:,2)=axes(:,3); axes(:,3)=vecta(:)
2277        cycle
2278      end if
2279 !    This case is observed when the three new vectors
2280 !    are pointing opposite to the three original vectors
2281 !    One takes their opposite, then switch to of them, then process
2282 !    them again in the loop
2283      if(sum(scprods(:,:))<tolsym)then
2284        axes(:,1)=-axes(:,1)
2285        vecta(:)=-axes(:,2)
2286        axes(:,2)=-axes(:,3)
2287        axes(:,3)=vecta(:)
2288        cycle
2289      end if
2290    end if
2291    exit
2292 !  Other cases might be coded ...
2293  end do
2294 
2295 !--------------------------------------------------------------------------
2296 
2297 !DEBUG
2298 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',&
2299 !&  rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
2300 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes  =',&
2301 !&  axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
2302 !ENDDEBUG
2303 
2304 !Compute the coordinates of rprimd in the system defined by axes(:,:)
2305  call matr3inv(axes,axesinvt)
2306  do ii=1,3
2307    coord(:,ii)=rprimd(1,ii)*axesinvt(1,:)+ &
2308 &   rprimd(2,ii)*axesinvt(2,:)+ &
2309 &   rprimd(3,ii)*axesinvt(3,:)
2310  end do
2311 
2312 !Check that the coordinates are integers, or half-integer in
2313 !the case there is a centering, and generate integer coordinates
2314  do ii=1,3
2315    do jj=1,3
2316      val=coord(ii,jj)*fact
2317      if(abs(val-nint(val))>fact*tolsym)then
2318        write(message,'(4a,a,3es18.10,a,a,3es18.10,a,a,3es18.10,a,a,i4)')&
2319 &       'One of the coordinates of rprimd in axes is non-integer,',ch10,&
2320 &       'or non-half-integer (if centering).',ch10,&
2321 &       'coord=',coord(:,1),ch10,&
2322 &       '      ',coord(:,2),ch10,&
2323 &       '      ',coord(:,3),ch10,&
2324 &       'fact=',fact
2325        MSG_BUG(message)
2326      end if
2327      icoord(ii,jj)=nint(val)
2328    end do
2329  end do
2330 
2331 !Store the bravais lattice characteristics
2332  bravais(1)=iholohedry
2333  bravais(2)=center
2334  bravais(3:5)=icoord(1:3,1)
2335  bravais(6:8)=icoord(1:3,2)
2336  bravais(9:11)=icoord(1:3,3)
2337 
2338 !--------------------------------------------------------------------------
2339 !Initialize the set of symmetries
2340 !Bravais lattices are always invariant under identity and inversion
2341 
2342 !Identity and inversion
2343  ptsymrel(:,:,1)=identity(:,:) ; ptsymrel(:,:,2)=-identity(:,:)
2344  nptsym=2
2345 
2346 !Keep this for IFCv70 compiler
2347  if(nptsym/=2)then
2348    write(message,'(a,a,a,a)')ch10,&
2349 &   ' symlatt : BUG -',ch10,&
2350 &   '  Crazy error, compiler bug '
2351    call wrtout(std_out,message,'COLL')
2352  end if
2353 
2354 !--------------------------------------------------------------------------
2355 !Initialize some generators
2356 !gen6 is defined in a coordinated system with gamma=120 degrees
2357  gen6(:,:)=0  ; gen6(3,3)=1  ; gen6(1,1)=1  ; gen6(1,2)=-1 ; gen6(2,1)=1
2358  gen3(:,:)=0  ; gen3(1,2)=1  ; gen3(2,3)=1  ; gen3(3,1)=1
2359  gen2xy(:,:)=0 ; gen2xy(2,1)=1 ; gen2xy(1,2)=1; gen2xy(3,3)=1
2360  gen2y(:,:)=0 ; gen2y(1,1)=-1; gen2y(2,2)=1 ; gen2y(3,3)=-1
2361  gen2z(:,:)=0 ; gen2z(1,1)=-1; gen2z(2,2)=-1; gen2z(3,3)=1
2362 
2363 !--------------------------------------------------------------------------
2364 
2365 !Define the generators for each holohedry (inversion is already included)
2366  if(iholohedry==6)then
2367    ngen=2
2368    gen(:,:,1)=gen2xy(:,:) ; order(1)=2
2369    gen(:,:,2)=gen6(:,:)   ; order(2)=6
2370  else if(iholohedry==5)then
2371    ngen=2
2372    gen(:,:,1)=gen2xy(:,:) ; order(1)=2
2373    gen(:,:,2)=gen3(:,:)   ; order(2)=3
2374  else
2375    gen(:,:,1)=gen2y(:,:)  ; order(1)=2
2376    gen(:,:,2)=gen2z(:,:)  ; order(2)=2
2377    gen(:,:,3)=gen2xy(:,:) ; order(3)=2
2378    gen(:,:,4)=gen3(:,:)   ; order(4)=3
2379    if(iholohedry<=4)ngen=iholohedry-1
2380    if(iholohedry==7)ngen=4
2381  end if
2382 
2383 !Build the point symmetry operations from generators, in the reduced system
2384 !of coordinates defined by axes(:,:)
2385  if(ngen/=0)then
2386    do igen=1,ngen
2387      do isym=1+nptsym,order(igen)*nptsym
2388        jsym=isym-nptsym
2389        do ii=1,3
2390          ptsymrel(:,ii,isym)=gen(:,1,igen)*ptsymrel(1,ii,jsym)+ &
2391 &         gen(:,2,igen)*ptsymrel(2,ii,jsym)+ &
2392 &         gen(:,3,igen)*ptsymrel(3,ii,jsym)
2393        end do
2394      end do
2395      nptsym=order(igen)*nptsym
2396 
2397    end do
2398  end if
2399 
2400 !--------------------------------------------------------------------------
2401 
2402 !Transform symmetry matrices in the system defined by rprimd
2403  call symrelrot(nptsym,axes,rprimd,ptsymrel,tolsym)
2404 
2405 end subroutine symlatt

m_symfind/symspgr [ Functions ]

[ Top ] [ m_symfind ] [ Functions ]

NAME

 symspgr

FUNCTION

 Using the type of each symmetry operation (found in symplanes.f and symaxes.f):
 proper symmetries 1,2,2_1,3,3_1,3_2,4,4_1,4_2,4_3,6,6_1,...6_5
 improper symmetries -1,m,a,b,c,d,n,g,-3,-4,-6 ,
 build an array with the number of such operations. then, call symlist.f to identify the space group.
 The identification is not unambiguous still ...

INPUTS

 bravais(11): bravais(1)=iholohedry
              bravais(2)=center
              bravais(3:11)=coordinates of rprimd in the axes
              of the conventional bravais lattice (*2 if center/=0)
 nsym=actual number of symmetries
 symrel(3,3,nsym)= nsym symmetry operations in real space in terms
   of primitive translations
 tnons(3,nsym)=nonsymmorphic translations for each symmetry (would
   be 0 0 0 each for a symmorphic space group)

OUTPUT

 spgroup=symmetry space group number

NOTES

 It is assumed that the symmetry operations will be entered in the
 symrel tnons arrays, for the PRIMITIVE cell. The matrix of transformation
 from the primitive cell to the conventional cell is described
 in the array "bravais" (see symlatt.F90).
 The present routine first make the transformation from the
 primitive coordinates to the conventional ones, then eventually
 generate additional symmetries, taking into account the
 centering translations.
 Then, the order and determinant of each symmetry operation
 is determined.

 For proper symmetries (rotations), the
 associated translation is also determined.
 However, left or right handed screw rotations are
 not (presently) distinguished, and will be attributed equally
 to left or right.

 For the detailed description of the labelling of the axes,
 see symaxes.f and symplanes.f

PARENTS

      symanal

CHILDREN

      spgdata,symcharac,symdet,symlist_bcc,symlist_fcc,symlist_others
      symlist_prim,symrelrot,wrtout,xred2xcart

SOURCE

 993 subroutine symspgr(bravais,nsym,spgroup,symrel,tnons,tolsym)
 994 
 995  use m_numeric_tools, only : OPERATOR(.x.)
 996 
 997 !This section has been created automatically by the script Abilint (TD).
 998 !Do not modify the following lines by hand.
 999 #undef ABI_FUNC
1000 #define ABI_FUNC 'symspgr'
1001 !End of the abilint section
1002 
1003  implicit none
1004 
1005 !Arguments ------------------------------------
1006 !scalars
1007  integer,intent(in) :: nsym
1008  integer,intent(out) :: spgroup
1009  real(dp),intent(in) :: tolsym
1010 !arrays
1011  integer,intent(in) :: bravais(11),symrel(3,3,nsym)
1012  real(dp),intent(inout) :: tnons(3,nsym)
1013 
1014 !Local variables-------------------------------
1015 !scalars
1016 ! logical,parameter :: verbose=.FALSE.
1017  integer :: additional_info,brvltt,center,direction=0,found,iholohedry,ii
1018  integer :: ishift,isym,jj,nshift,nsymconv,spgaxor,spgorig,sporder
1019  character(len=1) :: brvsb
1020  character(len=15) :: intsb,ptintsb,ptschsb,schsb
1021  character(len=35) :: intsbl
1022  character(len=500) :: message
1023  character(len=128) :: label
1024 !arrays
1025  integer :: n_axes(31),n_axest(31),prime(5),test_direction(3),symrel_uni(3,3)
1026  integer :: uniaxis(3),uniaxis_try(3)
1027  integer,allocatable :: determinant(:),symrelconv(:,:,:),t_axes(:)
1028  real(dp) :: axes(3,3),rprimdconv(3,3),trialt(3),vect(3,3)
1029  real(dp),allocatable :: shift(:,:),tnonsconv(:,:)
1030 
1031 !**************************************************************************
1032 
1033  DBG_ENTER("COLL")
1034 
1035 !Initialize brvltt, from bravais(2) and bravais(1)
1036  center=bravais(2)
1037  iholohedry=bravais(1)
1038  brvltt=1
1039  if(center==-1)brvltt=2  ! Inner centering
1040  if(center==-3)brvltt=3  ! Face centering
1041  if(center==1)brvltt=5  ! A-Face centering
1042  if(center==2)brvltt=6  ! B-Face centering
1043  if(center==3)brvltt=4  ! C-Face centering
1044  if(iholohedry==5)brvltt=7  ! Rhombohedral
1045 
1046 !Produce the symmetry operations, in the axis of the conventional cell
1047  nsymconv=nsym
1048  if(center/=0)nsymconv=2*nsymconv
1049  if(center==-3)nsymconv=4*nsym
1050  ABI_ALLOCATE(symrelconv,(3,3,nsymconv))
1051  ABI_ALLOCATE(tnonsconv,(3,nsymconv))
1052 
1053 !Produce symrel and tnons in conventional axes,
1054 !name them symrelconv and tnonsconv
1055  rprimdconv(:,1)=bravais(3:5)
1056  rprimdconv(:,2)=bravais(6:8)
1057  rprimdconv(:,3)=bravais(9:11)
1058 
1059  if(center/=0)rprimdconv(:,:)=rprimdconv(:,:)*half
1060 
1061  axes(:,:)=zero
1062  axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one
1063  symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym)
1064 !Note that the number of symmetry operations is still nsym
1065  call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym)
1066  call xred2xcart(nsym,rprimdconv,tnonsconv,tnons)
1067 !Gives the associated translation, with components in the
1068 !interval ]-0.5,0.5] .
1069  tnonsconv(:,1:nsym)=tnonsconv(:,1:nsym)-nint(tnonsconv(:,1:nsym)-tol6)
1070 
1071 !If the Bravais lattice is centered, duplicate or quadruplicate
1072 !the number of symmetry operations, using the Bravais
1073 !lattice shifts
1074  nshift=1
1075  if(center/=0)nshift=2
1076  if(center==-3)nshift=4
1077  ABI_ALLOCATE(shift,(3,nshift))
1078  shift(:,1)=zero
1079  if(center/=0 .and. center/=-3)then
1080    shift(:,2)=half
1081    if(center==1)shift(1,2)=zero
1082    if(center==2)shift(2,2)=zero
1083    if(center==3)shift(3,2)=zero
1084  else if(center==-3)then
1085    shift(:,2)=half ; shift(1,2)=zero
1086    shift(:,3)=half ; shift(2,3)=zero
1087    shift(:,4)=half ; shift(3,4)=zero
1088  end if ! center/=0 or -3
1089  if(nshift/=1)then
1090    do ishift=2,nshift
1091      symrelconv(:,:,(ishift-1)*nsym+1:ishift*nsym)=symrelconv(:,:,1:nsym)
1092      do isym=1,nsym
1093        tnonsconv(:,(ishift-1)*nsym+isym)=tnonsconv(:,isym)+shift(:,ishift)
1094      end do
1095    end do ! ishift
1096  end if ! nshift/=1
1097 
1098 !At this stage, all the symmetry operations are available,
1099 !expressed in the conventional axis, and also include
1100 !the Bravais lattive translations, and associated operations...
1101 
1102  n_axes(:)=0
1103 
1104  ABI_ALLOCATE(determinant,(nsymconv))
1105 
1106 !Get the determinant
1107  call symdet(determinant,nsymconv,symrelconv)
1108 
1109 !Get the order of each the symmetry operation, as well as the maximal order
1110 !Also, examine whether each symmetry operation is the inversion, or a root
1111 !of the inversion (like -3)
1112 !Decide which kind of point symmetry operation it is
1113 !Finally assign tnonsconv order and decide the space symmetry operation
1114 
1115  ABI_ALLOCATE(t_axes,(nsymconv))
1116 
1117  do isym=1,nsymconv
1118 
1119    call symcharac(center, determinant(isym), iholohedry, isym, label, &
1120    symrelconv(:,:,isym), tnonsconv(:,isym), t_axes(isym))
1121    if (t_axes(isym) == -1) then
1122      write(message, '(a,a,i3,a,3(a,3i4,a),a,3es22.12,a,a,3es22.12)' )ch10,&
1123 &     ' symspgr : problem with isym=',isym,ch10,&
1124 &     '  symrelconv(:,1,isym)=',symrelconv(:,1,isym),ch10,&
1125 &     '  symrelconv(:,2,isym)=',symrelconv(:,2,isym),ch10,&
1126 &     '  symrelconv(:,3,isym)=',symrelconv(:,3,isym),ch10,&
1127 &     '  tnonsconv(:,isym)=',tnonsconv(:,isym),ch10,&
1128 &     '  trialt(:)=',trialt(:)
1129      call wrtout(std_out,message,'COLL')
1130      write(message, '(a,i4,2a)' )&
1131 &     'The space symmetry operation number',isym,ch10,&
1132 &     'is not a (translated) root of unity'
1133      MSG_BUG(message)
1134    else if (t_axes(isym) == -2) then
1135      write(message, '(a,i0,a)' )'The symmetry operation number ',isym,' is not a root of unity'
1136      MSG_BUG(message)
1137    end if
1138 
1139    n_axes(t_axes(isym))=n_axes(t_axes(isym))+1
1140 
1141  end do ! isym=1,nsymconv
1142 
1143  if (sum(n_axes)-nsymconv/=0) then
1144    write(message, '(7a)' )&
1145 &   'Not all the symmetries have been recognized. ',ch10,&
1146 &   'This might be due either to an error in the input file',ch10,&
1147 &   'or to a BUG in ABINIT',ch10,&
1148 &   'Please contact the ABINIT group.'
1149    MSG_WARNING(message)
1150  end if
1151 
1152 !DEBUG
1153 !write(std_out,*)' symspgr : brvltt,nsymconv=',brvltt,nsymconv
1154 !write(std_out,*)' n_axes(1:10)=',n_axes(1:10)
1155 !write(std_out,*)' n_axes(11:20)=',n_axes(11:20)
1156 !write(std_out,*)' n_axes(21:31)=',n_axes(21:31)
1157 !ENDDEBUG
1158 
1159 !Treat cases in which the space group cannot be identified on the
1160 !basis of n_axes one need additional informations
1161  if(brvltt==1)then
1162 !  If the bravais lattice is primitive
1163    if(nsymconv==4)then
1164      n_axest=(/0,0,0,0,0,0,0,1,1,0,  0,0,0,0,0,2,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,0/)
1165      if(sum((n_axes-n_axest)**2)==0)then    ! Spgroup 27 (Pcc2) or 32 (Pba2)
1166        write(std_out,*)' symspgr: 27 or 32'
1167        additional_info=2
1168 !      Select binary axis
1169        do isym=1,nsymconv
1170          if(t_axes(isym)==8)then
1171 !          Find direction of binary axis
1172            if(symrelconv(1,1,isym)==1)direction=1
1173            if(symrelconv(2,2,isym)==1)direction=2
1174            if(symrelconv(3,3,isym)==1)direction=3
1175          end if
1176        end do
1177 !      Examine the projection of the translation vector of the a, b or c mirror planes
1178 !      onto the binary axis
1179        do isym=1,nsymconv
1180          if(t_axes(isym)==16)then
1181            if(abs(tnonsconv(direction,isym))>tol8)additional_info=1
1182          end if
1183        end do
1184      end if
1185    else if(nsymconv==8)then
1186      n_axest=(/0,0,0,0,1,0,0,1,1,0,  0,0,0,0,1,2,0,0,0,2,  0,0,0,0,0,0,0,0,0,0,0/)
1187      if(sum((n_axes-n_axest)**2)==0)then    ! Spgroup 55 (Pbam) or 57 (Pbcm)
1188        write(std_out,*)' symspgr: 55 or 57'
1189        additional_info=1
1190 !      Select mirror plane m
1191        do isym=1,nsymconv
1192          if(t_axes(isym)==15)then
1193 !          Find direction of mirror plane
1194            if(symrelconv(1,1,isym)==-1)direction=1
1195            if(symrelconv(2,2,isym)==-1)direction=2
1196            if(symrelconv(3,3,isym)==-1)direction=3
1197          end if
1198        end do
1199 !      Examine the projection of the translation vector of the a, b, or c mirror planes
1200 !      onto the binary axis
1201        do isym=1,nsymconv
1202          if(t_axes(isym)==16)then
1203            if(abs(tnonsconv(direction,isym))>tol8)additional_info=2
1204          end if
1205        end do
1206      end if
1207      n_axest=(/0,0,0,0,1,0,0,1,1,0,  0,0,0,0,0,2,0,1,0,2,  0,0,0,0,0,0,0,0,0,0,0/)
1208      if(sum((n_axes-n_axest)**2)==0)then    ! Spgroup 56 (Pccn) or 60 (Pbcn)
1209        write(std_out,*)' symspgr: 56 or 60'
1210        additional_info=1
1211 !      Select mirror plane n
1212        do isym=1,nsymconv
1213          if(t_axes(isym)==18)then
1214 !          Find direction of mirror plane
1215            if(symrelconv(1,1,isym)==-1)direction=1
1216            if(symrelconv(2,2,isym)==-1)direction=2
1217            if(symrelconv(3,3,isym)==-1)direction=3
1218          end if
1219        end do
1220 !      Examine the projection of the translation vector of the a, b, or c mirror planes
1221 !      onto the binary axis
1222        do isym=1,nsymconv
1223          if(t_axes(isym)==16)then
1224            if(abs(tnonsconv(direction,isym))<tol8)additional_info=2
1225          end if
1226        end do
1227      end if
1228    end if
1229  else if(brvltt==2)then
1230 !  In the few next lines, use additional_info as a flag
1231    additional_info=0
1232 !  If the bravais lattice is inner-centered
1233    if(nsymconv==8)then
1234 !    Test spgroup 23 (I222) or 24 (I2_{1}2_{1}2_{1})
1235      n_axest=(/0,0,0,0,0,0,1,1,3,0,  0,0,0,0,0,0,0,0,0,3,  0,0,0,0,0,0,0,0,0,0,0/)
1236      if(sum((n_axes-n_axest)**2)==0) additional_info=1
1237    else if(nsymconv==24)then
1238 !    Test spgroup 197 (I23) or 199 (I2_{1}3)
1239      n_axest=(/0,0,0,0,0,0,1,1,3,16, 0,0,0,0,0,0,0,0,0,3,  0,0,0,0,0,0,0,0,0,0,0/)
1240      if(sum((n_axes-n_axest)**2)==0) additional_info=1
1241    end if
1242    if(additional_info==1)then
1243      write(std_out,*)' symspgr: (23 or 24) or (197 or 199)'
1244 !    Select the three binary axes (they might be 2 or 2_1 !)
1245      test_direction(:)=0
1246      do isym=1,nsymconv
1247        if(t_axes(isym)==20)then
1248 !        Find direction of axis
1249          do direction=1,3
1250            if(symrelconv(direction,direction,isym)==1)then
1251              test_direction(direction)=1
1252              if(abs(tnonsconv(direction,isym))<tol8)then
1253                vect(:,direction)=tnonsconv(:,isym)
1254              else
1255                vect(:,direction)=tnonsconv(:,isym)+half
1256              end if
1257              vect(:,direction)=vect(:,direction)-nint(vect(:,direction)-tol8)
1258              vect(direction,direction)=zero
1259            end if
1260          end do ! direction=1,3
1261        end if ! if binary axis
1262      end do ! isym
1263      if(test_direction(1)/=1 .or. test_direction(2)/=1 .and. test_direction(3)/=1)then
1264        write(message, '(5a,3i4)' )&
1265 &       'For space groups 23, 24, 197 or 197, the three binary axes',ch10,&
1266 &       'are not equally partitioned along the x, y and z directions',ch10,&
1267 &       'test_direction(1:3)=',test_direction(:)
1268        MSG_BUG(message)
1269      end if
1270      additional_info=1
1271      if(abs(vect(1,2)-vect(1,3))>tol8 .or. &
1272 &     abs(vect(2,1)-vect(2,3))>tol8 .or. &
1273 &     abs(vect(3,1)-vect(3,2))>tol8) additional_info=2
1274    end if ! additional informations are needed
1275  end if ! brvltt==1
1276 
1277  if (brvltt==0 .or. brvltt==1) then ! Primitive
1278    call symlist_prim(additional_info,nsymconv,n_axes,spgroup)
1279  else if(brvltt==2)then
1280    call symlist_bcc(additional_info,nsymconv,n_axes,spgroup)
1281  else if(brvltt==3)then
1282    call symlist_fcc(nsymconv,n_axes,spgroup)
1283  else
1284    call symlist_others(brvltt,nsymconv,n_axes,spgroup)
1285  end if
1286 
1287  if(spgroup==0) then
1288    write(message, '(a,a,a,a,a)' )&
1289 &   'Could not find the space group.',ch10,&
1290 &   'This often happens when the user selects a restricted set of symmetries ',ch10,&
1291 &   'in the input file, instead of letting the code automatically find symmetries.'
1292    MSG_WARNING(message)
1293  end if
1294 
1295  spgorig=1 ; spgaxor=1
1296  call spgdata(brvsb,intsb,intsbl,ptintsb,ptschsb,schsb,spgaxor,spgroup,sporder,spgorig)
1297 
1298  if(spgroup/=0)then
1299    write(message, '(a,i4,2x,a,a,a,a,a)' ) ' symspgr : spgroup=',spgroup,trim(brvsb),trim(intsb),'   (=',trim(schsb),')'
1300    call wrtout(std_out,message,'COLL')
1301  end if
1302 
1303  if(bravais(1)==7)then
1304    write(message, '(a)' ) ' symspgr : optical characteristics = isotropic '
1305    call wrtout(std_out,message,'COLL')
1306  else if(bravais(1)==4 .or. bravais(1)==5 .or. bravais(1)==6)then
1307    write(message, '(a)' ) ' symspgr : optical characteristics = uniaxial '
1308    call wrtout(std_out,message,'COLL')
1309 !  Identify the first symmetry operation that is order 3, 4 or 6
1310    found=0
1311    do isym=1,nsym
1312 !    Proper rotations
1313      if( minval( abs( t_axes(isym)-(/10,12,14,22,23,24,25,26,27,28,29,30,31/) ))==0) then
1314        found=1 ; exit
1315 !   Improper symmetry operations
1316      else if( minval( abs( t_axes(isym)-(/1,2,3/) ))==0) then
1317        found=-1 ; exit
1318      end if
1319    end do
1320    if(found==-1 .or. found==1)then
1321      symrel_uni=symrel(:,:,isym)
1322      if(found==-1)symrel_uni=-symrel_uni
1323 !    Now, symrel_uni is a rotation of order 3, 4, 6, for which the axis must be identified
1324 !    It is actually the only eigenvector with eigenvalue 1. It can be found by cross products
1325 !    Subtract the unit matrix.
1326      do ii=1,3
1327        symrel_uni(ii,ii)=symrel_uni(ii,ii)-1
1328      end do
1329      found=0
1330      do ii=1,3
1331        jj=ii+1 ; if(jj==4)jj=1
1332 !      Cross product
1333        uniaxis=symrel_uni(ii,:).x.symrel_uni(jj,:)
1334        if(sum(uniaxis**2)/=0)then
1335          found=1 ; exit
1336        end if
1337      end do
1338      if(found==1)then
1339 !      Try to reduce the length, by an integer factor (try only primes 2, 3, 5, 7, 11)
1340        prime=(/2,3,5,7,11/)
1341        ii=1
1342        do while (ii<6)
1343          uniaxis_try=uniaxis/prime(ii)
1344          if(sum(abs(uniaxis_try*prime(ii)-uniaxis))==0)then
1345            uniaxis=uniaxis_try
1346          else
1347            ii=ii+1
1348          end if
1349        end do
1350        write(message, '(a,3i4)' ) ' Optical axis (in reduced coordinates, real space ) :',uniaxis
1351      end if
1352    end if
1353    if(found==0)then
1354      write(message, '(a)' ) ' However, the axis has not been found. Sorry for this.'
1355    end if
1356    call wrtout(std_out,message,'COLL')
1357  end if
1358 
1359  ABI_DEALLOCATE(determinant)
1360  ABI_DEALLOCATE(shift)
1361  ABI_DEALLOCATE(symrelconv)
1362  ABI_DEALLOCATE(tnonsconv)
1363  ABI_DEALLOCATE(t_axes)
1364 
1365  DBG_EXIT("COLL")
1366 
1367 end subroutine symspgr