Branch data 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 : 488 : subroutine d2frnl(becfrnl,cg,dtfil,dtset,dyfrnl,dyfr_cplex,dyfr_nondiag,efmasdeg,efmasval,eigen,eltfrnl,&
154 : 488 : & gsqcut,has_allddk,indsym,kg,mband_mem_rbz,mkmem_rbz,mgfftf,mpi_enreg,mpsang,my_natom,natom,nfftf,ngfft,ngfftf,npwarr,&
155 : 488 : & occ,paw_ij,pawang,pawbec,pawfgrtab,pawpiezo,pawrad,pawrhoij,pawtab,ph1d,ph1df,piezofrnl,psps,&
156 : 488 : & 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 : 488 : type(gs_hamiltonian_type) :: gs_ham
213 : : !arrays
214 : : integer :: ik_ddk(3),ddkfil(3)
215 : 976 : integer :: bands_treated_now(dtset%mband), band_procs(dtset%mband)
216 : : integer :: ndat,idat
217 : 488 : integer,allocatable :: dimlmn(:),kg_k(:,:),l_size_atm(:)
218 : 488 : 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 : 488 : complex(dp), allocatable :: ch2c_tmp(:)
222 : 488 : real(dp),allocatable :: becfrnl_tmp(:,:,:),becfrnlk(:,:,:),becij(:,:,:,:,:),cg_left(:,:)
223 : 488 : real(dp),allocatable :: occ_k(:),eig_k(:)
224 : 976 : real(dp),allocatable :: cwavef(:,:),ddk(:,:,:),ddkinpw(:,:,:),dyfrnlk(:,:)
225 : 488 : real(dp),allocatable :: elt_work(:,:),eltfrnlk(:,:),enlout_bec1(:),enlout_bec2(:),enlout_efmas(:)
226 : 488 : real(dp),allocatable :: enlout_piez1(:),enlout_piez2(:),enlout_phon(:),enlout_strs(:)
227 : 488 : real(dp),allocatable :: gh2c(:,:),gs2c(:,:)
228 : 488 : real(dp),allocatable :: kpg_k(:,:),mpibuf(:),nhat_dum(:,:),piezofrnlk(:,:),ph3d(:,:,:)
229 : 488 : real(dp),allocatable :: svectout(:,:),ylm_k(:,:),ylmgr_k(:,:,:)
230 : 488 : real(dp),allocatable,target :: ffnl(:,:,:,:),ffnl_str(:,:,:,:,:)
231 : : character(len=fnlen) :: fiwfddk(3)
232 : 488 : type(paw_ij_type),allocatable :: paw_ij_tmp(:)
233 : 488 : type(pawcprj_type),allocatable,target :: cwaveprj(:,:)
234 : 488 : type(pawfgrtab_type),allocatable :: pawfgrtab_tmp(:)
235 : 488 : type(pawrhoij_type),pointer :: pawrhoij_tot(:)
236 [ + + + + : 47824 : type(wfk_t) :: ddkfiles(3)
+ + + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - - + -
+ - + - ]
237 : : ! *************************************************************************
238 : :
239 : : DBG_ENTER("COLL")
240 : :
241 : 488 : call timab(159,1,tsec)
242 : : ABI_NVTX_START_RANGE(NVTX_D2FRNL)
243 : :
244 : 488 : write(msg,'(3a)')ch10,' ==> Calculation of the frozen part of the second order derivatives, this can take some time...',ch10
245 : 488 : call wrtout(std_out,msg,'COLL')
246 : :
247 : : !Set up parallelism
248 : 488 : spaceworld=mpi_enreg%comm_cell
249 : 488 : me=mpi_enreg%me_kpt
250 : 488 : master=0
251 : 488 : paral_atom=(my_natom/=natom)
252 : 488 : my_comm_atom=mpi_enreg%comm_atom
253 : 488 : my_atmtab=>mpi_enreg%my_atmtab
254 : 488 : ndat = dtset%bandpp
255 : :
256 : : !Compute gmet, gprimd and ucvol from rprimd
257 : 488 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
258 : :
259 : : !If needed, check for ddk files (used for effective charges)
260 [ + + - + ]: 488 : if (pawbec==1.or.pawpiezo==1) then
261 : 52 : ddkfil(:)=0
262 [ + + ]: 208 : do ii=1,3
263 : 156 : ddkcase=ii+natom*3
264 : 156 : call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk(ii))
265 : 156 : t_test = file_exists(fiwfddk(ii))
266 : : ! Trick needed to run Abinit test suite in netcdf mode.
267 [ - + + + ]: 156 : 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 : 237 : call wrtout(std_out,msg,'COLL')
271 : : end if
272 : 133 : if (t_test) ddkfil(ii)=20+ii ! Note the use of unit numbers 21, 22 and 23
273 : : end do
274 [ + + + + ]: 127 : has_ddk_file=(any(ddkfil(:)>0))
275 [ + + + + ]: 133 : has_allddk =(all(ddkfil(:)>0))
276 : : else
277 : 436 : has_ddk_file=.FALSE.
278 : 436 : has_allddk =.FALSE.
279 : : end if
280 : :
281 [ + + - + : 488 : if(pawbec==1.or.pawpiezo==1.and.has_ddk_file) then
- - ]
282 [ + + ]: 52 : if(.not.has_allddk) then
283 : 25 : write(msg,'(5a)')ch10,&
284 : 25 : & ' WARNING: All ddk perturbations are needed to compute',ch10,&
285 : 50 : & ' the frozen part of Born effective charges and/or piezoelectric tensor.',ch10
286 : 25 : 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 [ + + + + ]: 488 : need_becfr=(pawbec==1.and.has_ddk_file)
296 [ + + + + ]: 488 : need_piezofr=(pawpiezo==1.and.has_ddk_file)
297 : :
298 : : !Initialization of frozen non local array
299 [ + + ]: 488 : if(rfphon==1) then
300 [ + + + + : 32659 : dyfrnl(:,:,:,:,:)=zero
+ + + + +
+ ]
301 [ - + + - : 1353 : ABI_MALLOC(dyfrnlk,(6,natom))
+ - + - -
+ - + ]
302 : : end if
303 [ + + ]: 488 : if(rfstrs/=0)then
304 [ + + + + ]: 3283 : eltfrnl(:,:)=zero;enl=zero
305 [ + - - + : 185 : ABI_MALLOC(eltfrnlk,(6+3*natom,6))
+ - + - +
- + - - +
- + ]
306 : : end if
307 [ + + ]: 488 : if (need_becfr) then
308 [ + + + + : 948 : becfrnl(:,:,:)=zero
+ + ]
309 [ - + + - : 108 : ABI_MALLOC(becfrnlk,(3,natom,3))
+ - + - +
- - + -
+ ]
310 : : end if
311 [ + + ]: 488 : if (need_piezofr) then
312 [ + + + + ]: 396 : piezofrnl(:,:)=zero
313 [ - + ]: 18 : ABI_MALLOC(piezofrnlk,(6,3))
314 : : end if
315 : 488 : need_efmas=dtset%efmas>0
316 [ + + + - : 488 : 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 [ - + + - : 1464 : ABI_MALLOC(eig_k,(ndat))
+ - + - -
+ - + ]
324 [ - + + - : 976 : ABI_MALLOC(occ_k,(ndat))
- + - + ]
325 : :
326 : : !Common initialization
327 : 488 : bdtot_index=0;ibg=0;icg=0
328 [ + + ]: 488 : nsploop=dtset%nsppol;if (dtset%nspden==4) nsploop=4
329 : 488 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
330 : 488 : nfftot=ngfftf(1)*ngfftf(2)*ngfftf(3)
331 : :
332 : : !Common data for "nonlop" routine
333 : 488 : tim_nonlop=6
334 [ + + ]: 976 : 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 [ + + ]: 488 : idir_ffnl=0 ; if(need_efmas) idir_ffnl=4
338 : 488 : choice_phon=0;choice_strs=0
339 [ + + ]: 488 : if(rfphon==1)then
340 : 451 : shift_rhoij=0
341 : 451 : choice_phon=4
342 : 451 : nnlout_phon=max(1,6*natom)
343 [ - + + - : 1353 : ABI_MALLOC(enlout_phon,(nnlout_phon*ndat))
+ - + - -
+ - + -
+ ]
344 : : end if
345 [ + + ]: 488 : 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 [ + + ]: 488 : if (psps%usepaw==0) then
352 : 430 : paw_opt=0 ; cpopt=-1
353 : : else
354 : 58 : paw_opt=2 ; cpopt=1+2*usecprj
355 : : end if
356 [ + + ]: 488 : 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 [ + + ]: 488 : 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 : 461 : choice_bec2=0 ; choice_bec54=0
374 : 461 : nnlout_bec1=0
375 : : end if
376 [ + + ]: 488 : 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 : 488 : & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,gpu_option=dtset%gpu_option)
387 : :
388 : : !===== PAW specific section
389 [ + + ]: 488 : if (psps%usepaw==1) then
390 : :
391 : : ! Define several sizes & flags
392 : 58 : ncpgr=0;ngrhoij=0
393 [ + + ]: 58 : if(rfphon==1)then
394 : 45 : ncpgr=3;ngrhoij=3
395 : : end if
396 [ + + + + ]: 58 : if(rfphon==1.or.need_becfr)then
397 : 52 : ncpgr=6;ngrhoij=6
398 : : end if
399 [ + + + - ]: 58 : if(rfstrs/=0.and.rfphon==1)then
400 : 12 : ncpgr=9;ngrhoij=9
401 : : end if
402 [ + + + + ]: 58 : 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 [ + + ]: 58 : 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 [ + + ]: 58 : 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 : 32 : pawrhoij_tot => pawrhoij
445 : : end if
446 : :
447 : : ! Projected WF (cprj) and PAW occupancies (& gradients)
448 [ - + + - : 464 : ABI_MALLOC(cwaveprj,(natom,dtset%nspinor*ndat))
- + + - +
- + - - +
- + - + +
+ + + ]
449 : 58 : call pawcprj_alloc(cwaveprj,ncpgr,gs_ham%dimcprj)
450 [ + + ]: 220 : do iatom=1,natom
451 : 162 : sz2=pawrhoij_tot(iatom)%cplex_rhoij*pawrhoij_tot(iatom)%qphase*pawrhoij_tot(iatom)%lmn2_size
452 : 162 : sz3=pawrhoij_tot(iatom)%nspden
453 [ + + + - : 810 : ABI_MALLOC(pawrhoij_tot(iatom)%grhoij,(ngrhoij,sz2,sz3))
- + + - -
+ + - + -
+ + - + -
+ - + ]
454 : 162 : pawrhoij_tot(iatom)%ngrhoij=ngrhoij
455 [ + + + + : 66114 : pawrhoij_tot(iatom)%grhoij=zero
+ + ]
456 : : end do
457 : 58 : use_timerev=(dtset%kptopt>0.and.dtset%kptopt<3)
458 [ + + - + ]: 112 : 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 [ + + ]: 488 : 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 [ + + ]: 1012 : do isppol=1,dtset%nsppol
478 : :
479 : : ! Continue to initialize the Hamiltonian (PAW DIJ coefficients)
480 : 524 : call gs_ham%load_spin(isppol,with_nonlocal=.true.)
481 : :
482 : : ! Rewind (k+G) data if needed
483 : 524 : ikg=0
484 : :
485 : : ! Loop over k points
486 [ + + ]: 43579 : do ikpt=1,dtset%nkpt
487 : 42567 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
488 : 42567 : 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 : 42567 : & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
491 : :
492 : 42567 : istwf_k=dtset%istwfk(ikpt)
493 : 42567 : npw_k=npwarr(ikpt)
494 : 42567 : wtk_k=dtset%wtk(ikpt)
495 [ + + ]: 170268 : kpoint(:)=dtset%kptns(:,ikpt)
496 : :
497 : : ! Skip this k-point if not the proper processor
498 [ + + ]: 42567 : 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 [ + + ]: 31623 : 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 [ - + + - : 94869 : ABI_MALLOC(cwavef,(2,npw_k*dtset%nspinor*ndat))
+ - + - -
+ - + -
+ ]
525 [ + + ]: 31623 : if (need_becfr.or.need_piezofr) then
526 [ - + + - : 1856 : ABI_MALLOC(svectout,(2,npw_k*dtset%nspinor*ndat))
- + - + -
+ ]
527 : : end if
528 [ + + ]: 31623 : 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 [ - + + - : 126492 : ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
+ + + - +
- + + - +
- + - + ]
535 [ + + + + : 31623 : 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 [ - + - + ]: 25055 : ABI_MALLOC(ylmgr_k,(0,0,0))
539 : : end if
540 : :
541 [ - + + - : 94869 : ABI_MALLOC(kg_k,(3,npw_k))
+ - + - -
+ - + -
+ ]
542 [ + + + + ]: 19068583 : kg_k(:,:) = 0
543 : : !$OMP PARALLEL DO
544 [ + + ]: 4790863 : do ipw=1,npw_k
545 : 4759240 : kg_k(1,ipw)=kg(1,ipw+ikg)
546 : 4759240 : kg_k(2,ipw)=kg(2,ipw+ikg)
547 : 4790863 : kg_k(3,ipw)=kg(3,ipw+ikg)
548 : : end do
549 [ + + ]: 31623 : if (psps%useylm==1) then
550 : : !SOMP PARALLEL DO COLLAPSE(2)
551 [ + + ]: 57480 : do ilm=1,mpsang*mpsang
552 [ + + ]: 6751372 : do ipw=1,npw_k
553 : 6745718 : ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm)
554 : : end do
555 : : end do
556 [ + + + + : 5654 : 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 [ + + ]: 85594 : 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 [ - + ]: 31623 : cplex=2;if (istwf_k>1) cplex=1
569 : :
570 : : ! Compute (k+G) vectors (only if useylm=1)
571 : 31623 : nkpg=0
572 [ + + + + : 31623 : if (rfstrs/=0.or.need_efmas.or.pawpiezo==1) nkpg=3*dtset%nloalg(3)
+ + ]
573 [ - + + - : 126492 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
+ + + - +
- + + - +
- + - + ]
574 [ + + ]: 31623 : 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 [ + + ]: 31623 : 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 : 31623 : ider=0;dimffnl=1;
591 [ + + ]: 31623 : if(need_becfr) then
592 : 928 : ider=1;dimffnl=4
593 : : end if
594 [ + + + + ]: 31623 : if(rfstrs/=0.or.need_piezofr.or.need_efmas)then
595 : 6136 : ider=2;dimffnl=3+7*psps%useylm
596 : : end if
597 [ - + + - : 189738 : 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 : 31623 : & 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 [ + + ]: 31623 : 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 [ - + + - : 126492 : 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 : 31623 : & 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 [ - + ]: 31623 : 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 [ + + + + : 340799 : if(rfphon==1) dyfrnlk(:,:)=zero
+ + ]
641 [ + + ]: 31623 : if(rfstrs/=0)then
642 [ + + + + ]: 386576 : enlk=zero;eltfrnlk(:,:)=zero
643 : : end if
644 [ + + + + : 57447 : if (need_becfr) becfrnlk(:,:,:)=zero
+ + + + ]
645 [ + + + + : 42039 : if (need_piezofr) piezofrnlk(:,:)=zero
+ + ]
646 [ + + ]: 31623 : 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 : 31623 : iband_me = 0
665 [ - + - - : 210795 : do iband=1,nband_k,ndat
- + + + ]
666 : :
667 [ + + ]: 179172 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me)) cycle
668 : 179028 : iband_me = iband_me + ndat
669 : :
670 [ + + ]: 358056 : 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 [ + + + + ]: 122288841 : & iband_me *npw_k*dtset%nspinor+icg)
673 : :
674 : : ! Compute non-local contributions from n,k
675 [ + + + + ]: 193972 : if (psps%usepaw==1) eig_k(1:ndat)=eigen(iband+bdtot_index:iband+bdtot_index+ndat-1)
676 : :
677 : : ! === Dynamical matrix
678 [ + + ]: 179028 : if(rfphon==1) then
679 : : call nonlop(choice_phon,cpopt,cwaveprj,enlout_phon,gs_ham,idir,eig_k,mpi_enreg,ndat,&
680 : 162050 : & nnlout_phon,paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
681 : : ! Accumulate non-local contributions from n,k
682 [ + + ]: 324100 : do idat=1,ndat
683 : : dyfrnlk(:,:) = dyfrnlk(:,:) &
684 [ + + + + : 2490222 : & + 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 [ + + ]: 179028 : 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 [ + + + + ]: 179028 : if (psps%usepaw==1.and.(.not.need_efmas)) then
702 [ + + ]: 29536 : 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 [ - + ]: 29536 : & 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 [ + + ]: 179028 : 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 [ + + ]: 179028 : 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 [ + + ]: 179028 : 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 [ + + ]: 179028 : 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 [ + + ]: 210651 : 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 [ + + ]: 31623 : if(rfphon==1) then
901 [ + + ]: 71711 : do iatom=1,natom
902 [ + + ]: 44168 : ia=iatom;if (dyfr_nondiag==0) ia=1
903 : 44168 : dyfrnl(1,1,1,iatom,ia)=dyfrnl(1,1,1,iatom,ia)+wtk_k*dyfrnlk(1,iatom)
904 : 44168 : dyfrnl(1,2,2,iatom,ia)=dyfrnl(1,2,2,iatom,ia)+wtk_k*dyfrnlk(2,iatom)
905 : 44168 : dyfrnl(1,3,3,iatom,ia)=dyfrnl(1,3,3,iatom,ia)+wtk_k*dyfrnlk(3,iatom)
906 : 44168 : dyfrnl(1,2,3,iatom,ia)=dyfrnl(1,2,3,iatom,ia)+wtk_k*dyfrnlk(4,iatom)
907 : 44168 : dyfrnl(1,1,3,iatom,ia)=dyfrnl(1,1,3,iatom,ia)+wtk_k*dyfrnlk(5,iatom)
908 : 71711 : 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 [ + + ]: 31623 : if(rfstrs/=0) then
912 [ + + + + ]: 386576 : eltfrnl(:,:)=eltfrnl(:,:)+dtset%wtk(ikpt)*eltfrnlk(:,:)
913 : : end if
914 [ + + ]: 31623 : if(need_becfr) then
915 [ + + + + : 26752 : becfrnl(:,:,:)=becfrnl(:,:,:)+dtset%wtk(ikpt)*becfrnlk(:,:,:)
+ + ]
916 : : end if
917 [ + + ]: 31623 : if(need_piezofr) then
918 [ + + + + ]: 10912 : piezofrnl(:,:)=piezofrnl(:,:)+dtset%wtk(ikpt)*piezofrnlk(:,:)
919 : : end if
920 : : ! Increment indexes
921 : 31623 : bdtot_index=bdtot_index+nband_k
922 [ + - ]: 31623 : if (mkmem_rbz/=0) then
923 : 31623 : ibg=ibg+nband_k*dtset%nspinor
924 : 31623 : icg=icg+npw_k*dtset%nspinor*nband_me
925 : 31623 : 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 [ - + ]: 31623 : ABI_FREE(ffnl)
932 [ - + ]: 31623 : ABI_FREE(kpg_k)
933 [ - + ]: 31623 : ABI_FREE(ph3d)
934 [ - + ]: 31623 : ABI_FREE(ylm_k)
935 [ - + ]: 31623 : ABI_FREE(ylmgr_k)
936 [ - + ]: 31623 : ABI_FREE(cwavef)
937 [ - + ]: 31623 : ABI_FREE(kg_k)
938 [ + + ]: 31623 : 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 [ + + ]: 32147 : 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 [ + + ]: 488 : if(rfphon==1) then
958 [ - + ]: 451 : ABI_FREE(dyfrnlk)
959 [ - + ]: 451 : ABI_FREE(enlout_phon)
960 : : end if
961 [ + + ]: 488 : if(rfstrs/=0) then
962 [ - + ]: 37 : ABI_FREE(eltfrnlk)
963 [ - + ]: 37 : ABI_FREE(enlout_strs)
964 : : end if
965 [ + + ]: 488 : 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 [ + + ]: 488 : 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 [ + + ]: 488 : if(need_efmas) then
976 [ - + ]: 17 : ABI_FREE(enlout_efmas)
977 : : end if
978 [ + + ]: 488 : if (psps%usepaw==1) then
979 [ + + ]: 58 : if (need_becfr.or.need_piezofr) then
980 [ - + ]: 27 : ABI_FREE(becij)
981 : : end if
982 : 58 : call pawcprj_free(cwaveprj)
983 : : end if
984 [ + - + + : 658 : ABI_FREE(cwaveprj)
- + - + -
+ ]
985 : :
986 : : !Fill in lower triangle of matrixes
987 [ + + ]: 488 : if (rfphon==1) then
988 [ + + ]: 1393 : do iatom=1,natom
989 [ + + ]: 942 : ia=iatom;if (dyfr_nondiag==0) ia=1
990 : 942 : dyfrnl(1,3,2,iatom,ia)=dyfrnl(1,2,3,iatom,ia)
991 : 942 : dyfrnl(1,3,1,iatom,ia)=dyfrnl(1,1,3,iatom,ia)
992 : 1393 : dyfrnl(1,2,1,iatom,ia)=dyfrnl(1,1,2,iatom,ia)
993 : : end do
994 : : end if
995 [ + + ]: 488 : 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 : 488 : call timab(48,1,tsec)
1006 : : ! Accumulate dyfrnl
1007 [ + + ]: 488 : if(rfphon==1)then
1008 : 451 : call xmpi_sum(dyfrnl,spaceworld,ierr)
1009 : : end if
1010 : : ! Accumulate eltfrnl.
1011 [ + + ]: 488 : if(rfstrs/=0)then
1012 : 37 : call xmpi_sum(eltfrnl,spaceworld,ierr)
1013 : : end if
1014 : : ! Accumulate becfrnl
1015 [ + + ]: 488 : if (need_becfr) then
1016 : 27 : call xmpi_sum(becfrnl,spaceworld,ierr)
1017 : : end if
1018 : : ! Accumulate piezofrnl
1019 [ + + ]: 488 : if (need_piezofr) then
1020 : 18 : call xmpi_sum(piezofrnl,spaceworld,ierr)
1021 : : end if
1022 : :
1023 : : ! PAW: accumulate gradients of rhoij
1024 [ + + ]: 488 : if (psps%usepaw==1) then
1025 [ - + + - : 174 : ABI_MALLOC(dimlmn,(natom))
+ - + - -
+ - + ]
1026 [ + + ]: 220 : dimlmn(1:natom)=pawrhoij_tot(1:natom)%cplex_rhoij*pawrhoij_tot(1:natom)%qphase*pawrhoij_tot(1:natom)%lmn2_size
1027 [ + + ]: 220 : bufdim=ncpgr*sum(dimlmn)*nsploop
1028 [ + + + - : 174 : ABI_MALLOC(mpibuf,(bufdim))
+ - + + -
+ - + -
+ ]
1029 [ + + ]: 52216 : ii=0;mpibuf=zero
1030 [ + + ]: 220 : do iatom=1,natom
1031 [ + + ]: 392 : do isppol=1,nsploop
1032 [ + + ]: 1492 : do mu=1,ncpgr
1033 [ + + ]: 53316 : mpibuf(ii+1:ii+dimlmn(iatom))=pawrhoij_tot(iatom)%grhoij(mu,1:dimlmn(iatom),isppol)
1034 : 1330 : ii=ii+dimlmn(iatom)
1035 : : end do
1036 : : end do
1037 : : end do
1038 : 58 : call xmpi_sum(mpibuf,spaceworld,ierr)
1039 : 58 : ii=0
1040 [ + + ]: 220 : do iatom=1,natom
1041 [ + + ]: 392 : do isppol=1,nsploop
1042 [ + + ]: 1492 : do mu=1,ncpgr
1043 [ + + ]: 53316 : pawrhoij_tot(iatom)%grhoij(mu,1:dimlmn(iatom),isppol)=mpibuf(ii+1:ii+dimlmn(iatom))
1044 : 1330 : ii=ii+dimlmn(iatom)
1045 : : end do
1046 : : end do
1047 : : end do
1048 [ - + ]: 58 : ABI_FREE(mpibuf)
1049 [ - + ]: 116 : ABI_FREE(dimlmn)
1050 : : end if
1051 : 488 : call timab(48,2,tsec)
1052 : : end if
1053 : :
1054 : : !====== PAW: Additional steps
1055 [ + + ]: 488 : 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 [ + + + + ]: 58 : if (rfphon==1.and.rfstrs==0) then
1060 : 33 : 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 : 33 : & 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 [ - + ]: 58 : ABI_CHECK(nsploop/=4,'d2frnl: should we mix mx/my/mz when translating coordinates?')
1073 [ + + ]: 220 : do iatom=1,natom
1074 : 162 : cplx=pawrhoij_tot(iatom)%cplex_rhoij
1075 [ + + ]: 382 : do iq=1,pawrhoij_tot(iatom)%qphase
1076 [ - + ]: 162 : iq0=0;if (iq==2) iq0=cplx*pawrhoij_tot(iatom)%lmn2_size
1077 [ + + ]: 496 : do isppol=1,nsploop
1078 [ + + ]: 8778 : do klmn=1,pawrhoij_tot(iatom)%lmn2_size
1079 [ + + ]: 17788 : do ii=1,cplx
1080 [ + + - + ]: 9172 : if(rfphon==1.or.rfstrs/=0)then
1081 [ + + ]: 26168 : grhoij(1:3)=pawrhoij_tot(iatom)%grhoij(shift_rhoij+1:shift_rhoij+3,iq0+cplx*(klmn-1)+ii,isppol)
1082 [ + + ]: 26168 : do mu=1,3
1083 : : pawrhoij_tot(iatom)%grhoij(shift_rhoij+mu,iq0+cplx*(klmn-1)+ii,isppol)=gprimd(mu,1)*grhoij(1)&
1084 : 28798 : & +gprimd(mu,2)*grhoij(2)+gprimd(mu,3)*grhoij(3)
1085 : : end do
1086 : : end if
1087 [ + + ]: 17616 : 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 [ + + ]: 58 : 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 : 58 : dimnhat=0;optgr=0;optgr2=0;optstr=0;optstr2=0
1120 [ + + ]: 58 : if (rfphon==1) optgr2=1
1121 [ + + ]: 58 : if (rfstrs/=0) optstr2=1
1122 [ - + ]: 58 : 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 : 58 : & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
1127 [ - + ]: 58 : ABI_FREE(nhat_dum)
1128 : : end if !PAW
1129 : :
1130 : : !The indexing array atindx is used to reestablish the correct order of atoms
1131 [ + + ]: 488 : 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 [ + + ]: 488 : 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 [ + + ]: 488 : if (need_piezofr) then
1166 [ + + + + ]: 396 : piezofrnl=-piezofrnl
1167 : : end if
1168 : :
1169 [ - + ]: 488 : ABI_FREE(eig_k)
1170 [ - + ]: 488 : ABI_FREE(occ_k)
1171 : :
1172 : : !Close the ddk files
1173 [ + + ]: 1952 : do ii=1,3
1174 : 1952 : call ddkfiles(ii)%close()
1175 : : end do
1176 : :
1177 : : !Release now useless memory
1178 [ + + ]: 488 : if (psps%usepaw==1) then
1179 [ + + ]: 220 : do iatom=1,natom
1180 [ - + ]: 162 : ABI_FREE(pawrhoij_tot(iatom)%grhoij)
1181 : 220 : pawrhoij_tot(iatom)%ngrhoij=0
1182 : : end do
1183 [ + + ]: 58 : if (paral_atom) then
1184 : 26 : call pawrhoij_free(pawrhoij_tot)
1185 [ - + ]: 26 : ABI_FREE(pawrhoij_tot)
1186 : : end if
1187 : : end if
1188 : 488 : call gs_ham%free()
1189 : :
1190 : : ABI_NVTX_END_RANGE()
1191 : 488 : call timab(159,2,tsec)
1192 : :
1193 : 488 : write(msg,'(3a)')ch10,' ==> Calculation of the frozen part of the second order derivative done',ch10
1194 : 488 : call wrtout(std_out,msg,'COLL')
1195 : :
1196 : : DBG_EXIT("COLL")
1197 : :
1198 [ - + - + : 1464 : end subroutine d2frnl
- + - + -
+ - + - -
- - - - -
- - - - -
- - - - -
- - + - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
+ - + - +
- + - + ]
1199 : : !!***
1200 : :
1201 : : end module m_d2frnl
1202 : : !!***
|
