Commit d864b81b authored by Gregor Michalicek's avatar Gregor Michalicek

Introduced wrapper for core density IO and applied it to some files

Note: There is a change in this commit that may or may not fix
or produce a bug that might be there or not: The rhc array in the
subroutine signature of cored was missing the spin dimension that
it has in cdngen. This is now corrected and the array is now correctly
filled in cored. It is unclear if this changes the behavior of the
code at some place. It might have an effect in cdnovlp as this array
also goes there.
parent cedf5f25
......@@ -5,7 +5,7 @@ CONTAINS
& rho,DIMENSION,&
& sphhar,&
& vr,&
& qint,rhc,seig)
& qint,rhc,tec,seig)
! *******************************************************
! ***** set up the core densities for compounds. *****
......@@ -29,16 +29,19 @@ CONTAINS
REAL, INTENT (OUT) :: seig
! ..
! .. Array Arguments ..
REAL , INTENT (IN) :: vr(atoms%jmtd,atoms%ntype)
REAL, INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd)
REAL, INTENT (OUT) :: rhc(DIMENSION%msh,atoms%ntype),qint(atoms%ntype,DIMENSION%jspd)
REAL, INTENT(IN) :: vr(atoms%jmtd,atoms%ntype)
REAL, INTENT(INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd)
REAL, INTENT(INOUT) :: rhc(DIMENSION%msh,atoms%ntype,DIMENSION%jspd)
REAL, INTENT(INOUT) :: qint(atoms%ntype,DIMENSION%jspd)
REAL, INTENT(INOUT) :: tec(atoms%ntype,DIMENSION%jspd)
! ..
! .. Local Scalars ..
REAL e,fj,fl,fn,q,rad,rhos,rhs,sea,sume,t2,tec
REAL e,fj,fl,fn,q,rad,rhos,rhs,sea,sume,t2
REAL d,dxx,rn,rnot,z,t1,rr,r,lambd,c,bmu,weight
INTEGER i,j,jatom,jm,korb,n,ncmsh,nm,nm1,nst ,l,ierr
! ..
! .. Local Arrays ..
REAL rhcs(DIMENSION%msh),rhoc(DIMENSION%msh),rhoss(DIMENSION%msh),vrd(DIMENSION%msh),f(0:3)
REAL occ(DIMENSION%nstd),a(DIMENSION%msh),b(DIMENSION%msh),ain(DIMENSION%msh),ahelp(DIMENSION%msh)
REAL occ_h(DIMENSION%nstd,2)
......@@ -48,26 +51,18 @@ CONTAINS
! ..
c = c_light(1.0)
seig = 0.
IF (jspin.EQ.1) THEN
OPEN (17,file='cdnc',form='unformatted',status='unknown')
ENDIF
!
IF (input%frcor) THEN
IF (jspin.EQ.1) REWIND 17
DO n = 1,atoms%ntype
jm = atoms%jri(n)
rnot = atoms%rmsh(1,n) ; dxx = atoms%dx(n)
ncmsh = NINT( LOG( (atoms%rmt(n)+10.0)/rnot ) / dxx + 1 )
ncmsh = MIN( ncmsh, DIMENSION%msh )
! ---> read in core density
READ (17) (rhc(i,n),i=1,ncmsh)
! ---> update spherical charge density
DO i = 1,atoms%jri(n)
rhoc(i) = rhc(i,n)
rhoc(i) = rhc(i,n,jspin)
rho(i,0,n,jspin) = rho(i,0,n,jspin) + rhoc(i)/sfp_const
ENDDO
! ---> read in kinetic enrgy of the core
READ (17) tec
! ---> for total energy calculations, determine the sum of the
! ---> eigenvalues by requiring that the core kinetic energy
! ---> remains constant.
......@@ -76,13 +71,11 @@ CONTAINS
ENDDO
nm = atoms%jri(n)
CALL intgr3(rhoc,atoms%rmsh(1,n),atoms%dx(n),nm,rhos)
sea = tec + rhos
WRITE (16,FMT=8030) n,jspin,tec,sea
WRITE (6,FMT=8030) n,jspin,tec,sea
sea = tec(n,jspin) + rhos
WRITE (16,FMT=8030) n,jspin,tec(n,jspin),sea
WRITE (6,FMT=8030) n,jspin,tec(n,jspin),sea
seig = seig + atoms%neq(n)*sea
ENDDO
! ---> read in qint
READ (17) (qint(n,jspin),n=1,atoms%ntype)
RETURN
END IF
......@@ -172,20 +165,17 @@ CONTAINS
rho(j,0,jatom,jspin) = rho(j,0,jatom,jspin) + rhoc(j)/sfp_const
ENDDO
rhc(1:ncmsh,jatom) = rhoss(1:ncmsh) / input%jspins
rhc(ncmsh+1:DIMENSION%msh,jatom) = 0.0
rhc(1:ncmsh,jatom,jspin) = rhoss(1:ncmsh) / input%jspins
rhc(ncmsh+1:DIMENSION%msh,jatom,jspin) = 0.0
seig = seig + atoms%neq(jatom)*sume
! WRITE (17) (rhoc(i),i=1,nm)
WRITE (17) (rhc(i,jatom),i=1,ncmsh)
DO i = 1,nm
rhoc(i) = rhoc(i)*vr(i,jatom)/atoms%rmsh(i,jatom)
ENDDO
CALL intgr3(rhoc,atoms%rmsh(1,jatom),atoms%dx(jatom),nm,rhs)
tec = sume - rhs
WRITE (6,FMT=8030) jatom,jspin,tec,sume
WRITE (16,FMT=8030) jatom,jspin,tec,sume
WRITE (17) tec
tec(jatom,jspin) = sume - rhs
WRITE (6,FMT=8030) jatom,jspin,tec(jatom,jspin),sume
WRITE (16,FMT=8030) jatom,jspin,tec(jatom,jspin),sume
! ---> simpson integration
rad = atoms%rmt(jatom)
......@@ -206,7 +196,7 @@ CONTAINS
WRITE(attributes(1),'(i0)') jatom
WRITE(attributes(2),'(i0)') NINT(z)
WRITE(attributes(3),'(i0)') jspin
WRITE(attributes(4),'(f18.10)') tec
WRITE(attributes(4),'(f18.10)') tec(jatom,jspin)
WRITE(attributes(5),'(f18.10)') sume
WRITE(attributes(6),'(f9.6)') q/input%jspins
CALL openXMLElementForm('coreStates',(/'atomType ','atomicNumber ','spin ','kinEnergy ',&
......@@ -229,14 +219,9 @@ CONTAINS
CALL closeXMLElement('coreStates')
ENDDO
! qint=0.
WRITE (17) (qint(n,jspin),n=1,atoms%ntype)
!
IF (jspin.EQ.input%jspins) CLOSE (17)
RETURN
8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f8.6,5x,&
8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f9.6,5x,&
& 'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x,&
& 'weight')
8010 FORMAT (12x,2f5.0,f6.1,f10.4,f10.0)
......
......@@ -9,6 +9,7 @@ CONTAINS
USE m_etabinit
USE m_spratm
USE m_ccdnup
USE m_cdn_io
USE m_types
IMPLICIT NONE
......@@ -165,18 +166,7 @@ CONTAINS
END DO ! loop over atoms (jatom)
!
!---->store core charge densities to file.17
OPEN (17,file='cdnc',form='unformatted',status='unknown')
REWIND 17
DO jspin = 1,input%jspins
DO jatom = 1,atoms%ntype
WRITE (17) (rhcs(j,jatom,jspin),j=1,atoms%jri(jatom))
WRITE (17) tecs(jatom,jspin)
END DO
WRITE (17) (qints(jatom,jspin),jatom=1,atoms%ntype)
END DO
CLOSE (17)
!
RETURN
!----> store core charge densities
CALL writeCoreDensity(input,atoms,dimension,rhcs,tecs,qints)
END SUBROUTINE coredr
END MODULE m_coredr
......@@ -22,6 +22,7 @@ MODULE m_cdn_io
PRIVATE
PUBLIC readDensity, writeDensity, isDensityFilePresent
PUBLIC readCoreDensity, writeCoreDensity
PUBLIC CDN_INPUT_DEN_const, CDN_OUTPUT_DEN_const
PUBLIC CDN_ARCHIVE_TYPE_CDN1_const, CDN_ARCHIVE_TYPE_NOCO_const
PUBLIC CDN_ARCHIVE_TYPE_CDN_const
......@@ -345,6 +346,104 @@ MODULE m_cdn_io
END SUBROUTINE writeDensity
SUBROUTINE readCoreDensity(input,atoms,dimension,rhcs,tecs,qints)
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_input),INTENT(IN) :: input
TYPE(t_dimension),INTENT(IN) :: DIMENSION
REAL, INTENT(OUT) :: rhcs(atoms%jmtd,atoms%ntype,DIMENSION%jspd)
REAL, INTENT(OUT) :: tecs(atoms%ntype,DIMENSION%jspd)
REAL, INTENT(OUT) :: qints(atoms%ntype,DIMENSION%jspd)
INTEGER :: mode, iUnit, iSpin, iAtom, i
LOGICAL :: l_exist
CHARACTER(LEN=30) :: filename
CALL getMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
IF (l_exist) THEN
!load density from cdn.hdf and exit subroutine
RETURN
ELSE
WRITE(*,*) 'cdn.hdf file not found.'
WRITE(*,*) 'Falling back to stream access file cdn.str.'
mode = CDN_STREAM_MODE
END IF
END IF
IF(mode.EQ.CDN_STREAM_MODE) THEN
INQUIRE(FILE='cdn.str',EXIST=l_exist)
IF (l_exist) THEN
!load density from cdn.str and exit subroutine
RETURN
ELSE
WRITE(*,*) 'cdn.str file not found.'
WRITE(*,*) 'Falling back to direct access file cdnc.'
mode = CDN_DIRECT_MODE
END IF
END IF
IF (mode.EQ.CDN_DIRECT_MODE) THEN
filename = 'cdnc'
INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
IF(.NOT.l_exist) THEN
CALL juDFT_error("core charge density file "//TRIM(ADJUSTL(filename))//" missing",calledby ="readCoreDensity")
END IF
iUnit = 17
OPEN (iUnit,file=TRIM(ADJUSTL(filename)),form='unformatted',status='unknown')
DO iSpin = 1,input%jspins
DO iAtom = 1,atoms%ntype
READ (iUnit) (rhcs(i,iAtom,iSpin),i=1,atoms%jri(iAtom))
READ (iUnit) tecs(iAtom,iSpin)
END DO
READ (iUnit) (qints(iAtom,iSpin),iAtom=1,atoms%ntype)
END DO
CLOSE (iUnit)
END IF
END SUBROUTINE readCoreDensity
SUBROUTINE writeCoreDensity(input,atoms,dimension,rhcs,tecs,qints)
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_input),INTENT(IN) :: input
TYPE(t_dimension),INTENT(IN) :: DIMENSION
REAL, INTENT(IN) :: rhcs(atoms%jmtd,atoms%ntype,DIMENSION%jspd)
REAL, INTENT(IN) :: tecs(atoms%ntype,DIMENSION%jspd)
REAL, INTENT(IN) :: qints(atoms%ntype,DIMENSION%jspd)
INTEGER :: mode, iUnit, iSpin, iAtom, i
CALL getMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
! Write core density to cdn.hdf file
STOP 'CDN_HDF5_MODE not yet implemented!'
ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
! Write core density to cdn.str file
STOP 'CDN_STREAM_MODE not yet implemented!'
ELSE
iUnit = 17
OPEN (iUnit,file='cdnc',form='unformatted',status='unknown')
DO iSpin = 1,input%jspins
DO iAtom = 1,atoms%ntype
WRITE (iUnit) (rhcs(i,iAtom,iSpin),i=1,atoms%jri(iAtom))
WRITE (iUnit) tecs(iAtom,iSpin)
END DO
WRITE (iUnit) (qints(iAtom,iSpin),iAtom=1,atoms%ntype)
END DO
CLOSE (iUnit)
END IF
END SUBROUTINE writeCoreDensity
SUBROUTINE getMode(mode)
INTEGER, INTENT(OUT) :: mode
......
......@@ -67,6 +67,7 @@
! .. Local Arrays ..
REAL stdn(atoms%ntype,dimension%jspd),svdn(atoms%ntype,dimension%jspd),alpha_l(atoms%ntype),&
rh(dimension%msh,atoms%ntype,dimension%jspd),qint(atoms%ntype,dimension%jspd)
REAL tec(atoms%ntype,DIMENSION%jspd),rhTemp(dimension%msh,atoms%ntype,dimension%jspd)
REAL chmom(atoms%ntype,dimension%jspd),clmom(3,atoms%ntype,dimension%jspd)
INTEGER,ALLOCATABLE :: igq_fft(:)
REAL ,ALLOCATABLE :: vz(:,:,:),vr(:,:,:,:)
......@@ -196,58 +197,68 @@
OPEN (15,file='ecore',status='unknown', action='write',form='unformatted')
ENDIF
DO jspin = 1,input%jspins
rh = 0.0
tec = 0.0
qint = 0.0
IF (input%frcor) THEN
CALL readCoreDensity(input,atoms,dimension,rh,tec,qint)
END IF
DO jspin = 1,input%jspins
IF ((input%jspins.EQ.2).AND.(mpi%irank.EQ.0)) THEN
DO n = 1,atoms%ntype
svdn(n,jspin) = rho(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
enddo
END IF
END DO
END IF
!
! block 1 unnecessary for slicing: begin
IF (.NOT.sliceplot%slice) THEN
! ---> add in core density
IF (mpi%irank.EQ.0) THEN
CALL cored(input,jspin,atoms, rho,dimension, sphhar, vr(:,0,:,jspin), qint,rh,seig)
results%seigc = results%seigc + seig
IF (input%jspins.EQ.2) THEN
DO n = 1,atoms%ntype
stdn(n,jspin) = rho(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
enddo
ENDIF
ENDIF ! mpi%irank = 0
IF (mpi%irank.EQ.0) THEN
CALL cored(input,jspin,atoms, rho,dimension, sphhar, vr(:,0,:,jspin), qint,rh,tec,seig)
rhTemp(:,:,jspin) = rh(:,:,jspin)
results%seigc = results%seigc + seig
IF (input%jspins.EQ.2) THEN
DO n = 1,atoms%ntype
stdn(n,jspin) = rho(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
END DO
END IF
END IF ! mpi%irank = 0
! ---> add core tail charge to qpw
IF ((noco%l_noco).AND.(mpi%irank.EQ.0)) THEN
!---> pk non-collinear
!---> add the coretail-charge to the constant interstitial
!---> charge (star 0), taking into account the direction of
!---> magnetisation of this atom
IF (jspin .EQ. 2) THEN
DO ityp = 1,atoms%ntype
rhoint = (qint(ityp,1) + qint(ityp,2)) /cell%volint/input%jspins/2.0
momint = (qint(ityp,1) - qint(ityp,2)) /cell%volint/input%jspins/2.0
!---> rho_11
qpw(1,1) = qpw(1,1) + rhoint + momint*cos(noco%beta(ityp))
!---> rho_22
qpw(1,2) = qpw(1,2) + rhoint - momint*cos(noco%beta(ityp))
!---> real part rho_21
cdom(1) = cdom(1) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
!---> imaginary part rho_21
cdom(1) = cdom(1) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
ENDDO
ENDIF
!---> pk non-collinear
ELSEIF (input%ctail) THEN
CALL cdnovlp(mpi,&
IF ((noco%l_noco).AND.(mpi%irank.EQ.0)) THEN
!---> pk non-collinear
!---> add the coretail-charge to the constant interstitial
!---> charge (star 0), taking into account the direction of
!---> magnetisation of this atom
IF (jspin .EQ. 2) THEN
DO ityp = 1,atoms%ntype
rhoint = (qint(ityp,1) + qint(ityp,2)) /cell%volint/input%jspins/2.0
momint = (qint(ityp,1) - qint(ityp,2)) /cell%volint/input%jspins/2.0
!---> rho_11
qpw(1,1) = qpw(1,1) + rhoint + momint*cos(noco%beta(ityp))
!---> rho_22
qpw(1,2) = qpw(1,2) + rhoint - momint*cos(noco%beta(ityp))
!---> real part rho_21
cdom(1) = cdom(1) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
!---> imaginary part rho_21
cdom(1) = cdom(1) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
END DO
END IF
!---> pk non-collinear
ELSE IF (input%ctail) THEN
CALL cdnovlp(mpi,&
sphhar,stars,atoms,sym, dimension,vacuum, cell, input,oneD,l_st, jspin,rh, qpw,rhtxy,rho,rht)
ELSEIF (mpi%irank.EQ.0) THEN
ELSE IF (mpi%irank.EQ.0) THEN
DO ityp = 1,atoms%ntype
qpw(1,jspin) = qpw(1,jspin) + qint(ityp,jspin)/input%jspins/cell%volint
ENDDO
END IF
END DO
END IF
! block 1 unnecessary for slicing: end
END IF
!
ENDDO
END DO ! loop over spins
CALL writeCoreDensity(input,atoms,dimension,rhTemp,tec,qint)
IF ((input%gw.eq.1 .or. input%gw.eq.3).AND.(mpi%irank.EQ.0)) CLOSE(15)
ELSE
......
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