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