Line data Source code
1 : !!****m* ABINIT/m_dfptlw_loop
2 : !! NAME
3 : !! m_dfptlw_loop
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2022-2026 ABINIT group (MR)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_dfptlw_loop
22 :
23 : use defs_basis
24 : use defs_wvltypes
25 : use m_errors
26 : use m_abicore
27 : use m_hdr
28 : use m_nctk
29 : use m_wffile
30 : use m_wfk
31 : use m_dtset
32 : use m_dtfil
33 :
34 : use defs_datatypes, only : pseudopotential_type
35 : use defs_abitypes, only : MPI_type
36 : use m_time, only : timab
37 : use m_io_tools, only : file_exists,iomode_from_fname,get_unit
38 : use m_kg, only : getmpw,kpgio
39 : use m_inwffil, only : inwffil
40 : use m_fft, only : fourdp
41 : use m_ioarr, only : read_rhor
42 : use m_hamiltonian, only : gs_hamiltonian_type
43 : use m_pawdij, only : pawdij, pawdijfr, symdij
44 : use m_pawfgr, only : pawfgr_type
45 : use m_pawfgrtab, only : pawfgrtab_type
46 : use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags
47 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags, paw_ij_print
48 : use m_pawang, only : pawang_type
49 : use m_pawrad, only : pawrad_type
50 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_nullify, &
51 : & pawrhoij_io, pawrhoij_inquire_dim
52 : use m_paw_nhat, only : pawmknhat,pawnhatfr
53 : use m_paw_denpot, only : pawdenpot
54 : use m_pawtab, only : pawtab_type
55 : use m_rf2, only : rf2_getidir
56 : use m_initylmg, only : initylmg
57 : use m_atm2fft, only : dfpt_atm2fft
58 : use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxcccdq
59 : use m_dfpt_rhotov, only : dfpt_rhotov
60 : use m_mkcore, only : dfpt_mkcore
61 : use m_mklocl, only : dfpt_vlocal, vlocalstr,dfpt_vlocaldq,dfpt_vlocaldqdq,dfpt_vmetdqdq
62 : use m_dfptlw_pert, only : dfptlw_pert
63 : use m_dfpttd_berrycurv, only : dfpttd_berrycurv
64 : use m_dynmat, only : cart39
65 : use m_ebands , only : ebands_t
66 : use m_xmpi
67 :
68 : implicit none
69 :
70 : private
71 : !!***
72 :
73 : public :: dfptlw_loop
74 : !!***
75 : ! *************************************************************************
76 :
77 : contains
78 : !!***
79 :
80 : !!****f* ABINIT/m_dfptlw_loop/dfptlw_loop
81 : !! NAME
82 : !! dfptlw_loop
83 : !!
84 : !! FUNCTION
85 : !! Loop over two perturbations j1, j2 and a q gradient
86 : !!
87 : !! INPUTS
88 : !! atindx(natom)=index table for atoms (see gstate.f)
89 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions
90 : !! codvsn=code version
91 : !! d3e_pert1(mpert)=array with the i1pert cases to calculate
92 : !! d3e_pert2(mpert)=array with the i2pert cases to calculate
93 : !! dimffnl= third dimension of ffnl
94 : !! dtfil <type(datafiles_type)>=variables related to files
95 : !! dtset <type(dataset_type)>=all input variables for this dataset
96 : !! ffnl(dtset%mkmem,dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives
97 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
98 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
99 : !! kg(3,mpw*mkmem)=reduced planewave coordinates
100 : !! kxc(nfftf,nkxc)=exchange-correlation kernel
101 : !! just_timdisp= it .true. this is a calculation of time-dispersion only
102 : !! mband = maximum number of bands
103 : !! mgfft = maximum single fft dimension
104 : !! mkmem = Number of k points treated by this node.
105 : !! mk1mem = Number of k points for first-order WF treated by this node.
106 : !! mpert =maximum number of ipert
107 : !! mpi_enreg=MPI-parallelisation information
108 : !! mpw = maximum number of planewaves in basis sphere (large number)
109 : !! natom = number of atoms in unit cell
110 : !! nattyp(ntypat)= # atoms of each type.
111 : !! ncorespl(mqgrid_vl,2,ntypat)= spl interpolation of the pseudo core density for all atom types
112 : !! nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
113 : !! ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid (see NOTES in respfn.F90)
114 : !! nkpt = number of k points
115 : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
116 : !! nspinor = number of spinorial components of the wavefunctions
117 : !! nsppol = number of channels for spin-polarization (1 or 2)
118 : !! npwarr(nkpt) = array holding npw for each k point
119 : !! nylmgr=second dimension of ylmgr_k
120 : !! occ(mband*nkpt*nsppol) = occupation number for each band and k
121 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
122 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
123 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
124 : !! rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
125 : !! that have to be computed
126 : !! 1 -> element has to be computed explicitely
127 : !! -1 -> use symmetry operations to obtain the corresponding element
128 : !! rhog(2,nfftf)=array for Fourier transform of GS electron density
129 : !! rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3.
130 : !! rmet(3,3)=real space metric tensor in bohr**2
131 : !! rprimd(3,3)=dimensional primitive translations (bohr)
132 : !! ucvol = unit cell volume (bohr^3)
133 : !! useylmgr= if 1 use the derivative of spherical harmonics
134 : !! xred(3,natom) = reduced atomic coordinates
135 : !! ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics
136 : !! ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics
137 : !!
138 : !! OUTPUT
139 : !! blkflg(3,mpert,3,mpert) = flags for each element of the 3DTE
140 : !! (=1 if computed)
141 : !! d3etot(2,3,mpert,3,mpert,3,mpert) = third derivatives of the energy tensor
142 : !!
143 : !! SIDE EFFECTS
144 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
145 : !!
146 : !! SOURCE
147 :
148 17 : subroutine dfptlw_loop(atindx,blkflg,cg,codvsn,d3e_pert1,d3e_pert2,d3etot,dimffnl,dtfil,dtset,&
149 17 : & ffnl,gmet,gprimd,gsqcut,&
150 17 : & hdr,just_timdisp,kg,kxc,mband,mgfft,mkmem,mk1mem,&
151 17 : & mpert,mpi_enreg,mpw,natom,nattyp,ncorespl,ngfftf,nfftf,nkpt,nkxc,nspinor,nsppol,&
152 17 : & npwarr,nylmgr,occ,&
153 17 : & pawfgr,pawtab,ph1d,&
154 17 : & psps,rfpert,rhog,rhor,rmet,rprimd,ucvol,useylmgr,xred,ylm,ylmgr)
155 :
156 :
157 : !Arguments ------------------------------------
158 : !scalars
159 : integer,intent(in) :: dimffnl,mband,mgfft,mk1mem,mkmem,mpert,mpw,natom,nfftf
160 : integer,intent(in) :: nkpt,nkxc,nspinor,nsppol,nylmgr,useylmgr
161 : character(len=8),intent(in) :: codvsn
162 : logical,intent(in) :: just_timdisp
163 : real(dp),intent(in) :: gsqcut,ucvol
164 : type(MPI_type),intent(inout) :: mpi_enreg
165 : type(datafiles_type),intent(in) :: dtfil
166 : type(dataset_type),intent(in) :: dtset
167 : type(hdr_type),intent(inout) :: hdr
168 : type(pawfgr_type),intent(in) :: pawfgr
169 : type(pseudopotential_type),intent(in) :: psps
170 :
171 : !arrays
172 : integer,intent(in) :: atindx(natom),d3e_pert1(mpert),d3e_pert2(mpert)
173 : integer,intent(in) :: kg(3,mk1mem*mpw)
174 : integer,intent(in) :: nattyp(psps%ntypat),ngfftf(18),npwarr(nkpt)
175 : integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert)
176 : integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert)
177 : real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gmet(3,3)
178 : real(dp),intent(in) :: ffnl(mkmem,mpw,dimffnl,psps%lmnmax,psps%ntypat)
179 : real(dp),intent(in) :: gprimd(3,3),kxc(nfftf,nkxc)
180 : real(dp),intent(in) :: ncorespl(psps%mqgrid_vl,2,psps%ntypat)
181 : real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom)
182 : real(dp),intent(in) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden),rmet(3,3),rprimd(3,3)
183 : real(dp),intent(in) :: xred(3,natom)
184 : real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
185 : real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert)
186 : real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
187 : real(dp),intent(in) :: ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)
188 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
189 :
190 : !Local variables-------------------------------
191 : !scalars
192 : integer :: alpha,ask_accurate,bantot,beta,comm_cell,cplex
193 : integer :: delta,dkdk_index,formeig,gamma
194 : integer :: ia1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,idir_dkdk
195 : integer :: idq,ierr,ii,ikpt,ireadwf,istr,itypat,mcg1,me,mpsang
196 : integer :: mpw1, mpw1_mq
197 : integer :: n1,n2,n3,ndir,n1dq,n2dq,nhat1grdim,nfftotf,nspden,n3xccc
198 : integer :: optgeom,opthartdqdq,optorth,pawread
199 : integer :: pert1case,pert2case,pert3case,timrev,usexcnhat
200 : integer :: pert1case_mq,pert2case_mq
201 : real(dp) :: delad,delag,delbd,delbg,ecut,ecut_eff
202 : logical :: finite_q,kramers_deg,samepert
203 : character(len=500) :: message
204 : character(len=fnlen) :: fiden1i,fiwf1i,fiwf2i,fiwfddk,fiwfdkdk
205 : character(len=fnlen) :: fiwf1i_mq,fiwf2i_mq
206 17 : type(gs_hamiltonian_type) :: gs_hamkq
207 : type(wffile_type) :: wff1,wff2,wfft1,wfft2
208 187 : type(wfk_t) :: ddk_f,d2_dkdk_f,d2_dkdk_f2
209 17 : type(wvl_data) :: wvl
210 17 : type(hdr_type) :: hdr_den
211 17 : type(ebands_t) :: ebands_kq,ebands_kmq
212 17 : type(hdr_type) :: hdr1_pq,hdr1_mq
213 17 : type(hdr_type) :: hdr2_pq,hdr2_mq
214 : !arrays
215 : integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/)
216 17 : integer,allocatable :: kg1(:,:),kg1_mq(:,:)
217 17 : integer,allocatable :: npwar1(:),npwar1_mq(:),npwtot1(:),npwtot1_mq(:)
218 : real(dp) :: d3etot_td(2),d3etot_td_mq(2)
219 17 : real(dp),allocatable :: cg1(:,:),cg2(:,:)
220 17 : real(dp),allocatable :: cg1_mq(:,:),cg2_mq(:,:)
221 17 : real(dp),allocatable :: d3etot_t4(:,:),d3etot_t5(:,:),d3etot_tgeom(:,:),eigen1(:),eigen2(:)
222 17 : real(dp),allocatable :: eigenq(:), doccde_tmp(:)
223 17 : real(dp),allocatable :: eigen1_mq(:),eigen2_mq(:)
224 17 : real(dp),allocatable :: kpq(:,:),kmq(:,:)
225 17 : real(dp),allocatable :: nhat1(:,:)
226 17 : real(dp),allocatable :: rho1g1(:,:),rho1r1(:,:)
227 17 : real(dp),allocatable :: rho2g1(:,:),rho2r1(:,:)
228 17 : real(dp),allocatable :: t4_typeI(:,:,:,:,:,:),t4_typeII(:,:,:,:,:,:,:)
229 17 : real(dp),allocatable :: t5_typeI(:,:,:,:,:,:),t5_typeII(:,:,:,:,:,:,:)
230 17 : real(dp),allocatable :: tgeom_typeI(:,:,:,:,:,:),tgeom_typeII(:,:,:,:,:,:,:)
231 17 : real(dp),allocatable :: vhart1dqdq(:),vpsp1dqdq(:)
232 17 : real(dp),allocatable :: vpsp1_i1pertdq(:,:,:),vpsp1_i2pertdq(:,:,:)
233 17 : real(dp),allocatable :: vxccc1_i2pertdq(:,:,:)
234 17 : real(dp),allocatable :: vpsp1_i1pertdq_geom(:,:,:), vpsp1_i1pertdqdq(:,:,:)
235 17 : real(dp),allocatable :: vxc1dqdq(:),work(:)
236 17 : real(dp),allocatable :: xccc3d2(:),xccc3d2dq(:)
237 17 : type(pawrhoij_type),allocatable :: pawrhoij_read(:)
238 :
239 :
240 : ! *************************************************************************
241 :
242 : DBG_ENTER("COLL")
243 :
244 : !Init parallelism
245 17 : comm_cell=mpi_enreg%comm_cell
246 17 : me=mpi_enreg%me_kpt
247 :
248 : !Deactivate kramer's degeneracy for calculations that break TRS
249 : !(finite q in magnetic materials or finite omega)
250 17 : kramers_deg=.true.
251 17 : timrev = 1
252 17 : if (dtset%tim1rev==0) then
253 2 : kramers_deg=.false.
254 2 : timrev = 0
255 : end if
256 :
257 : !Special treatment for timdisp calculations at finite q
258 17 : finite_q=.false.
259 68 : if (sum(dtset%qptn(:)**2)>tol8) finite_q=.true.
260 :
261 : !Various initializations
262 17 : cplex = 2 - timrev
263 17 : nspden = dtset%nspden
264 17 : ecut=dtset%ecut
265 17 : ecut_eff = ecut*(dtset%dilatmx)**2
266 17 : mpsang = psps%mpsang
267 17 : optorth=1;if (psps%usepaw==1) optorth=0
268 17 : opthartdqdq=1
269 :
270 : !Set up some quantities required to calculate the Berry
271 : !curvature at finite q
272 51 : ABI_MALLOC(npwar1,(dtset%nkpt))
273 34 : ABI_MALLOC(npwtot1,(dtset%nkpt))
274 51 : ABI_MALLOC(kpq,(3,dtset%nkpt))
275 17 : if (.not.kramers_deg) then
276 4 : ABI_MALLOC(npwar1_mq,(dtset%nkpt))
277 4 : ABI_MALLOC(npwtot1_mq,(dtset%nkpt))
278 4 : ABI_MALLOC(kmq,(3,dtset%nkpt))
279 : end if
280 17 : if (finite_q.or..not.kramers_deg) then
281 :
282 : !Initialize k+q (and k-q) array
283 26 : do ikpt=1,dtset%nkpt
284 98 : kpq(:,ikpt)=dtset%qptn(:)+dtset%kptns(:,ikpt)
285 : end do
286 2 : if (.not.kramers_deg) then
287 26 : do ikpt=1,nkpt
288 98 : kmq(:,ikpt)=-dtset%qptn(:)+dtset%kptns(:,ikpt)
289 : end do
290 : end if
291 :
292 : !Compute maximum number of planewaves at k+q (and k-q)
293 2 : call getmpw(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kpq,mpi_enreg,mpw1,dtset%nkpt)
294 2 : if (.not.kramers_deg) then
295 2 : call getmpw(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kmq,mpi_enreg,mpw1_mq,dtset%nkpt)
296 : !number of plane waves at k+q and k-q should be in principle the same to reconstruct rhor1_pq (?)
297 2 : mpw1=max(mpw1,mpw1_mq)
298 2 : mpw1_mq=mpw1
299 : else
300 0 : mpw1_mq=0
301 : end if
302 :
303 : !Set up the basis sphere of planewaves at k+q (and k-q)
304 6 : ABI_MALLOC(kg1,(3,mpw1*mk1mem))
305 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg1,&
306 : & kpq,mk1mem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,mpw1,&
307 2 : & npwar1,npwtot1,dtset%nsppol)
308 2 : if (.not.kramers_deg) then
309 6 : ABI_MALLOC(kg1_mq,(3,mpw1_mq*mk1mem))
310 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg1_mq,&
311 : & kmq,mk1mem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,mpw1_mq,&
312 2 : & npwar1_mq,npwtot1_mq,dtset%nsppol)
313 : end if
314 :
315 : ! Initialize band structure datatype at k+q
316 26 : bantot=sum(dtset%nband(1:dtset%nkpt*dtset%nsppol))
317 6 : ABI_MALLOC(eigenq,(bantot))
318 6 : ABI_MALLOC(doccde_tmp,(dtset%mband*dtset%nkpt*dtset%nsppol))
319 1010 : eigenq(:)=zero
320 : call ebands_kq%init(bantot,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
321 : & doccde_tmp,eigenq,dtset%istwfk,kpq,dtset%nband,dtset%nkpt,npwar1,dtset%nsppol,&
322 : & dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ,dtset%wtk,&
323 : & dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
324 2 : & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
325 2 : if (.not.kramers_deg) then
326 1010 : eigenq(:)=zero
327 : call ebands_kmq%init(bantot,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
328 : & doccde_tmp,eigenq,dtset%istwfk,kmq,dtset%nband,dtset%nkpt,npwar1_mq,dtset%nsppol,&
329 : & dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ,dtset%wtk,&
330 : & dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
331 2 : & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
332 : end if
333 2 : ABI_FREE(eigenq)
334 4 : ABI_FREE(doccde_tmp)
335 :
336 : else
337 45 : ABI_MALLOC(kg1,(3,mpw*mk1mem))
338 4158 : kpq=dtset%kptns(:,1:dtset%nkpt)
339 15 : mpw1=mpw
340 658542 : kg1=kg
341 1062 : npwar1=npwarr
342 : end if
343 17 : mcg1=mpw1*nspinor*mband*mk1mem*nsppol
344 :
345 51 : ABI_MALLOC(cg1,(2,mcg1))
346 34 : ABI_MALLOC(cg2,(2,mcg1))
347 51 : ABI_MALLOC(eigen1,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
348 34 : ABI_MALLOC(eigen2,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
349 17 : if (.not.kramers_deg) then
350 4 : ABI_MALLOC(cg1_mq,(2,mcg1))
351 4 : ABI_MALLOC(cg2_mq,(2,mcg1))
352 4 : ABI_MALLOC(eigen1_mq,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
353 4 : ABI_MALLOC(eigen2_mq,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
354 : end if
355 68 : ABI_MALLOC(rho1r1,(cplex*nfftf,dtset%nspden))
356 51 : ABI_MALLOC(rho2r1,(cplex*nfftf,dtset%nspden))
357 51 : ABI_MALLOC(rho1g1,(2,nfftf))
358 34 : ABI_MALLOC(rho2g1,(2,nfftf))
359 :
360 17 : ask_accurate=1 ; formeig = 1 ; ireadwf = 1
361 17 : n1=ngfftf(1) ; n2=ngfftf(2) ; n3=ngfftf(3)
362 17 : nfftotf=n1*n2*n3
363 :
364 : !Allocations for type-I terms
365 119 : ABI_MALLOC(t4_typeII,(2,3,mpert,3,mpert,3,mpert))
366 6485593 : t4_typeII(:,:,:,:,:,:,:)=zero
367 17 : if (d3e_pert2(natom+3)==1.or.d3e_pert2(natom+4)==1) then
368 15 : ABI_MALLOC(t4_typeI,(2,3,mpert,3,3,3))
369 8234 : t4_typeI(:,:,:,:,:,:)=zero
370 : end if
371 68 : ABI_MALLOC(t5_typeII,(2,3,mpert,3,mpert,3,mpert))
372 6485593 : t5_typeII(:,:,:,:,:,:,:)=zero
373 17 : if (d3e_pert1(natom+3)==1.or.d3e_pert1(natom+4)==1) then
374 0 : ABI_MALLOC(t5_typeI,(2,3,mpert,3,3,3))
375 17 : t5_typeI(:,:,:,:,:,:)=zero
376 : end if
377 68 : ABI_MALLOC(tgeom_typeII,(2,3,mpert,3,mpert,3,mpert))
378 6485593 : tgeom_typeII(:,:,:,:,:,:,:)=zero
379 55 : if (any(d3e_pert1(1:natom)==1).and.(d3e_pert2(natom+3)==1.or.d3e_pert2(natom+4)==1)) then
380 15 : ABI_MALLOC(tgeom_typeI,(2,3,mpert,3,3,3))
381 8234 : tgeom_typeI(:,:,:,:,:,:)=zero
382 : end if
383 :
384 :
385 : !==== Initialize most of the Hamiltonian (and derivative) ====
386 : !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
387 : !2) Perform the setup needed for the non-local factors:
388 : !3) Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
389 : call gs_hamkq%init(psps,pawtab,dtset%nspinor,dtset%nsppol,nspden,dtset%natom,&
390 : & dtset%typat,xred,dtset%nfft,mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,&
391 17 : & gpu_option=dtset%gpu_option)
392 :
393 : !Specific allocations for strain-gradient perturbation
394 17 : if (dtset%lw_flexo==1.or.dtset%lw_flexo==2.or.dtset%lw_flexo==4) then
395 9 : ABI_MALLOC(vhart1dqdq,(2*nfftf))
396 6 : ABI_MALLOC(vpsp1dqdq,(2*nfftf))
397 6 : ABI_MALLOC(vxc1dqdq,(2*nfftf))
398 : end if
399 :
400 : !Allocate xccc (only for quadrupoles yet)
401 17 : if (dtset%lw_qdrpl==1.and.psps%n1xccc/=0) then
402 6 : ABI_MALLOC(xccc3d2,(cplex*nfftf))
403 : end if
404 :
405 : !This is necessary to deactivate paw options in the dfpt_rhotov routine
406 17 : ABI_MALLOC(pawrhoij_read,(0))
407 17 : usexcnhat=0
408 17 : n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf
409 17 : pawread=0
410 17 : nhat1grdim=0
411 68 : ABI_MALLOC(nhat1,(cplex*dtset%nfft,nspden))
412 223617 : nhat1=zero
413 :
414 :
415 17 : pert1case = 0 ; pert2case = 0 ; pert3case = 0
416 17 : pert1case_mq = 0 ; pert2case_mq = 0
417 :
418 259 : do i1pert = 1, mpert
419 985 : do i1dir = 1, 3
420 :
421 2623064 : if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then
422 :
423 55 : if (i1pert < natom + 10) then
424 52 : pert1case = i1dir + (i1pert-1)*3
425 3 : else if (i1pert > natom+11 .and. i1pert <= 2*natom+11) then
426 3 : pert1case = i1dir + i1pert*3
427 : end if
428 55 : call appdig(pert1case,dtfil%fnamewff1,fiwf1i)
429 :
430 55 : if (finite_q) then
431 : call hdr1_pq%init(ebands_kq,codvsn,dtset,pawtab,pert1case,psps,wvl%descr, &
432 4 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )
433 : else
434 51 : hdr1_pq=hdr
435 : end if
436 :
437 : call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
438 : & formeig,hdr1_pq,ireadwf,dtset%istwfk,kg1,kpq,dtset%localrdwf,&
439 : & dtset%mband,mcg1,dtset%mk1mem,mpi_enreg,mpw1,&
440 : & dtset%nband,dtset%ngfft,dtset%nkpt,npwar1,&
441 : & dtset%nsppol,dtset%nsym,&
442 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
443 55 : & dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl)
444 :
445 55 : if (ireadwf==1) then
446 55 : call WffClose (wff1,ierr)
447 : end if
448 :
449 55 : if (.not.kramers_deg) then
450 4 : if (finite_q) then
451 4 : pert1case_mq=pert1case+(2*dtset%natom+11)*3
452 : else
453 0 : pert1case_mq=pert1case
454 : endif
455 4 : call appdig(pert1case_mq,dtfil%fnamewff1,fiwf1i_mq)
456 :
457 4 : if (finite_q) then
458 : call hdr1_mq%init(ebands_kmq,codvsn,dtset,pawtab,pert1case_mq,psps,wvl%descr, &
459 4 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )
460 : else
461 0 : hdr1_mq=hdr
462 : end if
463 :
464 : call inwffil(ask_accurate,cg1_mq,dtset,dtset%ecut,ecut_eff,eigen1_mq,dtset%exchn2n3d,&
465 : & formeig,hdr1_mq,ireadwf,dtset%istwfk,kg1_mq,kmq,dtset%localrdwf,&
466 : & dtset%mband,mcg1,dtset%mk1mem,mpi_enreg,mpw1,&
467 : & dtset%nband,dtset%ngfft,dtset%nkpt,npwar1_mq,&
468 : & dtset%nsppol,dtset%nsym,&
469 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
470 4 : & dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i_mq,wvl)
471 :
472 4 : if (ireadwf==1) then
473 4 : call WffClose (wff1,ierr)
474 : end if
475 : end if
476 :
477 55 : if (.not.just_timdisp) then
478 43 : call read_1eig(eigen1,formeig,mband,nkpt,nsppol,fiwf1i)
479 :
480 1377834 : rho1r1(:,:) = zero; rho1g1(:,:) = zero
481 43 : if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
482 43 : call appdig(pert1case,dtfil%fildens1in,fiden1i)
483 :
484 : call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, psps%usepaw, mpi_enreg, rho1r1, &
485 43 : hdr_den, pawrhoij_read, comm_cell, check_hdr=hdr)
486 43 : call hdr_den%free()
487 : end if
488 :
489 : !Perform FFT rhor1 to rhog1
490 129 : ABI_MALLOC(work,(cplex*nfftf))
491 344480 : work(:)=rho1r1(:,1)
492 43 : call fourdp(cplex,rho1g1,work,-1,mpi_enreg,dtset%nfft,1,dtset%ngfft,0)
493 43 : ABI_FREE(work)
494 : end if
495 :
496 : !Allocate the first-order gradient local potential
497 55 : if (i1pert <= natom+3) then
498 49 : n1dq=1
499 196 : ABI_MALLOC(vpsp1_i1pertdq,(2*nfftf,dtset%nspden,n1dq))
500 6 : else if (i1pert == natom+4) then
501 0 : n1dq=2
502 0 : ABI_MALLOC(vpsp1_i1pertdq,(2*nfftf,dtset%nspden,n1dq))
503 : else
504 6 : n1dq=1
505 : end if
506 165 : ABI_MALLOC(d3etot_t5,(2,n1dq))
507 :
508 737 : do i2pert = 1, mpert
509 2783 : do i2dir = 1, 3
510 :
511 117016 : if ((maxval(rfpert(i1dir,i1pert,i2dir,i2pert,:,:))==1)) then
512 :
513 242 : if (i2pert < natom + 10) then
514 236 : pert2case = i2dir + (i2pert-1)*3
515 6 : else if (i2pert > natom+11 .and. i2pert <= 2*natom+11) then
516 6 : pert2case = i2dir + i2pert*3
517 : end if
518 :
519 242 : call appdig(pert2case,dtfil%fnamewff1,fiwf2i)
520 :
521 242 : if (finite_q) then
522 : call hdr2_pq%init(ebands_kq,codvsn,dtset,pawtab,pert2case,psps,wvl%descr, &
523 6 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )
524 : else
525 236 : hdr2_pq=hdr
526 : end if
527 :
528 : call inwffil(ask_accurate,cg2,dtset,dtset%ecut,ecut_eff,eigen2,dtset%exchn2n3d,&
529 : & formeig,hdr2_pq,ireadwf,dtset%istwfk,kg1,kpq,dtset%localrdwf,&
530 : & dtset%mband,mcg1,dtset%mk1mem,mpi_enreg,mpw1,&
531 : & dtset%nband,dtset%ngfft,dtset%nkpt,npwar1,&
532 : & dtset%nsppol,dtset%nsym,&
533 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
534 242 : & dtfil%unkg1,wff2,wfft2,dtfil%unwff2,fiwf2i,wvl)
535 :
536 242 : if (ireadwf==1) then
537 242 : call WffClose (wff2,ierr)
538 : end if
539 :
540 242 : if (.not.kramers_deg) then
541 6 : if (finite_q) then
542 6 : pert2case_mq=pert2case+(2*dtset%natom+11)*3
543 : else
544 0 : pert2case_mq=pert2case
545 : end if
546 6 : call appdig(pert2case_mq,dtfil%fnamewff1,fiwf2i_mq)
547 :
548 6 : if (finite_q) then
549 : call hdr2_mq%init(ebands_kmq,codvsn,dtset,pawtab,pert2case_mq,psps,wvl%descr, &
550 6 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )
551 : else
552 0 : hdr2_mq=hdr
553 : end if
554 :
555 : call inwffil(ask_accurate,cg2_mq,dtset,dtset%ecut,ecut_eff,eigen2_mq,dtset%exchn2n3d,&
556 : & formeig,hdr2_mq,ireadwf,dtset%istwfk,kg1_mq,kmq,dtset%localrdwf,&
557 : & dtset%mband,mcg1,dtset%mk1mem,mpi_enreg,mpw1,&
558 : & dtset%nband,dtset%ngfft,dtset%nkpt,npwar1_mq,&
559 : & dtset%nsppol,dtset%nsym,&
560 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
561 6 : & dtfil%unkg1,wff2,wfft2,dtfil%unwff2,fiwf2i_mq,wvl)
562 :
563 6 : if (ireadwf==1) then
564 6 : call WffClose (wff2,ierr)
565 : end if
566 : end if
567 :
568 242 : if (i1pert==i2pert.and.i1dir==i2dir) then
569 14 : samepert=.true.
570 : else
571 228 : samepert=.false.
572 : end if
573 :
574 242 : if (.not.just_timdisp) then
575 222 : call read_1eig(eigen2,formeig,mband,nkpt,nsppol,fiwf2i)
576 :
577 5990380 : rho2r1(:,:) = zero; rho2g1(:,:) = zero
578 222 : if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
579 222 : call appdig(pert2case,dtfil%fildens1in,fiden1i)
580 :
581 : call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, psps%usepaw, mpi_enreg, rho2r1, &
582 222 : hdr_den, pawrhoij_read, comm_cell, check_hdr=hdr)
583 222 : call hdr_den%free()
584 : end if
585 :
586 222 : if (.not.samepert) then
587 : !Perform FFT rhor1 to rhog1
588 654 : ABI_MALLOC(work,(cplex*nfftf))
589 1473753 : work(:)=rho2r1(:,1)
590 218 : call fourdp(cplex,rho2g1,work,-1,mpi_enreg,dtset%nfft,1,dtset%ngfft,0)
591 218 : ABI_FREE(work)
592 : end if !samepert
593 : end if
594 :
595 : !Allocate the first-order gradient local potential
596 242 : if (i2pert <= natom+3) then
597 149 : n2dq=1
598 596 : ABI_MALLOC(vpsp1_i2pertdq,(2*nfftf,dtset%nspden,n2dq))
599 596 : ABI_MALLOC(vxccc1_i2pertdq,(2*nfftf,dtset%nspden,n2dq))
600 447 : ABI_MALLOC(xccc3d2dq,(2*nfftf))
601 93 : else if (i2pert == natom+4) then
602 81 : n2dq=2
603 405 : ABI_MALLOC(vpsp1_i2pertdq,(2*nfftf,dtset%nspden,n2dq))
604 : else
605 12 : n2dq=1
606 : end if
607 726 : ABI_MALLOC(d3etot_t4,(2,n2dq))
608 484 : ABI_MALLOC(d3etot_tgeom,(2,n2dq))
609 :
610 : !Calculate the first-order gradient local potential that enters the geometric term
611 1210 : ABI_MALLOC(vpsp1_i1pertdq_geom,(2*nfftf,dtset%nspden,3))
612 242 : if (.not.just_timdisp) then
613 222 : if (i1pert <= natom .and. (i2pert == natom+3.or.i2pert == natom+4)) then
614 :
615 : !calculate the second of the two first-gradient directions
616 432 : do ii=1,3
617 : call dfpt_vlocaldq(atindx,2,gmet,gsqcut,i1dir,i1pert,mpi_enreg, &
618 : & psps%mqgrid_vl,dtset%natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3, &
619 546 : & ph1d,ii,psps%qgrid_vl,dtset%qptn,ucvol,psps%vlspl,vpsp1_i1pertdq_geom(:,1,ii))
620 : end do
621 : end if
622 : end if
623 :
624 : !Allocate the second-gradient array
625 1210 : ABI_MALLOC(vpsp1_i1pertdqdq,(2*nfftf,dtset%nspden,n2dq))
626 :
627 : !Calculate first-order pseudocore charge (still, only for quadrupoles)
628 242 : if (.not.just_timdisp) then
629 222 : if (i1pert==natom+2.and.i2pert<=natom.and.psps%n1xccc/=0.and.nkxc == 7) then
630 9 : if (psps%nc_xccc_gspace==1) then
631 0 : ndir=1
632 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i2dir,i2pert,&
633 : & mgfft,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,&
634 : & ph1d,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
635 0 : & atmrhor1=xccc3d2,optn_in=n3xccc/nfftf,optn2_in=1)
636 9 : else if (psps%nc_xccc_gspace==0) then
637 : call dfpt_mkcore(cplex,i2dir,i2pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,&
638 9 : & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d2,xred)
639 : end if
640 : end if
641 : end if
642 :
643 2898 : do i3pert = 1, mpert
644 10866 : do i3dir = 1, 3
645 :
646 10624 : if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
647 :
648 611 : blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
649 :
650 : !Calculate local potentials for first-order gradient Hamiltonians
651 : !gradient of i1pert:
652 611 : if (.not.just_timdisp) then
653 591 : if (i1pert<=natom) then
654 : !Get q-gradient of first-order local part of the pseudopotential
655 : call dfpt_vlocaldq(atindx,2,gmet,gsqcut,i1dir,i1pert,mpi_enreg, &
656 : & psps%mqgrid_vl,dtset%natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3, &
657 354 : & ph1d,i3dir,psps%qgrid_vl,dtset%qptn,ucvol,psps%vlspl,vpsp1_i1pertdq(:,1,1))
658 :
659 354 : if (i2pert == natom+3.or.i2pert == natom+4) then
660 324 : gamma=i3dir
661 810 : do idq= 1, n2dq
662 486 : if (i2pert==natom+3) then
663 162 : istr=i2dir
664 : else
665 324 : istr=idq*3+i2dir
666 : endif
667 486 : delta=idx(2*istr)
668 : call dfpt_vlocaldqdq(atindx,2,gs_hamkq%gmet,gsqcut,i1dir,i1pert,mpi_enreg, &
669 : & psps%mqgrid_vl,dtset%natom,&
670 : & nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3, &
671 : & ph1d,gamma,delta,psps%qgrid_vl,&
672 840 : & dtset%qptn,ucvol,psps%vlspl,vpsp1_i1pertdqdq(:,1,idq))
673 : end do
674 : end if
675 :
676 237 : else if (i1pert==natom+3.or.i1pert==natom+4) then
677 0 : istr=i1dir; if (i1pert==natom+4) istr=3+i1dir
678 : !Get 2nd q-gradient of first-order local part of the pseudopotential and of the Hartree
679 : !(and XC if GGA) contribution from ground state density
680 : call dfpt_vmetdqdq(2,gmet,gprimd,gsqcut,istr,i1pert,kxc,mpi_enreg, &
681 : & psps%mqgrid_vl,natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3,&
682 : & nkxc,nspden,opthartdqdq,ph1d,i3dir,psps%qgrid_vl,&
683 0 : & dtset%qptn,rhog,rhor,ucvol,psps%vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq)
684 0 : vpsp1_i1pertdq(:,1,1)=vhart1dqdq(:)+vpsp1dqdq(:)+vxc1dqdq(:)
685 0 : if (i1pert==natom+4) then
686 : !Here we need to calculate both extradiagonal shear-strains
687 : !because the second gradient of the metric perturbation is
688 : !type-I, i.e., non symmetric with respect to the
689 : !permutation of the strain indexes.
690 0 : istr=6+i1dir
691 : call dfpt_vmetdqdq(2,gmet,gprimd,gsqcut,istr,i1pert,kxc,mpi_enreg, &
692 : & psps%mqgrid_vl,natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3,&
693 : & nkxc,nspden,opthartdqdq,ph1d,i3dir,psps%qgrid_vl,&
694 0 : & dtset%qptn,rhog,rhor,ucvol,psps%vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq)
695 0 : vpsp1_i1pertdq(:,1,2)=vhart1dqdq(:)+vpsp1dqdq(:)+vxc1dqdq(:)
696 : end if
697 : end if
698 :
699 591 : if (.not.samepert) then
700 : !gradient of i2pert:
701 584 : if (i2pert<=natom) then
702 : !Get q-gradient of first-order local part of the pseudopotential
703 : call dfpt_vlocaldq(atindx,2,gmet,gsqcut,i2dir,i2pert,mpi_enreg, &
704 : & psps%mqgrid_vl,dtset%natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3, &
705 94 : & ph1d,i3dir,psps%qgrid_vl,dtset%qptn,ucvol,psps%vlspl,vpsp1_i2pertdq(:,1,1))
706 490 : else if (i2pert==natom+3.or.i2pert==natom+4) then
707 486 : istr=i2dir; if (i2pert==natom+4) istr=3+i2dir
708 : !Get 2nd q-gradient of first-order local part of the pseudopotential and of the Hartree
709 : !(and XC if GGA) contribution from ground state density
710 : call dfpt_vmetdqdq(2,gmet,gprimd,gsqcut,istr,i2pert,kxc,mpi_enreg, &
711 : & psps%mqgrid_vl,natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3,&
712 : & nkxc,nspden,opthartdqdq,ph1d,i3dir,psps%qgrid_vl,&
713 486 : & dtset%qptn,rhog,rhor,ucvol,psps%vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq)
714 3280986 : vpsp1_i2pertdq(:,1,1)=vhart1dqdq(:)+vpsp1dqdq(:)+vxc1dqdq(:)
715 486 : if (i2pert==natom+4) then
716 : !Here we need to calculate both extradiagonal shear-strains
717 : !because the second gradient of the metric perturbation is
718 : !type-I, i.e., non symmetric with respect to the
719 : !permutation of the strain indexes.
720 243 : istr=6+i2dir
721 : call dfpt_vmetdqdq(2,gmet,gprimd,gsqcut,istr,i2pert,kxc,mpi_enreg, &
722 : & psps%mqgrid_vl,natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3,&
723 : & nkxc,nspden,opthartdqdq,ph1d,i3dir,psps%qgrid_vl,&
724 243 : & dtset%qptn,rhog,rhor,ucvol,psps%vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq)
725 1640493 : vpsp1_i2pertdq(:,1,2)=vhart1dqdq(:)+vpsp1dqdq(:)+vxc1dqdq(:)
726 : end if
727 : end if
728 :
729 584 : if (i1pert==natom+2.and.i2pert<=natom.and.psps%n1xccc/=0) then
730 : !Get the q-gradient of the pseudocore density
731 : call dfpt_vlocaldq(atindx,2,gmet,gsqcut,i2dir,i2pert,mpi_enreg, &
732 : & psps%mqgrid_vl,dtset%natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3, &
733 27 : & ph1d,i3dir,psps%qgrid_vl,dtset%qptn,ucvol,ncorespl,xccc3d2dq,optnc=1)
734 :
735 : !Get the q-gradient of the first-order XC potential due to the pseudocore charge
736 : call dfpt_mkvxcccdq(cplex,i3dir,dtset%ixc,gprimd,kxc,mpi_enreg,nfftf,dtset%ngfft,&
737 27 : & nkxc,nspden,dtset%qptn,rprimd,vxccc1_i2pertdq,xccc3d2,xccc3d2dq)
738 :
739 : !Add this contribution to the gradient of the local PSP
740 244716 : vpsp1_i2pertdq= vpsp1_i2pertdq + vxccc1_i2pertdq
741 : end if
742 : end if !samepert
743 :
744 :
745 : !Prepare ddk wf file
746 591 : pert3case = i3dir + natom*3
747 591 : call appdig(pert3case,dtfil%fnamewffddk,fiwfddk)
748 : ! Checking the existence of data file
749 591 : if (.not. file_exists(fiwfddk)) then
750 : ! Trick needed to run Abinit test suite in netcdf mode.
751 591 : if (file_exists(nctk_ncify(fiwfddk))) then
752 591 : write(message,"(3a)")"- File: ",trim(fiwfddk),&
753 1182 : " does not exist but found netcdf file with similar name."
754 591 : call wrtout(std_out,message,'COLL')
755 1182 : fiwfddk = nctk_ncify(fiwfddk)
756 : end if
757 591 : if (.not. file_exists(fiwfddk)) then
758 0 : ABI_ERROR('Missing file: '//TRIM(fiwfddk))
759 : end if
760 : end if
761 591 : write(message,'(2a)')'-dfptlw_loop : read the ddk wavefunctions from file: ',trim(fiwfddk)
762 591 : call wrtout(std_out,message,'COLL')
763 : !call wrtout(ab_out,message,'COLL')
764 : !Note that the unit number for these files is 50,51,52 or 53 (dtfil%unddk=50)
765 591 : call ddk_f%open_read(fiwfddk,1,dtset%iomode,dtfil%unddk,mpi_enreg%comm_cell)
766 :
767 : !Prepare d2_dkdk wf file
768 : !For i1pert
769 591 : if (i1pert==natom+2) then
770 237 : call rf2_getidir(i1dir,i3dir,idir_dkdk)
771 : !if (idir_dkdk>6) idir_dkdk=idir_dkdk-3
772 237 : dkdk_index=idir_dkdk+(dtset%natom+6)*3
773 237 : call appdig(dkdk_index,dtfil%fnamewffdkdk,fiwfdkdk)
774 : !Check that d2_dkdk file exists and open it
775 237 : if (.not. file_exists(fiwfdkdk)) then
776 : ! Trick needed to run Abinit test suite in netcdf mode.
777 237 : if (file_exists(nctk_ncify(fiwfdkdk))) then
778 237 : write(message,"(3a)")"- File: ",trim(fiwfdkdk),&
779 474 : " does not exist but found netcdf file with similar name."
780 237 : call wrtout(std_out,message,'COLL')
781 474 : fiwfdkdk = nctk_ncify(fiwfdkdk)
782 : end if
783 237 : if (.not. file_exists(fiwfdkdk)) then
784 0 : ABI_ERROR('Missing file: '//TRIM(fiwfdkdk))
785 : end if
786 : end if
787 237 : write(message,'(2a)')'-dfptlw_loop : read the d2_dkdk wavefunctions from file: ',trim(fiwfdkdk)
788 237 : call wrtout(std_out,message,'COLL')
789 : !call wrtout(ab_out,message,'COLL')
790 237 : call d2_dkdk_f%open_read(fiwfdkdk,1,dtset%iomode,dtfil%unddk+1,mpi_enreg%comm_cell)
791 :
792 : end if
793 :
794 : !Prepare d2_dkdk wf file
795 : !For i2pert
796 591 : if (i2pert==natom+2.and..not.samepert) then
797 4 : call rf2_getidir(i2dir,i3dir,idir_dkdk)
798 : !if (idir_dkdk>6) idir_dkdk=idir_dkdk-3
799 4 : dkdk_index=idir_dkdk+(dtset%natom+6)*3
800 4 : call appdig(dkdk_index,dtfil%fnamewffdkdk,fiwfdkdk)
801 : !Check that d2_dkdk file exists and open it
802 4 : if (.not. file_exists(fiwfdkdk)) then
803 : ! Trick needed to run Abinit test suite in netcdf mode.
804 4 : if (file_exists(nctk_ncify(fiwfdkdk))) then
805 4 : write(message,"(3a)")"- File: ",trim(fiwfdkdk),&
806 8 : " does not exist but found netcdf file with similar name."
807 4 : call wrtout(std_out,message,'COLL')
808 8 : fiwfdkdk = nctk_ncify(fiwfdkdk)
809 : end if
810 4 : if (.not. file_exists(fiwfdkdk)) then
811 0 : ABI_ERROR('Missing file: '//TRIM(fiwfdkdk))
812 : end if
813 : end if
814 4 : write(message,'(2a)')'-dfptlw_loop : read the d2_dkdk wavefunctions from file: ',trim(fiwfdkdk)
815 4 : call wrtout(std_out,message,'COLL')
816 4 : call d2_dkdk_f2%open_read(fiwfdkdk,1,dtset%iomode,dtfil%unddk+2,mpi_enreg%comm_cell)
817 : end if
818 :
819 : !Perform the longwave DFPT part of the 3dte calculation
820 : call dfptlw_pert(cg,cg1,cg2,cplex,d3etot,d3etot_t4,d3etot_t5,d3etot_tgeom,dimffnl,dtset, &
821 : & eigen1,eigen2,ffnl,gmet,gs_hamkq,gsqcut,i1dir,&
822 : & i2dir,i3dir,i1pert,i2pert,i3pert,kg,kxc,mband,mkmem,mk1mem,mpert,mpi_enreg,&
823 : & mpsang,mpw1,natom,n1dq,n2dq,nfftf,ngfftf,nkpt,nkxc,nspden,nspinor,nsppol,npwarr,nylmgr,occ,&
824 : & pawfgr,psps,rho1g1,rho1r1,rho2r1,rmet,rprimd,samepert,ucvol,useylmgr,&
825 : & vpsp1_i1pertdq,vpsp1_i1pertdqdq,vpsp1_i1pertdq_geom,vpsp1_i2pertdq,&
826 591 : & ddk_f,d2_dkdk_f,d2_dkdk_f2,ylm,ylmgr)
827 :
828 : !close ddk file
829 591 : call ddk_f%close()
830 :
831 : !close d2_dkdk file (i1pert)
832 591 : if (i1pert==natom+2) call d2_dkdk_f%close()
833 :
834 : ! Close d2_dkdk file (i2pert)
835 591 : if (i2pert==natom+2.and..not.samepert) call d2_dkdk_f2%close()
836 :
837 : !Save the type-I terms
838 591 : if (i2pert==natom+3.or.i2pert==natom+4) then
839 486 : gamma=i3dir
840 1215 : do idq=1,n2dq
841 729 : if (i2pert==natom+3) then
842 243 : istr=i2dir
843 : else
844 486 : istr=idq*3+i2dir
845 : endif
846 729 : beta=idx(2*istr-1); delta=idx(2*istr)
847 2673 : t4_typeI(:,i1dir,i1pert,beta,delta,gamma)=d3etot_t4(:,idq)
848 : end do
849 : else
850 315 : t4_typeII(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=d3etot_t4(:,1)
851 : end if
852 :
853 591 : if (i1pert==natom+3.or.i1pert==natom+4) then
854 0 : gamma=i3dir
855 0 : do idq=1,n1dq
856 0 : if (i1pert==natom+3) then
857 0 : istr=i1dir
858 : else
859 0 : istr=idq*3+i1dir
860 : endif
861 0 : beta=idx(2*istr-1); delta=idx(2*istr)
862 0 : t5_typeI(:,i2dir,i2pert,beta,delta,gamma)=d3etot_t5(:,idq)
863 : end do
864 : else
865 1773 : t5_typeII(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=d3etot_t5(:,1)
866 : end if
867 :
868 :
869 591 : if (i1pert<=natom.and.(i2pert==natom+3.or.i2pert==natom+4)) then
870 324 : alpha=i1dir
871 324 : gamma=i3dir
872 810 : do idq=1,n2dq
873 486 : if (i2pert==natom+3) then
874 162 : istr=i2dir
875 : else
876 324 : istr=idq*3+i2dir
877 : endif
878 486 : beta=idx(2*istr-1); delta=idx(2*istr)
879 1458 : tgeom_typeI(:,i1dir,i1pert,beta,delta,gamma)=d3etot_tgeom(:,idq)
880 :
881 : !Incorporate here the G=0 contribution of the geometric term
882 486 : ia1=0
883 486 : itypat=0
884 972 : do ii=1,dtset%ntypat
885 486 : ia1=ia1+nattyp(ii)
886 972 : if (atindx(i1pert)<=ia1.and.itypat==0) itypat=ii
887 : end do
888 486 : delad=zero ; if (alpha==delta) delad=one
889 486 : delbd=zero ; if (beta==delta) delbd=one
890 486 : delag=zero ; if (alpha==gamma) delag=one
891 486 : delbg=zero ; if (beta==gamma) delbg=one
892 :
893 : tgeom_typeI(1,i1dir,i1pert,beta,delta,gamma)= &
894 : & tgeom_typeI(1,i1dir,i1pert,beta,delta,gamma) + &
895 810 : & pi*pi*rhog(1,1)*psps%vlspl(1,2,itypat)*(delag*delbd+delad*delbg)
896 : end do
897 : else
898 801 : tgeom_typeII(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=d3etot_tgeom(:,1)
899 : end if
900 : end if ! .not.just_timdisp
901 :
902 611 : if (i3pert==natom+9) then
903 :
904 20 : write(message,'(2a,3(a,i2,a,i1))') ch10,'TIMDISP : ',&
905 40 : ' perts : ',i1pert,'.',i1dir,' / ',i2pert,'.',i2dir,' / ',i3pert,'.',i3dir
906 20 : call wrtout(std_out,message,'COLL')
907 20 : call wrtout(ab_out,message,'COLL')
908 : !Perform the Berry curvature part of the time-disperion 3dte calculation
909 : call dfpttd_berrycurv(cg1,cg2,d3etot_td,dtset,&
910 : & mband,mk1mem,mpi_enreg,&
911 20 : & mpw1,nkpt,nspinor,nsppol,npwar1,occ)
912 20 : if (.not.kramers_deg) then
913 : call dfpttd_berrycurv(cg1_mq,cg2_mq,d3etot_td_mq,dtset,&
914 : & mband,mk1mem,mpi_enreg,&
915 6 : & mpw1,nkpt,nspinor,nsppol,npwar1_mq,occ)
916 : end if
917 :
918 : !Add the result to the big array
919 : if (kramers_deg) then
920 14 : d3etot(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=zero
921 14 : d3etot(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=-two*d3etot_td(2)
922 : else
923 6 : d3etot(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=-d3etot_td(1)+d3etot_td_mq(1)
924 6 : d3etot(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=-d3etot_td(2)-d3etot_td_mq(2)
925 : end if
926 :
927 : end if
928 :
929 : end if ! rfpert
930 : end do ! ir3dir
931 : end do ! ir3pert
932 :
933 242 : if (finite_q) then
934 6 : call hdr1_pq%free()
935 6 : call hdr2_pq%free()
936 6 : if (.not.kramers_deg) then
937 6 : call hdr1_mq%free()
938 6 : call hdr2_mq%free()
939 : end if
940 : end if
941 :
942 242 : ABI_SFREE(vpsp1_i2pertdq)
943 242 : ABI_SFREE(xccc3d2dq)
944 242 : ABI_SFREE(vxccc1_i2pertdq)
945 242 : ABI_FREE(vpsp1_i1pertdq_geom)
946 242 : ABI_FREE(vpsp1_i1pertdqdq)
947 242 : ABI_FREE(d3etot_t4)
948 242 : ABI_FREE(d3etot_tgeom)
949 :
950 : end if ! rfpert
951 : end do ! i2dir
952 : end do ! i2pert
953 :
954 55 : ABI_SFREE(vpsp1_i1pertdq)
955 55 : ABI_FREE(d3etot_t5)
956 :
957 : end if ! rfpert
958 : end do ! i1dir
959 : end do ! i1pert
960 :
961 : !More memory cleaning
962 17 : call gs_hamkq%free()
963 17 : ABI_FREE(kpq)
964 17 : ABI_FREE(kg1)
965 17 : ABI_FREE(npwar1)
966 17 : ABI_FREE(npwtot1)
967 17 : ABI_FREE(cg1)
968 17 : ABI_FREE(cg2)
969 17 : ABI_FREE(eigen1)
970 17 : ABI_FREE(eigen2)
971 17 : ABI_FREE(rho1r1)
972 17 : ABI_FREE(rho2r1)
973 17 : ABI_FREE(rho1g1)
974 17 : ABI_FREE(rho2g1)
975 17 : ABI_FREE(nhat1)
976 17 : ABI_FREE(pawrhoij_read)
977 17 : ABI_SFREE(xccc3d2)
978 :
979 17 : if (dtset%lw_flexo==1.or.dtset%lw_flexo==2.or.dtset%lw_flexo==4) then
980 3 : ABI_FREE(vhart1dqdq)
981 3 : ABI_FREE(vpsp1dqdq)
982 3 : ABI_FREE(vxc1dqdq)
983 : end if
984 :
985 17 : call ebands_kq%free()
986 17 : if (.not.kramers_deg) then
987 2 : call ebands_kmq%free()
988 2 : ABI_FREE(cg1_mq)
989 2 : ABI_FREE(cg2_mq)
990 2 : ABI_FREE(eigen1_mq)
991 2 : ABI_FREE(eigen2_mq)
992 2 : ABI_FREE(kg1_mq)
993 2 : ABI_FREE(npwar1_mq)
994 2 : ABI_FREE(npwtot1_mq)
995 2 : ABI_FREE(kmq)
996 : end if
997 : !Treatment of T4 and T5 terms that have a q-gradient of a rf Hamiltonian
998 : !they need to be converted to type-II for strain perturbation
999 17 : if (d3e_pert2(natom+3)==1.or.d3e_pert2(natom+4)==1) then
1000 3 : optgeom=0
1001 3 : call dfptlw_typeIproc(blkflg,gprimd,optgeom,mpert,natom,rfpert,rprimd,t4_typeI,t4_typeII)
1002 : end if
1003 :
1004 17 : if (d3e_pert1(natom+3)==1.or.d3e_pert1(natom+4)==1) then
1005 0 : optgeom=0
1006 0 : call dfptlw_typeIproc(blkflg,gprimd,optgeom,mpert,natom,rfpert,rprimd,t5_typeI,t5_typeII)
1007 : end if
1008 :
1009 : !Tgeom has to be converted to type-II
1010 : !To do it we need to convert the three involve indexes to cartessian. Then,
1011 : !after type-II conversion the q-gradient index is back converted to reduced.
1012 55 : if (any(d3e_pert1(1:natom)==1).and.(d3e_pert2(natom+3)==1.or.d3e_pert2(natom+4)==1)) then
1013 3 : optgeom=1
1014 3 : call dfptlw_typeIproc(blkflg,gprimd,optgeom,mpert,natom,rfpert,rprimd,tgeom_typeI,tgeom_typeII)
1015 : end if
1016 :
1017 : !Incorporate T4, T5 and Tgeom to d3etot
1018 : d3etot(:,:,:,:,:,:,:)= d3etot(:,:,:,:,:,:,:) + &
1019 : & t4_typeII(:,:,:,:,:,:,:) + &
1020 : & t5_typeII(:,:,:,:,:,:,:) + &
1021 6485593 : & tgeom_typeII(:,:,:,:,:,:,:)
1022 :
1023 : !Anounce end of spatial-dispersion calculation
1024 17 : write(message, '(a,a,a,a)' ) ch10,ch10,&
1025 34 : & ' -- Spatial-dispersion 3rd-order derivatives completed -- ',ch10
1026 17 : call wrtout(std_out,message,'COLL')
1027 17 : call wrtout(ab_out,message,'COLL')
1028 :
1029 : !Deallocations
1030 17 : ABI_FREE(t4_typeII)
1031 17 : ABI_FREE(t5_typeII)
1032 17 : ABI_FREE(tgeom_typeII)
1033 17 : ABI_SFREE(t4_typeI)
1034 17 : ABI_SFREE(t5_typeI)
1035 17 : ABI_SFREE(tgeom_typeI)
1036 :
1037 : DBG_EXIT("COLL")
1038 :
1039 17 : end subroutine dfptlw_loop
1040 : !!***
1041 :
1042 : !!****f* ABINIT/m_dfptlw_loop/dfptlw_typeIproc
1043 : !! NAME
1044 : !! dfptlw_typeIproc
1045 : !!
1046 : !! FUNCTION
1047 : !! Process type-I terms and convert them to type-II in
1048 : !! the d3etot mixed (reduced/cartessian) coordinates.
1049 : !!
1050 : !! COPYRIGHT
1051 : !! Copyright (C) 2022-2026 ABINIT group (MR)
1052 : !! This file is distributed under the terms of the
1053 : !! GNU General Public License, see ~abinit/COPYING
1054 : !! or http://www.gnu.org/copyleft/gpl.txt .
1055 : !!
1056 : !! INPUTS
1057 : !! blkflg(3,mpert,3,mpert) = flags for each element of the 3DTE
1058 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
1059 : !! optgeom= if 1 do special treatment for the geometric term
1060 : !! mpert =maximum number of ipert
1061 : !! natom = number of atoms in unit cell
1062 : !! rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
1063 : !! that have to be computed
1064 : !! 1 -> element has to be computed explicitely
1065 : !! -1 -> use symmetry operations to obtain the corresponding element
1066 : !! rprimd(3,3)=dimensional primitive translations (bohr)
1067 : !! t_typeI(2,3,mpert,3,3,3)= Input type-I tensor
1068 : !!
1069 : !! OUTPUT
1070 : !! t_typeII(2,3,mpert,3,mpert,3,mpert)= type-II tensor converted to the mixed
1071 : !! coordinates.
1072 : !!
1073 : !! SOURCE
1074 :
1075 12 : subroutine dfptlw_typeIproc(blkflg,gprimd,optgeom,mpert,natom,rfpert,rprimd,t_typeI,&
1076 6 : & t_typeII)
1077 :
1078 : use defs_basis
1079 : use m_errors
1080 : use m_profiling_abi
1081 : use m_dynmat, only : cart39
1082 :
1083 : !Arguments ------------------------------------
1084 : !scalars
1085 : integer,intent(in) :: mpert,natom,optgeom
1086 : !arrays
1087 : integer,intent(in) :: blkflg(3,mpert,3,mpert,3,mpert)
1088 : integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert)
1089 : real(dp),intent(in) :: gprimd(3,3),rprimd(3,3)
1090 : real(dp),intent(inout) :: t_typeI(2,3,mpert,3,3,3)
1091 : real(dp),intent(inout) :: t_typeII(2,3,mpert,3,mpert,3,mpert)
1092 :
1093 : !Local variables-------------------------------
1094 : !scalar
1095 : integer :: beta,delta,gamma,ii
1096 : integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,istr
1097 : real(dp) :: fac
1098 : !arrays
1099 : integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/)
1100 : integer :: flg1(3),flg2(3)
1101 : real(dp) :: vec1(3),vec2(3)
1102 : real(dp) :: rprimd_t(3,3),gprimd_t(3,3)
1103 :
1104 : ! *************************************************************************
1105 :
1106 : DBG_ENTER("COLL")
1107 :
1108 6 : if (optgeom==1) then
1109 : !Transform the metric perturbation direction
1110 : !(treat it as an atomic displacement)
1111 12 : flg1(:)=1
1112 9 : do i1pert=1,natom
1113 27 : do i1dir=1,3
1114 78 : do gamma=1,3
1115 180 : do ii=1,2
1116 486 : do delta=1,3
1117 1296 : do beta=1,3
1118 1296 : vec1(beta)=t_typeI(ii,i1dir,i1pert,beta,delta,gamma)
1119 : end do
1120 324 : call cart39(flg1,flg2,gprimd,i1pert,natom,rprimd,vec1,vec2)
1121 1404 : do beta=1,3
1122 1296 : t_typeI(ii,i1dir,i1pert,beta,delta,gamma)=vec2(beta)
1123 : end do
1124 : end do
1125 : end do
1126 : end do
1127 : end do
1128 : end do
1129 :
1130 : !Transform the second q-gradient direction
1131 : !(treat it as an electric field)
1132 9 : do i1pert=1,natom
1133 27 : do i1dir=1,3
1134 78 : do gamma=1,3
1135 180 : do ii=1,2
1136 486 : do beta=1,3
1137 1296 : do delta=1,3
1138 1296 : vec1(delta)=t_typeI(ii,i1dir,i1pert,beta,delta,gamma)
1139 : end do
1140 324 : call cart39(flg1,flg2,gprimd,natom+2,natom,rprimd,vec1,vec2)
1141 1404 : do delta=1,3
1142 1296 : t_typeI(ii,i1dir,i1pert,beta,delta,gamma)=vec2(delta)
1143 : end do
1144 : end do
1145 : end do
1146 : end do
1147 : end do
1148 : end do
1149 :
1150 : !Transform the first q-gradient direction
1151 : !(treat it as an electric field)
1152 9 : do i1pert=1,natom
1153 27 : do i1dir=1,3
1154 60 : do ii=1,2
1155 162 : do beta=1,3
1156 468 : do delta=1,3
1157 1296 : do gamma=1,3
1158 1296 : vec1(gamma)=t_typeI(ii,i1dir,i1pert,beta,delta,gamma)
1159 : end do
1160 324 : call cart39(flg1,flg2,gprimd,natom+2,natom,rprimd,vec1,vec2)
1161 1404 : do gamma=1,3
1162 1296 : t_typeI(ii,i1dir,i1pert,beta,delta,gamma)=vec2(gamma)
1163 : end do
1164 : end do
1165 : end do
1166 : end do
1167 : end do
1168 : end do
1169 :
1170 : end if
1171 :
1172 6 : fac=two_pi ** 2
1173 6 : i3pert= natom+8
1174 78 : rprimd_t=transpose(rprimd)
1175 78 : gprimd_t=transpose(gprimd)
1176 66 : do i1pert = 1, mpert
1177 246 : do i1dir = 1, 3
1178 223440 : if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then
1179 162 : do i2pert = natom+3, natom+4
1180 486 : do i2dir = 1, 3
1181 324 : istr=(i2pert-natom-3)*3+i2dir
1182 324 : beta=idx(2*istr-1); delta=idx(2*istr)
1183 1080 : do ii=1,2
1184 :
1185 : !Transform into type-II
1186 2592 : do i3dir=1,3
1187 1944 : gamma=i3dir
1188 : t_typeII(ii,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)= &
1189 : & t_typeI(ii,i1dir,i1pert,beta,delta,gamma) + &
1190 : & t_typeI(ii,i1dir,i1pert,delta,gamma,beta) - &
1191 2592 : & t_typeI(ii,i1dir,i1pert,gamma,beta,delta)
1192 : end do ! i3dir
1193 :
1194 : !Transform i3dir into reduced coordinates
1195 2592 : do i3dir=1,3
1196 1944 : vec1(i3dir)=t_typeII(ii,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1197 2592 : flg1(i3dir)=blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1198 : end do
1199 648 : call cart39(flg1,flg2,rprimd_t,natom+2,natom,gprimd_t,vec1,vec2)
1200 2916 : do i3dir=1,3
1201 2592 : t_typeII(ii,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=vec2(i3dir)*fac
1202 : end do
1203 :
1204 : end do ! ii
1205 : end do ! i2dir
1206 : end do ! i2pert
1207 : end if ! rfpert
1208 : end do ! i1dir
1209 : end do ! i1pert
1210 :
1211 : DBG_EXIT("COLL")
1212 :
1213 6 : end subroutine dfptlw_typeIproc
1214 : !!***
1215 :
1216 : !!****f* ABINIT/m_dfptlw_loop/read_1eig
1217 : !! NAME
1218 : !! read_1eig
1219 : !!
1220 : !! FUNCTION
1221 : !!
1222 : !! Reads all the first-order energies from a given _1WF file.
1223 : !! Data is read by the master and broadcasted.
1224 : !!
1225 : !! COPYRIGHT
1226 : !! Copyright (C) 2022-2026 ABINIT group (MR)
1227 : !! This file is distributed under the terms of the
1228 : !! GNU General Public License, see ~abinit/COPYING
1229 : !! or http://www.gnu.org/copyleft/gpl.txt .
1230 : !!
1231 : !! INPUTS
1232 : !! formeig option (format of the eigenvalues and occupations) :
1233 : !! 0 => ground-state format
1234 : !! 1 => respfn format
1235 : !! mband=maximum number of bands
1236 : !! nkpt= number of k points
1237 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1238 : !! wffnm=name (character data) of file for input wavefunctions.
1239 : !!
1240 : !! OUTPUT
1241 : !! eigen(2*mband*mband*nkpt*nsppol)=matrix of eigenvalues
1242 : !!
1243 : !! SOURCE
1244 :
1245 265 : subroutine read_1eig(eigen,formeig,mband,nkpt,nsppol,wffnm)
1246 :
1247 : !Arguments ------------------------------------
1248 : !scalars
1249 : integer,intent(in) :: formeig,mband,nkpt,nsppol
1250 : character(len=*),intent(inout) :: wffnm
1251 : !arrays
1252 : real(dp),intent(out) :: eigen((2*mband)**formeig*mband*nkpt*nsppol)
1253 :
1254 : !Local variables-------------------------------
1255 : !scalar
1256 : integer :: bd2tot,comm,ierr,ik_bz,iomode,isppol,master,my_rank
1257 1060 : type(wfk_t) :: Wfk1
1258 265 : type(hdr_type) :: hdr1
1259 : character(len=500) :: msg
1260 : !arrays
1261 265 : real(dp),allocatable :: eig_buffer(:)
1262 :
1263 : ! *************************************************************************
1264 :
1265 : DBG_ENTER("COLL")
1266 :
1267 265 : comm = xmpi_world
1268 265 : master = 0
1269 265 : my_rank = xmpi_comm_rank(comm)
1270 :
1271 : ! Master opens the 1WF file
1272 265 : if (my_rank == master) then
1273 265 : iomode = iomode_from_fname(wffnm)
1274 :
1275 : !Check that atdis file exists and open it
1276 265 : if (.not. file_exists(wffnm)) then
1277 : ! Trick needed to run Abinit test suite in netcdf mode.
1278 265 : if (file_exists(nctk_ncify(wffnm))) then
1279 265 : write(msg,"(3a)")"- File: ",trim(wffnm),&
1280 530 : " does not exist but found netcdf file with similar name."
1281 265 : call wrtout(std_out,msg,'COLL')
1282 530 : wffnm = nctk_ncify(wffnm)
1283 : end if
1284 265 : if (.not. file_exists(wffnm)) then
1285 0 : ABI_ERROR('Missing file: '//TRIM(wffnm))
1286 : end if
1287 : end if
1288 265 : write(msg,'(a,a)')'-open 1wf file :',trim(wffnm)
1289 265 : call wrtout(std_out,msg,'COLL')
1290 :
1291 265 : call wfk1%open_read(wffnm, formeig, iomode, get_unit(), xmpi_comm_self, Hdr_out=hdr1)
1292 : end if
1293 :
1294 : ! Master Broadcasts the header to all procs in comm
1295 265 : call hdr1%bcast(master, my_rank, comm)
1296 :
1297 : !Allocate buffer for MPI communicatio with max dimensions.
1298 795 : ABI_MALLOC(eig_buffer, ((2*mband)**formeig*mband*nsppol))
1299 :
1300 265 : bd2tot = 0
1301 530 : do isppol=1,nsppol
1302 30658 : do ik_bz=1,hdr1%nkpt
1303 :
1304 : !Master reads and broadcasts
1305 30128 : if (my_rank == master) then
1306 90384 : call Wfk1%read_band_block([1,mband], ik_bz, isppol, xmpio_single, eig_k=eig_buffer)
1307 : end if
1308 :
1309 30128 : call xmpi_bcast(eig_buffer, master, comm, ierr)
1310 :
1311 1049008 : eigen(1+bd2tot:2*mband**2+bd2tot)=eig_buffer(:)
1312 :
1313 : !Keep track of total number of bands
1314 60521 : bd2tot = bd2tot + 2*mband**2
1315 : end do
1316 : end do
1317 :
1318 265 : if (my_rank == master) call Wfk1%close()
1319 265 : call hdr1%free()
1320 :
1321 265 : ABI_FREE(eig_buffer)
1322 :
1323 : DBG_EXIT("COLL")
1324 :
1325 536 : end subroutine read_1eig
1326 : !!***
1327 :
1328 : end module m_dfptlw_loop
1329 : !!***
|
