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