PROGRESS  master
prg_implicit_fermi_mod.F90
Go to the documentation of this file.
1 ! The Implicit Recursive Fermi O(N) module.
2 !! \ingroup PROGRESS
3 !! \brief Here are subroutines implementing Niklasson's implicit recursive fermi dirac exact
4 !! density matrix purification algorithm.
5 !!
6 module prg_implicit_fermi_mod
7 
8  use omp_lib
9  use bml
10  use prg_normalize_mod
11  use prg_densitymatrix_mod
12  use prg_timer_mod
13  use prg_parallel_mod
14  use prg_ewald_mod
15 
16  implicit none
17 
18  private !Everything is private by default
19 
20  integer, parameter :: dp = kind(1.0d0)
21 
22  public :: prg_implicit_fermi
23  public :: prg_implicit_fermi_save_inverse
24  public :: prg_implicit_fermi_zero
25  public :: prg_test_density_matrix
26  public :: prg_implicit_fermi_response
27  public :: prg_implicit_fermi_first_order_response
28  public :: prg_finite_diff
29 
30 contains
31 
46  subroutine prg_implicit_fermi_save_inverse(Inv_bml, h_bml, p_bml, nsteps, nocc, &
47  mu, beta, occErrLimit, threshold, tol,SCF_IT, occiter, totns)
48 
49  implicit none
50 
51  type(bml_matrix_t), intent(in) :: h_bml
52  type(bml_matrix_t), intent(inout) :: p_bml, inv_bml(nsteps)
53  integer, intent(in) :: nsteps, scf_it
54  real(dp), intent(in) :: nocc, threshold
55  real(dp), intent(in) :: tol
56  real(dp), intent(in) :: occerrlimit, beta
57  real(dp), intent(inout) :: mu
58  integer, intent(inout) :: occiter, totns
59 
60  type(bml_matrix_t) :: w_bml, y_bml, d_bml, aux_bml, p2_bml, i_bml, ai_bml
61  real(dp) :: trdpdmu, trp0, occerr, alpha, newerr
62  real(dp) :: cnst, ofactor, mustep, preverr
63  real(dp), allocatable :: trace(:), gbnd(:)
64  character(20) :: bml_type
65  integer :: n, m, i, iter, muadj, prev, maxiter, nsiter, nsdm
66 
67  bml_type = bml_get_type(h_bml)
68  n = bml_get_n(h_bml)
69  m = bml_get_m(h_bml)
70 
71  allocate(trace(2))
72  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, p2_bml)
73  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, d_bml)
74  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, w_bml)
75  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, aux_bml)
76  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, y_bml)
77  call bml_identity_matrix(bml_type, bml_element_real, dp, n, m, i_bml)
78  call bml_identity_matrix(bml_type, bml_element_real, dp, n, m, ai_bml)
79 
80  occerr = 10.0_dp
81  newerr = 1000_dp
82  preverr = 1000_dp
83  alpha = 1.0_dp
84  prev = 1
85  iter = 0
86  maxiter = 30
87  cnst = beta/(1.0_dp*2**(nsteps+2))
88 
89  if (scf_it .eq. 1) then
90  alpha = 4.0_dp
91  ! Normalization
92  ! P0 = 0.5*I - cnst*(H0-mu0*I)
93  call bml_copy(h_bml, p_bml)
94  call prg_normalize_implicit_fermi(p_bml, cnst, mu)
95  ! Generate good starting guess for (2*(P2-P)+1)^-1 using conjugate gradient
96  call bml_multiply_x2(p_bml, p2_bml, threshold, trace)
97  ! Y = 2*(P2-P) + II
98  call bml_copy(p2_bml, y_bml)
99  call bml_add(y_bml, p_bml, 1.0_dp, -1.0_dp, threshold)
100  call bml_scale_add_identity(y_bml, 2.0_dp, 1.0_dp, threshold)
101  call prg_conjgrad(y_bml, ai_bml, i_bml, aux_bml, d_bml, w_bml, 0.0001_dp, threshold)
102  do i = 1, nsteps
103  call bml_copy(i_bml, inv_bml(i))
104  enddo
105  end if
106 
107  do while ((occerr .gt. occerrlimit .or. muadj .eq. 1) .and. iter < maxiter)
108  iter = iter + 1
109  muadj = 0
110  write(*,*) 'mu =', mu
111  ! Normalization
112  ! P0 = 0.5*I - cnst*(H0-mu0*I)
113  call bml_copy(h_bml, p_bml)
114  call prg_normalize_implicit_fermi(p_bml, cnst, mu)
115 
116  nsdm = 0
117  do i = 1, nsteps
118  call bml_multiply_x2(p_bml, p2_bml, threshold, trace)
119  ! Y = 2*(P2-P) + I
120  call bml_copy(p2_bml, y_bml)
121  call bml_add(y_bml, p_bml, 1.0_dp, -1.0_dp, threshold)
122  call bml_scale_add_identity(y_bml, 2.0_dp, 1.0_dp, threshold)
123  ! Find inverse ai = (2*(P2-P)+I)^-1
124  !call prg_conjgrad(y_bml, Inv_bml(i), I_bml, w_bml, d_bml, aux_bml, tol, threshold)
125  !call bml_copy(Inv_bml(i),ai_bml)
126  if (iter .eq. 1) then
127  call prg_conjgrad(y_bml, inv_bml(i), i_bml, aux_bml, d_bml, w_bml, 0.001_dp, threshold)
128  endif
129  call prg_newtonschulz(y_bml, inv_bml(i), d_bml, w_bml, aux_bml, i_bml, tol, threshold, nsiter)
130  nsdm = nsdm + nsiter
131  call bml_multiply(inv_bml(i), p2_bml, p_bml, 1.0_dp, 0.0_dp, threshold)
132  !call bml_copy(ai_bml, Inv_bml(i)) ! Save inverses for use in perturbation response calculation
133  enddo
134 
135  write(*,*) 'Number of Newton-Schulz iterations:',nsdm
136  totns = totns + nsdm
137 
138  trp0 = bml_trace(p_bml)
139  trdpdmu = beta*(trp0 - bml_sum_squares(p_bml)) ! sum p(i,j)**2
140  occerr = abs(trp0 - nocc)
141  write(*,*) 'occerr =', occerr
142 
143  ! If occupation error is too large, do bisection method
144  if (occerr > 1.0_dp) then
145  if (nocc-trp0 < 0.0_dp .and. prev .eq. -1) then
146  prev = -1
147  else if (nocc-trp0 > 0.0_dp .and. prev .eq. 1) then
148  prev = 1
149  else if (nocc-trp0 > 0.0_dp .and. prev .eq. -1) then
150  prev = 1
151  alpha = alpha/2
152  else
153  prev = -1
154  alpha = alpha/2
155  endif
156  mu = mu + prev*alpha
157  muadj = 1
158 
159  ! Otherwise do Newton
160  else if (occerr .gt. occerrlimit) then
161  mustep = (nocc -trp0)/trdpdmu
162  mu = mu + mustep
163  muadj = 1
164  preverr = occerr
165  end if
166  enddo
167 
168  if (iter .ge. maxiter) then
169  write(*,*) 'Could not converge chemical potential in prg_impplicit_fermi_save_inverse'
170  end if
171  ! Adjusting the occupation sometimes causes the perturbation calculation to not converge.
172  ! For now we recompute the DM one extra time if mu was adjusted.
173  !if (muadj .eq. 1) then
174  ! Adjust occupation
175  ! call bml_copy(p_bml, d_bml)
176  ! call bml_scale_add_identity(d_bml, -1.0_dp, 1.0_dp, threshold)
177  ! call bml_multiply(p_bml, d_bml, w_bml, 1.0_dp, 0.0_dp, threshold)
178  ! ofactor = ((nocc - trP0)/trdPdmu) * beta
179  ! call bml_add(p_bml, w_bml, 1.0_dp, ofactor, threshold)
180  !end if
181  occiter = occiter + iter
182  call bml_scale(2.0_dp,p_bml)
183  deallocate(trace)
184 
185  call bml_deallocate(p2_bml)
186  call bml_deallocate(w_bml)
187  call bml_deallocate(d_bml)
188  call bml_deallocate(y_bml)
189  call bml_deallocate(aux_bml)
190  call bml_deallocate(ai_bml)
191  call bml_deallocate(i_bml)
192 
193  end subroutine prg_implicit_fermi_save_inverse
194 
209  subroutine prg_implicit_fermi(h_bml, p_bml, nsteps, k, nocc, &
210  mu, beta, method, osteps, occErrLimit, threshold, tol)
211 
212  implicit none
213 
214  type(bml_matrix_t), intent(in) :: h_bml
215  type(bml_matrix_t), intent(inout) :: p_bml
216  integer, intent(in) :: osteps, nsteps, method, k
217  real(dp), intent(in) :: nocc, threshold
218  real(dp), intent(in) :: tol
219  real(dp), intent(in) :: occerrlimit, beta
220  real(dp), intent(inout) :: mu
221 
222  type(bml_matrix_t) :: w_bml, y_bml, d_bml, p2_bml, aux1_bml, aux2_bml, i_bml, ai_bml
223  real(dp) :: trdpdmu, trp0, occerr
224  real(dp) :: cnst, ofactor
225  real(dp), allocatable :: trace(:), gbnd(:)
226  character(20) :: bml_type
227  integer :: n, m, i, iter, exp_order
228 
229  bml_type = bml_get_type(h_bml)
230  n = bml_get_n(h_bml)
231  m = bml_get_m(h_bml)
232 
233  allocate(trace(2))
234  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, p2_bml)
235  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, d_bml)
236  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, w_bml)
237  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, y_bml)
238  if (k .ge. 2) then
239  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, aux1_bml)
240  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, aux2_bml)
241  endif
242  if (method .eq. 1) then
243  call bml_identity_matrix(bml_type, bml_element_real, dp, n, m, i_bml)
244  call bml_identity_matrix(bml_type, bml_element_real, dp, n, m, ai_bml)
245  endif
246 
247  occerr = 10000.0_dp
248  iter = 0
249  exp_order = k**nsteps
250  cnst = beta/(4*exp_order)
251 
252 
253  do while ((osteps .eq. 0 .and. occerr .gt. occerrlimit) .or. &
254  (osteps .gt. 0 .and. iter .lt. osteps))
255  iter = iter + 1
256 
257  ! Normalization
258  ! P0 = 0.5*II - cnst*(H0-mu0*II)
259  call bml_copy(h_bml, p_bml)
260  call prg_normalize_implicit_fermi(p_bml, cnst, mu)
261 
262  if (method .eq. 0) then
263  write(*,*) "Doing CG"
264  do i = 1, nsteps
265 
266  if (k .eq. 2) then
267  call bml_multiply_x2(p_bml, p2_bml, threshold, trace)
268 
269  ! Y = 2*(P2-P) + II
270  call bml_copy(p2_bml, y_bml)
271  call bml_add(y_bml, p_bml, 1.0_dp, -1.0_dp, threshold)
272  call bml_scale_add_identity(y_bml, 2.0_dp, 1.0_dp, threshold)
273  else
274  call prg_setup_linsys(p_bml, y_bml, p2_bml, d_bml, w_bml, aux1_bml, aux2_bml, k, threshold)
275  end if
276  call prg_conjgrad(y_bml, p_bml, p2_bml, d_bml, aux1_bml, w_bml, tol, threshold)
277  enddo
278  else
279  write(*,*) "Doing NS"
280  do i = 1, nsteps
281 
282  if (k .eq. 2) then
283  call bml_multiply_x2(p_bml, p2_bml, threshold, trace)
284 
285  ! Y = 2*(P2-P) + II
286  call bml_copy(p2_bml, y_bml)
287  call bml_add(y_bml, p_bml, 1.0_dp, -1.0_dp, threshold)
288  call bml_scale_add_identity(y_bml, 2.0_dp, 1.0_dp, threshold)
289  else
290  call prg_setup_linsys(p_bml, y_bml, p2_bml, d_bml, w_bml, aux1_bml, aux2_bml, k, threshold)
291  end if
292  if (i .eq. 1) then
293  call prg_conjgrad(y_bml, ai_bml, i_bml, aux1_bml, d_bml, w_bml, 0.9_dp, threshold)
294  end if
295  !call prg_newtonschulz(y_bml, ai_bml, d_bml, w_bml, aux1_bml, I_bml, tol, threshold)
296  call bml_multiply(ai_bml, p2_bml, p_bml, 1.0_dp, 0.0_dp, threshold)
297  enddo
298 
299  end if
300  trdpdmu = bml_trace(p_bml)
301  trp0 = trdpdmu
302  trdpdmu = trdpdmu - bml_sum_squares(p_bml) ! sum p(i,j)**2
303  trdpdmu = beta * trdpdmu
304  occerr = abs(trp0 - nocc)
305  if (occerr .gt. occerrlimit) then
306  mu = mu + (nocc - trp0)/trdpdmu
307  end if
308  write(*,*) "mu =", mu
309  enddo
310 
311  ! Adjust occupation
312  ! X = II-P0
313  call bml_copy(p_bml, d_bml)
314  call bml_scale_add_identity(d_bml, -1.0_dp, 1.0_dp, threshold)
315 
316  call bml_multiply(p_bml, d_bml, w_bml, 1.0_dp, 0.0_dp, threshold)
317  ofactor = ((nocc - trp0)/trdpdmu) * beta
318  !call bml_add(p_bml, w_bml, 1.0_dp, ofactor, threshold)
319  !call bml_print_matrix("P adjusted occupation",p_bml,0,10,0,10)
320 
321  deallocate(trace)
322 
323  call bml_deallocate(p2_bml)
324  call bml_deallocate(w_bml)
325  call bml_deallocate(d_bml)
326  call bml_deallocate(y_bml)
327  if (k .ge. 2) then
328  call bml_deallocate(aux1_bml)
329  call bml_deallocate(aux2_bml)
330  endif
331  if (method .eq. 1) then
332  call bml_deallocate(ai_bml)
333  call bml_deallocate(i_bml)
334  endif
335 
336  end subroutine prg_implicit_fermi
337 
347  subroutine prg_implicit_fermi_zero(h_bml, p_bml, nsteps, mu, method, threshold, tol)
348 
349  implicit none
350 
351  type(bml_matrix_t), intent(in) :: h_bml
352  type(bml_matrix_t), intent(inout) :: p_bml
353  integer, intent(in) :: nsteps, method
354  real(dp), intent(in) :: mu, threshold
355  real(dp), intent(inout), optional :: tol
356 
357  type(bml_matrix_t) :: w_bml, y_bml, c_bml, d_bml, p2_bml, aux1_bml, aux2_bml, i_bml, ai_bml
358  real(dp) :: cnst
359  real(dp), allocatable :: trace(:), gbnd(:)
360  character(20) :: bml_type
361  integer :: n, m, i
362 
363  bml_type = bml_get_type(h_bml)
364  n = bml_get_n(h_bml)
365  m = bml_get_m(h_bml)
366 
367  allocate(trace(2))
368  allocate(gbnd(2))
369  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, p2_bml)
370  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, d_bml)
371  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, w_bml)
372  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, y_bml)
373  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, c_bml)
374  if (method .eq. 1) then
375  call bml_identity_matrix(bml_type, bml_element_real, dp, n, m, i_bml)
376  call bml_identity_matrix(bml_type, bml_element_real, dp, n, m, ai_bml)
377  endif
378 
379  call bml_copy(h_bml, p_bml)
380  call bml_gershgorin(p_bml, gbnd)
381  cnst = 0.5*min(1/(mu-gbnd(1)),1/(gbnd(2)-mu))
382  call prg_normalize_implicit_fermi(p_bml, cnst, mu)
383 
384  if (method .eq. 0) then
385  write(*,*) "Doing CG"
386  do i = 1, nsteps
387 
388  call bml_multiply_x2(p_bml, p2_bml, threshold, trace)
389 
390  ! Y = 2*(P2-P) + II
391  call bml_copy(p2_bml, y_bml)
392  call bml_add(y_bml, p_bml, 1.0_dp, -1.0_dp, threshold)
393  call bml_scale_add_identity(y_bml, 2.0_dp, 1.0_dp, threshold)
394  call prg_conjgrad(y_bml, p_bml, p2_bml, d_bml, w_bml, c_bml, tol, threshold)
395  enddo
396  else
397  write(*,*) "Doing NS"
398  do i = 1, nsteps
399 
400  ! Y = 2*(P2-P) + II
401  call bml_copy(p2_bml, y_bml)
402  call bml_add(y_bml, p_bml, 1.0_dp, -1.0_dp, threshold)
403  call bml_scale_add_identity(y_bml, 2.0_dp, 1.0_dp, threshold)
404  if (i .eq. 1) then
405  call prg_conjgrad(y_bml, ai_bml, i_bml, c_bml, d_bml, w_bml, 0.9_dp, threshold)
406  end if
407  !call prg_newtonschulz(y_bml, ai_bml, d_bml, w_bml, c_bml, I_bml, tol, threshold)
408  call bml_multiply(ai_bml, p2_bml, p_bml, 1.0_dp, 0.0_dp, threshold)
409  enddo
410  endif
411 
412  deallocate(gbnd)
413  deallocate(trace)
414 
415  call bml_deallocate(p2_bml)
416  call bml_deallocate(w_bml)
417  call bml_deallocate(d_bml)
418  call bml_deallocate(y_bml)
419  call bml_deallocate(c_bml)
420  if (method .eq. 1) then
421  call bml_deallocate(ai_bml)
422  call bml_deallocate(i_bml)
423  endif
424 
425  end subroutine prg_implicit_fermi_zero
426 
427 
439  subroutine prg_implicit_fermi_first_order_response(H0_bml, H1_bml, P0_bml, P1_bml, &
440  Inv_bml, nsteps, mu0, beta, nocc, threshold)
441 
442  implicit none
443 
444  type(bml_matrix_t), intent(in) :: h0_bml, h1_bml, inv_bml(nsteps)
445  type(bml_matrix_t), intent(inout) :: p0_bml,p1_bml
446  real(dp), intent(in) :: mu0, threshold
447  real(dp) :: mu1
448  real(dp), intent(in) :: beta, nocc
449  integer, intent(in) :: nsteps
450  type(bml_matrix_t) :: b0_bml, b_bml, c_bml, c0_bml
451  character(20) :: bml_type
452  real(dp) :: p1_trace, dpdmu_trace, p1b_trace, mu1b, cnst
453  integer :: n, m, i, j, k
454 
455  bml_type = bml_get_type(h0_bml)
456  n = bml_get_n(h0_bml)
457  m = bml_get_m(h0_bml)
458 
459  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, b0_bml)
460  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, b_bml)
461  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, c_bml)
462  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, c0_bml)
463 
464  cnst = beta/(2**(2+nsteps))
465 
466  ! P0 = 0.5*II - cnst*(H0-mu0*II)
467  call bml_copy(h0_bml, p0_bml)
468  call prg_normalize_implicit_fermi(p0_bml, cnst, mu0)
469 
470  ! P1 = - cnst*H1
471  call bml_copy(h1_bml, p1_bml)
472  call bml_scale(-1.0_dp*cnst, p1_bml)
473  do i = 1, nsteps
474 
475  ! Calculate coefficient matrices
476  ! C0 = P0^2
477  call bml_multiply(p0_bml, p0_bml, c0_bml, 1.0_dp, 0.0_dp, threshold)
478  ! C = P0*P1+P1*P0, B = 2(P1 - C)
479  call bml_multiply(p0_bml, p1_bml, c_bml, 1.0_dp, 0.0_dp, threshold)
480  call bml_multiply(p1_bml, p0_bml, c_bml, 1.0_dp, 1.0_dp, threshold)
481  call bml_copy(p1_bml, b_bml)
482  call bml_add(b_bml, c_bml, 2.0_dp, -2.0_dp, threshold)
483  ! Get next P0
484  call bml_multiply(inv_bml(i), c0_bml, p0_bml, 1.0_dp, 0.0_dp, threshold)
485  ! Get next P1
486  ! C = P0*P1+P1*P0 + 2(P1 -P0*P1-P1*P0)*P0(i+1)
487  call bml_multiply(b_bml, p0_bml, c_bml, 1.0_dp, 1.0_dp, threshold)
488  call bml_multiply(inv_bml(i), c_bml, p1_bml, 1.0_dp, 0.0_dp, threshold)
489  enddo
490 
491  ! do i = 1, nsteps-1
492  ! D = A^-1*P0
493  ! call bml_multiply(Inv_bml(i), B0_bml, C0_bml, 1.0_dp, 0.0_dp, threshold)
494  ! call bml_multiply(C0_bml, B0_bml, B_bml, 1.0_dp, 0.0_dp, threshold)
495  ! B0 = A^-1*P0^2
496  ! call bml_copy(B_bml,B0_bml)
497  ! B = I + D -P0*D
498  ! call bml_add(B_bml, C0_bml, -1.0_dp, 1.0_dp, threshold)
499  ! call bml_scale_add_identity(B_bml, 1.0_dp, 1.0_dp, threshold)
500  ! P1 = 2D*P1(I+D-P0*D)
501  ! call bml_multiply(C0_bml, P1_bml, C_bml, 1.0_dp, 0.0_dp, threshold)
502  ! call bml_multiply(C_bml, B_bml, P1_bml, 2.0_dp, 0.0_dp, threshold)
503  ! enddo
504  ! call bml_multiply(B0_bml, P1_bml, C_bml, 2.0_dp, 0.0_dp, threshold)
505  ! call bml_copy(P1_bml, B_bml)
506  ! call bml_add(B_bml, C_bml, 2.0_dp, -2.0_dp, threshold)
507  ! Get next P1
508  ! call bml_multiply(B_bml, P0_bml, C_bml, 1.0_dp, 1.0_dp, threshold)
509  ! call bml_multiply(Inv_bml(i), C_bml, P1_bml, 1.0_dp, 0.0_dp, threshold)
510 
511 
512  ! dPdmu = beta*P0(I-P0)
513  call bml_copy(p0_bml, b_bml)
514  call bml_scale_add_identity(b_bml, -1.0_dp, 1.0_dp, threshold)
515  call bml_multiply(p0_bml, b_bml, c_bml, beta, 0.0_dp, threshold)
516  dpdmu_trace = bml_trace(c_bml)
517  p1_trace = bml_trace(p1_bml)
518  mu1 = - p1_trace/dpdmu_trace
519  if (abs(dpdmu_trace) > 1e-8) then
520  call bml_add(p1_bml,c_bml,1.0_dp,mu1,threshold)
521  endif
522 
523  call bml_deallocate(b_bml)
524  call bml_deallocate(b0_bml)
525  call bml_deallocate(c_bml)
526  call bml_deallocate(c0_bml)
527 
528 
529  end subroutine prg_implicit_fermi_first_order_response
530 
531 
547  subroutine prg_implicit_fermi_response(H0_bml, H1_bml, H2_bml, H3_bml, P0_bml, P1_bml, P2_bml, P3_bml, &
548  nsteps, mu0, mu, beta, nocc, occ_tol, lin_tol, order, threshold)
549 
550  implicit none
551 
552  type(bml_matrix_t), intent(in) :: h0_bml, h1_bml, h2_bml, h3_bml
553  type(bml_matrix_t), intent(inout) :: p0_bml, p1_bml, p2_bml, p3_bml
554  real(dp), intent(inout) :: mu0
555  real(dp), allocatable, intent(inout) :: mu(:)
556  real(dp), intent(in) :: beta, occ_tol, lin_tol, nocc
557  integer, intent(in) :: nsteps
558  type(bml_matrix_t) :: i_bml, tmp1_bml, tmp2_bml, tmp3_bml, c0_bml, t_bml, ti_bml
559  type(bml_matrix_t), allocatable :: b_bml(:), p_bml(:), c_bml(:), h_bml(:)
560  real(dp), allocatable :: p_trace(:), trace(:)
561  character(20) :: bml_type
562  real(dp) :: occ_err, p0_trace, pmu_trace, cnst, threshold, tol, lambda, h
563  integer :: n, m, order, i, j, k
564 
565  k = 0
566  occ_err = 10000.0
567  allocate(p_trace(order))
568  allocate(b_bml(order))
569  allocate(c_bml(order))
570  allocate(p_bml(order))
571  allocate(h_bml(order))
572 
573  bml_type = bml_get_type(h0_bml)
574  n = bml_get_n(h0_bml)
575  m = bml_get_m(h0_bml)
576 
577  do i = 1, order
578  mu(i) = 0.0_dp
579  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, b_bml(i))
580  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, c_bml(i))
581  end do
582 
583  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, tmp1_bml)
584  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, tmp3_bml)
585  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, tmp2_bml)
586  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, c0_bml)
587  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, t_bml)
588  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, ti_bml)
589  call bml_identity_matrix(bml_type, bml_element_real, dp, n, m, i_bml)
590 
591  h_bml(1) = h1_bml
592  p_bml(1) = p1_bml
593  if (order .gt. 1) then
594  h_bml(2) = h2_bml
595  p_bml(2) = p2_bml
596  end if
597  if (order .gt. 2) then
598  h_bml(3) = h3_bml
599  p_bml(3) = p3_bml
600  end if
601 
602  cnst = beta/(2**(2+nsteps))
603 
604  do while (occ_err .gt. occ_tol)
605  k = k + 1
606  ! P0 = 0.5*II - cnst*(H0-mu0*II)
607  call bml_copy(h0_bml, p0_bml)
608  call prg_normalize_implicit_fermi(p0_bml, cnst, mu0)
609 
610  ! P(j) = - cnst*(H(j)-mu(j)*II)
611  do j = 1, order
612  call bml_copy(h_bml(j), p_bml(j))
613  call prg_normalize_implicit_fermi(p_bml(j), cnst, mu(j))
614  call bml_scale_add_identity(p_bml(j), 1.0_dp, -0.5_dp, threshold)
615  end do
616 
617  do i = 1, nsteps
618 
619  ! Calculate coefficient matrices
620 
621  ! C0 = P0^2
622  call bml_multiply(p0_bml, p0_bml, c0_bml, 1.0_dp, 0.0_dp, threshold)
623  ! C1 = P0*P1+P1*P0, B1 = 2(P1 - C1)
624  call bml_multiply(p0_bml, p_bml(1), c_bml(1), 1.0_dp, 0.0_dp, threshold)
625  call bml_multiply(p_bml(1), p0_bml, c_bml(1), 1.0_dp, 1.0_dp, threshold)
626  call bml_copy(p_bml(1), b_bml(1))
627  call bml_add(b_bml(1), c_bml(1), 2.0_dp, -2.0_dp, threshold)
628  if (order > 1) then
629  ! C2 = P1^2 + P0*P2 + P2*P0, B2 = 2(P2 - C2)
630  call bml_multiply(p_bml(1), p_bml(1), c_bml(2), 1.0_dp, 0.0_dp, threshold)
631  call bml_multiply(p0_bml, p_bml(2), c_bml(2), 1.0_dp, 1.0_dp, threshold)
632  call bml_multiply(p_bml(2), p0_bml, c_bml(2), 1.0_dp, 1.0_dp, threshold)
633  call bml_copy(p_bml(2), b_bml(2))
634  call bml_add(b_bml(2), c_bml(2), 2.0_dp, -2.0_dp, threshold)
635  end if
636  if (order > 2) then
637  ! C3 = P1*P2 + P2+P1 + P0*P3 + P3*P0, B3 = 2(P3 - C3)
638  call bml_multiply(p_bml(1), p_bml(2), c_bml(3), 1.0_dp, 0.0_dp, threshold)
639  call bml_multiply(p_bml(2), p_bml(1), c_bml(3), 1.0_dp, 1.0_dp, threshold)
640  call bml_multiply(p0_bml, p_bml(3), c_bml(3), 1.0_dp, 1.0_dp, threshold)
641  call bml_multiply(p_bml(3), p0_bml, c_bml(3), 1.0_dp, 1.0_dp, threshold)
642  call bml_copy(p_bml(3), b_bml(3))
643  call bml_add(b_bml(3), c_bml(3), 2.0_dp, -2.0_dp, threshold)
644  endif
645  ! T = 2P0^2 - 2P0 + I
646  call bml_copy(c0_bml, t_bml)
647  call bml_add(t_bml, p0_bml, 1.0_dp, -1.0_dp, threshold)
648  call bml_scale_add_identity(t_bml, 2.0_dp, 1.0_dp, threshold)
649  ! Find T-inverse
650  if (i .eq. 1) then
651  call prg_conjgrad(t_bml, ti_bml, i_bml, tmp1_bml, tmp2_bml, tmp3_bml,0.01_dp, threshold)
652  call bml_identity_matrix(bml_type, bml_element_real, dp, n, m, i_bml)
653  end if
654  !call prg_newtonschulz(T_bml, Ti_bml, tmp1_bml, tmp2_bml, tmp3_bml, I_bml, lin_tol, threshold)
655  ! Get next P0
656  call bml_multiply(ti_bml, c0_bml, p0_bml, 1.0_dp, 0.0_dp, threshold)
657  ! Get next P1
658  call bml_multiply(b_bml(1), p0_bml, c_bml(1), 1.0_dp, 1.0_dp, threshold)
659  call bml_multiply(ti_bml, c_bml(1), p_bml(1), 1.0_dp, 0.0_dp, threshold)
660  if (order > 1) then
661  ! Get next P2
662  call bml_multiply(b_bml(2), p0_bml, c_bml(2), 1.0_dp, 1.0_dp, threshold)
663  call bml_multiply(b_bml(1), p_bml(1), c_bml(2), 1.0_dp, 1.0_dp, threshold)
664  call bml_multiply(ti_bml, c_bml(2), p_bml(2), 1.0_dp, 0.0_dp, threshold)
665  end if
666  if (order > 2) then
667  ! Get next P3
668  call bml_multiply(b_bml(3), p0_bml, c_bml(3), 1.0_dp, 1.0_dp, threshold)
669  call bml_multiply(b_bml(2), p_bml(1), c_bml(3), 1.0_dp, 1.0_dp, threshold)
670  call bml_multiply(b_bml(1), p_bml(2), c_bml(3), 1.0_dp, 1.0_dp, threshold)
671  call bml_multiply(ti_bml, c_bml(3), p_bml(3), 1.0_dp, 0.0_dp, threshold)
672  endif
673  enddo
674 
675  ! Pmu = beta*P0(I-P0)
676  call bml_copy(p0_bml, tmp1_bml)
677  call bml_scale_add_identity(tmp1_bml, -1.0_dp, 1.0_dp, threshold)
678  call bml_multiply(p0_bml, tmp1_bml, tmp2_bml, beta, 0.0_dp, threshold)
679 
680  pmu_trace = bml_trace(tmp2_bml)
681  p0_trace = bml_trace(p0_bml)
682  occ_err = abs(p0_trace-nocc)
683  mu0 = mu0 + (nocc - p0_trace)/pmu_trace
684  do i = 1, order
685  p_trace(i) = bml_trace(p_bml(i))
686  mu(i) = mu(i) - p_trace(i)/pmu_trace
687  occ_err = occ_err + abs(p_trace(i))
688  enddo
689 
690  write(*,*) "occ_err =", occ_err
691  if (k .gt. 50) then
692  write(*,*) "Chemical potential is not converging"
693  exit
694  endif
695 
696  enddo
697 
698  do i = 1, order
699  call bml_deallocate(b_bml(i))
700  call bml_deallocate(c_bml(i))
701  end do
702 
703  call bml_deallocate(i_bml)
704  call bml_deallocate(tmp1_bml)
705  call bml_deallocate(tmp2_bml)
706  call bml_deallocate(t_bml)
707  call bml_deallocate(ti_bml)
708  deallocate(p_trace)
709  deallocate(b_bml)
710  deallocate(c_bml)
711 
712  end subroutine prg_implicit_fermi_response
713 
715  !using finite differences.
716  !! \param H0_bml Input Hamiltonian matrix.
717  !! \param H_list Input List of one to third order Hamiltonian perturbations
718  !! \param mu0 Shifted chemical potential.
719  !! \param mu List of first to third order perturbations in the chemical
720  !! potential.
721  !! \param beta Input inverse temperature.
722  !! \param order Calculate response up to this order.
723  !! \param lambda Perturbation parameter
724  !! \param h Finite difference step size
725  !! \param threshold Threshold for matrix algebra.
726  subroutine prg_finite_diff(H0_bml, H_list, mu0, mu_list, beta, order, lambda, h, threshold)
727 
728  implicit none
729 
730  type(bml_matrix_t), intent(in) :: h0_bml
731  real(dp), intent(in) :: mu0
732  type(bml_matrix_t), allocatable, intent(in) :: h_list(:)
733  real(dp), allocatable, intent(in) :: mu_list(:)
734  real(dp), intent(in) :: lambda, beta, threshold, h
735  integer, intent(in) :: order
736  character(20) :: bml_type
737  real(dp) :: mu_1minus, mu_1plus, mu_2minus, mu_2plus, mu_central
738  real(dp) :: lambda_f, lambda_b, lambda_2f, lambda_2b
739  integer :: n, m, i
740  type(bml_matrix_t) :: d0_bml, d1minus_bml, d1plus_bml, d2plus_bml, d2minus_bml, tmp1_bml
741 
742  bml_type = bml_get_type(h0_bml)
743  n = bml_get_n(h0_bml)
744  m = bml_get_m(h0_bml)
745 
746  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, tmp1_bml)
747  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, d0_bml)
748  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, d1minus_bml)
749  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, d1plus_bml)
750  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, d2plus_bml)
751  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, d2minus_bml)
752 
753  lambda_f = lambda+h
754  lambda_b = lambda-h
755  lambda_2f = lambda+2*h
756  lambda_2b = lambda-2*h
757 
758  ! Calculate first order difference
759  ! P1 = (F(H0 + lambda_f*H1, mu0 + lambda_f*mu(1)) - F(H0 + lambda*H1, mu0 +
760  ! lambda*H1))/h
761  call bml_copy(h0_bml, tmp1_bml)
762  mu_central = mu0
763  do i = 1, order
764  call bml_add(tmp1_bml, h_list(i), 1.0_dp, lambda**i, threshold)
765  mu_central = mu_central + lambda**i*mu_list(i)
766  end do
767  call prg_get_density_matrix(tmp1_bml, d0_bml, beta, mu_central, threshold)
768 
769  call bml_copy(h0_bml, tmp1_bml)
770  mu_1plus = mu0
771  do i = 1, order
772  call bml_add(tmp1_bml, h_list(i), 1.0_dp, lambda_f**i, threshold)
773  mu_1plus = mu_1plus + lambda_f**i*mu_list(i)
774  end do
775  call prg_get_density_matrix(tmp1_bml, d1plus_bml, beta, mu_1plus, threshold)
776  call bml_copy(d1plus_bml, tmp1_bml)
777  call bml_add(tmp1_bml, d0_bml, 1.0_dp/h, -1.0_dp/h, threshold)
778 
779  call bml_scale(1000.0_dp, tmp1_bml)
780  call bml_print_matrix("Finite diff - Order 1 * 1000", tmp1_bml, 0,10,0,10)
781 
782  if (order .gt. 1) then
783  ! Calculate second order difference
784  call bml_copy(h0_bml, tmp1_bml)
785  mu_1minus = mu0
786  do i = 1, order
787  call bml_add(tmp1_bml, h_list(i), 1.0_dp, lambda_b**i, threshold)
788  mu_1minus = mu_1minus + lambda_b**i*mu_list(i)
789  end do
790  call prg_get_density_matrix(tmp1_bml, d1minus_bml, beta, mu_1minus, threshold)
791  call bml_copy(d0_bml, tmp1_bml)
792  call bml_add(tmp1_bml, d1minus_bml, -1.0_dp/(h*h), 0.5_dp/(h*h), threshold)
793  call bml_add(tmp1_bml, d1plus_bml, 1.0_dp, 0.5_dp/(h*h))
794 
795  call bml_scale(1000.0_dp, tmp1_bml)
796  call bml_print_matrix("Finite diff - Order 2 * 1000", tmp1_bml, 0,10,0,10)
797  end if
798 
799  if (order .gt. 2) then
800  ! Calculate third order difference
801  call bml_copy(h0_bml, tmp1_bml)
802  mu_2minus = mu0
803  do i = 1, order
804  call bml_add(tmp1_bml, h_list(i), 1.0_dp, lambda_2b**i, threshold)
805  mu_2minus = mu_2minus + lambda_2b**i*mu_list(i)
806  end do
807  call prg_get_density_matrix(tmp1_bml, d2minus_bml, beta, mu_2minus, threshold)
808  call bml_copy(h0_bml, tmp1_bml)
809  mu_2plus = mu0
810  do i = 1, order
811  call bml_add(tmp1_bml, h_list(i), 1.0_dp, lambda_2f**i, threshold)
812  mu_2plus = mu_2plus + lambda_2f**i*mu_list(i)
813  end do
814  call prg_get_density_matrix(tmp1_bml, d2plus_bml, beta, mu_2plus, threshold)
815  call bml_copy(d2plus_bml, tmp1_bml)
816  call bml_add(tmp1_bml, d2minus_bml, 1.0_dp/(12.0_dp*h**3), -1.0_dp/(12.0_dp*h**3), threshold)
817  call bml_add(tmp1_bml, d1plus_bml, 1.0_dp, -1.0/(6.0_dp*h**3), threshold)
818  call bml_add(tmp1_bml, d1minus_bml, 1.0_dp, 1.0/(6.0_dp*h**3), threshold)
819 
820  call bml_scale(1000.0_dp, tmp1_bml)
821  call bml_print_matrix("Finite diff - Order 3 * 1000", tmp1_bml, 0,10,0,10)
822  end if
823 
824  call bml_deallocate(tmp1_bml)
825  call bml_deallocate(d0_bml)
826  call bml_deallocate(d1minus_bml)
827  call bml_deallocate(d1plus_bml)
828  call bml_deallocate(d2plus_bml)
829  call bml_deallocate(d2minus_bml)
830 
831  end subroutine prg_finite_diff
832 
843  subroutine prg_setup_linsys(p_bml, A_bml, b_bml, p2_bml, y_bml, aux_bml, &
844  aux1_bml, k, threshold)
845 
846  implicit none
847 
848  type(bml_matrix_t), intent(inout) :: A_bml, b_bml, p2_bml, y_bml, aux_bml, aux1_bml
849  type(bml_matrix_t), intent(in) :: p_bml
850  real(dp), intent(in) :: threshold
851  integer, intent(in) :: k
852  character(20) :: bml_type
853  integer :: M, N, i
854 
855  if (k .eq. 2) then
856  call bml_multiply(p_bml, p_bml, b_bml, 1.0_dp, 0.0_dp, threshold)
857  call bml_copy(b_bml, a_bml)
858  call bml_add(a_bml, p_bml, 2.0_dp, -2.0_dp, threshold)
859  call bml_scale_add_identity(a_bml, 1.0_dp, 1.0_dp, threshold)
860 
861  else
862  call bml_multiply(p_bml, p_bml, p2_bml, 1.0_dp, 0.0_dp, threshold)
863  call bml_copy(p2_bml, y_bml)
864  call bml_add(y_bml, p_bml, 1.0_dp, -2.0_dp, threshold)
865  call bml_scale_add_identity(y_bml, 1.0_dp, 1.0_dp, threshold)
866 
867  call bml_copy(p2_bml, b_bml)
868  call bml_copy(y_bml, a_bml)
869  do i = 1,(k-2)/2
870  call bml_multiply(b_bml, p2_bml, aux_bml, 1.0_dp, 0.0_dp, threshold)
871  call bml_multiply(a_bml, y_bml, aux1_bml, 1.0_dp, 0.0_dp, threshold)
872  call bml_copy(aux_bml, b_bml)
873  call bml_copy(aux1_bml, a_bml)
874  enddo
875  call bml_add(a_bml, b_bml, 1.0_dp, 1.0_dp, threshold)
876  end if
877 
878  end subroutine prg_setup_linsys
879 
887 
888  subroutine prg_newtonschulz(a_bml, ai_bml, r_bml, tmp_bml, d_bml, I_bml, tol, threshold, num_iter)
889 
890  implicit none
891 
892  type(bml_matrix_t), intent(inout) :: ai_bml, r_bml, tmp_bml, d_bml
893  type(bml_matrix_t), intent(in) :: a_bml, I_bml
894  real(dp), intent(in) :: threshold, tol
895  integer, intent(out) :: num_iter
896  real(dp) :: err,prev_err,scaled_tol
897  integer :: i,N,N2
898 
899  n = bml_get_n(a_bml)
900  err = 100000.0
901  i = 0
902  do while(err > tol)
903  !write(*,*) 'iter = ', i
904  !write(*,*) 'ns error =', err
905  call bml_copy(ai_bml, tmp_bml)
906  call bml_multiply(a_bml, ai_bml, r_bml, 1.0_dp, 0.0_dp, threshold)
907  call bml_scale_add_identity(r_bml, -1.0_dp, 1.0_dp, threshold)
908  prev_err = err
909  err = bml_fnorm(r_bml)
910  !write(*,*) "err = ", err
911  !write(*,*) "prev_err = ", prev_err
912  if (prev_err+0.01 < err) then
913  write(*,*) 'NS did not converge, calling conjugate gradient'
914  call prg_conjgrad(a_bml, ai_bml, i_bml, r_bml, tmp_bml, d_bml, 0.00001_dp, threshold)
915  else
916  call bml_multiply(tmp_bml, r_bml, ai_bml, 1.0_dp, 1.0_dp, threshold)
917  endif
918  i = i + 1
919  enddo
920  num_iter = i
921  !write(*,*) "Number of NS iterations:", i
922  end subroutine prg_newtonschulz
923 
924  ! Preconditioned CG, preconditioner inverse diagonal of A
934  subroutine prg_pcg(A_bml, p_bml, p2_bml, d_bml, wtmp_bml, cg_tol, threshold)
935 
936  implicit none
937 
938  type(bml_matrix_t), intent(in) :: A_bml
939  type(bml_matrix_t), intent(inout) :: p_bml, p2_bml, d_bml, wtmp_bml
940  real(dp), intent(in) :: cg_tol, threshold
941 
942  type(bml_matrix_t) :: M_bml, z_bml
943  real(dp), allocatable :: diagonal(:)
944  real(dp) :: alpha, beta
945  character(20) :: bml_type
946  integer :: k,N,M
947  real(dp) :: r_norm_old, r_norm_new
948 
949  bml_type = bml_get_type(p_bml)
950  n = bml_get_n(p_bml)
951  m = bml_get_m(p_bml)
952 
953  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, z_bml)
954  call bml_zero_matrix(bml_type, bml_element_real, dp, n, m, m_bml)
955 
956  allocate(diagonal(n))
957  call bml_get_diagonal(a_bml, diagonal)
958  do k = 1,n
959  diagonal(k) = 1.0_dp/diagonal(k)
960  enddo
961  call bml_set_diagonal(m_bml, diagonal)
962 
963  call bml_multiply(a_bml, p_bml, p2_bml, -1.0_dp, 1.0_dp, threshold)
964  call bml_multiply(m_bml, p2_bml, z_bml, 1.0_dp, 0.0_dp, threshold)
965  r_norm_new = bml_trace_mult(z_bml, p2_bml)
966  call bml_copy(z_bml, d_bml)
967  k = 0
968 
969  do while (bml_sum_squares(p2_bml) .gt. cg_tol)
970 
971  write(*,*) "r_norm", bml_sum_squares(p2_bml)
972  k = k + 1
973  if (k .ne. 1) then
974  beta = r_norm_new/r_norm_old
975  call bml_add(d_bml, z_bml, beta, 1.0_dp, threshold)
976  endif
977 
978  call bml_multiply(a_bml, d_bml, wtmp_bml, 1.0_dp, 0.0_dp, threshold)
979  alpha = bml_trace_mult(p2_bml,z_bml)/bml_trace_mult(d_bml, wtmp_bml)
980  call bml_add(p_bml, d_bml, 1.0_dp, alpha, threshold)
981  call bml_add(p2_bml, wtmp_bml, 1.0_dp, -alpha, threshold)
982  call bml_multiply(m_bml, p2_bml, z_bml, 1.0_dp, 0.0_dp, threshold)
983  r_norm_old = r_norm_new
984  r_norm_new = bml_trace_mult(p2_bml,z_bml)
985  if (k .gt. 100) then
986  write(*,*) "PCG is not converging"
987  stop
988  endif
989  enddo
990  write(*,*) "Number of iterations:", k
991 
992  call bml_deallocate(z_bml)
993  call bml_deallocate(m_bml)
994  deallocate(diagonal)
995 
996  end subroutine prg_pcg
997 
1006  subroutine prg_conjgrad(A_bml, p_bml, p2_bml, tmp_bml, d_bml, w_bml, cg_tol, threshold)
1007 
1008  implicit none
1009 
1010  type(bml_matrix_t), intent(in) :: A_bml, p2_bml
1011  type(bml_matrix_t), intent(inout) :: p_bml, tmp_bml, d_bml, w_bml
1012  real(dp), intent(in) :: cg_tol, threshold
1013 
1014  real(dp) :: alpha, beta
1015  integer :: k
1016  real(dp) :: r_norm_old, r_norm_new
1017 
1018  call bml_copy(p2_bml,tmp_bml)
1019  call bml_multiply(a_bml, p_bml, tmp_bml, -1.0_dp, 1.0_dp, threshold)
1020  r_norm_new = bml_sum_squares(tmp_bml)
1021  k = 0
1022 
1023  do while (r_norm_new .gt. cg_tol)
1024 
1025  ! write(*,*) r_norm_new
1026  k = k + 1
1027  if (k .eq. 1) then
1028  call bml_copy(tmp_bml, d_bml)
1029  else
1030  beta = r_norm_new/r_norm_old
1031  call bml_add(d_bml, tmp_bml, beta, 1.0_dp, threshold)
1032  endif
1033 
1034  call bml_multiply(a_bml, d_bml, w_bml, 1.0_dp, 0.0_dp, threshold)
1035  alpha = r_norm_new/bml_trace_mult(d_bml, w_bml)
1036 
1037  call bml_add(p_bml, d_bml, 1.0_dp, alpha, threshold)
1038  call bml_add(tmp_bml, w_bml, 1.0_dp, -alpha, threshold)
1039  r_norm_old = r_norm_new
1040  r_norm_new = bml_sum_squares(tmp_bml)
1041  if (k .gt. 50) then
1042  write(*,*) "Conjugate gradient is not converging"
1043  stop
1044  endif
1045  enddo
1046  !write(*,*) "Number of CG-iterations:", k
1047 
1048  end subroutine prg_conjgrad
1049 
1056  subroutine prg_get_density_matrix(ham_bml, p_bml, beta, mu, threshold)
1057 
1058  implicit none
1059 
1060  type(bml_matrix_t), intent(in) :: ham_bml
1061  type(bml_matrix_t), intent(inout) :: p_bml
1062  real(dp), intent(in) :: beta, threshold
1063  real(dp), intent(in) :: mu
1064  character(20) :: bml_type
1065  integer :: N, M, i
1066  real(dp), allocatable :: eigenvalues(:)
1067  type(bml_matrix_t) :: eigenvectors_bml,occupation_bml,aux_bml,aux1_bml,i_bml
1068 
1069  bml_type = bml_get_type(p_bml)
1070  n = bml_get_n(p_bml)
1071  m = bml_get_m(p_bml)
1072 
1073  allocate(eigenvalues(n))
1074 
1075  call bml_zero_matrix(bml_type,bml_element_real,dp,n,m,eigenvectors_bml)
1076  call bml_zero_matrix(bml_type,bml_element_real,dp,n,m,occupation_bml)
1077  call bml_zero_matrix(bml_type,bml_element_real,dp,n,m,aux_bml)
1078  call bml_zero_matrix(bml_type,bml_element_real,dp,n,m,aux1_bml)
1079  call bml_identity_matrix(bml_type,bml_element_real,dp,n,m,i_bml)
1080 
1081  call bml_diagonalize(ham_bml,eigenvalues,eigenvectors_bml)
1082 
1083  do i=1,n
1084  eigenvalues(i) = fermi(eigenvalues(i),mu,beta)
1085  enddo
1086 
1087  call bml_set_diagonal(occupation_bml, eigenvalues)
1088  call bml_multiply(eigenvectors_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp,threshold)
1089  call bml_transpose_new(eigenvectors_bml, aux1_bml)
1090  call bml_multiply(aux_bml, aux1_bml, p_bml, 1.0_dp, 0.0_dp, threshold)
1091 
1092  call bml_deallocate(eigenvectors_bml)
1093  call bml_deallocate(occupation_bml)
1094  call bml_deallocate(aux_bml)
1095  call bml_deallocate(aux1_bml)
1096  call bml_deallocate(i_bml)
1097 
1098  deallocate(eigenvalues)
1099 
1100  end subroutine prg_get_density_matrix
1101 
1103  !potential
1104  !! \param ham_bml Input hamiltonian
1105  !! \param p_bml Output density matrix
1106  !! \param beta Inverse temperature
1107  !! \param mu Chemical potential
1108  !! \param nocc Number of occupied states
1109  !! \param osteps Outer loop steps to converge chemical potential
1110  !! \param occErrLimit Occupation error limit.
1111  !! \param threshold Threshold for matrix algebra
1112  subroutine prg_test_density_matrix(ham_bml, p_bml, beta, mu, nocc, osteps, occErrLimit, threshold)
1113 
1114  implicit none
1115 
1116  type(bml_matrix_t), intent(in) :: ham_bml
1117  type(bml_matrix_t), intent(inout) :: p_bml
1118  real(dp), intent(in) :: beta, nocc, occerrlimit, threshold
1119  real(dp), intent(inout) :: mu
1120  integer, intent(in) :: osteps
1121  character(20) :: bml_type
1122  integer :: n, m, i, iter
1123  real(dp) :: trdpdmu, trp0, ofactor, occerr
1124  real(dp), allocatable :: eigenvalues(:)
1125  type(bml_matrix_t) :: eigenvectors_bml,occupation_bml,aux_bml,aux1_bml,i_bml
1126 
1127  bml_type = bml_get_type(p_bml)
1128  n = bml_get_n(p_bml)
1129  m = bml_get_m(p_bml)
1130 
1131  allocate(eigenvalues(n))
1132 
1133  call bml_zero_matrix(bml_type,bml_element_real,dp,n,m,eigenvectors_bml)
1134  call bml_zero_matrix(bml_type,bml_element_real,dp,n,m,occupation_bml)
1135  call bml_zero_matrix(bml_type,bml_element_real,dp,n,m,aux_bml)
1136  call bml_zero_matrix(bml_type,bml_element_real,dp,n,m,aux1_bml)
1137  call bml_identity_matrix(bml_type,bml_element_real,dp,n,m,i_bml)
1138 
1139  occerr = 1000.0_dp
1140  iter = 0
1141 
1142  do while ((osteps .eq. 0 .and. occerr .gt. occerrlimit) .or. &
1143  (osteps .gt. 0 .and. iter .lt. osteps))
1144  iter = iter + 1
1145 
1146  call bml_diagonalize(ham_bml,eigenvalues,eigenvectors_bml)
1147 
1148  do i=1,n
1149  eigenvalues(i) = fermi(eigenvalues(i),mu,beta)
1150  enddo
1151 
1152  call bml_set_diagonal(occupation_bml, eigenvalues)
1153  call bml_multiply(eigenvectors_bml, occupation_bml, aux_bml, 1.0_dp, 0.0_dp,threshold)
1154  call bml_transpose_new(eigenvectors_bml, aux1_bml)
1155  call bml_multiply(aux_bml, aux1_bml, p_bml, 1.0_dp, 0.0_dp, threshold)
1156  call bml_print_matrix("test density",p_bml,0,10,0,10)
1157  trdpdmu = bml_trace(p_bml)
1158  trp0 = trdpdmu
1159  trdpdmu = trdpdmu - bml_sum_squares(p_bml) ! sum p(i,j)**2
1160  trdpdmu = beta * trdpdmu
1161  occerr = abs(trp0 - nocc)
1162  if (occerr .gt. occerrlimit) then
1163  mu = mu + (nocc - trp0)/trdpdmu
1164  end if
1165  !write(*,*) "mu = ", mu
1166  enddo
1167 
1168  ! Adjust occupation
1169  ! X = II-P0
1170  call bml_copy(p_bml, aux_bml)
1171  call bml_scale_add_identity(aux_bml, -1.0_dp, 1.0_dp, threshold)
1172 
1173  call bml_multiply(p_bml, aux_bml, aux1_bml, 1.0_dp, 0.0_dp, threshold)
1174  ofactor = ((nocc - trp0)/trdpdmu) * beta
1175  !call bml_add(p_bml, aux1_bml, 1.0_dp, ofactor, threshold)
1176  !call bml_print_matrix("Diagonalization - Adjusted occupation",p_bml,0,10,0,10)
1177 
1178  call bml_deallocate(eigenvectors_bml)
1179  call bml_deallocate(occupation_bml)
1180  call bml_deallocate(aux_bml)
1181  call bml_deallocate(aux1_bml)
1182  call bml_deallocate(i_bml)
1183 
1184  deallocate(eigenvalues)
1185 
1186  end subroutine prg_test_density_matrix
1187 
1192  real(dp) function fermi(e,mu,beta)
1193 
1194  real(dp), intent(in) :: e, mu, beta
1195 
1196  fermi = 1.0_dp/(1.0_dp+exp(beta*(e-mu)))
1197 
1198  end function fermi
1199 
1200 end module prg_implicit_fermi_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_ewald_mod
Definition: prg_ewald_mod.F90:2
prg_implicit_fermi_mod
Definition: prg_implicit_fermi_mod.F90:6
prg_implicit_fermi_mod::fermi
real(dp) function fermi(e, mu, beta)
Gives the Fermi distribution value for energy e.
Definition: prg_implicit_fermi_mod.F90:1193
prg_implicit_fermi_mod::prg_implicit_fermi
subroutine, public prg_implicit_fermi(h_bml, p_bml, nsteps, k, nocc, mu, beta, method, osteps, occErrLimit, threshold, tol)
Recursive Implicit Fermi Dirac for finite temperature.
Definition: prg_implicit_fermi_mod.F90:211
prg_implicit_fermi_mod::prg_implicit_fermi_save_inverse
subroutine, public prg_implicit_fermi_save_inverse(Inv_bml, h_bml, p_bml, nsteps, nocc, mu, beta, occErrLimit, threshold, tol, SCF_IT, occiter, totns)
Recursive Implicit Fermi Dirac for finite temperature.
Definition: prg_implicit_fermi_mod.F90:48
prg_implicit_fermi_mod::prg_pcg
subroutine prg_pcg(A_bml, p_bml, p2_bml, d_bml, wtmp_bml, cg_tol, threshold)
Solve the system AX = B with conjugate gradient.
Definition: prg_implicit_fermi_mod.F90:935
prg_implicit_fermi_mod::prg_implicit_fermi_first_order_response
subroutine, public prg_implicit_fermi_first_order_response(H0_bml, H1_bml, P0_bml, P1_bml, Inv_bml, nsteps, mu0, beta, nocc, threshold)
Calculate first order density matrix response to perturbations using Implicit Fermi Dirac.
Definition: prg_implicit_fermi_mod.F90:441
prg_implicit_fermi_mod::prg_get_density_matrix
subroutine prg_get_density_matrix(ham_bml, p_bml, beta, mu, threshold)
Calculate the density matrix with diagonalization.
Definition: prg_implicit_fermi_mod.F90:1057
prg_implicit_fermi_mod::prg_implicit_fermi_response
subroutine, public prg_implicit_fermi_response(H0_bml, H1_bml, H2_bml, H3_bml, P0_bml, P1_bml, P2_bml, P3_bml, nsteps, mu0, mu, beta, nocc, occ_tol, lin_tol, order, threshold)
Calculate density matrix response to perturbations using Implicit Fermi Dirac.
Definition: prg_implicit_fermi_mod.F90:549
prg_implicit_fermi_mod::prg_setup_linsys
subroutine prg_setup_linsys(p_bml, A_bml, b_bml, p2_bml, y_bml, aux_bml, aux1_bml, k, threshold)
Set up linear system for Implicit Fermi Dirac.
Definition: prg_implicit_fermi_mod.F90:845
prg_implicit_fermi_mod::prg_finite_diff
subroutine, public prg_finite_diff(H0_bml, H_list, mu0, mu_list, beta, order, lambda, h, threshold)
Calculate density matrix response from perturbations in the Hamiltonian.
Definition: prg_implicit_fermi_mod.F90:727
prg_implicit_fermi_mod::prg_test_density_matrix
subroutine, public prg_test_density_matrix(ham_bml, p_bml, beta, mu, nocc, osteps, occErrLimit, threshold)
Calculate the density matrix with diagonalization and converge chemical.
Definition: prg_implicit_fermi_mod.F90:1113
prg_implicit_fermi_mod::prg_conjgrad
subroutine prg_conjgrad(A_bml, p_bml, p2_bml, tmp_bml, d_bml, w_bml, cg_tol, threshold)
Solve the system AX = B with conjugate gradient.
Definition: prg_implicit_fermi_mod.F90:1007
prg_implicit_fermi_mod::prg_newtonschulz
subroutine prg_newtonschulz(a_bml, ai_bml, r_bml, tmp_bml, d_bml, I_bml, tol, threshold, num_iter)
Find the inverse of the matrix A with Newton-Schulz iteration.
Definition: prg_implicit_fermi_mod.F90:889
prg_implicit_fermi_mod::prg_implicit_fermi_zero
subroutine, public prg_implicit_fermi_zero(h_bml, p_bml, nsteps, mu, method, threshold, tol)
Recursive Implicit Fermi Dirac for zero temperature.
Definition: prg_implicit_fermi_mod.F90:348
prg_normalize_mod
The prg_normalize module.
Definition: prg_normalize_mod.F90:5
prg_normalize_mod::prg_normalize_implicit_fermi
subroutine, public prg_normalize_implicit_fermi(h_bml, cnst, mu)
Normalize a Hamiltonian matrix prior to running the implicit fermi dirac algorithm.
Definition: prg_normalize_mod.F90:84
prg_parallel_mod
The parallel module.
Definition: prg_parallel_mod.F90:6
prg_timer_mod
The timer module.
Definition: prg_timer_mod.F90:23