Commit 289c481a authored by Matthias Redies's avatar Matthias Redies

integrate kEDs and laplace

parent 885c4f84
......@@ -23,13 +23,13 @@ CONTAINS
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_sphhar), INTENT(IN):: sphhar
TYPE(t_noco), INTENT(INOUT) :: noco
TYPE(t_lapl), INTENT(in) :: is_inte, mt_inte(:)
TYPE(t_grid), INTENT(in) :: is_inte, mt_inte(:)
REAL, INTENT(out) :: q(input%jspins), qis(input%jspins), qmt(atoms%ntype,input%jspins),&
qvac(2,input%jspins), qtot, qistot
TYPE(t_potden) :: integrand
TYPE(t_lapl) :: is_inte_mut
TYPE(t_grid) :: is_inte_mut
INTEGER :: n
call init_pw_grid(xcpot, stars, sym, cell)
......@@ -39,12 +39,13 @@ CONTAINS
call integrand%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN)
allocate(integrand%pw_w, mold=integrand%pw)
call init_mt_grid(input%jspins, atoms, sphhar, xcpot, sym)
!put is in potden-basis
call pw_from_grid(xcpot, stars,.True., is_inte_mut%lapl, integrand%pw, integrand%pw_w)
call pw_from_grid(xcpot, stars,.True., is_inte_mut%grid, integrand%pw, integrand%pw_w)
!put mt in potden-basis
call init_mt_grid(input%jspins, atoms, sphhar, xcpot, sym)
do n = 1,atoms%ntype
call mt_from_grid(atoms,sphhar,n,input%jspins,mt_inte(n)%lapl,integrand%mt(:,0:,n,:))
call mt_from_grid(atoms,sphhar,n,input%jspins,mt_inte(n)%grid,integrand%mt(:,0:,n,:))
enddo
! integrate my integrand
......@@ -77,7 +78,6 @@ CONTAINS
REAL :: q2(vacuum%nmz), w, rht1(vacuum%nmzd,2,input%jspins)
COMPLEX :: x(stars%ng3)
CALL timestart("cdntot")
qtot = 0.0
qistot = 0.0
DO jsp = 1,input%jspins
......@@ -159,6 +159,7 @@ CONTAINS
INTEGER, ALLOCATABLE :: lengths(:,:)
CHARACTER(LEN=20) :: attributes(6), names(6)
CALL timestart("cdntot")
call integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, den, &
q, qis, qmt, qvac, qtot, qistot)
......
......@@ -104,6 +104,8 @@ CONTAINS
#endif
call integrate_lapl(xcpot, stars, atoms, sym, vacuum, input, cell, oneD, sphhar,noco)
call integrate_kED_schr(xcpot, stars, atoms, sym, vacuum, input, cell, oneD, sphhar,noco)
call integrate_kED_libxc(xcpot, stars, atoms, sym, vacuum, input, cell, oneD, sphhar,noco)
END SUBROUTINE vgen
SUBROUTINE integrate_lapl(xcpot, stars, atoms, sym, vacuum, input, cell, oneD, sphhar,noco)
......@@ -128,6 +130,7 @@ CONTAINS
xcpot%is_lapl, xcpot%mt_lapl, &
q, qis, qmt, qvac, qtot, qistot)
write (*,*) "------------------------"
write (*,*) "lapl integration:"
write (*,*) "q = "
write (*,*) q
......@@ -139,9 +142,103 @@ CONTAINS
write (*,*) qvac
write (*,*) "qtot = "
write (*,*) qtot
write (*,*) "q = "
write (*,*) "qistot = "
write (*,*) qistot
write (*,*) "------------------------"
END SUBROUTINE integrate_lapl
SUBROUTINE integrate_kED_schr(xcpot, stars, atoms, sym, vacuum, input, cell, oneD, sphhar,noco)
use m_cdntot
use m_types
implicit none
CLASS(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_input),INTENT(IN) :: input
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_sphhar), INTENT(IN):: sphhar
TYPE(t_noco), INTENT(INOUT) :: noco
REAL :: q(input%jspins), qis(input%jspins), qmt(atoms%ntype,input%jspins),&
qvac(2,input%jspins), qtot, qistot
if(allocated(xcpot%is_kED_schr%grid) .and. allocated(xcpot%mt_kED_schr)) then
call integrate_grid(xcpot, stars, atoms, sym, vacuum, input, cell, oneD, sphhar,noco,&
xcpot%is_kED_schr, xcpot%mt_kED_schr, &
q, qis, qmt, qvac, qtot, qistot)
write (*,*) "------------------------"
write (*,*) "kED_schr integration:"
write (*,*) "q = "
write (*,*) q
write (*,*) "qis = "
write (*,*) qis
write (*,*) "qmt = "
write (*,*) qmt
write (*,*) "qvac = "
write (*,*) qvac
write (*,*) "qtot = "
write (*,*) qtot
write (*,*) "qistot = "
write (*,*) qistot
write (*,*) "------------------------"
endif
END SUBROUTINE integrate_kED_schr
SUBROUTINE integrate_kED_libxc(xcpot, stars, atoms, sym, vacuum, input, cell, oneD, sphhar,noco)
use m_cdntot
use m_types
implicit none
CLASS(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_input),INTENT(IN) :: input
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_sphhar), INTENT(IN):: sphhar
TYPE(t_noco), INTENT(INOUT) :: noco
REAL :: q(input%jspins), qis(input%jspins), qmt(atoms%ntype,input%jspins),&
qvac(2,input%jspins), qtot, qistot
TYPE(t_grid) :: is_kED_libxc
TYPE(t_grid), allocatable :: mt_kED_libxc(:)
INTEGER :: i
if(allocated(xcpot%is_kED_schr%grid) .and. allocated(xcpot%mt_kED_schr)) then
allocate(mt_kED_libxc(size(xcpot%mt_kED_schr)))
is_kED_libXC%grid = xcpot%is_kED_schr%grid + 0.25 * xcpot%is_lapl%grid
do i = 1,atoms%ntype
mt_kED_libXC(i)%grid = xcpot%mt_kED_schr(i)%grid + 0.25 * xcpot%mt_lapl(i)%grid
enddo
call integrate_grid(xcpot, stars, atoms, sym, vacuum, input, cell, oneD, sphhar,noco,&
is_kED_libXC, mt_kED_libXC, &
q, qis, qmt, qvac, qtot, qistot)
write (*,*) "------------------------"
write (*,*) "kED_libXC integration:"
write (*,*) "q = "
write (*,*) q
write (*,*) "qis = "
write (*,*) qis
write (*,*) "qmt = "
write (*,*) qmt
write (*,*) "qvac = "
write (*,*) qvac
write (*,*) "qtot = "
write (*,*) qtot
write (*,*) "qistot = "
write (*,*) qistot
write (*,*) "------------------------"
endif
END SUBROUTINE integrate_kED_libxc
END MODULE m_vgen
......@@ -15,17 +15,17 @@ MODULE m_types_xcpot
use m_types_potden
IMPLICIT NONE
PRIVATE
PUBLIC :: t_xcpot,t_gradients, t_lapl
PUBLIC :: t_xcpot,t_gradients, t_grid
type t_lapl
real, allocatable :: lapl(:,:)
end type t_lapl
type t_grid
real, allocatable :: grid(:,:)
end type t_grid
TYPE,ABSTRACT :: t_xcpot
REAL :: gmaxxc
TYPE(t_potden) :: comparison_kinED_pw(3)
TYPE(t_lapl) :: is_lapl
TYPE(t_lapl), allocatable :: mt_lapl(:)
TYPE(t_grid) :: is_lapl, is_kED_schr
TYPE(t_grid), allocatable :: mt_lapl(:), mt_kED_schr(:)
CONTAINS
PROCEDURE :: vxc_is_LDA=>xcpot_vxc_is_LDA
PROCEDURE :: exc_is_LDA=>xcpot_exc_is_LDA
......
......@@ -27,22 +27,24 @@ CONTAINS
! generate nspd points on a sherical shell with radius 1.0
! angular mesh equidistant in phi,
! theta are zeros of the legendre polynomials
ALLOCATE (wt(atoms%nsp()), rx(3, atoms%nsp()), thet(atoms%nsp()))
CALL gaussp(atoms%lmaxd, rx, wt)
! generate the lattice harmonics on the angular mesh
ALLOCATE (ylh(atoms%nsp(), 0:sphhar%nlhd, sphhar%ntypsd))
IF (xcpot%needs_grad()) THEN
ALLOCATE (ylht, MOLD=ylh)
ALLOCATE (ylhtt, MOLD=ylh)
ALLOCATE (ylhf, MOLD=ylh)
ALLOCATE (ylhff, MOLD=ylh)
ALLOCATE (ylhtf, MOLD=ylh)
CALL lhglptg(sphhar, atoms, rx, atoms%nsp(), xcpot, sym, &
ylh, thet, ylht, ylhtt, ylhf, ylhff, ylhtf)
ELSE
CALL lhglpts(sphhar, atoms, rx, atoms%nsp(), sym, ylh)
END IF
if(.not. allocated(wt)) then
ALLOCATE (wt(atoms%nsp()), rx(3, atoms%nsp()), thet(atoms%nsp()))
CALL gaussp(atoms%lmaxd, rx, wt)
! generate the lattice harmonics on the angular mesh
ALLOCATE (ylh(atoms%nsp(), 0:sphhar%nlhd, sphhar%ntypsd))
IF (xcpot%needs_grad()) THEN
ALLOCATE (ylht, MOLD=ylh)
ALLOCATE (ylhtt, MOLD=ylh)
ALLOCATE (ylhf, MOLD=ylh)
ALLOCATE (ylhff, MOLD=ylh)
ALLOCATE (ylhtf, MOLD=ylh)
CALL lhglptg(sphhar, atoms, rx, atoms%nsp(), xcpot, sym, &
ylh, thet, ylht, ylhtt, ylhf, ylhff, ylhtf)
ELSE
CALL lhglpts(sphhar, atoms, rx, atoms%nsp(), sym, ylh)
END IF
ENDIF
END SUBROUTINE init_mt_grid
SUBROUTINE mt_to_grid(xcpot, jspins, atoms, sphhar, den_mt, n, grad, ch)
......
......@@ -58,8 +58,7 @@ CONTAINS
!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)
write (*,*) "set is_grad"
xcpot%is_lapl%lapl = grad%laplace
xcpot%is_lapl%grid = grad%laplace
ALLOCATE(v_xc,mold=rho)
ALLOCATE(v_x,mold=rho)
......@@ -108,6 +107,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)
xcpot%is_kED_schr%grid = kinED_rs
call save_npy("exc_pw.npy", e_xc(:,1))
ELSE
CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad)
......
......@@ -93,6 +93,7 @@ CONTAINS
ALLOCATE(ED_rs, mold=ch)
ALLOCATE(vTot_rs, mold=ch)
ALLOCATE(kinED_RS, mold=ch)
if(.not. allocated(xcpot%mt_kED_schr)) allocate(xcpot%mt_kED_schr(atoms%ntype))
ENDIF
CALL init_mt_grid(input%jspins,atoms,sphhar,xcpot,sym)
......@@ -114,7 +115,7 @@ CONTAINS
DO n = n_start,atoms%ntype,n_stride
CALL mt_to_grid(xcpot, input%jspins, atoms,sphhar,den%mt(:,0:,n,:),n,grad,ch)
write (*,*) "set grad for mt = ", n
xcpot%mt_lapl(n)%lapl = grad%laplace
xcpot%mt_lapl(n)%grid = grad%laplace
!
! calculate the ex.-cor. potential
......@@ -164,6 +165,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)
xcpot%mt_kED_schr(n)%grid = 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)
......
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