TABLE OF CONTENTS


ABINIT/m_symtk [ Modules ]

[ Top ] [ Modules ]

NAME

  m_symtk

FUNCTION

  Low-level tools related to symmetries

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (RC, XG, GMR, MG, JWZ)
  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_symtk
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32 
33  use m_numeric_tools,  only : isinteger, wrap2_pmhalf
34  use m_hide_lapack,    only : matrginv
35 
36  implicit none
37 
38  private

m_symtk/chkgrp [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 chkgrp

FUNCTION

 Checks that a set of input symmetries constitutes a group.

INPUTS

 nsym = number of symmetry operations
 symafm = (anti)ferromagnetic part of symmetry operations
 symrel = 3D matrix containg symmetry operations

OUTPUT

  ierr=Status error.

TODO

 SHOULD ALSO CHECK THE tnons !

PARENTS

      chkinp,gensymspgr,m_bz_mesh,m_esymm,m_sigmaph,setsym

CHILDREN

SOURCE

365 subroutine chkgrp(nsym,symafm,symrel,ierr)
366 
367 
368 !This section has been created automatically by the script Abilint (TD).
369 !Do not modify the following lines by hand.
370 #undef ABI_FUNC
371 #define ABI_FUNC 'chkgrp'
372 !End of the abilint section
373 
374  implicit none
375 
376 !Arguments ------------------------------------
377 !scalars
378  integer,intent(in) :: nsym
379  integer,intent(out) :: ierr
380 !arrays
381  integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
382 
383 !Local variables-------------------------------
384 !scalars
385  integer :: isym,jsym,ksym,symafmchk,testeq=1
386  logical :: found_inv
387  character(len=500) :: msg
388 !arrays
389  integer :: chk(3,3)
390 
391 ! *************************************************************************
392 
393 !DEBUG
394 !write(std_out,*)' chkgrp : enter'
395 !write(std_out,*)'     isym         symrel            symafm '
396 !do isym=1,nsym
397 !write(std_out,'(i3,a,9i3,a,i3)' )isym,'   ',symrel(:,:,isym),'   ',symafm(isym)
398 !end do
399 !ENDDEBUG
400 
401  ierr = 0
402 
403 !1) Identity must be the first symmetry.
404  if (ANY(symrel(:,:,1) /= identity_3d .or. symafm(1)/=1 )) then
405    MSG_WARNING("First operation must be the identity operator")
406    ierr = ierr+1
407  end if
408 !
409 !2) The inverse of each element must belong to the group.
410  do isym=1,nsym
411    call mati3inv(symrel(:,:,isym),chk)
412    chk = TRANSPOSE(chk)
413    found_inv = .FALSE.
414    do jsym=1,nsym
415      if ( ALL(symrel(:,:,jsym) == chk) .and. (symafm(jsym)*symafm(isym) == 1 )) then
416        found_inv = .TRUE.; EXIT
417      end if
418    end do
419 
420    if (.not.found_inv) then
421      write(msg,'(a,i0,2a)')&
422 &     "Cannot find the inverse of symmetry operation ",isym,ch10,&
423 &     "Input symmetries do not form a group "
424      MSG_WARNING(msg)
425      ierr = ierr+1
426    end if
427 
428  end do
429 !
430 !Closure relation under composition.
431  do isym=1,nsym
432    do jsym=1,nsym
433 !
434 !    Compute the product of the two symmetries
435      chk = MATMUL(symrel(:,:,jsym), symrel(:,:,isym))
436      symafmchk=symafm(jsym)*symafm(isym)
437 !
438 !    Check that product array is one of the original symmetries.
439      do ksym=1,nsym
440        testeq=1
441        if ( ANY(chk/=symrel(:,:,ksym) )) testeq=0
442 #if 0
443 !      FIXME this check make v4/t26 and v4/t27 fails.
444 !      The rotational part is in the group but with different magnetic part!
445        if (symafmchk/=symafm(ksym))testeq=0
446 #endif
447        if (testeq==1) exit ! The test is positive
448      end do
449 !
450      if(testeq==0) then ! The test is negative
451        write(msg, '(a,2i3,a,7a)' )&
452 &       'product of symmetries',isym,jsym,' is not in group.',ch10,&
453 &       'This indicates that the input symmetry elements',ch10,&
454 &       'do not possess closure under group composition.',ch10,&
455 &       'Action: check symrel, symafm and fix them.'
456        MSG_WARNING(msg)
457        ierr = ierr+1
458      end if
459 
460    end do ! jsym
461  end do ! isym
462 
463 end subroutine chkgrp

m_symtk/chkorthsy [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 chkorthsy

FUNCTION

 Check the orthogonality of the symmetry operations
 (lengths and absolute values of scalar products should be preserved)

INPUTS

 gprimd(3,3)=dimensional primitive transl. for reciprocal space (bohr**-1)
 rmet=Real space metric.
 nsym=actual number of symmetries
 rprimd(3,3)=dimensional primitive translations for real space (bohr)
 symrel(3,3,1:nsym)=symmetry operations in real space in terms of primitive translations

SIDE EFFECTS

 iexit= if 0 at input, will do the check, and stop if there is a problem, return 0 if no problem
        if 1 at input, will always input, return 0 if no problem, -1 if there is a problem,
                       also, suppresses printing of problem

PARENTS

      chkinp,ingeo,symmetrize_rprimd

CHILDREN

SOURCE

640 subroutine chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel)
641 
642 
643 !This section has been created automatically by the script Abilint (TD).
644 !Do not modify the following lines by hand.
645 #undef ABI_FUNC
646 #define ABI_FUNC 'chkorthsy'
647 !End of the abilint section
648 
649  implicit none
650 
651 !Arguments ------------------------------------
652 !scalars
653  integer,intent(in) :: nsym
654  integer,intent(inout) :: iexit
655 !arrays
656  integer,intent(in) :: symrel(3,3,nsym)
657  real(dp),intent(in) :: gprimd(3,3),rmet(3,3),rprimd(3,3)
658 
659 !Local variables-------------------------------
660 !scalars
661  integer :: ii,isym,jj
662  real(dp),parameter :: tol=2.0d-8
663  real(dp) :: residual,rmet2
664  character(len=500) :: message
665 !arrays
666  real(dp) :: prods(3,3),rmet_sym(3,3),rprimd_sym(3,3)
667 
668 ! *************************************************************************
669 
670  rmet2=zero
671  do ii=1,3
672    do jj=1,3
673      rmet2=rmet2+rmet(ii,jj)*2
674    end do
675  end do
676 
677 !Loop over all symmetry operations
678  do isym=1,nsym
679 
680 !  Compute symmetric of primitive vectors under point symmetry operations
681    do ii=1,3
682      rprimd_sym(:,ii)=symrel(1,ii,isym)*rprimd(:,1)+&
683 &     symrel(2,ii,isym)*rprimd(:,2)+&
684 &     symrel(3,ii,isym)*rprimd(:,3)
685    end do
686 
687 !  If the new lattice is the same as the original one,
688 !  the lengths and angles are preserved
689    do ii=1,3
690      rmet_sym(ii,:)=rprimd_sym(1,ii)*rprimd_sym(1,:)+&
691 &     rprimd_sym(2,ii)*rprimd_sym(2,:)+&
692 &     rprimd_sym(3,ii)*rprimd_sym(3,:)
693    end do
694 
695    residual=zero
696    do ii=1,3
697      do jj=1,3
698        residual=residual+(rmet_sym(ii,jj)-rmet(ii,jj))**2
699      end do
700    end do
701 
702    if(sqrt(residual)>tol*sqrt(rmet2))then
703      if(iexit==0)then
704        write(message, '(a,i5,a,a,a,a,a,es12.4,a,a,a,a,a,a,a)' )&
705 &       'The symmetry operation number ',isym,' does not preserve',ch10,&
706 &       'vector lengths and angles.',ch10,&
707 &       'The value of the residual is ',residual,'.',ch10,&
708 &       'Action: modify rprim, acell and/or symrel so that',ch10,&
709 &       'vector lengths and angles are preserved.',ch10,&
710 &       'Beware, the tolerance on symmetry operations is very small.'
711        MSG_ERROR(message)
712      else
713        iexit=-1
714      end if
715    end if
716 
717 !  Also, the scalar product of rprimd_sym and gprimd must give integer numbers
718    do ii=1,3
719      prods(ii,:)=rprimd_sym(1,ii)*gprimd(1,:)+ &
720 &     rprimd_sym(2,ii)*gprimd(2,:)+ &
721 &     rprimd_sym(3,ii)*gprimd(3,:)
722    end do
723 
724    do ii=1,3
725      do jj=1,3
726        residual=prods(ii,jj)-anint(prods(ii,jj))
727        if(abs(residual)>tol)then
728          if(iexit==0)then
729            write(message, '(a,i0,a,a,a,a,a,a,a)' )&
730 &           'The symmetry operation number ',isym,' generates',ch10,&
731 &           'a different lattice.',ch10,&
732 &           'Action: modify rprim, acell and/or symrel so that',ch10,&
733 &           'the lattice is preserved.'
734            MSG_ERROR(message)
735          else
736            iexit=-1
737          end if
738        end if
739      end do
740    end do
741 
742    if(iexit==-1) exit
743  end do ! isym
744 
745  if(iexit==1)iexit=0
746 
747 end subroutine chkorthsy

