(git:98357aa)
Loading...
Searching...
No Matches
qmmm_init.F
Go to the documentation of this file.
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! **************************************************************************************************
19 USE cell_types, ONLY: cell_type,&
30 USE cp_units, ONLY: cp_unit_from_cp2k,&
40 USE input_constants, ONLY: &
48 USE kinds, ONLY: default_string_length,&
49 dp
55 USE pw_env_types, ONLY: pw_env_type
61 USE qmmm_types_low, ONLY: add_set_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
92
93CONTAINS
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 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 TYPE(particle_type), DIMENSION(:), POINTER :: core_set, particle_set, shell_set
131 TYPE(shell_kind_type), POINTER :: shell_kind
132
133 NULLIFY (particle_set, my_kind, added_shells)
134 CALL cp_subsys_get(subsys=subsys, particles=particles, core_particles=core_particles, &
135 shell_particles=shell_particles)
136 particle_set => particles%els
137
138 IF (all(particle_set(:)%shell_index == 0)) THEN
139 shell_model = .false.
140 CALL create_add_shell_type(added_shells, ndim=0)
141 ELSE
142 shell_model = .true.
143 END IF
144
145 IF (shell_model) THEN
146 shell_set => shell_particles%els
147 core_set => core_particles%els
148 ishell = SIZE(shell_set)
149 CALL create_add_shell_type(added_shells, ndim=ishell)
150 added_shells%added_particles => shell_set
151 added_shells%added_cores => core_set
152 END IF
153
154 DO i = 1, SIZE(mm_atom_index)
155 indmm = mm_atom_index(i)
156 my_kind => particle_set(indmm)%atomic_kind
157 CALL get_atomic_kind(atomic_kind=my_kind, fist_potential=my_potential, &
158 shell_active=is_shell, shell=shell_kind)
159 CALL get_potential(potential=my_potential, &
160 qeff=qi, &
161 qmmm_radius=ri, &
162 qmmm_corr_radius=rc)
163 IF (ASSOCIATED(charges)) qi = charges(indmm)
164 mm_atom_chrg(i) = qi
165 mm_el_pot_radius(i) = ri
166 mm_el_pot_radius_corr(i) = rc
167 IF (is_shell) THEN
168 indshell = particle_set(indmm)%shell_index
169 IF (ASSOCIATED(shell_kind)) THEN
170 CALL get_shell(shell=shell_kind, &
171 charge_core=qcore, &
172 charge_shell=qshell)
173 mm_atom_chrg(i) = qcore
174 END IF
175 added_shells%mm_core_index(indshell) = indmm
176 added_shells%mm_core_chrg(indshell) = qshell
177 added_shells%mm_el_pot_radius(indshell) = ri*1.0_dp
178 added_shells%mm_el_pot_radius_corr(indshell) = rc*1.0_dp
179 END IF
180 END DO
181
182 IF (ASSOCIATED(mm_link_atoms)) THEN
183 DO ilink = 1, SIZE(mm_link_atoms)
184 DO i = 1, SIZE(mm_atom_index)
185 IF (mm_atom_index(i) == mm_link_atoms(ilink)) EXIT
186 END DO
187 indmm = mm_atom_index(i)
188 mm_atom_chrg(i) = mm_atom_chrg(i)*mm_link_scale_factor(ilink)
189 END DO
190 END IF
191
192 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 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 qtot = 0.0_dp
224 logger => cp_get_default_logger()
225 iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%QMMM_CHARGES", &
226 extension=".log")
227 IF (iw > 0) THEN
228 WRITE (iw, fmt="(/,T2,A)") repeat("-", 79)
229 WRITE (iw, fmt='(/5X,A)') "MM POINT CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
230 WRITE (iw, fmt="(/,T2,A)") repeat("-", 79)
231 DO i = 1, SIZE(mm_atom_index)
232 indmm = mm_atom_index(i)
233 qi = mm_atom_chrg(i)
234 qtot = qtot + qi
235 ri = mm_el_pot_radius(i)
236 rc = mm_el_pot_radius_corr(i)
237 IF (nocompatibility) THEN
238 WRITE (iw, '(5X,A9,T15,I5,T28,A8,T38,F12.6,T60,A8,T69,F12.6)') ' MM ATOM:', indmm, ' RADIUS:', ri, &
239 ' 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 ' MM ATOM:', indmm, ' RADIUS:', ri, ' CHARGE:', qi, 'CORR. RADIUS', rc
243 END IF
244 END DO
245 IF (added_charges%num_mm_atoms /= 0) THEN
246 WRITE (iw, fmt="(/,T2,A)") repeat("-", 79)
247 WRITE (iw, '(/5X,A)') "ADDED POINT CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
248 WRITE (iw, fmt="(/,T2,A)") repeat("-", 79)
249 DO i = 1, SIZE(added_charges%mm_atom_index)
250 indmm = added_charges%mm_atom_index(i)
251 qi = added_charges%mm_atom_chrg(i)
252 qtot = qtot + qi
253 ri = added_charges%mm_el_pot_radius(i)
254 ind1 = added_charges%add_env(i)%Index1
255 ind2 = added_charges%add_env(i)%Index2
256 IF (nocompatibility) THEN
257 WRITE (iw, '(5X,A9,I5,T25,A8,T35,F12.6,T50,A8,T59,F12.6,I5,I5)') 'MM POINT:', indmm, ' RADIUS:', ri, &
258 ' 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 'MM POINT:', indmm, ' RADIUS:', ri, ' CHARGE:', qi, ind1, ind2, 'CORR. RADIUS', rc
262 END IF
263 END DO
264
265 END IF
266
267 IF (shell_model) THEN
268 WRITE (iw, fmt="(/,T2,A)") repeat("-", 73)
269 WRITE (iw, '(/5X,A)') "ADDED SHELL CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
270 WRITE (iw, fmt="(/,T2,A)") repeat("-", 73)
271
272 DO i = 1, SIZE(added_shells%mm_core_index)
273 indmm = added_shells%mm_core_index(i)
274 qi = added_shells%mm_core_chrg(i)
275 qtot = qtot + qi
276 ri = added_shells%mm_el_pot_radius(i)
277 IF (nocompatibility) THEN
278 WRITE (iw, '(7X,A,I5,A8,F12.6,A8,F12.6,3F12.6)') 'SHELL:', indmm, ' RADIUS:', ri, &
279 ' CHARGE:', qi, added_shells%added_particles(i)%r
280 ELSE
281 WRITE (iw, '(7X,A,I5,A8,F12.6,A8,F12.6,A,F12.6)') 'SHELL:', indmm, ' RADIUS:', ri, &
282 ' CHARGE:', qi, ' CORR. RADIUS', rc
283 END IF
284
285 END DO
286
287 END IF
288
289 WRITE (iw, fmt="(/,T2,A)") repeat("-", 79)
290 WRITE (iw, '(/,T50,A,T69,F12.6)') ' TOTAL CHARGE:', qtot
291 WRITE (iw, fmt="(/,T2,A,/)") repeat("-", 79)
292 END IF
293 CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
294 "PRINT%QMMM_CHARGES")
295 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 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 logger => cp_get_default_logger()
314 iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%qmmm_link_info", extension=".log")
315 IF (iw > 0) THEN
316 IF (ASSOCIATED(qmmm_links)) THEN
317 WRITE (iw, fmt="(/,T2, A)") repeat("-", 73)
318 WRITE (iw, fmt="(/,T31,A)") " QM/MM LINKS "
319 WRITE (iw, fmt="(/,T2, A)") repeat("-", 73)
320 IF (ASSOCIATED(qmmm_links%imomm)) THEN
321 WRITE (iw, fmt="(/,T31,A)") " IMOMM TYPE LINK "
322 DO i = 1, SIZE(qmmm_links%imomm)
323 qm_index = qmmm_links%imomm(i)%link%qm_index
324 mm_index = qmmm_links%imomm(i)%link%mm_index
325 alpha = qmmm_links%imomm(i)%link%alpha
326 WRITE (iw, fmt="(T2,A,T20,A9,I8,1X,A9,I8,T62,A6,F12.6)") "TYPE: IMOMM", &
327 "QM INDEX:", qm_index, "MM INDEX:", mm_index, "ALPHA:", alpha
328 END DO
329 END IF
330 IF (ASSOCIATED(qmmm_links%pseudo)) THEN
331 WRITE (iw, fmt="(/,T31,A)") " PSEUDO TYPE LINK "
332 DO i = 1, SIZE(qmmm_links%pseudo)
333 qm_index = qmmm_links%pseudo(i)%link%qm_index
334 mm_index = qmmm_links%pseudo(i)%link%mm_index
335 WRITE (iw, fmt="(T2,A,T20,A9,I8,1X,A9,I8)") "TYPE: PSEUDO", &
336 "QM INDEX:", qm_index, "MM INDEX:", mm_index
337 END DO
338 END IF
339 WRITE (iw, fmt="(/,T2,A,/)") repeat("-", 73)
340 ELSE
341 WRITE (iw, fmt="(/,T2, A)") repeat("-", 73)
342 WRITE (iw, fmt="(/,T26,A)") " NO QM/MM LINKS DETECTED"
343 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 "PRINT%qmmm_link_info")
348 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: maxradius, maxradius2
378 TYPE(pw_env_type), POINTER :: pw_env
379
380 NULLIFY (maxradius, maxradius2, pw_env)
381
382 maxchrg = maxval(abs(mm_atom_chrg(:)))
383 CALL get_qs_env(qs_env, pw_env=pw_env)
384 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 qmmm_section=qmmm_section)
397
398 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 qmmm_section=qmmm_section)
411
412 SELECT CASE (qmmm_env_qm%qmmm_coupl_type)
414 DO i = 1, SIZE(maxradius)
415 maxradius(i) = max(maxradius(i), maxradius2(i))
416 END DO
417 END SELECT
418
419 IF (ASSOCIATED(maxradius2)) DEALLOCATE (maxradius2)
420 END IF
421
422 IF (qmmm_env_qm%added_shells%num_mm_atoms > 0) THEN
423
424 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 qmmm_section=qmmm_section)
437
438 SELECT CASE (qmmm_env_qm%qmmm_coupl_type)
440 DO i = 1, SIZE(maxradius)
441 maxradius(i) = max(maxradius(i), maxradius2(i))
442 END DO
443 END SELECT
444
445 IF (ASSOCIATED(maxradius2)) DEALLOCATE (maxradius2)
446
447 END IF
448
449 qmmm_env_qm%maxradius => maxradius
450
451 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 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 print_section=print_section)
479
480 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 print_section=print_section)
489 END IF
490
491 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 print_section=print_section)
500 END IF
501
502 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 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 IF (qmmm_env_qm%periodic) THEN
535
536 NULLIFY (dft_control)
537 CALL get_qs_env(qs_env, dft_control=dft_control)
538
539 IF (dft_control%qs_control%semi_empirical) THEN
540 cpabort("QM/MM periodic calculations not implemented for semi empirical methods")
541 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 para_env=para_env, qmmm_periodic=qmmm_periodic, print_section=print_section)
545 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 para_env=para_env, qmmm_periodic=qmmm_periodic, print_section=print_section)
549 ELSE
550
551 ! setup for GPW/GPAW
552 maxchrg = maxval(abs(mm_atom_chrg(:)))
553 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 ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
568
569 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 ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
584 END IF
585
586 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 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 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 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 INTEGER, DIMENSION(:), POINTER :: mm_link_atoms
645 LOGICAL :: explicit, is_mm, is_qm
646 REAL(kind=dp) :: tmp_radius, tmp_radius_c
647 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 NULLIFY (mm_link_atoms, tmp_sph_cut)
655 NULLIFY (image_charge_section)
656 qmmm_link = .false.
657
658 CALL section_vals_get(qmmm_section, explicit=explicit)
659 IF (explicit) THEN
660 CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type)
661 CALL section_vals_val_get(qmmm_section, "EPS_MM_RSPACE", r_val=eps_mm_rspace)
662 CALL section_vals_val_get(qmmm_section, "SPHERICAL_CUTOFF", r_vals=tmp_sph_cut)
663 cpassert(SIZE(tmp_sph_cut) == 2)
664 qmmm_env%spherical_cutoff = tmp_sph_cut
665 IF (qmmm_env%spherical_cutoff(1) <= 0.0_dp) THEN
666 qmmm_env%spherical_cutoff(2) = 0.0_dp
667 ELSE
668 IF (qmmm_env%spherical_cutoff(2) <= 0.0_dp) qmmm_env%spherical_cutoff(2) = epsilon(0.0_dp)
669 tmp_radius = qmmm_env%spherical_cutoff(1) - 20.0_dp*qmmm_env%spherical_cutoff(2)
670 IF (tmp_radius <= 0.0_dp) THEN
671 CALL cp_abort(__location__, &
672 "SPHERICAL_CUTOFF(1) > 20*SPHERICAL_CUTOFF(1)! Please correct parameters for "// &
673 "the Spherical Cutoff in order to satisfy the previous condition!")
674 END IF
675 END IF
676 !
677 ! Initialization of arrays and core_charge_radius...
678 !
679 tmp_radius = 0.0_dp
680 CALL cp_subsys_get(subsys=subsys_mm, atomic_kinds=atomic_kinds)
681 DO ikind = 1, SIZE(atomic_kinds%els)
682 atomic_kind => atomic_kinds%els(ikind)
683 CALL get_atomic_kind(atomic_kind=atomic_kind, &
684 fist_potential=fist_potential)
685 CALL set_potential(potential=fist_potential, &
686 qmmm_radius=tmp_radius, &
687 qmmm_corr_radius=tmp_radius)
688 CALL set_atomic_kind(atomic_kind=atomic_kind, &
689 fist_potential=fist_potential)
690 END DO
691 CALL setup_qm_atom_list(qmmm_section=qmmm_section, &
692 qm_atom_index=qm_atom_index, &
693 qm_atom_type=qm_atom_type, &
694 mm_link_atoms=mm_link_atoms, &
695 qmmm_link=qmmm_link)
696 !
697 ! MM_KINDS
698 !
699 mm_kinds => section_vals_get_subs_vals(qmmm_section, "MM_KIND")
700 CALL section_vals_get(mm_kinds, explicit=explicit, n_repetition=nkind)
701 !
702 ! Default
703 !
704 tmp_radius = cp_unit_to_cp2k(radius_qmmm_default, "angstrom")
705 set_radius_pot_0: DO ikindr = 1, SIZE(atomic_kinds%els)
706 atomic_kind => atomic_kinds%els(ikindr)
707 CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
708 CALL get_atomic_kind(atomic_kind=atomic_kind, &
709 fist_potential=fist_potential)
710 CALL set_potential(potential=fist_potential, qmmm_radius=tmp_radius, &
711 qmmm_corr_radius=tmp_radius)
712 CALL set_atomic_kind(atomic_kind=atomic_kind, &
713 fist_potential=fist_potential)
714 END DO set_radius_pot_0
715 !
716 ! If present overwrite the kind specified in input file...
717 !
718 IF (explicit) THEN
719 DO ikind = 1, nkind
720 CALL section_vals_val_get(mm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
721 c_val=mm_atom_kind)
722 CALL section_vals_val_get(mm_kinds, "RADIUS", i_rep_section=ikind, r_val=tmp_radius)
723 tmp_radius_c = tmp_radius
724 CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
725 IF (n_rep_val == 1) CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, &
726 r_val=tmp_radius_c)
727 set_radius_pot_1: DO ikindr = 1, SIZE(atomic_kinds%els)
728 atomic_kind => atomic_kinds%els(ikindr)
729 CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
730 is_qm = qmmm_ff_precond_only_qm(atmname)
731 IF (trim(mm_atom_kind) == atmname) THEN
732 CALL get_atomic_kind(atomic_kind=atomic_kind, &
733 fist_potential=fist_potential)
734 CALL set_potential(potential=fist_potential, &
735 qmmm_radius=tmp_radius, &
736 qmmm_corr_radius=tmp_radius_c)
737 CALL set_atomic_kind(atomic_kind=atomic_kind, &
738 fist_potential=fist_potential)
739 END IF
740 END DO set_radius_pot_1
741 END DO
742 END IF
743
744 !Image charge section
745
746 image_charge_section => section_vals_get_subs_vals(qmmm_section, "IMAGE_CHARGE")
747 CALL section_vals_get(image_charge_section, explicit=qmmm_env%image_charge)
748
749 ELSE
750 cpabort("QMMM section not present in input file!")
751 END IF
752 !
753 ! Build MM atoms list
754 !
755 size_mm_system = SIZE(subsys_mm%particles%els) - SIZE(qm_atom_index)
756 IF (qmmm_link .AND. ASSOCIATED(mm_link_atoms)) size_mm_system = size_mm_system + SIZE(mm_link_atoms)
757 ALLOCATE (mm_atom_index(size_mm_system))
758 icount = 0
759
760 DO i = 1, SIZE(subsys_mm%particles%els)
761 is_mm = .true.
762 IF (any(qm_atom_index == i)) THEN
763 is_mm = .false.
764 END IF
765 IF (ASSOCIATED(mm_link_atoms)) THEN
766 IF (any(mm_link_atoms == i) .AND. qmmm_link) is_mm = .true.
767 END IF
768 IF (is_mm) THEN
769 icount = icount + 1
770 IF (icount <= size_mm_system) mm_atom_index(icount) = i
771 END IF
772 END DO
773 cpassert(icount == size_mm_system)
774 IF (ASSOCIATED(mm_link_atoms)) THEN
775 DEALLOCATE (mm_link_atoms)
776 END IF
777
778 ! Build image charge atom list + set up variables
779 !
780 IF (qmmm_env%image_charge) THEN
781 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
782 explicit=explicit)
783 IF (explicit) qmmm_env%image_charge_pot%all_mm = .false.
784
785 IF (qmmm_env%image_charge_pot%all_mm) THEN
786 qmmm_env%image_charge_pot%image_mm_list => mm_atom_index
787 ELSE
788 CALL setup_image_atom_list(image_charge_section, qmmm_env, &
789 qm_atom_index, subsys_mm)
790 END IF
791
792 qmmm_env%image_charge_pot%particles_all => subsys_mm%particles%els
793
794 CALL section_vals_val_get(image_charge_section, "EXT_POTENTIAL", &
795 r_val=qmmm_env%image_charge_pot%V0)
796 CALL section_vals_val_get(image_charge_section, "WIDTH", &
797 r_val=qmmm_env%image_charge_pot%eta)
798 CALL section_vals_val_get(image_charge_section, "DETERM_COEFF", &
799 i_val=my_type)
800 SELECT CASE (my_type)
802 qmmm_env%image_charge_pot%coeff_iterative = .false.
803 CASE (do_qmmm_image_iter)
804 qmmm_env%image_charge_pot%coeff_iterative = .true.
805 END SELECT
806
807 CALL section_vals_val_get(image_charge_section, "RESTART_IMAGE_MATRIX", &
808 l_val=qmmm_env%image_charge_pot%image_restart)
809
810 CALL section_vals_val_get(image_charge_section, "IMAGE_MATRIX_METHOD", &
811 i_val=qmmm_env%image_charge_pot%image_matrix_method)
812
813 IF (qmmm_env%image_charge_pot%image_matrix_method == do_eri_mme) THEN
814 eri_mme_section => section_vals_get_subs_vals(image_charge_section, "ERI_MME")
815 CALL cp_eri_mme_init_read_input(eri_mme_section, qmmm_env%image_charge_pot%eri_mme_param)
816 CALL cp_eri_mme_set_params(qmmm_env%image_charge_pot%eri_mme_param, &
817 hmat=qm_cell_small%hmat, is_ortho=qm_cell_small%orthorhombic, &
818 zet_min=qmmm_env%image_charge_pot%eta, &
819 zet_max=qmmm_env%image_charge_pot%eta, &
820 l_max_zet=0, &
821 l_max=0, &
822 para_env=para_env)
823
824 END IF
825 END IF
826
827 END SUBROUTINE setup_qmmm_vars_qm
828
829! **************************************************************************************************
830!> \brief ...
831!> \param qmmm_section ...
832!> \param qmmm_env ...
833!> \param qm_atom_index ...
834!> \param mm_link_atoms ...
835!> \param mm_link_scale_factor ...
836!> \param fist_scale_charge_link ...
837!> \param qmmm_coupl_type ...
838!> \param qmmm_link ...
839!> \par History
840!> 12.2004 created [tlaino]
841!> \author Teodoro Laino
842! **************************************************************************************************
843 SUBROUTINE setup_qmmm_vars_mm(qmmm_section, qmmm_env, qm_atom_index, &
844 mm_link_atoms, mm_link_scale_factor, &
845 fist_scale_charge_link, qmmm_coupl_type, &
846 qmmm_link)
847 TYPE(section_vals_type), POINTER :: qmmm_section
848 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
849 INTEGER, DIMENSION(:), POINTER :: qm_atom_index, mm_link_atoms
850 REAL(kind=dp), DIMENSION(:), POINTER :: mm_link_scale_factor, &
851 fist_scale_charge_link
852 INTEGER, INTENT(OUT) :: qmmm_coupl_type
853 LOGICAL, INTENT(OUT) :: qmmm_link
854
855 LOGICAL :: explicit
856 TYPE(section_vals_type), POINTER :: qmmm_ff_section
857
858 NULLIFY (qmmm_ff_section)
859 qmmm_link = .false.
860 CALL section_vals_get(qmmm_section, explicit=explicit)
861 IF (explicit) THEN
862 CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type)
863 CALL setup_qm_atom_list(qmmm_section, qm_atom_index=qm_atom_index, qmmm_link=qmmm_link, &
864 mm_link_atoms=mm_link_atoms, mm_link_scale_factor=mm_link_scale_factor, &
865 fist_scale_charge_link=fist_scale_charge_link)
866 !
867 ! Do we want to use a different FF for the non-bonded QM/MM interactions?
868 !
869 qmmm_ff_section => section_vals_get_subs_vals(qmmm_section, "FORCEFIELD")
870 CALL section_vals_get(qmmm_ff_section, explicit=qmmm_env%use_qmmm_ff)
871 IF (qmmm_env%use_qmmm_ff) THEN
872 CALL section_vals_val_get(qmmm_ff_section, "MULTIPLE_POTENTIAL", &
873 l_val=qmmm_env%multiple_potential)
874 CALL read_qmmm_ff_section(qmmm_ff_section, qmmm_env%inp_info)
875 END IF
876 END IF
877 END SUBROUTINE setup_qmmm_vars_mm
878
879! **************************************************************************************************
880!> \brief reads information regarding the forcefield specific for the QM/MM
881!> interactions
882!> \param qmmm_ff_section ...
883!> \param inp_info ...
884!> \par History
885!> 12.2004 created [tlaino]
886!> \author Teodoro Laino
887! **************************************************************************************************
888 SUBROUTINE read_qmmm_ff_section(qmmm_ff_section, inp_info)
889 TYPE(section_vals_type), POINTER :: qmmm_ff_section
890 TYPE(input_info_type), POINTER :: inp_info
891
892 INTEGER :: n_gd, n_gp, n_lj, n_wl, np
893 TYPE(section_vals_type), POINTER :: gd_section, gp_section, lj_section, &
894 wl_section
895
896!
897! NONBONDED
898!
899
900 lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%LENNARD-JONES")
901 wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%WILLIAMS")
902 gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GOODWIN")
903 gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GENPOT")
904 CALL section_vals_get(lj_section, n_repetition=n_lj)
905 np = n_lj
906 IF (n_lj /= 0) THEN
907 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, lj_charmm=.true.)
908 CALL read_lj_section(inp_info%nonbonded, lj_section, start=0)
909 END IF
910 CALL section_vals_get(wl_section, n_repetition=n_wl)
911 np = n_lj + n_wl
912 IF (n_wl /= 0) THEN
913 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, williams=.true.)
914 CALL read_wl_section(inp_info%nonbonded, wl_section, start=n_lj)
915 END IF
916 CALL section_vals_get(gd_section, n_repetition=n_gd)
917 np = n_lj + n_wl + n_gd
918 IF (n_gd /= 0) THEN
919 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, goodwin=.true.)
920 CALL read_gd_section(inp_info%nonbonded, gd_section, start=n_lj + n_wl)
921 END IF
922 CALL section_vals_get(gp_section, n_repetition=n_gp)
923 np = n_lj + n_wl + n_gd + n_gp
924 IF (n_gp /= 0) THEN
925 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, gp=.true.)
926 CALL read_gp_section(inp_info%nonbonded, gp_section, start=n_lj + n_wl + n_gd)
927 END IF
928 !
929 ! NONBONDED14
930 !
931 lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%LENNARD-JONES")
932 wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%WILLIAMS")
933 gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GOODWIN")
934 gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GENPOT")
935 CALL section_vals_get(lj_section, n_repetition=n_lj)
936 np = n_lj
937 IF (n_lj /= 0) THEN
938 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, lj_charmm=.true.)
939 CALL read_lj_section(inp_info%nonbonded14, lj_section, start=0)
940 END IF
941 CALL section_vals_get(wl_section, n_repetition=n_wl)
942 np = n_lj + n_wl
943 IF (n_wl /= 0) THEN
944 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, williams=.true.)
945 CALL read_wl_section(inp_info%nonbonded14, wl_section, start=n_lj)
946 END IF
947 CALL section_vals_get(gd_section, n_repetition=n_gd)
948 np = n_lj + n_wl + n_gd
949 IF (n_gd /= 0) THEN
950 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, goodwin=.true.)
951 CALL read_gd_section(inp_info%nonbonded14, gd_section, start=n_lj + n_wl)
952 END IF
953 CALL section_vals_get(gp_section, n_repetition=n_gp)
954 np = n_lj + n_wl + n_gd + n_gp
955 IF (n_gp /= 0) THEN
956 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, gp=.true.)
957 CALL read_gp_section(inp_info%nonbonded14, gp_section, start=n_lj + n_wl + n_gd)
958 END IF
959
960 END SUBROUTINE read_qmmm_ff_section
961
962! **************************************************************************************************
963!> \brief ...
964!> \param qmmm_section ...
965!> \param qm_atom_index ...
966!> \param qm_atom_type ...
967!> \param mm_link_atoms ...
968!> \param mm_link_scale_factor ...
969!> \param qmmm_link ...
970!> \param fist_scale_charge_link ...
971!> \par History
972!> 12.2004 created [tlaino]
973!> \author Teodoro Laino
974! **************************************************************************************************
975 SUBROUTINE setup_qm_atom_list(qmmm_section, qm_atom_index, qm_atom_type, &
976 mm_link_atoms, mm_link_scale_factor, qmmm_link, fist_scale_charge_link)
977 TYPE(section_vals_type), POINTER :: qmmm_section
978 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: qm_atom_index
979 CHARACTER(len=default_string_length), &
980 DIMENSION(:), OPTIONAL, POINTER :: qm_atom_type
981 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: mm_link_atoms
982 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: mm_link_scale_factor
983 LOGICAL, INTENT(OUT), OPTIONAL :: qmmm_link
984 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: fist_scale_charge_link
985
986 CHARACTER(len=default_string_length) :: qm_atom_kind, qm_link_element
987 INTEGER :: ikind, k, link_involv_mm, link_type, &
988 mm_index, n_var, nkind, nlinks, &
989 num_qm_atom_tot
990 INTEGER, DIMENSION(:), POINTER :: mm_indexes
991 LOGICAL :: explicit
992 REAL(kind=dp) :: scale_f
993 TYPE(section_vals_type), POINTER :: qm_kinds, qmmm_links
994
995 num_qm_atom_tot = 0
996 link_involv_mm = 0
997 nlinks = 0
998 !
999 ! QM_KINDS
1000 !
1001 qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND")
1002 CALL section_vals_get(qm_kinds, n_repetition=nkind)
1003 DO ikind = 1, nkind
1004 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
1005 DO k = 1, n_var
1006 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
1007 i_vals=mm_indexes)
1008 num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
1009 END DO
1010 END DO
1011 !
1012 ! QM/MM LINKS
1013 !
1014 qmmm_links => section_vals_get_subs_vals(qmmm_section, "LINK")
1015 CALL section_vals_get(qmmm_links, explicit=explicit)
1016 IF (explicit) THEN
1017 qmmm_link = .true.
1018 CALL section_vals_get(qmmm_links, n_repetition=nlinks)
1019 ! Take care of the various link types
1020 DO ikind = 1, nlinks
1021 CALL section_vals_val_get(qmmm_links, "LINK_TYPE", i_rep_section=ikind, &
1022 i_val=link_type)
1023 SELECT CASE (link_type)
1024 CASE (do_qmmm_link_imomm)
1025 num_qm_atom_tot = num_qm_atom_tot + 1
1026 link_involv_mm = link_involv_mm + 1
1027 CASE (do_qmmm_link_pseudo)
1028 num_qm_atom_tot = num_qm_atom_tot + 1
1029 CASE (do_qmmm_link_gho)
1030 ! do nothing for the moment
1031 CASE DEFAULT
1032 cpabort("Unknown QM/MM link type")
1033 END SELECT
1034 END DO
1035 END IF
1036 IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) THEN
1037 ALLOCATE (mm_link_scale_factor(link_involv_mm))
1038 END IF
1039 IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) THEN
1040 ALLOCATE (fist_scale_charge_link(link_involv_mm))
1041 END IF
1042 IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) THEN
1043 ALLOCATE (mm_link_atoms(link_involv_mm))
1044 END IF
1045 IF (PRESENT(qm_atom_index)) ALLOCATE (qm_atom_index(num_qm_atom_tot))
1046 IF (PRESENT(qm_atom_type)) ALLOCATE (qm_atom_type(num_qm_atom_tot))
1047 IF (PRESENT(qm_atom_index)) qm_atom_index = 0
1048 IF (PRESENT(qm_atom_type)) qm_atom_type = " "
1049 num_qm_atom_tot = 1
1050 DO ikind = 1, nkind
1051 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
1052 DO k = 1, n_var
1053 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
1054 i_vals=mm_indexes)
1055 IF (PRESENT(qm_atom_index)) THEN
1056 qm_atom_index(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = mm_indexes(:)
1057 END IF
1058 IF (PRESENT(qm_atom_type)) THEN
1059 CALL section_vals_val_get(qm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
1060 c_val=qm_atom_kind)
1061 qm_atom_type(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = qm_atom_kind
1062 END IF
1063 num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
1064 END DO
1065 END DO
1066 IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) mm_link_scale_factor = 0.0_dp
1067 IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) fist_scale_charge_link = 0.0_dp
1068 IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) mm_link_atoms = 0
1069 IF (explicit) THEN
1070 DO ikind = 1, nlinks
1071 IF (PRESENT(qm_atom_type)) THEN
1072 CALL section_vals_val_get(qmmm_links, "QM_KIND", i_rep_section=ikind, c_val=qm_link_element)
1073 qm_atom_type(num_qm_atom_tot:num_qm_atom_tot) = trim(qm_link_element)//"_LINK"
1074 END IF
1075 IF (PRESENT(qm_atom_index)) THEN
1076 CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1077 cpassert(all(qm_atom_index /= mm_index))
1078 qm_atom_index(num_qm_atom_tot:num_qm_atom_tot) = mm_index
1079 num_qm_atom_tot = num_qm_atom_tot + 1
1080 END IF
1081 IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) THEN
1082 CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1083 mm_link_atoms(ikind) = mm_index
1084 END IF
1085 IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) THEN
1086 CALL section_vals_val_get(qmmm_links, "QMMM_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
1087 mm_link_scale_factor(ikind) = scale_f
1088 END IF
1089 IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) THEN
1090 CALL section_vals_val_get(qmmm_links, "FIST_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
1091 fist_scale_charge_link(ikind) = scale_f
1092 END IF
1093 END DO
1094 END IF
1095 cpassert(num_qm_atom_tot - 1 == SIZE(qm_atom_index))
1096
1097 END SUBROUTINE setup_qm_atom_list
1098
1099! **************************************************************************************************
1100!> \brief this routine sets up all variables to treat qmmm links
1101!> \param qmmm_section ...
1102!> \param qmmm_links ...
1103!> \param mm_el_pot_radius ...
1104!> \param mm_el_pot_radius_corr ...
1105!> \param mm_atom_index ...
1106!> \par History
1107!> 12.2004 created [tlaino]
1108!> \author Teodoro Laino
1109! **************************************************************************************************
1110 SUBROUTINE setup_qmmm_links(qmmm_section, qmmm_links, mm_el_pot_radius, mm_el_pot_radius_corr, &
1111 mm_atom_index)
1112 TYPE(section_vals_type), POINTER :: qmmm_section
1113 TYPE(qmmm_links_type), POINTER :: qmmm_links
1114 REAL(kind=dp), DIMENSION(:), POINTER :: mm_el_pot_radius, mm_el_pot_radius_corr
1115 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1116
1117 INTEGER :: ikind, link_type, mm_index, n_gho, &
1118 n_imomm, n_pseudo, n_rep_val, n_tot, &
1119 nlinks, qm_index
1120 INTEGER, DIMENSION(:), POINTER :: wrk_tmp
1121 REAL(kind=dp) :: alpha, my_radius
1122 TYPE(section_vals_type), POINTER :: qmmm_link_section
1123
1124 NULLIFY (wrk_tmp)
1125 n_imomm = 0
1126 n_gho = 0
1127 n_pseudo = 0
1128 qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
1129 CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
1130 cpassert(nlinks /= 0)
1131 DO ikind = 1, nlinks
1132 CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1133 IF (link_type == do_qmmm_link_imomm) n_imomm = n_imomm + 1
1134 IF (link_type == do_qmmm_link_gho) n_gho = n_gho + 1
1135 IF (link_type == do_qmmm_link_pseudo) n_pseudo = n_pseudo + 1
1136 END DO
1137 n_tot = n_imomm + n_gho + n_pseudo
1138 cpassert(n_tot /= 0)
1139 ALLOCATE (qmmm_links)
1140 NULLIFY (qmmm_links%imomm, &
1141 qmmm_links%pseudo)
1142 ! IMOMM
1143 IF (n_imomm /= 0) THEN
1144 ALLOCATE (qmmm_links%imomm(n_imomm))
1145 ALLOCATE (wrk_tmp(n_imomm))
1146 DO ikind = 1, n_imomm
1147 NULLIFY (qmmm_links%imomm(ikind)%link)
1148 ALLOCATE (qmmm_links%imomm(ikind)%link)
1149 END DO
1150 n_imomm = 0
1151 DO ikind = 1, nlinks
1152 CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1153 IF (link_type == do_qmmm_link_imomm) THEN
1154 n_imomm = n_imomm + 1
1155 CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
1156 CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1157 CALL section_vals_val_get(qmmm_link_section, "ALPHA_IMOMM", i_rep_section=ikind, r_val=alpha)
1158 CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
1159 qmmm_links%imomm(n_imomm)%link%qm_index = qm_index
1160 qmmm_links%imomm(n_imomm)%link%mm_index = mm_index
1161 qmmm_links%imomm(n_imomm)%link%alpha = alpha
1162 wrk_tmp(n_imomm) = mm_index
1163 IF (n_rep_val == 1) THEN
1164 CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, r_val=my_radius)
1165 WHERE (mm_atom_index == mm_index) mm_el_pot_radius = my_radius
1166 WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
1167 END IF
1168 CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
1169 IF (n_rep_val == 1) THEN
1170 CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, r_val=my_radius)
1171 WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
1172 END IF
1173 END IF
1174 END DO
1175 !
1176 ! Checking the link structure
1177 !
1178 DO ikind = 1, SIZE(wrk_tmp)
1179 IF (count(wrk_tmp == wrk_tmp(ikind)) > 1) THEN
1180 CALL cp_abort(__location__, &
1181 "In the IMOMM scheme no more than one QM atom can be bounded to the same "// &
1182 "MM atom. Multiple link MM atom not allowed. Check your link sections.")
1183 END IF
1184 END DO
1185 DEALLOCATE (wrk_tmp)
1186 END IF
1187 ! PSEUDO
1188 IF (n_pseudo /= 0) THEN
1189 ALLOCATE (qmmm_links%pseudo(n_pseudo))
1190 DO ikind = 1, n_pseudo
1191 NULLIFY (qmmm_links%pseudo(ikind)%link)
1192 ALLOCATE (qmmm_links%pseudo(ikind)%link)
1193 END DO
1194 n_pseudo = 0
1195 DO ikind = 1, nlinks
1196 CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1197 IF (link_type == do_qmmm_link_pseudo) THEN
1198 n_pseudo = n_pseudo + 1
1199 CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
1200 CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1201 qmmm_links%pseudo(n_pseudo)%link%qm_index = qm_index
1202 qmmm_links%pseudo(n_pseudo)%link%mm_index = mm_index
1203 END IF
1204 END DO
1205 END IF
1206 ! GHO
1207 IF (n_gho /= 0) THEN
1208 ! not yet implemented
1209 ! still to define : type, implementation into QS
1210 cpabort("QM/MM link with ghost atoms not yet implemented")
1211 END IF
1212 END SUBROUTINE setup_qmmm_links
1213
1214! **************************************************************************************************
1215!> \brief this routine sets up all variables to treat qmmm links
1216!> \param qmmm_section ...
1217!> \param move_mm_charges ...
1218!> \param add_mm_charges ...
1219!> \param mm_atom_chrg ...
1220!> \param mm_el_pot_radius ...
1221!> \param mm_el_pot_radius_corr ...
1222!> \param added_charges ...
1223!> \param mm_atom_index ...
1224!> \par History
1225!> 12.2004 created [tlaino]
1226!> \author Teodoro Laino
1227! **************************************************************************************************
1228 SUBROUTINE move_or_add_atoms(qmmm_section, move_mm_charges, add_mm_charges, &
1229 mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, &
1230 added_charges, mm_atom_index)
1231 TYPE(section_vals_type), POINTER :: qmmm_section
1232 LOGICAL, INTENT(OUT) :: move_mm_charges, add_mm_charges
1233 REAL(kind=dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, &
1234 mm_el_pot_radius_corr
1235 TYPE(add_set_type), POINTER :: added_charges
1236 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1237
1238 INTEGER :: i_add, icount, ikind, ind1, index1, &
1239 index2, n_add_tot, n_adds, n_move_tot, &
1240 n_moves, n_rep_val, nlinks
1241 LOGICAL :: explicit
1242 REAL(kind=dp) :: alpha, c_radius, charge, radius
1243 TYPE(section_vals_type), POINTER :: add_section, move_section, &
1244 qmmm_link_section
1245
1246 explicit = .false.
1247 move_mm_charges = .false.
1248 add_mm_charges = .false.
1249 NULLIFY (qmmm_link_section, move_section, add_section)
1250 qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
1251 CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
1252 cpassert(nlinks /= 0)
1253 icount = 0
1254 n_move_tot = 0
1255 n_add_tot = 0
1256 DO ikind = 1, nlinks
1257 move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
1258 i_rep_section=ikind)
1259 CALL section_vals_get(move_section, n_repetition=n_moves)
1260 add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", &
1261 i_rep_section=ikind)
1262 CALL section_vals_get(add_section, n_repetition=n_adds)
1263 n_move_tot = n_move_tot + n_moves
1264 n_add_tot = n_add_tot + n_adds
1265 END DO
1266 icount = n_move_tot + n_add_tot
1267 IF (n_add_tot /= 0) add_mm_charges = .true.
1268 IF (n_move_tot /= 0) move_mm_charges = .true.
1269 !
1270 ! create add_set_type
1271 !
1272 CALL create_add_set_type(added_charges, ndim=icount)
1273 !
1274 ! Fill in structures
1275 !
1276 icount = 0
1277 DO ikind = 1, nlinks
1278 move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
1279 i_rep_section=ikind)
1280 CALL section_vals_get(move_section, explicit=explicit, n_repetition=n_moves)
1281 !
1282 ! Moving charge atoms
1283 !
1284 IF (explicit) THEN
1285 DO i_add = 1, n_moves
1286 icount = icount + 1
1287 CALL section_vals_val_get(move_section, "ATOM_INDEX_1", i_val=index1, i_rep_section=i_add)
1288 CALL section_vals_val_get(move_section, "ATOM_INDEX_2", i_val=index2, i_rep_section=i_add)
1289 CALL section_vals_val_get(move_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
1290 CALL section_vals_val_get(move_section, "RADIUS", r_val=radius, i_rep_section=i_add)
1291 CALL section_vals_val_get(move_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
1292 c_radius = radius
1293 IF (n_rep_val == 1) THEN
1294 CALL section_vals_val_get(move_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add)
1295 END IF
1296
1297 CALL set_add_set_type(added_charges, icount, index1, index2, alpha, radius, c_radius, &
1298 mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, &
1299 mm_el_pot_radius_corr=mm_el_pot_radius_corr, &
1300 mm_atom_index=mm_atom_index, move=n_moves, ind1=ind1)
1301 END DO
1302 mm_atom_chrg(ind1) = 0.0_dp
1303 END IF
1304
1305 add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", &
1306 i_rep_section=ikind)
1307 CALL section_vals_get(add_section, explicit=explicit, n_repetition=n_adds)
1308 !
1309 ! Adding charge atoms
1310 !
1311 IF (explicit) THEN
1312 DO i_add = 1, n_adds
1313 icount = icount + 1
1314 CALL section_vals_val_get(add_section, "ATOM_INDEX_1", i_val=index1, i_rep_section=i_add)
1315 CALL section_vals_val_get(add_section, "ATOM_INDEX_2", i_val=index2, i_rep_section=i_add)
1316 CALL section_vals_val_get(add_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
1317 CALL section_vals_val_get(add_section, "RADIUS", r_val=radius, i_rep_section=i_add)
1318 CALL section_vals_val_get(add_section, "CHARGE", r_val=charge, i_rep_section=i_add)
1319 CALL section_vals_val_get(add_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
1320 c_radius = radius
1321 IF (n_rep_val == 1) THEN
1322 CALL section_vals_val_get(add_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add)
1323 END IF
1324
1325 CALL set_add_set_type(added_charges, icount, index1, index2, alpha, radius, c_radius, charge, &
1326 mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, &
1327 mm_el_pot_radius_corr=mm_el_pot_radius_corr, &
1328 mm_atom_index=mm_atom_index)
1329 END DO
1330 END IF
1331 END DO
1332
1333 END SUBROUTINE move_or_add_atoms
1334
1335! **************************************************************************************************
1336!> \brief this routine sets up all variables of the add_set_type type
1337!> \param added_charges ...
1338!> \param icount ...
1339!> \param Index1 ...
1340!> \param Index2 ...
1341!> \param alpha ...
1342!> \param radius ...
1343!> \param c_radius ...
1344!> \param charge ...
1345!> \param mm_atom_chrg ...
1346!> \param mm_el_pot_radius ...
1347!> \param mm_el_pot_radius_corr ...
1348!> \param mm_atom_index ...
1349!> \param move ...
1350!> \param ind1 ...
1351!> \par History
1352!> 12.2004 created [tlaino]
1353!> \author Teodoro Laino
1354! **************************************************************************************************
1355 SUBROUTINE set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, charge, &
1356 mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, mm_atom_index, move, ind1)
1357 TYPE(add_set_type), POINTER :: added_charges
1358 INTEGER, INTENT(IN) :: icount, index1, index2
1359 REAL(kind=dp), INTENT(IN) :: alpha, radius, c_radius
1360 REAL(kind=dp), INTENT(IN), OPTIONAL :: charge
1361 REAL(kind=dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, &
1362 mm_el_pot_radius_corr
1363 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1364 INTEGER, INTENT(in), OPTIONAL :: move
1365 INTEGER, INTENT(OUT), OPTIONAL :: ind1
1366
1367 INTEGER :: i, my_move
1368 REAL(kind=dp) :: my_c_radius, my_charge, my_radius
1369
1370 my_move = 0
1371 my_radius = radius
1372 my_c_radius = c_radius
1373 IF (PRESENT(charge)) my_charge = charge
1374 IF (PRESENT(move)) my_move = move
1375 i = 1
1376 getid: DO WHILE (i <= SIZE(mm_atom_index))
1377 IF (index1 == mm_atom_index(i)) EXIT getid
1378 i = i + 1
1379 END DO getid
1380 IF (PRESENT(ind1)) ind1 = i
1381 cpassert(i <= SIZE(mm_atom_index))
1382 IF (.NOT. PRESENT(charge)) my_charge = mm_atom_chrg(i)/real(my_move, kind=dp)
1383 IF (my_radius == 0.0_dp) my_radius = mm_el_pot_radius(i)
1384 IF (my_c_radius == 0.0_dp) my_c_radius = mm_el_pot_radius_corr(i)
1385
1386 added_charges%add_env(icount)%Index1 = index1
1387 added_charges%add_env(icount)%Index2 = index2
1388 added_charges%add_env(icount)%alpha = alpha
1389 added_charges%mm_atom_index(icount) = icount
1390 added_charges%mm_atom_chrg(icount) = my_charge
1391 added_charges%mm_el_pot_radius(icount) = my_radius
1392 added_charges%mm_el_pot_radius_corr(icount) = my_c_radius
1393 END SUBROUTINE set_add_set_type
1394
1395! **************************************************************************************************
1396!> \brief this routine sets up the origin of the MM cell respect to the
1397!> origin of the QM cell. The origin of the QM cell is assumed to be
1398!> in (0.0,0.0,0.0)...
1399!> \param qmmm_section ...
1400!> \param qmmm_env ...
1401!> \param qm_cell_small ...
1402!> \param dr ...
1403!> \par History
1404!> 02.2005 created [tlaino]
1405!> \author Teodoro Laino
1406! **************************************************************************************************
1407 SUBROUTINE setup_origin_mm_cell(qmmm_section, qmmm_env, qm_cell_small, &
1408 dr)
1409 TYPE(section_vals_type), POINTER :: qmmm_section
1410 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1411 TYPE(cell_type), POINTER :: qm_cell_small
1412 REAL(kind=dp), DIMENSION(3), INTENT(in) :: dr
1413
1414 LOGICAL :: center_grid
1415 REAL(kind=dp), DIMENSION(3) :: tmp
1416 REAL(kind=dp), DIMENSION(:), POINTER :: vec
1417
1418! This is the vector that corrects position to apply properly the PBC
1419
1420 tmp(1) = qm_cell_small%hmat(1, 1)
1421 tmp(2) = qm_cell_small%hmat(2, 2)
1422 tmp(3) = qm_cell_small%hmat(3, 3)
1423 cpassert(all(tmp > 0))
1424 qmmm_env%dOmmOqm = tmp/2.0_dp
1425 ! This is unit vector to translate the QM system in order to center it
1426 ! in QM cell
1427 CALL section_vals_val_get(qmmm_section, "CENTER_GRID", l_val=center_grid)
1428 IF (center_grid) THEN
1429 qmmm_env%utrasl = dr
1430 ELSE
1431 qmmm_env%utrasl = 1.0_dp
1432 END IF
1433 CALL section_vals_val_get(qmmm_section, "INITIAL_TRANSLATION_VECTOR", r_vals=vec)
1434 qmmm_env%transl_v = vec
1435 END SUBROUTINE setup_origin_mm_cell
1436
1437! **************************************************************************************************
1438!> \brief this routine sets up list of MM atoms carrying an image charge
1439!> \param image_charge_section ...
1440!> \param qmmm_env ...
1441!> \param qm_atom_index ...
1442!> \param subsys_mm ...
1443!> \par History
1444!> 02.2012 created
1445!> \author Dorothea Golze
1446! **************************************************************************************************
1447 SUBROUTINE setup_image_atom_list(image_charge_section, qmmm_env, &
1448 qm_atom_index, subsys_mm)
1449
1450 TYPE(section_vals_type), POINTER :: image_charge_section
1451 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1452 INTEGER, DIMENSION(:), POINTER :: qm_atom_index
1453 TYPE(cp_subsys_type), POINTER :: subsys_mm
1454
1455 INTEGER :: atom_a, atom_b, i, j, k, max_index, &
1456 n_var, num_const_atom, &
1457 num_image_mm_atom
1458 INTEGER, DIMENSION(:), POINTER :: mm_indexes
1459 LOGICAL :: fix_xyz, imageind_in_range
1460 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind
1461
1462 NULLIFY (mm_indexes, molecule_kind)
1463 imageind_in_range = .false.
1464 num_image_mm_atom = 0
1465 max_index = 0
1466
1467 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1468 n_rep_val=n_var)
1469 DO i = 1, n_var
1470 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1471 i_rep_val=i, i_vals=mm_indexes)
1472 num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
1473 END DO
1474
1475 ALLOCATE (qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom))
1476
1477 qmmm_env%image_charge_pot%image_mm_list = 0
1478 num_image_mm_atom = 1
1479
1480 DO i = 1, n_var
1481 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1482 i_rep_val=i, i_vals=mm_indexes)
1483 qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom:num_image_mm_atom &
1484 + SIZE(mm_indexes) - 1) = mm_indexes(:)
1485 num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
1486 END DO
1487
1488 ! checking, if in range, if list contains QM atoms or any atoms doubled
1489 num_image_mm_atom = num_image_mm_atom - 1
1490
1491 max_index = SIZE(subsys_mm%particles%els)
1492
1493 cpassert(SIZE(qmmm_env%image_charge_pot%image_mm_list) /= 0)
1494 imageind_in_range = (maxval(qmmm_env%image_charge_pot%image_mm_list) <= max_index) &
1495 .AND. (minval(qmmm_env%image_charge_pot%image_mm_list) > 0)
1496 cpassert(imageind_in_range)
1497
1498 DO i = 1, num_image_mm_atom
1499 atom_a = qmmm_env%image_charge_pot%image_mm_list(i)
1500 IF (any(qm_atom_index == atom_a)) THEN
1501 cpabort("Image atom list must only contain MM atoms")
1502 END IF
1503 DO j = i + 1, num_image_mm_atom
1504 atom_b = qmmm_env%image_charge_pot%image_mm_list(j)
1505 IF (atom_a == atom_b) THEN
1506 cpabort("There are atoms doubled in image list.")
1507 END IF
1508 END DO
1509 END DO
1510
1511 ! check if molecules in list carry constraints
1512 num_const_atom = 0
1513 fix_xyz = .true.
1514 IF (ASSOCIATED(subsys_mm%molecule_kinds)) THEN
1515 IF (ASSOCIATED(subsys_mm%molecule_kinds%els)) THEN
1516 molecule_kind => subsys_mm%molecule_kinds%els
1517 DO i = 1, SIZE(molecule_kind)
1518 IF (.NOT. ASSOCIATED(molecule_kind(i)%fixd_list)) EXIT
1519 IF (.NOT. fix_xyz) EXIT
1520 DO j = 1, SIZE(molecule_kind(i)%fixd_list)
1521 IF (.NOT. fix_xyz) EXIT
1522 DO k = 1, num_image_mm_atom
1523 atom_a = qmmm_env%image_charge_pot%image_mm_list(k)
1524 IF (atom_a == molecule_kind(i)%fixd_list(j)%fixd) THEN
1525 num_const_atom = num_const_atom + 1
1526 IF (molecule_kind(i)%fixd_list(j)%itype /= use_perd_xyz) THEN
1527 fix_xyz = .false.
1528 EXIT
1529 END IF
1530 END IF
1531 END DO
1532 END DO
1533 END DO
1534 END IF
1535 END IF
1536
1537 ! if all image atoms are constrained, calculate image matrix only
1538 ! once for the first MD or GEO_OPT step (for non-iterative case)
1539 IF (num_const_atom == num_image_mm_atom .AND. fix_xyz) THEN
1540 qmmm_env%image_charge_pot%state_image_matrix = calc_once
1541 ELSE
1542 qmmm_env%image_charge_pot%state_image_matrix = calc_always
1543 END IF
1544
1545 END SUBROUTINE setup_image_atom_list
1546
1547! **************************************************************************************************
1548!> \brief Print info on image charges
1549!> \param qmmm_env ...
1550!> \param qmmm_section ...
1551!> \par History
1552!> 03.2012 created
1553!> \author Dorothea Golze
1554! **************************************************************************************************
1555 SUBROUTINE print_image_charge_info(qmmm_env, qmmm_section)
1556
1557 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1558 TYPE(section_vals_type), POINTER :: qmmm_section
1559
1560 INTEGER :: iw
1561 REAL(kind=dp) :: eta, eta_conv, v0, v0_conv
1562 TYPE(cp_logger_type), POINTER :: logger
1563
1564 logger => cp_get_default_logger()
1565 iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%PROGRAM_RUN_INFO", &
1566 extension=".log")
1567 eta = qmmm_env%image_charge_pot%eta
1568 eta_conv = cp_unit_from_cp2k(eta, "angstrom", power=-2)
1569 v0 = qmmm_env%image_charge_pot%V0
1570 v0_conv = cp_unit_from_cp2k(v0, "volt")
1571
1572 IF (iw > 0) THEN
1573 WRITE (iw, fmt="(T25,A)") "IMAGE CHARGE PARAMETERS"
1574 WRITE (iw, fmt="(T25,A)") repeat("-", 23)
1575 WRITE (iw, fmt="(/)")
1576 WRITE (iw, fmt="(T2,A)") "INDEX OF MM ATOMS CARRYING AN IMAGE CHARGE:"
1577 WRITE (iw, fmt="(/)")
1578
1579 WRITE (iw, "(7X,10I6)") qmmm_env%image_charge_pot%image_mm_list
1580 WRITE (iw, fmt="(/)")
1581 WRITE (iw, "(T2,A52,T69,F12.8)") &
1582 "WIDTH OF GAUSSIAN CHARGE DISTRIBUTION [angstrom^-2]:", eta_conv
1583 WRITE (iw, "(T2,A26,T69,F12.8)") "EXTERNAL POTENTIAL [volt]:", v0_conv
1584 WRITE (iw, fmt="(/,T2,A,/)") repeat("-", 79)
1585 END IF
1586 CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
1587 "PRINT%PROGRAM_RUN_INFO")
1588
1589 END SUBROUTINE print_image_charge_info
1590
1591END 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_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, cell_ref, use_ref_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:1251
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:1222
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, nequip, bmhft, bmhftd, ipbv, buck4r, buckmo, gp, tersoff, siepmann, gal, gal21, tab, deepmd, ace)
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
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 ...
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:629
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:367
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:1231
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:114
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:1409
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:522
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:211
subroutine, public print_qmmm_links(qmmm_section, qmmm_links)
Print info on qm/mm links.
Definition qmmm_init.F:306
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:847
subroutine, public print_image_charge_info(qmmm_env, qmmm_section)
Print info on image charges.
Definition qmmm_init.F:1556
subroutine, public setup_qmmm_links(qmmm_section, qmmm_links, mm_el_pot_radius, mm_el_pot_radius_corr, mm_atom_index)
this routine sets up all variables to treat qmmm links
Definition qmmm_init.F:1112
subroutine, public qmmm_init_potential(qmmm_env_qm, mm_cell, added_charges, added_shells, print_section)
...
Definition qmmm_init.F:466
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, 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, mimic, 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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
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)
...
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represents a system: atoms, molecules, their pos,vel,...
stores all the informations relevant to an mpi environment
contained for different pw related things
parameters for core-shell model potentials