Atomic polarization tensors via numerical differentiation

Co-authored-by: Leo Decking <ldecking@mail.uni-paderborn.de>
This commit is contained in:
Hossam Elgabarty 2025-07-21 19:40:43 +02:00 committed by GitHub
parent ae167406ec
commit 716159d91c
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
10 changed files with 450 additions and 7 deletions

View file

@ -13,6 +13,9 @@ AM1
AMBER
Assisted Model Building and Energy Refinement
APT
Atomic Polarization Tensor
ASE
Atomic Simulation Environment

View file

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

View file

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

View file

@ -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", &

260
src/qs_apt_fdiff_methods.F Normal file
View file

@ -0,0 +1,260 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! 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

35
src/qs_apt_fdiff_types.F Normal file
View file

@ -0,0 +1,35 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! 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

View file

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

View file

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

View file

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

View file

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