Commit 01525471 authored by Daniel Wortmann's avatar Daniel Wortmann

Fixes to make sure occupations are used correctly in starting density generation

parent c2013add
......@@ -149,7 +149,7 @@ SUBROUTINE initParallelProcesses(atoms,vacuum,input,stars,sliceplot,banddos,&
ALLOCATE(atoms%lo1l(0:atoms%llod,atoms%ntype))
ALLOCATE(atoms%nlol(0:atoms%llod,atoms%ntype))
ALLOCATE(atoms%coreStateOccs(dimension%nstd,2,atoms%ntype))
ALLOCATE(atoms%coreStateOccs(dimension%nstd,2,atoms%ntype));atoms%coreStateOccs=0.0
ALLOCATE(atoms%coreStateNprnc(dimension%nstd,atoms%ntype))
ALLOCATE(atoms%coreStateKappa(dimension%nstd,atoms%ntype))
......
......@@ -87,6 +87,7 @@ CONTAINS
ENDIF
!-odim
ALLOCATE ( atoms%nz(atoms%ntype),atoms%relax(3,atoms%ntype),atoms%nlhtyp(atoms%ntype))
ALLOCATE (atoms%corestateoccs(1,2,atoms%ntype));atoms%corestateoccs=0.0
ALLOCATE ( sphhar%clnu(sphhar%memd,0:sphhar%nlhd,sphhar%ntypsd),stars%ustep(stars%ng3) )
ALLOCATE ( stars%ig(-stars%mx1:stars%mx1,-stars%mx2:stars%mx2,-stars%mx3:stars%mx3),stars%ig2(stars%ng3) )
ALLOCATE ( atoms%jri(atoms%ntype),stars%kv2(2,stars%ng2),stars%kv3(3,stars%ng3),sphhar%llh(0:sphhar%nlhd,sphhar%ntypsd) )
......
......@@ -1269,7 +1269,7 @@ SUBROUTINE r_inpXML(&
dimension%nstd = 29
ALLOCATE(atoms%coreStateOccs(dimension%nstd,2,atoms%ntype))
ALLOCATE(atoms%coreStateOccs(dimension%nstd,2,atoms%ntype));atoms%coreStateOccs=0.0
ALLOCATE(atoms%coreStateNprnc(dimension%nstd,atoms%ntype))
ALLOCATE(atoms%coreStateKappa(dimension%nstd,atoms%ntype))
......
......@@ -37,14 +37,6 @@
IF (lmax.GT.lmaxd) THEN
CALL juDFT_error("lmaxd too small in ylm4")
!$OMP MASTER
! WRITE(6,*) ' calling ylmnorm, lmax=',lmax,lmaxd
!--> first deallocate the array if it exists
IF ( allocated(ynorm) ) DEALLOCATE(ynorm)
ALLOCATE ( ynorm( (lmax+1)**2 ) ) ! allocate array
lmaxd = lmax
!CALL ylmnorm
!$OMP END MASTER
ENDIF
!---> calculate sin and cos of theta and phi
......
......@@ -61,7 +61,7 @@ SUBROUTINE stden(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,&
INTEGER lnum(DIMENSION%nstd,atoms%ntype),nst(atoms%ntype)
INTEGER jrc(atoms%ntype)
LOGICAL l_found(0:3),llo_found(atoms%nlod),l_enpara,l_st
REAL :: occ(SIZE(atoms%corestateoccs,1),2)
! Data statements
DATA del/1.e-6/
PARAMETER (l_st=.true.)
......@@ -113,12 +113,14 @@ SUBROUTINE stden(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,&
ELSE
bm = 0.
END IF
occ=atoms%coreStateOccs(:,:,n)
! check whether this atom has been done already
DO n1 = 1, n - 1
IF (ABS(z-atoms%zatom(n1)).GT.del) CYCLE
IF (ABS(r-atoms%rmt(n1)).GT.del) CYCLE
IF (ABS(h-atoms%dx(n1)).GT.del) CYCLE
IF (ABS(bm-atoms%bmu(n1)).GT.del) CYCLE
IF (ANY(ABS(occ(:,:)-atoms%coreStateOccs(:,:,n1))>del)) CYCLE
IF (jr.NE.atoms%jri(n1)) CYCLE
DO ispin = 1, input%jspins
DO i = 1,jrc(n) ! dimension%msh
......
......@@ -16,7 +16,6 @@ CONTAINS
& qpw,rho,&
& qlm)
#include"cpp_double.h"
USE m_types
IMPLICIT NONE
TYPE(t_mpi),INTENT(IN) :: mpi
......@@ -131,7 +130,7 @@ CONTAINS
SUBROUTINE pw_moments(mpi,stars,atoms,cell,&
& sym,oneD,&
& qpw,qlmp)
& qpw_in,qlmp_out)
!multipole moments of plane wave charge inside the spheres (q_{lm}^{Ii})
USE m_phasy1
USE m_sphbes
......@@ -145,20 +144,23 @@ CONTAINS
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_atoms),INTENT(IN) :: atoms
COMPLEX,INTENT(IN) :: qpw(:)
COMPLEX,INTENT(OUT):: qlmp(-atoms%lmaxd:,0:,:)
COMPLEX,INTENT(IN) :: qpw_in(:)
COMPLEX,INTENT(OUT):: qlmp_out(-atoms%lmaxd:,0:,:)
!locals
INTEGER:: n,k,l,ll1 ,lm,ierr(3),m
COMPLEX sk3i,cil,nqpw
COMPLEX pylm( (MAXVAL(atoms%lmax)+1)**2 ,atoms%ntype )
REAL::sk3r,rl3
REAL aj(0:MAXVAL(atoms%lmax)+1)
COMPLEX, ALLOCATABLE :: c_b(:)
COMPLEX, ALLOCATABLE :: qpw(:)
LOGICAL :: od
#ifdef CPP_MPI
INCLUDE 'mpif.h'
EXTERNAL MPI_REDUCE
#endif
COMPLEX qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
ALLOCATE(qpw(stars%ng3))
qpw=qpw_in(:stars%ng3)
qlmp= 0.0
IF (mpi%irank==0) THEN
!g eq 0 term : \sqrt{4 \pi}/3 R_i^3 \rho_I(0) \delta_{l,0}
......@@ -168,8 +170,7 @@ CONTAINS
ENDDO
ENDIF
#ifdef CPP_MPI
CALL MPI_BCAST(qpw,SIZE(qpw),CPP_MPI_COMPLEX,0,&
& mpi%mpi_comm,ierr)
CALL MPI_BCAST(qpw,SIZE(qpw),MPI_DOUBLE_COMPLEX,0, mpi%mpi_comm,ierr)
#endif
! g ne 0 terms : \sum_{K \= 0} 4 \pi i^l \rho_I(K) R_i^{l+3} \times
! j_{l+1} (KR_i) / KR_i \exp{iK\xi_i} Y^*_{lm} (K)
......@@ -199,7 +200,7 @@ CONTAINS
cil = aj(l+1)*sk3i*rl3
ll1 = l*(l+1) + 1
DO m = -l,l
lm = ll1 + m
lm = ll1 + m
qlmp(m,l,n) = qlmp(m,l,n) + cil*pylm(lm,n)
ENDDO
rl3 = rl3*atoms%rmt(n)
......@@ -207,15 +208,11 @@ CONTAINS
ENDDO ! n = 1, atoms%ntype
ENDDO ! k = 2, stars%ng3
!$OMP END PARALLEL DO
#ifdef CPP_MPI
n = SIZE(qlmp)
ALLOCATE(c_b(n))
CALL MPI_REDUCE(qlmp,c_b,n,CPP_MPI_COMPLEX,MPI_SUM,0,&
& mpi%mpi_comm,ierr)
IF (mpi%irank.EQ.0) THEN
qlmp=RESHAPE(c_b,(/SIZE(qlmp,1),SIZE(qlmp,2),SIZE(qlmp,3)/))
ENDIF
DEALLOCATE (c_b)
#ifdef CPP_MPII
PRINT *,"mpi",mpi%irank,qlmp(0,0,:)
CALL MPI_REDUCE(qlmp,qlmp_out,SIZE(qlmp),MPI_DOUBLE_COMPLEX,MPI_SUM,0,mpi%mpi_comm,ierr)
#else
qlmp_out=qlmp
#endif
END SUBROUTINE pw_moments
......
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