Branch data 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 : : !!***
|
