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 to calculate weights for correlation methods 8!> \par History 9!> 05.2019 Refactored from rpa_ri_gpw [Frederick Stein] 10! ************************************************************************************************** 11MODULE mp2_weights 12 USE cp_fm_types, ONLY: cp_fm_get_info,& 13 cp_fm_type 14 USE cp_para_env, ONLY: cp_para_env_create,& 15 cp_para_env_release 16 USE cp_para_types, ONLY: cp_para_env_type 17 USE kinds, ONLY: dp 18 USE kpoint_types, ONLY: get_kpoint_info,& 19 kpoint_env_type,& 20 kpoint_type 21 USE machine, ONLY: m_flush 22 USE mathconstants, ONLY: pi 23 USE message_passing, ONLY: mp_bcast,& 24 mp_comm_split_direct,& 25 mp_sum 26 USE minimax_exp, ONLY: get_exp_minimax_coeff 27 USE minimax_rpa, ONLY: get_rpa_minimax_coeff 28 USE mp2_types, ONLY: mp2_type 29 USE qs_environment_types, ONLY: get_qs_env,& 30 qs_environment_type 31 USE qs_mo_types, ONLY: get_mo_set,& 32 mo_set_type 33#include "./base/base_uses.f90" 34 35 IMPLICIT NONE 36 37 PRIVATE 38 39 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_weights' 40 41 PUBLIC :: get_minimax_weights, get_clenshaw_weights, test_least_square_ft 42 43CONTAINS 44 45! ************************************************************************************************** 46!> \brief ... 47!> \param para_env ... 48!> \param unit_nr ... 49!> \param homo ... 50!> \param Eigenval ... 51!> \param num_integ_points ... 52!> \param do_im_time ... 53!> \param do_ri_sos_laplace_mp2 ... 54!> \param do_print ... 55!> \param tau_tj ... 56!> \param tau_wj ... 57!> \param qs_env ... 58!> \param do_gw_im_time ... 59!> \param do_kpoints_cubic_RPA ... 60!> \param ext_scaling ... 61!> \param a_scaling ... 62!> \param e_fermi ... 63!> \param tj ... 64!> \param wj ... 65!> \param mp2_env ... 66!> \param weights_cos_tf_t_to_w ... 67!> \param weights_cos_tf_w_to_t ... 68!> \param weights_sin_tf_t_to_w ... 69!> \param homo_beta ... 70!> \param dimen_ia_beta ... 71!> \param Eigenval_beta ... 72! ************************************************************************************************** 73 SUBROUTINE get_minimax_weights(para_env, unit_nr, homo, Eigenval, num_integ_points, & 74 do_im_time, do_ri_sos_laplace_mp2, do_print, tau_tj, tau_wj, qs_env, do_gw_im_time, & 75 do_kpoints_cubic_RPA, ext_scaling, a_scaling, e_fermi, tj, wj, mp2_env, weights_cos_tf_t_to_w, & 76 weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, & 77 homo_beta, dimen_ia_beta, Eigenval_beta) 78 79 TYPE(cp_para_env_type), POINTER :: para_env 80 INTEGER, INTENT(IN) :: unit_nr, homo 81 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 82 INTEGER, INTENT(IN) :: num_integ_points 83 LOGICAL, INTENT(IN) :: do_im_time, do_ri_sos_laplace_mp2, & 84 do_print 85 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 86 INTENT(OUT) :: tau_tj, tau_wj 87 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env 88 LOGICAL, INTENT(IN), OPTIONAL :: do_gw_im_time, do_kpoints_cubic_RPA 89 REAL(KIND=dp), INTENT(IN), OPTIONAL :: ext_scaling 90 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: a_scaling, e_fermi 91 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 92 INTENT(OUT), OPTIONAL :: tj, wj 93 TYPE(mp2_type), OPTIONAL, POINTER :: mp2_env 94 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 95 INTENT(OUT), OPTIONAL :: weights_cos_tf_t_to_w, & 96 weights_cos_tf_w_to_t, & 97 weights_sin_tf_t_to_w 98 INTEGER, INTENT(IN), OPTIONAL :: homo_beta, dimen_ia_beta 99 REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: Eigenval_beta 100 101 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_minimax_weights', & 102 routineP = moduleN//':'//routineN 103 104 INTEGER :: handle, ierr, jquad, & 105 num_points_per_magnitude 106 LOGICAL :: my_do_kpoints, my_open_shell 107 REAL(KIND=dp) :: E_Range, Emax, Emax_beta, Emin, & 108 Emin_beta, max_error_min, scaling 109 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x_tw 110 111 CALL timeset(routineN, handle) 112 113 ! Test for spin unrestricted 114 my_open_shell = .FALSE. 115 IF (PRESENT(homo_beta) .AND. PRESENT(dimen_ia_beta) .AND. PRESENT(Eigenval_beta)) THEN 116 my_open_shell = .TRUE. 117 END IF 118 119 ! Test whether all necessary variables are available 120 my_do_kpoints = .FALSE. 121 IF (.NOT. do_ri_sos_laplace_mp2) THEN 122 IF (.NOT. (PRESENT(do_gw_im_time) .AND. PRESENT(do_kpoints_cubic_RPA) .AND. PRESENT(ext_scaling) .AND. PRESENT(a_scaling) & 123 .AND. PRESENT(e_fermi) .AND. PRESENT(tj) .AND. PRESENT(wj) .AND. PRESENT(weights_cos_tf_t_to_w) .AND. & 124 PRESENT(weights_cos_tf_w_to_t) .AND. PRESENT(weights_sin_tf_t_to_w) .AND. PRESENT(mp2_env))) THEN 125 CPABORT("Need more parameters without SOS-MP2") 126 END IF 127 my_do_kpoints = do_kpoints_cubic_RPA 128 END IF 129 130 IF (my_do_kpoints) THEN 131 132 CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, Emin, Emax, e_fermi) 133 134 E_Range = Emax/Emin 135 136 ELSE 137 138 Emin = Eigenval(homo + 1) - Eigenval(homo) 139 Emax = MAXVAL(Eigenval) - MINVAL(Eigenval) 140 IF (my_open_shell) THEN 141 IF (homo_beta > 0) THEN 142 Emin_beta = Eigenval_beta(homo_beta + 1) - Eigenval_beta(homo_beta) 143 Emax_beta = MAXVAL(Eigenval_beta) - MINVAL(Eigenval_beta) 144 Emin = MIN(Emin, Emin_beta) 145 Emax = MAX(Emax, Emax_beta) 146 END IF 147 END IF 148 E_Range = Emax/Emin 149 150 END IF 151 152 IF (.NOT. do_ri_sos_laplace_mp2) THEN 153 ALLOCATE (x_tw(2*num_integ_points)) 154 x_tw = 0.0_dp 155 ierr = 0 156 CALL get_rpa_minimax_coeff(num_integ_points, E_Range, x_tw, ierr) 157 158 ALLOCATE (tj(num_integ_points)) 159 tj = 0.0_dp 160 161 ALLOCATE (wj(num_integ_points)) 162 wj = 0.0_dp 163 164 DO jquad = 1, num_integ_points 165 tj(jquad) = x_tw(jquad) 166 wj(jquad) = x_tw(jquad + num_integ_points) 167 END DO 168 169 DEALLOCATE (x_tw) 170 171 IF (unit_nr > 0 .AND. do_print) THEN 172 WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & 173 "INTEG_INFO| Range for the minimax approximation:", E_Range 174 WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") "INTEG_INFO| Minimax parameters:", "Weights", "Abscissas" 175 DO jquad = 1, num_integ_points 176 WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") wj(jquad), tj(jquad) 177 END DO 178 CALL m_flush(unit_nr) 179 END IF 180 181 ! scale the minimax parameters 182 tj(:) = tj(:)*Emin 183 wj(:) = wj(:)*Emin 184 ELSE 185 ! When we perform SOS-MP2, we need an additional factor of 2 for the energies (compare with mp2_laplace.F) 186 ! We do not need weights etc. for the cosine transform 187 ! We do not scale Emax because it is not needed for SOS-MP2 188 Emin = Emin*2.0_dp 189 END IF 190 191 ! set up the minimax time grid 192 IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN 193 194 ALLOCATE (x_tw(2*num_integ_points)) 195 x_tw = 0.0_dp 196 197 CALL get_exp_minimax_coeff(num_integ_points, E_Range, x_tw) 198 199 ! For RPA we include already a factor of two (see later steps) 200 scaling = 2.0_dp 201 IF (do_ri_sos_laplace_mp2) scaling = 1.0_dp 202 203 ALLOCATE (tau_tj(0:num_integ_points)) 204 tau_tj = 0.0_dp 205 206 ALLOCATE (tau_wj(num_integ_points)) 207 tau_wj = 0.0_dp 208 209 DO jquad = 1, num_integ_points 210 tau_tj(jquad) = x_tw(jquad)/scaling 211 tau_wj(jquad) = x_tw(jquad + num_integ_points)/scaling 212 END DO 213 214 DEALLOCATE (x_tw) 215 216 IF (unit_nr > 0 .AND. do_print) THEN 217 WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & 218 "INTEG_INFO| Range for the minimax approximation:", E_Range 219 ! For testing the gap 220 WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & 221 "INTEG_INFO| Gap:", Emin 222 WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") & 223 "INTEG_INFO| Minimax parameters of the time grid:", "Weights", "Abscissas" 224 DO jquad = 1, num_integ_points 225 WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad) 226 END DO 227 CALL m_flush(unit_nr) 228 END IF 229 230 ! scale grid from [1,R] to [Emin,Emax] 231 tau_tj(:) = tau_tj(:)/Emin 232 tau_wj(:) = tau_wj(:)/Emin 233 234 IF (.NOT. do_ri_sos_laplace_mp2) THEN 235 ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points)) 236 weights_cos_tf_t_to_w = 0.0_dp 237 238 num_points_per_magnitude = mp2_env%ri_rpa_im_time%num_points_per_magnitude 239 CALL get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, tj, & 240 Emin, Emax, max_error_min, num_points_per_magnitude) 241 242 IF (do_gw_im_time) THEN 243 244 ! get the weights for the cosine transform W^c(iw) -> W^c(it) 245 ALLOCATE (weights_cos_tf_w_to_t(num_integ_points, num_integ_points)) 246 weights_cos_tf_w_to_t = 0.0_dp 247 248 CALL get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, tj, & 249 Emin, Emax, max_error_min, num_points_per_magnitude) 250 251 ! get the weights for the sine transform Sigma^sin(it) -> Sigma^sin(iw) (PRB 94, 165109 (2016), Eq. 71) 252 ALLOCATE (weights_sin_tf_t_to_w(num_integ_points, num_integ_points)) 253 weights_sin_tf_t_to_w = 0.0_dp 254 255 CALL get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, tj, & 256 Emin, Emax, max_error_min, num_points_per_magnitude) 257 258 IF (unit_nr > 0) THEN 259 WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") & 260 "INTEG_INFO| Maximum deviation of the imag. time fit:", max_error_min 261 END IF 262 END IF 263 264 END IF 265 266 END IF 267 268 CALL timestop(handle) 269 270 END SUBROUTINE get_minimax_weights 271 272! ************************************************************************************************** 273!> \brief ... 274!> \param para_env ... 275!> \param para_env_RPA ... 276!> \param unit_nr ... 277!> \param homo ... 278!> \param virtual ... 279!> \param Eigenval ... 280!> \param num_integ_points ... 281!> \param num_integ_group ... 282!> \param color_rpa_group ... 283!> \param fm_mat_S ... 284!> \param my_do_gw ... 285!> \param ext_scaling ... 286!> \param a_scaling ... 287!> \param tj ... 288!> \param wj ... 289!> \param homo_beta ... 290!> \param virtual_beta ... 291!> \param dimen_ia_beta ... 292!> \param Eigenval_beta ... 293!> \param fm_mat_S_beta ... 294! ************************************************************************************************** 295 SUBROUTINE get_clenshaw_weights(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, & 296 num_integ_group, color_rpa_group, fm_mat_S, my_do_gw, & 297 ext_scaling, a_scaling, tj, wj, & 298 homo_beta, virtual_beta, dimen_ia_beta, Eigenval_beta, fm_mat_S_beta) 299 300 TYPE(cp_para_env_type), POINTER :: para_env, para_env_RPA 301 INTEGER, INTENT(IN) :: unit_nr, homo, virtual 302 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 303 INTEGER, INTENT(IN) :: num_integ_points, num_integ_group, & 304 color_rpa_group 305 TYPE(cp_fm_type), POINTER :: fm_mat_S 306 LOGICAL, INTENT(IN) :: my_do_gw 307 REAL(KIND=dp), INTENT(IN) :: ext_scaling 308 REAL(KIND=dp), INTENT(OUT) :: a_scaling 309 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 310 INTENT(OUT) :: tj, wj 311 INTEGER, INTENT(IN), OPTIONAL :: homo_beta, virtual_beta, dimen_ia_beta 312 REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: Eigenval_beta 313 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_mat_S_beta 314 315 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_clenshaw_weights', & 316 routineP = moduleN//':'//routineN 317 318 INTEGER :: handle, jquad 319 LOGICAL :: my_open_shell 320 321 CALL timeset(routineN, handle) 322 323 my_open_shell = .FALSE. 324 IF (PRESENT(homo_beta) .AND. PRESENT(virtual_beta) .AND. PRESENT(dimen_ia_beta) .AND. PRESENT(Eigenval_beta) .AND. & 325 PRESENT(fm_mat_S_beta)) THEN 326 my_open_shell = .TRUE. 327 END IF 328 329 ! Now, start to prepare the different grid 330 ALLOCATE (tj(num_integ_points)) 331 tj = 0.0_dp 332 333 ALLOCATE (wj(num_integ_points)) 334 wj = 0.0_dp 335 336 DO jquad = 1, num_integ_points - 1 337 tj(jquad) = jquad*pi/(2.0_dp*num_integ_points) 338 wj(jquad) = pi/(num_integ_points*SIN(tj(jquad))**2) 339 END DO 340 tj(num_integ_points) = pi/2.0_dp 341 wj(num_integ_points) = pi/(2.0_dp*num_integ_points*SIN(tj(num_integ_points))**2) 342 343 a_scaling = 1.0_dp 344 IF (my_open_shell) THEN 345 CALL calc_scaling_factor(a_scaling, para_env, para_env_RPA, homo, virtual, Eigenval, & 346 num_integ_points, num_integ_group, color_rpa_group, & 347 tj, wj, fm_mat_S, & 348 homo_beta, virtual_beta, dimen_ia_beta, Eigenval_beta, fm_mat_S_beta) 349 ELSE 350 CALL calc_scaling_factor(a_scaling, para_env, para_env_RPA, homo, virtual, Eigenval, & 351 num_integ_points, num_integ_group, color_rpa_group, & 352 tj, wj, fm_mat_S) 353 END IF 354 355 ! for G0W0, we may set the scaling factor by hand 356 IF (my_do_gw .AND. ext_scaling > 0.0_dp) THEN 357 a_scaling = ext_scaling 358 END IF 359 360 IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.5)') 'INTEG_INFO| Scaling parameter:', a_scaling 361 362 wj(:) = wj(:)*a_scaling 363 364 CALL timestop(handle) 365 366 END SUBROUTINE get_clenshaw_weights 367 368! ************************************************************************************************** 369!> \brief ... 370!> \param a_scaling_ext ... 371!> \param para_env ... 372!> \param para_env_RPA ... 373!> \param homo ... 374!> \param virtual ... 375!> \param Eigenval ... 376!> \param num_integ_points ... 377!> \param num_integ_group ... 378!> \param color_rpa_group ... 379!> \param tj_ext ... 380!> \param wj_ext ... 381!> \param fm_mat_S ... 382!> \param homo_beta ... 383!> \param virtual_beta ... 384!> \param dimen_ia_beta ... 385!> \param Eigenval_beta ... 386!> \param fm_mat_S_beta ... 387! ************************************************************************************************** 388 SUBROUTINE calc_scaling_factor(a_scaling_ext, para_env, para_env_RPA, homo, virtual, Eigenval, & 389 num_integ_points, num_integ_group, color_rpa_group, & 390 tj_ext, wj_ext, fm_mat_S, & 391 homo_beta, virtual_beta, dimen_ia_beta, Eigenval_beta, fm_mat_S_beta) 392 REAL(KIND=dp), INTENT(INOUT) :: a_scaling_ext 393 TYPE(cp_para_env_type), POINTER :: para_env, para_env_RPA 394 INTEGER, INTENT(IN) :: homo, virtual 395 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 396 INTEGER, INTENT(IN) :: num_integ_points, num_integ_group, & 397 color_rpa_group 398 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 399 INTENT(IN) :: tj_ext, wj_ext 400 TYPE(cp_fm_type), POINTER :: fm_mat_S 401 INTEGER, INTENT(IN), OPTIONAL :: homo_beta, virtual_beta, dimen_ia_beta 402 REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: Eigenval_beta 403 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_mat_S_beta 404 405 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_scaling_factor', & 406 routineP = moduleN//':'//routineN 407 408 INTEGER :: handle, icycle, jquad, nrow_local, & 409 nrow_local_beta 410 LOGICAL :: my_open_shell 411 REAL(KIND=dp) :: a_high, a_low, a_scaling, conv_param, eps, first_deriv, left_term, & 412 right_term, right_term_ref, right_term_ref_beta, step 413 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cottj, D_ia, D_ia_beta, iaia_RI, & 414 iaia_RI_beta, M_ia, M_ia_beta 415 TYPE(cp_para_env_type), POINTER :: para_env_row, para_env_row_beta 416 417 CALL timeset(routineN, handle) 418 419 my_open_shell = .FALSE. 420 IF (PRESENT(homo_beta) .AND. & 421 PRESENT(virtual_beta) .AND. & 422 PRESENT(dimen_ia_beta) .AND. & 423 PRESENT(Eigenval_beta) .AND. & 424 PRESENT(fm_mat_S_beta)) my_open_shell = .TRUE. 425 426 eps = 1.0E-10_dp 427 428 ALLOCATE (cottj(num_integ_points)) 429 430 ! calculate the cotangent of the abscissa tj 431 DO jquad = 1, num_integ_points 432 cottj(jquad) = 1.0_dp/TAN(tj_ext(jquad)) 433 END DO 434 435 CALL calc_ia_ia_integrals(para_env_RPA, homo, virtual, nrow_local, right_term_ref, Eigenval, D_ia, iaia_RI, M_ia, fm_mat_S, & 436 para_env_row) 437 438 ! In the open shell case do point 1-2-3 for the beta spin 439 IF (my_open_shell) THEN 440 CALL calc_ia_ia_integrals(para_env_RPA, homo_beta, virtual_beta, nrow_local_beta, right_term_ref_beta, Eigenval_beta, & 441 D_ia_beta, iaia_RI_beta, M_ia_beta, fm_mat_S_beta, para_env_row_beta) 442 443 right_term_ref = right_term_ref + right_term_ref_beta 444 END IF 445 446 ! bcast the result 447 IF (para_env%mepos == 0) THEN 448 CALL mp_bcast(right_term_ref, 0, para_env%group) 449 ELSE 450 right_term_ref = 0.0_dp 451 CALL mp_bcast(right_term_ref, 0, para_env%group) 452 END IF 453 454 ! 5) start iteration for solving the non-linear equation by bisection 455 ! find limit, here step=0.5 seems a good compromise 456 conv_param = 100.0_dp*EPSILON(right_term_ref) 457 step = 0.5_dp 458 a_low = 0.0_dp 459 a_high = step 460 right_term = -right_term_ref 461 DO icycle = 1, num_integ_points*2 462 a_scaling = a_high 463 464 CALL calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, & 465 M_ia, cottj, wj_ext, D_ia, D_ia_beta, M_ia_beta, & 466 nrow_local, nrow_local_beta, num_integ_group, color_rpa_group, & 467 para_env, para_env_row, para_env_row_beta) 468 left_term = left_term/4.0_dp/pi*a_scaling 469 470 IF (ABS(left_term) > ABS(right_term) .OR. ABS(left_term + right_term) <= conv_param) EXIT 471 a_low = a_high 472 a_high = a_high + step 473 474 END DO 475 476 IF (ABS(left_term + right_term) >= conv_param) THEN 477 IF (a_scaling >= 2*num_integ_points*step) THEN 478 a_scaling = 1.0_dp 479 ELSE 480 481 DO icycle = 1, num_integ_points*2 482 a_scaling = (a_low + a_high)/2.0_dp 483 484 CALL calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, & 485 M_ia, cottj, wj_ext, D_ia, D_ia_beta, M_ia_beta, & 486 nrow_local, nrow_local_beta, num_integ_group, color_rpa_group, & 487 para_env, para_env_row, para_env_row_beta) 488 left_term = left_term/4.0_dp/pi*a_scaling 489 490 IF (ABS(left_term) > ABS(right_term)) THEN 491 a_high = a_scaling 492 ELSE 493 a_low = a_scaling 494 END IF 495 496 IF (ABS(a_high - a_low) < 1.0e-5_dp) EXIT 497 498 END DO 499 500 END IF 501 END IF 502 503 a_scaling_ext = a_scaling 504 CALL mp_bcast(a_scaling_ext, 0, para_env%group) 505 506 DEALLOCATE (cottj) 507 DEALLOCATE (iaia_RI) 508 DEALLOCATE (D_ia) 509 DEALLOCATE (M_ia) 510 CALL cp_para_env_release(para_env_row) 511 512 IF (my_open_shell) THEN 513 DEALLOCATE (iaia_RI_beta) 514 DEALLOCATE (D_ia_beta) 515 DEALLOCATE (M_ia_beta) 516 CALL cp_para_env_release(para_env_row_beta) 517 END IF 518 519 CALL timestop(handle) 520 521 END SUBROUTINE calc_scaling_factor 522 523! ************************************************************************************************** 524!> \brief ... 525!> \param para_env_RPA ... 526!> \param homo ... 527!> \param virtual ... 528!> \param nrow_local ... 529!> \param right_term_ref ... 530!> \param Eigenval ... 531!> \param D_ia ... 532!> \param iaia_RI ... 533!> \param M_ia ... 534!> \param fm_mat_S ... 535!> \param para_env_row ... 536! ************************************************************************************************** 537 SUBROUTINE calc_ia_ia_integrals(para_env_RPA, homo, virtual, nrow_local, right_term_ref, Eigenval, & 538 D_ia, iaia_RI, M_ia, fm_mat_S, para_env_row) 539 540 TYPE(cp_para_env_type), POINTER :: para_env_RPA 541 INTEGER, INTENT(IN) :: homo, virtual 542 INTEGER, INTENT(OUT) :: nrow_local 543 REAL(KIND=dp), INTENT(OUT) :: right_term_ref 544 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 545 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 546 INTENT(OUT) :: D_ia, iaia_RI, M_ia 547 TYPE(cp_fm_type), POINTER :: fm_mat_S 548 TYPE(cp_para_env_type), POINTER :: para_env_row 549 550 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_ia_ia_integrals', & 551 routineP = moduleN//':'//routineN 552 553 INTEGER :: avirt, color_col, color_row, comm_col, & 554 comm_row, handle, i_global, iiB, iocc, & 555 jjB, ncol_local 556 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 557 REAL(KIND=dp) :: eigen_diff 558 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: iaia_RI_dp 559 TYPE(cp_para_env_type), POINTER :: para_env_col 560 561 CALL timeset(routineN, handle) 562 563 ! calculate the (ia|ia) RI integrals 564 ! ---------------------------------- 565 ! 1) get info fm_mat_S 566 !XXX CALL cp_fm_to_fm(source=fm_mat_S,destination=fm_mat_G) 567 CALL cp_fm_get_info(matrix=fm_mat_S, & 568 nrow_local=nrow_local, & 569 ncol_local=ncol_local, & 570 row_indices=row_indices, & 571 col_indices=col_indices) 572 573 ! allocate the local buffer of iaia_RI integrals (dp kind) 574 ALLOCATE (iaia_RI_dp(nrow_local)) 575 iaia_RI_dp = 0.0_dp 576 577 ! 2) perform the local multiplication SUM_K (ia|K)*(ia|K) 578 DO jjB = 1, ncol_local 579 DO iiB = 1, nrow_local 580 iaia_RI_dp(iiB) = iaia_RI_dp(iiB) + fm_mat_S%local_data(iiB, jjB)*fm_mat_S%local_data(iiB, jjB) 581 END DO 582 END DO 583 584 ! 3) sum the result with the processes of the RPA_group having the same rows 585 ! _______K_______ _ 586 ! | | | | | | | 587 ! --> | 1 | 5 | 9 | 13| SUM --> | | 588 ! |___|__ |___|___| |_| 589 ! | | | | | | | 590 ! --> | 2 | 6 | 10| 14| SUM --> | | 591 ! ia |___|___|___|___| |_| (ia|ia)_RI 592 ! | | | | | | | 593 ! --> | 3 | 7 | 11| 15| SUM --> | | 594 ! |___|___|___|___| |_| 595 ! | | | | | | | 596 ! --> | 4 | 8 | 12| 16| SUM --> | | 597 ! |___|___|___|___| |_| 598 ! 599 600 color_row = fm_mat_S%matrix_struct%context%mepos(1) 601 CALL mp_comm_split_direct(para_env_RPA%group, comm_row, color_row) 602 NULLIFY (para_env_row) 603 CALL cp_para_env_create(para_env_row, comm_row) 604 605 CALL mp_sum(iaia_RI_dp, para_env_row%group) 606 607 ! convert the iaia_RI_dp into double-double precision 608 ALLOCATE (iaia_RI(nrow_local)) 609 DO iiB = 1, nrow_local 610 iaia_RI(iiB) = iaia_RI_dp(iiB) 611 END DO 612 DEALLOCATE (iaia_RI_dp) 613 614 ! 4) calculate the right hand term, D_ia is the matrix containing the 615 ! orbital energy differences, M_ia is the diagonal of the full RPA 'excitation' 616 ! matrix 617 ALLOCATE (D_ia(nrow_local)) 618 619 ALLOCATE (M_ia(nrow_local)) 620 621 DO iiB = 1, nrow_local 622 i_global = row_indices(iiB) 623 624 iocc = MAX(1, i_global - 1)/virtual + 1 625 avirt = i_global - (iocc - 1)*virtual 626 eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc) 627 628 D_ia(iiB) = eigen_diff 629 END DO 630 631 DO iiB = 1, nrow_local 632 M_ia(iiB) = D_ia(iiB)*D_ia(iiB) + 2.0_dp*D_ia(iiB)*iaia_RI(iiB) 633 END DO 634 635 right_term_ref = 0.0_dp 636 DO iiB = 1, nrow_local 637 right_term_ref = right_term_ref + (SQRT(M_ia(iiB)) - D_ia(iiB) - iaia_RI(iiB)) 638 END DO 639 right_term_ref = right_term_ref/2.0_dp 640 641 ! sum the result with the processes of the RPA_group having the same col 642 color_col = fm_mat_S%matrix_struct%context%mepos(2) 643 CALL mp_comm_split_direct(para_env_RPA%group, comm_col, color_col) 644 NULLIFY (para_env_col) 645 CALL cp_para_env_create(para_env_col, comm_col) 646 647 ! allocate communication array for columns 648 CALL mp_sum(right_term_ref, para_env_col%group) 649 650 CALL cp_para_env_release(para_env_col) 651 652 CALL timestop(handle) 653 654 END SUBROUTINE calc_ia_ia_integrals 655 656! ************************************************************************************************** 657!> \brief ... 658!> \param a_scaling ... 659!> \param left_term ... 660!> \param first_deriv ... 661!> \param num_integ_points ... 662!> \param my_open_shell ... 663!> \param M_ia ... 664!> \param cottj ... 665!> \param wj ... 666!> \param D_ia ... 667!> \param D_ia_beta ... 668!> \param M_ia_beta ... 669!> \param nrow_local ... 670!> \param nrow_local_beta ... 671!> \param num_integ_group ... 672!> \param color_rpa_group ... 673!> \param para_env ... 674!> \param para_env_row ... 675!> \param para_env_row_beta ... 676! ************************************************************************************************** 677 SUBROUTINE calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, & 678 M_ia, cottj, wj, D_ia, D_ia_beta, M_ia_beta, & 679 nrow_local, nrow_local_beta, num_integ_group, color_rpa_group, & 680 para_env, para_env_row, para_env_row_beta) 681 REAL(KIND=dp), INTENT(IN) :: a_scaling 682 REAL(KIND=dp), INTENT(INOUT) :: left_term, first_deriv 683 INTEGER, INTENT(IN) :: num_integ_points 684 LOGICAL, INTENT(IN) :: my_open_shell 685 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 686 INTENT(IN) :: M_ia, cottj, wj, D_ia, D_ia_beta, & 687 M_ia_beta 688 INTEGER, INTENT(IN) :: nrow_local, nrow_local_beta, & 689 num_integ_group, color_rpa_group 690 TYPE(cp_para_env_type), POINTER :: para_env, para_env_row, para_env_row_beta 691 692 INTEGER :: iiB, jquad 693 REAL(KIND=dp) :: first_deriv_beta, left_term_beta, omega 694 695 left_term = 0.0_dp 696 first_deriv = 0.0_dp 697 left_term_beta = 0.0_dp 698 first_deriv_beta = 0.0_dp 699 DO jquad = 1, num_integ_points 700 ! parallelize over integration points 701 IF (MODULO(jquad, num_integ_group) /= color_rpa_group) CYCLE 702 omega = a_scaling*cottj(jquad) 703 704 DO iiB = 1, nrow_local 705 ! parallelize over ia elements in the para_env_row group 706 IF (MODULO(iiB, para_env_row%num_pe) /= para_env_row%mepos) CYCLE 707 ! calculate left_term 708 left_term = left_term + wj(jquad)* & 709 (LOG(1.0_dp + (M_ia(iiB) - D_ia(iiB)**2)/(omega**2 + D_ia(iiB)**2)) - & 710 (M_ia(iiB) - D_ia(iiB)**2)/(omega**2 + D_ia(iiB)**2)) 711 first_deriv = first_deriv + wj(jquad)*cottj(jquad)**2* & 712 ((-M_ia(iiB) + D_ia(iiB)**2)**2/((omega**2 + D_ia(iiB)**2)**2*(omega**2 + M_ia(iiB)))) 713 END DO 714 715 IF (my_open_shell) THEN 716 DO iiB = 1, nrow_local_beta 717 ! parallelize over ia elements in the para_env_row group 718 IF (MODULO(iiB, para_env_row_beta%num_pe) /= para_env_row_beta%mepos) CYCLE 719 ! calculate left_term 720 left_term_beta = left_term_beta + wj(jquad)* & 721 (LOG(1.0_dp + (M_ia_beta(iiB) - D_ia_beta(iiB)**2)/(omega**2 + D_ia_beta(iiB)**2)) - & 722 (M_ia_beta(iiB) - D_ia_beta(iiB)**2)/(omega**2 + D_ia_beta(iiB)**2)) 723 first_deriv_beta = & 724 first_deriv_beta + wj(jquad)*cottj(jquad)**2* & 725 ((-M_ia_beta(iiB) + D_ia_beta(iiB)**2)**2/((omega**2 + D_ia_beta(iiB)**2)**2*(omega**2 + M_ia_beta(iiB)))) 726 END DO 727 END IF 728 729 END DO 730 731 ! sum the contribution from all proc, starting form the row group 732 CALL mp_sum(left_term, para_env%group) 733 CALL mp_sum(first_deriv, para_env%group) 734 735 IF (my_open_shell) THEN 736 CALL mp_sum(left_term_beta, para_env%group) 737 CALL mp_sum(first_deriv_beta, para_env%group) 738 739 left_term = left_term + left_term_beta 740 first_deriv = first_deriv + first_deriv_beta 741 END IF 742 743 END SUBROUTINE calculate_objfunc 744 745! ************************************************************************************************** 746!> \brief Calculate integration weights for the tau grid (in dependency of the omega node) 747!> \param num_integ_points ... 748!> \param tau_tj ... 749!> \param weights_cos_tf_t_to_w ... 750!> \param omega_tj ... 751!> \param E_min ... 752!> \param E_max ... 753!> \param max_error ... 754!> \param num_points_per_magnitude ... 755! ************************************************************************************************** 756 SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, omega_tj, & 757 E_min, E_max, max_error, num_points_per_magnitude) 758 759 INTEGER, INTENT(IN) :: num_integ_points 760 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 761 INTENT(IN) :: tau_tj 762 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 763 INTENT(INOUT) :: weights_cos_tf_t_to_w 764 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 765 INTENT(IN) :: omega_tj 766 REAL(KIND=dp), INTENT(IN) :: E_min, E_max 767 REAL(KIND=dp), INTENT(INOUT) :: max_error 768 INTEGER, INTENT(IN) :: num_points_per_magnitude 769 770 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_l_sq_wghts_cos_tf_t_to_w', & 771 routineP = moduleN//':'//routineN 772 773 INTEGER :: handle, iii, info, jjj, jquad, lwork, & 774 num_x_nodes 775 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork 776 REAL(KIND=dp) :: chi2_min_jquad, multiplicator, omega 777 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, tau_wj_work, vec_UTy, work, & 778 work_array, x_values, y_values 779 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, & 780 mat_SinvVSinvT, mat_U 781 782 CALL timeset(routineN, handle) 783 784 ! take num_points_per_magnitude points per 10-interval 785 num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude 786 787 ! take at least as many x points as integration points to have clear 788 ! input for the singular value decomposition 789 num_x_nodes = MAX(num_x_nodes, num_integ_points) 790 791 ALLOCATE (x_values(num_x_nodes)) 792 x_values = 0.0_dp 793 ALLOCATE (y_values(num_x_nodes)) 794 y_values = 0.0_dp 795 ALLOCATE (mat_A(num_x_nodes, num_integ_points)) 796 mat_A = 0.0_dp 797 ALLOCATE (tau_wj_work(num_integ_points)) 798 tau_wj_work = 0.0_dp 799 ALLOCATE (work_array(2*num_integ_points)) 800 work_array = 0.0_dp 801 ALLOCATE (sing_values(num_integ_points)) 802 sing_values = 0.0_dp 803 ALLOCATE (mat_U(num_x_nodes, num_x_nodes)) 804 mat_U = 0.0_dp 805 ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points)) 806 807 mat_SinvVSinvT = 0.0_dp 808 ! double the value nessary for 'A' to achieve good performance 809 lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes 810 ALLOCATE (work(lwork)) 811 work = 0.0_dp 812 ALLOCATE (iwork(8*num_integ_points)) 813 iwork = 0 814 ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes)) 815 mat_SinvVSinvSigma = 0.0_dp 816 ALLOCATE (vec_UTy(num_x_nodes)) 817 vec_UTy = 0.0_dp 818 819 max_error = 0.0_dp 820 821 ! loop over all omega frequency points 822 DO jquad = 1, num_integ_points 823 824 chi2_min_jquad = 100.0_dp 825 826 ! set the x-values logarithmically in the interval [Emin,Emax] 827 multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp)) 828 DO iii = 1, num_x_nodes 829 x_values(iii) = E_min*multiplicator**(iii - 1) 830 END DO 831 832 omega = omega_tj(jquad) 833 834 ! y=2x/(x^2+omega_k^2) 835 DO iii = 1, num_x_nodes 836 y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2 + omega**2) 837 END DO 838 839 ! calculate mat_A 840 DO jjj = 1, num_integ_points 841 DO iii = 1, num_x_nodes 842 mat_A(iii, jjj) = COS(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj)) 843 END DO 844 END DO 845 846 ! Singular value decomposition of mat_A 847 CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, & 848 mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info) 849 850 CPASSERT(info == 0) 851 852 ! integration weights = V Sigma U^T y 853 ! 1) V*Sigma 854 DO jjj = 1, num_integ_points 855 DO iii = 1, num_integ_points 856 mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)/sing_values(jjj) 857 END DO 858 END DO 859 860 ! 2) U^T y 861 CALL DGEMM('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, y_values, num_x_nodes, & 862 0.0_dp, vec_UTy, num_x_nodes) 863 864 ! 3) (V*Sigma) * (U^T y) 865 CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, vec_UTy, & 866 num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points) 867 868 weights_cos_tf_t_to_w(jquad, :) = tau_wj_work(:) 869 870 CALL calc_max_error_fit_tau_grid_with_cosine(max_error, omega, tau_tj, tau_wj_work, x_values, & 871 y_values, num_integ_points, num_x_nodes) 872 873 END DO ! jquad 874 875 DEALLOCATE (x_values, y_values, mat_A, tau_wj_work, work_array, sing_values, mat_U, mat_SinvVSinvT, & 876 work, iwork, mat_SinvVSinvSigma, vec_UTy) 877 878 CALL timestop(handle) 879 880 END SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w 881 882! ************************************************************************************************** 883!> \brief Calculate integration weights for the tau grid (in dependency of the omega node) 884!> \param num_integ_points ... 885!> \param tau_tj ... 886!> \param weights_sin_tf_t_to_w ... 887!> \param omega_tj ... 888!> \param E_min ... 889!> \param E_max ... 890!> \param max_error ... 891!> \param num_points_per_magnitude ... 892! ************************************************************************************************** 893 SUBROUTINE get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, omega_tj, & 894 E_min, E_max, max_error, num_points_per_magnitude) 895 896 INTEGER, INTENT(IN) :: num_integ_points 897 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 898 INTENT(IN) :: tau_tj 899 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 900 INTENT(INOUT) :: weights_sin_tf_t_to_w 901 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 902 INTENT(IN) :: omega_tj 903 REAL(KIND=dp), INTENT(IN) :: E_min, E_max 904 REAL(KIND=dp), INTENT(OUT) :: max_error 905 INTEGER, INTENT(IN) :: num_points_per_magnitude 906 907 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_l_sq_wghts_sin_tf_t_to_w', & 908 routineP = moduleN//':'//routineN 909 910 INTEGER :: handle, iii, info, jjj, jquad, lwork, & 911 num_x_nodes 912 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork 913 REAL(KIND=dp) :: chi2_min_jquad, multiplicator, omega 914 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, tau_wj_work, vec_UTy, work, & 915 work_array, x_values, y_values 916 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, & 917 mat_SinvVSinvT, mat_U 918 919 CALL timeset(routineN, handle) 920 921 ! take num_points_per_magnitude points per 10-interval 922 num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude 923 924 ! take at least as many x points as integration points to have clear 925 ! input for the singular value decomposition 926 num_x_nodes = MAX(num_x_nodes, num_integ_points) 927 928 ALLOCATE (x_values(num_x_nodes)) 929 x_values = 0.0_dp 930 ALLOCATE (y_values(num_x_nodes)) 931 y_values = 0.0_dp 932 ALLOCATE (mat_A(num_x_nodes, num_integ_points)) 933 mat_A = 0.0_dp 934 ALLOCATE (tau_wj_work(num_integ_points)) 935 tau_wj_work = 0.0_dp 936 ALLOCATE (work_array(2*num_integ_points)) 937 work_array = 0.0_dp 938 ALLOCATE (sing_values(num_integ_points)) 939 sing_values = 0.0_dp 940 ALLOCATE (mat_U(num_x_nodes, num_x_nodes)) 941 mat_U = 0.0_dp 942 ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points)) 943 944 mat_SinvVSinvT = 0.0_dp 945 ! double the value nessary for 'A' to achieve good performance 946 lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes 947 ALLOCATE (work(lwork)) 948 work = 0.0_dp 949 ALLOCATE (iwork(8*num_integ_points)) 950 iwork = 0 951 ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes)) 952 mat_SinvVSinvSigma = 0.0_dp 953 ALLOCATE (vec_UTy(num_x_nodes)) 954 vec_UTy = 0.0_dp 955 956 max_error = 0.0_dp 957 958 ! loop over all omega frequency points 959 DO jquad = 1, num_integ_points 960 961 chi2_min_jquad = 100.0_dp 962 963 ! set the x-values logarithmically in the interval [Emin,Emax] 964 multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp)) 965 DO iii = 1, num_x_nodes 966 x_values(iii) = E_min*multiplicator**(iii - 1) 967 END DO 968 969 omega = omega_tj(jquad) 970 971 ! y=2x/(x^2+omega_k^2) 972 DO iii = 1, num_x_nodes 973! y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2+omega**2) 974 y_values(iii) = 2.0_dp*omega/((x_values(iii))**2 + omega**2) 975 END DO 976 977 ! calculate mat_A 978 DO jjj = 1, num_integ_points 979 DO iii = 1, num_x_nodes 980 mat_A(iii, jjj) = SIN(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj)) 981 END DO 982 END DO 983 984 ! Singular value decomposition of mat_A 985 CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, & 986 mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info) 987 988 CPASSERT(info == 0) 989 990 ! integration weights = V Sigma U^T y 991 ! 1) V*Sigma 992 DO jjj = 1, num_integ_points 993 DO iii = 1, num_integ_points 994 mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)/sing_values(jjj) 995 END DO 996 END DO 997 998 ! 2) U^T y 999 CALL DGEMM('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, y_values, num_x_nodes, & 1000 0.0_dp, vec_UTy, num_x_nodes) 1001 1002 ! 3) (V*Sigma) * (U^T y) 1003 CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, vec_UTy, & 1004 num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points) 1005 1006 weights_sin_tf_t_to_w(jquad, :) = tau_wj_work(:) 1007 1008 CALL calc_max_error_fit_tau_grid_with_sine(max_error, omega, tau_tj, tau_wj_work, x_values, & 1009 y_values, num_integ_points, num_x_nodes) 1010 1011 END DO ! jquad 1012 1013 DEALLOCATE (x_values, y_values, mat_A, tau_wj_work, work_array, sing_values, mat_U, mat_SinvVSinvT, & 1014 work, iwork, mat_SinvVSinvSigma, vec_UTy) 1015 1016 CALL timestop(handle) 1017 1018 END SUBROUTINE get_l_sq_wghts_sin_tf_t_to_w 1019 1020! ************************************************************************************************** 1021!> \brief ... 1022!> \param max_error ... 1023!> \param omega ... 1024!> \param tau_tj ... 1025!> \param tau_wj_work ... 1026!> \param x_values ... 1027!> \param y_values ... 1028!> \param num_integ_points ... 1029!> \param num_x_nodes ... 1030! ************************************************************************************************** 1031 PURE SUBROUTINE calc_max_error_fit_tau_grid_with_cosine(max_error, omega, tau_tj, tau_wj_work, x_values, & 1032 y_values, num_integ_points, num_x_nodes) 1033 1034 REAL(KIND=dp), INTENT(INOUT) :: max_error, omega 1035 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1036 INTENT(IN) :: tau_tj, tau_wj_work, x_values, y_values 1037 INTEGER, INTENT(IN) :: num_integ_points, num_x_nodes 1038 1039 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_max_error_fit_tau_grid_with_cosine', & 1040 routineP = moduleN//':'//routineN 1041 1042 INTEGER :: kkk 1043 REAL(KIND=dp) :: func_val, func_val_temp, max_error_tmp 1044 1045 max_error_tmp = 0.0_dp 1046 1047 DO kkk = 1, num_x_nodes 1048 1049 func_val = 0.0_dp 1050 1051 CALL eval_fit_func_tau_grid_cosine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega) 1052 1053 IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN 1054 max_error_tmp = ABS(y_values(kkk) - func_val) 1055 func_val_temp = func_val 1056 END IF 1057 1058 END DO 1059 1060 IF (max_error_tmp > max_error) THEN 1061 1062 max_error = max_error_tmp 1063 1064 END IF 1065 1066 END SUBROUTINE calc_max_error_fit_tau_grid_with_cosine 1067 1068! ************************************************************************************************** 1069!> \brief Evaluate fit function when calculating tau grid for cosine transform 1070!> \param func_val ... 1071!> \param x_value ... 1072!> \param num_integ_points ... 1073!> \param tau_tj ... 1074!> \param tau_wj_work ... 1075!> \param omega ... 1076! ************************************************************************************************** 1077 PURE SUBROUTINE eval_fit_func_tau_grid_cosine(func_val, x_value, num_integ_points, tau_tj, tau_wj_work, omega) 1078 1079 REAL(KIND=dp), INTENT(OUT) :: func_val 1080 REAL(KIND=dp), INTENT(IN) :: x_value 1081 INTEGER, INTENT(IN) :: num_integ_points 1082 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1083 INTENT(IN) :: tau_tj, tau_wj_work 1084 REAL(KIND=dp), INTENT(IN) :: omega 1085 1086 CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_fit_func_tau_grid_cosine', & 1087 routineP = moduleN//':'//routineN 1088 1089 INTEGER :: iii 1090 1091 func_val = 0.0_dp 1092 1093 DO iii = 1, num_integ_points 1094 1095 ! calculate value of the fit function 1096 func_val = func_val + tau_wj_work(iii)*COS(omega*tau_tj(iii))*EXP(-x_value*tau_tj(iii)) 1097 1098 END DO 1099 1100 END SUBROUTINE eval_fit_func_tau_grid_cosine 1101 1102! ************************************************************************************************** 1103!> \brief Evaluate fit function when calculating tau grid for sine transform 1104!> \param func_val ... 1105!> \param x_value ... 1106!> \param num_integ_points ... 1107!> \param tau_tj ... 1108!> \param tau_wj_work ... 1109!> \param omega ... 1110! ************************************************************************************************** 1111 PURE SUBROUTINE eval_fit_func_tau_grid_sine(func_val, x_value, num_integ_points, tau_tj, tau_wj_work, omega) 1112 1113 REAL(KIND=dp), INTENT(INOUT) :: func_val 1114 REAL(KIND=dp), INTENT(IN) :: x_value 1115 INTEGER, INTENT(in) :: num_integ_points 1116 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1117 INTENT(IN) :: tau_tj, tau_wj_work 1118 REAL(KIND=dp), INTENT(IN) :: omega 1119 1120 CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_fit_func_tau_grid_sine', & 1121 routineP = moduleN//':'//routineN 1122 1123 INTEGER :: iii 1124 1125 func_val = 0.0_dp 1126 1127 DO iii = 1, num_integ_points 1128 1129 ! calculate value of the fit function 1130 func_val = func_val + tau_wj_work(iii)*SIN(omega*tau_tj(iii))*EXP(-x_value*tau_tj(iii)) 1131 1132 END DO 1133 1134 END SUBROUTINE eval_fit_func_tau_grid_sine 1135 1136! ************************************************************************************************** 1137!> \brief ... 1138!> \param max_error ... 1139!> \param omega ... 1140!> \param tau_tj ... 1141!> \param tau_wj_work ... 1142!> \param x_values ... 1143!> \param y_values ... 1144!> \param num_integ_points ... 1145!> \param num_x_nodes ... 1146! ************************************************************************************************** 1147 PURE SUBROUTINE calc_max_error_fit_tau_grid_with_sine(max_error, omega, tau_tj, tau_wj_work, x_values, & 1148 y_values, num_integ_points, num_x_nodes) 1149 1150 REAL(KIND=dp), INTENT(INOUT) :: max_error, omega 1151 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1152 INTENT(IN) :: tau_tj, tau_wj_work, x_values, y_values 1153 INTEGER, INTENT(IN) :: num_integ_points, num_x_nodes 1154 1155 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_max_error_fit_tau_grid_with_sine', & 1156 routineP = moduleN//':'//routineN 1157 1158 INTEGER :: kkk 1159 REAL(KIND=dp) :: func_val, func_val_temp, max_error_tmp 1160 1161 max_error_tmp = 0.0_dp 1162 1163 DO kkk = 1, num_x_nodes 1164 1165 func_val = 0.0_dp 1166 1167 CALL eval_fit_func_tau_grid_sine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega) 1168 1169 IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN 1170 max_error_tmp = ABS(y_values(kkk) - func_val) 1171 func_val_temp = func_val 1172 END IF 1173 1174 END DO 1175 1176 IF (max_error_tmp > max_error) THEN 1177 1178 max_error = max_error_tmp 1179 1180 END IF 1181 1182 END SUBROUTINE calc_max_error_fit_tau_grid_with_sine 1183 1184! ************************************************************************************************** 1185!> \brief test the singular value decomposition for the computation of integration weights for the 1186!> Fourier transform between time and frequency grid in cubic-scaling RPA 1187!> \param nR ... 1188!> \param iw ... 1189! ************************************************************************************************** 1190 SUBROUTINE test_least_square_ft(nR, iw) 1191 INTEGER, INTENT(IN) :: nR, iw 1192 1193 INTEGER :: ierr, iR, jquad, num_integ_points 1194 REAL(KIND=dp) :: max_error, multiplicator, Rc, Rc_max 1195 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tau_tj, tau_wj, tj, wj, x_tw 1196 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: weights_cos_tf_t_to_w 1197 1198 Rc_max = 1.0E+7 1199 1200 multiplicator = Rc_max**(1.0_dp/(REAL(nR, KIND=dp) - 1.0_dp)) 1201 1202 DO num_integ_points = 1, 20 1203 1204 ALLOCATE (x_tw(2*num_integ_points)) 1205 x_tw = 0.0_dp 1206 ALLOCATE (tau_tj(num_integ_points)) 1207 tau_tj = 0.0_dp 1208 ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points)) 1209 weights_cos_tf_t_to_w = 0.0_dp 1210 ALLOCATE (tau_wj(num_integ_points)) 1211 tau_wj = 0.0_dp 1212 ALLOCATE (tj(num_integ_points)) 1213 tj = 0.0_dp 1214 ALLOCATE (wj(num_integ_points)) 1215 wj = 0.0_dp 1216 1217 DO iR = 0, nR - 1 1218 1219 Rc = 2.0_dp*multiplicator**iR 1220 1221 ierr = 0 1222 CALL get_rpa_minimax_coeff(num_integ_points, Rc, x_tw, ierr, print_warning=.FALSE.) 1223 1224 DO jquad = 1, num_integ_points 1225 tj(jquad) = x_tw(jquad) 1226 wj(jquad) = x_tw(jquad + num_integ_points) 1227 END DO 1228 1229 x_tw = 0.0_dp 1230 1231 CALL get_exp_minimax_coeff(num_integ_points, Rc, x_tw) 1232 1233 DO jquad = 1, num_integ_points 1234 tau_tj(jquad) = x_tw(jquad)/2.0_dp 1235 tau_wj(jquad) = x_tw(jquad + num_integ_points)/2.0_dp 1236 END DO 1237 1238 CALL get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, tj, & 1239 1.0_dp, Rc, max_error, 200) 1240 1241 IF (iw > 0) THEN 1242 WRITE (iw, '(T2, I3, F12.1, ES12.3)') num_integ_points, Rc, max_error 1243 END IF 1244 1245 END DO 1246 1247 DEALLOCATE (x_tw, tau_tj, weights_cos_tf_t_to_w, tau_wj, wj, tj) 1248 1249 END DO 1250 1251 END SUBROUTINE test_least_square_ft 1252 1253! ************************************************************************************************** 1254!> \brief ... 1255!> \param num_integ_points ... 1256!> \param tau_tj ... 1257!> \param weights_cos_tf_w_to_t ... 1258!> \param omega_tj ... 1259!> \param E_min ... 1260!> \param E_max ... 1261!> \param max_error ... 1262!> \param num_points_per_magnitude ... 1263! ************************************************************************************************** 1264 SUBROUTINE get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, omega_tj, & 1265 E_min, E_max, max_error, num_points_per_magnitude) 1266 1267 INTEGER, INTENT(IN) :: num_integ_points 1268 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1269 INTENT(IN) :: tau_tj 1270 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 1271 INTENT(INOUT) :: weights_cos_tf_w_to_t 1272 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1273 INTENT(IN) :: omega_tj 1274 REAL(KIND=dp), INTENT(IN) :: E_min, E_max 1275 REAL(KIND=dp), INTENT(INOUT) :: max_error 1276 INTEGER, INTENT(IN) :: num_points_per_magnitude 1277 1278 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_l_sq_wghts_cos_tf_w_to_t', & 1279 routineP = moduleN//':'//routineN 1280 1281 INTEGER :: handle, iii, info, jjj, jquad, lwork, & 1282 num_x_nodes 1283 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork 1284 REAL(KIND=dp) :: chi2_min_jquad, multiplicator, omega, & 1285 tau, x_value 1286 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: omega_wj_work, sing_values, vec_UTy, & 1287 work, work_array, x_values, y_values 1288 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, & 1289 mat_SinvVSinvT, mat_U 1290 1291 CALL timeset(routineN, handle) 1292 1293 ! take num_points_per_magnitude points per 10-interval 1294 num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude 1295 1296 ! take at least as many x points as integration points to have clear 1297 ! input for the singular value decomposition 1298 num_x_nodes = MAX(num_x_nodes, num_integ_points) 1299 1300 ALLOCATE (x_values(num_x_nodes)) 1301 x_values = 0.0_dp 1302 ALLOCATE (y_values(num_x_nodes)) 1303 y_values = 0.0_dp 1304 ALLOCATE (mat_A(num_x_nodes, num_integ_points)) 1305 mat_A = 0.0_dp 1306 ALLOCATE (omega_wj_work(num_integ_points)) 1307 omega_wj_work = 0.0_dp 1308 ALLOCATE (work_array(2*num_integ_points)) 1309 work_array = 0.0_dp 1310 ALLOCATE (sing_values(num_integ_points)) 1311 sing_values = 0.0_dp 1312 ALLOCATE (mat_U(num_x_nodes, num_x_nodes)) 1313 mat_U = 0.0_dp 1314 ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points)) 1315 1316 mat_SinvVSinvT = 0.0_dp 1317 ! double the value nessary for 'A' to achieve good performance 1318 lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes 1319 ALLOCATE (work(lwork)) 1320 work = 0.0_dp 1321 ALLOCATE (iwork(8*num_integ_points)) 1322 iwork = 0 1323 ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes)) 1324 mat_SinvVSinvSigma = 0.0_dp 1325 ALLOCATE (vec_UTy(num_x_nodes)) 1326 vec_UTy = 0.0_dp 1327 1328 ! set the x-values logarithmically in the interval [Emin,Emax] 1329 multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp)) 1330 DO iii = 1, num_x_nodes 1331 x_values(iii) = E_min*multiplicator**(iii - 1) 1332 END DO 1333 1334 max_error = 0.0_dp 1335 1336 ! loop over all tau time points 1337 DO jquad = 1, num_integ_points 1338 1339 chi2_min_jquad = 100.0_dp 1340 1341 tau = tau_tj(jquad) 1342 1343 ! y=exp(-x*|tau_k|) 1344 DO iii = 1, num_x_nodes 1345 y_values(iii) = EXP(-x_values(iii)*tau) 1346 END DO 1347 1348 ! calculate mat_A 1349 DO jjj = 1, num_integ_points 1350 DO iii = 1, num_x_nodes 1351 omega = omega_tj(jjj) 1352 x_value = x_values(iii) 1353 mat_A(iii, jjj) = COS(tau*omega)*2.0_dp*x_value/(x_value**2 + omega**2) 1354 END DO 1355 END DO 1356 1357 ! Singular value decomposition of mat_A 1358 CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, & 1359 mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info) 1360 1361 CPASSERT(info == 0) 1362 1363 ! integration weights = V Sigma U^T y 1364 ! 1) V*Sigma 1365 DO jjj = 1, num_integ_points 1366 DO iii = 1, num_integ_points 1367 mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)/sing_values(jjj) 1368 END DO 1369 END DO 1370 1371 ! 2) U^T y 1372 CALL DGEMM('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, y_values, num_x_nodes, & 1373 0.0_dp, vec_UTy, num_x_nodes) 1374 1375 ! 3) (V*Sigma) * (U^T y) 1376 CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, vec_UTy, & 1377 num_x_nodes, 0.0_dp, omega_wj_work, num_integ_points) 1378 1379 weights_cos_tf_w_to_t(jquad, :) = omega_wj_work(:) 1380 1381 CALL calc_max_error_fit_omega_grid_with_cosine(max_error, tau, omega_tj, omega_wj_work, x_values, & 1382 y_values, num_integ_points, num_x_nodes) 1383 1384 END DO ! jquad 1385 1386 DEALLOCATE (x_values, y_values, mat_A, omega_wj_work, work_array, sing_values, mat_U, mat_SinvVSinvT, & 1387 work, iwork, mat_SinvVSinvSigma, vec_UTy) 1388 1389 CALL timestop(handle) 1390 1391 END SUBROUTINE get_l_sq_wghts_cos_tf_w_to_t 1392 1393! ************************************************************************************************** 1394!> \brief ... 1395!> \param max_error ... 1396!> \param tau ... 1397!> \param omega_tj ... 1398!> \param omega_wj_work ... 1399!> \param x_values ... 1400!> \param y_values ... 1401!> \param num_integ_points ... 1402!> \param num_x_nodes ... 1403! ************************************************************************************************** 1404 SUBROUTINE calc_max_error_fit_omega_grid_with_cosine(max_error, tau, omega_tj, omega_wj_work, x_values, & 1405 y_values, num_integ_points, num_x_nodes) 1406 1407 REAL(KIND=dp), INTENT(INOUT) :: max_error 1408 REAL(KIND=dp), INTENT(IN) :: tau 1409 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1410 INTENT(IN) :: omega_tj, omega_wj_work, x_values, & 1411 y_values 1412 INTEGER, INTENT(IN) :: num_integ_points, num_x_nodes 1413 1414 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_max_error_fit_omega_grid_with_cosine', & 1415 routineP = moduleN//':'//routineN 1416 1417 INTEGER :: handle, kkk 1418 REAL(KIND=dp) :: func_val, func_val_temp, max_error_tmp 1419 1420 CALL timeset(routineN, handle) 1421 1422 max_error_tmp = 0.0_dp 1423 1424 DO kkk = 1, num_x_nodes 1425 1426 func_val = 0.0_dp 1427 1428 CALL eval_fit_func_omega_grid_cosine(func_val, x_values(kkk), num_integ_points, omega_tj, omega_wj_work, tau) 1429 1430 IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN 1431 max_error_tmp = ABS(y_values(kkk) - func_val) 1432 func_val_temp = func_val 1433 END IF 1434 1435 END DO 1436 1437 IF (max_error_tmp > max_error) THEN 1438 1439 max_error = max_error_tmp 1440 1441 END IF 1442 1443 CALL timestop(handle) 1444 1445 END SUBROUTINE calc_max_error_fit_omega_grid_with_cosine 1446 1447! ************************************************************************************************** 1448!> \brief ... 1449!> \param func_val ... 1450!> \param x_value ... 1451!> \param num_integ_points ... 1452!> \param omega_tj ... 1453!> \param omega_wj_work ... 1454!> \param tau ... 1455! ************************************************************************************************** 1456 PURE SUBROUTINE eval_fit_func_omega_grid_cosine(func_val, x_value, num_integ_points, omega_tj, omega_wj_work, tau) 1457 REAL(KIND=dp), INTENT(OUT) :: func_val 1458 REAL(KIND=dp), INTENT(IN) :: x_value 1459 INTEGER, INTENT(IN) :: num_integ_points 1460 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1461 INTENT(IN) :: omega_tj, omega_wj_work 1462 REAL(KIND=dp), INTENT(IN) :: tau 1463 1464 CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_fit_func_omega_grid_cosine', & 1465 routineP = moduleN//':'//routineN 1466 1467 INTEGER :: iii 1468 REAL(KIND=dp) :: omega 1469 1470 func_val = 0.0_dp 1471 1472 DO iii = 1, num_integ_points 1473 1474 ! calculate value of the fit function 1475 omega = omega_tj(iii) 1476 func_val = func_val + omega_wj_work(iii)*COS(tau*omega)*2.0_dp*x_value/(x_value**2 + omega**2) 1477 1478 END DO 1479 1480 END SUBROUTINE eval_fit_func_omega_grid_cosine 1481 1482! ************************************************************************************************** 1483!> \brief ... 1484!> \param qs_env ... 1485!> \param para_env ... 1486!> \param gap ... 1487!> \param max_eig_diff ... 1488!> \param e_fermi ... 1489! ************************************************************************************************** 1490 SUBROUTINE gap_and_max_eig_diff_kpoints(qs_env, para_env, gap, max_eig_diff, e_fermi) 1491 1492 TYPE(qs_environment_type), POINTER :: qs_env 1493 TYPE(cp_para_env_type), POINTER :: para_env 1494 REAL(KIND=dp), INTENT(OUT) :: gap, max_eig_diff, e_fermi 1495 1496 CHARACTER(LEN=*), PARAMETER :: routineN = 'gap_and_max_eig_diff_kpoints', & 1497 routineP = moduleN//':'//routineN 1498 1499 INTEGER :: handle, homo, ikpgr, ispin, kplocal, & 1500 nmo, nspin 1501 INTEGER, DIMENSION(2) :: kp_range 1502 REAL(KIND=dp) :: e_homo, e_homo_temp, e_lumo, e_lumo_temp 1503 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_homo_array, e_lumo_array, & 1504 max_eig_diff_array 1505 REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues 1506 TYPE(kpoint_env_type), POINTER :: kp 1507 TYPE(kpoint_type), POINTER :: kpoint 1508 TYPE(mo_set_type), POINTER :: mo_set 1509 1510 CALL timeset(routineN, handle) 1511 1512 CALL get_qs_env(qs_env, & 1513 kpoints=kpoint) 1514 1515 mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)%mo_set 1516 CALL get_mo_set(mo_set, nmo=nmo) 1517 1518 CALL get_kpoint_info(kpoint, kp_range=kp_range) 1519 kplocal = kp_range(2) - kp_range(1) + 1 1520 1521 gap = 1000.0_dp 1522 max_eig_diff = 0.0_dp 1523 e_homo = -1000.0_dp 1524 e_lumo = 1000.0_dp 1525 1526 DO ikpgr = 1, kplocal 1527 kp => kpoint%kp_env(ikpgr)%kpoint_env 1528 nspin = SIZE(kp%mos, 2) 1529 DO ispin = 1, nspin 1530 mo_set => kp%mos(1, ispin)%mo_set 1531 CALL get_mo_set(mo_set, eigenvalues=eigenvalues, homo=homo) 1532 e_homo_temp = eigenvalues(homo) 1533 e_lumo_temp = eigenvalues(homo + 1) 1534 1535 IF (e_homo_temp > e_homo) e_homo = e_homo_temp 1536 IF (e_lumo_temp < e_lumo) e_lumo = e_lumo_temp 1537 IF (eigenvalues(nmo) - eigenvalues(1) > max_eig_diff) max_eig_diff = eigenvalues(nmo) - eigenvalues(1) 1538 1539 END DO 1540 END DO 1541 1542 ALLOCATE (e_homo_array(0:para_env%num_pe - 1)) 1543 e_homo_array = 0.0_dp 1544 e_homo_array(para_env%mepos) = e_homo 1545 1546 ALLOCATE (e_lumo_array(0:para_env%num_pe - 1)) 1547 e_lumo_array = 0.0_dp 1548 e_lumo_array(para_env%mepos) = e_lumo 1549 1550 ALLOCATE (max_eig_diff_array(0:para_env%num_pe - 1)) 1551 max_eig_diff_array = 0.0_dp 1552 max_eig_diff_array(para_env%mepos) = max_eig_diff 1553 1554 CALL mp_sum(e_homo_array, para_env%group) 1555 1556 CALL mp_sum(e_lumo_array, para_env%group) 1557 1558 CALL mp_sum(max_eig_diff_array, para_env%group) 1559 1560 gap = MINVAL(e_lumo_array) - MAXVAL(e_homo_array) 1561 1562 e_fermi = (MAXVAL(e_homo_array) + MINVAL(e_lumo_array))/2.0_dp 1563 1564 max_eig_diff = MAXVAL(max_eig_diff_array) 1565 1566 DEALLOCATE (max_eig_diff_array, e_homo_array, e_lumo_array) 1567 1568 CALL timestop(handle) 1569 1570 END SUBROUTINE 1571 1572END MODULE mp2_weights 1573