m_symtk/chkprimit [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 chkprimit

FUNCTION

 Check whether the cell is primitive or not.
 If chkprim/=0 and the cell is non-primitive, stops.

INPUTS

 chkprim= if non-zero, check that the unit cell is primitive.
 nsym=actual number of symmetries
 symafm(nsym)= (anti)ferromagnetic part of symmetry operations
 symrel(3,3,nsym)= nsym symmetry operations in real space in terms
   of primitive translations

OUTPUT

  multi=multiplicity of the unit cell

PARENTS

      symanal

CHILDREN

SOURCE

775 subroutine chkprimit(chkprim, multi, nsym, symafm, symrel)
776 
777 
778 !This section has been created automatically by the script Abilint (TD).
779 !Do not modify the following lines by hand.
780 #undef ABI_FUNC
781 #define ABI_FUNC 'chkprimit'
782 !End of the abilint section
783 
784  implicit none
785 
786 !Arguments ------------------------------------
787 !scalars
788  integer,intent(in) :: chkprim,nsym
789  integer,intent(out) :: multi
790 !arrays
791  integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
792 
793 !Local variables-------------------------------
794 !scalars
795  integer :: isym
796  character(len=500) :: message
797 
798 !**************************************************************************
799 
800 !Loop over each symmetry operation of the Bravais lattice
801 !Find whether it is the identity, or a pure translation,
802 !without change of sign of the spin
803  multi=0
804  do isym=1,nsym
805    if( abs(symrel(1,1,isym)-1)+&
806 &   abs(symrel(2,2,isym)-1)+&
807 &   abs(symrel(3,3,isym)-1)+&
808 &   abs(symrel(1,2,isym))+abs(symrel(2,1,isym))+&
809 &   abs(symrel(2,3,isym))+abs(symrel(3,2,isym))+&
810 &   abs(symrel(3,1,isym))+abs(symrel(1,3,isym))+&
811 &   abs(symafm(isym)-1) == 0 )then
812      multi=multi+1
813    end if
814  end do
815 
816 !Check whether the cell is primitive
817  if(multi>1)then
818    if(chkprim/=0)then
819      write(message,'(a,a,a,i0,a,a,a,a,a,a,a,a,a)')&
820 &     'According to the symmetry finder, the unit cell is',ch10,&
821 &     'NOT primitive. The multiplicity is ',multi,' .',ch10,&
822 &     'The use of non-primitive unit cells is allowed',ch10,&
823 &     'only when the input variable chkprim is 0.',ch10,&
824 &     'Action: either change your unit cell (rprim or angdeg),',ch10,&
825 &     'or set chkprim to 0.'
826      MSG_ERROR(message)
827    else
828      write(message,'(3a,i0,a,a,a)')&
829 &     'According to the symmetry finder, the unit cell is',ch10,&
830 &     'not primitive, with multiplicity= ',multi,'.',ch10,&
831 &     'This is allowed, as the input variable chkprim is 0.'
832      MSG_COMMENT(message)
833    end if
834  end if
835 
836 end subroutine chkprimit

m_symtk/holocell [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 holocell

FUNCTION

 Examine whether the trial conventional cell described by cell_base
 is coherent with the required holohedral group.
 Possibly enforce the holohedry and modify the basis vectors.
 Note: for iholohedry=4, the tetragonal axis is not required to be along the C axis.

INPUTS

  enforce= if 0, only check; if =1, enforce exactly the holohedry
  iholohedry=required holohegral group
  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
  tolsym=tolerance for the symmetry operations

OUTPUT

  foundc=1 if the basis vectors supports the required holohedry ; =0 otherwise

SIDE EFFECTS

  cell_base(3,3)=basis vectors of the conventional cell  (changed if enforce==1, otherwise unchanged)

PARENTS

      symlatt,symmetrize_rprimd

CHILDREN

SOURCE

1265 subroutine holocell(cell_base,enforce,foundc,iholohedry,tolsym)
1266 
1267 
1268 !This section has been created automatically by the script Abilint (TD).
1269 !Do not modify the following lines by hand.
1270 #undef ABI_FUNC
1271 #define ABI_FUNC 'holocell'
1272 !End of the abilint section
1273 
1274  implicit none
1275 
1276 !Arguments ------------------------------------
1277 !scalars
1278  integer,intent(in) :: enforce,iholohedry
1279  integer,intent(out) :: foundc
1280  real(dp),intent(in) :: tolsym
1281 !arrays
1282  real(dp),intent(inout) :: cell_base(3,3)
1283 
1284 !Local variables ------------------------------
1285 !scalars
1286  integer :: allequal,ii,jj,orth
1287  real(dp):: aa,reldiff,scprod1
1288  character(len=500) :: msg
1289 !arrays
1290  integer :: ang90(3),equal(3)
1291  real(dp) :: length(3),metric(3,3),norm(3),rbasis(3,3),rconv(3,3),rconv_new(3,3)
1292  real(dp) :: rnormalized(3,3),symmetrized_length(3)
1293 
1294 !**************************************************************************
1295 
1296  do ii=1,3
1297    metric(:,ii)=cell_base(1,:)*cell_base(1,ii)+&
1298 &   cell_base(2,:)*cell_base(2,ii)+&
1299 &   cell_base(3,:)*cell_base(3,ii)
1300  end do
1301 
1302 !Examine the angles and vector lengths
1303  ang90(:)=0
1304  if(metric(1,2)**2<tolsym**2*metric(1,1)*metric(2,2))ang90(3)=1
1305  if(metric(1,3)**2<tolsym**2*metric(1,1)*metric(3,3))ang90(2)=1
1306  if(metric(2,3)**2<tolsym**2*metric(2,2)*metric(3,3))ang90(1)=1
1307  orth=0
1308  if(ang90(1)==1 .and. ang90(2)==1 .and. ang90(3)==1) orth=1
1309  equal(:)=0
1310  if(abs(metric(1,1)-metric(2,2))<tolsym*half*(metric(1,1)+metric(2,2)))equal(3)=1
1311  if(abs(metric(1,1)-metric(3,3))<tolsym*half*(metric(1,1)+metric(3,3)))equal(2)=1
1312  if(abs(metric(2,2)-metric(3,3))<tolsym*half*(metric(2,2)+metric(3,3)))equal(1)=1
1313  allequal=0
1314  if(equal(1)==1 .and. equal(2)==1 .and. equal(3)==1) allequal=1
1315 
1316  foundc=0
1317  if(iholohedry==1)                                      foundc=1
1318  if(iholohedry==2 .and. ang90(1)+ang90(3)==2 )          foundc=1
1319  if(iholohedry==3 .and. orth==1)                        foundc=1
1320  if(iholohedry==4 .and. orth==1 .and.          &
1321 & (equal(3)==1 .or. equal(2)==1 .or. equal(1)==1) ) foundc=1
1322  if(iholohedry==5 .and. allequal==1 .and. &
1323 & (abs(metric(1,2)-metric(2,3))<tolsym*metric(2,2)) .and. &
1324 & (abs(metric(1,2)-metric(1,3))<tolsym*metric(1,1))         )      foundc=1
1325  if(iholohedry==6 .and. equal(3)==1 .and. &
1326 & ang90(1)==1 .and. ang90(2)==1 .and. &
1327 & (2*metric(1,2)-metric(1,1))<tolsym*metric(1,1) )      foundc=1
1328  if(iholohedry==7 .and. orth==1 .and. allequal==1)      foundc=1
1329 
1330 !-------------------------------------------------------------------------------------
1331 !Possibly enforce the holohedry (if it is to be enforced !)
1332 
1333  if(foundc==1.and.enforce==1.and.iholohedry/=1)then
1334 
1335 !  Copy the cell_base vectors, and possibly fix the tetragonal axis to be the c-axis
1336    if(iholohedry==4.and.equal(1)==1)then
1337      rconv(:,3)=cell_base(:,1) ; rconv(:,1)=cell_base(:,2) ; rconv(:,2)=cell_base(:,3)
1338    else if (iholohedry==4.and.equal(2)==1)then
1339      rconv(:,3)=cell_base(:,2) ; rconv(:,2)=cell_base(:,1) ; rconv(:,1)=cell_base(:,3)
1340    else
1341      rconv(:,:)=cell_base(:,:)
1342    end if
1343 
1344 !  Compute the length of the three conventional vectors
1345    length(1)=sqrt(sum(rconv(:,1)**2))
1346    length(2)=sqrt(sum(rconv(:,2)**2))
1347    length(3)=sqrt(sum(rconv(:,3)**2))
1348 
1349 !  Take care of the first conventional vector aligned with rbasis(:,3) (or aligned with the trigonal axis if rhombohedral)
1350 !  and choice of the first normalized direction
1351    if(iholohedry==5)then
1352      rbasis(:,3)=third*(rconv(:,1)+rconv(:,2)+rconv(:,3))
1353    else
1354      rbasis(:,3)=rconv(:,3)
1355    end if
1356    norm(3)=sqrt(sum(rbasis(:,3)**2))
1357    rnormalized(:,3)=rbasis(:,3)/norm(3)
1358 
1359 !  Projection of the first conventional vector perpendicular to rbasis(:,3)
1360 !  and choice of the first normalized direction
1361    scprod1=sum(rnormalized(:,3)*cell_base(:,1))
1362    rbasis(:,1)=rconv(:,1)-rnormalized(:,3)*scprod1
1363    norm(1)=sqrt(sum(rbasis(:,1)**2))
1364    rnormalized(:,1)=rbasis(:,1)/norm(1)
1365 
1366 !  Generation of the second vector, perpendicular to the third and first
1367    rnormalized(1,2)=rnormalized(2,3)*rnormalized(3,1)-rnormalized(3,3)*rnormalized(2,1)
1368    rnormalized(2,2)=rnormalized(3,3)*rnormalized(1,1)-rnormalized(1,3)*rnormalized(3,1)
1369    rnormalized(3,2)=rnormalized(1,3)*rnormalized(2,1)-rnormalized(2,3)*rnormalized(1,1)
1370 
1371 !  Compute the vectors of the conventional cell, on the basis of iholohedry
1372    if(iholohedry==2)then
1373      rconv_new(:,3)=rconv(:,3)
1374      rconv_new(:,1)=rconv(:,1)
1375      rconv_new(:,2)=rnormalized(:,2)*length(2) ! Now, the y axis is perpendicular to the two others, that have not been changed
1376    else if(iholohedry==3.or.iholohedry==4.or.iholohedry==7)then
1377      if(iholohedry==7)then
1378        symmetrized_length(1:3)=sum(length(:))*third
1379      else if(iholohedry==4)then
1380        symmetrized_length(3)=length(3)
1381        symmetrized_length(1:2)=half*(length(1)+length(2))
1382      else if(iholohedry==3)then
1383        symmetrized_length(:)=length(:)
1384      end if
1385      do ii=1,3
1386        rconv_new(:,ii)=rnormalized(:,ii)*symmetrized_length(ii)
1387      end do
1388    else if(iholohedry==5)then
1389 !    In the normalized basis, they have coordinates (a,0,c), and (-a/2,+-sqrt(3)/2*a,c)
1390 !    c is known, but a is computed from the knowledge of the average length of the initial vectors
1391      aa=sqrt(sum(length(:)**2)*third-norm(3)**2)
1392      rconv_new(:,1)=aa*rnormalized(:,1)+rbasis(:,3)
1393      rconv_new(:,2)=aa*half*(-rnormalized(:,1)+sqrt(three)*rnormalized(:,2))+rbasis(:,3)
1394      rconv_new(:,3)=aa*half*(-rnormalized(:,1)-sqrt(three)*rnormalized(:,2))+rbasis(:,3)
1395    else if(iholohedry==6)then
1396 
1397 !    In the normalized basis, they have coordinates (a,0,0), (-a/2,+-sqrt(3)/2*a,0), and (0,0,c)
1398 !    c is known, but a is computed from the knowledge of the average length of the initial vectors
1399      aa=half*(length(1)+length(2))
1400      rconv_new(:,1)=aa*rnormalized(:,1)
1401      rconv_new(:,2)=aa*half*(-rnormalized(:,1)+sqrt(three)*rnormalized(:,2))
1402      rconv_new(:,3)=rconv(:,3)
1403    end if
1404 
1405 !  Check whether the modification make sense
1406    do ii=1,3
1407      do jj=1,3
1408        reldiff=(rconv_new(ii,jj)-rconv(ii,jj))/length(jj)
1409 !      Allow for twice tolsym
1410        if(abs(reldiff)>two*tolsym)then
1411          write(msg,'(a,6(2a,3es14.6))')&
1412 &         'Failed rectification of lattice vectors to comply with Bravais lattice identification, modifs are too large',ch10,&
1413 &         '  rconv    =',rconv(:,1),ch10,&
1414 &         '            ',rconv(:,2),ch10,&
1415 &         '            ',rconv(:,3),ch10,&
1416 &         '  rconv_new=',rconv_new(:,1),ch10,&
1417 &         '            ',rconv_new(:,2),ch10,&
1418 &         '            ',rconv_new(:,3)
1419          MSG_ERROR_CLASS(msg, "TolSymError")
1420        end if
1421      end do
1422    end do
1423 
1424 !  Copy back the cell_base vectors
1425    if(iholohedry==4.and.equal(1)==1)then
1426      cell_base(:,3)=rconv_new(:,2) ; cell_base(:,2)=rconv_new(:,1) ; cell_base(:,1)=rconv_new(:,3)
1427    else if (iholohedry==4.and.equal(2)==1)then
1428      cell_base(:,3)=rconv_new(:,1) ; cell_base(:,1)=rconv_new(:,2) ; cell_base(:,2)=rconv_new(:,3)
1429    else
1430      cell_base(:,:)=rconv_new(:,:)
1431    end if
1432 
1433  end if
1434 
1435 end subroutine holocell

m_symtk/littlegroup_q [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 littlegroup_q

FUNCTION

 Determines the symmetry operations by which reciprocal vector q is preserved,
 modulo a primitive reciprocal lattice vector, and the time-reversal symmetry.

INPUTS

 nsym=number of space group symmetries
 qpt(3)= vector in reciprocal space
 symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
 [prtvol]=integer flag defining the verbosity of output. =0 if no output is provided.
 prtgkk= integer flag. If 1 provide output of electron-phonon "gkk" matrix elements, for further
     treatment by mrggkk utility or anaddb utility. If 0 no output is provided.

OUTPUT

 symq(4,2,nsym)= three first numbers define the G vector;
     fourth number is zero if the q-vector is not preserved, is 1 otherwise
     second index is one without time-reversal symmetry, two with time-reversal symmetry
 timrev=1 if the time-reversal symmetry preserves the wavevector,
   modulo a reciprocal lattice vector (in principle, see below).

NOTES

 The condition is:

    $q =  O  S(q) - G$

 with O being either the identity or the time reversal symmetry (= inversion in reciprocal space)
 and G being a primitive vector of the reciprocal lattice.
 If the time-reversal (alone) also preserves q, modulo a lattice vector, then timrev is set to 1, otherwise 0.

TODO

 timrev is put to 1 only for Gamma.
 Better handling should be provided in further version.

PARENTS

      get_npert_rbz,m_bz_mesh,m_ddb,m_dvdb,m_dynmat,m_esymm,m_gkk,m_phgamma
      m_sigmaph,memory_eval,read_gkk,respfn

CHILDREN

      wrap2_pmhalf,wrtout

SOURCE

 987 subroutine littlegroup_q(nsym,qpt,symq,symrec,symafm,timrev,prtvol,use_sym)
 988 
 989 
 990 !This section has been created automatically by the script Abilint (TD).
 991 !Do not modify the following lines by hand.
 992 #undef ABI_FUNC
 993 #define ABI_FUNC 'littlegroup_q'
 994 !End of the abilint section
 995 
 996  implicit none
 997 
 998 !Arguments -------------------------------
 999 !scalars
1000  integer,intent(in) :: nsym
1001  integer,intent(in),optional :: prtvol,use_sym
1002  integer,intent(out) :: timrev
1003 !arrays
1004  integer,intent(in) :: symrec(3,3,nsym)
1005  integer,intent(in) :: symafm(nsym)
1006  integer,intent(out) :: symq(4,2,nsym)
1007  real(dp),intent(in) :: qpt(3)
1008 
1009 !Local variables -------------------------
1010 !scalars
1011  integer :: ii,isign,isym,itirev,my_prtvol
1012  real(dp),parameter :: tol=2.d-8
1013  real(dp) :: reduce
1014  character(len=500) :: message
1015 !arrays
1016  real(dp) :: difq(3),qsym(3),shift(3)
1017 
1018 ! *********************************************************************
1019 
1020  my_prtvol=0 ; if (PRESENT(prtvol)) my_prtvol=prtvol
1021 
1022 ! Initialise the array symq
1023  symq = 0
1024 
1025  isym = symafm(1) ! just to fool abirules and use symafm for the moment
1026 
1027  do isym=1,nsym
1028 !   if (symafm(isym) /= 1) cycle ! skip afm symops
1029 ! TODO: check how much of the afm syms are coded in the rf part of the code. cf
1030 ! test v3 / 12
1031    do itirev=1,2
1032      isign=3-2*itirev  ! isign is 1 without time-reversal, -1 with time-reversal
1033 
1034 !    Get the symmetric of the vector
1035      do ii=1,3
1036        qsym(ii)=qpt(1)*isign*symrec(ii,1,isym)&
1037 &       +qpt(2)*isign*symrec(ii,2,isym)&
1038 &       +qpt(3)*isign*symrec(ii,3,isym)
1039      end do
1040 
1041 !    Get the difference between the symmetric and the original vector
1042 
1043      symq(4,itirev,isym)=1
1044      do ii=1,3
1045        difq(ii)=qsym(ii)-qpt(ii)
1046 !      Project modulo 1 in the interval ]-1/2,1/2] such that difq = reduce + shift
1047        call wrap2_pmhalf(difq(ii),reduce,shift(ii))
1048        if(abs(reduce)>tol)symq(4,itirev,isym)=0
1049      end do
1050 
1051 !    SP: When prtgkk is asked (GKK matrix element will be output), one has to
1052 !    disable symmetries. There is otherwise a jauge problem with the unperturbed
1053 !    and the perturbed wavefunctions. This leads to a +- 5% increase in computational
1054 !    cost but provide the correct GKKs (i.e. the same as without the use of
1055 !    symmerties.)
1056 
1057      if (PRESENT(use_sym)) then
1058        if (use_sym == 0) then
1059          symq(4,itirev,isym)=0
1060          symq(4,itirev,1)=1
1061        end if
1062      end if
1063 
1064 !    If the operation succeded, change shift from real(dp) to integer, then exit loop
1065      if(symq(4,itirev,isym)/=0)then
1066        if (my_prtvol>0) then
1067          if(itirev==1)write(message,'(a,i4,a)')' littlegroup_q : found symmetry',isym,' preserves q '
1068          if(itirev==2)write(message,'(a,i4,a)')' littlegroup_q : found symmetry ',isym,' + TimeReversal preserves q '
1069          call wrtout(std_out,message,'COLL')
1070        end if
1071 !      Uses the mathematical function NINT = nearest integer
1072        do ii=1,3
1073          symq(ii,itirev,isym)=nint(shift(ii))
1074        end do
1075      end if
1076 
1077    end do !itirev
1078  end do !isym
1079 
1080 !Test time-reversal symmetry
1081  timrev=1
1082  do ii=1,3
1083 !  Unfortunately, this version does not work yet ...
1084 !  call wrap2_pmhalf(2*qpt(ii),reduce,shift(ii))
1085 !  if(abs(reduce)>tol)timrev=0
1086 !  So, this is left ...
1087    if(abs(qpt(ii))>tol)timrev=0
1088  end do
1089 
1090  if(timrev==1.and.my_prtvol>0)then
1091    write(message, '(a,a,a)' )&
1092 &   ' littlegroup_q : able to use time-reversal symmetry. ',ch10,&
1093 &   '  (except for gamma, not yet able to use time-reversal symmetry)'
1094    call wrtout(std_out,message,'COLL')
1095  end if
1096 
1097 end subroutine littlegroup_q

m_symtk/mati3det [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 mati3det

FUNCTION

 Compute the determinant of a 3x3 matrix of INTEGER elements.

INPUTS

 mm = integer matrix

OUTPUT

 det = determinant of the matrix

PARENTS

      get_kpt_fullbz,getspinrot,m_ab7_symmetry,m_anaddb_dataset,symdet

CHILDREN

SOURCE

172 subroutine mati3det(mm, det)
173 
174 
175 !This section has been created automatically by the script Abilint (TD).
176 !Do not modify the following lines by hand.
177 #undef ABI_FUNC
178 #define ABI_FUNC 'mati3det'
179 !End of the abilint section
180 
181  implicit none
182 
183 !Arguments ------------------------------------
184 !arrays
185  integer,intent(in) :: mm(3,3)
186  integer,intent(out) :: det
187 
188 ! *************************************************************************
189  det=mm(1,1)*(mm(2,2) * mm(3,3) - mm(3,2) * mm(2,3)) &
190 & + mm(2,1)*(mm(3,2) * mm(1,3) - mm(1,2) * mm(3,3)) &
191 & + mm(3,1)*(mm(1,2) * mm(2,3) - mm(2,2) * mm(1,3))
192 
193 end subroutine mati3det

m_symtk/mati3inv [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 mati3inv

FUNCTION

 Invert and transpose orthogonal 3x3 matrix of INTEGER elements.

INPUTS

 mm = integer matrix to be inverted

OUTPUT

 mit = inverse of mm input matrix

NOTES

 Used for symmetry operations.
 This routine applies to ORTHOGONAL matrices only.
 Since these form a group, inverses are also integer
 arrays.  Returned array is TRANSPOSE of inverse, as needed.
 Note use of integer arithmetic.

PARENTS

      cg_rotate,chkgrp,classify_bands,debug_tools,dfpt_nstdy,get_full_kgrid
      get_npert_rbz,getkgrid,ingeo,m_ab7_symmetry,m_crystal,m_ddb,m_dvdb
      m_dynmat,m_fft_mesh,m_fock,m_ptgroups,m_tdep_sym,matpointsym
      memory_eval,optic,read_gkk,setsym,strainsym,thmeig,wfconv

CHILDREN

SOURCE

 97 subroutine mati3inv(mm, mit)
 98 
 99 
100 !This section has been created automatically by the script Abilint (TD).
101 !Do not modify the following lines by hand.
102 #undef ABI_FUNC
103 #define ABI_FUNC 'mati3inv'
104 !End of the abilint section
105 
106  implicit none
107 
108 !Arguments ------------------------------------
109 !arrays
110  integer,intent(in) :: mm(3,3)
111  integer,intent(out) :: mit(3,3)
112 
113 !Local variables-------------------------------
114 !scalars
115  integer :: dd
116  character(len=500) :: message
117 !arrays
118  integer :: tt(3,3)
119 
120 ! *************************************************************************
121 
122  tt(1,1) = mm(2,2) * mm(3,3) - mm(3,2) * mm(2,3)
123  tt(2,1) = mm(3,2) * mm(1,3) - mm(1,2) * mm(3,3)
124  tt(3,1) = mm(1,2) * mm(2,3) - mm(2,2) * mm(1,3)
125  tt(1,2) = mm(3,1) * mm(2,3) - mm(2,1) * mm(3,3)
126  tt(2,2) = mm(1,1) * mm(3,3) - mm(3,1) * mm(1,3)
127  tt(3,2) = mm(2,1) * mm(1,3) - mm(1,1) * mm(2,3)
128  tt(1,3) = mm(2,1) * mm(3,2) - mm(3,1) * mm(2,2)
129  tt(2,3) = mm(3,1) * mm(1,2) - mm(1,1) * mm(3,2)
130  tt(3,3) = mm(1,1) * mm(2,2) - mm(2,1) * mm(1,2)
131  dd  = mm(1,1) * tt(1,1) + mm(2,1) * tt(2,1) + mm(3,1) * tt(3,1)
132 
133 !Make sure matrix is not singular
134  if (dd/=0) then
135    mit(:,:)=tt(:,:)/dd
136  else
137    write(message, '(2a,2x,9i5,a)' )&
138 &   'Attempting to invert integer array',ch10,mm(:,:),'   ==> determinant is zero.'
139    MSG_BUG(message)
140  end if
141 
142 !If matrix is orthogonal, determinant must be 1 or -1
143  if (abs(dd)/=1) then
144    write(message, '(2a,2x,9i5,a)' )&
145 &   'Absolute value of determinant should be one',ch10,'but determinant= ',dd
146    MSG_BUG(message)
147  end if
148 
149 end subroutine mati3inv

m_symtk/matpointsym [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 matpointsym

FUNCTION

 For given order of point group, symmetrizes a 3x3 input matrix using the
 point symmetry of the input atom

INPUTS

 iatom=index of atom to symmetrize around
 natom=number of atoms in cell
 nsym=order of group
 rprimd(3,3)= real space primitive vectors
 symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
 tnons(3,nsym) = nonsymmorphic translations
 xred(3,natom)=locations of atoms in reduced coordinates

OUTPUT

SIDE EFFECTS

 mat3(3,3) = matrix to be symmetrized, in cartesian frame

PARENTS

      make_efg_el,make_efg_ion,make_efg_onsite

CHILDREN

      dgemm,dgemv,mati3inv,matrginv

SOURCE

1130 subroutine matpointsym(iatom,mat3,natom,nsym,rprimd,symrel,tnons,xred)
1131 
1132  use m_linalg_interfaces
1133 
1134 !This section has been created automatically by the script Abilint (TD).
1135 !Do not modify the following lines by hand.
1136 #undef ABI_FUNC
1137 #define ABI_FUNC 'matpointsym'
1138 !End of the abilint section
1139 
1140  implicit none
1141 
1142 !Arguments ------------------------------------
1143 !scalars
1144  integer,intent(in) :: iatom,natom,nsym
1145 !arrays
1146  integer,intent(in) :: symrel(3,3,nsym)
1147  real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym),xred(3,natom)
1148  real(dp),intent(inout) :: mat3(3,3)
1149 
1150 !Local variables-------------------------------
1151 !scalars
1152  integer :: cell_index,cell_indexp,ii,isym,nsym_point
1153  real(dp) :: xreddiff
1154 !arrays
1155  integer :: symrel_it(3,3)
1156  real(dp) :: mat3_tri(3,3),mat3_tri_sym(3,3),rprimd_inv(3,3),tmp_mat(3,3)
1157  real(dp) :: xredp(3)
1158 
1159 !**************************************************************************
1160 
1161 !copy rprimd input and construct inverse
1162  rprimd_inv = rprimd
1163  call matrginv(rprimd_inv,3,3)
1164 
1165 !transform input mat3 to triclinic frame with rprimd^{-1} * mat3 * rprimd
1166  call dgemm('N','N',3,3,3,one,rprimd_inv,3,mat3,3,zero,tmp_mat,3)
1167  call dgemm('N','N',3,3,3,one,tmp_mat,3,rprimd,3,zero,mat3_tri,3)
1168 
1169 !loop over symmetry elements to obtain symmetrized input matrix
1170  mat3_tri_sym = zero
1171  nsym_point = 0
1172  do isym = 1, nsym
1173 
1174 ! skip any nonsymmorphic symmetry elements, want to consider point elements only
1175    if(dot_product(tnons(:,isym),tnons(:,isym))>tol8) cycle
1176 
1177 ! for current symmetry element, find transformed reduced coordinates of target atom
1178 ! via xredp = symrel * xred
1179    call dgemv('N',3,3,one,dble(symrel(:,:,isym)),3,xred(:,iatom),1,zero,xredp,1)
1180 
1181 
1182 ! shift xredp into the same unit cell as xred, for comparison
1183 ! label cells as 0..1:0 1..2:1 2..3:2 and -1..0:-1 -2..-1:-2 and so forth
1184    do ii = 1, 3
1185 
1186      cell_index = int(xred(ii,iatom))
1187      if(xred(ii,iatom) < zero) cell_index = cell_index - 1
1188      cell_indexp = int(xredp(ii))
1189      if(xredp(ii) < zero) cell_indexp = cell_indexp - 1
1190 
1191      do while (cell_indexp < cell_index)
1192        xredp(ii) = xredp(ii)+one
1193        cell_indexp = cell_indexp + 1
1194      end do
1195      do while (cell_indexp > cell_index)
1196        xredp(ii) = xredp(ii)-one
1197        cell_indexp = cell_indexp - 1
1198      end do
1199 
1200    end do
1201 
1202 ! now compare xredp to xred
1203    xreddiff = dot_product(xredp-xred(:,iatom),xredp-xred(:,iatom))
1204 
1205    if (xreddiff < tol8) then
1206 
1207 !  accumulate symrel^{-1}*mat3_tri*symrel into mat3_tri_sym iff xredp = xred + L,
1208 !  where is a lattice vector, so symrel leaves the target atom invariant
1209 
1210 !  mati3inv gives the inverse transpose of symrel
1211      call mati3inv(symrel(:,:,isym),symrel_it)
1212      call dgemm('N','N',3,3,3,one,mat3_tri,3,dble(symrel(:,:,isym)),3,zero,tmp_mat,3)
1213      call dgemm('T','N',3,3,3,one,dble(symrel_it),3,tmp_mat,3,one,mat3_tri_sym,3)
1214      nsym_point = nsym_point + 1
1215    end if
1216 
1217  end do
1218 
1219 !normalize by number of point symmetry operations
1220  mat3_tri_sym = mat3_tri_sym/dble(nsym_point)
1221 
1222 !transform mat3_tri_sym to cartesian frame with rprimd * mat3_tri_sym * rprimd^{-1}
1223 
1224  call dgemm('N','N',3,3,3,one,mat3_tri_sym,3,rprimd_inv,3,zero,tmp_mat,3)
1225  call dgemm('N','N',3,3,3,one,rprimd,3,tmp_mat,3,zero,mat3,3)
1226 
1227 end subroutine matpointsym

m_symtk/matr3inv [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 matr3inv

FUNCTION

 Invert and transpose general 3x3 matrix of real*8 elements.

INPUTS

 aa = 3x3 matrix to be inverted

OUTPUT

 ait = inverse of aa input matrix

NOTES

 Returned array is TRANSPOSE of inverse, as needed to get g from r.

PARENTS

      berryphase,chkdilatmx,conducti_nc,ddb_hybrid,dfpt_mkvxc,dfpt_mkvxcstr
      dfpt_symph,electrooptic,ep_el_weights,ep_fs_weights,ep_ph_weights
      fock_getghc,get_kpt_fullbz,getkgrid,getspinrot,gstate,harmonic_thermo
      invars2,inwffil,m_cut3d,m_ddb,m_ddk,m_double_grid,m_dynmat
      m_effective_potential,m_esymm,m_ewald,m_fock,m_fstab,m_ifc,m_pimd
      m_psps,m_strain,m_supercell,m_tdep_latt,make_efg_el,make_efg_ion,metric
      mover,optic,outwant,pimd_langevin_npt,prtxf,relaxpol,respfn,smpbz
      stresssym,symbrav,symlatt,symmetrize_rprimd,symrelrot,symrhg,tddft
      testkgrid,thmeig,uderiv,xcart2xred,xfpack_x2vin

CHILDREN

SOURCE

227 subroutine matr3inv(aa, ait)
228 
229 
230 !This section has been created automatically by the script Abilint (TD).
231 !Do not modify the following lines by hand.
232 #undef ABI_FUNC
233 #define ABI_FUNC 'matr3inv'
234 !End of the abilint section
235 
236  implicit none
237 
238 !Arguments ------------------------------------
239 !arrays
240  real(dp),intent(in) :: aa(3,3)
241  real(dp),intent(out) :: ait(3,3)
242 
243 !Local variables-------------------------------
244 !scalars
245  real(dp) :: dd,det,t1,t2,t3
246  character(len=500) :: message
247 
248 ! *************************************************************************
249 
250  t1 = aa(2,2) * aa(3,3) - aa(3,2) * aa(2,3)
251  t2 = aa(3,2) * aa(1,3) - aa(1,2) * aa(3,3)
252  t3 = aa(1,2) * aa(2,3) - aa(2,2) * aa(1,3)
253  det  = aa(1,1) * t1 + aa(2,1) * t2 + aa(3,1) * t3
254 
255 !Make sure matrix is not singular
256  if (abs(det)>tol16) then
257    dd=one/det
258  else
259    write(message, '(2a,2x,9es16.8,a,a,es16.8,a)' )&
260 &   'Attempting to invert real(8) 3x3 array',ch10,aa(:,:),ch10,'   ==> determinant=',det,' is zero.'
261    MSG_BUG(message)
262  end if
263 
264  ait(1,1) = t1 * dd
265  ait(2,1) = t2 * dd
266  ait(3,1) = t3 * dd
267  ait(1,2) = (aa(3,1)*aa(2,3)-aa(2,1)*aa(3,3)) * dd
268  ait(2,2) = (aa(1,1)*aa(3,3)-aa(3,1)*aa(1,3)) * dd
269  ait(3,2) = (aa(2,1)*aa(1,3)-aa(1,1)*aa(2,3)) * dd
270  ait(1,3) = (aa(2,1)*aa(3,2)-aa(3,1)*aa(2,2)) * dd
271  ait(2,3) = (aa(3,1)*aa(1,2)-aa(1,1)*aa(3,2)) * dd
272  ait(3,3) = (aa(1,1)*aa(2,2)-aa(2,1)*aa(1,2)) * dd
273 
274 end subroutine matr3inv

m_symtk/print_symmetries [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 print_symmetries

FUNCTION

  Helper function to print the set of symmetries.

INPUTS

OUTPUT

PARENTS

      gensymspgr,hdr_vs_dtset,m_crystal

CHILDREN

      mati3inv,sg_multable

SOURCE

3136 subroutine print_symmetries(nsym, symrel, tnons, symafm, unit, mode_paral)
3137 
3138 
3139 !This section has been created automatically by the script Abilint (TD).
3140 !Do not modify the following lines by hand.
3141 #undef ABI_FUNC
3142 #define ABI_FUNC 'print_symmetries'
3143 !End of the abilint section
3144 
3145  implicit none
3146 
3147 !Arguments ------------------------------------
3148 !scalars
3149  integer,intent(in) :: nsym
3150  integer,optional,intent(in) :: unit
3151  character(len=4),optional,intent(in) :: mode_paral
3152 !arrays
3153  integer,intent(in) :: symrel(3,3,nsym),symafm(nsym)
3154  real(dp),intent(in) :: tnons(3,nsym)
3155 
3156 !Local variables-------------------------------
3157  integer :: my_unt,isym,isymin,isymend,ii,jj
3158  character(len=500) :: msg
3159  character(len=4) :: my_mode
3160 ! *********************************************************************
3161 
3162  my_unt =std_out; if (PRESENT(unit      )) my_unt =unit
3163  my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode=mode_paral
3164 
3165  !write(msg,'(2a)')ch10,' Rotations                           Translations     Symafm '
3166  !do isym=1,nsym
3167  ! write(msg,'(1x,3(3i3,1x),4x,3(f11.7,1x),6x,i2)')symrel(:,:,isym),tnons(:,isym),symafm(isym)
3168  ! call wrtout(my_unt,msg,my_mode)
3169  !end do
3170 
3171  write(msg,'(2a)')ch10,' Symmetry operations in real space (Rotation tnons AFM)'
3172  call wrtout(my_unt,msg,my_mode)
3173  do isymin=1,nsym,4
3174    isymend=isymin+3
3175    if (isymend>nsym) isymend=nsym
3176    do ii=1,3
3177      write(msg,'(4(3i3,f8.3,i3,3x))')((symrel(ii,jj,isym),jj=1,3),tnons(ii,isym),symafm(isym),isym=isymin,isymend)
3178      call wrtout(my_unt,msg,my_mode)
3179    end do
3180    write(msg,'(a)')ch10
3181    call wrtout(my_unt,msg,my_mode)
3182  end do
3183 
3184 end subroutine print_symmetries

m_symtk/sg_multable [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 sg_multable

FUNCTION

 Checks that a set of input symmetries constitutes a group.

INPUTS

 nsym=number of symmetry operations
 symafm(nsym)=(anti)ferromagnetic part of symmetry operations
 symrel(3,3,nsym)=symmetry operations in real space.
 tnons(3,nsym)=Fractional translations.

OUTPUT

  ierr=Status error. A non-zero value signals a failure.
  [multable(4,nsym,nsym)]= Optional output.
    multable(1,sym1,sym2) gives the index of the symmetry product S1 * S2 in the symrel array. 0 if not found.
    multable(2:4,sym1,sym2)= the lattice vector that has to added to the fractional translation
      of the operation of index multable(1,sym1,sym2) to obtain the fractional traslation of the product S1 * S2.
  [toinv(4,nsym)]= Optional output.
    toinv(1,sym1)=Gives the index of the inverse of the symmetry operation.
     S1 * S1^{-1} = {E, L} with E identity and L a lattice vector
    toinv(2:4,sym1)=The lattice vector L
      Note that toinv can be easily obtained from multable but sometimes we do not need the full table.

TODO

  This improved version should replace chkgrp.

PARENTS

      m_crystal,m_shirley

CHILDREN

SOURCE

501 subroutine sg_multable(nsym,symafm,symrel,tnons,tnons_tol,ierr,multable,toinv)
502 
503 
504 !This section has been created automatically by the script Abilint (TD).
505 !Do not modify the following lines by hand.
506 #undef ABI_FUNC
507 #define ABI_FUNC 'sg_multable'
508 !End of the abilint section
509 
510  implicit none
511 
512 !Arguments ------------------------------------
513 !scalars
514  integer,intent(in) :: nsym
515  integer,intent(out) :: ierr
516  real(dp),intent(in) :: tnons_tol
517 !arrays
518  integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
519  integer,optional,intent(out) :: multable(4,nsym,nsym)
520  integer,optional,intent(out) :: toinv(4,nsym)
521  real(dp),intent(in) :: tnons(3,nsym)
522 
523 !Local variables-------------------------------
524 !scalars
525  integer :: sym1,sym2,sym3,prd_symafm
526  logical :: found_inv,iseq
527  character(len=500) :: msg
528 !arrays
529  integer :: prd_symrel(3,3)
530  real(dp) :: prd_tnons(3)
531 
532 ! *************************************************************************
533 
534  ierr = 0
535 
536 !1) Identity must be the first symmetry. Do not check tnons, cell might not be primitive.
537  if (ANY(symrel(:,:,1) /= identity_3d .or. symafm(1)/=1 )) then
538    MSG_WARNING("First operation must be the identity operator")
539    ierr = ierr+1
540  end if
541 !
542 !2) The inverse of each element must belong to the group.
543  do sym1=1,nsym
544    found_inv = .FALSE.
545    do sym2=1,nsym
546      prd_symrel = MATMUL(symrel(:,:,sym1), symrel(:,:,sym2))
547      prd_tnons = tnons(:,sym1) + MATMUL(symrel(:,:,sym1),tnons(:,sym2))
548      prd_symafm = symafm(sym1)*symafm(sym2)
549      if ( ALL(prd_symrel == identity_3d) .and. isinteger(prd_tnons,tnons_tol) .and. prd_symafm == 1 ) then
550        found_inv = .TRUE.
551        if (PRESENT(toinv)) then
552          toinv(1,sym1) = sym2
553          toinv(2:4,sym1) = NINT(prd_tnons)
554        end if
555        EXIT
556      end if
557    end do
558 
559    if (.not.found_inv) then
560      write(msg,'(a,i0,2a)')&
561 &     "Cannot find the inverse of symmetry operation ",sym1,ch10,&
562 &     "Input symmetries do not form a group "
563      MSG_WARNING(msg)
564      ierr = ierr+1
565    end if
566  end do
567 !
568 !Check closure relation under composition and construct multiplication table.
569  do sym1=1,nsym
570    do sym2=1,nsym
571 !
572 !    Compute the product of the two symmetries. Convention {A,a} {B,b} = {AB, a+Ab}
573      prd_symrel = MATMUL(symrel(:,:,sym1), symrel(:,:,sym2))
574      prd_symafm = symafm(sym1)*symafm(sym2)
575      prd_tnons = tnons(:,sym1) + MATMUL(symrel(:,:,sym1),tnons(:,sym2))
576 !
577      iseq=.FALSE.
578      do sym3=1,nsym ! Check that product array is one of the original symmetries.
579        iseq = ( ALL(prd_symrel==symrel(:,:,sym3) )           .and. &
580 &       isinteger(prd_tnons-tnons(:,sym3),tnons_tol) .and. &
581 &       prd_symafm==symafm(sym3) )  ! Here v4/t26 and v4/t27 will fail.
582 !      The rotational part is in the group but with different magnetic part!
583 
584        if (iseq) then ! The test is positive
585          if (PRESENT(multable)) then
586            multable(1,sym1,sym2) = sym3
587            multable(2:4,sym1,sym2) = NINT(prd_tnons-tnons(:,sym3))
588          end if
589          EXIT
590        end if
591      end do
592 !
593      if (.not.iseq) then ! The test is negative
594        write(msg, '(a,2(i0,1x),a,7a)' )&
595 &       'Product of symmetries:',sym1,sym2,' is not in group.',ch10,&
596 &       'This indicates that the input symmetry elements',ch10,&
597 &       'do not possess closure under group composition.',ch10,&
598 &       'Action: check symrel, symafm and fix them.'
599        MSG_WARNING(msg)
600        ierr = ierr+1
601        if (PRESENT(multable)) then
602          multable(1,sym1,sym2) = 0
603          multable(2:4,sym1,sym2) = HUGE(0)
604        end if
605      end if
606 
607    end do ! sym2
608  end do ! sym1
609 
610 end subroutine sg_multable

m_symtk/smallprim [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 smallprim

FUNCTION

 Find the smallest possible primitive vectors for an input lattice
 This algorithm is not as restrictive as the conditions mentioned at p.740
 of the international tables for crystallography (1983).
 The final vectors form a right-handed basis, while their
 sign and ordering is chosen such as to maximize the overlap
 with the original vectors in order.

INPUTS

  rprimd(3,3)=primitive vectors

OUTPUT

  metmin(3,3)=metric for the new (minimal) primitive vectors
  minim(3,3)=minimal primitive translations

NOTES

 The routine might as well be defined without
 metmin as argument, but it is more convenient to have it

PARENTS

      getkgrid,symlatt,testkgrid

CHILDREN

      metric

SOURCE

2905 subroutine smallprim(metmin,minim,rprimd)
2906 
2907 
2908 !This section has been created automatically by the script Abilint (TD).
2909 !Do not modify the following lines by hand.
2910 #undef ABI_FUNC
2911 #define ABI_FUNC 'smallprim'
2912 !End of the abilint section
2913 
2914  implicit none
2915 
2916 !Arguments ------------------------------------
2917 !arrays
2918  real(dp),intent(in) :: rprimd(3,3)
2919  real(dp),intent(out) :: metmin(3,3),minim(3,3)
2920 
2921 !Local variables-------------------------------
2922 !scalars
2923  integer :: ia,ib,ii,itrial,minimal
2924  integer :: iiter, maxiter = 100000
2925  real(dp) :: determinant,length2,metsum
2926  character(len=500) :: message
2927 !arrays
2928  integer :: nvecta(3),nvectb(3)
2929  real(dp) :: rmet(3,3),scprod(3),tmpvect(3)
2930 
2931 !**************************************************************************
2932 
2933  !call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2934  rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
2935 
2936  nvecta(1)=2 ; nvectb(1)=3
2937  nvecta(2)=1 ; nvectb(2)=3
2938  nvecta(3)=1 ; nvectb(3)=2
2939 
2940  minim(:,:)=rprimd(:,:)
2941  metmin(:,:)=rmet(:,:)
2942 
2943 !DEBUG
2944 !write(std_out,*)' smallprim : starting values, rprim '
2945 !write(std_out,'(3f16.8)' )rprimd(:,1)
2946 !write(std_out,'(3f16.8)' )rprimd(:,2)
2947 !write(std_out,'(3f16.8)' )rprimd(:,3)
2948 !write(std_out,*)' smallprim : starting values, rmet '
2949 !write(std_out,'(3f16.8)' )rmet(:,1)
2950 !write(std_out,'(3f16.8)' )rmet(:,2)
2951 !write(std_out,'(3f16.8)' )rmet(:,3)
2952 !ENDDEBUG
2953 
2954 !Note this loop without index
2955  do iiter = 1, maxiter
2956 
2957 !  Will exit if minimal=1 is still valid after a trial
2958 !  to reduce the vectors of each of the three pairs
2959    minimal=1
2960 
2961    do itrial=1,3
2962 
2963      ia=nvecta(itrial) ; ib=nvectb(itrial)
2964 !    Make sure the scalar product is negative
2965      if(metmin(ia,ib)>tol8)then
2966        minim(:,ia)=-minim(:,ia)
2967        metmin(ia,ib)=-metmin(ia,ib) ; metmin(ib,ia)=-metmin(ib,ia)
2968        metmin(ia,itrial)=-metmin(ia,itrial)
2969        metmin(itrial,ia)=-metmin(itrial,ia)
2970      end if
2971 !    Compute the length of the sum vector
2972      length2=metmin(ia,ia)+2*metmin(ia,ib)+metmin(ib,ib)
2973 !    Replace the first vector by the sum vector if the latter is smaller
2974      if(length2/metmin(ia,ia) < one-tol8)then
2975        minim(:,ia)=minim(:,ia)+minim(:,ib)
2976        metmin(ia,ia)=length2
2977        metmin(ia,ib)=metmin(ia,ib)+metmin(ib,ib)
2978        metmin(ia,itrial)=metmin(ia,itrial)+metmin(ib,itrial)
2979        metmin(ib,ia)=metmin(ia,ib)
2980        metmin(itrial,ia)=metmin(ia,itrial)
2981        minimal=0
2982 !      Replace the second vector by the sum vector if the latter is smaller
2983      else if(length2/metmin(ib,ib) < one-tol8)then
2984        minim(:,ib)=minim(:,ia)+minim(:,ib)
2985        metmin(ib,ib)=length2
2986        metmin(ia,ib)=metmin(ia,ib)+metmin(ia,ia)
2987        metmin(itrial,ib)=metmin(itrial,ib)+metmin(itrial,ia)
2988        metmin(ib,ia)=metmin(ia,ib)
2989        metmin(ib,itrial)=metmin(itrial,ib)
2990        minimal=0
2991      end if
2992 
2993    end do
2994 
2995    if(minimal==1)exit
2996  end do
2997 
2998  if (iiter >= maxiter) then
2999    write(message,'(a,i0,a)') 'the loop has failed to find a set of minimal vectors in ',maxiter,' iterations.'
3000    MSG_BUG(message)
3001  end if
3002 
3003 !At this stage, the three vectors have angles between each other that are
3004 !comprised between 90 and 120 degrees. It might still be that minus the vector
3005 !that is the sum of the three vectors is smaller than the longest of these vectors
3006  do iiter = 1, maxiter
3007 
3008 !  Will exit if minimal=1 is still valid after a trial
3009 !  to replace each of the three vectors by minus the summ of the three vectors
3010    minimal=1
3011    metsum=sum(metmin(:,:))
3012    do itrial=1,3
3013      ia=nvecta(itrial) ; ib=nvectb(itrial)
3014      if(metmin(ia,ia)/metsum > one + tol8)then
3015        minim(:,ia)=-minim(:,1)-minim(:,2)-minim(:,3)
3016        metmin(ia,ib)=-sum(metmin(:,ib))
3017        metmin(ia,itrial)=-sum(metmin(:,itrial))
3018        metmin(ia,ia)=metsum
3019        metmin(ib,ia)=metmin(ia,ib)
3020        metmin(itrial,ia)=metmin(ia,itrial)
3021        minimal=0
3022      end if
3023    end do
3024 
3025    if(minimal==1)exit
3026 
3027  end do
3028 
3029  if (iiter >= maxiter) then
3030    write(message, '(a,i0,a)') 'the second loop has failed to find a set of minimal vectors in ',maxiter, 'iterations.'
3031    MSG_BUG(message)
3032  end if
3033 
3034 !DEBUG
3035 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
3036 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' minim =',minim(:,1),ch10,minim(:,2),ch10,minim(:,3)
3037 !ENDDEBUG
3038 
3039 !DEBUG
3040 !Change sign of the third vector if not right-handed basis
3041 !determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+&
3042 !&            minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+&
3043 !&            minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3))
3044 !write(std_out,*)' smallprim: determinant=',determinant
3045 !ENDDEBUG
3046 
3047 !Choose the first vector
3048 !Compute the scalar product of the three minimal vectors with the first original vector
3049  scprod(:)=zero
3050  do ii=1,3
3051    scprod(:)=scprod(:)+minim(ii,:)*rprimd(ii,1)
3052  end do
3053 !Determine the vector with the maximal absolute overlap
3054  itrial=1
3055  if(abs(scprod(2))>abs(scprod(1))+tol8)itrial=2
3056  if(abs(scprod(3))>abs(scprod(itrial))+tol8)itrial=3
3057 !Switch the vectors if needed
3058  if(itrial/=1)then
3059    tmpvect(:)=minim(:,1)
3060    minim(:,1)=minim(:,itrial)
3061    minim(:,itrial)=tmpvect(:)
3062  end if
3063 !Choose the sign
3064  if(scprod(itrial)<tol8)minim(:,1)=-minim(:,1)
3065 
3066 !DEBUG
3067 !Change sign of the third vector if not right-handed basis
3068 !determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+&
3069 !&            minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+&
3070 !&            minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3))
3071 !write(std_out,*)' smallprim: determinant=',determinant
3072 !ENDDEBUG
3073 
3074 !Choose the second vector
3075 !Compute the scalar product of the second and third minimal vectors with the second original vector
3076  scprod(2:3)=zero
3077  do ii=1,3
3078    scprod(2:3)=scprod(2:3)+minim(ii,2:3)*rprimd(ii,2)
3079  end do
3080 !Determine the vector with the maximal absolute overlap
3081  itrial=2
3082  if(abs(scprod(3))>abs(scprod(2))+tol8)itrial=3
3083 !Switch the vectors if needed
3084  if(itrial/=2)then
3085    tmpvect(:)=minim(:,2)
3086    minim(:,2)=minim(:,itrial)
3087    minim(:,itrial)=tmpvect(:)
3088  end if
3089 !Choose the sign
3090  if(scprod(itrial)<tol8)minim(:,2)=-minim(:,2)
3091 
3092 !Change sign of the third vector if not right-handed basis
3093  determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+&
3094 & minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+&
3095 & minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3))
3096  if(determinant<-tol8)minim(:,3)=-minim(:,3)
3097  if(abs(determinant)<tol8)then
3098    MSG_BUG('minim gives vanishing unit cell volume.')
3099  end if
3100 
3101 !Final computation of metmin
3102  do ii=1,3
3103    metmin(ii,:)=minim(1,ii)*minim(1,:)+&
3104 &   minim(2,ii)*minim(2,:)+&
3105 &   minim(3,ii)*minim(3,:)
3106  end do
3107 
3108 !DEBUG
3109 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
3110 !write(std_out,'(a,3es16.8,a,3es16.8,a,3es16.8)')' minim =',minim(:,1),ch10,minim(:,2),ch10,minim(:,3)
3111 !write(std_out,'(a,3es16.8,a,3es16.8,a,3es16.8)')' metmin =',metmin(:,1),ch10,metmin(:,2),ch10,metmin(:,3)
3112 !ENDDEBUG
3113 
3114 end subroutine smallprim

