Prepare for External method response forces (#3891)

This commit is contained in:
Juerg Hutter 2025-01-24 16:45:43 +01:00 committed by GitHub
parent e123cd09bd
commit 269bd3eb43
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
5 changed files with 241 additions and 61 deletions

View file

@ -112,10 +112,13 @@ MODULE ec_env_types
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz => Null(), matrix_wz => Null(), &
matrix_z => Null(), z_admm => Null()
! Harris (rhoout), and response density (rhoz) on grid
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rhoout_r => Null(), rhoz_r => Null()
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rhoout_r => Null(), &
rhoz_r => Null()
! potentials from input density
TYPE(pw_r3d_rs_type) :: vh_rspace = pw_r3d_rs_type()
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace => Null(), vtau_rspace => Null(), vadmm_rspace => Null()
TYPE(pw_r3d_rs_type) :: vh_rspace = pw_r3d_rs_type()
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace => Null(), &
vtau_rspace => Null(), &
vadmm_rspace => Null()
! efield
TYPE(efield_berry_type), POINTER :: efield => NULL()
! LS matrices and types
@ -123,8 +126,15 @@ MODULE ec_env_types
! Environment for Hartree-Fock exchange
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data => Null()
! ADMM XC environments
TYPE(section_vals_type), POINTER :: xc_section_primary => Null(), &
xc_section_aux => Null()
TYPE(section_vals_type), POINTER :: xc_section_primary => Null(), &
xc_section_aux => Null()
! External
CHARACTER(len=40) :: exresp_fn = ""
CHARACTER(len=40) :: exresult_fn = ""
LOGICAL :: do_error = .FALSE.
REAL(KIND=dp), DIMENSION(3, 3) :: rpv = 0.0_dp, rpverror = 0.0_dp
REAL(KIND=dp), DIMENSION(:, :), POINTER :: rf => NULL(), &
rferror => NULL()
END TYPE energy_correction_type
CONTAINS
@ -215,6 +225,13 @@ CONTAINS
CALL cp_fm_release(ec_env%mo_occ)
CALL cp_fm_release(ec_env%cpmos)
IF (ASSOCIATED(ec_env%rf)) THEN
DEALLOCATE (ec_env%rf)
END IF
IF (ASSOCIATED(ec_env%rferror)) THEN
DEALLOCATE (ec_env%rferror)
END IF
DEALLOCATE (ec_env)
END IF

View file

@ -202,6 +202,13 @@ CONTAINS
l_val=ec_env%debug_external)
! ADMM
CALL section_vals_val_get(ec_section, "ADMM", l_val=ec_env%do_ec_admm)
! EXTERNAL
CALL section_vals_val_get(ec_section, "EXTERNAL_RESPONSE_FILENAME", &
c_val=ec_env%exresp_fn)
CALL section_vals_val_get(ec_section, "EXTERNAL_RESULT_FILENAME", &
c_val=ec_env%exresult_fn)
CALL section_vals_val_get(ec_section, "ERROR_ESTIMATION", &
l_val=ec_env%do_error)
ec_env%do_skip = .FALSE.

View file

@ -154,10 +154,9 @@ CONTAINS
ec_env%cpmos => cpmos
ELSE
! get external information
!deb CALL ec_ext_interface(qs_env, ec_env%ex, mo_ref, cpmos, calculate_forces, unit_nr)
CALL ec_ext_interface(qs_env, ec_env, mo_ref, cpmos, calculate_forces, unit_nr)
ec_env%mo_occ => mo_ref
ec_env%cpmos => cpmos
CPABORT("EC EXT NYA")
END IF
IF (calculate_forces) THEN
@ -184,6 +183,57 @@ CONTAINS
! **************************************************************************************************
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param ec_env ...
!> \param mo_ref ...
!> \param cpmos ...
!> \param calculate_forces ...
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE ec_ext_interface(qs_env, ec_env, mo_ref, cpmos, calculate_forces, unit_nr)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_ref, cpmos
LOGICAL, INTENT(IN) :: calculate_forces
INTEGER, INTENT(IN) :: unit_nr
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_interface'
INTEGER :: handle, ispin, nao, nocc(2), nspins
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(mp_para_env_type), POINTER :: para_env
CALL timeset(routineN, handle)
MARK_USED(qs_env)
MARK_USED(ec_env)
MARK_USED(mo_ref)
MARK_USED(cpmos)
MARK_USED(calculate_forces)
MARK_USED(unit_nr)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
nspins = dft_control%nspins
CALL get_qs_env(qs_env, mos=mos)
nocc = 0
DO ispin = 1, nspins
CALL get_mo_set(mo_set=mos(ispin), nao=nao, homo=nocc(ispin))
END DO
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A)') "Read EXTERNAL Response from file: "//TRIM(ec_env%exresp_fn)
END IF
CALL timestop(handle)
END SUBROUTINE ec_ext_interface
! **************************************************************************************************
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...

View file

