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