Skip to content
GitLab
Commit c104ae6b authored by Gustav Bihlmayer's avatar Gustav Bihlmayer
Browse files

Changes in rw_symfile.f for some special film cases, modification in 

first_glance.f to read in l_gamma (w/o hybrid).
Updated atom_sym.f, crystal.f, inpgen.f90, lattice2.f, setab.f, and
struct_input.f for "Examples sorted by Bravais lattice" on web-page.
parent 50644063
Hide whitespace changes
Inline Side-by-side
init/first_glance.f
View file @ c104ae6b
......@@ -136,21 +136,26 @@ c
IF (.not.l_kpts) THEN
WRITE (6,*) 'No kpts-file exists, trying to generate it'
DO line = 1,6
READ (5,*,END=98,ERR=98)
READ (5,*,END=95,ERR=95)
ENDDO
READ (5,'(5x,i5,3(4x,i2),7x,l1)',END=97,ERR=97)
+ nkpt,nmop(1),nmop(2),nmop(3),l_gamma
GOTO 96
97 BACKSPACE (5)
READ (5,'(5x,i5,3(4x,i2))',END=98,ERR=98) nkpt
nmop(1) = 0 ; nmop(2) = 0 ; nmop(3) = 0
READ (5,'(5x,i5,3(4x,i2))',END=98,ERR=98)
+ nkpt,nmop(1),nmop(2),nmop(3)
l_gamma=.false.
GOTO 96
98 BACKSPACE (5)
READ (5,'(5x,i5)',END=95,ERR=95) nkpt
nmop(1) = 0 ; nmop(2) = 0 ; nmop(3) = 0 ; l_gamma=.false.
GOTO 96
98 WRITE (6,*) 'Since you did not provide a kpts-file, you'
95 WRITE (6,*) 'Since you did not provide a kpts-file, you'
WRITE (6,*) 'should give k-mesh information at the end of'
WRITE (6,*) 'the inp-file, at least how many k-points, e.g.'
WRITE (6,*) 'nkpt= 100 -- or give the divisions n in xyz:'
WRITE (6,*) 'nkpt= 36,nx= 6,ny= 6,nz= 8'
WRITE (6,*) 'nkpt= 36,nx= 6,ny= 6,nz= 8,gamma=F'
CLOSE (6)
STOP
ENDIF
......
inpgen/atom_sym.f
View file @ c104ae6b
......@@ -178,10 +178,8 @@
td = matmul( bs, tc )
tc = matmul( td, as )
mmrot(:,:,ng) = NINT(tc)
write(*,*) i_c, ttr(:,ng)
ttr2(:,1) = matmul( lmat(:,:,i_c), ttr(:,ng) )
ttr(:,ng) = ttr2(:,1)
write(*,*) mmrot(:,:,ng),ttr(:,ng)
ENDDO
ENDIF
mmrot2(:,:,1:ngen) = mmrot(:,:,2:ngen+1)
......@@ -262,7 +260,6 @@
mrot(:,:,n) = mmrot(:,:,n)
tau(:,n) = ttr(:,n)
index_op(n) = n
write(*,*) n,mrot(:,:,n),tau(:,n)
ENDDO
!---> check that the group is closed, etc.
......
inpgen/crystal.f
View file @ c104ae6b
......@@ -8,7 +8,7 @@
> dbgfh,errfh,outfh,dispfh,dispfn,
> cal_symm,cartesian,symor,oldfleur,
> natin,natmax,nop48,
> atomid,atompos,a1,a2,a3,aa,scale,noangles,i_c,
> atomid,atompos,a1,a2,a3,aa,scale,i_c,
< invs,zrfs,invs2,nop,nop2,
< ngen,mmrot,ttr,ntype,nat,nops,
< neq,ntyrep,zatom,natype,natrep,natmap,
......@@ -25,7 +25,7 @@
IMPLICIT NONE
!===> Arguments
LOGICAL, INTENT(IN) :: cal_symm,cartesian,oldfleur,noangles
LOGICAL, INTENT(IN) :: cal_symm,cartesian,oldfleur
INTEGER, INTENT(IN) :: ngen,natmax,nop48,i_c
INTEGER, INTENT(IN) :: dbgfh,errfh,outfh,dispfh ! file handles, mainly 6
REAL, INTENT(IN) :: aa
......@@ -99,7 +99,7 @@
!---> generate lattice matrices (everything in mod_lattice)
!
CALL setab(
> a1,a2,a3,aa,scale,noangles,
> a1,a2,a3,aa,scale,
< amat,bmat,aamat,bbmat,amatinv,omtil)
......
inpgen/inpgen.f90
View file @ c104ae6b
......@@ -22,7 +22,7 @@ PROGRAM inpgen
INTEGER natmax,nop48,nline,natin,ngen,i,j
INTEGER nops,no3,no2,na,numSpecies,i_c
INTEGER infh,errfh,bfh,warnfh,symfh,dbgfh,outfh,dispfh
LOGICAL cal_symm,checkinp,newSpecies,noangles
LOGICAL cal_symm,checkinp,newSpecies
LOGICAL cartesian,oldfleur,l_hyb ,inistop
REAL aa
......@@ -73,7 +73,7 @@ PROGRAM inpgen
& natmax,nop48,&
& nline,xl_buffer,buffer,&
& title,input%film,cal_symm,checkinp,sym%symor,&
& cartesian,oldfleur,a1,a2,a3,vacuum%dvac,aa,scale,noangles,i_c,&
& cartesian,oldfleur,a1,a2,a3,vacuum%dvac,aa,scale,i_c,&
& factor,natin,atomid,atompos,ngen,mmrot,ttr,&
& l_hyb,noco%l_soc,noco%l_ss,noco%theta,noco%phi,noco%qss,inistop)!keep
......@@ -113,7 +113,7 @@ PROGRAM inpgen
& dbgfh,errfh,outfh,dispfh,dispfn,&
& cal_symm,cartesian,sym%symor,input%film,&
& natin,natmax,nop48,&
& atomid,atompos,a1,a2,a3,aa,scale,noangles,i_c,&
& atomid,atompos,a1,a2,a3,aa,scale,i_c,&
& sym%invs,sym%zrfs,sym%invs2,sym%nop,sym%nop2,&
& ngen,mmrot,ttr,atoms%ntype,atoms%nat,nops,&
& atoms%neq,ntyrep,atoms%zatom,natype,natrep,natmap,&
......
inpgen/lattice2.f
View file @ c104ae6b
......@@ -8,7 +8,7 @@
CONTAINS
SUBROUTINE lattice2(
> buffer,xl_buffer,errfh,bfh,nline,
< a1,a2,a3,aa,scale,noangles,i_c,ios )
< a1,a2,a3,aa,scale,mat,i_c,ios )
USE m_constants
IMPLICIT NONE
......@@ -18,23 +18,22 @@
INTEGER, INTENT (IN) :: xl_buffer
CHARACTER(len=xl_buffer), INTENT(IN) :: buffer
REAL, INTENT (OUT) :: a1(3),a2(3),a3(3)
REAL, INTENT (OUT) :: aa
REAL, INTENT (OUT) :: scale(3)
REAL, INTENT (OUT) :: aa ! overall scaling constant
REAL, INTENT (OUT) :: scale(3),mat(3,3) ! for trigonal lattices
INTEGER, INTENT (OUT) :: i_c,ios
LOGICAL, INTENT (OUT) :: noangles
!==> Local Variables
CHARACTER(len=40) :: latsys
REAL :: a0
REAL :: a0,a_rho
REAL :: a,b,c
REAL :: alpha,beta,gamma
REAL :: c1(3),c2(3),c3(3)
REAL :: ar,br,cr,b1,b2,am(3,3)
REAL :: ca,cb,at
INTEGER :: i,j,err,i1,i2
LOGICAL :: noangles
REAL, PARAMETER :: eps = 1.0e-7
REAL, PARAMETER :: h = 0.5, o = 0.0 ,
REAL, PARAMETER :: eps = 1.0e-7,
& sqrt2 = 1.4142135623730950,
& sqrt3 = 1.7320508075688773,
& sqrt32 = 0.86602540378444,
......@@ -50,9 +49,9 @@
+ -0.5, 0.5, 0.5, ! 3: body-centered : I
& 0.5, -0.5, 0.5,
& 0.5, 0.5, -0.5,
+ h, msqrt32, o, ! 4: hexagonal-P : hP, hcp
& h, sqrt32, o,
& o, o, 1.0,
+ 0.5, msqrt32, 0.0, ! 4: hexagonal-P : hP, hcp
& 0.5, sqrt32, 0.0,
& 0.0, 0.0, 1.0,
+ 0.0, -1.0, 1.0, ! 5: hexagonal-R : hR, trigonal
& sqrt32, 0.5, 1.0,
& msqrt32, 0.5, 1.0,
......@@ -66,9 +65,6 @@
& 0.0, 0.5, 0.5,
& 0.0, -0.5, 0.5/
!===> 12: monoclinic-P (mP)
!===> 13: monoclinic-P (mS) (mA) (mB) (mC)
!===> namelists
NAMELIST /lattice/ latsys,a0,a,b,c,alpha,beta,gamma
......@@ -76,7 +72,7 @@
latsys = ' ' ; a0 = 0.0
a = 0.0 ; b = 0.0 ; c = 0.0
alpha = 0.0 ; beta = 0.0 ; gamma = 0.0
scale = 0.0
scale = 0.0 ; mat = 0.0
READ (bfh,lattice,err=911,end=911,iostat=ios)
......@@ -179,6 +175,7 @@
a3 = lmat(:,3,i_c)
IF ( a.NE.b ) err = 51
IF ( b.NE.c ) err = 51
IF ( alpha.EQ.0.0 .OR.
& alpha.NE.beta .OR. alpha.NE.gamma ) err = 52
......@@ -187,7 +184,8 @@
am(1,:) = at * am(1,:)
am(2,:) = at * am(2,:)
am(3,:) = cos(asin(at)) * am(3,:)
a1 = am(:,1) ; a2 = am(:,2) ; a3 = am(:,3)
a1 = am(:,1)*a ; a2 = am(:,2)*b ; a3 = am(:,3)*c
scale = 1.0
CALL angles( am )
......@@ -199,26 +197,35 @@
noangles=.false.
i_c = 5
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
IF ( a.NE.b .OR. a.NE.c ) err = 53
IF ( alpha.NE.beta .OR. alpha.NE.gamma ) err = 54
IF ( a.NE.b ) err = 53
IF ( alpha.NE.90 ) err = 54
IF ( alpha.NE.beta .OR. gamma.NE.120 ) err = 54
b1=(1.0 -sqrt(cos(ar)-cos(2.0*ar)))/cos(ar)
b2= 1.0 /sqrt( 2.0 + b1*b1 )
am = b2
DO i=1,3
am(i,i) = b1*b2
ENDDO
a1 = am(:,1)
a2 = am(:,2)
a3 = am(:,3)
mat(:,1) = a*lmat(:,1,4) ! to transfer atom
mat(:,2) = a*lmat(:,2,4) ! positions hex -> trig
mat(:,3) = c*lmat(:,3,4) ! in struct_inp.f
! transform hex -> rho
a_rho = sqrt( 3*a*a + c*c)/3.0
ar = acos( 1. - 9./(6.+2*(c/a)**2) )
at = sqrt( 2.0 / 3.0 * ( 1 - cos(ar) ) )
a = a_rho * at
b = a_rho * at
c = a_rho * cos( asin(at) )
ca = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3)
ca = ca/sqrt(a1(1)*a1(1)+a1(2)*a1(2)+a1(3)*a1(3))
ca = ca/sqrt(a2(1)*a2(1)+a2(2)*a2(2)+a2(3)*a2(3))
ca = acos(ca)*180/pi_const
am(:,1) = a1 ; am(:,2) = a2 ; am(:,3) = a3
am(1,:) = a*am(1,:)
am(2,:) = b*am(2,:)
am(3,:) = c*am(3,:)
a1 = am(:,1) ; a2 = am(:,2) ; a3 = am(:,3)
scale = 1.0
WRITE (6,*) 'dbg:trigonal angle =',ca
CALL angles( am )
!===> 6: tetragonal-P (tP) st
......@@ -335,9 +342,11 @@
+ ,calledby ="lattice2")
ENDIF
CALL brvmat ( alpha, beta, gamma, am )
a1 = am(:,1)
a2 = am(:,2)
a3 = am(:,3)
a1 = a*am(:,1)
a2 = b*am(:,2)
a3 = c*am(:,3)
scale = 1.0
CALL angles( am )
!===> 13: monoclinic-C (mC)
......@@ -432,10 +441,11 @@
noangles=.false.
i_c = 1
CALL brvmat ( alpha, beta, gamma, am )
a1 = am(:,1)
a2 = am(:,2)
a3 = am(:,3)
a1 = a*am(:,1)
a2 = b*am(:,2)
a3 = c*am(:,3)
CALL angles( am )
scale = 1.0
ELSE
WRITE (errfh,*)
......
inpgen/setab.f
View file @ c104ae6b
......@@ -5,7 +5,7 @@
!*********************************************************************
CONTAINS
SUBROUTINE setab(
> a1,a2,a3,aa,scale,noangles,
> a1,a2,a3,aa,scale,
< amat,bmat,aamat,bbmat,amatinv,omtil)
USE m_constants
......@@ -15,7 +15,6 @@
!==> Arguments
REAL, INTENT (IN) :: aa
REAL, INTENT (IN) :: a1(3),a2(3),a3(3),scale(3)
LOGICAL, INTENT (IN) :: noangles
REAL, INTENT (OUT) :: amat(3,3),bmat(3,3),amatinv(3,3)
REAL, INTENT (OUT) :: aamat(3,3),bbmat(3,3)
REAL, INTENT (OUT) :: omtil
......@@ -47,8 +46,6 @@
! matrices of lattice vectors in full Cartesian units
IF (noangles) THEN ! we scale the coordinate system in the case of cP, cI, cF, tI etc.
DO i=1,3
amat(i,1) = aa*scale(i)*a1(i)
amat(i,2) = aa*scale(i)*a2(i)
......@@ -67,29 +64,6 @@
amatinv(3,i) = (1.0/(aa*scale(i))) * b3(i)
ENDDO
ELSE ! we scale the lattice vectors in the case of mP, mA, mB, aP
DO i=1,3
amat(i,1) = aa*scale(1)*a1(i)
amat(i,2) = aa*scale(2)*a2(i)
amat(i,3) = aa*scale(3)*a3(i)
ENDDO
DO i=1,3
bmat(1,i) = (pi_const/(aa*scale(1))) * b1(i)
bmat(2,i) = (pi_const/(aa*scale(2))) * b2(i)
bmat(3,i) = (pi_const/(aa*scale(3))) * b3(i)
ENDDO
DO i=1,3
amatinv(1,i) = (1.0/(aa*scale(1))) * b1(i)
amatinv(2,i) = (1.0/(aa*scale(2))) * b2(i)
amatinv(3,i) = (1.0/(aa*scale(3))) * b3(i)
ENDDO
ENDIF
!---> check that amat and amatinv consistent
! (amat*amatinv should be identity)
......
inpgen/struct_input.f
View file @ c104ae6b
......@@ -9,7 +9,7 @@
> natmax,nop48,
X nline,xl_buffer,buffer,
< title,film,cal_symm,checkinp,symor,
< cartesian,oldfleur,a1,a2,a3,dvac,aa,scale,noangles,i_c,
< cartesian,oldfleur,a1,a2,a3,dvac,aa,scale,i_c,
< factor,natin,atomid,atompos,ngen,mmrot,ttr,
< l_hyb,l_soc,l_ss,theta,phi,qss,inistop)
......@@ -26,7 +26,7 @@
CHARACTER(len=xl_buffer) :: buffer
LOGICAL :: cal_symm, checkinp, symor, film
LOGICAL :: cartesian,oldfleur,inistop
LOGICAL, INTENT (OUT) :: l_hyb,l_soc,l_ss,noangles
LOGICAL, INTENT (OUT) :: l_hyb,l_soc,l_ss
INTEGER, INTENT (OUT) :: natin,i_c
INTEGER, INTENT (OUT) :: ngen
REAL, INTENT (OUT) :: aa,theta,phi
......@@ -46,7 +46,7 @@
!===> Local Variables
INTEGER :: n, ng, op, nbuffer, ios,nop2
REAL :: shift(3),rdummy(3,3),z_max,z_min
REAL :: shift(3),rdummy(3,3),z_max,z_min,mat(3,3),x(3)
LOGICAL :: oldfleurset,l_symfile,l_gen,hybrid
CHARACTER(len=10) :: chtmp
CHARACTER(len=3) :: ch_test
......@@ -140,7 +140,7 @@
IF ( buffer(1:8)=='&lattice' ) THEN
CALL lattice2(
> buffer,xl_buffer,errfh,bfh,nline,
< a1,a2,a3,aa,scale,noangles,i_c,ios )
< a1,a2,a3,aa,scale,mat,i_c,ios )
dvac = 0.00
IF ( ios.NE.0 ) THEN
WRITE (errfh,*)
......@@ -307,6 +307,20 @@
WRITE (warnfh,*)
ENDIF
IF (abs(mat(1,1)).GT.0.0000001) THEN ! transform hex->trig
CALL recip(a1,a2,a3,rdummy)
DO n = 1, abs(natin)
! CALL cotra0(atompos(1,n),x,mat)
x = matmul(mat,atompos(:,n))
write(*,'(3f10.5)') x(1:3)
write(*,'(3f10.5)') rdummy
! CALL cotra1(x,atompos(1,n),rdummy)
atompos(:,n) = matmul(rdummy,x)
write(*,'(3f10.5)') atompos(1:3,n)
ENDDO
i_c = 1
ENDIF
IF (film) THEN
z_max = MAXVAL( atompos(3,1:abs(natin)) ) ! check the outmost atomic position
......@@ -422,4 +436,32 @@
CALL juDFT_error("ERROR reading input",calledby="struct_input")
END SUBROUTINE struct_input
!-------------------------------------
SUBROUTINE recip(a1,a2,a3,b)
USE m_constants, ONLY : pimach
IMPLICIT NONE
REAL, INTENT (IN) :: a1(3),a2(3),a3(3)
REAL, INTENT (OUT):: b(3,3)
REAL volume
! volume (without scaling factor aa^3)
volume = a1(1)*a2(2)*a3(3) + a2(1)*a3(2)*a1(3) +
& a3(1)*a1(2)*a2(3) - a1(3)*a2(2)*a3(1) -
& a2(3)*a3(2)*a1(1) - a3(3)*a1(2)*a2(1)
! reciprocal lattice vectors in scaled Cartesian units
b(1,1) = (a2(2)*a3(3) - a2(3)*a3(2))
b(1,2) = (a2(3)*a3(1) - a2(1)*a3(3))
b(1,3) = (a2(1)*a3(2) - a2(2)*a3(1))
b(2,1) = (a3(2)*a1(3) - a3(3)*a1(2))
b(2,2) = (a3(3)*a1(1) - a3(1)*a1(3))
b(2,3) = (a3(1)*a1(2) - a3(2)*a1(1))
b(3,1) = (a1(2)*a2(3) - a1(3)*a2(2))
b(3,2) = (a1(3)*a2(1) - a1(1)*a2(3))
b(3,3) = (a1(1)*a2(2) - a1(2)*a2(1))
! b = 2.0*pimach()*b/volume
b = b/volume
END SUBROUTINE recip
END MODULE m_structinput
io/rw_symfile.f
View file @ c104ae6b
......@@ -23,7 +23,7 @@
LOGICAL, INTENT (INOUT) :: symor
!===> Variables
INTEGER i,j,n,ios,no2,no3,gen1
INTEGER i,j,n,ios,no2,no3,gen1,isrt(nop)
REAL t,d
LOGICAL ex,op
CHARACTER(len=3) :: type
......@@ -37,11 +37,28 @@
OPEN (symfh, file=sym2fn, status='unknown', err=911, iostat=ios)
WRITE (symfh,*) nop,nop2,symor,' ! nop,nop2,symor '
IF (nop == 2*nop2) THEN ! film-calculation
i = 1 ; j = nop2 + 1
DO n = 1, nop
IF (mrot(3,3,n) == 1) THEN
isrt(n) = i ; i = i + 1
ELSE
isrt(n) = j ; j = j + 1
ENDIF
ENDDO
ELSE
DO n = 1, nop
isrt(n) = n
ENDDO
ENDIF
DO n = 1, nop
WRITE (symfh,'(a1,i3)') '!', n
WRITE (symfh,'(3i5,5x,f10.5)')
& ((mrot(i,j,n),j=1,3),tau(i,n),i=1,3)
WRITE (symfh,'(3i5,5x,f10.5)')
& ((mrot(i,j,isrt(n)),j=1,3),tau(i,isrt(n)),i=1,3)
ENDDO
! WRITE (symfh,*) '! reciprocal lattice vectors'
! WRITE (symfh,'(3f25.15)') ((bmat(i,j),j=1,3),i=1,3)
......@@ -98,12 +115,7 @@
!
! Determine the kind of symmetry operation we have here
!
d = mrot(1,1,n)*mrot(2,2,n)*mrot(3,3,n) +
+ mrot(1,2,n)*mrot(2,3,n)*mrot(3,1,n) +
+ mrot(2,1,n)*mrot(3,2,n)*mrot(1,3,n) -
+ mrot(1,3,n)*mrot(2,2,n)*mrot(3,1,n) -
+ mrot(2,3,n)*mrot(3,2,n)*mrot(1,1,n) -
+ mrot(2,1,n)*mrot(1,2,n)*mrot(3,3,n)
d = det(mrot(:,:,n))
t = mrot(1,1,n) + mrot(2,2,n) + mrot(3,3,n)
IF (d.EQ.-1) THEN
......@@ -142,4 +154,15 @@
CALL juDFT_error("i/o ERROR",calledby="rw_symfile")
END SUBROUTINE rw_symfile
!--------------------------------------------------------------------
INTEGER FUNCTION det(m)
INTEGER m(3,3)
det = m(1,1)*m(2,2)*m(3,3) + m(1,2)*m(2,3)*m(3,1) +
+ m(2,1)*m(3,2)*m(1,3) - m(1,3)*m(2,2)*m(3,1) -
+ m(2,3)*m(3,2)*m(1,1) - m(2,1)*m(1,2)*m(3,3)
END FUNCTION det
!--------------------------------------------------------------------
END MODULE m_rwsymfile
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment