Commit 59c1cbfe authored by Gregor Michalicek's avatar Gregor Michalicek

Introduce oUnit to files in wannier directory

parent eaab3a04
......@@ -58,7 +58,7 @@ c Inclusion of Noco&&soc; tidied up version
c Frank Freimuth
c*********************************************************************
use m_intgr, only : intgr3
use m_constants, only: pimach, ImagUnit
use m_constants
use m_dwigner
use m_wann_tlmw
use m_wann_rad_twf
......@@ -181,15 +181,15 @@ c..reading the proj.1 / proj.2 / proj file
close (203)
if (ikpt.eq.1) then
write (6,*) 'Number of trial WFs:',nwfs
write (6,*)
write (oUnit,*) 'Number of trial WFs:',nwfs
write (oUnit,*)
do nwf = 1,nwfs
write (6,*) 'Twfs N:',nwf,' Atom N:',ind(nwf)
write (6,*) 'l=',lwf(nwf),' mr=',mrwf(nwf),' r=',rwf(nwf)
write (6,*) 'zona=',zona(nwf),' region=',regio(nwf),'*Rmt'
write (6,*) 'alpha=',alpha(nwf),
write (oUnit,*) 'Twfs N:',nwf,' Atom N:',ind(nwf)
write (oUnit,*) 'l=',lwf(nwf),' mr=',mrwf(nwf),' r=',rwf(nwf)
write (oUnit,*) 'zona=',zona(nwf),' region=',regio(nwf),'*Rmt'
write (oUnit,*) 'alpha=',alpha(nwf),
& ' beta=',beta(nwf),' gamma=',gamma(nwf)
write (6,*)
write (oUnit,*)
enddo
endif
......
......@@ -28,8 +28,12 @@ c host atoms.
c
c Frank Freimuth
c*********************************************
USE m_constants
use m_wann_readcenters
implicit none
integer, intent(in) :: natd
real, intent(in) :: pos(3,natd)
real, intent(in) :: omtil
......@@ -90,7 +94,7 @@ c--------reading the proj.1 / proj.2 / proj file
close(203)
endif
print*,"number of wannier functions= ",num_wann(jspin)
write(6,*)"number of wannier functions",num_wann(jspin)
write(oUnit,*)"number of wannier functions",num_wann(jspin)
allocate( wann_centers(3,num_wann(jspin)) )
call wann_readcenters(
> num_wann(jspin),jspin,
......
......@@ -20,7 +20,7 @@ c
c Frank Freimuth, February 2011
c*************************************************
use m_constants, only:pimach
use m_constants
use m_wann_read_umatrix
implicit none
......@@ -95,7 +95,7 @@ c real :: kpoints(3,nkpts)
jspins=jspins_in
if(l_soc)jspins=1
write(6,*)"nkpts=",nkpts
write(oUnit,*)"nkpts=",nkpts
c*****************************************************
c get num_bands and num_wann from the proj file
c*****************************************************
......@@ -114,8 +114,8 @@ c*****************************************************
endif
read (203,*) num_wann,num_bands
close (203)
write(6,*)'According to proj there are ',num_bands,' bands'
write(6,*)"and ",num_wann," wannier functions."
write(oUnit,*)'According to proj there are ',num_bands,' bands'
write(oUnit,*)"and ",num_wann," wannier functions."
c****************************************************************
c read in chk
......@@ -181,7 +181,7 @@ c Calculate matrix elements of SOC in the basis of
c rotated Bloch functions.
c****************************************************************
allocate( hsomtx2(num_wann,num_wann,3,nkpts) )
write(6,*)"calculate matrix elements of SOC commutator
write(oUnit,*)"calculate matrix elements of SOC commutator
&between wannier orbitals"
if(have_disentangled) then
......@@ -232,7 +232,7 @@ c****************************************************************
c************************************************************
c Calculate matrix elements in real space.
c***********************************************************
write(6,*)"calculate SOC-mat in rs"
write(oUnit,*)"calculate SOC-mat in rs"
allocate(hreal(num_wann,num_wann,3,rvecnum))
hreal=cmplx(0.0,0.0)
do rvecind=1,rvecnum
......
......@@ -22,7 +22,7 @@ c
c Frank Freimuth, February 2011
c*************************************************
use m_constants, only:pimach
use m_constants
use m_wann_read_umatrix
implicit none
......@@ -97,7 +97,7 @@ c*************************************************
tpi=2*pimach()
write(6,*)"nkpts=",nkpts
write(oUnit,*)"nkpts=",nkpts
c*****************************************************
c get num_bands and num_wann from the proj file
c*****************************************************
......@@ -116,8 +116,8 @@ c*****************************************************
endif
read (203,*) num_wann,num_bands
close (203)
write(6,*)'According to proj there are ',num_bands,' bands'
write(6,*)"and ",num_wann," wannier functions."
write(oUnit,*)'According to proj there are ',num_bands,' bands'
write(oUnit,*)"and ",num_wann," wannier functions."
c****************************************************************
c read in chk
......@@ -178,7 +178,7 @@ c****************************************************************
c fft5: transformation to Bloch-like
c****************************************************************
allocate( hsomtx2(dim1,dim2,num_wann,num_wann,nkpts) )
write(6,*)"fft5: transformation to Bloch-like"
write(oUnit,*)"fft5: transformation to Bloch-like"
if(have_disentangled) then
hsomtx2=0.0
......@@ -230,7 +230,7 @@ c****************************************************************
c************************************************************
c Calculate matrix elements in real space.
c***********************************************************
write(6,*)"transform to rs"
write(oUnit,*)"transform to rs"
allocate(hreal(dim1,dim2,num_wann,num_wann,rvecnum))
hreal=cmplx(0.0,0.0)
do rvecind=1,rvecnum
......
......@@ -16,7 +16,11 @@ c Read in the k-points from kpts/w90kpts file.
c
c Frank Freimuth
c********************************************************
USE m_constants
implicit none
TYPE(t_input), INTENT(IN) :: input
TYPE(t_kpts), INTENT(IN) :: kpts
logical,intent(in) :: l_bzsym,film
......@@ -35,7 +39,7 @@ c********************************************************
+ ="wann_get_kpts")
open(987,file='w90kpts',status='old',form='formatted')
read(987,*)nkpts, scale
write(6,*)"wann_get_kpts: nkpts=",nkpts
write(oUnit,*)"wann_get_kpts: nkpts=",nkpts
if(l_readkpts)then
IF(SIZE(kpoints,1)/=3) CALL juDFT_error("wann_get_kpts: 1"
+ ,calledby ="wann_get_kpts")
......@@ -48,7 +52,7 @@ c********************************************************
close(987)
else
nkpts = kpts%nkpt
write(6,*)"wann_get_kpts: nkpts=",nkpts
write(oUnit,*)"wann_get_kpts: nkpts=",nkpts
if(l_readkpts)then
do iter=1,nkpts
kpoints(:,iter) = kpts%bk(:,iter)
......@@ -58,7 +62,7 @@ c********************************************************
IF (l_readkpts) THEN
do iter=1,nkpts
write(6,*)kpoints(:,iter)
write(oUnit,*)kpoints(:,iter)
enddo
END IF
......
......@@ -16,7 +16,11 @@ c Monkhorst-Pack mesh.
c
c Frank Freimuth
c**************************************
USE m_constants
implicit none
integer,intent(in) :: nkpts
real,intent(in) :: kpoints(:,:)
integer,intent(out) :: num(:)
......@@ -48,8 +52,9 @@ c**************************************
num(dim)=(maxi-mini)/increm+1.01
endif
enddo
write(6,*)"wann_get_mp: determination of mp-grid parameters:"
write(6,*)"mp_1=",num(1),"mp_2=",num(2),"mp_3=",num(3)
write(oUnit,*)
+ "wann_get_mp: determination of mp-grid parameters:"
write(oUnit,*)"mp_1=",num(1),"mp_2=",num(2),"mp_3=",num(3)
IF(num(1)*num(2)*num(3)/=nkpts) CALL juDFT_error
+ ("mysterious kpoints",calledby ="wann_get_mp")
......
......@@ -9,7 +9,11 @@ c Read in the q-points from qpts file.
c
c
c********************************************************
USE m_constants
implicit none
logical,intent(in) :: l_bzsym,film
logical,intent(in) :: l_onedimens,l_readqpts
integer,intent(out) :: nqpts
......@@ -27,7 +31,7 @@ c********************************************************
+ ="wann_get_qpts")
open(987,file='w90qpts',status='old',form='formatted')
read(987,*)nqpts, scale
write(6,*)"wann_get_qpts: nqpts=",nqpts
write(oUnit,*)"wann_get_qpts: nqpts=",nqpts
if(l_readqpts)then
IF(SIZE(qpoints,1)/=3) CALL juDFT_error("wann_get_qpts: 1"
+ ,calledby ="wann_get_qpts")
......@@ -45,7 +49,7 @@ c********************************************************
+ ="wann_get_qpts")
open(987,file=param_file,status='old',form='formatted')
read(987,*)nqpts,scale
write(6,*)"wann_get_qpts: nqpts=",nqpts
write(oUnit,*)"wann_get_qpts: nqpts=",nqpts
if(l_readqpts)then
IF(SIZE(qpoints,1)/=3) CALL juDFT_error("wann_get_qpts: 1"
+ ,calledby ="wann_get_qpts")
......@@ -66,7 +70,7 @@ c********************************************************
qpoints(3,:)=0.0
endif
do iter=1,nqpts
write(6,*)qpoints(:,iter)
write(oUnit,*)qpoints(:,iter)
enddo
endif
......
......@@ -19,7 +19,7 @@ c in file WF1.chk produced by wannier90.
c
c Frank Freimuth
c****************************************************
use m_constants, only:pimach
use m_constants
use m_wann_read_umatrix
c use m_wann_get_mp
c use m_wann_get_kpts
......@@ -105,7 +105,7 @@ c real :: kpoints(3,nkpts)
if(l_soc)jspins=1
l_socham=.false.
write(6,*)"nkpts=",nkpts
write(oUnit,*)"nkpts=",nkpts
do jspin=1,jspins !spin loop
c*****************************************************
......@@ -127,8 +127,9 @@ c*****************************************************
endif
read (203,*) num_wann,num_bands
close (203)
write(6,*)'According to proj there are ',num_bands,' bands'
write(6,*)"and ",num_wann," wannier functions."
write(oUnit,*)'According to proj there are ',
+ num_bands,' bands'
write(oUnit,*)"and ",num_wann," wannier functions."
c****************************************************************
c read in chk
c****************************************************************
......@@ -150,7 +151,7 @@ c****************************************************************
if(l_soc.and.l_socmmn0)then
num_bands2=neigd
endif
write(6,*)"read in eig-file"
write(oUnit,*)"read in eig-file"
allocate(energy(num_bands2,num_kpts))
inquire(file=spin12(jspin)//'.eig',exist=l_umdat)
IF(.NOT.l_umdat) CALL juDFT_error
......@@ -166,8 +167,8 @@ c****************************************************************
minenerg=minval(energy(:,:))
maxenerg=maxval(energy(:,:))
write(6,*)"minenerg=",minenerg
write(6,*)"maxenerg=",maxenerg
write(oUnit,*)"minenerg=",minenerg
write(oUnit,*)"maxenerg=",maxenerg
c*********************************************
c Preparations for spin-orbit coupling
c*********************************************
......@@ -204,7 +205,7 @@ c****************************************************************
c calculate matrix elements of hamiltonian
c****************************************************************
write(6,*)"calculate matrix elements of hamiltonian
write(oUnit,*)"calculate matrix elements of hamiltonian
& between wannier orbitals"
......@@ -276,9 +277,9 @@ c****************************************************************
num_wann2=neigd
wann_shift=band_min(jspin)-1
endif
write(6,*)"num_wann2=",num_wann2
write(6,*)"num_wann=",num_wann
write(6,*)"wann_shift=",wann_shift
write(oUnit,*)"num_wann2=",num_wann2
write(oUnit,*)"num_wann=",num_wann
write(oUnit,*)"wann_shift=",wann_shift
allocate(hwannsoc(num_wann,num_wann,num_kpts,2,2))
hwannsoc=cmplx(0.0,0.0)
do kspin=1,2
......@@ -312,7 +313,7 @@ c****************************************************************
& hwannsoc(i,j,k,2,1)+hwannsoc(i,j,k,2,2)
eulav=eulav-hwann(i,j,k)
if(abs(eulav).gt.1.e-6)then
write(6,*)"soc-hop: something wrong:",eulav
write(oUnit,*)"soc-hop: something wrong:",eulav
endif
enddo
enddo
......@@ -331,7 +332,8 @@ c***************************************************************
if(repro_eig)then
write(6,*)"As a check, try to reproduce the old eigenvalues"
write(oUnit,*)
+ "As a check, try to reproduce the old eigenvalues"
c Note, that this check will give positive result,
c if u_matrix is unitary
......@@ -368,11 +370,11 @@ c if u_matrix is unitary
do j = 1,num_wann
eigw(i,j) = ei(j)
if(abs(eigw(i,j)-eigval2(j,i)).gt.0.0001)then
write(6,*)"found different eigenvalues:"
write(6,*)"kpt=",i
write(6,*)"band=",j
write(6,*)"eig=",eigval2(j,i)
write(6,*)"neweig=",eigw(i,j)
write(oUnit,*)"found different eigenvalues:"
write(oUnit,*)"kpt=",i
write(oUnit,*)"band=",j
write(oUnit,*)"eig=",eigval2(j,i)
write(oUnit,*)"neweig=",eigw(i,j)
endif
enddo
......@@ -391,7 +393,7 @@ c if u_matrix is unitary
c************************************************************
c Calculate hoppings.
c***********************************************************
write(6,*)"calculate hoppings"
write(oUnit,*)"calculate hoppings"
allocate( hreal(num_wann,num_wann,rvecnum) )
hreal=cmplx(0.0,0.0)
do rvecind=1,rvecnum
......
......@@ -14,7 +14,11 @@ c*****************************************************************
c Apply the symmetries to reduce the number of k-points.
c Frank Freimuth
c*****************************************************************
USE m_constants
implicit none
integer, intent(in) :: nop
logical,intent(in) :: film
real, intent(in) :: bmat(3,3)
......@@ -44,7 +48,7 @@ c*****************************************************************
inquire(file='onlysymor',exist=l_onlysymor)
bbmat=matmul(bmat,transpose(bmat))
write(6,*) "Apply the symmetries to w90kpts"
write(oUnit,*) "Apply the symmetries to w90kpts"
c**********************************************************
c read in kpoints from w90kpts file
c**********************************************************
......@@ -62,9 +66,9 @@ c**********************************************************
read(987,*)kpoints(:,iter)
enddo
close(987)
write(6,*) "kpoints read from w90kpts:"
write(oUnit,*) "kpoints read from w90kpts:"
do iter=1,nkpts
write(6,'(3f10.5)')kpoints(:,iter)/scale
write(oUnit,'(3f10.5)')kpoints(:,iter)/scale
enddo
c*********************************************************
......@@ -160,9 +164,9 @@ c****************************************************
bkpt(:)=kpoints(:,ikpt)
if (mapk(ikpt)==0)then
sumweights=sumweights+weight(ikpt)
write(6,*)"ikpt=",ikpt
write(6,*)"irreducible"
write(6,fmt='(a10,3f9.6)')"internal: ",bkpt(:)/scale
write(oUnit,*)"ikpt=",ikpt
write(oUnit,*)"irreducible"
write(oUnit,fmt='(a10,3f9.6)')"internal: ",bkpt(:)/scale
if (film.and..not.l_onedimens)then
write(119,'(3f10.5)')bkpt(1:2),weight(ikpt)
......@@ -171,32 +175,34 @@ c****************************************************
endif
elseif(mapkoper(ikpt).gt.0)then
write(6,*)"ikpt=",ikpt
write(6,*)"reducible"
write(6,*)"map=",mapk(ikpt)
write(oUnit,*)"ikpt=",ikpt
write(oUnit,*)"reducible"
write(oUnit,*)"map=",mapk(ikpt)
brot(:)=0.0
do k=1,3
brot(:)=brot(:)+
+ mrot(k,:,mapkoper(ikpt))*kpoints(k,mapk(ikpt))
enddo
write(6,'(a19,3f9.6)')"rotated internal: ",brot(:)/scale
write(oUnit,'(a19,3f9.6)')"rotated internal: ",
+ brot(:)/scale
elseif(mapkoper(ikpt).lt.0)then
write(6,*)"ikpt=",ikpt
write(6,*)"reducible"
write(6,*)"map=",mapk(ikpt)
write(oUnit,*)"ikpt=",ikpt
write(oUnit,*)"reducible"
write(oUnit,*)"map=",mapk(ikpt)
brot(:)=0.0
do k=1,3
brot(:)=brot(:)-
+ mrot(k,:,-mapkoper(ikpt))*kpoints(k,mapk(ikpt))
enddo
write(6,'(a19,3f9.6)')"rotated internal: ",brot(:)/scale
write(oUnit,'(a19,3f9.6)')"rotated internal: ",
+ brot(:)/scale
endif
enddo
close(119)
write(6,*)"reduznumk=",reduznumk
write(6,*)"nkpts=",nkpts
write(6,*)"sumweights=",sumweights
write(oUnit,*)"reduznumk=",reduznumk
write(oUnit,*)"nkpts=",nkpts
write(oUnit,*)"sumweights=",sumweights
IF(sumweights/=nkpts) CALL juDFT_error
+ ("sum of weights differs from nkpts",calledby
......
......@@ -15,7 +15,11 @@ c*****************************************************************
c Apply the symmetries to reduce the number of k-points.
c Frank Freimuth
c*****************************************************************
USE m_constants
implicit none
integer, intent(in) :: mhp(3)
integer, intent(in) :: nop
logical,intent(in) :: film
......@@ -48,7 +52,7 @@ c*****************************************************************
bbmat=matmul(bmat,transpose(bmat))
write(6,*) "Apply the symmetries to w90kpts"
write(oUnit,*) "Apply the symmetries to w90kpts"
c**************************************************************
c The array 'mhp' specifies the Monkhorst-Pack mesh.
c**************************************************************
......@@ -204,9 +208,9 @@ c close(117)
if (mapk(ikpt)==0)then
sumweights=sumweights+weight(ikpt)
write(6,*)"ikpt=",ikpt
write(6,*)"irreducible"
write(6,fmt='(a10,3f9.6)')"internal: ",kpoint(:)/scale
write(oUnit,*)"ikpt=",ikpt
write(oUnit,*)"irreducible"
write(oUnit,fmt='(a10,3f9.6)')"internal: ",kpoint(:)/scale
if (film.and..not.l_onedimens)then
write(119,'(3f10.5)')kpoint(1:2),weight(ikpt)
......@@ -215,34 +219,34 @@ c close(117)
endif
elseif(mapkoper(ikpt).gt.0)then
c write(6,*)"ikpt=",ikpt
c write(6,*)"reducible"
c write(6,*)"map=",mapk(ikpt)
c write(oUnit,*)"ikpt=",ikpt
c write(oUnit,*)"reducible"
c write(oUnit,*)"map=",mapk(ikpt)
c brot(:)=0.0
c do k=1,3
c brot(:)=brot(:)+
c + mrot(k,:,mapkoper(ikpt))*kpoints(k,mapk(ikpt))
c enddo
c write(6,'(a19,3f9.6)')"rotated internal: ",brot(:)/scale
c write(oUnit,'(a19,3f9.6)')"rotated internal: ",brot(:)/scale
elseif(mapkoper(ikpt).lt.0)then
c write(6,*)"ikpt=",ikpt
c write(6,*)"reducible"
c write(6,*)"map=",mapk(ikpt)
c write(oUnit,*)"ikpt=",ikpt
c write(oUnit,*)"reducible"
c write(oUnit,*)"map=",mapk(ikpt)
c brot(:)=0.0
c do k=1,3
c brot(:)=brot(:)-
c + mrot(k,:,-mapkoper(ikpt))*kpoints(k,mapk(ikpt))
c enddo
c write(6,'(a19,3f9.6)')"rotated internal: ",brot(:)/scale
c write(oUnit,'(a19,3f9.6)')"rotated internal: ",brot(:)/scale
endif
enddo !k1
enddo !k2
enddo !k3
close(119)
write(6,*)"reduznumk=",reduznumk
write(6,*)"nkpts=",nkpts
write(6,*)"sumweights=",sumweights
write(oUnit,*)"reduznumk=",reduznumk
write(oUnit,*)"nkpts=",nkpts
write(oUnit,*)"sumweights=",sumweights
IF(sumweights/=nkpts) CALL juDFT_error
+ ("sum of weights differs from nkpts",calledby
......
......@@ -22,10 +22,13 @@ c****************************************
> tau,
x bkpt,k1,k2,k3,
x zMat,nsfactor)
USE m_types
use m_constants
use m_inv3
use m_constants,only:pimach
implicit none
integer,intent(in) :: natd
integer,intent(in) :: nlod
integer,intent(in) :: llod
......@@ -162,8 +165,8 @@ c print*,testmat
k2(:,:) = k2(:,:) + shiftkpt(2)
k3(:,:) = k3(:,:) + shiftkpt(3)
write(6,*)"in wann_kptsrotate:"
write(6,*) "bkpt=",bkpt
write(oUnit,*)"in wann_kptsrotate:"
write(oUnit,*) "bkpt=",bkpt
arg = tpi*(
+ bkpt(1)*shiftnonsymm(1)+
......
......@@ -22,6 +22,8 @@ c******************************************************************
> invtab,mrot,l_onedimens,tau,
< pair_to_do,maptopair,kdiff,l_q,param_file)
USE m_constants
implicit none
TYPE(t_input),INTENT(IN) :: input
......@@ -207,14 +209,15 @@ c endif !gb=0
10 continue ! end of cycle by the k-points
if(l_p0)then
write(6,*)"pairs to calculate: ",num_pair
write(6,*)"maps by conjugation: ",num_conj
write(6,*)"maps by rotation:", num_rot
write(6,*)"num_pair+num_rot+num_conj:",num_pair+num_conj+num_rot
write(oUnit,*)"pairs to calculate: ",num_pair
write(oUnit,*)"maps by conjugation: ",num_conj
write(oUnit,*)"maps by rotation:", num_rot
write(oUnit,*)"num_pair+num_rot+num_conj:",
+ num_pair+num_conj+num_rot
if(.not.l_q) then
write(6,*)"fullnkpts*nntot:", fullnkpts*nntot
write(oUnit,*)"fullnkpts*nntot:", fullnkpts*nntot
else
write(6,*)"fullnqpts*nntot:", fullnkpts*nntot
write(oUnit,*)"fullnqpts*nntot:", fullnkpts*nntot
endif
endif !l_p0
......@@ -288,10 +291,13 @@ c*****************************************************************
enddo
enddo
end subroutine
SUBROUTINE close_pt(
> nops,mrot,
< mtable)
USE m_constants
IMPLICIT NONE
INTEGER, INTENT (IN) :: nops,mrot(3,3,nops)
......@@ -314,7 +320,7 @@ c*****************************************************************
IF ( map(i) .eq. 0 ) THEN
map(i) = k
ELSE
WRITE (6,'(" Symmetry error : multiple ops")')
WRITE (oUnit,'(" Symmetry error : multiple ops")')
CALL juDFT_error("close_pt: Multiple ops (Bravais)"
+ ,calledby ="wann_mmnk_symm")
ENDIF
......@@ -322,8 +328,8 @@ c*****************************************************************
ENDDO
IF (map(i).eq.0) THEN
WRITE (6,'(" Group not closed (Bravais lattice)")')
WRITE (6,'(" operation j=",i2," map=",12i4,:/,
WRITE (oUnit,'(" Group not closed (Bravais lattice)")')
WRITE (oUnit,'(" operation j=",i2," map=",12i4,:/,
& (21x,12i4))') j, map(1:nops)
CALL juDFT_error("close_pt: Not closed",calledby
+ ="wann_mmnk_symm")
......
......@@ -20,7 +20,7 @@ c*************************************************
> l_soc,band_min,band_max,