Line data Source code
1 : !!****m* ABINIT/m_respfn_driver
2 : !! NAME
3 : !! m_respfn_driver
4 : !!
5 : !! FUNCTION
6 : !! Subdriver for DFPT calculations.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (XG, DRH, MT, MKV, GA)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : ! nvtx related macro definition
23 : #include "nvtx_macros.h"
24 :
25 : module m_respfn_driver
26 :
27 : use defs_basis
28 : use defs_wvltypes
29 : use m_efmas_defs
30 : use m_abicore
31 : use m_xmpi
32 : use m_exit
33 : use m_wffile
34 : use m_errors
35 : use m_ebands
36 : use m_results_respfn
37 : use m_hdr
38 : use m_crystal
39 : use m_xcdata
40 : use m_dtset
41 : use m_dtfil
42 : use m_gemm_nonlop_projectors
43 :
44 : use defs_datatypes, only : pseudopotential_type
45 : use defs_abitypes, only : MPI_type
46 : use m_time, only : timab
47 : use m_fstrings, only : strcat, endswith
48 : use m_matrix, only : matr3inv
49 : use m_symtk, only : littlegroup_q, symmetrize_xred
50 : use m_fft, only : zerosym, fourdp
51 : use m_kpts, only : symkchk
52 : use m_geometry, only : irreducible_set_pert, symredcart, sylwtens
53 : use m_dynmat, only : chkph3, d2sym3, q0dy3_apply, q0dy3_calc, wings3, dfpt_phfrq, sytens, dfpt_prtph, &
54 : asria_calc, asria_corr, cart29, cart39, chneu9, dfpt_sydy
55 : use m_ddb, only : ddb_type
56 : use m_ddb_hdr, only : ddb_hdr_type
57 : use m_occ, only : newocc
58 : use m_efmas, only : efmasdeg_free_array, efmasval_free_array
59 : use m_wfk, only : wfk_read_eigenvalues, wfk_read_my_kptbands
60 : use m_ioarr, only : read_rhor
61 : use m_pawang, only : pawang_type
62 : use m_pawrad, only : pawrad_type
63 : use m_pawtab, only : pawtab_type, pawtab_get_lsize
64 : use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
65 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
66 : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
67 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, &
68 : pawrhoij_bcast, pawrhoij_nullify, pawrhoij_inquire_dim, &
69 : pawrhoij_print_rhoij, pawrhoij_io
70 : use m_pawdij, only : pawdij, symdij, pawdij_print_dij
71 : use m_paw_energies,only : paw_energies_type
72 : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
73 : use m_paw_finegrid,only : pawexpiqr
74 : use m_pawxc, only : pawxc_get_nkxc, pawxc_get_usekden
75 : use m_paw_dmft, only : paw_dmft_type
76 : use m_paw_sphharm, only : setsym_ylm
77 : use m_paw_nhat, only : nhatgrid,pawmknhat
78 : use m_paw_tools, only : chkpawovlp
79 : use m_paw_denpot, only : pawdenpot
80 : use m_paw_init, only : pawinit,paw_gencond
81 : use m_kg, only : getcut, getph, kpgio
82 : use m_eig2d, only : eig2tot, elph2_fanddw
83 : use m_inwffil, only : inwffil
84 : use m_spacepar, only : hartre, setsym
85 : use m_mkrho, only : mkrho
86 : use m_vdw_dftd2, only : vdw_dftd2
87 : use m_vdw_dftd3, only : vdw_dftd3
88 : use m_initylmg, only : initylmg
89 : use m_pspini, only : pspini
90 : use m_atm2fft, only : atm2fft
91 : use m_dfpt_loopert,only : dfpt_looppert, eigen_meandege
92 : use m_rhotoxc, only : rhotoxc
93 : use m_drivexc, only : check_kxc, xc_need_kden
94 : use m_xc_tb09, only : xc_tb09_update_c
95 : use m_mklocl, only : mklocl, mklocl_recipspace
96 : use m_common, only : setup1, prteigrs
97 : use m_fourier_interpol, only : transgrid
98 : use m_paral_atom, only : get_my_atmtab, free_my_atmtab
99 : use m_paw_occupancies, only : initrhoij
100 : use m_paw_correlations,only : pawpuxinit
101 : use m_mkcore, only : mkcore, dfpt_mkcore
102 : use m_dfpt_elt, only : dfpt_eltfrxc, dfpt_eltfrloc, dfpt_eltfrkin, dfpt_eltfrhar, elt_ewald, dfpt_ewald
103 : use m_d2frnl, only : d2frnl
104 : use m_pstat, only : pstat_proc
105 :
106 : #if defined HAVE_GPU
107 : use m_alloc_hamilt_gpu
108 : #endif
109 :
110 : #if defined(HAVE_GPU_MARKERS)
111 : use m_nvtx_data
112 : #endif
113 :
114 : implicit none
115 :
116 : private
117 : !!***
118 :
119 : public :: respfn
120 : !!***
121 :
122 : contains
123 : !!***
124 :
125 : !!****f* m_respfn_driver/respfn
126 : !! NAME
127 : !! respfn
128 : !!
129 : !! FUNCTION
130 : !! Primary routine for conducting DFT calculations of Response functions.
131 : !!
132 : !! INPUTS
133 : !! codvsn=code version
134 : !! cpui=initial cpu time
135 : !! dtfil <type(datafiles_type)>=variables related to files
136 : !! dtset <type(dataset_type)>=all input variables for this dataset
137 : !! | mband=maximum number of bands
138 : !! | mgfft=maximum single fft dimension
139 : !! | mkmem=Number of k points treated by this node
140 : !! | mpw=maximum number of planewaves in basis sphere (large number)
141 : !! | natom=number of atoms in unit cell
142 : !! | nfft=(effective) number of FFT grid points (for this processor)
143 : !! | nkpt=number of k points
144 : !! | nspden=number of spin-density components
145 : !! | nsppol=number of channels for spin-polarization (1 or 2)
146 : !! | nsym=number of symmetry elements in space group
147 : !! mkmems(3)=array containing the tree values of mkmem (see above) (k-GS, k+q-GS and RF)
148 : !! mpi_enreg=information about MPI parallelization
149 : !! npwtot(nkpt)=number of planewaves in basis and boundary at each k point
150 : !! xred(3,natom)=reduced dimensionless atomic coordinates
151 : !!
152 : !! OUTPUT
153 : !! etotal=total energy (sum of 7 or 8 contributions) (hartree)
154 : !!
155 : !! SIDE EFFECTS
156 : !! iexit=index of "exit" on first line of file (0 if not found)
157 : !! occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point
158 : !! Occupations number may have been read from a previous dataset...
159 : !! pawang <type(pawang_type)>=paw angular mesh and related data
160 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
161 : !! Some dimensions in pawrad have been set in driver.f
162 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
163 : !! Some dimensions in pawtab have been set in driver.f
164 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
165 : !! Before entering the first time in respfn, a significant part of psps
166 : !! has been initialized: the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,
167 : !! mpsso,mgrid,ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp
168 : !! All the remaining components of psps are to be initialized in the call
169 : !! to pspini. The next time the code enters respfn, psps might be identical
170 : !! to the one of the previous dtset, in which case, no reinitialisation
171 : !! is scheduled in pspini.f .
172 : !! results_respfn <type(results_respfn_type)>=stores some results of respfn calls
173 : !!
174 : !! NOTES
175 : !! USE OF FFT GRIDS:
176 : !! =================
177 : !! In case of PAW:
178 : !! ---------------
179 : !! Two FFT grids are used:
180 : !! - A "coarse" FFT grid (defined by ecut)
181 : !! for the application of the Hamiltonian on the plane waves basis.
182 : !! It is defined by nfft, ngfft, mgfft, ...
183 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ...
184 : !! are expressed on this grid.
185 : !! - A "fine" FFT grid (defined) by ecutdg)
186 : !! for the computation of the density inside PAW spheres.
187 : !! It is defined by nfftf, ngfftf, mgfftf, ...
188 : !! Total density, potentials, ...
189 : !! are expressed on this grid.
190 : !! In case of norm-conserving:
191 : !! ---------------------------
192 : !! - Only the usual FFT grid (defined by ecut) is used.
193 : !! It is defined by nfft, ngfft, mgfft, ...
194 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf)
195 : !! are set equal to (nfft,ngfft,mgfft) in that case.
196 : !!
197 : !! SOURCE
198 :
199 1003 : subroutine respfn(codvsn,cpui,dtfil,dtset,etotal,iexit,&
200 716 : & mkmems,mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,psps,results_respfn,xred)
201 :
202 : !Arguments ------------------------------------
203 : integer,intent(inout) :: iexit
204 : real(dp),intent(in) :: cpui
205 : real(dp),intent(inout) :: etotal !vz_i
206 : character(len=8),intent(in) :: codvsn
207 : type(MPI_type),intent(inout) :: mpi_enreg
208 : type(datafiles_type),intent(in) :: dtfil
209 : type(dataset_type),intent(in) :: dtset
210 : type(pawang_type),intent(inout) :: pawang
211 : type(pseudopotential_type),intent(inout) :: psps
212 : integer,intent(in) :: mkmems(3)
213 : integer,intent(inout) :: npwtot(dtset%nkpt)
214 : real(dp),intent(inout) :: xred(3,dtset%natom)
215 : real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
216 : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
217 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
218 : type(results_respfn_type),intent(inout) :: results_respfn
219 :
220 : !Local variables-------------------------------
221 : integer,parameter :: formeig=0,level=10
222 : integer,parameter :: response=1,syuse=0,master=0,cplex1=1
223 : integer :: nk3xc
224 : integer :: analyt,ask_accurate,asr,bantot,bdeigrf,chneut,coredens_method,coretau_method,cplex,cplex_rhoij
225 : !integer :: nkpt_eff, band_index, ikpt, isppol, nkpt_max, nband_k,
226 : integer :: dim_eig2nkq,dim_eigbrd,dyfr_cplex,dyfr_nondiag,gnt_option
227 : integer :: gscase,has_dijnd,has_diju,has_vhartree,has_kxc,iatom,iatom_tot,iband,idir,ider,ierr,ifft,ii,indx
228 : integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert
229 : integer :: initialized,ipert,ipert2,ireadwf0,iscf,iscf_eff,ispden,isym
230 : integer :: itypat,izero,mcg,me,mgfftf,mk1mem,mkqmem,mpert,mu
231 : integer :: my_natom,n1,natom,n3xccc,nfftf,nfftot,nfftotf,nhatdim,nhatgrdim
232 : integer :: nkpt_rbz,nkxc,nkxc1,nspden_rhoij,ntypat,nzlmopt,openexit
233 : integer :: optcut,option,optgr0,optgr1,optgr2,optorth,optrad
234 : integer :: optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv
235 : integer :: outd2,pawbec,pawpiezo,prtbbb,psp_gencond,qzero,rdwrpaw
236 : integer :: rfddk,rfelfd,rfphon,rfstrs,rf2_dkdk,rf2_dkde,rfmagn
237 : integer :: spaceworld,sumg0,sumg0_save,sz1,sz2,tim_mkrho,timrev,usecprj,usevdw,usevxctau,usevxctau_paw
238 : integer :: usexcnhat,use_sym,vloc_method,zero_by_symm
239 : logical :: has_full_piezo,has_allddk,is_dfpt=.true.,non_magnetic_xc
240 : logical :: paral_atom,qeq0,use_nhat_gga,call_pawinit
241 : real(dp) :: boxcut,compch_fft,compch_sph,cpus,ecore,ecut_eff,ecutdg_eff,ecutf
242 : real(dp) :: eei,eew,ehart,eii,ek,enl,entropy,bigexc,bigsxc
243 : real(dp) :: etot,evdw,fermie,fermih,gsqcut,gsqcut_eff,gsqcutc_eff,qphnrm,residm
244 : real(dp) :: ucvol,vxcavg,el_temp
245 : character(len=500) :: msg
246 716 : type(ebands_t) :: bstruct
247 716 : type(hdr_type) :: hdr,hdr_fine,hdr0,hdr_den
248 716 : type(ddb_type) :: ddb
249 37232 : type(ddb_hdr_type) :: ddb_hdr
250 716 : type(paw_dmft_type) :: paw_dmft
251 : type(paw_energies_type) :: paw_energies
252 : type(pawfgr_type) :: pawfgr
253 716 : type(wvl_data) :: wvl
254 : type(xcdata_type) :: xcdata
255 : integer :: ddkfil(3),ngfft(18),ngfftf(18),rfdir(3),rf2_dirs_from_rfpert_nl(3,3)
256 716 : integer,allocatable :: atindx(:),atindx1(:),blkflg(:,:,:,:),blkflgfrx1(:,:,:,:),blkflg1(:,:,:,:)
257 1432 : integer,allocatable :: blkflg2(:,:,:,:),carflg(:,:,:,:),clflg(:,:),indsym(:,:,:)
258 716 : integer,allocatable :: irrzon(:,:,:),kg(:,:),l_size_atm(:),nattyp(:),npwarr(:)
259 716 : integer,allocatable :: pertsy(:,:),rfpert(:)
260 716 : integer,allocatable :: rfpert_lw(:,:,:,:,:,:),rfpert_nl(:,:,:,:,:,:),symq(:,:,:),symrec(:,:,:)
261 716 : logical,allocatable :: distrb_flags(:,:,:)
262 : real(dp) :: dum_gauss(0),dum_dyfrn(0),dum_dyfrv(0),dum_eltfrxc(0)
263 : real(dp) :: dum_grn(0),dum_grv(0),dum_rhog(0),dum_vg(0)
264 : real(dp) :: dummy6(6),gmet(3,3),gmet_for_kg(3,3),gprimd(3,3),gprimd_for_kg(3,3),qphon(3)
265 : real(dp) :: dummy_in(0)
266 : real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0)
267 : real(dp) :: rmet(3,3),rprimd(3,3),rprimd_for_kg(3,3),strn_dummy6(6),strv_dummy6(6),tsec(2)
268 : real(dp),parameter :: k0(3)=(/zero,zero,zero/)
269 716 : real(dp),allocatable :: becfrnl(:,:,:),cg(:,:),d2bbb(:,:,:,:,:,:),d2cart(:,:,:,:,:)
270 716 : real(dp),allocatable :: d2cart_bbb(:,:,:,:,:,:),d2eig0(:,:,:,:,:)
271 716 : real(dp),allocatable :: d2k0(:,:,:,:,:),d2lo(:,:,:,:,:),d2loc0(:,:,:,:,:)
272 716 : real(dp),allocatable :: d2matr(:,:,:,:,:),d2nfr(:,:,:,:,:),d2nl(:,:,:,:,:),d2ovl(:,:,:,:,:)
273 716 : real(dp),allocatable :: d2nl0(:,:,:,:,:),d2nl1(:,:,:,:,:),d2tmp(:,:,:,:,:),d2vn(:,:,:,:,:)
274 716 : real(dp),allocatable :: displ(:),doccde(:)
275 716 : real(dp),allocatable :: dyew(:,:,:,:,:),dyewq0(:,:,:),dyfrlo(:,:,:),dyfrlo_indx(:,:,:)
276 716 : real(dp),allocatable :: dyfrnl(:,:,:,:,:),dyfrwf(:,:,:,:,:),dyfrx1(:,:,:,:,:),dyvdw(:,:,:,:,:)
277 716 : real(dp),allocatable :: dyfrx2(:,:,:),eigen0(:),eigval(:),eigvec(:)
278 716 : real(dp),allocatable :: eig2nkq(:,:,:,:,:,:,:),eigbrd(:,:,:,:,:,:,:)
279 716 : real(dp),allocatable :: eigen_fan(:),eigen_ddw(:),eigen_fanddw(:)
280 716 : real(dp),allocatable :: eigen_fan_mean(:),eigen_ddw_mean(:)
281 716 : real(dp),allocatable :: eltcore(:,:),elteew(:,:),eltfrhar(:,:),eltfrkin(:,:)
282 716 : real(dp),allocatable :: eltfrloc(:,:),eltfrnl(:,:),eltfrxc(:,:),eltvdw(:,:),grtn_indx(:,:)
283 716 : real(dp),allocatable :: grxc(:,:),kxc(:,:),nhat(:,:),nhatgr(:,:,:)
284 716 : real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phfrq(:),phnons(:,:,:),piezofrnl(:,:)
285 716 : real(dp),allocatable :: rhog(:,:),rhor(:,:),rhowfg(:,:),rhowfr(:,:)
286 716 : real(dp),allocatable :: symrel_cart(:,:,:),taug(:,:),taur(:,:)
287 716 : real(dp),allocatable :: vhartr(:),vpsp(:),vtrial(:,:)
288 716 : real(dp),allocatable :: vxc(:,:),vxctau(:,:,:),xccc3d(:),xcctau3d(:),ylm(:,:),ylmgr(:,:,:)
289 716 : real(dp),pointer :: eigenq_fine(:,:,:),eigen1_pert(:,:,:)
290 716 : real(dp),allocatable :: eigen0_pert(:),eigenq_pert(:),occ_rbz_pert(:)
291 716 : type(efmasdeg_type),allocatable :: efmasdeg(:)
292 716 : type(efmasval_type),allocatable :: efmasval(:,:)
293 716 : type(paw_an_type),allocatable :: paw_an(:)
294 716 : type(paw_ij_type),allocatable :: paw_ij(:)
295 716 : type(pawfgrtab_type),allocatable,save :: pawfgrtab(:)
296 716 : type(pawrhoij_type),allocatable :: pawrhoij(:),pawrhoij_read(:)
297 :
298 : ! ***********************************************************************
299 :
300 : DBG_ENTER("COLL")
301 :
302 716 : call timab(132,1,tsec)
303 716 : call timab(133,1,tsec)
304 : ABI_NVTX_START_RANGE(NVTX_RESPFN)
305 :
306 : !Some data for parallelism
307 :
308 716 : my_natom=mpi_enreg%my_natom
309 716 : paral_atom=(my_natom/=dtset%natom)
310 : !Define FFT grid(s) sizes (be careful !)
311 : !See NOTES in the comments at the beginning of this file.
312 716 : call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
313 :
314 : !Structured debugging if dtset%prtvol==-level
315 716 : if(dtset%prtvol==-level)then
316 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' respfn : enter , debug mode '
317 0 : call wrtout(std_out,msg)
318 : end if
319 :
320 : !Option input variables
321 716 : iscf=dtset%iscf
322 :
323 : !Get electronic temperature from dtset
324 716 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
325 :
326 : !Respfn input variables
327 2864 : asr=dtset%asr ; chneut=dtset%chneut ; rfdir(1:3)=dtset%rfdir(1:3)
328 716 : rfddk=dtset%rfddk ; rfelfd=dtset%rfelfd ; rfmagn=dtset%rfmagn
329 716 : rfphon=dtset%rfphon ; rfstrs=dtset%rfstrs
330 716 : rf2_dkdk=dtset%rf2_dkdk ; rf2_dkde=dtset%rf2_dkde
331 :
332 : !Scalar potential has no idir.
333 716 : if (rfmagn==3) then
334 1 : rfdir(:)= 0
335 1 : rfdir(1)= 1
336 : end if
337 :
338 716 : pawbec=0 ; if(psps%usepaw==1.and.(rfphon==1.or.(rfelfd==1.or.rfelfd==3))) pawbec=1
339 716 : pawpiezo=0; if(psps%usepaw==1.and.(rfstrs/=0.or.(rfelfd==1.or.rfelfd==3))) pawpiezo=1
340 : !AM 10152015 -- WARNING --- the full calculation of the piezoelectric tensor
341 : !from electric field perturbation is only available
342 : !if nsym/=1 (strain perturbation is not symmetrized):
343 716 : has_full_piezo=.False. ; if(pawpiezo==1.and.dtset%nsym==1) has_full_piezo=.True.
344 716 : usevdw=0;if (dtset%vdw_xc>=5.and.dtset%vdw_xc<=7) usevdw=1
345 : !mkmem variables (mkmem is already argument)
346 716 : mkqmem=mkmems(2) ; mk1mem=mkmems(3)
347 :
348 716 : ntypat=psps%ntypat
349 716 : natom=dtset%natom
350 2864 : nfftot=product(ngfft(1:3))
351 2864 : nfftotf=product(ngfftf(1:3))
352 :
353 : !LIKELY TO BE TAKEN AWAY
354 716 : initialized=0
355 716 : ek=zero ; ehart=zero ; bigexc=zero ; eei=zero ; enl=zero
356 716 : eii=zero ; eew=zero ; ecore=zero ; bigsxc=zero
357 :
358 : !Set up for iterations
359 : call setup1(dtset%acell_orig(1:3,1),bantot,dtset,&
360 : ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,&
361 : ngfftf,ngfft,dtset%nkpt,dtset%nsppol,&
362 716 : response,rmet,dtset%rprim_orig(1:3,1:3,1),rprimd,ucvol,psps%usepaw)
363 : !In some cases (e.g. getcell/=0), the plane wave vectors have
364 : ! to be generated from the original simulation cell
365 716 : rprimd_for_kg=rprimd
366 752 : if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=dtset%rprimd_orig(:,:,1)
367 716 : call matr3inv(rprimd_for_kg,gprimd_for_kg)
368 28640 : gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg)
369 :
370 : !Define the set of admitted perturbations
371 : ! Note that we have a global parameter (mpert=natom+MPERT_MAX)
372 : ! with MPERT_MAX=8, but we use a smaller value here.
373 716 : mpert=natom+7
374 716 : if (rf2_dkdk>0.or.rf2_dkde>0) mpert=natom+11
375 716 : if (dtset%rfmagn==2.or.dtset%prt1mag/=0) mpert=2*dtset%natom+11
376 :
377 : !Initialize the list of perturbations rfpert
378 2148 : ABI_MALLOC(rfpert,(mpert))
379 7425 : rfpert(:)=0
380 1486 : if(rfphon==1)rfpert(dtset%rfatpol(1):dtset%rfatpol(2))=1
381 :
382 716 : if(rfddk==1)rfpert(natom+1)=1
383 :
384 716 : if(rf2_dkdk/=0)rfpert(natom+10)=1
385 716 : if(rf2_dkde/=0)rfpert(natom+11)=1
386 :
387 716 : if(rfelfd==1.or.rfelfd==2)rfpert(natom+1)=1
388 716 : if(rfelfd==1.or.rfelfd==3)rfpert(natom+2)=1
389 :
390 716 : if(rfstrs==1.or.rfstrs==3)rfpert(natom+3)=1
391 716 : if(rfstrs==2.or.rfstrs==3)rfpert(natom+4)=1
392 :
393 716 : if(rfmagn==1) rfpert(natom+5)=1
394 719 : if(rfmagn==2) rfpert(natom+11+dtset%rfatpol(1):natom+11+dtset%rfatpol(2))=1
395 716 : if(rfmagn==3) rfpert(natom+6)=1
396 :
397 716 : qeq0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2<1.d-14)
398 :
399 : !Init spaceworld
400 716 : spaceworld=mpi_enreg%comm_cell
401 716 : me = xmpi_comm_rank(spaceworld)
402 :
403 : !Set up the basis sphere of planewaves
404 2148 : ABI_MALLOC(kg,(3,dtset%mpw*dtset%mkmem))
405 2148 : ABI_MALLOC(npwarr,(dtset%nkpt))
406 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet_for_kg,dtset%istwfk,kg,&
407 : & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw,npwarr,npwtot,&
408 716 : & dtset%nsppol)
409 :
410 : !Set up the Ylm for each k point
411 2864 : ABI_MALLOC(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm))
412 716 : if (rfstrs/=0.or.pawbec==1.or.pawpiezo==1.or.dtset%efmas>0) then
413 372 : ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,9,psps%mpsang*psps%mpsang*psps%useylm))
414 : else
415 1246 : ABI_MALLOC(ylmgr,(0,0,psps%useylm))
416 : end if
417 716 : if (psps%useylm==1) then
418 172 : option=0
419 172 : if (rfstrs/=0.or.pawbec==1.or.pawpiezo==1.or.dtset%efmas>0) option=2
420 : call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,&
421 172 : & npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr)
422 : end if
423 :
424 716 : call timab(133,2,tsec)
425 716 : call timab(134,1,tsec)
426 :
427 : !Open and read pseudopotential files
428 : call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,&
429 716 : & psps,rprimd,comm_mpi=mpi_enreg%comm_cell)
430 :
431 716 : call timab(134,2,tsec)
432 716 : call timab(135,1,tsec)
433 :
434 : !Initialize band structure datatype
435 716 : call bstruct%from_dtset(dtset, npwarr)
436 :
437 : !Initialize PAW atomic occupancies
438 716 : if (psps%usepaw==1) then
439 524 : ABI_MALLOC(pawrhoij,(my_natom))
440 112 : call pawrhoij_nullify(pawrhoij)
441 : call initrhoij(dtset%pawcpxocc,dtset%lexexch,dtset%lpawu, &
442 : & my_natom,natom,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,&
443 : & pawrhoij,dtset%pawspnorb,pawtab,cplex1,dtset%spinat,dtset%typat,&
444 112 : & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
445 : else
446 604 : ABI_MALLOC(pawrhoij,(0))
447 : end if
448 :
449 : !Initialize header
450 716 : gscase=0
451 : call hdr%init(bstruct,codvsn,dtset,pawtab,gscase,psps,wvl%descr, &
452 716 : comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
453 :
454 : !Update header, with evolving variables, when available
455 : !Here, rprimd, xred and occ are available
456 716 : etot=hdr%etot ; fermie=hdr%fermie ; fermih=hdr%fermih ; residm=hdr%residm
457 :
458 : !If parallelism over atom, hdr is distributed
459 : call hdr%update(bantot,etot,fermie,fermih,&
460 : residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1), &
461 716 : comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
462 :
463 : !Clean band structure datatype (should use it more in the future !)
464 716 : call bstruct%free()
465 :
466 : !Initialize wavefunction files and wavefunctions.
467 716 : ireadwf0=1
468 :
469 :
470 716 : mcg=dtset%mpw*dtset%nspinor*dtset%mband_mem*dtset%mkmem*dtset%nsppol
471 2148 : ABI_MALLOC(cg,(2,mcg))
472 : !ABI_MALLOC_OR_DIE(cg,(2,mcg), ierr)
473 :
474 2148 : ABI_MALLOC(eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol))
475 339337 : eigen0(:)=zero ; ask_accurate=1
476 716 : optorth=0
477 :
478 : ! Initialize the wave function type and read GS WFK
479 3580 : ABI_MALLOC(distrb_flags,(dtset%nkpt,dtset%mband,dtset%nsppol))
480 347942 : distrb_flags = (mpi_enreg%proc_distrb == mpi_enreg%me_kpt)
481 : call wfk_read_my_kptbands(dtfil%fnamewffk, distrb_flags, spaceworld, dtset%ecut*(dtset%dilatmx)**2, &
482 : & formeig, dtset%istwfk, dtset%kptns, mcg, dtset%mband, dtset%mband_mem,dtset%mkmem,dtset%mpw,&
483 : & dtset%natom, dtset%nkpt, npwarr, dtset%nspinor, dtset%nsppol, dtset%usepaw,&
484 716 : & cg, eigen=eigen0, pawrhoij=hdr%pawrhoij)
485 716 : ABI_FREE(distrb_flags)
486 :
487 716 : if (psps%usepaw==1 .and. ireadwf0==1) then
488 : ! if parallelism, pawrhoij is distributed, hdr%pawrhoij is not
489 : call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,&
490 112 : & mpi_atmtab=mpi_enreg%my_atmtab)
491 : end if
492 :
493 716 : call timab(135,2,tsec)
494 716 : call timab(136,1,tsec)
495 :
496 : ! Report on eigen0 values ! Should use prteigrs.F90
497 : !write(msg, '(a,a)' )
498 : !call wrtout(std_out,ch10//' respfn : eigen0 array')
499 : !nkpt_eff=dtset%nkpt
500 : !nkpt_max=50;if(xmpi_paral==1)nkpt_max=-1
501 : !if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. dtset%nkpt>nkpt_max ) nkpt_eff=nkpt_max
502 : !band_index=0
503 : !do isppol=1,dtset%nsppol
504 : ! do ikpt=1,dtset%nkpt
505 : ! nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
506 : ! if(ikpt<=nkpt_eff)then
507 : ! write(msg, '(a,i2,a,i5)' )' isppol=',isppol,', k point number',ikpt
508 : ! call wrtout(std_out,msg)
509 : ! do iband=1,nband_k,4
510 : ! write(msg, '(a,4es16.6)')' ',eigen0(iband+band_index:min(iband+3,nband_k)+band_index)
511 : ! call wrtout(std_out,msg)
512 : ! end do
513 : ! else if(ikpt==nkpt_eff+1)then
514 : ! write(msg,'(a,a)' )' respfn : prtvol=0, 1 or 2, stop printing eigen0.',ch10
515 : ! call wrtout(std_out,msg)
516 : ! end if
517 : ! band_index=band_index+nband_k
518 : ! end do
519 : !end do
520 :
521 : !Allocation for forces and atomic positions (should be taken away, also argument ... )
522 2148 : ABI_MALLOC(grxc,(3,natom))
523 :
524 : !Do symmetry stuff
525 2864 : ABI_MALLOC(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
526 2864 : ABI_MALLOC(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
527 2864 : ABI_MALLOC(indsym,(4,dtset%nsym,natom))
528 2148 : ABI_MALLOC(symrec,(3,3,dtset%nsym))
529 17611705 : irrzon=0;indsym=0;symrec=0;phnons=zero
530 : !If the density is to be computed by mkrho, need irrzon and phnons
531 716 : iscf_eff=0;if(dtset%getden==0)iscf_eff=1
532 : call setsym(indsym,irrzon,iscf_eff,natom,&
533 : & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
534 716 : & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred)
535 :
536 : !Symmetrize atomic coordinates over space group elements:
537 716 : call symmetrize_xred(natom,dtset%nsym,dtset%symrel,dtset%tnons,xred,indsym=indsym)
538 :
539 : !Examine the symmetries of the q wavevector
540 2148 : ABI_MALLOC(symq,(4,2,dtset%nsym))
541 : timrev=1
542 :
543 : ! By default use symmetries.
544 716 : use_sym = 1
545 716 : if (dtset%prtgkk == 1)then
546 22 : use_sym = 0
547 22 : call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol,use_sym=use_sym)
548 : else
549 694 : call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol)
550 : end if
551 :
552 : !Deactivate time-reversal symmetry for finite-omega calculations
553 716 : if (abs(dtset%rfomega)>tol10.or.dtset%tim1rev==0) timrev=0
554 :
555 : !Generate an index table of atoms, in order for them to be used
556 : !type after type.
557 2148 : ABI_MALLOC(atindx,(natom))
558 1432 : ABI_MALLOC(atindx1,(natom))
559 2148 : ABI_MALLOC(nattyp,(ntypat))
560 1789 : indx=1
561 1789 : do itypat=1,ntypat
562 1073 : nattyp(itypat)=0
563 4280 : do iatom=1,natom
564 3564 : if(dtset%typat(iatom)==itypat)then
565 1527 : atindx(iatom)=indx
566 1527 : atindx1(indx)=iatom
567 1527 : indx=indx+1
568 1527 : nattyp(itypat)=nattyp(itypat)+1
569 : end if
570 : end do
571 : end do
572 :
573 : !Here allocation of GPU for fft calculations
574 : #if defined HAVE_GPU
575 : if (dtset%gpu_option/=ABI_GPU_DISABLED) then
576 : call alloc_hamilt_gpu(atindx1,dtset,gprimd,mpi_enreg,nattyp,npwarr,0,psps,dtset%gpu_option)
577 : end if
578 : #endif
579 :
580 : !Compute structure factor phases for current atomic pos:
581 2148 : ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*natom))
582 2148 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*natom))
583 716 : call getph(atindx,natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
584 :
585 716 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
586 94 : call getph(atindx,natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
587 : else
588 493732 : ph1df(:,:)=ph1d(:,:)
589 : end if
590 :
591 : !Compute occupation numbers and fermi energy, in case occupation scheme is metallic.
592 1432 : ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
593 716 : if( dtset%occopt>=3.and.dtset%occopt<=9 ) then
594 :
595 : call newocc(doccde,eigen0,entropy,fermie,fermih,dtset%ivalence,&
596 : & dtset%spinmagntarget,dtset%mband,dtset%nband,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,&
597 : & dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,dtset%occopt,dtset%prtvol,&
598 166 : & dtset%tphysel,dtset%tsmear,dtset%wtk)
599 :
600 : ! Update fermie and occ
601 166 : etot=hdr%etot ; residm=hdr%residm
602 166 : call hdr%update(bantot,etot,fermie,fermih,residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1))
603 :
604 : else
605 : ! doccde is irrelevant in this case
606 190267 : doccde(:)=zero
607 : end if
608 :
609 : !Recompute first large sphere cut-off gsqcut, without taking into account dilatmx
610 716 : ecutf=dtset%ecut
611 716 : if (psps%usepaw==1) then
612 112 : ecutf=dtset%pawecutdg
613 112 : call wrtout(std_out,ch10//' FFT (fine) grid used in SCF cycle:')
614 : end if
615 :
616 716 : call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf)
617 :
618 : !PAW: 1- Initialize values for several arrays depending only on atomic data
619 : !2- Check overlap
620 : !3- Identify FFT points in spheres and compute g_l(r).Y_lm(r) (and exp(-i.q.r) if needed)
621 : !4- Allocate PAW specific arrays
622 : !5- Compute perturbed local potential inside spheres
623 : !6- Eventually open temporary storage files
624 716 : if(psps%usepaw==1) then
625 : ! 1-Initialize values for several arrays depending only on atomic data
626 112 : gnt_option=1
627 : if (dtset%pawxcdev==2.or.dtset%rfphon/=0.or.dtset%rfstrs/=0.or.dtset%rfelfd==1.or.&
628 112 : dtset%rfelfd==3.or.dtset%rf2_dkde==1) gnt_option=2
629 :
630 : ! Test if we have to call pawinit
631 112 : call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
632 :
633 112 : if (psp_gencond==1.or.call_pawinit) then
634 : ! Some gen-cond have to be added...
635 56 : call timab(553,1,tsec)
636 : call pawinit(dtset%effmass_free,gnt_option,zero,zero,dtset%pawlcutd,dtset%pawlmix,&
637 : & psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,&
638 56 : & pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero)
639 : call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,&
640 56 : & rprimd,symrec,pawang%zarot)
641 :
642 : ! Update internal values
643 56 : call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
644 :
645 56 : call timab(553,2,tsec)
646 : else
647 56 : if (pawtab(1)%has_kij ==1) pawtab(1:psps%ntypat)%has_kij =2
648 56 : if (pawtab(1)%has_nabla==1) pawtab(1:psps%ntypat)%has_nabla=2
649 : end if
650 310 : psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore)
651 112 : call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot)
652 : call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,&
653 : & is_dfpt,dtset%jpawu,dtset%lexexch,dtset%lpawu,dtset%nspinor,ntypat,dtset%optdcmagpawu,pawang,dtset%pawprtvol,pawrad,&
654 112 : & pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu)
655 112 : compch_fft=-1.d5;compch_sph=-1.d5
656 310 : usexcnhat=maxval(pawtab(:)%usexcnhat)
657 112 : usecprj=dtset%pawusecp
658 : ! 2-Check overlap
659 112 : call chkpawovlp(natom,psps%ntypat,dtset%pawovlp,pawtab,rmet,dtset%typat,xred)
660 : ! 3-Identify FFT points in spheres and compute g_l(r).Y_lm(r) and exp(-i.q.r)
661 524 : ABI_MALLOC(pawfgrtab,(my_natom))
662 112 : if (my_natom>0) then
663 : call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,&
664 102 : & mpi_atmtab=mpi_enreg%my_atmtab)
665 : call pawfgrtab_init(pawfgrtab,1,l_size_atm,pawrhoij(1)%nspden,dtset%typat,&
666 102 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
667 102 : ABI_FREE(l_size_atm)
668 : end if
669 112 : use_nhat_gga=(dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0)
670 112 : optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm
671 112 : if (use_nhat_gga) then
672 4 : optgr1=dtset%pawstgylm
673 4 : if (rfphon==1) optgr2=1
674 : end if
675 112 : if (rfphon==1.or.rfstrs/=0.or.rfelfd==3.or.rf2_dkde==1) then ! LB2016-11-28 : Why not rfelfd==1?
676 55 : if (optgr1==0) optgr1=dtset%pawstgylm
677 55 : if (optgr2==0) optgr2=dtset%pawstgylm
678 55 : if (optrad==0.and.(.not.qeq0.or.rfstrs/=0.or.rfelfd==3.or.rf2_dkde==1)) optrad=1 ! LB2016-11-28 : Why not rfelfd==1?
679 : end if
680 112 : if (rfelfd==1.or.rfelfd==3.or.rf2_dkde==1) then
681 22 : if (optgr1==0) optgr1=dtset%pawstgylm
682 : end if
683 : call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfftf,psps%ntypat,&
684 : & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,&
685 112 : & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab )
686 : ! Compute exp(iq.r) factors around the atoms
687 112 : if (.not.qeq0) then
688 28 : do iatom=1,my_natom
689 14 : iatom_tot=iatom; if (paral_atom) iatom_tot=mpi_enreg%my_atmtab(iatom)
690 14 : if (allocated(pawfgrtab(iatom)%expiqr)) then
691 14 : ABI_FREE(pawfgrtab(iatom)%expiqr)
692 : end if
693 42 : ABI_MALLOC(pawfgrtab(iatom)%expiqr,(2,pawfgrtab(iatom)%nfgd))
694 : call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,pawfgrtab(iatom)%nfgd,dtset%qptn,&
695 14 : & pawfgrtab(iatom)%rfgd,xred(:,iatom_tot))
696 28 : pawfgrtab(iatom)%expiqr_allocated=1
697 : end do
698 : end if
699 : ! 4-Allocate PAW specific arrays
700 524 : ABI_MALLOC(paw_an,(my_natom))
701 524 : ABI_MALLOC(paw_ij,(my_natom))
702 112 : call paw_an_nullify(paw_an)
703 112 : call paw_ij_nullify(paw_ij)
704 112 : has_kxc=0;nkxc1=0;cplex=1
705 1160 : has_dijnd=0; if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1
706 112 : has_diju=merge(0,1,dtset%usepawu==0)
707 112 : if (rfphon/=0.or.rfelfd==1.or.rfelfd==3.or.rfstrs/=0.or.rf2_dkde/=0) then
708 55 : has_kxc=1
709 55 : call pawxc_get_nkxc(nkxc1,dtset%nspden,dtset%xclevel)
710 : end if
711 112 : has_vhartree=0
712 112 : if(dtset%orbmag>0 .AND. dtset%pawspnorb > 0) has_vhartree=1
713 112 : usevxctau_paw=pawxc_get_usekden(dtset%ixc)
714 : call paw_an_init(paw_an,dtset%natom,dtset%ntypat,nkxc1,0,dtset%nspden,&
715 : & cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1,has_vhartree=has_vhartree,has_vxctau=usevxctau_paw,&
716 112 : & has_vxc_ex=1,has_kxc=has_kxc,mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
717 : call paw_ij_init(paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%pawspnorb,&
718 : & natom,dtset%ntypat,dtset%typat,pawtab,has_dij=1,has_dijhartree=1,has_dijnd=has_dijnd,&
719 : & has_dijso=1,has_dijU=has_diju,has_pawu_occ=1,has_exexch_pot=1,&
720 : & nucdipmom=dtset%nucdipmom,zora=dtset%zora,&
721 224 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
722 :
723 : else ! PAW vs NCPP
724 604 : usexcnhat=0;usecprj=0
725 604 : use_nhat_gga=.false.
726 604 : ABI_MALLOC(paw_an,(0))
727 604 : ABI_MALLOC(paw_ij,(0))
728 604 : ABI_MALLOC(pawfgrtab,(0))
729 : end if ! paw
730 :
731 2148 : ABI_MALLOC(rhog,(2,nfftf))
732 2864 : ABI_MALLOC(rhor,(nfftf,dtset%nspden))
733 2148 : ABI_MALLOC(taug,(2,nfftf*dtset%usekden))
734 2864 : ABI_MALLOC(taur,(nfftf,dtset%nspden*dtset%usekden))
735 :
736 : !>>> Initialize charge density
737 :
738 : ABI_NVTX_START_RANGE(NVTX_MKRHO)
739 716 : if (dtset%getden/=0.or.dtset%irdden/=0) then
740 : ! Choice 1: read charge density from a disk file and broadcast data
741 : ! This part is not compatible with MPI-FFT (note single_proc=.True. below)
742 86 : rdwrpaw=psps%usepaw ; if(ireadwf0/=0) rdwrpaw=0
743 : if (rdwrpaw/=0) then
744 : ABI_MALLOC(pawrhoij_read,(natom))
745 : call pawrhoij_nullify(pawrhoij_read)
746 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
747 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cpxocc=dtset%pawcpxocc)
748 : call pawrhoij_alloc(pawrhoij_read,cplex_rhoij,nspden_rhoij,dtset%nspinor,&
749 : & dtset%nsppol,dtset%typat,pawtab=pawtab)
750 : else
751 86 : ABI_MALLOC(pawrhoij_read,(0))
752 : end if
753 : ! Note MT july 2013: should we read rhoij from the density file?
754 : call read_rhor(dtfil%fildensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw,&
755 86 : & mpi_enreg,rhor,hdr_den,pawrhoij_read,spaceworld,check_hdr=hdr)
756 86 : etotal = hdr_den%etot
757 86 : call hdr_den%free()
758 : if (rdwrpaw/=0) then
759 : call pawrhoij_bcast(pawrhoij_read,hdr%pawrhoij,0,spaceworld)
760 : call pawrhoij_free(pawrhoij_read)
761 : end if
762 86 : ABI_FREE(pawrhoij_read)
763 : ! Compute up+down rho(G) by fft
764 86 : call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0)
765 :
766 : else
767 : ! Choice 2: obtain the charge density from read wfs
768 : ! Warning: in PAW, compensation density has to be added !
769 630 : tim_mkrho=4
770 630 : paw_dmft%use_dmft=0 ; paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented
771 630 : if (psps%usepaw==1) then
772 273 : ABI_MALLOC(rhowfg,(2,dtset%nfft))
773 364 : ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden))
774 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
775 91 : & rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
776 91 : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor)
777 91 : ABI_FREE(rhowfg)
778 91 : ABI_FREE(rhowfr)
779 : else
780 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
781 539 : & rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
782 : end if
783 :
784 : end if ! choice for charge density initialization
785 : ABI_NVTX_END_RANGE()
786 :
787 : !>>> Initialize kinetic energy density
788 716 : if (dtset%usekden==1) then
789 :
790 2 : if (dtset%getkden/=0.or.dtset%irdkden/=0) then
791 : ! Choice 1: read kinetic energy density from a disk file and broadcast data
792 : ! This part is not compatible with MPI-FFT (note single_proc=.True. below)
793 0 : rdwrpaw=0
794 0 : ABI_MALLOC(pawrhoij_read,(0))
795 : call read_rhor(dtfil%filkdensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw,&
796 0 : & mpi_enreg,taur,hdr_den,pawrhoij_read,spaceworld,check_hdr=hdr)
797 0 : call hdr_den%free()
798 0 : ABI_FREE(pawrhoij_read)
799 : ! Compute up+down tau(G) by fft
800 0 : call fourdp(1,taug,taur(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0)
801 :
802 : else
803 : ! Choice 2: obtain the kinetic energy density from read wfs
804 2 : tim_mkrho=4
805 2 : paw_dmft%use_dmft=0 ; paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented
806 2 : if (psps%usepaw==1) then
807 6 : ABI_MALLOC(rhowfg,(2,dtset%nfft))
808 8 : ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden))
809 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
810 2 : & rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
811 2 : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,taug,rhowfr,taur)
812 2 : ABI_FREE(rhowfg)
813 2 : ABI_FREE(rhowfr)
814 : else
815 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
816 0 : & taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
817 : end if
818 :
819 : end if ! choice for kinetic energy density initialization
820 : end if ! usekden
821 :
822 : !In PAW, compensation density has eventually to be added
823 716 : nhatgrdim=0;nhatdim=0
824 716 : ABI_MALLOC(nhatgr,(0,0,0))
825 716 : if (psps%usepaw==1.and. ((usexcnhat==0).or.(dtset%getden==0).or.dtset%xclevel==2)) then
826 106 : nhatdim=1
827 :
828 424 : ABI_MALLOC(nhat,(nfftf,dtset%nspden))
829 106 : call timab(558,1,tsec)
830 106 : nhatgrdim=0;if (dtset%xclevel==2.and.dtset%pawnhatxc>0) nhatgrdim=usexcnhat
831 106 : ider=2*nhatgrdim
832 106 : if (nhatgrdim>0) then
833 4 : ABI_FREE(nhatgr)
834 20 : ABI_MALLOC(nhatgr,(nfftf,dtset%nspden,3))
835 : end if
836 106 : izero=0;cplex=1;ipert=0;idir=0;qphon(:)=zero
837 : call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,natom,&
838 : & nfftf,ngfftf,nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,&
839 : & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, &
840 106 : & mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom)
841 106 : if (dtset%getden==0) then
842 3255350 : rhor(:,:)=rhor(:,:)+nhat(:,:)
843 91 : call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0)
844 : end if
845 106 : call timab(558,2,tsec)
846 : else
847 610 : ABI_MALLOC(nhat,(0,0))
848 : end if
849 :
850 : !The GS irrzon and phnons were only needed to symmetrize the GS density
851 716 : ABI_FREE(irrzon)
852 716 : ABI_FREE(phnons)
853 :
854 : !Will compute now the total potential
855 :
856 : !Compute local ionic pseudopotential vpsp and core electron density xccc3d
857 716 : n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf
858 2148 : ABI_MALLOC(xccc3d,(n3xccc))
859 2148 : ABI_MALLOC(vpsp,(nfftf))
860 716 : if(psps%usepaw==1) then
861 336 : ABI_MALLOC(xcctau3d,(n3xccc*dtset%usekden))
862 : else
863 604 : ABI_MALLOC(xcctau3d,(0))
864 : end if
865 :
866 : ! Handling GEMM nonlop use
867 : ! Not enabled by default for CPU and CUDA implementations
868 : ! Enabled if using OpenMP GPU offload
869 716 : gemm_nonlop_use_gemm = .false.
870 :
871 : ! OpenMP GPU offload case (GEMM nonlop used by default)
872 716 : if(dtset%gpu_option == ABI_GPU_OPENMP .or. dtset%use_gemm_nonlop == 1) then
873 0 : gemm_nonlop_use_gemm = .true.
874 0 : call init_gemm_nonlop(dtset%gpu_option)
875 : end if
876 :
877 716 : gemm_nonlop_is_distributed = .false.
878 716 : if(dtset%gpu_nl_distrib == 1) gemm_nonlop_is_distributed = .true.
879 :
880 : !Determine by which method the local ionic potential and/or
881 : ! the pseudo core charge density have to be computed
882 : !Local ionic potential:
883 : ! Method 1: PAW ; Method 2: Norm-conserving PP
884 716 : vloc_method=1;if (psps%usepaw==0) vloc_method=2
885 : !Pseudo core charge density:
886 : ! Method 1: PAW, nc_xccc_gspace ; Method 2: Norm-conserving PP
887 : coredens_method=1;if (psps%usepaw==0) coredens_method=2
888 716 : if (psps%nc_xccc_gspace==1) coredens_method=1
889 716 : if (psps%nc_xccc_gspace==0) coredens_method=2
890 :
891 : ! Core kinetic energy density method
892 716 : coretau_method = 0
893 716 : if (dtset%usekden==1.and.psps%usepaw==1) then
894 2 : coretau_method=1
895 2 : if (psps%nc_xccc_gspace==0) then
896 0 : coretau_method=2
897 : end if
898 : end if
899 :
900 :
901 : !Local ionic potential and/or pseudo core charge by method 1
902 716 : if (vloc_method==1.or.coredens_method==1) then
903 113 : call timab(562,1,tsec)
904 113 : optv=0;if (vloc_method==1) optv=1
905 113 : optn=0;if (coredens_method==1) optn=n3xccc/nfftf
906 113 : optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1
907 : call atm2fft(atindx1,xccc3d,vpsp,dum_dyfrn,dum_dyfrv,dum_eltfrxc,dum_gauss,gmet,gprimd,&
908 : & dum_grn,dum_grv,gsqcut,mgfftf,psps%mqgrid_vl,natom,nattyp,nfftf,ngfftf,&
909 : & ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1df,psps%qgrid_vl,&
910 113 : & dtset%qprtrb,dtset%rcut,dum_rhog,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dum_vg,dum_vg,dum_vg,dtset%vprtrb,psps%vlspl)
911 113 : call timab(562,2,tsec)
912 : end if
913 :
914 716 : if (coretau_method==1) then
915 2 : optv=0;optn=1
916 2 : optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=4
917 : call atm2fft(atindx1,xcctau3d,dummy_out6,dummy_out1,dummy_out2,dummy_out3,dummy_in,&
918 : & gmet,gprimd,dummy_out4,dummy_out5,gsqcut,mgfftf,psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,ntypat,&
919 : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1df,psps%qgrid_vl,dtset%qprtrb,&
920 : & dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,&
921 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
922 2 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
923 : end if
924 :
925 : !Local ionic potential by method 2
926 716 : if (vloc_method==2) then
927 604 : option=1
928 1812 : ABI_MALLOC(dyfrlo_indx,(3,3,natom))
929 1812 : ABI_MALLOC(grtn_indx,(3,natom))
930 : call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,&
931 : & grtn_indx,gsqcut,dummy6,mgfftf,mpi_enreg,natom,nattyp,&
932 : & nfftf,ngfftf,dtset%nspden,ntypat,option,pawtab,ph1df,psps,&
933 604 : & dtset%qprtrb,rhog,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred)
934 604 : ABI_FREE(dyfrlo_indx)
935 604 : ABI_FREE(grtn_indx)
936 : end if
937 :
938 : !Pseudo core electron density by method 2
939 716 : if (coredens_method==2.and.psps%n1xccc/=0) then
940 356 : option=1
941 1068 : ABI_MALLOC(dyfrx2,(3,3,natom))
942 356 : ABI_MALLOC(vxc,(0,0)) ! dummy
943 : call mkcore(dummy6,dyfrx2,grxc,mpi_enreg,natom,nfftf,dtset%nspden,ntypat,&
944 : & ngfftf(1),psps%n1xccc,ngfftf(2),ngfftf(3),option,rprimd,dtset%typat,ucvol,vxc,&
945 356 : & psps%xcccrc,psps%xccc1d,xccc3d,xred)
946 356 : ABI_FREE(dyfrx2)
947 356 : ABI_FREE(vxc) ! dummy
948 : end if
949 :
950 : !Set up hartree and xc potential. Compute kxc here.
951 1432 : ABI_MALLOC(vhartr,(nfftf))
952 :
953 : call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfftf,ngfftf,&
954 716 : &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr)
955 :
956 716 : option=2 ; nk3xc=1
957 716 : nkxc=2*min(dtset%nspden,2)-1;if(dtset%xclevel==2)nkxc=12*min(dtset%nspden,2)-5
958 716 : call check_kxc(dtset%ixc,dtset%optdriver)
959 2864 : ABI_MALLOC(kxc,(nfftf,nkxc))
960 2864 : ABI_MALLOC(vxc,(nfftf,dtset%nspden))
961 716 : usevxctau=merge(1,0,xc_need_kden(dtset%ixc))
962 3582 : ABI_MALLOC(vxctau,(nfftf,dtset%nspden,4*usevxctau))
963 :
964 716 : call xcdata_init(xcdata,dtset=dtset)
965 716 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
966 : !If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value
967 716 : if (dtset%xc_tb09_c>99._dp) then
968 : call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom,nfftf,ngfftf, &
969 : & nhat,psps%usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, &
970 : & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,psps%usepaw, &
971 716 : & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
972 : end if
973 :
974 : call rhotoxc(bigexc,bigsxc,kxc,mpi_enreg,nfftf,ngfftf,&
975 : & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,rhor,&
976 : & rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata,&
977 716 : & taur=taur,vhartr=vhartr,vxctau=vxctau,xcctau3d=xcctau3d)
978 :
979 : !Compute local + Hxc potential, and subtract mean potential.
980 2148 : ABI_MALLOC(vtrial,(nfftf,dtset%nspden))
981 1513 : do ispden=1,min(dtset%nspden,2)
982 8236460 : do ifft=1,nfftf
983 8235744 : vtrial(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft)
984 : end do
985 : end do
986 716 : if (dtset%nspden==4) then
987 60 : do ispden=3,4
988 234338 : do ifft=1,nfftf
989 234318 : vtrial(ifft,ispden)=vxc(ifft,ispden)
990 : end do
991 : end do
992 : end if
993 :
994 716 : ABI_FREE(vhartr)
995 :
996 716 : if(dtset%prtvol==-level) call wrtout(std_out,' respfn: ground-state density and potential set up.')
997 :
998 : !PAW: compute Dij quantities (psp strengths)
999 716 : if (psps%usepaw==1)then
1000 112 : cplex=1;ipert=0;option=1
1001 112 : nzlmopt=0;if (dtset%pawnzlm>0) nzlmopt=-1
1002 :
1003 : call pawdenpot(compch_sph,el_temp,gprimd,ipert,dtset%ixc,my_natom,natom,dtset%nspden,&
1004 : & ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_energies,paw_ij,pawang,&
1005 : & dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,&
1006 : & dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,xred,ucvol,psps%znuclpsp,dtset%spinaxis, &
1007 112 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1008 :
1009 112 : call timab(561,1,tsec)
1010 : call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,natom,nfftf,nfftotf,&
1011 : & dtset%nspden,ntypat,paw_an,paw_ij,pawang,pawfgrtab,dtset%pawprtvol,&
1012 : & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,k0,&
1013 : & dtset%spnorbscl,ucvol,dtset%cellcharge(1),vtrial,vxc,xred,dtset%znucl,&
1014 : & nucdipmom=dtset%nucdipmom,&
1015 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,&
1016 112 : & spinaxis=dtset%spinaxis)
1017 : call symdij(gprimd,indsym,ipert,my_natom,natom,dtset%nsym,ntypat,0,&
1018 : & paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,&
1019 112 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1020 112 : call timab(561,2,tsec)
1021 :
1022 : end if
1023 :
1024 : !-----2. Frozen-wavefunctions and Ewald(q=0) parts of 2DTE
1025 :
1026 716 : dyfr_nondiag=0;if (psps%usepaw==1.and.rfphon==1) dyfr_nondiag=1
1027 716 : dyfr_cplex=1;if (psps%usepaw==1.and.rfphon==1.and.(.not.qeq0)) dyfr_cplex=2
1028 716 : if (abs(dtset%rfomega)>tol10.or.dtset%tim1rev==0) dyfr_cplex=2
1029 3580 : ABI_MALLOC(dyew,(2,3,natom,3,natom))
1030 2148 : ABI_MALLOC(dyewq0,(3,3,natom))
1031 1432 : ABI_MALLOC(dyfrlo,(3,3,natom))
1032 1432 : ABI_MALLOC(dyfrx2,(3,3,natom))
1033 4296 : ABI_MALLOC(dyfrnl,(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag))
1034 2864 : ABI_MALLOC(dyfrwf,(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag))
1035 2864 : ABI_MALLOC(becfrnl,(3,natom,3*pawbec))
1036 2148 : ABI_MALLOC(piezofrnl,(6,3*pawpiezo))
1037 2864 : ABI_MALLOC(dyvdw,(2,3,natom,3,natom*usevdw))
1038 142154 : dyew(:,:,:,:,:)=zero
1039 20567 : dyewq0(:,:,:)=zero
1040 45981 : dyfrnl(:,:,:,:,:)=zero
1041 45981 : dyfrwf(:,:,:,:,:)=zero
1042 20567 : dyfrlo(:,:,:)=zero
1043 20567 : dyfrx2(:,:,:)=zero
1044 912 : if (usevdw==1) dyvdw(:,:,:,:,:)=zero
1045 2621 : if (pawbec==1) becfrnl(:,:,:)=zero
1046 1178 : if (pawpiezo==1) piezofrnl(:,:)=zero
1047 :
1048 716 : ABI_MALLOC(eltcore,(6,6))
1049 2864 : ABI_MALLOC(elteew,(6+3*natom,6))
1050 716 : ABI_MALLOC(eltfrhar,(6,6))
1051 1432 : ABI_MALLOC(eltfrnl,(6+3*natom,6))
1052 1432 : ABI_MALLOC(eltfrloc,(6+3*natom,6))
1053 716 : ABI_MALLOC(eltfrkin,(6,6))
1054 1432 : ABI_MALLOC(eltfrxc,(6+3*natom,6))
1055 2864 : ABI_MALLOC(eltvdw,(6+3*natom,6*usevdw))
1056 30788 : eltcore(:,:)=zero
1057 58274 : elteew(:,:)=zero
1058 58274 : eltfrnl(:,:)=zero
1059 58274 : eltfrloc(:,:)=zero
1060 30788 : eltfrkin(:,:)=zero
1061 30788 : eltfrhar(:,:)=zero
1062 58274 : eltfrxc(:,:)=zero
1063 914 : if (usevdw==1) eltvdw(:,:)=zero
1064 :
1065 : !Section common to all perturbations
1066 : !Compute the nonlocal part of the elastic tensor and/or dynamical matrix
1067 716 : if (rfstrs/=0.or.rfphon==1.or.dtset%efmas>0.or.pawbec==1.or.pawpiezo==1)then
1068 : call d2frnl(becfrnl,cg,dtfil,dtset,dyfrnl,dyfr_cplex,&
1069 : & dyfr_nondiag,efmasdeg,efmasval,eigen0,eltfrnl,gsqcut,has_allddk,indsym,kg,&
1070 : & dtset%mband_mem,dtset%mkmem,mgfftf,&
1071 : & mpi_enreg,psps%mpsang,my_natom,natom,nfftf,ngfft,ngfftf,&
1072 : & npwarr,occ,paw_ij,pawang,pawbec,pawfgrtab,pawpiezo,pawrad,&
1073 : & pawrhoij,pawtab,ph1d,ph1df,piezofrnl,psps,rprimd,rfphon,&
1074 487 : & rfstrs,symrec,vtrial,vxc,xred,ylm,ylmgr)
1075 : end if
1076 :
1077 : !No more need of these local derivatives
1078 716 : if (rfphon==1.and.psps%usepaw==1.and.(.not.use_nhat_gga)) then
1079 113 : do iatom=1,my_natom
1080 70 : if (allocated(pawfgrtab(iatom)%gylmgr2)) then
1081 70 : ABI_FREE(pawfgrtab(iatom)%gylmgr2)
1082 : end if
1083 113 : pawfgrtab(iatom)%gylmgr2_allocated=0
1084 : end do
1085 : end if
1086 :
1087 : !Section for the atomic displacement/electric field perturbations
1088 716 : if (rfphon==1) then
1089 :
1090 : ! Compute the local of the dynamical matrix
1091 : ! dyfrnl has not yet been symmetrized, but will be in the next routine
1092 : call dfpt_dyfro(atindx1,dyfrnl,dyfrlo,dyfrwf,dyfrx2,dyfr_cplex,dyfr_nondiag,&
1093 : & dtset,gmet,gprimd,gsqcut,indsym,mgfftf,mpi_enreg,psps%mqgrid_vl,&
1094 : & natom,nattyp, nfftf,ngfftf,dtset%nspden,dtset%nsym,ntypat,&
1095 : & psps%n1xccc,n3xccc,psps,pawtab,ph1df,psps%qgrid_vl,&
1096 : & dtset%qptn,rhog,rprimd,symq,symrec,dtset%typat,ucvol,&
1097 450 : & psps%usepaw,psps%vlspl,vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred)
1098 :
1099 : ! Compute Ewald (q=0) contribution
1100 450 : sumg0=0;qphon(:)=zero
1101 : call dfpt_ewald(dyew,gmet,gsqcut,dtset%icutcoul,my_natom,natom,ngfftf,dtset%nkpt,qphon,dtset%rcut,rmet,&
1102 : & rprimd,sumg0,dtset%typat,ucvol,dtset%vcutgeo,xred,psps%ziontypat,&
1103 450 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1104 450 : option=1
1105 450 : call q0dy3_calc(natom,dyewq0,dyew,option)
1106 : ! Calculate the DFT-D2 vdw part of the dynamical matrix
1107 450 : if (usevdw==1.and.dtset%vdw_xc==5) then
1108 : call vdw_dftd2(evdw,dtset%ixc,natom,ntypat,1,dtset%typat,rprimd,dtset%vdw_tol,&
1109 0 : & xred,psps%znuclpsp,dyn_vdw_dftd2=dyvdw,qphon=dtset%qptn)
1110 : end if
1111 : ! Calculate the DFT-D3/D3(BJ) vdw part of the dynamical matrix
1112 450 : if (usevdw==1.and.(dtset%vdw_xc==6.or.dtset%vdw_xc==7)) then
1113 : call vdw_dftd3(evdw,dtset%ixc,natom,ntypat,1,dtset%typat,rprimd,&
1114 : & dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,xred,psps%znuclpsp,&
1115 3 : & dyn_vdw_dftd3=dyvdw,qphon=dtset%qptn)
1116 : end if
1117 : !The frozen-wavefunction part of the dynamical matrix is now:
1118 : ! d2frnl: non-local contribution
1119 : ! dyfrlo: local contribution
1120 : ! dyfrx2: 2nd order xc core correction contribution
1121 : ! dyew : Ewald contribution
1122 : ! dyvdw : vdw DFT-D contribution
1123 : ! dyfrwf: all contributions
1124 : ! In case of PAW, it misses a term coming from the perturbed overlap operator
1125 : end if
1126 :
1127 : !Section for the strain perturbation
1128 716 : if(rfstrs/=0) then
1129 :
1130 : ABI_NVTX_START_RANGE(NVTX_DFPT_ELT)
1131 : ! Verify that k-point set has full space-group symmetry; otherwise exit
1132 37 : timrev=1
1133 37 : if (symkchk(dtset%kptns,dtset%nkpt,dtset%nsym,symrec,timrev,msg) /= 0) then
1134 0 : ABI_ERROR(msg)
1135 : end if
1136 :
1137 : ! Calculate the kinetic part of the elastic tensor
1138 : call dfpt_eltfrkin(cg,eltfrkin,dtset%ecut,dtset%ecutsm,dtset%effmass_free,&
1139 : & dtset%istwfk,kg,dtset%kptns,dtset%mband,dtset%mband_mem,dtset%mgfft,dtset%mkmem,mpi_enreg,&
1140 : & dtset%mpw,dtset%nband,dtset%nkpt,ngfft,npwarr,&
1141 37 : & dtset%nspinor,dtset%nsppol,occ,rprimd,dtset%wtk)
1142 :
1143 : ! Calculate the hartree part of the elastic tensor
1144 37 : call dfpt_eltfrhar(eltfrhar,rprimd,gsqcut,mpi_enreg,nfftf,ngfftf,rhog)
1145 :
1146 : ! Calculate the xc part of the elastic tensor
1147 : call dfpt_eltfrxc(atindx,dtset,eltfrxc,bigexc,gsqcut,kxc,mpi_enreg,mgfftf,&
1148 : & nattyp,nfftf,ngfftf,ngfftf,nhat,nkxc,n3xccc,pawtab,ph1df,psps,rhor,rprimd,&
1149 37 : & usexcnhat,vxc,xccc3d,xred)
1150 :
1151 : ! Calculate the local potential part of the elastic tensor
1152 : call dfpt_eltfrloc(atindx,eltfrloc,gmet,gprimd,gsqcut,mgfftf,mpi_enreg,psps%mqgrid_vl,&
1153 37 : & natom,nattyp,nfftf,ngfftf,ntypat,ph1df,psps%qgrid_vl,rhog,psps%vlspl)
1154 :
1155 : ! Calculate the Ewald part of the elastic tensor
1156 : call elt_ewald(elteew,gmet,gprimd,my_natom,natom,ntypat,rmet,rprimd,&
1157 : & dtset%typat,ucvol,xred,psps%ziontypat,&
1158 37 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1159 :
1160 : ! Calculate the DFT-D2 vdw part of the elastic tensor
1161 37 : if (usevdw==1.and.dtset%vdw_xc==5) then
1162 0 : option=1; if (rfphon==1) option=0
1163 : call vdw_dftd2(evdw,dtset%ixc,natom,ntypat,option,dtset%typat,rprimd,dtset%vdw_tol,&
1164 0 : & xred,psps%znuclpsp,elt_vdw_dftd2=eltvdw)
1165 : end if
1166 : ! Calculate the DFT-D3/D3(BJ) vdw part of the elastic tensor
1167 37 : if (usevdw==1.and.(dtset%vdw_xc==6.or.dtset%vdw_xc==7)) then
1168 1 : option=1; if (rfphon==1) option=0
1169 : call vdw_dftd3(evdw,dtset%ixc,natom,ntypat,option,dtset%typat,rprimd,&
1170 : & dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,xred,psps%znuclpsp,&
1171 1 : & elt_vdw_dftd3=eltvdw)
1172 : end if
1173 : ! Calculate the psp core energy part of elastic tensor (trivial)
1174 481 : eltcore(1:3,1:3)=ecore/ucvol
1175 :
1176 : !The frozen-wavefunction part of the elastic tensor is now:
1177 : ! eltfrnl: non-local contribution
1178 : ! eltfrloc: local contribution
1179 : ! eltfrkin: kinetic contribution
1180 : ! eltfrhar: Hartree contribution
1181 : ! eltfrx: XC contribution
1182 : ! eltcore: psps core contribution
1183 : ! elteew: Ewald contribution
1184 : ! eltvdw: vdw DFT-D contribution
1185 : ! In case of PAW, it misses a term coming from the perturbed overlap operator
1186 : ABI_NVTX_END_RANGE()
1187 : end if
1188 :
1189 716 : ABI_FREE(vpsp)
1190 716 : ABI_FREE(xccc3d)
1191 716 : if(allocated(xcctau3d)) then
1192 716 : ABI_FREE(xcctau3d)
1193 : end if
1194 :
1195 716 : if(dtset%prtvol==-level) call wrtout(std_out,' respfn: frozen wavef. and Ewald(q=0) part of 2DTE done.')
1196 :
1197 :
1198 716 : call timab(136,2,tsec)
1199 716 : call pstat_proc%print(_PSTAT_ARGS_)
1200 :
1201 : !-----3. Initialisation of 1st response, taking into account the q vector.
1202 :
1203 716 : call timab(137,1,tsec)
1204 :
1205 716 : write(msg,'(3a)')ch10,' ==> initialize data related to q vector <== ',ch10
1206 2148 : call wrtout([std_out, ab_out] ,msg)
1207 :
1208 2864 : qphon(:)=dtset%qptn(:)
1209 716 : sumg0=1
1210 :
1211 : !Treat the case of q=0 or q too close to 0
1212 716 : qzero=0
1213 716 : if(qeq0)then
1214 501 : qphon(:)=zero
1215 : write(msg,'(3a)')&
1216 501 : & ' respfn : the norm of the phonon wavelength (as input) was small (<1.d-7).',ch10,&
1217 1002 : & ' q has been set exactly to (0 0 0)'
1218 501 : call wrtout(std_out,msg)
1219 501 : sumg0=0
1220 501 : qzero=1
1221 : else
1222 215 : if(rfelfd/=0 .or. rfstrs/=0 .or. rfddk /= 0 .or. rf2_dkdk /= 0 .or. rf2_dkde /= 0) then
1223 : ! Temporarily, ...
1224 0 : write(msg, '(a,a,a,3es16.6,a,a,5(a,i2),a,a,a)' )ch10,&
1225 0 : & 'The treatment of non-zero wavevector q is restricted to phonons.',&
1226 0 : & 'However, the input normalized qpt is',qphon(:),',',ch10,&
1227 0 : & 'while rfelfd=',rfelfd,', rfddk=',rfddk,', rf2_dkdk=',rf2_dkdk,', rf2_dkde=',rf2_dkde,&
1228 0 : & ' and rfstrs=',rfstrs,'.',ch10,&
1229 0 : & 'Action: change qpt, or rfelfd, or rfstrs in the input file.'
1230 0 : ABI_ERROR(msg)
1231 : end if
1232 215 : if(chneut/=0)then
1233 165 : write(msg,'(2a)')ch10,' chneut/=0 not allowed with q/=0 => chneut reset to 0 locally.'
1234 165 : ABI_WARNING(msg)
1235 165 : chneut=0
1236 : end if
1237 215 : if(asr/=0)then
1238 162 : write(msg,'(2a)')ch10,' asr/=0 not allowed with q/=0 => asr reset to 0 locally.'
1239 162 : ABI_WARNING(msg)
1240 162 : asr=0
1241 : end if
1242 :
1243 : end if
1244 :
1245 : !Determine the symmetrical perturbations
1246 2148 : ABI_MALLOC(pertsy,(3,mpert))
1247 716 : call irreducible_set_pert(indsym,mpert,natom,dtset%nsym,pertsy,rfdir,rfpert,symq,symrec,dtset%symrel)
1248 :
1249 716 : write(msg,'(a)') ' The list of irreducible perturbations for this q vector is:'
1250 2148 : call wrtout([std_out, ab_out] ,msg)
1251 716 : ii=1
1252 7425 : do ipert=1,mpert
1253 27552 : do idir=1,3
1254 26836 : if(rfpert(ipert)==1.and.rfdir(idir)==1)then
1255 3142 : if( pertsy(idir,ipert)==1 )then
1256 2181 : write(msg, '(i5,a,i2,a,i4)' )ii,') idir=',idir,' ipert=',ipert
1257 6543 : call wrtout([std_out, ab_out] ,msg)
1258 2181 : ii=ii+1
1259 : end if
1260 : end if
1261 : end do
1262 : end do
1263 :
1264 : !test if the user left default rfdir 0 0 0
1265 716 : if (ii==1 .and. rf2_dkdk==0 .and. rf2_dkde==0) then
1266 0 : write(msg,'(5a)')ch10,&
1267 0 : & ' WARNING: no perturbations to be done at this q-point.',ch10,&
1268 0 : & ' You may have forgotten to set the rfdir or rfatpol variables. Continuing normally.',ch10
1269 0 : call wrtout(ab_out,msg)
1270 0 : ABI_WARNING(msg)
1271 : end if
1272 :
1273 716 : if (dtset%prepanl==1.and.(rf2_dkdk/=0 .or. rf2_dkde/=0)) then
1274 98 : ABI_MALLOC(rfpert_nl,(3,natom+2,3,natom+2,3,natom+2))
1275 35182 : rfpert_nl = 0
1276 560 : rfpert_nl(:,natom+2,:,natom+2,:,natom+2) = 1
1277 1190 : rfpert_nl(:,1:natom,:,natom+2,:,natom+2) = 1
1278 1148 : rfpert_nl(:,natom+2,:,1:natom,:,natom+2) = 1
1279 1134 : rfpert_nl(:,natom+2,:,natom+2,:,1:natom) = 1
1280 14 : call sytens(indsym,natom+2,natom,dtset%nsym,rfpert_nl,symrec,dtset%symrel)
1281 14 : write(msg, '(a,a,a,a,a)' ) ch10, &
1282 14 : & ' The list of irreducible elements of the Raman and non-linear',&
1283 28 : & ch10,' optical susceptibility tensors is:',ch10
1284 14 : call wrtout(std_out,msg)
1285 :
1286 14 : write(msg,'(12x,a)')'i1pert i1dir i2pert i2dir i3pert i3dir'
1287 14 : call wrtout(std_out,msg)
1288 14 : n1 = 0
1289 14 : rf2_dirs_from_rfpert_nl(:,:) = 0
1290 70 : do i1pert = 1, natom + 2
1291 238 : do i1dir = 1, 3
1292 896 : do i2pert = 1, natom + 2
1293 2856 : do i2dir = 1, 3
1294 10752 : do i3pert = 1, natom + 2
1295 34272 : do i3dir = 1,3
1296 32256 : if (rfpert_nl(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
1297 582 : n1 = n1 + 1
1298 582 : write(msg,'(2x,i4,a,6(5x,i3))') n1,')', i1pert,i1dir,i2pert,i2dir,i3pert,i3dir
1299 582 : call wrtout(std_out,msg)
1300 582 : if (i2pert==natom+2) then
1301 398 : if (i3pert==natom+2) then
1302 398 : rf2_dirs_from_rfpert_nl(i3dir,i2dir) = 1
1303 0 : else if (i1pert==natom+2) then
1304 0 : rf2_dirs_from_rfpert_nl(i1dir,i2dir) = 1
1305 : end if
1306 : end if
1307 : end if
1308 : end do
1309 : end do
1310 : end do
1311 : end do
1312 : end do
1313 : end do
1314 14 : write(msg,'(a,a)') ch10,ch10
1315 14 : call wrtout(std_out,msg)
1316 :
1317 14 : call wrtout(std_out,'rf2_dirs_from_rfpert_nl :')
1318 56 : do i1dir = 1, 3
1319 884 : do i2dir = 1, 3
1320 126 : write(msg,'(3(a,i1))') ' ',i1dir,' ',i2dir,' : ',rf2_dirs_from_rfpert_nl(i1dir,i2dir)
1321 168 : call wrtout(std_out,msg)
1322 : end do
1323 : end do
1324 : end if
1325 :
1326 : !For longwave calculation:
1327 : !Get symmetries in cartesian coordinates
1328 2148 : ABI_MALLOC(symrel_cart, (3, 3, dtset%nsym))
1329 17797 : do isym =1,dtset%nsym
1330 17081 : call symredcart(rprimd, gprimd, symrel_cart(:,:,isym), dtset%symrel(:,:,isym))
1331 : ! purify operations in cartesian coordinates.
1332 222769 : where (abs(symrel_cart(:,:,isym)) < tol14)
1333 : symrel_cart(:,:,isym) = zero
1334 : end where
1335 : end do
1336 :
1337 716 : if (dtset%prepalw/=0) then
1338 112 : ABI_MALLOC(rfpert_lw,(3,natom+8,3,natom+8,3,natom+8))
1339 782952 : rfpert_lw=0
1340 16 : if (dtset%prepalw==1) then
1341 948 : rfpert_lw(:,1:natom+2,:,1:natom,:,natom+8)=1
1342 948 : rfpert_lw(:,1:natom+2,:,natom+3:natom+4,:,natom+8)=1
1343 13 : else if (dtset%prepalw==2) then
1344 1101 : rfpert_lw(:,natom+2,:,1:natom,:,natom+8)=1
1345 1 : else if (dtset%prepalw==3) then
1346 0 : rfpert_lw(:,1:natom+2,:,1:natom,:,natom+8)=1
1347 1 : else if (dtset%prepalw==4) then
1348 40 : rfpert_lw(:,natom+2,:,natom+2,:,natom+8)=1
1349 : end if
1350 :
1351 : ! call sylwtens(indsym,natom+8,natom,dtset%nsym,rfpert_lw,symrec,dtset%symrel,symrel_cart)
1352 16 : call sylwtens(indsym,natom+8,natom,dtset%nsym,rfpert_lw,symrec,dtset%symrel)
1353 :
1354 16 : write(msg,'(7a)') ch10, ' The following reducible perturbations will also be ', ch10, &
1355 16 : & ' explicitly calculated for a correct subsequent ', ch10, &
1356 32 : & ' execution of the longwave driver:', ch10
1357 16 : call wrtout(ab_out,msg,'COLL')
1358 16 : call wrtout(std_out,msg,'COLL')
1359 186 : do i3pert = 1,natom+8
1360 696 : do i3dir = 1, 3
1361 6254 : do i2pert = 1, natom+8
1362 22806 : do i2dir = 1,3
1363 212286 : do i1pert = 1,natom+8
1364 776682 : do i1dir = 1, 3
1365 759960 : if (rfpert_lw(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
1366 930 : if (pertsy(i1dir,i1pert)==-1) then
1367 12 : pertsy(i1dir,i1pert)=1
1368 12 : write(msg,'(a,i2,a,i4)' )' idir=',i1dir,' ipert=',i1pert
1369 12 : call wrtout(ab_out,msg,'COLL')
1370 12 : call wrtout(std_out,msg,'COLL')
1371 : end if
1372 930 : if (pertsy(i2dir,i2pert)==-1) then
1373 16 : pertsy(i2dir,i2pert)=1
1374 16 : write(msg,'(a,i2,a,i4)' )' idir=',i2dir,' ipert=',i2pert
1375 16 : call wrtout(ab_out,msg,'COLL')
1376 16 : call wrtout(std_out,msg,'COLL')
1377 : end if
1378 : end if
1379 : end do
1380 : end do
1381 : end do
1382 : end do
1383 : end do
1384 : end do
1385 16 : write(msg,'(a,a)') ch10,ch10
1386 16 : call wrtout(std_out,msg,'COLL')
1387 16 : ABI_FREE(rfpert_lw)
1388 : end if
1389 :
1390 : !Contribution to the dynamical matrix from ion-ion energy
1391 716 : if(rfphon==1)then
1392 450 : sumg0_save= sumg0
1393 450 : if (dtset%icutcoul==55) sumg0= 0
1394 : call dfpt_ewald(dyew,gmet,gsqcut,dtset%icutcoul,my_natom,natom,ngfftf,dtset%nkpt,qphon,dtset%rcut,rmet,&
1395 : & rprimd,sumg0,dtset%typat,ucvol,dtset%vcutgeo,xred,psps%ziontypat, &
1396 450 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1397 450 : sumg0= sumg0_save
1398 450 : call q0dy3_apply(natom,dyewq0,dyew)
1399 : end if
1400 :
1401 : !1-order contribution of the xc core correction to the dynamical matrix
1402 3580 : ABI_MALLOC(dyfrx1,(2,3,natom,3,natom))
1403 142154 : dyfrx1(:,:,:,:,:)=zero
1404 716 : if(rfphon==1.and.psps%n1xccc/=0)then
1405 : ABI_NVTX_START_RANGE(NVTX_DFPT_DYXC)
1406 1435 : ABI_MALLOC(blkflgfrx1,(3,natom,3,natom))
1407 : !FR non-collinear magnetism
1408 287 : if (dtset%nspden==4) then
1409 : call dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,dtset%ixc,kxc,mgfftf,mpert,mpi_enreg,&
1410 : & psps%mqgrid_vl,natom,nfftf,ngfftf,nkxc,non_magnetic_xc,dtset%nspden,&
1411 : & ntypat,psps%n1xccc,psps,pawtab,ph1df,psps%qgrid_vl,qphon,&
1412 6 : & rfdir,rfpert,rprimd,timrev,dtset%typat,ucvol,psps%usepaw,psps%xcccrc,psps%xccc1d,xred,rhor=rhor,vxc=vxc,ixcrot=dtset%ixcrot)
1413 : else
1414 : call dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,dtset%ixc,kxc,mgfftf,mpert,mpi_enreg,&
1415 : & psps%mqgrid_vl,natom,nfftf,ngfftf,nkxc,non_magnetic_xc,dtset%nspden,&
1416 : & ntypat,psps%n1xccc,psps,pawtab,ph1df,psps%qgrid_vl,qphon,&
1417 281 : & rfdir,rfpert,rprimd,timrev,dtset%typat,ucvol,psps%usepaw,psps%xcccrc,psps%xccc1d,xred)
1418 : end if
1419 : ABI_NVTX_END_RANGE()
1420 : end if
1421 :
1422 : !Deallocate the arrays that were needed only for the frozen wavefunction part
1423 716 : ABI_FREE(ph1d)
1424 716 : ABI_FREE(ph1df)
1425 716 : ABI_FREE(cg)
1426 716 : ABI_FREE(kg)
1427 716 : ABI_FREE(npwarr)
1428 716 : if(xmpi_paral==1) then
1429 716 : ABI_FREE(mpi_enreg%proc_distrb)
1430 : end if
1431 :
1432 3580 : ABI_MALLOC(blkflg,(3,mpert,3,mpert))
1433 3580 : ABI_MALLOC(d2eig0,(2,3,mpert,3,mpert))
1434 2148 : ABI_MALLOC(d2k0,(2,3,mpert,3,mpert))
1435 2148 : ABI_MALLOC(d2lo,(2,3,mpert,3,mpert))
1436 2148 : ABI_MALLOC(d2loc0,(2,3,mpert,3,mpert))
1437 2148 : ABI_MALLOC(d2nfr,(2,3,mpert,3,mpert))
1438 2148 : ABI_MALLOC(d2nl,(2,3,mpert,3,mpert))
1439 2148 : ABI_MALLOC(d2nl0,(2,3,mpert,3,mpert))
1440 2148 : ABI_MALLOC(d2nl1,(2,3,mpert,3,mpert))
1441 2148 : ABI_MALLOC(d2vn,(2,3,mpert,3,mpert))
1442 2864 : ABI_MALLOC(d2ovl,(2,3,mpert,3,mpert*psps%usepaw))
1443 816828 : blkflg(:,:,:,:)=0
1444 4000768 : d2eig0(:,:,:,:,:)=zero ; d2k0(:,:,:,:,:)=zero
1445 4000768 : d2lo(:,:,:,:,:)=zero ; d2loc0(:,:,:,:,:)=zero
1446 4000768 : d2nfr(:,:,:,:,:)=zero ; d2nl(:,:,:,:,:)=zero
1447 4000768 : d2nl0(:,:,:,:,:)=zero ; d2nl1(:,:,:,:,:)=zero
1448 2000742 : d2vn(:,:,:,:,:)=zero
1449 335358 : if (psps%usepaw==1) d2ovl(:,:,:,:,:)=zero
1450 :
1451 716 : prtbbb=dtset%prtbbb
1452 3580 : ABI_MALLOC(d2bbb,(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb))
1453 2864 : ABI_MALLOC(d2cart_bbb,(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb))
1454 716 : if(prtbbb==1)then
1455 71716 : d2cart_bbb(:,:,:,:,:,:)=zero
1456 71716 : d2bbb(:,:,:,:,:,:)=zero
1457 : end if
1458 :
1459 716 : dim_eig2nkq = 0
1460 716 : if(dtset%ieig2rf /= 0) dim_eig2nkq = 1
1461 5728 : ABI_MALLOC(eig2nkq,(2,dtset%mband*dtset%nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq))
1462 716 : dim_eigbrd=0
1463 716 : if(dtset%ieig2rf /= 0 .and. dtset%smdelta>0 ) dim_eigbrd = 1
1464 4296 : ABI_MALLOC(eigbrd,(2,dtset%mband*dtset%nsppol,dtset%nkpt,3,natom,3,natom*dim_eigbrd))
1465 :
1466 716 : call timab(137,2,tsec)
1467 716 : call pstat_proc%print(_PSTAT_ARGS_)
1468 :
1469 :
1470 : !Check whether exiting was required by the user.
1471 : !If found then do not start minimization steps
1472 : !At this first call to exit_check, initialize cpus
1473 716 : cpus=dtset%cpus
1474 716 : if(abs(cpus)>1.0d-5)cpus=cpus+cpui
1475 716 : openexit=1; if(dtset%chkexit==0) openexit=0
1476 716 : call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit)
1477 :
1478 : !TEMPORARY: for testing purpose only
1479 : ! if (rfstrs/=0.and.dtset%usepaw==1) iexit=1
1480 :
1481 716 : if (iexit==0) then
1482 : ! #######################################################################
1483 57996 : write(msg,'(a,80a)')ch10,('=',mu=1,80)
1484 2148 : call wrtout([std_out, ab_out], msg)
1485 :
1486 716 : ddkfil(:)=0
1487 :
1488 : ! MGNAG: WHY THIS? all these variables are declared as optional pointers in dfpt_looppert!
1489 : ! but they are allocated here so why using pointers! Moreover OPTIONAL arguments MUST
1490 : ! be passed by keyword for better clarity and robustness!
1491 : ! People should learn how to program Fort90 before being allowed to change the code!
1492 : ! v5[26] crashes in dfpt_looppert
1493 : ! The best solution would be using a datatype to gather the results!
1494 1432 : ABI_MALLOC(clflg,(3,mpert))
1495 716 : if(dtset%ieig2rf > 0) then
1496 252 : ABI_MALLOC(eigen0_pert,(dtset%mband*dtset%nkpt*dtset%nsppol))
1497 168 : ABI_MALLOC(eigenq_pert,(dtset%mband*dtset%nkpt*dtset%nsppol))
1498 168 : ABI_MALLOC(occ_rbz_pert,(dtset%mband*dtset%nkpt*dtset%nsppol))
1499 : end if
1500 716 : if(dtset%efmas > 0) then
1501 51 : ABI_MALLOC(eigen0_pert,(dtset%mband*dtset%nkpt*dtset%nsppol))
1502 : end if
1503 : ! Note that kg, cg, eigen0, mpw and npwarr are NOT passed to dfpt_looppert :
1504 : ! they will be reinitialized for each perturbation, with an eventually
1505 : ! reduced set of k point, thanks to the use of symmetry operations.
1506 : call dfpt_looppert(atindx,blkflg,codvsn,cpus,dim_eigbrd,dim_eig2nkq,doccde,&
1507 : & ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,dyvdw,dyfr_cplex,dyfr_nondiag,&
1508 : & d2bbb,d2lo,d2nl,d2ovl,efmasdeg,efmasval,eigbrd,eig2nkq,&
1509 : & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
1510 : & etotal,fermie,iexit,indsym,kxc,&
1511 : & dtset%mkmem,mkqmem,mk1mem,mpert,mpi_enreg,my_natom,nattyp,&
1512 : & nfftf,nhat,dtset%nkpt,nkxc,dtset%nspden,dtset%nsym,occ,&
1513 : & paw_an,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
1514 : & pertsy,prtbbb,psps,rfpert,rf2_dirs_from_rfpert_nl,rhog,rhor,symq,symrec,timrev,&
1515 : & usecprj,usevdw,usevxctau,vtrial,vxc,vxcavg,vxctau,xred,clflg,occ_rbz_pert,eigen0_pert,eigenq_pert,&
1516 716 : & eigen1_pert,nkpt_rbz,eigenq_fine,hdr_fine,hdr0)
1517 :
1518 : ! #####################################################################
1519 : end if
1520 :
1521 716 : call timab(138,1,tsec)
1522 :
1523 57996 : write(msg, '(80a,a,a,a,a)' ) ('=',mu=1,80),ch10,ch10,&
1524 58712 : ' ---- first-order wavefunction calculations are completed ----',ch10
1525 2148 : call wrtout([std_out, ab_out], msg)
1526 :
1527 716 : ABI_FREE(vxc)
1528 716 : ABI_FREE(vxctau)
1529 :
1530 : ! Cleaning GEMM nonlop data
1531 716 : if(gemm_nonlop_use_gemm) then
1532 0 : call destroy_gemm_nonlop(dtset%gpu_option)
1533 0 : gemm_nonlop_use_gemm = .false.
1534 : end if
1535 :
1536 716 : if (dtset%prepanl==1.and.(rf2_dkdk/=0 .or. rf2_dkde/=0)) then
1537 14 : ABI_FREE(rfpert_nl)
1538 : end if
1539 :
1540 : !Output of the localization tensor
1541 716 : if ( rfpert(natom+1) /= 0 .and. (me == 0) .and. dtset%occopt<=2) then
1542 157 : call wrtloctens(blkflg,d2bbb,d2nl,dtset%mband,mpert,natom,dtset%prtbbb,rprimd,psps%usepaw)
1543 : end if
1544 :
1545 : !The perturbation natom+1 was only an auxiliary perturbation,
1546 : !needed to construct the electric field response, so its flag is now set to 0.
1547 : !rfpert(natom+1)=0
1548 :
1549 : !Were 2DTE computed ?
1550 716 : if(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0 .or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 .and. rfmagn==0)then
1551 :
1552 227 : write(msg,'(a,a)' )ch10,' respfn : d/dk was computed, but no 2DTE, so no DDB output.'
1553 681 : call wrtout([std_out, ab_out], msg)
1554 :
1555 : ! If 2DTE were computed, only one processor must output them and compute
1556 : ! frequencies.
1557 489 : else if(me==0)then
1558 :
1559 418 : write(msg,'(a,a)' )ch10,' ==> Compute Derivative Database <== '
1560 1254 : call wrtout([std_out, ab_out], msg)
1561 :
1562 : ! In the RESPFN code, dfpt_nstdy and stady3 were called here
1563 1140800 : d2nfr(:,:,:,:,:)=d2lo(:,:,:,:,:)+d2nl(:,:,:,:,:)
1564 89800 : if (psps%usepaw==1) d2nfr(:,:,:,:,:)=d2nfr(:,:,:,:,:)+d2ovl(:,:,:,:,:)
1565 :
1566 418 : zero_by_symm=1
1567 418 : if(dtset%rfmeth<0)zero_by_symm=0
1568 :
1569 : ! In case of bbb decomposition
1570 418 : if(prtbbb==1)then
1571 6 : ABI_MALLOC(blkflg1,(3,mpert,3,mpert))
1572 6 : ABI_MALLOC(blkflg2,(3,mpert,3,mpert))
1573 2018 : blkflg2(:,:,:,:) = blkflg(:,:,:,:)
1574 20 : do ipert = 1, mpert
1575 182 : do ipert2 = 1, mpert
1576 180 : if ((ipert /= natom + 2).and.(ipert>natom).and.(ipert2/=natom+2)) then
1577 1248 : blkflg2(:,ipert2,:,ipert) = 0
1578 : end if
1579 : end do
1580 : end do
1581 6 : ABI_MALLOC(d2tmp,(2,3,mpert,3,mpert))
1582 18 : do iband = 1,dtset%mband
1583 39472 : d2tmp(:,:,:,:,:)=zero
1584 16144 : blkflg1(:,:,:,:) = blkflg2(:,:,:,:)
1585 4480 : d2tmp(:,:,natom+2,:,:) = d2bbb(:,:,:,:,iband,iband)
1586 : call d2sym3(blkflg1,d2tmp,indsym,mpert,natom,dtset%nsym,qphon,symq,&
1587 16 : & symrec,dtset%symrel,timrev,zero_by_symm)
1588 4482 : d2bbb(:,:,:,:,iband,iband) = d2tmp(:,:,natom+2,:,:)
1589 : end do
1590 2 : ABI_FREE(blkflg1)
1591 2 : ABI_FREE(blkflg2)
1592 2 : ABI_FREE(d2tmp)
1593 : end if
1594 :
1595 : ! Complete the d2nfr matrix by symmetrization of the existing elements
1596 : !write(std_out,*)"blkflg before d2sym3: ", blkflg
1597 : call d2sym3(blkflg,d2nfr,indsym,mpert,natom,dtset%nsym,qphon,symq,symrec,&
1598 418 : & dtset%symrel,timrev,zero_by_symm,eta=dtset%rfeta)
1599 : !write(std_out,*)"blkflg after d2sym3: ", blkflg
1600 :
1601 418 : if(rfphon==1.and.psps%n1xccc/=0)then
1602 : ! Complete the dyfrx1 matrix by symmetrization of the existing elements
1603 : call d2sym3(blkflgfrx1,dyfrx1,indsym,natom,natom,dtset%nsym,qphon,symq,&
1604 239 : & symrec,dtset%symrel,timrev,zero_by_symm,eta=dtset%rfeta)
1605 : end if
1606 :
1607 : ! Note that d2sym3 usually complete the 2nd-order matrix
1608 : ! with elements that are zero by symmetry, automatically,
1609 : ! unless it has been explicitly asked not to do so.
1610 : ! blkflg is then set to 1 for these matrix elements, even if there has be no calculation.
1611 :
1612 : ! Add the frozen-wf (dyfrwf) part to the ewald part (dyew),
1613 : ! the part 1 of the frozen wf part of the xc core correction
1614 : ! (dyfrx1) and the non-frozen part (dynfr) to get the second-order
1615 : ! derivative matrix (d2matr), then
1616 : ! take account of the non-cartesian coordinates (d2cart).
1617 1254 : ABI_MALLOC(d2cart,(2,3,mpert,3,mpert))
1618 1254 : ABI_MALLOC(carflg,(3,mpert,3,mpert))
1619 1254 : ABI_MALLOC(d2matr,(2,3,mpert,3,mpert))
1620 418 : outd2=1
1621 : call dfpt_gatherdy(asr,becfrnl,dtset%berryopt,blkflg,carflg,&
1622 : & chneut,dyew,dyfrwf,dyfrx1,dyfr_cplex,dyfr_nondiag,dyvdw,d2bbb,d2cart,d2cart_bbb,d2matr,d2nfr,&
1623 : & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
1624 : & gprimd,dtset%mband,mpert,natom,ntypat,outd2,pawbec,pawpiezo,piezofrnl,dtset%prtbbb,&
1625 418 : & rfpert,rprimd,dtset%typat,ucvol,usevdw,psps%ziontypat)
1626 :
1627 : ! Output of the dynamical matrix (master only)
1628 : call dfpt_dyout(becfrnl,dtset%berryopt,blkflg,carflg,ddkfil,dyew,dyfrlo,&
1629 : & dyfrnl,dyfrx1,dyfrx2,dyfr_cplex,dyfr_nondiag,dyvdw,d2cart,d2cart_bbb,d2eig0,&
1630 : & d2k0,d2lo,d2loc0,d2matr,d2nl,d2nl0,d2nl1,d2ovl,d2vn,&
1631 : & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
1632 : & has_full_piezo,has_allddk,ab_out,dtset%mband,mpert,natom,ntypat,&
1633 : & outd2,pawbec,pawpiezo,piezofrnl,dtset%prtbbb,dtset%prtvol,qzero,&
1634 418 : & dtset%typat,rfdir,rfmagn,rfpert,rfphon,rfstrs,psps%usepaw,usevdw,psps%ziontypat)
1635 :
1636 :
1637 : ! Initialize ddb header object
1638 : call ddb_hdr%init(dtset,psps,pawtab, dscrpt=' Note : temporary (transfer) database ', &
1639 418 : nblok=1,xred=xred,occ=occ,ngfft=ngfft)
1640 :
1641 : ! Initialize ddb object
1642 418 : call ddb%init(dtset, nblok=1, mpert=mpert, with_d2E=.true.)
1643 :
1644 : ! Set the values for the 2nd order derivatives
1645 418 : call ddb%set_omega(iblok=1, omega=dtset%rfomega)
1646 418 : call ddb%set_qpt(iblok=1, qpt=qphon(1:3))
1647 418 : call ddb%set_d2matr(1, d2matr, blkflg)
1648 :
1649 : ! Output dynamical matrix
1650 418 : call ddb%write(ddb_hdr, dtfil%fnameabo_ddb)
1651 :
1652 : ! Deallocate ddb object
1653 418 : call ddb_hdr%free()
1654 418 : call ddb%free()
1655 :
1656 : ! In case of phonons, diagonalize the dynamical matrix
1657 418 : if(rfphon==1)then
1658 :
1659 : ! First, suppress the 'wings' elements,
1660 : ! for which the diagonal element is not known
1661 382 : call wings3(carflg,d2cart,mpert)
1662 :
1663 : ! Check the analyticity of the dynamical matrix
1664 382 : analyt=0
1665 382 : if (rfpert(natom+2)==0 .or. rfpert(natom+2)==2 .or. sumg0==1 ) analyt=1
1666 :
1667 : ! Diagonalize the analytic part
1668 1146 : ABI_MALLOC(displ,(2*3*natom*3*natom))
1669 1146 : ABI_MALLOC(eigval,(3*natom))
1670 764 : ABI_MALLOC(eigvec,(2*3*natom*3*natom))
1671 764 : ABI_MALLOC(phfrq,(3*natom))
1672 382 : qphnrm=one
1673 : call dfpt_phfrq(dtset%amu_orig(:,1),displ,d2cart,eigval,eigvec,indsym,mpert,&
1674 : & dtset%nsym,natom,dtset%nsym,ntypat,phfrq,qphnrm,qphon,&
1675 382 : & dtset%rprimd_orig(1:3,1:3,1),0,dtset%symrel,dtset%symafm,dtset%typat,ucvol)
1676 :
1677 : ! Print the phonon frequencies
1678 382 : call dfpt_prtph(displ,0,dtset%enunit,ab_out,natom,phfrq,qphnrm,qphon)
1679 :
1680 : ! Check the completeness of the dynamical matrix and eventually send a warning
1681 382 : call chkph3(carflg,0,mpert,natom)
1682 : end if ! end case of phonons
1683 : end if !end me == 0
1684 :
1685 : !Compute the other terms for AHC dynamic and AHC full
1686 716 : if (.not.(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0.or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 &
1687 : & .and. rfmagn==0)) then
1688 489 : if(rfphon==1) then ! AHC can only be computed in case of phonons
1689 :
1690 : ! Stuff for parallelism
1691 450 : if(master /= me) then
1692 204 : ABI_MALLOC(phfrq,(3*natom))
1693 204 : ABI_MALLOC(displ,(2*3*natom*3*natom))
1694 : end if
1695 450 : call xmpi_bcast (phfrq,master,mpi_enreg%comm_cell,ierr) !Broadcast phfrq and displ
1696 450 : call xmpi_bcast (displ,master,mpi_enreg%comm_cell,ierr) !to all processes
1697 :
1698 450 : if(dtset%ieig2rf == 3 .or. dtset%ieig2rf == 4 .or. dtset%ieig2rf == 5 ) then
1699 54 : bdeigrf = dtset%bdeigrf
1700 54 : if(dtset%bdeigrf == -1) bdeigrf = dtset%mband
1701 : ! if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
1702 : ! & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then
1703 : ! write(std_out,*)'Reading the dense grid WF file'
1704 : ! call wfk_read_eigenvalues(dtfil%fnameabi_wfkfine,eigenq_fine,hdr_fine,mpi_enreg%comm_world)
1705 : ! ABI_CHECK(SIZE(eigenq_fine,DIM=1)==Dtset%mband,"Size eigenq_fine != mband")
1706 : ! endif
1707 54 : if(dtset%kptopt==3 .or. dtset%kptopt==0 .or. dtset%nsym==1)then
1708 54 : write(std_out,*) 'Entering: eig2tot'
1709 54 : if(dtset%smdelta>0)then
1710 54 : if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
1711 : & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then
1712 : call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,&
1713 : & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,&
1714 : & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,&
1715 0 : & occ_rbz_pert,hdr0,eigbrd,eigenq_fine,hdr_fine)
1716 : else
1717 : call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,&
1718 : & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,&
1719 : & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,&
1720 54 : & occ_rbz_pert,hdr0,eigbrd)
1721 : end if
1722 : else
1723 0 : if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
1724 : & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then
1725 : call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,&
1726 : & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,&
1727 : & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,&
1728 0 : & occ_rbz_pert,hdr0,eigbrd,eigenq_fine,hdr_fine)
1729 : else
1730 : call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,&
1731 : & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,&
1732 : & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,&
1733 0 : & occ_rbz_pert,hdr0)
1734 : end if
1735 : end if
1736 54 : write(std_out,*) 'Leaving: eig2tot'
1737 : end if
1738 : end if
1739 450 : if (dtset%ieig2rf > 0) then
1740 84 : ABI_FREE(eigen0_pert)
1741 84 : ABI_FREE(eigenq_pert)
1742 84 : ABI_FREE(occ_rbz_pert)
1743 84 : ABI_FREE(eigen1_pert)
1744 84 : call hdr0%free()
1745 84 : if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
1746 : & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then
1747 3 : call hdr_fine%free()
1748 3 : ABI_FREE(eigenq_fine)
1749 : end if
1750 : end if ! ieig2rf == 3 or %ieig2rf == 4 or %ieig2rf == 5
1751 : end if ! rfphon==1
1752 : end if
1753 716 : if(dtset%efmas>0) then
1754 17 : ABI_FREE(eigen0_pert)
1755 17 : ABI_FREE(eigen1_pert)
1756 : end if
1757 :
1758 716 : ABI_FREE(doccde)
1759 :
1760 716 : if(me==0)then
1761 622 : if (.not.(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0 .or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 &
1762 : & .and. rfmagn==0) )then
1763 418 : if(rfphon==1)then
1764 : ! Compute and print the T=0 Fan, and possibly DDW contributions to the eigenenergies.
1765 382 : if(dtset%ieig2rf > 0) then
1766 5508 : write(msg, '(80a,9a)' ) ('=',mu=1,80),ch10,ch10,&
1767 68 : & ' ---- T=0 shift of eigenenergies due to electron-phonon interation at q ---- ',ch10,&
1768 68 : & ' Warning : the total shift must be computed through anaddb, ',ch10,&
1769 68 : & ' here, only the contribution of one q point is printed. ',ch10,&
1770 5576 : & ' Print first the electronic eigenvalues, then the q-dependent Fan shift of eigenvalues.'
1771 204 : call wrtout([std_out, ab_out], msg)
1772 :
1773 68 : if(qeq0)then
1774 21 : write(msg, '(a)' )' Phonons at gamma, also compute the Diagonal Debye-Waller shift of eigenvalues.'
1775 63 : call wrtout([std_out, ab_out], msg)
1776 : end if
1777 :
1778 68 : write(msg, '(a)' ) ' '
1779 204 : call wrtout([std_out, ab_out], msg)
1780 :
1781 : call prteigrs(eigen0,dtset%enunit,fermie,fermih,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,&
1782 : & dtset%mband,dtset%nband,dtset%nbdbuf,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,3,0,dtset%prtvol,&
1783 68 : & eigen0,zero,zero,dtset%wtk)
1784 :
1785 68 : write(msg, '(a)' ) ch10
1786 204 : call wrtout([std_out, ab_out], msg)
1787 :
1788 : ! Compute and print Fan contribution
1789 204 : ABI_MALLOC(eigen_fan,(dtset%mband*dtset%nkpt*dtset%nsppol))
1790 136 : ABI_MALLOC(eigen_fan_mean,(dtset%mband*dtset%nkpt*dtset%nsppol))
1791 : call elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_fan,gprimd,&
1792 68 : & dtset%mband,natom,dtset%nkpt,dtset%nsppol,1,phfrq,dtset%prtvol)
1793 68 : call eigen_meandege(eigen0,eigen_fan,eigen_fan_mean,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,2)
1794 : call prteigrs(eigen_fan_mean,dtset%enunit,fermie,fermih,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,&
1795 : & dtset%mband,dtset%nband,dtset%nbdbuf,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,5,0,dtset%prtvol,&
1796 68 : & eigen0,zero,zero,dtset%wtk)
1797 :
1798 68 : if(qeq0 .or. dtset%getgam_eig2nkq>0)then
1799 :
1800 29 : write(msg, '(a)' ) ch10
1801 87 : call wrtout([std_out, ab_out], msg)
1802 :
1803 : ! Compute and print Diagonal Debye-Waller contribution
1804 58 : ABI_MALLOC(eigen_ddw,(dtset%mband*dtset%nkpt*dtset%nsppol))
1805 58 : ABI_MALLOC(eigen_ddw_mean,(dtset%mband*dtset%nkpt*dtset%nsppol))
1806 29 : if(qeq0)then
1807 : call elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_ddw,gprimd,&
1808 21 : & dtset%mband,natom,dtset%nkpt,dtset%nsppol,2,phfrq,dtset%prtvol)
1809 21 : if(results_respfn%gam_jdtset == -dtset%jdtset)then
1810 2 : sz1=dtset%mband*dtset%nsppol
1811 2 : sz2=natom*dim_eig2nkq
1812 14 : ABI_MALLOC(results_respfn%gam_eig2nkq,(2,sz1,dtset%nkpt,3,natom,3,sz2))
1813 8210 : results_respfn%gam_eig2nkq(:,:,:,:,:,:,:)=eig2nkq(:,:,:,:,:,:,:)
1814 2 : results_respfn%gam_jdtset=dtset%jdtset
1815 : end if
1816 8 : else if(dtset%getgam_eig2nkq>0)then
1817 8 : if(results_respfn%gam_jdtset==dtset%getgam_eig2nkq)then
1818 : call elph2_fanddw(dim_eig2nkq,displ,results_respfn%gam_eig2nkq,eigen_ddw,&
1819 8 : & gprimd,dtset%mband,natom,dtset%nkpt,dtset%nsppol,2,phfrq,dtset%prtvol)
1820 : else
1821 : write(msg,'(a,i0,2a,i0,2a)')&
1822 0 : & 'results_respfn%gam_jdtset=',results_respfn%gam_jdtset,ch10,&
1823 0 : & 'dtset%getgam_eig2nkq=',dtset%getgam_eig2nkq,ch10,&
1824 0 : & 'So, it seems eig2nkq at gamma has not yet been computed, while it is needed now.'
1825 0 : ABI_BUG(msg)
1826 : end if
1827 : end if
1828 29 : call eigen_meandege(eigen0,eigen_ddw,eigen_ddw_mean,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,2)
1829 : call prteigrs(eigen_ddw_mean,dtset%enunit,fermie,fermih,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,&
1830 : & dtset%mband,dtset%nband,dtset%nbdbuf,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,6,0,dtset%prtvol,&
1831 29 : & eigen0,zero,zero,dtset%wtk)
1832 29 : write(msg, '(a)' ) ch10
1833 87 : call wrtout([std_out, ab_out], msg)
1834 :
1835 : ! Print sum of mean Fan and DDW
1836 87 : ABI_MALLOC(eigen_fanddw,(dtset%mband*dtset%nkpt*dtset%nsppol))
1837 3338 : eigen_fanddw=eigen_fan_mean+eigen_ddw_mean
1838 : call prteigrs(eigen_fanddw,dtset%enunit,fermie,fermih,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,&
1839 : & dtset%mband,dtset%nband,dtset%nbdbuf,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,7,0,dtset%prtvol,&
1840 29 : & eigen0,zero,zero,dtset%wtk)
1841 :
1842 29 : ABI_FREE(eigen_ddw)
1843 29 : ABI_FREE(eigen_ddw_mean)
1844 29 : ABI_FREE(eigen_fanddw)
1845 :
1846 : end if
1847 :
1848 68 : ABI_FREE(eigen_fan)
1849 68 : ABI_FREE(eigen_fan_mean)
1850 : end if
1851 :
1852 : ! In case of a non-analytical part,
1853 : ! get the phonon frequencies for three different directions (in cartesian coordinates)
1854 382 : if(analyt==0)then
1855 76 : qphnrm=zero
1856 304 : do idir=1,3
1857 : ! Need to know the corresponding dielectric constant
1858 304 : if(carflg(idir,natom+2,idir,natom+2)==1)then
1859 219 : qphon(:)=zero ; qphon(idir)=one
1860 : ! Get the phonon frequencies
1861 : call dfpt_phfrq(dtset%amu_orig(:,1),displ,d2cart,eigval,eigvec,indsym,mpert,&
1862 : & dtset%nsym,natom,dtset%nsym,ntypat,phfrq,qphnrm,qphon,&
1863 219 : & dtset%rprimd_orig(1:3,1:3,1),0,dtset%symrel,dtset%symafm,dtset%typat,ucvol)
1864 : ! Print the phonon frequencies
1865 219 : call dfpt_prtph(displ,0,dtset%enunit,ab_out,natom,phfrq,qphnrm,qphon)
1866 : ! Check the completeness of the dynamical matrix
1867 : ! and eventually send a warning
1868 219 : call chkph3(carflg,idir,mpert,natom)
1869 : end if
1870 : end do
1871 : if (idir < 4) then
1872 : qphon(idir)=zero
1873 : end if
1874 : end if
1875 :
1876 382 : ABI_FREE(displ)
1877 382 : ABI_FREE(eigval)
1878 382 : ABI_FREE(eigvec)
1879 382 : ABI_FREE(phfrq)
1880 : end if ! rfphon == 1
1881 418 : ABI_FREE(carflg)
1882 418 : ABI_FREE(d2cart)
1883 418 : ABI_FREE(d2matr)
1884 : end if ! End condition on if.not.
1885 : end if ! master node
1886 :
1887 : !Deallocate arrays
1888 716 : ABI_SFREE(displ)
1889 716 : ABI_SFREE(eigval)
1890 716 : ABI_SFREE(eigvec)
1891 716 : ABI_SFREE(phfrq)
1892 :
1893 716 : ABI_FREE(clflg)
1894 716 : ABI_FREE(atindx)
1895 716 : ABI_FREE(atindx1)
1896 716 : ABI_FREE(blkflg)
1897 716 : ABI_FREE(dyew)
1898 716 : ABI_FREE(dyewq0)
1899 716 : ABI_FREE(dyfrlo)
1900 716 : ABI_FREE(dyfrnl)
1901 716 : ABI_FREE(dyfrwf)
1902 716 : ABI_FREE(dyfrx1)
1903 716 : ABI_FREE(dyfrx2)
1904 716 : ABI_FREE(dyvdw)
1905 716 : ABI_FREE(d2bbb)
1906 716 : ABI_FREE(d2cart_bbb)
1907 716 : ABI_FREE(d2eig0)
1908 716 : ABI_FREE(d2k0)
1909 716 : ABI_FREE(d2lo)
1910 716 : ABI_FREE(d2loc0)
1911 716 : ABI_FREE(d2nfr)
1912 716 : ABI_FREE(d2nl)
1913 716 : ABI_FREE(d2nl0)
1914 716 : ABI_FREE(d2nl1)
1915 716 : ABI_FREE(d2ovl)
1916 716 : ABI_FREE(d2vn)
1917 716 : ABI_FREE(eigen0)
1918 716 : ABI_FREE(eig2nkq)
1919 716 : ABI_FREE(eigbrd)
1920 716 : ABI_FREE(eltcore)
1921 716 : ABI_FREE(elteew)
1922 716 : ABI_FREE(eltfrhar)
1923 716 : ABI_FREE(eltfrnl)
1924 716 : ABI_FREE(eltfrloc)
1925 716 : ABI_FREE(eltfrkin)
1926 716 : ABI_FREE(eltfrxc)
1927 716 : ABI_FREE(eltvdw)
1928 716 : ABI_FREE(becfrnl)
1929 716 : ABI_FREE(piezofrnl)
1930 716 : call efmasdeg_free_array(efmasdeg)
1931 716 : call efmasval_free_array(efmasval)
1932 716 : ABI_FREE(grxc)
1933 716 : ABI_FREE(indsym)
1934 716 : ABI_FREE(kxc)
1935 716 : ABI_FREE(nattyp)
1936 716 : ABI_FREE(pertsy)
1937 716 : ABI_FREE(rfpert)
1938 716 : ABI_FREE(rhog)
1939 716 : ABI_FREE(rhor)
1940 716 : ABI_FREE(taug)
1941 716 : ABI_FREE(taur)
1942 716 : ABI_FREE(symq)
1943 716 : ABI_FREE(symrec)
1944 716 : ABI_FREE(symrel_cart)
1945 716 : ABI_FREE(vtrial)
1946 716 : ABI_FREE(ylm)
1947 716 : ABI_FREE(ylmgr)
1948 716 : call pawfgr_destroy(pawfgr)
1949 716 : if (psps%usepaw==1) then
1950 112 : call pawrhoij_free(pawrhoij)
1951 112 : call paw_an_free(paw_an)
1952 112 : call paw_ij_free(paw_ij)
1953 112 : call pawfgrtab_free(pawfgrtab)
1954 : end if
1955 716 : ABI_FREE(nhat)
1956 716 : ABI_FREE(nhatgr)
1957 904 : ABI_FREE(pawrhoij)
1958 904 : ABI_FREE(paw_an)
1959 904 : ABI_FREE(paw_ij)
1960 904 : ABI_FREE(pawfgrtab)
1961 716 : if(rfphon==1.and.psps%n1xccc/=0)then
1962 287 : ABI_FREE(blkflgfrx1)
1963 : end if
1964 :
1965 : ! Clean the header
1966 716 : call hdr%free()
1967 :
1968 : !Clean GPU data
1969 : #if defined HAVE_GPU
1970 : if (dtset%gpu_option/=ABI_GPU_DISABLED) then
1971 : call dealloc_hamilt_gpu(0,dtset%gpu_option)
1972 : end if
1973 : #endif
1974 :
1975 : ABI_NVTX_END_RANGE()
1976 716 : call timab(138,2,tsec)
1977 716 : call timab(132,2,tsec)
1978 :
1979 : DBG_EXIT("COLL")
1980 :
1981 5012 : end subroutine respfn
1982 : !!***
1983 :
1984 : !!****f* m_respfn_driver/wrtloctens
1985 : !! NAME
1986 : !! wrtloctens
1987 : !!
1988 : !! FUNCTION
1989 : !! Output of the localisation tensor
1990 : !!
1991 : !! INPUTS
1992 : !! blkflg = flags for each element of the 2DTE (=1 if computed)
1993 : !! d2bbb = band by band decomposition of second order derivatives
1994 : !! d2nl = non-local contributions to the 2DTEs
1995 : !! mband = maximum number of bands
1996 : !! mpert = maximum number of ipert
1997 : !! natom = number of atoms in unit cell
1998 : !! prtbbb = if = 1, write the band by band decomposition of the localization tensor
1999 : !! rprimd = dimensional primitive translations for real space (bohr)
2000 : !! usepaw = flag for PAW
2001 : !!
2002 : !! OUTPUT
2003 : !! (only writing)
2004 : !!
2005 : !! TODO
2006 : !! The localization tensor cannot be defined in the metallic case. It should not be computed.
2007 : !!
2008 : !! SOURCE
2009 :
2010 157 : subroutine wrtloctens(blkflg,d2bbb,d2nl,mband,mpert,natom,prtbbb,rprimd,usepaw)
2011 :
2012 : !Arguments ------------------------------------
2013 : !scalars
2014 : integer,intent(in) :: mband,mpert,natom,prtbbb,usepaw
2015 : !arrays
2016 : integer,intent(in) :: blkflg(3,mpert,3,mpert)
2017 : real(dp),intent(in) :: rprimd(3,3)
2018 : real(dp),intent(inout) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb)
2019 : real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert)
2020 :
2021 : !Local variables ------------------------------
2022 : !scalars
2023 : integer :: flag,iband,idir,idir2,jband
2024 : character(len=500) :: msg
2025 : !arrays
2026 : real(dp) :: loctenscart(2,3,3)
2027 157 : real(dp),allocatable :: loctenscart_bbb(:,:,:,:,:)
2028 :
2029 : ! *********************************************************************
2030 :
2031 : !This feature is disabled in the PAW case
2032 157 : if (usepaw==1) return
2033 :
2034 119 : if(prtbbb==1)then
2035 8 : ABI_MALLOC(loctenscart_bbb,(2,3,3,mband,mband*prtbbb))
2036 3986 : loctenscart_bbb(:,:,:,:,:)=zero
2037 : end if
2038 :
2039 : !complete missing elements
2040 119 : flag = 0
2041 476 : do idir2 = 1,3
2042 1547 : do idir = 1,3
2043 :
2044 1428 : if (blkflg(idir2,natom+1,idir,natom+1)==0) then
2045 627 : if (blkflg(idir,natom+1,idir2,natom+1)==0) then
2046 : flag = 1
2047 : else
2048 195 : d2nl(1,idir2,natom+1,idir,natom+1) = d2nl(1,idir,natom+1,idir2,natom+1)
2049 195 : d2nl(2,idir2,natom+1,idir,natom+1) =-d2nl(2,idir,natom+1,idir2,natom+1)
2050 195 : if(prtbbb==1)then
2051 438 : d2bbb(1,idir2,idir,natom+1,:,:) = d2bbb(1,idir,idir2,natom+1,:,:)
2052 438 : d2bbb(2,idir2,idir,natom+1,:,:) =-d2bbb(2,idir,idir2,natom+1,:,:)
2053 : end if
2054 :
2055 : end if
2056 : end if
2057 :
2058 : end do ! idir=1,3
2059 : end do ! idir2=1,3
2060 :
2061 : !Transform the tensor to cartesian coordinates
2062 :
2063 6188 : loctenscart(1,:,:) = matmul(rprimd,d2nl(1,:,natom+1,:,natom+1))
2064 6188 : loctenscart(2,:,:) = matmul(rprimd,d2nl(2,:,natom+1,:,natom+1))
2065 :
2066 7616 : loctenscart(1,:,:) = matmul(loctenscart(1,:,:),transpose(rprimd))
2067 7616 : loctenscart(2,:,:) = matmul(loctenscart(2,:,:),transpose(rprimd))
2068 :
2069 3689 : loctenscart(:,:,:) = loctenscart(:,:,:)/(two_pi**2)
2070 :
2071 119 : if (prtbbb == 1) then
2072 :
2073 2 : write(msg,'(a,a)')ch10, ' Band by band decomposition of the localisation tensor (bohr^2)'
2074 6 : call wrtout([std_out, ab_out], msg)
2075 :
2076 18 : do iband = 1,mband
2077 146 : do jband = 1,mband
2078 :
2079 6656 : loctenscart_bbb(1,:,:,iband,jband) = matmul(rprimd,d2bbb(1,:,:,natom+1,iband,jband))
2080 6656 : loctenscart_bbb(2,:,:,iband,jband) = matmul(rprimd,d2bbb(2,:,:,natom+1,iband,jband))
2081 :
2082 8192 : loctenscart_bbb(1,:,:,iband,jband) = matmul(loctenscart_bbb(1,:,:,iband,jband),transpose(rprimd))
2083 8192 : loctenscart_bbb(2,:,:,iband,jband) = matmul(loctenscart_bbb(2,:,:,iband,jband),transpose(rprimd))
2084 :
2085 3968 : loctenscart_bbb(:,:,:,iband,jband) = loctenscart_bbb(:,:,:,iband,jband)/(two_pi**2)
2086 :
2087 :
2088 128 : write(msg,'(a,a,i5,a,i5,a)')ch10, &
2089 128 : ' Localisation tensor (bohr^2) for band ',iband,',',jband, &
2090 256 : ' in cartesian coordinates'
2091 384 : call wrtout([std_out, ab_out], msg)
2092 :
2093 128 : write(ab_out,*)' direction matrix element'
2094 128 : write(ab_out,*)' alpha beta real part imaginary part'
2095 528 : do idir2 = 1,3
2096 1664 : do idir = 1,3
2097 1152 : write(ab_out,'(5x,i1,8x,i1,3x,2f16.10)')idir2,idir,&
2098 1152 : loctenscart_bbb(1,idir2,idir,iband,jband),&
2099 2688 : loctenscart_bbb(2,idir2,idir,iband,jband)
2100 : end do
2101 : end do
2102 :
2103 : end do !jband
2104 : end do !iband
2105 :
2106 : end if !prtbbb
2107 :
2108 119 : if (usepaw==0) then
2109 119 : write(msg,'(a,a,a,a)')ch10, &
2110 119 : & ' Total localisation tensor (bohr^2) in cartesian coordinates',ch10,&
2111 238 : & ' WARNING : still subject to testing - especially symmetries.'
2112 : else
2113 0 : write(msg,'(a,a,a,a,a,a)')ch10, &
2114 0 : & ' Total localisation tensor (bohr^2) in cartesian coordinates',ch10,&
2115 0 : & ' WARNING : probably wrong for PAW (printing for testing purpose)',ch10,&
2116 0 : & ' WARNING : still subject to testing - especially symmetries.'
2117 : end if
2118 357 : call wrtout([std_out, ab_out], msg)
2119 :
2120 119 : write(ab_out,*)' direction matrix element'
2121 119 : write(ab_out,*)' alpha beta real part imaginary part'
2122 476 : do idir2 = 1,3
2123 1547 : do idir = 1,3
2124 1071 : write(ab_out,'(5x,i1,8x,i1,3x,2f16.10)')idir2,idir,&
2125 1071 : & loctenscart(1,idir2,idir),&
2126 2499 : & loctenscart(2,idir2,idir)
2127 : end do
2128 : end do
2129 :
2130 119 : if (flag == 1) then
2131 54 : write(msg,'(6a)')ch10,&
2132 54 : & ' WARNING : Localization tensor calculation (this does not apply to other properties).',ch10,&
2133 54 : & ' Not all d/dk perturbations were computed. So the localization tensor in reciprocal space is incomplete,',ch10,&
2134 108 : & ' and transformation to cartesian coordinates may be wrong. Check input variable rfdir.'
2135 162 : call wrtout([std_out, ab_out], msg)
2136 : end if
2137 :
2138 119 : ABI_SFREE(loctenscart_bbb)
2139 :
2140 : end subroutine wrtloctens
2141 : !!***
2142 :
2143 : !!****f* ABINIT/dfpt_dyout
2144 : !! NAME
2145 : !! dfpt_dyout
2146 : !!
2147 : !! FUNCTION
2148 : !! Output of all quantities related to the 2nd-order matrix:
2149 : !! Ewald part, local and non-local frozen wf part, core contributions,
2150 : !! local and non-local variational part, 2nd-order matrix itself, and, for the phonon part,
2151 : !! eigenfrequencies, in Hartree, meV and cm-1.
2152 : !! Also output unformatted 2nd-order matrix for later use in the Brillouin-zone interpolation
2153 : !!
2154 : !! INPUTS
2155 : !! becfrnl(3,natom,3*pawbec)=NL frozen contribution to Born Effective Charges (PAW only)
2156 : !! blkflg(3,mpert,3,mpert)= ( 1 if the element of the dynamical
2157 : !! matrix has been calculated ; 0 otherwise )
2158 : !! carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
2159 : !! 2DTE matrix has been calculated correctly ; 0 otherwise )
2160 : !! ddkfil(3)=components are 1 if corresponding d/dk file exists, otherwise 0
2161 : !! (in what follows, DYMX means dynamical matrix, and D2MX means 2nd-order matrix)
2162 : !! dyew(2,3,natom,3,natom)=Ewald part of the DYMX
2163 : !! dyfrlo(3,3,natom)=frozen wf local part of the DYMX
2164 : !! dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=frozen wf nonloc part of the DYMX
2165 : !! dyfrx1(2,3,natom,3,natom)=frozen wf nonlin. xc core corr.(1)
2166 : !! part of the DYMX
2167 : !! dyfrx2(3,3,natom)=frozen wf nonlin. xc core corr.(2) part of the DYMX
2168 : !! dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
2169 : !! dyfr_nondiag=1 if dyfrnl is non diagonal with respect to atoms; 0 otherwise
2170 : !! dyvdw(2,3,natom,3,natom*usevdw)=vdw DFT-D part of the dynamical matrix
2171 : !! d2cart(2,3,mpert,3,mpert)=D2MX in cartesian coordinates
2172 : !! d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)=
2173 : !! band by band decomposition of Born effective charges
2174 : !! (calculated from phonon-type perturbation) in cartesian coordinates
2175 : !! d2eig0(2,3,mpert,3,mpert)=0-order eigen. station. part of the D2MX
2176 : !! d2k0(2,3,mpert,3,mpert)=0-order kinet. station. part of the D2MX
2177 : !! d2lo(2,3,mpert,3,mpert)=nonstation. local part of the D2MX
2178 : !! d2loc0(2,3,mpert,3,mpert)=0-order loc station. part of the D2MX
2179 : !! d2matr(2,3,mpert,3,mpert)=D2MX in non-cartesian coordinates
2180 : !! d2nl(2,3,mpert,3,mpert)=nonstation. nonloc part of the D2MX
2181 : !! d2nl0(2,3,mpert,3,mpert)=0-order nonloc station. part of the D2MX
2182 : !! d2nl1(2,3,mpert,3,mpert)=1-order nonloc station. part of the D2MX
2183 : !! d2ovl(2,mpert,3,mpert*usepaw)=1st-order change of WF overlap contributions to the 2DTEs (PAW)
2184 : !! d2vn(2,3,mpert,3,mpert)=potential*dens station. part of the D2MX and without masses included)
2185 : !! eltcore(6,6)=core contribution to the elastic tensor
2186 : !! elteew(6+3*natom,6)=Ewald contribution to the elastic tsenor
2187 : !! eltfrhar(6,6)=hartree contribution to the elastic tensor
2188 : !! eltfrkin(6,6)=kinetic contribution to the elastic tensor
2189 : !! eltfrloc(6+3*natom,6)=local psp contribution to the elastic tensor
2190 : !! eltfrnl(6+3*natom,6)=non-local psp contribution to the elastic tensor
2191 : !! eltfrxc(6+3*natom,6)=exchange-correlation contribution to the elastic tensor
2192 : !! eltvdw(6+3*natom,6*usevdw)=vdw DFT-D part of the elastic tensor
2193 : !! has_full_piezo=the full calculation of the piezoelectric tensor from electric field perturbation
2194 : !! is only available if nsym=1 (strain perturbation is not symmetrized)
2195 : !! has_allddk= True if all ddk file are present on disk, so the effective charge or piezzo
2196 : !! electric tensor are correctly computed (PAW ONLY)
2197 : !! iout=unit number for long write-up
2198 : !! mband=maximum number of bands
2199 : !! mpert =maximum number of ipert
2200 : !! natom=number of atoms
2201 : !! ntypat=number of atom types
2202 : !! outd2=option for the output of the 2nd-order matrix :
2203 : !! if outd2=1, non-stationary part
2204 : !! if outd2=2, stationary part.
2205 : !! pawbec= flag for the computation of frozen part of Born Effective Charges (PAW only)
2206 : !! pawpiezo= flag for the computation of frozen part of Piezoelectric tensor (PAW only)
2207 : !! prtbbb=if 1, print the band-by-band decomposition
2208 : !! prtvol=print volume
2209 : !! qzero=1 if zero phonon wavevector
2210 : !! rfdir(3)=defines the directions for the perturbations
2211 : !! rfmagn= if 1 (2), there are response to uniform (local) spin fields
2212 : !! rfpert(mpert)=defines the perturbations
2213 : !! rfphon=if 1, there are phonon perturbations
2214 : !! rfstrs=if 1,2,3 there are strain perturbations
2215 : !! typat(natom)=integer label of each type of atom (1,2,...)
2216 : !! usepaw=1 if PAW, 0 otherwise
2217 : !! usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use
2218 : !! zion(ntypat)=charge corresponding to the atom type
2219 : !!
2220 : !! SIDE EFFECTS
2221 : !! d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)
2222 : !!
2223 : !! NOTES
2224 : !! This routine is called only by the processor me==0 .
2225 : !! In consequence, no use of message and wrtout routine.
2226 : !!
2227 : !! SOURCE
2228 :
2229 418 : subroutine dfpt_dyout(becfrnl,berryopt,blkflg,carflg,ddkfil,dyew,dyfrlo,dyfrnl,&
2230 418 : & dyfrx1,dyfrx2,dyfr_cplex,dyfr_nondiag,dyvdw,d2cart,d2cart_bbb,&
2231 418 : & d2eig0,d2k0,d2lo,d2loc0,d2matr,d2nl,d2nl0,d2nl1,d2ovl,d2vn,&
2232 418 : & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
2233 : & has_full_piezo,has_allddk,iout,mband,mpert,natom,ntypat,&
2234 418 : & outd2,pawbec,pawpiezo,piezofrnl,prtbbb,prtvol,qzero,typat,rfdir,&
2235 418 : & rfmagn,rfpert,rfphon,rfstrs,usepaw,usevdw,zion)
2236 :
2237 : !Arguments -------------------------------
2238 : !scalars
2239 : integer,intent(in) :: berryopt,dyfr_cplex,dyfr_nondiag,iout,mband,mpert
2240 : integer,intent(in) :: natom,ntypat,outd2,pawbec,pawpiezo,prtbbb,prtvol,qzero
2241 : integer, intent(in) :: rfmagn,rfphon,rfstrs,usepaw,usevdw
2242 : !arrays
2243 : integer,intent(in) :: blkflg(3,mpert,3,mpert),carflg(3,mpert,3,mpert)
2244 : integer,intent(in) :: ddkfil(3),rfdir(3),rfpert(mpert),typat(natom)
2245 : real(dp),intent(in) :: becfrnl(3,natom,3*pawbec)
2246 : real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert),d2eig0(2,3,mpert,3,mpert)
2247 : real(dp),intent(in) :: d2k0(2,3,mpert,3,mpert),d2lo(2,3,mpert,3,mpert)
2248 : real(dp),intent(in) :: d2loc0(2,3,mpert,3,mpert),d2matr(2,3,mpert,3,mpert)
2249 : real(dp),intent(in) :: d2nl(2,3,mpert,3,mpert),d2nl0(2,3,mpert,3,mpert)
2250 : real(dp),intent(in) :: d2nl1(2,3,mpert,3,mpert),d2ovl(2,3,mpert,3,mpert*usepaw)
2251 : real(dp),intent(in) :: d2vn(2,3,mpert,3,mpert)
2252 : real(dp),intent(in) :: dyew(2,3,natom,3,natom),dyfrlo(3,3,natom)
2253 : real(dp),intent(in) :: dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)
2254 : real(dp),intent(in) :: dyfrx1(2,3,natom,3,natom),dyfrx2(3,3,natom)
2255 : real(dp),intent(in) :: dyvdw(2,3,natom,3,natom*usevdw)
2256 : real(dp),intent(in) :: eltcore(6,6),elteew(6+3*natom,6)
2257 : real(dp),intent(in) :: eltfrhar(6,6),eltfrkin(6,6),eltfrloc(6+3*natom,6)
2258 : real(dp),intent(in) :: eltfrnl(6+3*natom,6),eltfrxc(6+3*natom,6)
2259 : real(dp),intent(in) :: eltvdw(6+3*natom,6*usevdw),piezofrnl(6,3*pawpiezo)
2260 : real(dp),intent(in) :: zion(ntypat)
2261 : real(dp),intent(inout) :: d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)
2262 : logical,intent(in) :: has_allddk,has_full_piezo
2263 :
2264 : !Local variables -------------------------
2265 : !scalars
2266 : integer :: iband,idir1,idir2,ii,ipert1,ipert2,jj,nelmts,nline
2267 : real(dp) :: zi,zr
2268 : !arrays
2269 : real(dp) :: delta(3,3)
2270 :
2271 : ! *********************************************************************
2272 :
2273 : ! GA: As much as I can tell, the option outd2 is always set to 1. This variable should be removed
2274 :
2275 : !Long print : includes detail of every part of the 2nd-order energy
2276 418 : if(prtvol>=10)then
2277 :
2278 : ! In case of phonon
2279 33 : if (rfphon==1)then
2280 :
2281 : ! write the Ewald part of the dynamical matrix
2282 28 : write(iout,*)' '
2283 28 : write(iout,*)' Ewald part of the dynamical matrix'
2284 28 : write(iout,*)' j1 j2 matrix element'
2285 28 : write(iout,*)' dir pert dir pert real part imaginary part'
2286 95 : do ipert1=1,natom
2287 296 : do idir1=1,3
2288 268 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2289 201 : write(iout,*)' '
2290 762 : do ipert2=1,natom
2291 2445 : do idir2=1,3
2292 2244 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2293 1395 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2294 1395 : dyew(1,idir1,ipert1,idir2,ipert2),&
2295 2790 : dyew(2,idir1,ipert1,idir2,ipert2)
2296 : end if
2297 : end do
2298 : end do
2299 : end if
2300 : end do
2301 : end do
2302 :
2303 : ! Now the local frozen wf part
2304 28 : write(iout,*)' '
2305 28 : write(iout,*)' Frozen wf local part of the dynamical matrix'
2306 28 : write(iout,*)' j1 j2 matrix element'
2307 28 : write(iout,*)' dir pert dir pert real part imaginary part'
2308 95 : do ipert1=1,natom
2309 296 : do idir1=1,3
2310 268 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2311 201 : write(iout,*)' '
2312 762 : do ipert2=1,natom
2313 2445 : do idir2=1,3
2314 2244 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2315 1395 : if(ipert1==ipert2)then
2316 459 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,dyfrlo(idir1,idir2,ipert2),zero
2317 : else
2318 936 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,zero,zero
2319 : end if
2320 : end if
2321 : end do
2322 : end do
2323 : end if
2324 : end do
2325 : end do
2326 :
2327 : ! Now the nonlo frozen wf part
2328 28 : write(iout,*)' '
2329 28 : write(iout,*)' Frozen wf non-local part of the dynamical matrix'
2330 28 : write(iout,*)' j1 j2 matrix element'
2331 28 : write(iout,*)' dir pert dir pert real part imaginary part'
2332 95 : do ipert1=1,natom
2333 296 : do idir1=1,3
2334 268 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2335 201 : write(iout,*)' '
2336 762 : do ipert2=1,natom
2337 2445 : do idir2=1,3
2338 2244 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2339 1395 : if(ipert1==ipert2.or.dyfr_nondiag==1)then
2340 1107 : if (dyfr_cplex==1) then
2341 1107 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2342 2214 : dyfrnl(1,idir1,idir2,ipert1,1+(ipert2-1)*dyfr_nondiag),zero
2343 : else
2344 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2345 0 : dyfrnl(:,idir1,idir2,ipert1,1+(ipert2-1)*dyfr_nondiag)
2346 : end if
2347 : else
2348 288 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,zero,zero
2349 : end if
2350 : end if
2351 : end do
2352 : end do
2353 : end if
2354 : end do
2355 : end do
2356 :
2357 : ! Now the nonlinear xc core correction(1) part
2358 28 : write(iout,*)' '
2359 28 : write(iout,*)' Frozen wf xc core (1) part',' of the dynamical matrix'
2360 28 : write(iout,*)' j1 j2 matrix element'
2361 28 : write(iout,*)' dir pert dir pert real part imaginary part'
2362 95 : do ipert1=1,natom
2363 296 : do idir1=1,3
2364 268 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2365 201 : write(iout,*)' '
2366 762 : do ipert2=1,natom
2367 2445 : do idir2=1,3
2368 2244 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2369 1395 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2370 1395 : dyfrx1(1,idir1,ipert1,idir2,ipert2),&
2371 2790 : dyfrx1(2,idir1,ipert1,idir2,ipert2)
2372 : end if
2373 : end do
2374 : end do
2375 : end if
2376 : end do
2377 : end do
2378 :
2379 : ! Now the nonlinear xc core correction(2) part
2380 28 : write(iout,*)' '
2381 28 : write(iout,*)' Frozen wf xc core (2) part',' of the dynamical matrix'
2382 28 : write(iout,*)' j1 j2 matrix element'
2383 28 : write(iout,*)' dir pert dir pert real part imaginary part'
2384 95 : do ipert1=1,natom
2385 296 : do idir1=1,3
2386 268 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2387 201 : write(iout,*)' '
2388 762 : do ipert2=1,natom
2389 2445 : do idir2=1,3
2390 2244 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2391 1395 : if(ipert1==ipert2)then
2392 459 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,dyfrx2(idir1,idir2,ipert2),zero
2393 : else
2394 936 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,zero,zero
2395 : end if
2396 : end if
2397 : end do
2398 : end do
2399 : end if
2400 : end do
2401 : end do
2402 :
2403 : ! Now the DFT-D vdw part of the dynamical matrix
2404 28 : if (usevdw==1) then
2405 0 : write(iout,*)' '
2406 0 : write(iout,*)' DFT-D van der Waals part of the dynamical matrix'
2407 0 : write(iout,*)' j1 j2 matrix element'
2408 0 : write(iout,*)' dir pert dir pert real part imaginary part'
2409 0 : do ipert1=1,natom
2410 0 : do idir1=1,3
2411 0 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2412 0 : write(iout,*)' '
2413 0 : do ipert2=1,natom
2414 0 : do idir2=1,3
2415 0 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2416 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2417 0 : dyvdw(1,idir1,ipert1,idir2,ipert2),&
2418 0 : dyvdw(2,idir1,ipert1,idir2,ipert2)
2419 : end if
2420 : end do
2421 : end do
2422 : end if
2423 : end do
2424 : end do
2425 : end if
2426 :
2427 : ! End of the phonon condition
2428 : end if
2429 :
2430 : ! In case of atom. strain/electric field perturbation (piezoelectric tensor)
2431 33 : if (pawpiezo==1.and.(rfpert(natom+2)==1.or.rfstrs/=0).and.outd2==1)then
2432 8 : write(iout,*)' '
2433 8 : write(iout,*)' Frozen wf part of the piezoelectric tensor'
2434 8 : write(iout,*)' j1 j2 matrix element'
2435 8 : write(iout,*)' dir pert dir pert real part imaginary part'
2436 8 : ipert1=natom+2
2437 32 : do idir1=1,3
2438 24 : write(iout,*)' '
2439 24 : ii=1
2440 80 : do ipert2=natom+3,natom+4
2441 216 : do idir2=1,3
2442 144 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,piezofrnl(ii,idir1),zero
2443 192 : ii=ii+1
2444 : end do
2445 : end do
2446 : end do
2447 : end if
2448 :
2449 : ! In case of atom. displ/electric field perturbation (Born Effective Charges)
2450 33 : if (pawbec==1.and.(rfpert(natom+2)==1.or.rfphon==1).and.outd2==1)then
2451 8 : write(iout,*)' '
2452 8 : write(iout,*)' Frozen wf part of the Born Effective Charges'
2453 8 : write(iout,*)' j1 j2 matrix element'
2454 8 : write(iout,*)' dir pert dir pert real part imaginary part'
2455 8 : ipert1 = natom+2
2456 34 : do ipert2=1,natom
2457 112 : do idir1=1,3
2458 78 : write(iout,*)' '
2459 338 : do idir2=1,3
2460 312 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,becfrnl(idir2,ipert2,idir1),zero
2461 : end do
2462 : end do
2463 : end do
2464 : end if
2465 :
2466 : ! In case of strain
2467 33 : if (rfstrs/=0)then
2468 :
2469 : ! Write the Ewald part of the elastic tensor
2470 15 : write(iout,*)' '
2471 15 : write(iout,*)' Ewald part of the elastic tensor in cartesian coordinates'
2472 15 : write(iout,*)' j1 j2 matrix element'
2473 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2474 45 : do ipert1=natom+3,natom+4
2475 135 : do idir1=1,3
2476 120 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2477 90 : ii=idir1+3*(ipert1-natom-3)
2478 90 : write(iout,*)' '
2479 270 : do ipert2=natom+3,natom+4
2480 810 : do idir2=1,3
2481 720 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2482 420 : jj=idir2+3*(ipert2-natom-3)
2483 420 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,elteew(ii,jj),zero
2484 : end if
2485 : end do
2486 : end do
2487 : end if
2488 : end do
2489 : end do
2490 :
2491 : ! Write the Ewald part of the internal strain coupling parameters
2492 15 : write(iout,*)' '
2493 15 : write(iout,*)' Ewald part of the internal strain coupling parameters'
2494 15 : write(iout,*)' (cartesian strain, reduced atomic coordinates)'
2495 15 : write(iout,*)' j1 j2 matrix element'
2496 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2497 58 : do ipert1=1,natom
2498 187 : do idir1=1,3
2499 172 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2500 129 : ii=idir1+6+3*(ipert1-1)
2501 129 : write(iout,*)' '
2502 387 : do ipert2=natom+3,natom+4
2503 1161 : do idir2=1,3
2504 1032 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2505 642 : jj=idir2+3*(ipert2-natom-3)
2506 642 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,elteew(ii,jj),zero
2507 : end if
2508 : end do
2509 : end do
2510 : end if
2511 : end do
2512 : end do
2513 :
2514 : ! Now the local frozen wf part
2515 15 : write(iout,*)' '
2516 15 : write(iout,*)' Frozen wf local part of the elastic tensor in cartesian coordinates'
2517 15 : write(iout,*)' j1 j2 matrix element'
2518 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2519 45 : do ipert1=natom+3,natom+4
2520 135 : do idir1=1,3
2521 120 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2522 90 : ii=idir1+3*(ipert1-natom-3)
2523 90 : write(iout,*)' '
2524 270 : do ipert2=natom+3,natom+4
2525 810 : do idir2=1,3
2526 720 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2527 420 : jj=idir2+3*(ipert2-natom-3)
2528 420 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2, eltfrloc(ii,jj),zero
2529 : end if
2530 : end do
2531 : end do
2532 : end if
2533 : end do
2534 : end do
2535 :
2536 15 : write(iout,*)' '
2537 15 : write(iout,*)' Frozen wf local part of the internal strain coupling parameters '
2538 15 : write(iout,*)' (cartesian strain, reduced atomic coordinates)'
2539 15 : write(iout,*)' j1 j2 matrix element'
2540 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2541 58 : do ipert1=1,natom
2542 187 : do idir1=1,3
2543 172 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2544 129 : ii=idir1+6+3*(ipert1-1)
2545 129 : write(iout,*)' '
2546 387 : do ipert2=natom+3,natom+4
2547 1161 : do idir2=1,3
2548 1032 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2549 642 : jj=idir2+3*(ipert2-natom-3)
2550 642 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2, eltfrloc(ii,jj),zero
2551 : end if
2552 : end do
2553 : end do
2554 : end if
2555 : end do
2556 : end do
2557 :
2558 : ! Now the nonlo frozen wf part
2559 15 : write(iout,*)' '
2560 15 : write(iout,*)' Frozen wf nonlocal part of the elastic tensor in cartesian coordinates'
2561 15 : write(iout,*)' j1 j2 matrix element'
2562 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2563 45 : do ipert1=natom+3,natom+4
2564 135 : do idir1=1,3
2565 120 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2566 90 : ii=idir1+3*(ipert1-natom-3)
2567 90 : write(iout,*)' '
2568 270 : do ipert2=natom+3,natom+4
2569 810 : do idir2=1,3
2570 720 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2571 420 : jj=idir2+3*(ipert2-natom-3)
2572 420 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2, eltfrnl(ii,jj),zero
2573 : end if
2574 : end do
2575 : end do
2576 : end if
2577 : end do
2578 : end do
2579 :
2580 15 : write(iout,*)' '
2581 15 : write(iout,*)' Frozen wf nonlocal part of the internal strain coupling parameters '
2582 15 : write(iout,*)' (cartesian strain, reduced atomic coordinates)'
2583 15 : write(iout,*)' j1 j2 matrix element'
2584 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2585 58 : do ipert1=1,natom
2586 187 : do idir1=1,3
2587 172 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2588 129 : ii=idir1+6+3*(ipert1-1)
2589 129 : write(iout,*)' '
2590 387 : do ipert2=natom+3,natom+4
2591 1161 : do idir2=1,3
2592 1032 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2593 642 : jj=idir2+3*(ipert2-natom-3)
2594 642 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,eltfrnl(ii,jj),zero
2595 : end if
2596 : end do
2597 : end do
2598 : end if
2599 : end do
2600 : end do
2601 :
2602 : ! Now the xc part
2603 15 : write(iout,*)' '
2604 15 : write(iout,*)' Frozen wf xc part of the elastic tensor in cartesian coordinates'
2605 15 : write(iout,*)' j1 j2 matrix element'
2606 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2607 45 : do ipert1=natom+3,natom+4
2608 135 : do idir1=1,3
2609 120 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2610 90 : ii=idir1+3*(ipert1-natom-3)
2611 90 : write(iout,*)' '
2612 270 : do ipert2=natom+3,natom+4
2613 810 : do idir2=1,3
2614 720 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2615 420 : jj=idir2+3*(ipert2-natom-3)
2616 420 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2, eltfrxc(ii,jj),zero
2617 : end if
2618 : end do
2619 : end do
2620 : end if
2621 : end do
2622 : end do
2623 :
2624 15 : write(iout,*)' '
2625 15 : write(iout,*)' Frozen wf xc part of the internal strain coupling parameters '
2626 15 : write(iout,*)' (cartesian strain, reduced atomic coordinates)'
2627 15 : write(iout,*)' j1 j2 matrix element'
2628 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2629 58 : do ipert1=1,natom
2630 187 : do idir1=1,3
2631 172 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2632 129 : ii=idir1+6+3*(ipert1-1)
2633 129 : write(iout,*)' '
2634 387 : do ipert2=natom+3,natom+4
2635 1161 : do idir2=1,3
2636 1032 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2637 642 : jj=idir2+3*(ipert2-natom-3)
2638 642 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,eltfrxc(ii,jj),zero
2639 : end if
2640 : end do
2641 : end do
2642 : end if
2643 : end do
2644 : end do
2645 :
2646 : ! Now the kinetic frozen wf part
2647 15 : write(iout,*)' '
2648 15 : write(iout,*)' Frozen wf kinetic part of the elastic tensor in cartesian coordinates'
2649 15 : write(iout,*)' j1 j2 matrix element'
2650 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2651 45 : do ipert1=natom+3,natom+4
2652 135 : do idir1=1,3
2653 120 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2654 90 : ii=idir1+3*(ipert1-natom-3)
2655 90 : write(iout,*)' '
2656 270 : do ipert2=natom+3,natom+4
2657 810 : do idir2=1,3
2658 720 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2659 420 : jj=idir2+3*(ipert2-natom-3)
2660 420 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,eltfrkin(ii,jj),zero
2661 : end if
2662 : end do
2663 : end do
2664 : end if
2665 : end do
2666 : end do
2667 :
2668 : ! Now the hartree frozen wf part
2669 15 : write(iout,*)' '
2670 15 : write(iout,*)' Frozen wf hartree part of the elastic tensor in cartesian coordinates'
2671 15 : write(iout,*)' j1 j2 matrix element'
2672 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2673 45 : do ipert1=natom+3,natom+4
2674 135 : do idir1=1,3
2675 120 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2676 90 : ii=idir1+3*(ipert1-natom-3)
2677 90 : write(iout,*)' '
2678 270 : do ipert2=natom+3,natom+4
2679 810 : do idir2=1,3
2680 720 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2681 420 : jj=idir2+3*(ipert2-natom-3)
2682 420 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2, eltfrhar(ii,jj),zero
2683 : end if
2684 : end do
2685 : end do
2686 : end if
2687 : end do
2688 : end do
2689 :
2690 : ! Now the psp core part
2691 15 : write(iout,*)' '
2692 15 : write(iout,*)' Psp core part of the elastic tensor in cartesian coordinates'
2693 15 : write(iout,*)' j1 j2 matrix element'
2694 15 : write(iout,*)' dir pert dir pert real part imaginary part'
2695 45 : do ipert1=natom+3,natom+4
2696 135 : do idir1=1,3
2697 120 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1) then
2698 90 : ii=idir1+3*(ipert1-natom-3)
2699 90 : write(iout,*)' '
2700 270 : do ipert2=natom+3,natom+4
2701 810 : do idir2=1,3
2702 720 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2703 420 : jj=idir2+3*(ipert2-natom-3)
2704 420 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,eltcore(ii,jj),zero
2705 : end if
2706 : end do
2707 : end do
2708 : end if
2709 : end do
2710 : end do
2711 :
2712 : ! Now the DFT-D vdw part
2713 15 : if (usevdw==1) then
2714 0 : write(iout,*)' '
2715 0 : write(iout,*)' DFT-D van der Waals part of the elastic tensor in cartesian coordinates'
2716 0 : write(iout,*)' j1 j2 matrix element'
2717 0 : write(iout,*)' dir pert dir pert real part imaginary part'
2718 0 : do ipert1=natom+3,natom+4
2719 0 : do idir1=1,3
2720 0 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1)then
2721 0 : ii=idir1+3*(ipert1-natom-3)
2722 0 : write(iout,*)' '
2723 0 : do ipert2=natom+3,natom+4
2724 0 : do idir2=1,3
2725 0 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2726 0 : jj=idir2+3*(ipert2-natom-3)
2727 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,eltvdw(ii,jj),zero
2728 : end if
2729 : end do
2730 : end do
2731 : end if
2732 : end do
2733 : end do
2734 :
2735 0 : write(iout,*)' '
2736 0 : write(iout,*)' DFT-D2 van der Waals part of the internal strain coupling parameters'
2737 0 : write(iout,*)' (cartesian strain, reduced atomic coordinates)'
2738 0 : write(iout,*)' j1 j2 matrix element'
2739 0 : write(iout,*)' dir pert dir pert real part imaginary part'
2740 0 : do ipert1=1,natom
2741 0 : do idir1=1,3
2742 0 : if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1 )then
2743 0 : ii=idir1+6+3*(ipert1-1)
2744 0 : write(iout,*)' '
2745 0 : do ipert2=natom+3,natom+4
2746 0 : do idir2=1,3
2747 0 : if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2748 0 : jj=idir2+3*(ipert2-natom-3)
2749 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,eltvdw(ii,jj),zero
2750 : end if
2751 : end do
2752 : end do
2753 : end if
2754 : end do
2755 : end do
2756 : end if ! usevdw
2757 :
2758 : end if ! strain condition
2759 :
2760 : ! Now the local nonstationary nonfrozenwf part
2761 33 : if (outd2==1)then
2762 33 : write(iout,*)' '
2763 33 : write(iout,*)' Non-stationary local part of the 2-order matrix'
2764 33 : write(iout,*)' j1 j2 matrix element'
2765 33 : write(iout,*)' dir pert dir pert real part imaginary part'
2766 342 : do ipert1=1,mpert
2767 1269 : do idir1=1,3
2768 : if ((ipert1<=natom .or.&
2769 927 : & (ipert1==natom+2.and.qzero==1.and.ddkfil(idir1)/=0)).or.&
2770 : & ((ipert1==natom+3.or.ipert1==natom+4).and.&
2771 309 : & (rfpert(natom+3)==1.or.rfpert(natom+4)==1)))then
2772 346 : write(iout,*)' '
2773 3748 : do ipert2=1,mpert
2774 14535 : do idir2=1,3
2775 13608 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2776 3668 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2777 3668 : & d2lo(1,idir1,ipert1,idir2,ipert2),&
2778 7336 : & d2lo(2,idir1,ipert1,idir2,ipert2)
2779 : end if
2780 : end do
2781 : end do
2782 : end if
2783 : end do
2784 : end do
2785 : end if
2786 :
2787 : ! Now the nonlocal nonstationary nonfrozenwf part
2788 33 : if (outd2==1)then
2789 33 : write(iout,*)' '
2790 33 : write(iout,*)' Non-stationary non-local part of the 2nd-order matrix'
2791 33 : write(iout,*)' j1 j2 matrix element'
2792 33 : write(iout,*)' dir pert dir pert real part imaginary part'
2793 342 : do ipert1=1,mpert
2794 1269 : do idir1=1,3
2795 : if ((ipert1<=natom .or.&
2796 927 : & (ipert1==natom+2.and.qzero==1.and.ddkfil(idir1)/=0)).or.&
2797 : & ((ipert1==natom+3.or.ipert1==natom+4).and.&
2798 309 : & (rfpert(natom+3)==1.or.rfpert(natom+4)==1)))then
2799 346 : write(iout,*)' '
2800 3748 : do ipert2=1,mpert
2801 14535 : do idir2=1,3
2802 13608 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2803 3668 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2804 3668 : & d2nl(1,idir1,ipert1,idir2,ipert2),&
2805 7336 : & d2nl(2,idir1,ipert1,idir2,ipert2)
2806 : end if
2807 : end do
2808 : end do
2809 : end if
2810 : end do
2811 : end do
2812 : end if
2813 :
2814 : ! Now the overlap change nonstationnary nonfrozenwf part (PAW only)
2815 33 : if (outd2==1.and.usepaw==1)then
2816 8 : write(iout,*)' '
2817 8 : write(iout,*)' PAW: Non-stationary WF-overlap part of the 2nd-order matrix'
2818 8 : write(iout,*)' j1 j2 matrix element'
2819 8 : write(iout,*)' dir pert dir pert real part imaginary part'
2820 90 : do ipert1=1,mpert
2821 336 : do idir1=1,3
2822 : if ((ipert1<=natom .or.&
2823 246 : & (ipert1==natom+2.and.qzero==1.and.ddkfil(idir1)/=0)).or.&
2824 : & ((ipert1==natom+3.or.ipert1==natom+4).and.&
2825 82 : & (rfpert(natom+3)==1.or.rfpert(natom+4)==1)))then
2826 138 : write(iout,*)' '
2827 1602 : do ipert2=1,mpert
2828 6102 : do idir2=1,3
2829 5856 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2830 2574 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2831 2574 : & d2ovl(1,idir1,ipert1,idir2,ipert2),&
2832 5148 : & d2ovl(2,idir1,ipert1,idir2,ipert2)
2833 : end if
2834 : end do
2835 : end do
2836 : end if
2837 : end do
2838 : end do
2839 : end if
2840 :
2841 : ! Now the 0-order local stationary nonfrozenwf part
2842 33 : if (outd2==2)then
2843 0 : write(iout,*)' '
2844 0 : write(iout,*)' Stationary 0-order local part of the 2nd-order matrix'
2845 0 : write(iout,*)' j1 j2 matrix element'
2846 0 : write(iout,*)' dir pert dir pert real part imaginary part'
2847 0 : do ipert1=1,mpert
2848 0 : do idir1=1,3
2849 0 : if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2850 0 : write(iout,*)' '
2851 0 : do ipert2=1,mpert
2852 0 : do idir2=1,3
2853 0 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2854 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2855 0 : & d2loc0(1,idir1,ipert1,idir2,ipert2),&
2856 0 : & d2loc0(2,idir1,ipert1,idir2,ipert2)
2857 : end if
2858 : end do
2859 : end do
2860 : end if
2861 : end do
2862 : end do
2863 : end if
2864 :
2865 : ! Now the stationary 0-order kinetic nonfrozenwf part
2866 33 : if (outd2==2)then
2867 0 : write(iout,*)' '
2868 0 : write(iout,*)' Stationary 0-order kinetic part of the 2nd-order matrix'
2869 0 : write(iout,*)' j1 j2 matrix element'
2870 0 : write(iout,*)' dir pert dir pert real part imaginary part'
2871 0 : do ipert1=1,mpert
2872 0 : do idir1=1,3
2873 0 : if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2874 0 : write(iout,*)' '
2875 0 : do ipert2=1,mpert
2876 0 : do idir2=1,3
2877 0 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2878 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2879 0 : & d2k0(1,idir1,ipert1,idir2,ipert2),&
2880 0 : & d2k0(2,idir1,ipert1,idir2,ipert2)
2881 : end if
2882 : end do
2883 : end do
2884 : end if
2885 : end do
2886 : end do
2887 : end if
2888 :
2889 : ! Now the stationary 0-order eigenvalue nonfrozenwf part
2890 33 : if (outd2==2)then
2891 0 : write(iout,*)' '
2892 0 : write(iout,*)' Stationary 0-order eigenvalue part of the' ,' 2nd-order matrix'
2893 0 : write(iout,*)' j1 j2 matrix element'
2894 0 : write(iout,*)' dir pert dir pert real part imaginary part'
2895 0 : do ipert1=1,mpert
2896 0 : do idir1=1,3
2897 0 : if(rfpert(ipert1)==1.and.rfdir(idir1)==1) then
2898 0 : write(iout,*)' '
2899 0 : do ipert2=1,mpert
2900 0 : do idir2=1,3
2901 0 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1) then
2902 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2903 0 : d2eig0(1,idir1,ipert1,idir2,ipert2),&
2904 0 : d2eig0(2,idir1,ipert1,idir2,ipert2)
2905 : end if
2906 : end do
2907 : end do
2908 : end if
2909 : end do
2910 : end do
2911 : end if
2912 :
2913 : ! Now the stationary potential-density nonfrozenwf part
2914 33 : if (outd2==2)then
2915 0 : write(iout,*)' '
2916 0 : write(iout,*)' Station. potential-density part of the ',&
2917 0 : & ' 2nd-order matrix'
2918 0 : write(iout,*)' (Note : include some xc core-correction) '
2919 0 : write(iout,*)' j1 j2 matrix element'
2920 0 : write(iout,*)' dir pert dir pert real part imaginary part'
2921 0 : do ipert1=1,mpert
2922 0 : do idir1=1,3
2923 0 : if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2924 0 : write(iout,*)' '
2925 0 : do ipert2=1,mpert
2926 0 : do idir2=1,3
2927 0 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2928 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2929 0 : & d2vn(1,idir1,ipert1,idir2,ipert2),&
2930 0 : & d2vn(2,idir1,ipert1,idir2,ipert2)
2931 : end if
2932 : end do
2933 : end do
2934 : end if
2935 : end do
2936 : end do
2937 : end if
2938 :
2939 : ! Now the stationary 0-order nonloc nonfrozenwf part
2940 33 : if (outd2==2)then
2941 0 : write(iout,*)' '
2942 0 : write(iout,*)' Stationary 0-order nonlocal part of the 2-order'&
2943 0 : & ,' matrix'
2944 0 : write(iout,*)' j1 j2 matrix element'
2945 0 : write(iout,*)' dir pert dir pert real part imaginary part'
2946 0 : do ipert1=1,mpert
2947 0 : do idir1=1,3
2948 0 : if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2949 0 : write(iout,*)' '
2950 0 : do ipert2=1,mpert
2951 0 : do idir2=1,3
2952 0 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2953 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2954 0 : & d2nl0(1,idir1,ipert1,idir2,ipert2),&
2955 0 : & d2nl0(2,idir1,ipert1,idir2,ipert2)
2956 : end if
2957 : end do
2958 : end do
2959 : end if
2960 : end do
2961 : end do
2962 : end if
2963 :
2964 : ! Now the stationary 1-order nonloc nonfrozenwf part
2965 33 : if (outd2==2)then
2966 0 : write(iout,*)' '
2967 0 : write(iout,*)' Stationary 1-order nonlocal part of the'&
2968 0 : & ,' 2nd-order matrix'
2969 0 : write(iout,*)' (or the ddk wf part of it, in case of',&
2970 0 : & ' an electric field perturbation )'
2971 0 : write(iout,*)' j1 j2 matrix element'
2972 0 : write(iout,*)' dir pert dir pert real part imaginary part'
2973 0 : do ipert1=1,mpert
2974 0 : do idir1=1,3
2975 0 : if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2976 0 : write(iout,*)' '
2977 0 : do ipert2=1,mpert
2978 0 : do idir2=1,3
2979 0 : if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2980 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2981 0 : & d2nl1(1,idir1,ipert1,idir2,ipert2),&
2982 0 : & d2nl1(2,idir1,ipert1,idir2,ipert2)
2983 : end if
2984 : end do
2985 : end do
2986 : end if
2987 : end do
2988 : end do
2989 : end if
2990 :
2991 : ! End of the long print out condition
2992 : end if
2993 :
2994 :
2995 : !Derivative database initialisation
2996 :
2997 : !Calculation of the number of elements
2998 418 : nelmts=0
2999 4271 : do ipert1=1,mpert
3000 15830 : do idir1=1,3
3001 127909 : do ipert2=1,mpert
3002 461547 : do idir2=1,3
3003 449988 : nelmts=nelmts+blkflg(idir1,ipert1,idir2,ipert2)
3004 : end do
3005 : end do
3006 : end do
3007 : end do
3008 :
3009 : !Now the whole 2nd-order matrix, but not in cartesian coordinates,
3010 : !and masses not included
3011 418 : write(iout,*)' '
3012 418 : write(iout,*)' 2nd-order matrix (non-cartesian coordinates,',' masses not included,'
3013 418 : write(iout,*)' asr not included )'
3014 418 : if(rfstrs/=0) then
3015 33 : write(iout,*)' cartesian coordinates for strain terms (1/ucvol factor '
3016 33 : write(iout,*)' for elastic tensor components not included) '
3017 : end if
3018 418 : write(iout,*)' j1 j2 matrix element'
3019 418 : write(iout,*)' dir pert dir pert real part imaginary part'
3020 418 : nline=1
3021 4271 : do ipert1=1,mpert
3022 15830 : do idir1=1,3
3023 11559 : if(nline/=0)write(iout,*)' '
3024 11559 : nline=0
3025 127909 : do ipert2=1,mpert
3026 461547 : do idir2=1,3
3027 449988 : if(blkflg(idir1,ipert1,idir2,ipert2)==1)then
3028 25555 : nline=nline+1
3029 25555 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3030 25555 : & d2matr(1,idir1,ipert1,idir2,ipert2),&
3031 51110 : & d2matr(2,idir1,ipert1,idir2,ipert2)
3032 : end if
3033 : end do
3034 : end do
3035 : end do
3036 : end do
3037 :
3038 : !Now the dynamical matrix
3039 418 : if(rfphon==1)then
3040 382 : write(iout,*)' '
3041 382 : write(iout,*)' Dynamical matrix, in cartesian coordinates,'
3042 382 : write(iout,*)' if specified in the inputs, asr has been imposed'
3043 382 : write(iout,*)' j1 j2 matrix element'
3044 382 : write(iout,*)' dir pert dir pert real part imaginary part'
3045 382 : nline=1
3046 1145 : do ipert1=1,natom
3047 3434 : do idir1=1,3
3048 2289 : if(nline/=0)write(iout,*)' '
3049 2289 : nline=0
3050 9421 : do ipert2=1,natom
3051 27765 : do idir2=1,3
3052 25476 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3053 14235 : nline=nline+1
3054 14235 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3055 14235 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3056 28470 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3057 : end if
3058 : end do
3059 : end do
3060 : end do
3061 : end do
3062 : end if
3063 :
3064 : !Now the dielectric tensor ! normal case
3065 418 : if(rfpert(natom+2)==1)then
3066 :
3067 92 : write(iout,*)' '
3068 92 : write(iout,*)' Dielectric tensor, in cartesian coordinates,'
3069 92 : write(iout,*)' j1 j2 matrix element'
3070 92 : write(iout,*)' dir pert dir pert real part imaginary part'
3071 92 : ipert1=natom+2
3072 92 : ipert2=natom+2
3073 92 : nline=1
3074 368 : do idir1=1,3
3075 276 : if(nline/=0)write(iout,*)' '
3076 276 : nline=0
3077 1196 : do idir2=1,3
3078 1104 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3079 776 : nline=nline+1
3080 776 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3081 776 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3082 1552 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3083 : end if
3084 : end do
3085 : end do
3086 :
3087 92 : if (prtbbb == 1) then
3088 :
3089 2 : delta(:,:) = zero
3090 2 : delta(1,1) = one ; delta(2,2) = one ; delta(3,3) = one
3091 :
3092 2 : write(iout,*)
3093 2 : write(iout,*)'Band by band decomposition of the dielectric tensor'
3094 2 : write(iout,*)' '
3095 :
3096 2 : write(iout,*)' Vacuum polarization'
3097 2 : write(iout,*)' j1 j2 matrix element'
3098 2 : write(iout,*)' dir pert dir pert real part imaginary part'
3099 2 : nline=1
3100 8 : do idir1=1,3
3101 6 : if(nline/=0)write(iout,*)' '
3102 6 : nline=0
3103 26 : do idir2=1,3
3104 18 : nline=nline+1
3105 18 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3106 42 : & delta(idir2,idir1),zero
3107 : end do
3108 : end do
3109 :
3110 18 : do iband = 1,mband
3111 16 : write(iout,*)' '
3112 16 : write(iout,*)' Dielectric tensor, in cartesian coordinates, for band',iband
3113 16 : write(iout,*)' j1 j2 matrix element'
3114 16 : write(iout,*)' dir pert dir pert real part imaginary part'
3115 16 : ipert1 = natom + 2
3116 16 : ipert2 = natom + 2
3117 16 : nline=1
3118 66 : do idir1=1,3
3119 48 : if(nline/=0)write(iout,*)' '
3120 48 : nline=0
3121 208 : do idir2=1,3
3122 192 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3123 : ! substract vacuum polarization
3124 144 : if (idir1 == idir2) then
3125 : d2cart_bbb(1,idir1,idir2,ipert2,iband,iband) = &
3126 48 : & d2cart_bbb(1,idir1,idir2,ipert2,iband,iband) - 1
3127 : end if
3128 144 : nline=nline+1
3129 144 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3130 144 : & d2cart_bbb(1,idir1,idir2,ipert2,iband,iband),&
3131 288 : & d2cart_bbb(2,idir1,idir2,ipert2,iband,iband)
3132 : end if
3133 : end do
3134 : end do
3135 : end do !iband
3136 :
3137 : end if !prtbbb
3138 :
3139 : end if ! end natom+2 dielectric output
3140 :
3141 : !Now the effective charges
3142 : !In case of the stationary calculation
3143 418 : if(outd2==2 .and. rfpert(natom+2)==1 .and.rfphon==1)then
3144 0 : write(iout,*)' '
3145 0 : write(iout,*)' Effective charges, in cartesian coordinates,'
3146 0 : write(iout,*)' if specified in the inputs, charge neutrality has been imposed'
3147 0 : write(iout,*)' j1 j2 matrix element'
3148 0 : write(iout,*)' dir pert dir pert real part imaginary part'
3149 0 : ipert1=natom+2
3150 0 : nline=1
3151 0 : do idir1=1,3
3152 0 : if(nline/=0)write(iout,*)' '
3153 0 : nline=0
3154 0 : do ipert2=1,natom
3155 0 : do idir2=1,3
3156 0 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3157 0 : nline=nline+1
3158 0 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3159 0 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3160 0 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3161 : end if
3162 : end do
3163 : end do
3164 : end do
3165 : end if
3166 :
3167 : !Now in case of the non-stationary calculation
3168 418 : if(outd2==1 .and. rfpert(natom+2)==1)then
3169 92 : write(iout,*)' '
3170 92 : if(usepaw==1.and..not.(has_allddk))then
3171 0 : write(iout,*)' Warning: Born effectives charges are not correctly computed'
3172 0 : write(iout,*)' you need all ddk perturbations!'
3173 : end if
3174 92 : write(iout,*)' Effective charges, in cartesian coordinates,'
3175 92 : write(iout,*)' (from electric field response) '
3176 92 : write(iout,*)' if specified in the inputs, charge neutrality has been imposed'
3177 92 : write(iout,*)' j1 j2 matrix element'
3178 92 : write(iout,*)' dir pert dir pert real part imaginary part'
3179 92 : ipert2=natom+2
3180 92 : nline=1
3181 368 : do idir2=1,3
3182 276 : if(nline/=0)write(iout,*)' '
3183 276 : nline=0
3184 1025 : do ipert1=1,natom
3185 2904 : do idir1=1,3
3186 2628 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3187 1765 : nline=nline+1
3188 1765 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3189 1765 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3190 3530 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3191 : end if
3192 : end do
3193 : end do
3194 : end do
3195 : end if
3196 :
3197 : if(outd2==1 .and. rfphon==1 .and. qzero==1&
3198 418 : & .and. ( (ddkfil(1)/=0.or.ddkfil(2)/=0.or.ddkfil(3)/=0) .or. &
3199 : & berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. berryopt==14 .or. berryopt==16 .or. berryopt==17 ) )then !!HONG need to test for fixed E and D
3200 90 : write(iout,*)' '
3201 90 : if(usepaw==1.and..not.(has_allddk))then
3202 0 : write(iout,*)' Warning: Born effectives charges are not correctly computed'
3203 0 : write(iout,*)' you need all ddk perturbations!'
3204 : end if
3205 90 : write(iout,*)' Effective charges, in cartesian coordinates,'
3206 90 : write(iout,*)' (from phonon response) '
3207 90 : write(iout,*)' if specified in the inputs, charge neutrality has been imposed'
3208 90 : write(iout,*)' j1 j2 matrix element'
3209 90 : write(iout,*)' dir pert dir pert real part imaginary part'
3210 90 : nline=1
3211 294 : do ipert2=1,natom
3212 906 : do idir2=1,3
3213 612 : if(nline/=0)write(iout,*)' '
3214 612 : nline=0
3215 612 : ipert1=natom+2
3216 2652 : do idir1=1,3
3217 2448 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3218 1738 : nline=nline+1
3219 1738 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3220 1738 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3221 3476 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3222 : end if
3223 : end do
3224 : end do
3225 : end do
3226 90 : write(iout,*)' '
3227 90 : write(iout,*)' '
3228 90 : write(iout,*)' '
3229 :
3230 90 : if (prtbbb == 1) then
3231 :
3232 2 : write(iout,*)'Band by band decomposition of the Born effective charges'
3233 2 : write(iout,*)' '
3234 2 : write(iout,*)'Ionic charges in cartesian coordinates'
3235 2 : write(iout,*)' j1 j2 matrix element'
3236 2 : write(iout,*)' dir pert dir pert real part imaginary part'
3237 2 : zr = zero
3238 2 : zi = zero
3239 6 : do ipert2=1,natom
3240 18 : do idir2=1,3
3241 12 : if(nline/=0)write(iout,*)' '
3242 12 : nline=0
3243 12 : ipert1=natom+2
3244 52 : do idir1=1,3
3245 36 : zr = zero
3246 36 : if (idir1 == idir2) then
3247 12 : zr = zion(typat(ipert2))
3248 : end if
3249 36 : nline=nline+1
3250 36 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3251 84 : & zr,zi
3252 : end do
3253 : end do
3254 : end do
3255 :
3256 18 : do iband = 1,mband
3257 16 : write(iout,*)' '
3258 16 : write(iout,*)' Effective charges, in cartesian coordinates, for band',iband
3259 16 : write(iout,*)' (from phonon response) '
3260 16 : write(iout,*)' if specified in the inputs, charge neutrality has been imposed'
3261 16 : write(iout,*)' j1 j2 matrix element'
3262 16 : write(iout,*)' dir pert dir pert real part imaginary part'
3263 16 : nline=1
3264 50 : do ipert2=1,natom
3265 144 : do idir2=1,3
3266 96 : if(nline/=0)write(iout,*)' '
3267 96 : nline=0
3268 96 : ipert1=natom+2
3269 416 : do idir1=1,3
3270 384 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3271 288 : nline=nline+1
3272 288 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3273 288 : & d2cart_bbb(1,idir1,idir2,ipert2,iband,iband),&
3274 576 : & d2cart_bbb(2,idir1,idir2,ipert2,iband,iband)
3275 : end if
3276 : end do
3277 : end do
3278 : end do
3279 : end do !iband
3280 : end if !prtbbb
3281 : end if ! end of print effective charges
3282 :
3283 : !Now the elastic tensor
3284 418 : if(rfstrs/=0) then
3285 33 : write(iout,*)' '
3286 33 : write(iout,*)' Rigid-atom elastic tensor , in cartesian coordinates,'
3287 33 : write(iout,*)' j1 j2 matrix element'
3288 33 : write(iout,*)' dir pert dir pert real part imaginary part'
3289 33 : nline=1
3290 99 : do ipert1=natom+3,natom+4
3291 297 : do idir1=1,3
3292 198 : if(nline/=0)write(iout,*)' '
3293 198 : nline=0
3294 660 : do ipert2=natom+3,natom+4
3295 1782 : do idir2=1,3
3296 1584 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3297 962 : nline=nline+1
3298 962 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3299 962 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3300 1924 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3301 : end if
3302 : end do
3303 : end do
3304 : end do
3305 : end do
3306 : end if
3307 :
3308 : !Now the internal strain coupling parameters
3309 418 : if(rfstrs/=0) then
3310 33 : write(iout,*)' '
3311 33 : write(iout,*)' Internal strain coupling parameters, in cartesian coordinates,'
3312 33 : write(iout,*)' zero average net force deriv. has been imposed '
3313 33 : write(iout,*)' j1 j2 matrix element'
3314 33 : write(iout,*)' dir pert dir pert real part imaginary part'
3315 33 : nline=1
3316 111 : do ipert1=1,natom
3317 345 : do idir1=1,3
3318 234 : if(nline/=0)write(iout,*)' '
3319 234 : nline=0
3320 780 : do ipert2=natom+3,natom+4
3321 2106 : do idir2=1,3
3322 1872 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3323 1101 : nline=nline+1
3324 1101 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3325 1101 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3326 2202 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3327 : end if
3328 : end do
3329 : end do
3330 : end do
3331 : end do
3332 : end if
3333 :
3334 : !Now the piezoelectric tensor
3335 418 : if(rfstrs/=0 .and. (ddkfil(1)/=0.or.ddkfil(2)/=0.or.ddkfil(3)/=0))then
3336 17 : write(iout,*)' '
3337 17 : if(usepaw==1.and..not.(has_allddk))then
3338 0 : write(iout,*)' Warning: Rigid-atom proper piezoelectric tensor is not correctly computed'
3339 0 : write(iout,*)' you need all ddk perturbations!'
3340 : end if
3341 17 : write(iout,*)' Rigid-atom proper piezoelectric tensor, in cartesian coordinates,'
3342 17 : write(iout,*)' (from strain response)'
3343 17 : write(iout,*)' j1 j2 matrix element'
3344 17 : write(iout,*)' dir pert dir pert real part imaginary part'
3345 17 : nline=1
3346 17 : ipert1=natom+2
3347 68 : do idir1=1,3
3348 51 : if(nline/=0)write(iout,*)' '
3349 51 : nline=0
3350 571 : do ipert2=natom+3,natom+4
3351 459 : do idir2=1,3
3352 408 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3353 265 : nline=nline+1
3354 265 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3355 265 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3356 530 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3357 : end if
3358 : end do
3359 : end do
3360 : end do
3361 : end if
3362 :
3363 : !Now the piezoelectric tensor
3364 : if(outd2==1 .and. (pawpiezo==1.and.rfpert(natom+2)==1)&
3365 418 : & .and. (ddkfil(1)/=0.or.ddkfil(2)/=0.or.ddkfil(3)/=0)) then
3366 11 : write(iout,*)' '
3367 11 : if(usepaw==1.and..not.(has_allddk))then
3368 0 : write(iout,*)' Warning: Rigid-atom proper piezoelectric tensor is not correctly computed'
3369 0 : write(iout,*)' you need all ddk perturbations!'
3370 : end if
3371 11 : if(usepaw==1.and..not.has_full_piezo)then
3372 10 : write(iout,*)' Warning: The rigid-atom proper piezoelectric tensor'
3373 10 : write(iout,*)' from electric field response requires nsym=1'
3374 : end if
3375 11 : if (has_full_piezo) then
3376 1 : write(iout,*)' Rigid-atom proper piezoelectric tensor, in cartesian coordinates,'
3377 1 : write(iout,*)' (from electric field response)'
3378 1 : write(iout,*)' j1 j2 matrix element'
3379 1 : write(iout,*)' dir pert dir pert real part imaginary part'
3380 1 : nline=1
3381 1 : ipert1=natom+2
3382 4 : do idir1=1,3
3383 3 : if(nline/=0)write(iout,*)' '
3384 3 : nline=0
3385 10 : do ipert2=natom+3,natom+4
3386 27 : do idir2=1,3
3387 24 : if(carflg(idir2,ipert2,idir1,ipert1)==1)then
3388 18 : nline=nline+1
3389 18 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir2,ipert2,idir1,ipert1,&
3390 18 : & d2cart(1,idir2,ipert2,idir1,ipert1),&
3391 36 : & d2cart(2,idir2,ipert2,idir1,ipert1)
3392 : end if
3393 : end do
3394 : end do
3395 : end do
3396 : end if
3397 : end if
3398 :
3399 : !Now the spin field quantities
3400 418 : if (rfmagn==1) then
3401 8 : write(iout,*)' '
3402 8 : write(iout,*)' Magnetic susceptibility, in cartesian coordinates'
3403 8 : write(iout,*)' (from uniform spin field response)'
3404 8 : write(iout,*)' j1 j2 matrix element'
3405 8 : write(iout,*)' dir pert dir pert real part imaginary part'
3406 8 : ipert1=natom+5
3407 8 : ipert2=natom+5
3408 8 : nline=1
3409 32 : do idir1=1,3
3410 24 : if (nline/=0) write(iout,*)' '
3411 24 : nline=0
3412 104 : do idir2=1,3
3413 96 : if (carflg(idir1,ipert1,idir2,ipert2)==1) then
3414 15 : nline=nline+1
3415 15 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3416 15 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3417 30 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3418 : end if
3419 : end do
3420 : end do
3421 : end if
3422 :
3423 418 : if (rfmagn==2) then
3424 3 : write(iout,*)' '
3425 3 : write(iout,*)' Local magnetic susceptibility, in cartesian coordinates'
3426 3 : write(iout,*)' (from local spin field response)'
3427 3 : write(iout,*)' j1 j2 matrix element'
3428 3 : write(iout,*)' dir pert dir pert real part imaginary part'
3429 3 : nline=1
3430 21 : do ipert1= natom+12,2*natom+11
3431 75 : do idir1=1,3
3432 54 : if(nline/=0)write(iout,*)' '
3433 54 : nline=0
3434 396 : do ipert2= natom+12,2*natom+11
3435 1350 : do idir2=1,3
3436 1296 : if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3437 105 : nline=nline+1
3438 105 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3439 105 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3440 210 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3441 : end if
3442 : end do
3443 : end do
3444 : end do
3445 : end do
3446 : end if
3447 :
3448 418 : if (rfmagn==1.and.rfpert(natom+2)==1) then
3449 2 : write(iout,*)' '
3450 2 : write(iout,*)' Magnetoelectric tensor, in cartesian coordinates'
3451 2 : write(iout,*)' (from magnetization induced by electric field)'
3452 2 : write(iout,*)' j1 j2 matrix element'
3453 2 : write(iout,*)' dir pert dir pert real part imaginary part'
3454 2 : ipert1=natom+5
3455 2 : ipert2=natom+2
3456 2 : nline=1
3457 8 : do idir1=1,3
3458 6 : if (nline/=0) write(iout,*)' '
3459 6 : nline=0
3460 26 : do idir2=1,3
3461 24 : if (carflg(idir1,ipert1,idir2,ipert2)==1) then
3462 10 : nline=nline+1
3463 10 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3464 10 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3465 20 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3466 : end if
3467 : end do
3468 : end do
3469 :
3470 2 : write(iout,*)' '
3471 2 : write(iout,*)' Magnetoelectric tensor, in cartesian coordinates'
3472 2 : write(iout,*)' (from polarization induced by Zeeman field)'
3473 2 : write(iout,*)' j1 j2 matrix element'
3474 2 : write(iout,*)' dir pert dir pert real part imaginary part'
3475 2 : ipert1=natom+2
3476 2 : ipert2=natom+5
3477 2 : nline=1
3478 8 : do idir1=1,3
3479 6 : if (nline/=0) write(iout,*)' '
3480 6 : nline=0
3481 26 : do idir2=1,3
3482 24 : if (carflg(idir1,ipert1,idir2,ipert2)==1) then
3483 10 : nline=nline+1
3484 10 : write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3485 10 : & d2cart(1,idir1,ipert1,idir2,ipert2),&
3486 20 : & d2cart(2,idir1,ipert1,idir2,ipert2)
3487 : end if
3488 : end do
3489 : end do
3490 :
3491 : end if
3492 :
3493 :
3494 418 : end subroutine dfpt_dyout
3495 : !!***
3496 :
3497 : !!****f* ABINIT/dfpt_gatherdy
3498 : !!
3499 : !! NAME
3500 : !! dfpt_gatherdy
3501 : !!
3502 : !! FUNCTION
3503 : !! Sum (gather) the different parts of the 2nd-order matrix,
3504 : !! to get the matrix of second-order derivatives (d2matr)
3505 : !! Then, generates the dynamical matrix, not including the masses,
3506 : !! but the correct non-cartesian coordinates ( => d2cart)
3507 : !!
3508 : !! INPUTS
3509 : !! asr= (0=> no acoustic sum rule [asr] imposed), (1 or 2=> asr is imposed) only for dynamical matrix at Gamma
3510 : !! becfrnl(3,natom,3*pawbec)=NL frozen contribution to Born Effective Charges (PAW only)
3511 : !! berryopt=option for berry phase treatment
3512 : !! blkflg(3,mpert,3,mpert)= ( 1 if the element of the dynamical
3513 : !! matrix has been calculated ; 0 otherwise )
3514 : !! chneut= (0=> no charge neutrality sum rule imposed), (1=> charge neutrality is imposed,
3515 : !! with equal repartition of the charge neutrality correction for the effective charges),
3516 : !! (2=> charge neutrality is imposed, with weighted repartition of the charge neutrality correction for the effective charges),
3517 : !! dyew(2,3,natom,3,natom)=Ewald part of the dyn.matrix
3518 : !! dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=frozen wf part of the dyn.matrix (except xc1)
3519 : !! dyfrx1(2,3,natom,3,natom)=xc core correction (1) part of the frozen-wf
3520 : !! part of the dynamical matrix.
3521 : !! dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
3522 : !! dyfr_nondiag=1 if dyfrwf is non diagonal with respect to atoms; 0 otherwise
3523 : !! dyvdw(2,3,natom,3,natom*usevdw)=vdw DFT-D part of the dynamical matrix
3524 : !! d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
3525 : !! second order derivatives
3526 : !! d2nfr(2,3,mpert,3,mpert)=non-frozen wf part of the 2nd-order matr
3527 : !! eltcore(6,6)=core contribution to the elastic tensor
3528 : !! elteew(6+3*natom,6)=Ewald contribution to the elastic tsenor
3529 : !! eltfrhar(6,6)=hartree contribution to the elastic tensor
3530 : !! eltfrkin(6,6)=kinetic contribution to the elastic tensor
3531 : !! eltfrloc(6+3*natom,6)=local psp contribution to the elastic tensor
3532 : !! eltfrnl(6+3*natom,6)=non-local psp contribution to the elastic tensor
3533 : !! eltfrxc(6+3*natom,6)=exchange-correlation contribution to the elastic tensor
3534 : !! eltvdw(6+3*natom,6*usevdw)=vdw DFT-D part of the elastic tensor
3535 : !! gprimd(3,3)=basis vector in the reciprocal space
3536 : !! mband=maximum number of bands
3537 : !! mpert =maximum number of ipert
3538 : !! natom=number of atoms in unit cell
3539 : !! ntypat=number of atom types
3540 : !! outd2=option for the output of the 2nd-order matrix :
3541 : !! if outd2=1, non-stationary part
3542 : !! if outd2=2, stationary part.
3543 : !! pawbec= flag for the computation of frozen part of Born Effective Charges (PAW only)
3544 : !! prtbbb=if 1, print the band-by-band decomposition, otherwise, prtbbb=0
3545 : !! rfpert(mpert)=define the perturbations
3546 : !! rprimd(3,3)=dimensional primitive translations (bohr)
3547 : !! typat(natom)=integer label of each type of atom (1,2,...)
3548 : !! ucvol=unit cell volume
3549 : !! usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use
3550 : !! zion(ntypat)=charge corresponding to the atom type
3551 : !!
3552 : !! OUTPUT
3553 : !! carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
3554 : !! 2DTE matrix has been calculated correctly ; 0 otherwise )
3555 : !! d2cart(2,3,mpert,3,mpert)=
3556 : !! dynamical matrix, effective charges, dielectric tensor,....
3557 : !! all in cartesian coordinates
3558 : !! d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)=
3559 : !! band by band decomposition of Born effective charges
3560 : !! (calculated from phonon-type perturbation) in cartesian coordinates
3561 : !! d2matr(2,3,mpert,3,mpert)=2nd-order matrix (masses non included,
3562 : !! no cartesian coordinates : simply second derivatives)
3563 : !!
3564 : !! SOURCE
3565 :
3566 418 : subroutine dfpt_gatherdy(asr,becfrnl,berryopt,blkflg,carflg,chneut,dyew,dyfrwf,dyfrx1,&
3567 418 : & dyfr_cplex,dyfr_nondiag,dyvdw,d2bbb,d2cart,d2cart_bbb,d2matr,d2nfr,&
3568 418 : & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
3569 418 : & gprimd,mband,mpert,natom,ntypat,outd2,pawbec,pawpiezo,piezofrnl,prtbbb,&
3570 418 : & rfpert,rprimd,typat,ucvol,usevdw,zion)
3571 :
3572 : !Arguments -------------------------------
3573 : !scalars
3574 : integer,intent(in) :: asr,berryopt,chneut,dyfr_cplex,dyfr_nondiag,mband,mpert,natom,ntypat,outd2
3575 : integer,intent(in) :: pawbec,pawpiezo,prtbbb,usevdw
3576 : real(dp),intent(in) :: ucvol
3577 : !arrays
3578 : integer,intent(in) :: rfpert(mpert),typat(natom)
3579 : integer,intent(inout) :: blkflg(3,mpert,3,mpert)
3580 : integer,intent(out) :: carflg(3,mpert,3,mpert)
3581 : real(dp),intent(in) :: becfrnl(3,natom,3*pawbec)
3582 : real(dp),intent(in) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb)
3583 : real(dp),intent(in) :: d2nfr(2,3,mpert,3,mpert),dyew(2,3,natom,3,natom)
3584 : real(dp),intent(in) :: dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)
3585 : real(dp),intent(in) :: dyfrx1(2,3,natom,3,natom),dyvdw(2,3,natom,3,natom*usevdw)
3586 : real(dp),intent(in) :: eltcore(6,6),elteew(6+3*natom,6)
3587 : real(dp),intent(in) :: eltfrhar(6,6),eltfrkin(6,6),eltfrloc(6+3*natom,6)
3588 : real(dp),intent(in) :: eltfrnl(6+3*natom,6),eltfrxc(6+3*natom,6)
3589 : real(dp),intent(in) :: eltvdw(6+3*natom,6*usevdw),gprimd(3,3)
3590 : real(dp),intent(in) :: piezofrnl(6,3*pawpiezo),rprimd(3,3),zion(ntypat)
3591 : real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert)
3592 : real(dp),intent(out) :: d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)
3593 : real(dp),intent(out) :: d2matr(2,3,mpert,3,mpert)
3594 :
3595 : !Local variables -------------------------
3596 : !scalars
3597 : integer :: iband,iblok,idir,idir1,idir2,ii,ipert,ipert1,ipert2
3598 : integer :: jj,nblok,selectz
3599 : character(len=500) :: msg
3600 : !arrays
3601 : integer :: flg1(3),flg2(3)
3602 : real(dp) :: vec1(3),vec2(3)
3603 : ! real(dp) :: ter(3,3) ! this variable appears commented out below
3604 418 : real(dp),allocatable :: d2tmp(:,:,:,:,:),d2work(:,:,:,:,:),elfrtot(:,:)
3605 :
3606 : ! *********************************************************************
3607 :
3608 :
3609 418 : if(outd2/=3)then
3610 :
3611 : ! Initialise the 2nd-derivative matrix
3612 1140800 : d2matr(:,:,:,:,:)=0.0_dp
3613 :
3614 : ! Copy first the non-frozen-part for all possible elements
3615 4271 : do ipert2=1,mpert
3616 15830 : do idir2=1,3
3617 127909 : do ipert1=1,mpert
3618 461547 : do idir1=1,3
3619 449988 : if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3620 76665 : do ii=1,2
3621 : d2matr(ii,idir1,ipert1,idir2,ipert2)=&
3622 76665 : & d2nfr(ii,idir1,ipert1,idir2,ipert2)
3623 : end do
3624 : end if
3625 : end do
3626 : end do
3627 : end do
3628 : end do
3629 :
3630 : ! For the dynamical matrix
3631 : ! Add the Ewald part and the xc1 part of the frozen-wf part
3632 : ! Add the vdw part (if any)
3633 1275 : do ipert2=1,natom
3634 3846 : do idir2=1,3
3635 10925 : do ipert1=1,natom
3636 32559 : do idir1=1,3
3637 29988 : if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3638 43662 : do ii=1,2
3639 : d2matr(ii,idir1,ipert1,idir2,ipert2)= &
3640 : & d2matr(ii,idir1,ipert1,idir2,ipert2)+ &
3641 : & dyew(ii,idir1,ipert1,idir2,ipert2) + &
3642 29108 : & dyfrx1(ii,idir1,ipert1,idir2,ipert2)
3643 43662 : if (usevdw==1) then
3644 : d2matr(ii,idir1,ipert1,idir2,ipert2)= &
3645 : & d2matr(ii,idir1,ipert1,idir2,ipert2)+ &
3646 108 : & dyvdw(ii,idir1,ipert1,idir2,ipert2)
3647 : end if
3648 : end do
3649 : end if
3650 : end do
3651 : end do
3652 : end do
3653 : end do
3654 :
3655 : ! For the dynamical matrix
3656 : ! Add the frozen-wavefunction part
3657 418 : if (dyfr_nondiag==0) then
3658 1176 : do ipert2=1,natom
3659 3534 : do idir2=1,3
3660 10218 : do idir1=1,3
3661 9432 : if( blkflg(idir1,ipert2,idir2,ipert2)==1 ) then
3662 : d2matr(1:dyfr_cplex,idir1,ipert2,idir2,ipert2)=&
3663 : & d2matr(1:dyfr_cplex,idir1,ipert2,idir2,ipert2)&
3664 10544 : & +dyfrwf(1:dyfr_cplex,idir1,idir2,ipert2,1)
3665 : end if
3666 : end do
3667 : end do
3668 : end do
3669 : else
3670 99 : do ipert2=1,natom
3671 348 : do ipert1=1,natom
3672 1067 : do idir2=1,3
3673 3237 : do idir1=1,3
3674 2988 : if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3675 : d2matr(1:dyfr_cplex,idir1,ipert1,idir2,ipert2)=&
3676 : & d2matr(1:dyfr_cplex,idir1,ipert1,idir2,ipert2)&
3677 4036 : & +dyfrwf(1:dyfr_cplex,idir1,idir2,ipert1,ipert2)
3678 : end if
3679 : end do
3680 : end do
3681 : end do
3682 : end do
3683 : end if
3684 :
3685 : ! Add the frozen-wavefunction part of Born Effective Charges
3686 418 : if (pawbec==1) then
3687 32 : ipert2=natom+2
3688 128 : do idir2=1,3 ! Direction of electric field
3689 365 : do ipert1=1,natom ! Atom
3690 1044 : do idir1=1,3 ! Direction of atom
3691 711 : if(blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3692 : d2matr(1,idir1,ipert1,idir2,ipert2)=&
3693 450 : & d2matr(1,idir1,ipert1,idir2,ipert2)+becfrnl(idir1,ipert1,idir2)
3694 : end if
3695 948 : if(blkflg(idir2,ipert2,idir1,ipert1)==1 ) then
3696 : d2matr(1,idir2,ipert2,idir1,ipert1)=&
3697 450 : & d2matr(1,idir2,ipert2,idir1,ipert1)+becfrnl(idir1,ipert1,idir2)
3698 : end if
3699 : end do
3700 : end do
3701 : end do
3702 : end if
3703 :
3704 : ! For piezoelectric tensor add a contribution in PAW case
3705 418 : if(pawpiezo==1) then
3706 15 : ipert2=natom+2
3707 60 : do idir2=1,3 ! Direction of electric field
3708 150 : do ipert1=natom+3,natom+4 ! Strain
3709 405 : do idir1=1,3
3710 270 : ii=idir1+3*(ipert1-natom-3)
3711 270 : if(blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3712 : d2matr(1,idir1,ipert1,idir2,ipert2)=&
3713 150 : d2matr(1,idir1,ipert1,idir2,ipert2)+piezofrnl(ii,idir2)
3714 : endif
3715 360 : if(blkflg(idir2,ipert2,idir1,ipert1)==1 ) then
3716 : d2matr(1,idir2,ipert2,idir1,ipert1)=&
3717 150 : d2matr(1,idir2,ipert2,idir1,ipert1)+piezofrnl(ii,idir2)
3718 : endif
3719 : end do
3720 : end do
3721 : end do
3722 : end if
3723 :
3724 : ! Internal strain and strain-strain
3725 :
3726 : ! Accumulate all frozen parts of the elastic tensor
3727 2090 : ABI_MALLOC(elfrtot,(6+3*natom,6))
3728 33400 : elfrtot(:,:)=elteew(:,:)+eltfrloc(:,:)+eltfrnl(:,:)+eltfrxc(:,:)
3729 17974 : elfrtot(1:6,1:6)=elfrtot(1:6,1:6)+eltcore(:,:)+eltfrhar(:,:)+eltfrkin(:,:)
3730 619 : if (usevdw==1) elfrtot(:,:)=elfrtot(:,:)+eltvdw(:,:)
3731 :
3732 1254 : do ipert2=natom+3,natom+4
3733 3762 : do idir2=1,3
3734 :
3735 : ! Internal strain components
3736 7650 : do ipert1=1,natom
3737 23076 : do idir1=1,3
3738 15426 : ii=idir1+6+3*(ipert1-1)
3739 15426 : jj=idir2+3*(ipert2-natom-3)
3740 15426 : if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3741 : d2matr(1,idir1,ipert1,idir2,ipert2)=&
3742 1101 : & d2matr(1,idir1,ipert1,idir2,ipert2)+elfrtot(ii,jj)
3743 : endif
3744 20568 : if( blkflg(idir2,ipert2,idir1,ipert1)==1 ) then
3745 : d2matr(1,idir2,ipert2,idir1,ipert1)=&
3746 1101 : & d2matr(1,idir2,ipert2,idir1,ipert1)+elfrtot(ii,jj)
3747 : endif
3748 : end do
3749 : end do
3750 :
3751 : ! Now, strain-strain 2nd derivatives
3752 8360 : do ipert1=natom+3,natom+4
3753 22572 : do idir1=1,3
3754 20064 : if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3755 962 : ii=idir1+3*(ipert1-natom-3)
3756 962 : jj=idir2+3*(ipert2-natom-3)
3757 : d2matr(1,idir1,ipert1,idir2,ipert2)=&
3758 962 : & d2matr(1,idir1,ipert1,idir2,ipert2)+elfrtot(ii,jj)
3759 : end if
3760 : end do
3761 : end do
3762 :
3763 : end do
3764 : end do
3765 418 : ABI_FREE(elfrtot)
3766 : ! End section for internal strain and strain-strain
3767 :
3768 : ! The second-order matrix has been computed.
3769 :
3770 : ! Filter now components smaller in absolute value than 1.0d-20,
3771 : ! for automatic testing reasons
3772 4271 : do ipert2=1,mpert
3773 15830 : do idir2=1,3
3774 127909 : do ipert1=1,mpert
3775 461547 : do idir1=1,3
3776 449988 : if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3777 76665 : do ii=1,2
3778 76665 : if( d2matr(ii,idir1,ipert1,idir2,ipert2)**2 < 1.0d-40)then
3779 17701 : d2matr(ii,idir1,ipert1,idir2,ipert2)=zero
3780 : end if
3781 : end do
3782 : end if
3783 : end do
3784 : end do
3785 : end do
3786 : end do
3787 :
3788 : ! Cartesian coordinates transformation
3789 418 : iblok=1 ; nblok=1
3790 :
3791 : ! In the case of finite electric field, the convention for the
3792 : ! direction of the electric field perturbation was NOT the usual convention ...
3793 : ! So, there is a transformation to the usual conventions
3794 : ! to be done first ...
3795 : if((berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. berryopt==14 .or. berryopt==16 .or. berryopt==17 ) &
3796 5434 : & .and. minval(abs(blkflg(:,natom+2,:,natom+2)))/=0)then !!HONG need to check for fixed D and E calculation
3797 13 : if(minval(abs(blkflg(:,natom+2,:,natom+2)-1))/=0)then
3798 : write(msg,'(5a)')&
3799 0 : & ' In case of finite electric field, and electric field perturbation,',ch10,&
3800 0 : & ' the three directions for the perturbations must be treated.',ch10,&
3801 0 : & ' Action : set idir to 1 1 1, or forget about finite electric field.'
3802 0 : ABI_ERROR(msg)
3803 : end if
3804 10 : do ipert=1,mpert
3805 37 : do idir=1,3
3806 90 : do ii=1,2
3807 216 : vec1(:)=d2matr(ii,idir,ipert,:,natom+2)
3808 216 : flg1(:)=blkflg(idir,ipert,:,natom+2)
3809 54 : call cart39(flg1,flg2,gprimd,1,1,rprimd,vec1,vec2)
3810 216 : d2matr(ii,idir,ipert,:,natom+2)=vec2(:)*two_pi
3811 243 : blkflg(idir,ipert,:,natom+2)=flg2(:)
3812 : end do
3813 : end do
3814 : end do
3815 10 : do ipert=1,mpert
3816 37 : do idir=1,3
3817 90 : do ii=1,2
3818 216 : vec1(:)=d2matr(ii,:,natom+2,idir,ipert)
3819 216 : flg1(:)=blkflg(:,natom+2,idir,ipert)
3820 54 : call cart39(flg1,flg2,gprimd,1,1,rprimd,vec1,vec2)
3821 216 : d2matr(ii,:,natom+2,idir,ipert)=vec2(:)*two_pi
3822 243 : blkflg(:,natom+2,idir,ipert)=flg2(:)
3823 : end do
3824 : end do
3825 : end do
3826 : ! Also to be done, a change of sign, that I do not understand (XG071110)
3827 : ! Perhaps due to d/dk replacing id/dk ? !
3828 31 : d2matr(:,:,natom+2,:,natom+2)=-d2matr(:,:,natom+2,:,natom+2)
3829 : end if
3830 :
3831 : call cart29(blkflg,d2matr,carflg,d2cart,&
3832 418 : & gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion)
3833 :
3834 : ! Band by band decomposition of the Born effective charges
3835 418 : if(prtbbb==1)then
3836 10 : ABI_MALLOC(d2work,(2,3,mpert,3,mpert))
3837 6 : ABI_MALLOC(d2tmp,(2,3,mpert,3,mpert))
3838 18 : do iband=1,mband
3839 39472 : d2work(:,:,:,:,:)=0.0_dp
3840 39472 : d2tmp(:,:,:,:,:)=0.0_dp
3841 4480 : d2work(:,:,natom+2,:,:) = d2bbb(:,:,:,:,iband,iband)
3842 : call cart29(blkflg,d2work,carflg,d2tmp,&
3843 16 : & gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion)
3844 :
3845 : ! Remove the ionic part
3846 48 : do ipert1=1,natom
3847 144 : do idir1=1,3
3848 : d2tmp(1,idir1,natom+2,idir1,ipert1) = &
3849 128 : & d2tmp(1,idir1,natom+2,idir1,ipert1) - zion(typat(ipert1))
3850 : end do
3851 : end do
3852 :
3853 4482 : d2cart_bbb(:,:,:,:,iband,iband)=d2tmp(:,:,natom+2,:,:)
3854 :
3855 : end do
3856 2 : ABI_FREE(d2tmp)
3857 2 : ABI_FREE(d2work)
3858 : end if ! prtbbb==1
3859 :
3860 : !
3861 : ! Now, the cartesian elements are ready for output
3862 : ! carflg give the information on their correctness
3863 : end if
3864 :
3865 : ! Imposition of the ASR on the analytical part of the DynMat
3866 : ! Assume that if asr/=0, the whole cartesian matrix is correct
3867 418 : if(asr/=0)then
3868 :
3869 785 : ABI_MALLOC(d2work,(2,3,mpert,3,mpert))
3870 157 : call asria_calc(asr,d2work,d2cart,mpert,natom)
3871 : ! The following line imposes ASR:
3872 157 : call asria_corr(asr,d2work,d2cart,mpert,natom)
3873 :
3874 157 : ABI_FREE(d2work)
3875 :
3876 : ! Imposition of the charge neutrality on the effective charges.
3877 157 : if(rfpert(natom+2)==1)then
3878 66 : selectz=0
3879 66 : call chneu9(chneut,d2cart,mpert,natom,ntypat,selectz,typat,zion)
3880 : end if
3881 :
3882 : end if
3883 :
3884 : !Additional operations on cartesian strain derivatives
3885 418 : if(rfpert(natom+3)==1 .or. rfpert(natom+4)==1) then
3886 : ! Impose zero-net-force condition on internal strain tensor
3887 99 : do ipert2=natom+3,natom+4
3888 297 : do idir2=1,3
3889 198 : vec1(:)=0.0_dp
3890 666 : do ipert1=1,natom
3891 2070 : do idir1=1,3
3892 1872 : if(carflg(idir1,ipert1,idir2,ipert2)==1) then
3893 1101 : vec1(idir1)=vec1(idir1)+d2cart(1,idir1,ipert1,idir2,ipert2)
3894 : end if
3895 : end do
3896 : end do
3897 792 : vec1(:)=vec1(:)/dble(natom)
3898 732 : do ipert1=1,natom
3899 2070 : do idir1=1,3
3900 1872 : if(carflg(idir1,ipert1,idir2,ipert2)==1) then
3901 : ! Note minus sign to convert gradients to forces
3902 : d2cart(1,idir1,ipert1,idir2,ipert2)=&
3903 1101 : & -(d2cart(1,idir1,ipert1,idir2,ipert2)-vec1(idir1))
3904 : end if
3905 : end do
3906 : end do
3907 : end do
3908 : end do
3909 : ! Divide strain 2nd derivative by ucvol to give elastic tensor
3910 99 : do ipert2=natom+3,natom+4
3911 682 : do idir2=1,3
3912 660 : do ipert1=natom+3,natom+4
3913 1782 : do idir1=1,3
3914 1584 : if(carflg(idir1,ipert1,idir2,ipert2)==1) then
3915 : d2cart(1,idir1,ipert1,idir2,ipert2)=&
3916 962 : & d2cart(1,idir1,ipert1,idir2,ipert2)/ucvol
3917 : end if
3918 : end do
3919 : end do
3920 : end do
3921 : end do
3922 : end if
3923 :
3924 : !calculate Born effective charges from electric field perturbation
3925 : !do ipert1=1,natom
3926 : !ter(:,:)=zero
3927 : !do idir1=1,3
3928 : !do ii=1,3
3929 : !do jj=1,3
3930 : !if(abs(gprimd(idir1,ii))>1.0d-10)then
3931 : !ter(idir1,ii)=ter(idir1,ii)+ d2nfr(1,idir1,natom+2,jj,ipert1)*gprimd(jj,ii)
3932 : !endif
3933 : !enddo
3934 : !enddo
3935 : !add zion to bec
3936 : !ter(idir1,idir1)=ter(idir1,idir1)+zion(typat(ipert1))
3937 : !enddo
3938 : !d2cart(1,:,ipert1,:,natom+2)=ter(:,:)
3939 : !enddo
3940 : !carflg(:,1:natom,:,natom+2)=1
3941 :
3942 : !Born effective charges from phonon perturbation
3943 : !do ipert1=1,natom
3944 : !ter(:,:)=zero
3945 : !do idir1=1,3
3946 : !do ii=1,3
3947 : !do jj=1,3
3948 : !if(abs(gprimd(idir1,ii))>1.0d-10)then
3949 : !ter(idir1,ii)=ter(idir1,ii)+ d2nfr(1,jj,ipert1,idir1,natom+2)*gprimd(jj,ii)
3950 : !endif
3951 : !enddo
3952 : !enddo
3953 : !! add zion to bec
3954 : !ter(idir1,idir1)=ter(idir1,idir1)+zion(typat(ipert1))
3955 : !enddo
3956 : !d2cart(1,:,natom+2,:,ipert1)=ter(:,:)
3957 : !enddo
3958 : !carflg(:,natom+2,:,1:natom)=1
3959 :
3960 :
3961 : !!Dielectric constant
3962 : !do ii=1,3
3963 : !do jj=1,3
3964 : !ter(ii,jj)=d2nfr(1,ii,natom+2,jj,natom+2)
3965 : !end do
3966 : !end do
3967 : !ter(:,:)=pi*four*ter(:,:)/ucvol
3968 : !
3969 : !do ii=1,3
3970 : !ter(ii,ii)=ter(ii,ii)+one
3971 : !end do
3972 : !d2cart(1,:,natom+2,:,natom+2)=ter(:,:)
3973 : !carflg(:,natom+2,:,1:natom+2)=1
3974 :
3975 : !DEBUG
3976 : !Debugging, but needed to make the stuff work on the IBM Dirac ? !
3977 : !write(std_out,*)' d2cart '
3978 : !ipert2=natom+2
3979 : !do idir2=1,3
3980 : !ipert1=natom+2
3981 : !do idir1=1,3
3982 : !write(std_out,'(5i4,2d20.10)' )idir1,ipert1,idir2,ipert2,&
3983 : !& carflg(idir1,ipert1,idir2,ipert2),&
3984 : !& d2cart(1,idir1,ipert1,idir2,ipert2),&
3985 : !& d2cart(2,idir1,ipert1,idir2,ipert2)
3986 : !end do
3987 : !end do
3988 : !ENDDEBUG
3989 :
3990 :
3991 418 : end subroutine dfpt_gatherdy
3992 : !!***
3993 :
3994 : !!****f* ABINIT/dfpt_dyfro
3995 : !! NAME
3996 : !! dfpt_dyfro
3997 : !!
3998 : !! FUNCTION
3999 : !! Compute the different parts of the frozen-wavefunction part of
4000 : !! the dynamical matrix, except the non-local one, computed previously.
4001 : !! Also (when installed) symmetrize the different part and their sum.
4002 : !!
4003 : !! INPUTS
4004 : !! atindx1(natom)=index table for atoms, inverse of atindx
4005 : !! dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
4006 : !! dyfr_nondiag=1 if dyfrnl and dyfrwf are non diagonal with respect to atoms; 0 otherwise
4007 : !! gmet(3,3)=reciprocal space metric (bohr^-2)
4008 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1)
4009 : !! gsqcut=cutoff on G^2 based on ecut
4010 : !! indsym(4,nsym,natom)=index showing transformation of atom labels
4011 : !! under symmetry operations (computed in symatm)
4012 : !! mgfft=maximum size of 1D FFTs
4013 : !! mpi_enreg=information about MPI parallelization
4014 : !! mqgrid=dimensioned number of q grid points for local psp spline
4015 : !! natom=number of atoms in unit cell
4016 : !! nattyp(ntypat)=number of atoms of each type
4017 : !! nfft=(effective) number of FFT grid points (for this processor)
4018 : !! ngfft(18)=contain all needed information about 3D FFT,
4019 : !! see ~abinit/doc/variables/vargs.htm#ngfft
4020 : !! nspden=number of spin-density components
4021 : !! nsym=number of symmetries in space group
4022 : !! ntypat=number of types of atoms
4023 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
4024 : !! n3xccc=dimension of the xccc3d array (0 or nfft).
4025 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
4026 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
4027 : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
4028 : !! qgrid(mqgrid)=q point array for local psp spline fits
4029 : !! qphon(3)=wavevector of the phonon
4030 : !! rhog(2,nfft)=electron density in G space
4031 : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
4032 : !! symq(4,2,nsym)=1 if symmetry preserves present qpoint. From littlegroup_q
4033 : !! symrec(3,3,nsym)=symmetries in reciprocal space
4034 : !! typat(natom)=integer type for each atom in cell
4035 : !! ucvol=unit cell volume (bohr**3).
4036 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
4037 : !! vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
4038 : !! vxc(nfft,nspden)=exchange-correlation potential (hartree) in real
4039 : !! space--only used when n1xccc/=0
4040 : !! xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
4041 : !! xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives,
4042 : !! for each type of atom, from psp
4043 : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
4044 : !! xred(3,natom)=reduced coordinates for atoms in unit cell
4045 : !!
4046 : !! OUTPUT
4047 : !! dyfrlo(3,3,natom)=frozen wavefunctions part of the dynamical matrix
4048 : !! (local only)
4049 : !! dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=
4050 : !! frozen wavefunctions part of the dynamical matrix
4051 : !! (local + non-local)
4052 : !! If NCPP, it depends on one atom
4053 : !! If PAW, it depends on two atoms
4054 : !! dyfrxc(3,3,natom)=frozen wavefunctions part of the dynamical matrix
4055 : !! (non-linear xc core correction)
4056 : !!
4057 : !! SIDE EFFECTS
4058 : !! Input/Output
4059 : !! dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=
4060 : !! frozen wavefunctions part of the dynamical matrix
4061 : !! (non-local only)
4062 : !! If NCPP, it depends on one atom
4063 : !! If PAW, it depends on two atoms
4064 : !!
4065 : !! SOURCE
4066 :
4067 450 : subroutine dfpt_dyfro(atindx1,dyfrnl,dyfrlo,dyfrwf,dyfrxc,dyfr_cplex,dyfr_nondiag,&
4068 450 : & dtset,gmet,gprimd,gsqcut,indsym,mgfft,mpi_enreg,mqgrid,natom,nattyp,&
4069 450 : & nfft,ngfft,nspden,nsym,ntypat,n1xccc,n3xccc,psps,pawtab,ph1d,qgrid,&
4070 450 : & qphon,rhog,rprimd,symq,symrec,typat,ucvol,usepaw,vlspl,vxc,&
4071 450 : & xcccrc,xccc1d,xccc3d,xred)
4072 :
4073 : !Arguments ------------------------------------
4074 : !scalars
4075 : integer,intent(in) :: dyfr_cplex,dyfr_nondiag,mgfft,mqgrid,n1xccc,n3xccc,natom,nfft,nspden
4076 : integer,intent(in) :: nsym,ntypat,usepaw
4077 : real(dp),intent(in) :: gsqcut,ucvol
4078 : type(dataset_type),intent(in) :: dtset
4079 : type(pseudopotential_type),intent(in) :: psps
4080 : type(MPI_type),intent(in) :: mpi_enreg
4081 : !arrays
4082 : integer,intent(in) :: atindx1(natom),indsym(4,nsym,natom),nattyp(ntypat)
4083 : integer,intent(in) :: ngfft(18),symq(4,2,nsym),symrec(3,3,nsym),typat(natom)
4084 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
4085 : real(dp),intent(in) :: qgrid(mqgrid),qphon(3),rhog(2,nfft)
4086 : real(dp),intent(in) :: vlspl(mqgrid,2,ntypat),vxc(nfft,nspden)
4087 : real(dp),intent(in) :: xccc1d(n1xccc,6,ntypat),xcccrc(ntypat),xred(3,natom)
4088 : real(dp),intent(inout) :: rprimd(3,3)
4089 : real(dp),intent(inout) :: dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag),xccc3d(n3xccc)
4090 : real(dp),intent(out) :: dyfrlo(3,3,natom),dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)
4091 : real(dp),intent(inout) :: dyfrxc(3,3,natom) !vz_i
4092 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
4093 :
4094 : !Local variables-------------------------------
4095 : !scalars
4096 : logical, parameter :: do_final_sym=.true.
4097 : integer :: iatom,jatom,n1,n2,n3,optatm,optdyfr,opteltfr,optgr,option
4098 : integer :: optn,optn2,optstr,optv
4099 : real(dp) :: eei
4100 : !arrays
4101 : integer :: qprtrb(3)
4102 : real(dp) :: dummy6(6),dum_strn(6),dum_strv(6)
4103 : real(dp) :: tsec(2),vprtrb(2)
4104 : real(dp) :: dum_atmrho(0),dum_atmvloc(0),dum_gauss(0),dum_grn(0),dum_grv(0),dum_eltfrxc(0)
4105 450 : real(dp),allocatable :: dyfrlo_tmp1(:,:,:),dyfrlo_tmp2(:,:,:,:,:),dyfrsym_tmp(:,:,:,:,:)
4106 450 : real(dp),allocatable :: gr_dum(:,:),v_dum(:),vxctotg(:,:)
4107 :
4108 : ! *************************************************************************
4109 :
4110 450 : if(nspden==4 .and. usepaw==1)then
4111 0 : ABI_WARNING('dfpt_dyfro : DFPT with nspden=4 works at the moment just for norm-conserving psp! (no paw support yet with nspden=4)')
4112 : end if
4113 :
4114 450 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
4115 :
4116 450 : if (usepaw==1 .or. psps%nc_xccc_gspace==1) then
4117 :
4118 : ! PAW or NC with nc_xccc_gspace: compute local psp and core charge contribs together
4119 : ! in reciprocal space
4120 : ! -----------------------------------------------------------------------
4121 45 : call timab(563,1,tsec)
4122 45 : if (n3xccc>0) then
4123 123 : ABI_MALLOC(v_dum,(nfft))
4124 123 : ABI_MALLOC(vxctotg,(2,nfft))
4125 853721 : v_dum(:)=vxc(:,1);if (nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2))
4126 41 : call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,1,ngfft,0)
4127 : call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),&
4128 41 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
4129 41 : ABI_FREE(v_dum)
4130 : else
4131 4 : ABI_MALLOC(vxctotg,(0,0))
4132 : end if
4133 45 : optatm=0;optdyfr=1;optgr=0;optstr=0;optv=1;optn=n3xccc/nfft;optn2=1;opteltfr=0
4134 : call atm2fft(atindx1,dum_atmrho,dum_atmvloc,dyfrxc,dyfrlo,dum_eltfrxc,&
4135 : & dum_gauss,gmet,gprimd,dum_grn,dum_grv,gsqcut,mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,&
4136 : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,qgrid,qprtrb,dtset%rcut,rhog,&
4137 45 : & rprimd,dum_strn,dum_strv,ucvol,usepaw,vxctotg,vxctotg,vxctotg,vprtrb,vlspl)
4138 45 : ABI_FREE(vxctotg)
4139 127 : if (n3xccc==0) dyfrxc=zero
4140 : else
4141 :
4142 : ! Norm-conserving: compute local psp contribution in reciprocal space
4143 : ! and core charge contribution in real space
4144 : ! -----------------------------------------------------------------------
4145 405 : option=4
4146 1215 : ABI_MALLOC(dyfrlo_tmp1,(3,3,natom))
4147 1215 : ABI_MALLOC(gr_dum,(3,natom))
4148 1215 : ABI_MALLOC(v_dum,(nfft))
4149 : call mklocl_recipspace(dyfrlo_tmp1,eei,gmet,gprimd,&
4150 : & gr_dum,gsqcut,dtset%icutcoul,dummy6,mgfft,mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,dtset%nkpt,&
4151 405 : & ntypat,option,ph1d,qgrid,qprtrb,dtset%rcut,rhog,rprimd,ucvol,dtset%vcutgeo,vlspl,vprtrb,v_dum)
4152 1209 : do iatom=1,natom
4153 : ! Reestablish correct order of atoms
4154 10857 : dyfrlo(1:3,1:3,atindx1(iatom))=dyfrlo_tmp1(1:3,1:3,iatom)
4155 : end do
4156 405 : ABI_FREE(dyfrlo_tmp1)
4157 405 : ABI_FREE(v_dum)
4158 405 : if(n1xccc/=0)then
4159 : call mkcore(dummy6,dyfrxc,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,&
4160 246 : & n1,n1xccc,n2,n3,option,rprimd,typat,ucvol,vxc,xcccrc,xccc1d,xccc3d,xred)
4161 : end if
4162 405 : ABI_FREE(gr_dum)
4163 : end if
4164 :
4165 : !Symmetrize dynamical matrix explicitly for given space group:
4166 :
4167 : !Symmetrize local part of the dynamical matrix dyfrlo:
4168 1350 : ABI_MALLOC(dyfrsym_tmp,(1,3,3,natom,1))
4169 900 : ABI_MALLOC(dyfrlo_tmp2,(1,3,3,natom,1))
4170 12644 : dyfrsym_tmp(1,:,:,:,1)=dyfrlo(:,:,:)
4171 450 : call dfpt_sydy(1,dyfrsym_tmp,indsym,natom,0,nsym,qphon,dyfrlo_tmp2,symq,symrec)
4172 12644 : dyfrlo(:,:,:)=dyfrlo_tmp2(1,:,:,:,1)
4173 : if (do_final_sym) then
4174 12644 : dyfrsym_tmp(1,:,:,:,1)=dyfrxc(:,:,:)
4175 450 : call dfpt_sydy(1,dyfrsym_tmp,indsym,natom,0,nsym,qphon,dyfrlo_tmp2,symq,symrec)
4176 12644 : dyfrxc(:,:,:)=dyfrlo_tmp2(1,:,:,:,1)
4177 : end if
4178 450 : ABI_FREE(dyfrsym_tmp)
4179 450 : ABI_FREE(dyfrlo_tmp2)
4180 :
4181 : !Symmetrize nonlocal part of the dynamical matrix dyfrnl:
4182 : !atindx1 is used to reestablish the correct order of atoms
4183 450 : if (dyfr_nondiag==0) then
4184 2030 : ABI_MALLOC(dyfrsym_tmp,(dyfr_cplex,3,3,natom,1))
4185 1212 : do iatom=1,natom
4186 18246 : dyfrsym_tmp(:,:,:,atindx1(iatom),1)=dyfrnl(:,:,:,iatom,1)
4187 : end do
4188 : else
4189 264 : ABI_MALLOC(dyfrsym_tmp,(dyfr_cplex,3,3,natom,natom))
4190 176 : do jatom=1,natom
4191 768 : do iatom=1,natom
4192 13606 : dyfrsym_tmp(:,:,:,atindx1(iatom),atindx1(jatom))=dyfrnl(:,:,:,iatom,jatom)
4193 : end do
4194 : end do
4195 : end if
4196 450 : call dfpt_sydy(dyfr_cplex,dyfrsym_tmp,indsym,natom,dyfr_nondiag,nsym,qphon,dyfrnl,symq,symrec)
4197 450 : ABI_FREE(dyfrsym_tmp)
4198 :
4199 : !Collect local, nl xc core, and non-local part
4200 : !of the frozen wf dynamical matrix.
4201 32302 : dyfrwf(:,:,:,:,:)=dyfrnl(:,:,:,:,:)
4202 450 : if (dyfr_nondiag==0) then
4203 10884 : dyfrwf(1,:,:,:,1)=dyfrwf(1,:,:,:,1)+dyfrlo(:,:,:)+dyfrxc(:,:,:)
4204 : else
4205 176 : do iatom=1,natom
4206 1760 : dyfrwf(1,:,:,iatom,iatom)=dyfrwf(1,:,:,iatom,iatom)+dyfrlo(:,:,iatom)+dyfrxc(:,:,iatom)
4207 : end do
4208 : end if
4209 :
4210 450 : end subroutine dfpt_dyfro
4211 : !!***
4212 :
4213 : !!****f* ABINIT/dfpt_dyxc1
4214 : !! NAME
4215 : !! dfpt_dyxc1
4216 : !!
4217 : !! FUNCTION
4218 : !! Compute 2nd-order non-linear xc core-correction (part1) to the dynamical matrix.
4219 : !! In case of derivative with respect to k or electric field perturbation,
4220 : !! the 1st-order local potential vanishes.
4221 : !!
4222 : !! INPUTS
4223 : !! atindx(natom)=index table for atoms
4224 : !! gmet(3,3)=metrix tensor in G space in Bohr**-2.
4225 : !! gsqcut=cutoff value on G**2 for sphere inside fft box.
4226 : !! ixc= choice of exchange-correlation scheme
4227 : !! kxc(nfft,nkxc)=first-order derivative of the xc potential
4228 : !! if (nkxc=1) LDA kxc(:,1)= d2Exc/drho2
4229 : !! if (nkxc=2) LDA kxc(:,1)= d2Exc/drho_up drho_up
4230 : !! kxc(:,2)= d2Exc/drho_up drho_dn
4231 : !! kxc(:,3)= d2Exc/drho_dn drho_dn
4232 : !! if (nkxc=7) GGA kxc(:,1)= d2Exc/drho2
4233 : !! kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
4234 : !! kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
4235 : !! kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
4236 : !! kxc(:,5)= gradx(rho)
4237 : !! kxc(:,6)= grady(rho)
4238 : !! kxc(:,7)= gradz(rho)
4239 : !! if (nkxc=19) spin-polarized GGA case (same as nkxc=7 with up and down components)
4240 : !! mgfft=maximum size of 1D FFTs
4241 : !! mpert=maximum number of ipert
4242 : !! mpi_enreg=information about MPI parallelization
4243 : !! mqgrid=number of grid pts in q array for f(q) spline.
4244 : !! natom=number of atoms in cell.
4245 : !! nfft=(effective) number of FFT grid points (for this processor)
4246 : !! ngfft(3)=fft grid dimensions.
4247 : !! nkxc=second dimension of the kxc array
4248 : !! (=1 for non-spin-polarized case, =3 for spin-polarized case)
4249 : !! nmxc= if true, handle density/potential as non-magnetic (even if it is)
4250 : !! nspden=number of spin-density components
4251 : !! ntypat=number of types of atoms in cell.
4252 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
4253 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
4254 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
4255 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
4256 : !! qgrid(mqgrid)=q grid for spline from 0 to qmax.
4257 : !! qphon(3)=wavevector of the phonon
4258 : !! rfdir(3)=array that define the directions of perturbations
4259 : !! rfpert(mpert)=array defining the type of perturbations that have to be computed
4260 : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
4261 : !! timrev=1 if time-reversal preserves the q wavevector; 0 otherwise.
4262 : !! typat(natom)=integer type for each atom in cell
4263 : !! ucvol=unit cell volume (bohr**3).
4264 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
4265 : !! xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
4266 : !! xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives,
4267 : !! for each type of atom, from psp
4268 : !! xred(3,natom)=fractional coordinates for atoms in unit cell
4269 : !!
4270 : !! OUTPUT! non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is)
4271 : !! blkflgfrx1(3,natom,3,natom)=flag to indicate whether an element has been computed or not
4272 : !! dyfrx1(2,3,natom,3,natom)=2nd-order non-linear xc
4273 : !! core-correction (part1) part of the dynamical matrix
4274 : !!
4275 : !! SOURCE
4276 :
4277 287 : subroutine dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,ixc,kxc,mgfft,mpert,mpi_enreg,mqgrid,&
4278 287 : & natom,nfft,ngfft,nkxc,nmxc,nspden,ntypat,n1xccc,psps,pawtab,&
4279 287 : & ph1d,qgrid,qphon,rfdir,rfpert,rprimd,timrev,typat,ucvol,usepaw,xcccrc,xccc1d,xred,rhor,vxc,ixcrot)
4280 :
4281 : use m_cgtools, only : dotprod_vn
4282 : use m_atm2fft, only : dfpt_atm2fft
4283 : use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll
4284 :
4285 : !Arguments ------------------------------------
4286 : !scalars
4287 : integer,intent(in) :: ixc,mgfft,mpert,mqgrid,n1xccc,natom,nfft,nkxc,nspden,ntypat
4288 : integer,intent(in) :: timrev,usepaw
4289 : integer,optional,intent(in) :: ixcrot
4290 : logical,intent(in) :: nmxc
4291 : real(dp),intent(in) :: gsqcut,ucvol
4292 : type(pseudopotential_type),intent(in) :: psps
4293 : type(MPI_type),intent(in) :: mpi_enreg
4294 : !arrays
4295 : integer,intent(in) :: atindx(natom),ngfft(18),rfdir(3),rfpert(mpert),typat(natom)
4296 : real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc)
4297 : real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid),qphon(3)
4298 : real(dp),intent(in) :: rprimd(3,3),xccc1d(n1xccc,6,ntypat),xcccrc(ntypat)
4299 : real(dp),intent(in) :: xred(3,natom)
4300 : integer,intent(out) :: blkflgfrx1(3,natom,3,natom)
4301 : real(dp),intent(out) :: dyfrx1(2,3,natom,3,natom)
4302 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
4303 : !optional
4304 : real(dp),optional,intent(in) :: rhor(nfft,nspden)
4305 : real(dp),optional,intent(in) :: vxc(nfft,nspden)
4306 :
4307 : !Local variables-------------------------------
4308 : !scalars
4309 : integer :: cplex,iat1,iatom1,iatom2,idir1,idir2,ierr,ifft,my_natom,comm_atom
4310 : integer :: n1,n2,n3,n3xccc,nfftot,option,upperdir,optnc
4311 : logical :: paral_atom
4312 : real(dp) :: valuei,valuer
4313 : !arrays
4314 287 : integer,pointer :: my_atmtab(:)
4315 : real(dp) :: tsec(2),gprimd_dummy(3,3)
4316 : real(dp) :: dum_nhat(0)
4317 287 : real(dp),allocatable :: rhor1(:,:),vxc10(:,:),xcccwk1(:),xcccwk2(:)
4318 : ! *********************************************************************
4319 :
4320 287 : call timab(182,1,tsec)
4321 :
4322 287 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
4323 287 : nfftot=n1*n2*n3
4324 :
4325 : !Set up parallelism over atoms
4326 287 : my_natom=mpi_enreg%my_natom
4327 287 : my_atmtab=>mpi_enreg%my_atmtab
4328 287 : comm_atom=mpi_enreg%comm_atom
4329 287 : paral_atom=(my_natom/=natom)
4330 :
4331 : !Zero out the output arrays :
4332 30735 : blkflgfrx1(:,:,:,:)=0
4333 72441 : dyfrx1(:,:,:,:,:)=zero
4334 :
4335 287 : cplex=2-timrev ; n3xccc=nfft
4336 1148 : ABI_MALLOC(vxc10,(cplex*nfft,nspden))
4337 :
4338 :
4339 : !Loop on the perturbation j1
4340 887 : do iat1=1,my_natom
4341 600 : iatom1=iat1; if(paral_atom)iatom1=my_atmtab(iat1)
4342 2687 : do idir1=1,3
4343 :
4344 : ! Compute the derivative of the core charge with respect to j1
4345 5400 : ABI_MALLOC(xcccwk1,(cplex*n3xccc))
4346 :
4347 : ! PAW or NC with nc_xccc_gspace: 1st-order core charge in reciprocal space
4348 1800 : if (usepaw==1 .or. psps%nc_xccc_gspace==1) then
4349 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd_dummy,gsqcut,idir1,iatom1,&
4350 : & mgfft,mqgrid,natom,1,nfft,ngfft,ntypat,ph1d,qgrid,&
4351 201 : & qphon,typat,ucvol,usepaw,xred,psps,pawtab,atmrhor1=xcccwk1,optn2_in=1)
4352 :
4353 : ! Norm-conserving psp: 1st-order core charge in real space
4354 : else
4355 : call dfpt_mkcore(cplex,idir1,iatom1,natom,ntypat,n1,n1xccc,&
4356 1599 : & n2,n3,qphon,rprimd,typat,ucvol,xcccrc,xccc1d,xcccwk1,xred)
4357 : end if
4358 :
4359 : ! Compute the corresponding potential
4360 1800 : option=0
4361 5400 : ABI_MALLOC(rhor1,(cplex*nfft,nspden))
4362 41821596 : rhor1=zero
4363 : !FR SPr EB Non-collinear magnetism
4364 1800 : if (nspden==4.and.present(rhor).and.present(vxc).and.present(ixcrot)) then
4365 90 : optnc=1
4366 : call dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,dum_nhat,0,dum_nhat,0,dum_nhat,0,nkxc,&
4367 90 : & nmxc,nspden,n3xccc,optnc,option,qphon,rhor,rhor1,rprimd,0,vxc,vxc10,xcccwk1,ixcrot=ixcrot)
4368 : else
4369 : call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,dum_nhat,0,dum_nhat,0,nkxc,&
4370 1710 : & nmxc,nspden,n3xccc,option,qphon,rhor1,rprimd,0,vxc10,xcccwk1)
4371 : end if
4372 1800 : ABI_FREE(rhor1)
4373 1800 : ABI_FREE(xcccwk1)
4374 :
4375 : ! vxc10 will couple with xcccwk2, that behaves like
4376 : ! a total density (ispden=1). Only the spin-up + spin-down
4377 : ! average of vxc10 is needed.
4378 1800 : if (nspden/=1)then
4379 1336926 : do ifft=1,cplex*nfft
4380 1336926 : vxc10(ifft,1)=(vxc10(ifft,1)+vxc10(ifft,2))*half
4381 : end do
4382 : end if
4383 :
4384 : ! Loop on the perturbation j2
4385 6261 : do iatom2=1,iatom1
4386 3861 : upperdir=3
4387 3861 : if(iatom1==iatom2)upperdir=idir1
4388 15444 : do idir2=1,upperdir
4389 9783 : if( (rfpert(iatom1)==1 .and. rfdir(idir1) == 1) .or. &
4390 3861 : & (rfpert(iatom2)==1 .and. rfdir(idir2) == 1) )then
4391 :
4392 : ! Compute the derivative of the core charge with respect to j2
4393 11472 : ABI_MALLOC(xcccwk2,(cplex*n3xccc))
4394 :
4395 : ! PAW or NC with nc_xccc_gspace: 1st-order core charge in reciprocal space
4396 5736 : if (usepaw==1 .or. psps%nc_xccc_gspace==1) then
4397 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd_dummy,gsqcut,idir2,iatom2,&
4398 : & mgfft,mqgrid,natom,1,nfft,ngfft,ntypat,ph1d,qgrid,&
4399 933 : & qphon,typat,ucvol,usepaw,xred,psps,pawtab,atmrhor1=xcccwk2,optn2_in=1)
4400 :
4401 : ! Norm-conserving psp: 1st-order core charge in real space
4402 : else
4403 : call dfpt_mkcore(cplex,idir2,iatom2,natom,ntypat,n1,n1xccc,&
4404 4803 : & n2,n3,qphon,rprimd,typat,ucvol,xcccrc,xccc1d,xcccwk2,xred)
4405 : end if
4406 :
4407 : ! Get the matrix element j1,j2
4408 :
4409 5736 : call dotprod_vn(cplex,xcccwk2,valuer,valuei,nfft,nfftot,1,2,vxc10,ucvol)
4410 :
4411 5736 : ABI_FREE(xcccwk2)
4412 :
4413 5736 : dyfrx1(1,idir1,iatom1,idir2,iatom2)= valuer
4414 5736 : dyfrx1(2,idir1,iatom1,idir2,iatom2)= valuei
4415 5736 : dyfrx1(1,idir2,iatom2,idir1,iatom1)= valuer
4416 5736 : dyfrx1(2,idir2,iatom2,idir1,iatom1)=-valuei
4417 5736 : blkflgfrx1(idir1,iatom1,idir2,iatom2)=1
4418 5736 : blkflgfrx1(idir2,iatom2,idir1,iatom1)=1
4419 : end if
4420 : end do
4421 : end do
4422 : end do
4423 : end do
4424 :
4425 287 : if (paral_atom) then
4426 22 : call timab(48,1,tsec)
4427 22 : call xmpi_sum(dyfrx1,comm_atom,ierr)
4428 22 : call xmpi_sum(blkflgfrx1,comm_atom,ierr)
4429 22 : call timab(48,2,tsec)
4430 : end if
4431 :
4432 287 : ABI_FREE(vxc10)
4433 :
4434 287 : call timab(182,2,tsec)
4435 :
4436 574 : end subroutine dfpt_dyxc1
4437 : !!***
4438 :
4439 : end module m_respfn_driver
4440 : !!***
|
