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