diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b856eb9b04..f23908c6cd 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -530,6 +530,7 @@ list( qs_condnum.F qs_core_energies.F qs_core_hamiltonian.F + qs_core_matrices.F qs_dcdr_ao.F qs_dcdr.F qs_dcdr_utils.F diff --git a/src/core_ppl.F b/src/core_ppl.F index 2ea6d34757..a0e425e3c7 100644 --- a/src/core_ppl.F +++ b/src/core_ppl.F @@ -168,6 +168,10 @@ CONTAINS do_dR = PRESENT(deltaR) doat = PRESENT(atcore) + IF ((calculate_forces .OR. doat) .AND. do_dR) THEN + CPABORT("core_ppl: incompatible options") + END IF + MARK_USED(int_8) ! Use internal integral routine for local ECP terms or use libgrrp diff --git a/src/core_ppnl.F b/src/core_ppnl.F index 1e6593be42..2e2f0a00a5 100644 --- a/src/core_ppnl.F +++ b/src/core_ppnl.F @@ -162,6 +162,9 @@ CONTAINS IF (PRESENT(deltaR)) do_dR = .TRUE. doat = .FALSE. IF (PRESENT(atcore)) doat = .TRUE. + IF ((calculate_forces .OR. doat) .AND. do_dR) THEN + CPABORT("core_ppl: incompatible options") + END IF IF (calculate_forces) THEN CALL timeset(routineN//"_forces", handle) diff --git a/src/ec_external.F b/src/ec_external.F index 0b9a1c2d8f..5b2bf7fc18 100644 --- a/src/ec_external.F +++ b/src/ec_external.F @@ -12,11 +12,7 @@ !> \author JGH ! ************************************************************************************************** MODULE ec_external - USE atomic_kind_types, ONLY: atomic_kind_type USE cell_types, ONLY: cell_type - USE core_ae, ONLY: build_core_ae - USE core_ppl, ONLY: build_core_ppl - USE core_ppnl, ONLY: build_core_ppnl USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_api, ONLY: dbcsr_add,& dbcsr_copy,& @@ -50,14 +46,13 @@ MODULE ec_external USE mathlib, ONLY: det_3x3 USE message_passing, ONLY: mp_para_env_type USE parallel_gemm_api, ONLY: parallel_gemm - USE particle_types, ONLY: particle_type USE physcon, ONLY: pascal USE qs_core_energies, ONLY: calculate_ptrace + USE qs_core_matrices, ONLY: core_matrices,& + kinetic_energy_matrix USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_force_types, ONLY: qs_force_type - USE qs_kind_types, ONLY: qs_kind_type - USE qs_kinetic, ONLY: build_kinetic_matrix USE qs_ks_types, ONLY: qs_ks_env_type USE qs_mo_types, ONLY: deallocate_mo_set,& get_mo_set,& @@ -361,23 +356,16 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_debug_force' INTEGER :: handle, iounit, nder, nimages - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index LOGICAL :: calculate_forces, debug_forces, & debug_stress, use_virial - REAL(KIND=dp) :: eps_ppnl, fconv - REAL(KIND=dp), DIMENSION(3) :: fodeb + REAL(KIND=dp) :: fconv REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot - TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm TYPE(dft_control_type), POINTER :: dft_control TYPE(mp_para_env_type), POINTER :: para_env TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl - TYPE(particle_type), DIMENSION(:), POINTER :: particle_set - TYPE(qs_force_type), DIMENSION(:), POINTER :: force - TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set - TYPE(qs_ks_env_type), POINTER :: ks_env + POINTER :: sab_orb TYPE(virial_type), POINTER :: virial CALL timeset(routineN, handle) @@ -389,12 +377,9 @@ CONTAINS calculate_forces = .TRUE. ! no k-points possible - NULLIFY (cell, dft_control, force, ks_env, para_env, virial) CALL get_qs_env(qs_env=qs_env, & cell=cell, & dft_control=dft_control, & - force=force, & - ks_env=ks_env, & para_env=para_env, & virial=virial) nimages = dft_control%nimages @@ -415,116 +400,22 @@ CONTAINS CPABORT("GAPW not implemented") END IF - ! get neighbor lists - NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl) - CALL get_qs_env(qs_env=qs_env, & - sab_orb=sab_orb, sac_ae=sac_ae, & - sac_ppl=sac_ppl, sap_ppnl=sap_ppnl) - ! initialize src matrix NULLIFY (scrm) CALL dbcsr_allocate_matrix_set(scrm, 1, 1) ALLOCATE (scrm(1, 1)%matrix) CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_p(1, 1)%matrix) + CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb) CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb) - nder = 1 - IF (SIZE(matrix_p, 1) == 2) THEN - CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, & - alpha_scalar=1.0_dp, beta_scalar=1.0_dp) - END IF - ! kinetic energy - IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ekinetic - CALL build_kinetic_matrix(ks_env, matrixkp_t=scrm, & - matrix_name="KINETIC ENERGY MATRIX", & - basis_type="ORB", & - sab_nl=sab_orb, calculate_forces=.TRUE., & - matrixkp_p=matrix_p) - IF (debug_forces) THEN - fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dT ", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ekinetic - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pout*dT', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - IF (SIZE(matrix_p, 1) == 2) THEN - CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, & - alpha_scalar=1.0_dp, beta_scalar=-1.0_dp) - END IF + CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm, matrix_p=matrix_p, & + calculate_forces=calculate_forces, & + debug_forces=debug_forces, debug_stress=debug_stress) - ! compute the ppl contribution to the core hamiltonian - NULLIFY (atomic_kind_set, particle_set, qs_kind_set) - CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, & - atomic_kind_set=atomic_kind_set) - - IF (ASSOCIATED(sac_ae)) THEN - IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%all_potential(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl - CALL build_core_ae(scrm, matrix_p, force, & - virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, & - nimages, cell_to_index, "ORB") - IF (calculate_forces .AND. debug_forces) THEN - fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_AE ", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ppl - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pout*dH_AE ', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - END IF - - IF (ASSOCIATED(sac_ppl)) THEN - IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl - CALL build_core_ppl(scrm, matrix_p, force, & - virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & - nimages, cell_to_index, "ORB") - IF (calculate_forces .AND. debug_forces) THEN - fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPL ", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ppl - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pout*dH_PPL', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - END IF - - ! compute the ppnl contribution to the core hamiltonian *** - eps_ppnl = dft_control%qs_control%eps_ppnl - IF (ASSOCIATED(sap_ppnl)) THEN - IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl - CALL build_core_ppnl(scrm, matrix_p, force, & - virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, & - sab_orb, sap_ppnl, eps_ppnl, & - nimages, cell_to_index, "ORB") - IF (calculate_forces .AND. debug_forces) THEN - fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPNL", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ppnl - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pout*dH_PPNL', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - END IF + nder = 1 + CALL core_matrices(qs_env, scrm, matrix_p, calculate_forces, nder, & + debug_forces=debug_forces, debug_stress=debug_stress) IF (debug_stress .AND. use_virial) THEN stdeb = fconv*(virial%pv_virial - sttot) diff --git a/src/ed_analysis.F b/src/ed_analysis.F index d7d96b9175..6d00a5216c 100644 --- a/src/ed_analysis.F +++ b/src/ed_analysis.F @@ -72,7 +72,7 @@ MODULE ed_analysis USE qs_core_energies, ONLY: calculate_ecore_alpha,& calculate_ecore_overlap,& calculate_ecore_self - USE qs_core_hamiltonian, ONLY: core_matrices + USE qs_core_matrices, ONLY: core_matrices USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot USE qs_dispersion_types, ONLY: qs_dispersion_type USE qs_energy_types, ONLY: qs_energy_type diff --git a/src/energy_corrections.F b/src/energy_corrections.F index 9385cbd84e..d0e87a070c 100644 --- a/src/energy_corrections.F +++ b/src/energy_corrections.F @@ -26,9 +26,6 @@ MODULE energy_corrections cite_reference USE cell_types, ONLY: cell_type,& pbc - USE core_ae, ONLY: build_core_ae - USE core_ppl, ONLY: build_core_ppl - USE core_ppnl, ONLY: build_core_ppnl USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_api, ONLY: & @@ -142,6 +139,8 @@ MODULE energy_corrections USE qs_collocate_density, ONLY: calculate_rho_elec USE qs_core_energies, ONLY: calculate_ecore_overlap,& calculate_ptrace + USE qs_core_matrices, ONLY: core_matrices,& + kinetic_energy_matrix USE qs_density_matrices, ONLY: calculate_density_matrix,& calculate_w_matrix USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot @@ -1322,23 +1321,21 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_core_hamiltonian' + CHARACTER(LEN=default_string_length) :: basis_type INTEGER :: handle, nder, nimages - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index LOGICAL :: calculate_forces, use_virial - REAL(KIND=dp) :: eps_ppnl TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(dft_control_type), POINTER :: dft_control TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl TYPE(particle_type), DIMENSION(:), POINTER :: particle_set - TYPE(qs_force_type), DIMENSION(:), POINTER :: force TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(qs_ks_env_type), POINTER :: ks_env - TYPE(virial_type), POINTER :: virial CALL timeset(routineN, handle) - NULLIFY (atomic_kind_set, cell_to_index, dft_control, ks_env, particle_set, qs_kind_set, virial) + NULLIFY (atomic_kind_set, dft_control, ks_env, particle_set, & + qs_kind_set) CALL get_qs_env(qs_env=qs_env, & atomic_kind_set=atomic_kind_set, & @@ -1369,16 +1366,18 @@ CONTAINS sac_ppl => ec_env%sac_ppl sap_ppnl => ec_env%sap_ppnl + basis_type = "HARRIS" + nder = 0 ! Overlap and kinetic energy matrices CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, & matrix_name="OVERLAP MATRIX", & - basis_type_a="HARRIS", & - basis_type_b="HARRIS", & + basis_type_a=basis_type, & + basis_type_b=basis_type, & sab_nl=sab_orb) CALL build_kinetic_matrix(ks_env, matrixkp_t=ec_env%matrix_t, & matrix_name="KINETIC ENERGY MATRIX", & - basis_type="HARRIS", & + basis_type=basis_type, & sab_nl=sab_orb) ! initialize H matrix @@ -1391,29 +1390,8 @@ CONTAINS CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, ec_env%matrix_t(1, 1)%matrix, & keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX") - ! Possible AE contributions (also ECP) - IF (ASSOCIATED(sac_ae)) THEN - CALL build_core_ae(ec_env%matrix_h, ec_env%matrix_p, force, & - virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, & - nimages, cell_to_index, "HARRIS") - END IF - ! compute the ppl contribution to the core hamiltonian - IF (ASSOCIATED(sac_ppl)) THEN - CALL build_core_ppl(ec_env%matrix_h, ec_env%matrix_p, force, & - virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & - nimages, cell_to_index, "HARRIS") - END IF - - ! compute the ppnl contribution to the core hamiltonian *** - eps_ppnl = dft_control%qs_control%eps_ppnl - IF (ASSOCIATED(sap_ppnl)) THEN - CALL build_core_ppnl(ec_env%matrix_h, ec_env%matrix_p, force, & - virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, & - eps_ppnl, nimages, cell_to_index, "HARRIS") - END IF + CALL core_matrices(qs_env, ec_env%matrix_h, ec_env%matrix_p, calculate_forces, nder, & + ec_env=ec_env, ec_env_matrices=.TRUE., basis_type=basis_type) ! External field (nonperiodic case) ec_env%efield_nuclear = 0.0_dp @@ -1660,24 +1638,21 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_core_hamiltonian_force' + CHARACTER(LEN=default_string_length) :: basis_type INTEGER :: handle, iounit, nder, nimages - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index LOGICAL :: calculate_forces, debug_forces, & debug_stress, use_virial - REAL(KIND=dp) :: eps_ppnl, fconv + REAL(KIND=dp) :: fconv REAL(KIND=dp), DIMENSION(3) :: fodeb REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot - TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm TYPE(dft_control_type), POINTER :: dft_control TYPE(mp_para_env_type), POINTER :: para_env TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl - TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + POINTER :: sab_orb TYPE(qs_force_type), DIMENSION(:), POINTER :: force - TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(virial_type), POINTER :: virial @@ -1695,6 +1670,8 @@ CONTAINS calculate_forces = .TRUE. + basis_type = "HARRIS" + ! no k-points possible NULLIFY (cell, dft_control, force, ks_env, para_env, virial) CALL get_qs_env(qs_env=qs_env, & @@ -1725,11 +1702,7 @@ CONTAINS END IF ! get neighbor lists, we need the full sab_orb list from the ec_env - NULLIFY (sab_orb, sac_ppl, sap_ppnl) sab_orb => ec_env%sab_orb - sac_ae => ec_env%sac_ae - sac_ppl => ec_env%sac_ppl - sap_ppnl => ec_env%sap_ppnl ! initialize src matrix NULLIFY (scrm) @@ -1740,8 +1713,6 @@ CONTAINS nder = 1 IF (SIZE(matrix_p, 1) == 2) THEN - CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, & - alpha_scalar=1.0_dp, beta_scalar=1.0_dp) CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, & alpha_scalar=1.0_dp, beta_scalar=1.0_dp) END IF @@ -1751,8 +1722,8 @@ CONTAINS IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, & matrix_name="OVERLAP MATRIX", & - basis_type_a="HARRIS", & - basis_type_b="HARRIS", & + basis_type_a=basis_type, & + basis_type_b=basis_type, & sab_nl=sab_orb, calculate_forces=.TRUE., & matrixkp_p=matrix_w) @@ -1760,103 +1731,23 @@ CONTAINS fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3) CALL para_env%sum(fodeb) IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wout*dS ", fodeb - fodeb(1:3) = force(1)%kinetic(1:3, 1) END IF IF (debug_stress .AND. use_virial) THEN stdeb = fconv*(virial%pv_overlap - stdeb) CALL para_env%sum(stdeb) IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & 'STRESS| Wout*dS', one_third_sum_diag(stdeb), det_3x3(stdeb) - stdeb = virial%pv_ekinetic - END IF - CALL build_kinetic_matrix(ks_env, matrixkp_t=scrm, & - matrix_name="KINETIC ENERGY MATRIX", & - basis_type="HARRIS", & - sab_nl=sab_orb, calculate_forces=.TRUE., & - matrixkp_p=matrix_p) - IF (debug_forces) THEN - fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dT ", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ekinetic - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pout*dT', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - IF (SIZE(matrix_p, 1) == 2) THEN - CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, & - alpha_scalar=1.0_dp, beta_scalar=-1.0_dp) END IF - ! compute the ppl contribution to the core hamiltonian - NULLIFY (atomic_kind_set, particle_set, qs_kind_set) - CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, & - atomic_kind_set=atomic_kind_set) + CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm, matrix_p=matrix_p, & + calculate_forces=.TRUE., & + sab_orb=sab_orb, & + basis_type=basis_type, & + debug_forces=debug_forces, debug_stress=debug_stress) - IF (ASSOCIATED(sac_ae)) THEN - IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%all_potential(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl - CALL build_core_ae(scrm, matrix_p, force, & - virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, & - nimages, cell_to_index, "HARRIS") - IF (calculate_forces .AND. debug_forces) THEN - fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dHae ", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ppl - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pout*dHae ', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - END IF - - IF (ASSOCIATED(sac_ppl)) THEN - IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl - CALL build_core_ppl(scrm, matrix_p, force, & - virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & - nimages, cell_to_index, "HARRIS") - IF (calculate_forces .AND. debug_forces) THEN - fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPL ", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ppl - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pout*dH_PPL', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - END IF - - ! compute the ppnl contribution to the core hamiltonian *** - eps_ppnl = dft_control%qs_control%eps_ppnl - IF (ASSOCIATED(sap_ppnl)) THEN - IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl - CALL build_core_ppnl(scrm, matrix_p, force, & - virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, & - sab_orb, sap_ppnl, eps_ppnl, & - nimages, cell_to_index, "HARRIS") - IF (calculate_forces .AND. debug_forces) THEN - fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPNL", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ppnl - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pout*dH_PPNL', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - END IF + CALL core_matrices(qs_env, scrm, matrix_p, calculate_forces, nder, & + ec_env=ec_env, ec_env_matrices=.FALSE., basis_type=basis_type, & + debug_forces=debug_forces, debug_stress=debug_stress) ! External field (nonperiodic case) ec_env%efield_nuclear = 0.0_dp diff --git a/src/mp2_cphf.F b/src/mp2_cphf.F index e07b6a8d8c..5dae2e540a 100644 --- a/src/mp2_cphf.F +++ b/src/mp2_cphf.F @@ -15,11 +15,7 @@ MODULE mp2_cphf USE admm_methods, ONLY: admm_projection_derivative USE admm_types, ONLY: admm_type,& get_admm_env - USE atomic_kind_types, ONLY: atomic_kind_type USE cell_types, ONLY: cell_type - USE core_ae, ONLY: build_core_ae - USE core_ppl, ONLY: build_core_ppl - USE core_ppnl, ONLY: build_core_ppnl USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_api, ONLY: dbcsr_add,& @@ -74,7 +70,6 @@ MODULE mp2_cphf USE mp2_types, ONLY: mp2_type,& ri_rpa_method_gpw USE parallel_gemm_api, ONLY: parallel_gemm - USE particle_types, ONLY: particle_type USE pw_env_types, ONLY: pw_env_get,& pw_env_type USE pw_methods, ONLY: pw_axpy,& @@ -90,6 +85,8 @@ MODULE mp2_cphf USE pw_types, ONLY: pw_c1d_gs_type,& pw_r3d_rs_type USE qs_2nd_kernel_ao, ONLY: apply_2nd_order_kernel + USE qs_core_matrices, ONLY: core_matrices,& + kinetic_energy_matrix USE qs_density_matrices, ONLY: calculate_whz_matrix USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot USE qs_dispersion_types, ONLY: qs_dispersion_type @@ -103,8 +100,6 @@ MODULE mp2_cphf zero_qs_force USE qs_integrate_potential, ONLY: integrate_v_core_rspace,& integrate_v_rspace - USE qs_kind_types, ONLY: qs_kind_type - USE qs_kinetic, ONLY: build_kinetic_matrix USE qs_ks_reference, ONLY: ks_ref_potential USE qs_ks_types, ONLY: qs_ks_env_type,& set_ks_env @@ -1291,13 +1286,11 @@ CONTAINS INTEGER :: alpha, beta, handle, idir, iounit, & ispin, nimages, nocc, nspins INTEGER, DIMENSION(3) :: comp - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index LOGICAL :: do_exx, do_hfx, use_virial REAL(KIND=dp) :: e_dummy, e_hartree, e_xc, ehartree, exc REAL(KIND=dp), DIMENSION(3) :: deb REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_virial TYPE(admm_type), POINTER :: admm_env - TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:), & @@ -1305,15 +1298,14 @@ CONTAINS TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_p_mp2, & matrix_p_mp2_admm, matrix_s, rho1, & rho_ao, rho_ao_aux, scrm - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp, scrm_kp + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, rho_ao_kp, scrm_kp TYPE(dft_control_type), POINTER :: dft_control TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data TYPE(linres_control_type), POINTER :: linres_control TYPE(mo_set_type), DIMENSION(:), POINTER :: mos TYPE(mp_para_env_type), POINTER :: para_env TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl - TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + POINTER :: sab_orb TYPE(pw_c1d_gs_type) :: pot_g, rho_tot_g, temp_pw_g TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: dvg TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_mp2_g @@ -1328,7 +1320,6 @@ CONTAINS TYPE(qs_dispersion_type), POINTER :: dispersion_env TYPE(qs_energy_type), POINTER :: energy TYPE(qs_force_type), DIMENSION(:), POINTER :: force - TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(qs_p_env_type), POINTER :: p_env TYPE(qs_rho_type), POINTER :: rho, rho_aux @@ -1373,7 +1364,6 @@ CONTAINS nimages = dft_control%nimages CPASSERT(nimages == 1) - NULLIFY (cell_to_index) p_env => qs_env%mp2_env%ri_grad%p_env @@ -1391,101 +1381,17 @@ CONTAINS CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp) END DO - IF (nspins == 2) THEN - CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, 1.0_dp, 1.0_dp) - END IF + ! kinetic energy + NULLIFY (scrm_kp) + matrix_p(1:nspins, 1:1) => rho_ao(1:nspins) + CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm_kp, matrix_p=matrix_p, & + calculate_forces=.TRUE., & + debug_forces=debug_forces, debug_stress=debug_forces) + CALL dbcsr_deallocate_matrix_set(scrm_kp) - ! Kinetic energy matrix - NULLIFY (scrm) - IF (debug_forces) THEN - deb(1:3) = force(1)%kinetic(1:3, 1) - IF (use_virial) e_dummy = third_tr(virial%pv_virial) - END IF - CALL build_kinetic_matrix(ks_env, matrix_t=scrm, & - matrix_name="KINETIC ENERGY MATRIX", & - basis_type="ORB", & - sab_nl=sab_orb, calculate_forces=.TRUE., & - matrix_p=rho_ao(1)%matrix) - IF (debug_forces) THEN - deb(1:3) = force(1)%kinetic(1:3, 1) - deb(1:3) - CALL para_env%sum(deb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT ", deb - IF (use_virial) THEN - e_dummy = third_tr(virial%pv_virial) - e_dummy - CALL para_env%sum(e_dummy) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: T ", e_dummy - END IF - END IF - CALL dbcsr_deallocate_matrix_set(scrm) - - IF (nspins == 2) THEN - CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, 1.0_dp, -1.0_dp) - END IF - - ! Add pseudo potential terms scrm_kp(1:nspins, 1:1) => matrix_ks(1:nspins) - CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, & - atomic_kind_set=atomic_kind_set, sac_ae=sac_ae, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl) - IF (ASSOCIATED(sac_ae)) THEN - IF (debug_forces) THEN - deb(1:3) = force(1)%all_potential(1:3, 1) - IF (use_virial) e_dummy = third_tr(virial%pv_virial) - END IF - CALL build_core_ae(scrm_kp, rho_ao_kp, force, & - virial, .TRUE., use_virial, 1, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, & - nimages, cell_to_index, basis_type="ORB") - IF (debug_forces) THEN - deb(1:3) = force(1)%all_potential(1:3, 1) - deb(1:3) - CALL para_env%sum(deb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae ", deb - IF (use_virial) THEN - e_dummy = third_tr(virial%pv_virial) - e_dummy - CALL para_env%sum(e_dummy) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hae ", e_dummy - END IF - END IF - END IF - IF (ASSOCIATED(sac_ppl)) THEN - IF (debug_forces) THEN - deb(1:3) = force(1)%gth_ppl(1:3, 1) - IF (use_virial) e_dummy = third_tr(virial%pv_virial) - END IF - CALL build_core_ppl(scrm_kp, rho_ao_kp, force, & - virial, .TRUE., use_virial, 1, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & - nimages=1, basis_type="ORB") - IF (debug_forces) THEN - deb(1:3) = force(1)%gth_ppl(1:3, 1) - deb(1:3) - CALL para_env%sum(deb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl ", deb - IF (use_virial) THEN - e_dummy = third_tr(virial%pv_virial) - e_dummy - CALL para_env%sum(e_dummy) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hppl ", e_dummy - END IF - END IF - END IF - IF (ASSOCIATED(sap_ppnl)) THEN - IF (debug_forces) THEN - deb(1:3) = force(1)%gth_ppnl(1:3, 1) - IF (use_virial) e_dummy = third_tr(virial%pv_virial) - END IF - CALL build_core_ppnl(scrm_kp, rho_ao_kp, force, & - virial, .TRUE., use_virial, 1, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, dft_control%qs_control%eps_ppnl, & - nimages=1, basis_type="ORB") - IF (debug_forces) THEN - deb(1:3) = force(1)%gth_ppnl(1:3, 1) - deb(1:3) - CALL para_env%sum(deb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl ", deb - IF (use_virial) THEN - e_dummy = third_tr(virial%pv_virial) - e_dummy - CALL para_env%sum(e_dummy) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hppnl ", e_dummy - END IF - END IF - END IF + CALL core_matrices(qs_env, scrm_kp, rho_ao_kp, .TRUE., 1, & + debug_forces=debug_forces, debug_stress=debug_forces) ! Get the different components of the KS potential NULLIFY (vxc_rspace, vtau_rspace, vadmm_rspace) diff --git a/src/qs_core_hamiltonian.F b/src/qs_core_hamiltonian.F index 02929b2a30..624de9b0ed 100644 --- a/src/qs_core_hamiltonian.F +++ b/src/qs_core_hamiltonian.F @@ -75,17 +75,17 @@ !> \author Matthias Krack (14.09.2000,21.03.02) ! ************************************************************************************************** MODULE qs_core_hamiltonian - USE atomic_kind_types, ONLY: atomic_kind_type,& - get_atomic_kind_set - USE core_ae, ONLY: build_core_ae - USE core_ppl, ONLY: build_core_ppl - USE core_ppnl, ONLY: build_core_ppnl + USE atomic_kind_types, ONLY: atomic_kind_type USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_control_types, ONLY: dft_control_type - USE cp_dbcsr_api, ONLY: & - dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_iterator_blocks_left, & - dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & - dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric + USE cp_dbcsr_api, ONLY: dbcsr_add,& + dbcsr_copy,& + dbcsr_create,& + dbcsr_distribution_type,& + dbcsr_p_type,& + dbcsr_set,& + dbcsr_type,& + dbcsr_type_antisymmetric USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,& dbcsr_deallocate_matrix_set @@ -98,28 +98,22 @@ MODULE qs_core_hamiltonian cp_print_key_should_output,& cp_print_key_unit_nr USE input_constants, ONLY: do_admm_purify_none,& - do_ppl_analytic,& - kg_tnadd_atomic,& - rel_none,& - rel_trans_atom + kg_tnadd_atomic USE input_section_types, ONLY: section_vals_val_get USE kg_environment_types, ONLY: kg_environment_type USE kg_tnadd_mat, ONLY: build_tnadd_mat USE kinds, ONLY: default_string_length,& dp - USE kpoint_types, ONLY: get_kpoint_info,& - kpoint_type - USE lri_environment_types, ONLY: lri_environment_type USE message_passing, ONLY: mp_para_env_type USE particle_types, ONLY: particle_type USE qs_condnum, ONLY: overlap_condnum + USE qs_core_matrices, ONLY: core_matrices,& + kinetic_energy_matrix USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type,& set_qs_env USE qs_force_types, ONLY: qs_force_type - USE qs_kind_types, ONLY: get_qs_kind,& - qs_kind_type - USE qs_kinetic, ONLY: build_kinetic_matrix + USE qs_kind_types, ONLY: qs_kind_type USE qs_ks_types, ONLY: get_ks_env,& qs_ks_env_type,& set_ks_env @@ -140,7 +134,7 @@ MODULE qs_core_hamiltonian CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_hamiltonian' - PUBLIC :: build_core_hamiltonian_matrix, core_matrices + PUBLIC :: build_core_hamiltonian_matrix PUBLIC :: dump_info_core_hamiltonian, qs_matrix_h_allocate_imag_from_real CONTAINS @@ -160,9 +154,8 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_hamiltonian_matrix' - INTEGER :: handle, ic, img, iw, nder, nders, & - nimages, nkind - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index + INTEGER :: handle, img, iw, nder, nders, nimages, & + nkind LOGICAL :: h_is_complex, norml1, norml2, ofdft, & use_arnoldi, use_virial REAL(KIND=dp) :: eps_filter, eps_fit @@ -175,7 +168,6 @@ CONTAINS matrix_w TYPE(dft_control_type), POINTER :: dft_control TYPE(kg_environment_type), POINTER :: kg_env - TYPE(kpoint_type), POINTER :: kpoints TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_orb, sap_oce TYPE(oce_matrix_type), POINTER :: oce @@ -273,21 +265,19 @@ CONTAINS alpha_scalar=1.0_dp, beta_scalar=1.0_dp) END DO END IF - ! S matrix - CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, & - matrix_name="OVERLAP MATRIX", & - basis_type_a="ORB", & - basis_type_b="ORB", & - sab_nl=sab_orb, calculate_forces=.TRUE., & - matrixkp_p=matrix_w) - ! T matrix - IF (.NOT. ofdft) & - CALL build_kinetic_matrix(ks_env, matrixkp_t=matrix_t, & - matrix_name="KINETIC ENERGY MATRIX", & - basis_type="ORB", & - sab_nl=sab_orb, calculate_forces=.TRUE., & - matrixkp_p=matrix_p, & - eps_filter=eps_filter) + ELSE + NULLIFY (matrix_p, matrix_w) + END IF + + ! S matrix + CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, & + matrix_name="OVERLAP MATRIX", & + basis_type_a="ORB", & + basis_type_b="ORB", & + sab_nl=sab_orb, calculate_forces=calculate_forces, & + matrixkp_p=matrix_w) + + IF (calculate_forces) THEN ! *** If LSD, then recover alpha density and beta density *** ! *** from the total density (1) and the spin density (2) *** ! *** The W matrix is neglected, since it will be destroyed *** @@ -300,22 +290,17 @@ CONTAINS alpha_scalar=-1.0_dp, beta_scalar=1.0_dp) END DO END IF - ELSE - ! S matrix - CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, & - matrix_name="OVERLAP MATRIX", & - basis_type_a="ORB", & - basis_type_b="ORB", & - sab_nl=sab_orb) - ! T matrix - IF (.NOT. ofdft) & - CALL build_kinetic_matrix(ks_env, matrixkp_t=matrix_t, & - matrix_name="KINETIC ENERGY MATRIX", & - basis_type="ORB", & - sab_nl=sab_orb, & - eps_filter=eps_filter) END IF + ! T matrix + CALL kinetic_energy_matrix(qs_env, matrixkp_t=matrix_t, & + matrix_p=matrix_p, & + matrix_name="KINETIC ENERGY MATRIX", & + basis_type="ORB", & + sab_orb=sab_orb, & + calculate_forces=calculate_forces, & + eps_filter=eps_filter) + ! (Re-)allocate H matrix based on overlap matrix CALL get_ks_env(ks_env, complex_ks=h_is_complex) CALL qs_matrix_h_allocate(qs_env, matrix_s(1, 1)%matrix, is_complex=h_is_complex) @@ -334,24 +319,6 @@ CONTAINS CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, & particle_set=particle_set) - IF (.NOT. ofdft) THEN - ! relativistic atomic correction to kinetic energy - IF (qs_env%rel_control%rel_method /= rel_none) THEN - IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN - IF (nimages == 1) THEN - ic = 1 - ELSE - CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) - CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) - ic = cell_to_index(0, 0, 0) - END IF - CALL build_atomic_relmat(matrix_h(1, 1)%matrix, atomic_kind_set, qs_kind_set) - ELSE - CPABORT("Relativistic corrections of this type are currently not implemented") - END IF - END IF - END IF - ! *** core and pseudopotentials CALL core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder) @@ -405,161 +372,6 @@ CONTAINS END SUBROUTINE build_core_hamiltonian_matrix -! ************************************************************************************************** -!> \brief ... -!> \param qs_env ... -!> \param matrix_h ... -!> \param matrix_p ... -!> \param calculate_forces ... -!> \param nder ... -!> \param atcore ... -! ************************************************************************************************** - SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, atcore) - - TYPE(qs_environment_type), POINTER :: qs_env - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p - LOGICAL, INTENT(IN) :: calculate_forces - INTEGER, INTENT(IN) :: nder - REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atcore - - INTEGER :: natom, nimages - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index - LOGICAL :: all_present, my_gt_nl, ppl_present, & - ppnl_present, use_virial - REAL(KIND=dp) :: eps_ppnl - TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set - TYPE(dft_control_type), POINTER :: dft_control - TYPE(kpoint_type), POINTER :: kpoints - TYPE(lri_environment_type), POINTER :: lri_env - TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl - TYPE(particle_type), DIMENSION(:), POINTER :: particle_set - TYPE(qs_force_type), DIMENSION(:), POINTER :: force - TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set - TYPE(qs_ks_env_type), POINTER :: ks_env - TYPE(virial_type), POINTER :: virial - - NULLIFY (dft_control) - CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control, natom=natom) - nimages = dft_control%nimages - IF (PRESENT(atcore)) THEN - CPASSERT(SIZE(atcore) >= natom) - END IF - - ! check whether a gauge transformed version of the non-local potential part has to be used - my_gt_nl = .FALSE. - IF (qs_env%run_rtp) THEN - CPASSERT(ASSOCIATED(dft_control%rtp_control)) - IF (dft_control%rtp_control%velocity_gauge) THEN - my_gt_nl = dft_control%rtp_control%nl_gauge_transform - END IF - END IF - - ! prepare for k-points - NULLIFY (cell_to_index) - IF (nimages > 1) THEN - CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) - CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) - dft_control%qs_control%do_ppl_method = do_ppl_analytic - END IF - ! force analytic ppl calculation for GAPW methods - IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN - dft_control%qs_control%do_ppl_method = do_ppl_analytic - END IF - - ! force - NULLIFY (force) - IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force) - ! virial - CALL get_qs_env(qs_env=qs_env, virial=virial) - use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) - - NULLIFY (qs_kind_set, atomic_kind_set, particle_set) - CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, & - particle_set=particle_set) - - NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl) - CALL get_qs_env(qs_env=qs_env, & - sab_orb=sab_orb, & - sac_ae=sac_ae, & - sac_ppl=sac_ppl, & - sap_ppnl=sap_ppnl) - - ! *** compute the nuclear attraction contribution to the core hamiltonian *** - all_present = ASSOCIATED(sac_ae) - IF (all_present) THEN - CALL build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, & - nimages, cell_to_index, "ORB", atcore=atcore) - END IF - - ! *** compute the ppl contribution to the core hamiltonian *** - ppl_present = ASSOCIATED(sac_ppl) - IF (ppl_present) THEN - IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN - IF (dft_control%qs_control%lrigpw) THEN - CALL get_qs_env(qs_env, lri_env=lri_env) - IF (lri_env%ppl_ri) THEN - IF (lri_env%exact_1c_terms) THEN - CPABORT("not implemented") - END IF - ELSE - CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & - nimages, cell_to_index, "ORB", atcore=atcore) - END IF - ELSE - CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & - nimages, cell_to_index, "ORB", atcore=atcore) - END IF - END IF - END IF - - ! *** compute the ppnl contribution to the core hamiltonian *** - eps_ppnl = dft_control%qs_control%eps_ppnl - ppnl_present = ASSOCIATED(sap_ppnl) - IF (ppnl_present) THEN - IF (.NOT. my_gt_nl) THEN - CALL build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, & - nimages, cell_to_index, "ORB", atcore=atcore) - END IF - END IF - - END SUBROUTINE core_matrices - -! ************************************************************************************************** -!> \brief Adds atomic blocks of relativistic correction for the kinetic energy -!> \param matrix_h ... -!> \param atomic_kind_set ... -!> \param qs_kind_set ... -! ************************************************************************************************** - SUBROUTINE build_atomic_relmat(matrix_h, atomic_kind_set, qs_kind_set) - TYPE(dbcsr_type), POINTER :: matrix_h - TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set - TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set - - INTEGER :: iatom, ikind, jatom - INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of - REAL(KIND=dp), DIMENSION(:, :), POINTER :: hblock, reltmat - TYPE(dbcsr_iterator_type) :: iter - - CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of) - - CALL dbcsr_iterator_start(iter, matrix_h) - DO WHILE (dbcsr_iterator_blocks_left(iter)) - CALL dbcsr_iterator_next_block(iter, iatom, jatom, hblock) - IF (iatom == jatom) THEN - ikind = kind_of(iatom) - CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat) - IF (ASSOCIATED(reltmat)) hblock = hblock + reltmat - END IF - END DO - CALL dbcsr_iterator_stop(iter) - - END SUBROUTINE build_atomic_relmat - ! ************************************************************************************************** !> \brief Possibly prints matrices after the construction of the Core !> Hamiltonian Matrix diff --git a/src/qs_core_matrices.F b/src/qs_core_matrices.F new file mode 100644 index 0000000000..6b18d8b09c --- /dev/null +++ b/src/qs_core_matrices.F @@ -0,0 +1,577 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2025 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Calculation of the core Hamiltonian integral matrix over +!> Cartesian Gaussian-type functions. +!> +!> Nuclear potential energy: +!> +!> a) Allelectron calculation: +!> +!> erfc(r) +!> = -Z* +!> r +!> +!> 1 - erf(r) +!> = -Z* +!> r +!> +!> 1 erf(r) +!> = -Z*( - ) +!> r r +!> +!> 1 +!> = -Z*( - N*) +!> r +!> +!> -Z +!> = + Z*N* +!> r +!> \_______/ \_____/ +!> | | +!> nuclear coulomb +!> +!> b) Pseudopotential calculation (Goedecker, Teter and Hutter; GTH): +!> +!> = +!> +!> = + +!> +!> = (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 + +!> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b> +!> +!> = *h(i,j)* +!> \par Literature +!> S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 1703 (1996) +!> C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B 58, 3641 (1998) +!> M. Krack and M. Parrinello, Phys. Chem. Chem. Phys. 2, 2105 (2000) +!> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986) +!> \par History +!> - Joost VandeVondele (April 2003) : added LSD forces +!> - Non-redundant calculation of the non-local part of the GTH PP +!> (22.05.2003,MK) +!> - New parallelization scheme (27.06.2003,MK) +!> - OpenMP version (07.12.2003,JGH) +!> - Binary search loop for VPPNL operators (09.01.2004,JGH,MK) +!> - Refactoring of pseudopotential and nuclear attraction integrals (25.02.2009,JGH) +!> - General refactoring (01.10.2010,JGH) +!> - Refactoring related to the new kinetic energy and overlap routines (07.2014,JGH) +!> - k-point functionality (07.2015,JGH) +!> - Refactor for PP functionality (08.2025,JGH) +!> \author Matthias Krack (14.09.2000,21.03.02) +! ************************************************************************************************** +MODULE qs_core_matrices + USE atomic_kind_types, ONLY: atomic_kind_type,& + get_atomic_kind_set + USE cell_types, ONLY: cell_type + USE core_ae, ONLY: build_core_ae + USE core_ppl, ONLY: build_core_ppl + USE core_ppnl, ONLY: build_core_ppnl + USE cp_control_types, ONLY: dft_control_type + USE cp_dbcsr_api, ONLY: dbcsr_add,& + dbcsr_iterator_blocks_left,& + dbcsr_iterator_next_block,& + dbcsr_iterator_start,& + dbcsr_iterator_stop,& + dbcsr_iterator_type,& + dbcsr_p_type,& + dbcsr_type + USE cp_log_handling, ONLY: cp_get_default_logger,& + cp_logger_get_default_unit_nr,& + cp_logger_type + USE ec_env_types, ONLY: energy_correction_type + USE input_constants, ONLY: do_ppl_analytic,& + rel_none,& + rel_trans_atom + USE kinds, ONLY: default_string_length,& + dp + USE kpoint_types, ONLY: get_kpoint_info,& + kpoint_type + USE lri_environment_types, ONLY: lri_environment_type + USE mathlib, ONLY: det_3x3 + USE message_passing, ONLY: mp_para_env_type + USE particle_types, ONLY: particle_type + USE physcon, ONLY: pascal + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type + USE qs_force_types, ONLY: qs_force_type + USE qs_kind_types, ONLY: get_qs_kind,& + qs_kind_type + USE qs_kinetic, ONLY: build_kinetic_matrix + USE qs_ks_types, ONLY: get_ks_env,& + qs_ks_env_type + USE qs_linres_types, ONLY: dcdr_env_type + USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type + USE virial_methods, ONLY: one_third_sum_diag + USE virial_types, ONLY: virial_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_matrices' + + PUBLIC :: core_matrices, kinetic_energy_matrix + +CONTAINS + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param matrix_h ... +!> \param matrix_p ... +!> \param calculate_forces ... +!> \param nder ... +!> \param ec_env ... +!> \param dcdr_env ... +!> \param ec_env_matrices ... +!> \param basis_type ... +!> \param debug_forces ... +!> \param debug_stress ... +!> \param atcore ... +! ************************************************************************************************** + SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, ec_env, dcdr_env, & + ec_env_matrices, basis_type, debug_forces, debug_stress, atcore) + + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p + LOGICAL, INTENT(IN) :: calculate_forces + INTEGER, INTENT(IN) :: nder + TYPE(energy_correction_type), OPTIONAL, POINTER :: ec_env + TYPE(dcdr_env_type), OPTIONAL :: dcdr_env + LOGICAL, INTENT(IN), OPTIONAL :: ec_env_matrices + CHARACTER(LEN=*), OPTIONAL :: basis_type + LOGICAL, INTENT(IN), OPTIONAL :: debug_forces, debug_stress + REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atcore + + CHARACTER(LEN=default_string_length) :: my_basis_type + INTEGER :: iounit, natom, nimages + INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index + LOGICAL :: all_present, debfor, debstr, my_gt_nl, & + ppl_present, ppnl_present, use_virial + REAL(KIND=dp) :: eps_ppnl, fconv + REAL(KIND=dp), DIMENSION(3) :: fodeb + REAL(KIND=dp), DIMENSION(3, 3) :: stdeb + REAL(KIND=dp), DIMENSION(:, :), POINTER :: deltaR + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cell_type), POINTER :: cell + TYPE(cp_logger_type), POINTER :: logger + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_hloc, matrix_ploc + TYPE(dft_control_type), POINTER :: dft_control + TYPE(kpoint_type), POINTER :: kpoints + TYPE(lri_environment_type), POINTER :: lri_env + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(virial_type), POINTER :: virial + + logger => cp_get_default_logger() + IF (logger%para_env%is_source()) THEN + iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + iounit = -1 + END IF + + NULLIFY (dft_control) + CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env) + + ! basis type + IF (PRESENT(basis_type)) THEN + my_basis_type = basis_type + ELSE + my_basis_type = "ORB" + END IF + + IF (PRESENT(dcdr_env) .AND. PRESENT(ec_env)) THEN + CPABORT("core_matrices: conflicting options") + END IF + + ! check size of atcore array + IF (PRESENT(atcore)) THEN + CALL get_qs_env(qs_env=qs_env, natom=natom) + CPASSERT(SIZE(atcore) >= natom) + END IF + + ! check whether a gauge transformed version of the non-local potential part has to be used + my_gt_nl = .FALSE. + IF (qs_env%run_rtp) THEN + CPASSERT(ASSOCIATED(dft_control%rtp_control)) + IF (dft_control%rtp_control%velocity_gauge) THEN + my_gt_nl = dft_control%rtp_control%nl_gauge_transform + END IF + END IF + + ! prepare for k-points + nimages = dft_control%nimages + NULLIFY (cell_to_index) + IF (nimages > 1) THEN + CALL get_qs_env(qs_env=qs_env, ks_env=ks_env) + CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) + CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) + dft_control%qs_control%do_ppl_method = do_ppl_analytic + END IF + + ! Possible dc/dr terms + IF (PRESENT(dcdr_env)) THEN + deltaR => dcdr_env%delta_basis_function + matrix_hloc(1:3, 1:1) => dcdr_env%matrix_hc(1:3) + ELSE + NULLIFY (deltaR) + matrix_hloc => matrix_h + END IF + matrix_ploc => matrix_p + + ! force analytic ppl calculation for GAPW methods + IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN + dft_control%qs_control%do_ppl_method = do_ppl_analytic + END IF + + ! force + NULLIFY (force) + IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force) + ! virial + CALL get_qs_env(qs_env=qs_env, virial=virial) + use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) + + !debug + debfor = .FALSE. + IF (PRESENT(debug_forces)) debfor = debug_forces + debfor = debfor .AND. calculate_forces + debstr = .FALSE. + IF (PRESENT(debug_stress)) debstr = debug_stress + debstr = debstr .AND. use_virial + IF (debstr) THEN + CALL get_qs_env(qs_env=qs_env, cell=cell) + fconv = 1.0E-9_dp*pascal/cell%deth + END IF + + NULLIFY (qs_kind_set, atomic_kind_set, particle_set) + CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, & + particle_set=particle_set) + + NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl) + IF (PRESENT(ec_env)) THEN + sab_orb => ec_env%sab_orb + sac_ae => ec_env%sac_ae + sac_ppl => ec_env%sac_ppl + sap_ppnl => ec_env%sap_ppnl + ELSE + CALL get_qs_env(qs_env=qs_env, & + sab_orb=sab_orb, & + sac_ae=sac_ae, & + sac_ppl=sac_ppl, & + sap_ppnl=sap_ppnl) + END IF + + IF (PRESENT(ec_env) .AND. PRESENT(ec_env_matrices)) THEN + IF (ec_env_matrices) THEN + matrix_hloc => ec_env%matrix_h + matrix_ploc => ec_env%matrix_p + END IF + END IF + + ! *** compute the nuclear attraction contribution to the core hamiltonian *** + all_present = ASSOCIATED(sac_ae) + IF (all_present) THEN + IF (PRESENT(dcdr_env)) THEN + CPABORT("ECP/AE functionality for dcdr missing") + END IF + IF (debfor) fodeb(1:3) = force(1)%all_potential(1:3, 1) + IF (debstr) stdeb = virial%pv_ppl + CALL build_core_ae(matrix_hloc, matrix_ploc, force, virial, calculate_forces, use_virial, nder, & + qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, & + nimages, cell_to_index, my_basis_type, atcore=atcore) + IF (debfor) THEN + fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae ", fodeb + END IF + IF (debstr) THEN + stdeb = fconv*(virial%pv_ppl - stdeb) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| P*dHae ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + END IF + + ! *** compute the ppl contribution to the core hamiltonian *** + ppl_present = ASSOCIATED(sac_ppl) + IF (ppl_present) THEN + IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN + CALL get_qs_env(qs_env, lri_env=lri_env) + IF (dft_control%qs_control%lrigpw .AND. lri_env%ppl_ri) THEN + IF (lri_env%exact_1c_terms) THEN + CPABORT("not implemented") + END IF + ELSE + IF (debfor) fodeb(1:3) = force(1)%gth_ppl(1:3, 1) + IF (debstr) stdeb = virial%pv_ppl + CALL build_core_ppl(matrix_hloc, matrix_ploc, force, & + virial, calculate_forces, use_virial, nder, & + qs_kind_set, atomic_kind_set, particle_set, & + sab_orb, sac_ppl, nimages, cell_to_index, my_basis_type, & + deltaR=deltaR, atcore=atcore) + IF (debfor) THEN + fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl ", fodeb + END IF + IF (debstr) THEN + stdeb = fconv*(virial%pv_ppl - stdeb) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| P*dHppl ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + END IF + END IF + END IF + + ! *** compute the ppnl contribution to the core hamiltonian *** + eps_ppnl = dft_control%qs_control%eps_ppnl + ppnl_present = ASSOCIATED(sap_ppnl) + IF (ppnl_present) THEN + IF (PRESENT(dcdr_env)) THEN + matrix_hloc(1:3, 1:1) => dcdr_env%matrix_ppnl_1(1:3) + END IF + IF (.NOT. my_gt_nl) THEN + IF (debfor) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) + IF (debstr) stdeb = virial%pv_ppnl + CALL build_core_ppnl(matrix_hloc, matrix_ploc, force, & + virial, calculate_forces, use_virial, nder, & + qs_kind_set, atomic_kind_set, particle_set, & + sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, my_basis_type, & + deltaR=deltaR, atcore=atcore) + IF (debfor) THEN + fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl ", fodeb + END IF + IF (debstr) THEN + stdeb = fconv*(virial%pv_ppnl - stdeb) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| P*dHppnl ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + END IF + END IF + + END SUBROUTINE core_matrices + +! ************************************************************************************************** +!> \brief Calculate kinetic energy matrix and possible relativistic correction +!> \param qs_env ... +!> \param matrixkp_t ... +!> \param matrix_t ... +!> \param matrix_p ... +!> \param matrix_name ... +!> \param calculate_forces ... +!> \param nderivative ... +!> \param sab_orb ... +!> \param eps_filter ... +!> \param basis_type ... +!> \param debug_forces ... +!> \param debug_stress ... +!> \author Creation (21.08.2025,JGH) +! ************************************************************************************************** + SUBROUTINE kinetic_energy_matrix(qs_env, matrixkp_t, matrix_t, matrix_p, & + matrix_name, calculate_forces, nderivative, & + sab_orb, eps_filter, basis_type, debug_forces, debug_stress) + + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, & + POINTER :: matrixkp_t + TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & + POINTER :: matrix_t + TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, & + POINTER :: matrix_p + CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name + LOGICAL, INTENT(IN) :: calculate_forces + INTEGER, INTENT(IN), OPTIONAL :: nderivative + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + OPTIONAL, POINTER :: sab_orb + REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter + CHARACTER(LEN=*), OPTIONAL :: basis_type + LOGICAL, INTENT(IN), OPTIONAL :: debug_forces, debug_stress + + CHARACTER(LEN=default_string_length) :: my_basis_type + INTEGER :: ic, img, iounit, nimages + INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index + LOGICAL :: debfor, debstr, use_virial + REAL(KIND=dp) :: fconv + REAL(KIND=dp), DIMENSION(3) :: fodeb + REAL(KIND=dp), DIMENSION(3, 3) :: stdeb + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cell_type), POINTER :: cell + TYPE(cp_logger_type), POINTER :: logger + TYPE(dft_control_type), POINTER :: dft_control + TYPE(kpoint_type), POINTER :: kpoints + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_nl + TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(virial_type), POINTER :: virial + + logger => cp_get_default_logger() + IF (logger%para_env%is_source()) THEN + iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + iounit = -1 + END IF + + CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env) + ! is this a orbital-free method calculation + IF (dft_control%qs_control%ofgpw) RETURN + + ! Matrix images (kp) + nimages = dft_control%nimages + + ! basis type + IF (PRESENT(basis_type)) THEN + my_basis_type = basis_type + ELSE + my_basis_type = "ORB" + END IF + + debfor = .FALSE. + IF (PRESENT(debug_forces)) debfor = debug_forces + debfor = debfor .AND. calculate_forces + + CALL get_qs_env(qs_env=qs_env, virial=virial) + use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) + + debstr = .FALSE. + IF (PRESENT(debug_stress)) debstr = debug_stress + debstr = debstr .AND. use_virial + IF (debstr) THEN + CALL get_qs_env(qs_env=qs_env, cell=cell) + fconv = 1.0E-9_dp*pascal/cell%deth + END IF + + NULLIFY (ks_env) + CALL get_qs_env(qs_env=qs_env, ks_env=ks_env) + NULLIFY (sab_nl) + IF (PRESENT(sab_orb)) THEN + sab_nl => sab_orb + ELSE + CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl) + END IF + + IF (calculate_forces) THEN + IF (SIZE(matrix_p, 1) == 2) THEN + DO img = 1, nimages + CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & + alpha_scalar=1.0_dp, beta_scalar=1.0_dp) + END DO + END IF + END IF + + IF (debfor) THEN + CALL get_qs_env(qs_env=qs_env, force=force) + fodeb(1:3) = force(1)%kinetic(1:3, 1) + END IF + IF (debstr) THEN + stdeb = virial%pv_ekinetic + END IF + + ! T matrix + CALL build_kinetic_matrix(ks_env, matrixkp_t=matrixkp_t, matrix_t=matrix_t, & + matrix_name=matrix_name, & + basis_type=my_basis_type, & + sab_nl=sab_nl, & + calculate_forces=calculate_forces, & + matrixkp_p=matrix_p, & + nderivative=nderivative, & + eps_filter=eps_filter) + + IF (debfor) THEN + fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT ", fodeb + END IF + IF (debstr) THEN + stdeb = fconv*(virial%pv_ekinetic - stdeb) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| P*dT', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + + IF (calculate_forces) THEN + IF (SIZE(matrix_p, 1) == 2) THEN + DO img = 1, nimages + CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & + alpha_scalar=1.0_dp, beta_scalar=-1.0_dp) + END DO + END IF + END IF + + ! relativistic atomic correction to kinetic energy + IF (qs_env%rel_control%rel_method /= rel_none) THEN + IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN + CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set) + IF (nimages == 1) THEN + ic = 1 + ELSE + CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) + CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) + ic = cell_to_index(0, 0, 0) + END IF + IF (my_basis_type /= "ORB") THEN + CPABORT("Basis mismatch for relativistic correction") + END IF + IF (PRESENT(matrixkp_t)) THEN + CALL build_atomic_relmat(matrixkp_t(1, ic)%matrix, atomic_kind_set, qs_kind_set) + ELSEIF (PRESENT(matrix_t)) THEN + CALL build_atomic_relmat(matrix_t(1)%matrix, atomic_kind_set, qs_kind_set) + END IF + ELSE + CPABORT("Relativistic corrections of this type are currently not implemented") + END IF + END IF + + END SUBROUTINE kinetic_energy_matrix + +! ************************************************************************************************** +!> \brief Adds atomic blocks of relativistic correction for the kinetic energy +!> \param matrix_t ... +!> \param atomic_kind_set ... +!> \param qs_kind_set ... +! ************************************************************************************************** + SUBROUTINE build_atomic_relmat(matrix_t, atomic_kind_set, qs_kind_set) + TYPE(dbcsr_type), POINTER :: matrix_t + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + + INTEGER :: iatom, ikind, jatom + INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of + REAL(KIND=dp), DIMENSION(:, :), POINTER :: reltmat, tblock + TYPE(dbcsr_iterator_type) :: iter + + CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of) + + CALL dbcsr_iterator_start(iter, matrix_t) + DO WHILE (dbcsr_iterator_blocks_left(iter)) + CALL dbcsr_iterator_next_block(iter, iatom, jatom, tblock) + IF (iatom == jatom) THEN + ikind = kind_of(iatom) + CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat) + IF (ASSOCIATED(reltmat)) tblock = tblock + reltmat + END IF + END DO + CALL dbcsr_iterator_stop(iter) + + END SUBROUTINE build_atomic_relmat + +END MODULE qs_core_matrices diff --git a/src/qs_dcdr_ao.F b/src/qs_dcdr_ao.F index f2b6963bda..f577c1406d 100644 --- a/src/qs_dcdr_ao.F +++ b/src/qs_dcdr_ao.F @@ -12,12 +12,8 @@ MODULE qs_dcdr_ao - USE atomic_kind_types, ONLY: atomic_kind_type USE basis_set_types, ONLY: gto_basis_set_p_type,& gto_basis_set_type - USE core_ae, ONLY: build_core_ae - USE core_ppl, ONLY: build_core_ppl - USE core_ppnl, ONLY: build_core_ppnl USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_api, ONLY: dbcsr_copy,& dbcsr_get_block_p,& @@ -31,14 +27,12 @@ MODULE qs_dcdr_ao cp_fm_type USE cp_log_handling, ONLY: cp_get_default_logger,& cp_logger_type - USE input_constants, ONLY: do_ppl_analytic USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type USE kinds, ONLY: default_string_length,& dp USE orbital_pointers, ONLY: ncoset USE parallel_gemm_api, ONLY: parallel_gemm - USE particle_types, ONLY: particle_type USE pw_env_types, ONLY: pw_env_get,& pw_env_type USE pw_methods, ONLY: pw_axpy,& @@ -54,10 +48,10 @@ MODULE qs_dcdr_ao pw_r3d_rs_type USE qs_collocate_density, ONLY: calculate_drho_core,& calculate_drho_elec_dR + USE qs_core_matrices, ONLY: core_matrices USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type - USE qs_force_types, ONLY: qs_force_type USE qs_integral_utils, ONLY: basis_set_list_setup,& get_memory_usage USE qs_integrate_potential, ONLY: integrate_v_dbasis,& @@ -80,7 +74,6 @@ MODULE qs_dcdr_ao qs_rho_release,& qs_rho_type USE qs_vxc, ONLY: qs_vxc_create - USE virial_types, ONLY: virial_type USE xc, ONLY: xc_calc_2nd_deriv,& xc_prep_2nd_deriv USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& @@ -403,95 +396,30 @@ CONTAINS CHARACTER(LEN=default_string_length) :: my_basis_type INTEGER :: handle, nder - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index - LOGICAL :: calculate_forces, failure, ppl_present, & - ppnl_present, use_virial - REAL(KIND=dp) :: eps_ppnl - REAL(KIND=dp), DIMENSION(:, :), POINTER :: deltaR - TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + LOGICAL :: calculate_forces TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_hc_pass, matrix_p_pass, & - matrix_ppnl_1_pass - TYPE(dft_control_type), POINTER :: dft_control - TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl - TYPE(particle_type), DIMENSION(:), POINTER :: particle_set - TYPE(qs_force_type), DIMENSION(:), POINTER :: force - TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p_pass TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(qs_rho_type), POINTER :: rho - TYPE(virial_type), POINTER :: virial CALL timeset(routineN, handle) - failure = .FALSE. - - NULLIFY (atomic_kind_set, qs_kind_set, ks_env, dft_control, particle_set, & - sab_orb, sac_ae, sac_ppl, sap_ppnl, virial, rho, rho_ao) - - CALL get_qs_env(qs_env=qs_env, & - atomic_kind_set=atomic_kind_set, & - qs_kind_set=qs_kind_set, & - ks_env=ks_env, & - dft_control=dft_control, & - particle_set=particle_set, & - sab_orb=sab_orb, & - sac_ae=sac_ae, & - sac_ppl=sac_ppl, & - sap_ppnl=sap_ppnl, & - virial=virial) + CALL get_qs_env(qs_env=qs_env, ks_env=ks_env) CALL get_ks_env(ks_env=ks_env, rho=rho) CALL qs_rho_get(rho, rho_ao=rho_ao) - deltaR => dcdr_env%delta_basis_function nder = 1 calculate_forces = .FALSE. my_basis_type = "ORB" - ! ECP/AE contribution to the core hamiltonian - IF (ASSOCIATED(sac_ae)) THEN - CPABORT("ECP/AE functionality in qs_dcdr_ao missing") - ! Missing feature: deltaR weighting factors of the derivatives wrt. nuclear positions - matrix_hc_pass(1:3, 1:1) => dcdr_env%matrix_hc(1:3) - matrix_p_pass(1:1, 1:1) => rho_ao(1:1) - CALL build_core_ae(matrix_h=matrix_hc_pass, matrix_p=matrix_p_pass, & - force=force, virial=virial, calculate_forces=calculate_forces, & - use_virial=use_virial, nder=nder, qs_kind_set=qs_kind_set, & - atomic_kind_set=atomic_kind_set, particle_set=particle_set, & - sab_orb=sab_orb, sac_ae=sac_ae, nimages=1, & - basis_type=my_basis_type, cell_to_index=cell_to_index) - END IF - ! *** compute the ppl contribution to the core hamiltonian *** - ppl_present = ASSOCIATED(sac_ppl) - IF (ppl_present) THEN - IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN - matrix_hc_pass(1:3, 1:1) => dcdr_env%matrix_hc(1:3) - matrix_p_pass(1:1, 1:1) => rho_ao(1:1) - CALL build_core_ppl(matrix_h=matrix_hc_pass, matrix_p=matrix_p_pass, & - force=force, virial=virial, calculate_forces=calculate_forces, & - use_virial=use_virial, nder=nder, qs_kind_set=qs_kind_set, & - atomic_kind_set=atomic_kind_set, particle_set=particle_set, & - sab_orb=sab_orb, sac_ppl=sac_ppl, basis_type=my_basis_type, & - nimages=1, cell_to_index=cell_to_index, deltaR=deltaR) - - END IF ! ppl_analytic - END IF ! ppl_present - - ! *** compute the ppnl contribution to the core hamiltonian *** - eps_ppnl = dft_control%qs_control%eps_ppnl - ppnl_present = ASSOCIATED(sap_ppnl) - IF (ppnl_present) THEN - matrix_ppnl_1_pass(1:3, 1:1) => dcdr_env%matrix_ppnl_1(1:3) - CALL build_core_ppnl(matrix_h=matrix_ppnl_1_pass, matrix_p=matrix_p_pass, force=force, virial=virial, & - calculate_forces=calculate_forces, use_virial=use_virial, nder=nder, & - qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, & - particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl, & - eps_ppnl=eps_ppnl, nimages=1, cell_to_index=cell_to_index, & - basis_type=my_basis_type, deltaR=deltaR) - END IF + NULLIFY (matrix_h) + matrix_p_pass(1:1, 1:1) => rho_ao(1:1) + CALL core_matrices(qs_env, matrix_h, matrix_p_pass, calculate_forces, nder, & + dcdr_env=dcdr_env) CALL timestop(handle) + END SUBROUTINE core_dR ! ************************************************************************************************** diff --git a/src/qs_dcdr_utils.F b/src/qs_dcdr_utils.F index beb5910c1e..84c7dba9a6 100644 --- a/src/qs_dcdr_utils.F +++ b/src/qs_dcdr_utils.F @@ -64,9 +64,9 @@ MODULE qs_dcdr_utils USE moments_utils, ONLY: get_reference_point USE parallel_gemm_api, ONLY: parallel_gemm USE particle_types, ONLY: particle_type + USE qs_core_matrices, ONLY: kinetic_energy_matrix USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type - USE qs_kinetic, ONLY: build_kinetic_matrix USE qs_ks_types, ONLY: qs_ks_env_type USE qs_linres_types, ONLY: dcdr_env_type,& linres_control_type @@ -666,6 +666,7 @@ CONTAINS TYPE(cp_fm_type), POINTER :: mo_coeff TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p TYPE(dft_control_type), POINTER :: dft_control TYPE(mo_set_type), DIMENSION(:), POINTER :: mos TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set @@ -723,11 +724,20 @@ CONTAINS sab_nl=sab_orb) NULLIFY (dcdr_env%matrix_t) - CALL build_kinetic_matrix(ks_env, matrix_t=dcdr_env%matrix_t, & - matrix_name="KINETIC ENERGY MATRIX", & - basis_type="ORB", & - sab_nl=sab_orb, nderivative=1, & - eps_filter=dft_control%qs_control%eps_filter_matrix) + NULLIFY (matrix_p) + CALL kinetic_energy_matrix(qs_env, matrix_t=dcdr_env%matrix_t, matrix_p=matrix_p, & + matrix_name="KINETIC ENERGY MATRIX", & + calculate_forces=.FALSE., & + basis_type="ORB", & + sab_orb=sab_orb, & + nderivative=1, & + eps_filter=dft_control%qs_control%eps_filter_matrix) + +! CALL build_kinetic_matrix(ks_env, matrix_t=dcdr_env%matrix_t, & +! matrix_name="KINETIC ENERGY MATRIX", & +! basis_type="ORB", & +! sab_nl=sab_orb, nderivative=1, & +! eps_filter=dft_control%qs_control%eps_filter_matrix) ! Get inputs CALL section_vals_val_get(dcdr_section, "DISTRIBUTED_ORIGIN", l_val=dcdr_env%distributed_origin) diff --git a/src/qs_energy_utils.F b/src/qs_energy_utils.F index 448acc8800..b89c935190 100644 --- a/src/qs_energy_utils.F +++ b/src/qs_energy_utils.F @@ -40,7 +40,7 @@ MODULE qs_energy_utils pw_scale USE pw_pool_types, ONLY: pw_pool_type USE pw_types, ONLY: pw_r3d_rs_type - USE qs_core_hamiltonian, ONLY: core_matrices + USE qs_core_matrices, ONLY: core_matrices USE qs_dispersion_types, ONLY: qs_dispersion_type USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env,& diff --git a/src/qs_kinetic.F b/src/qs_kinetic.F index 5bdb63bb07..c7f2f8113a 100644 --- a/src/qs_kinetic.F +++ b/src/qs_kinetic.F @@ -186,6 +186,7 @@ CONTAINS REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dab TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: k_block + TYPE(dbcsr_type), POINTER :: matp TYPE(dft_control_type), POINTER :: dft_control TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b @@ -260,7 +261,13 @@ CONTAINS IF (dokp) THEN CPASSERT(PRESENT(matrixkp_p)) ELSE - CPASSERT(PRESENT(matrix_p)) + IF (PRESENT(matrix_p)) THEN + matp => matrix_p + ELSEIF (PRESENT(matrixkp_p)) THEN + matp => matrixkp_p(1, 1)%matrix + ELSE + CPABORT("Missing density matrix") + END IF END IF END IF @@ -273,7 +280,7 @@ CONTAINS !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED (do_forces, ldsab, use_cell_mapping, do_symmetric, dokp,& !$OMP sab_nl, ncoset, maxder, nder, ndod, use_virial, matrix_t, matrixkp_t,& -!$OMP matrix_p, basis_set_list, atom_of_kind, cell_to_index, matrixkp_p, locks, natom) & +!$OMP matp, matrix_p, basis_set_list, atom_of_kind, cell_to_index, matrixkp_p, locks, natom) & !$OMP PRIVATE (k_block, kab, qab, pab, ikind, jkind, iatom, jatom, rab, cell, basis_set_a, basis_set_b, f, & !$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, & !$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, tab, & @@ -386,7 +393,9 @@ CONTAINS row=irow, col=icol, block=p_block, found=found) CPASSERT(found) ELSE - CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, & +!deb CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, & +!deb block=p_block, found=found) + CALL dbcsr_get_block_p(matrix=matp, row=irow, col=icol, & block=p_block, found=found) CPASSERT(found) END IF diff --git a/src/response_solver.F b/src/response_solver.F index 04409f2bf9..33af6ec231 100644 --- a/src/response_solver.F +++ b/src/response_solver.F @@ -21,9 +21,6 @@ MODULE response_solver USE atomic_kind_types, ONLY: atomic_kind_type,& get_atomic_kind USE cell_types, ONLY: cell_type - USE core_ae, ONLY: build_core_ae - USE core_ppl, ONLY: build_core_ppl - USE core_ppnl, ONLY: build_core_ppnl USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_api, ONLY: & @@ -98,6 +95,8 @@ MODULE response_solver pw_r3d_rs_type USE qs_2nd_kernel_ao, ONLY: build_dm_response USE qs_collocate_density, ONLY: calculate_rho_elec + USE qs_core_matrices, ONLY: core_matrices,& + kinetic_energy_matrix USE qs_density_matrices, ONLY: calculate_whz_matrix,& calculate_wz_matrix USE qs_energy_types, ONLY: qs_energy_type @@ -112,7 +111,6 @@ MODULE response_solver USE qs_kind_types, ONLY: get_qs_kind,& get_qs_kind_set,& qs_kind_type - USE qs_kinetic, ONLY: build_kinetic_matrix USE qs_ks_atom, ONLY: update_ks_atom USE qs_ks_methods, ONLY: calc_rho_tot_gspace USE qs_ks_types, ONLY: qs_ks_env_type @@ -834,13 +832,12 @@ CONTAINS CHARACTER(LEN=default_string_length) :: basis_type, unitstr INTEGER :: handle, iounit, ispin, mspin, myfun, & n_rep_hf, nao, nao_aux, natom, nder, & - nimages, nocc, nspins - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index + nocc, nspins LOGICAL :: debug_forces, debug_stress, distribute_fock_matrix, do_ex, do_hfx, gapw, gapw_xc, & hfx_treat_lsd_in_core, resp_only, s_mstruct_changed, use_virial REAL(KIND=dp) :: eh1, ehartree, ekin_mol, eps_filter, & - eps_ppnl, exc, exc_aux_fit, fconv, & - focc, hartree_gs, hartree_t + exc, exc_aux_fit, fconv, focc, & + hartree_gs, hartree_t REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot1, ftot2, ftot3 REAL(KIND=dp), DIMENSION(2) :: total_rho_gs, total_rho_t REAL(KIND=dp), DIMENSION(3) :: fodeb @@ -853,7 +850,7 @@ CONTAINS TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ht, matrix_pd, matrix_pza, & matrix_s, mpa, scrm TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, mhd, mhx, mhy, mhz, & - mpd, mpz + mpa2, mpd, mpz, scrm2 TYPE(dbcsr_type), POINTER :: dbwork TYPE(dft_control_type), POINTER :: dft_control TYPE(hartree_local_type), POINTER :: hartree_local_gs, hartree_local_t @@ -865,8 +862,7 @@ CONTAINS TYPE(mo_set_type), DIMENSION(:), POINTER :: mos TYPE(mp_para_env_type), POINTER :: para_env TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_aux_fit, sab_orb, sac_ae, sac_ppl, & - sap_ppnl + POINTER :: sab_aux_fit, sab_orb TYPE(oce_matrix_type), POINTER :: oce TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rho_tot_gspace_gs, rho_tot_gspace_t, & @@ -900,7 +896,6 @@ CONTAINS debug_stress = .FALSE. END IF - logger => cp_get_default_logger() logger => cp_get_default_logger() IF (logger%para_env%is_source()) THEN iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.) @@ -914,7 +909,7 @@ CONTAINS CPASSERT(PRESENT(p_env)) END IF - NULLIFY (ks_env, sab_orb, sac_ae, sac_ppl, sap_ppnl, virial) + NULLIFY (ks_env, sab_orb, virial) CALL get_qs_env(qs_env=qs_env, & cell=cell, & force=force, & @@ -922,9 +917,6 @@ CONTAINS dft_control=dft_control, & para_env=para_env, & sab_orb=sab_orb, & - sac_ae=sac_ae, & - sac_ppl=sac_ppl, & - sap_ppnl=sap_ppnl, & virial=virial) nspins = dft_control%nspins gapw = dft_control%qs_control%gapw @@ -975,39 +967,20 @@ CONTAINS ! ! START OF Tr(P+Z)Hcore ! - IF (nspins == 2) THEN - CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, 1.0_dp) - END IF - nimages = 1 ! Kinetic energy matrix - NULLIFY (scrm) - IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ekinetic - CALL build_kinetic_matrix(ks_env, matrix_t=scrm, & - matrix_name="KINETIC ENERGY MATRIX", & - basis_type="ORB", & - sab_nl=sab_orb, calculate_forces=.TRUE., & - matrix_p=mpa(1)%matrix) - IF (debug_forces) THEN - fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dT ", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ekinetic - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Kinetic energy', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - - IF (nspins == 2) THEN - CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, -1.0_dp) - END IF - CALL dbcsr_deallocate_matrix_set(scrm) + NULLIFY (scrm2) + mpa2(1:nspins, 1:1) => mpa(1:nspins) + CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm2, matrix_p=mpa2, & + matrix_name="KINETIC ENERGY MATRIX", & + basis_type="ORB", & + sab_orb=sab_orb, calculate_forces=.TRUE., & + debug_forces=debug_forces, debug_stress=debug_stress) + CALL dbcsr_deallocate_matrix_set(scrm2) ! Initialize a matrix scrm, later used for scratch purposes CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s) + NULLIFY (scrm) CALL dbcsr_allocate_matrix_set(scrm, nspins) DO ispin = 1, nspins ALLOCATE (scrm(ispin)%matrix) @@ -1019,7 +992,6 @@ CONTAINS CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, & atomic_kind_set=atomic_kind_set) - NULLIFY (cell_to_index) ALLOCATE (matrix_p(nspins, 1), matrix_h(nspins, 1)) DO ispin = 1, nspins matrix_p(ispin, 1)%matrix => mpa(ispin)%matrix @@ -1027,68 +999,10 @@ CONTAINS END DO matrix_h(1, 1)%matrix => scrm(1)%matrix - IF (ASSOCIATED(sac_ae)) THEN - nder = 1 - IF (debug_forces) fodeb(1:3) = force(1)%all_potential(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl - CALL build_core_ae(matrix_h, matrix_p, force, virial, .TRUE., use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, & - nimages, cell_to_index, "ORB") - IF (debug_forces) THEN - fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHae ", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ppl - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pz*dHae ', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - END IF + nder = 1 + CALL core_matrices(qs_env, matrix_h, matrix_p, .TRUE., nder, & + debug_forces=debug_forces, debug_stress=debug_stress) - IF (ASSOCIATED(sac_ppl)) THEN - nder = 1 - IF (debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl - CALL build_core_ppl(matrix_h, matrix_p, force, & - virial, .TRUE., use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, & - sab_orb, sac_ppl, nimages, cell_to_index, "ORB") - IF (debug_forces) THEN - fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHppl ", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ppl - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pz*dHppl ', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - END IF - eps_ppnl = dft_control%qs_control%eps_ppnl - IF (ASSOCIATED(sap_ppnl)) THEN - nder = 1 - IF (debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl - CALL build_core_ppnl(matrix_h, matrix_p, force, & - virial, .TRUE., use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, & - sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, "ORB") - IF (debug_forces) THEN - fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3) - CALL para_env%sum(fodeb) - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHppnl ", fodeb - END IF - IF (debug_stress .AND. use_virial) THEN - stdeb = fconv*(virial%pv_ppnl - stdeb) - CALL para_env%sum(stdeb) - IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & - 'STRESS| Pz*dHppnl ', one_third_sum_diag(stdeb), det_3x3(stdeb) - END IF - - END IF ! Kim-Gordon subsystem DFT ! Atomic potential for nonadditive kinetic energy contribution IF (dft_control%qs_control%do_kg) THEN diff --git a/src/rpa_im_time_force_methods.F b/src/rpa_im_time_force_methods.F index c583c831bc..08b9ca1766 100644 --- a/src/rpa_im_time_force_methods.F +++ b/src/rpa_im_time_force_methods.F @@ -22,9 +22,6 @@ MODULE rpa_im_time_force_methods cite_reference USE cell_types, ONLY: cell_type,& pbc - USE core_ae, ONLY: build_core_ae - USE core_ppl, ONLY: build_core_ppl - USE core_ppnl, ONLY: build_core_ppnl USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_api, ONLY: & @@ -118,6 +115,8 @@ MODULE rpa_im_time_force_methods pw_r3d_rs_type USE qs_collocate_density, ONLY: calculate_rho_elec,& collocate_function + USE qs_core_matrices, ONLY: core_matrices,& + kinetic_energy_matrix USE qs_density_matrices, ONLY: calculate_whz_matrix USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type,& @@ -129,7 +128,6 @@ MODULE rpa_im_time_force_methods integrate_v_rspace USE qs_interactions, ONLY: init_interaction_radii_orb_basis USE qs_kind_types, ONLY: qs_kind_type - USE qs_kinetic, ONLY: build_kinetic_matrix USE qs_ks_methods, ONLY: calc_rho_tot_gspace USE qs_ks_reference, ONLY: ks_ref_potential USE qs_ks_types, ONLY: set_ks_env @@ -3136,18 +3134,17 @@ CONTAINS INTEGER :: handle, i, idens, ispin, n_rep_hf, nao, & nao_aux, nder, nimages, nocc, nspins - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index LOGICAL :: do_exx, do_hfx, do_tau, do_tau_admm, & use_virial - REAL(dp) :: dummy_real1, dummy_real2, ehartree, & - eps_ppnl, exc, focc + REAL(dp) :: dummy_real1, dummy_real2, ehartree, exc, & + focc REAL(dp), DIMENSION(3, 3) :: h_stress, pv_loc TYPE(admm_type), POINTER :: admm_env TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: current_density, current_density_admm, & current_mat_h, matrix_p_mp2, matrix_p_mp2_admm, matrix_s, matrix_s_aux_fit, matrix_w, & rho_ao, rho_ao_aux, scrm, scrm_admm - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: dbcsr_work_h, dbcsr_work_p + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: dbcsr_work_h, dbcsr_work_p, mpa2 TYPE(dbcsr_type) :: dbcsr_work TYPE(dft_control_type), POINTER :: dft_control TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data @@ -3172,12 +3169,14 @@ CONTAINS TYPE(task_list_type), POINTER :: task_list_aux_fit TYPE(virial_type), POINTER :: virial - NULLIFY (scrm, rho, dft_control, matrix_p_mp2, matrix_s, matrix_p_mp2_admm, admm_env, sab_orb, & - cell_to_index, dbcsr_work_p, dbcsr_work_h, sac_ae, sac_ppl, sap_ppnl, force, virial, & - qs_kind_set, atomic_kind_set, particle_set, pw_env, poisson_env, auxbas_pw_pool, & - task_list_aux_fit, matrix_s_aux_fit, scrm_admm, rho_aux_fit, rho_ao_aux, x_data, & - hfx_section, xc_section, para_env, rhoz_g, rhoz_r, tauz_r, v_xc, v_xc_tau, & - vxc_rspace, vtau_rspace, vadmm_rspace, rho_ao, matrix_w) + NULLIFY (scrm, rho, dft_control, matrix_p_mp2, matrix_s, & + matrix_p_mp2_admm, admm_env, sab_orb, dbcsr_work_p, & + dbcsr_work_h, sac_ae, sac_ppl, sap_ppnl, force, virial, & + qs_kind_set, atomic_kind_set, particle_set, pw_env, poisson_env, & + auxbas_pw_pool, task_list_aux_fit, matrix_s_aux_fit, scrm_admm, & + rho_aux_fit, rho_ao_aux, x_data, hfx_section, xc_section, & + para_env, rhoz_g, rhoz_r, tauz_r, v_xc, v_xc_tau, vxc_rspace, & + vtau_rspace, vadmm_rspace, rho_ao, matrix_w) CALL timeset(routineN, handle) @@ -3202,13 +3201,11 @@ CONTAINS !The kinetic term (only response density) NULLIFY (scrm) - IF (nspins == 2) CALL dbcsr_add(matrix_p_mp2(1)%matrix, matrix_p_mp2(2)%matrix, 1.0_dp, 1.0_dp) - CALL build_kinetic_matrix(qs_env%ks_env, matrix_t=scrm, & - matrix_name="KINETIC ENERGY MATRIX", & - basis_type="ORB", & - sab_nl=sab_orb, calculate_forces=.TRUE., & - matrix_p=matrix_p_mp2(1)%matrix) - IF (nspins == 2) CALL dbcsr_add(matrix_p_mp2(1)%matrix, matrix_p_mp2(2)%matrix, 1.0_dp, -1.0_dp) + mpa2(1:nspins, 1:1) => matrix_p_mp2(1:nspins) + CALL kinetic_energy_matrix(qs_env, matrix_t=scrm, matrix_p=mpa2, & + matrix_name="KINETIC ENERGY MATRIX", & + basis_type="ORB", & + sab_orb=sab_orb, calculate_forces=.TRUE.) CALL dbcsr_deallocate_matrix_set(scrm) !The pseudo-potential terms (only reponse density) @@ -3222,34 +3219,14 @@ CONTAINS nder = 1 nimages = 1 - NULLIFY (cell_to_index) ALLOCATE (dbcsr_work_p(nspins, 1), dbcsr_work_h(nspins, 1)) DO ispin = 1, nspins dbcsr_work_p(ispin, 1)%matrix => matrix_p_mp2(ispin)%matrix dbcsr_work_h(ispin, 1)%matrix => scrm(ispin)%matrix END DO - IF (ASSOCIATED(sac_ae)) THEN - CALL build_core_ae(dbcsr_work_h, dbcsr_work_p, force, & - virial, .TRUE., use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, & - sab_orb, sac_ae, nimages, cell_to_index, "ORB") - END IF + CALL core_matrices(qs_env, dbcsr_work_h, dbcsr_work_p, .TRUE., nder) - IF (ASSOCIATED(sac_ppl)) THEN - CALL build_core_ppl(dbcsr_work_h, dbcsr_work_p, force, & - virial, .TRUE., use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, & - sab_orb, sac_ppl, nimages, cell_to_index, "ORB") - END IF - - IF (ASSOCIATED(sap_ppnl)) THEN - eps_ppnl = dft_control%qs_control%eps_ppnl - CALL build_core_ppnl(dbcsr_work_h, dbcsr_work_p, force, & - virial, .TRUE., use_virial, nder, & - qs_kind_set, atomic_kind_set, particle_set, & - sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, "ORB") - END IF DEALLOCATE (dbcsr_work_p, dbcsr_work_h) IF (use_virial) THEN