1! 2! Copyright (C) 2002-2011 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 10!=----------------------------------------------------------------------------=! 11 subroutine phfac_x( tau0, ei1, ei2, ei3, eigr) 12!=----------------------------------------------------------------------------=! 13 ! 14 ! this subroutine generates the complex matrices ei1, ei2, and ei3 15 ! used to compute the structure factor and forces on atoms : 16 ! ei1(n1,ia,is) = exp(-i*n1*b1*tau(ia,is)) -nr1<n1<nr1 17 ! and similar definitions for ei2 and ei3 ; and : 18 ! eigr(n,ia,is) = ei1*ei2*ei3 = exp(-i g*tau(ia,is)) 19 ! The value of n1,n2,n3 for a vector g is supplied by array mill 20 ! calculated in ggen . 21 ! 22 use kinds, only: DP 23 use control_flags, only: iverbosity 24 use io_global, only: stdout 25 use ions_base, only: nsp, na, nat 26 use cell_base, only: ainv, r_to_s 27 use fft_base, only: dfftp 28 use gvect, only: mill 29 use gvecw, only: ngw 30 use cp_interfaces, only: phfacs 31! 32 implicit none 33 real(DP) tau0(3,nat) 34! 35 complex(DP) ei1(-dfftp%nr1:dfftp%nr1,nat), ei2(-dfftp%nr2:dfftp%nr2,nat), & 36 & ei3(-dfftp%nr3:dfftp%nr3,nat), eigr(ngw,nat) 37! 38 integer :: i, isa 39 real(DP), allocatable :: taus(:,:) 40! 41 allocate( taus(3,nat) ) 42! 43 if( iverbosity > 2 ) then 44 WRITE( stdout,*) ' phfac: tau0 ' 45 WRITE( stdout,*) ( ( tau0(i,isa), i=1, 3 ), isa=1, nat ) 46 endif 47 CALL r_to_s( tau0, taus, nat, ainv ) 48 CALL phfacs( ei1, ei2, ei3, eigr, mill, taus, dfftp%nr1, dfftp%nr2, dfftp%nr3, nat ) 49 50 deallocate( taus ) 51! 52 return 53 end subroutine phfac_x 54 55 56!=----------------------------------------------------------------------------=! 57 SUBROUTINE phfacs_x( ei1, ei2, ei3, eigr, mill, taus, nr1, nr2, nr3, nat ) 58!=----------------------------------------------------------------------------=! 59 60 ! this routine computes the phase factors 61 ! 62 ! ei1(ix,ia) = exp ( -i ix G_1 dot R(ia)) 63 ! ei2(iy,ia) = exp ( -i iy G_2 dot R(ia)) 64 ! ei3(iz,ia) = exp ( -i iz G_3 dot R(ia)) 65 ! 66 ! eigr(ig,ia) = exp( -i G dot R(ia)) = 67 ! = ei1(ix,ia) * ei2(iy,ia) * ei3(iz,ia) 68 ! 69 ! G_1,G_2,G_3 = reciprocal lattice generators 70 ! 71 ! ia = index of ion 72 ! ig = index of G vector 73 ! ix,iy,iz = Miller indices 74 ! ---------------------------------------------- 75 76 USE kinds, ONLY: DP 77 USE constants, ONLY: tpi 78 79 IMPLICIT NONE 80 81 ! ... declare subroutine arguments 82 83 INTEGER, INTENT(IN) :: nat 84 INTEGER, INTENT(IN) :: nr1, nr2, nr3 85 COMPLEX(DP) :: ei1( -nr1 : nr1, nat ) 86 COMPLEX(DP) :: ei2( -nr2 : nr2, nat ) 87 COMPLEX(DP) :: ei3( -nr3 : nr3, nat ) 88 COMPLEX(DP) :: eigr( :, : ) 89 REAL(DP) :: taus( 3, nat ) 90 INTEGER :: mill( :, : ) 91 92 ! ... declare other variables 93 94 COMPLEX(DP) :: ctep1, ctep2, ctep3, ctem1, ctem2, ctem3 95 REAL(DP) :: ar1, ar2, ar3 96 INTEGER :: i, j, k, isa 97 INTEGER :: ig, ig1, ig2, ig3, ngw 98 99 ! ... --+ end of declarations +-- 100 101 if(nr1 < 3) call errore(' phfacs ',' nr1 too small ',1) 102 if(nr2 < 3) call errore(' phfacs ',' nr2 too small ',1) 103 if(nr3 < 3) call errore(' phfacs ',' nr3 too small ',1) 104 105 DO isa = 1, nat 106 107 ! ... Miller index = 0: exp(i 0 dot R(ia)) = 1 108 109 ei1( 0, isa ) = ( 1.d0, 0.d0 ) 110 ei2( 0, isa ) = ( 1.d0, 0.d0 ) 111 ei3( 0, isa ) = ( 1.d0, 0.d0 ) 112 113 ! ... let R_1,R_2,R_3 be the direct lattice generators, 114 ! ... G_1,G_2,G_3 the reciprocal lattice generators 115 ! ... by definition G_i dot R_j = 2 pi delta_{ij} 116 ! ... ionic coordinates are in units of R_1,R_2,R_3 117 ! ... then G_i dot R(ia) = 2 pi R(ia)_i 118 119 ar1 = tpi * taus(1,isa) ! G_1 dot R(ia) 120 ar2 = tpi * taus(2,isa) ! G_2 dot R(ia) 121 ar3 = tpi * taus(3,isa) ! G_3 dot R(ia) 122 123 ! ... Miller index = 1: exp(-i G_i dot R(ia)) 124 125 ctep1 = CMPLX( cos( ar1 ), -sin( ar1 ) ,kind=DP) 126 ctep2 = CMPLX( cos( ar2 ), -sin( ar2 ) ,kind=DP) 127 ctep3 = CMPLX( cos( ar3 ), -sin( ar3 ) ,kind=DP) 128 129 ! ... Miller index = -1: exp(-i G_im dot R(ia)) = exp(i G_i dot R(ia)) 130 131 ctem1 = CONJG(ctep1) 132 ctem2 = CONJG(ctep2) 133 ctem3 = CONJG(ctep3) 134 135 ! ... Miller index > 0: exp(i N G_i dot R(ia)) = 136 ! ... = exp(i G_i dot R(ia)) exp(i (N-1) G_i dot R(ia)) 137 ! ... Miller index < 0: exp(-i N G_i dot R(ia)) = 138 ! ... = exp(-i G_i dot R(ia)) exp(-i (N-1) G_i dot R(ia)) 139 140 DO i = 1, nr1 141 ei1( i, isa ) = ei1( i - 1, isa ) * ctep1 142 ei1( -i, isa ) = ei1( -i + 1, isa ) * ctem1 143 END DO 144 DO j = 1, nr2 145 ei2( j, isa ) = ei2( j - 1, isa ) * ctep2 146 ei2( -j, isa ) = ei2( -j + 1, isa ) * ctem2 147 END DO 148 DO k = 1, nr3 149 ei3( k, isa ) = ei3( k - 1, isa ) * ctep3 150 ei3( -k, isa ) = ei3( -k + 1, isa ) * ctem3 151 END DO 152 153 END DO 154 155 ngw = SIZE( eigr, 1 ) 156 IF( ngw > SIZE( mill, 2 ) ) THEN 157 CALL errore(' phfacs ',' inconsistent size for eigr ',ngw) 158 END IF 159 160 DO ig = 1, ngw 161 ig1 = mill( 1, ig ) 162 ig2 = mill( 2, ig ) 163 ig3 = mill( 3, ig ) 164 DO i = 1, nat 165 eigr( ig, i ) = ei1( ig1, i ) * ei2( ig2, i ) * ei3( ig3, i ) 166 END DO 167 END DO 168 169 RETURN 170 END SUBROUTINE phfacs_x 171 172 173!=----------------------------------------------------------------------------=! 174 SUBROUTINE strucf_x( sfac, ei1, ei2, ei3, mill, ngm ) 175!=----------------------------------------------------------------------------=! 176 177! this routine computes the structure factors 178! 179! sfac(ig,is) = (sum over ia) exp(i G dot R(ia)) = 180! (sum over ia) ei1(ix,ia) * ei2(iy,ia) * ei3(iz,ia) 181! 182! ei1(ix,ia) = exp (i ix G_1 dot R(ia)) 183! ei2(iy,ia) = exp (i iy G_2 dot R(ia)) 184! ei3(iz,ia) = exp (i iz G_3 dot R(ia)) 185! 186! G_1,G_2,G_3 = reciprocal lattice generators 187! 188! ia = index of ion (running over ions of species is) 189! ig = index of G vector 190! is = index of atomic species 191! ix,iy,iz = Miller indices of G vector 192 193 194 USE kinds, ONLY: DP 195 USE ions_base, ONLY: nat, na, nsp, ityp 196 use fft_base, only: dfftp 197 198 IMPLICIT NONE 199 200 ! ... declare subroutine arguments 201 ! 202 COMPLEX(DP) :: ei1( -dfftp%nr1 : dfftp%nr1, nat ) 203 COMPLEX(DP) :: ei2( -dfftp%nr2 : dfftp%nr2, nat ) 204 COMPLEX(DP) :: ei3( -dfftp%nr3 : dfftp%nr3, nat ) 205 INTEGER :: mill( :, : ) 206 INTEGER :: ngm 207 COMPLEX(DP), INTENT(OUT) :: sfac(:,:) 208 209 ! ... declare other variables 210 ! 211 INTEGER :: is, ig, ia, ig1, ig2, ig3 212 213 call start_clock( 'strucf' ) 214 215 sfac = (0.0d0, 0.0d0) 216 217 DO ia = 1, nat 218 is = ityp(ia) 219 DO ig = 1, ngm 220 ig1 = mill( 1, ig ) 221 ig2 = mill( 2, ig ) 222 ig3 = mill( 3, ig ) 223 sfac( ig, is ) = sfac( ig, is ) + ei1( ig1, ia ) * ei2( ig2, ia ) * ei3( ig3, ia ) 224 END DO 225 END DO 226 227 call stop_clock( 'strucf' ) 228 229 RETURN 230 END SUBROUTINE strucf_x 231 232 233