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