Skip to content

Commit

Permalink
Add ground state forces to TDDFPT forces in property calculation (#2386)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Nov 7, 2022
1 parent 0a7d28b commit a213519
Show file tree
Hide file tree
Showing 5 changed files with 335 additions and 176 deletions.
31 changes: 28 additions & 3 deletions src/qs_energy.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,19 +14,20 @@
MODULE qs_energy
USE almo_scf, ONLY: almo_entry_scf
USE cp_control_types, ONLY: dft_control_type
USE cp_external_control, ONLY: external_control
USE dm_ls_scf, ONLY: ls_scf
USE energy_corrections, ONLY: energy_correction
USE excited_states, ONLY: excited_state_energy
USE lri_environment_methods, ONLY: lri_print_stat
USE mp2, ONLY: mp2_main
USE qs_energy_init, ONLY: qs_energies_init
USE qs_energy_types, ONLY: qs_energy_type
USE qs_energy_utils, ONLY: qs_energies_compute_matrix_w,&
qs_energies_mp2,&
qs_energies_properties
USE qs_energy_utils, ONLY: qs_energies_properties
USE qs_environment_methods, ONLY: qs_env_rebuild_pw_env
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_ks_methods, ONLY: qs_ks_update_qs_env
USE qs_matrix_w, ONLY: qs_energies_compute_matrix_w
USE qs_scf, ONLY: scf
#include "./base/base_uses.f90"

Expand Down Expand Up @@ -127,4 +128,28 @@ SUBROUTINE qs_energies(qs_env, consistent_energies, calc_forces)

END SUBROUTINE qs_energies

! **************************************************************************************************
!> \brief Enters the mp2 part of cp2k
!> \param qs_env ...
!> \param calc_forces ...
! **************************************************************************************************

SUBROUTINE qs_energies_mp2(qs_env, calc_forces)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN) :: calc_forces

LOGICAL :: should_stop

! Compute MP2 energy

IF (ASSOCIATED(qs_env%mp2_env)) THEN

CALL external_control(should_stop, "MP2", target_time=qs_env%target_time, &
start_time=qs_env%start_time)

CALL mp2_main(qs_env=qs_env, calc_forces=calc_forces)
END IF

END SUBROUTINE qs_energies_mp2

END MODULE qs_energy
171 changes: 2 additions & 169 deletions src/qs_energy_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,47 +16,27 @@ MODULE qs_energy_utils
atprop_type
USE cp_control_types, ONLY: dft_control_type
USE cp_control_utils, ONLY: read_ddapc_section
USE cp_external_control, ONLY: external_control
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_p_type,&
cp_fm_release,&
cp_fm_type
USE dbcsr_api, ONLY: &
dbcsr_get_block_p, dbcsr_get_info, 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_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_type
USE et_coupling, ONLY: calc_et_coupling
USE et_coupling_proj, ONLY: calc_et_coupling_proj
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type
USE kinds, ONLY: dp
USE kpoint_methods, ONLY: kpoint_density_matrices,&
kpoint_density_transform
USE kpoint_types, ONLY: kpoint_type
USE mp2, ONLY: mp2_main
USE pw_env_types, ONLY: pw_env_type
USE pw_types, ONLY: pw_type
USE qs_density_matrices, ONLY: calculate_w_matrix,&
calculate_w_matrix_ot
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_integrate_potential, ONLY: integrate_v_core_rspace
USE qs_ks_methods, ONLY: qs_ks_update_qs_env
USE qs_linres_module, ONLY: linres_calculation_low
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_p_type,&
mo_set_type
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE qs_scf, ONLY: scf
USE qs_tddfpt2_methods, ONLY: tddfpt
USE scf_control_types, ONLY: scf_control_type
USE tip_scan_methods, ONLY: tip_scanning
USE xas_methods, ONLY: xas
USE xas_tdp_methods, ONLY: xas_tdp
Expand All @@ -70,133 +50,10 @@ MODULE qs_energy_utils

CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_energy_utils'

PUBLIC :: qs_energies_compute_matrix_w, &
qs_energies_properties, &
qs_energies_mp2
PUBLIC :: qs_energies_properties

CONTAINS

! **************************************************************************************************
!> \brief Refactoring of qs_energies_scf. Moves computation of matrix_w
!> into separate subroutine
!> \param qs_env ...
!> \param calc_forces ...
!> \par History
!> 05.2013 created [Florian Schiffmann]
! **************************************************************************************************

SUBROUTINE qs_energies_compute_matrix_w(qs_env, calc_forces)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN) :: calc_forces

CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_compute_matrix_w'

