Line data Source code
1 : !!****m* ABINIT/m_paw_mkaewf
2 : !! NAME
3 : !! m_paw_mkaewf
4 : !!
5 : !! FUNCTION
6 : !! Construct complete AE wave functions on the fine FFT grid adding onsite PAW corrections.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
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_paw_mkaewf
23 :
24 : use defs_basis
25 : use defs_wvltypes
26 : use m_abicore
27 : use m_xmpi
28 : use m_hide_blas
29 : use m_splines
30 : use m_errors
31 : use m_nctk
32 : use m_hdr
33 : use m_dtset
34 : use m_dtfil
35 : use netcdf
36 :
37 : use defs_abitypes, only : MPI_type
38 : use m_io_tools, only : flush_unit
39 : use m_numeric_tools, only : wrap2_zero_one
40 : use m_fftcore, only : sphereboundary
41 : use m_geometry, only : xcart2xred
42 : use m_crystal, only : crystal_t
43 : use m_ebands, only : ebands_t
44 : use m_pawrad, only : pawrad_type
45 : use m_pawtab, only : pawtab_type, pawtab_get_lsize
46 : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print
47 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free
48 : use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
49 : use m_paral_atom, only : get_my_atmtab, free_my_atmtab
50 : use m_paw_nhat, only : nhatgrid
51 : use m_mpinfo, only : proc_distrb_cycle
52 : use m_fft, only : fourwf
53 :
54 : implicit none
55 :
56 : private
57 :
58 : public :: pawmkaewf
59 :
60 : CONTAINS !========================================================================================
61 : !!***
62 :
63 : !!****f* m_paw_mkaewf/pawmkaewf
64 : !! NAME
65 : !! pawmkaewf
66 : !!
67 : !! FUNCTION
68 : !! Construct complete AE wave functions on the fine FFT grid adding onsite PAW corrections.
69 : !!
70 : !! INPUTS
71 : !! crystal<crystal_t>=Crystalline structure
72 : !! ebands<ebands_t>=Electronic energies
73 : !! dimcprj(natom)=array of dimensions of array cprj (not ordered)
74 : !! mband=maximum number of bands
75 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
76 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
77 : !! mkmem=number of k points treated by this node.
78 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
79 : !! [comm_atom]= MPI communicator over atoms
80 : !! mpw=maximum dimensioned size of npw.
81 : !! my_natom=number of atoms treated by current processor
82 : !! natom=number of atoms in cell
83 : !! ntypat=number of types of atoms in the cell
84 : !! nkpt=Total number of k-points
85 : !! nsppol=1 for unpolarized, 2 for spin-polarized
86 : !! unks=unit number for G vectors.
87 : !! nband(nkpt*nsppol)=Number of bands for each k-point and spin.
88 : !! istwfk(nkpt)=Storage mode at each k-point.
89 : !! Pawfgrtab(natom) <type(pawfgrtab_type)> : data about the fine grid around each atom
90 : !! Pawrad(ntypat) <type(pawrad_type)> : radial mesh data for each type of atom
91 : !! Pawtab(ntypat) <type(pawtab_type)> : PAW functions around each type of atom
92 : !! Dtfil <type(datafiles_type)>=variables related to files
93 : !! cg(2,mcg)=planewave coefficients of wavefunctions.
94 : !! Cprj(natom,nspinor*mband*mkmem*nsppol)=<p_lmn|Cnk> coefficients for each WF |Cnk>
95 : !! and each |p_lmn> non-local projector
96 : !! npwarr(nkpt)=Number of plane waves at each k-point
97 : !! ngfftf(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
98 : !! Note that ngfftf refers to the fine mesh.
99 : !! kg(3,mpw*mkmem)=reduced planewave coordinates
100 : !! Hdr<hdr_type>=the header of wf, den and pot files
101 : !! kpt(3,nkpt)=reduced coordinates of k points.
102 : !!
103 : !! OUTPUT
104 : !! ierr=Status error
105 : !! Main output is written on file (ETSF_IO file format).
106 : !!
107 : !! NOTES
108 : !! In PAW calculations, the pseudized wavefunction us represented
109 : !! on a relatively small plane wave basis set and is not normalized
110 : !! as it does not include the on-site PAW contributions which is described
111 : !! in terms of real spherical harmonics and radial functions.
112 : !! For post-processing and proper visualization, it is necessary
113 : !! to use the full electronic wave function, which is what this subroutine constructs.
114 : !! Specifically, it computes the pseudo part by doing an FFT from G- to r-space
115 : !! using the dense mesh defined by pawecutdg. The on-site PAW terms are also
116 : !! computed in real space inside each sphere and added to the pseudo part.
117 : !! Notice that this formula is expressed on the fine grid, and requires
118 : !! interpolating the PAW radial functions onto this grid, as well as calling
119 : !! initylmr in order to get the angular functions on the grid points.
120 : !!
121 : !! SOURCE
122 :
123 2 : subroutine pawmkaewf(Dtset,crystal,ebands,my_natom,mpw,mband,mcg,mcprj,nkpt,mkmem,nsppol,nband,&
124 2 : & istwfk,npwarr,kpt,ngfftf,kg,dimcprj,Pawfgrtab,Pawrad,Pawtab,&
125 1 : & Hdr,Dtfil,cg,Cprj,MPI_enreg,ierr,pseudo_norms,set_k,set_band , &
126 1 : & mpi_atmtab,comm_atom) ! Optional arguments
127 :
128 : !Arguments ------------------------------------
129 : !scalars
130 : integer,intent(in) :: my_natom,mband,mcg,mcprj,mkmem,mpw,nsppol,nkpt
131 : integer,intent(in),optional :: comm_atom,set_k,set_band
132 : integer,intent(out) :: ierr
133 : type(Datafiles_type),intent(in) :: Dtfil
134 : type(MPI_type),intent(in) :: MPI_enreg
135 : type(hdr_type),intent(inout) :: Hdr
136 : type(dataset_type),intent(in) :: Dtset
137 : type(crystal_t),intent(in) :: crystal
138 : type(ebands_t),intent(in) :: ebands
139 : !arrays
140 : integer,intent(in) :: nband(nkpt*nsppol),istwfk(nkpt),npwarr(nkpt),dimcprj(crystal%natom)
141 : integer,intent(in) :: ngfftf(18),kg(3,mpw*mkmem)
142 : integer,optional,target,intent(in) :: mpi_atmtab(:)
143 : real(dp),intent(in) :: cg(2,mcg)
144 : real(dp),intent(in) :: kpt(3,nkpt)
145 : real(dp),optional,intent(out) :: pseudo_norms(nsppol,nkpt,mband)
146 : type(pawfgrtab_type),intent(in) :: Pawfgrtab(my_natom)
147 : type(pawrad_type),intent(in) :: Pawrad(crystal%ntypat)
148 : type(pawtab_type),intent(in) :: Pawtab(crystal%ntypat)
149 : type(pawcprj_type),intent(in) :: Cprj(crystal%natom,mcprj)
150 :
151 : !Local variables-------------------------------
152 : !scalars
153 : integer,parameter :: tim_fourwf0=0,tim_rwwf0=0,master=0
154 : integer :: bdtot_index,iband,icg,mgfftf,paral_kgb
155 : integer :: iatom,iatom_tot,ifgd,ifftsph,ifft,itypat,ispinor,ipw,ndat,ii,i1,i2,i3
156 : integer :: jl,jm,jlmn,natom
157 : integer :: max_nfgd,nfgd,ln_size,lmn_size,my_comm_atom,option
158 : integer :: iorder_cprj,comm_cell,me_kpt,ibsp,ibg,isppol,ikpt,nband_k,cplex
159 : integer :: n1,n2,n3,n4,n5,n6,ikg,npwout,istwf_k,npw_k
160 : integer :: nfftot,nprocs,tmp_unt
161 : integer :: optcut,optgr0,optgr1,optgr2,optrad,start_band,start_kpt,stop_kpt,stop_band
162 : logical :: my_atmtab_allocated,paral_atom
163 : real(dp),parameter :: weight1=one
164 : real(dp) :: phj,tphj,re_p,im_p,norm,norm_rerr,max_rerr,imur,reur,arg
165 : character(len=500) :: msg
166 : character(len=nctk_slen) :: shape_str
167 : !arrays
168 1 : integer,allocatable :: l_size_atm(:)
169 1 : integer, pointer :: my_atmtab(:)
170 1 : integer,allocatable :: gbound(:,:),kg_k(:,:)
171 : real(dp) :: red(3),shift(3),rfft(3),kpoint(3),cp_fact(2)
172 1 : real(dp),allocatable :: r0shift(:,:,:),phk_atm(:,:,:),paw_compl(:)
173 1 : real(dp),allocatable :: buf_tmp(:,:,:),fofgin(:,:),fofgin_down(:,:),fofgout(:,:)
174 2 : real(dp),allocatable :: denpot(:,:,:),fofr(:,:,:,:),fofr_down(:,:,:,:),phkr(:,:)
175 1 : real(dp),allocatable :: ur_ae(:,:), ur_pw(:,:),ur_ae_onsite(:,:),ur_ps_onsite(:,:)
176 2 : real(dp),allocatable :: ur_mask(:),dummy_1d(:),rsph_red(:,:),rsph_cart(:,:)
177 1 : type(pawcprj_type),allocatable :: Cprj_k(:,:)
178 14 : type(pawfgrtab_type) :: local_pawfgrtab(my_natom)
179 1 : type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
180 : integer :: fform,ncerr,ncid,ae_ncid,pw_ncid,aeons_ncid,psons_ncid
181 : character(len=fnlen) :: fname
182 : ! ************************************************************************
183 :
184 : DBG_ENTER("COLL")
185 :
186 : !Init parallelism
187 1 : comm_cell = MPI_enreg%comm_cell; nprocs = xmpi_comm_size(comm_cell)
188 1 : me_kpt = MPI_enreg%me_kpt; paral_kgb=mpi_enreg%paral_kgb
189 :
190 : !Compatibility tests
191 1 : ABI_CHECK(mkmem/=0, "mkmem==0 not supported anymore!")
192 1 : ABI_CHECK(MPI_enreg%paral_kgb == 0, "paral_kgb/=0 not coded")
193 1 : ABI_CHECK(SIZE(dimcprj)>0, "dimcprj should be allocated")
194 1 : ABI_CHECK(mpi_enreg%paral_spinor==0, "parallelisation over spinors not implemented")
195 1 : ABI_CHECK(nprocs==1, "k spin parallelism not yet active")
196 1 : ABI_CHECK(dtset%nspinor==1, "nspinor == 2 is buggy")
197 :
198 1 : natom = crystal%natom
199 :
200 : !Set up parallelism over atoms
201 1 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
202 1 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
203 1 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
204 1 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
205 :
206 : !If collection of pseudo norms is enabled, make sure the array is initialised
207 19 : if (present(pseudo_norms)) pseudo_norms = zero
208 :
209 : !use a local copy of pawfgrtab to make sure we use the correction in the paw spheres
210 : !the usual pawfgrtab uses r_shape which may not be the same as r_paw
211 1 : if (my_natom>0) then
212 1 : if (paral_atom) then
213 0 : call pawtab_get_lsize(pawtab,l_size_atm,my_natom,Dtset%typat,mpi_atmtab=my_atmtab)
214 : call pawfgrtab_init(local_pawfgrtab,Pawfgrtab(1)%cplex,l_size_atm,Dtset%nspden,Dtset%typat,&
215 0 : & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
216 : else
217 1 : call pawtab_get_lsize(pawtab,l_size_atm,my_natom,Dtset%typat)
218 1 : call pawfgrtab_init(local_pawfgrtab,Pawfgrtab(1)%cplex,l_size_atm,Dtset%nspden,Dtset%typat)
219 : end if
220 1 : ABI_FREE(l_size_atm)
221 : end if
222 1 : optcut = 1 ! use rpaw to construct local_pawfgrtab
223 1 : optgr0 = 0; optgr1 = 0; optgr2 = 0 ! dont need gY terms locally
224 1 : optrad = 1 ! do store r-R
225 :
226 1 : if (paral_atom) then
227 : call nhatgrid(crystal%atindx1,crystal%gmet,my_natom,natom,crystal%nattyp,ngfftf,crystal%ntypat,&
228 : & optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,crystal%rprimd,Dtset%typat,crystal%ucvol,Hdr%xred,&
229 0 : & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
230 : else
231 : call nhatgrid(crystal%atindx1,crystal%gmet,my_natom,natom,crystal%nattyp,ngfftf,crystal%ntypat,&
232 1 : & optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,crystal%rprimd,Dtset%typat,crystal%ucvol,Hdr%xred)
233 : end if
234 : !now local_pawfgrtab is ready to use
235 :
236 2 : max_nfgd=MAXVAL(local_pawfgrtab(:)%nfgd) ! MAX no. of points in the fine grid for this PAW sphere
237 4 : ABI_MALLOC(r0shift,(3,max_nfgd,my_natom))
238 4 : ABI_MALLOC(phk_atm,(2,max_nfgd,my_natom))
239 :
240 2 : do iatom=1,my_natom
241 1 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
242 :
243 1 : nfgd=local_pawfgrtab(iatom)%nfgd ! no. of points in the fine grid for this PAW sphere
244 3 : ABI_MALLOC(rsph_red,(3,nfgd))
245 2 : ABI_MALLOC(rsph_cart,(3,nfgd))
246 1430 : do ifgd=1,nfgd
247 5717 : rsph_cart(:,ifgd) = local_pawfgrtab(iatom)%rfgd(:,ifgd) + crystal%xcart(:,iatom_tot)
248 : end do
249 1 : call xcart2xred(nfgd,crystal%rprimd,rsph_cart,rsph_red) ! we work in reduced coordinates.
250 1430 : do ifgd=1,nfgd
251 1429 : call wrap2_zero_one(rsph_red(1,ifgd),red(1),shift(1)) ! num = red + shift
252 1429 : call wrap2_zero_one(rsph_red(2,ifgd),red(2),shift(2))
253 1429 : call wrap2_zero_one(rsph_red(3,ifgd),red(3),shift(3))
254 5717 : r0shift(:,ifgd,iatom) = shift
255 : !if (ANY( ABS(shift) > tol12)) then
256 : ! ABI_WARNING("rmR_red is outside the first unit cell.")
257 : ! write(std_out,*)rsph_red(:,ifgd),shift
258 : !end if
259 : end do
260 1 : ABI_FREE(rsph_red)
261 2 : ABI_FREE(rsph_cart)
262 : end do
263 :
264 1 : if (.not.paral_atom .and. my_natom>0) then
265 : call pawfgrtab_print(local_pawfgrtab,natom=natom,unit=std_out,&
266 1 : & prtvol=Dtset%prtvol,mode_paral="COLL")
267 : end if
268 :
269 1 : ierr=0
270 :
271 : !FIXME check ordering in cprj and Eventually in external file
272 : !why is iorder_cprj not stored in the file for crosschecking purpose?
273 : !Here Im assuming cprj are not ordered!
274 1 : iorder_cprj=0
275 :
276 : !n4,n5,n6 are FFT dimensions, modified to avoid cache trashing
277 1 : n1=ngfftf(1); n2=ngfftf(2); n3=ngfftf(3)
278 1 : n4=ngfftf(4); n5=ngfftf(5); n6=ngfftf(6)
279 4 : nfftot=PRODUCT(ngfftf(1:3))
280 4 : mgfftf=MAXVAL(ngfftf(1:3))
281 :
282 3 : ABI_MALLOC(phkr,(2,nfftot))
283 4 : ABI_MALLOC(gbound,(2*mgfftf+8,2))
284 :
285 : !=== Initialize ETSF_IO files ===
286 : ! FIXME: nspinor == 2 is buggy
287 1 : if(dtset%pawprtwf==1) then
288 1 : fname = trim(dtfil%filnam_ds(4))//'_PAWAVES.nc'
289 1 : write(msg,'(2a)')' Opening file for AE PAW wave functions: ',trim(fname)
290 3 : call wrtout([std_out, ab_out], msg, 'PERS')
291 :
292 1 : if (xmpi_comm_rank(comm_cell) == master) then
293 1 : NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
294 :
295 1 : fform = 602
296 1 : NCF_CHECK(hdr%ncwrite(ncid, fform, nc_define=.True.))
297 :
298 : ! Define wavefunctions in real space on the dense FFT mesh
299 : ! Fortran layout:
300 : !real_space_wavefunctions: double 8d array with shape:
301 : ! [real_or_complex_wavefunctions]
302 : ! [number_of_grid_points_vector1][number_of_grid_points_vector2][number_of_grid_points_vector3]
303 : ! [number_of_spinor_components]
304 : ! [max_number_of_states][number_of_kpoints][number_of_spins]
305 :
306 : ncerr = nctk_def_dims(ncid, [ &
307 : nctkdim_t("real_or_complex_wavefunctions", 2), &
308 : nctkdim_t("number_of_grid_points_vector1", n1), &
309 : nctkdim_t("number_of_grid_points_vector2", n2), &
310 : nctkdim_t("number_of_grid_points_vector3", n3) &
311 5 : ], defmode=.True.)
312 1 : NCF_CHECK(ncerr)
313 :
314 : shape_str = "real_or_complex_wavefunctions, &
315 : & number_of_grid_points_vector1, number_of_grid_points_vector2, number_of_grid_points_vector3, &
316 : & number_of_spinor_components, &
317 1 : & max_number_of_states, number_of_kpoints, number_of_spins"
318 :
319 : ! Define wavefunctions in real space.
320 : ncerr = nctk_def_arrays(ncid, [&
321 : nctkarr_t('ur_ae', "dp", shape_str),&
322 : nctkarr_t('ur_pw', "dp", shape_str),&
323 : nctkarr_t('ur_ae_onsite', "dp", shape_str),&
324 : nctkarr_t('ur_ps_onsite', "dp", shape_str) &
325 5 : ], defmode=.True.)
326 1 : NCF_CHECK(ncerr)
327 :
328 : ! Complete the geometry information.
329 1 : NCF_CHECK(crystal%ncwrite(ncid))
330 1 : NCF_CHECK(ebands%ncwrite(ncid))
331 :
332 1 : NCF_CHECK(nf90_close(ncid))
333 : end if
334 :
335 1 : call xmpi_barrier(comm_cell)
336 :
337 : ! Reopen the file in parallel inside comm_cell
338 : ! Note that we use individual IO thus there's no need to handle idle processes
339 : ! if paral_kgb == 0 and nprocs > nkpt * nsppol
340 1 : NCF_CHECK(nctk_open_modify(ncid, fname, comm_cell))
341 1 : ae_ncid = nctk_idname(ncid, "ur_ae")
342 1 : pw_ncid = nctk_idname(ncid, "ur_pw")
343 1 : aeons_ncid = nctk_idname(ncid, "ur_ae_onsite")
344 1 : psons_ncid = nctk_idname(ncid, "ur_ps_onsite")
345 :
346 1 : NCF_CHECK(nctk_set_datamode(ncid))
347 : endif
348 :
349 : !Init structure storing phi_{nlm} and tphi_(nlm} on the dense FFT points located in the PAW spheres.
350 4 : ABI_MALLOC(Paw_onsite,(natom))
351 1 : if (paral_atom) then
352 : call paw_pwaves_lmn_init(Paw_onsite,my_natom,natom,crystal%ntypat,crystal%rprimd,crystal%xcart,&
353 0 : Pawtab,Pawrad,local_pawfgrtab, comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
354 : else
355 : call paw_pwaves_lmn_init(Paw_onsite,my_natom,natom,crystal%ntypat,crystal%rprimd,crystal%xcart,&
356 1 : Pawtab,Pawrad,local_pawfgrtab)
357 : end if
358 :
359 1 : bdtot_index=0; icg=0; ibg=0; norm_rerr=smallest_real
360 4 : ABI_MALLOC(paw_compl,(maxval(nband)))
361 7 : paw_compl=zero
362 :
363 : ! === Loop over spin ===
364 2 : do isppol=1,nsppol
365 1 : ikg=0; start_kpt=1; stop_kpt=nkpt
366 :
367 : ! Check if k-point was specified (only serial)
368 1 : if (present(set_k) .and. nprocs==1) then
369 1 : if (set_k/=0) then
370 1 : start_kpt = set_k
371 1 : stop_kpt = set_k
372 : !ABI_ERROR("set_k")
373 : end if
374 : end if
375 :
376 : ! === Loop over k points ===
377 3 : do ikpt=start_kpt,stop_kpt
378 4 : kpoint = kpt(:,ikpt)
379 1 : nband_k = nband(ikpt+(isppol-1)*nkpt)
380 1 : npw_k = npwarr(ikpt)
381 1 : istwf_k = istwfk(ikpt)
382 :
383 1 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) then
384 1 : bdtot_index=bdtot_index+nband_k
385 : !ABI_ERROR("cycle in seq!")
386 : cycle
387 : end if
388 :
389 21 : do i3=0,n3-1
390 20 : rfft(3)=DBLE(i3)/n3
391 381 : do i2=0,n2-1
392 360 : rfft(2)=DBLE(i2)/n2
393 6140 : do i1=0,n1-1
394 5760 : rfft(1)=DBLE(i1)/n1
395 5760 : ifft = 1 +i1 +i2*n1 +i3*n1*n2
396 23040 : phkr(1,ifft) = COS(two_pi*dot_product(kpoint,rfft))
397 23400 : phkr(2,ifft) = SIN(two_pi*dot_product(kpoint,rfft))
398 : end do
399 : end do
400 : end do
401 : ! phkr(1,:)=one; phkr(2,:)=zero
402 :
403 : ! Calculate the phase for the onsite PAW contributions.
404 2 : do iatom=1,my_natom
405 1 : nfgd=local_pawfgrtab(iatom)%nfgd ! no. of points in the fine grid for this PAW sphere
406 1431 : do ifgd=1,nfgd
407 5716 : arg = -two_pi* dot_product(r0shift(:,ifgd,iatom),kpoint)
408 1429 : phk_atm(1,ifgd,iatom) = COS(arg)
409 1430 : phk_atm(2,ifgd,iatom) = SIN(arg)
410 : end do
411 : end do
412 :
413 16 : ABI_MALLOC(Cprj_k,(natom,dtset%nspinor*nband_k))
414 1 : call pawcprj_alloc(Cprj_k,0,dimcprj)
415 :
416 : ! Extract cprj for this k-point.
417 1 : ibsp=0
418 7 : do iband=1,nband_k
419 13 : do ispinor=1,dtset%nspinor
420 6 : ibsp=ibsp+1
421 18 : do iatom=1,natom
422 246 : Cprj_k(iatom,ibsp)%cp(:,:)=Cprj(iatom,ibsp+ibg)%cp(:,:)
423 : end do
424 : end do
425 : end do
426 :
427 3 : ABI_MALLOC(kg_k,(3,npw_k))
428 :
429 : ! Extract G-vectors.
430 117 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
431 1 : call sphereboundary(gbound,istwf_k,kg_k,mgfftf,npw_k)
432 :
433 : ! If a single band is requested, neuter the loop (only serial)
434 1 : start_band = 1; stop_band = nband_k
435 1 : if (present(set_band).AND.nprocs==1) then
436 1 : if (set_band/=0) then
437 1 : start_band = set_band
438 1 : stop_band = set_band
439 : !ABI_ERROR("set_band")
440 : end if
441 : end if
442 :
443 : ! Loop over bands.
444 7 : do iband=start_band,stop_band
445 :
446 : ! Fourier transform on the real fft box of the smooth part.
447 6 : ndat=Dtset%nspinor
448 18 : ABI_MALLOC(fofgin,(2,npw_k*ndat))
449 30 : ABI_MALLOC(fofr,(2,n4,n5,n6*ndat))
450 :
451 180 : do ipw=1,npw_k*dtset%nspinor
452 528 : fofgin(:,ipw)=cg(:,ipw+(iband-1)*npw_k*dtset%nspinor+icg)
453 : end do
454 :
455 : ! Complex can be set to 0 with this option(0) of fourwf
456 6 : option=0; cplex=0; npwout=1
457 18 : ABI_MALLOC(denpot,(cplex*n4,n5,n6))
458 18 : ABI_MALLOC(fofgout,(2,npwout*ndat))
459 :
460 : call fourwf(cplex,denpot,fofgin(:,1:npw_k),fofgout,fofr(:,:,:,1:n6),gbound,gbound,istwf_k,kg_k,kg_k,&
461 : mgfftf,MPI_enreg,1,ngfftf,npw_k,npwout,n4,n5,n6,option,tim_fourwf0,weight1,weight1,&
462 6 : gpu_option=Dtset%gpu_option)
463 :
464 : ! Here I do not know if fourwf works in the case of spinors,
465 : ! It seems that not all fftalg option support ndata! should check!
466 : ! Do not forget to declare real(dp)::fofgin_down(:,:) to use the following statements
467 6 : if (Dtset%nspinor==2) then
468 0 : ABI_MALLOC(fofgin_down,(2,npw_k))
469 0 : ABI_MALLOC(fofr_down,(2,n4,n5,n6))
470 0 : fofgin_down(:,:)=fofgin(:,1+npw_k:2*npw_k)
471 : ! Complex can be set to 0 with this option(0) of fourwf
472 : ! cplex=1; option=1; npwout=1; ndat=1
473 : ! NOTE: fofr_down can NOT be replaced by fofr(:,:,:,n6+1:2*n6), or else
474 : ! the data in fofr(:,:,:,1:n6) will be the same with fofr(:,:,:,n6+1:2*n6)
475 : call fourwf(cplex,denpot,fofgin_down,fofgout,fofr_down,gbound,gbound,istwf_k,kg_k,kg_k,&
476 0 : mgfftf,MPI_enreg,1,ngfftf,npw_k,npwout,n4,n5,n6,option,tim_fourwf0,weight1,weight1)
477 0 : ABI_FREE(fofgin_down)
478 : end if
479 :
480 18 : ABI_MALLOC(ur_ae,(2,n1*n2*n3*ndat))
481 18 : ABI_MALLOC(ur_ae_onsite,(2,n1*n2*n3))
482 12 : ABI_MALLOC(ur_ps_onsite,(2,n1*n2*n3))
483 12 : ABI_MALLOC(ur_pw,(2,n1*n2*n3*ndat))
484 18 : ABI_MALLOC(ur_mask,(n1*n2*n3))
485 :
486 449286 : ur_ae=zero;ur_ae_onsite=zero;ur_ps_onsite=zero;ur_pw=zero;ur_mask=zero
487 :
488 : ! * Add phase e^{ikr} since it is contained in cprj.
489 126 : do i3=1,n3
490 2286 : do i2=1,n2
491 36840 : do i1=1,n1
492 34560 : ii = i1 + n1*(i2-1)+ n1*n2*(i3-1)
493 103680 : ur_pw(:,ii)=fofr(:,i1,i2,i3) ! Save pw part separately without the phase.
494 34560 : ur_ae(1,ii)= fofr(1,i1,i2,i3) * phkr(1,ii) - fofr(2,i1,i2,i3) * phkr(2,ii)
495 34560 : ur_ae(2,ii)= fofr(1,i1,i2,i3) * phkr(2,ii) + fofr(2,i1,i2,i3) * phkr(1,ii)
496 36720 : if(Dtset%nspinor==2) then
497 0 : ur_pw(:,ii+n1*n2*n3)=fofr_down(:,i1,i2,i3) ! Save pw part separately without the phase.
498 0 : ur_ae(1,ii+n1*n2*n3)= fofr_down(1,i1,i2,i3) * phkr(1,ii) - fofr_down(2,i1,i2,i3) * phkr(2,ii)
499 0 : ur_ae(2,ii+n1*n2*n3)= fofr_down(1,i1,i2,i3) * phkr(2,ii) + fofr_down(2,i1,i2,i3) * phkr(1,ii)
500 : end if
501 : end do
502 : end do
503 : end do
504 6 : ABI_FREE(fofr)
505 :
506 6 : if(Dtset%nspinor==2) then
507 0 : ABI_FREE(fofr_down)
508 : end if
509 :
510 : ! === Add onsite term on the augmented FFT mesh ===
511 12 : do iatom=1,my_natom
512 6 : itypat =local_pawfgrtab(iatom)%itypat
513 6 : lmn_size=Pawtab(itypat)%lmn_size
514 6 : ln_size =Pawtab(itypat)%basis_size ! no. of nl elements in PAW basis
515 6 : nfgd =local_pawfgrtab(iatom)%nfgd ! no. of points in the fine grid for this PAW sphere
516 :
517 6 : ibsp=(iband-1)*dtset%nspinor
518 18 : do ispinor=1,dtset%nspinor
519 6 : ibsp=ibsp+1
520 90 : do jlmn=1,lmn_size
521 78 : jl=Pawtab(itypat)%indlmn(1,jlmn)
522 78 : jm=Pawtab(itypat)%indlmn(2,jlmn)
523 78 : cp_fact(1) = Cprj_k(iatom,ibsp)%cp(1,jlmn) *sqrt(crystal%ucvol) ! Magic factor
524 78 : cp_fact(2) = Cprj_k(iatom,ibsp)%cp(2,jlmn) *sqrt(crystal%ucvol)
525 :
526 111546 : do ifgd=1,nfgd ! loop over fine grid points in current PAW sphere.
527 111462 : ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid
528 111462 : phj = Paw_onsite(iatom)% phi(ifgd,jlmn)
529 111462 : tphj = Paw_onsite(iatom)%tphi(ifgd,jlmn)
530 : ! old code
531 : !re_p = cp_fact(1); im_p = cp_fact(2)
532 : ! apply the phase
533 111462 : re_p = cp_fact(1) * phk_atm(1,ifgd,iatom) - cp_fact(2) * phk_atm(2,ifgd,iatom)
534 111462 : im_p = cp_fact(1) * phk_atm(2,ifgd,iatom) + cp_fact(2) * phk_atm(1,ifgd,iatom)
535 :
536 111462 : ur_ae(1,ifftsph+(ispinor-1)*nfftot) = ur_ae(1,ifftsph+(ispinor-1)*nfftot) + re_p * (phj-tphj)
537 111462 : ur_ae(2,ifftsph+(ispinor-1)*nfftot) = ur_ae(2,ifftsph+(ispinor-1)*nfftot) + im_p * (phj-tphj)
538 111462 : ur_ae_onsite(1,ifftsph) = ur_ae_onsite(1,ifftsph) + re_p * phj
539 111462 : ur_ae_onsite(2,ifftsph) = ur_ae_onsite(2,ifftsph) + im_p * phj
540 111462 : ur_ps_onsite(1,ifftsph) = ur_ps_onsite(1,ifftsph) + re_p * tphj
541 111462 : ur_ps_onsite(2,ifftsph) = ur_ps_onsite(2,ifftsph) + im_p * tphj
542 111540 : ur_mask(ifftsph) = one
543 : end do
544 :
545 : end do !jlmn
546 : end do !ispinor
547 : end do !iatom
548 :
549 6 : if (paral_atom) then
550 0 : ABI_MALLOC(buf_tmp,(2,n1*n2*n3,3))
551 0 : buf_tmp(:,:,1) = ur_ae
552 0 : buf_tmp(:,:,2) = ur_ae_onsite
553 0 : buf_tmp(:,:,3) = ur_ps_onsite
554 0 : call xmpi_sum(buf_tmp,my_comm_atom,ierr)
555 0 : ur_ae = buf_tmp(:,:,1)
556 0 : ur_ae_onsite= buf_tmp(:,:,2)
557 0 : ur_ps_onsite= buf_tmp(:,:,3)
558 0 : ABI_FREE(buf_tmp)
559 : end if
560 :
561 : ! * Remove the phase e^{ikr}, we store u(r).
562 126 : do i3=1,n3
563 2286 : do i2=1,n2
564 36840 : do i1=1,n1
565 34560 : ii = i1 + n1*(i2-1)+ n1*n2*(i3-1)
566 34560 : reur=ur_ae(1,ii)
567 34560 : imur=ur_ae(2,ii)
568 34560 : ur_ae(1,ii)= reur * phkr(1,ii) + imur * phkr(2,ii)
569 34560 : ur_ae(2,ii)= -reur * phkr(2,ii) + imur * phkr(1,ii)
570 34560 : if(Dtset%nspinor==2) then
571 0 : reur=ur_ae(1,ii+nfftot) ! Important!
572 0 : imur=ur_ae(2,ii+nfftot)
573 0 : ur_ae(1,ii+nfftot)= reur * phkr(1,ii) + imur * phkr(2,ii)
574 0 : ur_ae(2,ii+nfftot)= -reur * phkr(2,ii) + imur * phkr(1,ii)
575 : end if
576 34560 : reur=ur_ae_onsite(1,ii)
577 34560 : imur=ur_ae_onsite(2,ii)
578 34560 : ur_ae_onsite(1,ii)= reur * phkr(1,ii) + imur * phkr(2,ii)
579 34560 : ur_ae_onsite(2,ii)= -reur * phkr(2,ii) + imur * phkr(1,ii)
580 34560 : reur=ur_ps_onsite(1,ii)
581 34560 : imur=ur_ps_onsite(2,ii)
582 34560 : ur_ps_onsite(1,ii)= reur * phkr(1,ii) + imur * phkr(2,ii)
583 36720 : ur_ps_onsite(2,ii)= -reur * phkr(2,ii) + imur * phkr(1,ii)
584 : end do
585 : end do
586 : end do
587 :
588 6 : norm=zero
589 180 : do ii=1,npw_k*Dtset%nspinor
590 180 : norm=norm+fofgin(1,ii)**2+fofgin(2,ii)**2
591 : end do
592 6 : write(std_out,'(a,2i5,f22.16)',advance='no') 'ikpt,iband, norm (G,PSWF)=',ikpt,iband,norm
593 6 : norm=zero
594 34566 : do ifft=1,nfftot*Dtset%nspinor
595 34566 : norm = norm + ur_ae(1,ifft)**2+ur_ae(2,ifft)**2
596 : end do
597 6 : norm=norm/nfftot
598 6 : norm_rerr = MAX((ABS(norm-one))*100,norm_rerr)
599 6 : write(std_out,*)"norm (R,AEWF)= ",norm
600 6 : call flush_unit(std_out)
601 :
602 6 : norm=zero
603 34566 : do ifft=1, nfftot
604 34566 : norm=norm+(ur_pw(1,ifft)**2+ur_pw(2,ifft)**2)*ur_mask(ifft)
605 : enddo
606 34566 : do ifft=1, nfftot
607 : paw_compl(iband)=paw_compl(iband)+(ur_pw(1,ifft)*ur_ps_onsite(1,ifft)+ur_pw(2,ifft)*ur_ps_onsite(2,ifft))&
608 34566 : & /norm/nsppol*dtset%wtk(ikpt) !/(stop_kpt-start_kpt+1)
609 : enddo
610 :
611 : ! MS: Various testing and debugging options
612 6 : if (.TRUE..and.nprocs==1) then
613 6 : if (present(pseudo_norms)) then
614 : ! Check the supposedly zero overlap |\tilde{Psi_n}-\tilde{Psi_n^1}|^2
615 12 : ABI_MALLOC(dummy_1d,(n1*n2*n3))
616 34566 : dummy_1d = zero
617 6 : norm = zero
618 34566 : do ifft = 1, nfftot
619 : dummy_1d(ifft) = ((ur_pw(1,ifft)-ur_ps_onsite(1,ifft))**2 &
620 34560 : + (ur_pw(2,ifft)-ur_ps_onsite(2,ifft))**2) * ur_mask(ifft)
621 34566 : norm = norm + dummy_1d(ifft)
622 : end do
623 6 : norm = norm / nfftot
624 6 : pseudo_norms(isppol,ikpt,iband) = norm
625 6 : ABI_FREE(dummy_1d)
626 : end if
627 :
628 : else
629 : write(msg,'(5a)')&
630 0 : "The option to print PAW all-electron wavefunctions is on, but execution ",ch10,&
631 0 : "is in parallel on two or more processors. XcrysDen files with individual con-",ch10,&
632 0 : "tributions will not be written. In order to enable this you must run in serial."
633 0 : ABI_WARNING(msg)
634 : end if ! Check if serial run
635 :
636 6 : if(dtset%pawprtwf==1) then
637 : ncerr = nf90_put_var(ncid, ae_ncid, ur_ae, &
638 102 : start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1])
639 6 : NCF_CHECK(ncerr)
640 :
641 : ncerr = nf90_put_var(ncid, pw_ncid, ur_pw, &
642 102 : start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1])
643 6 : NCF_CHECK(ncerr)
644 :
645 : ncerr = nf90_put_var(ncid, aeons_ncid, ur_ae_onsite, &
646 102 : start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1])
647 6 : NCF_CHECK(ncerr)
648 :
649 : ncerr = nf90_put_var(ncid, psons_ncid, ur_ps_onsite, &
650 102 : start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1])
651 6 : NCF_CHECK(ncerr)
652 : endif
653 :
654 6 : ABI_FREE(ur_ae)
655 6 : ABI_FREE(ur_ae_onsite)
656 6 : ABI_FREE(ur_ps_onsite)
657 6 : ABI_FREE(ur_pw)
658 6 : ABI_FREE(ur_mask)
659 6 : ABI_FREE(fofgin)
660 6 : ABI_FREE(fofgout)
661 7 : ABI_FREE(denpot)
662 : end do !nband_k
663 :
664 1 : bdtot_index=bdtot_index+nband_k
665 :
666 1 : if (mkmem/=0) then
667 1 : ibg=ibg+dtset%nspinor*nband_k
668 1 : icg=icg+npw_k*dtset%nspinor*nband_k
669 1 : ikg=ikg+npw_k
670 : end if
671 :
672 1 : ABI_FREE(kg_k)
673 :
674 1 : call pawcprj_free(Cprj_k)
675 8 : ABI_FREE(Cprj_k)
676 :
677 : end do !ikpt
678 : end do !nsppol
679 :
680 1 : ABI_FREE(phkr)
681 1 : ABI_FREE(gbound)
682 :
683 : ! Free augmentation waves.
684 1 : call paw_pwaves_lmn_free(Paw_onsite)
685 2 : ABI_FREE(Paw_onsite)
686 :
687 : ! Maximum relative error over CPUs.
688 1 : call xmpi_max(norm_rerr,max_rerr,comm_cell,ierr)
689 1 : write(std_out,*)"max_rerr=",max_rerr
690 :
691 1 : if (max_rerr > ten) then
692 : write(msg,'(7a)')&
693 1 : "Inaccuracy on the normalization of the wave funtions exceeds 10%. ",ch10,&
694 1 : "Likely due to the use of a too coarse FFT mesh or unconverged wavefunctions. ",ch10,&
695 1 : "Numerical values inside the augmentation regions might be inaccurate. ",ch10,&
696 2 : "Action: increase pawecutdg in your input file. "
697 1 : ABI_COMMENT(msg)
698 : end if
699 :
700 :
701 1 : open(file=trim(dtfil%filnam_ds(4))//'_PAWCOMPL', newunit=tmp_unt,status='unknown',form='formatted')
702 8 : do iband=1,maxval(nband)
703 7 : write(tmp_unt,*) iband,paw_compl(iband)
704 : enddo
705 1 : close(tmp_unt)
706 :
707 1 : ABI_FREE(paw_compl)
708 1 : ABI_FREE(r0shift)
709 1 : ABI_FREE(phk_atm)
710 1 : call pawfgrtab_free(local_pawfgrtab)
711 :
712 : ! Destroy atom table used for parallelism
713 1 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
714 :
715 : DBG_EXIT("COLL")
716 :
717 5 : end subroutine pawmkaewf
718 : !!***
719 :
720 : end module m_paw_mkaewf
721 : !!***
|
