LCOV - code coverage report
Current view: top level - src - qs_force.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 80.6 % 319 257
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 3 3

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Quickstep force driver routine
      10              : !> \author MK (12.06.2002)
      11              : ! **************************************************************************************************
      12              : MODULE qs_force
      13              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      14              :                                               get_atomic_kind_set
      15              :    USE cp_control_types,                ONLY: dft_control_type
      16              :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      17              :                                               dbcsr_p_type,&
      18              :                                               dbcsr_set
      19              :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      20              :                                               dbcsr_deallocate_matrix_set
      21              :    USE cp_dbcsr_output,                 ONLY: cp_dbcsr_write_sparse_matrix
      22              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23              :                                               cp_logger_get_default_io_unit,&
      24              :                                               cp_logger_type
      25              :    USE cp_output_handling,              ONLY: cp_p_file,&
      26              :                                               cp_print_key_finished_output,&
      27              :                                               cp_print_key_should_output,&
      28              :                                               cp_print_key_unit_nr
      29              :    USE dft_plus_u,                      ONLY: plus_u
      30              :    USE ec_env_types,                    ONLY: energy_correction_type
      31              :    USE efield_utils,                    ONLY: calculate_ecore_efield,&
      32              :                                               efield_potential_lengh_gauge
      33              :    USE energy_corrections,              ONLY: energy_correction
      34              :    USE excited_states,                  ONLY: excited_state_energy
      35              :    USE hfx_exx,                         ONLY: calculate_exx
      36              :    USE input_constants,                 ONLY: ri_mp2_laplace,&
      37              :                                               ri_mp2_method_gpw,&
      38              :                                               ri_rpa_method_gpw
      39              :    USE input_section_types,             ONLY: section_vals_get,&
      40              :                                               section_vals_get_subs_vals,&
      41              :                                               section_vals_type,&
      42              :                                               section_vals_val_get
      43              :    USE kinds,                           ONLY: dp
      44              :    USE lri_environment_types,           ONLY: lri_environment_type
      45              :    USE message_passing,                 ONLY: mp_para_env_type
      46              :    USE mp2_cphf,                        ONLY: update_mp2_forces
      47              :    USE mulliken,                        ONLY: mulliken_restraint
      48              :    USE particle_types,                  ONLY: particle_type
      49              :    USE qs_core_energies,                ONLY: calculate_ecore_overlap,&
      50              :                                               calculate_ecore_self
      51              :    USE qs_core_hamiltonian,             ONLY: build_core_hamiltonian_matrix
      52              :    USE qs_dftb_dispersion,              ONLY: calculate_dftb_dispersion
      53              :    USE qs_dftb_matrices,                ONLY: build_dftb_matrices
      54              :    USE qs_energy,                       ONLY: qs_energies
      55              :    USE qs_energy_types,                 ONLY: qs_energy_type
      56              :    USE qs_environment_methods,          ONLY: qs_env_rebuild_pw_env
      57              :    USE qs_environment_types,            ONLY: get_qs_env,&
      58              :                                               qs_environment_type
      59              :    USE qs_external_potential,           ONLY: external_c_potential,&
      60              :                                               external_e_potential
      61              :    USE qs_force_types,                  ONLY: allocate_qs_force,&
      62              :                                               qs_force_type,&
      63              :                                               replicate_qs_force,&
      64              :                                               zero_qs_force
      65              :    USE qs_ks_methods,                   ONLY: qs_ks_update_qs_env
      66              :    USE qs_ks_types,                     ONLY: qs_ks_env_type,&
      67              :                                               set_ks_env
      68              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      69              :                                               qs_rho_type
      70              :    USE qs_scf_post_scf,                 ONLY: qs_scf_compute_properties
      71              :    USE qs_subsys_types,                 ONLY: qs_subsys_set,&
      72              :                                               qs_subsys_type
      73              :    USE ri_environment_methods,          ONLY: build_ri_matrices
      74              :    USE rt_propagation_forces,           ONLY: calc_c_mat_force,&
      75              :                                               rt_admm_force
      76              :    USE rt_propagation_velocity_gauge,   ONLY: velocity_gauge_ks_matrix,&
      77              :                                               velocity_gauge_nl_force
      78              :    USE se_core_core,                    ONLY: se_core_core_interaction
      79              :    USE se_core_matrix,                  ONLY: build_se_core_matrix
      80              :    USE tblite_interface,                ONLY: build_tblite_matrices
      81              :    USE virial_types,                    ONLY: symmetrize_virial,&
      82              :                                               virial_type
      83              :    USE xtb_matrices,                    ONLY: build_xtb_matrices
      84              : #include "./base/base_uses.f90"
      85              : 
      86              :    IMPLICIT NONE
      87              : 
      88              :    PRIVATE
      89              : 
      90              : ! *** Global parameters ***
      91              : 
      92              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_force'
      93              : 
      94              : ! *** Public subroutines ***
      95              : 
      96              :    PUBLIC :: qs_calc_energy_force
      97              : 
      98              : CONTAINS
      99              : 
     100              : ! **************************************************************************************************
     101              : !> \brief ...
     102              : !> \param qs_env ...
     103              : !> \param calc_force ...
     104              : !> \param consistent_energies ...
     105              : !> \param linres ...
     106              : ! **************************************************************************************************
     107        22765 :    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        22765 :       qs_env%linres_run = linres
     112        22765 :       IF (calc_force) THEN
     113        10581 :          CALL qs_forces(qs_env)
     114              :       ELSE
     115              :          CALL qs_energies(qs_env, calc_forces=.FALSE., &
     116        12184 :                           consistent_energies=consistent_energies)
     117              :       END IF
     118              : 
     119        22765 :    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        10581 :    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        10581 :       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        10581 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     143              :       TYPE(cp_logger_type), POINTER                      :: logger
     144        10581 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, matrix_w, rho_ao
     145        10581 :       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        10581 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     151              :       TYPE(qs_energy_type), POINTER                      :: energy
     152        10581 :       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        10581 :       CALL timeset(routineN, handle)
     160        10581 :       NULLIFY (logger)
     161        10581 :       logger => cp_get_default_logger()
     162              : 
     163              :       ! rebuild plane wave environment
     164        10581 :       CALL qs_env_rebuild_pw_env(qs_env)
     165              : 
     166              :       ! zero out the forces in particle set
     167        10581 :       CALL get_qs_env(qs_env, particle_set=particle_set)
     168        10581 :       natom = SIZE(particle_set)
     169        77750 :       DO iatom = 1, natom
     170       279257 :          particle_set(iatom)%f = 0.0_dp
     171              :       END DO
     172              : 
     173              :       ! get atom mapping
     174        10581 :       NULLIFY (atomic_kind_set)
     175        10581 :       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        10581 :                                kind_of=kind_of)
     179              : 
     180        10581 :       NULLIFY (force, subsys, dft_control)
     181              :       CALL get_qs_env(qs_env, &
     182              :                       force=force, &
     183              :                       subsys=subsys, &
     184        10581 :                       dft_control=dft_control)
     185        10581 :       IF (.NOT. ASSOCIATED(force)) THEN
     186              :          !   *** Allocate the force data structure ***
     187         2895 :          nkind = SIZE(atomic_kind_set)
     188         2895 :          CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
     189         2895 :          CALL allocate_qs_force(force, natom_of_kind)
     190         2895 :          DEALLOCATE (natom_of_kind)
     191         2895 :          CALL qs_subsys_set(subsys, force=force)
     192              :       END IF
     193        10581 :       CALL zero_qs_force(force)
     194              : 
     195              :       ! Check if CDFT potential is needed and save it until forces have been calculated
     196        10581 :       IF (dft_control%qs_control%cdft) THEN
     197          118 :          dft_control%qs_control%cdft_control%save_pot = .TRUE.
     198              :       END IF
     199              : 
     200              :       ! recalculate energy with forces
     201        10581 :       CALL qs_energies(qs_env, calc_forces=.TRUE.)
     202              : 
     203        10581 :       NULLIFY (para_env)
     204              :       CALL get_qs_env(qs_env, &
     205        10581 :                       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        10581 :       IF (qs_env%run_rtp) THEN
     210         1216 :          NULLIFY (matrix_w, matrix_s, ks_env)
     211              :          CALL get_qs_env(qs_env, &
     212              :                          ks_env=ks_env, &
     213              :                          matrix_w=matrix_w, &
     214         1216 :                          matrix_s=matrix_s)
     215         1216 :          CALL dbcsr_allocate_matrix_set(matrix_w, dft_control%nspins)
     216         2674 :          DO ispin = 1, dft_control%nspins
     217         1458 :             ALLOCATE (matrix_w(ispin)%matrix)
     218              :             CALL dbcsr_copy(matrix_w(ispin)%matrix, matrix_s(1)%matrix, &
     219         1458 :                             name="W MATRIX")
     220         2674 :             CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
     221              :          END DO
     222         1216 :          CALL set_ks_env(ks_env, matrix_w=matrix_w)
     223              : 
     224         1216 :          CALL calc_c_mat_force(qs_env)
     225         1216 :          IF (dft_control%do_admm) CALL rt_admm_force(qs_env)
     226         1216 :          IF (dft_control%rtp_control%velocity_gauge .AND. dft_control%rtp_control%nl_gauge_transform) &
     227           22 :             CALL velocity_gauge_nl_force(qs_env, particle_set)
     228              :       END IF
     229              :       ! from an eventual Mulliken restraint
     230        10581 :       IF (dft_control%qs_control%mulliken_restraint) THEN
     231            6 :          NULLIFY (matrix_w, matrix_s, rho)
     232              :          CALL get_qs_env(qs_env, &
     233              :                          matrix_w=matrix_w, &
     234              :                          matrix_s=matrix_s, &
     235            6 :                          rho=rho)
     236            6 :          NULLIFY (rho_ao)
     237            6 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     238              :          CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
     239            6 :                                  para_env, matrix_s(1)%matrix, rho_ao, w_matrix=matrix_w)
     240              :       END IF
     241              :       ! Add non-Pulay contribution of DFT+U to W matrix, since it has also to be
     242              :       ! digested with overlap matrix derivatives
     243        10581 :       IF (dft_control%dft_plus_u) THEN
     244           72 :          NULLIFY (matrix_w)
     245           72 :          CALL get_qs_env(qs_env, matrix_w=matrix_w)
     246           72 :          CALL plus_u(qs_env=qs_env, matrix_w=matrix_w)
     247              :       END IF
     248              : 
     249              :       ! Write W Matrix to output (if requested)
     250        10581 :       CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
     251        10581 :       IF (.NOT. has_unit_metric) THEN
     252         7579 :          NULLIFY (matrix_w_kp)
     253         7579 :          CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp)
     254         7579 :          nspin = SIZE(matrix_w_kp, 1)
     255        16162 :          DO ispin = 1, nspin
     256         8583 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     257         7579 :                                                  qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX"), cp_p_file)) THEN
     258              :                iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX", &
     259            8 :                                          extension=".Log")
     260            8 :                CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
     261            8 :                CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
     262            8 :                after = MIN(MAX(after, 1), 16)
     263           16 :                DO ic = 1, SIZE(matrix_w_kp, 2)
     264              :                   CALL cp_dbcsr_write_sparse_matrix(matrix_w_kp(ispin, ic)%matrix, 4, after, qs_env, &
     265           16 :                                                     para_env, output_unit=iw, omit_headers=omit_headers)
     266              :                END DO
     267              :                CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
     268            8 :                                                  "DFT%PRINT%AO_MATRICES/W_MATRIX")
     269              :             END IF
     270              :          END DO
     271              :       END IF
     272              : 
     273              :       ! Check if energy correction should be skipped
     274        10581 :       perform_ec = .FALSE.
     275        10581 :       IF (qs_env%energy_correction) THEN
     276          436 :          CALL get_qs_env(qs_env, ec_env=ec_env)
     277          436 :          IF (.NOT. ec_env%do_skip) perform_ec = .TRUE.
     278              :       END IF
     279              : 
     280              :       ! Compute core forces (also overwrites matrix_w)
     281        10581 :       IF (dft_control%qs_control%semi_empirical) THEN
     282              :          CALL build_se_core_matrix(qs_env=qs_env, para_env=para_env, &
     283         3002 :                                    calculate_forces=.TRUE.)
     284         3002 :          CALL se_core_core_interaction(qs_env, para_env, calculate_forces=.TRUE.)
     285         7579 :       ELSEIF (dft_control%qs_control%dftb) THEN
     286              :          CALL build_dftb_matrices(qs_env=qs_env, para_env=para_env, &
     287          724 :                                   calculate_forces=.TRUE.)
     288              :          CALL calculate_dftb_dispersion(qs_env=qs_env, para_env=para_env, &
     289          724 :                                         calculate_forces=.TRUE.)
     290         6855 :       ELSEIF (dft_control%qs_control%xtb) THEN
     291          568 :          IF (dft_control%qs_control%xtb_control%do_tblite) THEN
     292           28 :             CALL build_tblite_matrices(qs_env=qs_env, calculate_forces=.TRUE.)
     293              :          ELSE
     294          540 :             CALL build_xtb_matrices(qs_env=qs_env, calculate_forces=.TRUE.)
     295              :          END IF
     296         6287 :       ELSEIF (perform_ec) THEN
     297              :          ! Calculates core and grid based forces
     298          436 :          CALL energy_correction(qs_env, ec_init=.FALSE., calculate_forces=.TRUE.)
     299              :       ELSE
     300              :          ! Dispersion energy and forces are calculated in qs_energy?
     301         5851 :          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         5851 :          IF (qs_env%run_rtp) THEN
     304          804 :             IF (dft_control%apply_efield_field) &
     305          160 :                CALL efield_potential_lengh_gauge(qs_env)
     306          804 :             IF (dft_control%rtp_control%velocity_gauge) &
     307           22 :                CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.FALSE.)
     308              : 
     309              :          END IF
     310         5851 :          CALL calculate_ecore_self(qs_env)
     311         5851 :          CALL calculate_ecore_overlap(qs_env, para_env, calculate_forces=.TRUE.)
     312         5851 :          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         5851 :          CALL external_e_potential(qs_env)
     316         5851 :          IF (.NOT. dft_control%qs_control%gapw) THEN
     317         5419 :             CALL external_c_potential(qs_env, calculate_forces=.TRUE.)
     318              :          END IF
     319              :          ! RIGPW  matrices
     320         5851 :          IF (dft_control%qs_control%rigpw) THEN
     321            0 :             CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
     322            0 :             CALL build_ri_matrices(lri_env, qs_env, calculate_forces=.TRUE.)
     323              :          END IF
     324              :       END IF
     325              : 
     326              :       ! MP2 Code
     327        10581 :       IF (ASSOCIATED(qs_env%mp2_env)) THEN
     328          322 :          NULLIFY (energy)
     329          322 :          CALL get_qs_env(qs_env, energy=energy)
     330          322 :          CALL qs_scf_compute_properties(qs_env, wf_type='MP2   ', do_mp2=.TRUE.)
     331          322 :          CALL qs_ks_update_qs_env(qs_env, just_energy=.TRUE.)
     332          322 :          energy%total = energy%total + energy%mp2
     333              : 
     334              :          IF ((qs_env%mp2_env%method == ri_mp2_method_gpw .OR. qs_env%mp2_env%method == ri_mp2_laplace .OR. &
     335              :               qs_env%mp2_env%method == ri_rpa_method_gpw) &
     336          322 :              .AND. .NOT. qs_env%mp2_env%do_im_time) THEN
     337          272 :             CALL update_mp2_forces(qs_env)
     338              :          END IF
     339              : 
     340              :          !RPA EXX energy and forces
     341          322 :          IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
     342              :             do_exx = .FALSE.
     343           52 :             hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
     344           52 :             CALL section_vals_get(hfx_sections, explicit=do_exx)
     345           52 :             IF (do_exx) THEN
     346           26 :                do_gw = qs_env%mp2_env%ri_rpa%do_ri_g0w0
     347           26 :                do_admm = qs_env%mp2_env%ri_rpa%do_admm
     348           26 :                reuse_hfx = qs_env%mp2_env%ri_rpa%reuse_hfx
     349           26 :                do_im_time = qs_env%mp2_env%do_im_time
     350           26 :                output_unit = cp_logger_get_default_io_unit()
     351           26 :                dummy_real = 0.0_dp
     352              : 
     353              :                CALL calculate_exx(qs_env=qs_env, &
     354              :                                   unit_nr=output_unit, &
     355              :                                   hfx_sections=hfx_sections, &
     356              :                                   x_data=qs_env%mp2_env%ri_rpa%x_data, &
     357              :                                   do_gw=do_gw, &
     358              :                                   do_admm=do_admm, &
     359              :                                   calc_forces=.TRUE., &
     360              :                                   reuse_hfx=reuse_hfx, &
     361              :                                   do_im_time=do_im_time, &
     362              :                                   E_ex_from_GW=dummy_real, &
     363              :                                   E_admm_from_GW=dummy_real2, &
     364           26 :                                   t3=dummy_real)
     365              :             END IF
     366              :          END IF
     367        10259 :       ELSEIF (perform_ec) THEN
     368              :          ! energy correction forces postponed
     369         9823 :       ELSEIF (qs_env%harris_method) THEN
     370              :          ! Harris method forces already done in harris_energy_correction
     371              :       ELSE
     372              :          ! Compute grid-based forces
     373         9819 :          CALL qs_ks_update_qs_env(qs_env, calculate_forces=.TRUE.)
     374              :       END IF
     375              : 
     376              :       ! Excited state forces
     377        10581 :       CALL excited_state_energy(qs_env, calculate_forces=.TRUE.)
     378              : 
     379              :       ! replicate forces (get current pointer)
     380        10581 :       NULLIFY (force)
     381        10581 :       CALL get_qs_env(qs_env=qs_env, force=force)
     382        10581 :       CALL replicate_qs_force(force, para_env)
     383              : 
     384        77750 :       DO iatom = 1, natom
     385        67169 :          ikind = kind_of(iatom)
     386        67169 :          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       537352 :          force(ikind)%other(1:3, i) = -particle_set(iatom)%f(1:3) + force(ikind)%ch_pulay(1:3, i)
     393       268676 :          force(ikind)%total(1:3, i) = force(ikind)%total(1:3, i) + force(ikind)%other(1:3, i)
     394       547933 :          particle_set(iatom)%f = -force(ikind)%total(1:3, i)
     395              :       END DO
     396              : 
     397        10581 :       NULLIFY (virial, energy)
     398        10581 :       CALL get_qs_env(qs_env=qs_env, virial=virial, energy=energy)
     399        10581 :       IF (virial%pv_availability) THEN
     400          912 :          CALL para_env%sum(virial%pv_overlap)
     401          912 :          CALL para_env%sum(virial%pv_ekinetic)
     402          912 :          CALL para_env%sum(virial%pv_ppl)
     403          912 :          CALL para_env%sum(virial%pv_ppnl)
     404          912 :          CALL para_env%sum(virial%pv_ecore_overlap)
     405          912 :          CALL para_env%sum(virial%pv_ehartree)
     406          912 :          CALL para_env%sum(virial%pv_exc)
     407          912 :          CALL para_env%sum(virial%pv_exx)
     408          912 :          CALL para_env%sum(virial%pv_vdw)
     409          912 :          CALL para_env%sum(virial%pv_mp2)
     410          912 :          CALL para_env%sum(virial%pv_nlcc)
     411          912 :          CALL para_env%sum(virial%pv_gapw)
     412          912 :          CALL para_env%sum(virial%pv_lrigpw)
     413          912 :          CALL para_env%sum(virial%pv_virial)
     414          912 :          CALL symmetrize_virial(virial)
     415              :          ! Add the volume terms of the virial
     416          912 :          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          622 :             IF (perform_ec) THEN
     425          168 :                CALL get_qs_env(qs_env, ec_env=ec_env)
     426          168 :                energy%hartree = ec_env%ehartree
     427          168 :                energy%exc = ec_env%exc
     428          168 :                IF (dft_control%do_admm) THEN
     429           38 :                   energy%exc_aux_fit = ec_env%exc_aux_fit
     430              :                END IF
     431              :             END IF
     432         2488 :             DO i = 1, 3
     433              :                virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) &
     434         1866 :                                           - 2.0_dp*(energy%hartree + energy%sccs_pol)
     435              :                virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc &
     436         1866 :                                         - 2.0_dp*(energy%hartree + energy%sccs_pol)
     437         1866 :                virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc
     438         2488 :                IF (dft_control%do_admm) THEN
     439          198 :                   virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc_aux_fit
     440          198 :                   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        10581 :                                          extension=".Log")
     451        10581 :       print_section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%DERIVATIVES")
     452        10581 :       IF (dft_control%qs_control%semi_empirical) THEN
     453              :          CALL write_forces(force, atomic_kind_set, 2, output_unit=output_unit, &
     454         3002 :                            print_section=print_section)
     455         7579 :       ELSE IF (dft_control%qs_control%dftb) THEN
     456              :          CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, &
     457          724 :                            print_section=print_section)
     458         6855 :       ELSE IF (dft_control%qs_control%xtb) THEN
     459              :          CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, &
     460          568 :                            print_section=print_section)
     461         6287 :       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          504 :                            print_section=print_section)
     464              :       ELSE
     465              :          CALL write_forces(force, atomic_kind_set, 0, output_unit=output_unit, &
     466         5783 :                            print_section=print_section)
     467              :       END IF
     468              :       CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, &
     469        10581 :                                         "DFT%PRINT%DERIVATIVES")
     470              : 
     471              :       ! deallocate W Matrix:
     472        10581 :       NULLIFY (ks_env, matrix_w_kp)
     473              :       CALL get_qs_env(qs_env=qs_env, &
     474              :                       matrix_w_kp=matrix_w_kp, &
     475        10581 :                       ks_env=ks_env)
     476        10581 :       CALL dbcsr_deallocate_matrix_set(matrix_w_kp)
     477        10581 :       NULLIFY (matrix_w_kp)
     478        10581 :       CALL set_ks_env(ks_env, matrix_w_kp=matrix_w_kp)
     479              : 
     480        10581 :       DEALLOCATE (atom_of_kind, kind_of)
     481              : 
     482        10581 :       CALL timestop(handle)
     483              : 
     484        21162 :    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        10581 :    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        10581 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     512              :       REAL(KIND=dp), DIMENSION(3)                        :: grand_total
     513              : 
     514        10581 :       IF (output_unit > 0) THEN
     515              : 
     516          156 :          IF (.NOT. ASSOCIATED(qs_force)) THEN
     517              :             CALL cp_abort(__LOCATION__, &
     518              :                           "The qs_force pointer is not associated "// &
     519            0 :                           "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          156 :                                   kind_of=kind_of, natom=natom)
     524              : 
     525              :          ! Variable precision output of the forces
     526              :          CALL section_vals_val_get(print_section, "NDIGITS", &
     527          156 :                                    i_val=ndigits)
     528              : 
     529          156 :          fmtstr1 = "(/,/,T2,A,/,/,T3,A,T11,A,T23,A,T40,A1,2(  X,A1))"
     530          156 :          WRITE (UNIT=fmtstr1(41:42), FMT="(I2)") ndigits + 5
     531              : 
     532          156 :          fmtstr2 = "(/,(T2,I5,4X,I4,T18,A,T34,3F  .  ))"
     533          156 :          WRITE (UNIT=fmtstr2(32:33), FMT="(I2)") ndigits
     534          156 :          WRITE (UNIT=fmtstr2(29:30), FMT="(I2)") ndigits + 6
     535              : 
     536          156 :          fmtstr3 = "(/,T3,A,T34,3F  .  )"
     537          156 :          WRITE (UNIT=fmtstr3(18:19), FMT="(I2)") ndigits
     538          156 :          WRITE (UNIT=fmtstr3(15:16), FMT="(I2)") ndigits + 6
     539              : 
     540          156 :          fmtstr4 = "((T34,3F  .  ))"
     541          156 :          WRITE (UNIT=fmtstr4(12:13), FMT="(I2)") ndigits
     542          156 :          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          156 :             "FORCES [a.u.]", "Atom", "Kind", "Component", "X", "Y", "Z"
     548              : 
     549          156 :          grand_total(:) = 0.0_dp
     550              : 
     551          156 :          my_ftype = ftype
     552              : 
     553            0 :          SELECT CASE (my_ftype)
     554              :          CASE DEFAULT
     555            0 :             DO iatom = 1, natom
     556            0 :                ikind = kind_of(iatom)
     557            0 :                i = atom_of_kind(iatom)
     558              :                WRITE (UNIT=output_unit, FMT=fmtstr2) &
     559            0 :                   iatom, ikind, "         total", qs_force(ikind)%total(1:3, i)
     560            0 :                grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
     561              :             END DO
     562              :          CASE (0)
     563          476 :             DO iatom = 1, natom
     564          342 :                ikind = kind_of(iatom)
     565          342 :                i = atom_of_kind(iatom)
     566              :                WRITE (UNIT=output_unit, FMT=fmtstr2) &
     567         1368 :                   iatom, ikind, "       overlap", qs_force(ikind)%overlap(1:3, i), &
     568         1368 :                   iatom, ikind, "  overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
     569         1368 :                   iatom, ikind, "       kinetic", qs_force(ikind)%kinetic(1:3, i), &
     570         1368 :                   iatom, ikind, "       gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
     571         1368 :                   iatom, ikind, "      gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
     572         1368 :                   iatom, ikind, "      gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
     573         1368 :                   iatom, ikind, "  core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
     574         1368 :                   iatom, ikind, "      rho_core", qs_force(ikind)%rho_core(1:3, i), &
     575         1368 :                   iatom, ikind, "      rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
     576         1368 :                   iatom, ikind, "  rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
     577         1368 :                   iatom, ikind, "      ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
     578         1368 :                   iatom, ikind, "    dispersion", qs_force(ikind)%dispersion(1:3, i), &
     579         1368 :                   iatom, ikind, "           gCP", qs_force(ikind)%gcp(1:3, i), &
     580         1368 :                   iatom, ikind, "         other", qs_force(ikind)%other(1:3, i), &
     581         1368 :                   iatom, ikind, "       fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
     582         1368 :                   iatom, ikind, "     ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
     583         1368 :                   iatom, ikind, "        efield", qs_force(ikind)%efield(1:3, i), &
     584         1368 :                   iatom, ikind, "           eev", qs_force(ikind)%eev(1:3, i), &
     585         1368 :                   iatom, ikind, "   mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
     586         1710 :                   iatom, ikind, "         total", qs_force(ikind)%total(1:3, i)
     587         1502 :                grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
     588              :             END DO
     589              :          CASE (1)
     590            0 :             DO iatom = 1, natom
     591            0 :                ikind = kind_of(iatom)
     592            0 :                i = atom_of_kind(iatom)
     593              :                WRITE (UNIT=output_unit, FMT=fmtstr2) &
     594            0 :                   iatom, ikind, "       overlap", qs_force(ikind)%overlap(1:3, i), &
     595            0 :                   iatom, ikind, "  overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
     596            0 :                   iatom, ikind, "       kinetic", qs_force(ikind)%kinetic(1:3, i), &
     597            0 :                   iatom, ikind, "       gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
     598            0 :                   iatom, ikind, "      gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
     599            0 :                   iatom, ikind, "      gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
     600            0 :                   iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
     601            0 :                   iatom, ikind, "  core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
     602            0 :                   iatom, ikind, "      rho_core", qs_force(ikind)%rho_core(1:3, i), &
     603            0 :                   iatom, ikind, "      rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
     604            0 :                   iatom, ikind, "  rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
     605            0 :                   iatom, ikind, "     vhxc_atom", qs_force(ikind)%vhxc_atom(1:3, i), &
     606            0 :                   iatom, ikind, "   g0s_Vh_elec", qs_force(ikind)%g0s_Vh_elec(1:3, i), &
     607            0 :                   iatom, ikind, "      ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
     608            0 :                   iatom, ikind, "    dispersion", qs_force(ikind)%dispersion(1:3, i), &
     609            0 :                   iatom, ikind, "           gCP", qs_force(ikind)%gcp(1:3, i), &
     610            0 :                   iatom, ikind, "       fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
     611            0 :                   iatom, ikind, "     ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
     612            0 :                   iatom, ikind, "        efield", qs_force(ikind)%efield(1:3, i), &
     613            0 :                   iatom, ikind, "           eev", qs_force(ikind)%eev(1:3, i), &
     614            0 :                   iatom, ikind, "   mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
     615            0 :                   iatom, ikind, "         total", qs_force(ikind)%total(1:3, i)
     616            0 :                grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
     617              :             END DO
     618              :          CASE (2)
     619           75 :             DO iatom = 1, natom
     620           73 :                ikind = kind_of(iatom)
     621           73 :                i = atom_of_kind(iatom)
     622              :                WRITE (UNIT=output_unit, FMT=fmtstr2) &
     623          292 :                   iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
     624          292 :                   iatom, ikind, "      rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
     625          365 :                   iatom, ikind, "         total", qs_force(ikind)%total(1:3, i)
     626          294 :                grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
     627              :             END DO
     628              :          CASE (3)
     629            0 :             DO iatom = 1, natom
     630            0 :                ikind = kind_of(iatom)
     631            0 :                i = atom_of_kind(iatom)
     632              :                WRITE (UNIT=output_unit, FMT=fmtstr2) &
     633            0 :                   iatom, ikind, "        overlap", qs_force(ikind)%overlap(1:3, i), &
     634            0 :                   iatom, ikind, "overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
     635            0 :                   iatom, ikind, "        kinetic", qs_force(ikind)%kinetic(1:3, i), &
     636            0 :                   iatom, ikind, "        gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
     637            0 :                   iatom, ikind, "       gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
     638            0 :                   iatom, ikind, "       gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
     639            0 :                   iatom, ikind, "   core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
     640            0 :                   iatom, ikind, "       rho_core", qs_force(ikind)%rho_core(1:3, i), &
     641            0 :                   iatom, ikind, "       rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
     642            0 :                   iatom, ikind, "       ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
     643            0 :                   iatom, ikind, "        fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
     644            0 :                   iatom, ikind, "   mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
     645            0 :                   iatom, ikind, "          total", qs_force(ikind)%total(1:3, i)
     646            0 :                grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
     647              :             END DO
     648              :          CASE (4)
     649          160 :             DO iatom = 1, natom
     650          140 :                ikind = kind_of(iatom)
     651          140 :                i = atom_of_kind(iatom)
     652              :                WRITE (UNIT=output_unit, FMT=fmtstr2) &
     653          560 :                   iatom, ikind, "  all_potential", qs_force(ikind)%all_potential(1:3, i), &
     654          560 :                   iatom, ikind, "        overlap", qs_force(ikind)%overlap(1:3, i), &
     655          560 :                   iatom, ikind, "       rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
     656          560 :                   iatom, ikind, "      repulsive", qs_force(ikind)%repulsive(1:3, i), &
     657          560 :                   iatom, ikind, "     dispersion", qs_force(ikind)%dispersion(1:3, i), &
     658          560 :                   iatom, ikind, "        efield", qs_force(ikind)%efield(1:3, i), &
     659          560 :                   iatom, ikind, "     ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
     660          700 :                   iatom, ikind, "          total", qs_force(ikind)%total(1:3, i)
     661          580 :                grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
     662              :             END DO
     663              :          CASE (5)
     664          156 :             DO iatom = 1, natom
     665            0 :                ikind = kind_of(iatom)
     666            0 :                i = atom_of_kind(iatom)
     667              :                WRITE (UNIT=output_unit, FMT=fmtstr2) &
     668            0 :                   iatom, ikind, "       overlap", qs_force(ikind)%overlap(1:3, i), &
     669            0 :                   iatom, ikind, "       kinetic", qs_force(ikind)%kinetic(1:3, i), &
     670            0 :                   iatom, ikind, "      rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
     671            0 :                   iatom, ikind, "    dispersion", qs_force(ikind)%dispersion(1:3, i), &
     672            0 :                   iatom, ikind, " all potential", qs_force(ikind)%all_potential(1:3, i), &
     673            0 :                   iatom, ikind, "         other", qs_force(ikind)%other(1:3, i), &
     674            0 :                   iatom, ikind, "         total", qs_force(ikind)%total(1:3, i)
     675            0 :                grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
     676              :             END DO
     677              :          END SELECT
     678              : 
     679          156 :          WRITE (UNIT=output_unit, FMT=fmtstr3) "Sum of total", grand_total(1:3)
     680              : 
     681          156 :          DEALLOCATE (atom_of_kind)
     682          156 :          DEALLOCATE (kind_of)
     683              : 
     684              :       END IF
     685              : 
     686        10581 :    END SUBROUTINE write_forces
     687              : 
     688              : END MODULE qs_force
        

Generated by: LCOV version 2.0-1