Commit 0b4393da authored by Gregor Michalicek's avatar Gregor Michalicek

Introduce t_potden type to brysh1, brysh2

..also in this commit: slight cleanup in mix/metric.f90
parent 03e00312
......@@ -153,15 +153,11 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
!(in the spin polarized case the arrays sm and fsm consist of spin up and spin down densities)
CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
intfac,vacfac,inDen%pw,inDen%mt,inDen%vacz,inDen%vacxy,inDen%cdom,&
inDen%cdomvz,inDen%cdomvxy,inDen%mmpMat(-lmaxU_const,-lmaxU_const,1,1),&
nmap,nmaph,mapmt,mapvac,mapvac2,sm)
intfac,vacfac,inDen,nmap,nmaph,mapmt,mapvac,mapvac2,sm)
!put output charge density into array fsm
CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
intfac,vacfac,outDen%pw,outDen%mt,outDen%vacz,outDen%vacxy,outDen%cdom,&
outDen%cdomvz,outDen%cdomvxy,outDen%mmpMat(-lmaxU_const,-lmaxU_const,1,1),&
nmap,nmaph,mapmt,mapvac,mapvac2,fsm)
intfac,vacfac,outDen,nmap,nmaph,mapmt,mapvac,mapvac2,fsm)
!store fsm - sm the difference on fsm
DO imap = 1,nmap
......@@ -196,9 +192,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
inDen%cdomvxy = CMPLX(0.0,0.0)
inDen%mmpMat = CMPLX(0.0,0.0)
CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,inDen%mmpMat,oneD,&
inDen%pw,inDen%mt,inDen%vacz,inDen%vacxy,inDen%cdom,&
inDen%cdomvz,inDen%cdomvxy)
CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,oneD,inDen)
!calculate the distance of charge densities...
......
......@@ -8,7 +8,7 @@ MODULE m_brysh1
!******************************************************
CONTAINS
SUBROUTINE brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
intfac,vacfac,qpw,rho,rht,rhtxy,cdom,cdomvz,cdomvxy,n_mmp, nmap,nmaph,mapmt,mapvac,mapvac2,sout)
intfac,vacfac,den,nmap,nmaph,mapmt,mapvac,mapvac2,sout)
USE m_types
IMPLICIT NONE
......@@ -20,25 +20,18 @@ CONTAINS
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
! ..
! .. Scalar Arguments ..
REAL, INTENT (IN) :: intfac,vacfac
INTEGER, INTENT (OUT):: mapmt,mapvac,mapvac2,nmap,nmaph
! ..
! .. Array Arguments ..
!-odim
!+odim
COMPLEX, INTENT (IN) :: qpw(stars%ng3,input%jspins)
COMPLEX, INTENT (IN) :: cdomvz(vacuum%nmz,2),cdomvxy(vacuum%nmzxy,oneD%odi%n2d-1,2)
COMPLEX, INTENT (IN) :: cdom(stars%ng3),rhtxy(vacuum%nmzxy,oneD%odi%n2d-1,2,input%jspins)
COMPLEX, INTENT (IN) :: n_mmp(-3:3,-3:3,atoms%n_u,input%jspins)
REAL, INTENT (IN) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
REAL, INTENT (IN) :: rht(vacuum%nmz,2,input%jspins)
REAL, INTENT (OUT):: sout(:)
! ..
! .. Local Scalars ..
TYPE(t_potden),INTENT(IN) :: den
! Scalar Arguments
REAL, INTENT (IN) :: intfac,vacfac
INTEGER, INTENT (OUT) :: mapmt,mapvac,mapvac2,nmap,nmaph
! Array Arguments
REAL, INTENT (OUT) :: sout(:)
! Local Scalars
INTEGER i,iv,j,js,k,l,n,nall,na,nvaccoeff,nvaccoeff2,mapmtd
!
!---> put input into arrays sout
! in the spin polarized case the arrays consist of
! spin up and spin down densities
......@@ -47,12 +40,12 @@ CONTAINS
DO js = 1,input%jspins
DO i = 1,stars%ng3
j = j + 1
sout(j) = REAL(qpw(i,js))
sout(j) = REAL(den%pw(i,js))
END DO
IF (.NOT.sym%invs) THEN
DO i = 1,stars%ng3
j = j + 1
sout(j) = AIMAG(qpw(i,js))
sout(j) = AIMAG(den%pw(i,js))
END DO
ENDIF
mapmt=0
......@@ -62,7 +55,7 @@ CONTAINS
DO i = 1,atoms%jri(n)
mapmt = mapmt +1
j = j + 1
sout(j) = rho(i,l,n,js)
sout(j) = den%mt(i,l,n,js)
END DO
END DO
na = na + atoms%neq(n)
......@@ -73,13 +66,13 @@ CONTAINS
DO k = 1,vacuum%nmz
mapvac = mapvac + 1
j = j + 1
sout(j) = rht(k,iv,js)
sout(j) = den%vacz(k,iv,js)
END DO
DO k = 1,oneD%odi%nq2-1
DO i = 1,vacuum%nmzxy
mapvac = mapvac + 1
j = j + 1
sout(j) = REAL(rhtxy(i,k,iv,js))
sout(j) = REAL(den%vacxy(i,k,iv,js))
END DO
END DO
IF (.NOT.sym%invs2) THEN
......@@ -87,7 +80,7 @@ CONTAINS
DO i = 1,vacuum%nmzxy
mapvac = mapvac + 1
j = j + 1
sout(j) = AIMAG(rhtxy(i,k,iv,js))
sout(j) = AIMAG(den%vacxy(i,k,iv,js))
END DO
END DO
END IF
......@@ -101,24 +94,24 @@ CONTAINS
!---> off-diagonal part of the density matrix
DO i = 1,stars%ng3
j = j + 1
sout(j) = REAL(cdom(i))
sout(j) = REAL(den%cdom(i))
END DO
DO i = 1,stars%ng3
j = j + 1
sout(j) = AIMAG(cdom(i))
sout(j) = AIMAG(den%cdom(i))
END DO
IF (input%film) THEN
DO iv = 1,vacuum%nvac
DO k = 1,vacuum%nmz
mapvac2 = mapvac2 + 1
j = j + 1
sout(j) = REAL(cdomvz(k,iv))
sout(j) = REAL(den%cdomvz(k,iv))
END DO
DO k = 1,oneD%odi%nq2-1
DO i = 1,vacuum%nmzxy
mapvac2 = mapvac2 + 1
j = j + 1
sout(j) = REAL(cdomvxy(i,k,iv))
sout(j) = REAL(den%cdomvxy(i,k,iv))
END DO
END DO
END DO
......@@ -126,13 +119,13 @@ CONTAINS
DO k = 1,vacuum%nmz
mapvac2 = mapvac2 + 1
j = j + 1
sout(j) = AIMAG(cdomvz(k,iv))
sout(j) = AIMAG(den%cdomvz(k,iv))
END DO
DO k = 1,oneD%odi%nq2-1
DO i = 1,vacuum%nmzxy
mapvac2 = mapvac2 + 1
j = j + 1
sout(j) = AIMAG(cdomvxy(i,k,iv))
sout(j) = AIMAG(den%cdomvxy(i,k,iv))
END DO
END DO
END DO
......@@ -154,9 +147,9 @@ CONTAINS
DO k = -3, 3
DO i = -3, 3
j = j + 1
sout(j) = REAL(n_mmp(i,k,n,js))
sout(j) = REAL(den%mmpMat(i,k,n,js))
j = j + 1
sout(j) = AIMAG(n_mmp(i,k,n,js))
sout(j) = AIMAG(den%mmpMat(i,k,n,js))
ENDDO
ENDDO
ENDDO
......
......@@ -5,45 +5,37 @@ MODULE m_brysh2
! proper component of interstitial, m.t. and vacuum density
!******************************************************
CONTAINS
SUBROUTINE brysh2(&
& input,stars,atoms,sphhar,&
& noco,vacuum,&
& sym,s_in,&
& n_mmp,oneD,qpw,rho,rht,rhtxy,cdom,cdomvz,cdomvxy)
SUBROUTINE brysh2(input,stars,atoms,sphhar,noco,vacuum,&
sym,s_in,oneD,den)
USE m_types
IMPLICIT NONE
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
! ..
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden),INTENT(INOUT) :: den
REAL, INTENT (IN) :: s_in(:)
REAL, INTENT (OUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
REAL, INTENT (OUT) :: rht(vacuum%nmz,2,input%jspins)
COMPLEX, INTENT (OUT) :: qpw(stars%ng3,input%jspins),cdom(stars%ng3),cdomvz(vacuum%nmz,2)
COMPLEX, INTENT (OUT) :: rhtxy(vacuum%nmzxy,oneD%odi%n2d-1,2,input%jspins)
COMPLEX, INTENT (OUT) :: cdomvxy(vacuum%nmzxy,oneD%odi%n2d-1,2)
COMPLEX, INTENT (OUT) :: n_mmp(-3:3,-3:3,atoms%n_u,input%jspins)
! ..
! .. Local Scalars ..
! Local Scalars
INTEGER i,iv,j,js,k,l,n,na
!
j=0
DO js = 1,input%jspins
IF (sym%invs) THEN
DO i = 1,stars%ng3
j = j + 1
qpw(i,js) = CMPLX(s_in(j),0.0)
den%pw(i,js) = CMPLX(s_in(j),0.0)
END DO
ELSE
DO i = 1,stars%ng3
j = j + 1
qpw(i,js) = CMPLX(s_in(j),s_in(j+stars%ng3))
den%pw(i,js) = CMPLX(s_in(j),s_in(j+stars%ng3))
END DO
j = j + stars%ng3
ENDIF
......@@ -52,7 +44,7 @@ CONTAINS
DO l = 0,sphhar%nlh(atoms%ntypsy(na))
DO i = 1,atoms%jri(n)
j = j + 1
rho(i,l,n,js) = s_in(j)
den%mt(i,l,n,js) = s_in(j)
END DO
END DO
na = na + atoms%neq(n)
......@@ -61,19 +53,19 @@ CONTAINS
DO iv = 1,vacuum%nvac
DO k = 1,vacuum%nmz
j = j + 1
rht(k,iv,js) = s_in(j)
den%vacz(k,iv,js) = s_in(j)
END DO
DO k = 1,oneD%odi%nq2-1
DO i = 1,vacuum%nmzxy
j = j + 1
rhtxy(i,k,iv,js) = CMPLX(s_in(j),0.0)
den%vacxy(i,k,iv,js) = CMPLX(s_in(j),0.0)
END DO
END DO
IF (.NOT.sym%invs2) THEN
DO k = 1,oneD%odi%nq2-1
DO i = 1,vacuum%nmzxy
j = j + 1
rhtxy(i,k,iv,js) = rhtxy(i,k,iv,js) +CMPLX(0.0,s_in(j))
den%vacxy(i,k,iv,js) = den%vacxy(i,k,iv,js) + CMPLX(0.0,s_in(j))
END DO
END DO
END IF
......@@ -85,34 +77,34 @@ CONTAINS
!---> off-diagonal part of the density matrix
DO i = 1,stars%ng3
j = j + 1
cdom(i) = CMPLX(s_in(j),0.0)
den%cdom(i) = CMPLX(s_in(j),0.0)
END DO
DO i = 1,stars%ng3
j = j + 1
cdom(i) = cdom(i) + CMPLX(0.0,s_in(j))
den%cdom(i) = den%cdom(i) + CMPLX(0.0,s_in(j))
END DO
IF (input%film) THEN
DO iv = 1,vacuum%nvac
DO k = 1,vacuum%nmz
j = j + 1
cdomvz(k,iv) = CMPLX(s_in(j),0.0)
den%cdomvz(k,iv) = CMPLX(s_in(j),0.0)
END DO
DO k = 1,oneD%odi%nq2-1
DO i = 1,vacuum%nmzxy
j = j + 1
cdomvxy(i,k,iv) = CMPLX(s_in(j),0.0)
den%cdomvxy(i,k,iv) = CMPLX(s_in(j),0.0)
END DO
END DO
END DO
DO iv = 1,vacuum%nvac
DO k = 1,vacuum%nmz
j = j + 1
cdomvz(k,iv) = cdomvz(k,iv) + CMPLX(0.0,s_in(j))
den%cdomvz(k,iv) = den%cdomvz(k,iv) + CMPLX(0.0,s_in(j))
END DO
DO k = 1,oneD%odi%nq2-1
DO i = 1,vacuum%nmzxy
j = j + 1
cdomvxy(i,k,iv) = cdomvxy(i,k,iv)+ CMPLX(0.0,s_in(j))
den%cdomvxy(i,k,iv) = den%cdomvxy(i,k,iv)+ CMPLX(0.0,s_in(j))
END DO
END DO
END DO
......@@ -125,7 +117,7 @@ CONTAINS
DO k = -3, 3
DO i = -3, 3
j = j + 1
n_mmp(i,k,n,js) = CMPLX(s_in(j),s_in(j+1))
den%mmpMat(i,k,n,js) = CMPLX(s_in(j),s_in(j+1))
j = j + 1
ENDDO
ENDDO
......
......@@ -5,229 +5,211 @@ MODULE m_metric
! output vector sout
!*********************************************************
CONTAINS
SUBROUTINE metric(&
& cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
& mmap,nmaph,mapmt,mapvac2,s_in,sout,lpot)
SUBROUTINE metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
mmap,nmaph,mapmt,mapvac2,s_in,sout,lpot)
USE m_metrz0
USE m_convol
USE m_types
IMPLICIT NONE
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: mmap
INTEGER, INTENT (IN) :: mapmt,mapvac2
INTEGER, INTENT (IN) :: nmaph
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
! Scalar Arguments
INTEGER, INTENT (IN) :: mmap
INTEGER, INTENT (IN) :: mapmt,mapvac2
INTEGER, INTENT (IN) :: nmaph
LOGICAL, OPTIONAL,INTENT (IN) :: lpot !do we mix a potential??
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: s_in(mmap)
REAL, INTENT (OUT):: sout(mmap)
! ..
! .. Local Scalars ..
INTEGER :: imap,ivac,iz,j,js,k2,l,n,iv2c,iv2,na,ioff
REAL :: dvol,dxn,dxn2,dxn4,volnstr2
LOGICAL :: l_pot
! ..
! .. Local Arrays ..
! Array Arguments
REAL, INTENT (IN) :: s_in(mmap)
REAL, INTENT (OUT) :: sout(mmap)
! Local Scalars
INTEGER :: imap,ivac,iz,j,js,k2,l,n,iv2c,iv2,na,ioff
REAL :: dvol,dxn,dxn2,dxn4,volnstr2
LOGICAL :: l_pot
! Local Arrays
REAL, ALLOCATABLE :: g(:),wght(:)
COMPLEX, ALLOCATABLE :: ag3(:),fg3(:)
!
! calculate the coefficients of the g-matrix
! for the m.t. and vacuum region
!
IF (PRESENT(lpot)) THEN
l_pot=lpot
ELSE
l_pot=.FALSE.
ENDIF
! calculate the coefficients of the g-matrix
! for the m.t. and vacuum region
l_pot = .FALSE.
IF (PRESENT(lpot)) l_pot = lpot ! for potential mixing
ALLOCATE (g(mmap),wght(vacuum%nmzd),ag3(stars%ng3),fg3(stars%ng3))
g=0.0
IF (sym%invs) THEN
imap=stars%ng3
ELSE
imap=2*stars%ng3
ENDIF
iv2=1
IF (.NOT.sym%invs2) iv2=2
!
g = 0.0
imap = 2 * stars%ng3 ! complex values without invs
IF (sym%invs) imap = stars%ng3 ! only real values with invs
iv2 = 2
IF (sym%invs2) iv2 = 1
! metric for MT is r^2 dr/di di = r^3 dx ; divided by r^4 to
! compensate use of n(r) * r^2 in array rho
! simpson integration used, weights for first and last point: 1
! weights forthe rest alternating: 2 or 4
!
na = 1
DO n = 1,atoms%ntype
dxn = atoms%neq(n)*atoms%dx(n)/3.0e0
dxn2 =2.0e0 *dxn
dxn4 =4.0e0 *dxn
!
DO l = 0,sphhar%nlh(atoms%ntypsy(na))
DO n = 1, atoms%ntype
dxn = atoms%neq(n) * atoms%dx(n) / 3.0e0
dxn2 =2.0e0 * dxn
dxn4 =4.0e0 * dxn
DO l = 0, sphhar%nlh(atoms%ntypsy(na))
imap = imap + 1
g(imap) = dxn/atoms%rmsh(1,n)
g(imap) = dxn / atoms%rmsh(1,n)
IF (.NOT.l_pot) THEN
DO j = 2,atoms%jri(n)-1,2
DO j = 2, atoms%jri(n) - 1, 2
imap = imap + 2
g(imap-1) = dxn4/atoms%rmsh(j,n)
g(imap) = dxn2/atoms%rmsh(j+1,n)
ENDDO
g(imap-1) = dxn4 / atoms%rmsh(j,n)
g(imap) = dxn2 / atoms%rmsh(j+1,n)
END DO
! CHANGE JR 96/12/01
! take care when jri(n) is even
imap=imap+1-MOD(atoms%jri(n),2)
g(imap) = dxn/atoms%rmsh(atoms%jri(n),n)
imap = imap + 1 - MOD(atoms%jri(n),2)
g(imap) = dxn / atoms%rmsh(atoms%jri(n),n)
ELSE
!
! for the potential multiply by r^4
!
DO j = 2,atoms%jri(n)-1,2
! for the potential multiply by r^4
DO j = 2, atoms%jri(n) - 1, 2
imap = imap + 2
g(imap-1) = dxn4*atoms%rmsh(j,n)**3
g(imap) = dxn2*atoms%rmsh(j+1,n)**3
ENDDO
imap=imap+1-MOD(atoms%jri(n),2)
g(imap) = dxn*atoms%rmsh(atoms%jri(n),n)**3
ENDIF
ENDDO
g(imap-1) = dxn4 * atoms%rmsh(j,n)**3
g(imap) = dxn2 * atoms%rmsh(j+1,n)**3
END DO
imap = imap + 1 - MOD(atoms%jri(n),2)
g(imap) = dxn * atoms%rmsh(atoms%jri(n),n)**3
END IF
END DO
na = na + atoms%neq(n)
enddo
!
END DO
! vacuum contribution
!
IF (input%film) THEN
dvol = cell%area*vacuum%delz
! nvac=1 if (zrfs.or.invs)
IF ( vacuum%nvac.EQ.1 ) dvol = dvol + dvol
IF (oneD%odi%d1) THEN
dvol = cell%area*vacuum%delz
END IF
DO ivac = 1,vacuum%nvac
! G||=0 components
!
!---> use 7-point simpson integration in accordance to intgz0.f
! calculate weights for integration
! nvac=1 if (zrfs.or.invs)
IF (vacuum%nvac.EQ.1) dvol = dvol + dvol
IF (oneD%odi%d1) dvol = cell%area*vacuum%delz
DO ivac = 1, vacuum%nvac
! G||=0 components
!
! use 7-point simpson integration in accordance to intgz0.f
! calculate weights for integration
CALL metr_z0(vacuum%nmz,wght)
DO iz = 1,vacuum%nmz
DO iz = 1, vacuum%nmz
imap = imap + 1
IF (oneD%odi%d1) THEN
g(imap) = wght(iz)*dvol*(cell%z1+(iz-1)*vacuum%delz)
g(imap) = wght(iz) * dvol * (cell%z1+(iz-1)*vacuum%delz)
ELSE
g(imap) = wght(iz)*dvol
ENDIF
ENDDO
! G||.ne.0 components
! calculate weights for integration
g(imap) = wght(iz) * dvol
END IF
END DO
! G||.ne.0 components
!
! calculate weights for integration
CALL metr_z0(vacuum%nmzxy,wght)
DO iv2c=1,iv2
DO k2 = 1,oneD%odi%nq2-1
DO iv2c = 1, iv2
DO k2 = 1, oneD%odi%nq2 - 1
IF (oneD%odi%d1) THEN
DO iz = 1,vacuum%nmzxy
imap = imap + 1
g(imap) = wght(iz)*oneD%odi%nst2(k2)*&
& dvol*(cell%z1+(iz-1)*vacuum%delz)
ENDDO
g(imap) = wght(iz) * oneD%odi%nst2(k2) * dvol * (cell%z1+(iz-1)*vacuum%delz)
END DO
ELSE
volnstr2= dvol*stars%nstr2(k2)
DO iz = 1,vacuum%nmzxy
volnstr2 = dvol * stars%nstr2(k2)
DO iz = 1, vacuum%nmzxy
imap = imap + 1
g(imap) = wght(iz)*volnstr2
ENDDO
g(imap) = wght(iz) * volnstr2
END DO
END IF
ENDDO
ENDDO
enddo
END DO
END DO
END DO
END IF
!
! multiplicate the metric with the vector
!
DO js = 1,input%jspins
! map s_in on a complex help array ag3
! Interstitial region (metric here = step function)
! multiplicate the metric with the vector in real space
DO js = 1, input%jspins
! map s_in on a complex help array ag3
IF (sym%invs) THEN
DO imap = 1,stars%ng3
DO imap = 1, stars%ng3
ag3(imap) = CMPLX(s_in(imap+nmaph*(js-1)),0.0)
ENDDO
END DO
ELSE
DO imap = 1,stars%ng3
ag3(imap) = CMPLX(s_in(imap+nmaph*(js-1)),&
& s_in(imap+stars%ng3+nmaph*(js-1)))
ENDDO
DO imap = 1, stars%ng3
ag3(imap) = CMPLX(s_in(imap+nmaph*(js-1)),s_in(imap+stars%ng3+nmaph*(js-1)))
END DO
ENDIF
CALL convol(&
& stars,&
& fg3,&
& ag3,stars%ufft)
CALL convol(stars,fg3,ag3,stars%ufft)
IF (sym%invs) THEN
DO imap = 1,stars%ng3
sout(imap+nmaph*(js-1)) = cell%omtil*REAL(fg3(imap))
ENDDO
DO imap = stars%ng3+1,nmaph
sout(imap+nmaph*(js-1)) = g(imap)*s_in(imap+nmaph*(js-1))
ENDDO
DO imap = 1, stars%ng3
sout(imap+nmaph*(js-1)) = cell%omtil * REAL(fg3(imap))
END DO
DO imap = stars%ng3+1, nmaph
sout(imap+nmaph*(js-1)) = g(imap) * s_in(imap+nmaph*(js-1))
END DO
ELSE
DO imap = 1,stars%ng3
sout(imap+nmaph*(js-1)) = cell%omtil*REAL(fg3(imap))
sout(imap+stars%ng3+nmaph*(js-1)) = cell%omtil*AIMAG(fg3(imap))
ENDDO
DO imap = 2*stars%ng3+1,nmaph
sout(imap+nmaph*(js-1)) = g(imap)*s_in(imap+nmaph*(js-1))
ENDDO
ENDIF
enddo
DO imap = 1, stars%ng3
sout(imap+nmaph*(js-1)) = cell%omtil * REAL(fg3(imap))
sout(imap+stars%ng3+nmaph*(js-1)) = cell%omtil * AIMAG(fg3(imap))
END DO
DO imap = 2 * stars%ng3 + 1, nmaph
sout(imap+nmaph*(js-1)) = g(imap) * s_in(imap+nmaph*(js-1))
END DO
END IF
END DO
IF (noco%l_noco) THEN
DO imap = 1,stars%ng3
ag3(imap) = CMPLX(s_in(2*nmaph + imap),&
& s_in(2*nmaph + stars%ng3 + imap))
ENDDO
CALL convol(&
& stars,&
& fg3,&
& ag3,stars%ufft)
DO imap = 1,stars%ng3
sout(2*nmaph + imap) = cell%omtil*REAL(fg3(imap))
sout(2*nmaph + stars%ng3 + imap) = cell%omtil*AIMAG(fg3(imap))
ENDDO
DO imap = 1, stars%ng3
ag3(imap) = CMPLX(s_in(2*nmaph + imap),s_in(2*nmaph + stars%ng3 + imap))
END DO
CALL convol(stars,fg3,ag3,stars%ufft)
DO imap = 1, stars%ng3
sout(2*nmaph + imap) = cell%omtil * REAL(fg3(imap))
sout(2*nmaph + stars%ng3 + imap) = cell%omtil * AIMAG(fg3(imap))
END DO
IF (input%film) THEN
!---> js runs over the real and imaginary part of the vacuum density
!---> coefficients (not the spin).
IF (sym%invs) THEN
ioff = stars%ng3 + mapmt
ELSE
ioff = 2*stars%ng3 + mapmt
ENDIF
DO js = 1,2
DO imap = 1,mapvac2/2
! js runs over the real and imaginary part of the vacuum density
! coefficients (not the spin).
ioff = 2*stars%ng3 + mapmt ! (for complex values)
IF (sym%invs) ioff = stars%ng3 + mapmt ! (real values if sym%invs)
DO js = 1, 2
DO imap = 1, mapvac2 / 2
sout(2*nmaph + 2*stars%ng3 + mapvac2/2*(js-1) + imap) = &
& g(ioff + imap)&
& *s_in(2*nmaph + 2*stars%ng3 + mapvac2/2*(js-1) + imap)
ENDDO
ENDDO
ENDIF
ENDIF