LCOV - code coverage report
Current view: top level - src - qmmmx_force.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 90.1 % 71 64
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            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 Calculates QM/MM energy and forces with Force-Mixing
      10              : !> \par History
      11              : !>      2015 Factored out of force_env_methods.F
      12              : !> \author Ole Schuett
      13              : ! **************************************************************************************************
      14              : MODULE qmmmx_force
      15              :    USE cell_types,                      ONLY: cell_type
      16              :    USE cp_subsys_types,                 ONLY: cp_subsys_type
      17              :    USE fist_environment_types,          ONLY: fist_env_get
      18              :    USE input_constants,                 ONLY: do_fm_mom_conserv_QM,&
      19              :                                               do_fm_mom_conserv_buffer,&
      20              :                                               do_fm_mom_conserv_core,&
      21              :                                               do_fm_mom_conserv_equal_a,&
      22              :                                               do_fm_mom_conserv_equal_f,&
      23              :                                               do_fm_mom_conserv_none
      24              :    USE input_section_types,             ONLY: section_vals_type,&
      25              :                                               section_vals_val_get,&
      26              :                                               section_vals_val_set
      27              :    USE kinds,                           ONLY: default_string_length,&
      28              :                                               dp
      29              :    USE particle_types,                  ONLY: particle_type
      30              :    USE qmmm_force,                      ONLY: qmmm_calc_energy_force
      31              :    USE qmmm_types,                      ONLY: qmmm_env_get,&
      32              :                                               qmmm_env_type
      33              :    USE qmmm_types_low,                  ONLY: force_mixing_label_QM_core,&
      34              :                                               force_mixing_label_QM_dynamics,&
      35              :                                               force_mixing_label_buffer
      36              :    USE qmmm_util,                       ONLY: apply_qmmm_unwrap,&
      37              :                                               apply_qmmm_wrap
      38              :    USE qmmmx_types,                     ONLY: qmmmx_env_type
      39              :    USE qmmmx_util,                      ONLY: apply_qmmmx_translate
      40              :    USE qs_environment_types,            ONLY: get_qs_env
      41              : #include "./base/base_uses.f90"
      42              : 
      43              :    IMPLICIT NONE
      44              : 
      45              :    PRIVATE
      46              : 
      47              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmmx_force'
      48              : 
      49              :    PUBLIC :: qmmmx_calc_energy_force
      50              : 
      51              : CONTAINS
      52              : 
      53              : ! **************************************************************************************************
      54              : !> \brief calculates the qm/mm energy and forces
      55              : !> \param qmmmx_env ...
      56              : !> \param calc_force if also the forces should be calculated
      57              : !> \param consistent_energies ...
      58              : !> \param linres ...
      59              : !> \param require_consistent_energy_force ...
      60              : !> \par History
      61              : !>      05.2004 created [fawzi]
      62              : !> \author Fawzi Mohamed
      63              : ! **************************************************************************************************
      64          104 :    SUBROUTINE qmmmx_calc_energy_force(qmmmx_env, calc_force, consistent_energies, linres, &
      65              :                                       require_consistent_energy_force)
      66              :       TYPE(qmmmx_env_type), POINTER                      :: qmmmx_env
      67              :       LOGICAL, INTENT(IN)                                :: calc_force, consistent_energies, linres
      68              :       LOGICAL, INTENT(IN), OPTIONAL :: require_consistent_energy_force
      69              : 
      70              :       INTEGER                                            :: ip, mom_conserv_min_label, &
      71              :                                                             mom_conserv_n, mom_conserv_region, &
      72              :                                                             mom_conserv_type
      73           52 :       INTEGER, POINTER                                   :: cur_indices(:), cur_labels(:)
      74              :       REAL(dp)                                           :: delta_a(3), delta_f(3), &
      75              :                                                             mom_conserv_mass, total_f(3)
      76              :       TYPE(cp_subsys_type), POINTER                      :: subsys_primary, subsys_qmmm_core, &
      77              :                                                             subsys_qmmm_extended
      78           52 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles_primary, particles_qmmm_core, &
      79           52 :                                                             particles_qmmm_extended
      80              :       TYPE(section_vals_type), POINTER                   :: force_env_section
      81              : 
      82           52 :       IF (PRESENT(require_consistent_energy_force)) THEN
      83           48 :          IF (require_consistent_energy_force) &
      84              :             CALL cp_abort(__LOCATION__, &
      85            0 :                           "qmmmx_energy_and_forces got require_consistent_energy_force but force mixing is active. ")
      86              :       END IF
      87              : 
      88              :       ! Possibly translate the system
      89           52 :       CALL apply_qmmmx_translate(qmmmx_env)
      90              : 
      91              :       ! actual energy force calculation
      92           52 :       CALL qmmmx_calc_energy_force_low(qmmmx_env%ext, calc_force, consistent_energies, linres, "ext")
      93           52 :       CALL qmmmx_calc_energy_force_low(qmmmx_env%core, calc_force, consistent_energies, linres, "core")
      94              : 
      95              :       ! get forces from subsys of each sub force env
      96           52 :       CALL qmmm_env_get(qmmmx_env%core, subsys=subsys_qmmm_core)
      97           52 :       CALL qmmm_env_get(qmmmx_env%ext, subsys=subsys_qmmm_extended)
      98              : 
      99           52 :       CALL get_qs_env(qmmmx_env%ext%qs_env, input=force_env_section)
     100           52 :       CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%INDICES", i_vals=cur_indices)
     101           52 :       CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%LABELS", i_vals=cur_labels)
     102              : 
     103           52 :       particles_qmmm_extended => subsys_qmmm_extended%particles%els
     104           52 :       particles_qmmm_core => subsys_qmmm_core%particles%els
     105         1990 :       DO ip = 1, SIZE(cur_indices)
     106         1990 :          IF (cur_labels(ip) >= force_mixing_label_QM_dynamics) THEN ! this is a QM atom
     107              :             ! copy (QM) force from extended calculation
     108         5472 :             particles_qmmm_core(cur_indices(ip))%f = particles_qmmm_extended(cur_indices(ip))%f
     109              :          END IF
     110              :       END DO
     111              : 
     112              :       ! zero momentum
     113              :       CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%MOMENTUM_CONSERVATION_TYPE", &
     114           52 :                                 i_val=mom_conserv_type)
     115           52 :       IF (mom_conserv_type /= do_fm_mom_conserv_none) THEN
     116              :          CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%MOMENTUM_CONSERVATION_REGION", &
     117           52 :                                    i_val=mom_conserv_region)
     118              : 
     119           52 :          IF (mom_conserv_region == do_fm_mom_conserv_core) THEN
     120              :             mom_conserv_min_label = force_mixing_label_QM_core
     121              :          ELSEIF (mom_conserv_region == do_fm_mom_conserv_QM) THEN
     122              :             mom_conserv_min_label = force_mixing_label_QM_dynamics
     123              :          ELSEIF (mom_conserv_region == do_fm_mom_conserv_buffer) THEN
     124              :             mom_conserv_min_label = force_mixing_label_buffer
     125              :          ELSE
     126            0 :             CPABORT("Got unknown MOMENTUM_CONSERVATION_REGION (not CORE, QM, or BUFFER) !")
     127              :          END IF
     128              : 
     129           52 :          total_f = 0.0_dp
     130        99838 :          DO ip = 1, SIZE(particles_qmmm_core)
     131       399196 :             total_f(1:3) = total_f(1:3) + particles_qmmm_core(ip)%f(1:3)
     132              :          END DO
     133           52 :          IF (mom_conserv_type == do_fm_mom_conserv_equal_f) THEN
     134            0 :             mom_conserv_n = COUNT(cur_labels >= mom_conserv_min_label)
     135            0 :             delta_f = total_f/mom_conserv_n
     136            0 :             DO ip = 1, SIZE(cur_indices)
     137            0 :                IF (cur_labels(ip) >= mom_conserv_min_label) THEN
     138            0 :                   particles_qmmm_core(cur_indices(ip))%f = particles_qmmm_core(cur_indices(ip))%f - delta_f
     139              :                END IF
     140              :             END DO
     141           52 :          ELSE IF (mom_conserv_type == do_fm_mom_conserv_equal_a) THEN
     142           52 :             mom_conserv_mass = 0.0_dp
     143         1990 :             DO ip = 1, SIZE(cur_indices)
     144         1938 :                IF (cur_labels(ip) >= mom_conserv_min_label) &
     145          736 :                   mom_conserv_mass = mom_conserv_mass + particles_qmmm_core(cur_indices(ip))%atomic_kind%mass
     146              :             END DO
     147          208 :             delta_a = total_f/mom_conserv_mass
     148         1990 :             DO ip = 1, SIZE(cur_indices)
     149         1990 :                IF (cur_labels(ip) >= mom_conserv_min_label) THEN
     150              :                   particles_qmmm_core(cur_indices(ip))%f = particles_qmmm_core(cur_indices(ip))%f - &
     151         2736 :                                                            particles_qmmm_core(cur_indices(ip))%atomic_kind%mass*delta_a
     152              :                END IF
     153              :             END DO
     154              :          END IF
     155              :       END IF
     156              : 
     157           52 :       CALL qmmm_env_get(qmmmx_env%ext, subsys=subsys_primary)
     158           52 :       particles_primary => subsys_primary%particles%els
     159        99838 :       DO ip = 1, SIZE(particles_qmmm_core)
     160       798340 :          particles_primary(ip)%f = particles_qmmm_core(ip)%f
     161              :       END DO
     162              : 
     163           52 :    END SUBROUTINE qmmmx_calc_energy_force
     164              : 
     165              : ! **************************************************************************************************
     166              : !> \brief ...
     167              : !> \param qmmm_env ...
     168              : !> \param calc_force ...
     169              : !> \param consistent_energies ...
     170              : !> \param linres ...
     171              : !> \param label ...
     172              : ! **************************************************************************************************
     173          104 :    SUBROUTINE qmmmx_calc_energy_force_low(qmmm_env, calc_force, consistent_energies, linres, label)
     174              :       TYPE(qmmm_env_type), POINTER                       :: qmmm_env
     175              :       LOGICAL, INTENT(IN)                                :: calc_force, consistent_energies, linres
     176              :       CHARACTER(*)                                       :: label
     177              : 
     178              :       CHARACTER(default_string_length)                   :: new_restart_fn, new_restart_hist_fn, &
     179              :                                                             old_restart_fn, old_restart_hist_fn
     180          104 :       INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index
     181              :       LOGICAL                                            :: saved_do_translate
     182          104 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: saved_pos
     183              :       TYPE(cell_type), POINTER                           :: mm_cell
     184              :       TYPE(cp_subsys_type), POINTER                      :: subsys_mm, subsys_qm
     185              :       TYPE(section_vals_type), POINTER                   :: force_env_section
     186              : 
     187          104 :       NULLIFY (mm_cell, subsys_qm, subsys_mm, qm_atom_index)
     188              : 
     189          104 :       CALL get_qs_env(qmmm_env%qs_env, input=force_env_section)
     190              : 
     191              :       ! rewrite RESTART%FILENAME
     192              :       CALL section_vals_val_get(force_env_section, "DFT%SCF%PRINT%RESTART%FILENAME", &
     193          104 :                                 c_val=old_restart_fn)
     194          104 :       new_restart_fn = TRIM(old_restart_fn)//"-"//TRIM(label)
     195              :       CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART%FILENAME", &
     196          104 :                                 c_val=new_restart_fn)
     197              : 
     198              :       ! rewrite RESTART_HISTORY%FILENAME
     199              :       CALL section_vals_val_get(force_env_section, "DFT%SCF%PRINT%RESTART_HISTORY%FILENAME", &
     200          104 :                                 c_val=old_restart_hist_fn)
     201          104 :       new_restart_hist_fn = TRIM(old_restart_hist_fn)//"-"//TRIM(label)
     202              :       CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART_HISTORY%FILENAME", &
     203          104 :                                 c_val=new_restart_hist_fn)
     204              : 
     205              :       ! wrap positions before QM/MM calculation.
     206              :       ! Required if diffusion causes atoms outside of periodic box get added to QM
     207          104 :       CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
     208          104 :       CALL get_qs_env(qmmm_env%qs_env, cp_subsys=subsys_qm)
     209          104 :       qm_atom_index => qmmm_env%qm%qm_atom_index
     210          104 :       CALL apply_qmmm_wrap(subsys_mm, mm_cell, subsys_qm, qm_atom_index, saved_pos)
     211              : 
     212              :       ! Turn off box translation, it was already performed by apply_qmmmx_translate(),
     213              :       ! the particles coordinates will still be copied from MM to QM.
     214          104 :       saved_do_translate = qmmm_env%qm%do_translate
     215          104 :       qmmm_env%qm%do_translate = .FALSE.
     216              : 
     217              :       ! actual energy force calculation
     218          104 :       CALL qmmm_calc_energy_force(qmmm_env, calc_force, consistent_energies, linres)
     219              : 
     220              :       ! restore do_translate
     221          104 :       qmmm_env%qm%do_translate = saved_do_translate
     222              : 
     223              :       ! restore unwrapped positions
     224          104 :       CALL apply_qmmm_unwrap(subsys_mm, subsys_qm, qm_atom_index, saved_pos)
     225              : 
     226              :       ! restore RESTART filenames
     227              :       CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART%FILENAME", &
     228          104 :                                 c_val=old_restart_fn)
     229              :       CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART_HISTORY%FILENAME", &
     230          104 :                                 c_val=old_restart_hist_fn)
     231              : 
     232          208 :    END SUBROUTINE qmmmx_calc_energy_force_low
     233              : 
     234              : END MODULE qmmmx_force
        

Generated by: LCOV version 2.0-1