1!
2! Copyright (C) 2003 A. Smogunov
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC).
9!
10!
11SUBROUTINE poten(vppot,nrz,z)
12!
13! This subroutine computes the 2D Fourier components of the
14! local potential in each slab.
15!
16  USE constants, ONLY : tpi
17  USE cell_base, ONLY : at, bg
18  USE scf, only : vltot, v
19  USE noncollin_module, ONLY : noncolin, npol
20  USE cond
21  USE mp,               ONLY : mp_bcast
22  USE mp_world,         ONLY : world_comm
23  USE io_global,        ONLY : ionode_id
24  USE fft_scalar,       ONLY : cfft3d
25  USE fft_base,         ONLY : dfftp
26  USE scatter_mod,      ONLY : gather_grid
27
28  IMPLICIT NONE
29
30  INTEGER ::                                                &
31             i, j, ij, ijx, k, n, p, il, &
32             ix, jx, kx, nrz, info
33  INTEGER :: iis, jjs, is(4), js(4), ispin, nspin_eff
34  INTEGER, ALLOCATABLE :: ipiv(:)
35
36  REAL(DP), PARAMETER :: eps = 1.d-8
37  REAL(DP) :: arg, bet, z(nrz+1), zlen
38  REAL(DP), ALLOCATABLE :: gz(:), allv(:), auxr(:)
39
40  COMPLEX(DP), PARAMETER :: cim = (0.d0,1.d0)
41  COMPLEX(DP) :: caux, vppot(nrz,nrx*nry,npol,npol)
42  COMPLEX(DP), ALLOCATABLE :: aux(:), amat(:,:), amat0(:,:)
43  COMPLEX(DP), ALLOCATABLE :: vppot0(:,:,:,:)
44
45  CALL start_clock('poten')
46  ALLOCATE( ipiv( nrz ) )
47  ALLOCATE( gz( nrz ) )
48  ALLOCATE( aux( dfftp%nr1x*dfftp%nr2x*dfftp%nr3x ) )
49  ALLOCATE( auxr( dfftp%nnr ) )
50  ALLOCATE( amat( nrz, nrz ) )
51  ALLOCATE( amat0( nrz, nrz ) )
52
53
54  zlen = at(3,3)
55
56!
57!  Compute the Gz vectors in the z direction
58!
59  DO k = 1, nrz
60     il = k-1
61     IF (il.GT.nrz/2) il = il-nrz
62     gz(k) = il*bg(3,3)
63  ENDDO
64!
65! set up the matrix for the linear system
66!
67DO n=1,nrz
68   DO p=1,nrz
69      arg=gz(n)*z(p)*tpi
70      bet=gz(n)*(z(p+1)-z(p))*tpi
71      IF (ABS(gz(n)).GT.eps) THEN
72        caux=cim*(CMPLX(COS(bet),-SIN(bet),kind=DP)-(1.d0,0.d0))  &
73                                    /zlen/gz(n)/tpi
74      ELSE
75        caux=(z(p+1)-z(p))/zlen
76      ENDIF
77      amat0(n,p)=CMPLX(COS(arg),-SIN(arg),kind=DP)*caux
78   ENDDO
79ENDDO
80IF (noncolin) THEN
81   nspin_eff=4
82   ij=0
83   DO iis=1,2
84      DO jjs=1,2
85         ij=ij+1
86         is(ij)=iis
87         js(ij)=jjs
88      ENDDO
89   ENDDO
90ELSE
91   nspin_eff=1
92   is(1)=1
93   js(1)=1
94ENDIF
95!
96!     To form local potential on the real space mesh
97!
98!
99#if defined(__MPI)
100  allocate ( allv(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x) )
101#endif
102
103vppot = 0.d0
104DO ispin=1,nspin_eff
105   IF (noncolin) THEN
106      IF (ispin==1) THEN
107         auxr(:) = vltot(:)+v%of_r(:,1)
108      ELSE
109         auxr(:) = v%of_r(:,ispin)
110      ENDIF
111   ELSE
112      auxr(:) = vltot(:) + v%of_r(:,iofspin)
113   ENDIF
114!
115! To collect the potential from different CPUs
116!
117#if defined(__MPI)
118  call gather_grid ( dfftp, auxr, allv )
119  CALL mp_bcast( allv, ionode_id, world_comm )
120  aux(:) = CMPLX(allv(:), 0.d0,kind=DP)
121#else
122  aux(:) = CMPLX(auxr(:), 0.d0,kind=DP)
123#endif
124!
125!  To find FFT of the local potential
126!  (use serial FFT even in the parallel case)
127!
128  CALL cfft3d (aux,dfftp%nr1,dfftp%nr2,dfftp%nr3,dfftp%nr1x,dfftp%nr2x,dfftp%nr3x,1,-1)
129
130  DO i = 1, nrx
131    IF(i.GT.nrx/2+1) THEN
132        ix = dfftp%nr1-(nrx-i)
133    ELSE
134        ix = i
135    ENDIF
136    DO j = 1, nry
137      IF(j.GT.nry/2+1) THEN
138         jx = dfftp%nr2-(nry-j)
139      ELSE
140         jx = j
141      ENDIF
142      ij = i+(j-1)*nrx
143      ijx = ix+(jx-1)*dfftp%nr1x
144
145      DO k = 1, nrz
146        il = k-1
147        IF (il.GT.nrz/2) il = il-nrz
148        IF(il.LE.dfftp%nr3/2.AND.il.GE.-(dfftp%nr3-1)/2) THEN
149
150         IF(k.GT.nrz/2+1) THEN
151            kx = dfftp%nr3-(nrz-k)
152         ELSE
153            kx = k
154         ENDIF
155         vppot(k, ij, is(ispin), js(ispin)) = aux(ijx+(kx-1)*dfftp%nr1x*dfftp%nr2x)
156
157        ENDIF
158      ENDDO
159    ENDDO
160  ENDDO
161!
162! solve the linear system
163!
164  amat=amat0
165  CALL ZGESV(nrz, nrx*nry, amat, nrz, ipiv, vppot(1,1,is(ispin),js(ispin)),&
166                                             nrz, info)
167  CALL errore ('poten','info different from zero',ABS(info))
168ENDDO
169
170IF (noncolin) THEN
171   ALLOCATE( vppot0(nrz, nrx * nry, npol, npol) )
172   vppot0=vppot
173   vppot(:,:,1,1)=vppot0(:,:,1,1)+vppot0(:,:,2,2)
174   vppot(:,:,1,2)=vppot0(:,:,1,2)-(0.d0,1.d0)*vppot0(:,:,2,1)
175   vppot(:,:,2,1)=vppot0(:,:,1,2)+(0.d0,1.d0)*vppot0(:,:,2,1)
176   vppot(:,:,2,2)=vppot0(:,:,1,1)-vppot0(:,:,2,2)
177   DEALLOCATE( vppot0 )
178ENDIF
179
180!  do p = 1, nrz
181!    write(stdout,'(i5,2f12.6)') p, real(vppot(p,1,1,1)), imag(vppot(p,1,1,1))
182!  enddo
183!  stop
184
185  DEALLOCATE(ipiv)
186  DEALLOCATE(gz)
187  DEALLOCATE(aux)
188  DEALLOCATE(auxr)
189  DEALLOCATE(amat)
190  DEALLOCATE(amat0)
191#if defined(__MPI)
192  deallocate(allv)
193#endif
194
195  CALL stop_clock('poten')
196
197  RETURN
198END SUBROUTINE poten
199