Commit c2c0b513 authored by Gregor Michalicek's avatar Gregor Michalicek

Introduced core density IO wrappers to all relevant places

parent d864b81b
MODULE m_force_a4
CONTAINS
SUBROUTINE force_a4(atoms,sphhar,input,&
SUBROUTINE force_a4(atoms,sphhar,input,dimension,&
& vr,&
& force)
!
......@@ -13,10 +13,12 @@ CONTAINS
USE m_constants, ONLY : sfp_const
USE m_differentiate,ONLY: difcub
USE m_types
USE m_cdn_io
IMPLICIT NONE
TYPE(t_input),INTENT(IN) :: input
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_input),INTENT(IN) :: input
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_dimension),INTENT(IN) :: dimension
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
......@@ -28,7 +30,8 @@ CONTAINS
! ..
! .. Local Arrays ..
COMPLEX forc_a4(3),gv(3),ycomp1(3,-1:1)
REAL rhoaux(atoms%jmtd),rhoc(atoms%jmtd)
REAL rhoaux(atoms%jmtd),rhoc(atoms%jmtd,atoms%ntype,dimension%jspd)
REAL tec(atoms%ntype,dimension%jspd),qintc(atoms%ntype,dimension%jspd)
! ..
! .. Data statements ..
COMPLEX,PARAMETER:: czero=(0.000,0.000)
......@@ -47,28 +50,18 @@ CONTAINS
ycomp1(1,1) = CMPLX(-s23,0.0)
ycomp1(2,1) = CMPLX(0.0,-s23)
ycomp1(3,1) = czero
!
OPEN (17,file='cdnc',form='unformatted',status='old')
DO jsp = 1,input%jspins
IF (jsp.EQ.1) THEN
REWIND 17
END IF
! ---> read in core density
CALL readCoreDensity(input,atoms,dimension,rhoc,tec,qintc)
DO jsp = 1,input%jspins
na = 1
DO n = 1,atoms%ntype
! ---> read in core density
READ (17) (rhoc(ir),ir=1,atoms%jri(n))
! ---> core density is multiplied by r**2*sqrt(4*pi)
! ---> such that intrg0 gives correct core charge
!
! ---> read in kinetic enrgy of the core
READ (17)
!
DO n = 1,atoms%ntype
IF (atoms%l_geo(n)) THEN
DO i = 1,3
forc_a4(i) = czero
END DO
!
CALL intgr0(rhoc,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),qcore)
CALL intgr0(rhoc(:,n,jsp),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),qcore)
! write(16,1616) qcore
8000 FORMAT (' FORCE_A4: core charge=',1p,e16.8)
!
......@@ -97,22 +90,22 @@ CONTAINS
! difcub performs analytic derivative of Lagrangian of 3rd order
!
xi = atoms%rmsh(1,n)
rhoaux(1) = difcub(atoms%rmsh(1,n),vr(1,lh,n,jsp),xi)*rhoc(1)
rhoaux(1) = difcub(atoms%rmsh(1,n),vr(1,lh,n,jsp),xi)*rhoc(1,n,jsp)
DO ir = 2,atoms%jri(n) - 2
xi = atoms%rmsh(ir,n)
rhoaux(ir) = difcub(atoms%rmsh(ir-1,n),&
& vr(ir-1,lh,n,jsp),xi) * rhoc(ir)
& vr(ir-1,lh,n,jsp),xi) * rhoc(ir,n,jsp)
END DO
!
ir = atoms%jri(n) - 1
xi = atoms%rmsh(ir,n)
rhoaux(ir) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
& vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir)
& vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir,n,jsp)
!
ir = atoms%jri(n)
xi = atoms%rmsh(ir,n)
rhoaux(ir) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
& vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir)
& vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir,n,jsp)
CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
a4_1 = 0.5*w/sfp_const
!
......@@ -120,7 +113,7 @@ CONTAINS
! note: rhocore is already multiplied by r**2 and srt(4.*pi)
!
DO ir = 1,atoms%jri(n)
rhoaux(ir) = rhoc(ir)/atoms%rmsh(ir,n)*vr(ir,lh,n,jsp)
rhoaux(ir) = rhoc(ir,n,jsp)/atoms%rmsh(ir,n)*vr(ir,lh,n,jsp)
END DO
!
CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
......@@ -151,11 +144,7 @@ CONTAINS
ENDIF
na = na + atoms%neq(n)
ENDDO
READ (17)
! spin loop ends
ENDDO
!
CLOSE (17)
RETURN
END SUBROUTINE force_a4
END MODULE m_force_a4
......@@ -21,7 +21,8 @@ MODULE m_cdn_io
IMPLICIT NONE
PRIVATE
PUBLIC readDensity, writeDensity, isDensityFilePresent
PUBLIC readDensity, writeDensity
PUBLIC isDensityFilePresent, isCoreDensityPresent
PUBLIC readCoreDensity, writeCoreDensity
PUBLIC CDN_INPUT_DEN_const, CDN_OUTPUT_DEN_const
PUBLIC CDN_ARCHIVE_TYPE_CDN1_const, CDN_ARCHIVE_TYPE_NOCO_const
......@@ -494,4 +495,37 @@ MODULE m_cdn_io
isDensityFilePresent = l_exist
END FUNCTION isDensityFilePresent
LOGICAL FUNCTION isCoreDensityPresent()
LOGICAL :: l_exist
INTEGER :: mode
CALL getMode(mode)
IF (mode.EQ.CDN_HDF5_MODE) THEN
INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
IF(l_exist) THEN
STOP 'Not yet implemented!'
RETURN
END IF
END IF
IF ((mode.EQ.CDN_STREAM_MODE).OR.(mode.EQ.CDN_HDF5_MODE)) THEN
INQUIRE(FILE='cdn.str',EXIST=l_exist)
IF(l_exist) THEN
STOP 'Not yet implemented!'
RETURN
END IF
END IF
!cdnc should be enough for any mode...
INQUIRE(FILE='cdnc',EXIST=l_exist)
IF (l_exist) THEN
isCoreDensityPresent = l_exist
RETURN
END IF
isCoreDensityPresent = .FALSE.
END FUNCTION isCoreDensityPresent
END MODULE m_cdn_io
......@@ -744,7 +744,7 @@
CALL timestart('determination of total energy')
CALL totale(atoms,sphhar,stars,vacuum,&
CALL totale(atoms,sphhar,stars,vacuum,dimension,&
sym,input,noco,cell,oneD,xcpot,hybrid,it,results)
CALL timestop('determination of total energy')
......
......@@ -99,22 +99,22 @@ CONTAINS
IF (noco%l_noco) THEN
cdnfname = 'cdn'
CALL plotdop(&
& oneD,stars,vacuum,sphhar,atoms,&
& oneD,dimension,stars,vacuum,sphhar,atoms,&
& input,sym,cell,sliceplot,&
& noco%l_noco,cdnfname)
cdnfname = 'mdnx'
CALL plotdop(&
& oneD,stars,vacuum,sphhar,atoms,&
& oneD,dimension,stars,vacuum,sphhar,atoms,&
& input,sym,cell,sliceplot,&
& noco%l_noco,cdnfname)
cdnfname = 'mdny'
CALL plotdop(&
& oneD,stars,vacuum,sphhar,atoms,&
& oneD,dimension,stars,vacuum,sphhar,atoms,&
& input,sym,cell,sliceplot,&
& noco%l_noco,cdnfname)
cdnfname = 'mdnz'
CALL plotdop(&
& oneD,stars,vacuum,sphhar,atoms,&
& oneD,dimension,stars,vacuum,sphhar,atoms,&
& input,sym,cell,sliceplot,&
& noco%l_noco,cdnfname)
ELSE
......@@ -124,7 +124,7 @@ CONTAINS
cdnfname = 'cdn1'
ENDIF
CALL plotdop(&
& oneD,stars,vacuum,sphhar,atoms,&
& oneD,dimension,stars,vacuum,sphhar,atoms,&
& input,sym,cell,sliceplot,&
& noco%l_noco,cdnfname)
ENDIF
......
......@@ -5,7 +5,7 @@
!--------------------------------------------------------------------------------
MODULE m_totale
CONTAINS
SUBROUTINE totale(atoms,sphhar,stars,vacuum, &
SUBROUTINE totale(atoms,sphhar,stars,vacuum,dimension, &
sym,input,noco,cell,oneD, xcpot,hybrid, it,results)
!
! ***************************************************
......@@ -64,6 +64,7 @@ CONTAINS
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_dimension),INTENT(IN) :: dimension
! ..
! .. Scalar Arguments ..
INTEGER,INTENT (IN) :: it
......@@ -163,7 +164,7 @@ CONTAINS
CALL readPotential(stars,vacuum,atoms,sphhar,input,sym,POT_ARCHIVE_TYPE_TOT_const,&
iter,vr,vpw,vz,vxy)
CALL force_a4(atoms,sphhar,input, vr, results%force)
CALL force_a4(atoms,sphhar,input,dimension, vr, results%force)
!
ENDIF
!
......
......@@ -88,7 +88,8 @@ CONTAINS
COMPLEX, ALLOCATABLE :: alphm(:,:)
COMPLEX, ALLOCATABLE :: excpw(:),excxy(:,:,:),vpw_w(:,:),psq(:)
REAL, ALLOCATABLE :: vbar(:),af1(:),bf1(:),xp(:,:)
REAL, ALLOCATABLE :: rhoc(:),rhoc_vx(:)
REAL, ALLOCATABLE :: rhoc(:,:,:),rhoc_vx(:)
REAL, ALLOCATABLE :: tec(:,:), qintc(:,:)
!.....density
REAL, ALLOCATABLE :: rho(:,:,:,:),rht(:,:,:)
COMPLEX, ALLOCATABLE :: qpw(:,:),rhtxy(:,:,:,:)
......@@ -691,9 +692,10 @@ CONTAINS
!HF kinetic energy correction for core states
IF ( xcpot%icorr == icorr_pbe0 .OR. xcpot%icorr == icorr_hse .OR.&
xcpot%icorr == icorr_hf .OR. xcpot%icorr == icorr_vhse ) THEN
OPEN (17,file='cdnc',form='unformatted',status='unknown')
REWIND 17
ALLOCATE( rhoc(atoms%jmtd), rhoc_vx(atoms%jmtd) )
ALLOCATE(rhoc(atoms%jmtd,atoms%ntype,DIMENSION%jspd), rhoc_vx(atoms%jmtd))
ALLOCATE(tec(atoms%ntype,DIMENSION%jspd),qintc(atoms%ntype,DIMENSION%jspd))
CALL readCoreDensity(input,atoms,dimension,rhoc,tec,qintc)
DEALLOCATE(tec,qintc)
END IF
!HF
......@@ -712,12 +714,8 @@ CONTAINS
xcpot%icorr == icorr_hf .OR. xcpot%icorr == icorr_vhse ) THEN
nat = 1
DO n = 1,atoms%ntype
READ (17) ( rhoc(j), j = 1,atoms%jri(n) )
! Skip over parts in cdnc not used
READ (17) rdum
! calculate exchange correction to kinetic energy by core states
DO j = 1, atoms%jri(n)
rhoc_vx(j) = rhoc(j) * vxr(j,0,n,js) / sfp_const
rhoc_vx(j) = rhoc(j,n,js) * vxr(j,0,n,js) / sfp_const
END DO
CALL intgr3(rhoc_vx,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),dpdot)
IF( xcpot%icorr .EQ. icorr_pbe0 ) THEN
......@@ -732,7 +730,6 @@ CONTAINS
nat = nat + atoms%neq(n)
END DO
! Skip over parts in cdnc not used
READ (17) ( rdum, n = 1,atoms%ntype )
END IF
!HF
......@@ -741,7 +738,6 @@ CONTAINS
!HF
IF ( xcpot%icorr == icorr_pbe0 .OR. xcpot%icorr == icorr_hse .OR.&
xcpot%icorr == icorr_hf .OR. xcpot%icorr == icorr_vhse ) THEN
CLOSE ( 17 )
DEALLOCATE( rhoc, rhoc_vx )
END IF
!HF end kinetic energy correction
......
......@@ -40,7 +40,8 @@
LOGICAL n_exist
! ..
! .. Local Arrays ..
REAL rhoc(DIMENSION%msh,atoms%ntype)
REAL rhoc(atoms%jmtd,atoms%ntype,dimension%jspd)
REAL tec(atoms%ntype,dimension%jspd),qintc(atoms%ntype,dimension%jspd)
COMPLEX :: cdom(1),cdomvz(1,1),cdomvxy(1,1,1)
COMPLEX, ALLOCATABLE :: qpw(:,:),rhtxy(:,:,:,:)
REAL , ALLOCATABLE :: rho(:,:,:,:),rht(:,:,:)
......@@ -56,15 +57,8 @@
ALLOCATE ( rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins),qpw(stars%ng3,input%jspins) )
ALLOCATE ( rhtxy(vacuum%nmzxy,stars%ng2-1,2,input%jspins),rht(vacuum%nmz,2,input%jspins) )
!
OPEN (17,file='cdnc',form='unformatted',status='old')
DO n = 1,atoms%ntype
READ (17) (rhoc(j,n),j=1,atoms%jri(n))
READ (17) dummy
ENDDO
CLOSE (17)
input%jspins=1
CALL readCoreDensity(input,atoms,dimension,rhoc,tec,qintc)
CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
CDN_INPUT_DEN_const,0,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
input%jspins=2
......@@ -77,15 +71,15 @@
!
DO n = 1,atoms%ntype
DO j = 1,atoms%jri(n)
rho(j,0,n,1) = rho(j,0,n,1) - rhoc(j,n)/sfp
rho(j,0,n,1) = rho(j,0,n,1) - rhoc(j,n,1)/sfp
ENDDO
! WRITE (16,FMT='(8f10.4)') (rho(i,0,n,1),i=1,16)
CALL intgr3(rho(1,0,n,1),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),qval)
p = (atoms%bmu(n)+sfp*qval)/ (2.*sfp*qval)
pp = 1. - p
DO j = 1,atoms%jri(n)
rho(j,0,n,jsp_new) = pp*rho(j,0,n,1) + rhoc(j,n)/ (2.*sfp)
rho(j,0,n,1) = p*rho(j,0,n,1) + rhoc(j,n)/ (2.*sfp)
rho(j,0,n,jsp_new) = pp*rho(j,0,n,1) + rhoc(j,n,1)/ (2.*sfp)
rho(j,0,n,1) = p*rho(j,0,n,1) + rhoc(j,n,1)/ (2.*sfp)
ENDDO
DO lh = 1,sphhar%nlh(atoms%ntypsy(na))
DO j = 1,atoms%jri(n)
......
......@@ -24,7 +24,7 @@
PRIVATE
PUBLIC plotdop
CONTAINS
SUBROUTINE plotdop(oneD,&
SUBROUTINE plotdop(oneD,dimension,&
& stars,vacuum,sphhar,atoms,&
& input,sym,cell,sliceplot,&
& l_noco,cdnfname)
......@@ -36,6 +36,7 @@
! ..
! .. Scalar Arguments ..
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_dimension),INTENT(IN):: dimension
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_sphhar),INTENT(IN) :: sphhar
......@@ -69,7 +70,7 @@
oldform = .false.
INQUIRE(file ="plotin",exist = oldform)
IF ( oldform ) THEN
CALL priv_old_plot(oneD,&
CALL priv_old_plot(oneD,dimension,&
& stars,vacuum,sphhar,atoms,&
& input,sym,cell,sliceplot,&
& l_noco,cdnfname)
......@@ -290,17 +291,19 @@
!------------------------------------------
! The old subroutine from Fleur is here
!------------------------------------------
SUBROUTINE priv_old_plot(oneD,&
SUBROUTINE priv_old_plot(oneD,dimension,&
& stars,vacuum,sphhar,atoms,&
& input,sym,cell,sliceplot,&
& l_noco,cdnfname)
!
USE m_outcdn
USE m_loddop
USE m_cdn_io
IMPLICIT NONE
! ..
! .. Scalar Arguments ..
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_dimension),INTENT(IN):: dimension
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_sphhar),INTENT(IN) :: sphhar
......@@ -317,7 +320,7 @@
!+odim
! ..
! .. Local Scalars ..
REAL rx,ry,s,sl,sm,su,x,xm,y,ym,tec,qint,sfp,xdnout
REAL rx,ry,s,sl,sm,su,x,xm,y,ym,sfp,xdnout
INTEGER i,i1,i2,i3,ii3,imshx,imshy,ix,iy,j,jsp,na,nfile,nplo,&
& nplot,nq,nt,nplott,jm,jspin,iter,ii1,ii2
LOGICAL twodim
......@@ -326,7 +329,8 @@
COMPLEX qpw(stars%n3d,input%jspins),rhtxy(vacuum%nmzxyd,stars%n2d-1,2,input%jspins)
REAL rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins),rht(vacuum%nmzd,2,input%jspins)
REAL ptp(3),rngl(3),rngm(3),rngu(3),tl(3),tm(3)
REAL tu(3),vx1(3),vx2(3),tl_r(3),tm_r(3),tu_r(3),rhocc(atoms%jmtd)
REAL tu(3),vx1(3),vx2(3),tl_r(3),tm_r(3),tu_r(3),rhocc(atoms%jmtd,atoms%ntype,dimension%jspd)
REAL tec(atoms%ntype,dimension%jspd),qint(atoms%ntype,dimension%jspd)
REAL pt(3),a3(3)
REAL, ALLOCATABLE :: cdn(:,:)
CHARACTER(len=10), ALLOCATABLE :: plotname(:)
......@@ -342,23 +346,16 @@
& iter,rho,qpw,rht,rhtxy)
!
IF (input%score) THEN
OPEN (17,file='cdnc',form='unformatted',status='old')
REWIND 17
!
CALL readCoreDensity(input,atoms,dimension,rhocc,tec,qint)
DO jspin = 1,input%jspins
DO nt = 1,atoms%ntype
jm = atoms%jri(nt)
READ (17) (rhocc(i),i=1,jm)
DO i = 1,atoms%jri(nt)
rho(i,0,nt,jspin) = rho(i,0,nt,jspin) - rhocc(i)/sfp
rho(i,0,nt,jspin) = rho(i,0,nt,jspin) - rhocc(i,nt,jspin)/sfp
ENDDO
READ (17) tec
ENDDO
READ (17) qint
qpw(1,jspin) = qpw(1,jspin) - qint/cell%volint
qpw(1,jspin) = qpw(1,jspin) - qint(nt,jspin)/cell%volint
ENDDO
!
CLOSE (17)
END IF
OPEN (18,file='plotin')
READ (18,FMT='(7x,l1)') twodim
......
......@@ -17,6 +17,7 @@ CONTAINS
USE m_fft3d
USE m_qpwtonmt
USE m_convol
USE m_cdn_io
IMPLICIT NONE
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_atoms),INTENT(IN) :: atoms
......@@ -37,34 +38,31 @@ CONTAINS
!locals
TYPE(t_atoms) :: atoms_tmp
REAL :: e_vdW
REAL,ALLOCATABLE :: n_grid(:),v_grid(:),rhc(:)
REAL,ALLOCATABLE :: n_grid(:),v_grid(:),rhc(:,:,:)
COMPLEX,ALLOCATABLE:: vpw(:),psq(:)
INTEGER :: n,ncmsh,j,i
LOGICAL :: l_core,l_pot
REAL tec(atoms%ntype,dimension%jspd),qintc(atoms%ntype,dimension%jspd)
l_core=.FALSE. !try to subtract core charge?
ALLOCATE(n_grid(27*stars%k1d*stars%k2d*stars%k3d),v_grid(27*stars%k1d*stars%k2d*stars%k3d))
ALLOCATE(vpw(SIZE(qpw)))
ALLOCATE(psq(SIZE(qpw)),rhc(DIMENSION%msh))
ALLOCATE(psq(SIZE(qpw)),rhc(atoms%jmtd,atoms%ntype,dimension%jspd))
IF (l_core) INQUIRE(file="cdnc",exist=l_core)
IF (l_core) l_core = isCoreDensityPresent()
IF (l_core) THEN
WRITE(6,*) "VdW contribution without core charge"
WRITE(6,*) "VdW contribution without core charge"
! read the core charge
OPEN (17,file='cdnc',form='unformatted',status='unknown')
CALL readCoreDensity(input,atoms,dimension,rhc,tec,qintc)
DO j=1,input%jspins
DO n=1,atoms%ntype
ncmsh = NINT( LOG( (atoms%rmt(n)+10.0)/atoms%rmsh(1,n) ) / atoms%dx(n) + 1 )
ncmsh = MIN( ncmsh, DIMENSION%msh )
READ(17) (rhc(i),i=1,ncmsh)
rho(:,1,n) = rho(:,1,n) - rhc(:SIZE(rho,1))/(4. * SQRT( ATAN (1.) ))
READ (17) !kinetic energy not needed
rho(:,1,n) = rho(:,1,n) - rhc(:SIZE(rho,1),n,j)/(4. * SQRT( ATAN (1.) ))
ENDDO
READ(17)
ENDDO
CLOSE(17)
ENDIF
! Construct the pseudo charge
......
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