Commit 2300e5e0 authored by Daniel Wortmann's avatar Daniel Wortmann

Changes to hybrid code. Mostly more introduction of types, cleanup and bugfixes

parent 8ea58cfa
...@@ -341,7 +341,7 @@ CONTAINS ...@@ -341,7 +341,7 @@ CONTAINS
print *,"Wrong overlap matrix used, fix this later" print *,"Wrong overlap matrix used, fix this later"
call olap%from_packed(l_real,dimension%nbasfcn,hamovlp%b_r,hamovlp%b_c) call olap%from_packed(l_real,dimension%nbasfcn,hamovlp%b_r,hamovlp%b_c)
call write_olap(olap,nrec) call write_olap(olap,nrec)
PRINT *,"BASIS:",lapw%nv(jsp),atoms%nlotot
if (hybrid%l_addhf) CALL add_Vnonlocal(nk,hybrid,dimension, kpts,jsp,results,xcpot,hamovlp) if (hybrid%l_addhf) CALL add_Vnonlocal(nk,hybrid,dimension, kpts,jsp,results,xcpot,hamovlp)
......
...@@ -17,10 +17,11 @@ module m_types_rcmat ...@@ -17,10 +17,11 @@ module m_types_rcmat
END type t_mat END type t_mat
CONTAINS CONTAINS
SUBROUTINE t_mat_alloc(mat,l_real,matsize1,matsize2) SUBROUTINE t_mat_alloc(mat,l_real,matsize1,matsize2,init)
CLASS(t_mat) :: mat CLASS(t_mat) :: mat
LOGICAL,INTENT(IN),OPTIONAL:: l_real LOGICAL,INTENT(IN),OPTIONAL:: l_real
INTEGER,INTENT(IN),OPTIONAL:: matsize1,matsize2 INTEGER,INTENT(IN),OPTIONAL:: matsize1,matsize2
REAL,INTENT(IN),OPTIONAL :: init
INTEGER:: err INTEGER:: err
IF (present(l_real)) mat%l_real=l_real IF (present(l_real)) mat%l_real=l_real
...@@ -36,11 +37,11 @@ module m_types_rcmat ...@@ -36,11 +37,11 @@ module m_types_rcmat
IF (mat%l_real) THEN IF (mat%l_real) THEN
ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2),STAT=err) ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2),STAT=err)
ALLOCATE(mat%data_c(0,0)) ALLOCATE(mat%data_c(0,0))
mat%data_r=0.0 if (present(init)) mat%data_r=init
ELSE ELSE
ALLOCATE(mat%data_r(0,0)) ALLOCATE(mat%data_r(0,0))
ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2),STAT=err) ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2),STAT=err)
mat%data_c=0.0 IF (PRESENT(init)) mat%data_c=init
ENDIF ENDIF
IF (err>0) CALL judft_error("Allocation of memmory failed for mat datatype",hint="You probably run out of memory") IF (err>0) CALL judft_error("Allocation of memmory failed for mat datatype",hint="You probably run out of memory")
......
...@@ -90,12 +90,13 @@ MODULE m_add_vnonlocal ...@@ -90,12 +90,13 @@ MODULE m_add_vnonlocal
DO nn=1,n DO nn=1,n
IF (hamovlp%l_real) THEN IF (hamovlp%l_real) THEN
hamovlp%a_r(ic) = hamovlp%a_r(ic) - a_ex*v_x%data_r(n,nn) hamovlp%a_r(ic) = hamovlp%a_r(ic) - a_ex*v_x%data_r(n,nn)
WRITE(732,*) n,nn,v_x%data_r(n,nn)
ELSE ELSE
hamovlp%a_c(ic) = hamovlp%a_c(ic) - a_ex*v_x%data_c(n,nn) hamovlp%a_c(ic) = hamovlp%a_c(ic) - a_ex*v_x%data_c(n,nn)
ENDIF ENDIF
ENDDO ENDDO
END DO END DO
STOP "DEBUG"
! calculate HF energy ! calculate HF energy
IF( hybrid%l_calhf ) THEN IF( hybrid%l_calhf ) THEN
WRITE(6,'(A)') new_line('n')//new_line('n')//' ### '// ' diagonal HF exchange elements (eV) ###' WRITE(6,'(A)') new_line('n')//new_line('n')//' ### '// ' diagonal HF exchange elements (eV) ###'
......
...@@ -1150,7 +1150,7 @@ CONTAINS ...@@ -1150,7 +1150,7 @@ CONTAINS
DO ikpt = ikptmin,ikptmax DO ikpt = ikptmin,ikptmax
!calculate IR overlap-matrix !calculate IR overlap-matrix
call olapm%alloc(sym%invs,hybrid%ngptm(ikpt),hybrid%ngptm(ikpt)) CALL olapm%alloc(sym%invs,hybrid%ngptm(ikpt),hybrid%ngptm(ikpt),0.0)
CALL olap_pw(olapm,hybrid%gptm(:,hybrid%pgptm(:hybrid%ngptm(ikpt),ikpt)),hybrid%ngptm(ikpt), atoms,cell ) CALL olap_pw(olapm,hybrid%gptm(:,hybrid%pgptm(:hybrid%ngptm(ikpt),ikpt)),hybrid%ngptm(ikpt), atoms,cell )
...@@ -1494,7 +1494,17 @@ CONTAINS ...@@ -1494,7 +1494,17 @@ CONTAINS
# ifdef CPP_IRCOULOMBAPPROX # ifdef CPP_IRCOULOMBAPPROX
call write_coulomb_spm_r(ikpt,coulomb_mt1(:,:,:,:,1),coulomb_mt2_r(:,:,:,:,1),coulomb_mt3_r(:,:,:,1), coulomb_mtir_r(:,1)) call write_coulomb_spm_r(ikpt,coulomb_mt1(:,:,:,:,1),coulomb_mt2_r(:,:,:,:,1),coulomb_mt3_r(:,:,:,1), coulomb_mtir_r(:,1))
# else # else
call write_coulomb_spm_r(ikpt,coulomb_mt1(:,:,:,:,1),coulomb_mt2_r(:,:,:,:,1),coulomb_mt3_r(:,:,:,1), coulombp_mtir_r(:,1)) CALL write_coulomb_spm_r(ikpt,coulomb_mt1(:,:,:,:,1),coulomb_mt2_r(:,:,:,:,1),coulomb_mt3_r(:,:,:,1), coulombp_mtir_r(:,1))
!!$ print *,"DEBUG"
!!$ DO n1=1,SIZE(coulomb_mt1,1)
!!$ DO n2=1,SIZE(coulomb_mt1,2)
!!$ DO i=1,SIZE(coulomb_mt1,3)
!!$ DO j=1,SIZE(coulomb_mt1,4)
!!$ WRITE(732,*) n1,n2,i-1,j,coulomb_mt2_r(n1,n2,i-1,j,1)
!!$ ENDDO
!!$ ENDDO
!!$ ENDDO
!!$ ENDDO
# endif # endif
else else
# ifdef CPP_IRCOULOMBAPPROX # ifdef CPP_IRCOULOMBAPPROX
...@@ -1528,7 +1538,7 @@ CONTAINS ...@@ -1528,7 +1538,7 @@ CONTAINS
CALL cpu_TIME(time2) CALL cpu_TIME(time2)
WRITE(6,'(2X,A,F8.2,A)') '( Timing:',time2-time1,' )' WRITE(6,'(2X,A,F8.2,A)') '( Timing:',time2-time1,' )'
END IF END IF
CONTAINS CONTAINS
! Calculate body of Coulomb matrix at Gamma point: v_IJ = SUM(G) c^*_IG c_JG 4*pi/G**2 . ! Calculate body of Coulomb matrix at Gamma point: v_IJ = SUM(G) c^*_IG c_JG 4*pi/G**2 .
......
...@@ -343,13 +343,12 @@ endif ...@@ -343,13 +343,12 @@ endif
! from IBZ to current k-point ! from IBZ to current k-point
IF( kpts%bkp(ikpt0) .ne. ikpt0 ) THEN IF( kpts%bkp(ikpt0) .ne. ikpt0 ) THEN
CALL bra_trafo2(& CALL bra_trafo2(&
mat_ex%l_real,carr3_vv_r(:hybrid%nbasm(ikpt0),:,:),cprod_vv_r(:hybrid%nbasm(ikpt0),:,:),carr3_vv_c(:hybrid%nbasm(ikpt0),:,:),cprod_vv_c(:hybrid%nbasm(ikpt0),:,:),& mat_ex%l_real,carr3_vv_r(:hybrid%nbasm(ikpt0),:,:),cprod_vv_r(:hybrid%nbasm(ikpt0),:,:),&
hybrid%nbasm(ikpt0),psize,hybrid%nbands(nk),& carr3_vv_c(:hybrid%nbasm(ikpt0),:,:),cprod_vv_c(:hybrid%nbasm(ikpt0),:,:),&
ikpt0,kpts%bkp(ikpt0),kpts%bksym(ikpt0),sym,& hybrid%nbasm(ikpt0),psize,hybrid%nbands(nk),&
hybrid,kpts,cell,hybrid%maxlcutm1,atoms,& kpts%bkp(ikpt0),ikpt0,kpts%bksym(ikpt0),sym,&
hybrid%lcutm1,hybrid%nindxm1,hybrid%maxindxm1,hybrid%gptmd,& hybrid,kpts,cell,atoms,&
hybrid%nbasp,& phase_vv)
phase_vv)
IF (mat_ex%l_real) THEN IF (mat_ex%l_real) THEN
cprod_vv_r(:hybrid%nbasm(ikpt0),:,:) = carr3_vv_r(:hybrid%nbasm(ikpt0),:,:) cprod_vv_r(:hybrid%nbasm(ikpt0),:,:) = carr3_vv_r(:hybrid%nbasm(ikpt0),:,:)
...@@ -546,8 +545,12 @@ endif ...@@ -546,8 +545,12 @@ endif
ENDIF ENDIF
! write exch_vv in mat_ex ! write exch_vv in mat_ex
call mat_ex%alloc(matsize1=hybrid%nbands(nk)) CALL mat_ex%alloc(matsize1=hybrid%nbands(nk))
mat_ex%data_c=exch_vv IF (mat_ex%l_real) THEN
mat_ex%data_r=exch_vv
ELSE
mat_ex%data_c=exch_vv
END IF
END SUBROUTINE exchange_valence_hf END SUBROUTINE exchange_valence_hf
......
...@@ -194,6 +194,7 @@ MODULE m_hsfock ...@@ -194,6 +194,7 @@ MODULE m_hsfock
& noco,nsest,indx_sest,& & noco,nsest,indx_sest,&
& mpi,irank2,isize2,comm,& & mpi,irank2,isize2,comm,&
& ex) & ex)
DEALLOCATE ( rep_c ) DEALLOCATE ( rep_c )
CALL timestop("valence exchange calculation") CALL timestop("valence exchange calculation")
...@@ -267,7 +268,7 @@ MODULE m_hsfock ...@@ -267,7 +268,7 @@ MODULE m_hsfock
CALL timestop("time for performing T^-1*mat_ex*T^-1*") CALL timestop("time for performing T^-1*mat_ex*T^-1*")
CALL symmetrizeh(atoms,& CALL symmetrizeh(atoms,&
& kpts%bkf(:,nk),dimension,jsp,lapw,gpt,& & kpts%bkf(:,nk),dimension,jsp,lapw,gpt,&
& sym,hybdat%kveclo_eig,& & sym,hybdat%kveclo_eig,&
......
...@@ -549,9 +549,7 @@ ...@@ -549,9 +549,7 @@
SUBROUTINE bra_trafo2 (& SUBROUTINE bra_trafo2 (&
l_real,vecout_r,vecin_r,vecout_c,vecin_c,& l_real,vecout_r,vecin_r,vecout_c,vecin_c,&
dim,nobd,nbands,ikpt0,ikpt1,iop,sym,& dim,nobd,nbands,ikpt0,ikpt1,iop,sym,&
hybrid,kpts,cell,maxlcutm,atoms,& hybrid,kpts,cell,atoms,&
lcutm,nindxm,maxindxm,&
ngptmall,nbasp,&
phase) phase)
! ikpt0 :: parent of ikpt1 ! ikpt0 :: parent of ikpt1
...@@ -570,10 +568,9 @@ ...@@ -570,10 +568,9 @@
! - scalars - ! - scalars -
INTEGER,INTENT(IN) :: ikpt0,ikpt1,iop,dim,nobd,nbands INTEGER,INTENT(IN) :: ikpt0,ikpt1,iop,dim,nobd,nbands
INTEGER,INTENT(IN) :: maxlcutm ,maxindxm,nbasp
! - arrays - ! - arrays -
INTEGER,INTENT(IN) :: lcutm(atoms%ntype),nindxm(0:maxlcutm,atoms%ntype)
INTEGER,INTENT(IN) :: ngptmall
LOGICAL,INTENT(IN) :: l_real LOGICAL,INTENT(IN) :: l_real
...@@ -593,12 +590,12 @@ ...@@ -593,12 +590,12 @@
! - arrays - ! - arrays -
INTEGER :: rrot(3,3),invrot(3,3) INTEGER :: rrot(3,3),invrot(3,3)
INTEGER :: pnt(maxindxm,0:maxlcutm,atoms%nat) INTEGER :: pnt(hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%nat)
INTEGER :: g(3),g1(3) INTEGER :: g(3),g1(3)
REAL :: rkpt(3),rkpthlp(3),rtaual(3),trans(3) REAL :: rkpt(3),rkpthlp(3),rtaual(3),trans(3)
REAL :: arg REAL :: arg
COMPLEX :: dwgn(-maxlcutm:maxlcutm,& COMPLEX :: dwgn(-hybrid%maxlcutm1:hybrid%maxlcutm1,&
& -maxlcutm:maxlcutm,0:maxlcutm) & -hybrid%maxlcutm1:hybrid%maxlcutm1,0:hybrid%maxlcutm1)
! COMPLEX :: vecin1(dim,nobd,nbands),vecout1(dim,nobd,nbands) ! COMPLEX :: vecin1(dim,nobd,nbands),vecout1(dim,nobd,nbands)
COMPLEX, ALLOCATABLE :: vecin1(:,:,:),vecout1(:,:,:) COMPLEX, ALLOCATABLE :: vecin1(:,:,:),vecout1(:,:,:)
...@@ -608,15 +605,15 @@ ...@@ -608,15 +605,15 @@
& STOP 'bra_trafo2: error allocating vecin1 or vecout1' & STOP 'bra_trafo2: error allocating vecin1 or vecout1'
vecin1 = 0 ; vecout1 = 0 vecin1 = 0 ; vecout1 = 0
IF( maxlcutm .gt. atoms%lmaxd ) STOP 'bra_trafo2: maxlcutm > atoms%lmaxd' ! very improbable case IF( hybrid%maxlcutm1 .gt. atoms%lmaxd ) STOP 'bra_trafo2: maxlcutm > atoms%lmaxd' ! very improbable case
! transform back to unsymmetrized product basis in case of inversion symmetry ! transform back to unsymmetrized product basis in case of inversion symmetry
if (l_real) THEN if (l_real) THEN
vecin1 = vecin_r vecin1 = vecin_r
DO i=1,nbands DO i=1,nbands
DO j=1,nobd DO j=1,nobd
CALL desymmetrize(vecin1(:nbasp,j,i),nbasp,1,1,& CALL desymmetrize(vecin1(:hybrid%nbasp,j,i),hybrid%nbasp,1,1,&
atoms,lcutm,maxlcutm, nindxm,sym) atoms,hybrid%lcutm1,hybrid%maxlcutm1, hybrid%nindxm1,sym)
END DO END DO
END DO END DO
else else
...@@ -630,8 +627,8 @@ ...@@ -630,8 +627,8 @@
invrot = sym%mrot(:,:,sym%invtab(iop)) invrot = sym%mrot(:,:,sym%invtab(iop))
trans = sym%tau(:,iop) trans = sym%tau(:,iop)
dwgn (-maxlcutm:maxlcutm,-maxlcutm:maxlcutm,0:maxlcutm) & dwgn (-hybrid%maxlcutm1:hybrid%maxlcutm1,-hybrid%maxlcutm1:hybrid%maxlcutm1,0:hybrid%maxlcutm1) &
= hybrid%d_wgn2(-maxlcutm:maxlcutm,-maxlcutm:maxlcutm,0:maxlcutm,inviop) = hybrid%d_wgn2(-hybrid%maxlcutm1:hybrid%maxlcutm1,-hybrid%maxlcutm1:hybrid%maxlcutm1,0:hybrid%maxlcutm1,inviop)
ELSE ELSE
...@@ -641,8 +638,8 @@ ...@@ -641,8 +638,8 @@
invrot = sym%mrot(:,:,sym%invtab(iiop)) invrot = sym%mrot(:,:,sym%invtab(iiop))
trans = sym%tau(:,iiop) trans = sym%tau(:,iiop)
dwgn (-maxlcutm:maxlcutm,-maxlcutm:maxlcutm,0:maxlcutm)& dwgn (-hybrid%maxlcutm1:hybrid%maxlcutm1,-hybrid%maxlcutm1:hybrid%maxlcutm1,0:hybrid%maxlcutm1)&
= conjg( hybrid%d_wgn2(-maxlcutm:maxlcutm, -maxlcutm:maxlcutm, 0:maxlcutm,inviop)) = conjg( hybrid%d_wgn2(-hybrid%maxlcutm1:hybrid%maxlcutm1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1,inviop))
END IF END IF
...@@ -652,7 +649,7 @@ ...@@ -652,7 +649,7 @@
rkpt = modulo1(rkpt,kpts%nkpt3) rkpt = modulo1(rkpt,kpts%nkpt3)
g = nint(rkpthlp-rkpt) g = nint(rkpthlp-rkpt)
#ifdef CPP_DEBUG !#ifdef CPP_DEBUG
!test !test
DO i=1,kpts%nkpt DO i=1,kpts%nkpt
IF ( maxval( abs(rkpt - kpts%bk(:,i)) ) .le. 1E-06 ) THEN IF ( maxval( abs(rkpt - kpts%bk(:,i)) ) .le. 1E-06 ) THEN
...@@ -661,7 +658,7 @@ ...@@ -661,7 +658,7 @@
END IF END IF
END DO END DO
IF( nrkpt .ne. ikpt1 ) STOP 'bra_trafo2: rotation failed' IF( nrkpt .ne. ikpt1 ) STOP 'bra_trafo2: rotation failed'
#endif !#endif
! Define pointer to first mixed-basis functions (with m = -l) ! Define pointer to first mixed-basis functions (with m = -l)
i = 0 i = 0
...@@ -669,12 +666,12 @@ ...@@ -669,12 +666,12 @@
DO itype = 1,atoms%ntype DO itype = 1,atoms%ntype
DO ieq = 1,atoms%neq(itype) DO ieq = 1,atoms%neq(itype)
ic = ic + 1 ic = ic + 1
DO l = 0,lcutm(itype) DO l = 0,hybrid%lcutm1(itype)
DO n = 1,nindxm(l,itype) DO n = 1,hybrid%nindxm1(l,itype)
i = i + 1 i = i + 1
pnt(n,l,ic) = i pnt(n,l,ic) = i
END DO END DO
i = i + nindxm(l,itype) * 2*l i = i + hybrid%nindxm1(l,itype) * 2*l
END DO END DO
END DO END DO
END DO END DO
...@@ -692,8 +689,8 @@ ...@@ -692,8 +689,8 @@
cdum = cexp *exp(-img*tpi_const*dot_product(g,atoms%taual(:,rcent))) cdum = cexp *exp(-img*tpi_const*dot_product(g,atoms%taual(:,rcent)))
DO l = 0,lcutm(itype) DO l = 0,hybrid%lcutm1(itype)
nn = nindxm(l,itype) nn = hybrid%nindxm1(l,itype)
DO n = 1,nn DO n = 1,nn
i1 = pnt(n,l,ic) i1 = pnt(n,l,ic)
...@@ -729,13 +726,18 @@ ...@@ -729,13 +726,18 @@
WRITE(*,*) ikpt0,ikpt1,g1 WRITE(*,*) ikpt0,ikpt1,g1
WRITE(*,*) hybrid%ngptm(ikpt0),hybrid%ngptm(ikpt1) WRITE(*,*) hybrid%ngptm(ikpt0),hybrid%ngptm(ikpt1)
WRITE(*,*) WRITE(*,*)
WRITE(*,*) hybrid%gptm(:,igptp) WRITE(*,*) igptp,hybrid%gptm(:,igptp)
WRITE(*,*) g WRITE(*,*) g
WRITE(*,*) rrot
WRITE (*,*) "Failed tests:",g1
DO i=1,hybrid%ngptm(ikpt1)
WRITE(*,*) hybrid%gptm(:,hybrid%pgptm(i,ikpt1))
ENDDO
STOP 'bra_trafo2: G-point not found in G-point set.' STOP 'bra_trafo2: G-point not found in G-point set.'
END IF END IF
cdum = exp(img*tpi_const*dot_product(kpts%bkf(:,ikpt1)+g1,trans(:))) cdum = exp(img*tpi_const*dot_product(kpts%bkf(:,ikpt1)+g1,trans(:)))
vecout1(nbasp+igptm,:,:)= cdum * vecin1(nbasp+igptm2,:,:) vecout1(hybrid%nbasp+igptm,:,:)= cdum * vecin1(hybrid%nbasp+igptm2,:,:)
END DO END DO
DEALLOCATE ( vecin1 ) DEALLOCATE ( vecin1 )
...@@ -745,7 +747,7 @@ ...@@ -745,7 +747,7 @@
DO j=1,nobd DO j=1,nobd
CALL symmetrize(vecout1(:,j,i),dim,1,1,.false.,& CALL symmetrize(vecout1(:,j,i),dim,1,1,.false.,&
atoms,lcutm,maxlcutm, nindxm,sym) atoms,hybrid%lcutm1,hybrid%maxlcutm1, hybrid%nindxm1,sym)
CALL commonphase(phase(j,i),vecout1(:,j,i),dim) CALL commonphase(phase(j,i),vecout1(:,j,i),dim)
vecout1(:,j,i) = vecout1(:,j,i) / phase(j,i) vecout1(:,j,i) = vecout1(:,j,i) / phase(j,i)
......
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