1!
2! Copyright (C) 2012 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!---------------------------------------------------------------------------------
9MODULE bp
10!---------------------------------------------------------------------------------
11  !! Contains the variables needed for the Berry phase polarization calculation
12  !
13  USE kinds,   ONLY: DP
14  USE becmod,  ONLY: bec_type
15  !
16  SAVE
17  PRIVATE
18  PUBLIC :: lberry, lelfield, lorbm, gdir, nppstr, nberrycyc, evcel, evcelp, evcelm, &
19            fact_hepsi, bec_evcel, mapgp_global, mapgm_global, nppstr_3d, ion_pol,   &
20            el_pol, fc_pol, l_el_pol_old, el_pol_old, el_pol_acc, nx_el, l3dstring,  &
21            efield, efield_cart, efield_cry, transform_el, mapg_owner, phase_control
22  PUBLIC :: allocate_bp_efield, deallocate_bp_efield, bp_global_map
23  PUBLIC :: pdl_tot
24  !
25  !
26  LOGICAL :: lberry = .FALSE.
27  !! if .TRUE. calculate polarization using Berry phase
28  LOGICAL :: lelfield = .FALSE.
29  !! if .TRUE. finite electric field using Berry phase
30  LOGICAL :: lorbm = .FALSE.
31  !! if .TRUE. calculate orbital magnetization (Kubo terms)
32  INTEGER :: gdir
33  !! G-vector for polarization calculation
34  INTEGER :: nppstr
35  !! number of k-points (parallel vector)
36  INTEGER :: nberrycyc
37  !! number of cycles for convergence in electric field
38  !! without changing the selfconsistent charge
39  REAL(DP) :: efield
40  !! electric field intensity in a.u.
41  COMPLEX(DP), ALLOCATABLE , TARGET :: evcel(:,:)
42  !! wavefunctions for calculating the electric field operator
43  COMPLEX(DP), ALLOCATABLE , TARGET :: evcelm(:,:,:)
44  !! wavefunctions for  storing projectors for  electric field operator
45  COMPLEX(DP), ALLOCATABLE , TARGET :: evcelp(:,:,:)
46  !! wavefunctions for  storing projectors for  electric field operator
47  COMPLEX(DP), ALLOCATABLE, TARGET :: fact_hepsi(:,:)
48  !! factors for hermitean electric field operators
49  !COMPLEX(DP), ALLOCATABLE, TARGET :: bec_evcel(:,:)
50  !! for storing bec's factors with evcel
51  TYPE(bec_type) :: bec_evcel
52  INTEGER, ALLOCATABLE, TARGET :: mapgp_global(:,:)
53  !! map for G'= G+1 correspondence
54  INTEGER, ALLOCATABLE, TARGET :: mapgm_global(:,:)
55  !! map for G'= G-1 correspondence
56  REAL(DP) :: ion_pol(3)
57  !! the ionic polarization
58  REAL(DP) :: el_pol(3)
59  !! the electronic polarization
60  REAL(DP) :: fc_pol(3)
61  !! the prefactor for the electronic polarization
62  LOGICAL  :: l_el_pol_old
63  !! if true there is already stored a n older value for the polarization
64  !! neeeded for having correct polarization during MD
65  REAL(DP) :: el_pol_old(3)
66  !! the old  electronic polarization
67  REAL(DP) :: el_pol_acc(3)
68  !! accumulator for the electronic polarization
69  !
70  INTEGER :: nppstr_3d(3)
71  !! number of element of strings along the reciprocal directions
72  INTEGER, ALLOCATABLE :: nx_el(:,:)
73  !! index for string to k-point map, (nks*nspin,dir=3)
74  LOGICAL :: l3dstring
75  !! if true strings are on the 3 three directions
76  REAL(DP) :: efield_cart(3)
77  !! electric field vector in cartesian units
78  REAL(DP) :: efield_cry(3)
79  !! electric field vector in crystal units
80  REAL(DP) :: transform_el(3,3)
81  !! transformation matrix from cartesian coordinates to normed reciprocal
82  !! space
83  INTEGER, ALLOCATABLE :: mapg_owner(:,:)
84  REAL(DP) :: pdl_tot
85  !! the total phase calculated from bp_c_phase
86  INTEGER  :: phase_control
87  !! 0 no control, 1 write, 2 read
88  !
89 CONTAINS
90  !
91  !-------------------------------------------------------------------------
92  SUBROUTINE allocate_bp_efield()
93    !-----------------------------------------------------------------------
94    !! Allocate memory for the Berry's phase electric field.
95    ! NOTICE: should be allocated ONLY in parallel case, for gdir=1 or 2
96    !
97    USE gvect,  ONLY : ngm_g
98    !
99    IMPLICIT NONE
100    !
101    IF ( lberry .OR. lelfield .OR. lorbm ) THEN
102       ALLOCATE(mapgp_global(ngm_g,3))
103       ALLOCATE(mapgm_global(ngm_g,3))
104       ALLOCATE(mapg_owner(2,ngm_g))
105    ENDIF
106    !
107    l_el_pol_old = .FALSE.
108    el_pol_acc = 0.d0
109    !
110    !
111    RETURN
112    !
113 END SUBROUTINE allocate_bp_efield
114 !
115 !
116 !--------------------------------------------------------------------------
117 SUBROUTINE deallocate_bp_efield
118   !------------------------------------------------------------------------
119   !! Deallocate memory used in Berry's phase electric field calculation
120   !
121   IMPLICIT NONE
122   !
123   IF ( lberry .OR. lelfield .OR. lorbm ) THEN
124      IF ( ALLOCATED(mapgp_global) ) DEALLOCATE( mapgp_global )
125      IF ( ALLOCATED(mapgm_global) ) DEALLOCATE( mapgm_global )
126      IF ( ALLOCATED(nx_el)        ) DEALLOCATE( nx_el        )
127      IF ( ALLOCATED(mapg_owner)   ) DEALLOCATE( mapg_owner   )
128   ENDIF
129   !
130   !
131   RETURN
132   !
133 END SUBROUTINE deallocate_bp_efield
134 !
135 !
136 !-------------------------------------------------------------------------
137 SUBROUTINE bp_global_map
138   !-----------------------------------------------------------------------
139   !! This subroutine sets up the global correspondence map G+1 and G-1.
140   !
141   USE mp,                   ONLY : mp_sum
142   USE mp_images,            ONLY : me_image, intra_image_comm
143   USE gvect,                ONLY : ngm_g, g, ngm, ig_l2g
144   USE fft_base,             ONLY : dfftp
145   USE cell_base,            ONLY : at
146   !
147   IMPLICIT NONE
148   !
149   INTEGER :: ig, mk1,mk2,mk3, idir, imk(3)
150   INTEGER, ALLOCATABLE :: ln_g(:,:,:)
151   INTEGER, ALLOCATABLE :: g_ln(:,:)
152   !
153   IF ( .NOT.lberry .AND. .NOT. lelfield .AND. .NOT. lorbm ) RETURN
154   ! set up correspondence ln_g ix,iy,iz ---> global g index in
155   ! (for now...) coarse grid
156   ! and inverse realtion global g (coarse) to ix,iy,iz
157   !
158   ALLOCATE( ln_g(-dfftp%nr1:dfftp%nr1,-dfftp%nr2:dfftp%nr2,-dfftp%nr3:dfftp%nr3) )
159   ALLOCATE( g_ln(3,ngm_g) )
160   !
161   ln_g(:,:,:)=0!it means also not found
162   DO ig=1,ngm
163      mk1 = NINT( g(1,ig)*at(1,1) + g(2,ig)*at(2,1) + g(3,ig)*at(3,1) )
164      mk2 = NINT( g(1,ig)*at(1,2) + g(2,ig)*at(2,2) + g(3,ig)*at(3,2) )
165      mk3 = NINT( g(1,ig)*at(1,3) + g(2,ig)*at(2,3) + g(3,ig)*at(3,3) )
166      ln_g(mk1,mk2,mk3)=ig_l2g(ig)
167   ENDDO
168   !
169   CALL mp_sum( ln_g(:,:,:), intra_image_comm )
170   !
171   g_ln(:,:) = 0 !it means also not found
172   !
173   DO ig = 1, ngm
174      mk1 = NINT( g(1,ig)*at(1,1) + g(2,ig)*at(2,1) + g(3,ig)*at(3,1) )
175      mk2 = NINT( g(1,ig)*at(1,2) + g(2,ig)*at(2,2) + g(3,ig)*at(3,2) )
176      mk3 = NINT( g(1,ig)*at(1,3) + g(2,ig)*at(2,3) + g(3,ig)*at(3,3) )
177      g_ln(1,ig_l2g(ig)) = mk1
178      g_ln(2,ig_l2g(ig)) = mk2
179      g_ln(3,ig_l2g(ig)) = mk3
180   ENDDO
181   !
182   CALL mp_sum( g_ln(:,:), intra_image_comm )
183   !
184   !loop on direction
185    DO idir = 1, 3
186       ! for every g on global array find G+1 and G-1 and put on
187       DO ig = 1, ngm_g
188          imk(:) = g_ln(:,ig)
189          imk(idir) = imk(idir)+1
190          ! table array
191          mapgp_global(ig,idir) = ln_g(imk(1),imk(2),imk(3))
192          imk(idir) = imk(idir)-2
193          mapgm_global(ig,idir) = ln_g(imk(1),imk(2),imk(3))
194       ENDDO
195       !
196    ENDDO
197    !
198    mapg_owner = 0
199    DO ig = 1, ngm
200       mapg_owner(1,ig_l2g(ig)) = me_image + 1
201       mapg_owner(2,ig_l2g(ig)) = ig
202    ENDDO
203    !
204    CALL mp_sum( mapg_owner, intra_image_comm )
205    !
206    DEALLOCATE( ln_g, g_ln )
207    !
208    !
209    RETURN
210    !
211  END SUBROUTINE bp_global_map
212  !
213END MODULE bp
214