Commit e498c60c authored by Gregor Michalicek's avatar Gregor Michalicek

added Intraatomic magnetic dipole operator and multipole expansions

Note: The output unit for <T> is hbar/2. ...this has to be changed in the future.
Note2: Switches have to be implemented to toggle the computation of these quantities.
parent f72dc951
......@@ -134,7 +134,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
CALL denCoeffs%init(atoms,sphhar,jsp_start,jsp_end)
! The last entry in denCoeffsOffdiag%init is l_fmpl. It is meant as a switch to a plot of the full magnet.
! density without the atomic sphere approximation for the magnet. density. It is not completely implemented (lo's missing).
CALL denCoeffsOffdiag%init(atoms,noco,sphhar,noco%l_mtnocopot)
CALL denCoeffsOffdiag%init(atoms,noco,sphhar,noco%l_mtnocopot.OR.noco%l_mperp)
CALL force%init1(input,atoms)
CALL orb%init(atoms,noco,jsp_start,jsp_end)
......
......@@ -10,10 +10,13 @@ cdn_mt/abcof3.F90
cdn_mt/abcrot2.f90
cdn_mt/cdnmt.F90
cdn_mt/cdncore.F90
cdn_mt/denMultipoleExp.f90
cdn_mt/magMoms.f90
cdn_mt/magDiMom.f90
cdn_mt/orbMagMoms.f90
cdn_mt/orb_comp2.f90
cdn_mt/radfun.f90
cdn_mt/resMoms.f90
cdn_mt/rhomt.f90
cdn_mt/rhomt21.f90
cdn_mt/rhomtlo.f90
......
......@@ -38,7 +38,7 @@ CONTAINS
TYPE (t_denCoeffsOffdiag), INTENT(IN) :: denCoeffsOffdiag
! ..
! .. Local Scalars ..
INTEGER itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu
INTEGER itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu,llpb
INTEGER ilo,ilop,i
REAL s,wronk,sumlm,qmtt
COMPLEX cs
......@@ -63,11 +63,11 @@ CONTAINS
ENDIF
ENDIF
!$OMP PARALLEL DEFAULT(none) &
!$OMP SHARED(usdus,rho,moments,qmtl) &
!$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
!$OMP SHARED(orb,noco,denCoeffsOffdiag,jspd)&
!$OMP PRIVATE(itype,na,ispin,l,rho21,f,g,nodeu,noded,wronk,i,j,s,qmtllo,qmtt,nd,lh,lp,llp,cs)
! !$OMP PARALLEL DEFAULT(none) &
! !$OMP SHARED(usdus,rho,moments,qmtl) &
! !$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
! !$OMP SHARED(orb,noco,denCoeffsOffdiag,jspd)&
! !$OMP PRIVATE(itype,na,ispin,l,rho21,f,g,nodeu,noded,wronk,i,j,s,qmtllo,qmtt,nd,lh,lp,llp,llpb,cs)
IF (noco%l_mperp) THEN
ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jspd),g(atoms%jmtd,2,0:atoms%lmaxd,jspd) )
ELSE
......@@ -77,7 +77,7 @@ CONTAINS
qmtl = 0
!$OMP DO
! !$OMP DO
DO itype = 1,atoms%ntype
na = 1
DO i = 1, itype - 1
......@@ -88,11 +88,13 @@ CONTAINS
DO l = 0,atoms%lmax(itype)
CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vr(1,itype,ispin),atoms,&
f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
llp = (l* (l+1))/2 + l
DO j = 1,atoms%jri(itype)
s = denCoeffs%uu(l,itype,ispin)*( f(j,1,l,ispin)*f(j,1,l,ispin)+f(j,2,l,ispin)*f(j,2,l,ispin) )&
+ denCoeffs%dd(l,itype,ispin)*( g(j,1,l,ispin)*g(j,1,l,ispin)+g(j,2,l,ispin)*g(j,2,l,ispin) )&
+ 2*denCoeffs%du(l,itype,ispin)*( f(j,1,l,ispin)*g(j,1,l,ispin)+f(j,2,l,ispin)*g(j,2,l,ispin) )
rho(j,0,itype,ispin) = rho(j,0,itype,ispin)+ s/(atoms%neq(itype)*sfp_const)
moments%rhoLRes(j,0,llp,itype,ispin) = moments%rhoLRes(j,0,llp,itype,ispin)+ s/(atoms%neq(itype)*sfp_const)
ENDDO
ENDDO
......@@ -110,7 +112,7 @@ CONTAINS
denCoeffs%cclo(1,1,itype,ispin),denCoeffs%acnmt(0,1,1,itype,ispin),&
denCoeffs%bcnmt(0,1,1,itype,ispin),denCoeffs%ccnmt(1,1,1,itype,ispin),&
f(1,1,0,ispin),g(1,1,0,ispin),&
rho(:,0:,itype,ispin),qmtllo)
rho(:,0:,itype,ispin),qmtllo,moments%rhoLRes(:,0:,:,itype,ispin))
!---> l-decomposed density for each atom type
......@@ -146,6 +148,7 @@ CONTAINS
+ denCoeffs%dunmt(llp,lh,itype,ispin)*(g(j,1,l,ispin)*f(j,1,lp,ispin)&
+ g(j,2,l,ispin)*f(j,2,lp,ispin) )
rho(j,lh,itype,ispin) = rho(j,lh,itype,ispin)+ s/atoms%neq(itype)
moments%rhoLRes(j,lh,llp,itype,ispin) = moments%rhoLRes(j,lh,llp,itype,ispin)+ s/atoms%neq(itype)
ENDDO
ENDDO
ENDDO
......@@ -183,6 +186,7 @@ CONTAINS
!---> calculate off-diagonal part of the density matrix
!---> spherical component
DO l = 0,atoms%lmax(itype)
llp = (l* (l+1))/2 + l
DO j = 1,atoms%jri(itype)
cs = denCoeffsOffdiag%uu21(l,itype)*( f(j,1,l,2)*f(j,1,l,1) +f(j,2,l,2)*f(j,2,l,1) )&
+ denCoeffsOffdiag%ud21(l,itype)*( f(j,1,l,2)*g(j,1,l,1) +f(j,2,l,2)*g(j,2,l,1) )&
......@@ -192,6 +196,8 @@ CONTAINS
rho21=CONJG(cs)/(atoms%neq(itype)*sfp_const)
rho(j,0,itype,3)=rho(j,0,itype,3)+REAL(rho21)
rho(j,0,itype,4)=rho(j,0,itype,4)+imag(rho21)
moments%rhoLRes(j,0,llp,itype,3) = moments%rhoLRes(j,0,llp,itype,3)+ REAL(conjg(cs)/(atoms%neq(itype)*sfp_const))
moments%rhoLRes(j,0,llp,itype,4) = moments%rhoLRes(j,0,llp,itype,4)+ AIMAG(conjg(cs)/(atoms%neq(itype)*sfp_const))
ENDDO
ENDDO
......@@ -201,6 +207,7 @@ CONTAINS
DO l = 0,atoms%lmax(itype)
DO lp = 0,atoms%lmax(itype)
llp = lp*(atoms%lmax(itype)+1)+l+1
llpb = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
DO j = 1,atoms%jri(itype)
cs = denCoeffsOffdiag%uunmt21(llp,lh,itype)*(f(j,1,lp,2)*f(j,1,l,1)&
+ f(j,2,lp,2)*f(j,2,l,1) )+ denCoeffsOffdiag%udnmt21(llp,lh,itype)*(f(j,1,lp,2)*g(j,1,l,1)&
......@@ -211,6 +218,8 @@ CONTAINS
rho21=CONJG(cs)/atoms%neq(itype)
rho(j,lh,itype,3)=rho(j,lh,itype,3)+REAL(rho21)
rho(j,lh,itype,4)=rho(j,lh,itype,4)+imag(rho21)
moments%rhoLRes(j,lh,llpb,itype,3)= moments%rhoLRes(j,lh,llpb,itype,3) + REAL(conjg(cs)/atoms%neq(itype))
moments%rhoLRes(j,lh,llpb,itype,4)= moments%rhoLRes(j,lh,llpb,itype,4) + AIMAG(conjg(cs)/atoms%neq(itype))
ENDDO
ENDDO
ENDDO
......@@ -220,9 +229,9 @@ CONTAINS
ENDIF ! noco%l_mperp
ENDDO ! end of loop over atom types
!$OMP END DO
! !$OMP END DO
DEALLOCATE ( f,g)
!$OMP END PARALLEL
! !$OMP END PARALLEL
WRITE (6,FMT=8000)
8000 FORMAT (/,5x,'l-like charge',/,t6,'atom',t15,'s',t24,'p',&
......
MODULE m_denMultipoleExp
IMPLICIT NONE
CONTAINS
SUBROUTINE denMultipoleExp(input, mpi, atoms, sphhar, stars, sym, cell, oneD, den)
USE m_types
USE m_constants
USE m_mpmom
TYPE(t_input), INTENT(IN) :: input
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_oneD), INTENT(IN) :: oneD
TYPE(t_potden), INTENT(IN) :: den
type(t_potden) :: workDen
COMPLEX :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
IF(input%jspins == 2) THEN
IF(mpi%irank.EQ.0) THEN
WRITE(6,*) 'Multipole expansion for spin-up density:'
WRITE(6,*) '======================================='
END IF
qlm = CMPLX(0.0,0.0)
workDen = den
CALL mpmom(input,mpi,atoms,sphhar,stars,sym,cell,oneD,workDen%pw(1:,1),workDen%mt(:,0:,1:,1),POTDEN_TYPE_DEN,qlm,.FALSE.)
IF(mpi%irank.EQ.0) THEN
WRITE(6,*) '======================================='
END IF
IF(mpi%irank.EQ.0) THEN
WRITE(6,*) 'Multipole expansion for spin-down density:'
WRITE(6,*) '======================================='
END IF
qlm = CMPLX(0.0,0.0)
CALL mpmom(input,mpi,atoms,sphhar,stars,sym,cell,oneD,workDen%pw(1:,2),workDen%mt(:,0:,1:,2),POTDEN_TYPE_DEN,qlm,.FALSE.)
IF(mpi%irank.EQ.0) THEN
WRITE(6,*) '======================================='
END IF
END IF
IF(mpi%irank.EQ.0) THEN
WRITE(6,*) 'Multipole expansion for charge density:'
WRITE(6,*) '======================================='
END IF
qlm = CMPLX(0.0,0.0)
workDen = den
IF(input%jspins == 2) CALL workDen%SpinsToChargeAndMagnetisation()
CALL mpmom(input,mpi,atoms,sphhar,stars,sym,cell,oneD,workDen%pw(1:,1),workDen%mt(:,0:,1:,1),POTDEN_TYPE_DEN,qlm,.FALSE.)
IF(mpi%irank.EQ.0) THEN
WRITE(6,*) '======================================='
END IF
IF(input%jspins == 2) THEN
IF(mpi%irank.EQ.0) THEN
WRITE(6,*) 'Multipole expansion for magnetization density:'
WRITE(6,*) '======================================='
END IF
qlm = CMPLX(0.0,0.0)
CALL mpmom(input,mpi,atoms,sphhar,stars,sym,cell,oneD,workDen%pw(1:,2),workDen%mt(:,0:,1:,2),POTDEN_TYPE_DEN,qlm,.FALSE.)
IF(mpi%irank.EQ.0) THEN
WRITE(6,*) '======================================='
END IF
END IF
END SUBROUTINE denMultipoleExp
END MODULE m_denMultipoleExp
This diff is collapsed.
! Copyright (c) 2018 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_resMoms
CONTAINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This subroutine calculates and writes out intraatomic electric and magnetic dipole
! moments resolved with respect to their orbital (angular momentum) origins.
!
! GM'2018
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE resMoms(input,atoms,sphhar,noco,den,rhoLRes)
USE m_constants
USE m_types
USE m_juDFT
USE m_magDiMom
IMPLICIT NONE
TYPE(t_input), INTENT(IN) :: input
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), INTENT(IN) :: den
REAL, INTENT(IN) :: rhoLRes(:,0:,0:,:,:)
REAL, ALLOCATABLE :: rhoTemp(:,:,:,:)
REAL :: t_op(3,atoms%ntype), elecDip(3,atoms%ntype)
REAL :: res_T_op(3,atoms%ntype,0:(atoms%lmaxd*(atoms%lmaxd+1))/2+atoms%lmaxd)
REAL :: resElecDip(3,atoms%ntype,0:(atoms%lmaxd*(atoms%lmaxd+1))/2+atoms%lmaxd)
INTEGER :: iType, l, lp, llp
IF(input%jspins.EQ.1) RETURN
IF(.NOT.noco%l_noco) RETURN
t_op = 0.0
res_T_op = 0.0
elecDip = 0.0
resElecDip = 0.0
ALLOCATE(rhoTemp(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,4))
rhoTemp = 0.0
rhoTemp(:,:,:,1) = den%mt(:,:,:,1)
rhoTemp(:,:,:,2) = den%mt(:,:,:,2)
IF (noco%l_mperp) THEN
rhoTemp(:,:,:,3) = den%mt(:,:,:,3)
rhoTemp(:,:,:,4) = den%mt(:,:,:,4)
! WRITE(5000,'(f15.8)') den%mt(:,:,:,3)
! WRITE(5000,'(f15.8)') den%mt(:,:,:,4)
END IF
CALL magDiMom(input,atoms,sphhar,noco,noco%l_mperp,rhoTemp,t_op,elecDip)
DO l = 0, atoms%lmaxd
DO lp = 0, l
llp = (l* (l+1))/2 + lp
rhoTemp = 0.0
rhoTemp(:,:,:,1) = rhoLRes(:,:,llp,:,1)
rhoTemp(:,:,:,2) = rhoLRes(:,:,llp,:,2)
rhoTemp(:,:,:,3) = rhoLRes(:,:,llp,:,3)
rhoTemp(:,:,:,4) = rhoLRes(:,:,llp,:,4)
CALL magDiMom(input,atoms,sphhar,noco,noco%l_mperp,rhoTemp,res_T_op(:,:,llp),resElecDip(:,:,llp))
END DO
END DO
DO iType = 1, atoms%ntype
WRITE(6,*) 'Intraatomic electric and magnetic dipole moments for atom type ', iType,':'
WRITE(6,'(a)') ' lowL largeL p_x p_y p_z t_x t_y t_z'
WRITE(6,'(a,6f15.8)') 'Overall: ', elecDip(:,iType), t_op(:,iType)
DO l = 0, atoms%lmax(iType)
DO lp = 0, l
llp = (l* (l+1))/2 + lp
IF(ALL(ABS(res_T_op(:,iType,llp)).LT.1.0e-8).AND.&
ALL(ABS(resElecDip(:,iType,llp)).LT.1.0e-8)) CYCLE
WRITE(6,'(a,2i6,6f15.8)') ' ', lp, l, resElecDip(:,iType,llp),res_T_op(:,iType,llp)
END DO
END DO
END DO
END SUBROUTINE resMoms
END MODULE m_resMoms
......@@ -13,7 +13,7 @@ MODULE m_rhosphnlo
!***********************************************************************
CONTAINS
SUBROUTINE rhosphnlo(itype,atoms,sphhar, uloulopn,dulon,uulon,&
ello,vr, aclo,bclo,cclo,acnmt,bcnmt,ccnmt,f,g, rho,qmtllo)
ello,vr, aclo,bclo,cclo,acnmt,bcnmt,ccnmt,f,g, rho,qmtllo,rhoLRes)
USE m_constants, ONLY : c_light,sfp_const
USE m_radsra
......@@ -36,11 +36,12 @@ CONTAINS
REAL, INTENT (IN) :: f(atoms%jmtd,2,0:atoms%lmaxd),g(atoms%jmtd,2,0:atoms%lmaxd)
REAL, INTENT (INOUT) :: qmtllo(0:atoms%lmaxd)
REAL, INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd)
REAL, INTENT (INOUT) :: rhoLRes(atoms%jmtd,0:sphhar%nlhd,0:(atoms%lmaxd*(atoms%lmaxd+1))/2+atoms%lmaxd)
! ..
! .. Local Scalars ..
REAL dsdum,usdum ,c_1,c_2
INTEGER j,l,lh,lo,lop,lp,nodedum
REAL dus,ddn,c
INTEGER j,l,lh,lo,lop,lp,nodedum,llp
REAL dus,ddn,c,temp
! ..
! .. Local Arrays ..
REAL, ALLOCATABLE :: flo(:,:,:),glo(:,:)
......@@ -90,16 +91,21 @@ CONTAINS
DO lo = 1,atoms%nlo(itype)
l = atoms%llo(lo,itype)
llp = (l* (l+1))/2 + l
DO j = 1,atoms%jri(itype)
rho(j,0) = rho(j,0) + c_2 *&
(aclo(lo) * ( f(j,1,l)*flo(j,1,lo) +f(j,2,l)*flo(j,2,lo) ) +&
bclo(lo) * ( g(j,1,l)*flo(j,1,lo) +g(j,2,l)*flo(j,2,lo) ) )
temp = c_2 *&
(aclo(lo) * ( f(j,1,l)*flo(j,1,lo) +f(j,2,l)*flo(j,2,lo) ) +&
bclo(lo) * ( g(j,1,l)*flo(j,1,lo) +g(j,2,l)*flo(j,2,lo) ) )
rho(j,0) = rho(j,0) + temp
rhoLRes(j,0,llp) = rhoLRes(j,0,llp) + temp
END DO
DO lop = 1,atoms%nlo(itype)
IF (atoms%llo(lop,itype).EQ.l) THEN
DO j = 1,atoms%jri(itype)
rho(j,0) = rho(j,0) + c_2 * cclo(lop,lo) *&
temp = c_2 * cclo(lop,lo) *&
( flo(j,1,lop)*flo(j,1,lo) +flo(j,2,lop)*flo(j,2,lo) )
rho(j,0) = rho(j,0) + temp
rhoLRes(j,0,llp) = rhoLRes(j,0,llp) + temp
END DO
END IF
END DO
......@@ -111,18 +117,27 @@ CONTAINS
DO lh = 1,sphhar%nlh(atoms%ntypsy(atoms%nat))
DO lp = 0,atoms%lmax(itype)
DO lo = 1,atoms%nlo(itype)
l = atoms%llo(lo,itype)
llp = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
DO j = 1,atoms%jri(itype)
rho(j,lh) = rho(j,lh) + c_1 * (&
temp = c_1 * (&
acnmt(lp,lo,lh) * (f(j,1,lp)*flo(j,1,lo) +f(j,2,lp)*flo(j,2,lo) ) +&
bcnmt(lp,lo,lh) * (g(j,1,lp)*flo(j,1,lo) +g(j,2,lp)*flo(j,2,lo) ) )
rho(j,lh) = rho(j,lh) + temp
rhoLRes(j,lh,llp) = rhoLRes(j,lh,llp) + temp
END DO
END DO
END DO
DO lo = 1,atoms%nlo(itype)
l = atoms%llo(lo,itype)
DO lop = 1,atoms%nlo(itype)
lp = atoms%llo(lop,itype)
llp = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
DO j = 1,atoms%jri(itype)
rho(j,lh) = rho(j,lh) + c_1 * ccnmt(lop,lo,lh) *&
temp = c_1 * ccnmt(lop,lo,lh) *&
( flo(j,1,lop)*flo(j,1,lo) +flo(j,2,lop)*flo(j,2,lo) )
rho(j,lh) = rho(j,lh) + temp
rhoLRes(j,lh,llp) = rhoLRes(j,lh,llp) + temp
END DO
END DO
END DO
......
......@@ -1690,7 +1690,11 @@ MODULE m_cdnpot_io_hdf
dimsInt(:4)=(/jmtd,nlhd+1,ntype,input%jspins/)
CALL h5dopen_f(groupID, 'fr', frSetID, hdfError)
CALL io_write_real4(frSetID,(/1,1,1,1/),dimsInt(:4),den%mt)
! Note: The last dimension of den%mt (input%jspins) is temporary to
! avoid segmentation faults if l_mperp is set to true but there
! already is a data set with l_mperp=false. At the moment this is ok
! since the offdiagonal parts are never read.
CALL io_write_real4(frSetID,(/1,1,1,1/),dimsInt(:4),den%mt(:,0:,:,:input%jspins))
CALL h5dclose_f(frSetID, hdfError)
dimsInt(:3)=(/2,ng3,input%jspins/)
......@@ -1861,7 +1865,12 @@ MODULE m_cdnpot_io_hdf
CALL h5screate_simple_f(4,dims(:4),frSpaceID,hdfError)
CALL h5dcreate_f(groupID, "fr", H5T_NATIVE_DOUBLE, frSpaceID, frSetID, hdfError)
CALL h5sclose_f(frSpaceID,hdfError)
CALL io_write_real4(frSetID,(/1,1,1,1/),dimsInt(:4),den%mt)
! Note: The last dimension of den%mt (input%jspins) is temporary to
! avoid segmentation faults if l_mperp is set to true but there
! already is a data set with l_mperp=false. At the moment this is ok
! since the offdiagonal parts are never read.
CALL io_write_real4(frSetID,(/1,1,1,1/),dimsInt(:4),den%mt(:,0:,:,:input%jspins))
! CALL io_write_real4(frSetID,(/1,1,1,1/),dimsInt(:4),den%mt)
CALL h5dclose_f(frSetID, hdfError)
dims(:3)=(/2,ng3,input%jspins/)
......
......@@ -30,12 +30,14 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
USE m_xmlOutput
USE m_magMoms
USE m_orbMagMoms
USE m_resMoms
USE m_cdncore
USE m_doswrite
USE m_Ekwritesl
USE m_banddos_io
USE m_metagga
USE m_unfold_band_kpts
USE m_denMultipoleExp
#ifdef CPP_MPI
USE m_mpi_bc_potden
#endif
......@@ -96,7 +98,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
CALL regCharges%init(input,atoms)
CALL dos%init(input,atoms,dimension,kpts,vacuum)
CALL moments%init(input,atoms)
CALL moments%init(mpi,input,sphhar,atoms)
CALL mcd%init1(banddos,dimension,input,atoms,kpts)
CALL slab%init(banddos,dimension,atoms,cell,input,kpts)
CALL orbcomp%init(input,banddos,dimension,atoms,kpts)
......@@ -175,6 +177,11 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
call core_den%subPotDen(outDen, val_den)
CALL timestop("cdngen: cdncore")
IF(.TRUE.) CALL denMultipoleExp(input, mpi, atoms, sphhar, stars, sym, cell, oneD, outDen)
IF(mpi%irank.EQ.0) THEN
IF(.TRUE.) CALL resMoms(input,atoms,sphhar,noco,outDen,moments%rhoLRes) ! There should be a switch in the inp file for this
END IF
CALL enpara%calcOutParams(input,atoms,vacuum,regCharges)
IF (mpi%irank == 0) THEN
......
......@@ -200,10 +200,10 @@ CONTAINS
dimension%neigd2 = dimension%neigd
IF (noco%l_soc) dimension%neigd2 = dimension%neigd*2
IF (hybrid%l_hybrid) THEN
!HF
!$ num_threads = omp_get_max_threads()
!$ call omp_set_num_threads(1)
IF (hybrid%l_hybrid) THEN
SELECT TYPE(xcpot)
TYPE IS(t_xcpot_inbuild)
CALL calc_hybrid(eig_id,hybrid,kpts,atoms,input,DIMENSION,mpi,noco,&
......@@ -213,6 +213,7 @@ CONTAINS
call mixing_history_reset(mpi)
iter = 0
END IF
!$ call omp_set_num_threads(num_threads)
ENDIF
!RDMFT
IF(input%l_rdmft) THEN
......@@ -220,8 +221,6 @@ CONTAINS
END IF
CALL reset_eig(eig_id,noco%l_soc) ! This has to be placed after the calc_hybrid call but before eigen
!$ call omp_set_num_threads(num_threads)
!#endif
!!$ DO pc = 1, wann%nparampts
......
......@@ -29,6 +29,7 @@ math/fft3d.f90
math/fft_interface.F90
math/SphBessel.f90
math/DoubleFactorial.f90
math/LattHarmsSphHarmsConv.f90
math/ExpSave.f90
math/intgr.F90
math/ylm4.F90
......
MODULE m_lattHarmsSphHarmsConv
CONTAINS
SUBROUTINE lattHarmsRepToSphHarms(atoms, lattHarms, iType, funcLattHarms, funcSphHarms)
USE m_types
IMPLICIT NONE
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: lattHarms
INTEGER, INTENT(IN) :: iType
REAL, INTENT(IN) :: funcLattHarms(:,0:) ! (iR,iLH)
COMPLEX, INTENT(INOUT) :: funcSphHarms(:,:) ! (iR,lm)
INTEGER :: iAtom, iLH, ns, l, iM, m, lm, iR
iAtom = SUM(atoms%neq(:iType-1)) + 1
ns = atoms%ntypsy(iAtom)
funcSphHarms = CMPLX(0.0,0.0)
DO iLH = 0, lattHarms%nlh(ns)
l = lattHarms%llh(iLH,ns)
DO iM = 1, lattHarms%nmem(iLH,ns)
m = lattHarms%mlh(iM,iLH,ns)
lm = l*(l+1) + m + 1
DO iR = 1, atoms%jri(iType)
funcSphHarms(iR,lm) = funcSphHarms(iR,lm) + funcLattHarms(iR,iLH) * lattHarms%clnu(iM,iLH,ns)
END DO
END DO
END DO
END SUBROUTINE lattHarmsRepToSphHarms
SUBROUTINE sphHarmsRepToLattHarms()
IMPLICIT NONE
END SUBROUTINE sphHarmsRepToLattHarms
END MODULE m_lattHarmsSphHarmsConv
......@@ -253,7 +253,8 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
CALL regCharges%init(input,atoms)
CALL dos%init(input,atoms,dimension,kpts,vacuum)
CALL moments%init(input,atoms)
! CALL moments%init(input,atoms)
CALL moments%init(mpi,input,sphhar,atoms)
CALL overallDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
CALL overallVCoul%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTCOUL)
IF (ALLOCATED(vTot%pw_w)) DEALLOCATE (vTot%pw_w)
......
......@@ -105,6 +105,8 @@ PRIVATE
REAL, ALLOCATABLE :: stdn(:,:)
REAL, ALLOCATABLE :: svdn(:,:)
REAL, ALLOCATABLE :: rhoLRes(:,:,:,:,:)
CONTAINS
PROCEDURE,PASS :: init => moments_init
END TYPE t_moments
......@@ -391,14 +393,17 @@ SUBROUTINE mcd_init1(thisMCD,banddos,dimension,input,atoms,kpts)
END SUBROUTINE mcd_init1
SUBROUTINE moments_init(thisMoments,input,atoms)
SUBROUTINE moments_init(thisMoments,mpi,input,sphhar,atoms)
USE m_types_setup
USE m_types_mpi
IMPLICIT NONE
CLASS(t_moments), INTENT(INOUT) :: thisMoments
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_input), INTENT(IN) :: input
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_atoms), INTENT(IN) :: atoms
ALLOCATE(thisMoments%chmom(atoms%ntype,input%jspins))
......@@ -415,6 +420,11 @@ SUBROUTINE moments_init(thisMoments,input,atoms)
thisMoments%stdn = 0.0
thisMoments%svdn = 0.0
IF(mpi%irank.EQ.0) THEN
ALLOCATE(thisMoments%rhoLRes(atoms%jmtd,0:sphhar%nlhd,0:(atoms%lmaxd*(atoms%lmaxd+1))/2+atoms%lmaxd,atoms%ntype,4))
thisMoments%rhoLRes = 0.0
END IF
END SUBROUTINE moments_init
SUBROUTINE orbcomp_init(thisOrbcomp,input,banddos,dimension,atoms,kpts)
......
......@@ -274,7 +274,7 @@ CONTAINS
INTEGER,INTENT(IN) :: jspins, potden_type
CALL init_potden_simple(pd,stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,&
atoms%n_u,jspins,noco%l_noco,noco%l_mtnocopot,potden_type,&
atoms%n_u,jspins,noco%l_noco,noco%l_mtnocopot.OR.noco%l_mperp,potden_type,&
vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
END SUBROUTINE init_potden_types
......
......@@ -22,7 +22,7 @@ module m_mpmom
contains
subroutine mpmom( input, mpi, atoms, sphhar, stars, sym, cell, oneD, qpw, rho, potdenType, qlm )
subroutine mpmom( input, mpi, atoms, sphhar