PROGRESS  master
prg_xlbokernel_mod.F90
Go to the documentation of this file.
1 
5 module prg_xlbokernel_mod
6 
7  use omp_lib
8  use bml
9  use prg_normalize_mod
10  use prg_densitymatrix_mod
11  use prg_timer_mod
12  use prg_parallel_mod
13  use prg_ewald_mod
14  use prg_implicit_fermi_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_kernel_multirank
23  public :: prg_kernel_multirank_latte
24  public :: prg_kernel_matrix_multirank
25  public :: prg_full_kernel
26  public :: prg_full_kernel_latte
27 
28 contains
29 
30  subroutine invert(A,AI,N)
31 
32  implicit none
33  integer, parameter :: PREC = 8
34  integer, intent(in) :: N
35  real(PREC), intent(in) :: A(N,N)
36  real(PREC), intent(out) :: AI(N,N)
37  real(PREC) :: WORK(N+N*N)!, C(N,N)
38  integer :: LDA, LWORK, M, INFO, IPIV(N)
39  integer :: I,J,K
40 
41  external DGETRF
42  external DGETRI
43 
44  ai = a
45  lda = n
46  m = n
47  lwork = n+n*n
48 
49  call dgetrf(m, n, ai, lda, ipiv, info)
50  call dgetri(n, ai, n, ipiv, work, lwork, info)
51 
52  end subroutine invert
53 
84  subroutine prg_kernel_multirank_latte(KRes,KK0_bml,Res,FelTol,L,LMAX,NUMRANK,HO_bml,mu,beta,RXYZ,Box,Hubbard_U,Element_Pointer, &
85  Nr_atoms,HDIM,Max_Nr_Neigh,Coulomb_acc,nebcoul,totnebcoul,Hinxlist, &
86  S_bml,Z_bml,Nocc,Inv_bml,H1_bml,DO_bml,D1_bml,m_rec,threshold,Nr_elem,mdstep,update)
87 
88  !! Res = q[n] - n
89  !! KK0 is preconditioner
90  !! KRes rank-L approximation of (K0*J)^(-1)*K0*(q[n]-n) with (K0*J)^(-1) as in Eq. (41) in Ref. [*]
91 
92  implicit none
93  integer, intent(in) :: nr_atoms, hdim, max_nr_neigh,m_rec,nr_elem,numrank,mdstep,update
94  real(dp), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
95  real(dp), intent(in) :: coulomb_acc, feltol, beta, mu, nocc
96  real(dp), intent(in) :: rxyz(3,nr_atoms), box(3,3)
97  real(dp), intent(in) :: res(nr_atoms)
98  integer, intent(in) :: hinxlist(hdim),element_pointer(nr_atoms)
99  real(dp), intent(in) :: hubbard_u(nr_elem)
100  type(bml_matrix_t), intent(inout) :: ho_bml, s_bml, z_bml, inv_bml(m_rec), kk0_bml
101  real(dp) :: k0res(nr_atoms)
102  type(bml_matrix_t),intent(inout) :: h1_bml, do_bml,d1_bml
103  real(dp), intent(in) :: threshold
104  integer, intent(in) :: lmax
105  character(20) :: bml_type
106  integer, intent(in) :: totnebcoul(nr_atoms), nebcoul(4,max_nr_neigh,nr_atoms)
107  real(dp) :: coulomb_pot_real(nr_atoms), coulomb_pot(nr_atoms)
108  real(dp) :: coulomb_pot_real_i, coulomb_pot_k(nr_atoms),dq_dv(nr_atoms)
109  real(dp) :: coulomb_pot_real_dq_v(nr_atoms), coulomb_pot_dq_v(nr_atoms)
110  real(dp) :: dq_v(nr_atoms)
111  real(dp), intent(inout) :: kres(nr_atoms)
112  integer :: i,j,k,it,n,mn
113  integer, intent(out) :: l
114  real(dp) :: fel, proj_tmp, start, finish
115  real(dp) :: vi(nr_atoms,lmax), fi(nr_atoms,lmax), v(nr_atoms)
116  real(dp) :: dr(nr_atoms), identres(nr_atoms)
117  real(dp), allocatable :: o(:,:),m(:,:),row1(:),row2(:),row_na(:)
118  type(bml_matrix_t) :: kk0t_bml, k0res_bml, res_bml, zt_bml, t_bml, x_bml, y_bml
119 
120  call timer_prg_init()
121  bml_type = bml_get_type(ho_bml)
122  n = bml_get_n(ho_bml)
123  mn = bml_get_m(ho_bml)
124  allocate(row1(hdim)); allocate(row2(hdim)); allocate(row_na(nr_atoms))
125  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,kk0t_bml)
126  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,k0res_bml)
127  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,res_bml)
128  call bml_zero_matrix(bml_type,bml_element_real,dp,n,mn,zt_bml)
129  call bml_zero_matrix(bml_type,bml_element_real,dp,n,mn,x_bml)
130  call bml_zero_matrix(bml_type,bml_element_real,dp,n,mn,y_bml)
131 
132  call bml_transpose_new(z_bml,zt_bml)
133 
134  ! K0Res = KK0*Res temporary for matrix-vector multiplication
135  call bml_set_row(res_bml,1,res,1.0_dp*1e-10)
136  call bml_transpose_new(kk0_bml,kk0t_bml)
137  call bml_multiply(res_bml,kk0t_bml,k0res_bml,1.0_dp,0.0_dp,1.0_dp*1e-10)
138  call bml_get_row(k0res_bml,1,row_na)
139  k0res = row_na
140  dr = k0res
141 
142  write(*,*) 'resnorm =', norm2(res), 'kresnorm =', norm2(dr)
143  i = 0
144  fel = 1.d0
145  do while ((fel > feltol).AND.(i < (lmax)).AND.(i < numrank)) !! Fel = "Error" in Swedish
146  i = i + 1
147  !write(*,*) 'dr =', norm2(dr)
148  vi(:,i) = dr/norm2(dr)
149  do j = 1,i-1
150  vi(:,i) = vi(:,i) - dot_product(vi(:,i),vi(:,j))*vi(:,j) !! Orthogonalized v_i as in Eq. (42) Ref. [*]
151  enddo
152  vi(:,i) = vi(:,i)/norm2(vi(:,i))
153  v(:) = vi(:,i) ! v_i
154 
155  ! Compute H1 = H(v)
156  dq_v = v
157  call prg_timer_start(1)
158  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(J,Coulomb_Pot_Real_I)
159  do j = 1,nr_atoms
160  call ewald_real_space_latte(coulomb_pot_real_i,j,rxyz,box, &
161  dq_v,hubbard_u,element_pointer,nr_atoms,coulomb_acc,nebcoul,totnebcoul,hdim,max_nr_neigh,nr_elem)
162  coulomb_pot_real(j) = coulomb_pot_real_i
163  enddo
164  !$OMP END PARALLEL DO
165 
166  call ewald_k_space_latte(coulomb_pot_k,rxyz,box,dq_v,nr_atoms,coulomb_acc,max_nr_neigh)
167  coulomb_pot_dq_v = coulomb_pot_real+coulomb_pot_k
168  call prg_timer_stop(1,1)
169 
170 
171  call bml_deallocate(h1_bml)
172  call bml_zero_matrix(bml_type,bml_element_real,dp,hdim,mn,h1_bml)
173 
174  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(It,K)
175  do it = 1,nr_atoms-1
176  do k = hinxlist(it)+1,hinxlist(it+1)
177  row1(k) = hubbard_u(element_pointer(it))*dq_v(it) + coulomb_pot_dq_v(it)
178  enddo
179  enddo
180  !$OMP END PARALLEL DO
181  do k = hinxlist(nr_atoms)+1,hdim
182  row1(k) = hubbard_u(element_pointer(nr_atoms))*dq_v(nr_atoms) + coulomb_pot_dq_v(nr_atoms)
183  enddo
184 
185  ! H1 = 1/2(S*H1+H1*S)
186  call bml_set_diagonal(h1_bml,row1,threshold)
187  call bml_multiply(s_bml,h1_bml,x_bml,1.0_dp,0.0_dp,threshold)
188  call bml_multiply(h1_bml,s_bml,x_bml,0.5_dp,0.5_dp,threshold)
189  call bml_copy(x_bml,h1_bml)
190 
191  ! H1 = Z^T H1 Z
192  call bml_multiply(zt_bml,h1_bml,x_bml,1.0_dp,0.0_dp,threshold)
193  call bml_multiply(x_bml,z_bml,h1_bml,1.0_dp,0.0_dp,threshold)
194 
195  ! Compute D1 = F_FD(HO_bml + eps*H1_bml)/eps at eps = 0
196  call prg_implicit_fermi_first_order_response(ho_bml,h1_bml,do_bml,d1_bml,inv_bml, &
197  m_rec, mu, beta, real(nocc,dp), threshold)
198 
199  ! D1 = Z D1 Z^T
200  call bml_multiply(z_bml,d1_bml,y_bml,2.0_dp,0.0_dp,threshold)
201  call bml_multiply(y_bml,zt_bml,d1_bml,1.0_dp,0.0_dp,threshold)
202 
203  ! Compute dq/dv
204  call bml_multiply(d1_bml,s_bml,x_bml, 1.0_dp,0.0_dp,threshold)
205  call bml_get_diagonal(x_bml,row1)
206  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(It,K)
207  do it = 1, nr_atoms-1
208  dq_dv(it) = 0
209  do k = hinxlist(it)+1,hinxlist(it+1)
210  dq_dv(it) = dq_dv(it) + row1(k)
211  enddo
212  enddo
213  !$OMP END PARALLEL DO
214  dq_dv(nr_atoms) = 0
215  do k = hinxlist(nr_atoms)+1,hdim
216  dq_dv(nr_atoms) = dq_dv(nr_atoms) + row1(k)
217  enddo
218 
219  dr = dq_dv - v
220  ! fi = K0(dq_dv - v)
221  call bml_set_row(res_bml,1,dr,1.0_dp*1e-10)
222  call bml_multiply(res_bml,kk0t_bml,k0res_bml,1.0_dp,0.0_dp,1.0_dp*1e-10)
223  call bml_get_row(k0res_bml,1,row_na)
224  dr = row_na
225  fi(:,i) = dr
226 
227  l = i
228  allocate(o(l,l), m(l,l))
229  do k = 1,l
230  do j = 1,l
231  o(k,j) = dot_product(fi(:,k),fi(:,j)) ! O_KJ = < fv_i(K) | fv_i(J) > see below Eq. (31)
232  enddo
233  enddo
234  call invert(o,m,l) ! M = O^(-1)
235  identres = 0.d0*k0res
236  kres = 0.d0
237  do k = 1,l
238  do j = 1,l
239  proj_tmp = m(k,j)*dot_product(fi(:,j),k0res)
240  identres = identres + proj_tmp*fi(:,k)
241  kres = kres + proj_tmp*vi(:,k) !! KRes becomes the rank-L approximate of (K0*J)^(-1)*K0*(q[n]-n)
242  enddo
243  enddo
244  fel = norm2(identres-k0res)/norm2(identres) !! RELATIVE RESIDUAL ERROR ESTIMATE Eq. (48) Ref. [*]
245  write(*,*) '# I, L, Fel = ',i,l,fel !! Fel goes down with L
246  deallocate(o, m)
247 
248  enddo
249 
250  if (modulo(mdstep,update) .eq. 0) then
251  call prg_update_preconditioner(kk0_bml,vi(:,1),fi(:,1),nr_atoms,threshold)
252  endif
253 
254  deallocate(row1);deallocate(row2);deallocate(row_na)
255  call bml_deallocate(kk0t_bml)
256  call bml_deallocate(k0res_bml)
257  call bml_deallocate(res_bml)
258  call bml_deallocate(zt_bml)
259  call bml_deallocate(x_bml)
260  call bml_deallocate(y_bml)
261  call prg_timer_shutdown()
262 
263  end subroutine prg_kernel_multirank_latte
264 
265  subroutine prg_update_preconditioner(K0,v,fv,Nr_atoms,threshold)
266 
267  implicit none
268  type(bml_matrix_t), intent(inout) :: K0
269  real(dp), intent(in) :: v(Nr_atoms),fv(Nr_atoms),threshold
270  integer, intent(in) :: Nr_atoms
271  type(bml_matrix_t) :: v_bml,fv_bml,vt_bml,tmp1_bml,tmp2_bml
272  type(bml_matrix_t) :: ones_bml,onest_bml,K0_update,K0_T
273  real(dp) :: const,kthresh
274  real(dp) :: ones(Nr_atoms)
275  real(dp),allocatable :: row(:),norm(:)
276  integer :: I
277  character(20) :: bml_type
278 
279  bml_type = bml_get_type(k0)
280  kthresh = 1.0_dp*1e-5
281  allocate(row(nr_atoms))
282  allocate(norm(nr_atoms))
283 
284  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,v_bml)
285  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,vt_bml)
286  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,fv_bml)
287  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,tmp1_bml)
288  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,tmp2_bml)
289  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,ones_bml)
290  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,onest_bml)
291  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,k0_update)
292  call bml_zero_matrix(bml_type,bml_element_real,dp,nr_atoms,nr_atoms,k0_t)
293 
294  ones = 1.0_dp
295  call bml_set_row(onest_bml,1,ones,0.0_dp)
296  call bml_multiply(k0,ones_bml,tmp1_bml,1.0_dp,0.0_dp,threshold)
297  const = 1.0_dp/(bml_trace_mult(onest_bml,tmp1_bml) + 1)
298 
299  call bml_set_row(vt_bml,1,v,1.0_dp*1e-10)
300  call bml_set_row(fv_bml,1,fv,1.0_dp*1e-10)
301  call bml_transpose_new(vt_bml,v_bml)
302 
303  call bml_multiply(k0,v_bml,tmp1_bml,1.0_dp,0.0_dp,threshold)
304  call bml_copy(vt_bml,tmp2_bml)
305  call bml_multiply(fv_bml,k0,tmp2_bml,1.0_dp,-1.0_dp,threshold)
306  call bml_multiply(tmp1_bml,tmp2_bml,k0_update,const,0.0_dp,kthresh)
307  call bml_add(k0,k0_update,1.0_dp,-1.0_dp,kthresh)
308 
309  call bml_multiply(onest_bml,k0,tmp1_bml,1.0_dp,0.0_dp,0.0_dp)
310  call bml_get_row(tmp1_bml,1,norm)
311  call bml_transpose_new(k0,k0_t)
312  do i = 1,nr_atoms
313  call bml_get_row(k0_t,1,row)
314  call bml_set_row(k0_t,1,row/norm(i),1.0_dp*1e-10)
315  enddo
316  call bml_transpose_new(k0_t,k0)
317 
318  deallocate(row)
319  deallocate(norm)
320  call bml_deallocate(v_bml)
321  call bml_deallocate(vt_bml)
322  call bml_deallocate(fv_bml)
323  call bml_deallocate(tmp1_bml)
324  call bml_deallocate(tmp2_bml)
325  call bml_deallocate(onest_bml)
326  call bml_deallocate(ones_bml)
327  call bml_deallocate(k0_update)
328  call bml_deallocate(k0_t)
329 
330  end subroutine prg_update_preconditioner
331 
332 
333  ! Above routine but for development code
334  subroutine prg_kernel_multirank(KRes,KK0_bml,Res,FelTol,L,LMAX,HO_bml,mu,beta,RX,RY,RZ,LBox,Hubbard_U,Element_Type, &
335  Nr_atoms,HDIM,Max_Nr_Neigh,Coulomb_acc,TIMERATIO,nnRx,nnRy,nnRz,nrnnlist,nnType,H_INDEX_START,H_INDEX_END, &
336  S_bml,Z_bml,Nocc,Znuc,Inv_bml,H1_bml,X_bml,Y_bml,DO_bml,D1_bml,m_rec,threshold)
337 
338  !! Res = q[n] - n
339  !! KK0 is preconditioner
340  !! KRes rank-L approximation of (K0*J)^(-1)*K0*(q[n]-n) with (K0*J)^(-1) as in Eq. (41) in Ref. [*]
341 
342  implicit none
343  integer, intent(in) :: nr_atoms, hdim, nocc,max_nr_neigh,m_rec
344  real(dp), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
345  real(dp), intent(in) :: coulomb_acc, timeratio, feltol, beta, mu
346  real(dp), intent(in) :: rx(nr_atoms), ry(nr_atoms), rz(nr_atoms), lbox(3)
347  real(dp), intent(in) :: res(nr_atoms)
348  integer, intent(in) :: h_index_start(nr_atoms), h_index_end(nr_atoms)
349  real(dp), intent(in) :: znuc(nr_atoms), hubbard_u(nr_atoms)
350  type(bml_matrix_t), intent(in) :: ho_bml, s_bml, z_bml, inv_bml(m_rec), kk0_bml
351  real(dp) :: k0res(nr_atoms)
352  type(bml_matrix_t),intent(inout) :: h1_bml, x_bml,y_bml,do_bml,d1_bml
353  real(dp), intent(in) :: threshold
354  integer, intent(in) :: lmax
355  character(10), intent(in) :: element_type(nr_atoms)
356  integer, intent(in) :: nrnnlist(nr_atoms), nntype(nr_atoms,max_nr_neigh)
357  real(dp), intent(in) :: nnrx(nr_atoms,max_nr_neigh)
358  real(dp), intent(in) :: nnry(nr_atoms,max_nr_neigh), nnrz(nr_atoms,max_nr_neigh)
359  real(dp) :: coulomb_pot_real(nr_atoms), coulomb_pot(nr_atoms)
360  real(dp) :: coulomb_pot_real_i, coulomb_pot_k(nr_atoms),dq_dv(nr_atoms)
361  real(dp) :: coulomb_pot_real_dq_v(nr_atoms), coulomb_pot_dq_v(nr_atoms)
362  real(dp) :: coulomb_force_real_i(3),coulomb_force_k(3,nr_atoms)
363  real(dp) :: dq_v(nr_atoms)
364  real(dp), intent(inout) :: kres(nr_atoms)
365  integer :: i,j,k,it
366  integer, intent(out) :: l
367  real(dp) :: fel, proj_tmp
368  real(dp) :: vi(nr_atoms,lmax), fi(nr_atoms,lmax), v(nr_atoms)
369  real(dp) :: dr(nr_atoms), identres(nr_atoms)
370  real(dp), allocatable :: o(:,:),m(:,:),row1(:),row2(:),row_na(:)
371  type(bml_matrix_t) :: kk0t_bml, k0res_bml, res_bml
372 
373  allocate(row1(hdim)); allocate(row2(hdim)); allocate(row_na(nr_atoms))
374  call bml_zero_matrix("ellpack",bml_element_real,dp,nr_atoms,nr_atoms,kk0t_bml)
375  call bml_zero_matrix("ellpack",bml_element_real,dp,nr_atoms,nr_atoms,k0res_bml)
376  call bml_zero_matrix("ellpack",bml_element_real,dp,nr_atoms,nr_atoms,res_bml)
377 
378  ! K0Res = KK0*Res temporary for matrix-vector multiplication
379  call bml_transpose_new(kk0_bml,kk0t_bml)
380  call bml_set_row(res_bml,1,res,one*1e-14)
381  call bml_multiply(res_bml,kk0t_bml,k0res_bml,1.0_dp,0.0_dp,one*1e-14)
382  call bml_get_row(k0res_bml,1,row_na)
383  k0res = row_na
384  dr = k0res
385 
386  i = 0
387  fel = 1.d0
388  do while ((fel > feltol).AND.(i < (lmax))) !! Fel = "Error" in Swedish
389  i = i + 1
390  vi(:,i) = dr/norm2(dr)
391  do j = 1,i-1
392  vi(:,i) = vi(:,i) - dot_product(vi(:,i),vi(:,j))*vi(:,j) !! Orthogonalized v_i as in Eq. (42) Ref. [*]
393  enddo
394  vi(:,i) = vi(:,i)/norm2(vi(:,i))
395  v(:) = vi(:,i) ! v_i
396 
397  dq_v = v
398 
399  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(J,Coulomb_Pot_Real_I,Coulomb_Force_Real_I)
400  do j = 1,nr_atoms
401  call ewald_real_space(coulomb_pot_real_i,coulomb_force_real_i,j,rx,ry,rz,lbox, &
402  dq_v,hubbard_u,element_type,nr_atoms,coulomb_acc,timeratio,nnrx,nnry,nnrz,nrnnlist,nntype,hdim,max_nr_neigh)
403  coulomb_pot_real(j) = coulomb_pot_real_i
404  enddo
405  !$OMP END PARALLEL DO
406 
407  ! call Ewald_k_Space(Coulomb_Pot_k,Coulomb_Force_k,RX,RY,RZ,LBox,dq_v,Nr_atoms,Coulomb_acc,TIMERATIO,Max_Nr_Neigh)
408  ! Coulomb_Pot_dq_v = Coulomb_Pot_Real+Coulomb_Pot_k
409 
410  row1 = 0.0_dp
411  do j = 1,hdim
412  call bml_set_row(h1_bml,j,row1,threshold)
413  enddo
414 
415  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(J,K)
416  do j = 1,nr_atoms
417  do k = h_index_start(j),h_index_end(j)
418  row1(k) = hubbard_u(j)*dq_v(j) + coulomb_pot_dq_v(j)
419  enddo
420  enddo
421  !$OMP END PARALLEL DO
422 
423  call bml_set_diagonal(h1_bml,row1,threshold)
424  call bml_multiply(s_bml,h1_bml,d1_bml,1.0_dp,0.0_dp,threshold)
425  call bml_multiply(h1_bml,s_bml,d1_bml,0.5_dp,0.5_dp,threshold)
426  call bml_copy(d1_bml,h1_bml)
427 
428  call bml_transpose_new(z_bml,x_bml)
429  call bml_multiply(x_bml,h1_bml,d1_bml,1.0_dp,0.0_dp,threshold)
430  call bml_multiply(d1_bml,z_bml,h1_bml,1.0_dp,0.0_dp,threshold)
431 
432  call prg_implicit_fermi_first_order_response(ho_bml,h1_bml,do_bml,d1_bml,inv_bml, &
433  m_rec, mu, beta, real(nocc,dp), threshold)
434 
435  call bml_transpose_new(z_bml,x_bml)
436  call bml_multiply(z_bml,d1_bml,y_bml,2.0_dp,0.0_dp,threshold)
437  call bml_multiply(y_bml,x_bml,d1_bml,1.0_dp,0.0_dp,threshold)
438 
439  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(It,K,row1,row2)
440  do it = 1, nr_atoms
441  dq_dv(it) = 0
442  do k = h_index_start(it), h_index_end(it)
443  call bml_get_row(s_bml,k,row1)
444  call bml_get_row(d1_bml,k,row2)
445  dq_dv(it) = dq_dv(it) + dot_product(row1,row2)
446  enddo
447  enddo
448  !$OMP END PARALLEL DO
449 
450  dr = dq_dv - v
451  ! fi = K0(dq_dv - v)
452  call bml_set_row(res_bml,1,dr,threshold)
453  call bml_multiply(res_bml,kk0t_bml,k0res_bml,1.0_dp,0.0_dp,threshold)
454  call bml_get_row(k0res_bml,1,row_na)
455  dr = row_na
456  fi(:,i) = dr
457 
458  l = i
459  allocate(o(l,l), m(l,l))
460  do k = 1,l
461  do j = 1,l
462  o(k,j) = dot_product(fi(:,k),fi(:,j)) ! O_KJ = < fv_i(K) | fv_i(J) > see below Eq. (31)
463  enddo
464  enddo
465  call invert(o,m,l) ! M = O^(-1)
466  identres = 0.d0*k0res
467  kres = 0.d0
468  do k = 1,l
469  do j = 1,l
470  proj_tmp = m(k,j)*dot_product(fi(:,j),k0res)
471  identres = identres + proj_tmp*fi(:,k)
472  kres = kres + proj_tmp*vi(:,k) !! KRes becomes the rank-L approximate of (K0*J)^(-1)*K0*(q[n]-n)
473  enddo
474  enddo
475  fel = norm2(identres-k0res)/norm2(identres) !! RELATIVE RESIDUAL ERROR ESTIMATE Eq. (48) Ref. [*]
476  write(*,*) '# I, L, Fel = ',i,l,fel !! Fel goes down with L
477  deallocate(o, m)
478 
479  enddo
480 
481  deallocate(row1);deallocate(row2);deallocate(row_na)
482  call bml_deallocate(kk0t_bml)
483  call bml_deallocate(k0res_bml)
484  call bml_deallocate(res_bml)
485 
486  end subroutine prg_kernel_multirank
487 
511  subroutine prg_full_kernel_latte(KK,DO_bml,mu0,RXYZ,Box,Hubbard_U, &
512  Element_Pointer, Nr_atoms,HDIM, Max_Nr_Neigh,Coulomb_acc, &
513  Hinxlist,S_bml,Z_bml,Inv_bml,D1_bml,H1_bml,HO_bml, &
514  Nocc,m_rec,threshold,beta,Nr_elem,nebcoul,totnebcoul)
515 
516  use bml
517 
518  implicit none
519  integer, parameter :: prec = 8, dp = kind(1.0d0)
520  integer, intent(in) :: nr_atoms, nr_elem,hdim, max_nr_neigh,m_rec
521  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
522  real(prec), parameter :: kb = 8.61739d-5 ! eV/K, kB = 6.33366256e-6 Ry/K, kB = 3.166811429e-6 Ha/K
523  real(prec), intent(in) :: coulomb_acc, threshold,beta,nocc
524  real(prec) :: v(nr_atoms)
525  real(prec), intent(in) :: rxyz(3,nr_atoms),box(3,3)
526  integer, intent(in) :: hinxlist(hdim)
527  real(prec), intent(in) :: hubbard_u(nr_elem)
528  type(bml_matrix_t), intent(in) :: inv_bml(m_rec)
529  type(bml_matrix_t), intent(in) :: s_bml,z_bml,ho_bml
530  type(bml_matrix_t),intent(inout) :: h1_bml,do_bml,d1_bml
531  real(prec), intent(inout) :: mu0
532  integer, intent(in) :: element_pointer(nr_atoms),nebcoul(4,max_nr_neigh,nr_atoms),totnebcoul(nr_atoms)
533  real(prec) :: coulomb_pot_real(nr_atoms), coulomb_pot(nr_atoms)
534  real(prec) :: coulomb_pot_real_i, coulomb_pot_k(nr_atoms)
535  real(prec) :: coulomb_pot_real_dq_v(nr_atoms), coulomb_pot_dq_v(nr_atoms)
536  real(prec) :: dq_v(nr_atoms)
537  real(prec) :: dq_dv(nr_atoms), err,tol,start,ewaldk,ewaldr,ewaldkacc,ewaldracc,response,respacc
538  real(prec), intent(out) :: kk(nr_atoms,nr_atoms)
539  integer :: i,j,k, iter, mm,it,n,mn
540  real(prec),allocatable :: row1(:),row2(:),jj(:,:),e(:,:),ereal(:,:),ekspace(:,:)
541  type(bml_matrix_t) :: zt_bml, r_bml, jji_bml, jj_bml, tmp_bml, x_bml, y_bml
542  character(20) :: bml_type
543 
544  bml_type = bml_get_type(ho_bml)
545  n = bml_get_n(ho_bml)
546  mn = bml_get_m(ho_bml)
547  call bml_zero_matrix(bml_type,bml_element_real,dp,n,mn,zt_bml)
548  call bml_zero_matrix(bml_type,bml_element_real,dp,n,mn,x_bml)
549  call bml_zero_matrix(bml_type,bml_element_real,dp,n,mn,y_bml)
550  call bml_transpose_new(z_bml,zt_bml)
551  allocate(row1(hdim)); allocate(row2(hdim)); allocate(jj(nr_atoms,nr_atoms))
552  allocate(ereal(nr_atoms,nr_atoms)); allocate(ekspace(nr_atoms,nr_atoms))
553  allocate(e(nr_atoms,nr_atoms))
554  coulomb_pot_dq_v = zero
555  coulomb_pot_k = zero
556  dq_v = zero
557  jj = zero
558  kk = zero
559 
560  !call Ewald_Real_Space_Matrix_latte(EReal,RXYZ,Box,Hubbard_U,Element_Pointer, &
561  ! Nr_atoms,Coulomb_acc,Nebcoul,Totnebcoul,HDIM,Max_Nr_Neigh,Nr_Elem)
562  !call Ewald_k_Space_Matrix_latte(EKspace,RXYZ,Box,Nr_atoms,Coulomb_acc,Max_Nr_Neigh,Nebcoul,Totnebcoul)
563  !E = EReal + EKspace
564  write(*,*) 'Beginning response calculations'
565  ewaldracc = 0.0; ewaldkacc = 0.0; respacc = 0.0;
566  do j = 1,nr_atoms
567  call cpu_time(start)
568  write(*,*) 'J=',j
569  dq_v(j) = one
570  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,Coulomb_Pot_Real_I)
571  do i = 1,nr_atoms
572  call ewald_real_space_single_latte(coulomb_pot_real_i,i,rxyz,box,nr_elem, &
573  dq_v,j,hubbard_u,element_pointer,nr_atoms,coulomb_acc,hdim,max_nr_neigh)
574  coulomb_pot_real(i) = coulomb_pot_real_i
575  enddo
576  !$OMP END PARALLEL DO
577  call cpu_time(ewaldr)
578  ewaldracc = ewaldracc + ewaldr-start
579  call ewald_k_space_latte_single(coulomb_pot_k,j,rxyz,box,dq_v,nr_atoms,coulomb_acc)
580  coulomb_pot_dq_v = coulomb_pot_real+coulomb_pot_k
581  call cpu_time(ewaldk)
582  ewaldkacc = ewaldkacc + ewaldk-ewaldr
583 
584 
585  call bml_deallocate(h1_bml)
586  call bml_zero_matrix(bml_type,bml_element_real,dp,n,mn,h1_bml)
587 
588  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,K)
589  do i = 1,nr_atoms-1
590  do k = hinxlist(i)+1,hinxlist(i+1)
591  row1(k) = hubbard_u(element_pointer(i))*dq_v(i) + coulomb_pot_dq_v(i)
592  enddo
593  enddo
594  !$OMP END PARALLEL DO
595  do k = hinxlist(nr_atoms)+1,hdim
596  row1(k) = hubbard_u(element_pointer(nr_atoms))*dq_v(nr_atoms) + coulomb_pot_dq_v(nr_atoms)
597  enddo
598 
599  call bml_set_diagonal(h1_bml,row1,threshold)
600  call bml_multiply(s_bml,h1_bml,x_bml,1.0_dp,0.0_dp,threshold)
601  call bml_multiply(h1_bml,s_bml,x_bml,0.5_dp,0.5_dp,threshold)
602  call bml_copy(x_bml,h1_bml)
603 
604  call bml_multiply(zt_bml,h1_bml,x_bml,1.0_dp,0.0_dp,threshold)
605  call bml_multiply(x_bml,z_bml,h1_bml,1.0_dp,0.0_dp,threshold)
606 
607  call prg_implicit_fermi_first_order_response(ho_bml,h1_bml,do_bml,d1_bml,inv_bml, &
608  m_rec,mu0,beta,nocc,threshold)
609  call cpu_time(response)
610  respacc = respacc + response-ewaldk
611 
612  call bml_multiply(z_bml,d1_bml,y_bml,2.0_dp,0.0_dp,threshold)
613  call bml_multiply(y_bml,zt_bml,d1_bml,1.0_dp,0.0_dp,threshold)
614 
615  call bml_multiply(d1_bml,s_bml,x_bml,1.0_dp,0.0_dp,threshold)
616  call bml_get_diagonal(x_bml,row1)
617  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(It,K)
618  do it = 1, nr_atoms-1
619  dq_dv(it) = 0
620  do k = hinxlist(it)+1,hinxlist(it+1)
621  dq_dv(it) = dq_dv(it) + row1(k)
622  enddo
623  jj(it,j) = dq_dv(it)
624  enddo
625  !$OMP END PARALLEL DO
626  dq_dv(nr_atoms) = 0
627  do k = hinxlist(nr_atoms)+1,hdim
628  dq_dv(nr_atoms) = dq_dv(nr_atoms) + row1(k)
629  enddo
630  jj(nr_atoms,j) = dq_dv(nr_atoms)
631  dq_v = zero
632  enddo
633 
634  do i = 1,nr_atoms
635  jj(i,i) = jj(i,i) - one
636  enddo
637 
638  write(*,*) 'ewaldracc =', ewaldracc
639  write(*,*) 'ewaldkacc =', ewaldkacc
640  write(*,*) 'implcit response time =', respacc
641  call invert(jj,kk,nr_atoms)
642 
643  deallocate(row1); deallocate(row2); deallocate(jj)
644  deallocate(e); deallocate(ereal); deallocate(ekspace)
645  call bml_deallocate(zt_bml)
646  call bml_deallocate(x_bml)
647  call bml_deallocate(y_bml)
648 
649  end subroutine prg_full_kernel_latte
650 
679 
680  subroutine prg_full_kernel(KK,DO_bml,mu0,RX,RY,RZ,LBox,Hubbard_U,Element_Type, &
681  Nr_atoms,HDIM, Max_Nr_Neigh,Coulomb_acc,TIMERATIO,nnRx,nnRy,nnRz, &
682  nrnnlist,nnType,H_INDEX_START,H_INDEX_END,S_bml,Z_bml,Inv_bml,D1_bml,H1_bml,HO_bml,Y_bml,X_bml, &
683  Nocc,Znuc,m_rec,threshold,beta,diagonal)
684 
685  use bml
686 
687  implicit none
688  integer, parameter :: prec = 8, dp = kind(1.0d0)
689  integer, intent(in) :: nr_atoms, hdim, nocc, max_nr_neigh,m_rec
690  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
691  real(prec), parameter :: kb = 8.61739d-5 ! eV/K, kB = 6.33366256e-6 Ry/K, kB = 3.166811429e-6 Ha/K
692  real(prec), intent(in) :: coulomb_acc, timeratio,threshold,beta
693  real(prec) :: v(nr_atoms)
694  real(prec), intent(in) :: rx(nr_atoms), ry(nr_atoms), rz(nr_atoms), lbox(3)
695  integer, intent(in) :: h_index_start(nr_atoms), h_index_end(nr_atoms)
696  real(prec), intent(in) :: znuc(nr_atoms), hubbard_u(nr_atoms)
697  type(bml_matrix_t), intent(in) :: inv_bml(m_rec)
698  type(bml_matrix_t), intent(in) :: s_bml,z_bml,ho_bml
699  type(bml_matrix_t),intent(inout) :: h1_bml,do_bml,d1_bml,y_bml,x_bml
700  real(prec), intent(inout) :: mu0, diagonal(hdim)
701  character(10), intent(in) :: element_type(nr_atoms)
702  integer, intent(in) :: nrnnlist(nr_atoms), nntype(nr_atoms,max_nr_neigh)
703  real(prec), intent(in) :: nnrx(nr_atoms,max_nr_neigh)
704  real(prec), intent(in) :: nnry(nr_atoms,max_nr_neigh), nnrz(nr_atoms,max_nr_neigh)
705  real(prec) :: coulomb_pot_real(nr_atoms), coulomb_pot(nr_atoms)
706  real(prec) :: coulomb_pot_real_i, coulomb_pot_k(nr_atoms)
707  real(prec) :: coulomb_pot_real_dq_v(nr_atoms), coulomb_pot_dq_v(nr_atoms)
708  real(prec) :: coulomb_force_real_i(3),coulomb_force_k(3,nr_atoms)
709  real(prec) :: dq_v(nr_atoms)
710  real(prec) :: dq_dv(nr_atoms)
711  real(prec), intent(out) :: kk(nr_atoms,nr_atoms)
712  integer :: i,j,k, iter, mm,it
713  real(prec),allocatable :: row1(:),row2(:),jj(:,:)
714 
715  allocate(row1(hdim)); allocate(row2(hdim)); allocate(jj(nr_atoms,nr_atoms))
716  coulomb_pot_dq_v = zero
717  coulomb_pot_k = zero
718  dq_v = zero
719  jj = zero
720  kk = zero
721 
722  do j = 1,nr_atoms
723  dq_v(j) = one
724  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,Coulomb_Pot_Real_I)
725  do i = 1,nr_atoms
726  call ewald_real_space_single(coulomb_pot_real_i,coulomb_force_real_i,i,rx,ry,rz,lbox, &
727  dq_v,j,hubbard_u,element_type,nr_atoms,coulomb_acc,timeratio,hdim,max_nr_neigh)
728  coulomb_pot_real(i) = coulomb_pot_real_i
729  enddo
730  !$OMP END PARALLEL DO
731  ! call Ewald_k_Space(Coulomb_Pot_k,Coulomb_Force_k,RX,RY,RZ,LBox,dq_v,Nr_atoms,Coulomb_acc, &
732  ! TIMERATIO,Max_Nr_Neigh)
733  coulomb_pot_dq_v = coulomb_pot_real+coulomb_pot_k
734 
735  diagonal = 0.0_dp
736  do i = 1,hdim
737  call bml_set_row(h1_bml,i,diagonal,threshold)
738  enddo
739 
740  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,K)
741  do i = 1,nr_atoms
742  do k = h_index_start(i),h_index_end(i)
743  diagonal(k) = hubbard_u(i)*dq_v(i) + coulomb_pot_dq_v(i)
744  enddo
745  enddo
746  !$OMP END PARALLEL DO
747 
748  call bml_set_diagonal(h1_bml,diagonal,threshold)
749  call bml_multiply(s_bml,h1_bml,d1_bml,1.0_dp,0.0_dp,threshold)
750  call bml_multiply(h1_bml,s_bml,d1_bml,0.5_dp,0.5_dp,threshold)
751  call bml_copy(d1_bml,h1_bml)
752 
753  call bml_transpose_new(z_bml,x_bml)
754  call bml_multiply(x_bml,h1_bml,d1_bml,1.0_dp,0.0_dp,threshold)
755  call bml_multiply(d1_bml,z_bml,h1_bml,1.0_dp,0.0_dp,threshold)
756 
757  call prg_implicit_fermi_first_order_response(ho_bml,h1_bml,do_bml,d1_bml,inv_bml, &
758  m_rec,mu0,beta,real(nocc,prec),threshold)
759 
760  call bml_transpose_new(z_bml,x_bml)
761  call bml_multiply(z_bml,d1_bml,y_bml,2.0_dp,0.0_dp,threshold)
762  call bml_multiply(y_bml,x_bml,d1_bml,1.0_dp,0.0_dp,threshold)
763 
764  do it = 1, nr_atoms
765  dq_dv(it) = 0
766  do k = h_index_start(it), h_index_end(it)
767  call bml_get_row(s_bml,k,row1)
768  call bml_get_row(d1_bml,k,row2)
769  dq_dv(it) = dq_dv(it) + dot_product(row1,row2)
770  enddo
771  jj(it,j) = dq_dv(it)
772  enddo
773  dq_v = zero
774  enddo
775 
776  do i = 1,nr_atoms
777  jj(i,i) = jj(i,i) - one
778  enddo
779  call invert(jj,kk,nr_atoms)
780  deallocate(row1); deallocate(row2); deallocate(jj)
781 
782  end subroutine prg_full_kernel
783 
784 
785  ! Compute the low-rank kernel matrix. (For development code)
786  subroutine prg_kernel_matrix_multirank(KRes,KK0_bml,Res,FelTol,L,LMAX,HO_bml,mu,beta,RX,RY,RZ,LBox,Hubbard_U,Element_Type, &
787  Nr_atoms,HDIM,Max_Nr_Neigh,Coulomb_acc,TIMERATIO,nnRx,nnRy,nnRz,nrnnlist,nnType,H_INDEX_START,H_INDEX_END, &
788  S_bml,Z_bml,Nocc,Znuc,Inv_bml,H1_bml,X_bml,Y_bml,DO_bml,D1_bml,m_rec,threshold)
789 
790  !! Res = q[n] - n
791  !! KK0 is preconditioner
792  !! KRes rank-L approximation of (K0*J)^(-1)*K0*(q[n]-n) with (K0*J)^(-1) as in Eq. (41) in Ref. [*]
793 
794  implicit none
795  integer, intent(in) :: nr_atoms, hdim, max_nr_neigh,m_rec
796  real(dp), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
797  real(dp), intent(in) :: coulomb_acc, timeratio, feltol, beta, mu,nocc
798  real(dp), intent(in) :: rx(nr_atoms), ry(nr_atoms), rz(nr_atoms), lbox(3)
799  real(dp), intent(in) :: res(nr_atoms)
800  integer, intent(in) :: h_index_start(nr_atoms), h_index_end(nr_atoms)
801  real(dp), intent(in) :: znuc(nr_atoms), hubbard_u(nr_atoms)
802  type(bml_matrix_t), intent(in) :: ho_bml, s_bml, z_bml, inv_bml(m_rec)
803  real(dp) :: k0res(nr_atoms)
804  type(bml_matrix_t),intent(inout) :: h1_bml, x_bml,y_bml,do_bml,d1_bml
805  type(bml_matrix_t),intent(inout) :: kk0_bml
806  real(dp), intent(in) :: threshold
807  integer, intent(in) :: lmax
808  character(10), intent(in) :: element_type(nr_atoms)
809  integer, intent(in) :: nrnnlist(nr_atoms), nntype(nr_atoms,max_nr_neigh)
810  real(dp), intent(in) :: nnrx(nr_atoms,max_nr_neigh)
811  real(dp), intent(in) :: nnry(nr_atoms,max_nr_neigh), nnrz(nr_atoms,max_nr_neigh)
812  real(dp) :: coulomb_pot_real(nr_atoms), coulomb_pot(nr_atoms)
813  real(dp) :: coulomb_pot_real_i, coulomb_pot_k(nr_atoms),dq_dv(nr_atoms)
814  real(dp) :: coulomb_pot_real_dq_v(nr_atoms), coulomb_pot_dq_v(nr_atoms)
815  real(dp) :: coulomb_force_real_i(3),coulomb_force_k(3,nr_atoms)
816  real(dp) :: dq_v(nr_atoms)
817  real(dp), intent(out) :: kres(nr_atoms)
818  integer :: i,j,k,it,col,row
819  integer, intent(out) :: l
820  real(dp) :: fel, proj_tmp,elem
821  real(dp) :: vi(nr_atoms,lmax), fi(nr_atoms,lmax), v(nr_atoms)
822  real(dp) :: dr(nr_atoms), identres(nr_atoms)
823  real(dp), allocatable :: o(:,:),m(:,:),row1(:),row2(:)
824 
825  allocate(row1(hdim));allocate(row2(hdim));
826 
827  dr = res
828  i = 0
829  fel = 1.d0
830  do while ((fel > feltol).AND.(i < (lmax))) !! Fel = "Error" in Swedish
831  i = i + 1
832  vi(:,i) = dr/norm2(dr)
833  do j = 1,i-1
834  vi(:,i) = vi(:,i) - dot_product(vi(:,i),vi(:,j))*vi(:,j) !! Orthogonalized v_i as in Eq. (42) Ref. [*]
835  enddo
836  vi(:,i) = vi(:,i)/norm2(vi(:,i))
837  v(:) = vi(:,i) ! v_i
838 
839  dq_v = v
840 
841  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,Coulomb_Pot_Real_I,Coulomb_Force_Real_I)
842  do j = 1,nr_atoms
843  call ewald_real_space(coulomb_pot_real_i,coulomb_force_real_i,j,rx,ry,rz,lbox, &
844  dq_v,hubbard_u,element_type,nr_atoms,coulomb_acc,timeratio,nnrx,nnry,nnrz,nrnnlist,nntype,hdim,max_nr_neigh)
845  coulomb_pot_real(j) = coulomb_pot_real_i
846  enddo
847  !$OMP END PARALLEL DO
848 
849  ! call Ewald_k_Space(Coulomb_Pot_k,Coulomb_Force_k,RX,RY,RZ,LBox,dq_v,Nr_atoms,Coulomb_acc,TIMERATIO,Max_Nr_Neigh)
850  ! Coulomb_Pot_dq_v = Coulomb_Pot_Real+Coulomb_Pot_k
851 
852  row1 = 0.0_dp
853  do j = 1,hdim
854  call bml_set_row(h1_bml,j,row1,threshold)
855  enddo
856 
857  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(J,K)
858  do j = 1,nr_atoms
859  do k = h_index_start(j),h_index_end(j)
860  row1(k) = hubbard_u(j)*dq_v(j) + coulomb_pot_dq_v(j)
861  enddo
862  enddo
863  !$OMP END PARALLEL DO
864 
865  call bml_set_diagonal(h1_bml,row1,threshold)
866  call bml_multiply(s_bml,h1_bml,d1_bml,1.0_dp,0.0_dp,threshold)
867  call bml_multiply(h1_bml,s_bml,d1_bml,0.5_dp,0.5_dp,threshold)
868  call bml_copy(d1_bml,h1_bml)
869 
870  call bml_transpose_new(z_bml,x_bml)
871  call bml_multiply(x_bml,h1_bml,d1_bml,1.0_dp,0.0_dp,threshold)
872  call bml_multiply(d1_bml,z_bml,h1_bml,1.0_dp,0.0_dp,threshold)
873 
874  call prg_implicit_fermi_first_order_response(ho_bml,h1_bml,do_bml,d1_bml,inv_bml, &
875  m_rec, mu, beta, nocc, threshold)
876 
877  call bml_transpose_new(z_bml,x_bml)
878  call bml_multiply(z_bml,d1_bml,y_bml,2.0_dp,0.0_dp,threshold)
879  call bml_multiply(y_bml,x_bml,d1_bml,1.0_dp,0.0_dp,threshold)
880 
881  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(It,K,row1,row2)
882  do it = 1, nr_atoms
883  dq_dv(it) = 0
884  do k = h_index_start(it), h_index_end(it)
885  call bml_get_row(s_bml,k,row1)
886  call bml_get_row(d1_bml,k,row2)
887  dq_dv(it) = dq_dv(it) + dot_product(row1,row2)
888  enddo
889  enddo
890  !$OMP END PARALLEL DO
891 
892  dr = dq_dv - v
893  fi(:,i) = dr
894 
895  l = i
896  allocate(o(l,l), m(l,l))
897  do k = 1,l
898  do j = 1,l
899  o(k,j) = dot_product(fi(:,k),fi(:,j)) ! O_KJ = < fv_i(K) | fv_i(J) > see below Eq. (31)
900  enddo
901  enddo
902  call invert(o,m,l) ! M = O^(-1)
903  identres = 0.d0*res
904  kres = 0.d0
905  do k = 1,l
906  do j = 1,l
907  proj_tmp = m(k,j)*dot_product(fi(:,j),res)
908  identres = identres + proj_tmp*fi(:,k)
909  kres = kres + proj_tmp*vi(:,k) !! KRes becomes the rank-L approximate of (J)^(-1)*(q[n]-n)
910  enddo
911  enddo
912  fel = norm2(identres-res)/norm2(identres) !! RELATIVE RESIDUAL ERROR ESTIMATE Eq. (48) Ref. [*]
913  write(*,*) '# I, L, Fel = ',i,l,fel !! Fel goes down with L
914 
915  ! Does not work, need to normalize matrix
916  if ((fel > feltol).AND.(i < (lmax))) then
917  deallocate(o, m)
918  else
919  do row = 1,nr_atoms
920  do col = 1,nr_atoms
921  elem = 0.0
922  do k = 1,l
923  do j = 1,l
924  elem = elem + m(j,k)*vi(row,j)*fi(col,k)
925  enddo
926  enddo
927  if (abs(elem) > threshold) then
928  call bml_set_element(kk0_bml,row,col,elem)
929  endif
930  enddo
931  enddo
932  deallocate(o,m)
933  endif
934  enddo
935  deallocate(row1); deallocate(row2);
936 
937  end subroutine prg_kernel_matrix_multirank
938 
939 end module prg_xlbokernel_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_ewald_mod
Definition: prg_ewald_mod.F90:2
prg_ewald_mod::ewald_k_space_latte
subroutine, public ewald_k_space_latte(COULOMBV, RXYZ, Box, DELTAQ, Nr_atoms, COULACC, Max_Nr_Neigh)
Definition: prg_ewald_mod.F90:745
prg_ewald_mod::ewald_k_space_latte_single
subroutine, public ewald_k_space_latte_single(COULOMBV, J, RXYZ, Box, DELTAQ, Nr_atoms, COULACC)
Definition: prg_ewald_mod.F90:999
prg_ewald_mod::ewald_real_space_single_latte
subroutine, public ewald_real_space_single_latte(COULOMBV, I, RXYZ, Box, Nr_elem, DELTAQ, J, U, Element_Pointer, Nr_atoms, COULACC, HDIM, Max_Nr_Neigh)
Find Coulomb potential on site I from single charge at site J.
Definition: prg_ewald_mod.F90:30
prg_ewald_mod::ewald_real_space
subroutine, public ewald_real_space(COULOMBV, FCOUL, I, RX, RY, RZ, LBox, DELTAQ, U, Element_Type, Nr_atoms, COULACC, TIMERATIO, nnRx, nnRy, nnRz, nrnnlist, nnType, HDIM, Max_Nr_Neigh)
Definition: prg_ewald_mod.F90:623
prg_ewald_mod::ewald_real_space_latte
subroutine, public ewald_real_space_latte(COULOMBV, I, RXYZ, Box, DELTAQ, U, Element_Pointer, Nr_atoms, COULACC, nebcoul, totnebcoul, HDIM, Max_Nr_Neigh, Nr_Elem)
Definition: prg_ewald_mod.F90:389
prg_ewald_mod::ewald_real_space_single
subroutine, public ewald_real_space_single(COULOMBV, FCOUL, I, RX, RY, RZ, LBox, DELTAQ, J, U, Element_Type, Nr_atoms, COULACC, TIMERATIO, HDIM, Max_Nr_Neigh)
Definition: prg_ewald_mod.F90:142
prg_implicit_fermi_mod
Definition: prg_implicit_fermi_mod.F90:6
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_normalize_mod
The prg_normalize module.
Definition: prg_normalize_mod.F90:5
prg_parallel_mod
The parallel module.
Definition: prg_parallel_mod.F90:6
prg_timer_mod
The timer module.
Definition: prg_timer_mod.F90:23
prg_timer_mod::prg_timer_shutdown
subroutine, public prg_timer_shutdown()
Done with timers.
Definition: prg_timer_mod.F90:205
prg_timer_mod::timer_prg_init
subroutine, public timer_prg_init()
Initialize timers.
Definition: prg_timer_mod.F90:132
prg_timer_mod::prg_timer_start
subroutine, public prg_timer_start(itimer, tag)
Start Timing.
Definition: prg_timer_mod.F90:215
prg_timer_mod::prg_timer_stop
subroutine, public prg_timer_stop(itimer, verbose)
Stop timing.
Definition: prg_timer_mod.F90:233
prg_xlbokernel_mod
Pre-conditioned O(N) calculation of the kernel for XL-BOMD.
Definition: prg_xlbokernel_mod.F90:5
prg_xlbokernel_mod::prg_kernel_multirank
subroutine, public prg_kernel_multirank(KRes, KK0_bml, Res, FelTol, L, LMAX, HO_bml, mu, beta, RX, RY, RZ, LBox, Hubbard_U, Element_Type, Nr_atoms, HDIM, Max_Nr_Neigh, Coulomb_acc, TIMERATIO, nnRx, nnRy, nnRz, nrnnlist, nnType, H_INDEX_START, H_INDEX_END, S_bml, Z_bml, Nocc, Znuc, Inv_bml, H1_bml, X_bml, Y_bml, DO_bml, D1_bml, m_rec, threshold)
Definition: prg_xlbokernel_mod.F90:337
prg_xlbokernel_mod::prg_kernel_matrix_multirank
subroutine, public prg_kernel_matrix_multirank(KRes, KK0_bml, Res, FelTol, L, LMAX, HO_bml, mu, beta, RX, RY, RZ, LBox, Hubbard_U, Element_Type, Nr_atoms, HDIM, Max_Nr_Neigh, Coulomb_acc, TIMERATIO, nnRx, nnRy, nnRz, nrnnlist, nnType, H_INDEX_START, H_INDEX_END, S_bml, Z_bml, Nocc, Znuc, Inv_bml, H1_bml, X_bml, Y_bml, DO_bml, D1_bml, m_rec, threshold)
Definition: prg_xlbokernel_mod.F90:789
prg_xlbokernel_mod::prg_full_kernel
subroutine, public prg_full_kernel(KK, DO_bml, mu0, RX, RY, RZ, LBox, Hubbard_U, Element_Type, Nr_atoms, HDIM, Max_Nr_Neigh, Coulomb_acc, TIMERATIO, nnRx, nnRy, nnRz, nrnnlist, nnType, H_INDEX_START, H_INDEX_END, S_bml, Z_bml, Inv_bml, D1_bml, H1_bml, HO_bml, Y_bml, X_bml, Nocc, Znuc, m_rec, threshold, beta, diagonal)
Compute full inverse Jacobian of q[n]-n (for development code)
Definition: prg_xlbokernel_mod.F90:684
prg_xlbokernel_mod::prg_kernel_multirank_latte
subroutine, public prg_kernel_multirank_latte(KRes, KK0_bml, Res, FelTol, L, LMAX, NUMRANK, HO_bml, mu, beta, RXYZ, Box, Hubbard_U, Element_Pointer, Nr_atoms, HDIM, Max_Nr_Neigh, Coulomb_acc, nebcoul, totnebcoul, Hinxlist, S_bml, Z_bml, Nocc, Inv_bml, H1_bml, DO_bml, D1_bml, m_rec, threshold, Nr_elem, mdstep, update)
Compute low rank approximation of (K0*J)^(-1)*K0*(q[n]-n)(for LATTE)
Definition: prg_xlbokernel_mod.F90:87
prg_xlbokernel_mod::invert
subroutine invert(A, AI, N)
Definition: prg_xlbokernel_mod.F90:31
prg_xlbokernel_mod::prg_full_kernel_latte
subroutine, public prg_full_kernel_latte(KK, DO_bml, mu0, RXYZ, Box, Hubbard_U, Element_Pointer, Nr_atoms, HDIM, Max_Nr_Neigh, Coulomb_acc, Hinxlist, S_bml, Z_bml, Inv_bml, D1_bml, H1_bml, HO_bml, Nocc, m_rec, threshold, beta, Nr_elem, nebcoul, totnebcoul)
Compute full inverse Jacobian of q[n]-n (for LATTE)
Definition: prg_xlbokernel_mod.F90:515
prg_xlbokernel_mod::prg_update_preconditioner
subroutine prg_update_preconditioner(K0, v, fv, Nr_atoms, threshold)
Definition: prg_xlbokernel_mod.F90:266