Commit db1020c4 authored by Gregor Michalicek's avatar Gregor Michalicek

More adjustments to the Wannier code to integrate it into the new fleur

parent 5fc9552e
......@@ -15,15 +15,13 @@ c m.t. wavefunctions for each band and atom. c.l. fu
c ************************************************************
#include "cpp_double.h"
USE m_cotra, ONLY : cotra3
USE m_constants, ONLY : pimach
USE m_dotir, ONLY : dotirp
USE m_constants
USE m_setabc1locdn
USE m_sphbes
USE m_dsphbs
USE m_abclocdn
USE m_ylm
USE m_od_types, ONLY : od_inp, od_sym
USE m_types
IMPLICIT NONE
C ..
......@@ -186,7 +184,8 @@ c---> (jspin counts the local spin directions inside each MT)
fk(2) = bkpt(2) + k2(k,jspin) + qss2
fk(3) = bkpt(3) + k3(k,jspin) + qss3
ENDIF ! (l_ss)
s = dotirp(fk,fk,bbmat)
s = dot_product(fk,matmul(bbmat,fk))
! s = dotirp(fk,fk,bbmat)
s = sqrt(s)
r1 = rmt(n)*s
CALL sphbes(
......@@ -227,7 +226,8 @@ c ----> construct a and b coefficients
END IF
ENDDO
ENDDO
CALL cotra3(fkr,fkp,bmat)
fkp=MATMUL(fkr,bmat)
! CALL cotra3(fkr,fkp,bmat)
c ----> generate spherical harmonics
CALL ylm4(
> lmax(n),fkp,
......
module m_wann_uHu_commat
USE m_fleurenv
USE m_juDFT
implicit none
contains
......@@ -53,8 +53,8 @@
write(*,*)'dimension:',arr_len
write(*,*)'x?,y?,z? :',l_dim(1:3)
if(arr_len.le.3) call fleur_err("dimension<4",
> calledby='wann_gwf_commat')
if(arr_len.le.3) call juDFT_error("dimension<4",
> calledby='wann_gwf_commat')
nkqpts=nkpts*nqpts
tpi = 2.0*pimach()
......@@ -87,7 +87,7 @@
inquire(file='proj',exist=l_proj)
if(.not.l_proj) write(*,*)'missing: proj'
if(l_miss.or.(.not.l_exist).or.(.not.l_proj))
> call fleur_err("missing file(s) for HDWFs")
> call juDFT_error("missing file(s) for HDWFs")
! get number of bands and wfs from proj
open(405,file='proj',status='old')
......@@ -240,8 +240,8 @@ c write(*,*)'<k|q>',k,bpt(nk,k),q,bpt_q(nq2,q)
if(q.eq.bpt_q(nq2,q)) stop 'q.eq.bpt_q'
countkq=countkq+1
else
call fleur_err("problem k,q neighbors",
> calledby="wann_uHu_commat")
call juDFT_error("problem k,q neighbors",
> calledby="wann_uHu_commat")
endif
endif
enddo
......@@ -406,8 +406,8 @@ c > param_vec(:,qb2)+gb_q(:,nq2,q),latt_const_q,ovaux)
enddo
else
call fleur_err("problem k,q neighbors",
> calledby="wann_uHu_commat")
call juDFT_error("problem k,q neighbors",
> calledby="wann_uHu_commat")
endif
endif
enddo
......
module m_wann_uHu_commat
USE m_fleurenv
USE m_juDFT
implicit none
contains
......@@ -50,8 +50,8 @@
write(*,*)'dimension:',arr_len
write(*,*)'x?,y?,z? :',l_dim(1:3)
if(arr_len.le.3) call fleur_err("dimension<4",
> calledby='wann_gwf_commat')
if(arr_len.le.3) call juDFT_error("dimension<4",
> calledby='wann_gwf_commat')
nkqpts=nkpts*nqpts
tpi = 2.0*pimach()
......@@ -80,7 +80,7 @@
inquire(file='proj',exist=l_proj)
if(.not.l_proj) write(*,*)'missing: proj'
if(l_miss.or.(.not.l_exist).or.(.not.l_proj))
> call fleur_err("missing file(s) for HDWFs")
> call juDFT_error("missing file(s) for HDWFs")
! get number of bands and wfs from proj
open(405,file='proj',status='old')
......@@ -193,8 +193,8 @@ c write(*,'(a,8i8)')'<q|k>',q,k,qb,qb2,kb,kb2,nq,nk2
uHu(:,:,nn2,nn,ikqpt)=cmplx(0.,0.)
c write(*,'(a,8i8)')'<k|q>',q,k,qb,qb2,kb,kb2,nk,nq2
else
call fleur_err("problem k,q neighbors",
> calledby="wann_uHu_commat")
call juDFT_error("problem k,q neighbors",
> calledby="wann_uHu_commat")
endif
endif
enddo
......@@ -331,8 +331,8 @@ c > real(mq(j,i,nq)),aimag(mq(j,i,nq))
enddo
else
call fleur_err("problem k,q neighbors",
> calledby="wann_uHu_commat")
call juDFT_error("problem k,q neighbors",
> calledby="wann_uHu_commat")
endif
endif
enddo
......
......@@ -11,6 +11,7 @@ c*****************************************c
MODULE m_wann_uHu_soc_tlo
CONTAINS
SUBROUTINE wann_uHu_soc_tlo(
> input,atoms,
> ntypd,jmtd,lmaxd,jspd,
> ntype,dx,rmsh,jri,lmax,natd,
> lmd,irank,nlod,llod,
......@@ -23,17 +24,20 @@ c*****************************************c
USE m_intgr, ONLY : intgr0
USE m_sphbes
USE m_dotir
USE m_ylm
USE m_cotra
USE m_gaunt, ONLY: gaunt1
USE m_sointg
USE m_constants, ONLY : c_light
USE m_constants
USE m_types
IMPLICIT NONE
C .. Intrinsic Functions ..
INTRINSIC abs,cmplx,max,mod
TYPE(t_input),INTENT(IN) :: input
TYPE(t_atoms),INTENT(IN) :: atoms
C ..
C .. Scalar Arguments ..
INTEGER, INTENT (IN) :: ntypd,jmtd,lmaxd,jspd,jspins
......@@ -148,9 +152,7 @@ c****************************************************c
> +ello(lop,ntyp,1) + ello(lop,ntyp,jspins) )/4.
END IF
CALL sointg(
> e,vr(:,ntyp,:),v0,r0,dx(ntyp),jri(ntyp),jmtd,c,jspins,
< vso)
CALL sointg(ntyp,e,vr(:,ntyp,:),v0,atoms,input,vso)
! spin-averaging
if(l_spav)then
......@@ -210,9 +212,7 @@ c****************************************************c
> +ello(lop,ntyp,1) + ello(lop,ntyp,jspins) )/4.
END IF
CALL sointg(
> e,vr(:,ntyp,:),v0,r0,dx(ntyp),jri(ntyp),jmtd,c,jspins,
< vso)
CALL sointg(ntyp,e,vr(:,ntyp,:),v0,atoms,input,vso)
! spin-averaging
if(l_spav)then
......
......@@ -26,10 +26,9 @@ c****************************************c
> tuu,tud,tdu,tdd,tuulo,tulou,tdulo,tulod,tuloulo,
> kdiff,kdiff2,nntot,nntot2,uHu)
#include "cpp_double.h"
USE m_fleurenv
USE m_juDFT
use m_constants, only : pimach
use m_matmul , only : matmul3,matmul3r
use m_cotra
implicit none
......@@ -98,7 +97,7 @@ C .. intrinsic functions ..
do nene=1,nntot
if(all(abs(bpt(:)-kdiff(:,nene)).lt.1e-4)) exit
enddo
IF(nene==nntot+1) CALL fleur_err
IF(nene==nntot+1) CALL juDFT_error
+ ("cannot find matching nearest neighbor k+b1",calledby
+ ="wann_uHu_sph")
......@@ -107,7 +106,7 @@ C .. intrinsic functions ..
do nene2=1,nntot2
if(all(abs(bpt2(:)-kdiff2(:,nene2)).lt.1e-4)) exit
enddo
IF(nene2==nntot2+1) CALL fleur_err
IF(nene2==nntot2+1) CALL juDFT_error
+ ("cannot find matching nearest neighbor k+b2",calledby
+ ="wann_uHu_sph")
......
......@@ -7,7 +7,7 @@ c***************************************c
c J.-P. Hanke, Dec. 2015 c
c***************************************c
MODULE m_wann_uHu_vac
USE m_fleurenv
USE m_juDFT
CONTAINS
SUBROUTINE wann_uHu_vac(
......@@ -20,7 +20,6 @@ c***************************************c
#include "cpp_double.h"
use m_constants, only : pimach
use m_intgr, only : intgz0
use m_dotir
use m_d2fdz2cmplx
implicit none
......@@ -156,7 +155,7 @@ c write(*,*)'vz0',ivac,vz0(ivac),vz0_b(ivac)
write(*,*)n2,nv2d,'jspin',jspin
endif
IF (n2>nv2d) CALL fleur_err("wannier uHu vac",calledby
IF (n2>nv2d) CALL juDFT_error("wannier uHu vac",calledby
+ ="wann_uHu_vac")
kvac1(n2) = k1(k,jspin)
......@@ -181,7 +180,7 @@ c.. and for the b-point
write(*,*)n2_b,nv2d,'jspin_b',jspin_b
endif
IF (n2_b>nv2d) CALL fleur_err("wannier uHu vac",calledby
IF (n2_b>nv2d) CALL juDFT_error("wannier uHu vac",calledby
+ ="wann_uHu_vac")
kvac1_b(n2_b) = k1_b(k,jspin_b)
......@@ -236,7 +235,8 @@ c.. the body of the routine
v(1) = bkpt(1) + kvac1(ik) + qss1
v(2) = bkpt(2) + kvac2(ik) + qss2
v(3) = 0.
ev = evacp - 0.5*dotirp(v,v,bbmat)
ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
! ev = evacp - 0.5*dotirp(v,v,bbmat)
call vacuz(ev,vz(1,ivac,jspin2),vz0(ivac),nmz,delz,t(ik),
+ dt(ik),u(1,ik))
call vacudz(ev,vz(1,ivac,jspin2),vz0(ivac),nmz,delz,te(ik),
......@@ -289,7 +289,8 @@ c...now for the k+b point
v(1) = bkpt_b(1) + kvac1_b(ik) + qss1
v(2) = bkpt_b(2) + kvac2_b(ik) + qss2
v(3) = 0.
ev = evacp - 0.5*dotirp(v,v,bbmat)
ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
! ev = evacp - 0.5*dotirp(v,v,bbmat)
call vacuz(ev,vz(1,ivac,jspin2_b),vz0_b(ivac),nmz,delz,
+ t_b(ik),dt_b(ik),u_b(1,ik))
call vacudz(ev,vz(1,ivac,jspin2_b),vz0_b(ivac),nmz,delz,
......@@ -354,7 +355,7 @@ c i2 = -i2
endif
ind2 = ig2(ind3)
IF (ind2.EQ.0) CALL fleur_err("error in map2 for 2-d stars",
IF (ind2.EQ.0) CALL juDFT_error("error in map2 for 2-d stars",
> calledby="wann_uHu_vac")
ind2 = ind2 - 1
IF (ind2.NE.0) THEN ! warping components G=/=0
......@@ -421,13 +422,15 @@ c i2 = -i2
! determine |G+k+b1|^2 and |G'+k+b2|^2
s(1) = bkpt(1) + kvac1(l)
s(2) = bkpt(2) + kvac2(l)
s(3) = 0.0
rk = dotirp(s,s,bbmat)
s(3) = 0.0
rk = dot_product(s,matmul(bbmat,s))
! rk = dotirp(s,s,bbmat)
s(1) = bkpt_b(1) + kvac1_b(lp)
s(2) = bkpt_b(2) + kvac2_b(lp)
s(3) = 0.0
rk_b = dotirp(s,s,bbmat)
rk_b = dot_product(s,matmul(bbmat,s))
! rk_b = dotirp(s,s,bbmat)
c ! no kinetic energy if jspin =/= jspin_b
c if(jspin.ne.jspin_b) then
......
......@@ -9,7 +9,6 @@ c**************************c
< uHu_in,nkpt_loc,counts,displs,nnodes,
> l_unformatted,l_symcc,l_check)
use m_types
use m_fleurenv
use m_wann_uHu_symcheck
implicit none
......
module m_wann_get_qpts
USE m_fleurenv
USE m_juDFT
contains
subroutine wann_get_qpts(
> l_bzsym,film,l_onedimens,l_readqpts,
......@@ -23,33 +23,35 @@ c********************************************************
if(l_bzsym)then
inquire(file='w90qpts',exist=l_file)
IF(.NOT.l_file) CALL fleur_err("where is w90qpts?",calledby
IF(.NOT.l_file) CALL juDFT_error("where is w90qpts?",calledby
+ ="wann_get_qpts")
open(987,file='w90qpts',status='old',form='formatted')
read(987,*)nqpts, scale
write(6,*)"wann_get_qpts: nqpts=",nqpts
if(l_readqpts)then
IF(SIZE(qpoints,1)/=3) CALL fleur_err("wann_get_qpts: 1"
+ ,calledby ="wann_get_qpts")
IF(SIZE(qpoints,2)/=nqpts) CALL fleur_err("wann_get_qpts: 2"
IF(SIZE(qpoints,1)/=3) CALL juDFT_error("wann_get_qpts: 1"
+ ,calledby ="wann_get_qpts")
IF(SIZE(qpoints,2)/=nqpts)
+ CALL juDFT_error("wann_get_qpts: 2",
+ calledby ="wann_get_qpts")
do iter=1,nqpts
read(987,*)qpoints(:,iter)
enddo
endif
else
inquire(file=param_file,exist=l_file)
IF(.NOT.l_file) CALL fleur_err(
IF(.NOT.l_file) CALL juDFT_error(
> "where is "//trim(param_file)//"?",calledby
+ ="wann_get_qpts")
open(987,file=param_file,status='old',form='formatted')
read(987,*)nqpts,scale
write(6,*)"wann_get_qpts: nqpts=",nqpts
if(l_readqpts)then
IF(SIZE(qpoints,1)/=3) CALL fleur_err("wann_get_qpts: 1"
+ ,calledby ="wann_get_qpts")
IF(SIZE(qpoints,2)/=nqpts) CALL fleur_err("wann_get_qpts: 2"
IF(SIZE(qpoints,1)/=3) CALL juDFT_error("wann_get_qpts: 1"
+ ,calledby ="wann_get_qpts")
IF(SIZE(qpoints,2)/=nqpts)
+ CALL juDFT_error("wann_get_qpts: 2",
+ calledby ="wann_get_qpts")
do iter=1,nqpts
read(987,*)qpoints(:,iter)
enddo
......
......@@ -17,6 +17,7 @@ c***************************************************************
CONTAINS
SUBROUTINE wann_mmk0_od_vac(
> DIMENSION, oneD, vacuum, stars, cell,
> l_noco,nlotot,
> z1,nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n2d,n3d,
> ig,nmzxy,nmz,delz,ig2,n2d_1,
......@@ -30,8 +31,16 @@ c***************************************************************
use m_od_abvac
use m_cylbes
use m_dcylbs
USE m_types
implicit none
TYPE(t_dimension),INTENT(IN) :: DIMENSION
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
c .. scalar Arguments..
logical, intent (in) :: l_noco
integer, intent (in) :: nlotot
......@@ -119,10 +128,9 @@ c ..intrinsic functions..
qssbti(3,2) = + qss(3)/2.
DO ispin = 1,1 ! jspins
CALL od_abvac(
> z1,nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n2d,n3d,
> ig,odi%ig,tpi,qssbti(3,jspin),
> nmzxy,nmz,delz,ig2,odi%n2d,
> bbmat,wronk,evac,bkpt,odi%M,odi%mb,
> cell,vacuum,DIMENSION,stars,oneD,
> qssbti(3,jspin),odi%n2d,
> wronk,evac,bkpt,odi%M,odi%mb,
> vz,kvac3(1),nv2,
< uz(1,-vM),duz(1,-vM),u(1,1,-vM),udz(1,-vM),
< dudz(1,-vM),ddnv(1,-vM),ud(1,1,-vM))
......
MODULE m_wann_nocoplot
USE m_fleurenv
USE m_juDFT
c ++++++++++++++++++++++++++++++++++++++++++++++++
c + If noco=t, transform wave functions within +
c + mu-th MT back to the global frame to plot +
......@@ -7,19 +7,20 @@ c + to plot WFs in a meaningful way. +
c + +
c ++++++++++++++++++++++++++++++++++++++++++++++++
CONTAINS
SUBROUTINE wann_nocoplot(slice,nnne,amat,bmat
SUBROUTINE wann_nocoplot(atoms,slice,nnne,amat,bmat
> ,nkpts,odi,film,natd,ntypd,jmtd
> ,ntype,neq,pos,jri,rmsh,alph,beta
> ,nqpts,qss,z1,zatom)
! *****************************************************
USE m_od_types, ONLY : od_inp, od_sym
USE m_dotir
USE m_types
USE m_xsf_io
USE m_cotra, ONLY : cotra1
USE m_wann_gwf_tools, ONLY : get_index_q
USE m_constants, ONLY : pimach
USE m_constants
IMPLICIT NONE
TYPE(t_atoms),INTENT(IN):: atoms
LOGICAL, INTENT(IN) :: slice,film
INTEGER, INTENT(IN) :: nnne,nkpts,ntypd,jmtd,natd
INTEGER, INTENT(IN) :: ntype,neq(ntypd)
......@@ -59,7 +60,7 @@ c ++++++++++++++++++++++++++++++++++++++++++++++++
INQUIRE(file ="plot_inp",exist=twodim)
IF(.NOT.twodim) THEN
CALL fleur_err("Need the plot_inp from RNK generation!"
CALL juDFT_error("Need the plot_inp from RNK generation!"
> ,calledby="wann_nocoplot")
ENDIF
......@@ -68,7 +69,7 @@ c ++++++++++++++++++++++++++++++++++++++++++++++++
READ(18,'(i2,5x,l1)') nplot,xsf
IF (nplot.ge.2)
& CALL fleur_err
& CALL juDFT_error
+ ("plots one by one, please, this is not charge density"
+ ,calledby="wann_nocoplot")
......@@ -78,10 +79,10 @@ c ++++++++++++++++++++++++++++++++++++++++++++++++
zero = (/0.,0.,0./);filename="default"
READ(18,plot)
IF (twodim.AND.ANY(grid(1:2)<1))
+ CALL fleur_err("Illegal grid size in plot",calledby
+ CALL juDFT_error("Illegal grid size in plot",calledby
+ ="wann_nocoplot")
IF (.NOT.twodim.AND.ANY(grid<1))
+ CALL fleur_err("Illegal grid size in plot",calledby
+ CALL juDFT_error("Illegal grid size in plot",calledby
+ ="wann_nocoplot")
IF (twodim) grid(3) = 1
!calculate cartesian coordinates if needed
......@@ -130,17 +131,11 @@ c ++++++++++++++++++++++++++++++++++++++++++++++++
write (name3,24) ikpt,nbn,jspin
24 format (i5.5,'.',i3.3,'.absv.',i1,'.xsf')
OPEN(600+jspin,file=name1)
CALL xsf_WRITE_atoms(
> 600+jspin,film,odi%d1,amat,neq(:ntype),
> zatom(:ntype),pos)
CALL xsf_WRITE_atoms(600+jspin,atoms,film,odi%d1,amat)
OPEN(602+jspin,file=name2)
CALL xsf_WRITE_atoms(
> 602+jspin,film,odi%d1,amat,neq(:ntype),
> zatom(:ntype),pos)
CALL xsf_WRITE_atoms(602+jspin,atoms,film,odi%d1,amat)
OPEN(604+jspin,file=name3)
CALL xsf_WRITE_atoms(
> 604+jspin,film,odi%d1,amat,neq(:ntype),
> zatom(:ntype),pos)
CALL xsf_WRITE_atoms(604+jspin,atoms,film,odi%d1,amat)
CALL xsf_WRITE_header(600+jspin,twodim,filename,(vec1),
& (vec2),(vec3),zero
$ ,grid)
......@@ -195,7 +190,8 @@ c ++++++++++++++++++++++++++++++++++++++++++++++++
! transform wfs to global frame
DO iintsp=1,2
CALL cotra1(pt(1),rcc2,bmat)
pt = matmul(bmat,rcc2)/tpi_const
! CALL cotra1(pt(1),rcc2,bmat)
arg = -pi*real(2*iintsp-3)*( qss(1,iqpt)*rcc2(1)
> +qss(2,iqpt)*rcc2(2)
> +qss(3,iqpt)*rcc2(3))
......
......@@ -158,7 +158,7 @@ c... rotate into representative
60 CONTINUE
70 CONTINUE
c... switch back to cartesian units
x=matmul(amat,p)
x=matmul(amat,p)/tpi_const
END IF
DO 80 j = jri(n),2,-1
IF (sx.GE.rmsh(j,n)) GO TO 90
......
......@@ -41,7 +41,7 @@ c****************************************************************************
use m_od_types, only : od_inp, od_sym
use m_loddop
use m_constants, only : pimach
use m_wann_mmk0_od_vac
! use m_wann_mmk0_od_vac
use m_wann_mmk0_vac
use m_wann_mmk0_updown_sph
use m_wann_abinv
......
......@@ -8,7 +8,6 @@ c**********************************************************
c**********************************************************
use m_types
use m_fleurenv
implicit none
logical, intent(in) :: l_p0,l_unformatted
integer, intent(in) :: fullnkpts
......
Markdown is supported
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