m_symtk/symatm [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 symatm

FUNCTION

 For each symmetry operation, find the number of the position to
 which each atom is sent in the unit cell by the INVERSE of the
 symmetry operation inv(symrel); i.e. this is the atom which, when acted
 upon by the given symmetry element isym, gets transformed into atom iatom.
 This routine uses the fact that inv(symrel)=trans(symrec),
 the inverse of the symmetry operation expressed in the basis of real
 space primitive translations equals the transpose of the same symmetry
 operation expressed in the basis of reciprocal space primitive transl.
 $xred(nu,indsym(4,isym,ia))=symrec(mu,nu,isym)*(xred(mu,ia)-tnons(mu,isym))
 - transl(mu)$ where $transl$ is also a set of integers and
 where translation transl places coordinates within unit cell (note sign).
 Note that symrec is the set of arrays which are actually input here.
 These arrays have integer elements.
 tnons is the nonsymmorphic translation or else is zero.
 If nsym=1 (i.e. only the identity symmetry is present) then
 indsym merely takes each atom into itself.
 The array of integer translations "transl" gets included within array "indsym" as seen below.
 This routine has been improved using ideas of p. 649 of notes,
 implementing suggestion of Andrew Horsfield: replace search for
 equivalent atoms using direct primitive cell translations by
 use of dot product relation which must produce an integer.
 Relation: $[inv(S(i))*(x(a)-tnons(i)) - x(inv(S)(i,a))] = integer$
 where $S(i) =$ symmetry matrix in real space, tnons=nonsymmorphic translation
 (may be 0 0 0), and $x(inv(S)(i,a))$ is sought atom into which $x(a)$ gets
 rotated by $inv(S)$.  Integer gives primitive translation coordinates to get
 back to original unit cell.
 Equivalent to $S*t(b)+tnons-x(a)=another$ $integer$ for $x(b)=x(inv(S))$.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

 natom=number of atoms in cell.
 nsym=number of space group symmetries.
 symrec(3,3,nsym)=symmetries expressed in terms of their action on
                  reciprocal space primitive translations (integer).
 tnons(3,nsym)=nonsymmorphic translations for each symmetry (would
               be 0 0 0 each for a symmorphic space group)
 typat(natom)=integer identifying type of atom.
 xred(3,natom)=reduced coordinates of atoms in terms of real space
               primitive translations
 tolsym=tolerance for the symmetries

OUTPUT

 indsym(4,nsym,natom)=indirect indexing array described above: for each
                      isym,iatom, fourth element is label of atom into
                      which iatom is sent by INVERSE of symmetry operation
                      isym; first three elements are the primitive
                      translations which must be subtracted after the
                      transformation to get back to the original unit cell.

PARENTS

      get_npert_rbz,ingeo,initberry,initorbmag,m_ab7_symmetry,m_crystal,m_ddb
      m_polynomial_coeff,m_tdep_sym,setsym,thmeig

CHILDREN

      symchk,wrtout

SOURCE

1887 subroutine symatm(indsym,natom,nsym,symrec,tnons,tolsym,typat,xred)
1888 
1889 
1890 !This section has been created automatically by the script Abilint (TD).
1891 !Do not modify the following lines by hand.
1892 #undef ABI_FUNC
1893 #define ABI_FUNC 'symatm'
1894 !End of the abilint section
1895 
1896  implicit none
1897 
1898 !Arguments ------------------------------------
1899 !scalars
1900  integer,intent(in) :: natom,nsym
1901  real(dp), intent(in) :: tolsym
1902 !arrays
1903  integer,intent(in) :: symrec(3,3,nsym),typat(natom)
1904  integer,intent(out) :: indsym(4,nsym,natom)
1905  real(dp),intent(in) :: tnons(3,nsym),xred(3,natom)
1906 
1907 !Local variables-------------------------------
1908 !scalars
1909  integer :: eatom,errout,iatom,ii,isym,mu
1910  real(dp) :: difmax,err
1911  character(len=500) :: message
1912 !arrays
1913  integer :: transl(3)
1914  real(dp) :: difmin(3),tratom(3)
1915 
1916 ! *************************************************************************
1917 
1918  err=zero
1919  errout=0
1920 
1921  do isym=1,nsym
1922    do iatom=1,natom
1923 
1924      do mu=1,3 ! Apply inverse transformation to original coordinates. Note transpose of symrec.
1925        tratom(mu) = dble(symrec(1,mu,isym))*(xred(1,iatom)-tnons(1,isym))&
1926 &       +dble(symrec(2,mu,isym))*(xred(2,iatom)-tnons(2,isym))&
1927 &       +dble(symrec(3,mu,isym))*(xred(3,iatom)-tnons(3,isym))
1928      end do
1929 !
1930 !    Find symmetrically equivalent atom
1931      call symchk(difmin,eatom,natom,tratom,transl,typat(iatom),typat,xred)
1932 !
1933 !    Put information into array indsym: translations and label
1934      indsym(1,isym,iatom)=transl(1)
1935      indsym(2,isym,iatom)=transl(2)
1936      indsym(3,isym,iatom)=transl(3)
1937      indsym(4,isym,iatom)=eatom
1938 !
1939 !    Keep track of maximum difference between transformed coordinates and
1940 !    nearest "target" coordinate
1941      difmax=max(abs(difmin(1)),abs(difmin(2)),abs(difmin(3)))
1942      err=max(err,difmax)
1943 
1944      if (difmax>tolsym) then ! Print warnings if differences exceed tolerance
1945        write(message, '(3a,i3,a,i6,a,i3,a,a,3es12.4,3a)' )&
1946 &       'Trouble finding symmetrically equivalent atoms',ch10,&
1947 &       'Applying inv of symm number',isym,' to atom number',iatom,'  of typat',typat(iatom),ch10,&
1948 &       'gives tratom=',tratom(1:3),'.',ch10,&
1949 &       'This is further away from every atom in crystal than the allowed tolerance.'
1950        MSG_WARNING(message)
1951 
1952        write(message, '(a,3i3,a,a,3i3,a,a,3i3)' ) &
1953 &       '  The inverse symmetry matrix is',symrec(1,1:3,isym),ch10,&
1954 &       '                                ',symrec(2,1:3,isym),ch10,&
1955 &       '                                ',symrec(3,1:3,isym)
1956        call wrtout(std_out,message,'COLL')
1957        write(message, '(a,3f13.7)' )'  and the nonsymmorphic transl. tnons =',(tnons(mu,isym),mu=1,3)
1958 
1959        call wrtout(std_out,message,'COLL')
1960        write(message, '(a,1p,3e11.3,a,a,i5)' ) &
1961 &       '  The nearest coordinate differs by',difmin(1:3),ch10,&
1962 &       '  for indsym(nearest atom)=',indsym(4,isym,iatom)
1963        call wrtout(std_out,message,'COLL')
1964 !
1965 !      Use errout to reduce volume of error diagnostic output
1966        if (errout==0) then
1967          write(message,'(6a)') ch10,&
1968 &         '  This indicates that when symatm attempts to find atoms  symmetrically',ch10, &
1969 &         '  related to a given atom, the nearest candidate is further away than some',ch10,&
1970 &         '  tolerance.  Should check atomic coordinates and symmetry group input data.'
1971          call wrtout(std_out,message,'COLL')
1972          errout=1
1973        end if
1974 
1975      end if !difmax>tol
1976    end do !iatom
1977  end do !isym
1978 
1979  if (.FALSE.) then
1980    do iatom=1,natom
1981      write(message, '(a,i0,a)' )' symatm: atom number ',iatom,' is reached starting at atom'
1982      call wrtout(std_out,message,'COLL')
1983      do ii=1,(nsym-1)/24+1
1984        if(natom<100)then
1985          write(message, '(1x,24i3)' ) (indsym(4,isym,iatom),isym=1+(ii-1)*24,min(nsym,ii*24))
1986        else
1987          write(message, '(1x,24i6)' ) (indsym(4,isym,iatom),isym=1+(ii-1)*24,min(nsym,ii*24))
1988        end if
1989        call wrtout(std_out,message,'COLL')
1990      end do
1991    end do
1992  end if
1993 
1994  if (err>tolsym) then
1995    write(message, '(1x,a,1p,e14.5,a,e12.4)' )'symatm: maximum (delta t)=',err,' is larger than tol=',tolsym
1996    call wrtout(std_out,message,'COLL')
1997  end if
1998 
1999 !Stop execution if error is really big
2000  if (err>0.01d0) then
2001    write(message,'(5a)')&
2002 &   'Largest error (above) is so large (0.01) that either input atomic coordinates (xred)',ch10,&
2003 &   'are wrong or space group symmetry data is wrong.',ch10,&
2004 &   'Action: correct your input file.'
2005    MSG_ERROR(message)
2006  end if
2007 
2008 end subroutine symatm

m_symtk/symaxes [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 symaxes

FUNCTION

 Determines the type of symmetry operation, for
 the proper symmetries 2,2_1,3,3_1,3_2,4,4_1,4_2,4_3,6,6_1,...6_5

INPUTS

 center=type of bravais lattice 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
 iholohedry=type of holohedry
            iholohedry=1   triclinic      1bar
            iholohedry=2   monoclinic     2/m
            iholohedry=3   orthorhombic   mmm
            iholohedry=4   tetragonal     4/mmm
            iholohedry=5   trigonal       3bar m  (rhombohedral Bravais latt)
            iholohedry=6   hexagonal      6/mmm
            iholohedry=7   cubic          m3bar m
 isym=number of the symmetry operation that is currently analyzed
 isymrelconv=symrel matrix for the particular operation, in conv. axes
 ordersym=order of the symmetry operation
 tnons_order=order of the screw translation
 trialt(3)=screw translation associated with the symmetry operation
           in conventional axes (all components in the range ]-1/2,1/2] )

