Commit 0725c6e9 authored by Daniel Wortmann's avatar Daniel Wortmann

Starting to refactor potential generator, nothing compiles yet...

parent 294af201
......@@ -113,7 +113,7 @@ CONTAINS
ALLOCATE ( kpts%bk(3,kpts%nkpt),kpts%wtkpt(kpts%nkpt) )
ALLOCATE ( stars%pgfft(0:DIMENSION%nn3d-1),stars%pgfft2(0:DIMENSION%nn2d-1) )
ALLOCATE ( stars%ufft(0:27*stars%mx1*stars%mx2*stars%mx3-1) )
ALLOCATE ( atoms%bmu(atoms%ntype),atoms%vr0(atoms%ntype) )
ALLOCATE ( atoms%bmu(atoms%ntype) )
ALLOCATE ( atoms%l_geo(atoms%ntype) )
ALLOCATE ( atoms%nlo(atoms%ntype),atoms%llo(atoms%nlod,atoms%ntype) )
ALLOCATE ( atoms%lo1l(0:atoms%llod,atoms%ntype),atoms%nlol(0:atoms%llod,atoms%ntype),atoms%lapw_l(atoms%ntype) )
......@@ -140,7 +140,6 @@ CONTAINS
input%l_coreSpec = .FALSE.
atoms%vr0(:) = 0.0
results%force(:,:,:) = 0.0
IF(.NOT.juDFT_was_argument("-toXML")) THEN
......
......@@ -275,8 +275,6 @@ SUBROUTINE postprocessInput(mpi,input,sym,stars,atoms,vacuum,obsolete,kpts,&
dimension%msh = 0
ALLOCATE(atoms%rmsh(atoms%jmtd,atoms%ntype))
ALLOCATE(atoms%volmts(atoms%ntype))
ALLOCATE(atoms%vr0(atoms%ntype)) ! This should actually not be in the atoms type!
atoms%vr0(:) = 0.0
na = 0
DO iType = 1, atoms%ntype
l_vca = .FALSE.
......
......@@ -177,7 +177,7 @@ CONTAINS
WRITE (6,FMT=8045) zintn_r(n)
WRITE (16,FMT=8045) zintn_r(n)
CALL intgr3(mt(1,n),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),totz)
vmd(n) = atoms%rmt(n)*atoms%vr0(n)/sfp_const + atoms%zatom(n) - totz*sfp_const
vmd(n) = atoms%rmt(n)*vCoul%mt(atoms%jri(n),0,n,1)/sfp_const + atoms%zatom(n) - totz*sfp_const
vmd(n) = -atoms%neq(n)*atoms%zatom(n)*vmd(n)/ (2.*atoms%rmt(n))
WRITE (6,FMT=8050) n,vmd(n)
WRITE (16,FMT=8050) n,vmd(n)
......
......@@ -120,6 +120,10 @@ CONTAINS
! units: hartrees
!
WRITE (6,FMT=8000)
8000 FORMAT (/,/,t10,' p o t e n t i a l g e n e r a t o r',/)
CALL vTot%resetPotDen()
CALL vCoul%resetPotDen()
CALL vx%resetPotDen()
......@@ -135,18 +139,7 @@ CONTAINS
workDen = den
IF (mpi%irank == 0) THEN
!
! -- total = .false. and reap = .false. means, that we now calculate
! the potential from the output density
!
IF ((.NOT.input%total).AND.(.NOT.reap)) THEN
IF (noco%l_noco) THEN
CALL juDFT_error("vgen:1",calledby ="vgen")
ENDIF
END IF
WRITE (6,FMT=8000)
8000 FORMAT (/,/,t10,' p o t e n t i a l g e n e r a t o r',/)
vxpw_w = CMPLX(0.,0.)
! ---> perform spin summation of charge densities
! ---> for the calculation of the coulomb potentials
......@@ -160,231 +153,14 @@ CONTAINS
END IF
!
! ************** coulomb potential ***********************
IF (l_xyav) THEN
zat_l(:) = 0. ! for xy-averaged densities do not
ELSE ! include nuclear charge
zat_l(:) = atoms%zatom(:)
ENDIF
! ----> create pesudo-charge density coefficients
ENDIF ! (mpi%irank == 0)
#ifdef CPP_MPI
CALL MPI_BCAST(workDen%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
#endif
CALL timestart("psqpw")
CALL psqpw(mpi, atoms,sphhar,stars,vacuum, dimension,cell,input,sym,oneD,&
workDen%pw,workDen%mt,workDen%vacz,l_xyav, psq)
CALL timestop("psqpw")
IF (mpi%irank == 0) THEN
IF (l_xyav) THEN ! write out xy-averaged density & stop
CALL xy_av_den(stars,vacuum, cell,psq,workDen%vacz)
CALL juDFT_error("xy-averaged density calculated", calledby ="vgen")
ENDIF
! ------------------------------------------
IF (oneD%odi%d1) THEN
!-odim
CALL timestart("coulomb potential")
!---> generates the m=0,gz=0 component of the vacuum potential
CALL od_vvac(stars,vacuum,cell, psq,workDen%vacz, vTot%vacz)
!---> generation of the vacuum warped potential components and
!---> interstitial pw potential
!---> vvacxy_5.F is a symmetrized potential generator
CALL od_vvacis(oneD%odi%n2d,dimension,vacuum,oneD%odi%nq2,&
oneD%odi%kv,cell,oneD%odi%M,stars,oneD%odi%nst2,&
oneD, workDen%vacz,workDen%vacxy,psq,vTot%vacz,sym, vTot%vacxy,vTot%pw)
CALL timestop("coulomb potential")
CALL vgen_coulomb(1,mpi,DIMENSION,oneD,input,vacuum,sym,stars,cell,sphhar,atoms,den,vCoul,results)
!+odim
ELSEIF (input%film .AND. .NOT.oneD%odi%d1) THEN
! ----> potential in the vacuum region
!
CALL timestart("p vac")
CALL vvac(vacuum,stars, cell,sym,input, psq,workDen%vacz, vTot%vacz,rhobar,sig1dh,vz1dh)
CALL vvacis(stars,vacuum, sym,cell, psq, input, vTot%vacxy)
CALL vvacxy(stars,vacuum,cell,sym,input, workDen%vacxy, vTot%vacxy, alphm)
CALL timestop("p vac")
END IF
! ------------------------------------------
! ----> potential in the interstitial region
CALL timestart("p int")
WRITE (6,FMT=8010)
8010 FORMAT (/,5x,'coulomb potential in the interstitial region:')
IF (input%film .AND. .NOT.oneD%odi%d1) THEN
! -----> create v(z) for each 2-d reciprocal vector
ivfft = 3*stars%mx3
! ivfft = 2*mx3 - 1
ani = 1.0/REAL(ivfft)
DO irec2 = 1,stars%ng2
i = 0
DO i3 = 0,ivfft - 1
i = i + 1
z = cell%amat(3,3)*i3*ani
IF (z.GT.cell%amat(3,3)/2.) z = z - cell%amat(3,3)
vintcza = vintcz(stars,vacuum,cell,sym,input,&
z,irec2, psq,vTot%vacxy,vTot%vacz,rhobar,sig1dh,vz1dh,alphm)
af1(i) = REAL(vintcza)
bf1(i) = AIMAG(vintcza)
ENDDO
! z = (i_sm-1)*ani
! IF (z > 0.5) z = z - 1.0
! af1(i_sm) = af1(i_sm) + z * delta
! bf1(i_sm) = bf1(i_sm) + z * deltb
! ENDDO
! ENDIF
! --> 1-d fourier transform and store the coefficients in vTot%pw( ,1)
CALL cfft(af1,bf1,ivfft,ivfft,ivfft,-1)
! delta = ivfft * delta * 2 / fpi ! * amat(3,3)**2 * ani
i = 0
DO i3 = 0,ivfft - 1
k3 = i3
IF (k3 > FLOOR(ivfft/2.) ) k3 = k3 - ivfft
i = i + 1
IF ((k3.GE.-stars%mx3).AND.(k3.LE.stars%mx3)) THEN
irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
! IF ( (irec2 == 1).AND.(i3 > 0) ) THEN ! smooth potential
! corr = 2.0*mod(abs(k3),2) - 1.0
! bf1(i) = bf1(i) + delta * corr / k3
! ENDIF
! ----> only stars within g_max sphere (shz oct.97)
IF (irec3.NE.0) THEN
!
xint = CMPLX(af1(i),bf1(i))*ani
nzst1 = stars%nstr(irec3)/stars%nstr2(irec2)
vTot%pw(irec3,1) = vTot%pw(irec3,1) + xint/nzst1
END IF
ENDIF
ENDDO
ENDDO
ELSEIF (.NOT.input%film) THEN
vTot%pw(1,1) = CMPLX(0.0,0.0)
vTot%pw(2:stars%ng3,1)=fpi_const*psq(2:stars%ng3)/(stars%sk3(2:stars%ng3)*stars%sk3(2:stars%ng3))
END IF
CALL timestop("p int")
ENDIF ! mpi%irank == 0
! --------------------------------------------
! ---> potential in the muffin-tin spheres
CALL timestart("p vmts")
CALL vmts(mpi, stars,sphhar,atoms, sym,cell,oneD, vTot%pw,workDen%mt, vTot%mt)
! --------------------------------------------
CALL timestop("p vmts")
IF (mpi%irank == 0) THEN
! ---> check continuity of coulomb potential
IF (input%vchk) THEN
CALL timestart("checking")
CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,&
cell,vTot,1)
CALL timestop("checking")
END IF
!
!========TOTAL==============================================
!
! IF (l_xyav) THEN ! write out xy-averaged potential & stop
! CALL xy_av_den(
! > n3d,k3d,nq3,nmzd,nmz,dvac,delz,
! > area,ig2,kv3,amat,vTot%pw,vTot%vacz(1,1,1))
! CALL juDFT_error("xy-averaged potential calculated",calledby="vgen")
! ENDIF
IF (input%total) THEN
CALL timestart("int_nv")
!
! ---> AVERAGE COULOMB POTENTIAL ON THE SPHERE
! FOR CALCULATING THE MADELUNG TERM in totale.f
! r=Rmt
DO n=1,atoms%ntype
atoms%vr0(n)=vTot%mt(atoms%jri(n),0,n,1)
ENDDO
!
! CALCULATE THE INTEGRAL OF n*Vcoulomb
!
WRITE (6,FMT=8020)
WRITE (16,FMT=8020)
8020 FORMAT (/,10x,'density-coulomb potential integrals',/)
!
! interstitial first
!
! convolute ufft and pot: F(G) = \sum_(G') U(G - G') V(G')
!
CALL convol(stars, vpw_w, vTot%pw, stars%ufft)
!
IF (input%jspins.EQ.2) CALL CPP_BLAS_ccopy(stars%ng3,vpw_w(1,1),1,vpw_w(1,input%jspins),1)
!
results%te_vcoul = 0.0
CALL int_nv(stars,vacuum,atoms,sphhar, cell,sym,input,oneD,&
workDen%pw,vpw_w, workDen%vacxy,vTot%vacxy, workDen%vacz,vTot%vacz, workDen%mt,vTot%mt, results%te_vcoul)
WRITE (6,FMT=8030) results%te_vcoul
WRITE (16,FMT=8030) results%te_vcoul
8030 FORMAT (/,10x,'total density-coulomb potential integral :', t40,f20.10)
CALL timestop("int_nv")
INQUIRE(file='vdW_kernel_table',exist=l_vdw)
IF (l_vdw) THEN
CALL timestart("fleur_vdW")
! calculate vdW contribution to potential
CALL fleur_vdW(mpi,atoms,sphhar,stars, input,dimension,&
cell,sym,oneD,vacuum, workDen%pw(:,1),workDen%mt(:,:,:,1), vpw_w(:,1),vTot%mt(:,:,:,:))
CALL timestop("fleur_vdW")
ENDIF
END IF
!
!==========END TOTAL===================================================
!
! ----> reload the density for calculating vxc (for spin-pol. case)
!
IF (input%jspins.EQ.2) THEN
workDen = den
vTot%mt(:,0:,:,2) = vTot%mt(:,0:,:,1)
vTot%pw(:,2) = vTot%pw(:,1)
IF (input%film) THEN
vTot%vacxy(:,:,:,2) = vTot%vacxy(:,:,:,1)
vTot%vacz(:,:,2)=vTot%vacz(:,:,1)
END IF
END IF
IF (input%total) THEN
DO js = 1,input%jspins
! to enable a GW calculation,
vpw_w(1:stars%ng3,js)=vpw_w(1:stars%ng3,js)/stars%nstr(1:stars%ng3) ! the PW-coulomb part is not
! used otherwise anyway.
ENDDO
vCoul%iter = vTot%iter
vCoul%mt = vTot%mt
vCoul%pw(:,1:input%jspins) = vpw_w(:,1:input%jspins)
vCoul%vacz = vTot%vacz
vCoul%vacxy = vTot%vacxy
DO js = 1,input%jspins
DO i = 1,stars%ng3
vpw_w(i,js)=vpw_w(i,js)*stars%nstr(i)
ENDDO
ENDDO
END IF
IF (sliceplot%plpot) THEN
OPEN (11,file='potcoul_pl',form='unformatted',status='unknown')
CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
11, vTot%iter,vTot%mt,vTot%pw,vTot%vacz,vTot%vacxy)
CLOSE(11)
END IF
ENDIF !irank==0
CALL vCoul%copy_both_spin(vTot)
! ******** exchange correlation potential******************
!+ta
......
......@@ -9,7 +9,7 @@ MODULE m_types_potden
TYPE t_potden
INTEGER :: iter
INTEGER :: potdenType
COMPLEX,ALLOCATABLE :: pw(:,:)
COMPLEX,ALLOCATABLE :: pw(:,:),pw_w(:,:)
REAL,ALLOCATABLE :: mt(:,:,:,:)
REAL,ALLOCATABLE :: vacz(:,:,:)
COMPLEX,ALLOCATABLE :: vacxy(:,:,:,:)
......@@ -33,9 +33,24 @@ MODULE m_types_potden
PROCEDURE :: init_potden_simple
PROCEDURE :: resetpotden
GENERIC :: init=>init_potden_types,init_potden_simple
PROCEDURE :: copy_both_spin
END TYPE t_potden
CONTAINS
SUBROUTINE copy_both_spin(this,that)
IMPLICIT NONE
CLASS(t_potden),INTENT(IN) :: this
TYPE(t_potden),INTENT(INOUT) :: that
IF (SIZE(that%mt,4)==2) THEN
that%mt(:,0:,:,2)=this%mt(:,0:,:,1)
that%pw(:,2)=this%pw(:,1)
that%vacz(:,:,2)=this%vacz(:,:,1)
that%vacxy(:,:,:,2)=this%vacxy(:,:,:,1)
IF (ALLOCATED(that%pw_w).AND.ALLOCATED(this%pw_w)) that%pw_w(:,2)=this%pw_w(:,1)
END IF
END SUBROUTINE copy_both_spin
SUBROUTINE init_potden_types(pd,stars,atoms,sphhar,vacuum,noco,oneD,jspins,nocoExtraDim,potden_type)
USE m_judft
USE m_types_setup
......
......@@ -180,7 +180,6 @@ MODULE m_types_setup
TYPE(t_utype),ALLOCATABLE::lda_u(:)
INTEGER,ALLOCATABLE :: relax(:,:) !<(3,ntype)
INTEGER, ALLOCATABLE :: nflip(:) !<flip magnetisation of this atom
REAL,ALLOCATABLE:: vr0(:) !< Average Coulomb potential for atoms
END TYPE t_atoms
TYPE t_cell
......
......@@ -25,6 +25,7 @@ vgen/mpmom.F90
vgen/od_vvac.f90
vgen/od_vvacis.f90
vgen/pot_mod.f90
vgen/vgen_coulomb.F90
vgen/prp_xcfft_map.f90
vgen/psqpw.F90
vgen/rhodirgen.f90
......
......@@ -4,21 +4,16 @@ MODULE m_intnv
! and potential in the unit cell
! ************************************************
CONTAINS
SUBROUTINE int_nv(&
& stars,vacuum,atoms,sphhar,&
& cell,sym,input,oneD,&
& qpw,vpw_w,&
& rhtxy,vxy,&
& rht,vz,&
& rho,vr,&
& RESULT)
SUBROUTINE int_nv(ispin,stars,vacuum,atoms,sphhar,&
cell,sym,input,oneD,vpot,den,RESULT)
USE m_intgr, ONLY : intgr3,intgz0
USE m_types
IMPLICIT NONE
! ..
! .. Scalar Arguments ..
REAL RESULT
REAL RESULT
INTEGER,INTENT(IN) :: ispin
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_atoms),INTENT(IN) :: atoms
......@@ -27,14 +22,9 @@ CONTAINS
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_input),INTENT(IN) :: input
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_potden),INTENT(IN) :: vpot,den
! ..
! .. Array Arguments ..
COMPLEX qpw(stars%ng3),rhtxy(vacuum%nmzxyd,oneD%odi%n2d-1,2),&
& vpw_w(stars%ng3),vxy(vacuum%nmzxyd,oneD%odi%n2d-1,2)
REAL rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype),rht(vacuum%nmzd,2),&
& vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype),vz(vacuum%nmzd,2)
! ..
! .. Local Scalars ..
REAL dpdot,facv,tis,tmt,tvac,tvact
......@@ -51,7 +41,7 @@ CONTAINS
!
! -> warping has been moved to vgen and visxc resp. ...gustav
!
tis = cell%omtil * REAL( DOT_PRODUCT(vpw_w,qpw))
tis = cell%omtil * REAL( DOT_PRODUCT(vpot%pw_w(:,ispin),den%pw(:,ispin)))
WRITE (6,FMT=8020) tis
WRITE (16,FMT=8020) tis
......@@ -66,7 +56,7 @@ CONTAINS
DO n = 1,atoms%ntype
DO lh = 0,sphhar%nlh(atoms%ntypsy(nat))
DO j = 1,atoms%jri(n)
dpj(j) = rho(j,lh,n)*vr(j,lh,n)
dpj(j) = den%mt(j,lh,n,ispin)*vpot%mt(j,lh,n,ispin)
ENDDO
CALL intgr3(dpj,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),dpdot)
tmt = tmt + dpdot*atoms%neq(n)
......@@ -94,14 +84,14 @@ CONTAINS
dpz=0.0
DO ivac = 1,vacuum%nvac
DO ip = 1,vacuum%nmz
dpz(npz-ip) = rht(ip,ivac)*vz(ip,ivac)
dpz(npz-ip) = den%vacz(ip,ivac,ispin)*vpot%vacz(ip,ivac,ispin)
! ---> WARPING REGION
ENDDO
DO k2 = 2,stars%ng2
DO ip = 1,vacuum%nmzxy
dpz(npz-ip) = dpz(npz-ip) +&
& stars%nstr2(k2)*rhtxy(ip,k2-1,ivac)*&
& CONJG(vxy(ip,k2-1,ivac))
& stars%nstr2(k2)*den%vacxy(ip,k2-1,ivac,ispin)*&
& CONJG(vpot%vacxy(ip,k2-1,ivac,ispin))
ENDDO
ENDDO
CALL intgz0(dpz,vacuum%delz,vacuum%nmz,tvac,tail)
......@@ -121,15 +111,15 @@ CONTAINS
dpz=0.0
DO ip = 1,vacuum%nmz
dpz(npz-ip) = (cell%z1+vacuum%delz*(ip-1))*&
& rht(ip,vacuum%nvac)*vz(ip,vacuum%nvac)
& den%vacz(ip,vacuum%nvac,ispin)*vpot%vacz(ip,vacuum%nvac,ispin)
! ---> WARPING REGION
ENDDO
DO k2 = 2,oneD%odi%nq2
DO ip = 1,vacuum%nmzxy
dpz(npz-ip) = dpz(npz-ip)+&
& (cell%z1+vacuum%delz*(ip-1))*&
& rhtxy(ip,k2-1,vacuum%nvac)*&
& CONJG(vxy(ip,k2-1,vacuum%nvac))
& den%vacxy(ip,k2-1,vacuum%nvac,ispin)*&
& CONJG(vpot%vacxy(ip,k2-1,vacuum%nvac,ispin))
ENDDO
ENDDO
......
......@@ -28,8 +28,8 @@ CONTAINS
!
!
! .. Array Arguments ..
REAL, INTENT (IN) :: rho(:,0:,:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd)
COMPLEX, INTENT (IN) :: qpw(:,:) !(stars%ng3,dimension%jspd)
REAL, INTENT (IN) :: rho(:,0:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
COMPLEX, INTENT (IN) :: qpw(:) !(stars%ng3)
COMPLEX, INTENT (OUT):: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
!-odim
!+odim
......@@ -50,12 +50,12 @@ CONTAINS
!
IF (mpi%irank == 0) THEN
CALL mt_moments(atoms,sphhar,&
& rho(:,:,:,1),qlmo)
& rho(:,:,:),qlmo)
ENDIF ! mpi%irank == 0
CALL pw_moments(mpi,stars,atoms,cell,&
& sym,oneD,&
& qpw(:,1),qlmp)
& qpw(:),qlmp)
!
! eq.(15): \tilde q_(lm}^i = q_{lm}^i - q_{lm}^{Ii}
!
......
......@@ -29,8 +29,8 @@ CONTAINS
COMPLEX, INTENT (IN) :: psq(stars%ng3)
REAL, INTENT (IN) :: rht(:,:,:) !(vacuum%nmzd,2,dimension%jspd)
REAL, INTENT (OUT) :: vz(:,:,:) !(vacuum%nmzd,2,dimension%jspd)
REAL, INTENT (IN) :: rht(:,:) !(vacuum%nmzd,2)
REAL, INTENT (OUT) :: vz(:,:) !(vacuum%nmzd,2)
COMPLEX rhobar
INTEGER k1,k2,irec3,irec2,i,j,ivac,imz,imz1
......@@ -46,7 +46,7 @@ CONTAINS
f2(i) = 0.
f22(i) = 0.
DO ivac = 1,vacuum%nvac
vz(i,ivac,1) = 0.
vz(i,ivac) = 0.
END DO
END DO
......@@ -71,7 +71,7 @@ CONTAINS
!----> 1st equivalent way
DO i=1,vacuum%nmz
rht1(i) = fpi_const*(cell%z1+(i-1)*vacuum%delz)*rht(i,1,1)
rht1(i) = fpi_const*(cell%z1+(i-1)*vacuum%delz)*rht(i,1)
ENDDO
CALL qsf(vacuum%delz,rht1(1),f2(1),vacuum%nmz,1)
......@@ -89,7 +89,7 @@ CONTAINS
ENDDO
CALL qsf(vacuum%delz,f22(1),a,vacuum%nmz,0)
DO ivac =1,vacuum%nvac
vz(i,ivac,1) = -a(1)
vz(i,ivac) = -a(1)
ENDDO
ENDDO
!----> 2nd equivalent way (via the Green function)
......@@ -99,14 +99,14 @@ CONTAINS
DO imz1 = 1,vacuum%nmz
zp = cell%z1 + (imz1-1)*vacuum%delz
IF (imz1.LE.imz) THEN
rht1(imz1) = fpi_const*LOG(z)*zp*rht(imz1,1,1)
rht1(imz1) = fpi_const*LOG(z)*zp*rht(imz1,1)
ELSE
rht1(imz1) = fpi_const*LOG(zp)*zp*rht(imz1,1,1)
rht1(imz1) = fpi_const*LOG(zp)*zp*rht(imz1,1)
END IF
END DO
CALL qsf(vacuum%delz,rht1,a,vacuum%nmz,0)
vz(imz,1,1) = tpi_const*LOG(z)*(cell%z1*cell%z1)*rhobar - a(1)
vz(imz,1) = tpi_const*LOG(z)*(cell%z1*cell%z1)*rhobar - a(1)
END DO
RETURN
......
......@@ -63,11 +63,11 @@ CONTAINS
INTEGER, INTENT (IN) :: nstr2_1(n2d_1)
INTEGER, INTENT (IN) :: kv2_1(2,n2d_1)
COMPLEX, INTENT (INOUT) :: psq(stars%ng3)
REAL, INTENT (IN) :: vz(vacuum%nmzd,2,DIMENSION%jspd)
REAL, INTENT (IN) :: rht(vacuum%nmzd,2,DIMENSION%jspd)
COMPLEX, INTENT (IN) :: rhtxy(vacuum%nmzxyd,n2d_1-1,2,DIMENSION%jspd)
COMPLEX, INTENT (OUT):: vxy(vacuum%nmzxyd,n2d_1-1,2,DIMENSION%jspd)
COMPLEX, INTENT (OUT):: vpw(stars%ng3,DIMENSION%jspd)
REAL, INTENT (IN) :: vz(vacuum%nmzd,2)
REAL, INTENT (IN) :: rht(vacuum%nmzd,2)
COMPLEX, INTENT (IN) :: rhtxy(vacuum%nmzxyd,n2d_1-1,2)
COMPLEX, INTENT (OUT):: vxy(vacuum%nmzxyd,n2d_1-1,2)
COMPLEX, INTENT (OUT):: vpw(stars%ng3)
! local
INTEGER :: m
......@@ -174,19 +174,19 @@ CONTAINS
!----> vpw in the '1st aproximation' (V - tilde)
vpw(1,1) = CMPLX(0.,0.)
vpw(1) = CMPLX(0.,0.)
DO irec3 = 2,stars%ng3
g = stars%sk3(irec3)
vpw(irec3,1) = fpi_const*psq(irec3)/(g*g)
vpw(irec3) = fpi_const*psq(irec3)/(g*g)
ENDDO
DO irc1 = 2,nq2_1
DO i = 1,vacuum%nmzxy
vxy(i,irc1-1,1,1) = CMPLX(0.,0.)
vxy(i,irc1-1,1) = CMPLX(0.,0.)
END DO
END DO
......@@ -219,7 +219,7 @@ CONTAINS
irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
IF (irec3.NE.0) THEN
val(irc1) = val(irc1) +&
& (ic**m)*vpw(irec3,1)*EXP(-ic*&
& (ic**m)*vpw(irec3)*EXP(-ic*&
& m*phi)*fJJ(iabs(m),irec2)*&
& stars%nstr2(irec2)*mult
END IF
......@@ -268,7 +268,7 @@ CONTAINS
irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),gzi)
IF (irec3.NE.0) THEN
val_m(gzi,m) = val_m(gzi,m) +&
& (ic**m)*vpw(irec3,1)*EXP(-ic*&
& (ic**m)*vpw(irec3)*EXP(-ic*&
& m*phi)*fJJ(iabs(m),irec2)*&
& stars%nstr2(irec2)*mult
END IF
......@@ -290,12 +290,12 @@ CONTAINS
im = zf
q = zf - im
vis(ix,iy,1) = 0.5*(q-1.)*&
& (q-2.)*vz(im,1,1) -&
& q*(q-2.)*vz(im+1,1,1) +&
& 0.5*q*(q-1.)*vz(im+2,1,1)
& (q-2.)*vz(im,1) -&
& q*(q-2.)*vz(im+1,1) +&
& 0.5*q*(q-1.)*vz(im+2,1)
ELSE
vis(ix,iy,1) = &
& vz(1,1,1) - val(1) + tpi_const*&
& vz(1,1) - val(1) + tpi_const*&
& psq(1)*(cell%z1*cell%z1 - r*r)/2.
END IF
DO irc1 = 2,nq2_1
......@@ -398,7 +398,7 @@ CONTAINS
!----- this form of the density is just more easy to use
DO imz = 1,vacuum%nmzxy
rxy(imz) = rhtxy(imz,irec1(l)-1,1,1)
rxy(imz) = rhtxy(imz,irec1(l)-1,1)
END DO
!----- vacuum potential caused by the vacuum density
......@@ -435,7 +435,7 @@ CONTAINS
pint(:vacuum%nmzxy) = fact(:vacuum%nmzxy)*aa
vxy(:vacuum%nmzxy,irec1(l)-1,1,1) = pvac(:vacuum%nmzxy) + pint(:vacuum%nmzxy)
vxy(:vacuum%nmzxy,irec1(l)-1,1) = pvac(:vacuum%nmzxy) + pint(:vacuum%nmzxy)
!----- array val further is a boundary values of the
!----- potential V- \tilde \tilde which is created to compensate
......@@ -444,7 +444,7 @@ CONTAINS
!----- density, V - \tilde and V - \tilde\tilde are then added in
!----- order to obtain the real interstitial potential
val_help = vxy(1,irec1(l)-1,1,1) - val(irec1(l))
val_help = vxy(1,irec1(l)-1,1) - val(irec1(l))
!----- potential \tilde\tilde{V} is a solution of the Laplase equation
!----- in the interstitial with the boundary conditions val_0 and val_z
......@@ -471,7 +471,7 @@ CONTAINS
!------------------------------------------------------------->
DO imz = 1,vacuum%nmzxy
rxy(imz) = rhtxy(imz,irec1(l)-1,1,1)
rxy(imz) = rhtxy(imz,irec1(l)-1,1)
END DO
!----- vacuum potential caused by the vacuum density
......@@ -510,10 +510,10 @@ CONTAINS
DO imz = 1,vacuum%nmzxy
pint(imz) = fpi_const*a*KK(imz)
vxy(imz,irec1(l)-1,1,1) = pint(imz) + pvac(imz)
vxy(imz,irec1(l)-1,1) = pint(imz) + pvac(imz)
END DO
val_help = vxy(1,irec1(l)-1,1,1) - val(irec1(l))
val_help = vxy(1,irec1(l)-1,1) - val(irec1(l))
CALL visp5_z(&
& vacuum%nmzxyd,vacuum%nmzxyd,vacuum%delz,m,ivfft1,ivfft2,IIIR,&
......@@ -571,10 +571,10 @@ CONTAINS
irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
IF (irec3.NE.0) THEN
IF (irec2.EQ.1) THEN
fxy0 = REAL(vpw(irec3,1))
rhti = AIMAG(vpw(irec3,1))
fxy0 = REAL(vpw(irec3))
rhti = AIMAG(vpw(irec3))
ELSE
fxy(irec2-1) = vpw(irec3,1)
fxy(irec2-1) = vpw(irec3)
END IF
END IF
END DO
......@@ -668,15 +668,15 @@ CONTAINS
END DO ! gz -> Vpw(.,.,gz)
DO irec3 = 1,stars%ng3
vpw(irec3,1) = vpw_help(irec3)
vpw(irec3) = vpw_help(irec3)
!$$$ vpw(irec3,1) = vpw(irec3,1) + vpw_help(irec3)
END DO
DO irc1 = 2,nq2_1