1! 2! Copyright (C) 2001-2007 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 force_cc( forcecc ) 11 !---------------------------------------------------------------------- 12 !! Calculates the NLCC contribution to the force. 13 ! 14 USE kinds, ONLY : DP 15 USE constants, ONLY : tpi 16 USE atom, ONLY : rgrid, msh 17 USE uspp_param, ONLY : upf 18 USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau 19 USE cell_base, ONLY : alat, omega, tpiba, tpiba2 20 USE fft_base, ONLY : dfftp 21 USE fft_interfaces, ONLY : fwfft 22 USE gvect, ONLY : ngm, gstart, g, gg, ngl, gl, igtongl 23 USE ener, ONLY : etxc, vtxc 24 USE lsda_mod, ONLY : nspin 25 USE scf, ONLY : rho, rho_core, rhog_core 26 USE control_flags, ONLY : gamma_only 27 USE noncollin_module, ONLY : noncolin 28 USE wavefunctions, ONLY : psic 29 USE mp_bands, ONLY : intra_bgrp_comm 30 USE mp, ONLY : mp_sum 31 ! 32 IMPLICIT NONE 33 ! 34 REAL(DP) :: forcecc(3,nat) 35 !! output: the NLCC forces on atoms 36 ! 37 ! ... local variables 38 ! 39 INTEGER :: ig, ir, nt, na 40 ! counter on polarizations 41 ! counter on G vectors 42 ! counter on FFT grid points 43 ! counter on types of atoms 44 ! counter on atoms 45 REAL(DP), ALLOCATABLE :: vxc(:,:), rhocg(:) 46 ! exchange-correlation potential 47 ! radial fourier transform of rho core 48 REAL(DP) :: arg, fact 49 ! 50 ! 51 forcecc(:,:) = 0.d0 52 ! 53 IF ( ANY( upf(1:ntyp)%nlcc ) ) GOTO 15 54 RETURN 55 ! 5615 CONTINUE 57 IF (gamma_only) THEN 58 fact = 2.d0 59 ELSE 60 fact = 1.d0 61 ENDIF 62 ! 63 ! ... recalculate the exchange-correlation potential 64 ! 65 ALLOCATE( vxc(dfftp%nnr,nspin) ) 66 ! 67 CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, vxc ) 68 ! 69 psic = (0.0_DP,0.0_DP) 70 IF (nspin == 1 .OR. nspin == 4) THEN 71 DO ir = 1, dfftp%nnr 72 psic(ir) = vxc (ir,1) 73 ENDDO 74 ELSE 75 DO ir = 1, dfftp%nnr 76 psic (ir) = 0.5d0 * (vxc (ir, 1) + vxc (ir, 2) ) 77 ENDDO 78 ENDIF 79 ! 80 DEALLOCATE( vxc ) 81 ! 82 CALL fwfft( 'Rho', psic, dfftp ) 83 ! 84 ! ... psic contains now Vxc(G) 85 ! 86 ALLOCATE( rhocg(ngl) ) 87 ! 88 ! ... core correction term: sum on g of omega*ig*exp(-i*r_i*g)*n_core(g)*vxc 89 ! g = 0 term gives no contribution 90 ! 91 DO nt = 1, ntyp 92 IF ( upf(nt)%nlcc ) THEN 93 ! 94 CALL drhoc( ngl, gl, omega, tpiba2, msh(nt), rgrid(nt)%r, & 95 rgrid(nt)%rab, upf(nt)%rho_atc, rhocg ) 96!$omp parallel do private(arg) 97 DO na = 1, nat 98 IF (nt == ityp (na) ) THEN 99 DO ig = gstart, ngm 100 arg = (g(1,ig) * tau(1,na) + g (2, ig) * tau (2, na) & 101 + g(3,ig) * tau(3,na) ) * tpi 102 forcecc (1:3, na) = forcecc(1:3, na) + tpiba * omega * & 103 rhocg(igtongl(ig)) * CONJG(psic(dfftp%nl(ig) ) ) * & 104 CMPLX( SIN(arg), COS(arg), KIND=DP) * g(1:3,ig) * fact 105 ENDDO 106 ENDIF 107 ENDDO 108!$omp end parallel do 109 ENDIF 110 ENDDO 111 ! 112 CALL mp_sum( forcecc, intra_bgrp_comm ) 113 ! 114 DEALLOCATE( rhocg ) 115 ! 116 RETURN 117 ! 118END SUBROUTINE force_cc 119