OUTPUT

 label=a user friendly label for the rotation
 type_axis=type of the symmetry operation

NOTES

 It is assumed that the symmetry operations will be entered in the
 symrel tnonsconv arrays, for the CONVENTIONAL cell.
 For proper symmetries (rotations), the
 associated translation is determined.

 There is a subtlety with translations associated with rotations :
 all the rotations with axis
 parallel to the one analysed do not all have the
 same translation characteristics. This is clearly seen
 in the extended Hermann-Mauguin symbols, see the international
 table for crystallography, chapter 4.
 In the treatment that we adopt, one will distinguish
 the cases of primitive Bravais lattices, and centered
 bravais lattices. In the latter case, in the present routine,
 at the exception of the trigonal axis for the
 cubic system, we explicitely generate the correct ratio of different
 translations, so that their type can be explicitely assigned,
 without confusion. By contrast, for primitive lattices,
 the "tnons" that has been transmitted to the present routine
 might be one of the few possible translations vectors,
 nearly at random. We deal with this case by the explicit
 examination of the system classes, and the identification
 of such a possibility. In particular:
 (1) for the trigonal axis in the rhombohedral Bravais lattice,
 or in the cubic system, there is an equal number of 3, 3_1,
 and 3_2 axes parallel to each other, in a cell that
 is primitive (as well as conventional). In this particular case,
 in the present
 routine, all 3, 3_1 and 3_2 axes are assigned to be 3 axes,
 independently of the centering.
 (2) for the 4- or 6- axes, no confusion is possible :
 in the primitive cell, there is only one possible translation,
 while in the centered cells, the correct ratio of translation
 vectors will be generated
 (3) for the binary axes, there is no problem when the cell
 is centered, but there are problems
 (3a) for the tP Bravais lattice, for an axis in a tertiary direction,
 (see the description of the lattice symmetry directions
  table 2.4.1 of the international tables for crystallography),
  where the family of axes is made equally of 2 and 2_1 axis.
  In this case, we attribute the binary axis to the specific class
  of "tertiary 2-axis". We keep track of the 2 or 2_1
  characteristics of all other binary axes
 (3b) for the tI Bravais lattice, in all the directions,
  there is an equal number of 2 and 2_1 axes. We distinguish
  the primary and secondary family from the tertiary family.
 (3c) for the hP Bravais lattice, each binary axis can present
  no translation or be a screw axis (in the same direction).
  For primary axes, one need the "2" and "2_1" classification,
  while for secondary and tertiary axes, the associated
  translation vector will have not importance.
  However, one will need to distinguish secondary from
  tertiary, and these from primary axes.
  So, this is the most complicated case, for binary axes,
  with the following sets of binary axes : "2", "2_1",
  "secondary 2" and "tertiary 2".
 (3d) for the hR Bravais lattice, each binary axis can present
  no translation or be a screw axis (in the same direction).
  There is no distinction between tertiary axes and other, so that
  we simply assign a binary axis to "2-axis"
 (3e) for the cP lattice, the binary axes along tertiary directions
  can also have different translation vectors, while for the primary
  direction, there is no such ambiguity. So, we will attribute
  tertiary 2 axis to the "tertiary 2-axis" set (there are always 6),
  and attribute 2 and 2_1 primary axes to the corresponding sets.

