diff --git a/src/core_ae.F b/src/core_ae.F index dc6e826083..1563e406b7 100644 --- a/src/core_ae.F +++ b/src/core_ae.F @@ -78,11 +78,12 @@ CONTAINS !> \param sac_ae ... !> \param nimages ... !> \param cell_to_index ... +!> \param basis_type ... !> \param atcore ... ! ************************************************************************************************** SUBROUTINE 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, atcore) + nimages, cell_to_index, basis_type, atcore) TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p TYPE(qs_force_type), DIMENSION(:), POINTER :: force @@ -97,6 +98,7 @@ CONTAINS POINTER :: sab_orb, sac_ae INTEGER, INTENT(IN) :: nimages INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index + CHARACTER(LEN=*), INTENT(IN) :: basis_type REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), & OPTIONAL :: atcore @@ -168,7 +170,7 @@ CONTAINS ALLOCATE (basis_set_list(nkind)) DO ikind = 1, nkind - CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a) + CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type) IF (ASSOCIATED(basis_set_a)) THEN basis_set_list(ikind)%gto_basis_set => basis_set_a ELSE @@ -176,7 +178,7 @@ CONTAINS END IF END DO - CALL get_qs_kind_set(qs_kind_set, & + CALL get_qs_kind_set(qs_kind_set, basis_type=basis_type, & maxco=maxco, maxlgto=maxl, maxsgf=maxsgf, maxnset=maxnset) CALL init_orbital_pointers(maxl + nder + 1) ldsab = MAX(maxco, maxsgf) diff --git a/src/ec_env_types.F b/src/ec_env_types.F index 2ea73ca044..0d8690001f 100644 --- a/src/ec_env_types.F +++ b/src/ec_env_types.F @@ -88,7 +88,8 @@ MODULE ec_env_types TYPE(qs_force_type), DIMENSION(:), POINTER :: force => Null() ! full neighbor lists and corresponding task list TYPE(neighbor_list_set_p_type), & - DIMENSION(:), POINTER :: sab_orb => Null(), sac_ppl => Null(), sap_ppnl => Null() + DIMENSION(:), POINTER :: sab_orb => Null(), sac_ppl => Null(), & + sac_ae => Null(), sap_ppnl => Null() TYPE(task_list_type), POINTER :: task_list => Null() ! the XC function to be used for the correction, dispersion info TYPE(section_vals_type), POINTER :: xc_section => Null() @@ -155,6 +156,7 @@ CONTAINS ! neighbor lists CALL release_neighbor_list_sets(ec_env%sab_orb) CALL release_neighbor_list_sets(ec_env%sac_ppl) + CALL release_neighbor_list_sets(ec_env%sac_ae) CALL release_neighbor_list_sets(ec_env%sap_ppnl) ! forces IF (ASSOCIATED(ec_env%force)) CALL deallocate_qs_force(ec_env%force) diff --git a/src/ec_environment.F b/src/ec_environment.F index b963378fa5..8ee7335313 100644 --- a/src/ec_environment.F +++ b/src/ec_environment.F @@ -133,7 +133,7 @@ CONTAINS CALL timeset(routineN, handle) NULLIFY (atomic_kind_set, dispersion_env, ec_env%ls_env, para_env) - NULLIFY (ec_env%sab_orb, ec_env%sac_ppl, ec_env%sap_ppnl) + NULLIFY (ec_env%sab_orb, ec_env%sac_ae, ec_env%sac_ppl, ec_env%sap_ppnl) NULLIFY (ec_env%matrix_ks, ec_env%matrix_h, ec_env%matrix_s) NULLIFY (ec_env%matrix_t, ec_env%matrix_p, ec_env%matrix_w) NULLIFY (ec_env%task_list) diff --git a/src/ec_external.F b/src/ec_external.F index ac0e771873..0b9a1c2d8f 100644 --- a/src/ec_external.F +++ b/src/ec_external.F @@ -14,6 +14,7 @@ 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 @@ -120,6 +121,11 @@ CONTAINS unit_nr = -1 END IF + ec_env%etotal = 0.0_dp + ec_env%ex = 0.0_dp + ec_env%exc = 0.0_dp + ec_env%ehartree = 0.0_dp + CALL get_qs_env(qs_env, mos=mos) ALLOCATE (cpmos(nspins), mo_ref(nspins), mo_occ(nspins)) @@ -367,7 +373,7 @@ CONTAINS 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_ppl, sap_ppnl + 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 @@ -410,9 +416,10 @@ CONTAINS END IF ! get neighbor lists - NULLIFY (sab_orb, sac_ppl, sap_ppnl) + NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl) CALL get_qs_env(qs_env=qs_env, & - sab_orb=sab_orb, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl) + sab_orb=sab_orb, sac_ae=sac_ae, & + sac_ppl=sac_ppl, sap_ppnl=sap_ppnl) ! initialize src matrix NULLIFY (scrm) @@ -456,6 +463,26 @@ CONTAINS 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 diff --git a/src/energy_corrections.F b/src/energy_corrections.F index 6f75895a61..1831558eca 100644 --- a/src/energy_corrections.F +++ b/src/energy_corrections.F @@ -87,7 +87,8 @@ MODULE energy_corrections USE ec_external, ONLY: ec_ext_energy USE ec_methods, ONLY: create_kernel,& ec_mos_init - USE external_potential_types, ONLY: get_potential,& + USE external_potential_types, ONLY: all_potential_type,& + get_potential,& gth_potential_type,& sgp_potential_type USE hfx_exx, ONLY: add_exx_to_rhs,& @@ -147,7 +148,8 @@ MODULE energy_corrections USE qs_dispersion_types, ONLY: qs_dispersion_type USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env,& - qs_environment_type + qs_environment_type,& + set_qs_env USE qs_force_types, ONLY: qs_force_type,& total_qs_force,& zero_qs_force @@ -172,6 +174,10 @@ MODULE energy_corrections build_neighbor_lists,& local_atoms_type,& pair_radius_setup + USE qs_oce_methods, ONLY: build_oce_matrices + USE qs_oce_types, ONLY: allocate_oce_set,& + create_oce_set,& + oce_matrix_type USE qs_ot_eigensolver, ONLY: ot_eigensolver USE qs_overlap, ONLY: build_overlap_matrix USE qs_rho_types, ONLY: qs_rho_get,& @@ -338,12 +344,17 @@ CONTAINS LOGICAL, INTENT(IN) :: calculate_forces INTEGER, INTENT(IN) :: unit_nr - INTEGER :: ispin, nspins - LOGICAL :: debug_f - REAL(KIND=dp) :: exc + INTEGER :: ispin, nkind, nspins + LOGICAL :: debug_f, gapw, gapw_xc + REAL(KIND=dp) :: eps_fit, exc TYPE(dft_control_type), POINTER :: dft_control + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sap_oce + TYPE(oce_matrix_type), POINTER :: oce + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(pw_env_type), POINTER :: pw_env TYPE(pw_pool_type), POINTER :: auxbas_pw_pool + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set IF (ec_env%should_update) THEN CALL ec_build_neighborlist(qs_env, ec_env) @@ -443,12 +454,27 @@ CONTAINS CPABORT("unknown energy correction") END SELECT + CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) + nspins = dft_control%nspins + gapw = dft_control%qs_control%gapw + gapw_xc = dft_control%qs_control%gapw_xc + IF (gapw .OR. gapw_xc) THEN + CALL get_qs_env(qs_env=qs_env, nkind=nkind, & + qs_kind_set=qs_kind_set, particle_set=particle_set) + CALL get_qs_env(qs_env=qs_env, oce=oce, sap_oce=sap_oce) + CALL create_oce_set(oce) + CALL allocate_oce_set(oce, nkind) + eps_fit = dft_control%qs_control%gapw_control%eps_fit + CALL build_oce_matrices(oce%intac, .TRUE., 1, qs_kind_set, particle_set, & + sap_oce, eps_fit) + CALL set_qs_env(qs_env, oce=oce) + END IF + CALL response_calculation(qs_env, ec_env) ! Allocate response density on real space grid for use in properties ! Calculated in response_force - CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env) - nspins = dft_control%nspins + CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) CPASSERT(ASSOCIATED(pw_env)) CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) @@ -571,7 +597,7 @@ CONTAINS CHARACTER(LEN=default_string_length) :: unit_string INTEGER :: funit, ia, natom - REAL(KIND=dp) :: fconv + REAL(KIND=dp) :: evol, fconv REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell @@ -595,6 +621,11 @@ CONTAINS IF (virial%pv_availability .AND. (.NOT. virial%pv_numer)) THEN ec_env%rpv = virial%pv_virial - ec_env%rpv CALL para_env%sum(ec_env%rpv) + ! Volume terms + evol = ec_env%exc + ec_env%exc_aux_fit + 2.0_dp*ec_env%ehartree + ec_env%rpv(1, 1) = ec_env%rpv(1, 1) - evol + ec_env%rpv(2, 2) = ec_env%rpv(2, 2) - evol + ec_env%rpv(3, 3) = ec_env%rpv(3, 3) - evol END IF CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell) @@ -1334,6 +1365,7 @@ CONTAINS ! get neighbor lists, we need the full sab_orb list from the ec_env NULLIFY (sab_orb, sac_ae, 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 @@ -1361,13 +1393,10 @@ CONTAINS ! Possible AE contributions (also ECP) IF (ASSOCIATED(sac_ae)) THEN - CPABORT("ECP/AE not available for energy corrections") - ! missig code: sac_ae has to bee calculated and stored in ec_env - ! build_core_ae: needs functionality to set the basis type at input 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) + nimages, cell_to_index, "HARRIS") END IF ! compute the ppl contribution to the core hamiltonian IF (ASSOCIATED(sac_ppl)) THEN @@ -1645,7 +1674,7 @@ CONTAINS 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_ppl, sap_ppnl + 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 @@ -1690,12 +1719,15 @@ CONTAINS ! check for GAPW/GAPW_XC IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN - CPABORT("Harris functional for GAPW not implemented") + IF (ec_env%energy_functional == ec_functional_harris) THEN + CPABORT("Harris functional for GAPW not implemented") + END IF 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 @@ -1763,6 +1795,26 @@ CONTAINS 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, "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 @@ -2928,15 +2980,17 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_neighborlist' INTEGER :: handle, ikind, nkind, zat - LOGICAL :: gth_potential_present, & + LOGICAL :: all_potential_present, & + gth_potential_present, & sgp_potential_present, & skip_load_balance_distributed - LOGICAL, ALLOCATABLE, DIMENSION(:) :: default_present, orb_present, & - ppl_present, ppnl_present + LOGICAL, ALLOCATABLE, DIMENSION(:) :: all_present, default_present, & + orb_present, ppl_present, ppnl_present REAL(dp) :: subcells - REAL(dp), ALLOCATABLE, DIMENSION(:) :: c_radius, orb_radius, ppl_radius, & - ppnl_radius + REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_radius, c_radius, orb_radius, & + ppl_radius, ppnl_radius REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius + TYPE(all_potential_type), POINTER :: all_potential TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell TYPE(dft_control_type), POINTER :: dft_control @@ -2958,12 +3012,13 @@ CONTAINS CALL timeset(routineN, handle) CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set) - CALL get_qs_kind_set(qs_kind_set, gth_potential_present=gth_potential_present, & + CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present, & + gth_potential_present=gth_potential_present, & sgp_potential_present=sgp_potential_present) nkind = SIZE(qs_kind_set) ALLOCATE (c_radius(nkind), default_present(nkind)) - ALLOCATE (orb_radius(nkind), ppl_radius(nkind), ppnl_radius(nkind)) - ALLOCATE (orb_present(nkind), ppl_present(nkind), ppnl_present(nkind)) + ALLOCATE (orb_radius(nkind), all_radius(nkind), ppl_radius(nkind), ppnl_radius(nkind)) + ALLOCATE (orb_present(nkind), all_present(nkind), ppl_present(nkind), ppnl_present(nkind)) ALLOCATE (pair_radius(nkind, nkind)) ALLOCATE (atom2d(nkind)) @@ -2989,7 +3044,8 @@ CONTAINS orb_present(ikind) = .FALSE. orb_radius(ikind) = 0.0_dp END IF - CALL get_qs_kind(qs_kind, gth_potential=gth_potential, sgp_potential=sgp_potential) + CALL get_qs_kind(qs_kind, all_potential=all_potential, & + gth_potential=gth_potential, sgp_potential=sgp_potential) IF (gth_potential_present .OR. sgp_potential_present) THEN IF (ASSOCIATED(gth_potential)) THEN CALL get_potential(potential=gth_potential, & @@ -3010,6 +3066,20 @@ CONTAINS ppnl_radius(ikind) = 0.0_dp END IF END IF + ! Check the presence of an all electron potential or ERFC potential + IF (all_potential_present .OR. sgp_potential_present) THEN + all_present(ikind) = .FALSE. + all_radius(ikind) = 0.0_dp + IF (ASSOCIATED(all_potential)) THEN + all_present(ikind) = .TRUE. + CALL get_potential(potential=all_potential, core_charge_radius=all_radius(ikind)) + ELSE IF (ASSOCIATED(sgp_potential)) THEN + IF (sgp_potential%ecp_local) THEN + all_present(ikind) = .TRUE. + CALL get_potential(potential=sgp_potential, core_charge_radius=all_radius(ikind)) + END IF + END IF + END IF END DO CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells) @@ -3018,7 +3088,15 @@ CONTAINS CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius) CALL build_neighbor_lists(ec_env%sab_orb, particle_set, atom2d, cell, pair_radius, & subcells=subcells, nlname="sab_orb") - ! pseudopotential + ! pseudopotential/AE + IF (all_potential_present .OR. sgp_potential_present) THEN + IF (ANY(all_present)) THEN + CALL pair_radius_setup(orb_present, all_present, orb_radius, all_radius, pair_radius) + CALL build_neighbor_lists(ec_env%sac_ae, particle_set, atom2d, cell, pair_radius, & + subcells=subcells, operator_type="ABC", nlname="sac_ae") + END IF + END IF + IF (gth_potential_present .OR. sgp_potential_present) THEN IF (ANY(ppl_present)) THEN CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius) @@ -3062,8 +3140,8 @@ CONTAINS ! Release work storage CALL atom2d_cleanup(atom2d) DEALLOCATE (atom2d) - DEALLOCATE (orb_present, default_present, ppl_present, ppnl_present) - DEALLOCATE (orb_radius, ppl_radius, ppnl_radius, c_radius) + DEALLOCATE (orb_present, default_present, all_present, ppl_present, ppnl_present) + DEALLOCATE (orb_radius, all_radius, ppl_radius, ppnl_radius, c_radius) DEALLOCATE (pair_radius) ! Task list diff --git a/src/mp2_cphf.F b/src/mp2_cphf.F index eec03953fd..b37de4e3e5 100644 --- a/src/mp2_cphf.F +++ b/src/mp2_cphf.F @@ -1452,7 +1452,7 @@ CONTAINS 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) + 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) diff --git a/src/qs_core_hamiltonian.F b/src/qs_core_hamiltonian.F index 4c1a973d08..02929b2a30 100644 --- a/src/qs_core_hamiltonian.F +++ b/src/qs_core_hamiltonian.F @@ -490,7 +490,7 @@ CONTAINS 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, atcore=atcore) + nimages, cell_to_index, "ORB", atcore=atcore) END IF ! *** compute the ppl contribution to the core hamiltonian *** diff --git a/src/qs_dcdr_ao.F b/src/qs_dcdr_ao.F index 2b4231d76e..f2b6963bda 100644 --- a/src/qs_dcdr_ao.F +++ b/src/qs_dcdr_ao.F @@ -460,7 +460,7 @@ CONTAINS 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, & - cell_to_index=cell_to_index) + 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) diff --git a/src/response_solver.F b/src/response_solver.F index 2644aad461..04409f2bf9 100644 --- a/src/response_solver.F +++ b/src/response_solver.F @@ -831,7 +831,7 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'response_force' - CHARACTER(LEN=default_string_length) :: basis_type + 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 @@ -1033,7 +1033,7 @@ CONTAINS 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) + 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) @@ -2503,6 +2503,13 @@ CONTAINS CALL para_env%sum(stdeb) IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & 'STRESS| Stress Response ', one_third_sum_diag(stdeb), det_3x3(stdeb) + stdeb = fconv*(virial%pv_virial) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Total Stress ', one_third_sum_diag(stdeb), det_3x3(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,3(1X,ES19.11))") & + stdeb(1, 1), stdeb(2, 2), stdeb(3, 3) + unitstr = "bar" END IF IF (do_ex) THEN diff --git a/src/rpa_im_time_force_methods.F b/src/rpa_im_time_force_methods.F index eabf62de30..c583c831bc 100644 --- a/src/rpa_im_time_force_methods.F +++ b/src/rpa_im_time_force_methods.F @@ -3233,7 +3233,7 @@ CONTAINS 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) + sab_orb, sac_ae, nimages, cell_to_index, "ORB") END IF IF (ASSOCIATED(sac_ppl)) THEN