Line data Source code
1 : !{\src2tex{textfont=tt}}
2 : !!****m* ABINIT/m_dfptlw_wf
3 : !! NAME
4 : !! m_dfptlw_wf
5 : !!
6 : !! FUNCTION
7 : !! FIXME: add description.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2022-2026 ABINIT group (MR)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_dfptlw_wf
24 :
25 : use defs_basis
26 : use defs_abitypes
27 : use defs_datatypes
28 : use m_dtset
29 : use m_errors
30 : use m_profiling_abi
31 : use m_hamiltonian
32 : use m_cgtools
33 : use m_pawcprj
34 : use m_pawfgr
35 : use m_wfk
36 : use m_xmpi
37 : use m_getgh1c
38 : use m_mklocl
39 :
40 : use m_fstrings, only : itoa, sjoin
41 : use m_io_tools, only : file_exists
42 : use m_time, only : cwtime
43 : use m_kg, only : mkkpg
44 :
45 : implicit none
46 :
47 : public :: dfpt_1wf
48 :
49 : private
50 :
51 : ! *************************************************************************
52 :
53 : contains
54 : !!***
55 :
56 : !!****f* ABINIT/m_dfptlw_wf/dfpt_1wf
57 : !! NAME
58 : !! dfpt_1wf
59 : !!
60 : !! FUNCTION
61 : !! Compute the spin, band and kpt resolved contributions
62 : !! to the spatial-dispersion third-order energy derivatives
63 : !! that depend on first-order response functions.
64 : !!
65 : !! INPUTS
66 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k
67 : !! cplex: if 1, several magnitudes are REAL, if 2, COMPLEX
68 : !! ddk_f = wf files
69 : !! d2_dkdk_f = wf files
70 : !! dimffnl= third dimension of ffnl_k
71 : !! dtset <type(dataset_type)>=all input variables for this dataset
72 : !! eig1_k(2*nband_k**2)=1st-order eigenvalues at k for i1pert,i1dir
73 : !! eig2_k(2*nband_k**2)=1st-order eigenvalues at k for i2pert,i2dir
74 : !! ffnl_k(dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives for this k point
75 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
76 : !! cg1 = first derivative of cg with respect the perturbation i1pert
77 : !! cg2 = first derivative of cg with respect the perturbation i2pert
78 : !! icg=shift to be applied on the location of data in the array cg
79 : !! i1dir,i2dir,i3dir=directions of the corresponding perturbations
80 : !! i1pert,i2pert,i3pert = type of perturbation that has to be computed
81 : !! ikpt=number of the k-point
82 : !! isppol=1 for unpolarized, 2 for spin-polarized
83 : !! istwf_k=parameter that describes the storage of wfs
84 : !! kg_k(3,npw_k)=reduced planewave coordinates.
85 : !! kpt(3)=reduced coordinates of k point
86 : !! mkmem =number of k points treated by this node
87 : !! mpi_enreg=information about MPI parallelization
88 : !! mpw=maximum dimensioned size of npw or wfs at k
89 : !! natom= number of atoms in the unit cell
90 : !! natpert=number of atomic displacement perturbations
91 : !! nband_k=number of bands at this k point for that spin polarization
92 : !! nfft=(effective) number of FFT grid points (for this proc)
93 : !! ngfft(1:18)=integer array with FFT box dimensions and other
94 : !! npw_k=number of plane waves at this k point
95 : !! nspden=number of spin-density components
96 : !! nsppol=1 for unpolarized, 2 for spin-polarized
97 : !! nylmgr=second dimension of ylmgr_k
98 : !! occ_k(nband_k)=occupation number for each band (usually 2) for each k.
99 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
100 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
101 : !! rmet(3,3)=real space metric (bohr**2)
102 : !! rprimd(3,3) = dimensional primitive translations (bohr)
103 : !! samepert= .true. if i1pert=i2pert and i1dir=i2dir
104 : !! useylmgr= if 1 use the derivative of spherical harmonics
105 : !! vpsp1_i1pertdq(cplex*nfft,nspden,n1dq)= local potential of first-order
106 : !! gradient Hamiltonian for i1pert
107 : !! vpsp1_i2pertdq(cplex*nfft,nspden,n2dq)= local potential of first-order
108 : !! gradient Hamiltonian for i2pert
109 : !! wtk_k=weight assigned to the k point.
110 : !! ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics for the k point
111 : !! ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical
112 : !! harmonics for the k point
113 : !!
114 : !! OUTPUT
115 : !! d3etot_t(1-5)_k= stationary 1wf contributions to the third-order energy
116 : !! derivatives for kpt
117 : !!
118 : !! SOURCE
119 :
120 70860 : subroutine dfpt_1wf(cg,cg1,cg2,cplex,ddk_f,d2_dkdk_f,d2_dkdk_f2,&
121 : & d3etot_t1_k,d3etot_t2_k,d3etot_t3_k,&
122 70860 : & d3etot_t4_k,d3etot_t5_k,dimffnl,dtset,eig1_k,eig2_k,ffnl_k,gs_hamkq,icg,&
123 : & i1dir,i2dir,i3dir,i1pert,i2pert,ikpt,isppol,istwf_k,&
124 70860 : & kg_k,kpt,mkmem,mpi_enreg,mpw,natom,nband_k,&
125 70860 : & n1dq,n2dq,nfft,ngfft,npw_k,nspden,nsppol,nylmgr,occ_k,&
126 : & pawfgr,psps,rmet,rprimd,samepert,useylmgr,&
127 70860 : & vpsp1_i1pertdq,vpsp1_i2pertdq,&
128 70860 : & wtk_k,ylm_k,ylmgr_k)
129 :
130 : !Arguments ------------------------------------
131 : !scalars
132 : integer,intent(in) :: cplex,dimffnl,i1dir,i1pert,i2dir,i2pert,i3dir
133 : integer,intent(in) :: icg,ikpt,isppol,istwf_k
134 : integer,intent(in) :: mkmem,mpw,natom,nband_k,n1dq,n2dq,nfft
135 : integer,intent(in) :: npw_k,nspden,nsppol,nylmgr
136 : integer,intent(in) :: useylmgr
137 : real(dp),intent(in) :: wtk_k
138 : logical,intent(in) :: samepert
139 : type(dataset_type),intent(in) :: dtset
140 : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
141 : type(MPI_type),intent(in) :: mpi_enreg
142 : type(pseudopotential_type),intent(in) :: psps
143 : type(wfk_t),intent(inout) :: ddk_f,d2_dkdk_f, d2_dkdk_f2
144 : type(pawfgr_type),intent(in) :: pawfgr
145 :
146 : !arrays
147 : integer,intent(in) :: kg_k(3,npw_k),ngfft(18)
148 : real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol)
149 : real(dp),intent(in) :: cg1(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol)
150 : real(dp),intent(in) :: cg2(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol)
151 : real(dp),intent(out) :: d3etot_t1_k(2)
152 : real(dp),intent(out) :: d3etot_t2_k(2)
153 : real(dp),intent(out) :: d3etot_t3_k(2)
154 : real(dp),intent(out) :: d3etot_t4_k(2,n2dq)
155 : real(dp),intent(out) :: d3etot_t5_k(2,n1dq)
156 : real(dp),intent(in) :: eig1_k(2*nband_k**2),eig2_k(2*nband_k**2)
157 : real(dp),intent(in) :: ffnl_k(npw_k,dimffnl,psps%lmnmax,psps%ntypat)
158 : real(dp),intent(in) :: kpt(3),occ_k(nband_k)
159 : real(dp),intent(in) :: rmet(3,3),rprimd(3,3)
160 : real(dp),intent(in) :: vpsp1_i1pertdq(2*nfft,nspden,n1dq)
161 : real(dp),intent(in) :: vpsp1_i2pertdq(2*nfft,nspden,n2dq)
162 : real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)
163 : real(dp),intent(in) :: ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)
164 :
165 : !Local variables-------------------------------
166 : !scalars
167 : integer :: berryopt,dimffnlk,dimffnl1,iband,idir,idq,ii,jband,nkpg,nkpg1,nylmgrtmp
168 : integer :: offset_cgi,opt_gvnl1,optlocal,optnl,reuse_ffnlk,reuse_ffnl1,sij_opt
169 : integer :: size_wf,tim_getgh1c,usepaw,usevnl,useylmgr1
170 : real(dp) :: cprodi,cprodr,doti,dotr,dum_lambda,fac,tmpim,tmpre
171 : real(dp) :: cpu,wall,gflops
172 : logical :: with_nonlocal_i1pert,with_nonlocal_i2pert
173 70860 : type(rf_hamiltonian_type) :: rf_hamkq
174 :
175 : !arrays
176 70860 : real(dp),allocatable :: cg1_aux(:,:),cg1_ddk(:,:,:),cwave0i(:,:),cwave0j(:,:)
177 70860 : real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:)
178 70860 : real(dp),allocatable :: dkinpw(:),gv1c(:,:)
179 70860 : real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:)
180 70860 : real(dp),allocatable :: gvloc1dqc(:,:),gvnl1dqc(:,:)
181 : real(dp) :: cj_h1_ci(2),dum_grad_berry(1,1),dum_gs1(1,1),dum_gvnl1(1,1)
182 70860 : real(dp),allocatable :: kinpw1(:),kpg_k(:,:),kpg1_k(:,:)
183 70860 : real(dp),allocatable :: part_ylmgr_k(:,:,:),ph3d(:,:,:),ph3d1(:,:,:)
184 70860 : real(dp),allocatable :: dum_vlocal(:,:,:,:),vlocal1(:,:,:,:),dum_vpsp(:)
185 70860 : real(dp),allocatable :: vpsp1(:)
186 70860 : type(pawcprj_type),allocatable :: dum_cwaveprj(:,:)
187 :
188 : ! *************************************************************************
189 :
190 : DBG_ENTER("COLL")
191 :
192 : !Additional definitions
193 70860 : tim_getgh1c=0
194 70860 : useylmgr1=useylmgr
195 70860 : usepaw=dtset%usepaw
196 70860 : size_wf= dtset%nspinor*npw_k
197 70860 : with_nonlocal_i1pert=.true. ; if (i1pert==natom+2) with_nonlocal_i1pert=.false.
198 70860 : with_nonlocal_i2pert=.true. ; if (i2pert==natom+2) with_nonlocal_i2pert=.false.
199 70860 : reuse_ffnlk=1 ; if (dtset%ffnl_lw==1) reuse_ffnlk=0
200 70860 : reuse_ffnl1=1 ; if (dtset%ffnl_lw==1) reuse_ffnl1=0
201 :
202 : !Additional allocations
203 212580 : ABI_MALLOC(cwave0i,(2,size_wf))
204 141720 : ABI_MALLOC(cwave0j,(2,size_wf))
205 141720 : ABI_MALLOC(cwavef1,(2,size_wf))
206 141720 : ABI_MALLOC(cwavef2,(2,size_wf))
207 141720 : ABI_MALLOC(cg1_aux,(2,size_wf))
208 141720 : ABI_MALLOC(gv1c,(2,size_wf))
209 425160 : ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc))
210 212580 : ABI_MALLOC(dum_vpsp,(nfft))
211 425160 : ABI_MALLOC(dum_vlocal,(ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc))
212 212580 : ABI_MALLOC(vpsp1,(cplex*nfft))
213 70860 : ABI_MALLOC(dum_cwaveprj,(0,0))
214 283440 : ABI_MALLOC(part_ylmgr_k,(npw_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
215 60620160 : part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3,:)
216 :
217 : !------------------------------------T1------------------------------------------------
218 : !q1-gradient of gs Hamiltonian:
219 : ! < u_{i,k}^{\lambda1}} | \partial_{gamma} H^{(0)} | u_{i,k}^{\lambda2} >
220 : !--------------------------------------------------------------------------------------
221 :
222 70860 : call cwtime(cpu, wall, gflops, "start")
223 :
224 : !Specific definitions
225 70860 : d3etot_t1_k=zero
226 264293592 : vlocal1=zero
227 70860 : dum_lambda=zero
228 70860 : berryopt=0;optlocal=0;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0
229 :
230 : !Initialize rf Hamiltonian (the k-dependent part is prepared in getgh1c_setup)
231 : call rf_hamkq%init(cplex,gs_hamkq,natom+1,&
232 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
233 70860 : & mpi_spintab=mpi_enreg%my_isppoltab)
234 :
235 70860 : call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.)
236 :
237 : !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
238 70860 : if (dtset%ffnl_lw==0) then
239 193840 : ABI_MALLOC(ffnlk,(npw_k,0,psps%lmnmax,psps%ntypat))
240 242300 : ABI_MALLOC(ffnl1,(npw_k,2,psps%lmnmax,psps%ntypat))
241 75830651 : ffnl1(:,1,:,:)=ffnl_k(:,1,:,:)
242 75830651 : ffnl1(:,2,:,:)=ffnl_k(:,1+i3dir,:,:)
243 : end if
244 : call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In
245 : kpt,kpt,i3dir,natom+1,natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In
246 : npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,part_ylmgr_k,& ! In
247 : dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,& ! Out
248 70860 : reuse_ffnlk=reuse_ffnlk, reuse_ffnl1=reuse_ffnl1) ! Optional
249 :
250 : !LOOP OVER BANDS
251 356236 : do iband=1,nband_k
252 :
253 285376 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
254 :
255 : !Select bks wf1
256 285376 : offset_cgi = (iband-1)*size_wf+icg
257 90262492 : cwavef1(:,:)= cg1(:,1+offset_cgi:size_wf+offset_cgi)
258 90262492 : cwavef2(:,:)= cg2(:,1+offset_cgi:size_wf+offset_cgi)
259 :
260 : !Compute < g |\partial_{gamma} H^{(0)} | u_{i,k}^{\lambda2} >
261 : call getgh1c(berryopt,cwavef2,dum_cwaveprj,gv1c,dum_grad_berry,&
262 : & dum_gs1,gs_hamkq,dum_gvnl1,i3dir,natom+1,(/dum_lambda/),mpi_enreg,1,optlocal,&
263 570752 : & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
264 :
265 : !Apply the dot product with the ket wf (take into account occupation here)
266 : ! < u_{i,k}^{\lambda1}} | \partial_{gamma} H^{(0)} | u_{i,k}^{lambda2}} >
267 : call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef1,gv1c, &
268 285376 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
269 :
270 285376 : d3etot_t1_k(1)=d3etot_t1_k(1)+occ_k(iband)*dotr
271 356236 : d3etot_t1_k(2)=d3etot_t1_k(2)+occ_k(iband)*doti
272 :
273 : end do !iband
274 :
275 : !Clean rf_hamiltonian
276 70860 : call rf_hamkq%free()
277 :
278 : !Deallocations
279 70860 : ABI_FREE(kpg_k)
280 70860 : ABI_FREE(kpg1_k)
281 70860 : ABI_FREE(dkinpw)
282 70860 : ABI_FREE(kinpw1)
283 70860 : ABI_FREE(ffnlk)
284 70860 : ABI_FREE(ffnl1)
285 70860 : ABI_FREE(ph3d)
286 :
287 70860 : call cwtime(cpu, wall, gflops, "stop")
288 :
289 : !------------------------------------T2------------------------------------------------
290 : !q-gradient of CB projector x rf Hamiltonian lambda 2:
291 : ! < u_{i,k}^{\lambda1}} | \partial_{gamma} Q_k H^{\lambda2} | u_{i,k}^{(0)} >
292 : !--------------------------------------------------------------------------------------
293 :
294 : !Create array for ddk 1wf from file
295 283440 : ABI_MALLOC(cg1_ddk,(2,size_wf,nband_k))
296 356236 : do iband=1,nband_k
297 285376 : call ddk_f%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=cg1_aux)
298 90333352 : cg1_ddk(:,:,iband)=cg1_aux(:,:)
299 : end do
300 :
301 70860 : call cwtime(cpu, wall, gflops, "start")
302 :
303 : !For \lambda1=\lambda2 T2 is inferred from the cc of T3
304 70860 : if (.not.samepert) then
305 :
306 : !Specific definitions
307 70088 : d3etot_t2_k=zero
308 :
309 : !LOOP OVER BANDS
310 352296 : do iband=1,nband_k
311 :
312 282208 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
313 :
314 : !Select bks wfs
315 282208 : offset_cgi = (iband-1)*size_wf+icg
316 89198692 : cwavef1(:,:)= cg1(:,1+offset_cgi:size_wf+offset_cgi)
317 :
318 : !LOOP OVER BANDS
319 1519528 : do jband=1,nband_k
320 :
321 : !Select ddk wf1
322 388853008 : cg1_aux(:,:)=cg1_ddk(:,:,jband)
323 :
324 : !Load < u_{j,k}^{(0) | H^{\lambda2}+V^{\lambda2}} | u_{i,k}^{(0)} >
325 1167232 : ii=2*jband-1+(iband-1)*2*nband_k
326 1167232 : cj_h1_ci(1)=eig2_k(ii)
327 1167232 : cj_h1_ci(2)=eig2_k(ii+1)
328 :
329 : !Calculate: < u_{i,k}^{lambda1}} | u_{j,k}^{k_{\gamma}} >
330 : call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef1,cg1_aux, &
331 1167232 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
332 :
333 : !Calculate the contribution to T2
334 1167232 : cprodr=dotr*cj_h1_ci(1)-doti*cj_h1_ci(2)
335 1167232 : cprodi=dotr*cj_h1_ci(2)+doti*cj_h1_ci(1)
336 1167232 : d3etot_t2_k(1)=d3etot_t2_k(1)-cprodr*occ_k(iband)
337 1449440 : d3etot_t2_k(2)=d3etot_t2_k(2)-cprodi*occ_k(iband)
338 :
339 : end do !jband
340 :
341 : end do !iband
342 :
343 : end if !samepert
344 :
345 70860 : call cwtime(cpu, wall, gflops, "stop")
346 :
347 : !------------------------------------T3------------------------------------------------
348 : !rf Hamiltonian lambda 1 x q-gradient of CB projector
349 : ! < u_{i,k}^{(0) | (H^{\lambda1})^{\dagger} \partial_{gamma} Q_k | u_{i,k}^{\lambda2}} >
350 : !--------------------------------------------------------------------------------------
351 :
352 70860 : call cwtime(cpu, wall, gflops, "start")
353 : !Specific definitions
354 70860 : d3etot_t3_k=zero
355 :
356 : !LOOP OVER BANDS
357 356236 : do iband=1,nband_k
358 :
359 285376 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
360 :
361 : !Select bks wfs
362 285376 : offset_cgi = (iband-1)*size_wf+icg
363 90262492 : cwavef2(:,:)= cg2(:,1+offset_cgi:size_wf+offset_cgi)
364 :
365 : !LOOP OVER BANDS
366 1538060 : do jband=1,nband_k
367 :
368 : !Select ddk wf1
369 396115408 : cg1_aux(:,:)=cg1_ddk(:,:,jband)
370 :
371 : !Load (< u_{j,k}^{(0) | H^{\lambda1}+V^{\lambda1}} | u_{i,k}^{(0)} >)^*
372 1181824 : ii=2*jband-1+(iband-1)*2*nband_k
373 1181824 : cj_h1_ci(1)=eig1_k(ii)
374 1181824 : cj_h1_ci(2)=-eig1_k(ii+1)
375 :
376 : !Calculate: < u_{j,k}^{k_{\gamma}} | u_{i,k}^{lambda2}} >
377 : call dotprod_g(dotr,doti,istwf_k,size_wf,2,cg1_aux,cwavef2, &
378 1181824 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
379 :
380 : !Calculate the contribution to T3
381 1181824 : cprodr=dotr*cj_h1_ci(1)-doti*cj_h1_ci(2)
382 1181824 : cprodi=dotr*cj_h1_ci(2)+doti*cj_h1_ci(1)
383 1181824 : d3etot_t3_k(1)=d3etot_t3_k(1)-cprodr*occ_k(iband)
384 1467200 : d3etot_t3_k(2)=d3etot_t3_k(2)-cprodi*occ_k(iband)
385 :
386 : end do !jband
387 :
388 : end do !iband
389 :
390 70860 : ABI_FREE(cg1_ddk)
391 70860 : ABI_FREE(cg1_aux)
392 70860 : ABI_FREE(vpsp1)
393 70860 : ABI_FREE(vlocal1)
394 :
395 70860 : if (samepert) then
396 772 : d3etot_t2_k(1)=d3etot_t3_k(1)
397 772 : d3etot_t2_k(2)=-d3etot_t3_k(2)
398 : end if
399 :
400 70860 : call cwtime(cpu, wall, gflops, "stop")
401 : !------------------------------------T4------------------------------------------------
402 : !q-gradient of rf Hamiltonian lambda 2
403 : ! < u_{i,k}^{\lambda1} | H^{\lambda2}_{gamma} | u_{i,k}^{(0)} >
404 : !--------------------------------------------------------------------------------------
405 :
406 : !For \lambda1=\lambda2 T4 is inferred from the cc of T5
407 70860 : if (.not.samepert) then
408 :
409 : !Specific definitions and allocations
410 373664 : d3etot_t4_k=zero
411 70088 : optlocal=1;optnl=1
412 70088 : dimffnlk=0
413 70088 : if (i2pert/=natom+2) then
414 420432 : ABI_MALLOC(vlocal1,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc))
415 280288 : ABI_MALLOC(vpsp1,(2*nfft))
416 140144 : ABI_MALLOC(gvloc1dqc,(2,size_wf))
417 140144 : ABI_MALLOC(gvnl1dqc,(2,size_wf))
418 : end if
419 : if (i2pert<=natom) fac=-one
420 70088 : if (i2pert==natom+2) fac=one
421 70088 : if (i2pert==natom+3.or.i2pert==natom+4) fac=-half
422 70088 : if (i2pert<=natom) then
423 7864 : nylmgrtmp=3
424 7864 : dimffnlk=1
425 7864 : dimffnl1=2
426 62224 : else if (i2pert==natom+3.or.i2pert==natom+4) then
427 62208 : nylmgrtmp=nylmgr
428 62208 : dimffnl1=10
429 62208 : ABI_FREE(part_ylmgr_k)
430 311040 : ABI_MALLOC(part_ylmgr_k,(npw_k,nylmgrtmp,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
431 167843032 : part_ylmgr_k(:,:,:)=ylmgr_k(:,:,:)
432 : end if
433 :
434 : !Do loop to compute both extradiagonal shear-strain components
435 171280 : do idq=1,n2dq
436 :
437 101192 : call cwtime(cpu, wall, gflops, "start")
438 :
439 101192 : if (i2pert/=natom+2) then
440 101176 : idir=i2dir; if (i2pert==natom+4) idir=idq*3+i2dir
441 : !Initialize rf Hamiltonian (the k-dependent part is prepared in getgh1c_setup)
442 : call rf_hamkq%init(2,gs_hamkq,i2pert,&
443 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
444 101176 : & mpi_spintab=mpi_enreg%my_isppoltab)
445 :
446 : !Set up local potentials with proper dimensioning
447 : !and load the spin-dependent part of the Hamiltonians
448 694024304 : vpsp1=vpsp1_i2pertdq(:,isppol,idq)
449 : call rf_transgrid_and_pack(isppol,nspden,usepaw,2,nfft,nfft,ngfft,&
450 101176 : & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1)
451 : call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,&
452 101176 : & with_nonlocal=with_nonlocal_i2pert)
453 :
454 : !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
455 101176 : if (dtset%ffnl_lw==0) then
456 343320 : ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,psps%ntypat))
457 18878035 : if (dimffnlk==1) ffnlk(:,1,:,:)=ffnl_k(:,1,:,:)
458 411984 : ABI_MALLOC(ffnl1,(npw_k,dimffnl1,psps%lmnmax,psps%ntypat))
459 68664 : if (dimffnl1==2) then
460 18815827 : ffnl1(:,1,:,:)=ffnl_k(:,1,:,:)
461 18815827 : ffnl1(:,2,:,:)=ffnl_k(:,1+i3dir,:,:)
462 : else
463 838935144 : ffnl1(:,1:dimffnl1,:,:)=ffnl_k(:,1:dimffnl1,:,:)
464 : end if
465 : end if
466 :
467 : call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,i2pert,i3dir, &
468 : & dtset%natom,rmet,rprimd,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgrtmp,useylmgr1,kg_k, &
469 : & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1, &
470 101176 : & reuse_ffnlk=reuse_ffnlk,reuse_ffnl1=reuse_ffnl1)
471 :
472 : end if
473 :
474 : !LOOP OVER BANDS
475 507816 : do iband=1,nband_k
476 :
477 406624 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
478 :
479 : !Select bks wfs
480 406624 : offset_cgi = (iband-1)*size_wf+icg
481 126193012 : cwavef1(:,:)= cg1(:,1+offset_cgi:size_wf+offset_cgi)
482 :
483 : !Perturbation-specific part
484 406624 : if (i2pert==natom+2) then
485 : if (samepert) then
486 : ! Read from d2_dkdk_f
487 : call d2_dkdk_f%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=gv1c)
488 : else
489 : ! Read from d2_dkdk_f2
490 384 : call d2_dkdk_f2%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=gv1c)
491 : end if
492 : else
493 125591572 : cwave0i(:,:)= cg(:,1+offset_cgi:size_wf+offset_cgi)
494 :
495 : !Compute < g |H^{\lambda2}}_{\gamma} | u_{i,k}^{(0)} >
496 : call getgh1dqc(cwave0i,dum_cwaveprj,gv1c,gvloc1dqc,gvnl1dqc,gs_hamkq, &
497 406240 : & idir,i2pert,mpi_enreg,optlocal,optnl,i3dir,rf_hamkq)
498 : end if
499 :
500 : !Calculate: < u_{j,k}^{\lambda1} | |H^{\lambda2}}_{\gamma} | u_{i,k}^{(0)} >
501 : call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef1,gv1c, &
502 406624 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
503 :
504 : !Calculate the contribution to T4
505 406624 : d3etot_t4_k(1,idq)=d3etot_t4_k(1,idq)+dotr*occ_k(iband)
506 507816 : d3etot_t4_k(2,idq)=d3etot_t4_k(2,idq)+doti*occ_k(iband)
507 :
508 : end do !iband
509 :
510 101192 : if (i2pert/=natom+2) then
511 :
512 : !Clean the rf_hamiltonian
513 101176 : call rf_hamkq%free()
514 :
515 : !Deallocations
516 101176 : ABI_FREE(kpg_k)
517 101176 : ABI_FREE(kpg1_k)
518 101176 : ABI_FREE(dkinpw)
519 101176 : ABI_FREE(kinpw1)
520 101176 : ABI_FREE(ffnlk)
521 101176 : ABI_FREE(ffnl1)
522 101176 : ABI_FREE(ph3d)
523 :
524 : end if
525 :
526 101192 : call cwtime(cpu, wall, gflops, "stop")
527 :
528 : !Apply the perturbation-dependent prefactors on T4
529 101192 : tmpre=d3etot_t4_k(1,idq); tmpim=d3etot_t4_k(2,idq)
530 101192 : if (i2pert<=natom.or.i2pert==natom+2) then
531 7880 : d3etot_t4_k(1,idq)=-tmpim
532 7880 : d3etot_t4_k(2,idq)=tmpre
533 : end if
534 576048 : d3etot_t4_k(:,idq)=d3etot_t4_k(:,idq)*fac
535 :
536 : end do !idq
537 :
538 70088 : if (i2pert/=natom+2) then
539 70072 : ABI_FREE(gvloc1dqc)
540 70072 : ABI_FREE(gvnl1dqc)
541 70072 : ABI_FREE(vlocal1)
542 70072 : ABI_FREE(vpsp1)
543 : end if
544 :
545 : end if !samepert
546 :
547 : !------------------------------------T5------------------------------------------------
548 : !q-gradient of rf Hamiltonian lambda 1
549 : ! < u_{i,k}^{(0)} | (H^{\lambda1}_{gamma})^{\dagger} | u_{i,k}^{\lambda2} >
550 : !--------------------------------------------------------------------------------------
551 :
552 : !Specific definitions and allocations
553 283440 : d3etot_t5_k=zero
554 70860 : optlocal=1;optnl=1
555 70860 : dimffnlk=0
556 70860 : if (i1pert/=natom+2) then
557 271872 : ABI_MALLOC(vlocal1,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc))
558 181248 : ABI_MALLOC(vpsp1,(2*nfft))
559 90624 : ABI_MALLOC(gvloc1dqc,(2,size_wf))
560 90624 : ABI_MALLOC(gvnl1dqc,(2,size_wf))
561 : end if
562 70860 : if (i1pert<=natom) fac=-one
563 70860 : if (i1pert==natom+2) fac=one
564 70860 : if (i1pert==natom+3.or.i1pert==natom+4) fac=-half
565 70860 : if (i1pert<=natom) then
566 45312 : nylmgrtmp=3
567 45312 : dimffnlk=1
568 45312 : dimffnl1=2
569 45312 : ABI_FREE(part_ylmgr_k)
570 181248 : ABI_MALLOC(part_ylmgr_k,(npw_k,nylmgrtmp,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
571 40872840 : part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3,:)
572 25548 : else if (i1pert==natom+3.or.i1pert==natom+4) then
573 0 : nylmgrtmp=nylmgr
574 0 : dimffnl1=10
575 0 : ABI_FREE(part_ylmgr_k)
576 0 : ABI_MALLOC(part_ylmgr_k,(npw_k,nylmgrtmp,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
577 25548 : part_ylmgr_k(:,:,:)=ylmgr_k(:,:,:)
578 : end if
579 :
580 : !Do loop to compute both extradiagonal shear-strain components
581 141720 : do idq=1,n1dq
582 :
583 70860 : call cwtime(cpu, wall, gflops, "start")
584 70860 : if (i1pert/=natom+2) then
585 45312 : idir=i1dir; if (i1pert==natom+4) idir=idq*3+i1dir
586 : !Initialize rf Hamiltonian (the k-dependent part is prepared in getgh1c_setup)
587 : call rf_hamkq%init(2,gs_hamkq,i1pert,&
588 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
589 45312 : & mpi_spintab=mpi_enreg%my_isppoltab)
590 :
591 : !Set up local potentials with proper dimensioning
592 : !and load the spin-dependent part of the Hamiltonians
593 305946624 : vpsp1=vpsp1_i1pertdq(:,isppol,idq)
594 : call rf_transgrid_and_pack(isppol,nspden,usepaw,2,nfft,nfft,ngfft,&
595 45312 : & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1)
596 : call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,&
597 45312 : & with_nonlocal=with_nonlocal_i1pert)
598 :
599 : !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
600 45312 : if (dtset%ffnl_lw==0) then
601 151040 : ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,psps%ntypat))
602 40752008 : if (dimffnlk==1) ffnlk(:,1,:,:)=ffnl_k(:,1,:,:)
603 181248 : ABI_MALLOC(ffnl1,(npw_k,dimffnl1,psps%lmnmax,psps%ntypat))
604 30208 : if (dimffnl1==2) then
605 40752008 : ffnl1(:,1,:,:)=ffnl_k(:,1,:,:)
606 40752008 : ffnl1(:,2,:,:)=ffnl_k(:,1+i3dir,:,:)
607 : else
608 0 : ffnl1(:,1:dimffnl1,:,:)=ffnl_k(:,1:dimffnl1,:,:)
609 : end if
610 : end if
611 : call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,i1pert,i3dir, &
612 : & dtset%natom,rmet,rprimd,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgrtmp,useylmgr1,kg_k, &
613 : & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1, &
614 45312 : & reuse_ffnlk=reuse_ffnlk,reuse_ffnl1=reuse_ffnl1)
615 :
616 : end if
617 :
618 : !LOOP OVER BANDS
619 356236 : do iband=1,nband_k
620 :
621 285376 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
622 :
623 : !Select bks wfs
624 285376 : offset_cgi = (iband-1)*size_wf+icg
625 90262492 : cwavef2(:,:)= cg2(:,1+offset_cgi:size_wf+offset_cgi)
626 :
627 : !Perturbation-specific part
628 285376 : if (i1pert==natom+2) then
629 104128 : call d2_dkdk_f%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=gv1c)
630 : else
631 53892960 : cwave0i(:,:)= cg(:,1+offset_cgi:size_wf+offset_cgi)
632 :
633 : !Compute < g |H^{\lambda1}}_{\gamma} | u_{i,k}^{(0)} >
634 : call getgh1dqc(cwave0i,dum_cwaveprj,gv1c,gvloc1dqc,gvnl1dqc,gs_hamkq, &
635 181248 : & idir,i1pert,mpi_enreg,optlocal,optnl,i3dir,rf_hamkq)
636 : end if
637 :
638 : !Calculate: < u_{j,k}^{\lambda2} | |H^{\lambda1}}_{\gamma} | u_{i,k}^{(0)} >
639 : call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef2,gv1c, &
640 285376 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
641 :
642 : !Calculate the contribution to T5:
643 285376 : d3etot_t5_k(1,idq)=d3etot_t5_k(1,idq)+dotr*occ_k(iband)
644 356236 : d3etot_t5_k(2,idq)=d3etot_t5_k(2,idq)+doti*occ_k(iband)
645 :
646 : end do !iband
647 :
648 70860 : if (i1pert/=natom+2) then
649 :
650 : !Clean the rf_hamiltonian
651 45312 : call rf_hamkq%free()
652 :
653 : !Deallocations
654 45312 : ABI_FREE(kpg_k)
655 45312 : ABI_FREE(kpg1_k)
656 45312 : ABI_FREE(dkinpw)
657 45312 : ABI_FREE(kinpw1)
658 45312 : ABI_FREE(ffnlk)
659 45312 : ABI_FREE(ffnl1)
660 45312 : ABI_FREE(ph3d)
661 :
662 : end if
663 :
664 70860 : call cwtime(cpu, wall, gflops, "stop")
665 :
666 : !Apply the perturbation-dependent prefactors on T5
667 70860 : tmpre=d3etot_t5_k(1,idq); tmpim=d3etot_t5_k(2,idq)
668 70860 : if (i1pert<=natom.or.i1pert==natom+2) then
669 70860 : d3etot_t5_k(1,idq)=-tmpim
670 70860 : d3etot_t5_k(2,idq)=tmpre
671 : end if
672 212580 : d3etot_t5_k(:,idq)=d3etot_t5_k(:,idq)*fac
673 :
674 : !Apply now the conjugate complex:
675 : !(< u_{j,k}^{\lambda2} | |H^{\lambda1}}_{\gamma} | u_{i,k}^{(0)} >)*
676 70860 : tmpim=d3etot_t5_k(2,idq)
677 283440 : d3etot_t5_k(2,idq)=-tmpim
678 :
679 : end do !idq
680 :
681 :
682 70860 : if (i1pert/=natom+2) then
683 45312 : ABI_FREE(gvloc1dqc)
684 45312 : ABI_FREE(gvnl1dqc)
685 45312 : ABI_FREE(vlocal1)
686 45312 : ABI_FREE(vpsp1)
687 : end if
688 :
689 70860 : if (samepert) then
690 1544 : d3etot_t4_k(1,:)=d3etot_t5_k(1,:)
691 1544 : d3etot_t4_k(2,:)=-d3etot_t5_k(2,:)
692 : end if
693 :
694 : !Scale d3etot_k contributions by the kpt weight
695 212580 : d3etot_t1_k(:)=d3etot_t1_k(:)*wtk_k
696 212580 : d3etot_t2_k(:)=d3etot_t2_k(:)*wtk_k
697 212580 : d3etot_t3_k(:)=d3etot_t3_k(:)*wtk_k
698 376752 : d3etot_t4_k(:,:)=d3etot_t4_k(:,:)*wtk_k
699 283440 : d3etot_t5_k(:,:)=d3etot_t5_k(:,:)*wtk_k
700 :
701 : !Deallocations
702 70860 : ABI_FREE(cwave0i)
703 70860 : ABI_FREE(cwave0j)
704 70860 : ABI_FREE(cwavef1)
705 70860 : ABI_FREE(cwavef2)
706 70860 : ABI_FREE(gv1c)
707 70860 : ABI_FREE(dum_vpsp)
708 70860 : ABI_FREE(dum_vlocal)
709 70860 : ABI_FREE(dum_cwaveprj)
710 70860 : ABI_FREE(part_ylmgr_k)
711 :
712 :
713 : DBG_EXIT("COLL")
714 :
715 425160 : end subroutine dfpt_1wf
716 : !!***
717 :
718 : end module m_dfptlw_wf
719 : !!***
|