PARENTS

      symcharac

CHILDREN

      wrtout

SOURCE

2310 subroutine symaxes(center,iholohedry,isym,isymrelconv,label,ordersym,tnons_order,trialt,type_axis)
2311 
2312 
2313 !This section has been created automatically by the script Abilint (TD).
2314 !Do not modify the following lines by hand.
2315 #undef ABI_FUNC
2316 #define ABI_FUNC 'symaxes'
2317 !End of the abilint section
2318 
2319  implicit none
2320 
2321 !Arguments ------------------------------------
2322 !scalars
2323  integer,intent(in) :: center,iholohedry,isym,ordersym,tnons_order
2324  integer,intent(out) :: type_axis
2325  character(len=128),intent(out) :: label
2326 !arrays
2327  integer,intent(in) :: isymrelconv(3,3)
2328  real(dp),intent(in) :: trialt(3)
2329 
2330 !Local variables-------------------------------
2331 !scalars
2332  logical,parameter :: verbose=.FALSE.
2333  character(len=500) :: message
2334  integer :: direction,directiontype
2335  real(dp),parameter :: nzero=1.0d-6
2336 
2337 !**************************************************************************
2338 
2339 !write(std_out,*)' symaxes : enter, isym=',isym
2340 !write(std_out,*)' symaxes : iholohedry, ',iholohedry
2341 !write(std_out,*)' symaxes : center, ',center
2342 
2343  select case(ordersym)
2344 
2345  case(2)                       ! point symmetry 2
2346 !    Must characterize directiontype for cP, tP, tI, and hP Bravais lattices
2347    directiontype=1
2348    if( iholohedry==4 .or. iholohedry==7) then ! tP or cP Bravais lattices
2349      if(abs(isymrelconv(1,1))+ &
2350 &     abs(isymrelconv(2,2))+ &
2351 &     abs(isymrelconv(3,3))  ==1) directiontype=3
2352    else if(iholohedry==6)then   ! hP Bravais lattice
2353      if(sum(isymrelconv(:,:))/=-1 )directiontype=2
2354      if(sum(isymrelconv(:,:))==0 .or. sum(isymrelconv(:,:))==-3 )&
2355 &     directiontype=3
2356 !      directiontype=1 corresponds to a primary axis
2357 !      directiontype=2 corresponds to a tertiary axis
2358 !      directiontype=3 corresponds to a secondary axis
2359    end if
2360 
2361 !    DEBUG
2362 !    write(std_out,*)' directiontype=',directiontype
2363 !    write(std_out,'(a,3i6)' )' isymrelconv(1:3)=',isymrelconv(:,1)
2364 !    write(std_out,'(a,3i6)' )' isymrelconv(4:6)=',isymrelconv(:,2)
2365 !    write(std_out,'(a,3i6)' )' isymrelconv(7:9)=',isymrelconv(:,3)
2366 !    write(std_out,'(a,i)' )' tnons_order=',tnons_order
2367 !    ENDDEBUG
2368 
2369 !    Now, classify the 2 axes
2370    if(directiontype==2)then
2371      type_axis=4                 ! secondary 2  (only in the hP Bravais latt case)
2372      write(label,'(a)') 'a secondary 2-axis '
2373 
2374    else if(directiontype==3 .and. iholohedry==4)then
2375      type_axis=21                ! tertiary 2
2376      write(label,'(a)') 'a tertiary 2-axis '
2377    else if(directiontype==3 .and. &
2378 &     center==0 .and. (iholohedry==6.or.iholohedry==7) )then
2379      type_axis=21                ! tertiary 2
2380      write(label,'(a)') 'a tertiary 2-axis '
2381    else if(tnons_order==1 .or. (iholohedry==4 .and. center==-1) .or. &
2382 &     iholohedry==5)then
2383      type_axis=9                 ! 2
2384      write(label,'(a)') 'a 2-axis '
2385    else
2386      type_axis=20                ! 2_1
2387      write(label,'(a)') 'a 2_1-axis '
2388    end if
2389 
2390  case(3)                       ! point symmetry 3
2391    if(tnons_order==1)then
2392      type_axis=10                ! 3
2393      write(label,'(a)') 'a 3-axis '
2394    else if(iholohedry==5 .or. iholohedry==7)then
2395 !      This is a special situation : in the same family of parallel 3-axis,
2396 !      one will have an equal number of 3, 3_1 and 3_2 axes, so that
2397 !      it is non-sense to try to classify one of them.
2398      type_axis=10                ! 3, 3_1 or 3_2, undistinguishable
2399      write(label,'(a)') 'a 3, 3_1 or 3_2 axis '
2400    else
2401 !      DEBUG
2402 !      write(std_out,*)'isymrelconv=',isymrelconv(:,:)
2403 !      write(std_out,*)'trialt=',trialt(:)
2404 !      ENDDEBUG
2405 !      Must recognize 3_1 or 3_2
2406      if(isymrelconv(1,1)==0)then  ! 3+
2407        if(abs(trialt(3)-third)<nzero)type_axis=22   ! 3_1
2408        if(abs(trialt(3)+third)<nzero)type_axis=23   ! 3_2
2409      else if(isymrelconv(1,1)==-1)then  ! 3-
2410        if(abs(trialt(3)-third)<nzero)type_axis=23   ! 3_2
2411        if(abs(trialt(3)+third)<nzero)type_axis=22   ! 3_1
2412      end if
2413      write(label,'(a)') 'a 3_1 or 3_2-axis '
2414    end if
2415 
2416  case(4)                       ! point symmetry 4
2417    if(tnons_order==1)then
2418      type_axis=12                ! 4
2419      write(label,'(a)') 'a 4-axis '
2420    else if(tnons_order==2)then
2421      type_axis=25                ! 4_2
2422      write(label,'(a)') 'a 4_2-axis '
2423    else if(center/=0)then
2424      type_axis=24                ! 4_1 or 4_3
2425      write(label,'(a)') 'a 4_1 or 4_3-axis '
2426    else
2427 !      DEBUG
2428 !      write(std_out,*)'isymrelconv=',isymrelconv(:,:)
2429 !      write(std_out,*)'trialt=',trialt(:)
2430 !      ENDDEBUG
2431 !      Must recognize 4_1 or 4_3, along the three primary directions
2432      do direction=1,3
2433        if(isymrelconv(direction,direction)==1)then  !
2434          if( (direction==1 .and. isymrelconv(2,3)==-1) .or. &
2435 &         (direction==2 .and. isymrelconv(3,1)==-1) .or. &
2436 &         (direction==3 .and. isymrelconv(1,2)==-1)       )then ! 4+
2437            if(abs(trialt(direction)-quarter)<nzero)type_axis=24    ! 4_1
2438            if(abs(trialt(direction)+quarter)<nzero)type_axis=26    ! 4_3
2439          else if( (direction==1 .and. isymrelconv(2,3)==1) .or. &
2440 &           (direction==2 .and. isymrelconv(3,1)==1) .or. &
2441 &           (direction==3 .and. isymrelconv(1,2)==1)       )then ! 4-
2442            if(abs(trialt(direction)-quarter)<nzero)type_axis=26    ! 4_3
2443            if(abs(trialt(direction)+quarter)<nzero)type_axis=24    ! 4_1
2444          end if
2445        end if
2446      end do
2447      write(label,'(a)') 'a 4_1 or 4_3-axis '
2448    end if
2449 
2450  case(6)                       ! point symmetry 6
2451    if(tnons_order==1)then
2452      type_axis=14                ! 6
2453      write(label,'(a)') 'a 6-axis '
2454    else if(tnons_order==2)then
2455      type_axis=29                ! 6_3
2456      write(label,'(a)') 'a 6_3-axis '
2457    else if(tnons_order==3)then
2458 !      DEBUG
2459 !      write(std_out,*)'isymrelconv=',isymrelconv(:,:)
2460 !      write(std_out,*)'trialt=',trialt(:)
2461 !      ENDDEBUG
2462 !      Must recognize 6_2 or 6_4
2463      if(isymrelconv(1,1)==1)then  ! 6+
2464        if(abs(trialt(3)-third)<nzero)type_axis=28   ! 6_2
2465        if(abs(trialt(3)+third)<nzero)type_axis=30   ! 6_4
2466      else if(isymrelconv(1,1)==0)then  ! 6-
2467        if(abs(trialt(3)-third)<nzero)type_axis=30   ! 6_4
2468        if(abs(trialt(3)+third)<nzero)type_axis=28   ! 6_2
2469      end if
2470      write(label,'(a)') 'a 6_2 or 6_4-axis '
2471    else
2472 !      DEBUG
2473 !      write(std_out,*)'isymrelconv=',isymrelconv(:,:)
2474 !      write(std_out,*)'trialt=',trialt(:)
2475 !      ENDDEBUG
2476 !      Must recognize 6_1 or 6_5
2477      if(isymrelconv(1,1)==1)then  ! 6+
2478        if(abs(trialt(3)-sixth)<nzero)type_axis=27   ! 6_1
2479        if(abs(trialt(3)+sixth)<nzero)type_axis=31   ! 6_5
2480      else if(isymrelconv(1,1)==0)then  ! 6-
2481        if(abs(trialt(3)-sixth)<nzero)type_axis=31   ! 6_5
2482        if(abs(trialt(3)+sixth)<nzero)type_axis=27   ! 6_1
2483      end if
2484      write(label,'(a)') 'a 6_1 or 6_5-axis '
2485    end if
2486 
2487  end select
2488 
2489  if (verbose) then
2490    write(message,'(a,i3,a,a)')' symaxes : the symmetry operation no. ',isym,' is ', trim(label)
2491    call wrtout(std_out,message,'COLL')
2492  end if
2493 
2494 end subroutine symaxes

