1!
2! Copyright (C) 2001-2009 Quantum ESPRESSO group
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!
9!-----------------------------------------------------------------------
10SUBROUTINE do_shift_ew (alat, nat, ntyp, ityp, zv, delta_zv, at, bg, tau, &
11     omega, g, gg, ngm, gcutm, gstart, gamma_only, shift_ion)
12  !-----------------------------------------------------------------------
13  !
14  ! Calculates Ewald energy with both G- and R-space terms.
15  ! Determines optimal alpha. Should hopefully work for any structure.
16  !
17  !
18  USE kinds, ONLY : DP
19  USE constants, ONLY : tpi, e2
20  USE mp_pools,  ONLY : intra_pool_comm
21  USE mp, ONLY : mp_sum
22  IMPLICIT NONE
23  !
24  !   first the dummy variables
25  !
26
27  INTEGER :: nat, ntyp, ityp (nat), ngm, gstart
28  ! input: number of atoms in the unit cell
29  ! input: number of different types of atoms
30  ! input: the type of each atom
31  ! input: number of plane waves for G sum
32  ! input: first non-zero G vector
33
34  LOGICAL :: gamma_only
35
36  real(DP) :: tau (3, nat), g (3, ngm), gg (ngm), zv (ntyp), &
37       at (3, 3), bg (3, 3), omega, alat, gcutm, delta_zv(ntyp),  &
38       shift_ion(nat)
39  ! input: the positions of the atoms in the cell
40  ! input: the coordinates of G vectors
41  ! input: the square moduli of G vectors
42  ! input: the charge of each type of atoms
43  ! input: the direct lattice vectors
44  ! input: the reciprocal lattice vectors
45  ! input: the volume of the unit cell
46  ! input: lattice parameter
47  ! input: cut-off of g vectors
48  real(DP) :: ewald
49  ! output: the ewald energy
50  !
51  !    here the local variables
52  !
53  INTEGER, PARAMETER :: mxr = 50
54  ! the maximum number of R vectors included in r
55  INTEGER :: ng, nr, na, nb, nt, nrm, ipol
56  ! counter over reciprocal G vectors
57  ! counter over direct vectors
58  ! counter on atoms
59  ! counter on atoms
60  ! counter on atomic types
61  ! number of R vectors included in r sum
62  ! counter on polarization
63
64  real(DP) :: charge, tpiba2, ewaldg, ewaldr, dtau (3), alpha, &
65       r (3, mxr), r2 (mxr), rmax, rr, upperbound, fact, arg
66  ! total ionic charge in the cell
67  ! length in reciprocal space
68  ! ewald energy computed in reciprocal space
69  ! ewald energy computed in real space
70  ! the difference tau_s - tau_s'
71  ! alpha term in ewald sum
72  ! input of the rgen routine ( not used here )
73  ! the square modulus of R_j-tau_s-tau_s'
74  ! the maximum radius to consider real space sum
75  ! buffer variable
76  ! used to optimize alpha
77  COMPLEX(DP), ALLOCATABLE :: rhon(:)
78  real(DP), EXTERNAL :: qe_erfc
79
80  ALLOCATE (rhon(ngm))
81
82  shift_ion(:) = 0.d0
83
84  tpiba2 = (tpi / alat) **2
85  charge = 0.d0
86  DO na = 1, nat
87     charge = charge+zv (ityp (na) )
88  ENDDO
89  alpha = 2.9d0
90100 alpha = alpha - 0.1d0
91  !
92  ! choose alpha in order to have convergence in the sum over G
93  ! upperbound is a safe upper bound for the error in the sum over G
94  !
95  IF (alpha<=0.d0) CALL errore ('do_shift_ew', 'optimal alpha not found', 1)
96  upperbound = 2.d0 * charge**2 * sqrt (2.d0 * alpha / tpi) * qe_erfc ( &
97       sqrt (tpiba2 * gcutm / 4.d0 / alpha) )
98  IF (upperbound>1.0d-7) GOTO 100
99  !
100  ! G-space sum here.
101  ! Determine if this processor contains G=0 and set the constant term
102  !
103  IF (gstart==2) THEN
104     DO na =1,nat
105        shift_ion(na) = - charge * delta_zv(ityp(na)) /alpha/ 4.0d0
106     ENDDO
107  ENDIF
108  IF (gamma_only) THEN
109     fact = 2.d0
110  ELSE
111     fact = 1.d0
112  ENDIF
113  DO ng = gstart, ngm
114     rhon(ng) = (0.d0, 0.d0)
115     DO na =1, nat
116        arg = (g (1, ng) * tau (1, na) + &
117               g (2, ng) * tau (2, na) + &
118               g (3, ng) * tau (3, na) ) * tpi
119        rhon(ng) = rhon(ng) + zv (ityp(na)) * cmplx(cos (arg), -sin (arg),kind=DP)
120     ENDDO
121  ENDDO
122  DO na=1,nat
123     DO ng=gstart, ngm
124        arg = (g (1, ng) * tau (1, na) + g (2, ng) * tau (2, na) &
125             + g (3, ng) * tau (3, na) ) * tpi
126        shift_ion(na) = shift_ion(na) + fact * delta_zv(ityp(na)) * &
127                        conjg(rhon(ng)) * cmplx(cos (arg), -sin (arg),kind=DP) * &
128                        exp ( -gg(ng)*tpiba2/alpha/4.d0) / gg(ng)/tpiba2
129     ENDDO
130  ENDDO
131  shift_ion(:) = 2.d0 * tpi / omega * shift_ion(:)
132  !
133  !  Here add the other constant term
134  !
135  IF (gstart==2) THEN
136     DO na = 1, nat
137        shift_ion(na) = shift_ion(na) - &
138                        zv (ityp (na) ) * delta_zv(ityp(na)) *  &
139                        sqrt (8.d0/tpi*alpha)
140     ENDDO
141  ENDIF
142  !
143  ! R-space sum here (only for the processor that contains G=0)
144  !
145  IF (gstart==2) THEN
146     rmax = 4.d0 / sqrt (alpha) / alat
147     !
148     ! with this choice terms up to ZiZj*erfc(4) are counted (erfc(4)=2x10^-8
149     !
150     DO na = 1, nat
151        DO nb = 1, nat
152           DO ipol = 1, 3
153              dtau (ipol) = tau (ipol, na) - tau (ipol, nb)
154           ENDDO
155           !
156           ! generates nearest-neighbors shells
157           !
158           CALL rgen (dtau, rmax, mxr, at, bg, r, r2, nrm)
159           !
160           ! and sum to the real space part
161           !
162           DO nr = 1, nrm
163              rr = sqrt (r2 (nr) ) * alat
164              shift_ion(na) = shift_ion(na) + &
165                       delta_zv(ityp(na)) * zv (ityp (nb) ) *  &
166                       qe_erfc ( sqrt (alpha) * rr) / rr
167           ENDDO
168        ENDDO
169     ENDDO
170  ENDIF
171
172  shift_ion(:) = e2 * shift_ion(:)
173
174  CALL mp_sum ( shift_ion, intra_pool_comm )
175
176  DEALLOCATE (rhon)
177  RETURN
178END SUBROUTINE do_shift_ew
179
180