Commit ebc3d86f authored by Gregor Michalicek's avatar Gregor Michalicek

Beautify Julia

parent 426dfe3d
...@@ -3,13 +3,13 @@ ...@@ -3,13 +3,13 @@
! This file is part of FLEUR and available as free software under the conditions ! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail. ! of the MIT license as expressed in the LICENSE file in more detail.
!-------------------------------------------------------------------------------- !--------------------------------------------------------------------------------
MODULE m_julia
MODULE m_julia USE m_juDFT
use m_juDFT
CONTAINS CONTAINS
SUBROUTINE julia(&
& sym,cell,input,noco,banddos,& SUBROUTINE julia(sym,cell,input,noco,banddos,kpts,l_q,l_fillArrays)
& kpts,l_q,l_fillArrays)
!----------------------------------------------------------------------+ !----------------------------------------------------------------------+
! Generate a k-point file with approx. nkpt k-pts or a Monkhorst-Pack | ! Generate a k-point file with approx. nkpt k-pts or a Monkhorst-Pack |
! set with nmod(i) divisions in i=x,y,z direction. Interface to kptmop | ! set with nmod(i) divisions in i=x,y,z direction. Interface to kptmop |
...@@ -17,450 +17,411 @@ ...@@ -17,450 +17,411 @@
! G.B. 07/01 | ! G.B. 07/01 |
!----------------------------------------------------------------------+ !----------------------------------------------------------------------+
USE m_constants USE m_constants
USE m_bravais USE m_bravais
USE m_divi USE m_divi
USE m_brzone USE m_brzone
USE m_brzone2 USE m_brzone2
USE m_kptmop USE m_kptmop
USE m_kpttet USE m_kpttet
USE m_bandstr1 USE m_bandstr1
use m_types USE m_types
IMPLICIT NONE
TYPE(t_sym),INTENT(IN) :: sym IMPLICIT NONE
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_input),INTENT(IN) :: input TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_noco),INTENT(IN) :: noco TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_banddos),INTENT(IN) :: banddos TYPE(t_input), INTENT(IN) :: input
TYPE(t_kpts),INTENT(INOUT) :: kpts TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_banddos), INTENT(IN) :: banddos
LOGICAL, INTENT (IN) :: l_q, l_fillArrays TYPE(t_kpts), INTENT(INOUT) :: kpts
LOGICAL, INTENT (IN) :: l_q, l_fillArrays
INTEGER, PARAMETER :: nop48 = 48 INTEGER, PARAMETER :: nop48 = 48
INTEGER, PARAMETER :: mface = 51 INTEGER, PARAMETER :: mface = 51
INTEGER, PARAMETER :: mdir = 10 INTEGER, PARAMETER :: mdir = 10
INTEGER, PARAMETER :: nbsz = 3 INTEGER, PARAMETER :: nbsz = 3
INTEGER, PARAMETER :: ibfile = 6 INTEGER, PARAMETER :: ibfile = 6
INTEGER, PARAMETER :: nv48 = (2*nbsz+1)**3+48 INTEGER, PARAMETER :: nv48 = (2*nbsz+1)**3+48
INTEGER ndiv3 ! max. number of tetrahedrons (< 6*(kpts%nkpt+1) INTEGER ndiv3 ! max. number of tetrahedrons (< 6*(kpts%nkpt+1)
INTEGER ntet ! actual number of tetrahedrons INTEGER ntet ! actual number of tetrahedrons
REAL, ALLOCATABLE :: vkxyz(:,:) ! vector of kpoint generated; in cartesian representation
REAL, ALLOCATABLE :: wghtkp(:) ! associated with k-points for BZ integration REAL, ALLOCATABLE :: vkxyz(:,:) ! vector of kpoint generated; in cartesian representation
INTEGER, ALLOCATABLE :: ntetra(:,:) ! corners of the tetrahedrons REAL, ALLOCATABLE :: wghtkp(:) ! associated with k-points for BZ integration
REAL, ALLOCATABLE :: voltet(:) ! voulmes of the tetrahedrons INTEGER, ALLOCATABLE :: ntetra(:,:) ! corners of the tetrahedrons
REAL, ALLOCATABLE :: vktet(:,:) ! REAL, ALLOCATABLE :: voltet(:) ! voulmes of the tetrahedrons
REAL, ALLOCATABLE :: vktet(:,:)
REAL divis(4) ! Used to find more accurate representation of k-points
! vklmn(i,kpt)/divis(i) and weights as wght(kpt)/divis(4) REAL divis(4) ! Used to find more accurate representation of k-points
INTEGER nkstar ! number of stars for k-points generated in full stars ! vklmn(i,kpt)/divis(i) and weights as wght(kpt)/divis(4)
REAL bltv(3,3) ! cartesian Bravais lattice basis (a.u.) INTEGER nkstar ! number of stars for k-points generated in full stars
REAL rltv(3,3) ! reciprocal lattice basis (2\pi/a.u.) REAL bltv(3,3) ! cartesian Bravais lattice basis (a.u.)
REAL ccr(3,3,nop48) ! rotation matrices in cartesian repr. REAL rltv(3,3) ! reciprocal lattice basis (2\pi/a.u.)
REAL rlsymr(3,3,nop48) ! rotation matrices in reciprocal lattice basis representation REAL ccr(3,3,nop48) ! rotation matrices in cartesian repr.
REAL talfa(3,nop48) ! translation vector associated with (non-symmorphic) REAL rlsymr(3,3,nop48) ! rotation matrices in reciprocal lattice basis representation
! symmetry elements in Bravais lattice representation REAL talfa(3,nop48) ! translation vector associated with (non-symmorphic)
INTEGER ncorn,nedge,nface ! number of corners, faces and edges of the IBZ ! symmetry elements in Bravais lattice representation
REAL fnorm(3,mface) ! normal vector of the planes bordering the IBZ INTEGER ncorn,nedge,nface ! number of corners, faces and edges of the IBZ
REAL fdist(mface) ! distance vector of the planes bordering t IBZ REAL fnorm(3,mface) ! normal vector of the planes bordering the IBZ
REAL cpoint(3,mface) ! cartesian coordinates of corner points of IBZ REAL fdist(mface) ! distance vector of the planes bordering t IBZ
REAL xvec(3) ! arbitrary vector lying in the IBZ REAL cpoint(3,mface) ! cartesian coordinates of corner points of IBZ
REAL xvec(3) ! arbitrary vector lying in the IBZ
INTEGER idsyst ! crystal system identification in MDDFT programs
INTEGER idtype ! lattice type identification in MDDFT programs INTEGER idsyst ! crystal system identification in MDDFT programs
INTEGER idtype ! lattice type identification in MDDFT programs
INTEGER idimens ! number of dimensions for k-point set (2 or 3)
INTEGER nreg ! 1 kpoints in full BZ; 0 kpoints in irrBZ INTEGER idimens ! number of dimensions for k-point set (2 or 3)
INTEGER nfulst ! 1 kpoints ordered in full stars INTEGER nreg ! 1 kpoints in full BZ; 0 kpoints in irrBZ
! (meaningful only for nreg =1; full BZ) INTEGER nfulst ! 1 kpoints ordered in full stars
INTEGER nbound ! 0 no primary points on BZ boundary; ! (meaningful only for nreg =1; full BZ)
! 1 with boundary points (not for BZ integration!!!) INTEGER nbound ! 0 no primary points on BZ boundary;
INTEGER ikzero ! 0 no shift of k-points; ! 1 with boundary points (not for BZ integration!!!)
! 1 shift of k-points for better use of sym in irrBZ INTEGER ikzero ! 0 no shift of k-points;
REAL kzero(3) ! shifting vector to bring one k-point to or ! 1 shift of k-points for better use of sym in irrBZ
! away from (0,0,0) (for even/odd nkpt3) REAL kzero(3) ! shifting vector to bring one k-point to or
! away from (0,0,0) (for even/odd nkpt3)
INTEGER i,j,k,l,idiv,mkpt,addSym,nsym
INTEGER iofile,iokpt,kpri,ktest,kmidtet INTEGER i,j,k,l,idiv,mkpt,addSym,nsym
INTEGER idivis(3) INTEGER iofile,iokpt,kpri,ktest,kmidtet
LOGICAL random,trias INTEGER idivis(3)
REAL help(3),binv(3,3),rlsymr1(3,3),ccr1(3,3) LOGICAL random,trias
REAL help(3),binv(3,3),rlsymr1(3,3),ccr1(3,3)
random = .false. ! do not use random tetra-points
random = .false. ! do not use random tetra-points
!------------------------------------------------------------
! !------------------------------------------------------------
! idsyst idtype !
! ! idsyst idtype
! 1 cubic primitive !
! 2 tetragonal body centered ! 1 cubic primitive
! 3 orthorhombic face centered ! 2 tetragonal body centered
! 4 hexagonal A-face centered ! 3 orthorhombic face centered
! 5 trigonal B-face centered ! 4 hexagonal A-face centered
! 6 monoclinic C-face centered ! 5 trigonal B-face centered
! 7 triclinic ! 6 monoclinic C-face centered
! ! 7 triclinic
! ---> for 2 dimensions only the following Bravais lattices exist: !
! ! ---> for 2 dimensions only the following Bravais lattices exist:
! TYPE EQUIVALENT 3-DIM idsyst/idtype !
! square = p-tetragonal ( 1+2 axis ) 2/1 ! TYPE EQUIVALENT 3-DIM idsyst/idtype
! rectangular = p-orthorhomb ( 1+2 axis ) 3/1 ! square = p-tetragonal ( 1+2 axis ) 2/1
! centered rectangular = c-face-orthorhomb( 1+2 axis) 3/6 ! rectangular = p-orthorhomb ( 1+2 axis ) 3/1
! hexagonal = p-hexagonal ( 1+2 axis ) 4/1 ! centered rectangular = c-face-orthorhomb( 1+2 axis) 3/6
! oblique = p-monoclinic ( 1+2 axis ) 6/1 ! hexagonal = p-hexagonal ( 1+2 axis ) 4/1
! ! oblique = p-monoclinic ( 1+2 axis ) 6/1
!------------------------------------------------------------ !
!------------------------------------------------------------
IF(l_q) THEN
trias=input%tria IF(l_q) THEN
if (input%tria) call judft_error("tria=T not implemented for q-point generator",calledby='julia') trias=input%tria
!input%tria=.false. if (input%tria) call judft_error("tria=T not implemented for q-point generator",calledby='julia')
ENDIF !input%tria=.false.
ENDIF
IF (cell%latnam.EQ.'squ') THEN IF (cell%latnam.EQ.'squ') THEN
idsyst = 2 idsyst = 2
idtype = 1 idtype = 1
IF (.not.input%film) THEN IF (.not.input%film) THEN
IF (abs(cell%amat(1,1)-cell%amat(3,3)) < 0.0000001) THEN IF (abs(cell%amat(1,1)-cell%amat(3,3)) < 0.0000001) THEN
idsyst = 1 idsyst = 1
idtype = 1 idtype = 1
ENDIF END IF
ENDIF END IF
END IF
IF (cell%latnam.EQ.'p-r') THEN
idsyst = 3
idtype = 1
END IF
IF ((cell%latnam.EQ.'c-b').OR.(cell%latnam.EQ.'c-r')) THEN
idsyst = 3
idtype = 6
END IF
IF ((cell%latnam.EQ.'hex').OR.(cell%latnam.EQ.'hx3')) THEN
idsyst = 4
idtype = 1
END IF
IF (cell%latnam.EQ.'obl') THEN
idsyst = 6
idtype = 1
END IF END IF
IF (cell%latnam.EQ.'any') THEN IF (cell%latnam.EQ.'p-r') THEN
CALL bravais(& idsyst = 3
& cell%amat,& idtype = 1
& idsyst,idtype) END IF
ENDIF IF ((cell%latnam.EQ.'c-b').OR.(cell%latnam.EQ.'c-r')) THEN
nsym = sym%nop idsyst = 3
IF (input%film) nsym = sym%nop2 idtype = 6
! END IF
!-------------------- Want to make a Bandstructure ? -------- IF ((cell%latnam.EQ.'hex').OR.(cell%latnam.EQ.'hx3')) THEN
! idsyst = 4
IF (banddos%ndir == -4) THEN idtype = 1
CALL bandstr1(idsyst,idtype,cell%bmat,kpts,input,l_fillArrays,banddos) END IF
RETURN IF (cell%latnam.EQ.'obl') THEN
ENDIF idsyst = 6
! idtype = 1
!-------------------- Some variables we do not use ---------- END IF
! IF (cell%latnam.EQ.'any') THEN
iofile = 6 CALL bravais(cell%amat,idsyst,idtype)
iokpt = 6 END IF
kpri = 0 ! 3 nsym = sym%nop
ktest = 0 ! 5 IF (input%film) nsym = sym%nop2
kmidtet = 0
nreg = 0 ! Want to make a Bandstructure?
nfulst = 0 IF (banddos%ndir == -4) THEN
ikzero = 0 CALL bandstr1(idsyst,idtype,cell%bmat,kpts,input,l_fillArrays,banddos)
kzero(1) = 0.0 ; kzero(2) = 0.0 ; kzero(3) = 0.0 RETURN
nbound = 0 END IF
IF (input%tria) THEN
IF (input%film) nbound = 1 ! Some variables we do not use
! IF ((idsyst==1).AND.(idtype==1)) nbound = 1
! IF ((idsyst==2).AND.(idtype==1)) nbound = 1 iofile = 6
! IF ((idsyst==3).AND.(idtype==1)) nbound = 1 iokpt = 6
! IF ((idsyst==3).AND.(idtype==6)) nbound = 1 kpri = 0 ! 3
! IF ((idsyst==4).AND.(idtype==1)) nbound = 1 ktest = 0 ! 5
IF (nbound == 0) random = .true. kmidtet = 0
ENDIF nreg = 0
idimens = 3 nfulst = 0
IF (input%film) idimens = 2 ikzero = 0
! kzero(1) = 0.0
!--------------------- Lattice information ------------------ kzero(2) = 0.0
kzero(3) = 0.0
DO j = 1,3 nbound = 0
DO k = 1,3 IF (input%tria) THEN
bltv(j,k) = cell%amat(k,j) IF (input%film) nbound = 1
binv(j,k) = cell%bmat(k,j)/tpi_const ! IF ((idsyst==1).AND.(idtype==1)) nbound = 1
rltv(j,k) = cell%bmat(k,j) ! IF ((idsyst==2).AND.(idtype==1)) nbound = 1
DO i = 1,nsym ! IF ((idsyst==3).AND.(idtype==1)) nbound = 1
rlsymr(k,j,i) = real( sym%mrot(j,k,i) ) ! IF ((idsyst==3).AND.(idtype==6)) nbound = 1
ENDDO ! IF ((idsyst==4).AND.(idtype==1)) nbound = 1
ENDDO IF (nbound == 0) random = .true.
ENDDO END IF
idimens = 3
ccr = 0.0 IF (input%film) idimens = 2
DO i = 1,nsym
DO j = 1,3 ! Lattice information
talfa(j,i) = 0.0
DO k = 1,3 DO j = 1, 3
DO k = 1, 3
bltv(j,k) = cell%amat(k,j)
binv(j,k) = cell%bmat(k,j) / tpi_const
rltv(j,k) = cell%bmat(k,j)
DO i = 1,nsym
rlsymr(k,j,i) = real(sym%mrot(j,k,i))
END DO
END DO
END DO
ccr = 0.0
DO i = 1, nsym
DO j = 1, 3
talfa(j,i) = 0.0
DO k = 1, 3
talfa(j,i) = bltv(j,k) * sym%tau(k,i) talfa(j,i) = bltv(j,k) * sym%tau(k,i)
help(k) = 0.0 help(k) = 0.0
DO l = 1,3 DO l = 1, 3
help(k) = help(k) + rlsymr(l,k,i) * binv(j,l) help(k) = help(k) + rlsymr(l,k,i) * binv(j,l)
ENDDO END DO
ENDDO END DO
DO k = 1,3 DO k = 1, 3
ccr(j,k,i) = 0.0 ccr(j,k,i) = 0.0
DO l = 1,3 DO l = 1, 3
ccr(j,k,i) = ccr(j,k,i) + bltv(l,k) * help(l) ccr(j,k,i) = ccr(j,k,i) + bltv(l,k) * help(l)
ENDDO END DO
ENDDO END DO
ENDDO END DO
! write (*,'(3f12.6)') ((ccr(j,k,i),j=1,3),k=1,3) END DO
! write (*,*) DO i = 1, nsym
ENDDO rlsymr1(:,:) = rlsymr(:,:,i)
DO i = 1,nsym ccr1(:,:) = ccr(:,:,i)
rlsymr1(:,:) = rlsymr(:,:,i) DO j = 1, 3
ccr1(:,:) = ccr(:,:,i) DO k = 1, 3
DO j = 1,3
DO k = 1,3
rlsymr(k,j,i) = rlsymr1(j,k) rlsymr(k,j,i) = rlsymr1(j,k)
ccr(k,j,i) = ccr1(j,k) ccr(k,j,i) = ccr1(j,k)
ENDDO END DO
ENDDO END DO
ENDDO END DO
IF ((.not.noco%l_ss).AND.(.not.noco%l_soc).AND.(2*nsym<nop48)) THEN IF ((.not.noco%l_ss).AND.(.not.noco%l_soc).AND.(2*nsym<nop48)) THEN
IF ((input%film.AND.(.not.sym%invs2)).OR.((.not.input%film).AND.(.not.sym%invs))) THEN
IF ( (input%film.AND.(.not.sym%invs2)).OR.& addSym = 0
& ((.not.input%film).AND.(.not.sym%invs)) ) THEN ! Note: We have to add the negative of each symmetry operation
addSym = 0 ! to exploit time reversal symmetry. However, if the new
! Note: We have to add the negative of each symmetry operation ! symmetry operation is the identity matrix it is excluded.
! to exploit time reversal symmetry. However, if the new ! This is the case iff it is (-Id) + a translation vector.
! symmetry operation is the identity matrix it is excluded. DO i = 1, nsym
! This is the case iff it is (-Id) + a translation vector. ! This test assumes that ccr(:,:,1) is the identity matrix.
DO i = 1, nsym IF(.NOT.ALL(ABS(ccr(:,:,1)+ccr(:,:,i)).LT.10e-10) ) THEN
! This test assumes that ccr(:,:,1) is the identity matrix. ccr(:,:,nsym+addSym+1 ) = -ccr(:,:,i)
IF(.NOT.ALL(ABS(ccr(:,:,1)+ccr(:,:,i)).LT.10e-10) ) THEN rlsymr(:,:,nsym+addSym+1 ) = -rlsymr(:,:,i)
ccr(:,:,nsym+addSym+1 ) = -ccr(:,:,i) addSym = addSym + 1
rlsymr(:,:,nsym+addSym+1 ) = -rlsymr(:,:,i) END IF
addSym = addSym + 1 END DO
END IF nsym = nsym + addSym
END DO END IF
nsym = nsym + addSym END IF
ENDIF
! brzone and brzone2 find the corner-points, the edges, and the
ENDIF ! faces of the irreducible wedge of the brillouin zone (IBZ).
! In these subroutines many special cases can occur. Due to this the very
! brzone and brzone2 find the corner-points, the edges, and the ! sophisticated old routine brzone had a few bugs. The new routine
! faces of the irreducible wedge of the brillouin zone (IBZ). ! brzone2 was written with a different algorithm that is slightly slower
! In these subroutines many special cases can occur. Due to this the very ! but should be more stable. To make comparisons possible the old
! sophisticated old routine brzone had a few bugs. The new routine ! routine is only commented out. Both routines are directly
! brzone2 was written with a different algorithm that is slightly slower ! interchangable. GM, 2016.
! but should be more stable. To make comparisons possible the old
! routine is only commented out. Both routines are directly ! CALL brzone(rltv,nsym,ccr,mface,nbsz,nv48,cpoint,xvec,ncorn,nedge,nface,fnorm,fdist)
! interchangable. GM, 2016.
CALL brzone2(rltv,nsym,ccr,mface,nbsz,nv48,cpoint,xvec,ncorn,nedge,nface,fnorm,fdist)
! CALL brzone(&
! & rltv,nsym,ccr,mface,nbsz,nv48,& IF (input%tria.AND.random) THEN
! & cpoint,& ! Calculate the points for tetrahedron method
! & xvec,ncorn,nedge,nface,fnorm,fdist) mkpt = kpts%nkpt
ndiv3 = 6*(mkpt+1)
CALL brzone2(& ALLOCATE (vkxyz(3,mkpt),wghtkp(mkpt))
& rltv,nsym,ccr,mface,nbsz,nv48,& ALLOCATE (voltet(ndiv3),vktet(3,mkpt),ntetra(4,ndiv3))
& cpoint,& vkxyz = 0.0
& xvec,ncorn,nedge,nface,fnorm,fdist) CALL kpttet(iofile,ibfile,iokpt,kpri,ktest,kmidtet,mkpt,ndiv3,&
nreg,nfulst,rltv,cell%omtil,nsym,ccr,mdir,mface,&
ncorn,nface,fdist,fnorm,cpoint,voltet,ntetra,ntet,vktet,&
IF ( input%tria.AND.random ) THEN kpts%nkpt,divis,vkxyz,wghtkp)
! ELSE
! Calculate the points for tetrahedron method ! If just the total number of k-points is given, determine
! ! the divisions in each direction (nkpt3):
mkpt = kpts%nkpt
ndiv3 = 6*(mkpt+1) ! IF (tria) THEN
ALLOCATE (vkxyz(3,mkpt),wghtkp(mkpt) ) ! nkpt = nkpt/4
ALLOCATE ( voltet(ndiv3),vktet(3,mkpt),ntetra(4,ndiv3) ) ! nkpt3(:) = nkpt3(:) / 2
vkxyz = 0.0 ! END IF
CALL kpttet(& IF (sum(kpts%nkpt3).EQ.0) THEN
& iofile,ibfile,iokpt,& CALL divi(kpts%nkpt,cell%bmat,input%film,sym%nop,sym%nop2,kpts%nkpt3)
& kpri,ktest,kmidtet,mkpt,ndiv3,& END IF
& nreg,nfulst,rltv,cell%omtil,&
& nsym,ccr,mdir,mface,& ! Now calculate Monkhorst-Pack k-points:
& ncorn,nface,fdist,fnorm,cpoint,& IF (kpts%nkpt3(2).EQ.0) kpts%nkpt3(2) = kpts%nkpt3(1)
& voltet,ntetra,ntet,vktet,& IF ((.not.input%film).AND.(kpts%nkpt3(3).EQ.0)) kpts%nkpt3(3) = kpts%nkpt3(2)
& kpts%nkpt,& IF (nbound.EQ.1) THEN
& divis,vkxyz,wghtkp) mkpt = (2*kpts%nkpt3(1)+1)*(2*kpts%nkpt3(2)+1)
IF (.not.input%film) mkpt = mkpt*(2*kpts%nkpt3(3)+1)
ELSE ELSE
! mkpt = kpts%nkpt3(1)*kpts%nkpt3(2)
! If just the total number of k-points is given, determine IF (.not.input%film) mkpt = mkpt*kpts%nkpt3(3)
! the divisions in each direction (nkpt3): END IF
! ALLOCATE (vkxyz(3,mkpt),wghtkp(mkpt) )
! IF (tria) THEN vkxyz = 0.0
! nkpt = nkpt/4
! nkpt3(:) = nkpt3(:) / 2 CALL kptmop(iofile,iokpt,kpri,ktest,idsyst,idtype,kpts%nkpt3,ikzero,kzero,&
! ENDIF rltv,bltv,nreg,nfulst,nbound,idimens,xvec,fnorm,fdist,ncorn,nface,&
IF (sum(kpts%nkpt3).EQ.0) THEN nedge,cpoint,nsym,ccr,rlsymr,talfa,mkpt,mface,mdir,&
CALL divi(& kpts%nkpt,divis,vkxyz,nkstar,wghtkp)
& kpts%nkpt,cell%bmat,input%film,sym%nop,sym%nop2,& END IF
& kpts%nkpt3)
ENDIF idivis(1) = int(divis(1))
idivis(2) = int(divis(2))
! idivis(3) = int(divis(3))
! Now calculate Monkhorst-Pack k-points: idiv = lcm(3,idivis)
! IF (idiv.GE.200) idiv = 1
IF (kpts%nkpt3(2).EQ.0) kpts%nkpt3(2) = kpts%nkpt3(1) DO j=1,kpts%nkpt
IF ((.not.input%film).AND.(kpts%nkpt3(3).EQ.0)) kpts%nkpt3(3) = kpts%nkpt3(2) wghtkp(j) = wghtkp(j) * divis(4)
IF (nbound.EQ.1) THEN DO k = 1,3
mkpt = (2*kpts%nkpt3(1)+1)*(2*kpts%nkpt3(2)+1) help(k) = 0.0
IF (.not.input%film) mkpt = mkpt*(2*kpts%nkpt3(3)+1) DO l = 1,3
ELSE help(k) = help(k) + cell%amat(l,k) * vkxyz(l,j)
mkpt = kpts%nkpt3(1)*kpts%nkpt3(2) END DO
IF (.not.input%film) mkpt = mkpt*kpts%nkpt3(3) END DO
ENDIF DO i=1,3
ALLOCATE (vkxyz(3,mkpt),wghtkp(mkpt) ) vkxyz(i,j) = help(i) * idiv / tpi_const
vkxyz = 0.0 END DO
END DO
CALL kptmop(&
& iofile,iokpt,kpri,ktest,& ! if (l_q) write qpts file:
& idsyst,idtype,kpts%nkpt3,ikzero,kzero,& IF(l_q)THEN
& rltv,bltv,nreg,nfulst,nbound,idimens,& IF(input%film) THEN
& xvec,fnorm,fdist,ncorn,nface,nedge,cpoint,& CALL juDFT_error("For the case of input%film q-points generator not implemented!", calledby = "julia")
& nsym,ccr,rlsymr,talfa,mkpt,mface,mdir,& END IF
& kpts%nkpt,divis,vkxyz,nkstar,wghtkp)
ENDIF
!
idivis(1) = int(divis(1))
idivis(2) = int(divis(2))
idivis(3) = int(divis(3))
idiv = lcm(3,idivis)
! WRITE (*,'(2i5)') nkpt,idiv
IF (idiv.GE.200) idiv = 1
DO j=1,kpts%nkpt
! WRITE (*,'(4f10.5)') (vkxyz(i,j),i=1,3),wghtkp(j)
wghtkp(j) = wghtkp(j) * divis(4)
DO k = 1,3
help(k) = 0.0
DO l = 1,3
help(k) = help(k) + cell%amat(l,k) * vkxyz(l,j)
ENDDO
ENDDO
DO i=1,3
vkxyz(i,j) = help(i) * idiv / tpi_const
ENDDO
ENDDO
!
! if (l_q) write qpts file:
!
IF(l_q)THEN
IF(input%film) CALL juDFT_error("For the case of input%film q-points "//&
& "generator not implemented!",calledby ="julia")
OPEN(113,file='qpts',form='formatted',status='new') OPEN(113,file='qpts',form='formatted',status='new')
WRITE(113,'(i5)') kpts%nkpt+1 WRITE(113,'(i5)') kpts%nkpt+1
WRITE(113,8050) 0.,0.,0. WRITE(113,8050) 0.,0.,0.
DO j = 1, kpts%nkpt DO j = 1, kpts%nkpt
WRITE (113,FMT=8050) (vkxyz(i,j)/real(idiv),i=1,3) WRITE (113,FMT=8050) (vkxyz(i,j)/real(idiv),i=1,3)
ENDDO END DO
CLOSE(113) CLOSE(113)
!input%tria=trias !input%tria=trias
RETURN RETURN
ENDIF END IF
8050 FORMAT (2(f14.10,1x),f14.10) 8050 FORMAT (2(f14.10,1x),f14.10)
! ! write k-points file or write data into arrays
! write k-points file or write data into arrays IF (l_fillArrays) THEN
! IF (ALLOCATED(kpts%bk)) THEN
DEALLOCATE(kpts%bk)
IF (l_fillArrays) THEN END IF
IF (ALLOCATED(kpts%bk)) THEN IF (ALLOCATED(kpts%wtkpt)) THEN
DEALLOCATE(kpts%bk) DEALLOCATE(kpts%wtkpt)
END IF
ALLOCATE(kpts%bk(3,kpts%nkpt),kpts%wtkpt(kpts%nkpt))
IF (idiv.NE.0) kpts%posScale = REAL(idiv)
DO j = 1, kpts%nkpt
kpts%bk(1,j) = vkxyz(1,j)
kpts%bk(2,j) = vkxyz(2,j)
kpts%bk(3,j) = vkxyz(3,j)
kpts%wtkpt(j) = wghtkp(j)
END DO
IF (input%tria.AND.random) THEN
kpts%ntet = ntet
IF (ALLOCATED(kpts%ntetra)) THEN
DEALLOCATE(kpts%ntetra)
END IF END IF
IF (ALLOCATED(kpts%wtkpt)) THEN IF (ALLOCATED(kpts%voltet)) THEN
DEALLOCATE(kpts%wtkpt) DEALLOCATE(kpts%voltet)
END IF END IF
ALLOCATE(kpts%bk(3,kpts%nkpt),kpts%wtkpt(kpts%nkpt)) ALLOCATE(kpts%ntetra(4,kpts%ntet))
IF (idiv.NE.0) kpts%posScale = REAL(idiv) ALLOCATE(kpts%voltet(kpts%ntet))
DO j = 1, ntet
DO i = 1, 4
kpts%ntetra(i,j) = ntetra(i,j)
END DO
kpts%voltet(j) = ABS(voltet(j))
END DO
END IF
ELSE
OPEN (41,file='kpts',form='formatted',status='new')
IF (input%film) THEN
WRITE (41,FMT=8110) kpts%nkpt,real(idiv),.false.
DO j = kpts%nkpt, 1, -1
WRITE (41,FMT=8040) (vkxyz(i,j),i=1,2),wghtkp(j)
END DO
ELSE
WRITE (41,FMT=8100) kpts%nkpt,real(idiv)
DO j = 1, kpts%nkpt DO j = 1, kpts%nkpt
kpts%bk(1,j) = vkxyz(1,j) WRITE (41,FMT=8040) (vkxyz(i,j),i=1,3),wghtkp(j)
kpts%bk(2,j) = vkxyz(2,j)
kpts%bk(3,j) = vkxyz(3,j)
kpts%wtkpt(j) = wghtkp(j)