Commit fa58d844 authored by Gregor Michalicek's avatar Gregor Michalicek

Eliminate potden%cdomvxy

Note: This seems to have introduced some bugs. ...to be corrected soon.
parent f92042da
......@@ -2,8 +2,8 @@ MODULE m_cdnval
use m_juDFT
CONTAINS
SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,sliceplot,noco, input,banddos,cell,atoms,enpara,stars,&
vacuum,dimension,sphhar,sym,obsolete,igq_fft,vr,vz,oneD,coreSpecInput,n_mmp,results, qpw,rhtxy,&
rho,rht,cdom,cdomvz,cdomvxy,qvac,qvlay,qa21, chmom,clmom)
vacuum,dimension,sphhar,sym,obsolete,igq_fft,vr,vz,oneD,coreSpecInput,den,results,&
qvac,qvlay,qa21, chmom,clmom)
!
! ***********************************************************
! this subroutin is a modified version of cdnval.F.
......@@ -44,6 +44,7 @@ CONTAINS
! and bands
!***********************************************************************
!
USE m_constants
USE m_eig66_io,ONLY: write_dos
USE m_radfun
USE m_radflo
......@@ -111,26 +112,19 @@ CONTAINS
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_coreSpecInput),INTENT(IN) :: coreSpecInput
TYPE(t_potden),INTENT(INOUT) :: den
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: eig_id,jspin
! .. Array Arguments ..
COMPLEX, INTENT(INOUT) :: qpw(stars%ng3,dimension%jspd)
COMPLEX, INTENT(INOUT) :: rhtxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,dimension%jspd)
COMPLEX, INTENT(INOUT) :: cdom(stars%ng3)
COMPLEX, INTENT(INOUT) :: cdomvz(vacuum%nmzd,2)
COMPLEX, INTENT(INOUT) :: cdomvxy(vacuum%nmzxyd,oneD%odi%n2d-1,2)
COMPLEX, INTENT(INOUT) :: qa21(atoms%ntype)
INTEGER, INTENT (IN) :: igq_fft(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1)
REAL, INTENT (IN) :: vz(vacuum%nmzd,2)
REAL, INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd)
REAL, INTENT (OUT) :: chmom(atoms%ntype,dimension%jspd),clmom(3,atoms%ntype,dimension%jspd)
REAL, INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd)
REAL, INTENT (INOUT) :: rht(vacuum%nmzd,2,dimension%jspd)
REAL, INTENT (INOUT) :: qvac(dimension%neigd,2,kpts%nkpt,dimension%jspd)
REAL, INTENT (INOUT) :: qvlay(dimension%neigd,vacuum%layerd,2,kpts%nkpt,dimension%jspd)
COMPLEX, INTENT(INOUT) :: n_mmp(-3:3,-3:3,atoms%n_u)
#ifdef CPP_MPI
INCLUDE 'mpif.h'
......@@ -665,7 +659,7 @@ CONTAINS
IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
CALL timestart("cdnval: pwden")
CALL pwden(stars,kpts,banddos,oneD, input,mpi,noco,cell,atoms,sym,ikpt,&
jspin,lapw,noccbd,igq_fft,we, eig,qpw,cdom,qis,results%force,f_b8,zMat)
jspin,lapw,noccbd,igq_fft,we, eig,den,qis,results%force,f_b8,zMat)
CALL timestop("cdnval: pwden")
END IF
!+new
......@@ -687,7 +681,7 @@ CONTAINS
CALL timestart("cdnval: vacden")
CALL vacden(vacuum,dimension,stars,oneD, kpts,input, cell,atoms,noco,banddos,&
gvac1d,gvac2d, we,ikpt,jspin,vz,vz0, noccbd,lapw, evac,eig,&
rhtxy,rht,qvac,qvlay, qstars,cdomvz,cdomvxy,zMat)
den,qvac,qvlay, qstars,zMat)
CALL timestop("cdnval: vacden")
END IF
!---> perform Brillouin zone integration and summation over the
......@@ -742,7 +736,7 @@ CONTAINS
IF (atoms%n_u.GT.0) THEN
CALL n_mat(atoms,sym,noccbd,usdus,ispin,we, acof(:,0:,:,ispin),bcof(:,0:,:,ispin),&
ccof(-atoms%llod:,:,:,:,ispin), n_mmp)
ccof(-atoms%llod:,:,:,:,ispin), den%mmpMat(:,:,:,jspin))
END IF
!
!---> perform Brillouin zone integration and summation over the
......@@ -904,8 +898,8 @@ CONTAINS
CALL timestart("cdnval: mpi_col_den")
DO ispin = jsp_start,jsp_end
CALL mpi_col_den(mpi,sphhar,atoms,oneD,stars,vacuum,&
input,noco,l_fmpl,ispin,llpd, rhtxy(1,1,1,ispin),&
rht(1,1,ispin),qpw(1,ispin), ener(0,1,ispin),sqal(0,1,ispin),&
input,noco,l_fmpl,ispin,llpd, den%vacxy(1,1,1,ispin),&
den%vacz(1,1,ispin),den%pw(1,ispin), ener(0,1,ispin),sqal(0,1,ispin),&
results,svac(1,ispin),pvac(1,ispin),uu(0,1,ispin),&
dd(0,1,ispin),du(0,1,ispin),uunmt(0,1,1,ispin),ddnmt(0,1,1,ispin),&
udnmt(0,1,1,ispin),dunmt(0,1,1,ispin),sqlo(1,1,ispin),&
......@@ -914,7 +908,7 @@ CONTAINS
ccnmt(1,1,1,1,ispin),enerlo(1,1,ispin),&
orb(0,-atoms%lmaxd,1,ispin),orbl(1,-atoms%llod,1,ispin),&
orblo(1,1,-atoms%llod,1,ispin),mt21,lo21,uloulop21,&
uunmt21,ddnmt21,udnmt21,dunmt21,cdom,cdomvz,cdomvxy,n_mmp)
uunmt21,ddnmt21,udnmt21,dunmt21,den,den%mmpMat(:,:,:,jspin))
END DO
CALL timestop("cdnval: mpi_col_den")
#endif
......@@ -954,7 +948,7 @@ CONTAINS
orb,orbl,orblo,mt21,lo21,uloulopn21,uloulop21,&
uunmt21,ddnmt21,udnmt21,dunmt21,&
chmom,clmom,&
qa21,rho)
qa21,den%mt)
DO ispin = jsp_start,jsp_end
WRITE (6,*) 'Energy Parameters for spin:',ispin
......@@ -992,7 +986,7 @@ CONTAINS
xp,npd,0,0,ivac,1,ispin,.true.,dimension,atoms,&
sphhar,stars,sym,&
vacuum,cell,oneD,&
qpw,rho,rhtxy,rht)
den%pw,den%mt,den%vacxy,den%vacz)
END DO
ELSE IF (oneD%odi%d1) THEN
!-odim
......@@ -1002,7 +996,7 @@ CONTAINS
xp,npd,0,0,ivac,1,ispin,.true.,dimension,atoms,&
sphhar,stars,sym,&
vacuum,cell,oneD,&
qpw,rho,rhtxy,rht)
den%pw,den%mt,den%vacxy,den%vacz)
!+odim
END IF
!---> m.t. boundaries
......@@ -1013,7 +1007,7 @@ CONTAINS
xp,dimension%nspd,n,nat,0,-1,ispin,.true.,&
dimension,atoms,sphhar,stars,sym,&
vacuum,cell,oneD,&
qpw,rho,rhtxy,rht)
den%pw,den%mt,den%vacxy,den%vacz)
nat = nat + atoms%neq(n)
END DO
CALL timestop("cdnval: cdninf-stuff")
......@@ -1023,7 +1017,7 @@ CONTAINS
!---> forces of equ. A8 of Yu et al.
IF ((input%l_f)) THEN
CALL timestart("cdnval: force_a8")
CALL force_a8(input,atoms,sphhar, ispin, vr(:,:,:,ispin),rho,&
CALL force_a8(input,atoms,sphhar, ispin, vr(:,:,:,ispin),den%mt,&
f_a12,f_a21,f_b4,f_b8,results%force)
CALL timestop("cdnval: force_a8")
END IF
......
......@@ -18,6 +18,7 @@ CONTAINS
!
USE m_types
USE m_constants
IMPLICIT NONE
TYPE(t_usdus),INTENT(IN) :: usdus
TYPE(t_sym),INTENT(IN) :: sym
......@@ -31,7 +32,7 @@ CONTAINS
COMPLEX, INTENT (IN) :: acof(:,0:,:)!(nobd,0:atoms%lmaxd*(lmaxd+2) ,natd)
COMPLEX, INTENT (IN) :: bcof(:,0:,:)!(nobd,0:atoms%lmaxd*(lmaxd+2) ,natd)
COMPLEX, INTENT (IN) :: ccof(-atoms%llod:,:,:,:)!(-llod:llod,nobd,atoms%nlod,atoms%nat)
COMPLEX, INTENT (INOUT) :: n_mmp(-3:3,-3:3,atoms%n_u)
COMPLEX, INTENT (INOUT) :: n_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u)
! ..
! .. Local Scalars ..
COMPLEX c_0
......
......@@ -7,7 +7,7 @@
MODULE m_pwden
CONTAINS
SUBROUTINE pwden(stars,kpts,banddos,oneD, input,mpi,noco,cell,atoms,sym, &
ikpt,jspin,lapw,ne, igq_fft,we,eig, qpw,cdom, qis,forces,f_b8,zMat)
ikpt,jspin,lapw,ne, igq_fft,we,eig,den,qis,forces,f_b8,zMat)
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
! In this subroutine the star function expansion coefficients of
! the plane wave charge density is determined.
......@@ -31,24 +31,24 @@ CONTAINS
! Brillouine zone sampling
! FFT information
!
! OUTPUT: qpw(s)
! OUTPUT: den%pw(s)
! 1) using FFT
!
! 2) traditional method
!
! -1 ef
! qpw (g) = vol * sum{ sum{ sum{ sum{ w(k) * f(nu) *
! den%pw(g) = vol * sum{ sum{ sum{ sum{ w(k) * f(nu) *
! sp k nu g'
! *
! c(g'-g,nu,k) * c(g',nu,k) } } } }
! or :
! -1 ef
! qpw (g) = vol * sum{ sum{ sum{ sum{ w(k) * f(nu) *
! den%pw(g) = vol * sum{ sum{ sum{ sum{ w(k) * f(nu) *
! sp k nu g'
! *
! c(g',nu,k) * c(g'+g,nu,k) } } } }
!
! qpw(g) are actuall
! den%pw(g) are actuall
!
! the weights w(k) are normalized: sum{w(k)} = 1
! k -6
......@@ -63,8 +63,8 @@ CONTAINS
! In non-collinear calculations the density becomes a hermitian 2x2
! matrix. This subroutine generates this density matrix in the
! interstitial region. The diagonal elements of this matrix
! (n_11 & n_22) are stored in qpw, while the real and imaginary part
! of the off-diagonal element are store in cdom.
! (n_11 & n_22) are stored in den%pw, while the real and imaginary part
! of the off-diagonal element are store in den%cdom.
!
! Philipp Kurz 99/07
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
......@@ -80,18 +80,19 @@ CONTAINS
USE m_types
USE m_fft_interface
IMPLICIT NONE
TYPE(t_lapw),INTENT(IN) :: lapw
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_banddos),INTENT(IN) :: banddos
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_zMat),INTENT(IN) :: zMat
TYPE(t_lapw),INTENT(IN) :: lapw
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_banddos),INTENT(IN) :: banddos
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_zMat),INTENT(IN) :: zMat
TYPE(t_potden),INTENT(INOUT) :: den
INTEGER, INTENT (IN) :: igq_fft(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1)
REAL,INTENT(IN) :: we(:) !(nobd)
......@@ -100,8 +101,6 @@ CONTAINS
INTEGER,INTENT(IN):: ne
!-----> CHARGE DENSITY INFORMATION
INTEGER,INTENT(IN) :: ikpt,jspin
COMPLEX,INTENT(INOUT) :: qpw(:,:) !(stars%ng3,dimension%jspd)
COMPLEX,INTENT(INOUT) :: cdom(:)!(stars%ng3)
REAL,INTENT(OUT) :: qis(:,:,:) !(dimension%neigd,kpts%nkpt,dimension%jspd)
COMPLEX, INTENT (INOUT) :: f_b8(3,atoms%ntype)
REAL, INTENT (INOUT) :: forces(:,:,:) !(3,atoms%ntype,dimension%jspd)
......@@ -148,7 +147,7 @@ CONTAINS
! cv=z : wavefunction in g-space (reciprocal space)
! psir : wavefunction in r-space (real-space)
! cwk : complex work array: charge density in g-space (as stars)
! qpw : charge density stored as stars
! den%pw : charge density stored as stars
! trdchg: logical key, determines the mode of charge density
! calculation: false (default) : fft
! true : double sum over stars
......@@ -686,17 +685,17 @@ CONTAINS
ispin = jspin
IF (noco%l_noco) ispin = idens
DO istr = 1 , stars%ng3_fft
qpw(istr,ispin) = qpw(istr,ispin) + cwk(istr)
den%pw(istr,ispin) = den%pw(istr,ispin) + cwk(istr)
ENDDO
ELSE IF (idens.EQ.3) THEN
!---> add to off-diag. part of density matrix (only non-collinear)
DO istr = 1 , stars%ng3_fft
cdom(istr) = cdom(istr) + cwk(istr)
den%cdom(istr) = den%cdom(istr) + cwk(istr)
ENDDO
ELSE
!---> add to off-diag. part of density matrix (only non-collinear)
DO istr = 1 , stars%ng3_fft
cdom(istr) = cdom(istr) + CMPLX(0.0,1.0)*cwk(istr)
den%cdom(istr) = den%cdom(istr) + CMPLX(0.0,1.0)*cwk(istr)
ENDDO
ENDIF
......
......@@ -11,8 +11,8 @@ CONTAINS
gvac1,gvac2,&
we,ikpt,jspin,vz,vz0,&
ne,lapw,&
evac,eig,rhtxy,rht,qvac,qvlay,&
stcoeff,cdomvz,cdomvxy,zMat)
evac,eig,den,qvac,qvlay,&
stcoeff,zMat)
!***********************************************************************
! ****** change vacden(....,q) for vacuum density of states shz Jan.96
......@@ -21,8 +21,8 @@ CONTAINS
! In non-collinear calculations the density becomes a hermitian 2x2
! matrix. This subroutine generates this density matrix in the
! vacuum region. The diagonal elements of this matrix (n_11 & n_22)
! are store in rht and rhtxy, while the real and imaginary part
! of the off-diagonal element are store in cdomvz and cdomvxy.
! are store in den%vacz and den%vacxy, while the real and imaginary part
! of the off-diagonal element are store in den%cdomvz and den%vacxy(:,:,:,3).
!
! Philipp Kurz 99/07
!***********************************************************************
......@@ -33,13 +33,13 @@ CONTAINS
! vz : non-warping part of the vacuum potential (matrix)
! collinear : 2. index = ivac (# of vaccum)
! non-collinear: 2. index = ipot (comp. of pot. matr.)
! rht : non-warping part of the vacuum density matrix,
! den%vacz : non-warping part of the vacuum density matrix,
! diagonal elements n_11 and n_22
! rhtxy : warping part of the vacuum density matrix,
! den%vacxy: warping part of the vacuum density matrix,
! diagonal elements n_11 and n_22
! cdomvz : non-warping part of the vacuum density matrix,
! den%cdomvz: non-warping part of the vacuum density matrix,
! off-diagonal elements n_21
! cdomvxy : warping part of the vacuum density matrix,
! den%vacxy(:,:,:,3): warping part of the vacuum density matrix,
! off-diagonal elements n_21
!***********************************************************************
!
......@@ -65,6 +65,7 @@ CONTAINS
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_zMat),INTENT(IN) :: zMat
TYPE(t_potden),INTENT(INOUT) :: den
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jspin
INTEGER, INTENT (IN) :: ne
......@@ -73,12 +74,8 @@ CONTAINS
REAL,PARAMETER :: emax=2.0/hartree_to_ev_const
! .. Array Arguments ..
REAL, INTENT (IN) :: evac(2,DIMENSION%jspd)
COMPLEX, INTENT (INOUT):: rhtxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,DIMENSION%jspd)
REAL, INTENT (INOUT):: rht(vacuum%nmzd,2,DIMENSION%jspd)
REAL, INTENT (OUT) :: qvlay(DIMENSION%neigd,vacuum%layerd,2,kpts%nkpt,DIMENSION%jspd)
REAL qvac(DIMENSION%neigd,2,kpts%nkpt,DIMENSION%jspd),we(DIMENSION%neigd),vz(vacuum%nmzd,2),vz0(2)
COMPLEX, INTENT (INOUT):: cdomvz(vacuum%nmzd,2)
COMPLEX, INTENT (INOUT):: cdomvxy(vacuum%nmzxyd,oneD%odi%n2d-1,2)
!
! STM-Arguments
REAL, INTENT (IN) :: eig(DIMENSION%neigd)
......@@ -574,7 +571,7 @@ CONTAINS
uei = ue(jz,l,jspin)
ddui = ddu(jz,l)
dduei = ddue(jz,l)
rht(jz,ivac,jspin) = rht(jz,ivac,jspin) +&
den%vacz(jz,ivac,jspin) = den%vacz(jz,ivac,jspin) +&
REAL(aaee*ui*ui+bbee*uei*uei+&
(abee+baee)*ui*uei+aa*ddui*ddui+&
bb*dduei*dduei+(ab+ba)*ddui*dduei+&
......@@ -612,7 +609,7 @@ CONTAINS
DO jz = 1,vacuum%nmz
ui = u_1(jz,l,m,ispin)
uei = ue_1(jz,l,m,ispin)
rht(jz,ivac,ispin) = rht(jz,ivac,ispin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
den%vacz(jz,ivac,ispin) = den%vacz(jz,ivac,ispin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
END DO
END DO
END DO
......@@ -633,7 +630,7 @@ CONTAINS
DO jz = 1,vacuum%nmz
ui = u(jz,l,ispin)
uei = ue(jz,l,ispin)
rht(jz,ivac,ispin) = rht(jz,ivac,ispin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
den%vacz(jz,ivac,ispin) = den%vacz(jz,ivac,ispin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
ENDDO
ENDDO
END IF ! one-dimensional
......@@ -658,7 +655,7 @@ CONTAINS
DO jz = 1,vacuum%nmz
ui = u_1(jz,l,m,jspin)
uei = ue_1(jz,l,m,jspin)
rht(jz,ivac,jspin) = rht(jz,ivac,jspin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
den%vacz(jz,ivac,jspin) = den%vacz(jz,ivac,jspin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
END DO
END DO
END DO
......@@ -679,7 +676,7 @@ CONTAINS
DO jz = 1,vacuum%nmz
ui = u(jz,l,jspin)
uei = ue(jz,l,jspin)
rht(jz,ivac,jspin) = rht(jz,ivac,jspin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
den%vacz(jz,ivac,jspin) = den%vacz(jz,ivac,jspin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
ENDDO
END DO
END IF ! D1
......@@ -906,8 +903,8 @@ CONTAINS
+ abe*(ui*dduej+uj*dduei)+bae*(ddui*uej+dduj*uei)&
+ aa*ddui*dduj+bb*dduei*dduej+ba*ddui*dduej&
+ ab*dduei*dduj
rhtxy(jz,ind2-1,ivac,jspin) = rhtxy(jz,ind2-1,ivac,jspin) + t1*phs/stars%nstr2(ind2)
rhtxy(jz,ind2p-1,ivac,jspin)=rhtxy(jz,ind2p-1,ivac,jspin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
den%vacxy(jz,ind2-1,ivac,jspin) = den%vacxy(jz,ind2-1,ivac,jspin) + t1*phs/stars%nstr2(ind2)
den%vacxy(jz,ind2p-1,ivac,jspin)= den%vacxy(jz,ind2p-1,ivac,jspin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
ENDDO
ENDDO
END DO
......@@ -951,10 +948,10 @@ CONTAINS
uej = ue_1(jz,l1,m1,ispin)
t1 = aa*ui*uj + bb*uei*uej + ba*ui*uej + ab*uei*uj
IF (ind1.NE.0) THEN
rhtxy(jz,ind1-1,ivac,ispin) = rhtxy(jz,ind1-1,ivac,ispin) + t1/ oneD%odi%nst2(ind1)
den%vacxy(jz,ind1-1,ivac,ispin) = den%vacxy(jz,ind1-1,ivac,ispin) + t1/ oneD%odi%nst2(ind1)
END IF
IF (ind1p.NE.0) THEN
rhtxy(jz,ind1p-1,ivac,ispin) = rhtxy(jz,ind1p-1,ivac,ispin) + CONJG(t1)/ oneD%odi%nst2(ind1p)
den%vacxy(jz,ind1p-1,ivac,ispin) = den%vacxy(jz,ind1p-1,ivac,ispin) + CONJG(t1)/ oneD%odi%nst2(ind1p)
END IF
END DO xys1
......@@ -994,8 +991,8 @@ CONTAINS
uei = ue(jz,l,ispin)
uej = ue(jz,l1,ispin)
t1 = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
rhtxy(jz,ind2-1,ivac,ispin) = rhtxy(jz,ind2-1,ivac,ispin) + t1*phs/stars%nstr2(ind2)
rhtxy(jz,ind2p-1,ivac,ispin) = rhtxy(jz,ind2p-1,ivac,ispin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
den%vacxy(jz,ind2-1,ivac,ispin) = den%vacxy(jz,ind2-1,ivac,ispin) + t1*phs/stars%nstr2(ind2)
den%vacxy(jz,ind2p-1,ivac,ispin) = den%vacxy(jz,ind2p-1,ivac,ispin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
ENDDO
ENDDO
ENDDO
......@@ -1029,7 +1026,7 @@ CONTAINS
uj = u_1(jz,l1,m1,2)
uei = ue_1(jz,l,m,1)
uej = ue_1(jz,l1,m1,2)
cdomvz(jz,ivac) = cdomvz(jz,ivac) + aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
den%cdomvz(jz,ivac) = den%cdomvz(jz,ivac) + aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
END DO xys3
ELSE ! the warped part (ind1 > 1)
aa = CMPLX(0.,0.)
......@@ -1048,7 +1045,7 @@ CONTAINS
uei = ue_1(jz,l,m,1)
uej = ue_1(jz,l1,m1,2)
t1 = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
cdomvxy(jz,ind1-1,ivac) = cdomvxy(jz,ind1-1,ivac) + t1/oneD%odi%nst2(ind1)
den%vacxy(jz,ind1-1,ivac,3) = den%vacxy(jz,ind1-1,ivac,3) + t1/oneD%odi%nst2(ind1)
END DO xys2
END IF ! the non-warped (ind1 = 1)
END IF ! ind1.ne.0
......@@ -1086,7 +1083,7 @@ CONTAINS
ui2 = u(jz,l1,2)
uei = ue(jz,l,1)
uei2 = ue(jz,l1,2)
cdomvz(jz,ivac) = cdomvz(jz,ivac) + aa*ui2*ui + bb*uei2*uei + ab*ui2*uei + ba*uei2*ui
den%cdomvz(jz,ivac) = den%cdomvz(jz,ivac) + aa*ui2*ui + bb*uei2*uei + ab*ui2*uei + ba*uei2*ui
ENDDO
ELSE
!---> warping part
......@@ -1106,7 +1103,7 @@ CONTAINS
uei = ue(jz,l,1)
uej = ue(jz,l1,2)
t1 = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
cdomvxy(jz,ind2-1,ivac) = cdomvxy(jz, ind2-1,ivac) + t1*phs/stars%nstr2(ind2)
den%vacxy(jz,ind2-1,ivac,3) = den%vacxy(jz, ind2-1,ivac,3) + t1*phs/stars%nstr2(ind2)
ENDDO
ENDIF
ENDDO
......@@ -1144,10 +1141,10 @@ CONTAINS
uej = ue_1(jz,l1,m1,jspin)
t1 = aa*ui*uj + bb*uei*uej + ba*ui*uej + ab*uei*uj
IF (ind1.NE.0) THEN
rhtxy(jz,ind1-1,ivac,jspin) = rhtxy(jz,ind1-1,ivac,jspin) + t1/ oneD%odi%nst2(ind1)
den%vacxy(jz,ind1-1,ivac,jspin) = den%vacxy(jz,ind1-1,ivac,jspin) + t1/ oneD%odi%nst2(ind1)
END IF
IF (ind1p.NE.0) THEN
rhtxy(jz,ind1p-1,ivac,jspin) = rhtxy(jz,ind1p-1,ivac,jspin) + CONJG(t1)/ oneD%odi%nst2(ind1p)
den%vacxy(jz,ind1p-1,ivac,jspin) = den%vacxy(jz,ind1p-1,ivac,jspin) + CONJG(t1)/ oneD%odi%nst2(ind1p)
END IF
END DO xys
......@@ -1189,8 +1186,8 @@ CONTAINS
uei = ue(jz,l,jspin)
uej = ue(jz,l1,jspin)
t1 = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
rhtxy(jz,ind2-1,ivac,jspin) = rhtxy(jz,ind2-1, ivac,jspin) + t1*phs/stars%nstr2(ind2)
rhtxy(jz,ind2p-1,ivac,jspin) = rhtxy(jz,ind2p-1, ivac,jspin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
den%vacxy(jz,ind2-1,ivac,jspin) = den%vacxy(jz,ind2-1, ivac,jspin) + t1*phs/stars%nstr2(ind2)
den%vacxy(jz,ind2p-1,ivac,jspin) = den%vacxy(jz,ind2p-1, ivac,jspin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
ENDDO
ENDDO
END DO
......
......@@ -282,7 +282,7 @@ MODULE m_cdn_io
IF (datend == 0) THEN
IF (input%film) THEN
READ (iUnit) ((den%cdomvz(i,iVac),i=1,vacuum%nmz),iVac=1,vacuum%nvac)
READ (iUnit) (((den%cdomvxy(i,j-1,iVac),i=1,vacuum%nmzxy),j=2,oneD%odi%nq2), iVac=1,vacuum%nvac)
READ (iUnit) (((den%vacxy(i,j-1,iVac,3),i=1,vacuum%nmzxy),j=2,oneD%odi%nq2), iVac=1,vacuum%nvac)
END IF
ELSE
! (datend < 0) => no off-diagonal magnetisation stored
......@@ -294,14 +294,14 @@ MODULE m_cdn_io
den%cdom = CMPLX(0.0,0.0)
IF (input%film) THEN
den%cdomvz = CMPLX(0.0,0.0)
den%cdomvxy = CMPLX(0.0,0.0)
den%vacxy(:,:,:,3) = CMPLX(0.0,0.0)
END IF
END IF
ELSE IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
den%cdom = CMPLX(0.0,0.0)
IF (input%film) THEN
den%cdomvz = CMPLX(0.0,0.0)
den%cdomvxy = CMPLX(0.0,0.0)
den%vacxy(:,:,:,3) = CMPLX(0.0,0.0)
END IF
END IF
CLOSE(iUnit)
......@@ -598,7 +598,7 @@ MODULE m_cdn_io
WRITE (iUnit) (den%cdom(k),k=1,stars%ng3)
IF (input%film) THEN
WRITE (iUnit) ((den%cdomvz(i,iVac),i=1,vacuum%nmz),iVac=1,vacuum%nvac)
WRITE (iUnit) (((den%cdomvxy(i,j-1,iVac),i=1,vacuum%nmzxy),j=2,oneD%odi%nq2), iVac=1,vacuum%nvac)
WRITE (iUnit) (((den%vacxy(i,j-1,iVac,3),i=1,vacuum%nmzxy),j=2,oneD%odi%nq2), iVac=1,vacuum%nvac)
END IF
END IF
......
......@@ -1682,18 +1682,18 @@ MODULE m_cdnpot_io_hdf
dimsInt(:3)=(/2,ng3,input%jspins/)
CALL h5dopen_f(groupID, 'fpw', fpwSetID, hdfError)
CALL io_write_complex2(fpwSetID,(/-1,1,1/),dimsInt(:3),den%pw)
CALL io_write_complex2(fpwSetID,(/-1,1,1/),dimsInt(:3),den%pw(:,:input%jspins))
CALL h5dclose_f(fpwSetID, hdfError)
IF (l_film) THEN
dimsInt(:3)=(/nmzd,2,input%jspins/)
CALL h5dopen_f(groupID, 'fz', fzSetID, hdfError)
CALL io_write_real3(fzSetID,(/1,1,1/),dimsInt(:3),den%vacz)
CALL io_write_real3(fzSetID,(/1,1,1/),dimsInt(:3),den%vacz(:,:,:input%jspins))
CALL h5dclose_f(fzSetID, hdfError)
dimsInt(:5)=(/2,nmzxyd,ng2-1,2,input%jspins/)
CALL h5dopen_f(groupID, 'fzxy', fzxySetID, hdfError)
CALL io_write_complex4(fzxySetID,(/-1,1,1,1,1/),dimsInt(:5),den%vacxy)
CALL io_write_complex4(fzxySetID,(/-1,1,1,1,1/),dimsInt(:5),den%vacxy(:,:,:,:input%jspins))
CALL h5dclose_f(fzxySetID, hdfError)
END IF
......@@ -1713,7 +1713,7 @@ MODULE m_cdnpot_io_hdf
dimsInt(:4)=(/2,nmzxy,od_nq2-1,nvac/)
CALL h5dopen_f(groupID, 'cdomvxy', cdomvxySetID, hdfError)
CALL io_write_complex3(cdomvxySetID,(/-1,1,1,1/),dimsInt(:4),den%cdomvxy)
CALL io_write_complex3(cdomvxySetID,(/-1,1,1,1/),dimsInt(:4),den%vacxy(:,:,:,3))
CALL h5dclose_f(cdomvxySetID, hdfError)
END IF
END IF
......@@ -1745,7 +1745,7 @@ MODULE m_cdnpot_io_hdf
CALL h5screate_simple_f(3,dims(:3),fpwSpaceID,hdfError)
CALL h5dcreate_f(groupID, "fpw", H5T_NATIVE_DOUBLE, fpwSpaceID, fpwSetID, hdfError)
CALL h5sclose_f(fpwSpaceID,hdfError)
CALL io_write_complex2(fpwSetID,(/-1,1,1/),dimsInt(:3),den%pw)
CALL io_write_complex2(fpwSetID,(/-1,1,1/),dimsInt(:3),den%pw(:,:input%jspins))
CALL h5dclose_f(fpwSetID, hdfError)
IF (l_film) THEN
......@@ -1754,7 +1754,7 @@ MODULE m_cdnpot_io_hdf
CALL h5screate_simple_f(3,dims(:3),fzSpaceID,hdfError)
CALL h5dcreate_f(groupID, "fz", H5T_NATIVE_DOUBLE, fzSpaceID, fzSetID, hdfError)
CALL h5sclose_f(fzSpaceID,hdfError)
CALL io_write_real3(fzSetID,(/1,1,1/),dimsInt(:3),den%vacz)
CALL io_write_real3(fzSetID,(/1,1,1/),dimsInt(:3),den%vacz(:,:,:input%jspins))
CALL h5dclose_f(fzSetID, hdfError)
dims(:5)=(/2,nmzxyd,ng2-1,2,input%jspins/)
......@@ -1762,7 +1762,7 @@ MODULE m_cdnpot_io_hdf
CALL h5screate_simple_f(5,dims(:5),fzxySpaceID,hdfError)
CALL h5dcreate_f(groupID, "fzxy", H5T_NATIVE_DOUBLE, fzxySpaceID, fzxySetID, hdfError)
CALL h5sclose_f(fzxySpaceID,hdfError)
CALL io_write_complex4(fzxySetID,(/-1,1,1,1,1/),dimsInt(:5),den%vacxy)
CALL io_write_complex4(fzxySetID,(/-1,1,1,1,1/),dimsInt(:5),den%vacxy(:,:,:,:input%jspins))
CALL h5dclose_f(fzxySetID, hdfError)
END IF
......@@ -1791,7 +1791,7 @@ MODULE m_cdnpot_io_hdf
CALL h5screate_simple_f(4,dims(:4),cdomvxySpaceID,hdfError)
CALL h5dcreate_f(groupID, "cdomvxy", H5T_NATIVE_DOUBLE, cdomvxySpaceID, cdomvxySetID, hdfError)
CALL h5sclose_f(cdomvxySpaceID,hdfError)
CALL io_write_complex3(cdomvxySetID,(/-1,1,1,1/),dimsInt(:4),den%cdomvxy)
CALL io_write_complex3(cdomvxySetID,(/-1,1,1,1/),dimsInt(:4),den%vacxy(:,:,:,3))
CALL h5dclose_f(cdomvxySetID, hdfError)
END IF
END IF
......@@ -1842,7 +1842,7 @@ MODULE m_cdnpot_io_hdf
CALL h5screate_simple_f(3,dims(:3),fpwSpaceID,hdfError)
CALL h5dcreate_f(groupID, "fpw", H5T_NATIVE_DOUBLE, fpwSpaceID, fpwSetID, hdfError)
CALL h5sclose_f(fpwSpaceID,hdfError)
CALL io_write_complex2(fpwSetID,(/-1,1,1/),dimsInt(:3),den%pw)
CALL io_write_complex2