Line data Source code
1 : !!****m* ABINIT/m_longwave
2 : !! NAME
3 : !! m_longwave
4 : !!
5 : !! FUNCTION
6 : !! DFPT calculation of spatial and time dispersion properties
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2019-2026 ABINIT group (MR, MS)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! NOTES
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : module m_longwave
25 :
26 : use defs_basis
27 : use m_errors
28 : use m_xmpi
29 : use defs_datatypes
30 : use defs_abitypes, only : MPI_type
31 : use defs_wvltypes
32 : use m_dtset
33 : use m_dtfil
34 : use m_xcdata
35 : use m_hdr
36 : use m_ebands
37 : use m_wffile
38 :
39 : use m_pspini, only : pspini
40 : use m_common, only : setup1
41 : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
42 : use m_pawrhoij, only : pawrhoij_type
43 : use m_paw_dmft, only : paw_dmft_type
44 : use m_pawrad, only : pawrad_type
45 : use m_pawtab, only : pawtab_type
46 : use m_drivexc, only : check_kxc
47 : use m_rhotoxc, only : rhotoxc
48 : use m_ioarr, only : read_rhor
49 : use m_matrix, only : matr3inv
50 : use m_symtk, only : symmetrize_xred
51 : use m_kg, only : kpgio,getcut,getph
52 : use m_inwffil, only : inwffil
53 : use m_spacepar, only : setsym
54 : use m_mkrho, only : mkrho
55 : use m_fft, only : fourdp
56 : use m_ddb, only : ddb_type,lwcart
57 : use m_ddb_hdr, only : ddb_hdr_type
58 : use m_mkcore, only : mkcore
59 : use m_dfptlw_loop, only : dfptlw_loop
60 : use m_dfptlw_nv, only : dfptlw_nv
61 : use m_dfptlw_pert, only : preca_ffnl
62 : use m_initylmg, only : initylmg
63 : use m_geometry, only : symredcart, d3lwsym, sylwtens
64 : use m_atm2fft, only : atm2fft
65 :
66 : implicit none
67 :
68 : private
69 : !!***
70 :
71 : public :: longwave
72 : !!***
73 :
74 : ! *************************************************************************
75 :
76 : contains
77 : !!***
78 :
79 : !!****f* ABINIT/longwave
80 : !! NAME
81 : !! longwave
82 : !!
83 : !! FUNCTION
84 : !! Primary routine for conducting DFPT calculations of dispersion properties
85 : !!
86 : !! INPUTS
87 : !! codvsn = code version
88 : !! dtfil <type(datafiles_type)> = variables related to files
89 : !! dtset <type(dataset_type)> = all input variables for this dataset
90 : !! etotal = new total energy (no meaning at output)
91 : !! mpi_enreg=information about MPI pnarallelization
92 : !! occ(mband*nkpt*nsppol) = occupation number for each band and k
93 : !! xred(3,natom) = reduced atomic coordinates
94 : !!
95 : !! OUTPUT
96 : !! npwtot(nkpt) = total number of plane waves at each k point
97 : !!
98 : !! SIDE EFFECTS
99 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
100 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
101 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
102 : !!
103 : !! SOURCE
104 :
105 17 : subroutine longwave(codvsn,dtfil,dtset,etotal,mpi_enreg,npwtot,occ,&
106 17 : pawrad,pawtab,psps,xred)
107 :
108 : #ifdef FC_INTEL
109 : !DEC$ NOOPTIMIZE
110 : #endif
111 :
112 : !Arguments ------------------------------------
113 : !scalars
114 : real(dp),intent(inout) :: etotal
115 : character(len=8),intent(in) :: codvsn
116 : type(MPI_type),intent(inout) :: mpi_enreg
117 : type(datafiles_type),intent(in) :: dtfil
118 : type(dataset_type),intent(inout) :: dtset
119 : type(pseudopotential_type),intent(inout) :: psps
120 : !arrays
121 : integer,intent(out) :: npwtot(dtset%nkpt)
122 : real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),xred(3,dtset%natom)
123 : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
124 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
125 :
126 : !Local variables-------------------------------
127 : !scalars
128 : integer,parameter :: cplex1=1,formeig=0,response=1
129 : integer :: ask_accurate,bantot,dimffnl,dimffnl_i
130 : integer :: gscase,iatom,ierr,indx,ireadwf0,iscf_eff,itypat
131 : integer :: ider,idir0,idir
132 : integer :: i1dir,i1pert,i2dir,ii,i2pert,i3dir,i3pert
133 : integer :: me,mcg,mgfftf,natom,nfftf,nfftot,nfftotf,nhatdim,nhatgrdim
134 : ! integer :: isym
135 : integer :: mpert,my_natom,n1,nkxc,nk3xc,ntypat,n3xccc,nylmgr
136 : integer :: optatm,optdyfr,opteltfr,optgr,optstr,optv,optn,optn2
137 : integer :: option,optorth,psp_gencond,rdwrpaw,spaceworld,timrev,tim_mkrho
138 : integer :: usexcnhat,useylmgr
139 : real(dp) :: bigexc,bigsxc,boxcut,ecore,ecutdg_eff,ecut_eff,etot
140 : real(dp) :: fermie,fermih,gsqcut,gsqcut_eff,gsqcutc_eff,residm
141 : real(dp) :: ucvol,vxcavg
142 : logical :: non_magnetic_xc,just_timdisp
143 : ! logical :: has_strain
144 : character(len=500) :: msg
145 17 : type(ebands_t) :: bstruct
146 884 : type(ddb_hdr_type) :: ddb_hdr
147 17 : type(ddb_type) :: ddb
148 17 : type(paw_dmft_type) :: paw_dmft
149 : type(pawfgr_type) :: pawfgr
150 17 : type(hdr_type) :: hdr,hdr_den
151 : type(xcdata_type) :: xcdata
152 17 : type(wvl_data) :: wvl
153 : type(wffile_type) :: wffgs,wfftgs
154 : !arrays
155 : integer :: d3e_dir1(3),d3e_dir2(3),d3e_dir3(3)
156 : integer :: ngfft(18),ngfftf(18),perm(6)
157 : real(dp) :: dummy6(6),other_dummy6(6),gmet(3,3),gmet_for_kg(3,3),gprimd(3,3),gprimd_for_kg(3,3)
158 : real(dp) :: rmet(3,3),rprimd(3,3),rprimd_for_kg(3,3)
159 : real(dp) :: dum_gauss(0),dum_dyfrv(0),dum_eltfrxc(0)
160 : real(dp) :: dum_grn(0),dum_grv(0),dum_rhog(0),dum_vg(0)
161 17 : integer,allocatable :: atindx(:),atindx1(:)
162 17 : integer,allocatable :: blkflg(:,:,:,:,:,:),blkflg_car(:,:,:,:,:,:)
163 17 : integer,allocatable :: d3e_pert1(:),d3e_pert2(:),d3e_pert3(:)
164 17 : integer,allocatable :: indsym(:,:,:),irrzon(:,:,:),kg(:,:)
165 17 : integer,allocatable :: nattyp(:),npwarr(:),symrec(:,:,:)
166 17 : integer,allocatable :: rfpert(:,:,:,:,:,:)
167 17 : real(dp),allocatable :: cg(:,:)
168 17 : real(dp),allocatable :: d3etot(:,:,:,:,:,:,:),d3etot_car(:,:,:,:,:,:,:)
169 17 : real(dp),allocatable :: d3etot_nv(:,:,:,:,:,:,:),doccde(:)
170 17 : real(dp),allocatable :: eigen0(:),ffnl(:,:,:,:,:),ffnl_i(:,:,:,:,:)
171 17 : real(dp),allocatable :: grxc(:,:),kxc(:,:),vxc(:,:)
172 17 : real(dp),allocatable :: ncorespl(:,:,:),nhat(:,:),nhatgr(:,:,:)
173 17 : real(dp),allocatable :: phnons(:,:,:),ph1d(:,:),rhog(:,:),rhor(:,:),dummy_dyfrx2(:,:,:)
174 : ! real(dp),allocatable :: symrel_cart(:,:,:)
175 17 : real(dp),allocatable :: dummy_vpsp(:),work(:),xccc3d(:)
176 17 : real(dp),allocatable :: ylm(:,:),ylmgr(:,:,:)
177 17 : type(pawrhoij_type),allocatable :: pawrhoij(:),pawrhoij_read(:)
178 : ! *************************************************************************
179 :
180 : DBG_ENTER("COLL")
181 :
182 : !Just time dispersion calculation
183 17 : just_timdisp=.false.
184 : if (dtset%lw_flexo==0.and.dtset%lw_qdrpl==0.and.&
185 17 : & dtset%lw_natopt==0.and.dtset%timdisp==1) just_timdisp=.true.
186 :
187 : !Not valid for PAW
188 17 : if (psps%usepaw==1) then
189 0 : msg='This routine cannot be used for PAW!'
190 0 : ABI_BUG(msg)
191 : end if
192 :
193 : !Not valid for finite wave-vector perturbations
194 68 : if (.not.just_timdisp.and.sqrt(sum(dtset%qptn**2))>tol8) then
195 0 : msg='This routine can only be used at q /= 0 for time dispersion calculation '
196 0 : ABI_BUG(msg)
197 : end if
198 :
199 : !Only usable with spherical harmonics
200 17 : if (dtset%useylm/=1.and.(.not.just_timdisp.and.dtset%lw_qdrpl/=0.or.dtset%lw_flexo/=0)) then
201 0 : msg='This routine can only be used with useylm/=1 for the spatial dispersion given by lw_natopt=1'
202 0 : ABI_BUG(msg)
203 : end if
204 :
205 : !Not valid for spin-dependent calculations
206 17 : if (.not.just_timdisp.and.(dtset%nspinor/=1.or.dtset%nsppol/=1.or.dtset%nspden/=1)) then
207 0 : msg='This routine cannot be used for spin-dependent calculations of spatial dispersion'
208 0 : ABI_BUG(msg)
209 : end if
210 :
211 : !Not usable with core electron density corrections and flexo
212 17 : if (.not.just_timdisp.and.psps%n1xccc/=0.and.dtset%lw_flexo/=0) then
213 0 : msg='This routine cannot be used to calculate flexoelectric properties with n1xccc/=0'
214 0 : ABI_BUG(msg)
215 : end if
216 :
217 : !Only usable for insulators
218 : if(0>dtset%occopt .and. dtset%occopt>2)then
219 : msg='This routine cannot be used with metallic occupations'
220 : end if
221 :
222 : !Define some data
223 17 : ntypat=psps%ntypat
224 17 : natom=dtset%natom
225 17 : timrev=1
226 :
227 : !Init spaceworld
228 17 : spaceworld=mpi_enreg%comm_cell
229 17 : my_natom=mpi_enreg%my_natom
230 17 : me = xmpi_comm_rank(spaceworld)
231 :
232 : !Define FFT grid(s) sizes (be careful !)
233 : !See NOTES in the comments at the beginning of this file.
234 17 : call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
235 68 : nfftot=product(ngfft(1:3))
236 : nfftotf=product(ngfftf(1:3))
237 :
238 : !Set up for iterations
239 : call setup1(dtset%acell_orig(1:3,1),bantot,dtset,&
240 : & ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,&
241 : & ngfftf,ngfft,dtset%nkpt,dtset%nsppol,&
242 17 : & response,rmet,dtset%rprim_orig(1:3,1:3,1),rprimd,ucvol,psps%usepaw)
243 :
244 : !Define the set of admitted perturbations taking into account
245 : !the possible permutations
246 : ! -> natom+8 refers to ddq perturbation (spatial dispersion)
247 : ! -> natom+9 refers to ddw perturbation (time dispersion)
248 17 : if (dtset%timdisp==0) then
249 11 : mpert=natom+8
250 : else
251 6 : mpert=natom+9
252 6 : if (dtset%d3e_pert1_magn==2.or.dtset%d3e_pert2_magn==2) mpert=2*dtset%natom+11
253 : end if
254 119 : ABI_MALLOC(blkflg,(3,mpert,3,mpert,3,mpert))
255 119 : ABI_MALLOC(d3etot,(2,3,mpert,3,mpert,3,mpert))
256 68 : ABI_MALLOC(d3etot_nv,(2,3,mpert,3,mpert,3,mpert))
257 68 : ABI_MALLOC(rfpert,(3,mpert,3,mpert,3,mpert))
258 51 : ABI_MALLOC(d3e_pert1,(mpert))
259 34 : ABI_MALLOC(d3e_pert2,(mpert))
260 34 : ABI_MALLOC(d3e_pert3,(mpert))
261 2623081 : blkflg(:,:,:,:,:,:) = 0
262 6485593 : d3etot(:,:,:,:,:,:,:) = zero
263 6485593 : d3etot_nv(:,:,:,:,:,:,:) = zero
264 2623081 : rfpert(:,:,:,:,:,:) = 0
265 743 : d3e_pert1(:) = 0 ; d3e_pert2(:) = 0 ; d3e_pert3(:) = 0
266 :
267 : !Spatial-dispersion perturbations
268 17 : if (.not.just_timdisp) d3e_pert3(natom+8)=1
269 :
270 17 : if (dtset%lw_qdrpl==1) then
271 7 : d3e_pert1(natom+2)=1
272 24 : d3e_pert2(1:natom)=1
273 : end if
274 :
275 17 : if (dtset%lw_flexo==2.or.dtset%lw_flexo==1) then
276 3 : d3e_pert1(natom+2)=1
277 9 : d3e_pert2(natom+3:natom+4)=1
278 : end if
279 :
280 17 : if (dtset%lw_flexo==3.or.dtset%lw_flexo==1) then
281 9 : d3e_pert1(natom+2)=1 ; d3e_pert1(1:natom)=1
282 9 : d3e_pert2(1:natom)=1
283 : end if
284 :
285 17 : if (dtset%lw_flexo==4.or.dtset%lw_flexo==1) then
286 9 : d3e_pert1(1:natom)=1
287 9 : d3e_pert2(natom+3:natom+4)=1
288 : end if
289 :
290 17 : if (dtset%lw_natopt==1) then
291 1 : d3e_pert1(natom+2)=1
292 1 : d3e_pert2(natom+2)=1
293 : end if
294 :
295 17 : perm(:)=0
296 259 : do i1pert = 1, mpert
297 4183 : do i2pert = 1, mpert
298 75694 : do i3pert = 1, mpert
299 71528 : perm(1)=d3e_pert1(i1pert)*d3e_pert2(i2pert)*d3e_pert3(i3pert)
300 : ! perm(2)=d3e_pert1(i1pert)*d3e_pert2(i3pert)*d3e_pert3(i2pert)
301 : ! perm(3)=d3e_pert1(i2pert)*d3e_pert2(i1pert)*d3e_pert3(i3pert)
302 : ! perm(4)=d3e_pert1(i2pert)*d3e_pert2(i3pert)*d3e_pert3(i1pert)
303 : ! perm(5)=d3e_pert1(i3pert)*d3e_pert2(i2pert)*d3e_pert3(i1pert)
304 : ! perm(6)=d3e_pert1(i3pert)*d3e_pert2(i1pert)*d3e_pert3(i2pert)
305 506726 : if ( sum(perm(:)) > 0 ) rfpert(:,i1pert,:,i2pert,:,i3pert)=1
306 : end do
307 : end do
308 : end do
309 :
310 : !Time-dispersion perturbations
311 : ! d3e_pert1(:) = 0 ; d3e_pert2(:) = 0 ; d3e_pert3(:) = 0
312 : ! d3e_dir1(:) = 0 ; d3e_dir2(:) = 0 ; d3e_dir3(:) = 0
313 17 : if (dtset%timdisp==1) then
314 6 : d3e_pert3(natom+9)=1
315 6 : d3e_dir3(1)=1
316 6 : d3e_dir3(2)=0
317 6 : d3e_dir3(3)=0
318 :
319 : !Atomic displacements
320 10 : if (dtset%d3e_pert1_phon==1) d3e_pert1(dtset%d3e_pert1_atpol(1):dtset%d3e_pert1_atpol(2))=1
321 10 : if (dtset%d3e_pert2_phon==1) d3e_pert2(dtset%d3e_pert2_atpol(1):dtset%d3e_pert2_atpol(2))=1
322 :
323 : !Electric fields
324 6 : if (dtset%d3e_pert1_elfd/=0) d3e_pert1(natom+2)=1
325 6 : if (dtset%d3e_pert2_elfd/=0) d3e_pert2(natom+2)=1
326 :
327 : !Macroscopic hspinfields
328 6 : if (dtset%d3e_pert1_magn==1) d3e_pert1(natom+5)=1
329 6 : if (dtset%d3e_pert2_magn==1) d3e_pert2(natom+5)=1
330 :
331 : !Local hspinfields
332 9 : if (dtset%d3e_pert1_magn==2) d3e_pert1(natom+11+dtset%d3e_pert1_magat(1):natom+11+dtset%d3e_pert1_magat(2))=1
333 9 : if (dtset%d3e_pert2_magn==2) d3e_pert2(natom+11+dtset%d3e_pert2_magat(1):natom+11+dtset%d3e_pert2_magat(2))=1
334 :
335 6 : perm(:)=0
336 128 : do i1pert = 1, mpert
337 488 : d3e_dir1(:)=dtset%d3e_pert1_dir(:)
338 122 : if (i1pert>natom+11.and.i1pert<=2*natom+11) &
339 96 : & d3e_dir1(:)=dtset%d3e_pert1_magdir(:)
340 494 : do i1dir = 1, 3
341 8186 : do i2pert = 1, mpert
342 30792 : d3e_dir2(:)=dtset%d3e_pert2_dir(:)
343 7698 : if (i2pert>natom+11.and.i2pert<=2*natom+11) &
344 6624 : & d3e_dir2(:)=dtset%d3e_pert2_magdir(:)
345 31158 : do i2dir = 1, 3
346 53886 : do i3pert = natom+9, natom+9
347 115470 : do i3dir = 1, 3
348 : perm(1)=d3e_pert1(i1pert)*d3e_dir1(i1dir)* &
349 : & d3e_pert2(i2pert)*d3e_dir2(i2dir)* &
350 69282 : & d3e_pert3(i3pert)*d3e_dir3(i3dir)
351 508068 : if ( sum(perm(:)) > 0 ) rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=1
352 : end do
353 : end do
354 : end do
355 : end do
356 : end do
357 : end do
358 : end if
359 :
360 : !Do symmetry stuff
361 68 : ABI_MALLOC(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
362 68 : ABI_MALLOC(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
363 68 : ABI_MALLOC(indsym,(4,dtset%nsym,natom))
364 51 : ABI_MALLOC(symrec,(3,3,dtset%nsym))
365 715682 : irrzon=0;indsym=0;symrec=0;phnons=zero
366 : !If the density is to be computed by mkrho, need irrzon and phnons
367 17 : iscf_eff=0;if(dtset%getden==0)iscf_eff=1
368 : call setsym(indsym,irrzon,iscf_eff,natom,&
369 : & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
370 17 : & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred)
371 :
372 : !Symmetrize atomic coordinates over space group elements:
373 17 : call symmetrize_xred(natom,dtset%nsym,dtset%symrel,dtset%tnons,xred,indsym=indsym)
374 :
375 : ! Get symmetries in cartesian coordinates
376 : ! ABI_MALLOC(symrel_cart, (3, 3, dtset%nsym))
377 : ! do isym =1,dtset%nsym
378 : ! call symredcart(rprimd, gprimd, symrel_cart(:,:,isym), dtset%symrel(:,:,isym))
379 : ! ! purify operations in cartesian coordinates.
380 : ! where (abs(symrel_cart(:,:,isym)) < tol14)
381 : ! symrel_cart(:,:,isym) = zero
382 : ! end where
383 : ! end do
384 :
385 : ! call sylwtens(indsym,mpert,natom,dtset%nsym,rfpert,symrec,dtset%symrel,symrel_cart)
386 17 : call sylwtens(indsym,mpert,natom,dtset%nsym,rfpert,symrec,dtset%symrel)
387 :
388 17 : write(msg,'(a,a,a,a,a)') ch10, &
389 17 : & ' The list of irreducible elements of the spatial-dispersion third-order energy derivatives is: ', ch10,&
390 34 : & ' (in reduced coordinates except for strain and hspinfield perturbations) ', ch10
391 17 : call wrtout(ab_out,msg,'COLL')
392 17 : call wrtout(std_out,msg,'COLL')
393 :
394 17 : write(msg,'(12x,a)') 'i1dir i1pert i2dir i2pert i3dir i3pert'
395 17 : call wrtout(ab_out,msg,'COLL')
396 17 : call wrtout(std_out,msg,'COLL')
397 17 : n1 = 0
398 259 : do i3pert = 1, mpert
399 985 : do i3dir = 1, 3
400 12740 : do i2pert = 1, mpert
401 47814 : do i2dir = 1,3
402 690840 : do i1pert = 1, mpert
403 2610324 : do i1dir = 1, 3
404 2575008 : if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
405 611 : n1 = n1 + 1
406 611 : write(msg,'(2x,i4,a,6(5x,i3))') n1,')', &
407 1222 : & i1dir,i1pert,i2dir,i2pert,i3dir,i3pert
408 611 : call wrtout(ab_out,msg,'COLL')
409 611 : call wrtout(std_out,msg,'COLL')
410 1930645 : else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-2) then
411 263 : blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
412 263 : if (dtset%prtvol>=10) then
413 0 : n1 = n1 + 1
414 0 : write(msg,'(2x,i4,a,6(5x,i3),a)') n1,')', &
415 0 : & i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,' => must be zero, not computed'
416 0 : call wrtout(ab_out,msg,'COLL')
417 0 : call wrtout(std_out,msg,'COLL')
418 : end if
419 1930382 : else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-1) then
420 612 : if (dtset%prtvol>=10) then
421 0 : n1 = n1 + 1
422 0 : write(msg,'(2x,i4,a,6(5x,i3),a)') n1,')', &
423 0 : & i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,' => symmetric of another element, not computed'
424 0 : call wrtout(ab_out,msg,'COLL')
425 0 : call wrtout(std_out,msg,'COLL')
426 : end if
427 : end if
428 : end do
429 : end do
430 : end do
431 : end do
432 : end do
433 : end do
434 17 : write(msg,'(a,a)') ch10,ch10
435 17 : call wrtout(ab_out,msg,'COLL')
436 17 : call wrtout(std_out,msg,'COLL')
437 :
438 : !In some cases (e.g. getcell/=0), the plane wave vectors have
439 : !to be generated from the original simulation cell
440 17 : rprimd_for_kg=rprimd
441 17 : if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=dtset%rprimd_orig(:,:,1)
442 17 : call matr3inv(rprimd_for_kg,gprimd_for_kg)
443 680 : gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg)
444 :
445 : !Set up the basis sphere of planewaves
446 51 : ABI_MALLOC(kg,(3,dtset%mpw*dtset%mkmem))
447 51 : ABI_MALLOC(npwarr,(dtset%nkpt))
448 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet_for_kg,dtset%istwfk,kg,&
449 : & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw,npwarr,npwtot,&
450 17 : & dtset%nsppol)
451 :
452 : !Open and read pseudopotential files
453 : ecore=zero
454 : call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,&
455 17 : & psps,rprimd,comm_mpi=mpi_enreg%comm_cell)
456 :
457 : !Initialize band structure datatype
458 17 : call bstruct%from_dtset(dtset, npwarr)
459 :
460 : !Initialize PAW atomic occupancies to zero
461 17 : ABI_MALLOC(pawrhoij,(0))
462 :
463 : !Initialize header
464 17 : gscase=0
465 : call hdr%init(bstruct,codvsn,dtset,pawtab,gscase,psps,wvl%descr, &
466 17 : & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
467 :
468 : !Update header, with evolving variables, when available
469 : !Here, rprimd, xred and occ are available
470 17 : etot=hdr%etot ; fermie=hdr%fermie ; fermih=hdr%fermih ; residm=hdr%residm
471 :
472 : !If parallelism over atom, hdr is distributed
473 : call hdr%update(bantot,etot,fermie,fermih,&
474 : & residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1), &
475 17 : & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
476 :
477 : !Clean band structure datatype (should use it more in the future !)
478 17 : call bstruct%free()
479 :
480 : !Initialize wavefunction files and wavefunctions.
481 17 : mcg=dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
482 51 : ABI_STAT_MALLOC(cg,(2,mcg), ierr)
483 17 : ABI_CHECK(ierr==0, "out-of-memory in cg")
484 :
485 51 : ABI_MALLOC(eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol))
486 7121 : eigen0(:)=zero ; ask_accurate=1
487 17 : optorth=0
488 :
489 17 : ireadwf0=1
490 221 : hdr%rprimd=rprimd_for_kg ! We need the rprimd that was used to generate de G vectors
491 17 : if (.not.just_timdisp) then
492 : call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,&
493 : & formeig,hdr,ireadwf0,dtset%istwfk,kg,dtset%kptns,&
494 : & dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,&
495 : & dtset%nband,ngfft,dtset%nkpt,npwarr,dtset%nsppol,dtset%nsym,&
496 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
497 11 : & dtfil%unkg,wffgs,wfftgs,dtfil%unwffgs,dtfil%fnamewffk,wvl)
498 143 : hdr%rprimd=rprimd
499 :
500 : !Close wffgs, if it was ever opened (in inwffil)
501 11 : if (ireadwf0==1) then
502 11 : call WffClose(wffgs,ierr)
503 : end if
504 : end if
505 :
506 : !Generate an index table of atoms, in order for them to be used
507 : !type after type.
508 51 : ABI_MALLOC(atindx,(natom))
509 34 : ABI_MALLOC(atindx1,(natom))
510 51 : ABI_MALLOC(nattyp,(ntypat))
511 48 : indx=1
512 48 : do itypat=1,ntypat
513 31 : nattyp(itypat)=0
514 181 : do iatom=1,natom
515 164 : if(dtset%typat(iatom)==itypat)then
516 68 : atindx(iatom)=indx
517 68 : atindx1(indx)=iatom
518 68 : indx=indx+1
519 68 : nattyp(itypat)=nattyp(itypat)+1
520 : end if
521 : end do
522 : end do
523 :
524 : !Derivative of occupations is always zero for non metallic systems
525 51 : ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
526 7121 : doccde(:)=zero
527 :
528 : !Read ground-state charge density from diskfile in case getden /= 0
529 : !or compute it from wfs that were read previously : rhor
530 :
531 51 : ABI_MALLOC(rhog,(2,nfftf))
532 68 : ABI_MALLOC(rhor,(nfftf,dtset%nspden))
533 :
534 17 : if (.not.just_timdisp) then
535 11 : if (dtset%getden /= 0 .or. dtset%irdden /= 0) then
536 : ! Read rho1(r) from a disk file and broadcast data.
537 : ! This part is not compatible with MPI-FFT (note single_proc=.True. below)
538 :
539 0 : rdwrpaw=psps%usepaw
540 0 : ABI_MALLOC(pawrhoij_read,(0))
541 :
542 : !
543 : ! MT july 2013: Should we read rhoij from the density file ?
544 : call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor, &
545 0 : hdr_den, pawrhoij_read, spaceworld, check_hdr=hdr)
546 0 : etotal = hdr_den%etot; call hdr_den%free()
547 :
548 0 : ABI_FREE(pawrhoij_read)
549 :
550 : ! Compute up+down rho(G) by fft
551 0 : ABI_MALLOC(work,(nfftf))
552 0 : work(:)=rhor(:,1)
553 0 : call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,1,ngfftf,0)
554 0 : ABI_FREE(work)
555 : else
556 : ! Obtain the charge density from read wfs
557 : ! Be careful: in PAW, compensation density has to be added !
558 11 : tim_mkrho=4
559 11 : paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented
560 11 : paw_dmft%use_dmft=0 ! respfn with dmft not implemented
561 :
562 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
563 11 : & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
564 : end if ! getden
565 : ! ABI_FREE(cg)
566 : end if
567 :
568 : !Compute large sphere cut-off gsqcut
569 17 : call getcut(boxcut,dtset%ecut,gmet,gsqcut,dtset%iboxcut,std_out,dtset%qptn,dtset%ngfft)
570 :
571 : !Generate the 1-dimensional phases
572 51 : ABI_MALLOC(ph1d,(2,3*(2*mgfftf+1)*dtset%natom))
573 17 : call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),ph1d,xred)
574 :
575 : !Pseudo core electron density by method 2
576 : !TODO: The tasks to adapt the code to consider n3xccc in the long-wave
577 : !driver are under way.
578 17 : n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf
579 51 : ABI_MALLOC(xccc3d,(n3xccc))
580 17 : if (.not.just_timdisp) then
581 44 : ABI_MALLOC(ncorespl,(psps%mqgrid_vl,2,ntypat))
582 11 : if (psps%n1xccc/=0) then
583 6 : ABI_MALLOC(dummy_dyfrx2,(3,3,natom)) ! dummy
584 2 : if (psps%nc_xccc_gspace==1) then
585 0 : ABI_MALLOC(dummy_vpsp,(nfftf))
586 0 : optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optv=0;optn=n3xccc/nfftf;optn2=1
587 : call atm2fft(atindx1,xccc3d,dummy_vpsp,dummy_dyfrx2,dum_dyfrv,dum_eltfrxc,dum_gauss,gmet,gprimd,&
588 : & dum_grn,dum_grv,gsqcut,mgfftf,psps%mqgrid_vl,natom,nattyp,nfftf,ngfftf,&
589 : & ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,&
590 0 : & dtset%qprtrb,dtset%rcut,dum_rhog,rprimd,dummy6,other_dummy6,ucvol,psps%usepaw,dum_vg,dum_vg,dum_vg,dtset%vprtrb,psps%vlspl)
591 :
592 0 : ABI_FREE(dummy_vpsp)
593 : end if
594 2 : if (psps%nc_xccc_gspace==0) then
595 2 : option=1
596 2 : ABI_MALLOC(vxc,(0,0)) ! dummy
597 6 : ABI_MALLOC(grxc,(3,natom))
598 : call mkcore(dummy6,dummy_dyfrx2,grxc,mpi_enreg,natom,nfftf,dtset%nspden,ntypat,&
599 : & ngfftf(1),psps%n1xccc,ngfftf(2),ngfftf(3),option,rprimd,dtset%typat,ucvol,vxc,&
600 2 : & psps%xcccrc,psps%xccc1d,xccc3d,xred)
601 2 : ABI_FREE(vxc) ! dummy
602 2 : ABI_FREE(grxc) ! dummy
603 : end if
604 2 : ABI_FREE(dummy_dyfrx2) ! dummy
605 :
606 : !Write the spl interpolation of the pseudo core density for all atom types
607 6 : do itypat= 1, ntypat
608 24022 : ncorespl(:,:,itypat)= psps%nctab(itypat)%tcorespl(:,:)
609 : end do
610 : end if
611 : end if
612 :
613 : !Set up xc potential. Compute kxc here.
614 : !TODO: Iclude nonlinear core corrections (see m_respfn_driver.F90)
615 17 : option=2 ; nk3xc=1
616 17 : nkxc=2*min(dtset%nspden,2)-1;if(dtset%xclevel==2)nkxc=12*min(dtset%nspden,2)-5
617 17 : call check_kxc(dtset%ixc,dtset%optdriver)
618 68 : ABI_MALLOC(kxc,(nfftf,nkxc))
619 68 : ABI_MALLOC(vxc,(nfftf,dtset%nspden))
620 :
621 17 : nhatgrdim=0;nhatdim=0
622 17 : ABI_MALLOC(nhat,(0,0))
623 17 : ABI_MALLOC(nhatgr,(0,0,0))
624 17 : non_magnetic_xc=.false.
625 :
626 17 : usexcnhat=0
627 17 : bigexc=zero
628 17 : bigsxc=zero
629 :
630 17 : if (.not.just_timdisp) then
631 11 : call xcdata_init(xcdata,dtset=dtset)
632 : call rhotoxc(bigexc,bigsxc,kxc,mpi_enreg,nfftf,ngfftf,&
633 : & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,rhor,&
634 11 : & rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata)
635 : end if
636 :
637 17 : ABI_FREE(xccc3d)
638 :
639 : !Set up the spherical harmonics (Ylm) and gradients at each k point
640 17 : if (psps%useylm==1) then
641 11 : useylmgr=1; option=2 ; nylmgr=9
642 44 : ABI_MALLOC(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm))
643 55 : ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
644 11 : if (.not.just_timdisp) then
645 : call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,&
646 : & psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,npwarr,dtset%nsppol,option,&
647 11 : & rprimd,ylm,ylmgr)
648 : end if
649 : end if
650 :
651 : !Compute nonlocal form factors ffnl1, for all atoms and all k-points.
652 17 : if (dtset%ffnl_lw == 0) then
653 16 : if (dtset%lw_natopt==1) then
654 1 : ider=1;dimffnl=4;dimffnl_i=2
655 7 : ABI_MALLOC(ffnl,(dtset%mkmem,dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat))
656 7 : ABI_MALLOC(ffnl_i,(dtset%mkmem,dtset%mpw,dimffnl_i,psps%lmnmax,psps%ntypat))
657 4 : do idir=1, 3
658 3 : idir0=idir
659 3 : if (.not.just_timdisp) then
660 : call preca_ffnl(dimffnl_i,ffnl_i,gmet,gprimd,ider,idir0,kg, &
661 : & dtset%kptns,dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw, &
662 3 : & dtset%nkpt,npwarr,nylmgr,psps,rmet,useylmgr,ylm,ylmgr)
663 : end if
664 284157 : ffnl(:,:,1,:,:)=ffnl_i(:,:,1,:,:)
665 284158 : ffnl(:,:,1+idir,:,:)=ffnl_i(:,:,2,:,:)
666 : end do
667 1 : ABI_FREE(ffnl_i)
668 1 : if (psps%useylm==1) then
669 1 : useylmgr=0
670 1 : ABI_FREE(ylmgr)
671 3 : ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
672 : end if
673 : else
674 15 : if (dtset%lw_qdrpl==1.or.dtset%lw_flexo==3) ider=1; idir0=4; dimffnl=4
675 15 : if (dtset%lw_flexo==1.or.dtset%lw_flexo==2.or.dtset%lw_flexo==4) then
676 2 : ider=2; idir0=4; dimffnl=10
677 : end if
678 105 : ABI_MALLOC(ffnl,(dtset%mkmem,dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat))
679 15 : if (.not.just_timdisp) then
680 : call preca_ffnl(dimffnl,ffnl,gmet,gprimd,ider,idir0,kg, &
681 : & dtset%kptns,dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw, &
682 9 : & dtset%nkpt,npwarr,nylmgr,psps,rmet,useylmgr,ylm,ylmgr)
683 : end if
684 15 : useylmgr=0
685 15 : ABI_SFREE(ylmgr)
686 45 : ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
687 : end if
688 1 : else if (dtset%ffnl_lw == 1) then
689 1 : dimffnl=0
690 5 : ABI_MALLOC(ffnl,(dtset%mkmem,dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat))
691 : end if
692 :
693 : !############# SPATIAL-DISPERSION PROPERTIES CALCULATION ###########################
694 :
695 : !Anounce start of spatial-dispersion calculation
696 1394 : write(msg, '(a,80a,a,a,a)' ) ch10,('=',ii=1,80),ch10,&
697 1411 : & ' ==> Compute spatial-dispersion 3rd-order energy derivatives <== ',ch10
698 17 : call wrtout(std_out,msg,'COLL')
699 17 : call wrtout(ab_out,msg,'COLL')
700 :
701 17 : if (dtset%prtvol>=10) then
702 0 : write(msg,'(5a)') ' CAUTION: Individual contributions to the 3rd-order energy derivatives ',ch10, &
703 0 : & ' are not written in a unified form. Mixed cartesian/reduced coordinates ',ch10, &
704 0 : & ' and/or type-I/type-II forms are used.'
705 0 : call wrtout(std_out,msg,'COLL')
706 0 : call wrtout(ab_out,msg,'COLL')
707 : end if
708 :
709 : !Calculate the nonvariational Ewald terms
710 17 : if (dtset%lw_flexo==1.or.dtset%lw_flexo==3.or.dtset%lw_flexo==4) then
711 : call dfptlw_nv(d3etot_nv,dtset,gmet,gprimd,mpert,my_natom,rfpert,rmet,rprimd,ucvol,xred,psps%ziontypat, &
712 3 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
713 : end if
714 :
715 : !Main loop over the perturbations to calculate the stationary part
716 : call dfptlw_loop(atindx,blkflg,cg,codvsn,d3e_pert1,d3e_pert2,d3etot,dimffnl,dtfil,dtset,&
717 : & ffnl,gmet,gprimd,gsqcut,&
718 : & hdr,just_timdisp,kg,kxc,dtset%mband,dtset%mgfft,&
719 : & dtset%mkmem,dtset%mk1mem,mpert,mpi_enreg,dtset%mpw,natom,nattyp,ncorespl,ngfftf,nfftf,&
720 : & dtset%nkpt,nkxc,dtset%nspinor,dtset%nsppol,npwarr,nylmgr,occ,&
721 : & pawfgr,pawtab,ph1d,&
722 17 : & psps,rfpert,rhog,rhor,rmet,rprimd,ucvol,useylmgr,xred,ylm,ylmgr)
723 :
724 : !Merge stationay and nonvariational contributions
725 6485593 : d3etot(:,:,:,:,:,:,:)=d3etot(:,:,:,:,:,:,:) + d3etot_nv(:,:,:,:,:,:,:)
726 :
727 : !Real (imaginary) part of d3etot is zero for first (second) momentum derivatives
728 17 : if (dtset%kptopt /= 3) then
729 131 : do i3pert = 1, mpert
730 120 : if (i3pert == natom+9 ) cycle
731 491 : do i3dir = 1, 3
732 4554 : do i2pert = 1, mpert
733 16656 : do i2dir = 1,3
734 161286 : do i1pert = 1, mpert
735 592182 : do i1dir = 1, 3
736 579960 : if (blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) == 1) then
737 854 : if (i2pert /= natom+3 .and. i2pert /= natom+4) then
738 368 : d3etot(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero
739 : else
740 486 : d3etot(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero
741 : end if
742 : end if
743 : end do
744 : end do
745 : end do
746 : end do
747 : end do
748 : end do
749 : end if
750 :
751 :
752 : !Complete missing elements using symmetry operations
753 : ! has_strain=.false.
754 : ! if (dtset%lw_flexo==1.or.dtset%lw_flexo==2.or.dtset%lw_flexo==4) has_strain=.true.
755 : ! call d3lwsym(blkflg,d3etot,has_strain,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel,symrel_cart)
756 17 : call d3lwsym(blkflg,d3etot,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
757 :
758 : !Deallocate global proc_distrib
759 17 : if(xmpi_paral==1) then
760 17 : ABI_FREE(mpi_enreg%proc_distrb)
761 : end if
762 :
763 : ! Write the DDB file
764 : call ddb_hdr%init(dtset,psps,pawtab,&
765 : dscrpt=' Note : temporary (transfer) database ',&
766 17 : nblok=1,xred=xred,occ=occ)
767 :
768 17 : call ddb%init(dtset, 1, mpert, with_d3E=.true.)
769 :
770 17 : call ddb%set_omega(iblok=1, omega=dtset%rfomega)
771 :
772 17 : call ddb%set_qpt(iblok=1, qpt=dtset%qptn(1:3))
773 :
774 17 : call ddb%set_d3matr(1, d3etot, blkflg, lw=.true.)
775 :
776 17 : call ddb%write_txt(ddb_hdr, dtfil%fnameabo_ddb)
777 :
778 17 : call ddb_hdr%free()
779 17 : call ddb%free()
780 :
781 : !Calculate spatial-dispersion quantities in Cartesian coordinates and write
782 : !them in abi_out
783 119 : ABI_MALLOC(blkflg_car,(3,mpert,3,mpert,3,mpert))
784 119 : ABI_MALLOC(d3etot_car,(2,3,mpert,3,mpert,3,mpert))
785 17 : call lwcart(blkflg,blkflg_car,d3etot,d3etot_car,gprimd,mpert,natom,rprimd)
786 17 : if (me==0) then
787 17 : call dfptlw_out(blkflg_car,d3etot_car,dtset%lw_flexo,dtset%lw_qdrpl,dtset%lw_natopt,mpert,natom,dtset%timdisp,ucvol)
788 : end if
789 :
790 : !Deallocate arrays
791 17 : ABI_FREE(atindx)
792 17 : ABI_FREE(atindx1)
793 17 : ABI_FREE(blkflg)
794 17 : ABI_FREE(doccde)
795 17 : ABI_FREE(eigen0)
796 17 : ABI_FREE(cg)
797 17 : ABI_FREE(ph1d)
798 17 : ABI_SFREE(ffnl)
799 17 : ABI_FREE(indsym)
800 17 : ABI_FREE(irrzon)
801 17 : ABI_FREE(nattyp)
802 17 : ABI_FREE(kg)
803 17 : ABI_FREE(kxc)
804 17 : ABI_FREE(npwarr)
805 17 : ABI_FREE(phnons)
806 17 : ABI_FREE(rhog)
807 17 : ABI_FREE(rhor)
808 17 : ABI_FREE(symrec)
809 : ! ABI_FREE(symrel_cart)
810 17 : ABI_FREE(vxc)
811 17 : ABI_FREE(d3etot)
812 17 : ABI_FREE(d3etot_nv)
813 : ! ABI_FREE(pertsy)
814 17 : ABI_FREE(rfpert)
815 17 : ABI_FREE(d3e_pert1)
816 17 : ABI_FREE(d3e_pert2)
817 17 : ABI_FREE(d3e_pert3)
818 17 : ABI_SFREE(pawrhoij)
819 17 : ABI_SFREE(nhat)
820 17 : ABI_SFREE(nhatgr)
821 17 : ABI_SFREE(ylm)
822 17 : ABI_SFREE(ylmgr)
823 17 : ABI_SFREE(blkflg_car)
824 17 : ABI_SFREE(d3etot_car)
825 17 : ABI_SFREE(ncorespl)
826 17 : call pawfgr_destroy(pawfgr)
827 :
828 : ! Clean the header
829 17 : call hdr%free()
830 :
831 : DBG_EXIT("COLL")
832 :
833 34 : end subroutine longwave
834 : !!***
835 :
836 : !!****f* ABINIT/m_dfptlw_loop/dfptlw_out
837 : !! NAME
838 : !! dfptlw_out
839 : !!
840 : !! FUNCTION
841 : !! Write the relevant spatial-dispersion quantities in Cartesian coordinates
842 : !!
843 : !! COPYRIGHT
844 : !! Copyright (C) 2022-2026 ABINIT group (MR)
845 : !! This file is distributed under the terms of the
846 : !! GNU General Public License, see ~abinit/COPYING
847 : !! or http://www.gnu.org/copyleft/gpl.txt .
848 : !!
849 : !! INPUTS
850 : !! blkflg_car(3,mpert,3,mpert,3,mpert) =flags for each element of the 3DTE
851 : !! d3etot_car(2,3,mpert,3,mpert,3,mpert) =array with the cartesian thir-order derivatives
852 : !! lw_qdrpl= flag that activates quadrupoles calculation
853 : !! lw_flexo= flag that activates flexoelectric tensor calculation
854 : !! mpert =maximum number of ipert
855 : !! natom = number of atoms in unit cell
856 : !! timdisp = if 1 write derivatives with respect to omega
857 : !!
858 : !! OUTPUT
859 : !!
860 : !! SIDE EFFECTS
861 : !!
862 : !! SOURCE
863 :
864 17 : subroutine dfptlw_out(blkflg_car,d3etot_car,lw_flexo,lw_qdrpl,lw_natopt,mpert,natom,timdisp,ucvol)
865 :
866 : !Arguments ------------------------------------
867 : !scalars
868 : integer,intent(in) :: lw_flexo,lw_qdrpl,lw_natopt,mpert,natom,timdisp
869 : real(dp),intent(in) :: ucvol
870 : !arrays
871 : integer,intent(in) :: blkflg_car(3,mpert,3,mpert,3,mpert)
872 : real(dp),intent(in) :: d3etot_car(2,3,mpert,3,mpert,3,mpert)
873 :
874 : !Local variables-------------------------------
875 : !scalar
876 : integer :: beta,delta,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,istr
877 : !arrays
878 : integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/)
879 17 : real(dp),allocatable :: qdrp(:,:,:,:,:,:,:)
880 : real(dp) :: piezoci(2),piezofr(2),celastci(2)
881 :
882 : ! *************************************************************************
883 :
884 : DBG_ENTER("COLL")
885 :
886 17 : i3pert=natom+8
887 17 : if (lw_qdrpl==1.or.lw_flexo==3.or.lw_flexo==1) then
888 10 : write(ab_out,'(a)')' First real-space moment of the polarization response '
889 10 : write(ab_out,'(a)')' to an atomic displacementatom, in cartesian coordinates,'
890 10 : write(ab_out,'(a)')' (1/ucvol factor not included),'
891 10 : write(ab_out,'(a)')' efidir atom atdir qgrdir real part imaginary part'
892 10 : i1pert=natom+2
893 40 : do i3dir=1,3
894 120 : do i1dir=1,3
895 327 : do i2pert=1,natom
896 918 : do i2dir=1,3
897 828 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
898 621 : write(ab_out,'(4(i5,3x),2(1x,f20.10))') i1dir,i2pert,i2dir,i3dir, &
899 621 : & d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert), &
900 1242 : & -d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
901 : end if
902 : end do
903 : end do
904 : end do
905 40 : write(ab_out,*)' '
906 : end do
907 :
908 : !Calculate cuadrupoles (symmetrize i1dir/i3dir)
909 70 : ABI_MALLOC(qdrp,(2,3,mpert,3,mpert,3,mpert))
910 : i1pert=natom+2
911 33 : do i2pert=1,natom
912 102 : do i2dir=1,3
913 299 : do i1dir=1,3
914 414 : do i3dir=1,i1dir-1
915 414 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
916 : !real part
917 : qdrp(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=&
918 : & d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + &
919 207 : & d3etot_car(2,i3dir,i1pert,i2dir,i2pert,i1dir,i3pert)
920 :
921 : qdrp(1,i3dir,i1pert,i2dir,i2pert,i1dir,i3pert)=&
922 207 : & qdrp(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
923 :
924 : !imaginary part
925 : qdrp(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=&
926 : & -(d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + &
927 207 : & d3etot_car(1,i3dir,i1pert,i2dir,i2pert,i1dir,i3pert) )
928 :
929 : qdrp(2,i3dir,i1pert,i2dir,i2pert,i1dir,i3pert)=&
930 207 : & qdrp(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
931 : end if
932 : end do
933 276 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i1dir,i3pert)==1) then
934 : !real part
935 : qdrp(1,i1dir,i1pert,i2dir,i2pert,i1dir,i3pert)=&
936 207 : & two*d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i1dir,i3pert)
937 :
938 : !imaginary part
939 : qdrp(2,i1dir,i1pert,i2dir,i2pert,i1dir,i3pert)=&
940 207 : &-two*d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i1dir,i3pert)
941 : end if
942 : end do
943 : end do
944 : end do
945 :
946 10 : write(ab_out,'(a)')' Quadrupole tensor, in cartesian coordinates,'
947 10 : write(ab_out,'(a)')' efidir atom atdir qgrdir real part imaginary part'
948 40 : do i3dir=1,3
949 120 : do i1dir=1,3
950 327 : do i2pert=1,natom
951 918 : do i2dir=1,3
952 828 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
953 621 : write(ab_out,'(4(i5,3x),2(1x,f20.10))') i1dir,i2pert,i2dir,i3dir, &
954 621 : & qdrp(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert), &
955 1242 : & qdrp(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
956 : end if
957 : end do
958 : end do
959 : end do
960 40 : write(ab_out,*)' '
961 : end do
962 :
963 10 : write(ab_out,'(a)')' Electronic (clamped-ion) contribution to the piezoelectric tensor,'
964 10 : write(ab_out,'(a)')' in cartesian coordinates, (from sum rule of P^1 tensor)'
965 10 : write(ab_out,'(a)')' efidir atdir qgrdir real part imaginary part'
966 40 : do i3dir=1,3
967 120 : do i1dir=1,3
968 390 : do i2dir=1,3
969 270 : piezoci=zero
970 891 : do i2pert=1,natom
971 891 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
972 621 : piezoci(1)=piezoci(1)+d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
973 621 : piezoci(2)=piezoci(2)-d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
974 : end if
975 : end do
976 270 : piezoci(1)=-piezoci(1)/ucvol
977 270 : piezoci(2)=-piezoci(2)/ucvol
978 360 : write(ab_out,'(3(i5,3x),2(1x,f20.10))') i1dir,i2dir,i3dir,piezoci(1),piezoci(2)
979 : end do
980 : end do
981 40 : write(ab_out,*)' '
982 : end do
983 :
984 10 : write(ab_out,'(a)')' Electronic (clamped-ion) contribution to the piezoelectric tensor,'
985 10 : write(ab_out,'(a)')' in cartesian coordinates, (from sum rule of dynamic quadrupoles)'
986 10 : write(ab_out,'(a)')' efidir atdir qgrdir real part imaginary part'
987 40 : do i3dir=1,3
988 120 : do i1dir=1,3
989 390 : do i2dir=1,3
990 270 : piezoci=zero
991 891 : do i2pert=1,natom
992 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1.and. &
993 621 : blkflg_car(i1dir,i1pert,i3dir,i2pert,i2dir,i3pert)==1.and. &
994 270 : blkflg_car(i2dir,i1pert,i1dir,i2pert,i3dir,i3pert)==1) then
995 : piezoci(1)=piezoci(1)+qdrp(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) &
996 : & +qdrp(1,i1dir,i1pert,i3dir,i2pert,i2dir,i3pert) &
997 621 : & -qdrp(1,i2dir,i1pert,i1dir,i2pert,i3dir,i3pert)
998 : piezoci(2)=piezoci(2)+qdrp(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) &
999 : & +qdrp(2,i1dir,i1pert,i3dir,i2pert,i2dir,i3pert) &
1000 621 : & -qdrp(2,i2dir,i1pert,i1dir,i2pert,i3dir,i3pert)
1001 : end if
1002 : end do
1003 270 : piezoci(1)=-piezoci(1)/(two*ucvol)
1004 270 : piezoci(2)=-piezoci(2)/(two*ucvol)
1005 360 : write(ab_out,'(3(i5,3x),2(1x,f20.10))') i1dir,i2dir,i3dir,piezoci(1),piezoci(2)
1006 : end do
1007 : end do
1008 40 : write(ab_out,*)' '
1009 : end do
1010 10 : ABI_FREE(qdrp)
1011 : end if
1012 :
1013 17 : if (lw_flexo==2.or.lw_flexo==1) then
1014 3 : write(ab_out,'(a)')' Clamped-ion flexoelectric tensor (type-II), in cartesian coordinates,'
1015 3 : write(ab_out,'(a)')' efidir qgrdir strdir1 strdir2 real part imaginary part'
1016 3 : i1pert=natom+2
1017 12 : do i3dir=1,3
1018 30 : do i2pert=natom+3,natom+4
1019 72 : do i2dir=1,3
1020 54 : istr=(i2pert-natom-3)*3+i2dir
1021 54 : beta=idx(2*istr-1); delta=idx(2*istr)
1022 234 : do i1dir=1,3
1023 216 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
1024 162 : write(ab_out,'(4(i5,3x),2(1x,f20.10))') i1dir,i3dir,beta,delta, &
1025 162 : & d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/ucvol, &
1026 324 : & d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/ucvol
1027 : end if
1028 : end do
1029 : end do
1030 27 : write(ab_out,*)' '
1031 : end do
1032 : end do
1033 : end if
1034 :
1035 17 : if (lw_flexo==3.or.lw_flexo==1) then
1036 3 : write(ab_out,'(a)')' 1st real-space moment of IFCs, in cartesian coordinates,'
1037 3 : write(ab_out,'(a)')' iatdir iatom jatdir jatom qgrdir real part imaginary part'
1038 12 : do i3dir=1,3
1039 27 : do i1pert=1,natom
1040 81 : do i1dir=1,3
1041 180 : do i2pert=1,natom
1042 486 : do i2dir=1,3
1043 432 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
1044 324 : write(ab_out,'(5(i5,4x),2(1x,f20.10))') i1dir,i1pert,i2dir,i2pert,i3dir, &
1045 324 : & -d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),&
1046 648 : & d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1047 : end if
1048 : end do
1049 : end do
1050 : end do
1051 : end do
1052 12 : write(ab_out,*)' '
1053 : end do
1054 :
1055 3 : write(ab_out,'(a)')' Piezoelectric force-response tensor, in cartesian coordinates '
1056 3 : write(ab_out,'(a)')' (from sum rule of 1st moment of IFCs),'
1057 3 : write(ab_out,'(a)')' (for non-vanishing forces in the cell it lacks an improper contribution),'
1058 3 : write(ab_out,'(a)')' iatom iatddir jatddir qgrdir real part imaginary part'
1059 12 : do i3dir=1,3
1060 27 : do i1pert=1,natom
1061 81 : do i1dir=1,3
1062 234 : do i2dir=1,3
1063 162 : piezofr=zero
1064 486 : do i2pert=1,natom
1065 486 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
1066 324 : piezofr(1)=piezofr(1)-d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1067 324 : piezofr(2)=piezofr(2)+d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1068 : end if
1069 : end do
1070 162 : write(ab_out,'(4(i5,4x),2(1x,f20.10))') i1pert,i1dir,i2dir,i3dir, &
1071 378 : & piezofr(1), piezofr(2)
1072 : end do
1073 : end do
1074 : end do
1075 12 : write(ab_out,*)' '
1076 : end do
1077 : end if
1078 :
1079 17 : if (lw_flexo==4.or.lw_flexo==1) then
1080 3 : write(ab_out,'(a)')' Clamped-ion flexoelectric force-response tensor (type-II), in cartesian coordinates,'
1081 3 : write(ab_out,'(a)')' atom atdir qgrdir strdir1 strdir2 real part imaginary part'
1082 12 : do i3dir=1,3
1083 30 : do i1pert=1,natom
1084 81 : do i1dir=1,3
1085 162 : do i2pert=natom+3, natom+4
1086 486 : do i2dir=1,3
1087 324 : istr=(i2pert-natom-3)*3+i2dir
1088 324 : beta=idx(2*istr-1); delta=idx(2*istr)
1089 432 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
1090 324 : write(ab_out,'(5(i5,3x),2(1x,f20.10))') i1pert,i1dir,i3dir,beta,delta, &
1091 324 : & d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),&
1092 648 : & d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1093 : end if
1094 : end do
1095 : end do
1096 72 : write(ab_out,*)' '
1097 : end do
1098 : end do
1099 : end do
1100 :
1101 :
1102 3 : write(ab_out,'(a)')' Clamped-ion elastic tensor, in cartesian coordinates '
1103 3 : write(ab_out,'(a)')' (from sum rule of clamped-ion flexoelectric force-response tensor),'
1104 3 : write(ab_out,'(a)')' (for stressed cells it lacks an improper contribution),'
1105 3 : write(ab_out,'(a)')' atdir qgrdir strdir1 strdir2 real part imaginary part'
1106 12 : do i1dir=1,3
1107 30 : do i3dir=1,i1dir
1108 54 : do i2pert=natom+3, natom+4
1109 162 : do i2dir=1,3
1110 108 : istr=(i2pert-natom-3)*3+i2dir
1111 108 : beta=idx(2*istr-1); delta=idx(2*istr)
1112 108 : celastci=zero
1113 324 : do i1pert=1,natom
1114 324 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
1115 216 : celastci(1)=celastci(1)+d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1116 216 : celastci(2)=celastci(2)+d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1117 : end if
1118 : end do
1119 108 : write(ab_out,'(4(i5,3x),2(1x,f20.10))') i1dir,i3dir,beta,delta, &
1120 252 : & celastci(1)/ucvol,celastci(2)/ucvol
1121 : end do
1122 : end do
1123 27 : write(ab_out,*)' '
1124 : end do
1125 : end do
1126 : end if
1127 :
1128 17 : if (lw_natopt==1) then
1129 1 : write(ab_out,'(a)')' Natural optical activity tensor, in cartesian coordinates,'
1130 1 : write(ab_out,'(a)')' (1/ucvol factor not included),'
1131 1 : write(ab_out,'(a)')' efidir1 efidir2 qgrdir real part imaginary part'
1132 1 : i1pert=natom+2
1133 1 : i2pert=natom+2
1134 4 : do i3dir=1,3
1135 12 : do i1dir=1,3
1136 39 : do i2dir=1,3
1137 36 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
1138 27 : write(ab_out,'(3(i5,3x),2(1x,f20.10))') i1dir,i2dir,i3dir, &
1139 27 : & -four*pi*d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert), &
1140 54 : & four*pi*d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1141 : end if
1142 : end do
1143 : end do
1144 4 : write(ab_out,*)' '
1145 : end do
1146 : end if
1147 :
1148 17 : if (timdisp==1) then
1149 6 : i3pert=natom+9
1150 6 : i3dir=1
1151 6 : write(ab_out,'(a)')' Frequency derivative of second-order energies, in cartesian coordinates'
1152 6 : write(ab_out,'(a)')' j1 j2 matrix elementi'
1153 6 : write(ab_out,'(a)')' dir pert dir pert real part imaginary part'
1154 128 : do i1pert=1,mpert
1155 494 : do i1dir=1,3
1156 8186 : do i2pert=1,mpert
1157 31158 : do i2dir=1,3
1158 30792 : if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
1159 30 : write(ab_out,'(2(i4,i5),2(1x,f20.10))')i1dir,i1pert,i2dir,i2pert,&
1160 30 : & d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert), &
1161 60 : & d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1162 : end if
1163 : end do
1164 : end do
1165 : end do
1166 : end do
1167 : end if
1168 :
1169 : DBG_EXIT("COLL")
1170 :
1171 17 : end subroutine dfptlw_out
1172 : !!***
1173 :
1174 : end module m_longwave
1175 : !!***
|
