(git:97501a3)
Loading...
Searching...
No Matches
qs_apt_fdiff_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Atomic Polarization Tensor calculation by dF/d(E-field) finite differences
10!> \author Leo Decking, Hossam Elgabarty
11! **************************************************************************************************
12
13MODULE qs_apt_fdiff_methods
14 USE cp_control_types, ONLY: dft_control_type
15 USE cp_log_handling, ONLY: cp_add_default_logger,&
16 cp_get_default_logger,&
17 cp_logger_create,&
18 cp_logger_get_default_io_unit,&
19 cp_logger_release,&
20 cp_logger_set,&
21 cp_logger_type,&
22 cp_rm_default_logger
23 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
24 cp_print_key_unit_nr
25 USE cp_subsys_types, ONLY: cp_subsys_get,&
26 cp_subsys_type
27 USE force_env_types, ONLY: force_env_get,&
28 force_env_type
29 USE input_section_types, ONLY: section_vals_get_subs_vals,&
30 section_vals_type,&
31 section_vals_val_get
32 USE kinds, ONLY: dp
33 USE message_passing, ONLY: mp_para_env_type
34 USE particle_list_types, ONLY: particle_list_type
35 USE qs_apt_fdiff_types, ONLY: apt_fdiff_point_type,&
36 apt_fdiff_points_type
37 USE qs_environment_types, ONLY: get_qs_env,&
38 qs_environment_type
39 USE qs_force, ONLY: qs_calc_energy_force
40 USE qs_subsys_types, ONLY: qs_subsys_type
41#include "./base/base_uses.f90"
42
43 IMPLICIT NONE
44 PRIVATE
45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_apt_fdiff_methods'
46 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
47 PUBLIC :: apt_fdiff
48
49CONTAINS
50
51! **************************************************************************************************
52!> \brief Perform the 2PNT symmetric difference
53!> \param apt_fdiff_points ...
54!> \param ap_tensor ...
55!> \param natoms ...
56!> \author Leo Decking
57! **************************************************************************************************
58 SUBROUTINE fdiff_2pnt(apt_fdiff_points, ap_tensor, natoms)
59 TYPE(apt_fdiff_points_type) :: apt_fdiff_points
60 REAL(kind=dp), DIMENSION(:, :, :) :: ap_tensor
61 INTEGER, INTENT(IN) :: natoms
62
63 INTEGER :: i, j, n
64
65 DO j = 1, 3 ! axis force
66 DO i = 1, 3 ! axis field
67 DO n = 1, natoms
68 ap_tensor(n, i, j) = (apt_fdiff_points%point_field(i, 1)%forces(n, j) - &
69 apt_fdiff_points%point_field(i, 2)%forces(n, j)) &
70 /(2*apt_fdiff_points%field_strength)
71 END DO
72 END DO
73 END DO
74
75 END SUBROUTINE fdiff_2pnt
76
77! **************************************************************************************************
78!> \brief ...
79!> \param apt_fdiff_point ...
80!> \param particles ...
81!> \author Leo Decking
82! **************************************************************************************************
83 SUBROUTINE get_forces(apt_fdiff_point, particles)
84 TYPE(apt_fdiff_point_type) :: apt_fdiff_point
85 TYPE(particle_list_type), POINTER :: particles
86
87 INTEGER :: i
88
89 cpassert(ASSOCIATED(particles))
90
91 DO i = 1, particles%n_els
92 apt_fdiff_point%forces(i, 1:3) = particles%els(i)%f(1:3)
93 END DO
94
95 END SUBROUTINE get_forces
96
97! **************************************************************************************************
98!> \brief Calculate Atomic Polarization Tensors by dF/d(E-field) finite differences
99!> \param force_env ...
100!> \author Leo Decking, Hossam Elgabarty
101! **************************************************************************************************
102 SUBROUTINE apt_fdiff(force_env)
103
104 TYPE(force_env_type), POINTER :: force_env
105
106 CHARACTER(LEN=*), PARAMETER :: routinen = 'apt_fdiff'
107
108 INTEGER :: apt_log, fd_method, handle, i, j, &
109 log_unit, n, natoms, output_fdiff_scf
110 REAL(kind=dp) :: born_sum
111 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ap_tensor
112 TYPE(apt_fdiff_points_type) :: apt_fdiff_points
113 TYPE(cp_logger_type), POINTER :: logger, logger_apt
114 TYPE(cp_subsys_type), POINTER :: cp_subsys
115 TYPE(dft_control_type), POINTER :: dft_control
116 TYPE(mp_para_env_type), POINTER :: para_env
117 TYPE(particle_list_type), POINTER :: particles
118 TYPE(qs_environment_type), POINTER :: qs_env
119 TYPE(qs_subsys_type), POINTER :: subsys
120 TYPE(section_vals_type), POINTER :: dcdr_section
121
122 CALL timeset(routinen, handle)
123
124 NULLIFY (qs_env, logger, dcdr_section, logger_apt, para_env, dft_control, subsys)
125 NULLIFY (particles)
126
127 cpassert(ASSOCIATED(force_env))
128
129 CALL force_env_get(force_env, subsys=cp_subsys, qs_env=qs_env, para_env=para_env)
130
131 CALL cp_subsys_get(cp_subsys, particles=particles)
132 cpassert(ASSOCIATED(particles))
133 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)!, para_env=para_env)
134
135 natoms = particles%n_els
136 IF (dft_control%apply_period_efield .OR. dft_control%apply_efield) THEN
137 cpabort("APT calculation not available in the presence of an external electric field")
138 END IF
139
140 logger => cp_get_default_logger()
141
142 log_unit = cp_logger_get_default_io_unit(logger)
143
144 dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
145
146 apt_log = cp_print_key_unit_nr(logger, dcdr_section, "PRINT%APT", &
147 extension=".data", middle_name="apt", log_filename=.false., &
148 file_position="APPEND", file_status="UNKNOWN")
149
150 output_fdiff_scf = cp_print_key_unit_nr(logger, dcdr_section, "PRINT%APT", &
151 extension=".scfLog", middle_name="apt", log_filename=.false., &
152 file_position="APPEND", file_status="UNKNOWN")
153
154 CALL cp_logger_create(logger_apt, para_env=para_env, print_level=0, &
155 default_global_unit_nr=output_fdiff_scf, &
156 close_global_unit_on_dealloc=.true.)
157
158 CALL cp_logger_set(logger_apt, global_filename="APT_localLog")
159
160 CALL cp_add_default_logger(logger_apt)
161
162 IF (output_fdiff_scf > 0) THEN
163 WRITE (output_fdiff_scf, '(T2,A)') &
164 '!----------------------------------------------------------------------------!'
165 WRITE (output_fdiff_scf, '(/,T2,A)') "SCF log for finite difference steps"
166 WRITE (output_fdiff_scf, '(/,T2,A)') &
167 '!----------------------------------------------------------------------------!'
168 END IF
169
170 IF (log_unit > 0) THEN
171 WRITE (log_unit, '(/,T2,A)') &
172 '!----------------------------------------------------------------------------!'
173 WRITE (log_unit, '(/,T10,A)') "Computing Atomic polarization tensors using finite differences"
174 WRITE (log_unit, '(T2,A)') " "
175 END IF
176
177 CALL section_vals_val_get(dcdr_section, "APT_FD_DE", r_val=apt_fdiff_points%field_strength)
178 CALL section_vals_val_get(dcdr_section, "APT_FD_METHOD", i_val=fd_method)
179
180 dft_control%apply_period_efield = .true.
181 ALLOCATE (dft_control%period_efield)
182 dft_control%period_efield%displacement_field = .false.
183
184 DO i = 1, 3
185 dft_control%period_efield%polarisation(1:3) = (/0.0_dp, 0.0_dp, 0.0_dp/)
186 dft_control%period_efield%polarisation(i) = 1.0_dp
187
188 IF (log_unit > 0) THEN
189 WRITE (log_unit, '(T2,A)') " "
190 WRITE (log_unit, "(T2,A)") "Computing forces under efield in direction +/-"//achar(i + 119)
191 END IF
192
193 DO j = 1, 2 ! 1 -> positive, 2 -> negative
194 IF (j == 1) THEN
195 dft_control%period_efield%strength = apt_fdiff_points%field_strength
196 ELSE
197 dft_control%period_efield%strength = -1.0*apt_fdiff_points%field_strength
198 END IF
199
200 CALL qs_calc_energy_force(qs_env, calc_force=.true., consistent_energies=.true., linres=.false.)
201
202 ALLOCATE (apt_fdiff_points%point_field(i, j)%forces(natoms, 1:3))
203 CALL get_forces(apt_fdiff_points%point_field(i, j), particles)
204
205 END DO
206 END DO
207
208 IF (output_fdiff_scf > 0) THEN
209 WRITE (output_fdiff_scf, '(/,T2,A)') &
210 '!----------------------------------------------------------------------------!'
211 WRITE (output_fdiff_scf, '(/,T2,A)') "Finite differences done!"
212 WRITE (output_fdiff_scf, '(/,T2,A)') &
213 '!----------------------------------------------------------------------------!'
214 END IF
215
216 CALL cp_print_key_finished_output(output_fdiff_scf, logger_apt, dcdr_section, "PRINT%APT")
217 CALL cp_logger_release(logger_apt)
218 CALL cp_rm_default_logger()
219
220 ALLOCATE (ap_tensor(natoms, 3, 3))
221 CALL fdiff_2pnt(apt_fdiff_points, ap_tensor, natoms)
222
223 !Print
224 born_sum = 0.0_dp
225 IF (apt_log > 0) THEN
226 DO n = 1, natoms
227 born_sum = born_sum + (ap_tensor(n, 1, 1) + ap_tensor(n, 1, 1) + ap_tensor(n, 1, 1))/3.0
228 WRITE (apt_log, "(I6, A6, F20.10)") n, particles%els(n)%atomic_kind%element_symbol, &
229 (ap_tensor(n, 1, 1) + ap_tensor(n, 2, 2) + ap_tensor(n, 3, 3))/3.0
230 WRITE (apt_log, '(F20.10,F20.10,F20.10)') &
231 ap_tensor(n, 1, 1), ap_tensor(n, 1, 2), ap_tensor(n, 1, 3)
232 WRITE (apt_log, '(F20.10,F20.10,F20.10)') &
233 ap_tensor(n, 2, 1), ap_tensor(n, 2, 2), ap_tensor(n, 2, 3)
234 WRITE (apt_log, '(F20.10,F20.10,F20.10)') &
235 ap_tensor(n, 3, 1), ap_tensor(n, 3, 2), ap_tensor(n, 3, 3)
236 END DO
237 WRITE (apt_log, '(/,A20, F20.10)') "Sum of Born charges:", born_sum
238 WRITE (log_unit, '(/,A30, F20.10)') "Checksum (Acoustic Sum Rule):", born_sum
239 END IF
240 DEALLOCATE (ap_tensor)
241 DEALLOCATE (dft_control%period_efield)
242
243 CALL cp_print_key_finished_output(apt_log, logger, dcdr_section, "PRINT%APT")
244
245 dft_control%apply_period_efield = .false.
246 qs_env%linres_run = .false.
247
248 IF (log_unit > 0) THEN
249 WRITE (log_unit, '(T2,A)') " "
250 WRITE (log_unit, '(T2,A)') " "
251 WRITE (log_unit, '(T22,A)') "APT calculation Done!"
252 WRITE (log_unit, '(/,T2,A)') &
253 '!----------------------------------------------------------------------------!'
254 END IF
255
256 CALL timestop(handle)
257
258 END SUBROUTINE apt_fdiff
259
260END MODULE qs_apt_fdiff_methods
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_set
subroutine, public cp_logger_set(logger, local_filename, global_filename)
sets various attributes of the given logger
Definition cp_log_handling.F:717
cp_log_handling::cp_rm_default_logger
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
Definition cp_log_handling.F:213
cp_log_handling::cp_logger_release
subroutine, public cp_logger_release(logger)
releases this logger
Definition cp_log_handling.F:423
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_log_handling::cp_logger_create
subroutine, public cp_logger_create(logger, para_env, print_level, default_global_unit_nr, default_local_unit_nr, global_filename, local_filename, close_global_unit_on_dealloc, iter_info, close_local_unit_on_dealloc, suffix, template_logger)
initializes a logger
Definition cp_log_handling.F:285
cp_log_handling::cp_add_default_logger
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
Definition cp_log_handling.F:188
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:345
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
qs_apt_fdiff_methods
Atomic Polarization Tensor calculation by dF/d(E-field) finite differences.
Definition qs_apt_fdiff_methods.F:13
qs_apt_fdiff_methods::apt_fdiff
subroutine, public apt_fdiff(force_env)
Calculate Atomic Polarization Tensors by dF/d(E-field) finite differences.
Definition qs_apt_fdiff_methods.F:103
qs_apt_fdiff_types
Atomic Polarization Tensor calculation by dF/d(E-field) finite differences.
Definition qs_apt_fdiff_types.F:13
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:520
qs_force
Quickstep force driver routine.
Definition qs_force.F:12
qs_force::qs_calc_energy_force
subroutine, public qs_calc_energy_force(qs_env, calc_force, consistent_energies, linres)
...
Definition qs_force.F:108
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
cp_control_types::dft_control_type
Definition cp_control_types.F:599
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
qs_apt_fdiff_types::apt_fdiff_point_type
Definition qs_apt_fdiff_types.F:25
qs_apt_fdiff_types::apt_fdiff_points_type
Definition qs_apt_fdiff_types.F:29
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:219
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56