Commit 37c58a36 authored by Matthias Redies's avatar Matthias Redies

pre-dev commit

parent 3b4872d9
......@@ -21,4 +21,27 @@ contains
output = input(front:back)
end function strip
function str2int(str) result(int)
implicit none
character(len=*), intent(in) :: str
integer :: int
integer :: stat
read (str, *, iostat=stat) int
if (stat /= 0) then
write (*, *) "str reading failed", str
stop 9
endif
end function str2int
function int2str(num) result(ret_str)
implicit none
integer, intent(in) :: num
character(len=:), allocatable :: ret_str
allocate(character(100) :: ret_str)
write (ret_str,*) num
ret_str = strip(ret_str)
end function int2str
end module m_juDFT_string
......@@ -91,25 +91,25 @@ CONTAINS
eigSum = results%seigscv + results%seigc
results%tote = eigSum
WRITE (6,FMT=8010) results%tote
8010 FORMAT (/,10x,'sum of eigenvalues =',t40,f20.10)
8010 FORMAT (/,10x,'sum of eigenvalues =',t40,ES20.10)
!
! ---> add contribution of coulomb potential
!
results%tote = results%tote + 0.5e0*results%te_vcoul
WRITE (6,FMT=8020) results%te_vcoul
8020 FORMAT (/,10x,'density-coulomb potential integral =',t40,f20.10)
8020 FORMAT (/,10x,'density-coulomb potential integral =',t40,ES20.10)
!
! ---> add contribution of effective potential
!
results%tote = results%tote - results%te_veff
WRITE (6,FMT=8030) results%te_veff
8030 FORMAT (/,10x,'density-effective potential integral =',t40,f20.10)
8030 FORMAT (/,10x,'density-effective potential integral =',t40,ES20.10)
!
! ---> add contribution of exchange-correlation energy
!
results%tote = results%tote + results%te_exc
WRITE (6,FMT=8040) results%te_exc
8040 FORMAT (/,10x,'charge density-ex.-corr.energy density integral=', t40,f20.10)
8040 FORMAT (/,10x,'charge density-ex.-corr.energy density integral=', t40,ES20.10)
!
! ---> Fock exchange contribution
!
......@@ -122,8 +122,8 @@ CONTAINS
ENDIF
WRITE (6,FMT=8100) 0.5e0*results%te_hfex%valence
WRITE (6,FMT=8101) 0.5e0*results%te_hfex%core
8100 FORMAT (/,10x,'Fock-exchange energy (valence)=',t40,f20.10)
8101 FORMAT (10x,'Fock-exchange energy (core)=',t40,f20.10)
8100 FORMAT (/,10x,'Fock-exchange energy (valence)=',t40,ES20.10)
8101 FORMAT (10x,'Fock-exchange energy (core)=',t40,ES20.10)
! ----> VM terms
......@@ -201,7 +201,7 @@ CONTAINS
WRITE ( 6,FMT=8081) results%tote-0.5e0*results%ts
END IF
WRITE(attributes(1),'(f20.10)') results%tote
WRITE(attributes(1),'(ES20.10)') results%tote
WRITE(attributes(2),'(a)') 'Htr'
WRITE(attributes(3),'(a)') 'HF'
IF (hybrid%l_calhf) THEN
......@@ -231,19 +231,19 @@ CONTAINS
CALL writeXMLElementFormPoly('freeEnergy',(/'value'/),(/results%tote-results%ts/),reshape((/38,20/),(/1,2/)))
CALL writeXMLElementFormPoly('extrapolationTo0K',(/'value'/),(/results%tote-0.5e0*results%ts/),reshape((/31,20/),(/1,2/)))
CALL closeXMLElement('totalEnergy')
8060 FORMAT (/,/,' ----> input%total energy=',t40,f20.10,' htr')
8061 FORMAT (/,/,' ----> HF input%total energy=',t40,f20.10,' htr')
8050 FORMAT (/,10x,'Madelung term for atom type:',i3,t40,f20.10)
8045 FORMAT (/,10x,'el.-nucl. inter. diff. m.t.',t40,f20.10)
8065 FORMAT (/,/,' ----> (input%tkb*entropy) TS=',t40,f20.10,' htr')
8066 FORMAT (/,/,' ----> HF (input%tkb*entropy) TS=',t40,f20.10,' htr')
8070 FORMAT (/,/,' ----> free energy=',t40,f20.10,' htr')
8071 FORMAT (/,/,' ----> HF free energy=',t40,f20.10,' htr')
8060 FORMAT (/,/,' ----> input%total energy=',t40,ES20.10,' htr')
8061 FORMAT (/,/,' ----> HF input%total energy=',t40,ES20.10,' htr')
8050 FORMAT (/,10x,'Madelung term for atom type:',i3,t40,ES20.10)
8045 FORMAT (/,10x,'el.-nucl. inter. diff. m.t.',t40,ES20.10)
8065 FORMAT (/,/,' ----> (input%tkb*entropy) TS=',t40,ES20.10,' htr')
8066 FORMAT (/,/,' ----> HF (input%tkb*entropy) TS=',t40,ES20.10,' htr')
8070 FORMAT (/,/,' ----> free energy=',t40,ES20.10,' htr')
8071 FORMAT (/,/,' ----> HF free energy=',t40,ES20.10,' htr')
8080 FORMAT (/,/,' extrapolation for T->0',&
/,' ----> input%total electron energy=',t40,f20.10,' htr')
/,' ----> input%total electron energy=',t40,ES20.10,' htr')
8081 FORMAT (/,/,' extrapolation for T->0',&
/,' ----> HF input%total electron energy=',t40,f20.10,' htr')
8090 FORMAT (/,/,' ----> correction for lda+U =',t40,f20.10,' htr')
/,' ----> HF input%total electron energy=',t40,ES20.10,' htr')
8090 FORMAT (/,/,' ----> correction for lda+U =',t40,ES20.10,' htr')
END SUBROUTINE totale
END MODULE m_totale
......@@ -261,6 +261,7 @@ CONTAINS
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinEnergyDen_KS)
use m_npy
use m_constants
IMPLICIT NONE
CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
INTEGER, INTENT (IN) :: jspins
......@@ -282,7 +283,7 @@ CONTAINS
! tau = 0.5 * sum[|grad phi_i(r)|²]
! see eq (3) in https://doi.org/10.1063/1.1565316
REAL, ALLOCATABLE :: kinEnergyDen_libXC(:,:)
REAL, ALLOCATABLE :: kinEnergyDen_libXC(:,:), pkzb_ratio(:,:), pkzb_zaehler(:,:), pkzb_nenner(:,:)
IF (xcpot%exc_is_gga()) THEN
IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a GGA potential without providing derivatives")
......@@ -300,13 +301,40 @@ CONTAINS
ELSEIF(xcpot%exc_is_MetaGGA()) THEN
IF(PRESENT(kinEnergyDen_KS)) THEN
! apply correction in eq (4) in https://doi.org/10.1063/1.1565316
kinEnergyDen_libXC = transpose(kinEnergyDen_KS + 0.25 * grad%laplace)
if(any(kinEnergyDen_libXC < 0.0)) write (*,*) "kED negative", &
minval(kinEnergyDen_libxc)
write (filename, '("kED_libxc_", I0.6, ".npy")') size(kinEnergyDen_libxc, dim=2)
!kinEnergyDen_libXC = transpose(kinEnergyDen_KS + 0.25 * grad%laplace)
!where(kinEnergyDen_libXC < 1d-5) kinEnergyDen_libXC = 1d-5
write (*,*) "apply tf approx. shapes: "
write (*,*) "shape(rh) = ", shape(rh)
write (*,*) "shape(grad%sigma) = ", shape(transpose(grad%sigma))
write (*,*) "shape(grad%lapl) = ", shape(grad%laplace)
kinEnergyDen_libXC = 0.3 * (3.0*pi_const**2)**(2./3.) * rh**(5./3.) &
+ 1.0/72.0 * transpose(abs(grad%sigma))/rh &
+ 1.0/6.0 * grad%laplace
pkzb_zaehler = (1./8. * transpose(abs(grad%sigma))/rh)**2
pkzb_nenner = kinEnergyDen_libxc**2
pkzb_ratio = pkzb_zaehler/pkzb_nenner
write (*,*) "pkzb ratio:"
write (*,*) "min = ", minval(pkzb_ratio)
write (*,*) "max = ", maxval(pkzb_ratio)
write (filename, '("kED_libxc_", I0.6, ".npy")') size(kinEnergyDen_libxc, dim=1)
call save_npy(filename, transpose(kinEnergyDen_libxc))
write (filename, '("sigma_", I0.6, ".npy")') size(grad%sigma, dim=2)
call save_npy(filename, grad%sigma)
write (filename, '("pkzb_zaehler_", I0.6, ".npy")') size(kinEnergyDen_libxc, dim=1)
call save_npy(filename, transpose(pkzb_zaehler))
write (filename, '("pkzb_nenner_", I0.6, ".npy")') size(kinEnergyDen_libxc, dim=1)
call save_npy(filename, transpose(pkzb_nenner))
write (filename, '("pkzb_ratio_", I0.6, ".npy")') size(kinEnergyDen_libxc, dim=1)
call save_npy(filename, transpose(pkzb_ratio))
exc = 0.0
excc = 0.0
call xc_f03_mgga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), grad%sigma, &
transpose(grad%laplace), kinEnergyDen_libXC, exc)
......
......@@ -44,7 +44,7 @@ CONTAINS
tis = cell%omtil * REAL( DOT_PRODUCT(vpot%pw_w(:stars%ng3,ispin),den%pw(:stars%ng3,ispin)))
WRITE (6,FMT=8020) tis
8020 FORMAT (/,10x,'interstitial :',t40,f20.10)
8020 FORMAT (/,10x,'interstitial :',t40,ES20.10)
RESULT = RESULT + tis
!
......@@ -63,7 +63,7 @@ CONTAINS
nat = nat + atoms%neq(n)
ENDDO
WRITE (6,FMT=8030) tmt
8030 FORMAT (/,10x,'muffin tin spheres :',t40,f20.10)
8030 FORMAT (/,10x,'muffin tin spheres :',t40,ES20.10)
RESULT = RESULT + tmt
!
! *********** VACUUM REGION**************
......
......@@ -167,8 +167,7 @@ contains
! if( abs( real( psq(1) ) ) * cell%omtil < 0.01 ) vCoul%pw(1,ispin) = 0.0
! there is a better option now using qfix in mix
else
vCoul%pw(1,ispin) = cmplx(0.0,0.0) ! that's my line
write (*,*) "shifting vtot by", vCoul%pw(1,ispin)
vCoul%pw(1,ispin) = cmplx(0.0,0.0)
vCoul%pw(2:stars%ng3,ispin) = fpi_const * psq(2:stars%ng3) / stars%sk3(2:stars%ng3) ** 2
end if
end if
......
......@@ -105,6 +105,7 @@ CONTAINS
IF(ALLOCATED(EnergyDen%pw) .AND. xcpot%exc_is_MetaGGA()) THEN
CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad, kinED_rs)
call save_npy("exc_pw.npy", e_xc(:,1))
ELSE
CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad)
ENDIF
......
......@@ -24,6 +24,7 @@ MODULE m_vmt_xc
CONTAINS
SUBROUTINE vmt_xc(DIMENSION,mpi,sphhar,atoms,&
den,xcpot,input,sym, obsolete,EnergyDen,vTot,vx,exc)
use m_npy
#include"cpp_double.h"
use m_libxc_postprocess_gga
......@@ -31,6 +32,7 @@ CONTAINS
USE m_types_xcpot_inbuild
USE m_types
USE m_metagga
USE m_juDFT_string
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
......@@ -106,7 +108,7 @@ CONTAINS
n_start=1
n_stride=1
#endif
call save_npy("rmsh.npy", atoms%rmsh)
DO n = n_start,atoms%ntype,n_stride
CALL mt_to_grid(xcpot, input%jspins, atoms,sphhar,den%mt(:,0:,n,:),nsp,n,grad,ch)
!
......@@ -129,7 +131,7 @@ CONTAINS
ENDIF
!Add postprocessing for libxc
IF (l_libxc.AND.xcpot%needs_grad()) CALL libxc_postprocess_gga_mt(xcpot,atoms,sphhar,n,v_xc,grad)
IF (l_libxc.AND.xcpot%needs_grad()) CALL libxc_postprocess_gga_mt(xcpot,atoms,sphhar,n,v_xc,grad, atom_num=n)
CALL mt_from_grid(atoms,sphhar,nsp,n,input%jspins,v_xc,vTot%mt(:,0:,n,:))
CALL mt_from_grid(atoms,sphhar,nsp,n,input%jspins,v_x,vx%mt(:,0:,n,:))
......@@ -147,8 +149,7 @@ CONTAINS
ENDDO
CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, vTot_tmp%mt(:,0:,n,:), &
nsp, n, tmp_grad, vTot_rs)
CALL calc_kinEnergyDen(ED_rs, vTot_rs, ch, kinED_rs, .False.)
CALL calc_kinEnergyDen(ED_rs, vTot_rs, ch, kinED_rs, is_pw=.False., nsp=nsp, atm_idx=n)
ENDIF
IF (ALLOCATED(exc%mt)) THEN
......@@ -158,6 +159,7 @@ CONTAINS
IF(perform_MetaGGA) THEN
CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),e_xc(:nsp*atoms%jri(n),1),grad, kinED_rs)
call save_npy("exc_mt.npy", e_xc(:,1))
ELSE
CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),e_xc(:nsp*atoms%jri(n),1),grad)
ENDIF
......@@ -187,18 +189,4 @@ CONTAINS
!
RETURN
END SUBROUTINE vmt_xc
function get_radial_line(den, line_idx, nsp) result(line)
implicit none
real, intent(in) :: den(:, :)
integer, intent(in) :: line_idx, nsp
real, allocatable :: line(:, :)
integer :: num_elem
write (*, *) "size(den(:,1)) =", size(den(:, 1))
write (*, *) "nsp =", nsp
write (*, *) "size(den(:,1))/nsp =", size(den(:, 1))/nsp
end function
END MODULE m_vmt_xc
......@@ -6,9 +6,12 @@
MODULE m_libxc_postprocess_gga
CONTAINS
SUBROUTINE libxc_postprocess_gga_mt(xcpot,atoms,sphhar,n,v_xc,grad)
SUBROUTINE libxc_postprocess_gga_mt(xcpot,atoms,sphhar,n,v_xc,grad, atom_num)
USE m_mt_tofrom_grid
USE m_types
USE m_npy
use m_judft_string
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_atoms),INTENT(IN) :: atoms
......@@ -16,10 +19,12 @@ CONTAINS
INTEGER,INTENT(IN) :: n
REAL,INTENT(INOUT) :: v_xc(:,:)
TYPE(t_gradients),INTENT(IN):: grad
INTEGER, OPTIONAL :: atom_num
INTEGER :: nsp,n_sigma,i
REAL,ALLOCATABLE:: vsigma(:,:),vsigma_mt(:,:,:)
TYPE(t_gradients)::grad_vsigma
character(len=:), allocatable :: fname
n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !Number of contracted gradients in libxc 1 for non-spin-polarized, 3 otherwise
nsp=SIZE(v_xc,1) !no of points
......@@ -33,6 +38,11 @@ CONTAINS
ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma))
CALL mt_to_grid(xcpot,n_sigma,atoms,sphhar,vsigma_mt,nsp/atoms%jmtd,n,grad=grad_vsigma)
fname = merge("mt=" //int2str(atom_num) // "_lapl.npy","mt_lapl.npy", present(atom_num))
call save_npy(fname, grad%laplace)
fname = merge("mt="// int2str(atom_num) // "_grad.npy", "mt_grad.npy", present(atom_num))
call save_npy(fname, grad%gr)
CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc)
END SUBROUTINE libxc_postprocess_gga_mt
......@@ -62,8 +72,8 @@ CONTAINS
ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma))
CALL pw_to_grid(xcpot,n_sigma,.false.,stars,cell,vsigma_g,grad_vsigma)
call save_npy("lapl.npy", grad%laplace)
call save_npy("grad.npy", grad%gr)
call save_npy("pw_lapl.npy", grad%laplace)
call save_npy("pw_grad.npy", grad%gr)
CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc)
END SUBROUTINE libxc_postprocess_gga_pw
......
......@@ -12,14 +12,16 @@ MODULE m_metagga
end type t_RS_potden
CONTAINS
SUBROUTINE calc_kinEnergyDen(EnergyDen_rs, vTot_rs, den_rs, kinEnergyDen_RS, is_pw)
SUBROUTINE calc_kinEnergyDen(EnergyDen_rs, vTot_rs, den_rs, kinEnergyDen_RS, is_pw, nsp, atm_idx)
USE m_juDFT_stop
USE m_juDFT_string
use m_npy
!use m_cdngen
IMPLICIT NONE
REAL, INTENT(in) :: den_RS(:,:), EnergyDen_RS(:,:), vTot_RS(:,:)
REAL, INTENT(inout), allocatable :: kinEnergyDen_RS(:,:)
LOGICAL, INTENT(IN), optional :: is_pw
INTEGER, INTENT(IN), optional :: nsp, atm_idx
#ifdef CPP_LIBXC
REAL, PARAMETER :: eps = 1e-15
......@@ -27,20 +29,27 @@ CONTAINS
kinEnergyDen_RS = EnergyDen_RS - vTot_RS * den_RS
if(is_pw) then
call save_npy("EnergyDen_RS.npy", EnergyDen_RS)
call save_npy("vTot_RS.npy", vTot_RS)
call save_npy("den_RS.npy", den_RS)
call save_npy("vTot_RSxdenRS.npy", vtot_RS*den_RS)
call save_npy("kinED_pw_schroeway_precut.npy",kinEnergyDen_RS)
endif
if(is_pw) then
call save_npy("kinED_pw_schroeway.npy",kinEnergyDen_RS)
!write (*,*) "read new"
!open(unit=69, file="kin_ED_pwway.dat")
!read(69,*) kinEnergyDen_RS
!close(69)
call save_npy("pw_EnergyDen_RS.npy", EnergyDen_RS)
call save_npy("pw_vTot_RS.npy", vTot_RS)
call save_npy("pw_den_RS.npy", den_RS)
call save_npy("pw_vTot_RSxdenRS.npy", vtot_RS*den_RS)
call save_npy("pw_kinED_pw_schroeway.npy",kinEnergyDen_RS)
else
if(present(nsp) .and. present(atm_idx)) then
!call save_npy("mt=" // int2str(atm_idx) // "_EnergyDen_RS.npy", &
!get_radial_line(EnergyDen_RS,1,nsp))
!call save_npy("mt=" // int2str(atm_idx) // "_vTot_RS.npy", &
!get_radial_line(vTot_RS,1,nsp))
!call save_npy("mt=" // int2str(atm_idx) // "_den_RS.npy", &
!get_radial_line(den_RS,1,nsp))
!call save_npy("mt=" // int2str(atm_idx) // "_vTot_RSxdenRS.npy", &
!get_radial_line(vtot_RS*den_RS,1,nsp))
!call save_npy("mt=" // int2str(atm_idx) // "_kinED_schroeway.npy", &
!get_radial_line(kinEnergyDen_RS,1,nsp))
call save_npy("mt=" // int2str(atm_idx) // "_den_RS.npy", den_RS)
else
write (*,*) "something not present"
endif
endif
#else
CALL juDFT_error("MetaGGA require LibXC",hint="compile Fleur with LibXC (e.g. by giving '-external libxc' to ./configure")
......@@ -223,7 +232,7 @@ CONTAINS
INTEGER :: ikpt,jsp_start,jsp_end,ispin,jsp
INTEGER :: iErr,nbands,noccbd,iType
INTEGER :: skip_t,skip_tt,nStart,nEnd,nbasfcn
LOGICAL :: l_orbcomprot, l_real, l_dosNdir, l_corespec
LOGICAL :: l_orbcomprot, l_real, l_dosNdir
! Local Arrays
REAL, ALLOCATABLE :: we(:)
......@@ -279,8 +288,6 @@ CONTAINS
END DO
IF (noco%l_mperp) CALL denCoeffsOffdiag%addRadFunScalarProducts(atoms,f,g,flo,iType)
IF (banddos%l_mcd) CALL mcd_init(atoms,input,dimension,vTot%mt(:,0,:,:),g,f,mcd,iType,jspin)
IF (l_coreSpec) CALL corespec_rme(atoms,input,iType,dimension%nstd,input%jspins,jspin,results%ef,&
dimension%msh,vTot%mt(:,0,:,:),f,g)
END DO
DEALLOCATE (f,g,flo)
......@@ -381,4 +388,23 @@ CONTAINS
res = matmul(transpose(cell%bmat), vec)
end function internal_to_rez
function get_radial_line(den, line_idx, nsp) result(line)
use m_npy
implicit none
real, intent(in) :: den(:, :)
integer, intent(in) :: line_idx, nsp
real, allocatable :: line(:, :)
integer :: num_elem
integer, allocatable :: mask(:)
allocate(mask(size(den, dim=1)))
mask = 0
mask(line_idx:size(mask):nsp) = 1
call save_npy("line_mask.npy", mask)
num_elem = size(den(:, 1))/nsp
line = den(line_idx:size(den(:,1)):nsp,:)
end function
END MODULE m_metagga
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