Line data Source code
1 : !{\src2tex{textfont=tt}}
2 : !!****m* ABINIT/m_dfptlw_nv
3 : !! NAME
4 : !! m_dfptlw_nv
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_nv
24 :
25 : use defs_basis
26 : use defs_abitypes
27 : use defs_datatypes
28 : use m_abicore
29 : use m_xmpi
30 : use m_errors
31 : use m_mpinfo
32 : use m_dtset
33 : use m_hamiltonian
34 : use m_cgtools
35 : use m_wfk
36 : use m_xmpi
37 : use m_getgh1c
38 : use m_mklocl
39 : use m_pawcprj
40 : use m_pawfgr
41 :
42 : use m_dfpt_elt, only : dfpt_ewalddq, dfpt_ewalddqdq
43 : use m_kg, only : mkkpg
44 : use m_dynmat, only : cart39
45 :
46 : implicit none
47 :
48 : private
49 : !!***
50 :
51 : public :: dfptlw_nv
52 : public :: dfptlw_geom
53 : !!***
54 :
55 : ! *************************************************************************
56 :
57 : contains
58 : !!***
59 :
60 : !!****f* ABINIT/m_dfptlw_nv/dfptlw_nv
61 : !! NAME
62 : !! dfptlw_nv
63 : !!
64 : !! FUNCTION
65 : !! This routine calculates the nonvariational Ewald contributions to the
66 : !! spatial-dispersion third-order energy derivatives.
67 : !!
68 : !! INPUTS
69 : !! dtset <type(dataset_type)>=all input variables for this dataset
70 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
71 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
72 : !! mpert=maximum number of ipert
73 : !! my_natom=number of atoms treated by current processor
74 : !! rmet(3,3)=metric tensor in real space (length units squared)
75 : !! rprimd(3,3)=dimensional primitive translations (bohr)
76 : !! rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
77 : !! that have to be computed
78 : !! 1 -> element has to be computed explicitely
79 : !! -1 -> use symmetry operations to obtain the corresponding element
80 : !! ucvol=unit cell volume in (whatever length scale units)**3
81 : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
82 : !! zion(ntypat)=charge on each type of atom (real number)
83 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
84 : !! comm_atom=--optional-- MPI communicator over atoms
85 : !!
86 : !! OUTPUT
87 : !! d3etot_nv(2,3,mpert,3,mpert,3,mpert)= array with the nonvariational
88 : !! contributions of d3etot
89 : !!
90 : !! SOURCE
91 :
92 6 : subroutine dfptlw_nv(d3etot_nv,dtset,gmet,gprimd,mpert,my_natom,rfpert,rmet,rprimd,ucvol,xred,zion, &
93 3 : & mpi_atmtab,comm_atom ) ! optional arguments (parallelism))
94 :
95 : !Arguments ------------------------------------
96 : !scalars
97 : integer , intent(in) :: mpert,my_natom
98 : real(dp) :: ucvol
99 : type(dataset_type),intent(in) :: dtset
100 : integer,optional,intent(in) :: comm_atom
101 :
102 : !arrays
103 : integer,optional,target,intent(in) :: mpi_atmtab(:)
104 : integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert)
105 : real(dp), intent(inout) :: d3etot_nv(2,3,mpert,3,mpert,3,mpert)
106 : real(dp), intent(in) :: gmet(3,3),rmet(3,3),xred(3,dtset%natom),zion(*)
107 : real(dp), intent(in) :: gprimd(3,3),rprimd(3,3)
108 :
109 : !Local variables-------------------------------
110 : !scalars
111 : integer :: alpha,beta,delta,gamma,i1dir,i2dir,i3dir,ii,i1pert,i2pert,i3pert,istr,natom,sumg0
112 : real(dp) :: fac,tmpim,tmpre
113 : character(len=500) :: msg
114 :
115 : !arrays
116 : integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/)
117 : integer :: flg1(3),flg2(3)
118 3 : real(dp),allocatable :: dyewdq(:,:,:,:,:,:),dyewdqdq(:,:,:,:,:,:)
119 3 : real(dp),allocatable :: dyewdqdq_tII(:,:,:,:,:,:)
120 : real(dp) :: qphon(3),vec1(3),vec2(3)
121 : real(dp) :: rprimd_t(3,3),gprimd_t(3,3)
122 :
123 : ! *************************************************************************
124 :
125 : DBG_ENTER("COLL")
126 :
127 : !Initialiations
128 3 : natom=dtset%natom
129 :
130 3 : if (dtset%lw_flexo==1.or.dtset%lw_flexo==3) then
131 :
132 : !1st q-gradient of Ewald contribution to the IFCs
133 18 : ABI_MALLOC(dyewdq,(2,3,natom,3,natom,3))
134 3 : sumg0=0;qphon(:)=zero
135 : call dfpt_ewalddq(dyewdq,gmet,my_natom,natom,qphon,rmet,sumg0,dtset%typat,ucvol,xred,zion,&
136 6 : & mpi_atmtab=mpi_atmtab,comm_atom=comm_atom)
137 :
138 3 : i3pert=natom+8
139 9 : do i1pert=1,natom
140 27 : do i1dir=1,3
141 60 : do i2pert=1,natom
142 162 : do i2dir=1,3
143 468 : do i3dir=1,3
144 324 : tmpre=dyewdq(1,i1dir,i1pert,i2dir,i2pert,i3dir)
145 324 : tmpim=dyewdq(2,i1dir,i1pert,i2dir,i2pert,i3dir)
146 324 : if (abs(tmpre)>=tol8) d3etot_nv(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)= tmpre
147 432 : if (abs(tmpim)>=tol8) d3etot_nv(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)= tmpim
148 : end do
149 : end do
150 : end do
151 : end do
152 : end do
153 3 : ABI_FREE(dyewdq)
154 :
155 : end if
156 :
157 3 : if (dtset%lw_flexo==1.or.dtset%lw_flexo==4) then
158 :
159 : !2nd q-gradient of Ewald contribution to the IFCs
160 18 : ABI_MALLOC(dyewdqdq,(2,3,natom,3,3,3))
161 6 : ABI_MALLOC(dyewdqdq_tII,(2,3,natom,3,3,3))
162 3 : sumg0=1;qphon(:)=zero
163 : call dfpt_ewalddqdq(dyewdqdq,gmet,my_natom,natom,qphon,rmet,sumg0,dtset%typat,ucvol,xred,zion,&
164 6 : & mpi_atmtab=mpi_atmtab,comm_atom=comm_atom)
165 :
166 :
167 : !Convert the indexes labelling the strain perturbation into cartesian coordinates
168 : !Transform the metric perturbation direction
169 : !(treat it as an atomic displacement)
170 12 : flg1(:)=1
171 9 : do i1pert=1,natom
172 27 : do i1dir=1,3
173 78 : do gamma=1,3
174 180 : do ii=1,2
175 486 : do delta=1,3
176 1296 : do beta=1,3
177 1296 : vec1(beta)=dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma)
178 : end do
179 324 : call cart39(flg1,flg2,gprimd,i1pert,natom,rprimd,vec1,vec2)
180 1404 : do beta=1,3
181 1296 : dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma)=vec2(beta)
182 : end do
183 : end do
184 : end do
185 : end do
186 : end do
187 : end do
188 :
189 : !Transform the second q-gradient direction
190 : !(treat it as an electric field)
191 9 : do i1pert=1,natom
192 27 : do i1dir=1,3
193 78 : do gamma=1,3
194 180 : do ii=1,2
195 486 : do beta=1,3
196 1296 : do delta=1,3
197 1296 : vec1(delta)=dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma)
198 : end do
199 324 : call cart39(flg1,flg2,gprimd,natom+2,natom,rprimd,vec1,vec2)
200 1404 : do delta=1,3
201 1296 : dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma)=vec2(delta)
202 : end do
203 : end do
204 : end do
205 : end do
206 : end do
207 : end do
208 :
209 : !Transform the first q-gradient direction
210 : !(treat it as an electric field)
211 9 : do i1pert=1,natom
212 27 : do i1dir=1,3
213 60 : do ii=1,2
214 162 : do beta=1,3
215 468 : do delta=1,3
216 1296 : do gamma=1,3
217 1296 : vec1(gamma)=dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma)
218 : end do
219 324 : call cart39(flg1,flg2,gprimd,natom+2,natom,rprimd,vec1,vec2)
220 1404 : do gamma=1,3
221 1296 : dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma)=vec2(gamma)
222 : end do
223 : end do
224 : end do
225 : end do
226 : end do
227 : end do
228 :
229 : !Convert to a type-II quantity
230 1740 : dyewdqdq_tII(:,:,:,:,:,:)=zero
231 9 : do i1pert=1,natom
232 27 : do alpha=1,3
233 78 : do gamma=1,3
234 234 : do beta=1,3
235 702 : do delta=1,3
236 : dyewdqdq_tII(:,alpha,i1pert,gamma,beta,delta)= &
237 : & dyewdqdq(:,alpha,i1pert,beta,delta,gamma) + &
238 : & dyewdqdq(:,alpha,i1pert,delta,gamma,beta) - &
239 1620 : & dyewdqdq(:,alpha,i1pert,gamma,beta,delta)
240 : end do
241 : end do
242 : end do
243 : end do
244 : end do
245 3 : ABI_FREE(dyewdqdq)
246 :
247 : !Transform back the first q-gradient direction to reduced coordinates
248 : !(treat it as an electric field)
249 39 : rprimd_t=transpose(rprimd)
250 39 : gprimd_t=transpose(gprimd)
251 3 : fac=two_pi**2
252 9 : do i1pert=1,natom
253 27 : do alpha=1,3
254 60 : do ii=1,2
255 162 : do delta=1,3
256 468 : do beta=1,3
257 1296 : do gamma=1,3
258 1296 : vec1(gamma)=dyewdqdq_tII(ii,alpha,i1pert,gamma,beta,delta)
259 : end do
260 324 : call cart39(flg1,flg2,rprimd_t,natom+2,natom,gprimd_t,vec1,vec2)
261 1404 : do gamma=1,3
262 1296 : dyewdqdq_tII(ii,alpha,i1pert,gamma,beta,delta)=vec2(gamma)*fac
263 : end do
264 : end do
265 : end do
266 : end do
267 : end do
268 : end do
269 :
270 3 : i3pert=natom+8
271 9 : do i1pert=1,natom
272 27 : do i1dir=1,3
273 60 : do i2pert=natom+3,natom+4
274 162 : do i2dir=1,3
275 108 : istr=(i2pert-natom-3)*3+i2dir
276 108 : beta=idx(2*istr-1); delta=idx(2*istr)
277 468 : do i3dir=1,3
278 324 : tmpre=dyewdqdq_tII(1,i1dir,i1pert,i3dir,beta,delta)
279 324 : tmpim=dyewdqdq_tII(2,i1dir,i1pert,i3dir,beta,delta)
280 324 : if (abs(tmpre)>=tol8) d3etot_nv(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)= half*tmpre
281 432 : if (abs(tmpim)>=tol8) d3etot_nv(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)= half*tmpim
282 : end do
283 : end do
284 : end do
285 : end do
286 : end do
287 3 : ABI_FREE(dyewdqdq_tII)
288 :
289 : end if
290 :
291 : !Print results
292 3 : if (dtset%prtvol>=10) then
293 0 : write(msg,'(3a)') ch10,' LONGWAVE NONVARIATIONAL EWALD D3ETOT: ',ch10
294 0 : call wrtout(std_out,msg,'COLL')
295 0 : call wrtout(ab_out,msg,'COLL')
296 0 : do i1pert=1,mpert
297 0 : do i1dir=1,3
298 0 : do i2pert=1,mpert
299 0 : do i2dir=1,3
300 0 : do i3pert=1,mpert
301 0 : do i3dir=1,3
302 0 : if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
303 0 : tmpre=d3etot_nv(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
304 0 : tmpim=d3etot_nv(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
305 0 : if (abs(tmpre)>zero.or.abs(tmpim)>zero) then
306 : write(msg,'(3(a,i2,a,i1),2f18.8)') &
307 0 : ' perts : ',i1pert,'.',i1dir,' / ',i2pert,'.',i2dir,' / ',i3pert,'.',i3dir,&
308 0 : & tmpre, tmpim
309 0 : call wrtout(std_out,msg,'COLL')
310 0 : call wrtout(ab_out,msg,'COLL')
311 : end if
312 : end if
313 : end do
314 : end do
315 : end do
316 : end do
317 : end do
318 : end do
319 0 : write(msg,'(a)') ch10
320 0 : call wrtout(std_out,msg,'COLL')
321 0 : call wrtout(ab_out,msg,'COLL')
322 : end if
323 :
324 : DBG_EXIT("COLL")
325 :
326 3 : end subroutine dfptlw_nv
327 : !!***
328 :
329 : !!****f* ABINIT/dfptlw_geom
330 : !! NAME
331 : !! dfptlw_geom
332 : !!
333 : !! FUNCTION
334 : !! This routine computes the nonvariational geometric contribution to the
335 : !! third-order energy derivative of the flexoelectric force-response tensor.
336 : !!
337 : !! INPUTS
338 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k
339 : !! cplex: if 1, several magnitudes are REAL, if 2, COMPLEX
340 : !! dimffnl= third dimension of ffnl_k
341 : !! dtset <type(dataset_type)>=all input variables for this dataset
342 : !! ffnl_k(dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives for this k point
343 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
344 : !! icg=shift to be applied on the location of data in the array cg
345 : !! i1dir,i2dir,i3dir=directions of the corresponding perturbations
346 : !! i1pert,i2pert = type of perturbation that has to be computed
347 : !! ikpt=number of the k-point
348 : !! isppol=1 for unpolarized, 2 for spin-polarized
349 : !! istwf_k=parameter that describes the storage of wfs
350 : !! kg_k(3,npw_k)=reduced planewave coordinates.
351 : !! kpt(3)=reduced coordinates of k point
352 : !! natom= number of atoms in the cell
353 : !! mkmem =number of k points treated by this node
354 : !! mpi_enreg=information about MPI parallelization
355 : !! mpw=maximum dimensioned size of npw or wfs at k
356 : !! natpert=number of atomic displacement perturbations
357 : !! nband_k=number of bands at this k point for that spin polarization
358 : !! n2dq= second dimension of d3etot_tgeom_k
359 : !! nfft=(effective) number of FFT grid points (for this proc)
360 : !! ngfft(1:18)=integer array with FFT box dimensions and other
361 : !! npw_k=number of plane waves at this k point
362 : !! nspden=number of spin-density components
363 : !! nsppol=1 for unpolarized, 2 for spin-polarized
364 : !! nylmgr=second dimension of ylmgr_k
365 : !! occ_k(nband_k)=occupation number for each band (usually 2) for each k.
366 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
367 : !! rmet(3,3)=real space metric (bohr**2)
368 : !! rprimd(3,3) = dimensional primitive translations (bohr)
369 : !! vpsp1_i1pertdqdq(cplex*nfft,nspden,n2dq)= local potential of second-order
370 : !! gradient Hamiltonian for i1pert
371 : !! vpsp1_i1pertdq_geom(cplex*nfft,nspden,3)= local potential of first-order
372 : !! gradient Hamiltonian for i1pert wrt i3dir and i2dir
373 : !! useylmgr= if 1 use the derivative of spherical harmonics
374 : !! wtk_k=weight assigned to the k point.
375 : !! ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics for the k point
376 : !! ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical
377 : !! harmonics for the k point
378 : !!
379 : !! OUTPUT
380 : !! d3etot_tgeom_k(2,n2dq)= nonvariational geometric contribution to d3etot for
381 : ! this kpt.
382 : !!
383 : !! SOURCE
384 :
385 41472 : subroutine dfptlw_geom(cg,d3etot_tgeom_k,dimffnl,dtset,ffnl_k, &
386 : & gs_hamkq,icg, &
387 : & i1dir,i2dir,i3dir,i1pert,i2pert,ikpt, &
388 41472 : & isppol,istwf_k,kg_k,kpt,mkmem,mpi_enreg,natom,mpw,nband_k,n2dq,nfft, &
389 41472 : & ngfft,npw_k,nspden,nsppol,nylmgr,occ_k, &
390 41472 : & psps,rmet,rprimd,useylmgr,vpsp1_i1pertdqdq,vpsp1_i1pertdq_geom,wtk_k,ylm_k,ylmgr_k)
391 :
392 : !Arguments ------------------------------------
393 : !scalars
394 : integer,intent(in) :: dimffnl,icg,ikpt,isppol,istwf_k
395 : integer,intent(in) :: i1dir,i1pert,i2dir,i2pert,i3dir
396 : integer,intent(in) :: natom,mkmem,mpw,nband_k,nfft
397 : integer,intent(in) :: npw_k,n2dq,nspden,nsppol,nylmgr
398 : integer,intent(in) :: useylmgr
399 : real(dp),intent(in) :: wtk_k
400 : type(dataset_type),intent(in) :: dtset
401 : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
402 : type(MPI_type),intent(in) :: mpi_enreg
403 : type(pseudopotential_type),intent(in) :: psps
404 :
405 : !arrays
406 : integer,intent(in) :: kg_k(3,npw_k),ngfft(18)
407 : real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol)
408 : real(dp),intent(in) :: ffnl_k(npw_k,dimffnl,psps%lmnmax,psps%ntypat)
409 : real(dp),intent(in) :: kpt(3),occ_k(nband_k)
410 : real(dp),intent(in) :: rmet(3,3),rprimd(3,3)
411 : real(dp),intent(in) :: vpsp1_i1pertdqdq(2*nfft,nspden,n2dq)
412 : real(dp),intent(in) :: vpsp1_i1pertdq_geom(2*nfft,nspden,3)
413 : real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)
414 : real(dp),intent(in) :: ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)
415 : real(dp),intent(out) :: d3etot_tgeom_k(2,n2dq)
416 :
417 : !Local variables-------------------------------
418 : !scalars
419 : integer :: beta,delta,dimffnlk,dimffnl1,gamma,iband,idq,ii,ipw,istr,nkpg,nkpg1,nylmgrpart
420 : integer :: optlocal,optnl,q1dir,q2dir,reuse_ffnlk,reuse_ffnl1,tim_getgh1c,useylmgr1
421 : real(dp) :: doti,dotr
422 : type(pawfgr_type) :: pawfgr
423 41472 : type(rf_hamiltonian_type) :: rf_hamkq
424 :
425 : !arrays
426 : integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/)
427 : real(dp) :: q1dirs(2),q2dirs(2)
428 41472 : real(dp),allocatable :: cwave0i(:,:)
429 41472 : real(dp),allocatable :: dkinpw(:)
430 41472 : real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:)
431 41472 : real(dp),allocatable :: gh1dqc(:,:),gh1dqpkc(:,:),gvloc1dqc(:,:),gvnl1dqc(:,:)
432 41472 : real(dp),allocatable :: kinpw1(:),kpg_k(:,:),kpg1_k(:,:),kpg_pk(:,:),ph3d(:,:,:),ph3d1(:,:,:)
433 41472 : real(dp),allocatable :: dum_vlocal(:,:,:,:),vlocal1dq(:,:,:,:), dum_vpsp(:)
434 41472 : real(dp),allocatable :: vpsp1dq(:),part_ylmgr_k(:,:,:)
435 41472 : type(pawcprj_type),allocatable :: dum_cwaveprj(:,:)
436 :
437 : ! *************************************************************************
438 :
439 : DBG_ENTER("COLL")
440 :
441 : !Definitions
442 41472 : tim_getgh1c=0
443 41472 : useylmgr1=useylmgr;optlocal=1;optnl=1
444 41472 : nylmgrpart=3
445 41472 : nkpg=3
446 228096 : d3etot_tgeom_k(:,:)=zero
447 41472 : reuse_ffnlk=1 ; if (dtset%ffnl_lw==1) reuse_ffnlk=0
448 41472 : reuse_ffnl1=1 ; if (dtset%ffnl_lw==1) reuse_ffnl1=0
449 :
450 : !Allocations
451 124416 : ABI_MALLOC(cwave0i,(2,npw_k*dtset%nspinor))
452 124416 : ABI_MALLOC(dum_vpsp,(nfft))
453 248832 : ABI_MALLOC(dum_vlocal,(ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc))
454 41472 : ABI_MALLOC(dum_cwaveprj,(0,0))
455 165888 : ABI_MALLOC(vpsp1dq,(2*nfft))
456 248832 : ABI_MALLOC(vlocal1dq,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc))
457 82944 : ABI_MALLOC(gh1dqc,(2,npw_k*dtset%nspinor))
458 82944 : ABI_MALLOC(gvloc1dqc,(2,npw_k*dtset%nspinor))
459 82944 : ABI_MALLOC(gvnl1dqc,(2,npw_k*dtset%nspinor))
460 165888 : ABI_MALLOC(part_ylmgr_k,(npw_k,3, psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
461 37409040 : part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3,:)
462 82944 : ABI_MALLOC(gh1dqpkc,(2,npw_k*dtset%nspinor))
463 165888 : ABI_MALLOC(kpg_pk,(npw_k,nkpg))
464 :
465 : !Generate k+G vectors
466 41472 : call mkkpg(kg_k,kpg_pk,kpt,nkpg,npw_k)
467 :
468 : !Since this is a type-I term, it has to be done for both up and down
469 : !extradiagonal shear strains
470 41472 : gamma=i3dir
471 103680 : do idq=1, n2dq
472 62208 : if (i2pert==natom+3) then
473 20736 : istr=i2dir
474 : else
475 41472 : istr=idq*3+i2dir
476 : endif
477 62208 : beta=idx(2*istr-1); delta=idx(2*istr)
478 :
479 : !-----------------------------------------------------------------------------------------------
480 : ! q1-gradient of atomic displacement 1st order hamiltonian:
481 : ! < u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}_{\q1dir} \delta_{\beta\q2dir}| u_{i,k}^{(0)} >
482 : !-----------------------------------------------------------------------------------------------
483 62208 : dimffnlk=1
484 62208 : dimffnl1=2
485 186624 : q1dirs=(/gamma,delta/)
486 186624 : q2dirs=(/delta,gamma/)
487 186624 : do ii=1,2
488 124416 : q1dir=q1dirs(ii)
489 124416 : q2dir=q2dirs(ii)
490 :
491 186624 : if (beta==q2dir) then
492 :
493 : !Get q-gradient of first-order local part of the pseudopotential
494 : ! call dfpt_vlocaldq(atindx,2,gs_hamkq%gmet,gsqcut,i1dir,i1pert,mpi_enreg, &
495 : ! & psps%mqgrid_vl,dtset%natom,&
496 : ! & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), &
497 : ! & ph1d,q1dir,psps%qgrid_vl,&
498 : ! & dtset%qptn,ucvol,psps%vlspl,vpsp1dq)
499 : ! write(300,*) vpsp1dq(:)-vpsp1_i1pertdq_geom(:,isppol,q1dir)
500 :
501 :
502 : !Set up q-gradient of local potential vlocal1dq with proper dimensioning
503 279977472 : vpsp1dq(:)=vpsp1_i1pertdq_geom(:,isppol,q1dir)
504 : call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,nfft,ngfft,&
505 41472 : & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dq,dum_vlocal,vlocal1dq)
506 :
507 : !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup)
508 : call rf_hamkq%init(2,gs_hamkq,i1pert,&
509 41472 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
510 41472 : call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dq,with_nonlocal=.true.)
511 :
512 : !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
513 41472 : if (dtset%ffnl_lw==0) then
514 138240 : ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,psps%ntypat))
515 37298448 : ffnlk(:,1,:,:)=ffnl_k(:,1,:,:)
516 138240 : ABI_MALLOC(ffnl1,(npw_k,dimffnl1,psps%lmnmax,psps%ntypat))
517 37298448 : ffnl1(:,1,:,:)=ffnl_k(:,1,:,:)
518 37298448 : ffnl1(:,2,:,:)=ffnl_k(:,1+q1dir,:,:)
519 : end if
520 : call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,i1dir,i1pert,q1dir, &
521 : & dtset%natom,rmet,rprimd,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgrpart,useylmgr1,kg_k, &
522 : & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,&
523 41472 : & ph3d,ph3d1,reuse_ffnlk=reuse_ffnlk,reuse_ffnl1=reuse_ffnl1)
524 :
525 : !LOOP OVER BANDS
526 207360 : do iband=1,nband_k
527 :
528 165888 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
529 :
530 : !Read ket ground-state wavefunctions
531 49325760 : cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg)
532 :
533 : !Compute < g |H^{\tau_{\kappa\alpha}}_{\q1dir} | u_{i,k}^{(0)} >
534 : call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq, &
535 165888 : & i1dir,i1pert,mpi_enreg,optlocal,optnl,q1dir,rf_hamkq)
536 :
537 : !Calculate:
538 : !<u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}}_{\q1dir} | u_{i,k}^{(0)} >
539 : call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0i,gh1dqc, &
540 165888 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
541 :
542 : !Take into account the two pi factor from the term
543 : !(\hat{p}_{k\beta + \frac{q_{\beta}}{2}}) appearing before the double q-derivation
544 : !Take also into account here the -i factor and the complex conjugate
545 165888 : d3etot_tgeom_k(1,idq)=d3etot_tgeom_k(1,idq)-occ_k(iband)*half*doti*two_pi
546 207360 : d3etot_tgeom_k(2,idq)=d3etot_tgeom_k(2,idq)-occ_k(iband)*half*dotr*two_pi
547 :
548 : end do !iband
549 :
550 : !Clean the rf_hamiltonian
551 41472 : call rf_hamkq%free()
552 :
553 : !Deallocations
554 41472 : ABI_FREE(kpg_k)
555 41472 : ABI_FREE(kpg1_k)
556 41472 : ABI_FREE(dkinpw)
557 41472 : ABI_FREE(kinpw1)
558 41472 : ABI_FREE(ffnlk)
559 41472 : ABI_FREE(ffnl1)
560 41472 : ABI_FREE(ph3d)
561 :
562 : end if
563 :
564 : end do !ii
565 :
566 : !-----------------------------------------------------------------------------------------------
567 : ! 2nd q-gradient of atomic displacement 1st order hamiltonian * momentum operator :
568 : ! <u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}}_{\gamma\delta} (k+G)_{\beta} | u_{i,k}^{(0)} >
569 : !-----------------------------------------------------------------------------------------------
570 :
571 : !Get q-gradient of first-order local part of the pseudopotential
572 : ! call dfpt_vlocaldqdq(atindx,2,gs_hamkq%gmet,gsqcut,i1dir,i1pert,mpi_enreg, &
573 : ! & psps%mqgrid_vl,dtset%natom,&
574 : ! & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), &
575 : ! & ph1d,gamma,delta,psps%qgrid_vl,&
576 : ! & dtset%qptn,ucvol,psps%vlspl,vpsp1dq)
577 :
578 : !Set up q-gradient of local potential vlocal1dq with proper dimensioning
579 419966208 : vpsp1dq(:)=vpsp1_i1pertdqdq(:,isppol,idq)
580 : call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,dtset%nfft,dtset%ngfft,&
581 62208 : & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dq,dum_vlocal,vlocal1dq)
582 :
583 : !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup)
584 : call rf_hamkq%init(2,gs_hamkq,i1pert,&
585 62208 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
586 62208 : call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dq,with_nonlocal=.true.)
587 :
588 : !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
589 62208 : if (dtset%ffnl_lw==0) then
590 41472 : dimffnlk=1
591 207360 : ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,psps%ntypat))
592 55947672 : ffnlk(:,1,:,:)=ffnl_k(:,1,:,:)
593 41472 : dimffnl1=10
594 207360 : ABI_MALLOC(ffnl1,(npw_k,dimffnl1,psps%lmnmax,psps%ntypat))
595 559290096 : ffnl1(:,1:dimffnl1,:,:)=ffnl_k(:,1:dimffnl1,:,:)
596 : end if
597 : call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,i1dir,i1pert,gamma, &
598 : & dtset%natom,rmet,rprimd,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgr,useylmgr1,kg_k, &
599 : & ylm_k,kg_k,ylm_k,ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1, &
600 62208 : & reuse_ffnlk=reuse_ffnlk,reuse_ffnl1=reuse_ffnl1,qdir2=delta)
601 :
602 : !LOOP OVER BANDS
603 311040 : do iband=1,nband_k
604 :
605 248832 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
606 :
607 : !Read ket ground-state wavefunctions
608 73988640 : cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg)
609 :
610 : !Compute < g |H^{\tau_{\kappa\alpha}}_{\gamma\delta} | u_{i,k}^{(0)} >
611 : call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq, &
612 248832 : & i1dir,i1pert,mpi_enreg,optlocal,optnl,gamma,rf_hamkq,qdir2=delta)
613 :
614 :
615 24828768 : do ipw=1,npw_k
616 73988640 : gh1dqpkc(:,ipw)=gh1dqc(:,ipw)*kpg_pk(ipw,beta)
617 : end do
618 :
619 : call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0i,gh1dqpkc, &
620 248832 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
621 :
622 : !Take into account the two pi factor from the term
623 : !(\hat{p}_{k\beta + \frac{q_{\beta}}{2}}) appearing before the double q-derivation
624 248832 : d3etot_tgeom_k(1,idq)=d3etot_tgeom_k(1,idq)-occ_k(iband)*doti*two_pi
625 311040 : d3etot_tgeom_k(2,idq)=d3etot_tgeom_k(2,idq)-occ_k(iband)*dotr*two_pi
626 :
627 : end do !iband
628 :
629 : !Clean the rf_hamiltonian
630 62208 : call rf_hamkq%free()
631 :
632 : !Deallocations
633 62208 : ABI_FREE(kpg_k)
634 62208 : ABI_FREE(kpg1_k)
635 62208 : ABI_FREE(dkinpw)
636 62208 : ABI_FREE(kinpw1)
637 62208 : ABI_FREE(ffnlk)
638 62208 : ABI_FREE(ffnl1)
639 103680 : ABI_FREE(ph3d)
640 :
641 : end do !idq
642 :
643 : !scale by the k-point weight
644 228096 : d3etot_tgeom_k(:,:)=d3etot_tgeom_k(:,:)*wtk_k
645 :
646 : !Deallocations
647 41472 : ABI_FREE(dum_cwaveprj)
648 41472 : ABI_FREE(gh1dqc)
649 41472 : ABI_FREE(gh1dqpkc)
650 41472 : ABI_FREE(gvloc1dqc)
651 41472 : ABI_FREE(gvnl1dqc)
652 41472 : ABI_FREE(vpsp1dq)
653 41472 : ABI_FREE(vlocal1dq)
654 41472 : ABI_FREE(dum_vpsp)
655 41472 : ABI_FREE(dum_vlocal)
656 41472 : ABI_FREE(kpg_pk)
657 41472 : ABI_FREE(cwave0i)
658 41472 : ABI_FREE(part_ylmgr_k)
659 :
660 : DBG_EXIT("COLL")
661 :
662 41472 : end subroutine dfptlw_geom
663 : !!***
664 : end module m_dfptlw_nv
665 : !!***
|
