Line data Source code
1 : !!****m* ABINIT/m_nonlinear
2 : !! NAME
3 : !! m_nonlinear
4 : !!
5 : !! FUNCTION
6 : !! DFT calculations of non linear response functions.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2002-2026 ABINIT group (MVeithen,MB,LB)
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 : module m_nonlinear
23 :
24 : use defs_basis
25 : use defs_wvltypes
26 : use m_wffile
27 : use m_errors
28 : use m_abicore
29 : use m_xmpi
30 : use m_hdr
31 : use m_ebands
32 : use m_xcdata
33 : use m_dtset
34 : use m_dtfil
35 :
36 : use defs_datatypes, only : pseudopotential_type
37 : use defs_abitypes, only : MPI_type
38 : use m_fstrings, only : sjoin, itoa
39 : use m_time, only : timab
40 : use m_symtk, only : symmetrize_xred, littlegroup_q
41 : use m_dynmat, only : d3sym, sytens
42 : use m_ddb, only : ddb_type, nlopt
43 : use m_ddb_hdr, only : ddb_hdr_type
44 : use m_ioarr, only : read_rhor
45 : use m_kg, only : getcut, kpgio, getph
46 : use m_fft, only : fourdp
47 : use m_kpts, only : getkgrid
48 : use m_inwffil, only : inwffil
49 : use m_spacepar, only : hartre, setsym
50 : use m_pawfgr, only : pawfgr_type,pawfgr_init, pawfgr_destroy
51 : use m_pawang, only : pawang_type, pawang_init, pawang_free
52 : use m_pawrad, only : pawrad_type
53 : use m_pawtab, only : pawtab_type,pawtab_get_lsize
54 : use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_print
55 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_print
56 : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
57 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, &
58 : pawrhoij_bcast, pawrhoij_nullify, pawrhoij_inquire_dim
59 : use m_paw_energies,only : paw_energies_type
60 : use m_pawdij, only : pawdij, symdij
61 : use m_paw_finegrid,only : pawexpiqr
62 : use m_pawxc, only : pawxc_get_nkxc
63 : use m_paw_dmft, only : paw_dmft_type
64 : use m_paw_sphharm, only : setsym_ylm
65 : use m_paw_nhat, only : nhatgrid,pawmknhat
66 : use m_paw_denpot, only : pawdenpot
67 : use m_paw_init, only : pawinit,paw_gencond
68 : use m_paw_tools, only : chkpawovlp
69 : use m_mkrho, only : mkrho
70 : use m_getshell, only : getshell
71 : use m_pspini, only : pspini
72 : use m_atm2fft, only : atm2fft
73 : use m_rhotoxc, only : rhotoxc
74 : use m_drivexc, only : check_kxc
75 : use m_mpinfo, only : proc_distrb_cycle
76 : use m_mklocl, only : mklocl
77 : use m_common, only : setup1
78 : use m_fourier_interpol, only : transgrid
79 : use m_paw_occupancies, only : initrhoij
80 : use m_paw_correlations, only : pawpuxinit
81 : use m_mkcore, only : mkcore
82 : use m_pead_nl_loop, only : pead_nl_loop
83 : use m_dfptnl_loop, only : dfptnl_loop
84 :
85 : implicit none
86 :
87 : private
88 : !!***
89 :
90 : public :: nonlinear
91 : public :: dfptnl_doutput ! Write the matrix of third-order derivatives to the output file
92 : !!***
93 :
94 : contains
95 : !!***
96 :
97 : !!****f* ABINIT/nonlinear
98 : !! NAME
99 : !! nonlinear
100 : !!
101 : !! FUNCTION
102 : !! Primary routine for conducting DFT calculations of non linear response functions.
103 : !!
104 : !! INPUTS
105 : !! codvsn = code version
106 : !! dtfil <type(datafiles_type)> = variables related to files
107 : !! dtset <type(dataset_type)> = all input variables for this dataset
108 : !! etotal = new total energy (no meaning at output)
109 : !! mpi_enreg=information about MPI pnarallelization
110 : !! occ(mband*nkpt*nsppol) = occupation number for each band and k
111 : !! xred(3,natom) = reduced atomic coordinates
112 : !!
113 : !! OUTPUT
114 : !!
115 : !! npwtot(nkpt) = total number of plane waves at each k point
116 : !!
117 : !! SIDE EFFECTS
118 : !! pawang <type(pawang_type)>=paw angular mesh and related data
119 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
120 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
121 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
122 : !!
123 : !! NOTES
124 : !! USE OF FFT GRIDS:
125 : !! =================
126 : !! In case of PAW:
127 : !! ---------------
128 : !! Two FFT grids are used:
129 : !! - A "coarse" FFT grid (defined by ecut)
130 : !! for the application of the Hamiltonian on the plane waves basis.
131 : !! It is defined by nfft, ngfft, mgfft, ...
132 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ...
133 : !! are expressed on this grid.
134 : !! - A "fine" FFT grid (defined) by ecutdg)
135 : !! for the computation of the density inside PAW spheres.
136 : !! It is defined by nfftf, ngfftf, mgfftf, ...
137 : !! Total density, potentials, ...
138 : !! are expressed on this grid.
139 : !! In case of norm-conserving:
140 : !! ---------------------------
141 : !! - Only the usual FFT grid (defined by ecut) is used.
142 : !! It is defined by nfft, ngfft, mgfft, ...
143 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf)
144 : !! are set equal to (nfft,ngfft,mgfft) in that case.
145 : !!
146 : !! SOURCE
147 :
148 28 : subroutine nonlinear(codvsn,dtfil,dtset,etotal,mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,psps,xred)
149 :
150 : !Arguments ------------------------------------
151 : !scalars
152 : real(dp),intent(inout) :: etotal
153 : character(len=8),intent(in) :: codvsn
154 : type(MPI_type),intent(inout) :: mpi_enreg
155 : type(datafiles_type),intent(in) :: dtfil
156 : type(dataset_type),intent(inout) :: dtset
157 : type(pawang_type),intent(inout) :: pawang
158 : type(pseudopotential_type),intent(inout) :: psps
159 : !arrays
160 : integer,intent(out) :: npwtot(dtset%nkpt)
161 : real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),xred(3,dtset%natom)
162 : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
163 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
164 :
165 : !Local variables-------------------------------
166 : !scalars
167 : logical :: paral_atom,call_pawinit,qeq0
168 : integer,parameter :: level=50,formeig=0,response=1,cplex1=1
169 : integer :: ask_accurate,band_index,bantot,cplex,cplex_rhoij,dum_nshiftk,flag,gnt_option,gscase
170 : integer :: has_dijnd,has_diju,has_kxc,has_k3xc
171 : integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert
172 : integer :: iatom,indx,iband,ider,idir,ierr,ifft,ikpt,ipert,isppol
173 : integer :: ireadwf0,iscf_eff,ispden,itypat,izero,mcg,me,mgfftf,mkmem_max,mpert,my_natom
174 : integer :: n1,n3xccc,natom,nband_k,nfftf,nfftot,nfftotf,nhatdim,nhatgrdim
175 : integer :: nkpt_eff,nkpt_max,nkpt3,nkxc,nkxc1,nk3xc,nk3xc1,nneigh,ntypat,nsym1,nspden_rhoij,nzlmopt
176 : integer :: optcut,optgr0,optgr1,optgr2,optrad,option,optorth
177 : integer :: optatm,optdyfr,opteltfr,optgr,optstr,optv,optn,optn2
178 : integer :: psp_gencond,pead,qphase_rhoij,rdwr,rdwrpaw,spaceworld,tim_mkrho,timrev
179 : integer :: use_sym,usecprj,usexcnhat
180 : logical :: is_dfpt=.true.,nmxc
181 : real(dp),parameter :: k0(3)=(/zero,zero,zero/)
182 : real(dp) :: boxcut,compch_fft,compch_sph,ecore,ecut_eff,ecutdg_eff,ecutf
183 : real(dp) :: eei,bigexc,bigsxc,etot,fermie,fermih
184 : real(dp) :: gsqcut,gsqcut_eff,gsqcutc_eff
185 : real(dp) :: rdum,residm,ucvol,vxcavg,el_temp
186 : character(len=500) :: msg
187 : character(len=30) :: small_msg
188 : character(len=fnlen) :: dscrpt
189 28 : type(pawang_type) :: pawang1
190 : type(paw_energies_type) :: paw_energies
191 28 : type(ebands_t) :: bstruct
192 28 : type(hdr_type) :: hdr,hdr_den
193 1456 : type(ddb_hdr_type) :: ddb_hdr
194 28 : type(ddb_type) :: ddb
195 : type(wffile_type) :: wffgs,wfftgs
196 28 : type(wvl_data) :: wvl
197 : type(xcdata_type) :: xcdata
198 : !arrays
199 : integer :: dum_kptrlatt(3,3),dum_vacuum(3),ngfft(18),ngfftf(18),perm(6),ii,theunit
200 28 : integer,allocatable :: atindx(:),atindx1(:),blkflg(:,:,:,:,:,:),carflg(:,:,:,:,:,:),cgindex(:,:)
201 28 : integer,allocatable :: flg_tmp(:,:,:,:,:,:)
202 28 : integer,allocatable :: d3e_pert1(:),d3e_pert2(:),d3e_pert3(:)
203 28 : integer,allocatable :: indsym(:,:,:),indsy1(:,:,:),irrzon(:,:,:),irrzon1(:,:,:)
204 28 : integer,allocatable :: kg(:,:),kneigh(:,:),kg_neigh(:,:,:)
205 28 : integer,allocatable :: kptindex(:,:),l_size_atm(:)
206 28 : integer,allocatable :: npwarr(:),nattyp(:),pwind(:,:,:),rfpert(:,:,:,:,:,:)
207 28 : integer,allocatable :: symq(:,:,:),symrec(:,:,:),symaf1(:),symrc1(:,:,:),symrl1(:,:,:)
208 : real(dp) :: dum_gauss(0),dum_dyfrn(0),dum_dyfrv(0),dum_eltfrxc(0)
209 : real(dp) :: dum_grn(0),dum_grv(0),dum_rhog(0),dum_vg(0)
210 : real(dp) :: dum_shiftk(3,MAX_NSHIFTK),dummy6(6),other_dummy6(6),gmet(3,3),gprimd(3,3)
211 : real(dp) :: qphon(3),rmet(3,3),rprimd(3,3),tsec(2)
212 28 : real(dp),allocatable :: cg(:,:),d3cart(:,:,:,:,:,:,:)
213 56 : real(dp),allocatable :: d3etot(:,:,:,:,:,:,:),dum_kptns(:,:)
214 : ! We need all these arrays instead of one because in Fortran the maximum number of dimensions is 7...
215 28 : real(dp),allocatable :: d3e_1(:,:,:,:,:,:,:),d3cart_1(:,:,:,:,:,:,:)
216 28 : real(dp),allocatable :: d3e_2(:,:,:,:,:,:,:),d3cart_2(:,:,:,:,:,:,:)
217 28 : real(dp),allocatable :: d3e_3(:,:,:,:,:,:,:),d3cart_3(:,:,:,:,:,:,:)
218 28 : real(dp),allocatable :: d3e_4(:,:,:,:,:,:,:),d3cart_4(:,:,:,:,:,:,:)
219 28 : real(dp),allocatable :: d3e_5(:,:,:,:,:,:,:),d3cart_5(:,:,:,:,:,:,:)
220 28 : real(dp),allocatable :: d3e_6(:,:,:,:,:,:,:),d3cart_6(:,:,:,:,:,:,:)
221 28 : real(dp),allocatable :: d3e_7(:,:,:,:,:,:,:),d3cart_7(:,:,:,:,:,:,:)
222 28 : real(dp),allocatable :: d3e_8(:,:,:,:,:,:,:),d3cart_8(:,:,:,:,:,:,:)
223 28 : real(dp),allocatable :: d3e_9(:,:,:,:,:,:,:),d3cart_9(:,:,:,:,:,:,:)
224 28 : real(dp),allocatable :: dum_wtk(:),dyfrlo_indx(:,:,:),dyfrx2(:,:,:),eigen0(:)
225 28 : real(dp),allocatable :: grtn_indx(:,:),grxc(:,:),k3xc(:,:),kpt3(:,:),kxc(:,:)
226 28 : real(dp),allocatable :: mvwtk(:,:),nhat(:,:),nhatgr(:,:,:),ph1d(:,:),ph1df(:,:),phnons(:,:,:),phnons1(:,:,:)
227 28 : real(dp),allocatable :: rhog(:,:),rhor(:,:),rhowfg(:,:),rhowfr(:,:),tnons1(:,:)
228 28 : real(dp),allocatable :: vhartr(:),vpsp(:),vtrial(:,:),vxc(:,:),work(:),xccc3d(:)
229 : type(pawfgr_type) :: pawfgr
230 28 : type(pawrhoij_type),allocatable :: pawrhoij(:),pawrhoij_read(:)
231 28 : type(pawfgrtab_type),allocatable,save :: pawfgrtab(:)
232 28 : type(paw_an_type),allocatable :: paw_an(:)
233 28 : type(paw_ij_type),allocatable :: paw_ij(:)
234 28 : type(paw_dmft_type) :: paw_dmft
235 :
236 : ! ***********************************************************************
237 :
238 : DBG_ENTER("COLL")
239 :
240 28 : call timab(501,1,tsec)
241 :
242 : !Structured debugging if dtset%prtvol==-level
243 28 : if(dtset%prtvol==-level)then
244 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' nonlinear : enter , debug mode '
245 0 : call wrtout(std_out,msg)
246 : end if
247 :
248 : !Check if the perturbations asked in the input file can be computed
249 :
250 : if (((dtset%d3e_pert1_phon == 1).and.(dtset%d3e_pert2_phon == 1)).or. &
251 28 : & ((dtset%d3e_pert1_phon == 1).and.(dtset%d3e_pert3_phon == 1)).or. &
252 : & ((dtset%d3e_pert2_phon == 1).and.(dtset%d3e_pert3_phon == 1))) then
253 : write(msg,'(7a)')&
254 0 : & 'You have asked for a third-order derivative with respect to',ch10,&
255 0 : & '2 or more atomic displacements.',ch10,&
256 0 : & 'This is not allowed yet.',ch10,&
257 0 : & 'Action : change d3e_pert1_phon, d3e_pert2_phon or d3e_pert3_phon in your input file.'
258 0 : ABI_ERROR(msg)
259 : end if
260 :
261 : !Computation of third order derivatives from PEAD (pead=1) or full DPFT formalism (pead=0):
262 28 : pead = dtset%usepead
263 28 : if (pead==0) then
264 7 : write(msg, '(2a)' ) ch10,'NONLINEAR : PEAD=0, full DFPT computation of third order derivatives'
265 21 : call wrtout([ab_out, std_out], msg)
266 : end if
267 :
268 : !Some data for parallelism
269 28 : nkpt_max=50;if(xmpi_paral==1)nkpt_max=-1
270 28 : my_natom=mpi_enreg%my_natom
271 28 : paral_atom=(my_natom/=dtset%natom)
272 28 : if (paral_atom) then
273 0 : ABI_BUG(" Nonlinear routine is not available yet with parallelization over atoms...")
274 : end if
275 :
276 : !Init spaceworld
277 28 : spaceworld=mpi_enreg%comm_cell
278 28 : me = xmpi_comm_rank(spaceworld)
279 :
280 : !Define FFT grid(s) sizes (be careful !)
281 : !See NOTES in the comments at the beginning of this file.
282 28 : call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
283 :
284 28 : ntypat=psps%ntypat
285 28 : natom=dtset%natom
286 112 : nfftot=product(ngfft(1:3))
287 112 : nfftotf=product(ngfftf(1:3))
288 :
289 : !Define the set of admitted perturbations taking into account
290 : !the possible permutations
291 28 : mpert=natom+6
292 196 : ABI_MALLOC(blkflg,(3,mpert,3,mpert,3,mpert))
293 112 : ABI_MALLOC(carflg,(3,mpert,3,mpert,3,mpert))
294 112 : ABI_MALLOC(rfpert,(3,mpert,3,mpert,3,mpert))
295 84 : ABI_MALLOC(d3e_pert1,(mpert))
296 56 : ABI_MALLOC(d3e_pert2,(mpert))
297 56 : ABI_MALLOC(d3e_pert3,(mpert))
298 196 : ABI_MALLOC(d3etot,(2,3,mpert,3,mpert,3,mpert))
299 112 : ABI_MALLOC(d3cart,(2,3,mpert,3,mpert,3,mpert))
300 28 : if (pead==0) then
301 28 : ABI_MALLOC(d3e_1,(2,3,mpert,3,mpert,3,mpert))
302 28 : ABI_MALLOC(d3e_2,(2,3,mpert,3,mpert,3,mpert))
303 28 : ABI_MALLOC(d3e_3,(2,3,mpert,3,mpert,3,mpert))
304 28 : ABI_MALLOC(d3e_4,(2,3,mpert,3,mpert,3,mpert))
305 28 : ABI_MALLOC(d3e_5,(2,3,mpert,3,mpert,3,mpert))
306 28 : ABI_MALLOC(d3e_6,(2,3,mpert,3,mpert,3,mpert))
307 28 : ABI_MALLOC(d3e_7,(2,3,mpert,3,mpert,3,mpert))
308 28 : ABI_MALLOC(d3e_8,(2,3,mpert,3,mpert,3,mpert))
309 28 : ABI_MALLOC(d3e_9,(2,3,mpert,3,mpert,3,mpert))
310 328167 : d3e_1(:,:,:,:,:,:,:) = 0_dp
311 328167 : d3e_2(:,:,:,:,:,:,:) = 0_dp
312 328167 : d3e_3(:,:,:,:,:,:,:) = 0_dp
313 328167 : d3e_4(:,:,:,:,:,:,:) = 0_dp
314 328167 : d3e_5(:,:,:,:,:,:,:) = 0_dp
315 328167 : d3e_6(:,:,:,:,:,:,:) = 0_dp
316 328167 : d3e_7(:,:,:,:,:,:,:) = 0_dp
317 328167 : d3e_8(:,:,:,:,:,:,:) = 0_dp
318 328167 : d3e_9(:,:,:,:,:,:,:) = 0_dp
319 7 : if (dtset%nonlinear_info>0) then
320 4 : ABI_MALLOC(flg_tmp,(3,mpert,3,mpert,3,mpert))
321 4 : ABI_MALLOC(d3cart_1,(2,3,mpert,3,mpert,3,mpert))
322 4 : ABI_MALLOC(d3cart_2,(2,3,mpert,3,mpert,3,mpert))
323 4 : ABI_MALLOC(d3cart_3,(2,3,mpert,3,mpert,3,mpert))
324 4 : ABI_MALLOC(d3cart_4,(2,3,mpert,3,mpert,3,mpert))
325 4 : ABI_MALLOC(d3cart_5,(2,3,mpert,3,mpert,3,mpert))
326 4 : ABI_MALLOC(d3cart_6,(2,3,mpert,3,mpert,3,mpert))
327 4 : ABI_MALLOC(d3cart_7,(2,3,mpert,3,mpert,3,mpert))
328 4 : ABI_MALLOC(d3cart_8,(2,3,mpert,3,mpert,3,mpert))
329 31 : ABI_MALLOC(d3cart_9,(2,3,mpert,3,mpert,3,mpert))
330 : end if
331 : end if
332 538524 : blkflg(:,:,:,:,:,:) = 0
333 1312668 : d3etot(:,:,:,:,:,:,:) = 0_dp
334 538524 : rfpert(:,:,:,:,:,:) = 0
335 700 : d3e_pert1(:) = 0 ; d3e_pert2(:) = 0 ; d3e_pert3(:) = 0
336 :
337 74 : if (dtset%d3e_pert1_phon==1) d3e_pert1(dtset%d3e_pert1_atpol(1):dtset%d3e_pert1_atpol(2))=1
338 29 : if (dtset%d3e_pert2_phon==1) d3e_pert2(dtset%d3e_pert2_atpol(1):dtset%d3e_pert2_atpol(2))=1
339 29 : if (dtset%d3e_pert3_phon==1) d3e_pert3(dtset%d3e_pert3_atpol(1):dtset%d3e_pert3_atpol(2))=1
340 28 : if (dtset%d3e_pert1_elfd/=0) d3e_pert1(natom+2)=1
341 28 : if (dtset%d3e_pert2_elfd/=0) d3e_pert2(natom+2)=1
342 28 : if (dtset%d3e_pert3_elfd/=0) d3e_pert3(natom+2)=1
343 :
344 252 : do i1pert = 1, mpert
345 924 : do i1dir = 1, 3
346 6272 : do i2pert = 1, mpert
347 22176 : do i2dir = 1, 3
348 150528 : do i3pert = 1, mpert
349 532224 : do i3dir = 1, 3
350 : perm(1) = &
351 : & d3e_pert1(i1pert)*dtset%d3e_pert1_dir(i1dir) &
352 : & *d3e_pert2(i2pert)*dtset%d3e_pert2_dir(i2dir) &
353 387072 : & *d3e_pert3(i3pert)*dtset%d3e_pert3_dir(i3dir)
354 : perm(2) = &
355 : & d3e_pert1(i1pert)*dtset%d3e_pert1_dir(i1dir) &
356 : & *d3e_pert2(i3pert)*dtset%d3e_pert2_dir(i3dir) &
357 387072 : & *d3e_pert3(i2pert)*dtset%d3e_pert3_dir(i2dir)
358 : perm(3) = &
359 : & d3e_pert1(i2pert)*dtset%d3e_pert1_dir(i2dir) &
360 : & *d3e_pert2(i1pert)*dtset%d3e_pert2_dir(i1dir) &
361 387072 : & *d3e_pert3(i3pert)*dtset%d3e_pert3_dir(i3dir)
362 : perm(4) = &
363 : & d3e_pert1(i2pert)*dtset%d3e_pert1_dir(i2dir) &
364 : & *d3e_pert2(i3pert)*dtset%d3e_pert2_dir(i3dir) &
365 387072 : & *d3e_pert3(i1pert)*dtset%d3e_pert3_dir(i1dir)
366 : perm(5) = &
367 : & d3e_pert1(i3pert)*dtset%d3e_pert1_dir(i3dir) &
368 : & *d3e_pert2(i2pert)*dtset%d3e_pert2_dir(i2dir) &
369 387072 : & *d3e_pert3(i1pert)*dtset%d3e_pert3_dir(i1dir)
370 : perm(6) = &
371 : & d3e_pert1(i3pert)*dtset%d3e_pert1_dir(i3dir) &
372 : & *d3e_pert2(i1pert)*dtset%d3e_pert2_dir(i1dir) &
373 387072 : & *d3e_pert3(i2pert)*dtset%d3e_pert3_dir(i2dir)
374 2838528 : if (sum(perm(:)) > 0) rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
375 : end do
376 : end do
377 : end do
378 : end do
379 : end do
380 : end do
381 :
382 : ! call timab(134,2,tsec)
383 : ! call timab(135,1,tsec)
384 :
385 : !Do symmetry stuff
386 112 : ABI_MALLOC(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
387 112 : ABI_MALLOC(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
388 112 : ABI_MALLOC(indsym,(4,dtset%nsym,natom))
389 84 : ABI_MALLOC(symrec,(3,3,dtset%nsym))
390 399078 : irrzon=0;indsym=0;symrec=0;phnons=zero
391 : !If the density is to be computed by mkrho, need irrzon and phnons
392 28 : iscf_eff=0;if(dtset%getden==0)iscf_eff=1
393 : call setsym(indsym,irrzon,iscf_eff,natom,&
394 : & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
395 28 : & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred)
396 :
397 : !Symmetrize atomic coordinates over space group elements:
398 28 : call symmetrize_xred(natom,dtset%nsym,dtset%symrel,dtset%tnons,xred,indsym=indsym)
399 :
400 28 : call sytens(indsym,mpert,natom,dtset%nsym,rfpert,symrec,dtset%symrel)
401 :
402 28 : write(msg, '(a,a,a,a,a)' ) ch10, &
403 28 : & ' The list of irreducible elements of the Raman and non-linear',&
404 56 : & ch10,' optical susceptibility tensors is:',ch10
405 28 : call wrtout(ab_out,msg)
406 28 : call wrtout(std_out,msg)
407 :
408 : write(msg,'(12x,a)')&
409 28 : & 'i1pert i1dir i2pert i2dir i3pert i3dir'
410 28 : call wrtout(ab_out,msg)
411 28 : call wrtout(std_out,msg)
412 28 : n1 = 0
413 140 : do i1pert = 1, natom + 2
414 476 : do i1dir = 1, 3
415 1792 : do i2pert = 1, natom + 2
416 5712 : do i2dir = 1,3
417 21504 : do i3pert = 1, natom + 2
418 68544 : do i3dir = 1, 3
419 64512 : if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
420 698 : n1 = n1 + 1
421 698 : write(msg,'(2x,i4,a,6(5x,i3))') n1,')', &
422 1396 : & i1pert,i1dir,i2pert,i2dir,i3pert,i3dir
423 698 : call wrtout(ab_out,msg)
424 698 : call wrtout(std_out,msg)
425 47686 : else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-2) then
426 1220 : blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
427 1220 : if (dtset%nonlinear_info>0) then
428 : ! n1 = n1 + 1
429 0 : write(msg,'(2x,i4,a,6(5x,i3),a)') n1,')', &
430 0 : & i1pert,i1dir,i2pert,i2dir,i3pert,i3dir,' => must be zero, not computed'
431 0 : call wrtout(ab_out,msg)
432 0 : call wrtout(std_out,msg)
433 : end if
434 46466 : else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-1) then
435 1894 : if (dtset%nonlinear_info>0) then
436 : ! n1 = n1 + 1
437 81 : write(msg,'(2x,i4,a,6(5x,i3),a)') n1,')', &
438 162 : & i1pert,i1dir,i2pert,i2dir,i3pert,i3dir,' => symmetric of an other element, not computed'
439 81 : call wrtout(ab_out,msg)
440 81 : call wrtout(std_out,msg)
441 : end if
442 : end if
443 : end do
444 : end do
445 : end do
446 : end do
447 : end do
448 : end do
449 28 : write(msg,'(a,a)') ch10,ch10
450 28 : call wrtout(ab_out,msg)
451 28 : call wrtout(std_out,msg)
452 :
453 : ! For abipy :
454 28 : if (dtset%paral_rf == -1) then
455 0 : write(std_out,'(a)')"--- !IrredPerts"
456 0 : write(std_out,'(a)')'# List of irreducible perturbations for nonlinear'
457 0 : write(std_out,'(a)')'irred_perts:'
458 :
459 0 : n1 = 0
460 0 : do i1pert = 1, natom + 2
461 0 : do i1dir = 1, 3
462 0 : do i2pert = 1, natom + 2
463 0 : do i2dir = 1, 3
464 0 : do i3pert = 1, natom + 2
465 0 : do i3dir = 1,3
466 0 : if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
467 0 : n1 = n1 + 1
468 0 : write(std_out,'(a,i0)')" - i1pert: ",i1pert
469 0 : write(std_out,'(a,i0)')" i1dir: ",i1dir
470 0 : write(std_out,'(a,i0)')" i2pert: ",i2pert
471 0 : write(std_out,'(a,i0)')" i2dir: ",i2dir
472 0 : write(std_out,'(a,i0)')" i3pert: ",i3pert
473 0 : write(std_out,'(a,i0)')" i3dir: ",i3dir
474 : end if
475 : end do
476 : end do
477 : end do
478 : end do
479 : end do
480 : end do
481 0 : write(std_out,'(a)')"..."
482 0 : ABI_ERROR_NODUMP("aborting now")
483 : end if
484 :
485 : !Set up for iterations
486 : call setup1(dtset%acell_orig(1:3,1),bantot,dtset,&
487 : ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,&
488 : ngfftf,ngfft,dtset%nkpt,dtset%nsppol,&
489 28 : response,rmet,dtset%rprim_orig(1:3,1:3,1),rprimd,ucvol,psps%usepaw)
490 :
491 : !Set up the basis sphere of planewaves
492 84 : ABI_MALLOC(kg,(3,dtset%mpw*dtset%mkmem))
493 84 : ABI_MALLOC(npwarr,(dtset%nkpt))
494 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg,&
495 : & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw,npwarr,npwtot,&
496 28 : & dtset%nsppol)
497 :
498 : !Recompute first large sphere cut-off gsqcut, without taking into account dilatmx
499 28 : ecutf=dtset%ecut
500 28 : if (psps%usepaw==1) then
501 4 : ecutf=dtset%pawecutdg
502 4 : call wrtout(std_out,ch10//' FFT (fine) grid used in SCF cycle:')
503 : end if
504 28 : call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf)
505 :
506 : !Open and read pseudopotential files
507 : ecore = 0_dp
508 : call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,&
509 28 : & psps,rprimd,comm_mpi=mpi_enreg%comm_cell)
510 :
511 : !Initialize band structure datatype
512 28 : call bstruct%from_dtset(dtset, npwarr)
513 :
514 : !Initialize PAW atomic occupancies
515 28 : if (psps%usepaw==1) then
516 20 : ABI_MALLOC(pawrhoij,(my_natom))
517 4 : call pawrhoij_nullify(pawrhoij)
518 : call initrhoij(dtset%pawcpxocc,dtset%lexexch,dtset%lpawu, &
519 : & my_natom,natom,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,&
520 : & pawrhoij,dtset%pawspnorb,pawtab,cplex1,dtset%spinat,dtset%typat,&
521 4 : & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
522 : else
523 24 : ABI_MALLOC(pawrhoij,(0))
524 : end if
525 :
526 : !Initialize header
527 28 : gscase=0
528 : call hdr%init(bstruct,codvsn,dtset,pawtab,gscase,psps,wvl%descr, &
529 28 : & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
530 :
531 : !Update header, with evolving variables, when available
532 : !Here, rprimd, xred and occ are available
533 28 : etot=hdr%etot ; fermie=hdr%fermie ; fermih=hdr%fermih ; residm=hdr%residm
534 :
535 : !If parallelism over atom, hdr is distributed
536 : call hdr%update(bantot,etot,fermie,fermih,&
537 : residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1), &
538 28 : comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
539 :
540 : !Clean band structure datatype (should use it more in the future !)
541 28 : call bstruct%free()
542 :
543 : !Initialize wavefunction files and wavefunctions.
544 28 : ireadwf0=1
545 :
546 28 : mcg=dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
547 84 : ABI_MALLOC_OR_DIE(cg,(2,mcg), ierr)
548 :
549 84 : ABI_MALLOC(eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol))
550 5398 : eigen0(:)=zero ; ask_accurate=1
551 28 : optorth=0
552 :
553 : call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,&
554 : & formeig,hdr,ireadwf0,dtset%istwfk,kg,dtset%kptns,&
555 : & dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,&
556 : & dtset%nband,ngfft,dtset%nkpt,npwarr,dtset%nsppol,dtset%nsym,&
557 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
558 28 : & dtfil%unkg,wffgs,wfftgs,dtfil%unwffgs,dtfil%fnamewffk,wvl)
559 :
560 : !Close wffgs, if it was ever opened (in inwffil)
561 28 : if (ireadwf0==1) then
562 28 : call WffClose(wffgs,ierr)
563 : end if
564 :
565 28 : if (psps%usepaw==1.and.ireadwf0==1) then
566 : ! if parallelism, pawrhoij is distributed, hdr%pawrhoij is not
567 : call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,&
568 4 : & mpi_atmtab=mpi_enreg%my_atmtab)
569 : end if
570 :
571 : ! call timab(135,2,tsec)
572 : ! call timab(136,1,tsec)
573 :
574 : !Report on eigen0 values ! Should use prteigrs.F90
575 28 : write(msg, '(a,a)' )
576 28 : call wrtout(std_out,ch10//' respfn : eigen0 array')
577 28 : nkpt_eff=dtset%nkpt
578 28 : if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. dtset%nkpt>nkpt_max ) nkpt_eff=nkpt_max
579 28 : band_index=0
580 63 : do isppol=1,dtset%nsppol
581 1295 : do ikpt=1,dtset%nkpt
582 1232 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
583 1232 : if(ikpt<=nkpt_eff)then
584 0 : write(msg, '(a,i2,a,i5)' )' isppol=',isppol,', k point number',ikpt
585 0 : call wrtout(std_out,msg)
586 0 : do iband=1,nband_k,4
587 0 : write(msg, '(a,4es16.6)')' ',eigen0(iband+band_index:min(iband+3,nband_k)+band_index)
588 0 : call wrtout(std_out,msg)
589 : end do
590 1232 : else if(ikpt==nkpt_eff+1)then
591 0 : write(msg,'(a,a)' )' respfn : prtvol=0, 1 or 2, stop printing eigen0.',ch10
592 0 : call wrtout(std_out,msg)
593 : end if
594 1267 : band_index=band_index+nband_k
595 : end do
596 : end do
597 :
598 : !Allocation for forces and atomic positions (should be taken away, also argument ... )
599 84 : ABI_MALLOC(grxc,(3,natom))
600 :
601 : !Examine the symmetries of the q wavevector
602 84 : ABI_MALLOC(symq,(4,2,dtset%nsym))
603 : timrev=1
604 :
605 : ! By default use symmetries.
606 28 : use_sym = 1
607 28 : if (dtset%prtgkk == 1)then
608 0 : use_sym = 0
609 0 : call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol,use_sym=use_sym)
610 : else
611 28 : call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol)
612 : end if
613 :
614 :
615 :
616 : !Generate an index table of atoms, in order for them to be used
617 : !type after type.
618 84 : ABI_MALLOC(atindx,(natom))
619 56 : ABI_MALLOC(atindx1,(natom))
620 84 : ABI_MALLOC(nattyp,(ntypat))
621 77 : indx=1
622 77 : do itypat=1,ntypat
623 49 : nattyp(itypat)=0
624 175 : do iatom=1,natom
625 147 : if(dtset%typat(iatom)==itypat)then
626 56 : atindx(iatom)=indx
627 56 : atindx1(indx)=iatom
628 56 : indx=indx+1
629 56 : nattyp(itypat)=nattyp(itypat)+1
630 : end if
631 : end do
632 : end do
633 :
634 : !Compute structure factor phases for current atomic pos:
635 84 : ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*natom))
636 84 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*natom))
637 28 : call getph(atindx,natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
638 :
639 28 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
640 4 : call getph(atindx,natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
641 : else
642 15432 : ph1df(:,:)=ph1d(:,:)
643 : end if
644 :
645 28 : qeq0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2<1.d-14)
646 28 : if (.not.qeq0) then
647 0 : ABI_BUG('NONLINEAR with dtset%qptn!=0 is not implemented yet')
648 : end if
649 :
650 : !PAW: 1- Initialize values for several arrays depending only on atomic data
651 : !2- Check overlap
652 : !3- Identify FFT points in spheres and compute g_l(r).Y_lm(r) (and exp(-i.q.r) if needed)
653 : !4- Allocate PAW specific arrays
654 : !5- Compute perturbed local potential inside spheres
655 : !6- Eventually open temporary storage files
656 28 : if(psps%usepaw==1) then
657 : ! 1-Initialize values for several arrays depending only on atomic data
658 :
659 4 : gnt_option=2
660 :
661 : ! Test if we have to call pawinit
662 4 : call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
663 :
664 4 : if (psp_gencond==1.or.call_pawinit) then
665 : ! Some gen-cond have to be added...
666 0 : call timab(553,1,tsec)
667 : call pawinit(dtset%effmass_free,gnt_option,zero,zero,dtset%pawlcutd,dtset%pawlmix,&
668 : & psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,&
669 0 : & pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero)
670 : call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,&
671 0 : & rprimd,symrec,pawang%zarot)
672 :
673 : ! Update internal values
674 0 : call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
675 :
676 0 : call timab(553,2,tsec)
677 : else
678 4 : if (pawtab(1)%has_kij ==1) pawtab(1:psps%ntypat)%has_kij =2
679 4 : if (pawtab(1)%has_nabla==1) pawtab(1:psps%ntypat)%has_nabla=2
680 : end if
681 12 : psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore)
682 4 : call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot)
683 : call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,&
684 : & is_dfpt,dtset%jpawu,dtset%lexexch,dtset%lpawu,dtset%nspinor,ntypat,dtset%optdcmagpawu,pawang,dtset%pawprtvol,pawrad,&
685 4 : & pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu)
686 4 : compch_fft=-1.d5;compch_sph=-1.d5
687 12 : usexcnhat=maxval(pawtab(:)%usexcnhat)
688 :
689 : ! Note: many derivatives of cprj are needed and used a few times only, so for simplicity the
690 : ! computation of all needed derivatives will be done on-the-fly.
691 4 : usecprj=0
692 :
693 : ! 2-Check overlap
694 4 : call chkpawovlp(natom,psps%ntypat,dtset%pawovlp,pawtab,rmet,dtset%typat,xred)
695 : ! 3-Identify FFT points in spheres and compute g_l(r).Y_lm(r) and exp(-i.q.r)
696 20 : ABI_MALLOC(pawfgrtab,(my_natom))
697 4 : if (my_natom>0) then
698 : call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,&
699 4 : & mpi_atmtab=mpi_enreg%my_atmtab)
700 : call pawfgrtab_init(pawfgrtab,1,l_size_atm,pawrhoij(1)%nspden,dtset%typat,&
701 4 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
702 4 : ABI_FREE(l_size_atm)
703 : end if
704 4 : optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm
705 4 : optgr1=dtset%pawstgylm
706 4 : optgr2=dtset%pawstgylm
707 : call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfftf,psps%ntypat,&
708 : & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,&
709 4 : & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab )
710 20 : ABI_MALLOC(paw_an,(my_natom))
711 20 : ABI_MALLOC(paw_ij,(my_natom))
712 4 : call paw_an_nullify(paw_an)
713 4 : call paw_ij_nullify(paw_ij)
714 4 : has_kxc=0;nkxc1=0;cplex=1
715 36 : has_dijnd=0;if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1
716 4 : has_diju=merge(0,1,dtset%usepawu==0)
717 4 : has_kxc=1;nkxc1=2*dtset%nspden-1 ! LDA only
718 4 : call pawxc_get_nkxc(nkxc1,dtset%nspden,dtset%xclevel)
719 4 : has_k3xc=1; nk3xc1=3*min(dtset%nspden,2)-2 ! LDA only
720 : call paw_an_init(paw_an,dtset%natom,dtset%ntypat,nkxc1,nk3xc1,dtset%nspden,cplex,dtset%pawxcdev,&
721 : & dtset%typat,pawang,pawtab,has_vxc=1,has_vxc_ex=1,has_kxc=has_kxc,has_k3xc=has_k3xc,&
722 4 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
723 : call paw_ij_init(paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%pawspnorb,&
724 : & natom,dtset%ntypat,dtset%typat,pawtab,has_dij=1,has_dijhartree=1,has_dijnd=has_dijnd,&
725 : & has_dijso=1,has_dijU=has_diju,has_pawu_occ=1,has_exexch_pot=1,nucdipmom=dtset%nucdipmom,&
726 12 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
727 : else ! PAW vs NCPP
728 24 : usexcnhat=0;usecprj=0
729 24 : ABI_MALLOC(paw_an,(0))
730 24 : ABI_MALLOC(paw_ij,(0))
731 24 : ABI_MALLOC(pawfgrtab,(0))
732 : end if
733 :
734 84 : ABI_MALLOC(rhog,(2,nfftf))
735 112 : ABI_MALLOC(rhor,(nfftf,dtset%nspden))
736 :
737 : !Read ground-state charge density from diskfile in case getden /= 0
738 : !or compute it from wfs that were read previously : rhor as well as rhog
739 :
740 28 : if (dtset%getden /= 0 .or. dtset%irdden /= 0) then
741 : ! Read rho1(r) from a disk file and broadcast data.
742 : ! This part is not compatible with MPI-FFT (note single_proc=.True. below)
743 :
744 26 : rdwr=1;rdwrpaw=psps%usepaw;if(ireadwf0/=0) rdwrpaw=0
745 : if (rdwrpaw/=0) then
746 : ABI_MALLOC(pawrhoij_read,(natom))
747 : call pawrhoij_nullify(pawrhoij_read)
748 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,&
749 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc)
750 : call pawrhoij_alloc(pawrhoij_read,cplex_rhoij,nspden_rhoij,dtset%nspinor,&
751 : & dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab)
752 : else
753 26 : ABI_MALLOC(pawrhoij_read,(0))
754 : end if
755 :
756 : ! MT july 2013: Should we read rhoij from the density file ?
757 : call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor, &
758 26 : hdr_den, pawrhoij_read, spaceworld, check_hdr=hdr)
759 26 : call hdr_den%free()
760 :
761 : if (rdwrpaw/=0) then
762 : call pawrhoij_bcast(pawrhoij_read,hdr%pawrhoij,0,spaceworld)
763 : call pawrhoij_free(pawrhoij_read)
764 : end if
765 26 : ABI_FREE(pawrhoij_read)
766 :
767 : ! Compute up+down rho(G) by fft
768 78 : ABI_MALLOC(work,(nfftf))
769 163138 : work(:)=rhor(:,1)
770 26 : call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,1,ngfftf,0)
771 26 : ABI_FREE(work)
772 :
773 : else
774 2 : izero=0
775 : ! Obtain the charge density from read wfs
776 : ! Be careful: in PAW, compensation density has to be added !
777 2 : tim_mkrho=4
778 2 : paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented
779 2 : paw_dmft%use_dmft=0 ! respfn with dmft not implemented
780 2 : if (psps%usepaw==1) then
781 0 : ABI_MALLOC(rhowfg,(2,dtset%nfft))
782 0 : ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden))
783 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
784 0 : & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
785 0 : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor)
786 0 : ABI_FREE(rhowfg)
787 0 : ABI_FREE(rhowfr)
788 : else
789 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
790 2 : & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
791 : end if
792 : end if ! getden
793 :
794 : !In PAW, compensation density has eventually to be added
795 28 : nhatgrdim=0;nhatdim=0
796 28 : ABI_MALLOC(nhatgr,(0,0,0))
797 28 : if (psps%usepaw==1.and. ((usexcnhat==0).or.(dtset%getden==0).or.dtset%xclevel==2)) then
798 4 : nhatdim=1
799 16 : ABI_MALLOC(nhat,(nfftf,dtset%nspden))
800 4 : call timab(558,1,tsec)
801 4 : nhatgrdim=0;if (dtset%xclevel==2.and.dtset%pawnhatxc>0) nhatgrdim=usexcnhat
802 4 : ider=2*nhatgrdim
803 4 : if (nhatgrdim>0) then
804 0 : ABI_FREE(nhatgr)
805 0 : ABI_MALLOC(nhatgr,(nfftf,dtset%nspden,3))
806 : end if
807 4 : izero=0;cplex=1;ipert=0;idir=0;qphon(:)=zero
808 : call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,natom,&
809 : & nfftf,ngfftf,nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,&
810 : & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, &
811 4 : & mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom)
812 4 : if (dtset%getden==0) then
813 0 : rhor(:,:)=rhor(:,:)+nhat(:,:)
814 0 : call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0)
815 : end if
816 4 : call timab(558,2,tsec)
817 : else
818 24 : ABI_MALLOC(nhat,(0,0))
819 : end if
820 :
821 : !The GS irrzon and phnons were only needed to symmetrize the GS density
822 28 : ABI_FREE(irrzon)
823 28 : ABI_FREE(phnons)
824 :
825 : !!jmb 2012 write(std_out,'(a)')' ' ! needed to make ibm6_xlf12 pass tests. No idea why this works. JWZ 5 Sept 2011
826 : !!Will compute now the total potential
827 :
828 : !Compute local ionic pseudopotential vpsp and core electron density xccc3d:
829 28 : n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf
830 84 : ABI_MALLOC(xccc3d,(n3xccc))
831 84 : ABI_MALLOC(vpsp,(nfftf))
832 :
833 28 : eei = zero
834 28 : if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
835 : ! PAW or NC with nc_xccc_gspace: compute Vloc and core charge together in reciprocal space
836 4 : call timab(562,1,tsec)
837 4 : optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optv=1;optn=n3xccc/nfftf;optn2=1
838 : call atm2fft(atindx1,xccc3d,vpsp,dum_dyfrn,dum_dyfrv,dum_eltfrxc,dum_gauss,gmet,gprimd,&
839 : & dum_grn,dum_grv,gsqcut,mgfftf,psps%mqgrid_vl,natom,nattyp,nfftf,ngfftf,&
840 : & ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1df,psps%qgrid_vl,&
841 4 : & dtset%qprtrb,dtset%rcut,dum_rhog,rprimd,dummy6,other_dummy6,ucvol,psps%usepaw,dum_vg,dum_vg,dum_vg,dtset%vprtrb,psps%vlspl)
842 4 : call timab(562,2,tsec)
843 : else
844 : ! Norm-cons.: compute Vloc in reciprocal space and core charge in real space
845 24 : option=1
846 72 : ABI_MALLOC(dyfrlo_indx,(3,3,natom))
847 72 : ABI_MALLOC(grtn_indx,(3,natom))
848 : call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,&
849 : & grtn_indx,gsqcut,dummy6,mgfftf,mpi_enreg,natom,nattyp,&
850 : & nfftf,ngfftf,dtset%nspden,ntypat,option,pawtab,ph1df,psps,&
851 24 : & dtset%qprtrb,rhog,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred)
852 24 : ABI_FREE(dyfrlo_indx)
853 24 : ABI_FREE(grtn_indx)
854 24 : if (psps%n1xccc/=0) then
855 18 : ABI_MALLOC(dyfrx2,(3,3,natom))
856 9 : ABI_MALLOC(vxc,(0,0)) ! dummy
857 : call mkcore(dummy6,dyfrx2,grxc,mpi_enreg,natom,nfftf,dtset%nspden,ntypat,&
858 : & ngfftf(1),psps%n1xccc,ngfftf(2),ngfftf(3),option,rprimd,dtset%typat,ucvol,vxc,&
859 9 : & psps%xcccrc,psps%xccc1d,xccc3d,xred)
860 9 : ABI_FREE(dyfrx2)
861 9 : ABI_FREE(vxc) ! dummy
862 : end if
863 : end if
864 :
865 : !Set up hartree and xc potential. Compute kxc here.
866 56 : ABI_MALLOC(vhartr,(nfftf))
867 : call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfftf,ngfftf,&
868 28 : &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr)
869 :
870 28 : option=3
871 28 : nkxc=2*dtset%nspden-1 ! LDA
872 28 : if(dtset%xclevel==2.and.dtset%nspden==1) nkxc=7 ! non-polarized GGA
873 28 : if(dtset%xclevel==2.and.dtset%nspden==2) nkxc=19 ! polarized GGA
874 28 : nk3xc=3*dtset%nspden-2
875 28 : call check_kxc(dtset%ixc,dtset%optdriver,check_k3xc=.true.)
876 112 : ABI_MALLOC(kxc,(nfftf,nkxc))
877 112 : ABI_MALLOC(k3xc,(nfftf,nk3xc))
878 112 : ABI_MALLOC(vxc,(nfftf,dtset%nspden))
879 :
880 28 : call xcdata_init(xcdata,dtset=dtset)
881 28 : nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
882 : call rhotoxc(bigexc,bigsxc,kxc,mpi_enreg,nfftf,ngfftf,&
883 : & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc,option,rhor,&
884 28 : & rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata,k3xc=k3xc,vhartr=vhartr)
885 :
886 : !Get electronic temperature from dtset
887 28 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
888 :
889 : !Compute local + Hxc potential, and subtract mean potential.
890 112 : ABI_MALLOC(vtrial,(nfftf,dtset%nspden))
891 63 : do ispden=1,min(dtset%nspden,2)
892 215999 : do ifft=1,nfftf
893 215971 : vtrial(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft)
894 : end do
895 : end do
896 28 : if (dtset%nspden==4) then
897 0 : do ispden=3,4
898 0 : do ifft=1,nfftf
899 0 : vtrial(ifft,ispden)=vxc(ifft,ispden)
900 : end do
901 : end do
902 : end if
903 28 : ABI_FREE(vpsp)
904 28 : ABI_FREE(vhartr)
905 :
906 28 : if(dtset%prtvol==-level)then
907 0 : call wrtout(std_out,' nonlinear : ground-state density and potential set up.')
908 : end if
909 :
910 : !PAW: compute Dij quantities (psp strengths)
911 28 : if (psps%usepaw==1)then
912 4 : cplex=1;ipert=0;option=1
913 4 : nzlmopt=0;if (dtset%pawnzlm>0) nzlmopt=-1
914 : call pawdenpot(compch_sph,el_temp,gprimd,ipert,dtset%ixc,my_natom,natom,dtset%nspden,&
915 : & ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_energies,paw_ij,pawang,&
916 : & dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,&
917 : & dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,xred,ucvol,psps%znuclpsp,dtset%spinaxis, &
918 4 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
919 :
920 4 : call timab(561,1,tsec)
921 : call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,natom,nfftf,nfftotf,&
922 : & dtset%nspden,ntypat,paw_an,paw_ij,pawang,pawfgrtab,dtset%pawprtvol,&
923 : & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,k0,&
924 : & dtset%spnorbscl,ucvol,dtset%cellcharge(1),vtrial,vxc,xred,dtset%znucl,&
925 : & nucdipmom=dtset%nucdipmom,&
926 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,&
927 4 : & spinaxis=dtset%spinaxis)
928 : call symdij(gprimd,indsym,ipert,my_natom,natom,dtset%nsym,ntypat,0,&
929 : & paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,&
930 4 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
931 4 : call timab(561,2,tsec)
932 : end if
933 :
934 28 : ABI_FREE(xccc3d)
935 :
936 : ! Determine the subset of symmetry operations (nsym1 operations)
937 : ! that leaves the perturbation invariant, and initialize corresponding arrays
938 : ! symaf1, symrl1, tnons1 (and pawang1%zarot, if PAW)..
939 28 : nsym1 = 1
940 : ! symaf1_tmp(1) = 1
941 : ! symrl1_tmp(:,:,1) = dtset%symrel(:,:,1)
942 : ! tnons1_tmp(:,1) = 0_dp
943 84 : ABI_MALLOC(indsy1,(4,nsym1,dtset%natom))
944 28 : ABI_MALLOC(symrc1,(3,3,nsym1))
945 28 : ABI_MALLOC(symaf1,(nsym1))
946 28 : ABI_MALLOC(symrl1,(3,3,nsym1))
947 28 : ABI_MALLOC(tnons1,(3,nsym1))
948 28 : symaf1(1)= 1 !symaf1_tmp(1:nsym1)
949 364 : symrl1(:,:,1)= dtset%symrel(:,:,1) !symrl1_tmp(:,:,1:nsym1)
950 112 : tnons1(:,1)= 0_dp !tnons1_tmp(:,1:nsym1)
951 : ! ABI_FREE(symaf1_tmp)
952 : ! ABI_FREE(symrl1_tmp)
953 : ! ABI_FREE(tnons1_tmp)
954 :
955 : ! Set up corresponding symmetry data
956 112 : ABI_MALLOC(irrzon1,(dtset%nfft**(1-1/nsym1),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
957 112 : ABI_MALLOC(phnons1,(2,dtset%nfft**(1-1/nsym1),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
958 : call setsym(indsy1,irrzon1,1,dtset%natom,dtset%nfft,dtset%ngfft,dtset%nspden,dtset%nsppol,&
959 28 : & nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred)
960 28 : if (psps%usepaw==1) then
961 : ! Allocate/initialize only zarot in pawang1 datastructure
962 4 : call pawang_init(pawang1,0,0,pawang%l_max-1,0,0,nsym1,0,0,0,0)
963 4 : call setsym_ylm(gprimd,pawang1%l_max-1,pawang1%nsym,0,rprimd,symrc1,pawang1%zarot)
964 : end if
965 :
966 28 : if (pead/=0) then
967 : ! Initialize finite difference calculation of the ddk
968 :
969 21 : nkpt3 = 0
970 :
971 : ! Prepare first call to getkgrid (obtain number of k points in FBZ)
972 273 : dum_kptrlatt(:,:) = dtset%kptrlatt(:,:)
973 21 : dum_nshiftk = dtset%nshiftk
974 21 : ABI_CHECK(dum_nshiftk <= MAX_NSHIFTK, sjoin("dum_nshiftk must be <= ", itoa(MAX_NSHIFTK)))
975 21 : dum_shiftk(:,:) = zero
976 105 : dum_shiftk(:,1:dtset%nshiftk) = dtset%shiftk(:,1:dtset%nshiftk)
977 21 : dum_vacuum(:) = 0
978 :
979 21 : ABI_MALLOC(dum_kptns,(3,0))
980 21 : ABI_MALLOC(dum_wtk,(0))
981 : call getkgrid(0,0,dtset%iscf,dum_kptns,3,dum_kptrlatt,&
982 : & rdum,dtset%nsym,0,nkpt3,dum_nshiftk,dtset%nsym,&
983 : & rprimd,dum_shiftk,dtset%symafm,dtset%symrel,&
984 21 : & dum_vacuum,dum_wtk)
985 21 : ABI_FREE(dum_kptns)
986 21 : ABI_FREE(dum_wtk)
987 :
988 : ! write(std_out,*) 'nonlinear : nkpt, nkpt3 = ',dtset%nkpt,nkpt3
989 : !call flush(6)
990 : !jmb : malloc() problem with gcc461_openmpi under max2 : change order of allocations works ?!?
991 : !allocate(kneigh(30,nkpt),kg_neigh(30,nkpt,3),mvwtk(30,nkpt))
992 84 : ABI_MALLOC(kg_neigh,(30,dtset%nkpt,3))
993 63 : ABI_MALLOC(mvwtk,(30,dtset%nkpt))
994 63 : ABI_MALLOC(kneigh,(30,dtset%nkpt))
995 :
996 63 : ABI_MALLOC(kptindex,(2,nkpt3))
997 63 : ABI_MALLOC(kpt3,(3,nkpt3))
998 :
999 : call getshell(gmet,kneigh,kg_neigh,kptindex,dtset%kptopt,&
1000 : & dtset%kptrlatt,dtset%kptns,kpt3,dtset%mkmem,mkmem_max,mvwtk,&
1001 21 : & dtset%nkpt,nkpt3,nneigh,dtset%nshiftk,rmet,rprimd,dtset%shiftk,dtset%wtk, mpi_enreg%comm_cell)
1002 :
1003 105 : ABI_MALLOC(pwind,(dtset%mpw,nneigh,dtset%mkmem))
1004 84 : ABI_MALLOC(cgindex,(dtset%nkpt,dtset%nsppol))
1005 105 : ABI_MALLOC(mpi_enreg%kpt_loc2ibz_sp,(0:mpi_enreg%nproc-1,1:mkmem_max, 1:2))
1006 63 : ABI_MALLOC(mpi_enreg%mkmem,(0:mpi_enreg%nproc-1))
1007 :
1008 : call initmv(cgindex,dtset,gmet,kg,kneigh,kg_neigh,kptindex,&
1009 : & kpt3,dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nband,dtset%nkpt,&
1010 21 : & nkpt3,nneigh,npwarr,dtset%nsppol,occ,pwind)
1011 :
1012 : call pead_nl_loop(blkflg,cg,cgindex,dtfil,dtset,d3etot,gmet,gprimd,gsqcut,&
1013 : & hdr,kg,kneigh,kg_neigh,kptindex,kpt3,kxc,k3xc,dtset%mband,dtset%mgfft,&
1014 : & dtset%mkmem,mkmem_max,dtset%mk1mem,mpert,mpi_enreg,dtset%mpw,mvwtk,natom,nfftf,&
1015 : & dtset%nkpt,nkpt3,nkxc,nk3xc,nneigh,dtset%nspinor,dtset%nsppol,npwarr,occ,psps,pwind,&
1016 21 : & rfpert,rprimd,ucvol,xred)
1017 :
1018 : else ! pead=0 in this case
1019 :
1020 : call dfptnl_loop(atindx,blkflg,cg,dtfil,dtset,d3etot,eigen0,gmet,gprimd,gsqcut,&
1021 : & hdr,kg,kxc,k3xc,dtset%mband,dtset%mgfft,mgfftf,&
1022 : & dtset%mkmem,dtset%mk1mem,mpert,mpi_enreg,dtset%mpw,natom,nattyp,ngfftf,nfftf,nhat,&
1023 : & dtset%nkpt,nkxc,nk3xc,dtset%nspinor,dtset%nsppol,npwarr,occ,&
1024 : & paw_an,paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
1025 : & ph1d,ph1df,psps,rfpert,rhog,rhor,rprimd,ucvol,usecprj,vtrial,vxc,xred,&
1026 : & nsym1,indsy1,symaf1,symrc1,&
1027 7 : & d3e_1,d3e_2,d3e_3,d3e_4,d3e_5,d3e_6,d3e_7,d3e_8,d3e_9)
1028 :
1029 : !Complete missing elements using symmetry operations
1030 :
1031 7 : if (dtset%nonlinear_info>0) then
1032 19234 : flg_tmp = blkflg
1033 1 : call d3sym(flg_tmp,d3e_1,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1034 19234 : flg_tmp = blkflg
1035 1 : call d3sym(flg_tmp,d3e_2,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1036 19234 : flg_tmp = blkflg
1037 1 : call d3sym(flg_tmp,d3e_3,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1038 19234 : flg_tmp = blkflg
1039 1 : call d3sym(flg_tmp,d3e_4,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1040 19234 : flg_tmp = blkflg
1041 1 : call d3sym(flg_tmp,d3e_5,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1042 19234 : flg_tmp = blkflg
1043 1 : call d3sym(flg_tmp,d3e_6,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1044 19234 : flg_tmp = blkflg
1045 1 : call d3sym(flg_tmp,d3e_7,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1046 19234 : flg_tmp = blkflg
1047 1 : call d3sym(flg_tmp,d3e_8,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1048 19234 : flg_tmp = blkflg
1049 1 : call d3sym(flg_tmp,d3e_9,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1050 : end if
1051 :
1052 : end if ! end pead/=0
1053 :
1054 28 : write(msg,'(a,a,a)')ch10,&
1055 56 : & ' --- Third order energy calculation completed --- ',ch10
1056 28 : call wrtout(ab_out,msg)
1057 :
1058 :
1059 : !Complete missing elements using symmetry operations
1060 28 : call d3sym(blkflg,d3etot,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1061 :
1062 28 : if (mpi_enreg%me == 0) then
1063 :
1064 : ! Write 3rd order derivatives in the output file
1065 24 : call dfptnl_doutput(blkflg,d3etot,mpert)
1066 :
1067 : ! Write the DDB file
1068 24 : dscrpt=' Note : temporary (transfer) database '
1069 24 : call ddb_hdr%init(dtset,psps,pawtab,dscrpt,1,xred=xred,occ=occ)
1070 :
1071 24 : call ddb%init(dtset, 1, mpert, with_d3E=.true.)
1072 :
1073 24 : call ddb%set_d3matr(1, d3etot, blkflg)
1074 :
1075 24 : call ddb%write(ddb_hdr, dtfil%fnameabo_ddb)
1076 :
1077 24 : call ddb_hdr%free()
1078 24 : call ddb%free()
1079 :
1080 : ! Compute tensors related to third-order derivatives
1081 24 : call nlopt(blkflg,carflg,d3etot,d3cart,gprimd,mpert,natom,rprimd,ucvol)
1082 : ! Note that the imaginary part is not transformed into cartesian coordinates
1083 24 : if (pead==0.and.(dtset%nonlinear_info>0)) then
1084 1 : call nlopt(blkflg,flg_tmp,d3e_1,d3cart_1,gprimd,mpert,natom,rprimd,ucvol)
1085 1 : call nlopt(blkflg,flg_tmp,d3e_2,d3cart_2,gprimd,mpert,natom,rprimd,ucvol)
1086 1 : call nlopt(blkflg,flg_tmp,d3e_3,d3cart_3,gprimd,mpert,natom,rprimd,ucvol)
1087 1 : call nlopt(blkflg,flg_tmp,d3e_4,d3cart_4,gprimd,mpert,natom,rprimd,ucvol)
1088 1 : call nlopt(blkflg,flg_tmp,d3e_5,d3cart_5,gprimd,mpert,natom,rprimd,ucvol)
1089 1 : call nlopt(blkflg,flg_tmp,d3e_6,d3cart_6,gprimd,mpert,natom,rprimd,ucvol)
1090 1 : call nlopt(blkflg,flg_tmp,d3e_7,d3cart_7,gprimd,mpert,natom,rprimd,ucvol)
1091 1 : call nlopt(blkflg,flg_tmp,d3e_8,d3cart_8,gprimd,mpert,natom,rprimd,ucvol)
1092 1 : call nlopt(blkflg,flg_tmp,d3e_9,d3cart_9,gprimd,mpert,natom,rprimd,ucvol)
1093 : end if
1094 :
1095 24 : if ((d3e_pert1(natom+2)==1).and.(d3e_pert2(natom+2)==1).and. &
1096 : & (d3e_pert3(natom+2)==1)) then
1097 :
1098 24 : flag = 1
1099 24 : i1pert = natom+2
1100 :
1101 : d3cart(:,:,i1pert,:,i1pert,:,i1pert) = &
1102 2256 : & d3cart(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
1103 :
1104 24 : write(ab_out,*)ch10
1105 24 : write(ab_out,*)' Non-linear optical susceptibility tensor d (pm/V)'
1106 24 : write(ab_out,*)' in cartesian coordinates'
1107 24 : write(ab_out,*)' i1dir i2dir i3dir d'
1108 :
1109 96 : do i1dir = 1, 3
1110 312 : do i2dir = 1, 3
1111 936 : do i3dir = 1, 3
1112 648 : write(ab_out,'(3(5x,i2),5x,f16.9)') i1dir,i2dir,i3dir,&
1113 1296 : & d3cart(1,i1dir,i1pert,i2dir,i1pert,i3dir,i1pert)
1114 648 : if ((blkflg(i1dir,i1pert,i2dir,i1pert,i3dir,i1pert)/=1).or.&
1115 427 : & (carflg(i1dir,i1pert,i2dir,i1pert,i3dir,i1pert)/=1)) flag = 0
1116 : end do
1117 : end do
1118 : end do
1119 :
1120 24 : if (pead==0.and.(dtset%nonlinear_info>0)) then
1121 :
1122 : d3cart_1(:,:,i1pert,:,i1pert,:,i1pert) = &
1123 94 : & d3cart_1(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
1124 : d3cart_2(:,:,i1pert,:,i1pert,:,i1pert) = &
1125 94 : & d3cart_2(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
1126 : d3cart_8(:,:,i1pert,:,i1pert,:,i1pert) = &
1127 94 : & d3cart_8(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
1128 : d3cart_9(:,:,i1pert,:,i1pert,:,i1pert) = &
1129 94 : & d3cart_9(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
1130 :
1131 1 : theunit = ab_out
1132 :
1133 1 : write(small_msg,'(a)') ' ** Total :'
1134 1 : call print_chi2(d3cart,small_msg,theunit)
1135 :
1136 1 : write(small_msg,'(a)') ' ** sum_psi1H1psi1 :'
1137 1 : call print_chi2(d3cart_1,small_msg,theunit)
1138 :
1139 1 : write(small_msg,'(a)') ' ** sum_lambda1psi1psi1 :'
1140 1 : call print_chi2(d3cart_2,small_msg,theunit)
1141 :
1142 1 : write(small_msg,'(a)') ' ** exc3 :'
1143 1 : call print_chi2(d3cart_8,small_msg,theunit)
1144 :
1145 1 : write(small_msg,'(a)') ' ** exc3_paw :'
1146 1 : call print_chi2(d3cart_9,small_msg,theunit)
1147 :
1148 : end if ! nonlinear_info > 0
1149 :
1150 24 : if (flag == 0) then
1151 8 : write(msg,'(a,a,a,a,a,a)')ch10,&
1152 8 : & ' dfptnl_doutput: WARNING -',ch10,&
1153 8 : & ' matrix of third-order energies incomplete,',ch10,&
1154 16 : & ' non-linear optical coefficients may be wrong, check input variables rfatpol and rfdir.'
1155 8 : call wrtout(ab_out,msg)
1156 8 : call wrtout(std_out,msg)
1157 : end if
1158 :
1159 : end if ! d3e_pert1,d3e_pert2,d3e_pert3
1160 :
1161 : if (((maxval(d3e_pert1(1:natom))/=0).and.(d3e_pert2(natom+2)/=0).and. &
1162 : & (d3e_pert3(natom+2)/=0)).or.&
1163 : ((maxval(d3e_pert2(1:natom))/=0).and.(d3e_pert1(natom+2)/=0).and. &
1164 168 : & (d3e_pert3(natom+2)/=0)).or.&
1165 : ((maxval(d3e_pert3(1:natom))/=0).and.(d3e_pert2(natom+2)/=0).and. &
1166 : & (d3e_pert1(natom+2)/=0))) then
1167 : ! Perform a check if all relevant elements are available
1168 :
1169 24 : flag = 1
1170 72 : do i1pert = 1, natom
1171 216 : do i1dir = 1, 3
1172 624 : do i2dir = 1, 3
1173 1872 : do i3dir = 1, 3
1174 : if ((blkflg(i1dir,i1pert,i2dir,natom+2,i3dir,natom+2) /= 1).or.&
1175 1296 : (blkflg(i1dir,natom+2,i2dir,i1pert,i3dir,natom+2) /= 1).or.&
1176 418 : (blkflg(i1dir,natom+2,i2dir,natom+2,i3dir,i1pert) /= 1)) flag = 0
1177 : if ((carflg(i1dir,i1pert,i2dir,natom+2,i3dir,natom+2) /= 1).or.&
1178 1296 : (carflg(i1dir,natom+2,i2dir,i1pert,i3dir,natom+2) /= 1).or.&
1179 859 : (carflg(i1dir,natom+2,i2dir,natom+2,i3dir,i1pert) /= 1)) flag = 0
1180 : end do
1181 : end do
1182 : end do
1183 : end do
1184 :
1185 24 : write(ab_out,*)ch10
1186 24 : write(ab_out,*)' First-order change in the electronic dielectric '
1187 24 : write(ab_out,*)' susceptibility tensor (Bohr^-1)'
1188 24 : write(ab_out,*)' induced by an atomic displacement'
1189 24 : write(ab_out,*)' atom displacement'
1190 :
1191 72 : do i1pert = 1,natom
1192 192 : do i1dir = 1,3
1193 144 : write(ab_out,'(1x,i4,9x,i2,3(3x,f16.9))')i1pert,i1dir,&
1194 288 : & d3cart(1,i1dir,i1pert,1,natom+2,:,natom+2)
1195 : write(ab_out,'(16x,3(3x,f16.9))')&
1196 144 : & d3cart(1,i1dir,i1pert,2,natom+2,:,natom+2)
1197 : write(ab_out,'(16x,3(3x,f16.9))')&
1198 192 : & d3cart(1,i1dir,i1pert,3,natom+2,:,natom+2)
1199 : end do
1200 72 : write(ab_out,*)
1201 : end do
1202 :
1203 24 : if (flag == 0) then
1204 8 : write(msg,'(a,a,a,a,a,a)')ch10,&
1205 8 : & ' dfptnl_doutput: WARNING -',ch10,&
1206 8 : & ' matrix of third-order energies incomplete,',ch10,&
1207 16 : & ' changes in the dielectric susceptibility may be wrong, check input variables rfatpol and rfdir.'
1208 8 : call wrtout(ab_out,msg)
1209 8 : call wrtout(std_out,msg)
1210 : end if
1211 :
1212 24 : if (pead==0.and.(dtset%nonlinear_info>0)) then
1213 1 : theunit = ab_out
1214 :
1215 1 : write(small_msg,'(a)') ' ** Total :'
1216 1 : call print_dchidtau(d3cart,small_msg,theunit)
1217 :
1218 1 : write(small_msg,'(a)') ' ** sum_psi1H1psi1 :'
1219 1 : call print_dchidtau(d3cart_1,small_msg,theunit)
1220 :
1221 1 : write(small_msg,'(a)') ' ** sum_lambda1psi1psi1 :'
1222 1 : call print_dchidtau(d3cart_2,small_msg,theunit)
1223 :
1224 1 : write(small_msg,'(a)') ' ** sum_lambda1psi0S1psi1 :'
1225 1 : call print_dchidtau(d3cart_3,small_msg,theunit)
1226 :
1227 1 : write(small_msg,'(a)') ' ** sum_psi0H2psi1a :'
1228 1 : call print_dchidtau(d3cart_4,small_msg,theunit)
1229 :
1230 1 : write(small_msg,'(a)') ' ** sum_psi0H2psi1b :'
1231 1 : call print_dchidtau(d3cart_5,small_msg,theunit)
1232 :
1233 1 : write(small_msg,'(a)') ' ** eHxc21_paw :'
1234 1 : call print_dchidtau(d3cart_6,small_msg,theunit)
1235 :
1236 1 : write(small_msg,'(a)') ' ** eHxc21_nhat :'
1237 1 : call print_dchidtau(d3cart_7,small_msg,theunit)
1238 :
1239 1 : write(small_msg,'(a)') ' ** exc3 :'
1240 1 : call print_dchidtau(d3cart_8,small_msg,theunit)
1241 :
1242 1 : write(small_msg,'(a)') ' ** exc3_paw :'
1243 1 : call print_dchidtau(d3cart_9,small_msg,theunit)
1244 :
1245 : end if ! nonlinear_info > 0
1246 :
1247 : end if ! d3e_pert1,d3e_pert2,d3e_pert3
1248 : end if ! mpi_enreg%me
1249 :
1250 : ! TO OPTIMIZE DEALLOCATION !
1251 11 : if (pead/=0) then
1252 21 : ABI_FREE(cgindex)
1253 21 : ABI_FREE(kg_neigh)
1254 21 : ABI_FREE(kneigh)
1255 21 : ABI_FREE(kptindex)
1256 21 : ABI_FREE(kpt3)
1257 21 : ABI_FREE(mpi_enreg%kpt_loc2ibz_sp)
1258 21 : ABI_FREE(mpi_enreg%mkmem)
1259 21 : ABI_FREE(mvwtk)
1260 21 : ABI_FREE(pwind)
1261 : else
1262 7 : if (dtset%nonlinear_info>0) then
1263 1 : ABI_FREE(d3cart_1)
1264 1 : ABI_FREE(d3cart_2)
1265 1 : ABI_FREE(d3cart_3)
1266 1 : ABI_FREE(d3cart_4)
1267 1 : ABI_FREE(d3cart_5)
1268 1 : ABI_FREE(d3cart_6)
1269 1 : ABI_FREE(d3cart_7)
1270 1 : ABI_FREE(d3cart_8)
1271 1 : ABI_FREE(d3cart_9)
1272 1 : ABI_FREE(flg_tmp)
1273 : end if
1274 7 : ABI_FREE(d3e_1)
1275 7 : ABI_FREE(d3e_2)
1276 7 : ABI_FREE(d3e_3)
1277 7 : ABI_FREE(d3e_4)
1278 7 : ABI_FREE(d3e_5)
1279 7 : ABI_FREE(d3e_6)
1280 7 : ABI_FREE(d3e_7)
1281 7 : ABI_FREE(d3e_8)
1282 7 : ABI_FREE(d3e_9)
1283 : end if
1284 28 : ABI_FREE(atindx)
1285 28 : ABI_FREE(atindx1)
1286 28 : ABI_FREE(blkflg)
1287 28 : ABI_FREE(carflg)
1288 28 : ABI_FREE(cg)
1289 28 : ABI_FREE(d3cart)
1290 28 : ABI_FREE(d3etot)
1291 28 : ABI_FREE(d3e_pert1)
1292 28 : ABI_FREE(d3e_pert2)
1293 28 : ABI_FREE(d3e_pert3)
1294 28 : ABI_FREE(eigen0)
1295 28 : ABI_FREE(rhog)
1296 28 : ABI_FREE(rhor)
1297 28 : ABI_FREE(nhat)
1298 28 : ABI_FREE(nhatgr)
1299 28 : ABI_FREE(rfpert)
1300 28 : ABI_FREE(grxc)
1301 28 : ABI_FREE(kg)
1302 28 : ABI_FREE(kxc)
1303 28 : ABI_FREE(k3xc)
1304 28 : ABI_FREE(indsym)
1305 28 : ABI_FREE(indsy1)
1306 28 : ABI_FREE(nattyp)
1307 28 : ABI_FREE(npwarr)
1308 28 : ABI_FREE(symrec)
1309 28 : ABI_FREE(symrc1)
1310 28 : ABI_FREE(symaf1)
1311 28 : ABI_FREE(symrl1)
1312 28 : ABI_FREE(tnons1)
1313 28 : ABI_FREE(irrzon1)
1314 28 : ABI_FREE(phnons1)
1315 28 : ABI_FREE(symq)
1316 28 : ABI_FREE(ph1d)
1317 28 : ABI_FREE(ph1df)
1318 28 : ABI_FREE(vtrial)
1319 28 : ABI_FREE(vxc)
1320 28 : call pawfgr_destroy(pawfgr)
1321 28 : if (psps%usepaw==1) then
1322 4 : call pawang_free(pawang1)
1323 4 : call pawrhoij_free(pawrhoij)
1324 4 : call paw_an_free(paw_an)
1325 4 : call paw_ij_free(paw_ij)
1326 4 : call pawfgrtab_free(pawfgrtab)
1327 : end if
1328 36 : ABI_FREE(pawrhoij)
1329 36 : ABI_FREE(paw_an)
1330 36 : ABI_FREE(paw_ij)
1331 36 : ABI_FREE(pawfgrtab)
1332 :
1333 : ! Clean the header
1334 28 : call hdr%free()
1335 :
1336 : !As the etotal energy has no meaning here, we set it to zero
1337 : !(to avoid meaningless side-effects when comparing ouputs...)
1338 28 : etotal = zero
1339 :
1340 140 : call timab(501,2,tsec)
1341 :
1342 : DBG_EXIT("COLL")
1343 :
1344 : contains
1345 : !!***
1346 :
1347 : !!****f* nonlinear/print_chi2
1348 : !! NAME
1349 : !! print_chi2
1350 : !!
1351 : !! FUNCTION
1352 : !! Print a third derivative tensor. Used only in nonlinear
1353 : !!
1354 : !! INPUTS
1355 : !! d3cart0 = the tensor to print
1356 : !! msg = a short message printed before the tensor
1357 : !! theunit = unit where the tensor is written
1358 : !!
1359 : !! SOURCE
1360 :
1361 5 : subroutine print_chi2(d3cart0,msg,theunit)
1362 :
1363 : integer,intent(in) :: theunit
1364 : character(len=30) :: msg
1365 : real(dp) :: elem1,elem2
1366 : real(dp),intent(in) :: d3cart0(2,3,mpert,3,mpert,3,mpert)
1367 : ! *************************************************************************
1368 :
1369 5 : write(theunit,'(2a)') ch10,msg
1370 20 : do i1dir = 1, 3
1371 65 : do i2dir = 1, 3
1372 195 : do i3dir = 1, 3
1373 135 : elem1 = d3cart0(1,i1dir,natom+2,i2dir,natom+2,i3dir,natom+2)
1374 135 : elem2 = d3cart0(2,i1dir,natom+2,i2dir,natom+2,i3dir,natom+2)
1375 180 : write(theunit,'(3(5x,i2),5x,f16.9,2x,f16.9)') i1dir,i2dir,i3dir,elem1,elem2
1376 : end do
1377 : end do
1378 : end do
1379 :
1380 5 : end subroutine print_chi2
1381 : !!***
1382 :
1383 : !!****f* nonlinear/print_dchidtau
1384 : !! NAME
1385 : !! print_dchidtau
1386 : !!
1387 : !! FUNCTION
1388 : !! Print a third derivative tensor. Used only in nonlinear
1389 : !!
1390 : !! INPUTS
1391 : !! d3cart0 = the tensor to print
1392 : !! msg = a short message printed before the tensor
1393 : !! theunit = unit where the tensor is written
1394 : !!
1395 : !! OUTPUT
1396 : !!
1397 : !! SIDE EFFECTS
1398 : !!
1399 : !! SOURCE
1400 :
1401 10 : subroutine print_dchidtau(d3cart0,msg,theunit)
1402 :
1403 : integer,intent(in) :: theunit
1404 : character(len=30) :: msg
1405 : real(dp),intent(in) :: d3cart0(2,3,mpert,3,mpert,3,mpert)
1406 : ! *************************************************************************
1407 :
1408 10 : write(theunit,'(a)') msg
1409 30 : do i1pert = 1,natom
1410 90 : do i1dir = 1,3
1411 60 : write(theunit,'(1x,i4,9x,i2,3(3x,f16.9),3(3x,f16.9))')i1pert,i1dir,&
1412 120 : d3cart0(1,i1dir,i1pert,1,natom+2,:,natom+2),d3cart0(2,i1dir,i1pert,1,natom+2,:,natom+2)
1413 : write(theunit,'(16x,3(3x,f16.9),3(3x,f16.9))')&
1414 60 : d3cart0(1,i1dir,i1pert,2,natom+2,:,natom+2),d3cart0(2,i1dir,i1pert,2,natom+2,:,natom+2)
1415 : write(theunit,'(16x,3(3x,f16.9),3(3x,f16.9))')&
1416 80 : d3cart0(1,i1dir,i1pert,3,natom+2,:,natom+2),d3cart0(2,i1dir,i1pert,3,natom+2,:,natom+2)
1417 : end do
1418 : end do
1419 :
1420 10 : end subroutine print_dchidtau
1421 : !!***
1422 :
1423 : end subroutine nonlinear
1424 : !!***
1425 :
1426 : !!****f* ABINIT/initmv
1427 : !! NAME
1428 : !! initmv
1429 : !!
1430 : !! FUNCTION
1431 : !! Initialize finite difference calculation of the ddk im dfptnl_mv.f
1432 : !!
1433 : !! INPUTS
1434 : !! dtset <type(dataset_type)> = all input variables in this dataset
1435 : !! gmet(3,3) = reciprocal space metric tensor in bohr**-2
1436 : !! kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere
1437 : !! kneigh(30,nkpt2) = index of the neighbours of each k-point
1438 : !! kg_neigh(30,nkpt2,3) = necessary to construct the vector joining a k-point
1439 : !! to its nearest neighbour in case of a single k-point,
1440 : !! a line of k-points or a plane of k-points.
1441 : !! See getshell.F90 for details
1442 : !! kptindex(2,nkpt3)= index of the k-points in the reduced BZ
1443 : !! related to a k-point in the full BZ
1444 : !! kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ
1445 : !! mband = maximum number of bands
1446 : !! mkmem = number of k points which can fit in memory
1447 : !! mpi_enreg = information about MPI parallelization
1448 : !! mpw = maximum number of plane waves
1449 : !! nband(nkpt*nsppol)=number of bands at each k point, for each polarization
1450 : !! nkpt2 = number of k-points in the reduced BZ
1451 : !! nkpt3 = number of k-points in the full BZ
1452 : !! nneigh = total number of neighbours required to evaluate the finite
1453 : !! difference formula
1454 : !! npwarr(nkpt2)=number of planewaves at each k point
1455 : !! nsppol = number of spin polarizations
1456 : !! occ(mband*nkpt*nsppol) = occupation number for each band for each k
1457 : !!
1458 : !! OUTPUT
1459 : !! cgindex(nkpt2,nsppol) = for each k-point, cgindex tores the location of the WF in the cg array
1460 : !! me = index of the current processor
1461 : !! ineigh = index of a neighbour
1462 : !! ikpt_loc = index of the iteration on ikpt on the current processor
1463 : !! ikpt_rbz = index of a k-point in the reduced BZ
1464 : !! pwind(mpw,nneigh,mkmem) = array used to compute the overlap matrix smat between k-points
1465 : !! (see initberry.f for more explanations)
1466 : !!
1467 : !! SOURCE
1468 :
1469 21 : subroutine initmv(cgindex,dtset,gmet,kg,kneigh,kg_neigh,kptindex,&
1470 21 : & kpt3,mband,mkmem,mpi_enreg,mpw,nband,nkpt2,&
1471 21 : & nkpt3,nneigh,npwarr,nsppol,occ,pwind)
1472 :
1473 : !Arguments ------------------------------------
1474 : !scalars
1475 : integer,intent(in) :: mband,mkmem,mpw,nkpt2,nkpt3,nneigh,nsppol
1476 : type(MPI_type),intent(inout) :: mpi_enreg
1477 : type(dataset_type),intent(in) :: dtset
1478 : !arrays
1479 : integer,intent(in) :: kg(3,mpw*mkmem),kneigh(30,nkpt2),kg_neigh(30,nkpt2,3)
1480 : integer,intent(in) :: nband(nkpt2*nsppol),npwarr(nkpt2),kptindex(2,nkpt3)
1481 : integer,intent(out) :: cgindex(nkpt2,nsppol),pwind(mpw,nneigh,mkmem)
1482 : real(dp),intent(in) :: gmet(3,3),kpt3(3,nkpt3),occ(mband*nkpt2*nsppol)
1483 :
1484 : !Local variables-------------------------------
1485 : !scalars
1486 : integer :: flag,iband,icg,ierr,ikg,ikg1,ikpt,ikpt2,ikpt_loc,ikpt_rbz
1487 : integer :: index,ineigh,ipw,isppol,jpw,nband_k,mband_occ,mband_occ_k,npw_k
1488 : integer :: npw_k1,orig,spaceComm
1489 : real(dp) :: ecut_eff,sdeg
1490 : character(len=500) :: msg
1491 : !arrays
1492 : integer :: dg(3)
1493 21 : integer,allocatable :: kg1(:,:),kg1_k(:,:),npwar1(:),npwtot(:)
1494 : real(dp) :: dk(3),dk_(3)
1495 21 : real(dp),allocatable :: kpt1(:,:)
1496 :
1497 : !************************************************************************
1498 :
1499 : if (xmpi_paral== 1) then
1500 21 : spaceComm=mpi_enreg%comm_cell
1501 2567 : mpi_enreg%kpt_loc2ibz_sp(:,:,:) = 0
1502 56 : mpi_enreg%mkmem(:) = 0
1503 : end if
1504 :
1505 21 : ecut_eff = dtset%ecut*(dtset%dilatmx)**2
1506 63 : ABI_MALLOC(kg1_k,(3,mpw))
1507 63 : ABI_MALLOC(kg1,(3,mkmem*mpw))
1508 63 : ABI_MALLOC(kpt1,(3,nkpt2))
1509 63 : ABI_MALLOC(npwar1,(nkpt2))
1510 42 : ABI_MALLOC(npwtot,(nkpt2))
1511 29997 : kg1_k(:,:) = 0
1512 328593 : pwind(:,:,:) = 0
1513 717 : cgindex(:,:) = 0
1514 :
1515 : !Compute the number of occupied bands.
1516 : !Check that it is the same for every k-point and that
1517 : !nband(ikpt) is equal to this value
1518 :
1519 21 : if (nsppol == 1) then
1520 : sdeg = two
1521 3 : else if (nsppol == 2) then
1522 3 : sdeg = one
1523 : end if
1524 :
1525 : !DEBUG
1526 : !write(std_out,*)' list of nband '
1527 : !do isppol = 1, nsppol
1528 : !do ikpt = 1, nkpt2
1529 : !nband_k = nband(ikpt + (isppol - 1)*nkpt2)
1530 : !write(std_out,*)' isppol, ikpt, nband_k=',isppol, ikpt, nband_k
1531 : !end do
1532 : !end do
1533 : !ENDDEBUG
1534 :
1535 21 : index = 0
1536 45 : do isppol = 1, nsppol
1537 717 : do ikpt = 1, nkpt2
1538 :
1539 672 : mband_occ_k = 0
1540 672 : nband_k = nband(ikpt + (isppol - 1)*nkpt2)
1541 :
1542 3280 : do iband = 1, nband_k
1543 2608 : index = index + 1
1544 3280 : if (abs(occ(index) - sdeg) < tol8) mband_occ_k = mband_occ_k + 1
1545 : end do
1546 :
1547 672 : if (nband_k /= mband_occ_k) then
1548 : write(msg,'(a,a,a)')&
1549 0 : & ' In a non-linear response calculation, nband must be equal ',ch10,&
1550 0 : & ' to the number of valence bands.'
1551 0 : ABI_ERROR(msg)
1552 : end if
1553 :
1554 : ! Note that the number of bands can be different for spin up and spin down
1555 696 : if (ikpt > 1) then
1556 648 : if (mband_occ /= mband_occ_k) then
1557 0 : ABI_ERROR('The number of valence bands is not the same for every k-point')
1558 : end if
1559 : else
1560 : mband_occ = mband_occ_k
1561 : end if
1562 :
1563 : end do ! close loop over ikpt
1564 : end do ! close loop over isppol
1565 :
1566 : !Find the location of each wavefunction
1567 :
1568 21 : icg = 0
1569 45 : do isppol = 1, nsppol
1570 717 : do ikpt = 1, nkpt2
1571 : ! fab: inserted the shift due to the spin...
1572 672 : nband_k = dtset%nband(ikpt+(isppol - 1)*nkpt2)
1573 672 : npw_k = npwarr(ikpt)
1574 :
1575 672 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me)) cycle
1576 :
1577 608 : cgindex(ikpt,isppol) = icg
1578 696 : icg = icg + dtset%nspinor*npw_k*nband_k
1579 :
1580 : end do
1581 : end do
1582 :
1583 :
1584 : !Build pwind
1585 :
1586 147 : do ineigh = 1, nneigh
1587 :
1588 4074 : do ikpt = 1, nkpt2
1589 3948 : ikpt2 = kneigh(ineigh,ikpt)
1590 3948 : ikpt_rbz = kptindex(1,ikpt2) ! index of the k-point in the reduced BZ
1591 15918 : kpt1(:,ikpt) = dtset%kptns(:,ikpt_rbz)
1592 : end do
1593 :
1594 : ! Set up the basis sphere of plane waves at kpt1
1595 1297782 : kg1(:,:) = 0
1596 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg1,&
1597 : & kpt1,mkmem,dtset%nband,nkpt2,'PERS',mpi_enreg,mpw,&
1598 126 : & npwar1,npwtot,dtset%nsppol)
1599 :
1600 126 : ikg = 0 ; ikg1 = 0 ; ikpt_loc = 0
1601 :
1602 126 : if(dtset%nsppol/=1)then
1603 18 : if(mpi_enreg%nproc/=1)then
1604 0 : ABI_ERROR('At present, non-linear response calculations for spin-polarized system cannot be done in parallel.')
1605 : else
1606 18 : isppol=1
1607 : end if
1608 : else
1609 108 : isppol=1
1610 : end if
1611 :
1612 4095 : do ikpt = 1, nkpt2
1613 :
1614 3948 : nband_k = dtset%nband(ikpt+(isppol - 1)*nkpt2)
1615 3948 : ikpt2 = kneigh(ineigh,ikpt)
1616 3948 : ikpt_rbz = kptindex(1,ikpt2) ! index of the k-point in the reduced BZ
1617 :
1618 3948 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,-1,mpi_enreg%me)) cycle
1619 :
1620 3564 : ikpt_loc = ikpt_loc + 1
1621 :
1622 3564 : mpi_enreg%kpt_loc2ibz_sp(mpi_enreg%me, ikpt_loc, 1) = ikpt
1623 :
1624 3564 : flag = 0
1625 3564 : npw_k = npwarr(ikpt)
1626 3564 : npw_k1 = npwarr(ikpt_rbz)
1627 14256 : dk_(:) = kpt3(:,ikpt2) - dtset%kptns(:,ikpt)
1628 14256 : dk(:) = dk_(:) - nint(dk_(:)) + real(kg_neigh(ineigh,ikpt,:),dp)
1629 14256 : dg(:) = nint(dk(:) - dk_(:))
1630 :
1631 :
1632 3564 : if (kptindex(2,ikpt2) == 0) then
1633 1071030 : kg1_k(:,1:npw_k1) = kg1(:,ikg1+1:ikg1+npw_k1)
1634 2854 : if (dg(1)==0.and.dg(2)==0.and.dg(3)==0) flag = 1
1635 : else
1636 184494 : kg1_k(:,1:npw_k1) = -1*kg1(:,ikg1+1:ikg1+npw_k1)
1637 : end if
1638 :
1639 3564 : orig = 1
1640 316554 : do ipw = 1, npw_k
1641 61901662 : do jpw = orig, npw_k1
1642 :
1643 : if ((kg(1,ikg + ipw) == kg1_k(1,jpw) - dg(1)).and. &
1644 61875792 : & (kg(2,ikg + ipw) == kg1_k(2,jpw) - dg(2)).and. &
1645 22306 : & (kg(3,ikg + ipw) == kg1_k(3,jpw) - dg(3))) then
1646 :
1647 290684 : pwind(ipw,ineigh,ikpt_loc) = jpw
1648 290684 : if (flag == 1) orig = jpw + 1
1649 : exit
1650 :
1651 : end if
1652 :
1653 : end do
1654 : end do
1655 :
1656 3564 : ikg = ikg + npw_k
1657 4074 : ikg1 = ikg1 + npw_k1
1658 :
1659 : end do ! close loop over k-points
1660 : end do ! close loop over ineigh
1661 :
1662 21 : mpi_enreg%mkmem(mpi_enreg%me) = mkmem
1663 :
1664 21 : call xmpi_sum(mpi_enreg%kpt_loc2ibz_sp,spaceComm,ierr)
1665 21 : call xmpi_sum(mpi_enreg%mkmem,spaceComm,ierr)
1666 :
1667 21 : ABI_FREE(kg1)
1668 21 : ABI_FREE(kg1_k)
1669 21 : ABI_FREE(kpt1)
1670 21 : ABI_FREE(npwar1)
1671 21 : ABI_FREE(npwtot)
1672 :
1673 21 : end subroutine initmv
1674 : !!***
1675 :
1676 : !----------------------------------------------------------------------
1677 :
1678 : !!****f* m_nonlinear/dfptnl_doutput
1679 : !! NAME
1680 : !! dfptnl_doutput
1681 : !!
1682 : !! FUNCTION
1683 : !! Write the matrix of third-order derivatives to the output file
1684 : !!
1685 : !! INPUTS
1686 : !! blkflg(3,mpert,3,mpert,3,mpert)= ( 1 if the element of the 3dte
1687 : !! has been calculated ; 0 otherwise )
1688 : !! d3(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
1689 : !! mpert =maximum number of ipert
1690 : !! natom=Number of atoms
1691 : !! ntypat=Number of type of atoms
1692 : !! unddb = unit number for DDB output
1693 : !!
1694 : !! NOTES
1695 : !! d3 holds the third-order derivatives before computing
1696 : !! the permutations of the perturbations.
1697 : !!
1698 : !! SOURCE
1699 :
1700 24 : subroutine dfptnl_doutput(blkflg,d3,mpert)
1701 :
1702 : !Arguments -------------------------------
1703 : !scalars
1704 : integer,intent(in) :: mpert
1705 : !arrays
1706 : integer,intent(in) :: blkflg(3,mpert,3,mpert,3,mpert)
1707 : real(dp),intent(in) :: d3(2,3,mpert,3,mpert,3,mpert)
1708 :
1709 : !Local variables -------------------------
1710 : !scalars
1711 : integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert
1712 : character(len=500) :: msg
1713 :
1714 : !*************************************************************************
1715 :
1716 : ! Write blok of third-order derivatives to ouput file
1717 :
1718 24 : write(msg,'(a,a,a,a,a)')ch10,&
1719 24 : ' Matrix of third-order derivatives (reduced coordinates)',ch10,&
1720 48 : ' before computing the permutations of the perturbations',ch10
1721 24 : call wrtout(ab_out,msg)
1722 :
1723 24 : write(ab_out,*)' j1 j2 j3 matrix element'
1724 24 : write(ab_out,*)' dir pert dir pert dir pert real part imaginary part'
1725 :
1726 216 : do i1pert=1,mpert
1727 792 : do i1dir=1,3
1728 5376 : do i2pert=1,mpert
1729 19008 : do i2dir=1,3
1730 129024 : do i3pert=1,mpert
1731 456192 : do i3dir=1,3
1732 :
1733 442368 : if (blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=0) then
1734 :
1735 : write(ab_out,'(3(i4,i5),2f22.10)')&
1736 3080 : i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,&
1737 6160 : d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
1738 : end if
1739 :
1740 : end do
1741 : end do
1742 : end do
1743 : end do
1744 : end do
1745 : end do
1746 :
1747 24 : end subroutine dfptnl_doutput
1748 : !!***
1749 :
1750 : end module m_nonlinear
1751 : !!***
|
