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