Commit 5cafd10b authored by Matthias Redies's avatar Matthias Redies

pre monday merge

parent 1ddd915b
......@@ -21,4 +21,7 @@ juDFT/string.f90
juDFT/time.F90
juDFT/args.F90
juDFT/sysinfo.F90
# just for debugging purposes
juDFT/npy.F90
juDFT/endian_swap.f90
)
......@@ -230,6 +230,7 @@ subroutine save_kinED(xcpot, input, noco, stars, cell, sym)
use m_pw_tofrom_grid
use m_judft_stop
use m_metagga, only: give_stats
use m_npy
implicit none
CLASS(t_xcpot),INTENT(IN) :: xcpot
......@@ -247,10 +248,6 @@ subroutine save_kinED(xcpot, input, noco, stars, cell, sym)
do dim_idx = 1,3
call pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, cell, &
xcpot%comparison_kinED_pw(dim_idx)%pw, grad, tmp)
if(any(tmp > 1e5)) then
write (*,*) "tmp > 1e5 for", dim_idx
call give_stats(tmp, "tmp(dim_idx)")
endif
if(.not. allocated(kinED)) then
allocate(kinED, mold=tmp)
kinED = 0.0
......@@ -276,6 +273,7 @@ subroutine save_kinED(xcpot, input, noco, stars, cell, sym)
call finish_pw_grid()
write (*,*) "kED shape =", shape(kinED)
call save_npy("kin_ED_pwway.npy", kinED)
open(unit=69, file="kin_ED_pwway.dat")
write (69,'(ES17.10)') kinED
close(69)
......
......@@ -111,6 +111,7 @@ CONTAINS
ALLOCATE( mx(0:ifftd-1),my(0:ifftd-1),magmom(0:ifftd-1))
ENDIF
END IF
IF (PRESENT(rho)) THEN
!Put den_pw on grid and store into rho(:,1:2)
DO js=1,jspins
......
......@@ -167,7 +167,8 @@ 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)
vCoul%pw(1,ispin) = cmplx(-2.0,0.0) ! that's my line
write (*,*) "shifting vtot by", vCoul%pw(1,ispin)
vCoul%pw(2:stars%ng3,ispin) = fpi_const * psq(2:stars%ng3) / stars%sk3(2:stars%ng3) ** 2
end if
end if
......
......@@ -10,6 +10,7 @@
!--------------------------------------------------------------------------------
MODULE m_vis_xc
USE m_juDFT
use m_convol
! ******************************************************
! subroutine generates the exchange-correlation potential
! in the interstitial region c.l.fu
......@@ -34,6 +35,7 @@ CONTAINS
USE m_types_xcpot_libxc
USE m_libxc_postprocess_gga
USE m_metagga
USE m_npy
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
......@@ -42,18 +44,20 @@ CONTAINS
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_potden),INTENT(IN) :: den, EnergyDen
TYPE(t_potden),INTENT(IN) :: den, EnergyDen
TYPE(t_potden),INTENT(INOUT) :: vTot,vx,exc
TYPE(t_gradients) :: grad, tmp_grad
REAL, ALLOCATABLE :: rho(:,:), ED_rs(:,:), vTot_rs(:,:), kinED_rs(:,:)
REAL, ALLOCATABLE :: rho_conv(:,:), ED_conv(:,:), vTot_conv(:,:)
COMPLEX, ALLOCATABLE :: den_pw_w(:,:), EnergyDen_pw_w(:,:)
REAL, ALLOCATABLE :: v_x(:,:),v_xc(:,:),e_xc(:,:)
INTEGER :: jspin
CALL init_pw_grid(xcpot,stars,sym,cell)
!Put the charge on the grid, in GGA case also calculate gradients
CALL pw_to_grid(xcpot,input%jspins,noco%l_noco,stars,cell,den%pw,grad,rho)
call give_stats(rho, "rho")
ALLOCATE(v_xc,mold=rho)
ALLOCATE(v_x,mold=rho)
......@@ -71,11 +75,32 @@ CONTAINS
! use updated vTot for exc calculation
IF(ALLOCATED(EnergyDen%pw) .AND. xcpot%exc_is_MetaGGA()) THEN
allocate(den_pw_w, mold=vTot%pw_w)
allocate(EnergyDen_pw_w, mold=vTot%pw_w)
do jspin = 1,input%jspins
call convol(stars, vTot%pw_w(:,jspin), vTot%pw, stars%ufft)
call convol(stars, den_pw_w(:,jspin), den%pw, stars%ufft)
call convol(stars, EnergyDen_pw_w(:,jspin), EnergyDen%pw, stars%ufft)
enddo
CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
cell, vTot%pw_w, tmp_grad, vTot_conv)
CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
cell, den_pw_w, tmp_grad, rho_conv)
CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
cell, EnergyDen_pw_w, tmp_grad, ED_conv)
call save_npy("den_conv.npy", rho_conv)
call save_npy("EnergyDen_conv.npy", ED_conv)
call save_npy("vTot_conv.npy", vTot_conv)
call save_npy("kinED_schr_conv_precut.npy", ED_conv - vTot_conv*rho_conv)
CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
cell, EnergyDen%pw, tmp_grad, ED_rs)
CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
cell, vTot%pw, tmp_grad, vTot_rs)
CALL calc_kinEnergyDen(ED_rs, vTot_rs, rho, kinED_rs)
CALL calc_kinEnergyDen(ED_rs, vTot_rs, rho, kinED_rs, .True.)
ENDIF
!calculate the ex.-cor energy density
......
......@@ -148,7 +148,7 @@ CONTAINS
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)
CALL calc_kinEnergyDen(ED_rs, vTot_rs, ch, kinED_rs, .False.)
ENDIF
IF (ALLOCATED(exc%mt)) THEN
......
......@@ -12,36 +12,35 @@ MODULE m_metagga
end type t_RS_potden
CONTAINS
SUBROUTINE calc_kinEnergyDen(EnergyDen_rs, vTot_rs, den_rs, kinEnergyDen_RS)
SUBROUTINE calc_kinEnergyDen(EnergyDen_rs, vTot_rs, den_rs, kinEnergyDen_RS, is_pw)
USE m_juDFT_stop
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
#ifdef CPP_LIBXC
REAL, PARAMETER :: eps = 1e-15
!implicit allocation
kinEnergyDen_RS = EnergyDen_RS - vTot_RS * den_RS
if(all(shape(kinEnergyDen_RS) == [6144,1])) then
call give_stats(EnergyDen_RS, array_name="ED_RS")
call give_stats(vTot_RS, array_name="vTot_RS")
call give_stats(den_RS, array_name="den_RS")
call give_stats(kinEnergyDen_RS, array_name="kin_ED_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
write (*,*) ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
if(any(kinEnergyDen_RS < eps)) then
write (6,*) " lowest kinetic energy density cutoff = ", minval(kinEnergyDen_RS)
kinEnergyDen_RS = max(kinEnergyDen_RS, eps)
endif
if(all(shape(kinEnergyDen_RS) == [6144,1])) then
write (*,*) "kinED shape:", shape(kinEnergyDen_RS)
write (*,*) "write old"
open(unit=69, file="kinED_pw_schroeway.dat")
write (69,'(ES17.10)') kinEnergyDen_RS
close(69)
if(is_pw) then
call save_npy("kinED_pw_schroeway.npy",kinEnergyDen_RS)
write (*,*) "read new"
open(unit=69, file="kin_ED_pwway.dat")
......@@ -66,6 +65,17 @@ CONTAINS
write (*,*) "#############################"
END SUBROUTINE give_stats
SUBROUTINE dump_array(array, filename)
implicit none
real, intent(in) :: array(:,:)
character(len=*), intent(in) :: filename
open(69, file=filename)
write (69,'(ES17.10)') array
close(69)
END SUBROUTINE
SUBROUTINE calc_EnergyDen(eig_id, mpi, kpts, noco, input, banddos, cell, atoms, enpara, stars, &
vacuum, DIMENSION, sphhar, sym, vTot, oneD, results, EnergyDen)
......
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