Commit 99cd6819 authored by Alexander Neukirchen's avatar Alexander Neukirchen

This is the big one. divergence.f90 is written and building it works. Needs to...

This is the big one. divergence.f90 is written and building it works. Needs to be tested in program flow though.
parent 7315bd26
......@@ -48,7 +48,7 @@ vgen/b_field.F90
vgen/write_xcstuff.f90
vgen/xy_av_den.f90
vgen/VYukawaFilm.f90
vgen/mt_grad_on_grid.F90
vgen/divergence.f90
)
#vdW Stuff
set(fleur_F90 ${fleur_F90}
......
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! Copyright (c) 2019 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_pw_grad_on_grid
MODULE m_divergence
USE m_types
PRIVATE
REAL,PARAMETER:: d_15=1.e-15
PUBLIC :: mt_div, pw_div, divergence
CONTAINS
SUBROUTINE pw_do_div(cell,stars,n,atoms,sphhar,sym,xcB,div)
SUBROUTINE mt_div(jspins,n,atoms,sphhar,sym,xcB,div)
!-----------------------------------------------------------------------------!
!By use of the cartesian components of a field and its cartesian derivatives !
!in the interstitial space and vacuum: !
!By use of the cartesian components of a field, its radial/angular derivati- !
!ves in the muffin tin at each spherical grid point and the corresponding an- !
!gles: !
! !
!Make the divergence of said field in real space and store it as a source !
!density, again represented by pw-coefficients in a potden. !
!density, again represented by mt-coefficients in a potden. !
! !
!Code by A. Neukirchen, September 2019 !
!-----------------------------------------------------------------------------!
USE m_grdrsis
USE m_pw_tofrom_grid
USE m_mt_tofrom_grid
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
INTEGER,INTENT(IN) :: jspins
LOGICAL,INTENT(IN) :: l_noco
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
INTEGER, INTENT(IN) :: jspins, n
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_potden), dimension(3), INTENT(INOUT) :: xcB
TYPE(t_potden), INTENT(INOUT) :: div
TYPE(t_gradients) :: gradx, grady, gradz
REAL, ALLOCATABLE :: div_temp(:, :)
REAL, ALLOCATABLE :: thet(:), phi(:)
REAL :: r,th,ph
INTEGER :: jr, k, nsp, kt
ALLOCATE (gradx%gr(3,stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft,jspins),grady%gr(3,atoms%jri(n)*nsp,jspins),gradz%gr(3,atoms%jri(n)*nsp,jspins))
nsp = atoms%nsp()
ALLOCATE (gradx%gr(3,atoms%jri(n)*nsp,jspins),grady%gr(3,atoms%jri(n)*nsp,jspins),gradz%gr(3,atoms%jri(n)*nsp,jspins))
ALLOCATE (div_temp(atoms%jri(n)*nsp,jspins))
CALL grdrsis(,cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,gradx)
CALL mt_do_grad(xcpot, jspins, n, atoms, sphhar, sym, xcB(2)%mt(:,0:,n,:), grady)
CALL mt_do_grad(xcpot, jspins, n, atoms, sphhar, sym, xcB(3)%mt(:,0:,n,:), gradz)
CALL init_mt_grid(jspins, atoms, sphhar, .TRUE., sym)
CALL mt_to_grid(.TRUE., jspins, atoms, sphhar, xcB(1)%mt(:,0:,n,:), n, gradx)
CALL mt_to_grid(.TRUE., jspins, atoms, sphhar, xcB(2)%mt(:,0:,n,:), n, grady)
CALL mt_to_grid(.TRUE., jspins, atoms, sphhar, xcB(3)%mt(:,0:,n,:), n, gradz)
kt = 0
DO jr = 1, atoms%jri(n)
......@@ -57,11 +62,82 @@ CONTAINS
kt = kt+nsp
ENDDO ! jr
CALL pw_from_grid(xcpot,stars,l_pw_w,v_in,v_out_pw,v_out_pw_w)
CALL mt_from_grid(atoms, sphhar, n, jspins, div_temp, div%mt(:,0:,n,:))
DEALLOCATE (ylh, wt, rx, thet, phi, ylht, ylhtt, ylhf, ylhff, ylhtf)
CALL finish_mt_grid
END SUBROUTINE mt_div
SUBROUTINE pw_div(ifftxc3,jspins,stars,cell,noco,sym,xcB,div)
!-----------------------------------------------------------------------------!
!By use of the cartesian components of a field, its radial/angular derivati- !
!ves in the muffin tin at each spherical grid point and the corresponding an- !
!gles: !
! !
!Make the divergence of said field in real space and store it as a source !
!density, again represented by mt-coefficients in a potden. !
! !
!Code by A. Neukirchen, September 2019 !
!-----------------------------------------------------------------------------!
USE m_pw_tofrom_grid
IMPLICIT NONE
END SUBROUTINE pw_do_div
END MODULE m_pw_grad_on_grid
INTEGER, INTENT(IN) :: jspins, ifftxc3
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_potden), dimension(3), INTENT(INOUT) :: xcB
TYPE(t_potden), INTENT(INOUT) :: div
TYPE(t_gradients) :: gradx, grady, gradz
REAL, ALLOCATABLE :: div_temp(:, :)
INTEGER :: i, nsp
nsp = 3*ifftxc3
ALLOCATE (gradx%gr(3,nsp,jspins),grady%gr(3,nsp,jspins),gradz%gr(3,nsp,jspins))
ALLOCATE (div_temp(nsp,jspins))
CALL init_pw_grid(.TRUE.,stars,sym,cell)
CALL pw_to_grid(.TRUE.,jspins,noco%l_noco,stars,cell,xcB(1)%pw,gradx)
CALL pw_to_grid(.TRUE.,jspins,noco%l_noco,stars,cell,xcB(2)%pw,grady)
CALL pw_to_grid(.TRUE.,jspins,noco%l_noco,stars,cell,xcB(3)%pw,gradz)
DO i = 1, nsp
div_temp(i,1)=gradx%gr(1,i,1)+grady%gr(2,i,1)+gradz%gr(3,i,1)
ENDDO ! i
CALL pw_from_grid(.TRUE.,stars,.TRUE.,div_temp,div%pw,div%pw_w)
CALL finish_pw_grid()
END SUBROUTINE pw_div
SUBROUTINE divergence(jspins,n,ifftxc3,atoms,sphhar,sym,stars,cell,vacuum,noco,xcB,div)
USE m_types
IMPLICIT NONE
INTEGER, INTENT(IN) :: jspins, n, ifftxc3
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_potden), dimension(3), INTENT(INOUT) :: xcB
TYPE(t_potden), INTENT(OUT) :: div
CALL div%init(stars,atoms,sphhar,vacuum,noco,jspins,1001)
CALL mt_div(jspins,n,atoms,sphhar,sym,xcB,div)
CALL pw_div(ifftxc3,jspins,stars,cell,noco,sym,xcB,div)
END SUBROUTINE divergence
END MODULE m_divergence
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_mt_grad_on_grid
USE m_types
PRIVATE
REAL,PARAMETER:: d_15=1.e-15
INTEGER, PARAMETER :: ndvgrd = 6
REAL, ALLOCATABLE :: ylh(:, :, :), ylht(:, :, :), ylhtt(:, :, :)
REAL, ALLOCATABLE :: ylhf(:, :, :), ylhff(:, :, :), ylhtf(:, :, :)
REAL, ALLOCATABLE :: wt(:), rx(:, :), thet(:), phi(:)
PUBLIC :: mt_do_grad, mt_do_div
CONTAINS
!SUBROUTINE init_mt_grad_grid()
!END SUBROUTINE init_mt_grad_grid
SUBROUTINE mt_do_grad(xcpot, jspins, n, atoms, sphhar, sym, den_mt, grad, ch)
!-----------------------------------------------------------------------------!
!By use of the radial functions/spherical harmonics representation of a field !
!component in the Muffin Tins: !
! !
!Make the gradient of said field in real space and store it. !
! !
!Based on mt_tofrom_grid.F90, A. Neukirchen, September 2019. !
!-----------------------------------------------------------------------------!
USE m_gaussp
USE m_lhglptg
USE m_grdchlh
USE m_mkgylm
IMPLICIT NONE
CLASS(t_xcpot), INTENT(IN) :: xcpot
INTEGER, INTENT(IN) :: jspins, n
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_sym), INTENT(IN) :: sym
REAL, INTENT(IN) :: den_mt(:, 0:, :)
TYPE(t_gradients), INTENT(INOUT):: grad
REAL, INTENT(OUT), OPTIONAL :: ch(:, :)
REAL, ALLOCATABLE :: chlh(:, :, :), chlhdr(:, :, :), chlhdrr(:, :, :)
REAL, ALLOCATABLE :: ch_tmp(:, :)
REAL, ALLOCATABLE :: chdr(:, :), chdt(:, :), chdf(:, :)
REAL, ALLOCATABLE :: chdrr(:, :), chdtt(:, :), chdff(:, :)
REAL, ALLOCATABLE :: chdrt(:, :), chdrf(:, :), chdtf(:, :)
INTEGER:: nd, lh, js, jr, kt, k, nsp
!Generate nspd points on a sherical shell with radius 1.0:
!The angular mesh is equidistant in phi,
!theta are the zeros of the legendre polynomials
ALLOCATE (wt(atoms%nsp()), rx(3, atoms%nsp()), thet(atoms%nsp()), phi(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))
!Generate the derivatives of the lattice harmonics:
ALLOCATE (ylht, MOLD=ylh)
ALLOCATE (ylhtt, MOLD=ylh)
ALLOCATE (ylhf, MOLD=ylh)
ALLOCATE (ylhff, MOLD=ylh)
ALLOCATE (ylhtf, MOLD=ylh)
!Calculate the required lattice harmonics and their derivatives:
CALL lhglptg(sphhar, atoms, rx, atoms%nsp(), xcpot%needs_grad(), sym, &
ylh, thet, phi, ylht, ylhtt, ylhf, ylhff, ylhtf)
nd = atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1)
nsp = atoms%nsp()
ALLOCATE (ch_tmp(nsp, jspins))
ALLOCATE (chdr(nsp, jspins), chdt(nsp, jspins), chdf(nsp, jspins), &
chdrr(nsp, jspins),chdtt(nsp, jspins), chdff(nsp, jspins), &
chdtf(nsp, jspins), chdrt(nsp, jspins), chdrf(nsp, jspins))
ALLOCATE (chlh(atoms%jmtd, 0:sphhar%nlhd, jspins))
ALLOCATE (chlhdr(atoms%jmtd, 0:sphhar%nlhd, jspins))
ALLOCATE (chlhdrr(atoms%jmtd, 0:sphhar%nlhd, jspins))
!Calculates the derivatives of the charge density in spherical
!coordinates as follows [mt are the muffin tin coefficients]:
!
!chlh=pw/r**2.
!chlhdr=d(chlh)/dr, chlhdrr=dd(chlh)/drr [a]
!rho=sum(chlh*ylh)
! --> d(rho)/dr=sum(chlhdr*ylh) for radial derivative
! --> d(rho)/d[angle]=sum(chlh*ylh[angle]) for the angular ones
! --> higher and mixed derivatives follow as one would expect
!First get chlh, chlhdr and chlhdrr for all grid points [a]:
DO lh = 0, sphhar%nlh(nd)
DO js = 1, jspins
DO jr = 1, atoms%jri(n)
chlh(jr, lh, js) = den_mt(jr, lh, js)/(atoms%rmsh(jr, n)*atoms%rmsh(jr, n))
ENDDO ! jr
CALL grdchlh(1, 1, atoms%jri(n), atoms%dx(n), atoms%rmsh(1, n), chlh(1, lh, js), &
ndvgrd, chlhdr(1, lh, js), chlhdrr(1, lh, js))
ENDDO ! js
ENDDO ! lh
kt = 0
DO jr = 1, atoms%jri(n)
!On extended grid for all jr:
!Generate rho and all its derivatives on the jr-th sphere:
ch_tmp(:, :) = 0.0
chdr(:, :) = 0.0 ! d(ch)/dr
chdt(:, :) = 0.0 ! d(ch)/dtheta
chdf(:, :) = 0.0 ! d(ch)/dfai
chdrr(:, :) = 0.0 ! dd(ch)/drr
chdtt(:, :) = 0.0 ! dd(ch)/dtt
chdff(:, :) = 0.0 ! dd(ch)/dff
chdtf(:, :) = 0.0 ! dd(ch)/dtf
chdrt(:, :) = 0.0 ! d(d(ch)/dr)dt
chdrf(:, :) = 0.0 ! d(d(ch)/dr)df
!Calculate [b] and sum [c] the density/derivatives on the angular mesh:
DO js = 1, jspins
DO lh = 0, sphhar%nlh(nd)
DO k = 1, nsp
ch_tmp(k, js) = ch_tmp(k, js) + ylh(k, lh, nd)*chlh(jr, lh, js)
chdr(k, js) = chdr(k, js) + ylh(k, lh, nd)*chlhdr(jr, lh, js)
chdrr(k, js) = chdrr(k, js) + ylh(k, lh, nd)*chlhdrr(jr, lh, js)
chdrt(k, js) = chdrt(k, js) + ylht(k, lh, nd)*chlhdr(jr, lh, js)
chdrf(k, js) = chdrf(k, js) + ylhf(k, lh, nd)*chlhdr(jr, lh, js)
chdt(k, js) = chdt(k, js) + ylht(k, lh, nd)*chlh(jr, lh, js)
chdf(k, js) = chdf(k, js) + ylhf(k, lh, nd)*chlh(jr, lh, js)
chdtt(k, js) = chdtt(k, js) + ylhtt(k, lh, nd)*chlh(jr, lh, js)
chdff(k, js) = chdff(k, js) + ylhff(k, lh, nd)*chlh(jr, lh, js)
chdtf(k, js) = chdtf(k, js) + ylhtf(k, lh, nd)*chlh(jr, lh, js)
ENDDO ! k
ENDDO ! lh
ENDDO ! js
!Combine the derivatives and put them into grad; grad%gr will contain the
!first derivatives in r, theta, phi on an angular mesh:
CALL mkgylm(jspins, atoms%rmsh(jr, n), thet, nsp, ch_tmp, chdr, chdt, &
chdf, chdrr, chdtt, chdff, chdtf, chdrt, chdrf, grad, kt)
!If requested also pipe out the real space density:
IF (PRESENT(ch)) THEN
WHERE (ABS(ch_tmp) < d_15) ch_tmp = d_15
ch(kt + 1:kt + nsp, :) = ch_tmp(:nsp, :)
ENDIF
kt = kt + nsp
END DO ! jr
END SUBROUTINE mt_do_grad
SUBROUTINE mt_do_div(xcpot,jspins,n,atoms,sphhar,sym,xcB,div)
!-----------------------------------------------------------------------------!
!By use of the cartesian components of a field, its radial/angular derivati- !
!ves in the muffin tin at each spherical grid point and the corresponding an- !
!gles: !
! !
!Make the divergence of said field in real space and store it as a source !
!density, again represented by mt-coefficients in a potden. !
! !
!Code by A. Neukirchen, September 2019 !
!-----------------------------------------------------------------------------!
USE m_mt_tofrom_grid
IMPLICIT NONE
CLASS(t_xcpot), INTENT(IN) :: xcpot
INTEGER, INTENT(IN) :: jspins, n
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_potden), dimension(3), INTENT(INOUT) :: xcB
TYPE(t_potden), INTENT(INOUT) :: div
TYPE(t_gradients) :: gradx, grady, gradz
REAL, ALLOCATABLE :: div_temp(:, :)
REAL :: r,th,ph
INTEGER :: jr, k, nsp, kt
ALLOCATE (gradx%gr(3,atoms%jri(n)*nsp,jspins),grady%gr(3,atoms%jri(n)*nsp,jspins),gradz%gr(3,atoms%jri(n)*nsp,jspins))
ALLOCATE (div_temp(atoms%jri(n)*nsp,jspins))
CALL mt_do_grad(xcpot, jspins, n, atoms, sphhar, sym, xcB(1)%mt(:,0:,n,:), gradx)
CALL mt_do_grad(xcpot, jspins, n, atoms, sphhar, sym, xcB(2)%mt(:,0:,n,:), grady)
CALL mt_do_grad(xcpot, jspins, n, atoms, sphhar, sym, xcB(3)%mt(:,0:,n,:), gradz)
kt = 0
DO jr = 1, atoms%jri(n)
r=atoms%rmsh(jr, n)
DO k = 1, nsp
th = thet(k)
ph = phi(k)
div_temp(kt+nsp,1) = (SIN(th)*COS(ph)*gradx%gr(1,kt+nsp,jspins) + SIN(th)*SIN(ph)*grady%gr(1,kt+nsp,jspins) + COS(th)*gradz%gr(1,kt+nsp,jspins))&
+(COS(th)*COS(ph)*gradx%gr(2,kt+nsp,jspins) + COS(th)*SIN(ph)*grady%gr(2,kt+nsp,jspins) - SIN(th)*gradz%gr(2,kt+nsp,jspins))/r&
-(SIN(ph)*gradx%gr(3,kt+nsp,jspins) + COS(ph)*grady%gr(3,kt+nsp,jspins))/(r*SIN(th))
ENDDO ! k
kt = kt+nsp
ENDDO ! jr
CALL mt_from_grid(atoms, sphhar, n, jspins, div_temp, div%mt(:,0:,n,:))
DEALLOCATE (ylh, wt, rx, thet, phi, ylht, ylhtt, ylhf, ylhff, ylhtf)
END SUBROUTINE mt_do_div
END MODULE m_mt_grad_on_grid
......@@ -217,23 +217,10 @@ CONTAINS
ENDDO !jdm
ENDDO !idm
END IF
IF (PRESENT(xcpot)) THEN
CALL xcpot%alloc_gradients(ifftxc3,jspins,grad)
END IF
!!!!!!!THIS IS A QUICKFIX! TO BE REMOVED ASAP!!!!!!!!!
!A. Neukirchen 25.09.19
!IF (ALLOCATED(grad%agrt)) THEN
!DEALLOCATE(grad%agrt,grad%agru,grad%agrd)
!DEALLOCATE(grad%g2ru,grad%g2rd,grad%gggrt)
!DEALLOCATE(grad%gggru,grad%gzgr,grad%g2rt)
!DEALLOCATE(grad%gggrd,grad%grgru,grad%grgrd)
!ENDIF
!ALLOCATE(grad%agrt(ifftxc3),grad%agru(ifftxc3),grad%agrd(ifftxc3))
!ALLOCATE(grad%g2ru(ifftxc3),grad%g2rd(ifftxc3),grad%gggrt(ifftxc3))
!ALLOCATE(grad%gggru(ifftxc3),grad%gzgr(ifftxc3),grad%g2rt(ifftxc3))
!ALLOCATE(grad%gggrd(ifftxc3),grad%grgru(ifftxc3),grad%grgrd(ifftxc3))
!!!!!!!!!
!
! calculate the quantities such as abs(grad(rho)),.. used in
......
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