TABLE OF CONTENTS
- ABINIT/m_symtk
- m_symtk/chkgrp
- m_symtk/chkorthsy
- m_symtk/chkprimit
- m_symtk/holocell
- m_symtk/littlegroup_q
- m_symtk/mati3det
- m_symtk/mati3inv
- m_symtk/matpointsym
- m_symtk/matr3inv
- m_symtk/print_symmetries
- m_symtk/sg_multable
- m_symtk/smallprim
- m_symtk/symatm
- m_symtk/symaxes
- m_symtk/symcharac
- m_symtk/symchk
- m_symtk/symdet
- m_symtk/symmetrize_rprimd
- m_symtk/symmetrize_tnons
- m_symtk/symmetrize_xred
- m_symtk/symplanes
- m_symtk/symrelrot
ABINIT/m_symtk [ Modules ]
NAME
m_symtk
FUNCTION
Low-level tools related to symmetries
COPYRIGHT
Copyright (C) 1998-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_symtk 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 28 use m_numeric_tools, only : isinteger, wrap2_pmhalf 29 use m_hide_lapack, only : matrginv 30 31 implicit none 32 33 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 !
SOURCE
286 subroutine chkgrp(nsym, symafm, symrel, ierr) 287 288 !Arguments ------------------------------------ 289 !scalars 290 integer,intent(in) :: nsym 291 integer,intent(out) :: ierr 292 !arrays 293 integer,intent(in) :: symafm(nsym),symrel(3,3,nsym) 294 295 !Local variables------------------------------- 296 !scalars 297 integer :: isym, jsym, ksym, symafmchk, testeq, print_warning 298 logical :: found_inv 299 character(len=500) :: msg 300 !arrays 301 integer :: chk(3,3) 302 303 ! ************************************************************************* 304 305 !DEBUG 306 !write(std_out,*)' chkgrp : enter' 307 !write(std_out,*)' isym symrel symafm ' 308 !do isym=1,nsym 309 !write(std_out,'(i3,a,9i3,a,i3)' )isym,' ',symrel(:,:,isym),' ',symafm(isym) 310 !end do 311 !ENDDEBUG 312 313 ierr = 0 314 print_warning = 1 315 316 ! 1) Identity must be the first symmetry. 317 if (ANY(symrel(:,:,1) /= identity_3d .or. symafm(1)/=1 )) then 318 ABI_WARNING("First operation must be the identity operator") 319 ierr = ierr + 1 320 end if 321 322 ! 2) The inverse of each element must belong to the group. 323 do isym=1,nsym 324 call mati3inv(symrel(:,:,isym), chk) 325 chk = transpose(chk) 326 found_inv = .FALSE. 327 do jsym=1,nsym 328 if (all(symrel(:,:,jsym) == chk) .and. (symafm(jsym) * symafm(isym) == 1)) then 329 found_inv = .TRUE.; EXIT 330 end if 331 end do 332 333 if (.not. found_inv) then 334 write(msg,'(a,i0,2a)')& 335 "Cannot find the inverse of symmetry operation ",isym,ch10,& 336 "Input symmetries do not form a group!" 337 ABI_WARNING(msg) 338 ierr = ierr + 1 339 end if 340 end do 341 342 ! Check closure under composition. 343 do isym=1,nsym 344 do jsym=1,nsym 345 346 ! Compute the product of the two symmetries 347 chk = MATMUL(symrel(:,:,jsym), symrel(:,:,isym)) 348 symafmchk = symafm(jsym) * symafm(isym) 349 350 ! Check that product array is one of the original symmetries. 351 do ksym=1,nsym 352 testeq = 1 353 if ( ANY(chk/=symrel(:,:,ksym) )) testeq = 0 354 #if 0 355 ! FIXME this check make v4/t26 and v4/t27 fails. 356 ! The rotational part is in the group but with different magnetic part! 357 if (symafmchk /= symafm(ksym)) testeq=0 358 #endif 359 if (testeq==1) exit ! The test is positive 360 end do 361 362 if (testeq == 0 .and. print_warning == 1) then 363 ! The test is negative 364 write(msg, '(a,2i3,a,9a)' )& 365 'Product of symmetries',isym,jsym,' is not in group.',ch10,& 366 'This indicates that the input symmetry elements',ch10,& 367 'do not possess closure under group composition.',ch10,& 368 'ABINIT might stop with an ERROR after trying to correct and making a few more checks.',ch10,& 369 'Action: check symrel, symafm and possibly atomic positions, and fix them.' 370 ABI_WARNING(msg) 371 ierr = ierr + 1 372 print_warning = 0 373 end if 374 375 end do ! jsym 376 end do ! isym 377 378 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 tolsym=defines the tolerance on the orthogonality, after multiplication by 2.
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 output, return 0 if no problem, -1 if there is a problem,
also, suppresses printing of problem
SOURCE
744 subroutine chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tolsym) 745 746 !Arguments ------------------------------------ 747 !scalars 748 integer,intent(in) :: nsym 749 integer,intent(inout) :: iexit 750 real(dp),intent(in) :: tolsym 751 !arrays 752 integer,intent(in) :: symrel(3,3,nsym) 753 real(dp),intent(in) :: gprimd(3,3),rmet(3,3),rprimd(3,3) 754 755 !Local variables------------------------------- 756 !scalars 757 integer :: ii,isym,jj 758 real(dp) :: residual,rmet2 759 character(len=500) :: msg 760 !arrays 761 real(dp) :: prods(3,3),rmet_sym(3,3),rprimd_sym(3,3) 762 763 ! ************************************************************************* 764 765 !DEBUG 766 !write(std_out,'(a,i3)') ' chkorthsy : enter, iexit= ',iexit 767 !write(std_out,'(a,i3)') ' nsym=',nsym 768 !do isym=1,nsym 769 ! write(std_out,'(9i4)')symrel(:,:,isym) 770 !enddo 771 !write(std_out, '(a)') ' Matrix rprimd :' 772 !do ii=1,3 773 ! write(std_out, '(3es16.8)')rprimd(:,ii) 774 !enddo 775 !write(std_out, '(a)') ' Matrix rmet :' 776 !do ii=1,3 777 ! write(std_out, '(3es16.8)')rmet(:,ii) 778 !enddo 779 !ENDDEBUG 780 781 rmet2=zero 782 do ii=1,3 783 do jj=1,3 784 rmet2=rmet2+rmet(ii,jj)**2 785 end do 786 end do 787 788 !Loop over all symmetry operations 789 do isym=1,nsym 790 791 !DEBUG 792 !write(std_out,'(a,a,i4)') ch10,' Check for isym=',isym 793 !ENDDEBUG 794 795 ! Compute symmetric of primitive vectors under point symmetry operations 796 do ii=1,3 797 rprimd_sym(:,ii)=symrel(1,ii,isym)*rprimd(:,1)+& 798 symrel(2,ii,isym)*rprimd(:,2)+& 799 symrel(3,ii,isym)*rprimd(:,3) 800 end do 801 802 ! If the new lattice is the same as the original one, 803 ! the lengths and angles are preserved 804 do ii=1,3 805 rmet_sym(ii,:)=rprimd_sym(1,ii)*rprimd_sym(1,:)+& 806 rprimd_sym(2,ii)*rprimd_sym(2,:)+& 807 rprimd_sym(3,ii)*rprimd_sym(3,:) 808 end do 809 810 residual=zero 811 do ii=1,3 812 do jj=1,3 813 residual=residual+(rmet_sym(ii,jj)-rmet(ii,jj))**2 814 end do 815 end do 816 817 if(sqrt(residual) > four*tolsym*sqrt(rmet2))then 818 if(iexit==0)then 819 write(std_out, '(a)') ' Matrix rprimd :' 820 do ii=1,3 821 write(std_out, '(3es16.8)')rprimd(:,ii) 822 enddo 823 write(std_out, '(a)') ' Matrix rmet :' 824 do ii=1,3 825 write(std_out, '(3es16.8)')rmet(:,ii) 826 enddo 827 write(std_out, '(a)') ' Matrix rprimd_sym :' 828 do ii=1,3 829 write(std_out, '(3es16.8)')rprimd_sym(:,ii) 830 enddo 831 write(std_out, '(a)') ' Matrix rmet_sym :' 832 do ii=1,3 833 write(std_out, '(3es16.8)')rmet_sym(:,ii) 834 enddo 835 write(std_out, '(a)') ' Matrix rmet_sym-rmet :' 836 do ii=1,3 837 write(std_out, '(3es16.8)')(rmet_sym(:,ii)-rmet(:,ii)) 838 enddo 839 write(msg, '(a,i0,5a,es12.4,a,es12.4,6a)' )& 840 'The symmetry operation number ',isym,' does not preserve',ch10,& 841 'vector lengths and angles.',ch10,& 842 'The value of the square root of residual is: ',sqrt(residual),& 843 & ' that is greater than threshold:', four*tolsym*sqrt(rmet2),ch10,& 844 'Action: modify rprim, acell and/or symrel so that',ch10,& 845 'vector lengths and angles are preserved.',ch10,& 846 'Beware, the tolerance on symmetry operations is very small.' 847 ABI_ERROR(msg) 848 else 849 iexit=-1 850 end if 851 end if 852 853 ! Also, the scalar product of rprimd_sym and gprimd must give integer numbers 854 do ii=1,3 855 prods(ii,:)=rprimd_sym(1,ii)*gprimd(1,:)+ & 856 rprimd_sym(2,ii)*gprimd(2,:)+ & 857 rprimd_sym(3,ii)*gprimd(3,:) 858 end do 859 860 do ii=1,3 861 do jj=1,3 862 residual=prods(ii,jj)-anint(prods(ii,jj)) 863 if(abs(residual)>two*tolsym)then 864 if(iexit==0)then 865 write(msg, '(a,i0,5a,es12.4,a,es12.4,4a)' )& 866 'The symmetry operation number ',isym,' generates',ch10,& 867 'a different lattice.',ch10,& 868 'The value of the residual is: ',residual, 'that is greater than the threshold:', two*tolsym, ch10,& 869 'Action: modify rprim, acell and/or symrel so that',ch10,& 870 'the lattice is preserved.' 871 ABI_ERROR(msg) 872 else 873 iexit=-1 874 end if 875 end if 876 end do 877 end do 878 879 if(iexit==-1) exit 880 end do ! isym 881 882 if(iexit==1)iexit=0 883 884 !DEBUG 885 !write(std_out,'(a)') ' chkorthsy : exit ' 886 !ENDDEBUG 887 888 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 translation(nsym)= (optional) set to 1 if the symetry operation is a pure translation
SOURCE
912 subroutine chkprimit(chkprim, multi, nsym, symafm, symrel, is_translation) 913 914 !Arguments ------------------------------------ 915 !scalars 916 integer,intent(in) :: chkprim,nsym 917 integer,intent(out) :: multi 918 !arrays 919 integer,intent(in) :: symafm(nsym),symrel(3,3,nsym) 920 integer,intent(out),optional :: is_translation(nsym) 921 922 !Local variables------------------------------- 923 !scalars 924 integer :: isym 925 character(len=500) :: msg 926 927 !************************************************************************** 928 929 if(present(is_translation))then 930 is_translation(:)=0 931 endif 932 933 !Loop over each symmetry operation of the Bravais lattice 934 !Find whether it is the identity, or a pure translation, 935 !without change of sign of the spin 936 multi=0 937 do isym=1,nsym 938 if( abs(symrel(1,1,isym)-1)+& 939 & abs(symrel(2,2,isym)-1)+& 940 & abs(symrel(3,3,isym)-1)+& 941 & abs(symrel(1,2,isym))+abs(symrel(2,1,isym))+& 942 & abs(symrel(2,3,isym))+abs(symrel(3,2,isym))+& 943 & abs(symrel(3,1,isym))+abs(symrel(1,3,isym))+& 944 & abs(symafm(isym)-1) == 0 )then 945 multi=multi+1 946 if(present(is_translation))then 947 is_translation(isym)=1 948 endif 949 end if 950 end do 951 952 !Check whether the cell is primitive 953 if(multi>1)then 954 if(chkprim>0)then 955 write(msg,'(a,a,a,i0,a,a,a,a,a,a,a,a,a)')& 956 'According to the symmetry finder, the unit cell is',ch10,& 957 'NOT primitive. The multiplicity is ',multi,' .',ch10,& 958 'The use of non-primitive unit cells is allowed',ch10,& 959 'only when the current chkprim is 0.',ch10,& 960 'Action: either change your unit cell (rprim or angdeg),',ch10,& 961 'or set chkprim to 0.' 962 ABI_ERROR(msg) 963 else if(chkprim==0)then 964 write(msg,'(3a,i0,a,a,a)')& 965 'According to the symmetry finder, the unit cell is',ch10,& 966 'not primitive, with multiplicity= ',multi,'.',ch10,& 967 'This is allowed, as the current chkprim is 0.' 968 ABI_COMMENT(msg) 969 end if 970 end if 971 972 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 (uses its absolute value, since when the multiplicity of the cell is more than one, the sign of iholohedry is changed). 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)
SOURCE
1371 subroutine holocell(cell_base,enforce,foundc,iholohedry,tolsym) 1372 1373 !Arguments ------------------------------------ 1374 !scalars 1375 integer,intent(in) :: enforce,iholohedry 1376 integer,intent(out) :: foundc 1377 real(dp),intent(in) :: tolsym 1378 !arrays 1379 real(dp),intent(inout) :: cell_base(3,3) 1380 1381 !Local variables ------------------------------ 1382 !scalars 1383 integer :: allequal,ii,orth 1384 real(dp):: aa,scprod1 1385 character(len=500) :: msg 1386 !arrays 1387 integer :: ang90(3),equal(3) 1388 real(dp) :: length(3),metric(3,3),norm(3),rbasis(3,3),rconv(3,3),rconv_new(3,3) 1389 real(dp) :: rnormalized(3,3),symmetrized_length(3) 1390 1391 !************************************************************************** 1392 1393 if(abs(iholohedry)<1 .or. abs(iholohedry)>7)then 1394 write(msg, '(a,i0)' )& 1395 & 'Abs(iholohedry) should be between 1 and 7, while iholohedry=',iholohedry 1396 ABI_BUG(msg) 1397 end if 1398 1399 do ii=1,3 1400 metric(:,ii)=cell_base(1,:)*cell_base(1,ii)+& 1401 & cell_base(2,:)*cell_base(2,ii)+& 1402 & cell_base(3,:)*cell_base(3,ii) 1403 end do 1404 1405 !Examine the angles and vector lengths 1406 ang90(:)=0 1407 if(metric(1,2)**2<tolsym**2*metric(1,1)*metric(2,2))ang90(3)=1 1408 if(metric(1,3)**2<tolsym**2*metric(1,1)*metric(3,3))ang90(2)=1 1409 if(metric(2,3)**2<tolsym**2*metric(2,2)*metric(3,3))ang90(1)=1 1410 orth=0 1411 if(ang90(1)==1 .and. ang90(2)==1 .and. ang90(3)==1) orth=1 1412 equal(:)=0 1413 if(abs(metric(1,1)-metric(2,2))<tolsym*half*(metric(1,1)+metric(2,2)))equal(3)=1 1414 if(abs(metric(1,1)-metric(3,3))<tolsym*half*(metric(1,1)+metric(3,3)))equal(2)=1 1415 if(abs(metric(2,2)-metric(3,3))<tolsym*half*(metric(2,2)+metric(3,3)))equal(1)=1 1416 allequal=0 1417 if(equal(1)==1 .and. equal(2)==1 .and. equal(3)==1) allequal=1 1418 1419 !DEBUG 1420 !write(std_out,*)' holocell : enforce, iholohedry=',enforce, iholohedry 1421 !write(std_out,*)' holocell : ang90=',ang90 1422 !write(std_out,*)' holocell : equal=',equal 1423 !!write(std_out,*)' holocell : tolsym=',tolsym 1424 !!write(std_out,*)' holocell : metric(1,1)=',metric(1,1) 1425 !!write(std_out,*)' holocell : metric(1,2)=',metric(1,2) 1426 !ENDDEBUG 1427 1428 foundc=0 1429 if(abs(iholohedry)==1) foundc=1 1430 if(abs(iholohedry)==2 .and. ang90(1)+ang90(3)==2 ) foundc=1 1431 if(abs(iholohedry)==3 .and. orth==1) foundc=1 1432 if(abs(iholohedry)==4 .and. orth==1 .and. & 1433 & (equal(3)==1 .or. equal(2)==1 .or. equal(1)==1) ) foundc=1 1434 if(abs(iholohedry)==5 .and. allequal==1 .and. & 1435 & (abs(metric(1,2)-metric(2,3))<tolsym*metric(2,2)) .and. & 1436 & (abs(metric(1,2)-metric(1,3))<tolsym*metric(1,1)) ) foundc=1 1437 if(abs(iholohedry)==6 .and. equal(3)==1 .and. & 1438 & ang90(1)==1 .and. ang90(2)==1 .and. & 1439 & abs(2*metric(1,2)+metric(1,1))<tolsym*metric(1,1) ) foundc=1 1440 if(abs(iholohedry)==7 .and. orth==1 .and. allequal==1) foundc=1 1441 1442 !DEBUG 1443 !write(std_out, '(a,2i4)' )' holocell : foundc, enforce=',foundc,enforce 1444 !ENDDEBUG 1445 1446 !------------------------------------------------------------------------------------- 1447 !Possibly enforce the holohedry (if it is to be enforced !) 1448 1449 if(foundc==0.and.enforce==1.and.abs(iholohedry)/=1)then 1450 1451 ! Copy the cell_base vectors, and possibly fix the tetragonal axis to be the c-axis 1452 ! XG20201016 WARNING : in principle, one should NOT use the 'equal' information, since precisely this enforcement 1453 ! has the aim to reinstall the symmetries while they are broken !! 1454 if(abs(iholohedry)==4.and.equal(1)==1)then 1455 rconv(:,3)=cell_base(:,1) ; rconv(:,1)=cell_base(:,2) ; rconv(:,2)=cell_base(:,3) 1456 else if (abs(iholohedry)==4.and.equal(2)==1)then 1457 rconv(:,3)=cell_base(:,2) ; rconv(:,2)=cell_base(:,1) ; rconv(:,1)=cell_base(:,3) 1458 else 1459 rconv(:,:)=cell_base(:,:) 1460 end if 1461 1462 ! Compute the length of the three conventional vectors 1463 length(1)=sqrt(sum(rconv(:,1)**2)) 1464 length(2)=sqrt(sum(rconv(:,2)**2)) 1465 length(3)=sqrt(sum(rconv(:,3)**2)) 1466 1467 ! Take care of the first conventional vector aligned with rbasis(:,3) (or aligned with the trigonal axis if rhombohedral) 1468 ! and choice of the first normalized direction 1469 if(abs(iholohedry)==5)then 1470 rbasis(:,3)=third*(rconv(:,1)+rconv(:,2)+rconv(:,3)) 1471 else 1472 rbasis(:,3)=rconv(:,3) 1473 end if 1474 norm(3)=sqrt(sum(rbasis(:,3)**2)) 1475 rnormalized(:,3)=rbasis(:,3)/norm(3) 1476 1477 ! Projection of the first conventional vector perpendicular to rbasis(:,3) 1478 ! and choice of the first normalized direction 1479 scprod1=sum(rnormalized(:,3)*rconv(:,1)) 1480 rbasis(:,1)=rconv(:,1)-rnormalized(:,3)*scprod1 1481 norm(1)=sqrt(sum(rbasis(:,1)**2)) 1482 rnormalized(:,1)=rbasis(:,1)/norm(1) 1483 1484 ! Generation of the second vector, perpendicular to the third and first 1485 rnormalized(1,2)=rnormalized(2,3)*rnormalized(3,1)-rnormalized(3,3)*rnormalized(2,1) 1486 rnormalized(2,2)=rnormalized(3,3)*rnormalized(1,1)-rnormalized(1,3)*rnormalized(3,1) 1487 rnormalized(3,2)=rnormalized(1,3)*rnormalized(2,1)-rnormalized(2,3)*rnormalized(1,1) 1488 1489 ! Compute the vectors of the conventional cell, on the basis of iholohedry 1490 if(abs(iholohedry)==2)then 1491 rconv_new(:,3)=rconv(:,3) 1492 rconv_new(:,1)=rconv(:,1) 1493 rconv_new(:,2)=rnormalized(:,2)*length(2) ! Now, the y axis is perpendicular to the two others, that have not been changed 1494 else if(abs(iholohedry)==3.or.abs(iholohedry)==4.or.abs(iholohedry)==7)then 1495 if(abs(iholohedry)==7)then 1496 symmetrized_length(1:3)=sum(length(:))*third 1497 else if(abs(iholohedry)==4)then 1498 symmetrized_length(3)=length(3) 1499 symmetrized_length(1:2)=half*(length(1)+length(2)) 1500 else if(abs(iholohedry)==3)then 1501 symmetrized_length(:)=length(:) 1502 end if 1503 do ii=1,3 1504 rconv_new(:,ii)=rnormalized(:,ii)*symmetrized_length(ii) 1505 end do 1506 else if(abs(iholohedry)==5)then 1507 ! In the normalized basis, they have coordinates (a,0,c), and (-a/2,+-sqrt(3)/2*a,c) 1508 ! c is known, but a is computed from the knowledge of the average length of the initial vectors 1509 aa=sqrt(sum(length(:)**2)*third-norm(3)**2) 1510 rconv_new(:,1)=aa*rnormalized(:,1)+rbasis(:,3) 1511 rconv_new(:,2)=aa*half*(-rnormalized(:,1)+sqrt(three)*rnormalized(:,2))+rbasis(:,3) 1512 rconv_new(:,3)=aa*half*(-rnormalized(:,1)-sqrt(three)*rnormalized(:,2))+rbasis(:,3) 1513 else if(abs(iholohedry)==6)then 1514 1515 ! In the normalized basis, they have coordinates (a,0,0), (-a/2,+-sqrt(3)/2*a,0), and (0,0,c) 1516 ! c is known, but a is computed from the knowledge of the average length of the initial vectors 1517 aa=half*(length(1)+length(2)) 1518 rconv_new(:,1)=aa*rnormalized(:,1) 1519 rconv_new(:,2)=aa*half*(-rnormalized(:,1)+sqrt(three)*rnormalized(:,2)) 1520 rconv_new(:,3)=rconv(:,3) 1521 end if 1522 1523 ! Copy back the cell_base vectors 1524 if(abs(iholohedry)==4.and.equal(1)==1)then 1525 cell_base(:,3)=rconv_new(:,2) ; cell_base(:,2)=rconv_new(:,1) ; cell_base(:,1)=rconv_new(:,3) 1526 else if (abs(iholohedry)==4.and.equal(2)==1)then 1527 cell_base(:,3)=rconv_new(:,1) ; cell_base(:,1)=rconv_new(:,2) ; cell_base(:,2)=rconv_new(:,3) 1528 else 1529 cell_base(:,:)=rconv_new(:,:) 1530 end if 1531 1532 end if 1533 1534 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, 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.
SOURCE
1123 subroutine littlegroup_q(nsym,qpt,symq,symrec,symafm,timrev,prtvol,use_sym) 1124 1125 !Arguments ------------------------------- 1126 !scalars 1127 integer,intent(in) :: nsym 1128 integer,intent(in),optional :: prtvol,use_sym 1129 integer,intent(out) :: timrev 1130 !arrays 1131 integer,intent(in) :: symrec(3,3,nsym), symafm(nsym) 1132 integer,intent(out) :: symq(4,2,nsym) 1133 real(dp),intent(in) :: qpt(3) 1134 1135 !Local variables ------------------------- 1136 !scalars 1137 integer :: ii,isign,isym,itirev,my_prtvol 1138 real(dp),parameter :: tol=2.d-8 1139 real(dp) :: reduce 1140 character(len=500) :: msg 1141 !arrays 1142 real(dp) :: difq(3),qsym(3),shift(3) 1143 1144 ! ********************************************************************* 1145 1146 my_prtvol=0; if (PRESENT(prtvol)) my_prtvol=prtvol 1147 1148 ! Initialise the array symq 1149 symq = 0 1150 1151 isym = symafm(1) ! just to fool abirules and use symafm for the moment 1152 1153 do isym=1,nsym 1154 ! if (symafm(isym) /= 1) cycle ! skip afm symops 1155 ! TODO: check how much of the afm syms are coded in the rf part of the code. cf 1156 ! test v3 / 12 1157 do itirev=1,2 1158 isign=3-2*itirev ! isign is 1 without time-reversal, -1 with time-reversal 1159 1160 ! Get the symmetric of the vector 1161 do ii=1,3 1162 qsym(ii)=qpt(1)*isign*symrec(ii,1,isym)& 1163 +qpt(2)*isign*symrec(ii,2,isym)& 1164 +qpt(3)*isign*symrec(ii,3,isym) 1165 end do 1166 1167 ! Get the difference between the symmetric and the original vector 1168 symq(4,itirev,isym)=1 1169 do ii=1,3 1170 difq(ii)=qsym(ii)-qpt(ii) 1171 ! Project modulo 1 in the interval ]-1/2,1/2] such that difq = reduce + shift 1172 call wrap2_pmhalf(difq(ii),reduce,shift(ii)) 1173 if(abs(reduce)>tol)symq(4,itirev,isym)=0 1174 end do 1175 1176 ! SP: When prtgkk is asked (GKK matrix element will be output), one has to 1177 ! disable symmetries. There is otherwise a gauge problem with the unperturbed 1178 ! and the perturbed wavefunctions. This leads to a +- 5% increase in computational 1179 ! cost but provide the correct GKKs (i.e. the same as without the use of symmetries.) 1180 1181 if (PRESENT(use_sym)) then 1182 if (use_sym == 0) then 1183 symq(4,itirev,isym)=0 1184 symq(4,itirev,1)=1 1185 end if 1186 end if 1187 1188 ! If the operation succeded, change shift from real(dp) to integer, then exit loop 1189 if(symq(4,itirev,isym)/=0)then 1190 if (my_prtvol>0) then 1191 if(itirev==1)write(msg,'(a,i4,a)')' littlegroup_q : found symmetry',isym,' preserves q ' 1192 if(itirev==2)write(msg,'(a,i4,a)')' littlegroup_q : found symmetry ',isym,' + TimeReversal preserves q ' 1193 call wrtout(std_out,msg) 1194 end if 1195 ! Uses the mathematical function NINT = nearest integer 1196 do ii=1,3 1197 symq(ii,itirev,isym)=nint(shift(ii)) 1198 end do 1199 end if 1200 1201 end do !itirev 1202 end do !isym 1203 1204 ! Test time-reversal symmetry 1205 timrev=1 1206 do ii=1,3 1207 ! Unfortunately, this version does not work yet ... 1208 ! call wrap2_pmhalf(2*qpt(ii),reduce,shift(ii)) 1209 ! if(abs(reduce)>tol)timrev=0 1210 ! So, this is left ... 1211 if(abs(qpt(ii))>tol)timrev=0 1212 end do 1213 1214 if(timrev==1.and.my_prtvol>0)then 1215 write(msg, '(3a)' )& 1216 ' littlegroup_q: able to use time-reversal symmetry. ',ch10,& 1217 ' (except for gamma, not yet able to use time-reversal symmetry)' 1218 call wrtout(std_out,msg) 1219 end if 1220 1221 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
SOURCE
144 subroutine mati3det(mm, det) 145 146 !Arguments ------------------------------------ 147 !arrays 148 integer,intent(in) :: mm(3,3) 149 integer,intent(out) :: det 150 151 ! ************************************************************************* 152 det=mm(1,1)*(mm(2,2) * mm(3,3) - mm(3,2) * mm(2,3)) & 153 + mm(2,1)*(mm(3,2) * mm(1,3) - mm(1,2) * mm(3,3)) & 154 + mm(3,1)*(mm(1,2) * mm(2,3) - mm(2,2) * mm(1,3)) 155 156 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.
SOURCE
85 subroutine mati3inv(mm, mit) 86 87 !Arguments ------------------------------------ 88 !arrays 89 integer,intent(in) :: mm(3,3) 90 integer,intent(out) :: mit(3,3) 91 92 !Local variables------------------------------- 93 !scalars 94 integer :: dd 95 character(len=500) :: msg 96 !arrays 97 integer :: tt(3,3) 98 99 ! ************************************************************************* 100 101 tt(1,1) = mm(2,2) * mm(3,3) - mm(3,2) * mm(2,3) 102 tt(2,1) = mm(3,2) * mm(1,3) - mm(1,2) * mm(3,3) 103 tt(3,1) = mm(1,2) * mm(2,3) - mm(2,2) * mm(1,3) 104 tt(1,2) = mm(3,1) * mm(2,3) - mm(2,1) * mm(3,3) 105 tt(2,2) = mm(1,1) * mm(3,3) - mm(3,1) * mm(1,3) 106 tt(3,2) = mm(2,1) * mm(1,3) - mm(1,1) * mm(2,3) 107 tt(1,3) = mm(2,1) * mm(3,2) - mm(3,1) * mm(2,2) 108 tt(2,3) = mm(3,1) * mm(1,2) - mm(1,1) * mm(3,2) 109 tt(3,3) = mm(1,1) * mm(2,2) - mm(2,1) * mm(1,2) 110 dd = mm(1,1) * tt(1,1) + mm(2,1) * tt(2,1) + mm(3,1) * tt(3,1) 111 112 ! Make sure matrix is not singular 113 if (dd /= 0) then 114 mit(:,:)=tt(:,:)/dd 115 else 116 write(msg, '(2a,2x,9(i0,1x),a)' )'Attempting to invert integer array',ch10,mm,' ==> determinant is zero.' 117 ABI_ERROR(msg) 118 end if 119 120 ! If matrix is orthogonal, determinant must be 1 or -1 121 if (abs(dd) /= 1) then 122 write(msg, '(3a,i0)' )'Absolute value of determinant should be one',ch10,'but determinant= ',dd 123 ABI_ERROR(msg) 124 end if 125 126 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
SOURCE
1248 subroutine matpointsym(iatom,mat3,natom,nsym,rprimd,symrel,tnons,xred) 1249 1250 use m_linalg_interfaces 1251 1252 !Arguments ------------------------------------ 1253 !scalars 1254 integer,intent(in) :: iatom,natom,nsym 1255 !arrays 1256 integer,intent(in) :: symrel(3,3,nsym) 1257 real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym),xred(3,natom) 1258 real(dp),intent(inout) :: mat3(3,3) 1259 1260 !Local variables------------------------------- 1261 !scalars 1262 integer :: cell_index,cell_indexp,ii,isym,nsym_point 1263 real(dp) :: xreddiff 1264 !arrays 1265 integer :: symrel_it(3,3) 1266 real(dp) :: mat3_tri(3,3),mat3_tri_sym(3,3),rprimd_inv(3,3),tmp_mat(3,3) 1267 real(dp) :: xredp(3) 1268 1269 !************************************************************************** 1270 1271 !copy rprimd input and construct inverse 1272 rprimd_inv = rprimd 1273 call matrginv(rprimd_inv,3,3) 1274 1275 !transform input mat3 to triclinic frame with rprimd^{-1} * mat3 * rprimd 1276 call dgemm('N','N',3,3,3,one,rprimd_inv,3,mat3,3,zero,tmp_mat,3) 1277 call dgemm('N','N',3,3,3,one,tmp_mat,3,rprimd,3,zero,mat3_tri,3) 1278 1279 !loop over symmetry elements to obtain symmetrized input matrix 1280 mat3_tri_sym = zero 1281 nsym_point = 0 1282 do isym = 1, nsym 1283 1284 ! skip any nonsymmorphic symmetry elements, want to consider point elements only 1285 if(dot_product(tnons(:,isym),tnons(:,isym))>tol8) cycle 1286 1287 ! for current symmetry element, find transformed reduced coordinates of target atom 1288 ! via xredp = symrel * xred 1289 call dgemv('N',3,3,one,dble(symrel(:,:,isym)),3,xred(:,iatom),1,zero,xredp,1) 1290 1291 1292 ! shift xredp into the same unit cell as xred, for comparison 1293 ! label cells as 0..1:0 1..2:1 2..3:2 and -1..0:-1 -2..-1:-2 and so forth 1294 do ii = 1, 3 1295 1296 cell_index = int(xred(ii,iatom)) 1297 if(xred(ii,iatom) < zero) cell_index = cell_index - 1 1298 cell_indexp = int(xredp(ii)) 1299 if(xredp(ii) < zero) cell_indexp = cell_indexp - 1 1300 1301 do while (cell_indexp < cell_index) 1302 xredp(ii) = xredp(ii)+one 1303 cell_indexp = cell_indexp + 1 1304 end do 1305 do while (cell_indexp > cell_index) 1306 xredp(ii) = xredp(ii)-one 1307 cell_indexp = cell_indexp - 1 1308 end do 1309 1310 end do 1311 1312 ! now compare xredp to xred 1313 xreddiff = dot_product(xredp-xred(:,iatom),xredp-xred(:,iatom)) 1314 1315 if (xreddiff < tol8) then 1316 1317 ! accumulate symrel^{-1}*mat3_tri*symrel into mat3_tri_sym iff xredp = xred + L, 1318 ! where is a lattice vector, so symrel leaves the target atom invariant 1319 1320 ! mati3inv gives the inverse transpose of symrel 1321 call mati3inv(symrel(:,:,isym),symrel_it) 1322 call dgemm('N','N',3,3,3,one,mat3_tri,3,dble(symrel(:,:,isym)),3,zero,tmp_mat,3) 1323 call dgemm('T','N',3,3,3,one,dble(symrel_it),3,tmp_mat,3,one,mat3_tri_sym,3) 1324 nsym_point = nsym_point + 1 1325 end if 1326 1327 end do 1328 1329 !normalize by number of point symmetry operations 1330 mat3_tri_sym = mat3_tri_sym/dble(nsym_point) 1331 1332 !transform mat3_tri_sym to cartesian frame with rprimd * mat3_tri_sym * rprimd^{-1} 1333 1334 call dgemm('N','N',3,3,3,one,mat3_tri_sym,3,rprimd_inv,3,zero,tmp_mat,3) 1335 call dgemm('N','N',3,3,3,one,rprimd,3,tmp_mat,3,zero,mat3,3) 1336 1337 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.
SOURCE
177 subroutine matr3inv(aa, ait) 178 179 !Arguments ------------------------------------ 180 !arrays 181 real(dp),intent(in) :: aa(3,3) 182 real(dp),intent(out) :: ait(3,3) 183 184 !Local variables------------------------------- 185 !scalars 186 real(dp) :: dd,det,t1,t2,t3 187 character(len=500) :: msg 188 189 ! ************************************************************************* 190 191 t1 = aa(2,2) * aa(3,3) - aa(3,2) * aa(2,3) 192 t2 = aa(3,2) * aa(1,3) - aa(1,2) * aa(3,3) 193 t3 = aa(1,2) * aa(2,3) - aa(2,2) * aa(1,3) 194 det = aa(1,1) * t1 + aa(2,1) * t2 + aa(3,1) * t3 195 196 !Make sure matrix is not singular 197 if (abs(det)>tol16) then 198 dd=one/det 199 else 200 write(msg, '(2a,2x,9es16.8,a,a,es16.8,a)' )& 201 'Attempting to invert real(8) 3x3 array',ch10,aa(:,:),ch10,' ==> determinant=',det,' is zero.' 202 ABI_BUG(msg) 203 end if 204 205 ait(1,1) = t1 * dd 206 ait(2,1) = t2 * dd 207 ait(3,1) = t3 * dd 208 ait(1,2) = (aa(3,1)*aa(2,3)-aa(2,1)*aa(3,3)) * dd 209 ait(2,2) = (aa(1,1)*aa(3,3)-aa(3,1)*aa(1,3)) * dd 210 ait(3,2) = (aa(2,1)*aa(1,3)-aa(1,1)*aa(2,3)) * dd 211 ait(1,3) = (aa(2,1)*aa(3,2)-aa(3,1)*aa(2,2)) * dd 212 ait(2,3) = (aa(3,1)*aa(1,2)-aa(1,1)*aa(3,2)) * dd 213 ait(3,3) = (aa(1,1)*aa(2,2)-aa(2,1)*aa(1,2)) * dd 214 215 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
SOURCE
3421 subroutine print_symmetries(nsym, symrel, tnons, symafm, unit, mode_paral) 3422 3423 !Arguments ------------------------------------ 3424 !scalars 3425 integer,intent(in) :: nsym 3426 integer,optional,intent(in) :: unit 3427 character(len=4),optional,intent(in) :: mode_paral 3428 !arrays 3429 integer,intent(in) :: symrel(3,3,nsym),symafm(nsym) 3430 real(dp),intent(in) :: tnons(3,nsym) 3431 3432 !Local variables------------------------------- 3433 integer :: my_unt,isym,isymin,isymend,ii,jj 3434 character(len=500) :: msg 3435 character(len=4) :: my_mode 3436 ! ********************************************************************* 3437 3438 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 3439 my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode=mode_paral 3440 3441 !write(msg,'(2a)')ch10,' Rotations Translations Symafm ' 3442 !do isym=1,nsym 3443 ! write(msg,'(1x,3(3i3,1x),4x,3(f11.7,1x),6x,i2)')symrel(:,:,isym),tnons(:,isym),symafm(isym) 3444 ! call wrtout(my_unt,msg,my_mode) 3445 !end do 3446 3447 write(msg,'(2a)')ch10,' Symmetry operations in real space (Rotation tnons AFM)' 3448 call wrtout(my_unt,msg,my_mode) 3449 do isymin=1,nsym,4 3450 isymend=isymin+3 3451 if (isymend>nsym) isymend=nsym 3452 do ii=1,3 3453 write(msg,'(4(3i3,f11.6,i3,3x))')((symrel(ii,jj,isym),jj=1,3),tnons(ii,isym),symafm(isym),isym=isymin,isymend) 3454 call wrtout(my_unt,msg,my_mode) 3455 end do 3456 write(msg,'(a)')ch10 3457 call wrtout(my_unt,msg,my_mode) 3458 end do 3459 3460 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. Treat reasonably well large set of symmetries, where pure translations are present. The translations are optional. This allows to test symrec.
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) [optional]=Fractional translations. tnons_tol [optional]= tolerance on the match for tnons
OUTPUT
ierr=Status error. A non-zero value signals 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 translation 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 the identity and L a real-space 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.
SOURCE
414 subroutine sg_multable(nsym, symafm, symrel, ierr, & 415 & tnons, tnons_tol, multable, toinv) ! optional 416 417 !Arguments ------------------------------------ 418 !scalars 419 integer,intent(in) :: nsym 420 integer,intent(out) :: ierr 421 real(dp),optional,intent(in) :: tnons_tol 422 !arrays 423 integer,intent(in) :: symafm(nsym),symrel(3,3,nsym) 424 integer,optional,intent(out) :: multable(4,nsym,nsym), toinv(4,nsym) 425 real(dp),optional,intent(in) :: tnons(3,nsym) 426 427 !Local variables------------------------------- 428 !scalars 429 integer :: echo,found,ilist_symrel,nptsymm,prd_symafm,prd_ptsymm,ptsymm1,ptsymm2,ptsymm3 430 integer :: sym1,sym2,sym3 431 !integer :: isym 432 real(dp) :: tnons_tol_ 433 logical :: found_inv,iseq 434 character(len=500) :: msg 435 !arrays 436 integer :: nlist_symrel(48),prd_symrel(3,3),ptmultable(48,48),ptsymrel(3,3,48) 437 integer,allocatable :: ptsymm(:),list_symrel(:,:) 438 real(dp) :: prd_tnons(3) 439 real(dp),allocatable :: tnons_(:,:) 440 441 ! ************************************************************************* 442 443 !DEBUG 444 !write(std_out,*)' m_symtk%sg_multable : enter, nsym= ',nsym 445 !ENDDEBUG 446 447 ierr = 0 448 449 ABI_MALLOC(tnons_,(3,nsym)) 450 if(present(tnons))then 451 tnons_=tnons 452 else 453 tnons_=zero 454 endif 455 if(present(tnons_tol))then 456 tnons_tol_=tnons_tol 457 else 458 tnons_tol_=tol5 459 endif 460 461 ! 1) Identity must be the first symmetry. Do not check if tnons_ == 0 as cell might not be primitive. 462 if (any(symrel(:,:,1) /= identity_3d .or. symafm(1) /= 1)) then 463 ABI_WARNING("First operation must be the identity operator") 464 ierr = ierr + 1 465 end if 466 467 ! 2) The inverse of each element must belong to the group. 468 echo = 1 469 do sym1=1,nsym 470 found_inv = .FALSE. 471 do sym2=1,nsym 472 prd_symrel = matmul(symrel(:,:,sym1), symrel(:,:,sym2)) 473 prd_tnons = tnons_(:,sym1) + matmul(symrel(:,:,sym1), tnons_(:,sym2)) 474 prd_symafm = symafm(sym1)*symafm(sym2) 475 if ( all(prd_symrel == identity_3d) .and. isinteger(prd_tnons, tnons_tol_) .and. prd_symafm == 1 ) then 476 found_inv = .TRUE. 477 if (present(toinv)) then 478 toinv(1, sym1) = sym2; toinv(2:4, sym1) = nint(prd_tnons) 479 end if 480 exit 481 end if 482 end do 483 484 if (.not. found_inv) then 485 if(echo == 1) then 486 write(msg,'(a,i0,2a)')& 487 "Cannot find the inverse of symmetry operation ",sym1,ch10,& 488 "Input symmetries do not form a group " 489 ABI_WARNING(msg) 490 echo = 0 491 endif 492 ierr = ierr + 1 493 exit 494 end if 495 end do 496 497 ! 3) 498 !In order to avoid potential cubic scaling with number of atoms, in exotic cases, with large prefactor, 499 !set up lookup table for the point symmetry part of the symmetry operations. 500 !Still cubic, but with a reduced prefactor. To fully eliminate cubic scaling, should 501 !set up lookup table for the tnons_ as well. 502 503 ABI_MALLOC(list_symrel,(nsym,48)) 504 ABI_MALLOC(ptsymm,(nsym)) 505 506 nlist_symrel(:)=0 507 ! Initialize with the first symmetry operation 508 ptsymrel(1:3,1:3,1)=symrel(:,:,1) 509 ptsymm(1)=1 510 nptsymm=1 511 list_symrel(1,1)=1 512 nlist_symrel(1)=1 513 !If more than one symmetry operation, then loop on the other ones, find whether the ptsymm has already been found, 514 !or create one new item in the list 515 if(nsym/=1)then 516 do sym1=2,nsym 517 found=0 518 do ptsymm2=1,nptsymm 519 if(all(symrel(:,:,list_symrel(1,ptsymm2)) == symrel(:,:,sym1)))then 520 ptsymm(sym1)=ptsymm2 ; found=1 521 nlist_symrel(ptsymm2)=nlist_symrel(ptsymm2)+1 522 list_symrel(nlist_symrel(ptsymm2),ptsymm2)=sym1 523 cycle 524 endif 525 enddo 526 if(found==0)then 527 nptsymm=nptsymm+1 528 !DEBUG 529 ! write(std_out,*)' current value of nptsymm, sym1=',nptsymm, sym1 530 !ENDDEBUG 531 ptsymrel(1:3,1:3,nptsymm)=symrel(:,:,sym1) 532 ptsymm(sym1)=nptsymm 533 nlist_symrel(nptsymm)=1 534 list_symrel(1,nptsymm)=sym1 535 endif 536 enddo 537 endif 538 539 !Check that each point symmetry is associated to the same number of translations 540 if(nptsymm/=1)then 541 do ptsymm1=1,nptsymm 542 if(nlist_symrel(ptsymm1)/=nlist_symrel(1))then 543 write(msg, '(9a)' )& 544 & 'The number of translations (and possibly symafm) associated to the same symrel',ch10,& 545 & 'is not the same for all point symmetries',ch10,& 546 & 'This indicates that the input symmetry elements',ch10,& 547 & 'do not possess closure under group composition.',ch10,& 548 & 'Action: check symrel, symafm and fix them.' 549 ABI_WARNING(msg) 550 echo = 0 551 ierr = ierr + 1 552 if (present(multable)) then 553 multable(1,:,:) = 0; multable(2:4,:,:) = huge(0) 554 end if 555 exit 556 endif 557 enddo 558 endif 559 560 !DEBUG 561 ! write(std_out,*)' final value of nptsymm=',nptsymm 562 !ENDDEBUG 563 564 ! 4) 565 !Check closure under composition and construct multiplication table of ptsymrel 566 echo = 1 567 do ptsymm1=1,nptsymm 568 sym1=list_symrel(1,ptsymm1) 569 do ptsymm2=1,nptsymm 570 sym2=list_symrel(1,ptsymm2) 571 ! Compute the product of the two symmetries. 572 prd_symrel = matmul(symrel(:,:,sym1), symrel(:,:,sym2)) 573 ! Check that product array is one of the original point symmetries. 574 iseq= .false. 575 do ptsymm3=1,nptsymm 576 iseq= all(prd_symrel == symrel(:,:,list_symrel(1,ptsymm3) )) 577 if(iseq)then 578 ptmultable(ptsymm1,ptsymm2) = ptsymm3 579 exit 580 endif 581 end do 582 if (.not. iseq .and. echo == 1) then 583 if (echo == 1)then 584 ! The test is negative 585 prd_symafm = symafm(sym1) * symafm(sym2) 586 prd_tnons = tnons_(:, sym1) + matmul(symrel(:,:,sym1), tnons_(:,sym2)) 587 write(msg, '(a,2(i0,1x),2a,3i3,f11.6,i3,a,2(3i3,f11.6,a),5a)' )& 588 'Product of symmetries:',sym1,sym2,' is not in group.',ch10,& 589 prd_symrel(1,1:3),prd_tnons(1),prd_symafm,ch10,& 590 prd_symrel(2,1:3),prd_tnons(2),ch10,& 591 prd_symrel(3,1:3),prd_tnons(3),ch10,& 592 'This indicates that the input symmetry elements',ch10,& 593 'do not possess closure under group composition.',ch10,& 594 'Action: check symrel, symafm and fix them.' 595 ABI_WARNING(msg) 596 echo = 0 597 endif 598 ierr = ierr + 1 599 if (present(multable)) then 600 multable(1, sym1, sym2) = 0; multable(2:4, sym1, sym2) = huge(0) 601 end if 602 exit 603 end if 604 605 end do ! ptsymm2 606 607 !DEBUG 608 ! write(std_out,*)' ptmultable for ptsymm1=',ptsymm1,' by batch of 16 values ' 609 ! write(std_out,'(16i3)')ptmultable(ptsymm1,1:16) 610 ! write(std_out,'(16i3)')ptmultable(ptsymm1,17:32) 611 ! write(std_out,'(16i3)')ptmultable(ptsymm1,33:48) 612 !ENDDEBUG 613 614 if (echo == 0) exit 615 end do ! ptsymm1 616 617 ! 5) 618 ! Check closure under composition and construct multiplication table. 619 ! However, does this only if the ptgroup has been successfull. 620 if(echo/=0 .and. ierr==0)then 621 do sym1=1,nsym 622 ptsymm1=ptsymm(sym1) 623 do sym2=1,nsym 624 ptsymm2=ptsymm(sym2) 625 626 !The equal number of translations for each point symmetry has been checked earlier. 627 !If the full table is not requested, it is now sufficient to check that 628 !the product of all symmetry operations sym1 with a pure translation (ptsymm=1), or with one of the instances 629 !for each point symmetries is indeed present in the table. 630 !This is done to save CPU time when the number of symmetry operations is bigger than 384. 631 632 if (nsym>384 .and. .not.(present(multable))) then 633 if(ptsymm2/=1 .and. sym2/=list_symrel(1,ptsymm2))then 634 cycle 635 endif 636 end if 637 638 !DEBUG 639 !if(ptsymm2<1 .or. ptsymm2>48)then 640 !write(std_out,*)' sym1,sym2,ptsymm1,ptsymm2=',sym1,sym2,ptsymm1,ptsymm2 641 !endif 642 !ENDDEBUG 643 644 ! Compute the product of the two symmetries. Convention {A,a} {B,b} = {AB, a + Ab} 645 ! prd_symrel = matmul(symrel(:,:,sym1), symrel(:,:,sym2)) 646 prd_ptsymm=ptmultable(ptsymm1,ptsymm2) 647 prd_symrel=ptsymrel(:,:,prd_ptsymm) 648 prd_symafm = symafm(sym1) * symafm(sym2) 649 prd_tnons = tnons_(:, sym1) + matmul(symrel(:,:,sym1), tnons_(:,sym2)) 650 !DEBUG 651 !write(std_out,*)' prd_ptsymm,prdsymrel=',prd_ptsymm,prd_symrel 652 !DEBUG 653 654 ! Check that product array is one of the original symmetries. 655 ! Only explore those symmetries that have a symrel that is the product of the two symrel of sym1 and sym2. 656 iseq = .False. 657 do ilist_symrel=1,nlist_symrel(prd_ptsymm) 658 sym3=list_symrel(ilist_symrel,prd_ptsymm) 659 iseq = isinteger(prd_tnons(1) - tnons_(1,sym3), tnons_tol_) 660 if(iseq)then 661 iseq = isinteger(prd_tnons(2) - tnons_(2,sym3), tnons_tol_) 662 if(iseq)then 663 iseq = isinteger(prd_tnons(3) - tnons_(3,sym3), tnons_tol_) 664 if(iseq)then 665 iseq = (prd_symafm == symafm(sym3)) 666 if(iseq)then 667 ! The test is positive 668 if (present(multable)) then 669 multable(1,sym1,sym2) = sym3; multable(2:4,sym1,sym2) = nint(prd_tnons - tnons_(:,sym3)) 670 end if 671 exit 672 endif 673 endif 674 endif 675 endif 676 end do 677 if (.not. iseq .and. echo == 1) then 678 if (echo == 1)then 679 ! The test is negative 680 write(msg, '(a,2(i0,1x),2a,3i3,f11.6,i3,a,2(3i3,f11.6,a),5a)' )& 681 'Product of symmetries:',sym1,sym2,' is not in group.',ch10,& 682 prd_symrel(1,1:3),prd_tnons(1),prd_symafm,ch10,& 683 prd_symrel(2,1:3),prd_tnons(2),ch10,& 684 prd_symrel(3,1:3),prd_tnons(3),ch10,& 685 'This indicates that the input symmetry elements',ch10,& 686 'do not possess closure under group composition.',ch10,& 687 'Action: check symrel, symafm and fix them.' 688 ABI_WARNING(msg) 689 echo = 0 690 endif 691 ierr = ierr + 1 692 if (present(multable)) then 693 multable(1, sym1, sym2) = 0; multable(2:4, sym1, sym2) = huge(0) 694 end if 695 exit 696 end if 697 end do ! sym2 698 if (echo == 0) exit 699 end do ! sym1 700 else 701 if (present(multable)) then 702 do sym1=1,nsym 703 do sym2=1,nsym 704 multable(1, sym1, sym2) = 0; multable(2:4, sym1, sym2) = huge(0) 705 enddo 706 enddo 707 endif 708 endif 709 710 ABI_FREE(list_symrel) 711 ABI_FREE(ptsymm) 712 ABI_FREE(tnons_) 713 714 !DEBUG 715 !write(std_out,*)' m_symtk%sg_multable : exit ' 716 !ENDDEBUG 717 718 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
SOURCE
3153 subroutine smallprim(metmin,minim,rprimd) 3154 3155 !Arguments ------------------------------------ 3156 !arrays 3157 real(dp),intent(in) :: rprimd(3,3) 3158 real(dp),intent(out) :: metmin(3,3),minim(3,3) 3159 3160 !Local variables------------------------------- 3161 !scalars 3162 integer :: ia,ib,ii,ilong,itrial,minimal 3163 integer :: iiter, maxiter = 100000 3164 real(dp) :: determinant,length2,metsum 3165 character(len=500) :: msg 3166 !arrays 3167 integer :: nvecta(3),nvectb(3) 3168 real(dp) :: rmet(3,3),scprod(3),tmpvect(3) 3169 3170 !************************************************************************** 3171 3172 !call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 3173 rmet = MATMUL(TRANSPOSE(rprimd),rprimd) 3174 3175 nvecta(1)=2 ; nvectb(1)=3 3176 nvecta(2)=1 ; nvectb(2)=3 3177 nvecta(3)=1 ; nvectb(3)=2 3178 3179 minim(:,:)=rprimd(:,:) 3180 metmin(:,:)=rmet(:,:) 3181 3182 !DEBUG 3183 !write(std_out,*)' smallprim : starting values, rprim ' 3184 !write(std_out,'(3f16.8)' )rprimd(:,1) 3185 !write(std_out,'(3f16.8)' )rprimd(:,2) 3186 !write(std_out,'(3f16.8)' )rprimd(:,3) 3187 !write(std_out,*)' smallprim : starting values, rmet ' 3188 !write(std_out,'(3f16.8)' )rmet(:,1) 3189 !write(std_out,'(3f16.8)' )rmet(:,2) 3190 !write(std_out,'(3f16.8)' )rmet(:,3) 3191 !call flush(std_out) 3192 !ENDDEBUG 3193 3194 !Note this loop without index 3195 do iiter = 1, maxiter 3196 3197 ! Will exit if minimal=1 is still valid after a trial 3198 ! to reduce the vectors of each of the three pairs 3199 minimal=1 3200 3201 do itrial=1,3 3202 3203 ia=nvecta(itrial) ; ib=nvectb(itrial) 3204 ! Make sure the scalar product is negative 3205 if(metmin(ia,ib)>tol8)then 3206 minim(:,ia)=-minim(:,ia) 3207 metmin(ia,ib)=-metmin(ia,ib) ; metmin(ib,ia)=-metmin(ib,ia) 3208 metmin(ia,itrial)=-metmin(ia,itrial) 3209 metmin(itrial,ia)=-metmin(itrial,ia) 3210 end if 3211 ! Compute the length of the sum vector 3212 length2=metmin(ia,ia)+2*metmin(ia,ib)+metmin(ib,ib) 3213 ! Replace the first vector by the sum vector if the latter is smaller 3214 if(length2/metmin(ia,ia) < one-tol8)then 3215 minim(:,ia)=minim(:,ia)+minim(:,ib) 3216 metmin(ia,ia)=length2 3217 metmin(ia,ib)=metmin(ia,ib)+metmin(ib,ib) 3218 metmin(ia,itrial)=metmin(ia,itrial)+metmin(ib,itrial) 3219 metmin(ib,ia)=metmin(ia,ib) 3220 metmin(itrial,ia)=metmin(ia,itrial) 3221 minimal=0 3222 ! Replace the second vector by the sum vector if the latter is smaller 3223 else if(length2/metmin(ib,ib) < one-tol8)then 3224 minim(:,ib)=minim(:,ia)+minim(:,ib) 3225 metmin(ib,ib)=length2 3226 metmin(ia,ib)=metmin(ia,ib)+metmin(ia,ia) 3227 metmin(itrial,ib)=metmin(itrial,ib)+metmin(itrial,ia) 3228 metmin(ib,ia)=metmin(ia,ib) 3229 metmin(ib,itrial)=metmin(itrial,ib) 3230 minimal=0 3231 end if 3232 3233 end do 3234 3235 if(minimal==1)exit 3236 end do 3237 3238 if (iiter >= maxiter) then 3239 write(msg,'(a,i0,a)') 'the loop has failed to find a set of minimal vectors in ',maxiter,' iterations.' 3240 ABI_BUG(msg) 3241 end if 3242 3243 !DEBUG 3244 !write(std_out,*)' smallprim : after pair optimization ' 3245 !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' minim =',ch10,minim(:,1),ch10,minim(:,2),ch10,minim(:,3) 3246 !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' metmin =',ch10,metmin(:,1),ch10,metmin(:,2),ch10,metmin(:,3) 3247 !write(std_out,*)' smallprim : will start triplet optimization ',ch10 3248 !call flush(std_out) 3249 !ENDDEBUG 3250 3251 !At this stage, the three vectors have angles between each other that are 3252 !comprised between 90 and 120 degrees. It might still be that minus the vector 3253 !that is the sum of the three vectors is smaller than the longest of these vectors 3254 do iiter = 1, maxiter 3255 3256 ! Will exit if minimal=1 is still valid after a trial 3257 ! to replace the longest of the three vectors by one of the triplet sum of the three vectors, with plus or minus sign 3258 3259 ! Find longest of the three vectors 3260 ilong=1 3261 if( metmin(2,2)/metmin(1,1) > one + tol8 )ilong=2 3262 if( metmin(3,3)/metmin(ilong,ilong) > one + tol8)ilong=3 3263 3264 ! Try combination with all same signs 3265 minimal=1 3266 metsum=sum(metmin(:,:)) 3267 itrial=0 3268 if( metsum/metmin(ilong,ilong) < one - tol8)then 3269 ! Better combination indeed ... 3270 minim(:,ilong)=minim(:,1)+minim(:,2)+minim(:,3) 3271 metmin=MATMUL(TRANSPOSE(minim),minim) 3272 minimal=0 3273 else 3274 ! Try combinations with sign of itrial different from others 3275 metsum=two*(metmin(1,1)+metmin(2,2)+metmin(3,3))-metsum 3276 do itrial=1,3 3277 ia=nvecta(itrial) ; ib=nvectb(itrial) 3278 if( (metsum+four*metmin(ia,ib))/metmin(ilong,ilong) < one - tol8)then 3279 ! Better combination indeed ... 3280 metsum=metsum+four*metmin(ia,ib) 3281 minim(:,ilong)=-minim(:,itrial)+minim(:,ia)+minim(:,ib) 3282 metmin=MATMUL(TRANSPOSE(minim),minim) 3283 minimal=0 3284 exit 3285 endif 3286 enddo 3287 endif 3288 3289 !DEBUG 3290 !write(std_out,*)' smallprim : triplet optimization, iiter,ilong,itrial= ',iiter,ilong,itrial 3291 !write(std_out,*)' smallprim : predict met for the new vector=',metsum 3292 !call flush(std_out) 3293 !ENDDEBUG 3294 3295 ! do itrial=1,3 3296 ! ia=nvecta(itrial) ; ib=nvectb(itrial) 3297 ! if(metmin(ia,ia)/metsum > one + tol8)then 3298 ! minim(:,ia)=-minim(:,1)-minim(:,2)-minim(:,3) 3299 ! metmin(ia,ib)=-sum(metmin(:,ib)) 3300 ! metmin(ia,itrial)=-sum(metmin(:,itrial)) 3301 ! metmin(ia,ia)=metsum 3302 ! metmin(ib,ia)=metmin(ia,ib) 3303 ! metmin(itrial,ia)=metmin(ia,itrial) 3304 ! minimal=0 3305 ! end if 3306 ! end do 3307 3308 !DEBUG 3309 !write(std_out,*)' smallprim : found better primitive vector using triplets, itrial= ',itrial 3310 !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' minim =',ch10,minim(:,1),ch10,minim(:,2),ch10,minim(:,3) 3311 !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' metmin =',ch10,metmin(:,1),ch10,metmin(:,2),ch10,metmin(:,3) 3312 !write(std_out,*)' smallprim : will continue triplet optimization ',ch10 3313 !call flush(std_out) 3314 !ENDDEBUG 3315 3316 if(minimal==1)exit 3317 3318 end do 3319 3320 if (iiter >= maxiter) then 3321 write(msg, '(a,i0,a)') 'the second loop has failed to find a set of minimal vectors in ',maxiter, 'iterations.' 3322 ABI_BUG(msg) 3323 end if 3324 3325 !DEBUG 3326 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3) 3327 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' minim =',minim(:,1),ch10,minim(:,2),ch10,minim(:,3) 3328 !ENDDEBUG 3329 3330 !DEBUG 3331 !Change sign of the third vector if not right-handed basis 3332 !determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+& 3333 !& minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+& 3334 !& minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3)) 3335 !write(std_out,*)' smallprim: determinant=',determinant 3336 !ENDDEBUG 3337 3338 !Choose the first vector 3339 !Compute the scalar product of the three minimal vectors with the first original vector 3340 scprod(:)=zero 3341 do ii=1,3 3342 scprod(:)=scprod(:)+minim(ii,:)*rprimd(ii,1) 3343 end do 3344 !Determine the vector with the maximal absolute overlap 3345 itrial=1 3346 if(abs(scprod(2))>abs(scprod(1))+tol8)itrial=2 3347 if(abs(scprod(3))>abs(scprod(itrial))+tol8)itrial=3 3348 !Switch the vectors if needed 3349 if(itrial/=1)then 3350 tmpvect(:)=minim(:,1) 3351 minim(:,1)=minim(:,itrial) 3352 minim(:,itrial)=tmpvect(:) 3353 end if 3354 !Choose the sign 3355 if(scprod(itrial)<tol8)minim(:,1)=-minim(:,1) 3356 3357 !DEBUG 3358 !Change sign of the third vector if not right-handed basis 3359 !determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+& 3360 !& minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+& 3361 !& minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3)) 3362 !write(std_out,*)' smallprim: determinant=',determinant 3363 !ENDDEBUG 3364 3365 !Choose the second vector 3366 !Compute the scalar product of the second and third minimal vectors with the second original vector 3367 scprod(2:3)=zero 3368 do ii=1,3 3369 scprod(2:3)=scprod(2:3)+minim(ii,2:3)*rprimd(ii,2) 3370 end do 3371 !Determine the vector with the maximal absolute overlap 3372 itrial=2 3373 if(abs(scprod(3))>abs(scprod(2))+tol8)itrial=3 3374 !Switch the vectors if needed 3375 if(itrial/=2)then 3376 tmpvect(:)=minim(:,2) 3377 minim(:,2)=minim(:,itrial) 3378 minim(:,itrial)=tmpvect(:) 3379 end if 3380 !Choose the sign 3381 if(scprod(itrial)<tol8)minim(:,2)=-minim(:,2) 3382 3383 !Change sign of the third vector if not right-handed basis 3384 determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+& 3385 & minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+& 3386 & minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3)) 3387 if(determinant<-tol8)minim(:,3)=-minim(:,3) 3388 if(abs(determinant)<tol8)then 3389 ABI_BUG('minim gives vanishing unit cell volume.') 3390 end if 3391 3392 !Final computation of metmin 3393 do ii=1,3 3394 metmin(ii,:)=minim(1,ii)*minim(1,:)+ minim(2,ii)*minim(2,:)+ minim(3,ii)*minim(3,:) 3395 end do 3396 3397 !DEBUG 3398 !write(std_out,'(2a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',ch10,rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3) 3399 !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' minim =',ch10,minim(:,1),ch10,minim(:,2),ch10,minim(:,3) 3400 !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' metmin =',ch10,metmin(:,1),ch10,metmin(:,2),ch10,metmin(:,3) 3401 !write(std_out,'(a)')' smallprim : exit ' 3402 !call flush(std_out) 3403 !ENDDEBUG 3404 3405 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) is the 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-2024 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 [print_indsym]: Print indsym table to std_out if the number of atoms is smaller that print_indsym Default: -1 i.e. no output is provided.
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.
SOURCE
2179 subroutine symatm(indsym, natom, nsym, symrec, tnons, tolsym, typat, xred, print_indsym) 2180 2181 !Arguments ------------------------------------ 2182 !scalars 2183 integer,intent(in) :: natom,nsym 2184 integer,optional,intent(in) :: print_indsym 2185 real(dp), intent(in) :: tolsym 2186 !arrays 2187 integer,intent(in) :: symrec(3,3,nsym),typat(natom) 2188 integer,intent(out) :: indsym(4,nsym,natom) 2189 real(dp),intent(in) :: tnons(3,nsym),xred(3,natom) 2190 2191 !Local variables------------------------------- 2192 !scalars 2193 integer :: eatom,errout,iatom,ii,isym,mu,print_indsym_ 2194 real(dp) :: difmax,err 2195 character(len=500) :: msg 2196 !arrays 2197 integer :: transl(3) 2198 real(dp) :: difmin(3),tratom(3) 2199 2200 ! ************************************************************************* 2201 2202 !DEBUG 2203 !write(std_out,'(a,i4,es12.4)')' symatm : enter, nsym,tolsym=',nsym,tolsym 2204 !write(std_out,'(a,es12.4)')' symatm : xred=' 2205 !do ii=1,natom 2206 ! write(std_out,'(i4,3es18.10)')ii,xred(1:3,ii) 2207 !enddo 2208 !write(std_out,'(a,es12.4)')' symatm : isym,symrec,tnons=' 2209 !do isym=1,nsym 2210 ! write(std_out,'(i6,9i4,3es18.10)')isym,symrec(:,:,isym),tnons(1:3,isym) 2211 !enddo 2212 !ENDDEBUG 2213 2214 err=zero 2215 errout=0 2216 2217 do isym=1,nsym 2218 do iatom=1,natom 2219 2220 do mu=1,3 ! Apply inverse transformation to original coordinates. Note transpose of symrec. 2221 tratom(mu) = dble(symrec(1,mu,isym))*(xred(1,iatom)-tnons(1,isym))& 2222 & +dble(symrec(2,mu,isym))*(xred(2,iatom)-tnons(2,isym))& 2223 & +dble(symrec(3,mu,isym))*(xred(3,iatom)-tnons(3,isym)) 2224 end do 2225 ! 2226 ! Find symmetrically equivalent atom 2227 call symchk(difmin,eatom,natom,tratom,transl,typat(iatom),typat,xred) 2228 ! 2229 ! Put information into array indsym: translations and label 2230 indsym(1,isym,iatom)=transl(1) 2231 indsym(2,isym,iatom)=transl(2) 2232 indsym(3,isym,iatom)=transl(3) 2233 indsym(4,isym,iatom)=eatom 2234 ! 2235 ! Keep track of maximum difference between transformed coordinates and 2236 ! nearest "target" coordinate 2237 difmax=max(abs(difmin(1)),abs(difmin(2)),abs(difmin(3))) 2238 err=max(err,difmax) 2239 2240 if(errout==3)then 2241 write(msg, '(a)' )& 2242 ' Suppress warning about finding symmetrically equivalent atoms, as mentioned already three times.' 2243 ABI_WARNING(msg) 2244 errout=errout+1 2245 endif 2246 2247 if (difmax>tolsym .and. errout<3) then ! Print warnings if differences exceed tolerance 2248 write(msg, '(3a,i3,a,i6,a,i3,a,a,3f18.12,3a,es12.4)' )& 2249 ' Trouble finding symmetrically equivalent atoms',ch10,& 2250 ' Applying inv of symm number',isym,' to atom number',iatom,' of typat',typat(iatom),ch10,& 2251 ' gives tratom=',tratom(1:3),'.',ch10,& 2252 ' This is further away from every atom in crystal than the allowed tolerance, tolsym=',tolsym 2253 ABI_WARNING(msg) 2254 2255 write(msg, '(a,3i3,a,a,3i3,a,a,3i3)' ) & 2256 ' The inverse symmetry matrix is',symrec(1,1:3,isym),ch10,& 2257 ' ',symrec(2,1:3,isym),ch10,& 2258 ' ',symrec(3,1:3,isym) 2259 call wrtout(std_out,msg) 2260 write(msg, '(a,3f18.12)' )' and the nonsymmorphic transl. tnons =',(tnons(mu,isym),mu=1,3) 2261 2262 call wrtout(std_out,msg) 2263 write(msg, '(a,1p,3es12.4,a,a,i5)' ) & 2264 ' The nearest coordinate differs by',difmin(1:3),ch10,& 2265 ' for indsym(nearest atom)=',indsym(4,isym,iatom) 2266 call wrtout(std_out,msg) 2267 ! 2268 ! Use errout to reduce volume of error diagnostic output 2269 if (errout==0) then 2270 write(msg,'(6a)') ch10,& 2271 ' This indicates that when symatm attempts to find atoms symmetrically',ch10, & 2272 ' related to a given atom, the nearest candidate is further away than some',ch10,& 2273 ' tolerance. Should check atomic coordinates and symmetry group input data.' 2274 call wrtout(std_out,msg) 2275 end if 2276 errout=errout+1 2277 2278 end if !difmax>tol 2279 end do !iatom 2280 end do !isym 2281 2282 ! MG: Do not change this behaviour. symatm is called many times in the EPH code in which we have tons of q-points 2283 ! and it's really annoying to see this output repeated over and over again. 2284 ! If you need to print the indsym table at the beginning of the calculation, find the call to symatm 2285 ! and pass the optional argument print_indsym_ or use `abitk crystal_print FILE --prtvol 1` 2286 print_indsym_ = -1; if (present(print_indsym)) print_indsym_ = print_indsym 2287 if (natom <= print_indsym_) then 2288 do iatom=1,natom 2289 write(msg, '(a,i0,a)' )' symatm: atom number ',iatom,' is reached starting at atom' 2290 call wrtout(std_out,msg) 2291 do ii=1,(nsym-1)/24+1 2292 if(natom<100)then 2293 write(msg, '(1x,24i3)' ) (indsym(4,isym,iatom),isym=1+(ii-1)*24,min(nsym,ii*24)) 2294 else 2295 write(msg, '(1x,24i6)' ) (indsym(4,isym,iatom),isym=1+(ii-1)*24,min(nsym,ii*24)) 2296 end if 2297 call wrtout(std_out,msg) 2298 end do 2299 end do 2300 end if 2301 2302 if (err>tolsym) then 2303 write(msg, '(1x,a,1p,e14.5,a,e12.4)' )'symatm: maximum (delta t)=',err,' is larger than tol=',tolsym 2304 ABI_WARNING(msg) 2305 end if 2306 2307 !Stop execution if error is really big 2308 if (err>0.01d0) then 2309 write(msg,'(5a)')& 2310 'Largest error (above) is so large (0.01) that either input atomic coordinates (xred)',ch10,& 2311 'are wrong or space group symmetry data is wrong.',ch10,& 2312 'Action: correct your input file.' 2313 ABI_ERROR(msg) 2314 end if 2315 2316 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.
SOURCE
2597 subroutine symaxes(center,iholohedry,isym,isymrelconv,label,ordersym,tnons_order,trialt,type_axis) 2598 2599 !Arguments ------------------------------------ 2600 !scalars 2601 integer,intent(in) :: center,iholohedry,isym,ordersym,tnons_order 2602 integer,intent(out) :: type_axis 2603 character(len=128),intent(out) :: label 2604 !arrays 2605 integer,intent(in) :: isymrelconv(3,3) 2606 real(dp),intent(in) :: trialt(3) 2607 2608 !Local variables------------------------------- 2609 !scalars 2610 logical,parameter :: verbose=.FALSE. 2611 character(len=500) :: msg 2612 integer :: direction,directiontype 2613 real(dp),parameter :: nzero=1.0d-6 2614 2615 !************************************************************************** 2616 2617 !write(std_out,*)' symaxes : enter, isym=',isym 2618 !write(std_out,*)' symaxes : iholohedry, ',iholohedry 2619 !write(std_out,*)' symaxes : center, ',center 2620 2621 select case(ordersym) 2622 2623 case(2) ! point symmetry 2 2624 ! Must characterize directiontype for cP, tP, tI, and hP Bravais lattices 2625 directiontype=1 2626 if( iholohedry==4 .or. iholohedry==7) then ! tP or cP Bravais lattices 2627 if(abs(isymrelconv(1,1))+ & 2628 & abs(isymrelconv(2,2))+ & 2629 & abs(isymrelconv(3,3)) ==1) directiontype=3 2630 else if(iholohedry==6)then ! hP Bravais lattice 2631 if(sum(isymrelconv(:,:))/=-1 )directiontype=2 2632 if(sum(isymrelconv(:,:))==0 .or. sum(isymrelconv(:,:))==-3 )& 2633 & directiontype=3 2634 ! directiontype=1 corresponds to a primary axis 2635 ! directiontype=2 corresponds to a tertiary axis 2636 ! directiontype=3 corresponds to a secondary axis 2637 end if 2638 2639 ! DEBUG 2640 ! write(std_out,*)' directiontype=',directiontype 2641 ! write(std_out,'(a,3i6)' )' isymrelconv(1:3)=',isymrelconv(:,1) 2642 ! write(std_out,'(a,3i6)' )' isymrelconv(4:6)=',isymrelconv(:,2) 2643 ! write(std_out,'(a,3i6)' )' isymrelconv(7:9)=',isymrelconv(:,3) 2644 ! write(std_out,'(a,i)' )' tnons_order=',tnons_order 2645 ! ENDDEBUG 2646 2647 ! Now, classify the 2 axes 2648 if(directiontype==2)then 2649 type_axis=4 ! secondary 2 (only in the hP Bravais latt case) 2650 write(label,'(a)') 'a secondary 2-axis ' 2651 2652 else if(directiontype==3 .and. iholohedry==4)then 2653 type_axis=21 ! tertiary 2 2654 write(label,'(a)') 'a tertiary 2-axis ' 2655 else if(directiontype==3 .and. & 2656 & center==0 .and. (iholohedry==6.or.iholohedry==7) )then 2657 type_axis=21 ! tertiary 2 2658 write(label,'(a)') 'a tertiary 2-axis ' 2659 else if(tnons_order==1 .or. (iholohedry==4 .and. center==-1) .or. & 2660 & iholohedry==5)then 2661 type_axis=9 ! 2 2662 write(label,'(a)') 'a 2-axis ' 2663 else 2664 type_axis=20 ! 2_1 2665 write(label,'(a)') 'a 2_1-axis ' 2666 end if 2667 2668 case(3) ! point symmetry 3 2669 if(tnons_order==1)then 2670 type_axis=10 ! 3 2671 write(label,'(a)') 'a 3-axis ' 2672 else if(iholohedry==5 .or. iholohedry==7)then 2673 ! This is a special situation : in the same family of parallel 3-axis, 2674 ! one will have an equal number of 3, 3_1 and 3_2 axes, so that 2675 ! it is non-sense to try to classify one of them. 2676 type_axis=10 ! 3, 3_1 or 3_2, undistinguishable 2677 write(label,'(a)') 'a 3, 3_1 or 3_2 axis ' 2678 else 2679 ! DEBUG 2680 ! write(std_out,*)'isymrelconv=',isymrelconv(:,:) 2681 ! write(std_out,*)'trialt=',trialt(:) 2682 ! ENDDEBUG 2683 ! Must recognize 3_1 or 3_2 2684 if(isymrelconv(1,1)==0)then ! 3+ 2685 if(abs(trialt(3)-third)<nzero)type_axis=22 ! 3_1 2686 if(abs(trialt(3)+third)<nzero)type_axis=23 ! 3_2 2687 else if(isymrelconv(1,1)==-1)then ! 3- 2688 if(abs(trialt(3)-third)<nzero)type_axis=23 ! 3_2 2689 if(abs(trialt(3)+third)<nzero)type_axis=22 ! 3_1 2690 end if 2691 write(label,'(a)') 'a 3_1 or 3_2-axis ' 2692 end if 2693 2694 case(4) ! point symmetry 4 2695 if(tnons_order==1)then 2696 type_axis=12 ! 4 2697 write(label,'(a)') 'a 4-axis ' 2698 else if(tnons_order==2)then 2699 type_axis=25 ! 4_2 2700 write(label,'(a)') 'a 4_2-axis ' 2701 else if(center/=0)then 2702 type_axis=24 ! 4_1 or 4_3 2703 write(label,'(a)') 'a 4_1 or 4_3-axis ' 2704 else 2705 ! DEBUG 2706 ! write(std_out,*)'isymrelconv=',isymrelconv(:,:) 2707 ! write(std_out,*)'trialt=',trialt(:) 2708 ! ENDDEBUG 2709 ! Must recognize 4_1 or 4_3, along the three primary directions 2710 do direction=1,3 2711 if(isymrelconv(direction,direction)==1)then ! 2712 if( (direction==1 .and. isymrelconv(2,3)==-1) .or. & 2713 & (direction==2 .and. isymrelconv(3,1)==-1) .or. & 2714 & (direction==3 .and. isymrelconv(1,2)==-1) )then ! 4+ 2715 if(abs(trialt(direction)-quarter)<nzero)type_axis=24 ! 4_1 2716 if(abs(trialt(direction)+quarter)<nzero)type_axis=26 ! 4_3 2717 else if( (direction==1 .and. isymrelconv(2,3)==1) .or. & 2718 & (direction==2 .and. isymrelconv(3,1)==1) .or. & 2719 & (direction==3 .and. isymrelconv(1,2)==1) )then ! 4- 2720 if(abs(trialt(direction)-quarter)<nzero)type_axis=26 ! 4_3 2721 if(abs(trialt(direction)+quarter)<nzero)type_axis=24 ! 4_1 2722 end if 2723 end if 2724 end do 2725 write(label,'(a)') 'a 4_1 or 4_3-axis ' 2726 end if 2727 2728 case(6) ! point symmetry 6 2729 if(tnons_order==1)then 2730 type_axis=14 ! 6 2731 write(label,'(a)') 'a 6-axis ' 2732 else if(tnons_order==2)then 2733 type_axis=29 ! 6_3 2734 write(label,'(a)') 'a 6_3-axis ' 2735 else if(tnons_order==3)then 2736 !write(std_out,*)'isymrelconv=',isymrelconv(:,:) 2737 !write(std_out,*)'trialt=',trialt(:) 2738 !Must recognize 6_2 or 6_4 2739 if(isymrelconv(1,1)==1)then ! 6+ 2740 if(abs(trialt(3)-third)<nzero)type_axis=28 ! 6_2 2741 if(abs(trialt(3)+third)<nzero)type_axis=30 ! 6_4 2742 else if(isymrelconv(1,1)==0)then ! 6- 2743 if(abs(trialt(3)-third)<nzero)type_axis=30 ! 6_4 2744 if(abs(trialt(3)+third)<nzero)type_axis=28 ! 6_2 2745 end if 2746 write(label,'(a)') 'a 6_2 or 6_4-axis ' 2747 else 2748 !write(std_out,*)'isymrelconv=',isymrelconv(:,:) 2749 !write(std_out,*)'trialt=',trialt(:) 2750 !Must recognize 6_1 or 6_5 2751 if(isymrelconv(1,1)==1)then ! 6+ 2752 if(abs(trialt(3)-sixth)<nzero)type_axis=27 ! 6_1 2753 if(abs(trialt(3)+sixth)<nzero)type_axis=31 ! 6_5 2754 else if(isymrelconv(1,1)==0)then ! 6- 2755 if(abs(trialt(3)-sixth)<nzero)type_axis=31 ! 6_5 2756 if(abs(trialt(3)+sixth)<nzero)type_axis=27 ! 6_1 2757 end if 2758 write(label,'(a)') 'a 6_1 or 6_5-axis ' 2759 end if 2760 2761 end select 2762 2763 if (verbose) then 2764 write(msg,'(a,i3,a,a)')' symaxes : the symmetry operation no. ',isym,' is ', trim(label) 2765 call wrtout(std_out,msg) 2766 end if 2767 2768 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
SOURCE
2341 subroutine symcharac(center, determinant, iholohedry, isym, label, symrel, tnons, type_axis) 2342 2343 !Arguments ------------------------------------ 2344 !scalars 2345 integer, intent(in) :: determinant, center, iholohedry, isym 2346 integer, intent(out) :: type_axis 2347 character(len = 128), intent(out) :: label 2348 !arrays 2349 integer,intent(in) :: symrel(3,3) 2350 real(dp),intent(in) :: tnons(3) 2351 2352 !Local variables------------------------------- 2353 !scalars 2354 logical,parameter :: verbose=.FALSE. 2355 integer :: tnons_order, identified, ii, order, iorder 2356 character(len=500) :: msg 2357 !arrays 2358 integer :: identity(3,3),matrix(3,3),trial(3,3) 2359 real(dp) :: reduced(3),trialt(3) 2360 2361 !************************************************************************** 2362 2363 identity(:,:)=0 2364 identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1 2365 trial(:,:)=identity(:,:) 2366 matrix(:,:)=symrel(:,:) 2367 2368 order=0 2369 do iorder=1,6 2370 trial=matmul(matrix,trial) 2371 if(sum((trial-identity)**2)==0)then 2372 order=iorder 2373 exit 2374 end if 2375 if(sum((trial+identity)**2)==0)then 2376 order=iorder 2377 exit 2378 end if 2379 end do 2380 2381 if(order==0)then 2382 type_axis = -2 2383 return 2384 end if 2385 2386 !Determination of the characteristics of proper symmetries (rotations) 2387 if(determinant==1)then 2388 2389 ! Determine the translation vector associated to the rotations 2390 ! and its order : apply the symmetry operation 2391 ! then analyse the resulting vector. 2392 identified=0 2393 trialt(:)=zero 2394 do ii=1,order 2395 trialt(:)=matmul(symrel(:,:),trialt(:))+tnons(:) 2396 end do 2397 ! Gives the associated translation, with components in the 2398 ! interval [-0.5,0.5] . 2399 reduced(:)=trialt(:)-nint(trialt(:)-tol6) 2400 2401 if(sum(abs(reduced(:)))<tol6)identified=1 2402 if( (center==1 .or. center==-3) .and. & 2403 & sum(abs(reduced(:)-(/zero,half,half/)))<tol6 )identified=2 2404 if( (center==2 .or. center==-3) .and. & 2405 & sum(abs(reduced(:)-(/half,zero,half/)))<tol6 )identified=3 2406 if( (center==3 .or. center==-3) .and. & 2407 & sum(abs(reduced(:)-(/half,half,zero/)))<tol6 )identified=4 2408 if(center==-1.and. sum(abs(reduced(:)-(/half,half,half/)))<tol6 )identified=5 2409 2410 ! If the symmetry operation has not been identified, there is a problem ... 2411 if(identified==0) then 2412 type_axis = -1 2413 return 2414 end if 2415 2416 ! Compute the translation vector associated with one rotation 2417 trialt(:)=trialt(:)/order 2418 trialt(:)=trialt(:)-nint(trialt(:)-tol6) 2419 2420 ! Analyse the resulting vector. 2421 identified=0 2422 do ii=1,order 2423 reduced(:)=ii*trialt(:)-nint(ii*trialt(:)-tol6) 2424 if(sum(abs(reduced(:)))<tol6)identified=1 2425 if( (center==1 .or. center==-3) .and. & 2426 & sum(abs(reduced(:)-(/zero,half,half/)))<tol6 )identified=2 2427 if( (center==2 .or. center==-3) .and. & 2428 & sum(abs(reduced(:)-(/half,zero,half/)))<tol6 )identified=3 2429 if( (center==3 .or. center==-3) .and. & 2430 & sum(abs(reduced(:)-(/half,half,zero/)))<tol6 )identified=4 2431 if(center==-1.and. sum(abs(reduced(:)-(/half,half,half/)))<tol6 )identified=5 2432 2433 if(identified/=0)then 2434 tnons_order=ii 2435 exit 2436 end if 2437 end do ! ii 2438 2439 ! Determinant (here=+1, as we are dealing with proper symmetry operations), 2440 ! order, tnons_order and identified are enough to 2441 ! determine the kind of symmetry operation 2442 2443 select case(order) 2444 case(1) ! point symmetry 1 2445 if(identified==1) then 2446 type_axis=8 ! 1 2447 write(label,'(a)') 'the identity' 2448 else 2449 type_axis=7 ! t 2450 write(label,'(a)') 'a pure translation ' 2451 end if 2452 2453 if (verbose) then 2454 write(msg,'(a,i3,2a)')' symspgr : the symmetry operation no. ',isym,' is ',trim(label) 2455 call wrtout(std_out,msg) 2456 end if 2457 2458 case(2,3,4,6) ! point symmetry 2,3,4,6 - rotations 2459 call symaxes(center,iholohedry,isym,symrel,label,order,tnons_order,trialt,type_axis) 2460 end select 2461 2462 else if (determinant==-1)then 2463 2464 ! Now, take care of the improper symmetry operations. 2465 ! Their treatment is relatively easy, except for the mirror planes 2466 select case(order) 2467 case(1) ! point symmetry 1 2468 type_axis=5 ! -1 2469 write(label,'(a)') 'an inversion' 2470 case(2) ! point symmetry 2 - planes 2471 call symplanes(center,iholohedry,isym,symrel,tnons,label,type_axis) 2472 case(3) ! point symmetry 3 2473 type_axis=3 ! -3 2474 write(label,'(a)') 'a -3 axis ' 2475 case(4) ! point symmetry 1 2476 type_axis=2 ! -4 2477 write(label,'(a)') 'a -4 axis ' 2478 case(6) ! point symmetry 1 2479 type_axis=1 ! -6 2480 write(label,'(a)') 'a -6 axis ' 2481 end select 2482 2483 if (order /= 2 .and. verbose) then 2484 write(msg,'(a,i3,2a)')' symspgr : the symmetry operation no. ',isym,' is ',trim(label) 2485 call wrtout(std_out,msg) 2486 end if 2487 2488 end if ! determinant==1 or -1 2489 2490 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)
SOURCE
2030 subroutine symchk(difmin,eatom,natom,tratom,transl,trtypat,typat,xred) 2031 2032 !Arguments ------------------------------------ 2033 !scalars 2034 integer,intent(in) :: natom,trtypat 2035 integer,intent(out) :: eatom 2036 !arrays 2037 integer,intent(in) :: typat(natom) 2038 integer,intent(out) :: transl(3) 2039 real(dp),intent(in) :: tratom(3),xred(3,natom) 2040 real(dp),intent(out) :: difmin(3) 2041 2042 !Local variables------------------------------- 2043 !scalars 2044 integer :: iatom,jatom,trans1,trans2,trans3 2045 real(dp) :: test,test1,test2,test3,testmn 2046 2047 ! ************************************************************************* 2048 2049 !DEBUG 2050 ! write(std_out,'(a,a,i4,3f18.12)') ch10,' symchk : enter, trtypat,tratom=',trtypat,tratom 2051 !ENDDEBUG 2052 2053 !Start testmn out at large value 2054 testmn=1000000.d0 2055 2056 !Loop through atoms-- 2057 !when types agree, check for agreement after primitive translation 2058 jatom=1 2059 do iatom=1,natom 2060 if (trtypat/=typat(iatom)) cycle 2061 2062 ! Check all three components 2063 test1=tratom(1)-xred(1,iatom) 2064 test2=tratom(2)-xred(2,iatom) 2065 test3=tratom(3)-xred(3,iatom) 2066 ! Find nearest integer part of difference 2067 trans1=nint(test1) 2068 trans2=nint(test2) 2069 trans3=nint(test3) 2070 ! Check whether, after translation, they agree 2071 test1=test1-dble(trans1) 2072 test2=test2-dble(trans2) 2073 test3=test3-dble(trans3) 2074 test=abs(test1)+abs(test2)+abs(test3) 2075 if (test<tol10) then 2076 ! Note that abs() is not taken here 2077 difmin(1)=test1 2078 difmin(2)=test2 2079 difmin(3)=test3 2080 jatom=iatom 2081 transl(1)=trans1 2082 transl(2)=trans2 2083 transl(3)=trans3 2084 ! Break out of loop when agreement is within tolerance 2085 exit 2086 else 2087 ! Keep track of smallest difference if greater than tol10 2088 if (test<testmn) then 2089 testmn=test 2090 ! Note that abs() is not taken here 2091 difmin(1)=test1 2092 difmin(2)=test2 2093 difmin(3)=test3 2094 jatom=iatom 2095 transl(1)=trans1 2096 transl(2)=trans2 2097 transl(3)=trans3 2098 end if 2099 end if 2100 2101 ! End loop over iatom. Note a "cycle" and an "exit" inside the loop 2102 end do 2103 2104 eatom=jatom 2105 2106 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
SOURCE
235 subroutine symdet(determinant, nsym, sym) 236 237 !Arguments ------------------------------------ 238 !scalars 239 integer,intent(in) :: nsym 240 !arrays 241 integer,intent(in) :: sym(3,3,nsym) 242 integer,intent(out) :: determinant(nsym) 243 244 !Local variables------------------------------- 245 !scalars 246 integer :: det,isym 247 character(len=500) :: msg 248 249 ! ************************************************************************* 250 251 do isym=1,nsym 252 call mati3det(sym(:,:,isym),det) 253 determinant(isym)=det 254 if (abs(det)/=1) then 255 write(msg,'(2(a,i0), 5a)')& 256 'Abs(determinant) for symmetry number ',isym,' is ',det,' .',ch10,& 257 'For a legitimate symmetry, abs(determinant) must be 1.',ch10,& 258 'Action: check your symmetry operations (symrel) in input file.' 259 ABI_ERROR(msg) 260 end if 261 end do 262 263 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)
SOURCE
1561 subroutine symmetrize_rprimd(bravais,nsym,rprimd,symrel,tolsym) 1562 1563 !Arguments ------------------------------------ 1564 !scalars 1565 integer,intent(in) :: nsym 1566 real(dp),intent(in) :: tolsym 1567 !arrays 1568 integer,intent(in) :: bravais(11),symrel(3,3,nsym) 1569 real(dp),intent(inout) :: rprimd(3,3) 1570 1571 !Local variables------------------------------- 1572 !scalars 1573 integer :: foundc,iexit,ii,jj 1574 real(dp):: rprimd_maxabs 1575 !character(len=500) :: msg 1576 !arrays 1577 real(dp):: aa(3,3),ait(3,3),cell_base(3,3),gprimd(3,3),rmet(3,3),rprimd_new(3,3) 1578 1579 ! ************************************************************************* 1580 1581 !DEBUG 1582 !write(std_out,'(a)') ' symmetrize_rprimd : enter ' 1583 !ENDDEBUG 1584 1585 !Build the conventional cell basis vectors in cartesian coordinates 1586 aa(:,1)=bravais(3:5) 1587 aa(:,2)=bravais(6:8) 1588 aa(:,3)=bravais(9:11) 1589 !Inverse transpose 1590 call matr3inv(aa,ait) 1591 do ii=1,3 1592 cell_base(:,ii)=ait(ii,1)*rprimd(:,1)+ait(ii,2)*rprimd(:,2)+ait(ii,3)*rprimd(:,3) 1593 end do 1594 1595 !DEBUG 1596 !write(std_out,'(a)') ' before holocell, cell_base =' 1597 !do ii=1,3 1598 ! write(std_out,'(3es16.8)') cell_base(:,ii) 1599 !enddo 1600 !ENDDEBUG 1601 1602 !Enforce the proper holohedry on the conventional cell vectors. 1603 call holocell(cell_base,1,foundc,bravais(1),tolsym) 1604 1605 !DEBUG 1606 !write(std_out,'(a)') ' after holocell, cell_base =' 1607 !do ii=1,3 1608 ! write(std_out,'(3es16.8)') cell_base(:,ii) 1609 !enddo 1610 !ENDDEBUG 1611 1612 !Reconstruct the dimensional primitive vectors 1613 do ii=1,3 1614 rprimd_new(:,ii)=aa(1,ii)*cell_base(:,1)+aa(2,ii)*cell_base(:,2)+aa(3,ii)*cell_base(:,3) 1615 end do 1616 1617 !Suppress meaningless values 1618 rprimd_maxabs=maxval(abs(rprimd_new)) 1619 do ii=1,3 1620 do jj=1,3 1621 if(abs(rprimd(ii,jj))<tol12*rprimd_maxabs)rprimd(ii,jj)=zero 1622 enddo 1623 enddo 1624 1625 rprimd(:,:)=rprimd_new(:,:) 1626 1627 !Check whether the symmetry operations are consistent with the lattice vectors 1628 rmet = MATMUL(TRANSPOSE(rprimd), rprimd) 1629 call matr3inv(rprimd, gprimd) 1630 !call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1631 iexit=0 1632 1633 call chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tolsym) 1634 1635 !DEBUG 1636 !write(std_out,'(a)') ' symmetrize_rprimd : exit ' 1637 !ENDDEBUG 1638 1639 end subroutine symmetrize_rprimd
m_symtk/symmetrize_tnons [ Functions ]
[ Top ] [ m_symtk ] [ Functions ]
NAME
symmetrize_tnons
FUNCTION
Given the order of a symmetry operation, make sure that tnons is such that applying "order" times the symmetry operation generate the unity operation, accurately.
INPUTS
nsym=actual number of symmetries symrel(3,3,1:nsym)=symmetry operations in real space in terms of primitive translations tolsym=tolerance that was used to determine the symmetry operations
SIDE EFFECTS
tnons(3,1:nsym)= non-symmorphic translation vectors
SOURCE
1661 subroutine symmetrize_tnons(nsym,symrel,tnons,tolsym) 1662 1663 !Arguments ------------------------------------ 1664 !scalars 1665 integer,intent(in) :: nsym 1666 real(dp),intent(in) :: tolsym 1667 !arrays 1668 integer,intent(in) :: symrel(3,3,nsym) 1669 real(dp),intent(inout) :: tnons(3,nsym) 1670 1671 !Local variables------------------------------- 1672 !scalars 1673 integer :: iorder,isym,order 1674 !character(len=500) :: msg 1675 !arrays 1676 integer :: symrel_mult(3,3) 1677 integer :: unitmat(3,3) 1678 real(dp):: tnons_mult(3) 1679 1680 ! ************************************************************************* 1681 1682 !DEBUG 1683 !write(std_out,'(a)') ' symmetrize_tnons : enter ' 1684 !ENDDEBUG 1685 1686 unitmat=0 1687 unitmat(1,1)=1 ; unitmat(2,2)=1 ; unitmat(3,3)=1 1688 1689 do isym=1,nsym 1690 1691 !DEBUG 1692 !write(std_out,'(a,i4,9i3,3es16.6)') ' isym,symrel,tnons=',isym,symrel(:,:,isym),tnons(:,isym) 1693 !ENDDEBUG 1694 1695 symrel_mult(:,:)=symrel(:,:,isym) ; tnons_mult(:)=tnons(:,isym) 1696 order=0 1697 !Determine the order of the operation 1698 do iorder=1,48 1699 symrel_mult(:,:)=matmul(symrel(:,:,isym),symrel_mult(:,:)) 1700 tnons_mult(:)=matmul(symrel(:,:,isym),tnons_mult(:))+tnons(:,isym) 1701 if(sum(abs(symrel_mult-unitmat))==0)then 1702 if(abs(tnons_mult(1)-nint(tnons_mult(1)))<tolsym*iorder .and. & 1703 & abs(tnons_mult(2)-nint(tnons_mult(2)))<tolsym*iorder .and. & 1704 & abs(tnons_mult(3)-nint(tnons_mult(3)))<tolsym*iorder)then 1705 !The order has been found 1706 order=iorder+1 1707 !Now, adjust the tnons vector, in order to obtain the exact identity 1708 !operation at order "order" 1709 tnons_mult(:)=(tnons_mult(:)-nint(tnons_mult(:)))/(dble(order)) 1710 if(abs(tnons_mult(1))>1.00001e-8) tnons(1,isym)=tnons(1,isym)-tnons_mult(1) 1711 if(abs(tnons_mult(2))>1.00001e-8) tnons(2,isym)=tnons(2,isym)-tnons_mult(2) 1712 if(abs(tnons_mult(3))>1.00001e-8) tnons(3,isym)=tnons(3,isym)-tnons_mult(3) 1713 exit 1714 endif 1715 endif 1716 1717 enddo ! iorder 1718 1719 if(order==0)then 1720 ABI_BUG("Was unable to find order of operation") 1721 endif 1722 enddo 1723 1724 !DEBUG 1725 !write(std_out,'(a)') ' symmetrize_tnons : exit ' 1726 !ENDDEBUG 1727 1728 end subroutine symmetrize_tnons
m_symtk/symmetrize_xred [ Functions ]
[ Top ] [ m_symtk ] [ Functions ]
NAME
symmetrize_xred
FUNCTION
Symmetrize atomic coordinates. Two tasks can be executed : A. If optional argument indsym is present. Using input symmetry matrices symrel which are expressed in terms of the basis of real space primitive translations (array elements are integers), use indsym to make all corresponding atoms coordinate fullfil exactly symmetry operations. 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). B. If optional argument tolsym AND tnons_new are defined. Might also adjust xred in order for tnons to be aligned with the FFT grids. This will deliver new tnons_new. NOTE : Actually, should make two separate routines !
INPUTS
indsym(4,nsym,natom)=(optional) 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 tolsym=(optional) tolerance on symmetries. When defined, one will try to align the symmetry operations with the FFT grid, if the modification is less than tolsym. Take tolsym equal to 1 to deliver possibly large changes of xred, giving suggestions of xred modifications, to be proposed to users.
OUTPUT
tnons_new(3,nsym)=(optional)nonsymmorphic translations for symmetries
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 fixed_mismatch=(optional) At input, needs to be present for tnons_new to be computed At output : 1 if there is a mismatch and this mismatch has been fixed, 0 otherwise mismatch_fft_tnons=(optional) At input, needs to be present for tnons_new to be computed Atd output : non-zero if there is a mismatch between the fft grid and the tnons, gives the number of the first symmetry operation for which there is such a mismatch. Zero otherwise.
SOURCE
1783 subroutine symmetrize_xred(natom,nsym,symrel,tnons,xred,fixed_mismatch,indsym,mismatch_fft_tnons,tnons_new,tolsym) 1784 1785 !Arguments ------------------------------------ 1786 !scalars 1787 integer,intent(in) :: natom,nsym 1788 integer,intent(out),optional :: fixed_mismatch,mismatch_fft_tnons 1789 !arrays 1790 integer,intent(in),optional :: indsym(4,nsym,natom) 1791 integer,intent(in) :: symrel(3,3,nsym) 1792 real(dp),intent(in) :: tnons(3,nsym) 1793 real(dp),intent(in),optional :: tolsym 1794 real(dp),intent(out),optional :: tnons_new(3,nsym) 1795 real(dp),intent(inout) :: xred(3,natom) 1796 1797 !Local variables------------------------------- 1798 !scalars 1799 integer :: iatom,ib,ii,info,irank,isym,isym2 1800 integer :: jj,mismatch_fft_tnons_current 1801 real(dp) :: diff 1802 logical :: dissimilar 1803 !arrays 1804 real(dp) :: delta(4),fc(3),mat(3,3),mult(4)=(/eight,nine,ten,three*four/) 1805 real(dp) :: sgval(3),tsum(3),tt(3),work(15),xredshift(3,1) 1806 real(dp),allocatable :: xredsym(:,:) 1807 real(dp) :: transl(3) ! translation vector 1808 1809 ! ************************************************************************* 1810 ! 1811 !Check whether group contains more than identity; 1812 !if not then simply return after possible copying 1813 if(present(tnons_new))then 1814 tnons_new(:,1:nsym)=tnons(:,1:nsym) 1815 endif 1816 if(present(fixed_mismatch))fixed_mismatch=0 1817 if(present(mismatch_fft_tnons))mismatch_fft_tnons=0 1818 1819 if (nsym>1) then 1820 1821 !DEBUG 1822 ! write(std_out,*) 1823 ! write(std_out,'(a,i4)') 'symmetrize_xred: enter, nsym=',nsym 1824 ! do iatom=1,natom 1825 ! write(std_out,'(a,i4,3es16.6)') 'iatom,xred=',iatom,xred(:,iatom) 1826 ! enddo 1827 ! do isym=1,nsym 1828 ! write(std_out,'(a,i4,9i3,3es16.6)') 'isym,symrel,tnons',isym,symrel(:,:,isym),tnons(:,isym) 1829 ! enddo 1830 ! write(std_out,*)' present(tnons_new),present(tolsym)=',present(tnons_new),present(tolsym) 1831 ! write(std_out,*) 1832 !ENDDEBUG 1833 1834 ABI_MALLOC(xredsym,(3,natom)) 1835 xredsym(:,:)=xred(:,1:natom) 1836 1837 if(present(indsym))then 1838 1839 ! Loop over atoms to determine new, symmetrized positions. 1840 do iatom=1,natom 1841 tsum(:)=0.0d0 1842 ! 1843 ! Loop over symmetries 1844 do isym=1,nsym 1845 ! atom ib is atom into which iatom is rotated by inverse of 1846 ! symmetry isym (inverse of symrel(mu,nu,isym)) 1847 ib=indsym(4,isym,iatom) 1848 ! Find the reduced coordinates after translation=t(indsym)+transl 1849 fc(:)=xred(:,ib)+dble(indsym(1:3,isym,iatom)) 1850 ! Compute [S * (x(indsym)+transl) ] + tnonsymmorphic 1851 tt(:)=dble(symrel(:,1,isym))*fc(1)+& 1852 & dble(symrel(:,2,isym))*fc(2)+& 1853 & dble(symrel(:,3,isym))*fc(3)+ tnons(:,isym) 1854 1855 ! Average over nominally equivalent atomic positions 1856 tsum(:)=tsum(:)+tt(:) 1857 end do ! isym 1858 ! 1859 ! Set symmetrized result to sum over number of terms 1860 xredsym(:,iatom)=tsum(:)/dble(nsym) 1861 1862 end do ! iatom 1863 endif ! present(indsym) 1864 1865 !DEBUG 1866 !do iatom=1,natom 1867 ! write(std_out,'(a,i4,3es20.10)') 'iatom,xredsym=',iatom,xredsym(:,iatom) 1868 !enddo 1869 !ENDDEBUG 1870 1871 ! Loop over symmetry operations to determine possibly new tnons, as well as symmetrized positions. 1872 if(present(tolsym) .and. present(tnons_new) .and. present(fixed_mismatch) .and. present(mismatch_fft_tnons))then 1873 fixed_mismatch=0 1874 mismatch_fft_tnons=0 1875 !The use of tolsym here is only to favor 0.5 over -0.5 1876 tnons_new(:,:)=tnons(:,:)-nint(tnons(:,:)-tolsym) 1877 do isym=1,nsym 1878 mismatch_fft_tnons_current=0 1879 do ii=1,3 1880 delta(:)=tnons(ii,isym)*mult(:) 1881 delta(:)=delta(:)-nint(delta(:)) 1882 ! Is there is a mismatch between FFT and isym for all multipliers ? 1883 if( all(abs(delta(:))>tol8*mult(:)) ) mismatch_fft_tnons_current=1 1884 enddo 1885 !Declare the first symmetry operation that induces a problem 1886 if(mismatch_fft_tnons_current>0 .and. (mismatch_fft_tnons==0)) mismatch_fft_tnons=isym 1887 1888 ! However, also try to propose a solution. 1889 if(mismatch_fft_tnons_current==1)then 1890 ! Compute the pseudo-inverse of symrel-1, then multiply tnons 1891 mat(:,:)=zero; mat(1,1)=one; mat(2,2)=one; mat(3,3)=one 1892 ! This is symrel-1 1893 mat(:,:)=symrel(:,:,isym)-mat(:,:) 1894 do ii=1,3 1895 !Select the smallest modification tnons 1896 delta(:)=tnons(ii,isym)*mult(:) 1897 delta(:)=(delta(:)-nint(delta(:)))/mult(:) 1898 xredshift(ii,1)=delta(1) 1899 do jj=2,4 1900 if(abs(delta(jj))<abs(xredshift(ii,1)))xredshift(ii,1)=delta(jj) 1901 enddo 1902 enddo 1903 call dgelss(3,3,1,mat,3,xredshift(:,1),3,sgval,tol5,irank,work,15,info) 1904 1905 ! xredshift(:,1) is now the tentative shift, to be tested for all symmetries 1906 if( all(abs(xredshift(:,1))<tolsym) )then 1907 fixed_mismatch=1 1908 do isym2=1, nsym 1909 tnons_new(:,isym2)=tnons(:,isym2)+xredshift(:,1)-matmul(symrel(:,:,isym2),xredshift(:,1)) 1910 do ii=1,3 1911 !tnons might now be slighly non-zero. Set to zero such values 1912 if(abs(tnons_new(ii,isym2))<tol6**2)tnons_new(ii,isym2)=zero 1913 delta(:)=tnons_new(ii,isym2)*mult(:) 1914 delta(:)=delta(:)-nint(delta(:)) 1915 ! Is the mismatch between FFT and symmetries still present for all the multipliers ?? 1916 if( all(abs(delta(:))>tol8*mult(:)) ) fixed_mismatch=0 1917 enddo 1918 if(fixed_mismatch==0)exit 1919 enddo 1920 endif 1921 if(fixed_mismatch==1)exit 1922 endif ! mismatch_fft_tnons_current==1 1923 end do ! isym 1924 if(mismatch_fft_tnons/=0)then 1925 if(fixed_mismatch==1)then 1926 do iatom=1,natom 1927 xredsym(:,iatom)=xredsym(:,iatom)+xredshift(:,1) 1928 enddo 1929 endif 1930 endif 1931 1932 !DEBUG 1933 ! write(std_out,*) ' ' 1934 ! write(std_out,*) ' mismatch_fft_tnons, fixed_mismatch=',mismatch_fft_tnons, fixed_mismatch 1935 ! do iatom=1,natom 1936 ! write(std_out,'(a,i4,3es20.10)') ' iatom,xredsym=',iatom,xredsym(:,iatom) 1937 ! enddo 1938 !ENDDEBUG 1939 1940 endif ! present(tolsym) .and. present(tnons_new) 1941 1942 ! -------------------------------------------------------------- 1943 ! Will update the atomic positions only if it is worth to do so. 1944 1945 transl(:)=xredsym(:,1)-nint(xredsym(:,1)) 1946 1947 ! Compute the smallest translation to an integer 1948 do jj=2,natom 1949 do ii=1,3 1950 diff=xredsym(ii,jj)-nint(xredsym(ii,jj)) 1951 if (diff<transl(ii)) transl(ii)=diff 1952 end do 1953 end do 1954 1955 ! Test if the translation on each direction is small 1956 ! Tolerance 1E-13 1957 do ii=1,3 1958 if (abs(transl(ii))>1e-13) transl(ii)=0.0 1959 end do 1960 1961 ! Execute translation 1962 do jj=1,natom 1963 do ii=1,3 1964 xredsym(ii,jj)=xredsym(ii,jj)-transl(ii) 1965 end do 1966 end do 1967 1968 ! Test if xredsym is too similar to xred 1969 ! Tolerance 1E-15 1970 dissimilar=.FALSE. 1971 do jj=1,natom 1972 do ii=1,3 1973 if (abs(xredsym(ii,jj)-xred(ii,jj))>1E-15) dissimilar=.TRUE. 1974 end do 1975 end do 1976 1977 if (dissimilar) xred(:,:)=xredsym(:,:) 1978 ABI_FREE(xredsym) 1979 1980 ! End condition of nsym/=1 1981 end if 1982 1983 !DEBUG 1984 !write(std_out,*) 1985 !write(std_out,'(a)') 'symmetrize_xred : exit' 1986 !do iatom=1,natom 1987 ! write(std_out,'(a,i4,3es20.10)') 'iatom,xred=',iatom,xred(:,iatom) 1988 !enddo 1989 !if(present(tnons_new))then 1990 ! do isym=1,nsym 1991 ! write(std_out,'(a,i4,9i3,3es20.10)') 'isym,symrel,tnons_new',isym,symrel(:,:,isym),tnons_new(:,isym) 1992 ! enddo 1993 !endif 1994 !write(std_out,*) 1995 !ENDDEBUG 1996 1997 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.
SOURCE
2849 subroutine symplanes(center,iholohedry,isym,isymrelconv,itnonsconv,label,type_axis) 2850 2851 !Arguments ------------------------------------ 2852 !scalars 2853 integer,intent(in) :: center,iholohedry,isym 2854 integer,intent(out) :: type_axis 2855 character(len = 128), intent(out) :: label 2856 !arrays 2857 integer,intent(in) :: isymrelconv(3,3) 2858 real(dp),intent(in) :: itnonsconv(3) 2859 2860 !Local variables------------------------------- 2861 !scalars 2862 logical,parameter :: verbose=.FALSE. 2863 character(len=500) :: msg 2864 integer :: directiontype,sum_elements 2865 real(dp),parameter :: nzero=1.0d-6 2866 !arrays 2867 integer :: identity(3,3),mirrormxy(3,3),mirrormyz(3,3),mirrormzx(3,3) 2868 integer :: mirrorx(3,3),mirrorxy(3,3),mirrory(3,3),mirroryz(3,3),mirrorz(3,3) 2869 integer :: mirrorzx(3,3) 2870 real(dp) :: trialt(3) 2871 ! real(dp) :: itnonsconv2(3),trialt2(3) 2872 2873 !************************************************************************** 2874 2875 !write(std_out,*)' symplanes : enter' 2876 !write(std_out,*)' center,iholohedry,isym,isymrelconv,itnonsconv=',center,iholohedry,isym,isymrelconv,itnonsconv 2877 2878 identity(:,:)=0 2879 identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1 2880 2881 !Will be a mirror plane, but one must characterize 2882 !(1) the type of plane (primary, secondary or tertiary) 2883 !(2) the gliding vector. One now defines a few matrices. 2884 mirrorx(:,:)=identity(:,:) ; mirrorx(1,1)=-1 2885 mirrory(:,:)=identity(:,:) ; mirrory(2,2)=-1 2886 mirrorz(:,:)=identity(:,:) ; mirrorz(3,3)=-1 2887 mirrorxy(:,:)=0 ; mirrorxy(1,2)=1 ; mirrorxy(2,1)=1 ; mirrorxy(3,3)=1 2888 mirrorzx(:,:)=0 ; mirrorzx(1,3)=1 ; mirrorzx(3,1)=1 ; mirrorzx(2,2)=1 2889 mirroryz(:,:)=0 ; mirroryz(2,3)=1 ; mirroryz(3,2)=1 ; mirroryz(1,1)=1 2890 mirrormxy(:,:)=0 ; mirrormxy(1,2)=-1 ; mirrormxy(2,1)=-1 ; mirrormxy(3,3)=1 2891 mirrormzx(:,:)=0 ; mirrormzx(1,3)=-1 ; mirrormzx(3,1)=-1 ; mirrormzx(2,2)=1 2892 mirrormyz(:,:)=0 ; mirrormyz(2,3)=-1 ; mirrormyz(3,2)=-1 ; mirrormyz(1,1)=1 2893 2894 !Determine the type of plane. At the end, 2895 !directiontype=1 will correspond to a primary axis (or equivalent 2896 !axes for orthorhombic) 2897 !directiontype=2 will correspond to a secondary axis 2898 !directiontype=3 will correspond to a tertiary axis 2899 !See table 2.4.1, 11.2 and 11.3 of the international tables for crystallography 2900 directiontype=0 2901 !The sum of elements of the matrices allow to characterize them 2902 sum_elements=sum(isymrelconv(:,:)) 2903 2904 if(sum_elements==1)then 2905 ! The mirror plane perpendicular to the c axis is always primary 2906 if( sum(abs(isymrelconv(:,:)-mirrorz(:,:)))==0 )then 2907 directiontype=1 2908 ! All the other planes with a symrel matrix whose sum of elements is 1 2909 ! are a or b planes. They are primary or 2910 ! secondary planes, depending the holohedry. 2911 else if(sum(isymrelconv(:,:))==1)then 2912 if( iholohedry==2 .or. iholohedry==3 .or. iholohedry==7 )then 2913 directiontype=1 2914 else if(iholohedry==4 .or. iholohedry==6)then 2915 directiontype=2 2916 end if 2917 end if 2918 end if 2919 2920 !All the planes with a symrel matrix whose sum of elements 2921 !is 2 are secondary planes (table 11.3). 2922 if( sum_elements==2 ) directiontype=2 2923 2924 !The planes with a symrel matrix whose sum of elements 2925 !is 3 or 0 are tertiary planes 2926 if( sum_elements==3 .or. sum_elements==0 )directiontype=3 2927 2928 !One is left with sum_elements=-1, tertiary for tetragonal 2929 !or cubic, secondary for hexagonal 2930 if( sum_elements==-1)then 2931 if(iholohedry==4 .or. iholohedry==7)directiontype=3 2932 if(iholohedry==6)directiontype=2 2933 end if 2934 2935 2936 !Now, determine the gliding vector 2937 !First, apply the symmetry operation 2938 !to itnonsconv, in order to get the translation vector 2939 !under the application of twice the symmetry operation 2940 trialt(:)=matmul(isymrelconv(:,:),itnonsconv(:)) +itnonsconv(:) 2941 !Get the translation associated with one application, 2942 !and force its components to be in the interval ]-0.5,0.5] . 2943 trialt(:)=trialt(:)*half 2944 trialt(:)=trialt(:)-nint(trialt(:)-nzero) 2945 2946 !If there is a glide vector for the initial choice of itnonsconv, 2947 !it might be that it disappears if itnonsconv is translated by a 2948 !lattice vector of the conventional cell 2949 !if(trialt(1)**2+trialt(2)**2+trialt(3)**2>tol5)then 2950 !do ii=1,3 2951 !itnonsconv2(:)=itnonsconv(:) 2952 !itnonsconv2(ii)=itnonsconv(ii)+one 2953 !trialt2(:)=matmul(isymrelconv(:,:),itnonsconv2(:)) +itnonsconv2(:) 2954 !trialt2(:)=trialt2(:)*half 2955 !trialt2(:)=trialt2(:)-nint(trialt2(:)-nzero) 2956 !if(trialt2(1)**2+trialt2(2)**2+trialt2(3)**2<tol5)then 2957 !trialt(:)=trialt2(:) 2958 !endif 2959 !enddo 2960 !endif 2961 2962 write(msg,'(a)') ' symplanes...' 2963 2964 !Must use the convention of table 1.3 of the international 2965 !tables for crystallography, see also pp 788 and 789. 2966 !Often, one needs to specialize the selection according 2967 !to the Bravais lattice or the system. 2968 2969 if(sum(abs(trialt(:)))<nzero .and. iholohedry/=6)then 2970 type_axis=15 ! m 2971 write(label,'(a)') 'a mirror plane' 2972 else if(iholohedry==4 .and. center==0)then ! primitive tetragonal 2973 2974 if(directiontype==1)then 2975 type_axis=18 ! primary n 2976 write(label,'(a)') 'a primary n plane' 2977 else if(directiontype==2)then 2978 if(sum(abs(trialt(:)-(/half,zero,zero/)))<nzero .or. sum(abs(trialt(:)-(/zero,half,zero/)))<nzero)then 2979 type_axis=16 ! secondary a or b 2980 write(label,'(a)') 'a secondary a or b plane' 2981 else if(sum(abs(trialt(:)-(/zero,zero,half/)))<nzero)then 2982 type_axis=17 ! secondary c 2983 write(label,'(a)') 'a secondary c plane' 2984 else 2985 type_axis=18 ! secondary n 2986 write(label,'(a)') 'a secondary n plane' 2987 end if ! directiontype==2 2988 else if(directiontype==3)then 2989 if( abs(trialt(3))<nzero )then 2990 type_axis=15 ! tertiary m 2991 write(label,'(a)') 'a tertiary m plane' 2992 else if( abs(trialt(3)-half)<nzero )then 2993 type_axis=19 ! tertiary c 2994 write(label,'(a)') 'a tertiary c plane' 2995 end if 2996 end if 2997 2998 else if(iholohedry==4 .and. center==-1)then ! inner tetragonal 2999 3000 if(directiontype==1 .or. directiontype==2)then 3001 if(sum(abs(trialt(:)-(/half,zero,zero/)))<nzero .or. & 3002 sum(abs(trialt(:)-(/zero,half,zero/)))<nzero .or. & 3003 sum(abs(trialt(:)-(/zero,zero,half/)))<nzero )then 3004 type_axis=16 ! a, b, or c 3005 write(label,'(a)') 'an a, b or c plane' 3006 else if(sum(abs(trialt(:)-(/half,half,zero/)))<nzero .or. & 3007 sum(abs(trialt(:)-(/zero,half,half/)))<nzero .or. & 3008 sum(abs(trialt(:)-(/half,zero,half/)))<nzero )then 3009 type_axis=15 ! n plane, equivalent to m 3010 write(label,'(a)') 'a m plane' 3011 end if ! directiontype==1 or 2 3012 else if(directiontype==3)then 3013 if( abs(trialt(3))<nzero .or. abs(trialt(3)-half)<nzero )then 3014 type_axis=15 ! tertiary c, equivalent to m 3015 write(label,'(a)') 'a tertiary m plane' 3016 else 3017 type_axis=17 ! tertiary d 3018 write(label,'(a)') 'a tertiary d plane' 3019 end if 3020 end if 3021 3022 else if(iholohedry==5)then ! hR 3023 3024 if( abs(sum(abs(trialt(:)))-one) < nzero) then 3025 type_axis=15 ! secondary m 3026 write(label,'(a)') 'a secondary m plane' 3027 else if( abs(sum(abs(trialt(:)))-half) < nzero .or. abs(sum(abs(trialt(:)))-three*half) < nzero )then 3028 type_axis=16 ! secondary c 3029 write(label,'(a)') 'a secondary c plane' 3030 end if 3031 3032 else if(iholohedry==6)then ! hP 3033 3034 if(directiontype==1)then 3035 if( abs(trialt(3)) <nzero )then 3036 type_axis=15 ! primary m 3037 write(label,'(a)') 'a primary m plane' 3038 end if 3039 else if(directiontype==2)then 3040 if( abs(trialt(3)) <nzero )then 3041 type_axis=15 ! secondary m 3042 write(label,'(a)') 'a secondary m plane' 3043 else if( abs(trialt(3)-half) < nzero ) then 3044 type_axis=16 ! secondary c 3045 write(label,'(a)') 'a secondary c plane' 3046 end if 3047 else if(directiontype==3)then 3048 if( abs(trialt(3)) <nzero )then 3049 type_axis=17 ! tertiary m 3050 write(label,'(a)') 'a tertiary m plane' 3051 else if( abs(trialt(3)-half) < nzero ) then 3052 type_axis=18 ! tertiary c 3053 write(label,'(a)') 'a tertiary c plane' 3054 end if 3055 end if ! directiontype 3056 3057 ! else if(iholohedry==7 .and. center==0)then ! cP 3058 else if(iholohedry==7)then ! cP 3059 3060 if(directiontype==1)then 3061 if((sum(abs(isymrelconv(:,:)-mirrorx(:,:)))==0 .and. & 3062 sum(abs(two*abs(trialt(:))-(/zero,half,half/)))<nzero ).or. & 3063 (sum(abs(isymrelconv(:,:)-mirrory(:,:)))==0 .and. & 3064 sum(abs(two*abs(trialt(:))-(/half,zero,half/)))<nzero ).or. & 3065 (sum(abs(isymrelconv(:,:)-mirrorz(:,:)))==0 .and. & 3066 sum(abs(two*abs(trialt(:))-(/half,half,zero/)))<nzero ) ) then 3067 type_axis=17 ! d 3068 write(label,'(a)') 'a d plane' 3069 else 3070 type_axis=18 ! primary n 3071 write(label,'(a)') 'a primary n plane' 3072 end if 3073 else if(directiontype==3)then 3074 if(sum(abs(two*abs(trialt(:))-(/half,half,half/)))<nzero )then 3075 type_axis=17 ! d 3076 write(label,'(a)') 'a d plane' 3077 else if( abs(sum(abs(trialt(:)))-half) < nzero .or. abs(sum(abs(trialt(:)))-three*half) < nzero ) then 3078 type_axis=18 ! tertiary n 3079 write(label,'(a)') 'a tertiary n plane' 3080 else if( abs(sum(abs(trialt(:)))-one) < nzero )then 3081 type_axis=15 ! tertiary m 3082 write(label,'(a)') 'a tertiary m plane' 3083 end if 3084 end if 3085 3086 ! Now, treat all other cases (including other centered Bravais lattices) 3087 else if(sum(abs(trialt(:)-(/half,zero,zero/)))<nzero .or. & 3088 sum(abs(trialt(:)-(/zero,half,zero/)))<nzero .or. & 3089 sum(abs(trialt(:)-(/zero,zero,half/)))<nzero )then 3090 type_axis=16 ! a, b or c 3091 write(label,'(a)') 'an a,b, or c plane' 3092 else if( (directiontype==1 .or. directiontype==2) .and. & 3093 (sum(abs(trialt(:)-(/half,half,zero/)))<nzero .or. & 3094 sum(abs(trialt(:)-(/zero,half,half/)))<nzero .or. & 3095 sum(abs(trialt(:)-(/half,zero,half/)))<nzero ) )then 3096 type_axis=18 ! n 3097 write(label,'(a)') 'an n plane' 3098 else if( directiontype==3 .and. sum(abs(trialt(:)-(/half,half,half/)))<nzero )then 3099 type_axis=18 ! n 3100 write(label,'(a)') 'an n plane' 3101 else if((sum(abs(isymrelconv(:,:)-mirrorx(:,:)))==0 .and. & 3102 sum(abs(two*abs(trialt(:))-(/zero,half,half/)))<nzero ).or. & 3103 (sum(abs(isymrelconv(:,:)-mirrory(:,:)))==0 .and. & 3104 sum(abs(two*abs(trialt(:))-(/half,zero,half/)))<nzero ).or. & 3105 (sum(abs(isymrelconv(:,:)-mirrorz(:,:)))==0 .and. & 3106 sum(abs(two*abs(trialt(:))-(/half,half,zero/)))<nzero ) ) then 3107 type_axis=17 ! d 3108 write(label,'(a)') 'a d plane' 3109 else if( directiontype==3 .and. sum(abs(two*abs(trialt(:))-(/half,half,half/)))<nzero)then 3110 type_axis=17 ! d 3111 write(label,'(a)') 'a d plane' 3112 else 3113 type_axis=19 ! g (all other planes with 3114 ! unconventional glide vector) 3115 write(label,'(a)') 'a g plane' 3116 end if 3117 3118 if (verbose) then 3119 write(msg,'(a,i3,a,a)')' symplanes : the symmetry operation no. ',isym,' is ', trim(label) 3120 call wrtout(std_out,msg) 3121 end if 3122 3123 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
ierr= (at input) if present, will deal with error code outside of the routine.
(at output) return 0 if no problem, 1 otherwise
symrel(3,3,nsym)=symmetry operations in real space in terms
of primitive translations rprimd at input and rprimd_new at output
SOURCE
1000 subroutine symrelrot(nsym, rprimd, rprimd_new, symrel, tolsym, ierr) 1001 1002 !Arguments ------------------------------------ 1003 !scalars 1004 integer,intent(in) :: nsym 1005 integer,intent(inout),optional :: ierr 1006 real(dp),intent(in) :: tolsym 1007 !arrays 1008 integer,intent(inout) :: symrel(3,3,nsym) 1009 real(dp),intent(in) :: rprimd(3,3),rprimd_new(3,3) 1010 1011 !Local variables------------------------------- 1012 !scalars 1013 integer :: ierr_,ii,isym,jj 1014 real(dp) :: val 1015 character(len=500) :: msg 1016 !arrays 1017 integer :: symrel_tmp(3,3,nsym) 1018 real(dp) :: coord(3,3),coordinvt(3,3),matr1(3,3),matr2(3,3),rprimd_invt(3,3) 1019 1020 !************************************************************************** 1021 1022 ierr_=0 1023 1024 !Compute the coordinates of rprimd_new in the system defined by rprimd(:,:) 1025 call matr3inv(rprimd,rprimd_invt) 1026 do ii=1,3 1027 coord(:,ii)=rprimd_new(1,ii)*rprimd_invt(1,:)+ & 1028 & rprimd_new(2,ii)*rprimd_invt(2,:)+ & 1029 & rprimd_new(3,ii)*rprimd_invt(3,:) 1030 end do 1031 1032 !Transform symmetry matrices in the system defined by rprimd_new 1033 call matr3inv(coord,coordinvt) 1034 do isym=1,nsym 1035 do ii=1,3 1036 matr1(:,ii)=symrel(:,1,isym)*coord(1,ii)+& 1037 & symrel(:,2,isym)*coord(2,ii)+& 1038 & symrel(:,3,isym)*coord(3,ii) 1039 end do 1040 do ii=1,3 1041 matr2(:,ii)=coordinvt(1,:)*matr1(1,ii)+& 1042 & coordinvt(2,:)*matr1(2,ii)+& 1043 & coordinvt(3,:)*matr1(3,ii) 1044 end do 1045 1046 !DEBUG 1047 ! write(std_out, '(a,10i4)')' symrelrot : isym, symrel=',isym,symrel(:,:,isym) 1048 ! write(std_out, '(a,9es16.6)')' transformed to ', matr2(:,:) 1049 !ENDDEBUG 1050 1051 ! Check that the new symmetry matrices are made of integers, and store them 1052 do ii=1,3 1053 do jj=1,3 1054 val=matr2(ii,jj) 1055 ! Need to allow for ten times tolsym, in case of centered Bravais lattices (but do it for all lattices ...) 1056 if(abs(val-nint(val))>ten*tolsym)then 1057 ierr_=1 1058 if(.not.(present(ierr))) then 1059 write(msg,'(2a,a,i3,a,a,3es14.6,a,a,3es14.6,a,a,3es14.6)')& 1060 'One of the components of symrel is non-integer within 10*tolsym,',ch10,& 1061 ' for isym=',isym,ch10,& 1062 ' symrel=',matr2(:,1),ch10,& 1063 ' ',matr2(:,2),ch10,& 1064 ' ',matr2(:,3) 1065 ABI_ERROR_CLASS(msg, "TolSymError") 1066 endif 1067 end if 1068 symrel_tmp(ii,jj,isym)=nint(val) 1069 end do 1070 end do 1071 end do ! isym 1072 1073 if(ierr_==0)then 1074 ! Upgrade symrel only if there is no error 1075 symrel(:,:,:)=symrel_tmp(:,:,:) 1076 endif 1077 1078 if(present(ierr))then 1079 ierr=ierr_ 1080 endif 1081 1082 end subroutine symrelrot