Commit 4e2627f6 authored by Alexander Neukirchen's avatar Alexander Neukirchen

some small on getting div to run; extended it to the next steps; debugging tomorrow

parent 3d474af6
......@@ -38,6 +38,8 @@ MODULE m_constants
INTEGER, PARAMETER :: PLOT_POT_COU=9
INTEGER, PARAMETER :: PLOT_POT_VXC=10
INTEGER, PARAMETER :: PLOT_SPECIAL=14
INTEGER, PARAMETER :: PLOT_SPECIAL2=15
INTEGER, PARAMETER :: PLOT_SPECIAL3=16
......
......@@ -65,6 +65,7 @@ CONTAINS
USE m_ylm
USE m_metagga
USE m_divergence
! USE m_gradfromgrid
USE m_plot
#ifdef CPP_MPI
USE m_mpi_bc_potden
......@@ -98,9 +99,9 @@ CONTAINS
TYPE(t_mpi) :: mpi
TYPE(t_coreSpecInput) :: coreSpecInput
TYPE(t_wann) :: wann
TYPE(t_potden) :: vTot, vx, vCoul, vTemp, vxcForPlotting
TYPE(t_potden) :: vTot, vx, vCoul, vTemp, vxcForPlotting, vDiv
TYPE(t_potden) :: inDen, outDen, EnergyDen, divB
TYPE(t_potden), dimension(3):: xcB
TYPE(t_potden), dimension(3):: xcB, graddiv
CLASS(t_xcpot), ALLOCATABLE :: xcpot
CLASS(t_forcetheo), ALLOCATABLE :: forcetheo
......@@ -500,16 +501,38 @@ CONTAINS
! DIVERGENCE
CALL divB%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
ALLOCATE(divB%pw_w,mold=divB%pw)
! CALL divB%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
! ALLOCATE(divB%pw_w,mold=divB%pw)
DO i=1,atoms%ntype
CALL divergence(input%jspins,i,stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft,atoms,sphhar,sym,stars,cell,vacuum,noco,xcB,divB)
END DO
IF ((sliceplot%iplot.NE.0).AND.(mpi%irank==0) ) THEN
CALL makeplots(mpi,sym,stars,vacuum,atoms,sphhar,input,cell,oneD,noco,sliceplot,divB,PLOT_SPECIAL)
END IF
! DO i=1,atoms%ntype
! CALL divergence(input%jspins,i,stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft,atoms,sphhar,sym,stars,cell,vacuum,noco,xcB,divB)
! END DO
! CALL vDiv%init(stars,atoms,sphhar,vacuum,noco,1,POTDEN_TYPE_POTCOUL)
! ALLOCATE(vDiv%pw_w(SIZE(vDiv%pw,1),size(vDiv%pw,2)))
! vDiv%pw_w = CMPLX(0.0,0.0)
! CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,divB,vDiv)
! DO i=1,3
! CALL graddiv(i)%init(stars,atoms,sphhar,vacuum,noco,1,POTDEN_TYPE_DEN)
! ALLOCATE(graddiv(i)%pw_w,mold=graddiv(i)%pw)
! ENDDO
! DO i=1,atoms%ntype
! CALL gradfromgrid(input%jspins,i,stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft,atoms,sphhar,sym,stars,cell,vacuum,noco,vDiv,graddiv)
! END DO
! DO i=1,atoms%ntype
! xcB(i)=xcB(i)+graddiv/(4.0*pi_const)
! END DO
! IF ((sliceplot%iplot.NE.0).AND.(mpi%irank==0) ) THEN
! CALL makeplots(mpi,sym,stars,vacuum,atoms,sphhar,input,cell,oneD,noco,sliceplot,divB,PLOT_SPECIAL)
! DO i=1,3
! CALL makeplots(mpi,sym,stars,vacuum,atoms,sphhar,input,cell,oneD,noco,sliceplot,graddiv(i),PLOT_SPECIAL2)
! CALL makeplots(mpi,sym,stars,vacuum,atoms,sphhar,input,cell,oneD,noco,sliceplot,xcB(i),PLOT_SPECIAL3)
! END IF
CALL add_usage_data("Iterations",iter)
......
......@@ -879,7 +879,7 @@ CONTAINS
END IF
!Plotting the divergence of B_xc_vec. identifier: 14
!Plotting the divergence of B_xc_vec. identifier: 14 (Pot: 15, B_corrected: 16)
!No core subtraction done!
! --> Additive term for iplot: 16384
IF (plot_const.EQ.14) THEN
......@@ -887,6 +887,23 @@ CONTAINS
score = .FALSE.
potnorm = .FALSE.
CALL savxsf(potnorm,oneD,stars,vacuum,sphhar,atoms,input,sym,cell,sliceplot,noco,score,denName,denmat)
END IF
IF (plot_const.EQ.15) THEN
denName = 'divPot'
score = .FALSE.
potnorm = .FALSE.
CALL savxsf(potnorm,oneD,stars,vacuum,sphhar,atoms,input,sym,cell,sliceplot,noco,score,denName,denmat)
END IF
IF (plot_const.EQ.16) THEN
denName = 'xcBcorr'
score = .FALSE.
potnorm = .FALSE.
CALL savxsf(potnorm,oneD,stars,vacuum,sphhar,atoms,input,sym,cell,sliceplot,noco,score,denName,denmat)
END IF
......
......@@ -98,17 +98,37 @@ CONTAINS
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 xcpot%alloc_gradients(ifftxc3,1,gradx)
! CALL xcpot%alloc_gradients(ifftxc3,2,gradx)
! CALL xcpot%alloc_gradients(ifftxc3,3,gradx)
ALLOCATE (gradx%gr(3,nsp,1))
ALLOCATE (grady%gr(3,nsp,1))
ALLOCATE (gradz%gr(3,nsp,1))
! ALLOCATE (gradx%agrt(nsp),grady%agrt(nsp),gradz%agrt(nsp))
! ALLOCATE (gradx%agru(nsp),grady%agru(nsp),gradz%agru(nsp))
! ALLOCATE (gradx%agrd(nsp),grady%agrd(nsp),gradz%agrd(nsp))
! ALLOCATE (gradx%gggrt(nsp),grady%gggrt(nsp),gradz%gggrt(nsp))
! ALLOCATE (gradx%gggru(nsp),grady%gggru(nsp),gradz%gggru(nsp))
! ALLOCATE (gradx%gggrd(nsp),grady%gggrd(nsp),gradz%gggrd(nsp))
! ALLOCATE (gradx%gzgr(nsp),grady%gzgr(nsp),gradz%gzgr(nsp))
! ALLOCATE (gradx%g2rt(nsp),grady%g2rt(nsp),gradz%g2rt(nsp))
! ALLOCATE (gradx%g2ru(nsp),grady%g2ru(nsp),gradz%g2ru(nsp))
! ALLOCATE (gradx%g2rd(nsp),grady%g2rd(nsp),gradz%g2rd(nsp))
ALLOCATE (div_temp(ifftxc3,1))
CALL init_pw_grid(.TRUE.,stars,sym,cell)
CALL pw_to_grid(.TRUE.,1,noco%l_noco,stars,cell,xcB(1)%pw,gradx)
CALL pw_to_grid(.TRUE.,1,noco%l_noco,stars,cell,xcB(2)%pw,grady)
CALL pw_to_grid(.TRUE.,1,noco%l_noco,stars,cell,xcB(3)%pw,gradz)
CALL pw_to_grid(.TRUE.,1,.FALSE.,stars,cell,xcB(1)%pw,gradx)
CALL pw_to_grid(.TRUE.,1,.FALSE.,stars,cell,xcB(2)%pw,grady)
CALL pw_to_grid(.TRUE.,1,.FALSE.,stars,cell,xcB(3)%pw,gradz)
DO i = 1, nsp
! DEALLOCATE(gradx%agrt,gradx%agru,gradx%agrd,gradx%gggrt,gradx%gggru,gradx%gggrd,gradx%gzgr,gradx%g2rt,gradx%g2ru,gradx%g2rd)
! DEALLOCATE(grady%agrt,grady%agru,grady%agrd,grady%gggrt,grady%gggru,grady%gggrd,grady%gzgr,grady%g2rt,grady%g2ru,grady%g2rd)
! DEALLOCATE(gradz%agrt,gradz%agru,gradz%agrd,gradz%gggrt,gradz%gggru,gradz%gggrd,gradz%gzgr,gradz%g2rt,gradz%g2ru,gradz%g2rd)
print *, ifftxc3
DO i = 1, ifftxc3
div_temp(i,1)=gradx%gr(1,i,1)+grady%gr(2,i,1)+gradz%gr(3,i,1)
ENDDO ! i
......
......@@ -67,18 +67,20 @@ CONTAINS
IF (ANY(SHAPE(vl).NE.SHAPE(dvx))) CALL judft_error("Gradients for internal GGA called with inconsistent sizes",hint="This is a bug")
DO i = 1,size(grad%agrt)
grad%agrt(i) = 0.0
grad%agru(i) = 0.0
grad%agrd(i) = 0.0
grad%gggrt(i) = 0.0
grad%gggru(i) = 0.0
grad%gggrd(i) = 0.0
grad%gzgr(i) = 0.0
grad%g2rt(i) = 0.0
grad%g2ru(i) = 0.0
grad%g2rd(i) = 0.0
ENDDO
IF(ALLOCATED(grad%agrt)) THEN
DO i = 1,size(grad%agrt)
grad%agrt(i) = 0.0
grad%agru(i) = 0.0
grad%agrd(i) = 0.0
grad%gggrt(i) = 0.0
grad%gggru(i) = 0.0
grad%gggrd(i) = 0.0
grad%gzgr(i) = 0.0
grad%g2rt(i) = 0.0
grad%g2ru(i) = 0.0
grad%g2rd(i) = 0.0
ENDDO
END IF
IF (jspins.eq.1) THEN
......@@ -118,7 +120,7 @@ CONTAINS
dvyzt = dvyzu + dvyzd
dvzxt = dvzxu + dvzxd
dvxyt = dvxyu + dvxyd
IF(ALLOCATED(grad%agrt)) THEN
! agr: abs(grad(ro)), t,u,d for total, up and down.
grad%agrt(i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
......@@ -157,7 +159,7 @@ CONTAINS
grad%g2ru(i) = dvxxu + dvyyu + dvzzu
grad%g2rd(i) = dvxxd + dvyyd + dvzzd
END IF
ENDDO
ELSE
......@@ -197,7 +199,7 @@ CONTAINS
dvyzt = dvyzu + dvyzd
dvzxt = dvzxu + dvzxd
dvxyt = dvxyu + dvxyd
IF(ALLOCATED(grad%agrt)) THEN
!c agr: abs(grad(ro)), t,u,d for total, up and down.
grad%agrt(i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
......@@ -235,7 +237,7 @@ CONTAINS
grad%g2rt(i) = dvxxt + dvyyt + dvzzt
grad%g2ru(i) = dvxxu + dvyyu + dvzzu
grad%g2rd(i) = dvxxd + dvyyd + dvzzd
END IF
ENDDO
ENDIF
......
......@@ -231,7 +231,9 @@ CONTAINS
rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
ELSE
!Dummy rho (only possible if grad is used for libxc mode)
CALL mkgxyz3 (RESHAPE((/0.0/),(/1,1/)),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
!CALL mkgxyz3 (RESHAPE((/0.0/),(/1,1/)),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
! rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
CALL mkgxyz3 (0*rhd1(0:,:,1),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
END IF
......
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