Commit 00b6e7d7 authored by Matthias Redies's avatar Matthias Redies

Merge branch 'develop' into hybrid_purge

parents 647681ef 36edbcaa
......@@ -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)
......
......@@ -36,7 +36,7 @@ CONTAINS
! .. Local Scalars ..
INTEGER iri
REAL b_xavh,scale,b_con_outx,b_con_outy,mx,my,mz,&
& alphh,betah,mz_tmp,mx_mix,my_mix,mz_mix
& alphh,betah,mz_tmp,mx_mix,my_mix,mz_mix,absmag
REAL rho11,rho22, alphdiff
COMPLEX rho21
! ..
......@@ -56,11 +56,12 @@ CONTAINS
mx = 2*REAL(qa21(itype))
my = 2*AIMAG(qa21(itype))
mz = chmom(itype,1) - chmom(itype,2)
WRITE (6,8025) mx,my
absmag=SQRT(mx*mx+my*my+mz*mz)
WRITE (6,8025) mx,my,mz,absmag
!---> determine the polar angles of the moment vector in the local frame
CALL pol_angle(mx,my,mz,betah,alphh)
WRITE (6,8026) betah,alphh
8025 FORMAT(2x,'--> local frame: ','mx=',f9.5,' my=',f9.5)
8025 FORMAT(2x,'--> local frame: ','mx=',f9.5,' my=',f9.5,' mz=',f9.5,' |m|=',f9.5)
8026 FORMAT(2x,'-->',10x,' delta beta=',f9.5,&
& ' delta alpha=',f9.5)
......
......@@ -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
......@@ -24,5 +27,7 @@ cdn_mt/rhosphnlo.f90
cdn_mt/setabc1locdn.f90
cdn_mt/setabc1locdn1.f90
cdn_mt/calcDenCoeffs.f90
cdn_mt/magnMomFromDen.f90
cdn_mt/alignSpinAxisMagn.f90
)
!--------------------------------------------------------------------------------
! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and avhttps://gcc.gnu.org/onlinedocs/gfortran/SQRT.htmlailable as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!------------------------------------------------------------------------------
! This routine allows to rotate the cdn in a way that the direction of magnetization aligns with the direction of the spin quantization axis.
! This routine also allows to reverse the rotation by using the angles stored in atoms (phi_mt_avg,theta_mt_avg) which are generated by the
! routine magnMomFromDen.
!
! Robin Hilgers, Nov '19
MODULE m_alignSpinAxisMagn
USE m_magnMomFromDen
USE m_types
USE m_flipcdn
USE m_constants
USE m_polangle
CONTAINS
SUBROUTINE rotateMagnetToSpinAxis(vacuum,sphhar,stars&
,sym,oneD,cell,noco,input,atoms,den)
TYPE(t_input), INTENT(INOUT) :: input
TYPE(t_atoms), INTENT(INOUT) :: atoms
TYPE(t_noco), INTENT(INOUT) :: noco
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_potden), INTENT(INOUT) :: den
REAL :: moments(atoms%ntype,3)
REAL :: phiTemp(atoms%ntype),thetaTemp(atoms%ntype),pi
INTEGER :: i
pi=pimach()
!!TEMP
! REAL :: x,y,z
phiTemp=noco%alph
thetaTemp=noco%beta
CALL magnMomFromDen(input,atoms,noco,den,moments)
! DO i=1, atoms%ntype
! IF (abs(atoms%theta_mt_avg(i)).LE. 0.001) THEN
! atoms%phi_mt_avg(i)=0.0
! atoms%theta_mt_avg(i)=0.0
! END IF
! END DO
!write(*,*) "mx1"
!write(*,*) moments(1,1)
!write(*,*) "mz1"
!write(*,*) moments(1,3)
!write(*,*) "mx2"
!write(*,*) moments(2,1)
!write(*,*) "mz2"
!write(*,*) moments(2,3)
CALL flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell,-atoms%phi_mt_avg,-atoms%theta_mt_avg,den)
!write (*,*)"mx my mz"
!CALL sphericaltocart(SQRT(moments(1,1)**2+moments(1,2)**2+moments(1,3)**2),thetaTemp(1),phiTemp(1),x,y,z)
!write(*,*) x,y,z
!CALL sphericaltocart(SQRT(moments(2,1)**2+moments(2,2)**2+moments(2,3)**2),thetaTemp(2),phiTemp(2),x,y,z)
!write(*,*) x,y,z
!write(*,*) "atoms%phi_mt_avg"
!!write(*,*) atoms%phi_mt_avg
!write(*,*) "atoms%theta_mt_avg"
!write(*,*) atoms%theta_mt_avg
noco%alph=mod(atoms%phi_mt_avg+phiTemp,2*pimach())
noco%beta=mod(atoms%theta_mt_avg+thetaTemp,2*pimach())
atoms%phi_mt_avg=noco%alph
atoms%theta_mt_avg=noco%beta
DO i=1, atoms%ntype
IF(noco%alph(i)<0) noco%alph(i)=noco%alph(i)+2*pi
IF(noco%beta(i)<0) THEN
noco%beta(i)=-noco%beta(i)
noco%alph=noco%alph+pi
END IF
IF(noco%beta(i)>pi) THEN
noco%beta(i)=pi-mod(noco%beta(i),pi)
noco%alph(i)=noco%alph(i)+pi
END IF
noco%alph=mod(noco%alph,2*pi)
End Do
write(*,*) "Noco Phi"
write(*,*) noco%alph
write(*,*) "Noco Theta"
write(*,*) noco%beta
END SUBROUTINE rotateMagnetToSpinAxis
SUBROUTINE rotateMagnetFromSpinAxis(noco,vacuum,sphhar,stars&
,sym,oneD,cell,input,atoms,den,inDen)
TYPE(t_input), INTENT(INOUT) :: input
TYPE(t_atoms), INTENT(INOUT) :: atoms
TYPE(t_noco), INTENT(INOUT) :: noco
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_potden), INTENT(INOUT) :: den, inDen
CALL flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell,atoms%phi_mt_avg,atoms%theta_mt_avg,den)
CALL flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell,atoms%phi_mt_avg,atoms%theta_mt_avg,inDen)
atoms%flipSpinPhi=0
atoms%flipSpinTheta=0
noco%alph=0
noco%beta=0
END SUBROUTINE rotateMagnetFromSpinAxis
END MODULE m_alignSpinAxisMagn
......@@ -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) )&
......@@ -191,7 +195,9 @@ CONTAINS
!rho21(j,0,itype) = rho21(j,0,itype)+ conjg(cs)/(atoms%neq(itype)*sfp_const)
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)
rho(j,0,itype,4)=rho(j,0,itype,4)+aimag(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)&
......@@ -210,7 +217,9 @@ CONTAINS
!rho21(j,lh,itype)= rho21(j,lh,itype)+ CONJG(cs)/atoms%neq(itype)
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)
rho(j,lh,itype,4)=rho(j,lh,itype,4)+aimag(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
! 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_magDiMom
CONTAINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This subroutine calculates intraatomic magnetic dipole moments.
!
! GM'2018
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE magDiMom(input,atoms,sphhar,noco,l_fmpl2,rho,magDipoles,elecDipoles)
USE m_constants
USE m_types
USE m_juDFT
USE m_rotdenmat
USE m_lattHarmsSphHarmsConv
USE m_gaunt
USE m_intgr
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
REAL, INTENT(IN) :: rho(:,0:,:,:) ! if l_fmpl last dimension is 4, otherwise 2.
LOGICAL, INTENT(IN) :: l_fmpl2
REAL, INTENT(INOUT) :: magDipoles(:,:)
REAL, INTENT(INOUT) :: elecDipoles(:,:)
REAL, ALLOCATABLE :: inRho(:,:,:,:)
COMPLEX, ALLOCATABLE :: rhoSphHarms(:,:,:,:), rhoTempSphHarms(:,:,:,:)
COMPLEX, ALLOCATABLE :: rhoSphHarmsR(:,:,:)
INTEGER :: iType, ilh, l, m, lm, lp, mp, lmp, i
REAL :: theta, phi, cdn11, cdn22
REAL :: magDipole(3), myCharge, elecDipole(3)
REAL :: constA
COMPLEX :: gauntA, gauntB, gauntC
COMPLEX :: cdn21
IF(input%jspins.EQ.1) RETURN
IF(.NOT.noco%l_noco) RETURN
!---> calculate the charge and magnetization density in the muffin tins
ALLOCATE(inRho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,4))
DO iType = 1,atoms%ntype
IF (.NOT.l_fmpl2) THEN
theta = noco%beta(iType)
phi = noco%alph(iType)
inRho(:,:,iType,1) = rho(:,:,iType,1) + rho(:,:,iType,2)
inRho(:,:,iType,2) = rho(:,:,iType,1) - rho(:,:,iType,2)
inRho(:,:,iType,3) = inRho(:,:,iType,2) * SIN(phi)*SIN(theta)
inRho(:,:,iType,4) = inRho(:,:,iType,2) * COS(theta)
inRho(:,:,iType,2) = inRho(:,:,iType,2) * COS(phi)*SIN(theta)
ELSE
DO ilh = 0,sphhar%nlh(atoms%ntypsy(iType))
DO i = 1,atoms%jri(iType)
cdn11 = rho(i,ilh,iType,1)
cdn22 = rho(i,ilh,iType,2)
cdn21 = CMPLX(rho(i,ilh,iType,3),rho(i,ilh,iType,4))
CALL rot_den_mat(noco%alph(iType),noco%beta(iType),cdn11,cdn22,cdn21)
inRho(i,ilh,iType,1) = cdn11 + cdn22
inRho(i,ilh,iType,2) = 2.0 * REAL(cdn21)
! Note: The minus sign in the following line is temporary to adjust for differences in the offdiagonal
! part of the density between this fleur version and ancient (v0.26) fleur.
inRho(i,ilh,iType,3) = -2.0 * AIMAG(cdn21)
inRho(i,ilh,iType,4) = cdn11 - cdn22
END DO
END DO
END IF
END DO
ALLOCATE (rhoSphHarms(atoms%jmtd,(atoms%lmaxd+1)**2,atoms%ntype,4))
ALLOCATE (rhoTempSphHarms(atoms%jmtd,(atoms%lmaxd+1)**2,atoms%ntype,4))
ALLOCATE (rhoSphHarmsR(atoms%jmtd,(atoms%lmaxd+1)**2,atoms%ntype))
rhoSphHarms = CMPLX(0.0,0.0)
DO i = 1, 4
DO iType = 1, atoms%ntype
CALL lattHarmsRepToSphHarms(atoms,sphhar,iType,inRho(:,0:,iType,i),rhoSphHarms(:,:,iType,i))
END DO
END DO
! electric dipole moment (start)
DO iType = 1, atoms%ntype
DO i = 1, atoms%jri(iType)
rhoSphHarmsR(i,:,iType) = rhoSphHarms(i,:,iType,1) * atoms%rmsh(i,iType)
END DO
END DO
constA = SQRT(2.0*pi_const/3.0)
rhoTempSphHarms = CMPLX(0.0,0.0)
DO iType = 1, atoms%ntype
DO lp = 0, MIN(2,atoms%lmax(iType))
DO mp = -lp, lp
DO l = 0, MIN(2,atoms%lmax(iType))
DO m = -l, l
!note 1: For refinement maybe I could make use of selection rules.
!note 2: ls for r^\hat is 1, ms is -1..1.
gauntA = gaunt1(lp,1,l,mp,-1,m,atoms%lmaxd)
gauntB = gaunt1(lp,1,l,mp,0,m,atoms%lmaxd)
gauntC = gaunt1(lp,1,l,mp,1,m,atoms%lmaxd)
lmp = lp*(lp+1) + mp + 1
lm = l*(l+1) + m + 1
DO i = 1, atoms%jri(iType)
rhoTempSphHarms(i,lmp,iType,2) = rhoTempSphHarms(i,lmp,iType,2) +&
constA*(gauntA-gauntC)*rhoSphHarmsR(i,lm,iType)
rhoTempSphHarms(i,lmp,iType,3) = rhoTempSphHarms(i,lmp,iType,3) +&
constA*CMPLX(0.0,1.0)*(gauntA+gauntC)*rhoSphHarmsR(i,lm,iType)
rhoTempSphHarms(i,lmp,iType,4) = rhoTempSphHarms(i,lmp,iType,4) +&
constA*SQRT(2.0)*gauntB*rhoSphHarmsR(i,lm,iType)
END DO
END DO
END DO
END DO
END DO
END DO
elecDipole = 0.0
DO iType = 1, atoms%ntype
CALL intgr3(REAL(rhoTempSphHarms(:,1,iType,2)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),elecDipole(1))
CALL intgr3(REAL(rhoTempSphHarms(:,1,iType,3)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),elecDipole(2))
CALL intgr3(REAL(rhoTempSphHarms(:,1,iType,4)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),elecDipole(3))
elecDipoles(:,iType) = elecDipole(:) * sfp_const
END DO
! electric dipole moment (end)
! WRITE(7534,*) '===================================================='
DO iType = 1, atoms%ntype
DO l = 0, 2
DO m = -l, l
magDipole(:) = 0.0
myCharge = 0.0
lm = l*(l+1) + m + 1
CALL intgr3(REAL(rhoSphHarms(:,lm,iType,1)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),myCharge)
CALL intgr3(REAL(rhoSphHarms(:,lm,iType,2)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),magDipole(1))
CALL intgr3(REAL(rhoSphHarms(:,lm,iType,3)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),magDipole(2))
CALL intgr3(REAL(rhoSphHarms(:,lm,iType,4)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),magDipole(3))
! WRITE(7534,'(3i7,4e24.8)') iType, l, m, myCharge, magDipole(:)
END DO
END DO
rhoSphHarms(:,:,iType,1) = CMPLX(0.0,0.0)
END DO
constA = SQRT(2.0*pi_const/3.0)
DO iType = 1, atoms%ntype
DO lp = 0, MIN(2,atoms%lmax(iType))
DO mp = -lp, lp
DO l = 0, MIN(2,atoms%lmax(iType))
DO m = -l, l
!note 1: For refinement maybe I could make use of selection rules.
!note 2: ls for r^\hat is 1, ms is -1..1.
gauntA = gaunt1(lp,1,l,mp,-1,m,atoms%lmaxd)
gauntB = gaunt1(lp,1,l,mp,0,m,atoms%lmaxd)
gauntC = gaunt1(lp,1,l,mp,1,m,atoms%lmaxd)
lmp = lp*(lp+1) + mp + 1
lm = l*(l+1) + m + 1
DO i = 1, atoms%jri(iType)
rhoSphHarms(i,lmp,iType,1) = rhoSphHarms(i,lmp,iType,1) +&
constA*(gauntA-gauntC)*rhoSphHarms(i,lm,iType,2)
rhoSphHarms(i,lmp,iType,1) = rhoSphHarms(i,lmp,iType,1) +&