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