1!
2! Copyright (C) 2001-2008 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 drhodvus (irr, imode0, dvscfin, npe)
11  !-----------------------------------------------------------------------
12  !
13  !    This subroutine calculates the term of the dynamical matrix
14  !    which comes from the interaction of the change of the self
15  !    consistent potential with the change of the charge density
16  !    at fixed wavefunctions.
17  !    See Eq.B36 of PRB 64, 235118 (2001).
18  !    In the PAW case this part of the dynamical matrix has an additional
19  !    contribution.
20  !
21  !
22  USE kinds,     ONLY : DP
23  USE ions_base, ONLY : nat, ntyp=>nsp, ityp
24  USE cell_base, ONLY : omega
25  USE ions_base, ONLY : nat
26  USE fft_base,  ONLY : dfftp
27  USE uspp,      ONLY : okvan
28  USE io_global, ONLY : stdout
29  USE buffers,   ONLY : get_buffer
30  USE uspp_param, ONLY : upf, nh
31  USE paw_variables, ONLY : okpaw
32  USE noncollin_module, ONLY : nspin_mag
33
34  USE modes,     ONLY : npert, nirr, u
35  USE dynmat,    ONLY : dyn, dyn_rec
36  USE phus,      ONLY : becsumort
37  USE units_ph,  ONLY : iudrhous, lrdrhous
38
39  USE mp_pools,  ONLY : inter_pool_comm
40  USE mp_bands,  ONLY : intra_bgrp_comm
41  USE mp,        ONLY : mp_sum
42
43  USE lrus,      ONLY : int3_paw
44  implicit none
45
46  integer :: irr, imode0, npe
47  ! input: the irreducible representation
48  ! input: starting position of this represe
49  ! input: the number of perturbations
50
51  complex(DP) :: dvscfin (dfftp%nnr, nspin_mag, npe)
52  ! input: the change of V_Hxc
53
54  integer :: ipert, irr1, mode0, mu, is, nu_i, nu_j, nrtot, &
55             ih, jh, ijh, na, nt
56  ! counters
57  ! mode0: starting position of the represention
58  ! nrtot: the total number of mesh points
59
60  complex(DP) :: dyn1 (3 * nat, 3 * nat)
61  ! the dynamical matrix
62  complex(DP), allocatable ::  drhous (:,:)
63  ! the change of the charge
64  complex(DP), external :: zdotc
65
66  if (.not.okvan) then
67     dyn_rec=(0.0_DP,0.0_DP)
68     return
69  endif
70  call start_clock ('drhodvus')
71  allocate (drhous ( dfftp%nnr, nspin_mag))
72  dyn1 (:,:) = (0.d0, 0.d0)
73  nrtot = dfftp%nr1 * dfftp%nr2 * dfftp%nr3
74  mode0 = 0
75  do irr1 = 1, nirr
76     do ipert = 1, npert (irr1)
77        nu_j = mode0 + ipert
78        call get_buffer (drhous,  lrdrhous, iudrhous, nu_j)
79        do mu = 1, npert (irr)
80           nu_i = imode0 + mu
81           ! FIXME : use zgemm instead of zdotc
82           dyn1 (nu_i, nu_j) = dyn1 (nu_i, nu_j) + &
83                   zdotc (dfftp%nnr*nspin_mag,dvscfin(1,1,mu),1,drhous, 1) &
84                   * omega / DBLE (nrtot)
85        enddo
86     enddo
87     mode0 = mode0 + npert (irr1)
88  enddo
89  deallocate (drhous)
90  !
91  ! collect contributions from all pools (sum over k-points)
92  !
93  call mp_sum ( dyn1, inter_pool_comm )
94  call mp_sum ( dyn1, intra_bgrp_comm )
95!
96!  PAW contribution: this part of the dynamical matrix is present only
97!  with PAW. PAW and US dynamical matrices differ only at this point.
98!
99  IF (okpaw) THEN
100     mode0 = 0
101     do irr1 = 1, nirr
102        do ipert = 1, npert (irr1)
103           nu_j = mode0 + ipert
104           do mu = 1, npert (irr)
105              nu_i = imode0 + mu
106              do nt=1,ntyp
107                 if (upf(nt)%tpawp) then
108                    ijh=0
109                    do ih=1,nh(nt)
110                       do jh=ih,nh(nt)
111                          ijh=ijh+1
112                          do na=1,nat
113                             if (ityp(na)==nt) then
114                                do is = 1, nspin_mag
115                                   dyn1(nu_i,nu_j)=dyn1(nu_i,nu_j)+ &
116                                       CONJG(int3_paw(ih,jh,na,is,mu))* &
117                                       becsumort(ijh,na,is,nu_j)
118                                enddo
119                             endif
120                          enddo
121                       enddo
122                    enddo
123                 endif
124              enddo
125           enddo
126        enddo
127        mode0 = mode0 + npert (irr1)
128      enddo
129   endif
130  !       WRITE( stdout,*) 'drhodvus dyn1, dyn'
131  !       call tra_write_matrix('drhodvus dyn1',dyn1,u,nat)
132  !       call tra_write_matrix('drhodvus dyn',dyn,u,nat)
133  !       call stop_ph(.true.)
134  dyn (:,:) = dyn (:,:) + dyn1 (:,:)
135  dyn_rec(:,:) = dyn1(:,:)
136  call stop_clock ('drhodvus')
137  return
138
139end subroutine drhodvus
140