TABLE OF CONTENTS
- ABINIT/irreducible_set_pert
- ABINIT/m_geometry
- m_geometry/acrossb
- m_geometry/bonds_lgth_angles
- m_geometry/chkdilatmx
- m_geometry/chkrprimd
- m_geometry/det3r
- m_geometry/dist2
- m_geometry/fcart2gred
- m_geometry/fixsym
- m_geometry/getspinrot
- m_geometry/gred2fcart
- m_geometry/ioniondist
- m_geometry/littlegroup_pert
- m_geometry/metric
- m_geometry/mkradim
- m_geometry/mkrdim
- m_geometry/normv_int_vector
- m_geometry/normv_int_vector_array
- m_geometry/normv_rdp_vector
- m_geometry/normv_rdp_vector_array
- m_geometry/phdispl_cart2red
- m_geometry/randomcellpos
- m_geometry/remove_inversion
- m_geometry/rotmat
- m_geometry/shellstruct
- m_geometry/spinrot_cmat
- m_geometry/strainsym
- m_geometry/strconv
- m_geometry/stress_voigt_to_mat
- m_geometry/stresssym
- m_geometry/symredcart
- m_geometry/vdotw_rc_vector
- m_geometry/vdotw_rr_vector
- m_geometry/wedge_basis
- m_geometry/wedge_product
- m_geometry/wigner_seitz
- m_geometry/xcart2xred
- m_geometry/xred2xcart
- m_symtk/reduce2primitive
ABINIT/irreducible_set_pert [ Functions ]
NAME
irreducible_set_pert
FUNCTION
Determines a set of perturbations that form a basis in that, using symmetry, they can be used to generate all other perturbations that are asked to be calculated (target).
INPUTS
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. mpert =maximum number of iper natom= number of atoms nsym=number of space group symmetries rfdir(3)=direction for the perturbations rfpert(mpert)=information on the perturbations symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space) symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space) symq(4,2,nsym)= (integer) three first numbers define the G vector; fourth number is 0 if the q-vector is not preserved, is 1 otherwise second index is one without time-reversal symmetry, two with time-reversal symmetry
OUTPUT
pertsy(3,mpert)= the target perturbation is described by the two last indices (idir, and ipert),
the value is 0, 1 or -1, see notes.
NOTES
Output will be in the pertsy array, 0 for non-target perturbations 1 for basis perturbations -1 for perturbations that can be found from basis perturbations
SOURCE
3475 subroutine irreducible_set_pert(indsym,mpert,natom,nsym,pertsy,rfdir,rfpert,symq,symrec,symrel) 3476 3477 !Arguments ------------------------------- 3478 !scalars 3479 integer,intent(in) :: mpert,natom,nsym 3480 !arrays 3481 integer,intent(in) :: indsym(4,nsym,natom),rfdir(3),rfpert(mpert) 3482 integer,intent(in) :: symq(4,2,nsym),symrec(3,3,nsym),symrel(3,3,nsym) 3483 integer,intent(out) :: pertsy(3,mpert) 3484 3485 !Local variables ------------------------- 3486 !scalars 3487 integer :: found,idir1,idisy1,ii,ipert1,ipesy1,isign,isym,itirev,jj 3488 !arrays 3489 integer :: sym1(3,3) 3490 3491 ! ********************************************************************* 3492 3493 !Zero pertsy 3494 pertsy(:,:)=0 3495 3496 do ipert1=1,mpert 3497 do idir1=1,3 3498 if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then 3499 ! write(std_out,*)' for candidate idir =',idir1,' ipert = ',ipert1 3500 3501 ! Loop on all symmetries, including time-reversal 3502 do isym=1,nsym 3503 do itirev=1,2 3504 isign=3-2*itirev 3505 3506 if(symq(4,itirev,isym)/=0)then 3507 3508 found=1 3509 3510 ! Here select the symmetric of ipert1 3511 if(ipert1<=natom)then 3512 ipesy1=indsym(4,isym,ipert1) 3513 do ii=1,3 3514 do jj=1,3 3515 sym1(ii,jj)=symrec(ii,jj,isym) 3516 end do 3517 end do 3518 else if(ipert1==(natom+2) .or. ipert1==(natom+6))then 3519 ipesy1=ipert1 3520 do ii=1,3 3521 do jj=1,3 3522 sym1(ii,jj)=symrel(ii,jj,isym) 3523 end do 3524 end do 3525 else 3526 found=0 3527 end if 3528 3529 ! Now that a symmetric perturbation has been obtained, 3530 ! including the expression of the symmetry matrix, see 3531 ! if the symmetric perturbations are available 3532 if( found==1 ) then 3533 3534 do idisy1=1,3 3535 if(sym1(idir1,idisy1)/=0)then 3536 if(pertsy(idisy1,ipesy1)==0)then 3537 found=0 3538 exit 3539 end if 3540 end if 3541 end do 3542 end if 3543 3544 ! Now, if still found, then it is a symmetric 3545 ! of some linear combination of existing perturbations 3546 if(found==1)then 3547 3548 ! DEBUG 3549 ! write(std_out,*)' all found ! isym, isign= ',isym,isign 3550 ! write(std_out,1010)((sym1(ii,jj),ii=1,3),jj=1,3) 3551 ! write(std_out,1010)((sym2(ii,jj),ii=1,3),jj=1,3) 3552 ! write(std_out,*)sumr,sumi 3553 ! 1010 format(9i4) 3554 ! ENDDEBUG 3555 3556 pertsy(idir1,ipert1)=-1 3557 exit ! Exit loop on symmetry operations 3558 3559 end if 3560 3561 end if ! End loop on all symmetries + time-reversal 3562 end do 3563 end do 3564 3565 ! Now that all symmetries have been examined, 3566 ! if still not symmetric of a linear combination 3567 ! of basis perturbations, then it is a basis perturbation 3568 if(pertsy(idir1,ipert1)/=-1) pertsy(idir1,ipert1)=1 3569 ! write(std_out,'(a,3i5)' ) ' irreducible_set_pert :',idir1,ipert1,pertsy(idir1,ipert1) 3570 3571 end if ! End big loop on all elements 3572 end do 3573 end do 3574 3575 end subroutine irreducible_set_pert
ABINIT/m_geometry [ Modules ]
NAME
m_geometry
FUNCTION
This module contains basic tools to operate on vectors expressed in reduced coordinates.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG, MT, FJ, TRangel, DCA, XG, AHR, DJA, DRH) 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_geometry 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_atomdata 28 use m_sort 29 30 use m_io_tools, only : open_file 31 use m_numeric_tools, only : uniformrandom, isinteger, set2unit 32 use m_symtk, only : mati3inv, mati3det, matr3inv, symdet 33 use m_hide_lapack, only : matr3eigval 34 use m_pptools, only : prmat 35 use m_numeric_tools, only : wrap2_pmhalf 36 use m_hide_lapack, only : matrginv 37 38 implicit none 39 40 private 41 42 public :: normv ! Norm of vector(s) in reduced coordinates either in real or reciprocal space. 43 public :: vdotw ! Scalar product between two reduced vectors either in real or reciprocal space. 44 public :: acrossb ! Cross product of two 3-vectors. 45 public :: wigner_seitz ! Find the grid of points falling inside the Wigner-Seitz cell. 46 public :: phdispl_cart2red ! Calculate the displacement vectors for all branches in reduced coordinates. 47 public :: getspinrot ! Compute the components of the spinor rotation matrix 48 public :: spinrot_cmat ! Construct 2x2 complex matrix representing rotation operator in spin-space. 49 public :: rotmat ! Finds the rotation matrix. 50 public :: fixsym ! Check that iatfix does not break symmetry. 51 public :: det3r ! Compute determinant of a 3x3 real matrix 52 public :: metric ! Compute metric matrices. 53 public :: mkradim ! Make rprim and acell from rprimd 54 public :: mkrdim ! Make rprimd from acell from rprim 55 public :: chkrprimd ! Test if {rprim,acell,rprimd} are consistent 56 public :: chkdilatmx ! check if dilatation of unit cell is consistent with initial G-sphere 57 public :: xcart2xred ! From cart coords to reduced 58 public :: xred2xcart ! From reduced coords to cart. 59 public :: gred2fcart ! Convert reduced gradients into cartesian forces 60 public :: fcart2gred ! Convert cartesian forces into reduced gradients 61 public :: bonds_lgth_angles ! Write GEO file 62 public :: randomcellpos ! Creates unit cell with random atomic positions. 63 public :: ioniondist ! Compute ion-ion distances 64 public :: dist2 ! Calculates the distance of v1 and v2 in a crystal by repeating the unit cell 65 public :: shellstruct ! Calculates shell structure (multiplicities, radii) 66 public :: remove_inversion ! Remove the inversion symmetry and improper rotations 67 public :: reduce2primitive ! Find real space primitive vectors from non-primitive ones and set of translations 68 public :: symredcart ! Convert a symmetry operation from reduced coordinates (integers) to cart coords (reals) 69 public :: strainsym ! Symmetrize the strain tensor. 70 public :: stresssym ! Symmetrize the stress tensor. 71 public :: stress_voigt_to_mat! Build 3x3 symmetric stress tensor from stress vector in Voigt notation. 72 public :: strconv ! Convert from symmetric storage mode in reduced coords to cart coords. 73 public :: littlegroup_pert ! Determines the set of symmetries that leaves a perturbation invariant. 74 public :: irreducible_set_pert ! Determines a set of perturbations that form a basis 75 public :: wedge_basis ! compute rprimd x gprimd vectors needed for generalized cross product 76 public :: wedge_product ! compute wedge product given wedge basis 77 78 interface normv 79 module procedure normv_rdp_vector 80 module procedure normv_int_vector 81 !module procedure normv_int_vector_array ! WARNING for the time being, do not use these 2 procedures, 82 !module procedure normv_rdp_vector_array ! sunstudio12 is not able to resolve which sub should be called. 83 end interface normv 84 85 interface vdotw 86 module procedure vdotw_rr_vector 87 module procedure vdotw_rc_vector 88 end interface vdotw 89 90 CONTAINS !===========================================================
m_geometry/acrossb [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
acrossb
FUNCTION
Calculates the cross product of two 3-vectors
INPUTS
a(3): real(dp) vector b(3): real(dp) vector
OUTPUT
c(3): real(dp) vector = a X b
SOURCE
396 subroutine acrossb(a,b,c) 397 398 !Arguments --------------------------------------------- 399 !arrays 400 real(dp),intent(in) :: a(3),b(3) 401 real(dp),intent(out) :: c(3) 402 403 ! ********************************************************************* 404 405 c(1) = a(2)*b(3) - a(3)*b(2) 406 c(2) = -a(1)*b(3) + a(3)*b(1) 407 c(3) = a(1)*b(2) - b(1)*a(2) 408 409 end subroutine acrossb
m_geometry/bonds_lgth_angles [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
bonds_lgth_angles
FUNCTION
From list of coordinates and primitive translations, output a list of bonds lengths and bond angles.
INPUTS
coordn = maximum coordination number to be taken into account fnameabo_app_geo=name of file for _GEO data natom = number of atoms in unit cell ntypat = number of types of atoms in unit cell. rprimd(3,3) = real space dimensional primitive translations (bohr) typat(natom) = type integer for each atom in cell znucl(ntypat)= real(dp), atomic number of atom type xred(3,natom)= reduced coordinates of atoms
OUTPUT
data written in file fnameabo_app_geo
NOTES
The tolerance tol8 aims at giving a machine-independent ordering. (this trick is used in bonds.f, listkk.f, prtrhomxmn.f and rsiaf9.f)
SOURCE
1798 subroutine bonds_lgth_angles(coordn,fnameabo_app_geo,natom,ntypat,rprimd,typat,xred,znucl) 1799 1800 !Arguments ------------------------------------ 1801 !scalars 1802 integer,intent(in) :: coordn,natom,ntypat 1803 character(len=*),intent(in) :: fnameabo_app_geo 1804 !arrays 1805 integer,intent(in) :: typat(natom) 1806 real(dp),intent(in) :: rprimd(3,3),znucl(ntypat) 1807 real(dp),intent(inout) :: xred(3,natom) 1808 1809 !Local variables------------------------------- 1810 !scalars 1811 integer :: done,ia,ib,ic,ii,ineighb,jneighb,mneighb,mu,ndig,nu,t1,t2,t3,tmax,temp_unit 1812 real(dp) :: adotb,asq,bsq,co,length,sq,thdeg 1813 !real(dp)u1,u2,u3,v1,v2,v3 1814 character(len=500) :: msg 1815 type(atomdata_t) :: atom 1816 !arrays 1817 integer,allocatable :: list_neighb(:,:,:) 1818 real(dp) :: bab(3),bac(3),dif(3),rmet(3,3) 1819 real(dp),allocatable :: sqrlength(:),xcart(:,:) 1820 character(len=8),allocatable :: iden(:) 1821 1822 ! ************************************************************************* 1823 1824 !Initialize the file 1825 write(msg, '(a,a,a)' )' bonds_lgth_angles : about to open file ',trim(fnameabo_app_geo),ch10 1826 call wrtout(std_out,msg,'COLL'); call wrtout(ab_out,msg,'COLL') 1827 1828 if (open_file(fnameabo_app_geo,msg,newunit=temp_unit,status='unknown',form='formatted') /= 0) then 1829 ABI_ERROR(msg) 1830 end if 1831 rewind(temp_unit) 1832 1833 write(msg, '(a,a)' ) ch10,' ABINIT package : GEO file ' 1834 call wrtout(temp_unit,msg,'COLL') 1835 1836 !Compute maximum number of neighbors is the neighbor list, 1837 !from the indicative coordination number 1838 !Note : the following formula includes next nearest neighbors, but not others 1839 mneighb=1+coordn+coordn*(coordn-1) 1840 1841 write(msg, '(a,a,i2,a,a,i4,a,a,a,i4,a)' ) ch10,& 1842 & ' Maximal coordination number, as estimated by the user : ',coordn,ch10,& 1843 & ' giving a maximum of ',coordn*coordn,& 1844 & ' nearest neighbors and next nearest neighbors, ',ch10,& 1845 & ' and ',(coordn*(coordn-1))/2,& 1846 & ' distinct angles between nearest neighbors' 1847 call wrtout(temp_unit,msg,'COLL') 1848 1849 !Compute metric tensor in real space rmet 1850 do nu=1,3 1851 do mu=1,3 1852 rmet(mu,nu)=rprimd(1,mu)*rprimd(1,nu)+& 1853 & rprimd(2,mu)*rprimd(2,nu)+& 1854 & rprimd(3,mu)*rprimd(3,nu) 1855 end do 1856 end do 1857 1858 write(msg, '(a,a)' )ch10,' Primitive vectors of the periodic cell (bohr)' 1859 call wrtout(temp_unit,msg,'COLL') 1860 do nu=1,3 1861 write(msg, '(1x,a,i1,a,3f10.5)' ) ' R(',nu,')=',rprimd(:,nu) 1862 call wrtout(temp_unit,msg,'COLL') 1863 end do 1864 1865 write(msg, '(a,a)' ) ch10,& 1866 & ' Atom list Reduced coordinates Cartesian coordinates (bohr)' 1867 call wrtout(temp_unit,msg,'COLL') 1868 1869 !Set up a list of character identifiers for all atoms : iden(ia) 1870 ABI_MALLOC(iden,(natom)) 1871 iden(:)=' ' 1872 do ia=1,natom 1873 ndig=int(log10(dble(ia)+0.5d0))+1 1874 call atomdata_from_znucl(atom,znucl(typat(ia))) 1875 if(ndig==1) write(iden(ia), '(a,a,i1,a)' ) atom%symbol,'(',ia,') ' 1876 if(ndig==2) write(iden(ia), '(a,a,i2,a)' ) atom%symbol,'(',ia,') ' 1877 if(ndig==3) write(iden(ia), '(a,a,i3,a)' ) atom%symbol,'(',ia,') ' 1878 if(ndig==4) write(iden(ia), '(a,a,i4,a)' ) atom%symbol,'(',ia,')' 1879 if(ndig>4)then 1880 close(temp_unit) 1881 write(msg, '(a,i8,a,a)' )& 1882 & 'bonds_lgth_angles cannot handle more than 9999 atoms, while natom=',natom,ch10,& 1883 & 'Action: decrease natom, or contact ABINIT group.' 1884 ABI_BUG(msg) 1885 end if 1886 end do 1887 1888 !Compute cartesian coordinates, and print reduced and cartesian coordinates 1889 !then print coordinates in angstrom, with the format neede for xmol 1890 ABI_MALLOC(xcart,(3,natom)) 1891 call xred2xcart(natom,rprimd,xcart,xred) 1892 1893 do ia=1,natom 1894 write(msg, '(a,a,3f10.5,a,3f10.5)' ) & 1895 & ' ',iden(ia),(xred(ii,ia)+tol10,ii=1,3),& 1896 & ' ',(xcart(ii,ia)+tol10,ii=1,3) 1897 call wrtout(temp_unit,msg,'COLL') 1898 end do 1899 1900 write(msg, '(a,a,a,a,i4,a)' )ch10,& 1901 & ' XMOL data : natom, followed by cartesian coordinates in Angstrom',& 1902 & ch10,ch10,natom,ch10 1903 call wrtout(temp_unit,msg,'COLL') 1904 1905 do ia=1,natom 1906 call atomdata_from_znucl(atom,znucl(typat(ia))) 1907 write(msg, '(a,a,3f10.5)' )' ',atom%symbol,xcart(1:3,ia)*Bohr_Ang 1908 call wrtout(temp_unit,msg,'COLL') 1909 end do 1910 1911 ABI_FREE(xcart) 1912 1913 ABI_MALLOC(list_neighb,(0:mneighb+1,4,2)) 1914 ABI_MALLOC(sqrlength,(0:mneighb+1)) 1915 1916 !Compute list of neighbors 1917 do ia=1,natom 1918 1919 write(msg, '(a,a,a,a,a,a,a,a,a)' ) ch10,'===========',& 1920 & '=====================================================================',& 1921 & ch10,' ',iden(ia),ch10,ch10,' Bond lengths ' 1922 call wrtout(temp_unit,msg,'COLL') 1923 1924 ! Search other atoms for bonds, but must proceed 1925 ! in such a way to consider a search box sufficiently large, 1926 ! so increase the size of the search box until the 1927 ! final bond length list do not change 1928 do tmax=0,5 1929 1930 ! Set initial list of neighbors to zero, 1931 ! and initial square of bond lengths to a very large number. 1932 ! Note that the dimension is larger than neighb to ease 1933 ! the later sorting : neighbors 0 and neighb+1 are non-existent, while 1934 ! neighbor 1 will be the atom itself ... 1935 list_neighb(0:mneighb+1,1:4,1)=0 1936 sqrlength(1:mneighb+1)=huge(zero) 1937 sqrlength(0)=-1.0d0 1938 1939 ! Here search on all atoms inside the box defined by tmax 1940 do ib=1,natom 1941 do t3=-tmax,tmax 1942 do t2=-tmax,tmax 1943 do t1=-tmax,tmax 1944 dif(1)=xred(1,ia)-(xred(1,ib)+dble(t1)) 1945 dif(2)=xred(2,ia)-(xred(2,ib)+dble(t2)) 1946 dif(3)=xred(3,ia)-(xred(3,ib)+dble(t3)) 1947 sq=rsdot(dif(1),dif(2),dif(3),dif(1),dif(2),dif(3),rmet) 1948 1949 ! Insert the atom at the proper place in the neighbor list. 1950 do ineighb=mneighb,0,-1 1951 ! Note the tolerance 1952 if(sq+tol8>sqrlength(ineighb))then 1953 sqrlength(ineighb+1)=sq 1954 list_neighb(ineighb+1,1,1)=ib 1955 list_neighb(ineighb+1,2,1)=t1 1956 list_neighb(ineighb+1,3,1)=t2 1957 list_neighb(ineighb+1,4,1)=t3 1958 ! DEBUG 1959 ! if(ineighb/=mneighb)then 1960 ! write(std_out,*)' ' 1961 ! do ii=1,mneighb 1962 ! write(std_out,*)ii,sqrlength(ii) 1963 ! end do 1964 ! end if 1965 ! ENDDEBUG 1966 exit 1967 else 1968 sqrlength(ineighb+1)=sqrlength(ineighb) 1969 list_neighb(ineighb+1,1:4,1)=list_neighb(ineighb,1:4,1) 1970 end if 1971 end do 1972 1973 end do 1974 end do 1975 end do 1976 ! end ib loop: 1977 end do 1978 1979 ! Now, check that the box defined by tmax was large enough : 1980 ! require the present and old lists to be the same 1981 done=0 1982 1983 if(tmax>0)then 1984 done=1 1985 do ineighb=1,mneighb 1986 ! DEBUG 1987 ! write(std_out,'(5i5,f12.5)' )ineighb,list_neighb(ineighb,1:4,1),& 1988 ! & sqrlength(ineighb) 1989 ! write(std_out,'(5i5)' )ineighb,list_neighb(ineighb,1:4,2) 1990 ! ENDDEBUG 1991 if( list_neighb(ineighb,1,1)/=list_neighb(ineighb,1,2) .or. & 1992 & list_neighb(ineighb,2,1)/=list_neighb(ineighb,2,2) .or. & 1993 & list_neighb(ineighb,3,1)/=list_neighb(ineighb,3,2) .or. & 1994 & list_neighb(ineighb,4,1)/=list_neighb(ineighb,4,2) )then 1995 done=0 1996 end if 1997 end do 1998 end if 1999 2000 ! If done==1, then one can exit the loop : the correct list of 2001 ! neighbors is contained in list_neighb(1:neighb,1:4,1), 2002 ! with the first neighbor being the atom itself 2003 if(done==1)exit 2004 2005 ! If the work is not done, while tmax==5, then there is a problem . 2006 if(tmax==5)then 2007 close(temp_unit) 2008 write(msg, '(2a)' )& 2009 & 'Did not succeed to generate a reliable list of bonds ',& 2010 & 'since tmax is exceeded.' 2011 ABI_BUG(msg) 2012 end if 2013 2014 ! Copy the new list into the old list. 2015 list_neighb(1:mneighb,1:4,2)=list_neighb(1:mneighb,1:4,1) 2016 2017 ! Loop on tmax (note that there are exit instruction inside the loop) 2018 end do 2019 2020 2021 2022 ! Output the bond list 2023 do ineighb=2,mneighb 2024 ib=list_neighb(ineighb,1,1) 2025 length=sqrt(sqrlength(ineighb)) 2026 write(msg, '(a,a,a,a,3i2,t27,a,f10.5,a,f9.5,a)' )& 2027 & ' ',trim(iden(ia)),' - ',trim(iden(ib)),& 2028 & list_neighb(ineighb,2:4,1),'bond length is ',& 2029 & length,' bohr ( or ',Bohr_Ang*length,' Angst.)' 2030 call wrtout(temp_unit,msg,'COLL') 2031 end do 2032 2033 ! Output the angle list 2034 if(coordn>1)then 2035 2036 write(msg, '(a,a)' ) ch10,' Bond angles ' 2037 call wrtout(temp_unit,msg,'COLL') 2038 2039 do ineighb=2,coordn 2040 do jneighb=ineighb+1,coordn+1 2041 2042 ib=list_neighb(ineighb,1,1) 2043 ic=list_neighb(jneighb,1,1) 2044 do mu=1,3 2045 bab(mu)=xred(mu,ib)+dble(list_neighb(ineighb,1+mu,1))-xred(mu,ia) 2046 bac(mu)=xred(mu,ic)+dble(list_neighb(jneighb,1+mu,1))-xred(mu,ia) 2047 end do 2048 asq=rsdot(bab(1),bab(2),bab(3),bab(1),bab(2),bab(3),rmet) 2049 bsq=rsdot(bac(1),bac(2),bac(3),bac(1),bac(2),bac(3),rmet) 2050 adotb=rsdot(bab(1),bab(2),bab(3),bac(1),bac(2),bac(3),rmet) 2051 co=adotb/sqrt(asq*bsq) 2052 if( abs(co)-1.0d0 >= 0.0d0 )then 2053 if( abs(co)-1.0d0 <= 1.0d-12 )then 2054 ! Allows for a small numerical inaccuracy 2055 thdeg=0.0d0 2056 if(co < 0.0d0) thdeg=180.0d0 2057 else 2058 ABI_BUG('the evaluation of the angle is wrong.') 2059 end if 2060 else 2061 thdeg=acos(co)*180.d0*piinv 2062 end if 2063 2064 write(msg, '(a,a,3i2,a,a,a,a,3i2,t44,a,f13.5,a)' )& 2065 & ' ',trim(iden(ib)),list_neighb(ineighb,2:4,1),' - ',& 2066 & trim(iden(ia)),' - ',trim(iden(ic)),& 2067 & list_neighb(jneighb,2:4,1),'bond angle is ',thdeg,' degrees ' 2068 call wrtout(temp_unit,msg,'COLL') 2069 end do 2070 end do 2071 2072 end if 2073 end do ! End big ia loop: 2074 2075 ABI_FREE(iden) 2076 ABI_FREE(list_neighb) 2077 ABI_FREE(sqrlength) 2078 2079 close(temp_unit) 2080 2081 contains 2082 2083 function rsdot(u1,u2,u3,v1,v2,v3,rmet) 2084 2085 real(dp) :: rsdot 2086 real(dp),intent(in) :: u1,u2,u3,v1,v2,v3 2087 real(dp),intent(in) :: rmet(3,3) 2088 rsdot=rmet(1,1)*u1*v1+rmet(2,1)*u2*v1+& 2089 & rmet(3,1)*u3*v1+rmet(1,2)*u1*v2+rmet(2,2)*u2*v2+& 2090 & rmet(3,2)*u3*v2+rmet(1,3)*u1*v3+rmet(2,3)*u2*v3+rmet(3,3)*u3*v3 2091 end function rsdot 2092 2093 end subroutine bonds_lgth_angles
m_geometry/chkdilatmx [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
chkdilatmx
FUNCTION
Check whether the new rprimd does not give a too large number of plane waves, compared to the one booked for rprimd, taking into account the maximal dilatation dilatmx. Actually check whether the new Fermi sphere is inside the old one, dilated.
INPUTS
chkdilatmx_ = if 1, will prevent to have any vector outside the Fermi sphere, possibly by rescaling (three times at most), and then stopping the execution if 0, simply send a warning, but continues execution dilatmx = maximal dilatation factor (usually the input variable) rprimd = new primitive vectors rprimd_orig = original primitive vectors (usually the input variable)
OUTPUT
dilatmx_errmsg=Emptry string if calculation can continue.
If the calculation cannot continue, dilatmx_errmsg will contain
the message that should be reported in the output file.
Client code should handle a possible problem with the following test:
if (LEN_TRIM(dilatmx_errmsg) then
dump dilatmx_errmsg to the main output file.
handle_error
end if
SOURCE
1430 subroutine chkdilatmx(chkdilatmx_,dilatmx,rprimd,rprimd_orig,dilatmx_errmsg) 1431 1432 !Arguments ------------------------------------ 1433 !scalars 1434 integer,intent(in) :: chkdilatmx_ 1435 real(dp),intent(in) :: dilatmx 1436 character(len=500),intent(out) :: dilatmx_errmsg 1437 !arrays 1438 real(dp),intent(inout) :: rprimd(3,3) 1439 real(dp),intent(in) :: rprimd_orig(3,3) 1440 1441 !Local variables------------------------------- 1442 !scalars 1443 integer :: ii,jj,mu 1444 real(dp) :: alpha,dilatmx_new 1445 !arrays 1446 real(dp) :: eigval(3),gprimd_orig(3,3),met(3,3),old_to_new(3,3) 1447 character(len=500) :: msg 1448 1449 ! ************************************************************************* 1450 1451 !Generates gprimd 1452 call matr3inv(rprimd_orig,gprimd_orig) 1453 1454 !Find the matrix that transform an original xcart to xred, then to the new xcart 1455 do mu=1,3 1456 old_to_new(mu,:)=rprimd(mu,1)*gprimd_orig(:,1)+& 1457 & rprimd(mu,2)*gprimd_orig(:,2)+& 1458 & rprimd(mu,3)*gprimd_orig(:,3) 1459 end do 1460 1461 !The largest increase in length will be obtained thanks 1462 !to the diagonalization of the corresponding metric matrix : 1463 !it is the square root of its largest eigenvalue. 1464 do ii=1,3 1465 do jj=1,3 1466 met(ii,jj)=old_to_new(1,ii)*old_to_new(1,jj)+& 1467 & old_to_new(2,ii)*old_to_new(2,jj)+& 1468 & old_to_new(3,ii)*old_to_new(3,jj) 1469 end do 1470 end do 1471 1472 call matr3eigval(eigval,met) 1473 1474 dilatmx_new=sqrt(maxval(eigval(:))) 1475 1476 dilatmx_errmsg = "" 1477 if(dilatmx_new>dilatmx+tol6)then 1478 1479 ! MJV 2014 07 22: correct rprim to maximum jump allowed by dilatmx 1480 ! XG 20171011 : eigenvalues of "old_to_old" tensor are of course the unity ! 1481 1482 if(chkdilatmx_/=0)then 1483 alpha = (dilatmx - one) / (dilatmx_new - one) 1484 ! for safety, only 90 percent of max jump 1485 alpha = 0.9_dp * alpha 1486 1487 rprimd = alpha * rprimd + (one - alpha) * rprimd_orig 1488 1489 write(dilatmx_errmsg,'(3a,es16.6,4a,es16.6,2a,es16.6,a)')& 1490 'The new primitive vectors rprimd (an evolving quantity)',ch10,& 1491 'are too large with respect to the old rprimd and the accompanying dilatmx: ',dilatmx,ch10,& 1492 'This large change of unit cell parameters is not allowed by the present value of dilatmx.',ch10,& 1493 'An adequate value would have been dilatmx_new= ',dilatmx_new,ch10,& 1494 'Calculation continues with limited jump, by rescaling the projected move by the factor: ',alpha,'.' 1495 else 1496 write(msg, '(3a,es16.6,2a,es16.6,2a)' )& 1497 'The new primitive vectors rprimd (an evolving quantity)',ch10,& 1498 'are too large, given the initial rprimd and the accompanying dilatmx: ',dilatmx,ch10,& 1499 'An adequate value would have been dilatmx_new= ',dilatmx_new,ch10,& 1500 'As chkdilatmx=0, assume experienced user. Execution will continue.' 1501 ABI_WARNING(msg) 1502 end if 1503 1504 end if 1505 1506 end subroutine chkdilatmx
m_geometry/chkrprimd [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
chkrprimd
FUNCTION
Test if {rprim,acell,rprimd} are consistent
It means that rprimd can be reconstructed from the rprim and acell
Output a message if is not the case
INPUTS
OUTPUT
(only writing)
SOURCE
1349 subroutine chkrprimd(acell,rprim,rprimd,iout) 1350 1351 !Arguments ------------------------------------ 1352 !scalars 1353 integer,intent(in) :: iout 1354 !arrays 1355 real(dp),intent(in) :: rprim(3,3) 1356 real(dp),intent(in) :: rprimd(3,3) 1357 real(dp),intent(in) :: acell(3) 1358 1359 !Local variables------------------------------- 1360 !scalars 1361 integer :: ii,jj 1362 !arrays 1363 real(dp) :: rprimd_test(3,3) 1364 logical :: equal 1365 1366 ! *********************************************************** 1367 1368 !########################################################### 1369 !### 1. Compute rprimd from rprim and acell 1370 do ii=1,3 1371 do jj=1,3 1372 rprimd_test(ii,jj)=rprim(ii,jj)*acell(jj) 1373 end do 1374 end do 1375 1376 1377 !########################################################### 1378 !### 2. Compare rprimd and rprimd_test 1379 1380 equal=.TRUE. 1381 do ii=1,3 1382 do jj=1,3 1383 if (abs(rprimd_test(ii,jj)-rprimd(ii,jj))>1.E-12) then 1384 equal=.FALSE. 1385 end if 1386 end do 1387 end do 1388 1389 if (equal)then 1390 write(iout,*) 'chkrprimd: rprimd is consistent' 1391 else 1392 write(iout,*) 'chkrprimd: rprimd is NOT consistent ERROR' 1393 end if 1394 1395 end subroutine chkrprimd
m_geometry/det3r [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
det3r
FUNCTION
Compute determinant of a 3x3 real matrix
SOURCE
1156 pure real(dp) function det3r(rprimd) 1157 1158 !Arguments ------------------------------------ 1159 real(dp),intent(in) :: rprimd(3,3) 1160 1161 ! ************************************************************************* 1162 1163 ! Compute unit cell volume 1164 det3r = rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+& 1165 rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+& 1166 rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3)) 1167 1168 end function det3r
m_geometry/dist2 [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
dist2
FUNCTION
Calculates the distance of v1 and v2 in a crystal by repeating the unit cell
INPUTS
v1,v2
rprimd: dimensions of the unit cell. if not given 1,0,0/0,1,0/0,0,1 is assumed
option: 0 v1, v2 given in cartesian coordinates (default)
1 v1,v2 given in reduced coordinates
-1 v1 and v2 are supposed equal, and the routine returns the length of the smallest Bravais lattice vector
OUTPUT
dist2
SOURCE
2597 function dist2(v1,v2,rprimd,option) 2598 2599 !Arguments ------------------------------------ 2600 !scalars 2601 integer,intent(in),optional :: option 2602 real(dp) :: dist2 2603 !arrays 2604 real(dp),intent(in),optional :: rprimd(3,3) 2605 real(dp),intent(in) :: v1(3),v2(3) 2606 2607 !Local variables------------------------------- 2608 !scalars 2609 integer :: i1,i2,i3,opt,s1,s2,s3 2610 real(dp):: min2,norm2,ucvol 2611 !arrays 2612 integer :: limits(3) 2613 real(dp) :: corner(3),dred(3),dtot(3),dv(3),dwrap(3),sh(3) 2614 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 2615 real(dp) :: vprimd(3,3) 2616 2617 ! ************************************************************************* 2618 2619 if (.not.PRESENT(rprimd)) then 2620 vprimd=reshape((/1,0,0, 0,1,0, 0,0,1/),(/3,3/)) 2621 else 2622 vprimd=rprimd 2623 end if 2624 2625 call metric(gmet,gprimd,-1,rmet,vprimd,ucvol) 2626 2627 dv(:)=v2(:)-v1(:) 2628 2629 !If in cartesian coordinates, need to be transformed to reduced coordinates. 2630 opt=0 2631 if(present(option))then 2632 opt=option 2633 end if 2634 if(opt==0)then 2635 dred(:)=gprimd(1,:)*dv(1)+gprimd(2,:)*dv(2)+gprimd(3,:)*dv(3) 2636 else if(opt==1)then 2637 dred(:)=dv(:) 2638 else if(opt==-1)then 2639 dred(:)=zero 2640 end if 2641 2642 !Wrap in the ]-1/2,1/2] interval 2643 call wrap2_pmhalf(dred(1),dwrap(1),sh(1)) 2644 call wrap2_pmhalf(dred(2),dwrap(2),sh(2)) 2645 call wrap2_pmhalf(dred(3),dwrap(3),sh(3)) 2646 2647 !Compute the limits of the parallelipipedic box that contains the Wigner-Seitz cell 2648 !The reduced coordinates of the corners of the Wigner-Seitz cell are computed (multiplied by two) 2649 !Then, the maximal values of these reduced coordinates are stored. 2650 limits(:)=0 2651 do s1=-1,1,2 2652 do s2=-1,1,2 2653 do s3=-1,1,2 2654 corner(:)=gmet(:,1)*s1*rmet(1,1)+gmet(:,2)*s2*rmet(2,2)+gmet(:,3)*s3*rmet(3,3) 2655 limits(1)=max(limits(1),ceiling(abs(corner(1))+tol14)) 2656 limits(2)=max(limits(2),ceiling(abs(corner(2))+tol14)) 2657 limits(3)=max(limits(3),ceiling(abs(corner(3))+tol14)) 2658 end do 2659 end do 2660 end do 2661 2662 !Use all relevant primitive real space lattice vectors to find the minimal difference vector 2663 min2=huge(zero) 2664 do i1=-limits(1),limits(1) 2665 dtot(1)=dwrap(1)+i1 2666 do i2=-limits(2),limits(2) 2667 dtot(2)=dwrap(2)+i2 2668 do i3=-limits(3),limits(3) 2669 if(opt/=-1.or.i1/=0.or.i2/=0.or.i3/=0)then 2670 dtot(3)=dwrap(3)+i3 2671 norm2=dtot(1)*rmet(1,1)*dtot(1)+dtot(2)*rmet(2,2)*dtot(2)+dtot(3)*rmet(3,3)*dtot(3)+& 2672 & 2*(dtot(1)*rmet(1,2)*dtot(2)+dtot(2)*rmet(2,3)*dtot(3)+dtot(3)*rmet(3,1)*dtot(1)) 2673 min2=min(norm2,min2) 2674 endif 2675 end do 2676 end do 2677 end do 2678 dist2=sqrt(min2) 2679 2680 end function dist2
m_geometry/fcart2gred [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
fcart2gred
FUNCTION
Convert cartesian forces into reduced forces
INPUTS
fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr) natom=Number of atoms in the unitary cell rprimd(3,3)=dimensional primitive
OUTPUT
gred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
NOTES
Unlike gred, fcart has been corrected by enforcing the translational symmetry, namely that the sum of force on all atoms is zero.
SOURCE
1734 subroutine fcart2gred(fcart,gred,rprimd,natom) 1735 1736 !Arguments ------------------------------------ 1737 !scalars 1738 integer,intent(in) :: natom 1739 !arrays 1740 real(dp),intent(in) :: fcart(3,natom) 1741 real(dp),intent(out) :: gred(3,natom) 1742 real(dp),intent(in) :: rprimd(3,3) 1743 1744 !Local variables------------------------------- 1745 !scalars 1746 integer :: iatom,mu 1747 1748 ! ************************************************************************* 1749 1750 !MT, april 2012: the coding was not consistent with gred2fcart 1751 do iatom=1,natom 1752 do mu=1,3 1753 gred(mu,iatom)= - (rprimd(1,mu)*fcart(1,iatom)+& 1754 & rprimd(2,mu)*fcart(2,iatom)+& 1755 & rprimd(3,mu)*fcart(3,iatom)) 1756 end do 1757 end do 1758 1759 !Previous version 1760 !do iatom=1,natom 1761 !do mu=1,3 1762 !gred(mu,iatom)= - (rprimd(mu,1)*fcart(1,iatom)+& 1763 !& rprimd(mu,2)*fcart(2,iatom)+& 1764 !& rprimd(mu,3)*fcart(3,iatom)) 1765 !end do 1766 !end do 1767 1768 end subroutine fcart2gred
m_geometry/fixsym [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
fixsym
FUNCTION
Using input indsym which tells which atoms are related by symmetry, check that iatfix consistently fixes (freezes) all atoms which are related by symmetry, i.e. that iatfix does not break symmetry.
INPUTS
iatfix(3,natom)=integer array with 1 in every position for which the atom is to be kept fixed NOTE that this is not the input data structure for iatfix but it is the internal data structure used through most of the subroutines indsym(4,nsym,natom)=indirect indexing array for symmetrically related atoms; 4th element is label of symmetrically related atom natom=number of atoms nsym=number of symmetries (should be > 1 when this is called)
OUTPUT
(only checking) NOTE Stops execution with an error message if iatfix breaks symmetry.
SOURCE
1109 subroutine fixsym(iatfix,indsym,natom,nsym) 1110 1111 !Arguments ------------------------------------ 1112 !scalars 1113 integer,intent(in) :: natom,nsym 1114 !arrays 1115 integer,intent(in) :: iatfix(3,natom),indsym(4,nsym,natom) 1116 1117 !Local variables------------------------------- 1118 !scalars 1119 integer :: iatom,isym,jatom 1120 character(len=500) :: msg 1121 1122 ! ************************************************************************* 1123 1124 if (nsym > 1) then 1125 do iatom=1,natom 1126 do isym=1,nsym 1127 ! jatom is the label of a symmetrically related atom 1128 jatom=indsym(4,isym,iatom) 1129 ! Thus the atoms jatom and iatom must be fixed along the same directions 1130 if (iatfix(1,jatom) /= iatfix(1,iatom) .or. & 1131 iatfix(2,jatom) /= iatfix(2,iatom) .or. & 1132 iatfix(3,jatom) /= iatfix(3,iatom)) then 1133 write(msg, '(a,i0,a,i0,7a)' )& 1134 'Atom number: ',jatom,' is symmetrically equivalent to atom number: ',iatom,',',ch10,& 1135 'but according to iatfix, iatfixx, iatfixy and iatfixz, they',ch10,& 1136 'are not fixed along the same directions, which is forbidden.',ch10,& 1137 'Action: modify either the symmetry or iatfix(x,y,z) and resubmit.' 1138 ABI_ERROR(msg) 1139 end if 1140 end do 1141 end do 1142 end if 1143 1144 end subroutine fixsym
m_geometry/getspinrot [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
getspinrot
FUNCTION
From the symmetry matrix symrel expressed in the coordinate system rprimd, compute the components of the spinor rotation matrix
INPUTS
rprimd(3,3)=dimensional primitive translations for real space (bohr) symrel(3,3)=symmetry operation in real space in terms of primitive translations rprimd
OUTPUT
spinrot(4)=components of the spinor rotation matrix: spinrot(1)=$\cos \phi / 2$ spinrot(2)=$\sin \phi / 2 \times u_x$ spinrot(3)=$\sin \phi / 2 \times u_y$ spinrot(4)=$\sin \phi / 2 \times u_z$ where $\phi$ is the angle of rotation, and $(u_x,u_y,u_z)$ is the normalized direction of the rotation axis
NOTES
Only the proper part of the symmetry operation is taken into account: pure rotations, while the inversion part is taken away, if present. The whole collection of symmetry matrices is call symrel(3,3,nsym) symrel1 contains just one of those matrices symrel1(3,3)
SOURCE
782 subroutine getspinrot(rprimd, spinrot, symrel) 783 784 !Arguments ------------------------------------ 785 !arrays 786 integer,intent(in) :: symrel(3,3) 787 real(dp),intent(in) :: rprimd(3,3) 788 real(dp),intent(out) :: spinrot(4) 789 790 !Local variables------------------------------- 791 !scalars 792 integer :: det,ii 793 real(dp) :: cos_phi,norminv,phi,scprod,sin_phi 794 !character(len=500) :: msg 795 !arrays 796 integer :: identity(3,3),symrel1(3,3) 797 real(dp) :: axis(3),coord(3,3),coordinvt(3,3),matr1(3,3),matr2(3,3) 798 real(dp) :: rprimd_invt(3,3),vecta(3),vectb(3),vectc(3) 799 800 !************************************************************************** 801 802 symrel1(:,:) = symrel(:,:) 803 804 ! Compute determinant of the matrix 805 call mati3det(symrel1, det) 806 807 ! Produce a rotation from an improper symmetry 808 if (det==-1) symrel1(:,:) = -symrel1(:,:) 809 810 ! Test the possibility of the unit matrix 811 identity(:,:)=0; identity(1,1)=1; identity(2,2)=1; identity(3,3)=1 812 813 if (sum((symrel1(:,:) - identity(:,:))**2)/=0) then 814 815 ! Transform symmetry matrix in the system defined by rprimd 816 call matr3inv(rprimd, rprimd_invt) 817 do ii=1,3 818 coord(:,ii)=rprimd_invt(ii,:) 819 end do 820 call matr3inv(coord,coordinvt) 821 do ii=1,3 822 matr1(:,ii) = symrel1(:,1)*coord(1,ii) + symrel1(:,2)*coord(2,ii) + symrel1(:,3)*coord(3,ii) 823 end do 824 do ii=1,3 825 matr2(:,ii) = coordinvt(1,:)*matr1(1,ii) + coordinvt(2,:)*matr1(2,ii) + coordinvt(3,:)*matr1(3,ii) 826 end do 827 828 ! Find the eigenvector with unit eigenvalue of the 829 ! rotation matrix in cartesian coordinate, matr2 830 831 matr1(:,:)=matr2(:,:) 832 matr1(1,1)=matr1(1,1)-one 833 matr1(2,2)=matr1(2,2)-one 834 matr1(3,3)=matr1(3,3)-one 835 836 ! Compute the axis of rotation and the cos and sin of rotation angle 837 if(matr1(1,1)**2 + matr1(2,1)**2 + matr1(3,1)**2 < tol8 )then 838 ! The first direction is the axis 839 axis(1)=one ; axis(2)=zero ; axis(3)=zero 840 cos_phi=matr2(2,2) 841 sin_phi=matr2(3,2) 842 else if(matr1(1,2)**2 + matr1(2,2)**2 + matr1(3,2)**2 < tol8 )then 843 ! The second direction is the axis 844 axis(1)=zero ; axis(2)=one ; axis(3)=zero 845 cos_phi=matr2(3,3) 846 sin_phi=matr2(1,3) 847 else 848 ! In this case, try use the first and second vector to build the 849 ! rotation axis: compute their cross product 850 axis(1)=matr1(2,1)*matr1(3,2)-matr1(2,2)*matr1(3,1) 851 axis(2)=matr1(3,1)*matr1(1,2)-matr1(3,2)*matr1(1,1) 852 axis(3)=matr1(1,1)*matr1(2,2)-matr1(1,2)*matr1(2,1) 853 ! Then, try to normalize it 854 scprod=sum(axis(:)**2) 855 if(scprod<tol8)then 856 ! The first and second vectors were linearly dependent 857 ! Thus, use the first and third vectors 858 axis(1)=matr1(2,1)*matr1(3,3)-matr1(2,3)*matr1(3,1) 859 axis(2)=matr1(3,1)*matr1(1,3)-matr1(3,3)*matr1(1,1) 860 axis(3)=matr1(1,1)*matr1(2,3)-matr1(1,3)*matr1(2,1) 861 ! Normalize the vector 862 scprod=sum(axis(:)**2) 863 if(scprod < tol8)then 864 ABI_BUG('Cannot find the rotation axis.') 865 end if 866 end if 867 norminv=one/sqrt(scprod) 868 axis(:)=axis(:)*norminv 869 870 ! Project the axis vector out of the first unit vector, 871 ! and renormalize the projected vector 872 ! (the first vector cannot be the axis, as tested before) 873 vecta(1)=one-axis(1)**2 874 vecta(2)=-axis(1)*axis(2) 875 vecta(3)=-axis(1)*axis(3) 876 scprod=sum(vecta(:)**2) 877 norminv=one/sqrt(scprod) 878 vecta(:)=vecta(:)*norminv 879 ! Rotate the vector A, to get vector B 880 vectb(:)=matr2(:,1)*vecta(1)+matr2(:,2)*vecta(2)+matr2(:,3)*vecta(3) 881 ! Get dot product of vectors A and B, giving cos of the rotation angle 882 cos_phi=sum(vecta(:)*vectb(:)) 883 ! Compute the cross product of the axis and vector A 884 vectc(1)=axis(2)*vecta(3)-axis(3)*vecta(2) 885 vectc(2)=axis(3)*vecta(1)-axis(1)*vecta(3) 886 vectc(3)=axis(1)*vecta(2)-axis(2)*vecta(1) 887 ! Get dot product of vectors B and C, giving sin of the rotation angle 888 sin_phi=sum(vectb(:)*vectc(:)) 889 end if 890 891 ! Get the rotation angle, then the parameters of the spinor rotation 892 ! Here, treat possible inaccurate values of the cosine of phi 893 if(cos_phi> one-tol8 )cos_phi= one-tol8 894 if(cos_phi<-(one-tol8))cos_phi=-(one-tol8) 895 phi=acos(cos_phi) 896 if(sin_phi<zero)phi=-phi 897 ! Rectify the angle, such that its absolute values corresponds to 180, 120, 90, 60, or 0 degrees 898 phi=(nint(six*phi/pi))/six*pi 899 ! Compute components of the spinor matrix 900 spinrot(1)=cos(half*phi) 901 spinrot(2)=axis(1)*sin(half*phi) 902 spinrot(3)=axis(2)*sin(half*phi) 903 spinrot(4)=axis(3)*sin(half*phi) 904 905 else 906 907 ! Here, the case of the unit matrix 908 axis(:)=zero 909 phi=zero 910 spinrot(1)=one 911 spinrot(2)=zero 912 spinrot(3)=zero 913 spinrot(4)=zero 914 915 end if ! the case of the identity matrix 916 917 !DEBUG 918 !write(std_out,*)' getspinrot :' 919 !write(std_out,*)' symre =',symrel(:,:) 920 !write(std_out,*)' symrel1 =',symrel1(:,:) 921 !write(std_out,*)' rprimd =',rprimd(:,:) 922 !write(std_out,*)' matr2 =',matr2(:,:) 923 !write(std_out,*)' matr1 =',matr1(:,:) 924 !write(std_out,*)' phi (degree)=',phi*180._dp/pi 925 !write(std_out,'(a,3d16.6)' )' axis=',axis(:) 926 !write(std_out,*)' vecta=',vecta(:) 927 !stop 928 !ENDDEBUG 929 930 end subroutine getspinrot
m_geometry/gred2fcart [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
gred2fcart
FUNCTION
Convert reduced forces into cartesian forces
INPUTS
gred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
natom=Number of atoms in the unitary cell
Favgz_null=TRUE if the average cartesian force has to be set to zero
FALSE if it is set to zero only in x,y directions (not z)
gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
OUTPUT
fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
NOTES
Unlike gred, fcart has been corrected by enforcing
the translational symmetry, namely that the sum of force
on all atoms is zero (except is a slab is used)
SOURCE
1670 subroutine gred2fcart(favg,Favgz_null,fcart,gred,gprimd,natom) 1671 1672 !Arguments ------------------------------------ 1673 !scalars 1674 integer,intent(in) :: natom 1675 logical :: Favgz_null 1676 !arrays 1677 real(dp),intent(out) :: fcart(3,natom) 1678 real(dp),intent(in) :: gred(3,natom) 1679 real(dp),intent(in) :: gprimd(3,3) 1680 real(dp),intent(out) :: favg(3) 1681 1682 !Local variables------------------------------- 1683 !scalars 1684 integer :: iatom,mu 1685 1686 ! ************************************************************************* 1687 1688 !Note conversion to cartesian coordinates (bohr) AND 1689 !negation to make a force out of a gradient 1690 favg(:)=zero 1691 do iatom=1,natom 1692 do mu=1,3 1693 fcart(mu,iatom)= - (gprimd(mu,1)*gred(1,iatom)+& 1694 & gprimd(mu,2)*gred(2,iatom)+& 1695 & gprimd(mu,3)*gred(3,iatom)) 1696 favg(mu)=favg(mu)+fcart(mu,iatom) 1697 end do 1698 end do 1699 1700 !Subtract off average force from each force component 1701 !to avoid spurious drifting of atoms across cell. 1702 favg(:)=favg(:)/dble(natom) 1703 if(.not.Favgz_null) favg(3)=zero 1704 do iatom=1,natom 1705 fcart(:,iatom)=fcart(:,iatom)-favg(:) 1706 end do 1707 1708 end subroutine gred2fcart
m_geometry/ioniondist [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
ioniondist
FUNCTION
Compute ion-ion distances
INPUTS
natom= number of atoms in unit cell rprimd(3,3)=dimensional real space primitive translations (bohr) xred(3,natom)=dimensionless reduced coordinates of atoms inm(natom,natom)=index (m,n) of the atom option= 1 output ion-ion distances / 2 output ordering of ion-ion distances / 3 output variables in varlist according to ion-ion distances * magnetic ordering magv magnetic ordering of atoms given als 1 and -1, if not given fm is assumed varlist=List of variables magv(natom)= magnetic ordering of atoms atp=atom on which the perturbation was done
OUTPUT
SOURCE
2475 subroutine ioniondist(natom,rprimd,xred,inm,option,varlist,magv,atp,prtvol) 2476 2477 !Arguments ------------------------------------ 2478 !scalars 2479 integer,intent(in) :: natom,option 2480 integer,intent(in),optional :: atp !atom on which the perturbation was done 2481 !arrays 2482 real(dp),intent(in) :: rprimd(3,3) 2483 real(dp),intent(in) :: xred(3,natom) 2484 real(dp),intent(out) :: inm(natom,natom) 2485 integer,intent(in),optional :: magv(natom) 2486 real(dp),intent(in),optional :: varlist(natom) 2487 integer,intent(in),optional :: prtvol 2488 2489 !Local variables------------------------------- 2490 !scalars 2491 integer :: iatom,jatom,katom,kdum,atpp,prtvoll 2492 !character(len=500) :: msg 2493 !arrays 2494 integer :: interq(natom) 2495 real(dp) :: hxcart(3,natom),distm(natom,natom) 2496 real(dp) :: magvv(natom) 2497 2498 ! ************************************************************************* 2499 2500 hxcart=matmul(rprimd,xred) 2501 interq=(/(iatom,iatom=1,natom)/) 2502 inm=0 2503 2504 if (present(magv)) then 2505 magvv=magv 2506 else 2507 magvv=(/ (1, iatom=1,natom) /) 2508 end if 2509 2510 if (present(atp)) then 2511 atpp=atp 2512 else 2513 atpp=1 2514 end if 2515 2516 if (present(prtvol)) then 2517 prtvoll=prtvol 2518 else 2519 prtvoll=1 2520 end if 2521 2522 if (option==3.and.(.not.present(varlist))) then 2523 call wrtout(std_out,'ioniondist error: option=3 but no variable list provided for symmetrization','COLL') 2524 return 2525 end if 2526 2527 !call wrtout(std_out,' ioniondist start ','COLL') 2528 2529 distm=0 2530 katom=atpp-1 2531 do iatom=1,natom 2532 katom=katom+1 2533 if (katom > natom) katom=1 2534 distm(iatom,iatom)=0 2535 do jatom=iatom,natom 2536 distm(iatom,jatom)=dist2(xred(:,katom),xred(:,jatom),rprimd,1)*magvv(katom)*magvv(jatom) 2537 distm(jatom,iatom)=distm(iatom,jatom) 2538 end do 2539 end do 2540 2541 if (prtvoll>=3) then 2542 call wrtout(std_out,'ioniondist: ionic distances:','COLL') 2543 call prmat(distm,natom,natom,natom,std_out) 2544 end if 2545 2546 distm=anint(distm*10000_dp)/10000_dp ! rounding needed else distm(iatom,jatom)/= distm(1,kdum) sometimes fails 2547 2548 do iatom=1,natom 2549 if (option==1) then 2550 inm(iatom,:)=distm(iatom,:) 2551 else 2552 do jatom=iatom,natom 2553 kdum=1 2554 do while ( (kdum <= natom) .and. (distm(iatom,jatom)/= distm(1,kdum)) ) 2555 kdum=kdum+1 2556 end do 2557 if (option==2) then 2558 inm(iatom,jatom)=interq(kdum) 2559 else if (option==3) then 2560 inm(iatom,jatom)=varlist(kdum) 2561 end if 2562 inm(jatom,iatom)=inm(iatom,jatom) 2563 end do 2564 end if 2565 end do 2566 2567 if (prtvoll==2) then 2568 call wrtout(std_out,'ioniondist: symmetrized matrix:','COLL') 2569 call prmat(distm,1,natom,natom,std_out) 2570 else if (prtvoll>=3) then 2571 call wrtout(std_out,'ioniondist: symmetrized matrix:','COLL') 2572 call prmat(distm,natom,natom,natom,std_out) 2573 end if 2574 2575 end subroutine ioniondist
m_geometry/littlegroup_pert [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
littlegroup_pert
FUNCTION
If syuse==0 and abs(rfmeth)==2, determines the set of symmetries that leaves a perturbation invariant. (Actually, all symmetries that leaves a q-wavevector invariant should be used to reduce the number of k-points for all perturbations. Unfortunately, one has to take into account the sign reversal of the perturbation under the symmetry operations, which makes GS routines not usable for the respfn code. The intermediate choice was to select only those that keep also the perturbation invariant. Note that the wavevector of the perturbation must also be invariant, a translation vector in real space is NOT allowed ).
INPUTS
gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1) idir=direction of the perturbation indsym(4,nsym,natom)=indirect indexing of atom labels--see subroutine symatm for definition (if nsym>1) iout=if non-zero, output on unit iout ipert=characteristics of the perturbation natom= number of atoms nsym=number of space group symmetries rfmeth = 1 or -1 if non-stationary block 2 or -2 if stationary block 3 or -3 if third order derivatives positive if symmetries are used to set elements to zero whenever possible, negative to prevent this to happen. symq(4,2,nsym)= Table computed by littlegroup_q. three first numbers define the G vector; fourth number is zero if the q-vector is not preserved, is 1 otherwise second index is one without time-reversal symmetry, two with time-reversal symmetry symafm(nsym)=(anti)ferromagnetic part of the symmetry operations symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space) symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space) syuse= flag to use the symmetries or not. If 0 usei it, if 1 do not use it. tnons(3,nsym)=nonsymmorphic translations of space group in terms of real space primitive translations (may be 0) [unit]=By default the routine writes to std_out and this is very annoying if we are inside a big loop. Use unit=dev_null or a negative integer to disable writing.
OUTPUT
nsym1 =number of space group symmetries that leaves the perturbation invariant symaf1(nsym1)=(anti)ferromagnetic part of the corresponding symmetry operations symrl1(3,3,nsym1)=corresponding 3x3 matrices of the group symmetries (real space) tnons1(3,nsym1)=corresponding nonsymmorphic translations of space group in terms of real space primitive translations (may be 0)!!
SOURCE
3297 subroutine littlegroup_pert(gprimd,idir,indsym,iout,ipert,natom,nsym,nsym1, & 3298 & rfmeth,symafm,symaf1,symq,symrec,symrel,symrl1,syuse,tnons,tnons1, & 3299 & unit) ! Optional 3300 3301 !Arguments ------------------------------- 3302 !scalars 3303 integer,intent(in) :: idir,iout,ipert,natom,nsym,rfmeth,syuse 3304 integer,intent(in),optional :: unit 3305 integer,intent(out) :: nsym1 3306 !arrays 3307 integer,intent(in) :: indsym(4,nsym,natom),symafm(nsym),symq(4,2,nsym) 3308 integer,intent(in) :: symrec(3,3,nsym),symrel(3,3,nsym) 3309 integer,intent(out) :: symaf1(nsym),symrl1(3,3,nsym) 3310 real(dp),intent(in) :: gprimd(3,3),tnons(3,nsym) 3311 real(dp),intent(out) :: tnons1(3,nsym) 3312 3313 !Local variables ------------------------- 3314 !scalars 3315 integer :: idir1,ii,istr,isym,jj,nsym_test,tok,ount 3316 character(len=500) :: msg 3317 !arrays 3318 integer :: sym_test(3,3,2) 3319 real(dp) :: str_test(6) 3320 3321 ! ********************************************************************* 3322 3323 ount = std_out; if (present(unit)) ount = unit 3324 3325 nsym1=0 3326 if((ipert==natom+3 .or. ipert==natom+4) .and. syuse==0 .and. abs(rfmeth)==2) then 3327 ! Strain perturbation section 3328 ! Use ground state routine which symmetrizes cartesian stress as a quick 3329 ! and dirty test for the invariance of the strain (ipert,idir) under 3330 ! each candidate symmetry 3331 ! I am presently assuming that translations are acceptable because I dont 3332 ! see why not. 3333 3334 istr=3*(ipert-natom-3)+idir 3335 nsym_test=2 3336 ! Store identity as first element for test 3337 sym_test(:,:,1)=0 3338 sym_test(1,1,1)=1; sym_test(2,2,1)=1; sym_test(3,3,1)=1 3339 do isym=1,nsym 3340 sym_test(:,:,2)=symrec(:,:,isym) 3341 str_test(:)=0.0_dp 3342 str_test(istr)=1.0_dp 3343 call stresssym(gprimd,nsym_test,str_test,sym_test) 3344 if(abs(str_test(istr)-1.0_dp)<tol8)then 3345 ! The test has been successful ! 3346 nsym1=nsym1+1 3347 symaf1(nsym1)=symafm(isym) 3348 do ii=1,3 3349 tnons1(ii,nsym1)=tnons(ii,isym) 3350 do jj=1,3 3351 symrl1(ii,jj,nsym1)=symrel(ii,jj,isym) 3352 end do 3353 end do 3354 end if 3355 end do 3356 3357 else if(ipert>natom .or. syuse/=0 .or. abs(rfmeth)/=2)then 3358 3359 ! Not yet coded for d/dk or electric field perturbations 3360 nsym1=1 3361 do ii=1,3 3362 tnons1(ii,1)=0._dp 3363 symaf1(1)=1 3364 do jj=1,3 3365 symrl1(ii,jj,1)=0 3366 if(ii==jj)symrl1(ii,jj,1)=1 3367 end do 3368 end do 3369 3370 else 3371 3372 do isym=1,nsym 3373 ! Check that the symmetry operation preserves the wavevector 3374 ! (a translation is NOT allowed) 3375 if(symq(4,1,isym)==1 .and.& 3376 & symq(1,1,isym)==0 .and.& 3377 & symq(2,1,isym)==0 .and.& 3378 & symq(3,1,isym)==0 )then 3379 ! Check that the symmetry operation preserves the atom 3380 if(ipert==indsym(4,isym,ipert))then 3381 ! Check if the direction is preserved 3382 tok=1 3383 do idir1=1,3 3384 if((idir1==idir.and.symrec(idir,idir1,isym)/=1) .or.& 3385 & (idir1/=idir.and.symrec(idir,idir1,isym)/=0))then 3386 tok=0 3387 end if 3388 end do 3389 if(tok==1)then 3390 ! All the tests have been successful ! 3391 nsym1=nsym1+1 3392 symaf1(nsym1)=symafm(isym) 3393 do ii=1,3 3394 tnons1(ii,nsym1)=tnons(ii,isym) 3395 do jj=1,3 3396 symrl1(ii,jj,nsym1)=symrel(ii,jj,isym) 3397 end do 3398 end do 3399 end if 3400 3401 end if 3402 end if 3403 end do 3404 end if 3405 3406 if (nsym1<1) then 3407 write(msg,'(a,i0,a)')' The number of selected symmetries should be > 0, while it is nsym= ',nsym1,'.' 3408 ABI_BUG(msg) 3409 end if 3410 3411 if (nsym1 /= 1) then 3412 if (iout /= ount .and. iout > 0) then 3413 write(msg,'(a,i5,a)')' Found ',nsym1,' symmetries that leave the perturbation invariant.' 3414 call wrtout(iout,msg,'COLL') 3415 end if 3416 write(msg,'(a,i5,a)')' littlegroup_pert: found ',nsym1,' symmetries that leave the perturbation invariant: ' 3417 call wrtout(ount,msg,'COLL') 3418 else 3419 if (iout /= ount .and. iout > 0) then 3420 write(msg,'(a,a)')' The set of symmetries contains',' only one element for this perturbation.' 3421 call wrtout(iout,msg,'COLL') 3422 end if 3423 write(msg,'(a)')' littlegroup_pert: only one element in the set of symmetries for this perturbation:' 3424 call wrtout(ount,msg,'COLL') 3425 end if 3426 3427 if (ount > 0) then 3428 do isym=1,nsym1 3429 write(msg, '(9i4)' )((symrl1(ii,jj,isym),ii=1,3),jj=1,3) 3430 call wrtout(ount,msg,'COLL') 3431 end do 3432 end if 3433 3434 end subroutine littlegroup_pert
m_geometry/metric [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
metric
FUNCTION
Compute first dimensional primitive translation vectors in reciprocal space gprimd from rprimd, and eventually writes out. Then, computes metrics for real and recip space rmet and gmet using length dimensional primitive translation vectors in columns of rprimd(3,3) and gprimd(3,3). gprimd is the inverse transpose of rprimd. i.e. $ rmet_{i,j}= \sum_k ( rprimd_{k,i}*rprimd_{k,j} ) $ $ gmet_{i,j}= \sum_k ( gprimd_{k,i}*gprimd_{k,j} ) $ Also computes unit cell volume ucvol in $\textrm{bohr}^3$
INPUTS
rprimd(3,3)=dimensional primitive translations for real space (bohr) iout=unit number of output file. If iout<0, do not write output.
OUTPUT
gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$). gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) rmet(3,3)=real space metric ($\textrm{bohr}^{2}$). ucvol=unit cell volume ($\textrm{bohr}^{3}$).
SOURCE
1197 subroutine metric(gmet, gprimd, iout, rmet, rprimd, ucvol) 1198 1199 !Arguments ------------------------------------ 1200 !scalars 1201 integer,intent(in) :: iout 1202 real(dp),intent(out) :: ucvol 1203 !arrays 1204 real(dp),intent(in) :: rprimd(3,3) 1205 real(dp),intent(out) :: gmet(3,3),gprimd(3,3),rmet(3,3) 1206 1207 !Local variables------------------------------- 1208 !scalars 1209 integer :: nu 1210 character(len=500) :: msg 1211 !arrays 1212 real(dp) :: angle(3) 1213 1214 ! ************************************************************************* 1215 1216 ! Compute unit cell volume 1217 ucvol=rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+& 1218 rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+& 1219 rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3)) 1220 !ucvol = det3r(rprimd) 1221 1222 ! Check that the input primitive translations are not linearly dependent (and none is zero); i.e. ucvol~=0 1223 ! Also ask that the mixed product is positive. 1224 if (abs(ucvol)<tol12) then 1225 !write(std_out,*)"rprimd",rprimd,"ucvol",ucvol 1226 write(msg,'(5a)')& 1227 'Input rprim and acell gives vanishing unit cell volume.',ch10,& 1228 'This indicates linear dependency between primitive lattice vectors',ch10,& 1229 'Action: correct either rprim or acell in input file.' 1230 ABI_ERROR(msg) 1231 end if 1232 if (ucvol<zero)then 1233 write(msg,'(2a,3(a,3es16.6,a),7a)')& 1234 'Current rprimd gives negative (R1 x R2) . R3 . ',ch10,& 1235 'Rprimd =',rprimd(:,1),ch10,& 1236 ' ',rprimd(:,2),ch10,& 1237 ' ',rprimd(:,3),ch10,& 1238 'Action: if the cell size and shape are fixed (optcell==0),',ch10,& 1239 ' exchange two of the input rprim vectors;',ch10,& 1240 ' if you are optimizing the cell size and shape (optcell/=0),',ch10,& 1241 ' maybe the move was too large, and you might try to decrease strprecon.' 1242 ABI_ERROR(msg) 1243 end if 1244 1245 ! Generate gprimd 1246 call matr3inv(rprimd,gprimd) 1247 1248 ! Write out rprimd, gprimd and ucvol 1249 if (iout>=0) then 1250 write(msg,'(2a)')' Real(R)+Recip(G) ','space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):' 1251 call wrtout(iout,msg) 1252 do nu=1,3 1253 write(msg, '(1x,a,i1,a,3f11.7,2x,a,i1,a,3f11.7)' ) & 1254 'R(',nu,')=',rprimd(:,nu)+tol10,& 1255 'G(',nu,')=',gprimd(:,nu)+tol10 1256 call wrtout(iout,msg) 1257 end do 1258 write(msg,'(a,1p,e15.7,a)') ' Unit cell volume ucvol=',ucvol+tol10,' bohr^3' 1259 call wrtout(iout,msg,'COLL') 1260 call wrtout(std_out,msg,'COLL') 1261 end if 1262 1263 ! Compute real space metric. 1264 rmet = MATMUL(TRANSPOSE(rprimd),rprimd) 1265 1266 ! Compute reciprocal space metric. 1267 gmet = MATMUL(TRANSPOSE(gprimd),gprimd) 1268 1269 ! Write out the angles 1270 if (iout>=0) then 1271 angle(1)=acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))/two_pi*360.0d0 1272 angle(2)=acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))/two_pi*360.0d0 1273 angle(3)=acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))/two_pi*360.0d0 1274 write(msg, '(a,3es16.8,a)' )' Angles (23,13,12)=',angle(1:3),' degrees' 1275 call wrtout(iout,msg,'COLL') 1276 call wrtout(std_out,msg,'COLL') 1277 end if 1278 1279 end subroutine metric
m_geometry/mkradim [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
mkradim
FUNCTION
Not so trivial subroutine to make dimensionless real space primitive translations rprim(3,3) from dimensional rprimd(3). also make acell(3).
INPUTS
rprimd(3,3)=dimensional real space primitive translations (bohr) where: rprimd(i,j)=rprim(i,j)*acell(j)
OUTPUT
acell(3)=unit cell length scales (bohr) rprim(3,3)=dimensionless real space primitive translations
SOURCE
1301 subroutine mkradim(acell,rprim,rprimd) 1302 1303 !Arguments ------------------------------------ 1304 !arrays 1305 real(dp),intent(out) :: acell(3),rprim(3,3) 1306 real(dp),intent(in) :: rprimd(3,3) 1307 1308 !Local variables------------------------------- 1309 !scalars 1310 integer :: ii,jj 1311 real(dp) :: rprim_maxabs 1312 1313 ! ************************************************************************* 1314 1315 !Use a representation based on normalised rprim vectors 1316 do ii=1,3 1317 acell(ii)=sqrt(rprimd(1,ii)**2+rprimd(2,ii)**2+rprimd(3,ii)**2) 1318 rprim(:,ii)=rprimd(:,ii)/acell(ii) 1319 end do 1320 1321 !Suppress meaningless values 1322 rprim_maxabs=maxval(abs(rprim)) 1323 do ii=1,3 1324 do jj=1,3 1325 if(abs(rprim(ii,jj))<tol12*rprim_maxabs)rprim(ii,jj)=zero 1326 enddo 1327 enddo 1328 1329 end subroutine mkradim
m_geometry/mkrdim [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
mkrdim
FUNCTION
Trivial subroutine to make dimensional real space primitive translations from length scales acell(3) and dimensionless translations rprim(3,3).
INPUTS
acell(3)=unit cell length scales (bohr) rprim(3,3)=dimensionless real space primitive translations
OUTPUT
rprimd(3,3)=dimensional real space primitive translations (bohr) where: rprimd(i,j)=rprim(i,j)*acell(j)
SOURCE
1528 subroutine mkrdim(acell,rprim,rprimd) 1529 1530 !Arguments ------------------------------------ 1531 !arrays 1532 real(dp),intent(in) :: acell(3),rprim(3,3) 1533 real(dp),intent(out) :: rprimd(3,3) 1534 1535 !Local variables------------------------------- 1536 !scalars 1537 integer :: ii,jj 1538 1539 ! ************************************************************************* 1540 1541 do ii=1,3 1542 do jj=1,3 1543 rprimd(ii,jj)=rprim(ii,jj)*acell(jj) 1544 end do 1545 end do 1546 1547 end subroutine mkrdim
m_geometry/normv_int_vector [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
normv_int_vector
FUNCTION
Returns the norm of an integer 3D vector expressed in reduced coordinates. either in real or reciprocal space. In the later case the factor 2pi has to be included, due to the conventions used in abinit to define the reciprocal lattice.
INPUTS
OUTPUT
SOURCE
159 function normv_int_vector(xv, met, space) result(res) 160 161 !Arguments ------------------------------------ 162 !scalars 163 real(dp) :: res 164 character(len=1),intent(in) :: space 165 !arrays 166 real(dp),intent(in) :: met(3,3) 167 integer,intent(in) :: xv(3) 168 169 ! ************************************************************************* 170 171 res = ( xv(1)*met(1,1)*xv(1) + xv(2)*met(2,2)*xv(2) + xv(3)*met(3,3)*xv(3) & 172 & +two*( xv(1)*met(1,2)*xv(2) + xv(1)*met(1,3)*xv(3) + xv(2)*met(2,3)*xv(3)) ) 173 174 select case (space) 175 case ('r','R') 176 res=SQRT(res) 177 case ('g','G') 178 res=two_pi*SQRT(res) 179 case default 180 ABI_BUG('Wrong value for space') 181 end select 182 183 end function normv_int_vector
m_geometry/normv_int_vector_array [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
normv_int_vector_array
FUNCTION
Returns the norm of an array of integer 3D vectors expressed in reduced coordinates. either in real or reciprocal space. In the later case the factor 2pi has to be included, due to the conventions used in abinit to define the reciprocal lattice.
INPUTS
OUTPUT
SOURCE
203 function normv_int_vector_array(xv,met,space) result(res) 204 205 !Arguments ------------------------------------ 206 !scalars 207 character(len=1),intent(in) :: space 208 !arrays 209 real(dp),intent(in) :: met(3,3) 210 integer,intent(in) :: xv(:,:) 211 !this awful trick is needed to avoid problems with abilint 212 real(dp) :: res(SIZE(xv(1,:))) 213 !real(dp) :: res(SIZE(xv,DIM=2)) 214 215 ! ************************************************************************* 216 217 res(:) = ( xv(1,:)*met(1,1)*xv(1,:) + xv(2,:)*met(2,2)*xv(2,:) + xv(3,:)*met(3,3)*xv(3,:) & 218 +two*(xv(1,:)*met(1,2)*xv(2,:) + xv(1,:)*met(1,3)*xv(3,:) + xv(2,:)*met(2,3)*xv(3,:)) ) 219 220 select case (space) 221 case ('r','R') 222 res(:)=SQRT(res(:)) 223 case ('g','G') 224 res(:)=two_pi*SQRT(res(:)) 225 case default 226 ABI_BUG('Wrong value for space') 227 end select 228 229 end function normv_int_vector_array
m_geometry/normv_rdp_vector [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
normv_rdp_vector
FUNCTION
Compute the norm of a vector expressed in reduced coordinates using the metric met. The result is multiplied by 2pi in case of a vector in reciprocal space to take into account the correct normalisation of the reciprocal lattice vectors
INPUTS
xv(3)=Vector in reduced coordinates met(3,3)=Metric tensor space=Character defining whether we are working in real (r|R) or reciprocal space (g|G)
OUTPUT
normv_rdp_vector=norm of xv
NOTES
The routine is able to deal both with a single vector as well as arrays of vectors. Versions for integer and real vectors are provided.
SOURCE
116 function normv_rdp_vector(xv,met,space) result(res) 117 118 !Arguments ------------------------------------ 119 !scalars 120 real(dp) :: res 121 character(len=1),intent(in) :: space 122 !arrays 123 real(dp),intent(in) :: met(3,3),xv(3) 124 125 ! ************************************************************************* 126 127 res = (xv(1)*met(1,1)*xv(1) + xv(2)*met(2,2)*xv(2) + xv(3)*met(3,3)*xv(3) & 128 & +two*(xv(1)*met(1,2)*xv(2) + xv(1)*met(1,3)*xv(3) + xv(2)*met(2,3)*xv(3)) ) 129 130 select case (space) 131 case ('r','R') 132 res=SQRT(res) 133 case ('g','G') 134 res=two_pi*SQRT(res) 135 case default 136 ABI_BUG('Wrong value for space') 137 end select 138 139 end function normv_rdp_vector
m_geometry/normv_rdp_vector_array [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
normv_rdp_vector_array
FUNCTION
Returns the norm of an array of real 3D vectors expressed in reduced coordinates. either in real or reciprocal space. In the later case the factor 2pi has to be included, due to the conventions used in abinit to define the reciprocal lattice.
INPUTS
OUTPUT
SOURCE
249 function normv_rdp_vector_array(xv,met,space) result(res) 250 251 !Arguments ------------------------------------ 252 !scalars 253 character(len=1),intent(in) :: space 254 !arrays 255 real(dp),intent(in) :: met(3,3) 256 real(dp),intent(in) :: xv(:,:) 257 !this awful trick is needed to avoid problems with abilint 258 real(dp) :: res(SIZE(xv(1,:))) 259 !real(dp) :: res(SIZE(xv,DIM=2)) 260 261 ! ************************************************************************* 262 263 res(:) = ( xv(1,:)*met(1,1)*xv(1,:) + xv(2,:)*met(2,2)*xv(2,:) + xv(3,:)*met(3,3)*xv(3,:) & 264 & +two*(xv(1,:)*met(1,2)*xv(2,:) + xv(1,:)*met(1,3)*xv(3,:) + xv(2,:)*met(2,3)*xv(3,:)) ) 265 266 select case (space) 267 case ('r','R') 268 res(:)=SQRT(res(:)) 269 case ('g','G') 270 res(:)=two_pi*SQRT(res) 271 case default 272 ABI_BUG('Wrong value for space') 273 end select 274 275 end function normv_rdp_vector_array
m_geometry/phdispl_cart2red [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
phdispl_cart2red
FUNCTION
Calculates the displacement vectors for all branches in reduced coordinates. $ displ_red = displ_cart \cdot gprimd $ for each phonon branch.
INPUTS
natom=Number of atoms. gprimd(3,3)=Dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) displ_cart(2,3*natom,3*natom)=Phonon displacement in Cartesian coordinates.
OUTPUT
displ_red(2,3*natom,3*natom)=Phonon displacement in reduded coordinates.
SOURCE
703 subroutine phdispl_cart2red(natom, gprimd, displ_cart, displ_red) 704 705 !Arguments ------------------------------------ 706 !scalars 707 integer,intent(in) :: natom 708 !arrays 709 real(dp),intent(in) :: gprimd(3,3) 710 real(dp),intent(in) :: displ_cart(2,3*natom,3*natom) 711 real(dp),intent(out) :: displ_red(2,3*natom,3*natom) 712 713 !Local variables------------------------- 714 !scalars 715 integer :: nbranch,jbranch,iatom,idir,ibranch,kdir,k1 716 717 ! ************************************************************************* 718 719 displ_red = zero 720 721 nbranch=3*natom 722 723 do jbranch=1,nbranch 724 ! 725 do iatom=1,natom 726 do idir=1,3 727 ibranch=idir+3*(iatom-1) 728 do kdir=1,3 729 k1 = kdir+3*(iatom-1) 730 ! WARNING: could be non-transpose of rprimd matrix : to be checked. 731 ! 23 june 2004: rprimd becomes gprimd 732 ! could be gprim and then multiply by acell... 733 ! Nope, checked and ok with gprimd 24 jun 2004 734 displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + & 735 & gprimd(kdir,idir) * displ_cart(1,k1,jbranch) 736 737 displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + & 738 & gprimd(kdir,idir) * displ_cart(2,k1,jbranch) 739 740 end do !kdir 741 end do !idir 742 end do !iatom 743 ! 744 end do !jbranch 745 746 end subroutine phdispl_cart2red
m_geometry/randomcellpos [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
randomcellpos
FUNCTION
This subroutine creates a unit cell with random atomic positions. It is assumed that the cell parameters are given and fixed. Several methods are used to generate the cell.
INPUTS
natom=number of atoms
npsp=number of pseudopotentials (needed for the dimension of znucl)
ntypat=number of type of atoms
random_atpos=input variable
0 no generation of random atomic potision
1 completely random atomic potisions
2 random atomic positions, avoiding too close atoms
(prevent coming closer than a fraction of the sum of covalent radii)
3 same than 2 but also generates the rprim and acell randomly
within some given ranges (angles between 50 and 130)
ratsph(1:ntypat)=radius of the atomic sphere
rprimd(3,3)=dimensional primitive translations in real space (bohr)
typat(1:natom)= input variable giving the type of each atom
znucl(1:npsp)=nuclear number of atom as specified in psp file
OUTPUT
xred(3,natom)=reduced dimensionless atomic coordinates
SIDE EFFECTS
NOTES
SOURCE
2130 subroutine randomcellpos(natom,npsp,ntypat,random_atpos,ratsph,rprim,rprimd,typat,xred,znucl,acell) 2131 2132 !Arguments ------------------------------------ 2133 !scalars 2134 integer,intent(in) :: natom,npsp,ntypat,random_atpos 2135 !arrays 2136 integer, intent(in) :: typat(natom) 2137 real(dp),intent(in) :: ratsph(ntypat) 2138 real(dp), intent(inout) :: rprim(3,3) 2139 real(dp), intent(inout) :: rprimd(3,3) 2140 real(dp), intent(inout) :: xred(3,natom) 2141 real(dp), intent(in) :: znucl(npsp) 2142 real(dp), intent(inout) :: acell(3) 2143 2144 !Local variables------------------------------- 2145 integer :: iatom=0,ii,idum=-20 2146 real(dp) :: rij(3), rijd(3), radiuscovi, radiuscovj, dist, rati, ratj, angdeg(3) 2147 real(dp) :: cosang,aa,cc,a2 2148 character(len=500) :: msg 2149 type(atomdata_t) :: atom 2150 2151 ! ************************************************************************* 2152 2153 !DEBUG 2154 !For the time being, print rprimd to keep it as an argument, in spite of abirule checking. 2155 !write (std_out,*) ' randomcellpos : enter' 2156 !write(std_out,*)' rprimd=',rprimd 2157 !write(std_out,*)' znucl=',znucl 2158 !write(std_out,*)' typat=',typat 2159 !write(std_out,*)' random_atpos=',random_atpos 2160 !ENDDEBUG 2161 2162 if(random_atpos==2 .and. npsp/=ntypat)then 2163 write(msg, '(a,i5,2a,i5,a,i5,4a)' )& 2164 & 'Input variable random_atpos= ',random_atpos,ch10,& 2165 & 'However, the number of pseudopotentials ',npsp,', is not equal to the number of type of atoms ',ntypat,ch10,& 2166 & 'The use of alchemical mixing cannot be combined with the constraint based on the mixing of covalent radii.',ch10,& 2167 & 'Action: switch to another value of random_atpos.' 2168 ABI_ERROR(msg) 2169 end if 2170 2171 !random_atpos = 0 Default value, no random initialisation 2172 !random_atpos = 1 Fully random (Is it really useful ???) 2173 !random_atpos = 2 Random, but the sum of the two covalent radii is 2174 !less than the interatomic distance 2175 !random_atpos = 3 Random, but the sum of the two (other type of) 2176 !radii is less than the interatomic distance 2177 !random_atpos = 4 Random, but the sum of the two pseudopotential 2178 !radii is less than the interatomic distance 2179 !random_atpos = 5 Random, but the interatomic distance must be bigger 2180 !than the sum of 2181 !some input variable (well, instead of defining a new variable, why 2182 !not use ratsph ?) 2183 !Right now we are not using a factor for the tested distance.. something to be done, after a new variable has been defined 2184 2185 if (random_atpos /= 0) then 2186 select case (random_atpos) 2187 case (1) 2188 do ii=1,natom 2189 xred(1,ii)=uniformrandom(idum) 2190 xred(2,ii)=uniformrandom(idum) 2191 xred(3,ii)=uniformrandom(idum) 2192 end do 2193 case (2) 2194 iatom=0 2195 do 2196 iatom=iatom+1 2197 xred(1,iatom)=uniformrandom(idum) 2198 xred(2,iatom)=uniformrandom(idum) 2199 xred(3,iatom)=uniformrandom(idum) 2200 call atomdata_from_znucl(atom,znucl(typat(iatom))) 2201 radiuscovi = atom%rcov 2202 do ii=1,iatom-1 2203 rij=xred(:,iatom)-xred(:,ii) 2204 ! periodic boundary conditions 2205 rij = rij - 0.5 2206 rij = rij - anint (rij) 2207 ! coming back to cube between (0,1) 2208 rij = rij + 0.5 2209 ! convert reduced coordinates to cartesian coordinates 2210 call xred2xcart(1,rprimd,rijd,rij) 2211 dist=dot_product(rijd,rijd) 2212 call atomdata_from_znucl(atom,znucl(typat(ii))) 2213 radiuscovj = atom%rcov 2214 if (dist<(radiuscovj+radiuscovi)) then 2215 iatom = iatom -1 2216 EXIT 2217 end if 2218 end do 2219 if (iatom>=natom) EXIT 2220 end do 2221 case(3) 2222 iatom=0 2223 do 2224 iatom=iatom+1 2225 xred(1,iatom)=uniformrandom(idum) 2226 xred(2,iatom)=uniformrandom(idum) 2227 xred(3,iatom)=uniformrandom(idum) 2228 call atomdata_from_znucl(atom,znucl(typat(iatom))) 2229 radiuscovi = atom%rcov 2230 do ii=1,iatom-1 2231 rij=xred(:,iatom)-xred(:,ii) 2232 ! periodic boundary conditions 2233 rij = rij - 0.5 2234 rij = rij - anint (rij) 2235 ! coming back to cube between (0,1) 2236 rij = rij + 0.5 2237 ! convert reduced coordinates to cartesian coordinates 2238 call xred2xcart(1,rprimd,rijd,rij) 2239 dist=dot_product(rijd,rijd) 2240 call atomdata_from_znucl(atom,znucl(typat(ii))) 2241 radiuscovj = atom%rcov 2242 if (dist<(radiuscovj+radiuscovi)) then 2243 iatom = iatom -1 2244 EXIT 2245 end if 2246 end do 2247 if (iatom>=natom) EXIT 2248 end do 2249 do ii=1,3 2250 ! generates cells with angles between 60 and 120 degrees 2251 angdeg(ii)=60_dp+uniformrandom(idum)*60.0_dp 2252 end do 2253 if (angdeg(1)+angdeg(2)+angdeg(3)>360._dp) then 2254 angdeg(3)=360._dp-angdeg(1)-angdeg(2) 2255 end if 2256 ! check if angles are between the limits and create rprim 2257 if( abs(angdeg(1)-angdeg(2))<tol12 .and. & 2258 & abs(angdeg(2)-angdeg(3))<tol12 .and. & 2259 & abs(angdeg(1)-90._dp)+abs(angdeg(2)-90._dp)+abs(angdeg(3)-90._dp)>tol12 )then 2260 ! Treat the case of equal angles (except all right angles) : 2261 ! generates trigonal symmetry wrt third axis 2262 cosang=cos(pi*angdeg(1)/180.0_dp) 2263 a2=2.0_dp/3.0_dp*(1.0_dp-cosang) 2264 aa=sqrt(a2) 2265 cc=sqrt(1.0_dp-a2) 2266 rprim(1,1)=aa ; rprim(2,1)=0.0_dp ; rprim(3,1)=cc 2267 rprim(1,2)=-0.5_dp*aa ; rprim(2,2)= sqrt(3.0_dp)*0.5_dp*aa ; rprim(3,2)=cc 2268 rprim(1,3)=-0.5_dp*aa ; rprim(2,3)=-sqrt(3.0_dp)*0.5_dp*aa ; rprim(3,3)=cc 2269 ! DEBUG 2270 ! write(std_out,*)' ingeo : angdeg=',angdeg(1:3) 2271 ! write(std_out,*)' ingeo : aa,cc=',aa,cc 2272 ! ENDDEBUG 2273 else 2274 ! Treat all the other cases 2275 rprim(:,:)=0.0_dp 2276 rprim(1,1)=1.0_dp 2277 rprim(1,2)=cos(pi*angdeg(3)/180.0_dp) 2278 rprim(2,2)=sin(pi*angdeg(3)/180.0_dp) 2279 rprim(1,3)=cos(pi*angdeg(2)/180.0_dp) 2280 rprim(2,3)=(cos(pi*angdeg(1)/180.0_dp)-rprim(1,2)*rprim(1,3))/rprim(2,2) 2281 rprim(3,3)=sqrt(1.0_dp-rprim(1,3)**2-rprim(2,3)**2) 2282 end if 2283 ! generate acell 2284 aa=zero 2285 do ii=1,npsp 2286 aa=znucl(ii) 2287 end do 2288 do ii=1,3 2289 acell(ii)=aa+uniformrandom(idum)*4.0 2290 end do 2291 call mkrdim(acell,rprim,rprimd) 2292 case(4) 2293 write(std_out,*) 'Not implemented yet' 2294 case(5) 2295 iatom=0 2296 do 2297 iatom=iatom+1 2298 xred(1,iatom)=uniformrandom(idum) 2299 xred(2,iatom)=uniformrandom(idum) 2300 xred(3,iatom)=uniformrandom(idum) 2301 rati=ratsph(typat(iatom)) 2302 do ii=1,iatom-1 2303 ratj=ratsph(typat(ii)) 2304 ! apply periodic boundary conditions 2305 rij=(xred(:,iatom)-xred(:,ii))-0.5 2306 rij = rij - ANINT ( rij ) 2307 rij = rij + 0.5 2308 call xred2xcart(natom,rprimd,rijd,rij) 2309 dist=dot_product(rijd,rijd) 2310 if (dist<(rati+ratj)) EXIT 2311 end do 2312 if (iatom==natom) EXIT 2313 if (ii<(iatom-1)) iatom=iatom-1 2314 end do 2315 end select 2316 end if 2317 2318 end subroutine randomcellpos
m_geometry/remove_inversion [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
remove_inversion
FUNCTION
Remove the inversion symmetry from a symmetry set as well all the improper rotations (if present)
INPUTS
nsym=initial number of symmetries symrel(3,3,nsym)=Initial set of symmetry operarations in real space tnons(3,nsym)=Initial fractional translations
OUTPUT
nsym_out=Number of symmetries in the set without improper rotation symrel_out(:,:) [pointer] = output symmetries without improper rotations tnons_out(:) [pointer] = fractional translations associated to symrel_out pinv=-1 if the inversion has been removed, 1 otherwise
NOTES
Note the use of pointers, memory is allocated inside the procedure and passed back to the caller. Thus memory deallocation is relegated to the caller. To be on the safe side the pointers should be nullified before entering.
SOURCE
2709 subroutine remove_inversion(nsym,symrel,tnons,nsym_out,symrel_out,tnons_out,pinv) 2710 2711 !Arguments ------------------------------------ 2712 !scalars 2713 integer,intent(in) :: nsym 2714 integer,intent(out) :: nsym_out,pinv 2715 !arrays 2716 integer,intent(in) :: symrel(3,3,nsym) 2717 integer,pointer :: symrel_out(:,:,:) 2718 real(dp),intent(in) :: tnons(3,nsym) 2719 real(dp),pointer :: tnons_out(:,:) 2720 2721 !Local variables------------------------------- 2722 !scalars 2723 integer :: is,is2,is_discarded,is_inv,is_retained,nsym2 2724 logical :: found 2725 character(len=500) :: msg 2726 !arrays 2727 integer :: determinant(nsym),inversion(3,3),symrel2(3,3,nsym) 2728 real(dp) :: dtnons(3),tnons2(3,nsym) 2729 2730 ! ********************************************************************* 2731 2732 ABI_WARNING('Removing inversion related symmetrie from initial set') 2733 2734 ! Find the occurence of the inversion symmetry. 2735 call set2unit(inversion) ; inversion=-inversion 2736 2737 is_inv=0; found=.FALSE. 2738 do while (is_inv<nsym .and. .not.found) 2739 is_inv=is_inv+1; found=ALL(symrel(:,:,is_inv)==inversion) 2740 end do 2741 if (found) then 2742 write(msg,'(a,i3)')' The inversion is symmetry operation no. ',is_inv 2743 else 2744 write(msg,'(a)')' The inversion was not found in the symmetries list.' 2745 end if 2746 call wrtout(std_out,msg,'COLL') 2747 2748 ! Find the symmetries that are related through the inversion symmetry 2749 call symdet(determinant,nsym,symrel) 2750 nsym2=0 2751 do is=1,nsym-1 2752 do is2=is+1,nsym 2753 2754 dtnons(:)=tnons(:,is2)-tnons(:,is)-tnons(:,is_inv) 2755 found=ALL(symrel(:,:,is)==-symrel(:,:,is2)).and.isinteger(dtnons,tol8) 2756 2757 if (found) then 2758 nsym2=nsym2+1 2759 ! Retain symmetries with positive determinant 2760 if (ALL(tnons(:,is2)<tol8).and.ALL(tnons(:,is)<tol8)) then 2761 is_retained=is2 ; is_discarded=is 2762 if (determinant(is)==1) then 2763 is_retained=is ; is_discarded=is2 2764 end if 2765 else if (ALL(tnons(:,is2)<tol8)) then 2766 is_retained=is2 ; is_discarded=is 2767 else 2768 is_retained=is ; is_discarded=is2 2769 end if 2770 2771 symrel2(:,:,nsym2)=symrel(:,:,is_retained) 2772 tnons2 (:,nsym2)=tnons (:,is_retained) 2773 write(msg,'(a,i3,a,i3,3a,i3,a)')& 2774 & ' Symmetry operations no. ',is,' and no. ',is2,& 2775 & ' are related through the inversion.',ch10,& 2776 & ' Symmetry operation no. ',is_discarded,' will be suppressed.' 2777 call wrtout(std_out,msg,'COLL') 2778 end if ! found 2779 2780 end do !is2 2781 end do !is 2782 2783 if (nsym2/=(nsym/2).or.nsym==1) then 2784 call wrtout(std_out, ' Program uses the original set of symmetries ', 'COLL') 2785 nsym_out=nsym 2786 ABI_MALLOC(symrel_out,(3,3,nsym)) 2787 ABI_MALLOC(tnons_out,(3,nsym)) 2788 symrel_out(:,:,:)=symrel(:,:,1:nsym) 2789 tnons_out(:,:)=tnons(:,1:nsym) 2790 pinv=1 2791 else 2792 write(msg,'(a)')' Inversion related operations have been suppressed from symmetries list.' 2793 call wrtout(std_out,msg,'COLL') 2794 nsym_out=nsym2 2795 ABI_MALLOC(symrel_out,(3,3,nsym2)) 2796 ABI_MALLOC(tnons_out,(3,nsym2)) 2797 symrel_out(:,:,:)=symrel2(:,:,1:nsym2) 2798 tnons_out(:,:)=tnons(:,1:nsym2) 2799 pinv=-1 2800 end if 2801 2802 end subroutine remove_inversion
m_geometry/rotmat [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
rotmat
FUNCTION
Finds the rotation matrix.
INPUTS
xaxis(3)= vectors defining the x axis zaxis(3)= vectors defining the z axis
OUTPUT
inversion_flag = flag that indicates that an inversion operation on the coordinate system should be done umat(3,3)= matrix that rotates the x=(1 0 0) and z=(0 0 1) to the new values defined in xaxis and zaxis
NOTES
Here I set that the axe x is originally at the 1 0 0 direction and z is originally 0 0 1. So calling rotmat(x',z') will find the rotation matrix for the case in which we rotate the x and z axes from their default values to x' and z'.
SOURCE
1012 subroutine rotmat(xaxis,zaxis,inversion_flag,umat) 1013 1014 !Arguments ------------------------------------ 1015 !scalars 1016 integer,intent(out) :: inversion_flag 1017 !arrays 1018 real(dp),intent(in) :: xaxis(3),zaxis(3) 1019 real(dp),intent(out) :: umat(3,3) 1020 1021 !Local variables------------------------------- 1022 !scalars 1023 real(dp) :: cosine,xmod,zmod 1024 character(len=500) :: msg 1025 !arrays 1026 real(dp) :: yaxis(3) 1027 1028 ! ************************************************************************* 1029 1030 xmod = sqrt(xaxis(1)**2 + xaxis(2)**2 + xaxis(3)**2) 1031 zmod = sqrt(zaxis(1)**2 + zaxis(2)**2 + zaxis(3)**2) 1032 1033 if(xmod < 1.d-8)then 1034 write(msg,'(a,a,a,i6)')& 1035 & 'The module of the xaxis should be greater than 1.d-8,',ch10,& 1036 & 'however, |xaxis|=',xmod 1037 ABI_BUG(msg) 1038 end if 1039 1040 if(zmod < 1.d-8)then 1041 write(msg,'(a,a,a,i6)')& 1042 & 'The module of the zaxis should be greater than 1.d-8,',ch10,& 1043 & 'however, |zaxis|=',zmod 1044 ABI_ERROR(msg) 1045 end if 1046 1047 !verify that both axis are perpendicular 1048 cosine = (xaxis(1)*zaxis(1) + xaxis(2)*zaxis(2) & 1049 & + xaxis(3)*zaxis(3))/(xmod*zmod) 1050 1051 if(abs(cosine) > 1.d-8)then 1052 write(msg,'(a,a,a,i6)')& 1053 & 'xaxis and zaxis should be perpendicular,',ch10,& 1054 & 'however, cosine=',cosine 1055 ABI_BUG(msg) 1056 end if 1057 1058 !new y axis as cross product 1059 yaxis(1) = (zaxis(2)*xaxis(3) - xaxis(2)*zaxis(3))/(xmod*zmod) 1060 yaxis(2) = (zaxis(3)*xaxis(1) - xaxis(3)*zaxis(1))/(xmod*zmod) 1061 yaxis(3) = (zaxis(1)*xaxis(2) - xaxis(1)*zaxis(2))/(xmod*zmod) 1062 1063 !hack to allow inversion operation on coordinate transformation 1064 !uses unlikely large but legal values of proj_x and/or proj_z 1065 !to flag inversion 1066 inversion_flag=0 1067 if(xmod>10._dp .or. zmod>10._dp) then 1068 inversion_flag=1 1069 write(msg, '(4a)' )& 1070 & 'inversion operation will be appended to axis transformation',ch10,& 1071 & 'Action: If you did not intend this, make |z|<10 and |x|<10 ',ch10 1072 call wrtout(std_out,msg) 1073 end if 1074 1075 umat(1,:) = xaxis(:)/xmod 1076 umat(2,:) = yaxis(:) 1077 umat(3,:) = zaxis(:)/zmod 1078 1079 end subroutine rotmat
m_geometry/shellstruct [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
shellstruct
FUNCTION
Calculates shell structure (multiplicities, radii) of an atomic configuration
INPUTS
natom=number of atoms in unit cell xred=reduced coordinates of atoms rprimd=unit cell vectors magv = magnetic ordering of atoms given as 1 and -1, if not given fm is assumed atp = atom on which the perturbation was done
OUTPUT
sdisv(nat)= distance of each shell to central atom (only the first nsh entries are relevant) nsh= number of shells mult(nat) = number of atoms on shell (only the first nsh entries are relevant)
SOURCE
2342 subroutine shellstruct(xred,rprimd,natom,magv,distv,smult,sdisv,nsh,atp,prtvol) 2343 2344 !Arguments ------------------------------------ 2345 !scalars 2346 integer,intent(in) :: natom 2347 integer,intent(in),optional :: atp 2348 integer,intent(in),optional :: prtvol 2349 integer,intent(out) :: nsh 2350 !arrays 2351 real(dp),intent(in) :: rprimd(3,3) 2352 real(dp),intent(in) :: xred(3,natom) 2353 integer,intent(out) :: smult(natom) 2354 integer,intent(in),optional :: magv(natom) 2355 real(dp),intent(out) :: sdisv(natom) 2356 real(dp),intent(out) :: distv(natom) 2357 2358 !Local variables------------------------------- 2359 !scalars 2360 integer :: iatom,atpp,ish,prtvoll 2361 character(len=500) :: msg 2362 real(dp),parameter :: rndfact=10000_dp 2363 !arrays 2364 integer :: iperm(natom),jperm(natom) 2365 real(dp) :: distvh(natom,natom) 2366 real(dp) :: magvv(natom) 2367 2368 ! ************************************************************************* 2369 2370 if (present(magv)) then 2371 magvv=magv 2372 else 2373 magvv=(/ (1, iatom=1,natom) /) 2374 end if 2375 2376 if (present(atp)) then 2377 atpp=atp 2378 else 2379 atpp=1 2380 end if 2381 2382 if (present(prtvol)) then 2383 prtvoll=prtvol 2384 else 2385 prtvoll=1 2386 end if 2387 2388 !DEBUB 2389 write(std_out,*)'shellstruct start' 2390 !END DEBUG 2391 2392 !Calculate ionic distances 2393 call ioniondist(natom,rprimd,xred,distvh,1,magv=int(magvv),atp=atpp) 2394 distv=distvh(1,:) 2395 2396 if (prtvol>2) then 2397 write(std_out,'(a)')' shellstruct ionic distances in cell (distv) : ' 2398 call prmat(distv(1:natom),1,natom,1,std_out) 2399 end if 2400 2401 iperm=(/ (iatom, iatom=1,natom ) /) 2402 jperm=iperm 2403 distv=anint(distv*rndfact)/rndfact 2404 !Sort distances 2405 call sort_dp(natom,distv,iperm,10d-5) 2406 call sort_int(natom,iperm,jperm) 2407 2408 smult=0 2409 sdisv=dot_product(rprimd(1,:),rprimd(1,:))+dot_product(rprimd(2,:),rprimd(2,:))+dot_product(rprimd(3,:),rprimd(3,:)) 2410 2411 nsh=1 2412 smult(1)=1 2413 sdisv(1)=distv(1) 2414 2415 do iatom=2,natom 2416 do ish=1,natom 2417 if (distv(iatom)>sdisv(ish)) then 2418 cycle 2419 else if (distv(iatom)==sdisv(ish)) then 2420 smult(ish)=smult(ish)+1 2421 exit 2422 else if (distv(iatom)<sdisv(ish)) then 2423 smult(ish+1:natom)=smult(ish:natom-1) 2424 sdisv(ish+1:natom)=sdisv(ish:natom-1) 2425 smult(ish)=1 2426 sdisv(ish)=distv(iatom) 2427 nsh=nsh+1 2428 exit 2429 end if 2430 end do 2431 end do 2432 2433 distv=(/ ( distv(jperm(iatom)),iatom=1,natom ) /) 2434 2435 if (prtvoll>2) then 2436 write(msg,'(a,i4,a)')' shellstruct found ',nsh,' shells at distances (sdisv) ' 2437 call wrtout(std_out,msg,'COLL') 2438 call prmat(sdisv(1:nsh),1,nsh,1,std_out) 2439 write(msg,fmt='(a,150i4)')' and multiplicities (smult) ', smult(1:nsh) 2440 call wrtout(std_out,msg,'COLL') 2441 end if 2442 2443 !DEBUB 2444 write(std_out,*)'shellstruct leave' 2445 !END DEBUG 2446 2447 end subroutine shellstruct
m_geometry/spinrot_cmat [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
spinrot_cmat
FUNCTION
Construct 2x2 complex matrix representing the rotation operator in spin-space.
INPUTS
spinrot(4)=components of the spinor rotation matrix computed by getspinrot
OUTPUT
spinrot(2,2)=Rotation matrix (complex array)
SOURCE
948 pure function spinrot_cmat(spinrot) 949 950 !Arguments ------------------------------------ 951 real(dp),intent(in) :: spinrot(4) 952 complex(dpc) :: spinrot_cmat(2,2) 953 954 ! ************************************************************************* 955 956 ! Build rotation matrix from spinrot: 957 ! 958 ! ( cos(phi/2) + i n_z sin(phi/2), (+n_y + i n_x) sin(phi/2) ) 959 ! ( (-n_y + i n_x) sin(phi/2) , cos(phi/2) - i n_z sin(phi/2) ) 960 961 ! spinrot(1)=cos(half*phi) 962 ! spinrot(2)=axis(1)*sin(half*phi) 963 ! spinrot(3)=axis(2)*sin(half*phi) 964 ! spinrot(4)=axis(3)*sin(half*phi) 965 966 ! Rotation in spinor space (same equations as in wfconv) 967 ! TODO: Be careful here as wfconv uses symrel^T to map k-points (listkk) 968 ! thus the inverse of the corresponding symrec. 969 ! This may explain why all the terms with sin(phi/2) change sign (phi --> -phi) 970 971 spinrot_cmat(1,1) = spinrot(1) + j_dpc*spinrot(4) 972 spinrot_cmat(1,2) = spinrot(3) + j_dpc*spinrot(2) 973 spinrot_cmat(2,1) =-spinrot(3) + j_dpc*spinrot(2) 974 spinrot_cmat(2,2) = spinrot(1) - j_dpc*spinrot(4) 975 976 ! My equation 977 !spinrot_cmat(1,1) = spinrot(1) - j_dpc*spinrot(4) 978 !spinrot_cmat(1,2) =-spinrot(3) - j_dpc*spinrot(2) 979 !spinrot_cmat(2,1) = spinrot(3) - j_dpc*spinrot(2) 980 !spinrot_cmat(2,2) = spinrot(1) + j_dpc*spinrot(4) 981 982 end function spinrot_cmat
m_geometry/strainsym [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
strainsym
FUNCTION
For given order of point group, symmetrizes the strain tensor, then produce primitive vectors based on the symmetrized strain.
INPUTS
nsym=order of group. rprimd(3,3)= primitive vectors, to be symmetrized rprimd0(3,3)= reference primitive vectors, already symmetrized symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
OUTPUT
rprimd_symm(3,3)= symmetrized primitive vectors
SOURCE
2982 subroutine strainsym(nsym,rprimd0,rprimd,rprimd_symm,symrel) 2983 2984 use m_linalg_interfaces 2985 2986 !Arguments ------------------------------------ 2987 !scalars 2988 integer,intent(in) :: nsym 2989 !arrays 2990 integer,intent(in) :: symrel(3,3,nsym) 2991 real(dp),intent(in) :: rprimd(3,3),rprimd0(3,3) 2992 real(dp),intent(out) :: rprimd_symm(3,3) 2993 2994 !Local variables------------------------------- 2995 !scalars 2996 integer :: isym 2997 !arrays 2998 integer :: symrel_it(3,3) 2999 real(dp) :: rprimd0_inv(3,3),strain(3,3),strain_symm(3,3),tmp_mat(3,3) 3000 3001 !************************************************************************** 3002 3003 !copy initial rprimd input and construct inverse 3004 rprimd0_inv = rprimd0 3005 call matrginv(rprimd0_inv,3,3) 3006 3007 !define strain as rprimd = strain * rprimd0 (in cartesian frame) 3008 !so strain = rprimd * rprimd0^{-1} 3009 !transform to triclinic frame with rprimd0^{-1} * strain * rprimd0 3010 !giving strain as rprimd0^{-1} * rprimd 3011 call dgemm('N','N',3,3,3,one,rprimd0_inv,3,rprimd,3,zero,strain,3) 3012 3013 !loop over symmetry elements to obtain symmetrized strain matrix 3014 strain_symm = zero 3015 do isym = 1, nsym 3016 3017 ! this loop accumulates symrel^{-1}*strain*symrel into strain_symm 3018 3019 ! mati3inv gives the inverse transpose of symrel 3020 call mati3inv(symrel(:,:,isym),symrel_it) 3021 call dgemm('N','N',3,3,3,one,strain,3,dble(symrel(:,:,isym)),3,zero,tmp_mat,3) 3022 call dgemm('T','N',3,3,3,one,dble(symrel_it),3,tmp_mat,3,one,strain_symm,3) 3023 3024 end do 3025 3026 !normalize by number of symmetry operations 3027 strain_symm = strain_symm/dble(nsym) 3028 3029 !this step is equivalent to r_new = r_old * strain * r_old^{-1} * r_old, 3030 !that is, convert strain back to cartesian frame and then multipy by r_old, 3031 !to get the r_new primitive vectors 3032 3033 call dgemm('N','N',3,3,3,one,rprimd0,3,strain_symm,3,zero,rprimd_symm,3) 3034 3035 end subroutine strainsym
m_geometry/strconv [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
strconv
FUNCTION
If original gprimd is input, convert from symmetric storage mode 3x3 tensor in reduced coordinates "frac" to symmetric storage mode symmetric tensor in cartesian coordinates "cart".
INPUTS
frac(6)=3x3 tensor in symmetric storage mode, reduced coordinates gprimd(3,3)=reciprocal space dimensional primitive translations (bohr^-1)
OUTPUT
cart(6)=symmetric storage mode for symmetric 3x3 tensor in cartesian coords.
NOTES
$cart(i,j)=G(i,a) G(j,b) frac(a,b)$ "Symmetric" storage mode for 3x3 tensor is 6 element array with elements 11, 22, 33, 32, 31, and 21. "cart" may be same array as "frac". If rprimd transpose is input instead of gprimd, then convert tensor in cartesian coordinates to reduced coordinates
SOURCE
3199 subroutine strconv(frac,gprimd,cart) 3200 3201 !Arguments ------------------------------------ 3202 !arrays 3203 real(dp),intent(in) :: frac(6),gprimd(3,3) 3204 real(dp),intent(inout) :: cart(6) ! alias of frac !vz_i 3205 3206 !Local variables------------------------------- 3207 !scalars 3208 integer :: ii,jj 3209 !arrays 3210 real(dp) :: work1(3,3),work2(3,3) 3211 3212 ! ************************************************************************* 3213 3214 work1(1,1)=frac(1) 3215 work1(2,2)=frac(2) 3216 work1(3,3)=frac(3) 3217 work1(3,2)=frac(4) ; work1(2,3)=frac(4) 3218 work1(3,1)=frac(5) ; work1(1,3)=frac(5) 3219 work1(2,1)=frac(6) ; work1(1,2)=frac(6) 3220 3221 ! TODO: these are matmuls, replace or get BLAS 3222 ! work2 = work1 * gprimd^T 3223 do ii=1,3 3224 work2(:,ii)=zero 3225 do jj=1,3 3226 work2(:,ii)=work2(:,ii)+gprimd(ii,jj)*work1(:,jj) 3227 end do 3228 end do 3229 3230 ! work1 = gprimd * work2 = gprimd * input * gprimd^T 3231 do ii=1,3 3232 work1(ii,:)=zero 3233 do jj=1,3 3234 work1(ii,:)=work1(ii,:)+gprimd(ii,jj)*work2(jj,:) 3235 end do 3236 end do 3237 3238 cart(1)=work1(1,1) 3239 cart(2)=work1(2,2) 3240 cart(3)=work1(3,3) 3241 cart(4)=work1(2,3) 3242 cart(5)=work1(1,3) 3243 cart(6)=work1(1,2) 3244 3245 end subroutine strconv
m_geometry/stress_voigt_to_mat [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
stress_voigt_to_mat
FUNCTION
Build 3x3 symmetric stress tensor from stress vector in Voigt notation.
INPUTS
OUTPUT
SOURCE
3155 subroutine stress_voigt_to_mat(stress6, stress_mat) 3156 3157 real(dp),intent(in) :: stress6(6) 3158 real(dp),intent(out) :: stress_mat(3,3) 3159 3160 stress_mat(1,1) = stress6(1) 3161 stress_mat(2,2) = stress6(2) 3162 stress_mat(3,3) = stress6(3) 3163 stress_mat(2,3) = stress6(4) 3164 stress_mat(3,2) = stress6(4) 3165 stress_mat(1,3) = stress6(5) 3166 stress_mat(3,1) = stress6(5) 3167 stress_mat(1,2) = stress6(6) 3168 stress_mat(2,1) = stress6(6) 3169 3170 end subroutine stress_voigt_to_mat
m_geometry/stresssym [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
stresssym
FUNCTION
For given order of point group, symmetrizes the stress tensor, in symmetrized storage mode and cartesian coordinates, using input 3x3 symmetry operators in reduced coordinates. symmetrized tensor replaces input tensor.
INPUTS
gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1) nsym=order of group. sym(3,3,nsym)=symmetry operators (usually symrec=expressed in terms of action on reciprocal lattice primitive translations); integers.
SIDE EFFECTS
stress(6)=stress tensor, in cartesian coordinates, in symmetric storage mode
SOURCE
3059 subroutine stresssym(gprimd,nsym,stress,sym) 3060 3061 !Arguments ------------------------------------ 3062 !scalars 3063 integer,intent(in) :: nsym 3064 !arrays 3065 integer,intent(in) :: sym(3,3,nsym) 3066 real(dp),intent(in) :: gprimd(3,3) 3067 real(dp),intent(inout) :: stress(6) 3068 3069 !Local variables------------------------------- 3070 !scalars 3071 integer :: ii,isym,mu,nu 3072 real(dp) :: summ,tmp 3073 !arrays 3074 real(dp) :: rprimd(3,3),rprimdt(3,3),strfrac(6),tensor(3,3),tt(3,3) 3075 3076 !************************************************************************* 3077 3078 !Obtain matrix of real space dimensional primitive translations 3079 !(inverse tranpose of gprimd), and its transpose 3080 call matr3inv(gprimd,rprimd) 3081 rprimdt=transpose(rprimd) 3082 3083 !Compute stress tensor in reduced coordinates 3084 ! strfrac = rprimd^T * stress * rprimd 3085 call strconv(stress,rprimdt,strfrac) 3086 3087 !Switch to full storage mode 3088 tensor(1,1)=strfrac(1) 3089 tensor(2,2)=strfrac(2) 3090 tensor(3,3)=strfrac(3) 3091 tensor(3,2)=strfrac(4) 3092 tensor(3,1)=strfrac(5) 3093 tensor(2,1)=strfrac(6) 3094 tensor(2,3)=tensor(3,2) 3095 tensor(1,3)=tensor(3,1) 3096 tensor(1,2)=tensor(2,1) 3097 3098 ! these loops are useless - trivial action: 3099 ! tt = tensor / dble(nsym) 3100 ! tensor = zero 3101 do nu=1,3 3102 do mu=1,3 3103 tt(mu,nu)=tensor(mu,nu)/dble(nsym) 3104 tensor(mu,nu)=0.0_dp 3105 end do 3106 end do 3107 3108 !loop over all symmetry operations: 3109 ! tensor = symrec * tt * symrec^T = symrec * rprimd^T * input * rprimd symrec^T 3110 ! TODO: this should be replaced by a little BLAS call or two 3111 do isym=1,nsym 3112 do mu=1,3 3113 do nu=1,3 3114 summ=0._dp 3115 do ii=1,3 3116 tmp=tt(ii,1)*sym(nu,1,isym)+tt(ii,2)*sym(nu,2,isym)+& 3117 & tt(ii,3)*sym(nu,3,isym) 3118 summ=summ+sym(mu,ii,isym)*tmp 3119 end do 3120 tensor(mu,nu)=tensor(mu,nu)+summ 3121 end do 3122 end do 3123 end do 3124 3125 !Switch back to symmetric storage mode 3126 strfrac(1)=tensor(1,1) 3127 strfrac(2)=tensor(2,2) 3128 strfrac(3)=tensor(3,3) 3129 strfrac(4)=tensor(3,2) 3130 strfrac(5)=tensor(3,1) 3131 strfrac(6)=tensor(2,1) 3132 3133 !Convert back stress tensor (symmetrized) in cartesian coordinates 3134 ! stress = gprimd * symrec * rprimd^T * input * rprimd symrec^T * gprimd^T 3135 ! symrec_cart = gprimd * symrec * rprimd^T 3136 ! sym_cart = symrec_cart^-1 ^T = rprimd * sym * gprimd^T 3137 call strconv(strfrac,gprimd,stress) 3138 3139 end subroutine stresssym
m_geometry/symredcart [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
symredcart
FUNCTION
Convert a symmetry operation from reduced coordinates (integers) to cartesian coordinates (reals). Can operate in real or reciprocal space
INPUTS
symred(3,3)=symmetry matrice in reduced coordinates (integers) (real or reciprocal space) aprim(3,3)=real or reciprocal space dimensional primitive translations (see below) bprim(3,3)=real or reciprocal space dimensional primitive translations (see below)
OUTPUT
symcart(3,3)=symmetry matrice in cartesian coordinates (reals)
NOTES
When aprim=rprimd and bprim=gprimd, the routine operates in real space (on a real space symmetry) When aprim=gprimd and bprim=rprimd, the routine operates in reciprocal space (on a real space symmetry)
SOURCE
2920 subroutine symredcart(aprim,bprim,symcart,symred) 2921 2922 !Arguments ------------------------------------ 2923 !arrays 2924 integer,intent(in) :: symred(3,3) 2925 real(dp),intent(in) :: aprim(3,3),bprim(3,3) 2926 real(dp),intent(out) :: symcart(3,3) 2927 2928 !Local variables------------------------------- 2929 !scalars 2930 integer :: ii,jj,kk 2931 real(dp) :: symtmp 2932 !arrays 2933 real(dp) :: work(3,3) 2934 2935 ! ************************************************************************* 2936 2937 work=zero 2938 do kk=1,3 2939 do jj=1,3 2940 symtmp=dble(symred(jj,kk)) 2941 do ii=1,3 2942 work(ii,jj)=work(ii,jj)+bprim(ii,kk)*symtmp 2943 end do 2944 end do 2945 end do 2946 2947 ! work = bprim * symred^T 2948 2949 symcart=zero 2950 do kk=1,3 2951 do jj=1,3 2952 symtmp=work(jj,kk) 2953 do ii=1,3 2954 ! symcart = aprim * work^T = aprim * symred * bprim^T 2955 symcart(ii,jj)=symcart(ii,jj)+aprim(ii,kk)*symtmp 2956 end do 2957 end do 2958 end do 2959 2960 end subroutine symredcart
m_geometry/vdotw_rc_vector [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
vdotw_rc_vector
FUNCTION
Compute the scalar product between two vectors expressed in reduced coordinates First vector is real, the second one is complex. The result is multiplied by (2pi)**2 in case of vectors in reciprocal space to take into account the correct normalisation of the reciprocal lattice vectors
INPUTS
xv(3),xw(3)=Vectors in reduced coordinates met(3,3)=Metric tensor space=Character defining whether we are working in real (r) or reciprocal space (g)
OUTPUT
res=complex scalar product of xv and xw
SOURCE
349 complex(dp) function vdotw_rc_vector(xv, xw, met, space) result(res) 350 351 !Arguments ------------------------------------ 352 character(len=1),intent(in) :: space 353 !arrays 354 real(dp),intent(in) :: met(3,3),xv(3) 355 complex(dpc),intent(in) :: xw(3) 356 357 ! ************************************************************************* 358 359 res = ( met(1,1)* xv(1)*xw(1) & 360 +met(2,2)* xv(2)*xw(2) & 361 +met(3,3)* xv(3)*xw(3) & 362 +met(1,2)*(xv(1)*xw(2) + xv(2)*xw(1)) & 363 +met(1,3)*(xv(1)*xw(3) + xv(3)*xw(1)) & 364 +met(2,3)*(xv(2)*xw(3) + xv(3)*xw(2)) ) 365 366 select case (space) 367 case ('r', 'R') 368 return 369 case ('g', 'G') 370 res= res * (two_pi**2) 371 case default 372 ABI_BUG('Wrong value for space') 373 end select 374 375 end function vdotw_rc_vector
m_geometry/vdotw_rr_vector [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
vdotw_rr_vector
FUNCTION
Compute the scalar product between two vectors expressed in reduced coordinates The result is multiplied by (2pi)**2 in case of vectors in reciprocal space to take into account the correct normalisation of the reciprocal lattice vectors
INPUTS
xv(3),xw(3)=Vectors in reduced coordinates met(3,3)=Metric tensor space=Character defining whether we are working in real (r) or reciprocal space (g)
OUTPUT
res=scalar product of xv and xw
SOURCE
299 real(dp) function vdotw_rr_vector(xv,xw,met,space) result(res) 300 301 !Arguments ------------------------------------ 302 character(len=1),intent(in) :: space 303 !arrays 304 real(dp),intent(in) :: met(3,3),xv(3),xw(3) 305 306 ! ************************************************************************* 307 308 res = ( met(1,1)* xv(1)*xw(1) & 309 +met(2,2)* xv(2)*xw(2) & 310 +met(3,3)* xv(3)*xw(3) & 311 +met(1,2)*(xv(1)*xw(2) + xv(2)*xw(1)) & 312 +met(1,3)*(xv(1)*xw(3) + xv(3)*xw(1)) & 313 +met(2,3)*(xv(2)*xw(3) + xv(3)*xw(2)) ) 314 315 select case (space) 316 case ('r', 'R') 317 return 318 case ('g', 'G') 319 res= res * (two_pi**2) 320 case default 321 ABI_BUG('Wrong value for space') 322 end select 323 324 end function vdotw_rr_vector
m_geometry/wedge_basis [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
wedge_basis
FUNCTION
Calculates the basis vectors a ^ a* for a in rprimd and a* in gprimd, needed for some generalized cross products
INPUTS
rprimd(3,3) : real(dp) matrix gprimd(3,3) : real(dp) matrix normalize,optional : whether to normalize the output vectors
OUTPUT
wedge(3,3,3) : 9 basis vectors of rprimd ^ gprimd
SOURCE
430 subroutine wedge_basis(gprimd,rprimd,wedge,normalize) 431 432 !Arguments --------------------------------------------- 433 ! scalars 434 logical,optional,intent(in) :: normalize 435 !arrays 436 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3) 437 real(dp),intent(out) :: wedge(3,3,3) 438 439 ! local 440 !scalars 441 integer :: igprimd, irprimd 442 real(dp) :: nfac 443 logical :: nvec 444 445 ! ********************************************************************* 446 447 if(present(normalize)) then 448 nvec = normalize 449 else 450 nvec = .FALSE. 451 end if 452 453 do irprimd = 1, 3 454 do igprimd = 1, 3 455 wedge(1,irprimd,igprimd) = rprimd(2,irprimd)*gprimd(3,igprimd) - rprimd(3,irprimd)*gprimd(2,igprimd) 456 wedge(2,irprimd,igprimd) = rprimd(3,irprimd)*gprimd(1,igprimd) - rprimd(1,irprimd)*gprimd(3,igprimd) 457 wedge(3,irprimd,igprimd) = rprimd(1,irprimd)*gprimd(2,igprimd) - rprimd(2,irprimd)*gprimd(1,igprimd) 458 end do 459 end do 460 461 if (nvec) then 462 do irprimd = 1, 3 463 do igprimd = 1, 3 464 if(any(abs(wedge(1:3,irprimd,igprimd)).GT.tol8)) then 465 nfac = SQRT(DOT_PRODUCT(wedge(1:3,irprimd,igprimd),wedge(1:3,irprimd,igprimd))) 466 wedge(1:3,irprimd,igprimd) = wedge(1:3,irprimd,igprimd)/nfac 467 end if 468 end do 469 end do 470 end if 471 472 end subroutine wedge_basis
m_geometry/wedge_product [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
wedge_product
FUNCTION
Calculates the wedge product u^w, given the wedge product basis a^b typically u=(u1 a + u2 b + u3 c) and w = (w1 a* + w2 b* + w3 c*)
INPUTS
u(3) :: real(dp) input vector v(3) :: real(dp) input vector wedgebasis(3,3,3) :: real(dp) input matrix
OUTPUT
produv(3) :: real(dp) output vector
SOURCE
493 subroutine wedge_product(produv,u,v,wedgebasis) 494 495 !Arguments --------------------------------------------- 496 !arrays 497 real(dp),intent(in) :: u(3),v(3),wedgebasis(3,3,3) 498 real(dp),intent(out) :: produv(3) 499 500 ! local 501 !scalars 502 integer :: igprimd, ii, irprimd 503 504 ! ********************************************************************* 505 506 produv(:) = zero 507 do irprimd = 1, 3 508 do igprimd = 1, 3 509 do ii = 1, 3 510 produv(ii) = produv(ii) + u(irprimd)*v(igprimd)*wedgebasis(ii,irprimd,igprimd) 511 end do 512 end do 513 end do 514 515 end subroutine wedge_product
m_geometry/wigner_seitz [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
wigner_seitz
FUNCTION
Calculates a grid of points that falls inside of (and eventually on the surface of) the Wigner-Seitz supercell centered on the origin of the B lattice with primitive translations nmonkh(1)*a_1+nmonkh(2)*a_2+nmonkh(3)*a_3. Subroutine taken from the Wannier90 code. Modified by MG to fulfil abinit coding rules. API slightly changed wrt the wannier90 version.
COPYRIGHT
Copyright (C) 2007 Jonathan Yates, Arash Mostofi, Young-Su Lee, Nicola Marzari, Ivo Souza, David Vanderbilt. This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
center(3)=The Wigner-Seitz cell is centered on this point in reduced coordinates. rmet(3,3)=Real space metric ($\textrm{bohr}^{2}$). kptrlatt(3)=Values defining the supercell. prtvol=If different from 0 print out the points falling inside the W-S cell and the correponding weights. lmax(3)=see Notes below.
OUTPUT
npts=number of points falling inside the Wigner-Seitz cell irvec(3,npts)=Reduced coordinated of the points inside the W-S cell ndegen(npts)=Weigths associated to each point.
SIDE EFFECTS
irvec and ndegen are are allocated with the correct size inside the routine and returned to the caller.
NOTES
The Wannier functions live in a supercell of the real space unit cell. This supercell is mp_grid unit cells long in each direction The algorithm loops over grid points r on a unit cell that is 8 times larger than this primitive supercell. One of these points is in the W-S cell if it is closer to center(:) than any of the other points R where R are the translation vectors of the supercell. In the end npts contains the total number of grid points that have been found in the Wigner-Seitz cell The number of lattice vectors R along each direction of the supercell is defined by lmax.
SOURCE
564 subroutine wigner_seitz(center, lmax, kptrlatt, rmet, npts, irvec, ndegen, prtvol) 565 566 !Arguments ------------------------------------ 567 !scalars 568 integer,optional,intent(in) :: prtvol 569 integer,intent(out) :: npts 570 !arrays 571 integer,intent(in) :: kptrlatt(3,3),lmax(3) 572 integer,allocatable,intent(out) :: irvec(:,:),ndegen(:) 573 real(dp),intent(in) :: center(3),rmet(3,3) 574 575 !Local variables------------------------------- 576 !scalars 577 integer :: in1,in2,in3,l1,l2,l3,ii,icount,n1,n2,n3 578 integer :: l0,l1_max,l2_max,l3_max,nl,verbose,mm1,mm2,mm3 579 real(dp) :: tot,dist_min 580 real(dp),parameter :: TOL_DIST=tol7 581 character(len=500) :: msg 582 !arrays 583 real(dp) :: diff(3) 584 real(dp),allocatable :: dist(:),swap2(:,:),swap1(:) 585 586 ! ************************************************************************* 587 588 verbose = 0; if (PRESENT(prtvol)) verbose = prtvol 589 590 if (kptrlatt(1,2) /= 0 .or. kptrlatt(2,1) /= 0 .or. & 591 kptrlatt(1,3) /= 0 .or. kptrlatt(3,1) /= 0 .or. & 592 kptrlatt(2,3) /= 0 .or. kptrlatt(3,2) /= 0 ) then 593 ABI_ERROR('Off-diagonal elements of kptrlatt must be zero') 594 end if 595 596 n1 = kptrlatt(1,1); n2 = kptrlatt(2,2); n3 = kptrlatt(3,3) 597 l1_max = lmax(1); l2_max = lmax(2); l3_max = lmax(3) 598 599 nl=(2*l1_max+1)*(2*l2_max+1)*(2*l3_max+1) 600 l0=1+l1_max*(1+(2*l2_max+1)**2+(2*l3_max+1)) ! Index of the origin. 601 ABI_MALLOC(dist,(nl)) 602 603 ! Allocate with maximum size 604 mm1 = 2 * n1 + 1 605 mm2 = 2 * n2 + 1 606 mm3 = 2 * n3 + 1 607 ABI_MALLOC(irvec, (3, mm1*mm2*mm3)) 608 ABI_MALLOC(ndegen, (mm1*mm2*mm3)) 609 610 npts = 0 611 do in1=-n1,n1 612 do in2=-n2,n2 613 do in3=-n3,n3 614 615 ! Loop over the nl points R. R=0 corresponds to l1=l2=l3=1, or icount=l0 616 icount = 0 617 do l1=-l1_max,l1_max 618 do l2=-l2_max,l2_max 619 do l3=-l3_max,l3_max 620 ! Calculate |r - R -r0|^2. 621 diff(1) = in1 - l1 * n1 - center(1) 622 diff(2) = in2 - l2 * n2 - center(2) 623 diff(3) = in3 - l3 * n3 - center(3) 624 icount = icount+1 625 dist(icount) = DOT_PRODUCT(diff, MATMUL(rmet, diff)) 626 end do 627 end do 628 end do 629 630 dist_min = MINVAL(dist) 631 632 if (ABS(dist(l0) - dist_min) < TOL_DIST) then 633 npts = npts + 1 634 ndegen (npts) = 0 635 do ii=1,nl 636 if (ABS(dist(ii) - dist_min) < TOL_DIST) ndegen(npts) = ndegen(npts) + 1 637 end do 638 irvec(:, npts) = [in1, in2, in3] 639 end if 640 end do !in3 641 end do !in2 642 end do !in1 643 644 if (verbose >= 1) then 645 write(msg,'(a,i0)')' lattice points in Wigner-Seitz supercell: ',npts 646 call wrtout(std_out,msg,'COLL') 647 do ii=1,npts 648 write(msg,'(a,3(i3),a,i4)')' vector ', irvec(:,ii),' degeneracy: ', ndegen(ii) 649 call wrtout(std_out,msg,'COLL') 650 end do 651 end if 652 653 ! Check the "sum rule" 654 tot = zero 655 do ii=1,npts 656 tot = tot + one/ndegen(ii) 657 end do 658 if (ABS(tot-(n1*n2*n3)) > tol8) then 659 write(msg,'(a,es16.8,a,i0)')'Something wrong in the generation of WS mesh: tot ',tot,' /= ',n1*n2*n3 660 ABI_ERROR(msg) 661 end if 662 663 ABI_FREE(dist) 664 665 ! Reallocate the output with correct size. 666 ABI_MALLOC(swap1,(npts)) 667 swap1(:)=ndegen(1:npts) 668 ABI_FREE(ndegen) 669 ABI_MALLOC(ndegen,(npts)) 670 ndegen=swap1 671 ABI_FREE(swap1) 672 673 ABI_MALLOC(swap2,(3,npts)) 674 swap2(:,:)=irvec(1:3,1:npts) 675 ABI_FREE(irvec) 676 ABI_MALLOC(irvec,(3,npts)) 677 irvec=swap2 678 ABI_FREE(swap2) 679 680 end subroutine wigner_seitz
m_geometry/xcart2xred [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
xcart2xred
FUNCTION
Convert from cartesian coordinates xcart(3,natom) in bohr to
dimensionless reduced coordinates xred(3,natom) by using
xred(mu,ia)=gprimd(1,mu)*xcart(1,ia)
+gprimd(2,mu)*xcart(2,ia)
+gprimd(3,mu)*xcart(3,ia)
where gprimd is the inverse of rprimd
Note that the reverse operation is done by xred2xcart
INPUTS
natom=number of atoms in unit cell rprimd(3,3)=dimensional real space primitive translations (bohr) xcart(3,natom)=cartesian coordinates of atoms (bohr)
OUTPUT
xred(3,natom)=dimensionless reduced coordinates of atoms
SOURCE
1573 subroutine xcart2xred(natom,rprimd,xcart,xred) 1574 1575 !Arguments ------------------------------------ 1576 !scalars 1577 integer,intent(in) :: natom 1578 !arrays 1579 real(dp),intent(in) :: rprimd(3,3),xcart(3,natom) 1580 real(dp),intent(out) :: xred(3,natom) 1581 1582 !Local variables------------------------------- 1583 !scalars 1584 integer :: iatom,mu 1585 !arrays 1586 real(dp) :: gprimd(3,3) 1587 1588 ! ************************************************************************* 1589 1590 call matr3inv(rprimd,gprimd) 1591 do iatom=1,natom 1592 do mu=1,3 1593 xred(mu,iatom)= gprimd(1,mu)*xcart(1,iatom)+gprimd(2,mu)*xcart(2,iatom)+gprimd(3,mu)*xcart(3,iatom) 1594 end do 1595 end do 1596 1597 end subroutine xcart2xred
m_geometry/xred2xcart [ Functions ]
[ Top ] [ m_geometry ] [ Functions ]
NAME
xred2xcart
FUNCTION
Convert from dimensionless reduced coordinates xred(3,natom)
to cartesian coordinates xcart(3,natom) in bohr by using
xcart(mu,ia)=rprimd(mu,1)*xred(1,ia)
+rprimd(mu,2)*xred(2,ia)
+rprimd(mu,3)*xred(3,ia)
Note that the reverse operation is done by xcart2xred.F90
INPUTS
natom=number of atoms in unit cell rprimd(3,3)=dimensional real space primitive translations (bohr) xred(3,natom)=dimensionless reduced coordinates of atoms
OUTPUT
xcart(3,natom)=cartesian coordinates of atoms (bohr)
SOURCE
1622 subroutine xred2xcart(natom,rprimd,xcart,xred) 1623 1624 !Arguments ------------------------------------ 1625 !scalars 1626 integer,intent(in) :: natom 1627 !arrays 1628 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 1629 real(dp),intent(out) :: xcart(3,natom) 1630 1631 !Local variables------------------------------- 1632 !scalars 1633 integer :: iatom,mu 1634 1635 ! ************************************************************************* 1636 1637 do iatom=1,natom 1638 do mu=1,3 1639 xcart(mu,iatom)=rprimd(mu,1)*xred(1,iatom)+rprimd(mu,2)*xred(2,iatom)+rprimd(mu,3)*xred(3,iatom) 1640 end do 1641 end do 1642 1643 end subroutine xred2xcart
m_symtk/reduce2primitive [ Functions ]
[ Top ] [ m_symtk ] [ Functions ]
NAME
reduce2primitive
FUNCTION
Find real space primitive vectors from non-primitive ones and the set of non-integer translations that leave the system invariant
INPUTS
ntranslat=number of translations rprimd(3,3)=dimensional non-primitive vectors in real space (bohr) tolsym=tolerance for the symmetry operations translations(3,ntranslat)=translation vectors, in reduced coordinates
OUTPUT
rprimd_primitive(3,3)=dimensional primitive vectors in real space (bohr)
SOURCE
2824 subroutine reduce2primitive(ntranslat, rprimd, rprimd_primitive, tolsym, translations) 2825 2826 !Arguments ------------------------------------ 2827 !scalars 2828 integer,intent(in) :: ntranslat 2829 real(dp),intent(in) :: tolsym 2830 !arrays 2831 real(dp),intent(in) :: rprimd(3,3),translations(3,ntranslat) 2832 real(dp),intent(out) :: rprimd_primitive(3,3) 2833 2834 2835 !Local variables------------------------------- 2836 !scalars 2837 integer :: idir,itentative,itrans,replace 2838 character(len=500) :: msg 2839 !arrays 2840 real(dp) :: trans_cart(3,ntranslat),trans_red(3,ntranslat) 2841 2842 !************************************************************************** 2843 2844 !These translations should form the primitive lattice when combined with the non-primitive vectors. 2845 !Each translation, in reduced coordinates, should be constituted of rational numbers. 2846 !They should pave the non-primitive cell homogeneously. The issue is to replace 2847 !at least one (or more) of the non-primitive vectors by one (or more) selected translation vectors among the list. 2848 !All translation vectors should be an integer linear combination of the vectors of the new basis. 2849 2850 rprimd_primitive(:,:)=rprimd(:,:) 2851 2852 !First, the reduced coordinates of translation vectors are transferred to the [0,1[ interval 2853 trans_red(:,1:ntranslat)=translations(:,1:ntranslat)-nint(translations(:,1:ntranslat)-tolsym) 2854 2855 !Then, one of the translation vectors with the smallest non-zero first coordinate will replace the first vector, if any. 2856 !Similarly for the three directions. 2857 do idir=1,3 2858 replace=0 2859 do itrans=1,ntranslat 2860 if(trans_red(idir,itrans)>tolsym)then 2861 if(replace==0)then 2862 replace=1 ; itentative=itrans 2863 else 2864 if(trans_red(idir,itentative)>trans_red(idir,itrans)+tolsym)then 2865 itentative=itrans 2866 endif 2867 endif 2868 endif 2869 enddo 2870 if(replace==1)then 2871 ! Change the trans vectors to cartesian coordinates, using "old" rprimd 2872 call xred2xcart(ntranslat,rprimd_primitive,trans_cart,trans_red) 2873 ! Replace rprimd vector with index idir by the selected trans_cart vector 2874 rprimd_primitive(:,idir)=trans_cart(:,itentative) 2875 ! Change the translation vectors to new reduced coordinates using updated rprimd_primitive 2876 call xcart2xred(ntranslat,rprimd_primitive,trans_cart,trans_red) 2877 !Transfer to the [0,1[ interval 2878 trans_red(:,1:ntranslat)=trans_red(:,1:ntranslat)-nint(trans_red(:,1:ntranslat)-tolsym) 2879 endif 2880 enddo ! idir 2881 2882 !Now, check that all translation vectors have zero reduced coordinates. 2883 do itrans=1,ntranslat 2884 do idir=1,3 2885 if (abs(trans_red(idir,itrans))>tolsym) then 2886 write(msg,'(5a)')& 2887 'Did not succeed to find primitive cell from non-primitive one.',ch10,& 2888 'Indeed, there remains a non-vanishing pure translation after reduction.',ch10,& 2889 'Action: this is a bug, contact ABINIT group. Then, use a primitive cell in your input file.' 2890 ABI_ERROR(msg) 2891 end if 2892 enddo 2893 enddo 2894 2895 end subroutine reduce2primitive