1!
2! Copyright (C) 2002-2005 FPMD-CPV groups
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!
8MODULE efield_module
9
10  USE kinds, ONLY : DP
11
12  IMPLICIT NONE
13  SAVE
14
15  logical      :: tefield  = .FALSE.
16  logical      :: tefield2 = .FALSE.
17  integer      :: epol     = 3 !direction electric field
18  real(kind=DP) :: efield   = 0.d0 !intensity electric field
19  real(kind=DP)  :: efield2 =0.d0
20  real(kind=DP)    evalue!strenght of electric field
21  real(kind=DP)  evalue2
22  integer epol2,ipolp2
23  integer         ipolp  !direction of electric field
24
25  real(kind=DP) :: pberryel = 0.0d0, pberryion = 0.0d0
26  real(kind=DP) :: pberryel2 = 0.0d0, pberryion2 = 0.0d0
27
28!***
29!***  Berry phase
30!***
31      integer, allocatable:: ctable(:,:,:)!correspondence tables for diff. polarization
32      integer, allocatable:: ctabin(:,:,:)!inverse correspondence table
33      complex(DP), allocatable:: qmat(:,:)!inverse of matrix Q, for Barry's phase
34      complex(DP), allocatable:: gqq(:,:,:,:)!factors int beta_Ri^*beta_Rj exp(iGr)dr
35      complex(DP), allocatable:: gqqm(:,:,:,:)! the same with exp(-iGr)
36      complex(DP), allocatable:: gqq0(:,:,:,:)!factors int beta_Ri^*beta_Rj exp(iGr)dr, at Gamma
37      complex(DP), allocatable:: gqqm0(:,:,:,:)! the same with exp(-iGr), at Gamma
38      complex(DP), allocatable:: df(:)
39      integer, allocatable:: ctable2(:,:,:)!correspondence tables for diff. polarization
40      integer, allocatable:: ctabin2(:,:,:)!inverse correspondence table
41      complex(DP), allocatable:: qmat2(:,:)!inverse of matrix Q, for Barry's phase
42      complex(DP), allocatable:: gqq2(:,:,:,:)!factors int beta_Ri^*beta_Rj exp(iGr)dr
43      complex(DP), allocatable:: gqqm2(:,:,:,:)! the same with exp(-iGr)
44      complex(DP), allocatable:: gqq02(:,:,:,:)!factors int beta_Ri^*beta_Rj exp(iGr)dr, at Gamma
45      complex(DP), allocatable:: gqqm02(:,:,:,:)! the same with exp(-iGr), at Gamma
46      complex(DP) detq
47      complex(DP) detq2
48      real(DP) cdzp(3),cdzm(3), cdz0(3)!centers of ionic charges
49
50!for parallelization for direcions 1 and 2
51      integer :: n_g_missing_p(2)!number of g vector with correspondence G-->G+1 is missing
52      integer :: n_g_missing_m(2)!number of g vector with correspondence G-->G-1 is missing
53      integer, allocatable :: whose_is_g(:) !correspondence G(plane waves, global) ---> processor
54      integer, allocatable :: ctable_missing_1(:,:,:)!correspondence G(plane waves local)--> array for mpi_alltoall
55                                              !n_g_missing*nproc
56      integer, allocatable :: ctable_missing_rev_1(:,:,:)!missing_g --> G (plane waves local)
57      integer, allocatable :: ctable_missing_2(:,:,:)!correspondence G(plane waves local)--> array for mpi_alltoall
58                                              !n_g_missing*nproc
59      integer, allocatable :: ctable_missing_rev_2(:,:,:)!missing_g --> G (plane waves local)
60      integer, allocatable :: ctabin_missing_1(:,:,:)!correspondence G(plane waves local)--> array for mpi_alltoall
61                                              !n_g_missing*nproc
62      integer, allocatable :: ctabin_missing_rev_1(:,:,:)!missing_g --> G (plane waves local)
63      integer, allocatable :: ctabin_missing_2(:,:,:)!correspondence G(plane waves local)--> array for mpi_alltoall
64                                              !n_g_missing*nproc
65      integer, allocatable :: ctabin_missing_rev_2(:,:,:)!missing_g --> G (plane waves local)
66
67CONTAINS
68
69
70  SUBROUTINE efield_init( epol_ , efield_ )
71    USE kinds, ONLY: DP
72    REAL(DP), INTENT(IN) :: efield_
73    INTEGER, INTENT(IN)    :: epol_
74    epol   = epol_
75    efield = efield_
76    RETURN
77  END SUBROUTINE efield_init
78
79  SUBROUTINE efield_info( )
80    USE io_global, ONLY: ionode,stdout
81    if(ionode) write (stdout,401) epol, efield
82
83401   format (/4x,'====================================='                       &
84     &        /4x,'|  BERRY PHASE ELECTRIC FIELD 1        '                       &
85     &        /4x,'====================================='                       &
86     &        /4x,'| direction    =',i10,'            '                         &
87     &        /4x,'| intensity    =',f10.5,' a.u.     '                         &
88     &        /4x,'=====================================')
89
90    RETURN
91  END SUBROUTINE efield_info
92
93
94  SUBROUTINE efield_berry_setup( eigr, tau0 )
95    USE io_global, ONLY: ionode,stdout
96    IMPLICIT NONE
97    COMPLEX(DP), INTENT(IN)  :: eigr(:,:)
98    REAL(DP), INTENT(IN)  :: tau0(:,:)
99    if(ionode) write(stdout,'("Initialize Berry phase electric field")')
100    ipolp = epol
101    evalue = efield
102!set up for parallel calculations
103
104#if defined(__MPI)
105    call find_whose_is_g
106    call gtable_missing
107    call gtable_missing_inv
108#endif
109
110    call gtable(ipolp,ctable(1,1,ipolp))
111    call gtablein(ipolp,ctabin(1,1,ipolp))
112    call qqberry2(gqq0,gqqm0,ipolp)!for Vanderbilt pps
113    call qqupdate(eigr,gqqm0,gqq,gqqm,ipolp)
114    !the following line was to keep the center of charge fixed
115    !when performing molecular dynamics in the presence of an electric
116    !field
117    !call cofcharge(tau0,cdz0)
118
119    RETURN
120  END SUBROUTINE efield_berry_setup
121
122
123  SUBROUTINE efield_update( eigr )
124    IMPLICIT NONE
125    COMPLEX(DP), INTENT(IN)  :: eigr(:,:)
126    call qqupdate(eigr,gqqm0,gqq,gqqm,ipolp)
127    RETURN
128  END SUBROUTINE efield_update
129
130
131  SUBROUTINE allocate_efield( ngw, ngw_g, nx, nhx, nax, nsp )
132    IMPLICIT NONE
133    INTEGER, INTENT(IN) :: ngw, ngw_g, nx, nhx, nax, nsp
134      allocate( ctable(ngw,2,3))
135      allocate( ctabin(ngw,2,3))
136      allocate( qmat(nx,nx))
137      allocate( gqq(nhx,nhx,nax,nsp))
138      allocate( gqqm(nhx,nhx,nax,nsp))
139      allocate( df(ngw))
140      allocate( gqq0(nhx,nhx,nax,nsp))
141      allocate( gqqm0(nhx,nhx,nax,nsp))
142      allocate( whose_is_g(ngw_g))
143
144    RETURN
145  END SUBROUTINE allocate_efield
146
147
148  SUBROUTINE deallocate_efield( )
149    IMPLICIT NONE
150    IF( allocated( ctable ) )  deallocate( ctable )
151    IF( allocated( ctabin ) ) deallocate( ctabin )
152    IF( allocated( qmat ) ) deallocate( qmat )
153    IF( allocated( gqq ) ) deallocate( gqq )
154    IF( allocated( gqqm ) ) deallocate( gqqm )
155    IF( allocated( df ) ) deallocate( df )
156    IF( allocated( gqq0 ) ) deallocate( gqq0 )
157    IF( allocated( gqqm0 ) )  deallocate( gqqm0 )
158    IF( allocated( whose_is_g) ) deallocate(whose_is_g)
159    IF( allocated( ctable_missing_1) ) deallocate( ctable_missing_1)
160    IF( allocated( ctable_missing_2) ) deallocate( ctable_missing_2)
161    IF( allocated( ctable_missing_rev_1) ) deallocate( ctable_missing_rev_1)
162    IF( allocated( ctable_missing_rev_1) ) deallocate( ctable_missing_rev_2)
163    IF( allocated( ctabin_missing_1) ) deallocate( ctabin_missing_1)
164    IF( allocated( ctabin_missing_2) ) deallocate( ctabin_missing_2)
165    IF( allocated( ctabin_missing_rev_1) ) deallocate( ctabin_missing_rev_1)
166    IF( allocated( ctabin_missing_rev_1) ) deallocate( ctabin_missing_rev_2)
167    RETURN
168  END SUBROUTINE deallocate_efield
169
170
171  SUBROUTINE berry_energy( enb, enbi, bec, cm, fion )
172    USE ions_positions, ONLY: tau0
173    USE control_flags, ONLY: tfor, tprnfor
174    IMPLICIT NONE
175    real(DP), intent(out) :: enb, enbi
176    real(DP) :: bec(:,:)
177    real(DP) :: fion(:,:)
178    complex(DP) :: cm(:,:)
179    real(dp), external :: enberry
180
181    call qmatrixd(cm,bec,ctable(1,1,ipolp),gqq,qmat,detq,ipolp)
182    enb =  enberry( detq, ipolp )
183    call berryion(tau0,fion,tfor.or.tprnfor,ipolp,evalue,enbi)
184    pberryel=enb
185    pberryion=enbi
186    enb=enb*evalue
187    enbi=enbi*evalue
188  END SUBROUTINE berry_energy
189
190
191  SUBROUTINE dforce_efield (bec,i,cm,c2,c3,rhos)
192    USE uspp, ONLY: betae => vkb, deeq
193    USE gvecw, ONLY: ngw
194    IMPLICIT NONE
195    complex(DP), intent(out) :: c2(:), c3(:)
196    complex(DP), intent(in) :: cm(:,:)
197    REAL(DP) :: rhos(:,:)
198    real(DP) :: bec(:,:)
199    integer :: i
200    integer :: ig
201    call dforceb (cm, i, betae, ipolp, bec ,ctabin(1,1,ipolp), gqq, gqqm, qmat, deeq, df)
202    do ig=1,ngw
203      c2(ig)=c2(ig)+evalue*df(ig)
204    enddo
205    call dforceb (cm, i+1, betae, ipolp, bec ,ctabin(1,1,ipolp), gqq, gqqm, qmat, deeq, df)
206    do ig=1,ngw
207      c3(ig)=c3(ig)+evalue*df(ig)
208    enddo
209  END SUBROUTINE dforce_efield
210
211 SUBROUTINE efield_init2( epol_ , efield_ )
212    USE kinds, ONLY: DP
213    REAL(DP), INTENT(IN) :: efield_
214    INTEGER, INTENT(IN)    :: epol_
215    epol2   = epol_
216    efield2 = efield_
217    RETURN
218  END SUBROUTINE efield_init2
219
220  SUBROUTINE efield_info2( )
221    USE io_global, ONLY: ionode,stdout
222    if(ionode) write (stdout,402) epol2, efield2
223
224402   format (/4x,'====================================='                       &
225     &        /4x,'|  BERRY PHASE ELECTRIC FIELD 2        '                       &
226     &        /4x,'====================================='                       &
227     &        /4x,'| direction    =',i10,'            '                         &
228     &        /4x,'| intensity    =',f10.5,' a.u.     '                         &
229     &        /4x,'=====================================')
230
231    RETURN
232  END SUBROUTINE efield_info2
233
234
235  SUBROUTINE efield_berry_setup2( eigr, tau0 )
236    USE io_global, ONLY: ionode,stdout
237    IMPLICIT NONE
238    COMPLEX(DP), INTENT(IN)  :: eigr(:,:)
239    REAL(DP), INTENT(IN)  :: tau0(:,:)
240    if(ionode) write(stdout,'("Initialize Berry phase electric field")')
241    ipolp2 = epol2
242    evalue2 = efield2
243    call gtable(ipolp2,ctable2(1,1,ipolp2))
244    call gtablein(ipolp2,ctabin2(1,1,ipolp2))
245    call qqberry2(gqq02,gqqm02,ipolp2)!for Vanderbilt pps
246    call qqupdate(eigr,gqqm02,gqq2,gqqm2,ipolp2)
247    !the following line was to keep the center of charge fixed
248    !when performing molecular dynamics in the presence of an electric
249    !field
250    !call cofcharge(tau0,cdz0)
251    RETURN
252  END SUBROUTINE efield_berry_setup2
253
254
255  SUBROUTINE efield_update2( eigr )
256    IMPLICIT NONE
257    COMPLEX(DP), INTENT(IN)  :: eigr(:,:)
258    call qqupdate(eigr,gqqm02,gqq2,gqqm2,ipolp2)
259    RETURN
260  END SUBROUTINE efield_update2
261
262
263  SUBROUTINE allocate_efield2( ngw, nx, nhx, nax, nsp )
264    IMPLICIT NONE
265    INTEGER, INTENT(IN) :: ngw, nx, nhx, nax, nsp
266      allocate( ctable2(ngw,2,3))
267      allocate( ctabin2(ngw,2,3))
268      allocate( qmat2(nx,nx))
269      allocate( gqq2(nhx,nhx,nax,nsp))
270      allocate( gqqm2(nhx,nhx,nax,nsp))
271      allocate( gqq02(nhx,nhx,nax,nsp))
272      allocate( gqqm02(nhx,nhx,nax,nsp))
273    RETURN
274  END SUBROUTINE allocate_efield2
275
276
277  SUBROUTINE deallocate_efield2( )
278    IMPLICIT NONE
279    IF( allocated( ctable2 ) )  deallocate( ctable2 )
280    IF( allocated( ctabin2 ) ) deallocate( ctabin2 )
281    IF( allocated( qmat2 ) ) deallocate( qmat2 )
282    IF( allocated( gqq2 ) ) deallocate( gqq2 )
283    IF( allocated( gqqm2 ) ) deallocate( gqqm2 )
284    IF( allocated( gqq02 ) ) deallocate( gqq02 )
285    IF( allocated( gqqm02 ) )  deallocate( gqqm02 )
286    RETURN
287  END SUBROUTINE deallocate_efield2
288
289
290  SUBROUTINE berry_energy2( enb, enbi, bec, cm, fion )
291    USE ions_positions, ONLY: tau0
292    USE control_flags, ONLY: tfor, tprnfor
293    IMPLICIT NONE
294    real(DP), intent(out) :: enb, enbi
295    real(DP) :: bec(:,:)
296    real(DP) :: fion(:,:)
297    complex(DP) :: cm(:,:)
298    real(dp), external :: enberry
299
300    call qmatrixd(cm,bec,ctable2(1,1,ipolp2),gqq2,qmat2,detq2,ipolp2)
301    enb =  enberry( detq2, ipolp2 )
302    call berryion(tau0,fion,tfor.or.tprnfor,ipolp2,evalue2,enbi)
303    pberryel2=enb
304    pberryion2=enbi
305    enb=enb*evalue2
306    enbi=enbi*evalue2
307  END SUBROUTINE berry_energy2
308
309
310  SUBROUTINE dforce_efield2 (bec,i,cm,c2,c3,rhos)
311    USE uspp, ONLY: betae => vkb, deeq
312    USE gvecw, ONLY: ngw
313    IMPLICIT NONE
314    complex(DP), intent(out) :: c2(:), c3(:)
315    complex(DP), intent(in) :: cm(:,:)
316    REAL(DP) :: rhos(:,:)
317    real(DP) :: bec(:,:)
318    integer :: i
319    integer :: ig
320    call dforceb (cm, i, betae, ipolp2, bec ,ctabin2(1,1,ipolp2), gqq2, gqqm2, qmat2, deeq, df)
321    do ig=1,ngw
322      c2(ig)=c2(ig)+evalue2*df(ig)
323    enddo
324    call dforceb (cm, i+1, betae, ipolp2, bec ,ctabin2(1,1,ipolp2), gqq2, gqqm2, qmat2, deeq, df)
325    do ig=1,ngw
326      c3(ig)=c3(ig)+evalue2*df(ig)
327    enddo
328  END SUBROUTINE dforce_efield2
329
330END MODULE efield_module
331