1! 2! Copyright (C) 2018 Quantum ESPRESSO Foundation 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! Various routines computing gradient and similar quantities via FFT 10!-------------------------------------------------------------------- 11! FIXME: there is a dependency upon "cell_base" via variable tpiba 12! (2\pi/a) that maybe should be taken out from here? 13!-------------------------------------------------------------------- 14SUBROUTINE external_gradient( a, grada ) 15!-------------------------------------------------------------------- 16 ! 17 ! Interface for computing gradients of a real function in real space, 18 ! to be called by an external module 19 ! 20 USE kinds, ONLY : DP 21 USE fft_base, ONLY : dfftp 22 USE gvect, ONLY : g 23 ! 24 IMPLICIT NONE 25 ! 26 REAL( DP ), INTENT(IN) :: a( dfftp%nnr ) 27 REAL( DP ), INTENT(OUT) :: grada( 3, dfftp%nnr ) 28 29! A in real space, grad(A) in real space 30 CALL fft_gradient_r2r( dfftp, a, g, grada ) 31 32 RETURN 33 34END SUBROUTINE external_gradient 35!---------------------------------------------------------------------------- 36SUBROUTINE fft_gradient_r2r( dfft, a, g, ga ) 37 !---------------------------------------------------------------------------- 38 ! 39 ! ... Calculates ga = \grad a 40 ! ... input : dfft FFT descriptor 41 ! ... a(:) a real function on the real-space FFT grid 42 ! ... g(3,:) G-vectors, in 2\pi/a units 43 ! ... output: ga(3,:) \grad a, real, on the real-space FFT grid 44 ! 45 USE kinds, ONLY : DP 46 USE cell_base, ONLY : tpiba 47 USE fft_interfaces,ONLY : fwfft, invfft 48 USE fft_types, ONLY : fft_type_descriptor 49 ! 50 IMPLICIT NONE 51 ! 52 TYPE(fft_type_descriptor),INTENT(IN) :: dfft 53 REAL(DP), INTENT(IN) :: a(dfft%nnr), g(3,dfft%ngm) 54 REAL(DP), INTENT(OUT) :: ga(3,dfft%nnr) 55 ! 56 INTEGER :: ipol 57 COMPLEX(DP), ALLOCATABLE :: aux(:), gaux(:) 58 ! 59 ALLOCATE( aux( dfft%nnr ) ) 60 ALLOCATE( gaux( dfft%nnr ) ) 61 ! 62 aux = CMPLX( a(:), 0.0_dp, kind=DP) 63 ! 64 ! ... bring a(r) to G-space, a(G) ... 65 ! 66 CALL fwfft ('Rho', aux, dfft) 67 ! 68 ! ... multiply by (iG) to get (\grad_ipol a)(G) ... 69 ! 70 DO ipol = 1, 3 71 ! 72 gaux(:) = (0.0_dp, 0.0_dp) 73 ! 74 gaux(dfft%nl(:)) = g(ipol,:) * CMPLX( -AIMAG( aux(dfft%nl(:)) ), & 75 REAL( aux(dfft%nl(:)) ), kind=DP) 76 ! 77 IF ( dfft%lgamma ) THEN 78 ! 79 gaux(dfft%nlm(:)) = CMPLX( REAL( gaux(dfft%nl(:)) ), & 80 -AIMAG( gaux(dfft%nl(:)) ), kind=DP) 81 ! 82 END IF 83 ! 84 ! ... bring back to R-space, (\grad_ipol a)(r) ... 85 ! 86 CALL invfft ('Rho', gaux, dfft) 87 ! 88 ! ...and add the factor 2\pi/a missing in the definition of G 89 ! 90 ga(ipol,:) = tpiba * DBLE( gaux(:) ) 91 ! 92 END DO 93 ! 94 DEALLOCATE( gaux ) 95 DEALLOCATE( aux ) 96 ! 97 RETURN 98 ! 99END SUBROUTINE fft_gradient_r2r 100! 101!-------------------------------------------------------------------- 102SUBROUTINE fft_qgradient (dfft, a, xq, g, ga) 103 !-------------------------------------------------------------------- 104 ! 105 ! Like fft_gradient_r2r, for complex arrays having a e^{iqr} behavior 106 ! ... input : dfft FFT descriptor 107 ! ... a(:) a complex function on the real-space FFT grid 108 ! ... xq(3) q-vector, in 2\pi/a units 109 ! ... g(3,:) G-vectors, in 2\pi/a units 110 ! ... output: ga(3,:) \grad a, complex, on the real-space FFT grid 111 ! 112 USE kinds, ONLY: dp 113 USE cell_base, ONLY: tpiba 114 USE fft_types, ONLY : fft_type_descriptor 115 USE fft_interfaces, ONLY: fwfft, invfft 116 ! 117 IMPLICIT NONE 118 ! 119 TYPE(fft_type_descriptor),INTENT(IN) :: dfft 120 ! 121 COMPLEX(DP), INTENT(IN) :: a(dfft%nnr) 122 REAL(DP), INTENT(IN):: xq(3), g(3,dfft%ngm) 123 COMPLEX(DP), INTENT(OUT) :: ga(3,dfft%nnr) 124 125 INTEGER :: n, ipol 126 COMPLEX(DP), ALLOCATABLE :: aux(:), gaux(:) 127 128 ALLOCATE (gaux(dfft%nnr)) 129 ALLOCATE (aux (dfft%nnr)) 130 131 ! bring a(r) to G-space, a(G) ... 132 aux (:) = a(:) 133 134 CALL fwfft ('Rho', aux, dfft) 135 ! multiply by i(q+G) to get (\grad_ipol a)(q+G) ... 136 DO ipol = 1, 3 137 gaux (:) = (0.0_dp, 0.0_dp) 138 DO n = 1, dfft%ngm 139 gaux(dfft%nl(n)) = CMPLX( 0.0_dp, xq (ipol) + g(ipol,n), kind=DP ) * & 140 aux (dfft%nl(n)) 141 IF ( dfft%lgamma ) gaux(dfft%nlm(n)) = CONJG( gaux (dfft%nl(n)) ) 142 END DO 143 ! bring back to R-space, (\grad_ipol a)(r) ... 144 145 CALL invfft ('Rho', gaux, dfft) 146 ! ...and add the factor 2\pi/a missing in the definition of q+G 147 DO n = 1, dfft%nnr 148 ga (ipol,n) = gaux (n) * tpiba 149 END DO 150 END DO 151 152 DEALLOCATE (aux) 153 DEALLOCATE (gaux) 154 155 RETURN 156 157END SUBROUTINE fft_qgradient 158! 159!---------------------------------------------------------------------------- 160SUBROUTINE fft_gradient_g2r( dfft, a, g, ga ) 161 !---------------------------------------------------------------------------- 162 ! 163 ! ... Calculates ga = \grad a - like fft_gradient with a(G) instead of a(r) 164 ! ... input : dfft FFT descriptor 165 ! ... a(:) a(G), a complex function in G-space 166 ! ... g(3,:) G-vectors, in 2\pi/a units 167 ! ... output: ga(3,:) \grad a, real, on the real-space FFT grid 168 ! 169 USE cell_base, ONLY : tpiba 170 USE kinds, ONLY : DP 171 USE fft_interfaces,ONLY : invfft 172 USE fft_types, ONLY : fft_type_descriptor 173 ! 174 IMPLICIT NONE 175 ! 176 TYPE(fft_type_descriptor),INTENT(IN) :: dfft 177 COMPLEX(DP), INTENT(IN) :: a(dfft%ngm) 178 REAL(DP), INTENT(IN) :: g(3,dfft%ngm) 179 REAL(DP), INTENT(OUT) :: ga(3,dfft%nnr) 180 ! 181 INTEGER :: ipol, n 182 COMPLEX(DP), ALLOCATABLE :: gaux(:) 183 ! 184 ! 185 ALLOCATE( gaux( dfft%nnr ) ) 186 ga(:,:) = 0.D0 187 ! 188 IF ( dfft%lgamma) THEN 189 ! 190 ! ... Gamma tricks: perform 2 FFT's in a single shot 191 ! x and y 192 ipol = 1 193 gaux(:) = (0.0_dp,0.0_dp) 194 ! 195 ! ... multiply a(G) by iG to get the gradient in real space 196 ! 197 DO n = 1, dfft%ngm 198 gaux(dfft%nl (n)) = CMPLX( 0.0_dp, g(ipol, n), kind=DP )* a(n) - & 199 g(ipol+1,n) * a(n) 200 gaux(dfft%nlm(n)) = CMPLX( 0.0_dp,-g(ipol, n), kind=DP )*CONJG(a(n)) +& 201 g(ipol+1,n) * CONJG(a(n)) 202 ENDDO 203 ! 204 ! ... bring back to R-space, (\grad_ipol a)(r) ... 205 ! 206 CALL invfft ('Rho', gaux, dfft) 207 ! 208 ! ... bring back to R-space, (\grad_ipol a)(r) 209 ! ... add the factor 2\pi/a missing in the definition of q+G 210 ! 211 DO n = 1, dfft%nnr 212 ga (ipol , n) = REAL( gaux(n) ) * tpiba 213 ga (ipol+1, n) = AIMAG( gaux(n) ) * tpiba 214 ENDDO 215 ! z 216 ipol = 3 217 gaux(:) = (0.0_dp,0.0_dp) 218 ! 219 ! ... multiply a(G) by iG to get the gradient in real space 220 ! 221 gaux(dfft%nl (:)) = g(ipol,:) * CMPLX( -AIMAG(a(:)), REAL(a(:)), kind=DP) 222 gaux(dfft%nlm(:)) = CONJG( gaux(dfft%nl(:)) ) 223 ! 224 ! ... bring back to R-space, (\grad_ipol a)(r) ... 225 ! 226 CALL invfft ('Rho', gaux, dfft) 227 ! 228 ! ...and add the factor 2\pi/a missing in the definition of G 229 ! 230 ga(ipol,:) = tpiba * REAL( gaux(:) ) 231 ! 232 ELSE 233 ! 234 DO ipol = 1, 3 235 ! 236 gaux(:) = (0.0_dp,0.0_dp) 237 ! 238 ! ... multiply a(G) by iG to get the gradient in real space 239 ! 240 gaux(dfft%nl(:)) = g(ipol,:) * CMPLX( -AIMAG(a(:)), REAL(a(:)), kind=DP) 241 ! 242 ! ... bring back to R-space, (\grad_ipol a)(r) ... 243 ! 244 CALL invfft ('Rho', gaux, dfft) 245 ! 246 ! ...and add the factor 2\pi/a missing in the definition of G 247 ! 248 ga(ipol,:) = tpiba * REAL( gaux(:) ) 249 ! 250 END DO 251 ! 252 END IF 253 ! 254 DEALLOCATE( gaux ) 255 ! 256 RETURN 257 ! 258END SUBROUTINE fft_gradient_g2r 259 260!---------------------------------------------------------------------------- 261SUBROUTINE fft_graddot( dfft, a, g, da ) 262 !---------------------------------------------------------------------------- 263 ! 264 ! ... Calculates da = \sum_i \grad_i a_i in R-space 265 ! ... input : dfft FFT descriptor 266 ! ... a(3,:) a real function on the real-space FFT grid 267 ! ... g(3,:) G-vectors, in 2\pi/a units 268 ! ... output: ga(:) \sum_i \grad_i a_i, real, on the real-space FFT grid 269 ! 270 USE cell_base, ONLY : tpiba 271 USE kinds, ONLY : DP 272 USE fft_interfaces,ONLY : fwfft, invfft 273 USE fft_types, ONLY : fft_type_descriptor 274 ! 275 IMPLICIT NONE 276 ! 277 TYPE(fft_type_descriptor),INTENT(IN) :: dfft 278 REAL(DP), INTENT(IN) :: a(3,dfft%nnr), g(3,dfft%ngm) 279 REAL(DP), INTENT(OUT) :: da(dfft%nnr) 280 ! 281 INTEGER :: n, ipol 282 COMPLEX(DP), ALLOCATABLE :: aux(:), gaux(:) 283 COMPLEX(DP) :: fp, fm, aux1, aux2 284 ! 285 ALLOCATE( aux(dfft%nnr) ) 286 ALLOCATE( gaux(dfft%nnr) ) 287 ! 288 gaux(:) = (0.0_dp,0.0_dp) 289 ! 290 IF ( dfft%lgamma ) THEN 291 ! 292 ! Gamma tricks: perform 2 FFT's in a single shot 293 ! x and y 294 ipol = 1 295 aux(:) = CMPLX( a(ipol,:), a(ipol+1,:), kind=DP) 296 ! 297 ! ... bring a(ipol,r) to G-space, a(G) ... 298 ! 299 CALL fwfft ('Rho', aux, dfft) 300 ! 301 ! ... multiply by iG to get the gradient in G-space 302 ! 303 DO n = 1, dfft%ngm 304 ! 305 fp = (aux(dfft%nl(n)) + aux (dfft%nlm(n)))*0.5_dp 306 fm = (aux(dfft%nl(n)) - aux (dfft%nlm(n)))*0.5_dp 307 aux1 = CMPLX( REAL(fp), AIMAG(fm), kind=DP) 308 aux2 = CMPLX(AIMAG(fp), -REAL(fm), kind=DP) 309 gaux (dfft%nl(n)) = & 310 CMPLX(0.0_dp, g(ipol ,n),kind=DP) * aux1 + & 311 CMPLX(0.0_dp, g(ipol+1,n),kind=DP) * aux2 312 ENDDO 313 ! z 314 ipol = 3 315 aux(:) = CMPLX( a(ipol,:), 0.0_dp, kind=DP) 316 ! 317 ! ... bring a(ipol,r) to G-space, a(G) ... 318 ! 319 CALL fwfft ('Rho', aux, dfft) 320 ! 321 ! ... multiply by iG to get the gradient in G-space 322 ! ... fill both gaux(G) and gaux(-G) = gaux*(G) 323 ! 324 DO n = 1, dfft%ngm 325 gaux(dfft%nl(n)) = gaux(dfft%nl(n)) + g(ipol,n) * & 326 CMPLX( -AIMAG( aux(dfft%nl(n)) ), & 327 REAL( aux(dfft%nl(n)) ), kind=DP) 328 gaux(dfft%nlm(n)) = CONJG( gaux(dfft%nl(n)) ) 329 END DO 330 ! 331 ELSE 332 ! 333 DO ipol = 1, 3 334 ! 335 aux = CMPLX( a(ipol,:), 0.0_dp, kind=DP) 336 ! 337 ! ... bring a(ipol,r) to G-space, a(G) ... 338 ! 339 CALL fwfft ('Rho', aux, dfft) 340 ! 341 ! ... multiply by iG to get the gradient in G-space 342 ! 343 DO n = 1, dfft%ngm 344 gaux(dfft%nl(n)) = gaux(dfft%nl(n)) + g(ipol,n) * & 345 CMPLX( -AIMAG( aux(dfft%nl(n)) ), & 346 REAL( aux(dfft%nl(n)) ), kind=DP) 347 END DO 348 ! 349 END DO 350 ! 351 END IF 352 ! 353 ! ... bring back to R-space, (\grad_ipol a)(r) ... 354 ! 355 CALL invfft ('Rho', gaux, dfft) 356 ! 357 ! ... add the factor 2\pi/a missing in the definition of G and sum 358 ! 359 da(:) = tpiba * REAL( gaux(:) ) 360 ! 361 DEALLOCATE( aux, gaux ) 362 ! 363 RETURN 364 ! 365END SUBROUTINE fft_graddot 366 367!-------------------------------------------------------------------- 368SUBROUTINE fft_qgraddot ( dfft, a, xq, g, da) 369 !-------------------------------------------------------------------- 370 ! 371 ! Like fft_graddot, for complex arrays having a e^{iqr} dependency 372 ! ... input : dfft FFT descriptor 373 ! ... a(3,:) a complex function on the real-space FFT grid 374 ! ... xq(3) q-vector, in 2\pi/a units 375 ! ... g(3,:) G-vectors, in 2\pi/a units 376 ! ... output: ga(:) \sum_i \grad_i a_i, complex, on the real-space FFT grid 377 ! 378 USE kinds, ONLY : DP 379 USE cell_base, ONLY : tpiba 380 USE fft_interfaces, ONLY : fwfft, invfft 381 USE fft_types, ONLY : fft_type_descriptor 382 ! 383 IMPLICIT NONE 384 ! 385 TYPE(fft_type_descriptor),INTENT(IN) :: dfft 386 COMPLEX(DP), INTENT(IN) :: a(3,dfft%nnr) 387 REAL(DP), INTENT(IN) :: xq(3), g(3,dfft%ngm) 388 COMPLEX(DP), INTENT(OUT) :: da(dfft%nnr) 389 390 INTEGER :: n, ipol 391 COMPLEX(DP), allocatable :: aux (:) 392 393 ALLOCATE (aux (dfft%nnr)) 394 da(:) = (0.0_dp, 0.0_dp) 395 DO ipol = 1, 3 396 ! copy a(ipol,r) to a complex array... 397 DO n = 1, dfft%nnr 398 aux (n) = a (ipol, n) 399 END DO 400 ! bring a(ipol,r) to G-space, a(G) ... 401 CALL fwfft ('Rho', aux, dfft) 402 ! multiply by i(q+G) to get (\grad_ipol a)(q+G) ... 403 DO n = 1, dfft%ngm 404 da (dfft%nl(n)) = da (dfft%nl(n)) + & 405 CMPLX(0.0_dp, xq (ipol) + g (ipol, n),kind=DP) * aux(dfft%nl(n)) 406 END DO 407 END DO 408 IF ( dfft%lgamma ) THEN 409 DO n = 1, dfft%ngm 410 da (dfft%nlm(n)) = CONJG( da (dfft%nl(n)) ) 411 END DO 412 END IF 413 ! bring back to R-space, (\grad_ipol a)(r) ... 414 CALL invfft ('Rho', da, dfft) 415 ! ...add the factor 2\pi/a missing in the definition of q+G 416 da (:) = da (:) * tpiba 417 418 DEALLOCATE(aux) 419 420 RETURN 421 422END SUBROUTINE fft_qgraddot 423 424!-------------------------------------------------------------------- 425! Routines computing laplacian via FFT 426!-------------------------------------------------------------------- 427 428!-------------------------------------------------------------------- 429SUBROUTINE external_laplacian( a, lapla ) 430!-------------------------------------------------------------------- 431 ! 432 ! Interface for computing laplacian in real space, to be called by 433 ! an external module 434 ! 435 USE kinds, ONLY : DP 436 USE fft_base, ONLY : dfftp 437 USE gvect, ONLY : gg 438 ! 439 IMPLICIT NONE 440 ! 441 REAL( DP ), INTENT(IN) :: a( dfftp%nnr ) 442 REAL( DP ), INTENT(OUT) :: lapla( dfftp%nnr ) 443 444! A in real space, lapl(A) in real space 445 CALL fft_laplacian( dfftp, a, gg, lapla ) 446 447 RETURN 448 449END SUBROUTINE external_laplacian 450 451!-------------------------------------------------------------------- 452SUBROUTINE fft_laplacian( dfft, a, gg, lapla ) 453!-------------------------------------------------------------------- 454 ! 455 ! ... Calculates lapla = laplacian(a) 456 ! ... input : dfft FFT descriptor 457 ! ... a(:) a real function on the real-space FFT grid 458 ! ... gg(:) square modules of G-vectors, in (2\pi/a)^2 units 459 ! ... output: lapla(:) \nabla^2 a, real, on the real-space FFT grid 460 ! 461 USE kinds, ONLY : DP 462 USE cell_base, ONLY : tpiba2 463 USE fft_types, ONLY : fft_type_descriptor 464 USE fft_interfaces,ONLY : fwfft, invfft 465 ! 466 IMPLICIT NONE 467 ! 468 TYPE(fft_type_descriptor),INTENT(IN) :: dfft 469 REAL(DP), INTENT(IN) :: a(dfft%nnr), gg(dfft%ngm) 470 REAL(DP), INTENT(OUT) :: lapla(dfft%nnr) 471 ! 472 INTEGER :: ig 473 COMPLEX(DP), ALLOCATABLE :: aux(:), laux(:) 474 ! 475 ! 476 ALLOCATE( aux( dfft%nnr ) ) 477 ALLOCATE( laux( dfft%nnr ) ) 478 ! 479 aux = CMPLX( a(:), 0.0_dp, kind=DP) 480 ! 481 ! ... bring a(r) to G-space, a(G) ... 482 ! 483 CALL fwfft ('Rho', aux, dfft) 484 ! 485 ! ... Compute the laplacian 486 ! 487 laux(:) = (0.0_dp, 0.0_dp) 488 ! 489 DO ig = 1, dfft%ngm 490 ! 491 laux(dfft%nl(ig)) = -gg(ig)*aux(dfft%nl(ig)) 492 ! 493 END DO 494 ! 495 IF ( dfft%lgamma ) THEN 496 ! 497 laux(dfft%nlm(:)) = CMPLX( REAL(laux(dfft%nl(:)) ), & 498 -AIMAG(laux(dfft%nl(:)) ), kind=DP) 499 ! 500 ENDIF 501 ! 502 ! ... bring back to R-space, (\lapl a)(r) ... 503 ! 504 CALL invfft ('Rho', laux, dfft) 505 ! 506 ! ... add the missing factor (2\pi/a)^2 in G 507 ! 508 lapla = tpiba2 * REAL( laux ) 509 ! 510 DEALLOCATE( laux ) 511 DEALLOCATE( aux ) 512 ! 513 RETURN 514 ! 515END SUBROUTINE fft_laplacian 516! 517!-------------------------------------------------------------------- 518! Routines computing hessian via FFT 519!-------------------------------------------------------------------- 520! 521!---------------------------------------------------------------------- 522SUBROUTINE fft_hessian_g2r ( dfft, a, g, ha ) 523!---------------------------------------------------------------------- 524 ! 525 ! ... Calculates ha = hessian(a) 526 ! ... input : dfft FFT descriptor 527 ! ... a(:) a real function on the real-space FFT grid 528 ! ... g(3,:) G-vectors, in (2\pi/a)^2 units 529 ! ... output: ha(6,:) hessian(a), real, on the real-space FFT grid 530 ! ... lower-packed matrix indeces 1-6 correspond to: 531 ! ... 1 = xx, 2 = yx, 3 = yy, 4 = zx, 5 = zy, 6 = zz 532 ! 533 USE kinds, ONLY : DP 534 USE cell_base, ONLY : tpiba 535 USE fft_types, ONLY : fft_type_descriptor 536 USE fft_interfaces,ONLY : fwfft, invfft 537 USE fft_helper_subroutines, ONLY: fftx_oned2threed 538 ! 539 IMPLICIT NONE 540 ! 541 TYPE(fft_type_descriptor),INTENT(IN) :: dfft 542 REAL(DP), INTENT(IN) :: g(3,dfft%ngm) 543 COMPLEX(DP), INTENT(IN) :: a(dfft%ngm) 544 REAL(DP), INTENT(OUT) :: ha( 6, dfft%nnr ) 545 ! 546 INTEGER :: ig, ir 547 COMPLEX(DP), ALLOCATABLE :: aux(:), haux(:,:) 548 ! 549 IF ( .NOT. dfft%lgamma ) CALL errore ('fft_hessian_g2r',& 550 'only gamma case is implemented',1) 551 ALLOCATE ( aux(dfft%nnr)) 552 ALLOCATE (haux(dfft%ngm,2)) 553 ! xx, yx 554 DO ig=1,dfft%ngm 555 haux(ig,1) = -tpiba**2*g(1,ig)**2 *a(ig) 556 haux(ig,2) = -tpiba**2*g(1,ig)*g(2,ig)*a(ig) 557 END DO 558 CALL fftx_oned2threed( dfft, aux, haux(:,1), haux(:,2) ) 559 CALL invfft('Rho', aux, dfft) 560 DO ir=1,dfft%nnr 561 ha(1,ir) = DBLE(aux(ir)) 562 ha(2,ir) =AIMAG(aux(ir)) 563 END DO 564 ! yy, zx 565 DO ig=1,dfft%ngm 566 haux(ig,1) = -tpiba**2*g(2,ig)**2 *a(ig) 567 haux(ig,2) = -tpiba**2*g(1,ig)*g(3,ig)*a(ig) 568 END DO 569 CALL fftx_oned2threed( dfft, aux, haux(:,1), haux(:,2) ) 570 CALL invfft('Rho', aux, dfft) 571 DO ir=1,dfft%nnr 572 ha(3,ir) = DBLE(aux(ir)) 573 ha(4,ir) =AIMAG(aux(ir)) 574 END DO 575 ! zy, zz 576 DO ig=1,dfft%ngm 577 haux(ig,1) = -tpiba**2*g(2,ig)*g(3,ig)*a(ig) 578 haux(ig,2) = -tpiba**2*g(3,ig)**2 *a(ig) 579 END DO 580 CALL fftx_oned2threed( dfft, aux, haux(:,1), haux(:,2) ) 581 CALL invfft('Rho', aux, dfft) 582 DO ir=1,dfft%nnr 583 ha(5,ir) = DBLE(aux(ir)) 584 ha(6,ir) =AIMAG(aux(ir)) 585 END DO 586 ! 587 DEALLOCATE(aux) 588 DEALLOCATE(haux) 589 590END SUBROUTINE fft_hessian_g2r 591!-------------------------------------------------------------------- 592SUBROUTINE fft_hessian( dfft, a, g, ga, ha ) 593!-------------------------------------------------------------------- 594 ! 595 ! ... Calculates ga = \grad a and ha = hessian(a) 596 ! ... input : dfft FFT descriptor 597 ! ... a(:) a real function on the real-space FFT grid 598 ! ... g(3,:) G-vectors, in (2\pi/a)^2 units 599 ! ... output: ga(3,:) \grad a, real, on the real-space FFT grid 600 ! ... ha(3,3,:) hessian(a), real, on the real-space FFT grid 601 ! 602 USE kinds, ONLY : DP 603 USE cell_base, ONLY : tpiba 604 USE fft_types, ONLY : fft_type_descriptor 605 USE fft_interfaces,ONLY : fwfft, invfft 606 ! 607 IMPLICIT NONE 608 ! 609 TYPE(fft_type_descriptor),INTENT(IN) :: dfft 610 REAL(DP), INTENT(IN) :: a(dfft%nnr), g(3,dfft%ngm) 611 REAL(DP), INTENT(OUT) :: ga( 3, dfft%nnr ) 612 REAL(DP), INTENT(OUT) :: ha( 3, 3, dfft%nnr ) 613 ! 614 INTEGER :: ipol, jpol 615 COMPLEX(DP), ALLOCATABLE :: aux(:), gaux(:), haux(:) 616 ! 617 ! 618 ALLOCATE( aux( dfft%nnr ) ) 619 ALLOCATE( gaux( dfft%nnr ) ) 620 ALLOCATE( haux( dfft%nnr ) ) 621 ! 622 aux = CMPLX( a(:), 0.0_dp, kind=DP) 623 ! 624 ! ... bring a(r) to G-space, a(G) ... 625 ! 626 CALL fwfft ('Rho', aux, dfft) 627 ! 628 ! ... multiply by (iG) to get (\grad_ipol a)(G) ... 629 ! 630 DO ipol = 1, 3 631 ! 632 gaux(:) = (0.0_dp,0.0_dp) 633 ! 634 gaux(dfft%nl(:)) = g(ipol,:) * CMPLX( -AIMAG( aux(dfft%nl(:)) ), & 635 REAL( aux(dfft%nl(:)) ), kind=DP ) 636 ! 637 IF ( dfft%lgamma ) THEN 638 ! 639 gaux(dfft%nlm(:)) = CMPLX( REAL( gaux(dfft%nl(:)) ), & 640 -AIMAG( gaux(dfft%nl(:)) ), kind=DP) 641 ! 642 END IF 643 ! 644 ! ... bring back to R-space, (\grad_ipol a)(r) ... 645 ! 646 CALL invfft ('Rho', gaux, dfft) 647 ! 648 ! ...and add the factor 2\pi/a missing in the definition of G 649 ! 650 ga(ipol,:) = tpiba * REAL( gaux(:) ) 651 ! 652 ! ... compute the second derivatives 653 ! 654 DO jpol = 1, ipol 655 ! 656 haux(:) = (0.0_dp,0.0_dp) 657 ! 658 haux(dfft%nl(:)) = - g(ipol,:) * g(jpol,:) * & 659 CMPLX( REAL( aux(dfft%nl(:)) ), & 660 AIMAG( aux(dfft%nl(:)) ), kind=DP) 661 ! 662 IF ( dfft%lgamma ) THEN 663 ! 664 haux(dfft%nlm(:)) = CMPLX( REAL( haux(dfft%nl(:)) ), & 665 -AIMAG( haux(dfft%nl(:)) ), kind=DP) 666 ! 667 END IF 668 ! 669 ! ... bring back to R-space, (\grad_ipol a)(r) ... 670 ! 671 CALL invfft ('Rho', haux, dfft) 672 ! 673 ! ...and add the factor 2\pi/a missing in the definition of G 674 ! 675 ha(ipol, jpol, :) = tpiba * tpiba * REAL( haux(:) ) 676 ! 677 ha(jpol, ipol, :) = ha(ipol, jpol, :) 678 ! 679 END DO 680 ! 681 END DO 682 ! 683 DEALLOCATE( haux ) 684 DEALLOCATE( gaux ) 685 DEALLOCATE( aux ) 686 ! 687 RETURN 688 ! 689END SUBROUTINE fft_hessian 690!-------------------------------------------------------------------- 691SUBROUTINE external_hessian( a, grada, hessa ) 692!-------------------------------------------------------------------- 693 ! 694 ! Interface for computing hessian in real space, to be called by 695 ! an external module 696 ! 697 USE kinds, ONLY : DP 698 USE fft_base, ONLY : dfftp 699 USE gvect, ONLY : g 700 ! 701 IMPLICIT NONE 702 ! 703 REAL( DP ), INTENT(IN) :: a( dfftp%nnr ) 704 REAL( DP ), INTENT(OUT) :: grada( 3, dfftp%nnr ) 705 REAL( DP ), INTENT(OUT) :: hessa( 3, 3, dfftp%nnr ) 706 707! A in real space, grad(A) and hess(A) in real space 708 CALL fft_hessian( dfftp, a, g, grada, hessa ) 709 710 RETURN 711 712END SUBROUTINE external_hessian 713 714!-------------------------------------------------------------------- 715