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