1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Quickstep force driver routine 8!> \author MK (12.06.2002) 9! ************************************************************************************************** 10MODULE qs_force 11 USE admm_methods, ONLY: calc_aux_mo_derivs_none,& 12 calc_mixed_overlap_force 13 USE admm_types, ONLY: admm_type 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind_set 16 USE cp_control_types, ONLY: dft_control_type 17 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,& 18 dbcsr_deallocate_matrix_set 19 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 20 USE cp_fm_types, ONLY: cp_fm_type 21 USE cp_log_handling, ONLY: cp_get_default_logger,& 22 cp_logger_type 23 USE cp_output_handling, ONLY: cp_p_file,& 24 cp_print_key_finished_output,& 25 cp_print_key_should_output,& 26 cp_print_key_unit_nr 27 USE cp_para_types, ONLY: cp_para_env_type 28 USE dbcsr_api, ONLY: dbcsr_add,& 29 dbcsr_copy,& 30 dbcsr_p_type,& 31 dbcsr_set 32 USE dft_plus_u, ONLY: plus_u 33 USE efield_utils, ONLY: calculate_ecore_efield 34 USE input_constants, ONLY: do_admm_purify_none 35 USE input_section_types, ONLY: section_vals_get_subs_vals,& 36 section_vals_type,& 37 section_vals_val_get 38 USE kinds, ONLY: dp 39 USE lri_environment_types, ONLY: lri_environment_type 40 USE message_passing, ONLY: mp_sum 41 USE mulliken, ONLY: mulliken_restraint 42 USE particle_types, ONLY: particle_type 43 USE qs_core_energies, ONLY: calculate_ecore_overlap,& 44 calculate_ecore_self 45 USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix 46 USE qs_dftb_dispersion, ONLY: calculate_dftb_dispersion 47 USE qs_dftb_matrices, ONLY: build_dftb_matrices 48 USE qs_energy, ONLY: qs_energies 49 USE qs_energy_types, ONLY: qs_energy_type 50 USE qs_environment_methods, ONLY: qs_env_rebuild_pw_env 51 USE qs_environment_types, ONLY: get_qs_env,& 52 qs_environment_type,& 53 set_qs_env 54 USE qs_external_potential, ONLY: external_c_potential,& 55 external_e_potential 56 USE qs_force_types, ONLY: allocate_qs_force,& 57 qs_force_type,& 58 replicate_qs_force,& 59 zero_qs_force 60 USE qs_ks_methods, ONLY: qs_ks_update_qs_env 61 USE qs_ks_types, ONLY: qs_ks_did_change,& 62 qs_ks_env_type,& 63 set_ks_env 64 USE qs_mo_types, ONLY: get_mo_set,& 65 mo_set_p_type,& 66 mo_set_type 67 USE qs_rho_methods, ONLY: qs_rho_update_rho 68 USE qs_rho_types, ONLY: qs_rho_get,& 69 qs_rho_type 70 USE qs_scf_post_scf, ONLY: qs_scf_compute_properties 71 USE qs_subsys_types, ONLY: qs_subsys_set,& 72 qs_subsys_type 73 USE ri_environment_methods, ONLY: build_ri_matrices 74 USE rt_propagation_forces, ONLY: calc_c_mat_force,& 75 rt_admm_force 76 USE se_core_core, ONLY: se_core_core_interaction 77 USE se_core_matrix, ONLY: build_se_core_matrix 78 USE virial_types, ONLY: virial_type 79 USE xtb_matrices, ONLY: build_xtb_matrices 80#include "./base/base_uses.f90" 81 82 IMPLICIT NONE 83 84 PRIVATE 85 86! *** Global parameters *** 87 88 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_force' 89 90! *** Public subroutines *** 91 92 PUBLIC :: qs_calc_energy_force 93 94CONTAINS 95 96! ************************************************************************************************** 97!> \brief ... 98!> \param qs_env ... 99!> \param calc_force ... 100!> \param consistent_energies ... 101!> \param linres ... 102! ************************************************************************************************** 103 SUBROUTINE qs_calc_energy_force(qs_env, calc_force, consistent_energies, linres) 104 TYPE(qs_environment_type), POINTER :: qs_env 105 LOGICAL :: calc_force, consistent_energies, linres 106 107 CHARACTER(len=*), PARAMETER :: routineN = 'qs_calc_energy_force', & 108 routineP = moduleN//':'//routineN 109 110 qs_env%linres_run = linres 111 CALL set_qs_env(qs_env) 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 CALL get_qs_env(qs_env) 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', routineP = moduleN//':'//routineN 133 134 INTEGER :: after, dir, 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 :: has_unit_metric, omit_headers 139 TYPE(admm_type), POINTER :: admm_env 140 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 141 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 142 TYPE(cp_logger_type), POINTER :: logger 143 TYPE(cp_para_env_type), POINTER :: para_env 144 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit, matrix_p_mp2, matrix_s, & 145 matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, matrix_w, matrix_w_mp2, rho_ao 146 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_w_kp 147 TYPE(dft_control_type), POINTER :: dft_control 148 TYPE(lri_environment_type), POINTER :: lri_env 149 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_aux_fit 150 TYPE(mo_set_type), POINTER :: mo_set 151 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 152 TYPE(qs_energy_type), POINTER :: energy 153 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 154 TYPE(qs_ks_env_type), POINTER :: ks_env 155 TYPE(qs_rho_type), POINTER :: rho 156 TYPE(qs_subsys_type), POINTER :: subsys 157 TYPE(section_vals_type), POINTER :: print_section 158 TYPE(virial_type), POINTER :: virial 159 160 CALL timeset(routineN, handle) 161 NULLIFY (logger) 162 logger => cp_get_default_logger() 163 164 ! rebuild plane wave environment 165 CALL qs_env_rebuild_pw_env(qs_env) 166 167 ! zero out the forces in particle set 168 CALL get_qs_env(qs_env, particle_set=particle_set) 169 natom = SIZE(particle_set) 170 DO iatom = 1, natom 171 particle_set(iatom)%f = 0.0_dp 172 END DO 173 174 ! get atom mapping 175 NULLIFY (atomic_kind_set) 176 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set) 177 ALLOCATE (atom_of_kind(natom), kind_of(natom)) 178 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 179 atom_of_kind=atom_of_kind, & 180 kind_of=kind_of) 181 182 NULLIFY (force, subsys, dft_control) 183 CALL get_qs_env(qs_env, & 184 force=force, & 185 subsys=subsys, & 186 dft_control=dft_control) 187 IF (.NOT. ASSOCIATED(force)) THEN 188 ! *** Allocate the force data structure *** 189 nkind = SIZE(atomic_kind_set) 190 ALLOCATE (natom_of_kind(nkind)) 191 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 192 natom_of_kind=natom_of_kind) 193 CALL allocate_qs_force(force, natom_of_kind) 194 DEALLOCATE (natom_of_kind) 195 CALL qs_subsys_set(subsys, force=force) 196 END IF 197 CALL zero_qs_force(force) 198 199 ! Check if CDFT potential is needed and save it until forces have been calculated 200 IF (dft_control%qs_control%cdft) THEN 201 dft_control%qs_control%cdft_control%save_pot = .TRUE. 202 END IF 203 204 ! Set parameter for P screening with MP2 205 IF (ASSOCIATED(qs_env%mp2_env)) qs_env%mp2_env%not_last_hfx = .TRUE. 206 207 ! recalculate energy with forces 208 CALL qs_energies(qs_env, calc_forces=.TRUE.) 209 210 NULLIFY (para_env) 211 CALL get_qs_env(qs_env, & 212 para_env=para_env) 213 214 ! Now we handle some special cases 215 ! Maybe some of these would be better dealt with in qs_energies? 216 IF (qs_env%run_rtp) THEN 217 NULLIFY (matrix_w, matrix_s, ks_env) 218 CALL get_qs_env(qs_env, & 219 ks_env=ks_env, & 220 matrix_w=matrix_w, & 221 matrix_s=matrix_s) 222 CALL dbcsr_allocate_matrix_set(matrix_w, dft_control%nspins) 223 DO ispin = 1, dft_control%nspins 224 ALLOCATE (matrix_w(ispin)%matrix) 225 CALL dbcsr_copy(matrix_w(ispin)%matrix, matrix_s(1)%matrix, & 226 name="W MATRIX") 227 CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp) 228 END DO 229 CALL set_ks_env(ks_env, matrix_w=matrix_w) 230 231 CALL calc_c_mat_force(qs_env) 232 IF (dft_control%do_admm) CALL rt_admm_force(qs_env) 233 END IF 234 ! from an eventual Mulliken restraint 235 IF (dft_control%qs_control%mulliken_restraint) THEN 236 NULLIFY (matrix_w, matrix_s, rho) 237 CALL get_qs_env(qs_env, & 238 matrix_w=matrix_w, & 239 matrix_s=matrix_s, & 240 rho=rho) 241 NULLIFY (rho_ao) 242 CALL qs_rho_get(rho, rho_ao=rho_ao) 243 CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, & 244 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 IF (dft_control%dft_plus_u) THEN 249 NULLIFY (matrix_w) 250 CALL get_qs_env(qs_env, matrix_w=matrix_w) 251 CALL plus_u(qs_env=qs_env, matrix_w=matrix_w) 252 END IF 253 254 ! Write W Matrix to output (if requested) 255 CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric) 256 IF (.NOT. has_unit_metric) THEN 257 NULLIFY (matrix_w_kp) 258 CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp) 259 nspin = SIZE(matrix_w_kp, 1) 260 DO ispin = 1, nspin 261 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 262 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 extension=".Log") 265 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 266 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers) 267 after = MIN(MAX(after, 1), 16) 268 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 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 "DFT%PRINT%AO_MATRICES/W_MATRIX") 274 END IF 275 END DO 276 ENDIF 277 278 ! Compute core forces (also overwrites matrix_w) 279 IF (dft_control%qs_control%semi_empirical) THEN 280 CALL build_se_core_matrix(qs_env=qs_env, para_env=para_env, & 281 calculate_forces=.TRUE.) 282 CALL se_core_core_interaction(qs_env, para_env, calculate_forces=.TRUE.) 283 ELSEIF (dft_control%qs_control%dftb) THEN 284 CALL build_dftb_matrices(qs_env=qs_env, para_env=para_env, & 285 calculate_forces=.TRUE.) 286 CALL calculate_dftb_dispersion(qs_env=qs_env, para_env=para_env, & 287 calculate_forces=.TRUE.) 288 ELSEIF (dft_control%qs_control%xtb) THEN 289 CALL build_xtb_matrices(qs_env=qs_env, para_env=para_env, & 290 calculate_forces=.TRUE.) 291 ! Dispersion energy and forces are calculated in qs_energy? 292 ELSE 293 CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.TRUE.) 294 CALL calculate_ecore_self(qs_env) 295 CALL calculate_ecore_overlap(qs_env, para_env, calculate_forces=.TRUE.) 296 CALL calculate_ecore_efield(qs_env, calculate_forces=.TRUE.) 297 !swap external_e_potential before external_c_potential, to ensure 298 !that external potential on grid is loaded before calculating energy of cores 299 CALL external_e_potential(qs_env) 300 IF (.NOT. dft_control%qs_control%gapw) THEN 301 CALL external_c_potential(qs_env, calculate_forces=.TRUE.) 302 END IF 303 ! RIGPW matrices 304 IF (dft_control%qs_control%rigpw) THEN 305 CALL get_qs_env(qs_env=qs_env, lri_env=lri_env) 306 CALL build_ri_matrices(lri_env, qs_env, calculate_forces=.TRUE.) 307 ENDIF 308 END IF 309 310 ! Compute grid-based forces 311 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.TRUE.) 312 313 ! MP2 Code 314 IF (ASSOCIATED(qs_env%mp2_env)) THEN 315 NULLIFY (matrix_p_mp2, matrix_w_mp2, rho, ks_env, energy) 316 CALL get_qs_env(qs_env, & 317 matrix_p_mp2=matrix_p_mp2, & 318 matrix_w_mp2=matrix_w_mp2, & 319 ks_env=ks_env, & 320 rho=rho, & 321 energy=energy) 322 NULLIFY (rho_ao) 323 CALL qs_rho_get(rho, rho_ao=rho_ao) 324 325 ! with MP2 we have to recalculate the SCF energy with the 326 ! correct density 327 DO ispin = 1, dft_control%nspins 328 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp) 329 END DO 330 CALL qs_rho_update_rho(rho, qs_env=qs_env) 331 CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.) 332 CALL qs_ks_update_qs_env(qs_env, just_energy=.TRUE.) 333 energy%total = energy%total + energy%mp2 334 335 ! Compute MP2 properties 336 ! Get the HF+MP2 density 337 DO ispin = 1, dft_control%nspins 338 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp) 339 END DO 340 CALL qs_rho_update_rho(rho, qs_env=qs_env) 341 CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.) 342 CALL qs_scf_compute_properties(qs_env, wf_type='MP2 ') 343 ! Get everything back 344 DO ispin = 1, dft_control%nspins 345 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp) 346 END DO 347 CALL qs_rho_update_rho(rho, qs_env=qs_env) 348 CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.) 349 350 ! deallocate mp2_W 351 CALL dbcsr_deallocate_matrix_set(matrix_w_mp2) 352 CALL set_ks_env(ks_env, matrix_w_mp2=Null()) 353 354 END IF 355 356 ! Add forces resulting from wavefunction fitting 357 IF (dft_control%do_admm_dm) THEN 358 CPABORT("Forces with ADMM DM methods not implemented") 359 END IF 360 IF (dft_control%do_admm_mo .AND. .NOT. qs_env%run_rtp) THEN 361 NULLIFY (matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, matrix_ks_aux_fit, & 362 mos_aux_fit, mos, admm_env) 363 CALL get_qs_env(qs_env=qs_env, & 364 matrix_s_aux_fit=matrix_s_aux_fit, & 365 matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, & 366 matrix_ks_aux_fit=matrix_ks_aux_fit, & 367 mos_aux_fit=mos_aux_fit, & 368 mos=mos, & 369 admm_env=admm_env) 370 DO ispin = 1, dft_control%nspins 371 mo_set => mos(ispin)%mo_set 372 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff) 373 ! if no purification we need to calculate the H matrix for forces 374 IF (admm_env%purification_method == do_admm_purify_none) THEN 375 CALL get_mo_set(mo_set=mos_aux_fit(ispin)%mo_set, mo_coeff=mo_coeff_aux_fit) 376 CALL calc_aux_mo_derivs_none(ispin, qs_env%admm_env, mo_set, & 377 mo_coeff_aux_fit, matrix_ks_aux_fit) 378 END IF 379 END DO 380 CALL calc_mixed_overlap_force(qs_env) 381 END IF 382 383 ! *** replicate forces *** 384 CALL replicate_qs_force(force, para_env) 385 386 DO iatom = 1, natom 387 ikind = kind_of(iatom) 388 i = atom_of_kind(iatom) 389 ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 390 ! the force is - dE/dR, what is called force is actually the gradient 391 ! Things should have the right name 392 ! The minus sign below is a hack 393 ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 394 force(ikind)%other(1:3, i) = -particle_set(iatom)%f(1:3) + force(ikind)%ch_pulay(1:3, i) 395 force(ikind)%total(1:3, i) = force(ikind)%total(1:3, i) + force(ikind)%other(1:3, i) 396 particle_set(iatom)%f = -force(ikind)%total(1:3, i) 397 END DO 398 399 NULLIFY (virial, energy) 400 CALL get_qs_env(qs_env=qs_env, virial=virial, energy=energy) 401 ! *** distribute virial *** 402 IF (virial%pv_availability) THEN 403 CALL mp_sum(virial%pv_virial, para_env%group) 404 ! *** add the volume terms of the virial *** 405 IF ((.NOT. virial%pv_numer) .AND. & 406 (.NOT. (dft_control%qs_control%dftb .OR. & 407 dft_control%qs_control%xtb .OR. & 408 dft_control%qs_control%semi_empirical))) THEN 409 DO dir = 1, 3 410 virial%pv_virial(dir, dir) = virial%pv_virial(dir, dir) - energy%exc & 411 - 2.0_dp*energy%hartree 412 IF (dft_control%do_admm) THEN 413 virial%pv_virial(dir, dir) = virial%pv_virial(dir, dir) - energy%exc_aux_fit 414 END IF 415 ! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve. 416 ! The sign in pw_poisson_solve is correct for FIST, but not for QS. 417 ! There should be a more elegant solution to that... 418 END DO 419 END IF 420 END IF 421 422 output_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%DERIVATIVES", & 423 extension=".Log") 424 print_section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%DERIVATIVES") 425 IF (dft_control%qs_control%semi_empirical) THEN 426 CALL write_forces(force, atomic_kind_set, 2, output_unit=output_unit, & 427 print_section=print_section) 428 ELSE IF (dft_control%qs_control%dftb) THEN 429 CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, & 430 print_section=print_section) 431 ELSE IF (dft_control%qs_control%xtb) THEN 432 CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, & 433 print_section=print_section) 434 ELSE IF (dft_control%qs_control%gapw) THEN 435 CALL write_forces(force, atomic_kind_set, 1, output_unit=output_unit, & 436 print_section=print_section) 437 ELSE 438 CALL write_forces(force, atomic_kind_set, 0, output_unit=output_unit, & 439 print_section=print_section) 440 END IF 441 CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, & 442 "DFT%PRINT%DERIVATIVES") 443 444 ! deallocate W Matrix: 445 NULLIFY (ks_env, matrix_w_kp) 446 CALL get_qs_env(qs_env=qs_env, & 447 matrix_w_kp=matrix_w_kp, & 448 ks_env=ks_env) 449 CALL dbcsr_deallocate_matrix_set(matrix_w_kp) 450 CALL set_ks_env(ks_env, matrix_w=Null(), matrix_w_kp=Null()) 451 452 DEALLOCATE (atom_of_kind, kind_of) 453 454 CALL timestop(handle) 455 456 END SUBROUTINE qs_forces 457 458! ************************************************************************************************** 459!> \brief Write a Quickstep force data structure to output unit 460!> \param qs_force ... 461!> \param atomic_kind_set ... 462!> \param ftype ... 463!> \param output_unit ... 464!> \param print_section ... 465!> \date 05.06.2002 466!> \author MK 467!> \version 1.0 468! ************************************************************************************************** 469 SUBROUTINE write_forces(qs_force, atomic_kind_set, ftype, output_unit, & 470 print_section) 471 472 TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force 473 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 474 INTEGER, INTENT(IN) :: ftype, output_unit 475 TYPE(section_vals_type), POINTER :: print_section 476 477 CHARACTER(len=*), PARAMETER :: routineN = 'write_forces', routineP = moduleN//':'//routineN 478 479 CHARACTER(LEN=13) :: fmtstr5 480 CHARACTER(LEN=15) :: fmtstr4 481 CHARACTER(LEN=20) :: fmtstr3 482 CHARACTER(LEN=35) :: fmtstr2 483 CHARACTER(LEN=48) :: fmtstr1 484 INTEGER :: i, iatom, ikind, my_ftype, natom, ndigits 485 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of 486 REAL(KIND=dp), DIMENSION(3) :: grand_total 487 488 IF (output_unit > 0) THEN 489 490 IF (.NOT. ASSOCIATED(qs_force)) THEN 491 CALL cp_abort(__LOCATION__, & 492 "The qs_force pointer is not associated "// & 493 "and cannot be printed") 494 END IF 495 496 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 497 natom=natom) 498 ALLOCATE (atom_of_kind(natom)) 499 ALLOCATE (kind_of(natom)) 500 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 501 atom_of_kind=atom_of_kind, & 502 kind_of=kind_of) 503 504 ! Variable precision output of the forces 505 CALL section_vals_val_get(print_section, "NDIGITS", & 506 i_val=ndigits) 507 508 fmtstr1 = "(/,/,T2,A,/,/,T3,A,T11,A,T23,A,T40,A1,2( X,A1))" 509 WRITE (UNIT=fmtstr1(41:42), FMT="(I2)") ndigits + 5 510 511 fmtstr2 = "(/,(T2,I5,4X,I4,T18,A,T34,3F . ))" 512 WRITE (UNIT=fmtstr2(32:33), FMT="(I2)") ndigits 513 WRITE (UNIT=fmtstr2(29:30), FMT="(I2)") ndigits + 6 514 515 fmtstr3 = "(/,T3,A,T34,3F . )" 516 WRITE (UNIT=fmtstr3(18:19), FMT="(I2)") ndigits 517 WRITE (UNIT=fmtstr3(15:16), FMT="(I2)") ndigits + 6 518 519 fmtstr4 = "((T34,3F . ))" 520 WRITE (UNIT=fmtstr4(12:13), FMT="(I2)") ndigits 521 WRITE (UNIT=fmtstr4(9:10), FMT="(I2)") ndigits + 6 522 523 fmtstr5 = "(/T2,A//T3,A)" 524 525 WRITE (UNIT=output_unit, FMT=fmtstr1) & 526 "FORCES [a.u.]", "Atom", "Kind", "Component", "X", "Y", "Z" 527 528 grand_total(:) = 0.0_dp 529 530 my_ftype = ftype 531 532 SELECT CASE (my_ftype) 533 CASE DEFAULT 534 DO iatom = 1, natom 535 ikind = kind_of(iatom) 536 i = atom_of_kind(iatom) 537 WRITE (UNIT=output_unit, FMT=fmtstr2) & 538 iatom, ikind, " total", qs_force(ikind)%total(1:3, i) 539 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i) 540 END DO 541 CASE (0) 542 DO iatom = 1, natom 543 ikind = kind_of(iatom) 544 i = atom_of_kind(iatom) 545 WRITE (UNIT=output_unit, FMT=fmtstr2) & 546 iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), & 547 iatom, ikind, " overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), & 548 iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), & 549 iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), & 550 iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), & 551 iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), & 552 iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), & 553 iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), & 554 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), & 555 iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), & 556 iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), & 557 iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), & 558 iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), & 559 iatom, ikind, " other", qs_force(ikind)%other(1:3, i), & 560 iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), & 561 iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), & 562 iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), & 563 iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), & 564 iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), & 565 iatom, ikind, " mp2_sep", qs_force(ikind)%mp2_sep(1:3, i), & 566 iatom, ikind, " total", qs_force(ikind)%total(1:3, i) 567 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i) 568 END DO 569 CASE (1) 570 DO iatom = 1, natom 571 ikind = kind_of(iatom) 572 i = atom_of_kind(iatom) 573 WRITE (UNIT=output_unit, FMT=fmtstr2) & 574 iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), & 575 iatom, ikind, " overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), & 576 iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), & 577 iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), & 578 iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), & 579 iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), & 580 iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), & 581 iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), & 582 iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), & 583 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), & 584 iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), & 585 iatom, ikind, " vhxc_atom", qs_force(ikind)%vhxc_atom(1:3, i), & 586 iatom, ikind, " g0s_Vh_elec", qs_force(ikind)%g0s_Vh_elec(1:3, i), & 587 iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), & 588 iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), & 589 iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), & 590 iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), & 591 iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), & 592 iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), & 593 iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), & 594 iatom, ikind, " mp2_sep", qs_force(ikind)%mp2_sep(1:3, i), & 595 iatom, ikind, " total", qs_force(ikind)%total(1:3, i) 596 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i) 597 END DO 598 CASE (2) 599 DO iatom = 1, natom 600 ikind = kind_of(iatom) 601 i = atom_of_kind(iatom) 602 WRITE (UNIT=output_unit, FMT=fmtstr2) & 603 iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), & 604 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), & 605 iatom, ikind, " total", qs_force(ikind)%total(1:3, i) 606 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i) 607 END DO 608 CASE (3) 609 DO iatom = 1, natom 610 ikind = kind_of(iatom) 611 i = atom_of_kind(iatom) 612 WRITE (UNIT=output_unit, FMT=fmtstr2) & 613 iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), & 614 iatom, ikind, "overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), & 615 iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), & 616 iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), & 617 iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), & 618 iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), & 619 iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), & 620 iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), & 621 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), & 622 iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), & 623 iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), & 624 iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), & 625 iatom, ikind, " mp2_sep", qs_force(ikind)%mp2_sep(1:3, i), & 626 iatom, ikind, " total", qs_force(ikind)%total(1:3, i) 627 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i) 628 END DO 629 CASE (4) 630 DO iatom = 1, natom 631 ikind = kind_of(iatom) 632 i = atom_of_kind(iatom) 633 WRITE (UNIT=output_unit, FMT=fmtstr2) & 634 iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), & 635 iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), & 636 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), & 637 iatom, ikind, " repulsive", qs_force(ikind)%repulsive(1:3, i), & 638 iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), & 639 iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), & 640 iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), & 641 iatom, ikind, " total", qs_force(ikind)%total(1:3, i) 642 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i) 643 END DO 644 CASE (5) 645 DO iatom = 1, natom 646 ikind = kind_of(iatom) 647 i = atom_of_kind(iatom) 648 WRITE (UNIT=output_unit, FMT=fmtstr2) & 649 iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), & 650 iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), & 651 iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), & 652 iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), & 653 iatom, ikind, " all potential", qs_force(ikind)%all_potential(1:3, i), & 654 iatom, ikind, " other", qs_force(ikind)%other(1:3, i), & 655 iatom, ikind, " total", qs_force(ikind)%total(1:3, i) 656 grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i) 657 END DO 658 END SELECT 659 660 WRITE (UNIT=output_unit, FMT=fmtstr3) "Sum of total", grand_total(1:3) 661 662 DEALLOCATE (atom_of_kind) 663 DEALLOCATE (kind_of) 664 665 END IF 666 667 END SUBROUTINE write_forces 668 669END MODULE qs_force 670