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