1!----------------------------------------------------------------------- 2!------- DRIVERS FOR DERIVATIVES OF XC POTENTIAL (GGA CASE) ------------ 3!----------------------------------------------------------------------- 4! 5!--------------------------------------------------------------------- 6SUBROUTINE dgcxc( length, sp, r_in, g_in, dvxc_rr, dvxc_sr, dvxc_ss ) 7 !--------------------------------------------------------------------- 8 !! Wrapper routine. Calls dgcx-driver routines from internal libraries 9 !! or from the external libxc, depending on the input choice. 10 ! 11 USE constants, ONLY: e2 12 USE kinds, ONLY: DP 13 USE funct, ONLY: get_igcx, get_igcc, is_libxc 14 USE xc_gga, ONLY: gcxc, gcx_spin 15#if defined(__LIBXC) 16#include "xc_version.h" 17 USE xc_f03_lib_m 18#endif 19 ! 20 IMPLICIT NONE 21 ! 22 INTEGER, INTENT(IN) :: length 23 !! length of the I/O arrays 24 INTEGER, INTENT(IN) :: sp 25 !! number of spin components 26 REAL(DP), INTENT(IN) :: r_in(length,sp) 27 !! charge density 28 REAL(DP), INTENT(IN) :: g_in(length,3,sp) 29 !! gradient 30 REAL(DP), INTENT(OUT) :: dvxc_rr(length,sp,sp), dvxc_sr(length,sp,sp), & 31 dvxc_ss(length,sp,sp) 32 ! 33 ! ... local variables 34 ! 35 REAL(DP), ALLOCATABLE :: vrrx(:,:), vsrx(:,:), vssx(:,:) 36 REAL(DP), ALLOCATABLE :: vrrc(:,:), vsrc(:,:), vssc(:), vrzc(:,:) 37 ! 38#if defined(__LIBXC) 39 TYPE(xc_f03_func_t) :: xc_func 40 TYPE(xc_f03_func_info_t) :: xc_info1, xc_info2 41 INTEGER :: fkind 42 REAL(DP), ALLOCATABLE :: rho_lbxc(:) 43 REAL(DP), ALLOCATABLE :: v2rho2_x(:), v2rhosigma_x(:), v2sigma2_x(:) 44 REAL(DP), ALLOCATABLE :: v2rho2_c(:), v2rhosigma_c(:), v2sigma2_c(:) 45#if (XC_MAJOR_VERSION > 4) 46 INTEGER(8) :: lengthxc 47#else 48 INTEGER :: lengthxc 49#endif 50#endif 51 ! 52 INTEGER :: igcx, igcc 53 INTEGER :: k, ir, length_lxc, length_dlxc 54 REAL(DP) :: rht, zeta 55 REAL(DP), ALLOCATABLE :: sigma(:) 56 REAL(DP), PARAMETER :: small = 1.E-10_DP, rho_trash = 0.5_DP 57 REAL(DP), PARAMETER :: epsr=1.0d-6, epsg=1.0d-10 58 ! 59 igcx = get_igcx() 60 igcc = get_igcc() 61 ! 62#if defined(__LIBXC) 63 ! 64 fkind = -1 65 lengthxc = length 66 ! 67 IF ( (is_libxc(3) .OR. igcx==0) .AND. (is_libxc(4) .OR. igcc==0)) THEN 68 ! 69 length_lxc = length*sp 70 length_dlxc = length 71 IF (sp == 2) length_dlxc = length*3 72 ! 73 ALLOCATE( rho_lbxc(length_lxc), sigma(length_dlxc) ) 74 ! 75 ! ... set libxc input 76 SELECT CASE( sp ) 77 CASE( 1 ) 78 ! 79 DO k = 1, length 80 rho_lbxc(k) = r_in(k,1) 81 sigma(k) = g_in(k,1,1)**2 + g_in(k,2,1)**2 + g_in(k,3,1)**2 82 ENDDO 83 ! 84 CASE( 2 ) 85 ! 86 DO k = 1, length 87 rho_lbxc(2*k-1) = r_in(k,1) 88 rho_lbxc(2*k) = r_in(k,2) 89 ! 90 sigma(3*k-2) = g_in(k,1,1)**2 + g_in(k,2,1)**2 + g_in(k,3,1)**2 91 sigma(3*k-1) = g_in(k,1,1) * g_in(k,1,2) + g_in(k,2,1) * g_in(k,2,2) + & 92 g_in(k,3,1) * g_in(k,3,2) 93 sigma(3*k) = g_in(k,1,2)**2 + g_in(k,2,2)**2 + g_in(k,3,2)**2 94 ENDDO 95 ! 96 CASE( 4 ) 97 ! 98 CALL errore( 'dgcxc', 'The derivative of the xc potential with libxc & 99 &is not available for noncollinear case', 1 ) 100 ! 101 CASE DEFAULT 102 ! 103 CALL errore( 'dgcxc', 'Wrong number of spin dimensions', 2 ) 104 ! 105 END SELECT 106 ! 107 ALLOCATE( v2rho2_x(length_dlxc), v2rhosigma_x(length_dlxc*sp), v2sigma2_x(length_dlxc*sp) ) 108 ALLOCATE( v2rho2_c(length_dlxc), v2rhosigma_c(length_dlxc*sp), v2sigma2_c(length_dlxc*sp) ) 109 ! 110 ! ... DERIVATIVE FOR EXCHANGE 111 v2rho2_x = 0.d0 ; v2rhosigma_x = 0.d0 ; v2sigma2_x = 0.d0 112 IF (igcx /= 0) THEN 113 CALL xc_f03_func_init( xc_func, igcx, sp ) 114 xc_info1 = xc_f03_func_get_info( xc_func ) 115 fkind = xc_f03_func_info_get_kind( xc_info1 ) 116 CALL xc_f03_gga_fxc( xc_func, lengthxc, rho_lbxc(1), sigma(1), v2rho2_x(1), v2rhosigma_x(1), v2sigma2_x(1) ) 117 CALL xc_f03_func_end( xc_func ) 118 ENDIF 119 ! 120 ! ... DERIVATIVE FOR CORRELATION 121 v2rho2_c = 0.d0 ; v2rhosigma_c = 0.d0 ; v2sigma2_c = 0.d0 122 IF (igcc /= 0) THEN 123 CALL xc_f03_func_init( xc_func, igcc, sp ) 124 xc_info2 = xc_f03_func_get_info( xc_func ) 125 CALL xc_f03_gga_fxc( xc_func, lengthxc, rho_lbxc(1), sigma(1), v2rho2_c(1), v2rhosigma_c(1), v2sigma2_c(1) ) 126 CALL xc_f03_func_end( xc_func ) 127 ENDIF 128 ! 129 dvxc_rr = 0.d0 130 dvxc_sr = 0.d0 131 dvxc_ss = 0.d0 132 ! 133 IF ( sp == 1 ) THEN 134 ! 135 dvxc_rr(:,1,1) = e2 * (v2rho2_x(:) + v2rho2_c(:)) 136 dvxc_sr(:,1,1) = e2 * (v2rhosigma_x(:) + v2rhosigma_c(:))*2.d0 137 dvxc_ss(:,1,1) = e2 * (v2sigma2_x(:) + v2sigma2_c(:))*4.d0 138 ! 139 ELSEIF ( sp == 2 ) THEN 140 ! 141 DO k = 1, length 142 rht = r_in(k,1) + r_in(k,2) 143 IF (rht > epsr) THEN 144 dvxc_rr(k,1,1) = e2 * (v2rho2_x(3*k-2) + v2rho2_c(3*k-2)) 145 dvxc_rr(k,1,2) = e2 * (v2rho2_x(3*k-1) + v2rho2_c(3*k-1)) 146 dvxc_rr(k,2,1) = e2 * (v2rho2_x(3*k-1) + v2rho2_c(3*k-1)) 147 dvxc_rr(k,2,2) = e2 * (v2rho2_x(3*k) + v2rho2_c(3*k)) 148 ENDIF 149 ! 150 dvxc_sr(k,1,1) = e2 * (v2rhosigma_x(6*k-5) + v2rhosigma_c(6*k-5))*2.d0 151 dvxc_ss(k,1,1) = e2 * (v2sigma2_x(6*k-5) + v2sigma2_c(6*k))*4.d0 152 IF ( fkind==XC_EXCHANGE_CORRELATION ) THEN 153 dvxc_sr(k,1,2) = e2 * v2rhosigma_x(6*k-4) 154 dvxc_sr(k,2,1) = e2 * v2rhosigma_x(6*k-1) 155 dvxc_ss(k,1,2) = e2 * v2sigma2_x(6*k-2) 156 dvxc_ss(k,2,1) = e2 * v2sigma2_x(6*k-2) 157 ELSE 158 dvxc_sr(k,1,2) = e2 * v2rhosigma_c(6*k-4) 159 dvxc_sr(k,2,1) = e2 * v2rhosigma_c(6*k-1) 160 dvxc_ss(k,1,2) = e2 * v2sigma2_c(6*k-2) 161 dvxc_ss(k,2,1) = e2 * v2sigma2_c(6*k-2) 162 ENDIF 163 dvxc_sr(k,2,2) = e2 * (v2rhosigma_x(6*k) + v2rhosigma_c(6*k))*2.d0 164 dvxc_ss(k,2,2) = e2 * (v2sigma2_x(6*k) + v2sigma2_c(6*k))*4.d0 165 ! 166 ENDDO 167 ! 168 ENDIF 169 ! 170 DEALLOCATE( rho_lbxc, sigma ) 171 DEALLOCATE( v2rho2_x, v2rhosigma_x, v2sigma2_x ) 172 DEALLOCATE( v2rho2_c, v2rhosigma_c, v2sigma2_c ) 173 ! 174 ELSEIF ((.NOT.is_libxc(3)) .AND. (.NOT.is_libxc(4))) THEN 175 ! 176 ALLOCATE( vrrx(length,sp), vsrx(length,sp), vssx(length,sp) ) 177 ALLOCATE( vrrc(length,sp), vsrc(length,sp), vssc(length) ) 178 ! 179 IF ( sp == 1 ) THEN 180 ! 181 ALLOCATE( sigma(length) ) 182 sigma(:) = g_in(:,1,1)**2 + g_in(:,2,1)**2 + g_in(:,3,1)**2 183 CALL dgcxc_unpol( length, r_in, sigma, vrrx, vsrx, vssx, vrrc, vsrc, vssc ) 184 DEALLOCATE( sigma ) 185 ! 186 dvxc_rr(:,1,1) = e2 * (vrrx(:,1) + vrrc(:,1)) 187 dvxc_sr(:,1,1) = e2 * (vsrx(:,1) + vsrc(:,1)) 188 dvxc_ss(:,1,1) = e2 * (vssx(:,1) + vssc(:) ) 189 ! 190 ELSEIF ( sp == 2 ) THEN 191 ! 192 ALLOCATE( vrzc(length,sp) ) 193 ! 194 CALL dgcxc_spin( length, r_in, g_in, vrrx, vsrx, vssx, vrrc, vsrc, vssc, vrzc ) 195 ! 196 DO k = 1, length 197 ! 198 rht = r_in(k,1) + r_in(k,2) 199 IF (rht > epsr) THEN 200 zeta = (r_in(k,1) - r_in(k,2)) / rht 201 ! 202 dvxc_rr(k,1,1) = e2 * (vrrx(k,1) + vrrc(k,1) + vrzc(k,1) * (1.d0 - zeta) / rht) 203 dvxc_rr(k,1,2) = e2 * (vrrc(k,1) - vrzc(k,1) * (1.d0 + zeta) / rht) 204 dvxc_rr(k,2,1) = e2 * (vrrc(k,2) + vrzc(k,2) * (1.d0 - zeta) / rht) 205 dvxc_rr(k,2,2) = e2 * (vrrx(k,2) + vrrc(k,2) - vrzc(k,2) * (1.d0 + zeta) / rht) 206 ENDIF 207 ! 208 dvxc_sr(k,1,1) = e2 * (vsrx(k,1) + vsrc(k,1)) 209 dvxc_sr(k,1,2) = e2 * vsrc(k,1) 210 dvxc_sr(k,2,1) = e2 * vsrc(k,2) 211 dvxc_sr(k,2,2) = e2 * (vsrx(k,2) + vsrc(k,2)) 212 ! 213 dvxc_ss(k,1,1) = e2 * (vssx(k,1) + vssc(k)) 214 dvxc_ss(k,1,2) = e2 * vssc(k) 215 dvxc_ss(k,2,1) = e2 * vssc(k) 216 dvxc_ss(k,2,2) = e2 * (vssx(k,2) + vssc(k)) 217 ! 218 ENDDO 219 ! 220 DEALLOCATE( vrzc ) 221 ! 222 ENDIF 223 ! 224 DEALLOCATE( vrrx, vsrx, vssx ) 225 DEALLOCATE( vrrc, vsrc, vssc ) 226 ! 227 ELSE 228 ! 229 CALL errore( 'dgcxc', 'Derivatives of exchange and correlation terms, & 230 & at present, must be both qe or both libxc.', 3 ) 231 ! 232 ENDIF 233 ! 234#else 235 ! 236 ALLOCATE( vrrx(length,sp), vsrx(length,sp), vssx(length,sp) ) 237 ALLOCATE( vrrc(length,sp), vsrc(length,sp), vssc(length) ) 238 ! 239 SELECT CASE( sp ) 240 CASE( 1 ) 241 ! 242 ALLOCATE( sigma(length) ) 243 sigma(:) = g_in(:,1,1)**2 + g_in(:,2,1)**2 + g_in(:,3,1)**2 244 CALL dgcxc_unpol( length, r_in, sigma, vrrx, vsrx, vssx, vrrc, vsrc, vssc ) 245 DEALLOCATE( sigma ) 246 ! 247 dvxc_rr(:,1,1) = e2 * (vrrx(:,1) + vrrc(:,1)) 248 dvxc_sr(:,1,1) = e2 * (vsrx(:,1) + vsrc(:,1)) 249 dvxc_ss(:,1,1) = e2 * (vssx(:,1) + vssc(:) ) 250 ! 251 CASE( 2 ) 252 ! 253 ALLOCATE( vrzc(length,sp) ) 254 ! 255 CALL dgcxc_spin( length, r_in, g_in, vrrx, vsrx, vssx, vrrc, vsrc, vssc, vrzc ) 256 ! 257 DO k = 1, length 258 ! 259 rht = r_in(k,1) + r_in(k,2) 260 IF (rht > epsr) THEN 261 zeta = (r_in(k,1) - r_in(k,2)) / rht 262 ! 263 dvxc_rr(k,1,1) = e2 * (vrrx(k,1) + vrrc(k,1) + vrzc(k,1) * (1.d0 - zeta) / rht) 264 dvxc_rr(k,1,2) = e2 * (vrrc(k,1) - vrzc(k,1) * (1.d0 + zeta) / rht) 265 dvxc_rr(k,2,1) = e2 * (vrrc(k,2) + vrzc(k,2) * (1.d0 - zeta) / rht) 266 dvxc_rr(k,2,2) = e2 * (vrrx(k,2) + vrrc(k,2) - vrzc(k,2) * (1.d0 + zeta) / rht) 267 ENDIF 268 ! 269 dvxc_sr(k,1,1) = e2 * (vsrx(k,1) + vsrc(k,1)) 270 dvxc_sr(k,1,2) = e2 * vsrc(k,1) 271 dvxc_sr(k,2,1) = e2 * vsrc(k,2) 272 dvxc_sr(k,2,2) = e2 * (vsrx(k,2) + vsrc(k,2)) 273 ! 274 dvxc_ss(k,1,1) = e2 * (vssx(k,1) + vssc(k)) 275 dvxc_ss(k,1,2) = e2 * vssc(k) 276 dvxc_ss(k,2,1) = e2 * vssc(k) 277 dvxc_ss(k,2,2) = e2 * (vssx(k,2) + vssc(k)) 278 ENDDO 279 ! 280 DEALLOCATE( vrzc ) 281 ! 282 CASE DEFAULT 283 ! 284 CALL errore( 'dgcxc', 'Wrong ns input', 4 ) 285 ! 286 END SELECT 287 ! 288 DEALLOCATE( vrrx, vsrx, vssx ) 289 DEALLOCATE( vrrc, vsrc, vssc ) 290 ! 291#endif 292 ! 293 ! 294 RETURN 295 ! 296END SUBROUTINE 297! 298! 299!--------------------------------------------------------------------------- 300SUBROUTINE dgcxc_unpol( length, r_in, s2_in, vrrx, vsrx, vssx, vrrc, vsrc, vssc ) 301 !------------------------------------------------------------------------- 302 !! This routine computes the derivative of the exchange and correlation 303 !! potentials. 304 ! 305 USE kinds, ONLY: DP 306 USE xc_gga, ONLY: gcxc 307 ! 308 IMPLICIT NONE 309 ! 310 INTEGER, INTENT(IN) :: length 311 REAL(DP), INTENT(IN), DIMENSION(length) :: r_in, s2_in 312 REAL(DP), INTENT(OUT), DIMENSION(length) :: vrrx, vsrx, vssx 313 REAL(DP), INTENT(OUT), DIMENSION(length) :: vrrc, vsrc, vssc 314 ! 315 ! ... local variables 316 ! 317 INTEGER :: i1, i2, i3, i4, f1, f2, f3, f4 318 REAL(DP), DIMENSION(length) :: dr, s, ds 319 REAL(DP), DIMENSION(4*length) :: raux, s2aux 320 REAL(DP), ALLOCATABLE :: v1x(:), v2x(:), v1c(:), v2c(:) 321 REAL(DP), ALLOCATABLE :: sx(:), sc(:) 322 REAL(DP), PARAMETER :: small = 1.E-30_DP 323 ! 324 ALLOCATE( v1x(4*length), v2x(4*length), sx(4*length) ) 325 ALLOCATE( v1c(4*length), v2c(4*length), sc(4*length) ) 326 ! 327 i1 = 1 ; f1 = length !4 blocks: [ rho+dr , grho2 ] 328 i2 = f1+1 ; f2 = 2*length ! [ rho-dr , grho2 ] 329 i3 = f2+1 ; f3 = 3*length ! [ rho , (grho+ds)^2 ] 330 i4 = f3+1 ; f4 = 4*length ! [ rho , (grho-ds)^2 ] 331 ! 332 s = SQRT(s2_in) 333 dr = MIN(1.d-4, 1.d-2*r_in) 334 ds = MIN(1.d-4, 1.d-2*s) 335 ! 336 raux(i1:f1) = r_in+dr ; s2aux(i1:f1) = s2_in 337 raux(i2:f2) = r_in-dr ; s2aux(i2:f2) = s2_in 338 raux(i3:f3) = r_in ; s2aux(i3:f3) = (s+ds)**2 339 raux(i4:f4) = r_in ; s2aux(i4:f4) = (s-ds)**2 340 ! 341 CALL gcxc( length*4, raux, s2aux, sx, sc, v1x, v2x, v1c, v2c ) 342 ! 343 ! ... to avoid NaN in the next operations 344 WHERE( r_in<=small .OR. s2_in<=small ) 345 dr = 1._DP ; ds = 1._DP ; s = 1._DP 346 END WHERE 347 ! 348 vrrx = 0.5_DP * (v1x(i1:f1) - v1x(i2:f2)) / dr 349 vrrc = 0.5_DP * (v1c(i1:f1) - v1c(i2:f2)) / dr 350 ! 351 vsrx = 0.25_DP * ((v2x(i1:f1) - v2x(i2:f2)) / dr + & 352 (v1x(i3:f3) - v1x(i4:f4)) / ds / s) 353 vsrc = 0.25_DP * ((v2c(i1:f1) - v2c(i2:f2)) / dr + & 354 (v1c(i3:f3) - v1c(i4:f4)) / ds / s) 355 ! 356 vssx = 0.5_DP * (v2x(i3:f3) - v2x(i4:f4)) / ds / s 357 vssc = 0.5_DP * (v2c(i3:f3) - v2c(i4:f4)) / ds / s 358 ! 359 DEALLOCATE( v1x, v2x, sx ) 360 DEALLOCATE( v1c, v2c, sc ) 361 ! 362 RETURN 363 ! 364END SUBROUTINE dgcxc_unpol 365! 366! 367!-------------------------------------------------------------------------- 368SUBROUTINE dgcxc_spin( length, r_in, g_in, vrrx, vrsx, vssx, vrrc, vrsc, & 369 vssc, vrzc ) 370 !------------------------------------------------------------------------ 371 !! This routine computes the derivative of the exchange and correlation 372 !! potentials in the spin-polarized case. 373 ! 374 USE xc_gga, ONLY: gcx_spin, gcc_spin 375 USE kinds, ONLY: DP 376 ! 377 IMPLICIT NONE 378 ! 379 INTEGER, INTENT(IN) :: length 380 REAL(DP), INTENT(IN), DIMENSION(length,2) :: r_in 381 REAL(DP), INTENT(IN), DIMENSION(length,3,2) :: g_in 382 ! input: the charges and the gradient 383 REAL(DP), INTENT(OUT), DIMENSION(length,2) :: vrrx, vrsx, vssx 384 REAL(DP), INTENT(OUT), DIMENSION(length,2) :: vrrc, vrsc, vrzc 385 REAL(DP), INTENT(OUT), DIMENSION(length) :: vssc 386 ! output: derivatives of the exchange and of the correlation 387 ! 388 ! ... local variables 389 ! 390 INTEGER :: i1, i2, i3, i4, i5, i6, i7, i8 391 INTEGER :: f1, f2, f3, f4, f5, f6, f7, f8 392 ! block delimiters 393 REAL(DP), DIMENSION(length,2) :: r, s, s2 394 REAL(DP), DIMENSION(length,2) :: drup, drdw, dsup, dsdw 395 ! deltas for rho and gradient 396 REAL(DP), ALLOCATABLE :: sx(:), v1x(:,:), v2x(:,:) 397 ! exchange energy and potentials for each block 398 REAL(DP), ALLOCATABLE :: sc(:), v1c(:,:), v2c(:) 399 ! correlation energy and potentials for each block 400 REAL(DP), DIMENSION(length) :: rt, zeta, st, s2t 401 ! rho tot, zeta, gradient, square tot gradient 402 REAL(DP), DIMENSION(length) :: dr, ds, dz 403 ! deltas for rho tot, gradient and zeta 404 REAL(DP), DIMENSION(length,2) :: null_v 405 ! used to set output values to zero when input values 406 ! are too small (e.g. rho<eps) 407 ! 408 REAL(DP), ALLOCATABLE :: raux(:,:), s2aux(:,:) 409 REAL(DP), ALLOCATABLE :: rtaux(:), s2taux(:), zetaux(:) 410 ! auxiliary arrays for gcx- and gcc- routines input 411 ! 412 REAL(DP), PARAMETER :: eps = 1.D-6 413 REAL(DP), PARAMETER :: rho_trash = 0.4_DP, zeta_trash = 0.2_DP, & 414 s2_trash = 0.1_DP 415 ! 416 vrrx = 0.0_DP ; vrsx = 0.0_DP ; vssx = 0.0_DP 417 vrrc = 0.0_DP ; vrsc = 0.0_DP ; vrzc = 0.0_DP 418 vssc = 0.0_DP 419 ! 420 ! ... EXCHANGE 421 ! 422 i1 = 1 ; f1 = length !8 blocks(x2): [ rho+drup , grho2 ] 423 i2 = f1+1 ; f2 = 2*length ! [ rho-drup , grho2 ] 424 i3 = f2+1 ; f3 = 3*length ! [ rho , (grho+dsup)^2 ] 425 i4 = f3+1 ; f4 = 4*length ! [ rho , (grho-dsup)^2 ] 426 i5 = f4+1 ; f5 = 5*length ! [ rho+drdw , grho2 ] 427 i6 = f5+1 ; f6 = 6*length ! [ rho-drdw , grho2 ] 428 i7 = f6+1 ; f7 = 7*length ! [ rho , (grho+dsdw)^2 ] 429 i8 = f7+1 ; f8 = 8*length ! [ rho , (grho-dsdw)^2 ] 430 ! 431 ALLOCATE( raux(length*8,2), s2aux(length*8,2) ) 432 ALLOCATE( sx(length*8), v1x(length*8,2), v2x(length*8,2) ) 433 ! 434 s2(:,1) = g_in(:,1,1)**2 + g_in(:,2,1)**2 + g_in(:,3,1)**2 !up 435 s2(:,2) = g_in(:,1,2)**2 + g_in(:,2,2)**2 + g_in(:,3,2)**2 !down 436 s = SQRT(s2) 437 r = r_in 438 ! 439 null_v = 1.0_DP 440 ! 441 ! ... thresholds 442 WHERE ( r_in(:,1)<=eps .OR. s(:,1)<=eps ) 443 r(:,1) = rho_trash 444 s2(:,1) = s2_trash ; s(:,1) = SQRT(s2_trash) 445 null_v(:,1) = 0.0_DP 446 END WHERE 447 ! 448 WHERE ( r_in(:,2)<=eps .OR. s(:,2)<=eps ) 449 r(:,2) = rho_trash 450 s2(:,2) = s2_trash ; s(:,2) = SQRT(s2_trash) 451 null_v(:,2) = 0.0_DP 452 END WHERE 453 ! 454 drup = 0.0_DP ; drdw = 0.0_DP 455 dsup = 0.0_DP ; dsdw = 0.0_DP 456 ! 457 drup(:,1) = MIN(1.D-4, 1.D-2*r(:,1)) ; dsup(:,1) = MIN(1.D-4, 1.D-2*s(:,1)) 458 drdw(:,2) = MIN(1.D-4, 1.D-2*r(:,2)) ; dsdw(:,2) = MIN(1.D-4, 1.D-2*s(:,2)) 459 ! 460 ! ... up 461 raux(i1:f1,:) = r+drup ; s2aux(i1:f1,:) = s2 462 raux(i2:f2,:) = r-drup ; s2aux(i2:f2,:) = s2 463 raux(i3:f3,:) = r ; s2aux(i3:f3,:) = (s+dsup)**2 464 raux(i4:f4,:) = r ; s2aux(i4:f4,:) = (s-dsup)**2 465 ! ... down 466 raux(i5:f5,:) = r+drdw ; s2aux(i5:f5,:) = s2 467 raux(i6:f6,:) = r-drdw ; s2aux(i6:f6,:) = s2 468 raux(i7:f7,:) = r ; s2aux(i7:f7,:) = (s+dsdw)**2 469 raux(i8:f8,:) = r ; s2aux(i8:f8,:) = (s-dsdw)**2 470 ! 471 ! 472 CALL gcx_spin( length*8, raux, s2aux, sx, v1x, v2x ) 473 ! 474 ! ... up 475 vrrx(:,1) = 0.5_DP * (v1x(i1:f1,1) - v1x(i2:f2,1)) / drup(:,1) 476 vrsx(:,1) = 0.25_DP * ((v2x(i1:f1,1) - v2x(i2:f2,1)) / drup(:,1) + & 477 (v1x(i3:f3,1) - v1x(i4:f4,1)) / dsup(:,1) / s(:,1)) 478 vssx(:,1) = 0.5_DP * (v2x(i3:f3,1) - v2x(i4:f4,1)) / dsup(:,1) / s(:,1) 479 ! ... down 480 vrrx(:,2) = 0.5_DP * (v1x(i5:f5,2) - v1x(i6:f6,2)) / drdw(:,2) 481 vrsx(:,2) = 0.25_DP * ((v2x(i5:f5,2) - v2x(i6:f6,2)) / drdw(:,2) + & 482 (v1x(i7:f7,2) - v1x(i8:f8,2)) / dsdw(:,2) / s(:,2)) 483 vssx(:,2) = 0.5_DP * (v2x(i7:f7,2) - v2x(i8:f8,2)) / dsdw(:,2) / s(:,2) 484 ! 485 vrrx(:,1) = vrrx(:,1)*null_v(:,1) ; vrrx(:,2) = vrrx(:,2)*null_v(:,2) 486 vrsx(:,1) = vrsx(:,1)*null_v(:,1) ; vrsx(:,2) = vrsx(:,2)*null_v(:,2) 487 vssx(:,1) = vssx(:,1)*null_v(:,1) ; vssx(:,2) = vssx(:,2)*null_v(:,2) 488 ! 489 DEALLOCATE( raux, s2aux ) 490 DEALLOCATE( sx, v1x, v2x ) 491 ! 492 ! ... CORRELATION 493 ! 494 i1 = 1 ; f1 = length !6 blocks: [ rt+dr , s2t , zeta ] 495 i2 = f1+1 ; f2 = 2*length ! [ rt-dr , s2t , zeta ] 496 i3 = f2+1 ; f3 = 3*length ! [ rt , (st+ds)^2 , zeta ] 497 i4 = f3+1 ; f4 = 4*length ! [ rt , (st-ds)^2 , zeta ] 498 i5 = f4+1 ; f5 = 5*length ! [ rt , grho2 , zeta+dz ] 499 i6 = f5+1 ; f6 = 6*length ! [ rt , grho2 , zeta-dz ] 500 ! 501 ALLOCATE( rtaux(length*6), s2taux(length*6), zetaux(length*6) ) 502 ALLOCATE( v1c(length*6,2), v2c(length*6), sc(length*6) ) 503 ! 504 rt(:) = r_in(:,1) + r_in(:,2) 505 ! 506 null_v = 1.0_DP 507 ! 508 WHERE (rt > eps) 509 zeta = (r_in(:,1) - r_in(:,2)) / rt(:) 510 ELSEWHERE 511 zeta = zeta_trash 512 null_v(:,1) = 0.0_DP 513 END WHERE 514 ! 515 s2t = (g_in(:,1,1) + g_in(:,1,2))**2 + & 516 (g_in(:,2,1) + g_in(:,2,2))**2 + & 517 (g_in(:,3,1) + g_in(:,3,2))**2 518 st = SQRT(s2t) 519 ! 520 WHERE (rt<eps .OR. ABS(zeta)>1._DP .OR. st<eps) 521 rt(:) = rho_trash 522 s2t(:) = s2_trash ; st = SQRT(s2_trash) 523 zeta(:) = zeta_trash 524 null_v(:,1) = 0.0_DP 525 END WHERE 526 ! 527 dr = MIN(1.D-4, 1.D-2 * rt) 528 ds = MIN(1.D-4, 1.D-2 * st) 529 !dz = MIN(1.D-4, 1.D-2 * ABS(zeta) ) 530 dz = 1.D-6 531 ! 532 ! ... If zeta is too close to +-1 the derivative is evaluated at a 533 ! slightly smaller value. 534 zeta = SIGN( MIN(ABS(zeta), (1.0_DP - 2.0_DP*dz)), zeta ) 535 ! 536 rtaux(i1:f1) = rt+dr ; s2taux(i1:f1) = s2t ; zetaux(i1:f1) = zeta 537 rtaux(i2:f2) = rt-dr ; s2taux(i2:f2) = s2t ; zetaux(i2:f2) = zeta 538 rtaux(i3:f3) = rt ; s2taux(i3:f3) = (st+ds)**2 ; zetaux(i3:f3) = zeta 539 rtaux(i4:f4) = rt ; s2taux(i4:f4) = (st-ds)**2 ; zetaux(i4:f4) = zeta 540 rtaux(i5:f5) = rt ; s2taux(i5:f5) = s2t ; zetaux(i5:f5) = zeta+dz 541 rtaux(i6:f6) = rt ; s2taux(i6:f6) = s2t ; zetaux(i6:f6) = zeta-dz 542 ! 543 CALL gcc_spin( length*6, rtaux, zetaux, s2taux, sc, v1c, v2c ) 544 ! 545 vrrc(:,1) = 0.5_DP * (v1c(i1:f1,1) - v1c(i2:f2,1)) / dr * null_v(:,1) 546 vrrc(:,2) = 0.5_DP * (v1c(i1:f1,2) - v1c(i2:f2,2)) / dr * null_v(:,1) 547 vrsc(:,1) = 0.5_DP * (v1c(i3:f3,1) - v1c(i4:f4,1)) / ds/st * null_v(:,1) 548 vrsc(:,2) = 0.5_DP * (v1c(i3:f3,2) - v1c(i4:f4,2)) / ds/st * null_v(:,1) 549 vssc(:) = 0.5_DP * (v2c(i3:f3) - v2c(i4:f4) ) / ds/st * null_v(:,1) 550 vrzc(:,1) = 0.5_DP * (v1c(i5:f5,1) - v1c(i6:f6,1)) / dz * null_v(:,1) 551 vrzc(:,2) = 0.5_DP * (v1c(i5:f5,2) - v1c(i6:f6,2)) / dz * null_v(:,1) 552 ! 553 RETURN 554 ! 555END SUBROUTINE dgcxc_spin 556! 557! 558!----------------------------------------------------------------------- 559SUBROUTINE d3gcxc( r, s2, vrrrx, vsrrx, vssrx, vsssx, & 560 vrrrc, vsrrc, vssrc, vsssc ) 561 !----------------------------------------------------------------------- 562 ! wat20101006: Calculates all derivatives of the exchange (x) and 563 ! correlation (c) potential in third order. 564 ! of the Exc. 565 ! 566 ! input: r = rho, s2=|\nabla rho|^2 567 ! definition: E_xc = \int ( f_x(r,s2) + f_c(r,s2) ) dr 568 ! output: vrrrx = d^3(f_x)/d(r)^3 569 ! vsrrx = d^3(f_x)/d(|\nabla r|)d(r)^2 / |\nabla r| 570 ! vssrx = d/d(|\nabla r|) [ & 571 ! d^2(f_x)/d(|\nabla r|)d(r) / |\nabla r| ] & 572 ! / |\nabla r| 573 ! vsssx = d/d(|\nabla r|) [ & 574 ! d/d(|\nabla r|) [ & 575 ! d(f_x)/d(|\nabla r|) / |\nabla r| ] & 576 ! / |\nabla r| ] & 577 ! / |\nabla r| 578 ! same for (c) 579 ! 580 USE kinds, ONLY : DP 581 ! 582 IMPLICIT NONE 583 ! 584 REAL(DP) :: r, s2 585 REAL(DP) :: vrrrx, vsrrx, vssrx, vsssx, vrrrc, vsrrc, vssrc, vsssc 586 ! 587 ! ... local variables 588 ! 589 REAL(DP) :: dr, s, ds 590 REAL(DP), DIMENSION(4) :: raux, s2aux 591 REAL(DP), DIMENSION(4) :: vrrx_rs, vsrx_rs, vssx_rs, vrrc_rs, & 592 vsrc_rs, vssc_rs 593 ! 594 s = SQRT(s2) 595 dr = MIN(1.d-4, 1.d-2 * r) 596 ds = MIN(1.d-4, 1.d-2 * s) 597 ! 598 raux(1) = r+dr ; s2aux(1) = s2 !4 blocks: [ rho+dr , grho2 ] 599 raux(2) = r-dr ; s2aux(2) = s2 ! [ rho-dr , grho2 ] 600 raux(3) = r ; s2aux(3) = (s+ds)**2 ! [ rho , (grho+ds)^2 ] 601 raux(4) = r ; s2aux(4) = (s-ds)**2 ! [ rho , (grho-ds)^2 ] 602 ! 603 CALL dgcxc_unpol( 4, raux, s2aux, vrrx_rs, vsrx_rs, vssx_rs, vrrc_rs, vsrc_rs, vssc_rs ) 604 ! 605 vrrrx = 0.5d0 * (vrrx_rs(1) - vrrx_rs(2)) / dr 606 vsrrx = 0.25d0 * ((vsrx_rs(1) - vsrx_rs(2)) / dr & 607 +(vrrx_rs(3) - vrrx_rs(4)) / ds / s) 608 vssrx = 0.25d0 * ((vssx_rs(1) - vssx_rs(2)) / dr & 609 +(vsrx_rs(3) - vsrx_rs(4)) / ds / s) 610 vsssx = 0.5d0 * (vssx_rs(3) - vssx_rs(4)) / ds / s 611 ! 612 vrrrc = 0.5d0 * (vrrc_rs(1) - vrrc_rs(2)) / dr 613 vsrrc = 0.25d0 * ((vsrc_rs(1) - vsrc_rs(2)) / dr & 614 +(vrrc_rs(3) - vrrc_rs(4)) / ds / s) 615 vssrc = 0.25d0 * ((vssc_rs(1) - vssc_rs(2)) / dr & 616 +(vsrc_rs(3) - vsrc_rs(4)) / ds / s) 617 vsssc = 0.5d0 * (vssc_rs(3) - vssc_rs(4)) / ds / s 618 ! 619 RETURN 620 ! 621END SUBROUTINE d3gcxc 622