Commit 0cf0bbf4 authored by Daniel Wortmann's avatar Daniel Wortmann

Bugfix in noco-vacuum part. Should hopefully fix #233

parent b066adc9
......@@ -29,7 +29,7 @@ CONTAINS
TYPE(t_lapw),INTENT(IN) :: lapw
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_potden),INTENT(IN) :: v
CLASS(t_mat),INTENT(INOUT) :: hmat(:,:),smat(:,:)
CLASS(t_mat),INTENT(INOUT) :: hmat(:,:),smat(:,:)
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jsp
......@@ -40,14 +40,14 @@ CONTAINS
! .. Local Scalars ..
COMPLEX hij,sij,apw_lo,c_1
REAL d2,gz,sign,th,wronk
INTEGER i,i2,ii,jj,ik,j,jk,k,jspin,ipot,ii0,i0
INTEGER i,i2,ii,jj,ik,j,jk,k,jspin,ii0,i0
INTEGER ivac,irec,imz,igvm2,igvm2i,s1,s2
INTEGER jspin1,jspin2,jmax,jsp_start,jsp_end
INTEGER i_start,nc,nc_0
! ..
! .. Local Arrays ..
INTEGER:: nv2(input%jspins)
INTEGER kvac1(DIMENSION%nv2d,input%jspins),kvac2(DIMENSION%nv2d,input%jspins)
INTEGER kvac(2,DIMENSION%nv2d,input%jspins)
INTEGER map2(DIMENSION%nvd,input%jspins)
COMPLEX tddv(DIMENSION%nv2d,DIMENSION%nv2d),tduv(DIMENSION%nv2d,DIMENSION%nv2d)
COMPLEX tudv(DIMENSION%nv2d,DIMENSION%nv2d),tuuv(DIMENSION%nv2d,DIMENSION%nv2d)
......@@ -66,16 +66,14 @@ CONTAINS
nv2(jspin) = 0
k_loop:DO k = 1,lapw%nv(jspin)
DO j = 1,nv2(jspin)
IF (lapw%k1(k,jspin).EQ.kvac1(j,jspin)&
.AND. lapw%k2(k,jspin).EQ.kvac2(j,jspin)) THEN
IF (all(lapw%gvec(1:2,k,jspin)==kvac(1:2,j,jspin))) THEN
map2(k,jspin) = j
CYCLE k_loop
END IF
ENDDO
nv2(jspin) = nv2(jspin) + 1
IF (nv2(jspin)>DIMENSION%nv2d) CALL juDFT_error("hsvac:dimension%nv2d",calledby ="hsvac")
kvac1(nv2(jspin),jspin) = lapw%k1(k,jspin)
kvac2(nv2(jspin),jspin) = lapw%k2(k,jspin)
kvac(1:2,nv2(jspin),jspin) = lapw%gvec(1:2,k,jspin)
map2(k,jspin) = nv2(jspin)
ENDDO k_loop
ENDDO
......@@ -86,20 +84,16 @@ CONTAINS
s1=MIN(SIZE(hmat,1),jspin1) !in colinear case s1=1
DO jspin2=MERGE(1,jsp,noco%l_noco),MERGE(2,jsp,noco%l_noco) !loop over global spin
s2=MIN(SIZE(hmat,1),jspin2) !in colinear case s2=1
ipot=3
IF (jspin1==jspin2) ipot=jspin1
!---> get the wavefunctions and set up the tuuv, etc matrices
jspin=jsp
CALL vacfun(&
vacuum,DIMENSION,stars,&
jsp,input,noco,jspin1,jspin2,&
sym, cell,ivac,evac(1,1),lapw%bkpt,v%vacxy(:,:,ivac,ipot),v%vacz(:,:,:),kvac1,kvac2,nv2,&
tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
CALL vacfun(&
vacuum,stars,&
input,noco,jspin1,jspin2,&
sym, cell,ivac,evac,lapw%bkpt,v%vacxy,v%vacz,kvac,nv2,&
tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
!
!---> generate a and b coeffficients
!
IF (noco%l_noco) THEN
DO jspin = 1,input%jspins
DO jspin = MIN(jspin1,jspin2),MAX(jspin1,jspin2)
DO k = 1,lapw%nv(jspin)
gz = sign*cell%bmat(3,3)*lapw%k3(k,jspin)
i2 = map2(k,jspin)
......@@ -109,16 +103,6 @@ CONTAINS
b(k,jspin) = c_1 * CMPLX(duz(i2,jspin), gz* uz(i2,jspin) )
ENDDO
ENDDO
ELSE
DO k = 1,lapw%nv(jsp)
gz = sign*cell%bmat(3,3)*lapw%gvec(3,k,jsp)
i2 = map2(k,jsp)
th = gz*cell%z1
c_1 = CMPLX( COS(th), SIN(th) )/ (d2*wronk)
a(k,jsp) = - c_1 * CMPLX(dudz(i2,jsp), gz*udz(i2,jsp) )
b(k,jsp) = c_1 * CMPLX(duz(i2,jsp), gz* uz(i2,jsp) )
ENDDO
ENDIF
!---> update hamiltonian and overlap matrices
IF (jspin1==jspin2) THEN
DO i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size
......@@ -153,30 +137,25 @@ CONTAINS
!Diagonal term of Overlapp matrix, Hamiltonian later
sij = CONJG(a(i,jspin2))*a(i,jspin2) + CONJG(b(i,jspin2))*b(i,jspin2)*ddnv(ik,jspin2)
IF (hmat(1,1)%l_real) THEN
smat(s1,s2)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij)
smat(s2,s1)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij)
ELSE
smat(s1,s2)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij
smat(s2,s1)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij
ENDIF
ENDDO
ENDIF
!---> hamiltonian update
DO i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size
DO i = mpi%n_rank+1,lapw%nv(jspin1),mpi%n_size
i0=(i-1)/mpi%n_size+1 !local column index
ik = map2(i,jspin2)
DO j = 1,MIN(i,lapw%nv(jspin1))
jk = map2(j,jspin1)
IF (jspin2>jspin1) THEN
hij = CONJG(CONJG(a(j,jspin1))* (tuuv(jk,ik)*a(i,jspin2) +tudv(jk,ik)*b(i,jspin2))&
+ CONJG(b(j,jspin1))* (tddv(jk,ik)*b(i,jspin2) +tduv(jk,ik)*a(i,jspin2)))
ELSE
hij = CONJG(a(i,jspin2))* (tuuv(ik,jk)*a(j,jspin1) +tudv(ik,jk)*b(j,jspin1))&
+ CONJG(b(i,jspin2))* (tddv(ik,jk)*b(j,jspin1) +tduv(ik,jk)*a(j,jspin1))
ENDIF
ik = map2(i,jspin1)
DO j = 1,MERGE(i,lapw%nv(jspin2),jspin1==jspin2)
jk = map2(j,jspin2)
hij = CONJG(a(i,jspin1))* (tuuv(ik,jk)*a(j,jspin2) +tudv(ik,jk)*b(j,jspin2))&
+ CONJG(b(i,jspin1))* (tddv(ik,jk)*b(j,jspin2) +tduv(ik,jk)*a(j,jspin2))
IF (hmat(1,1)%l_real) THEN
hmat(s1,s2)%data_r(j,i0) = hmat(s1,s2)%data_r(j,i0) + REAL(hij)
hmat(s2,s1)%data_r(j,i0) = hmat(s2,s1)%data_r(j,i0) + REAL(hij)
ELSE
hmat(s1,s2)%data_c(j,i0) = hmat(s1,s2)%data_c(j,i0) + hij
hmat(s2,s1)%data_c(j,i0) = hmat(s2,s1)%data_c(j,i0) + hij
ENDIF
ENDDO
ENDDO
......
......@@ -2,8 +2,8 @@ MODULE m_vacfun
use m_juDFT
CONTAINS
SUBROUTINE vacfun(&
vacuum,DIMENSION,stars, jsp,input,noco,jsp1,jsp2,&
sym, cell,ivac,evac,bkpt, vxy,vz,kvac1,kvac2,nv2,&
vacuum,stars,input,noco,jspin1,jspin2,&
sym, cell,ivac,evac,bkpt, vxy,vz,kvac,nv2,&
tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
!*********************************************************************
! determines the necessary values and derivatives on the vacuum
......@@ -19,7 +19,6 @@ CONTAINS
USE m_types
IMPLICIT NONE
TYPE(t_dimension),INTENT(IN) :: dimension
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
......@@ -28,31 +27,35 @@ CONTAINS
TYPE(t_cell),INTENT(IN) :: cell
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jsp ,ivac,jsp1,jsp2
INTEGER, INTENT (IN) :: ivac,jspin1,jspin2
REAL, INTENT (OUT) :: wronk
! ..
! .. Array Arguments ..
INTEGER, INTENT (IN) :: nv2(input%jspins)
INTEGER, INTENT (IN) :: kvac1(dimension%nv2d,input%jspins),kvac2(dimension%nv2d,input%jspins)
COMPLEX, INTENT (IN) :: vxy(vacuum%nmzxyd,stars%ng2-1)
COMPLEX, INTENT (OUT):: tddv(dimension%nv2d,dimension%nv2d),tduv(dimension%nv2d,dimension%nv2d)
COMPLEX, INTENT (OUT):: tudv(dimension%nv2d,dimension%nv2d),tuuv(dimension%nv2d,dimension%nv2d)
REAL, INTENT (IN) :: vz(vacuum%nmzd,2,4) ,evac(2,input%jspins)
INTEGER, INTENT (IN) :: nv2(:)!(input%jspins)
INTEGER, INTENT (IN) :: kvac(:,:,:)!(2,dimension%nv2d,input%jspins)
COMPLEX, INTENT (IN) :: vxy(:,:,:,:) !(vacuum%nmzxyd,stars%ng2-1,nvac,:)
COMPLEX, INTENT (OUT):: tddv(:,:),tduv(:,:)!(dimension%nv2d,dimension%nv2d)
COMPLEX, INTENT (OUT):: tudv(:,:),tuuv(:,:)!(dimension%nv2d,dimension%nv2d)
REAL,ALLOCATABLE,INTENT (IN) :: vz(:,:,:) !(vacuum%nmzd,2,4) ,
REAL, INTENT (IN) :: evac(:,:)!(2,input%jspins)
REAL, INTENT (IN) :: bkpt(3)
REAL, INTENT (OUT):: udz(dimension%nv2d,input%jspins),uz(dimension%nv2d,input%jspins)
REAL, INTENT (OUT):: dudz(dimension%nv2d,input%jspins),duz(dimension%nv2d,input%jspins)
REAL, INTENT (OUT):: ddnv(dimension%nv2d,input%jspins)
REAL, INTENT (OUT):: udz(:,:),uz(:,:)!(dimension%nv2d,input%jspins)
REAL, INTENT (OUT):: dudz(:,:),duz(:,:)!(dimension%nv2d,input%jspins)
REAL, INTENT (OUT):: ddnv(:,:)!(dimension%nv2d,input%jspins)
! ..
! .. Local Scalars ..
REAL ev,scale,xv,yv,vzero
REAL ev,scale,xv,yv,vzero,fac
COMPLEX phase
INTEGER i,i1,i2,i3,ik,ind2,ind3,jk,np1,jspin
INTEGER i,i1,i2,i3,ik,ind2,ind3,jk,np1,jspin,ipot
LOGICAL tail
! ..
! .. Local Arrays ..
REAL u(vacuum%nmzd,dimension%nv2d,input%jspins),ud(vacuum%nmzd,dimension%nv2d,input%jspins)
REAL u(vacuum%nmzd,size(duz,1),input%jspins),ud(vacuum%nmzd,size(duz,1),input%jspins)
REAL v(3),x(vacuum%nmzd), qssbti(2,2)
! ..
fac=MERGE(1.0,-1.0,jspin1<=jspin2)
ipot=MERGE(jspin1,3,jspin1==jspin2)
tuuv=0.0;tudv=0.0;tddv=0.0;tduv=0.0
udz=0.0;duz=0.0;ddnv=0.0;udz=0.;uz=0.
tail = .true.
......@@ -60,15 +63,12 @@ CONTAINS
!---> wronksian for the schrodinger equation given by an identity
wronk = 2.0
!---> setup the spin-spiral q-vector
qssbti(1,1) = - noco%qss(1)/2
qssbti(2,1) = - noco%qss(2)/2
qssbti(1,2) = + noco%qss(1)/2
qssbti(2,2) = + noco%qss(2)/2
qssbti(1:2,1) = - noco%qss(1:2)/2
qssbti(1:2,2) = + noco%qss(1:2)/2
!---> generate basis functions for each 2-d k+g
DO jspin = 1,input%jspins
DO jspin = MIN(jspin1,jspin2),MAX(jspin1,jspin2)
DO ik = 1,nv2(jspin)
v(1) = bkpt(1) + kvac1(ik,jspin) + qssbti(1,jspin)
v(2) = bkpt(2) + kvac2(ik,jspin) + qssbti(2,jspin)
v(1:2) = bkpt(1:2) + kvac(:,ik,jspin) + qssbti(1:2,jspin)
v(3) = 0.0
ev = evac(ivac,jspin) - 0.5*dot_product(v,matmul(v,cell%bbmat))
vzero = vz(vacuum%nmzd,ivac,jspin)
......@@ -80,19 +80,19 @@ CONTAINS
!---> make sure the solutions satisfy the wronksian
scale = wronk/ (udz(ik,jspin)*duz(ik,jspin)-&
& dudz(ik,jspin)*uz(ik,jspin))
udz(ik,jspin) = scale*udz(ik,jspin)
udz(ik,jspin) = scale*udz(ik,jspin)
dudz(ik,jspin) = scale*dudz(ik,jspin)
ddnv(ik,jspin) = scale*ddnv(ik,jspin)
ud(:,ik,jspin) = scale*ud(:,ik,jspin)
enddo
ENDDO
!---> set up the tuuv, etc. matrices
DO ik = 1,nv2(jsp1)
DO jk = 1,nv2(jsp2)
DO ik = 1,nv2(jspin1)
DO jk = 1,nv2(jspin2)
!---> determine the warping component of the potential
i1 = kvac1(ik,jsp1) - kvac1(jk,jsp2)
i2 = kvac2(ik,jsp1) - kvac2(jk,jsp2)
i1 = kvac(1,ik,jspin1) - kvac(1,jk,jspin2)
i2 = kvac(2,ik,jspin1) - kvac(2,jk,jspin2)
i3 = 0
ind3 = stars%ig(i1,i2,i3)
IF (ind3.EQ.0) CYCLE
......@@ -114,98 +114,96 @@ CONTAINS
!---> tuuv
DO i = 1,vacuum%nmzxy
x(np1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*real(vxy(i,ind2))
x(np1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*REAL(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
DO i = 1,vacuum%nmzxy
x(np1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*aimag(vxy(i,ind2))
x(np1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
tuuv(ik,jk) = phase*cmplx(xv,yv)
!---> tddv
DO i = 1,vacuum%nmzxy
x(np1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*real(vxy(i,ind2))
x(np1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*REAL(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
DO i = 1,vacuum%nmzxy
x(np1-i) =ud(i,ik,jsp1)*ud(i,jk,jsp2)*aimag(vxy(i,ind2))
x(np1-i) =ud(i,ik,jspin1)*ud(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
tddv(ik,jk) = phase*cmplx(xv,yv)
!---> tudv
DO i = 1,vacuum%nmzxy
x(np1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*real(vxy(i,ind2))
x(np1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*real(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
DO i = 1,vacuum%nmzxy
x(np1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*aimag(vxy(i,ind2))
x(np1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
tudv(ik,jk) = phase*cmplx(xv,yv)
!---> tduv
DO i = 1,vacuum%nmzxy
x(np1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*real(vxy(i,ind2))
x(np1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*real(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
DO i = 1,vacuum%nmzxy
x(np1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*aimag(vxy(i,ind2))
x(np1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
tduv(ik,jk) = phase*cmplx(xv,yv)
ELSE
!---> diagonal (film muffin-tin) terms
IF (jsp1==jsp2) THEN
tuuv(ik,ik) = cmplx(evac(ivac,jsp1),0.0)
tddv(ik,ik) = cmplx(evac(ivac,jsp1)*ddnv(ik,jsp1),0.0)
IF (jspin1==jspin2) THEN
tuuv(ik,ik) = cmplx(evac(ivac,jspin1),0.0)
tddv(ik,ik) = cmplx(evac(ivac,jspin1)*ddnv(ik,jspin1),0.0)
tudv(ik,ik) = cmplx(0.5,0.0)
tduv(ik,ik) = cmplx(0.5,0.0)
ELSE
!---> tuuv
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,3)
x(vacuum%nmz+1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*vz(i,ivac,3)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,4)
x(vacuum%nmz+1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*fac*vz(i,ivac,4)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
tuuv(ik,jk) = cmplx(xv,yv)
!---> tddv
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,3)
x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*vz(i,ivac,3)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,4)
x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*fac*vz(i,ivac,4)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
tddv(ik,jk) = cmplx(xv,yv)
!---> tudv
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,3)
x(vacuum%nmz+1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*vz(i,ivac,3)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,4)
x(vacuum%nmz+1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*fac*vz(i,ivac,4)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
tudv(ik,jk) = cmplx(xv,yv)
!---> tduv
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,3)
x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*vz(i,ivac,3)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,4)
x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*fac*vz(i,ivac,4)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
tduv(ik,jk) = cmplx(xv,yv)
......
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