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