Commit 5874820a authored by Daniel Wortmann's avatar Daniel Wortmann
Browse files

Merge remote-tracking branch 'origin/develop' into strgn

parents ff5e7d6a 3f84f373
......@@ -443,10 +443,10 @@ CONTAINS
!
!=====> calculate the fourier transform of the core-pseudocharge
IF (l_f2) THEN
CALL ft_of_CorePseudocharge(fmpi,atoms,mshc,alpha,tol_14,rh, &
CALL ft_of_CorePseudocharge(fmpi,input,atoms,mshc,alpha,tol_14,rh, &
acoff,stars,method2,rat,cell,oneD,sym,qpwc,jspin,l_f2,vpw,ffonat,force_a4_mt_loc)
ELSE
CALL ft_of_CorePseudocharge(fmpi,atoms,mshc,alpha,tol_14,rh, &
CALL ft_of_CorePseudocharge(fmpi,input,atoms,mshc,alpha,tol_14,rh, &
acoff,stars,method2,rat,cell,oneD,sym,qpwc,jspin,l_f2)
END IF
......@@ -687,7 +687,7 @@ CONTAINS
! INTERNAL SUBROUTINES
!***********************************************************************
subroutine ft_of_CorePseudocharge(fmpi,atoms,mshc,alpha,&
subroutine ft_of_CorePseudocharge(fmpi,input,atoms,mshc,alpha,&
tol_14,rh,acoff,stars,method2,rat,cell,oneD,sym,qpwc,jspin,l_f2,vpw,ffonat,force_a4_mt_loc)
!=====> calculate the fourier transform of the core-pseudocharge
......@@ -700,7 +700,7 @@ CONTAINS
USE m_types
type(t_mpi) ,intent(in) :: fmpi
TYPE(t_input), INTENT(in) ::input
type(t_atoms) ,intent(in) :: atoms
integer ,intent(in) :: mshc(atoms%ntype),jspin
real ,intent(in) :: alpha(atoms%ntype), tol_14
......@@ -719,6 +719,7 @@ CONTAINS
! ..Local variables
integer nat1, n, n_out_p, k
INTEGER :: reducedStarsCutoff ! This is introduced to avoid numerical instabilities.
complex czero
! ..Local arrays
......@@ -730,11 +731,12 @@ CONTAINS
#endif
czero = (0.0,0.0)
#ifdef CPP_MPI
DO k = 1 , stars%ng3
DO k = 1, stars%ng3
qpwc_loc(k) = czero
ENDDO
#endif
DO k = 1 , stars%ng3
DO k = 1, stars%ng3
IF (stars%sk3(k).LE.3.0*input%rkmax) reducedStarsCutoff = k ! The factor 3.0 is arbitrary. One could try going down to 2.0.
qpwc(k) = czero
ENDDO
......@@ -749,7 +751,7 @@ CONTAINS
! (1) Form factor for each atom type
CALL FormFactor_forAtomType(atoms%msh,method2,n_out_p,&
CALL FormFactor_forAtomType(atoms%msh,method2,n_out_p,reducedStarsCutoff,&
atoms%rmt(n),atoms%jri(n),atoms%dx(n),mshc(n),rat(:,n), &
rh(:,n),alpha(n),stars,cell,acoff(n),qf)
......@@ -763,11 +765,11 @@ CONTAINS
END IF
IF (l_f2) THEN
CALL StructureConst_forAtom(nat1,stars,oneD,sym,&
CALL StructureConst_forAtom(nat1,stars,oneD,sym,reducedStarsCutoff,&
atoms%neq(n),atoms%nat,atoms%taual,&
cell,qf,qpwc_at,jspin,l_f2,n,vpw,ffonat)
ELSE
CALL StructureConst_forAtom(nat1,stars,oneD,sym,&
CALL StructureConst_forAtom(nat1,stars,oneD,sym,reducedStarsCutoff,&
atoms%neq(n),atoms%nat,atoms%taual,&
cell,qf,qpwc_at,jspin,l_f2,n)
END IF
......@@ -796,7 +798,7 @@ CONTAINS
END IF
end subroutine ft_of_CorePseudocharge
SUBROUTINE StructureConst_forAtom(nat1,stars,oneD,sym,&
SUBROUTINE StructureConst_forAtom(nat1,stars,oneD,sym,reducedStarsCutoff,&
neq,natd,taual,cell,qf,qpwc_at,jspin,l_f2,n,vpw,ffonat)
! Calculates the structure constant for each atom of atom type
......@@ -809,6 +811,7 @@ CONTAINS
type(t_stars), intent(in) :: stars
type(t_oneD), intent(in) :: oneD
type(t_sym), intent(in) :: sym
INTEGER, INTENT(IN) :: reducedStarsCutoff
integer, intent(in) :: neq,natd, jspin, n
real, intent(in) :: taual(3,natd)
type(t_cell), intent(in) :: cell
......@@ -829,9 +832,7 @@ CONTAINS
complex phaso(oneD%ods%nop), kcmplx(3)
czero = (0.0,0.0)
DO k = 1 , stars%ng3
qpwc_at(k) = czero
END DO
qpwc_at(:) = czero
! first G=0
k=1
......@@ -842,11 +843,11 @@ CONTAINS
force_mt_loc=0.0
force_is_loc=cmplx(0.0,0.0)
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP SHARED(stars,oneD,sym,neq,natd,nat1,taual,cell,qf,qpwc_at,l_f2,ffonat,n,jspin,vpw) &
!$OMP SHARED(stars,oneD,sym,reducedStarsCutoff,neq,natd,nat1,taual,cell,qf,qpwc_at,l_f2,ffonat,n,jspin,vpw) &
!$OMP FIRSTPRIVATE(czero) &
!$OMP PRIVATE(k,kr,phas,nat2,nat,sf,j,x,kro,phaso,kcmplx,phase) &
!$OMP REDUCTION(+:force_mt_loc,force_is_loc)
DO k = 2,stars%ng3
DO k = 2,reducedStarsCutoff
IF (.NOT.oneD%odi%d1) THEN
CALL spgrot(sym%nop, sym%symor, sym%mrot, sym%tau, sym%invtab, &
stars%kv3(:,k), kr, phas)
......@@ -912,7 +913,7 @@ CONTAINS
END IF
END SUBROUTINE StructureConst_forAtom
SUBROUTINE FormFactor_forAtomType(msh, method2, n_out_p, rmt, jri, dx, &
SUBROUTINE FormFactor_forAtomType(msh, method2, n_out_p, reducedStarsCutoff, rmt, jri, dx, &
mshc, rat, rh, alpha, stars, cell, acoff, &
qf)
......@@ -923,6 +924,7 @@ CONTAINS
integer ,intent(in) :: msh,method2, n_out_p
INTEGER, INTENT(IN) :: reducedStarsCutoff
real ,intent(in) :: rmt
integer ,intent(in) :: jri
real ,intent(in) :: dx
......@@ -944,9 +946,7 @@ CONTAINS
real rhohelp(msh)
zero = 0.0
DO k = 1,stars%ng3
qf(k) = 0.0
END DO
qf(:) = 0.0
tail = .FALSE.
f11 = tpi_const * rmt * rh(jri) / alpha
......@@ -954,11 +954,11 @@ CONTAINS
ar = SQRT( alpha ) * rmt
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP SHARED(stars,f11,f12,ar,method2,n_out_p,jri,rat,rh,dx,tail) &
!$OMP SHARED(stars,f11,f12,ar,method2,n_out_p,reducedStarsCutoff,jri,rat,rh,dx,tail) &
!$OMP SHARED(alpha,cell,mshc,rmt,qf) &
!$OMP FIRSTPRIVATE(zero) &
!$OMP PRIVATE(k,g,ai,qfin,ir,j,rhohelp,qfout,gr,a4,alpha3)
DO k = 1,stars%ng3
DO k = 1, reducedStarsCutoff
g = stars%sk3(k)
! first G=0
IF ( k.EQ.1 ) THEN
......
......@@ -297,7 +297,7 @@ SUBROUTINE cdnval(eig_id, fmpi,kpts,jspin,noco,nococonv,input,banddos,cell,atoms
IF (noco%l_soc) CALL orbmom(atoms,noccbd,we,ispin,eigVecCoeffs,orb)
IF (input%l_f) THEN
CALL tlmplm%init(atoms,input%jspins,.FALSE.)
CALL tlmplm_cholesky(sphhar,atoms,sym,noco,nococonv,enpara,ispin,fmpi,vTot,vtot,den,input,hub1inp,hub1data,tlmplm,usdus,0.0,.FALSE.)
CALL tlmplm_cholesky(sphhar,atoms,sym,noco,nococonv,enpara,ispin,fmpi,vTot,vtot,den,input,hub1inp,hub1data,tlmplm,usdus,0.0)
CALL force%addContribsA21A12(input,atoms,sym,cell,oneD,enpara,&
usdus,tlmplm,vtot,eigVecCoeffs,noccbd,ispin,eig,we,results,jsp_start,jspin,nbasfcn,zMat,lapw,sphhar,lapw%gvec(1,:,:),lapw%gvec(2,:,:),lapw%gvec(3,:,:),bkpt)
ENDIF
......
......@@ -133,7 +133,7 @@ CONTAINS
CALL setabc1lo(atoms,iType,usdus,jspin,alo1,blo1,clo1)
! generate the spinors (chi)
IF(noco%l_noco) ccchi=nococonv%chi(itype)
IF(noco%l_noco) ccchi=conjg(nococonv%umat(itype))
nintsp = 1
......@@ -212,7 +212,7 @@ CONTAINS
!$acc host_data use_device(work_c,abCoeffs,abTemp)
CALL zgemm_acc("T","C",ne,2*abSize,nvmax,CMPLX(1.0,0.0),work_c,MAXVAL(lapw%nv),abCoeffs,2*atoms%lmaxd*(atoms%lmaxd+2)+2,CMPLX(0.0,0.0),abTemp,acof_size)
CALL zgemm_acc("T","T",ne,2*abSize,nvmax,CMPLX(1.0,0.0),work_c,MAXVAL(lapw%nv),abCoeffs,2*atoms%lmaxd*(atoms%lmaxd+2)+2,CMPLX(0.0,0.0),abTemp,acof_size)
!$acc end host_data
!stop "DEBUG"
......@@ -302,8 +302,8 @@ CONTAINS
lm = ll1 + m
lmp = ll1 - m
inv_f = (-1)**(l-m)
acof(:,lmp,jatom)=acof(:,lmp,jatom)+inv_f*matmul(abCoeffs(lm+1,:),work_c(:nvmax,:))
bcof(:,lmp,jatom)=bcof(:,lmp,jatom)+inv_f*matmul(abCoeffs(lm+1+abSize,:),work_c(:nvmax,:))
acof(:,lmp,jatom)=acof(:,lmp,jatom)+inv_f*matmul(CONJG(abCoeffs(lm+1,:)),work_c(:nvmax,:)) !TODO: Is this conjugation costly?
bcof(:,lmp,jatom)=bcof(:,lmp,jatom)+inv_f*matmul(CONJG(abCoeffs(lm+1+abSize,:)),work_c(:nvmax,:)) !TODO: Is this conjugation costly?
!CALL zaxpy(ne,c_1,workTrans_c(:,iLAPW),1, acof(:,lmp,jatom),1)
!CALL zaxpy(ne,c_2,workTrans_c(:,iLAPW),1, bcof(:,lmp,jatom),1)
END DO
......@@ -336,10 +336,10 @@ CONTAINS
fgpl(:,iLAPW) = MATMUL(fgr,cell%bmat)
ENDDO
helpMat_c = abCoeffs(1+abSize:,:)
helpMat_c = CONJG(abCoeffs(1+abSize:,:)) !TODO: Is this conjugation costly?
workTrans_cf = 0.0
CALL zgemm("N","C",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),s2h_e,ne,abCoeffs,size(abcoeffs,1),CMPLX(1.0,0.0),force%e1cof(:,:,iAtom),ne)
CALL zgemm("N","T",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),s2h_e,ne,abCoeffs,size(abcoeffs,1),CMPLX(1.0,0.0),force%e1cof(:,:,iAtom),ne)
CALL zgemm("N","C",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),s2h_e,ne,helpMat_c,size(helpMat_c,1),CMPLX(1.0,0.0),force%e2cof(:,:,iAtom),ne)
DO i =1,3
IF (zmat%l_real) THEN
......@@ -351,7 +351,7 @@ CONTAINS
workTrans_cf(:,iLAPW) = workTrans_c(:,iLAPW) * fgpl(i,iLAPW)
ENDDO
ENDIF
CALL zgemm("N","C",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),workTrans_cf,ne,abCoeffs,size(abCoeffs,1),CMPLX(0.0,0.0),helpMat_force,ne)
CALL zgemm("N","T",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),workTrans_cf,ne,abCoeffs,size(abCoeffs,1),CMPLX(0.0,0.0),helpMat_force,ne)
force%aveccof(i,:,:,iAtom) = force%aveccof(i,:,:,iAtom) + helpMat_force(:,:)
CALL zgemm("N","C",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),workTrans_cf,ne,helpMat_c,size(helpMat_c,1),CMPLX(0.0,0.0),helpMat_force,ne)
force%bveccof(i,:,:,iAtom) = force%bveccof(i,:,:,iAtom) + helpMat_force(:,:)
......@@ -363,8 +363,8 @@ CONTAINS
ll1 = l* (l+1)
DO m = -l,l
lm = ll1 + m
c_1 = CONJG(abCoeffs(lm+1,iLAPW))
c_2 = CONJG(abCoeffs(lm+1+abSize,iLAPW))
c_1 = abCoeffs(lm+1,iLAPW)
c_2 = abCoeffs(lm+1+abSize,iLAPW)
jatom = sym%invsatnr(iAtom)
lmp = ll1 - m
inv_f = (-1)**(l-m)
......
......@@ -36,6 +36,13 @@ if (CMAKE_Fortran_COMPILER_ID MATCHES "Intel")
else()
set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -check arg_temp_created,assume,bounds,contiguous,format,output_conversion,pointers,stack,uninit -traceback -O0 -g -check uninit -check pointers -DCPP_DEBUG")
endif()
elseif (CMAKE_Fortran_COMPILER_ID MATCHES "NVHPC")
message("NVHPC Fortran detected")
set(FLEUR_PRECISION_OPTION "-Mr8;-Mr8intrinsics")
set(FLEUR_OPENMP_FLAG "-mp")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -tp=zen2 -mp -O1 -g ")
set(FLEUR_COMPILE_OPTIONS -mavx2 -Mlre -Mautoinline -Mpre -Mvect=simd -Mcache_align -Mflushz -O2 -g)
set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -C -traceback -O0 -g -Mchkstk -Mchkptr -Ktrap=fp -DCPP_DEBUG")
elseif (CMAKE_Fortran_COMPILER_ID MATCHES "PGI")
message("PGI Fortran detected")
set(FLEUR_PRECISION_OPTION "-Mr8;-Mr8intrinsics")
......
......@@ -30,6 +30,7 @@ eigen/tlmplm.F90
eigen/tlo.f90
eigen/vacfun.F90
eigen/eigen_redist_matrix.f90
eigen/hs_int_direct.F90
)
#if (FLEUR_USE_GPU)
# set(fleur_F90 ${fleur_F90}
......
......@@ -20,9 +20,18 @@ CONTAINS
!> 4. writing (saving) of eigenvectors
!>
!>@author D. Wortmann
!
! Modifications done to use this with DFPT phonons:
! a) We need additional MT-integrals from mt_setup that cover the potential variation V1.
! b) The eigenvalues are to be evaluated for k+q, not k.
! c) Additionally, load in the occupied states for k without q.
! d) The work isn't done once the eigenvectors and eigenvalues are found. There is post-
! processing to gain the perturbed eigenvalues and eigenvectors.
! e) The latter are the actual output for the routine when used for DFPT. They are saved
! the same way as the eigenvalues before, but for a shifted eig_id.
SUBROUTINE eigen(fi,fmpi,stars,sphhar,xcpot,&
enpara,nococonv,mpdata,hybdat,&
iter,eig_id,results,inden,v,vx,hub1data,nvfull,GbasVec_eig,z0)
iter,eig_id,results,inden,v,vx,hub1data,nvfull,GbasVec_eig,bqpt,dfpt_eig_id,starsq,v1real,v1imag)
#include"cpp_double.h"
USE m_types
......@@ -34,13 +43,14 @@ CONTAINS
USE m_mt_setup
USE m_util
!USE m_icorrkeys
USE m_eig66_io, ONLY : open_eig, write_eig, read_eig
USE m_eig66_io, ONLY : write_eig, read_eig
USE m_xmlOutput
USE m_symmetrize_matrix
USE m_unfold_band_kpts !used for unfolding bands
USE m_types_mpimat
use m_store_load_hybrid
USE m_dfpt_tlmplm
IMPLICIT NONE
......@@ -60,20 +70,25 @@ CONTAINS
TYPE(t_potden), INTENT(IN) :: vx
TYPE(t_potden),INTENT(IN) :: v
INTEGER, OPTIONAL, ALLOCATABLE, INTENT(OUT) :: nvfull(:, :), GbasVec_eig(:, :, :, :)
COMPLEX, OPTIONAL, ALLOCATABLE, INTENT(OUT) :: z0(:, :, :, :)
! EXTERNAL MPI_BCAST !only used by band_unfolding to broadcast the gvec
! Scalar Arguments
INTEGER,INTENT(IN) :: iter
INTEGER,INTENT(IN) :: eig_id
INTEGER, OPTIONAL, ALLOCATABLE, INTENT(OUT) :: nvfull(:, :), GbasVec_eig(:, :, :, :)
REAL, OPTIONAL, INTENT(IN) :: bqpt
INTEGER, OPTIONAL, INTENT(IN) :: dfpt_eig_id
TYPE(t_stars), OPTIONAL, INTENT(IN) :: starsq
TYPE(t_potden), OPTIONAL, INTENT(IN) :: v1real, v1imag
! Local Scalars
INTEGER jsp,nk,nred,ne_all,ne_found,neigd2
INTEGER ne, nk_i,n_size,n_rank
INTEGER isp,i,j,err
LOGICAL l_wu,l_file,l_real,l_zref
LOGICAL l_real,l_zref
INTEGER :: solver=0
! Local Arrays
INTEGER :: ierr
......@@ -82,33 +97,44 @@ CONTAINS
COMPLEX :: unfoldingBuffer(SIZE(results%unfolding_weights,1),fi%kpts%nkpt,fi%input%jspins) ! needed for unfolding bandstructure fmpi case
INTEGER, ALLOCATABLE :: nvBuffer(:,:), nvBufferTemp(:,:), nvfullBuffer(:,:), GbasVecBuffer(:, :, :, :)
COMPLEX, ALLOCATABLE :: z0Buffer(:,:,:,:)
REAL, ALLOCATABLE :: bkpt(:)
REAL, ALLOCATABLE :: eig(:), eigBuffer(:,:,:)
INTEGER :: jsp_m, i_kpt_m, i_m
INTEGER :: maxspin
TYPE(t_tlmplm) :: td
TYPE(t_usdus) :: ud
TYPE(t_tlmplm) :: td, tdV1, tdmod
TYPE(t_usdus) :: ud, uddummy
TYPE(t_lapw) :: lapw
CLASS(t_mat), ALLOCATABLE :: zMat
CLASS(t_mat), ALLOCATABLE :: hmat,smat
CLASS(t_mat), ALLOCATABLE :: smat_unfold !used for unfolding bandstructure
TYPE(t_kpts) :: kqpts ! basically kpts, but with q added onto each one.
TYPE(t_hub1data) :: hub1datadummy
! Variables for HF or fi%hybinp functional calculation
INTEGER :: comm(fi%kpts%nkpt),irank2(fi%kpts%nkpt),isize2(fi%kpts%nkpt), dealloc_stat
character(len=300) :: errmsg
real :: alpha_hybrid
LOGICAL :: l_dfpteigen
l_dfpteigen = PRESENT(bqpt)
kqpts = fi%kpts
! Modify this from kpts only in DFPT case.
IF (l_dfpteigen) THEN
DO nk_i = 1, fi%kpts%nkpt
kqpts%bk(3, nk_i) = kqpts%bk(3, nk_i) + bqpt
END DO
END IF
call ud%init(fi%atoms,fi%input%jspins)
call uddummy%init(fi%atoms,fi%input%jspins)
ALLOCATE(eig(fi%input%neig))
ALLOCATE(bkpt(3))
ALLOCATE(eigBuffer(fi%input%neig,fi%kpts%nkpt,fi%input%jspins))
ALLOCATE(nvBuffer(fi%kpts%nkpt,MERGE(1,fi%input%jspins,fi%noco%l_noco)),nvBufferTemp(fi%kpts%nkpt,MERGE(1,fi%input%jspins,fi%noco%l_noco)))
! check if z-reflection trick can be used
l_zref=(fi%sym%zrfs.AND.(SUM(ABS(fi%kpts%bk(3,:fi%kpts%nkpt))).LT.1e-9).AND..NOT.fi%noco%l_noco)
l_zref=(fi%sym%zrfs.AND.(SUM(ABS(kqpts%bk(3,:fi%kpts%nkpt))).LT.1e-9).AND..NOT.fi%noco%l_noco)
IF (fmpi%n_size > 1) l_zref = .FALSE.
!IF (fmpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber '/),(/iter,v%iter/),&
......@@ -119,7 +145,13 @@ CONTAINS
! set up k-point independent t(l'm',lm) matrices
alpha_hybrid = MERGE(xcpot%get_exchange_weight(),0.0,hybdat%l_subvxc)
CALL mt_setup(fi%atoms,fi%sym,sphhar,fi%input,fi%noco,nococonv,enpara,fi%hub1inp,hub1data,inden,v,vx,fmpi,results,td,ud,alpha_hybrid,.FALSE.)
CALL mt_setup(fi%atoms,fi%sym,sphhar,fi%input,fi%noco,nococonv,enpara,fi%hub1inp,hub1data,inden,v,vx,fmpi,td,ud,alpha_hybrid)
! Get matrix elements of perturbed potential and modified H/S in DFPT case.
IF (l_dfpteigen) THEN
hub1datadummy = hub1data
CALL dfpt_tlmplm(fi%atoms,fi%sym,sphhar,fi%input,fi%noco,enpara,fi%hub1inp,hub1data,v,fmpi,tdV1,v1real,v1imag)
CALL mt_setup(fi%atoms,fi%sym,sphhar,fi%input,fi%noco,nococonv,enpara,fi%hub1inp,hub1datadummy,inden,v,vx,fmpi,tdmod,uddummy,alpha_hybrid,.TRUE.)
END IF
neigBuffer = 0
results%neig = 0
......@@ -138,17 +170,14 @@ CONTAINS
GbasVec_eig = 0
ALLOCATE(GbasVecBuffer(3, fi%input%neig, fi%kpts%nkpt, MERGE(1,fi%input%jspins,fi%noco%l_noco)))
GbasVecBuffer = 0
ALLOCATE(z0(fi%input%neig, fi%input%neig, fi%kpts%nkpt, fi%input%jspins))
z0 = CMPLX(0.0, 0.0)
ALLOCATE(z0Buffer(fi%input%neig, fi%input%neig, fi%kpts%nkpt, fi%input%jspins))
z0Buffer = CMPLX(0.0, 0.0)
END IF
DO jsp = 1, MERGE(1,fi%input%jspins,fi%noco%l_noco)
k_loop:DO nk_i = 1,size(fmpi%k_list)
nk=fmpi%k_list(nk_i)
! Set up lapw list
CALL lapw%init(fi%input,fi%noco,nococonv, fi%kpts, fi%atoms, fi%sym, nk, fi%cell, l_zref, fmpi)
CALL lapw%init(fi%input,fi%noco,nococonv, kqpts, fi%atoms, fi%sym, nk, fi%cell, l_zref, fmpi)
call timestart("Setup of H&S matrices")
CALL eigen_hssetup(jsp,fmpi,fi,mpdata,results,vx,xcpot,enpara,nococonv,stars,sphhar,hybdat,ud,td,v,lapw,nk,smat,hmat)
CALL timestop("Setup of H&S matrices")
......@@ -160,13 +189,12 @@ CONTAINS
GbasVecBuffer(:, :lapw%nv(jsp), nk, jsp) = lapw%gvec(:, :lapw%nv(jsp), jsp)
END IF
l_wu=.FALSE.
ne_all=fi%input%neig
IF(ne_all < 0) ne_all = lapw%nmat
IF(ne_all > lapw%nmat) ne_all = lapw%nmat
!Try to symmetrize matrix
CALL symmetrize_matrix(fmpi,fi%noco,fi%kpts,nk,hmat,smat)
CALL symmetrize_matrix(fmpi,fi%noco,kqpts,nk,hmat,smat)
IF (fi%banddos%unfoldband .AND. (.NOT. fi%noco%l_soc)) THEN
select type(smat)
......@@ -175,19 +203,21 @@ CONTAINS
select type(smat_unfold)
type is (t_mat)
smat_unfold=smat
smat_unfold%data_c=CONJG(smat%data_c)
end select
type is (t_mpimat)
allocate(t_mpimat::smat_unfold)
select type(smat_unfold)
type is (t_mpimat)
smat_unfold=smat
smat_unfold%data_c=CONJG(smat%data_c)
end select
end select
END IF
! Solve generalized eigenvalue problem.
! ne_all ... number of eigenpairs searched (and found) on this node
! on fi%input, overall number of eigenpairs searched,
! on input, overall number of eigenpairs searched,
! on output, local number of eigenpairs found
! eig ...... all eigenvalues, output
! zMat ..... local eigenvectors, output
......@@ -202,15 +232,6 @@ CONTAINS
if(dealloc_stat /= 0) call juDFT_error("deallocate failed for hmat or smat",&
hint=errmsg, calledby="eigen.F90")
! TODO: Should this go before or after the complex conjugation some lines below?
IF (PRESENT(nvfull)) THEN
IF (zMat%l_real) THEN
z0Buffer(:lapw%nv(jsp), :lapw%nv(jsp), nk, jsp) = CMPLX(1.0,0.0)*zMat%data_r(:lapw%nv(jsp), :lapw%nv(jsp))
ELSE
z0Buffer(:lapw%nv(jsp), :lapw%nv(jsp), nk, jsp) = zMat%data_c(:lapw%nv(jsp), :lapw%nv(jsp))
END IF
END IF
! Output results
CALL timestart("EV output")
ne_found=ne_all
......@@ -220,9 +241,6 @@ CONTAINS
CALL MPI_ALLREDUCE(ne_found,ne_all,1,MPI_INTEGER,MPI_SUM, fmpi%diag_sub_comm,ierr)
ne_all=MIN(fi%input%neig,ne_all)
#endif
IF (.NOT.zMat%l_real) THEN
zMat%data_c(:lapw%nmat,:ne_found) = CONJG(zMat%data_c(:lapw%nmat,:ne_found))
END IF
endif
IF (fmpi%n_rank == 0) THEN
......@@ -234,12 +252,22 @@ CONTAINS
#else
n_rank = 0; n_size=1;
#endif
CALL write_eig(eig_id, nk,jsp,ne_found,ne_all,&
eig(:ne_all),n_start=n_size,n_end=n_rank,zMat=zMat)
!IF (.NOT.l_dfpteigen) THEN
CALL write_eig(eig_id, nk,jsp,ne_found,ne_all,&
eig(:ne_all),n_start=n_size,n_end=n_rank,zMat=zMat)
!ELSE
! CALL dfpt_eigen()
! RETURN
!END IF
eigBuffer(:ne_all,nk,jsp) = eig(:ne_all)
ELSE
if (fmpi%pe_diag) CALL write_eig(eig_id, nk,jsp,ne_found,&
n_start=fmpi%n_size,n_end=fmpi%n_rank,zMat=zMat)
!IF (.NOT.l_dfpteigen) THEN
if (fmpi%pe_diag) CALL write_eig(eig_id, nk,jsp,ne_found,&
n_start=fmpi%n_size,n_end=fmpi%n_rank,zMat=zMat)
!ELSE
! if (fmpi%pe_diag) CALL dfpt_eigen()
! RETURN
!END IF
ENDIF
neigBuffer(nk,jsp) = ne_found
#if defined(CPP_MPI)
......@@ -277,9 +305,6 @@ CONTAINS
IF (PRESENT(nvfull)) THEN
CALL MPI_ALLREDUCE(nvfullBuffer(:,:),nvfull(:,:),size(nvfullBuffer),MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,ierr)
CALL MPI_ALLREDUCE(GbasVecBuffer(:,:,:,:),GbasVec_eig(:,:,:,:),size(GbasVecBuffer),MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,ierr)
!CALL MPI_ALLREDUCE(z0Buffer(:,:,:,:),z0(:,:,:,:),size(z0),MPI_COMPLEX,MPI_SUM,MPI_COMM_WORLD,ierr)
! TODO: This last broadcast is surely faulty.
z0(:,:,:,:) = z0Buffer(:,:,:,:)
END IF
CALL MPI_BARRIER(fmpi%MPI_COMM,ierr)
#else
......@@ -290,7 +315,6 @@ CONTAINS
IF (PRESENT(nvfull)) THEN
nvfull(:,:) = nvfullBuffer(:,:)
GbasVec_eig(:,:,:,:) = GbasVecBuffer(:,:,:,:)
z0(:,:,:,:) = z0Buffer(:,:,:,:)
END IF
#endif
......
......@@ -70,7 +70,7 @@ CONTAINS
!Generate interstitial part of Hamiltonian
ALLOCATE(vpw_wTemp(SIZE(v%pw_w,1),SIZE(v%pw_w,2)))
vpw_wTemp = merge(v%pw_w - xcpot%get_exchange_weight() * vx%pw_w, v%pw_w, hybdat%l_subvxc)
CALL hs_int(fi%input, fi%noco, stars, lapw, fmpi, fi%cell, isp, vpw_wTemp, smat, hmat)
CALL hs_int(fi%input, fi%noco, nococonv, stars, lapw, fmpi, fi%cell%bbmat, isp, vpw_wTemp, smat, hmat)
DEALLOCATE(vpw_wTemp)
CALL timestop("Interstitial part")
......@@ -103,8 +103,9 @@ CONTAINS
if (any(shape(smat) /= 1)) then
call judft_error("Hybrid doesn't do noco.")
end if
smat(1,1)%data_c = CONJG(smat(1,1)%data_c)
CALL write_eig(hybdat%eig_id, nk, isp, smat=smat(1, 1), n_start=fmpi%n_size, n_end=fmpi%n_rank)
smat(1,1)%data_c = CONJG(smat(1,1)%data_c)
END IF
IF (fi%hybinp%l_hybrid) THEN
......
This diff is collapsed.
......@@ -6,83 +6,77 @@
MODULE m_hs_int
CONTAINS
!Subroutine to construct the interstitial Hamiltonian and overlap matrix
SUBROUTINE hs_int(input,noco,stars,lapw,fmpi,cell,isp,vpw,&
smat,hmat)
USE m_types
IMPLICIT NONE
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_lapw),INTENT(IN) :: lapw
TYPE(t_mpi),INTENT(IN) :: fmpi
INTEGER,INTENT(IN) :: isp
COMPLEX,INTENT(IN) :: vpw(:,:)
CLASS(t_mat),INTENT(INOUT) :: smat(:,:),hmat(:,:)
!Subroutine to construct the interstitial Hamiltonian and overlap matrix
SUBROUTINE hs_int(input, noco, nococonv, stars, lapw, fmpi, bbmat, isp, vpw, &
& smat, hmat)
! Control subroutine for the calculation of Hamiltonian/overlap matrix
! elements in the interstitial. The spin logic and case selections are
! found here while the actual calculation loop is one layer deeper in
! hs_int_direct.
!
! <\Phi'|H/S|\Phi>
!
! All primed indices are to be understood as corresponding to the Bra
! basis function, e.g. iSpinPr is the left and iSpin the right hand spin.
USE m_types
USE m_hs_int_direct
INTEGER :: ispin,jspin !spin indices
INTEGER :: i,j,ii(3),iispin,jjspin,i0
INTEGER :: in
COMPLEX :: th,ts,phase
REAL :: b1(3),b2(3),r2
IMPLICIT NONE
IF (noco%l_noco.AND.isp==2) RETURN !was done already
DO ispin=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
iispin=MIN(ispin,SIZE(smat,1))
DO jspin=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
jjspin=MIN(jspin,SIZE(smat,1))
!$OMP PARALLEL DO SCHEDULE(dynamic) DEFAULT(none) &
!$OMP SHARED(fmpi,lapw,stars,input,cell,vpw) &