Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief 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