Line data Source code
1 : !!****m* ABINIT/m_berryphase_new !! NAME
2 : !! m_berryphase_new
3 : !!
4 : !! FUNCTION
5 : !!
6 : !! COPYRIGHT
7 : !! Copyright (C) 2003-2026 ABINIT group (MVeithen)
8 : !! This file is distributed under the terms of the
9 : !! GNU General Public License, see ~abinit/COPYING
10 : !! or http://www.gnu.org/copyleft/gpl.txt .
11 : !!
12 : !! SOURCE
13 :
14 : #if defined HAVE_CONFIG_H
15 : #include "config.h"
16 : #endif
17 :
18 : #include "abi_common.h"
19 :
20 : module m_berryphase_new
21 :
22 : use defs_basis
23 : use defs_wvltypes
24 : use m_efield
25 : use m_errors
26 : use m_abicore
27 : use m_xmpi
28 : use m_hdr
29 : use m_dtset
30 : use m_dtfil
31 :
32 : use defs_datatypes, only : pseudopotential_type
33 : use defs_abitypes, only : MPI_type
34 : use m_berrytk, only : smatrix, polcart
35 : use m_cgprj, only : ctocprj
36 : use m_fftcore, only : kpgsph
37 : use m_geometry, only : xred2xcart, metric
38 : use m_io_tools, only : open_file
39 : use m_iowf, only : outwf, outresid
40 : use m_kg, only : getph
41 : use m_kpts, only : listkk, smpbz
42 : use m_mpinfo, only : proc_distrb_cycle
43 : use m_numeric_tools,only : rhophi
44 : use m_pawang, only : pawang_type
45 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_allgather, &
46 : pawcprj_put, pawcprj_copy, pawcprj_mpi_recv, &
47 : pawcprj_mpi_send, pawcprj_free, pawcprj_getdim, pawcprj_symkn
48 : use m_paw_dfpt, only : dsdr_k_paw
49 : use m_paw_efield, only : pawpolev
50 : use m_pawrad, only : pawrad_type
51 : use m_pawrhoij, only : pawrhoij_type
52 : use m_paw_sphharm, only : setsym_ylm
53 : use m_pawtab, only : pawtab_type
54 : use m_paw_overlap, only : expibi,qijb_kk,smatrix_k_paw
55 : use m_symtk, only : symatm
56 : use m_time, only : timab
57 :
58 : implicit none
59 :
60 : private
61 : !!***
62 :
63 : public :: berryphase_new
64 : public :: prtefield
65 : public :: init_e_field_vars
66 : public :: initberry
67 : public :: update_e_field_vars
68 : !!***
69 :
70 : contains
71 : !!***
72 :
73 : !!****f* ABINIT/berryphase_new
74 : !! NAME
75 : !! berryphase_new
76 : !!
77 : !! FUNCTION
78 : !! This routine computes the Berry Phase polarization
79 : !! and the finite difference expression of the ddk.
80 : !! See for example Na Sai et al., PRB 66, 104108 (2002) [[cite:Sai2002]]
81 : !!
82 : !! INPUTS
83 : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
84 : !! cg(2,mcg)=planewave coefficients of wavefunctions
85 : !! cprj(natom,mcprj*usecrpj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
86 : !! dtfil <type(datafiles_type)>=variables related to files
87 : !! dtset <type(dataset_type)>=all input variables in this dataset
88 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
89 : !! gprimd(3,3)=reciprocal space dimensional primitive translations
90 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
91 : !! indlmn(6,lmnmax,ntypat)
92 : !! array giving l,m,n,lm,ln,spin for i=ln (if useylm=0)
93 : !! or i=lmn (if useylm=1)
94 : !! kg(3,mpw*mkmem)=reduced planewave coordinates
95 : !! lmnmax If useylm=0, max number of (l,m,n) comp. over all type of psps (lnproj)
96 : !! If useylm=1, max number of (l,n) comp. over all type of psps (lmnproj)
97 : !! mband=maximum number of bands
98 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
99 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
100 : !! mkmem=number of k points treated by this node
101 : !! mpi_enreg=information about MPI parallelization
102 : !! mpw=maximum dimensioned size of npw
103 : !! my_natom=number of atoms treated by current processor
104 : !! natom=number of atoms in cell
105 : !! nkpt=number of k points
106 : !! npwarr(nkpt)=number of planewaves in basis at this k point
107 : !! nsppol=1 for unpolarized, 2 for spin-polarized
108 : !! ntypat=number of types of atoms in unit cell
109 : !! nkpt=number of k-points
110 : !! calc_pol_ddk = 1: compute Berryphase polarization
111 : !! 2: compute finite difference expression of the ddk
112 : !! 3: compute polarization & ddk
113 : !! pawrhoij(natom*usepaw) <type(pawrhoij_type)> atomic occupancies
114 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
115 : !! pwind(pwind_alloc,2,3) = array used to compute
116 : !! the overlap matrix smat between k-points (see initberry.f)
117 : !! pwind_alloc = first dimension of pwind
118 : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
119 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
120 : !! typat(natom)=type integer for each atom in cell
121 : !! ucvol=unit cell volume in bohr**3.
122 : !! unit_out= unit for output of the results (usually the .out file of ABINIT)
123 : !! The option unit_out = 0 is allowed. In this case, no information is written
124 : !! to the output file but only to the log file.
125 : !! usecprj=1 if cprj datastructure has been allocated
126 : !! usepaw= 1: use paw framework. 0:do not use paw.
127 : !! xred(3,natom)=reduced atomic coordinates
128 : !! zion(ntypat)=valence charge of each type of atom
129 : !!
130 : !! OUTPUT
131 : !! ptot(3) = total polarization including correction for jumps
132 : !! red_ptot(3) = total polarization including correction for jumps reduced units
133 : !! pel(3) = reduced coordinates of the electronic polarization (a. u.)
134 : !! pelev(3)= expectation value polarization term (PAW only) in cartesian coordinates (already contained in pel)
135 : !! pion(3)= reduced coordinates of the ionic polarization (a. u.)
136 : !!
137 : !! SIDE EFFECTS
138 : !! Input/Output
139 : !! dtefield <type(efield_type)> = variables related to Berry phase
140 : !! and electric field calculations (see initberry.f).
141 : !! In case berryopt = 4, the overlap matrices computed
142 : !! in this routine are stored in dtefield%smat in order
143 : !! to be used in the electric field calculation.
144 : !!
145 : !! TODO
146 : !! - Use the analytical relation between the overlap matrices
147 : !! S(k,k+dk) and S(k+dk,k) to avoid to recompute them when ifor = 2.
148 : !!
149 : !! NOTES
150 : !! - pel and pion do not take into account the factor 1/ucvol
151 : !! - In case of a ddk calculation, the eigenvalues are not computed.
152 : !! - The ddk computed by this routine should not be used to
153 : !! compute the electronic dielectric tensor.
154 : !!
155 : !! SOURCE
156 :
157 277 : subroutine berryphase_new(atindx1,cg,cg13,cprj,dtefield,dtfil,dtset,psps,&
158 277 : & gprimd,hdr,indlmn,kg,lmnmax,mband,mcg,mcg13,mcprj,&
159 277 : & mkmem,mpi_enreg,mpw,my_natom,natom,npwarr,nsppol,ntypat,&
160 277 : & nkpt,calc_pol_ddk,pawrhoij,pawtab,pel,pelev,pion,ptot,red_ptot,pwind,& !!REC
161 277 : & pwind_alloc,pwnsfac,&
162 277 : & rprimd,save_cg13,typat,ucvol,unit_out,usecprj,usepaw,xred,zion)
163 :
164 : !Arguments ------------------------------------
165 : integer, intent(in) :: lmnmax,mband,mcg,mcg13,mcprj,mkmem,mpw,my_natom,natom,nkpt
166 : integer, intent(in) :: nsppol,ntypat,calc_pol_ddk
167 : integer, intent(in) :: pwind_alloc,unit_out,usecprj,usepaw
168 : real(dp), intent(in) :: ucvol
169 : logical,intent(in) :: save_cg13
170 : type(MPI_type), intent(in) :: mpi_enreg
171 : type(datafiles_type), intent(in) :: dtfil
172 : type(dataset_type), intent(in) :: dtset
173 : type(pseudopotential_type),intent(in) :: psps
174 : type(efield_type), intent(inout) :: dtefield
175 : type(hdr_type), intent(inout) :: hdr
176 : !arrays
177 : integer, intent(in) :: atindx1(natom),indlmn(6,lmnmax,ntypat),kg(3,mpw*mkmem)
178 : integer, intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3)
179 : integer, intent(in) :: typat(natom)
180 : real(dp), intent(in) :: cg(2,mcg),gprimd(3,3)
181 : real(dp), intent(in) :: pwnsfac(2,pwind_alloc)
182 : real(dp), intent(in) :: rprimd(3,3),zion(ntypat)
183 : real(dp), intent(inout) :: xred(3,natom)
184 : real(dp), intent(out) :: pel(3),pelev(3),pion(3)
185 : real(dp), intent(out) :: ptot(3),red_ptot(3) !!REC
186 : real(dp),intent(out) :: cg13(2,mcg13,3)
187 : type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*usepaw)
188 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
189 : type(pawcprj_type),intent(in) :: cprj(natom,mcprj*usecprj)
190 :
191 : !Local variables -------------------------
192 : integer :: count,count1,dest,fdir,unt
193 : integer :: iatom,iband,icg,icg1,idir,idum,ikpt1i_sp
194 : integer :: ierr,ifor,ikg,ikpt,ikpt1,ikpt_loc!,ikpt2,ikpt2i,npw_k2, itrs
195 : integer :: icp1, icp2,icpgr_offset,iproc
196 : !! integer :: ii ! appears commented out below in a debug section
197 : integer :: inibz,ikpt1i
198 : integer :: isppol,istr,itypat,jband,jkpt,jkstr,jsppol
199 : integer :: det_inv_smat, det_smat, inv_smat
200 : integer :: maxbd,mcg1_k
201 : integer :: minbd,my_nspinor,nband_k,ncpgr,nfor,npw_k1,ntotcp,n2dim,nproc,pertcase
202 : integer :: response,shiftbd,source,spaceComm,tag
203 : integer :: jj,jstr,kk,ineigh_str
204 554 : integer :: istep,jstep,kpt_mark(dtefield%fnkpt),nkstr,nstr,iunmark,berrystep
205 : integer :: jkpt2, jkpt2i, me
206 : real(dp) :: det_mod,dkinv,dphase,dtm_real,dtm_imag,fac,gmod,phase0
207 : real(dp) :: pol,polbtot,polion,politot,poltot,rho
208 : logical :: calc_epaw3_force,calc_epaw3_stress,efield_flag
209 : integer :: polflag, ddkflag
210 : !!REC start
211 : integer :: jump
212 : real(dp),save :: pol0(3)
213 : logical, save :: first=.true.
214 : logical :: lexist
215 : !!REC end
216 : real(dp) :: dphase_new,dphase_init
217 : character(len=fnlen) :: fiwf1o
218 : character(len=500) :: msg
219 : type(wvl_wf_type) :: wfs
220 277 : type(wvl_internal_type) :: wvl
221 277 : integer,allocatable :: dimlmn(:),ikpt1_recv(:), sflag_k(:)!,pwind_k(:)
222 277 : integer,allocatable :: ikpt3(:), ikpt3i(:), sflag_k_mult(:,:),nattyp_dum(:),npw_k3(:)
223 277 : integer,allocatable :: idxkstr_mult(:,:), pwind_k_mult(:,:),itrs_mult(:)
224 : real(dp) :: det_average(2),dk(3),dtm_k(2),gpard(3),pel_cart(3),pion_cart(3)
225 831 : real(dp) :: polb(nsppol),ptot_cart(3),rel_string(2),xcart(3,natom)
226 : real(dp) :: delta_str(2),dist_,dstr(2)
227 277 : real(dp),allocatable :: buffer(:,:),buffer1(:),buffer2(:)
228 277 : real(dp),allocatable :: cg1(:,:),cg1_k(:,:),cgq(:,:)
229 277 : real(dp),allocatable :: det_string(:,:),dsdr(:,:,:,:,:),dsdr_sum(:,:,:),dsds_sum(:,:,:)
230 277 : real(dp),allocatable :: eig_dum(:),epawf3_str(:,:,:),epaws3_str(:,:,:)
231 277 : real(dp),allocatable :: occ_dum(:),polberry(:),resid(:),pwnsfac_k(:,:)
232 277 : real(dp),allocatable :: smat_inv(:,:,:),smat_k(:,:,:),smat_k_paw(:,:,:)
233 277 : real(dp),allocatable :: str_flag(:)
234 : ! real(dp),allocatable :: dist_str(:,:),det_string_test(:,:)
235 277 : real(dp),allocatable :: dtm_mult(:,:,:), coef(:,:), polb_mult(:,:)
236 277 : type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kb(:,:),cprj_buf(:,:),cprj_gat(:,:)
237 277 : type(pawcprj_type),allocatable :: cprj_fkn(:,:),cprj_ikn(:,:)
238 : ! integer :: bband,bbs,bra_start,bra_end,ilmn,ipw,ispinor,jlmn,kband,kbs,ket_start,ket_end,klmn,npw_k
239 : ! integer :: nspinor,spnipw,spnshft
240 : ! real(dp) :: err_ovlp,mag_ovlp,max_err_ovlp, ovlp_r, ovlp_i, paw_r, paw_i
241 : ! real(dp) :: tot_r, tot_i
242 : ! real(dp),allocatable :: bra(:,:),ket(:,:)
243 : ! complex(dp) :: cpb,cpk,cterm
244 : ! ***********************************************************************
245 :
246 : !DEBUG
247 : !write(std_out,*)' berryphase_new : enter'
248 : !do ii=1,pwind_alloc
249 : !write(std_out,*)ii,pwnsfac(:,ii)
250 : !end do
251 : !stop
252 : !ENDDEBUG
253 :
254 : !Init MPI
255 277 : spaceComm=mpi_enreg%comm_cell
256 277 : nproc=xmpi_comm_size(spaceComm)
257 277 : me=mpi_enreg%me_kpt
258 :
259 277 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
260 :
261 277 : polflag = 1
262 277 : ddkflag = 1
263 277 : if (calc_pol_ddk == 1) then
264 272 : ddkflag = 0
265 5 : else if (calc_pol_ddk == 2) then
266 277 : polflag = 0
267 : end if
268 :
269 : !allocate(pwind_k(mpw))
270 831 : ABI_MALLOC(pwnsfac_k,(4,mpw))
271 831 : ABI_MALLOC(sflag_k,(dtefield%mband_occ))
272 : !pwind_k(:) = 0
273 47890 : pwnsfac_k(1,:) = 1.0_dp ! bra real
274 47890 : pwnsfac_k(2,:) = 0.0_dp ! bra imag
275 47890 : pwnsfac_k(3,:) = 1.0_dp ! ket real
276 47890 : pwnsfac_k(4,:) = 0.0_dp ! ket imag
277 :
278 19685 : if (maxval(dtset%istwfk(:)) /= 1) then
279 : write(msg, '(3a)' )&
280 0 : 'This routine does not work yet with istwfk /= 1.',ch10,&
281 0 : 'This should have been tested previously ...'
282 0 : ABI_BUG(msg)
283 : end if
284 :
285 277 : if (usepaw == 1 .and. usecprj /= 1) then
286 0 : ABI_BUG('PAW calculation but cprj datastructure has not been allocated !')
287 : end if
288 :
289 277 : if (save_cg13 .AND. (ddkflag /= 1) ) then
290 0 : ABI_BUG('cg13 output requested but ddkflag not set')
291 : end if
292 :
293 277 : if (save_cg13 .AND. (mcg13 /= mcg) ) then
294 0 : ABI_BUG('cg13 output requested but mcg13 /= mcg')
295 : end if
296 :
297 1108 : if (save_cg13 .AND. (size(cg13) /= 2*mcg13*3 ) ) then
298 0 : ABI_BUG('cg13 output requested but cg13 size incorrect')
299 : end if
300 :
301 : ! useful flags for various efield possibilities
302 : efield_flag = (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7 .or. &
303 277 : dtset%berryopt ==14 .or. dtset%berryopt ==16 .or. dtset%berryopt ==17 )
304 216 : calc_epaw3_force = ( efield_flag .and. dtset%optforces /= 0 .and. usepaw == 1 )
305 216 : calc_epaw3_stress = ( efield_flag .and. dtset%optstress /= 0 .and. usepaw == 1 )
306 :
307 277 : mcg1_k = mpw*mband*dtset%nspinor
308 277 : shiftbd = 1
309 277 : if (ddkflag==1) then
310 15 : ABI_MALLOC(cg1,(2,mcg))
311 15 : ABI_MALLOC(eig_dum,(2*mband*mband*nkpt*nsppol))
312 15 : ABI_MALLOC(occ_dum,(mband*nkpt*nsppol))
313 25605 : eig_dum(:) = zero
314 2693 : occ_dum(:) = dtefield%sdeg
315 : end if
316 :
317 : !initialize variable tied to multiple step computation
318 277 : berrystep=dtset%berrystep
319 831 : ABI_MALLOC(ikpt3,(berrystep))
320 554 : ABI_MALLOC(ikpt3i,(berrystep))
321 1108 : ABI_MALLOC(sflag_k_mult,(dtefield%mband_occ,berrystep))
322 554 : ABI_MALLOC(npw_k3,(berrystep))
323 1108 : ABI_MALLOC(pwind_k_mult,(mpw,berrystep))
324 554 : ABI_MALLOC(itrs_mult,(berrystep))
325 1108 : ABI_MALLOC(coef,(berrystep,berrystep))
326 1108 : ABI_MALLOC(polb_mult,(nsppol,berrystep))
327 : !coefficient for berryphase computation
328 915 : coef(:,:) = 0.0_dp
329 566 : do jstep = 1, berrystep
330 289 : coef(jstep,1) = 1.d0/real(jstep*jstep,dp)
331 566 : if(jstep/=1)coef(jstep,1)=coef(jstep,1)/real(1-jstep*jstep,dp)
332 : end do
333 289 : do istep = 2, berrystep
334 349 : do jstep = 1, berrystep
335 60 : coef(jstep, istep) = real(istep*istep,dp)*coef(jstep,istep-1)
336 72 : if(jstep /= istep)coef(jstep, istep)=coef(jstep,istep)/real(istep*istep-jstep*jstep,dp)
337 : end do
338 : end do
339 : !the berryphase using the strings of steps dk, 2*dk, ..., istep*dk is :
340 : !coef(1,istep)*berryphase(dk) + coef(2,istep)*berryphase(2*dk) + ... + coef(istep,istep)*berryphase(istep*dk)
341 : !DEBUG
342 : !write(std_out,*)'coef, sum coef'
343 : !do istep=1,step
344 : !write(std_out,*)coef(:,istep), sum(coef(1:istep,istep))
345 : !end do
346 : !ENDDEBUG
347 :
348 : !allocate(dtm(2,dtefield%fnkpt*nsppol))
349 1108 : ABI_MALLOC(dtm_mult,(2,dtefield%fnkpt*nsppol,berrystep))
350 831 : ABI_MALLOC(cg1_k,(2,mcg1_k))
351 :
352 277 : if (usepaw == 1) then ! cprj allocation
353 99 : ncpgr = cprj(1,1)%ncpgr
354 99 : if ( calc_epaw3_force ) then
355 312 : ABI_MALLOC(dsdr_sum,(natom,3,dtefield%fnkpt*nsppol))
356 312 : ABI_MALLOC(epawf3_str,(natom,3,3))
357 : end if
358 99 : if ( calc_epaw3_stress ) then
359 312 : ABI_MALLOC(dsds_sum,(natom,6,dtefield%fnkpt*nsppol))
360 312 : ABI_MALLOC(epaws3_str,(natom,3,6))
361 : end if
362 297 : ABI_MALLOC(dimlmn,(natom))
363 99 : call pawcprj_getdim(dimlmn,natom,nattyp_dum,ntypat,typat,pawtab,'R')
364 1584 : ABI_MALLOC(cprj_k,(natom,dtefield%nspinor*mband))
365 1584 : ABI_MALLOC(cprj_kb,(natom,dtefield%nspinor*mband))
366 3096 : ABI_MALLOC(cprj_gat,(natom,nproc*dtefield%nspinor*mband))
367 99 : call pawcprj_alloc(cprj_k,ncpgr,dimlmn)
368 99 : call pawcprj_alloc(cprj_kb,ncpgr,dimlmn)
369 99 : call pawcprj_alloc(cprj_gat,ncpgr,dimlmn)
370 99 : if (dtset%kptopt /= 3) then
371 1520 : ABI_MALLOC(cprj_ikn,(natom,dtefield%nspinor*mband))
372 1520 : ABI_MALLOC(cprj_fkn,(natom,dtefield%nspinor*mband))
373 95 : call pawcprj_alloc(cprj_ikn,ncpgr,dimlmn)
374 95 : call pawcprj_alloc(cprj_fkn,ncpgr,dimlmn)
375 : end if
376 :
377 99 : n2dim = dtefield%nspinor*mband
378 : ntotcp = n2dim*SUM(dimlmn(:))
379 99 : if (nproc>1) then
380 864 : ABI_MALLOC(cprj_buf,(natom,dtefield%nspinor*mband))
381 54 : call pawcprj_alloc(cprj_buf,ncpgr,dimlmn)
382 : end if
383 :
384 99 : if ( efield_flag ) then
385 78 : write(msg,'(2a,i5,2a)')ch10,' nkpt = ',nkpt,ch10,' copy cprj to dtefield%cprj '
386 78 : call wrtout(std_out,msg,'COLL')
387 :
388 156 : do isppol = 1, nsppol
389 :
390 78 : ikpt_loc = 0
391 78 : ikpt1 = 0
392 8426 : do while (ikpt_loc < mkmem)
393 :
394 8270 : if (ikpt_loc < mkmem) ikpt1 = ikpt1 + 1
395 8270 : if ((ikpt1 > nkpt).and.(ikpt_loc < mkmem)) exit
396 8270 : nband_k = dtset%nband(ikpt1)
397 :
398 8270 : if ( (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt1,1,nband_k,isppol,me)) .and. (ikpt_loc <= mkmem) ) cycle
399 :
400 5006 : ikpt_loc = ikpt_loc + 1
401 :
402 15018 : ABI_MALLOC(ikpt1_recv,(nproc))
403 5006 : call xmpi_allgather(ikpt1,ikpt1_recv,spaceComm,ierr)
404 : call pawcprj_get(atindx1,cprj_k,cprj,natom,1,(ikpt_loc-1)*nband_k*my_nspinor,ikpt1,0,isppol,mband,&
405 : & mkmem,natom,nband_k,nband_k,my_nspinor,nsppol,0,&
406 5006 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
407 5006 : call pawcprj_mpi_allgather(cprj_k,cprj_gat,natom,n2dim,1,dimlmn,ncpgr,nproc,spaceComm,ierr,rank_ordered=.true.)
408 16540 : do iproc = 1, nproc
409 11534 : icp2=nband_k*(iproc-1)*my_nspinor
410 : call pawcprj_get(atindx1,cprj_k,cprj_gat,natom,1,icp2,ikpt1,0,isppol,mband,&
411 : & nproc,natom,nband_k,nband_k,my_nspinor,1,0,&
412 11534 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
413 11534 : icp1 = nband_k*(ikpt1_recv(iproc)-1)*my_nspinor
414 : call pawcprj_put(atindx1,cprj_k,dtefield%cprj,natom,1,icp1,ikpt1,0,isppol,&
415 : & mband,dtefield%fnkpt,natom,nband_k,nband_k,dimlmn,my_nspinor,nsppol,0,&
416 16540 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
417 : end do
418 10090 : ABI_FREE(ikpt1_recv)
419 :
420 : end do ! close loop over k-points
421 : end do ! end loop over nsppol
422 : end if ! end check on efield
423 : end if
424 :
425 : !!=======================================
426 : !! code to test orthonormality of cg_k
427 : !!=======================================
428 : !
429 : !ikpt = 3
430 : !npw_k = npwarr(ikpt)
431 : !isppol = 1
432 : !nband_k = dtefield%mband_occ
433 : !ABI_MALLOC(bra,(2,npw_k*my_nspinor))
434 : !ABI_MALLOC(ket,(2,npw_k*my_nspinor))
435 : !max_err_ovlp=0.0
436 : !call pawcprj_get(atindx1,cprj_k,cprj,natom,1,dtefield%cprjindex(ikpt,isppol),ikpt,0,isppol,mband,&
437 : !& mkmem,natom,nband_k,nband_k,my_nspinor,nsppol,0)
438 : !do bband = 1, nband_k
439 : !bra_start = dtefield%cgindex(ikpt,nsppol)+1+(bband-1)*npw_k*my_nspinor
440 : !bra_end = bra_start + npw_k*my_nspinor - 1
441 : !bra(1:2,1:npw_k*my_nspinor) = cg(1:2,bra_start:bra_end)
442 : !do kband = 1, nband_k
443 : !ket_start = dtefield%cgindex(ikpt,nsppol)+1+(kband-1)*npw_k*my_nspinor
444 : !ket_end = ket_start + npw_k*my_nspinor - 1
445 : !ket(1:2,1:npw_k*my_nspinor) = cg(1:2,ket_start:ket_end)
446 : !
447 : !tot_r = 0.0; tot_i = 0.0
448 : !do ispinor = 1, my_nspinor
449 : !ovlp_r = 0.0; ovlp_i = 0.0
450 : !spnshft = (ispinor-1)*npw_k
451 : !do ipw = 1, npw_k
452 : !spnipw = ipw + spnshft
453 : !ovlp_r = ovlp_r + bra(1,spnipw)*ket(1,spnipw)+bra(2,spnipw)*ket(2,spnipw)
454 : !ovlp_i = ovlp_i - bra(2,spnipw)*ket(1,spnipw)+bra(1,spnipw)*ket(2,spnipw)
455 : !end do ! end loop over ipw
456 : !paw_r = 0.0; paw_i = 0.0
457 : !do iatom = 1, natom
458 : !itypat = typat(iatom)
459 : !do ilmn = 1, dtefield%lmn_size(itypat)
460 : !do jlmn = 1, dtefield%lmn_size(itypat)
461 : !klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn)
462 : !bbs = my_nspinor*(bband-1)+ispinor
463 : !kbs = my_nspinor*(kband-1)+ispinor
464 : !cpb=cmplx(cprj_k(iatom,bbs)%cp(1,ilmn),cprj_k(iatom,bbs)%cp(2,ilmn))
465 : !cpk=cmplx(cprj_k(iatom,kbs)%cp(1,jlmn),cprj_k(iatom,kbs)%cp(2,jlmn))
466 : !cterm = conjg(cpb)*pawtab(itypat)%sij(klmn)*cpk
467 : !paw_r = paw_r + real(cterm)
468 : !paw_i = paw_i + aimag(cterm)
469 : !end do ! end loop over jlmn
470 : !end do ! end loop over ilmn
471 : !end do ! end loop over iatom
472 : !tot_r = tot_r + ovlp_r + paw_r
473 : !tot_i = tot_i + ovlp_i + paw_i
474 : !end do ! end loop over ispinor
475 : !
476 : !! write(std_out,'(a,2i4,2es16.8)')' JWZ Debug: berryphase_new bband kband ovlp : ',&
477 : !!& bband,kband,tot_r,tot_i
478 : !mag_ovlp = tot_r*tot_r + tot_i*tot_i
479 : !if(bband==kband) then
480 : !err_ovlp=abs(mag_ovlp-1.0)
481 : !else
482 : !err_ovlp=abs(mag_ovlp)
483 : !end if
484 : !max_err_ovlp=MAX(max_err_ovlp,err_ovlp)
485 : !end do ! end loop over kband
486 : !end do ! end loop over bband
487 : !write(std_out,'(a,i4,es16.8)')' JWZ Debug: berrphase_new ikpt ovlp err : ',&
488 : !& ikpt,max_err_ovlp
489 : !ABI_FREE(bra)
490 : !ABI_FREE(ket)
491 : !
492 : !!=========================================
493 : !! end code to test orthonormality of cg_k
494 : !!=========================================
495 :
496 277 : pel(:) = zero ; pelev(:)=zero ; pion(:) = zero ; ptot(:)=zero ; red_ptot(:)=zero
497 :
498 277 : minbd = 1 ; maxbd = dtefield%mband_occ
499 :
500 2617 : if(calc_epaw3_force) dtefield%epawf3(:,:,:) = zero
501 4957 : if(calc_epaw3_stress) dtefield%epaws3(:,:,:) = zero
502 :
503 1108 : do idir = 1, 3
504 :
505 : ! dtm(:,:) = zero
506 1494735 : dtm_mult(:,:,:) = zero
507 3462591 : if (calc_epaw3_force) dsdr_sum(:,:,:) = zero
508 6578175 : if (calc_epaw3_stress) dsds_sum(:,:,:) = zero
509 :
510 831 : if (dtset%rfdir(idir) /= 1) cycle
511 :
512 101784 : if (abs(dtefield%efield_dot(idir)) < tol12) dtefield%sflag(:,:,:,idir) = 0
513 :
514 : ! calculate vector steps in k space
515 3256 : dk(:) = dtefield%dkvecs(:,idir)
516 3256 : gpard(:) = dk(1)*gprimd(:,1) + dk(2)*gprimd(:,2) + dk(3)*gprimd(:,3)
517 3256 : gmod = sqrt(dot_product(gpard,gpard))
518 :
519 814 : write(msg,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,&
520 814 : & ' Computing the polarization (Berry phase) for reciprocal vector:',ch10,&
521 814 : & dk(:),' (in reduced coordinates)',ch10,&
522 1628 : & gpard(1:3),' (in cartesian coordinates - atomic units)'
523 814 : call wrtout(std_out,msg,'COLL')
524 814 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
525 :
526 : write(msg,'(a,i5,a,a,i5)')&
527 814 : & ' Number of strings: ',dtefield%nstr(idir),ch10,&
528 1628 : & ' Number of k points in string:', dtefield%nkstr(idir)
529 814 : call wrtout(std_out,msg,'COLL')
530 814 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
531 :
532 : ! Check whether the polarization or the ddk must be computed
533 :
534 : ! nfor = 1 : to compute P, I only need the WF at k + dk
535 : ! nfor = 2 : to compute the ddk, I need the WF at k + dk and k - dk
536 : ! dkinv : +-1/2dk
537 :
538 :
539 : ! default for polarization
540 814 : nfor = 1
541 814 : if (ddkflag == 1) then
542 13 : nfor = 2
543 : end if
544 :
545 814 : if (ddkflag == 1) then
546 :
547 3478093 : cg1(:,:) = zero
548 13 : dkinv = one/(two*dk(idir))
549 :
550 13 : write(msg,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,&
551 13 : & ' Computing the ddk (Berry phase) for reciprocal vector:',ch10,&
552 13 : & dk(:),' (in reduced coordinates)',ch10,&
553 26 : & gpard(1:3),' (in cartesian coordinates - atomic units)'
554 13 : call wrtout(std_out,msg,'COLL')
555 13 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
556 : end if
557 :
558 : ! From smatrix routine: det_inv_smat = type of calculation
559 : ! 1 : compute inverse of the overlap matrix
560 : ! 10 : compute determinant of the overlap matrix
561 : ! 11 : compute determinant and inverse of the overlap matrix
562 814 : inv_smat = 0
563 814 : det_smat = 0
564 :
565 : ! for ddk need inverse matrix
566 814 : if (ddkflag == 1) then
567 13 : inv_smat = 1
568 : end if
569 :
570 : ! if polarization is requested need smat determinant as well
571 814 : if (polflag == 1) then
572 804 : det_smat = 1
573 : end if
574 :
575 : ! electric fields with PAW also needs S_inverse for forces and stresses, even just for polarization
576 814 : if (calc_epaw3_force .or. calc_epaw3_stress) then
577 234 : inv_smat = 1
578 : end if
579 :
580 814 : det_inv_smat = 10*det_smat + inv_smat
581 :
582 : !--------------------------------------------------------------------
583 : ! for each dk we require, calculate the smatrix, derivatives etc...
584 : !--------------------------------------------------------------------
585 1641 : do ifor = 1, nfor
586 :
587 827 : if (ifor == 2) then
588 52 : dk(:) = -1_dp*dk(:)
589 : ! only the inverse of the overlap matrix is required on second pass, speeds things up a bit
590 13 : det_inv_smat = 1
591 13 : dkinv = -1_dp*dkinv
592 : end if
593 :
594 :
595 : ! Compute the determinant and/or the inverse of the overlap matrix
596 : ! for each pair of k-points < u_nk | u_nk+dk >
597 :
598 827 : icg = 0 ; icg1 = 0
599 3308 : ABI_MALLOC(smat_k,(2,dtefield%mband_occ,dtefield%mband_occ))
600 2481 : ABI_MALLOC(smat_inv,(2,dtefield%mband_occ,dtefield%mband_occ))
601 3308 : ABI_MALLOC(smat_k_paw,(2,usepaw*dtefield%mband_occ,usepaw*dtefield%mband_occ))
602 827 : if (calc_epaw3_force .or. calc_epaw3_stress) then ! dsdr needed for forces and stresses in electric field with PAW
603 1404 : ABI_MALLOC(dsdr,(2,natom,ncpgr,usepaw*dtefield%mband_occ,usepaw*dtefield%mband_occ))
604 240786 : dsdr = zero
605 : end if
606 :
607 :
608 : ! Loop on the values of ikpt_loc and ikpt1 :
609 : ! ikpt1 is incremented one by one, and number the k points in the FBZ
610 : ! ikpt1i refer to the k point numbering in the IBZ
611 : ! ikpt_loc differs from ikpt1 only in the parallel case, and gives
612 : ! the index of the k point in the FBZ, in the set treated by the present processor
613 : ! NOTE : in order to allow synchronisation, ikpt_loc contain information about
614 : ! ikpt AND ISPPOL !
615 : ! It means that the following loop is equivalent to a double loop :
616 : ! do isppol = 1, nsppol
617 : ! do ikpt1 = 1, dtefield%fmkmem
618 : !
619 279913 : do ikpt_loc = 1, dtefield%fmkmem_max*nsppol
620 :
621 279086 : ikpt1=mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc,1)
622 279086 : isppol=mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc,2)
623 :
624 : ! if this k and spin are for me do it
625 279086 : if (ikpt1 > 0 .and. isppol > 0) then
626 :
627 270122 : ikpt1i = dtefield%indkk_f2ibz(ikpt1,1)
628 270122 : nband_k = dtset%nband(ikpt1i + (isppol-1)*dtset%nkpt)
629 :
630 : ! DEBUG
631 : ! Please keep this debugging feature
632 : ! write(std_out,'(a,5i4)' )' berryphase_new : ikpt_loc,ikpt1,isppol,idir,ifor=',ikpt_loc,ikpt1,isppol,idir,ifor
633 : ! ENDDEBUG
634 :
635 270122 : inibz=0
636 : if (dtset%kptns(1,ikpt1i) == dtefield%fkptns(1,ikpt1) .and. &
637 270122 : & dtset%kptns(2,ikpt1i) == dtefield%fkptns(2,ikpt1) .and. &
638 270122 : & dtset%kptns(3,ikpt1i) == dtefield%fkptns(3,ikpt1)) inibz=1
639 :
640 270122 : ikg = dtefield%fkgindex(ikpt1)
641 : ! ikpt2 = dtefield%ikpt_dk(ikpt1,ifor,idir)
642 : ! ikpt2i = dtefield%indkk_f2ibz(ikpt2,1)
643 :
644 : ! ikpt3(istep) : index of kpt1 + istep*dk in the FBZ
645 : ! ikpt3i(istep) : index of kpt1 + istep*dk in the IBZ
646 270122 : ikpt3(1) = dtefield%ikpt_dk(ikpt1,ifor,idir)
647 270122 : ikpt3i(1) = dtefield%indkk_f2ibz(ikpt3(1),1)
648 270890 : do istep = 1, berrystep-1
649 768 : ikpt3(istep+1) = dtefield%ikpt_dk(ikpt3(istep),ifor,idir)
650 270890 : ikpt3i(istep+1) = dtefield%indkk_f2ibz(ikpt3(istep+1),1)
651 : end do
652 :
653 : ! itrs = 0
654 : ! if (dtefield%indkk_f2ibz(ikpt1,6) == 1 ) itrs = itrs + 1
655 : ! if (dtefield%indkk_f2ibz(ikpt2,6) == 1 ) itrs = itrs + 10
656 :
657 541012 : itrs_mult(:)=0
658 398668 : if (dtefield%indkk_f2ibz(ikpt1,6) == 1 ) itrs_mult(:) = itrs_mult(:) + 1
659 541012 : do istep=1,berrystep
660 541012 : if (dtefield%indkk_f2ibz(ikpt3(istep),6) == 1 ) itrs_mult(istep) = itrs_mult(istep) + 10
661 : end do
662 :
663 270122 : npw_k1 = npwarr(ikpt1i)
664 : ! npw_k2 = npwarr(ikpt2i)
665 :
666 541012 : do istep = 1, berrystep
667 541012 : npw_k3(istep)=npwarr(ikpt3i(istep))
668 : end do
669 :
670 : ! ji: the loop is over the FBZ, but sflag and smat only apply to the IBZ
671 270122 : if ( efield_flag .and. inibz == 1) then !!HONG
672 23124 : ikpt1i_sp=ikpt1i+(isppol-1)*dtset%nkpt
673 5750172 : smat_k(:,:,:) = dtefield%smat(:,:,:,ikpt1i_sp,ifor,idir)
674 : else
675 25136406 : smat_k(:,:,:) = zero
676 : end if
677 :
678 : ! pwind_k(1:npw_k1) = pwind(ikg+1:ikg+npw_k1,ifor,idir)
679 31128776 : pwnsfac_k(1,1:npw_k1) = pwnsfac(1,ikg+1:ikg+npw_k1)
680 31128776 : pwnsfac_k(2,1:npw_k1) = pwnsfac(2,ikg+1:ikg+npw_k1)
681 :
682 : ! the array needed to compute the overlap matrix between k and k+istep*dk (with multiple steps)
683 : ! the 0-case (no corresponding pw in k and k+dk) could be handled better (k+2*dk could have a corresponding pw ?)
684 31128776 : pwind_k_mult(1:npw_k1,1)=pwind(ikg+1:ikg+npw_k1,ifor,idir)
685 270890 : do istep = 1, berrystep-1
686 327410 : do jj=1, npw_k1
687 57288 : if(pwind_k_mult(jj,istep)/=0)then
688 45528 : pwind_k_mult(jj,istep+1) = pwind(dtefield%fkgindex(ikpt3(istep))+pwind_k_mult(jj,istep),ifor,idir)
689 : else
690 10992 : pwind_k_mult(jj,istep+1) = 0
691 : end if
692 : end do
693 : end do
694 :
695 : ! DEBUG
696 : ! write(std_out,*)' berryphase_new : dtset%berryopt,inibz,ikpt1i,isppol,dtset%nkpt,ifor,idir', &
697 : ! & dtset%berryopt,inibz,ikpt1i,isppol,dtset%nkpt,ifor,idir
698 : ! write(std_out,'(a,4i4)' )' berryphase_new : sflag_k(:)=',sflag_k(:)
699 : ! ENDDEBUG
700 :
701 270122 : if ( efield_flag .and. inibz == 1) then !!HONG
702 23124 : ikpt1i_sp=ikpt1i+(isppol-1)*dtset%nkpt
703 158220 : sflag_k(:) = dtefield%sflag(:,ikpt1i_sp,ifor,idir)
704 : else
705 1348806 : sflag_k(:) = 0
706 : end if
707 :
708 270122 : if (usepaw == 1) then
709 169440 : icp1=dtefield%cprjindex(ikpt1i,isppol)
710 : call pawcprj_get(atindx1,cprj_k,cprj,natom,1,icp1,ikpt1i,0,isppol,&
711 : & mband,mkmem,natom,dtefield%mband_occ,dtefield%mband_occ,&
712 : & my_nspinor,nsppol,0,mpicomm=mpi_enreg%comm_kpt,&
713 169440 : & proc_distrb=mpi_enreg%proc_distrb)
714 :
715 169440 : if ( ikpt1i /= ikpt1 ) then
716 168720 : call pawcprj_copy(cprj_k,cprj_ikn)
717 : call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,indlmn,&
718 : & dtefield%indkk_f2ibz(ikpt1,2),dtefield%indkk_f2ibz(ikpt1,6),&
719 : & dtefield%fkptns(:,dtefield%i2fbz(ikpt1i)),&
720 : & dtefield%lmax,dtefield%lmnmax,mband,natom,dtefield%mband_occ,my_nspinor,&
721 168720 : & dtefield%nsym,ntypat,typat,dtefield%zarot)
722 168720 : call pawcprj_copy(cprj_fkn,cprj_k)
723 : end if
724 :
725 : end if ! end if usepaw
726 :
727 : ! DEBUG
728 : ! write(std_out,'(a,4i4)' )' berryphase_new : sflag_k(:)=',sflag_k(:)
729 : ! ENDDEBUG
730 :
731 : ! DEBUG
732 : ! write(std_out,'(a,7i4)')'me, idir,ifor, ikpt_loc, ikpt1, isppol = ',&
733 : ! & me,idir,ifor,ikpt_loc,ikpt1,isppol
734 : ! write(std_out,'(a,10i3)')'pwind_k(1:10) = ',pwind_k(1:10)
735 : ! ENDDEBUG
736 :
737 541012 : do istep=1,berrystep
738 1780988 : sflag_k_mult(:,istep) = sflag_k(:)
739 : end do
740 :
741 : end if ! end check that ikpt1 > 0 and isppol > 0
742 :
743 : ! --------------------------------------------------------------------------------
744 : ! Communication
745 : ! --------------------------------------------------------------------------------
746 :
747 558940 : do istep=1,berrystep
748 :
749 : ! if(ikpt_loc <= nsppol*dtefield%fmkmem) then
750 279854 : if (ikpt1 > 0 .and. isppol > 0) then ! I currently have a true kpt to use
751 :
752 270890 : count = npw_k3(istep)*my_nspinor*nband_k
753 812670 : ABI_MALLOC(cgq,(2,count))
754 470315054 : cgq = zero
755 270890 : source = me
756 270890 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt3i(istep),1,nband_k,isppol,me)) then
757 : ! I need the datas from someone else
758 42376 : source = mpi_enreg%proc_distrb(ikpt3i(istep),1,isppol)
759 : end if
760 : else
761 8964 : source = -1 ! I do not have a kpt to use
762 : end if
763 :
764 806080 : do dest = 0, nproc-1
765 :
766 806080 : if ((dest==me) .and. (ikpt1>0) .and. (isppol>0)) then
767 : ! I am destination and I have something to do
768 : ! if (mpi_enreg%paral_compil_kpt == 1) write(std_out,*) &
769 : ! & 'coucou 2, mpi_enreg%proc_distrb(ikpt3i(istep),1:nband_k,isppol) : ', &
770 : ! & mpi_enreg%proc_distrb(ikpt3i(istep),1:nband_k,isppol)
771 : ! write(std_out,*)'ikpt3i(istep) ', ikpt3i(istep)
772 : ! write(std_out,*)'nband_k ',nband_k
773 : ! write(std_out,*)'isppol ', isppol
774 : ! write(std_out,*)'mpi_enreg%proc_distrb',mpi_enreg%proc_distrb
775 :
776 270890 : if (source == me) then
777 : ! I am destination and source
778 : ! DEBUG
779 : ! write(std_out,*)'copying ... '
780 : ! write(std_out,*)'me: ',me, 'ikpt3i(istep) ', ikpt3i(istep), 'isppol ', isppol
781 : ! ENDDEBUG
782 :
783 : ! pwnsfac
784 228514 : idum = dtefield%fkgindex(ikpt3(istep))
785 27627668 : pwnsfac_k(3,1:npw_k3(istep)) = pwnsfac(1,idum+1:idum+npw_k3(istep))
786 27627668 : pwnsfac_k(4,1:npw_k3(istep)) = pwnsfac(2,idum+1:idum+npw_k3(istep))
787 :
788 : ! cgq (and cprj)
789 228514 : icg1 = dtefield%cgindex(ikpt3i(istep),isppol)
790 :
791 228514 : if (usepaw == 1) then
792 128556 : icp2=dtefield%cprjindex(ikpt3i(istep),isppol)
793 : call pawcprj_get(atindx1,cprj_kb,cprj,natom,1,icp2,ikpt3i(istep),0,isppol,&
794 : & mband,mkmem,natom,dtefield%mband_occ,dtefield%mband_occ,my_nspinor,&
795 128556 : & nsppol,0,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
796 : end if
797 :
798 405257968 : cgq(:,1:count) = cg(:,icg1+1:icg1+count)
799 : ! if (usepaw == 1) then
800 : ! call pawcprj_copy(cprj_buf,cprj_kb)
801 : ! end if
802 :
803 : ! if ((source /= me)) then
804 : else
805 : ! I am the destination but not the source -> receive
806 : ! DEBUG
807 : ! write(std_out,'(a)')'receiving ...'
808 : ! write(std_out,'(a,i4,a,i4,a,i4,a,i4)')'me: ',me, 'source ', source,'ikpt3i(istep) ', ikpt3i(istep), 'isppol ', isppol
809 : ! ENDDEBUG
810 :
811 : ! receive pwnsfac
812 127128 : ABI_MALLOC(buffer,(2,npw_k3(istep)))
813 42376 : tag = ikpt3(istep) + (isppol - 1)*dtefield%fnkpt
814 42376 : call xmpi_recv(buffer,source,tag,spaceComm,ierr)
815 3558396 : pwnsfac_k(3,1:npw_k3(istep)) = buffer(1,1:npw_k3(istep))
816 3558396 : pwnsfac_k(4,1:npw_k3(istep)) = buffer(2,1:npw_k3(istep))
817 42376 : ABI_FREE(buffer)
818 :
819 : ! receive cgq (and cprj)
820 42376 : tag = ikpt3i(istep) + (isppol - 1)*nkpt
821 42376 : call xmpi_recv(cgq,source,tag,spaceComm,ierr)
822 :
823 84752 : if (usepaw == 1) then
824 40884 : call pawcprj_mpi_recv(natom,n2dim,dimlmn,ncpgr,cprj_kb,source,spaceComm,ierr)
825 : end if
826 :
827 : end if
828 :
829 255336 : else if (dest /= me) then
830 :
831 : ! jkpt is the kpt which is being treated by dest
832 : ! jsppol is his isppol
833 246372 : jkpt = mpi_enreg%kpt_loc2fbz_sp(dest, ikpt_loc,1)
834 246372 : jsppol = mpi_enreg%kpt_loc2fbz_sp(dest, ikpt_loc,2)
835 :
836 246372 : if (jkpt > 0 .and. jsppol > 0) then ! dest is treating a true kpt
837 :
838 227328 : jkpt2 = dtefield%ikpt_dk(jkpt,ifor,idir)
839 227328 : jkpt2i = dtefield%indkk_f2ibz(jkpt2,1)
840 :
841 : ! check if I am his source
842 227328 : if((mpi_enreg%proc_distrb(jkpt2i,1,jsppol) == me)) then
843 : ! I know something about jkpt3i and I must send it
844 : ! DEBUG
845 : ! write(std_out,'(a)')'sending ...'
846 : ! write(std_out,'(a,i4,a,i4,a,i4,a,i4)')'dest: ',dest,' me: ',me,&
847 : ! & ' jkpt2i ',jkpt2i,' jsppol: ',jsppol
848 : ! ENDDEBUG
849 :
850 : ! pwnsfac
851 42376 : tag = jkpt2 + (jsppol - 1)*dtefield%fnkpt
852 42376 : count1 = npwarr(jkpt2i)
853 127128 : ABI_MALLOC(buffer,(2,count1))
854 42376 : idum = dtefield%fkgindex(jkpt2)
855 3558396 : buffer(1,1:count1) = pwnsfac(1,idum+1:idum+count1)
856 3558396 : buffer(2,1:count1) = pwnsfac(2,idum+1:idum+count1)
857 42376 : call xmpi_send(buffer,dest,tag,spaceComm,ierr)
858 42376 : ABI_FREE(buffer)
859 :
860 : ! cgq (and cprj)
861 42376 : icg1 = dtefield%cgindex(jkpt2i,jsppol)
862 :
863 42376 : if (usepaw == 1) then
864 40884 : icp2=dtefield%cprjindex(jkpt2i,jsppol)
865 : call pawcprj_get(atindx1,cprj_buf,cprj,natom,1,icp2,jkpt2i,0,jsppol,&
866 : & mband,mkmem,natom,dtefield%mband_occ,dtefield%mband_occ,&
867 : & my_nspinor,nsppol,0,mpicomm=mpi_enreg%comm_kpt,&
868 40884 : & proc_distrb=mpi_enreg%proc_distrb)
869 : end if
870 :
871 42376 : tag = jkpt2i + (jsppol - 1)*nkpt
872 42376 : count1 = npwarr(jkpt2i)*my_nspinor*nband_k
873 127128 : ABI_MALLOC(buffer,(2,count1))
874 65057086 : buffer(:,1:count1) = cg(:,icg1+1:icg1+count1)
875 42376 : call xmpi_send(buffer,dest,tag,spaceComm,ierr)
876 42376 : ABI_FREE(buffer)
877 :
878 127128 : if (usepaw == 1 ) then
879 40884 : call pawcprj_mpi_send(natom,n2dim,dimlmn,ncpgr,cprj_buf,dest,spaceComm,ierr)
880 : end if
881 :
882 : end if ! end check that I am his source
883 : end if ! end check that jkpt > 0 and jsppol > 0
884 : end if ! end if statements on dest == me or dest /= me
885 : end do ! end loop over dest = 0, nproc - 1
886 :
887 558940 : if (ikpt1 > 0 .and. isppol > 0) then ! if I am treating a kpt, compute the smatrix
888 :
889 270890 : if (usepaw == 1) then
890 169440 : if (ikpt3(istep) /= ikpt3i(istep)) then ! cprj_kb refers to ikpt3i(istep), must compute ikpt3(istep) value
891 168720 : call pawcprj_copy(cprj_kb,cprj_ikn)
892 :
893 : call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,indlmn,&
894 : & dtefield%indkk_f2ibz(ikpt3(istep),2),dtefield%indkk_f2ibz(ikpt3(istep),6),&
895 : & dtefield%fkptns(:,dtefield%i2fbz(ikpt3i(istep))),&
896 : & dtefield%lmax,dtefield%lmnmax,mband,natom,&
897 : & dtefield%mband_occ,my_nspinor,dtefield%nsym,ntypat,typat,&
898 168720 : & dtefield%zarot)
899 168720 : call pawcprj_copy(cprj_fkn,cprj_kb)
900 : end if
901 169440 : call smatrix_k_paw(cprj_k,cprj_kb,dtefield,idir,ifor,mband,natom,smat_k_paw,typat)
902 : ! write(std_out,'(a,5i4)')' JWZ berryphase_new : ikpt_loc,ikpt1,ikpt1i,ikpt2,ikpt2i ',ikpt_loc,ikpt1,ikpt1i,ikpt3(istep),ikpt3i(istep)
903 : ! call smatrix_k0_paw(atindx1,cprj_k,cprj_k,dtefield,ikpt1i,idir,ifor,&
904 : ! & mband,mpi_enreg,natom,ntypat,pawtab,smat_k_paw,typat)
905 169440 : if (calc_epaw3_force .or. calc_epaw3_stress) then
906 149568 : call dsdr_k_paw(cprj_k,cprj_kb,dsdr,dtefield,idir,ifor,mband,natom,ncpgr,typat)
907 : end if
908 : end if
909 :
910 270890 : icg1 = 0
911 270890 : icg = dtefield%cgindex(ikpt1i,isppol)
912 : ! DEBUG
913 : ! if(istep<=2)then
914 : ! if(ikpt1==1)then
915 : ! write(std_out,'(a,2i4,3e15.4)')'istep ikpt3, kpt, cgq', istep, ikpt3(istep), dtefield%fkptns(:,ikpt3(istep))
916 : ! write(std_out,*) cgq
917 : ! write(std_out,*)
918 : ! end if
919 : ! end if
920 : ! ENDDEBUG
921 : call smatrix(cg,cgq,cg1_k,ddkflag,dtm_k,icg,icg1,itrs_mult(istep),det_inv_smat,maxbd,&
922 : & mcg,count,mcg1_k,minbd,&
923 : & mpw,dtefield%mband_occ,dtefield%nband_occ(isppol),&
924 : & npw_k1,npw_k3(istep),my_nspinor,pwind_k_mult(:,istep),pwnsfac_k,sflag_k_mult(:,istep),&
925 270890 : & shiftbd,smat_inv,smat_k,smat_k_paw,usepaw)
926 :
927 : ! in finite electric field case with paw must save additional F3 term in forces
928 270890 : if(calc_epaw3_force) then
929 : ! when ncpgr = 3, gradients are wrt to atom displacements
930 : ! but when ncpgr = 9, first 6 gradients are wrt strains, last three are displacements
931 149568 : icpgr_offset = 0
932 149568 : if (ncpgr == 9) icpgr_offset = 6
933 448704 : do iatom = 1, natom
934 1346112 : do fdir = 1, 3
935 897408 : dsdr_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = zero
936 4786176 : do iband = 1, dtefield%nband_occ(isppol)
937 18845568 : do jband = 1, dtefield%nband_occ(isppol)
938 : ! collect Im{Trace{S^{-1}.dS/dR}} for this k point
939 : dsdr_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = &
940 : & dsdr_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) + &
941 : & smat_inv(2,iband,jband)*dsdr(1,iatom,icpgr_offset+fdir,jband,iband) + &
942 17948160 : & smat_inv(1,iband,jband)*dsdr(2,iatom,icpgr_offset+fdir,jband,iband)
943 : end do ! end sum over jband
944 : end do ! end sum over iband
945 : end do ! end sum over fdir
946 : end do ! end sum over iatom
947 : end if ! end check on calc_epaw3_force
948 :
949 : ! in finite electric field case with paw must save additional F3 term in stress
950 : ! note that when strains are present they are always saved before forces
951 : ! therefore no need for icpgr_offset in this case
952 270890 : if(calc_epaw3_stress) then
953 448704 : do iatom = 1, natom
954 2243520 : do fdir = 1, 6
955 1794816 : dsds_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = zero
956 9273216 : do iband = 1, dtefield%nband_occ(isppol)
957 37691136 : do jband = 1, dtefield%nband_occ(isppol)
958 : ! collect Im{Trace{S^{-1}.dS/de}} for this k point
959 : dsds_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = &
960 : & dsds_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) + &
961 : & smat_inv(2,iband,jband)*dsdr(1,iatom,fdir,jband,iband) + &
962 35896320 : & smat_inv(1,iband,jband)*dsdr(2,iatom,fdir,jband,iband)
963 : end do ! end sum over jband
964 : end do ! end sum over iband
965 : end do ! end sum over fdir
966 : end do ! end sum over iatom
967 : end if ! end check on calc_epaw3_stress
968 :
969 270890 : if ((det_inv_smat == 10).or.(det_inv_smat == 11)) then
970 :
971 264906 : if (sqrt(dtm_k(1)*dtm_k(1) + dtm_k(2)*dtm_k(2)) < tol12) then
972 : ! EB: the ABI_BUG has been replaced here by what is done in 67_common/m_cgwf.F90
973 : ! This avoid the code to stop for phonons under E-field too.
974 : ! TODO: Since the same is done in m_cgwf.F90 and in m_berryphase_new.F90,
975 : ! rationalization should be done with one single module.
976 : write(msg,'(a,i5,a,a,a)')&
977 0 : & ' For k-point #',ikpt1,',',ch10,&
978 0 : & ' the determinant of the overlap matrix is found to be 0. Fixing...'
979 : ! Try this:
980 0 : write(std_out,*)msg,dtm_k(1:2)
981 0 : if(abs(dtm_k(1))<=1d-12)dtm_k(1)=1d-12
982 0 : if(abs(dtm_k(2))<=1d-12)dtm_k(2)=1d-12
983 0 : write(std_out,*)' Changing to:',dtm_k(1:2)
984 : ! ABI_BUG(msg)
985 : end if
986 :
987 264906 : dtm_mult(1,ikpt1+(isppol-1)*dtefield%fnkpt,istep) = dtm_k(1)
988 264906 : dtm_mult(2,ikpt1+(isppol-1)*dtefield%fnkpt,istep) = dtm_k(2)
989 :
990 : end if
991 :
992 270890 : if ( efield_flag .and. inibz == 1 .and. istep == 1) then !!HONG
993 23124 : ikpt1i_sp=ikpt1i+(isppol-1)*dtset%nkpt
994 : dtefield%smat(:,:,:,ikpt1i_sp,ifor,idir) = &
995 5750172 : & smat_k(:,:,:)
996 : dtefield%sflag(:,ikpt1i_sp,ifor,idir) = &
997 158220 : & sflag_k_mult(:,1)
998 : end if
999 :
1000 : ! for IBZ k-points and first step, add
1001 270890 : if ((ddkflag==1 .and.((det_inv_smat == 1).or.(det_inv_smat == 11))) .and. inibz == 1 .and. istep == 1) then
1002 : cg1(:,icg + 1: icg + npw_k1*my_nspinor*nband_k) = &
1003 : cg1(:,icg + 1:icg + npw_k1*my_nspinor*nband_k) + &
1004 6652640 : dkinv*cg1_k(:,1:npw_k1*my_nspinor*nband_k)
1005 : end if
1006 :
1007 270890 : ABI_FREE(cgq)
1008 :
1009 : end if ! end if ikpt1 > 0 and isppol > 0
1010 :
1011 : end do ! end loop over istep
1012 :
1013 : ! if (ikpt_loc <= dtefield%fmkmem) sflag_k(:) = sflag_k_mult(:,1)
1014 1516817 : if (ikpt1 > 0) sflag_k(:) = sflag_k_mult(:,1)
1015 :
1016 : end do ! close loop over ikpt_loc (k-points, isppol)
1017 :
1018 827 : ABI_FREE(smat_inv)
1019 827 : ABI_FREE(smat_k)
1020 827 : ABI_FREE(smat_k_paw)
1021 1641 : if (calc_epaw3_force .or. calc_epaw3_stress) then
1022 234 : ABI_FREE(dsdr)
1023 : end if
1024 :
1025 : end do ! close loop over ifor
1026 :
1027 : ! MPI communicate stuff between everyone
1028 814 : if (nproc>1) then
1029 306 : count = 2*dtefield%fnkpt*nsppol*berrystep
1030 918 : ABI_MALLOC(buffer1,(count))
1031 612 : ABI_MALLOC(buffer2,(count))
1032 612 : buffer1(:) = reshape(dtm_mult(:,:,:),(/count/))
1033 306 : call xmpi_sum(buffer1,buffer2,count,spaceComm,ierr)
1034 1224 : dtm_mult(:,:,:) = reshape(buffer2(:),(/2,dtefield%fnkpt*nsppol,berrystep/))
1035 306 : ABI_FREE(buffer1)
1036 306 : ABI_FREE(buffer2)
1037 306 : if (calc_epaw3_force) then
1038 144 : count = natom*3*dtefield%fnkpt*nsppol
1039 432 : ABI_MALLOC(buffer1,(count))
1040 288 : ABI_MALLOC(buffer2,(count))
1041 288 : buffer1(:) = reshape(dsdr_sum(:,:,:),(/count/))
1042 144 : call xmpi_sum(buffer1,buffer2,count,spaceComm,ierr)
1043 576 : dsdr_sum(:,:,:) = reshape(buffer2(:),(/natom,3,dtefield%fnkpt*nsppol/))
1044 144 : ABI_FREE(buffer1)
1045 288 : ABI_FREE(buffer2)
1046 : end if
1047 612 : if (calc_epaw3_stress) then
1048 144 : count = natom*6*dtefield%fnkpt*nsppol
1049 432 : ABI_MALLOC(buffer1,(count))
1050 288 : ABI_MALLOC(buffer2,(count))
1051 288 : buffer1(:) = reshape(dsds_sum(:,:,:),(/count/))
1052 144 : call xmpi_sum(buffer1,buffer2,count,spaceComm,ierr)
1053 576 : dsds_sum(:,:,:) = reshape(buffer2(:),(/natom,6,dtefield%fnkpt*nsppol/))
1054 144 : ABI_FREE(buffer1)
1055 288 : ABI_FREE(buffer2)
1056 : end if
1057 : end if ! if parallel
1058 :
1059 : ! DEBUG
1060 : ! write(std_out,*)
1061 : ! write(std_out,*)'istep = 1, nsppol =',nsppol
1062 : ! istep=1
1063 : ! isppol=1
1064 : ! do jkpt = 1, dtefield%fnkpt
1065 : ! write(std_out,'(a,i4,3e15.4,2e15.4)')'jkpt, kpt, dtm_mult(:,kpt,1)', jkpt, dtefield%fkptns(:,jkpt), dtm_mult(:,jkpt+(isppol-1)*dtefield%fnkpt,istep)
1066 : ! end do
1067 : ! write(std_out,*)
1068 : ! write(std_out,*) "istep = 2"
1069 : ! if(berrystep>=2)then
1070 : ! istep=2
1071 : ! isppol=1
1072 : ! do jkpt = 1, dtefield%fnkpt
1073 : ! write(std_out,'(a,i4,3e15.4,2e15.4)')'jkpt, kpt, dtm_mult(:,kpt,2)', jkpt, dtefield%fkptns(:,jkpt), dtm_mult(:,jkpt+(isppol-1)*dtefield%fnkpt,istep)
1074 : ! end do
1075 : ! end if
1076 : ! ENDDEBUG
1077 :
1078 : ! ===========================================================================
1079 : ! in DDK case everything has been calculated above from finite difference
1080 : ! Now write the ddk WF to a file
1081 : ! ===========================================================================
1082 :
1083 814 : if (ddkflag == 1) then
1084 :
1085 13 : pertcase = idir + 3*natom
1086 13 : response = 1
1087 13 : call appdig(pertcase,dtfil%fnameabo_1wf,fiwf1o)
1088 39 : ABI_MALLOC(resid,(mband*nkpt*nsppol))
1089 7053 : resid(:) = zero
1090 :
1091 13 : call outresid(dtset,dtset%kptns,mband,dtset%nband,nkpt,nsppol,resid)
1092 :
1093 : call outwf(cg1,dtset,psps,eig_dum,fiwf1o,hdr,kg,dtset%kptns,&
1094 : & mband,mcg,mkmem,mpi_enreg,mpw,natom,dtset%nband,&
1095 : & nkpt,npwarr,nsppol,&
1096 13 : & occ_dum,response,dtfil%unwff2,wfs,wvl)
1097 :
1098 13 : if (save_cg13) then
1099 172803 : cg13(1:2,1:mcg13,idir) = cg1(1:2,1:mcg)
1100 : end if
1101 :
1102 13 : ABI_FREE(resid)
1103 : end if ! ddkflag == 1
1104 : ! end of ddk part for this idir
1105 :
1106 :
1107 : ! ===========================================================================
1108 : ! Compute the Berry phase polarization
1109 : ! ===========================================================================
1110 :
1111 1091 : if (polflag == 1) then
1112 :
1113 : ! Compute the electronic Berry phase
1114 :
1115 2436 : polb_mult(:,:)=zero
1116 1620 : do istep = 1,berrystep
1117 :
1118 816 : if(berrystep==1) then
1119 801 : write(msg,'(a,a)')ch10,' Compute the electronic contribution to polarization'
1120 801 : call wrtout(std_out,msg,'COLL')
1121 : else
1122 15 : write(msg,'(a,a,i4,a)')ch10,&
1123 30 : ' Compute the electronic contribution to polarization for a step of istep=',istep,'*dk'
1124 15 : call wrtout(std_out,msg,'COLL')
1125 : end if
1126 :
1127 816 : if(istep /= 1) then
1128 : ! construct the strings for a step of istep*dk
1129 : ! string length
1130 12 : istr=1
1131 12 : nkstr=1
1132 12 : ikpt1=1
1133 33 : do ikpt=1,dtefield%fnkpt
1134 153 : do jstep = 1,istep
1135 153 : ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir)
1136 : end do
1137 33 : if (ikpt1 == 1) exit
1138 33 : nkstr = nkstr + 1
1139 : end do
1140 : ! Check that the string length is a divisor of nkpt
1141 12 : if(mod(dtefield%fnkpt,nkstr) /= 0) then
1142 : write(msg,'(a,i5,a,i5,a,i7)')&
1143 0 : & ' For istep = ', istep,&
1144 0 : & ' The string length = ',nkstr,&
1145 0 : & ', is not a divisor of fnkpt =',dtefield%fnkpt
1146 0 : ABI_BUG(msg)
1147 : end if
1148 12 : nstr = dtefield%fnkpt/nkstr
1149 :
1150 : write(msg,'(a,i1,a,i2,a,i3,a,i6)')&
1151 12 : & ' berryphase_new: for direction ',idir, ' and istep ', istep, ', nkstr = ',nkstr,', nstr = ',nstr
1152 12 : call wrtout(std_out,msg,'COLL')
1153 12 : call wrtout(ab_out,msg,'COLL')
1154 :
1155 48 : ABI_MALLOC(idxkstr_mult,(nkstr,nstr))
1156 780 : iunmark = 1
1157 780 : kpt_mark(:)=0
1158 396 : do istr=1,nstr
1159 913 : do while(kpt_mark(iunmark) /= 0)
1160 529 : iunmark = iunmark + 1
1161 : end do
1162 384 : idxkstr_mult(1,istr) = iunmark
1163 384 : kpt_mark(iunmark)=1
1164 :
1165 384 : ikpt1 = idxkstr_mult(1,istr)
1166 780 : do jkstr=2, nkstr
1167 1728 : do jstep = 1, istep
1168 1728 : ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir)
1169 : end do
1170 384 : idxkstr_mult(jkstr,istr) = ikpt1
1171 768 : kpt_mark(ikpt1) = 1
1172 : end do
1173 : end do
1174 : else
1175 804 : nstr = dtefield%nstr(idir)
1176 804 : nkstr = dtefield%nkstr(idir)
1177 3216 : ABI_MALLOC(idxkstr_mult,(nkstr,nstr))
1178 562293 : idxkstr_mult(:,:) = dtefield%idxkstr(1:nkstr,1:nstr,idir)
1179 : end if
1180 : ! DEBUG
1181 : ! do istr=1,nstr
1182 : ! write(std_out,*)'string ', idxkstr_mult(:,istr)
1183 : ! end do
1184 : ! ENDBEBUG
1185 :
1186 2448 : ABI_MALLOC(det_string,(2,nstr))
1187 2448 : ABI_MALLOC(polberry,(nstr))
1188 816 : write(msg,'(a,10x,a,10x,a)')ch10,'istr','polberry(istr)'
1189 816 : call wrtout(std_out,msg,'COLL')
1190 :
1191 816 : polbtot = zero
1192 1632 : do isppol = 1, nsppol
1193 :
1194 141630 : det_string(1,:) = one ; det_string(2,:) = zero
1195 2448 : dtm_k(:) = one
1196 816 : det_average(:) = zero
1197 :
1198 :
1199 71223 : do istr = 1, nstr
1200 :
1201 1392327 : if(calc_epaw3_force) epawf3_str(:,:,:) = zero
1202 2714247 : if(calc_epaw3_stress) epaws3_str(:,:,:) = zero
1203 :
1204 562641 : do jkstr = 1, nkstr
1205 :
1206 492234 : ikpt=idxkstr_mult(jkstr,istr)
1207 :
1208 492234 : dtm_real=dtm_mult(1,ikpt+(isppol-1)*dtefield%fnkpt,istep)
1209 492234 : dtm_imag=dtm_mult(2,ikpt+(isppol-1)*dtefield%fnkpt,istep)
1210 :
1211 : dtm_k(1) = det_string(1,istr)*dtm_real - &
1212 492234 : & det_string(2,istr)*dtm_imag
1213 : dtm_k(2) = det_string(1,istr)*dtm_imag + &
1214 492234 : & det_string(2,istr)*dtm_real
1215 1476702 : det_string(1:2,istr) = dtm_k(1:2)
1216 : ! DEBUG
1217 : ! write(std_out,'(a,i4,3e15.4,2e15.4)')'ikpt, kpt, dtm', ikpt, dtefield%fkptns(:,ikpt), dtm_k
1218 : ! ENDDEBUG
1219 :
1220 492234 : if(calc_epaw3_force) then
1221 1038528 : do iatom = 1, natom
1222 3115584 : do fdir = 1, 3
1223 : epawf3_str(iatom,idir,fdir) = epawf3_str(iatom,idir,fdir) + &
1224 2769408 : & dsdr_sum(iatom,fdir,ikpt+(isppol-1)*dtefield%fnkpt)
1225 : end do ! end loop over fdir
1226 : end do ! end loop over natom
1227 : end if ! end check on calc_epaw3_force
1228 562641 : if(calc_epaw3_stress) then
1229 1038528 : do iatom = 1, natom
1230 5192640 : do fdir = 1, 6
1231 : epaws3_str(iatom,idir,fdir) = epaws3_str(iatom,idir,fdir) + &
1232 4846464 : & dsds_sum(iatom,fdir,ikpt+(isppol-1)*dtefield%fnkpt)
1233 : end do ! end loop over fdir
1234 : end do ! end loop over natom
1235 : end if ! end check on calc_epaw3_stress
1236 :
1237 : end do
1238 :
1239 70407 : if(calc_epaw3_force) then
1240 132192 : do iatom = 1, natom
1241 396576 : do fdir = 1, 3
1242 : dtefield%epawf3(iatom,idir,fdir) = dtefield%epawf3(iatom,idir,fdir) + &
1243 352512 : & epawf3_str(iatom,idir,fdir)
1244 : end do ! end loop over fdir
1245 : end do ! end loop over natom
1246 : end if ! end check on calc_epaw3_force
1247 70407 : if(calc_epaw3_stress) then
1248 132192 : do iatom = 1, natom
1249 660960 : do fdir = 1, 6
1250 : dtefield%epaws3(iatom,idir,fdir) = dtefield%epaws3(iatom,idir,fdir) + &
1251 616896 : & epaws3_str(iatom,idir,fdir)
1252 : end do ! end loop over fdir
1253 : end do ! end loop over natom
1254 : end if ! end check on calc_epaw3_stress
1255 :
1256 : det_average(:) = det_average(:) + &
1257 212037 : & det_string(:,istr)/dble(nstr)
1258 :
1259 : end do
1260 :
1261 :
1262 : ! correction to obtain a smooth logarithm of the determinant
1263 1632 : ABI_MALLOC(str_flag,(nstr))
1264 : ! DEBUG
1265 : ! since we don't have any case of non-nul Chern number,
1266 : ! we must change the det_string value "by brute force" if we want debug this
1267 : ! allocate(det_string_test(2,dtefield%nstr(idir)))
1268 : ! det_string_test(:,:)=det_string(:,:)
1269 : ! kk=0
1270 : ! det_string(1,1)=cos(2._dp*Pi*real(kk,dp)/four)
1271 : ! det_string(2,1)=sin(2._dp*Pi*real(kk,dp)/four)
1272 : ! jj=dtefield%str_neigh(1,1,idir)
1273 : ! ll=dtefield%str_neigh(2,1,idir)
1274 : ! do while (jj/=1)
1275 : ! kk=kk+1
1276 : ! det_string(1,jj)=cos(2._dp*Pi*real(kk,dp)/four)
1277 : ! det_string(2,jj)=sin(2._dp*Pi*real(kk,dp)/four)
1278 : ! det_string(1,ll)=cos(-2._dp*Pi*real(kk,dp)/four)
1279 : ! det_string(2,ll)=sin(-2._dp*Pi*real(kk,dp)/four)
1280 : ! jj=dtefield%str_neigh(1,jj,idir)
1281 : ! ll=dtefield%str_neigh(2,ll,idir)
1282 : ! enddo
1283 : ! ENDDEBUG
1284 816 : if (istep==1) then
1285 2412 : do ineigh_str = 1,2
1286 141654 : str_flag(:)=0
1287 : delta_str(:) = &
1288 1608 : & dtefield%coord_str(:,dtefield%str_neigh(ineigh_str,1,idir),idir) - dtefield%coord_str(:,1,idir)
1289 : dstr(:)= delta_str(:) - nint(delta_str(:)) - real(dtefield%strg_neigh(ineigh_str,1,:,idir),dp)
1290 : dist_=0._dp
1291 : do kk = 1,2
1292 : do jj = 1,2
1293 : dist_ = dist_ + dstr(kk)*dtefield%gmet_str(kk,jj,idir)*dstr(jj)
1294 : end do
1295 : end do
1296 : dist_=sqrt(dist_)
1297 142458 : do istr = 1,dtefield%nstr(idir)
1298 141654 : if(str_flag(istr)==0)then
1299 : ! write(std_out,*)'new string'
1300 12339 : str_flag(istr)=1
1301 12339 : call rhophi(det_string(:,istr),dphase,rho)
1302 : ! write(std_out,'(i4,e15.4,e15.4,e15.4)')istr, det_string(:,istr),dphase
1303 12339 : dphase_init=dphase
1304 12339 : jstr = dtefield%str_neigh(ineigh_str,istr,idir)
1305 140046 : do while (istr/=jstr)
1306 127707 : str_flag(jstr)=1
1307 127707 : call rhophi(det_string(:,jstr),dphase_new,rho)
1308 127707 : jj=nint((dphase_new-dphase)/(2._dp*Pi))
1309 : ! DEBUG
1310 : ! write(std_out,'(i4,e15.4,e15.4,e15.4,e15.4,i4)')jstr, det_string(:,jstr),dphase_new,dphase_new-dphase,jj
1311 : ! ENDDEBUG
1312 127707 : dphase_new=dphase_new-two*Pi*real(jj,dp)
1313 127707 : if(jj/=0)then
1314 840 : write(msg,'(6a)') ch10,&
1315 840 : & ' berryphase_new : WARNING -',ch10,&
1316 840 : & ' the berry phase has some huge variation in the space of strings of k-points',ch10,&
1317 1680 : & ' ABINIT is trying to correct the berry phase, but it is highly experimental'
1318 840 : call wrtout(std_out,msg,'PERS')
1319 : end if
1320 : ! if(jj/=0)write(std_out,'(i4,e15.4,e15.4,e15.4,e15.4)')jstr, det_string(:,jstr),dphase_new,dphase_new-dphase
1321 127707 : dphase=dphase_new
1322 127707 : jstr=dtefield%str_neigh(ineigh_str,jstr,idir)
1323 : end do
1324 : ! write(std_out,*)dphase_init, dphase, (dphase-dphase_init)/(2._dp*Pi),nint((dphase-dphase_init)/(2._dp*Pi))
1325 : end if
1326 : end do
1327 : end do
1328 : end if
1329 816 : ABI_FREE(str_flag)
1330 : ! DEBUG
1331 : ! deallocate(dist_str)
1332 : ! det_string(:,:)=det_string_test(:,:)
1333 : ! deallocate(det_string_test)
1334 : ! ENDDEBUG
1335 :
1336 : ! First berry phase that corresponds to det_average
1337 : ! phase0 = atan2(det_average(2),det_average(1))
1338 816 : call rhophi(det_average,phase0,rho)
1339 816 : det_mod = det_average(1)**2+det_average(2)**2
1340 :
1341 : ! Then berry phase that corresponds to each string relative to the average
1342 71223 : do istr = 1, nstr
1343 :
1344 : rel_string(1) = (det_string(1,istr)*det_average(1) + &
1345 70407 : det_string(2,istr)*det_average(2))/det_mod
1346 : rel_string(2) = (det_string(2,istr)*det_average(1) - &
1347 70407 : det_string(1,istr)*det_average(2))/det_mod
1348 : ! dphase = atan2(rel_string(2),rel_string(1))
1349 70407 : call rhophi(rel_string,dphase,rho)
1350 70407 : polberry(istr) = dtefield%sdeg*(phase0 + dphase)/two_pi
1351 70407 : polb_mult(isppol,istep) = polb_mult(isppol,istep) + polberry(istr)/(istep*dtefield%nstr(idir))
1352 70407 : polb(isppol) = zero
1353 141774 : do jstep=1, istep
1354 141774 : polb(isppol)=polb(isppol)+coef(jstep,istep)*polb_mult(isppol,jstep)
1355 : end do
1356 :
1357 70407 : write(msg,'(10x,i6,7x,e16.9)')istr,polberry(istr)
1358 141630 : call wrtout(std_out,msg,'COLL')
1359 :
1360 : end do
1361 :
1362 816 : if(berrystep>1)then
1363 : write(msg,'(9x,a,7x,e16.9,1x,a,i4,a,i4,a)')&
1364 15 : & 'total',polb_mult(isppol,istep),'(isppol=',isppol,', istep=',istep,')'!,ch10
1365 15 : call wrtout(std_out,msg,'COLL')
1366 :
1367 : write(msg,'(3x,a,7x,e16.9,1x,a,i4,a,i4,a,a)')&
1368 15 : & '+correction',polb(isppol),'(isppol=',isppol,', istep=1..',istep,')',ch10
1369 15 : call wrtout(std_out,msg,'COLL')
1370 :
1371 : else
1372 :
1373 : write(msg,'(9x,a,7x,e16.9,1x,a,i4,a)')&
1374 801 : & 'total',polb_mult(isppol,istep),'(isppol=',isppol,')'!,ch10
1375 801 : call wrtout(std_out,msg,'COLL')
1376 : end if
1377 :
1378 2448 : polbtot = polbtot + polb(isppol)
1379 :
1380 : end do ! isppol
1381 :
1382 : ! Fold into interval [-1,1]
1383 816 : polbtot = polbtot - 2_dp*nint(polbtot/2_dp)
1384 :
1385 816 : ABI_FREE(det_string)
1386 816 : ABI_FREE(polberry)
1387 :
1388 : ! ==========================================================================
1389 :
1390 : ! Compute the ionic Berry phase
1391 :
1392 816 : call xred2xcart(natom,rprimd,xcart,xred)
1393 816 : politot = zero
1394 816 : write(msg,'(a)')' Compute the ionic contributions'
1395 816 : call wrtout(std_out,msg,'COLL')
1396 :
1397 816 : write(msg,'(a,2x,a,2x,a,15x,a)')ch10,'itom', 'itypat', 'polion'
1398 816 : call wrtout(std_out,msg,'COLL')
1399 :
1400 3804 : do iatom = 1, natom
1401 2988 : itypat = typat(iatom)
1402 :
1403 : ! The ionic phase can be computed much easier
1404 2988 : polion = zion(itypat)*xred(idir,iatom)
1405 :
1406 : ! Fold into interval (-1,1)
1407 2988 : polion = polion - 2_dp*nint(polion/2_dp)
1408 2988 : politot = politot + polion
1409 2988 : write(msg,'(2x,i2,5x,i2,10x,e16.9)') iatom,itypat,polion
1410 3804 : call wrtout(std_out,msg,'COLL')
1411 : end do
1412 :
1413 : ! Fold into interval [-1,1] again
1414 816 : politot = politot - 2_dp*nint(politot/2_dp)
1415 816 : pion(idir) = politot
1416 :
1417 816 : write(msg,'(9x,a,7x,es19.9)') 'total',politot
1418 816 : call wrtout(std_out,msg,'COLL')
1419 :
1420 :
1421 : ! ==========================================================================
1422 :
1423 : ! Compute the total polarization
1424 :
1425 816 : poltot = politot + polbtot
1426 :
1427 816 : if (berrystep==1)then
1428 801 : write(msg,'(a,a)')ch10,' Summary of the results'
1429 801 : call wrtout(std_out,msg,'COLL')
1430 801 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
1431 : else
1432 15 : write(msg,'(a,a,i4)')ch10,' Summary of the results for istep =',istep
1433 15 : call wrtout(std_out,msg,'COLL')
1434 15 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
1435 : end if
1436 :
1437 816 : write(msg,'(a,es19.9)')' Electronic Berry phase ' ,polbtot
1438 816 : call wrtout(std_out,msg,'COLL')
1439 816 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
1440 :
1441 816 : write(msg,'(a,es19.9)')' Ionic phase ', politot
1442 816 : call wrtout(std_out,msg,'COLL')
1443 816 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
1444 :
1445 816 : write(msg,'(a,es19.9)')' Total phase ', poltot
1446 816 : call wrtout(std_out,msg,'COLL')
1447 816 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
1448 :
1449 : ! REC start
1450 816 : if(abs(dtset%polcen(idir))>tol8)then
1451 0 : poltot = poltot-dtset%polcen(idir)
1452 : write(msg,'(a,f15.10)') &
1453 0 : & ' Translating Polarization by P0 for centrosymmetric cell: ',dtset%polcen(idir)
1454 0 : call wrtout(std_out,msg,'COLL')
1455 0 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
1456 : end if
1457 : ! REC end
1458 :
1459 816 : poltot = poltot - 2.0_dp*nint(poltot/2._dp)
1460 816 : write(msg,'(a,es19.9)')' Remapping in [-1,1] ', poltot
1461 816 : call wrtout(std_out,msg,'COLL')
1462 816 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
1463 :
1464 : ! ! REC and HONG
1465 : ! =====================================================================================
1466 : ! Polarization branch control (start)
1467 : ! -------------------------------------------------------------------------------------
1468 : ! berrysav == 0, for non fixed D/d calculation, polarizaion is in [-1,1],done above
1469 : ! for fixed D/d calculation, choose polarization to minimize internal
1470 : ! energy, or minimize |red_efiled|. (red_dfield=red_efiled+red_ptot)
1471 : ! (d=e+p, as (26) of Stengel, Suppl.) [[cite:Stengel2009]]
1472 : ! This is default value.
1473 : !
1474 : ! berrysav == 1, keep the polarization on the same branch, which saved in file POLSAVE
1475 : ! ======================================================================================
1476 :
1477 : ! for fixed D/d calculation, choose polarization to minimize internal energy, or to minimize reduced electric field |red_efield|
1478 816 : if(dtset%berrysav ==0 .and. (dtset%berryopt == 6 .or. dtset%berryopt == 7 .or. &
1479 : & dtset%berryopt == 16 .or. dtset%berryopt == 17)) then
1480 :
1481 0 : jump=-nint(dtset%red_dfield(idir) - poltot) ! red_efield = red_dfield - poltot
1482 :
1483 0 : if(jump /= 0)then
1484 : write(msg,'(a,i1,a,es19.9,a,i2)') &
1485 0 : & ' P(',idir,') Shifted polarization branch to minimize red_efield &
1486 0 : & k from ',poltot, ' by ',jump
1487 0 : call wrtout(std_out,msg,'COLL')
1488 0 : if (unit_out /= 0) then
1489 0 : call wrtout(unit_out,msg,'COLL')
1490 : end if
1491 0 : poltot=poltot-jump
1492 : end if
1493 0 : pol0(idir)=poltot
1494 :
1495 : end if
1496 :
1497 :
1498 : ! keep the polarization on the same branch.
1499 816 : if (dtset%berrysav == 1) then
1500 :
1501 : ! use saved polarization to keep on same branch
1502 0 : inquire(file='POLSAVE',exist=lexist)
1503 0 : if(lexist)then
1504 0 : if(idir==1)then
1505 0 : if(mpi_enreg%me==0)then
1506 0 : if (open_file('POLSAVE',msg,newunit=unt,status='OLD') /= 0) then
1507 0 : ABI_ERROR(msg)
1508 : end if
1509 0 : read(unt,*)pol0
1510 0 : write(msg,'(a,3f20.12)')'Reading old polarization:',pol0
1511 0 : call wrtout(std_out,msg,'COLL')
1512 0 : if (unit_out /= 0) then
1513 0 : call wrtout(unit_out,msg,'COLL')
1514 : end if
1515 0 : close(unt)
1516 : end if
1517 0 : call xmpi_bcast(pol0,0,spaceComm,ierr)
1518 : end if
1519 : else
1520 0 : pol0(idir)=poltot
1521 : end if
1522 0 : jump=nint(poltot-pol0(idir))
1523 0 : if(jump /= 0)then
1524 : write(msg,'(a,i1,a,es19.9,a,i2)') &
1525 0 : & ' P(',idir,') jumped to new branch. Shifting bac&
1526 0 : & k from ',poltot, ' by ',jump
1527 0 : call wrtout(std_out,msg,'COLL')
1528 0 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
1529 0 : poltot=poltot-jump
1530 : end if
1531 :
1532 0 : pol0(idir)=poltot
1533 :
1534 : end if
1535 :
1536 : ! =====================================================================================
1537 : ! Polarization branch control (end)
1538 : ! =====================================================================================
1539 :
1540 :
1541 : ! Transform the phase into a polarization
1542 816 : fac = 1._dp/(gmod*dtefield%nkstr(idir))
1543 : ! !REC fac = fac/ucvol
1544 : ! !REC pol = fac*poltot
1545 816 : red_ptot(idir)=poltot !!REC
1546 816 : pol = fac*red_ptot(idir)/ucvol !!REC
1547 816 : ptot(idir)=red_ptot(idir)/ucvol !!REC
1548 816 : write(msg,'(a,a,es19.9,a,a,a,es19.9,a,a)')ch10,&
1549 816 : & ' Polarization ', pol,' (a.u. of charge)/bohr^2',ch10,&
1550 816 : & ' Polarization ', pol*(e_Cb)/(Bohr_Ang*1d-10)**2,&
1551 1632 : & ' C/m^2',ch10
1552 816 : call wrtout(std_out,msg,'COLL')
1553 816 : if (unit_out /= 0) call wrtout(unit_out,msg,'COLL')
1554 :
1555 1620 : ABI_FREE(idxkstr_mult)
1556 :
1557 : end do !istep
1558 804 : pel(idir) = polbtot
1559 : end if ! if calculate polarization polflag==1
1560 : end do ! Close loop over idir
1561 :
1562 :
1563 : !!REC start
1564 277 : if (dtset%berrysav == 1) then
1565 0 : if(mpi_enreg%me==0)then
1566 0 : if (open_file('POLSAVE',msg,newunit=unt,status='UNKNOWN') /= 0) then
1567 0 : ABI_ERROR(msg)
1568 : end if
1569 0 : write(unt,'(3F20.12)') pol0
1570 0 : close(unt)
1571 : end if
1572 0 : first=.false.
1573 : end if
1574 : !!REC end
1575 :
1576 : !-------------------------------------------------
1577 : ! Compute polarization in cartesian coordinates
1578 : !-------------------------------------------------
1579 : !if (all(dtset%rfdir(:) == 1) .AND. (calc_pol_ddk .NE. 2)) then
1580 1086 : if (all(dtset%rfdir(:) == 1)) then
1581 :
1582 267 : if(usepaw.ne.1) then
1583 168 : pelev=zero
1584 : else
1585 99 : call pawpolev(my_natom,natom,ntypat,pawrhoij,pawtab,pelev,comm_atom=mpi_enreg%comm_atom)
1586 : ! note that in the PAW case, the pelev contribution is already
1587 : ! implicitly included in the electronic polarization, from the
1588 : ! discretized derivative operator. In the NCPP case no such
1589 : ! terms exist anyway. Actually in the PAW formulation
1590 : ! such terms are included to all orders, unlike in USPP where only
1591 : ! zeroth and first-order terms are. In USPP the first-order term
1592 : ! is pelev. Here we compute pelev separately only for reporting
1593 : ! purposes in polcart, it is not added into pel or used in the the
1594 : ! PAW finite field code in make_grad_berry.F90
1595 : ! 13 June 2012 J Zwanziger
1596 : end if
1597 267 : call polcart(red_ptot,pel,pel_cart,pelev,pion,pion_cart,3,ptot_cart,rprimd,ucvol,unit_out,usepaw)
1598 : end if
1599 :
1600 277 : ABI_FREE(pwnsfac_k)
1601 277 : ABI_FREE(sflag_k)
1602 277 : ABI_FREE(cg1_k)
1603 277 : if (ddkflag == 1) then
1604 5 : ABI_FREE(cg1)
1605 5 : ABI_FREE(eig_dum)
1606 5 : ABI_FREE(occ_dum)
1607 : end if
1608 :
1609 277 : if (usepaw == 1) then
1610 99 : ABI_FREE(dimlmn)
1611 99 : call pawcprj_free(cprj_k)
1612 99 : call pawcprj_free(cprj_kb)
1613 99 : call pawcprj_free(cprj_gat)
1614 891 : ABI_FREE(cprj_k)
1615 891 : ABI_FREE(cprj_kb)
1616 1899 : ABI_FREE(cprj_gat)
1617 99 : if (dtset%kptopt /= 3) then
1618 95 : call pawcprj_free(cprj_ikn)
1619 95 : call pawcprj_free(cprj_fkn)
1620 855 : ABI_FREE(cprj_ikn)
1621 855 : ABI_FREE(cprj_fkn)
1622 : end if
1623 99 : if (calc_epaw3_force) then
1624 78 : ABI_FREE(dsdr_sum)
1625 78 : ABI_FREE(epawf3_str)
1626 : end if
1627 99 : if (calc_epaw3_stress) then
1628 78 : ABI_FREE(dsds_sum)
1629 78 : ABI_FREE(epaws3_str)
1630 : end if
1631 :
1632 99 : if (nproc>1) then
1633 54 : call pawcprj_free(cprj_buf)
1634 486 : ABI_FREE(cprj_buf)
1635 : end if
1636 : end if
1637 :
1638 277 : ABI_FREE(ikpt3)
1639 277 : ABI_FREE(ikpt3i)
1640 277 : ABI_FREE(sflag_k_mult)
1641 277 : ABI_FREE(npw_k3)
1642 277 : ABI_FREE(pwind_k_mult)
1643 277 : ABI_FREE(itrs_mult)
1644 277 : ABI_FREE(coef)
1645 277 : ABI_FREE(polb_mult)
1646 277 : ABI_FREE(dtm_mult)
1647 :
1648 : !write(std_out,*)'berryphase_new exit'
1649 :
1650 277 : end subroutine berryphase_new
1651 : !!***
1652 :
1653 : !!****f* ABINIT/update_e_field_vars
1654 : !! NAME
1655 : !! update_e_field_vars
1656 : !!
1657 : !! FUNCTION
1658 : !! This routine updates E field variables
1659 : !!
1660 : !! INPUTS
1661 : !! atindx(natom)=index table for atoms, inverse of atindx (see gstate.f)
1662 : !! atindx1(natom)=index table for atoms (see gstate.f)
1663 : !! cg(2,mcg)=planewave coefficients of wavefunctions
1664 : !! dimcprj(usepaw*natom)=lmn_size for each atom
1665 : !! dtfil <type(datafiles_type)>=variables related to files
1666 : !! gmet(3,3)=metric in reciprocal space
1667 : !! gprimd(3,3)=reciprocal space dimensional primitive translations
1668 : !! idir = determines directions for derivatives computed in ctocprj (0 for all)
1669 : !! kg(3,mpw*mkmem)=reduced planewave coordinates
1670 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
1671 : !! mkmem=number of k points treated by this node.
1672 : !! mpw=maximum dimensioned size of npw
1673 : !! my_natom=number of atoms treated by current processor
1674 : !! natom=number of atoms in cell
1675 : !! nattyp(ntypat)=number of atoms of each type
1676 : !! ngfft(18)=contain all needed information about 3D FFT, see ~ABINIT/Infos/vargs.htm#ngfft
1677 : !! nkpt=number of k-points
1678 : !! npwarr(nkpt)=number of planewaves in basis at this k point
1679 : !! ntypat=number of types of atoms in unit cell
1680 : !! pawrhoij(natom*usepaw) <type(pawrhoij_type)> atomic occupancies
1681 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
1682 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
1683 : !! pwind(pwind_alloc,2,3) = array used to compute
1684 : !! the overlap matrix smat between k-points (see initberry.f)
1685 : !! pwind_alloc = first dimension of pwind
1686 : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
1687 : !! rmet(3,3)=metric in real space
1688 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1689 : !! scfcv_level= 0 if calling before scf loop, 1 if during
1690 : !! scfcv_quit=signals whether calling during scf quit (see scfcv.F90)
1691 : !! scfcv_step=istep value of loop counter from scfcv.F90
1692 : !! ucvol=unit cell volume in bohr**3.
1693 : !! unit_out= unit for output of the results (usually the .out file of ABINIT)
1694 : !! The option unit_out = 0 is allowed. In this case, no information is written
1695 : !! to the output file but only to the log file.
1696 : !! usepaw= 1: use paw framework. 0:do not use paw.
1697 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for
1698 : !! each G and k point
1699 : !! ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real
1700 : !! spherical harmonics
1701 : !!
1702 : !! OUTPUT
1703 : !! efield_old_cart(3)=updating cartesian values of efield (used in berryopt 6,16,17)
1704 : !! pel_cg(3)=electronic polarization
1705 : !! pelev(3)=leading order PAW contribution in pel_cg (for reporting purposes only)
1706 : !! pion(3)=ionic part of polarization
1707 : !! ptot(3)=total polarization
1708 : !! red_efield2=updating efield used in berryopt 16,17
1709 : !! red_efield2_old=updating efield used in berryopt 16.17
1710 : !! red_ptot=updating efield used in berryopt 16.17
1711 : !!
1712 : !! SIDE EFFECTS
1713 : !! Input/Output
1714 : !! dtset <type(dataset_type)>=all input variables in this dataset
1715 : !! dtefield <type(efield_type)> = efield variables
1716 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
1717 : !! mpi_enreg=information about MPI parallelization
1718 : !! ptot_cart(3)=total polarization in cartesian coordinates
1719 : !! xred(3,natom)=reduced atomic coordinates
1720 : !!
1721 : !! SOURCE
1722 :
1723 52323 : subroutine update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,&
1724 52323 : & efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,&
1725 52323 : & mkmem,mpi_enreg,mpw,my_natom,natom,nattyp,ngfft,nkpt,npwarr,ntypat,&
1726 52323 : & pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,pwind,&
1727 52323 : & pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,&
1728 : & scfcv_level,scfcv_quit,scfcv_step,ucvol,unit_out,&
1729 52323 : & usepaw,xred,ylm,ylmgr)
1730 :
1731 : !Arguments ------------------------------------
1732 : integer, intent(in) :: idir,mcg,mkmem,mpw,my_natom,natom,nkpt,ntypat
1733 : integer, intent(in) :: pwind_alloc,scfcv_level,scfcv_quit,scfcv_step,unit_out,usepaw
1734 : real(dp), intent(in) :: ucvol
1735 : type(datafiles_type), intent(in) :: dtfil
1736 : type(pseudopotential_type),intent(in) :: psps
1737 : type(dataset_type), intent(inout) :: dtset
1738 : type(efield_type), intent(inout) :: dtefield
1739 : type(hdr_type), intent(inout) :: hdr
1740 : type(MPI_type), intent(inout) :: mpi_enreg
1741 : !arrays
1742 : integer, intent(in) :: atindx(natom),atindx1(natom),dimcprj(usepaw*natom)
1743 : integer, intent(in) :: kg(3,mpw*mkmem),nattyp(ntypat)
1744 : integer, intent(in) :: ngfft(18),npwarr(nkpt),pwind(pwind_alloc,2,3)
1745 : real(dp), intent(in) :: cg(2,mcg),gmet(3,3),gprimd(3,3)
1746 : real(dp), intent(in) :: pwnsfac(2,pwind_alloc)
1747 : real(dp), intent(in) :: rmet(3,3),rprimd(3,3)
1748 : real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
1749 : real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
1750 : real(dp), intent(inout) :: ptot_cart(3),xred(3,natom),efield_old_cart(3) !vz_i
1751 : real(dp), intent(out) :: pel_cg(3),pelev(3),pion(3) !vz_i
1752 : real(dp), intent(inout) :: red_efield2(3),red_efield2_old(3) !vz_i
1753 : real(dp), intent(out) :: ptot(3),red_ptot(3) !vz_i
1754 : type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*usepaw)
1755 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
1756 :
1757 : !Local variables -------------------------
1758 : !scalars
1759 : character(len=500) :: msg
1760 : integer :: ctocprj_choice,iatom,ii,iorder_cprj,mcg13,mcprj,my_nspinor,ncpgr
1761 : integer :: optberry,usecprj
1762 : logical :: calc_epaw3_force, calc_epaw3_stress, efield, save_cg13
1763 : !arrays
1764 : real(dp) :: efield_test_cart(3),red_efield1(3)
1765 52323 : real(dp),allocatable :: cg13(:,:,:),ph1d(:,:)
1766 52323 : type(pawcprj_type),allocatable :: cprj(:,:)
1767 :
1768 : ! *************************************************************************
1769 :
1770 52323 : efield = .false.
1771 :
1772 : if ( dtset%berryopt == 4 .or. &
1773 : & dtset%berryopt == 6 .or. &
1774 : & dtset%berryopt == 7 .or. &
1775 : & dtset%berryopt ==14 .or. &
1776 52323 : & dtset%berryopt ==16 .or. &
1777 : & dtset%berryopt ==17 ) efield = .true.
1778 191 : calc_epaw3_force = ( efield .and. dtset%optforces /= 0 .and. usepaw == 1 )
1779 52323 : calc_epaw3_stress = ( efield .and. dtset%optstress /= 0 .and. usepaw == 1 )
1780 :
1781 52323 : usecprj=1; if (psps%usepaw==0) usecprj = 0
1782 52323 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
1783 52323 : mcprj=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol/mpi_enreg%nproc_band
1784 : !Just to be sure
1785 52323 : if (efield.and.mpi_enreg%nproc_band/=1) then
1786 0 : msg='Berry phase calculations are not designed for band parallelization!'
1787 0 : ABI_BUG(msg)
1788 : end if
1789 :
1790 52323 : ncpgr = 0
1791 52323 : ctocprj_choice = 1 ! no derivs
1792 52323 : if ( efield .and. psps%usepaw == 1) then
1793 52884 : ABI_MALLOC(cprj,(dtset%natom,mcprj))
1794 : ! finite electric field may need gradients for forces, stress
1795 69 : if (calc_epaw3_force .and. .not. calc_epaw3_stress) then
1796 0 : ncpgr = 3; ctocprj_choice = 2 ! derivs w.r.t. position
1797 69 : else if (.not. calc_epaw3_force .and. calc_epaw3_stress) then
1798 0 : ncpgr = 6; ctocprj_choice = 3 ! derivs w.r.t strain
1799 69 : else if (calc_epaw3_force .and. calc_epaw3_stress) then
1800 69 : ncpgr = 9; ctocprj_choice = 23 ! derivs w.r.t. position and strain
1801 : end if
1802 69 : call pawcprj_alloc(cprj,ncpgr,dimcprj)
1803 69 : iatom=0 ; iorder_cprj=1 ! retain ordering of input list
1804 : ! all arguments to ctocprj are defined already except ph1d, do that here
1805 207 : ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
1806 69 : call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
1807 : call ctocprj(atindx,cg,ctocprj_choice,cprj,gmet,gprimd,iatom,idir,iorder_cprj,&
1808 : & dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,&
1809 : & mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,nattyp,dtset%nband,&
1810 : & dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,&
1811 : & dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,&
1812 69 : & dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr)
1813 69 : ABI_FREE(ph1d)
1814 : else
1815 52254 : ABI_MALLOC(cprj,(0,0))
1816 : end if ! end update of cprj
1817 :
1818 52323 : if ( efield ) then ! compute polarization and if necessary store cprj in efield
1819 191 : optberry=1
1820 191 : pel_cg(:) = zero;pelev=zero
1821 :
1822 191 : mcg13 = 0
1823 191 : save_cg13 = .FALSE.
1824 191 : ABI_MALLOC(cg13,(2,mcg13,3))
1825 : call berryphase_new(atindx1,cg,cg13,cprj,dtefield,dtfil,dtset,psps,gprimd,hdr,psps%indlmn,kg,&
1826 : & psps%lmnmax,dtset%mband,mcg,mcg13,mcprj,dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,&
1827 : & dtset%natom,npwarr,dtset%nsppol,psps%ntypat,dtset%nkpt,optberry,pawrhoij,pawtab,&
1828 : & pel_cg,pelev,pion,ptot,red_ptot,pwind,&
1829 : & pwind_alloc,pwnsfac,rprimd,save_cg13,dtset%typat,ucvol,&
1830 191 : & unit_out,usecprj,psps%usepaw,xred,psps%ziontypat)
1831 191 : ABI_FREE(cg13)
1832 :
1833 764 : dtefield%red_ptot1(:)=red_ptot(:)
1834 :
1835 : end if ! end compute polarization and store cprj for efield
1836 :
1837 191 : if (efield .and. (scfcv_level == 0) ) then ! do this before scfcv loop
1838 :
1839 100 : efield_old_cart(:)=dtset%efield(:) !!HONG
1840 :
1841 : ! save this value in order to print the final value of real electric field, comparing with the desired red_fieldbar
1842 100 : dtefield%efield2(:)=dtset%efield(:)
1843 :
1844 25 : if ( dtset%berryopt ==16 .or. dtset%berryopt ==17) then !!HONG
1845 0 : do ii=1,3
1846 0 : red_efield2(ii)=zero
1847 0 : red_efield2_old(ii) =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii))
1848 : end do
1849 : end if
1850 :
1851 25 : if (dtset%berryopt == 14 .and. scfcv_quit /=1) then
1852 : ! ! Convert polarization to cartesian coords
1853 :
1854 0 : ptot_cart(:)=zero
1855 0 : do ii = 1,3
1856 : ptot_cart(ii)=rprimd(ii,1)*red_ptot(1) + rprimd(ii,2)*red_ptot(2) + &
1857 0 : & rprimd(ii,3)*red_ptot(3)
1858 : end do
1859 0 : ptot_cart(:)=ptot_cart(:)/ucvol
1860 :
1861 0 : do ii=1,3
1862 0 : dtefield%efield_dot(ii) = dot_product(dtset%efield(:),rprimd(:,ii))
1863 : end do
1864 :
1865 : ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar
1866 0 : write(msg,'(a,a)') ch10, 'scfcv: Constant reduced ebar-field:'
1867 :
1868 0 : call wrtout(std_out,msg,'COLL')
1869 0 : call prtefield(dtset,dtefield,std_out,rprimd)
1870 :
1871 0 : if(dtset%prtvol>=10)then
1872 0 : call wrtout(ab_out,msg,'COLL')
1873 0 : call prtefield(dtset,dtefield,ab_out,rprimd)
1874 : end if
1875 :
1876 : ! updating E field
1877 0 : do ii =1,3 ! desired E field
1878 : efield_test_cart(ii)=gprimd(ii,1)*dtset%red_efieldbar(1) + &
1879 0 : & gprimd(ii,2)*dtset%red_efieldbar(2)+gprimd(ii,3)*dtset%red_efieldbar(3)
1880 : end do
1881 :
1882 : ! if not convergence well, need to add some code here to make sure efield_test_cart(:) not change much
1883 0 : dtset%efield(:) = efield_test_cart(:)
1884 :
1885 : end if ! berryopt ==14
1886 :
1887 : end if ! end efield .and. scfcv_level 0 tasks
1888 :
1889 : !!!
1890 : !!! Various printing and update steps for the different efield options
1891 : !!!
1892 :
1893 191 : if (efield .and. (scfcv_level == 1) ) then ! do this each scf step
1894 :
1895 166 : if (dtset%prtvol >= 10)then
1896 0 : write(msg,'(6(a),3(e16.9,2x),a,a,3(e16.9,2x))')ch10,&
1897 0 : & ' scfcv: New value of the polarization:',ch10,&
1898 0 : & ' (reduced coordinates, a. u.)',ch10,&
1899 0 : & ' Electronic berry phase: ', (pel_cg(ii), ii = 1, 3)
1900 0 : call wrtout(ab_out,msg,'COLL')
1901 0 : call wrtout(std_out,msg,'COLL')
1902 0 : if(psps%usepaw==1) then
1903 : write(msg,'(a,3(e16.9,2x))')&
1904 0 : & ' ...includes PAW on-site term: ', (pelev(ii), ii = 1, 3)
1905 0 : call wrtout(ab_out,msg,'COLL')
1906 0 : call wrtout(std_out,msg,'COLL')
1907 : end if
1908 : write(msg,'(a,3(e16.9,2x),a,a,3(e16.9,2x))')&
1909 0 : & ' Ionic: ', (pion(ii), ii = 1, 3), ch10, &
1910 0 : & ' Total: ', (red_ptot(ii), ii = 1, 3) !!REC
1911 0 : call wrtout(ab_out,msg,'COLL')
1912 0 : call wrtout(std_out,msg,'COLL')
1913 : end if ! end prtvol >= 10 output
1914 :
1915 166 : ptot_cart(:)=zero
1916 664 : do ii = 1,3
1917 : ptot_cart(ii)=rprimd(ii,1)*red_ptot(1) + rprimd(ii,2)*red_ptot(2) + &
1918 664 : & rprimd(ii,3)*red_ptot(3)
1919 : end do
1920 664 : ptot_cart(:)=ptot_cart(:)/ucvol
1921 :
1922 : ! !===================================================================================================
1923 : ! ! OUTPUT for fixed E
1924 : ! !===================================================================================================
1925 :
1926 166 : if (dtset%berryopt == 4) then
1927 :
1928 : ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar
1929 166 : write(msg,'(a,a)') ch10, 'scfcv: Constant unreduced E-field:'
1930 166 : call wrtout(std_out,msg,'COLL')
1931 166 : call prtefield(dtset,dtefield,std_out,rprimd)
1932 166 : if(dtset%prtvol>=10)then
1933 0 : call wrtout(ab_out,msg,'COLL')
1934 0 : call prtefield(dtset,dtefield,ab_out,rprimd)
1935 : end if
1936 : end if ! end berryopt 4 output
1937 :
1938 : ! =====================================================================================
1939 : ! ! fixed D calculation
1940 : ! !====================================================================================
1941 166 : if (dtset%berryopt == 6) then
1942 0 : if (scfcv_step > 1) then
1943 :
1944 : ! ! update efield taking damping into account dfield is in cartesian in dtset structure (contains input value)
1945 : ! ! same goes for efield - update the dtset%efield value
1946 : efield_test_cart(:)=dtset%ddamp*(dtset%dfield(:)-4.0d0*pi*ptot_cart(:))+&
1947 0 : & (1.0d0-dtset%ddamp)*efield_old_cart(:)
1948 :
1949 : ! ! test whether change in efield in any direction exceed maxestep, if so, set the
1950 : ! ! change to maxestep instead ! need optimized !
1951 0 : do ii = 1,3
1952 :
1953 0 : if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then
1954 :
1955 0 : write(std_out,'(a,a,i5)') "JH - "," E-field component:",ii
1956 0 : write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), &
1957 0 : & ", E(n-1)=",efield_old_cart(ii), ", E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), &
1958 0 : & ", maxestep=",dtset%maxestep
1959 :
1960 :
1961 0 : if (efield_test_cart(ii) > efield_old_cart(ii)) then
1962 0 : efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep)
1963 : else
1964 0 : efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep)
1965 : end if
1966 : end if
1967 : end do
1968 :
1969 0 : dtset%efield(:) = efield_test_cart(:)
1970 :
1971 : ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar
1972 0 : write(msg,'(a,a)') ch10, 'scfcv: Constant unreduced D-field - updating E-field:'
1973 0 : call wrtout(std_out,msg,'COLL')
1974 0 : call prtefield(dtset,dtefield,std_out,rprimd)
1975 0 : if(dtset%prtvol>=10)then
1976 0 : call wrtout(ab_out,msg,'COLL')
1977 0 : call prtefield(dtset,dtefield,ab_out,rprimd)
1978 : end if
1979 :
1980 : ! ! need to update dtset%efield_dot(:) with new value
1981 0 : dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1))
1982 0 : dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2))
1983 0 : dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3))
1984 :
1985 : else
1986 :
1987 0 : write(msg,'(a,a)') ch10, 'scfcv: Constant unreduced D-field - Pre E-field:'
1988 0 : call wrtout(std_out,msg,'COLL')
1989 0 : call prtefield(dtset,dtefield,std_out,rprimd)
1990 0 : if(dtset%prtvol>=10)then
1991 0 : call wrtout(ab_out,msg,'COLL')
1992 0 : call prtefield(dtset,dtefield,ab_out,rprimd)
1993 : end if
1994 :
1995 : end if ! scfcv_step >1
1996 :
1997 0 : efield_old_cart(:)=dtset%efield(:)
1998 : end if ! berryopt ==6
1999 : ! !===================================================================================================
2000 : ! ! fixed reduced d calculation
2001 : ! !===================================================================================================
2002 166 : if (dtset%berryopt == 16) then
2003 :
2004 0 : if (scfcv_step > 1) then
2005 : ! ! update efield taking damping into account reduced red_dfield
2006 : ! red_efield2 is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009) [[cite:Stengel2009]]
2007 :
2008 0 : red_efield2(:)=dtset%ddamp*(dtset%red_dfield(:)-red_ptot(:))+ (1.0d0-dtset%ddamp)*red_efield2_old(:)
2009 :
2010 : ! to calculate unreduced E
2011 0 : efield_test_cart(:)=(4*pi/ucvol)*(rprimd(:,1)*red_efield2(1)+rprimd(:,2)*red_efield2(2)+rprimd(:,3)*red_efield2(3))
2012 :
2013 : ! ! test whether change in efield in any direction exceed maxestep, if so, set the
2014 : ! ! change to maxestep instead ! need optimized !
2015 0 : do ii = 1,3
2016 :
2017 0 : if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then
2018 :
2019 0 : write(std_out,'(a,a,i5)') "JH - "," E-field component:",ii
2020 0 : write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), &
2021 0 : & ", E(n-1)=",efield_old_cart(ii), ", E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), &
2022 0 : & ", maxestep=",dtset%maxestep
2023 :
2024 0 : if (efield_test_cart(ii) > efield_old_cart(ii)) then
2025 0 : efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep)
2026 : else
2027 0 : efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep)
2028 : end if
2029 : end if
2030 : end do
2031 :
2032 0 : dtset%efield(:) = efield_test_cart(:)
2033 :
2034 : ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar
2035 0 : write(msg,'(a,a)') ch10, 'scfcv: Constant reduced d-field - updating E-field:'
2036 0 : call wrtout(std_out,msg,'COLL')
2037 0 : call prtefield(dtset,dtefield,std_out,rprimd)
2038 0 : if(dtset%prtvol>=10)then
2039 0 : call wrtout(ab_out,msg,'COLL')
2040 0 : call prtefield(dtset,dtefield,ab_out,rprimd)
2041 : end if
2042 :
2043 : ! ! need to update dtset%efield_dot(:) with new value
2044 : ! ! This needs to be deleted when efield_dot is deleted
2045 0 : dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1))
2046 0 : dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2))
2047 0 : dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3))
2048 :
2049 : else
2050 :
2051 0 : write(msg,'(a,a)') ch10, 'scfcv: Constant reduced d-field - Pre E-field:'
2052 0 : call wrtout(std_out,msg,'COLL')
2053 0 : call prtefield(dtset,dtefield,std_out,rprimd)
2054 0 : if(dtset%prtvol>=10)then
2055 0 : call wrtout(ab_out,msg,'COLL')
2056 0 : call prtefield(dtset,dtefield,ab_out,rprimd)
2057 : end if
2058 :
2059 : end if ! scfcv_step > 1
2060 :
2061 0 : efield_old_cart(:)=dtset%efield(:)
2062 0 : red_efield2_old(:)=red_efield2(:)
2063 : end if ! berryopt ==16
2064 :
2065 :
2066 : ! !===================================================================================================
2067 : ! ! fixed reduced d and ebar calculation (mixed BC)
2068 : ! !===================================================================================================
2069 166 : if (dtset%berryopt == 17) then
2070 :
2071 0 : if (scfcv_step > 1) then
2072 : ! ! update efield taking damping into account reduced red_dfield
2073 : ! red_efield1 and red_efield2 is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009) [[cite:Stengel1999]]
2074 : ! red_efield1 for fixed ebar, red_efield2 for fixed d calculation
2075 :
2076 : ! save this value in order to print the final value of real electric field, comparing with the desired red_fieldbar
2077 0 : dtefield%efield2(:)=dtset%efield(:)
2078 :
2079 : ! write(*,'(a,3i4)') "jfielddir=", (dtset%jfielddir(ii),ii=1,3)
2080 :
2081 0 : do ii=1,3
2082 0 : if (dtset%jfielddir(ii) ==2 ) then ! direction under fixed d
2083 0 : dtset%red_efieldbar(ii) = dot_product(dtset%efield(:),rprimd(:,ii)) ! update ebar which is not fixed
2084 0 : dtefield%efield_dot(ii) = dot_product(dtset%efield(:),rprimd(:,ii))
2085 : red_efield2(ii)=dtset%ddamp*(dtset%red_dfield(ii) - red_ptot(ii)) + &
2086 0 : & (1.0d0-dtset%ddamp)*red_efield2_old(ii) ! d(ii) is fixed, update e(ii) may need ddamping here
2087 :
2088 : ! write(msg,'(a,a,i5,a,i5)') ch10, 'direction ', ii,' for fixed d, value is (2) ', dtset%jfielddir(ii)
2089 : ! call wrtout(ab_out,msg,'COLL')
2090 : ! call wrtout(std_out,msg,'COLL')
2091 :
2092 0 : else if (dtset%jfielddir(ii) ==1 ) then ! direction under fixed ebar
2093 0 : red_efield2(ii)= (ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii)) ! update e which is not fixed
2094 0 : dtset%red_dfield(ii)=red_ptot(ii) + (ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii)) ! update d
2095 :
2096 : ! write(msg,'(a,a,i5,a,i5)') ch10, 'direction ', ii,' for fixed ebar, value is (1) ', dtset%jfielddir(ii)
2097 : ! call wrtout(ab_out,msg,'COLL')
2098 : ! call wrtout(std_out,msg,'COLL')
2099 :
2100 : end if
2101 : end do
2102 :
2103 0 : do ii=1,3
2104 0 : red_efield1(ii) =(ucvol/(4*pi))*dot_product(dtset%red_efieldbar(:),gmet(:,ii))
2105 : end do
2106 :
2107 :
2108 0 : dtset%red_efield(:)=(red_efield1(:) + red_efield2(:))/2.0d0 ! average reduced efield,
2109 : ! one is from fixed ebar part,
2110 : ! the other is from fixed d part.
2111 : ! This may need to be optimized !!
2112 :
2113 0 : write(msg,'(a,a,a,a,3(es16.9,2x),a)') ch10, 'Reduced efield from fixed ebar:', ch10, &
2114 0 : & ' e: ', (red_efield1(ii),ii=1,3), ch10
2115 :
2116 : ! call wrtout(ab_out,msg,'COLL')
2117 0 : call wrtout(std_out,msg,'COLL')
2118 :
2119 0 : write(msg,'(a,a,a,a,3(es16.9,2x),a)') ch10, 'Reduced efield from fixed d:', ch10, &
2120 0 : & ' e: ', (red_efield2(ii),ii=1,3), ch10
2121 :
2122 : ! call wrtout(ab_out,msg,'COLL')
2123 0 : call wrtout(std_out,msg,'COLL')
2124 :
2125 0 : write(msg,'(a,a,a,a,3(es16.9,2x),a)') ch10, 'Average reduced efield:', ch10, &
2126 0 : & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10
2127 :
2128 : ! call wrtout(ab_out,msg,'COLL')
2129 0 : call wrtout(std_out,msg,'COLL')
2130 :
2131 : ! to calculate unreduced E
2132 0 : do ii=1,3
2133 0 : efield_test_cart(ii) = (4*pi/ucvol)* dot_product(dtset%red_efield(:),rprimd(:,ii))
2134 : end do
2135 :
2136 : ! ! test whether change in efield in any direction exceed maxestep, if so, set the
2137 : ! ! change to maxestep instead ! need optimized !
2138 0 : do ii = 1,3
2139 0 : if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then
2140 :
2141 0 : write(std_out,'(a,a,i5)') "JH - "," E-field component:",ii
2142 0 : write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), &
2143 0 : & ", E(n-1)=",efield_old_cart(ii), ", E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), &
2144 0 : & ", maxestep=",dtset%maxestep
2145 :
2146 0 : if (efield_test_cart(ii) > efield_old_cart(ii)) then
2147 0 : efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep)
2148 : else
2149 0 : efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep)
2150 : end if
2151 : end if
2152 : end do
2153 :
2154 0 : dtset%efield(:) = efield_test_cart(:)
2155 :
2156 : ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar
2157 0 : write(msg,'(a,a)') ch10, 'scfcv: Constant reduced ebar and d-field - updating E-field:'
2158 0 : call wrtout(std_out,msg,'COLL')
2159 0 : call prtefield(dtset,dtefield,std_out,rprimd)
2160 0 : if(dtset%prtvol>=10)then
2161 0 : call wrtout(ab_out,msg,'COLL')
2162 0 : call prtefield(dtset,dtefield,ab_out,rprimd)
2163 : end if
2164 :
2165 :
2166 : ! ! need to update dtset%efield_dot(:) with new value
2167 : ! ! This needs to be deleted when efield_dot is deleted
2168 0 : dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1))
2169 0 : dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2))
2170 0 : dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3))
2171 :
2172 : else
2173 :
2174 0 : write(msg,'(a,a)') ch10, 'scfcv: Constant reduced ebar and d-field - Pre E-field:'
2175 0 : call wrtout(std_out,msg,'COLL')
2176 0 : call prtefield(dtset,dtefield,std_out,rprimd)
2177 0 : if(dtset%prtvol>=10)then
2178 0 : call wrtout(ab_out,msg,'COLL')
2179 0 : call prtefield(dtset,dtefield,ab_out,rprimd)
2180 : end if
2181 :
2182 : end if ! scfcv_step > 1
2183 :
2184 0 : efield_old_cart(:)=dtset%efield(:)
2185 0 : red_efield2_old(:)=red_efield2(:)
2186 :
2187 : end if ! berryopt ==17
2188 :
2189 : end if ! end efield .and. scfcv_level 1 tasks
2190 :
2191 : !deallocate cprj
2192 191 : if ( efield .and. psps%usepaw == 1) then
2193 69 : call pawcprj_free(cprj)
2194 : end if
2195 87395 : ABI_FREE(cprj)
2196 :
2197 52323 : end subroutine update_e_field_vars
2198 : !!***
2199 :
2200 : !!****f* ABINIT/prtefield
2201 : !!
2202 : !! NAME
2203 : !! prtefield
2204 : !!
2205 : !! FUNCTION
2206 : !! Print components of electric field, displacement field and polarization in nice format
2207 : !!
2208 : !! INPUTS
2209 : !! dtset <type(dataset_type)>=all input variables in this dataset
2210 : !! | berryopt
2211 : !! | efield
2212 : !! | dfield
2213 : !! | red_efield
2214 : !! | red_efieldbar
2215 : !! | red_dfield
2216 : !! dtefield <type(efield_type)>
2217 : !! | efield2
2218 : !! | red_ptot1
2219 : !! iunit = unit number to which the data is printed
2220 : !! rprimd
2221 : !!
2222 : !! OUTPUT
2223 : !! (only writing)
2224 : !!
2225 : !! SOURCE
2226 :
2227 216 : subroutine prtefield(dtset,dtefield,iunit,rprimd)
2228 :
2229 : !Arguments ------------------------------------
2230 : integer :: iunit
2231 : real(dp),intent(in) :: rprimd(3,3)
2232 : type(efield_type),intent(in) :: dtefield
2233 : type(dataset_type),intent(inout) :: dtset
2234 :
2235 : !Local variables-------------------------------
2236 : ! Do not modify the length of this string
2237 : !scalars
2238 : integer :: idir,ii
2239 : character(len=1500) :: msg
2240 : character(len=7) :: flag_field(3)
2241 :
2242 : real(dp) :: ucvol
2243 : ! arrays
2244 : real(dp) :: ptot_cart(3),gmet(3,3),gprimd(3,3),rmet(3,3),red_pbar(3),red_dbar(3),red_dfieldbar(3)
2245 : real(dp) :: red_efieldbar_lc(3),red_efield_lc(3)
2246 :
2247 :
2248 : ! *************************************************************************
2249 :
2250 : !DEBUG
2251 : !write(iout, '(a)') ' prtefield : enter '
2252 : !ENDDEBUG
2253 :
2254 : !write here
2255 :
2256 216 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2257 :
2258 :
2259 216 : ptot_cart(:)=zero
2260 864 : do idir = 1,3
2261 : ptot_cart(idir)=rprimd(idir,1)*dtefield%red_ptot1(1) + rprimd(idir,2)*dtefield%red_ptot1(2) + &
2262 864 : & rprimd(idir,3)*dtefield%red_ptot1(3)
2263 : end do
2264 864 : ptot_cart(:)=ptot_cart(:)/ucvol
2265 :
2266 216 : if (dtset%berryopt == 4) then
2267 :
2268 : ! to calculate e Eq.(25)
2269 864 : do idir=1,3
2270 2808 : dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,idir))
2271 : end do
2272 :
2273 : ! to calculate ebar Eq.(25)
2274 864 : do idir=1,3
2275 2808 : dtset%red_efieldbar(idir) =dot_product(dtset%efield(:),rprimd(:,idir))
2276 : end do
2277 :
2278 :
2279 : ! to calculate pbar
2280 864 : do idir=1,3
2281 2808 : red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2282 : end do
2283 :
2284 : !MGNAG: This msg is too long and causes
2285 : ! Runtime Error: wrtout_cpp.f90, line 893: Buffer overflow on output
2286 : ! with NAG in test seq_tsv6_125 where we write to std_out!
2287 : ! I cannot change the RECLEN of std_out!
2288 :
2289 216 : write(msg,'(a,a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ' (a. u.)', ch10,&
2290 1080 : & ' E: ', (dtset%efield(ii), ii=1,3), ch10, &
2291 1296 : & ' P: ', (ptot_cart(ii), ii=1,3)
2292 216 : call wrtout(iunit,msg,'COLL')
2293 :
2294 216 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10,&
2295 1080 : & ' ebar: ', (dtset%red_efieldbar(ii),ii=1,3), ch10, & !!HONG need to change
2296 1296 : & ' pbar: ', (red_pbar(ii),ii=1,3)
2297 216 : call wrtout(iunit,msg,'COLL')
2298 :
2299 216 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10,&
2300 1080 : & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10, &
2301 1944 : & ' p: ', (dtefield%red_ptot1(ii), ii=1,3)
2302 216 : call wrtout(iunit,msg,'COLL')
2303 :
2304 216 : write(msg,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10,&
2305 216 : & ' (S.I.), that is V/m for E, and C/m^2 for P', ch10, &
2306 1080 : & '- E: ', (dtset%efield(ii)*(Ha_J/(e_Cb*Bohr_Ang*1d-10)), ii=1,3), ch10, & !(Ha_J/(e_Cb*Bohr_Ang*1d-10))= 5.14220652*1d+11
2307 1944 : & ' P: ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3),ch10
2308 216 : call wrtout(iunit,msg,'COLL')
2309 :
2310 : end if ! berryopt ==4
2311 :
2312 :
2313 216 : if (dtset%berryopt == 6) then
2314 0 : do idir=1,3
2315 0 : dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,idir))
2316 : end do
2317 :
2318 : ! to calculate ebar !! Need to be changed
2319 0 : do idir=1,3
2320 0 : dtset%red_efieldbar(idir) = dot_product(dtset%efield(:),rprimd(:,idir))
2321 : end do
2322 :
2323 : ! to calculate red_pbar
2324 0 : do idir=1,3
2325 0 : red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2326 : end do
2327 :
2328 : ! to calculate red_dbar
2329 0 : do idir=1,3
2330 0 : red_dbar(idir) = dot_product(dtset%dfield(:),rprimd(:,idir))
2331 : end do
2332 :
2333 : ! to calculate d
2334 0 : do idir=1,3
2335 0 : dtset%red_dfield(idir) =(ucvol/(4*pi))*dot_product(dtset%dfield(:),gprimd(:,idir))
2336 : end do
2337 :
2338 0 : write(msg,'(a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ' (a. u.)', ch10,&
2339 0 : & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10, &
2340 0 : & ' p: ', (dtefield%red_ptot1(ii), ii=1,3), ch10, &
2341 0 : & ' d: ', (dtset%red_dfield(ii),ii = 1, 3), ch10, &
2342 0 : & ' e + p: ', (dtset%red_efield(ii)+dtefield%red_ptot1(ii),ii=1,3)
2343 0 : call wrtout(iunit,msg,'COLL')
2344 :
2345 0 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10,&
2346 0 : & ' ebar: ', (dtset%red_efieldbar(ii),ii=1,3), ch10, & !!HONG need to change
2347 0 : & ' pbar: ', (red_pbar(ii),ii=1,3), ch10, &
2348 0 : & ' dbar: ', (red_dbar(ii),ii=1,3), ch10, &
2349 0 : & ' eba+pba: ', (dtset%red_efieldbar(ii)+red_pbar(ii),ii=1,3)
2350 0 : call wrtout(iunit,msg,'COLL')
2351 :
2352 0 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2353 0 : & ' E: ', (dtset%efield(ii), ii=1,3), ch10, &
2354 0 : & ' P: ', (ptot_cart(ii), ii=1,3), ch10, &
2355 0 : & ' D: ', (dtset%dfield(ii),ii = 1, 3), ch10, &
2356 0 : & 'E+4*pi*P: ', (dtset%efield(ii)+4.0d0*pi*ptot_cart(ii),ii=1,3)
2357 0 : call wrtout(iunit,msg,'COLL')
2358 :
2359 0 : write(msg,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10,&
2360 0 : & ' (S.I.), that is V/m for E, and C/m^2 for P and D', ch10, &
2361 0 : & '- E: ', (dtset%efield(ii)*(Ha_J/(e_Cb*Bohr_Ang*1d-10)), ii=1,3), ch10,& !(Ha_J/(e_Cb*Bohr_Ang*1d-10))= 5.14220652*1d+11
2362 0 : & ' P: ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3), ch10,&
2363 0 : & ' D: ', ((1.0d0/(4*pi))*dtset%dfield(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii = 1, 3),ch10,&
2364 0 : & 'eps0*E+P: ', (dtset%efield(ii)*eps0*(Ha_J/(e_Cb*Bohr_Ang*1d-10))+ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii=1,3),ch10
2365 : ! eps0*(Ha_J/(e_Cb*Bohr_Ang*1d-10))=8.854187817620*5.14220652*1d-1
2366 0 : call wrtout(iunit,msg,'COLL')
2367 :
2368 : !MGNAG Runtime Error: wrtout_cpp.f90, line 896: Buffer overflow on output
2369 :
2370 : end if ! berryopt ==6
2371 :
2372 :
2373 216 : if (dtset%berryopt == 14) then
2374 :
2375 0 : do idir=1,3 ! ebar local
2376 0 : red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir))
2377 : end do
2378 :
2379 :
2380 0 : do idir=1,3
2381 0 : dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir))
2382 : end do
2383 :
2384 : ! to calculate pbar
2385 0 : do idir=1,3
2386 0 : red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2387 : end do
2388 :
2389 0 : write(msg,'(a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ' (a. u.)', ch10,&
2390 0 : & ' ebar0: ', (dtset%red_efieldbar(ii),ii=1,3), ch10, &
2391 0 : & ' ebar: ', (red_efieldbar_lc(ii),ii=1,3), ch10, &
2392 0 : & ' pbar: ', (red_pbar(ii),ii=1,3)
2393 0 : call wrtout(iunit,msg,'COLL')
2394 :
2395 0 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2396 0 : & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10, &
2397 0 : & ' p: ', (dtefield%red_ptot1(ii), ii=1,3)
2398 0 : call wrtout(iunit,msg,'COLL')
2399 :
2400 0 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2401 0 : & ' E: ', (dtefield%efield2(ii), ii=1,3), ch10, &
2402 0 : & ' P: ', (ptot_cart(ii), ii=1,3)
2403 0 : call wrtout(iunit,msg,'COLL')
2404 :
2405 0 : write(msg,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10, &
2406 0 : & ' (S.I.), that is V/m for E, and C/m^2 for P', ch10, &
2407 0 : & '- E: ', (dtefield%efield2(ii)*(Ha_J/(e_Cb*Bohr_Ang*1d-10)), ii=1,3), ch10, & !(Ha_J/(e_Cb*Bohr_Ang*1d-10))= 5.14220652*1d+11
2408 0 : & ' P: ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3),ch10
2409 0 : call wrtout(iunit,msg,'COLL')
2410 :
2411 :
2412 : end if ! berryopt ==14
2413 :
2414 :
2415 216 : if (dtset%berryopt == 16) then
2416 :
2417 : ! to calculate e Eq.(25)
2418 0 : do idir=1,3
2419 0 : dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,idir))
2420 : end do
2421 :
2422 : ! to calculate ebar
2423 0 : do idir=1,3
2424 0 : dtset%red_efieldbar(idir) = dot_product(dtset%efield(:),rprimd(:,idir))
2425 : end do
2426 :
2427 : ! to calculate pbar
2428 0 : do idir=1,3
2429 0 : red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2430 : end do
2431 :
2432 : ! to calculate dbar
2433 0 : do idir=1,3
2434 0 : red_dfieldbar(idir) = (4*pi/ucvol)*dot_product(dtset%red_dfield(:),rmet(:,idir))
2435 : end do
2436 :
2437 : ! to calculate D
2438 0 : do idir=1,3
2439 0 : dtset%dfield(idir) =(4*pi/ucvol)*dot_product(dtset%red_dfield(:),rprimd(:,idir))
2440 : end do
2441 :
2442 0 : write(msg,'(a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ' (a. u.)', ch10,&
2443 0 : & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10, &
2444 0 : & ' p: ', (dtefield%red_ptot1(ii), ii=1,3), ch10, &
2445 0 : & ' d: ', (dtset%red_dfield(ii),ii = 1, 3), ch10, &
2446 0 : & ' e + p: ', (dtset%red_efield(ii)+dtefield%red_ptot1(ii),ii=1,3)
2447 0 : call wrtout(iunit,msg,'COLL')
2448 :
2449 0 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2450 0 : & ' ebar: ', (dtset%red_efieldbar(ii),ii=1,3), ch10, &
2451 0 : & ' pbar: ', (red_pbar(ii),ii=1,3), ch10, &
2452 0 : & ' dbar: ', (red_dfieldbar(ii),ii=1,3), ch10, &
2453 0 : & ' eba+pba: ', (dtset%red_efieldbar(ii)+red_pbar(ii),ii=1,3)
2454 0 : call wrtout(iunit,msg,'COLL')
2455 :
2456 0 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2457 0 : & ' E: ', (dtset%efield(ii), ii=1,3), ch10, &
2458 0 : & ' P: ', (ptot_cart(ii), ii=1,3), ch10, &
2459 0 : & ' D: ', (dtset%dfield(ii),ii = 1, 3), ch10, &
2460 0 : & 'E+4*pi*P: ', (dtset%efield(ii)+4.0d0*pi*ptot_cart(ii),ii=1,3)
2461 0 : call wrtout(iunit,msg,'COLL')
2462 :
2463 0 : write(msg,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10, &
2464 0 : & ' (S.I.), that is V/m for E, and C/m^2 for P and D', ch10, &
2465 0 : & '- E: ', (dtset%efield(ii)*(Ha_J/(e_Cb*Bohr_Ang*1d-10)), ii=1,3), ch10, & !(Ha_J/(e_Cb*Bohr_Ang*1d-10))= 5.14220652*1d+11
2466 0 : & ' P: ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3), ch10, &
2467 0 : & ' D: ', ((1.0d0/(4*pi))*dtset%dfield(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii = 1, 3), ch10, &
2468 0 : & 'eps0*E+P: ', (dtset%efield(ii)*eps0*(Ha_J/(e_Cb*Bohr_Ang*1d-10))+ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii=1,3),ch10
2469 0 : call wrtout(iunit,msg,'COLL')
2470 :
2471 : end if ! berryopt ==16
2472 :
2473 216 : if (dtset%berryopt == 17) then
2474 :
2475 0 : do idir=1,3 ! ebar local
2476 0 : red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir))
2477 : end do
2478 :
2479 0 : do idir=1,3
2480 0 : dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir))
2481 : end do
2482 :
2483 : ! to calculate pbar
2484 0 : do idir=1,3
2485 0 : red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2486 : end do
2487 :
2488 :
2489 : ! do idir=1,3
2490 : ! if (dtset%rfdir(idir)==1) then ! fixed ebar
2491 : ! red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir)) ! local efieldbar
2492 : ! dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir))
2493 : ! red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2494 : ! dtset%dfield(idir)=dtefield%efield2(idir)+4*pi*ptot_cart(idir)
2495 : ! dtset%red_dfield(idir)=dtset%red_efield+dtefield%red_ptot1(idir)
2496 : ! dtset%red_dfieldbar(idir)=red_efieldbar_lc(idir)+red_pbar(idir)
2497 : ! E_lc(idir)=dtefield%efield2(idir)
2498 : ! e_lc(idir)=red_efieldbar_lc(idir)
2499 : ! ebar_lc(idir)=dtset%red_efieldbar(idir)
2500 : ! else if (dtset%rfdir(idir)==2) then ! fixed d
2501 : ! dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir))
2502 : ! dtset%red_efieldbar(idir) = dot_product(dtefield%efield2(:),rprimd(:,idir))
2503 : ! red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2504 : ! red_dfieldbar(idir) = (4*pi/ucvol)*dot_product(dtset%red_dfield(:),rmet(:,idir))
2505 : ! dtset%dfield(idir) =(4*pi/ucvol)*dot_product(dtset%red_dfield(:),rprimd(:,idir))
2506 : ! E_lc(idir)=dtefield%efield2(idir)
2507 : ! e_lc(idir)=dtset%red_efield(idir)
2508 : ! ebar_lc(idir)=dtset%red_efieldbar(idir)
2509 : ! end if
2510 : ! enddo
2511 :
2512 :
2513 :
2514 0 : do idir=1,3
2515 0 : red_efield_lc(idir)= (ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir))
2516 0 : red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir)) ! local efieldbar
2517 0 : red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2518 0 : red_dfieldbar(idir) = (4*pi/ucvol)*dot_product(dtset%red_dfield(:),rmet(:,idir))
2519 0 : dtset%dfield(idir) =(4*pi/ucvol)*dot_product(dtset%red_dfield(:),rprimd(:,idir))
2520 :
2521 : end do
2522 :
2523 0 : do idir=1,3
2524 0 : if(dtset%jfielddir(idir)==1) then
2525 0 : flag_field(idir)="E-field"
2526 : else
2527 0 : flag_field(idir)="D-field"
2528 : end if
2529 : end do
2530 :
2531 : write(msg,'(a,a,a,6x,a,11x,a,11x,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') &
2532 0 : & ' (a. u.)', ch10,&
2533 0 : & ' ', (flag_field(ii),ii=1,3),ch10, &
2534 0 : & ' ebar0: ', (dtset%red_efieldbar(ii),ii=1,3), ch10, &
2535 0 : & ' ebar: ', (red_efieldbar_lc(ii),ii=1,3), ch10, &
2536 0 : & ' d: ', (dtset%red_dfield(ii),ii = 1, 3), ch10, &
2537 0 : & ' e + p: ', (red_efield_lc(ii)+dtefield%red_ptot1(ii),ii=1,3)
2538 0 : call wrtout(iunit,msg,'COLL')
2539 :
2540 0 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2541 0 : & ' e: ', (red_efield_lc(ii),ii=1,3), ch10, &
2542 0 : & ' p: ', (dtefield%red_ptot1(ii), ii=1,3), ch10, &
2543 0 : & ' d: ', (dtset%red_dfield(ii),ii = 1, 3), ch10, &
2544 0 : & ' e + p: ', (red_efield_lc(ii)+dtefield%red_ptot1(ii),ii=1,3)
2545 0 : call wrtout(iunit,msg,'COLL')
2546 :
2547 0 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2548 0 : & ' ebar: ', (red_efieldbar_lc(ii),ii=1,3), ch10, &
2549 0 : & ' pbar: ', (red_pbar(ii),ii=1,3), ch10, &
2550 0 : & ' dbar: ', (red_dfieldbar(ii),ii=1,3), ch10, &
2551 0 : & ' eba+pba: ', (red_efieldbar_lc(ii)+red_pbar(ii),ii=1,3)
2552 0 : call wrtout(iunit,msg,'COLL')
2553 :
2554 0 : write(msg,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2555 0 : & ' E: ', (dtefield%efield2(ii), ii=1,3), ch10, &
2556 0 : & ' P: ', (ptot_cart(ii), ii=1,3), ch10, &
2557 0 : & ' D: ', (dtset%dfield(ii),ii = 1, 3), ch10, &
2558 0 : & 'E+4*pi*P: ', (dtset%efield(ii)+4.0d0*pi*ptot_cart(ii),ii=1,3)
2559 0 : call wrtout(iunit,msg,'COLL')
2560 :
2561 0 : write(msg,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10, &
2562 0 : & ' (S.I.), that is V/m for E, and C/m^2 for P and D', ch10, &
2563 0 : & ' E: ', (dtefield%efield2(ii)*(Ha_J/(e_Cb*Bohr_Ang*1d-10)), ii=1,3), ch10,& !(Ha_J/(e_Cb*Bohr_Ang*1d-10))= 5.14220652*1d+11
2564 0 : & ' P: ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3), ch10, &
2565 0 : & ' D: ', ((1.0d0/(4*pi))*dtset%dfield(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii = 1, 3), ch10, &
2566 0 : & 'eps0*E+P: ', (dtefield%efield2(ii)*eps0*(Ha_J/(e_Cb*Bohr_Ang*1d-10))+ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii=1,3),ch10
2567 0 : call wrtout(iunit,msg,'COLL')
2568 :
2569 : end if ! berryopt ==17
2570 :
2571 216 : end subroutine prtefield
2572 : !!***
2573 :
2574 : !!****f* ABINIT/init_e_field_vars
2575 : !! NAME
2576 : !! init_e_field_vars
2577 : !!
2578 : !! FUNCTION
2579 : !! Initialization of variables and data structures used in polarization
2580 : !! calculations
2581 : !!
2582 : !! INPUTS
2583 : !! dtset <type(dataset_type)> = all input variables in this dataset
2584 : !! gmet(3,3) = reciprocal space metric tensor in bohr**-2
2585 : !! gprimd(3,3) = primitive translations in recip space
2586 : !! kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere
2587 : !! mpi_enreg=information about MPI parallelization
2588 : !! npwarr(nkpt) = number of planewaves in basis and boundary at this k point
2589 : !! occ(mband*nkpt*nsppol) = occup number for each band at each k point
2590 : !! pawang <type(pawang_type)>=paw angular mesh and related data
2591 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
2592 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
2593 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
2594 : !! rprimd(3,3) = dimensional primitive vectors
2595 : !! symrec(3,3,nsym) = symmetries in reciprocal space in terms of
2596 : !! reciprocal space primitive translations
2597 : !! xred(3,natom) = location of atoms in reduced units
2598 : !!
2599 : !! OUTPUT
2600 : !! dtefield <type(efield_type)> :: initialized polarization variables
2601 : !! pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat
2602 : !! between k-points k and k +- dk where dk is
2603 : !! parallel to the direction idir
2604 : !! pwind_alloc = first dimension of pwind and pwnsfac
2605 : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
2606 : !!
2607 : !! SIDE EFFECTS
2608 : !!
2609 : !! TO DO
2610 : !!
2611 : !! NOTES
2612 : !!
2613 : !! SOURCE
2614 :
2615 5629 : subroutine init_e_field_vars(dtefield,dtset,gmet,gprimd,kg,&
2616 5629 : & mpi_enreg,npwarr,occ,pawang,pawrad,pawtab,psps,&
2617 5629 : & pwind,pwind_alloc,pwnsfac,rprimd,symrec,xred)
2618 :
2619 : !Arguments ------------------------------------
2620 : !scalars
2621 : integer,intent(out) :: pwind_alloc
2622 : type(MPI_type),intent(inout) :: mpi_enreg
2623 : type(dataset_type),intent(inout) :: dtset
2624 : type(efield_type),intent(inout) :: dtefield !vz_i needs efield2
2625 : type(pawang_type),intent(in) :: pawang
2626 : type(pseudopotential_type),intent(in) :: psps
2627 : !arrays
2628 : integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
2629 : integer,intent(in) :: symrec(3,3,dtset%nsym)
2630 : integer,pointer :: pwind(:,:,:)
2631 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
2632 : real(dp),intent(in) :: rprimd(3,3),xred(3,dtset%natom)
2633 : real(dp),pointer :: pwnsfac(:,:)
2634 : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw)
2635 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
2636 :
2637 : !Local variables-------------------------------
2638 : logical :: initfield
2639 : !scalars
2640 :
2641 : ! *************************************************************************
2642 :
2643 5629 : initfield = .false.
2644 :
2645 : !initialization
2646 5629 : dtefield%has_qijb = 0
2647 5629 : dtefield%has_epawf3 = 0
2648 5629 : dtefield%has_epaws3 = 0
2649 5629 : dtefield%has_expibi = 0
2650 5629 : dtefield%has_rij = 0
2651 5629 : dtefield%usecprj = 0
2652 5629 : dtefield%berryopt = 0
2653 :
2654 : if ((dtset%berryopt < 0).or.(dtset%berryopt == 4) .or. (dtset%berryopt == 6) .or.(dtset%berryopt == 7) .or. &
2655 5629 : & (dtset%berryopt == 14) .or.(dtset%berryopt == 16) .or.(dtset%berryopt == 17)) then
2656 82 : nullify(pwind,pwnsfac)
2657 : call initberry(dtefield,dtset,gmet,gprimd,kg,&
2658 : & dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw,&
2659 : & dtset%natom,dtset%nkpt,npwarr,dtset%nsppol,&
2660 : & dtset%nsym,dtset%ntypat,occ,pawang,pawrad,pawtab,&
2661 : & psps,pwind,pwind_alloc,pwnsfac,rprimd,symrec,&
2662 82 : & dtset%typat,psps%usepaw,xred)
2663 : initfield = .true.
2664 : end if
2665 :
2666 5547 : if (.not. initfield .and. dtset%orbmag == 0) then
2667 : ! initorbmag.F90 also allocates pwind and pwnsfac
2668 5547 : pwind_alloc = 1
2669 5547 : ABI_MALLOC(pwind,(pwind_alloc,2,3))
2670 5547 : ABI_MALLOC(pwnsfac,(2,pwind_alloc))
2671 88752 : pwind(:,:,:)=0
2672 27735 : pwnsfac(:,:)=zero
2673 : end if
2674 :
2675 5629 : end subroutine init_e_field_vars
2676 : !!***
2677 :
2678 : !!****f* ABINIT/initberry
2679 : !! NAME
2680 : !! initberry
2681 : !!
2682 : !! FUNCTION
2683 : !! Initialization of Berryphase calculation of the polarization, the
2684 : !! ddk and the response of an insulator to a homogenous electric field.
2685 : !!
2686 : !! INPUTS
2687 : !! dtset <type(dataset_type)> = all input variables in this dataset
2688 : !! gmet(3,3) = reciprocal space metric tensor in bohr**-2
2689 : !! gprimd(3,3) = primitive translations in recip space
2690 : !! kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere
2691 : !! mband = maximum number of bands
2692 : !! mkmem = maximum number of k-points in core memory
2693 : !! mpw = maximum number of plane waves
2694 : !! natom = number of atoms in unit cell
2695 : !! nkpt = number of k points
2696 : !! npwarr(nkpt) = number of planewaves in basis and boundary at this k point
2697 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
2698 : !! nsym = number of symmetry operations
2699 : !! ntypat = number of types of atoms in unit cell
2700 : !! occ(mband*nkpt*nsppol) = occup number for each band at each k point
2701 : !! pawang <type(pawang_type)>=paw angular mesh and related data
2702 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
2703 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
2704 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
2705 : !! rprimd(3,3) = dimensional primitive vectors
2706 : !! symrec(3,3,nsym) = symmetries in reciprocal space in terms of
2707 : !! reciprocal space primitive translations
2708 : !! typat = typat(natom) list of atom types
2709 : !! usepaw = flag for PAW (1 PAW, 0 NCPP)
2710 : !! xred(3,natom) = location of atoms in reduced units
2711 : !!
2712 : !! OUTPUT
2713 : !! dtefield <type(efield_type)> = variables related to Berry phase
2714 : !! calculations
2715 : !! pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat
2716 : !! between k-points k and k +- dk where dk is
2717 : !! parallel to the direction idir
2718 : !! jpw = pwind(ipw,ifor,idir)
2719 : !! * ipw = index of plane wave vector G for a given k-point k
2720 : !! * ifor = 1: k + dk
2721 : !! 2: k - dk
2722 : !! * idir = direction of the polarization/ddk calculation [dk(idir)
2723 : !! is the only non-zero element of dk(:)]
2724 : !! * jpw = index of plane wave vector G (+dG) at k +- dk
2725 : !! where dG is a shift of one reciprocal lattice vector
2726 : !! (required to close the strings of k-points using the
2727 : !! periodic gauge condition)
2728 : !! In case a G-vector of the basis sphere of plane waves at k
2729 : !! does not belong to the basis sphere of plane waves at k+dk, jpw = 0.
2730 : !! pwind_alloc = first dimension of pwind and pwnsfac
2731 : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
2732 : !!
2733 : !! SIDE EFFECTS
2734 : !! mpi_enreg = information about MPI parallelization
2735 : !! kptdstrb(nproc,nneighbour,fmkmem_max*nsppol) : Array required
2736 : !! by berryphase_new.f for MPI // over k-points. Defined
2737 : !! for k-points in the fBZ
2738 : !! but for k-points in the iBZ. Used by vtorho.f
2739 : !! nproc = number of cpus
2740 : !! nneighbour = number of neighbours for each k-point (= 6)
2741 : !!
2742 : !! SOURCE
2743 :
2744 82 : subroutine initberry(dtefield,dtset,gmet,gprimd,kg,mband,&
2745 82 : mkmem,mpi_enreg,mpw,natom,nkpt,npwarr,nsppol,&
2746 82 : nsym,ntypat,occ,pawang,pawrad,pawtab,psps,&
2747 : pwind,pwind_alloc,pwnsfac,&
2748 82 : rprimd,symrec,typat,usepaw,xred)
2749 :
2750 : !Arguments ------------------------------------
2751 : !scalars
2752 : integer,intent(in) :: mband,mkmem,mpw,natom,nkpt,nsppol,nsym,ntypat,usepaw
2753 : integer,intent(out) :: pwind_alloc
2754 : type(MPI_type),intent(inout) :: mpi_enreg
2755 : type(dataset_type),intent(inout) :: dtset
2756 : type(efield_type),intent(inout) :: dtefield !vz_i
2757 : type(pawang_type),intent(in) :: pawang
2758 : type(pseudopotential_type),intent(in) :: psps
2759 : !arrays
2760 : integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt)
2761 : integer,intent(in) :: symrec(3,3,nsym),typat(natom)
2762 : integer,pointer :: pwind(:,:,:)
2763 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(mband*nkpt*nsppol)
2764 : real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
2765 : real(dp),pointer :: pwnsfac(:,:)
2766 : type(pawrad_type),intent(in) :: pawrad(ntypat)
2767 : type(pawtab_type),intent(in) :: pawtab(ntypat)
2768 :
2769 : !Local variables-------------------------------
2770 : !scalars
2771 : integer :: exchn2n3d,flag,flag_kpt,fnkpt_computed,iband,icg,icprj
2772 : integer :: idir,idum,idum1,ierr,ifor,ikg,ikg1,ikpt,ikpt1,ikpt1f
2773 : integer :: ikpt1i,ikpt2,ikpt_loc,ikptf,ikpti,ikstr,index,ineigh,ipw,ipwnsfac
2774 : integer :: isppol,istr,istwf_k,isym,isym1,itrs,itypat,iunmark,jpw,klmn,lmax,lmn2_size_max
2775 : integer :: me,me_g0,mkmem_,my_nspinor,nband_k,mband_occ_k,ncpgr,nkstr,nproc,npw_k,npw_k1,spaceComm
2776 : integer :: option, brav, mkpt, nkptlatt
2777 : integer :: jstr,ii,jj,isign
2778 : integer :: dk_flag, coord1, coord2
2779 : integer :: mult
2780 : real(dp) :: c1,ecut_eff,eg,eg_ev,rdum,diffk1,diffk2,diffk3
2781 : real(dp) :: dist_, max_dist, last_dist, dist,kpt_shifted1,kpt_shifted2,kpt_shifted3
2782 : real(dp) :: gprimdlc(3,3),rmetllc(3,3),gmetllc(3,3),ucvol_local
2783 : ! gprimd(3,3) = inverse of rprimd
2784 : ! rmetlcl(3,3)=real-space metric (same as rmet in metric.F90)
2785 : ! gmetlcl(3,3)= same as gmet in metric.F90
2786 : ! ucvol = volume of the unit cell in Bohr**3
2787 : character(len=500) :: msg
2788 : logical :: calc_epaw3_force,calc_epaw3_stress,fieldflag
2789 : !arrays
2790 : integer :: dg(3),iadum(3),iadum1(3),neigh(6)
2791 82 : integer,allocatable :: dimlmn(:),kg1_k(:,:),kpt_mark(:),nattyp_dum(:)
2792 : real(dp) :: diffk(3),dk(3),dum33(3,3),eg_dir(3)
2793 : real(dp) :: kpt1(3)
2794 : real(dp) :: delta_str3(2), dstr(2),dk_str(2,2,3)
2795 : real(dp) :: tsec(2)
2796 82 : real(dp),allocatable :: calc_expibi(:,:),calc_qijb(:,:,:),spkpt(:,:)
2797 : ! *************************************************************************
2798 :
2799 : DBG_ENTER("COLL")
2800 :
2801 82 : call timab(1001,1,tsec)
2802 82 : call timab(1002,1,tsec)
2803 :
2804 82 : spaceComm=mpi_enreg%comm_cell
2805 82 : nproc=xmpi_comm_size(spaceComm)
2806 82 : me=xmpi_comm_rank(spaceComm)
2807 :
2808 : !save the current value of berryopt
2809 82 : dtefield%berryopt = dtset%berryopt
2810 : !save the current value of nspinor
2811 82 : dtefield%nspinor = dtset%nspinor
2812 :
2813 : !----------------------------------------------------------------------------
2814 : !-------------------- Obtain k-point grid in the full BZ --------------------
2815 : !----------------------------------------------------------------------------
2816 :
2817 82 : if(dtset%kptopt==1 .or. dtset%kptopt==2 .or. dtset%kptopt==4)then
2818 : ! Compute the number of k points in the G-space unit cell
2819 : nkptlatt=dtset%kptrlatt(1,1)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,3) &
2820 : & +dtset%kptrlatt(1,2)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,1) &
2821 : & +dtset%kptrlatt(1,3)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,2) &
2822 : & -dtset%kptrlatt(1,2)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,3) &
2823 : & -dtset%kptrlatt(1,3)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,1) &
2824 74 : & -dtset%kptrlatt(1,1)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,2)
2825 :
2826 : ! Call smpbz to obtain the list of k-point in the full BZ - without symmetry reduction
2827 74 : option = 0
2828 74 : brav = 1
2829 74 : mkpt=nkptlatt*dtset%nshiftk
2830 222 : ABI_MALLOC(spkpt,(3,mkpt))
2831 74 : call smpbz(1,ab_out,dtset%kptrlatt,mkpt,fnkpt_computed,dtset%nshiftk,option,dtset%shiftk,spkpt)
2832 74 : dtefield%fnkpt = fnkpt_computed
2833 222 : ABI_MALLOC(dtefield%fkptns,(3,dtefield%fnkpt))
2834 167966 : dtefield%fkptns(:,:)=spkpt(:,1:dtefield%fnkpt)
2835 74 : ABI_FREE(spkpt)
2836 8 : else if(dtset%kptopt==3.or.dtset%kptopt==0)then
2837 8 : dtefield%fnkpt=nkpt
2838 24 : ABI_MALLOC(dtefield%fkptns,(3,dtefield%fnkpt))
2839 1544 : dtefield%fkptns(1:3,1:dtefield%fnkpt)=dtset%kpt(1:3,1:dtefield%fnkpt)
2840 8 : if(dtset%kptopt==0)then
2841 0 : write(msg,'(10a)') ch10,&
2842 0 : & ' initberry : WARNING -',ch10,&
2843 0 : & ' you have defined manually the k-point grid with kptopt = 0',ch10,&
2844 0 : & ' the berry phase calculation works only with a regular k-points grid,',ch10,&
2845 0 : & ' abinit doesn''t check if your grid is regular...'
2846 0 : call wrtout(std_out,msg,'PERS')
2847 : end if
2848 : end if
2849 :
2850 : !call listkk to get mapping from FBZ to IBZ
2851 : rdum=1.0d-5 ! cutoff distance to decide when two k points match
2852 246 : ABI_MALLOC(dtefield%indkk_f2ibz,(dtefield%fnkpt,6))
2853 :
2854 82 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
2855 :
2856 : !ji: The following may need modification in the future
2857 : !**** no spin-polarization doubling ; allow use of time reversal symmetry ****
2858 :
2859 : !Here is original call
2860 : !
2861 : !call listkk(rdum,gmet,dtefield%indkk_f2ibz,dtset%kptns,dtefield%fkptns,nkpt,&
2862 : !& dtefield%fnkpt,dtset%nsym,1,dtset%symafm,dtset%symrel,1, spaceComm)
2863 :
2864 82 : call timab(1002,2,tsec)
2865 82 : call timab(1003,1,tsec)
2866 :
2867 : call listkk(rdum,gmet,dtefield%indkk_f2ibz,dtset%kptns,dtefield%fkptns,nkpt,&
2868 82 : & dtefield%fnkpt,dtset%nsym,1,dtset%symafm,symrec,1, spaceComm, use_symrec=.True.)
2869 :
2870 82 : call timab(1003,2,tsec)
2871 82 : call timab(1004,1,tsec)
2872 :
2873 : !Construct i2fbz and f2ibz
2874 246 : ABI_MALLOC(dtefield%i2fbz,(nkpt))
2875 82 : idum=0
2876 42439 : do ikpt=1,dtefield%fnkpt
2877 : if (dtefield%indkk_f2ibz(ikpt,2)==1 .and. &
2878 169428 : & dtefield%indkk_f2ibz(ikpt,6) == 0 .and. &
2879 82 : & maxval(abs(dtefield%indkk_f2ibz(ikpt,3:5))) == 0 ) then
2880 4219 : dtefield%i2fbz(dtefield%indkk_f2ibz(ikpt,1))=ikpt
2881 4219 : idum=idum+1
2882 : end if
2883 : end do
2884 82 : if (idum/=nkpt)then
2885 0 : ABI_ERROR('Found wrong number of k-points in IBZ')
2886 : end if
2887 :
2888 : !set flags for fields, forces, stresses
2889 : fieldflag = ( (dtset%berryopt== 4) .or. (dtset%berryopt== 6) .or. (dtset%berryopt== 7) &
2890 82 : & .or. (dtset%berryopt==14) .or. (dtset%berryopt==16) .or. (dtset%berryopt==17) )
2891 : ! following two flags activates computation of projector gradient contributions to force and
2892 : ! stress in finite field PAW calculations
2893 25 : calc_epaw3_force = (fieldflag .and. usepaw == 1 .and. dtset%optforces /= 0)
2894 25 : calc_epaw3_stress = (fieldflag .and. usepaw == 1 .and. dtset%optstress /= 0)
2895 :
2896 :
2897 :
2898 : !----------------------------------------------------------------------------
2899 : !------------- Allocate PAW space if necessary ------------------------------
2900 : !----------------------------------------------------------------------------
2901 :
2902 82 : if (usepaw == 1) then
2903 :
2904 30 : dtefield%usepaw = usepaw
2905 30 : dtefield%natom = natom
2906 30 : dtefield%my_natom = mpi_enreg%my_natom
2907 :
2908 90 : ABI_MALLOC(dtefield%lmn_size,(ntypat))
2909 60 : ABI_MALLOC(dtefield%lmn2_size,(ntypat))
2910 90 : do itypat = 1, ntypat
2911 60 : dtefield%lmn_size(itypat) = pawtab(itypat)%lmn_size
2912 90 : dtefield%lmn2_size(itypat) = pawtab(itypat)%lmn2_size
2913 : end do
2914 :
2915 30 : lmn2_size_max = psps%lmnmax*(psps%lmnmax+1)/2
2916 30 : dtefield%lmn2max = lmn2_size_max
2917 :
2918 : ! expibi and qijb_kk are NOT parallelized over atoms
2919 : ! this may change in the future (JZwanziger 18 March 2014)
2920 150 : ABI_MALLOC(dtefield%qijb_kk,(2,lmn2_size_max,dtefield%natom,3))
2921 120 : ABI_MALLOC(dtefield%expibi,(2,dtefield%natom,3))
2922 30 : dtefield%has_expibi = 1
2923 30 : dtefield%has_qijb = 1
2924 :
2925 30 : if ( fieldflag .and. dtefield%has_rij==0) then
2926 9 : lmn2_size_max = psps%lmnmax*(psps%lmnmax+1)/2
2927 45 : ABI_MALLOC(dtefield%rij,(lmn2_size_max,ntypat,3))
2928 9 : dtefield%has_rij = 1
2929 : end if
2930 :
2931 : ! additional F3-type force term for finite electric field with PAW. Same term
2932 : ! might also apply for other displacement-type field calculations, but not sure yet
2933 : ! JZwanziger 4 April 2014
2934 30 : if ( calc_epaw3_force ) then
2935 36 : ABI_MALLOC(dtefield%epawf3,(dtefield%natom,3,3))
2936 9 : dtefield%has_epawf3 = 1
2937 : end if
2938 30 : if ( calc_epaw3_stress ) then
2939 36 : ABI_MALLOC(dtefield%epaws3,(dtefield%natom,3,6))
2940 9 : dtefield%has_epaws3 = 1
2941 : end if
2942 :
2943 30 : ncpgr = 0
2944 30 : if ( fieldflag .and. dtefield%usecprj == 0) then
2945 27 : ABI_MALLOC(dimlmn,(natom))
2946 9 : call pawcprj_getdim(dimlmn,natom,nattyp_dum,ntypat,typat,pawtab,'R')
2947 : ! allocate space for cprj at kpts in BZ (IBZ or FBZ)
2948 17292 : ABI_MALLOC(dtefield%cprj,(natom, mband*dtset%nspinor*dtset%nkpt*nsppol))
2949 : ! write(std_out,*) "initberry alloc of cprj ", shape(dtefield%cprj)
2950 9 : if (calc_epaw3_force .and. .not. calc_epaw3_stress) ncpgr = 3
2951 9 : if (.not. calc_epaw3_force .and. calc_epaw3_stress) ncpgr = 6
2952 9 : if (calc_epaw3_force .and. calc_epaw3_stress) ncpgr = 9
2953 9 : call pawcprj_alloc(dtefield%cprj,ncpgr,dimlmn)
2954 9 : dtefield%usecprj = 1
2955 9 : ABI_FREE(dimlmn)
2956 : end if
2957 :
2958 120 : ABI_MALLOC(dtefield%cprjindex,(nkpt,nsppol))
2959 2138 : dtefield%cprjindex(:,:) = 0
2960 :
2961 30 : if (dtset%kptopt /= 3) then
2962 104 : ABI_MALLOC(dtefield%atom_indsym,(4,nsym,natom))
2963 26 : call symatm(dtefield%atom_indsym,natom,nsym,symrec,dtset%tnons,tol8,typat,xred)
2964 26 : lmax = psps%mpsang - 1
2965 182 : ABI_MALLOC(dtefield%zarot,(2*lmax+1,2*lmax+1,lmax+1,nsym))
2966 26 : call setsym_ylm(gprimd,lmax,nsym,1,rprimd,symrec,dtefield%zarot)
2967 26 : dtefield%nsym = nsym
2968 26 : dtefield%lmax = lmax
2969 26 : dtefield%lmnmax = psps%lmnmax
2970 : end if
2971 :
2972 : end if
2973 :
2974 : !------------------------------------------------------------------------------
2975 : !------------------- Compute variables related to MPI // ----------------------
2976 : !------------------------------------------------------------------------------
2977 :
2978 :
2979 :
2980 :
2981 82 : if (nproc==1) then
2982 58 : dtefield%fmkmem = dtefield%fnkpt
2983 58 : dtefield%fmkmem_max = dtefield%fnkpt
2984 58 : dtefield%mkmem_max = nkpt
2985 : else
2986 24 : dtefield%fmkmem = 0
2987 25368 : do ikpt = 1, dtefield%fnkpt
2988 25344 : ikpti = dtefield%indkk_f2ibz(ikpt,1)
2989 25344 : nband_k = dtset%nband(ikpti)
2990 25344 : if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,-1,me))) &
2991 8472 : & dtefield%fmkmem = dtefield%fmkmem + 1
2992 : end do
2993 : ! Maximum value of mkmem and fmkmem
2994 24 : call xmpi_max(dtefield%fmkmem,dtefield%fmkmem_max,spaceComm,ierr)
2995 : ! I have to use the dummy variable mkmem_ because
2996 : ! mkmem is declared as intent(in) while the first
2997 : ! argument of xmpi_max must be intent(inout)
2998 24 : mkmem_ = mkmem
2999 24 : call xmpi_max(mkmem_,dtefield%mkmem_max,spaceComm,ierr)
3000 : end if
3001 :
3002 410 : ABI_MALLOC(mpi_enreg%kpt_loc2fbz_sp,(0:nproc-1,1:dtefield%fmkmem_max*nsppol, 1:2))
3003 410 : ABI_MALLOC(mpi_enreg%kpt_loc2ibz_sp,(0:nproc-1,1:dtefield%mkmem_max*nsppol, 1:2))
3004 328 : ABI_MALLOC(mpi_enreg%kptdstrb,(nproc,6,dtefield%fmkmem_max*nsppol*2))
3005 246 : ABI_MALLOC(mpi_enreg%mkmem,(0:nproc-1))
3006 142154 : mpi_enreg%kpt_loc2fbz_sp(:,:,:) = 0
3007 14790 : mpi_enreg%kpt_loc2ibz_sp(:,:,:) = 0
3008 903924 : mpi_enreg%kptdstrb(:,:,:) = 0
3009 220 : mpi_enreg%mkmem(:) = 0
3010 :
3011 82 : if (fieldflag) then
3012 75 : ABI_MALLOC(dtefield%cgqindex,(3,6,nkpt*nsppol))
3013 50 : ABI_MALLOC(dtefield%nneigh,(nkpt))
3014 47449 : dtefield%cgqindex(:,:,:) = 0 ; dtefield%nneigh(:) = 0
3015 : end if
3016 :
3017 82 : pwind_alloc = mpw*dtefield%fmkmem_max
3018 328 : ABI_MALLOC(pwind,(pwind_alloc,2,3))
3019 246 : ABI_MALLOC(pwnsfac,(2,pwind_alloc))
3020 :
3021 : !------------------------------------------------------------------------------
3022 : !---------------------- Compute efield_type variables -------------------------
3023 : !------------------------------------------------------------------------------
3024 :
3025 : !Initialization of efield_type variables
3026 82 : mult=dtset%useria+1
3027 328 : dtefield%efield_dot(:) = zero
3028 1066 : dtefield%dkvecs(:,:) = zero
3029 82 : dtefield%maxnstr = 0 ; dtefield%maxnkstr = 0
3030 574 : dtefield%nstr(:) = 0 ; dtefield%nkstr(:) = 0
3031 328 : ABI_MALLOC(dtefield%ikpt_dk,(dtefield%fnkpt,2,3))
3032 328 : ABI_MALLOC(dtefield%cgindex,(nkpt,nsppol))
3033 164 : ABI_MALLOC(dtefield%kgindex,(nkpt))
3034 246 : ABI_MALLOC(dtefield%fkgindex,(dtefield%fnkpt))
3035 254962 : dtefield%ikpt_dk(:,:,:) = 0
3036 4383 : dtefield%cgindex(:,:) = 0
3037 82 : dtefield%mband_occ = 0
3038 246 : ABI_MALLOC(dtefield%nband_occ,(nsppol))
3039 4301 : dtefield%kgindex(:) = 0
3040 42439 : dtefield%fkgindex(:) = 0
3041 :
3042 82 : if (fieldflag) then
3043 100 : dtset%rfdir(1:3) = 1
3044 : end if
3045 :
3046 :
3047 : !Compute spin degeneracy
3048 82 : if (nsppol == 1 .and. dtset%nspinor == 1) then
3049 76 : dtefield%sdeg = two
3050 6 : else if (nsppol == 2 .or. my_nspinor == 2) then
3051 6 : dtefield%sdeg = one
3052 : end if
3053 :
3054 : !Compute the number of occupied bands and check that
3055 : !it is the same for each k-point
3056 :
3057 82 : index = 0
3058 164 : do isppol = 1, nsppol
3059 82 : dtefield%nband_occ(isppol) = 0
3060 4383 : do ikpt = 1, nkpt
3061 :
3062 4219 : mband_occ_k = 0
3063 4219 : nband_k = dtset%nband(ikpt + (isppol - 1)*nkpt)
3064 :
3065 28373 : do iband = 1, nband_k
3066 24154 : index = index + 1
3067 28373 : if (abs(occ(index) - dtefield%sdeg) < tol8) mband_occ_k = mband_occ_k + 1
3068 : end do
3069 :
3070 4219 : if (fieldflag) then
3071 1824 : if (nband_k /= mband_occ_k) then
3072 : write(msg,'(a,a,a)')&
3073 0 : & ' In a finite electric field, nband must be equal ',ch10,&
3074 0 : & ' to the number of valence bands.'
3075 0 : ABI_ERROR(msg)
3076 : end if
3077 : end if
3078 :
3079 4301 : if (ikpt > 1) then
3080 4137 : if (dtefield%nband_occ(isppol) /= mband_occ_k) then
3081 0 : ABI_ERROR("The number of valence bands is not the same for every k-point of present spin channel")
3082 : end if
3083 : else
3084 82 : dtefield%mband_occ = max(dtefield%mband_occ, mband_occ_k)
3085 82 : dtefield%nband_occ(isppol) = mband_occ_k
3086 : end if
3087 :
3088 : end do ! close loop over ikpt
3089 : end do ! close loop over isppol
3090 :
3091 82 : if (fieldflag) then
3092 175 : ABI_MALLOC(dtefield%smat,(2,dtefield%mband_occ,dtefield%mband_occ,nkpt*nsppol,2,3))
3093 :
3094 3513802 : dtefield%smat(:,:,:,:,:,:) = zero
3095 : end if
3096 :
3097 492 : ABI_MALLOC(dtefield%sflag,(dtefield%mband_occ,nkpt*nsppol,2,3))
3098 170950 : dtefield%sflag(:,:,:,:) = 0
3099 :
3100 : !Compute the location of each wavefunction
3101 :
3102 82 : icg = 0
3103 82 : icprj = 0
3104 : !ikg = 0
3105 164 : do isppol = 1, nsppol
3106 4383 : do ikpt = 1, nkpt
3107 :
3108 4219 : nband_k = dtset%nband(ikpt + (isppol-1)*nkpt)
3109 :
3110 4219 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
3111 :
3112 3043 : dtefield%cgindex(ikpt,isppol) = icg
3113 3043 : npw_k = npwarr(ikpt)
3114 3043 : icg = icg + npw_k*dtefield%nspinor*nband_k
3115 :
3116 3125 : if (usepaw == 1) then
3117 1022 : dtefield%cprjindex(ikpt,isppol) = icprj
3118 1022 : icprj = icprj + dtefield%nspinor*nband_k
3119 : end if
3120 :
3121 : end do
3122 : end do
3123 :
3124 82 : ikg = 0
3125 4301 : do ikpt = 1, nkpt
3126 4219 : if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,me)).and.&
3127 : & (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,nsppol,me))) cycle
3128 :
3129 3043 : npw_k = npwarr(ikpt)
3130 3043 : dtefield%kgindex(ikpt) = ikg
3131 4301 : ikg = ikg + npw_k
3132 : end do
3133 :
3134 : !Need to use dtset%red_efieldbar in the whole code
3135 : !Compute the reciprocal lattice coordinates of the electric field
3136 82 : if (fieldflag) then
3137 :
3138 25 : call metric(gmetllc,gprimdlc,-1,rmetllc,rprimd,ucvol_local)
3139 :
3140 25 : if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7) then
3141 :
3142 100 : do ii=1,3
3143 300 : dtset%red_efieldbar(ii) = dot_product(dtset%efield(:),rprimd(:,ii))
3144 100 : dtefield%efield_dot(ii) = dtset%red_efieldbar(ii)
3145 : end do
3146 :
3147 : ! dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1))
3148 : ! dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2))
3149 : ! dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3))
3150 :
3151 25 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3152 25 : & ' initberry: Reduced electric field (ebar)',ch10,&
3153 50 : & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10
3154 25 : call wrtout(std_out,msg,'COLL')
3155 :
3156 : end if
3157 :
3158 25 : if (dtset%berryopt == 6 .or. dtset%berryopt ==7 ) then
3159 :
3160 0 : do ii=1,3
3161 0 : dtset%red_dfield(ii)= (dot_product(dtset%dfield(:),gprimdlc(:,ii)))*ucvol_local/(4.d0*pi)
3162 : end do
3163 :
3164 0 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3165 0 : & ' initberry: Reduced electric displacement field',ch10,&
3166 0 : & ' red_dfield(1:3) = ',dtset%red_dfield(1:3),ch10
3167 0 : call wrtout(std_out,msg,'COLL')
3168 :
3169 : end if
3170 :
3171 :
3172 25 : if ( dtset%berryopt == 14 ) then
3173 : ! transfer to unreduced electric field.
3174 0 : do idir=1,3
3175 0 : dtset%efield(idir)= dot_product(dtset%red_efieldbar(:),gprimdlc(:,idir))
3176 0 : dtefield%efield_dot(idir) = dtset%red_efieldbar(idir)
3177 : ! dtefield%efield2(idir)=dtset%red_efieldbar(idir)
3178 : end do
3179 :
3180 : ! dtefield%efield_dot(1) = dtset%red_efieldbar(1)
3181 : ! dtefield%efield_dot(2) = dtset%red_efieldbar(2)
3182 : ! dtefield%efield_dot(3) = dtset%red_efieldbar(3)
3183 :
3184 0 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3185 0 : & ' initberry: Unreduced electric field (a.u.)',ch10,&
3186 0 : & ' efield(1:3) = ',dtset%efield(1:3),ch10
3187 0 : call wrtout(std_out,msg,'COLL')
3188 :
3189 0 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3190 0 : & ' initberry: Reduced electric field (ebar)',ch10,&
3191 0 : & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10
3192 0 : call wrtout(std_out,msg,'COLL')
3193 :
3194 : end if
3195 :
3196 :
3197 25 : if ( dtset%berryopt == 16 ) then
3198 :
3199 : ! to calculate D
3200 0 : do ii=1,3
3201 0 : dtset%dfield(ii) =(4*pi/ucvol_local)*dot_product(dtset%red_dfield(:),rprimd(:,ii))
3202 : end do
3203 :
3204 0 : do idir=1,3
3205 0 : dtset%efield(idir)= (4*pi/ucvol_local)*dot_product(dtset%red_efield(:),rprimd(:,idir))
3206 : end do
3207 :
3208 0 : do idir=1,3
3209 0 : dtset%red_efieldbar(idir)= (4*pi/ucvol_local)*dot_product(dtset%red_efield(:),rmetllc(:,idir))
3210 0 : dtefield%efield_dot(idir) = dtset%red_efieldbar(idir)
3211 : end do
3212 :
3213 : ! dtefield%efield_dot(1) = dtset%red_efieldbar(1)
3214 : ! dtefield%efield_dot(2) = dtset%red_efieldbar(2)
3215 : ! dtefield%efield_dot(3) = dtset%red_efieldbar(3)
3216 :
3217 :
3218 0 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3219 0 : & ' initberry: Unreduced electric displacement field (a.u.)',ch10,&
3220 0 : & ' dfield(1:3) = ',dtset%dfield(1:3),ch10
3221 0 : call wrtout(std_out,msg,'COLL')
3222 :
3223 0 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3224 0 : & ' initberry: Unreduced electric field (a.u.)',ch10,&
3225 0 : & ' efield(1:3) = ',dtset%efield(1:3),ch10
3226 0 : call wrtout(std_out,msg,'COLL')
3227 :
3228 0 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3229 0 : & ' initberry: Reduced electric field (ebar)',ch10,&
3230 0 : & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10
3231 0 : call wrtout(std_out,msg,'COLL')
3232 :
3233 : end if
3234 :
3235 25 : if ( dtset%berryopt ==17) then
3236 :
3237 : ! to calculate D
3238 :
3239 0 : do idir=1,3
3240 0 : dtset%efield(idir)= dot_product(dtset%red_efieldbar(:),gprimdlc(:,idir)) ! from ebar
3241 0 : dtset%dfield(idir) =(4*pi/ucvol_local)*dot_product(dtset%red_dfield(:),rprimd(:,idir))
3242 : ! dtset%red_efield(idir) = (ucvol_local/(4*pi))*dot_product(dtset%red_efieldbar(:),gmetllc(:,idir))
3243 0 : dtefield%efield_dot(idir) = dtset%red_efieldbar(idir)
3244 : end do
3245 :
3246 0 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3247 0 : & ' initberry: Reduced electric field (ebar)',ch10,&
3248 0 : & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10
3249 0 : call wrtout(std_out,msg,'COLL')
3250 :
3251 :
3252 0 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3253 0 : & ' initberry: Unreduced electric field (a.u.)',ch10,&
3254 0 : & ' efield(1:3) = ',dtset%efield(1:3),ch10
3255 0 : call wrtout(std_out,msg,'COLL')
3256 :
3257 0 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3258 0 : & ' initberry: Reduced electric displacement field (a.u.)',ch10,&
3259 0 : & ' red_dfield(1:3) = ',dtset%red_dfield(1:3),ch10
3260 0 : call wrtout(std_out,msg,'COLL')
3261 :
3262 0 : write(msg,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3263 0 : & ' initberry: Unreduced electric displacement field (a.u.)',ch10,&
3264 0 : & ' dfield(1:3) = ',dtset%dfield(1:3),ch10
3265 0 : call wrtout(std_out,msg,'COLL')
3266 :
3267 :
3268 : end if
3269 :
3270 :
3271 :
3272 : end if
3273 :
3274 82 : call timab(1004,2,tsec)
3275 :
3276 : !------------------------------------------------------------------------------
3277 : !---------------------- Compute dk --------------------------------------------
3278 : !------------------------------------------------------------------------------
3279 :
3280 82 : call timab(1005,1,tsec)
3281 :
3282 328 : do idir = 1, 3
3283 :
3284 246 : if (dtset%rfdir(idir) == 1) then
3285 :
3286 : ! Compute dk(:), the vector between a k-point and its nearest
3287 : ! neighbour along the direction idir
3288 :
3289 229 : dk(:) = zero
3290 229 : dk(idir) = 1._dp ! 1 mean there is no other k-point un the direction idir
3291 126106 : do ikpt = 2, dtefield%fnkpt
3292 503508 : diffk(:) = abs(dtefield%fkptns(:,ikpt) - dtefield%fkptns(:,1))
3293 125877 : if ((diffk(1) < dk(1)+tol8).and.(diffk(2) < dk(2)+tol8).and.&
3294 597 : & (diffk(3) < dk(3)+tol8)) dk(:) = diffk(:)
3295 : end do
3296 916 : dtefield%dkvecs(:,idir) = dk(:)
3297 : ! DEBUG
3298 : ! write(std_out,*)' initberry : idir, dk', idir, dk
3299 : ! ENDDEBUG
3300 :
3301 : ! For each k point, find k_prim such that k_prim= k + dk mod(G)
3302 : ! where G is a vector of the reciprocal lattice
3303 :
3304 126335 : do ikpt = 1, dtefield%fnkpt
3305 :
3306 : ! First k+dk, then k-dk
3307 252441 : do isign=-1,1,2
3308 252212 : kpt_shifted1=dtefield%fkptns(1,ikpt)- isign*dk(1)
3309 252212 : kpt_shifted2=dtefield%fkptns(2,ikpt)- isign*dk(2)
3310 252212 : kpt_shifted3=dtefield%fkptns(3,ikpt)- isign*dk(3)
3311 : ! Note that this is still a order fnkpt**2 algorithm.
3312 : ! It is possible to implement a order fnkpt algorithm, see listkk.F90.
3313 205852374 : do ikpt1 = 1, dtefield%fnkpt
3314 205726268 : diffk1=dtefield%fkptns(1,ikpt1) - kpt_shifted1
3315 205726268 : if(abs(diffk1-nint(diffk1))>tol8)cycle
3316 13642852 : diffk2=dtefield%fkptns(2,ikpt1) - kpt_shifted2
3317 13642852 : if(abs(diffk2-nint(diffk2))>tol8)cycle
3318 1036356 : diffk3=dtefield%fkptns(3,ikpt1) - kpt_shifted3
3319 1036356 : if(abs(diffk3-nint(diffk3))>tol8)cycle
3320 252212 : dtefield%ikpt_dk(ikpt,(isign+3)/2,idir) = ikpt1
3321 205726268 : exit
3322 : end do ! ikpt1
3323 : end do ! isign
3324 :
3325 : ! OLD CODING
3326 : ! First: k + dk
3327 : ! do ikpt1 = 1, dtefield%fnkpt
3328 : ! diffk(:) = abs(dtefield%fkptns(:,ikpt1) - &
3329 : ! & dtefield%fkptns(:,ikpt) - dk(:))
3330 : ! if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
3331 : ! dtefield%ikpt_dk(ikpt,1,idir) = ikpt1
3332 : ! exit
3333 : ! end if
3334 : ! end do
3335 :
3336 : ! Second: k - dk
3337 : ! do ikpt1 = 1, dtefield%fnkpt
3338 : ! diffk(:) = abs(dtefield%fkptns(:,ikpt1) - &
3339 : ! & dtefield%fkptns(:,ikpt) + dk(:))
3340 : ! if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
3341 : ! dtefield%ikpt_dk(ikpt,2,idir) = ikpt1
3342 : ! exit
3343 : ! end if
3344 : ! end do
3345 :
3346 : end do ! ikpt
3347 :
3348 : ! Find the string length, starting from k point 1
3349 : ! (all strings must have the same number of points)
3350 :
3351 229 : nkstr = 1
3352 229 : ikpt1 = 1
3353 985 : do ikpt = 1, dtefield%fnkpt
3354 985 : ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir)
3355 985 : if (ikpt1 == 1) exit
3356 985 : nkstr = nkstr + 1
3357 : end do
3358 :
3359 : ! Check that the string length is a divisor of nkpt
3360 229 : if(mod(dtefield%fnkpt,nkstr) /= 0) then
3361 : write(msg,'(a,i5,a,i7)')&
3362 0 : & ' The string length = ',nkstr,&
3363 0 : & ', is not a divisor of fnkpt =',dtefield%fnkpt
3364 0 : ABI_BUG(msg)
3365 : end if
3366 :
3367 229 : dtefield%nkstr(idir) = nkstr
3368 229 : dtefield%nstr(idir) = dtefield%fnkpt/nkstr
3369 :
3370 : end if ! dtset%rfdir(idir) == 1
3371 :
3372 : write(msg,'(a,i1,a,i3,a,i6)')&
3373 246 : & ' initberry: for direction ',idir,', nkstr = ',dtefield%nkstr(idir),&
3374 492 : & ', nstr = ',dtefield%nstr(idir)
3375 246 : call wrtout(std_out,msg,'COLL')
3376 328 : call wrtout(ab_out,msg,'COLL')
3377 :
3378 : end do ! close loop over idir
3379 :
3380 82 : call timab(1005,2,tsec)
3381 82 : call timab(1006,1,tsec)
3382 :
3383 328 : dtefield%maxnstr = maxval(dtefield%nstr(:))
3384 328 : dtefield%maxnkstr = maxval(dtefield%nkstr(:))
3385 410 : ABI_MALLOC(dtefield%idxkstr,(dtefield%maxnkstr,dtefield%maxnstr,3))
3386 146290 : dtefield%idxkstr(:,:,:) = 0
3387 :
3388 : !for the geometry of the string space :
3389 328 : ABI_MALLOC(dtefield%coord_str,(2,dtefield%maxnstr,3))
3390 328 : ABI_MALLOC(dtefield%str_neigh,(-2:2,dtefield%maxnstr,3))
3391 410 : ABI_MALLOC(dtefield%strg_neigh,(-2:2,dtefield%maxnstr,2,3))
3392 56893 : dtefield%coord_str(:,:,:) = 0.d0
3393 113458 : dtefield%str_neigh(:,:,:)=0
3394 227080 : dtefield%strg_neigh(:,:,:,:)=0
3395 1804 : dtefield%gmet_str(:,:,:)=0.d0
3396 :
3397 : !------------------------------------------------------------------------------
3398 : !---------------------- Build the strings -------------------------------------
3399 : !------------------------------------------------------------------------------
3400 :
3401 246 : ABI_MALLOC(kpt_mark,(dtefield%fnkpt))
3402 328 : do idir = 1, 3
3403 :
3404 328 : if (dtset%rfdir(idir) == 1) then
3405 :
3406 126335 : iunmark = 1
3407 126335 : kpt_mark(:) = 0
3408 18836 : do istr = 1, dtefield%nstr(idir)
3409 :
3410 123606 : do while(kpt_mark(iunmark) /= 0)
3411 104999 : iunmark = iunmark + 1
3412 : end do
3413 18607 : dtefield%idxkstr(1,istr,idir) = iunmark
3414 18607 : kpt_mark(iunmark) = 1
3415 126335 : do ikstr = 2, dtefield%nkstr(idir)
3416 107499 : ikpt1 = dtefield%idxkstr(ikstr-1,istr,idir)
3417 107499 : ikpt2 = dtefield%ikpt_dk(ikpt1,1,idir)
3418 107499 : dtefield%idxkstr(ikstr,istr,idir) = ikpt2
3419 126106 : kpt_mark(ikpt2) = 1
3420 : end do
3421 :
3422 : end do ! istr
3423 :
3424 : ! compute distance between strings
3425 : ! compute the metric matrix of the strings space in the direction idir
3426 916 : do ii = 1,3
3427 2977 : do jj = 1,3
3428 2061 : if (ii<idir.and.jj<idir) dtefield%gmet_str(ii ,jj ,idir) = &
3429 377 : & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir)
3430 2061 : if (ii<idir.and.jj>idir) dtefield%gmet_str(ii ,jj-1,idir) = &
3431 77 : & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir)
3432 2061 : if (ii>idir.and.jj<idir) dtefield%gmet_str(ii-1,jj ,idir) = &
3433 77 : & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir)
3434 2061 : if (ii>idir.and.jj>idir) dtefield%gmet_str(ii-1,jj-1,idir) = &
3435 1072 : & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir)
3436 : end do
3437 : end do
3438 : ! DEBUG
3439 : ! write(std_out,*)'gmet'
3440 : ! do ii=1,3
3441 : ! write(std_out,*)gmet(ii,:)
3442 : ! end do
3443 : ! write(std_out,*)'gmet_str'
3444 : ! do ii=1,2
3445 : ! write(std_out,*)dtefield%gmet_str(ii,:,idir)
3446 : ! end do
3447 : ! ENDDEBUG
3448 18836 : do istr = 1, dtefield%nstr(idir)
3449 74657 : do ii = 1,3
3450 55821 : if (ii<idir) dtefield%coord_str(ii,istr,idir)=dtefield%fkptns(ii,dtefield%idxkstr(1,istr,idir))
3451 74428 : if (ii>idir) dtefield%coord_str(ii-1,istr,idir)=dtefield%fkptns(ii,dtefield%idxkstr(1,istr,idir))
3452 : end do
3453 : end do
3454 :
3455 : ! the following is very similar to getshell
3456 : dist_ = 0._dp
3457 687 : do ii = 1,2
3458 687 : dist_ = dist_ + dtefield%gmet_str(ii,ii,idir)
3459 : end do
3460 229 : max_dist = 2._dp * dist_ * 2._dp
3461 :
3462 1603 : dk_str(:,:,idir) = 0._dp
3463 : last_dist = 0._dp
3464 : ! ishell = 0
3465 : ! dtefield%str_neigh(:,:,:) = 0
3466 : dk_flag = 0
3467 521 : do while (dk_flag /= 2)
3468 : ! Advance shell counter
3469 : ! ishell = ishell + 1
3470 :
3471 : ! Search the smallest distance between two strings
3472 : dist = max_dist
3473 20588 : do istr = 1,dtefield%nstr(idir)
3474 60888 : delta_str3(:) = dtefield%coord_str(:,1,idir) - dtefield%coord_str(:,istr,idir)
3475 81476 : do coord1 = -1,1 !two loop to search also on the border of the BZ
3476 263848 : do coord2 = -1,1
3477 547992 : dist_ = 0._dp
3478 547992 : dstr(:) = delta_str3(:) - nint(delta_str3(:))
3479 182664 : dstr(1) = dstr(1) + real(coord1,dp)
3480 182664 : dstr(2) = dstr(2) + real(coord2,dp)
3481 547992 : do ii = 1,2
3482 1278648 : do jj = 1,2
3483 1095984 : dist_ = dist_ + dstr(ii)*dtefield%gmet_str(ii,jj,idir)*dstr(jj)
3484 : end do
3485 : end do
3486 243552 : if ((dist_ < dist).and.(dist_ - last_dist > tol8)) then
3487 182664 : dist = dist_
3488 : end if
3489 : end do
3490 : end do
3491 : end do
3492 :
3493 5920 : last_dist = dist
3494 :
3495 : ! search the connecting vectors for that distance
3496 6149 : do istr = 1,dtefield%nstr(idir)
3497 17571 : delta_str3(:) = dtefield%coord_str(:,istr,idir) - dtefield%coord_str(:,1,idir)
3498 22970 : do coord1 = -1,1
3499 68902 : do coord2 = -1,1
3500 155367 : dist_ = 0._dp
3501 155367 : dstr(:) = delta_str3(:) - nint(delta_str3(:))
3502 51789 : dstr(1) = dstr(1) + real(coord1,dp)
3503 51789 : dstr(2) = dstr(2) + real(coord2,dp)
3504 155367 : do ii = 1,2
3505 362523 : do jj = 1,2
3506 310734 : dist_ = dist_ + dstr(ii)*dtefield%gmet_str(ii,jj,idir)*dstr(jj)
3507 : end do
3508 : end do
3509 68902 : if (abs(dist_ - dist) < tol8) then
3510 728 : if (dk_flag == 0) then
3511 687 : dk_str(:,1,idir) = dstr(:)
3512 : dk_flag = 1
3513 : ! DEBUG
3514 : ! write(std_out,'(a,i4,2e15.4)')'1st connect', istr, dstr
3515 : ! ENDDEBUG
3516 499 : elseif (dk_str(1,1,idir)*dstr(2)-dk_str(2,1,idir)*dstr(1) > tol8) then
3517 687 : dk_str(:,2,idir) = dstr(:)
3518 : dk_flag = 2
3519 : ! DEBUG
3520 : ! write(std_out,'(a,i4,2e15.4)')'2nd connect', istr, dstr
3521 : ! ENDDEBUG
3522 : exit
3523 : end if
3524 : end if
3525 : end do
3526 22970 : if (dk_flag == 2) exit
3527 : end do
3528 5920 : if (dk_flag == 2) exit
3529 : end do
3530 :
3531 : end do ! do while
3532 :
3533 : ! search the two neighbours for each string
3534 18836 : do istr = 1,dtefield%nstr(idir)
3535 18607 : dtefield%str_neigh(0,istr,idir) = istr
3536 55821 : dtefield%strg_neigh(0,istr,:,idir) = 0
3537 3615495 : do jstr = 1,dtefield%nstr(idir)
3538 10789977 : delta_str3(:) = dtefield%coord_str(:,jstr,idir) - dtefield%coord_str(:,istr,idir)
3539 14405243 : do coord1 = -1,1
3540 46756567 : do coord2 = -1,1
3541 97109793 : dist_ = 0._dp
3542 97109793 : dstr(:) = delta_str3(:) - nint(delta_str3(:))
3543 32369931 : dstr(1) = dstr(1) + real(coord1,dp)
3544 32369931 : dstr(2) = dstr(2) + real(coord2,dp)
3545 107899770 : do ii = 1,2
3546 226589517 : if (sum(abs(dstr(:)-dk_str(:,ii,idir)))<tol8) then
3547 37214 : dtefield%str_neigh(ii,istr,idir) = jstr
3548 37214 : dtefield%strg_neigh(ii,istr,1,idir) = coord1
3549 37214 : dtefield%strg_neigh(ii,istr,2,idir) = coord2
3550 194107944 : elseif (sum(abs(dstr(:)+dk_str(:,ii,idir)))<tol8) then
3551 37214 : dtefield%str_neigh(-ii,istr,idir) = jstr
3552 37214 : dtefield%strg_neigh(-ii,istr,1,idir) = coord1
3553 37214 : dtefield%strg_neigh(-ii,istr,2,idir) = coord2
3554 : end if
3555 : end do
3556 : end do
3557 : end do
3558 : end do
3559 : end do
3560 :
3561 : ! DEBUG
3562 : ! write(std_out,'(a,e15.4,e15.4,e15.4,e15.4)')'dk_str',dk_str(1,1,idir),dk_str(2,1,idir),dk_str(1,2,idir),dk_str(2,2,idir)
3563 : ! write(std_out,*)'istr, neigh1, strg(1,:), neigh2, strg(2,:),neigh-1, strg(-1,:), neigh-2, strg(-2,:)'
3564 : ! do istr=1,dtefield%nstr(idir)
3565 : ! write(std_out,'(13i4)')istr, &
3566 : ! & dtefield%str_neigh(1,istr,idir), dtefield%strg_neigh(1,istr,:,idir),&
3567 : ! & dtefield%str_neigh(2,istr,idir), dtefield%strg_neigh(2,istr,:,idir),&
3568 : ! & dtefield%str_neigh(-1,istr,idir), dtefield%strg_neigh(-1,istr,:,idir),&
3569 : ! & dtefield%str_neigh(-2,istr,idir), dtefield%strg_neigh(-2,istr,:,idir)
3570 : ! end do
3571 : ! ENDDEBUG
3572 :
3573 :
3574 : end if ! rfdir(idir) == 1
3575 :
3576 : end do ! close loop over idir
3577 :
3578 82 : ABI_FREE(kpt_mark)
3579 :
3580 82 : call timab(1006,2,tsec)
3581 82 : call timab(1007,1,tsec)
3582 :
3583 : !------------------------------------------------------------------------------
3584 : !------------ Compute PAW on-site terms if necessary --------------------------
3585 : !------------------------------------------------------------------------------
3586 :
3587 82 : if (usepaw == 1 .and. dtefield%has_expibi == 1) then
3588 90 : ABI_MALLOC(calc_expibi,(2,natom))
3589 120 : do idir = 1, 3
3590 360 : dk = dtefield%dkvecs(1:3,idir)
3591 630 : calc_expibi = zero
3592 90 : call expibi(calc_expibi,dk,natom,xred)
3593 660 : dtefield%expibi(1:2,1:natom,idir) = calc_expibi
3594 : end do
3595 : ! call expibi(dtefield%expibi,dtefield%dkvecs,natom,xred)
3596 30 : dtefield%has_expibi = 2
3597 30 : ABI_FREE(calc_expibi)
3598 : end if
3599 :
3600 82 : if (usepaw == 1 .and. dtefield%has_qijb == 1) then
3601 120 : ABI_MALLOC(calc_qijb,(2,dtefield%lmn2max,natom))
3602 :
3603 120 : do idir = 1, 3
3604 360 : dk = dtefield%dkvecs(1:3,idir)
3605 19710 : calc_qijb = zero
3606 : call qijb_kk(calc_qijb,dk,dtefield%expibi(1:2,1:natom,idir),&
3607 90 : & gprimd,dtefield%lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat)
3608 19740 : dtefield%qijb_kk(1:2,1:dtefield%lmn2max,1:natom,idir) = calc_qijb
3609 : ! call qijb_kk(dtefield%qijb_kk,dtefield%dkvecs,dtefield%expibi,&
3610 : ! & gprimd,dtefield%lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat)
3611 : end do
3612 30 : dtefield%has_qijb = 2
3613 30 : ABI_FREE(calc_qijb)
3614 : end if
3615 :
3616 82 : if (usepaw == 1 .and. dtefield%has_rij == 1) then
3617 9 : c1=sqrt(four_pi/three)
3618 27 : do itypat = 1, ntypat
3619 675 : do klmn = 1, pawtab(itypat)%lmn2_size
3620 648 : dtefield%rij(klmn,itypat,1) = c1*pawtab(itypat)%qijl(4,klmn) ! S_{1,1} ~ x
3621 648 : dtefield%rij(klmn,itypat,2) = c1*pawtab(itypat)%qijl(2,klmn) ! S_{1,-1} ~ y
3622 666 : dtefield%rij(klmn,itypat,3) = c1*pawtab(itypat)%qijl(3,klmn) ! S_{1,0} ~ z
3623 : end do ! end loop over klmn
3624 : end do ! end loop over itypat
3625 9 : dtefield%has_rij = 2
3626 : end if !
3627 :
3628 82 : call timab(1007,2,tsec)
3629 82 : call timab(1008,1,tsec)
3630 :
3631 : !------------------------------------------------------------------------------
3632 : !------------ Build the array pwind that is needed to compute the -------------
3633 : !------------ overlap matrices at k +- dk -------------
3634 : !------------------------------------------------------------------------------
3635 :
3636 82 : ecut_eff = dtset%ecut*(dtset%dilatmx)**2
3637 82 : exchn2n3d = 0 ; istwf_k = 1 ; ikg1 = 0
3638 21373282 : pwind(:,:,:) = 0
3639 3562159 : pwnsfac(1,:) = 1.0_dp
3640 3562159 : pwnsfac(2,:) = 0.0_dp
3641 246 : ABI_MALLOC(kg1_k,(3,mpw))
3642 :
3643 82 : ipwnsfac = 0
3644 :
3645 328 : do idir = 1, 3
3646 :
3647 328 : if (dtset%rfdir(idir) == 1) then
3648 :
3649 916 : dk(:) = dtefield%dkvecs(:,idir)
3650 :
3651 687 : do ifor = 1, 2
3652 :
3653 1145 : if (ifor == 2) dk(:) = -1._dp*dk(:)
3654 :
3655 : ! Build pwind and kgindex
3656 : ! NOTE: The array kgindex is important for parallel execution.
3657 : ! In case nsppol = 2, it may happent that a particular processor
3658 : ! treats k-points at different spin polarizations.
3659 : ! In this case, it is not possible to address the elements of
3660 : ! pwind correctly without making use of the kgindex array.
3661 :
3662 458 : ikg = 0 ; ikpt_loc = 0 ; isppol = 1
3663 252670 : do ikpt = 1, dtefield%fnkpt
3664 :
3665 252212 : ikpti = dtefield%indkk_f2ibz(ikpt,1)
3666 252212 : nband_k = dtset%nband(ikpti)
3667 252212 : ikpt1f = dtefield%ikpt_dk(ikpt,ifor,idir)
3668 252212 : ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1)
3669 :
3670 252212 : if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,1,me)).and.&
3671 : & (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,nsppol,me))) cycle
3672 :
3673 20479080 : ikpt_loc = ikpt_loc + 1
3674 :
3675 : ! Build basis sphere of plane waves for the nearest neighbour of
3676 : ! the k-point (important for MPI //)
3677 :
3678 81463812 : kg1_k(:,:) = 0
3679 603344 : kpt1(:) = dtset%kptns(:,ikpt1i)
3680 : call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt,istwf_k,kg1_k,kpt1,&
3681 150836 : & 1,mpi_enreg,mpw,npw_k1)
3682 150836 : me_g0=mpi_enreg%me_g0
3683 :
3684 :
3685 : ! ji: fkgindex is defined here !
3686 150836 : dtefield%fkgindex(ikpt) = ikg
3687 :
3688 : !
3689 : ! Deal with symmetry transformations
3690 : !
3691 :
3692 : ! bra k-point k(b) and IBZ k-point kIBZ(b) related by
3693 : ! k(b) = alpha(b) S(b)^t kIBZ(b) + G(b)
3694 : ! where alpha(b), S(b) and G(b) are given by indkk_f2ibz
3695 : !
3696 : ! For the ket k-point:
3697 : ! k(k) = alpha(k) S(k)^t kIBZ(k) + G(k) - GBZ(k)
3698 : ! where GBZ(k) takes k(k) to the BZ
3699 : !
3700 :
3701 150836 : isym = dtefield%indkk_f2ibz(ikpt,2)
3702 150836 : isym1 = dtefield%indkk_f2ibz(ikpt1f,2)
3703 :
3704 : ! Construct transformed G vector that enters the matching condition:
3705 : ! alpha(k) S(k)^{t,-1} ( -G(b) - GBZ(k) + G(k) )
3706 :
3707 : dg(:) = -dtefield%indkk_f2ibz(ikpt,3:5) &
3708 : & -nint(-dtefield%fkptns(:,ikpt) - dk(:) - tol10 + &
3709 : & dtefield%fkptns(:,ikpt1f)) &
3710 603344 : & +dtefield%indkk_f2ibz(ikpt1f,3:5)
3711 :
3712 : ! old code
3713 : ! iadum(:)=0
3714 : ! do idum=1,3
3715 : ! iadum(:)=iadum(:)+ symrec(:,idum,isym1)*dg(idum)
3716 : ! end do
3717 :
3718 : ! new code
3719 150836 : iadum(:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),dg(:))
3720 :
3721 150836 : dg(:) = iadum(:)
3722 :
3723 361304 : if ( dtefield%indkk_f2ibz(ikpt1f,6) == 1 ) dg(:) = -dg(:)
3724 :
3725 : ! Construct S(k)^{t,-1} S(b)^{t}
3726 :
3727 14329420 : dum33(:,:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),symrec(:,:,isym))
3728 :
3729 : ! Construct alpha(k) alpha(b)
3730 :
3731 150836 : if (dtefield%indkk_f2ibz(ikpt,6) == dtefield%indkk_f2ibz(ikpt1f,6)) then
3732 : itrs=0
3733 : else
3734 47684 : itrs=1
3735 : end if
3736 :
3737 :
3738 150836 : npw_k = npwarr(ikpti)
3739 : ! npw_k1 = npwarr(ikpt1i)
3740 :
3741 : ! loop over bra G vectors
3742 19603392 : do ipw = 1, npw_k
3743 :
3744 : ! NOTE: the bra G vector is taken for the sym-related IBZ k point,
3745 : ! not for the FBZ k point
3746 77810224 : iadum(:) = kg(:,dtefield%kgindex(ikpti) + ipw)
3747 :
3748 : ! Store non-symmorphic operation phase factor exp[i2\pi \alpha G \cdot t]
3749 :
3750 19452556 : if ( ipwnsfac == 0 ) then
3751 : ! old code
3752 3346561 : rdum=0.0_dp
3753 13386244 : do idum=1,3
3754 13386244 : rdum=rdum+dble(iadum(idum))*dtset%tnons(idum,isym)
3755 : end do
3756 3346561 : rdum=two_pi*rdum
3757 3346561 : if ( dtefield%indkk_f2ibz(ikpt,6) == 1 ) rdum=-rdum
3758 3346561 : pwnsfac(1,ikg+ipw) = cos(rdum)
3759 3346561 : pwnsfac(2,ikg+ipw) = sin(rdum)
3760 : !
3761 : ! new code
3762 : ! rdum = DOT_PRODUCT(dble(iadum(:)),dtset%tnons(:,isym))
3763 : ! rdum= two_pi*rdum
3764 : ! if ( dtefield%indkk_f2ibz(ikpt,6) == 1 ) rdum=-rdum
3765 : ! pwnsfac(1,ikg+ipw) = cos(rdum)
3766 : ! pwnsfac(2,ikg+ipw) = sin(rdum)
3767 :
3768 : end if
3769 :
3770 : ! to determine r.l.v. matchings, we transformed the bra vector
3771 : ! Rotation
3772 19452556 : iadum1(:)=0
3773 77810224 : do idum1=1,3
3774 252883228 : iadum1(:)=iadum1(:)+dum33(:,idum1)*iadum(idum1)
3775 : end do
3776 19452556 : iadum(:)=iadum1(:)
3777 : ! Time reversal
3778 38253460 : if (itrs==1) iadum(:)=-iadum(:)
3779 : ! Translation
3780 77810224 : iadum(:) = iadum(:) + dg(:)
3781 :
3782 2672045062 : do jpw = 1, npw_k1
3783 10684269576 : iadum1(1:3) = kg1_k(1:3,jpw)
3784 : if ( (iadum(1) == iadum1(1)).and. &
3785 2671067394 : & (iadum(2) == iadum1(2)).and. &
3786 826832 : & (iadum(3) == iadum1(3)) ) then
3787 18625724 : pwind(ikg + ipw,ifor,idir) = jpw
3788 : ! write(std_out,'(a,2x,3i4,2x,i4)') 'Found !:',iadum1(:),jpw
3789 18625724 : exit
3790 : end if
3791 : end do
3792 : end do
3793 :
3794 151294 : ikg = ikg + npw_k
3795 :
3796 : end do ! close loop over ikpt
3797 :
3798 687 : ipwnsfac = 1
3799 :
3800 : end do ! close loop over ifor
3801 :
3802 : end if ! rfdir(idir) == 1
3803 :
3804 : end do ! close loop over idir
3805 :
3806 :
3807 82 : call timab(1008,2,tsec)
3808 82 : call timab(1009,1,tsec)
3809 :
3810 : !Build mpi_enreg%kptdstrb
3811 : !array required to communicate the WFs between cpus in berryphase_new.f
3812 : !(MPI // over k-points)
3813 82 : if (nproc>1) then
3814 96 : do idir = 1, 3
3815 96 : if (dtset%rfdir(idir) == 1) then
3816 216 : do ifor = 1, 2
3817 :
3818 144 : ikpt_loc = 0
3819 360 : do isppol = 1, nsppol
3820 :
3821 152352 : do ikpt = 1, dtefield%fnkpt
3822 :
3823 152064 : ikpti = dtefield%indkk_f2ibz(ikpt,1)
3824 152064 : nband_k = dtset%nband(ikpti)
3825 152064 : ikpt1f = dtefield%ikpt_dk(ikpt,ifor,idir)
3826 152064 : ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1)
3827 :
3828 152064 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle
3829 :
3830 50688 : ikpt_loc = ikpt_loc + 1
3831 : mpi_enreg%kptdstrb(me + 1,ifor+2*(idir-1),ikpt_loc) = &
3832 50688 : & ikpt1i + (isppol - 1)*nkpt
3833 :
3834 : mpi_enreg%kptdstrb(me+1,ifor+2*(idir-1),&
3835 : & ikpt_loc+dtefield%fmkmem_max*nsppol) = &
3836 152208 : & ikpt1f + (isppol - 1)*dtefield%fnkpt
3837 :
3838 : end do ! ikpt
3839 : end do ! isppol
3840 : end do ! ifor
3841 : end if ! dtset%rfdir(idir) == 1
3842 : end do ! idir
3843 : end if ! nproc>1
3844 :
3845 : !build mpi_enreg%kpt_loc2fbz_sp
3846 82 : ikpt_loc = 0
3847 164 : do isppol = 1, nsppol
3848 42521 : do ikpt = 1, dtefield%fnkpt
3849 :
3850 42357 : ikpti = dtefield%indkk_f2ibz(ikpt,1)
3851 42357 : nband_k = dtset%nband(ikpti)
3852 :
3853 42357 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle
3854 :
3855 25461 : ikpt_loc = ikpt_loc + 1
3856 :
3857 25461 : mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 1) = ikpt
3858 42439 : mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 2) = isppol
3859 :
3860 : end do
3861 : end do
3862 :
3863 :
3864 : !parallel case only :
3865 : !build mpi_enreg%kpt_loc2ibz_sp, dtefield%cgqindex and dtefield%nneigh
3866 82 : if ((fieldflag).and.(nproc>1)) then
3867 12 : ikpt_loc = 0
3868 24 : do isppol = 1, nsppol
3869 1368 : do ikpt = 1, nkpt
3870 :
3871 1344 : ikptf = dtefield%i2fbz(ikpt)
3872 1344 : nband_k = dtset%nband(ikpti)
3873 :
3874 1344 : neigh(:) = 0 ; icg = 0 ; ikg = 0 ; flag_kpt = 0; icprj = 0
3875 5376 : do idir=1, 3
3876 :
3877 : ! skip idir values for which efield_dot(idir) = 0
3878 4032 : if (abs(dtefield%efield_dot(idir)) < tol12) cycle
3879 :
3880 13440 : do ifor = 1, 2
3881 :
3882 8064 : flag = 0
3883 :
3884 8064 : ikpt1f = dtefield%ikpt_dk(ikptf,ifor,idir)
3885 8064 : ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1)
3886 :
3887 8064 : dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ikg
3888 8064 : ikg = ikg + npwarr(ikpt1i)
3889 :
3890 : ! check if this neighbour is also a previous neighbour
3891 33456 : do ineigh = 1, (ifor+2*(idir-1))
3892 33456 : if (neigh(ineigh) == ikpt1i) then
3893 912 : flag = 1
3894 912 : dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ineigh
3895 : dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = &
3896 912 : & dtefield%cgqindex(2,ineigh,ikpt+(isppol-1)*nkpt)
3897 : exit
3898 : end if
3899 : end do
3900 : ! create the cgqindex of the neighbour if necessary
3901 4032 : if (flag == 0) then
3902 7152 : neigh(ifor+2*(idir-1)) = ikpt1i
3903 : dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = &
3904 7152 : & ifor+2*(idir-1)
3905 7152 : dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = icg
3906 7152 : if (isppol == 1) dtefield%nneigh(ikpt) = dtefield%nneigh(ikpt) + 1
3907 7152 : icg = icg + npwarr(ikpt1i)*dtefield%nspinor*nband_k
3908 : end if
3909 : end do !ifor
3910 : end do !idir
3911 :
3912 1356 : if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))) then
3913 : ! ikpt is one of my kpt_loc
3914 448 : ikpt_loc = ikpt_loc + 1
3915 448 : mpi_enreg%kpt_loc2ibz_sp(me, ikpt_loc, 1) = ikpt
3916 448 : mpi_enreg%kpt_loc2ibz_sp(me, ikpt_loc, 2) = isppol
3917 : end if
3918 :
3919 : end do !ikpt
3920 : end do !isppol
3921 : end if !nproc>1
3922 :
3923 : !should be temporary
3924 : !unassigned mpi_enreg%kpt_loc2fbz_sp are empty ; inform other cpu (there are better ways...)
3925 82 : mpi_enreg%mkmem(me) = mkmem
3926 : !do ii=ikpt_loc+1,dtefield%fmkmem_max
3927 : !mpi_enreg%kpt_loc2fbz_sp(me, ii, 1) = -1
3928 : !end do
3929 :
3930 :
3931 : !(same as mpi_enreg%kptdstrb but for k-points in the iBZ),
3932 : !dtefield%cgqindex and dtefield%nneigh
3933 :
3934 82 : if ((fieldflag).and.(nproc>1)) then
3935 :
3936 12 : ikpt_loc = 1
3937 24 : do isppol = 1, nsppol
3938 1368 : do ikpt = 1, nkpt
3939 :
3940 1344 : nband_k = dtset%nband(ikpt)
3941 1344 : ikptf = dtefield%i2fbz(ikpt)
3942 :
3943 1344 : neigh(:) = 0 ; icg = 0 ; ikg = 0 ; flag_kpt = 0; icprj = 0
3944 5376 : do idir = 1, 3
3945 :
3946 : ! Skip idir values for which efield_dot(idir) = 0
3947 4032 : if (abs(dtefield%efield_dot(idir)) < tol12 .and. (fieldflag)) cycle
3948 :
3949 1344 : do ifor = 1, 2
3950 :
3951 : ikpt1f = dtefield%ikpt_dk(ikptf,ifor,idir)
3952 : ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1)
3953 :
3954 : ! dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ikg
3955 : ikg = ikg + npwarr(ikpt1i)
3956 :
3957 : flag = 0
3958 : do ineigh = 1, (ifor+2*(idir-1))
3959 : if (neigh(ineigh) == ikpt1i) then
3960 : flag = 1
3961 : ! dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ineigh
3962 : ! dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = &
3963 : ! & dtefield%cgqindex(2,ineigh,ikpt+(isppol-1)*nkpt)
3964 : exit
3965 : end if
3966 : end do
3967 : if (flag == 0) then
3968 : ! neigh(ifor+2*(idir-1)) = ikpt1i
3969 : ! dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = &
3970 : ! & ifor+2*(idir-1)
3971 : ! dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = icg
3972 : ! if (isppol == 1) dtefield%nneigh(ikpt) = dtefield%nneigh(ikpt) + 1
3973 : ! icg = icg + npwarr(ikpt1i)*dtset%nspinor*nband_k
3974 : end if
3975 :
3976 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
3977 :
3978 : flag_kpt = 1
3979 :
3980 : ! MVeithen: the if condition allows to avoid that the same wavefunction
3981 : ! is send several times to a particular cpu
3982 :
3983 : end do ! ifor
3984 : end do ! idir
3985 :
3986 12 : if (flag_kpt == 1) ikpt_loc = ikpt_loc + 1
3987 :
3988 : end do ! ikpt
3989 : end do ! isppol
3990 :
3991 : end if ! fieldflag
3992 :
3993 82 : call xmpi_sum(mpi_enreg%kptdstrb,spaceComm,ierr)
3994 82 : call xmpi_sum(mpi_enreg%kpt_loc2fbz_sp,spaceComm,ierr)
3995 82 : if (fieldflag) then
3996 25 : call xmpi_sum(mpi_enreg%kpt_loc2ibz_sp,spaceComm,ierr)
3997 25 : call xmpi_sum(mpi_enreg%mkmem,spaceComm,ierr)
3998 : end if
3999 :
4000 : !------------------------------------------------------------------------------
4001 : !------------------------ Estimate critical field -----------------------------
4002 : !------------------------------------------------------------------------------
4003 :
4004 : !Compute the minimal value of the bandgap required to be below
4005 : !the critical field as defined by the relation
4006 : !| E_i*a_i | < E_g/n_i
4007 :
4008 : if (fieldflag) then
4009 :
4010 100 : do idir = 1, 3
4011 : ! eg_dir(idir) = abs(dtefield%efield_dot(idir))*dtefield%nkstr(idir)
4012 100 : eg_dir(idir) = abs(dtset%red_efieldbar(idir))*dtefield%nkstr(idir)
4013 : end do
4014 :
4015 :
4016 125 : eg = maxval(eg_dir)
4017 25 : eg_ev = eg*Ha_eV
4018 :
4019 25 : if (dtset%optcell ==0 .and. (dtset%berryopt == 4 .or. dtset%berryopt == 14)) then
4020 25 : write(msg,'(a,a,a,a,a,a,a,a,f7.2,a,a)')ch10,&
4021 25 : & ' initberry: COMMENT - ',ch10,&
4022 25 : & ' As a rough estimate,',ch10,&
4023 25 : & ' to be below the critical field, the bandgap of your system',ch10,&
4024 50 : & ' should be larger than ',eg_ev,' eV.',ch10
4025 25 : call wrtout(ab_out,msg,'COLL')
4026 25 : call wrtout(std_out,msg,'COLL')
4027 :
4028 : else
4029 :
4030 0 : write(msg,'(a,a,a,a,a,a,a)') ch10,&
4031 0 : & ' initberry: COMMENT - ',ch10,&
4032 0 : & ' The estimation of critical electric field should be checked after calculation.',ch10,&
4033 0 : & ' It is printed out just after total energy.' ,ch10
4034 :
4035 0 : call wrtout(ab_out,msg,'COLL')
4036 0 : call wrtout(std_out,msg,'COLL')
4037 : end if
4038 : end if
4039 :
4040 82 : ABI_FREE(kg1_k)
4041 :
4042 82 : call timab(1009,2,tsec)
4043 82 : call timab(1001,2,tsec)
4044 :
4045 : DBG_EXIT("COLL")
4046 :
4047 164 : end subroutine initberry
4048 : !!***
4049 :
4050 150836 : end module m_berryphase_new
4051 : !!***
|
