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 Routines used for force-mixing QM/MM calculations
10 : !> \par History
11 : !> 2.2012 created [noam]
12 : !> \author Noam Bernstein
13 : ! **************************************************************************************************
14 : MODULE qmmmx_util
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind,&
17 : set_atomic_kind
18 : USE cell_types, ONLY: cell_copy,&
19 : cell_type,&
20 : pbc
21 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
22 : cp_logger_get_default_unit_nr
23 : USE cp_subsys_types, ONLY: cp_subsys_get,&
24 : cp_subsys_type
25 : USE fist_neighbor_list_types, ONLY: fist_neighbor_deallocate,&
26 : fist_neighbor_type
27 : USE fist_neighbor_lists, ONLY: build_fist_neighbor_lists
28 : USE input_section_types, ONLY: &
29 : section_vals_add_values, section_vals_duplicate, section_vals_get, &
30 : section_vals_get_subs_vals, section_vals_get_subs_vals3, section_vals_release, &
31 : section_vals_remove_values, section_vals_set_subs_vals, section_vals_type, &
32 : section_vals_val_get, section_vals_val_set, section_vals_write
33 : USE kinds, ONLY: default_string_length,&
34 : dp
35 : USE memory_utilities, ONLY: reallocate
36 : USE message_passing, ONLY: mp_para_env_type
37 : USE molecule_list_types, ONLY: molecule_list_type
38 : USE molecule_types, ONLY: molecule_type
39 : USE particle_list_types, ONLY: particle_list_type
40 : USE particle_types, ONLY: particle_type
41 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
42 : USE qmmm_types, ONLY: qmmm_env_get
43 : USE qmmm_types_low, ONLY: force_mixing_label_QM_core,&
44 : force_mixing_label_QM_core_list,&
45 : force_mixing_label_QM_dynamics,&
46 : force_mixing_label_QM_dynamics_list,&
47 : force_mixing_label_buffer,&
48 : force_mixing_label_buffer_list,&
49 : force_mixing_label_none,&
50 : force_mixing_label_termination
51 : USE qmmm_util, ONLY: apply_qmmm_translate
52 : USE qmmmx_types, ONLY: qmmmx_env_type
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 : PRIVATE
57 :
58 : LOGICAL, PRIVATE :: debug_this_module = .FALSE.
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmmx_util'
60 :
61 : PUBLIC :: setup_force_mixing_qmmm_sections, update_force_mixing_labels, &
62 : apply_qmmmx_translate
63 :
64 : CONTAINS
65 :
66 : ! **************************************************************************************************
67 : !> \brief Apply translation to the full system in order to center the QM
68 : !> system into the QM box
69 : !> \param qmmmx_env ...
70 : !> \par History
71 : !> 08.2007 created [tlaino] - Zurich University
72 : !> \author Teodoro Laino
73 : ! **************************************************************************************************
74 52 : SUBROUTINE apply_qmmmx_translate(qmmmx_env)
75 : TYPE(qmmmx_env_type), POINTER :: qmmmx_env
76 :
77 : INTEGER :: ip
78 : TYPE(cell_type), POINTER :: cell_core, cell_extended
79 : TYPE(cp_subsys_type), POINTER :: subsys_core, subsys_extended
80 52 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_core, particles_extended
81 :
82 52 : NULLIFY (cell_core, cell_extended)
83 52 : NULLIFY (subsys_core, subsys_extended)
84 52 : NULLIFY (particles_core, particles_extended)
85 :
86 : ! want to center extended, and make core consistent with that
87 52 : CALL apply_qmmm_translate(qmmmx_env%ext)
88 :
89 : ! translate core fist particles
90 52 : CALL qmmm_env_get(qmmmx_env%ext, subsys=subsys_extended)
91 52 : CALL cp_subsys_get(subsys_extended, cell=cell_extended)
92 52 : CALL qmmm_env_get(qmmmx_env%core, subsys=subsys_core)
93 52 : CALL cp_subsys_get(subsys_core, cell=cell_core)
94 52 : particles_extended => subsys_extended%particles%els
95 52 : particles_core => subsys_core%particles%els
96 99838 : DO ip = 1, SIZE(particles_extended)
97 798340 : particles_core(ip)%r = particles_extended(ip)%r
98 : END DO
99 52 : CALL cell_copy(cell_extended, cell_core)
100 :
101 : ! The core QM particles will be updated the regular call
102 : ! to apply_qmmm_translate() from within qmmm_calc_energy_force()
103 :
104 52 : END SUBROUTINE apply_qmmmx_translate
105 :
106 : ! **************************************************************************************************
107 : !> \brief ...
108 : !> \param subsys ...
109 : !> \param qmmm_section ...
110 : !> \param labels_changed ...
111 : !> \par History
112 : !> 02.2012 created [noam]
113 : !> \author Noam Bernstein
114 : ! **************************************************************************************************
115 56 : SUBROUTINE update_force_mixing_labels(subsys, qmmm_section, labels_changed)
116 : TYPE(cp_subsys_type), POINTER :: subsys
117 : TYPE(section_vals_type), POINTER :: qmmm_section
118 : LOGICAL, OPTIONAL :: labels_changed
119 :
120 56 : CHARACTER(LEN=default_string_length), POINTER :: adaptive_exclude_molecules(:)
121 : INTEGER :: i_rep_section, i_rep_val, ip, max_n_qm, n_new, n_rep_exclude, n_rep_section, &
122 : n_rep_val, natoms, output_unit, QM_extended_seed_min_label_val
123 56 : INTEGER, ALLOCATABLE :: new_full_labels(:), orig_full_labels(:)
124 56 : INTEGER, POINTER :: broken_bonds(:), cur_indices(:), &
125 56 : cur_labels(:), mm_index_entry(:), &
126 56 : new_indices(:), new_labels(:)
127 : LOGICAL :: explicit, QM_extended_seed_is_core_list
128 56 : REAL(dp), ALLOCATABLE :: nearest_dist(:)
129 56 : REAL(dp), POINTER :: r_buf(:), r_core(:), r_qm(:)
130 : TYPE(cell_type), POINTER :: cell
131 : TYPE(fist_neighbor_type), POINTER :: nlist
132 : TYPE(molecule_list_type), POINTER :: molecules
133 56 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
134 : TYPE(mp_para_env_type), POINTER :: para_env
135 : TYPE(particle_list_type), POINTER :: particles
136 56 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
137 : TYPE(section_vals_type), POINTER :: force_mixing_section, &
138 : non_adaptive_section, qm_kind_section, &
139 : restart_section
140 :
141 112 : output_unit = cp_logger_get_default_io_unit()
142 :
143 56 : IF (debug_this_module .AND. output_unit > 0) WRITE (output_unit, *) "BOB starting update_force_mixing_labels"
144 : ! get cur indices, labels
145 56 : force_mixing_section => section_vals_get_subs_vals3(qmmm_section, "FORCE_MIXING")
146 56 : CALL get_force_mixing_indices(force_mixing_section, cur_indices, cur_labels)
147 56 : IF (debug_this_module .AND. output_unit > 0) WRITE (output_unit, *) "BOB got cur_indices ", SIZE(cur_indices)
148 56 : IF (debug_this_module .AND. output_unit > 0) WRITE (output_unit, *) "BOB got cur_labels ", SIZE(cur_labels)
149 :
150 : ! read from input
151 : ![NB] breakable bonds will come from here, too
152 56 : NULLIFY (r_core, r_qm, r_buf, adaptive_exclude_molecules, broken_bonds)
153 56 : CALL section_vals_val_get(force_mixing_section, "R_CORE", r_vals=r_core)
154 56 : CALL section_vals_val_get(force_mixing_section, "R_QM", r_vals=r_qm)
155 : CALL section_vals_val_get(force_mixing_section, "QM_EXTENDED_SEED_IS_ONLY_CORE_LIST", &
156 56 : l_val=QM_extended_seed_is_core_list)
157 56 : CALL section_vals_val_get(force_mixing_section, "R_BUF", r_vals=r_buf)
158 56 : CALL section_vals_val_get(force_mixing_section, "MAX_N_QM", i_val=max_n_qm)
159 :
160 56 : CALL section_vals_val_get(force_mixing_section, "ADAPTIVE_EXCLUDE_MOLECULES", n_rep_val=n_rep_exclude)
161 56 : IF (n_rep_exclude > 0) THEN
162 12 : CALL section_vals_val_get(force_mixing_section, "ADAPTIVE_EXCLUDE_MOLECULES", c_vals=adaptive_exclude_molecules)
163 : END IF
164 : ![NB] need to read real list from input
165 : ! should be 2xN_bb integer arrays, with (1,:) indices of inside atoms, and (2,:) indices of outside atoms
166 : ! maybe also breakable_bond_types, with _atomic numbers_ of inside/outside atoms?
167 : ! separate lists for core/buffer?
168 :
169 : ! get particles, molecules
170 56 : NULLIFY (particles, molecules)
171 56 : CALL cp_subsys_get(subsys=subsys, particles=particles, molecules=molecules)
172 56 : particle_set => particles%els
173 56 : molecule_set => molecules%els
174 :
175 56 : natoms = SIZE(particle_set)
176 :
177 : ! initialize new indices, labels, and new_full_labels
178 56 : NULLIFY (new_indices, new_labels)
179 56 : CALL reallocate(new_indices, 1, SIZE(cur_indices))
180 56 : CALL reallocate(new_labels, 1, SIZE(cur_labels))
181 1484 : new_indices = 0
182 1484 : new_labels = force_mixing_label_none
183 168 : ALLOCATE (new_full_labels(natoms))
184 99902 : new_full_labels = force_mixing_label_none
185 :
186 : ! neighbor list for various hysteretic distance calls
187 56 : NULLIFY (cell)
188 56 : CALL cp_subsys_get(subsys, cell=cell)
189 56 : NULLIFY (nlist)
190 56 : CALL make_neighbor_list(force_mixing_section, subsys, cell, MAX(r_core(2), r_qm(2), r_buf(2)), nlist)
191 :
192 : ! create labels for core_list from QM_KIND
193 56 : NULLIFY (mm_index_entry)
194 56 : qm_kind_section => section_vals_get_subs_vals3(qmmm_section, "QM_KIND")
195 56 : CALL section_vals_get(qm_kind_section, n_repetition=n_rep_section)
196 56 : n_new = 0
197 174 : DO i_rep_section = 1, n_rep_section
198 118 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, n_rep_val=n_rep_val)
199 304 : DO i_rep_val = 1, n_rep_val
200 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, i_rep_val=i_rep_val, &
201 130 : i_vals=mm_index_entry)
202 422 : DO ip = 1, SIZE(mm_index_entry)
203 : CALL add_new_label(mm_index_entry(ip), force_mixing_label_QM_core_list, n_new, new_indices, new_labels, &
204 304 : new_full_labels, max_n_qm)
205 : END DO ! ip
206 : END DO ! i_rep_val
207 : END DO ! i_rep_section
208 :
209 56 : IF (debug_this_module .AND. output_unit > 0) THEN
210 0 : WRITE (output_unit, *) "BOB core_list new_indices ", new_indices(1:n_new)
211 0 : WRITE (output_unit, *) "BOB core_list new_labels ", new_labels(1:n_new)
212 : END IF
213 :
214 : ! create labels for non adaptive QM and buffer regions from *_NON_ADAPTIVE&QM_KIND sections
215 : non_adaptive_section => section_vals_get_subs_vals(qmmm_section, "FORCE_MIXING%QM_NON_ADAPTIVE", &
216 56 : can_return_null=.TRUE.)
217 56 : IF (ASSOCIATED(non_adaptive_section)) THEN
218 56 : qm_kind_section => section_vals_get_subs_vals3(non_adaptive_section, "QM_KIND")
219 56 : CALL section_vals_get(qm_kind_section, n_repetition=n_rep_section)
220 86 : DO i_rep_section = 1, n_rep_section
221 30 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, n_rep_val=n_rep_val)
222 164 : DO i_rep_val = 1, n_rep_val
223 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, i_rep_val=i_rep_val, &
224 78 : i_vals=mm_index_entry)
225 186 : DO ip = 1, SIZE(mm_index_entry)
226 : CALL add_new_label(mm_index_entry(ip), force_mixing_label_QM_dynamics_list, n_new, new_indices, new_labels, &
227 156 : new_full_labels, max_n_qm)
228 : END DO ! ip
229 : END DO ! i_rep_val
230 : END DO ! i_rep_section
231 : END IF
232 56 : IF (debug_this_module .AND. output_unit > 0) THEN
233 0 : WRITE (output_unit, *) "BOB core_list + non adaptive QM new_indices ", new_indices(1:n_new)
234 0 : WRITE (output_unit, *) "BOB core_list + non adaptive QM new_labels ", new_labels(1:n_new)
235 : END IF
236 : non_adaptive_section => section_vals_get_subs_vals(qmmm_section, "FORCE_MIXING%BUFFER_NON_ADAPTIVE", &
237 56 : can_return_null=.TRUE.)
238 56 : IF (ASSOCIATED(non_adaptive_section)) THEN
239 56 : qm_kind_section => section_vals_get_subs_vals3(non_adaptive_section, "QM_KIND")
240 56 : CALL section_vals_get(qm_kind_section, n_repetition=n_rep_section)
241 92 : DO i_rep_section = 1, n_rep_section
242 36 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, n_rep_val=n_rep_val)
243 692 : DO i_rep_val = 1, n_rep_val
244 : CALL section_vals_val_get(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section, i_rep_val=i_rep_val, &
245 600 : i_vals=mm_index_entry)
246 1236 : DO ip = 1, SIZE(mm_index_entry)
247 : CALL add_new_label(mm_index_entry(ip), force_mixing_label_buffer_list, n_new, new_indices, new_labels, &
248 1200 : new_full_labels, max_n_qm)
249 : END DO ! ip
250 : END DO ! i_rep_val
251 : END DO ! i_rep_section
252 : END IF
253 :
254 56 : IF (debug_this_module .AND. output_unit > 0) THEN
255 0 : WRITE (output_unit, *) "BOB core_list + non adaptive QM+buffer new_indices ", new_indices(1:n_new)
256 0 : WRITE (output_unit, *) "BOB core_list + non adaptive QM+buffer new_labels ", new_labels(1:n_new)
257 : END IF
258 :
259 : ! allocate and initialize full atom set labels for hysteretic loops
260 168 : ALLOCATE (nearest_dist(natoms))
261 :
262 : ! orig_full_labels is full array (natoms) with orig labels
263 112 : ALLOCATE (orig_full_labels(natoms))
264 99902 : orig_full_labels = force_mixing_label_none
265 1484 : orig_full_labels(cur_indices(:)) = cur_labels(:)
266 :
267 : ! hysteretically set QM core from QM_core_list and radii, whole molecule
268 : ![NB] need to replace all the whole molecule stuff with pad to breakable bonds. not quite done
269 : ! (need intra molecule bond info, which isn't available for QM molecules yet)
270 :
271 : ! add core using hysteretic selection(core_list, r_core) + unbreakable bonds
272 : CALL add_layer_hysteretically( &
273 : nlist, particle_set, cell, nearest_dist, &
274 : orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
275 : force_mixing_label_QM_core_list, force_mixing_label_QM_core_list, force_mixing_label_QM_core, r_core, &
276 56 : max_n_qm, adaptive_exclude_molecules, molecule_set, broken_bonds)
277 : ![NB] should actually pass this back for making link sections?
278 56 : DEALLOCATE (broken_bonds)
279 :
280 56 : IF (debug_this_module .AND. output_unit > 0) THEN
281 0 : WRITE (output_unit, *) "BOB core new_indices ", new_indices(1:n_new)
282 0 : WRITE (output_unit, *) "BOB core new_labels ", new_labels(1:n_new)
283 : END IF
284 :
285 : ![NB] need more sophisticated QM extended, buffer rules
286 :
287 : ! add QM using hysteretic selection (core_list, r_qm) + unbreakable bonds
288 56 : IF (debug_this_module .AND. output_unit > 0) &
289 0 : WRITE (output_unit, *) "BOB QM_extended_seed_is_core_list ", QM_extended_seed_is_core_list
290 56 : IF (QM_extended_seed_is_core_list) THEN
291 0 : QM_extended_seed_min_label_val = force_mixing_label_QM_core_list
292 : ELSE ! QM region seed is all of core, not just core list + unbreakable bonds
293 56 : QM_extended_seed_min_label_val = force_mixing_label_QM_core
294 : END IF
295 : CALL add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, &
296 : orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
297 : QM_extended_seed_min_label_val, force_mixing_label_QM_core_list, &
298 : force_mixing_label_QM_dynamics, r_qm, &
299 56 : max_n_qm, adaptive_exclude_molecules, molecule_set)
300 :
301 56 : IF (debug_this_module .AND. output_unit > 0) THEN
302 0 : WRITE (output_unit, *) "BOB extended new_indices ", new_indices(1:n_new)
303 0 : WRITE (output_unit, *) "BOB extended new_labels ", new_labels(1:n_new)
304 : END IF
305 :
306 : ! add buffer using hysteretic selection (>= QM extended, r_buf) + unbreakable bonds
307 : CALL add_layer_hysteretically( &
308 : nlist, particle_set, cell, nearest_dist, &
309 : orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
310 : force_mixing_label_QM_dynamics, force_mixing_label_QM_core_list, force_mixing_label_buffer, r_buf, &
311 56 : max_n_qm, adaptive_exclude_molecules, molecule_set, broken_bonds)
312 : ![NB] should actually pass this back for making link sections?
313 56 : DEALLOCATE (broken_bonds)
314 :
315 56 : IF (debug_this_module .AND. output_unit > 0) THEN
316 0 : WRITE (output_unit, *) "BOB buffer new_indices ", new_indices(1:n_new)
317 0 : WRITE (output_unit, *) "BOB buffer new_labels ", new_labels(1:n_new)
318 : END IF
319 :
320 56 : DEALLOCATE (nearest_dist)
321 :
322 65568 : IF (PRESENT(labels_changed)) labels_changed = ANY(new_full_labels /= orig_full_labels)
323 :
324 56 : DEALLOCATE (orig_full_labels)
325 56 : DEALLOCATE (new_full_labels)
326 :
327 : ! reduce new indices, labels to actually used size
328 56 : CALL reallocate(new_indices, 1, n_new)
329 56 : CALL reallocate(new_labels, 1, n_new)
330 :
331 : ! save info in input structure
332 56 : restart_section => section_vals_get_subs_vals(qmmm_section, "FORCE_MIXING%RESTART_INFO")
333 56 : CALL section_vals_get(restart_section, explicit=explicit)
334 56 : IF (explicit) CALL section_vals_remove_values(restart_section)
335 56 : CALL section_vals_val_set(restart_section, "INDICES", i_vals_ptr=new_indices)
336 56 : CALL section_vals_val_set(restart_section, "LABELS", i_vals_ptr=new_labels)
337 :
338 56 : DEALLOCATE (cur_indices, cur_labels)
339 56 : CALL fist_neighbor_deallocate(nlist)
340 :
341 : ![NB] perhap be controlled by some &PRINT section?
342 56 : CALL cp_subsys_get(subsys, para_env=para_env)
343 56 : IF (para_env%is_source() .AND. output_unit > 0) THEN
344 : WRITE (unit=output_unit, fmt='(A,A,I6,A,I5,A,I5,A,I5)') &
345 28 : "QMMM FORCE MIXING final count (not including links): ", &
346 1021 : " N_QM core_list ", COUNT(new_labels == force_mixing_label_QM_core_list), &
347 1021 : " N_QM core ", COUNT(new_labels == force_mixing_label_QM_core), &
348 28 : " N_QM extended ", COUNT(new_labels == force_mixing_label_QM_dynamics .OR. &
349 1021 : new_labels == force_mixing_label_QM_dynamics_list), &
350 28 : " N_QM buffered ", COUNT(new_labels == force_mixing_label_buffer .OR. &
351 1049 : new_labels == force_mixing_label_buffer_list)
352 : END IF
353 :
354 280 : END SUBROUTINE update_force_mixing_labels
355 :
356 : ! **************************************************************************************************
357 : !> \brief ...
358 : !> \param ip ...
359 : !> \param label ...
360 : !> \param n_new ...
361 : !> \param new_indices ...
362 : !> \param new_labels ...
363 : !> \param new_full_labels ...
364 : !> \param max_n_qm ...
365 : ! **************************************************************************************************
366 1986 : SUBROUTINE add_new_label(ip, label, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
367 : INTEGER :: ip, label, n_new
368 : INTEGER, POINTER :: new_indices(:), new_labels(:)
369 : INTEGER :: new_full_labels(:), max_n_qm
370 :
371 : INTEGER :: i, old_index
372 :
373 1986 : IF (new_full_labels(ip) > force_mixing_label_none) THEN ! already marked, just change mark
374 0 : old_index = -1
375 0 : DO i = 1, n_new
376 0 : IF (new_indices(i) == ip) THEN
377 : old_index = i
378 : EXIT
379 : END IF
380 : END DO
381 0 : IF (old_index <= 0) &
382 : CALL cp_abort(__LOCATION__, &
383 : "add_new_label found atom with a label "// &
384 0 : "already set, but not in new_indices array")
385 0 : new_labels(old_index) = label
386 : ELSE
387 1986 : n_new = n_new + 1
388 1986 : IF (n_new > max_n_qm) &
389 : CALL cp_abort(__LOCATION__, &
390 : "add_new_label tried to add more atoms "// &
391 0 : "than allowed by &FORCE_MIXING&MAX_N_QM!")
392 1986 : IF (n_new > SIZE(new_indices)) CALL reallocate(new_indices, 1, n_new + 9)
393 1986 : IF (n_new > SIZE(new_labels)) CALL reallocate(new_labels, 1, n_new + 9)
394 1986 : new_indices(n_new) = ip
395 1986 : new_labels(n_new) = label
396 : END IF
397 1986 : new_full_labels(ip) = label
398 1986 : END SUBROUTINE add_new_label
399 :
400 : ! **************************************************************************************************
401 : !> \brief ...
402 : !> \param nlist ...
403 : !> \param particle_set ...
404 : !> \param cell ...
405 : !> \param nearest_dist ...
406 : !> \param orig_full_labels ...
407 : !> \param new_full_labels ...
408 : !> \param n_new ...
409 : !> \param new_indices ...
410 : !> \param new_labels ...
411 : !> \param seed_min_label_val ...
412 : !> \param seed_max_label_val ...
413 : !> \param set_label_val ...
414 : !> \param r_inout ...
415 : !> \param max_n_qm ...
416 : !> \param adaptive_exclude_molecules ...
417 : !> \param molecule_set ...
418 : !> \param broken_bonds ...
419 : ! **************************************************************************************************
420 168 : SUBROUTINE add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, &
421 168 : orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
422 : seed_min_label_val, seed_max_label_val, set_label_val, r_inout, max_n_qm, &
423 : adaptive_exclude_molecules, molecule_set, broken_bonds)
424 : TYPE(fist_neighbor_type), POINTER :: nlist
425 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
426 : TYPE(cell_type), POINTER :: cell
427 : REAL(dp) :: nearest_dist(:)
428 : INTEGER :: orig_full_labels(:), new_full_labels(:), &
429 : n_new
430 : INTEGER, POINTER :: new_indices(:), new_labels(:)
431 : INTEGER :: seed_min_label_val, seed_max_label_val, &
432 : set_label_val
433 : REAL(dp) :: r_inout(2)
434 : INTEGER :: max_n_qm
435 : CHARACTER(len=*), POINTER :: adaptive_exclude_molecules(:)
436 : TYPE(molecule_type), DIMENSION(:), OPTIONAL, &
437 : POINTER :: molecule_set
438 : INTEGER, OPTIONAL, POINTER :: broken_bonds(:)
439 :
440 : INTEGER :: i_ind, im, im_exclude, ip, ipair, &
441 : ipairkind, j_ind, output_unit
442 : LOGICAL :: adaptive_exclude, i_in_new_seed, i_outside_new_seed, j_in_new_seed, &
443 : j_outside_new_seed, molec_in_inner, molec_in_outer
444 : REAL(dp) :: r_ij(3), r_ij_mag
445 :
446 168 : output_unit = cp_logger_get_default_unit_nr()
447 :
448 168 : IF (debug_this_module .AND. output_unit > 0) WRITE (output_unit, *) "BOB adding hysteretically seed ", &
449 0 : seed_min_label_val, seed_max_label_val, " set ", set_label_val, " r ", r_inout
450 : ! calculate nearest dist from each atom outside of new seed to nearest atom inside of new seed
451 299706 : nearest_dist = HUGE(1.0_dp)
452 : ! loop over pairs of all kinds in random order
453 4704 : DO ipairkind = 1, SIZE(nlist%neighbor_kind_pairs)
454 2621184 : DO ipair = 1, nlist%neighbor_kind_pairs(ipairkind)%npairs
455 :
456 2616480 : i_ind = nlist%neighbor_kind_pairs(ipairkind)%list(1, ipair)
457 2616480 : j_ind = nlist%neighbor_kind_pairs(ipairkind)%list(2, ipair)
458 :
459 2616480 : i_in_new_seed = (new_full_labels(i_ind) >= seed_min_label_val .AND. new_full_labels(i_ind) <= seed_max_label_val)
460 2616480 : i_outside_new_seed = (new_full_labels(i_ind) < seed_min_label_val)
461 2616480 : j_in_new_seed = (new_full_labels(j_ind) >= seed_min_label_val .AND. new_full_labels(j_ind) <= seed_max_label_val)
462 2616480 : j_outside_new_seed = (new_full_labels(j_ind) < seed_min_label_val)
463 :
464 2621016 : IF ((i_in_new_seed .AND. j_outside_new_seed) .OR. (j_in_new_seed .AND. i_outside_new_seed)) THEN
465 22936 : r_ij = pbc(particle_set(i_ind)%r - particle_set(j_ind)%r, cell)
466 22936 : r_ij_mag = SQRT(SUM(r_ij**2))
467 5734 : IF (i_in_new_seed .AND. j_outside_new_seed .AND. (r_ij_mag < nearest_dist(j_ind))) THEN
468 1694 : nearest_dist(j_ind) = r_ij_mag
469 : END IF
470 5734 : IF (j_in_new_seed .AND. i_outside_new_seed .AND. (r_ij_mag < nearest_dist(i_ind))) THEN
471 1660 : nearest_dist(i_ind) = r_ij_mag
472 : END IF
473 : END IF
474 :
475 : END DO
476 : END DO
477 :
478 : ![NB] this is whole molecule. Should be replaced with labeling of individual atoms +
479 : ! pad_to_breakable_bonds (below), but QM molecule bond information isn't available yet
480 90558 : DO im = 1, SIZE(molecule_set)
481 : ! molecule_set(im)%first_atom,molecule_set(im)%last_atom
482 90390 : IF (ASSOCIATED(adaptive_exclude_molecules)) THEN
483 89730 : adaptive_exclude = .FALSE.
484 179460 : DO im_exclude = 1, SIZE(adaptive_exclude_molecules)
485 177120 : IF (TRIM(molecule_set(im)%molecule_kind%name) == TRIM(adaptive_exclude_molecules(im_exclude)) .OR. &
486 : TRIM(molecule_set(im)%molecule_kind%name) == '_QM_'//TRIM(adaptive_exclude_molecules(im_exclude))) &
487 92070 : adaptive_exclude = .TRUE.
488 : END DO
489 89730 : IF (adaptive_exclude) CYCLE
490 : END IF
491 350854 : molec_in_inner = ANY(nearest_dist(molecule_set(im)%first_atom:molecule_set(im)%last_atom) <= r_inout(1))
492 350504 : molec_in_outer = ANY(nearest_dist(molecule_set(im)%first_atom:molecule_set(im)%last_atom) <= r_inout(2))
493 88218 : IF (molec_in_inner) THEN
494 1480 : DO ip = molecule_set(im)%first_atom, molecule_set(im)%last_atom
495 : ! labels are being rebuild from scratch, so never overwrite new label that's higher level
496 1110 : IF (new_full_labels(ip) < set_label_val) &
497 1480 : CALL add_new_label(ip, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
498 : END DO
499 87680 : ELSE IF (molec_in_outer) THEN
500 520 : IF (ANY(orig_full_labels(molecule_set(im)%first_atom:molecule_set(im)%last_atom) >= set_label_val)) THEN
501 32 : DO ip = molecule_set(im)%first_atom, molecule_set(im)%last_atom
502 : ! labels are being rebuild from scratch, so never overwrite new label that's higher level
503 24 : IF (new_full_labels(ip) < set_label_val) &
504 32 : CALL add_new_label(ip, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
505 : END DO
506 : END IF
507 : END IF
508 : END DO
509 168 : IF (PRESENT(broken_bonds)) CALL reallocate(broken_bonds, 1, 0)
510 :
511 168 : END SUBROUTINE add_layer_hysteretically
512 :
513 : ! **************************************************************************************************
514 : !> \brief ...
515 : !> \param force_mixing_section ...
516 : !> \param subsys ...
517 : !> \param cell ...
518 : !> \param r_max ...
519 : !> \param nlist ...
520 : ! **************************************************************************************************
521 56 : SUBROUTINE make_neighbor_list(force_mixing_section, subsys, cell, r_max, nlist)
522 : TYPE(section_vals_type), POINTER :: force_mixing_section
523 : TYPE(cp_subsys_type), POINTER :: subsys
524 : TYPE(cell_type), POINTER :: cell
525 : REAL(dp) :: r_max
526 : TYPE(fist_neighbor_type), POINTER :: nlist
527 :
528 : CHARACTER(LEN=default_string_length) :: kind_name
529 56 : CHARACTER(LEN=default_string_length), POINTER :: kind_name_a(:)
530 : INTEGER :: ik
531 : LOGICAL :: skip_kind
532 56 : REAL(dp), ALLOCATABLE :: r_max_a(:, :), r_minsq_a(:, :)
533 : TYPE(atomic_kind_type), POINTER :: atomic_kind
534 :
535 224 : ALLOCATE (r_max_a(SIZE(subsys%atomic_kinds%els), SIZE(subsys%atomic_kinds%els)))
536 168 : ALLOCATE (r_minsq_a(SIZE(subsys%atomic_kinds%els), SIZE(subsys%atomic_kinds%els)))
537 472612 : r_max_a = r_max
538 472612 : r_minsq_a = EPSILON(1.0_dp)
539 :
540 : ! save kind names
541 168 : ALLOCATE (kind_name_a(SIZE(subsys%atomic_kinds%els)))
542 2012 : DO ik = 1, SIZE(subsys%atomic_kinds%els)
543 1956 : atomic_kind => subsys%atomic_kinds%els(ik)
544 1956 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=kind_name)
545 2012 : kind_name_a(ik) = kind_name
546 : END DO
547 :
548 : ! overwrite kind names so that none are QM, and so excluding QM-QM interactions
549 : ! (which is not what we want) will not happen
550 2012 : DO ik = 1, SIZE(subsys%atomic_kinds%els)
551 1956 : atomic_kind => subsys%atomic_kinds%els(ik)
552 1956 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=kind_name)
553 : ! when atom is QM atom, kind_name is replaced with original
554 : ! mm kind name, and return status is logical .TRUE.
555 1956 : skip_kind = qmmm_ff_precond_only_qm(kind_name)
556 2012 : CALL set_atomic_kind(atomic_kind=atomic_kind, name=kind_name)
557 : END DO
558 :
559 56 : NULLIFY (nlist)
560 : CALL build_fist_neighbor_lists(subsys%atomic_kinds%els, subsys%particles%els, &
561 : cell=cell, r_max=r_max_a, r_minsq=r_minsq_a, &
562 : ei_scale14=1.0_dp, vdw_scale14=1.0_dp, nonbonded=nlist, &
563 : para_env=subsys%para_env, build_from_scratch=.TRUE., geo_check=.FALSE., &
564 56 : mm_section=force_mixing_section)
565 :
566 56 : DEALLOCATE (r_max_a, r_minsq_a)
567 :
568 : ! restore kind names
569 2012 : DO ik = 1, SIZE(subsys%atomic_kinds%els)
570 2012 : CALL set_atomic_kind(atomic_kind=atomic_kind, name=kind_name_a(ik))
571 : END DO
572 56 : DEALLOCATE (kind_name_a)
573 :
574 56 : END SUBROUTINE make_neighbor_list
575 :
576 : ! **************************************************************************************************
577 : !> \brief ...
578 : !> \param subsys ...
579 : !> \param qmmm_section ...
580 : !> \param qmmm_core_section ...
581 : !> \param qmmm_extended_section ...
582 : !> \par History
583 : !> 02.2012 created [noam]
584 : !> \author Noam Bernstein
585 : ! **************************************************************************************************
586 30 : SUBROUTINE setup_force_mixing_qmmm_sections(subsys, qmmm_section, qmmm_core_section, qmmm_extended_section)
587 : TYPE(cp_subsys_type), POINTER :: subsys
588 : TYPE(section_vals_type), POINTER :: qmmm_section, qmmm_core_section, &
589 : qmmm_extended_section
590 :
591 30 : CHARACTER(len=default_string_length), POINTER :: elem_mapping(:, :), elem_mapping_entry(:)
592 : INTEGER :: delta_charge, i_rep_section_core, i_rep_section_extended, i_rep_val_core, &
593 : i_rep_val_extended, ielem, ip, n_elements, output_unit
594 30 : INTEGER, POINTER :: cur_indices(:), cur_labels(:)
595 : LOGICAL :: mapped, new_element_core, &
596 : new_element_extended
597 30 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
598 : TYPE(section_vals_type), POINTER :: buffer_non_adaptive_section, &
599 : dup_link_section, &
600 : force_mixing_section, link_section, &
601 : qm_kind_section
602 :
603 30 : NULLIFY (qmmm_core_section, qmmm_extended_section)
604 60 : output_unit = cp_logger_get_default_unit_nr()
605 :
606 : ! create new qmmm sections for core and extended
607 30 : CALL section_vals_duplicate(qmmm_section, qmmm_core_section)
608 30 : CALL section_vals_duplicate(qmmm_section, qmmm_extended_section)
609 :
610 : ! remove LINKs (specified by user for core) from extended
611 30 : link_section => section_vals_get_subs_vals(qmmm_extended_section, "LINK", can_return_null=.TRUE.)
612 30 : IF (ASSOCIATED(link_section)) THEN
613 30 : CALL section_vals_remove_values(link_section)
614 : END IF
615 : ! for LINKs to be added to extended
616 : buffer_non_adaptive_section => section_vals_get_subs_vals(qmmm_extended_section, "FORCE_MIXING%BUFFER_NON_ADAPTIVE", &
617 30 : can_return_null=.TRUE.)
618 30 : link_section => section_vals_get_subs_vals(buffer_non_adaptive_section, "LINK", can_return_null=.TRUE.)
619 30 : IF (ASSOCIATED(link_section)) THEN
620 30 : NULLIFY (dup_link_section)
621 30 : CALL section_vals_duplicate(link_section, dup_link_section)
622 30 : CALL section_vals_set_subs_vals(qmmm_extended_section, "LINK", dup_link_section)
623 30 : CALL section_vals_release(dup_link_section)
624 : END IF
625 :
626 30 : IF (debug_this_module .AND. output_unit > 0) THEN
627 0 : link_section => section_vals_get_subs_vals(qmmm_core_section, "LINK", can_return_null=.TRUE.)
628 0 : WRITE (output_unit, *) "core section has LINKs ", ASSOCIATED(link_section)
629 0 : CALL section_vals_write(link_section, unit_nr=6)
630 0 : link_section => section_vals_get_subs_vals(qmmm_extended_section, "LINK", can_return_null=.TRUE.)
631 0 : WRITE (output_unit, *) "extended section has LINKs ", ASSOCIATED(link_section)
632 0 : CALL section_vals_write(link_section, unit_nr=6)
633 : END IF
634 :
635 30 : force_mixing_section => section_vals_get_subs_vals(qmmm_section, "FORCE_MIXING")
636 :
637 : ! get QM_KIND_ELEMENT_MAPPING
638 30 : CALL section_vals_val_get(force_mixing_section, "QM_KIND_ELEMENT_MAPPING", n_rep_val=n_elements)
639 90 : ALLOCATE (elem_mapping(2, n_elements))
640 102 : DO ielem = 1, n_elements
641 72 : CALL section_vals_val_get(force_mixing_section, "QM_KIND_ELEMENT_MAPPING", i_rep_val=ielem, c_vals=elem_mapping_entry)
642 390 : elem_mapping(1:2, ielem) = elem_mapping_entry(1:2)
643 : END DO
644 :
645 : ! get CUR_INDICES, CUR_LABELS
646 30 : CALL get_force_mixing_indices(force_mixing_section, cur_indices, cur_labels)
647 30 : IF (SIZE(cur_indices) <= 0) &
648 0 : CPABORT("cur_indices is empty, found no QM atoms")
649 :
650 30 : IF (debug_this_module .AND. output_unit > 0) THEN
651 0 : WRITE (output_unit, *) "cur_indices ", cur_indices
652 0 : WRITE (output_unit, *) "cur_labels ", cur_labels
653 : END IF
654 :
655 : ! loop through elements and atoms, and set up new QM_KIND sections
656 30 : particles => subsys%particles%els
657 :
658 1014 : DO ip = 1, SIZE(cur_indices)
659 984 : IF (cur_labels(ip) > force_mixing_label_none .AND. cur_labels(ip) < force_mixing_label_QM_core_list .AND. &
660 30 : cur_labels(ip) /= force_mixing_label_termination) THEN
661 892 : mapped = .FALSE.
662 1472 : DO ielem = 1, n_elements
663 1472 : IF (TRIM(particles(cur_indices(ip))%atomic_kind%element_symbol) == TRIM(elem_mapping(1, ielem))) THEN
664 : mapped = .TRUE.
665 : EXIT
666 : END IF
667 : END DO
668 892 : IF (.NOT. mapped) &
669 : CALL cp_abort(__LOCATION__, &
670 : "Force-mixing failed to find QM_KIND mapping for atom of type "// &
671 : TRIM(particles(cur_indices(ip))%atomic_kind%element_symbol)// &
672 0 : "! ")
673 : END IF
674 : END DO
675 :
676 : ! pre-existing QM_KIND section specifies list of core atom
677 30 : qm_kind_section => section_vals_get_subs_vals3(qmmm_section, "QM_KIND")
678 30 : CALL section_vals_get(qm_kind_section, n_repetition=i_rep_section_core)
679 30 : IF (i_rep_section_core <= 0) &
680 : CALL cp_abort(__LOCATION__, &
681 : "Force-mixing QM didn't find any QM_KIND sections, "// &
682 0 : "so no core specified!")
683 30 : i_rep_section_extended = i_rep_section_core
684 102 : DO ielem = 1, n_elements
685 : new_element_core = .TRUE.
686 : new_element_extended = .TRUE.
687 3522 : DO ip = 1, SIZE(cur_indices) ! particles with label
688 3420 : IF (TRIM(particles(cur_indices(ip))%atomic_kind%element_symbol) /= TRIM(elem_mapping(1, ielem))) CYCLE
689 : ! extended
690 : ! if current particle is some sort of QM atom, and not in core list
691 : ! (those the user gave explicit QM_KIND sections for), and not a
692 : ! termination atom, need to make a QM_KIND section for it
693 : IF (cur_labels(ip) > force_mixing_label_none .AND. &
694 984 : cur_labels(ip) /= force_mixing_label_QM_core_list .AND. &
695 : cur_labels(ip) /= force_mixing_label_termination) THEN
696 892 : qm_kind_section => section_vals_get_subs_vals3(qmmm_extended_section, "QM_KIND")
697 892 : IF (new_element_extended) THEN ! add new QM_KIND section for this element
698 72 : i_rep_section_extended = i_rep_section_extended + 1
699 72 : CALL section_vals_add_values(qm_kind_section)
700 : CALL section_vals_val_set(qm_kind_section, "_SECTION_PARAMETERS_", i_rep_section=i_rep_section_extended, &
701 72 : c_val=elem_mapping(2, ielem))
702 72 : i_rep_val_extended = 0
703 72 : new_element_extended = .FALSE.
704 : END IF
705 892 : i_rep_val_extended = i_rep_val_extended + 1
706 : CALL section_vals_val_set(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section_extended, &
707 892 : i_rep_val=i_rep_val_extended, i_val=cur_indices(ip))
708 : END IF ! is a non-termination QM atom
709 :
710 : ! core
711 : ! if current particle is a core QM atom, and not in core list (those the user
712 : ! gave explicit QM_KIND sections for, need to make a QM_KIND section for it
713 1056 : IF (cur_labels(ip) == force_mixing_label_QM_core) THEN
714 78 : qm_kind_section => section_vals_get_subs_vals3(qmmm_core_section, "QM_KIND")
715 78 : IF (new_element_core) THEN ! add new QM_KIND section for this element
716 48 : i_rep_section_core = i_rep_section_core + 1
717 48 : CALL section_vals_add_values(qm_kind_section)
718 : CALL section_vals_val_set(qm_kind_section, "_SECTION_PARAMETERS_", i_rep_section=i_rep_section_core, &
719 48 : c_val=elem_mapping(2, ielem))
720 48 : i_rep_val_core = 0
721 48 : new_element_core = .FALSE.
722 : END IF
723 78 : i_rep_val_core = i_rep_val_core + 1
724 : CALL section_vals_val_set(qm_kind_section, "MM_INDEX", i_rep_section=i_rep_section_core, &
725 78 : i_rep_val=i_rep_val_core, i_val=cur_indices(ip))
726 : END IF ! is a non-termination QM atom
727 :
728 : END DO ! atom index ip
729 : END DO ! element index ielem
730 :
731 30 : CALL section_vals_val_get(force_mixing_section, "EXTENDED_DELTA_CHARGE", i_val=delta_charge)
732 30 : CALL section_vals_val_set(qmmm_extended_section, "DELTA_CHARGE", i_val=delta_charge)
733 :
734 : ![NB] check
735 30 : DEALLOCATE (elem_mapping, cur_indices, cur_labels)
736 :
737 30 : IF (debug_this_module .AND. output_unit > 0) THEN
738 0 : WRITE (output_unit, *) "qmmm_core_section"
739 0 : CALL section_vals_write(qmmm_core_section, unit_nr=6)
740 0 : WRITE (output_unit, *) "qmmm_extended_section"
741 0 : CALL section_vals_write(qmmm_extended_section, unit_nr=6)
742 : END IF
743 :
744 120 : END SUBROUTINE setup_force_mixing_qmmm_sections
745 :
746 : ! **************************************************************************************************
747 : !> \brief ...
748 : !> \param force_mixing_section ...
749 : !> \param indices ...
750 : !> \param labels ...
751 : ! **************************************************************************************************
752 86 : SUBROUTINE get_force_mixing_indices(force_mixing_section, indices, labels)
753 : TYPE(section_vals_type), POINTER :: force_mixing_section
754 : INTEGER, POINTER :: indices(:), labels(:)
755 :
756 : INTEGER :: i_rep_val, n_indices, n_labels, n_reps
757 86 : INTEGER, POINTER :: indices_entry(:), labels_entry(:)
758 : LOGICAL :: explicit
759 : TYPE(section_vals_type), POINTER :: restart_section
760 :
761 86 : NULLIFY (indices, labels)
762 172 : restart_section => section_vals_get_subs_vals(force_mixing_section, "RESTART_INFO")
763 86 : CALL section_vals_get(restart_section, explicit=explicit)
764 86 : IF (.NOT. explicit) THEN ! no old indices, labels, return empty arrays
765 8 : ALLOCATE (indices(0))
766 8 : ALLOCATE (labels(0))
767 8 : RETURN
768 : END IF
769 :
770 : ![NB] maybe switch to reallocatable array
771 78 : CALL section_vals_val_get(restart_section, "INDICES", n_rep_val=n_reps)
772 78 : n_indices = 0
773 156 : DO i_rep_val = 1, n_reps
774 : CALL section_vals_val_get(restart_section, "INDICES", &
775 78 : i_rep_val=i_rep_val, i_vals=indices_entry)
776 156 : n_indices = n_indices + SIZE(indices_entry)
777 : END DO
778 234 : ALLOCATE (indices(n_indices))
779 78 : n_indices = 0
780 156 : DO i_rep_val = 1, n_reps
781 : CALL section_vals_val_get(restart_section, "INDICES", &
782 78 : i_rep_val=i_rep_val, i_vals=indices_entry)
783 4902 : indices(n_indices + 1:n_indices + SIZE(indices_entry)) = indices_entry
784 156 : n_indices = n_indices + SIZE(indices_entry)
785 : END DO
786 :
787 78 : CALL section_vals_val_get(restart_section, "LABELS", n_rep_val=n_reps)
788 78 : n_labels = 0
789 156 : DO i_rep_val = 1, n_reps
790 : CALL section_vals_val_get(restart_section, "LABELS", &
791 78 : i_rep_val=i_rep_val, i_vals=labels_entry)
792 156 : n_labels = n_labels + SIZE(labels_entry)
793 : END DO
794 234 : ALLOCATE (labels(n_labels))
795 78 : n_labels = 0
796 156 : DO i_rep_val = 1, n_reps
797 : CALL section_vals_val_get(restart_section, "LABELS", &
798 78 : i_rep_val=i_rep_val, i_vals=labels_entry)
799 4902 : labels(n_labels + 1:n_labels + SIZE(labels_entry)) = labels_entry
800 156 : n_labels = n_labels + SIZE(labels_entry)
801 : END DO
802 :
803 78 : IF (n_indices /= n_labels) &
804 0 : CPABORT("got unequal numbers of force_mixing indices and labels!")
805 242 : END SUBROUTINE get_force_mixing_indices
806 :
807 : END MODULE qmmmx_util
|
