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