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