Branch data 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 : 52830 : subroutine update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,&
1724 : 52830 : & efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,&
1725 : 52830 : & mkmem,mpi_enreg,mpw,my_natom,natom,nattyp,ngfft,nkpt,npwarr,ntypat,&
1726 : 52830 : & pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,pwind,&
1727 : 52830 : & pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,&
1728 : : & scfcv_level,scfcv_quit,scfcv_step,ucvol,unit_out,&
1729 : 52830 : & 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 : 52830 : real(dp),allocatable :: cg13(:,:,:),ph1d(:,:)
1766 : 52830 : type(pawcprj_type),allocatable :: cprj(:,:)
1767 : :
1768 : : ! *************************************************************************
1769 : :
1770 : 52830 : 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 [ + + ]: 52830 : & dtset%berryopt ==16 .or. &
1777 : : & dtset%berryopt ==17 ) efield = .true.
1778 [ + - + + ]: 191 : calc_epaw3_force = ( efield .and. dtset%optforces /= 0 .and. usepaw == 1 )
1779 [ + + + - : 52830 : calc_epaw3_stress = ( efield .and. dtset%optstress /= 0 .and. usepaw == 1 )
+ + ]
1780 : :
1781 [ + + ]: 52830 : usecprj=1; if (psps%usepaw==0) usecprj = 0
1782 : 52830 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
1783 : 52830 : mcprj=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol/mpi_enreg%nproc_band
1784 : : !Just to be sure
1785 [ + + - + ]: 52830 : 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 : 52830 : ncpgr = 0
1791 : 52830 : ctocprj_choice = 1 ! no derivs
1792 [ + + + + ]: 52830 : 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 [ - + ]: 52761 : ABI_MALLOC(cprj,(0,0))
1816 : : end if ! end update of cprj
1817 : :
1818 [ + + ]: 52830 : 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 [ + - + + : 87902 : ABI_FREE(cprj)
- + - + -
+ ]
2196 : :
2197 [ - + ]: 52830 : 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 : 5662 : subroutine init_e_field_vars(dtefield,dtset,gmet,gprimd,kg,&
2616 : 5662 : & mpi_enreg,npwarr,occ,pawang,pawrad,pawtab,psps,&
2617 : 5662 : & 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 : 5662 : initfield = .false.
2644 : :
2645 : : !initialization
2646 : 5662 : dtefield%has_qijb = 0
2647 : 5662 : dtefield%has_epawf3 = 0
2648 : 5662 : dtefield%has_epaws3 = 0
2649 : 5662 : dtefield%has_expibi = 0
2650 : 5662 : dtefield%has_rij = 0
2651 : 5662 : dtefield%usecprj = 0
2652 : 5662 : dtefield%berryopt = 0
2653 : :
2654 : : if ((dtset%berryopt < 0).or.(dtset%berryopt == 4) .or. (dtset%berryopt == 6) .or.(dtset%berryopt == 7) .or. &
2655 [ + + ]: 5662 : & (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 [ + - ]: 5580 : if (.not. initfield .and. dtset%orbmag == 0) then
2667 : : ! initorbmag.F90 also allocates pwind and pwnsfac
2668 : 5580 : pwind_alloc = 1
2669 [ - + ]: 5580 : ABI_MALLOC(pwind,(pwind_alloc,2,3))
2670 [ - + ]: 5580 : ABI_MALLOC(pwnsfac,(2,pwind_alloc))
2671 [ + + + + : 89280 : pwind(:,:,:)=0
+ + ]
2672 [ + + + + ]: 27900 : pwnsfac(:,:)=zero
2673 : : end if
2674 : :
2675 : 5662 : 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 : : !!***
|
