diff --git a/docs/acronyms.md b/docs/acronyms.md index 9f712fdd1f..dcb1c6275b 100644 --- a/docs/acronyms.md +++ b/docs/acronyms.md @@ -13,6 +13,9 @@ AM1 AMBER Assisted Model Building and Energy Refinement +APT + Atomic Polarization Tensor + ASE Atomic Simulation Environment diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b367221172..f3f2c38342 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -511,6 +511,8 @@ list( qcschema.F hdf5_wrapper.F sockets_interface.F + qs_apt_fdiff_methods.F + qs_apt_fdiff_types.F qs_atomic_block.F qs_band_structure.F qs_basis_gradient.F diff --git a/src/force_env_methods.F b/src/force_env_methods.F index 03bdb153f2..9597562a8c 100644 --- a/src/force_env_methods.F +++ b/src/force_env_methods.F @@ -39,7 +39,9 @@ MODULE force_env_methods cp_logger_type,& cp_rm_default_logger,& cp_to_string - USE cp_output_handling, ONLY: cp_print_key_finished_output,& + USE cp_output_handling, ONLY: cp_p_file,& + cp_print_key_finished_output,& + cp_print_key_should_output,& cp_print_key_unit_nr,& low_print_level USE cp_result_methods, ONLY: cp_results_erase,& @@ -143,6 +145,7 @@ MODULE force_env_methods USE qmmm_util, ONLY: apply_qmmm_translate USE qmmmx_force, ONLY: qmmmx_calc_energy_force USE qmmmx_types, ONLY: qmmmx_env_type + USE qs_apt_fdiff_methods, ONLY: apt_fdiff USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type,& @@ -205,8 +208,8 @@ CONTAINS nfixed_atoms_total, nkind, & output_unit, print_forces, print_grrm, & print_scine - LOGICAL :: calculate_forces, calculate_stress_tensor, energy_consistency, eval_ef, & - linres_run, my_skip, print_components + LOGICAL :: calculate_forces, calculate_stress_tensor, do_apt_fd, energy_consistency, & + eval_ef, linres_run, my_skip, print_components REAL(KIND=dp) :: checksum, e_entropy, e_gap, e_pot, & fconv, sum_energy REAL(KIND=dp), DIMENSION(3) :: grand_total_force, total_force @@ -219,6 +222,7 @@ CONTAINS TYPE(molecule_kind_type), POINTER :: molecule_kind TYPE(particle_list_type), POINTER :: core_particles, particles, & shell_particles + TYPE(section_vals_type), POINTER :: print_key TYPE(virial_type), POINTER :: virial NULLIFY (logger, virial, subsys, atprop_env, cell) @@ -310,6 +314,19 @@ CONTAINS END IF END IF + ! In case requested, compute the APT numerically + do_apt_FD = .FALSE. + IF (force_env%in_use == use_qs_force) THEN + CALL section_vals_val_get(force_env%qs_env%input, "PROPERTIES%LINRES%DCDR%APT_FD", l_val=do_apt_FD) + IF (do_apt_FD) THEN + print_key => section_vals_get_subs_vals(force_env%qs_env%input, & + subsection_name="PROPERTIES%LINRES%DCDR%PRINT%APT") + IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN + CALL apt_fdiff(force_env) + END IF + END IF + END IF + !sample peak memory CALL m_memory() diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index da81b0b2c9..d91a13a09b 100644 --- a/src/input_cp2k_properties_dft.F +++ b/src/input_cp2k_properties_dft.F @@ -414,6 +414,39 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="APT_FD", & + description="Use numerical differentiation to compute the APT, "// & + "switches off the calculation of dcdr analytical derivatives. "// & + "Requires RUN_TYPE = ENERGY_FORCE or MD.", & + usage="APT_FD T", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="APT_FD_DE", & + description="Electric field strength (atomic units) to use for finite differences", & + repeats=.FALSE., & + n_var=1, & + type_of_var=real_t, & + default_r_val=0.0003_dp, & + usage="APT_FD_DE 1.0E-4") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="APT_FD_METHOD", & + description="Numerical differentiation method", & + usage="APT_FD_METHOD FD", & + default_i_val=1, & + !enum_c_vals=s2a("FD", "2PNT"), & + enum_c_vals=s2a("2PNT"), & + !enum_desc=s2a("Forward differences.", & + ! "Symmetric two-point differences."), & + enum_desc=s2a("Symmetric two-point differences."), & + !enum_i_vals=(/0, 1/)) + enum_i_vals=(/1/)) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + NULLIFY (subsection) CALL section_create(subsection, __LOCATION__, name="PRINT", & description="print results of the magnetic dipole moment calculation", & diff --git a/src/qs_apt_fdiff_methods.F b/src/qs_apt_fdiff_methods.F new file mode 100644 index 0000000000..11b3a4ecac --- /dev/null +++ b/src/qs_apt_fdiff_methods.F @@ -0,0 +1,260 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2025 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Atomic Polarization Tensor calculation by dF/d(E-field) finite differences +!> \author Leo Decking, Hossam Elgabarty +! ************************************************************************************************** + +MODULE qs_apt_fdiff_methods + USE cp_control_types, ONLY: dft_control_type + USE cp_log_handling, ONLY: cp_add_default_logger,& + cp_get_default_logger,& + cp_logger_create,& + cp_logger_get_default_io_unit,& + cp_logger_release,& + cp_logger_set,& + cp_logger_type,& + cp_rm_default_logger + USE cp_output_handling, ONLY: cp_print_key_finished_output,& + cp_print_key_unit_nr + USE cp_subsys_types, ONLY: cp_subsys_get,& + cp_subsys_type + USE force_env_types, ONLY: force_env_get,& + force_env_type + USE input_section_types, ONLY: section_vals_get_subs_vals,& + section_vals_type,& + section_vals_val_get + USE kinds, ONLY: dp + USE message_passing, ONLY: mp_para_env_type + USE particle_list_types, ONLY: particle_list_type + USE qs_apt_fdiff_types, ONLY: apt_fdiff_point_type,& + apt_fdiff_points_type + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type + USE qs_force, ONLY: qs_calc_energy_force + USE qs_subsys_types, ONLY: qs_subsys_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + PRIVATE + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_apt_fdiff_methods' + LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. + PUBLIC :: apt_fdiff + +CONTAINS + +! ************************************************************************************************** +!> \brief Perform the 2PNT symmetric difference +!> \param apt_fdiff_points ... +!> \param ap_tensor ... +!> \param natoms ... +!> \author Leo Decking +! ************************************************************************************************** + SUBROUTINE fdiff_2pnt(apt_fdiff_points, ap_tensor, natoms) + TYPE(apt_fdiff_points_type) :: apt_fdiff_points + REAL(kind=dp), DIMENSION(:, :, :) :: ap_tensor + INTEGER, INTENT(IN) :: natoms + + INTEGER :: i, j, n + + DO j = 1, 3 ! axis force + DO i = 1, 3 ! axis field + DO n = 1, natoms + ap_tensor(n, i, j) = (apt_fdiff_points%point_field(i, 1)%forces(n, j) - & + apt_fdiff_points%point_field(i, 2)%forces(n, j)) & + /(2*apt_fdiff_points%field_strength) + END DO + END DO + END DO + + END SUBROUTINE fdiff_2pnt + +! ************************************************************************************************** +!> \brief ... +!> \param apt_fdiff_point ... +!> \param particles ... +!> \author Leo Decking +! ************************************************************************************************** + SUBROUTINE get_forces(apt_fdiff_point, particles) + TYPE(apt_fdiff_point_type) :: apt_fdiff_point + TYPE(particle_list_type), POINTER :: particles + + INTEGER :: i + + CPASSERT(ASSOCIATED(particles)) + + DO i = 1, particles%n_els + apt_fdiff_point%forces(i, 1:3) = particles%els(i)%f(1:3) + END DO + + END SUBROUTINE get_forces + +! ************************************************************************************************** +!> \brief Calculate Atomic Polarization Tensors by dF/d(E-field) finite differences +!> \param force_env ... +!> \author Leo Decking, Hossam Elgabarty +! ************************************************************************************************** + SUBROUTINE apt_fdiff(force_env) + + TYPE(force_env_type), POINTER :: force_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'apt_fdiff' + + INTEGER :: apt_log, fd_method, handle, i, j, & + log_unit, n, natoms, output_fdiff_scf + REAL(kind=dp) :: born_sum + REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ap_tensor + TYPE(apt_fdiff_points_type) :: apt_fdiff_points + TYPE(cp_logger_type), POINTER :: logger, logger_apt + TYPE(cp_subsys_type), POINTER :: cp_subsys + TYPE(dft_control_type), POINTER :: dft_control + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(particle_list_type), POINTER :: particles + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(qs_subsys_type), POINTER :: subsys + TYPE(section_vals_type), POINTER :: dcdr_section + + CALL timeset(routineN, handle) + + NULLIFY (qs_env, logger, dcdr_section, logger_apt, para_env, dft_control, subsys) + NULLIFY (particles) + + CPASSERT(ASSOCIATED(force_env)) + + CALL force_env_get(force_env, subsys=cp_subsys, qs_env=qs_env, para_env=para_env) + + CALL cp_subsys_get(cp_subsys, particles=particles) + CPASSERT(ASSOCIATED(particles)) + CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)!, para_env=para_env) + + natoms = particles%n_els + IF (dft_control%apply_period_efield .OR. dft_control%apply_efield) THEN + CPABORT("APT calculation not available in the presence of an external electric field") + END IF + + logger => cp_get_default_logger() + + log_unit = cp_logger_get_default_io_unit(logger) + + dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR") + + apt_log = cp_print_key_unit_nr(logger, dcdr_section, "PRINT%APT", & + extension=".data", middle_name="apt", log_filename=.FALSE., & + file_position="APPEND", file_status="UNKNOWN") + + output_fdiff_scf = cp_print_key_unit_nr(logger, dcdr_section, "PRINT%APT", & + extension=".scfLog", middle_name="apt", log_filename=.FALSE., & + file_position="APPEND", file_status="UNKNOWN") + + CALL cp_logger_create(logger_apt, para_env=para_env, print_level=0, & + default_global_unit_nr=output_fdiff_scf, & + close_global_unit_on_dealloc=.TRUE.) + + CALL cp_logger_set(logger_apt, global_filename="APT_localLog") + + CALL cp_add_default_logger(logger_apt) + + IF (output_fdiff_scf > 0) THEN + WRITE (output_fdiff_scf, '(T2,A)') & + '!----------------------------------------------------------------------------!' + WRITE (output_fdiff_scf, '(/,T2,A)') "SCF log for finite difference steps" + WRITE (output_fdiff_scf, '(/,T2,A)') & + '!----------------------------------------------------------------------------!' + END IF + + IF (log_unit > 0) THEN + WRITE (log_unit, '(/,T2,A)') & + '!----------------------------------------------------------------------------!' + WRITE (log_unit, '(/,T10,A)') "Computing Atomic polarization tensors using finite differences" + WRITE (log_unit, '(T2,A)') " " + END IF + + CALL section_vals_val_get(dcdr_section, "APT_FD_DE", r_val=apt_fdiff_points%field_strength) + CALL section_vals_val_get(dcdr_section, "APT_FD_METHOD", i_val=fd_method) + + dft_control%apply_period_efield = .TRUE. + ALLOCATE (dft_control%period_efield) + dft_control%period_efield%displacement_field = .FALSE. + + DO i = 1, 3 + dft_control%period_efield%polarisation(1:3) = (/0.0_dp, 0.0_dp, 0.0_dp/) + dft_control%period_efield%polarisation(i) = 1.0_dp + + IF (log_unit > 0) THEN + WRITE (log_unit, '(T2,A)') " " + WRITE (log_unit, "(T2,A)") "Computing forces under efield in direction +/-"//ACHAR(i + 119) + END IF + + DO j = 1, 2 ! 1 -> positive, 2 -> negative + IF (j == 1) THEN + dft_control%period_efield%strength = apt_fdiff_points%field_strength + ELSE + dft_control%period_efield%strength = -1.0*apt_fdiff_points%field_strength + END IF + + CALL qs_calc_energy_force(qs_env, calc_force=.TRUE., consistent_energies=.TRUE., linres=.FALSE.) + + ALLOCATE (apt_fdiff_points%point_field(i, j)%forces(natoms, 1:3)) + CALL get_forces(apt_fdiff_points%point_field(i, j), particles) + + END DO + END DO + + IF (output_fdiff_scf > 0) THEN + WRITE (output_fdiff_scf, '(/,T2,A)') & + '!----------------------------------------------------------------------------!' + WRITE (output_fdiff_scf, '(/,T2,A)') "Finite differences done!" + WRITE (output_fdiff_scf, '(/,T2,A)') & + '!----------------------------------------------------------------------------!' + END IF + + CALL cp_print_key_finished_output(output_fdiff_scf, logger_apt, dcdr_section, "PRINT%APT") + CALL cp_logger_release(logger_apt) + CALL cp_rm_default_logger() + + ALLOCATE (ap_tensor(natoms, 3, 3)) + CALL fdiff_2pnt(apt_fdiff_points, ap_tensor, natoms) + + !Print + born_sum = 0.0_dp + IF (apt_log > 0) THEN + DO n = 1, natoms + born_sum = born_sum + (ap_tensor(n, 1, 1) + ap_tensor(n, 1, 1) + ap_tensor(n, 1, 1))/3.0 + WRITE (apt_log, "(I6, A6, F20.10)") n, particles%els(n)%atomic_kind%element_symbol, & + (ap_tensor(n, 1, 1) + ap_tensor(n, 2, 2) + ap_tensor(n, 3, 3))/3.0 + WRITE (apt_log, '(F20.10,F20.10,F20.10)') & + ap_tensor(n, 1, 1), ap_tensor(n, 1, 2), ap_tensor(n, 1, 3) + WRITE (apt_log, '(F20.10,F20.10,F20.10)') & + ap_tensor(n, 2, 1), ap_tensor(n, 2, 2), ap_tensor(n, 2, 3) + WRITE (apt_log, '(F20.10,F20.10,F20.10)') & + ap_tensor(n, 3, 1), ap_tensor(n, 3, 2), ap_tensor(n, 3, 3) + END DO + WRITE (apt_log, '(/,A20, F20.10)') "Sum of Born charges:", born_sum + WRITE (log_unit, '(/,A30, F20.10)') "Checksum (Acoustic Sum Rule):", born_sum + END IF + DEALLOCATE (ap_tensor) + DEALLOCATE (dft_control%period_efield) + + CALL cp_print_key_finished_output(apt_log, logger, dcdr_section, "PRINT%APT") + + dft_control%apply_period_efield = .FALSE. + qs_env%linres_run = .FALSE. + + IF (log_unit > 0) THEN + WRITE (log_unit, '(T2,A)') " " + WRITE (log_unit, '(T2,A)') " " + WRITE (log_unit, '(T22,A)') "APT calculation Done!" + WRITE (log_unit, '(/,T2,A)') & + '!----------------------------------------------------------------------------!' + END IF + + CALL timestop(handle) + + END SUBROUTINE apt_fdiff + +END MODULE qs_apt_fdiff_methods diff --git a/src/qs_apt_fdiff_types.F b/src/qs_apt_fdiff_types.F new file mode 100644 index 0000000000..180c75b44a --- /dev/null +++ b/src/qs_apt_fdiff_types.F @@ -0,0 +1,35 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2025 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Atomic Polarization Tensor calculation by dF/d(E-field) finite differences +!> \author Leo Decking, Hossam Elgabarty +! ************************************************************************************************** + +MODULE qs_apt_fdiff_types + + USE kinds, ONLY: dp + !USE physcon, ONLY: debye +#include "./base/base_uses.f90" + + IMPLICIT NONE + PRIVATE + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_apt_fdiff_types' + LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. + PUBLIC :: apt_fdiff_point_type, apt_fdiff_points_type + + TYPE apt_fdiff_point_type + REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: forces + END TYPE apt_fdiff_point_type + + TYPE apt_fdiff_points_type + !TYPE(apt_fdiff_point_type) :: point_no_field ! Only for forward differences + TYPE(apt_fdiff_point_type), DIMENSION(1:3, 1:2) :: point_field + REAL(kind=dp) :: field_strength + END TYPE apt_fdiff_points_type + +END MODULE qs_apt_fdiff_types diff --git a/src/qs_linres_module.F b/src/qs_linres_module.F index 4d1ac8fa85..02716e6ba6 100644 --- a/src/qs_linres_module.F +++ b/src/qs_linres_module.F @@ -357,9 +357,9 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_calculation_low' INTEGER :: every_n_step, handle, iounit - LOGICAL :: dcdr_present, epr_present, issc_present, & - lr_calculation, nmr_present, & - polar_present, vcd_present + LOGICAL :: dcdr_present, do_apt_fd, epr_present, & + issc_present, lr_calculation, & + nmr_present, polar_present, vcd_present TYPE(cp_logger_type), POINTER :: logger TYPE(dft_control_type), POINTER :: dft_control TYPE(linres_control_type), POINTER :: linres_control @@ -374,13 +374,21 @@ CONTAINS issc_present = .FALSE. polar_present = .FALSE. dcdr_present = .FALSE. + do_apt_fd = .FALSE. NULLIFY (dft_control, linres_control, logger, prop_section, lr_section) - logger => cp_get_default_logger() lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") CALL section_vals_get(lr_section, explicit=lr_calculation) + CALL section_vals_val_get(lr_section, "DCDR%APT_FD", explicit=do_apt_fd) + IF (do_apt_fd) THEN + CALL timestop(handle) + RETURN + END IF + + logger => cp_get_default_logger() + CALL section_vals_val_get(lr_section, "EVERY_N_STEP", i_val=every_n_step) IF (lr_calculation .AND. MODULO(qs_env%sim_step, every_n_step) == 0) THEN diff --git a/tests/QS/regtest-dcdr/TEST_FILES.toml b/tests/QS/regtest-dcdr/TEST_FILES.toml index 77b528ca35..e4f9daaed0 100644 --- a/tests/QS/regtest-dcdr/TEST_FILES.toml +++ b/tests/QS/regtest-dcdr/TEST_FILES.toml @@ -12,4 +12,5 @@ "h2o_apt_pbc.inp" = [{matcher="M096", tol=1e-06, ref=-0.938707}] "h2o_apt_pbc_loc.inp" = [{matcher="M096", tol=1e-06, ref=-0.938770}] "h2o_aat.inp" = [{matcher="M100", tol=1e-06, ref=-0.857427}] +"h2o_apt_fdiff.inp" = [{matcher="M127", tol=1e-04, ref=0.0034319}] #EOF diff --git a/tests/QS/regtest-dcdr/h2o_apt_fdiff.inp b/tests/QS/regtest-dcdr/h2o_apt_fdiff.inp new file mode 100644 index 0000000000..0ed05806b4 --- /dev/null +++ b/tests/QS/regtest-dcdr/h2o_apt_fdiff.inp @@ -0,0 +1,82 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT wat-dim-apt-fdiff + RUN_TYPE ENERGY + #RUN_TYPE MD +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + CHARGE 0 + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 200 + NGRIDS 2 + &END MGRID + &POISSON + PERIODIC XYZ + POISSON_SOLVER PERIODIC + &END POISSON + &QS + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-3 + MAX_SCF 50 + SCF_GUESS ATOMIC + &OT + PRECONDITIONER FULL_SINGLE_INVERSE + &END OT + &END SCF + &XC + &XC_FUNCTIONAL LDA + &END XC_FUNCTIONAL + &END XC + &END DFT + &PROPERTIES + &LINRES + &DCDR + APT_FD + # &PRINT + # &APT + # &EACH + # MD 5 + # &END EACH + # &END APT + # &END PRINT + &END DCDR + &END LINRES + &END PROPERTIES + &SUBSYS + &CELL + ABC [angstrom] 8.0 5.0 5.0 + PERIODIC XYZ + &END CELL + &COORD + O 7.0414782577 8.1944693132 4.7203167161 + O 6.4136423281 6.9680225843 7.2942345321 + H 7.0024272528 9.1522284981 4.7945554450 + H 6.3623068415 7.9578974691 4.0821913660 + H 7.1357809996 6.4074215266 7.5863285535 + H 6.7003933205 7.3362896088 6.4455943873 + &END COORD + &KIND DEFAULT + BASIS_SET SZV-GTH + POTENTIAL GTH-LDA + &END KIND + &END SUBSYS +&END FORCE_EVAL + +# &MOTION +# &MD +# STEPS 20 +# TIMESTEP 0.5 +# &THERMOSTAT +# &CSVR +# TIMECON 5 +# &END CSVR +# &END +# &END MD +# &END MOTION \ No newline at end of file diff --git a/tests/matchers.py b/tests/matchers.py index 57b95ad917..ee13a7730f 100644 --- a/tests/matchers.py +++ b/tests/matchers.py @@ -257,4 +257,6 @@ registry["RTBSE_GXAC_H2_pol"] = GenericMatcher( ) registry["M126"] = GenericMatcher(r" # Total charge ", col=5) + +registry["M127"] = GenericMatcher(r"Checksum (Acoustic Sum Rule):", col=5) # EOF