1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 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