Commit c30b1989 authored by Matthias Redies's avatar Matthias Redies

save lapl for integration

parent 62bc4f46
......@@ -8,10 +8,12 @@ CONTAINS
is_inte, mt_inte, &
q, qis, qmt, qvac, qtot, qistot)
USE m_pw_tofrom_grid
USE m_mt_tofrom_grid
USE m_types
USE m_constants
!USE m_types_xcpot
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
CLASS(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym),INTENT(IN) :: sym
......@@ -21,27 +23,34 @@ CONTAINS
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_sphhar), INTENT(IN):: sphhar
TYPE(t_noco), INTENT(INOUT) :: noco
REAL, intent(in) :: is_inte(:,:), mt_inte(:,:)
TYPE(t_lapl), 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
INTEGER :: n
call init_pw_grid(xcpot, stars, sym, cell)
is_inte_mut = is_inte
!allocate potden type
call integrand%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN)
allocate(integrand%pw_w, mold=integrand%pw)
call init_mt_grid()
!put is in potden-basis
call pw_from_grid(xcpot, stars,.True., is_inte, integrand%pw, integrand%pw_w)
call pw_from_grid(xcpot, stars,.True., is_inte_mut%lapl, integrand%pw, integrand%pw_w)
!put mt in potden-basis
do n = 1,atoms%ntype
call mt_from_grid(atoms,sphhar,n,input%jspins,mt_inte,integrand%mt(:,0:,n,:))
call mt_from_grid(atoms,sphhar,n,input%jspins,mt_inte(n)%lapl,integrand%mt(:,0:,n,:))
enddo
! integrate my integrand
call integrate_cdn(stars, atoms, sym, vacuum, input, cell, oneD, integrand,&
q, qis, qmt, qvac, qtot, qistot)
call finish_pw_grid()
END SUBROUTINE integrate_grid
SUBROUTINE integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, integrand, &
......
......@@ -13,7 +13,7 @@ if (${CMAKE_Fortran_COMPILER_ID} MATCHES "Intel")
endif()
set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -xHost -O2 -g")
if (${CMAKE_Fortran_COMPILER_VERSION} VERSION_LESS "19.0.0.0")
set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -C -traceback -O0 -g -ftrapuv -check uninit -check pointers -DCPP_DEBUG")
set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -C -traceback -O0 -g -ftrapuv -check uninit -check pointers -DCPP_DEBUG -warn=all")
else()
set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -CB -traceback -O0 -g -ftrapuv -check uninit -check pointers -DCPP_DEBUG")
endif()
......@@ -44,5 +44,5 @@ elseif(${CMAKE_Fortran_COMPILER_ID} MATCHES "GNU")
endif()
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-none -fopenmp -fdefault-real-8 -Wno-missing-include-dirs")
set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -O2 -g")
set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -fdump-core -Wall -Wextra -Warray-temporaries -fbacktrace -fcheck=all -finit-real=nan -O0 -g -DCPP_DEBUG")
set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -fdump-core -Wall -Wextra -Wno-array-temporaries -fbacktrace -fcheck=all -finit-real=nan -O0 -g -DCPP_DEBUG")
endif()
......@@ -82,8 +82,8 @@ MODULE m_constants
! Hartree and Rydbergs changed by fac = 1.0 or 2.0
REAL, INTENT (IN) :: fac
c_light = 137.0359895e0 * fac
!c_light = 1e6*fac
!c_light = 137.0359895e0 * fac
c_light = 1e6*fac
END FUNCTION c_light
END MODULE m_constants
......@@ -35,7 +35,7 @@ CONTAINS
IMPLICIT NONE
TYPE(t_results), INTENT(INOUT) :: results
CLASS(t_xcpot), INTENT(IN) :: xcpot
CLASS(t_xcpot), INTENT(INOUT) :: xcpot
TYPE(t_hybrid), INTENT(IN) :: hybrid
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: dimension
......@@ -45,7 +45,7 @@ CONTAINS
TYPE(t_input), INTENT(IN) :: input
TYPE(t_field), INTENT(INOUT) :: field !efield can be modified
TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_noco), INTENT(INOUT) :: noco
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_cell), INTENT(IN) :: cell
......@@ -103,6 +103,45 @@ CONTAINS
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vx)
#endif
call integrate_lapl(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)
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
call integrate_grid(xcpot, stars, atoms, sym, vacuum, input, cell, oneD, sphhar,noco,&
xcpot%is_lapl, xcpot%mt_lapl, &
q, qis, qmt, qvac, qtot, qistot)
write (*,*) "lapl integration:"
write (*,*) "q = "
write (*,*) q
write (*,*) "qis = "
write (*,*) qis
write (*,*) "qmt = "
write (*,*) qmt
write (*,*) "qvac = "
write (*,*) qvac
write (*,*) "qtot = "
write (*,*) qtot
write (*,*) "q = "
write (*,*) qistot
END SUBROUTINE integrate_lapl
END MODULE m_vgen
......@@ -15,11 +15,17 @@ MODULE m_types_xcpot
use m_types_potden
IMPLICIT NONE
PRIVATE
PUBLIC :: t_xcpot,t_gradients
PUBLIC :: t_xcpot,t_gradients, t_lapl
type t_lapl
real, allocatable :: lapl(:,:)
end type t_lapl
TYPE,ABSTRACT :: t_xcpot
REAL :: gmaxxc
TYPE(t_potden) :: comparison_kinED_pw(3)
TYPE(t_lapl) :: is_lapl
TYPE(t_lapl), allocatable :: mt_lapl(:)
CONTAINS
PROCEDURE :: vxc_is_LDA=>xcpot_vxc_is_LDA
PROCEDURE :: exc_is_LDA=>xcpot_exc_is_LDA
......
......@@ -233,6 +233,9 @@ CONTAINS
!libxc uses the spin as a first index, hence we have to transpose....
ALLOCATE(vxc_tmp(SIZE(vxc,2),SIZE(vxc,1)));vxc_tmp=0.0
ALLOCATE(vx_tmp(SIZE(vx,2),SIZE(vx,1)));vx_tmp=0.0
write (*,*) "vxc_tmp shape = ", shape(vxc_tmp)
write (*,*) "vx_tmp shape = ", shape(vx_tmp)
IF (xcpot%needs_grad()) THEN
IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
ALLOCATE(vsigma,mold=grad%vsigma)
......
......@@ -248,7 +248,7 @@ CONTAINS
REAL,INTENT(INOUT) :: v_in(0:,:)
LOGICAL,INTENT(in) :: l_pw_w
COMPLEX,INTENT(INOUT) :: v_out_pw(:,:)
COMPLEX,INTENT(INOUT),OPTIONAL::v_out_pw_w(:,:)
COMPLEX,INTENT(INOUT),OPTIONAL:: v_out_pw_w(:,:)
INTEGER :: js,k,i
......
......@@ -36,7 +36,7 @@ CONTAINS
IMPLICIT NONE
CLASS(t_xcpot), INTENT(IN) :: xcpot
CLASS(t_xcpot), INTENT(INOUT) :: xcpot
TYPE(t_hybrid), INTENT(IN) :: hybrid
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: dimension
......
......@@ -38,7 +38,7 @@ CONTAINS
USE m_npy
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
CLASS(t_xcpot),INTENT(INOUT) :: xcpot
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
......@@ -58,6 +58,9 @@ 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
ALLOCATE(v_xc,mold=rho)
ALLOCATE(v_x,mold=rho)
......
......@@ -36,7 +36,7 @@ CONTAINS
USE m_juDFT_string
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
CLASS(t_xcpot),INTENT(INOUT) :: xcpot
TYPE(t_dimension),INTENT(IN) :: dimension
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_obsolete),INTENT(IN) :: obsolete
......@@ -110,8 +110,12 @@ CONTAINS
n_stride=1
#endif
call save_npy("rmsh.npy", atoms%rmsh)
if(.not. allocated(xcpot%mt_lapl)) allocate(xcpot%mt_lapl(atoms%ntype))
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
!
! calculate the ex.-cor. potential
CALL xcpot%get_vxc(input%jspins,ch(:nsp*atoms%jri(n),:),v_xc(:nsp*atoms%jri(n),:),v_x(:nsp*atoms%jri(n),:),grad)
......
......@@ -38,10 +38,18 @@ CONTAINS
ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma))
CALL mt_to_grid(xcpot,n_sigma,atoms,sphhar,vsigma_mt,n,grad=grad_vsigma)
fname = merge("mt=" //int2str(atom_num) // "_lapl.npy","mt_lapl.npy", present(atom_num))
if(present(atom_num)) then
fname = "mt=" //int2str(atom_num) // "_lapl.npy"
else
fname = "mt_lapl.npy"
endif
call save_npy(fname, grad%laplace)
fname = merge("mt="// int2str(atom_num) // "_grad.npy", "mt_grad.npy", present(atom_num))
if(present(atom_num)) then
fname = "mt="// int2str(atom_num) // "_grad.npy"
else
fname = "mt_grad.npy"
endif
call save_npy(fname, grad%gr)
CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc)
END SUBROUTINE libxc_postprocess_gga_mt
......
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