(git:ccc2433)
qmmmx_force.F
Go to the documentation of this file.
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  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  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  TYPE(particle_type), DIMENSION(:), POINTER :: particles_primary, particles_qmmm_core, &
79  particles_qmmm_extended
80  TYPE(section_vals_type), POINTER :: force_env_section
81 
82  IF (PRESENT(require_consistent_energy_force)) THEN
83  IF (require_consistent_energy_force) &
84  CALL cp_abort(__location__, &
85  "qmmmx_energy_and_forces got require_consistent_energy_force but force mixing is active. ")
86  END IF
87 
88  ! Possibly translate the system
89  CALL apply_qmmmx_translate(qmmmx_env)
90 
91  ! actual energy force calculation
92  CALL qmmmx_calc_energy_force_low(qmmmx_env%ext, calc_force, consistent_energies, linres, "ext")
93  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  CALL qmmm_env_get(qmmmx_env%core, subsys=subsys_qmmm_core)
97  CALL qmmm_env_get(qmmmx_env%ext, subsys=subsys_qmmm_extended)
98 
99  CALL get_qs_env(qmmmx_env%ext%qs_env, input=force_env_section)
100  CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%INDICES", i_vals=cur_indices)
101  CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%LABELS", i_vals=cur_labels)
102 
103  particles_qmmm_extended => subsys_qmmm_extended%particles%els
104  particles_qmmm_core => subsys_qmmm_core%particles%els
105  DO ip = 1, SIZE(cur_indices)
106  IF (cur_labels(ip) >= force_mixing_label_qm_dynamics) THEN ! this is a QM atom
107  ! copy (QM) force from extended calculation
108  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  i_val=mom_conserv_type)
115  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  i_val=mom_conserv_region)
118 
119  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  cpabort("Got unknown MOMENTUM_CONSERVATION_REGION (not CORE, QM, or BUFFER) !")
127  END IF
128 
129  total_f = 0.0_dp
130  DO ip = 1, SIZE(particles_qmmm_core)
131  total_f(1:3) = total_f(1:3) + particles_qmmm_core(ip)%f(1:3)
132  END DO
133  IF (mom_conserv_type == do_fm_mom_conserv_equal_f) THEN
134  mom_conserv_n = count(cur_labels >= mom_conserv_min_label)
135  delta_f = total_f/mom_conserv_n
136  DO ip = 1, SIZE(cur_indices)
137  IF (cur_labels(ip) >= mom_conserv_min_label) THEN
138  particles_qmmm_core(cur_indices(ip))%f = particles_qmmm_core(cur_indices(ip))%f - delta_f
139  END IF
140  END DO
141  ELSE IF (mom_conserv_type == do_fm_mom_conserv_equal_a) THEN
142  mom_conserv_mass = 0.0_dp
143  DO ip = 1, SIZE(cur_indices)
144  IF (cur_labels(ip) >= mom_conserv_min_label) &
145  mom_conserv_mass = mom_conserv_mass + particles_qmmm_core(cur_indices(ip))%atomic_kind%mass
146  END DO
147  delta_a = total_f/mom_conserv_mass
148  DO ip = 1, SIZE(cur_indices)
149  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  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  CALL qmmm_env_get(qmmmx_env%ext, subsys=subsys_primary)
158  particles_primary => subsys_primary%particles%els
159  DO ip = 1, SIZE(particles_qmmm_core)
160  particles_primary(ip)%f = particles_qmmm_core(ip)%f
161  END DO
162 
163  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  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  INTEGER, DIMENSION(:), POINTER :: qm_atom_index
181  LOGICAL :: saved_do_translate
182  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  NULLIFY (mm_cell, subsys_qm, subsys_mm, qm_atom_index)
188 
189  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  c_val=old_restart_fn)
194  new_restart_fn = trim(old_restart_fn)//"-"//trim(label)
195  CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART%FILENAME", &
196  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  c_val=old_restart_hist_fn)
201  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  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  CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
208  CALL get_qs_env(qmmm_env%qs_env, cp_subsys=subsys_qm)
209  qm_atom_index => qmmm_env%qm%qm_atom_index
210  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  saved_do_translate = qmmm_env%qm%do_translate
215  qmmm_env%qm%do_translate = .false.
216 
217  ! actual energy force calculation
218  CALL qmmm_calc_energy_force(qmmm_env, calc_force, consistent_energies, linres)
219 
220  ! restore do_translate
221  qmmm_env%qm%do_translate = saved_do_translate
222 
223  ! restore unwrapped positions
224  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  c_val=old_restart_fn)
229  CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART_HISTORY%FILENAME", &
230  c_val=old_restart_hist_fn)
231 
232  END SUBROUTINE qmmmx_calc_energy_force_low
233 
234 END MODULE qmmmx_force
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition: cp_subsys_types.F:16
fist_environment_types
Definition: fist_environment_types.F:15
fist_environment_types::fist_env_get
subroutine, public fist_env_get(fist_env, atomic_kind_set, particle_set, ewald_pw, local_particles, local_molecules, molecule_kind_set, molecule_set, cell, cell_ref, ewald_env, fist_nonbond_env, thermo, para_env, subsys, qmmm, qmmm_env, input, shell_model, shell_model_ad, shell_particle_set, core_particle_set, multipoles, results, exclusions, efield)
Purpose: Get the FIST environment.
Definition: fist_environment_types.F:140
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_fm_mom_conserv_none
integer, parameter, public do_fm_mom_conserv_none
Definition: input_constants.F:404
input_constants::do_fm_mom_conserv_buffer
integer, parameter, public do_fm_mom_conserv_buffer
Definition: input_constants.F:408
input_constants::do_fm_mom_conserv_equal_f
integer, parameter, public do_fm_mom_conserv_equal_f
Definition: input_constants.F:404
input_constants::do_fm_mom_conserv_qm
integer, parameter, public do_fm_mom_conserv_qm
Definition: input_constants.F:408
input_constants::do_fm_mom_conserv_equal_a
integer, parameter, public do_fm_mom_conserv_equal_a
Definition: input_constants.F:404
input_constants::do_fm_mom_conserv_core
integer, parameter, public do_fm_mom_conserv_core
Definition: input_constants.F:408
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition: input_section_types.F:1216
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qmmm_force
Calculates QM/MM energy and forces.
Definition: qmmm_force.F:14
qmmm_force::qmmm_calc_energy_force
subroutine, public qmmm_calc_energy_force(qmmm_env, calc_force, consistent_energies, linres)
calculates the qm/mm energy and forces
Definition: qmmm_force.F:76
qmmm_types_low
Definition: qmmm_types_low.F:13
qmmm_types_low::force_mixing_label_qm_core
integer, parameter, public force_mixing_label_qm_core
Definition: qmmm_types_low.F:44
qmmm_types_low::force_mixing_label_buffer
integer, parameter, public force_mixing_label_buffer
Definition: qmmm_types_low.F:44
qmmm_types_low::force_mixing_label_qm_dynamics
integer, parameter, public force_mixing_label_qm_dynamics
Definition: qmmm_types_low.F:44
qmmm_types
Basic container type for QM/MM.
Definition: qmmm_types.F:12
qmmm_types::qmmm_env_get
subroutine, public qmmm_env_get(qmmm_env, subsys, potential_energy, kinetic_energy)
...
Definition: qmmm_types.F:50
qmmm_util
Definition: qmmm_util.F:13
qmmm_util::apply_qmmm_unwrap
subroutine, public apply_qmmm_unwrap(subsys_mm, subsys_qm, qm_atom_index, saved_pos)
...
Definition: qmmm_util.F:344
qmmm_util::apply_qmmm_wrap
subroutine, public apply_qmmm_wrap(subsys_mm, mm_cell, subsys_qm, qm_atom_index, saved_pos)
wrap positions (with mm periodicity)
Definition: qmmm_util.F:308
qmmmx_force
Calculates QM/MM energy and forces with Force-Mixing.
Definition: qmmmx_force.F:14
qmmmx_force::qmmmx_calc_energy_force
subroutine, public qmmmx_calc_energy_force(qmmmx_env, calc_force, consistent_energies, linres, require_consistent_energy_force)
calculates the qm/mm energy and forces
Definition: qmmmx_force.F:66
qmmmx_types
Basic container type for QM/MM with force mixing.
Definition: qmmmx_types.F:12
qmmmx_util
Routines used for force-mixing QM/MM calculations.
Definition: qmmmx_util.F:14
qmmmx_util::apply_qmmmx_translate
subroutine, public apply_qmmmx_translate(qmmmx_env)
Apply translation to the full system in order to center the QM system into the QM box.
Definition: qmmmx_util.F:75
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502