m_symtk/symcharac [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 symcharac

FUNCTION

 Get the type of axis for the symmetry.

INPUTS

 center=bravais(2)
 determinant=the value of the determinant of sym
 iholohedry=bravais(1)
 isym=number of the symmetry operation that is currently analyzed
 order=the order of the symmetry
 symrel(3,3)= the symmetry matrix
 tnons(3)=nonsymmorphic translations

OUTPUT

 label=a human readable text for the characteristic of the symmetry
 type_axis=an identifier for the type of symmetry

PARENTS

      m_ab7_symmetry,symspgr

CHILDREN

      symaxes,symplanes,wrtout

SOURCE

2039 subroutine symcharac(center, determinant, iholohedry, isym, label, symrel, tnons, type_axis)
2040 
2041 
2042 !This section has been created automatically by the script Abilint (TD).
2043 !Do not modify the following lines by hand.
2044 #undef ABI_FUNC
2045 #define ABI_FUNC 'symcharac'
2046 !End of the abilint section
2047 
2048  implicit none
2049 
2050 !Arguments ------------------------------------
2051 !scalars
2052  integer, intent(in) :: determinant, center, iholohedry, isym
2053  integer, intent(out) :: type_axis
2054  character(len = 128), intent(out) :: label
2055  !arrays
2056  integer,intent(in) :: symrel(3,3)
2057  real(dp),intent(in) :: tnons(3)
2058 
2059  !Local variables-------------------------------
2060  !scalars
2061  logical,parameter :: verbose=.FALSE.
2062  integer :: tnons_order, identified, ii, order, iorder
2063  character(len=500) :: message
2064  !arrays
2065  integer :: identity(3,3),matrix(3,3),trial(3,3)
2066  real(dp) :: reduced(3),trialt(3)
2067 
2068  !**************************************************************************
2069 
2070  identity(:,:)=0
2071  identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
2072  trial(:,:)=identity(:,:)
2073  matrix(:,:)=symrel(:,:)
2074 
2075  order=0
2076  do iorder=1,6
2077    trial=matmul(matrix,trial)
2078    if(sum((trial-identity)**2)==0)then
2079      order=iorder
2080      exit
2081    end if
2082    if(sum((trial+identity)**2)==0)then
2083      order=iorder
2084      exit
2085    end if
2086  end do
2087 
2088  if(order==0)then
2089    type_axis = -2
2090    return
2091  end if
2092 
2093 !Determination of the characteristics of proper symmetries (rotations)
2094  if(determinant==1)then
2095 
2096 !  Determine the translation vector associated to the rotations
2097 !  and its order : apply the symmetry operation
2098 !  then analyse the resulting vector.
2099    identified=0
2100    trialt(:)=zero
2101    do ii=1,order
2102      trialt(:)=matmul(symrel(:,:),trialt(:))+tnons(:)
2103    end do
2104 !  Gives the associated translation, with components in the
2105 !  interval [-0.5,0.5] .
2106    reduced(:)=trialt(:)-nint(trialt(:)-tol6)
2107 
2108    if(sum(abs(reduced(:)))<tol6)identified=1
2109    if( (center==1 .or. center==-3) .and. &
2110 &   sum(abs(reduced(:)-(/zero,half,half/)))<tol6 )identified=2
2111    if( (center==2 .or. center==-3) .and. &
2112 &   sum(abs(reduced(:)-(/half,zero,half/)))<tol6 )identified=3
2113    if( (center==3 .or. center==-3) .and. &
2114 &   sum(abs(reduced(:)-(/half,half,zero/)))<tol6 )identified=4
2115    if(center==-1.and. sum(abs(reduced(:)-(/half,half,half/)))<tol6 )identified=5
2116 
2117 !  If the symmetry operation has not been identified, there is a problem ...
2118    if(identified==0) then
2119      type_axis = -1
2120      return
2121    end if
2122 
2123 !  Compute the translation vector associated with one rotation
2124    trialt(:)=trialt(:)/order
2125    trialt(:)=trialt(:)-nint(trialt(:)-tol6)
2126 
2127 !  Analyse the resulting vector.
2128    identified=0
2129    do ii=1,order
2130      reduced(:)=ii*trialt(:)-nint(ii*trialt(:)-tol6)
2131      if(sum(abs(reduced(:)))<tol6)identified=1
2132      if( (center==1 .or. center==-3) .and. &
2133 &     sum(abs(reduced(:)-(/zero,half,half/)))<tol6 )identified=2
2134      if( (center==2 .or. center==-3) .and. &
2135 &     sum(abs(reduced(:)-(/half,zero,half/)))<tol6 )identified=3
2136      if( (center==3 .or. center==-3) .and. &
2137 &     sum(abs(reduced(:)-(/half,half,zero/)))<tol6 )identified=4
2138      if(center==-1.and. sum(abs(reduced(:)-(/half,half,half/)))<tol6 )identified=5
2139 
2140      if(identified/=0)then
2141        tnons_order=ii
2142        exit
2143      end if
2144    end do ! ii
2145 
2146 !  Determinant (here=+1, as we are dealing with proper symmetry operations),
2147 !  order, tnons_order and identified are enough to
2148 !  determine the kind of symmetry operation
2149 
2150    select case(order)
2151    case(1)                       ! point symmetry 1
2152      if(identified==1) then
2153        type_axis=8                 ! 1
2154        write(label,'(a)') 'the identity'
2155      else
2156        type_axis=7                 ! t
2157        write(label,'(a)') 'a pure translation '
2158      end if
2159 
2160      if (verbose) then
2161        write(message,'(a,i3,2a)')' symspgr : the symmetry operation no. ',isym,' is ',trim(label)
2162        call wrtout(std_out,message,'COLL')
2163      end if
2164 
2165    case(2,3,4,6)                 ! point symmetry 2,3,4,6 - rotations
2166      call symaxes(center,iholohedry,isym,symrel,label,order,tnons_order,trialt,type_axis)
2167    end select
2168 
2169  else if (determinant==-1)then
2170 
2171 !  Now, take care of the improper symmetry operations.
2172 !  Their treatment is relatively easy, except for the mirror planes
2173    select case(order)
2174    case(1)                       ! point symmetry 1
2175      type_axis=5                  ! -1
2176      write(label,'(a)') 'an inversion'
2177    case(2)                       ! point symmetry 2 - planes
2178      call symplanes(center,iholohedry,isym,symrel,tnons,label,type_axis)
2179    case(3)                       ! point symmetry 3
2180      type_axis=3                  ! -3
2181      write(label,'(a)') 'a -3 axis '
2182    case(4)                       ! point symmetry 1
2183      type_axis=2                  ! -4
2184      write(label,'(a)') 'a -4 axis '
2185    case(6)                       ! point symmetry 1
2186      type_axis=1                  ! -6
2187      write(label,'(a)') 'a -6 axis '
2188    end select
2189 
2190    if (order /= 2 .and. verbose) then
2191      write(message,'(a,i3,2a)')' symspgr : the symmetry operation no. ',isym,' is ',trim(label)
2192      call wrtout(std_out,message,'COLL')
2193    end if
2194 
2195  end if ! determinant==1 or -1
2196 
2197 end subroutine symcharac

m_symtk/symchk [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 symchk

FUNCTION

 Symmetry checker for atomic coordinates.
 Checks for translated atomic coordinate tratom(3) to agree
 with some coordinate xred(3,iatom) where atomic types agree too.
 All coordinates are "reduced", i.e. given in terms of primitive
 reciprocal translations.

INPUTS

 natom=number of atoms in unit cell
 tratom(3)=reduced coordinates for a single atom which presumably
   result from the application of a symmetry operation to an atomic
   coordinate
 trtypat=type of atom (integer) translated to tratom
 typat(natom)=types of all atoms in unit cell (integer)
 xred(3,natom)=reduced coordinates for all atoms in unit cell

OUTPUT

 difmin(3)=minimum difference between apparently equivalent atoms
   (give value separately for each coordinate)--note that value
   may be NEGATIVE so take abs later if needed
 eatom=atom label of atom which is SAME as tratom to within a primitive
   cell translation ("equivalent atom")
 transl(3)=primitive cell translation to make iatom same as tratom (integers)

PARENTS

      m_polynomial_coeff,symatm

CHILDREN

SOURCE

1732 subroutine symchk(difmin,eatom,natom,tratom,transl,trtypat,typat,xred)
1733 
1734 
1735 !This section has been created automatically by the script Abilint (TD).
1736 !Do not modify the following lines by hand.
1737 #undef ABI_FUNC
1738 #define ABI_FUNC 'symchk'
1739 !End of the abilint section
1740 
1741  implicit none
1742 
1743 !Arguments ------------------------------------
1744 !scalars
1745  integer,intent(in) :: natom,trtypat
1746  integer,intent(out) :: eatom
1747 !arrays
1748  integer,intent(in) :: typat(natom)
1749  integer,intent(out) :: transl(3)
1750  real(dp),intent(in) :: tratom(3),xred(3,natom)
1751  real(dp),intent(out) :: difmin(3)
1752 
1753 !Local variables-------------------------------
1754 !scalars
1755  integer :: iatom,jatom,trans1,trans2,trans3
1756  real(dp) :: test,test1,test2,test3,testmn
1757 
1758 ! *************************************************************************
1759 
1760 !Start testmn out at large value
1761  testmn=1000000.d0
1762 
1763 !Loop through atoms--
1764 !when types agree, check for agreement after primitive translation
1765  jatom=1
1766  do iatom=1,natom
1767    if (trtypat/=typat(iatom)) cycle
1768 
1769 !  Check all three components
1770    test1=tratom(1)-xred(1,iatom)
1771    test2=tratom(2)-xred(2,iatom)
1772    test3=tratom(3)-xred(3,iatom)
1773 !  Find nearest integer part of difference
1774    trans1=nint(test1)
1775    trans2=nint(test2)
1776    trans3=nint(test3)
1777 !  Check whether, after translation, they agree
1778    test1=test1-dble(trans1)
1779    test2=test2-dble(trans2)
1780    test3=test3-dble(trans3)
1781 
1782    test=abs(test1)+abs(test2)+abs(test3)
1783    if (test<tol10) then
1784 !    Note that abs() is not taken here
1785      difmin(1)=test1
1786      difmin(2)=test2
1787      difmin(3)=test3
1788      jatom=iatom
1789      transl(1)=trans1
1790      transl(2)=trans2
1791      transl(3)=trans3
1792 !    Break out of loop when agreement is within tolerance
1793      exit
1794    else
1795 !    Keep track of smallest difference if greater than tol10
1796      if (test<testmn) then
1797        testmn=test
1798 !      Note that abs() is not taken here
1799        difmin(1)=test1
1800        difmin(2)=test2
1801        difmin(3)=test3
1802        jatom=iatom
1803        transl(1)=trans1
1804        transl(2)=trans2
1805        transl(3)=trans3
1806      end if
1807    end if
1808 
1809 !  End loop over iatom. Note a "cycle" and an "exit" inside the loop
1810  end do
1811 
1812  eatom=jatom
1813 
1814 end subroutine symchk

m_symtk/symdet [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 symdet

FUNCTION

 Compute determinant of each input symmetry matrix sym(3,3,i)
 and check that the determinant is always +/- 1.  Integer arithmetic.

INPUTS

 nsym=number of symmetry operations
 sym(3,3,nsym)=integer symmetry array

OUTPUT

 determinant(nsym)=determinant of each symmetry operation

PARENTS

      remove_inversion,setsym,symptgroup,symspgr

CHILDREN

      mati3det

SOURCE

300 subroutine symdet(determinant, nsym, sym)
301 
302 
303 !This section has been created automatically by the script Abilint (TD).
304 !Do not modify the following lines by hand.
305 #undef ABI_FUNC
306 #define ABI_FUNC 'symdet'
307 !End of the abilint section
308 
309  implicit none
310 
311 !Arguments ------------------------------------
312 !scalars
313  integer,intent(in) :: nsym
314 !arrays
315  integer,intent(in) :: sym(3,3,nsym)
316  integer,intent(out) :: determinant(nsym)
317 
318 !Local variables-------------------------------
319 !scalars
320  integer :: det,isym
321  character(len=500) :: message
322 
323 ! *************************************************************************
324 
325  do isym=1,nsym
326    call mati3det(sym(:,:,isym),det)
327    determinant(isym)=det
328    if (abs(det)/=1) then
329      write(message,'(a,i5,a,i10,a,a,a,a,a)')&
330 &     'Abs(determinant) for symmetry number',isym,' is',det,' .',ch10,&
331 &     'For a legitimate symmetry, abs(determinant) must be 1.',ch10,&
332 &     'Action: check your symmetry operations (symrel) in input file.'
333      MSG_ERROR(message)
334    end if
335  end do
336 
337 end subroutine symdet

m_symtk/symmetrize_rprimd [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 symmetrize_rprimd

FUNCTION

 Supposing the input rprimd does not preserve the length and angles
 following the symmetries, will generates a new set rprimd,
 on the basis of the expected characteristics of the conventional cell,
 as specified in bravais(:)

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,1:nsym)=symmetry operations in real space in terms of primitive translations
 tolsym=tolerance for the symmetry operations (only for checking purposes, the new set rprimd will
     be coherent with the symmetry operations at a much accurate level).

SIDE EFFECTS

 rprimd(3,3)=dimensional primitive translations for real space (bohr)

PARENTS

      ingeo

CHILDREN

      chkorthsy,holocell,matr3inv,metric

SOURCE

1469 subroutine symmetrize_rprimd(bravais,nsym,rprimd,symrel,tolsym)
1470 
1471 
1472 !This section has been created automatically by the script Abilint (TD).
1473 !Do not modify the following lines by hand.
1474 #undef ABI_FUNC
1475 #define ABI_FUNC 'symmetrize_rprimd'
1476 !End of the abilint section
1477 
1478  implicit none
1479 
1480 !Arguments ------------------------------------
1481 !scalars
1482  integer,intent(in) :: nsym
1483  real(dp),intent(in) :: tolsym
1484 !arrays
1485  integer,intent(in) :: bravais(11),symrel(3,3,nsym)
1486  real(dp),intent(inout) :: rprimd(3,3)
1487 
1488 !Local variables-------------------------------
1489 !scalars
1490  integer :: foundc,iexit,ii,jj
1491  real(dp):: reldiff
1492  character(len=500) :: msg
1493 !arrays
1494  real(dp):: aa(3,3),ait(3,3),cell_base(3,3),gprimd(3,3),rmet(3,3),rprimd_new(3,3)
1495 
1496 ! *************************************************************************
1497 
1498 !Build the conventional cell basis vectors in cartesian coordinates
1499  aa(:,1)=bravais(3:5)
1500  aa(:,2)=bravais(6:8)
1501  aa(:,3)=bravais(9:11)
1502 !Inverse transpose
1503  call matr3inv(aa,ait)
1504  do ii=1,3
1505    cell_base(:,ii)=ait(ii,1)*rprimd(:,1)+ait(ii,2)*rprimd(:,2)+ait(ii,3)*rprimd(:,3)
1506  end do
1507 
1508 !Enforce the proper holohedry on the conventional cell vectors.
1509  call holocell(cell_base,1,foundc,bravais(1),tolsym)
1510 
1511 !Reconstruct the dimensional primitive vectors
1512  do ii=1,3
1513    rprimd_new(:,ii)=aa(1,ii)*cell_base(:,1)+aa(2,ii)*cell_base(:,2)+aa(3,ii)*cell_base(:,3)
1514  end do
1515 
1516 !Check whether the modification make sense
1517  do ii=1,3
1518    do jj=1,3
1519      reldiff=(rprimd_new(ii,jj)-rprimd(ii,jj))/sqrt(sum(rprimd(:,jj)**2))
1520 !    Allow for twice tolsym
1521      if(abs(reldiff)>two*tolsym)then
1522        write(msg,'(a,6(2a,3es14.6))')&
1523 &       'Failed rectification of lattice vectors to comply with Bravais lattice identification, modifs are too large',ch10,&
1524 &       '  rprimd    =',rprimd(:,1),ch10,&
1525 &       '             ',rprimd(:,2),ch10,&
1526 &       '             ',rprimd(:,3),ch10,&
1527 &       '  rprimd_new=',rprimd_new(:,1),ch10,&
1528 &       '             ',rprimd_new(:,2),ch10,&
1529 &       '             ',rprimd_new(:,3)
1530        MSG_ERROR_CLASS(msg, "TolSymError")
1531      end if
1532    end do
1533  end do
1534 
1535  rprimd(:,:)=rprimd_new(:,:)
1536 
1537 !Check whether the symmetry operations are consistent with the lattice vectors
1538  rmet = MATMUL(TRANSPOSE(rprimd), rprimd)
1539  call matr3inv(rprimd, gprimd)
1540  !call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1541  iexit=1
1542 
1543  call chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel)
1544 
1545 end subroutine symmetrize_rprimd

m_symtk/symmetrize_xred [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 symmetrize_xred

FUNCTION

 Symmetrize atomic coordinates using input symmetry matrices symrel
 which are expressed in terms of the basis of real space primitive
 translations (array elements are integers).
 Input array indsym(4,isym,iatom) gives label of atom into which iatom
 is rotated by INVERSE of symmetry element isym and also gives primitive
 translation to get back to unit cell.
 This version uses improvement in algorithm suggested by Andrew
 Horsfield (see symatm.f).

INPUTS

 indsym(4,nsym,natom)=indirect indexing array giving label of atom
   into which iatom is rotated by symmetry element isym
 natom=number of atoms
 nsym=number of symmetries in group
 symrel(3,3,nsym)=symmetry matrices in terms of real space
   primitive translations
 tnons(3,nsym)=nonsymmorphic translations for symmetries

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/Output
 xred(3,natom)=
  (input) atomic coordinates in terms of real space translations
  (output) symmetrized atomic coordinates in terms
    of real space translations

PARENTS

      ingeo,mover,nonlinear,respfn,scfcv

CHILDREN

SOURCE

1588 subroutine symmetrize_xred(indsym,natom,nsym,symrel,tnons,xred)
1589 
1590 
1591 !This section has been created automatically by the script Abilint (TD).
1592 !Do not modify the following lines by hand.
1593 #undef ABI_FUNC
1594 #define ABI_FUNC 'symmetrize_xred'
1595 !End of the abilint section
1596 
1597  implicit none
1598 
1599 !Arguments ------------------------------------
1600 !scalars
1601  integer,intent(in) :: natom,nsym
1602 !arrays
1603  integer,intent(in) :: indsym(4,nsym,natom),symrel(3,3,nsym)
1604  real(dp),intent(in) :: tnons(3,nsym)
1605  real(dp),intent(inout) :: xred(3,natom)
1606 
1607 !Local variables-------------------------------
1608 !scalars
1609  integer  :: iatom,ib,isym
1610  integer  :: ii,jj
1611  real(dp) :: fc1,fc2,fc3
1612  real(dp) :: diff
1613  logical  :: dissimilar
1614 !arrays
1615  real(dp) :: tsum(3),tt(3)
1616  real(dp),allocatable :: xredsym(:,:)
1617  real(dp) :: transl(3) ! translation vector
1618 
1619 ! *************************************************************************
1620 !
1621 !Check whether group contains more than identity;
1622 !if not then simply return
1623  if (nsym>1) then
1624 
1625 !  loop over atoms
1626    ABI_ALLOCATE(xredsym,(3,natom))
1627    do iatom=1,natom
1628      tsum(1)=0.0d0
1629      tsum(2)=0.0d0
1630      tsum(3)=0.0d0
1631 !
1632 !    loop over symmetries
1633      do isym=1,nsym
1634 !      atom ib is atom into which iatom is rotated by inverse of
1635 !      symmetry isym (inverse of symrel(mu,nu,isym))
1636        ib=indsym(4,isym,iatom)
1637 !      Find the reduced coordinates after translation=t(indsym)+transl
1638        fc1=xred(1,ib)+dble(indsym(1,isym,iatom))
1639        fc2=xred(2,ib)+dble(indsym(2,isym,iatom))
1640        fc3=xred(3,ib)+dble(indsym(3,isym,iatom))
1641 !      Compute [S * (x(indsym)+transl) ] + tnonsymmorphic
1642        tt(:)=dble(symrel(:,1,isym))*fc1+&
1643 &       dble(symrel(:,2,isym))*fc2+&
1644 &       dble(symrel(:,3,isym))*fc3+ tnons(:,isym)
1645 
1646 !      Average over nominally equivalent atomic positions
1647        tsum(:)=tsum(:)+tt(:)
1648      end do
1649 !
1650 !    Set symmetrized result to sum over number of terms
1651      xredsym(:,iatom)=tsum(:)/dble(nsym)
1652 
1653 !    End loop over iatom
1654    end do
1655 
1656    transl(:)=xredsym(:,1)-nint(xredsym(:,1))
1657 
1658 !  Compute the smallest translation to an integer
1659    do jj=2,natom
1660      do ii=1,3
1661        diff=xredsym(ii,jj)-nint(xredsym(ii,jj))
1662        if (diff<transl(ii)) transl(ii)=diff
1663      end do
1664    end do
1665 
1666 !  Test if the translation on each direction is small
1667 !  Tolerance 1E-13
1668    do ii=1,3
1669      if (abs(transl(ii))>1e-13) transl(ii)=0.0
1670    end do
1671 
1672 !  Execute translation
1673    do jj=1,natom
1674      do ii=1,3
1675        xredsym(ii,jj)=xredsym(ii,jj)-transl(ii)
1676      end do
1677    end do
1678 
1679 !  Test if xredsym is too similar to xred
1680 !  Tolerance 1E-15
1681    dissimilar=.FALSE.
1682    do jj=1,natom
1683      do ii=1,3
1684        if (abs(xredsym(ii,jj)-xred(ii,jj))>1E-15) dissimilar=.TRUE.
1685      end do
1686    end do
1687 
1688    if (dissimilar) xred(:,:)=xredsym(:,:)
1689    ABI_DEALLOCATE(xredsym)
1690 
1691 !  End condition of nsym/=1
1692  end if
1693 
1694 end subroutine symmetrize_xred

m_symtk/symplanes [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 symplanes

FUNCTION

 Determines the type of symmetry mirror planes: m,a,b,c,d,n,g.
 This is used (see symlist.f) to identify the space group.

INPUTS

 center=type of bravais lattice 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
 iholohedry=type of holohedry
   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
 isym=number of the symmetry operation that is currently analyzed
 isymrelconv=symrel matrix for the particular operation, in conv. coord.
 itnonsconv=tnons vector for the particular operation, in conv. coord

OUTPUT

 label=user friendly label of the plane
 type_axis=type of the symmetry operation

NOTES

 One follows the
 conventions explained in table 1.3 of the international tables for
 crystallography. In the case of the rhombohedral system,
 one takes into account the first footnote of this table 1.3 .
 In general, we will assign the different symmetries to
 the following numbers :  m -> 15 , (a, b or c) -> 16,
  d -> 17, n -> 18 , g -> 19
 However, there is the same problem as for binary axes,
 namely, for parallel mirror planes, one can find different
 translation vectors, and these might be found at random,
 depending on the input tnons.
 (1) In the tP case, one will distinguish tertiary
  mirror plane, for which it is important to know whether they are
  m or c (for tertiary planes in tP, g is equivalent to m and n is equivalent to c).
  On the other hand, it is important to distinguish among
  primary and secondary mirror planes, those that are m,(a or b),c, or n.
  To summarize, the number of the symmetry will be :
  m (primary, secondary or tertiary) -> 15 ,
  secondary (a or b) -> 16, secondary c -> 17,
  primary or secondary n -> 18 , tertiary c -> 19
 (2) In the tI case, one will distinguish tertiary
  mirror plane, for which it is important to know whether they are
  m or d (for tertiary planes in tI, c is equivalent to m.
  On the other hand, it is important to distinguish among
  primary and secondary mirror planes, those that are m (equivalent to n),
  or a,b or c.
  To summarize, the number of the symmetry will be :
  m (primary, secondary, tertiary) -> 15 ,
  a,b or c (primary or secondary) -> 16, tertiary d -> 17
 (3) For hP and hR, a m plane is always coupled to a a or b plane,
  while a c plane is always coupled to an n plane. On the other
  hand, it is important to distinguish between primary or secondary
  mirror planes, and tertiary mirror planes. So we will keep the
  following sets : m non-tertiary (that includes a or b non-tertiary) -> 15,
  c non-tertiary (that includes n non-tertiary) -> 16,
  m tertiary (that includes a or b non-tertiary) -> 17,
  c tertiary (that includes n non-tertiary) -> 18.
  For hR, all mirror planes are secondary.
 (4) For the cP lattice, in the same spirit, one can see that
  the tertiary m and g mirror planes are to be classified as "m" -> 15,
  while n, a and c are to be classified as "n" -> 18. There is no need
  to distinguish between primary, secondary or tertiary axes.

PARENTS

      symcharac

CHILDREN

      wrtout

SOURCE

2581 subroutine symplanes(center,iholohedry,isym,isymrelconv,itnonsconv,label,type_axis)
2582 
2583 
2584 !This section has been created automatically by the script Abilint (TD).
2585 !Do not modify the following lines by hand.
2586 #undef ABI_FUNC
2587 #define ABI_FUNC 'symplanes'
2588 !End of the abilint section
2589 
2590  implicit none
2591 
2592 !Arguments ------------------------------------
2593 !scalars
2594  integer,intent(in) :: center,iholohedry,isym
2595  integer,intent(out) :: type_axis
2596  character(len = 128), intent(out) :: label
2597 !arrays
2598  integer,intent(in) :: isymrelconv(3,3)
2599  real(dp),intent(in) :: itnonsconv(3)
2600 
2601 !Local variables-------------------------------
2602 !scalars
2603  logical,parameter :: verbose=.FALSE.
2604  character(len=500) :: message
2605  integer :: directiontype,sum_elements
2606  real(dp),parameter :: nzero=1.0d-6
2607 !arrays
2608  integer :: identity(3,3),mirrormxy(3,3),mirrormyz(3,3),mirrormzx(3,3)
2609  integer :: mirrorx(3,3),mirrorxy(3,3),mirrory(3,3),mirroryz(3,3),mirrorz(3,3)
2610  integer :: mirrorzx(3,3)
2611  real(dp) :: trialt(3)
2612 ! real(dp) :: itnonsconv2(3),trialt2(3)
2613 
2614 !**************************************************************************
2615 
2616 !write(std_out,*)' symplanes : enter'
2617 !write(std_out,*)' center,iholohedry,isym,isymrelconv,itnonsconv=',center,iholohedry,isym,isymrelconv,itnonsconv
2618 
2619  identity(:,:)=0
2620  identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
2621 
2622 !Will be a mirror plane, but one must characterize
2623 !(1) the type of plane (primary, secondary or tertiary)
2624 !(2) the gliding vector. One now defines a few matrices.
2625  mirrorx(:,:)=identity(:,:) ; mirrorx(1,1)=-1
2626  mirrory(:,:)=identity(:,:) ; mirrory(2,2)=-1
2627  mirrorz(:,:)=identity(:,:) ; mirrorz(3,3)=-1
2628  mirrorxy(:,:)=0 ; mirrorxy(1,2)=1 ; mirrorxy(2,1)=1 ; mirrorxy(3,3)=1
2629  mirrorzx(:,:)=0 ; mirrorzx(1,3)=1 ; mirrorzx(3,1)=1 ; mirrorzx(2,2)=1
2630  mirroryz(:,:)=0 ; mirroryz(2,3)=1 ; mirroryz(3,2)=1 ; mirroryz(1,1)=1
2631  mirrormxy(:,:)=0 ; mirrormxy(1,2)=-1 ; mirrormxy(2,1)=-1 ; mirrormxy(3,3)=1
2632  mirrormzx(:,:)=0 ; mirrormzx(1,3)=-1 ; mirrormzx(3,1)=-1 ; mirrormzx(2,2)=1
2633  mirrormyz(:,:)=0 ; mirrormyz(2,3)=-1 ; mirrormyz(3,2)=-1 ; mirrormyz(1,1)=1
2634 
2635 !Determine the type of plane. At the end,
2636 !directiontype=1 will correspond to a primary axis (or equivalent
2637 !axes for orthorhombic)
2638 !directiontype=2 will correspond to a secondary axis
2639 !directiontype=3 will correspond to a tertiary axis
2640 !See table 2.4.1, 11.2 and 11.3 of the international tables for crystallography
2641  directiontype=0
2642 !The sum of elements of the matrices allow to characterize them
2643  sum_elements=sum(isymrelconv(:,:))
2644 
2645  if(sum_elements==1)then
2646 !  The mirror plane perpendicular to the c axis is always primary
2647    if( sum(abs(isymrelconv(:,:)-mirrorz(:,:)))==0 )then
2648      directiontype=1
2649 !    All the other planes with a symrel matrix whose sum of elements is 1
2650 !    are a or b planes. They are primary or
2651 !    secondary planes, depending the holohedry.
2652    else if(sum(isymrelconv(:,:))==1)then
2653      if( iholohedry==2 .or. iholohedry==3 .or. iholohedry==7 )then
2654        directiontype=1
2655      else if(iholohedry==4 .or. iholohedry==6)then
2656        directiontype=2
2657      end if
2658    end if
2659  end if
2660 
2661 !All the planes with a symrel matrix whose sum of elements
2662 !is 2 are secondary planes (table 11.3).
2663  if( sum_elements==2 ) directiontype=2
2664 
2665 !The planes with a symrel matrix whose sum of elements
2666 !is 3 or 0 are tertiary planes
2667  if( sum_elements==3 .or. sum_elements==0 )directiontype=3
2668 
2669 !One is left with sum_elements=-1, tertiary for tetragonal
2670 !or cubic, secondary for hexagonal
2671  if( sum_elements==-1)then
2672    if(iholohedry==4 .or. iholohedry==7)directiontype=3
2673    if(iholohedry==6)directiontype=2
2674  end if
2675 
2676 
2677 !Now, determine the gliding vector
2678 !First, apply the symmetry operation
2679 !to itnonsconv, in order to get the translation vector
2680 !under the application of twice the symmetry operation
2681  trialt(:)=matmul(isymrelconv(:,:),itnonsconv(:)) +itnonsconv(:)
2682 !Get the translation associated with one application,
2683 !and force its components to be in the interval ]-0.5,0.5] .
2684  trialt(:)=trialt(:)*half
2685  trialt(:)=trialt(:)-nint(trialt(:)-nzero)
2686 
2687 !If there is a glide vector for the initial choice of itnonsconv,
2688 !it might be that it disappears if itnonsconv is translated by a
2689 !lattice vector of the conventional cell
2690 !if(trialt(1)**2+trialt(2)**2+trialt(3)**2>tol5)then
2691 !do ii=1,3
2692 !itnonsconv2(:)=itnonsconv(:)
2693 !itnonsconv2(ii)=itnonsconv(ii)+one
2694 !trialt2(:)=matmul(isymrelconv(:,:),itnonsconv2(:)) +itnonsconv2(:)
2695 !trialt2(:)=trialt2(:)*half
2696 !trialt2(:)=trialt2(:)-nint(trialt2(:)-nzero)
2697 !if(trialt2(1)**2+trialt2(2)**2+trialt2(3)**2<tol5)then
2698 !trialt(:)=trialt2(:)
2699 !endif
2700 !enddo
2701 !endif
2702 
2703  write(message,'(a)') ' symplanes...'
2704 
2705 !Must use the convention of table 1.3 of the international
2706 !tables for crystallography, see also pp 788 and 789.
2707 !Often, one needs to specialize the selection according
2708 !to the Bravais lattice or the system.
2709 
2710  if(sum(abs(trialt(:)))<nzero .and. iholohedry/=6)then
2711    type_axis=15  ! m
2712    write(label,'(a)') 'a mirror plane'
2713  else if(iholohedry==4 .and. center==0)then    ! primitive tetragonal
2714 
2715    if(directiontype==1)then
2716      type_axis=18  ! primary n
2717      write(label,'(a)') 'a primary n plane'
2718    else if(directiontype==2)then
2719      if(sum(abs(trialt(:)-(/half,zero,zero/)))<nzero .or. &
2720 &     sum(abs(trialt(:)-(/zero,half,zero/)))<nzero       )then
2721        type_axis=16  ! secondary a or b
2722        write(label,'(a)') 'a secondary a or b plane'
2723      else if(sum(abs(trialt(:)-(/zero,zero,half/)))<nzero)then
2724        type_axis=17    ! secondary c
2725        write(label,'(a)') 'a secondary c plane'
2726      else
2727        type_axis=18    ! secondary n
2728        write(label,'(a)') 'a secondary n plane'
2729      end if ! directiontype==2
2730    else if(directiontype==3)then
2731      if( abs(trialt(3))<nzero )then
2732        type_axis=15    ! tertiary m
2733        write(label,'(a)') 'a tertiary m plane'
2734      else if( abs(trialt(3)-half)<nzero )then
2735        type_axis=19    ! tertiary c
2736        write(label,'(a)') 'a tertiary c plane'
2737      end if
2738    end if
2739 
2740  else if(iholohedry==4 .and. center==-1)then    ! inner tetragonal
2741 
2742    if(directiontype==1 .or. directiontype==2)then
2743      if(sum(abs(trialt(:)-(/half,zero,zero/)))<nzero .or. &
2744 &     sum(abs(trialt(:)-(/zero,half,zero/)))<nzero .or. &
2745 &     sum(abs(trialt(:)-(/zero,zero,half/)))<nzero      )then
2746        type_axis=16    ! a, b, or c
2747        write(label,'(a)') 'an a, b or c plane'
2748      else if(sum(abs(trialt(:)-(/half,half,zero/)))<nzero .or. &
2749 &       sum(abs(trialt(:)-(/zero,half,half/)))<nzero .or. &
2750 &       sum(abs(trialt(:)-(/half,zero,half/)))<nzero       )then
2751        type_axis=15    ! n plane, equivalent to m
2752        write(label,'(a)') 'a m plane'
2753      end if ! directiontype==1 or 2
2754    else if(directiontype==3)then
2755      if( abs(trialt(3))<nzero .or. abs(trialt(3)-half)<nzero )then
2756        type_axis=15    ! tertiary c, equivalent to m
2757        write(label,'(a)') 'a tertiary m plane'
2758      else
2759        type_axis=17    ! tertiary d
2760        write(label,'(a)') 'a tertiary d plane'
2761      end if
2762    end if
2763 
2764  else if(iholohedry==5)then    ! hR
2765 
2766    if( abs(sum(abs(trialt(:)))-one) < nzero) then
2767      type_axis=15    ! secondary m
2768      write(label,'(a)') 'a secondary m plane'
2769    else if( abs(sum(abs(trialt(:)))-half) < nzero .or. &
2770 &     abs(sum(abs(trialt(:)))-three*half) < nzero )then
2771      type_axis=16    ! secondary c
2772      write(label,'(a)') 'a secondary c plane'
2773    end if
2774 
2775  else if(iholohedry==6)then    ! hP
2776 
2777    if(directiontype==1)then
2778      if( abs(trialt(3)) <nzero )then
2779        type_axis=15    ! primary m
2780        write(label,'(a)') 'a primary m plane'
2781      end if
2782    else if(directiontype==2)then
2783      if( abs(trialt(3)) <nzero )then
2784        type_axis=15    ! secondary m
2785        write(label,'(a)') 'a secondary m plane'
2786      else if( abs(trialt(3)-half) < nzero ) then
2787        type_axis=16    ! secondary c
2788        write(label,'(a)') 'a secondary c plane'
2789      end if
2790    else if(directiontype==3)then
2791      if( abs(trialt(3)) <nzero )then
2792        type_axis=17    ! tertiary m
2793        write(label,'(a)') 'a tertiary m plane'
2794      else if( abs(trialt(3)-half) < nzero ) then
2795        type_axis=18    ! tertiary c
2796        write(label,'(a)') 'a tertiary c plane'
2797      end if
2798    end if ! directiontype
2799 
2800 !  else if(iholohedry==7 .and. center==0)then    ! cP
2801  else if(iholohedry==7)then    ! cP
2802 
2803    if(directiontype==1)then
2804      if((sum(abs(isymrelconv(:,:)-mirrorx(:,:)))==0 .and.  &
2805 &     sum(abs(two*abs(trialt(:))-(/zero,half,half/)))<nzero   ).or. &
2806 &     (sum(abs(isymrelconv(:,:)-mirrory(:,:)))==0 .and.  &
2807 &     sum(abs(two*abs(trialt(:))-(/half,zero,half/)))<nzero   ).or. &
2808 &     (sum(abs(isymrelconv(:,:)-mirrorz(:,:)))==0 .and.  &
2809 &     sum(abs(two*abs(trialt(:))-(/half,half,zero/)))<nzero   )    ) then
2810        type_axis=17     ! d
2811        write(label,'(a)') 'a d plane'
2812      else
2813        type_axis=18    ! primary n
2814        write(label,'(a)') 'a primary n plane'
2815      end if
2816    else if(directiontype==3)then
2817      if(sum(abs(two*abs(trialt(:))-(/half,half,half/)))<nzero       )then
2818        type_axis=17     ! d
2819        write(label,'(a)') 'a d plane'
2820      else if( abs(sum(abs(trialt(:)))-half) < nzero .or. &
2821 &       abs(sum(abs(trialt(:)))-three*half) < nzero ) then
2822        type_axis=18    ! tertiary n
2823        write(label,'(a)') 'a tertiary n plane'
2824      else if( abs(sum(abs(trialt(:)))-one) < nzero )then
2825        type_axis=15    ! tertiary m
2826        write(label,'(a)') 'a tertiary m plane'
2827      end if
2828    end if
2829 
2830 !  Now, treat all other cases (including other centered Bravais lattices)
2831  else if( sum(abs(trialt(:)-(/half,zero,zero/)))<nzero .or. &
2832 &   sum(abs(trialt(:)-(/zero,half,zero/)))<nzero .or. &
2833 &   sum(abs(trialt(:)-(/zero,zero,half/)))<nzero       )then
2834    type_axis=16     ! a, b or c
2835    write(label,'(a)') 'an a,b, or c plane'
2836  else if( (directiontype==1 .or. directiontype==2) .and. &
2837 &   (sum(abs(trialt(:)-(/half,half,zero/)))<nzero .or. &
2838 &   sum(abs(trialt(:)-(/zero,half,half/)))<nzero .or. &
2839 &   sum(abs(trialt(:)-(/half,zero,half/)))<nzero     ) )then
2840    type_axis=18     ! n
2841    write(label,'(a)') 'an n plane'
2842  else if( directiontype==3 .and. &
2843 &   sum(abs(trialt(:)-(/half,half,half/)))<nzero )then
2844    type_axis=18     ! n
2845    write(label,'(a)') 'an n plane'
2846  else if((sum(abs(isymrelconv(:,:)-mirrorx(:,:)))==0 .and.  &
2847 &   sum(abs(two*abs(trialt(:))-(/zero,half,half/)))<nzero   ).or. &
2848 &   (sum(abs(isymrelconv(:,:)-mirrory(:,:)))==0 .and.  &
2849 &   sum(abs(two*abs(trialt(:))-(/half,zero,half/)))<nzero   ).or. &
2850 &   (sum(abs(isymrelconv(:,:)-mirrorz(:,:)))==0 .and.  &
2851 &   sum(abs(two*abs(trialt(:))-(/half,half,zero/)))<nzero   )    ) then
2852    type_axis=17     ! d
2853    write(label,'(a)') 'a d plane'
2854  else if( directiontype==3 .and. &
2855 &   sum(abs(two*abs(trialt(:))-(/half,half,half/)))<nzero       )then
2856    type_axis=17     ! d
2857    write(label,'(a)') 'a d plane'
2858  else
2859    type_axis=19     ! g (all other planes with
2860 !  unconventional glide vector)
2861    write(label,'(a)') 'a g plane'
2862  end if
2863 
2864  if (verbose) then
2865    write(message,'(a,i3,a,a)')' symplanes : the symmetry operation no. ',isym,' is ', trim(label)
2866    call wrtout(std_out,message,'COLL')
2867  end if
2868 
2869 end subroutine symplanes

m_symtk/symrelrot [ Functions ]

[ Top ] [ m_symtk ] [ Functions ]

NAME

 symrelrot

FUNCTION

 Transform the symmetry matrices symrel expressed in the coordinate system rprimd,
 to symmetry matrices symrel expressed in the new coordinate system rprimd_new

INPUTS

 nsym=number of symmetries
 rprimd(3,3)=dimensional primitive translations for real space (bohr)
 rprimd_new(3,3)=new dimensional primitive translations for real space (bohr)

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/Output
 symrel(3,3,nsym)=symmetry operations in real space in terms
 of primitive translations rprimd at input and rprimd_new at output

PARENTS

      ingeo,m_esymm,symbrav,symlatt,symspgr

CHILDREN

      matr3inv

SOURCE

868 subroutine symrelrot(nsym,rprimd,rprimd_new,symrel,tolsym)
869 
870 
871 !This section has been created automatically by the script Abilint (TD).
872 !Do not modify the following lines by hand.
873 #undef ABI_FUNC
874 #define ABI_FUNC 'symrelrot'
875 !End of the abilint section
876 
877  implicit none
878 
879 !Arguments ------------------------------------
880 !scalars
881  integer,intent(in) :: nsym
882  real(dp),intent(in) :: tolsym
883 !arrays
884  integer,intent(inout) :: symrel(3,3,nsym)
885  real(dp),intent(in) :: rprimd(3,3),rprimd_new(3,3)
886 
887 !Local variables-------------------------------
888 !scalars
889  integer :: ii,isym,jj
890  real(dp) :: val
891  character(len=500) :: message
892 !arrays
893  real(dp) :: coord(3,3),coordinvt(3,3),matr1(3,3),matr2(3,3),rprimd_invt(3,3)
894 
895 !**************************************************************************
896 
897 !Compute the coordinates of rprimd_new in the system defined by rprimd(:,:)
898  call matr3inv(rprimd,rprimd_invt)
899  do ii=1,3
900    coord(:,ii)=rprimd_new(1,ii)*rprimd_invt(1,:)+ &
901 &   rprimd_new(2,ii)*rprimd_invt(2,:)+ &
902 &   rprimd_new(3,ii)*rprimd_invt(3,:)
903  end do
904 
905 !Transform symmetry matrices in the system defined by rprimd_new
906  call matr3inv(coord,coordinvt)
907  do isym=1,nsym
908    do ii=1,3
909      matr1(:,ii)=symrel(:,1,isym)*coord(1,ii)+&
910 &     symrel(:,2,isym)*coord(2,ii)+&
911 &     symrel(:,3,isym)*coord(3,ii)
912    end do
913    do ii=1,3
914      matr2(:,ii)=coordinvt(1,:)*matr1(1,ii)+&
915 &     coordinvt(2,:)*matr1(2,ii)+&
916 &     coordinvt(3,:)*matr1(3,ii)
917    end do
918 
919 !  Check that the new symmetry matrices are made of integers, and store them
920    do ii=1,3
921      do jj=1,3
922        val=matr2(ii,jj)
923 !      Need to allow for twice tolsym, in case of centered Bravais lattices (but do it for all lattices ...)
924        if(abs(val-nint(val))>two*tolsym)then
925          write(message,'(2a,a,i3,a,a,3es14.6,a,a,3es14.6,a,a,3es14.6)')&
926 &         'One of the components of symrel is non-integer,',ch10,&
927 &         '  for isym=',isym,ch10,&
928 &         '  symrel=',matr2(:,1),ch10,&
929 &         '         ',matr2(:,2),ch10,&
930 &         '         ',matr2(:,3)
931          MSG_ERROR_CLASS(message, "TolSymError")
932        end if
933        symrel(ii,jj,isym)=nint(val)
934      end do
935    end do
936 !  End loop isym
937  end do
938 
939 end subroutine symrelrot