@ -42,6 +42,8 @@ MODULE energy_corrections
cp_dbcsr_sm_fm_multiply,&
dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
USE cp_files, ONLY: close_file,&
open_file
USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
cp_fm_triangular_invert
USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
@ -97,10 +99,12 @@ MODULE energy_corrections
vdw_pairpot_dftd3bj, xc_vdw_fun_pairpot
USE input_section_types, ONLY: section_get_ival,&
section_get_lval,&
section_vals_duplicate,&
section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
section_vals_val_get,&
section_vals_val_set
USE kinds, ONLY: default_path_length,&
default_string_length,&
dp
@ -180,11 +184,13 @@ MODULE energy_corrections
USE task_list_methods, ONLY: generate_qs_task_list
USE task_list_types, ONLY: allocate_task_list,&
deallocate_task_list
USE trexio_utils, ONLY: write_trexio
USE virial_methods, ONLY: one_third_sum_diag,&
write_stress_tensor,&
write_stress_tensor_components
USE virial_types, ONLY: symmetrize_virial,&
virial_type
virial_type,&
zero_virial
USE voronoi_interface, ONLY: entry_voronoi_or_bqb
#include "./base/base_uses.f90"
@ -225,6 +231,7 @@ CONTAINS
TYPE(energy_correction_type), POINTER :: ec_env
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_force_type), DIMENSION(:), POINTER :: ks_force
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
@ -271,8 +278,9 @@ CONTAINS
WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 25), &
" Energy Correction Forces ", REPEAT("-", 26), "!"
END IF
CALL get_qs_env(qs_env, force=ks_force)
CALL get_qs_env(qs_env, force=ks_force, virial=virial)
CALL zero_qs_force(ks_force)
CALL zero_virial(virial, reset=.FALSE.)
ELSE
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 29), &
@ -405,6 +413,11 @@ CONTAINS
ec_env%matrix_s, &
ec_env%matrix_w)
CALL ec_build_ks_matrix_force(qs_env, ec_env)
IF (ec_env%debug_external) THEN
CALL write_response_interface(qs_env, ec_env)
CALL init_response_deriv(qs_env, ec_env)
END IF
CASE (ec_functional_dc)
! Prepare Density-corrected DFT (DC-DFT) calculation
@ -416,10 +429,15 @@ CONTAINS
ec_env%matrix_s, &
ec_env%matrix_w)
CALL ec_dc_build_ks_matrix_force(qs_env, ec_env)
IF (ec_env%debug_external) THEN
CALL write_response_interface(qs_env, ec_env)
CALL init_response_deriv(qs_env, ec_env)
END IF
CASE (ec_functional_ext)
CALL ec_ext_energy(qs_env, ec_env, calculate_forces=.TRUE.)
CALL init_response_deriv(qs_env, ec_env)
CASE DEFAULT
CPABORT("unknown energy correction")
@ -455,6 +473,8 @@ CONTAINS
p_env=ec_env%p_env, &
debug=debug_f)
CALL output_response_deriv(qs_env, ec_env, unit_nr)
CALL ec_properties(qs_env, ec_env)
! Deallocate Harris density and response density on grid
@ -486,6 +506,118 @@ CONTAINS
END SUBROUTINE energy_correction_low
! **************************************************************************************************
!> \brief Output response information to TREXIO file (for testing external method read in)
!> \param qs_env ...
!> \param ec_env ...
!> \author JHU
! **************************************************************************************************
SUBROUTINE write_response_interface(qs_env, ec_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
TYPE(section_vals_type), POINTER :: section, trexio_section
section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%TREXIO")
NULLIFY (trexio_section)
CALL section_vals_duplicate(section, trexio_section)
CALL section_vals_val_set(trexio_section, "FILENAME", c_val=ec_env%exresp_fn)
CALL section_vals_val_set(trexio_section, "CARTESIAN", l_val=.FALSE.)
CALL write_trexio(qs_env, trexio_section)
END SUBROUTINE write_response_interface
! **************************************************************************************************
!> \brief Initialize arrays for response derivatives
!> \param qs_env ...
!> \param ec_env ...
!> \author JHU
! **************************************************************************************************
SUBROUTINE init_response_deriv(qs_env, ec_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
INTEGER :: natom
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(virial_type), POINTER :: virial
CALL get_qs_env(qs_env, natom=natom)
ALLOCATE (ec_env%rf(3, natom), ec_env%rferror(3, natom))
ec_env%rf = 0.0_dp
ec_env%rferror = 0.0_dp
ec_env%rpv = 0.0_dp
ec_env%rpverror = 0.0_dp
CALL get_qs_env(qs_env, force=force, virial=virial)
CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
CALL total_qs_force(ec_env%rf, force, atomic_kind_set)
IF (virial%pv_availability .AND. (.NOT. virial%pv_numer)) THEN
ec_env%rpv = virial%pv_virial
END IF
END SUBROUTINE init_response_deriv
! **************************************************************************************************
!> \brief Write the reponse forces to file
!> \param qs_env ...
!> \param ec_env ...
!> \param unit_nr ...
!> \author JHU
! **************************************************************************************************
SUBROUTINE output_response_deriv(qs_env, ec_env, unit_nr)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
INTEGER, INTENT(IN) :: unit_nr
INTEGER :: funit, ia, natom
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(virial_type), POINTER :: virial
IF (ASSOCIATED(ec_env%rf)) THEN
CALL get_qs_env(qs_env, natom=natom)
ALLOCATE (ftot(3, natom))
ftot = 0.0_dp
CALL get_qs_env(qs_env, force=force, virial=virial)
CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
CALL total_qs_force(ftot, force, atomic_kind_set)
ec_env%rf(1:3, 1:natom) = ftot(1:3, 1:natom) - ec_env%rf(1:3, 1:natom)
DEALLOCATE (ftot)
IF (virial%pv_availability .AND. (.NOT. virial%pv_numer)) THEN
ec_env%rpv = virial%pv_virial - ec_env%rpv
END IF
CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell)
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A)') "Write EXTERNAL Response Derivative: "//TRIM(ec_env%exresult_fn)
CALL open_file(ec_env%exresult_fn, file_status="REPLACE", file_form="FORMATTED", &
file_action="WRITE", unit_number=funit)
WRITE (funit, *) " COORDINATES // RESPONSE FORCES // ERRORS "
DO ia = 1, natom
WRITE (funit, "(3(3F15.8,5x))") particle_set(ia)%r(1:3), &
ec_env%rf(1:3, ia), ec_env%rferror(1:3, ia)
END DO
WRITE (funit, *)
WRITE (funit, *) " CELL // RESPONSE PRESSURE // ERRORS "
DO ia = 1, 3
WRITE (funit, "(3F15.8,5x,3G15.8,5x,3G15.8)") cell%hmat(ia, 1:3), &
ec_env%rpv(ia, 1:3), ec_env%rpverror(ia, 1:3)
END DO
CALL close_file(funit)
END IF
END IF
END SUBROUTINE output_response_deriv
! **************************************************************************************************
!> \brief Calculates the traces of the core matrices and the density matrix.
!> \param qs_env ...
@ -2438,57 +2570,6 @@ CONTAINS
CPASSERT(.FALSE.)
END SELECT
! OUtput density available now
! HFX contribution to Harris functional and energy
! Can't calculate this earlier, cause ec_env%matrix_p doesnt exist yet
!------------------------------------------------------------------------------
! ! Exact exchange contribution (hybrid functionals)
! ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION")
! ec_hfx_sections => section_vals_get_subs_vals(ec_section, "XC%HF")
! CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
!
! IF (do_ec_hfx) THEN
!
! ! Check what works
! IF (dft_control%do_admm) THEN
! CALL cp_warn(__LOCATION__, "Energy correction with hybrid functional does not use ADMM.")
! END IF
!
! adiabatic_rescaling_section => section_vals_get_subs_vals(ec_section, "XC%ADIABATIC_RESCALING")
! CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
! IF (do_adiabatic_rescaling) THEN
! CALL cp_abort(__LOCATION__, "Adiabatic rescaling NYI for energy correction")
! END IF
! CALL section_vals_val_get(ec_hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core)
! IF (hfx_treat_lsd_in_core) THEN
! CALL cp_abort(__LOCATION__, "HFX_TREAT_LSD_IN_CORE NYI for energy correction")
! END IF
!
! ! Exchange matrix
! IF (ASSOCIATED(ec_env%matrix_x)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_x)
! CALL dbcsr_allocate_matrix_set(ec_env%matrix_x, nspins)
! DO ispin = 1, nspins
! headline = "EXCHANGE MATRIX"
! ALLOCATE (ec_env%matrix_x(ispin)%matrix)
! CALL dbcsr_create(ec_env%matrix_x(ispin)%matrix, name=TRIM(headline), &
! template=ec_env%matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
! CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_x(ispin)%matrix, ec_env%sab_orb)
! CALL dbcsr_set(ec_env%matrix_x(ispin)%matrix, 0.0_dp)
! END DO
!
! ! Get exact exchange energy (fraction) and its contribution to the EC hamiltonian
! should_update=.TRUE.
! ks_mat => ec_env%matrix_ks(:,1)
! CALL ec_hfx_contributions(qs_env, ks_mat, matrix_p, &
! ec_hfx_sections, ec_env%x_data, use_virial, &
! should_update, calculate_forces, matrix_x = ec_env%matrix_x, ex = ec_env%ex)
!
! END IF
!------------------------------------------------------------------------------
IF (ec_env%mao) THEN
CALL mao_release_matrices(ec_env, ksmat, smat, pmat, wmat)
END IF

View file

@ -361,6 +361,31 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EXTERNAL_RESPONSE_FILENAME", &
description="Name of the file that contains response information.", &
usage="EXTERNAL_RESPONSE_FILENAME <FILENAME>", &
default_c_val="TREXIO")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EXTERNAL_RESULT_FILENAME", &
description="Name of the file that contains results from external response calculation.", &
usage="EXTERNAL_RESULT_FILENAME <FILENAME>", &
default_c_val="CP2K_EXRESP.result")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create( &
keyword, __LOCATION__, &
name="ERROR_ESTIMATION", &
description="Perform an error estimation for the response forces/stress. "// &
"Requires error estimates for the RHS of the response equation from input. ", &
usage="ERROR_ESTIMATION", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_ec_section
! **************************************************************************************************