Commit 1dd2bbda authored by Gregor Michalicek's avatar Gregor Michalicek

First steps to re-implement the hybrid functionals

parent 0ab711bf
......@@ -65,59 +65,53 @@ CONTAINS
#ifdef CPP_MPI
INCLUDE 'mpif.h'
#endif
! ..
! .. Scalar Arguments ..
! Scalar Arguments
INTEGER,INTENT(IN) :: iter
INTEGER,INTENT(IN) :: eig_id
! ..
!-odim
!+odim
! ..
! .. Local Scalars ..
! Local Scalars
INTEGER jsp,nk,nred,ne_all,ne_found
INTEGER ne,lh0
INTEGER isp,i,j,err
LOGICAL l_wu,l_file,l_real,l_zref
! ..
! .. Local Arrays ..
INTEGER, PARAMETER :: lmaxb=3
! Local Arrays
INTEGER :: ierr(3)
INTEGER :: neigBuffer(kpts%nkpt,input%jspins)
INTEGER, PARAMETER :: lmaxb = 3
REAL, ALLOCATABLE :: bkpt(:)
REAL, ALLOCATABLE :: eig(:)
COMPLEX, ALLOCATABLE :: vs_mmp(:,:,:,:)
TYPE(t_tlmplm) :: td
TYPE(t_usdus) :: ud
TYPE(t_lapw) :: lapw
CLASS(t_Mat),ALLOCATABLE :: zMat
CLASS(t_mat),ALLOCATABLE :: hmat,smat
TYPE(T_mat) :: olap
!
INTEGER nn,n
INTEGER ierr(3)
INTEGER :: neigBuffer(kpts%nkpt,input%jspins)
! .. variables for HF or hybrid functional calculation ..
INTEGER :: comm(kpts%nkpt),irank2(kpts%nkpt),isize2(kpts%nkpt)
TYPE(t_tlmplm) :: td
TYPE(t_usdus) :: ud
TYPE(t_lapw) :: lapw
CLASS(t_mat), ALLOCATABLE :: zMat
CLASS(t_mat), ALLOCATABLE :: hmat,smat
! Variables for HF or hybrid functional calculation
INTEGER :: comm(kpts%nkpt),irank2(kpts%nkpt),isize2(kpts%nkpt)
call ud%init(atoms,DIMENSION%jspd)
ALLOCATE (eig(DIMENSION%neigd),bkpt(3))
l_real=sym%invs.AND..NOT.noco%l_noco
!check if z-reflection trick can be used
! check if z-reflection trick can be used
l_zref=(sym%zrfs.AND.(SUM(ABS(kpts%bk(3,:kpts%nkpt))).LT.1e-9).AND..NOT.noco%l_noco)
IF (mpi%n_size > 1) l_zref = .FALSE.
#ifdef CPP_MPI
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,v)
#endif
!IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber '/),(/iter,v%iter/),&
! RESHAPE((/19,13,5,5/),(/2,2/)))
!---> set up and solve the eigenvalue problem
!---> loop over spins
!---> set up k-point independent t(l'm',lm) matrices
! Set up and solve the eigenvalue problem
! loop over spins
! set up k-point independent t(l'm',lm) matrices
CALL mt_setup(atoms,sym,sphhar,input,noco,enpara,inden,v,mpi,results,DIMENSION,td,ud)
neigBuffer = 0
......@@ -127,43 +121,38 @@ CONTAINS
DO jsp = 1,MERGE(1,input%jspins,noco%l_noco)
k_loop:DO nk = mpi%n_start,kpts%nkpt,mpi%n_stride
!
!---> set up lapw list
! Set up lapw list
CALL lapw%init(input,noco, kpts,atoms,sym,nk,cell,l_zref, mpi)
call timestart("Setup of H&S matrices")
CALL eigen_hssetup(jsp,mpi,DIMENSION,hybrid,enpara,input,vacuum,noco,sym,&
stars,cell,sphhar,atoms,ud,td,v,lapw,l_real,smat,hmat)
stars,cell,sphhar,atoms,ud,td,v,lapw,l_real,smat,hmat)
CALL timestop("Setup of H&S matrices")
IF( hybrid%l_hybrid ) THEN
!write overlap matrix b to direct access file olap
IF(hybrid%l_hybrid) THEN
! Write overlap matrix smat to direct access file olap
print *,"Wrong overlap matrix used, fix this later"
!call write_olap(smat,nrec)
CALL write_olap(smat,(jsp-1)*kpts%nkpt+nk) ! Note: At this moment this only works without MPI parallelization
PRINT *,"TODO"
! STOP "TODO"
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_subvxc ) THEN
PRINT *,"BASIS:", lapw%nv(jsp), atoms%nlotot
IF (hybrid%l_addhf) CALL add_Vnonlocal(nk,hybrid,dimension,kpts,jsp,results,xcpot,hmat)
IF(hybrid%l_subvxc) THEN
!CALL subvxc(lapw,kpts%bk(:,nk),DIMENSION,input,jsp,v%mt(:,0,:,:),atoms,ud,hybrid,enpara%el0,enpara%ello0,&
! sym, atoms%nlotot,kveclo, cell,sphhar, stars, xcpot,mpi,&
! oneD, hamovlp,vx)
! sym,atoms%nlotot,kveclo,cell,sphhar,stars,xcpot,mpi,oneD,hamovlp,vx)
END IF
END IF ! hybrid%l_hybrid
l_wu=.FALSE.
ne_all=DIMENSION%neigd
if (allocated(zmat)) deallocate(zmat)
CALL eigen_diag(hmat,smat,nk,jsp,iter,ne_all,eig,zMat)
DEALLOCATE(hmat,smat)
!
!---> output results
!
! Output results
CALL timestart("EV output")
#if defined(CPP_MPI)
!Collect number of all eigenvalues
! Collect number of all eigenvalues
ne_found=ne_all
CALL MPI_ALLREDUCE(ne_found,ne_all,1,MPI_INTEGER,MPI_SUM, mpi%sub_comm,ierr)
ne_all=MIN(DIMENSION%neigd,ne_all)
......@@ -172,12 +161,12 @@ CONTAINS
#endif
IF (.NOT.l_real) THEN
zMat%data_c(:lapw%nmat,:ne_found) = CONJG(zMat%data_c(:lapw%nmat,:ne_found))
ENDIF
END IF
CALL write_eig(eig_id, nk,jsp,ne_found,ne_all,&
eig(:ne_found),n_start=mpi%n_size,n_end=mpi%n_rank,zmat=zMat)
neigBuffer(nk,jsp) = ne_found
#if defined(CPP_MPI)
!RMA synchronization
! RMA synchronization
CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
#endif
CALL timestop("EV output")
......
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_add_vnonlocal
USE m_judft
USE m_judft
! c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c
! This module is the driver routine for the calculation of the Hartree c
! Fock exchange term by using the mixed basis set. c
......@@ -34,108 +40,95 @@ MODULE m_add_vnonlocal
! c
! M.Betzinger (09/07) c
! c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c
CONTAINS
CONTAINS
SUBROUTINE add_vnonlocal(nk,hybrid,dimension, kpts,jsp,results,xcpot,hamovlp)
SUBROUTINE add_vnonlocal(nk,hybrid,dimension,kpts,jsp,results,xcpot,hmat)
USE m_symm_hf ,ONLY: symm_hf
USE m_util ,ONLY: intgrf,intgrf_init
USE m_symm_hf, ONLY: symm_hf
USE m_util, ONLY: intgrf,intgrf_init
USE m_exchange_valence_hf
USE m_exchange_core
USE m_symmetrizeh
USE m_wrapper
USE m_hsefunctional ,ONLY: exchange_vccvHSE,exchange_ccccHSE
USE m_hsefunctional, ONLY: exchange_vccvHSE,exchange_ccccHSE
USE m_types
USE m_io_hybrid
IMPLICIT NONE
TYPE(t_results),INTENT(INOUT) :: results
TYPE(t_xcpot_inbuild),INTENT(IN) :: xcpot
TYPE(t_dimension),INTENT(IN) :: dimension
TYPE(t_hybrid),INTENT(INOUT) :: hybrid
TYPE(t_kpts),INTENT(IN) :: kpts
INTEGER,INTENT(IN) :: jsp
INTEGER,INTENT(IN) :: nk
TYPE(t_hamovlp),INTENT(INOUT)::hamovlp
! - local scalars -
INTEGER :: ic,n,nn
INTEGER :: iband
IMPLICIT NONE
REAL :: a_ex
TYPE(t_results), INTENT(INOUT) :: results
CLASS(t_xcpot), INTENT(IN) :: xcpot
TYPE(t_dimension), INTENT(IN) :: dimension
TYPE(t_hybrid), INTENT(INOUT) :: hybrid
TYPE(t_kpts), INTENT(IN) :: kpts
TYPE(t_mat), INTENT(INOUT) :: hmat
INTEGER, INTENT(IN) :: jsp
INTEGER, INTENT(IN) :: nk
! local scalars
INTEGER :: n,nn,iband
REAL :: a_ex
TYPE(t_mat) :: olap,tmp,v_x,z
COMPLEX :: exch(dimension%neigd,dimension%neigd)
!
COMPLEX :: exch(dimension%neigd,dimension%neigd)
! initialize weighting factor for HF exchange part
!
a_ex=xcpot%get_exchange_weight()
a_ex=xcpot%get_exchange_weight()
v_x%l_real=hamovlp%l_real
v_x%l_real=hmat%l_real
v_x%matsize1=dimension%nbasfcn
CALL read_v_x(v_x,kpts%nkpt*(jsp-1) + nk)
! add non-local x-potential to the hamiltonian a (in packed storage)
ic=1
DO n=1,v_x%matsize1
DO nn=1,n
IF (hamovlp%l_real) THEN
hamovlp%a_r(ic) = hamovlp%a_r(ic) - a_ex*v_x%data_r(n,nn)
ELSE
hamovlp%a_c(ic) = hamovlp%a_c(ic) - a_ex*v_x%data_c(n,nn)
ENDIF
ic=ic+1
ENDDO
END DO
! calculate HF energy
IF( hybrid%l_calhf ) THEN
WRITE(6,'(A)') new_line('n')//new_line('n')//' ### '// ' diagonal HF exchange elements (eV) ###'
CALL read_v_x(v_x,kpts%nkpt*(jsp-1) + nk)
! add non-local x-potential to the hamiltonian hmat
DO n=1,v_x%matsize1
DO nn=1,n
IF (hmat%l_real) THEN
hmat%data_r(n,nn) = hmat%data_r(n,nn) - a_ex*v_x%data_r(n,nn)
ELSE
hmat%data_c(n,nn) = hmat%data_c(n,nn) - a_ex*v_x%data_c(n,nn)
ENDIF
END DO
END DO
! calculate HF energy
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') // ' k-point '// 'band tail pole total(valence+core)'
END IF
WRITE(6,'(A)') new_line('n') // ' k-point '// 'band tail pole total(valence+core)'
END IF
! read in lower triangle part of overlap matrix from direct acces file olap
call olap%alloc(hamovlp%l_real,dimension%nbasfcn)
call read_olap(olap, kpts%nkpt*(jsp-1) + nk)
if (.not.olap%l_real) olap%data_c=conjg(olap%data_c)
call z%alloc(olap%l_real,dimension%nbasfcn,dimension%neigd)
! read in lower triangle part of overlap matrix from direct acces file olap
CALL olap%alloc(hmat%l_real,dimension%nbasfcn)
CALL read_olap(olap,kpts%nkpt*(jsp-1) + nk)
IF (.NOT.olap%l_real) olap%data_c=conjg(olap%data_c)
call read_z(z,nk) !what about spin?
CALL z%alloc(olap%l_real,dimension%nbasfcn,dimension%neigd)
CALL read_z(z,nk) !what about spin?
WRITE(*,*) 'What about spin (in add_Vnonlocal)?'
! calculate exchange contribution of current k-point nk to total energy (te_hfex)
! in the case of a spin-unpolarized calculation the factor 2 is added in eigen.F90
if (.not.v_x%l_real) v_x%data_c=conjg(v_x%data_c)
exch = 0
z%matsize1=MIN(z%matsize1,v_x%matsize2)
call v_x%multiply(z,tmp)
DO iband = 1,hybrid%nbands(nk)
if (z%l_real) THEN
exch(iband,iband) = dot_product(z%data_r(:z%matsize1,iband),tmp%data_r(:,iband))
else
exch(iband,iband) = dot_product(z%data_r(:z%matsize1,iband),tmp%data_r(:,iband))
endif
IF( iband .le. hybrid%nobd(nk) ) THEN
results%te_hfex%valence = results%te_hfex%valence -a_ex*results%w_iks(iband,nk,jsp)*exch(iband,iband)
END IF
IF(hybrid%l_calhf) THEN
WRITE(6, '( '' ('',F5.3,'','',F5.3,'','',F5.3,'')'',I4,4X,3F15.5)')&
& kpts%bkf(:,nk),iband, (REAL(exch(iband,iband))-hybrid%div_vv(iband,nk,jsp))*(-27.211608),&
& hybrid%div_vv(iband,nk,jsp)*(-27.211608),REAL(exch(iband,iband))*(-27.211608)
END IF
END DO
END SUBROUTINE add_vnonlocal
! calculate exchange contribution of current k-point nk to total energy (te_hfex)
! in the case of a spin-unpolarized calculation the factor 2 is added in eigen.F90
IF (.NOT.v_x%l_real) v_x%data_c=conjg(v_x%data_c)
exch = 0
z%matsize1=MIN(z%matsize1,v_x%matsize2)
CALL v_x%multiply(z,tmp)
DO iband = 1, hybrid%nbands(nk)
IF (z%l_real) THEN
exch(iband,iband) = dot_product(z%data_r(:z%matsize1,iband),tmp%data_r(:,iband))
ELSE
exch(iband,iband) = dot_product(z%data_r(:z%matsize1,iband),tmp%data_r(:,iband))
END IF
IF(iband.LE.hybrid%nobd(nk)) THEN
results%te_hfex%valence = results%te_hfex%valence -a_ex*results%w_iks(iband,nk,jsp)*exch(iband,iband)
END IF
IF(hybrid%l_calhf) THEN
WRITE(6, '( '' ('',F5.3,'','',F5.3,'','',F5.3,'')'',I4,4X,3F15.5)')&
kpts%bkf(:,nk),iband, (REAL(exch(iband,iband))-hybrid%div_vv(iband,nk,jsp))*(-27.211608),&
hybrid%div_vv(iband,nk,jsp)*(-27.211608),REAL(exch(iband,iband))*(-27.211608)
END IF
END DO
END SUBROUTINE add_vnonlocal
END MODULE m_add_vnonlocal
END MODULE m_add_vnonlocal
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_subvxc
CONTAINS
SUBROUTINE subvxc(lapw,bk, DIMENSION,input,jsp,vr0, atoms,usdus, hybrid, el,ello,sym,&
cell, sphhar, stars,xcpot,mpi,oneD,hamovlp,vx)
USE m_intgr, ONLY : intgr3
USE m_constants
USE m_gaunt, ONLY : gaunt1
USE m_wrapper
USE m_loddop
USE m_radflo
USE m_radfun
USE m_abcof3
USE m_types
IMPLICIT NONE
TYPE(t_xcpot_inbuild),INTENT(IN) :: xcpot
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_dimension),INTENT(IN) :: DIMENSION
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_hybrid),INTENT(INOUT) :: hybrid
TYPE(t_input),INTENT(IN) :: input
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_lapw),INTENT(IN) :: lapw
TYPE(t_usdus),INTENT(INOUT) :: usdus
TYPE(t_potden),INTENT(IN) :: vx
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jsp
! .. Array Arguments ..
REAL, INTENT (IN) :: vr0(atoms%jmtd,atoms%ntype,DIMENSION%jspd) ! just for radial functions
REAL, INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd)
REAL, INTENT (IN) :: ello(atoms%nlod,atoms%ntype,DIMENSION%jspd)
REAL, INTENT (IN) :: bk(3)
TYPE(t_hamovlp),INTENT(INOUT)::hamovlp
! .. Local Scalars ..
INTEGER :: ic,indx,m,ig1,ig2
INTEGER :: nlharm,nnbas,typsym,lm
INTEGER :: noded,nodeu
INTEGER :: nbasf0
INTEGER :: i,j,l,ll,l1,l2 ,m1,m2 ,j1,j2
INTEGER :: ok,p1,p2,lh,mh,pp1,pp2
INTEGER :: igrid,itype,ilharm,istar
INTEGER :: ineq,iatom,ilo,ilop,ieq,icentry
INTEGER :: ikvecat,ikvecprevat,invsfct,ikvec,ikvecp
INTEGER :: lp,mp,pp
REAL :: a_ex
REAL :: wronk
COMPLEX :: rc,rr
! .. Local Arrays ..
INTEGER :: gg(3)
INTEGER :: pointer_lo(atoms%nlod,atoms%ntype)
REAL :: integ(0:sphhar%nlhd,hybrid%maxindx,0:atoms%lmaxd,hybrid%maxindx,0:atoms%lmaxd)
REAL :: grid(atoms%jmtd)
REAL :: vr(atoms%jmtd,0:sphhar%nlhd)
REAL :: f(atoms%jmtd,2,0:atoms%lmaxd),g(atoms%jmtd,2,0:atoms%lmaxd)
REAL :: flo(atoms%jmtd,2,atoms%nlod)
REAL :: uuilon(atoms%nlod,atoms%ntype),duilon(atoms%nlod,atoms%ntype)
REAL :: ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype)
REAL :: bas1(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype),&
bas2(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype)
COMPLEX :: vpw(stars%ng3)
COMPLEX :: vxc(hamovlp%matsize)
COMPLEX :: vrmat(hybrid%maxlmindx,hybrid%maxlmindx)
COMPLEX :: carr(hybrid%maxlmindx,DIMENSION%nvd),carr1(DIMENSION%nvd,DIMENSION%nvd)
COMPLEX ,ALLOCATABLE :: ahlp(:,:,:),bhlp(:,:,:)
COMPLEX, ALLOCATABLE :: bascof(:,:,:)
COMPLEX :: bascof_lo(3,-atoms%llod:atoms%llod,4*atoms%llod+2,atoms%nlod, atoms%nat)
CALL timestart("subvxc")
vxc=0
! calculate radial functions
hybrid%nindx = 2
DO itype = 1,atoms%ntype
!
!---> generate the radial basis-functions for each l
!
WRITE(6,'(a,i3,a)') new_LINE('n')//new_LINE('n')//' wavefunction parameters for atom type',itype,':'
WRITE(6,'(31x,a,32x,a)') 'radial function','energy derivative'
WRITE(6,'(a)') ' l energy value '//&
'derivative nodes value derivative nodes norm wronskian'
DO l = 0,atoms%lmax(itype)
CALL radfun(l,itype,jsp,el(l,itype,jsp),vr0(1,itype,jsp),atoms,&
f(1,1,l),g(1,1,l),usdus,&
nodeu,noded,wronk)
WRITE (6,FMT=8010) l,el(l,itype,jsp),usdus%us(l,itype,jsp),&
usdus%dus(l,itype,jsp),nodeu,usdus%uds(l,itype,jsp),usdus%duds(l,itype,jsp),noded,&
usdus%ddn(l,itype,jsp),wronk
END DO
! 8000 FORMAT (1x,/,/,' wavefunction parameters for atom type',i3,':',
! + /,t32,'radial function',t79,'energy derivative',/,t3,
! + 'l',t8,'energy',t26,'value',t39,'derivative',t53,
! + 'nodes',t68,'value',t81,'derivative',t95,'nodes',t107,
! + 'norm',t119,'wronskian')
8010 FORMAT (i3,f10.5,2 (5x,1p,2e16.7,i5),1p,2e16.7)
bas1(:,1,:,itype)=f(:,1,:)
bas1(:,2,:,itype)=g(:,1,:)
bas2(:,1,:,itype)=f(:,2,:)
bas2(:,2,:,itype)=g(:,2,:)
!
!---> generate the extra radial basis-functions for the local orbitals,
!---> if there are any.
!
IF (atoms%nlo(itype).GE.1) THEN
CALL radflo(&
atoms,itype,jsp,ello(1,1,jsp),vr0(1,itype,jsp),&
f,g,mpi,&
usdus,&
uuilon,duilon,ulouilopn,flo,.TRUE.)
DO i=1,atoms%nlo(itype)
hybrid%nindx(atoms%llo(i,itype),itype) = hybrid%nindx(atoms%llo(i,itype),itype) + 1
pointer_lo(i,itype) = hybrid%nindx(atoms%llo(i,itype),itype)
bas1(:,hybrid%nindx(atoms%llo(i,itype),itype),atoms%llo(i,itype),itype)=&
flo(:,1,i)
bas2(:,hybrid%nindx(atoms%llo(i,itype),itype),atoms%llo(i,itype),itype)=&
flo(:,2,i)
END DO
END IF
END DO
! compute APW coefficients
! calculate bascof
ALLOCATE( ahlp(DIMENSION%nvd,0:DIMENSION%lmd,atoms%nat),bhlp(DIMENSION%nvd,0:DIMENSION%lmd,atoms%nat),stat=ok)
IF( ok .NE. 0 ) STOP 'subvxc: error in allocation of ahlp/bhlp'
CALL abcof3( input,atoms,sym,jsp,cell, bk,lapw,&
usdus,oneD,ahlp,bhlp,bascof_lo)
ALLOCATE( bascof(DIMENSION%nvd,2*(DIMENSION%lmd+1),atoms%nat), stat=ok )
IF( ok .NE. 0 ) STOP 'subvxc: error in allocation of bascof'
bascof = 0
ic = 0
DO itype=1,atoms%ntype
DO ieq=1,atoms%neq(itype)
ic = ic + 1
indx = 0
DO l=0,atoms%lmax(itype)
ll = l*(l+1)
DO M=-l,l
lm=ll+M
DO i=1,2
indx = indx + 1
IF( i .EQ. 1) THEN
bascof(:,indx,ic) = ahlp(:,lm,ic)
ELSE IF( i .EQ. 2 ) THEN
bascof(:,indx,ic) = bhlp(:,lm,ic)
END IF
END DO
END DO
END DO
END DO
END DO
DEALLOCATE( ahlp,bhlp )
! Loop over atom types
iatom = 0
DO itype = 1,atoms%ntype
typsym = atoms%ntypsy( SUM(atoms%neq(:itype-1))+1 )
nlharm = sphhar%nlh(typsym)
! Calculate vxc = vtot - vcoul
DO l=0,nlharm
DO i=1,atoms%jri(itype)
IF(l.EQ.0) THEN
! vr(i,0)= vrtot(i,0,itype)*sfp/rmsh(i,itype) - vrcou(i,0,itype,jsp)
vr(i,0)= vx%mt(i,0,itype,jsp)*sfp_const/atoms%rmsh(i,itype) !
ELSE ! vxc = vtot - vcoul
! vr(i,l)= vrtot(i,l,itype)-vrcou(i,l,itype,jsp)
vr(i,l)= vx%mt(i,l,itype,jsp)
END IF
END DO
END DO
! Calculate MT contribution to vxc matrix elements
! Precompute auxiliary radial integrals
DO ilharm = 0,nlharm
i = 0
DO l1 = 0,atoms%lmax(itype)
DO p1 = 1,2
i = i + 1
j = 0
DO l2 = 0,atoms%lmax(itype)
DO p2 = 1,2
j = j + 1
IF( j .LE. i) THEN
DO igrid = 1,atoms%jri(itype)
grid(igrid)=vr(igrid,ilharm)*(bas1(igrid,p1,l1,itype)*bas1(igrid,p2,l2,itype)+ bas2(igrid,p1,l1,itype)*bas2(igrid,p2,l2,itype) )
END DO
CALL intgr3(grid,atoms%rmsh(:,itype),atoms%dx(itype),atoms%jri(itype),integ(ilharm,p1,l1,p2,l2) ) ! numerical integration
integ(ilharm,p2,l2,p1,l1)=integ(ilharm,p1,l1,p2,l2)
END IF
END DO
END DO
END DO
END DO
END DO
! Calculate muffin tin contribution to vxc matrix
vrmat=0
j1=0
DO l1 = 0,atoms%lmax(itype) ! loop: left basis function
DO m1 = -l1,l1
DO p1 = 1,2
j1 = j1+1
j2 = 0
DO l2 = 0,atoms%lmax(itype) ! loop: right basis function
DO m2 = -l2,l2
DO p2 = 1,2
j2 = j2+1
rr = 0
DO ilharm = 0,nlharm ! loop: lattice harmonics of vxc
l = sphhar%llh(ilharm,typsym)
DO i = 1,sphhar%nmem(ilharm,typsym)
M = sphhar%mlh(i,ilharm,typsym)
rc = sphhar%clnu(i,ilharm,typsym)* gaunt1(l1,l,l2,m1,M,m2,atoms%lmaxd)
rr = rr+integ(ilharm,p1,l1,p2,l2)*rc
END DO
END DO
rc = CMPLX(0,1)**(l2-l1) ! adjusts to a/b/ccof-scaling
vrmat(j1,j2) = rr*rc
END DO
END DO
END DO
END DO
END DO
END DO
nnbas = j1
! ! Project on bascof
DO ineq = 1,atoms%neq(itype)
iatom = iatom+1
carr (:nnbas,:lapw%nv(jsp)) = CONJG(MATMUL(vrmat(:nnbas,:nnbas), TRANSPOSE(bascof(:lapw%nv(jsp),:nnbas,iatom)) ))
carr1(:lapw%nv(jsp),:lapw%nv(jsp)) = MATMUL(bascof(:lapw%nv(jsp),:nnbas,iatom),carr(:nnbas,:lapw%nv(jsp)) )
ic = 0
DO j = 1,lapw%nv(jsp)
! carr(:nnbas) = matmul(vrmat(:nnbas,:nnbas),
! + bascof(j,:nnbas,iatom) )
DO i = 1,j
ic = ic + 1
vxc(ic) = vxc(ic) + carr1(i,j)
! vxc(ic) = vxc(ic) + conjg(dotprod ( bascof(i,:nnbas,iatom),
! + carr(:nnbas) ))
END DO
END DO
END DO
END DO ! End loop over atom types
! ---------------------------------------------------------------
! Calculate plane wave contribution
DO i=1,stars%ng3
vpw(i)= vx%pw(i,jsp)
! vpw(i)=vpwtot(i)-vpwcou(i,jsp)
END DO
! Calculate vxc-matrix, left basis function (ig1)
! right basis function (ig2)
ic = 0
DO ig1=1,lapw%nv(jsp)
DO ig2=1,ig1
ic = ic + 1
gg(1)=lapw%k1(ig1,jsp)-lapw%k1(ig2,jsp)
gg(2)=lapw%k2(ig1,jsp)-lapw%k2(ig2,jsp)
gg(3)=lapw%k3(ig1,jsp)-lapw%k3(ig2,jsp)
istar=stars%ig(gg(1),gg(2),gg(3))
IF(istar.NE.0) THEN
vxc(ic)= vxc(ic) + stars%rgphs(gg(1),gg(2),gg(3))*vpw(istar)
ELSE
IF ( mpi%irank == 0 ) WRITE(6,'(A,/6I5)') 'Warning: Gi-Gj not in any star:',&
lapw%k1(ig1,jsp),lapw%k2(ig1,jsp),lapw%k3(ig1,jsp),&
lapw%k1(ig2,jsp),lapw%k2(ig2,jsp),lapw%k3(ig2,jsp)
ENDIF
ENDDO
ENDDO
!
! -------------------------------------------------------------------
! Calculate local orbital contribution