(git:374b731)
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-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Initialize a QM/MM calculation
10!> \par History
11!> 5.2004 created [fawzi]
12!> \author Fawzi Mohamed
13! **************************************************************************************************
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 para_env=para_env, &
563 compatibility=qmmm_env_qm%compatibility, &
564 qmmm_periodic=qmmm_periodic, &
565 print_section=print_section, &
566 eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
567 maxchrg=maxchrg, &
568 ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, &
569 ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
570
571 IF (qmmm_env_qm%move_mm_charges .OR. qmmm_env_qm%add_mm_charges) THEN
572
573 CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
574 per_potentials=added_charges%per_potentials, &
575 potentials=added_charges%potentials, &
576 pgfs=added_charges%pgfs, &
577 qm_cell_small=qm_cell_small, &
578 mm_cell=mm_cell, &
579 para_env=para_env, &
580 compatibility=qmmm_env_qm%compatibility, &
581 qmmm_periodic=qmmm_periodic, &
582 print_section=print_section, &
583 eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
584 maxchrg=maxchrg, &
585 ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, &
586 ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
587 END IF
588
589 IF (qmmm_env_qm%added_shells%num_mm_atoms .GT. 0) THEN
590
591 CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
592 per_potentials=added_shells%per_potentials, &
593 potentials=added_shells%potentials, &
594 pgfs=added_shells%pgfs, &
595 qm_cell_small=qm_cell_small, &
596 mm_cell=mm_cell, &
597 para_env=para_env, &
598 compatibility=qmmm_env_qm%compatibility, &
599 qmmm_periodic=qmmm_periodic, &
600 print_section=print_section, &
601 eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
602 maxchrg=maxchrg, &
603 ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, &
604 ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
605 END IF
606
607 END IF
608
609 END IF
610
611 END SUBROUTINE qmmm_init_periodic_potential
612
613! **************************************************************************************************
614!> \brief ...
615!> \param qmmm_section ...
616!> \param qmmm_env ...
617!> \param subsys_mm ...
618!> \param qm_atom_type ...
619!> \param qm_atom_index ...
620!> \param mm_atom_index ...
621!> \param qm_cell_small ...
622!> \param qmmm_coupl_type ...
623!> \param eps_mm_rspace ...
624!> \param qmmm_link ...
625!> \param para_env ...
626!> \par History
627!> 11.2004 created [tlaino]
628!> \author Teodoro Laino
629! **************************************************************************************************
630 SUBROUTINE setup_qmmm_vars_qm(qmmm_section, qmmm_env, subsys_mm, qm_atom_type, &
631 qm_atom_index, mm_atom_index, qm_cell_small, qmmm_coupl_type, eps_mm_rspace, &
632 qmmm_link, para_env)
633 TYPE(section_vals_type), POINTER :: qmmm_section
634 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
635 TYPE(cp_subsys_type), POINTER :: subsys_mm
636 CHARACTER(len=default_string_length), &
637 DIMENSION(:), POINTER :: qm_atom_type
638 INTEGER, DIMENSION(:), POINTER :: qm_atom_index, mm_atom_index
639 TYPE(cell_type), POINTER :: qm_cell_small
640 INTEGER, INTENT(OUT) :: qmmm_coupl_type
641 REAL(kind=dp), INTENT(OUT) :: eps_mm_rspace
642 LOGICAL, INTENT(OUT) :: qmmm_link
643 TYPE(mp_para_env_type), POINTER :: para_env
644
645 CHARACTER(len=default_string_length) :: atmname, mm_atom_kind
646 INTEGER :: i, icount, ikind, ikindr, my_type, &
647 n_rep_val, nkind, size_mm_system
648 INTEGER, DIMENSION(:), POINTER :: mm_link_atoms
649 LOGICAL :: explicit, is_mm, is_qm
650 REAL(kind=dp) :: tmp_radius, tmp_radius_c
651 REAL(kind=dp), DIMENSION(:), POINTER :: tmp_sph_cut
652 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
653 TYPE(atomic_kind_type), POINTER :: atomic_kind
654 TYPE(fist_potential_type), POINTER :: fist_potential
655 TYPE(section_vals_type), POINTER :: cell_section, eri_mme_section, &
656 image_charge_section, mm_kinds
657
658 NULLIFY (mm_link_atoms, cell_section, tmp_sph_cut)
659 NULLIFY (image_charge_section)
660 qmmm_link = .false.
661
662 CALL section_vals_get(qmmm_section, explicit=explicit)
663 IF (explicit) THEN
664 !
665 ! QM_CELL
666 !
667 cell_section => section_vals_get_subs_vals(qmmm_section, "CELL")
668 CALL read_cell(qm_cell_small, qm_cell_small, cell_section=cell_section, &
669 check_for_ref=.false., para_env=para_env)
670 qm_cell_small%tag = "CELL_QM"
671 CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type)
672 CALL section_vals_val_get(qmmm_section, "EPS_MM_RSPACE", r_val=eps_mm_rspace)
673 CALL section_vals_val_get(qmmm_section, "SPHERICAL_CUTOFF", r_vals=tmp_sph_cut)
674 cpassert(SIZE(tmp_sph_cut) == 2)
675 qmmm_env%spherical_cutoff = tmp_sph_cut
676 IF (qmmm_env%spherical_cutoff(1) <= 0.0_dp) THEN
677 qmmm_env%spherical_cutoff(2) = 0.0_dp
678 ELSE
679 IF (qmmm_env%spherical_cutoff(2) <= 0.0_dp) qmmm_env%spherical_cutoff(2) = epsilon(0.0_dp)
680 tmp_radius = qmmm_env%spherical_cutoff(1) - 20.0_dp*qmmm_env%spherical_cutoff(2)
681 IF (tmp_radius <= 0.0_dp) &
682 CALL cp_abort(__location__, &
683 "SPHERICAL_CUTOFF(1) > 20*SPHERICAL_CUTOFF(1)! Please correct parameters for "// &
684 "the Spherical Cutoff in order to satisfy the previous condition!")
685 END IF
686 !
687 ! Initialization of arrays and core_charge_radius...
688 !
689 tmp_radius = 0.0_dp
690 CALL cp_subsys_get(subsys=subsys_mm, atomic_kinds=atomic_kinds)
691 DO ikind = 1, SIZE(atomic_kinds%els)
692 atomic_kind => atomic_kinds%els(ikind)
693 CALL get_atomic_kind(atomic_kind=atomic_kind, &
694 fist_potential=fist_potential)
695 CALL set_potential(potential=fist_potential, &
696 qmmm_radius=tmp_radius, &
697 qmmm_corr_radius=tmp_radius)
698 CALL set_atomic_kind(atomic_kind=atomic_kind, &
699 fist_potential=fist_potential)
700 END DO
701 CALL setup_qm_atom_list(qmmm_section=qmmm_section, &
702 qm_atom_index=qm_atom_index, &
703 qm_atom_type=qm_atom_type, &
704 mm_link_atoms=mm_link_atoms, &
705 qmmm_link=qmmm_link)
706 !
707 ! MM_KINDS
708 !
709 mm_kinds => section_vals_get_subs_vals(qmmm_section, "MM_KIND")
710 CALL section_vals_get(mm_kinds, explicit=explicit, n_repetition=nkind)
711 !
712 ! Default
713 !
714 tmp_radius = cp_unit_to_cp2k(radius_qmmm_default, "angstrom")
715 set_radius_pot_0: DO ikindr = 1, SIZE(atomic_kinds%els)
716 atomic_kind => atomic_kinds%els(ikindr)
717 CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
718 CALL get_atomic_kind(atomic_kind=atomic_kind, &
719 fist_potential=fist_potential)
720 CALL set_potential(potential=fist_potential, qmmm_radius=tmp_radius, &
721 qmmm_corr_radius=tmp_radius)
722 CALL set_atomic_kind(atomic_kind=atomic_kind, &
723 fist_potential=fist_potential)
724 END DO set_radius_pot_0
725 !
726 ! If present overwrite the kind specified in input file...
727 !
728 IF (explicit) THEN
729 DO ikind = 1, nkind
730 CALL section_vals_val_get(mm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
731 c_val=mm_atom_kind)
732 CALL section_vals_val_get(mm_kinds, "RADIUS", i_rep_section=ikind, r_val=tmp_radius)
733 tmp_radius_c = tmp_radius
734 CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
735 IF (n_rep_val == 1) CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, &
736 r_val=tmp_radius_c)
737 set_radius_pot_1: DO ikindr = 1, SIZE(atomic_kinds%els)
738 atomic_kind => atomic_kinds%els(ikindr)
739 CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
740 is_qm = qmmm_ff_precond_only_qm(atmname)
741 IF (trim(mm_atom_kind) == atmname) THEN
742 CALL get_atomic_kind(atomic_kind=atomic_kind, &
743 fist_potential=fist_potential)
744 CALL set_potential(potential=fist_potential, &
745 qmmm_radius=tmp_radius, &
746 qmmm_corr_radius=tmp_radius_c)
747 CALL set_atomic_kind(atomic_kind=atomic_kind, &
748 fist_potential=fist_potential)
749 END IF
750 END DO set_radius_pot_1
751 END DO
752 END IF
753
754 !Image charge section
755
756 image_charge_section => section_vals_get_subs_vals(qmmm_section, "IMAGE_CHARGE")
757 CALL section_vals_get(image_charge_section, explicit=qmmm_env%image_charge)
758
759 ELSE
760 cpabort("QMMM section not present in input file!")
761 END IF
762 !
763 ! Build MM atoms list
764 !
765 size_mm_system = SIZE(subsys_mm%particles%els) - SIZE(qm_atom_index)
766 IF (qmmm_link .AND. ASSOCIATED(mm_link_atoms)) size_mm_system = size_mm_system + SIZE(mm_link_atoms)
767 ALLOCATE (mm_atom_index(size_mm_system))
768 icount = 0
769
770 DO i = 1, SIZE(subsys_mm%particles%els)
771 is_mm = .true.
772 IF (any(qm_atom_index == i)) THEN
773 is_mm = .false.
774 END IF
775 IF (ASSOCIATED(mm_link_atoms)) THEN
776 IF (any(mm_link_atoms == i) .AND. qmmm_link) is_mm = .true.
777 END IF
778 IF (is_mm) THEN
779 icount = icount + 1
780 IF (icount <= size_mm_system) mm_atom_index(icount) = i
781 END IF
782 END DO
783 cpassert(icount == size_mm_system)
784 IF (ASSOCIATED(mm_link_atoms)) THEN
785 DEALLOCATE (mm_link_atoms)
786 END IF
787
788 ! Build image charge atom list + set up variables
789 !
790 IF (qmmm_env%image_charge) THEN
791 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
792 explicit=explicit)
793 IF (explicit) qmmm_env%image_charge_pot%all_mm = .false.
794
795 IF (qmmm_env%image_charge_pot%all_mm) THEN
796 qmmm_env%image_charge_pot%image_mm_list => mm_atom_index
797 ELSE
798 CALL setup_image_atom_list(image_charge_section, qmmm_env, &
799 qm_atom_index, subsys_mm)
800 END IF
801
802 qmmm_env%image_charge_pot%particles_all => subsys_mm%particles%els
803
804 CALL section_vals_val_get(image_charge_section, "EXT_POTENTIAL", &
805 r_val=qmmm_env%image_charge_pot%V0)
806 CALL section_vals_val_get(image_charge_section, "WIDTH", &
807 r_val=qmmm_env%image_charge_pot%eta)
808 CALL section_vals_val_get(image_charge_section, "DETERM_COEFF", &
809 i_val=my_type)
810 SELECT CASE (my_type)
812 qmmm_env%image_charge_pot%coeff_iterative = .false.
813 CASE (do_qmmm_image_iter)
814 qmmm_env%image_charge_pot%coeff_iterative = .true.
815 END SELECT
816
817 CALL section_vals_val_get(image_charge_section, "RESTART_IMAGE_MATRIX", &
818 l_val=qmmm_env%image_charge_pot%image_restart)
819
820 CALL section_vals_val_get(image_charge_section, "IMAGE_MATRIX_METHOD", &
821 i_val=qmmm_env%image_charge_pot%image_matrix_method)
822
823 IF (qmmm_env%image_charge_pot%image_matrix_method .EQ. do_eri_mme) THEN
824 eri_mme_section => section_vals_get_subs_vals(image_charge_section, "ERI_MME")
825 CALL cp_eri_mme_init_read_input(eri_mme_section, qmmm_env%image_charge_pot%eri_mme_param)
826 CALL cp_eri_mme_set_params(qmmm_env%image_charge_pot%eri_mme_param, &
827 hmat=qm_cell_small%hmat, is_ortho=qm_cell_small%orthorhombic, &
828 zet_min=qmmm_env%image_charge_pot%eta, &
829 zet_max=qmmm_env%image_charge_pot%eta, &
830 l_max_zet=0, &
831 l_max=0, &
832 para_env=para_env)
833
834 END IF
835 END IF
836
837 END SUBROUTINE setup_qmmm_vars_qm
838
839! **************************************************************************************************
840!> \brief ...
841!> \param qmmm_section ...
842!> \param qmmm_env ...
843!> \param qm_atom_index ...
844!> \param mm_link_atoms ...
845!> \param mm_link_scale_factor ...
846!> \param fist_scale_charge_link ...
847!> \param qmmm_coupl_type ...
848!> \param qmmm_link ...
849!> \par History
850!> 12.2004 created [tlaino]
851!> \author Teodoro Laino
852! **************************************************************************************************
853 SUBROUTINE setup_qmmm_vars_mm(qmmm_section, qmmm_env, qm_atom_index, &
854 mm_link_atoms, mm_link_scale_factor, &
855 fist_scale_charge_link, qmmm_coupl_type, &
856 qmmm_link)
857 TYPE(section_vals_type), POINTER :: qmmm_section
858 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
859 INTEGER, DIMENSION(:), POINTER :: qm_atom_index, mm_link_atoms
860 REAL(kind=dp), DIMENSION(:), POINTER :: mm_link_scale_factor, &
861 fist_scale_charge_link
862 INTEGER, INTENT(OUT) :: qmmm_coupl_type
863 LOGICAL, INTENT(OUT) :: qmmm_link
864
865 LOGICAL :: explicit
866 TYPE(section_vals_type), POINTER :: qmmm_ff_section
867
868 NULLIFY (qmmm_ff_section)
869 qmmm_link = .false.
870 CALL section_vals_get(qmmm_section, explicit=explicit)
871 IF (explicit) THEN
872 CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type)
873 CALL setup_qm_atom_list(qmmm_section, qm_atom_index=qm_atom_index, qmmm_link=qmmm_link, &
874 mm_link_atoms=mm_link_atoms, mm_link_scale_factor=mm_link_scale_factor, &
875 fist_scale_charge_link=fist_scale_charge_link)
876 !
877 ! Do we want to use a different FF for the non-bonded QM/MM interactions?
878 !
879 qmmm_ff_section => section_vals_get_subs_vals(qmmm_section, "FORCEFIELD")
880 CALL section_vals_get(qmmm_ff_section, explicit=qmmm_env%use_qmmm_ff)
881 IF (qmmm_env%use_qmmm_ff) THEN
882 CALL section_vals_val_get(qmmm_ff_section, "MULTIPLE_POTENTIAL", &
883 l_val=qmmm_env%multiple_potential)
884 CALL read_qmmm_ff_section(qmmm_ff_section, qmmm_env%inp_info)
885 END IF
886 END IF
887 END SUBROUTINE setup_qmmm_vars_mm
888
889! **************************************************************************************************
890!> \brief reads information regarding the forcefield specific for the QM/MM
891!> interactions
892!> \param qmmm_ff_section ...
893!> \param inp_info ...
894!> \par History
895!> 12.2004 created [tlaino]
896!> \author Teodoro Laino
897! **************************************************************************************************
898 SUBROUTINE read_qmmm_ff_section(qmmm_ff_section, inp_info)
899 TYPE(section_vals_type), POINTER :: qmmm_ff_section
900 TYPE(input_info_type), POINTER :: inp_info
901
902 INTEGER :: n_gd, n_gp, n_lj, n_wl, np
903 TYPE(section_vals_type), POINTER :: gd_section, gp_section, lj_section, &
904 wl_section
905
906!
907! NONBONDED
908!
909
910 lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%LENNARD-JONES")
911 wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%WILLIAMS")
912 gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GOODWIN")
913 gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GENPOT")
914 CALL section_vals_get(lj_section, n_repetition=n_lj)
915 np = n_lj
916 IF (n_lj /= 0) THEN
917 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, lj_charmm=.true.)
918 CALL read_lj_section(inp_info%nonbonded, lj_section, start=0)
919 END IF
920 CALL section_vals_get(wl_section, n_repetition=n_wl)
921 np = n_lj + n_wl
922 IF (n_wl /= 0) THEN
923 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, williams=.true.)
924 CALL read_wl_section(inp_info%nonbonded, wl_section, start=n_lj)
925 END IF
926 CALL section_vals_get(gd_section, n_repetition=n_gd)
927 np = n_lj + n_wl + n_gd
928 IF (n_gd /= 0) THEN
929 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, goodwin=.true.)
930 CALL read_gd_section(inp_info%nonbonded, gd_section, start=n_lj + n_wl)
931 END IF
932 CALL section_vals_get(gp_section, n_repetition=n_gp)
933 np = n_lj + n_wl + n_gd + n_gp
934 IF (n_gp /= 0) THEN
935 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, gp=.true.)
936 CALL read_gp_section(inp_info%nonbonded, gp_section, start=n_lj + n_wl + n_gd)
937 END IF
938 !
939 ! NONBONDED14
940 !
941 lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%LENNARD-JONES")
942 wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%WILLIAMS")
943 gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GOODWIN")
944 gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GENPOT")
945 CALL section_vals_get(lj_section, n_repetition=n_lj)
946 np = n_lj
947 IF (n_lj /= 0) THEN
948 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, lj_charmm=.true.)
949 CALL read_lj_section(inp_info%nonbonded14, lj_section, start=0)
950 END IF
951 CALL section_vals_get(wl_section, n_repetition=n_wl)
952 np = n_lj + n_wl
953 IF (n_wl /= 0) THEN
954 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, williams=.true.)
955 CALL read_wl_section(inp_info%nonbonded14, wl_section, start=n_lj)
956 END IF
957 CALL section_vals_get(gd_section, n_repetition=n_gd)
958 np = n_lj + n_wl + n_gd
959 IF (n_gd /= 0) THEN
960 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, goodwin=.true.)
961 CALL read_gd_section(inp_info%nonbonded14, gd_section, start=n_lj + n_wl)
962 END IF
963 CALL section_vals_get(gp_section, n_repetition=n_gp)
964 np = n_lj + n_wl + n_gd + n_gp
965 IF (n_gp /= 0) THEN
966 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, gp=.true.)
967 CALL read_gp_section(inp_info%nonbonded14, gp_section, start=n_lj + n_wl + n_gd)
968 END IF
969
970 END SUBROUTINE read_qmmm_ff_section
971
972! **************************************************************************************************
973!> \brief ...
974!> \param qmmm_section ...
975!> \param qm_atom_index ...
976!> \param qm_atom_type ...
977!> \param mm_link_atoms ...
978!> \param mm_link_scale_factor ...
979!> \param qmmm_link ...
980!> \param fist_scale_charge_link ...
981!> \par History
982!> 12.2004 created [tlaino]
983!> \author Teodoro Laino
984! **************************************************************************************************
985 SUBROUTINE setup_qm_atom_list(qmmm_section, qm_atom_index, qm_atom_type, &
986 mm_link_atoms, mm_link_scale_factor, qmmm_link, fist_scale_charge_link)
987 TYPE(section_vals_type), POINTER :: qmmm_section
988 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: qm_atom_index
989 CHARACTER(len=default_string_length), &
990 DIMENSION(:), OPTIONAL, POINTER :: qm_atom_type
991 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: mm_link_atoms
992 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: mm_link_scale_factor
993 LOGICAL, INTENT(OUT), OPTIONAL :: qmmm_link
994 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: fist_scale_charge_link
995
996 CHARACTER(len=default_string_length) :: qm_atom_kind, qm_link_element
997 INTEGER :: ikind, k, link_involv_mm, link_type, &
998 mm_index, n_var, nkind, nlinks, &
999 num_qm_atom_tot
1000 INTEGER, DIMENSION(:), POINTER :: mm_indexes
1001 LOGICAL :: explicit
1002 REAL(kind=dp) :: scale_f
1003 TYPE(section_vals_type), POINTER :: qm_kinds, qmmm_links
1004
1005 num_qm_atom_tot = 0
1006 link_involv_mm = 0
1007 nlinks = 0
1008 !
1009 ! QM_KINDS
1010 !
1011 qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND")
1012 CALL section_vals_get(qm_kinds, n_repetition=nkind)
1013 DO ikind = 1, nkind
1014 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
1015 DO k = 1, n_var
1016 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
1017 i_vals=mm_indexes)
1018 num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
1019 END DO
1020 END DO
1021 !
1022 ! QM/MM LINKS
1023 !
1024 qmmm_links => section_vals_get_subs_vals(qmmm_section, "LINK")
1025 CALL section_vals_get(qmmm_links, explicit=explicit)
1026 IF (explicit) THEN
1027 qmmm_link = .true.
1028 CALL section_vals_get(qmmm_links, n_repetition=nlinks)
1029 ! Take care of the various link types
1030 DO ikind = 1, nlinks
1031 CALL section_vals_val_get(qmmm_links, "LINK_TYPE", i_rep_section=ikind, &
1032 i_val=link_type)
1033 SELECT CASE (link_type)
1034 CASE (do_qmmm_link_imomm)
1035 num_qm_atom_tot = num_qm_atom_tot + 1
1036 link_involv_mm = link_involv_mm + 1
1037 CASE (do_qmmm_link_pseudo)
1038 num_qm_atom_tot = num_qm_atom_tot + 1
1039 CASE (do_qmmm_link_gho)
1040 ! do nothing for the moment
1041 CASE DEFAULT
1042 cpabort("")
1043 END SELECT
1044 END DO
1045 END IF
1046 IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) &
1047 ALLOCATE (mm_link_scale_factor(link_involv_mm))
1048 IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) &
1049 ALLOCATE (fist_scale_charge_link(link_involv_mm))
1050 IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) &
1051 ALLOCATE (mm_link_atoms(link_involv_mm))
1052 IF (PRESENT(qm_atom_index)) ALLOCATE (qm_atom_index(num_qm_atom_tot))
1053 IF (PRESENT(qm_atom_type)) ALLOCATE (qm_atom_type(num_qm_atom_tot))
1054 IF (PRESENT(qm_atom_index)) qm_atom_index = 0
1055 IF (PRESENT(qm_atom_type)) qm_atom_type = " "
1056 num_qm_atom_tot = 1
1057 DO ikind = 1, nkind
1058 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
1059 DO k = 1, n_var
1060 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
1061 i_vals=mm_indexes)
1062 IF (PRESENT(qm_atom_index)) THEN
1063 qm_atom_index(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = mm_indexes(:)
1064 END IF
1065 IF (PRESENT(qm_atom_type)) THEN
1066 CALL section_vals_val_get(qm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
1067 c_val=qm_atom_kind)
1068 qm_atom_type(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = qm_atom_kind
1069 END IF
1070 num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
1071 END DO
1072 END DO
1073 IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) mm_link_scale_factor = 0.0_dp
1074 IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) fist_scale_charge_link = 0.0_dp
1075 IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) mm_link_atoms = 0
1076 IF (explicit) THEN
1077 DO ikind = 1, nlinks
1078 IF (PRESENT(qm_atom_type)) THEN
1079 CALL section_vals_val_get(qmmm_links, "QM_KIND", i_rep_section=ikind, c_val=qm_link_element)
1080 qm_atom_type(num_qm_atom_tot:num_qm_atom_tot) = trim(qm_link_element)//"_LINK"
1081 END IF
1082 IF (PRESENT(qm_atom_index)) THEN
1083 CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1084 cpassert(all(qm_atom_index /= mm_index))
1085 qm_atom_index(num_qm_atom_tot:num_qm_atom_tot) = mm_index
1086 num_qm_atom_tot = num_qm_atom_tot + 1
1087 END IF
1088 IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) THEN
1089 CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1090 mm_link_atoms(ikind) = mm_index
1091 END IF
1092 IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) THEN
1093 CALL section_vals_val_get(qmmm_links, "QMMM_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
1094 mm_link_scale_factor(ikind) = scale_f
1095 END IF
1096 IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) THEN
1097 CALL section_vals_val_get(qmmm_links, "FIST_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
1098 fist_scale_charge_link(ikind) = scale_f
1099 END IF
1100 END DO
1101 END IF
1102 cpassert(num_qm_atom_tot - 1 == SIZE(qm_atom_index))
1103
1104 END SUBROUTINE setup_qm_atom_list
1105
1106! **************************************************************************************************
1107!> \brief this routine sets up all variables to treat qmmm links
1108!> \param qmmm_section ...
1109!> \param qmmm_links ...
1110!> \param mm_el_pot_radius ...
1111!> \param mm_el_pot_radius_corr ...
1112!> \param mm_atom_index ...
1113!> \param iw ...
1114!> \par History
1115!> 12.2004 created [tlaino]
1116!> \author Teodoro Laino
1117! **************************************************************************************************
1118 SUBROUTINE setup_qmmm_links(qmmm_section, qmmm_links, mm_el_pot_radius, mm_el_pot_radius_corr, &
1119 mm_atom_index, iw)
1120 TYPE(section_vals_type), POINTER :: qmmm_section
1121 TYPE(qmmm_links_type), POINTER :: qmmm_links
1122 REAL(kind=dp), DIMENSION(:), POINTER :: mm_el_pot_radius, mm_el_pot_radius_corr
1123 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1124 INTEGER, INTENT(IN) :: iw
1125
1126 INTEGER :: ikind, link_type, mm_index, n_gho, &
1127 n_imomm, n_pseudo, n_rep_val, n_tot, &
1128 nlinks, qm_index
1129 INTEGER, DIMENSION(:), POINTER :: wrk_tmp
1130 REAL(kind=dp) :: alpha, my_radius
1131 TYPE(section_vals_type), POINTER :: qmmm_link_section
1132
1133 NULLIFY (wrk_tmp)
1134 n_imomm = 0
1135 n_gho = 0
1136 n_pseudo = 0
1137 qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
1138 CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
1139 cpassert(nlinks /= 0)
1140 DO ikind = 1, nlinks
1141 CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1142 IF (link_type == do_qmmm_link_imomm) n_imomm = n_imomm + 1
1143 IF (link_type == do_qmmm_link_gho) n_gho = n_gho + 1
1144 IF (link_type == do_qmmm_link_pseudo) n_pseudo = n_pseudo + 1
1145 END DO
1146 n_tot = n_imomm + n_gho + n_pseudo
1147 cpassert(n_tot /= 0)
1148 ALLOCATE (qmmm_links)
1149 NULLIFY (qmmm_links%imomm, &
1150 qmmm_links%pseudo)
1151 ! IMOMM
1152 IF (n_imomm /= 0) THEN
1153 ALLOCATE (qmmm_links%imomm(n_imomm))
1154 ALLOCATE (wrk_tmp(n_imomm))
1155 DO ikind = 1, n_imomm
1156 NULLIFY (qmmm_links%imomm(ikind)%link)
1157 ALLOCATE (qmmm_links%imomm(ikind)%link)
1158 END DO
1159 n_imomm = 0
1160 DO ikind = 1, nlinks
1161 CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1162 IF (link_type == do_qmmm_link_imomm) THEN
1163 n_imomm = n_imomm + 1
1164 CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
1165 CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1166 CALL section_vals_val_get(qmmm_link_section, "ALPHA_IMOMM", i_rep_section=ikind, r_val=alpha)
1167 CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
1168 qmmm_links%imomm(n_imomm)%link%qm_index = qm_index
1169 qmmm_links%imomm(n_imomm)%link%mm_index = mm_index
1170 qmmm_links%imomm(n_imomm)%link%alpha = alpha
1171 wrk_tmp(n_imomm) = mm_index
1172 IF (n_rep_val == 1) THEN
1173 CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, r_val=my_radius)
1174 WHERE (mm_atom_index == mm_index) mm_el_pot_radius = my_radius
1175 WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
1176 END IF
1177 CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
1178 IF (n_rep_val == 1) THEN
1179 CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, r_val=my_radius)
1180 WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
1181 END IF
1182 END IF
1183 END DO
1184 !
1185 ! Checking the link structure
1186 !
1187 DO ikind = 1, SIZE(wrk_tmp)
1188 IF (count(wrk_tmp == wrk_tmp(ikind)) > 1) THEN
1189 WRITE (iw, '(/A)') "In the IMOMM scheme no more than one QM atom can be bounded to the same MM atom."
1190 WRITE (iw, '(A)') "Multiple link MM atom not allowed. Check your link sections."
1191 cpabort("")
1192 END IF
1193 END DO
1194 DEALLOCATE (wrk_tmp)
1195 END IF
1196 ! PSEUDO
1197 IF (n_pseudo /= 0) THEN
1198 ALLOCATE (qmmm_links%pseudo(n_pseudo))
1199 DO ikind = 1, n_pseudo
1200 NULLIFY (qmmm_links%pseudo(ikind)%link)
1201 ALLOCATE (qmmm_links%pseudo(ikind)%link)
1202 END DO
1203 n_pseudo = 0
1204 DO ikind = 1, nlinks
1205 CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
1206 IF (link_type == do_qmmm_link_pseudo) THEN
1207 n_pseudo = n_pseudo + 1
1208 CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
1209 CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
1210 qmmm_links%pseudo(n_pseudo)%link%qm_index = qm_index
1211 qmmm_links%pseudo(n_pseudo)%link%mm_index = mm_index
1212 END IF
1213 END DO
1214 END IF
1215 ! GHO
1216 IF (n_gho /= 0) THEN
1217 ! not yet implemented
1218 ! still to define : type, implementation into QS
1219 cpabort("")
1220 END IF
1221 END SUBROUTINE setup_qmmm_links
1222
1223! **************************************************************************************************
1224!> \brief this routine sets up all variables to treat qmmm links
1225!> \param qmmm_section ...
1226!> \param move_mm_charges ...
1227!> \param add_mm_charges ...
1228!> \param mm_atom_chrg ...
1229!> \param mm_el_pot_radius ...
1230!> \param mm_el_pot_radius_corr ...
1231!> \param added_charges ...
1232!> \param mm_atom_index ...
1233!> \par History
1234!> 12.2004 created [tlaino]
1235!> \author Teodoro Laino
1236! **************************************************************************************************
1237 SUBROUTINE move_or_add_atoms(qmmm_section, move_mm_charges, add_mm_charges, &
1238 mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, &
1239 added_charges, mm_atom_index)
1240 TYPE(section_vals_type), POINTER :: qmmm_section
1241 LOGICAL, INTENT(OUT) :: move_mm_charges, add_mm_charges
1242 REAL(kind=dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, &
1243 mm_el_pot_radius_corr
1244 TYPE(add_set_type), POINTER :: added_charges
1245 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1246
1247 INTEGER :: i_add, icount, ikind, ind1, index1, &
1248 index2, n_add_tot, n_adds, n_move_tot, &
1249 n_moves, n_rep_val, nlinks
1250 LOGICAL :: explicit
1251 REAL(kind=dp) :: alpha, c_radius, charge, radius
1252 TYPE(section_vals_type), POINTER :: add_section, move_section, &
1253 qmmm_link_section
1254
1255 explicit = .false.
1256 move_mm_charges = .false.
1257 add_mm_charges = .false.
1258 NULLIFY (qmmm_link_section, move_section, add_section)
1259 qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
1260 CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
1261 cpassert(nlinks /= 0)
1262 icount = 0
1263 n_move_tot = 0
1264 n_add_tot = 0
1265 DO ikind = 1, nlinks
1266 move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
1267 i_rep_section=ikind)
1268 CALL section_vals_get(move_section, n_repetition=n_moves)
1269 add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", &
1270 i_rep_section=ikind)
1271 CALL section_vals_get(add_section, n_repetition=n_adds)
1272 n_move_tot = n_move_tot + n_moves
1273 n_add_tot = n_add_tot + n_adds
1274 END DO
1275 icount = n_move_tot + n_add_tot
1276 IF (n_add_tot /= 0) add_mm_charges = .true.
1277 IF (n_move_tot /= 0) move_mm_charges = .true.
1278 !
1279 ! create add_set_type
1280 !
1281 CALL create_add_set_type(added_charges, ndim=icount)
1282 !
1283 ! Fill in structures
1284 !
1285 icount = 0
1286 DO ikind = 1, nlinks
1287 move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
1288 i_rep_section=ikind)
1289 CALL section_vals_get(move_section, explicit=explicit, n_repetition=n_moves)
1290 !
1291 ! Moving charge atoms
1292 !
1293 IF (explicit) THEN
1294 DO i_add = 1, n_moves
1295 icount = icount + 1
1296 CALL section_vals_val_get(move_section, "ATOM_INDEX_1", i_val=index1, i_rep_section=i_add)
1297 CALL section_vals_val_get(move_section, "ATOM_INDEX_2", i_val=index2, i_rep_section=i_add)
1298 CALL section_vals_val_get(move_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
1299 CALL section_vals_val_get(move_section, "RADIUS", r_val=radius, i_rep_section=i_add)
1300 CALL section_vals_val_get(move_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
1301 c_radius = radius
1302 IF (n_rep_val == 1) &
1303 CALL section_vals_val_get(move_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add)
1304
1305 CALL set_add_set_type(added_charges, icount, index1, index2, alpha, radius, c_radius, &
1306 mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, &
1307 mm_el_pot_radius_corr=mm_el_pot_radius_corr, &
1308 mm_atom_index=mm_atom_index, move=n_moves, ind1=ind1)
1309 END DO
1310 mm_atom_chrg(ind1) = 0.0_dp
1311 END IF
1312
1313 add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", &
1314 i_rep_section=ikind)
1315 CALL section_vals_get(add_section, explicit=explicit, n_repetition=n_adds)
1316 !
1317 ! Adding charge atoms
1318 !
1319 IF (explicit) THEN
1320 DO i_add = 1, n_adds
1321 icount = icount + 1
1322 CALL section_vals_val_get(add_section, "ATOM_INDEX_1", i_val=index1, i_rep_section=i_add)
1323 CALL section_vals_val_get(add_section, "ATOM_INDEX_2", i_val=index2, i_rep_section=i_add)
1324 CALL section_vals_val_get(add_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
1325 CALL section_vals_val_get(add_section, "RADIUS", r_val=radius, i_rep_section=i_add)
1326 CALL section_vals_val_get(add_section, "CHARGE", r_val=charge, i_rep_section=i_add)
1327 CALL section_vals_val_get(add_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
1328 c_radius = radius
1329 IF (n_rep_val == 1) &
1330 CALL section_vals_val_get(add_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add)
1331
1332 CALL set_add_set_type(added_charges, icount, index1, index2, alpha, radius, c_radius, charge, &
1333 mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, &
1334 mm_el_pot_radius_corr=mm_el_pot_radius_corr, &
1335 mm_atom_index=mm_atom_index)
1336 END DO
1337 END IF
1338 END DO
1339
1340 END SUBROUTINE move_or_add_atoms
1341
1342! **************************************************************************************************
1343!> \brief this routine sets up all variables of the add_set_type type
1344!> \param added_charges ...
1345!> \param icount ...
1346!> \param Index1 ...
1347!> \param Index2 ...
1348!> \param alpha ...
1349!> \param radius ...
1350!> \param c_radius ...
1351!> \param charge ...
1352!> \param mm_atom_chrg ...
1353!> \param mm_el_pot_radius ...
1354!> \param mm_el_pot_radius_corr ...
1355!> \param mm_atom_index ...
1356!> \param move ...
1357!> \param ind1 ...
1358!> \par History
1359!> 12.2004 created [tlaino]
1360!> \author Teodoro Laino
1361! **************************************************************************************************
1362 SUBROUTINE set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, charge, &
1363 mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, mm_atom_index, move, ind1)
1364 TYPE(add_set_type), POINTER :: added_charges
1365 INTEGER, INTENT(IN) :: icount, index1, index2
1366 REAL(kind=dp), INTENT(IN) :: alpha, radius, c_radius
1367 REAL(kind=dp), INTENT(IN), OPTIONAL :: charge
1368 REAL(kind=dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, &
1369 mm_el_pot_radius_corr
1370 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1371 INTEGER, INTENT(in), OPTIONAL :: move
1372 INTEGER, INTENT(OUT), OPTIONAL :: ind1
1373
1374 INTEGER :: i, my_move
1375 REAL(kind=dp) :: my_c_radius, my_charge, my_radius
1376
1377 my_move = 0
1378 my_radius = radius
1379 my_c_radius = c_radius
1380 IF (PRESENT(charge)) my_charge = charge
1381 IF (PRESENT(move)) my_move = move
1382 i = 1
1383 getid: DO WHILE (i <= SIZE(mm_atom_index))
1384 IF (index1 == mm_atom_index(i)) EXIT getid
1385 i = i + 1
1386 END DO getid
1387 IF (PRESENT(ind1)) ind1 = i
1388 cpassert(i <= SIZE(mm_atom_index))
1389 IF (.NOT. PRESENT(charge)) my_charge = mm_atom_chrg(i)/real(my_move, kind=dp)
1390 IF (my_radius == 0.0_dp) my_radius = mm_el_pot_radius(i)
1391 IF (my_c_radius == 0.0_dp) my_c_radius = mm_el_pot_radius_corr(i)
1392
1393 added_charges%add_env(icount)%Index1 = index1
1394 added_charges%add_env(icount)%Index2 = index2
1395 added_charges%add_env(icount)%alpha = alpha
1396 added_charges%mm_atom_index(icount) = icount
1397 added_charges%mm_atom_chrg(icount) = my_charge
1398 added_charges%mm_el_pot_radius(icount) = my_radius
1399 added_charges%mm_el_pot_radius_corr(icount) = my_c_radius
1400 END SUBROUTINE set_add_set_type
1401
1402! **************************************************************************************************
1403!> \brief this routine sets up the origin of the MM cell respect to the
1404!> origin of the QM cell. The origin of the QM cell is assumed to be
1405!> in (0.0,0.0,0.0)...
1406!> \param qmmm_section ...
1407!> \param qmmm_env ...
1408!> \param qm_cell_small ...
1409!> \param dr ...
1410!> \par History
1411!> 02.2005 created [tlaino]
1412!> \author Teodoro Laino
1413! **************************************************************************************************
1414 SUBROUTINE setup_origin_mm_cell(qmmm_section, qmmm_env, qm_cell_small, &
1415 dr)
1416 TYPE(section_vals_type), POINTER :: qmmm_section
1417 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1418 TYPE(cell_type), POINTER :: qm_cell_small
1419 REAL(kind=dp), DIMENSION(3), INTENT(in) :: dr
1420
1421 LOGICAL :: center_grid
1422 REAL(kind=dp), DIMENSION(3) :: tmp
1423 REAL(kind=dp), DIMENSION(:), POINTER :: vec
1424
1425! This is the vector that corrects position to apply properly the PBC
1426
1427 tmp(1) = qm_cell_small%hmat(1, 1)
1428 tmp(2) = qm_cell_small%hmat(2, 2)
1429 tmp(3) = qm_cell_small%hmat(3, 3)
1430 cpassert(all(tmp > 0))
1431 qmmm_env%dOmmOqm = tmp/2.0_dp
1432 ! This is unit vector to translate the QM system in order to center it
1433 ! in QM cell
1434 CALL section_vals_val_get(qmmm_section, "CENTER_GRID", l_val=center_grid)
1435 IF (center_grid) THEN
1436 qmmm_env%utrasl = dr
1437 ELSE
1438 qmmm_env%utrasl = 1.0_dp
1439 END IF
1440 CALL section_vals_val_get(qmmm_section, "INITIAL_TRANSLATION_VECTOR", r_vals=vec)
1441 qmmm_env%transl_v = vec
1442 END SUBROUTINE setup_origin_mm_cell
1443
1444! **************************************************************************************************
1445!> \brief this routine sets up list of MM atoms carrying an image charge
1446!> \param image_charge_section ...
1447!> \param qmmm_env ...
1448!> \param qm_atom_index ...
1449!> \param subsys_mm ...
1450!> \par History
1451!> 02.2012 created
1452!> \author Dorothea Golze
1453! **************************************************************************************************
1454 SUBROUTINE setup_image_atom_list(image_charge_section, qmmm_env, &
1455 qm_atom_index, subsys_mm)
1456
1457 TYPE(section_vals_type), POINTER :: image_charge_section
1458 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1459 INTEGER, DIMENSION(:), POINTER :: qm_atom_index
1460 TYPE(cp_subsys_type), POINTER :: subsys_mm
1461
1462 INTEGER :: atom_a, atom_b, i, j, k, max_index, &
1463 n_var, num_const_atom, &
1464 num_image_mm_atom
1465 INTEGER, DIMENSION(:), POINTER :: mm_indexes
1466 LOGICAL :: fix_xyz, imageind_in_range
1467 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind
1468
1469 NULLIFY (mm_indexes, molecule_kind)
1470 imageind_in_range = .false.
1471 num_image_mm_atom = 0
1472 max_index = 0
1473
1474 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1475 n_rep_val=n_var)
1476 DO i = 1, n_var
1477 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1478 i_rep_val=i, i_vals=mm_indexes)
1479 num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
1480 END DO
1481
1482 ALLOCATE (qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom))
1483
1484 qmmm_env%image_charge_pot%image_mm_list = 0
1485 num_image_mm_atom = 1
1486
1487 DO i = 1, n_var
1488 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
1489 i_rep_val=i, i_vals=mm_indexes)
1490 qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom:num_image_mm_atom &
1491 + SIZE(mm_indexes) - 1) = mm_indexes(:)
1492 num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
1493 END DO
1494
1495 ! checking, if in range, if list contains QM atoms or any atoms doubled
1496 num_image_mm_atom = num_image_mm_atom - 1
1497
1498 max_index = SIZE(subsys_mm%particles%els)
1499
1500 cpassert(SIZE(qmmm_env%image_charge_pot%image_mm_list) /= 0)
1501 imageind_in_range = (maxval(qmmm_env%image_charge_pot%image_mm_list) <= max_index) &
1502 .AND. (minval(qmmm_env%image_charge_pot%image_mm_list) > 0)
1503 cpassert(imageind_in_range)
1504
1505 DO i = 1, num_image_mm_atom
1506 atom_a = qmmm_env%image_charge_pot%image_mm_list(i)
1507 IF (any(qm_atom_index == atom_a)) THEN
1508 cpabort("Image atom list must only contain MM atoms")
1509 END IF
1510 DO j = i + 1, num_image_mm_atom
1511 atom_b = qmmm_env%image_charge_pot%image_mm_list(j)
1512 IF (atom_a == atom_b) &
1513 cpabort("There are atoms doubled in image list.")
1514 END DO
1515 END DO
1516
1517 ! check if molecules in list carry constraints
1518 num_const_atom = 0
1519 fix_xyz = .true.
1520 IF (ASSOCIATED(subsys_mm%molecule_kinds)) THEN
1521 IF (ASSOCIATED(subsys_mm%molecule_kinds%els)) THEN
1522 molecule_kind => subsys_mm%molecule_kinds%els
1523 DO i = 1, SIZE(molecule_kind)
1524 IF (.NOT. ASSOCIATED(molecule_kind(i)%fixd_list)) EXIT
1525 IF (.NOT. fix_xyz) EXIT
1526 DO j = 1, SIZE(molecule_kind(i)%fixd_list)
1527 IF (.NOT. fix_xyz) EXIT
1528 DO k = 1, num_image_mm_atom
1529 atom_a = qmmm_env%image_charge_pot%image_mm_list(k)
1530 IF (atom_a == molecule_kind(i)%fixd_list(j)%fixd) THEN
1531 num_const_atom = num_const_atom + 1
1532 IF (molecule_kind(i)%fixd_list(j)%itype /= use_perd_xyz) THEN
1533 fix_xyz = .false.
1534 EXIT
1535 END IF
1536 END IF
1537 END DO
1538 END DO
1539 END DO
1540 END IF
1541 END IF
1542
1543 ! if all image atoms are constrained, calculate image matrix only
1544 ! once for the first MD or GEO_OPT step (for non-iterative case)
1545 IF (num_const_atom == num_image_mm_atom .AND. fix_xyz) THEN
1546 qmmm_env%image_charge_pot%state_image_matrix = calc_once
1547 ELSE
1548 qmmm_env%image_charge_pot%state_image_matrix = calc_always
1549 END IF
1550
1551 END SUBROUTINE setup_image_atom_list
1552
1553! **************************************************************************************************
1554!> \brief Print info on image charges
1555!> \param qmmm_env ...
1556!> \param qmmm_section ...
1557!> \par History
1558!> 03.2012 created
1559!> \author Dorothea Golze
1560! **************************************************************************************************
1561 SUBROUTINE print_image_charge_info(qmmm_env, qmmm_section)
1562
1563 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1564 TYPE(section_vals_type), POINTER :: qmmm_section
1565
1566 INTEGER :: iw
1567 REAL(kind=dp) :: eta, eta_conv, v0, v0_conv
1568 TYPE(cp_logger_type), POINTER :: logger
1569
1570 logger => cp_get_default_logger()
1571 iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%PROGRAM_RUN_INFO", &
1572 extension=".log")
1573 eta = qmmm_env%image_charge_pot%eta
1574 eta_conv = cp_unit_from_cp2k(eta, "angstrom", power=-2)
1575 v0 = qmmm_env%image_charge_pot%V0
1576 v0_conv = cp_unit_from_cp2k(v0, "volt")
1577
1578 IF (iw > 0) THEN
1579 WRITE (iw, fmt="(T25,A)") "IMAGE CHARGE PARAMETERS"
1580 WRITE (iw, fmt="(T25,A)") repeat("-", 23)
1581 WRITE (iw, fmt="(/)")
1582 WRITE (iw, fmt="(T2,A)") "INDEX OF MM ATOMS CARRYING AN IMAGE CHARGE:"
1583 WRITE (iw, fmt="(/)")
1584
1585 WRITE (iw, "(7X,10I6)") qmmm_env%image_charge_pot%image_mm_list
1586 WRITE (iw, fmt="(/)")
1587 WRITE (iw, "(T2,A52,T69,F12.8)") &
1588 "WIDTH OF GAUSSIAN CHARGE DISTRIBUTION [angstrom^-2]:", eta_conv
1589 WRITE (iw, "(T2,A26,T69,F12.8)") "EXTERNAL POTENTIAL [volt]:", v0_conv
1590 WRITE (iw, fmt="(/,T2,A,/)") repeat("-", 79)
1591 END IF
1592 CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
1593 "PRINT%PROGRAM_RUN_INFO")
1594
1595 END SUBROUTINE print_image_charge_info
1596
1597END 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:633
subroutine, public setup_qmmm_links(qmmm_section, qmmm_links, mm_el_pot_radius, mm_el_pot_radius_corr, mm_atom_index, iw)
this routine sets up all variables to treat qmmm links
Definition qmmm_init.F:1120
subroutine, public qmmm_init_gaussian_type(qmmm_env_qm, para_env, mm_atom_chrg, qs_env, added_charges, added_shells, print_section, qmmm_section)
...
Definition qmmm_init.F:368
subroutine, public move_or_add_atoms(qmmm_section, move_mm_charges, add_mm_charges, mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, added_charges, mm_atom_index)
this routine sets up all variables to treat qmmm links
Definition qmmm_init.F:1240
subroutine, public assign_mm_charges_and_radius(subsys, charges, mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, mm_atom_index, mm_link_atoms, mm_link_scale_factor, added_shells, shell_model)
Assigns charges and radius to evaluate the MM electrostatic potential.
Definition qmmm_init.F:115
subroutine, public setup_origin_mm_cell(qmmm_section, qmmm_env, qm_cell_small, dr)
this routine sets up the origin of the MM cell respect to the origin of the QM cell....
Definition qmmm_init.F:1416
subroutine, public qmmm_init_periodic_potential(qmmm_env_qm, qm_cell_small, mm_cell, para_env, qs_env, added_charges, added_shells, qmmm_periodic, print_section, mm_atom_chrg)
...
Definition qmmm_init.F:523
subroutine, public print_qmmm_charges(mm_atom_index, mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, added_charges, added_shells, qmmm_section, nocompatibility, shell_model)
Print info on charges generating the qmmm potential..
Definition qmmm_init.F:212
subroutine, public print_qmmm_links(qmmm_section, qmmm_links)
Print info on qm/mm links.
Definition qmmm_init.F:307
subroutine, public setup_qmmm_vars_mm(qmmm_section, qmmm_env, qm_atom_index, mm_link_atoms, mm_link_scale_factor, fist_scale_charge_link, qmmm_coupl_type, qmmm_link)
...
Definition qmmm_init.F:857
subroutine, public print_image_charge_info(qmmm_env, qmmm_section)
Print info on image charges.
Definition qmmm_init.F:1562
subroutine, public qmmm_init_potential(qmmm_env_qm, mm_cell, added_charges, added_shells, print_section)
...
Definition qmmm_init.F:467
Setting up the potential for QM/MM periodic boundary conditions calculations.
subroutine, public qmmm_per_potential_init(qmmm_coupl_type, per_potentials, potentials, pgfs, qm_cell_small, mm_cell, para_env, compatibility, qmmm_periodic, print_section, eps_mm_rspace, maxchrg, ncp, ncpl)
Initialize the QMMM potential stored on vector, according the qmmm_coupl_type.
subroutine, public qmmm_ewald_potential_init(ewald_env, ewald_pw, qmmm_coupl_type, mm_cell, para_env, qmmm_periodic, print_section)
Initialize the QMMM Ewald potential needed for QM-QM Coupling using point charges.
subroutine, public create_add_shell_type(added_shells, ndim)
creates the add_shell_type structure
subroutine, public create_add_set_type(added_charges, ndim)
creates the add_set_type structure
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
elemental subroutine, public get_shell(shell, charge, charge_core, charge_shell, mass_core, mass_shell, k2_spring, k4_spring, max_dist, shell_cutoff)
...
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