diff --git a/src/ec_env_types.F b/src/ec_env_types.F index 500ee859d8..f9e9c996af 100644 --- a/src/ec_env_types.F +++ b/src/ec_env_types.F @@ -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 diff --git a/src/ec_environment.F b/src/ec_environment.F index a380c35f8b..b963378fa5 100644 --- a/src/ec_environment.F +++ b/src/ec_environment.F @@ -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. diff --git a/src/ec_external.F b/src/ec_external.F index 5e6378c35f..c6ae525e1f 100644 --- a/src/ec_external.F +++ b/src/ec_external.F @@ -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 ... diff --git a/src/energy_corrections.F b/src/energy_corrections.F index 7d5a6d3114..06fcdda801 100644 --- a/src/energy_corrections.F +++ b/src/energy_corrections.F @@ -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 diff --git a/src/input_cp2k_ec.F b/src/input_cp2k_ec.F index 24c4067680..6adc516a7b 100644 --- a/src/input_cp2k_ec.F +++ b/src/input_cp2k_ec.F @@ -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 ", & + 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 ", & + 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 ! **************************************************************************************************