C ------------------------------------------------------------------
C PROGRAM Group V1.1
C SPACE GROUP FINDER FOR PERIODIC STRUCTURES, BASED ON THE IDEA
C THAT A GIVEN BRAVAIS LATTICE CAN ONLY SUPPORT A FINITE NUMBER
C OF ROTATION OPERATIONS, WHICH ARE THEN VERIFIED ONE BY ONE.
C FROM THESE OPERATIONS WHICH SHOULD COMMUTE WITH THE SINGLE
C ELECTRON HAMILTONIAN*, ONE CAN COUNT THE ENERGY DEGENERACY OF
C A REGULARLY MESHED FIRST BZ, THUS IN TOTAL ENERGY SUMMATION**
C ONE ONLY NEEDS TO EVALUATE FOR A SMALLER SET OF K-POINTS.
C
C AUG. 4, 1998
C ORIGINAL PROGRAM BY AMES LAB
C MINOR REVISIONS BY JU LI (MIT)
C ------------------------------------------------------------------
PROGRAM Group
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (MMA=10,MML=2,EPS=1D-4)
C **********************************************************
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MML: MAXIMUM NUMBER OF SHELLS TO FULLY ENCLOSE THE CENTER
C warning: IF PARAMETERS MMA AND EPS ARE CHANGED HERE,
C OTHER PLACES MUST CHANGE TOO.
C **********************************************************
DIMENSION H(3,3), G(3,3), HH(3,3), RO(3,3)
DIMENSION TAU(MMA,3),S(MMA,3),SS(MMA,3),SSS(MMA,3)
DIMENSION X(3),XX(3,(2*MML+1)**3,3)
DIMENSION NTYPE(MMA),NT(3)
C MAXIMUM ORDER OF SPACE GROUP = 48
DIMENSION ROT(48,3,3),ROTATION(48,3,3),TRANSLATION(48,3)
CHARACTER *50 BUF, NAME
DATA LP_GROUP /27/
LOGICAL UNITARY, EQUIVALENT
print *, 'Reading instructions...'
print *, ' '
C Read file header: "Name:", NAME, "Lattice:"
READ (*, '(A50,/,A50,/,A50)') BUF,NAME,BUF
C Read in Bravais lattice vectors in certain "characteristic"
C length unit L. Notice that a vector is a row in the matrix,
C and so (x, y, z) = (s1, s2, s3) * A.
READ *, ((H(I,J),J=1,3),I=1,3)
C Form the metric matrix: HH = HH+
DO I=1,3
DO J=1,3
HH(I,J) = 0.D0
DO K=1,3
HH(I,J) = HH(I,J) + H(I,K)*H(J,K)
ENDDO
ENDDO
ENDDO
C Calculate the inverse matrix: G = H^(-1)
CALL MATINV (H,G,VOLUME)
C Now G(:,m) is the reciprocal to lattice vector A(m,:)
READ (*, '(A50)') BUF
C "Number of atoms in the unit cell:"
READ *, NPA
IF (NPA.GT.MMA) THEN
PRINT *, 'Number of atoms in unit cell greater than', MMA
STOP
ENDIF
READ (*, '(A50)') BUF
C "Type ( X Y Z ):"
READ *, (NTYPE(N),(TAU(N,I),I=1,3),N=1,NPA)
C Read in atom type and real space coordinates (in L)
C Transform real space coordinates into coefficients
C of lattice vectors (reduced coordinates)
DO N = 1,NPA
DO I = 1,3
S(N,I) = G(1,I)*TAU(N,1)+G(2,I)*TAU(N,2)+G(3,I)*TAU(N,3)
ENDDO
ENDDO
C ---------------------------------------------------
C Find all possible rotational parts and then verify
C whether they could form space group operations.
C ---------------------------------------------------
C Step 1: select lattice vector transformations which
C preserve length for any of the Bravais lattice vectors:
DO NX = 1,3
C Pick lattice vector NX
DO J = 1,3
IF (J.EQ.NX) THEN
X(J)= 1.D0
ELSE
X(J) = 0.D0
ENDIF
ENDDO
C Calculate its norm
DDV = PRODUCT(X, X, HH)
C Go over surrounding shells
IC = 0
DO 10 N1 = -MML,MML
DO 10 N2 = -MML,MML
DO 10 N3 = -MML,MML
X(1) = N1
X(2) = N2
X(3) = N3
DD = PRODUCT(X, X, HH)
C If norm not equal, look for new ones;
IF (ABS(DD-DDV).GT.EPS) GOTO 10
C else, increase the index
IC = IC+1
C and save this record for vector NX
XX(NX,IC,1) = N1
XX(NX,IC,2) = N2
XX(NX,IC,3) = N3
10 CONTINUE
C total number of records for vector NX
NT(NX) = IC
WRITE (*, '(i2," lattice vectors has the same",
A " length as Bravais vector ",i1)') NT(NX), NX
ENDDO
C Step 2: select unitary transformation matrices
C among NT(1)*NT(2)*NT(3) candidates:
NTRANS = 0
DO 20 IX = 1,NT(1)
DO 20 IY = 1,NT(2)
DO 20 IZ = 1,NT(3)
DO J = 1,3
RO(1,J) = XX(1,IX,J)
RO(2,J) = XX(2,IY,J)
RO(3,J) = XX(3,IZ,J)
ENDDO
IF (.NOT.UNITARY(RO,HH)) GOTO 20
NTRANS = NTRANS + 1
DO I=1,3
DO J=1,3
ROT(NTRANS,I,J) = RO(I,J)
ENDDO
ENDDO
20 CONTINUE
WRITE (*, '(/," Among the", i5, " possible candidates,",
A " we find ", i2, " orthogonal matrices.")')
A NT(1)*NT(2)*NT(3), NTRANS
C Number of group operations:
NGROUP = 0
C Number of non-symorphic operations:
NONSYM = 0
C Loop over rotation
DO 30 NG = 1,NTRANS
DO N = 1,NPA
C Operate on the atoms in the cluster
DO J = 1,3
SS(N,J) = 0.D0
DO K=1,3
SS(N,J) = SS(N,J) + S(N,K)*ROT(NG,K,J)
ENDDO
ENDDO
ENDDO
DO M = 1,NPA
C For all non-primitive translations S(M,:) - SS(1,:)
DO N = 1,NPA
DO K = 1,3
SSS(N,K) = SS(N,K) + S(M,K) - SS(1,K)
ENDDO
ENDDO
C If the two structures are equivalent:
IF (EQUIVALENT(NPA,NTYPE,SSS,NTYPE,S)) GOTO 40
ENDDO
C No, none of the translations are good
GOTO 30
C We have found a valid group operation
40 NGROUP = NGROUP+1
DO K = 1,3
TRANSLATION(NGROUP,K) =
A S(M,K)-SS(1,K)-NINT(S(M,K)-SS(1,K))
DO J = 1,3
ROTATION(NGROUP,K,J) = ROT(NG,K,J)
ENDDO
ENDDO
C If it is also an non-symmorphic operation:
IF (ABS(TRANSLATION(NGROUP,1)).GT.EPS.OR.
A ABS(TRANSLATION(NGROUP,2)).GT.EPS.OR.
A ABS(TRANSLATION(NGROUP,3)).GT.EPS) NONSYM = NONSYM+1
30 CONTINUE
write (*, '(/, " Among them, we found ", i2,
A " symmorphic operations,", /, 18x, "and ", i2,
A " non-symmorphic operations.",/)') NGROUP-NONSYM,NONSYM
print *,'The group operations will be saved in file \"group-op\",'
print *, 'both in terms of lattice vector coefficients (left),'
print *,
A 'and in Cartesian frame and characteristic length L (right).'
print *, ' '
C Write group operations in file "group-op", in both
C reduced and real coordinates:
OPEN (UNIT=LP_GROUP, STATUS='UNKNOWN', FORM='FORMATTED',
A FILE="group-op")
WRITE (LP_GROUP,
A '("Number of group operations found = ",I2)') NGROUP
WRITE (LP_GROUP,
A '("Number of non-symmorphic group operations =",I2)') NONSYM
DO L = 1,NGROUP
DO I = 1,3
X(I) = 0.
DO J = 1,3
X(I) = X(I) + TRANSLATION(L,J)*H(J,I)
RO(I,J) = 0.D0
DO K1=1,3
DO K2=1,3
RO(I,J) = RO(I,J) + G(I,K1)*ROTATION(L,K1,K2)*H(K2,J)
ENDDO
ENDDO
ENDDO
ENDDO
WRITE (LP_GROUP,'(/,"No.",i2)') L
WRITE (LP_GROUP,'(" |",3f9.5," | |",3f9.5," |")')
A ((ROTATION(L,I,J),J=1,3),(RO(I,J),J=1,3),I=1,3)
WRITE (LP_GROUP,'(3x,"+ [",3f9.5," ]",3x,"+ [",3f9.5," ]")')
A (TRANSLATION(L,J),J=1,3), (X(J),J=1,3)
ENDDO
CLOSE (LP_GROUP)
print *,
A 'Generating IBZ k-points (in reciprocal vector coeff.)...'
print *, ' '
C Generate IBZ k-points from a regularly meshed BZ:
READ (*, '(/,A50)') BUF
C "Number of IBZ k-point files:"
READ *, NFILE
DO I = 1,NFILE
READ (*, '(/,A50)') BUF
C "Mesh density + shift:"
READ *, NX,NY,NZ,SX,SY,SZ
READ (*, '(A50)') BUF
C "filename:"
READ (*, '(A50)') BUF
NRK = IBZ_K_POINTS(G,NX,NY,NZ,SX,SY,SZ,NGROUP,ROTATION,BUF)
PRINT *, NRK, ' k-points saved in ', BUF(1:INDEX(BUF,' ')-1)
ENDDO
STOP
END
C GENERATE IBZ K-POINT GRID FROM A REGULARLY MESHED
C BZ: NX, NY, NZ ARE THE NUMBER OF K-POINTS IN THREE
C DIRECTIONS. SX, SY, SZ SHIFT THE GRID FROM THE ORIGIN.
INTEGER FUNCTION IBZ_K_POINTS
A (G,NX,NY,NZ,SX,SY,SZ,NGROUP,ROTATION,FILENAME)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (MAXK=40000,EPS=1D-4)
DIMENSION G(3,3),ROTATION(48,3,3)
DIMENSION XK(3,MAXK),RKTRAN(3),KDEGEN(MAXK)
CHARACTER *50 FILENAME
DATA LP_KPTS /28/
IF (NX*NY*NZ.GT.MAXK) THEN
PRINT *, 'IBZ_K_POINTS(): increase NMAX'
STOP
ENDIF
C A uniform grid of k points:
NK = 0
DO I = 1,NX
DO J = 1,NY
DO K = 1,NZ
NK = NK+1
XK(1,NK) = (I-1+SX)/NX
XK(2,NK) = (J-1+SY)/NY
XK(3,NK) = (K-1+SZ)/NZ
KDEGEN(NK) = 1
ENDDO
ENDDO
ENDDO
NRK = 0
DO 10 I = 1,NK
C Only RK will be recorded:
IF (KDEGEN(I).EQ.0) GOTO 10
NRK = NRK+1
DO 20 L = 1,NGROUP
C Operate on RK with all symmetry operations
DO J = 1,3
RKTRAN(J) = ROTATION(L,J,1)*XK(1,I)
A + ROTATION(L,J,2)*XK(2,I)
A + ROTATION(L,J,3)*XK(3,I)
C Translate to first BZ
RKTRAN(J) = RKTRAN(J) - NINT(RKTRAN(J))
C Make it [0,1) for comparison:
IF (RKTRAN(J).LT.0) RKTRAN(J) = RKTRAN(J)+1.D0
ENDDO
C Check for upward k-points equivalence
DO 30 K = I+1,NK
IF (KDEGEN(K).EQ.0) GOTO 30
C Both the rotated k-vector and its inverse
C (time-reversal symmetry if no spin-polarization)
IF ( (ABS(RKTRAN(1)-XK(1,K)).LT.EPS.AND.
A ABS(RKTRAN(2)-XK(2,K)).LT.EPS.AND.
A ABS(RKTRAN(3)-XK(3,K)).LT.EPS) .OR.
A (ABS(1.D0-RKTRAN(1)-XK(1,K)).LT.EPS.AND.
A ABS(1.D0-RKTRAN(2)-XK(2,K)).LT.EPS.AND.
A ABS(1.D0-RKTRAN(3)-XK(3,K)).LT.EPS) ) THEN
KDEGEN(I) = KDEGEN(I)+1
KDEGEN(K) = 0
ENDIF
30 CONTINUE
20 CONTINUE
10 CONTINUE
C THE K-POINTS WILL BE STORED IN RECIPROCAL LATTICE
C VECTOR COEFFICIENTS [0,1)
OPEN (UNIT=LP_KPTS, STATUS='UNKNOWN', FORM='FORMATTED',
A FILE=FILENAME(1:INDEX(FILENAME,' ')-1))
WRITE (LP_KPTS, *) NRK
DO N = 1,NK
IF (KDEGEN(N).GT.0) THEN
C RKTRAN(1) = 2.D0*(XK(1,N)*G(1,1)+XK(2,N)*G(1,2)+XK(3,N)*G(1,3))
C RKTRAN(2) = 2.D0*(XK(1,N)*G(2,1)+XK(2,N)*G(2,2)+XK(3,N)*G(2,3))
C RKTRAN(3) = 2.D0*(XK(1,N)*G(3,1)+XK(2,N)*G(3,2)+XK(3,N)*G(3,3))
WRITE (LP_KPTS,'(4f17.14)')
A XK(1,N),XK(2,N),XK(3,N),dble(KDEGEN(N))/NK
C A RKTRAN(1),RKTRAN(2),RKTRAN(3),dble(KDEGEN(N))/NK
ENDIF
ENDDO
CLOSE(LP_KPTS)
IBZ_K_POINTS = NRK
RETURN
END
C SUBROUTINE TO EVALUATE THE INNER PRODUCT OF
C TWO VECTORS UNDER ARBITRARY METRIC.
DOUBLE PRECISION FUNCTION PRODUCT (X, Y, HH)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION X(3), HH(3,3), Y(3)
PRODUCT = 0.D0
DO I = 1,3
DO J = 1,3
PRODUCT = PRODUCT + X(I) * HH(I,J) * Y(J)
ENDDO
ENDDO
RETURN
END
C CHECK IF THE TRANSFORMATION R, EXPRESSED FROM AND INTO
C LATTICE VECTOR COEFFICIENTS, IS UNITARY: RHH+R+ = HH+.
LOGICAL FUNCTION UNITARY (R, HH)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (EPS=1D-4)
DIMENSION R(3,3),HH(3,3),AR(3,3)
DO I=1,3
DO J=1,3
AR(I,J) = 0.D0
DO K1=1,3
DO K2=1,3
AR(I,J) = AR(I,J) + R(I,K1) * HH(K1,K2) * R(J,K2)
ENDDO
ENDDO
ENDDO
ENDDO
DO I=1,3
DO J=1,3
IF (ABS(HH(I,J)-AR(I,J)).GT.EPS) THEN
UNITARY = .FALSE.
RETURN
ENDIF
ENDDO
ENDDO
UNITARY = .TRUE.
RETURN
END
C Compare two periodic group of atoms to see if they
C are structurally equivalent, indices notwithstanding.
LOGICAL FUNCTION EQUIVALENT (NPA,NTYPE1,S1,NTYPE2,S2)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (MMA=10,EPS=1D-4)
DIMENSION S1(MMA,3),S2(MMA,3)
DIMENSION NTYPE1(MMA),NTYPE2(MMA)
C Assume that no atoms in each group are close within EPS
DO 10 M = 1,NPA
DO N = 1,NPA
IF (NTYPE1(M).EQ.NTYPE2(N).AND.
A ABS(S1(M,1)-S2(N,1)-NINT(S1(M,1)-S2(N,1))).LT.EPS.AND.
A ABS(S1(M,2)-S2(N,2)-NINT(S1(M,2)-S2(N,2))).LT.EPS.AND.
A ABS(S1(M,3)-S2(N,3)-NINT(S1(M,3)-S2(N,3))).LT.EPS) GOTO 10
ENDDO
EQUIVALENT = .FALSE.
RETURN
10 CONTINUE
EQUIVALENT = .TRUE.
RETURN
END
SUBROUTINE MATINV(A,B,C)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION A(3,3),B(3,3)
D11=A(2,2)*A(3,3)-A(2,3)*A(3,2)
D22=A(3,3)*A(1,1)-A(3,1)*A(1,3)
D33=A(1,1)*A(2,2)-A(1,2)*A(2,1)
D12=A(2,3)*A(3,1)-A(2,1)*A(3,3)
D23=A(3,1)*A(1,2)-A(3,2)*A(1,1)
D31=A(1,2)*A(2,3)-A(1,3)*A(2,2)
D13=A(2,1)*A(3,2)-A(3,1)*A(2,2)
D21=A(3,2)*A(1,3)-A(1,2)*A(3,3)
D32=A(1,3)*A(2,1)-A(2,3)*A(1,1)
C=A(1,1)*D11+A(1,2)*D12+A(1,3)*D13
B(1,1)=D11/C
B(2,2)=D22/C
B(3,3)=D33/C
B(1,2)=D21/C
B(2,3)=D32/C
B(3,1)=D13/C
B(2,1)=D12/C
B(3,2)=D23/C
B(1,3)=D31/C
RETURN
END
C ------------------------------------------------------------------
C Convention: row vectors, and {R|t}r = rR + t
C * All rotational and translational operations commute with
C the total Hamiltonian H(r1,r2,..R1,R2..) where {R1,R2..} are
C nuclear coordinates. However if O{R1,R2..} = {R1,R2..}, then
C OH=HO even if O operates on {r1,r2,..} only. Furthermore, one
C assume the ground state charge density rho(r) to be the identity
C representation of the symmetry group, which is self-consistent
C but in principle not necessary.
C ** Although IBZ is enough for total energy summation, it can
C lead to fictitious results in LDOS calculation if the symmetry
C operations contains non-primitive translations, such as for the
C two Si atoms in 2H SiC, though physically equivalent, may have
C different IBZ summed charges. The reason is simply although
C O\psi(k)=\psi(Ok) has the same energy as \psi(k), the charge
C density on atom i in \psi(k) can shift to an equivalent but
C different site after O\psi(k): so although \sum_{BZ}E(k) =
C N\sum_{IBZ}E(k) and \sum_{BZ}\rho(k,r) has the full symmetry of
C the crystal, \sum_{IBZ}\rho(k,r) may not. This problem has a
C simple solution: one shows that charges belonging to one
C class of atoms (two Si above) never goes other classes for
C any k, so one just needs to symmetrize within each class.
C ------------------------------------------------------------------