Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Quickstep force driver routine
10 : !> \author MK (12.06.2002)
11 : ! **************************************************************************************************
12 : MODULE qs_force
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind_set
15 : USE cell_types, ONLY: cell_type
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
18 : dbcsr_p_type,&
19 : dbcsr_set
20 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
21 : dbcsr_deallocate_matrix_set
22 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_get_default_io_unit,&
25 : cp_logger_type
26 : USE cp_output_handling, ONLY: cp_p_file,&
27 : cp_print_key_finished_output,&
28 : cp_print_key_should_output,&
29 : cp_print_key_unit_nr
30 : USE dft_plus_u, ONLY: plus_u
31 : USE ec_env_types, ONLY: energy_correction_type
32 : USE efield_utils, ONLY: calculate_ecore_efield,&
33 : efield_potential_lengh_gauge
34 : USE energy_corrections, ONLY: energy_correction
35 : USE excited_states, ONLY: excited_state_energy
36 : USE hfx_exx, ONLY: calculate_exx
37 : USE input_constants, ONLY: &
38 : do_method_gapw, do_method_gapw_xc, do_method_gpw, do_method_lrigpw, do_method_ofgpw, &
39 : do_method_rigpw, plus_u_lowdin, ri_mp2_laplace, ri_mp2_method_gpw, ri_rpa_method_gpw, &
40 : smeagol_runtype_emtransport
41 : USE input_section_types, ONLY: section_vals_get,&
42 : section_vals_get_subs_vals,&
43 : section_vals_type,&
44 : section_vals_val_get
45 : USE kinds, ONLY: dp
46 : USE lri_environment_types, ONLY: lri_environment_type
47 : USE message_passing, ONLY: mp_para_env_type
48 : USE mp2_cphf, ONLY: update_mp2_forces
49 : USE mulliken, ONLY: mulliken_restraint
50 : USE particle_types, ONLY: particle_type
51 : USE qs_core_energies, ONLY: calculate_ecore_overlap,&
52 : calculate_ecore_self
53 : USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix
54 : USE qs_dftb_dispersion, ONLY: calculate_dftb_dispersion
55 : USE qs_dftb_matrices, ONLY: build_dftb_matrices
56 : USE qs_energy, ONLY: qs_energies
57 : USE qs_energy_types, ONLY: qs_energy_type
58 : USE qs_environment_methods, ONLY: qs_env_rebuild_pw_env
59 : USE qs_environment_types, ONLY: get_qs_env,&
60 : qs_environment_type
61 : USE qs_external_potential, ONLY: external_c_potential,&
62 : external_e_potential
63 : USE qs_force_types, ONLY: allocate_qs_force,&
64 : qs_force_type,&
65 : replicate_qs_force,&
66 : zero_qs_force
67 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
68 : USE qs_ks_types, ONLY: qs_ks_env_type,&
69 : set_ks_env
70 : USE qs_rho_types, ONLY: qs_rho_get,&
71 : qs_rho_type
72 : USE qs_scf_post_scf, ONLY: qs_scf_compute_properties
73 : USE qs_subsys_types, ONLY: qs_subsys_set,&
74 : qs_subsys_type
75 : USE ri_environment_methods, ONLY: build_ri_matrices
76 : USE rt_propagation_forces, ONLY: calc_c_mat_force,&
77 : rt_admm_force
78 : USE rt_propagation_velocity_gauge, ONLY: velocity_gauge_ks_matrix,&
79 : velocity_gauge_nl_force
80 : USE se_core_core, ONLY: se_core_core_interaction
81 : USE se_core_matrix, ONLY: build_se_core_matrix
82 : USE tblite_interface, ONLY: build_tblite_matrices
83 : USE tblite_reference_cli, ONLY: tb_reference_cli_compare
84 : USE virial_types, ONLY: project_virial_to_periodic_subspace,&
85 : symmetrize_virial,&
86 : virial_type
87 : USE xtb_matrices, ONLY: build_xtb_matrices
88 : #include "./base/base_uses.f90"
89 :
90 : IMPLICIT NONE
91 :
92 : PRIVATE
93 :
94 : ! *** Global parameters ***
95 :
96 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_force'
97 :
98 : ! *** Public subroutines ***
99 :
100 : PUBLIC :: qs_calc_energy_force
101 :
102 : CONTAINS
103 :
104 : ! **************************************************************************************************
105 : !> \brief ...
106 : !> \param qs_env ...
107 : !> \param calc_force ...
108 : !> \param consistent_energies ...
109 : !> \param linres ...
110 : ! **************************************************************************************************
111 30008 : SUBROUTINE qs_calc_energy_force(qs_env, calc_force, consistent_energies, linres)
112 : TYPE(qs_environment_type), POINTER :: qs_env
113 : LOGICAL :: calc_force, consistent_energies, linres
114 :
115 30008 : qs_env%linres_run = linres
116 30008 : IF (calc_force) THEN
117 11439 : CALL qs_forces(qs_env)
118 : ELSE
119 : CALL qs_energies(qs_env, calc_forces=.FALSE., &
120 18569 : consistent_energies=consistent_energies)
121 : END IF
122 :
123 30007 : END SUBROUTINE qs_calc_energy_force
124 :
125 : ! **************************************************************************************************
126 : !> \brief Calculate the Quickstep forces.
127 : !> \param qs_env ...
128 : !> \date 29.10.2002
129 : !> \author MK
130 : !> \version 1.0
131 : ! **************************************************************************************************
132 11439 : SUBROUTINE qs_forces(qs_env)
133 :
134 : TYPE(qs_environment_type), POINTER :: qs_env
135 :
136 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_forces'
137 :
138 : INTEGER :: after, handle, i, iatom, ic, ikind, &
139 : ispin, iw, natom, nkind, nspin, &
140 : output_unit
141 11439 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind
142 : LOGICAL :: do_admm, do_exx, do_gw, do_im_time, &
143 : do_kpoints, has_unit_metric, &
144 : omit_headers, perform_ec, reuse_hfx
145 : REAL(dp) :: dummy_real, dummy_real2(2)
146 11439 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
147 : TYPE(cell_type), POINTER :: cell
148 : TYPE(cp_logger_type), POINTER :: logger
149 11439 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_w, rho_ao
150 11439 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_w_kp
151 : TYPE(dft_control_type), POINTER :: dft_control
152 : TYPE(energy_correction_type), POINTER :: ec_env
153 : TYPE(lri_environment_type), POINTER :: lri_env
154 : TYPE(mp_para_env_type), POINTER :: para_env
155 11439 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
156 : TYPE(qs_energy_type), POINTER :: energy
157 11439 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
158 : TYPE(qs_ks_env_type), POINTER :: ks_env
159 : TYPE(qs_rho_type), POINTER :: rho
160 : TYPE(qs_subsys_type), POINTER :: subsys
161 : TYPE(section_vals_type), POINTER :: hfx_sections, print_section
162 : TYPE(virial_type), POINTER :: virial
163 :
164 11439 : CALL timeset(routineN, handle)
165 11439 : NULLIFY (logger)
166 11439 : logger => cp_get_default_logger()
167 :
168 : ! rebuild plane wave environment
169 11439 : CALL qs_env_rebuild_pw_env(qs_env)
170 :
171 : ! zero out the forces in particle set
172 11439 : CALL get_qs_env(qs_env, particle_set=particle_set)
173 11439 : natom = SIZE(particle_set)
174 82488 : DO iatom = 1, natom
175 295635 : particle_set(iatom)%f = 0.0_dp
176 : END DO
177 :
178 : ! get atom mapping
179 11439 : NULLIFY (atomic_kind_set)
180 11439 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
181 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
182 : atom_of_kind=atom_of_kind, &
183 11439 : kind_of=kind_of)
184 :
185 11439 : NULLIFY (force, subsys, dft_control)
186 : CALL get_qs_env(qs_env, &
187 : force=force, &
188 : subsys=subsys, &
189 : do_kpoints=do_kpoints, &
190 11439 : dft_control=dft_control)
191 11439 : IF (do_kpoints .AND. dft_control%dft_plus_u .AND. dft_control%plus_u_method_id == plus_u_lowdin) THEN
192 8 : IF (ASSOCIATED(qs_env%mp2_env) .OR. qs_env%run_rtp .OR. &
193 : (dft_control%smeagol_control%smeagol_enabled .AND. &
194 : dft_control%smeagol_control%run_type == smeagol_runtype_emtransport)) THEN
195 0 : CPABORT("Lowdin DFT+U forces with k-points require the ground-state DFT W matrix")
196 : END IF
197 : END IF
198 11439 : IF (.NOT. ASSOCIATED(force)) THEN
199 : ! *** Allocate the force data structure ***
200 3589 : nkind = SIZE(atomic_kind_set)
201 3589 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
202 3589 : CALL allocate_qs_force(force, natom_of_kind)
203 3589 : DEALLOCATE (natom_of_kind)
204 3589 : CALL qs_subsys_set(subsys, force=force)
205 : END IF
206 11439 : CALL zero_qs_force(force)
207 :
208 : ! Check if CDFT potential is needed and save it until forces have been calculated
209 11439 : IF (dft_control%qs_control%cdft) THEN
210 126 : dft_control%qs_control%cdft_control%save_pot = .TRUE.
211 : END IF
212 :
213 : ! recalculate energy and the response vector for the Z-vector linear equation system if calc_force = .true.
214 11439 : CALL qs_energies(qs_env, calc_forces=.TRUE.)
215 :
216 11439 : NULLIFY (para_env)
217 : CALL get_qs_env(qs_env, &
218 11439 : para_env=para_env)
219 :
220 : ! Now we handle some special cases
221 : ! Maybe some of these would be better dealt with in qs_energies?
222 11439 : IF (qs_env%run_rtp) THEN
223 1228 : NULLIFY (matrix_w, matrix_s, ks_env)
224 : CALL get_qs_env(qs_env, &
225 : ks_env=ks_env, &
226 : matrix_w=matrix_w, &
227 1228 : matrix_s=matrix_s)
228 1228 : CALL dbcsr_allocate_matrix_set(matrix_w, dft_control%nspins)
229 2708 : DO ispin = 1, dft_control%nspins
230 1480 : ALLOCATE (matrix_w(ispin)%matrix)
231 : CALL dbcsr_copy(matrix_w(ispin)%matrix, matrix_s(1)%matrix, &
232 1480 : name="W MATRIX")
233 2708 : CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
234 : END DO
235 1228 : CALL set_ks_env(ks_env, matrix_w=matrix_w)
236 :
237 1228 : CALL calc_c_mat_force(qs_env)
238 1228 : IF (dft_control%do_admm) CALL rt_admm_force(qs_env)
239 1228 : IF (dft_control%rtp_control%velocity_gauge .AND. dft_control%rtp_control%nl_gauge_transform) THEN
240 32 : CALL velocity_gauge_nl_force(qs_env, particle_set)
241 : END IF
242 : END IF
243 : ! from an eventual Mulliken restraint
244 11439 : IF (dft_control%qs_control%mulliken_restraint) THEN
245 6 : NULLIFY (matrix_w, matrix_s, rho)
246 : CALL get_qs_env(qs_env, &
247 : matrix_w=matrix_w, &
248 : matrix_s=matrix_s, &
249 6 : rho=rho)
250 6 : NULLIFY (rho_ao)
251 6 : CALL qs_rho_get(rho, rho_ao=rho_ao)
252 : CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
253 6 : para_env, matrix_s(1)%matrix, rho_ao, w_matrix=matrix_w)
254 : END IF
255 : ! Add non-Pulay contribution of DFT+U to W matrix, since it has also to be
256 : ! digested with overlap matrix derivatives
257 : ! Lowdin k-point derivatives were already added before transforming W(k).
258 11439 : IF (dft_control%dft_plus_u .AND. .NOT. (do_kpoints .AND. dft_control%plus_u_method_id == plus_u_lowdin)) THEN
259 92 : NULLIFY (matrix_w_kp)
260 92 : CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp)
261 92 : CALL plus_u(qs_env=qs_env, matrix_w=matrix_w_kp)
262 : END IF
263 :
264 : ! Write W Matrix to output (if requested)
265 11439 : CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
266 11439 : IF (.NOT. has_unit_metric) THEN
267 8415 : NULLIFY (matrix_w_kp)
268 8415 : CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp)
269 8415 : nspin = SIZE(matrix_w_kp, 1)
270 17926 : DO ispin = 1, nspin
271 9511 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
272 8415 : qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX"), cp_p_file)) THEN
273 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX", &
274 8 : extension=".Log")
275 8 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
276 8 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
277 8 : after = MIN(MAX(after, 1), 16)
278 16 : DO ic = 1, SIZE(matrix_w_kp, 2)
279 : CALL cp_dbcsr_write_sparse_matrix(matrix_w_kp(ispin, ic)%matrix, 4, after, qs_env, &
280 16 : para_env, output_unit=iw, omit_headers=omit_headers)
281 : END DO
282 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
283 8 : "DFT%PRINT%AO_MATRICES/W_MATRIX")
284 : END IF
285 : END DO
286 : END IF
287 :
288 : ! Check if energy correction should be skipped
289 11439 : perform_ec = .FALSE.
290 11439 : IF (qs_env%energy_correction) THEN
291 498 : CALL get_qs_env(qs_env, ec_env=ec_env)
292 498 : IF (.NOT. ec_env%do_skip) perform_ec = .TRUE.
293 : END IF
294 :
295 : ! Compute core forces (also overwrites matrix_w)
296 11439 : IF (dft_control%qs_control%semi_empirical) THEN
297 : CALL build_se_core_matrix(qs_env=qs_env, para_env=para_env, &
298 3024 : calculate_forces=.TRUE.)
299 3024 : CALL se_core_core_interaction(qs_env, para_env, calculate_forces=.TRUE.)
300 8415 : ELSE IF (dft_control%qs_control%dftb) THEN
301 : CALL build_dftb_matrices(qs_env=qs_env, para_env=para_env, &
302 788 : calculate_forces=.TRUE.)
303 : CALL calculate_dftb_dispersion(qs_env=qs_env, para_env=para_env, &
304 788 : calculate_forces=.TRUE.)
305 7627 : ELSE IF (dft_control%qs_control%xtb) THEN
306 782 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
307 154 : CALL build_tblite_matrices(qs_env=qs_env, calculate_forces=.TRUE.)
308 : ELSE
309 628 : CALL build_xtb_matrices(qs_env=qs_env, calculate_forces=.TRUE.)
310 : END IF
311 6845 : ELSE IF (perform_ec) THEN
312 : ! Calculates core and grid based forces
313 498 : CALL energy_correction(qs_env, ec_init=.FALSE., calculate_forces=.TRUE.)
314 : ELSE
315 : ! Dispersion energy and forces are calculated in qs_energy?
316 6347 : CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.TRUE.)
317 : ! The above line reset the core H, which should be re-updated in case a TD field is applied:
318 6347 : IF (qs_env%run_rtp) THEN
319 824 : IF (dft_control%apply_efield_field) THEN
320 160 : CALL efield_potential_lengh_gauge(qs_env)
321 : END IF
322 824 : IF (dft_control%rtp_control%velocity_gauge) THEN
323 32 : CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.FALSE.)
324 : END IF
325 :
326 : END IF
327 6347 : CALL calculate_ecore_self(qs_env)
328 6347 : CALL calculate_ecore_overlap(qs_env, para_env, calculate_forces=.TRUE.)
329 6347 : CALL calculate_ecore_efield(qs_env, calculate_forces=.TRUE.)
330 : !swap external_e_potential before external_c_potential, to ensure
331 : !that external potential on grid is loaded before calculating energy of cores
332 6347 : CALL external_e_potential(qs_env)
333 6347 : IF (.NOT. dft_control%qs_control%gapw) THEN
334 5681 : CALL external_c_potential(qs_env, calculate_forces=.TRUE.)
335 : END IF
336 : ! RIGPW matrices
337 6347 : IF (dft_control%qs_control%rigpw) THEN
338 2 : CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
339 2 : CALL build_ri_matrices(lri_env, qs_env, calculate_forces=.TRUE.)
340 : END IF
341 : END IF
342 :
343 : ! MP2 Code
344 11439 : IF (ASSOCIATED(qs_env%mp2_env)) THEN
345 322 : NULLIFY (energy)
346 322 : CALL get_qs_env(qs_env, energy=energy)
347 322 : CALL qs_scf_compute_properties(qs_env, wf_type='MP2 ', do_mp2=.TRUE.)
348 322 : CALL qs_ks_update_qs_env(qs_env, just_energy=.TRUE.)
349 322 : energy%total = energy%total + energy%mp2
350 :
351 : IF ((qs_env%mp2_env%method == ri_mp2_method_gpw .OR. qs_env%mp2_env%method == ri_mp2_laplace .OR. &
352 : qs_env%mp2_env%method == ri_rpa_method_gpw) &
353 322 : .AND. .NOT. qs_env%mp2_env%do_im_time) THEN
354 272 : CALL update_mp2_forces(qs_env)
355 : END IF
356 :
357 : !RPA EXX energy and forces
358 322 : IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
359 : do_exx = .FALSE.
360 52 : hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
361 52 : CALL section_vals_get(hfx_sections, explicit=do_exx)
362 52 : IF (do_exx) THEN
363 26 : do_gw = qs_env%mp2_env%ri_rpa%do_ri_g0w0
364 26 : do_admm = qs_env%mp2_env%ri_rpa%do_admm
365 26 : reuse_hfx = qs_env%mp2_env%ri_rpa%reuse_hfx
366 26 : do_im_time = qs_env%mp2_env%do_im_time
367 26 : output_unit = cp_logger_get_default_io_unit()
368 26 : dummy_real = 0.0_dp
369 :
370 : CALL calculate_exx(qs_env=qs_env, &
371 : unit_nr=output_unit, &
372 : hfx_sections=hfx_sections, &
373 : x_data=qs_env%mp2_env%ri_rpa%x_data, &
374 : do_gw=do_gw, &
375 : do_admm=do_admm, &
376 : calc_forces=.TRUE., &
377 : reuse_hfx=reuse_hfx, &
378 : do_im_time=do_im_time, &
379 : E_ex_from_GW=dummy_real, &
380 : E_admm_from_GW=dummy_real2, &
381 26 : t3=dummy_real)
382 : END IF
383 : END IF
384 11117 : ELSE IF (perform_ec) THEN
385 : ! energy correction forces postponed
386 10619 : ELSE IF (qs_env%harris_method) THEN
387 : ! Harris method forces already done in harris_energy_correction
388 : ELSE
389 : ! Compute grid-based forces
390 10613 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.TRUE.)
391 : END IF
392 :
393 : ! Excited state forces
394 : ! Solve the response linear equation system for the Z-vector method
395 : ! and calculate remaining terms of the force
396 11439 : CALL excited_state_energy(qs_env, calculate_forces=.TRUE.)
397 :
398 : ! replicate forces (get current pointer)
399 11439 : NULLIFY (force)
400 11439 : CALL get_qs_env(qs_env=qs_env, force=force)
401 11439 : CALL replicate_qs_force(force, para_env)
402 :
403 82488 : DO iatom = 1, natom
404 71049 : ikind = kind_of(iatom)
405 71049 : i = atom_of_kind(iatom)
406 : ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
407 : ! the force is - dE/dR, what is called force is actually the gradient
408 : ! Things should have the right name
409 : ! The minus sign below is a hack
410 : ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
411 568392 : force(ikind)%other(1:3, i) = -particle_set(iatom)%f(1:3) + force(ikind)%ch_pulay(1:3, i)
412 284196 : force(ikind)%total(1:3, i) = force(ikind)%total(1:3, i) + force(ikind)%other(1:3, i)
413 579831 : particle_set(iatom)%f = -force(ikind)%total(1:3, i)
414 : END DO
415 :
416 11439 : NULLIFY (cell, virial, energy)
417 11439 : CALL get_qs_env(qs_env=qs_env, cell=cell, virial=virial, energy=energy)
418 11439 : IF (virial%pv_availability) THEN
419 1264 : CALL para_env%sum(virial%pv_overlap)
420 1264 : CALL para_env%sum(virial%pv_ekinetic)
421 1264 : CALL para_env%sum(virial%pv_ppl)
422 1264 : CALL para_env%sum(virial%pv_ppnl)
423 1264 : CALL para_env%sum(virial%pv_ecore_overlap)
424 1264 : CALL para_env%sum(virial%pv_ehartree)
425 1264 : CALL para_env%sum(virial%pv_exc)
426 1264 : CALL para_env%sum(virial%pv_exx)
427 1264 : CALL para_env%sum(virial%pv_vdw)
428 1264 : CALL para_env%sum(virial%pv_mp2)
429 1264 : CALL para_env%sum(virial%pv_nlcc)
430 1264 : CALL para_env%sum(virial%pv_gapw)
431 1264 : CALL para_env%sum(virial%pv_lrigpw)
432 1264 : CALL para_env%sum(virial%pv_virial)
433 1264 : CALL symmetrize_virial(virial)
434 : ! Add the volume terms of the virial
435 1264 : IF ((.NOT. virial%pv_numer) .AND. &
436 : (.NOT. (dft_control%qs_control%dftb .OR. &
437 : dft_control%qs_control%xtb .OR. &
438 : dft_control%qs_control%semi_empirical))) THEN
439 :
440 : ! Harris energy correction requires volume terms from
441 : ! 1) Harris functional contribution, and
442 : ! 2) Linear Response solver
443 762 : IF (perform_ec) THEN
444 168 : CALL get_qs_env(qs_env, ec_env=ec_env)
445 168 : energy%hartree = ec_env%ehartree
446 168 : energy%exc = ec_env%exc
447 168 : IF (dft_control%do_admm) THEN
448 38 : energy%exc_aux_fit = ec_env%exc_aux_fit
449 : END IF
450 : END IF
451 3048 : DO i = 1, 3
452 : virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) &
453 2286 : - 2.0_dp*(energy%hartree + energy%sccs_pol)
454 : virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc &
455 2286 : - 2.0_dp*(energy%hartree + energy%sccs_pol)
456 2286 : virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc
457 3048 : IF (dft_control%do_admm) THEN
458 222 : virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc_aux_fit
459 222 : virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc_aux_fit
460 : END IF
461 : ! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve.
462 : ! The sign in pw_poisson_solve is correct for FIST, but not for QS.
463 : ! There should be a more elegant solution to that ...
464 : END DO
465 : END IF
466 5056 : IF ((.NOT. virial%pv_numer) .AND. COUNT(cell%perd /= 0) == 2) THEN
467 48 : SELECT CASE (dft_control%qs_control%method_id)
468 : CASE (do_method_gapw, do_method_gapw_xc, do_method_gpw, &
469 : do_method_lrigpw, do_method_ofgpw, do_method_rigpw)
470 32 : CALL project_virial_to_periodic_subspace(virial, cell%perd)
471 : END SELECT
472 : END IF
473 : END IF
474 :
475 11439 : IF (dft_control%qs_control%xtb .AND. dft_control%qs_control%xtb_control%do_tblite) THEN
476 154 : CALL tb_reference_cli_compare(qs_env)
477 : END IF
478 :
479 : output_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%DERIVATIVES", &
480 11439 : extension=".Log")
481 11439 : print_section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%DERIVATIVES")
482 11439 : IF (dft_control%qs_control%semi_empirical) THEN
483 : CALL write_forces(force, atomic_kind_set, 2, output_unit=output_unit, &
484 3024 : print_section=print_section)
485 8415 : ELSE IF (dft_control%qs_control%dftb) THEN
486 : CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, &
487 788 : print_section=print_section)
488 7627 : ELSE IF (dft_control%qs_control%xtb) THEN
489 : CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, &
490 782 : print_section=print_section)
491 6845 : ELSE IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
492 : CALL write_forces(force, atomic_kind_set, 1, output_unit=output_unit, &
493 840 : print_section=print_section)
494 : ELSE
495 : CALL write_forces(force, atomic_kind_set, 0, output_unit=output_unit, &
496 6005 : print_section=print_section)
497 : END IF
498 : CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, &
499 11439 : "DFT%PRINT%DERIVATIVES")
500 :
501 : ! deallocate W Matrix:
502 11439 : NULLIFY (ks_env, matrix_w_kp)
503 : CALL get_qs_env(qs_env=qs_env, &
504 : matrix_w_kp=matrix_w_kp, &
505 11439 : ks_env=ks_env)
506 11439 : CALL dbcsr_deallocate_matrix_set(matrix_w_kp)
507 11439 : NULLIFY (matrix_w_kp)
508 11439 : CALL set_ks_env(ks_env, matrix_w_kp=matrix_w_kp)
509 :
510 11439 : DEALLOCATE (atom_of_kind, kind_of)
511 :
512 11439 : CALL timestop(handle)
513 :
514 22878 : END SUBROUTINE qs_forces
515 :
516 : ! **************************************************************************************************
517 : !> \brief Write a Quickstep force data structure to output unit
518 : !> \param qs_force ...
519 : !> \param atomic_kind_set ...
520 : !> \param ftype ...
521 : !> \param output_unit ...
522 : !> \param print_section ...
523 : !> \date 05.06.2002
524 : !> \author MK
525 : !> \version 1.0
526 : ! **************************************************************************************************
527 11439 : SUBROUTINE write_forces(qs_force, atomic_kind_set, ftype, output_unit, &
528 : print_section)
529 :
530 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
531 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
532 : INTEGER, INTENT(IN) :: ftype, output_unit
533 : TYPE(section_vals_type), POINTER :: print_section
534 :
535 : CHARACTER(LEN=13) :: fmtstr5
536 : CHARACTER(LEN=15) :: fmtstr4
537 : CHARACTER(LEN=20) :: fmtstr3
538 : CHARACTER(LEN=35) :: fmtstr2
539 : CHARACTER(LEN=48) :: fmtstr1
540 : INTEGER :: i, iatom, ikind, my_ftype, natom, ndigits
541 11439 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
542 : REAL(KIND=dp), DIMENSION(3) :: grand_total
543 :
544 11439 : IF (output_unit > 0) THEN
545 :
546 180 : IF (.NOT. ASSOCIATED(qs_force)) THEN
547 : CALL cp_abort(__LOCATION__, &
548 : "The qs_force pointer is not associated "// &
549 0 : "and cannot be printed")
550 : END IF
551 :
552 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, &
553 180 : kind_of=kind_of, natom=natom)
554 :
555 : ! Variable precision output of the forces
556 : CALL section_vals_val_get(print_section, "NDIGITS", &
557 180 : i_val=ndigits)
558 :
559 180 : fmtstr1 = "(/,/,T2,A,/,/,T3,A,T11,A,T23,A,T40,A1,2( X,A1))"
560 180 : WRITE (UNIT=fmtstr1(41:42), FMT="(I2)") ndigits + 5
561 :
562 180 : fmtstr2 = "(/,(T2,I5,4X,I4,T18,A,T34,3F . ))"
563 180 : WRITE (UNIT=fmtstr2(32:33), FMT="(I2)") ndigits
564 180 : WRITE (UNIT=fmtstr2(29:30), FMT="(I2)") ndigits + 6
565 :
566 180 : fmtstr3 = "(/,T3,A,T34,3F . )"
567 180 : WRITE (UNIT=fmtstr3(18:19), FMT="(I2)") ndigits
568 180 : WRITE (UNIT=fmtstr3(15:16), FMT="(I2)") ndigits + 6
569 :
570 180 : fmtstr4 = "((T34,3F . ))"
571 180 : WRITE (UNIT=fmtstr4(12:13), FMT="(I2)") ndigits
572 180 : WRITE (UNIT=fmtstr4(9:10), FMT="(I2)") ndigits + 6
573 :
574 : fmtstr5 = "(/T2,A//T3,A)"
575 :
576 : WRITE (UNIT=output_unit, FMT=fmtstr1) &
577 180 : "FORCES [a.u.]", "Atom", "Kind", "Component", "X", "Y", "Z"
578 :
579 180 : grand_total(:) = 0.0_dp
580 :
581 180 : my_ftype = ftype
582 :
583 0 : SELECT CASE (my_ftype)
584 : CASE DEFAULT
585 0 : DO iatom = 1, natom
586 0 : ikind = kind_of(iatom)
587 0 : i = atom_of_kind(iatom)
588 : WRITE (UNIT=output_unit, FMT=fmtstr2) &
589 0 : iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
590 0 : grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
591 : END DO
592 : CASE (0)
593 476 : DO iatom = 1, natom
594 342 : ikind = kind_of(iatom)
595 342 : i = atom_of_kind(iatom)
596 : WRITE (UNIT=output_unit, FMT=fmtstr2) &
597 1368 : iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
598 1368 : iatom, ikind, " overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
599 1368 : iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
600 1368 : iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
601 1368 : iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
602 1368 : iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
603 1368 : iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
604 1368 : iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
605 1368 : iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
606 1368 : iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
607 1368 : iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
608 1368 : iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
609 1368 : iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), &
610 1368 : iatom, ikind, " other", qs_force(ikind)%other(1:3, i), &
611 1368 : iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
612 1368 : iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
613 1368 : iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
614 1368 : iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), &
615 1368 : iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
616 1368 : iatom, ikind, " tensorial_u", qs_force(ikind)%tensorial_u(1:3, i), &
617 1710 : iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
618 1502 : grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
619 : END DO
620 : CASE (1)
621 72 : DO iatom = 1, natom
622 52 : ikind = kind_of(iatom)
623 52 : i = atom_of_kind(iatom)
624 : WRITE (UNIT=output_unit, FMT=fmtstr2) &
625 208 : iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
626 208 : iatom, ikind, " overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
627 208 : iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
628 208 : iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
629 208 : iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
630 208 : iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
631 208 : iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
632 208 : iatom, ikind, "cneo_potential", qs_force(ikind)%cneo_potential(1:3, i), &
633 208 : iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
634 208 : iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
635 208 : iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
636 208 : iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
637 208 : iatom, ikind, " rho_cneo_nuc", qs_force(ikind)%rho_cneo_nuc(1:3, i), &
638 208 : iatom, ikind, " vhxc_atom", qs_force(ikind)%vhxc_atom(1:3, i), &
639 208 : iatom, ikind, " g0s_Vh_elec", qs_force(ikind)%g0s_Vh_elec(1:3, i), &
640 208 : iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
641 208 : iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
642 208 : iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), &
643 208 : iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
644 208 : iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
645 208 : iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
646 208 : iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), &
647 208 : iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
648 208 : iatom, ikind, " tensorial_u", qs_force(ikind)%tensorial_u(1:3, i), &
649 260 : iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
650 228 : grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
651 : END DO
652 : CASE (2)
653 75 : DO iatom = 1, natom
654 73 : ikind = kind_of(iatom)
655 73 : i = atom_of_kind(iatom)
656 : WRITE (UNIT=output_unit, FMT=fmtstr2) &
657 292 : iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
658 292 : iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
659 365 : iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
660 294 : grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
661 : END DO
662 : CASE (3)
663 0 : DO iatom = 1, natom
664 0 : ikind = kind_of(iatom)
665 0 : i = atom_of_kind(iatom)
666 : WRITE (UNIT=output_unit, FMT=fmtstr2) &
667 0 : iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
668 0 : iatom, ikind, "overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
669 0 : iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
670 0 : iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
671 0 : iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
672 0 : iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
673 0 : iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
674 0 : iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
675 0 : iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
676 0 : iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
677 0 : iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
678 0 : iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
679 0 : iatom, ikind, " tensorial_u", qs_force(ikind)%tensorial_u(1:3, i), &
680 0 : iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
681 0 : grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
682 : END DO
683 : CASE (4)
684 188 : DO iatom = 1, natom
685 164 : ikind = kind_of(iatom)
686 164 : i = atom_of_kind(iatom)
687 : WRITE (UNIT=output_unit, FMT=fmtstr2) &
688 656 : iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
689 656 : iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
690 656 : iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
691 656 : iatom, ikind, " repulsive", qs_force(ikind)%repulsive(1:3, i), &
692 656 : iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
693 656 : iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
694 656 : iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
695 820 : iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
696 680 : grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
697 : END DO
698 : CASE (5)
699 180 : DO iatom = 1, natom
700 0 : ikind = kind_of(iatom)
701 0 : i = atom_of_kind(iatom)
702 : WRITE (UNIT=output_unit, FMT=fmtstr2) &
703 0 : iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
704 0 : iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
705 0 : iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
706 0 : iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
707 0 : iatom, ikind, " all potential", qs_force(ikind)%all_potential(1:3, i), &
708 0 : iatom, ikind, " other", qs_force(ikind)%other(1:3, i), &
709 0 : iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
710 0 : grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
711 : END DO
712 : END SELECT
713 :
714 180 : WRITE (UNIT=output_unit, FMT=fmtstr3) "Sum of total", grand_total(1:3)
715 :
716 180 : DEALLOCATE (atom_of_kind)
717 180 : DEALLOCATE (kind_of)
718 :
719 : END IF
720 :
721 11439 : END SUBROUTINE write_forces
722 :
723 : END MODULE qs_force
|
