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 Initialize a QM/MM calculation
10 : !> \par History
11 : !> 5.2004 created [fawzi]
12 : !> \author Fawzi Mohamed
13 : ! **************************************************************************************************
14 : MODULE qmmm_init
15 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind,&
18 : set_atomic_kind
19 : USE cell_methods, ONLY: read_cell
20 : USE cell_types, ONLY: cell_type,&
21 : use_perd_xyz
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_eri_mme_interface, ONLY: cp_eri_mme_init_read_input,&
24 : cp_eri_mme_set_params
25 : USE cp_log_handling, ONLY: cp_get_default_logger,&
26 : cp_logger_type
27 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
28 : cp_print_key_unit_nr
29 : USE cp_subsys_types, ONLY: cp_subsys_get,&
30 : cp_subsys_type
31 : USE cp_units, ONLY: cp_unit_from_cp2k,&
32 : cp_unit_to_cp2k
33 : USE external_potential_types, ONLY: fist_potential_type,&
34 : get_potential,&
35 : set_potential
36 : USE force_field_types, ONLY: input_info_type
37 : USE force_fields_input, ONLY: read_gd_section,&
38 : read_gp_section,&
39 : read_lj_section,&
40 : read_wl_section
41 : USE input_constants, ONLY: &
42 : RADIUS_QMMM_DEFAULT, calc_always, calc_once, do_eri_mme, do_qmmm_gauss, &
43 : do_qmmm_image_calcmatrix, do_qmmm_image_iter, do_qmmm_link_gho, do_qmmm_link_imomm, &
44 : do_qmmm_link_pseudo, do_qmmm_pcharge, do_qmmm_swave
45 : USE input_section_types, ONLY: section_vals_get,&
46 : section_vals_get_subs_vals,&
47 : section_vals_type,&
48 : section_vals_val_get
49 : USE kinds, ONLY: default_string_length,&
50 : dp
51 : USE message_passing, ONLY: mp_para_env_type
52 : USE molecule_kind_types, ONLY: molecule_kind_type
53 : USE pair_potential_types, ONLY: pair_potential_reallocate
54 : USE particle_list_types, ONLY: particle_list_type
55 : USE particle_types, ONLY: particle_type
56 : USE pw_env_types, ONLY: pw_env_type
57 : USE qmmm_elpot, ONLY: qmmm_potential_init
58 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
59 : USE qmmm_gaussian_init, ONLY: qmmm_gaussian_initialize
60 : USE qmmm_per_elpot, ONLY: qmmm_ewald_potential_init,&
61 : qmmm_per_potential_init
62 : USE qmmm_types_low, ONLY: add_set_type,&
63 : add_shell_type,&
64 : create_add_set_type,&
65 : create_add_shell_type,&
66 : qmmm_env_mm_type,&
67 : qmmm_env_qm_type,&
68 : qmmm_links_type
69 : USE qs_environment_types, ONLY: get_qs_env,&
70 : qs_environment_type
71 : USE shell_potential_types, ONLY: get_shell,&
72 : shell_kind_type
73 : #include "./base/base_uses.f90"
74 :
75 : IMPLICIT NONE
76 : PRIVATE
77 :
78 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
79 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_init'
80 :
81 : PUBLIC :: assign_mm_charges_and_radius, &
82 : print_qmmm_charges, &
83 : print_qmmm_links, &
84 : print_image_charge_info, &
85 : qmmm_init_gaussian_type, &
86 : qmmm_init_potential, &
87 : qmmm_init_periodic_potential, &
88 : setup_qmmm_vars_qm, &
89 : setup_qmmm_vars_mm, &
90 : setup_qmmm_links, &
91 : move_or_add_atoms, &
92 : setup_origin_mm_cell
93 :
94 : CONTAINS
95 :
96 : ! **************************************************************************************************
97 : !> \brief Assigns charges and radius to evaluate the MM electrostatic potential
98 : !> \param subsys the subsys containing the MM charges
99 : !> \param charges ...
100 : !> \param mm_atom_chrg ...
101 : !> \param mm_el_pot_radius ...
102 : !> \param mm_el_pot_radius_corr ...
103 : !> \param mm_atom_index ...
104 : !> \param mm_link_atoms ...
105 : !> \param mm_link_scale_factor ...
106 : !> \param added_shells ...
107 : !> \param shell_model ...
108 : !> \par History
109 : !> 06.2004 created [tlaino]
110 : !> \author Teodoro Laino
111 : ! **************************************************************************************************
112 394 : SUBROUTINE assign_mm_charges_and_radius(subsys, charges, mm_atom_chrg, mm_el_pot_radius, &
113 : mm_el_pot_radius_corr, mm_atom_index, mm_link_atoms, &
114 : mm_link_scale_factor, added_shells, shell_model)
115 : TYPE(cp_subsys_type), POINTER :: subsys
116 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
117 : REAL(dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, &
118 : mm_el_pot_radius_corr
119 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index, mm_link_atoms
120 : REAL(dp), DIMENSION(:), POINTER :: mm_link_scale_factor
121 : TYPE(add_shell_type), OPTIONAL, POINTER :: added_shells
122 : LOGICAL :: shell_model
123 :
124 : INTEGER :: I, ilink, IndMM, IndShell, ishell
125 : LOGICAL :: is_shell
126 : REAL(dp) :: qcore, qi, qshell, rc, ri
127 : TYPE(atomic_kind_type), POINTER :: my_kind
128 : TYPE(fist_potential_type), POINTER :: my_potential
129 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
130 : shell_particles
131 394 : TYPE(particle_type), DIMENSION(:), POINTER :: core_set, particle_set, shell_set
132 : TYPE(shell_kind_type), POINTER :: shell_kind
133 :
134 394 : NULLIFY (particle_set, my_kind, added_shells)
135 : CALL cp_subsys_get(subsys=subsys, particles=particles, core_particles=core_particles, &
136 394 : shell_particles=shell_particles)
137 394 : particle_set => particles%els
138 :
139 190762 : IF (ALL(particle_set(:)%shell_index == 0)) THEN
140 392 : shell_model = .FALSE.
141 392 : CALL create_add_shell_type(added_shells, ndim=0)
142 : ELSE
143 2 : shell_model = .TRUE.
144 : END IF
145 :
146 394 : IF (shell_model) THEN
147 2 : shell_set => shell_particles%els
148 2 : core_set => core_particles%els
149 2 : ishell = SIZE(shell_set)
150 2 : CALL create_add_shell_type(added_shells, ndim=ishell)
151 2 : added_shells%added_particles => shell_set
152 2 : added_shells%added_cores => core_set
153 : END IF
154 :
155 188174 : DO I = 1, SIZE(mm_atom_index)
156 187780 : IndMM = mm_atom_index(I)
157 187780 : my_kind => particle_set(IndMM)%atomic_kind
158 : CALL get_atomic_kind(atomic_kind=my_kind, fist_potential=my_potential, &
159 187780 : shell_active=is_shell, shell=shell_kind)
160 : CALL get_potential(potential=my_potential, &
161 : qeff=qi, &
162 : qmmm_radius=ri, &
163 187780 : qmmm_corr_radius=rc)
164 187780 : IF (ASSOCIATED(charges)) qi = charges(IndMM)
165 187780 : mm_atom_chrg(I) = qi
166 187780 : mm_el_pot_radius(I) = ri
167 187780 : mm_el_pot_radius_corr(I) = rc
168 375954 : IF (is_shell) THEN
169 56 : IndShell = particle_set(IndMM)%shell_index
170 56 : IF (ASSOCIATED(shell_kind)) THEN
171 : CALL get_shell(shell=shell_kind, &
172 : charge_core=qcore, &
173 56 : charge_shell=qshell)
174 56 : mm_atom_chrg(I) = qcore
175 : END IF
176 56 : added_shells%mm_core_index(IndShell) = IndMM
177 56 : added_shells%mm_core_chrg(IndShell) = qshell
178 56 : added_shells%mm_el_pot_radius(Indshell) = ri*1.0_dp
179 56 : added_shells%mm_el_pot_radius_corr(Indshell) = rc*1.0_dp
180 : END IF
181 : END DO
182 :
183 394 : IF (ASSOCIATED(mm_link_atoms)) THEN
184 256 : DO ilink = 1, SIZE(mm_link_atoms)
185 40710 : DO i = 1, SIZE(mm_atom_index)
186 40710 : IF (mm_atom_index(i) == mm_link_atoms(ilink)) EXIT
187 : END DO
188 194 : IndMM = mm_atom_index(I)
189 256 : mm_atom_chrg(i) = mm_atom_chrg(i)*mm_link_scale_factor(ilink)
190 : END DO
191 : END IF
192 :
193 394 : END SUBROUTINE assign_mm_charges_and_radius
194 :
195 : ! **************************************************************************************************
196 : !> \brief Print info on charges generating the qmmm potential..
197 : !> \param mm_atom_index ...
198 : !> \param mm_atom_chrg ...
199 : !> \param mm_el_pot_radius ...
200 : !> \param mm_el_pot_radius_corr ...
201 : !> \param added_charges ...
202 : !> \param added_shells ...
203 : !> \param qmmm_section ...
204 : !> \param nocompatibility ...
205 : !> \param shell_model ...
206 : !> \par History
207 : !> 01.2005 created [tlaino]
208 : !> \author Teodoro Laino
209 : ! **************************************************************************************************
210 394 : SUBROUTINE print_qmmm_charges(mm_atom_index, mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, &
211 : added_charges, added_shells, qmmm_section, nocompatibility, shell_model)
212 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
213 : REAL(dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, &
214 : mm_el_pot_radius_corr
215 : TYPE(add_set_type), POINTER :: added_charges
216 : TYPE(add_shell_type), POINTER :: added_shells
217 : TYPE(section_vals_type), POINTER :: qmmm_section
218 : LOGICAL, INTENT(IN) :: nocompatibility, shell_model
219 :
220 : INTEGER :: I, ind1, ind2, IndMM, iw
221 : REAL(KIND=dp) :: qi, qtot, rc, ri
222 : TYPE(cp_logger_type), POINTER :: logger
223 :
224 394 : qtot = 0.0_dp
225 394 : logger => cp_get_default_logger()
226 : iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%QMMM_CHARGES", &
227 394 : extension=".log")
228 394 : IF (iw > 0) THEN
229 169 : WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
230 169 : WRITE (iw, FMT='(/5X,A)') "MM POINT CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
231 169 : WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
232 23623 : DO I = 1, SIZE(mm_atom_index)
233 23454 : IndMM = mm_atom_index(I)
234 23454 : qi = mm_atom_chrg(I)
235 23454 : qtot = qtot + qi
236 23454 : ri = mm_el_pot_radius(I)
237 23454 : rc = mm_el_pot_radius_corr(I)
238 23623 : IF (nocompatibility) THEN
239 2177 : WRITE (iw, '(5X,A9,T15,I5,T28,A8,T38,F12.6,T60,A8,T69,F12.6)') ' MM ATOM:', IndMM, ' RADIUS:', ri, &
240 4354 : ' CHARGE:', qi
241 : ELSE
242 : WRITE (iw, '(5X,A9,T15,I5,T28,A8,T38,F12.6,T60,A8,T69,F12.6,/,T56,A12,T69,F12.6)') &
243 21277 : ' MM ATOM:', IndMM, ' RADIUS:', ri, ' CHARGE:', qi, 'CORR. RADIUS', rc
244 : END IF
245 : END DO
246 169 : IF (added_charges%num_mm_atoms /= 0) THEN
247 4 : WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
248 4 : WRITE (iw, '(/5X,A)') "ADDED POINT CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
249 4 : WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
250 18 : DO I = 1, SIZE(added_charges%mm_atom_index)
251 14 : IndMM = added_charges%mm_atom_index(I)
252 14 : qi = added_charges%mm_atom_chrg(I)
253 14 : qtot = qtot + qi
254 14 : ri = added_charges%mm_el_pot_radius(I)
255 14 : ind1 = added_charges%add_env(I)%Index1
256 14 : ind2 = added_charges%add_env(I)%Index2
257 18 : IF (nocompatibility) THEN
258 14 : WRITE (iw, '(5X,A9,I5,T25,A8,T35,F12.6,T50,A8,T59,F12.6,I5,I5)') 'MM POINT:', IndMM, ' RADIUS:', ri, &
259 28 : ' CHARGE:', qi, ind1, ind2
260 : ELSE
261 : WRITE (iw, '(5X,A9,I5,T25,A8,T35,F12.6,T50,A8,T59,F12.6,I5,I5,/,T56,A12,T69,F12.6)') &
262 0 : 'MM POINT:', IndMM, ' RADIUS:', ri, ' CHARGE:', qi, ind1, ind2, 'CORR. RADIUS', rc
263 : END IF
264 : END DO
265 :
266 : END IF
267 :
268 169 : IF (shell_model) THEN
269 1 : WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 73)
270 1 : WRITE (iw, '(/5X,A)') "ADDED SHELL CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
271 1 : WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 73)
272 :
273 29 : DO I = 1, SIZE(added_shells%mm_core_index)
274 28 : IndMM = added_shells%mm_core_index(I)
275 28 : qi = added_shells%mm_core_chrg(I)
276 28 : qtot = qtot + qi
277 28 : ri = added_shells%mm_el_pot_radius(I)
278 29 : IF (nocompatibility) THEN
279 28 : WRITE (iw, '(7X,A,I5,A8,F12.6,A8,F12.6,3F12.6)') 'SHELL:', IndMM, ' RADIUS:', ri, &
280 140 : ' CHARGE:', qi, added_shells%added_particles(I)%r
281 : ELSE
282 0 : WRITE (iw, '(7X,A,I5,A8,F12.6,A8,F12.6,A,F12.6)') 'SHELL:', IndMM, ' RADIUS:', ri, &
283 0 : ' CHARGE:', qi, ' CORR. RADIUS', rc
284 : END IF
285 :
286 : END DO
287 :
288 : END IF
289 :
290 169 : WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
291 169 : WRITE (iw, '(/,T50,A,T69,F12.6)') ' TOTAL CHARGE:', qtot
292 169 : WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 79)
293 : END IF
294 : CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
295 394 : "PRINT%QMMM_CHARGES")
296 394 : END SUBROUTINE print_qmmm_charges
297 :
298 : ! **************************************************************************************************
299 : !> \brief Print info on qm/mm links
300 : !> \param qmmm_section ...
301 : !> \param qmmm_links ...
302 : !> \par History
303 : !> 01.2005 created [tlaino]
304 : !> \author Teodoro Laino
305 : ! **************************************************************************************************
306 64 : SUBROUTINE print_qmmm_links(qmmm_section, qmmm_links)
307 : TYPE(section_vals_type), POINTER :: qmmm_section
308 : TYPE(qmmm_links_type), POINTER :: qmmm_links
309 :
310 : INTEGER :: i, iw, mm_index, qm_index
311 : REAL(KIND=dp) :: alpha
312 : TYPE(cp_logger_type), POINTER :: logger
313 :
314 64 : logger => cp_get_default_logger()
315 64 : iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%qmmm_link_info", extension=".log")
316 64 : IF (iw > 0) THEN
317 21 : IF (ASSOCIATED(qmmm_links)) THEN
318 21 : WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
319 21 : WRITE (iw, FMT="(/,T31,A)") " QM/MM LINKS "
320 21 : WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
321 21 : IF (ASSOCIATED(qmmm_links%imomm)) THEN
322 20 : WRITE (iw, FMT="(/,T31,A)") " IMOMM TYPE LINK "
323 56 : DO I = 1, SIZE(qmmm_links%imomm)
324 36 : qm_index = qmmm_links%imomm(I)%link%qm_index
325 36 : mm_index = qmmm_links%imomm(I)%link%mm_index
326 36 : alpha = qmmm_links%imomm(I)%link%alpha
327 36 : WRITE (iw, FMT="(T2,A,T20,A9,I8,1X,A9,I8,T62,A6,F12.6)") "TYPE: IMOMM", &
328 92 : "QM INDEX:", qm_index, "MM INDEX:", mm_index, "ALPHA:", alpha
329 : END DO
330 : END IF
331 21 : IF (ASSOCIATED(qmmm_links%pseudo)) THEN
332 1 : WRITE (iw, FMT="(/,T31,A)") " PSEUDO TYPE LINK "
333 3 : DO I = 1, SIZE(qmmm_links%pseudo)
334 2 : qm_index = qmmm_links%pseudo(I)%link%qm_index
335 2 : mm_index = qmmm_links%pseudo(I)%link%mm_index
336 2 : WRITE (iw, FMT="(T2,A,T20,A9,I8,1X,A9,I8)") "TYPE: PSEUDO", &
337 5 : "QM INDEX:", qm_index, "MM INDEX:", mm_index
338 : END DO
339 : END IF
340 21 : WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 73)
341 : ELSE
342 0 : WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
343 0 : WRITE (iw, FMT="(/,T26,A)") " NO QM/MM LINKS DETECTED"
344 0 : WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
345 : END IF
346 : END IF
347 : CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
348 64 : "PRINT%qmmm_link_info")
349 64 : END SUBROUTINE print_qmmm_links
350 :
351 : ! **************************************************************************************************
352 : !> \brief ...
353 : !> \param qmmm_env_qm ...
354 : !> \param para_env ...
355 : !> \param mm_atom_chrg ...
356 : !> \param qs_env ...
357 : !> \param added_charges ...
358 : !> \param added_shells ...
359 : !> \param print_section ...
360 : !> \param qmmm_section ...
361 : !> \par History
362 : !> 1.2005 created [tlaino]
363 : !> \author Teodoro Laino
364 : ! **************************************************************************************************
365 394 : SUBROUTINE qmmm_init_gaussian_type(qmmm_env_qm, para_env, &
366 : mm_atom_chrg, qs_env, added_charges, added_shells, &
367 : print_section, qmmm_section)
368 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm
369 : TYPE(mp_para_env_type), POINTER :: para_env
370 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_atom_chrg
371 : TYPE(qs_environment_type), POINTER :: qs_env
372 : TYPE(add_set_type), POINTER :: added_charges
373 : TYPE(add_shell_type), POINTER :: added_shells
374 : TYPE(section_vals_type), POINTER :: print_section, qmmm_section
375 :
376 : INTEGER :: i
377 : REAL(KIND=dp) :: maxchrg
378 394 : REAL(KIND=dp), DIMENSION(:), POINTER :: maxradius, maxradius2
379 : TYPE(pw_env_type), POINTER :: pw_env
380 :
381 394 : NULLIFY (maxradius, maxradius2, pw_env)
382 :
383 188568 : maxchrg = MAXVAL(ABS(mm_atom_chrg(:)))
384 394 : CALL get_qs_env(qs_env, pw_env=pw_env)
385 410 : IF (qmmm_env_qm%add_mm_charges) maxchrg = MAX(maxchrg, MAXVAL(ABS(added_charges%mm_atom_chrg(:))))
386 : CALL qmmm_gaussian_initialize(qmmm_gaussian_fns=qmmm_env_qm%pgfs, &
387 : para_env=para_env, &
388 : pw_env=pw_env, &
389 : mm_el_pot_radius=qmmm_env_qm%mm_el_pot_radius, &
390 : mm_el_pot_radius_corr=qmmm_env_qm%mm_el_pot_radius_corr, &
391 : qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
392 : eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
393 : maxradius=maxradius, &
394 : maxchrg=maxchrg, &
395 : compatibility=qmmm_env_qm%compatibility, &
396 : print_section=print_section, &
397 394 : qmmm_section=qmmm_section)
398 :
399 394 : IF (qmmm_env_qm%move_mm_charges .OR. qmmm_env_qm%add_mm_charges) THEN
400 : CALL qmmm_gaussian_initialize(qmmm_gaussian_fns=added_charges%pgfs, &
401 : para_env=para_env, &
402 : pw_env=pw_env, &
403 : mm_el_pot_radius=added_charges%mm_el_pot_radius, &
404 : mm_el_pot_radius_corr=added_charges%mm_el_pot_radius_corr, &
405 : qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
406 : eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
407 : maxradius=maxradius2, &
408 : maxchrg=maxchrg, &
409 : compatibility=qmmm_env_qm%compatibility, &
410 : print_section=print_section, &
411 8 : qmmm_section=qmmm_section)
412 :
413 16 : SELECT CASE (qmmm_env_qm%qmmm_coupl_type)
414 : CASE (do_qmmm_gauss, do_qmmm_swave, do_qmmm_pcharge)
415 40 : DO i = 1, SIZE(maxradius)
416 40 : maxradius(i) = MAX(maxradius(i), maxradius2(i))
417 : END DO
418 : END SELECT
419 :
420 8 : IF (ASSOCIATED(maxradius2)) DEALLOCATE (maxradius2)
421 : END IF
422 :
423 394 : IF (qmmm_env_qm%added_shells%num_mm_atoms > 0) THEN
424 :
425 60 : maxchrg = MAXVAL(ABS(added_shells%mm_core_chrg(:)))
426 : CALL qmmm_gaussian_initialize(qmmm_gaussian_fns=added_shells%pgfs, &
427 : para_env=para_env, &
428 : pw_env=pw_env, &
429 : mm_el_pot_radius=added_shells%mm_el_pot_radius, &
430 : mm_el_pot_radius_corr=added_shells%mm_el_pot_radius_corr, &
431 : qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
432 : eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
433 : maxradius=maxradius2, &
434 : maxchrg=maxchrg, &
435 : compatibility=qmmm_env_qm%compatibility, &
436 : print_section=print_section, &
437 2 : qmmm_section=qmmm_section)
438 :
439 4 : SELECT CASE (qmmm_env_qm%qmmm_coupl_type)
440 : CASE (do_qmmm_gauss, do_qmmm_swave, do_qmmm_pcharge)
441 10 : DO i = 1, SIZE(maxradius)
442 10 : maxradius(i) = MAX(maxradius(i), maxradius2(i))
443 : END DO
444 : END SELECT
445 :
446 2 : IF (ASSOCIATED(maxradius2)) DEALLOCATE (maxradius2)
447 :
448 : END IF
449 :
450 394 : qmmm_env_qm%maxradius => maxradius
451 :
452 394 : END SUBROUTINE qmmm_init_gaussian_type
453 :
454 : ! **************************************************************************************************
455 : !> \brief ...
456 : !> \param qmmm_env_qm ...
457 : !> \param mm_cell ...
458 : !> \param added_charges ...
459 : !> \param added_shells ...
460 : !> \param print_section ...
461 : !> \par History
462 : !> 1.2005 created [tlaino]
463 : !> \author Teodoro Laino
464 : ! **************************************************************************************************
465 394 : SUBROUTINE qmmm_init_potential(qmmm_env_qm, mm_cell, &
466 : added_charges, added_shells, print_section)
467 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm
468 : TYPE(cell_type), POINTER :: mm_cell
469 : TYPE(add_set_type), POINTER :: added_charges
470 : TYPE(add_shell_type), POINTER :: added_shells
471 : TYPE(section_vals_type), POINTER :: print_section
472 :
473 : CALL qmmm_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
474 : mm_el_pot_radius=qmmm_env_qm%mm_el_pot_radius, &
475 : potentials=qmmm_env_qm%potentials, &
476 : pgfs=qmmm_env_qm%pgfs, &
477 : mm_cell=mm_cell, &
478 : compatibility=qmmm_env_qm%compatibility, &
479 394 : print_section=print_section)
480 :
481 394 : IF (qmmm_env_qm%move_mm_charges .OR. qmmm_env_qm%add_mm_charges) THEN
482 :
483 : CALL qmmm_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
484 : mm_el_pot_radius=added_charges%mm_el_pot_radius, &
485 : potentials=added_charges%potentials, &
486 : pgfs=added_charges%pgfs, &
487 : mm_cell=mm_cell, &
488 : compatibility=qmmm_env_qm%compatibility, &
489 8 : print_section=print_section)
490 : END IF
491 :
492 394 : IF (qmmm_env_qm%added_shells%num_mm_atoms > 0) THEN
493 :
494 : CALL qmmm_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
495 : mm_el_pot_radius=added_shells%mm_el_pot_radius, &
496 : potentials=added_shells%potentials, &
497 : pgfs=added_shells%pgfs, &
498 : mm_cell=mm_cell, &
499 : compatibility=qmmm_env_qm%compatibility, &
500 2 : print_section=print_section)
501 : END IF
502 :
503 394 : END SUBROUTINE qmmm_init_potential
504 :
505 : ! **************************************************************************************************
506 : !> \brief ...
507 : !> \param qmmm_env_qm ...
508 : !> \param qm_cell_small ...
509 : !> \param mm_cell ...
510 : !> \param para_env ...
511 : !> \param qs_env ...
512 : !> \param added_charges ...
513 : !> \param added_shells ...
514 : !> \param qmmm_periodic ...
515 : !> \param print_section ...
516 : !> \param mm_atom_chrg ...
517 : !> \par History
518 : !> 7.2005 created [tlaino]
519 : !> \author Teodoro Laino
520 : ! **************************************************************************************************
521 394 : SUBROUTINE qmmm_init_periodic_potential(qmmm_env_qm, qm_cell_small, mm_cell, para_env, qs_env, &
522 : added_charges, added_shells, qmmm_periodic, print_section, mm_atom_chrg)
523 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm
524 : TYPE(cell_type), POINTER :: qm_cell_small, mm_cell
525 : TYPE(mp_para_env_type), POINTER :: para_env
526 : TYPE(qs_environment_type), POINTER :: qs_env
527 : TYPE(add_set_type), POINTER :: added_charges
528 : TYPE(add_shell_type), POINTER :: added_shells
529 : TYPE(section_vals_type), POINTER :: qmmm_periodic, print_section
530 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_atom_chrg
531 :
532 : REAL(KIND=dp) :: maxchrg
533 : TYPE(dft_control_type), POINTER :: dft_control
534 :
535 394 : IF (qmmm_env_qm%periodic) THEN
536 :
537 46 : NULLIFY (dft_control)
538 46 : CALL get_qs_env(qs_env, dft_control=dft_control)
539 :
540 46 : IF (dft_control%qs_control%semi_empirical) THEN
541 0 : CPABORT("QM/MM periodic calculations not implemented for semi empirical methods")
542 46 : ELSE IF (dft_control%qs_control%dftb) THEN
543 : CALL qmmm_ewald_potential_init(qmmm_env_qm%ewald_env, qmmm_env_qm%ewald_pw, &
544 : qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, mm_cell=mm_cell, &
545 4 : para_env=para_env, qmmm_periodic=qmmm_periodic, print_section=print_section)
546 42 : ELSE IF (dft_control%qs_control%xtb) THEN
547 : CALL qmmm_ewald_potential_init(qmmm_env_qm%ewald_env, qmmm_env_qm%ewald_pw, &
548 : qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, mm_cell=mm_cell, &
549 4 : para_env=para_env, qmmm_periodic=qmmm_periodic, print_section=print_section)
550 : ELSE
551 :
552 : ! setup for GPW/GPAW
553 1560 : maxchrg = MAXVAL(ABS(mm_atom_chrg(:)))
554 38 : IF (qmmm_env_qm%add_mm_charges) maxchrg = MAX(maxchrg, MAXVAL(ABS(added_charges%mm_atom_chrg(:))))
555 :
556 : CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
557 : per_potentials=qmmm_env_qm%per_potentials, &
558 : potentials=qmmm_env_qm%potentials, &
559 : pgfs=qmmm_env_qm%pgfs, &
560 : qm_cell_small=qm_cell_small, &
561 : mm_cell=mm_cell, &
562 : compatibility=qmmm_env_qm%compatibility, &
563 : qmmm_periodic=qmmm_periodic, &
564 : print_section=print_section, &
565 : eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
566 : maxchrg=maxchrg, &
567 : ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, &
568 38 : ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
569 :
570 38 : IF (qmmm_env_qm%move_mm_charges .OR. qmmm_env_qm%add_mm_charges) THEN
571 :
572 : CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
573 : per_potentials=added_charges%per_potentials, &
574 : potentials=added_charges%potentials, &
575 : pgfs=added_charges%pgfs, &
576 : qm_cell_small=qm_cell_small, &
577 : mm_cell=mm_cell, &
578 : compatibility=qmmm_env_qm%compatibility, &
579 : qmmm_periodic=qmmm_periodic, &
580 : print_section=print_section, &
581 : eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
582 : maxchrg=maxchrg, &
583 : ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, &
584 0 : ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
585 : END IF
586 :
587 38 : IF (qmmm_env_qm%added_shells%num_mm_atoms > 0) THEN
588 :
589 : CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
590 : per_potentials=added_shells%per_potentials, &
591 : potentials=added_shells%potentials, &
592 : pgfs=added_shells%pgfs, &
593 : qm_cell_small=qm_cell_small, &
594 : mm_cell=mm_cell, &
595 : compatibility=qmmm_env_qm%compatibility, &
596 : qmmm_periodic=qmmm_periodic, &
597 : print_section=print_section, &
598 : eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
599 : maxchrg=maxchrg, &
600 : ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, &
601 2 : ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
602 : END IF
603 :
604 : END IF
605 :
606 : END IF
607 :
608 394 : END SUBROUTINE qmmm_init_periodic_potential
609 :
610 : ! **************************************************************************************************
611 : !> \brief ...
612 : !> \param qmmm_section ...
613 : !> \param qmmm_env ...
614 : !> \param subsys_mm ...
615 : !> \param qm_atom_type ...
616 : !> \param qm_atom_index ...
617 : !> \param mm_atom_index ...
618 : !> \param qm_cell_small ...
619 : !> \param qmmm_coupl_type ...
620 : !> \param eps_mm_rspace ...
621 : !> \param qmmm_link ...
622 : !> \param para_env ...
623 : !> \par History
624 : !> 11.2004 created [tlaino]
625 : !> \author Teodoro Laino
626 : ! **************************************************************************************************
627 394 : SUBROUTINE setup_qmmm_vars_qm(qmmm_section, qmmm_env, subsys_mm, qm_atom_type, &
628 : qm_atom_index, mm_atom_index, qm_cell_small, qmmm_coupl_type, eps_mm_rspace, &
629 : qmmm_link, para_env)
630 : TYPE(section_vals_type), POINTER :: qmmm_section
631 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
632 : TYPE(cp_subsys_type), POINTER :: subsys_mm
633 : CHARACTER(len=default_string_length), &
634 : DIMENSION(:), POINTER :: qm_atom_type
635 : INTEGER, DIMENSION(:), POINTER :: qm_atom_index, mm_atom_index
636 : TYPE(cell_type), POINTER :: qm_cell_small
637 : INTEGER, INTENT(OUT) :: qmmm_coupl_type
638 : REAL(KIND=dp), INTENT(OUT) :: eps_mm_rspace
639 : LOGICAL, INTENT(OUT) :: qmmm_link
640 : TYPE(mp_para_env_type), POINTER :: para_env
641 :
642 : CHARACTER(len=default_string_length) :: atmname, mm_atom_kind
643 : INTEGER :: i, icount, ikind, ikindr, my_type, &
644 : n_rep_val, nkind, size_mm_system
645 394 : INTEGER, DIMENSION(:), POINTER :: mm_link_atoms
646 : LOGICAL :: explicit, is_mm, is_qm
647 : REAL(KIND=dp) :: tmp_radius, tmp_radius_c
648 394 : REAL(KIND=dp), DIMENSION(:), POINTER :: tmp_sph_cut
649 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
650 : TYPE(atomic_kind_type), POINTER :: atomic_kind
651 : TYPE(fist_potential_type), POINTER :: fist_potential
652 : TYPE(section_vals_type), POINTER :: cell_section, eri_mme_section, &
653 : image_charge_section, mm_kinds
654 :
655 394 : NULLIFY (mm_link_atoms, cell_section, tmp_sph_cut)
656 394 : NULLIFY (image_charge_section)
657 394 : qmmm_link = .FALSE.
658 :
659 394 : CALL section_vals_get(qmmm_section, explicit=explicit)
660 394 : IF (explicit) THEN
661 : !
662 : ! QM_CELL
663 : !
664 394 : cell_section => section_vals_get_subs_vals(qmmm_section, "CELL")
665 : CALL read_cell(qm_cell_small, qm_cell_small, cell_section=cell_section, &
666 394 : check_for_ref=.FALSE., para_env=para_env)
667 394 : qm_cell_small%tag = "CELL_QM"
668 394 : CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type)
669 394 : CALL section_vals_val_get(qmmm_section, "EPS_MM_RSPACE", r_val=eps_mm_rspace)
670 394 : CALL section_vals_val_get(qmmm_section, "SPHERICAL_CUTOFF", r_vals=tmp_sph_cut)
671 394 : CPASSERT(SIZE(tmp_sph_cut) == 2)
672 1970 : qmmm_env%spherical_cutoff = tmp_sph_cut
673 394 : IF (qmmm_env%spherical_cutoff(1) <= 0.0_dp) THEN
674 392 : qmmm_env%spherical_cutoff(2) = 0.0_dp
675 : ELSE
676 2 : IF (qmmm_env%spherical_cutoff(2) <= 0.0_dp) qmmm_env%spherical_cutoff(2) = EPSILON(0.0_dp)
677 2 : tmp_radius = qmmm_env%spherical_cutoff(1) - 20.0_dp*qmmm_env%spherical_cutoff(2)
678 2 : IF (tmp_radius <= 0.0_dp) &
679 : CALL cp_abort(__LOCATION__, &
680 : "SPHERICAL_CUTOFF(1) > 20*SPHERICAL_CUTOFF(1)! Please correct parameters for "// &
681 0 : "the Spherical Cutoff in order to satisfy the previous condition!")
682 : END IF
683 : !
684 : ! Initialization of arrays and core_charge_radius...
685 : !
686 394 : tmp_radius = 0.0_dp
687 394 : CALL cp_subsys_get(subsys=subsys_mm, atomic_kinds=atomic_kinds)
688 4064 : DO Ikind = 1, SIZE(atomic_kinds%els)
689 3670 : atomic_kind => atomic_kinds%els(Ikind)
690 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
691 3670 : fist_potential=fist_potential)
692 : CALL set_potential(potential=fist_potential, &
693 : qmmm_radius=tmp_radius, &
694 3670 : qmmm_corr_radius=tmp_radius)
695 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
696 4064 : fist_potential=fist_potential)
697 : END DO
698 : CALL setup_qm_atom_list(qmmm_section=qmmm_section, &
699 : qm_atom_index=qm_atom_index, &
700 : qm_atom_type=qm_atom_type, &
701 : mm_link_atoms=mm_link_atoms, &
702 394 : qmmm_link=qmmm_link)
703 : !
704 : ! MM_KINDS
705 : !
706 394 : mm_kinds => section_vals_get_subs_vals(qmmm_section, "MM_KIND")
707 394 : CALL section_vals_get(mm_kinds, explicit=explicit, n_repetition=nkind)
708 : !
709 : ! Default
710 : !
711 394 : tmp_radius = cp_unit_to_cp2k(RADIUS_QMMM_DEFAULT, "angstrom")
712 4064 : Set_Radius_Pot_0: DO IkindR = 1, SIZE(atomic_kinds%els)
713 3670 : atomic_kind => atomic_kinds%els(IkindR)
714 3670 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
715 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
716 3670 : fist_potential=fist_potential)
717 : CALL set_potential(potential=fist_potential, qmmm_radius=tmp_radius, &
718 3670 : qmmm_corr_radius=tmp_radius)
719 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
720 4064 : fist_potential=fist_potential)
721 : END DO Set_Radius_Pot_0
722 : !
723 : ! If present overwrite the kind specified in input file...
724 : !
725 394 : IF (explicit) THEN
726 738 : DO ikind = 1, nkind
727 : CALL section_vals_val_get(mm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
728 504 : c_val=mm_atom_kind)
729 504 : CALL section_vals_val_get(mm_kinds, "RADIUS", i_rep_section=ikind, r_val=tmp_radius)
730 504 : tmp_radius_c = tmp_radius
731 504 : CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
732 504 : IF (n_rep_val == 1) CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, &
733 2 : r_val=tmp_radius_c)
734 7254 : Set_Radius_Pot_1: DO IkindR = 1, SIZE(atomic_kinds%els)
735 6012 : atomic_kind => atomic_kinds%els(IkindR)
736 6012 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
737 6012 : is_qm = qmmm_ff_precond_only_qm(atmname)
738 6516 : IF (TRIM(mm_atom_kind) == atmname) THEN
739 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
740 870 : fist_potential=fist_potential)
741 : CALL set_potential(potential=fist_potential, &
742 : qmmm_radius=tmp_radius, &
743 870 : qmmm_corr_radius=tmp_radius_c)
744 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
745 870 : fist_potential=fist_potential)
746 : END IF
747 : END DO Set_Radius_Pot_1
748 : END DO
749 : END IF
750 :
751 : !Image charge section
752 :
753 394 : image_charge_section => section_vals_get_subs_vals(qmmm_section, "IMAGE_CHARGE")
754 394 : CALL section_vals_get(image_charge_section, explicit=qmmm_env%image_charge)
755 :
756 : ELSE
757 0 : CPABORT("QMMM section not present in input file!")
758 : END IF
759 : !
760 : ! Build MM atoms list
761 : !
762 394 : size_mm_system = SIZE(subsys_mm%particles%els) - SIZE(qm_atom_index)
763 394 : IF (qmmm_link .AND. ASSOCIATED(mm_link_atoms)) size_mm_system = size_mm_system + SIZE(mm_link_atoms)
764 1180 : ALLOCATE (mm_atom_index(size_mm_system))
765 394 : icount = 0
766 :
767 190872 : DO i = 1, SIZE(subsys_mm%particles%els)
768 190478 : is_mm = .TRUE.
769 7211572 : IF (ANY(qm_atom_index == i)) THEN
770 2892 : is_mm = .FALSE.
771 : END IF
772 190478 : IF (ASSOCIATED(mm_link_atoms)) THEN
773 790244 : IF (ANY(mm_link_atoms == i) .AND. qmmm_link) is_mm = .TRUE.
774 : END IF
775 190678 : IF (is_mm) THEN
776 187780 : icount = icount + 1
777 187780 : IF (icount <= size_mm_system) mm_atom_index(icount) = i
778 : END IF
779 : END DO
780 394 : CPASSERT(icount == size_mm_system)
781 394 : IF (ASSOCIATED(mm_link_atoms)) THEN
782 62 : DEALLOCATE (mm_link_atoms)
783 : END IF
784 :
785 : ! Build image charge atom list + set up variables
786 : !
787 394 : IF (qmmm_env%image_charge) THEN
788 : CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
789 10 : explicit=explicit)
790 10 : IF (explicit) qmmm_env%image_charge_pot%all_mm = .FALSE.
791 :
792 10 : IF (qmmm_env%image_charge_pot%all_mm) THEN
793 0 : qmmm_env%image_charge_pot%image_mm_list => mm_atom_index
794 : ELSE
795 : CALL setup_image_atom_list(image_charge_section, qmmm_env, &
796 10 : qm_atom_index, subsys_mm)
797 : END IF
798 :
799 10 : qmmm_env%image_charge_pot%particles_all => subsys_mm%particles%els
800 :
801 : CALL section_vals_val_get(image_charge_section, "EXT_POTENTIAL", &
802 10 : r_val=qmmm_env%image_charge_pot%V0)
803 : CALL section_vals_val_get(image_charge_section, "WIDTH", &
804 10 : r_val=qmmm_env%image_charge_pot%eta)
805 : CALL section_vals_val_get(image_charge_section, "DETERM_COEFF", &
806 10 : i_val=my_type)
807 8 : SELECT CASE (my_type)
808 : CASE (do_qmmm_image_calcmatrix)
809 8 : qmmm_env%image_charge_pot%coeff_iterative = .FALSE.
810 : CASE (do_qmmm_image_iter)
811 10 : qmmm_env%image_charge_pot%coeff_iterative = .TRUE.
812 : END SELECT
813 :
814 : CALL section_vals_val_get(image_charge_section, "RESTART_IMAGE_MATRIX", &
815 10 : l_val=qmmm_env%image_charge_pot%image_restart)
816 :
817 : CALL section_vals_val_get(image_charge_section, "IMAGE_MATRIX_METHOD", &
818 10 : i_val=qmmm_env%image_charge_pot%image_matrix_method)
819 :
820 10 : IF (qmmm_env%image_charge_pot%image_matrix_method == do_eri_mme) THEN
821 8 : eri_mme_section => section_vals_get_subs_vals(image_charge_section, "ERI_MME")
822 8 : CALL cp_eri_mme_init_read_input(eri_mme_section, qmmm_env%image_charge_pot%eri_mme_param)
823 : CALL cp_eri_mme_set_params(qmmm_env%image_charge_pot%eri_mme_param, &
824 : hmat=qm_cell_small%hmat, is_ortho=qm_cell_small%orthorhombic, &
825 : zet_min=qmmm_env%image_charge_pot%eta, &
826 : zet_max=qmmm_env%image_charge_pot%eta, &
827 : l_max_zet=0, &
828 : l_max=0, &
829 8 : para_env=para_env)
830 :
831 : END IF
832 : END IF
833 :
834 394 : END SUBROUTINE setup_qmmm_vars_qm
835 :
836 : ! **************************************************************************************************
837 : !> \brief ...
838 : !> \param qmmm_section ...
839 : !> \param qmmm_env ...
840 : !> \param qm_atom_index ...
841 : !> \param mm_link_atoms ...
842 : !> \param mm_link_scale_factor ...
843 : !> \param fist_scale_charge_link ...
844 : !> \param qmmm_coupl_type ...
845 : !> \param qmmm_link ...
846 : !> \par History
847 : !> 12.2004 created [tlaino]
848 : !> \author Teodoro Laino
849 : ! **************************************************************************************************
850 394 : SUBROUTINE setup_qmmm_vars_mm(qmmm_section, qmmm_env, qm_atom_index, &
851 : mm_link_atoms, mm_link_scale_factor, &
852 : fist_scale_charge_link, qmmm_coupl_type, &
853 : qmmm_link)
854 : TYPE(section_vals_type), POINTER :: qmmm_section
855 : TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
856 : INTEGER, DIMENSION(:), POINTER :: qm_atom_index, mm_link_atoms
857 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_link_scale_factor, &
858 : fist_scale_charge_link
859 : INTEGER, INTENT(OUT) :: qmmm_coupl_type
860 : LOGICAL, INTENT(OUT) :: qmmm_link
861 :
862 : LOGICAL :: explicit
863 : TYPE(section_vals_type), POINTER :: qmmm_ff_section
864 :
865 394 : NULLIFY (qmmm_ff_section)
866 394 : qmmm_link = .FALSE.
867 394 : CALL section_vals_get(qmmm_section, explicit=explicit)
868 394 : IF (explicit) THEN
869 394 : CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type)
870 : CALL setup_qm_atom_list(qmmm_section, qm_atom_index=qm_atom_index, qmmm_link=qmmm_link, &
871 : mm_link_atoms=mm_link_atoms, mm_link_scale_factor=mm_link_scale_factor, &
872 394 : fist_scale_charge_link=fist_scale_charge_link)
873 : !
874 : ! Do we want to use a different FF for the non-bonded QM/MM interactions?
875 : !
876 394 : qmmm_ff_section => section_vals_get_subs_vals(qmmm_section, "FORCEFIELD")
877 394 : CALL section_vals_get(qmmm_ff_section, explicit=qmmm_env%use_qmmm_ff)
878 394 : IF (qmmm_env%use_qmmm_ff) THEN
879 : CALL section_vals_val_get(qmmm_ff_section, "MULTIPLE_POTENTIAL", &
880 20 : l_val=qmmm_env%multiple_potential)
881 20 : CALL read_qmmm_ff_section(qmmm_ff_section, qmmm_env%inp_info)
882 : END IF
883 : END IF
884 394 : END SUBROUTINE setup_qmmm_vars_mm
885 :
886 : ! **************************************************************************************************
887 : !> \brief reads information regarding the forcefield specific for the QM/MM
888 : !> interactions
889 : !> \param qmmm_ff_section ...
890 : !> \param inp_info ...
891 : !> \par History
892 : !> 12.2004 created [tlaino]
893 : !> \author Teodoro Laino
894 : ! **************************************************************************************************
895 180 : SUBROUTINE read_qmmm_ff_section(qmmm_ff_section, inp_info)
896 : TYPE(section_vals_type), POINTER :: qmmm_ff_section
897 : TYPE(input_info_type), POINTER :: inp_info
898 :
899 : INTEGER :: n_gd, n_gp, n_lj, n_wl, np
900 : TYPE(section_vals_type), POINTER :: gd_section, gp_section, lj_section, &
901 : wl_section
902 :
903 : !
904 : ! NONBONDED
905 : !
906 :
907 20 : lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%LENNARD-JONES")
908 20 : wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%WILLIAMS")
909 20 : gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GOODWIN")
910 20 : gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GENPOT")
911 20 : CALL section_vals_get(lj_section, n_repetition=n_lj)
912 20 : np = n_lj
913 20 : IF (n_lj /= 0) THEN
914 18 : CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, lj_charmm=.TRUE.)
915 18 : CALL read_lj_section(inp_info%nonbonded, lj_section, start=0)
916 : END IF
917 20 : CALL section_vals_get(wl_section, n_repetition=n_wl)
918 20 : np = n_lj + n_wl
919 20 : IF (n_wl /= 0) THEN
920 2 : CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, williams=.TRUE.)
921 2 : CALL read_wl_section(inp_info%nonbonded, wl_section, start=n_lj)
922 : END IF
923 20 : CALL section_vals_get(gd_section, n_repetition=n_gd)
924 20 : np = n_lj + n_wl + n_gd
925 20 : IF (n_gd /= 0) THEN
926 0 : CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, goodwin=.TRUE.)
927 0 : CALL read_gd_section(inp_info%nonbonded, gd_section, start=n_lj + n_wl)
928 : END IF
929 20 : CALL section_vals_get(gp_section, n_repetition=n_gp)
930 20 : np = n_lj + n_wl + n_gd + n_gp
931 20 : IF (n_gp /= 0) THEN
932 0 : CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, gp=.TRUE.)
933 0 : CALL read_gp_section(inp_info%nonbonded, gp_section, start=n_lj + n_wl + n_gd)
934 : END IF
935 : !
936 : ! NONBONDED14
937 : !
938 20 : lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%LENNARD-JONES")
939 20 : wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%WILLIAMS")
940 20 : gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GOODWIN")
941 20 : gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GENPOT")
942 20 : CALL section_vals_get(lj_section, n_repetition=n_lj)
943 20 : np = n_lj
944 20 : IF (n_lj /= 0) THEN
945 0 : CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, lj_charmm=.TRUE.)
946 0 : CALL read_lj_section(inp_info%nonbonded14, lj_section, start=0)
947 : END IF
948 20 : CALL section_vals_get(wl_section, n_repetition=n_wl)
949 20 : np = n_lj + n_wl
950 20 : IF (n_wl /= 0) THEN
951 0 : CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, williams=.TRUE.)
952 0 : CALL read_wl_section(inp_info%nonbonded14, wl_section, start=n_lj)
953 : END IF
954 20 : CALL section_vals_get(gd_section, n_repetition=n_gd)
955 20 : np = n_lj + n_wl + n_gd
956 20 : IF (n_gd /= 0) THEN
957 0 : CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, goodwin=.TRUE.)
958 0 : CALL read_gd_section(inp_info%nonbonded14, gd_section, start=n_lj + n_wl)
959 : END IF
960 20 : CALL section_vals_get(gp_section, n_repetition=n_gp)
961 20 : np = n_lj + n_wl + n_gd + n_gp
962 20 : IF (n_gp /= 0) THEN
963 0 : CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, gp=.TRUE.)
964 0 : CALL read_gp_section(inp_info%nonbonded14, gp_section, start=n_lj + n_wl + n_gd)
965 : END IF
966 :
967 20 : END SUBROUTINE read_qmmm_ff_section
968 :
969 : ! **************************************************************************************************
970 : !> \brief ...
971 : !> \param qmmm_section ...
972 : !> \param qm_atom_index ...
973 : !> \param qm_atom_type ...
974 : !> \param mm_link_atoms ...
975 : !> \param mm_link_scale_factor ...
976 : !> \param qmmm_link ...
977 : !> \param fist_scale_charge_link ...
978 : !> \par History
979 : !> 12.2004 created [tlaino]
980 : !> \author Teodoro Laino
981 : ! **************************************************************************************************
982 2364 : SUBROUTINE setup_qm_atom_list(qmmm_section, qm_atom_index, qm_atom_type, &
983 : mm_link_atoms, mm_link_scale_factor, qmmm_link, fist_scale_charge_link)
984 : TYPE(section_vals_type), POINTER :: qmmm_section
985 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: qm_atom_index
986 : CHARACTER(len=default_string_length), &
987 : DIMENSION(:), OPTIONAL, POINTER :: qm_atom_type
988 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: mm_link_atoms
989 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: mm_link_scale_factor
990 : LOGICAL, INTENT(OUT), OPTIONAL :: qmmm_link
991 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: fist_scale_charge_link
992 :
993 : CHARACTER(len=default_string_length) :: qm_atom_kind, qm_link_element
994 : INTEGER :: ikind, k, link_involv_mm, link_type, &
995 : mm_index, n_var, nkind, nlinks, &
996 : num_qm_atom_tot
997 788 : INTEGER, DIMENSION(:), POINTER :: mm_indexes
998 : LOGICAL :: explicit
999 : REAL(KIND=dp) :: scale_f
1000 : TYPE(section_vals_type), POINTER :: qm_kinds, qmmm_links
1001 :
1002 788 : num_qm_atom_tot = 0
1003 788 : link_involv_mm = 0
1004 788 : nlinks = 0
1005 : !
1006 : ! QM_KINDS
1007 : !
1008 1576 : qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND")
1009 788 : CALL section_vals_get(qm_kinds, n_repetition=nkind)
1010 2596 : DO ikind = 1, nkind
1011 1808 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
1012 6336 : DO k = 1, n_var
1013 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
1014 3740 : i_vals=mm_indexes)
1015 5548 : num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
1016 : END DO
1017 : END DO
1018 : !
1019 : ! QM/MM LINKS
1020 : !
1021 788 : qmmm_links => section_vals_get_subs_vals(qmmm_section, "LINK")
1022 788 : CALL section_vals_get(qmmm_links, explicit=explicit)
1023 788 : IF (explicit) THEN
1024 128 : qmmm_link = .TRUE.
1025 128 : CALL section_vals_get(qmmm_links, n_repetition=nlinks)
1026 : ! Take care of the various link types
1027 524 : DO ikind = 1, nlinks
1028 : CALL section_vals_val_get(qmmm_links, "LINK_TYPE", i_rep_section=ikind, &
1029 396 : i_val=link_type)
1030 524 : SELECT CASE (link_type)
1031 : CASE (do_qmmm_link_imomm)
1032 388 : num_qm_atom_tot = num_qm_atom_tot + 1
1033 388 : link_involv_mm = link_involv_mm + 1
1034 : CASE (do_qmmm_link_pseudo)
1035 8 : num_qm_atom_tot = num_qm_atom_tot + 1
1036 : CASE (do_qmmm_link_gho)
1037 : ! do nothing for the moment
1038 : CASE DEFAULT
1039 396 : CPABORT("")
1040 : END SELECT
1041 : END DO
1042 : END IF
1043 788 : IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) &
1044 186 : ALLOCATE (mm_link_scale_factor(link_involv_mm))
1045 788 : IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) &
1046 186 : ALLOCATE (fist_scale_charge_link(link_involv_mm))
1047 788 : IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) &
1048 372 : ALLOCATE (mm_link_atoms(link_involv_mm))
1049 2364 : IF (PRESENT(qm_atom_index)) ALLOCATE (qm_atom_index(num_qm_atom_tot))
1050 1576 : IF (PRESENT(qm_atom_type)) ALLOCATE (qm_atom_type(num_qm_atom_tot))
1051 6572 : IF (PRESENT(qm_atom_index)) qm_atom_index = 0
1052 3680 : IF (PRESENT(qm_atom_type)) qm_atom_type = " "
1053 788 : num_qm_atom_tot = 1
1054 2596 : DO ikind = 1, nkind
1055 1808 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
1056 6336 : DO k = 1, n_var
1057 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
1058 3740 : i_vals=mm_indexes)
1059 3740 : IF (PRESENT(qm_atom_index)) THEN
1060 14516 : qm_atom_index(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = mm_indexes(:)
1061 : END IF
1062 3740 : IF (PRESENT(qm_atom_type)) THEN
1063 : CALL section_vals_val_get(qm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
1064 1870 : c_val=qm_atom_kind)
1065 4564 : qm_atom_type(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = qm_atom_kind
1066 : END IF
1067 5548 : num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
1068 : END DO
1069 : END DO
1070 982 : IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) mm_link_scale_factor = 0.0_dp
1071 982 : IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) fist_scale_charge_link = 0.0_dp
1072 1176 : IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) mm_link_atoms = 0
1073 788 : IF (explicit) THEN
1074 524 : DO ikind = 1, nlinks
1075 396 : IF (PRESENT(qm_atom_type)) THEN
1076 198 : CALL section_vals_val_get(qmmm_links, "QM_KIND", i_rep_section=ikind, c_val=qm_link_element)
1077 396 : qm_atom_type(num_qm_atom_tot:num_qm_atom_tot) = TRIM(qm_link_element)//"_LINK"
1078 : END IF
1079 396 : IF (PRESENT(qm_atom_index)) THEN
1080 396 : CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1081 23508 : CPASSERT(ALL(qm_atom_index /= mm_index))
1082 792 : qm_atom_index(num_qm_atom_tot:num_qm_atom_tot) = mm_index
1083 396 : num_qm_atom_tot = num_qm_atom_tot + 1
1084 : END IF
1085 396 : IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) THEN
1086 388 : CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1087 388 : mm_link_atoms(ikind) = mm_index
1088 : END IF
1089 396 : IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) THEN
1090 194 : CALL section_vals_val_get(qmmm_links, "QMMM_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
1091 194 : mm_link_scale_factor(ikind) = scale_f
1092 : END IF
1093 524 : IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) THEN
1094 194 : CALL section_vals_val_get(qmmm_links, "FIST_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
1095 194 : fist_scale_charge_link(ikind) = scale_f
1096 : END IF
1097 : END DO
1098 : END IF
1099 788 : CPASSERT(num_qm_atom_tot - 1 == SIZE(qm_atom_index))
1100 :
1101 788 : END SUBROUTINE setup_qm_atom_list
1102 :
1103 : ! **************************************************************************************************
1104 : !> \brief this routine sets up all variables to treat qmmm links
1105 : !> \param qmmm_section ...
1106 : !> \param qmmm_links ...
1107 : !> \param mm_el_pot_radius ...
1108 : !> \param mm_el_pot_radius_corr ...
1109 : !> \param mm_atom_index ...
1110 : !> \param iw ...
1111 : !> \par History
1112 : !> 12.2004 created [tlaino]
1113 : !> \author Teodoro Laino
1114 : ! **************************************************************************************************
1115 128 : SUBROUTINE setup_qmmm_links(qmmm_section, qmmm_links, mm_el_pot_radius, mm_el_pot_radius_corr, &
1116 : mm_atom_index, iw)
1117 : TYPE(section_vals_type), POINTER :: qmmm_section
1118 : TYPE(qmmm_links_type), POINTER :: qmmm_links
1119 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_el_pot_radius, mm_el_pot_radius_corr
1120 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1121 : INTEGER, INTENT(IN) :: iw
1122 :
1123 : INTEGER :: ikind, link_type, mm_index, n_gho, &
1124 : n_imomm, n_pseudo, n_rep_val, n_tot, &
1125 : nlinks, qm_index
1126 64 : INTEGER, DIMENSION(:), POINTER :: wrk_tmp
1127 : REAL(KIND=dp) :: alpha, my_radius
1128 : TYPE(section_vals_type), POINTER :: qmmm_link_section
1129 :
1130 64 : NULLIFY (wrk_tmp)
1131 64 : n_imomm = 0
1132 64 : n_gho = 0
1133 64 : n_pseudo = 0
1134 128 : qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
1135 64 : CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
1136 64 : CPASSERT(nlinks /= 0)
1137 262 : DO ikind = 1, nlinks
1138 198 : CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1139 198 : IF (link_type == do_qmmm_link_imomm) n_imomm = n_imomm + 1
1140 198 : IF (link_type == do_qmmm_link_gho) n_gho = n_gho + 1
1141 460 : IF (link_type == do_qmmm_link_pseudo) n_pseudo = n_pseudo + 1
1142 : END DO
1143 64 : n_tot = n_imomm + n_gho + n_pseudo
1144 64 : CPASSERT(n_tot /= 0)
1145 64 : ALLOCATE (qmmm_links)
1146 : NULLIFY (qmmm_links%imomm, &
1147 : qmmm_links%pseudo)
1148 : ! IMOMM
1149 64 : IF (n_imomm /= 0) THEN
1150 380 : ALLOCATE (qmmm_links%imomm(n_imomm))
1151 186 : ALLOCATE (wrk_tmp(n_imomm))
1152 256 : DO ikind = 1, n_imomm
1153 194 : NULLIFY (qmmm_links%imomm(ikind)%link)
1154 256 : ALLOCATE (qmmm_links%imomm(ikind)%link)
1155 : END DO
1156 62 : n_imomm = 0
1157 256 : DO ikind = 1, nlinks
1158 194 : CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1159 256 : IF (link_type == do_qmmm_link_imomm) THEN
1160 194 : n_imomm = n_imomm + 1
1161 194 : CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
1162 194 : CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1163 194 : CALL section_vals_val_get(qmmm_link_section, "ALPHA_IMOMM", i_rep_section=ikind, r_val=alpha)
1164 194 : CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
1165 194 : qmmm_links%imomm(n_imomm)%link%qm_index = qm_index
1166 194 : qmmm_links%imomm(n_imomm)%link%mm_index = mm_index
1167 194 : qmmm_links%imomm(n_imomm)%link%alpha = alpha
1168 194 : wrk_tmp(n_imomm) = mm_index
1169 194 : IF (n_rep_val == 1) THEN
1170 64 : CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, r_val=my_radius)
1171 996 : WHERE (mm_atom_index == mm_index) mm_el_pot_radius = my_radius
1172 996 : WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
1173 : END IF
1174 194 : CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
1175 194 : IF (n_rep_val == 1) THEN
1176 0 : CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, r_val=my_radius)
1177 0 : WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
1178 : END IF
1179 : END IF
1180 : END DO
1181 : !
1182 : ! Checking the link structure
1183 : !
1184 256 : DO ikind = 1, SIZE(wrk_tmp)
1185 3514 : IF (COUNT(wrk_tmp == wrk_tmp(ikind)) > 1) THEN
1186 0 : WRITE (iw, '(/A)') "In the IMOMM scheme no more than one QM atom can be bounded to the same MM atom."
1187 0 : WRITE (iw, '(A)') "Multiple link MM atom not allowed. Check your link sections."
1188 0 : CPABORT("")
1189 : END IF
1190 : END DO
1191 62 : DEALLOCATE (wrk_tmp)
1192 : END IF
1193 : ! PSEUDO
1194 64 : IF (n_pseudo /= 0) THEN
1195 10 : ALLOCATE (qmmm_links%pseudo(n_pseudo))
1196 6 : DO ikind = 1, n_pseudo
1197 4 : NULLIFY (qmmm_links%pseudo(ikind)%link)
1198 6 : ALLOCATE (qmmm_links%pseudo(ikind)%link)
1199 : END DO
1200 2 : n_pseudo = 0
1201 6 : DO ikind = 1, nlinks
1202 4 : CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1203 6 : IF (link_type == do_qmmm_link_pseudo) THEN
1204 4 : n_pseudo = n_pseudo + 1
1205 4 : CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
1206 4 : CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1207 4 : qmmm_links%pseudo(n_pseudo)%link%qm_index = qm_index
1208 4 : qmmm_links%pseudo(n_pseudo)%link%mm_index = mm_index
1209 : END IF
1210 : END DO
1211 : END IF
1212 : ! GHO
1213 64 : IF (n_gho /= 0) THEN
1214 : ! not yet implemented
1215 : ! still to define : type, implementation into QS
1216 0 : CPABORT("")
1217 : END IF
1218 64 : END SUBROUTINE setup_qmmm_links
1219 :
1220 : ! **************************************************************************************************
1221 : !> \brief this routine sets up all variables to treat qmmm links
1222 : !> \param qmmm_section ...
1223 : !> \param move_mm_charges ...
1224 : !> \param add_mm_charges ...
1225 : !> \param mm_atom_chrg ...
1226 : !> \param mm_el_pot_radius ...
1227 : !> \param mm_el_pot_radius_corr ...
1228 : !> \param added_charges ...
1229 : !> \param mm_atom_index ...
1230 : !> \par History
1231 : !> 12.2004 created [tlaino]
1232 : !> \author Teodoro Laino
1233 : ! **************************************************************************************************
1234 128 : SUBROUTINE move_or_add_atoms(qmmm_section, move_mm_charges, add_mm_charges, &
1235 : mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, &
1236 : added_charges, mm_atom_index)
1237 : TYPE(section_vals_type), POINTER :: qmmm_section
1238 : LOGICAL, INTENT(OUT) :: move_mm_charges, add_mm_charges
1239 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, &
1240 : mm_el_pot_radius_corr
1241 : TYPE(add_set_type), POINTER :: added_charges
1242 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1243 :
1244 : INTEGER :: i_add, icount, ikind, ind1, Index1, &
1245 : Index2, n_add_tot, n_adds, n_move_tot, &
1246 : n_moves, n_rep_val, nlinks
1247 : LOGICAL :: explicit
1248 : REAL(KIND=dp) :: alpha, c_radius, charge, radius
1249 : TYPE(section_vals_type), POINTER :: add_section, move_section, &
1250 : qmmm_link_section
1251 :
1252 64 : explicit = .FALSE.
1253 64 : move_mm_charges = .FALSE.
1254 64 : add_mm_charges = .FALSE.
1255 64 : NULLIFY (qmmm_link_section, move_section, add_section)
1256 64 : qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
1257 64 : CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
1258 64 : CPASSERT(nlinks /= 0)
1259 : icount = 0
1260 64 : n_move_tot = 0
1261 64 : n_add_tot = 0
1262 262 : DO ikind = 1, nlinks
1263 : move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
1264 198 : i_rep_section=ikind)
1265 198 : CALL section_vals_get(move_section, n_repetition=n_moves)
1266 : add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", &
1267 198 : i_rep_section=ikind)
1268 198 : CALL section_vals_get(add_section, n_repetition=n_adds)
1269 198 : n_move_tot = n_move_tot + n_moves
1270 460 : n_add_tot = n_add_tot + n_adds
1271 : END DO
1272 64 : icount = n_move_tot + n_add_tot
1273 64 : IF (n_add_tot /= 0) add_mm_charges = .TRUE.
1274 64 : IF (n_move_tot /= 0) move_mm_charges = .TRUE.
1275 : !
1276 : ! create add_set_type
1277 : !
1278 64 : CALL create_add_set_type(added_charges, ndim=icount)
1279 : !
1280 : ! Fill in structures
1281 : !
1282 64 : icount = 0
1283 262 : DO ikind = 1, nlinks
1284 : move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
1285 198 : i_rep_section=ikind)
1286 198 : CALL section_vals_get(move_section, explicit=explicit, n_repetition=n_moves)
1287 : !
1288 : ! Moving charge atoms
1289 : !
1290 198 : IF (explicit) THEN
1291 36 : DO i_add = 1, n_moves
1292 26 : icount = icount + 1
1293 26 : CALL section_vals_val_get(move_section, "ATOM_INDEX_1", i_val=Index1, i_rep_section=i_add)
1294 26 : CALL section_vals_val_get(move_section, "ATOM_INDEX_2", i_val=Index2, i_rep_section=i_add)
1295 26 : CALL section_vals_val_get(move_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
1296 26 : CALL section_vals_val_get(move_section, "RADIUS", r_val=radius, i_rep_section=i_add)
1297 26 : CALL section_vals_val_get(move_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
1298 26 : c_radius = radius
1299 26 : IF (n_rep_val == 1) &
1300 26 : CALL section_vals_val_get(move_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add)
1301 :
1302 : CALL set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, &
1303 : mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, &
1304 : mm_el_pot_radius_corr=mm_el_pot_radius_corr, &
1305 62 : mm_atom_index=mm_atom_index, move=n_moves, Ind1=ind1)
1306 : END DO
1307 10 : mm_atom_chrg(ind1) = 0.0_dp
1308 : END IF
1309 :
1310 : add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", &
1311 198 : i_rep_section=ikind)
1312 198 : CALL section_vals_get(add_section, explicit=explicit, n_repetition=n_adds)
1313 : !
1314 : ! Adding charge atoms
1315 : !
1316 460 : IF (explicit) THEN
1317 4 : DO i_add = 1, n_adds
1318 2 : icount = icount + 1
1319 2 : CALL section_vals_val_get(add_section, "ATOM_INDEX_1", i_val=Index1, i_rep_section=i_add)
1320 2 : CALL section_vals_val_get(add_section, "ATOM_INDEX_2", i_val=Index2, i_rep_section=i_add)
1321 2 : CALL section_vals_val_get(add_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
1322 2 : CALL section_vals_val_get(add_section, "RADIUS", r_val=radius, i_rep_section=i_add)
1323 2 : CALL section_vals_val_get(add_section, "CHARGE", r_val=charge, i_rep_section=i_add)
1324 2 : CALL section_vals_val_get(add_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
1325 2 : c_radius = radius
1326 2 : IF (n_rep_val == 1) &
1327 2 : CALL section_vals_val_get(add_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add)
1328 :
1329 : CALL set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, charge, &
1330 : mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, &
1331 : mm_el_pot_radius_corr=mm_el_pot_radius_corr, &
1332 6 : mm_atom_index=mm_atom_index)
1333 : END DO
1334 : END IF
1335 : END DO
1336 :
1337 64 : END SUBROUTINE move_or_add_atoms
1338 :
1339 : ! **************************************************************************************************
1340 : !> \brief this routine sets up all variables of the add_set_type type
1341 : !> \param added_charges ...
1342 : !> \param icount ...
1343 : !> \param Index1 ...
1344 : !> \param Index2 ...
1345 : !> \param alpha ...
1346 : !> \param radius ...
1347 : !> \param c_radius ...
1348 : !> \param charge ...
1349 : !> \param mm_atom_chrg ...
1350 : !> \param mm_el_pot_radius ...
1351 : !> \param mm_el_pot_radius_corr ...
1352 : !> \param mm_atom_index ...
1353 : !> \param move ...
1354 : !> \param ind1 ...
1355 : !> \par History
1356 : !> 12.2004 created [tlaino]
1357 : !> \author Teodoro Laino
1358 : ! **************************************************************************************************
1359 28 : SUBROUTINE set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, charge, &
1360 : mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, mm_atom_index, move, ind1)
1361 : TYPE(add_set_type), POINTER :: added_charges
1362 : INTEGER, INTENT(IN) :: icount, Index1, Index2
1363 : REAL(KIND=dp), INTENT(IN) :: alpha, radius, c_radius
1364 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: charge
1365 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, &
1366 : mm_el_pot_radius_corr
1367 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1368 : INTEGER, INTENT(in), OPTIONAL :: move
1369 : INTEGER, INTENT(OUT), OPTIONAL :: ind1
1370 :
1371 : INTEGER :: i, my_move
1372 : REAL(KIND=dp) :: my_c_radius, my_charge, my_radius
1373 :
1374 28 : my_move = 0
1375 28 : my_radius = radius
1376 28 : my_c_radius = c_radius
1377 28 : IF (PRESENT(charge)) my_charge = charge
1378 28 : IF (PRESENT(move)) my_move = move
1379 28 : i = 1
1380 60 : GetId: DO WHILE (i <= SIZE(mm_atom_index))
1381 60 : IF (Index1 == mm_atom_index(i)) EXIT GetId
1382 60 : i = i + 1
1383 : END DO GetId
1384 28 : IF (PRESENT(ind1)) ind1 = i
1385 28 : CPASSERT(i <= SIZE(mm_atom_index))
1386 28 : IF (.NOT. PRESENT(charge)) my_charge = mm_atom_chrg(i)/REAL(my_move, KIND=dp)
1387 28 : IF (my_radius == 0.0_dp) my_radius = mm_el_pot_radius(i)
1388 28 : IF (my_c_radius == 0.0_dp) my_c_radius = mm_el_pot_radius_corr(i)
1389 :
1390 28 : added_charges%add_env(icount)%Index1 = Index1
1391 28 : added_charges%add_env(icount)%Index2 = Index2
1392 28 : added_charges%add_env(icount)%alpha = alpha
1393 28 : added_charges%mm_atom_index(icount) = icount
1394 28 : added_charges%mm_atom_chrg(icount) = my_charge
1395 28 : added_charges%mm_el_pot_radius(icount) = my_radius
1396 28 : added_charges%mm_el_pot_radius_corr(icount) = my_c_radius
1397 28 : END SUBROUTINE set_add_set_type
1398 :
1399 : ! **************************************************************************************************
1400 : !> \brief this routine sets up the origin of the MM cell respect to the
1401 : !> origin of the QM cell. The origin of the QM cell is assumed to be
1402 : !> in (0.0,0.0,0.0)...
1403 : !> \param qmmm_section ...
1404 : !> \param qmmm_env ...
1405 : !> \param qm_cell_small ...
1406 : !> \param dr ...
1407 : !> \par History
1408 : !> 02.2005 created [tlaino]
1409 : !> \author Teodoro Laino
1410 : ! **************************************************************************************************
1411 788 : SUBROUTINE setup_origin_mm_cell(qmmm_section, qmmm_env, qm_cell_small, &
1412 : dr)
1413 : TYPE(section_vals_type), POINTER :: qmmm_section
1414 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1415 : TYPE(cell_type), POINTER :: qm_cell_small
1416 : REAL(KIND=dp), DIMENSION(3), INTENT(in) :: dr
1417 :
1418 : LOGICAL :: center_grid
1419 : REAL(KIND=dp), DIMENSION(3) :: tmp
1420 394 : REAL(KINd=dp), DIMENSION(:), POINTER :: vec
1421 :
1422 : ! This is the vector that corrects position to apply properly the PBC
1423 :
1424 394 : tmp(1) = qm_cell_small%hmat(1, 1)
1425 394 : tmp(2) = qm_cell_small%hmat(2, 2)
1426 394 : tmp(3) = qm_cell_small%hmat(3, 3)
1427 1576 : CPASSERT(ALL(tmp > 0))
1428 1576 : qmmm_env%dOmmOqm = tmp/2.0_dp
1429 : ! This is unit vector to translate the QM system in order to center it
1430 : ! in QM cell
1431 394 : CALL section_vals_val_get(qmmm_section, "CENTER_GRID", l_val=center_grid)
1432 394 : IF (center_grid) THEN
1433 72 : qmmm_env%utrasl = dr
1434 : ELSE
1435 1504 : qmmm_env%utrasl = 1.0_dp
1436 : END IF
1437 394 : CALL section_vals_val_get(qmmm_section, "INITIAL_TRANSLATION_VECTOR", r_vals=vec)
1438 3152 : qmmm_env%transl_v = vec
1439 394 : END SUBROUTINE setup_origin_mm_cell
1440 :
1441 : ! **************************************************************************************************
1442 : !> \brief this routine sets up list of MM atoms carrying an image charge
1443 : !> \param image_charge_section ...
1444 : !> \param qmmm_env ...
1445 : !> \param qm_atom_index ...
1446 : !> \param subsys_mm ...
1447 : !> \par History
1448 : !> 02.2012 created
1449 : !> \author Dorothea Golze
1450 : ! **************************************************************************************************
1451 10 : SUBROUTINE setup_image_atom_list(image_charge_section, qmmm_env, &
1452 : qm_atom_index, subsys_mm)
1453 :
1454 : TYPE(section_vals_type), POINTER :: image_charge_section
1455 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1456 : INTEGER, DIMENSION(:), POINTER :: qm_atom_index
1457 : TYPE(cp_subsys_type), POINTER :: subsys_mm
1458 :
1459 : INTEGER :: atom_a, atom_b, i, j, k, max_index, &
1460 : n_var, num_const_atom, &
1461 : num_image_mm_atom
1462 10 : INTEGER, DIMENSION(:), POINTER :: mm_indexes
1463 : LOGICAL :: fix_xyz, imageind_in_range
1464 10 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind
1465 :
1466 10 : NULLIFY (mm_indexes, molecule_kind)
1467 10 : imageind_in_range = .FALSE.
1468 10 : num_image_mm_atom = 0
1469 10 : max_index = 0
1470 :
1471 : CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1472 10 : n_rep_val=n_var)
1473 20 : DO i = 1, n_var
1474 : CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1475 10 : i_rep_val=i, i_vals=mm_indexes)
1476 20 : num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
1477 : END DO
1478 :
1479 30 : ALLOCATE (qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom))
1480 :
1481 30 : qmmm_env%image_charge_pot%image_mm_list = 0
1482 10 : num_image_mm_atom = 1
1483 :
1484 20 : DO i = 1, n_var
1485 : CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1486 10 : i_rep_val=i, i_vals=mm_indexes)
1487 : qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom:num_image_mm_atom &
1488 60 : + SIZE(mm_indexes) - 1) = mm_indexes(:)
1489 20 : num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
1490 : END DO
1491 :
1492 : ! checking, if in range, if list contains QM atoms or any atoms doubled
1493 10 : num_image_mm_atom = num_image_mm_atom - 1
1494 :
1495 10 : max_index = SIZE(subsys_mm%particles%els)
1496 :
1497 10 : CPASSERT(SIZE(qmmm_env%image_charge_pot%image_mm_list) /= 0)
1498 : imageind_in_range = (MAXVAL(qmmm_env%image_charge_pot%image_mm_list) <= max_index) &
1499 60 : .AND. (MINVAL(qmmm_env%image_charge_pot%image_mm_list) > 0)
1500 10 : CPASSERT(imageind_in_range)
1501 :
1502 30 : DO i = 1, num_image_mm_atom
1503 20 : atom_a = qmmm_env%image_charge_pot%image_mm_list(i)
1504 60 : IF (ANY(qm_atom_index == atom_a)) THEN
1505 0 : CPABORT("Image atom list must only contain MM atoms")
1506 : END IF
1507 40 : DO j = i + 1, num_image_mm_atom
1508 10 : atom_b = qmmm_env%image_charge_pot%image_mm_list(j)
1509 10 : IF (atom_a == atom_b) &
1510 20 : CPABORT("There are atoms doubled in image list.")
1511 : END DO
1512 : END DO
1513 :
1514 : ! check if molecules in list carry constraints
1515 10 : num_const_atom = 0
1516 10 : fix_xyz = .TRUE.
1517 10 : IF (ASSOCIATED(subsys_mm%molecule_kinds)) THEN
1518 10 : IF (ASSOCIATED(subsys_mm%molecule_kinds%els)) THEN
1519 10 : molecule_kind => subsys_mm%molecule_kinds%els
1520 78 : DO i = 1, SIZE(molecule_kind)
1521 76 : IF (.NOT. ASSOCIATED(molecule_kind(i)%fixd_list)) EXIT
1522 68 : IF (.NOT. fix_xyz) EXIT
1523 82 : DO j = 1, SIZE(molecule_kind(i)%fixd_list)
1524 4 : IF (.NOT. fix_xyz) EXIT
1525 80 : DO k = 1, num_image_mm_atom
1526 8 : atom_a = qmmm_env%image_charge_pot%image_mm_list(k)
1527 12 : IF (atom_a == molecule_kind(i)%fixd_list(j)%fixd) THEN
1528 4 : num_const_atom = num_const_atom + 1
1529 4 : IF (molecule_kind(i)%fixd_list(j)%itype /= use_perd_xyz) THEN
1530 : fix_xyz = .FALSE.
1531 : EXIT
1532 : END IF
1533 : END IF
1534 : END DO
1535 : END DO
1536 : END DO
1537 : END IF
1538 : END IF
1539 :
1540 : ! if all image atoms are constrained, calculate image matrix only
1541 : ! once for the first MD or GEO_OPT step (for non-iterative case)
1542 10 : IF (num_const_atom == num_image_mm_atom .AND. fix_xyz) THEN
1543 2 : qmmm_env%image_charge_pot%state_image_matrix = calc_once
1544 : ELSE
1545 8 : qmmm_env%image_charge_pot%state_image_matrix = calc_always
1546 : END IF
1547 :
1548 10 : END SUBROUTINE setup_image_atom_list
1549 :
1550 : ! **************************************************************************************************
1551 : !> \brief Print info on image charges
1552 : !> \param qmmm_env ...
1553 : !> \param qmmm_section ...
1554 : !> \par History
1555 : !> 03.2012 created
1556 : !> \author Dorothea Golze
1557 : ! **************************************************************************************************
1558 10 : SUBROUTINE print_image_charge_info(qmmm_env, qmmm_section)
1559 :
1560 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1561 : TYPE(section_vals_type), POINTER :: qmmm_section
1562 :
1563 : INTEGER :: iw
1564 : REAL(KIND=dp) :: eta, eta_conv, V0, V0_conv
1565 : TYPE(cp_logger_type), POINTER :: logger
1566 :
1567 10 : logger => cp_get_default_logger()
1568 : iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%PROGRAM_RUN_INFO", &
1569 10 : extension=".log")
1570 10 : eta = qmmm_env%image_charge_pot%eta
1571 10 : eta_conv = cp_unit_from_cp2k(eta, "angstrom", power=-2)
1572 10 : V0 = qmmm_env%image_charge_pot%V0
1573 10 : V0_conv = cp_unit_from_cp2k(V0, "volt")
1574 :
1575 10 : IF (iw > 0) THEN
1576 5 : WRITE (iw, FMT="(T25,A)") "IMAGE CHARGE PARAMETERS"
1577 5 : WRITE (iw, FMT="(T25,A)") REPEAT("-", 23)
1578 5 : WRITE (iw, FMT="(/)")
1579 5 : WRITE (iw, FMT="(T2,A)") "INDEX OF MM ATOMS CARRYING AN IMAGE CHARGE:"
1580 5 : WRITE (iw, FMT="(/)")
1581 :
1582 15 : WRITE (iw, "(7X,10I6)") qmmm_env%image_charge_pot%image_mm_list
1583 5 : WRITE (iw, FMT="(/)")
1584 : WRITE (iw, "(T2,A52,T69,F12.8)") &
1585 5 : "WIDTH OF GAUSSIAN CHARGE DISTRIBUTION [angstrom^-2]:", eta_conv
1586 5 : WRITE (iw, "(T2,A26,T69,F12.8)") "EXTERNAL POTENTIAL [volt]:", V0_conv
1587 5 : WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 79)
1588 : END IF
1589 : CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
1590 10 : "PRINT%PROGRAM_RUN_INFO")
1591 :
1592 10 : END SUBROUTINE print_image_charge_info
1593 :
1594 : END MODULE qmmm_init
|
