Branch data Line data Source code
1 : : !!****m* ABINIT/m_extfpmd
2 : : !! NAME
3 : : !! m_extfpmd
4 : : !!
5 : : !! FUNCTION
6 : : !! This module provides routines to run computations at very high temperature
7 : : !! with reduced number of bands. High energy orbitals are represented as
8 : : !! pure plane waves.
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2018-2026 ABINIT group (A. Blanchet)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !!
16 : : !! TODO
17 : : !! 1) Add contribution to conductivity.
18 : : !! 2) Smooth the contributions.
19 : : !!
20 : : !! SOURCE
21 : :
22 : : #if defined HAVE_CONFIG_H
23 : : #include "config.h"
24 : : #endif
25 : :
26 : : #include "abi_common.h"
27 : :
28 : : module m_extfpmd
29 : : use defs_basis
30 : : use defs_abitypes
31 : : use m_io_tools
32 : : use m_errors
33 : : use m_geometry
34 : : use m_special_funcs
35 : : use m_specialmsg
36 : : use m_xmpi
37 : : use m_dtset, only : dataset_type
38 : : use m_energies, only : energies_type
39 : : use m_gsphere, only : getkpgnorm
40 : : use m_kg, only : mkkin,kpgio
41 : : use m_mpinfo, only : ptabs_fourdp,proc_distrb_cycle,copy_mpi_enreg,destroy_mpi_enreg
42 : : use m_numeric_tools, only : simpson,simpson_int
43 : : use m_spacepar, only : meanvalue_g
44 : :
45 : : implicit none
46 : : public :: extfpmd_dos,extfpmd_e_fg,extfpmd_i_fg,extfpmd_chkinp
47 : : !!***
48 : :
49 : : !----------------------------------------------------------------------
50 : :
51 : : !!****t* m_extfpmd/extfpmd_type
52 : : !! NAME
53 : : !! extfpmd_type
54 : : !!
55 : : !! FUNCTION
56 : : !! Store extfpmd functions and parameters.
57 : : !!
58 : : !! SOURCE
59 : : type,public :: extfpmd_type
60 : : logical :: truecg
61 : : integer :: bcut,mband,nbcut,nbdbuf,nfftf,nspden,version,pawsph
62 : : real(dp) :: ebcut,edc_kinetic,e_kinetic,entropy
63 : : real(dp) :: nelect,eshift,ucvol,el_temp,bandshift,eshift_paw
64 : : real(dp) :: nelect_res, nelect_respc
65 : : real(dp),allocatable :: vtrial(:,:)
66 : : real(dp),allocatable :: nelectarr(:,:)
67 : : real(dp),allocatable :: bandshiftk(:)
68 : : type(MPI_type) :: mpi_enreg
69 : : contains
70 : : procedure :: compute_e_kinetic
71 : : procedure :: compute_entropy
72 : : procedure :: compute_nelect
73 : : procedure :: compute_eshift
74 : : procedure :: init
75 : : procedure :: destroy
76 : : end type extfpmd_type
77 : : !!***
78 : :
79 : : contains
80 : :
81 : : !!****f* ABINIT/m_extfpmd/init
82 : : !! NAME
83 : : !! init
84 : : !!
85 : : !! FUNCTION
86 : : !! Initialize extfpmd_type object, memory allocation of arrays...
87 : : !!
88 : : !! INPUTS
89 : : !! this=extfpmd_type object concerned
90 : : !! mband=maximum number of bands
91 : : !! extfpmd_eshift=pre-defined extfpmd energy shift
92 : : !! nbcut=number of states used to average the constant potential value
93 : : !! nbdbuf=Number of bands in the buffer to converge scf cycle with extfpmd models
94 : : !! nfftf=number of FFT fine grid points
95 : : !! nspden=number of spin-density components
96 : : !! nsppol=number of independent spin WF components
97 : : !! nkpt=number of k-points
98 : : !! occopt=option for occupancies
99 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
100 : : !! tphysel="physical" electronic temperature with FD occupations
101 : : !! tsmear=smearing energy or temperature (if metal)
102 : : !! version=extfpmd implementation version
103 : : !! mpi_enreg=information about MPI parallelization
104 : : !! extfpmd_mband=number of extfpmd bands
105 : : !!
106 : : !! OUTPUT
107 : : !! this=extfpmd_type object concerned
108 : : !!
109 : : !! SOURCE
110 : 7 : subroutine init(this,mband,extfpmd_eshift,nbcut,nbdbuf,nfftf,nspden,&
111 : : & nsppol,nkpt,occopt,rprimd,tphysel,tsmear,version,mpi_enreg,extfpmd_mband,pawsph)
112 : : ! Arguments -------------------------------
113 : : ! Scalars
114 : : class(extfpmd_type),intent(inout) :: this
115 : : integer,intent(in) :: mband,nbcut,nbdbuf,nfftf,nspden,pawsph
116 : : integer,intent(in) :: nsppol,nkpt,version,extfpmd_mband,occopt
117 : : real(dp),intent(in) :: extfpmd_eshift,tphysel,tsmear
118 : : type(MPI_type),intent(in) :: mpi_enreg
119 : : ! Arrays
120 : : real(dp),intent(in) :: rprimd(3,3)
121 : :
122 : : ! Local variables -------------------------
123 : : ! Arrays
124 : : real(dp) :: gprimd(3,3),rmet(3,3),gmet(3,3)
125 : :
126 : : ! *********************************************************************
127 : :
128 : 7 : this%bcut=mband-nbdbuf
129 : 7 : this%nbcut=nbcut
130 : 7 : this%mband=extfpmd_mband
131 : 7 : this%nbdbuf=nbdbuf
132 : 7 : this%version=version
133 [ - + + - : 28 : ABI_MALLOC(this%vtrial,(nfftf,nspden))
- + + - +
- + - - +
- + - + ]
134 [ + + + + ]: 117902 : this%vtrial(:,:)=zero
135 : 7 : this%nfftf=nfftf
136 : 7 : this%nspden=nspden
137 : 7 : this%ebcut=zero
138 : 7 : this%edc_kinetic=zero
139 : 7 : this%e_kinetic=zero
140 : 7 : this%entropy=zero
141 : 7 : this%nelect=zero
142 : 7 : this%nelect_res=zero
143 : 7 : this%nelect_respc=zero
144 : 7 : this%bandshift=zero
145 [ - + + - : 21 : ABI_MALLOC(this%bandshiftk,(nkpt*nsppol))
+ - + - -
+ - + -
+ ]
146 [ + + ]: 20 : this%bandshiftk(:)=zero
147 : 7 : this%eshift=extfpmd_eshift
148 : 7 : this%eshift_paw=zero
149 : 7 : call metric(gmet,gprimd,-1,rmet,rprimd,this%ucvol)
150 [ - + - - : 7 : this%el_temp=merge(tphysel,tsmear,tphysel>tol8.and.occopt/=3.and.occopt/=9)
- - ]
151 : 7 : this%pawsph=pawsph
152 : :
153 [ + + ]: 7 : if(this%version==5) then
154 : : ! Make a copy of mpi_enreg in order to cycle.
155 : 1 : call copy_mpi_enreg(mpi_enreg,this%mpi_enreg)
156 : : end if
157 : :
158 : 7 : end subroutine init
159 : : !!***
160 : :
161 : : !!****f* ABINIT/m_extfpmd/destroy
162 : : !! NAME
163 : : !! destroy
164 : : !!
165 : : !! FUNCTION
166 : : !! Destroy extfpmd_type object, memory deallocation of arrays...
167 : : !!
168 : : !! INPUTS
169 : : !! this=extfpmd_type object concerned
170 : : !!
171 : : !! OUTPUT
172 : : !! this=extfpmd_type object concerned
173 : : !!
174 : : !! SOURCE
175 : 7 : subroutine destroy(this)
176 : :
177 : : ! Arguments -------------------------------
178 : : ! Scalars
179 : : class(extfpmd_type),intent(inout) :: this
180 : :
181 : : ! *********************************************************************
182 : :
183 [ + + ]: 7 : if(this%version==5) then
184 : 1 : call destroy_mpi_enreg(this%mpi_enreg)
185 : : end if
186 : :
187 [ + + + + ]: 117902 : this%vtrial(:,:)=zero
188 [ - + ]: 7 : ABI_FREE(this%vtrial)
189 [ - + ]: 7 : if(allocated(this%nelectarr)) then
190 [ # # # # ]: 0 : this%nelectarr(:,:)=zero
191 : 0 : ABI_FREE(this%nelectarr)
192 : : end if
193 [ + + ]: 20 : this%bandshiftk(:)=zero
194 [ - + ]: 7 : ABI_FREE(this%bandshiftk)
195 : 7 : this%nfftf=0
196 : 7 : this%nspden=0
197 : 7 : this%bcut=0
198 : 7 : this%mband=0
199 : 7 : this%nbcut=0
200 : 7 : this%nbdbuf=0
201 : 7 : this%version=1
202 : 7 : this%ebcut=zero
203 : 7 : this%edc_kinetic=zero
204 : 7 : this%e_kinetic=zero
205 : 7 : this%entropy=zero
206 : 7 : this%nelect=zero
207 : 7 : this%nelect_res=zero
208 : 7 : this%nelect_respc=zero
209 : 7 : this%bandshift=zero
210 : 7 : this%eshift=zero
211 : 7 : this%eshift_paw=zero
212 : 7 : this%ucvol=zero
213 : 7 : this%el_temp=zero
214 : 7 : this%pawsph=0
215 : 7 : end subroutine destroy
216 : : !!***
217 : :
218 : : !!****f* ABINIT/m_extfpmd/compute_eshift
219 : : !! NAME
220 : : !! compute_eshift
221 : : !!
222 : : !! FUNCTION
223 : : !! Computes the energy shift factor $U_0$ corresponding to constant
224 : : !! potential contribution.
225 : : !!
226 : : !! INPUTS
227 : : !! this=extfpmd_type object concerned
228 : : !! eigen(mband*nkpt*nsppol)=eigenvalues (hartree)
229 : : !! eknk(mband*nkpt*nsppol)=kinetic energies (hartree)
230 : : !! mband=maximum number of bands
231 : : !! nband(nkpt*nsppol)=desired number of bands at each k point
232 : : !! nfftf=number of FFT fine grid points
233 : : !! nkpt=number of k points
234 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
235 : : !! nspden=number of spin-density components
236 : : !! wtk(nkpt)=k point weights
237 : : !! vtrial(nfftf,nspden)=GS potential on the fine grid (Hartree)
238 : : !!
239 : : !! OUTPUT
240 : : !! this=extfpmd_type object concerned
241 : : !!
242 : : !! SOURCE
243 : 90 : subroutine compute_eshift(this,eigen,eknk,mband,nband,nfftf,nkpt,nsppol,nspden,wtk,vtrial)
244 : : ! Arguments -------------------------------
245 : : ! Scalars
246 : : class(extfpmd_type),intent(inout) :: this
247 : : integer,intent(in) :: mband,nfftf,nkpt,nsppol,nspden
248 : : ! Arrays
249 : : integer,intent(in) :: nband(nkpt*nsppol)
250 : : real(dp),intent(in) :: eigen(mband*nkpt*nsppol)
251 : : real(dp),intent(in) :: eknk(mband*nkpt*nsppol)
252 : : real(dp),intent(in) :: wtk(nkpt)
253 : : real(dp),intent(in) :: vtrial(nfftf,nspden)
254 : :
255 : : ! Local variables -------------------------
256 : : ! Scalars
257 : : integer :: band_index,ii,ikpt,isppol,nband_k
258 : :
259 : : ! *********************************************************************
260 [ + - + - : 993422 : this%vtrial=vtrial
- + - - +
+ + + ]
261 : :
262 [ + + ]: 90 : if(this%version==2) then
263 : : ! Computes U_0^{HEG} from the difference between
264 : : ! eigenvalues and Fermi gas energies, averaged
265 : : ! over lasts nbcut bands.
266 : 11 : this%eshift=zero
267 : 11 : band_index=0
268 [ + + ]: 22 : do isppol=1,nsppol
269 [ + + ]: 44 : do ikpt=1,nkpt
270 : 22 : nband_k=nband(ikpt+(isppol-1)*nkpt)
271 [ + + ]: 242 : do ii=nband_k-this%nbdbuf-this%nbcut+1,nband_k-this%nbdbuf
272 : : this%eshift=this%eshift+&
273 : 242 : & wtk(ikpt)*(eigen(band_index+ii)-extfpmd_e_fg(dble(ii),this%ucvol))
274 : : end do
275 : 33 : band_index=band_index+nband_k
276 : : end do
277 : : end do
278 : 11 : this%eshift=this%eshift/this%nbcut
279 [ + + ]: 79 : else if(this%version==3) then
280 : : ! Computes U_0^K from the difference between
281 : : ! eigenvalues and kinetic energies, averaged
282 : : ! over lasts nbcut bands.
283 : 11 : this%eshift=zero
284 : 11 : band_index=0
285 [ + + ]: 22 : do isppol=1,nsppol
286 [ + + ]: 44 : do ikpt=1,nkpt
287 : 22 : nband_k=nband(ikpt+(isppol-1)*nkpt)
288 [ + + ]: 242 : do ii=nband_k-this%nbdbuf-this%nbcut+1,nband_k-this%nbdbuf
289 : : this%eshift=this%eshift+&
290 : 242 : & wtk(ikpt)*(eigen(band_index+ii)-eknk(band_index+ii))
291 : : end do
292 : 33 : band_index=band_index+nband_k
293 : : end do
294 : : end do
295 : 11 : this%eshift=this%eshift/this%nbcut
296 : : else
297 : : ! Computes U_0 from the sum of local
298 : : ! potentials (vtrial), averaging over all space.
299 : : ! Simplest and most precise way to evaluate U_0.
300 [ + + + + ]: 903176 : this%eshift=sum(this%vtrial)/(nfftf*nspden)
301 : : end if
302 : :
303 [ - + ]: 90 : if(this%pawsph==2) then
304 : 0 : this%eshift=this%eshift+this%eshift_paw
305 : : endif
306 : :
307 : : ! Get extended FPMD band energy cutoff
308 : 90 : this%ebcut=zero
309 : 90 : this%bandshift=zero
310 [ + + ]: 263 : this%bandshiftk(:)=zero
311 : : band_index=0
312 [ + + ]: 180 : do isppol=1,nsppol
313 [ + + ]: 353 : do ikpt=1,nkpt
314 : 173 : nband_k=nband(ikpt+(isppol-1)*nkpt)
315 : 173 : this%ebcut=this%ebcut+wtk(ikpt)*eigen(band_index+nband_k-this%nbdbuf)/nsppol
316 : : this%bandshift=this%bandshift+wtk(ikpt)*&
317 : 173 : & (extfpmd_i_fg(eigen(band_index+nband_k-this%nbdbuf)-this%eshift,this%ucvol)-(nband_k-this%nbdbuf))/nsppol
318 : 173 : this%bandshiftk(ikpt+(isppol-1)*nkpt)=extfpmd_i_fg(eigen(band_index+nband_k-this%nbdbuf)-this%eshift,this%ucvol)-(nband_k-this%nbdbuf)
319 : 263 : band_index=band_index+nband_k
320 : : end do
321 : : end do
322 : 90 : end subroutine compute_eshift
323 : : !!***eigen(band_index+nband_k-this%nbdbuf)
324 : :
325 : : !!****f* ABINIT/m_extfpmd/compute_nelect
326 : : !! NAME
327 : : !! compute_nelect
328 : : !!
329 : : !! FUNCTION
330 : : !! Computes the value of the integral corresponding to the missing
331 : : !! free electrons contribution after band cut, with an order 1/2
332 : : !! incomplete Fermi-Dirac integral.
333 : : !!
334 : : !! INPUTS
335 : : !! this=extfpmd_type object concerned
336 : : !! fermie=chemical potential (Hartree)
337 : : !! nband(nkpt*nsppol)=desired number of bands at each k point
338 : : !! nelect=number of electrons per unit cell
339 : : !! nkpt=number of k points
340 : : !! nspinor=number of spinor components
341 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
342 : : !! wtk(nkpt)=k point weights
343 : : !!
344 : : !! OUTPUT
345 : : !! this=extfpmd_type object concerned
346 : : !! nelect=number of electrons per unit cell
347 : : !!
348 : : !! SOURCE
349 : 4731 : subroutine compute_nelect(this,fermie,nband,nelect,nkpt,nspinor,nsppol,wtk)
350 : : ! Arguments -------------------------------
351 : : ! Scalars
352 : : integer,intent(in) :: nkpt,nsppol,nspinor
353 : : real(dp),intent(in) :: fermie
354 : : real(dp),intent(inout) :: nelect
355 : : class(extfpmd_type),intent(inout) :: this
356 : : ! Arrays
357 : : integer,intent(in) :: nband(nkpt*nsppol)
358 : : real(dp),intent(in) :: wtk(nkpt)
359 : :
360 : : ! Local variables -------------------------
361 : : ! Scalars
362 : : integer :: ifft,ispden,isppol,ikpt,iband,nband_k,ierr
363 : : real(dp) :: factor,gamma,xcut,fn,maxocc,nelect_tmp
364 : : ! Arrays
365 : 4731 : real(dp),allocatable :: gamma_hybrid_tf(:,:)
366 : 4731 : real(dp),allocatable :: xcut_hybrid_tf(:,:)
367 : :
368 : : ! *********************************************************************
369 : :
370 : 4731 : maxocc=two/(nsppol*nspinor)
371 : 4731 : factor=dsqrt(two)/(PI*PI)*this%ucvol*this%el_temp**(1.5)
372 : 4731 : gamma=(fermie-this%eshift)/this%el_temp
373 : 4731 : nelect_tmp=zero
374 : :
375 : : ! Computes extfpmd contribution to nelect integrating
376 : : ! over accessible states from bcut to infinity with
377 : : ! order 1/2 incomplete Fermi-Dirac integral.
378 [ + + ]: 4731 : if(this%version==2.or.this%version==4) then
379 : 1151 : xcut=extfpmd_e_fg(one*this%bcut+this%bandshift,this%ucvol)/this%el_temp
380 [ - + ]: 1151 : if(one*this%bcut+this%bandshift.lt.zero) xcut=zero
381 : 1151 : nelect=nelect+factor*djp12(xcut,gamma)
382 : : end if
383 : :
384 : : ! Computes extfpmd contribution to nelect integrating
385 : : ! over energy from ebcut to infinity with order 1/2
386 : : ! incomplete Fermi-Dirac integral.
387 [ + + ]: 4731 : if(this%version==1.or.this%version==3) then
388 : 2998 : xcut=(this%ebcut-this%eshift)/this%el_temp
389 [ - + ]: 2998 : if(this%ebcut.lt.this%eshift) xcut=zero
390 : 2998 : nelect=nelect+factor*djp12(xcut,gamma)
391 : : end if
392 : :
393 : : ! Computes extfpmd contribution to nelect summing
394 : : ! over accessible states from bcut to mband, with
395 : : ! integer band numbers. Total number of bands
396 : : ! is controlled with the input variable extfpmd_nband.
397 [ + + ]: 4731 : if(this%version==5) then
398 [ + + ]: 1164 : do isppol=1,nsppol
399 [ + + ]: 2328 : do ikpt=1,nkpt
400 : 1164 : nband_k=nband(ikpt+(isppol-1)*nkpt)
401 [ + - ]: 1164 : if(proc_distrb_cycle(this%mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,this%mpi_enreg%me_kpt)) cycle
402 [ + + ]: 9479034 : do iband=nband_k-this%nbdbuf+1,this%mband
403 : 9477288 : fn=fermi_dirac(extfpmd_e_fg(one*iband+this%bandshiftk(ikpt+(isppol-1)*nkpt),this%ucvol)+this%eshift,fermie,this%el_temp)
404 : 9478452 : nelect_tmp=nelect_tmp+wtk(ikpt)*maxocc*fn
405 : : end do
406 : : end do
407 : : end do
408 : 582 : call xmpi_sum(nelect_tmp,this%mpi_enreg%comm_kpt,ierr)
409 : 582 : nelect=nelect+nelect_tmp
410 : : end if
411 : :
412 : : ! Computes extfpmd contribution to nelect using a sum
413 : : ! of Fermi gas contributions for each point of the fftf grid.
414 : : ! Warning: This is not yet operational. Work in progress.
415 [ - + ]: 4731 : if(this%version==10) then
416 [ # # # # : 0 : ABI_MALLOC(gamma_hybrid_tf,(this%nfftf,this%nspden))
# # # # #
# # # # #
# # ]
417 [ # # # # : 0 : ABI_MALLOC(xcut_hybrid_tf,(this%nfftf,this%nspden))
# # # # #
# # # ]
418 [ # # ]: 0 : if(.not.allocated(this%nelectarr)) then
419 [ # # # # : 0 : ABI_MALLOC(this%nelectarr,(this%nfftf,this%nspden))
# # # # #
# # # #
# ]
420 [ # # # # ]: 0 : this%nelectarr(:,:)=zero
421 : : end if
422 [ # # # # ]: 0 : gamma_hybrid_tf(:,:)=(fermie-this%vtrial(:,:))/this%el_temp
423 [ # # # # ]: 0 : xcut_hybrid_tf(:,:)=(this%ebcut-this%vtrial(:,:))/this%el_temp
424 [ # # # # : 0 : if(ANY(this%ebcut.lt.this%vtrial(:,:))) xcut_hybrid_tf(:,:)=zero
# # # # #
# # # ]
425 : :
426 : : !$OMP PARALLEL DO
427 [ # # ]: 0 : do ifft=1,this%nfftf
428 [ # # ]: 0 : do ispden=1,this%nspden
429 : 0 : this%nelectarr(ifft,ispden)=factor*djp12(xcut_hybrid_tf(ifft,ispden),gamma_hybrid_tf(ifft,ispden))
430 : : end do
431 : : end do
432 : : !$OMP END PARALLEL DO
433 : :
434 [ # # # # ]: 0 : nelect=nelect+sum(this%nelectarr)/(this%nfftf*this%nspden)
435 [ # # # # ]: 0 : gamma_hybrid_tf(:,:)=zero
436 [ # # # # ]: 0 : xcut_hybrid_tf(:,:)=zero
437 [ # # ]: 0 : ABI_FREE(gamma_hybrid_tf)
438 [ # # ]: 0 : ABI_FREE(xcut_hybrid_tf)
439 : : end if
440 : 4731 : end subroutine compute_nelect
441 : : !!***
442 : :
443 : : !!****f* ABINIT/m_extfpmd/compute_e_kinetic
444 : : !! NAME
445 : : !! compute_e_kinetic
446 : : !!
447 : : !! FUNCTION
448 : : !! Computes the value of the integral corresponding to the missing
449 : : !! kinetic energy contribution of free electrons after band cut,
450 : : !! with an order 3/2 incomplete Fermi-Dirac integral.
451 : : !!
452 : : !! INPUTS
453 : : !! this=extfpmd_type object concerned
454 : : !! fermie=chemical potential (Hartree)
455 : : !! nkpt=number of k points
456 : : !! nspinor=number of spinor components
457 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
458 : : !! nband(nkpt*nsppol)=desired number of bands at each k point
459 : : !! wtk(nkpt)=k point weights
460 : : !!
461 : : !! OUTPUT
462 : : !! this=extfpmd_type object concerned
463 : : !!
464 : : !! SOURCE
465 : 90 : subroutine compute_e_kinetic(this,fermie,nkpt,nspinor,nsppol,nband,wtk)
466 : : ! Arguments -------------------------------
467 : : ! Scalars
468 : : integer,intent(in) :: nkpt,nspinor,nsppol
469 : : class(extfpmd_type),intent(inout) :: this
470 : : real(dp),intent(in) :: fermie
471 : : ! Arrays
472 : : integer,intent(in) :: nband(nkpt*nsppol)
473 : : real(dp),intent(in) :: wtk(nkpt)
474 : :
475 : : ! Local variables -------------------------
476 : : ! Scalars
477 : : logical :: cut_warn=.false.
478 : : integer :: ikpt,isppol,nband_k,iband,ierr,ifft,ispden
479 : : real(dp) :: factor,gamma,xcut,dotr
480 : : real(dp) :: e_kinetic_hybrid_tf,maxocc,fn
481 : : character(len=500) :: msg
482 : : ! Arrays
483 : 90 : real(dp),allocatable :: gamma_hybrid_tf(:,:),xcut_hybrid_tf(:,:)
484 : :
485 : : ! *********************************************************************
486 : :
487 : 90 : dotr=zero
488 : 90 : maxocc=two/(nsppol*nspinor)
489 : 90 : this%e_kinetic=zero
490 : 90 : factor=dsqrt(two)/(PI*PI)*this%ucvol*this%el_temp**(2.5)
491 : 90 : gamma=(fermie-this%eshift)/this%el_temp
492 : :
493 : : ! Computes extfpmd contribution to kinetic energy integrating
494 : : ! over accessible states from bcut to infinity with
495 : : ! order 3/2 incomplete Fermi-Dirac integral.
496 [ + + ]: 90 : if(this%version==2.or.this%version==4) then
497 : 22 : xcut=extfpmd_e_fg(one*this%bcut+this%bandshift,this%ucvol)/this%el_temp
498 [ - + ]: 22 : if(one*this%bcut+this%bandshift.lt.zero) then
499 : 0 : cut_warn=.true.
500 : 0 : xcut=zero
501 : : end if
502 : 22 : this%e_kinetic=this%e_kinetic+factor*djp32(xcut,gamma)
503 : : end if
504 : :
505 : : ! Computes extfpmd contribution to kinetic energy integrating
506 : : ! over energy from ebcut to infinity with order 3/2
507 : : ! incomplete Fermi-Dirac integral.
508 [ + + ]: 90 : if(this%version==1.or.this%version==3) then
509 : 57 : xcut=(this%ebcut-this%eshift)/this%el_temp
510 [ - + ]: 57 : if(this%ebcut.lt.this%eshift) then
511 : 0 : cut_warn=.true.
512 : 0 : xcut=zero
513 : : end if
514 : 57 : this%e_kinetic=this%e_kinetic+factor*djp32(xcut,gamma)
515 : : end if
516 : :
517 : : ! Computes extfpmd contribution to kinetic energy summing
518 : : ! over accessible states from bcut to mband, with
519 : : ! integer band numbers. Total number of bands
520 : : ! is controlled with the input variable extfpmd_nband.
521 [ + + ]: 90 : if(this%version==5) then
522 [ + + ]: 22 : do isppol=1,nsppol
523 [ + + ]: 44 : do ikpt=1,nkpt
524 : 22 : nband_k=nband(ikpt+(isppol-1)*nkpt)
525 [ + - ]: 22 : if(proc_distrb_cycle(this%mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,this%mpi_enreg%me_kpt)) cycle
526 [ - + ]: 22 : if(one*this%bcut+this%bandshiftk(ikpt+(isppol-1)*nkpt).lt.zero) then
527 : 0 : cut_warn=.true.
528 : : end if
529 [ + + ]: 179157 : do iband=nband_k-this%nbdbuf+1,this%mband
530 : 179124 : dotr=extfpmd_e_fg(one*iband+this%bandshiftk(ikpt+(isppol-1)*nkpt),this%ucvol)
531 : 179124 : fn=fermi_dirac(dotr+this%eshift,fermie,this%el_temp)
532 : 179146 : this%e_kinetic=this%e_kinetic+wtk(ikpt)*maxocc*fn*dotr
533 : : end do
534 : : end do
535 : : end do
536 : 11 : call xmpi_sum(this%e_kinetic,this%mpi_enreg%comm_kpt,ierr)
537 : : end if
538 : :
539 : : ! Computes extfpmd contribution to kinetic energy using a sum
540 : : ! of Fermi gas contributions for each point of the fftf grid.
541 : : ! Warning: This is not yet operational. Work in progress.
542 [ - + ]: 90 : if(this%version==10) then
543 [ # # # # : 0 : ABI_MALLOC(gamma_hybrid_tf,(this%nfftf,this%nspden))
# # # # #
# # # # #
# # ]
544 [ # # # # : 0 : ABI_MALLOC(xcut_hybrid_tf,(this%nfftf,this%nspden))
# # # # #
# # # ]
545 [ # # # # ]: 0 : gamma_hybrid_tf(:,:)=(fermie-this%vtrial(:,:))/this%el_temp
546 [ # # # # ]: 0 : xcut_hybrid_tf(:,:)=(this%ebcut-this%vtrial(:,:))/this%el_temp
547 : : e_kinetic_hybrid_tf=zero
548 : :
549 : : !$OMP PARALLEL DO REDUCTION (+:e_kinetic_hybrid_tf)
550 [ # # ]: 0 : do ifft=1,this%nfftf
551 [ # # ]: 0 : do ispden=1,this%nspden
552 : : e_kinetic_hybrid_tf=e_kinetic_hybrid_tf+factor*djp32(xcut_hybrid_tf(ifft,ispden),gamma_hybrid_tf(ifft,ispden))/&
553 : 0 : & (this%nfftf*this%nspden)
554 : : end do
555 : : end do
556 : : !$OMP END PARALLEL DO
557 : :
558 : 0 : this%e_kinetic=e_kinetic_hybrid_tf
559 [ # # # # ]: 0 : gamma_hybrid_tf(:,:)=zero
560 [ # # # # ]: 0 : xcut_hybrid_tf(:,:)=zero
561 [ # # ]: 0 : ABI_FREE(gamma_hybrid_tf)
562 [ # # ]: 0 : ABI_FREE(xcut_hybrid_tf)
563 : : end if
564 : :
565 : : ! Computes the double counting term from the eshift, and
566 : : ! from the contributions to the kinetic energy and
567 : : ! the number of electrons
568 [ - + ]: 90 : if(this%version==10) then
569 [ # # # # ]: 0 : this%edc_kinetic=this%e_kinetic+sum(this%nelectarr(:,:)*this%vtrial(:,:)/(this%nfftf*this%nspden))
570 : : else
571 : 90 : this%edc_kinetic=this%e_kinetic+this%nelect*this%eshift
572 : : end if
573 : :
574 [ + + ]: 90 : if(this%pawsph==1) then
575 : 7 : this%edc_kinetic=this%edc_kinetic+this%nelect*this%eshift_paw
576 : : endif
577 : :
578 [ - + ]: 90 : if(cut_warn) then
579 : : write(msg,'(11a)')&
580 : 0 : & 'Extended FPMD could not properly compute the contribution to the energy.',ch10,&
581 : 0 : & 'This can be due to a too low number of bands in the calculation.',ch10,&
582 : 0 : & 'This can also happen when restarting from a previous calculation.',ch10,&
583 : 0 : & 'Poor prediction of the electron density based on forces may results in this error.',ch10,&
584 : 0 : & 'Action: slightly increase nband if the electron density is supposed to be converged.',ch10,&
585 : 0 : & 'Otherwise: wait for the density to be converged.'
586 : 0 : ABI_WARNING(msg)
587 : : end if
588 : 90 : end subroutine compute_e_kinetic
589 : : !!***
590 : :
591 : : !!****f* ABINIT/m_extfpmd/compute_entropy
592 : : !! NAME
593 : : !! compute_entropy
594 : : !!
595 : : !! FUNCTION
596 : : !! Computes the value of the integral corresponding to the missing
597 : : !! entropy contribution of free electrons after band cut using
598 : : !! incomplete Fermi-Dirac integrals.
599 : : !!
600 : : !! INPUTS
601 : : !! this=extfpmd_type object concerned
602 : : !! fermie=chemical potential (Hartree)
603 : : !! nkpt=number of k points
604 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
605 : : !! nspinor=number of spinor components
606 : : !! wtk(nkpt)=k point weights
607 : : !! nband(nkpt*nsppol)=desired number of bands at each k point
608 : : !!
609 : : !! OUTPUT
610 : : !! this=extfpmd_type object concerned
611 : : !! entropy_extfpmd=extfpmd contribution to the entropy
612 : : !!
613 : : !! SOURCE
614 : 90 : subroutine compute_entropy(this,entropy_extfpmd,fermie,nkpt,nsppol,nspinor,wtk,nband)
615 : : ! Arguments -------------------------------
616 : : ! Scalars
617 : : class(extfpmd_type),intent(inout) :: this
618 : : integer,intent(in) :: nkpt,nsppol,nspinor
619 : : real(dp),intent(in) :: fermie
620 : : real(dp),intent(out) :: entropy_extfpmd
621 : : ! Arrays
622 : : integer,intent(in) :: nband(nkpt*nsppol)
623 : : real(dp),intent(in) :: wtk(nkpt)
624 : :
625 : : ! Local variables -------------------------
626 : : ! Scalars
627 : : integer :: ii,ifft,ispden,isppol,ikpt,iband,nband_k,ierr,nom
628 : : real(dp) :: ix,step,factor,fn,gamma,maxocc
629 : : ! Arrays
630 : 90 : real(dp),dimension(:),allocatable :: valuesent
631 : 90 : real(dp),dimension(:,:),allocatable :: gamma_hybrid_tf
632 : 90 : real(dp),dimension(:,:),allocatable :: step_hybrid_tf
633 : :
634 : : ! *********************************************************************
635 : 90 : maxocc=two/(nsppol*nspinor)
636 : 90 : this%entropy=zero
637 : 90 : factor=dsqrt(two)/(PI*PI)*this%ucvol*this%el_temp**(2.5)
638 : 90 : gamma=(fermie-this%eshift)/this%el_temp
639 [ + - - + : 270 : ABI_MALLOC(valuesent,(this%bcut+1))
+ - + - +
- - + -
+ ]
640 : :
641 : : ! Computes extfpmd contribution to the entropy integrating
642 : : ! over accessible states with Fermi-Dirac complete integrals and
643 : : ! substracting 0 to bcut contribution with numeric integration.
644 [ + + ]: 90 : if(this%version==2.or.this%version==4) then
645 : 22 : step=(dble(this%bcut)+this%bandshift)/(this%bcut)
646 : : !$OMP PARALLEL DO PRIVATE(fn,ix) SHARED(valuesent)
647 [ + + ]: 1144 : do ii=1,this%bcut+1
648 : 1122 : ix=(dble(ii)-one)*step
649 : 1122 : fn=fermi_dirac(extfpmd_e_fg(ix,this%ucvol)+this%eshift,fermie,this%el_temp)
650 [ + - + - ]: 1144 : if(fn>tol16.and.(one-fn)>tol16) then
651 : 1122 : valuesent(ii)=-maxocc*(fn*log(fn)+(one-fn)*log(one-fn))
652 : : else
653 : 0 : valuesent(ii)=zero
654 : : end if
655 : : end do
656 : : !$OMP END PARALLEL DO
657 : :
658 : : ! We need at least 6 elements in valuesent to call simpson function.
659 [ + - ]: 22 : if(size(valuesent)>=6) then
660 : : this%entropy=5./3.*factor*dip32(gamma)/this%el_temp-&
661 : : & gamma*factor*dip12(gamma)/this%el_temp-&
662 : 22 : simpson(step,valuesent)
663 : : end if
664 : : end if
665 : :
666 : : ! Computes extfpmd contribution to the entropy integrating
667 : : ! over energy with Fermi-Dirac complete integrals and
668 : : ! substracting 0 to bcut contribution with numeric integration.
669 [ + + ]: 90 : if(this%version==1.or.this%version==3) then
670 [ - + ]: 57 : ABI_FREE(valuesent)
671 : 57 : nom=10000
672 [ - + ]: 57 : ABI_MALLOC(valuesent,(nom+1))
673 : 57 : step=(this%ebcut-this%eshift)/(nom)
674 : : !$OMP PARALLEL DO PRIVATE(fn,ix) SHARED(valuesent)
675 [ + + ]: 570114 : do ii=1,nom+1
676 : 570057 : ix=this%eshift+(dble(ii)-one)*step
677 : 570057 : fn=fermi_dirac(ix,fermie,this%el_temp)
678 [ + - + - ]: 570114 : if(fn>tol16.and.(one-fn)>tol16) then
679 : : valuesent(ii)=-(fn*log(fn)+(one-fn)*log(one-fn))*&
680 : 570057 : & extfpmd_dos(ix,this%eshift,this%ucvol)
681 : : else
682 : 0 : valuesent(ii)=zero
683 : : end if
684 : : end do
685 : : !$OMP END PARALLEL DO
686 : :
687 : : ! We need at least 6 elements in valuesent to call simpson function.
688 : : if(size(valuesent)>=6) then
689 : : this%entropy=5./3.*factor*dip32(gamma)/this%el_temp-&
690 : 57 : & gamma*factor*dip12(gamma)/this%el_temp-simpson(step,valuesent)
691 : : end if
692 : : end if
693 : :
694 : : ! Computes extfpmd contribution to the entropy summing
695 : : ! over accessible states from bcut to mband, with
696 : : ! integer band numbers. Total number of bands
697 : : ! is controlled with the input variable extfpmd_nband.
698 [ + + ]: 90 : if(this%version==5) then
699 [ + + ]: 22 : do isppol=1,nsppol
700 [ + + ]: 44 : do ikpt=1,nkpt
701 : 22 : nband_k=nband(ikpt+(isppol-1)*nkpt)
702 [ + - ]: 22 : if(proc_distrb_cycle(this%mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,this%mpi_enreg%me_kpt)) cycle
703 [ + + ]: 179157 : do iband=nband_k-this%nbdbuf+1,this%mband
704 : 179124 : fn=fermi_dirac(extfpmd_e_fg(one*iband+this%bandshiftk(ikpt+(isppol-1)*nkpt),this%ucvol)+this%eshift,fermie,this%el_temp)
705 : 179146 : this%entropy=this%entropy-wtk(ikpt)*maxocc*(fn*log(fn)+(one-fn)*log(one-fn))/nsppol
706 : : end do
707 : : end do
708 : : end do
709 : 11 : call xmpi_sum(this%entropy,this%mpi_enreg%comm_kpt,ierr)
710 : : end if
711 : :
712 : : ! Computes extfpmd contribution to the entropy using a sum
713 : : ! of Fermi gas contributions for each point of the fftf grid,
714 : : ! as we do for version=1 and version=2.
715 : : ! Warning: This is not yet operational. Work in progress.
716 [ - + ]: 90 : if(this%version==10) then
717 [ # # # # : 0 : ABI_MALLOC(gamma_hybrid_tf,(this%nfftf,this%nspden))
# # # # #
# # # # #
# # ]
718 [ # # # # : 0 : ABI_MALLOC(step_hybrid_tf,(this%nfftf,this%nspden))
# # # # #
# # # ]
719 [ # # # # ]: 0 : gamma_hybrid_tf(:,:)=(fermie-this%vtrial(:,:))/this%el_temp
720 [ # # # # ]: 0 : step_hybrid_tf(:,:)=(this%ebcut-this%vtrial(:,:))/(this%bcut)
721 : 0 : this%entropy=zero
722 : :
723 [ # # ]: 0 : do ifft=1,this%nfftf
724 [ # # ]: 0 : do ispden=1,this%nspden
725 : : !$OMP PARALLEL DO PRIVATE(ix,fn) SHARED(valuesent)
726 [ # # ]: 0 : do ii=1,this%bcut+1
727 : 0 : ix=this%vtrial(ifft,ispden)+(dble(ii)-one)*step_hybrid_tf(ifft,ispden)
728 : 0 : fn=fermi_dirac(ix,fermie,this%el_temp)
729 [ # # # # ]: 0 : if(fn>tol16.and.(one-fn)>tol16) then
730 : : valuesent(ii)=-(fn*log(fn)+(one-fn)*log(one-fn))*&
731 : 0 : & extfpmd_dos(ix,this%vtrial(ifft,ispden),this%ucvol)
732 : : else
733 : 0 : valuesent(ii)=zero
734 : : end if
735 : : end do
736 : : !$OMP END PARALLEL DO
737 : :
738 : : ! We need at least 6 elements in valuesent to call simpson function.
739 [ # # ]: 0 : if(size(valuesent)>=6) then
740 : : this%entropy=this%entropy+(5./3.*factor*dip32(gamma_hybrid_tf(ifft,ispden))/this%el_temp-&
741 : : & gamma_hybrid_tf(ifft,ispden)*factor*dip12(gamma_hybrid_tf(ifft,ispden))/this%el_temp-&
742 : 0 : & simpson(step_hybrid_tf(ifft,ispden),valuesent))/(this%nfftf*this%nspden)
743 : : end if
744 : : end do
745 : : end do
746 : :
747 [ # # # # ]: 0 : gamma_hybrid_tf(:,:)=zero
748 [ # # # # ]: 0 : step_hybrid_tf(:,:)=zero
749 [ # # ]: 0 : ABI_FREE(step_hybrid_tf)
750 [ # # ]: 0 : ABI_FREE(gamma_hybrid_tf)
751 : : end if
752 [ - + ]: 90 : ABI_FREE(valuesent)
753 : 90 : entropy_extfpmd=this%entropy
754 : 90 : end subroutine compute_entropy
755 : : !!***
756 : :
757 : : !----------------------------------------------------------------------
758 : :
759 : : !!****f* ABINIT/m_extfpmd/extfpmd_dos
760 : : !! NAME
761 : : !! extfpmd_dos
762 : : !!
763 : : !! FUNCTION
764 : : !! Returns the free particle density of states for a given energy.
765 : : !!
766 : : !! INPUTS
767 : : !! energy=get the value of the free particle density of states at this energy
768 : : !! eshift=energy shift factor
769 : : !! ucvol=unit cell volume (bohr^3)
770 : : !!
771 : : !! OUTPUT
772 : : !! extfpmd_dos=value of free particle density of states at given energy
773 : : !!
774 : : !! SOURCE
775 : 570057 : function extfpmd_dos(energy,eshift,ucvol)
776 : : ! Arguments -------------------------------
777 : : ! Scalars
778 : : real(dp),intent(in) :: energy,eshift,ucvol
779 : : real(dp) :: extfpmd_dos
780 : :
781 : : ! *********************************************************************
782 : :
783 : 570057 : extfpmd_dos=dsqrt(two)*ucvol*dsqrt(energy-eshift)/(PI*PI)
784 : 570057 : end function extfpmd_dos
785 : : !!***
786 : :
787 : : !!****f* ABINIT/m_extfpmd/extfpmd_e_fg
788 : : !! NAME
789 : : !! extfpmd_e_fg
790 : : !!
791 : : !! FUNCTION
792 : : !! Returns the energy of the Fermi gas for a given number of
793 : : !! accessible states.
794 : : !!
795 : : !! INPUTS
796 : : !! iband=number of accessible states
797 : : !! ucvol=unit cell volume (bohr^3)
798 : : !!
799 : : !! OUTPUT
800 : : !! extfpmd_e_fg=energy of homogeneous electron gas for a given number of accessible states
801 : : !!
802 : : !! SOURCE
803 : 9838051 : function extfpmd_e_fg(iband,ucvol)
804 : : ! Arguments -------------------------------
805 : : ! Scalars
806 : : real(dp),intent(in) :: iband,ucvol
807 : : real(dp) :: extfpmd_e_fg
808 : :
809 : : ! *********************************************************************
810 : :
811 : 9838051 : extfpmd_e_fg=half*(iband*six*PI*PI/ucvol)**(two/three)
812 : 9838051 : end function extfpmd_e_fg
813 : :
814 : : !!***
815 : : !!****f* ABINIT/m_extfpmd/extfpmd_i_fg
816 : : !! NAME
817 : : !! extfpmd_i_fg
818 : : !!
819 : : !! FUNCTION
820 : : !! Returns the number of doubly occupied orbitals of the
821 : : !! for a Fermi Gas for a given kinetic energy.
822 : : !!
823 : : !! INPUTS
824 : : !! ekin=kinetic energy
825 : : !! ucvol=unit cell volume (bohr^3)
826 : : !!
827 : : !! OUTPUT
828 : : !! extfpmd_i_fg=number of doubly occupied states of the Fermi gas
829 : : !!
830 : : !! SOURCE
831 : 346 : function extfpmd_i_fg(ekin,ucvol)
832 : : ! Arguments -------------------------------
833 : : ! Scalars
834 : : real(dp),intent(in) :: ekin,ucvol
835 : : real(dp) :: extfpmd_i_fg
836 : :
837 : : ! *********************************************************************
838 : :
839 : 346 : extfpmd_i_fg=(two*ekin)**(three/two)*ucvol/(six*PI*PI)
840 : 346 : end function extfpmd_i_fg
841 : : !!***
842 : :
843 : : !!***
844 : : !!****f* ABINIT/m_extfpmd/extfpmd_chkinp
845 : : !! NAME
846 : : !! extfpmd_chkinp
847 : : !!
848 : : !! FUNCTION
849 : : !! Returns true if extfpmd input variables are valid and false otherwise.
850 : : !!
851 : : !! INPUTS
852 : : !! dtset=<type datafiles_type>contains all input variables.
853 : : !!
854 : : !! OUTPUT
855 : : !! extfpmd_chkinp=allocate extfpmd object or not
856 : : !!
857 : : !! SOURCE
858 : 7 : function extfpmd_chkinp(dtset)
859 : : ! Arguments -------------------------------
860 : : ! Scalars
861 : : class(dataset_type),intent(in) :: dtset
862 : : logical :: extfpmd_chkinp
863 : : ! Local variables -------------------------
864 : : ! Scalars
865 : : character(len=500) :: msg
866 : :
867 : : ! *********************************************************************
868 : :
869 : 7 : extfpmd_chkinp=.false.
870 [ - + ]: 7 : if(.not.(dtset%occopt>=3.and.dtset%occopt<=9)) then
871 : 0 : write(msg,'(3a)') "ExtFPMD routines need metallic occupation option.",ch10,&
872 : 0 : & "Action: Set occopt input variable to a value >= 3 and <= 9."
873 : 0 : ABI_ERROR(msg)
874 [ + + - + ]: 7 : else if((dtset%useextfpmd==2.or.dtset%useextfpmd==3).and.(dtset%extfpmd_nbcut>dtset%mband)) then
875 : 0 : write(msg,'(3a,i0,a,i0,3a)') "Not enough bands to activate ExtFPMD routines.",ch10,&
876 : 0 : & "extfpmd_nbcut = ",dtset%extfpmd_nbcut," should be less than or equal to nband = ",dtset%mband,".",ch10,&
877 : 0 : & "Action: Increase nband or decrease extfpmd_nbcut."
878 : 0 : ABI_ERROR(msg)
879 [ + + - + ]: 7 : else if((dtset%useextfpmd==2.or.dtset%useextfpmd==3).and.(dtset%extfpmd_nbdbuf+dtset%extfpmd_nbcut>dtset%mband)) then
880 : 0 : write(msg,'(a,i0,a,i0,a,i0,2a,i0,a)') "(extfpmd_nbdbuf = ",dtset%extfpmd_nbdbuf," + extfpmd_nbcut = ",&
881 : 0 : & dtset%extfpmd_nbcut,") = ",dtset%extfpmd_nbdbuf+dtset%extfpmd_nbcut,ch10,&
882 : 0 : & "should be less than or equal to nband = ",dtset%mband,"."
883 : 0 : ABI_ERROR(msg)
884 [ - + ]: 7 : else if(dtset%extfpmd_nbdbuf>dtset%mband) then
885 : 0 : write(msg,'(a,i0,a,i0,a)') "extfpmd_nbdbuf = ",dtset%extfpmd_nbdbuf,&
886 : 0 : & " should be less than or equal to nband = ",dtset%mband,"."
887 : 0 : ABI_ERROR(msg)
888 [ + + - + ]: 7 : else if((dtset%useextfpmd==5.or.dtset%useextfpmd==11).and.(dtset%extfpmd_nband<=dtset%mband)) then
889 : 0 : write(msg,'(3a,i0,a,i0,3a)') "Not enough bands to activate ExtFPMD routines.",ch10,&
890 : 0 : & "extfpmd_nband = ",dtset%extfpmd_nband," should be strictly greater than nband = ",dtset%mband,".",ch10,&
891 : 0 : & "Action: Increase extfpmd_nband or decrease nband."
892 : 0 : ABI_ERROR(msg)
893 : : else
894 : : extfpmd_chkinp=.true.
895 : : end if
896 : 7 : end function extfpmd_chkinp
897 : : !!***
898 : :
899 : : !!***
900 : : !!****f* ABINIT/m_extfpmd/extfpmd_prterr
901 : : !! NAME
902 : : !! extfpmd_prterr
903 : : !!
904 : : !! SUBROUTINE
905 : : !! Print extfpmd error
906 : : !!
907 : : !! INPUTS
908 : : !!
909 : : !! OUTPUT
910 : : !!
911 : : !! SOURCE
912 : 0 : subroutine extfpmd_err(this,eigen,mband,nband,nkpt,nsppol,wtk,fname)
913 : : ! Arguments -------------------------------
914 : : ! Scalars
915 : : class(extfpmd_type),intent(inout) :: this
916 : : character(len=*),intent(in) :: fname
917 : : integer,intent(in) :: mband,nkpt,nsppol
918 : : ! Arrays
919 : : integer,intent(in) :: nband(nkpt*nsppol)
920 : : real(dp),intent(in) :: eigen(mband*nkpt*nsppol)
921 : : real(dp),intent(in) :: wtk(nkpt)
922 : :
923 : : ! Local variables -------------------------
924 : : ! Scalars
925 : : integer :: band_index,ii,ikpt,isppol,nband_k,tmp_unt
926 [ # # ]: 0 : real(dp) :: err(minval(nband)-this%nbdbuf)
927 : :
928 : : ! *********************************************************************
929 [ # # ]: 0 : err=zero
930 : : band_index=0
931 [ # # ]: 0 : do isppol=1,nsppol
932 [ # # ]: 0 : do ikpt=1,nkpt
933 : 0 : nband_k=nband(ikpt+(isppol-1)*nkpt)
934 [ # # ]: 0 : do ii=1,nband_k-this%nbdbuf
935 : : err(ii)=err(ii)+&
936 : : & wtk(ikpt)*(eigen(band_index+ii)-this%eshift-extfpmd_e_fg(dble(ii),this%ucvol))&
937 : 0 : & /nsppol/extfpmd_e_fg(dble(ii),this%ucvol)
938 : : end do
939 : 0 : band_index=band_index+nband_k
940 : : end do
941 : : end do
942 : 0 : open(file=fname,newunit=tmp_unt,status='unknown',form='formatted')
943 [ # # # # ]: 0 : do ii=1,minval(nband)-this%nbdbuf
944 : 0 : write(tmp_unt,*) ii,abs(err(ii))
945 : : enddo
946 : 0 : close(tmp_unt)
947 : 0 : end subroutine extfpmd_err
948 : : !!***
949 : :
950 [ # # # # : 0 : end module m_extfpmd
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
951 : : !!***
|
