Commit 62689124 authored by Daniel Wortmann's avatar Daniel Wortmann

First attempt for bugfix for noco. Still not working correctly.

parent 510cc289
......@@ -229,7 +229,7 @@ CONTAINS
IF (noco%l_noco) THEN
CALL timestart("gen. spin-up and -down density")
CALL rhodirgen(DIMENSION,sym,stars,atoms,sphhar,&
vacuum,22,cell,input,noco,oneD)
vacuum,cell,input,noco,oneD,inDen)
CALL timestop("gen. spin-up and -down density")
ENDIF
!---> pk non-collinear
......@@ -307,7 +307,7 @@ CONTAINS
!---> and down potentials and direction of the magnetic field
IF (noco%l_noco) THEN
CALL timestart("generation of potential-matrix")
CALL vmatgen(stars, atoms,sphhar,vacuum,sym,input,oneD,8,22,26)
CALL vmatgen(stars, atoms,sphhar,vacuum,sym,input,oneD,inDen,vTot)
CALL timestop("generation of potential-matrix")
ENDIF
!---> end pk non-collinear
......
......@@ -659,8 +659,9 @@ CONTAINS
IF (input%total) THEN
IF (noco%l_noco) THEN ! load workDen%pw,workDen%vacz,workDen%vacxyxy from 'cdn'-file
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
0,fermiEnergyTemp,l_qfix,workDen)
PRINT *,"TODO: check that read is not needed in vgen"
! CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
! 0,fermiEnergyTemp,l_qfix,workDen)
ENDIF
!
! CALCULATE THE INTEGRAL OF n1*Veff1 + n2*Veff2
......
......@@ -9,4 +9,5 @@ types/types_lapw.F90
types/types_tlmplm.F90
types/types_misc.F90
types/types_mpimat.F90
types/types_potden.f90
)
......@@ -11,4 +11,5 @@ MODULE m_types
USE m_types_mpi
USE m_types_tlmplm
USE m_types_misc
USE m_types_potden
END MODULE m_types
......@@ -863,30 +863,7 @@ MODULE m_types_misc
END TYPE t_wann
TYPE t_potden
INTEGER :: iter
INTEGER :: potdenType
COMPLEX,ALLOCATABLE :: pw(:,:)
REAL,ALLOCATABLE :: mt(:,:,:,:)
REAL,ALLOCATABLE :: vacz(:,:,:)
COMPLEX,ALLOCATABLE :: vacxy(:,:,:,:)
! For density only (noco case)
COMPLEX, ALLOCATABLE :: cdom(:)
COMPLEX, ALLOCATABLE :: cdomvz(:,:)
COMPLEX, ALLOCATABLE :: cdomvxy(:,:,:)
! For density matrix and associated potential matrix
COMPLEX, ALLOCATABLE :: mmpMat(:,:,:,:)
!this type contains two init routines that should be used to allocate
!memory. You can either specify the datatypes or give the dimensions as integers
!See implementation below!
CONTAINS
PROCEDURE :: init_potden_types
PROCEDURE :: init_potden_simple
PROCEDURE :: resetpotden
GENERIC :: init=>init_potden_types,init_potden_simple
END TYPE t_potden
CONTAINS
SUBROUTINE usdus_init(ud,atoms,jsp)
USE m_judft
......@@ -910,86 +887,6 @@ CONTAINS
IF (ANY(err>0)) CALL judft_error("Not enough memory allocating usdus datatype")
END SUBROUTINE usdus_init
SUBROUTINE init_potden_types(pd,stars,atoms,sphhar,vacuum,noco,oneD,jsp,nocoExtraDim,potden_type)
USE m_judft
IMPLICIT NONE
CLASS(t_potden),INTENT(OUT):: pd
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_sphhar),INTENT(IN):: sphhar
TYPE(t_vacuum),INTENT(IN):: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_oneD),INTENT(IN) :: oneD
INTEGER,INTENT(IN) :: jsp, potden_type
LOGICAL,INTENT(IN) :: nocoExtraDim
CALL init_potden_simple(pd,stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,&
atoms%n_u,noco%l_noco,jsp,nocoExtraDim,potden_type,&
vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
END SUBROUTINE init_potden_types
SUBROUTINE init_potden_simple(pd,ng3,jmtd,nlhd,ntype,n_u,l_noco,jsp,nocoExtraDim,potden_type,nmzd,nmzxyd,n2d)
USE m_constants
USE m_judft
IMPLICIT NONE
CLASS(t_potden),INTENT(OUT) :: pd
INTEGER,INTENT(IN) :: ng3,jmtd,nlhd,ntype,n_u,jsp,potden_type
LOGICAL,INTENT(IN) :: l_noco,nocoExtraDim
INTEGER,INTENT(IN) :: nmzd,nmzxyd,n2d
INTEGER:: err(4)
err=0
pd%iter=0
pd%potdenType=potden_type
IF(ALLOCATED(pd%pw)) DEALLOCATE (pd%pw)
IF(ALLOCATED(pd%mt)) DEALLOCATE (pd%mt)
IF(ALLOCATED(pd%vacz)) DEALLOCATE (pd%vacz)
IF(ALLOCATED(pd%vacxy)) DEALLOCATE (pd%vacxy)
IF(ALLOCATED(pd%cdom)) DEALLOCATE (pd%cdom)
IF(ALLOCATED(pd%cdomvz)) DEALLOCATE (pd%cdomvz)
IF(ALLOCATED(pd%cdomvxy)) DEALLOCATE (pd%cdomvxy)
IF(ALLOCATED(pd%mmpMat)) DEALLOCATE (pd%mmpMat)
ALLOCATE (pd%pw(ng3,jsp),stat=err(1))
ALLOCATE (pd%mt(jmtd,0:nlhd,ntype,jsp),stat=err(2))
ALLOCATE (pd%vacz(nmzd,2,MERGE(4,jsp,nocoExtraDim)),stat=err(3))
ALLOCATE (pd%vacxy(nmzxyd,n2d-1,2,jsp),stat=err(4))
IF (l_noco) THEN
ALLOCATE (pd%cdom(ng3))
ALLOCATE (pd%cdomvz(nmzd,2))
ALLOCATE (pd%cdomvxy(nmzxyd,n2d-1,2))
ELSE
ALLOCATE (pd%cdom(1))
ALLOCATE (pd%cdomvz(1,1),pd%cdomvxy(1,1,1))
END IF
ALLOCATE (pd%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,MAX(1,n_u),jsp))
IF (ANY(err>0)) CALL judft_error("Not enough memory allocating potential or density")
pd%pw=CMPLX(0.0,0.0)
pd%mt=0.0
pd%vacz=0.0
pd%vacxy=CMPLX(0.0,0.0)
pd%cdom = CMPLX(0.0,0.0)
pd%cdomvz = CMPLX(0.0,0.0)
pd%cdomvxy = CMPLX(0.0,0.0)
pd%mmpMat = CMPLX(0.0,0.0)
END SUBROUTINE init_potden_simple
SUBROUTINE resetPotDen(pd)
IMPLICIT NONE
CLASS(t_potden),INTENT(INOUT) :: pd
pd%pw=CMPLX(0.0,0.0)
pd%mt=0.0
pd%vacz=0.0
pd%vacxy=CMPLX(0.0,0.0)
pd%cdom = CMPLX(0.0,0.0)
pd%cdomvz = CMPLX(0.0,0.0)
pd%cdomvxy = CMPLX(0.0,0.0)
pd%mmpMat = CMPLX(0.0,0.0)
END SUBROUTINE resetPotDen
END MODULE m_types_misc
......@@ -26,7 +26,7 @@ MODULE m_rhodirgen
!**********************************************************************
CONTAINS
SUBROUTINE rhodirgen(DIMENSION,sym,stars,atoms,sphhar,vacuum,&
ndomfile,cell,input,noco,oneD)
cell,input,noco,oneD,den)
!******** ABBREVIATIONS ***********************************************
! ifft3 : size of the 3d real space mesh
......@@ -42,8 +42,6 @@ CONTAINS
!**********************************************************************
USE m_constants
USE m_cdn_io
USE m_qfix
USE m_fft2d
USE m_fft3d
USE m_types
......@@ -59,74 +57,51 @@ CONTAINS
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: ndomfile
! ..
! ..
!-odim
!+odim
! Local type instances
TYPE(t_potden) :: den
TYPE(t_potden),INTENT(INOUT) :: den
! .. Local Scalars ..
INTEGER iden,jspin,ivac,ifft2,ifft3
INTEGER imz,ityp,iri,ilh,imesh,iq2,iq3
REAL zero,rho_11,rho_22,rho_21r,rho_21i,rhotot,magmom,phi
REAL rho_up,rho_down,mx,my,mz,eps,pi,fix,vz_r,vz_i,rziw,theta
REAL fermiEnergyTemp
COMPLEX czero
CHARACTER*8 dop,iop,name(10)
LOGICAL l_qfix
REAL rho_11,rho_22,rho_21r,rho_21i,rhotot,magmom,phi
REAL rho_up,rho_down,mx,my,mz,eps,vz_r,vz_i,rziw,theta
! ..
! .. Local Arrays ..
REAL, ALLOCATABLE :: rz(:,:,:)
REAL, ALLOCATABLE :: rvacxy(:,:,:,:),ris(:,:),fftwork(:)
! ..
zero = 0.0 ; czero = CMPLX(0.0,0.0)
eps = 1.0e-20
CALL den%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,.FALSE.,POTDEN_TYPE_DEN)
ALLOCATE (ris(0:27*stars%mx1*stars%mx2*stars%mx3-1,4),fftwork(0:27*stars%mx1*stars%mx2*stars%mx3-1),&
rz(vacuum%nmzd,2,2))
!
!---> initialize arrays for the density matrix
!
den%cdom(:) = czero
den%cdomvz(:,:) = czero
den%cdomvxy(:,:,:) = czero
IF (input%film) THEN
rz(:,:,:) = zero
ENDIF
ifft3 = 27*stars%mx1*stars%mx2*stars%mx3
ifft2 = 9*stars%mx1*stars%mx2
IF (oneD%odi%d1) ifft2 = 9*stars%mx3*oneD%odi%M
IF (input%film) ALLOCATE(rvacxy(0:ifft2-1,vacuum%nmzxyd,2,4))
IF (input%film) THEN
ifft2 = 9*stars%mx1*stars%mx2
IF (oneD%odi%d1) ifft2 = 9*stars%mx3*oneD%odi%M
ELSE
ifft2=0
END IF
IF (input%jspins .NE. 2) THEN
WRITE (6,*) 'This is the non-collinear version of the flapw-'
WRITE (6,*) 'program. It can only perform spin-polarized'
WRITE (6,*) 'calculations.'
CALL juDFT_error("jspins not equal 2",calledby="rhodirgen")
IF (ALLOCATED(den%phi_pw)) THEN
DEALLOCATE(den%phi_pw,den%phi_vacz,den%phi_vacxy)
DEALLOCATE(den%theta_pw,den%theta_vacz,den%theta_vacxy)
ENDIF
ALLOCATE(den%phi_pw(ifft3),den%theta_pw(ifft3))
ALLOCATE(den%phi_vacz(vacuum%nmzd,2),den%theta_vacz(vacuum%nmzd,2))
ALLOCATE(den%phi_vacxy(ifft2,vacuum%nmzxyd,2),den%theta_vacxy(ifft2,vacuum%nmzxyd,2))
!---> reload the density matrix from file rhomat_inp
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_NOCO_const,CDN_INPUT_DEN_const,&
0,fermiEnergyTemp,l_qfix,den)
CALL qfix(stars,atoms,sym,vacuum,sphhar,input,cell,oneD,&
den%pw,den%vacxy,den%mt,den%vacz,.FALSE.,.FALSE.,fix)
ALLOCATE (ris(ifft3,4),fftwork(ifft3))
!---> fouriertransform the diagonal part of the density matrix
!---> in the interstitial, den%pw, to real space (ris)
DO iden = 1,2
CALL fft3d(ris(0,iden),fftwork,den%pw(1,iden),stars,+1)
CALL fft3d(ris(:,iden),fftwork,den%pw(:,iden),stars,+1)
ENDDO
!---> fouriertransform the off-diagonal part of the density matrix
CALL fft3d(ris(0,3),ris(0,4),den%cdom(1),stars,+1)
CALL fft3d(ris(:,3),ris(:,4),den%cdom(:),stars,+1)
!test
! DO iden=1,4
......@@ -136,7 +111,7 @@ CONTAINS
!test
!---> calculate the charge and magnetization density on the
!---> real space mesh
DO imesh = 0,ifft3-1
DO imesh = 1,ifft3
rho_11 = ris(imesh,1)
rho_22 = ris(imesh,2)
rho_21r = ris(imesh,3)
......@@ -178,117 +153,72 @@ CONTAINS
! write(36,'(i4,2f12.6)') mod(imesh,33),rho_11,rho_22
ris(imesh,1) = rho_up
ris(imesh,2) = rho_down
ris(imesh,3) = theta
ris(imesh,4) = phi
den%theta_pw(imesh) = theta
den%phi_pw(imesh) = phi
ENDDO
!test
! DO iden=1,4
! write(*,*)'iden=',iden
! write(*,8500)(ris(imesh,iden),imesh=0,ifft3-1)
! 8500 format(10e13.5)
! enddo
!test
!---> Fouriertransform the density matrix back to reciprocal space
DO jspin = 1,input%jspins
fftwork=0.0
CALL fft3d(ris(0,jspin),fftwork,den%pw(1,jspin),stars,-1)
CALL fft3d(ris(:,jspin),fftwork,den%pw(:,jspin),stars,-1)
ENDDO
IF (.NOT.input%film) RETURN
!Now the vacuum part starts
ALLOCATE(rvacxy(ifft2,vacuum%nmzxyd,2,4))
ALLOCATE (rz(vacuum%nmzd,2,2))
!---> fouriertransform the diagonal part of the density matrix
!---> in the vacuum, rz & rxy, to real space (rvacxy)
IF (input%film) THEN
DO iden = 1,2
DO ivac = 1,vacuum%nvac
DO imz = 1,vacuum%nmzxyd
rziw = 0.0
IF (oneD%odi%d1) THEN
call judft_error("oneD not implemented",calledby="rhodirgen")
!CALL fft2d(oneD%k3,odi%M,odi%n2d,rvacxy(0,imz,ivac,iden),fftwork,&
! den%vacz(imz,ivac,iden),rziw,den%vacxy(imz,1,ivac,iden),&
! vacuum,odi%nq2,odi%kimax2,1,&
! & %igf,odl%pgf,odi%nst2)
ELSE
CALL fft2d(stars,rvacxy(0,imz,ivac,iden),fftwork,&
den%vacz(imz,ivac,iden),rziw,den%vacxy(imz,1,ivac,iden),&
vacuum%nmzxyd,1)
ENDIF
ENDDO
ENDDO
ENDDO
!---> fouriertransform the off-diagonal part of the density matrix
DO iden = 1,2
DO ivac = 1,vacuum%nvac
DO imz = 1,vacuum%nmzxyd
rziw = 0.0
vz_r = REAL(den%cdomvz(imz,ivac))
vz_i = AIMAG(den%cdomvz(imz,ivac))
IF (oneD%odi%d1) THEN
call judft_error("oneD not implemented",calledby="rhodirgen")
!CALL fft2d(oneD%k3,odi%M,odi%n2d,&
! rvacxy(0,imz,ivac,3),rvacxy(0,imz,ivac,4),&
! vz_r,vz_i,den%cdomvxy(imz,1,ivac),&
CALL judft_error("oneD not implemented",calledby="rhodirgen")
!CALL fft2d(oneD%k3,odi%M,odi%n2d,rvacxy(0,imz,ivac,iden),fftwork,&
! den%vacz(imz,ivac,iden),rziw,den%vacxy(imz,1,ivac,iden),&
! vacuum,odi%nq2,odi%kimax2,1,&
! & %igf,odl%pgf,odi%nst2)
! & %igf,odl%pgf,odi%nst2)
ELSE
CALL fft2d(stars,rvacxy(0,imz,ivac,3),rvacxy(0,imz,ivac,4),&
vz_r,vz_i,den%cdomvxy(imz,1,ivac),vacuum%nmzxyd,1)
CALL fft2d(stars,rvacxy(:,imz,ivac,iden),fftwork,&
den%vacz(imz,ivac,iden),rziw,den%vacxy(imz,1,ivac,iden),&
vacuum%nmzxyd,1)
ENDIF
ENDDO
ENDDO
ENDDO
!---> fouriertransform the off-diagonal part of the density matrix
DO ivac = 1,vacuum%nvac
DO imz = 1,vacuum%nmzxyd
rziw = 0.0
vz_r = REAL(den%cdomvz(imz,ivac))
vz_i = AIMAG(den%cdomvz(imz,ivac))
IF (oneD%odi%d1) THEN
CALL judft_error("oneD not implemented",calledby="rhodirgen")
!CALL fft2d(oneD%k3,odi%M,odi%n2d,&
! rvacxy(0,imz,ivac,3),rvacxy(0,imz,ivac,4),&
! vz_r,vz_i,den%cdomvxy(imz,1,ivac),&
! vacuum,odi%nq2,odi%kimax2,1,&
! & %igf,odl%pgf,odi%nst2)
ELSE
CALL fft2d(stars,rvacxy(:,imz,ivac,3),rvacxy(:,imz,ivac,4),&
vz_r,vz_i,den%cdomvxy(imz,1,ivac),vacuum%nmzxyd,1)
ENDIF
ENDDO
ENDDO
!---> calculate the four components of the matrix potential on
!---> real space mesh
DO ivac = 1,vacuum%nvac
DO imz = 1,vacuum%nmzxyd
DO imesh = 0,ifft2-1
rho_11 = rvacxy(imesh,imz,ivac,1)
rho_22 = rvacxy(imesh,imz,ivac,2)
rho_21r = rvacxy(imesh,imz,ivac,3)
rho_21i = rvacxy(imesh,imz,ivac,4)
mx = 2*rho_21r
my = -2*rho_21i
mz = (rho_11-rho_22)
magmom = SQRT(mx**2 + my**2 + mz**2)
rhotot = rho_11 + rho_22
rho_up = (rhotot + magmom)/2
rho_down= (rhotot - magmom)/2
IF (ABS(mz) .LE. eps) THEN
theta = pi_const/2
ELSEIF (mz .GE. 0.0) THEN
theta = ATAN(SQRT(mx**2 + my**2)/mz)
ELSE
theta = ATAN(SQRT(mx**2 + my**2)/mz) + pi_const
ENDIF
IF (ABS(mx) .LE. eps) THEN
IF (ABS(my) .LE. eps) THEN
phi = 0.0
ELSEIF (my .GE. 0.0) THEN
phi = pi_const/2
ELSE
phi = -pi_const/2
ENDIF
ELSEIF (mx .GE. 0.0) THEN
phi = ATAN(my/mx)
ELSE
IF (my .GE. 0.0) THEN
phi = ATAN(my/mx) + pi_const
ELSE
phi = ATAN(my/mx) - pi_const
ENDIF
ENDIF
rvacxy(imesh,imz,ivac,1) = rho_up
rvacxy(imesh,imz,ivac,2) = rho_down
rvacxy(imesh,imz,ivac,3) = theta
rvacxy(imesh,imz,ivac,4) = phi
ENDDO
ENDDO
DO imz = vacuum%nmzxyd+1,vacuum%nmzd
rho_11 = den%vacz(imz,ivac,1)
rho_22 = den%vacz(imz,ivac,2)
rho_21r = REAL(den%cdomvz(imz,ivac))
rho_21i = AIMAG(den%cdomvz(imz,ivac))
!---> calculate the four components of the matrix potential on
!---> real space mesh
DO ivac = 1,vacuum%nvac
DO imz = 1,vacuum%nmzxyd
DO imesh = 1,ifft2
rho_11 = rvacxy(imesh,imz,ivac,1)
rho_22 = rvacxy(imesh,imz,ivac,2)
rho_21r = rvacxy(imesh,imz,ivac,3)
rho_21i = rvacxy(imesh,imz,ivac,4)
mx = 2*rho_21r
my = -2*rho_21i
mz = (rho_11-rho_22)
......@@ -323,52 +253,78 @@ CONTAINS
ENDIF
ENDIF
den%vacz(imz,ivac,1) = rho_up
den%vacz(imz,ivac,2) = rho_down
rz(imz,ivac,1) = theta
rz(imz,ivac,2) = phi
rvacxy(imesh,imz,ivac,1) = rho_up
rvacxy(imesh,imz,ivac,2) = rho_down
den%theta_vacxy(imesh,imz,ivac) = theta
den%phi_vacxy(imesh,imz,ivac) = phi
ENDDO
ENDDO
!---> Fouriertransform the matrix potential back to reciprocal space
DO jspin = 1,input%jspins
DO ivac = 1,vacuum%nvac
DO imz = 1,vacuum%nmzxyd
fftwork=0.0
IF (oneD%odi%d1) THEN
call judft_error("oneD not implemented",calledby="rhodirgen")
!CALL fft2d(oneD%k3,odi%M,odi%n2d,&
! rvacxy(0,imz,ivac,jspin),fftwork,&
! den%vacz(imz,ivac,jspin),rziw,den%vacxy(imz,1,ivac,jspin),&
! vacuum,odi%nq2,odi%kimax2,-1,&
! & %igf,odl%pgf,odi%nst2)
ELSE
CALL fft2d(stars,rvacxy(0,imz,ivac,jspin),fftwork,&
den%vacz(imz,ivac,jspin),rziw,den%vacxy(imz,1,ivac,jspin),&
vacuum%nmzxyd,-1)
END IF
ENDDO
ENDDO
ENDDO
ENDIF
!---> ndomfile is the dirofmag-file
OPEN (ndomfile,FILE='dirofmag',FORM='unformatted',STATUS='unknown')
WRITE (ndomfile) (ris(imesh,3),imesh=0,ifft3-1)
WRITE (ndomfile) (ris(imesh,4),imesh=0,ifft3-1)
IF (input%film) THEN
WRITE (ndomfile) ((rz(imz,ivac,1),imz=vacuum%nmzxyd+1,vacuum%nmzd),ivac=1,vacuum%nvac)
WRITE (ndomfile) ((rz(imz,ivac,2),imz=vacuum%nmzxyd+1,vacuum%nmzd),ivac=1,vacuum%nvac)
WRITE (ndomfile) (((rvacxy(imesh,imz,ivac,3),imesh=0,ifft2-1),imz=1,vacuum%nmzxyd),ivac=1,vacuum%nvac)
WRITE (ndomfile) (((rvacxy(imesh,imz,ivac,4),imesh=0,ifft2-1),imz=1,vacuum%nmzxyd),ivac=1,vacuum%nvac)
ENDIF
CLOSE (ndomfile)
DO imz = vacuum%nmzxyd+1,vacuum%nmzd
rho_11 = den%vacz(imz,ivac,1)
rho_22 = den%vacz(imz,ivac,2)
rho_21r = REAL(den%cdomvz(imz,ivac))
rho_21i = AIMAG(den%cdomvz(imz,ivac))
mx = 2*rho_21r
my = -2*rho_21i
mz = (rho_11-rho_22)
magmom = SQRT(mx**2 + my**2 + mz**2)
rhotot = rho_11 + rho_22
rho_up = (rhotot + magmom)/2
rho_down= (rhotot - magmom)/2
!---> write spin-up and -down density on file cdn
IF (ABS(mz) .LE. eps) THEN
theta = pi_const/2
ELSEIF (mz .GE. 0.0) THEN
theta = ATAN(SQRT(mx**2 + my**2)/mz)
ELSE
theta = ATAN(SQRT(mx**2 + my**2)/mz) + pi_const
ENDIF
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
0,-1.0,0.0,.FALSE.,den)
IF (ABS(mx) .LE. eps) THEN
IF (ABS(my) .LE. eps) THEN
phi = 0.0
ELSEIF (my .GE. 0.0) THEN
phi = pi_const/2
ELSE
phi = -pi_const/2
ENDIF
ELSEIF (mx .GE. 0.0) THEN
phi = ATAN(my/mx)
ELSE
IF (my .GE. 0.0) THEN
phi = ATAN(my/mx) + pi_const
ELSE
phi = ATAN(my/mx) - pi_const
ENDIF
ENDIF
DEALLOCATE (ris,fftwork,rz)
IF (input%film) DEALLOCATE(rvacxy)
den%vacz(imz,ivac,1) = rho_up
den%vacz(imz,ivac,2) = rho_down
den%theta_vacz(imz,ivac) = theta
den%phi_vacz(imz,ivac) = phi
ENDDO
ENDDO
!---> Fouriertransform the matrix potential back to reciprocal space
DO jspin = 1,input%jspins
DO ivac = 1,vacuum%nvac
DO imz = 1,vacuum%nmzxyd
fftwork=0.0
IF (oneD%odi%d1) THEN
call judft_error("oneD not implemented",calledby="rhodirgen")
!CALL fft2d(oneD%k3,odi%M,odi%n2d,&
! rvacxy(0,imz,ivac,jspin),fftwork,&
! den%vacz(imz,ivac,jspin),rziw,den%vacxy(imz,1,ivac,jspin),&
! vacuum,odi%nq2,odi%kimax2,-1,&
! & %igf,odl%pgf,odi%nst2)
ELSE
CALL fft2d(stars,rvacxy(:,imz,ivac,jspin),fftwork,&
den%vacz(imz,ivac,jspin),rziw,den%vacxy(imz,1,ivac,jspin),&
vacuum%nmzxyd,-1)
END IF
ENDDO
ENDDO
ENDDO
RETURN
END SUBROUTINE rhodirgen
END MODULE m_rhodirgen
......@@ -23,11 +23,7 @@ MODULE m_vmatgen
! Philipp Kurz 99/11/01
!**********************************************************************
CONTAINS
SUBROUTINE vmatgen(&
& stars,&
& atoms,sphhar,vacuum,&
& sym,input,oneD,&
& nu,ndomfile,npotmatfile)
SUBROUTINE vmatgen(stars, atoms,sphhar,vacuum, sym,input,oneD, den,v)
!******** ABBREVIATIONS ***********************************************
! ifft3 : size of the 3d real space mesh
......@@ -41,7 +37,6 @@ CONTAINS
! all stored on real space mesh
!**********************************************************************
USE m_loddop
USE m_fft2d
USE m_fft3d
USE m_types
......@@ -53,102 +48,39 @@ CONTAINS
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden),INTENT(IN) :: den
TYPE(t_potden),INTENT(INOUT):: v
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: nu,ndomfile,npotmatfile
! ..
! .. Local Scalars ..