Line data Source code
1 : !!****m* ABINIT/m_dfptnl_pert
2 : !! NAME
3 : !! m_dfptnl_pert
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group ()
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_dfptnl_pert
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_wffile
27 : use m_wfk
28 : use m_xmpi
29 : use m_hamiltonian
30 : use m_errors
31 : use m_rf2
32 : use m_kg
33 : use m_dtset
34 : use m_dtfil
35 :
36 : use defs_datatypes, only : pseudopotential_type
37 : use defs_abitypes, only : MPI_type
38 : use m_cgtools, only : dotprod_g,sqnorm_g,dotprod_vn
39 : use m_pawang, only : pawang_type
40 : use m_pawfgrtab, only : pawfgrtab_type
41 : use m_pawrad, only : pawrad_type
42 : use m_pawtab, only : pawtab_type
43 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free
44 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify,paw_ij_reset_flags
45 : use m_pawdij, only : pawdijfr
46 : use m_pawfgr, only : pawfgr_type
47 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc , pawrhoij_nullify, pawrhoij_free,&
48 : & pawrhoij_init_unpacked, pawrhoij_mpisum_unpacked, pawrhoij_inquire_dim
49 : use m_paw_an, only : paw_an_type
50 : use m_paw_mkrho, only : pawmkrho
51 : use m_paw_nhat, only : pawnhatfr
52 : use m_paw_dfpt, only : pawdfptenergy
53 : use m_paw_dfptnl, only : paw_dfptnl_accrhoij,paw_dfptnl_energy
54 : use m_initylmg, only : initylmg
55 : use m_mkffnl, only : mkffnl
56 : use m_getghc, only : getgsc
57 : use m_getgh1c, only : rf_transgrid_and_pack
58 : use m_mpinfo, only : proc_distrb_cycle
59 : use m_nonlop, only : nonlop
60 : use m_fourier_interpol, only : transgrid
61 : use m_cgprj, only : getcprj
62 :
63 : implicit none
64 :
65 : private
66 : !!***
67 :
68 : public :: dfptnl_pert
69 : public :: dfptnl_exc3
70 : !!***
71 :
72 : contains
73 : !!***
74 :
75 : !!****f* ABINIT/dfptnl_pert
76 : !! NAME
77 : !! dfptnl_pert
78 : !!
79 : !! FUNCTION
80 : !! Compute the linear response part to the 3dte. The main inputs are :
81 : !! - GS WFs and Hamiltonian (cg,gs_hamkq)
82 : !! - 1st-order WFs for three perturbations i1pert/i1dir,i2pert/i2dir,i3pert/i3dir (cg1,cg2,cg3)
83 : !! - 1st-order potentials for i2pert (vhartr1_i2pert,vtrial1_i2pert,vxc1_i2pert)
84 : !! - 1st-order WFs DDK,DDE and 2nd-order WF DKDE (ddk_f)
85 : !!
86 : !! COPYRIGHT
87 : !! Copyright (C) 2018-2026 ABINIT group (LB)
88 : !! This file is distributed under the terms of the
89 : !! GNU General Public License, see ~abinit/COPYING
90 : !! or http://www.gnu.org/copyleft/gpl.txt .
91 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
92 : !!
93 : !! INPUTS
94 : !! atindx(natom)=index table for atoms (see gstate.f)
95 : !! cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = array for planewave
96 : !! coefficients of wavefunctions
97 : !! cg1 = first derivative of cg with respect the perturbation i1pert
98 : !! cg2 = first derivative of cg with respect the perturbation i2pert
99 : !! cg3 = first derivative of cg with respect the perturbation i3pert
100 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
101 : !! if 2, COMPLEX
102 : !! dtfil <type(datafiles_type)>=variables related to files
103 : !! dtset <type(dataset_type)>=all input variables for this dataset
104 : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
105 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
106 : !! k3xc(nfftf,nk3xc)=third-order exchange-correlation kernel
107 : !! indsy1(4,nsym1,natom)=indirect indexing array for atom labels
108 : !! i1dir,i2dir,i3dir=directions of the corresponding perturbations
109 : !! i1pert,i2pert,i3pert = type of perturbation that has to be computed
110 : !! kg(3,mpw*mkmem_rbz)=reduced planewave coordinates
111 : !! mband = maximum number of bands
112 : !! mgfft=maximum size of 1D FFTs
113 : !! mkmem_rbz = maximum number of k points which can fit in core memory
114 : !! mk1mem = maximum number of k points for first-order WF
115 : !! which can fit in core memory
116 : !! mpert =maximum number of ipert
117 : !! mpi_enreg=MPI-parallelisation information
118 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
119 : !! mpw = maximum number of planewaves in basis sphere (large number)
120 : !! natom = number of atoms in unit cell
121 : !! nattyp(ntypat)= # atoms of each type.
122 : !! nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
123 : !! nfftotf=total number of real space fine grid points
124 : !! ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid (see NOTES in respfn.F90)
125 : !! nkpt = number of k points
126 : !! nk3xc=second dimension of the array k3xc
127 : !! nspden = number of spin-density components
128 : !! nspinor = number of spinorial components of the wavefunctions
129 : !! nsppol = number of channels for spin-polarization (1 or 2)
130 : !! nsym1=number of symmetry elements in space group consistent with the perturbation
131 : !! npwarr(nkpt) = array holding npw for each k point
132 : !! occ(mband*nkpt*nsppol) = occupation number for each band and k
133 : !! pawang <type(pawang_type)>=paw angular mesh and related data
134 : !! pawang1 <type(pawang_type)>=pawang datastructure containing only the symmetries preserving the perturbation
135 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
136 : !! pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
137 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
138 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
139 : !! pawrhoij0(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
140 : !! pawrhoij1_i1pert(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data (i1pert)
141 : !! pawrhoij1_i2pert(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data (i2pert)
142 : !! pawrhoij1_i3pert(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data (i3pert)
143 : !! paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh
144 : !! paw_an1_i2pert(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh (i2pert)
145 : !! paw_ij1_i2pert(natom) <type(paw_ij_type)>=1st-order paw arrays given on (i,j) channels (i2pert)
146 : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
147 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
148 : !! rho1r1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i1pert)
149 : !! rho1r2(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i2pert)
150 : !! rho1r3(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i3pert)
151 : !! rprimd(3,3) = dimensional primitive translations (bohr)
152 : !! symaf1(nsym)=(anti)ferromagnetic part of symmetry operations
153 : !! symrc1(3,3,nsym)=symmetries of group in terms of operations on reciprocal space primitive translations
154 : !! ucvol=volume of the unit cell
155 : !! vtrial(nfftf,nspden)=GS Vtrial(r).
156 : !! vhartr1_i2pert(cplex*nfftf,nspden)=firs-order hartree potential
157 : !! vtrial1_i2pert(cplex*nfft,nspden)=firs-order local potential
158 : !! vxc1_i2pert(cplex*nfft,nspden)=firs-order exchange-correlation potential
159 : !! ddk_f = wf files
160 : !! xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc (i1pert)
161 : !! xccc3d2(cplex*n3xccc)=3D change in core charge density, see n3xccc (i2pert)
162 : !! xccc3d3(cplex*n3xccc)=3D change in core charge density, see n3xccc (i3pert)
163 : !! xred(3,natom) = reduced atomic coordinates
164 : !!
165 : !! OUTPUT
166 : !! d3etot(2,3,mpert,3,mpert,3,mpert) = third derivatives of the energy tensor
167 : !! = \sum_{i=1}^9 d3etot_i
168 : !! d3etot_1(2,3,mpert,3,mpert,3,mpert) = 1st term of d3etot
169 : !! d3etot_2(2,3,mpert,3,mpert,3,mpert) = 2nd term of d3etot
170 : !! d3etot_3(2,3,mpert,3,mpert,3,mpert) = 3rd term of d3etot
171 : !! d3etot_4(2,3,mpert,3,mpert,3,mpert) = 4th term of d3etot
172 : !! d3etot_5(2,3,mpert,3,mpert,3,mpert) = 5th term of d3etot
173 : !! d3etot_6(2,3,mpert,3,mpert,3,mpert) = 6th term of d3etot
174 : !! d3etot_7(2,3,mpert,3,mpert,3,mpert) = 7th term of d3etot
175 : !! d3etot_8(2,3,mpert,3,mpert,3,mpert) = 8th term of d3etot
176 : !! d3etot_9(2,3,mpert,3,mpert,3,mpert) = 9th term of d3etot
177 : !!
178 : !! SIDE EFFECTS
179 : !! TO DO!
180 : !!
181 : !! SOURCE
182 :
183 291 : subroutine dfptnl_pert(atindx,cg,cg1,cg2,cg3,cplex,dtfil,dtset,d3etot,eigen0,gs_hamkq,k3xc,indsy1,i1dir,i2dir,i3dir,&
184 291 : & i1pert,i2pert,i3pert,kg,mband,mgfft,mkmem_rbz,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,nattyp,nfftf,nfftotf,ngfftf,nkpt,nk3xc,&
185 291 : & nspden,nspinor,nsppol,nsym1,npwarr,occ,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawtab,&
186 291 : & pawrhoij0,pawrhoij1_i1pert,pawrhoij1_i2pert,pawrhoij1_i3pert,&
187 291 : & paw_an0,paw_an1_i2pert,paw_ij1_i2pert,ph1d,psps,rho1r1,rho2r1,rho3r1,rprimd,symaf1,symrc1,&
188 291 : & ucvol,vtrial,vhartr1_i2pert,vtrial1_i2pert,vxc1_i2pert,ddk_f,xccc3d1,xccc3d2,xccc3d3,xred,&
189 291 : & d3etot_1,d3etot_2,d3etot_3,d3etot_4,d3etot_5,d3etot_6,d3etot_7,d3etot_8,d3etot_9)
190 :
191 : !Arguments ------------------------------------
192 : !scalars
193 : integer,intent(in) :: cplex,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband,mgfft
194 : integer,intent(in) :: mk1mem,mkmem_rbz,mpert,mpsang,mpw,natom,nfftf,nfftotf,nkpt,nspden,nsym1
195 : integer,intent(in) :: nk3xc,nspinor,nsppol
196 : real(dp),intent(in) :: ucvol
197 : type(MPI_type),intent(inout) :: mpi_enreg
198 : type(datafiles_type),intent(in) :: dtfil
199 : type(dataset_type),intent(in) :: dtset
200 : type(pseudopotential_type),intent(in) :: psps
201 : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
202 : type(pawang_type),intent(inout) :: pawang,pawang1
203 : type(pawfgr_type),intent(in) :: pawfgr
204 : type(wfk_t),intent(inout) :: ddk_f(5)
205 :
206 : !arrays
207 : integer,intent(in) :: atindx(natom),kg(3,mpw*mkmem_rbz),nattyp(psps%ntypat),ngfftf(18),npwarr(nkpt)
208 : integer,intent(in) :: indsy1(4,nsym1,dtset%natom),symaf1(nsym1),symrc1(3,3,nsym1)
209 : real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol)
210 : real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol)
211 : real(dp),intent(in) :: cg2(2,mpw*nspinor*mband*mk1mem*nsppol)
212 : real(dp),intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol)
213 : real(dp),intent(in) :: eigen0(dtset%mband*dtset%nkpt*dtset%nsppol)
214 : real(dp),intent(in) :: k3xc(nfftf,nk3xc)
215 : real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom)
216 : real(dp),intent(in) :: rho1r1(cplex*nfftf,dtset%nspden),rho2r1(cplex*nfftf,dtset%nspden)
217 : real(dp),intent(in) :: rho3r1(cplex*nfftf,dtset%nspden),rprimd(3,3)
218 : real(dp),intent(in) :: vtrial(cplex*nfftf,nspden)
219 : real(dp),intent(in) :: xccc3d1(cplex*nfftf),xccc3d2(cplex*nfftf),xccc3d3(cplex*nfftf),xred(3,natom)
220 : real(dp),intent(in) :: vxc1_i2pert(cplex*nfftf,nspden),vhartr1_i2pert(cplex*nfftf)
221 : real(dp),intent(inout) :: vtrial1_i2pert(cplex*nfftf,nspden),d3etot(2,3,mpert,3,mpert,3,mpert)
222 : real(dp),intent(inout) :: d3etot_1(2,3,mpert,3,mpert,3,mpert)
223 : real(dp),intent(inout) :: d3etot_2(2,3,mpert,3,mpert,3,mpert)
224 : real(dp),intent(inout) :: d3etot_3(2,3,mpert,3,mpert,3,mpert)
225 : real(dp),intent(inout) :: d3etot_4(2,3,mpert,3,mpert,3,mpert)
226 : real(dp),intent(inout) :: d3etot_5(2,3,mpert,3,mpert,3,mpert)
227 : real(dp),intent(inout) :: d3etot_6(2,3,mpert,3,mpert,3,mpert)
228 : real(dp),intent(inout) :: d3etot_7(2,3,mpert,3,mpert,3,mpert)
229 : real(dp),intent(inout) :: d3etot_8(2,3,mpert,3,mpert,3,mpert)
230 : real(dp),intent(inout) :: d3etot_9(2,3,mpert,3,mpert,3,mpert)
231 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(natom*psps%usepaw)
232 : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
233 : type(pawrhoij_type),intent(in) :: pawrhoij0(natom*psps%usepaw)
234 : type(pawrhoij_type),intent(in),target :: pawrhoij1_i1pert(natom*psps%usepaw)
235 : type(pawrhoij_type),intent(in) :: pawrhoij1_i2pert(natom*psps%usepaw)
236 : type(pawrhoij_type),intent(in),target :: pawrhoij1_i3pert(natom*psps%usepaw)
237 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
238 : type(paw_an_type),intent(in) :: paw_an0(natom*psps%usepaw)
239 : type(paw_an_type),intent(inout) :: paw_an1_i2pert(natom*psps%usepaw)
240 : type(paw_ij_type),intent(inout) :: paw_ij1_i2pert(natom*psps%usepaw)
241 :
242 : !Local variables-------------------------------
243 : !scalars
244 : logical :: has_cprj_jband,compute_conjugate,compute_rho21
245 : integer,parameter :: level=52,tim_nonlop=0
246 : integer :: bandtot,choice,counter,cplex_cprj,cplex_loc,cplex_rhoij,cpopt,dimffnl1,iband,icg0,ider,ierr
247 : integer :: idir0,idir_getgh2c,idir_phon,idir_elfd,ipert_phon,ipert_elfd
248 : integer :: ia,iatm,ibg,ii,ikg,ikg1,ikpt,ilm,isppol,istwf_k,jband
249 : integer :: me,n1,n2,n3,n4,n5,n6,nband_k,nkpg,nkpg1,nnlout,nsp,nspden_rhoij,npert_phon,npw_k,npw1_k,nzlmopt
250 : integer :: offset_cgi,offset_cgj,offset_eig0,option,paw_opt,qphase_rhoij,debug_mode
251 : integer :: signs,size_wf,size_cprj,spaceComm,typat_ipert_phon,usepaw,useylmgr1
252 : real(dp) :: arg,dot1i,dot1r,dot2i,dot2r,doti,dotr,e3tot,lagi,lagi_paw,lagr,lagr_paw
253 : real(dp) :: sumi,sum_psi1H1psi1,sum_psi1H1psi1_i
254 : real(dp) :: sum_lambda1psi1psi1,sum_lambda1psi1psi1_i
255 : real(dp) :: sum_psi0H2psi1a,sum_psi0H2psi1a_i,sum_psi0H2psi1b,sum_psi0H2psi1b_i
256 : real(dp) :: sum_lambda1psi0S1psi1,sum_lambda1psi0S1psi1_i
257 : character(len=1000) :: msg
258 : !arrays
259 291 : integer,allocatable :: kg_k(:,:),kg1_k(:,:)
260 582 : real(dp) :: buffer(10),eHxc21_paw(2),eHxc21_nhat(2),exc3(2),exc3_paw(2),kpt(3),eig0_k(mband)
261 : real(dp) :: enlout1(2),enlout2(2)
262 : real(dp) :: rmet(3,3),wtk_k
263 291 : real(dp),allocatable :: cgi(:,:),cgj(:,:),cg_jband(:,:,:),cwavef1(:,:),cwavef2(:,:),cwavef3(:,:),dkinpw(:)
264 291 : real(dp),allocatable :: eig1_k_i2pert(:),eig1_k_stored(:)
265 291 : real(dp),allocatable :: chi_ij(:,:,:,:),cwave_right(:,:),cwave_left(:,:),dudk(:,:),dudkde(:,:),dummy_array(:),dummy_array2(:,:)
266 291 : real(dp),allocatable :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:)
267 291 : real(dp),allocatable :: h_cwave(:,:),iddk(:,:),kinpw1(:),kpg_k(:,:),kpg1_k(:,:),nhat21(:,:),occ_k(:)
268 291 : real(dp),allocatable :: phkxred(:,:),ph3d(:,:,:),s_cwave(:,:)
269 291 : real(dp),allocatable :: vlocal(:,:,:,:),vlocal1_i2pert(:,:,:,:),v_i2pert(:,:),wfraug(:,:,:,:)
270 291 : real(dp),allocatable :: ylm(:,:),ylm1(:,:),ylmgr(:,:,:),ylmgr1(:,:,:)
271 291 : real(dp),allocatable :: ylm_k(:,:),ylm1_k(:,:),ylmgr1_k(:,:,:)
272 291 : type(pawcprj_type),allocatable :: cwaveprj0(:,:),cwaveprj1(:,:)
273 291 : type(pawcprj_type),target :: cprj_empty(0,0)
274 291 : type(pawcprj_type),allocatable,target :: cprj_jband(:,:)
275 291 : type(pawrhoij_type),allocatable,target :: pawrhoij21(:)
276 291 : type(pawrhoij_type),pointer :: pawrhoij21_unsym(:),pawrhoij11(:)
277 291 : type(paw_ij_type),allocatable :: paw_ij_tmp(:)
278 291 : type(rf_hamiltonian_type) :: rf_hamkq_i2pert
279 :
280 : !***********************************************************************
281 :
282 : DBG_ENTER("COLL")
283 :
284 : ABI_UNUSED(dtfil%ireadwf)
285 :
286 :
287 291 : me = mpi_enreg%me
288 291 : spaceComm=mpi_enreg%comm_cell
289 :
290 291 : npert_phon = 0
291 138 : if(i1pert<=natom) npert_phon = npert_phon + 1
292 291 : if(i2pert<=natom) npert_phon = npert_phon + 1
293 291 : if(i3pert<=natom) npert_phon = npert_phon + 1
294 291 : if (npert_phon>1) then
295 0 : ABI_ERROR("dfptnl_pert is available with at most one phonon perturbation. Change your input!")
296 : end if
297 :
298 291 : usepaw = psps%usepaw
299 291 : size_cprj = nspinor
300 :
301 291 : call rf_hamkq_i2pert%init(cplex,gs_hamkq,i2pert,paw_ij1=paw_ij1_i2pert,has_e1kbsc=.true.)
302 :
303 291 : ABI_MALLOC(dummy_array,(0))
304 291 : ABI_MALLOC(dummy_array2,(0,0))
305 :
306 : !Acivate computation of rho^(2:1) and related energy derivatives if needed
307 291 : compute_rho21 = .false.
308 291 : if (usepaw==1.and.npert_phon==1.and.(i1pert<=natom.or.i3pert<=natom)) then ! so i2pert==natom+2
309 120 : compute_rho21 = .true.
310 120 : if (i1pert<=natom) then
311 120 : ipert_phon = i1pert
312 120 : idir_phon = i1dir
313 120 : ipert_elfd = i3pert
314 120 : idir_elfd = i3dir
315 120 : pawrhoij11 => pawrhoij1_i3pert
316 0 : else if (i3pert<=natom) then
317 0 : ipert_phon = i3pert
318 0 : idir_phon = i3dir
319 0 : ipert_elfd = i1pert
320 0 : idir_elfd = i1dir
321 0 : pawrhoij11 => pawrhoij1_i1pert
322 : end if
323 600 : ABI_MALLOC(pawrhoij21,(natom))
324 120 : call pawrhoij_nullify(pawrhoij21)
325 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,&
326 120 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc)
327 : call pawrhoij_alloc(pawrhoij21,cplex_rhoij,nspden_rhoij,nspinor,dtset%nsppol,dtset%typat,&
328 120 : & qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
329 840 : ABI_MALLOC(cwaveprj0,(natom,size_cprj))
330 720 : ABI_MALLOC(cwaveprj1,(natom,size_cprj))
331 120 : call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj)
332 120 : call pawcprj_alloc(cwaveprj1,1,gs_hamkq%dimcprj)
333 : ! if (paral_atom) then
334 : ! ABI_MALLOC(pawrhoij1_unsym,(natom))
335 : ! call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,&
336 : !& nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc)
337 : ! call pawrhoij_alloc(pawrhoij1_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,&
338 : !& dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1)
339 : ! else
340 120 : pawrhoij21_unsym => pawrhoij21
341 120 : call pawrhoij_init_unpacked(pawrhoij21_unsym)
342 : ! end if
343 : ! Compute phkxred :
344 360 : ABI_MALLOC(phkxred,(2,natom))
345 360 : do ia=1,natom
346 240 : iatm=min(atindx(ia),natom)
347 240 : arg=two_pi*(kpt(1)*xred(1,ia)+kpt(2)*xred(2,ia)+kpt(3)*xred(3,ia))
348 360 : phkxred(1,iatm)=cos(arg);phkxred(2,iatm)=sin(arg)
349 : end do
350 600 : ABI_MALLOC(paw_ij_tmp,(natom))
351 120 : call paw_ij_nullify(paw_ij_tmp)
352 120 : cplex_loc=1;nsp=1 ! Force nsppol/nspden to 1 because Dij^(1) due to electric field is spin-independent
353 : call paw_ij_init(paw_ij_tmp,cplex_loc,dtset%nspinor,nsp,nsp,dtset%pawspnorb,natom,psps%ntypat,&
354 120 : & dtset%typat,pawtab,has_dijfr=1,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
355 120 : call paw_ij_reset_flags(paw_ij_tmp,all=.True.)
356 : call pawdijfr(gs_hamkq%gprimd,idir_elfd,ipert_elfd,natom,natom,nfftf,ngfftf,nsp,nsp,psps%ntypat,&
357 : & 1,paw_ij_tmp,pawang,pawfgrtab,pawrad,pawtab,cplex_loc,&
358 : & (/zero,zero,zero/),rprimd,ucvol,dummy_array2,dummy_array2,dummy_array2,xred,&
359 120 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
360 600 : ABI_MALLOC(chi_ij,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,dtset%nspinor**2,cplex_loc))
361 120 : call pawdij2e1kb(paw_ij_tmp,1,mpi_enreg%comm_atom,mpi_enreg%my_atmtab,e1kbfr=chi_ij)
362 120 : call paw_ij_free(paw_ij_tmp)
363 480 : ABI_FREE(paw_ij_tmp)
364 : else
365 171 : ABI_MALLOC(chi_ij,(0,0,0,0))
366 171 : ABI_MALLOC(phkxred,(0,0))
367 171 : ABI_MALLOC(pawrhoij21,(0))
368 171 : pawrhoij21_unsym => pawrhoij21
369 171 : ABI_MALLOC(cwaveprj0,(0,0))
370 171 : ABI_MALLOC(cwaveprj1,(0,0))
371 : end if
372 :
373 291 : nnlout = 0
374 291 : bandtot = 0
375 291 : icg0 = 0
376 291 : option = 2
377 291 : n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
378 291 : n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
379 :
380 1746 : ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
381 1746 : ABI_MALLOC(vlocal1_i2pert,(cplex*n4,n5,n6,gs_hamkq%nvloc))
382 :
383 1455 : ABI_MALLOC(wfraug,(2,n4,n5,n6))
384 :
385 11640 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
386 :
387 291 : sumi = zero
388 :
389 : !Set up the Ylm for each k point
390 1164 : ABI_MALLOC(ylm,(dtset%mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
391 1164 : ABI_MALLOC(ylmgr,(dtset%mpw*mkmem_rbz,9,psps%mpsang*psps%mpsang*psps%useylm))
392 291 : if (psps%useylm==1) then
393 : option=2
394 : call initylmg(gs_hamkq%gprimd,kg,dtset%kptns,mkmem_rbz,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,&
395 246 : dtset%nkpt,npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr)
396 : end if
397 :
398 : !Set up the spherical harmonics (Ylm) at k+q
399 291 : useylmgr1=0; option=0
400 291 : if (psps%useylm==1.and. &
401 : & (i2pert==natom+1.or.i2pert==natom+3.or.i2pert==natom+4.or.(usepaw==1.and.i2pert==natom+2))) then
402 166 : useylmgr1=1; option=1
403 : end if
404 873 : ABI_MALLOC(ylm1,(dtset%mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
405 1164 : ABI_MALLOC(ylmgr1,(dtset%mpw*mkmem_rbz,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
406 : !To change the following when q/=0
407 291 : if (psps%useylm==1) then
408 : call initylmg(gs_hamkq%gprimd,kg,dtset%kptns,mkmem_rbz,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,&
409 246 : dtset%nkpt,npwarr,dtset%nsppol,option,rprimd,ylm1,ylmgr1)
410 : end if
411 :
412 291 : debug_mode = 0
413 291 : if (dtset%nonlinear_info>3.or.dtset%nonlinear_info==2) debug_mode = 1
414 :
415 : !Real parts
416 291 : sum_psi1H1psi1 = zero
417 291 : sum_lambda1psi1psi1 = zero
418 291 : sum_lambda1psi0S1psi1 = zero
419 291 : sum_psi0H2psi1a = zero
420 291 : sum_psi0H2psi1b = zero
421 : !Imaginary parts
422 291 : sum_psi1H1psi1_i = zero
423 291 : sum_lambda1psi1psi1_i = zero
424 291 : sum_lambda1psi0S1psi1_i = zero
425 291 : sum_psi0H2psi1a_i = zero
426 291 : sum_psi0H2psi1b_i = zero
427 :
428 291 : compute_conjugate = .false.
429 : !We have to compute < u^(ip1) | H^(ip2) | u^(ip3) >
430 : !For some cases, we want to apply H^(ip2) on < u^(ip1) |, not on | u^(ip3) > (see below)
431 291 : if (i3pert<=natom) then ! As npert_phon<=1, this implies that i1pert=natom+2 and i2pert=natom+2
432 0 : compute_conjugate = .true.
433 : end if
434 :
435 : !Loop over spins
436 828 : do isppol = 1, nsppol
437 :
438 : ! Set up local potential vlocal1 with proper dimensioning, from vtrial1
439 : ! Same thing for vlocal from vtrial Also take into account the spin.
440 : call rf_transgrid_and_pack(isppol,nspden,usepaw,cplex,nfftf,dtset%nfft,dtset%ngfft,&
441 537 : & gs_hamkq%nvloc,pawfgr,mpi_enreg,vtrial,vtrial1_i2pert,vlocal,vlocal1_i2pert)
442 :
443 : ! Continue to initialize the Hamiltonian
444 537 : call gs_hamkq%load_spin(isppol,vlocal=vlocal,with_nonlocal=.true.)
445 537 : call rf_hamkq_i2pert%load_spin(isppol,vlocal1=vlocal1_i2pert,with_nonlocal=.true.)
446 :
447 : ! Loop over k-points
448 :
449 537 : ikg = 0
450 537 : ikg1 = 0
451 :
452 12204 : do ikpt = 1, nkpt
453 :
454 11376 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,mpi_enreg%me)) then
455 : cycle ! Skip the rest of the k-point loop
456 : end if
457 :
458 11376 : counter = 100*ikpt
459 :
460 11376 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
461 11376 : npw_k = npwarr(ikpt)
462 11376 : npw1_k = npw_k ! To change for q/=0
463 11376 : istwf_k = dtset%istwfk(ikpt)
464 34128 : ABI_MALLOC(occ_k,(nband_k))
465 88368 : occ_k(:) = occ(1+bandtot:nband_k+bandtot)
466 11376 : wtk_k = dtset%wtk(ikpt)
467 :
468 11376 : size_wf = nspinor*npw_k
469 :
470 45504 : kpt(:) = dtset%kptns(:,ikpt)
471 :
472 34128 : ABI_MALLOC(cwavef1,(2,npw_k*nspinor))
473 22752 : ABI_MALLOC(cwavef3,(2,npw_k*nspinor))
474 11376 : if (compute_rho21) then
475 4032 : ABI_MALLOC(cwavef2,(2,npw_k*nspinor))
476 : end if
477 :
478 34128 : ABI_MALLOC(kg_k,(3,npw_k))
479 22752 : ABI_MALLOC(kg1_k,(3,npw1_k))
480 45504 : ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
481 34128 : ABI_MALLOC(ylm1_k,(npw1_k,mpsang*mpsang*psps%useylm))
482 45504 : ABI_MALLOC(ylmgr1_k,(npw1_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
483 :
484 : ! Get (k+G) wave vectors and associated spherical harmonics
485 3297408 : kg_k(:,1:npw_k) = kg(:,1+ikg:npw_k+ikg)
486 11376 : if (psps%useylm==1) then
487 20880 : do ilm=1,mpsang*mpsang
488 1595232 : ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
489 : end do
490 : end if
491 :
492 : ! Get (k+q+G) wave vectors and associated spherical harmonics
493 3297408 : kg1_k(:,1:npw1_k)=kg(:,1+ikg1:npw1_k+ikg1) ! To change for q/=0
494 11376 : if (psps%useylm==1) then
495 20880 : do ilm=1,psps%mpsang*psps%mpsang
496 1595232 : ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm)
497 : end do
498 4176 : if (useylmgr1==1) then
499 14160 : do ilm=1,psps%mpsang*psps%mpsang
500 48144 : do ii=1,3
501 3248304 : ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm)
502 : end do
503 : end do
504 : end if
505 : end if
506 :
507 : ! Compute (k+G) vectors
508 11376 : nkpg=0;if(i2pert>=1.and.i2pert<=natom) nkpg=3*dtset%nloalg(3)
509 45504 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
510 11376 : if (nkpg>0) then
511 1344 : call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k)
512 : end if
513 :
514 : ! Compute (k+q+G) vectors
515 11376 : nkpg1=0;if(i2pert>=1.and.i2pert<=natom) nkpg1=3*dtset%nloalg(3)
516 45504 : ABI_MALLOC(kpg1_k,(npw1_k,nkpg1))
517 11376 : if (nkpg1>0) then
518 1344 : call mkkpg(kg1_k,kpg1_k,kpt,nkpg1,npw1_k)
519 : end if
520 :
521 : ! ===== Preparation of the non-local contributions
522 :
523 : ! Compute nonlocal form factors ffnl1 at (k+q+G)
524 : !-- Atomic displacement perturbation
525 11376 : if (i2pert<=natom) then
526 3264 : ider=0;idir0=0
527 : !-- Electric field perturbation
528 8112 : else if (i2pert==natom+2) then
529 8112 : if (psps%usepaw==1) then
530 2832 : ider=1;idir0=i2dir
531 : else
532 5280 : ider=0;idir0=0
533 : end if
534 : end if
535 11376 : if (compute_rho21) then ! compute_rho21 implies i2pert==natom+2
536 2016 : ider=1; idir0=4
537 : end if
538 :
539 : ! Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms
540 11376 : dimffnl1=1+ider
541 11376 : if (ider==1.and.(idir0==0.or.idir0==4)) dimffnl1=2+2*psps%useylm
542 : ! if (ider==2.and.idir0==4) dimffnl1=3+7*psps%useylm
543 68256 : ABI_MALLOC(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat))
544 : call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,ider,idir0,&
545 : & psps%indlmn,kg1_k,kpg1_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,&
546 11376 : & npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k)
547 :
548 : ! Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms
549 11376 : if (compute_rho21.and.debug_mode/=0) then
550 4320 : ABI_MALLOC(ffnl1_test,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat))
551 864 : idir0 = 0 ! for nonlop with signs = 1
552 : call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1_test,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,ider,idir0,&
553 : & psps%indlmn,kg1_k,kpg1_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,&
554 864 : & npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k)
555 : end if
556 :
557 : ! ===== Preparation of the kinetic contributions
558 :
559 : ! Note that not all these arrays should be allocated in the general case when wtk_k vanishes
560 :
561 : ! Compute (1/2) (2 Pi)**2 (k+q+G)**2:
562 34128 : ABI_MALLOC(kinpw1,(npw1_k))
563 832884 : kinpw1(:)=zero
564 11376 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gs_hamkq%gmet,kg1_k,kinpw1,kpt,npw1_k,0,0)
565 :
566 22752 : ABI_MALLOC(dkinpw,(npw_k)) ! 1st derivative (1st direction)
567 832884 : dkinpw(:)=zero
568 :
569 : !===== Load the k/k+q dependent parts of the Hamiltonian
570 :
571 : ! Load k-dependent part in the Hamiltonian datastructure
572 45504 : ABI_MALLOC(ph3d,(2,npw_k,gs_hamkq%matblk))
573 : call gs_hamkq%load_k(kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,&
574 11376 : & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.)
575 : call gs_hamkq%load_k(ffnl_k=ffnl1,kpt_k=kpt,npw_k=npw1_k,istwf_k=istwf_k,&
576 11376 : & kinpw_k=kinpw1,kg_k=kg1_k,kpg_k=kpg1_k,compute_gbound=.true.)
577 : ! end if
578 :
579 : ! Load k-dependent part in the 1st-order Hamiltonian datastructure
580 11376 : call rf_hamkq_i2pert%load_k(npw_k=npw_k,dkinpw_k=dkinpw)
581 :
582 34128 : ABI_MALLOC_OR_DIE(dudk, (2,nband_k*size_wf), ierr)
583 22752 : ABI_MALLOC_OR_DIE(dudkde,(2,nband_k*size_wf), ierr)
584 34128 : ABI_MALLOC_OR_DIE(eig1_k_i2pert,(2*nband_k), ierr)
585 34128 : ABI_MALLOC(eig1_k_stored,(2*nband_k**2))
586 22752 : ABI_MALLOC(cgi,(2,size_wf))
587 22752 : ABI_MALLOC(cwave_right,(2,size_wf))
588 22752 : ABI_MALLOC(cwave_left,(2,size_wf))
589 :
590 : ! **************************************************************************************************
591 : ! Read dudk and dudkde
592 : ! **************************************************************************************************
593 :
594 88368 : do iband = 1,nband_k
595 88368 : if (occ_k(iband)>tol10) then
596 :
597 : ! Read dude file
598 45504 : call ddk_f(1)%read_bks(iband, ikpt, isppol, xmpio_single, cg_bks=cwave_right,eig1_bks=eig1_k_i2pert)
599 : ! Copy eig1_k_i2pert in "eig1_k_stored"
600 661440 : eig1_k_stored(1+(iband-1)*2*nband_k:2*nband_k+(iband-1)*2*nband_k)=eig1_k_i2pert(:)
601 :
602 45504 : if (i2pert==natom+2) then
603 : ! Read dudk file
604 32448 : call ddk_f(2)%read_bks(iband, ikpt, isppol, xmpio_single, cg_bks=cwave_right,eig1_bks=eig1_k_i2pert)
605 32448 : offset_cgi = (iband-1)*size_wf+icg0
606 7001136 : cgi(:,:) = cg(:,1+offset_cgi:size_wf+offset_cgi)
607 : ! Copy cwave_right in "dudk"
608 7001136 : dudk(:,1+(iband-1)*size_wf:iband*size_wf)=cwave_right(:,:)
609 :
610 : ! Read dudkde file
611 32448 : call ddk_f(3)%read_bks(iband, ikpt, isppol, xmpio_single, cg_bks=cwave_right,eig1_bks=eig1_k_i2pert)
612 32448 : offset_cgi = (iband-1)*size_wf+icg0
613 7001136 : cgi(:,:) = cg(:,1+offset_cgi:size_wf+offset_cgi)
614 : ! Copy cwave_right in "dudkde"
615 7001136 : dudkde(:,1+(iband-1)*size_wf:iband*size_wf)=cwave_right(:,:)
616 : end if
617 :
618 : end if
619 : end do
620 :
621 22752 : ABI_MALLOC(cgj,(2,size_wf))
622 22752 : ABI_MALLOC(iddk,(2,size_wf))
623 :
624 11376 : offset_eig0 = mband*(ikpt-1+nkpt*(isppol-1))
625 88368 : eig0_k(:) = eigen0(1+offset_eig0:mband+offset_eig0)
626 :
627 22752 : ABI_MALLOC_OR_DIE(h_cwave,(2,size_wf), ierr)
628 22752 : ABI_MALLOC_OR_DIE(s_cwave,(2,size_wf), ierr)
629 :
630 : ! Allocate work spaces when debug_mode is activated
631 11376 : has_cprj_jband=.false.
632 11376 : if (debug_mode/=0) then ! Only for test purposes
633 6912 : ABI_MALLOC(cg_jband,(2,size_wf*nband_k,2))
634 5864832 : cg_jband(:,:,1) = cg(:,1+icg0:size_wf*nband_k+icg0)
635 1728 : if (i2pert==natom+2) then ! Note the multiplication by "i"
636 1304064 : cg_jband(1,:,2) = -dudk(2,1:size_wf*nband_k)
637 1304064 : cg_jband(2,:,2) = dudk(1,1:size_wf*nband_k)
638 : end if
639 1728 : if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
640 0 : ABI_MALLOC(cprj_jband,(natom,size_cprj*nband_k))
641 0 : has_cprj_jband=.true.
642 : else
643 3456 : ABI_MALLOC(cprj_jband,(natom,0))
644 : end if
645 : else
646 9648 : ABI_MALLOC(cg_jband,(2,0,2))
647 19296 : ABI_MALLOC(cprj_jband,(natom,0))
648 : end if
649 :
650 : ! Loop over bands
651 88368 : do jband = 1,nband_k
652 : ! Skip bands not treated by current proc
653 76992 : if((mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me)) cycle
654 :
655 88368 : if (occ_k(jband)>tol10) then
656 :
657 : ! tol_test = tol8
658 45504 : offset_cgj = (jband-1)*size_wf+icg0
659 9903600 : cgj(:,:) = cg(:,1+offset_cgj:size_wf+offset_cgj)
660 :
661 : ! **************************************************************************************************
662 : ! Compute < u^(1) | ( H^(1) - eps^(0) S^(1) ) | u^(1) >
663 : ! **************************************************************************************************
664 :
665 661440 : eig1_k_i2pert(:) = eig1_k_stored(1+(jband-1)*2*nband_k:jband*2*nband_k)
666 9903600 : cwavef1(:,:) = cg1(:,1+offset_cgj:size_wf+offset_cgj)
667 9903600 : cwavef3(:,:) = cg3(:,1+offset_cgj:size_wf+offset_cgj)
668 45504 : if (i2pert==natom+2) then ! Note the multiplication by i
669 2355344 : iddk(1,:) = -dudkde(2,1+(jband-1)*size_wf:jband*size_wf)
670 2355344 : iddk(2,:) = dudkde(1,1+(jband-1)*size_wf:jband*size_wf)
671 : else
672 2902464 : iddk(:,:) = zero
673 : end if
674 9903600 : cwave_right(:,:) = cwavef3(:,:)
675 9903600 : cwave_left(:,:) = cwavef1(:,:)
676 45504 : if (compute_conjugate) then
677 0 : cwave_right(:,:) = cwavef1(:,:)
678 0 : cwave_left(:,:) = cwavef3(:,:)
679 : end if
680 :
681 : ! Compute : < u^(ip1) | ( H^(ip2) - eps^(0) S^(ip2) ) | u^(ip3) >
682 : ! or : < u^(ip3) | ( H^(ip2) - eps^(0) S^(ip2) ) | u^(ip1) >
683 : call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_right,cprj_empty,h_cwave,s_cwave,eig0_k,eig1_k_i2pert,&
684 : & jband,gs_hamkq,iddk,i2dir,i2pert,ikpt,isppol,mkmem_rbz,mpi_enreg,nband_k,nsppol,&
685 45504 : & debug_mode,dtset%prtvol,rf_hamkq_i2pert,size_cprj,size_wf)
686 45504 : call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_left,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
687 :
688 45504 : if (usepaw==1.and.i2pert/=natom+2) then ! S^(1) is zero for ipert=natom+2
689 5376 : call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cwave_left,s_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
690 5376 : dotr = dotr - eig0_k(jband)*dot2r
691 5376 : doti = doti - eig0_k(jband)*dot2i
692 : end if
693 45504 : if (compute_conjugate) doti = -doti
694 :
695 : ! **************************************************************************************************
696 : ! Compute sum_i Lambda_ij^(1) < u_j^(1) | u_i^(1)>
697 : ! **************************************************************************************************
698 :
699 : ! eig1_k_i2pert(:) = eig1_k_stored(1+(jband-1)*2*nband_k:jband*2*nband_k)
700 45504 : lagr = zero ; lagi = zero
701 45504 : lagr_paw = zero ; lagi_paw = zero
702 :
703 353472 : do iband = 1, nband_k
704 353472 : if(occ_k(iband)>tol10) then
705 :
706 182016 : offset_cgi = (iband-1)*size_wf+icg0
707 39614400 : cwavef3(:,:) = cg3(:,1+offset_cgi:size_wf+offset_cgi)
708 :
709 182016 : if(debug_mode/=0) then
710 691200 : eig1_k_i2pert(:) = eig1_k_stored(1+(jband-1)*2*nband_k:jband*2*nband_k)
711 : end if
712 :
713 : ! Get Lambda_ij^(1) = < u_i^(0) | H^(1) - eps^(0) S^(1) | u_j^(0) > (see dfpt_cgwf.F90)
714 182016 : dot1r = eig1_k_i2pert(2*iband-1)
715 182016 : dot1i = eig1_k_i2pert(2*iband )
716 :
717 : ! Compute < u_j^(1) | S^(0) | u_i^(1) >
718 182016 : if (usepaw==1) then
719 66816 : ibg = 0
720 : call getgsc(cwavef3,cprj_empty,gs_hamkq,s_cwave,ibg,0,0,ikpt,isppol,&
721 66816 : & size_wf,size_cprj,size_wf,mpi_enreg,1,natom,-1,npw_k,nspinor,select_k=KPRIME_H_KPRIME)
722 : else
723 20655360 : s_cwave(:,:) = cwavef3(:,:)
724 : end if
725 :
726 182016 : call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cwavef1,s_cwave,mpi_enreg%me_g0, mpi_enreg%comm_spinorfft)
727 182016 : lagr = lagr + dot1r*dot2r - dot1i*dot2i
728 182016 : lagi = lagi + dot1r*dot2i + dot1i*dot2r
729 :
730 : ! Compute < u_j^(0) | S^(1) | u_i^(1) >
731 182016 : if (usepaw==1.and.i1pert<=natom) then ! S^(1) is zero for ipert=natom+2
732 :
733 32256 : cpopt=-1+5*gs_hamkq%usecprj ; choice=2 ; signs=2 ; paw_opt=3
734 : call nonlop(choice,cpopt,cprj_empty,dummy_array,gs_hamkq,i1dir,(/zero/),mpi_enreg,1,nnlout,&
735 32256 : & paw_opt,signs,s_cwave,tim_nonlop,cwavef3,dummy_array2,iatom_only=i1pert)
736 32256 : call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cgj,s_cwave,mpi_enreg%me_g0, mpi_enreg%comm_spinorfft)
737 32256 : lagr_paw = lagr_paw + dot1r*dot2r - dot1i*dot2i
738 32256 : lagi_paw = lagi_paw + dot1r*dot2i + dot1i*dot2r
739 :
740 : end if
741 :
742 : ! Compute < u_j^(1) | S^(1) | u_i^(0) >
743 182016 : if (usepaw==1.and.i3pert<=natom) then ! S^(1) is zero for ipert=natom+2
744 :
745 0 : cgi(:,:) = cg(:,1+offset_cgi:size_wf+offset_cgi)
746 0 : cpopt=-1+5*gs_hamkq%usecprj ; choice=2 ; signs=2 ; paw_opt=3
747 : call nonlop(choice,cpopt,cprj_empty,dummy_array,gs_hamkq,i3dir,(/zero/),mpi_enreg,1,nnlout,&
748 0 : & paw_opt,signs,s_cwave,tim_nonlop,cgi,dummy_array2,iatom_only=i3pert)
749 0 : call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cwavef1,s_cwave,mpi_enreg%me_g0, mpi_enreg%comm_spinorfft)
750 0 : lagr_paw = lagr_paw + dot1r*dot2r - dot1i*dot2i
751 0 : lagi_paw = lagi_paw + dot1r*dot2i + dot1i*dot2r
752 :
753 : end if
754 :
755 : end if
756 : end do ! iband
757 :
758 : ! **************************************************************************************************
759 : ! Sum all band_by_band contributions
760 : ! **************************************************************************************************
761 :
762 : ! Real part
763 45504 : sum_psi1H1psi1 = sum_psi1H1psi1 + dtset%wtk(ikpt)*occ_k(jband)*dotr
764 45504 : sum_lambda1psi1psi1 = sum_lambda1psi1psi1 - dtset%wtk(ikpt)*occ_k(jband)*lagr
765 45504 : sum_lambda1psi0S1psi1 = sum_lambda1psi0S1psi1 - dtset%wtk(ikpt)*occ_k(jband)*lagr_paw
766 :
767 : ! Imaginary part
768 45504 : sum_psi1H1psi1_i = sum_psi1H1psi1_i + dtset%wtk(ikpt)*occ_k(jband)*doti
769 45504 : sum_lambda1psi1psi1_i = sum_lambda1psi1psi1_i - dtset%wtk(ikpt)*occ_k(jband)*lagi
770 45504 : sum_lambda1psi0S1psi1_i = sum_lambda1psi0S1psi1_i - dtset%wtk(ikpt)*occ_k(jband)*lagi_paw
771 :
772 : ! **************************************************************************************************
773 : ! If compute_rho21 : accumulate rhoij and compute term with H_KV^(2)
774 : ! **************************************************************************************************
775 :
776 45504 : if (compute_rho21) then
777 :
778 8064 : if (i1pert<=natom) then ! If true, i3pert==natom+2
779 2296224 : cwave_right = cg3(:,1+offset_cgj:size_wf+offset_cgj)
780 0 : else if (i3pert<=natom) then ! If true, i1pert==natom+2
781 0 : cwave_right = cg1(:,1+offset_cgj:size_wf+offset_cgj)
782 : end if
783 8064 : choice = 2
784 8064 : cpopt = 0
785 : call getcprj(choice,cpopt,cgj,cwaveprj0,&
786 : & ffnl1,idir_phon,psps%indlmn,gs_hamkq%istwf_k,kg_k,kpg_k,kpt,psps%lmnmax,&
787 : & mgfft,mpi_enreg,1,natom,nattyp,dtset%ngfft,dtset%nloalg,&
788 8064 : & npw_k,nspinor,psps%ntypat,phkxred,ph1d,ph3d,ucvol,psps%useylm)
789 : call getcprj(choice,cpopt,cwave_right,cwaveprj1,&
790 : & ffnl1,idir_phon,psps%indlmn,gs_hamkq%istwf_k,kg_k,kpg_k,kpt,psps%lmnmax,&
791 : & mgfft,mpi_enreg,1,natom,nattyp,dtset%ngfft,dtset%nloalg,&
792 8064 : & npw_k,nspinor,psps%ntypat,phkxred,ph1d,ph3d,ucvol,psps%useylm)
793 :
794 8064 : cplex_cprj=2;if (gs_hamkq%istwf_k>1) cplex_cprj=1
795 : call paw_dfptnl_accrhoij(atindx,cplex_cprj,cwaveprj0,cwaveprj0,cwaveprj1,cwaveprj1,i1pert,i3pert,isppol,natom,natom,&
796 8064 : & nspinor,occ_k(jband),pawrhoij21_unsym,wtk_k)
797 : !& comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
798 :
799 : ! Compute < psi^(0) | H_KV^(pert1pert3) | psi^(pert2) > + < psi^(pert2) | H_KV^(pert1pert3) | psi^(0) >
800 2288160 : cwavef2(:,:) = cg2(:,1+offset_cgj:size_wf+offset_cgj)
801 :
802 : ! Read dkk file (for tests only)
803 8064 : if (debug_mode/=0) then
804 3456 : if(idir_elfd==i2dir) then
805 1152 : call ddk_f(2)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
806 : else
807 2304 : call ddk_f(4)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
808 : end if
809 329184 : cg_jband(1,1+size_wf*(jband-1):size_wf*jband,2) = -iddk(2,1:size_wf)
810 329184 : cg_jband(2,1+size_wf*(jband-1):size_wf*jband,2) = iddk(1,1:size_wf)
811 : end if
812 : ! Read dkde file
813 8064 : if(idir_elfd==i2dir) then
814 2688 : call ddk_f(3)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
815 : else
816 5376 : call ddk_f(5)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
817 : end if
818 2296224 : s_cwave = iddk
819 768096 : iddk(1,:) = -s_cwave(2,:)
820 768096 : iddk(2,:) = s_cwave(1,:)
821 8064 : call rf2_getidir(idir_phon,idir_elfd,idir_getgh2c)
822 : call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwavef2,cprj_empty,s_cwave,dummy_array2,eig0_k,eig1_k_i2pert,&
823 : & jband,gs_hamkq,iddk,idir_getgh2c,ipert_phon+natom+11,ikpt,isppol,mkmem_rbz,mpi_enreg,nband_k,nsppol,&
824 8064 : & debug_mode,dtset%prtvol,rf_hamkq_i2pert,size_cprj,size_wf,enl=chi_ij,ffnl1=ffnl1,ffnl1_test=ffnl1_test)
825 : call dotprod_g(enlout1(1),enlout1(2),gs_hamkq%istwf_k,npw_k*nspinor,2,cgj,s_cwave,&
826 8064 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
827 8064 : sum_psi0H2psi1a = sum_psi0H2psi1a + half*dtset%wtk(ikpt)*occ_k(jband)*enlout1(1) ! be careful : factor 0.5
828 8064 : sum_psi0H2psi1a_i = sum_psi0H2psi1a_i + half*dtset%wtk(ikpt)*occ_k(jband)*enlout1(2) ! be careful : factor 0.5
829 :
830 : ! Read ddk file
831 8064 : if(idir_elfd==i2dir) then
832 2688 : call ddk_f(2)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
833 : else
834 5376 : call ddk_f(4)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
835 : end if
836 2296224 : s_cwave = iddk
837 768096 : iddk(1,:) = -s_cwave(2,:)
838 768096 : iddk(2,:) = s_cwave(1,:)
839 : call rf2_apply_hamiltonian(cg_jband,cprj_jband,cgj,cprj_empty,s_cwave,dummy_array2,eig0_k,eig1_k_i2pert,&
840 : & jband,gs_hamkq,iddk,idir_getgh2c,ipert_phon+natom+11,ikpt,isppol,mkmem_rbz,mpi_enreg,nband_k,nsppol,&
841 8064 : & debug_mode,dtset%prtvol,rf_hamkq_i2pert,size_cprj,size_wf,enl=chi_ij,ffnl1=ffnl1,ffnl1_test=ffnl1_test)
842 : call dotprod_g(enlout2(1),enlout2(2),gs_hamkq%istwf_k,npw_k*nspinor,2,cwavef2,s_cwave,&
843 8064 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
844 8064 : sum_psi0H2psi1b = sum_psi0H2psi1b + half*dtset%wtk(ikpt)*occ_k(jband)*enlout2(1) ! be careful : factor 0.5
845 8064 : sum_psi0H2psi1b_i = sum_psi0H2psi1b_i + half*dtset%wtk(ikpt)*occ_k(jband)*enlout2(2) ! be careful : factor 0.5
846 :
847 : end if ! end if compute_rho21
848 :
849 : end if
850 : end do ! end loop over jband
851 :
852 : ! **************************************************************************************************
853 : ! END OF BAND LOOP
854 : ! **************************************************************************************************
855 :
856 11376 : ABI_FREE(cgi)
857 11376 : ABI_FREE(cgj)
858 11376 : ABI_FREE(iddk)
859 11376 : ABI_FREE(h_cwave)
860 11376 : ABI_FREE(s_cwave)
861 :
862 11376 : ABI_FREE(cwave_right)
863 11376 : ABI_FREE(cwave_left)
864 11376 : ABI_FREE(eig1_k_i2pert)
865 11376 : ABI_FREE(eig1_k_stored)
866 :
867 11376 : bandtot = bandtot + nband_k
868 11376 : icg0 = icg0 + npw_k*nspinor*nband_k
869 11376 : ikg = ikg + npw_k
870 11376 : ikg1 = ikg1 + npw1_k
871 :
872 11376 : ABI_FREE(cwavef1)
873 11376 : if (compute_rho21) then
874 2016 : ABI_FREE(cwavef2)
875 : end if
876 11376 : ABI_FREE(cwavef3)
877 11376 : ABI_FREE(dkinpw)
878 11376 : ABI_FREE(kg_k)
879 11376 : ABI_FREE(kg1_k)
880 11376 : ABI_FREE(kinpw1)
881 11376 : ABI_FREE(dudk)
882 11376 : ABI_FREE(dudkde)
883 11376 : ABI_FREE(ylm_k)
884 11376 : ABI_FREE(ylm1_k)
885 11376 : ABI_FREE(ylmgr1_k)
886 11376 : ABI_FREE(ffnl1)
887 11376 : if (compute_rho21.and.debug_mode/=0) then
888 864 : ABI_FREE(ffnl1_test)
889 : end if
890 11376 : ABI_FREE(kpg_k)
891 11376 : ABI_FREE(kpg1_k)
892 11376 : ABI_FREE(cg_jband)
893 11376 : ABI_FREE(occ_k)
894 11376 : ABI_FREE(ph3d)
895 11376 : if (has_cprj_jband) call pawcprj_free(cprj_jband)
896 11913 : ABI_FREE(cprj_jband)
897 :
898 : end do ! end loop over k-points
899 :
900 : end do ! end loop over spins
901 :
902 291 : call rf_hamkq_i2pert%free()
903 :
904 : ! **************************************************************************************************
905 : ! GATHER BAND-BY-BAND AND XC CONTRIBUTIONS
906 : ! **************************************************************************************************
907 :
908 : if (xmpi_paral == 1) then
909 :
910 : ! Real parts
911 291 : buffer(1) = sum_psi1H1psi1
912 291 : buffer(2) = sum_lambda1psi1psi1
913 291 : buffer(3) = sum_lambda1psi0S1psi1
914 291 : buffer(4) = sum_psi0H2psi1a
915 291 : buffer(5) = sum_psi0H2psi1b
916 :
917 : ! Imaginary parts
918 291 : buffer(6) = sum_psi1H1psi1_i
919 291 : buffer(7) = sum_lambda1psi1psi1_i
920 291 : buffer(8) = sum_lambda1psi0S1psi1_i
921 291 : buffer(9) = sum_psi0H2psi1a_i
922 291 : buffer(10) = sum_psi0H2psi1b_i
923 :
924 291 : call xmpi_sum(buffer,spaceComm,ierr)
925 :
926 : ! Real parts
927 291 : sum_psi1H1psi1 = buffer(1)
928 291 : sum_lambda1psi1psi1 = buffer(2)
929 291 : sum_lambda1psi0S1psi1 = buffer(3)
930 291 : sum_psi0H2psi1a = buffer(4)
931 291 : sum_psi0H2psi1b = buffer(5)
932 :
933 : ! Imaginary parts
934 291 : sum_psi1H1psi1_i = buffer(6)
935 291 : sum_lambda1psi1psi1_i = buffer(7)
936 291 : sum_lambda1psi0S1psi1_i = buffer(8)
937 291 : sum_psi0H2psi1a_i = buffer(9)
938 291 : sum_psi0H2psi1b_i = buffer(10)
939 :
940 : ! Accumulate PAW occupancies
941 291 : if (compute_rho21) then
942 120 : call pawrhoij_mpisum_unpacked(pawrhoij21_unsym,spaceComm)
943 : end if
944 :
945 : end if
946 :
947 : ! **************************************************************************************************
948 : ! Compute E_xc^(3) (NOTE : E_H^(3) = 0)
949 : ! **************************************************************************************************
950 :
951 291 : call dfptnl_exc3(cplex,exc3,k3xc,mpi_enreg,nk3xc,nfftf,nfftotf,nspden,rho1r1,rho2r1,rho3r1,ucvol,xccc3d1,xccc3d2,xccc3d3)
952 :
953 291 : exc3_paw = zero
954 291 : if (usepaw==1) then
955 :
956 : call paw_dfptnl_energy(exc3_paw,dtset%ixc,natom,natom,psps%ntypat,paw_an0,pawang,dtset%pawprtvol,pawrad,&
957 246 : & pawrhoij1_i1pert,pawrhoij1_i2pert,pawrhoij1_i3pert,pawtab,dtset%pawxcdev,mpi_enreg%my_atmtab,mpi_enreg%comm_atom)
958 :
959 : end if
960 :
961 : ! **************************************************************************************************
962 : ! Compute E_Hxc^(2:1)
963 : ! **************************************************************************************************
964 :
965 291 : eHxc21_paw = zero
966 291 : eHxc21_nhat = zero
967 291 : if (compute_rho21) then
968 :
969 120 : if (pawfgr%nfft/=nfftf) then
970 0 : write(msg,'(2(a,i10))') 'pawfgr%nfft/=nfftf : pawfgr%nfft=',pawfgr%nfft,' nfftf = ',nfftf
971 0 : ABI_ERROR(msg)
972 : end if
973 :
974 : call pawnhatfr(0,idir_phon,ipert_phon,natom,dtset%natom,nspden,psps%ntypat,&
975 120 : & pawang,pawfgrtab,pawrhoij11,pawtab,rprimd)
976 :
977 480 : ABI_MALLOC(nhat21,(cplex*nfftf,nspden))
978 : call pawmkrho(0,arg,cplex,gs_hamkq%gprimd,idir_phon,indsy1,ipert_phon,mpi_enreg,&
979 : & natom,natom,nspden,nsym1,psps%ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
980 : & dtset%pawprtvol,pawrhoij21,pawrhoij21_unsym,pawtab,dtset%qptn,dummy_array2,dummy_array2,&
981 : & dummy_array2,rprimd,symaf1,symrc1,dtset%typat,ucvol,dtset%usewvl,xred,&
982 120 : & pawang_sym=pawang1,pawnhat=nhat21,pawrhoij0=pawrhoij0)
983 :
984 : ! if (paral_atom) then
985 : ! call pawrhoij_free(pawrhoij21_unsym)
986 : ! ABI_FREE(pawrhoij21_unsym)
987 : ! end if
988 120 : nzlmopt = 0
989 : call pawdfptenergy(eHxc21_paw,i2pert,ipert_phon,dtset%ixc,natom,dtset%natom,dtset%ntypat,&
990 : & nzlmopt,nzlmopt,paw_an0,paw_an1_i2pert,paw_ij1_i2pert,pawang,dtset%pawprtvol,&
991 120 : & pawrad,pawrhoij1_i2pert,pawrhoij21,pawtab,dtset%pawxcdev,dtset%xclevel)
992 : !& mpi_atmtab=my_atmtab,comm_atom=my_comm_atom
993 :
994 360 : ABI_MALLOC(v_i2pert,(cplex*nfftf,nspden))
995 699960 : v_i2pert(:,1) = vhartr1_i2pert(:)
996 120 : if(nspden>1) then
997 664962 : v_i2pert(:,2) = vhartr1_i2pert(:)
998 : end if
999 :
1000 120 : typat_ipert_phon = dtset%typat(ipert_phon)
1001 120 : if (debug_mode/=0) then
1002 54 : write(msg,'(2(a,i6))') ' DFPTNL_PERT : pawtab(',typat_ipert_phon,')%usexcnhat = ',pawtab(ipert_phon)%usexcnhat
1003 54 : call wrtout(std_out,msg,'COLL')
1004 : end if
1005 120 : if (pawtab(typat_ipert_phon)%usexcnhat>0) then
1006 0 : v_i2pert(:,:) = v_i2pert(:,:) + vxc1_i2pert(:,:)
1007 : end if
1008 :
1009 : call dotprod_vn(cplex,nhat21,eHxc21_nhat(1),eHxc21_nhat(2),nfftf,nfftotf,nspden,2,v_i2pert,&
1010 120 : & ucvol,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1011 :
1012 120 : ABI_FREE(v_i2pert)
1013 120 : ABI_FREE(nhat21)
1014 :
1015 : end if
1016 :
1017 : ! **************************************************************************************************
1018 : ! ALL TERMS HAVE BEEN COMPUTED
1019 : ! **************************************************************************************************
1020 :
1021 : !Real part
1022 291 : e3tot = sum_psi1H1psi1 + sum_lambda1psi1psi1 + sum_lambda1psi0S1psi1 + sum_psi0H2psi1a + sum_psi0H2psi1b
1023 291 : e3tot = e3tot + half * (eHxc21_paw(1)+eHxc21_nhat(1)) + sixth * (exc3(1) + exc3_paw(1))
1024 :
1025 : !Imaginary part
1026 291 : sumi = sum_psi1H1psi1_i + sum_lambda1psi1psi1_i + sum_lambda1psi0S1psi1_i + sum_psi0H2psi1a_i + sum_psi0H2psi1b_i
1027 291 : sumi = sumi + half * (eHxc21_paw(2)+eHxc21_nhat(2)) + sixth * (exc3(2) + exc3_paw(2))
1028 :
1029 : !Real parts
1030 291 : d3etot(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = e3tot
1031 291 : d3etot_1(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi1H1psi1
1032 291 : d3etot_2(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi1psi1
1033 291 : d3etot_3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi0S1psi1
1034 291 : d3etot_4(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1a
1035 291 : d3etot_5(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1b
1036 291 : d3etot_6(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_paw(1)
1037 291 : d3etot_7(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_nhat(1)
1038 291 : d3etot_8(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3(1)
1039 291 : d3etot_9(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3_paw(1)
1040 : !Imaginary parts
1041 291 : d3etot(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi
1042 291 : d3etot_1(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi1H1psi1_i
1043 291 : d3etot_2(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi1psi1_i
1044 291 : d3etot_3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi0S1psi1_i
1045 291 : d3etot_4(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1a_i
1046 291 : d3etot_5(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1b_i
1047 291 : d3etot_6(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_paw(2)
1048 291 : d3etot_7(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_nhat(2)
1049 291 : d3etot_8(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3(2)
1050 291 : d3etot_9(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3_paw(2)
1051 :
1052 : !Before printing, set small contributions to zero
1053 : !Real parts
1054 291 : if (abs(sum_psi1H1psi1) <tol8) sum_psi1H1psi1 = zero
1055 291 : if (abs(sum_lambda1psi1psi1) <tol8) sum_lambda1psi1psi1 = zero
1056 291 : if (abs(sum_lambda1psi0S1psi1)<tol8) sum_lambda1psi0S1psi1 = zero
1057 291 : if (abs(sum_psi0H2psi1a) <tol8) sum_psi0H2psi1a = zero
1058 291 : if (abs(sum_psi0H2psi1b) <tol8) sum_psi0H2psi1b = zero
1059 291 : if (abs(eHxc21_paw(1)) <tol8) eHxc21_paw(1) = zero
1060 291 : if (abs(eHxc21_nhat(1)) <tol8) eHxc21_nhat(1) = zero
1061 291 : if (abs(exc3(1)) <tol8) exc3(1) = zero
1062 291 : if (abs(exc3_paw(1)) <tol8) exc3_paw(1) = zero
1063 291 : if (abs(e3tot) <tol8) e3tot = zero
1064 :
1065 : !Imaginary parts
1066 291 : if (abs(sum_psi1H1psi1_i) <tol8) sum_psi1H1psi1_i = zero
1067 291 : if (abs(sum_lambda1psi1psi1_i) <tol8) sum_lambda1psi1psi1_i = zero
1068 291 : if (abs(sum_lambda1psi0S1psi1_i)<tol8) sum_lambda1psi0S1psi1_i = zero
1069 291 : if (abs(sum_psi0H2psi1a_i) <tol8) sum_psi0H2psi1a_i = zero
1070 291 : if (abs(sum_psi0H2psi1b_i) <tol8) sum_psi0H2psi1b_i = zero
1071 291 : if (abs(eHxc21_paw(2)) <tol8) eHxc21_paw(2) = zero
1072 291 : if (abs(eHxc21_nhat(2)) <tol8) eHxc21_nhat(2) = zero
1073 291 : if (abs(exc3(2)) <tol8) exc3(2) = zero
1074 291 : if (abs(exc3_paw(2)) <tol8) exc3_paw(2) = zero
1075 291 : if (abs(sumi) <tol8) sumi = zero
1076 :
1077 291 : write(msg,'(2a,3(a,i2,a,i1))') ch10,'NONLINEAR : ',&
1078 582 : ' perts : ',i1pert,'.',i1dir,' / ',i2pert,'.',i2dir,' / ',i3pert,'.',i3dir
1079 291 : call wrtout(std_out,msg,'COLL')
1080 291 : call wrtout(ab_out,msg,'COLL')
1081 291 : if (dtset%nonlinear_info>0) then
1082 : write(msg,'(10(a,2(a,f18.8)),a)') &
1083 108 : ch10,' sum_psi1H1psi1 = ',sum_psi1H1psi1, ',',sum_psi1H1psi1_i,&
1084 108 : ch10,' sum_lambda1psi1psi1 = ',sum_lambda1psi1psi1, ',',sum_lambda1psi1psi1_i,&
1085 108 : ch10,' sum_lambda1psi0S1psi1 = ',sum_lambda1psi0S1psi1, ',',sum_lambda1psi0S1psi1_i,&
1086 108 : ch10,' sum_psi0H2psi1a = ',sum_psi0H2psi1a, ',',sum_psi0H2psi1a_i,&
1087 108 : ch10,' sum_psi0H2psi1b = ',sum_psi0H2psi1b, ',',sum_psi0H2psi1b_i,&
1088 108 : ch10,' eHxc21_paw/2 = ',half*eHxc21_paw(1), ',',half*eHxc21_paw(2),&
1089 108 : ch10,' eHxc21_nhat/2 = ',half*eHxc21_nhat(1), ',',half*eHxc21_nhat(2),&
1090 108 : ch10,' exc3/6 = ',sixth*exc3(1), ',',sixth*exc3(2),&
1091 108 : ch10,' exc3_paw/6 = ',sixth*exc3_paw(1), ',',sixth*exc3_paw(2),&
1092 216 : ch10,' >>>>>>>>>>>>>>> e3tot = ',e3tot, ',',sumi,ch10
1093 108 : call wrtout(std_out,msg,'COLL')
1094 108 : call wrtout(ab_out,msg,'COLL')
1095 : end if
1096 :
1097 291 : if (compute_rho21) then
1098 120 : call pawcprj_free(cwaveprj0)
1099 120 : call pawcprj_free(cwaveprj1)
1100 120 : call pawrhoij_free(pawrhoij21)
1101 : end if
1102 291 : ABI_FREE(chi_ij)
1103 291 : ABI_FREE(dummy_array)
1104 291 : ABI_FREE(dummy_array2)
1105 291 : ABI_FREE(phkxred)
1106 531 : ABI_FREE(cwaveprj0)
1107 531 : ABI_FREE(cwaveprj1)
1108 531 : ABI_FREE(pawrhoij21)
1109 291 : ABI_FREE(ylm)
1110 291 : ABI_FREE(ylm1)
1111 291 : ABI_FREE(ylmgr)
1112 291 : ABI_FREE(ylmgr1)
1113 291 : ABI_FREE(vlocal)
1114 291 : ABI_FREE(vlocal1_i2pert)
1115 291 : ABI_FREE(wfraug)
1116 :
1117 : DBG_EXIT("COLL")
1118 :
1119 582 : end subroutine dfptnl_pert
1120 : !!***
1121 :
1122 : !!****f* ABINIT/dfptnl_exc3
1123 : !! NAME
1124 : !! dfptnl_exc3
1125 : !!
1126 : !! FUNCTION
1127 : !! Compute the third-order xc energy.
1128 : !!
1129 : !! Take into account the contribution of the term
1130 : !!$
1131 : !! \frac{d}{d \lambda}
1132 : !! \frac{\delta^2 E_{Hxc}}{\delta n(r) \delta n(r\prim)}
1133 : !!$
1134 : !! (seventh term of Eq. (110) of X. Gonze, PRA 52, 1096 (1995) [[cite:Gonze1995]]).
1135 : !!
1136 : !! The following is essentially the 4th and the 3rd terms of PRB 71,125107 [[cite:Veithen2005]].
1137 : !!
1138 : !! However, we treat the nonlinear xc core correction in a slightly different way,
1139 : !! in order to keep the symmetry between pert1,pert2 and pert3. It helps for debugging.
1140 : !!
1141 : !! Namely, here we consider Exc as a functional depending on the TOTAL density (core+valence),
1142 : !! so it does not depend explicitely on the perturbation, and the term given above is always zero.
1143 : !! The "lost" terms are recovered adding the derivative of the core densities for EVERY perturbations,
1144 : !! and not only to 'pert2', as in the first case.
1145 : !!
1146 : !! COPYRIGHT
1147 : !! Copyright (C) 2020-2026 ABINIT group (LB)
1148 : !! This file is distributed under the terms of the
1149 : !! GNU General Public License, see ~abinit/COPYING
1150 : !! or http://www.gnu.org/copyleft/gpl.txt .
1151 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
1152 : !!
1153 : !! INPUTS
1154 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
1155 : !! if 2, COMPLEX
1156 : !! k3xc(nfftf,nk3xc)=third-order exchange-correlation kernel
1157 : !! mpi_enreg=MPI-parallelisation information
1158 : !! nk3xc=second dimension of the array k3xc
1159 : !! nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
1160 : !! nfftotf=total number of real space fine grid points
1161 : !! nspden = number of spin-density components
1162 : !! rho1r1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i1pert)
1163 : !! rho1r2(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i2pert)
1164 : !! rho1r3(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i3pert)
1165 : !! ucvol=volume of the unit cell
1166 : !! xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc (i1pert)
1167 : !! xccc3d2(cplex*n3xccc)=3D change in core charge density, see n3xccc (i2pert)
1168 : !! xccc3d3(cplex*n3xccc)=3D change in core charge density, see n3xccc (i3pert)
1169 : !!
1170 : !! OUTPUT
1171 : !! exc3(2)=real and imaginray part of the exchange correlation energy term
1172 : !!
1173 : !! SIDE EFFECTS
1174 : !!
1175 : !! SOURCE
1176 :
1177 698 : subroutine dfptnl_exc3(cplex,exc3,k3xc,mpi_enreg,nk3xc,nfftf,nfftotf,nspden,rho1r1,rho2r1,rho3r1,ucvol,xccc3d1,xccc3d2,xccc3d3)
1178 :
1179 : !Arguments ------------------------------------
1180 : !scalars
1181 : integer,intent(in) :: cplex,nk3xc,nfftf,nfftotf,nspden
1182 : real(dp),intent(in) :: ucvol
1183 : type(MPI_type),intent(inout) :: mpi_enreg
1184 :
1185 : !arrays
1186 : real(dp),intent(in) :: k3xc(nfftf,nk3xc)
1187 : real(dp),intent(in) :: rho1r1(cplex*nfftf,nspden),rho2r1(cplex*nfftf,nspden),rho3r1(cplex*nfftf,nspden)
1188 : real(dp),intent(in) :: xccc3d1(cplex*nfftf),xccc3d2(cplex*nfftf),xccc3d3(cplex*nfftf)
1189 : real(dp),intent(out) :: exc3(2)
1190 :
1191 : !Local variables-------------------------------
1192 : !scalars
1193 : integer :: ifft,ifft_im,ifft_re
1194 : real(dp) :: rho2ur,rho2ui,rho2dr,rho2di,rho3ur,rho3ui,rho3dr,rho3di
1195 : ! character(len=1000) :: msg
1196 : !arrays
1197 698 : real(dp),allocatable :: rho1r1_tot(:,:),xc_tmp(:,:)
1198 :
1199 : !***********************************************************************
1200 :
1201 : DBG_ENTER("COLL")
1202 :
1203 2792 : ABI_MALLOC(xc_tmp,(cplex*nfftf,nspden))
1204 2094 : ABI_MALLOC(rho1r1_tot,(cplex*nfftf,nspden))
1205 698 : if (nspden==1)then
1206 :
1207 338 : if (cplex==1) then
1208 1969730 : do ifft=1,nfftf
1209 1969392 : rho1r1_tot(ifft,1) = rho1r1(ifft,1) + xccc3d1(ifft)
1210 1969392 : rho2ur = rho2r1(ifft,1)+xccc3d2(ifft)
1211 1969392 : rho3ur = rho3r1(ifft,1)+xccc3d3(ifft)
1212 1969730 : xc_tmp(ifft,1)= k3xc(ifft,1)*rho2ur*rho3ur
1213 : end do
1214 : else
1215 0 : do ifft=1,nfftf
1216 0 : ifft_re = 2*ifft-1 ! Real part
1217 0 : ifft_im = 2*ifft ! Imaginary part
1218 0 : rho1r1_tot(ifft_re,1) = rho1r1(ifft_re,1) + xccc3d1(ifft_re)
1219 0 : rho1r1_tot(ifft_im,1) = rho1r1(ifft_im,1) + xccc3d1(ifft_im)
1220 0 : rho2ur = rho2r1(ifft_re,1)+xccc3d2(ifft_re)
1221 0 : rho2ui = rho2r1(ifft_im,1)+xccc3d2(ifft_im)
1222 0 : rho3ur = rho3r1(ifft_re,1)+xccc3d3(ifft_re)
1223 0 : rho3ui = rho3r1(ifft_im,1)+xccc3d3(ifft_im)
1224 0 : xc_tmp(ifft_re,1)= k3xc(ifft,1)*(rho2ur*rho3ur-rho2ui*rho3ui)
1225 0 : xc_tmp(ifft_im,1)= k3xc(ifft,1)*(rho2ur*rho3ui+rho2ui*rho3ur)
1226 : end do
1227 : end if
1228 :
1229 360 : else if (nspden==2) then
1230 :
1231 : ! Remember : rhor(...,1) = total density ( n_up(r) + n_down(r) )
1232 : ! rhor(...,2) = up density ( n_up(r) )
1233 : ! But : pot(...,1) = up potential ( v_up(r) )
1234 : ! pot(...,2) = down potential ( v_down(r) )
1235 :
1236 360 : if (cplex==1) then
1237 2109024 : do ifft=1,nfftf
1238 2108664 : rho1r1_tot(ifft,1) = rho1r1(ifft,1) + xccc3d1(ifft) ! 1 tot
1239 2108664 : rho1r1_tot(ifft,2) = rho1r1(ifft,2) + xccc3d1(ifft)*half ! 1 up
1240 2108664 : rho2ur = rho2r1(ifft,2) + xccc3d2(ifft)*half ! 2 up
1241 2108664 : rho2dr = rho2r1(ifft,1)-rho2r1(ifft,2) + xccc3d2(ifft)*half ! 2 down
1242 2108664 : rho3ur = rho3r1(ifft,2) + xccc3d3(ifft)*half ! 3 up
1243 2108664 : rho3dr = rho3r1(ifft,1)-rho3r1(ifft,2) + xccc3d3(ifft)*half ! 3 down
1244 : ! uuu uud
1245 : xc_tmp(ifft,1)= k3xc(ifft,1)*rho2ur*rho3ur + k3xc(ifft,2)*rho2ur*rho3dr + &
1246 : ! udu udd
1247 2108664 : & k3xc(ifft,2)*rho2dr*rho3ur + k3xc(ifft,3)*rho2dr*rho3dr
1248 : ! duu dud
1249 : xc_tmp(ifft,2)= k3xc(ifft,2)*rho2ur*rho3ur + k3xc(ifft,3)*rho2ur*rho3dr + &
1250 : ! ddu ddd
1251 2109024 : & k3xc(ifft,3)*rho2dr*rho3ur + k3xc(ifft,4)*rho2dr*rho3dr
1252 : end do
1253 :
1254 : else ! cplex = 2
1255 :
1256 0 : do ifft=1,nfftf
1257 0 : ifft_re = 2*ifft-1 ! Real part
1258 0 : ifft_im = 2*ifft ! Imaginary part
1259 0 : rho1r1_tot(ifft_re,1) = rho1r1(ifft_re,1) + xccc3d1(ifft_re) ! 1 tot re
1260 0 : rho1r1_tot(ifft_im,1) = rho1r1(ifft_im,1) + xccc3d1(ifft_im) ! 1 tot im
1261 0 : rho1r1_tot(ifft_re,2) = rho1r1(ifft_re,2) + xccc3d1(ifft_re)*half ! 1 up re
1262 0 : rho1r1_tot(ifft_im,2) = rho1r1(ifft_im,2) + xccc3d1(ifft_im)*half ! 1 up im
1263 0 : rho2ur = rho2r1(ifft_re,2) + xccc3d2(ifft_re)*half ! 2 up re
1264 0 : rho2ur = rho2r1(ifft_im,2) + xccc3d2(ifft_im)*half ! 2 up im
1265 0 : rho2dr = rho2r1(ifft_re,1)-rho2r1(ifft_re,2) + xccc3d2(ifft_re)*half ! 2 down re
1266 0 : rho2dr = rho2r1(ifft_im,1)-rho2r1(ifft_im,2) + xccc3d2(ifft_im)*half ! 2 down im
1267 0 : rho3ur = rho3r1(ifft_re,2) + xccc3d3(ifft_re)*half ! 3 up re
1268 0 : rho3ur = rho3r1(ifft_im,2) + xccc3d3(ifft_im)*half ! 3 up im
1269 0 : rho3dr = rho3r1(ifft_re,1)-rho3r1(ifft_re,2) + xccc3d3(ifft_re)*half ! 3 down re
1270 0 : rho3dr = rho3r1(ifft_im,1)-rho3r1(ifft_im,2) + xccc3d3(ifft_im)*half ! 3 down im
1271 : ! Real part:
1272 : ! uuu uud
1273 : xc_tmp(ifft_re,1)= k3xc(ifft,1)*(rho2ur*rho3ur-rho2ui*rho3ui) + k3xc(ifft,2)*(rho2ur*rho3dr-rho2ui*rho3di) + &
1274 : ! udu udd
1275 0 : & k3xc(ifft,2)*(rho2dr*rho3ur-rho2di*rho3ui) + k3xc(ifft,3)*(rho2dr*rho3dr-rho2di*rho3di)
1276 : ! duu dud
1277 : xc_tmp(ifft_re,2)= k3xc(ifft,2)*(rho2ur*rho3ur-rho2ui*rho3ui) + k3xc(ifft,3)*(rho2ur*rho3dr-rho2ui*rho3di) + &
1278 : ! ddu ddd
1279 0 : & k3xc(ifft,3)*(rho2dr*rho3ur-rho2di*rho3ui) + k3xc(ifft,4)*(rho2dr*rho3dr-rho2di*rho3di)
1280 : ! Imaginary part:
1281 : ! uuu uud
1282 : xc_tmp(ifft_im,1)= k3xc(ifft,1)*(rho2ur*rho3ui+rho2ui*rho3ur) + k3xc(ifft,2)*(rho2ur*rho3di+rho2ui*rho3dr) + &
1283 : ! udu udd
1284 0 : & k3xc(ifft,2)*(rho2dr*rho3ui+rho2di*rho3ur) + k3xc(ifft,3)*(rho2dr*rho3di+rho2di*rho3dr)
1285 : ! duu dud
1286 : xc_tmp(ifft_im,2)= k3xc(ifft,2)*(rho2ur*rho3ui+rho2ui*rho3ur) + k3xc(ifft,3)*(rho2ur*rho3di+rho2ui*rho3dr) + &
1287 : ! ddu ddd
1288 0 : & k3xc(ifft,3)*(rho2dr*rho3ui+rho2di*rho3ur) + k3xc(ifft,4)*(rho2dr*rho3di+rho2di*rho3dr)
1289 : end do
1290 :
1291 : end if
1292 :
1293 : else
1294 0 : ABI_BUG('DFPTNL_PERT is implemented only for nspden=1 or 2')
1295 : end if
1296 :
1297 698 : call dotprod_vn(cplex,rho1r1_tot,exc3(1),exc3(2),nfftf,nfftotf,nspden,2,xc_tmp,ucvol,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1298 698 : ABI_FREE(xc_tmp)
1299 698 : ABI_FREE(rho1r1_tot)
1300 :
1301 : DBG_EXIT("COLL")
1302 :
1303 698 : end subroutine dfptnl_exc3
1304 : !!***
1305 :
1306 : end module m_dfptnl_pert
1307 : !!***
|
