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