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