Line data Source code
1 : !!****m* ABINIT/m_d2frnl
2 : !! NAME
3 : !! m_d2frnl
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GM, AR, MB, MT, AM)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : ! nvtx related macro definition
23 : #include "nvtx_macros.h"
24 :
25 : module m_d2frnl
26 :
27 : use defs_basis
28 : use m_abicore
29 : use m_xmpi
30 : use m_mpinfo
31 : use m_errors
32 : use m_cgtools
33 : use m_nctk
34 : use m_hamiltonian
35 : use m_efmas_defs
36 : use m_wfk
37 : use m_dtset
38 : use m_dtfil
39 :
40 :
41 : use defs_datatypes, only : pseudopotential_type
42 : use defs_abitypes, only : MPI_type
43 : use m_time, only : timab
44 : use m_geometry, only : metric, strconv
45 : use m_efmas, only : check_degeneracies
46 : use m_io_tools, only : file_exists
47 : use m_hdr, only : hdr_skip
48 : use m_pawang, only : pawang_type
49 : use m_pawrad, only : pawrad_type
50 : use m_pawtab, only : pawtab_type,pawtab_get_lsize
51 : use m_pawfgrtab,only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
52 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags
53 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_free, pawrhoij_gather, &
54 : pawrhoij_nullify, pawrhoij_symrhoij
55 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_copy, pawcprj_free
56 : use m_pawdij, only : pawdijfr
57 : use m_paw_dfpt, only : pawgrnl
58 : use m_kg, only : mkkin, mkkpg
59 : use m_mkffnl, only : mkffnl
60 : use m_nonlop, only : nonlop
61 : use m_paw_occupancies, only : pawaccrhoij
62 : use m_gemm_nonlop_projectors, only : set_gemm_nonlop_ikpt, gemm_nonlop_use_gemm
63 :
64 : #if defined(HAVE_GPU_MARKERS)
65 : use m_nvtx_data
66 : #endif
67 :
68 : implicit none
69 :
70 : private
71 : !!***
72 :
73 : public :: d2frnl
74 : !!***
75 :
76 : contains
77 : !!***
78 :
79 : !!****f* ABINIT/d2frnl
80 : !! NAME
81 : !! d2frnl
82 : !!
83 : !! FUNCTION
84 : !! Compute the frozen-wavefunction non-local contribution for response functions
85 : !! (strain and/or phonon)
86 : !!
87 : !! INPUTS
88 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)=<G|Cnk>=Fourier coefficients of WF
89 : !! dtfil <type(datafiles_type)>=variables related to files
90 : !! dtset <type(dataset_type)>=all input variables for this dataset
91 : !! dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
92 : !! dyfr_nondiag=1 if dyfrnl is non diagonal with respect to atoms; 0 otherwise
93 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
94 : !! gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere
95 : !! has_allddk= True if all ddk file are present on disk
96 : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
97 : !! kg(3,mpw*mkmem)=work array for coordinates of G vectors in basis
98 : !! primitive translations
99 : !! mgfftf=maximum size of 1D FFTs for the fine FFT grid (PAW)
100 : !! mpi_enreg=information about MPI parallelization
101 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
102 : !! my_natom=number of atoms treated by current processor
103 : !! natom=number of atoms in unit cell
104 : !! nfftf= -PAW ONLY- number of FFT grid points for the fine grid
105 : !! (nfftf=nfft for norm-conserving potential runs)
106 : !! ngfft(18)=contain all needed information about 3D FFT,
107 : !! see ~abinit/doc/variables/vargs.htm#ngfft
108 : !! ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid
109 : !! (ngs_rbzfftf=ngfft for norm-conserving potential runs)
110 : !! npwarr(nkpt)=number of planewaves at each k point, and boundary
111 : !! ntypat=integer specification of atom type (1, 2, ...)
112 : !! occ(mband*nkpt*nsppol)=occupation numbers of bands (usually 2) at each k point
113 : !! rfphon=1 if non local contribution of dynamical matrix have to be computed
114 : !! rfstrs!=0 if non local contribution of elastic tensor have to be computed
115 : !! paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
116 : !! pawang <type(pawang_type)>=paw angular mesh and related data
117 : !! pawbec= flag for the computation of Born Effective Charge within PAW ; set to 1 if yes
118 : !! pawpiezo= flag for the computation of piezoelectric tensor within PAW ; set to 1 if yes
119 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
120 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
121 : !! ph1d(2,3*(2*mgfft+1)*natom)=phase information related to structure factor
122 : !! ph1df(2,3*(2*mgfftf+1)*natom)=phase information related to structure factor on the fine FFT grid (PAW)
123 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
124 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
125 : !! symrec(3,3,nsym)=symmetries in reciprocal space (dimensionless)
126 : !! vtrial(nfftf,nspden)=total potential (Hartree+XC+loc)
127 : !! vxc(nfftf,nspden)=XC potential
128 : !! xred(3,natom)=reduced coordinates of atoms (dimensionless)
129 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
130 : !! ylmgr(mpw*mkmem,9,mpsang*mpsang*useylm)= gradients of real spherical harmonics for each G and k point
131 : !!
132 : !! OUTPUT
133 : !! becfrnl(3,natom,3*pawbec)=NL frozen contribution to Born Effective Charges (PAW only)
134 : !! (3,natom) = derivative wr to the displ. of one atom in one direction
135 : !! (3) = derivative wr to electric field in one direction
136 : !! piezofrnl(3,6*pawpiezo)=NL frozen contribution to piezoelectric tensor (PAW only)
137 : !! dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=
138 : !! non-symmetrized non-local contribution to the dynamical matrix
139 : !! If NCPP, it depends on one atom
140 : !! If PAW, it depends on two atoms
141 : !! eltfrnl(6+3*natom,6)=non-symmetrized non-local contribution to the
142 : !! elastic tensor
143 : !!
144 : !! SIDE EFFECTS
145 : !! ===== if psps%usepaw==1
146 : !! pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
147 : !! pawfgrtab(:)%gylmgr2 are deallocated here
148 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
149 : !! (gradients of rhoij for each atom with respect to atomic positions are computed here)
150 : !!
151 : !! SOURCE
152 :
153 487 : subroutine d2frnl(becfrnl,cg,dtfil,dtset,dyfrnl,dyfr_cplex,dyfr_nondiag,efmasdeg,efmasval,eigen,eltfrnl,&
154 487 : & gsqcut,has_allddk,indsym,kg,mband_mem_rbz,mkmem_rbz,mgfftf,mpi_enreg,mpsang,my_natom,natom,nfftf,ngfft,ngfftf,npwarr,&
155 487 : & occ,paw_ij,pawang,pawbec,pawfgrtab,pawpiezo,pawrad,pawrhoij,pawtab,ph1d,ph1df,piezofrnl,psps,&
156 487 : & rprimd,rfphon,rfstrs,symrec,vtrial,vxc,xred,ylm,ylmgr)
157 :
158 : !Arguments ------------------------------------
159 : !scalars
160 : integer,intent(in) :: dyfr_cplex,dyfr_nondiag,mgfftf,mpsang,my_natom,natom
161 : integer,intent(in) :: nfftf,pawbec,pawpiezo,rfphon,rfstrs
162 : integer,intent(in) :: mkmem_rbz,mband_mem_rbz
163 : real(dp),intent(in) :: gsqcut
164 : type(MPI_type),intent(in) :: mpi_enreg
165 : type(datafiles_type),intent(in) :: dtfil
166 : type(dataset_type),intent(in) :: dtset
167 : type(pawang_type),intent(in) :: pawang
168 : type(pseudopotential_type),intent(in) :: psps
169 : !arrays
170 : integer,intent(in) :: indsym(4,dtset%nsym,natom),kg(3,dtset%mpw*mkmem_rbz)
171 : integer,intent(in) :: ngfft(18),ngfftf(18),npwarr(dtset%nkpt)
172 : integer,intent(in) :: symrec(3,3,dtset%nsym)
173 : real(dp),intent(in) :: cg(2,dtset%mpw*dtset%nspinor*mband_mem_rbz*mkmem_rbz*dtset%nsppol)
174 : real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
175 : real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
176 : real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*natom)
177 : real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*natom),rprimd(3,3)
178 : real(dp),intent(in) :: vxc(nfftf,dtset%nspden),xred(3,natom)
179 : real(dp),intent(in) :: ylm(dtset%mpw*mkmem_rbz,mpsang*mpsang*psps%useylm)
180 : real(dp),intent(in) :: ylmgr(dtset%mpw*mkmem_rbz,9,mpsang*mpsang*psps%useylm)
181 : real(dp),intent(in),target :: vtrial(nfftf,dtset%nspden)
182 : real(dp),intent(out) :: becfrnl(3,natom,3*pawbec),piezofrnl(6,3*pawpiezo)
183 : real(dp),intent(out) :: dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)
184 : real(dp),intent(out) :: eltfrnl(6+3*natom,6)
185 : logical,intent(inout):: has_allddk
186 : type(efmasdeg_type),allocatable,intent(out):: efmasdeg(:)
187 : type(efmasval_type),allocatable,intent(out):: efmasval(:,:)
188 : type(paw_ij_type),intent(in) :: paw_ij(my_natom)
189 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
190 : type(pawrad_type),intent(in) :: pawrad(psps%ntypat)
191 : type(pawrhoij_type),intent(inout),target :: pawrhoij(my_natom*psps%usepaw)
192 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat)
193 :
194 : !Local variables-------------------------------
195 : !scalars
196 : integer,parameter :: formeig1=1,usecprj=0
197 : integer :: bandmin,bandmax,bdtot_index,bufdim
198 : integer :: choice_bec2,choice_bec54,choice_efmas,choice_phon,choice_strs,choice_piez3,choice_piez55
199 : integer :: cplex,cplx,cpopt,cpopt_bec,ddkcase,deg_dim
200 : integer :: dimffnl,dimffnl_str,dimnhat,ia,iatom,iashift,iband,jband,ibg,icg,icplx,ideg,ider,idir
201 : integer :: ider_str,idir_ffnl,idir_str,ielt,ieltx,ierr,ii,ikg,ikpt,ilm,ipw,iq,iq0
202 : integer :: ispinor,isppol,istwf_k,isub,itypat,jj,jsub,klmn,master,me,mu
203 : integer :: my_comm_atom,n1,n2,n3,nband_k,ncpgr,nfftot,ngrhoij,nkpg,nnlout_bec1,nnlout_bec2,nnlout_efmas
204 : integer :: nnlout_piez1,nnlout_piez2,nnlout_phon,nnlout_strs,npw_,npw_k,nsp,nsploop,nu
205 : integer :: optgr,optgr2,option,option_rhoij,optstr,optstr2,paw_opt,paw_opt_1,paw_opt_3,paw_opt_efmas
206 : integer :: shift_rhoij,signs,signs_field,spaceworld,sz2,sz3,tim_nonlop
207 : integer :: iband_, iband_me, jband_me, nband_me
208 : real(dp) :: arg,enl,enlk,ucvol,wtk_k
209 : logical :: has_ddk_file,need_becfr,need_efmas,need_piezofr,paral_atom,t_test
210 : logical :: use_timerev,use_zeromag
211 : character(len=500) :: msg
212 487 : type(gs_hamiltonian_type) :: gs_ham
213 : !arrays
214 : integer :: ik_ddk(3),ddkfil(3)
215 974 : integer :: bands_treated_now(dtset%mband), band_procs(dtset%mband)
216 : integer :: ndat,idat
217 487 : integer,allocatable :: dimlmn(:),kg_k(:,:),l_size_atm(:)
218 487 : integer,pointer :: my_atmtab(:)
219 : real(dp) :: dotprod(2),dummy(0),gmet(3,3),gprimd(3,3),grhoij(3),kpoint(3),nonlop_dum(1,1)
220 : real(dp) :: rmet(3,3),tsec(2)
221 487 : complex(dp), allocatable :: ch2c_tmp(:)
222 487 : real(dp),allocatable :: becfrnl_tmp(:,:,:),becfrnlk(:,:,:),becij(:,:,:,:,:),cg_left(:,:)
223 487 : real(dp),allocatable :: occ_k(:),eig_k(:)
224 974 : real(dp),allocatable :: cwavef(:,:),ddk(:,:,:),ddkinpw(:,:,:),dyfrnlk(:,:)
225 487 : real(dp),allocatable :: elt_work(:,:),eltfrnlk(:,:),enlout_bec1(:),enlout_bec2(:),enlout_efmas(:)
226 487 : real(dp),allocatable :: enlout_piez1(:),enlout_piez2(:),enlout_phon(:),enlout_strs(:)
227 487 : real(dp),allocatable :: gh2c(:,:),gs2c(:,:)
228 487 : real(dp),allocatable :: kpg_k(:,:),mpibuf(:),nhat_dum(:,:),piezofrnlk(:,:),ph3d(:,:,:)
229 487 : real(dp),allocatable :: svectout(:,:),ylm_k(:,:),ylmgr_k(:,:,:)
230 487 : real(dp),allocatable,target :: ffnl(:,:,:,:),ffnl_str(:,:,:,:,:)
231 : character(len=fnlen) :: fiwfddk(3)
232 487 : type(paw_ij_type),allocatable :: paw_ij_tmp(:)
233 487 : type(pawcprj_type),allocatable,target :: cwaveprj(:,:)
234 487 : type(pawfgrtab_type),allocatable :: pawfgrtab_tmp(:)
235 487 : type(pawrhoij_type),pointer :: pawrhoij_tot(:)
236 47726 : type(wfk_t) :: ddkfiles(3)
237 : ! *************************************************************************
238 :
239 : DBG_ENTER("COLL")
240 :
241 487 : call timab(159,1,tsec)
242 : ABI_NVTX_START_RANGE(NVTX_D2FRNL)
243 :
244 487 : write(msg,'(3a)')ch10,' ==> Calculation of the frozen part of the second order derivatives, this can take some time...',ch10
245 487 : call wrtout(std_out,msg,'COLL')
246 :
247 : !Set up parallelism
248 487 : spaceworld=mpi_enreg%comm_cell
249 487 : me=mpi_enreg%me_kpt
250 487 : master=0
251 487 : paral_atom=(my_natom/=natom)
252 487 : my_comm_atom=mpi_enreg%comm_atom
253 487 : my_atmtab=>mpi_enreg%my_atmtab
254 487 : ndat = dtset%bandpp
255 :
256 : !Compute gmet, gprimd and ucvol from rprimd
257 487 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
258 :
259 : !If needed, check for ddk files (used for effective charges)
260 487 : if (pawbec==1.or.pawpiezo==1) then
261 51 : ddkfil(:)=0
262 204 : do ii=1,3
263 153 : ddkcase=ii+natom*3
264 153 : call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk(ii))
265 153 : t_test = file_exists(fiwfddk(ii))
266 : ! Trick needed to run Abinit test suite in netcdf mode.
267 153 : if (.not. t_test .and. file_exists(nctk_ncify(fiwfddk(ii)))) then
268 81 : t_test = .True.; fiwfddk(ii) = nctk_ncify(fiwfddk(ii))
269 81 : write(msg,"(3a)")"- File: ",trim(fiwfddk(ii))," does not exist but found netcdf file with similar name."
270 234 : call wrtout(std_out,msg,'COLL')
271 : end if
272 132 : if (t_test) ddkfil(ii)=20+ii ! Note the use of unit numbers 21, 22 and 23
273 : end do
274 123 : has_ddk_file=(any(ddkfil(:)>0))
275 132 : has_allddk =(all(ddkfil(:)>0))
276 : else
277 436 : has_ddk_file=.FALSE.
278 436 : has_allddk =.FALSE.
279 : end if
280 :
281 487 : if(pawbec==1.or.pawpiezo==1.and.has_ddk_file) then
282 51 : if(.not.has_allddk) then
283 24 : write(msg,'(5a)')ch10,&
284 24 : & ' WARNING: All ddk perturbations are needed to compute',ch10,&
285 48 : & ' the frozen part of Born effective charges and/or piezoelectric tensor.',ch10
286 24 : call wrtout(std_out,msg,'COLL')
287 : else
288 27 : write(msg,'(5a)')ch10,&
289 27 : & ' All ddk perturbations are available.',ch10,&
290 54 : & ' The frozen part of Born effective charges and/or piezoelectric tensor will be computed',ch10
291 27 : call wrtout(std_out,msg,'COLL')
292 : end if
293 : end if
294 :
295 487 : need_becfr=(pawbec==1.and.has_ddk_file)
296 487 : need_piezofr=(pawpiezo==1.and.has_ddk_file)
297 :
298 : !Initialization of frozen non local array
299 487 : if(rfphon==1) then
300 32302 : dyfrnl(:,:,:,:,:)=zero
301 1350 : ABI_MALLOC(dyfrnlk,(6,natom))
302 : end if
303 487 : if(rfstrs/=0)then
304 3283 : eltfrnl(:,:)=zero;enl=zero
305 185 : ABI_MALLOC(eltfrnlk,(6+3*natom,6))
306 : end if
307 487 : if (need_becfr) then
308 948 : becfrnl(:,:,:)=zero
309 108 : ABI_MALLOC(becfrnlk,(3,natom,3))
310 : end if
311 487 : if (need_piezofr) then
312 396 : piezofrnl(:,:)=zero
313 18 : ABI_MALLOC(piezofrnlk,(6,3))
314 : end if
315 487 : need_efmas=dtset%efmas>0
316 487 : if(need_efmas.and.(rfphon==1.or.rfstrs/=0.or.need_becfr.or.need_piezofr)) then
317 0 : write(msg,'(5a)')ch10,&
318 0 : & ' ERROR: Efmas calculation is incompatible with phonons, elastic tensor, Born effective charges,',ch10,&
319 0 : & ' and piezoelectric tensor calculations. Please revise your input.',ch10
320 0 : ABI_ERROR(msg)
321 : end if
322 :
323 1461 : ABI_MALLOC(eig_k,(ndat))
324 974 : ABI_MALLOC(occ_k,(ndat))
325 :
326 : !Common initialization
327 487 : bdtot_index=0;ibg=0;icg=0
328 487 : nsploop=dtset%nsppol;if (dtset%nspden==4) nsploop=4
329 487 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
330 487 : nfftot=ngfftf(1)*ngfftf(2)*ngfftf(3)
331 :
332 : !Common data for "nonlop" routine
333 487 : tim_nonlop=6
334 974 : signs=1 ; signs_field = 2 ; eig_k=zero ; idir=0
335 : ! ffnl are in cartesian coordinates for EFMAS (idir==4),
336 : ! in contrast to reduced coordinates for the other responses (idir==0).
337 487 : idir_ffnl=0 ; if(need_efmas) idir_ffnl=4
338 487 : choice_phon=0;choice_strs=0
339 487 : if(rfphon==1)then
340 450 : shift_rhoij=0
341 450 : choice_phon=4
342 450 : nnlout_phon=max(1,6*natom)
343 1350 : ABI_MALLOC(enlout_phon,(nnlout_phon*ndat))
344 : end if
345 487 : if(rfstrs/=0)then
346 37 : shift_rhoij=6
347 37 : choice_strs=6
348 37 : nnlout_strs=6*(3*natom+6)
349 111 : ABI_MALLOC(enlout_strs,(nnlout_strs*ndat))
350 : end if
351 487 : if (psps%usepaw==0) then
352 430 : paw_opt=0 ; cpopt=-1
353 : else
354 57 : paw_opt=2 ; cpopt=1+2*usecprj
355 : end if
356 487 : if(need_piezofr)then
357 18 : choice_piez3 = 3
358 18 : choice_piez55 = 55
359 18 : nnlout_piez1 = 6
360 18 : nnlout_piez2 = 36
361 18 : paw_opt_1 = 1
362 18 : paw_opt_3 = 3
363 54 : ABI_MALLOC(enlout_piez1,(nnlout_piez1*ndat))
364 54 : ABI_MALLOC(enlout_piez2,(nnlout_piez2*ndat))
365 : end if
366 487 : if (need_becfr) then
367 27 : choice_bec2=2 ; choice_bec54=54
368 27 : nnlout_bec1=max(1,3*natom) ; nnlout_bec2=max(1,18*natom);
369 27 : paw_opt_1=1 ; paw_opt_3=3 ; cpopt_bec=-1
370 81 : ABI_MALLOC(enlout_bec1,(nnlout_bec1*ndat))
371 81 : ABI_MALLOC(enlout_bec2,(nnlout_bec2*ndat))
372 : else
373 460 : choice_bec2=0 ; choice_bec54=0
374 460 : nnlout_bec1=0
375 : end if
376 487 : if(need_efmas) then
377 17 : ABI_MALLOC(enlout_efmas,(0))
378 75 : ABI_MALLOC(efmasdeg,(dtset%nkpt))
379 430 : ABI_MALLOC(efmasval,(dtset%mband,dtset%nkpt))
380 : end if
381 :
382 : !Initialize Hamiltonian (k-independent terms)
383 : call gs_ham%init(psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,natom,&
384 : & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
385 : & paw_ij=paw_ij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
386 487 : & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,gpu_option=dtset%gpu_option)
387 :
388 : !===== PAW specific section
389 487 : if (psps%usepaw==1) then
390 :
391 : ! Define several sizes & flags
392 57 : ncpgr=0;ngrhoij=0
393 57 : if(rfphon==1)then
394 44 : ncpgr=3;ngrhoij=3
395 : end if
396 57 : if(rfphon==1.or.need_becfr)then
397 51 : ncpgr=6;ngrhoij=6
398 : end if
399 57 : if(rfstrs/=0.and.rfphon==1)then
400 12 : ncpgr=9;ngrhoij=9
401 : end if
402 57 : if(rfstrs/=0.or.need_piezofr)then
403 22 : ncpgr=9;ngrhoij=9
404 : end if
405 :
406 : ! If PAW and Born Eff. Charges, one has to compute some additional data:
407 : ! For each atom and for electric field direction k:
408 : ! becij(k)=<Phi_i|r_k-R_k|Phi_j>-<tPhi_i|r_k-R_k|tPhi_j> + sij.R_k
409 57 : if (need_becfr.or.need_piezofr) then
410 162 : ABI_MALLOC(becij,(gs_ham%dimekb1,gs_ham%dimekb2,dtset%nspinor**2,1,3))
411 11868 : becij=zero
412 123 : ABI_MALLOC(paw_ij_tmp,(my_natom))
413 123 : ABI_MALLOC(pawfgrtab_tmp,(my_natom))
414 27 : call paw_ij_nullify(paw_ij_tmp)
415 27 : cplex=1;nsp=1 ! Force nsppol/nspden to 1 because Dij^(1) due to electric field is spin-independent
416 : call paw_ij_init(paw_ij_tmp,cplex,dtset%nspinor,nsp,nsp,dtset%pawspnorb,natom,psps%ntypat,&
417 27 : & dtset%typat,pawtab,has_dijfr=1,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab )
418 27 : call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,mpi_atmtab=my_atmtab)
419 : call pawfgrtab_init(pawfgrtab_tmp,1,l_size_atm,dtset%nspden,dtset%typat,&
420 27 : & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
421 27 : ABI_FREE(l_size_atm)
422 108 : do ii=1,3 ! Loop over direction of electric field
423 81 : call paw_ij_reset_flags(paw_ij_tmp,all=.True.)
424 : call pawdijfr(gprimd,ii,natom+2,my_natom,natom,nfftf,ngfftf,nsp,nsp,psps%ntypat,&
425 : & 0,paw_ij_tmp,pawang,pawfgrtab_tmp,pawrad,pawtab,cplex,&
426 : & (/zero,zero,zero/),rprimd,ucvol,vtrial,vtrial,vxc,xred,&
427 81 : & comm_atom=my_comm_atom, mpi_atmtab=my_atmtab ) ! vtrial not used here
428 189 : do isppol=1,dtset%nspinor**2
429 162 : call pawdij2e1kb(paw_ij_tmp(:),nsp,my_comm_atom,e1kbfr=becij(:,:,:,:,ii),mpi_atmtab=my_atmtab)
430 : end do
431 : end do
432 27 : call paw_ij_free(paw_ij_tmp)
433 27 : call pawfgrtab_free(pawfgrtab_tmp)
434 69 : ABI_FREE(paw_ij_tmp)
435 69 : ABI_FREE(pawfgrtab_tmp)
436 : end if
437 :
438 : ! PAW occupancies: need to communicate when paral atom is activated
439 57 : if (paral_atom) then
440 170 : ABI_MALLOC(pawrhoij_tot,(natom))
441 26 : call pawrhoij_nullify(pawrhoij_tot)
442 26 : call pawrhoij_gather(pawrhoij,pawrhoij_tot,-1,my_comm_atom)
443 : else
444 31 : pawrhoij_tot => pawrhoij
445 : end if
446 :
447 : ! Projected WF (cprj) and PAW occupancies (& gradients)
448 455 : ABI_MALLOC(cwaveprj,(natom,dtset%nspinor*ndat))
449 57 : call pawcprj_alloc(cwaveprj,ncpgr,gs_ham%dimcprj)
450 215 : do iatom=1,natom
451 158 : sz2=pawrhoij_tot(iatom)%cplex_rhoij*pawrhoij_tot(iatom)%qphase*pawrhoij_tot(iatom)%lmn2_size
452 158 : sz3=pawrhoij_tot(iatom)%nspden
453 790 : ABI_MALLOC(pawrhoij_tot(iatom)%grhoij,(ngrhoij,sz2,sz3))
454 158 : pawrhoij_tot(iatom)%ngrhoij=ngrhoij
455 62545 : pawrhoij_tot(iatom)%grhoij=zero
456 : end do
457 57 : use_timerev=(dtset%kptopt>0.and.dtset%kptopt<3)
458 110 : use_zeromag=(pawrhoij_tot(1)%nspden==4.and.dtset%nspden==1)
459 :
460 : else
461 430 : ABI_MALLOC(cwaveprj,(0,0))
462 : end if !PAW
463 :
464 : !If needed, manage ddk files
465 : !Open ddk WF file(s) in sequential mode
466 487 : if (need_becfr.or.need_piezofr) then
467 108 : do ii=1,3 ! Loop over elect. field directions
468 108 : if (ddkfil(ii)/=0) then
469 81 : write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk(ii))
470 243 : call wrtout([std_out, ab_out], msg)
471 81 : call ddkfiles(ii)%open_read(fiwfddk(ii),formeig1,dtset%iomode,ddkfil(ii), xmpi_comm_self)
472 : end if
473 : end do
474 : end if
475 :
476 : !LOOP OVER SPINS
477 1009 : do isppol=1,dtset%nsppol
478 :
479 : ! Continue to initialize the Hamiltonian (PAW DIJ coefficients)
480 522 : call gs_ham%load_spin(isppol,with_nonlocal=.true.)
481 :
482 : ! Rewind (k+G) data if needed
483 522 : ikg=0
484 :
485 : ! Loop over k points
486 43560 : do ikpt=1,dtset%nkpt
487 42551 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
488 42551 : nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me)
489 : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,dtset%mband,&
490 42551 : & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
491 :
492 42551 : istwf_k=dtset%istwfk(ikpt)
493 42551 : npw_k=npwarr(ikpt)
494 42551 : wtk_k=dtset%wtk(ikpt)
495 170204 : kpoint(:)=dtset%kptns(:,ikpt)
496 :
497 : ! Skip this k-point if not the proper processor
498 42551 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
499 10944 : bdtot_index=bdtot_index+nband_k
500 10944 : cycle
501 : end if
502 : ABI_NVTX_START_RANGE(NVTX_D2FRNL_KPT)
503 :
504 : ! If needed, manage ddk files
505 31607 : if (need_becfr.or.need_piezofr) then
506 3712 : do ii=1,3 ! Loop over elect. field directions
507 3712 : if (ddkfil(ii)/=0)then
508 : ! Number of k points to skip in the full set of k pointsp
509 2784 : ik_ddk(ii) = ddkfiles(ii)%findk(kpoint)
510 2784 : ABI_CHECK(ik_ddk(ii) /= -1, "Cannot find k-point in DDK")
511 2784 : npw_ = ddkfiles(ii)%hdr%npwarr(ik_ddk(ii))
512 2784 : if (npw_/=npw_k) then
513 : write(msg, '(a,i0,a,i0,a,i0,a,a,i0,a,a,i0)')&
514 0 : 'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',ii,ch10,&
515 0 : 'the number of plane waves in the ddk file is equal to', npw_,ch10,&
516 0 : 'while it should be ',npw_k
517 0 : ABI_ERROR(msg)
518 : end if
519 :
520 : end if
521 : end do
522 : end if
523 :
524 94821 : ABI_MALLOC(cwavef,(2,npw_k*dtset%nspinor*ndat))
525 31607 : if (need_becfr.or.need_piezofr) then
526 1856 : ABI_MALLOC(svectout,(2,npw_k*dtset%nspinor*ndat))
527 : end if
528 31607 : if (need_efmas) then
529 72 : ABI_MALLOC(cg_left,(2,npw_k*dtset%nspinor))
530 48 : ABI_MALLOC(gh2c,(2,npw_k*dtset%nspinor))
531 48 : ABI_MALLOC(gs2c,(2,npw_k*dtset%nspinor))
532 : end if
533 :
534 126428 : ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
535 31607 : if(rfstrs/=0.or.need_becfr.or.need_piezofr.or.need_efmas)then
536 26272 : ABI_MALLOC(ylmgr_k,(npw_k,9,mpsang*mpsang*psps%useylm))
537 : else
538 25039 : ABI_MALLOC(ylmgr_k,(0,0,0))
539 : end if
540 :
541 94821 : ABI_MALLOC(kg_k,(3,npw_k))
542 19037663 : kg_k(:,:) = 0
543 : !$OMP PARALLEL DO
544 4783121 : do ipw=1,npw_k
545 4751514 : kg_k(1,ipw)=kg(1,ipw+ikg)
546 4751514 : kg_k(2,ipw)=kg(2,ipw+ikg)
547 4783121 : kg_k(3,ipw)=kg(3,ipw+ikg)
548 : end do
549 31607 : if (psps%useylm==1) then
550 : !SOMP PARALLEL DO COLLAPSE(2)
551 57320 : do ilm=1,mpsang*mpsang
552 6681678 : do ipw=1,npw_k
553 6676040 : ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm)
554 : end do
555 : end do
556 5638 : if(rfstrs/=0.or.need_becfr.or.need_piezofr.or.need_efmas)then
557 : !SOMP PARALLEL DO COLLAPSE(3)
558 9410 : do ilm=1,mpsang*mpsang
559 85578 : do ii=1,9
560 8434829 : do ipw=1,npw_k
561 8426835 : ylmgr_k(ipw,ii,ilm)=ylmgr(ipw+ikg,ii,ilm)
562 : end do
563 : end do
564 : end do
565 : end if
566 : end if
567 :
568 31607 : cplex=2;if (istwf_k>1) cplex=1
569 :
570 : ! Compute (k+G) vectors (only if useylm=1)
571 31607 : nkpg=0
572 31607 : if (rfstrs/=0.or.need_efmas.or.pawpiezo==1) nkpg=3*dtset%nloalg(3)
573 126428 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
574 31607 : if (nkpg>0) then
575 586 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
576 : end if
577 :
578 : !EFMAS: Compute second order derivatives w/r to k for all direction for this k-point.
579 31607 : if (need_efmas) then
580 96 : ABI_MALLOC(ddkinpw,(npw_k,3,3))
581 96 : do mu=1,3
582 312 : do nu=1,3
583 : ! call d2kpg(ddkinpw(:,mu,nu),dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,mu,nu,kg_k,kpoint,npw_k)
584 288 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw(:,mu,nu),kpoint,npw_k,mu,nu)
585 : end do
586 : end do
587 : end if
588 :
589 : ! Compute nonlocal form factors ffnl at all (k+G):
590 31607 : ider=0;dimffnl=1;
591 31607 : if(need_becfr) then
592 928 : ider=1;dimffnl=4
593 : end if
594 31607 : if(rfstrs/=0.or.need_piezofr.or.need_efmas)then
595 6136 : ider=2;dimffnl=3+7*psps%useylm
596 : end if
597 189642 : ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
598 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
599 : & gmet,gprimd,ider,idir_ffnl,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
600 : & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,&
601 31607 : & psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
602 : #ifdef HAVE_OPENMP_OFFLOAD
603 : !$OMP TARGET ENTER DATA MAP(to:ffnl) IF(dtset%gpu_option==ABI_GPU_OPENMP)
604 : #endif
605 :
606 : ! For piezoelectric tensor need additional ffnl derivatives
607 31607 : if(need_piezofr)then
608 496 : ider_str=1 ; dimffnl_str=2
609 3472 : ABI_MALLOC(ffnl_str,(npw_k,dimffnl_str,psps%lmnmax,psps%ntypat,6))
610 3472 : do mu=1,6 !loop over strain
611 2976 : idir_str=-mu
612 : call mkffnl(psps%dimekb,dimffnl_str,psps%ekb,ffnl_str(:,:,:,:,mu),&
613 : & psps%ffspl,gmet,gprimd,ider_str,idir_str,psps%indlmn,kg_k,kpg_k,&
614 : & kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,&
615 3472 : & psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
616 : end do
617 : #ifdef HAVE_OPENMP_OFFLOAD
618 : !$OMP TARGET ENTER DATA MAP(to:ffnl_str) IF(dtset%gpu_option==ABI_GPU_OPENMP)
619 : #endif
620 : end if
621 :
622 : ! Load k-dependent part in the Hamiltonian datastructure
623 126428 : ABI_MALLOC(ph3d,(2,npw_k,gs_ham%matblk))
624 : call gs_ham%load_k(kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,&
625 31607 : & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_ph3d=.true.)
626 : #ifdef HAVE_OPENMP_OFFLOAD
627 : !$OMP TARGET ENTER DATA MAP(to:ph3d) IF(dtset%gpu_option==ABI_GPU_OPENMP)
628 : #endif
629 :
630 :
631 : ! Setup gemm_nonlop
632 31607 : if (gemm_nonlop_use_gemm) then
633 : call set_gemm_nonlop_ikpt(ikpt,gs_ham%npw_fft_k,gs_ham%istwf_k,gs_ham%indlmn,&
634 0 : & gs_ham%ntypat,gs_ham%nattyp,gs_ham%gpu_option)
635 : end if ! gemm_nonlop_use_gemm
636 :
637 :
638 :
639 : ! Initialize contributions from current k point
640 340335 : if(rfphon==1) dyfrnlk(:,:)=zero
641 31607 : if(rfstrs/=0)then
642 386576 : enlk=zero;eltfrnlk(:,:)=zero
643 : end if
644 57431 : if (need_becfr) becfrnlk(:,:,:)=zero
645 42023 : if (need_piezofr) piezofrnlk(:,:)=zero
646 31607 : if(need_efmas) then
647 : call check_degeneracies(efmasdeg(ikpt),dtset%efmas_bands(:,ikpt),nband_k,eigen(bdtot_index+1:bdtot_index+nband_k), &
648 24 : & dtset%efmas_deg_tol)
649 208 : do ideg=1,efmasdeg(ikpt)%ndegs
650 208 : if( efmasdeg(ikpt)%deg_range(1) <= ideg .and. ideg <= efmasdeg(ikpt)%deg_range(2) ) then
651 42 : deg_dim=efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
652 168 : ABI_MALLOC(efmasval(ideg,ikpt)%ch2c,(3,3,deg_dim,deg_dim))
653 126 : ABI_MALLOC(efmasval(ideg,ikpt)%eig2_diag,(3,3,deg_dim,deg_dim))
654 2734 : efmasval(ideg,ikpt)%ch2c=zero
655 2734 : efmasval(ideg,ikpt)%eig2_diag=zero
656 : else
657 142 : ABI_MALLOC(efmasval(ideg,ikpt)%ch2c,(0,0,0,0))
658 142 : ABI_MALLOC(efmasval(ideg,ikpt)%eig2_diag,(0,0,0,0))
659 : end if
660 : end do
661 : end if
662 :
663 : ! Loop over bands
664 31607 : iband_me = 0
665 210299 : do iband=1,nband_k,ndat
666 :
667 178692 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me)) cycle
668 178548 : iband_me = iband_me + ndat
669 :
670 357096 : occ_k(1:ndat)=occ(iband+bdtot_index:iband+bdtot_index+ndat-1)
671 : cwavef(:,1:npw_k*dtset%nspinor*ndat) = cg(:,1+(iband_me-ndat)*npw_k*dtset%nspinor+icg: &
672 121593021 : & iband_me *npw_k*dtset%nspinor+icg)
673 :
674 : ! Compute non-local contributions from n,k
675 193012 : if (psps%usepaw==1) eig_k(1:ndat)=eigen(iband+bdtot_index:iband+bdtot_index+ndat-1)
676 :
677 : ! === Dynamical matrix
678 178548 : if(rfphon==1) then
679 : call nonlop(choice_phon,cpopt,cwaveprj,enlout_phon,gs_ham,idir,eig_k,mpi_enreg,ndat,&
680 161570 : & nnlout_phon,paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
681 : ! Accumulate non-local contributions from n,k
682 323140 : do idat=1,ndat
683 : dyfrnlk(:,:) = dyfrnlk(:,:) &
684 2474862 : & + occ_k(idat) * reshape(enlout_phon(1+(idat-1)*nnlout_phon:idat*nnlout_phon),(/6,natom/))
685 : end do
686 : end if
687 :
688 : ! === Elastic tensor
689 178548 : if(rfstrs/=0) then
690 : call nonlop(choice_strs,cpopt,cwaveprj,enlout_strs,gs_ham,idir,(/eig_k/),mpi_enreg,ndat,&
691 70296 : & nnlout_strs,paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
692 : ! Accumulate non-local contribut ions from n,k
693 46864 : do idat=1,ndat
694 : eltfrnlk(:,:) = eltfrnlk(:,:) &
695 1699088 : & + occ_k(idat) * reshape(enlout_strs(1+(idat-1)*nnlout_strs:idat*nnlout_strs),(/3*natom+6,6/))
696 : end do
697 : end if !endo if strs
698 :
699 : ! PAW: accumulate gradients of rhoij
700 : !EFMAS: Bug with efmas currently; to be looked into...
701 178548 : if (psps%usepaw==1.and.(.not.need_efmas)) then
702 28576 : do idat=1,ndat
703 : call pawaccrhoij(gs_ham%atindx,cplex,cwaveprj(:,1+(idat-1)*dtset%nspinor:idat*dtset%nspinor),cwaveprj(:,1+(idat-1)*dtset%nspinor:idat*dtset%nspinor),0,isppol,natom,&
704 28576 : & natom,dtset%nspinor,occ_k(idat),3,pawrhoij_tot,use_timerev,use_zeromag,wtk_k)
705 : end do
706 : end if
707 :
708 : ! PAW: Compute frozen contribution to piezo electric tensor
709 178548 : if (need_piezofr) then
710 10304 : do ii=1,3 ! Loop over elect. field directions
711 : call nonlop(choice_piez3,cpopt,cwaveprj,enlout_piez1,gs_ham,0,(/zero/),mpi_enreg,ndat,&
712 7728 : & nnlout_piez1,paw_opt_1,signs,nonlop_dum,tim_nonlop,cwavef,cwavef,enl=becij(:,:,:,:,ii))
713 18032 : do idat=1,ndat
714 : piezofrnlk(:,ii) = piezofrnlk(:,ii) &
715 61824 : & + occ_k(idat) * enlout_piez1(1+(idat-1)*nnlout_piez1:idat*nnlout_piez1)
716 : end do
717 : end do !end do ii
718 : end if
719 :
720 : ! PAW: Compute frozen contribution to Born Effective Charges
721 178548 : if (need_becfr) then
722 17920 : do ii=1,3 ! Loop over elect. field directions
723 : call nonlop(choice_bec2,cpopt,cwaveprj,enlout_bec1,gs_ham,0,(/zero/),mpi_enreg,ndat,&
724 13440 : & nnlout_bec1,paw_opt_1,signs,nonlop_dum,tim_nonlop,cwavef,cwavef,enl=becij(:,:,:,:,ii))
725 31360 : do idat=1,ndat
726 : becfrnlk(:,:,ii) = becfrnlk(:,:,ii) &
727 168960 : & + occ_k(idat) * reshape(enlout_bec1(1+(idat-1)*nnlout_bec1:idat*nnlout_bec1),(/3,natom/))
728 : end do
729 : end do !end do ii
730 : end if
731 :
732 178548 : if (need_becfr.or.need_piezofr) then
733 17920 : do ii=1,3 ! Loop over elect. field directions
734 : ! Not able to compute if ipert=(Elect. field) and no ddk WF file
735 13440 : if (ddkfil(ii)==0) cycle
736 : ! Read ddk wave function
737 53760 : ABI_MALLOC(ddk,(2,npw_k*dtset%nspinor,ndat))
738 : if (ddkfil(ii)/=0) then
739 26880 : do idat=1,ndat
740 13440 : call ddkfiles(ii)%read_bks(iband+(idat-1), ik_ddk(ii), isppol, xmpio_single, cg_bks=ddk(:,:,idat))
741 : ! Multiply ddk by +i
742 1739136 : do jj=1,npw_k*dtset%nspinor
743 1712256 : arg=ddk(1,jj,idat)
744 1725696 : ddk(1,jj,idat)=-ddk(2,jj,idat);ddk(2,jj,idat)=arg
745 : end do
746 : end do
747 : else
748 : ddk=zero
749 : end if
750 :
751 13440 : if(need_becfr)then
752 53760 : do mu=1,3 !loop over atom direction
753 140160 : do iatom=1,natom !Loop over atom
754 86400 : ia=gs_ham%atindx(iatom)
755 : call nonlop(choice_bec2,cpopt_bec,cwaveprj,enlout_bec1,gs_ham,mu,(/zero/),&
756 : & mpi_enreg,ndat,nnlout_bec1,paw_opt_3,signs_field,svectout,tim_nonlop,&
757 86400 : & cwavef,cwavef,iatom_only=iatom)
758 213120 : do idat=1,ndat
759 : call dotprod_g(dotprod(1),dotprod(2),istwf_k,npw_k*dtset%nspinor,2,&
760 : & svectout(:,1+(idat-1)*npw_k*dtset%nspinor:idat*npw_k*dtset%nspinor),ddk(:,:,idat),&
761 86400 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
762 172800 : becfrnlk(mu,ia,ii)=becfrnlk(mu,ia,ii)+occ_k(idat)*dotprod(1)
763 : end do
764 : end do
765 : end do
766 : end if
767 :
768 13440 : if(need_piezofr)then
769 54096 : do mu=1,6 !loop over strain
770 46368 : call gs_ham%load_k(ffnl_k=ffnl_str(:,:,:,:,mu))
771 : call nonlop(choice_piez3,cpopt,cwaveprj,enlout_piez1,gs_ham,mu,(/zero/),mpi_enreg,ndat,&
772 46368 : & nnlout_piez1,paw_opt_3,signs_field,svectout,tim_nonlop,cwavef,svectout)
773 100464 : do idat=1,ndat
774 : call dotprod_g(dotprod(1),dotprod(2),istwf_k,npw_k*dtset%nspinor,2,&
775 : & svectout(:,1+(idat-1)*npw_k*dtset%nspinor:idat*npw_k*dtset%nspinor),ddk(:,:,idat),&
776 46368 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
777 92736 : piezofrnlk(mu,ii)=piezofrnlk(mu,ii)+occ_k(idat)*dotprod(1)
778 : end do
779 : end do
780 7728 : call gs_ham%load_k(ffnl_k=ffnl)
781 : end if
782 :
783 17920 : ABI_FREE(ddk)
784 : end do ! End loop ddk file
785 : end if
786 :
787 4480 : if(need_piezofr)then
788 95312 : enlout_piez2 = zero
789 : call nonlop(choice_piez55,cpopt,cwaveprj,enlout_piez2,gs_ham,0,(/zero/),mpi_enreg,ndat,&
790 2576 : & nnlout_piez2,paw_opt_3,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
791 : ! Multiply enlout by +i
792 2576 : iashift = 1
793 5152 : do idat=1,ndat
794 20608 : do mu=1,6 ! strain
795 64400 : do nu=1,3 ! k
796 46368 : piezofrnlk(mu,nu)=piezofrnlk(mu,nu)-occ_k(idat)*(enlout_piez2(iashift+1)) ! Real part
797 : ! piezofrnlk(mu,nu)=piezofrnlk(mu,nu)+occ_k(idat)*(enlout_piez2(iashift ))! Imaginary part
798 61824 : iashift = iashift + 2
799 : end do
800 : end do
801 : end do
802 : end if
803 :
804 178548 : if(need_becfr)then
805 : call nonlop(choice_bec54,cpopt,cwaveprj,enlout_bec2,gs_ham,0,(/zero/),mpi_enreg,ndat,&
806 4480 : & nnlout_bec2,paw_opt_3,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
807 : ! Multiply enlout by +i
808 4480 : iashift = 1
809 8960 : do idat=1,ndat
810 18560 : do iatom=1,natom ! atm
811 42880 : do mu=1,3 ! atm pos.
812 124800 : do nu=1,3 ! k
813 86400 : becfrnlk(mu,iatom,nu)=becfrnlk(mu,iatom,nu)-occ_k(idat)*(enlout_bec2(iashift+1)) ! Real part
814 : ! becfrnlk(mu,iatom,nu)=becfrnlk(mu,iatom,nu)+occ_k(idat)*(enlout_bec2(iashift ))! Imaginary part
815 115200 : iashift = iashift + 2
816 : end do
817 : end do
818 : end do
819 : end do
820 : end if
821 :
822 210155 : if(need_efmas) then
823 338 : bandmin=efmasdeg(ikpt)%degs_bounds(1, efmasdeg(ikpt)%deg_range(1) )
824 338 : bandmax=efmasdeg(ikpt)%degs_bounds(2, efmasdeg(ikpt)%deg_range(2) )
825 :
826 338 : choice_efmas=8; signs=2
827 338 : cpopt=-1 !To prevent re-use of stored dgxdt, which are not for all direction required for EFMAS.
828 338 : paw_opt_efmas=0; if(psps%usepaw/=0) paw_opt_efmas=4 !To get both gh2c and gs2c
829 338 : nnlout_efmas=0; tim_nonlop=0 ! No tim_nonlop for efmas, currently.
830 :
831 676 : do idat=1,ndat
832 : ! find list of iband which are running now:
833 5604 : bands_treated_now = 0
834 338 : bands_treated_now(iband+idat-1) = 1
835 338 : call xmpi_sum(bands_treated_now,mpi_enreg%comm_band,ierr)
836 :
837 : ! for all iband running right now
838 1905 : do iband_ = bandmin, bandmax
839 1229 : if (bands_treated_now(iband_) == 0) cycle
840 :
841 654 : do mu=1,3
842 2177 : do nu=1,3
843 : ! if I have iband_ prepare things
844 711 : if (iband_ == iband+idat-1) then
845 711 : idir=3*(mu-1)+nu !xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9, (xyz,xyz)=(mu,nu)
846 2450529 : gh2c=zero; gs2c=zero
847 : call nonlop(choice_efmas,cpopt,cwaveprj,enlout_efmas,gs_ham,idir,(/eig_k/),mpi_enreg,&
848 2133 : 1,nnlout_efmas,paw_opt_efmas,signs,gs2c,tim_nonlop,cwavef(:,1+(idat-1)*npw_k*dtset%nspinor:idat*npw_k*dtset%nspinor),gh2c)
849 : !DEBUG
850 : ! gh2c=zero; gs2c=zero
851 : !ENDDEBUG
852 1692 : do ispinor=1,dtset%nspinor
853 981 : ii = 1+(ispinor-1)*npw_k
854 3654 : do icplx=1,2
855 : gh2c(icplx,ii:ispinor*npw_k) = gh2c(icplx,ii:ispinor*npw_k) + &
856 819549 : & ddkinpw(1:npw_k,mu,nu)*cwavef(icplx,ii+(idat-1)*npw_k:ispinor*npw_k+(idat-1)*npw_k)
857 : end do
858 : end do
859 1226331 : gh2c = gh2c - eig_k(idat)*gs2c
860 : !DEBUG
861 : ! gh2c=zero; gs2c=zero
862 : !ENDDEBUG
863 : end if
864 711 : ideg = efmasdeg(ikpt)%ideg(iband+idat-1)
865 2133 : ABI_MALLOC( ch2c_tmp, (size(efmasval(ideg,ikpt)%ch2c, dim=3)) )
866 :
867 : ! share gh2c
868 711 : call xmpi_bcast(gh2c, band_procs(iband+idat-1), mpi_enreg%comm_band,ierr)
869 :
870 711 : jband_me = 0
871 5661 : do jband=1,efmasdeg(ikpt)%degs_bounds(2,ideg)
872 : ! jband treated on current proc?
873 4950 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,jband,jband,isppol,me)) cycle
874 : ! if so, indexing of the bands in my cg array
875 4950 : jband_me = jband_me + 1
876 : ! if we do not need to treat it for efmas, skip
877 4950 : if (jband < efmasdeg(ikpt)%degs_bounds(1,ideg)) cycle
878 :
879 2881710 : cg_left(:,1:npw_k*dtset%nspinor) = cg(:,1+(jband_me-1)*npw_k*dtset%nspinor+icg:jband_me*npw_k*dtset%nspinor+icg)
880 1809 : dotprod=0
881 : call dotprod_g(dotprod(1),dotprod(2),istwf_k,npw_k*dtset%nspinor,2,cg_left,gh2c,mpi_enreg%me_g0,&
882 1809 : mpi_enreg%comm_spinorfft)
883 1809 : isub = iband-efmasdeg(ikpt)%degs_bounds(1,ideg)+1
884 1809 : jsub = jband-efmasdeg(ikpt)%degs_bounds(1,ideg)+1
885 :
886 5661 : ch2c_tmp(jsub)=cmplx(dotprod(1),dotprod(2),kind=dp)
887 : end do ! jband
888 : !mpi_sum ch2c_tmp to get all jband contribs
889 711 : call xmpi_sum(ch2c_tmp,mpi_enreg%comm_band,ierr)
890 2520 : efmasval(ideg,ikpt)%ch2c(mu,nu,:,isub)=ch2c_tmp(:)
891 2370 : ABI_FREE( ch2c_tmp )
892 : end do ! nu
893 : end do ! mu
894 : end do ! iband_
895 : end do ! idat
896 : end if ! need_efmas
897 :
898 : end do ! End of loop on bands
899 :
900 31607 : if(rfphon==1) then
901 71631 : do iatom=1,natom
902 44104 : ia=iatom;if (dyfr_nondiag==0) ia=1
903 44104 : dyfrnl(1,1,1,iatom,ia)=dyfrnl(1,1,1,iatom,ia)+wtk_k*dyfrnlk(1,iatom)
904 44104 : dyfrnl(1,2,2,iatom,ia)=dyfrnl(1,2,2,iatom,ia)+wtk_k*dyfrnlk(2,iatom)
905 44104 : dyfrnl(1,3,3,iatom,ia)=dyfrnl(1,3,3,iatom,ia)+wtk_k*dyfrnlk(3,iatom)
906 44104 : dyfrnl(1,2,3,iatom,ia)=dyfrnl(1,2,3,iatom,ia)+wtk_k*dyfrnlk(4,iatom)
907 44104 : dyfrnl(1,1,3,iatom,ia)=dyfrnl(1,1,3,iatom,ia)+wtk_k*dyfrnlk(5,iatom)
908 71631 : dyfrnl(1,1,2,iatom,ia)=dyfrnl(1,1,2,iatom,ia)+wtk_k*dyfrnlk(6,iatom)
909 : end do
910 : end if ! end if rfphon
911 31607 : if(rfstrs/=0) then
912 386576 : eltfrnl(:,:)=eltfrnl(:,:)+dtset%wtk(ikpt)*eltfrnlk(:,:)
913 : end if
914 31607 : if(need_becfr) then
915 26752 : becfrnl(:,:,:)=becfrnl(:,:,:)+dtset%wtk(ikpt)*becfrnlk(:,:,:)
916 : end if
917 31607 : if(need_piezofr) then
918 10912 : piezofrnl(:,:)=piezofrnl(:,:)+dtset%wtk(ikpt)*piezofrnlk(:,:)
919 : end if
920 : ! Increment indexes
921 31607 : bdtot_index=bdtot_index+nband_k
922 31607 : if (mkmem_rbz/=0) then
923 31607 : ibg=ibg+nband_k*dtset%nspinor
924 31607 : icg=icg+npw_k*dtset%nspinor*nband_me
925 31607 : ikg=ikg+npw_k
926 : end if
927 :
928 : #ifdef HAVE_OPENMP_OFFLOAD
929 : !$OMP TARGET EXIT DATA MAP(delete:ffnl,ph3d) IF(dtset%gpu_option==ABI_GPU_OPENMP)
930 : #endif
931 31607 : ABI_FREE(ffnl)
932 31607 : ABI_FREE(kpg_k)
933 31607 : ABI_FREE(ph3d)
934 31607 : ABI_FREE(ylm_k)
935 31607 : ABI_FREE(ylmgr_k)
936 31607 : ABI_FREE(cwavef)
937 31607 : ABI_FREE(kg_k)
938 31607 : if (need_becfr.or.need_piezofr) then
939 928 : ABI_FREE(svectout)
940 : end if
941 928 : if (need_piezofr) then
942 : #ifdef HAVE_OPENMP_OFFLOAD
943 : !$OMP TARGET EXIT DATA MAP(delete:ffnl_str) IF(dtset%gpu_option==ABI_GPU_OPENMP)
944 : #endif
945 496 : ABI_FREE(ffnl_str)
946 : end if
947 32129 : if (need_efmas) then
948 24 : ABI_FREE(ddkinpw)
949 24 : ABI_FREE(cg_left)
950 24 : ABI_FREE(gh2c)
951 24 : ABI_FREE(gs2c)
952 : end if
953 :
954 : ABI_NVTX_END_RANGE()
955 : end do ! End loops on isppol and ikpt
956 : end do
957 487 : if(rfphon==1) then
958 450 : ABI_FREE(dyfrnlk)
959 450 : ABI_FREE(enlout_phon)
960 : end if
961 487 : if(rfstrs/=0) then
962 37 : ABI_FREE(eltfrnlk)
963 37 : ABI_FREE(enlout_strs)
964 : end if
965 487 : if (need_becfr) then
966 27 : ABI_FREE(becfrnlk)
967 27 : ABI_FREE(enlout_bec1)
968 27 : ABI_FREE(enlout_bec2)
969 : end if
970 487 : if(need_piezofr)then
971 18 : ABI_FREE(enlout_piez1)
972 18 : ABI_FREE(enlout_piez2)
973 18 : ABI_FREE(piezofrnlk)
974 : end if
975 487 : if(need_efmas) then
976 17 : ABI_FREE(enlout_efmas)
977 : end if
978 487 : if (psps%usepaw==1) then
979 57 : if (need_becfr.or.need_piezofr) then
980 27 : ABI_FREE(becij)
981 : end if
982 57 : call pawcprj_free(cwaveprj)
983 : end if
984 653 : ABI_FREE(cwaveprj)
985 :
986 : !Fill in lower triangle of matrixes
987 487 : if (rfphon==1) then
988 1388 : do iatom=1,natom
989 938 : ia=iatom;if (dyfr_nondiag==0) ia=1
990 938 : dyfrnl(1,3,2,iatom,ia)=dyfrnl(1,2,3,iatom,ia)
991 938 : dyfrnl(1,3,1,iatom,ia)=dyfrnl(1,1,3,iatom,ia)
992 1388 : dyfrnl(1,2,1,iatom,ia)=dyfrnl(1,1,2,iatom,ia)
993 : end do
994 : end if
995 487 : if(rfstrs/=0)then
996 222 : do jj=2,6
997 777 : do ii=1,jj-1
998 740 : eltfrnl(jj,ii)=eltfrnl(ii,jj)
999 : end do
1000 : end do
1001 : end if
1002 :
1003 : !Parallel case: accumulate (n,k) contributions
1004 : if (xmpi_paral==1) then
1005 487 : call timab(48,1,tsec)
1006 : ! Accumulate dyfrnl
1007 487 : if(rfphon==1)then
1008 450 : call xmpi_sum(dyfrnl,spaceworld,ierr)
1009 : end if
1010 : ! Accumulate eltfrnl.
1011 487 : if(rfstrs/=0)then
1012 37 : call xmpi_sum(eltfrnl,spaceworld,ierr)
1013 : end if
1014 : ! Accumulate becfrnl
1015 487 : if (need_becfr) then
1016 27 : call xmpi_sum(becfrnl,spaceworld,ierr)
1017 : end if
1018 : ! Accumulate piezofrnl
1019 487 : if (need_piezofr) then
1020 18 : call xmpi_sum(piezofrnl,spaceworld,ierr)
1021 : end if
1022 :
1023 : ! PAW: accumulate gradients of rhoij
1024 487 : if (psps%usepaw==1) then
1025 171 : ABI_MALLOC(dimlmn,(natom))
1026 215 : dimlmn(1:natom)=pawrhoij_tot(1:natom)%cplex_rhoij*pawrhoij_tot(1:natom)%qphase*pawrhoij_tot(1:natom)%lmn2_size
1027 215 : bufdim=ncpgr*sum(dimlmn)*nsploop
1028 171 : ABI_MALLOC(mpibuf,(bufdim))
1029 49167 : ii=0;mpibuf=zero
1030 215 : do iatom=1,natom
1031 379 : do isppol=1,nsploop
1032 1432 : do mu=1,ncpgr
1033 50220 : mpibuf(ii+1:ii+dimlmn(iatom))=pawrhoij_tot(iatom)%grhoij(mu,1:dimlmn(iatom),isppol)
1034 1274 : ii=ii+dimlmn(iatom)
1035 : end do
1036 : end do
1037 : end do
1038 57 : call xmpi_sum(mpibuf,spaceworld,ierr)
1039 57 : ii=0
1040 215 : do iatom=1,natom
1041 379 : do isppol=1,nsploop
1042 1432 : do mu=1,ncpgr
1043 50220 : pawrhoij_tot(iatom)%grhoij(mu,1:dimlmn(iatom),isppol)=mpibuf(ii+1:ii+dimlmn(iatom))
1044 1274 : ii=ii+dimlmn(iatom)
1045 : end do
1046 : end do
1047 : end do
1048 57 : ABI_FREE(mpibuf)
1049 114 : ABI_FREE(dimlmn)
1050 : end if
1051 487 : call timab(48,2,tsec)
1052 : end if
1053 :
1054 : !====== PAW: Additional steps
1055 487 : if (psps%usepaw==1) then
1056 :
1057 : ! Symmetrize rhoij gradients and transfer to cartesian (reciprocal space) coord.
1058 : ! This symetrization is necessary in the antiferromagnetic case...
1059 57 : if (rfphon==1.and.rfstrs==0) then
1060 32 : option_rhoij=2;option=0
1061 : call pawrhoij_symrhoij(pawrhoij_tot,pawrhoij_tot,option_rhoij,gprimd,indsym,0,natom,dtset%nsym,&
1062 : & psps%ntypat,option,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,dtset%typat,&
1063 32 : & comm_atom=my_comm_atom, mpi_atmtab=my_atmtab)
1064 25 : else if (rfphon==1.and.rfstrs==1) then
1065 0 : option_rhoij=23;option=0
1066 : call pawrhoij_symrhoij(pawrhoij_tot,pawrhoij_tot,option_rhoij,gprimd,indsym,0,natom,dtset%nsym,&
1067 : & psps%ntypat,option,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,dtset%typat,&
1068 0 : & comm_atom=my_comm_atom, mpi_atmtab=my_atmtab)
1069 : end if
1070 :
1071 : ! Translate coordinates
1072 57 : ABI_CHECK(nsploop/=4,'d2frnl: should we mix mx/my/mz when translating coordinates?')
1073 215 : do iatom=1,natom
1074 158 : cplx=pawrhoij_tot(iatom)%cplex_rhoij
1075 373 : do iq=1,pawrhoij_tot(iatom)%qphase
1076 158 : iq0=0;if (iq==2) iq0=cplx*pawrhoij_tot(iatom)%lmn2_size
1077 480 : do isppol=1,nsploop
1078 8258 : do klmn=1,pawrhoij_tot(iatom)%lmn2_size
1079 16764 : do ii=1,cplx
1080 8664 : if(rfphon==1.or.rfstrs/=0)then
1081 24136 : grhoij(1:3)=pawrhoij_tot(iatom)%grhoij(shift_rhoij+1:shift_rhoij+3,iq0+cplx*(klmn-1)+ii,isppol)
1082 24136 : do mu=1,3
1083 : pawrhoij_tot(iatom)%grhoij(shift_rhoij+mu,iq0+cplx*(klmn-1)+ii,isppol)=gprimd(mu,1)*grhoij(1)&
1084 26766 : & +gprimd(mu,2)*grhoij(2)+gprimd(mu,3)*grhoij(3)
1085 : end do
1086 : end if
1087 16600 : if(rfstrs/=0)then
1088 : call strconv(pawrhoij_tot(iatom)%grhoij(1:6,iq0+cplx*(klmn-1)+ii,isppol),gprimd,&
1089 1512 : & pawrhoij_tot(iatom)%grhoij(1:6,iq0+cplx*(klmn-1)+ii,isppol))
1090 : end if
1091 : end do
1092 : end do
1093 : end do
1094 : end do
1095 : end do
1096 :
1097 : ! In case of elastic tensor computation, add diagonal contribution:
1098 : ! -delta_{alphabeta} rhoi_{ij} to drhoij/d_eps
1099 57 : if(rfstrs/=0)then
1100 54 : do iatom=1,natom
1101 42 : cplx=pawrhoij_tot(iatom)%cplex_rhoij
1102 96 : do iq=1,pawrhoij_tot(iatom)%qphase
1103 42 : iq0=0;if (iq==2) iq0=cplx*pawrhoij_tot(iatom)%lmn2_size
1104 126 : do isppol=1,nsploop
1105 748 : do nu=1,pawrhoij_tot(iatom)%nrhoijsel
1106 664 : klmn=pawrhoij_tot(iatom)%rhoijselect(nu)
1107 1370 : do ii=1,cplx
1108 : pawrhoij_tot(iatom)%grhoij(1:3,iq0+cplx*(klmn-1)+ii,isppol)= &
1109 : & pawrhoij_tot(iatom)%grhoij(1:3,iq0+cplx*(klmn-1)+ii,isppol)&
1110 3320 : & -pawrhoij_tot(iatom)%rhoijp(iq0+cplx*(nu-1)+ii,isppol)
1111 : end do
1112 : end do
1113 : end do
1114 : end do
1115 : end do
1116 : end if
1117 :
1118 : ! Add gradients due to Dij derivatives to dynamical matrix/stress tensor
1119 57 : dimnhat=0;optgr=0;optgr2=0;optstr=0;optstr2=0
1120 57 : if (rfphon==1) optgr2=1
1121 57 : if (rfstrs/=0) optstr2=1
1122 57 : ABI_MALLOC(nhat_dum,(1,0))
1123 : call pawgrnl(gs_ham%atindx1,dimnhat,dyfrnl,dyfr_cplex,eltfrnl,dummy,gsqcut,mgfftf,my_natom,natom,&
1124 : & gs_ham%nattyp,nfftf,ngfftf,nhat_dum,dummy,dtset%nspden,dtset%nsym,psps%ntypat,optgr,optgr2,optstr,optstr2,&
1125 : & pawang,pawfgrtab,pawrhoij_tot,pawtab,ph1df,psps,dtset%qptn,rprimd,symrec,dtset%typat,ucvol,vtrial,vxc,xred,&
1126 57 : & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
1127 57 : ABI_FREE(nhat_dum)
1128 : end if !PAW
1129 :
1130 : !The indexing array atindx is used to reestablish the correct order of atoms
1131 487 : if (rfstrs/=0)then
1132 148 : ABI_MALLOC(elt_work,(6+3*natom,6))
1133 1591 : elt_work(1:6,1:6)=eltfrnl(1:6,1:6)
1134 131 : do ia=1,natom
1135 94 : ielt=7+3*(ia-1)
1136 94 : ieltx=7+3*(gs_ham%atindx(ia)-1)
1137 2387 : elt_work(ielt:ielt+2,1:6)=eltfrnl(ieltx:ieltx+2,1:6)
1138 : end do
1139 3283 : eltfrnl(:,:)=elt_work(:,:)
1140 37 : ABI_FREE(elt_work)
1141 : end if
1142 :
1143 : !Born Effective Charges and PAW:
1144 : !1-Re-order atoms -- 2-Add diagonal contribution from rhoij
1145 : !3-Multiply by -1 because that the effective charges
1146 : ! are minus the second derivatives of the energy
1147 487 : if (need_becfr) then
1148 108 : ABI_MALLOC(becfrnl_tmp,(3,natom,3))
1149 975 : becfrnl_tmp=-becfrnl
1150 97 : do ia=1,natom ! Atom (sorted by type)
1151 70 : iatom=gs_ham%atindx1(ia) ! Atom (not sorted)
1152 70 : itypat=dtset%typat(iatom)
1153 307 : do ii=1,3 ! Direction of electric field
1154 910 : do jj=1,3 ! Direction of atom
1155 840 : becfrnl(jj,iatom,ii)=becfrnl_tmp(jj,ia,ii)
1156 : end do
1157 : end do
1158 : end do
1159 27 : ABI_FREE(becfrnl_tmp)
1160 : end if
1161 :
1162 : !Piezoelectric Tensor
1163 : !-Multiply by -1 because that the piezoelectric tensor
1164 : ! are minus the second derivatives of the energy
1165 487 : if (need_piezofr) then
1166 396 : piezofrnl=-piezofrnl
1167 : end if
1168 :
1169 487 : ABI_FREE(eig_k)
1170 487 : ABI_FREE(occ_k)
1171 :
1172 : !Close the ddk files
1173 1948 : do ii=1,3
1174 1948 : call ddkfiles(ii)%close()
1175 : end do
1176 :
1177 : !Release now useless memory
1178 487 : if (psps%usepaw==1) then
1179 215 : do iatom=1,natom
1180 158 : ABI_FREE(pawrhoij_tot(iatom)%grhoij)
1181 215 : pawrhoij_tot(iatom)%ngrhoij=0
1182 : end do
1183 57 : if (paral_atom) then
1184 26 : call pawrhoij_free(pawrhoij_tot)
1185 26 : ABI_FREE(pawrhoij_tot)
1186 : end if
1187 : end if
1188 487 : call gs_ham%free()
1189 :
1190 : ABI_NVTX_END_RANGE()
1191 487 : call timab(159,2,tsec)
1192 :
1193 487 : write(msg,'(3a)')ch10,' ==> Calculation of the frozen part of the second order derivative done',ch10
1194 487 : call wrtout(std_out,msg,'COLL')
1195 :
1196 : DBG_EXIT("COLL")
1197 :
1198 1461 : end subroutine d2frnl
1199 : !!***
1200 :
1201 : end module m_d2frnl
1202 : !!***
|
