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