Line data Source code
1 : !!****m* ABINIT/m_pead_nl_loop
2 : !! NAME
3 : !! m_pead_nl_loop
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2002-2026 ABINIT group (MVeithen,MB)
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_pead_nl_loop
23 :
24 : use defs_basis
25 : use defs_wvltypes
26 : use m_wffile
27 : use m_abicore
28 : use m_xmpi
29 : use m_hdr
30 : use m_dtset
31 : use m_dtfil
32 : USE_MPI
33 :
34 : use defs_datatypes, only : pseudopotential_type
35 : use defs_abitypes, only : MPI_type
36 : use m_time, only : timab
37 : use m_kg, only : getph, mkkpg
38 : use m_cgtools, only : dotprod_vn, dotprod_g
39 : use m_fft, only : fourdp, fftpac, fourwf
40 : use m_ioarr, only : read_rhor
41 : use m_pawtab, only : pawtab_type
42 : use m_pawrhoij, only : pawrhoij_type
43 : use m_pawcprj, only : pawcprj_type
44 : use m_inwffil, only : inwffil
45 : use m_spacepar, only : hartre
46 : use m_initylmg, only : initylmg
47 : use m_dfpt_mkvxc, only : dfpt_mkvxc
48 : use m_mkcore, only : dfpt_mkcore
49 : use m_mklocl, only : dfpt_vlocal
50 : use m_hamiltonian,only : gs_hamiltonian_type
51 : use m_mkffnl, only : mkffnl
52 : use m_mpinfo, only : proc_distrb_cycle
53 : use m_nonlop, only : nonlop
54 : use m_dfptnl_pert, only : dfptnl_exc3
55 :
56 : implicit none
57 :
58 : private
59 :
60 : #if defined HAVE_MPI1
61 : include 'mpif.h'
62 : #endif
63 : !!***
64 :
65 : public :: pead_nl_loop
66 : !!***
67 :
68 : contains
69 : !!***
70 :
71 : !!****f* ABINIT/pead_nl_loop
72 : !! NAME
73 : !! pead_nl_loop
74 : !!
75 : !! FUNCTION
76 : !! Loop over the perturbations j1, j2 and j3
77 : !!
78 : !! INPUTS
79 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions
80 : !! cgindex(nkpt,nsppol) = for each k-point, cgindex tores the location
81 : !! of the WF in the cg array
82 : !! dtfil <type(datafiles_type)>=variables related to files
83 : !! dtset <type(dataset_type)>=all input variables for this dataset
84 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
85 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
86 : !! gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double
87 : !! that of the basis sphere--appropriate for charge density rho(G),
88 : !! Hartree potential, and pseudopotentials
89 : !! kg(3,mpw*mkmem)=reduced planewave coordinates
90 : !! kneigh(30,nkpt) = index of the neighbours of each k-point
91 : !! kg_neigh(30,nkpt,3) = necessary to construct the vector joining a k-point
92 : !! to its nearest neighbour in case of a single k-point,
93 : !! a line of k-points or a plane of k-points.
94 : !! kptindex(2,nkpt3)= index of the k-points in the reduced BZ
95 : !! related to a k-point in the full BZ
96 : !! kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ
97 : !! kxc(nfft,nkxc)=exchange-correlation kernel
98 : !! k3xc(nfft,nk3xc)=third-order exchange-correlation kernel
99 : !! mband = maximum number of bands
100 : !! mgfft = maximum single fft dimension
101 : !! mkmem = Number of k points treated by this node.
102 : !! mkmem_max = maximal number of k-points on each processor (MPI //)
103 : !! mk1mem = Number of k points for first-order WF treated by this node.
104 : !! mpert =maximum number of ipert
105 : !! mpi_enreg=MPI-parallelisation information
106 : !! mpw = maximum number of planewaves in basis sphere (large number)
107 : !! mvwtk(30,nkpt) = weights to compute the finite difference ddk
108 : !! natom = number of atoms in unit cell
109 : !! nfft = (effective) number of FFT grid points (for this processor)
110 : !! nkpt = number of k points
111 : !! nkpt3 = number of k-points in the full BZ
112 : !! nkxc=second dimension of the array kxc, see rhohxc.f for a description
113 : !! nneigh = total number of neighbours required to evaluate the finite
114 : !! difference formula
115 : !! nspinor = number of spinorial components of the wavefunctions
116 : !! nsppol = number of channels for spin-polarization (1 or 2)
117 : !! npwarr(nkpt) = array holding npw for each k point
118 : !! occ(mband*nkpt*nsppol) = occupation number for each band and k
119 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
120 : !! pwind(mpw,nneigh,mkmem) = array used to compute the overlap matrix smat
121 : !! between k-points
122 : !! rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
123 : !! that have to be computed
124 : !! 1 -> element has to be computed explicitely
125 : !! -1 -> use symmetry operations to obtain the corresponding element
126 : !! rprimd(3,3)=dimensional primitive translations (bohr)
127 : !! ucvol = unit cell volume (bohr^3)
128 : !! xred(3,natom) = reduced atomic coordinates
129 : !!
130 : !! OUTPUT
131 : !! blkflg(3,mpert,3,mpert) = flags for each element of the 3DTE
132 : !! (=1 if computed)
133 : !! d3lo(2,3,mpert,3,mpert,3,mpert) = matrix of the 3DTEs
134 : !!
135 : !! SIDE EFFECTS
136 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
137 : !!
138 : !! SOURCE
139 :
140 21 : subroutine pead_nl_loop(blkflg,cg,cgindex,dtfil,dtset,d3lo,&
141 : & gmet,gprimd,gsqcut, &
142 21 : & hdr,kg,kneigh,kg_neigh,kptindex,kpt3,kxc,k3xc,mband,mgfft,mkmem,mkmem_max,mk1mem,&
143 21 : & mpert,mpi_enreg,mpw,mvwtk,natom,nfft,nkpt,nkpt3,nkxc,nk3xc,nneigh,nspinor,nsppol,&
144 21 : & npwarr,occ,psps,pwind,&
145 21 : & rfpert,rprimd,ucvol,xred)
146 :
147 : !Arguments ------------------------------------
148 : !scalars
149 : integer,intent(in) :: mband,mgfft,mk1mem,mkmem,mkmem_max,mpert,mpw,natom,nfft
150 : integer,intent(in) :: nk3xc,nkpt,nkpt3,nkxc,nneigh,nspinor,nsppol
151 : real(dp),intent(in) :: gsqcut,ucvol
152 : type(MPI_type),intent(inout) :: mpi_enreg
153 : type(datafiles_type),intent(in) :: dtfil
154 : type(dataset_type),intent(inout) :: dtset
155 : type(hdr_type),intent(inout) :: hdr
156 : type(pseudopotential_type),intent(in) :: psps
157 : !arrays
158 : integer,intent(in) :: cgindex(nkpt,nsppol),kg(3,mk1mem*mpw),kneigh(30,nkpt)
159 : integer,intent(in) :: kg_neigh(30,nkpt,3)
160 : integer,intent(in) :: kptindex(2,nkpt3),npwarr(nkpt),pwind(mpw,nneigh,mkmem)
161 : integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert)
162 : integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) !vz_i
163 : real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gmet(3,3)
164 : real(dp),intent(in) :: gprimd(3,3),k3xc(nfft,nk3xc),kpt3(3,nkpt3)
165 : real(dp),intent(in) :: kxc(nfft,nkxc),mvwtk(30,nkpt),rprimd(3,3)
166 : real(dp),intent(in) :: xred(3,natom)
167 : real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
168 : real(dp),intent(inout) :: d3lo(2,3,mpert,3,mpert,3,mpert) !vz_i
169 :
170 : !Local variables-------------------------------
171 : !scalars
172 : integer,parameter :: level=51
173 : integer :: ask_accurate,counter,cplex,formeig,i1dir
174 : integer :: i1pert,i2dir,i2pert,i3dir,i3pert,iatom,ierr,index,ir
175 : integer :: ireadwf,itypat,mcg,mpsang,n1,n2,n3,n3xccc,nfftot,nspden,option,optorth
176 : integer :: pert1case,pert2case,pert3case,rdwrpaw,timrev,comm_cell
177 : logical :: nmxc
178 : real(dp) :: ecut_eff,exc3(2)
179 : character(len=500) :: message
180 : character(len=fnlen) :: fiden1i,fiwf1i,fiwf3i
181 : type(wffile_type) :: wff1,wff2,wfft1,wfft2
182 21 : type(wvl_data) :: wvl
183 21 : type(hdr_type) :: hdr_den
184 : !arrays
185 21 : integer,allocatable :: atindx(:),atindx1(:),nattyp(:)
186 : real(dp) :: d3_berry(2,3),rho_dum(1),tsec(2),ylmgr_dum(1)
187 21 : real(dp),allocatable :: cg1(:,:),cg3(:,:),eigen1(:),ph1d(:,:),rho1r1(:,:)
188 21 : real(dp),allocatable :: rho2g1(:,:),rho2r1(:,:),rho3r1(:,:),vhartr1(:)
189 21 : real(dp),allocatable :: vpsp1(:),vtrial1(:,:),vxc1(:,:),work(:)
190 21 : real(dp),allocatable :: xccc3d1(:),xccc3d2(:),xccc3d3(:),ylm(:,:,:)
191 21 : type(pawrhoij_type),allocatable :: rhoij_dum(:)
192 :
193 : ! ***********************************************************************
194 :
195 21 : call timab(502,1,tsec)
196 :
197 21 : comm_cell = mpi_enreg%comm_cell
198 :
199 21 : timrev = 1
200 21 : cplex = 2 - timrev
201 21 : nspden = dtset%nspden
202 21 : ecut_eff = (dtset%ecut)*(dtset%dilatmx)**2
203 21 : mpsang = psps%mpsang
204 21 : optorth=1;if (psps%usepaw==1) optorth=0
205 :
206 63 : ABI_MALLOC(cg1,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
207 42 : ABI_MALLOC(cg3,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
208 63 : ABI_MALLOC(eigen1,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
209 84 : ABI_MALLOC(rho1r1,(cplex*nfft,dtset%nspden))
210 63 : ABI_MALLOC(rho2r1,(cplex*nfft,dtset%nspden))
211 63 : ABI_MALLOC(rho2g1,(2,nfft))
212 63 : ABI_MALLOC(rho3r1,(cplex*nfft,dtset%nspden))
213 84 : ABI_MALLOC(ylm,(2,dtset%mpw*dtset%mkmem,mpsang*mpsang*psps%useylm))
214 :
215 21 : ask_accurate=1 ; formeig = 1 ; ireadwf = 1
216 21 : n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
217 21 : nfftot=n1*n2*n3
218 :
219 : !Generate an index table of atoms, in order for them to be used
220 : !type after type.
221 63 : ABI_MALLOC(atindx,(natom))
222 42 : ABI_MALLOC(atindx1,(natom))
223 63 : ABI_MALLOC(nattyp,(psps%ntypat))
224 56 : index=1
225 56 : do itypat=1,psps%ntypat
226 35 : nattyp(itypat)=0
227 126 : do iatom=1,natom
228 105 : if(dtset%typat(iatom)==itypat)then
229 42 : atindx(iatom)=index
230 42 : atindx1(index)=iatom
231 42 : index=index+1
232 42 : nattyp(itypat)=nattyp(itypat)+1
233 : end if
234 : end do
235 : end do
236 :
237 : !Generate the 1-dimensional phases
238 63 : ABI_MALLOC(ph1d,(2,3*(2*mgfft+1)*natom))
239 21 : call getph(atindx,natom,n1,n2,n3,ph1d,xred)
240 :
241 : !Set up the Ylm for each k point
242 21 : if (psps%useylm==1) then
243 : call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,&
244 : & dtset%mpw,dtset%nband,dtset%nkpt,&
245 0 : & npwarr,dtset%nsppol,0,rprimd,ylm,ylmgr_dum)
246 : end if
247 :
248 63 : ABI_MALLOC(vpsp1,(cplex*nfft))
249 42 : ABI_MALLOC(xccc3d1,(cplex*nfft))
250 42 : ABI_MALLOC(xccc3d2,(cplex*nfft))
251 42 : ABI_MALLOC(xccc3d3,(cplex*nfft))
252 42 : ABI_MALLOC(vhartr1,(cplex*nfft))
253 84 : ABI_MALLOC(vxc1,(cplex*nfft,dtset%nspden))
254 63 : ABI_MALLOC(vtrial1,(cplex*nfft,dtset%nspden))
255 :
256 : !Loop over the perturbations j1, j2, j3
257 :
258 21 : pert1case = 0 ; pert2case = 0 ; pert3case = 0
259 :
260 189 : do i1pert = 1, mpert
261 693 : do i1dir = 1, 3
262 :
263 403872 : if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then
264 :
265 70 : pert1case = i1dir + (i1pert-1)*3
266 70 : counter = pert1case
267 70 : call appdig(pert1case,dtfil%fnamewff1,fiwf1i)
268 :
269 70 : mcg=mpw*nspinor*mband*mkmem*nsppol
270 : call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
271 : & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
272 : & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
273 : & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
274 : & dtset%nsppol,dtset%nsym,&
275 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
276 70 : & dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl)
277 :
278 70 : if (ireadwf==1) then
279 70 : call WffClose (wff1,ierr)
280 : end if
281 :
282 690345 : rho1r1(:,:) = 0._dp
283 70 : if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
284 70 : rdwrpaw=0
285 70 : call appdig(pert1case,dtfil%fildens1in,fiden1i)
286 :
287 : call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho1r1, &
288 70 : hdr_den, rhoij_dum, comm_cell, check_hdr=hdr)
289 70 : call hdr_den%free()
290 : end if
291 :
292 586966 : xccc3d1(:) = 0._dp
293 70 : if ((psps%n1xccc/=0).and.(i1pert <= natom)) then
294 : call dfpt_mkcore(cplex,i1dir,i1pert,natom,psps%ntypat,n1,psps%n1xccc,&
295 : & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,&
296 22 : & psps%xcccrc,psps%xccc1d,xccc3d1,xred)
297 : end if ! psps%n1xccc/=0
298 :
299 630 : do i3pert = 1, mpert
300 2310 : do i3dir = 1, 3
301 :
302 56000 : if ((maxval(rfpert(i1dir,i1pert,:,:,i3dir,i3pert))==1)) then
303 :
304 145 : pert3case = i3dir + (i3pert-1)*3
305 145 : counter = 100*pert3case + pert1case
306 145 : call appdig(pert3case,dtfil%fnamewff1,fiwf3i)
307 :
308 145 : mcg=mpw*nspinor*mband*mkmem*nsppol
309 : call inwffil(ask_accurate,cg3,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
310 : & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
311 : & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
312 : & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
313 : & dtset%nsppol,dtset%nsym,&
314 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
315 : & dtfil%unkg1,wff2,wfft2,dtfil%unwff2,&
316 145 : & fiwf3i,wvl)
317 145 : if (ireadwf==1) then
318 145 : call WffClose (wff2,ierr)
319 : end if
320 :
321 1272062 : rho3r1(:,:) = 0._dp
322 145 : if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
323 145 : rdwrpaw=0
324 145 : call appdig(pert3case,dtfil%fildens1in,fiden1i)
325 :
326 : call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho3r1, &
327 145 : hdr_den, rhoij_dum, comm_cell, check_hdr=hdr)
328 145 : call hdr_den%free()
329 : end if
330 :
331 1076433 : xccc3d3(:) = 0._dp
332 145 : if ((psps%n1xccc/=0).and.(i3pert <= natom)) then
333 : call dfpt_mkcore(cplex,i3dir,i3pert,natom,psps%ntypat,n1,psps%n1xccc,&
334 : & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,&
335 0 : & psps%xcccrc,psps%xccc1d,xccc3d3,xred)
336 : end if ! psps%n1xccc/=0
337 :
338 1305 : do i2pert = 1, mpert
339 :
340 : ! In case of electric field perturbation, evaluate the ddk
341 : ! using the finite difference expression of
342 : ! Marzari and Vanderbilt PRB 56, 12847 (1997) [[cite:Marzari1997]].
343 :
344 1160 : d3_berry(:,:) = 0._dp
345 :
346 4640 : if ((i2pert==dtset%natom+2).and.&
347 : & (maxval(rfpert(i1dir,i1pert,:,i2pert,i3dir,i3pert)) == 1)) then
348 :
349 145 : call timab(511,1,tsec)
350 : call pead_nl_mv(cg,cgindex,cg1,cg3,dtset,dtfil,d3_berry,gmet,&
351 : & i1pert,i3pert,i1dir,i3dir,&
352 : & kneigh,kg_neigh,kptindex,kpt3,mband,mkmem,mkmem_max,mk1mem,&
353 145 : & mpi_enreg,mpw,mvwtk,natom,nkpt,nkpt3,nneigh,npwarr,nspinor,nsppol,pwind)
354 145 : call timab(511,2,tsec)
355 :
356 : end if
357 :
358 1160 : if (mpi_enreg%me == 0) then
359 :
360 3744 : if(sum(rfpert(i1dir,i1pert,:,i2pert,i3dir,i3pert))>0)then
361 120 : write(message,'(a,a,a,a,a,a)')ch10,ch10,&
362 240 : & ' Decomposition of the third-order energy for the set of perturbations',ch10
363 120 : call wrtout(std_out,message,'COLL')
364 120 : call wrtout(ab_out,message,'COLL')
365 120 : if (i1pert < natom + 1) then
366 : write(message,'(a,i3,a,i3)') &
367 56 : & ' j1 : displacement of atom ',i1pert,' along direction ', i1dir
368 : end if
369 120 : if (i1pert == dtset%natom + 2) then
370 64 : write(message,'(a,i4)')' j1 : homogeneous electric field along direction ',i1dir
371 : end if
372 120 : call wrtout(std_out,message,'COLL')
373 120 : call wrtout(ab_out,message,'COLL')
374 120 : if (i3pert < natom + 1) then
375 : write(message,'(a,i3,a,i3,a)') &
376 0 : & ' j3 : displacement of atom ',i3pert,' along direction ', i3dir,ch10
377 : end if
378 120 : if (i3pert == dtset%natom + 2) then
379 120 : write(message,'(a,i4,a)')' j3 : homogeneous electric field along direction ',i3dir,ch10
380 : end if
381 120 : call wrtout(std_out,message,'COLL')
382 120 : call wrtout(ab_out,message,'COLL')
383 : end if
384 :
385 : end if ! mpi_enreg%me == 0
386 :
387 4785 : do i2dir = 1, 3
388 :
389 4640 : if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
390 407 : pert2case = i2dir + (i2pert-1)*3
391 :
392 407 : blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
393 :
394 : ! Read the first-order densities from disk-files
395 10999383 : rho2r1(:,:) = 0._dp ; rho2g1(:,:) = 0._dp
396 :
397 407 : if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
398 407 : rdwrpaw=0
399 407 : call appdig(pert2case,dtfil%fildens1in,fiden1i)
400 :
401 : call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho2r1, &
402 407 : hdr_den, rhoij_dum, comm_cell, check_hdr=hdr)
403 407 : call hdr_den%free()
404 :
405 : ! Compute up+down rho1(G) by fft
406 814 : ABI_MALLOC(work,(cplex*nfft))
407 2566031 : work(:)=rho2r1(:,1)
408 407 : call fourdp(cplex,rho2g1,work,-1,mpi_enreg,nfft,1,dtset%ngfft,0)
409 407 : ABI_FREE(work)
410 :
411 : end if
412 :
413 : ! Compute first-order local potentials
414 : ! (hartree, xc and pseudopotential)
415 :
416 407 : n3xccc=0; if(psps%n1xccc/=0)n3xccc=nfft
417 5132062 : xccc3d2(:)=0._dp ; vpsp1(:)=0._dp
418 :
419 407 : if (i2pert <= natom) then
420 :
421 : call dfpt_vlocal(atindx,cplex,gmet,gsqcut,dtset%icutcoul,i2dir,i2pert,mpi_enreg,psps%mqgrid_vl,natom,&
422 : & nattyp,nfft,dtset%ngfft,nkpt,psps%ntypat,n1,n2,n3,ph1d,psps%qgrid_vl,&
423 126 : & dtset%qptn,dtset%rcut,rprimd,ucvol,dtset%vcutgeo,psps%vlspl,vpsp1,xred)
424 :
425 126 : if (psps%n1xccc/=0) then
426 : call dfpt_mkcore(cplex,i2dir,i2pert,natom,psps%ntypat,n1,psps%n1xccc,&
427 : & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,&
428 90 : & psps%xcccrc,psps%xccc1d,xccc3d2,xred)
429 : end if ! psps%n1xccc/=0
430 :
431 : end if ! i2pert <= natom
432 :
433 : call hartre(cplex,gsqcut,3,0,mpi_enreg,nfft,dtset%ngfft,dtset%nkpt,&
434 407 : &dtset%rcut,rho2g1,rprimd,dtset%vcutgeo,vhartr1)
435 407 : option=1 ; nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
436 : call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfft,dtset%ngfft,&
437 : & rho_dum,0,rho_dum,0,nkxc,nmxc,dtset%nspden,n3xccc,option,&
438 407 : & dtset%qptn,rho2r1,rprimd,0,vxc1,xccc3d2)
439 :
440 407 : if(dtset%nsppol==1)then
441 : if(cplex==1)then
442 1830365 : do ir=1,nfft
443 1830365 : vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1)
444 : end do
445 : else
446 : do ir=1,nfft
447 : vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1)
448 : vtrial1(2*ir ,1)=vpsp1(2*ir )+vhartr1(2*ir )+vxc1(2*ir ,1)
449 : end do
450 : end if
451 : else
452 : if(cplex==1)then
453 735666 : do ir=1,nfft
454 735552 : vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1)
455 735666 : vtrial1(ir,2)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,2)
456 : end do
457 : else
458 : ! fab: I think there was an error in the definition of vtrial1(2*ir-1,2); I have corrected it...
459 : do ir=1,nfft
460 : vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1)
461 : vtrial1(2*ir ,1)=vpsp1(2*ir )+vhartr1(2*ir )+vxc1(2*ir ,1)
462 : vtrial1(2*ir-1,2)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1 ,2)
463 : vtrial1(2*ir ,2)=vpsp1(2*ir )+vhartr1(2*ir )+vxc1(2*ir ,2)
464 : end do
465 : end if
466 : end if
467 :
468 : ! Compute the third-order xc energy
469 : call dfptnl_exc3(cplex,exc3,k3xc,mpi_enreg,nk3xc,nfft,nfftot,dtset%nspden,&
470 407 : & rho1r1,rho2r1,rho3r1,ucvol,xccc3d1,xccc3d2,xccc3d3)
471 :
472 : ! Perform DFPT part of the 3dte calculation
473 :
474 407 : call timab(512,1,tsec)
475 : call pead_nl_resp(cg,cg1,cg3,cplex,dtfil,dtset,d3lo,i1dir,i2dir,i3dir,i1pert,i2pert,i3pert,&
476 : & kg,mband,mgfft,mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,nfft,nkpt,nspden,&
477 407 : & nspinor,nsppol,npwarr,occ,ph1d,psps,rprimd,vtrial1,xred,ylm)
478 407 : call timab(512,2,tsec)
479 :
480 :
481 : ! Describe the perturbation and write out the result
482 407 : if (mpi_enreg%me == 0) then
483 347 : if (i2pert < natom + 1) then
484 : write(message,'(a,i3,a,i3)') &
485 110 : & ' j2 : displacement of atom ',i2pert,&
486 220 : & ' along direction ', i2dir
487 : end if
488 347 : if (i2pert == dtset%natom + 2) then
489 : write(message,'(a,i4)') &
490 237 : & ' j2 : homogeneous electric field along direction ',&
491 474 : & i2dir
492 : end if
493 347 : call wrtout(std_out,message,'COLL')
494 347 : call wrtout(ab_out,message,'COLL')
495 347 : write(ab_out,'(20x,a,13x,a)')'real part','imaginary part'
496 347 : write(ab_out,'(5x,a2,1x,f22.10,3x,f22.10)')'xc',exc3(1)*sixth,zero
497 347 : if (i2pert == natom + 2) then
498 237 : write(ab_out,'(5x,a3,f22.10,3x,f22.10)')'ddk',&
499 474 : & d3_berry(1,i2dir),d3_berry(2,i2dir)
500 : end if
501 347 : write(ab_out,'(5x,a3,f22.10,3x,f22.10)')'dft',&
502 347 : & d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),&
503 694 : & d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
504 347 : write(ab_out,*)
505 347 : write(std_out,'(18x,a,11x,a)')'real part','imaginary part'
506 347 : write(std_out,'(5x,a2,1x,f20.10,3x,f20.10)')'xc',exc3(1)*sixth,zero
507 347 : write(std_out,'(5x,a3,f22.10,3x,f22.10)')'ddk',&
508 694 : & d3_berry(1,i2dir),d3_berry(2,i2dir)
509 347 : write(std_out,'(5x,a3,f22.10,3x,f22.10)')'dft',&
510 347 : & d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),&
511 694 : & d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
512 347 : write(std_out,*)
513 : end if ! mpi_enreg%me == 0
514 :
515 : d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = &
516 407 : & d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + exc3(1)*sixth
517 : d3lo(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = &
518 1221 : & d3lo(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + d3_berry(:,i2dir)
519 :
520 : end if !rfpert
521 : end do !i2dir
522 : end do ! i2pert
523 :
524 : end if ! rfpert
525 : end do ! i3dir
526 : end do ! i3pert
527 :
528 : end if ! rfpert
529 : end do ! i1dir
530 : end do ! i1pert
531 :
532 :
533 21 : ABI_FREE(cg1)
534 21 : ABI_FREE(cg3)
535 21 : ABI_FREE(eigen1)
536 21 : ABI_FREE(rho1r1)
537 21 : ABI_FREE(rho2r1)
538 21 : ABI_FREE(rho2g1)
539 21 : ABI_FREE(rho3r1)
540 21 : ABI_FREE(atindx1)
541 21 : ABI_FREE(atindx)
542 21 : ABI_FREE(nattyp)
543 21 : ABI_FREE(ph1d)
544 21 : ABI_FREE(ylm)
545 21 : ABI_FREE(vtrial1)
546 21 : ABI_FREE(vxc1)
547 21 : ABI_FREE(vhartr1)
548 21 : ABI_FREE(vpsp1)
549 21 : ABI_FREE(xccc3d1)
550 21 : ABI_FREE(xccc3d2)
551 21 : ABI_FREE(xccc3d3)
552 :
553 21 : call timab(502,2,tsec)
554 :
555 21 : end subroutine pead_nl_loop
556 : !!***
557 :
558 : !!****f* ABINIT/pead_nl_resp
559 : !! NAME
560 : !! pead_nl_resp
561 : !!
562 : !! FUNCTION
563 : !! Compute the linear response part to the 3dte
564 : !!
565 : !! INPUTS
566 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave
567 : !! coefficients of wavefunctions
568 : !! cg1 = first-order wavefunction relative to the perturbations i1pert
569 : !! cg3 = first-order wavefunction relative to the perturbations i3pert
570 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
571 : !! if 2, COMPLEX
572 : !! dtfil <type(datafiles_type)>=variables related to files
573 : !! dtset <type(dataset_type)>=all input variables for this dataset
574 : !! i1dir,i2dir,i3dir=directions of the corresponding perturbations
575 : !! i1pert,i2pert,i3pert = type of perturbation that has to be computed
576 : !! kg(3,mpw*mkmem)=reduced planewave coordinates
577 : !! mband = maximum number of bands
578 : !! mgfft=maximum size of 1D FFTs
579 : !! mkmem = maximum number of k points which can fit in core memory
580 : !! mk1mem = maximum number of k points for first-order WF
581 : !! which can fit in core memory
582 : !! mpert =maximum number of ipert
583 : !! mpi_enreg=MPI-parallelisation information
584 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
585 : !! mpw = maximum number of planewaves in basis sphere (large number)
586 : !! natom = number of atoms in unit cell
587 : !! nfft = (effective) number of FFT grid points (for this processor)
588 : !! nkpt = number of k points
589 : !! nspden = number of spin-density components
590 : !! nspinor = number of spinorial components of the wavefunctions
591 : !! nsppol = number of channels for spin-polarization (1 or 2)
592 : !! npwarr(nkpt) = array holding npw for each k point
593 : !! occ(mband*nkpt*nsppol) = occupation number for each band and k
594 : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
595 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
596 : !! rprimd(3,3) = dimensional primitive translations (bohr)
597 : !! vtrial1(cplex*nfft,nspden)=firs-order local potential
598 : !! xred(3,natom) = reduced atomic coordinates
599 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= spherical harmonics for
600 : !! each G and k point
601 : !!
602 : !! OUTPUT
603 : !! d3lo(2,3,mpert,3,mpert,3,mpert) = matrix of the 3DTEs
604 : !!
605 : !! SOURCE
606 :
607 407 : subroutine pead_nl_resp(cg,cg1,cg3,cplex,dtfil,dtset,d3lo,&
608 : & i1dir,i2dir,i3dir,i1pert,i2pert,i3pert,&
609 407 : & kg,mband,mgfft,mkmem,mk1mem,&
610 : & mpert,mpi_enreg,mpsang,mpw,natom,nfft,nkpt,nspden,nspinor,nsppol,&
611 407 : & npwarr,occ,ph1d,psps,rprimd,vtrial1,xred,ylm)
612 :
613 : !Arguments ------------------------------------
614 : !scalars
615 : integer,intent(in) :: cplex,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband,mgfft
616 : integer,intent(in) :: mk1mem,mkmem,mpert,mpsang,mpw,natom,nfft,nkpt,nspden
617 : integer,intent(in) :: nspinor,nsppol
618 : type(MPI_type),intent(in) :: mpi_enreg
619 : type(datafiles_type),intent(in) :: dtfil
620 : type(dataset_type),intent(in) :: dtset
621 : type(pseudopotential_type),intent(in) :: psps
622 : !arrays
623 : integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt)
624 : real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
625 : real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol)
626 : real(dp),intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol)
627 : real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom),rprimd(3,3)
628 : real(dp),intent(in) :: xred(3,natom),ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
629 : real(dp),intent(inout) :: vtrial1(cplex*nfft,nspden)
630 : real(dp),intent(inout) :: d3lo(2,3,mpert,3,mpert,3,mpert)
631 :
632 : !Local variables-------------------------------
633 : !scalars
634 : integer,parameter :: level=52
635 : integer :: bantot,choice,counter,cpopt,dimffnl,iband,icg0,ider,ierr
636 : integer :: ii,ikg,ikpt,ilm,ipw,isppol,istwf_k,jband,jj
637 : integer :: me,n1,n2,n3,n4,n5,n6,nband_k,nkpg,nnlout,npw_k
638 : integer :: option,paw_opt,signs,spaceComm,tim_fourwf,tim_nonlop
639 : real(dp) :: dot1i,dot1r,dot2i,dot2r,doti,dotr,lagi,lagr,sumi,sumr,weight
640 407 : type(gs_hamiltonian_type) :: gs_hamk
641 : !arrays
642 407 : integer,allocatable :: kg_k(:,:)
643 : real(dp) :: buffer(2),enlout(3),kpq(3),kpt(3)
644 : real(dp) :: dum_svectout(1,1),dum(1),rmet(3,3),ylmgr_dum(1,1,1)
645 407 : real(dp),allocatable :: cwave0(:,:),cwavef3(:,:),ffnlk(:,:,:,:)
646 407 : real(dp),allocatable :: gh0(:,:),gh1(:,:),gvnl(:,:),kpg_k(:,:)
647 407 : real(dp),allocatable :: vlocal1(:,:,:),wfraug(:,:,:,:),ylm_k(:,:)
648 2849 : type(pawcprj_type) :: cprj_dum(1,1)
649 407 : type(pawtab_type) :: pawtab_dum(0)
650 :
651 : !***********************************************************************
652 :
653 : ABI_UNUSED(dtfil%ireadwf)
654 :
655 407 : me = mpi_enreg%me
656 407 : spaceComm=mpi_enreg%comm_cell
657 :
658 407 : bantot = 0
659 407 : icg0 = 0
660 407 : option = 2
661 407 : n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
662 407 : n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
663 :
664 2035 : ABI_MALLOC(vlocal1,(cplex*n4,n5,n6))
665 2035 : ABI_MALLOC(wfraug,(2,n4,n5,n6))
666 :
667 : !Initialize Hamiltonian (k-independent terms) - NCPP only
668 : call gs_hamk%init(psps,pawtab_dum,nspinor,nsppol,nspden,natom,&
669 : & dtset%typat,xred,nfft,mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,&
670 407 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
671 : !& paw_ij=paw_ij)
672 16280 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
673 :
674 407 : sumr = zero ; sumi = zero
675 :
676 : !Loop over spins
677 :
678 928 : do isppol = 1, nsppol
679 :
680 521 : call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,vtrial1,vlocal1,option)
681 :
682 : ! Loop over k-points
683 :
684 521 : ikg = 0
685 11130 : do ikpt = 1, nkpt
686 :
687 10202 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me))cycle
688 :
689 9242 : counter = 100*ikpt
690 :
691 9242 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
692 9242 : npw_k = npwarr(ikpt)
693 9242 : istwf_k = dtset%istwfk(ikpt)
694 :
695 36968 : kpt(:) = dtset%kptns(:,ikpt)
696 : kpq(:) = dtset%kptns(:,ikpt) ! In case of non zero q, kpt = kpt + q
697 :
698 27726 : ABI_MALLOC(cwave0,(2,npw_k*dtset%nspinor))
699 18484 : ABI_MALLOC(cwavef3,(2,npw_k*dtset%nspinor))
700 18484 : ABI_MALLOC(gh0,(2,npw_k*dtset%nspinor))
701 18484 : ABI_MALLOC(gvnl,(2,npw_k*dtset%nspinor))
702 18484 : ABI_MALLOC(gh1,(2,npw_k*dtset%nspinor))
703 :
704 27726 : ABI_MALLOC(kg_k,(3,npw_k))
705 36968 : ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
706 3846766 : kg_k(:,1:npw_k) = kg(:,1+ikg:npw_k+ikg)
707 9242 : if (psps%useylm==1) then
708 0 : do ilm=1,mpsang*mpsang
709 0 : ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
710 : end do
711 : end if
712 :
713 : ! Compute (k+G) and (k+q+G) vectors (only if useylm=1)
714 9242 : nkpg=0;if (i2pert<natom+1) nkpg=3*dtset%nloalg(3)
715 36968 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
716 9242 : if (nkpg>0) then
717 0 : call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k)
718 : end if
719 :
720 : ! Compute nonlocal form factors ffnl at (k+G), for all atoms
721 9242 : dimffnl=1
722 46210 : ABI_MALLOC(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
723 9242 : if (i2pert<natom+1) then
724 2569 : ider=0
725 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
726 : & ider,ider,psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,&
727 : & psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,&
728 2569 : & psps%usepaw,psps%useylm,ylm_k,ylmgr_dum)
729 : end if
730 :
731 : ! Load k-dependent part in the Hamiltonian datastructure
732 : call gs_hamk%load_k(kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,&
733 9242 : & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnlk,compute_gbound=.true.)
734 : ! Load k+q-dependent part in the Hamiltonian datastructure
735 : ! call load_kprime_hamiltonian... !! To be activated when q/=0
736 :
737 : ! Loop over bands
738 :
739 45970 : do iband = 1,nband_k
740 :
741 : cwave0(:,:)=cg(:,1+(iband - 1)*npw_k*dtset%nspinor+icg0:&
742 11068340 : & iband*npw_k*dtset%nspinor+icg0)
743 : cwavef3(:,:)=cg3(:,1+(iband-1)*npw_k*dtset%nspinor+icg0:&
744 11068340 : & iband*npw_k*dtset%nspinor+icg0)
745 :
746 : ! Compute vtrial1 | cwafef3 >
747 36728 : tim_fourwf = 0 ; weight = one
748 : call fourwf(cplex,vlocal1,cwavef3,gh1,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
749 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,n4,n5,n6,option,&
750 : & tim_fourwf,weight,weight,&
751 36728 : & gpu_option=dtset%gpu_option)
752 :
753 : ! In case i2pert = phonon-type perturbation
754 : ! add first-order change in the nonlocal potential
755 36728 : if (i2pert<natom+1) then
756 10196 : signs=2 ; choice=2 ; nnlout=3 ; tim_nonlop = 0 ; paw_opt=0 ; cpopt=-1
757 : call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,i2dir,dum,mpi_enreg,1,nnlout,paw_opt,&
758 10196 : & signs,dum_svectout,tim_nonlop,cwavef3,gvnl,iatom_only=i2pert)
759 3364556 : gh1(:,:) = gh1(:,:) + gvnl(:,:)
760 : end if
761 :
762 36728 : ii = (iband-1)*npw_k*dtset%nspinor + icg0
763 36728 : call dotprod_g(dotr,doti,istwf_k,npw_k,2,cg1(:,ii+1:ii+npw_k),gh1,mpi_enreg%me_g0,xmpi_comm_self)
764 :
765 : ! Compute vtrial1 | cwave0 >
766 : tim_fourwf = 0 ; weight = one
767 : call fourwf(cplex,vlocal1,cwave0,gh0,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
768 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,n4,n5,n6,option,&
769 36728 : & tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
770 :
771 : ! In case i2pert = phonon-type perturbation
772 : ! add first-order change in the nonlocal potential
773 36728 : if (i2pert<natom+1) then
774 10196 : signs=2 ; choice=2 ; nnlout=3 ; tim_nonlop = 0 ; paw_opt=0 ; cpopt=-1
775 : call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,i2dir,dum,mpi_enreg,1,nnlout,paw_opt,&
776 10196 : & signs,dum_svectout,tim_nonlop,cwave0,gvnl,iatom_only=i2pert)
777 3364556 : gh0(:,:) = gh0(:,:) + gvnl(:,:)
778 : end if
779 :
780 : ! Compute the dft contribution to the Lagrange multiplier
781 : ! cwavef3 and cwave0 have been transferred to gh1 and gh0
782 : ! these vectors will be used to store the wavefunctions of band iband
783 : ! cg1 and gh0 contain the wavefunctions of band jband
784 :
785 36728 : lagr = zero ; lagi = zero
786 183250 : do jband = 1, nband_k
787 :
788 146522 : ii = (jband - 1)*npw_k*dtset%nspinor + icg0
789 146522 : jj = (iband - 1)*npw_k*dtset%nspinor + icg0
790 :
791 : ! dot1r and dot1i contain < u_mk | v^(1) | u_nk >
792 : ! dot2r and dot2i contain < u_nk^(1) | u_mk^(1) >
793 : ! m -> jband and n -> iband
794 :
795 146522 : dot1r = zero ; dot1i = zero
796 146522 : dot2r = zero ; dot2i = zero
797 14594818 : do ipw = 1, npw_k
798 14448296 : ii = ii + 1 ; jj = jj + 1
799 14448296 : dot1r = dot1r + cg(1,ii)*gh0(1,ipw) + cg(2,ii)*gh0(2,ipw)
800 14448296 : dot1i = dot1i + cg(1,ii)*gh0(2,ipw) - cg(2,ii)*gh0(1,ipw)
801 : dot2r = dot2r + cg1(1,jj)*cg3(1,ii) + &
802 14448296 : & cg1(2,jj)*cg3(2,ii)
803 : dot2i = dot2i + cg1(1,jj)*cg3(2,ii) - &
804 14594818 : & cg1(2,jj)*cg3(1,ii)
805 : end do ! ipw
806 :
807 146522 : lagr = lagr + dot1r*dot2r - dot1i*dot2i
808 183250 : lagi = lagi + dot1r*dot2i + dot1i*dot2r
809 :
810 : end do ! jband
811 :
812 : sumr = sumr + &
813 36728 : & dtset%wtk(ikpt)*occ(bantot+iband)*(dotr-lagr)
814 : sumi = sumi + &
815 45970 : & dtset%wtk(ikpt)*occ(bantot+iband)*(doti-lagi)
816 :
817 : end do ! end loop over bands
818 :
819 9242 : bantot = bantot + nband_k
820 9242 : icg0 = icg0 + npw_k*dtset%nspinor*nband_k
821 9242 : ikg = ikg + npw_k
822 :
823 9242 : ABI_FREE(cwave0)
824 9242 : ABI_FREE(cwavef3)
825 9242 : ABI_FREE(gh0)
826 9242 : ABI_FREE(gh1)
827 9242 : ABI_FREE(gvnl)
828 9242 : ABI_FREE(kg_k)
829 9242 : ABI_FREE(ylm_k)
830 9242 : ABI_FREE(ffnlk)
831 10723 : ABI_FREE(kpg_k)
832 :
833 : end do ! end loop over k-points
834 :
835 : end do ! end loop over spins
836 :
837 : if (xmpi_paral == 1) then
838 407 : buffer(1) = sumr ; buffer(2) = sumi
839 407 : call xmpi_sum(buffer,spaceComm,ierr)
840 407 : sumr = buffer(1) ; sumi = buffer(2)
841 : end if
842 :
843 :
844 407 : d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr
845 : !d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi
846 :
847 : !In some cases, the imaginary part is /= 0 because of the
848 : !use of time reversal symmetry
849 407 : d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero
850 :
851 407 : call gs_hamk%free()
852 :
853 407 : ABI_FREE(vlocal1)
854 407 : ABI_FREE(wfraug)
855 :
856 407 : end subroutine pead_nl_resp
857 : !!***
858 :
859 : !!****f* ABINIT/pead_nl_mv
860 : !! NAME
861 : !! pead_nl_mv
862 : !!
863 : !! FUNCTION
864 : !! Compute the finite difference expression of the k-point derivative
865 : !! using the PEAD formulation of the third-order energy
866 : !! (see Nunes and Gonze PRB 63, 155107 (2001) [[cite:Nunes2001]] Eq. 102)
867 : !! and the finite difference formula of Marzari and Vanderbilt
868 : !! (see Marzari and Vanderbilt, PRB 56, 12847 (1997) [[cite:Marzari1997]], Appendix B)
869 : !!
870 : !! INPUTS
871 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions
872 : !! cgindex(nkpt2,nsppol) = for each k-point, cgindex stores the location of the WF in the cg array
873 : !! cg1 = first-order wavefunction relative to the perturbations i1pert
874 : !! cg3 = first-order wavefunction relative to the perturbations i3pert
875 : !! dtfil <type(datafiles_type)>=variables related to files
876 : !! dtset <type(dataset_type)>=all input variables for this dataset
877 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
878 : !! i1pert,i3pert = type of perturbation that has to be computed
879 : !! i1dir,i3dir=directions of the corresponding perturbations
880 : !! kneigh(30,nkpt2) = index of the neighbours of each k-point
881 : !! kg_neigh(30,nkpt2,3) = necessary to construct the vector joining a k-point
882 : !! to its nearest neighbour in case of a single k-point,
883 : !! a line of k-points or a plane of k-points.
884 : !! See getshell.F90 for details
885 : !! kptindex(2,nkpt3)= index of the k-points in the reduced BZ
886 : !! related to a k-point in the full BZ
887 : !! kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ
888 : !! mband = maximum number of bands
889 : !! mkmem = maximum number of k points which can fit in core memory
890 : !! mkmem_max = maximal number of k-points on each processor (MPI //)
891 : !! mk1mem = maximum number of k points for first-order WF which can fit in core memory
892 : !! mpi_enreg=MPI-parallelisation information
893 : !! mpw = maximum number of planewaves in basis sphere (large number)
894 : !! mvwtk(30,nkpt) = weights to compute the finite difference ddk
895 : !! natom = number of atoms in unit cell
896 : !! nkpt2 = number of k-points in the reduced part of the BZ
897 : !! nkpt2 = nkpt/2 in case of time-reversal symmetry (kptopt = 2)
898 : !! nkpt3 = number of k-points in the full BZ
899 : !! nneigh = total number of neighbours required to evaluate the finite difference formula
900 : !! npwarr(nkpt) = array holding npw for each k point
901 : !! nspinor = number of spinorial components of the wavefunctions
902 : !! nsppol = number of channels for spin-polarization (1 or 2)
903 : !! pwind(mpw,nneigh,mkmem) = array used to compute the overlap matrix smat between k-points
904 : !!
905 : !! OUTPUT
906 : !! d3_berry(2,3) = Berry-phase part of the third-order energy
907 : !!
908 : !! SIDE EFFECTS
909 : !! mpi_enreg=MPI-parallelisation information
910 : !!
911 : !! NOTES
912 : !! For a given set of values of i1pert,i3pert,i1dir and
913 : !! i3dir, the routine computes the k-point derivatives for
914 : !! 12dir = 1,2,3
915 : !!
916 : !! SOURCE
917 :
918 145 : subroutine pead_nl_mv(cg,cgindex,cg1,cg3,dtset,dtfil,d3_berry,gmet,&
919 145 : & i1pert,i3pert,i1dir,i3dir,kneigh,kg_neigh,kptindex,&
920 145 : & kpt3,mband,mkmem,mkmem_max,mk1mem,mpi_enreg,&
921 145 : & mpw,mvwtk,natom,nkpt2,nkpt3,nneigh,npwarr,nspinor,&
922 145 : & nsppol,pwind)
923 :
924 : use m_hide_lapack, only : dzgedi, dzgefa
925 :
926 : !Arguments ------------------------------------
927 : !
928 : !--- Arguments : integer scalars
929 : integer, intent(in) :: i1dir,i1pert,i3dir,i3pert,mband,mk1mem
930 : integer, intent(in) :: mkmem,mkmem_max,mpw,natom
931 : integer, intent(in) :: nkpt2,nkpt3,nneigh,nspinor,nsppol
932 : !
933 : !--- Arguments : integer arrays
934 : integer, intent(in) :: cgindex(nkpt2,nsppol)
935 : integer, intent(in) :: kneigh(30,nkpt2),kg_neigh(30,nkpt2,3),kptindex(2,nkpt3)
936 : integer, intent(in) :: npwarr(nkpt2),pwind(mpw,nneigh,mkmem)
937 : !
938 : !--- Arguments : real(dp) scalars
939 : !
940 : !--- Arguments : real(dp) arrays
941 : real(dp), intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
942 : real(dp), intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol)
943 : real(dp), intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol)
944 : real(dp), intent(in) :: gmet(3,3),kpt3(3,nkpt3)
945 : real(dp), intent(in) :: mvwtk(30,nkpt2)
946 : real(dp), intent(out) :: d3_berry(2,3)
947 : !
948 : !--- Arguments : structured datatypes
949 : type(MPI_type), intent(in) :: mpi_enreg
950 : type(datafiles_type), intent(in) :: dtfil
951 : type(dataset_type), intent(in) :: dtset
952 :
953 : !Local variables-------------------------------
954 : !
955 : !---- Local variables : integer scalars
956 : integer :: count,counter,count1,iband,icg
957 : integer :: ierr,ii,ikpt,ikpt_loc,ikpt2
958 : integer :: ikpt_rbz,ineigh,info,ipw,isppol,jband,jcg,jj,jkpt,job,jpw, jkpt2, jkpt_rbz
959 : integer :: lband,lpband,nband_occ,npw_k,npw_k1,my_source,his_source,dest,tag
960 : integer :: spaceComm
961 : integer,parameter :: level=52
962 : integer :: bdtot_index
963 : !
964 : !---- Local variables : integer arrays
965 145 : integer,allocatable :: ipvt(:)
966 145 : integer, allocatable :: bd_index(:,:)
967 : !
968 : !---- Local variables : real(dp) scalars
969 : real(dp) :: dotnegi,dotnegr,dotposi,dotposr
970 : ! real(dp) :: c1,c2 ! appear commented out below
971 : !
972 : !---- Local variables : real(dp) arrays
973 : real(dp) :: d3_aux(2,3),det(2,2),dk(3),dk_(3)
974 : real(dp) :: z1(2),z2(2)
975 145 : real(dp),allocatable :: buffer(:,:),cgq(:,:),cg1q(:,:),cg3q(:,:)
976 145 : real(dp),allocatable :: qmat(:,:,:),s13mat(:,:,:),s1mat(:,:,:),s3mat(:,:,:)
977 145 : real(dp),allocatable :: smat(:,:,:),zgwork(:,:)
978 : !
979 : !---- Local variables : character variables
980 : character(len=500) :: message
981 : !
982 : !---- Local variables : structured datatypes
983 :
984 :
985 : #if defined HAVE_MPI
986 : integer :: status1(MPI_STATUS_SIZE)
987 145 : spaceComm=mpi_enreg%comm_cell
988 : #endif
989 :
990 : ABI_UNUSED(dtfil%ireadwf)
991 :
992 : ! ***********************************************************************
993 :
994 145 : write(message,'(8a)') ch10,&
995 145 : & ' pead_nl_mv : finite difference expression of the k-point derivative',ch10,&
996 145 : & ' is performed using the PEAD formulation of ',&
997 145 : & 'the third-order energy',ch10,&
998 290 : & ' (see Nunes and Gonze PRB 63, 155107 (2001) [[cite:Nunes2001]] Eq. 102)',ch10
999 : !call wrtout(ab_out,message,'COLL')
1000 145 : call wrtout(std_out, message,'COLL')
1001 :
1002 :
1003 : !fab: I think that the following restriction must be eliminated:
1004 : !isppol = 1
1005 :
1006 145 : ikpt_loc = 0
1007 145 : d3_aux(:,:) = 0_dp
1008 :
1009 580 : ABI_MALLOC(s13mat,(2,mband,mband))
1010 435 : ABI_MALLOC(smat,(2,mband,mband))
1011 435 : ABI_MALLOC(s1mat,(2,mband,mband))
1012 435 : ABI_MALLOC(qmat,(2,mband,mband))
1013 435 : ABI_MALLOC(ipvt,(mband))
1014 435 : ABI_MALLOC(s3mat,(2,mband,mband))
1015 435 : ABI_MALLOC(zgwork,(2,mband))
1016 580 : ABI_MALLOC(bd_index, (nkpt2, nsppol))
1017 :
1018 318 : bdtot_index = 0
1019 318 : do isppol = 1, nsppol
1020 4815 : do ikpt_rbz = 1, nkpt2
1021 4497 : bd_index(ikpt_rbz,isppol) = bdtot_index
1022 4670 : bdtot_index = bdtot_index + dtset%nband(ikpt_rbz+nkpt2*(isppol-1))
1023 : end do
1024 : end do
1025 :
1026 : !fab: I think here I have to add the loop over spin
1027 :
1028 318 : do isppol = 1, nsppol
1029 :
1030 : ! Loop over k-points
1031 : ! COMMENT: Every processor has to make mkmem_max iterations
1032 : ! even if mkmem < mkemem_max. This is due to the fact
1033 : ! that it still has to communicate its wavefunctions
1034 : ! to other processors even if it has no more overlap
1035 : ! matrices to compute.
1036 :
1037 173 : ikpt_loc = 0 ; ikpt = 0
1038 :
1039 4591 : do while (ikpt_loc < mkmem_max)
1040 :
1041 4273 : if (ikpt_loc < mkmem) ikpt = ikpt + 1
1042 :
1043 : if (xmpi_paral == 1) then
1044 : ! if ((minval(abs(mpi_enreg%proc_distrb(ikpt,1:mband,1:dtset%nsppol) &
1045 : ! & - mpi_enreg%me)) /= 0).and.(ikpt_loc < mkmem)) cycle
1046 4273 : if(ikpt>nkpt2)then
1047 : ikpt_loc=mkmem_max
1048 : cycle
1049 : end if
1050 4273 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me)) then
1051 224 : if(ikpt==nkpt2) ikpt_loc=mkmem_max
1052 : cycle
1053 : end if
1054 : end if
1055 :
1056 4049 : ikpt_loc = ikpt_loc + 1
1057 4049 : npw_k = npwarr(ikpt)
1058 4049 : counter = 100*ikpt
1059 :
1060 4049 : ii = cgindex(ikpt,isppol)
1061 :
1062 : ! Loop on the neighbours
1063 :
1064 28516 : do ineigh = 1,nneigh
1065 :
1066 1271622 : s13mat(:,:,:) = zero
1067 1271622 : smat(:,:,:) = zero
1068 1271622 : s1mat(:,:,:) = zero
1069 1271622 : s3mat(:,:,:) = zero
1070 1271622 : qmat(:,:,:) = zero
1071 :
1072 24294 : ikpt2 = kneigh(ineigh,ikpt)
1073 24294 : ikpt_rbz = kptindex(1,ikpt2) ! index of the k-point in the reduced BZ
1074 24294 : jj = cgindex(ikpt_rbz,isppol)
1075 : ! previous fixed value for nband_k now called nband_occ:
1076 : !nband_occ = dtset%nband(ikpt_rbz+nkpt2*(isppol-1))
1077 : ! TODO: check if all these bands are occupied in nsppol = 2 case
1078 24294 : nband_occ = 0
1079 120510 : do iband = 1, dtset%nband(ikpt_rbz+nkpt2*(isppol-1))
1080 : !Note, only one image is allowed here (or occ_orig should be the same or all images)
1081 120510 : if (dtset%occ_orig(bd_index(ikpt_rbz,isppol) + iband,1) > tol10) nband_occ = nband_occ + 1
1082 : end do
1083 24294 : npw_k1 = npwarr(ikpt_rbz)
1084 97176 : dk_(:) = kpt3(:,ikpt2) - dtset%kptns(:,ikpt)
1085 97176 : dk(:) = dk_(:) - nint(dk_(:)) + real(kg_neigh(ineigh,ikpt,:),dp)
1086 :
1087 24294 : count = nspinor*mband*npw_k1
1088 72882 : ABI_MALLOC(cgq,(2,count))
1089 48588 : ABI_MALLOC(cg1q,(2,count))
1090 48588 : ABI_MALLOC(cg3q,(2,count))
1091 :
1092 : #if defined HAVE_MPI
1093 :
1094 24294 : my_source = mpi_enreg%proc_distrb(ikpt_rbz,1,1)
1095 :
1096 : ! do dest = 0, mpi_enreg%nproc-1
1097 31836066 : do dest = 0, maxval(mpi_enreg%proc_distrb(1:nkpt2,1:mband,1:dtset%nsppol))
1098 :
1099 51276 : if ((dest==mpi_enreg%me).and.(ikpt_loc <= mkmem)) then
1100 : ! I am dest and have something to do
1101 :
1102 24294 : if (my_source == mpi_enreg%me) then
1103 : ! I am destination and source
1104 23510 : jcg = cgindex(ikpt_rbz,isppol)
1105 :
1106 24919310 : cgq(:,1:count) = cg(:,jcg+1:jcg+count)
1107 24919310 : cg1q(:,1:count) = cg1(:,jcg+1:jcg+count)
1108 24919310 : cg3q(:,1:count) = cg3(:,jcg+1:jcg+count)
1109 :
1110 : else
1111 : ! I am the destination but not the source -> receive
1112 :
1113 784 : tag = ikpt_rbz
1114 :
1115 2352 : ABI_MALLOC(buffer,(2,3*count))
1116 :
1117 784 : call MPI_RECV(buffer,2*3*count,MPI_DOUBLE_PRECISION,my_source,tag,spaceComm,status1,ierr)
1118 :
1119 141904 : cgq(:,1:count) = buffer(:,1:count)
1120 141904 : cg1q(:,1:count) = buffer(:,count+1:2*count)
1121 141904 : cg3q(:,1:count) = buffer(:,2*count+1:3*count)
1122 784 : ABI_FREE(buffer)
1123 :
1124 : end if
1125 :
1126 2688 : else if (ikpt_loc <= mpi_enreg%mkmem(dest)) then ! dest != me and the dest has a k-point to treat
1127 :
1128 2688 : jkpt=mpi_enreg%kpt_loc2ibz_sp(dest, ikpt_loc,1)
1129 2688 : jkpt2 = kneigh(ineigh,jkpt)
1130 2688 : jkpt_rbz = kptindex(1,jkpt2) ! index of the k-point in the reduced BZ
1131 :
1132 2688 : his_source = mpi_enreg%proc_distrb(jkpt_rbz,1,1)
1133 :
1134 2688 : if (his_source == mpi_enreg%me) then
1135 :
1136 784 : jcg = cgindex(jkpt_rbz,isppol)
1137 :
1138 784 : tag = jkpt_rbz
1139 784 : count1 = npwarr(jkpt_rbz)*mband*nspinor
1140 2352 : ABI_MALLOC(buffer,(2,3*count1))
1141 141904 : buffer(:,1:count1) = cg(:,jcg+1:jcg+count1)
1142 141904 : buffer(:,count1+1:2*count1) = cg1(:,jcg+1:jcg+count1)
1143 141904 : buffer(:,2*count1+1:3*count1) = cg3(:,jcg+1:jcg+count1)
1144 :
1145 784 : call MPI_SEND(buffer,2*3*count1,MPI_DOUBLE_PRECISION,dest,tag,spaceComm,ierr)
1146 :
1147 784 : ABI_FREE(buffer)
1148 :
1149 : end if
1150 :
1151 : end if
1152 :
1153 : end do
1154 : !
1155 : ! do jkpt = 1, nkpt2
1156 : !
1157 : ! if ((jkpt == ikpt_rbz).and.(source /= mpi_enreg%me).and.&
1158 : ! & (ikpt_loc <= mkmem)) then
1159 : !
1160 : ! tag = jkpt
1161 : !
1162 : ! allocate(buffer(2,3*count))
1163 : ! call MPI_RECV(buffer,2*3*count,MPI_DOUBLE_PRECISION,&
1164 : ! source,tag,spaceComm,status1,ierr)
1165 : !
1166 : ! cgq(:,1:count) = buffer(:,1:count)
1167 : ! cg1q(:,1:count) = buffer(:,count+1:2*count)
1168 : ! cg3q(:,1:count) = buffer(:,2*count+1:3*count)
1169 : ! deallocate(buffer)
1170 : !
1171 : ! end if
1172 : !
1173 : ! ! ----------------------------------------------------------------------------
1174 : ! ! --------------- Here: send the WF to all the cpus that need it -------------
1175 : ! ! ----------------------------------------------------------------------------
1176 : !
1177 : ! do dest = 1, mpi_enreg%nproc
1178 : !
1179 : ! if ((minval(abs(mpi_enreg%proc_distrb(jkpt,1:mband,1:dtset%nsppol) &
1180 : ! & - mpi_enreg%me)) == 0).and.&
1181 : ! & (mpi_enreg%kptdstrb(dest,ineigh,ikpt_loc) == jkpt)) then
1182 : !
1183 : !
1184 : !
1185 : ! jcg = cgindex(jkpt,isppol)
1186 : !
1187 : ! if (((dest-1) == mpi_enreg%me)) then
1188 : !
1189 : ! cgq(:,1:count) = cg(:,jcg+1:jcg+count)
1190 : ! cg1q(:,1:count) = cg1(:,jcg+1:jcg+count)
1191 : ! cg3q(:,1:count) = cg3(:,jcg+1:jcg+count)
1192 : !
1193 : ! else
1194 : !
1195 : ! tag = jkpt
1196 : ! count1 = npwarr(jkpt)*mband*nspinor
1197 : ! allocate(buffer(2,3*count1))
1198 : ! buffer(:,1:count1) = cg(:,jcg+1:jcg+count1)
1199 : ! buffer(:,count1+1:2*count1) = cg1(:,jcg+1:jcg+count1)
1200 : ! buffer(:,2*count1+1:3*count1) = cg3(:,jcg+1:jcg+count1)
1201 : !
1202 : ! call MPI_SEND(buffer,2*3*count1,MPI_DOUBLE_PRECISION,(dest-1),tag,spaceComm,ierr)
1203 : !
1204 : ! deallocate(buffer)
1205 : !
1206 : ! end if
1207 : !
1208 : ! end if
1209 : !
1210 : ! end do ! loop over dest
1211 : !
1212 : ! end do ! loop over jkpt
1213 :
1214 24294 : if (ikpt_loc > mkmem) then
1215 0 : ABI_FREE(cgq)
1216 0 : ABI_FREE(cg1q)
1217 0 : ABI_FREE(cg3q)
1218 0 : cycle
1219 : end if
1220 :
1221 : #else
1222 : ! no // over k-points
1223 :
1224 : cgq(:,1:count) = cg(:,jj+1:jj+count)
1225 : cg1q(:,1:count) = cg1(:,jj+1:jj+count)
1226 : cg3q(:,1:count) = cg3(:,jj+1:jj+count)
1227 :
1228 : #endif
1229 :
1230 : ! Compute overlap matrices
1231 :
1232 24294 : if (kptindex(2,ikpt2) == 0) then ! no time-reversal symmetry
1233 :
1234 1827044 : do ipw = 1, npw_k
1235 :
1236 1807602 : jpw = pwind(ipw,ineigh,ikpt_loc)
1237 1827044 : if (jpw /= 0) then
1238 :
1239 7912824 : do iband = 1, nband_occ
1240 31954392 : do jband = 1, nband_occ
1241 :
1242 24041568 : icg = ii + (iband-1)*npw_k + ipw
1243 24041568 : jcg = (jband-1)*npw_k1 + jpw
1244 :
1245 : smat(1,iband,jband) = smat(1,iband,jband) + &
1246 24041568 : & cg(1,icg)*cgq(1,jcg) + cg(2,icg)*cgq(2,jcg)
1247 : smat(2,iband,jband) = smat(2,iband,jband) + &
1248 24041568 : & cg(1,icg)*cgq(2,jcg) - cg(2,icg)*cgq(1,jcg)
1249 :
1250 : s13mat(1,iband,jband) = s13mat(1,iband,jband) + &
1251 24041568 : & cg1(1,icg)*cg3q(1,jcg) + cg1(2,icg)*cg3q(2,jcg)
1252 : s13mat(2,iband,jband) = s13mat(2,iband,jband) + &
1253 24041568 : & cg1(1,icg)*cg3q(2,jcg) - cg1(2,icg)*cg3q(1,jcg)
1254 :
1255 : s1mat(1,iband,jband) = s1mat(1,iband,jband) + &
1256 : & cg1(1,icg)*cgq(1,jcg) + cg1(2,icg)*cgq(2,jcg) + &
1257 24041568 : & cg(1,icg)*cg1q(1,jcg) + cg(2,icg)*cg1q(2,jcg)
1258 : s1mat(2,iband,jband) = s1mat(2,iband,jband) + &
1259 : & cg1(1,icg)*cgq(2,jcg) - cg1(2,icg)*cgq(1,jcg) + &
1260 24041568 : & cg(1,icg)*cg1q(2,jcg) - cg(2,icg)*cg1q(1,jcg)
1261 :
1262 : s3mat(1,iband,jband) = s3mat(1,iband,jband) + &
1263 : & cg3(1,icg)*cgq(1,jcg) + cg3(2,icg)*cgq(2,jcg) + &
1264 24041568 : & cg(1,icg)*cg3q(1,jcg) + cg(2,icg)*cg3q(2,jcg)
1265 : s3mat(2,iband,jband) = s3mat(2,iband,jband) + &
1266 : & cg3(1,icg)*cgq(2,jcg) - cg3(2,icg)*cgq(1,jcg) + &
1267 30266304 : & cg(1,icg)*cg3q(2,jcg) - cg(2,icg)*cg3q(1,jcg)
1268 :
1269 : end do
1270 : end do
1271 :
1272 : end if
1273 :
1274 : end do ! ipw
1275 :
1276 : else ! use time-reversal symmetry
1277 :
1278 423940 : do ipw = 1,npw_k
1279 :
1280 419088 : jpw = pwind(ipw,ineigh,ikpt_loc)
1281 423940 : if (jpw /= 0) then
1282 :
1283 1901116 : do iband = 1, nband_occ
1284 7881324 : do jband = 1, nband_occ
1285 :
1286 5980208 : icg = ii + (iband-1)*npw_k + ipw
1287 5980208 : jcg = (jband-1)*npw_k1 + jpw
1288 :
1289 : smat(1,iband,jband) = smat(1,iband,jband) + &
1290 5980208 : & cg(1,icg)*cgq(1,jcg) - cg(2,icg)*cgq(2,jcg)
1291 : smat(2,iband,jband) = smat(2,iband,jband) - &
1292 5980208 : & cg(1,icg)*cgq(2,jcg) - cg(2,icg)*cgq(1,jcg)
1293 :
1294 : s13mat(1,iband,jband) = s13mat(1,iband,jband) + &
1295 5980208 : & cg1(1,icg)*cg3q(1,jcg) - cg1(2,icg)*cg3q(2,jcg)
1296 : s13mat(2,iband,jband) = s13mat(2,iband,jband) - &
1297 5980208 : & cg1(1,icg)*cg3q(2,jcg) - cg1(2,icg)*cg3q(1,jcg)
1298 :
1299 : s1mat(1,iband,jband) = s1mat(1,iband,jband) + &
1300 : & cg1(1,icg)*cgq(1,jcg) - cg1(2,icg)*cgq(2,jcg) + &
1301 5980208 : & cg(1,icg)*cg1q(1,jcg) - cg(2,icg)*cg1q(2,jcg)
1302 : s1mat(2,iband,jband) = s1mat(2,iband,jband) - &
1303 : & cg1(1,icg)*cgq(2,jcg) - cg1(2,icg)*cgq(1,jcg) - &
1304 5980208 : & cg(1,icg)*cg1q(2,jcg) - cg(2,icg)*cg1q(1,jcg)
1305 :
1306 : s3mat(1,iband,jband) = s3mat(1,iband,jband) + &
1307 : & cg3(1,icg)*cgq(1,jcg) - cg3(2,icg)*cgq(2,jcg) + &
1308 5980208 : & cg(1,icg)*cg3q(1,jcg) - cg(2,icg)*cg3q(2,jcg)
1309 : s3mat(2,iband,jband) = s3mat(2,iband,jband) - &
1310 : & cg3(1,icg)*cgq(2,jcg) - cg3(2,icg)*cgq(1,jcg) - &
1311 7492576 : & cg(1,icg)*cg3q(2,jcg) - cg(2,icg)*cg3q(1,jcg)
1312 :
1313 : end do
1314 : end do
1315 :
1316 : end if
1317 :
1318 : end do ! ipw
1319 :
1320 : end if
1321 :
1322 24294 : ABI_FREE(cgq)
1323 24294 : ABI_FREE(cg1q)
1324 24294 : ABI_FREE(cg3q)
1325 :
1326 : ! Compute qmat, the inverse of smat
1327 :
1328 24294 : job = 1 ! compute inverse only
1329 1271622 : qmat(:,:,:) = smat(:,:,:)
1330 :
1331 24294 : call dzgefa(qmat,mband,nband_occ,ipvt,info)
1332 24294 : call dzgedi(qmat,mband,nband_occ,ipvt,det,zgwork,job)
1333 :
1334 : ! DEBUG
1335 : ! write(100,*)
1336 : ! write(100,*)'ikpt = ',ikpt,'ineigh = ',ineigh
1337 : ! do iband = 1,nband_occ
1338 : ! do jband = 1,nband_occ
1339 : ! c1 = 0_dp ; c2 = 0_dp
1340 : ! do lband = 1,nband_occ
1341 : ! c1 = c1 + smat(1,iband,lband)*qmat(1,lband,jband) - &
1342 : ! & smat(2,iband,lband)*qmat(2,lband,jband)
1343 : ! c2 = c2 + smat(1,iband,lband)*qmat(2,lband,jband) + &
1344 : ! & smat(2,iband,lband)*qmat(1,lband,jband)
1345 : ! end do
1346 : ! write(100,'(2(2x,i2),2(2x,f16.9))')iband,jband,&
1347 : ! & c1,c2
1348 : ! end do
1349 : ! end do
1350 : ! ENDDEBUG
1351 :
1352 :
1353 :
1354 : ! Accumulate sum over bands
1355 :
1356 24294 : dotposr = 0_dp ; dotposi = 0_dp
1357 24294 : dotnegr = 0_dp ; dotnegi = 0_dp
1358 120510 : do iband = 1, nband_occ
1359 503814 : do jband = 1, nband_occ
1360 :
1361 : dotposr = dotposr + &
1362 : & s13mat(1,iband,jband)*qmat(1,jband,iband) - &
1363 383304 : & s13mat(2,iband,jband)*qmat(2,jband,iband)
1364 : dotposi = dotposi + &
1365 : & s13mat(1,iband,jband)*qmat(2,jband,iband) + &
1366 383304 : & s13mat(2,iband,jband)*qmat(1,jband,iband)
1367 :
1368 :
1369 2009976 : do lband = 1, nband_occ
1370 8030424 : do lpband= 1, nband_occ
1371 :
1372 : z1(1) = s1mat(1,iband,jband)*qmat(1,jband,lband) - &
1373 6116664 : & s1mat(2,iband,jband)*qmat(2,jband,lband)
1374 : z1(2) = s1mat(1,iband,jband)*qmat(2,jband,lband) + &
1375 6116664 : & s1mat(2,iband,jband)*qmat(1,jband,lband)
1376 :
1377 : z2(1) = s3mat(1,lband,lpband)*qmat(1,lpband,iband) - &
1378 6116664 : & s3mat(2,lband,lpband)*qmat(2,lpband,iband)
1379 : z2(2) = s3mat(1,lband,lpband)*qmat(2,lpband,iband) + &
1380 6116664 : & s3mat(2,lband,lpband)*qmat(1,lpband,iband)
1381 :
1382 : dotnegr = dotnegr + &
1383 6116664 : & z1(1)*z2(1) - z1(2)*z2(2)
1384 : dotnegi = dotnegi + &
1385 7647120 : & z1(1)*z2(2) + z1(2)*z2(1)
1386 :
1387 : end do ! lpband
1388 : end do ! lband
1389 :
1390 : end do ! jband
1391 : end do ! iband
1392 :
1393 : d3_aux(1,:) = d3_aux(1,:) + &
1394 97176 : & dk(:)*mvwtk(ineigh,ikpt)*dtset%wtk(ikpt)*(2_dp*dotposr-dotnegr)
1395 : d3_aux(2,:) = d3_aux(2,:) + &
1396 101225 : & dk(:)*mvwtk(ineigh,ikpt)*dtset%wtk(ikpt)*(2_dp*dotposi-dotnegi)
1397 :
1398 : end do ! End loop over neighbours
1399 :
1400 :
1401 : end do ! End loop over k-points
1402 :
1403 : end do ! fab: end loop over spin
1404 :
1405 :
1406 :
1407 :
1408 145 : call xmpi_sum(d3_aux,spaceComm,ierr)
1409 :
1410 :
1411 145 : ABI_FREE(s13mat)
1412 145 : ABI_FREE(smat)
1413 145 : ABI_FREE(s1mat)
1414 145 : ABI_FREE(qmat)
1415 145 : ABI_FREE(ipvt)
1416 145 : ABI_FREE(s3mat)
1417 145 : ABI_FREE(zgwork)
1418 145 : ABI_FREE(bd_index)
1419 :
1420 :
1421 : !fab: I think that in the following we have to make a distinction:
1422 : !for the spin unpolarized case we leave the PEAD expression as it is, while
1423 : !in the spin polarized case we have simply to divide by 2
1424 : !(see eq.19 di PRB 63,155107 [[cite:Nunes2001]], eq. 7 di PRB 71,125107 [[cite:Veithen2005]]
1425 : ! and eq 13 di PRB 71, 125107 [[cite:Veithen2005]] ...
1426 : !in this latter equation the 2 must be simply replaced by the sum over the spin components...
1427 : !and indeed we have inserted the loop over the spin,
1428 : !but there was a factor 2 already present in the routine due to spin degenracy that had to be removed)
1429 :
1430 :
1431 145 : if (nsppol==1) then
1432 :
1433 : ! Take minus the imaginary part
1434 :
1435 468 : d3_berry(1,:) = -1_dp*d3_aux(2,:)
1436 468 : d3_berry(2,:) = d3_aux(1,:)
1437 :
1438 468 : d3_berry(2,:) = 0_dp
1439 :
1440 : else
1441 :
1442 112 : d3_berry(1,:) = -1_dp*d3_aux(2,:)/2._dp
1443 112 : d3_berry(2,:) = d3_aux(1,:)/2._dp
1444 :
1445 112 : d3_berry(2,:) = 0_dp/2._dp
1446 :
1447 : end if
1448 :
1449 : !DEBUG
1450 : !write(100,*)'pead_nl_mv.f : d3_berry'
1451 : !write(100,*)'Perturbation',i1dir,i3dir
1452 : !write(100,*)
1453 : !write(100,*)'before transformation'
1454 : !write(100,*)'real part'
1455 : !write(100,'(3(2x,f20.9))')d3_berry(1,:)
1456 : !write(100,*)
1457 : !write(100,*)'imaginary part'
1458 : !write(100,'(3(2x,f20.9))')d3_berry(2,:)
1459 : !write(100,*)
1460 : !write(100,*)'after transformation'
1461 : !ENDDEBUG
1462 :
1463 : !Compute the projection on the basis vectors of
1464 : !reciprocal space
1465 :
1466 145 : d3_aux(:,:) = 0_dp
1467 580 : do ii = 1,3
1468 1885 : do jj = 1,3
1469 4350 : d3_aux(:,ii) = d3_aux(:,ii) + gmet(ii,jj)*d3_berry(:,jj)
1470 : end do
1471 : end do
1472 145 : d3_berry(:,:) = d3_aux(:,:)
1473 :
1474 : !Write out the berryphase part of the third order energy
1475 :
1476 145 : if (mpi_enreg%me == 0) then
1477 :
1478 117 : write(message,'(a,a,a)')ch10,' Berryphase part of the third-order energy:',ch10
1479 117 : call wrtout(std_out, message,'COLL')
1480 :
1481 117 : if (i1pert < natom + 1) then
1482 : write(message,'(a,i3,a,i3)')&
1483 76 : & ' j1: Displacement of atom ',i1pert,&
1484 152 : & ' along direction ',i1dir
1485 41 : else if (i1pert == natom + 2) then
1486 : write(message,'(a,i3)')&
1487 41 : & ' j1: homogenous electric field along direction ',i1dir
1488 : end if
1489 117 : call wrtout(std_out, message,'COLL')
1490 :
1491 : write(message,'(a)')&
1492 117 : & ' j2: k-point derivative along direction i2dir '
1493 117 : call wrtout(std_out, message,'COLL')
1494 :
1495 117 : if (i3pert < natom + 1) then
1496 : write(message,'(a,i3,a,i3,a)')&
1497 0 : & ' j3: Displacement of atom ',i3pert,&
1498 0 : & ' along direction ',i3dir,ch10
1499 117 : else if (i3pert == natom + 2) then
1500 : write(message,'(a,i3,a)')&
1501 117 : & ' j3: homogenous electric field along direction ',i3dir,ch10
1502 : end if
1503 117 : call wrtout(std_out, message,'COLL')
1504 :
1505 : ! write(ab_out,'(5x,a5,8x,a9,5x,a14)')'i2dir','real part','imaginary part'
1506 117 : write(std_out,'(5x,a5,8x,a9,5x,a14)')'i2dir','real part','imaginary part'
1507 468 : do ii = 1,3
1508 351 : write(std_out,'(7x,i1,3x,f16.9,3x,f16.9)')ii,&
1509 702 : & d3_berry(1,ii),d3_berry(2,ii)
1510 351 : write(std_out,'(7x,i1,3x,f16.9,3x,f16.9)')ii,&
1511 819 : & d3_berry(1,ii),d3_berry(2,ii)
1512 : end do
1513 :
1514 : end if ! mpi_enreg%me == 0
1515 :
1516 : !DEBUG
1517 : !write(100,*)'real part'
1518 : !write(100,'(3(2x,f20.9))')d3_berry(1,:)
1519 : !write(100,*)
1520 : !write(100,*)'imaginary part'
1521 : !write(100,'(3(2x,f20.9))')d3_berry(2,:)
1522 : !ENDDEBUG
1523 :
1524 145 : end subroutine pead_nl_mv
1525 : !!***
1526 :
1527 : end module m_pead_nl_loop
1528 : !!***
|
