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