Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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) THEN
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 : END IF
88 :
89 : ! Possibly translate the system
90 52 : CALL apply_qmmmx_translate(qmmmx_env)
91 :
92 : ! actual energy force calculation
93 52 : CALL qmmmx_calc_energy_force_low(qmmmx_env%ext, calc_force, consistent_energies, linres, "ext")
94 52 : CALL qmmmx_calc_energy_force_low(qmmmx_env%core, calc_force, consistent_energies, linres, "core")
95 :
96 : ! get forces from subsys of each sub force env
97 52 : CALL qmmm_env_get(qmmmx_env%core, subsys=subsys_qmmm_core)
98 52 : CALL qmmm_env_get(qmmmx_env%ext, subsys=subsys_qmmm_extended)
99 :
100 52 : CALL get_qs_env(qmmmx_env%ext%qs_env, input=force_env_section)
101 52 : CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%INDICES", i_vals=cur_indices)
102 52 : CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%LABELS", i_vals=cur_labels)
103 :
104 52 : particles_qmmm_extended => subsys_qmmm_extended%particles%els
105 52 : particles_qmmm_core => subsys_qmmm_core%particles%els
106 1990 : DO ip = 1, SIZE(cur_indices)
107 1990 : IF (cur_labels(ip) >= force_mixing_label_QM_dynamics) THEN ! this is a QM atom
108 : ! copy (QM) force from extended calculation
109 5472 : particles_qmmm_core(cur_indices(ip))%f = particles_qmmm_extended(cur_indices(ip))%f
110 : END IF
111 : END DO
112 :
113 : ! zero momentum
114 : CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%MOMENTUM_CONSERVATION_TYPE", &
115 52 : i_val=mom_conserv_type)
116 52 : IF (mom_conserv_type /= do_fm_mom_conserv_none) THEN
117 : CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%MOMENTUM_CONSERVATION_REGION", &
118 52 : i_val=mom_conserv_region)
119 :
120 52 : IF (mom_conserv_region == do_fm_mom_conserv_core) THEN
121 : mom_conserv_min_label = force_mixing_label_QM_core
122 : ELSE IF (mom_conserv_region == do_fm_mom_conserv_QM) THEN
123 : mom_conserv_min_label = force_mixing_label_QM_dynamics
124 : ELSE IF (mom_conserv_region == do_fm_mom_conserv_buffer) THEN
125 : mom_conserv_min_label = force_mixing_label_buffer
126 : ELSE
127 0 : CPABORT("Got unknown MOMENTUM_CONSERVATION_REGION (not CORE, QM, or BUFFER) !")
128 : END IF
129 :
130 52 : total_f = 0.0_dp
131 99838 : DO ip = 1, SIZE(particles_qmmm_core)
132 399196 : total_f(1:3) = total_f(1:3) + particles_qmmm_core(ip)%f(1:3)
133 : END DO
134 52 : IF (mom_conserv_type == do_fm_mom_conserv_equal_f) THEN
135 0 : mom_conserv_n = COUNT(cur_labels >= mom_conserv_min_label)
136 0 : delta_f = total_f/mom_conserv_n
137 0 : DO ip = 1, SIZE(cur_indices)
138 0 : IF (cur_labels(ip) >= mom_conserv_min_label) THEN
139 0 : particles_qmmm_core(cur_indices(ip))%f = particles_qmmm_core(cur_indices(ip))%f - delta_f
140 : END IF
141 : END DO
142 52 : ELSE IF (mom_conserv_type == do_fm_mom_conserv_equal_a) THEN
143 52 : mom_conserv_mass = 0.0_dp
144 1990 : DO ip = 1, SIZE(cur_indices)
145 1990 : IF (cur_labels(ip) >= mom_conserv_min_label) THEN
146 684 : mom_conserv_mass = mom_conserv_mass + particles_qmmm_core(cur_indices(ip))%atomic_kind%mass
147 : END IF
148 : END DO
149 208 : delta_a = total_f/mom_conserv_mass
150 1990 : DO ip = 1, SIZE(cur_indices)
151 1990 : IF (cur_labels(ip) >= mom_conserv_min_label) THEN
152 : particles_qmmm_core(cur_indices(ip))%f = particles_qmmm_core(cur_indices(ip))%f - &
153 2736 : particles_qmmm_core(cur_indices(ip))%atomic_kind%mass*delta_a
154 : END IF
155 : END DO
156 : END IF
157 : END IF
158 :
159 52 : CALL qmmm_env_get(qmmmx_env%ext, subsys=subsys_primary)
160 52 : particles_primary => subsys_primary%particles%els
161 99838 : DO ip = 1, SIZE(particles_qmmm_core)
162 798340 : particles_primary(ip)%f = particles_qmmm_core(ip)%f
163 : END DO
164 :
165 52 : END SUBROUTINE qmmmx_calc_energy_force
166 :
167 : ! **************************************************************************************************
168 : !> \brief ...
169 : !> \param qmmm_env ...
170 : !> \param calc_force ...
171 : !> \param consistent_energies ...
172 : !> \param linres ...
173 : !> \param label ...
174 : ! **************************************************************************************************
175 104 : SUBROUTINE qmmmx_calc_energy_force_low(qmmm_env, calc_force, consistent_energies, linres, label)
176 : TYPE(qmmm_env_type), POINTER :: qmmm_env
177 : LOGICAL, INTENT(IN) :: calc_force, consistent_energies, linres
178 : CHARACTER(*) :: label
179 :
180 : CHARACTER(default_string_length) :: new_restart_fn, new_restart_hist_fn, &
181 : old_restart_fn, old_restart_hist_fn
182 104 : INTEGER, DIMENSION(:), POINTER :: qm_atom_index
183 : LOGICAL :: saved_do_translate
184 104 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: saved_pos
185 : TYPE(cell_type), POINTER :: mm_cell
186 : TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm
187 : TYPE(section_vals_type), POINTER :: force_env_section
188 :
189 104 : NULLIFY (mm_cell, subsys_qm, subsys_mm, qm_atom_index)
190 :
191 104 : CALL get_qs_env(qmmm_env%qs_env, input=force_env_section)
192 :
193 : ! rewrite RESTART%FILENAME
194 : CALL section_vals_val_get(force_env_section, "DFT%SCF%PRINT%RESTART%FILENAME", &
195 104 : c_val=old_restart_fn)
196 104 : new_restart_fn = TRIM(old_restart_fn)//"-"//TRIM(label)
197 : CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART%FILENAME", &
198 104 : c_val=new_restart_fn)
199 :
200 : ! rewrite RESTART_HISTORY%FILENAME
201 : CALL section_vals_val_get(force_env_section, "DFT%SCF%PRINT%RESTART_HISTORY%FILENAME", &
202 104 : c_val=old_restart_hist_fn)
203 104 : new_restart_hist_fn = TRIM(old_restart_hist_fn)//"-"//TRIM(label)
204 : CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART_HISTORY%FILENAME", &
205 104 : c_val=new_restart_hist_fn)
206 :
207 : ! wrap positions before QM/MM calculation.
208 : ! Required if diffusion causes atoms outside of periodic box get added to QM
209 104 : CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
210 104 : CALL get_qs_env(qmmm_env%qs_env, cp_subsys=subsys_qm)
211 104 : qm_atom_index => qmmm_env%qm%qm_atom_index
212 104 : CALL apply_qmmm_wrap(subsys_mm, mm_cell, subsys_qm, qm_atom_index, saved_pos)
213 :
214 : ! Turn off box translation, it was already performed by apply_qmmmx_translate(),
215 : ! the particles coordinates will still be copied from MM to QM.
216 104 : saved_do_translate = qmmm_env%qm%do_translate
217 104 : qmmm_env%qm%do_translate = .FALSE.
218 :
219 : ! actual energy force calculation
220 104 : CALL qmmm_calc_energy_force(qmmm_env, calc_force, consistent_energies, linres)
221 :
222 : ! restore do_translate
223 104 : qmmm_env%qm%do_translate = saved_do_translate
224 :
225 : ! restore unwrapped positions
226 104 : CALL apply_qmmm_unwrap(subsys_mm, subsys_qm, qm_atom_index, saved_pos)
227 :
228 : ! restore RESTART filenames
229 : CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART%FILENAME", &
230 104 : c_val=old_restart_fn)
231 : CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART_HISTORY%FILENAME", &
232 104 : c_val=old_restart_hist_fn)
233 :
234 208 : END SUBROUTINE qmmmx_calc_energy_force_low
235 :
236 : END MODULE qmmmx_force
|
