comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

ECP/bug fix for energy correction (#4275)

Browse source
This commit is contained in:
Juerg Hutter 2025-07-03 13:49:09 +02:00 • committed by GitHub
parent 93bc24a11a
commit eb49f8a800
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
10 changed files with 156 additions and 40 deletions
Download patch file Download diff file Expand all files Collapse all files

8
src/core_ae.F
View file

@ -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)

4
src/ec_env_types.F
View file

@ -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)

2
src/ec_environment.F
View file

@ -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)

33
src/ec_external.F
View file

@ -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

130
src/energy_corrections.F
View file

@ -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

2
src/mp2_cphf.F
View file

@ -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)

2
src/qs_core_hamiltonian.F
View file

@ -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 ***

2
src/qs_dcdr_ao.F
View file

@ -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)

11
src/response_solver.F
View file

@ -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

2
src/rpa_im_time_force_methods.F
View file

@ -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