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