Line data Source code
1 : !{\src2tex{textfont=tt}}
2 : !!****m* ABINIT/m_dfptlw_pert
3 : !! NAME
4 : !! m_dfptlw_pert
5 : !!
6 : !! FUNCTION
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2022-2026 ABINIT group (MR)
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 : module m_dfptlw_pert
23 :
24 : use defs_basis
25 : use defs_abitypes
26 : use defs_datatypes
27 : use m_dtset
28 : use m_dtfil
29 : use m_errors
30 : use m_profiling_abi
31 : use m_hamiltonian
32 : use m_pawcprj
33 : use m_pawfgr
34 : use m_wfk
35 : use m_xmpi
36 : use m_getgh1c
37 : use m_mklocl
38 : use m_initylmg, only : initylmg
39 : use m_fstrings, only : itoa, sjoin
40 : use m_io_tools, only : file_exists
41 : use m_time, only : cwtime
42 : use m_kg, only : mkkpg
43 : use m_mpinfo, only : proc_distrb_cycle
44 : use m_dfptlw_wf
45 : use m_dfpt_mkvxc, only : dfpt_mkvxcggadq
46 : use m_dfptlw_nv, only : dfptlw_geom
47 : use m_spacepar, only : hartredq
48 : use m_cgtools, only : dotprod_vn
49 : use m_mkffnl, only : mkffnl
50 :
51 : implicit none
52 :
53 : public :: dfptlw_pert
54 : public :: preca_ffnl
55 :
56 : private
57 :
58 : ! *************************************************************************
59 :
60 : contains
61 : !!***
62 :
63 : !!****f* ABINIT/m_dfptlw_pert/dfptlw_pert
64 : !! NAME
65 : !! dfptlw_pert
66 : !!
67 : !! FUNCTION
68 : !! Compute first-order response function contributions to the spatial-dispersion
69 : !! 3rd order energy derivatives of the longwave driver.
70 : !! The main inputs are :
71 : !! - GS WFs and Hamiltonian (cg,gs_hamkq)
72 : !! - 1st-order WFs for two perturbations i1pert/i1dir,i2pert/i2dir (cg1,cg2)
73 : !! - 1st-order Local+SCF potentials for i1pert and i2pert
74 : !! - 1st-order WFs DDK and 2nd-order WF D2_DKDK (d2_dkdk_f)
75 : !!
76 : !! INPUTS
77 : !! cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = array for planewave
78 : !! coefficients of wavefunctions
79 : !! cg1 = first derivative of cg with respect the perturbation i1pert
80 : !! cg2 = first derivative of cg with respect the perturbation i2pert
81 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
82 : !! if 2, COMPLEX
83 : !! dimffnl= third dimension of ffnl
84 : !! dtset <type(dataset_type)>=all input variables for this dataset
85 : !! eigen1(2*mband*mband*nkpt*nsppol)=1st-order eigenvalues for i1pert,i1dir (hartree)
86 : !! eigen2(2*mband*mband*nkpt*nsppol)=1st-order eigenvalues for i2pert,i2dir (hartree)
87 : !! ffnl(dtset%mkmem,dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives
88 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
89 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
90 : !! gsqcut=large sphere cut-off
91 : !! i1dir,i2dir,i3dir=directions of the corresponding perturbations
92 : !! i1pert,i2pert,i3pert = type of perturbation that has to be computed
93 : !! kg(3,mpw*mkmem_rbz)=reduced planewave coordinates
94 : !! kxc(nfft,nkxc)=exchange and correlation kernel
95 : !! mband = maximum number of bands
96 : !! mkmem_rbz = maximum number of k points which can fit in core memory
97 : !! mk1mem = maximum number of k points for first-order WF
98 : !! which can fit in core memory
99 : !! mpert =maximum number of ipert
100 : !! mpi_enreg=MPI-parallelisation information
101 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
102 : !! mpw = maximum number of planewaves in basis sphere (large number)
103 : !! natom = number of atoms in unit cell
104 : !! n1dq= third dimension of vlocal1_i1pertdq
105 : !! n2dq= third dimension of vlocal1_i2pertdq
106 : !! nfft= number of FFT grid points (for this proc)
107 : !! ngfft(1:18)=integer array with FFT box dimensions and other
108 : !! nkpt = number of k points
109 : !! nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed.
110 : !! nspden = number of spin-density components
111 : !! nspinor = number of spinorial components of the wavefunctions
112 : !! nsppol = number of channels for spin-polarization (1 or 2)
113 : !! npwarr(nkpt) = array holding npw for each k point
114 : !! nylmgr=second dimension of ylmgr_k
115 : !! occ(mband*nkpt*nsppol) = occupation number for each band and k
116 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
117 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
118 : !! rho1g1(2,nfft)=G-space RF electron density in electrons/bohr**3 (i1pert)
119 : !! rho1r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i1pert)
120 : !! rho2r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i2pert)
121 : !! rmet(3,3)=real space metric tensor in bohr**2
122 : !! rprimd(3,3) = dimensional primitive translations (bohr)
123 : !! samepert= .true. if i1pert=i2pert and i1dir=i2dir
124 : !! ucvol=volume of the unit cell
125 : !! useylmgr= if 1 use the derivative of spherical harmonics
126 : !! vpsp1_i1pertdq(cplex*nfft,nspden,n1dq)= local potential of first-order
127 : !! gradient Hamiltonian for i1pert
128 : !! vpsp1_i1pertdq(cplex*nfft,nspden,n1dq)= local potential of second-order
129 : !! gradient Hamiltonian for i1pert
130 : !! vpsp1_i1pertdq_geom(cplex*nfft,nspden,3)= local potential of first-order
131 : !! gradient Hamiltonian for i1pert wrp to i3dir and i2dir
132 : !! vpsp1_i2pertdq(cplex*nfft,nspden,n2dq)= local potential of first-order
133 : !! gradient Hamiltonian for i2pert
134 : !! ddk_f = wf files
135 : !! d2_dkdk_f = wf files
136 : !! ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics
137 : !! ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics
138 : !!
139 : !! OUTPUT
140 : !! d3etot(2,3,mpert,3,mpert,3,mpert) = third derivatives of the energy tensor
141 : !! d3etot_t4(2,n2dq)= t4 term which might need to be converted to type-II
142 : !! d3etot_t5(2,n1dq)= t5 term which might need to be converted to type-II
143 : !! d3etot_tgeom(2,2)= Geometric term which needs to be converted to type-II
144 : !!
145 : !! SIDE EFFECTS
146 : !! TO DO!
147 : !!
148 : !! SOURCE
149 :
150 591 : subroutine dfptlw_pert(cg,cg1,cg2,cplex,d3etot,d3etot_t4,d3etot_t5,d3etot_tgeom,&
151 591 : & dimffnl,dtset,eigen1,eigen2,ffnl,gmet,gs_hamkq,gsqcut,i1dir,i2dir,i3dir,&
152 591 : & i1pert,i2pert,i3pert,kg,kxc,mband,mkmem_rbz,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,&
153 : & n1dq,n2dq,nfft,ngfft,nkpt,nkxc,&
154 591 : & nspden,nspinor,nsppol,npwarr,nylmgr,occ,pawfgr,psps,rho1g1,rho1r1,rho2r1,rmet,rprimd,samepert,&
155 591 : & ucvol,useylmgr,vpsp1_i1pertdq,vpsp1_i1pertdqdq,vpsp1_i1pertdq_geom,vpsp1_i2pertdq,ddk_f,d2_dkdk_f,d2_dkdk_f2,ylm,ylmgr)
156 :
157 : !Arguments ------------------------------------
158 : !scalars
159 : integer,intent(in) :: cplex,dimffnl,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband
160 : integer,intent(in) :: mk1mem,mkmem_rbz,mpert,mpsang,mpw,natom,n1dq,n2dq,nfft,nkpt,nkxc,nspden
161 : integer,intent(in) :: nspinor,nsppol,nylmgr,useylmgr
162 : real(dp),intent(in) :: gsqcut,ucvol
163 : logical,intent(in) :: samepert
164 : type(MPI_type),intent(inout) :: mpi_enreg
165 : type(dataset_type),intent(in) :: dtset
166 : type(pseudopotential_type),intent(in) :: psps
167 : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
168 : type(pawfgr_type),intent(in) :: pawfgr
169 : type(wfk_t),intent(inout) :: ddk_f,d2_dkdk_f, d2_dkdk_f2
170 :
171 : !arrays
172 : integer,intent(in) :: kg(3,mpw*mkmem_rbz),ngfft(18),npwarr(nkpt)
173 : real(dp),intent(in) :: eigen1(2*mband*mband*nkpt*nsppol)
174 : real(dp),intent(in) :: eigen2(2*mband*mband*nkpt*nsppol)
175 : real(dp),intent(in) :: ffnl(mkmem_rbz,mpw,dimffnl,psps%lmnmax,psps%ntypat)
176 : real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol)
177 : real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol)
178 : real(dp),intent(in) :: cg2(2,mpw*nspinor*mband*mk1mem*nsppol)
179 : real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc)
180 : real(dp),intent(in) :: occ(mband*nkpt*nsppol)
181 : real(dp),intent(in) :: rho1g1(2,nfft),rho1r1(cplex*nfft,dtset%nspden)
182 : real(dp),intent(in) :: rho2r1(cplex*nfft,dtset%nspden)
183 : real(dp),intent(in) :: rmet(3,3),rprimd(3,3)
184 : real(dp),intent(in) :: vpsp1_i1pertdq(2*nfft,nspden,n1dq)
185 : real(dp),intent(in) :: vpsp1_i1pertdqdq(2*nfft,nspden,n2dq)
186 : real(dp),intent(in) :: vpsp1_i1pertdq_geom(2*nfft,nspden,3)
187 : real(dp),intent(in) :: vpsp1_i2pertdq(2*nfft,nspden,n2dq)
188 : real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert)
189 : real(dp),intent(out) :: d3etot_t4(2,n2dq),d3etot_t5(2,n1dq)
190 : real(dp),intent(out) :: d3etot_tgeom(2,n2dq)
191 : real(dp),intent(in) :: ylm(mpw*mk1mem,psps%mpsang*psps%mpsang*psps%useylm)
192 : real(dp),intent(in) :: ylmgr(mpw*mk1mem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)
193 :
194 : !Variables ------------------------------------
195 : !scalars
196 : integer :: bandtot,bd2tot,icg,idq,ierr,ii,ikc,ikg,ikpt,ilm,isppol,istwf_k,me,n1,n2,n3,n4,n5,n6
197 : integer :: nband_k,npw_k,spaceworld,tim_getgh1c
198 : integer :: usepaw
199 : real(dp) :: tmpim,tmpre,wtk_k
200 : real(dp) :: cpu, wall, gflops
201 : character(len=1000) :: msg
202 : logical :: with_nonlocal_i1pert, with_nonlocal_i2pert
203 : !arrays
204 591 : integer,allocatable :: kg_k(:,:)
205 : real(dp) :: d3etot_t1(2),d3etot_t1_k(2)
206 : real(dp) :: d3etot_t2(2),d3etot_t2_k(2)
207 : real(dp) :: d3etot_t3(2),d3etot_t3_k(2)
208 1182 : real(dp) :: d3etot_t4_k(2,n2dq)
209 1182 : real(dp) :: d3etot_t5_k(2,n1dq)
210 1182 : real(dp) :: d3etot_tgeom_k(2,n2dq)
211 : real(dp) :: d3etot_telec(2)
212 : real(dp) :: e3tot(2),kpt(3)
213 591 : real(dp),allocatable :: eig1_k(:),eig2_k(:),occ_k(:)
214 591 : real(dp),allocatable :: ylm_k(:,:),ylmgr_k(:,:,:)
215 591 : real(dp),allocatable :: ffnl_k(:,:,:,:)
216 :
217 : ! *************************************************************************
218 :
219 : DBG_ENTER("COLL")
220 :
221 591 : write(msg,'(2a,3(a,i2,a,i1))') ch10,'LONGWAVE : ',&
222 1182 : ' perts : ',i1pert,'.',i1dir,' / ',i2pert,'.',i2dir,' / ',i3pert,'.',i3dir
223 591 : call wrtout(std_out,msg,'COLL')
224 591 : call wrtout(ab_out,msg,'COLL')
225 :
226 : !Init parallelism
227 591 : spaceworld=mpi_enreg%comm_cell
228 591 : me=mpi_enreg%me_kpt
229 :
230 : !Additional definitions
231 591 : tim_getgh1c=0
232 591 : usepaw=dtset%usepaw
233 591 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
234 591 : n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6)
235 591 : with_nonlocal_i1pert=.true. ; if (i1pert==natom+2) with_nonlocal_i1pert=.false.
236 591 : with_nonlocal_i2pert=.true. ; if (i2pert==natom+2) with_nonlocal_i2pert=.false.
237 :
238 : !Initialize d3etot parts
239 591 : d3etot_t1=zero
240 591 : d3etot_t2=zero
241 591 : d3etot_t3=zero
242 3093 : d3etot_t4=zero
243 2364 : d3etot_t5=zero
244 : d3etot_telec=zero
245 3093 : d3etot_tgeom=zero
246 :
247 : !Calculate the electrostatic contribution
248 : call lw_elecstic(cplex,d3etot_telec,gmet,gs_hamkq%gprimd,gsqcut,&
249 591 : & i3dir,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,rho1g1,rho1r1,rho2r1,ucvol)
250 :
251 : !Loop over spins
252 591 : bandtot = 0
253 591 : bd2tot = 0
254 591 : icg=0
255 1182 : do isppol = 1, nsppol
256 :
257 : ! Loop over k-points
258 591 : ikg = 0
259 591 : ikc = 0
260 72042 : do ikpt = 1, nkpt
261 :
262 70860 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
263 70860 : npw_k = npwarr(ikpt)
264 70860 : istwf_k = dtset%istwfk(ikpt)
265 :
266 70860 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,mpi_enreg%me)) then
267 0 : bandtot = bandtot + nband_k
268 0 : bd2tot = bd2tot + 2*nband_k**2
269 0 : cycle ! Skip the rest of the k-point loop
270 : end if
271 70860 : ikc= ikc + 1
272 :
273 212580 : ABI_MALLOC(occ_k,(nband_k))
274 356236 : occ_k(:) = occ(1+bandtot:nband_k+bandtot)
275 70860 : wtk_k = dtset%wtk(ikpt)
276 283440 : kpt(:) = dtset%kptns(:,ikpt)
277 :
278 212580 : ABI_MALLOC(eig1_k,(2*nband_k**2))
279 141720 : ABI_MALLOC(eig2_k,(2*nband_k**2))
280 212580 : ABI_MALLOC(kg_k,(3,npw_k))
281 283440 : ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
282 354300 : ABI_MALLOC(ylmgr_k,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
283 425160 : ABI_MALLOC(ffnl_k,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
284 :
285 : !Get plane-wave vectors and related data at k
286 29521220 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
287 70860 : if (dtset%ffnl_lw==1) then
288 22400 : if (psps%useylm==1) then
289 224000 : do ilm=1,psps%mpsang*psps%mpsang
290 20138300 : ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
291 : end do
292 22400 : if (useylmgr==1) then
293 224000 : do ilm=1,psps%mpsang*psps%mpsang
294 2038400 : do ii=1,nylmgr
295 181244700 : ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm)
296 : end do
297 : end do
298 : end if
299 : end if
300 48460 : else if (dtset%ffnl_lw==0) then
301 665791688 : ffnl_k(1:npw_k,:,:,:)=ffnl(ikc,1:npw_k,:,:,:)
302 : end if
303 :
304 : !Get matrix elements for uniform perturbations
305 2434508 : eig1_k(:)=eigen1(1+bd2tot:2*nband_k**2+bd2tot)
306 2434508 : eig2_k(:)=eigen2(1+bd2tot:2*nband_k**2+bd2tot)
307 :
308 : !Compute the stationary terms of d3etot depending on response functions
309 : call dfpt_1wf(cg,cg1,cg2,cplex,ddk_f,d2_dkdk_f,d2_dkdk_f2,d3etot_t1_k,d3etot_t2_k,d3etot_t3_k,&
310 : & d3etot_t4_k,d3etot_t5_k,dimffnl,dtset,eig1_k,eig2_k,ffnl_k,gs_hamkq,icg,&
311 : & i1dir,i2dir,i3dir,i1pert,i2pert,ikpt,isppol,istwf_k,&
312 : & kg_k,kpt,mkmem_rbz,mpi_enreg,mpw,natom,nband_k,&
313 : & n1dq,n2dq,nfft,ngfft,npw_k,nspden,nsppol,nylmgr,occ_k,&
314 : & pawfgr,psps,rmet,rprimd,samepert,useylmgr,&
315 : & vpsp1_i1pertdq,vpsp1_i2pertdq,&
316 70860 : & wtk_k,ylm_k,ylmgr_k)
317 :
318 : ! Add the contribution from each k-point.
319 212580 : d3etot_t1=d3etot_t1 + d3etot_t1_k
320 212580 : d3etot_t2=d3etot_t2 + d3etot_t2_k
321 212580 : d3etot_t3=d3etot_t3 + d3etot_t3_k
322 376752 : d3etot_t4=d3etot_t4 + d3etot_t4_k
323 283440 : d3etot_t5=d3etot_t5 + d3etot_t5_k
324 :
325 : !Compute the nonvariational geometric term
326 70860 : call cwtime(cpu, wall, gflops, "start")
327 70860 : if (i1pert<=natom.and.(i2pert==natom+3.or.i2pert==natom+4)) then
328 : call dfptlw_geom(cg,d3etot_tgeom_k,dimffnl,dtset, &
329 : & ffnl_k,gs_hamkq,icg, &
330 : & i1dir,i2dir,i3dir,i1pert,i2pert,ikpt, &
331 : & isppol,istwf_k,kg_k,kpt,mkmem_rbz,mpi_enreg,natom,mpw,nband_k,n2dq,nfft, &
332 : & ngfft,npw_k,nspden,nsppol,nylmgr,occ_k, &
333 41472 : & psps,rmet,rprimd,useylmgr,vpsp1_i1pertdqdq,vpsp1_i1pertdq_geom,wtk_k,ylm_k,ylmgr_k)
334 :
335 : !Add the contribution from each k-point
336 257484 : d3etot_tgeom=d3etot_tgeom + d3etot_tgeom_k
337 : end if
338 70860 : call cwtime(cpu, wall, gflops, "stop")
339 :
340 : ! Keep track of total number of bands
341 70860 : bandtot = bandtot + nband_k
342 70860 : bd2tot = bd2tot + 2*nband_k**2
343 :
344 : ! Shift arrays memory
345 70860 : icg=icg+npw_k*dtset%nspinor*nband_k
346 70860 : ikg=ikg+npw_k
347 :
348 70860 : ABI_FREE(eig1_k)
349 70860 : ABI_FREE(eig2_k)
350 70860 : ABI_FREE(occ_k)
351 70860 : ABI_FREE(kg_k)
352 70860 : ABI_FREE(ylm_k)
353 70860 : ABI_FREE(ylmgr_k)
354 213171 : ABI_FREE(ffnl_k)
355 :
356 : end do !ikpt
357 :
358 : end do !isppol
359 :
360 :
361 : !=== MPI communications ==================
362 : if (xmpi_paral==1) then
363 591 : call xmpi_sum(d3etot_t1,spaceworld,ierr)
364 591 : call xmpi_sum(d3etot_t2,spaceworld,ierr)
365 591 : call xmpi_sum(d3etot_t3,spaceworld,ierr)
366 591 : call xmpi_sum(d3etot_t4,spaceworld,ierr)
367 591 : call xmpi_sum(d3etot_t5,spaceworld,ierr)
368 591 : call xmpi_sum(d3etot_tgeom,spaceworld,ierr)
369 : end if
370 :
371 : !Apply +i or -i in case of strain perturbation.
372 591 : if (i1pert==natom+3.or.i1pert==natom+4) then
373 0 : tmpre=d3etot_telec(1);tmpim=d3etot_telec(2) ; d3etot_telec(1)=tmpim;d3etot_telec(2)=-tmpre
374 0 : tmpre=d3etot_t1(1);tmpim=d3etot_t1(2) ; d3etot_t1(1)=tmpim;d3etot_t1(2)=-tmpre
375 0 : tmpre=d3etot_t2(1);tmpim=d3etot_t2(2) ; d3etot_t2(1)=tmpim;d3etot_t2(2)=-tmpre
376 0 : tmpre=d3etot_t3(1);tmpim=d3etot_t3(2) ; d3etot_t3(1)=tmpim;d3etot_t3(2)=-tmpre
377 0 : do idq=1,n2dq
378 0 : tmpre=d3etot_t4(1,idq);tmpim=d3etot_t4(2,idq) ; d3etot_t4(1,idq)=tmpim;d3etot_t4(2,idq)=-tmpre
379 : end do
380 0 : do idq=1,n1dq
381 591 : tmpre=d3etot_t5(1,idq);tmpim=d3etot_t5(2,idq) ; d3etot_t5(1,idq)=tmpim;d3etot_t5(2,idq)=-tmpre
382 : end do
383 : end if
384 591 : if (i2pert==natom+3.or.i2pert==natom+4) then
385 486 : tmpre=d3etot_telec(1);tmpim=d3etot_telec(2) ; d3etot_telec(1)=-tmpim;d3etot_telec(2)=tmpre
386 486 : tmpre=d3etot_t1(1);tmpim=d3etot_t1(2) ; d3etot_t1(1)=-tmpim;d3etot_t1(2)=tmpre
387 486 : tmpre=d3etot_t2(1);tmpim=d3etot_t2(2) ; d3etot_t2(1)=-tmpim;d3etot_t2(2)=tmpre
388 486 : tmpre=d3etot_t3(1);tmpim=d3etot_t3(2) ; d3etot_t3(1)=-tmpim;d3etot_t3(2)=tmpre
389 1215 : do idq=1,n2dq
390 729 : tmpre=d3etot_t4(1,idq);tmpim=d3etot_t4(2,idq) ; d3etot_t4(1,idq)=-tmpim;d3etot_t4(2,idq)=tmpre
391 1215 : if (i1pert<=natom) then
392 486 : tmpre=d3etot_tgeom(1,idq);tmpim=d3etot_tgeom(2,idq) ; d3etot_tgeom(1,idq)=-tmpim;d3etot_tgeom(2,idq)=tmpre
393 : end if
394 : end do
395 972 : do idq=1,n1dq
396 1077 : tmpre=d3etot_t5(1,idq);tmpim=d3etot_t5(2,idq) ; d3etot_t5(1,idq)=-tmpim;d3etot_t5(2,idq)=tmpre
397 : end do
398 : end if
399 :
400 : !Join all the contributions in e3tot except t4 and t5 which may need to be
401 : !converted to type-II in case of strain perturbation.
402 : !Apply here the two factor to the stationary wf1 contributions
403 : !(see PRB 105, 064101 (2022))
404 1773 : d3etot_t1(:)=two*d3etot_t1(:)
405 1773 : d3etot_t2(:)=two*d3etot_t2(:)
406 1773 : d3etot_t3(:)=two*d3etot_t3(:)
407 3093 : d3etot_t4(:,:)=two*d3etot_t4(:,:)
408 2364 : d3etot_t5(:,:)=two*d3etot_t5(:,:)
409 1773 : e3tot(:)=d3etot_t1(:)+d3etot_t2(:)+d3etot_t3(:)+d3etot_telec(:)
410 :
411 :
412 : !Before printing, set small contributions to zero
413 : !Real parts
414 591 : if (abs(d3etot_t1(1))<tol8) d3etot_t1(1)= zero
415 591 : if (abs(d3etot_t2(1))<tol8) d3etot_t2(1)= zero
416 591 : if (abs(d3etot_t3(1))<tol8) d3etot_t3(1)= zero
417 1425 : do idq=1,n2dq
418 834 : if (abs(d3etot_t4(1,idq))<tol8) d3etot_t4(1,idq)= zero
419 1425 : if (abs(d3etot_tgeom(1,idq))<tol8) d3etot_tgeom(1,idq)= zero
420 : end do
421 1182 : do idq=1,n1dq
422 1182 : if (abs(d3etot_t5(1,idq))<tol8) d3etot_t5(1,idq)= zero
423 : end do
424 591 : if (abs(d3etot_telec(1))<tol8) d3etot_telec(1)= zero
425 591 : if (abs(e3tot(1)) <tol8) e3tot(1)= zero
426 :
427 : !Imaginary parts
428 591 : if (abs(d3etot_t1(2))<tol8) d3etot_t1(2)= zero
429 591 : if (abs(d3etot_t2(2))<tol8) d3etot_t2(2)= zero
430 591 : if (abs(d3etot_t3(2))<tol8) d3etot_t3(2)= zero
431 1425 : do idq=1,n2dq
432 834 : if (abs(d3etot_t4(2,idq))<tol8) d3etot_t4(2,idq)= zero
433 1425 : if (abs(d3etot_tgeom(2,idq))<tol8) d3etot_tgeom(2,idq)= zero
434 : end do
435 1182 : do idq=1,n1dq
436 1182 : if (abs(d3etot_t5(2,idq))<tol8) d3etot_t5(2,idq)= zero
437 : end do
438 591 : if (abs(d3etot_telec(2))<tol8) d3etot_telec(2)= zero
439 591 : if (abs(e3tot(2)) <tol8) e3tot(2)= zero
440 :
441 591 : if (dtset%prtvol>=10) then
442 : write(msg,'(4(a,2(a,f18.8)),a)') &
443 0 : ch10,' d3etot_telec = ',d3etot_telec(1), ',',d3etot_telec(2),&
444 0 : ch10,' d3etot_t1 = ',d3etot_t1(1), ',',d3etot_t1(2),&
445 0 : ch10,' d3etot_t2 = ',d3etot_t2(1), ',',d3etot_t2(2),&
446 0 : ch10,' d3etot_t3 = ',d3etot_t3(1), ',',d3etot_t3(2)
447 0 : call wrtout(std_out,msg,'COLL')
448 0 : call wrtout(ab_out,msg,'COLL')
449 0 : if (n2dq==1) then
450 : write(msg,'(2(a,f18.8))') &
451 0 : ' d3etot_t4 = ',d3etot_t4(1,1), ',',d3etot_t4(2,1)
452 0 : else if (n2dq==2) then
453 : write(msg,'(2(2(a,f18.8),a))') &
454 0 : ' d3etot_t4(dw shear) = ',d3etot_t4(1,1), ',',d3etot_t4(2,1),ch10,&
455 0 : ' d3etot_t4(up shear) = ',d3etot_t4(1,2), ',',d3etot_t4(2,2)
456 : end if
457 0 : call wrtout(std_out,msg,'COLL')
458 0 : call wrtout(ab_out,msg,'COLL')
459 0 : if (n1dq==1) then
460 : write(msg,'(2(a,f18.8))') &
461 0 : ' d3etot_t5 = ',d3etot_t5(1,1), ',',d3etot_t5(2,1)
462 0 : else if (n1dq==2) then
463 : write(msg,'(2(2(a,f18.8),a))') &
464 0 : ' d3etot_t5(dw shear) = ',d3etot_t5(1,1), ',',d3etot_t5(2,1),ch10,&
465 0 : ' d3etot_t5(up shear) = ',d3etot_t5(1,2), ',',d3etot_t5(2,2)
466 : end if
467 0 : call wrtout(std_out,msg,'COLL')
468 0 : call wrtout(ab_out,msg,'COLL')
469 0 : if (i1pert<=natom.and.(i2pert==natom+3.or.i2pert==natom+4)) then
470 0 : if (n2dq==1) then
471 : write(msg,'(2(a,f18.8))') &
472 0 : ' d3etot_tgeom = ',d3etot_tgeom(1,1), ',',d3etot_tgeom(2,1)
473 0 : else if (n2dq==2) then
474 : write(msg,'(2(2(a,f18.8),a))') &
475 0 : 'd3etot_tgeom(dw shear) = ',d3etot_tgeom(1,1), ',',d3etot_tgeom(2,1),ch10,&
476 0 : 'd3etot_tgeom(up shear) = ',d3etot_tgeom(1,2), ',',d3etot_tgeom(2,2)
477 : end if
478 0 : call wrtout(std_out,msg,'COLL')
479 0 : call wrtout(ab_out,msg,'COLL')
480 : end if
481 : end if
482 :
483 1773 : d3etot(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=e3tot(:)
484 :
485 : !Deallocations
486 :
487 : DBG_EXIT("COLL")
488 :
489 591 : end subroutine dfptlw_pert
490 : !!***
491 :
492 : !!****f* ABINIT/m_dfptlw_pert/lw_elecstic
493 : !! NAME
494 : !! lw_elecstic
495 : !!
496 : !! FUNCTION
497 : !! This routine calculates the electrostatic term of the spatial-dispersion
498 : !! third-order energy derivative for a couple of perturbations and a gradient
499 : !! direction.
500 : !!
501 : !! INPUTS
502 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
503 : !! if 2, COMPLEX
504 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
505 : !! gprimd(3,3)=reciprocal space dimensional primitive translations
506 : !! gsqcut=large sphere cut-off
507 : !! i3dir= directions of the 3th perturbations
508 : !! kxc(nfft,nkxc)=exchange and correlation kernel
509 : !! mpi_enreg=information about MPI parallelization
510 : !! nfft= number of FFT grid points (for this proc)
511 : !! ngfft(1:18)=integer array with FFT box dimensions and other
512 : !! nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed.
513 : !! nspden = number of spin-density components
514 : !! rho1g1(2,nfft)=G-space RF electron density in electrons/bohr**3 (i1pert)
515 : !! rho1r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i1pert)
516 : !! rho2r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i2pert)
517 : !! ucvol=volume of the unit cell
518 : !!
519 : !! OUTPUT
520 : !! d3etot_telec(2)= Electrostatic term of the third-order energy derivative
521 : !!
522 : !! SOURCE
523 :
524 591 : subroutine lw_elecstic(cplex,d3etot_telec,gmet,gprimd,gsqcut,&
525 591 : & i3dir,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,rho1g1,rho1r1,rho2r1,ucvol)
526 :
527 : !Arguments ------------------------------------
528 : integer,intent(in) :: cplex,i3dir
529 : integer,intent(in) :: nfft,nkxc,nspden
530 : real(dp),intent(in) :: gsqcut,ucvol
531 : type(MPI_type),intent(inout) :: mpi_enreg
532 : !arrays
533 : integer,intent(in) :: ngfft(18)
534 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
535 : real(dp),intent(in) :: rho1g1(2,nfft),rho1r1(cplex*nfft,nspden)
536 : real(dp),intent(in) :: rho2r1(cplex*nfft,nspden),kxc(nfft,nkxc)
537 : real(dp),intent(out) :: d3etot_telec(2)
538 :
539 : !Local variables-------------------------------
540 : !scalars
541 : integer :: ii,jj,nfftot,qcar
542 : real(dp) :: doti,dotr
543 :
544 : !arrays
545 591 : real(dp),allocatable :: rhor1_cplx(:,:)
546 591 : real(dp),allocatable :: vxc1dq(:,:),vxc1dq_car(:,:,:),vqgradhart(:)
547 :
548 : ! *************************************************************************
549 :
550 : DBG_ENTER("COLL")
551 :
552 : !If GGA xc first calculate the Cartesian q gradient of the xc potential
553 591 : if (nkxc == 7) then
554 995 : ABI_MALLOC(vxc1dq,(2*nfft,nspden))
555 995 : ABI_MALLOC(vxc1dq_car,(2*nfft,nspden,3))
556 796 : do qcar=1,3
557 597 : call dfpt_mkvxcggadq(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,qcar,rho1r1,vxc1dq)
558 4134967 : vxc1dq_car(:,:,qcar)=vxc1dq(:,:)
559 : end do
560 : end if
561 :
562 : !Calculate the q gradient of the Hartree potential
563 2364 : ABI_MALLOC(vqgradhart,(2*nfft))
564 591 : call hartredq(2,gmet,gsqcut,mpi_enreg,nfft,ngfft,i3dir,rho1g1,vqgradhart)
565 :
566 : !If GGA convert the gradient of xc kernel to reduced coordinates and incorporate it to the Hartree part
567 591 : if (nkxc == 7) then
568 1378256 : vxc1dq=zero
569 796 : do qcar=1,3
570 4134967 : vxc1dq(:,:)=vxc1dq(:,:) + gprimd(qcar,i3dir) * vxc1dq_car(:,:,qcar)
571 : end do
572 1378057 : vqgradhart(:)=vqgradhart(:)+vxc1dq(:,1)
573 199 : ABI_FREE(vxc1dq_car)
574 : end if
575 :
576 : !Calculate the electrostatic energy term with the i2pert density response
577 : !I need a complex density for the dotprod_vn
578 2364 : ABI_MALLOC(rhor1_cplx,(2*nfft,nspden))
579 5568228 : rhor1_cplx=zero
580 2784114 : do ii=1,nfft
581 2783523 : jj=ii*2
582 5567637 : rhor1_cplx(jj-1,:)=rho2r1(ii,:)
583 : end do
584 :
585 591 : nfftot=ngfft(1)*ngfft(2)*ngfft(3)
586 591 : call dotprod_vn(2,rhor1_cplx,dotr,doti,nfft,nfftot,nspden,2,vqgradhart,ucvol)
587 :
588 591 : d3etot_telec(1)=dotr
589 591 : d3etot_telec(2)=doti
590 :
591 : !Deallocations
592 591 : ABI_SFREE(vxc1dq)
593 591 : ABI_FREE(vqgradhart)
594 591 : ABI_FREE(rhor1_cplx)
595 :
596 : DBG_EXIT("COLL")
597 :
598 591 : end subroutine lw_elecstic
599 : !!***
600 :
601 : !!****f* ABINIT/m_dfptlw_pert/preca_ffnl
602 : !! NAME
603 : !! preca_ffnl
604 : !!
605 : !! FUNCTION
606 : !! Calculates the nonlocal form factors and derivatives for all the atoms
607 : !! and k points.
608 : !!
609 : !! INPUTS
610 : !! dimffnl= second dimension of ffnl
611 : !! gmet(3,3)= reciprocal-space metric tensor
612 : !! gprimd(3,3)= dimensional reciprocal space primitive translations (b^-1)
613 : !! ider= if 1 first order derivatives of ffnl are calculated
614 : !! if 2 first and second order derivatives of ffnl are calculated
615 : !! idir0= variable that controls the way in which the derivatives of ffnl are
616 : !! calculated and saved
617 : !! kg(3,mpw)=integer coordinates of G vectors in basis sphere
618 : !! kptns(3,nkpt)=k points in terms of reciprocal translations
619 : !! mband= masimum number of bands
620 : !! mkmem= maximum number of k points which can fit in core memory
621 : !! mpi_enreg=information about MPI parallelization
622 : !! mpw = maximum number of planewaves in basis sphere (large number)
623 : !! nkpt = number of k point
624 : !! npwarr(nkpt)=array holding npw for each k point
625 : !! nylmgr=second dimension of ylmgr
626 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
627 : !! rmet(3,3)= real-space metric tensor
628 : !! useylmgr= if 1 use the derivative of spherical harmonics
629 : !! ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics
630 : !! ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics
631 : !!
632 : !! OUTPUT
633 : !! ffnl(mkmem,npw_k,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives
634 : !!
635 : !! SIDE EFFECTS
636 : !!
637 : !! NOTES
638 : !!
639 : !! SOURCE
640 :
641 12 : subroutine preca_ffnl(dimffnl,ffnl,gmet,gprimd,ider,idir0,kg,kptns,mband,mkmem,mpi_enreg,mpw,nkpt, &
642 12 : & npwarr,nylmgr,psps,rmet,useylmgr,ylm,ylmgr)
643 :
644 : !Arguments ------------------------------------
645 : !scalars
646 : integer , intent(in) :: dimffnl,ider,idir0,mband,mkmem,mpw,nkpt,nylmgr,useylmgr
647 : type(pseudopotential_type),intent(in) :: psps
648 : type(MPI_type),intent(in) :: mpi_enreg
649 : !arrays
650 : integer,intent(in) :: kg(3,mpw*mkmem)
651 : integer,intent(in) :: npwarr(nkpt)
652 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kptns(3,nkpt),rmet(3,3)
653 : real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
654 : real(dp),intent(in) :: ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)
655 : real(dp),intent(out) :: ffnl(mkmem,mpw,dimffnl,psps%lmnmax,psps%ntypat)
656 :
657 : !Local variables-------------------------------
658 : !scalars
659 : integer :: ii,ikc,ikg,ikpt,ilm,nkpg,npw_k
660 : !character(len=500) :: msg
661 : !arrays
662 12 : integer,allocatable :: kg_k(:,:)
663 : real(dp) :: kpt(3)
664 12 : real(dp),allocatable :: ffnl_k(:,:,:,:),kpg_k(:,:)
665 12 : real(dp),allocatable :: ylm_k(:,:),ylmgr_k(:,:,:),ylmgr_k_part(:,:,:)
666 :
667 : ! *************************************************************************
668 :
669 : DBG_ENTER("COLL")
670 :
671 : !Loop over k-points
672 12 : ikg=0
673 12 : ikc=0
674 876 : do ikpt = 1, nkpt
675 :
676 864 : npw_k = npwarr(ikpt)
677 864 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,1,mpi_enreg%me)) then
678 : cycle ! Skip the rest of the k-point loop
679 : end if
680 864 : ikc= ikc + 1
681 :
682 2592 : ABI_MALLOC(kg_k,(3,npw_k))
683 3456 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
684 4320 : ABI_MALLOC(ylmgr_k,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
685 :
686 3456 : kpt(:)= kptns(:,ikpt)
687 :
688 : !Get plane-wave vectors and related data at k
689 572632 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
690 864 : if (psps%useylm==1) then
691 9496 : do ilm=1,psps%mpsang*psps%mpsang
692 1404265 : ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
693 : end do
694 864 : if (useylmgr==1) then
695 9496 : do ilm=1,psps%mpsang*psps%mpsang
696 87184 : do ii=1,nylmgr
697 12639241 : ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm)
698 : end do
699 : end do
700 : end if
701 : end if
702 :
703 864 : if (dimffnl==2.or.dimffnl==4) then
704 2432 : ABI_MALLOC(ylmgr_k_part,(npw_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
705 3527451 : ylmgr_k_part(:,:,:)=ylmgr_k(:,1:3,:)
706 256 : else if (dimffnl==10) then
707 1024 : ABI_MALLOC(ylmgr_k_part,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
708 2071624 : ylmgr_k_part(:,:,:)=ylmgr_k(:,:,:)
709 : else
710 0 : ABI_BUG('wrong size for ffnl via dimffnl!')
711 : end if
712 :
713 :
714 864 : nkpg=0
715 1728 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
716 5184 : ABI_MALLOC(ffnl_k,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
717 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl_k,psps%ffspl,gmet,gprimd,ider,idir0,&
718 : & psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
719 864 : & npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k_part)
720 :
721 12905392 : ffnl(ikc,1:npw_k,:,:,:)=ffnl_k(1:npw_k,:,:,:)
722 :
723 864 : ABI_FREE(kg_k)
724 864 : ABI_FREE(ylm_k)
725 864 : ABI_FREE(ylmgr_k)
726 864 : ABI_FREE(ylmgr_k_part)
727 864 : ABI_FREE(ffnl_k)
728 864 : ABI_FREE(kpg_k)
729 :
730 : !Shift arrays memory
731 876 : ikg=ikg+npw_k
732 :
733 : end do
734 :
735 : DBG_EXIT("COLL")
736 :
737 12 : end subroutine preca_ffnl
738 : !!***
739 :
740 : end module m_dfptlw_pert
741 : !!***
|
