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 188568 : maxchrg = MAXVAL(ABS(mm_atom_chrg(:)))
383 394 : CALL get_qs_env(qs_env, pw_env=pw_env)
384 410 : 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 60 : 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 1560 : 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) &
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 : !
676 : ! Initialization of arrays and core_charge_radius...
677 : !
678 394 : tmp_radius = 0.0_dp
679 394 : CALL cp_subsys_get(subsys=subsys_mm, atomic_kinds=atomic_kinds)
680 4064 : DO Ikind = 1, SIZE(atomic_kinds%els)
681 3670 : atomic_kind => atomic_kinds%els(Ikind)
682 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
683 3670 : fist_potential=fist_potential)
684 : CALL set_potential(potential=fist_potential, &
685 : qmmm_radius=tmp_radius, &
686 3670 : qmmm_corr_radius=tmp_radius)
687 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
688 4064 : fist_potential=fist_potential)
689 : END DO
690 : CALL setup_qm_atom_list(qmmm_section=qmmm_section, &
691 : qm_atom_index=qm_atom_index, &
692 : qm_atom_type=qm_atom_type, &
693 : mm_link_atoms=mm_link_atoms, &
694 394 : qmmm_link=qmmm_link)
695 : !
696 : ! MM_KINDS
697 : !
698 394 : mm_kinds => section_vals_get_subs_vals(qmmm_section, "MM_KIND")
699 394 : CALL section_vals_get(mm_kinds, explicit=explicit, n_repetition=nkind)
700 : !
701 : ! Default
702 : !
703 394 : tmp_radius = cp_unit_to_cp2k(RADIUS_QMMM_DEFAULT, "angstrom")
704 4064 : Set_Radius_Pot_0: DO IkindR = 1, SIZE(atomic_kinds%els)
705 3670 : atomic_kind => atomic_kinds%els(IkindR)
706 3670 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
707 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
708 3670 : fist_potential=fist_potential)
709 : CALL set_potential(potential=fist_potential, qmmm_radius=tmp_radius, &
710 3670 : qmmm_corr_radius=tmp_radius)
711 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
712 4064 : fist_potential=fist_potential)
713 : END DO Set_Radius_Pot_0
714 : !
715 : ! If present overwrite the kind specified in input file...
716 : !
717 394 : IF (explicit) THEN
718 738 : DO ikind = 1, nkind
719 : CALL section_vals_val_get(mm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
720 504 : c_val=mm_atom_kind)
721 504 : CALL section_vals_val_get(mm_kinds, "RADIUS", i_rep_section=ikind, r_val=tmp_radius)
722 504 : tmp_radius_c = tmp_radius
723 504 : CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
724 504 : IF (n_rep_val == 1) CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, &
725 2 : r_val=tmp_radius_c)
726 7254 : Set_Radius_Pot_1: DO IkindR = 1, SIZE(atomic_kinds%els)
727 6012 : atomic_kind => atomic_kinds%els(IkindR)
728 6012 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
729 6012 : is_qm = qmmm_ff_precond_only_qm(atmname)
730 6516 : IF (TRIM(mm_atom_kind) == atmname) THEN
731 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
732 870 : fist_potential=fist_potential)
733 : CALL set_potential(potential=fist_potential, &
734 : qmmm_radius=tmp_radius, &
735 870 : qmmm_corr_radius=tmp_radius_c)
736 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
737 870 : fist_potential=fist_potential)
738 : END IF
739 : END DO Set_Radius_Pot_1
740 : END DO
741 : END IF
742 :
743 : !Image charge section
744 :
745 394 : image_charge_section => section_vals_get_subs_vals(qmmm_section, "IMAGE_CHARGE")
746 394 : CALL section_vals_get(image_charge_section, explicit=qmmm_env%image_charge)
747 :
748 : ELSE
749 0 : CPABORT("QMMM section not present in input file!")
750 : END IF
751 : !
752 : ! Build MM atoms list
753 : !
754 394 : size_mm_system = SIZE(subsys_mm%particles%els) - SIZE(qm_atom_index)
755 394 : IF (qmmm_link .AND. ASSOCIATED(mm_link_atoms)) size_mm_system = size_mm_system + SIZE(mm_link_atoms)
756 1180 : ALLOCATE (mm_atom_index(size_mm_system))
757 394 : icount = 0
758 :
759 190872 : DO i = 1, SIZE(subsys_mm%particles%els)
760 190478 : is_mm = .TRUE.
761 7211572 : IF (ANY(qm_atom_index == i)) THEN
762 2892 : is_mm = .FALSE.
763 : END IF
764 190478 : IF (ASSOCIATED(mm_link_atoms)) THEN
765 790244 : IF (ANY(mm_link_atoms == i) .AND. qmmm_link) is_mm = .TRUE.
766 : END IF
767 190678 : IF (is_mm) THEN
768 187780 : icount = icount + 1
769 187780 : IF (icount <= size_mm_system) mm_atom_index(icount) = i
770 : END IF
771 : END DO
772 394 : CPASSERT(icount == size_mm_system)
773 394 : IF (ASSOCIATED(mm_link_atoms)) THEN
774 62 : DEALLOCATE (mm_link_atoms)
775 : END IF
776 :
777 : ! Build image charge atom list + set up variables
778 : !
779 394 : IF (qmmm_env%image_charge) THEN
780 : CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
781 10 : explicit=explicit)
782 10 : IF (explicit) qmmm_env%image_charge_pot%all_mm = .FALSE.
783 :
784 10 : IF (qmmm_env%image_charge_pot%all_mm) THEN
785 0 : qmmm_env%image_charge_pot%image_mm_list => mm_atom_index
786 : ELSE
787 : CALL setup_image_atom_list(image_charge_section, qmmm_env, &
788 10 : qm_atom_index, subsys_mm)
789 : END IF
790 :
791 10 : qmmm_env%image_charge_pot%particles_all => subsys_mm%particles%els
792 :
793 : CALL section_vals_val_get(image_charge_section, "EXT_POTENTIAL", &
794 10 : r_val=qmmm_env%image_charge_pot%V0)
795 : CALL section_vals_val_get(image_charge_section, "WIDTH", &
796 10 : r_val=qmmm_env%image_charge_pot%eta)
797 : CALL section_vals_val_get(image_charge_section, "DETERM_COEFF", &
798 10 : i_val=my_type)
799 8 : SELECT CASE (my_type)
800 : CASE (do_qmmm_image_calcmatrix)
801 8 : qmmm_env%image_charge_pot%coeff_iterative = .FALSE.
802 : CASE (do_qmmm_image_iter)
803 10 : qmmm_env%image_charge_pot%coeff_iterative = .TRUE.
804 : END SELECT
805 :
806 : CALL section_vals_val_get(image_charge_section, "RESTART_IMAGE_MATRIX", &
807 10 : l_val=qmmm_env%image_charge_pot%image_restart)
808 :
809 : CALL section_vals_val_get(image_charge_section, "IMAGE_MATRIX_METHOD", &
810 10 : i_val=qmmm_env%image_charge_pot%image_matrix_method)
811 :
812 10 : IF (qmmm_env%image_charge_pot%image_matrix_method == do_eri_mme) THEN
813 8 : eri_mme_section => section_vals_get_subs_vals(image_charge_section, "ERI_MME")
814 8 : CALL cp_eri_mme_init_read_input(eri_mme_section, qmmm_env%image_charge_pot%eri_mme_param)
815 : CALL cp_eri_mme_set_params(qmmm_env%image_charge_pot%eri_mme_param, &
816 : hmat=qm_cell_small%hmat, is_ortho=qm_cell_small%orthorhombic, &
817 : zet_min=qmmm_env%image_charge_pot%eta, &
818 : zet_max=qmmm_env%image_charge_pot%eta, &
819 : l_max_zet=0, &
820 : l_max=0, &
821 8 : para_env=para_env)
822 :
823 : END IF
824 : END IF
825 :
826 394 : END SUBROUTINE setup_qmmm_vars_qm
827 :
828 : ! **************************************************************************************************
829 : !> \brief ...
830 : !> \param qmmm_section ...
831 : !> \param qmmm_env ...
832 : !> \param qm_atom_index ...
833 : !> \param mm_link_atoms ...
834 : !> \param mm_link_scale_factor ...
835 : !> \param fist_scale_charge_link ...
836 : !> \param qmmm_coupl_type ...
837 : !> \param qmmm_link ...
838 : !> \par History
839 : !> 12.2004 created [tlaino]
840 : !> \author Teodoro Laino
841 : ! **************************************************************************************************
842 394 : SUBROUTINE setup_qmmm_vars_mm(qmmm_section, qmmm_env, qm_atom_index, &
843 : mm_link_atoms, mm_link_scale_factor, &
844 : fist_scale_charge_link, qmmm_coupl_type, &
845 : qmmm_link)
846 : TYPE(section_vals_type), POINTER :: qmmm_section
847 : TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
848 : INTEGER, DIMENSION(:), POINTER :: qm_atom_index, mm_link_atoms
849 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_link_scale_factor, &
850 : fist_scale_charge_link
851 : INTEGER, INTENT(OUT) :: qmmm_coupl_type
852 : LOGICAL, INTENT(OUT) :: qmmm_link
853 :
854 : LOGICAL :: explicit
855 : TYPE(section_vals_type), POINTER :: qmmm_ff_section
856 :
857 394 : NULLIFY (qmmm_ff_section)
858 394 : qmmm_link = .FALSE.
859 394 : CALL section_vals_get(qmmm_section, explicit=explicit)
860 394 : IF (explicit) THEN
861 394 : CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type)
862 : CALL setup_qm_atom_list(qmmm_section, qm_atom_index=qm_atom_index, qmmm_link=qmmm_link, &
863 : mm_link_atoms=mm_link_atoms, mm_link_scale_factor=mm_link_scale_factor, &
864 394 : fist_scale_charge_link=fist_scale_charge_link)
865 : !
866 : ! Do we want to use a different FF for the non-bonded QM/MM interactions?
867 : !
868 394 : qmmm_ff_section => section_vals_get_subs_vals(qmmm_section, "FORCEFIELD")
869 394 : CALL section_vals_get(qmmm_ff_section, explicit=qmmm_env%use_qmmm_ff)
870 394 : IF (qmmm_env%use_qmmm_ff) THEN
871 : CALL section_vals_val_get(qmmm_ff_section, "MULTIPLE_POTENTIAL", &
872 20 : l_val=qmmm_env%multiple_potential)
873 20 : CALL read_qmmm_ff_section(qmmm_ff_section, qmmm_env%inp_info)
874 : END IF
875 : END IF
876 394 : END SUBROUTINE setup_qmmm_vars_mm
877 :
878 : ! **************************************************************************************************
879 : !> \brief reads information regarding the forcefield specific for the QM/MM
880 : !> interactions
881 : !> \param qmmm_ff_section ...
882 : !> \param inp_info ...
883 : !> \par History
884 : !> 12.2004 created [tlaino]
885 : !> \author Teodoro Laino
886 : ! **************************************************************************************************
887 180 : SUBROUTINE read_qmmm_ff_section(qmmm_ff_section, inp_info)
888 : TYPE(section_vals_type), POINTER :: qmmm_ff_section
889 : TYPE(input_info_type), POINTER :: inp_info
890 :
891 : INTEGER :: n_gd, n_gp, n_lj, n_wl, np
892 : TYPE(section_vals_type), POINTER :: gd_section, gp_section, lj_section, &
893 : wl_section
894 :
895 : !
896 : ! NONBONDED
897 : !
898 :
899 20 : lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%LENNARD-JONES")
900 20 : wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%WILLIAMS")
901 20 : gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GOODWIN")
902 20 : gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GENPOT")
903 20 : CALL section_vals_get(lj_section, n_repetition=n_lj)
904 20 : np = n_lj
905 20 : IF (n_lj /= 0) THEN
906 18 : CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, lj_charmm=.TRUE.)
907 18 : CALL read_lj_section(inp_info%nonbonded, lj_section, start=0)
908 : END IF
909 20 : CALL section_vals_get(wl_section, n_repetition=n_wl)
910 20 : np = n_lj + n_wl
911 20 : IF (n_wl /= 0) THEN
912 2 : CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, williams=.TRUE.)
913 2 : CALL read_wl_section(inp_info%nonbonded, wl_section, start=n_lj)
914 : END IF
915 20 : CALL section_vals_get(gd_section, n_repetition=n_gd)
916 20 : np = n_lj + n_wl + n_gd
917 20 : IF (n_gd /= 0) THEN
918 0 : CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, goodwin=.TRUE.)
919 0 : CALL read_gd_section(inp_info%nonbonded, gd_section, start=n_lj + n_wl)
920 : END IF
921 20 : CALL section_vals_get(gp_section, n_repetition=n_gp)
922 20 : np = n_lj + n_wl + n_gd + n_gp
923 20 : IF (n_gp /= 0) THEN
924 0 : CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, gp=.TRUE.)
925 0 : CALL read_gp_section(inp_info%nonbonded, gp_section, start=n_lj + n_wl + n_gd)
926 : END IF
927 : !
928 : ! NONBONDED14
929 : !
930 20 : lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%LENNARD-JONES")
931 20 : wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%WILLIAMS")
932 20 : gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GOODWIN")
933 20 : gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GENPOT")
934 20 : CALL section_vals_get(lj_section, n_repetition=n_lj)
935 20 : np = n_lj
936 20 : IF (n_lj /= 0) THEN
937 0 : CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, lj_charmm=.TRUE.)
938 0 : CALL read_lj_section(inp_info%nonbonded14, lj_section, start=0)
939 : END IF
940 20 : CALL section_vals_get(wl_section, n_repetition=n_wl)
941 20 : np = n_lj + n_wl
942 20 : IF (n_wl /= 0) THEN
943 0 : CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, williams=.TRUE.)
944 0 : CALL read_wl_section(inp_info%nonbonded14, wl_section, start=n_lj)
945 : END IF
946 20 : CALL section_vals_get(gd_section, n_repetition=n_gd)
947 20 : np = n_lj + n_wl + n_gd
948 20 : IF (n_gd /= 0) THEN
949 0 : CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, goodwin=.TRUE.)
950 0 : CALL read_gd_section(inp_info%nonbonded14, gd_section, start=n_lj + n_wl)
951 : END IF
952 20 : CALL section_vals_get(gp_section, n_repetition=n_gp)
953 20 : np = n_lj + n_wl + n_gd + n_gp
954 20 : IF (n_gp /= 0) THEN
955 0 : CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, gp=.TRUE.)
956 0 : CALL read_gp_section(inp_info%nonbonded14, gp_section, start=n_lj + n_wl + n_gd)
957 : END IF
958 :
959 20 : END SUBROUTINE read_qmmm_ff_section
960 :
961 : ! **************************************************************************************************
962 : !> \brief ...
963 : !> \param qmmm_section ...
964 : !> \param qm_atom_index ...
965 : !> \param qm_atom_type ...
966 : !> \param mm_link_atoms ...
967 : !> \param mm_link_scale_factor ...
968 : !> \param qmmm_link ...
969 : !> \param fist_scale_charge_link ...
970 : !> \par History
971 : !> 12.2004 created [tlaino]
972 : !> \author Teodoro Laino
973 : ! **************************************************************************************************
974 2364 : SUBROUTINE setup_qm_atom_list(qmmm_section, qm_atom_index, qm_atom_type, &
975 : mm_link_atoms, mm_link_scale_factor, qmmm_link, fist_scale_charge_link)
976 : TYPE(section_vals_type), POINTER :: qmmm_section
977 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: qm_atom_index
978 : CHARACTER(len=default_string_length), &
979 : DIMENSION(:), OPTIONAL, POINTER :: qm_atom_type
980 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: mm_link_atoms
981 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: mm_link_scale_factor
982 : LOGICAL, INTENT(OUT), OPTIONAL :: qmmm_link
983 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: fist_scale_charge_link
984 :
985 : CHARACTER(len=default_string_length) :: qm_atom_kind, qm_link_element
986 : INTEGER :: ikind, k, link_involv_mm, link_type, &
987 : mm_index, n_var, nkind, nlinks, &
988 : num_qm_atom_tot
989 788 : INTEGER, DIMENSION(:), POINTER :: mm_indexes
990 : LOGICAL :: explicit
991 : REAL(KIND=dp) :: scale_f
992 : TYPE(section_vals_type), POINTER :: qm_kinds, qmmm_links
993 :
994 788 : num_qm_atom_tot = 0
995 788 : link_involv_mm = 0
996 788 : nlinks = 0
997 : !
998 : ! QM_KINDS
999 : !
1000 1576 : qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND")
1001 788 : CALL section_vals_get(qm_kinds, n_repetition=nkind)
1002 2596 : DO ikind = 1, nkind
1003 1808 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
1004 6336 : DO k = 1, n_var
1005 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
1006 3740 : i_vals=mm_indexes)
1007 5548 : num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
1008 : END DO
1009 : END DO
1010 : !
1011 : ! QM/MM LINKS
1012 : !
1013 788 : qmmm_links => section_vals_get_subs_vals(qmmm_section, "LINK")
1014 788 : CALL section_vals_get(qmmm_links, explicit=explicit)
1015 788 : IF (explicit) THEN
1016 128 : qmmm_link = .TRUE.
1017 128 : CALL section_vals_get(qmmm_links, n_repetition=nlinks)
1018 : ! Take care of the various link types
1019 524 : DO ikind = 1, nlinks
1020 : CALL section_vals_val_get(qmmm_links, "LINK_TYPE", i_rep_section=ikind, &
1021 396 : i_val=link_type)
1022 524 : SELECT CASE (link_type)
1023 : CASE (do_qmmm_link_imomm)
1024 388 : num_qm_atom_tot = num_qm_atom_tot + 1
1025 388 : link_involv_mm = link_involv_mm + 1
1026 : CASE (do_qmmm_link_pseudo)
1027 8 : num_qm_atom_tot = num_qm_atom_tot + 1
1028 : CASE (do_qmmm_link_gho)
1029 : ! do nothing for the moment
1030 : CASE DEFAULT
1031 396 : CPABORT("")
1032 : END SELECT
1033 : END DO
1034 : END IF
1035 788 : IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) &
1036 186 : ALLOCATE (mm_link_scale_factor(link_involv_mm))
1037 788 : IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) &
1038 186 : ALLOCATE (fist_scale_charge_link(link_involv_mm))
1039 788 : IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) &
1040 372 : ALLOCATE (mm_link_atoms(link_involv_mm))
1041 2364 : IF (PRESENT(qm_atom_index)) ALLOCATE (qm_atom_index(num_qm_atom_tot))
1042 1576 : IF (PRESENT(qm_atom_type)) ALLOCATE (qm_atom_type(num_qm_atom_tot))
1043 6572 : IF (PRESENT(qm_atom_index)) qm_atom_index = 0
1044 3680 : IF (PRESENT(qm_atom_type)) qm_atom_type = " "
1045 788 : num_qm_atom_tot = 1
1046 2596 : DO ikind = 1, nkind
1047 1808 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
1048 6336 : DO k = 1, n_var
1049 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
1050 3740 : i_vals=mm_indexes)
1051 3740 : IF (PRESENT(qm_atom_index)) THEN
1052 14516 : qm_atom_index(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = mm_indexes(:)
1053 : END IF
1054 3740 : IF (PRESENT(qm_atom_type)) THEN
1055 : CALL section_vals_val_get(qm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
1056 1870 : c_val=qm_atom_kind)
1057 4564 : qm_atom_type(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = qm_atom_kind
1058 : END IF
1059 5548 : num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
1060 : END DO
1061 : END DO
1062 982 : IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) mm_link_scale_factor = 0.0_dp
1063 982 : IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) fist_scale_charge_link = 0.0_dp
1064 1176 : IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) mm_link_atoms = 0
1065 788 : IF (explicit) THEN
1066 524 : DO ikind = 1, nlinks
1067 396 : IF (PRESENT(qm_atom_type)) THEN
1068 198 : CALL section_vals_val_get(qmmm_links, "QM_KIND", i_rep_section=ikind, c_val=qm_link_element)
1069 396 : qm_atom_type(num_qm_atom_tot:num_qm_atom_tot) = TRIM(qm_link_element)//"_LINK"
1070 : END IF
1071 396 : IF (PRESENT(qm_atom_index)) THEN
1072 396 : CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1073 23508 : CPASSERT(ALL(qm_atom_index /= mm_index))
1074 792 : qm_atom_index(num_qm_atom_tot:num_qm_atom_tot) = mm_index
1075 396 : num_qm_atom_tot = num_qm_atom_tot + 1
1076 : END IF
1077 396 : IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) THEN
1078 388 : CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1079 388 : mm_link_atoms(ikind) = mm_index
1080 : END IF
1081 396 : IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) THEN
1082 194 : CALL section_vals_val_get(qmmm_links, "QMMM_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
1083 194 : mm_link_scale_factor(ikind) = scale_f
1084 : END IF
1085 524 : IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) THEN
1086 194 : CALL section_vals_val_get(qmmm_links, "FIST_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
1087 194 : fist_scale_charge_link(ikind) = scale_f
1088 : END IF
1089 : END DO
1090 : END IF
1091 788 : CPASSERT(num_qm_atom_tot - 1 == SIZE(qm_atom_index))
1092 :
1093 788 : END SUBROUTINE setup_qm_atom_list
1094 :
1095 : ! **************************************************************************************************
1096 : !> \brief this routine sets up all variables to treat qmmm links
1097 : !> \param qmmm_section ...
1098 : !> \param qmmm_links ...
1099 : !> \param mm_el_pot_radius ...
1100 : !> \param mm_el_pot_radius_corr ...
1101 : !> \param mm_atom_index ...
1102 : !> \param iw ...
1103 : !> \par History
1104 : !> 12.2004 created [tlaino]
1105 : !> \author Teodoro Laino
1106 : ! **************************************************************************************************
1107 128 : SUBROUTINE setup_qmmm_links(qmmm_section, qmmm_links, mm_el_pot_radius, mm_el_pot_radius_corr, &
1108 : mm_atom_index, iw)
1109 : TYPE(section_vals_type), POINTER :: qmmm_section
1110 : TYPE(qmmm_links_type), POINTER :: qmmm_links
1111 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_el_pot_radius, mm_el_pot_radius_corr
1112 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1113 : INTEGER, INTENT(IN) :: iw
1114 :
1115 : INTEGER :: ikind, link_type, mm_index, n_gho, &
1116 : n_imomm, n_pseudo, n_rep_val, n_tot, &
1117 : nlinks, qm_index
1118 64 : INTEGER, DIMENSION(:), POINTER :: wrk_tmp
1119 : REAL(KIND=dp) :: alpha, my_radius
1120 : TYPE(section_vals_type), POINTER :: qmmm_link_section
1121 :
1122 64 : NULLIFY (wrk_tmp)
1123 64 : n_imomm = 0
1124 64 : n_gho = 0
1125 64 : n_pseudo = 0
1126 128 : qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
1127 64 : CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
1128 64 : CPASSERT(nlinks /= 0)
1129 262 : DO ikind = 1, nlinks
1130 198 : CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1131 198 : IF (link_type == do_qmmm_link_imomm) n_imomm = n_imomm + 1
1132 198 : IF (link_type == do_qmmm_link_gho) n_gho = n_gho + 1
1133 460 : IF (link_type == do_qmmm_link_pseudo) n_pseudo = n_pseudo + 1
1134 : END DO
1135 64 : n_tot = n_imomm + n_gho + n_pseudo
1136 64 : CPASSERT(n_tot /= 0)
1137 64 : ALLOCATE (qmmm_links)
1138 : NULLIFY (qmmm_links%imomm, &
1139 : qmmm_links%pseudo)
1140 : ! IMOMM
1141 64 : IF (n_imomm /= 0) THEN
1142 380 : ALLOCATE (qmmm_links%imomm(n_imomm))
1143 186 : ALLOCATE (wrk_tmp(n_imomm))
1144 256 : DO ikind = 1, n_imomm
1145 194 : NULLIFY (qmmm_links%imomm(ikind)%link)
1146 256 : ALLOCATE (qmmm_links%imomm(ikind)%link)
1147 : END DO
1148 62 : n_imomm = 0
1149 256 : DO ikind = 1, nlinks
1150 194 : CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1151 256 : IF (link_type == do_qmmm_link_imomm) THEN
1152 194 : n_imomm = n_imomm + 1
1153 194 : CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
1154 194 : CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1155 194 : CALL section_vals_val_get(qmmm_link_section, "ALPHA_IMOMM", i_rep_section=ikind, r_val=alpha)
1156 194 : CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
1157 194 : qmmm_links%imomm(n_imomm)%link%qm_index = qm_index
1158 194 : qmmm_links%imomm(n_imomm)%link%mm_index = mm_index
1159 194 : qmmm_links%imomm(n_imomm)%link%alpha = alpha
1160 194 : wrk_tmp(n_imomm) = mm_index
1161 194 : IF (n_rep_val == 1) THEN
1162 64 : CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, r_val=my_radius)
1163 996 : WHERE (mm_atom_index == mm_index) mm_el_pot_radius = my_radius
1164 996 : WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
1165 : END IF
1166 194 : CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
1167 194 : IF (n_rep_val == 1) THEN
1168 0 : CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, r_val=my_radius)
1169 0 : WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
1170 : END IF
1171 : END IF
1172 : END DO
1173 : !
1174 : ! Checking the link structure
1175 : !
1176 256 : DO ikind = 1, SIZE(wrk_tmp)
1177 3514 : IF (COUNT(wrk_tmp == wrk_tmp(ikind)) > 1) THEN
1178 0 : WRITE (iw, '(/A)') "In the IMOMM scheme no more than one QM atom can be bounded to the same MM atom."
1179 0 : WRITE (iw, '(A)') "Multiple link MM atom not allowed. Check your link sections."
1180 0 : CPABORT("")
1181 : END IF
1182 : END DO
1183 62 : DEALLOCATE (wrk_tmp)
1184 : END IF
1185 : ! PSEUDO
1186 64 : IF (n_pseudo /= 0) THEN
1187 10 : ALLOCATE (qmmm_links%pseudo(n_pseudo))
1188 6 : DO ikind = 1, n_pseudo
1189 4 : NULLIFY (qmmm_links%pseudo(ikind)%link)
1190 6 : ALLOCATE (qmmm_links%pseudo(ikind)%link)
1191 : END DO
1192 2 : n_pseudo = 0
1193 6 : DO ikind = 1, nlinks
1194 4 : CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1195 6 : IF (link_type == do_qmmm_link_pseudo) THEN
1196 4 : n_pseudo = n_pseudo + 1
1197 4 : CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
1198 4 : CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1199 4 : qmmm_links%pseudo(n_pseudo)%link%qm_index = qm_index
1200 4 : qmmm_links%pseudo(n_pseudo)%link%mm_index = mm_index
1201 : END IF
1202 : END DO
1203 : END IF
1204 : ! GHO
1205 64 : IF (n_gho /= 0) THEN
1206 : ! not yet implemented
1207 : ! still to define : type, implementation into QS
1208 0 : CPABORT("")
1209 : END IF
1210 64 : END SUBROUTINE setup_qmmm_links
1211 :
1212 : ! **************************************************************************************************
1213 : !> \brief this routine sets up all variables to treat qmmm links
1214 : !> \param qmmm_section ...
1215 : !> \param move_mm_charges ...
1216 : !> \param add_mm_charges ...
1217 : !> \param mm_atom_chrg ...
1218 : !> \param mm_el_pot_radius ...
1219 : !> \param mm_el_pot_radius_corr ...
1220 : !> \param added_charges ...
1221 : !> \param mm_atom_index ...
1222 : !> \par History
1223 : !> 12.2004 created [tlaino]
1224 : !> \author Teodoro Laino
1225 : ! **************************************************************************************************
1226 128 : SUBROUTINE move_or_add_atoms(qmmm_section, move_mm_charges, add_mm_charges, &
1227 : mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, &
1228 : added_charges, mm_atom_index)
1229 : TYPE(section_vals_type), POINTER :: qmmm_section
1230 : LOGICAL, INTENT(OUT) :: move_mm_charges, add_mm_charges
1231 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, &
1232 : mm_el_pot_radius_corr
1233 : TYPE(add_set_type), POINTER :: added_charges
1234 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1235 :
1236 : INTEGER :: i_add, icount, ikind, ind1, Index1, &
1237 : Index2, n_add_tot, n_adds, n_move_tot, &
1238 : n_moves, n_rep_val, nlinks
1239 : LOGICAL :: explicit
1240 : REAL(KIND=dp) :: alpha, c_radius, charge, radius
1241 : TYPE(section_vals_type), POINTER :: add_section, move_section, &
1242 : qmmm_link_section
1243 :
1244 64 : explicit = .FALSE.
1245 64 : move_mm_charges = .FALSE.
1246 64 : add_mm_charges = .FALSE.
1247 64 : NULLIFY (qmmm_link_section, move_section, add_section)
1248 64 : qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
1249 64 : CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
1250 64 : CPASSERT(nlinks /= 0)
1251 : icount = 0
1252 64 : n_move_tot = 0
1253 64 : n_add_tot = 0
1254 262 : DO ikind = 1, nlinks
1255 : move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
1256 198 : i_rep_section=ikind)
1257 198 : CALL section_vals_get(move_section, n_repetition=n_moves)
1258 : add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", &
1259 198 : i_rep_section=ikind)
1260 198 : CALL section_vals_get(add_section, n_repetition=n_adds)
1261 198 : n_move_tot = n_move_tot + n_moves
1262 460 : n_add_tot = n_add_tot + n_adds
1263 : END DO
1264 64 : icount = n_move_tot + n_add_tot
1265 64 : IF (n_add_tot /= 0) add_mm_charges = .TRUE.
1266 64 : IF (n_move_tot /= 0) move_mm_charges = .TRUE.
1267 : !
1268 : ! create add_set_type
1269 : !
1270 64 : CALL create_add_set_type(added_charges, ndim=icount)
1271 : !
1272 : ! Fill in structures
1273 : !
1274 64 : icount = 0
1275 262 : DO ikind = 1, nlinks
1276 : move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
1277 198 : i_rep_section=ikind)
1278 198 : CALL section_vals_get(move_section, explicit=explicit, n_repetition=n_moves)
1279 : !
1280 : ! Moving charge atoms
1281 : !
1282 198 : IF (explicit) THEN
1283 36 : DO i_add = 1, n_moves
1284 26 : icount = icount + 1
1285 26 : CALL section_vals_val_get(move_section, "ATOM_INDEX_1", i_val=Index1, i_rep_section=i_add)
1286 26 : CALL section_vals_val_get(move_section, "ATOM_INDEX_2", i_val=Index2, i_rep_section=i_add)
1287 26 : CALL section_vals_val_get(move_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
1288 26 : CALL section_vals_val_get(move_section, "RADIUS", r_val=radius, i_rep_section=i_add)
1289 26 : CALL section_vals_val_get(move_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
1290 26 : c_radius = radius
1291 26 : IF (n_rep_val == 1) &
1292 26 : CALL section_vals_val_get(move_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add)
1293 :
1294 : CALL set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, &
1295 : mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, &
1296 : mm_el_pot_radius_corr=mm_el_pot_radius_corr, &
1297 62 : mm_atom_index=mm_atom_index, move=n_moves, Ind1=ind1)
1298 : END DO
1299 10 : mm_atom_chrg(ind1) = 0.0_dp
1300 : END IF
1301 :
1302 : add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", &
1303 198 : i_rep_section=ikind)
1304 198 : CALL section_vals_get(add_section, explicit=explicit, n_repetition=n_adds)
1305 : !
1306 : ! Adding charge atoms
1307 : !
1308 460 : IF (explicit) THEN
1309 4 : DO i_add = 1, n_adds
1310 2 : icount = icount + 1
1311 2 : CALL section_vals_val_get(add_section, "ATOM_INDEX_1", i_val=Index1, i_rep_section=i_add)
1312 2 : CALL section_vals_val_get(add_section, "ATOM_INDEX_2", i_val=Index2, i_rep_section=i_add)
1313 2 : CALL section_vals_val_get(add_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
1314 2 : CALL section_vals_val_get(add_section, "RADIUS", r_val=radius, i_rep_section=i_add)
1315 2 : CALL section_vals_val_get(add_section, "CHARGE", r_val=charge, i_rep_section=i_add)
1316 2 : CALL section_vals_val_get(add_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
1317 2 : c_radius = radius
1318 2 : IF (n_rep_val == 1) &
1319 2 : CALL section_vals_val_get(add_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add)
1320 :
1321 : CALL set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, charge, &
1322 : mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, &
1323 : mm_el_pot_radius_corr=mm_el_pot_radius_corr, &
1324 6 : mm_atom_index=mm_atom_index)
1325 : END DO
1326 : END IF
1327 : END DO
1328 :
1329 64 : END SUBROUTINE move_or_add_atoms
1330 :
1331 : ! **************************************************************************************************
1332 : !> \brief this routine sets up all variables of the add_set_type type
1333 : !> \param added_charges ...
1334 : !> \param icount ...
1335 : !> \param Index1 ...
1336 : !> \param Index2 ...
1337 : !> \param alpha ...
1338 : !> \param radius ...
1339 : !> \param c_radius ...
1340 : !> \param charge ...
1341 : !> \param mm_atom_chrg ...
1342 : !> \param mm_el_pot_radius ...
1343 : !> \param mm_el_pot_radius_corr ...
1344 : !> \param mm_atom_index ...
1345 : !> \param move ...
1346 : !> \param ind1 ...
1347 : !> \par History
1348 : !> 12.2004 created [tlaino]
1349 : !> \author Teodoro Laino
1350 : ! **************************************************************************************************
1351 28 : SUBROUTINE set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, charge, &
1352 : mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, mm_atom_index, move, ind1)
1353 : TYPE(add_set_type), POINTER :: added_charges
1354 : INTEGER, INTENT(IN) :: icount, Index1, Index2
1355 : REAL(KIND=dp), INTENT(IN) :: alpha, radius, c_radius
1356 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: charge
1357 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, &
1358 : mm_el_pot_radius_corr
1359 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1360 : INTEGER, INTENT(in), OPTIONAL :: move
1361 : INTEGER, INTENT(OUT), OPTIONAL :: ind1
1362 :
1363 : INTEGER :: i, my_move
1364 : REAL(KIND=dp) :: my_c_radius, my_charge, my_radius
1365 :
1366 28 : my_move = 0
1367 28 : my_radius = radius
1368 28 : my_c_radius = c_radius
1369 28 : IF (PRESENT(charge)) my_charge = charge
1370 28 : IF (PRESENT(move)) my_move = move
1371 28 : i = 1
1372 60 : GetId: DO WHILE (i <= SIZE(mm_atom_index))
1373 60 : IF (Index1 == mm_atom_index(i)) EXIT GetId
1374 60 : i = i + 1
1375 : END DO GetId
1376 28 : IF (PRESENT(ind1)) ind1 = i
1377 28 : CPASSERT(i <= SIZE(mm_atom_index))
1378 28 : IF (.NOT. PRESENT(charge)) my_charge = mm_atom_chrg(i)/REAL(my_move, KIND=dp)
1379 28 : IF (my_radius == 0.0_dp) my_radius = mm_el_pot_radius(i)
1380 28 : IF (my_c_radius == 0.0_dp) my_c_radius = mm_el_pot_radius_corr(i)
1381 :
1382 28 : added_charges%add_env(icount)%Index1 = Index1
1383 28 : added_charges%add_env(icount)%Index2 = Index2
1384 28 : added_charges%add_env(icount)%alpha = alpha
1385 28 : added_charges%mm_atom_index(icount) = icount
1386 28 : added_charges%mm_atom_chrg(icount) = my_charge
1387 28 : added_charges%mm_el_pot_radius(icount) = my_radius
1388 28 : added_charges%mm_el_pot_radius_corr(icount) = my_c_radius
1389 28 : END SUBROUTINE set_add_set_type
1390 :
1391 : ! **************************************************************************************************
1392 : !> \brief this routine sets up the origin of the MM cell respect to the
1393 : !> origin of the QM cell. The origin of the QM cell is assumed to be
1394 : !> in (0.0,0.0,0.0)...
1395 : !> \param qmmm_section ...
1396 : !> \param qmmm_env ...
1397 : !> \param qm_cell_small ...
1398 : !> \param dr ...
1399 : !> \par History
1400 : !> 02.2005 created [tlaino]
1401 : !> \author Teodoro Laino
1402 : ! **************************************************************************************************
1403 788 : SUBROUTINE setup_origin_mm_cell(qmmm_section, qmmm_env, qm_cell_small, &
1404 : dr)
1405 : TYPE(section_vals_type), POINTER :: qmmm_section
1406 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1407 : TYPE(cell_type), POINTER :: qm_cell_small
1408 : REAL(KIND=dp), DIMENSION(3), INTENT(in) :: dr
1409 :
1410 : LOGICAL :: center_grid
1411 : REAL(KIND=dp), DIMENSION(3) :: tmp
1412 394 : REAL(KINd=dp), DIMENSION(:), POINTER :: vec
1413 :
1414 : ! This is the vector that corrects position to apply properly the PBC
1415 :
1416 394 : tmp(1) = qm_cell_small%hmat(1, 1)
1417 394 : tmp(2) = qm_cell_small%hmat(2, 2)
1418 394 : tmp(3) = qm_cell_small%hmat(3, 3)
1419 1576 : CPASSERT(ALL(tmp > 0))
1420 1576 : qmmm_env%dOmmOqm = tmp/2.0_dp
1421 : ! This is unit vector to translate the QM system in order to center it
1422 : ! in QM cell
1423 394 : CALL section_vals_val_get(qmmm_section, "CENTER_GRID", l_val=center_grid)
1424 394 : IF (center_grid) THEN
1425 72 : qmmm_env%utrasl = dr
1426 : ELSE
1427 1504 : qmmm_env%utrasl = 1.0_dp
1428 : END IF
1429 394 : CALL section_vals_val_get(qmmm_section, "INITIAL_TRANSLATION_VECTOR", r_vals=vec)
1430 3152 : qmmm_env%transl_v = vec
1431 394 : END SUBROUTINE setup_origin_mm_cell
1432 :
1433 : ! **************************************************************************************************
1434 : !> \brief this routine sets up list of MM atoms carrying an image charge
1435 : !> \param image_charge_section ...
1436 : !> \param qmmm_env ...
1437 : !> \param qm_atom_index ...
1438 : !> \param subsys_mm ...
1439 : !> \par History
1440 : !> 02.2012 created
1441 : !> \author Dorothea Golze
1442 : ! **************************************************************************************************
1443 10 : SUBROUTINE setup_image_atom_list(image_charge_section, qmmm_env, &
1444 : qm_atom_index, subsys_mm)
1445 :
1446 : TYPE(section_vals_type), POINTER :: image_charge_section
1447 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1448 : INTEGER, DIMENSION(:), POINTER :: qm_atom_index
1449 : TYPE(cp_subsys_type), POINTER :: subsys_mm
1450 :
1451 : INTEGER :: atom_a, atom_b, i, j, k, max_index, &
1452 : n_var, num_const_atom, &
1453 : num_image_mm_atom
1454 10 : INTEGER, DIMENSION(:), POINTER :: mm_indexes
1455 : LOGICAL :: fix_xyz, imageind_in_range
1456 10 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind
1457 :
1458 10 : NULLIFY (mm_indexes, molecule_kind)
1459 10 : imageind_in_range = .FALSE.
1460 10 : num_image_mm_atom = 0
1461 10 : max_index = 0
1462 :
1463 : CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1464 10 : n_rep_val=n_var)
1465 20 : DO i = 1, n_var
1466 : CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1467 10 : i_rep_val=i, i_vals=mm_indexes)
1468 20 : num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
1469 : END DO
1470 :
1471 30 : ALLOCATE (qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom))
1472 :
1473 30 : qmmm_env%image_charge_pot%image_mm_list = 0
1474 10 : num_image_mm_atom = 1
1475 :
1476 20 : DO i = 1, n_var
1477 : CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1478 10 : i_rep_val=i, i_vals=mm_indexes)
1479 : qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom:num_image_mm_atom &
1480 60 : + SIZE(mm_indexes) - 1) = mm_indexes(:)
1481 20 : num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
1482 : END DO
1483 :
1484 : ! checking, if in range, if list contains QM atoms or any atoms doubled
1485 10 : num_image_mm_atom = num_image_mm_atom - 1
1486 :
1487 10 : max_index = SIZE(subsys_mm%particles%els)
1488 :
1489 10 : CPASSERT(SIZE(qmmm_env%image_charge_pot%image_mm_list) /= 0)
1490 : imageind_in_range = (MAXVAL(qmmm_env%image_charge_pot%image_mm_list) <= max_index) &
1491 60 : .AND. (MINVAL(qmmm_env%image_charge_pot%image_mm_list) > 0)
1492 10 : CPASSERT(imageind_in_range)
1493 :
1494 30 : DO i = 1, num_image_mm_atom
1495 20 : atom_a = qmmm_env%image_charge_pot%image_mm_list(i)
1496 60 : IF (ANY(qm_atom_index == atom_a)) THEN
1497 0 : CPABORT("Image atom list must only contain MM atoms")
1498 : END IF
1499 40 : DO j = i + 1, num_image_mm_atom
1500 10 : atom_b = qmmm_env%image_charge_pot%image_mm_list(j)
1501 10 : IF (atom_a == atom_b) &
1502 20 : CPABORT("There are atoms doubled in image list.")
1503 : END DO
1504 : END DO
1505 :
1506 : ! check if molecules in list carry constraints
1507 10 : num_const_atom = 0
1508 10 : fix_xyz = .TRUE.
1509 10 : IF (ASSOCIATED(subsys_mm%molecule_kinds)) THEN
1510 10 : IF (ASSOCIATED(subsys_mm%molecule_kinds%els)) THEN
1511 10 : molecule_kind => subsys_mm%molecule_kinds%els
1512 78 : DO i = 1, SIZE(molecule_kind)
1513 76 : IF (.NOT. ASSOCIATED(molecule_kind(i)%fixd_list)) EXIT
1514 68 : IF (.NOT. fix_xyz) EXIT
1515 82 : DO j = 1, SIZE(molecule_kind(i)%fixd_list)
1516 4 : IF (.NOT. fix_xyz) EXIT
1517 80 : DO k = 1, num_image_mm_atom
1518 8 : atom_a = qmmm_env%image_charge_pot%image_mm_list(k)
1519 12 : IF (atom_a == molecule_kind(i)%fixd_list(j)%fixd) THEN
1520 4 : num_const_atom = num_const_atom + 1
1521 4 : IF (molecule_kind(i)%fixd_list(j)%itype /= use_perd_xyz) THEN
1522 : fix_xyz = .FALSE.
1523 : EXIT
1524 : END IF
1525 : END IF
1526 : END DO
1527 : END DO
1528 : END DO
1529 : END IF
1530 : END IF
1531 :
1532 : ! if all image atoms are constrained, calculate image matrix only
1533 : ! once for the first MD or GEO_OPT step (for non-iterative case)
1534 10 : IF (num_const_atom == num_image_mm_atom .AND. fix_xyz) THEN
1535 2 : qmmm_env%image_charge_pot%state_image_matrix = calc_once
1536 : ELSE
1537 8 : qmmm_env%image_charge_pot%state_image_matrix = calc_always
1538 : END IF
1539 :
1540 10 : END SUBROUTINE setup_image_atom_list
1541 :
1542 : ! **************************************************************************************************
1543 : !> \brief Print info on image charges
1544 : !> \param qmmm_env ...
1545 : !> \param qmmm_section ...
1546 : !> \par History
1547 : !> 03.2012 created
1548 : !> \author Dorothea Golze
1549 : ! **************************************************************************************************
1550 10 : SUBROUTINE print_image_charge_info(qmmm_env, qmmm_section)
1551 :
1552 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1553 : TYPE(section_vals_type), POINTER :: qmmm_section
1554 :
1555 : INTEGER :: iw
1556 : REAL(KIND=dp) :: eta, eta_conv, V0, V0_conv
1557 : TYPE(cp_logger_type), POINTER :: logger
1558 :
1559 10 : logger => cp_get_default_logger()
1560 : iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%PROGRAM_RUN_INFO", &
1561 10 : extension=".log")
1562 10 : eta = qmmm_env%image_charge_pot%eta
1563 10 : eta_conv = cp_unit_from_cp2k(eta, "angstrom", power=-2)
1564 10 : V0 = qmmm_env%image_charge_pot%V0
1565 10 : V0_conv = cp_unit_from_cp2k(V0, "volt")
1566 :
1567 10 : IF (iw > 0) THEN
1568 5 : WRITE (iw, FMT="(T25,A)") "IMAGE CHARGE PARAMETERS"
1569 5 : WRITE (iw, FMT="(T25,A)") REPEAT("-", 23)
1570 5 : WRITE (iw, FMT="(/)")
1571 5 : WRITE (iw, FMT="(T2,A)") "INDEX OF MM ATOMS CARRYING AN IMAGE CHARGE:"
1572 5 : WRITE (iw, FMT="(/)")
1573 :
1574 15 : WRITE (iw, "(7X,10I6)") qmmm_env%image_charge_pot%image_mm_list
1575 5 : WRITE (iw, FMT="(/)")
1576 : WRITE (iw, "(T2,A52,T69,F12.8)") &
1577 5 : "WIDTH OF GAUSSIAN CHARGE DISTRIBUTION [angstrom^-2]:", eta_conv
1578 5 : WRITE (iw, "(T2,A26,T69,F12.8)") "EXTERNAL POTENTIAL [volt]:", V0_conv
1579 5 : WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 79)
1580 : END IF
1581 : CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
1582 10 : "PRINT%PROGRAM_RUN_INFO")
1583 :
1584 10 : END SUBROUTINE print_image_charge_info
1585 :
1586 : END MODULE qmmm_init
|
