1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Routines using linear scaling chebyshev methods
8!> \par History
9!>       2012.10 created [Jinwoong Cha]
10!> \author Jinwoong Cha
11! **************************************************************************************************
12MODULE dm_ls_chebyshev
13   USE arnoldi_api,                     ONLY: arnoldi_extremal
14   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
15                                              cp_logger_get_default_unit_nr,&
16                                              cp_logger_type
17   USE cp_output_handling,              ONLY: cp_p_file,&
18                                              cp_print_key_finished_output,&
19                                              cp_print_key_should_output,&
20                                              cp_print_key_unit_nr
21   USE dbcsr_api,                       ONLY: &
22        dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_frobenius_norm, &
23        dbcsr_get_info, dbcsr_get_occupation, dbcsr_multiply, dbcsr_release, dbcsr_scale, &
24        dbcsr_set, dbcsr_trace, dbcsr_type, dbcsr_type_no_symmetry
25   USE dm_ls_scf_qs,                    ONLY: write_matrix_to_cube
26   USE dm_ls_scf_types,                 ONLY: ls_scf_env_type
27   USE input_section_types,             ONLY: section_get_ivals,&
28                                              section_vals_val_get
29   USE kinds,                           ONLY: default_string_length,&
30                                              dp
31   USE machine,                         ONLY: m_flush,&
32                                              m_walltime
33   USE mathconstants,                   ONLY: pi
34   USE qs_environment_types,            ONLY: qs_environment_type
35#include "./base/base_uses.f90"
36
37   IMPLICIT NONE
38
39   PRIVATE
40
41   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_chebyshev'
42
43   PUBLIC :: compute_chebyshev
44
45CONTAINS
46
47! **************************************************************************************************
48!> \brief compute chebyshev polynomials up to order n for a given value of x
49!> \param value ...
50!> \param x ...
51!> \param n ...
52!> \par History
53!>       2012.11 created [Jinwoong Cha]
54!> \author Jinwoong Cha
55! **************************************************************************************************
56   SUBROUTINE chebyshev_poly(value, x, n)
57      REAL(KIND=dp), INTENT(OUT)                         :: value
58      REAL(KIND=dp), INTENT(IN)                          :: x
59      INTEGER, INTENT(IN)                                :: n
60
61!polynomial values
62!number of chev polynomials
63
64      value = COS((n - 1)*ACOS(x))
65
66   END SUBROUTINE chebyshev_poly
67
68! **************************************************************************************************
69!> \brief kernel for chebyshev polynomials expansion (Jackson kernel)
70!> \param value ...
71!> \param n ...
72!> \param nc ...
73!> \par History
74!>       2012.11 created [Jinwoong Cha]
75!> \author Jinwoong Cha
76! **************************************************************************************************
77   SUBROUTINE kernel(value, n, nc)
78      REAL(KIND=dp), INTENT(OUT)                         :: value
79      INTEGER, INTENT(IN)                                :: n, nc
80
81!kernel at n
82!n-1 order of chebyshev polynomials
83!number of total chebyshev polynomials
84!Kernel define
85
86      value = 1.0_dp/(nc + 1.0_dp)*((nc - (n - 1) + 1.0_dp)* &
87                                    COS(pi*(n - 1)/(nc + 1.0_dp)) + SIN(pi*(n - 1)/(nc + 1.0_dp))*1.0_dp/TAN(pi/(nc + 1.0_dp)))
88
89   END SUBROUTINE kernel
90
91! **************************************************************************************************
92!> \brief compute properties based on chebyshev expansion
93!> \param qs_env ...
94!> \param ls_scf_env ...
95!> \par History
96!>       2012.10 created [Jinwoong Cha]
97!> \author Jinwoong Cha
98! **************************************************************************************************
99   SUBROUTINE compute_chebyshev(qs_env, ls_scf_env)
100      TYPE(qs_environment_type), POINTER                 :: qs_env
101      TYPE(ls_scf_env_type)                              :: ls_scf_env
102
103      CHARACTER(len=*), PARAMETER :: routineN = 'compute_chebyshev', &
104         routineP = moduleN//':'//routineN
105      REAL(KIND=dp), PARAMETER                           :: scale_evals = 1.01_dp
106
107      CHARACTER(LEN=30)                                  :: middle_name
108      CHARACTER(LEN=default_string_length)               :: title
109      INTEGER                                            :: handle, icheb, igrid, iinte, ispin, &
110                                                            iwindow, n_gridpoint_dos, ncheb, &
111                                                            ninte, Nrows, nwindow, unit_cube, &
112                                                            unit_dos, unit_nr
113      LOGICAL                                            :: converged, write_cubes
114      REAL(KIND=dp) :: chev_T, chev_T_dos, dummy1, final, frob_matrix, initial, interval_a, &
115         interval_b, max_ev, min_ev, occ, orbital_occ, summa, t1, t2
116      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: chev_E, chev_Es_dos, dos, dummy2, ev1, &
117                                                            ev2, kernel_g, mu, sev1, sev2, trace_dm
118      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: aitchev_T, E_inte, gdensity, sqrt_vec
119      REAL(KIND=dp), DIMENSION(:), POINTER               :: tmp_r
120      TYPE(cp_logger_type), POINTER                      :: logger
121      TYPE(dbcsr_type)                                   :: matrix_dummy1, matrix_F, matrix_tmp1, &
122                                                            matrix_tmp2, matrix_tmp3
123      TYPE(dbcsr_type), DIMENSION(:), POINTER            :: matrix_dummy2
124
125      IF (.NOT. ls_scf_env%chebyshev%compute_chebyshev) RETURN
126
127      CALL timeset(routineN, handle)
128
129      ! get a useful output_unit
130      logger => cp_get_default_logger()
131      IF (logger%para_env%ionode) THEN
132         unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
133      ELSE
134         unit_nr = -1
135      ENDIF
136
137      ncheb = ls_scf_env%chebyshev%n_chebyshev
138      ninte = 2*ncheb
139      n_gridpoint_dos = ls_scf_env%chebyshev%n_gridpoint_dos
140
141      write_cubes = BTEST(cp_print_key_should_output(logger%iter_info, ls_scf_env%chebyshev%print_key_cube), cp_p_file)
142      IF (write_cubes) THEN
143         IF (ASSOCIATED(ls_scf_env%chebyshev%min_energy)) DEALLOCATE (ls_scf_env%chebyshev%min_energy)
144         CALL section_vals_val_get(ls_scf_env%chebyshev%print_key_cube, "MIN_ENERGY", r_vals=tmp_r)
145         ALLOCATE (ls_scf_env%chebyshev%min_energy(SIZE(tmp_r)))
146         ls_scf_env%chebyshev%min_energy = tmp_r
147
148         IF (ASSOCIATED(ls_scf_env%chebyshev%max_energy)) DEALLOCATE (ls_scf_env%chebyshev%max_energy)
149         CALL section_vals_val_get(ls_scf_env%chebyshev%print_key_cube, "MAX_ENERGY", r_vals=tmp_r)
150         ALLOCATE (ls_scf_env%chebyshev%max_energy(SIZE(tmp_r)))
151         ls_scf_env%chebyshev%max_energy = tmp_r
152
153         nwindow = SIZE(ls_scf_env%chebyshev%min_energy)
154      ELSE
155         nwindow = 0
156      ENDIF
157
158      ALLOCATE (ev1(1:nwindow))
159      ALLOCATE (ev2(1:nwindow))
160      ALLOCATE (sev1(1:nwindow))
161      ALLOCATE (sev2(1:nwindow))
162      ALLOCATE (trace_dm(1:nwindow))
163      ALLOCATE (matrix_dummy2(1:nwindow))
164
165      DO iwindow = 1, nwindow
166         ev1(iwindow) = ls_scf_env%chebyshev%min_energy(iwindow)
167         ev2(iwindow) = ls_scf_env%chebyshev%max_energy(iwindow)
168      END DO
169
170      IF (unit_nr > 0) THEN
171         WRITE (unit_nr, '()')
172         WRITE (unit_nr, '(T2,A)') "STARTING CHEBYSHEV CALCULATION"
173      ENDIF
174
175      ! create 3 temporary matrices
176      CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
177      CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
178      CALL dbcsr_create(matrix_tmp3, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
179      CALL dbcsr_create(matrix_F, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
180      CALL dbcsr_create(matrix_dummy1, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
181
182      DO iwindow = 1, nwindow
183         CALL dbcsr_create(matrix_dummy2(iwindow), template=ls_scf_env%matrix_s, &
184                           matrix_type=dbcsr_type_no_symmetry)
185      END DO
186
187      DO ispin = 1, SIZE(ls_scf_env%matrix_ks)
188         ! create matrix_F=inv(sqrt(S))*H*inv(sqrt(S))
189         CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_ks(ispin), &
190                             0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
191         CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, &
192                             0.0_dp, matrix_F, filter_eps=ls_scf_env%eps_filter)
193
194         ! find largest and smallest eigenvalues
195         CALL arnoldi_extremal(matrix_F, max_ev, min_ev, converged=converged, max_iter=ls_scf_env%max_iter_lanczos, &
196                               threshold=ls_scf_env%eps_lanczos) !Lanczos algorithm to calculate eigenvalue
197         IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,2F16.8,A,L2)') &
198            "smallest largest eigenvalue", min_ev, max_ev, " converged ", converged
199         IF (nwindow > 0) THEN
200            IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,1000F16.8)') "requested interval-min_energy", ev1(:)
201            IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,1000F16.8)') "requested interval-max_energy", ev2(:)
202         ENDIF
203         interval_a = (max_ev - min_ev)*scale_evals/2
204         interval_b = (max_ev + min_ev)/2
205
206         sev1(:) = (ev1(:) - interval_b)/interval_a !scaled ev1 vector
207         sev2(:) = (ev2(:) - interval_b)/interval_a !scaled ev2 vector
208
209         !chebyshev domain,pi*sqrt(1-x^2) vector construction and chebyshev polynomials for integration (for g(E))
210         ALLOCATE (E_inte(1:ninte + 1, 1:nwindow))
211         ALLOCATE (sqrt_vec(1:ninte + 1, 1:nwindow))
212
213         DO iwindow = 1, nwindow
214            DO iinte = 1, ninte + 1
215               E_inte(iinte, iwindow) = sev1(iwindow) + ((sev2(iwindow) - sev1(iwindow))/ninte)*(iinte - 1)
216               sqrt_vec(iinte, iwindow) = pi*SQRT(1.0_dp - E_inte(iinte, iwindow)*E_inte(iinte, iwindow))
217            END DO
218         END DO
219
220         !integral.. (identical to the coefficient for g(E))
221
222         ALLOCATE (aitchev_T(1:ncheb, 1:nwindow)) !after intergral. =>ainte
223
224         DO iwindow = 1, nwindow
225            DO icheb = 1, ncheb
226               CALL chebyshev_poly(initial, E_inte(1, iwindow), icheb)
227               CALL chebyshev_poly(final, E_inte(1, iwindow), icheb)
228          summa = (sev2(iwindow) - sev1(iwindow))/(2.0_dp*ninte)*(initial/sqrt_vec(1, iwindow) + final/sqrt_vec(ninte + 1, iwindow))
229               DO iinte = 2, ninte
230                  CALL chebyshev_poly(chev_T, E_inte(iinte, iwindow), icheb)
231                  summa = summa + ((sev2(iwindow) - sev1(iwindow))/ninte)*(chev_T/sqrt_vec(iinte, iwindow))
232               END DO
233               aitchev_T(icheb, iwindow) = summa
234               summa = 0
235            END DO
236         END DO
237
238         ! scale the matrix to get evals in the interval -1,1
239         CALL dbcsr_add_on_diag(matrix_F, -interval_b)
240         CALL dbcsr_scale(matrix_F, 1/interval_a)
241
242         ! compute chebyshev matrix recursion
243         CALL dbcsr_get_info(matrix=matrix_F, nfullrows_total=Nrows) !get information about a matrix
244         CALL dbcsr_set(matrix_dummy1, 0.0_dp) !empty matrix creation(for density matrix)
245
246         DO iwindow = 1, nwindow
247            CALL dbcsr_set(matrix_dummy2(iwindow), 0.0_dp) !empty matrix creation(for density matrix)
248         END DO
249
250         ALLOCATE (mu(1:ncheb))
251         ALLOCATE (kernel_g(1:ncheb))
252         CALL kernel(kernel_g(1), 1, ncheb)
253         CALL kernel(kernel_g(2), 2, ncheb)
254
255         CALL dbcsr_set(matrix_tmp1, 0.0_dp) !matrix creation
256         CALL dbcsr_add_on_diag(matrix_tmp1, 1.0_dp) !add a only number to diagonal elements
257         CALL dbcsr_trace(matrix_tmp1, trace=mu(1))
258         CALL dbcsr_copy(matrix_tmp2, matrix_F) !make matrix_tmp2 = matrix_F
259         CALL dbcsr_trace(matrix_tmp2, trace=mu(2))
260
261         DO iwindow = 1, nwindow
262            CALL dbcsr_copy(matrix_dummy1, matrix_tmp1)
263            CALL dbcsr_copy(matrix_dummy2(iwindow), matrix_tmp2) !matrix_dummy2=
264            CALL dbcsr_scale(matrix_dummy1, kernel_g(1)*aitchev_T(1, iwindow)) !first term of chebyshev poly(matrix)
265            CALL dbcsr_scale(matrix_dummy2(iwindow), 2.0_dp*kernel_g(2)*aitchev_T(2, iwindow)) !second term of chebyshev poly(matrix)
266
267            CALL dbcsr_add(matrix_dummy2(iwindow), matrix_dummy1, 1.0_dp, 1.0_dp)
268         END DO
269
270         DO icheb = 2, ncheb - 1
271            t1 = m_walltime()
272            CALL dbcsr_multiply("N", "N", 2.0_dp, matrix_F, matrix_tmp2, &
273                                -1.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter) !matrix multiplication(Recursion)
274            CALL dbcsr_copy(matrix_tmp3, matrix_tmp1)
275            CALL dbcsr_copy(matrix_tmp1, matrix_tmp2)
276            CALL dbcsr_copy(matrix_tmp2, matrix_tmp3)
277            CALL dbcsr_trace(matrix_tmp2, trace=mu(icheb + 1)) !icheb+1 th coefficient
278            CALL kernel(kernel_g(icheb + 1), icheb + 1, ncheb)
279
280            DO iwindow = 1, nwindow
281
282               CALL dbcsr_copy(matrix_dummy1, matrix_tmp2)
283               CALL dbcsr_scale(matrix_dummy1, 2.0_dp*kernel_g(icheb + 1)*aitchev_T(icheb + 1, iwindow)) !second term of chebyshev poly(matrix)
284               CALL dbcsr_add(matrix_dummy2(iwindow), matrix_dummy1, 1.0_dp, 1.0_dp)
285               CALL dbcsr_trace(matrix_dummy2(iwindow), trace=trace_dm(iwindow)) !icheb+1 th coefficient
286
287            END DO
288
289            occ = dbcsr_get_occupation(matrix_tmp1)
290            t2 = m_walltime()
291            IF (unit_nr > 0 .AND. MOD(icheb, 20) == 0) THEN
292               CALL m_flush(unit_nr)
293               IF (nwindow > 0) THEN
294                  WRITE (unit_nr, '(T2,A,I5,1X,A,1X,F8.3,1X,A,1X,F8.6,1X,A,1X,1000F16.8)') &
295                     "Iter.", icheb, "time=", t2 - t1, "occ=", occ, "traces=", trace_dm(:)
296               ELSE
297                  WRITE (unit_nr, '(T2,A,I5,1X,A,1X,F8.3,1X,A,1X,F8.6)') &
298                     "Iter.", icheb, "time=", t2 - t1, "occ=", occ
299               ENDIF
300            ENDIF
301         ENDDO
302
303         DO iwindow = 1, nwindow
304            IF (SIZE(ls_scf_env%matrix_ks) == 1) THEN
305               orbital_occ = 2.0_dp
306            ELSE
307               orbital_occ = 1.0_dp
308            ENDIF
309            CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, matrix_dummy2(iwindow), &
310                                0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
311            CALL dbcsr_multiply("N", "N", orbital_occ, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, &
312                                0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
313            CALL dbcsr_copy(matrix_dummy2(iwindow), matrix_tmp2)
314
315            ! look at the difference with the density matrix from the ls routines
316            IF (.FALSE.) THEN
317               CALL dbcsr_copy(matrix_tmp1, matrix_tmp2)
318               CALL dbcsr_add(matrix_tmp1, ls_scf_env%matrix_p(ispin), 1.0_dp, -1.0_dp) !comparison
319               frob_matrix = dbcsr_frobenius_norm(matrix_tmp1)
320               IF (unit_nr > 0) WRITE (unit_nr, *) "Difference between Chebyshev DM and LS DM", frob_matrix
321            ENDIF
322         END DO
323
324         write_cubes = BTEST(cp_print_key_should_output(logger%iter_info, &
325                                                        ls_scf_env%chebyshev%print_key_cube), cp_p_file)
326         IF (write_cubes) THEN
327            DO iwindow = 1, nwindow
328               WRITE (middle_name, "(A,I0)") "E_DENSITY_WINDOW_", iwindow
329               WRITE (title, "(A,1X,F16.8,1X,A,1X,F16.8)") "Energy range : ", ev1(iwindow), "to", ev2(iwindow)
330               unit_cube = cp_print_key_unit_nr(logger, ls_scf_env%chebyshev%print_key_cube, &
331                                                "", extension=".cube", & !added 01/22/2012
332                                                middle_name=TRIM(middle_name), log_filename=.FALSE.)
333               CALL write_matrix_to_cube(qs_env, ls_scf_env, matrix_dummy2(iwindow), unit_cube, title, &
334                                         section_get_ivals(ls_scf_env%chebyshev%print_key_cube, "STRIDE"))
335               CALL cp_print_key_finished_output(unit_cube, logger, ls_scf_env%chebyshev%print_key_cube, "")
336            END DO
337         ENDIF
338
339      ENDDO
340
341      ! Chebyshev expansion with calculated coefficient
342      ! grid construction and rescaling (by J)
343      unit_dos = cp_print_key_unit_nr(logger, ls_scf_env%chebyshev%print_key_dos, "", extension=".xy", &
344                                      middle_name="DOS", log_filename=.FALSE.)
345
346      IF (unit_dos > 0) THEN
347         ALLOCATE (dos(1:n_gridpoint_dos))
348         ALLOCATE (gdensity(1:n_gridpoint_dos, 1:nwindow))
349         ALLOCATE (chev_E(1:n_gridpoint_dos))
350         ALLOCATE (chev_Es_dos(1:n_gridpoint_dos))
351         ALLOCATE (dummy2(1:nwindow))
352         DO igrid = 1, n_gridpoint_dos
353            chev_E(igrid) = min_ev + (igrid - 1)*(max_ev - min_ev)/(n_gridpoint_dos - 1)
354            chev_Es_dos(igrid) = (chev_E(igrid) - interval_b)/interval_a
355         END DO
356         DO igrid = 1, n_gridpoint_dos
357            dummy1 = 0.0_dp !summation of polynomials
358            dummy2(:) = 0.0_dp !summation of polynomials
359            DO icheb = 2, ncheb
360               CALL chebyshev_poly(chev_T_dos, chev_Es_dos(igrid), icheb)
361               dummy1 = dummy1 + kernel_g(icheb)*mu(icheb)*chev_T_dos
362               DO iwindow = 1, nwindow
363                  dummy2(iwindow) = dummy2(iwindow) + kernel_g(icheb)*aitchev_T(icheb, iwindow)*chev_T_dos
364               END DO
365            END DO
366            dos(igrid) = 1.0_dp/(interval_a*Nrows* &
367                                 (pi*SQRT(1.0_dp - chev_Es_dos(igrid)*chev_Es_dos(igrid))))*(kernel_g(1)*mu(1) + 2.0_dp*dummy1)
368            DO iwindow = 1, nwindow
369               gdensity(igrid, iwindow) = kernel_g(1)*aitchev_T(1, iwindow) + 2.0_dp*dummy2(iwindow)
370            END DO
371            WRITE (unit_dos, '(1000F16.8)') chev_E(igrid), dos(igrid), gdensity(igrid, :)
372         END DO
373         DEALLOCATE (chev_Es_dos, chev_E, dos, gdensity)
374      ENDIF
375      CALL cp_print_key_finished_output(unit_dos, logger, ls_scf_env%chebyshev%print_key_dos, "")
376
377      ! free the matrices
378      CALL dbcsr_release(matrix_tmp1)
379      CALL dbcsr_release(matrix_tmp2)
380      CALL dbcsr_release(matrix_tmp3)
381      CALL dbcsr_release(matrix_F)
382      CALL dbcsr_release(matrix_dummy1)
383
384      DO iwindow = 1, nwindow
385         CALL dbcsr_release(matrix_dummy2(iwindow))
386      END DO
387
388      DEALLOCATE (ev1, ev2, sev1, sev2, matrix_dummy2)
389
390      !Need deallocation
391      DEALLOCATE (mu, kernel_g, aitchev_T, E_inte, sqrt_vec)
392
393      IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') "ENDING CHEBYSHEV CALCULATION"
394
395      CALL timestop(handle)
396
397   END SUBROUTINE compute_chebyshev
398
399END MODULE dm_ls_chebyshev
400