Commit 14d7c1c0 authored by Daniel Wortmann's avatar Daniel Wortmann

Bugfix in SS part of Hamiltonian setup. Fixes #194

parent 0f321e18
......@@ -60,8 +60,8 @@ CONTAINS
ENDIF
DO i=1,nspins
DO j=1,nspins
CALL smat(i,j)%init(l_real,lapw%nv(i)+atoms%nlotot,lapw%nv(i)+atoms%nlotot,mpi%sub_comm,.false.)
CALL hmat(i,j)%init(l_real,lapw%nv(i)+atoms%nlotot,lapw%nv(i)+atoms%nlotot,mpi%sub_comm,.false.)
CALL smat(i,j)%init(l_real,lapw%nv(i)+atoms%nlotot,lapw%nv(j)+atoms%nlotot,mpi%sub_comm,.false.)
CALL hmat(i,j)%init(l_real,lapw%nv(i)+atoms%nlotot,lapw%nv(j)+atoms%nlotot,mpi%sub_comm,.false.)
ENDDO
ENDDO
......
......@@ -46,8 +46,8 @@ CONTAINS
CALL mat(2,2)%free()
!Now collect off-diagonal parts
CALL mat(2,1)%add_transpose(mat(1,2))
CALL mat_final%copy(mat(2,1),1,lapw%nv(1)+atoms%nlotot+1)
CALL mat(1,2)%add_transpose(mat(2,1))
CALL mat_final%copy(mat(1,2),1,lapw%nv(1)+atoms%nlotot+1)
CALL mat(1,2)%free()
CALL mat(2,1)%free()
......
......@@ -42,18 +42,19 @@ CONTAINS
DO i = mpi%n_rank+1,lapw%nv(ispin),mpi%n_size
i0=(i-1)/mpi%n_size+1
!---> loop over (k+g)
DO j = 1,i
ii = lapw%gvec(:,i,jspin) - lapw%gvec(:,j,ispin)
DO j = 1,MIN(i,lapw%nv(jspin))
ii = lapw%gvec(:,i,ispin) - lapw%gvec(:,j,jspin)
IF (ispin==1.AND.jspin==2) THEN
ii=-1*ii
in = stars%ig(ii(1),ii(2),ii(3))
IF (in.EQ.0) CYCLE
th = stars%rgphs(ii(1),ii(2),ii(3))*vpw(in,3)
th = stars%rgphs(ii(1),ii(2),ii(3))*CONJG(vpw(in,3))
ts=0.0
ELSEIF(ispin==2.and.jspin==1) THEN
ii = -1*ii
! ii = -1*ii
in = stars%ig(ii(1),ii(2),ii(3))
IF (in.EQ.0) CYCLE
th = stars%rgphs(ii(1),ii(2),ii(3))*CONJG(vpw(in,3))
th = stars%rgphs(ii(1),ii(2),ii(3))*vpw(in,3)
ts=0.0
ELSE
!--> determine index and phase factor
......
......@@ -40,8 +40,8 @@ CONTAINS
real,allocatable :: gkrot(:,:)
LOGICAL :: l_apw
ALLOCATE(c_ph(lapw%nv(1),MERGE(2,1,noco%l_ss)))
ALLOCATE(gkrot(3,lapw%nv(1)))
ALLOCATE(c_ph(maxval(lapw%nv),MERGE(2,1,noco%l_ss)))
ALLOCATE(gkrot(3,MAXVAL(lapw%nv)))
lmax=MERGE(atoms%lnonsph(n),atoms%lmax(n),l_nonsph)
......
......@@ -62,9 +62,9 @@ CONTAINS
COMPLEX,ALLOCATABLE:: ab(:,:),ab1(:,:),ab2(:,:)
real :: rchi
ALLOCATE(ab(MAXVAL(lapw%nv),2*atoms%lmaxd*(atoms%lmaxd+2)+2),ab1(lapw%nv(iintsp),2*atoms%lmaxd*(atoms%lmaxd+2)+2))
ALLOCATE(ab(MAXVAL(lapw%nv),2*atoms%lmaxd*(atoms%lmaxd+2)+2),ab1(lapw%nv(jintsp),2*atoms%lmaxd*(atoms%lmaxd+2)+2))
IF (iintsp.NE.jintsp) ALLOCATE(ab2(lapw%nv(jintsp),2*atoms%lmaxd*(atoms%lmaxd+2)+2))
IF (iintsp.NE.jintsp) ALLOCATE(ab2(lapw%nv(iintsp),2*atoms%lmaxd*(atoms%lmaxd+2)+2))
IF (hmat%l_real) THEN
IF (ANY(SHAPE(hmat%data_c)/=SHAPE(hmat%data_r))) THEN
......@@ -79,18 +79,18 @@ CONTAINS
IF ((atoms%invsat(na)==0) .OR. (atoms%invsat(na)==1)) THEN
rchi=MERGE(REAL(chi),REAL(chi)*2,(atoms%invsat(na)==0))
CALL hsmt_ab(sym,atoms,noco,isp,iintsp,n,na,cell,lapw,fj,gj,ab,ab_size,.TRUE.)
CALL hsmt_ab(sym,atoms,noco,isp,jintsp,n,na,cell,lapw,fj,gj,ab,ab_size,.TRUE.)
!Calculate Hamiltonian
CALL zgemm("N","N",lapw%nv(iintsp),ab_size,ab_size,CMPLX(1.0,0.0),ab,SIZE(ab,1),td%h_loc(0:,0:,n,isp),SIZE(td%h_loc,1),CMPLX(0.,0.),ab1,SIZE(ab1,1))
CALL zgemm("N","N",lapw%nv(jintsp),ab_size,ab_size,CMPLX(1.0,0.0),ab,SIZE(ab,1),td%h_loc(0:,0:,n,isp),SIZE(td%h_loc,1),CMPLX(0.,0.),ab1,SIZE(ab1,1))
!ab1=MATMUL(ab(:lapw%nv(iintsp),:ab_size),td%h_loc(:ab_size,:ab_size,n,isp))
IF (iintsp==jintsp) THEN
CALL ZHERK("U","N",lapw%nv(iintsp),ab_size,Rchi,CONJG(ab1),SIZE(ab1,1),1.0,hmat%data_c,SIZE(hmat%data_c,1))
ELSE !here the l_ss off-diagonal part starts
!Second set of ab is needed
CALL hsmt_ab(sym,atoms,noco,isp,jintsp,n,na,cell,lapw,fj,gj,ab,ab_size,.TRUE.)
CALL zgemm("N","N",lapw%nv(jintsp),ab_size,ab_size,CMPLX(1.0,0.0),ab,SIZE(ab,1),td%h_loc(:,:,n,isp),SIZE(td%h_loc,1),CMPLX(0.,0.),ab2,SIZE(ab2,1))
CALL hsmt_ab(sym,atoms,noco,isp,iintsp,n,na,cell,lapw,fj,gj,ab,ab_size,.TRUE.)
CALL zgemm("N","N",lapw%nv(iintsp),ab_size,ab_size,CMPLX(1.0,0.0),ab,SIZE(ab,1),td%h_loc(:,:,n,isp),SIZE(td%h_loc,1),CMPLX(0.,0.),ab2,SIZE(ab2,1))
!Multiply for Hamiltonian
CALL zgemm("N","T",lapw%nv(jintsp),lapw%nv(iintsp),ab_size,chi,conjg(ab2),SIZE(ab2,1),ab1,SIZE(ab1,1),CMPLX(1.0,0.0),hmat%data_c,SIZE(hmat%data_c,1))
CALL zgemm("N","T",lapw%nv(iintsp),lapw%nv(jintsp),ab_size,chi,conjg(ab2),SIZE(ab2,1),ab1,SIZE(ab1,1),CMPLX(1.0,0.0),hmat%data_c,SIZE(hmat%data_c,1))
ENDIF
ENDIF
END DO
......@@ -130,7 +130,7 @@ CONTAINS
COMPLEX,ALLOCATABLE:: ab(:,:),ab1(:,:),ab_select(:,:)
real :: rchi
ALLOCATE(ab(MAXVAL(lapw%nv),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2),ab1(lapw%nv(iintsp),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2),ab_select(lapw%num_local_cols(jintsp),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2))
ALLOCATE(ab(MAXVAL(lapw%nv),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2),ab1(lapw%nv(jintsp),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2),ab_select(lapw%num_local_cols(jintsp),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2))
!IF (iintsp.NE.jintsp) ALLOCATE(ab_select1(lapw%num_local_cols(jintsp),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2))
......@@ -147,20 +147,20 @@ CONTAINS
IF ((atoms%invsat(na)==0) .OR. (atoms%invsat(na)==1)) THEN
rchi=MERGE(REAL(chi),REAL(chi)*2,(atoms%invsat(na)==0))
CALL hsmt_ab(sym,atoms,noco,isp,iintsp,n,na,cell,lapw,fj,gj,ab,ab_size,.TRUE.)
CALL hsmt_ab(sym,atoms,noco,isp,jintsp,n,na,cell,lapw,fj,gj,ab,ab_size,.TRUE.)
!Calculate Hamiltonian
CALL zgemm("N","N",lapw%nv(iintsp),ab_size,ab_size,CMPLX(1.0,0.0),ab,SIZE(ab,1),td%h_loc(0:,0:,n,isp),SIZE(td%h_loc,1),CMPLX(0.,0.),ab1,SIZE(ab1,1))
CALL zgemm("N","N",lapw%nv(jintsp),ab_size,ab_size,CMPLX(1.0,0.0),ab,SIZE(ab,1),td%h_loc(0:,0:,n,isp),SIZE(td%h_loc,1),CMPLX(0.,0.),ab1,SIZE(ab1,1))
!Cut out of ab1 only the needed elements here
ab_select=ab1(mpi%n_rank+1:lapw%nv(iintsp):mpi%n_size,:)
ab_select=ab1(mpi%n_rank+1:lapw%nv(jintsp):mpi%n_size,:)
IF (iintsp==jintsp) THEN
CALL zgemm("N","T",lapw%nv(iintsp),lapw%num_local_cols(iintsp),ab_size,CMPLX(rchi,0.0),CONJG(ab1),SIZE(ab1,1),ab_select,lapw%num_local_cols(iintsp),CMPLX(1.,0.0),hmat%data_c,SIZE(hmat%data_c,1))
ELSE
!Second set of ab is needed
CALL hsmt_ab(sym,atoms,noco,isp,jintsp,n,na,cell,lapw,fj,gj,ab,ab_size,.TRUE.)
CALL hsmt_ab(sym,atoms,noco,isp,iintsp,n,na,cell,lapw,fj,gj,ab,ab_size,.TRUE.)
CALL zgemm("N","N",lapw%nv(iintsp),ab_size,ab_size,CMPLX(1.0,0.0),ab,SIZE(ab,1),td%h_loc(:,:,n,isp),SIZE(td%h_loc,1),CMPLX(0.,0.),ab1,SIZE(ab1,1))
!Multiply for Hamiltonian
CALL zgemm("N","t",lapw%nv(iintsp),lapw%num_local_cols(iintsp),ab_size,chi,conjg(ab1),SIZE(ab1,1),ab_select,lapw%num_local_cols(iintsp),CMPLX(1.,0.0),hmat%data_c,SIZE(hmat%data_c,1))
CALL zgemm("N","T",lapw%nv(iintsp),lapw%num_local_cols(jintsp),ab_size,chi,conjg(ab1),SIZE(ab1,1),ab_select,lapw%num_local_cols(jintsp),CMPLX(1.,0.0),hmat%data_c,SIZE(hmat%data_c,1))
ENDIF
ENDIF
END DO
......
......@@ -63,15 +63,15 @@ CONTAINS
ALLOCATE(plegend(MAXVAL(lapw%nv),0:atoms%lmaxd))
plegend=0.0
plegend(:,0)=1.0
qssbti=MERGE(- noco%qss/2,+ noco%qss/2,iintsp.EQ.1)
qssbtj=MERGE(- noco%qss/2,+ noco%qss/2,jintsp.EQ.1)
qssbti=MERGE(- noco%qss/2,+ noco%qss/2,jintsp.EQ.1)
qssbtj=MERGE(- noco%qss/2,+ noco%qss/2,iintsp.EQ.1)
!$OMP DO SCHEDULE(DYNAMIC,1)
DO ki = mpi%n_rank+1, lapw%nv(iintsp), mpi%n_size
DO ki = mpi%n_rank+1, lapw%nv(jintsp), mpi%n_size
kii=(ki-1)/mpi%n_size+1
ski = lapw%gvec(:,ki,iintsp) + qssbti
ski = lapw%gvec(:,ki,jintsp) + qssbti
!---> legendre polynomials
DO kj = 1,ki
plegend(kj,1) = DOT_PRODUCT(lapw%gk(:,kj,jintsp),lapw%gk(:,ki,iintsp))
plegend(kj,1) = DOT_PRODUCT(lapw%gk(:,kj,iintsp),lapw%gk(:,ki,jintsp))
END DO
DO l = 1,atoms%lmax(n) - 1
plegend(:ki,l+1) = fleg1(l)*plegend(:ki,1)*plegend(:ki,l) - fleg2(l)*plegend(:ki,l-1)
......@@ -82,16 +82,16 @@ CONTAINS
tnn = tpi_const*atoms%taual(:,nn)
DO kj = 1,ki
cph(kj) = cph(kj) +&
CMPLX(COS(DOT_PRODUCT(ski-lapw%gvec(:,kj,jintsp)-qssbtj,tnn)),&
SIN(DOT_PRODUCT(lapw%gvec(:,kj,jintsp)+qssbtj-ski,tnn)))
IF (iintsp.NE.jintsp) cph(kj)=CONJG(cph(kj))
CMPLX(COS(DOT_PRODUCT(ski-lapw%gvec(:,kj,iintsp)-qssbtj,tnn)),&
SIN(DOT_PRODUCT(lapw%gvec(:,kj,iintsp)+qssbtj-ski,tnn)))
! IF (iintsp.NE.jintsp) cph(kj)=CONJG(cph(kj))
END DO
END DO
!---> update overlap and l-diagonal hamiltonian matrix
DO l = 0,atoms%lmax(n)
fjkiln = fj(ki,l,iintsp)
gjkiln = gj(ki,l,iintsp)
fjkiln = fj(ki,l,jintsp)
gjkiln = gj(ki,l,jintsp)
!
IF (input%l_useapw) THEN
w1 = 0.5 * ( usdus%uds(l,n,isp)*usdus%dus(l,n,isp) + &
......@@ -108,11 +108,11 @@ CONTAINS
IF (smat%l_real) THEN
DO kj = 1,ki
fct = plegend(kj,l)*fl2p1(l)*&
( fjkiln*fj(kj,l,jintsp) + gjkiln*gj(kj,l,jintsp)*ddnln )
( fjkiln*fj(kj,l,iintsp) + gjkiln*gj(kj,l,iintsp)*ddnln )
smat%data_r(kj,kii)=smat%data_r(kj,kii)+REAL(cph(kj))*fct
hmat%data_r(kj,kii)=hmat%data_r(kj,kii) + REAL(cph(kj)) * &
( fct * elall + plegend(kj,l) * fl2p1bt(l) *&
( fjkiln*gj(kj,l,jintsp) + gjkiln*fj(kj,l,jintsp) ) )
( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) ) )
!+APW
IF (input%l_useapw) THEN
apw1 = REAL(cph(kj)) * plegend(kj,l) * &
......@@ -122,13 +122,13 @@ CONTAINS
!-APW
ENDDO
ELSE
DO kj = 1,ki
DO kj = 1,MIN(ki,lapw%nv(iintsp))
fct = chi*plegend(kj,l)*fl2p1(l)*&
( fjkiln*fj(kj,l,jintsp) + gjkiln*gj(kj,l,jintsp)*ddnln )
( fjkiln*fj(kj,l,iintsp) + gjkiln*gj(kj,l,iintsp)*ddnln )
smat%data_c(kj,kii)=smat%data_c(kj,kii) + cph(kj)*fct
hmat%data_c(kj,kii)=hmat%data_c(kj,kii) + cph(kj) * ( fct*elall &
+ chi*plegend(kj,l)*fl2p1bt(l) * ( fjkiln*gj(kj,l,jintsp) + gjkiln*fj(kj,l,jintsp) ) )
+ chi*plegend(kj,l)*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) ) )
IF (input%l_useapw) THEN
capw1 = cph(kj)*plegend(kj,l)&
......
......@@ -40,15 +40,18 @@ MODULE m_types_mat
SUBROUTINE t_mat_add_transpose(mat,mat1)
CLASS(t_mat),INTENT(INOUT)::mat,mat1
INTEGER::i,j
IF ((mat%matsize1.NE.mat1%matsize2).OR. &
(mat%matsize2.NE.mat1%matsize1)) &
CALL judft_error("Matrix sizes missmatch in add_transpose")
IF (mat%l_real.AND.mat1%l_real) THEN
DO i=1,mat%matsize1
DO j=i+1,mat%matsize2
DO i=1,mat%matsize2
DO j=i+1,mat%matsize1
mat%data_r(j,i)=mat1%data_r(i,j)
ENDDO
ENDDO
ELSEIF((.NOT.mat%l_real).AND.(.NOT.mat1%l_real)) THEN
DO i=1,mat%matsize1
DO j=i+1,mat%matsize2
DO i=1,mat%matsize2
DO j=i+1,mat%matsize1
mat%data_c(j,i)=CONJG(mat1%data_c(i,j))
ENDDO
ENDDO
......
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