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