Line data Source code
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 :
13 : MODULE 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 :
49 : CONTAINS
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 2 : 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 8 : DO j = 1, 3 ! axis force
66 26 : DO i = 1, 3 ! axis field
67 132 : 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 126 : /(2*apt_fdiff_points%field_strength)
71 : END DO
72 : END DO
73 : END DO
74 :
75 2 : END SUBROUTINE fdiff_2pnt
76 :
77 : ! **************************************************************************************************
78 : !> \brief ...
79 : !> \param apt_fdiff_point ...
80 : !> \param particles ...
81 : !> \author Leo Decking
82 : ! **************************************************************************************************
83 12 : 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 12 : CPASSERT(ASSOCIATED(particles))
90 :
91 84 : DO i = 1, particles%n_els
92 300 : apt_fdiff_point%forces(i, 1:3) = particles%els(i)%f(1:3)
93 : END DO
94 :
95 12 : 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 2 : 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 2 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ap_tensor
112 18 : 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 2 : CALL timeset(routineN, handle)
123 :
124 2 : NULLIFY (qs_env, logger, dcdr_section, logger_apt, para_env, dft_control, subsys)
125 2 : NULLIFY (particles)
126 :
127 2 : CPASSERT(ASSOCIATED(force_env))
128 :
129 2 : CALL force_env_get(force_env, subsys=cp_subsys, qs_env=qs_env, para_env=para_env)
130 :
131 2 : CALL cp_subsys_get(cp_subsys, particles=particles)
132 2 : CPASSERT(ASSOCIATED(particles))
133 2 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)!, para_env=para_env)
134 :
135 2 : natoms = particles%n_els
136 2 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield) THEN
137 0 : CPABORT("APT calculation not available in the presence of an external electric field")
138 : END IF
139 :
140 2 : logger => cp_get_default_logger()
141 :
142 2 : log_unit = cp_logger_get_default_io_unit(logger)
143 :
144 2 : 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 2 : 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 2 : 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 2 : close_global_unit_on_dealloc=.TRUE.)
157 :
158 2 : CALL cp_logger_set(logger_apt, global_filename="APT_localLog")
159 :
160 2 : CALL cp_add_default_logger(logger_apt)
161 :
162 2 : IF (output_fdiff_scf > 0) THEN
163 : WRITE (output_fdiff_scf, '(T2,A)') &
164 1 : '!----------------------------------------------------------------------------!'
165 1 : WRITE (output_fdiff_scf, '(/,T2,A)') "SCF log for finite difference steps"
166 : WRITE (output_fdiff_scf, '(/,T2,A)') &
167 1 : '!----------------------------------------------------------------------------!'
168 : END IF
169 :
170 2 : IF (log_unit > 0) THEN
171 : WRITE (log_unit, '(/,T2,A)') &
172 1 : '!----------------------------------------------------------------------------!'
173 1 : WRITE (log_unit, '(/,T10,A)') "Computing Atomic polarization tensors using finite differences"
174 1 : WRITE (log_unit, '(T2,A)') " "
175 : END IF
176 :
177 2 : CALL section_vals_val_get(dcdr_section, "APT_FD_DE", r_val=apt_fdiff_points%field_strength)
178 2 : CALL section_vals_val_get(dcdr_section, "APT_FD_METHOD", i_val=fd_method)
179 :
180 2 : dft_control%apply_period_efield = .TRUE.
181 18 : ALLOCATE (dft_control%period_efield)
182 2 : dft_control%period_efield%displacement_field = .FALSE.
183 :
184 8 : DO i = 1, 3
185 24 : dft_control%period_efield%polarisation(1:3) = [0.0_dp, 0.0_dp, 0.0_dp]
186 6 : dft_control%period_efield%polarisation(i) = 1.0_dp
187 :
188 6 : IF (log_unit > 0) THEN
189 3 : WRITE (log_unit, '(T2,A)') " "
190 3 : WRITE (log_unit, "(T2,A)") "Computing forces under efield in direction +/-"//ACHAR(i + 119)
191 : END IF
192 :
193 20 : DO j = 1, 2 ! 1 -> positive, 2 -> negative
194 12 : IF (j == 1) THEN
195 6 : dft_control%period_efield%strength = apt_fdiff_points%field_strength
196 : ELSE
197 6 : dft_control%period_efield%strength = -1.0*apt_fdiff_points%field_strength
198 : END IF
199 :
200 12 : CALL qs_calc_energy_force(qs_env, calc_force=.TRUE., consistent_energies=.TRUE., linres=.FALSE.)
201 :
202 36 : ALLOCATE (apt_fdiff_points%point_field(i, j)%forces(natoms, 1:3))
203 18 : CALL get_forces(apt_fdiff_points%point_field(i, j), particles)
204 :
205 : END DO
206 : END DO
207 :
208 2 : IF (output_fdiff_scf > 0) THEN
209 : WRITE (output_fdiff_scf, '(/,T2,A)') &
210 1 : '!----------------------------------------------------------------------------!'
211 1 : WRITE (output_fdiff_scf, '(/,T2,A)') "Finite differences done!"
212 : WRITE (output_fdiff_scf, '(/,T2,A)') &
213 1 : '!----------------------------------------------------------------------------!'
214 : END IF
215 :
216 2 : CALL cp_print_key_finished_output(output_fdiff_scf, logger_apt, dcdr_section, "PRINT%APT")
217 2 : CALL cp_logger_release(logger_apt)
218 2 : CALL cp_rm_default_logger()
219 :
220 8 : ALLOCATE (ap_tensor(natoms, 3, 3))
221 2 : CALL fdiff_2pnt(apt_fdiff_points, ap_tensor, natoms)
222 :
223 : !Print
224 2 : born_sum = 0.0_dp
225 2 : IF (apt_log > 0) THEN
226 7 : DO n = 1, natoms
227 6 : born_sum = born_sum + (ap_tensor(n, 1, 1) + ap_tensor(n, 1, 1) + ap_tensor(n, 1, 1))/3.0
228 6 : WRITE (apt_log, "(I6, A6, F20.10)") n, particles%els(n)%atomic_kind%element_symbol, &
229 12 : (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 6 : 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 6 : 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 7 : ap_tensor(n, 3, 1), ap_tensor(n, 3, 2), ap_tensor(n, 3, 3)
236 : END DO
237 1 : WRITE (apt_log, '(/,A20, F20.10)') "Sum of Born charges:", born_sum
238 1 : WRITE (log_unit, '(/,A30, F20.10)') "Checksum (Acoustic Sum Rule):", born_sum
239 : END IF
240 2 : DEALLOCATE (ap_tensor)
241 2 : DEALLOCATE (dft_control%period_efield)
242 :
243 2 : CALL cp_print_key_finished_output(apt_log, logger, dcdr_section, "PRINT%APT")
244 :
245 2 : dft_control%apply_period_efield = .FALSE.
246 2 : qs_env%linres_run = .FALSE.
247 :
248 2 : IF (log_unit > 0) THEN
249 1 : WRITE (log_unit, '(T2,A)') " "
250 1 : WRITE (log_unit, '(T2,A)') " "
251 1 : WRITE (log_unit, '(T22,A)') "APT calculation Done!"
252 : WRITE (log_unit, '(/,T2,A)') &
253 1 : '!----------------------------------------------------------------------------!'
254 : END IF
255 :
256 2 : CALL timestop(handle)
257 :
258 18 : END SUBROUTINE apt_fdiff
259 :
260 : END MODULE qs_apt_fdiff_methods
|
