PROGRESS  master
prg_densitymatrix_mod.F90
Go to the documentation of this file.
1 
5 module prg_densitymatrix_mod
6 
7  use bml
8  use iso_c_binding
9  use prg_parallel_mod
10 
11 #ifdef USE_NVTX
12  use prg_nvtx_mod
13 #endif
14 
15  implicit none
16 
17  private !Everything is private by default
18 
19  integer, parameter :: dp = kind(1.0d0)
20 
21  public :: prg_build_density_t0, prg_check_idempotency, prg_get_eigenvalues
22  public :: prg_get_flevel, prg_build_density_t, prg_build_atomic_density
23  public :: prg_build_density_t_fermi, prg_get_flevel_nt, prg_build_density_t_fulldata
24  public :: prg_build_density_t_ed, prg_toeigenspace, prg_tocanonicalspace
25  public :: canon_dm_prt, prg_get_evalsdvalsevects, prg_build_density_fromevalsandevects
26 
27 contains
28 
41  subroutine prg_build_density_t0(ham_bml, rho_bml, threshold, bndfil, eigenvalues_out)
42 
43  character(20) :: bml_type
44  integer :: i, norb
45  real(8), intent(in) :: bndfil, threshold
46  real(dp) :: nocc
47  real(dp), allocatable :: eigenvalues(:)
48  real(dp), allocatable, optional, intent(out) :: eigenvalues_out(:)
49  type(bml_matrix_t) :: aux1_bml, aux2_bml, eigenvectors_bml
50  type(bml_matrix_t), intent(in) :: ham_bml
51  type(bml_matrix_t), intent(inout) :: rho_bml
52  character(20) :: bml_dmode
53 
54  if (getnranks().gt.1)then
55  if (printrank() .eq. 1) print*,'prg_build_density_T0: BML_DMODE_DISTRIBUTED'
56  bml_dmode = bml_dmode_distributed
57  else
58  print*,'prg_build_density_T0: BML_DMODE_SEQUENTIAL'
59  bml_dmode = bml_dmode_sequential
60  endif
61 
62  norb = bml_get_n(ham_bml)
63  bml_type = bml_get_deep_type(ham_bml)
64 
65  allocate(eigenvalues(norb))
66 
67  call bml_noinit_matrix(bml_type,bml_element_real,dp,norb,norb,eigenvectors_bml,bml_dmode)
68 
69  call bml_diagonalize(ham_bml,eigenvalues,eigenvectors_bml)
70  if(present(eigenvalues_out))then
71  if(allocated(eigenvalues_out))deallocate(eigenvalues_out)
72  allocate(eigenvalues_out(norb))
73  eigenvalues_out = eigenvalues
74  endif
75 
76  nocc = norb*bndfil
77 
78  do i=1,norb !Reusing eigenvalues to apply the theta function.
79  if(real(i)-nocc < 0.0001_dp) then
80  eigenvalues(i) = 2.0_dp
81  elseif(abs(real(i)-real(nocc)) < 0.0001_dp) then
82  eigenvalues(i) = 2.0_dp
83  else
84  eigenvalues(i) = 0.0_dp
85  endif
86  enddo
87  if(abs(nocc - int(nocc)) > 0.01_dp)then
88  eigenvalues(int(nocc)+1) = 1.0_dp
89  endif
90 
91  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml,bml_dmode)
92  call bml_set_diagonal(aux1_bml, eigenvalues) !eps(i,i) = eps(i)
93 
94  deallocate(eigenvalues)
95 
96  call bml_noinit_matrix(bml_type,bml_element_real,dp,norb,norb,aux2_bml,bml_dmode)
97  call bml_multiply(eigenvectors_bml, aux1_bml, aux2_bml, 1.0_dp, 0.0_dp, threshold)
98 
99  call bml_transpose_new(eigenvectors_bml, aux1_bml)
100  call bml_deallocate(eigenvectors_bml)
101 
102  call bml_multiply(aux2_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
103  call bml_deallocate(aux2_bml)
104  call bml_deallocate(aux1_bml)
105 
106  end subroutine prg_build_density_t0
107 
122  subroutine prg_build_density_t(ham_bml, rho_bml, threshold, bndfil, kbt, ef, eigenvalues_out)
123 
124  character(20) :: bml_type
125  integer :: i, norb
126  real(dp), intent(in) :: bndfil, threshold, kbt
127  real(dp), intent(inout) :: ef
128  real(dp) :: nocc, fleveltol, mlsi, efold
129  real(dp), allocatable :: eigenvalues(:)
130  real(dp), allocatable, optional, intent(out) :: eigenvalues_out(:)
131  type(bml_matrix_t) :: aux1_bml, aux2_bml, eigenvectors_bml
132  type(bml_matrix_t), intent(in) :: ham_bml
133  type(bml_matrix_t), intent(inout) :: rho_bml
134  logical :: err
135 
136  if (printrank() .eq. 1) then
137  write(*,*)"In get_density_t ..."
138  endif
139 
140  norb = bml_get_n(ham_bml)
141  bml_type = bml_get_deep_type(ham_bml)
142 
143  allocate(eigenvalues(norb))
144  call bml_noinit_matrix(bml_type,bml_element_real,dp,norb,norb,eigenvectors_bml)
145 
146  call bml_diagonalize(ham_bml,eigenvalues,eigenvectors_bml)
147 
148  if(present(eigenvalues_out))then
149  if(allocated(eigenvalues_out))deallocate(eigenvalues_out)
150  allocate(eigenvalues_out(norb))
151  eigenvalues_out = eigenvalues
152  endif
153 
154  fleveltol = 1.0e-11
155 
156  efold = ef
157  call prg_get_flevel_nt(eigenvalues,kbt,bndfil,fleveltol,ef,err)
158  if(err)then
159  ef = efold
160  call prg_get_flevel(eigenvalues,kbt,bndfil,fleveltol,ef,err)
161  endif
162  if(err)then
163  write(*,*)"WARNING: Ef/Chemical potential search failed. We'll use the previous one to proceed"
164  ef = efold
165  endif
166 
167  nocc = norb*bndfil
168 
169  do i=1,norb !Reusing eigenvalues to apply the theta function.
170  eigenvalues(i) = 2.0_dp*fermi(eigenvalues(i),ef,kbt)
171  enddo
172 
173  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml)
174  call bml_set_diagonal(aux1_bml, eigenvalues) !eps(i,i) = eps(i)
175 
176  deallocate(eigenvalues)
177 
178  call bml_noinit_matrix(bml_type,bml_element_real,dp,norb,norb,aux2_bml)
179  call bml_multiply(eigenvectors_bml, aux1_bml, aux2_bml, 1.0_dp, 0.0_dp, threshold)
180 
181  call bml_transpose_new(eigenvectors_bml, aux1_bml)
182  call bml_deallocate(eigenvectors_bml)
183 
184  call bml_multiply(aux2_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
185  call bml_deallocate(aux2_bml)
186  call bml_deallocate(aux1_bml)
187 
188  end subroutine prg_build_density_t
189 
190 
207  subroutine prg_build_density_t_fulldata(ham_bml, rho_bml, threshold, bndfil, kbt, ef, eigenvalues_out&
208  &, evects_bml, fvals)
209 
210  character(20) :: bml_type
211  integer :: i, norb
212  real(dp), intent(in) :: bndfil, threshold, kbt
213  real(dp), intent(inout) :: ef
214  real(dp) :: nocc, fleveltol, efold
215  real(dp), allocatable :: eigenvalues(:)
216  real(dp), allocatable, intent(inout) :: eigenvalues_out(:), fvals(:)
217  type(bml_matrix_t) :: aux1_bml, aux_bml, occupation_bml
218  type(bml_matrix_t), intent(in) :: ham_bml
219  type(bml_matrix_t), intent(inout) :: rho_bml, evects_bml
220  logical :: err
221 
222  if (printrank() .eq. 1) then
223  write(*,*)"In get_density_t_fulldata ..."
224  endif
225 
226  norb = bml_get_n(ham_bml)
227  bml_type = bml_get_deep_type(ham_bml)
228 
229  allocate(eigenvalues(norb))
230  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,evects_bml)
231 
232  call bml_diagonalize(ham_bml,eigenvalues,evects_bml)
233 
234  if(.not.allocated(eigenvalues_out))then
235  allocate(eigenvalues_out(norb))
236  allocate(fvals(norb))
237  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,evects_bml)
238  endif
239 
240  eigenvalues_out = eigenvalues
241 
242  fleveltol = 1.0e-11
243  fvals = 0.0_dp
244 
245  call prg_get_flevel_nt(eigenvalues,kbt,bndfil,fleveltol,ef,err)
246  if(err)call prg_get_flevel(eigenvalues,kbt,bndfil,fleveltol,ef,err)
247  if(err)then
248  write(*,*)"WARNING: Ef/Chemical potential search failed. We'll use the previous one to proceed"
249  ef = efold
250  endif
251 
252  nocc = norb*bndfil
253 
254  do i=1,norb !Aapply Fermi function.
255  fvals(i) = 2.0_dp*fermi(eigenvalues(i),ef,kbt)
256  enddo
257 
258  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,occupation_bml)
259  call bml_set_diagonal(occupation_bml, fvals) !eps(i,i) = eps(i)
260 
261  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux_bml)
262  call bml_multiply(evects_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp, threshold)
263  call bml_deallocate(occupation_bml)
264 
265  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml)
266  call bml_transpose_new(evects_bml, aux1_bml)
267 
268  call bml_multiply(aux_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
269  call bml_deallocate(aux_bml)
270  call bml_deallocate(aux1_bml)
271 
272  end subroutine prg_build_density_t_fulldata
273 
274 
292  subroutine prg_build_density_t_ed(ham_bml, rho_bml, evects_bml, threshold, bndfil, kbt, ef, evals&
293  &, dvals, hindex, llsize, verbose)
294 
295  character(20) :: bml_type
296  integer, allocatable, intent(in) :: hindex(:,:)
297  integer :: i, j, jj, k, kk, l, norb, norbcore
298  integer, intent(in) :: llsize, verbose
299  real(dp), intent(in) :: bndfil, threshold, kbt
300  real(dp), intent(inout) :: ef
301  real(dp) :: nocc
302  real(dp), allocatable :: eigenvalues(:), row(:),fvals(:),ham(:,:),evects(:,:)
303  real(dp), allocatable, intent(inout) :: evals(:), dvals(:)
304  type(bml_matrix_t) :: aux1_bml, aux_bml, occupation_bml
305  type(bml_matrix_t), intent(in) :: ham_bml
306  type(bml_matrix_t), intent(inout) :: rho_bml, evects_bml
307 
308  if (printrank() .eq. 1 .and. verbose >= 1) then
309  write(*,*)"In get_density_t_efd ..."
310  endif
311 
312  norb = bml_get_n(ham_bml)
313  bml_type = bml_get_type(ham_bml)
314 
315  if(allocated(evals))then
316  deallocate(evals)
317  deallocate(dvals)
318  endif
319  allocate(evals(norb))
320  allocate(dvals(norb))
321  allocate(fvals(norb))
322  allocate(ham(norb,norb))
323  allocate(evects(norb,norb))
324  call bml_export_to_dense(ham_bml,ham)
325  if(bml_get_n(evects_bml) < 0)then
326  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,evects_bml)
327  endif
328  call eig(ham,evects,evals,'V',norb)
329  call bml_import_from_dense(bml_type,evects,evects_bml,0.0_dp,norb)
330  !call bml_diagonalize(ham_bml,evals,evects_bml)
331  write(*,*)"ham",ham(1,1),ham(1,2),ham(1,3)
332  deallocate(evects)
333  deallocate(ham)
334  call bml_print_matrix("evects",evects_bml,0,4,0,4)
335  write(*,*)"Evals",evals(1),evals(2),evals(3)
336  nocc = norb*bndfil
337 
338  do i=1,norb !Apply Fermi function.
339  fvals(i) = 2.0_dp*fermi(evals(i),ef,kbt)
340  enddo
341 
342  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,occupation_bml)
343  call bml_set_diagonal(occupation_bml, fvals) !eps(i,i) = eps(i)
344 
345  deallocate(fvals)
346 
347  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux_bml)
348  call bml_multiply(evects_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp, threshold)
349  call bml_deallocate(occupation_bml)
350 
351  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml)
352  call bml_transpose_new(evects_bml, aux1_bml)
353 
354  allocate(row(norb))
355  dvals = 0.0_dp
356  do i = 1,norb
357  call bml_get_row(aux1_bml,i,row)
358  do k = 1,llsize
359  do l = hindex(1,k),hindex(2,k)
360  dvals(i) = dvals(i) + row(l)**2
361  enddo
362  enddo
363  enddo
364  deallocate(row)
365 
366  call bml_multiply(aux_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
367 
368  call bml_deallocate(aux_bml)
369  call bml_deallocate(aux1_bml)
370 
371  end subroutine prg_build_density_t_ed
372 
373 
386  subroutine prg_get_evalsdvalsevects(ham_bml, threshold, hindex,&
387  & llsize, evals, dvals, evects_bml, verbose)
388 
389  character(20) :: bml_type
390  integer, allocatable, intent(in) :: hindex(:,:)
391  integer :: i, k, l, norb, nats
392  integer, intent(in) :: llsize, verbose
393  real(dp), intent(in) :: threshold
394  real(dp), allocatable :: row(:),ham(:,:),evects(:,:),aux(:,:)
395  real(dp), allocatable, intent(inout) :: evals(:), dvals(:)
396  type(bml_matrix_t), intent(in) :: ham_bml
397  type(bml_matrix_t), intent(inout) :: evects_bml
398  type(bml_matrix_t) :: aux1_bml
399 #ifdef USE_OFFLOAD
400  type(c_ptr) :: evects_bml_c_ptr
401  integer :: ld
402  real(c_double), pointer :: evects_bml_ptr(:,:)
403  real(dp) :: this_dval
404 #endif
405 
406  if (printrank() .eq. 1 .and. verbose >= 1) then
407  write(*,*)"In get_evalsDvalsEvects ..."
408  endif
409 
410  norb = bml_get_n(ham_bml)
411  bml_type = bml_get_type(ham_bml)
412 
413  if(.not.allocated(evals))allocate(evals(norb))
414  if(.not.allocated(dvals))allocate(dvals(norb))
415  !allocate(ham(norb,norb))
416  !allocate(evects(norb,norb))
417  !call bml_export_to_dense(ham_bml,ham)
418  if(.not. bml_allocated(evects_bml))then
419  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,evects_bml)
420  endif
421  !call Eig(ham,evects,evals,'V',norb)
422  !call bml_import_from_dense(bml_type,evects,evects_bml,0.0_dp,norb)
423  call bml_diagonalize(ham_bml,evals,evects_bml)
424  !deallocate(evects)
425  !deallocate(ham)
426  dvals = 0.0_dp
427 #ifdef USE_OFFLOAD
428  evects_bml_c_ptr = bml_get_data_ptr_dense(evects_bml)
429  ld = bml_get_ld_dense(evects_bml)
430  nats = size(hindex,dim=2)
431  call c_f_pointer(evects_bml_c_ptr,evects_bml_ptr,shape=[ld,norb])
432  !$acc enter data copyin(dvals(1:norb),hindex(1:2,1:nats))
433  !$acc parallel loop deviceptr(evects_bml_ptr) &
434  !$acc present(hindex,dvals) &
435  !$acc private(i,k,l,this_dval)
436  do i = 1,norb
437  this_dval = 0.0_dp
438  !$acc loop reduction(+:this_dval)
439  do k = 1,llsize
440  do l = hindex(1,k),hindex(2,k)
441  this_dval = this_dval + evects_bml_ptr(i,l)*evects_bml_ptr(i,l)
442  enddo
443  enddo
444  dvals(i) = this_dval
445  enddo
446  !$acc exit data copyout(dvals(1:norb)) delete(hindex(1:2,1:nats))
447 #else
448 ! call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml)
449 ! call bml_transpose_new(evects_bml, aux1_bml)
450 ! allocate(aux(norb,norb))
451 ! call bml_export_to_dense(aux1_bml,aux)
452 ! call bml_deallocate(aux1_bml)
453  call bml_export_to_dense(evects_bml,aux)
454  allocate(row(norb))
455  dvals = 0.0_dp
456  do i = 1,norb
457  !call bml_get_row(aux1_bml,i,row)
458  row(:) = aux(:,i)
459  do k = 1,llsize
460  do l = hindex(1,k),hindex(2,k)
461  dvals(i) = dvals(i) + row(l)**2
462  enddo
463  enddo
464  enddo
465  deallocate(row)
466  deallocate(aux)
467 #endif
468 
469  end subroutine prg_get_evalsdvalsevects
470 
471 
481  subroutine prg_build_density_fromevalsandevects(evects_bml, evals, rho_bml, threshold,&
482  & bndfil, kbt, ef, verbose)
483 
484  character(20) :: bml_type
485  integer :: i, j, norb
486  integer, intent(in) :: verbose
487  real(dp), intent(in) :: bndfil, threshold, kbt
488  real(dp), intent(inout) :: ef
489  real(dp) :: nocc
490  real(dp), intent(in) :: evals(*)
491  real(dp), allocatable :: fvals(:)
492  type(bml_matrix_t) :: aux_bml, occupation_bml
493  type(bml_matrix_t), intent(inout) :: evects_bml
494  type(bml_matrix_t), intent(inout) :: rho_bml
495 #ifdef USE_OFFLOAD
496  type(c_ptr) :: evects_bml_c_ptr, aux_bml_c_ptr
497  integer :: ld
498  real(c_double), pointer :: evects_bml_ptr(:,:), aux_bml_ptr(:,:)
499 #endif
500 
501  if (printrank() .eq. 1 .and. verbose >= 1) then
502  write(*,*)"In get_density_fromEvalsAndEvects ..."
503  endif
504 
505  norb = bml_get_n(evects_bml)
506  bml_type = bml_get_type(evects_bml)
507 
508  allocate(fvals(norb))
509 
510  nocc = norb*bndfil
511 
512  !$omp parallel do simd shared(fvals,evals,ef,kbt)
513  do i=1,norb !Apply Fermi function.
514  fvals(i) = 2.0_dp*fermi(evals(i),ef,kbt)
515  enddo
516 
517 #ifdef USE_OFFLOAD
518 
519  call bml_noinit_matrix(bml_type,bml_element_real,dp,norb,norb,aux_bml)
520 
521  evects_bml_c_ptr = bml_get_data_ptr_dense(evects_bml)
522  aux_bml_c_ptr = bml_get_data_ptr_dense(aux_bml)
523  ld = bml_get_ld_dense(evects_bml)
524  call c_f_pointer(evects_bml_c_ptr,evects_bml_ptr,shape=[ld,norb])
525  call c_f_pointer(aux_bml_c_ptr,aux_bml_ptr,shape=[ld,norb])
526  !$acc enter data copyin(fvals(:))
527  !$acc parallel loop gang deviceptr(aux_bml_ptr,evects_bml_ptr) present(fvals)
528  do i = 1,norb
529  !$acc loop vector
530  do j = 1,norb
531  aux_bml_ptr(j,i) = evects_bml_ptr(i,j)*fvals(i)
532  enddo
533  enddo
534  !$acc exit data delete(fvals(:))
535 #else
536  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,occupation_bml)
537  call bml_set_diagonal(occupation_bml, fvals) !eps(i,i) = eps(i)
538 
539  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux_bml)
540  call bml_multiply(evects_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp,threshold)
541  call bml_deallocate(occupation_bml)
542 
543  call bml_transpose(aux_bml)
544 #endif
545 
546  deallocate(fvals)
547 
548  call bml_multiply(evects_bml, aux_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
549 
550  call bml_deallocate(aux_bml)
551 
552  end subroutine prg_build_density_fromevalsandevects
553 
554 
555 
566  subroutine prg_build_density_t_fermi(ham_bml, rho_bml, threshold, kbt, ef, verbose, drho)
567 
568  character(20) :: bml_type
569  integer :: i, norb, mdim
570  integer, optional, intent(in) :: verbose
571  real(dp), intent(in) :: threshold, kbt
572  real(dp), intent(in) :: ef
573  real(dp) :: fleveltol
574  real(dp), allocatable :: eigenvalues(:)
575  real(dp), allocatable :: nocc(:)
576  type(bml_matrix_t) :: aux1_bml, aux_bml, eigenvectors_bml, occupation_bml
577  type(bml_matrix_t), intent(in) :: ham_bml
578  type(bml_matrix_t), intent(inout) :: rho_bml
579  real(dp), optional, intent(inout) :: drho ! gradient of density w.r.t. ef
580 
581  if (printrank() .eq. 1) then
582  if(present(verbose))then
583  if(verbose >= 1)then
584  write(*,*)"In get_density_t_Fermi ..."
585  write(*,*)"Ef = ",ef
586  endif
587  endif
588  endif
589 
590  norb = bml_get_n(ham_bml)
591  mdim = bml_get_m(ham_bml)
592  bml_type = bml_get_deep_type(ham_bml)
593 
594  allocate(eigenvalues(norb),nocc(norb))
595  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,eigenvectors_bml)
596 
597  call bml_diagonalize(ham_bml,eigenvalues,eigenvectors_bml)
598 
599  do i=1,norb !Reusing eigenvalues to apply the theta function.
600  nocc(i) = 2.0_dp*fermi(eigenvalues(i),ef,kbt)
601  enddo
602 
603  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,occupation_bml)
604  call bml_set_diagonal(occupation_bml, nocc) !eps(i,i) = eps(i)
605 
606  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux_bml)
607  call bml_multiply(eigenvectors_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp, threshold)
608  call bml_deallocate(occupation_bml)
609 
610  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux1_bml)
611  call bml_transpose_new(eigenvectors_bml, aux1_bml)
612 
613  call bml_multiply(aux_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
614 
615  if (present(drho)) then
616  do i=1,norb !Reusing eigenvalues to apply the theta function.
617  nocc(i) = 2.0_dp * dfermi(eigenvalues(i),ef,kbt)
618  enddo
619 
620  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,occupation_bml)
621  call bml_set_diagonal(occupation_bml, nocc) !eps(i,i) = eps(i)
622 
623  call bml_multiply(eigenvectors_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp, threshold)
624 
625  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux1_bml)
626  call bml_transpose_new(eigenvectors_bml, aux1_bml)
627 
628  call bml_multiply(aux_bml, aux1_bml, occupation_bml, 1.0_dp, 0.0_dp, threshold)
629  drho = bml_trace(occupation_bml)
630  call bml_deallocate(occupation_bml)
631  endif
632 
633  deallocate(nocc)
634  deallocate(eigenvalues)
635  call bml_deallocate(eigenvectors_bml)
636  call bml_deallocate(aux_bml)
637  call bml_deallocate(aux1_bml)
638 
639  end subroutine prg_build_density_t_fermi
640 
650  subroutine prg_build_atomic_density(rhoat_bml,numel,hindex,spindex,norb,bml_type)
651 
652  character(len=*), intent(in) :: bml_type
653  integer :: i, j, index, n_orb, nats
654  integer, intent(in) :: hindex(:,:), norb, spindex(:)
655  real(dp) :: occ
656  real(dp), allocatable :: d_atomic(:)
657  real(dp), intent(in) :: numel(:)
658  type(bml_matrix_t), intent(inout) :: rhoat_bml
659 #ifdef USE_OFFLOAD
660  type(c_ptr) :: rhoat_bml_c_ptr
661  integer :: ld
662  real(c_double), pointer :: rhoat_bml_ptr(:,:)
663 #else
664  real(dp), allocatable :: rhoat(:)
665 #endif
666 
667  nats = size(hindex,dim=2)
668 
669 #ifdef USE_OFFLOAD
670  rhoat_bml_c_ptr = bml_get_data_ptr_dense(rhoat_bml)
671  ld = bml_get_ld_dense(rhoat_bml)
672 
673  call c_f_pointer(rhoat_bml_c_ptr,rhoat_bml_ptr,shape=[ld,norb])
674 
675  !$acc enter data copyin(hindex(:,:),spindex(:),numel(:))
676 
677  !$acc parallel loop gang collapse(2) deviceptr(rhoat_bml_ptr)
678  do i = 1,norb
679  do j = 1,norb
680  rhoat_bml_ptr(i,j) = 0.0_dp
681  enddo
682  enddo
683 
684  !$acc parallel loop gang deviceptr(rhoat_bml_ptr) &
685  !$acc present(hindex,spindex,numel) &
686  !$acc private(i,index,n_orb,occ)
687 
688  do i = 1,nats
689  index = hindex(1,i) - 1
690  n_orb = hindex(2,i)-hindex(1,i) + 1;
691  if(n_orb == 1)then
692  index = index + 1;
693  rhoat_bml_ptr(index,index) = numel(spindex(i));
694  else
695  if(numel(spindex(i)) <= 2)then
696  index = index + 1;
697  rhoat_bml_ptr(index,index) = numel(spindex(i));
698  index = index + 1;
699  rhoat_bml_ptr(index,index) = 0.0_dp;
700  index = index + 1;
701  rhoat_bml_ptr(index,index) = 0.0_dp;
702  index = index + 1;
703  rhoat_bml_ptr(index,index) = 0.0_dp;
704  else
705  index = index + 1;
706  rhoat_bml_ptr(index,index) = 2.0_dp;
707 
708  index = index + 1;
709  occ = (numel(spindex(i))-2.0_dp)/3.0_dp;
710  rhoat_bml_ptr(index,index) = occ;
711  index = index + 1;
712  rhoat_bml_ptr(index,index) = occ;
713  index = index + 1;
714  rhoat_bml_ptr(index,index) = occ;
715  endif
716  endif
717  enddo
718 
719  !$acc exit data delete(hindex(:,:),spindex(:),numel(:))
720 
721 #else
722  allocate(rhoat(norb))
723 
724  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,rhoat_bml)
725 
726  index = 0;
727 
728  do i = 1,nats
729  n_orb = hindex(2,i)-hindex(1,i) + 1;
730  if(n_orb == 1)then
731  index = index + 1;
732  rhoat(index) = numel(spindex(i));
733  else
734  if(numel(spindex(i)) <= 2)then
735  index = index + 1;
736  rhoat(index) = numel(spindex(i));
737  index = index + 1;
738  rhoat(index) = 0.0_dp;
739  index = index + 1;
740  rhoat(index) = 0.0_dp;
741  index = index + 1;
742  rhoat(index) = 0.0_dp;
743  else
744  index = index + 1;
745  rhoat(index) = 2.0_dp;
746 
747  index = index + 1;
748  occ = (numel(spindex(i))-2.0_dp)/3.0_dp;
749  rhoat(index) = occ;
750  index = index + 1;
751  rhoat(index) = occ;
752  index = index + 1;
753  rhoat(index) = occ;
754  endif
755  endif
756  enddo
757 
758  call bml_set_diagonal(rhoat_bml,rhoat,0.0_dp)
759 
760  deallocate(rhoat)
761 #endif
762 
763  end subroutine prg_build_atomic_density
764 
775  subroutine prg_get_flevel(eigenvalues,kbt,bndfil,tol,Ef,err)
776 
777  integer :: i, j, k, m
778  integer :: norb
779  real(dp) :: ft1, ft2, kb, prod
780  real(dp) :: t, step, tol, nel
781  real(dp), intent(in) :: bndfil, eigenvalues(:), kbt
782  real(dp), intent(inout) :: ef
783  logical, intent(inout) :: err
784 
785  norb = size(eigenvalues,dim=1)
786  nel = bndfil*2.0_dp*norb
787  ef=eigenvalues(1)
788  step=abs((eigenvalues(norb)-eigenvalues(1)))
789  ft1=0.0_dp
790  ft2=0.0_dp
791  prod=0.0_dp
792 
793  !Sum of the occupations
794  do i=1,norb
795  ft1 = ft1 + 2.0_dp*fermi(eigenvalues(i),ef,kbt)
796  enddo
797  ft1=ft1-nel
798 
799  err = .false.
800  do m=1,1000001
801 
802  if(m.gt.1000000)then
803  err = .true.
804  write(*,*) "WARNING: Bisection method in prg_get_flevel not converging ..."
805  endif
806 
807  if(abs(ft1).lt.tol)then !tolerance control
808  return
809  endif
810 
811  ef = ef + step
812 
813  ft2=0.0_dp
814 
815  !New sum of the occupations
816  do i=1,norb
817  ft2 = ft2 + 2.0_dp*fermi(eigenvalues(i),ef,kbt)
818  enddo
819 
820  ft2=ft2-nel
821 
822  !Product to see the change in sign.
823  prod = ft2*ft1
824 
825  if(prod.lt.0)then
826  ef=ef-step
827  step=step/2.0_dp !If the root is inside we shorten the step.
828  else
829  ft1=ft2 !If not, Ef moves forward.
830  endif
831 
832  enddo
833 
834  end subroutine prg_get_flevel
835 
847  subroutine prg_get_flevel_nt(eigenvalues,kbt,bndfil,tol,ef,err,verbose)
848 
849  integer :: i, m
850  integer :: norb
851  real(dp) :: ef0, f1, f2, nel
852  real(dp) :: step
853  real(dp), intent(in) :: bndfil, kbt
854  real(dp), intent(in) :: tol, eigenvalues(:)
855  real(dp), intent(inout) :: ef
856  integer, optional, intent(in) :: verbose
857  logical, intent(inout) :: err
858 
859  norb = size(eigenvalues, dim=1)
860  nel = bndfil*2.0_dp*real(norb,dp)
861  ef0 = ef
862  f1 = 0.0_dp
863  f2 = 0.0_dp
864  step = 0.1_dp
865  err = .false.
866 
867  !$omp parallel do default(none) private(i) &
868  !$omp shared(eigenvalues,kbt,ef,norb) &
869  !$omp reduction(+:f1)
870  do i=1,norb
871  f1 = f1 + 2.0_dp*fermi(eigenvalues(i),ef,kbt)
872  enddo
873  !$omp end parallel do
874 
875  f1=f1-nel
876  ef = ef0 + step
877 
878  f2 = 0.0_dp
879 
880  !$omp parallel do default(none) private(i) &
881  !$omp shared(eigenvalues,kbt,ef,norb) &
882  !$omp reduction(+:f2)
883  do i=1,norb
884  f2 = f2 + 2.0_dp*fermi(eigenvalues(i),ef,kbt)
885  enddo
886  !$omp end parallel do
887 
888  f2=f2-nel
889  ef0 = ef
890  if(abs(f2 - f1) < 1.0d-5)then
891  err = .true.
892  return
893  endif
894 
895  ef = -f2*step/(f2-f1) + ef0
896  f1 = f2
897  step = ef - ef0
898 
899  do m = 1,101
900  if(m.gt.100)then
901  write(*,*) "WARNING: Newton method in prg_get_flevel_nt is not converging ..."
902  err = .true.
903  ef = ef0
904  exit
905  endif
906 
907  !New sum of the occupations
908  f2 = 0.0_dp
909  !$omp parallel do default(none) private(i) &
910  !$omp shared(eigenvalues,ef,kbt,norb) &
911  !$omp reduction(+:f2)
912  do i=1,norb
913  f2 = f2 + 2.0_dp*fermi(eigenvalues(i),ef,kbt)
914  enddo
915  !$omp end parallel do
916 
917  f2=f2-nel
918  ef0 = ef
919  ef = -f2*step/(f2-f1) + ef0
920  f1 = f2
921  step = ef - ef0
922  if(abs(f1).lt.tol)then !tolerance control
923  return
924  endif
925  enddo
926 
927  end subroutine prg_get_flevel_nt
928 
929 
935  subroutine prg_get_eigenvalues(ham_bml,eigenvalues,verbose)
936 
937  character(20) :: bml_type
938  integer :: i, norb
939  integer, intent(in) :: verbose
940  real(dp) :: nocc
941  real(dp), allocatable :: aux(:,:)
942  real(dp), allocatable, intent(inout) :: eigenvalues(:)
943  type(bml_matrix_t) :: aux_bml, eigenvectors_bml
944  type(bml_matrix_t), intent(in) :: ham_bml
945 
946  norb = bml_get_n(ham_bml)
947  bml_type = bml_get_deep_type(ham_bml)
948 
949  if(verbose.ge.1)write(*,*)"In prg_get_eigenvalues ..."
950  if(verbose.ge.1)write(*,*)"Number of states =",norb
951 
952  if(.not.allocated(eigenvalues))allocate(eigenvalues(norb))
953 
954  !Ensure dense type to diagonalize
955  if(trim(bml_type).ne."dense")then
956  allocate(aux(norb,norb))
957  call bml_export_to_dense(ham_bml,aux)
958  call bml_zero_matrix(bml_matrix_dense,bml_element_real,dp,norb,norb,aux_bml)
959  call bml_import_from_dense(bml_matrix_dense,aux,aux_bml,0.0_dp,norb)
960  deallocate(aux)
961  else
962  call bml_copy_new(ham_bml,aux_bml)
963  endif
964 
965  call bml_zero_matrix(bml_matrix_dense,bml_element_real,dp,norb,norb,eigenvectors_bml)
966 
967  call bml_diagonalize(aux_bml,eigenvalues,eigenvectors_bml)
968 
969  call bml_deallocate(eigenvectors_bml)
970  call bml_deallocate(aux_bml)
971 
972  end subroutine prg_get_eigenvalues
973 
974 
980  subroutine prg_check_idempotency(mat_bml, threshold, idempotency)
981 
982  character(20) :: bml_type
983  integer :: n, i, j
984  real(dp), intent(in) :: threshold
985  real(dp), intent(out) :: idempotency
986  type(bml_matrix_t) :: aux_bml
987  type(bml_matrix_t), intent(in) :: mat_bml
988 
989  call bml_copy_new(mat_bml, aux_bml)
990 
991  call bml_multiply(mat_bml, mat_bml, aux_bml, 1.0_dp, 0.0_dp, threshold) !A^2
992  call bml_add_deprecated(1.0_dp,aux_bml,-1.0_dp,mat_bml, threshold) !A^2 - A
993 
994  idempotency = bml_fnorm(aux_bml)
995 
996  call bml_deallocate(aux_bml)
997 
998  end subroutine prg_check_idempotency
999 
1004  real(dp) function fermi(e,ef,kbt)
1005 
1006  real(dp), intent(in) :: e, ef, kbt
1007 
1008  if ((e-ef)/kbt > 100.0_dp) then
1009  fermi = 0.0_dp
1010  else
1011  fermi = 1.0_dp/(1.0_dp+exp((e-ef)/(kbt)))
1012  endif
1013 
1014  end function fermi
1015 
1020  real(dp) function dfermi(e,ef,kbt)
1021 
1022  real(dp), intent(in) :: e, ef, kbt
1023 
1024  dfermi = 0.0_dp
1025  if ((e-ef)/kbt > 100.0_dp) then
1026  if (abs(e-ef)<0.001_dp) then
1027  dfermi = - 1000.0_dp
1028  endif
1029  else
1030  dfermi = 1.0_dp/(1.0_dp+exp((e-ef)/(kbt)))
1031  dfermi = dfermi * dfermi * exp((e-ef)/(kbt)) * 1.0_dp/kbt
1032  endif
1033 
1034  end function dfermi
1035 
1046  subroutine prg_toeigenspace(mat_bml,matEig_bml,evects_bml,threshold,verbose)
1047 
1048  integer, optional, intent(in) :: verbose
1049  type(bml_matrix_t), intent(in) :: mat_bml, evects_bml
1050  type(bml_matrix_t), intent(inout) :: mateig_bml
1051  type(bml_matrix_t) :: aux_bml
1052  real(dp), intent(in) :: threshold
1053  integer :: hdim, mdim
1054  character(20) :: bml_type
1055 
1056  if(verbose.eq.1) write(*,*)"In prg_toEigenspace ..."
1057 
1058  hdim= bml_get_n(mat_bml)
1059  mdim= bml_get_m(mat_bml)
1060  bml_type = bml_get_type(mat_bml)
1061  !Allocate bml's
1062  if(bml_get_n(mateig_bml) < 0)then
1063  call bml_zero_matrix(bml_type,bml_element_real,dp,hdim ,mdim,mateig_bml, &
1064  & bml_get_distribution_mode(mat_bml))
1065  endif
1066  !Do the operations in bml
1067  call bml_transpose_new(evects_bml, aux_bml)
1068 
1069  call bml_multiply(aux_bml, mat_bml, mateig_bml, 1.0_dp, 0.0_dp,threshold) !U^t*O
1070 
1071  call bml_multiply(mateig_bml, evects_bml, aux_bml, 1.0_dp, 0.0_dp,threshold) !U^t*O * U
1072 
1073  call bml_copy(aux_bml, mateig_bml)
1074 
1075  call bml_deallocate(aux_bml)
1076 
1077  end subroutine prg_toeigenspace
1078 
1082  !transformations
1083  !! constructed from the eigenvectors of a Hamiltonian
1084  !! \param mat_bml Operator to be transformed
1085  !! \param matEig_bml Output operator after transformation
1086  !! \param evects_bml Eigenvectors matrix (U)
1087  !! \param threshold Threshold parameter
1088  !! \verbose Verbosity level
1089  !!
1090  subroutine prg_tocanonicalspace(mat_bml,matCan_bml,evects_bml,threshold,verbose)
1091 
1092  integer, optional, intent(in) :: verbose
1093  type(bml_matrix_t), intent(in) :: mat_bml, evects_bml
1094  type(bml_matrix_t), intent(inout) :: matcan_bml
1095  type(bml_matrix_t) :: aux_bml
1096  real(dp), intent(in) :: threshold
1097  integer :: hdim, mdim
1098  character(20) :: bml_type
1099 
1100  if(verbose.eq.1) write(*,*)"In prg_toCacninicalspace ..."
1101 
1102  hdim= bml_get_n(mat_bml)
1103  mdim= bml_get_m(mat_bml)
1104  bml_type = bml_get_type(mat_bml)
1105 
1106  !Allocate bml's
1107  if(bml_get_n(matcan_bml) < 0)then
1108  call bml_zero_matrix(bml_type,bml_element_real,dp,hdim ,mdim,matcan_bml,&
1109  & bml_get_distribution_mode(mat_bml))
1110  endif
1111 
1112  !Do the operations in bml
1113 
1114  call bml_transpose_new(evects_bml, aux_bml)
1115  call bml_multiply(mat_bml, aux_bml, matcan_bml, 1.0_dp, 0.0_dp,threshold) !O*U^t
1116 
1117  call bml_multiply(evects_bml,matcan_bml, aux_bml, 1.0_dp, 0.0_dp,threshold) !U*O * U^t
1118 
1119  call bml_copy(aux_bml, matcan_bml)
1120 
1121  call bml_deallocate(aux_bml)
1122 
1123  end subroutine prg_tocanonicalspace
1124 
1125 
1126  subroutine canon_dm_prt(P1,H1,Nocc,T,Q,e,mu0,m,HDIM)
1127 
1128  implicit none
1129  integer, parameter :: prec = 8
1130  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
1131  integer, intent(in) :: hdim, m ! Nocc = Number of occupied orbitals, m
1132  real(prec), intent(in) :: h1(hdim,hdim), q(hdim,hdim), e(hdim), nocc !
1133  real(prec), intent(out) :: p1(hdim,hdim) ! Density matrix response derivative
1134  real(prec), intent(in) :: t, mu0 ! Electronic temperature and chemicalpotential
1135  real(prec) :: x(hdim,hdim), dx1(hdim,hdim), y(hdim,hdim) !Temporary matrices
1136  real(prec) :: h_0(hdim), p_0(hdim), dpdmu(hdim), p_02(hdim),id0(hdim)
1137  real(prec) :: beta, cnst, kb, mu1, sumdpdmu
1138  integer :: i, j, k, norbscore
1139 
1140  kb = 8.61739e-5 ! (eV/K)
1141  beta = 1.d0/(kb*t) ! Temp in Kelvin
1142  h_0 = e ! Diagonal Hamiltonian H0 respresented in the
1143  cnst = beta/(1.d0*2**(m+2)) ! Scaling constant
1144  p_0 = 0.5d0 + cnst*(h_0-mu0) ! Initialization for P0 represented in
1145 
1146  ! call MMult(ONE,Q,H1,ZERO,X,'T','N',HDIM) ! Main cost are the
1147  ! transformation
1148  ! call MMult(ONE,X,Q,ZERO,Y,'N','N',HDIM) ! to the eigenbasis of H1
1149  ! P1 = -cnst*Y !(set mu1 = 0 for simplicity) ! Initialization of DM response
1150  ! (not diagonal in Q)
1151  p1 = -cnst*h1 !(set mu1 = 0 for simplicity) ! Initialization of DM response
1152 
1153 
1154  do i = 1,m ! Loop over m recursion steps
1155  p_02 = p_0*p_0
1156  do j = 1,hdim
1157  do k = 1,hdim
1158  dx1(k,j) = p_0(k)*p1(k,j) + p1(k,j)*p_0(j)
1159  enddo
1160  enddo
1161  id0 = 1.d0/(2.d0*(p_02-p_0)+1.d0)
1162  p_0 = id0*p_02
1163  do j = 1,hdim
1164  do k = 1,hdim
1165  p1(k,j) = id0(k)*(dx1(k,j) + 2.d0*(p1(k,j)-dx1(k,j))*p_0(j))
1166  enddo
1167  enddo
1168  enddo
1169  dpdmu = beta*p_0*(1.d0-p_0)
1170  mu1 = 0.d0
1171  sumdpdmu = 0.d0
1172  !do i = 1,HDIM
1173  do i = 1,norbscore
1174  mu1 = mu1 + p1(i,i)
1175  sumdpdmu = sumdpdmu + dpdmu(i)
1176  enddo
1177  write(*,*)"mu1",mu1
1178 
1179  !mu1 = -mu1/SUM(dPdmu)
1180  mu1 = -mu1/sumdpdmu
1181  do i = 1,hdim
1182  p1(i,i) = p1(i,i) + mu1*dpdmu(i) ! Trace correction by adding (dP/dmu)*(dmu/dH1) to dP/dH1
1183  enddo
1184 
1185  ! call MMult(ONE,Q,P1,ZERO,X,'N','N',HDIM) ! and back transformation of P1,
1186  ! with a total of
1187  ! call MMult(ONE,X,Q,ZERO,P1,'N','T',HDIM) ! 4 matrix-matrix multiplication
1188  ! dominates the cost
1189 
1190  end subroutine canon_dm_prt
1191 
1192  subroutine eig(A,Q,ee,TYPE,HDIM)
1193 
1194  implicit none
1195  integer, parameter :: PREC = 8
1196  integer, intent(in) :: HDIM
1197  integer :: INFO
1198  integer :: DIAG_LWORK
1199  real(PREC) :: DIAG_WORK(1+6*HDIM+2*HDIM*HDIM)
1200  integer :: DIAG_LIWORK
1201  integer :: DIAG_IWORK(3+5*HDIM)
1202  real(PREC), intent(in) :: A(HDIM,HDIM)
1203  real(PREC), intent(out) :: Q(HDIM,HDIM), ee(HDIM)
1204  real(PREC) :: EVECS(HDIM,HDIM), EVALS(HDIM)
1205  character(LEN=1), parameter :: JOBZ = "V", uplo = "U"
1206  character(1), intent(in) :: TYPE ! 'N'(for eigenvaules only) or 'V'(otherwise)
1207 
1208  diag_lwork = 1+6*hdim+2*hdim*hdim
1209  diag_liwork = 3 + 5*hdim
1210 
1211  evecs = a
1212  CALL dsyevd(jobz, uplo, hdim, evecs, hdim, evals, &
1213  diag_work, diag_lwork, diag_iwork, diag_liwork, info)
1214 
1215  q = evecs
1216  ee = evals
1217 
1218  end subroutine eig
1219 
1220 
1221 end module prg_densitymatrix_mod
prg_densitymatrix_mod
Module to obtain the density matrix by diagonalizing an orthogonalized Hamiltonian.
Definition: prg_densitymatrix_mod.F90:5
prg_densitymatrix_mod::dp
integer, parameter dp
Definition: prg_densitymatrix_mod.F90:19
prg_densitymatrix_mod::fermi
real(dp) function fermi(e, ef, kbt)
Gives the Fermi distribution value for energy e.
Definition: prg_densitymatrix_mod.F90:1005
prg_densitymatrix_mod::prg_get_flevel_nt
subroutine, public prg_get_flevel_nt(eigenvalues, kbt, bndfil, tol, ef, err, verbose)
Routine to compute the Fermi level given a set of eigenvalues and a temperature. It applies the Newto...
Definition: prg_densitymatrix_mod.F90:848
prg_densitymatrix_mod::prg_build_density_t
subroutine, public prg_build_density_t(ham_bml, rho_bml, threshold, bndfil, kbt, ef, eigenvalues_out)
Builds the density matrix from for electronic temperature T. Where, is the matrix eigenvector and ...
Definition: prg_densitymatrix_mod.F90:123
prg_densitymatrix_mod::prg_build_density_t_fermi
subroutine, public prg_build_density_t_fermi(ham_bml, rho_bml, threshold, kbt, ef, verbose, drho)
Builds the density matrix from for electronic temperature T. Where, is the matrix eigenvector and ...
Definition: prg_densitymatrix_mod.F90:567
prg_densitymatrix_mod::prg_build_density_t0
subroutine, public prg_build_density_t0(ham_bml, rho_bml, threshold, bndfil, eigenvalues_out)
Builds the density matrix from for zero electronic temperature. Where, is the matrix eigenvector a...
Definition: prg_densitymatrix_mod.F90:42
prg_densitymatrix_mod::prg_build_density_t_ed
subroutine, public prg_build_density_t_ed(ham_bml, rho_bml, evects_bml, threshold, bndfil, kbt, ef, evals, dvals, hindex, llsize, verbose)
Builds the density matrix from for electronic temperature T. Where, is the matrix eigenvector and ...
Definition: prg_densitymatrix_mod.F90:294
prg_densitymatrix_mod::prg_get_eigenvalues
subroutine, public prg_get_eigenvalues(ham_bml, eigenvalues, verbose)
Gets the eigenvalues of the Orthogonalized Hamiltonian.
Definition: prg_densitymatrix_mod.F90:936
prg_densitymatrix_mod::prg_tocanonicalspace
subroutine, public prg_tocanonicalspace(mat_bml, matCan_bml, evects_bml, threshold, verbose)
Change an operator into the eigenspace representation.
Definition: prg_densitymatrix_mod.F90:1091
prg_densitymatrix_mod::prg_check_idempotency
subroutine, public prg_check_idempotency(mat_bml, threshold, idempotency)
To check the idempotency error of a matrix. This is calculated as the Frobenius norm of .
Definition: prg_densitymatrix_mod.F90:981
prg_densitymatrix_mod::prg_build_atomic_density
subroutine, public prg_build_atomic_density(rhoat_bml, numel, hindex, spindex, norb, bml_type)
Builds the atomic density matrix. Where, is the number of electrons for orbital i.
Definition: prg_densitymatrix_mod.F90:651
prg_densitymatrix_mod::dfermi
real(dp) function dfermi(e, ef, kbt)
Gives the gradient of Fermi distribution w.r.t. ef.
Definition: prg_densitymatrix_mod.F90:1021
prg_densitymatrix_mod::prg_build_density_t_fulldata
subroutine, public prg_build_density_t_fulldata(ham_bml, rho_bml, threshold, bndfil, kbt, ef, eigenvalues_out, evects_bml, fvals)
Builds the density matrix from for electronic temperature T. Where, is the matrix eigenvector and ...
Definition: prg_densitymatrix_mod.F90:209
prg_densitymatrix_mod::canon_dm_prt
subroutine, public canon_dm_prt(P1, H1, Nocc, T, Q, e, mu0, m, HDIM)
Definition: prg_densitymatrix_mod.F90:1127
prg_densitymatrix_mod::prg_get_evalsdvalsevects
subroutine, public prg_get_evalsdvalsevects(ham_bml, threshold, hindex, llsize, evals, dvals, evects_bml, verbose)
Gets the eigenvalues and eigenvectors and the core contribution to each eigenvalue (dvals) .
Definition: prg_densitymatrix_mod.F90:388
prg_densitymatrix_mod::prg_toeigenspace
subroutine, public prg_toeigenspace(mat_bml, matEig_bml, evects_bml, threshold, verbose)
Change an operator into the eigenspace representation.
Definition: prg_densitymatrix_mod.F90:1047
prg_densitymatrix_mod::eig
subroutine eig(A, Q, ee, TYPE, HDIM)
Definition: prg_densitymatrix_mod.F90:1193
prg_densitymatrix_mod::prg_build_density_fromevalsandevects
subroutine, public prg_build_density_fromevalsandevects(evects_bml, evals, rho_bml, threshold, bndfil, kbt, ef, verbose)
Builds the density matrix from the evects and evals for electronic temperature T.
Definition: prg_densitymatrix_mod.F90:483
prg_densitymatrix_mod::prg_get_flevel
subroutine, public prg_get_flevel(eigenvalues, kbt, bndfil, tol, Ef, err)
Routine to compute the Fermi level given a set of eigenvalues and a temperature. It applies the Bisec...
Definition: prg_densitymatrix_mod.F90:776
prg_nvtx_mod
Definition: prg_nvtx_mod.F90:1
prg_parallel_mod
The parallel module.
Definition: prg_parallel_mod.F90:6
prg_parallel_mod::printrank
integer function, public printrank()
Definition: prg_parallel_mod.F90:113
prg_parallel_mod::getnranks
integer function, public getnranks()
Definition: prg_parallel_mod.F90:90