(git:b195825)
qs_force.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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,&
15  USE cp_control_types, ONLY: dft_control_type
21  cp_logger_type
22  USE cp_output_handling, ONLY: cp_p_file,&
26  USE dbcsr_api, ONLY: dbcsr_copy,&
27  dbcsr_p_type,&
28  dbcsr_set
29  USE dft_plus_u, ONLY: plus_u
30  USE ec_env_types, ONLY: energy_correction_type
35  USE hfx_exx, ONLY: calculate_exx
36  USE input_constants, ONLY: ri_mp2_laplace,&
41  section_vals_type,&
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
54  USE qs_energy, ONLY: qs_energies
55  USE qs_energy_types, ONLY: qs_energy_type
57  USE qs_environment_types, ONLY: get_qs_env,&
58  qs_environment_type
62  qs_force_type,&
66  USE qs_ks_types, ONLY: qs_ks_env_type,&
68  USE qs_rho_types, ONLY: qs_rho_get,&
69  qs_rho_type
71  USE qs_subsys_types, ONLY: qs_subsys_set,&
72  qs_subsys_type
80  USE virial_types, ONLY: symmetrize_virial,&
81  virial_type
83 #include "./base/base_uses.f90"
84 
85  IMPLICIT NONE
86 
87  PRIVATE
88 
89 ! *** Global parameters ***
90 
91  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_force'
92 
93 ! *** Public subroutines ***
94 
95  PUBLIC :: qs_calc_energy_force
96 
97 CONTAINS
98 
99 ! **************************************************************************************************
100 !> \brief ...
101 !> \param qs_env ...
102 !> \param calc_force ...
103 !> \param consistent_energies ...
104 !> \param linres ...
105 ! **************************************************************************************************
106  SUBROUTINE qs_calc_energy_force(qs_env, calc_force, consistent_energies, linres)
107  TYPE(qs_environment_type), POINTER :: qs_env
108  LOGICAL :: calc_force, consistent_energies, linres
109 
110  qs_env%linres_run = linres
111  IF (calc_force) THEN
112  CALL qs_forces(qs_env)
113  ELSE
114  CALL qs_energies(qs_env, calc_forces=.false., &
115  consistent_energies=consistent_energies)
116  END IF
117 
118  END SUBROUTINE qs_calc_energy_force
119 
120 ! **************************************************************************************************
121 !> \brief Calculate the Quickstep forces.
122 !> \param qs_env ...
123 !> \date 29.10.2002
124 !> \author MK
125 !> \version 1.0
126 ! **************************************************************************************************
127  SUBROUTINE qs_forces(qs_env)
128 
129  TYPE(qs_environment_type), POINTER :: qs_env
130 
131  CHARACTER(len=*), PARAMETER :: routinen = 'qs_forces'
132 
133  INTEGER :: after, handle, i, iatom, ic, ikind, &
134  ispin, iw, natom, nkind, nspin, &
135  output_unit
136  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind
137  LOGICAL :: do_admm, do_exx, do_gw, do_im_time, &
138  has_unit_metric, omit_headers, &
139  perform_ec, reuse_hfx
140  REAL(dp) :: dummy_real, dummy_real2(2)
141  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
142  TYPE(cp_logger_type), POINTER :: logger
143  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_w, rho_ao
144  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_w_kp
145  TYPE(dft_control_type), POINTER :: dft_control
146  TYPE(energy_correction_type), POINTER :: ec_env
147  TYPE(lri_environment_type), POINTER :: lri_env
148  TYPE(mp_para_env_type), POINTER :: para_env
149  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
150  TYPE(qs_energy_type), POINTER :: energy
151  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
152  TYPE(qs_ks_env_type), POINTER :: ks_env
153  TYPE(qs_rho_type), POINTER :: rho
154  TYPE(qs_subsys_type), POINTER :: subsys
155  TYPE(section_vals_type), POINTER :: hfx_sections, print_section
156  TYPE(virial_type), POINTER :: virial
157 
158  CALL timeset(routinen, handle)
159  NULLIFY (logger)
160  logger => cp_get_default_logger()
161 
162  ! rebuild plane wave environment
163  CALL qs_env_rebuild_pw_env(qs_env)
164 
165  ! zero out the forces in particle set
166  CALL get_qs_env(qs_env, particle_set=particle_set)
167  natom = SIZE(particle_set)
168  DO iatom = 1, natom
169  particle_set(iatom)%f = 0.0_dp
170  END DO
171 
172  ! get atom mapping
173  NULLIFY (atomic_kind_set)
174  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
175  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
176  atom_of_kind=atom_of_kind, &
177  kind_of=kind_of)
178 
179  NULLIFY (force, subsys, dft_control)
180  CALL get_qs_env(qs_env, &
181  force=force, &
182  subsys=subsys, &
183  dft_control=dft_control)
184  IF (.NOT. ASSOCIATED(force)) THEN
185  ! *** Allocate the force data structure ***
186  nkind = SIZE(atomic_kind_set)
187  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
188  CALL allocate_qs_force(force, natom_of_kind)
189  DEALLOCATE (natom_of_kind)
190  CALL qs_subsys_set(subsys, force=force)
191  END IF
192  CALL zero_qs_force(force)
193 
194  ! Check if CDFT potential is needed and save it until forces have been calculated
195  IF (dft_control%qs_control%cdft) THEN
196  dft_control%qs_control%cdft_control%save_pot = .true.
197  END IF
198 
199  ! recalculate energy with forces
200  CALL qs_energies(qs_env, calc_forces=.true.)
201 
202  NULLIFY (para_env)
203  CALL get_qs_env(qs_env, &
204  para_env=para_env)
205 
206  ! Now we handle some special cases
207  ! Maybe some of these would be better dealt with in qs_energies?
208  IF (qs_env%run_rtp) THEN
209  NULLIFY (matrix_w, matrix_s, ks_env)
210  CALL get_qs_env(qs_env, &
211  ks_env=ks_env, &
212  matrix_w=matrix_w, &
213  matrix_s=matrix_s)
214  CALL dbcsr_allocate_matrix_set(matrix_w, dft_control%nspins)
215  DO ispin = 1, dft_control%nspins
216  ALLOCATE (matrix_w(ispin)%matrix)
217  CALL dbcsr_copy(matrix_w(ispin)%matrix, matrix_s(1)%matrix, &
218  name="W MATRIX")
219  CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
220  END DO
221  CALL set_ks_env(ks_env, matrix_w=matrix_w)
222 
223  CALL calc_c_mat_force(qs_env)
224  IF (dft_control%do_admm) CALL rt_admm_force(qs_env)
225  IF (dft_control%rtp_control%velocity_gauge .AND. dft_control%rtp_control%nl_gauge_transform) &
226  CALL velocity_gauge_nl_force(qs_env, particle_set)
227  END IF
228  ! from an eventual Mulliken restraint
229  IF (dft_control%qs_control%mulliken_restraint) THEN
230  NULLIFY (matrix_w, matrix_s, rho)
231  CALL get_qs_env(qs_env, &
232  matrix_w=matrix_w, &
233  matrix_s=matrix_s, &
234  rho=rho)
235  NULLIFY (rho_ao)
236  CALL qs_rho_get(rho, rho_ao=rho_ao)
237  CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
238  para_env, matrix_s(1)%matrix, rho_ao, w_matrix=matrix_w)
239  END IF
240  ! Add non-Pulay contribution of DFT+U to W matrix, since it has also to be
241  ! digested with overlap matrix derivatives
242  IF (dft_control%dft_plus_u) THEN
243  NULLIFY (matrix_w)
244  CALL get_qs_env(qs_env, matrix_w=matrix_w)
245  CALL plus_u(qs_env=qs_env, matrix_w=matrix_w)
246  END IF
247 
248  ! Write W Matrix to output (if requested)
249  CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
250  IF (.NOT. has_unit_metric) THEN
251  NULLIFY (matrix_w_kp)
252  CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp)
253  nspin = SIZE(matrix_w_kp, 1)
254  DO ispin = 1, nspin
255  IF (btest(cp_print_key_should_output(logger%iter_info, &
256  qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX"), cp_p_file)) THEN
257  iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX", &
258  extension=".Log")
259  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
260  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
261  after = min(max(after, 1), 16)
262  DO ic = 1, SIZE(matrix_w_kp, 2)
263  CALL cp_dbcsr_write_sparse_matrix(matrix_w_kp(ispin, ic)%matrix, 4, after, qs_env, &
264  para_env, output_unit=iw, omit_headers=omit_headers)
265  END DO
266  CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
267  "DFT%PRINT%AO_MATRICES/W_MATRIX")
268  END IF
269  END DO
270  END IF
271 
272  ! Check if energy correction should be skipped
273  perform_ec = .false.
274  IF (qs_env%energy_correction) THEN
275  CALL get_qs_env(qs_env, ec_env=ec_env)
276  IF (.NOT. ec_env%do_skip) perform_ec = .true.
277  END IF
278 
279  ! Compute core forces (also overwrites matrix_w)
280  IF (dft_control%qs_control%semi_empirical) THEN
281  CALL build_se_core_matrix(qs_env=qs_env, para_env=para_env, &
282  calculate_forces=.true.)
283  CALL se_core_core_interaction(qs_env, para_env, calculate_forces=.true.)
284  ELSEIF (dft_control%qs_control%dftb) THEN
285  CALL build_dftb_matrices(qs_env=qs_env, para_env=para_env, &
286  calculate_forces=.true.)
287  CALL calculate_dftb_dispersion(qs_env=qs_env, para_env=para_env, &
288  calculate_forces=.true.)
289  ELSEIF (dft_control%qs_control%xtb) THEN
290  CALL build_xtb_matrices(qs_env=qs_env, para_env=para_env, &
291  calculate_forces=.true.)
292  ELSEIF (perform_ec) THEN
293  ! Calculates core and grid based forces
294  CALL energy_correction(qs_env, ec_init=.false., calculate_forces=.true.)
295  ELSE
296  ! Dispersion energy and forces are calculated in qs_energy?
297  CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.true.)
298  ! The above line reset the core H, which should be re-updated in case a TD field is applied:
299  IF (qs_env%run_rtp) THEN
300  IF (dft_control%apply_efield_field) &
301  CALL efield_potential_lengh_gauge(qs_env)
302  IF (dft_control%rtp_control%velocity_gauge) &
303  CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.false.)
304 
305  END IF
306  CALL calculate_ecore_self(qs_env)
307  CALL calculate_ecore_overlap(qs_env, para_env, calculate_forces=.true.)
308  CALL calculate_ecore_efield(qs_env, calculate_forces=.true.)
309  !swap external_e_potential before external_c_potential, to ensure
310  !that external potential on grid is loaded before calculating energy of cores
311  CALL external_e_potential(qs_env)
312  IF (.NOT. dft_control%qs_control%gapw) THEN
313  CALL external_c_potential(qs_env, calculate_forces=.true.)
314  END IF
315  ! RIGPW matrices
316  IF (dft_control%qs_control%rigpw) THEN
317  CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
318  CALL build_ri_matrices(lri_env, qs_env, calculate_forces=.true.)
319  END IF
320  END IF
321 
322  ! MP2 Code
323  IF (ASSOCIATED(qs_env%mp2_env)) THEN
324  NULLIFY (energy)
325  CALL get_qs_env(qs_env, energy=energy)
326  CALL qs_scf_compute_properties(qs_env, wf_type='MP2 ', do_mp2=.true.)
327  CALL qs_ks_update_qs_env(qs_env, just_energy=.true.)
328  energy%total = energy%total + energy%mp2
329 
330  IF ((qs_env%mp2_env%method == ri_mp2_method_gpw .OR. qs_env%mp2_env%method == ri_mp2_laplace .OR. &
331  qs_env%mp2_env%method == ri_rpa_method_gpw) &
332  .AND. .NOT. qs_env%mp2_env%do_im_time) THEN
333  CALL update_mp2_forces(qs_env)
334  END IF
335 
336  !RPA EXX energy and forces
337  IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
338  do_exx = .false.
339  hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
340  CALL section_vals_get(hfx_sections, explicit=do_exx)
341  IF (do_exx) THEN
342  do_gw = qs_env%mp2_env%ri_rpa%do_ri_g0w0
343  do_admm = qs_env%mp2_env%ri_rpa%do_admm
344  reuse_hfx = qs_env%mp2_env%ri_rpa%reuse_hfx
345  do_im_time = qs_env%mp2_env%do_im_time
346  output_unit = cp_logger_get_default_io_unit()
347  dummy_real = 0.0_dp
348 
349  CALL calculate_exx(qs_env=qs_env, &
350  unit_nr=output_unit, &
351  hfx_sections=hfx_sections, &
352  x_data=qs_env%mp2_env%ri_rpa%x_data, &
353  do_gw=do_gw, &
354  do_admm=do_admm, &
355  calc_forces=.true., &
356  reuse_hfx=reuse_hfx, &
357  do_im_time=do_im_time, &
358  e_ex_from_gw=dummy_real, &
359  e_admm_from_gw=dummy_real2, &
360  t3=dummy_real)
361  END IF
362  END IF
363  ELSEIF (.NOT. perform_ec) THEN
364  ! Compute grid-based forces
365  CALL qs_ks_update_qs_env(qs_env, calculate_forces=.true.)
366  END IF
367 
368  ! Excited state forces
369  CALL excited_state_energy(qs_env, calculate_forces=.true.)
370 
371  ! replicate forces (get current pointer)
372  NULLIFY (force)
373  CALL get_qs_env(qs_env=qs_env, force=force)
374  CALL replicate_qs_force(force, para_env)
375 
376  DO iatom = 1, natom
377  ikind = kind_of(iatom)
378  i = atom_of_kind(iatom)
379  ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
380  ! the force is - dE/dR, what is called force is actually the gradient
381  ! Things should have the right name
382  ! The minus sign below is a hack
383  ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
384  force(ikind)%other(1:3, i) = -particle_set(iatom)%f(1:3) + force(ikind)%ch_pulay(1:3, i)
385  force(ikind)%total(1:3, i) = force(ikind)%total(1:3, i) + force(ikind)%other(1:3, i)
386  particle_set(iatom)%f = -force(ikind)%total(1:3, i)
387  END DO
388 
389  NULLIFY (virial, energy)
390  CALL get_qs_env(qs_env=qs_env, virial=virial, energy=energy)
391  IF (virial%pv_availability) THEN
392  CALL para_env%sum(virial%pv_overlap)
393  CALL para_env%sum(virial%pv_ekinetic)
394  CALL para_env%sum(virial%pv_ppl)
395  CALL para_env%sum(virial%pv_ppnl)
396  CALL para_env%sum(virial%pv_ecore_overlap)
397  CALL para_env%sum(virial%pv_ehartree)
398  CALL para_env%sum(virial%pv_exc)
399  CALL para_env%sum(virial%pv_exx)
400  CALL para_env%sum(virial%pv_vdw)
401  CALL para_env%sum(virial%pv_mp2)
402  CALL para_env%sum(virial%pv_nlcc)
403  CALL para_env%sum(virial%pv_gapw)
404  CALL para_env%sum(virial%pv_lrigpw)
405  CALL para_env%sum(virial%pv_virial)
406  CALL symmetrize_virial(virial)
407  ! Add the volume terms of the virial
408  IF ((.NOT. virial%pv_numer) .AND. &
409  (.NOT. (dft_control%qs_control%dftb .OR. &
410  dft_control%qs_control%xtb .OR. &
411  dft_control%qs_control%semi_empirical))) THEN
412 
413  ! Harris energy correction requires volume terms from
414  ! 1) Harris functional contribution, and
415  ! 2) Linear Response solver
416  IF (perform_ec) THEN
417  CALL get_qs_env(qs_env, ec_env=ec_env)
418  energy%hartree = ec_env%ehartree
419  energy%exc = ec_env%exc
420  IF (dft_control%do_admm) THEN
421  energy%exc_aux_fit = ec_env%exc_aux_fit
422  END IF
423  END IF
424  DO i = 1, 3
425  virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) &
426  - 2.0_dp*(energy%hartree + energy%sccs_pol)
427  virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc &
428  - 2.0_dp*(energy%hartree + energy%sccs_pol)
429  virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc
430  IF (dft_control%do_admm) THEN
431  virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc_aux_fit
432  virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc_aux_fit
433  END IF
434  ! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve.
435  ! The sign in pw_poisson_solve is correct for FIST, but not for QS.
436  ! There should be a more elegant solution to that ...
437  END DO
438  END IF
439  END IF
440 
441  output_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%DERIVATIVES", &
442  extension=".Log")
443  print_section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%DERIVATIVES")
444  IF (dft_control%qs_control%semi_empirical) THEN
445  CALL write_forces(force, atomic_kind_set, 2, output_unit=output_unit, &
446  print_section=print_section)
447  ELSE IF (dft_control%qs_control%dftb) THEN
448  CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, &
449  print_section=print_section)
450  ELSE IF (dft_control%qs_control%xtb) THEN
451  CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, &
452  print_section=print_section)
453  ELSE IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
454  CALL write_forces(force, atomic_kind_set, 1, output_unit=output_unit, &
455  print_section=print_section)
456  ELSE
457  CALL write_forces(force, atomic_kind_set, 0, output_unit=output_unit, &
458  print_section=print_section)
459  END IF
460  CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, &
461  "DFT%PRINT%DERIVATIVES")
462 
463  ! deallocate W Matrix:
464  NULLIFY (ks_env, matrix_w_kp)
465  CALL get_qs_env(qs_env=qs_env, &
466  matrix_w_kp=matrix_w_kp, &
467  ks_env=ks_env)
468  CALL dbcsr_deallocate_matrix_set(matrix_w_kp)
469  NULLIFY (matrix_w_kp)
470  CALL set_ks_env(ks_env, matrix_w_kp=matrix_w_kp)
471 
472  DEALLOCATE (atom_of_kind, kind_of)
473 
474  CALL timestop(handle)
475 
476  END SUBROUTINE qs_forces
477 
478 ! **************************************************************************************************
479 !> \brief Write a Quickstep force data structure to output unit
480 !> \param qs_force ...
481 !> \param atomic_kind_set ...
482 !> \param ftype ...
483 !> \param output_unit ...
484 !> \param print_section ...
485 !> \date 05.06.2002
486 !> \author MK
487 !> \version 1.0
488 ! **************************************************************************************************
489  SUBROUTINE write_forces(qs_force, atomic_kind_set, ftype, output_unit, &
490  print_section)
491 
492  TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
493  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
494  INTEGER, INTENT(IN) :: ftype, output_unit
495  TYPE(section_vals_type), POINTER :: print_section
496 
497  CHARACTER(LEN=13) :: fmtstr5
498  CHARACTER(LEN=15) :: fmtstr4
499  CHARACTER(LEN=20) :: fmtstr3
500  CHARACTER(LEN=35) :: fmtstr2
501  CHARACTER(LEN=48) :: fmtstr1
502  INTEGER :: i, iatom, ikind, my_ftype, natom, ndigits
503  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
504  REAL(kind=dp), DIMENSION(3) :: grand_total
505 
506  IF (output_unit > 0) THEN
507 
508  IF (.NOT. ASSOCIATED(qs_force)) THEN
509  CALL cp_abort(__location__, &
510  "The qs_force pointer is not associated "// &
511  "and cannot be printed")
512  END IF
513 
514  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, &
515  kind_of=kind_of, natom=natom)
516 
517  ! Variable precision output of the forces
518  CALL section_vals_val_get(print_section, "NDIGITS", &
519  i_val=ndigits)
520 
521  fmtstr1 = "(/,/,T2,A,/,/,T3,A,T11,A,T23,A,T40,A1,2( X,A1))"
522  WRITE (unit=fmtstr1(41:42), fmt="(I2)") ndigits + 5
523 
524  fmtstr2 = "(/,(T2,I5,4X,I4,T18,A,T34,3F . ))"
525  WRITE (unit=fmtstr2(32:33), fmt="(I2)") ndigits
526  WRITE (unit=fmtstr2(29:30), fmt="(I2)") ndigits + 6
527 
528  fmtstr3 = "(/,T3,A,T34,3F . )"
529  WRITE (unit=fmtstr3(18:19), fmt="(I2)") ndigits
530  WRITE (unit=fmtstr3(15:16), fmt="(I2)") ndigits + 6
531 
532  fmtstr4 = "((T34,3F . ))"
533  WRITE (unit=fmtstr4(12:13), fmt="(I2)") ndigits
534  WRITE (unit=fmtstr4(9:10), fmt="(I2)") ndigits + 6
535 
536  fmtstr5 = "(/T2,A//T3,A)"
537 
538  WRITE (unit=output_unit, fmt=fmtstr1) &
539  "FORCES [a.u.]", "Atom", "Kind", "Component", "X", "Y", "Z"
540 
541  grand_total(:) = 0.0_dp
542 
543  my_ftype = ftype
544 
545  SELECT CASE (my_ftype)
546  CASE DEFAULT
547  DO iatom = 1, natom
548  ikind = kind_of(iatom)
549  i = atom_of_kind(iatom)
550  WRITE (unit=output_unit, fmt=fmtstr2) &
551  iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
552  grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
553  END DO
554  CASE (0)
555  DO iatom = 1, natom
556  ikind = kind_of(iatom)
557  i = atom_of_kind(iatom)
558  WRITE (unit=output_unit, fmt=fmtstr2) &
559  iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
560  iatom, ikind, " overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
561  iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
562  iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
563  iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
564  iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
565  iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
566  iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
567  iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
568  iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
569  iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
570  iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
571  iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), &
572  iatom, ikind, " other", qs_force(ikind)%other(1:3, i), &
573  iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
574  iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
575  iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
576  iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), &
577  iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
578  iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
579  grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
580  END DO
581  CASE (1)
582  DO iatom = 1, natom
583  ikind = kind_of(iatom)
584  i = atom_of_kind(iatom)
585  WRITE (unit=output_unit, fmt=fmtstr2) &
586  iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
587  iatom, ikind, " overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
588  iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
589  iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
590  iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
591  iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
592  iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
593  iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
594  iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
595  iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
596  iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
597  iatom, ikind, " vhxc_atom", qs_force(ikind)%vhxc_atom(1:3, i), &
598  iatom, ikind, " g0s_Vh_elec", qs_force(ikind)%g0s_Vh_elec(1:3, i), &
599  iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
600  iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
601  iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), &
602  iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
603  iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
604  iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
605  iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), &
606  iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
607  iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
608  grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
609  END DO
610  CASE (2)
611  DO iatom = 1, natom
612  ikind = kind_of(iatom)
613  i = atom_of_kind(iatom)
614  WRITE (unit=output_unit, fmt=fmtstr2) &
615  iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
616  iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
617  iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
618  grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
619  END DO
620  CASE (3)
621  DO iatom = 1, natom
622  ikind = kind_of(iatom)
623  i = atom_of_kind(iatom)
624  WRITE (unit=output_unit, fmt=fmtstr2) &
625  iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
626  iatom, ikind, "overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
627  iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
628  iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
629  iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
630  iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
631  iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
632  iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
633  iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
634  iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
635  iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
636  iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
637  iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
638  grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
639  END DO
640  CASE (4)
641  DO iatom = 1, natom
642  ikind = kind_of(iatom)
643  i = atom_of_kind(iatom)
644  WRITE (unit=output_unit, fmt=fmtstr2) &
645  iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
646  iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
647  iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
648  iatom, ikind, " repulsive", qs_force(ikind)%repulsive(1:3, i), &
649  iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
650  iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
651  iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
652  iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
653  grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
654  END DO
655  CASE (5)
656  DO iatom = 1, natom
657  ikind = kind_of(iatom)
658  i = atom_of_kind(iatom)
659  WRITE (unit=output_unit, fmt=fmtstr2) &
660  iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
661  iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
662  iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
663  iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
664  iatom, ikind, " all potential", qs_force(ikind)%all_potential(1:3, i), &
665  iatom, ikind, " other", qs_force(ikind)%other(1:3, i), &
666  iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
667  grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
668  END DO
669  END SELECT
670 
671  WRITE (unit=output_unit, fmt=fmtstr3) "Sum of total", grand_total(1:3)
672 
673  DEALLOCATE (atom_of_kind)
674  DEALLOCATE (kind_of)
675 
676  END IF
677 
678  END SUBROUTINE write_forces
679 
680 END MODULE qs_force
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
DBCSR output in CP2K.
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Add the DFT+U contribution to the Hamiltonian matrix.
Definition: dft_plus_u.F:18
subroutine, public plus_u(qs_env, matrix_h, matrix_w)
Add the DFT+U contribution to the Hamiltonian matrix. Wrapper routine for all "+U" methods.
Definition: dft_plus_u.F:98
Types needed for a for a Energy Correction.
Definition: ec_env_types.F:14
all routins needed for a nonperiodic electric field
Definition: efield_utils.F:12
subroutine, public calculate_ecore_efield(qs_env, calculate_forces)
Computes the force and the energy due to a efield on the cores Note: In the velocity gauge,...
Definition: efield_utils.F:186
subroutine, public efield_potential_lengh_gauge(qs_env)
Replace the original implementation of the electric-electronic interaction in the length gauge....
Definition: efield_utils.F:65
Routines for an energy correction on top of a Kohn-Sham calculation.
subroutine, public energy_correction(qs_env, ec_init, calculate_forces)
Energy Correction to a Kohn-Sham simulation Available energy corrections: (1) Harris energy functiona...
Routines for total energy and forces of excited states.
subroutine, public excited_state_energy(qs_env, calculate_forces)
Excited state energy and forces.
Routines to calculate EXX in RPA and energy correction methods.
Definition: hfx_exx.F:16
subroutine, public calculate_exx(qs_env, unit_nr, hfx_sections, x_data, do_gw, do_admm, calc_forces, reuse_hfx, do_im_time, E_ex_from_GW, E_admm_from_GW, t3)
...
Definition: hfx_exx.F:106
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public ri_rpa_method_gpw
integer, parameter, public ri_mp2_method_gpw
integer, parameter, public ri_mp2_laplace
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Interface to the message passing library MPI.
Routines to calculate CPHF like update and solve Z-vector equation for MP2 gradients (only GPW)
Definition: mp2_cphf.F:14
subroutine, public update_mp2_forces(qs_env)
...
Definition: mp2_cphf.F:1301
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition: mulliken.F:13
subroutine, public mulliken_restraint(mulliken_restraint_control, para_env, s_matrix, p_matrix, energy, order_p, ks_matrix, w_matrix)
computes the energy and density matrix derivate of a constraint on the mulliken charges
Definition: mulliken.F:73
Define the data structure for the particle information.
Calculation of the energies concerning the core charge distribution.
subroutine, public calculate_ecore_self(qs_env, E_self_core, atecc)
Calculate the self energy of the core charge distribution.
subroutine, public calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular, E_overlap_core, atecc)
Calculate the overlap energy of the core charge distribution.
Calculation of the core Hamiltonian integral matrix <a|H|b> over Cartesian Gaussian-type functions.
subroutine, public build_core_hamiltonian_matrix(qs_env, calculate_forces)
Cosntruction of the QS Core Hamiltonian Matrix.
Calculation of dispersion in DFTB.
subroutine, public calculate_dftb_dispersion(qs_env, para_env, calculate_forces)
...
Calculation of Overlap and Hamiltonian matrices in DFTB.
subroutine, public build_dftb_matrices(qs_env, para_env, calculate_forces)
...
Perform a QUICKSTEP wavefunction optimization (single point)
Definition: qs_energy.F:14
subroutine, public qs_energies(qs_env, consistent_energies, calc_forces)
Driver routine for QUICKSTEP single point wavefunction optimization.
Definition: qs_energy.F:65
qs_environment methods that use many other modules
subroutine, public qs_env_rebuild_pw_env(qs_env)
rebuilds the pw_env in the given qs_env, allocating it if necessary
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Routines to handle an external electrostatic field The external field can be generic and is provided ...
subroutine, public external_c_potential(qs_env, calculate_forces)
Computes the force and the energy due to the external potential on the cores.
subroutine, public external_e_potential(qs_env)
Computes the external potential on the grid.
subroutine, public replicate_qs_force(qs_force, para_env)
Replicate and sum up the force.
subroutine, public zero_qs_force(qs_force)
Initialize a Quickstep force data structure.
subroutine, public allocate_qs_force(qs_force, natom_of_kind)
Allocate a Quickstep force data structure.
Quickstep force driver routine.
Definition: qs_force.F:12
subroutine, public qs_calc_energy_force(qs_env, calc_force, consistent_energies, linres)
...
Definition: qs_force.F:107
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition: qs_ks_types.F:567
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
Utility routines for qs_scf.
subroutine, public qs_scf_compute_properties(qs_env, wf_type, do_mp2)
computes properties for a given hamilonian using the current wfn
types that represent a quickstep subsys
subroutine, public qs_subsys_set(subsys, cp_subsys, local_particles, local_molecules, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, nelectron_total, nelectron_spin)
...
Calculates integral matrices for RIGPW method.
subroutine, public build_ri_matrices(lri_env, qs_env, calculate_forces)
creates and initializes an lri_env
Routines needed for EMD.
subroutine, public rt_admm_force(qs_env)
...
subroutine, public calc_c_mat_force(qs_env)
calculates the three additional force contributions needed in EMD P_imag*C , P_imag*B*S^-1*S_der ,...
Routines to perform the RTP in the velocity gauge.
subroutine, public velocity_gauge_nl_force(qs_env, particle_set)
Calculate the force associated to non-local pseudo potential in the velocity gauge.
subroutine, public velocity_gauge_ks_matrix(qs_env, subtract_nl_term)
...
Split and build its own idependent core_core SE interaction module.
Definition: se_core_core.F:14
subroutine, public se_core_core_interaction(qs_env, para_env, calculate_forces)
Evaluates the core-core interactions for NDDO methods.
Definition: se_core_core.F:79
Calculation of the Hamiltonian integral matrix <a|H|b> for semi-empirical methods.
subroutine, public build_se_core_matrix(qs_env, para_env, calculate_forces)
...
subroutine, public symmetrize_virial(virial)
Symmetrize the virial components.
Definition: virial_types.F:67
Calculation of Overlap and Hamiltonian matrices in xTB Reference: Stefan Grimme, Christoph Bannwarth,...
Definition: xtb_matrices.F:15
subroutine, public build_xtb_matrices(qs_env, para_env, calculate_forces)
...
Definition: xtb_matrices.F:134