Commit 5c6ae6b5 authored by Alexander Neukirchen's avatar Alexander Neukirchen

Playing around with boundaries and cutting components.

parent 9b23d5db
......@@ -517,20 +517,20 @@ CONTAINS
! Test: Build a field, for which the theoretical divergence etc. are known and
! compare with the result of the routine.
CALL builddivtest(stars,atoms,sphhar,vacuum,sym,cell,1,testDen)
!CALL makeBxc(stars,atoms,sphhar,vacuum,input,noco,vTot,testDen)
!CALL matrixsplit(stars, atoms, sphhar, vacuum, input, noco, 1.0, inDen, dummyDen, testDen(1), testDen(2), testDen(3))
!CALL builddivtest(stars,atoms,sphhar,vacuum,sym,cell,1,testDen)
!CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,inDen,1.0,testDen)
CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,vtot,2.0,testDen)
!CALL checkplotinp()
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDen ', dummyDen, testDen(1), testDen(2), testDen(3))
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDen ', testDen(1), testDen(1), testDen(2), testDen(3))
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDeny ', testDen(2))
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDenz ', testDen(3))
!CALL sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,testDen)
DO i=1,3
CALL testGrad(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(testGrad(i)%pw_w,mold=testGrad(i)%pw)
ENDDO
CALL divpotgrad(stars,atoms,sphhar,vacuum,sym,cell,noco,testDen(2),testGrad)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testGrad ', testGrad(1), testGrad(1), testGrad(2), testGrad(3))
CALL sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,testDen)
!DO i=1,3
!CALL testGrad(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
!ALLOCATE(testGrad(i)%pw_w,mold=testGrad(i)%pw)
!ENDDO
!CALL divpotgrad(stars,atoms,sphhar,vacuum,sym,cell,noco,testDen(3),testGrad)
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testGrad ', testGrad(1), testGrad(1), testGrad(2), testGrad(3))
CALL add_usage_data("Iterations",iter)
......
......@@ -80,7 +80,7 @@ CONTAINS
!sum up both spins in den into workden
CALL den%sum_both_spin(workden)
CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,workden,vCoul,results)
CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,.FALSE.,workden,vCoul,results)
CALL vCoul%copy_both_spin(vTot)
vCoul%mt(:,:,:,input%jspins)=vCoul%mt(:,:,:,1)
......
......@@ -58,7 +58,7 @@ CONTAINS
#endif
IF ( .NOT. input%film ) THEN
CALL vgen_coulomb( 1, mpi, DIMENSION, oneD, input, field, vacuum, sym, stars, cell, &
sphhar, atoms, resDen, vYukawa )
sphhar, atoms, .FALSE., resDen, vYukawa )
ELSE
if( mpi%irank == 0 ) then
call resDenMod%init( stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN )
......
......@@ -1096,10 +1096,10 @@ CONTAINS
.FALSE., .TRUE., 'phiDiv ', phi)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .TRUE., 'gradPhiDiv ', divGrx, divGrx, divGry, divGrz)
.FALSE., .FALSE., 'gradPhiDiv ', divGrx, divGrx, divGry, divGrz)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .TRUE., 'bCorrected ', xcBmodx, xcBmodx, xcBmody, xcBmodz)
.FALSE., .FALSE., 'bCorrected ', xcBmodx, xcBmodx, xcBmody, xcBmodz)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .FALSE., 'divCorrected ', div2)
......
......@@ -401,7 +401,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
CALL overallVCoul%resetPotDen()
ALLOCATE(overallVCoul%pw_w(size(overallVCoul%pw,1),size(overallVCoul%pw,2)))
overallVCoul%pw_w(:,:) = 0.0
CALL vgen_coulomb(1,mpi,DIMENSION,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,overallDen,overallVCoul)
CALL vgen_coulomb(1,mpi,DIMENSION,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,.FALSE.,overallDen,overallVCoul)
CALL convol(stars,overallVCoul%pw_w(:,1),overallVCoul%pw(:,1),stars%ufft) ! Is there a problem with a second spin?!
#ifdef CPP_MPI
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallVCoul)
......
......@@ -9,7 +9,7 @@ MODULE m_divergence
PUBLIC :: mt_div, pw_div, divergence, mt_grad, pw_grad, divpotgrad
CONTAINS
SUBROUTINE mt_div(n,atoms,sphhar,sym,bxc,div)
SUBROUTINE mt_div(atoms,sphhar,sym,bxc,div)
USE m_mt_tofrom_grid
!--------------------------------------------------------------------------
......@@ -25,68 +25,72 @@ CONTAINS
IMPLICIT NONE
INTEGER, INTENT(IN) :: 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(IN) :: bxc
TYPE(t_potden), DIMENSION(3), INTENT(INOUT) :: bxc
TYPE(t_potden), INTENT(INOUT) :: div
TYPE(t_gradients) :: gradx, grady, gradz
TYPE(t_gradients) :: gradx, grady, gradz
REAL, ALLOCATABLE :: thet(:), phi(:), div_temp(:, :), div_temp2(:, :)
REAL, ALLOCATABLE :: thet(:), phi(:), div_temp(:, :)
REAL :: r, th, ph, eps
INTEGER :: jr, k, nsp, kt, i, lh, lhmax
INTEGER :: jr, k, nsp, kt, i, lh, lhmax, n
nsp = atoms%nsp()
lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
eps=1.e-10
eps=1.e-12
ALLOCATE (gradx%gr(3,atoms%jri(n)*nsp,1),grady%gr(3,atoms%jri(n)*nsp,1), &
gradz%gr(3,atoms%jri(n)*nsp,1))
ALLOCATE (div_temp(atoms%jri(n)*nsp,1))
ALLOCATE (div_temp2(atoms%jri(n)*nsp,1))
ALLOCATE (thet(atoms%nsp()),phi(atoms%nsp()))
CALL init_mt_grid(1, atoms, sphhar, .TRUE., sym, thet, phi)
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, bxc(1)%mt(:,0:,n,:), n, gradx)
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, bxc(2)%mt(:,0:,n,:), n, grady)
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, bxc(3)%mt(:,0:,n,:), 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+k,1) = (SIN(th)*COS(ph)*gradx%gr(1,kt+k,1) + SIN(th)*SIN(ph)*grady%gr(1,kt+k,1) + COS(th)*gradz%gr(1,kt+k,1))&
+(COS(th)*COS(ph)*gradx%gr(2,kt+k,1) + COS(th)*SIN(ph)*grady%gr(2,kt+k,1) - SIN(th)*gradz%gr(2,kt+k,1))/r&
-(SIN(ph)*gradx%gr(3,kt+k,1) - COS(ph)*grady%gr(3,kt+k,1))/(r*SIN(th))
ENDDO ! k
kt = kt+nsp
ENDDO ! jr
CALL mt_from_grid(atoms, sphhar, n, 1, div_temp, div%mt(:,0:,n,:))
DO jr = 1, atoms%jri(n)
DO lh=0, lhmax
IF (ABS(div%mt(jr,lh,n,1))<eps) THEN
div%mt(jr,lh,n,:)=0.0
END IF
END DO
DO i=1,3
bxc(i)%mt(:,1:,:,:)=0.0
END DO
kt = 0
DO jr = 1, atoms%jri(n)
r =atoms%rmsh(jr, n)
div%mt(jr,0:,n,:) = div%mt(jr,0:,n,:)*r**2
kt = kt+nsp
ENDDO ! jr
CALL mt_to_grid(.FALSE., 1, atoms, sphhar, div%mt(:,0:,n,:), n, gradx, div_temp2)
DO n = 1, atoms%ntype
ALLOCATE (gradx%gr(3,atoms%jri(n)*nsp,1),grady%gr(3,atoms%jri(n)*nsp,1), &
gradz%gr(3,atoms%jri(n)*nsp,1))
ALLOCATE (div_temp(atoms%jri(n)*nsp,1))
div_temp=0.0
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, bxc(1)%mt(:,0:,n,:), n, gradx)
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, bxc(2)%mt(:,0:,n,:), n, grady)
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, bxc(3)%mt(:,0:,n,:), 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+k,1) = (SIN(th)*COS(ph)*gradx%gr(1,kt+k,1) + SIN(th)*SIN(ph)*grady%gr(1,kt+k,1) + COS(th)*gradz%gr(1,kt+k,1))&
+(COS(th)*COS(ph)*gradx%gr(2,kt+k,1) + COS(th)*SIN(ph)*grady%gr(2,kt+k,1) - SIN(th)*gradz%gr(2,kt+k,1))/r&
-(SIN(ph)*gradx%gr(3,kt+k,1) - COS(ph)*grady%gr(3,kt+k,1))/(r*SIN(th))
END DO !k
kt = kt+nsp
END DO !jr
CALL mt_from_grid(atoms, sphhar, n, 1, div_temp, div%mt(:,0:,n,:))
DEALLOCATE(gradx%gr,grady%gr,gradz%gr,div_temp)
END DO !n
!DO n = 1, atoms%ntype
! lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
! DO jr = 1, atoms%jri(n)
! DO lh=0, lhmax
! IF (ABS(div%mt(jr,lh,n,1))<eps) THEN
! IF (lh/=1) THEN
! div%mt(jr,lh,n,:)=0.0
! END IF
! END DO
! END DO
!END DO
DO n = 1, atoms%ntype
lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
DO lh = 0, lhmax
div%mt(:,lh,n,1) = div%mt(:,lh,n,1)*atoms%rmsh(:, n)**2
ENDDO !lh
ENDDO !n
CALL finish_mt_grid
......@@ -155,15 +159,10 @@ CONTAINS
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), DIMENSION(3), INTENT(IN) :: bxc
TYPE(t_potden), DIMENSION(3), INTENT(INOUT) :: bxc
TYPE(t_potden), INTENT(INOUT) :: div
INTEGER :: n
DO n=1,atoms%ntype
CALL mt_div(n,atoms,sphhar,sym,bxc,div)
END DO
CALL mt_div(atoms,sphhar,sym,bxc,div)
CALL pw_div(stars,sym,cell,noco,bxc,div)
END SUBROUTINE divergence
......@@ -191,6 +190,7 @@ CONTAINS
TYPE(t_potden), INTENT(IN) :: den
TYPE(t_potden), dimension(3), INTENT(INOUT) :: gradphi
TYPE(t_potden) :: denloc
TYPE(t_gradients) :: grad
REAL, ALLOCATABLE :: thet(:), phi(:), grad_temp(:, :, :)
......@@ -205,9 +205,15 @@ CONTAINS
ALLOCATE (grad_temp(atoms%jri(n)*nsp,1,3))
ALLOCATE (thet(atoms%nsp()),phi(atoms%nsp()))
denloc=den
DO lh=0, lhmax
denloc%mt(:,lh,n,1) = denloc%mt(:,lh,n,1)*atoms%rmsh(:, n)**2
END DO ! lh
CALL init_mt_grid(1, atoms, sphhar, .TRUE., sym, thet, phi)
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, den%mt(:,0:,n,:), n, grad)
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, denloc%mt(:,0:,n,:), n, grad)
kt = 0
DO jr = 1, atoms%jri(n)
......@@ -226,6 +232,9 @@ CONTAINS
CALL mt_from_grid(atoms, sphhar, n, 1, grad_temp(:,:,i), gradphi(i)%mt(:,0:,n,:))
DO lh=0, lhmax
gradphi(i)%mt(:,lh,n,1) = gradphi(i)%mt(:,lh,n,1)*atoms%rmsh(:, n)**2
IF ((sphhar%llh(lh,1)/=0).AND.(sphhar%llh(lh,1)/=2)) THEN
gradphi(i)%mt(:,lh,n,1) = 0.0
END IF
END DO ! lh
END DO ! i
......
......@@ -11,7 +11,7 @@ module m_vgen_coulomb
contains
subroutine vgen_coulomb( ispin, mpi, dimension, oneD, input, field, vacuum, sym, stars, &
cell, sphhar, atoms, den, vCoul, results )
cell, sphhar, atoms, dosf, den, vCoul, results )
!----------------------------------------------------------------------------
! FLAPW potential generator
!----------------------------------------------------------------------------
......@@ -50,6 +50,7 @@ contains
type(t_cell), intent(in) :: cell
type(t_sphhar), intent(in) :: sphhar
type(t_atoms), intent(in) :: atoms
LOGICAL, INTENT(IN) :: dosf
type(t_potden), intent(in) :: den
type(t_potden), intent(inout) :: vCoul
type(t_results), intent(inout), optional :: results
......@@ -174,7 +175,9 @@ contains
call timestop("interstitial")
end if ! mpi%irank == 0
if (dosf) then
vCoul%pw(:,:)=0.0
end do
! MUFFIN-TIN POTENTIAL
call timestart( "MT-spheres" )
......
......@@ -74,8 +74,6 @@ contains
#endif
integer :: OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM
! SPHERE BOUNDARY CONTRIBUTION to the coefficients calculated from the values
! of the interstitial Coulomb / Yukawa potential on the sphere boundary
......
......@@ -17,20 +17,22 @@ MODULE m_xcBfield
! consistent source-free fields.
!-----------------------------------------------------------------------------
PUBLIC :: makeBxc, sourcefree, builddivtest
PUBLIC :: makeBxc, sourcefree, builddivtest ! TODO: Build a routine to pack A_vec
! back into the matrix correctly.
CONTAINS
SUBROUTINE makeBxc(stars,atoms,sphhar,vacuum,input,noco,vTot,bxc)
SUBROUTINE makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denmat,factor,aVec)
! Contructs the exchange-correlation magnetic field from the total poten-
! tial matrix. Only used for the implementation of source free fields and
! therefore only applicable in a fully non-collinear description of magne-
! tism.
! tial matrix or the magnetic density for the density matrix. Only used for
! the implementation of source free fields and therefore only applicable in
! a fully non-collinear description of magnetism.
!
! Assumes the following form for the potential matrix:
! V_mat = V*Id_(2x2) + sigma_vec*B_vec
! Assumes the following form for the density/potential matrix:
! rho_mat = n*Id_(2x2) + conj(sigma_vec)*m_vec
! V_mat = V*Id_(2x2) + sigma_vec*B_vec
!
! B_vec is saved as a density type with an additional r^2-factor.
! A_vec is saved as a density type with an additional r^2-factor.
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_atoms), INTENT(IN) :: atoms
......@@ -38,8 +40,9 @@ CONTAINS
TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_input), INTENT(IN) :: input
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), INTENT(IN) :: vTot
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: bxc
TYPE(t_potden), INTENT(IN) :: denmat
REAL, INTENT(IN) :: factor
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: aVec
TYPE(t_potden), DIMENSION(4) :: dummyDen
INTEGER :: i, itype, ir
......@@ -54,28 +57,31 @@ CONTAINS
ALLOCATE(dummyDen(i)%pw_w,mold=dummyDen(i)%pw)
ENDDO
CALL matrixsplit(stars,atoms,sphhar,vacuum,input,noco,2.0,vtot, &
CALL matrixsplit(stars,atoms,sphhar,vacuum,input,noco,factor,denmat, &
dummyDen(1),dummyDen(2),dummyDen(3),dummyDen(4))
r2=1.0
! Initialize and fill the magnetic field.
DO itype=1,atoms%ntype
DO ir=1,atoms%jri(itype)
r2=atoms%rmsh(ir,itype)**2
DO i=1,3
CALL bxc(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd, &
atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE., &
POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(bxc(i)%pw_w,mold=bxc(i)%pw)
bxc(i)%mt(:,0:,:,:) = dummyDen(i+1)%mt(:,0:,:,:)
bxc(i)%pw(1:,:) = dummyDen(i+1)%pw(1:,:)
bxc(i)%vacz(1:,1:,:) = dummyDen(i+1)%vacz(1:,1:,:)
bxc(i)%vacxy(1:,1:,1:,:) = dummyDen(i+1)%vacxy(1:,1:,1:,:)
ENDDO
END DO
END DO
! Initialize and fill the vector field.
DO i=1,3
CALL aVec(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd, &
atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE., &
POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(aVec(i)%pw_w,mold=aVec(i)%pw)
DO itype=1,atoms%ntype
DO ir=1,atoms%jri(itype)
IF (factor==2.0) THEN
r2=atoms%rmsh(ir,itype)**2
END IF
aVec(i)%mt(ir,0:,itype,:) = dummyDen(i+1)%mt(ir,0:,itype,:)*r2
END DO !ir
END DO !itype
aVec(i)%pw(1:,:) = dummyDen(i+1)%pw(1:,:)
aVec(i)%vacz(1:,1:,:) = dummyDen(i+1)%vacz(1:,1:,:)
aVec(i)%vacxy(1:,1:,1:,:) = dummyDen(i+1)%vacxy(1:,1:,1:,:)
END DO !i
END SUBROUTINE makeBxc
END SUBROUTINE makeVectorField
SUBROUTINE sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,bxc)
USE m_vgen_coulomb
......@@ -87,6 +93,9 @@ CONTAINS
! b) Solve the Poisson equation (nabla_vec*nabla_vec)phi=-4*pi*d.
! c) Construct an auxiliary vector field C_vec=(nabla_vec phi)/(4*pi).
! d) Build A_vec_sf=A_vec+C_vec, which is source free by construction.
!
! Note: A_vec is assumed to be a density with an additional factor of r^2.
! A field of the same form will also be calculated.
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: dimension
......@@ -100,12 +109,12 @@ CONTAINS
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), DIMENSION(3), INTENT(IN) :: bxc
TYPE(t_potden), DIMENSION(3), INTENT(INOUT) :: bxc
TYPE(t_atoms) :: atloc
TYPE(t_potden) :: div,phi,checkdiv
TYPE(t_potden), DIMENSION(3) :: cvec, corrB
INTEGER :: n, jr
INTEGER :: n, jr, lh, lhmax
CALL div%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype, &
atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_DEN, &
......@@ -113,33 +122,41 @@ CONTAINS
ALLOCATE(div%pw_w,mold=div%pw)
CALL divergence(stars,atoms,sphhar,vacuum,sym,cell,noco,bxc,div)
div%mt(:,2:,:,:)=0.0
div%mt(:,0,:,:)=0.0
atloc=atoms
atloc%zatom=0.0
atloc%zatom=0.0 !Local atoms variable with no charges; needed for the potential generation.
eps=1.e-10
CALL phi%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTCOUL,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(phi%pw_w(SIZE(phi%pw,1),size(phi%pw,2)))
phi%pw_w = CMPLX(0.0,0.0)
CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atloc,div,phi)
DO n=1, atoms%ntype
CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atloc,.TRUE.,div,phi)
DO n=1,atoms%ntype
lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
DO lh=0, lhmax
phi%mt(:,lh,n,1) = phi%mt(:,lh,n,1)*atoms%rmsh(:, n)**2
DO jr = 1, atoms%jri(n)
DO lh=0, lhmax
!IF (ABS(phi%mt(jr,lh,n,1))<eps) THEN
IF (lh/=1) THEN
phi%mt(jr,lh,n,:)=0.0
END IF
END DO
END DO
END DO
DO i=1,3
CALL cvec(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
CALL cvec(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_DEN,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(cvec(i)%pw_w,mold=cvec(i)%pw)
ENDDO
CALL divpotgrad(stars,atoms,sphhar,vacuum,sym,cell,noco,phi,cvec)
DO i=1,3
CALL corrB(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
CALL corrB(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_DEN,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(corrB(i)%pw_w,mold=corrB(i)%pw)
CALL corrB(i)%addPotDen(bxc(i),cvec(i))
ENDDO
......@@ -151,9 +168,9 @@ CONTAINS
CALL divergence(stars,atoms,sphhar,vacuum,sym,cell,noco,corrB,checkdiv)
!CALL plotBtest(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
! noco, div, phi, cvec(1), cvec(2), cvec(3), &
! corrB(1), corrB(2), corrB(3), checkdiv)
CALL plotBtest(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
noco, div, phi, cvec(1), cvec(2), cvec(3), &
corrB(1), corrB(2), corrB(3), checkdiv)
END SUBROUTINE sourcefree
......@@ -226,9 +243,9 @@ CONTAINS
!print *, th/pi_const
!print *, ph/pi_const
!(r/atoms%rmt(n))
A_temp(kt+k,1,1)=(r**2)*r*SIN(th)*COS(ph)
A_temp(kt+k,2,1)=(r**2)*r*SIN(th)*SIN(ph)
A_temp(kt+k,3,1)=(r**2)*r*COS(th)
A_temp(kt+k,1,1)=(r**2)*r**2!*SIN(th)*COS(ph)
A_temp(kt+k,2,1)=(r**2)*r**2!*SIN(th)*SIN(ph)
A_temp(kt+k,3,1)=(r**2)*r**2!*COS(th)
ENDDO ! k
kt = kt + nsp
END DO ! ir
......@@ -286,10 +303,6 @@ CONTAINS
END DO
END DO !z-loop
!DO i=1,ifftxc3
!A_temp(i,:,:)=0.5
!END DO
DO i=1,3
CALL pw_from_grid(.TRUE.,stars,.TRUE.,A_temp(:,i,:),Avec(i)%pw,Avec(i)%pw_w)
END DO
......
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