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 prg_parallel_mod
9 
10  implicit none
11 
12  private !Everything is private by default
13 
14  integer, parameter :: dp = kind(1.0d0)
15 
16  public :: prg_build_density_t0, prg_check_idempotency, prg_get_eigenvalues
17  public :: prg_get_flevel, prg_build_density_t, prg_build_atomic_density
18  public :: prg_build_density_t_fermi, prg_get_flevel_nt, prg_build_density_t_fulldata
19  public :: prg_build_density_t_ed, prg_toeigenspace, prg_tocanonicalspace
20  public :: canon_dm_prt, prg_get_evalsdvalsevects, prg_build_density_fromevalsandevects
21 
22 contains
23 
36  subroutine prg_build_density_t0(ham_bml, rho_bml, threshold, bndfil, eigenvalues_out)
37 
38  character(20) :: bml_type
39  integer :: i, norb
40  real(8), intent(in) :: bndfil, threshold
41  real(dp) :: nocc
42  real(dp), allocatable :: eigenvalues(:)
43  real(dp), allocatable, optional, intent(out) :: eigenvalues_out(:)
44  type(bml_matrix_t) :: aux1_bml, aux2_bml, eigenvectors_bml
45  type(bml_matrix_t), intent(in) :: ham_bml
46  type(bml_matrix_t), intent(inout) :: rho_bml
47  character(20) :: bml_dmode
48 
49  if (getnranks().gt.1)then
50  if (printrank() .eq. 1) print*,'prg_build_density_T0: BML_DMODE_DISTRIBUTED'
51  bml_dmode = bml_dmode_distributed
52  else
53  print*,'prg_build_density_T0: BML_DMODE_SEQUENTIAL'
54  bml_dmode = bml_dmode_sequential
55  endif
56 
57  norb = bml_get_n(ham_bml)
58  bml_type = bml_get_deep_type(ham_bml)
59 
60  allocate(eigenvalues(norb))
61 
62  call bml_noinit_matrix(bml_type,bml_element_real,dp,norb,norb,eigenvectors_bml,bml_dmode)
63 
64  call bml_diagonalize(ham_bml,eigenvalues,eigenvectors_bml)
65  if(present(eigenvalues_out))then
66  if(allocated(eigenvalues_out))deallocate(eigenvalues_out)
67  allocate(eigenvalues_out(norb))
68  eigenvalues_out = eigenvalues
69  endif
70 
71  nocc = norb*bndfil
72 
73  do i=1,norb !Reusing eigenvalues to apply the theta function.
74  if(real(i)-nocc < 0.0001_dp) then
75  eigenvalues(i) = 2.0_dp
76  elseif(abs(real(i)-real(nocc)) < 0.0001_dp) then
77  eigenvalues(i) = 2.0_dp
78  else
79  eigenvalues(i) = 0.0_dp
80  endif
81  enddo
82  if(abs(nocc - int(nocc)) > 0.01_dp)then
83  eigenvalues(int(nocc)+1) = 1.0_dp
84  endif
85 
86  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml,bml_dmode)
87  call bml_set_diagonal(aux1_bml, eigenvalues) !eps(i,i) = eps(i)
88 
89  deallocate(eigenvalues)
90 
91  call bml_noinit_matrix(bml_type,bml_element_real,dp,norb,norb,aux2_bml,bml_dmode)
92  call bml_multiply(eigenvectors_bml, aux1_bml, aux2_bml, 1.0_dp, 0.0_dp, threshold)
93 
94  call bml_transpose(eigenvectors_bml, aux1_bml)
95  call bml_deallocate(eigenvectors_bml)
96 
97  call bml_multiply(aux2_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
98  call bml_deallocate(aux2_bml)
99  call bml_deallocate(aux1_bml)
100 
101  end subroutine prg_build_density_t0
102 
117  subroutine prg_build_density_t(ham_bml, rho_bml, threshold, bndfil, kbt, ef, eigenvalues_out)
118 
119  character(20) :: bml_type
120  integer :: i, norb
121  real(dp), intent(in) :: bndfil, threshold, kbt
122  real(dp), intent(inout) :: ef
123  real(dp) :: nocc, fleveltol, mlsi, efold
124  real(dp), allocatable :: eigenvalues(:)
125  real(dp), allocatable, optional, intent(out) :: eigenvalues_out(:)
126  type(bml_matrix_t) :: aux1_bml, aux2_bml, eigenvectors_bml
127  type(bml_matrix_t), intent(in) :: ham_bml
128  type(bml_matrix_t), intent(inout) :: rho_bml
129  logical :: err
130 
131  if (printrank() .eq. 1) then
132  write(*,*)"In get_density_t ..."
133  endif
134 
135  norb = bml_get_n(ham_bml)
136  bml_type = bml_get_deep_type(ham_bml)
137 
138  allocate(eigenvalues(norb))
139  call bml_noinit_matrix(bml_type,bml_element_real,dp,norb,norb,eigenvectors_bml)
140 
141  call bml_diagonalize(ham_bml,eigenvalues,eigenvectors_bml)
142 
143  if(present(eigenvalues_out))then
144  if(allocated(eigenvalues_out))deallocate(eigenvalues_out)
145  allocate(eigenvalues_out(norb))
146  eigenvalues_out = eigenvalues
147  endif
148 
149  fleveltol = 1.0e-11
150 
151  efold = ef
152  call prg_get_flevel_nt(eigenvalues,kbt,bndfil,fleveltol,ef,err)
153  if(err)then
154  ef = efold
155  call prg_get_flevel(eigenvalues,kbt,bndfil,fleveltol,ef,err)
156  endif
157  if(err)then
158  write(*,*)"WARNING: Ef/Chemical potential search failed. We'll use the previous one to proceed"
159  ef = efold
160  endif
161 
162  nocc = norb*bndfil
163 
164  do i=1,norb !Reusing eigenvalues to apply the theta function.
165  eigenvalues(i) = 2.0_dp*fermi(eigenvalues(i),ef,kbt)
166  enddo
167 
168  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml)
169  call bml_set_diagonal(aux1_bml, eigenvalues) !eps(i,i) = eps(i)
170 
171  deallocate(eigenvalues)
172 
173  call bml_noinit_matrix(bml_type,bml_element_real,dp,norb,norb,aux2_bml)
174  call bml_multiply(eigenvectors_bml, aux1_bml, aux2_bml, 1.0_dp, 0.0_dp, threshold)
175 
176  call bml_transpose(eigenvectors_bml, aux1_bml)
177  call bml_deallocate(eigenvectors_bml)
178 
179  call bml_multiply(aux2_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
180  call bml_deallocate(aux2_bml)
181  call bml_deallocate(aux1_bml)
182 
183  end subroutine prg_build_density_t
184 
185 
202  subroutine prg_build_density_t_fulldata(ham_bml, rho_bml, threshold, bndfil, kbt, ef, eigenvalues_out&
203  &, evects_bml, fvals)
204 
205  character(20) :: bml_type
206  integer :: i, norb
207  real(dp), intent(in) :: bndfil, threshold, kbt
208  real(dp), intent(inout) :: ef
209  real(dp) :: nocc, fleveltol, efold
210  real(dp), allocatable :: eigenvalues(:)
211  real(dp), allocatable, intent(inout) :: eigenvalues_out(:), fvals(:)
212  type(bml_matrix_t) :: aux1_bml, aux_bml, occupation_bml
213  type(bml_matrix_t), intent(in) :: ham_bml
214  type(bml_matrix_t), intent(inout) :: rho_bml, evects_bml
215  logical :: err
216 
217  if (printrank() .eq. 1) then
218  write(*,*)"In get_density_t_fulldata ..."
219  endif
220 
221  norb = bml_get_n(ham_bml)
222  bml_type = bml_get_deep_type(ham_bml)
223 
224  allocate(eigenvalues(norb))
225  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,evects_bml)
226 
227  call bml_diagonalize(ham_bml,eigenvalues,evects_bml)
228 
229  if(.not.allocated(eigenvalues_out))then
230  allocate(eigenvalues_out(norb))
231  allocate(fvals(norb))
232  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,evects_bml)
233  endif
234 
235  eigenvalues_out = eigenvalues
236 
237  fleveltol = 1.0e-11
238  fvals = 0.0_dp
239 
240  call prg_get_flevel_nt(eigenvalues,kbt,bndfil,fleveltol,ef,err)
241  if(err)call prg_get_flevel(eigenvalues,kbt,bndfil,fleveltol,ef,err)
242  if(err)then
243  write(*,*)"WARNING: Ef/Chemical potential search failed. We'll use the previous one to proceed"
244  ef = efold
245  endif
246 
247  nocc = norb*bndfil
248 
249  do i=1,norb !Aapply Fermi function.
250  fvals(i) = 2.0_dp*fermi(eigenvalues(i),ef,kbt)
251  enddo
252 
253  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,occupation_bml)
254  call bml_set_diagonal(occupation_bml, fvals) !eps(i,i) = eps(i)
255 
256  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux_bml)
257  call bml_multiply(evects_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp, threshold)
258  call bml_deallocate(occupation_bml)
259 
260  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml)
261  call bml_transpose(evects_bml, aux1_bml)
262 
263  call bml_multiply(aux_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
264  call bml_deallocate(aux_bml)
265  call bml_deallocate(aux1_bml)
266 
267  end subroutine prg_build_density_t_fulldata
268 
269 
287  subroutine prg_build_density_t_ed(ham_bml, rho_bml, evects_bml, threshold, bndfil, kbt, ef, evals&
288  &, dvals, hindex, llsize, verbose)
289 
290  character(20) :: bml_type
291  integer, allocatable, intent(in) :: hindex(:,:)
292  integer :: i, j, jj, k, kk, l, norb, norbcore
293  integer, intent(in) :: llsize, verbose
294  real(dp), intent(in) :: bndfil, threshold, kbt
295  real(dp), intent(inout) :: ef
296  real(dp) :: nocc
297  real(dp), allocatable :: eigenvalues(:), row(:),fvals(:),ham(:,:),evects(:,:)
298  real(dp), allocatable, intent(inout) :: evals(:), dvals(:)
299  type(bml_matrix_t) :: aux1_bml, aux_bml, occupation_bml
300  type(bml_matrix_t), intent(in) :: ham_bml
301  type(bml_matrix_t), intent(inout) :: rho_bml, evects_bml
302 
303  if (printrank() .eq. 1 .and. verbose >= 1) then
304  write(*,*)"In get_density_t_efd ..."
305  endif
306 
307  norb = bml_get_n(ham_bml)
308  bml_type = bml_get_type(ham_bml)
309 
310  if(allocated(evals))then
311  deallocate(evals)
312  deallocate(dvals)
313  endif
314  allocate(evals(norb))
315  allocate(dvals(norb))
316  allocate(fvals(norb))
317  allocate(ham(norb,norb))
318  allocate(evects(norb,norb))
319  call bml_export_to_dense(ham_bml,ham)
320  if(bml_get_n(evects_bml) < 0)then
321  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,evects_bml)
322  endif
323  call eig(ham,evects,evals,'V',norb)
324  call bml_import_from_dense(bml_type,evects,evects_bml,0.0_dp,norb)
325  !call bml_diagonalize(ham_bml,evals,evects_bml)
326  write(*,*)"ham",ham(1,1),ham(1,2),ham(1,3)
327  deallocate(evects)
328  deallocate(ham)
329  call bml_print_matrix("evects",evects_bml,0,4,0,4)
330  write(*,*)"Evals",evals(1),evals(2),evals(3)
331  nocc = norb*bndfil
332 
333  do i=1,norb !Apply Fermi function.
334  fvals(i) = 2.0_dp*fermi(evals(i),ef,kbt)
335  enddo
336 
337  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,occupation_bml)
338  call bml_set_diagonal(occupation_bml, fvals) !eps(i,i) = eps(i)
339 
340  deallocate(fvals)
341 
342  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux_bml)
343  call bml_multiply(evects_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp, threshold)
344  call bml_deallocate(occupation_bml)
345 
346  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml)
347  call bml_transpose(evects_bml, aux1_bml)
348 
349  allocate(row(norb))
350  dvals = 0.0_dp
351  do i = 1,norb
352  call bml_get_row(aux1_bml,i,row)
353  do k = 1,llsize
354  do l = hindex(1,k),hindex(2,k)
355  dvals(i) = dvals(i) + row(l)**2
356  enddo
357  enddo
358  enddo
359  deallocate(row)
360 
361  call bml_multiply(aux_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
362 
363  call bml_deallocate(aux_bml)
364  call bml_deallocate(aux1_bml)
365 
366  end subroutine prg_build_density_t_ed
367 
368 
381  subroutine prg_get_evalsdvalsevects(ham_bml, threshold, hindex,&
382  & llsize, evals, dvals, evects_bml, verbose)
383 
384  character(20) :: bml_type
385  integer, allocatable, intent(in) :: hindex(:,:)
386  integer :: i, k, l, norb
387  integer, intent(in) :: llsize, verbose
388  real(dp), intent(in) :: threshold
389  real(dp), allocatable :: row(:),ham(:,:),evects(:,:),aux(:,:)
390  real(dp), allocatable, intent(inout) :: evals(:), dvals(:)
391  type(bml_matrix_t), intent(in) :: ham_bml
392  type(bml_matrix_t), intent(inout) :: evects_bml
393  type(bml_matrix_t) :: aux1_bml
394 
395  if (printrank() .eq. 1 .and. verbose >= 1) then
396  write(*,*)"In get_evalsDvalsEvects ..."
397  endif
398 
399  norb = bml_get_n(ham_bml)
400  bml_type = bml_get_type(ham_bml)
401 
402  if(.not.allocated(evals))allocate(evals(norb))
403  if(.not.allocated(dvals))allocate(dvals(norb))
404  !allocate(ham(norb,norb))
405  !allocate(evects(norb,norb))
406  !call bml_export_to_dense(ham_bml,ham)
407  if(.not. bml_allocated(evects_bml))then
408  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,evects_bml)
409  endif
410  !call Eig(ham,evects,evals,'V',norb)
411  !call bml_import_from_dense(bml_type,evects,evects_bml,0.0_dp,norb)
412  call bml_diagonalize(ham_bml,evals,evects_bml)
413  !deallocate(evects)
414  !deallocate(ham)
415 
416  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml)
417  call bml_transpose(evects_bml, aux1_bml)
418  allocate(aux(norb,norb))
419  call bml_export_to_dense(aux1_bml,aux)
420  call bml_deallocate(aux1_bml)
421 
422  allocate(row(norb))
423  dvals = 0.0_dp
424  do i = 1,norb
425  !call bml_get_row(aux1_bml,i,row)
426  row(:) = aux(i,:)
427  do k = 1,llsize
428  do l = hindex(1,k),hindex(2,k)
429  dvals(i) = dvals(i) + row(l)**2
430  enddo
431  enddo
432  enddo
433  deallocate(row)
434  deallocate(aux)
435 
436  end subroutine prg_get_evalsdvalsevects
437 
438 
448  subroutine prg_build_density_fromevalsandevects(evects_bml, evals, rho_bml, threshold,&
449  & bndfil, kbt, ef, verbose)
450 
451  character(20) :: bml_type
452  integer :: i, norb
453  integer, intent(in) :: verbose
454  real(dp), intent(in) :: bndfil, threshold, kbt
455  real(dp), intent(inout) :: ef
456  real(dp) :: nocc
457  real(dp), intent(in) :: evals(*)
458  real(dp), allocatable :: fvals(:)
459  type(bml_matrix_t) :: aux1_bml, aux_bml, occupation_bml
460  type(bml_matrix_t), intent(in) :: evects_bml
461  type(bml_matrix_t), intent(inout) :: rho_bml
462 
463  if (printrank() .eq. 1 .and. verbose >= 1) then
464  write(*,*)"In get_density_fromEvalsAndEvects ..."
465  endif
466 
467  norb = bml_get_n(evects_bml)
468  bml_type = bml_get_type(evects_bml)
469 
470  allocate(fvals(norb))
471 
472  nocc = norb*bndfil
473 
474  do i=1,norb !Apply Fermi function.
475  fvals(i) = 2.0_dp*fermi(evals(i),ef,kbt)
476  enddo
477 
478  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,occupation_bml)
479  call bml_set_diagonal(occupation_bml, fvals) !eps(i,i) = eps(i)
480 
481  deallocate(fvals)
482 
483  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux_bml)
484  call bml_multiply(evects_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp,threshold)
485  call bml_deallocate(occupation_bml)
486 
487  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml)
488  call bml_transpose(evects_bml, aux1_bml)
489 
490  call bml_multiply(aux_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
491 
492  call bml_deallocate(aux_bml)
493  call bml_deallocate(aux1_bml)
494 
495  end subroutine prg_build_density_fromevalsandevects
496 
497 
498 
509  subroutine prg_build_density_t_fermi(ham_bml, rho_bml, threshold, kbt, ef, verbose, drho)
510 
511  character(20) :: bml_type
512  integer :: i, norb, mdim
513  integer, optional, intent(in) :: verbose
514  real(dp), intent(in) :: threshold, kbt
515  real(dp), intent(in) :: ef
516  real(dp) :: fleveltol
517  real(dp), allocatable :: eigenvalues(:)
518  real(dp), allocatable :: nocc(:)
519  type(bml_matrix_t) :: aux1_bml, aux_bml, eigenvectors_bml, occupation_bml
520  type(bml_matrix_t), intent(in) :: ham_bml
521  type(bml_matrix_t), intent(inout) :: rho_bml
522  real(dp), optional, intent(inout) :: drho ! gradient of density w.r.t. ef
523 
524  if (printrank() .eq. 1) then
525  if(present(verbose))then
526  if(verbose >= 1)then
527  write(*,*)"In get_density_t_Fermi ..."
528  write(*,*)"Ef = ",ef
529  endif
530  endif
531  endif
532 
533  norb = bml_get_n(ham_bml)
534  mdim = bml_get_m(ham_bml)
535  bml_type = bml_get_deep_type(ham_bml)
536 
537  allocate(eigenvalues(norb),nocc(norb))
538  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,eigenvectors_bml)
539 
540  call bml_diagonalize(ham_bml,eigenvalues,eigenvectors_bml)
541 
542  do i=1,norb !Reusing eigenvalues to apply the theta function.
543  nocc(i) = 2.0_dp*fermi(eigenvalues(i),ef,kbt)
544  enddo
545 
546  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,occupation_bml)
547  call bml_set_diagonal(occupation_bml, nocc) !eps(i,i) = eps(i)
548 
549  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux_bml)
550  call bml_multiply(eigenvectors_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp, threshold)
551  call bml_deallocate(occupation_bml)
552 
553  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux1_bml)
554  call bml_transpose(eigenvectors_bml, aux1_bml)
555 
556  call bml_multiply(aux_bml, aux1_bml, rho_bml, 1.0_dp, 0.0_dp, threshold)
557 
558  if (present(drho)) then
559  do i=1,norb !Reusing eigenvalues to apply the theta function.
560  nocc(i) = 2.0_dp * dfermi(eigenvalues(i),ef,kbt)
561  enddo
562 
563  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,occupation_bml)
564  call bml_set_diagonal(occupation_bml, nocc) !eps(i,i) = eps(i)
565 
566  call bml_multiply(eigenvectors_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp, threshold)
567 
568  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux1_bml)
569  call bml_transpose(eigenvectors_bml, aux1_bml)
570 
571  call bml_multiply(aux_bml, aux1_bml, occupation_bml, 1.0_dp, 0.0_dp, threshold)
572  drho = bml_trace(occupation_bml)
573  call bml_deallocate(occupation_bml)
574  endif
575 
576  deallocate(nocc)
577  deallocate(eigenvalues)
578  call bml_deallocate(eigenvectors_bml)
579  call bml_deallocate(aux_bml)
580  call bml_deallocate(aux1_bml)
581 
582  end subroutine prg_build_density_t_fermi
583 
593  subroutine prg_build_atomic_density(rhoat_bml,numel,hindex,spindex,norb,bml_type)
594 
595  character(len=*), intent(in) :: bml_type
596  integer :: i, index, n_orb, nats
597  integer, intent(in) :: hindex(:,:), norb, spindex(:)
598  real(dp) :: occ
599  real(dp), allocatable :: d_atomic(:), rhoat(:)
600  real(dp), intent(in) :: numel(:)
601  type(bml_matrix_t), intent(inout) :: rhoat_bml
602 
603  nats = size(hindex,dim=2)
604 
605  allocate(rhoat(norb))
606 
607  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,rhoat_bml)
608 
609  index = 0;
610 
611  do i = 1,nats
612  n_orb = hindex(2,i)-hindex(1,i) + 1;
613  if(n_orb == 1)then
614  index = index + 1;
615  rhoat(index) = numel(spindex(i));
616  else
617  if(numel(spindex(i)) <= 2)then
618  index = index + 1;
619  rhoat(index) = numel(spindex(i));
620  index = index + 1;
621  rhoat(index) = 0.0_dp;
622  index = index + 1;
623  rhoat(index) = 0.0_dp;
624  index = index + 1;
625  rhoat(index) = 0.0_dp;
626  else
627  index = index + 1;
628  rhoat(index) = 2.0_dp;
629 
630  index = index + 1;
631  occ = (numel(spindex(i))-2.0_dp)/3.0_dp;
632  rhoat(index) = occ;
633  index = index + 1;
634  rhoat(index) = occ;
635  index = index + 1;
636  rhoat(index) = occ;
637  endif
638  endif
639  enddo
640 
641  call bml_set_diagonal(rhoat_bml,rhoat,0.0_dp)
642 
643  deallocate(rhoat)
644 
645  end subroutine prg_build_atomic_density
646 
657  subroutine prg_get_flevel(eigenvalues,kbt,bndfil,tol,Ef,err)
658 
659  integer :: i, j, k, m
660  integer :: norb
661  real(dp) :: ft1, ft2, kb, prod
662  real(dp) :: t, step, tol, nel
663  real(dp), intent(in) :: bndfil, eigenvalues(:), kbt
664  real(dp), intent(inout) :: ef
665  logical, intent(inout) :: err
666 
667  norb = size(eigenvalues,dim=1)
668  nel = bndfil*2.0_dp*norb
669  ef=eigenvalues(1)
670  step=abs((eigenvalues(norb)-eigenvalues(1)))
671  ft1=0.0_dp
672  ft2=0.0_dp
673  prod=0.0_dp
674 
675  !Sum of the occupations
676  do i=1,norb
677  ft1 = ft1 + 2.0_dp*fermi(eigenvalues(i),ef,kbt)
678  enddo
679  ft1=ft1-nel
680 
681  err = .false.
682  do m=1,1000001
683 
684  if(m.gt.1000000)then
685  err = .true.
686  write(*,*) "WARNING: Bisection method in prg_get_flevel not converging ..."
687  endif
688 
689  if(abs(ft1).lt.tol)then !tolerance control
690  return
691  endif
692 
693  ef = ef + step
694 
695  ft2=0.0_dp
696 
697  !New sum of the occupations
698  do i=1,norb
699  ft2 = ft2 + 2.0_dp*fermi(eigenvalues(i),ef,kbt)
700  enddo
701 
702  ft2=ft2-nel
703 
704  !Product to see the change in sign.
705  prod = ft2*ft1
706 
707  if(prod.lt.0)then
708  ef=ef-step
709  step=step/2.0_dp !If the root is inside we shorten the step.
710  else
711  ft1=ft2 !If not, Ef moves forward.
712  endif
713 
714  enddo
715 
716  end subroutine prg_get_flevel
717 
729  subroutine prg_get_flevel_nt(eigenvalues,kbt,bndfil,tol,ef,err,verbose)
730 
731  integer :: i, m
732  integer :: norb
733  real(dp) :: ef0, f1, f2, nel
734  real(dp) :: step
735  real(dp), intent(in) :: bndfil, kbt
736  real(dp), intent(in) :: tol, eigenvalues(:)
737  real(dp), intent(inout) :: ef
738  integer, optional, intent(in) :: verbose
739  logical, intent(inout) :: err
740 
741  norb = size(eigenvalues, dim=1)
742  nel = bndfil*2.0_dp*real(norb,dp)
743  ef0 = ef
744  f1 = 0.0_dp
745  f2 = 0.0_dp
746  step = 0.1_dp
747  err = .false.
748 
749  !$omp parallel do default(none) private(i) &
750  !$omp shared(eigenvalues,kbt,ef,norb) &
751  !$omp reduction(+:f1)
752  do i=1,norb
753  f1 = f1 + 2.0_dp*fermi(eigenvalues(i),ef,kbt)
754  enddo
755  !$omp end parallel do
756 
757  f1=f1-nel
758  ef = ef0 + step
759 
760  f2 = 0.0_dp
761 
762  !$omp parallel do default(none) private(i) &
763  !$omp shared(eigenvalues,kbt,ef,norb) &
764  !$omp reduction(+:f2)
765  do i=1,norb
766  f2 = f2 + 2.0_dp*fermi(eigenvalues(i),ef,kbt)
767  enddo
768  !$omp end parallel do
769 
770  f2=f2-nel
771  ef0 = ef
772  if(abs(f2 - f1) < 1.0d-5)then
773  err = .true.
774  return
775  endif
776 
777  ef = -f2*step/(f2-f1) + ef0
778  f1 = f2
779  step = ef - ef0
780 
781  do m = 1,101
782  if(m.gt.100)then
783  write(*,*) "WARNING: Newton method in prg_get_flevel_nt is not converging ..."
784  err = .true.
785  ef = ef0
786  exit
787  endif
788 
789  !New sum of the occupations
790  f2 = 0.0_dp
791  !$omp parallel do default(none) private(i) &
792  !$omp shared(eigenvalues,ef,kbt,norb) &
793  !$omp reduction(+:f2)
794  do i=1,norb
795  f2 = f2 + 2.0_dp*fermi(eigenvalues(i),ef,kbt)
796  enddo
797  !$omp end parallel do
798 
799  f2=f2-nel
800  ef0 = ef
801  ef = -f2*step/(f2-f1) + ef0
802  f1 = f2
803  step = ef - ef0
804  if(abs(f1).lt.tol)then !tolerance control
805  return
806  endif
807  enddo
808 
809  end subroutine prg_get_flevel_nt
810 
811 
817  subroutine prg_get_eigenvalues(ham_bml,eigenvalues,verbose)
818 
819  character(20) :: bml_type
820  integer :: i, norb
821  integer, intent(in) :: verbose
822  real(dp) :: nocc
823  real(dp), allocatable :: aux(:,:)
824  real(dp), allocatable, intent(inout) :: eigenvalues(:)
825  type(bml_matrix_t) :: aux_bml, eigenvectors_bml
826  type(bml_matrix_t), intent(in) :: ham_bml
827 
828  norb = bml_get_n(ham_bml)
829  bml_type = bml_get_deep_type(ham_bml)
830 
831  if(verbose.ge.1)write(*,*)"In prg_get_eigenvalues ..."
832  if(verbose.ge.1)write(*,*)"Number of states =",norb
833 
834  if(.not.allocated(eigenvalues))allocate(eigenvalues(norb))
835 
836  !Ensure dense type to diagonalize
837  if(trim(bml_type).ne."dense")then
838  allocate(aux(norb,norb))
839  call bml_export_to_dense(ham_bml,aux)
840  call bml_zero_matrix(bml_matrix_dense,bml_element_real,dp,norb,norb,aux_bml)
841  call bml_import_from_dense(bml_matrix_dense,aux,aux_bml,0.0_dp,norb)
842  deallocate(aux)
843  else
844  call bml_copy_new(ham_bml,aux_bml)
845  endif
846 
847  call bml_zero_matrix(bml_matrix_dense,bml_element_real,dp,norb,norb,eigenvectors_bml)
848 
849  call bml_diagonalize(aux_bml,eigenvalues,eigenvectors_bml)
850 
851  call bml_deallocate(eigenvectors_bml)
852  call bml_deallocate(aux_bml)
853 
854  end subroutine prg_get_eigenvalues
855 
856 
862  subroutine prg_check_idempotency(mat_bml, threshold, idempotency)
863 
864  character(20) :: bml_type
865  integer :: n, i, j
866  real(dp), intent(in) :: threshold
867  real(dp), intent(out) :: idempotency
868  type(bml_matrix_t) :: aux_bml
869  type(bml_matrix_t), intent(in) :: mat_bml
870 
871  call bml_copy_new(mat_bml, aux_bml)
872 
873  call bml_multiply(mat_bml, mat_bml, aux_bml, 1.0_dp, 0.0_dp, threshold) !A^2
874  call bml_add_deprecated(1.0_dp,aux_bml,-1.0_dp,mat_bml, threshold) !A^2 - A
875 
876  idempotency = bml_fnorm(aux_bml)
877 
878  call bml_deallocate(aux_bml)
879 
880  end subroutine prg_check_idempotency
881 
886  real(dp) function fermi(e,ef,kbt)
887 
888  real(dp), intent(in) :: e, ef, kbt
889 
890  if ((e-ef)/kbt > 100.0_dp) then
891  fermi = 0.0_dp
892  else
893  fermi = 1.0_dp/(1.0_dp+exp((e-ef)/(kbt)))
894  endif
895 
896  end function fermi
897 
902  real(dp) function dfermi(e,ef,kbt)
903 
904  real(dp), intent(in) :: e, ef, kbt
905 
906  dfermi = 0.0_dp
907  if ((e-ef)/kbt > 100.0_dp) then
908  if (abs(e-ef)<0.001_dp) then
909  dfermi = - 1000.0_dp
910  endif
911  else
912  dfermi = 1.0_dp/(1.0_dp+exp((e-ef)/(kbt)))
913  dfermi = dfermi * dfermi * exp((e-ef)/(kbt)) * 1.0_dp/kbt
914  endif
915 
916  end function dfermi
917 
928  subroutine prg_toeigenspace(mat_bml,matEig_bml,evects_bml,threshold,verbose)
929 
930  integer, optional, intent(in) :: verbose
931  type(bml_matrix_t), intent(in) :: mat_bml, evects_bml
932  type(bml_matrix_t), intent(inout) :: mateig_bml
933  type(bml_matrix_t) :: aux_bml
934  real(dp), intent(in) :: threshold
935  integer :: hdim, mdim
936  character(20) :: bml_type
937 
938  if(verbose.eq.1) write(*,*)"In prg_toEigenspace ..."
939 
940  hdim= bml_get_n(mat_bml)
941  mdim= bml_get_m(mat_bml)
942  bml_type = bml_get_type(mat_bml)
943  !Allocate bml's
944  if(bml_get_n(mateig_bml) < 0)then
945  call bml_zero_matrix(bml_type,bml_element_real,dp,hdim ,mdim,mateig_bml, &
946  & bml_get_distribution_mode(mat_bml))
947  endif
948  !Do the operations in bml
949  call bml_transpose(evects_bml, aux_bml)
950 
951  call bml_multiply(aux_bml, mat_bml, mateig_bml, 1.0_dp, 0.0_dp,threshold) !U^t*O
952 
953  call bml_multiply(mateig_bml, evects_bml, aux_bml, 1.0_dp, 0.0_dp,threshold) !U^t*O * U
954 
955  call bml_copy(aux_bml, mateig_bml)
956 
957  call bml_deallocate(aux_bml)
958 
959  end subroutine prg_toeigenspace
960 
964  !transformations
965  !! constructed from the eigenvectors of a Hamiltonian
966  !! \param mat_bml Operator to be transformed
967  !! \param matEig_bml Output operator after transformation
968  !! \param evects_bml Eigenvectors matrix (U)
969  !! \param threshold Threshold parameter
970  !! \verbose Verbosity level
971  !!
972  subroutine prg_tocanonicalspace(mat_bml,matCan_bml,evects_bml,threshold,verbose)
973 
974  integer, optional, intent(in) :: verbose
975  type(bml_matrix_t), intent(in) :: mat_bml, evects_bml
976  type(bml_matrix_t), intent(inout) :: matcan_bml
977  type(bml_matrix_t) :: aux_bml
978  real(dp), intent(in) :: threshold
979  integer :: hdim, mdim
980  character(20) :: bml_type
981 
982  if(verbose.eq.1) write(*,*)"In prg_toCacninicalspace ..."
983 
984  hdim= bml_get_n(mat_bml)
985  mdim= bml_get_m(mat_bml)
986  bml_type = bml_get_type(mat_bml)
987 
988  !Allocate bml's
989  if(bml_get_n(matcan_bml) < 0)then
990  call bml_zero_matrix(bml_type,bml_element_real,dp,hdim ,mdim,matcan_bml,&
991  & bml_get_distribution_mode(mat_bml))
992  endif
993 
994  !Do the operations in bml
995 
996  call bml_transpose(evects_bml, aux_bml)
997  call bml_multiply(mat_bml, aux_bml, matcan_bml, 1.0_dp, 0.0_dp,threshold) !O*U^t
998 
999  call bml_multiply(evects_bml,matcan_bml, aux_bml, 1.0_dp, 0.0_dp,threshold) !U*O * U^t
1000 
1001  call bml_copy(aux_bml, matcan_bml)
1002 
1003  call bml_deallocate(aux_bml)
1004 
1005  end subroutine prg_tocanonicalspace
1006 
1007 
1008  subroutine canon_dm_prt(P1,H1,Nocc,T,Q,e,mu0,m,HDIM)
1009 
1010  implicit none
1011  integer, parameter :: prec = 8
1012  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
1013  integer, intent(in) :: hdim, m ! Nocc = Number of occupied orbitals, m
1014  real(prec), intent(in) :: h1(hdim,hdim), q(hdim,hdim), e(hdim), nocc !
1015  real(prec), intent(out) :: p1(hdim,hdim) ! Density matrix response derivative
1016  real(prec), intent(in) :: t, mu0 ! Electronic temperature and chemicalpotential
1017  real(prec) :: x(hdim,hdim), dx1(hdim,hdim), y(hdim,hdim) !Temporary matrices
1018  real(prec) :: h_0(hdim), p_0(hdim), dpdmu(hdim), p_02(hdim),id0(hdim)
1019  real(prec) :: beta, cnst, kb, mu1, sumdpdmu
1020  integer :: i, j, k, norbscore
1021 
1022  kb = 8.61739e-5 ! (eV/K)
1023  beta = 1.d0/(kb*t) ! Temp in Kelvin
1024  h_0 = e ! Diagonal Hamiltonian H0 respresented in the
1025  cnst = beta/(1.d0*2**(m+2)) ! Scaling constant
1026  p_0 = 0.5d0 + cnst*(h_0-mu0) ! Initialization for P0 represented in
1027 
1028  ! call MMult(ONE,Q,H1,ZERO,X,'T','N',HDIM) ! Main cost are the
1029  ! transformation
1030  ! call MMult(ONE,X,Q,ZERO,Y,'N','N',HDIM) ! to the eigenbasis of H1
1031  ! P1 = -cnst*Y !(set mu1 = 0 for simplicity) ! Initialization of DM response
1032  ! (not diagonal in Q)
1033  p1 = -cnst*h1 !(set mu1 = 0 for simplicity) ! Initialization of DM response
1034 
1035 
1036  do i = 1,m ! Loop over m recursion steps
1037  p_02 = p_0*p_0
1038  do j = 1,hdim
1039  do k = 1,hdim
1040  dx1(k,j) = p_0(k)*p1(k,j) + p1(k,j)*p_0(j)
1041  enddo
1042  enddo
1043  id0 = 1.d0/(2.d0*(p_02-p_0)+1.d0)
1044  p_0 = id0*p_02
1045  do j = 1,hdim
1046  do k = 1,hdim
1047  p1(k,j) = id0(k)*(dx1(k,j) + 2.d0*(p1(k,j)-dx1(k,j))*p_0(j))
1048  enddo
1049  enddo
1050  enddo
1051  dpdmu = beta*p_0*(1.d0-p_0)
1052  mu1 = 0.d0
1053  sumdpdmu = 0.d0
1054  !do i = 1,HDIM
1055  do i = 1,norbscore
1056  mu1 = mu1 + p1(i,i)
1057  sumdpdmu = sumdpdmu + dpdmu(i)
1058  enddo
1059  write(*,*)"mu1",mu1
1060 
1061  !mu1 = -mu1/SUM(dPdmu)
1062  mu1 = -mu1/sumdpdmu
1063  do i = 1,hdim
1064  p1(i,i) = p1(i,i) + mu1*dpdmu(i) ! Trace correction by adding (dP/dmu)*(dmu/dH1) to dP/dH1
1065  enddo
1066 
1067  ! call MMult(ONE,Q,P1,ZERO,X,'N','N',HDIM) ! and back transformation of P1,
1068  ! with a total of
1069  ! call MMult(ONE,X,Q,ZERO,P1,'N','T',HDIM) ! 4 matrix-matrix multiplication
1070  ! dominates the cost
1071 
1072  end subroutine canon_dm_prt
1073 
1074  subroutine eig(A,Q,ee,TYPE,HDIM)
1075 
1076  implicit none
1077  integer, parameter :: PREC = 8
1078  integer, intent(in) :: HDIM
1079  integer :: INFO
1080  integer :: DIAG_LWORK
1081  real(PREC) :: DIAG_WORK(1+6*HDIM+2*HDIM*HDIM)
1082  integer :: DIAG_LIWORK
1083  integer :: DIAG_IWORK(3+5*HDIM)
1084  real(PREC), intent(in) :: A(HDIM,HDIM)
1085  real(PREC), intent(out) :: Q(HDIM,HDIM), ee(HDIM)
1086  real(PREC) :: EVECS(HDIM,HDIM), EVALS(HDIM)
1087  character(LEN=1), parameter :: JOBZ = "V", uplo = "U"
1088  character(1), intent(in) :: TYPE ! 'N'(for eigenvaules only) or 'V'(otherwise)
1089 
1090  diag_lwork = 1+6*hdim+2*hdim*hdim
1091  diag_liwork = 3 + 5*hdim
1092 
1093  evecs = a
1094  CALL dsyevd(jobz, uplo, hdim, evecs, hdim, evals, &
1095  diag_work, diag_lwork, diag_iwork, diag_liwork, info)
1096 
1097  q = evecs
1098  ee = evals
1099 
1100  end subroutine eig
1101 
1102 
1103 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:14
prg_densitymatrix_mod::fermi
real(dp) function fermi(e, ef, kbt)
Gives the Fermi distribution value for energy e.
Definition: prg_densitymatrix_mod.F90:887
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:730
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:118
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:510
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:37
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:289
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:818
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:973
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:863
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:594
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:903
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:204
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:1009
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:383
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:929
prg_densitymatrix_mod::eig
subroutine eig(A, Q, ee, TYPE, HDIM)
Definition: prg_densitymatrix_mod.F90:1075
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:450
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:658
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:111
prg_parallel_mod::getnranks
integer function, public getnranks()
Definition: prg_parallel_mod.F90:88