findsymlat.f90

Go to the documentation of this file.

! Copyright (C) 2002-2007 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU Lesser General Public
! License. See the file COPYING for license details.

!BOP
! !ROUTINE: findsymlat
! !INTERFACE:
subroutine findsymlat
! !USES:
use modmain
use modtddft
! !DESCRIPTION:
!   Finds the point group symmetries which leave the Bravais lattice invariant.
!   Let $A$ be the matrix consisting of the lattice vectors in columns, then
!   $$ g=A^{\rm T}A $$
!   is the metric tensor. Any $3\times 3$ matrix $S$ with elements $-1$, 0 or 1
!   is a point group symmetry of the lattice if $\det(S)$ is $-1$ or 1, and
!   $$ S^{\rm T}gS=g. $$
!   The first matrix in the set returned is the identity.
!
! !REVISION HISTORY:
!   Created January 2003 (JKD)
!   Removed arguments and simplified, April 2007 (JKD)
!EOP
!BOC
implicit none
! local variables
integer md,sym(3,3),its,i,j
integer i11,i12,i13,i21,i22,i23,i31,i32,i33
real(8) s(3,3),g(3,3),sgs(3,3),sc(3,3)
real(8) c(3,3),v(3),t1
! determine metric tensor
call r3mtm(avec,avec,g)
! loop over all possible symmetry matrices
nsymlat=0
do i11=-1,1; do i12=-1,1; do i13=-1,1
do i21=-1,1; do i22=-1,1; do i23=-1,1
do i31=-1,1; do i32=-1,1; do i33=-1,1
  sym(1,1)=i11; sym(1,2)=i12; sym(1,3)=i13
  sym(2,1)=i21; sym(2,2)=i22; sym(2,3)=i23
  sym(3,1)=i31; sym(3,2)=i32; sym(3,3)=i33
! determinant of matrix
  md=i3mdet(sym)
! matrix should be orthogonal
  if (abs(md) /= 1) goto 10
! check invariance of metric tensor
  s(1:3,1:3)=dble(sym(1:3,1:3))
  call r3mtm(s,g,c)
  call r3mm(c,s,sgs)
  do j=1,3
    do i=1,3
      if (abs(sgs(i,j)-g(i,j)) > epslat) goto 10
    end do
  end do
! check invariance of spin-spiral q-vector if required
  if (spinsprl) then
    call r3mtv(s,vqlss,v)
    t1=abs(vqlss(1)-v(1))+abs(vqlss(2)-v(2))+abs(vqlss(3)-v(3))
    if (t1 > epslat) goto 10
  end if
! check invariance of electric field if required
  if (tefield) then
    call r3mv(s,efieldl,v)
    t1=abs(efieldl(1)-v(1))+abs(efieldl(2)-v(2))+abs(efieldl(3)-v(3))
    if (t1 > epslat) goto 10
  end if
! check invariance of static A-field if required
  if (tafield) then
    call r3mv(s,afieldl,v)
    t1=abs(afieldl(1)-v(1))+abs(afieldl(2)-v(2))+abs(afieldl(3)-v(3))
    if (t1 > epslat) goto 10
  end if
! check invariance of time-dependent A-field if required
  if (tafieldt) then
! symmetry matrix in Cartesian coordinates
    call r3mm(s,ainv,c)
    call r3mm(avec,c,sc)
    do its=1,ntimes
      call r3mv(sc,afieldt(:,its),v)
      t1=abs(afieldt(1,its)-v(1)) &
        +abs(afieldt(2,its)-v(2)) &
        +abs(afieldt(3,its)-v(3))
      if (t1 > epslat) goto 10
    end do
  end if
  nsymlat=nsymlat+1
  if (nsymlat > 48) then
    write(*,*)
    write(*,'("Error(findsymlat): more than 48 symmetries found")')
    write(*,'(" (lattice vectors may be linearly dependent)")')
    write(*,*)
    stop
  end if
! store the symmetry and determinant in global arrays
  symlat(1:3,1:3,nsymlat)=sym(1:3,1:3)
  symlatd(nsymlat)=md
10 continue
end do; end do; end do
end do; end do; end do
end do; end do; end do
if (nsymlat == 0) then
  write(*,*)
  write(*,'("Error(findsymlat): no symmetries found")')
  write(*,*)
  stop
end if
! make the first symmetry the identity
do i=1,nsymlat
  if ((symlat(1,1,i) == 1).and.(symlat(1,2,i) == 0).and.(symlat(1,3,i) == 0) &
 .and.(symlat(2,1,i) == 0).and.(symlat(2,2,i) == 1).and.(symlat(2,3,i) == 0) &
 .and.(symlat(3,1,i) == 0).and.(symlat(3,2,i) == 0).and.(symlat(3,3,i) == 1)) &
  then
    sym(1:3,1:3)=symlat(1:3,1:3,1)
    symlat(1:3,1:3,1)=symlat(1:3,1:3,i)
    symlat(1:3,1:3,i)=sym(1:3,1:3)
    md=symlatd(1)
    symlatd(1)=symlatd(i)
    symlatd(i)=md
    exit
  end if
end do
! index to the inverse of each operation
do i=1,nsymlat
  call i3minv(symlat(:,:,i),sym)
  do j=1,nsymlat
    if (all(symlat(1:3,1:3,j) == sym(1:3,1:3))) then
      isymlat(i)=j
      goto 30
    end if
  end do
  write(*,*)
  write(*,'("Error(findsymlat): inverse operation not found")')
  write(*,'(" for lattice symmetry ",I0)') i
  write(*,*)
  stop
30 continue
end do
! determine the lattice symmetries in Cartesian coordinates
do i=1,nsymlat
  s(1:3,1:3)=dble(symlat(1:3,1:3,i))
  call r3mm(s,ainv,c)
  call r3mm(avec,c,symlatc(:,:,i))
end do

contains

pure integer function i3mdet(a)
implicit none
! arguments
integer, intent(in) :: a(3,3)
! determinant of integer matrix
i3mdet=a(1,1)*(a(2,2)*a(3,3)-a(3,2)*a(2,3)) &
      +a(2,1)*(a(3,2)*a(1,3)-a(1,2)*a(3,3)) &
      +a(3,1)*(a(1,2)*a(2,3)-a(2,2)*a(1,3))
end function

pure subroutine i3minv(a,b)
implicit none
! arguments
integer, intent(in) :: a(3,3)
integer, intent(out) :: b(3,3)
! local variables
integer md
md=i3mdet(a)
b(1,1)=md*(a(2,2)*a(3,3)-a(2,3)*a(3,2))
b(2,1)=md*(a(2,3)*a(3,1)-a(2,1)*a(3,3))
b(3,1)=md*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
b(1,2)=md*(a(1,3)*a(3,2)-a(1,2)*a(3,3))
b(2,2)=md*(a(1,1)*a(3,3)-a(1,3)*a(3,1))
b(3,2)=md*(a(1,2)*a(3,1)-a(1,1)*a(3,2))
b(1,3)=md*(a(1,2)*a(2,3)-a(1,3)*a(2,2))
b(2,3)=md*(a(1,3)*a(2,1)-a(1,1)*a(2,3))
b(3,3)=md*(a(1,1)*a(2,2)-a(1,2)*a(2,1))
end subroutine

end subroutine
!EOC