INTEGER :: handle, is, ispin, nao, nspin
LOGICAL :: do_kpoints, has_unit_metric
TYPE(cp_fm_p_type), DIMENSION(2) :: fmwork
TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_w, &
mo_derivs, rho_ao
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, matrix_w_kp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kpoint_type), POINTER :: kpoints
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos
TYPE(mo_set_type), POINTER :: mo_set
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_nl
TYPE(qs_rho_type), POINTER :: rho
TYPE(scf_control_type), POINTER :: scf_control

CALL timeset(routineN, handle)

! if calculate forces, time to compute the w matrix
CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)

IF (calc_forces .AND. .NOT. has_unit_metric) THEN
CALL get_qs_env(qs_env, do_kpoints=do_kpoints)

IF (do_kpoints) THEN

CALL get_qs_env(qs_env, &
matrix_w_kp=matrix_w_kp, &
matrix_s_kp=matrix_s_kp, &
sab_orb=sab_nl, &
mos=mos, &
kpoints=kpoints)

CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff, nao=nao)
CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, nrow_global=nao, ncol_global=nao, &
template_fmstruct=mo_coeff%matrix_struct)

DO is = 1, SIZE(fmwork)
ALLOCATE (fmwork(is)%matrix)
CALL cp_fm_create(fmwork(is)%matrix, matrix_struct=ao_ao_fmstruct)
END DO
CALL cp_fm_struct_release(ao_ao_fmstruct)

! energy weighted density matrices in k-space
CALL kpoint_density_matrices(kpoints, energy_weighted=.TRUE.)
! energy weighted density matrices in real space
CALL kpoint_density_transform(kpoints, matrix_w_kp, .TRUE., &
matrix_s_kp(1, 1)%matrix, sab_nl, fmwork)

DO is = 1, SIZE(fmwork)
CALL cp_fm_release(fmwork(is)%matrix)
DEALLOCATE (fmwork(is)%matrix)
END DO

ELSE

NULLIFY (dft_control, rho_ao)
CALL get_qs_env(qs_env, &
matrix_w=matrix_w, &
matrix_ks=matrix_ks, &
matrix_s=matrix_s, &
mo_derivs=mo_derivs, &
scf_control=scf_control, &
mos=mos, &
rho=rho, &
dft_control=dft_control)

CALL qs_rho_get(rho, rho_ao=rho_ao)

nspin = SIZE(mos)
DO ispin = 1, nspin
mo_set => mos(ispin)%mo_set
IF (dft_control%roks) THEN
IF (scf_control%use_ot) THEN
IF (ispin > 1) THEN
! not very elegant, indeed ...
CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
ELSE
CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, &
matrix_w(ispin)%matrix, matrix_s(1)%matrix)
END IF
ELSE
CALL calculate_w_matrix(mo_set=mo_set, &
matrix_ks=matrix_ks(ispin)%matrix, &
matrix_p=rho_ao(ispin)%matrix, &
matrix_w=matrix_w(ispin)%matrix)
END IF
ELSE
IF (scf_control%use_ot) THEN
CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, &
matrix_w(ispin)%matrix, matrix_s(1)%matrix)
ELSE
CALL calculate_w_matrix(mo_set, matrix_w(ispin)%matrix)
END IF
END IF
END DO

END IF

END IF

CALL timestop(handle)

END SUBROUTINE qs_energies_compute_matrix_w

! **************************************************************************************************
!> \brief Refactoring of qs_energies_scf. Moves computation of properties
!> into separate subroutine
Expand Down Expand Up @@ -394,28 +251,4 @@ SUBROUTINE atom_trace(amat, bmat, factor, atrace)

END SUBROUTINE atom_trace

! **************************************************************************************************
!> \brief Enters the mp2 part of cp2k
!> \param qs_env ...
!> \param calc_forces ...
! **************************************************************************************************

SUBROUTINE qs_energies_mp2(qs_env, calc_forces)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN) :: calc_forces

LOGICAL :: should_stop

! Compute MP2 energy

IF (ASSOCIATED(qs_env%mp2_env)) THEN

CALL external_control(should_stop, "MP2", target_time=qs_env%target_time, &
start_time=qs_env%start_time)

CALL mp2_main(qs_env=qs_env, calc_forces=calc_forces)
END IF

END SUBROUTINE qs_energies_mp2

END MODULE qs_energy_utils
2 changes: 1 addition & 1 deletion src/qs_force.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ SUBROUTINE qs_calc_energy_force(qs_env, calc_force, consistent_energies, linres)
CALL qs_energies(qs_env, calc_forces=.FALSE., &
consistent_energies=consistent_energies)
END IF
CALL get_qs_env(qs_env)

END SUBROUTINE qs_calc_energy_force

! **************************************************************************************************
Expand Down

0 comments on commit a213519

Please sign in to comment.