1! 2! Copyright (C) 2010 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! ... original code written by Giovanni Cantele and Paolo Cazzato 9! ... adapted to work in the parallel case by Carlo Sbraccia 10! ... code for the calculation of the vacuum level written by Carlo Sbraccia 11! ... code ported from PW to CP by Federico Zipoli 12! 13SUBROUTINE makov_payne(etot) 14! 15! CP Modules 16 USE kinds, ONLY : DP 17 USE ions_base, ONLY : nat, zv, ityp 18 USE ions_positions, ONLY : tau0 19 USE io_global, ONLY : stdout, ionode, ionode_id 20 USE constants, ONLY : pi, autoev, au_debye 21 USE cp_main_variables, ONLY : rhor 22 USE electrons_base, ONLY : nspin 23 USE cell_base, ONLY : at, bg, omega, alat, ibrav 24 USE parallel_include 25 USE gvecw , ONLY : ngw 26 USE fft_base, ONLY : dfftp 27#if defined __MPI 28 USE mp_global, ONLY : me_bgrp, nproc_bgrp, intra_bgrp_comm 29 USE mp, ONLY : mp_barrier 30 USE mp_world, ONLY : world_comm 31#endif 32! 33IMPLICIT NONE 34INTEGER :: nspecie 35INTEGER :: i,j,k,l,m,n,ip 36INTEGER, ALLOCATABLE, DIMENSION(:) :: zvv 37REAL(DP), ALLOCATABLE, DIMENSION(:,:,:) :: rhof 38REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: r 39REAL(DP) :: h(3,3),volumetto,a(3,3) 40REAL(DP) :: usunx,usuny,usunz,R0(3),qq,aa,bb 41REAL(DP) :: charge, charge_ion, charge_el 42REAL(DP) :: dipole(3), dipole_ion(3), dipole_el(3) 43REAL(DP) :: quadrupole, quadrupole_ion, quadrupole_el 44REAL(DP) :: corr1, corr2, etot 45REAL(DP), ALLOCATABLE, DIMENSION(:) :: rgx,rgy,rgz 46INTEGER :: ir, is 47INTEGER :: ierr 48#if defined __MPI 49INTEGER :: proc 50INTEGER, ALLOCATABLE:: displs(:), recvcount(:) 51#endif 52REAL(KIND=DP), ALLOCATABLE:: rhodist1(:) 53REAL(KIND=DP), ALLOCATABLE:: rhodist2(:) 54! 55 IF(ibrav.NE.1)THEN 56 WRITE(*,*)" WARNING Makov-Payne implemented in CP only when ibrav=1 " 57 RETURN 58 ENDIF 59! 60 usunx=1.0D0/DBLE(dfftp%nr1x) 61 usuny=1.0D0/DBLE(dfftp%nr2x) 62 usunz=1.0D0/DBLE(dfftp%nr3x) 63 ALLOCATE ( r(nat,3),rhof(dfftp%nr1x,dfftp%nr2x,dfftp%nr3x),& 64 & rgx(dfftp%nr1x),rgy(dfftp%nr2x),rgz(dfftp%nr3x),zvv(nat) ) 65 ! 66 DO i=1,nat 67 zvv(i)=zv(ityp(i)) 68 DO j=1,3 69 r(i,j)=tau0(j,i) 70 ENDDO 71 ENDDO 72 ! 73 ip=0 74 rhof=0.0D0 75! 76!-------------------------------------------------------------------- 77ALLOCATE(rhodist1(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x)) 78IF (nspin.EQ.2) ALLOCATE(rhodist2(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x)) 79#if defined __MPI 80ALLOCATE( displs( nproc_bgrp ), recvcount( nproc_bgrp ) ) 81! 82 do proc=1,nproc_bgrp 83 recvcount(proc) = dfftp%nnp * ( dfftp%nr3p(proc) ) 84 if (proc.eq.1) then 85 displs(proc)=0 86 else 87 displs(proc)=displs(proc-1) + recvcount(proc-1) 88 end if 89 end do 90! 91! gather the charge density on the first node 92! 93 call mp_barrier( world_comm ) 94 call mpi_gatherv( rhor(1,1), recvcount(me_bgrp+1), MPI_DOUBLE_PRECISION,& 95 & rhodist1,recvcount, displs, MPI_DOUBLE_PRECISION,& 96 & ionode_id, intra_bgrp_comm, ierr) 97 call errore('mpi_gatherv','ierr<>0',ierr) 98! 99IF(nspin .eq. 2)THEN 100 call mp_barrier( world_comm ) 101 call mpi_gatherv( rhor(1,2), recvcount(me_bgrp+1), MPI_DOUBLE_PRECISION, & 102 & rhodist2,recvcount, displs, MPI_DOUBLE_PRECISION, & 103 & ionode_id, intra_bgrp_comm, ierr) 104 call errore('mpi_gatherv','ierr<>0',ierr) 105ENDIF 106#else 107 rhodist1=rhor(:,1) 108 IF(nspin .eq. 2) rhodist2=rhor(:,2) 109#endif 110! 111#if defined __MPI 112IF ( ionode ) THEN 113#endif 114 DO k = 1, dfftp%nr3x 115 DO j = 1, dfftp%nr2x 116 DO i = 1, dfftp%nr1x 117 ip=ip+1 118 IF (nspin == 1 )rhof(i,j,k)=rhodist1(ip) 119 IF (nspin == 2 )rhof(i,j,k)=rhodist1(ip)+rhodist2(ip) 120 ENDDO 121 ENDDO 122 ENDDO 123 ip=0 124 DO i=1,dfftp%nr1x 125 rgx(i)=DBLE(i-1)*usunx*alat 126 ENDDO 127 DO i=1,dfftp%nr2x 128 rgy(i)=DBLE(i-1)*usuny*alat 129 ENDDO 130 DO i=1,dfftp%nr3x 131 rgz(i)=DBLE(i-1)*usunz*alat 132 ENDDO 133 ! 134 !---------------------------------------------------------- 135 ! 136 ! center of charge of the ions 137 ! 138 R0=0.0D0 139 DO i=1,nat 140 R0(1)=R0(1)+zvv(i)*r(i,1) 141 R0(2)=R0(2)+zvv(i)*r(i,2) 142 R0(3)=R0(3)+zvv(i)*r(i,3) 143 ENDDO 144 R0=R0/SUM(zvv(1:nat)) 145 ! 146 ! shift of the ions (no PBC) 147 ! 148 DO i=1,nat 149 r(i,1)=(r(i,1)-R0(1)) 150 r(i,2)=(r(i,2)-R0(2)) 151 r(i,3)=(r(i,3)-R0(3)) 152 ENDDO 153 ! 154 ! shift of the electon density 155 ! 156 DO i=1,dfftp%nr1x 157 rgx(i)=(rgx(i)-R0(1))-alat*anint( (rgx(i)-R0(1))/alat ) 158 ENDDO 159 DO i=1,dfftp%nr2x 160 rgy(i)=(rgy(i)-R0(2))-alat*anint( (rgy(i)-R0(2))/alat ) 161 ENDDO 162 DO i=1,dfftp%nr3x 163 rgz(i)=(rgz(i)-R0(3))-alat*anint( (rgz(i)-R0(3))/alat ) 164 ENDDO 165 166 ! 167 ! ions 168 ! 169 charge_ion = SUM(zvv(1:nat)) 170 dipole_ion = 0.D0 171 quadrupole_ion = 0.D0 172 173 DO i = 1, nat 174 DO j = 1, 3 175 dipole_ion(j) = dipole_ion(j) + zvv(i)*r(i,j) 176 quadrupole_ion = quadrupole_ion + zvv(i)*(r(i,j))**2 177 ENDDO 178 ENDDO 179 180 ! 181 ! electrons 182 ! 183 charge_el = 0.0D0 184 dipole_el = 0.0D0 185 quadrupole_el = 0.0D0 186 187 DO i = 1, dfftp%nr1x 188 DO j = 1, dfftp%nr2x 189 DO k = 1, dfftp%nr3x 190 191 charge_el = charge_el + rhof(i,j,k) 192 dipole_el(1) = dipole_el(1) + rgx(i)*rhof(i,j,k) 193 dipole_el(2) = dipole_el(2) + rgy(j)*rhof(i,j,k) 194 dipole_el(3) = dipole_el(3) + rgz(k)*rhof(i,j,k) 195 196 quadrupole_el = quadrupole_el + rhof(i,j,k) * & 197 & ( (rgx(i))**2 + (rgy(j))**2 + (rgz(k))**2 ) 198 199 ENDDO 200 ENDDO 201 ENDDO 202 charge_el=charge_el*alat**3/DBLE(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x) 203 dipole_el=dipole_el*alat**3/DBLE(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x) 204 quadrupole_el=quadrupole_el*alat**3/DBLE(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x) 205 206 ! ... compute ionic+electronic total charge, dipole and quadrupole moments 207 ! 208 charge = -charge_el + charge_ion 209 dipole = -dipole_el + dipole_ion 210 quadrupole = -quadrupole_el + quadrupole_ion 211 ! 212 ! 213 WRITE( stdout, * )"total charge of the system ",charge 214 WRITE( stdout, '(/5X,"charge density inside the ", & 215 & "Wigner-Seitz cell:",3F14.8," el.")' ) charge_el 216 ! 217 WRITE( stdout, & 218 '(/5X,"reference position (R0):",5X,3F14.8," bohr")' ) R0(:) 219 ! 220 ! ... A positive dipole goes from the - charge to the + charge. 221 ! 222 WRITE( stdout, '(/5X,"Dipole moments (with respect to x0):")' ) 223 WRITE( stdout, '( 5X,"Elect",3F10.4," au, ", 3F10.4," Debye")' ) & 224 (-dipole_el(ip), ip = 1, 3), (-dipole_el(ip)*au_debye, ip = 1, 3 ) 225 WRITE( stdout, '( 5X,"Ionic",3F10.4," au, ", 3F10.4," Debye")' ) & 226 ( dipole_ion(ip),ip = 1, 3), ( dipole_ion(ip)*au_debye,ip = 1, 3 ) 227 WRITE( stdout, '( 5X,"Total",3F10.4," au, ", 3F10.4," Debye")' ) & 228 ( dipole(ip), ip = 1, 3), ( dipole(ip)*au_debye, ip = 1, 3 ) 229 ! 230 ! ... print the electronic, ionic and total quadrupole moments 231 ! 232 WRITE( stdout, '(/5X,"Electrons quadrupole moment",F20.8," a.u.")' ) & 233 -quadrupole_el 234 WRITE( stdout, '( 5X," Ions quadrupole moment",F20.8," a.u.")' ) & 235 quadrupole_ion 236 WRITE( stdout, '( 5X," Total quadrupole moment",F20.8," a.u.")' ) & 237 quadrupole 238 ! 239 ! ... Makov-Payne correction, PRB 51, 43014 (1995) 240 ! ... Note that Eq. 15 has the wrong sign for the quadrupole term 241 ! 242 ! 1 / 2 Ry -> a.u. 243 corr1 = - 2.8373D0 / alat * charge**2 / 2.0D0 244 ! 245 aa = quadrupole 246 bb = dipole(1)**2 + dipole(2)**2 + dipole(3)**2 247 ! 248 corr2 = ( 2.D0 / 3.D0 * pi )*( charge*aa - bb ) / alat**3 249 ! 250 ! ... print the Makov-Payne correction 251 ! 252 WRITE( stdout, '(/,5X,"********* MAKOV-PAYNE CORRECTION *********")' ) 253 ! 254 WRITE( stdout,'(/5X,"Makov-Payne correction ",F14.8," a.u. = ",F6.3, & 255 & " eV (1st order, 1/a0)")' ) -corr1, -corr1*autoev 256 WRITE( stdout,'( 5X," ",F14.8," a.u. = ",F6.3, & 257 & " eV (2nd order, 1/a0^3)")' ) -corr2, -corr2*autoev 258 WRITE( stdout,'( 5X," ",F14.8," a.u. = ",F6.3, & 259 & " eV (total)")' ) -corr1-corr2, (-corr1-corr2)*autoev 260 ! 261 WRITE( stdout,'(/5X,"corrected Total energy = ",F14.8," a.u.")' ) & 262 etot - corr1 - corr2 263! 264#if defined __MPI 265ENDIF ! ionode 266#endif 267! 268IF ( ALLOCATED( rhodist1 ) ) DEALLOCATE( rhodist1 ) 269IF ( ALLOCATED( rhodist2 ) ) DEALLOCATE( rhodist2 ) 270#if defined __MPI 271IF ( ALLOCATED( displs ) ) DEALLOCATE( displs ) 272IF ( ALLOCATED( recvcount ) ) DEALLOCATE( recvcount ) 273#endif 274IF ( ALLOCATED( r ) ) DEALLOCATE( r ) 275IF ( ALLOCATED( rgx ) ) DEALLOCATE( rgx ) 276IF ( ALLOCATED( rgy ) ) DEALLOCATE( rgy ) 277IF ( ALLOCATED( rgz ) ) DEALLOCATE( rgz ) 278IF ( ALLOCATED( zvv ) ) DEALLOCATE( zvv ) 279IF ( ALLOCATED( rhof ) ) DEALLOCATE( rhof ) 280! 281RETURN 282END 283