Line data Source code
1 : !!****m* ABINIT/m_paw_dfptnl
2 : !! NAME
3 : !! m_paw_dfptnl
4 : !!
5 : !! FUNCTION
6 : !! This module contains several routines used to compute PAW contributions to a 3rd-order energy
7 : !! or 2nd-order PAW occupancies.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2018-2026 ABINIT group (LB)
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_paw_dfptnl
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_errors
28 : use m_xmpi
29 :
30 : use m_pawang, only : pawang_type
31 : use m_pawrad, only : pawrad_type,simp_gen
32 : use m_pawtab, only : pawtab_type
33 : use m_paw_an, only : paw_an_type
34 : use m_pawrhoij, only : pawrhoij_type
35 : use m_pawcprj, only : pawcprj_type
36 : use m_paw_denpot, only : pawdensities
37 : use m_paral_atom, only : get_my_atmtab, free_my_atmtab
38 :
39 : implicit none
40 :
41 : private
42 :
43 : !public procedures.
44 : public :: paw_dfptnl_energy ! Compute the XC PAW on-site contributions to a 3rd-order energy
45 : public :: paw_dfptnl_xc ! Compute a contribution of the 3rd-derivative of XC energy of ONE PAW sphere
46 : public :: paw_dfptnl_accrhoij ! Accumulate the 2nd order PAW quantities rhoij^(2)
47 :
48 : CONTAINS !========================================================================================
49 : !!***
50 :
51 : !----------------------------------------------------------------------
52 :
53 : !!****f* m_paw_dfptnl/paw_dfptnl_energy
54 : !! NAME
55 : !! paw_dfptnl_energy
56 : !!
57 : !! FUNCTION
58 : !! Compute the XC PAW on-site contributions to a 3rd-order energy.
59 : !! It is equal to:
60 : !! E_onsite= \sum_at [ E_at(kxc,rho1,rho2,rho3) - E_at(tkxc,trho1,trho2,trho3) ]
61 : !! where E_at(...) is computed in paw_dfptnl_xc.F90.
62 : !! The atomic densities are computed from pawrhoij1,pawrhoij2 and pawrhoij3.
63 : !! This routine is similar to pawdfptenergy.F90 but is implemented independently
64 : !! in order to not overload the original routine.
65 : !! LDA ONLY - USE THE DENSITY OVER A WHOLE SPHERICAL GRID (r,theta,phi)
66 : !!
67 : !! INPUTS
68 : !! ixc= choice of exchange-correlation scheme
69 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
70 : !! comm_atom=--optional-- MPI communicator over atoms
71 : !! my_natom=number of atoms treated by current processor
72 : !! natom=total number of atoms in cell
73 : !! ntypat=number of types of atoms in unit cell.
74 : !! paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh
75 : !! paw_an1(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh
76 : !! This corresponds to (j1) perturbation
77 : !! pawang <type(pawang_type)>=paw angular mesh and related data
78 : !! pawprtvol=control print volume and debugging output for PAW
79 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
80 : !! pawrhoij_1-2-3(natom) <type(pawrhoij_type)>= paw rhoij 1st-order occupancies
81 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
82 : !! pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
83 : !!
84 : !! OUTPUT
85 : !! d3exc= real and imaginary parts of the contribution to the third derivative of the total energy
86 : !!
87 : !! SIDE EFFECTS
88 : !!
89 : !! SOURCE
90 :
91 246 : subroutine paw_dfptnl_energy(d3exc,ixc,my_natom,natom,ntypat,&
92 246 : & paw_an0,pawang,pawprtvol,pawrad,&
93 246 : & pawrhoij_1,pawrhoij_2,pawrhoij_3,&
94 246 : & pawtab,pawxcdev,&
95 246 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
96 :
97 : !Arguments ---------------------------------------------
98 : !scalars
99 : integer,intent(in) :: ixc,my_natom,natom,ntypat
100 : integer,intent(in) :: pawprtvol,pawxcdev
101 : integer,optional,intent(in) :: comm_atom
102 : type(pawang_type),intent(in) :: pawang
103 : !arrays
104 : integer,optional,target,intent(in) :: mpi_atmtab(:)
105 : real(dp),intent(out) :: d3exc(2)
106 : type(paw_an_type),intent(in) :: paw_an0(my_natom)
107 : type(pawrad_type),intent(in) :: pawrad(ntypat)
108 : type(pawrhoij_type),intent(in) :: pawrhoij_1(my_natom)
109 : type(pawrhoij_type),intent(in) :: pawrhoij_2(my_natom)
110 : type(pawrhoij_type),intent(in) :: pawrhoij_3(my_natom)
111 : type(pawtab_type),intent(in) :: pawtab(ntypat)
112 :
113 : !Local variables ---------------------------------------
114 : !scalars
115 : integer :: cplex_1,cplex_2,cplex_3,iatom,iatom_tot,itypat
116 : integer :: lm_size_all,mesh_size,my_comm_atom,npts,nspden,nzlmopt
117 : integer :: opt_compch,usecore,usetcore,usexcnhat
118 : logical :: my_atmtab_allocated,paral_atom
119 : real(dp) :: compch,d3exc1_iat(2)
120 : character(len=500) :: msg
121 : !arrays
122 246 : integer,pointer :: my_atmtab(:)
123 246 : logical,allocatable :: lmselect_1(:),lmselect_2(:),lmselect_3(:),lmselect_tmp(:)
124 246 : real(dp),allocatable :: nhat1_1(:,:,:),rho1_1(:,:,:),trho1_1(:,:,:)
125 246 : real(dp),allocatable :: nhat1_2(:,:,:),rho1_2(:,:,:),trho1_2(:,:,:)
126 246 : real(dp),allocatable :: nhat1_3(:,:,:),rho1_3(:,:,:),trho1_3(:,:,:)
127 :
128 : ! *************************************************************************
129 :
130 : DBG_ENTER("COLL")
131 :
132 246 : nzlmopt = 0 ! compute all LM-moments of the density and use all LM-moments
133 :
134 246 : if (pawxcdev/=0) then
135 0 : msg="paw_dfptnl_energy is not implemented for pawxcdev/=0"
136 0 : ABI_BUG(msg)
137 : end if
138 246 : if (my_natom>0) then
139 246 : if (pawrhoij_1(1)%qphase/=1.or.pawrhoij_2(1)%qphase/=1.or.pawrhoij_3(1)%qphase/=1) then
140 0 : msg="paw_dfptnl_energy not supposed to be called with q/=0!"
141 0 : ABI_BUG(msg)
142 : end if
143 : end if
144 :
145 : !Set up parallelism over atoms
146 246 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
147 246 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
148 246 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
149 246 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
150 :
151 : !!Various inits
152 246 : opt_compch=0; !optvxc=1;optexc=3
153 246 : usecore=0;usetcore=0 ! This is true for phonons and Efield pert.
154 738 : usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
155 :
156 246 : npts=pawang%angl_size
157 :
158 246 : d3exc = zero
159 :
160 : !================ Loop on atomic sites =======================
161 738 : do iatom=1,my_natom
162 492 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
163 :
164 492 : itypat=pawrhoij_1(iatom)%itypat
165 492 : mesh_size=pawtab(itypat)%mesh_size
166 492 : nspden=pawrhoij_1(iatom)%nspden
167 492 : cplex_1=pawrhoij_1(iatom)%cplex_rhoij
168 492 : cplex_2=pawrhoij_2(iatom)%cplex_rhoij
169 492 : cplex_3=pawrhoij_3(iatom)%cplex_rhoij
170 492 : lm_size_all=paw_an0(iatom)%lm_size
171 :
172 1476 : ABI_MALLOC(lmselect_tmp,(lm_size_all))
173 4920 : lmselect_tmp(:)=.true.
174 :
175 : ! Compute on-site 1st-order densities (pert1)
176 984 : ABI_MALLOC(lmselect_1,(lm_size_all))
177 4920 : lmselect_1(:)=paw_an0(iatom)%lmselect(:)
178 2460 : ABI_MALLOC(rho1_1,(cplex_1*mesh_size,lm_size_all,nspden))
179 1968 : ABI_MALLOC(trho1_1,(cplex_1*mesh_size,lm_size_all,nspden))
180 2460 : ABI_MALLOC(nhat1_1,(cplex_1*mesh_size,lm_size_all,nspden*usexcnhat))
181 : call pawdensities(compch,cplex_1,iatom_tot,lmselect_tmp,lmselect_1,&
182 : & lm_size_all,nhat1_1,nspden,nzlmopt,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,&
183 492 : & pawrad(itypat),pawrhoij_1(iatom),pawtab(itypat),rho1_1,trho1_1)
184 : ! Compute on-site 1st-order densities (pert2)
185 984 : ABI_MALLOC(lmselect_2,(lm_size_all))
186 4920 : lmselect_2(:)=paw_an0(iatom)%lmselect(:)
187 2460 : ABI_MALLOC(rho1_2,(cplex_2*mesh_size,lm_size_all,nspden))
188 1968 : ABI_MALLOC(trho1_2,(cplex_2*mesh_size,lm_size_all,nspden))
189 2460 : ABI_MALLOC(nhat1_2,(cplex_2*mesh_size,lm_size_all,nspden*usexcnhat))
190 : call pawdensities(compch,cplex_2,iatom_tot,lmselect_tmp,lmselect_2,&
191 : & lm_size_all,nhat1_2,nspden,nzlmopt,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,&
192 492 : & pawrad(itypat),pawrhoij_2(iatom),pawtab(itypat),rho1_2,trho1_2)
193 : ! Compute on-site 1st-order densities (pert3)
194 984 : ABI_MALLOC(lmselect_3,(lm_size_all))
195 4920 : lmselect_3(:)=paw_an0(iatom)%lmselect(:)
196 2460 : ABI_MALLOC(rho1_3,(cplex_3*mesh_size,lm_size_all,nspden))
197 1968 : ABI_MALLOC(trho1_3,(cplex_3*mesh_size,lm_size_all,nspden))
198 2460 : ABI_MALLOC(nhat1_3,(cplex_3*mesh_size,lm_size_all,nspden*usexcnhat))
199 : call pawdensities(compch,cplex_3,iatom_tot,lmselect_tmp,lmselect_3,&
200 : & lm_size_all,nhat1_3,nspden,nzlmopt,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,&
201 492 : & pawrad(itypat),pawrhoij_3(iatom),pawtab(itypat),rho1_3,trho1_3)
202 492 : ABI_FREE(lmselect_tmp)
203 :
204 : call paw_dfptnl_xc(cplex_1,cplex_2,cplex_3,d3exc1_iat,ixc,paw_an0(iatom)%k3xc1,lm_size_all,&
205 : & lmselect_1,lmselect_2,lmselect_3,nhat1_1,nhat1_2,nhat1_3,&
206 : & paw_an0(iatom)%nk3xc1,mesh_size,nspden,pawang,pawrad(itypat),&
207 492 : & rho1_1,rho1_2,rho1_3,0)
208 1476 : d3exc = d3exc + d3exc1_iat
209 :
210 : call paw_dfptnl_xc(cplex_1,cplex_2,cplex_3,d3exc1_iat,ixc,paw_an0(iatom)%k3xct1,lm_size_all,&
211 : & lmselect_1,lmselect_2,lmselect_3,nhat1_1,nhat1_2,nhat1_3,&
212 : & paw_an0(iatom)%nk3xc1,mesh_size,nspden,pawang,pawrad(itypat),&
213 492 : & trho1_1,trho1_2,trho1_3,usexcnhat)
214 1476 : d3exc = d3exc - d3exc1_iat
215 :
216 492 : ABI_FREE(lmselect_1)
217 492 : ABI_FREE(lmselect_2)
218 492 : ABI_FREE(lmselect_3)
219 492 : ABI_FREE(nhat1_1)
220 492 : ABI_FREE(nhat1_2)
221 492 : ABI_FREE(nhat1_3)
222 492 : ABI_FREE(rho1_1)
223 492 : ABI_FREE(rho1_2)
224 492 : ABI_FREE(rho1_3)
225 492 : ABI_FREE(trho1_1)
226 492 : ABI_FREE(trho1_2)
227 738 : ABI_FREE(trho1_3)
228 :
229 : ! ================ End loop oon atomic sites =======================
230 : end do
231 :
232 : !!Reduction in case of parallelism
233 : ! if (paral_atom) then
234 : ! call xmpi_sum(delta_energy,my_comm_atom,ierr)
235 : ! end if
236 :
237 : !!Destroy atom table used for parallelism
238 : ! call free_my_atmtab(my_atmtab,my_atmtab_allocated)
239 :
240 : ! call timab(567,2,tsec)
241 :
242 : DBG_EXIT("COLL")
243 :
244 246 : end subroutine paw_dfptnl_energy
245 : !!***
246 :
247 : !----------------------------------------------------------------------
248 :
249 : !!****f* m_paw_dfptnl/paw_dfptnl_xc
250 : !! NAME
251 : !! paw_dfptnl_xc
252 : !!
253 : !! FUNCTION
254 : !! Compute a contribution of the third derivative of XC energy of ONE PAW sphere Om_a.
255 : !! It is equal to:
256 : !! E_at(kxc,rho1,rho2,rho3) = Int_{Om_a} dr kxc(r) * rho1(r) * rho2(r) * rho3(r)
257 : !! where kxc,rho1,rho2 and rho3 are inputs.
258 : !! This routine is similar to m_pawxc.F90:pawxc_dfpt(...) but is implemented independently
259 : !! in order to not overload the original routine.
260 : !! LDA ONLY - USE THE DENSITY OVER A WHOLE SPHERICAL GRID (r,theta,phi)
261 : !!
262 : !! INPUTS
263 : !! cplex_1-2-3= if 1, 1st-order densities are REAL, if 2, COMPLEX
264 : !! d3exc1_iat=third-order derivative to compute
265 : !! ixc= choice of exchange-correlation scheme
266 : !! kxc(nrad,pawang%angl_size,nkxc)=GS xc kernel
267 : !! lm_size=size of density array rhor (see below)
268 : !! lmselect1-2-3(lm_size)=select the non-zero LM-moments of input density rhor1-2-3
269 : !! nhat1-2-3(cplex_den*nrad,lm_size,nspden)=first-order change of compensation density
270 : !! (total in 1st half and spin-up in 2nd half if nspden=2)
271 : !! nkxc=second dimension of the kxc array
272 : !! nrad=size of radial mesh for densities/potentials (might be different from pawrad%mesh_size)
273 : !! nspden=number of spin-density components
274 : !! option=0 compute both 2nd-order XC energy and 1st-order potential
275 : !! 1 compute only 1st-order XC potential
276 : !! 2 compute only 2nd-order XC energy, XC potential is temporary computed here
277 : !! 3 compute only 2nd-order XC energy, XC potential is input in vxc1(:)
278 : !! pawang <type(pawang_type)>=paw angular mesh and related data
279 : !! pawrad <type(pawrad_type)>=paw radial mesh and related data
280 : !! rhor1-2-3(cplex_den*nrad,lm_size,nspden)=first-order change of density
281 : !! usexcnhat= 0 if compensation density does not have to be used
282 : !! 1 if compensation density has to be used in d2Exc only
283 : !!
284 : !! OUTPUT
285 : !! d3exc1_iat = E_at(kxc,rho1,rho2,rho3) (see FUNCTION above)
286 : !!
287 : !! SIDE EFFECTS
288 : !!
289 : !! SOURCE
290 :
291 984 : subroutine paw_dfptnl_xc(cplex_1,cplex_2,cplex_3,d3exc1_iat,ixc,kxc,lm_size,lmselect1,lmselect2,lmselect3,&
292 984 : & nhat1,nhat2,nhat3,nkxc,nrad,nspden,pawang,pawrad,rhor1,rhor2,rhor3,usexcnhat)
293 :
294 : !Arguments ------------------------------------
295 : !scalars
296 : integer,intent(in) :: cplex_1,cplex_2,cplex_3,ixc,lm_size,nkxc,nrad,nspden,usexcnhat
297 : type(pawang_type),intent(in) :: pawang
298 : type(pawrad_type),intent(in) :: pawrad
299 : !arrays
300 : logical,intent(in) :: lmselect1(lm_size),lmselect2(lm_size),lmselect3(lm_size)
301 : real(dp),intent(out) :: d3exc1_iat(2)
302 : real(dp),intent(in) :: kxc(nrad,pawang%angl_size,nkxc)
303 : real(dp),intent(in) :: nhat1(cplex_1*nrad,lm_size,nspden*((usexcnhat+1)/2))
304 : real(dp),intent(in) :: nhat2(cplex_2*nrad,lm_size,nspden*((usexcnhat+1)/2))
305 : real(dp),intent(in) :: nhat3(cplex_3*nrad,lm_size,nspden*((usexcnhat+1)/2))
306 : real(dp),intent(in) :: rhor1(cplex_1*nrad,lm_size,nspden)
307 : real(dp),intent(in) :: rhor2(cplex_2*nrad,lm_size,nspden)
308 : real(dp),intent(in) :: rhor3(cplex_3*nrad,lm_size,nspden)
309 :
310 : !Local variables-------------------------------
311 : !scalars
312 : integer :: ii,ilm,ipts,ispden,lm_size_eff,npts
313 : real(dp) :: d3exc1_int,rho1u,rho1d,rho2u,rho2d,rho3u,rho3d
314 : character(len=500) :: msg
315 : !arrays
316 : ! real(dp) :: tsec(2)
317 984 : real(dp),allocatable :: ff(:),rho1arr(:,:),rho2arr(:,:),rho3arr(:,:)
318 :
319 : ! *************************************************************************
320 :
321 : !----------------------------------------------------------------------
322 : !----- Initializations
323 : !----------------------------------------------------------------------
324 :
325 984 : npts=pawang%angl_size
326 984 : lm_size_eff=min(lm_size,pawang%ylm_size)
327 :
328 984 : d3exc1_iat(:) = zero
329 :
330 : !Special case: no XC applied
331 984 : if (ixc==0) then
332 0 : msg='Note that no xc is applied (ixc=0). Returning'
333 0 : ABI_WARNING(msg)
334 : return
335 : end if
336 :
337 3936 : ABI_MALLOC(rho1arr,(cplex_1*nrad,nspden))
338 3936 : ABI_MALLOC(rho2arr,(cplex_2*nrad,nspden))
339 3936 : ABI_MALLOC(rho3arr,(cplex_3*nrad,nspden))
340 :
341 : !Restriction : all cplex must be 1
342 984 : if (cplex_1/=1.or.cplex_2/=1.or.cplex_3/=1) then
343 0 : msg='All cplex must be one (for the moment...)'
344 0 : ABI_BUG(msg)
345 : end if
346 : !Restriction : nspden must be 1
347 : ! if (nkxc>1) then
348 : ! msg='nkxc must be one (<=> nspden=1) (for the moment...)'
349 : ! ABI_BUG(msg)
350 : ! end if
351 :
352 2952 : ABI_MALLOC(ff,(nrad))
353 :
354 : !!----------------------------------------------------------------------
355 : !!----- Loop on the angular part and inits
356 : !!----------------------------------------------------------------------
357 :
358 : !Do loop on the angular part (theta,phi)
359 154488 : do ipts=1,npts
360 :
361 : ! Copy the input 1st-order density for this (theta,phi) - PERT1
362 366558192 : rho1arr(:,:)=zero
363 153504 : if (usexcnhat==0) then
364 451152 : do ispden=1,nspden
365 3129984 : do ilm=1,lm_size_eff
366 2678832 : if (lmselect1(ilm)) rho1arr(:,ispden)=rho1arr(:,ispden) &
367 3297939840 : & +rhor1(:,ilm,ispden)*pawang%ylmr(ilm,ipts)
368 : end do
369 : end do
370 : else
371 0 : do ispden=1,nspden
372 0 : do ilm=1,lm_size_eff
373 0 : if (lmselect1(ilm)) rho1arr(:,ispden)=rho1arr(:,ispden) &
374 0 : & +(rhor1(:,ilm,ispden)+nhat1(:,ilm,ispden))*pawang%ylmr(ilm,ipts)
375 : end do
376 : end do
377 : end if
378 :
379 : ! Copy the input 1st-order density for this (theta,phi) - PERT2
380 366558192 : rho2arr(:,:)=zero
381 153504 : if (usexcnhat==0) then
382 451152 : do ispden=1,nspden
383 3129984 : do ilm=1,lm_size_eff
384 2678832 : if (lmselect2(ilm)) rho2arr(:,ispden)=rho2arr(:,ispden) &
385 3297939840 : & +rhor2(:,ilm,ispden)*pawang%ylmr(ilm,ipts)
386 : end do
387 : end do
388 : else
389 0 : do ispden=1,nspden
390 0 : do ilm=1,lm_size_eff
391 0 : if (lmselect2(ilm)) rho2arr(:,ispden)=rho2arr(:,ispden) &
392 0 : & +(rhor2(:,ilm,ispden)+nhat2(:,ilm,ispden))*pawang%ylmr(ilm,ipts)
393 : end do
394 : end do
395 : end if
396 :
397 : ! Copy the input 1st-order density for this (theta,phi) - PERT3
398 366558192 : rho3arr(:,:)=zero
399 153504 : if (usexcnhat==0) then
400 451152 : do ispden=1,nspden
401 3129984 : do ilm=1,lm_size_eff
402 2678832 : if (lmselect3(ilm)) rho3arr(:,ispden)=rho3arr(:,ispden) &
403 3297939840 : & +rhor3(:,ilm,ispden)*pawang%ylmr(ilm,ipts)
404 : end do
405 : end do
406 : else
407 0 : do ispden=1,nspden
408 0 : do ilm=1,lm_size_eff
409 0 : if (lmselect3(ilm)) rho3arr(:,ispden)=rho3arr(:,ispden) &
410 0 : & +(rhor3(:,ilm,ispden)+nhat3(:,ilm,ispden))*pawang%ylmr(ilm,ipts)
411 : end do
412 : end do
413 : end if
414 :
415 : ! ----------------------------------------------------------------------
416 : ! ----- Accumulate and store 3nd-order change of XC energy
417 : ! ----------------------------------------------------------------------
418 :
419 153504 : if (cplex_1==1.and.cplex_2==1.and.cplex_3==1) then ! all cplex are 1 :
420 153504 : if (nspden==1) then
421 11522160 : ff(:)=kxc(:,ipts,1)*rho1arr(:,1)*rho2arr(:,1)*rho3arr(:,1)
422 144144 : else if (nspden==2) then
423 177441264 : do ii=1,nrad
424 177297120 : rho1u=rho1arr(ii,2)
425 177297120 : rho1d=rho1arr(ii,1)-rho1arr(ii,2)
426 177297120 : rho2u=rho2arr(ii,2)
427 177297120 : rho2d=rho2arr(ii,1)-rho2arr(ii,2)
428 177297120 : rho3u=rho3arr(ii,2)
429 177297120 : rho3d=rho3arr(ii,1)-rho3arr(ii,2)
430 : ff(ii)=&
431 : ! uuu uud
432 : & kxc(ii,ipts,1)*rho1u*rho2u*rho3u + kxc(ii,ipts,2)*rho1u*rho2u*rho3d + &
433 : ! udu udd
434 : & kxc(ii,ipts,2)*rho1u*rho2d*rho3u + kxc(ii,ipts,3)*rho1u*rho2d*rho3d + &
435 : ! duu dud
436 : & kxc(ii,ipts,2)*rho1d*rho2u*rho3u + kxc(ii,ipts,3)*rho1d*rho2u*rho3d + &
437 : ! ddu ddd
438 177441264 : & kxc(ii,ipts,3)*rho1d*rho2d*rho3u + kxc(ii,ipts,4)*rho1d*rho2d*rho3d
439 : end do
440 : end if
441 : end if
442 :
443 188963424 : ff(1:nrad)=ff(1:nrad)*pawrad%rad(1:nrad)**2
444 153504 : call simp_gen(d3exc1_int,ff,pawrad)
445 154488 : d3exc1_iat(1)=d3exc1_iat(1)+d3exc1_int*pawang%angwgth(ipts)
446 :
447 : ! ----- End of the loop on npts (angular part)
448 : end do
449 :
450 2952 : d3exc1_iat = d3exc1_iat*four_pi
451 :
452 984 : ABI_FREE(ff)
453 984 : ABI_FREE(rho1arr)
454 984 : ABI_FREE(rho2arr)
455 984 : ABI_FREE(rho3arr)
456 :
457 : end subroutine paw_dfptnl_xc
458 : !!***
459 :
460 : !----------------------------------------------------------------------
461 :
462 : !!****f* ABINIT/paw_dfptnl_accrhoij
463 : !!
464 : !! NAME
465 : !! paw_dfptnl_accrhoij
466 : !!
467 : !! FUNCTION
468 : !! Accumulate the 2nd order PAW quantities rhoij^(2) (augmentation occupancies)
469 : !! This routine is similar to pawaccrhoij.F90 but is implemented independently
470 : !! in order to not overload the original routine.
471 : !!
472 : !! INPUTS
473 : !! atindx(natom)=index table for atoms (sorted-->random), inverse of atindx.
474 : !! cplex: if 1, WFs (or 1st-order WFs) are REAL, if 2, COMPLEX
475 : !! cwaveprj0_pert1(natom,nspinor) = wave function at given n,k projected with non-local projectors:
476 : !! cwaveprj0%cp =<p_i|
477 : !! cwaveprj0%dcp(1)=<p_i^(pert1)|Cnk>
478 : !! cwaveprj0_pert2(natom,nspinor) = wave function at given n,k projected with non-local projectors:
479 : !! cwaveprj0%cp =<p_i|Cnk>
480 : !! cwaveprj0%dcp(1)=<p_i^(pert1)|Cnk>
481 : !! cwaveprj1_pert12(natom,nspinor)= 1st order wave function at given n,k projected with non-local projectors:
482 : !! cwaveprj1%cp =<p_i|Cnk^(pert1)>
483 : !! cwaveprj1%dcp(1)=<p_i^(pert2)|Cnk^(pert1)>
484 : !! cwaveprj1_pert21(natom,nspinor)= 1st order wave function at given n,k projected with non-local projectors:
485 : !! cwaveprj1%cp =<p_i|Cnk^(pert2)>
486 : !! cwaveprj1%dcp(1)=<p_i^(pert1)|Cnk^(pert2)>
487 : !! ipert1=index of the first perturbation
488 : !! ipert2=index of the second perturbation
489 : !! isppol=index of current spin component
490 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
491 : !! comm_atom=--optional-- MPI communicator over atoms
492 : !! my_natom=number of atoms treated by current processor
493 : !! natom=number of atoms in cell
494 : !! nspinor=number of spinorial components (on current proc)
495 : !! occ_k=occupation number for current band n,k
496 : !! wtk_k=weight assigned to current k-point
497 : !!
498 : !! SIDE EFFECTS
499 : !! pawrhoij(natom) <type(pawrhoij_type)>= 2-nd order paw rhoij occupancies and related data
500 : !! On output, has been updated with the contribution of current n,k
501 : !! pawrhoij(:)%rhoij_(lmn2_size,nspden) (non symetrized)
502 : !!
503 : !! SOURCE
504 :
505 16128 : subroutine paw_dfptnl_accrhoij(atindx,cplex,cwaveprj0_pert1,cwaveprj0_pert2,&
506 8064 : & cwaveprj1_pert12,cwaveprj1_pert21,ipert1,ipert2,isppol,my_natom,natom,&
507 8064 : & nspinor,occ_k,pawrhoij,wtk_k,&
508 8064 : & comm_atom,mpi_atmtab ) ! optional (parallelism)
509 :
510 : !Arguments ---------------------------------------------
511 : !scalars
512 : integer,intent(in) :: cplex,ipert1,ipert2,isppol,my_natom,natom,nspinor
513 : integer,optional,intent(in) :: comm_atom
514 : real(dp),intent(in) :: occ_k,wtk_k
515 : !arrays
516 : integer,intent(in) :: atindx(natom)
517 : integer,optional,target,intent(in) :: mpi_atmtab(:)
518 : type(pawcprj_type),intent(in) :: cwaveprj0_pert1(natom,nspinor),cwaveprj0_pert2(natom,nspinor)
519 : type(pawcprj_type),intent(in) :: cwaveprj1_pert12(natom,nspinor),cwaveprj1_pert21(natom,nspinor)
520 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
521 :
522 : !Local variables ---------------------------------------
523 : !scalars
524 : integer :: cplex_rhoij,iatm,iatom,iatom1,ilmn,iplex,j0lmn,jlmn,klmn,klmn_im,klmn_re
525 : integer :: my_comm_atom,ncpgr
526 : logical :: compute_impart,compute_impart_cplex
527 : logical :: my_atmtab_allocated,paral_atom
528 : real(dp) :: ro11_im,ro11_re,weight
529 : character(len=500) :: message
530 : !arrays
531 8064 : integer,pointer :: my_atmtab(:)
532 16128 : real(dp) :: cpi0(2,nspinor),d1cpi0(2,nspinor),d2cpi0(2,nspinor)
533 16128 : real(dp) :: cpj0(2,nspinor),d1cpj0(2,nspinor),d2cpj0(2,nspinor)
534 16128 : real(dp) :: cpi1(2,nspinor),d2cpi1(2,nspinor)
535 16128 : real(dp) :: cpj1(2,nspinor),d2cpj1(2,nspinor)
536 16128 : real(dp) :: cpi2(2,nspinor),d1cpi2(2,nspinor)
537 8064 : real(dp) :: cpj2(2,nspinor),d1cpj2(2,nspinor)
538 :
539 : ! ***********************************************************************
540 :
541 : DBG_ENTER("COLL")
542 :
543 8064 : if (my_natom==0) return
544 :
545 8064 : ncpgr=1
546 8064 : if (ipert1<0.or.ipert1>natom+2.or.ipert2<0.or.ipert2>natom+2) then
547 0 : message = 'paw_dfptnl_accrhoij: Necessary conditions on ipert1 or ipert2: 0<=ipert<=natom+2'
548 0 : ABI_BUG(message)
549 : end if
550 8064 : if (pawrhoij(1)%qphase/=1) then
551 0 : message="paw_dfptnl_accrhoij not supposed to be called with q/=0!"
552 0 : ABI_BUG(message)
553 : end if
554 :
555 : !Set up parallelism over atoms
556 8064 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
557 8064 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
558 8064 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
559 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
560 8064 : & my_natom_ref=my_natom)
561 :
562 8064 : weight=wtk_k*occ_k
563 8064 : if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1) weight=half*weight
564 :
565 : !! ==================================================================
566 : !! === Accumulate (n,k) contribution to partial 2nd-order rhoij ===
567 : !! ==================================================================
568 :
569 8064 : compute_impart=(pawrhoij(1)%cplex_rhoij==2)
570 8064 : compute_impart_cplex=((pawrhoij(1)%cplex_rhoij==2).and.(cplex==2))
571 :
572 : ! NOT USED FOR PAWRHO21! => only for PAWRHO2 (full second derivative)
573 : !!Accumulate : < Psi^(pert1) | p_i^(0) > < p_j^(0) | Psi^(pert2) >
574 : !! + < Psi^(pert2) | p_i^(0) > < p_j^(0) | Psi^(pert1) >
575 : ! if (nspinor==1) then
576 : ! do iatom=1,my_natom
577 : ! iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
578 : ! iatm=atindx(iatom1)
579 : ! cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
580 : ! do jlmn=1,pawrhoij(iatom)%lmn_size
581 : ! j0lmn=jlmn*(jlmn-1)/2
582 : ! cpj1(1:2,1)=cwaveprj1_pert12(iatm,1)%cp(1:2,jlmn) ! < p_j^(0) | Psi^(pert1) >
583 : ! cpj2(1:2,1)=cwaveprj1_pert21(iatm,1)%cp(1:2,jlmn) ! < p_j^(0) | Psi^(pert2) >
584 : ! do ilmn=1,jlmn
585 : ! klmn=j0lmn+ilmn
586 : ! klmn_re=cplex_rhoij*(klmn-1)+1
587 : ! cpi1(1:2,1)=cwaveprj1_pert12(iatm,1)%cp(1:2,ilmn) ! < p_i^(0) | Psi^(pert1) >
588 : ! cpi2(1:2,1)=cwaveprj1_pert21(iatm,1)%cp(1:2,ilmn) ! < p_i^(0) | Psi^(pert2) >
589 : ! ro11_re=zero
590 : ! do iplex=1,cplex
591 : ! ro11_re=ro11_re+cpi1(iplex,1)*cpj2(iplex,1)+cpj1(iplex,1)*cpi2(iplex,1)
592 : ! end do
593 : ! pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
594 : ! if (compute_impart_cplex) then
595 : ! klmn_im=klmn_re+1
596 : ! ro11_im= cpi1(1,1)*cpj2(2,1)-cpi1(2,1)*cpj2(1,1)
597 : ! ro11_im=ro11_im+cpj1(1,1)*cpi2(2,1)-cpj1(2,1)*cpi2(1,1)
598 : ! pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
599 : ! end if
600 : ! end do
601 : ! end do
602 : ! end do
603 : ! else ! nspinor=2
604 : ! ABI_BUG("paw_dfptnl_accrhoij is not implemented for nspinor=2")
605 : ! end if
606 :
607 : !Accumulate : < Psi^(pert1) | p_i^(pert2) > < p_j^(0) | Psi^(0) >
608 : ! + < Psi^(pert1) | p_i^(0) > < p_j^(pert2) | Psi^(0) >
609 : ! + < Psi^(0) | p_i^(pert2) > < p_j^(0) | Psi^(pert1) >
610 : ! + < Psi^(0) | p_i^(0) > < p_j^(pert2) | Psi^(pert1) >
611 8064 : if (ipert2>0.and.ipert2<=natom) then
612 0 : if (nspinor==1) then
613 0 : do iatom=1,my_natom
614 0 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
615 0 : iatm=atindx(iatom1)
616 0 : if (iatom/=ipert2) cycle ! To move atom "ipert2" does not change projectors of other atoms
617 0 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
618 0 : do jlmn=1,pawrhoij(iatom)%lmn_size
619 0 : j0lmn=jlmn*(jlmn-1)/2
620 0 : cpj0(1:2,1) =cwaveprj0_pert2 (iatm,1)% cp(1:2, jlmn) ! < p_j^(0) | Psi^(0) >
621 0 : d2cpj0(1:2,1)=cwaveprj0_pert2 (iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert2) | Psi^(0) >
622 0 : cpj1(1:2,1) =cwaveprj1_pert12(iatm,1)% cp(1:2, jlmn) ! < p_j^(0) | Psi^(pert1) >
623 0 : d2cpj1(1:2,1)=cwaveprj1_pert12(iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert2) | Psi^(pert1) >
624 0 : do ilmn=1,jlmn
625 0 : klmn=j0lmn+ilmn
626 0 : klmn_re=cplex_rhoij*(klmn-1)+1
627 0 : cpi0(1:2,1) =cwaveprj0_pert2 (iatm,1)% cp(1:2, ilmn) ! < p_i^(0) | Psi^(0) >
628 0 : d2cpi0(1:2,1)=cwaveprj0_pert2 (iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert2) | Psi^(0) >
629 0 : cpi1(1:2,1) =cwaveprj1_pert12(iatm,1)% cp(1:2, ilmn) ! < p_i^(0) | Psi^(pert1) >
630 0 : d2cpi1(1:2,1)=cwaveprj1_pert12(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert2) | Psi^(pert1) >
631 0 : ro11_re=zero
632 0 : do iplex=1,cplex
633 0 : ro11_re=ro11_re+d2cpi1(iplex,1)* cpj0(iplex,1)
634 0 : ro11_re=ro11_re+ cpi1(iplex,1)*d2cpj0(iplex,1)
635 0 : ro11_re=ro11_re+d2cpi0(iplex,1)* cpj1(iplex,1)
636 0 : ro11_re=ro11_re+ cpi0(iplex,1)*d2cpj1(iplex,1)
637 : end do
638 0 : pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
639 0 : if (compute_impart_cplex) then
640 0 : klmn_im=klmn_re+1
641 0 : ro11_im= d2cpi1(1,1)* cpj0(2,1)-d2cpi1(2,1)* cpj0(1,1)
642 0 : ro11_im=ro11_im+ cpi1(1,1)*d2cpj0(2,1)- cpi1(2,1)*d2cpj0(1,1)
643 0 : ro11_im=ro11_im+d2cpi0(1,1)* cpj1(2,1)-d2cpi0(2,1)* cpj1(1,1)
644 0 : ro11_im=ro11_im+ cpi0(1,1)*d2cpj1(2,1)- cpi0(2,1)*d2cpj1(1,1)
645 0 : pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
646 : end if
647 : end do
648 : end do
649 : end do
650 : else ! nspinor=2
651 0 : ABI_BUG("paw_dfptnl_accrhoij is not implemented for nspinor=2")
652 : end if
653 : end if
654 :
655 : !Accumulate : < Psi^(pert2) | p_i^(pert1) > < p_j^(0) | Psi^(0) >
656 : ! + < Psi^(pert2) | p_i^(0) > < p_j^(pert1) | Psi^(0) >
657 : ! + < Psi^(0) | p_i^(pert1) > < p_j^(0) | Psi^(pert2) >
658 : ! + < Psi^(0) | p_i^(0) > < p_j^(pert1) | Psi^(pert2) >
659 8064 : if (ipert1>0.and.ipert1<=natom) then
660 8064 : if (nspinor==1) then
661 24192 : do iatom=1,my_natom
662 16128 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
663 16128 : iatm=atindx(iatom1)
664 16128 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
665 16128 : if (iatom/=ipert1) cycle ! To move atom "ipert1" does not change projectors of other atoms
666 80640 : do jlmn=1,pawrhoij(iatom)%lmn_size
667 64512 : j0lmn=jlmn*(jlmn-1)/2
668 193536 : cpj0(1:2,1) =cwaveprj0_pert1 (iatm,1)% cp(1:2, jlmn) ! < p_j^(0) | Psi^(0) >
669 193536 : d1cpj0(1:2,1)=cwaveprj0_pert1 (iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert1) | Psi^(0) >
670 193536 : cpj2(1:2,1) =cwaveprj1_pert21(iatm,1)% cp(1:2, jlmn) ! < p_j^(0) | Psi^(pert2) >
671 193536 : d1cpj2(1:2,1)=cwaveprj1_pert21(iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert1) | Psi^(pert2) >
672 370944 : do ilmn=1,jlmn
673 290304 : klmn=j0lmn+ilmn
674 290304 : klmn_re=cplex_rhoij*(klmn-1)+1
675 870912 : cpi0(1:2,1) =cwaveprj0_pert1 (iatm,1)% cp(1:2, ilmn) ! < p_i^(0) | Psi^(0) >
676 870912 : d1cpi0(1:2,1)=cwaveprj0_pert1 (iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert1) | Psi^(0) >
677 870912 : cpi2(1:2,1) =cwaveprj1_pert21(iatm,1)% cp(1:2, ilmn) ! < p_i^(0) | Psi^(pert2) >
678 870912 : d1cpi2(1:2,1)=cwaveprj1_pert21(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert1) | Psi^(pert2) >
679 290304 : ro11_re=zero
680 870912 : do iplex=1,cplex
681 580608 : ro11_re=ro11_re+d1cpi2(iplex,1)* cpj0(iplex,1)
682 580608 : ro11_re=ro11_re+ cpi2(iplex,1)*d1cpj0(iplex,1)
683 580608 : ro11_re=ro11_re+d1cpi0(iplex,1)* cpj2(iplex,1)
684 870912 : ro11_re=ro11_re+ cpi0(iplex,1)*d1cpj2(iplex,1)
685 : end do
686 290304 : pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
687 354816 : if (compute_impart_cplex) then
688 0 : klmn_im=klmn_re+1
689 0 : ro11_im= d1cpi2(1,1)* cpj0(2,1)-d1cpi2(2,1)* cpj0(1,1)
690 0 : ro11_im=ro11_im+ cpi2(1,1)*d1cpj0(2,1)- cpi2(2,1)*d1cpj0(1,1)
691 0 : ro11_im=ro11_im+d1cpi0(1,1)* cpj2(2,1)-d1cpi0(2,1)* cpj2(1,1)
692 0 : ro11_im=ro11_im+ cpi0(1,1)*d1cpj2(2,1)- cpi0(2,1)*d1cpj2(1,1)
693 0 : pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
694 : end if
695 : end do
696 : end do
697 : end do
698 : else ! nspinor=2
699 0 : message="paw_dfptnl_accrhoij is not implemented for nspinor=2!"
700 0 : ABI_BUG(message)
701 : end if
702 : end if
703 : ! End
704 :
705 : !Accumulate : < Psi^(0) | p_i^(pert1) > < p_j^(pert2) | Psi^(0) >
706 : ! + < Psi^(0) | p_i^(pert2) > < p_j^(pert1) | Psi^(0) >
707 8064 : if (ipert1>0.and.ipert1<=natom.and.ipert2>0.and.ipert2<=natom) then
708 0 : if (nspinor==1) then
709 0 : do iatom=1,my_natom
710 0 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
711 0 : iatm=atindx(iatom1)
712 0 : if (iatom/=ipert1.or.iatom/=ipert2) cycle ! To move atom "ipert" does not change projectors of other atoms
713 0 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
714 0 : do jlmn=1,pawrhoij(iatom)%lmn_size
715 0 : j0lmn=jlmn*(jlmn-1)/2
716 0 : d1cpj0(1:2,1)=cwaveprj0_pert1(iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert1) | Psi^(0) >
717 0 : d2cpj0(1:2,1)=cwaveprj0_pert2(iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert2) | Psi^(0) >
718 0 : do ilmn=1,jlmn
719 0 : klmn=j0lmn+ilmn
720 0 : klmn_re=cplex_rhoij*(klmn-1)+1
721 0 : d1cpi0(1:2,1)=cwaveprj0_pert1(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert1) | Psi^(0) >
722 0 : d2cpi0(1:2,1)=cwaveprj0_pert2(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert2) | Psi^(0) >
723 0 : ro11_re=zero
724 0 : do iplex=1,cplex
725 0 : ro11_re=ro11_re+d1cpi0(iplex,1)*d2cpj0(iplex,1)+d2cpi0(iplex,1)*d1cpj0(iplex,1)
726 : end do
727 0 : pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
728 0 : if (compute_impart_cplex) then
729 0 : klmn_im=klmn_re+1
730 0 : ro11_im= d1cpi0(1,1)*d2cpj0(2,1)-d1cpi0(2,1)*d2cpj0(1,1)
731 0 : ro11_im=ro11_im+d2cpi0(1,1)*d1cpj0(2,1)-d2cpi0(2,1)*d1cpj0(1,1)
732 0 : pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
733 : end if
734 : end do
735 : end do
736 : end do
737 : else ! nspinor=2
738 0 : message="paw_dfptnl_accrhoij is not implemented for nspinor=2!"
739 0 : ABI_BUG(message)
740 : end if
741 : end if
742 :
743 : !Destroy atom table used for parallelism
744 8064 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
745 :
746 : DBG_EXIT("COLL")
747 :
748 8064 : end subroutine paw_dfptnl_accrhoij
749 : !!***
750 :
751 : !----------------------------------------------------------------------
752 :
753 : END MODULE m_paw_dfptnl
754 : !!***
|
