(git:d657d38)
Loading...
Searching...
No Matches
qs_force.F
Go to the documentation of this file.
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! **************************************************************************************************
16 USE cp_dbcsr_api, ONLY: dbcsr_copy,&
25 USE cp_output_handling, ONLY: cp_p_file,&
29 USE dft_plus_u, ONLY: plus_u
35 USE hfx_exx, ONLY: calculate_exx
43 USE kinds, ONLY: dp
54 USE qs_energy, ONLY: qs_energies
66 USE qs_ks_types, ONLY: qs_ks_env_type,&
68 USE qs_rho_types, ONLY: qs_rho_get,&
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
98CONTAINS
99
100! **************************************************************************************************
101!> \brief ...
102!> \param qs_env ...
103!> \param calc_force ...
104!> \param consistent_energies ...
105!> \param linres ...
106! **************************************************************************************************
107 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 qs_env%linres_run = linres
112 IF (calc_force) THEN
113 CALL qs_forces(qs_env)
114 ELSE
115 CALL qs_energies(qs_env, calc_forces=.false., &
116 consistent_energies=consistent_energies)
117 END IF
118
119 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 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 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 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
143 TYPE(cp_logger_type), POINTER :: logger
144 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_w, rho_ao
145 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
151 TYPE(qs_energy_type), POINTER :: energy
152 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 CALL timeset(routinen, handle)
160 NULLIFY (logger)
161 logger => cp_get_default_logger()
162
163 ! rebuild plane wave environment
164 CALL qs_env_rebuild_pw_env(qs_env)
165
166 ! zero out the forces in particle set
167 CALL get_qs_env(qs_env, particle_set=particle_set)
168 natom = SIZE(particle_set)
169 DO iatom = 1, natom
170 particle_set(iatom)%f = 0.0_dp
171 END DO
172
173 ! get atom mapping
174 NULLIFY (atomic_kind_set)
175 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 kind_of=kind_of)
179
180 NULLIFY (force, subsys, dft_control)
181 CALL get_qs_env(qs_env, &
182 force=force, &
183 subsys=subsys, &
184 dft_control=dft_control)
185 IF (.NOT. ASSOCIATED(force)) THEN
186 ! *** Allocate the force data structure ***
187 nkind = SIZE(atomic_kind_set)
188 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
189 CALL allocate_qs_force(force, natom_of_kind)
190 DEALLOCATE (natom_of_kind)
191 CALL qs_subsys_set(subsys, force=force)
192 END IF
193 CALL zero_qs_force(force)
194
195 ! Check if CDFT potential is needed and save it until forces have been calculated
196 IF (dft_control%qs_control%cdft) THEN
197 dft_control%qs_control%cdft_control%save_pot = .true.
198 END IF
199
200 ! recalculate energy with forces
201 CALL qs_energies(qs_env, calc_forces=.true.)
202
203 NULLIFY (para_env)
204 CALL get_qs_env(qs_env, &
205 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 IF (qs_env%run_rtp) THEN
210 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 matrix_s=matrix_s)
215 CALL dbcsr_allocate_matrix_set(matrix_w, dft_control%nspins)
216 DO ispin = 1, dft_control%nspins
217 ALLOCATE (matrix_w(ispin)%matrix)
218 CALL dbcsr_copy(matrix_w(ispin)%matrix, matrix_s(1)%matrix, &
219 name="W MATRIX")
220 CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
221 END DO
222 CALL set_ks_env(ks_env, matrix_w=matrix_w)
223
224 CALL calc_c_mat_force(qs_env)
225 IF (dft_control%do_admm) CALL rt_admm_force(qs_env)
226 IF (dft_control%rtp_control%velocity_gauge .AND. dft_control%rtp_control%nl_gauge_transform) &
227 CALL velocity_gauge_nl_force(qs_env, particle_set)
228 END IF
229 ! from an eventual Mulliken restraint
230 IF (dft_control%qs_control%mulliken_restraint) THEN
231 NULLIFY (matrix_w, matrix_s, rho)
232 CALL get_qs_env(qs_env, &
233 matrix_w=matrix_w, &
234 matrix_s=matrix_s, &
235 rho=rho)
236 NULLIFY (rho_ao)
237 CALL qs_rho_get(rho, rho_ao=rho_ao)
238 CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
239 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 IF (dft_control%dft_plus_u) THEN
244 NULLIFY (matrix_w)
245 CALL get_qs_env(qs_env, matrix_w=matrix_w)
246 CALL plus_u(qs_env=qs_env, matrix_w=matrix_w)
247 END IF
248
249 ! Write W Matrix to output (if requested)
250 CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
251 IF (.NOT. has_unit_metric) THEN
252 NULLIFY (matrix_w_kp)
253 CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp)
254 nspin = SIZE(matrix_w_kp, 1)
255 DO ispin = 1, nspin
256 IF (btest(cp_print_key_should_output(logger%iter_info, &
257 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 extension=".Log")
260 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
261 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
262 after = min(max(after, 1), 16)
263 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 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 "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 perform_ec = .false.
275 IF (qs_env%energy_correction) THEN
276 CALL get_qs_env(qs_env, ec_env=ec_env)
277 IF (.NOT. ec_env%do_skip) perform_ec = .true.
278 END IF
279
280 ! Compute core forces (also overwrites matrix_w)
281 IF (dft_control%qs_control%semi_empirical) THEN
282 CALL build_se_core_matrix(qs_env=qs_env, para_env=para_env, &
283 calculate_forces=.true.)
284 CALL se_core_core_interaction(qs_env, para_env, calculate_forces=.true.)
285 ELSEIF (dft_control%qs_control%dftb) THEN
286 CALL build_dftb_matrices(qs_env=qs_env, para_env=para_env, &
287 calculate_forces=.true.)
288 CALL calculate_dftb_dispersion(qs_env=qs_env, para_env=para_env, &
289 calculate_forces=.true.)
290 ELSEIF (dft_control%qs_control%xtb) THEN
291 IF (dft_control%qs_control%xtb_control%do_tblite) THEN
292 CALL build_tblite_matrices(qs_env=qs_env, calculate_forces=.true.)
293 ELSE
294 CALL build_xtb_matrices(qs_env=qs_env, calculate_forces=.true.)
295 END IF
296 ELSEIF (perform_ec) THEN
297 ! Calculates core and grid based forces
298 CALL energy_correction(qs_env, ec_init=.false., calculate_forces=.true.)
299 ELSE
300 ! Dispersion energy and forces are calculated in qs_energy?
301 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 IF (qs_env%run_rtp) THEN
304 IF (dft_control%apply_efield_field) &
306 IF (dft_control%rtp_control%velocity_gauge) &
307 CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.false.)
308
309 END IF
310 CALL calculate_ecore_self(qs_env)
311 CALL calculate_ecore_overlap(qs_env, para_env, calculate_forces=.true.)
312 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 CALL external_e_potential(qs_env)
316 IF (.NOT. dft_control%qs_control%gapw) THEN
317 CALL external_c_potential(qs_env, calculate_forces=.true.)
318 END IF
319 ! RIGPW matrices
320 IF (dft_control%qs_control%rigpw) THEN
321 CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
322 CALL build_ri_matrices(lri_env, qs_env, calculate_forces=.true.)
323 END IF
324 END IF
325
326 ! MP2 Code
327 IF (ASSOCIATED(qs_env%mp2_env)) THEN
328 NULLIFY (energy)
329 CALL get_qs_env(qs_env, energy=energy)
330 CALL qs_scf_compute_properties(qs_env, wf_type='MP2 ', do_mp2=.true.)
331 CALL qs_ks_update_qs_env(qs_env, just_energy=.true.)
332 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 .AND. .NOT. qs_env%mp2_env%do_im_time) THEN
337 CALL update_mp2_forces(qs_env)
338 END IF
339
340 !RPA EXX energy and forces
341 IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
342 do_exx = .false.
343 hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
344 CALL section_vals_get(hfx_sections, explicit=do_exx)
345 IF (do_exx) THEN
346 do_gw = qs_env%mp2_env%ri_rpa%do_ri_g0w0
347 do_admm = qs_env%mp2_env%ri_rpa%do_admm
348 reuse_hfx = qs_env%mp2_env%ri_rpa%reuse_hfx
349 do_im_time = qs_env%mp2_env%do_im_time
350 output_unit = cp_logger_get_default_io_unit()
351 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 t3=dummy_real)
365 END IF
366 END IF
367 ELSEIF (perform_ec) THEN
368 ! energy correction forces postponed
369 ELSEIF (qs_env%harris_method) THEN
370 ! Harris method forces already done in harris_energy_correction
371 ELSE
372 ! Compute grid-based forces
373 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.true.)
374 END IF
375
376 ! Excited state forces
377 CALL excited_state_energy(qs_env, calculate_forces=.true.)
378
379 ! replicate forces (get current pointer)
380 NULLIFY (force)
381 CALL get_qs_env(qs_env=qs_env, force=force)
382 CALL replicate_qs_force(force, para_env)
383
384 DO iatom = 1, natom
385 ikind = kind_of(iatom)
386 i = atom_of_kind(iatom)
387 ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
388 ! the force is - dE/dR, what is called force is actually the gradient
389 ! Things should have the right name
390 ! The minus sign below is a hack
391 ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
392 force(ikind)%other(1:3, i) = -particle_set(iatom)%f(1:3) + force(ikind)%ch_pulay(1:3, i)
393 force(ikind)%total(1:3, i) = force(ikind)%total(1:3, i) + force(ikind)%other(1:3, i)
394 particle_set(iatom)%f = -force(ikind)%total(1:3, i)
395 END DO
396
397 NULLIFY (virial, energy)
398 CALL get_qs_env(qs_env=qs_env, virial=virial, energy=energy)
399 IF (virial%pv_availability) THEN
400 CALL para_env%sum(virial%pv_overlap)
401 CALL para_env%sum(virial%pv_ekinetic)
402 CALL para_env%sum(virial%pv_ppl)
403 CALL para_env%sum(virial%pv_ppnl)
404 CALL para_env%sum(virial%pv_ecore_overlap)
405 CALL para_env%sum(virial%pv_ehartree)
406 CALL para_env%sum(virial%pv_exc)
407 CALL para_env%sum(virial%pv_exx)
408 CALL para_env%sum(virial%pv_vdw)
409 CALL para_env%sum(virial%pv_mp2)
410 CALL para_env%sum(virial%pv_nlcc)
411 CALL para_env%sum(virial%pv_gapw)
412 CALL para_env%sum(virial%pv_lrigpw)
413 CALL para_env%sum(virial%pv_virial)
414 CALL symmetrize_virial(virial)
415 ! Add the volume terms of the virial
416 IF ((.NOT. virial%pv_numer) .AND. &
417 (.NOT. (dft_control%qs_control%dftb .OR. &
418 dft_control%qs_control%xtb .OR. &
419 dft_control%qs_control%semi_empirical))) THEN
420
421 ! Harris energy correction requires volume terms from
422 ! 1) Harris functional contribution, and
423 ! 2) Linear Response solver
424 IF (perform_ec) THEN
425 CALL get_qs_env(qs_env, ec_env=ec_env)
426 energy%hartree = ec_env%ehartree
427 energy%exc = ec_env%exc
428 IF (dft_control%do_admm) THEN
429 energy%exc_aux_fit = ec_env%exc_aux_fit
430 END IF
431 END IF
432 DO i = 1, 3
433 virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) &
434 - 2.0_dp*(energy%hartree + energy%sccs_pol)
435 virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc &
436 - 2.0_dp*(energy%hartree + energy%sccs_pol)
437 virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc
438 IF (dft_control%do_admm) THEN
439 virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc_aux_fit
440 virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc_aux_fit
441 END IF
442 ! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve.
443 ! The sign in pw_poisson_solve is correct for FIST, but not for QS.
444 ! There should be a more elegant solution to that ...
445 END DO
446 END IF
447 END IF
448
449 output_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%DERIVATIVES", &
450 extension=".Log")
451 print_section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%DERIVATIVES")
452 IF (dft_control%qs_control%semi_empirical) THEN
453 CALL write_forces(force, atomic_kind_set, 2, output_unit=output_unit, &
454 print_section=print_section)
455 ELSE IF (dft_control%qs_control%dftb) THEN
456 CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, &
457 print_section=print_section)
458 ELSE IF (dft_control%qs_control%xtb) THEN
459 CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, &
460 print_section=print_section)
461 ELSE IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
462 CALL write_forces(force, atomic_kind_set, 1, output_unit=output_unit, &
463 print_section=print_section)
464 ELSE
465 CALL write_forces(force, atomic_kind_set, 0, output_unit=output_unit, &
466 print_section=print_section)
467 END IF
468 CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, &
469 "DFT%PRINT%DERIVATIVES")
470
471 ! deallocate W Matrix:
472 NULLIFY (ks_env, matrix_w_kp)
473 CALL get_qs_env(qs_env=qs_env, &
474 matrix_w_kp=matrix_w_kp, &
475 ks_env=ks_env)
476 CALL dbcsr_deallocate_matrix_set(matrix_w_kp)
477 NULLIFY (matrix_w_kp)
478 CALL set_ks_env(ks_env, matrix_w_kp=matrix_w_kp)
479
480 DEALLOCATE (atom_of_kind, kind_of)
481
482 CALL timestop(handle)
483
484 END SUBROUTINE qs_forces
485
486! **************************************************************************************************
487!> \brief Write a Quickstep force data structure to output unit
488!> \param qs_force ...
489!> \param atomic_kind_set ...
490!> \param ftype ...
491!> \param output_unit ...
492!> \param print_section ...
493!> \date 05.06.2002
494!> \author MK
495!> \version 1.0
496! **************************************************************************************************
497 SUBROUTINE write_forces(qs_force, atomic_kind_set, ftype, output_unit, &
498 print_section)
499
500 TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
501 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
502 INTEGER, INTENT(IN) :: ftype, output_unit
503 TYPE(section_vals_type), POINTER :: print_section
504
505 CHARACTER(LEN=13) :: fmtstr5
506 CHARACTER(LEN=15) :: fmtstr4
507 CHARACTER(LEN=20) :: fmtstr3
508 CHARACTER(LEN=35) :: fmtstr2
509 CHARACTER(LEN=48) :: fmtstr1
510 INTEGER :: i, iatom, ikind, my_ftype, natom, ndigits
511 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
512 REAL(kind=dp), DIMENSION(3) :: grand_total
513
514 IF (output_unit > 0) THEN
515
516 IF (.NOT. ASSOCIATED(qs_force)) THEN
517 CALL cp_abort(__location__, &
518 "The qs_force pointer is not associated "// &
519 "and cannot be printed")
520 END IF
521
522 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, &
523 kind_of=kind_of, natom=natom)
524
525 ! Variable precision output of the forces
526 CALL section_vals_val_get(print_section, "NDIGITS", &
527 i_val=ndigits)
528
529 fmtstr1 = "(/,/,T2,A,/,/,T3,A,T11,A,T23,A,T40,A1,2( X,A1))"
530 WRITE (unit=fmtstr1(41:42), fmt="(I2)") ndigits + 5
531
532 fmtstr2 = "(/,(T2,I5,4X,I4,T18,A,T34,3F . ))"
533 WRITE (unit=fmtstr2(32:33), fmt="(I2)") ndigits
534 WRITE (unit=fmtstr2(29:30), fmt="(I2)") ndigits + 6
535
536 fmtstr3 = "(/,T3,A,T34,3F . )"
537 WRITE (unit=fmtstr3(18:19), fmt="(I2)") ndigits
538 WRITE (unit=fmtstr3(15:16), fmt="(I2)") ndigits + 6
539
540 fmtstr4 = "((T34,3F . ))"
541 WRITE (unit=fmtstr4(12:13), fmt="(I2)") ndigits
542 WRITE (unit=fmtstr4(9:10), fmt="(I2)") ndigits + 6
543
544 fmtstr5 = "(/T2,A//T3,A)"
545
546 WRITE (unit=output_unit, fmt=fmtstr1) &
547 "FORCES [a.u.]", "Atom", "Kind", "Component", "X", "Y", "Z"
548
549 grand_total(:) = 0.0_dp
550
551 my_ftype = ftype
552
553 SELECT CASE (my_ftype)
554 CASE DEFAULT
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, " total", qs_force(ikind)%total(1:3, i)
560 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
561 END DO
562 CASE (0)
563 DO iatom = 1, natom
564 ikind = kind_of(iatom)
565 i = atom_of_kind(iatom)
566 WRITE (unit=output_unit, fmt=fmtstr2) &
567 iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
568 iatom, ikind, " overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
569 iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
570 iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
571 iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
572 iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
573 iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
574 iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
575 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
576 iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
577 iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
578 iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
579 iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), &
580 iatom, ikind, " other", qs_force(ikind)%other(1:3, i), &
581 iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
582 iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
583 iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
584 iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), &
585 iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
586 iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
587 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
588 END DO
589 CASE (1)
590 DO iatom = 1, natom
591 ikind = kind_of(iatom)
592 i = atom_of_kind(iatom)
593 WRITE (unit=output_unit, fmt=fmtstr2) &
594 iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
595 iatom, ikind, " overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
596 iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
597 iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
598 iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
599 iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
600 iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
601 iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
602 iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
603 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
604 iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
605 iatom, ikind, " vhxc_atom", qs_force(ikind)%vhxc_atom(1:3, i), &
606 iatom, ikind, " g0s_Vh_elec", qs_force(ikind)%g0s_Vh_elec(1:3, i), &
607 iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
608 iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
609 iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), &
610 iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
611 iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
612 iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
613 iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), &
614 iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
615 iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
616 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
617 END DO
618 CASE (2)
619 DO iatom = 1, natom
620 ikind = kind_of(iatom)
621 i = atom_of_kind(iatom)
622 WRITE (unit=output_unit, fmt=fmtstr2) &
623 iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
624 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
625 iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
626 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
627 END DO
628 CASE (3)
629 DO iatom = 1, natom
630 ikind = kind_of(iatom)
631 i = atom_of_kind(iatom)
632 WRITE (unit=output_unit, fmt=fmtstr2) &
633 iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
634 iatom, ikind, "overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
635 iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
636 iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
637 iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
638 iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
639 iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
640 iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
641 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
642 iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
643 iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
644 iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
645 iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
646 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
647 END DO
648 CASE (4)
649 DO iatom = 1, natom
650 ikind = kind_of(iatom)
651 i = atom_of_kind(iatom)
652 WRITE (unit=output_unit, fmt=fmtstr2) &
653 iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
654 iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
655 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
656 iatom, ikind, " repulsive", qs_force(ikind)%repulsive(1:3, i), &
657 iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
658 iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
659 iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
660 iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
661 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
662 END DO
663 CASE (5)
664 DO iatom = 1, natom
665 ikind = kind_of(iatom)
666 i = atom_of_kind(iatom)
667 WRITE (unit=output_unit, fmt=fmtstr2) &
668 iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
669 iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
670 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
671 iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
672 iatom, ikind, " all potential", qs_force(ikind)%all_potential(1:3, i), &
673 iatom, ikind, " other", qs_force(ikind)%other(1:3, i), &
674 iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
675 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
676 END DO
677 END SELECT
678
679 WRITE (unit=output_unit, fmt=fmtstr3) "Sum of total", grand_total(1:3)
680
681 DEALLOCATE (atom_of_kind)
682 DEALLOCATE (kind_of)
683
684 END IF
685
686 END SUBROUTINE write_forces
687
688END 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...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_set(matrix, alpha)
...
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:99
Types needed for a for a Energy Correction.
all routins needed for a nonperiodic electric field
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,...
subroutine, public efield_potential_lengh_gauge(qs_env)
Replace the original implementation of the electric-electronic interaction in the length gauge....
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:1305
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_overlap(qs_env, para_env, calculate_forces, molecular, e_overlap_core, atecc)
Calculate the overlap energy of the core charge distribution.
subroutine, public calculate_ecore_self(qs_env, e_self_core, atecc)
Calculate the self 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:70
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_pp, 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, harris_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, eeq, rhs, tb_tblite)
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:108
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
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_pp, 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)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
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...
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.
subroutine, public se_core_core_interaction(qs_env, para_env, calculate_forces)
Evaluates the core-core interactions for NDDO methods.
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)
...
interface to tblite
subroutine, public build_tblite_matrices(qs_env, calculate_forces)
...
subroutine, public symmetrize_virial(virial)
Symmetrize the virial components.
Calculation of Overlap and Hamiltonian matrices in xTB Reference: Stefan Grimme, Christoph Bannwarth,...
subroutine, public build_xtb_matrices(qs_env, calculate_forces)
...
Provides all information about an atomic kind.
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Contains information on the energy correction functional for KG.
stores all the informations relevant to an mpi environment
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.