LCOV - code coverage report
Current view: top level - src - qmmm_init.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 594 624 95.2 %
Date: 2024-03-29 07:50:05 Functions: 16 16 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Initialize a QM/MM calculation
      10             : !> \par History
      11             : !>      5.2004 created [fawzi]
      12             : !> \author Fawzi Mohamed
      13             : ! **************************************************************************************************
      14             : MODULE qmmm_init
      15             :    USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
      16             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      17             :                                               get_atomic_kind,&
      18             :                                               set_atomic_kind
      19             :    USE cell_methods,                    ONLY: read_cell
      20             :    USE cell_types,                      ONLY: cell_type,&
      21             :                                               use_perd_xyz
      22             :    USE cp_control_types,                ONLY: dft_control_type
      23             :    USE cp_eri_mme_interface,            ONLY: cp_eri_mme_init_read_input,&
      24             :                                               cp_eri_mme_set_params
      25             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      26             :                                               cp_logger_type
      27             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      28             :                                               cp_print_key_unit_nr
      29             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      30             :                                               cp_subsys_type
      31             :    USE cp_units,                        ONLY: cp_unit_from_cp2k,&
      32             :                                               cp_unit_to_cp2k
      33             :    USE external_potential_types,        ONLY: fist_potential_type,&
      34             :                                               get_potential,&
      35             :                                               set_potential
      36             :    USE force_field_types,               ONLY: input_info_type
      37             :    USE force_fields_input,              ONLY: read_gd_section,&
      38             :                                               read_gp_section,&
      39             :                                               read_lj_section,&
      40             :                                               read_wl_section
      41             :    USE input_constants,                 ONLY: &
      42             :         RADIUS_QMMM_DEFAULT, calc_always, calc_once, do_eri_mme, do_qmmm_gauss, &
      43             :         do_qmmm_image_calcmatrix, do_qmmm_image_iter, do_qmmm_link_gho, do_qmmm_link_imomm, &
      44             :         do_qmmm_link_pseudo, do_qmmm_pcharge, do_qmmm_swave
      45             :    USE input_section_types,             ONLY: section_vals_get,&
      46             :                                               section_vals_get_subs_vals,&
      47             :                                               section_vals_type,&
      48             :                                               section_vals_val_get
      49             :    USE kinds,                           ONLY: default_string_length,&
      50             :                                               dp
      51             :    USE message_passing,                 ONLY: mp_para_env_type
      52             :    USE molecule_kind_types,             ONLY: molecule_kind_type
      53             :    USE pair_potential_types,            ONLY: pair_potential_reallocate
      54             :    USE particle_list_types,             ONLY: particle_list_type
      55             :    USE particle_types,                  ONLY: particle_type
      56             :    USE pw_env_types,                    ONLY: pw_env_type
      57             :    USE qmmm_elpot,                      ONLY: qmmm_potential_init
      58             :    USE qmmm_ff_fist,                    ONLY: qmmm_ff_precond_only_qm
      59             :    USE qmmm_gaussian_init,              ONLY: qmmm_gaussian_initialize
      60             :    USE qmmm_per_elpot,                  ONLY: qmmm_ewald_potential_init,&
      61             :                                               qmmm_per_potential_init
      62             :    USE qmmm_types_low,                  ONLY: add_set_type,&
      63             :                                               add_shell_type,&
      64             :                                               create_add_set_type,&
      65             :                                               create_add_shell_type,&
      66             :                                               qmmm_env_mm_type,&
      67             :                                               qmmm_env_qm_type,&
      68             :                                               qmmm_links_type
      69             :    USE qs_environment_types,            ONLY: get_qs_env,&
      70             :                                               qs_environment_type
      71             :    USE shell_potential_types,           ONLY: get_shell,&
      72             :                                               shell_kind_type
      73             : #include "./base/base_uses.f90"
      74             : 
      75             :    IMPLICIT NONE
      76             :    PRIVATE
      77             : 
      78             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      79             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_init'
      80             : 
      81             :    PUBLIC :: assign_mm_charges_and_radius, &
      82             :              print_qmmm_charges, &
      83             :              print_qmmm_links, &
      84             :              print_image_charge_info, &
      85             :              qmmm_init_gaussian_type, &
      86             :              qmmm_init_potential, &
      87             :              qmmm_init_periodic_potential, &
      88             :              setup_qmmm_vars_qm, &
      89             :              setup_qmmm_vars_mm, &
      90             :              setup_qmmm_links, &
      91             :              move_or_add_atoms, &
      92             :              setup_origin_mm_cell
      93             : 
      94             : CONTAINS
      95             : 
      96             : ! **************************************************************************************************
      97             : !> \brief Assigns charges and radius to evaluate the MM electrostatic potential
      98             : !> \param subsys the subsys containing the MM charges
      99             : !> \param charges ...
     100             : !> \param mm_atom_chrg ...
     101             : !> \param mm_el_pot_radius ...
     102             : !> \param mm_el_pot_radius_corr ...
     103             : !> \param mm_atom_index ...
     104             : !> \param mm_link_atoms ...
     105             : !> \param mm_link_scale_factor ...
     106             : !> \param added_shells ...
     107             : !> \param shell_model ...
     108             : !> \par History
     109             : !>      06.2004 created [tlaino]
     110             : !> \author Teodoro Laino
     111             : ! **************************************************************************************************
     112         394 :    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         394 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_set, particle_set, shell_set
     132             :       TYPE(shell_kind_type), POINTER                     :: shell_kind
     133             : 
     134         394 :       NULLIFY (particle_set, my_kind, added_shells)
     135             :       CALL cp_subsys_get(subsys=subsys, particles=particles, core_particles=core_particles, &
     136         394 :                          shell_particles=shell_particles)
     137         394 :       particle_set => particles%els
     138             : 
     139      190762 :       IF (ALL(particle_set(:)%shell_index .EQ. 0)) THEN
     140         392 :          shell_model = .FALSE.
     141         392 :          CALL create_add_shell_type(added_shells, ndim=0)
     142             :       ELSE
     143           2 :          shell_model = .TRUE.
     144             :       END IF
     145             : 
     146         394 :       IF (shell_model) THEN
     147           2 :          shell_set => shell_particles%els
     148           2 :          core_set => core_particles%els
     149           2 :          ishell = SIZE(shell_set)
     150           2 :          CALL create_add_shell_type(added_shells, ndim=ishell)
     151           2 :          added_shells%added_particles => shell_set
     152           2 :          added_shells%added_cores => core_set
     153             :       END IF
     154             : 
     155      188174 :       DO I = 1, SIZE(mm_atom_index)
     156      187780 :          IndMM = mm_atom_index(I)
     157      187780 :          my_kind => particle_set(IndMM)%atomic_kind
     158             :          CALL get_atomic_kind(atomic_kind=my_kind, fist_potential=my_potential, &
     159      187780 :                               shell_active=is_shell, shell=shell_kind)
     160             :          CALL get_potential(potential=my_potential, &
     161             :                             qeff=qi, &
     162             :                             qmmm_radius=ri, &
     163      187780 :                             qmmm_corr_radius=rc)
     164      187780 :          IF (ASSOCIATED(charges)) qi = charges(IndMM)
     165      187780 :          mm_atom_chrg(I) = qi
     166      187780 :          mm_el_pot_radius(I) = ri
     167      187780 :          mm_el_pot_radius_corr(I) = rc
     168      375954 :          IF (is_shell) THEN
     169          56 :             IndShell = particle_set(IndMM)%shell_index
     170          56 :             IF (ASSOCIATED(shell_kind)) THEN
     171             :                CALL get_shell(shell=shell_kind, &
     172             :                               charge_core=qcore, &
     173          56 :                               charge_shell=qshell)
     174          56 :                mm_atom_chrg(I) = qcore
     175             :             END IF
     176          56 :             added_shells%mm_core_index(IndShell) = IndMM
     177          56 :             added_shells%mm_core_chrg(IndShell) = qshell
     178          56 :             added_shells%mm_el_pot_radius(Indshell) = ri*1.0_dp
     179          56 :             added_shells%mm_el_pot_radius_corr(Indshell) = rc*1.0_dp
     180             :          END IF
     181             :       END DO
     182             : 
     183         394 :       IF (ASSOCIATED(mm_link_atoms)) THEN
     184         256 :          DO ilink = 1, SIZE(mm_link_atoms)
     185       40710 :             DO i = 1, SIZE(mm_atom_index)
     186       40710 :                IF (mm_atom_index(i) == mm_link_atoms(ilink)) EXIT
     187             :             END DO
     188         194 :             IndMM = mm_atom_index(I)
     189         256 :             mm_atom_chrg(i) = mm_atom_chrg(i)*mm_link_scale_factor(ilink)
     190             :          END DO
     191             :       END IF
     192             : 
     193         394 :    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         394 :    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         394 :       qtot = 0.0_dp
     225         394 :       logger => cp_get_default_logger()
     226             :       iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%QMMM_CHARGES", &
     227         394 :                                 extension=".log")
     228         394 :       IF (iw > 0) THEN
     229         169 :          WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
     230         169 :          WRITE (iw, FMT='(/5X,A)') "MM    POINT CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
     231         169 :          WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
     232       23623 :          DO I = 1, SIZE(mm_atom_index)
     233       23454 :             IndMM = mm_atom_index(I)
     234       23454 :             qi = mm_atom_chrg(I)
     235       23454 :             qtot = qtot + qi
     236       23454 :             ri = mm_el_pot_radius(I)
     237       23454 :             rc = mm_el_pot_radius_corr(I)
     238       23623 :             IF (nocompatibility) THEN
     239        2177 :                WRITE (iw, '(5X,A9,T15,I5,T28,A8,T38,F12.6,T60,A8,T69,F12.6)') ' MM ATOM:', IndMM, ' RADIUS:', ri, &
     240        4354 :                   ' 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       21277 :                   ' MM ATOM:', IndMM, ' RADIUS:', ri, ' CHARGE:', qi, 'CORR. RADIUS', rc
     244             :             END IF
     245             :          END DO
     246         169 :          IF (added_charges%num_mm_atoms /= 0) THEN
     247           4 :             WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
     248           4 :             WRITE (iw, '(/5X,A)') "ADDED POINT CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
     249           4 :             WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
     250          18 :             DO I = 1, SIZE(added_charges%mm_atom_index)
     251          14 :                IndMM = added_charges%mm_atom_index(I)
     252          14 :                qi = added_charges%mm_atom_chrg(I)
     253          14 :                qtot = qtot + qi
     254          14 :                ri = added_charges%mm_el_pot_radius(I)
     255          14 :                ind1 = added_charges%add_env(I)%Index1
     256          14 :                ind2 = added_charges%add_env(I)%Index2
     257          18 :                IF (nocompatibility) THEN
     258          14 :                   WRITE (iw, '(5X,A9,I5,T25,A8,T35,F12.6,T50,A8,T59,F12.6,I5,I5)') 'MM POINT:', IndMM, ' RADIUS:', ri, &
     259          28 :                      ' 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           0 :                      'MM POINT:', IndMM, ' RADIUS:', ri, ' CHARGE:', qi, ind1, ind2, 'CORR. RADIUS', rc
     263             :                END IF
     264             :             END DO
     265             : 
     266             :          END IF
     267             : 
     268         169 :          IF (shell_model) THEN
     269           1 :             WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 73)
     270           1 :             WRITE (iw, '(/5X,A)') "ADDED SHELL CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
     271           1 :             WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 73)
     272             : 
     273          29 :             DO I = 1, SIZE(added_shells%mm_core_index)
     274          28 :                IndMM = added_shells%mm_core_index(I)
     275          28 :                qi = added_shells%mm_core_chrg(I)
     276          28 :                qtot = qtot + qi
     277          28 :                ri = added_shells%mm_el_pot_radius(I)
     278          29 :                IF (nocompatibility) THEN
     279          28 :                   WRITE (iw, '(7X,A,I5,A8,F12.6,A8,F12.6,3F12.6)') 'SHELL:', IndMM, ' RADIUS:', ri, &
     280         140 :                      ' CHARGE:', qi, added_shells%added_particles(I)%r
     281             :                ELSE
     282           0 :                   WRITE (iw, '(7X,A,I5,A8,F12.6,A8,F12.6,A,F12.6)') 'SHELL:', IndMM, ' RADIUS:', ri, &
     283           0 :                      ' CHARGE:', qi, ' CORR. RADIUS', rc
     284             :                END IF
     285             : 
     286             :             END DO
     287             : 
     288             :          END IF
     289             : 
     290         169 :          WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
     291         169 :          WRITE (iw, '(/,T50,A,T69,F12.6)') ' TOTAL CHARGE:', qtot
     292         169 :          WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 79)
     293             :       END IF
     294             :       CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
     295         394 :                                         "PRINT%QMMM_CHARGES")
     296         394 :    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          64 :    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          64 :       logger => cp_get_default_logger()
     315          64 :       iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%qmmm_link_info", extension=".log")
     316          64 :       IF (iw > 0) THEN
     317          21 :          IF (ASSOCIATED(qmmm_links)) THEN
     318          21 :             WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
     319          21 :             WRITE (iw, FMT="(/,T31,A)") " QM/MM LINKS "
     320          21 :             WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
     321          21 :             IF (ASSOCIATED(qmmm_links%imomm)) THEN
     322          20 :                WRITE (iw, FMT="(/,T31,A)") " IMOMM TYPE LINK "
     323          56 :                DO I = 1, SIZE(qmmm_links%imomm)
     324          36 :                   qm_index = qmmm_links%imomm(I)%link%qm_index
     325          36 :                   mm_index = qmmm_links%imomm(I)%link%mm_index
     326          36 :                   alpha = qmmm_links%imomm(I)%link%alpha
     327          36 :                   WRITE (iw, FMT="(T2,A,T20,A9,I8,1X,A9,I8,T62,A6,F12.6)") "TYPE: IMOMM", &
     328          92 :                      "QM INDEX:", qm_index, "MM INDEX:", mm_index, "ALPHA:", alpha
     329             :                END DO
     330             :             END IF
     331          21 :             IF (ASSOCIATED(qmmm_links%pseudo)) THEN
     332           1 :                WRITE (iw, FMT="(/,T31,A)") " PSEUDO TYPE LINK "
     333           3 :                DO I = 1, SIZE(qmmm_links%pseudo)
     334           2 :                   qm_index = qmmm_links%pseudo(I)%link%qm_index
     335           2 :                   mm_index = qmmm_links%pseudo(I)%link%mm_index
     336           2 :                   WRITE (iw, FMT="(T2,A,T20,A9,I8,1X,A9,I8)") "TYPE: PSEUDO", &
     337           5 :                      "QM INDEX:", qm_index, "MM INDEX:", mm_index
     338             :                END DO
     339             :             END IF
     340          21 :             WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 73)
     341             :          ELSE
     342           0 :             WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
     343           0 :             WRITE (iw, FMT="(/,T26,A)") " NO QM/MM LINKS DETECTED"
     344           0 :             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          64 :                                         "PRINT%qmmm_link_info")
     349          64 :    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         394 :    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         394 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: maxradius, maxradius2
     379             :       TYPE(pw_env_type), POINTER                         :: pw_env
     380             : 
     381         394 :       NULLIFY (maxradius, maxradius2, pw_env)
     382             : 
     383      188568 :       maxchrg = MAXVAL(ABS(mm_atom_chrg(:)))
     384         394 :       CALL get_qs_env(qs_env, pw_env=pw_env)
     385         410 :       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         394 :                                     qmmm_section=qmmm_section)
     398             : 
     399         394 :       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           8 :                                        qmmm_section=qmmm_section)
     412             : 
     413          16 :          SELECT CASE (qmmm_env_qm%qmmm_coupl_type)
     414             :          CASE (do_qmmm_gauss, do_qmmm_swave, do_qmmm_pcharge)
     415          40 :             DO i = 1, SIZE(maxradius)
     416          40 :                maxradius(i) = MAX(maxradius(i), maxradius2(i))
     417             :             END DO
     418             :          END SELECT
     419             : 
     420           8 :          IF (ASSOCIATED(maxradius2)) DEALLOCATE (maxradius2)
     421             :       END IF
     422             : 
     423         394 :       IF (qmmm_env_qm%added_shells%num_mm_atoms .GT. 0) THEN
     424             : 
     425          60 :          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           2 :                                        qmmm_section=qmmm_section)
     438             : 
     439           4 :          SELECT CASE (qmmm_env_qm%qmmm_coupl_type)
     440             :          CASE (do_qmmm_gauss, do_qmmm_swave, do_qmmm_pcharge)
     441          10 :             DO i = 1, SIZE(maxradius)
     442          10 :                maxradius(i) = MAX(maxradius(i), maxradius2(i))
     443             :             END DO
     444             :          END SELECT
     445             : 
     446           2 :          IF (ASSOCIATED(maxradius2)) DEALLOCATE (maxradius2)
     447             : 
     448             :       END IF
     449             : 
     450         394 :       qmmm_env_qm%maxradius => maxradius
     451             : 
     452         394 :    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         394 :    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         394 :                                print_section=print_section)
     480             : 
     481         394 :       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           8 :                                   print_section=print_section)
     490             :       END IF
     491             : 
     492         394 :       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           2 :                                   print_section=print_section)
     501             :       END IF
     502             : 
     503         394 :    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         394 :    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         394 :       IF (qmmm_env_qm%periodic) THEN
     536             : 
     537          46 :          NULLIFY (dft_control)
     538          46 :          CALL get_qs_env(qs_env, dft_control=dft_control)
     539             : 
     540          46 :          IF (dft_control%qs_control%semi_empirical) THEN
     541           0 :             CPABORT("QM/MM periodic calculations not implemented for semi empirical methods")
     542          46 :          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           4 :                                            para_env=para_env, qmmm_periodic=qmmm_periodic, print_section=print_section)
     546          42 :          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           4 :                                            para_env=para_env, qmmm_periodic=qmmm_periodic, print_section=print_section)
     550             :          ELSE
     551             : 
     552             :             ! setup for GPW/GPAW
     553        1560 :             maxchrg = MAXVAL(ABS(mm_atom_chrg(:)))
     554          38 :             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          38 :                                          ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
     570             : 
     571          38 :             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           0 :                                             ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
     587             :             END IF
     588             : 
     589          38 :             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           2 :                                             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         394 :    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         394 :    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         394 :       INTEGER, DIMENSION(:), POINTER                     :: mm_link_atoms
     649             :       LOGICAL                                            :: explicit, is_mm, is_qm
     650             :       REAL(KIND=dp)                                      :: tmp_radius, tmp_radius_c
     651         394 :       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         394 :       NULLIFY (mm_link_atoms, cell_section, tmp_sph_cut)
     659         394 :       NULLIFY (image_charge_section)
     660         394 :       qmmm_link = .FALSE.
     661             : 
     662         394 :       CALL section_vals_get(qmmm_section, explicit=explicit)
     663         394 :       IF (explicit) THEN
     664             :          !
     665             :          ! QM_CELL
     666             :          !
     667         394 :          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         394 :                         check_for_ref=.FALSE., para_env=para_env)
     670         394 :          qm_cell_small%tag = "CELL_QM"
     671         394 :          CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type)
     672         394 :          CALL section_vals_val_get(qmmm_section, "EPS_MM_RSPACE", r_val=eps_mm_rspace)
     673         394 :          CALL section_vals_val_get(qmmm_section, "SPHERICAL_CUTOFF", r_vals=tmp_sph_cut)
     674         394 :          CPASSERT(SIZE(tmp_sph_cut) == 2)
     675        1970 :          qmmm_env%spherical_cutoff = tmp_sph_cut
     676         394 :          IF (qmmm_env%spherical_cutoff(1) <= 0.0_dp) THEN
     677         392 :             qmmm_env%spherical_cutoff(2) = 0.0_dp
     678             :          ELSE
     679           2 :             IF (qmmm_env%spherical_cutoff(2) <= 0.0_dp) qmmm_env%spherical_cutoff(2) = EPSILON(0.0_dp)
     680           2 :             tmp_radius = qmmm_env%spherical_cutoff(1) - 20.0_dp*qmmm_env%spherical_cutoff(2)
     681           2 :             IF (tmp_radius <= 0.0_dp) &
     682             :                CALL cp_abort(__LOCATION__, &
     683             :                              "SPHERICAL_CUTOFF(1) > 20*SPHERICAL_CUTOFF(1)! Please correct parameters for "// &
     684           0 :                              "the Spherical Cutoff in order to satisfy the previous condition!")
     685             :          END IF
     686             :          !
     687             :          ! Initialization of arrays and core_charge_radius...
     688             :          !
     689         394 :          tmp_radius = 0.0_dp
     690         394 :          CALL cp_subsys_get(subsys=subsys_mm, atomic_kinds=atomic_kinds)
     691        4064 :          DO Ikind = 1, SIZE(atomic_kinds%els)
     692        3670 :             atomic_kind => atomic_kinds%els(Ikind)
     693             :             CALL get_atomic_kind(atomic_kind=atomic_kind, &
     694        3670 :                                  fist_potential=fist_potential)
     695             :             CALL set_potential(potential=fist_potential, &
     696             :                                qmmm_radius=tmp_radius, &
     697        3670 :                                qmmm_corr_radius=tmp_radius)
     698             :             CALL set_atomic_kind(atomic_kind=atomic_kind, &
     699        4064 :                                  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         394 :                                  qmmm_link=qmmm_link)
     706             :          !
     707             :          ! MM_KINDS
     708             :          !
     709         394 :          mm_kinds => section_vals_get_subs_vals(qmmm_section, "MM_KIND")
     710         394 :          CALL section_vals_get(mm_kinds, explicit=explicit, n_repetition=nkind)
     711             :          !
     712             :          ! Default
     713             :          !
     714         394 :          tmp_radius = cp_unit_to_cp2k(RADIUS_QMMM_DEFAULT, "angstrom")
     715        4064 :          Set_Radius_Pot_0: DO IkindR = 1, SIZE(atomic_kinds%els)
     716        3670 :             atomic_kind => atomic_kinds%els(IkindR)
     717        3670 :             CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
     718             :             CALL get_atomic_kind(atomic_kind=atomic_kind, &
     719        3670 :                                  fist_potential=fist_potential)
     720             :             CALL set_potential(potential=fist_potential, qmmm_radius=tmp_radius, &
     721        3670 :                                qmmm_corr_radius=tmp_radius)
     722             :             CALL set_atomic_kind(atomic_kind=atomic_kind, &
     723        4064 :                                  fist_potential=fist_potential)
     724             :          END DO Set_Radius_Pot_0
     725             :          !
     726             :          ! If present overwrite the kind specified in input file...
     727             :          !
     728         394 :          IF (explicit) THEN
     729         738 :             DO ikind = 1, nkind
     730             :                CALL section_vals_val_get(mm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
     731         504 :                                          c_val=mm_atom_kind)
     732         504 :                CALL section_vals_val_get(mm_kinds, "RADIUS", i_rep_section=ikind, r_val=tmp_radius)
     733         504 :                tmp_radius_c = tmp_radius
     734         504 :                CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
     735         504 :                IF (n_rep_val == 1) CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, &
     736           2 :                                                              r_val=tmp_radius_c)
     737        7254 :                Set_Radius_Pot_1: DO IkindR = 1, SIZE(atomic_kinds%els)
     738        6012 :                   atomic_kind => atomic_kinds%els(IkindR)
     739        6012 :                   CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
     740        6012 :                   is_qm = qmmm_ff_precond_only_qm(atmname)
     741        6516 :                   IF (TRIM(mm_atom_kind) == atmname) THEN
     742             :                      CALL get_atomic_kind(atomic_kind=atomic_kind, &
     743         870 :                                           fist_potential=fist_potential)
     744             :                      CALL set_potential(potential=fist_potential, &
     745             :                                         qmmm_radius=tmp_radius, &
     746         870 :                                         qmmm_corr_radius=tmp_radius_c)
     747             :                      CALL set_atomic_kind(atomic_kind=atomic_kind, &
     748         870 :                                           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         394 :          image_charge_section => section_vals_get_subs_vals(qmmm_section, "IMAGE_CHARGE")
     757         394 :          CALL section_vals_get(image_charge_section, explicit=qmmm_env%image_charge)
     758             : 
     759             :       ELSE
     760           0 :          CPABORT("QMMM section not present in input file!")
     761             :       END IF
     762             :       !
     763             :       ! Build MM atoms list
     764             :       !
     765         394 :       size_mm_system = SIZE(subsys_mm%particles%els) - SIZE(qm_atom_index)
     766         394 :       IF (qmmm_link .AND. ASSOCIATED(mm_link_atoms)) size_mm_system = size_mm_system + SIZE(mm_link_atoms)
     767        1180 :       ALLOCATE (mm_atom_index(size_mm_system))
     768         394 :       icount = 0
     769             : 
     770      190872 :       DO i = 1, SIZE(subsys_mm%particles%els)
     771      190478 :          is_mm = .TRUE.
     772     7211572 :          IF (ANY(qm_atom_index == i)) THEN
     773        2892 :             is_mm = .FALSE.
     774             :          END IF
     775      190478 :          IF (ASSOCIATED(mm_link_atoms)) THEN
     776      790244 :             IF (ANY(mm_link_atoms == i) .AND. qmmm_link) is_mm = .TRUE.
     777             :          END IF
     778      190678 :          IF (is_mm) THEN
     779      187780 :             icount = icount + 1
     780      187780 :             IF (icount <= size_mm_system) mm_atom_index(icount) = i
     781             :          END IF
     782             :       END DO
     783         394 :       CPASSERT(icount == size_mm_system)
     784         394 :       IF (ASSOCIATED(mm_link_atoms)) THEN
     785          62 :          DEALLOCATE (mm_link_atoms)
     786             :       END IF
     787             : 
     788             :       ! Build image charge atom list + set up variables
     789             :       !
     790         394 :       IF (qmmm_env%image_charge) THEN
     791             :          CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
     792          10 :                                    explicit=explicit)
     793          10 :          IF (explicit) qmmm_env%image_charge_pot%all_mm = .FALSE.
     794             : 
     795          10 :          IF (qmmm_env%image_charge_pot%all_mm) THEN
     796           0 :             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          10 :                                        qm_atom_index, subsys_mm)
     800             :          END IF
     801             : 
     802          10 :          qmmm_env%image_charge_pot%particles_all => subsys_mm%particles%els
     803             : 
     804             :          CALL section_vals_val_get(image_charge_section, "EXT_POTENTIAL", &
     805          10 :                                    r_val=qmmm_env%image_charge_pot%V0)
     806             :          CALL section_vals_val_get(image_charge_section, "WIDTH", &
     807          10 :                                    r_val=qmmm_env%image_charge_pot%eta)
     808             :          CALL section_vals_val_get(image_charge_section, "DETERM_COEFF", &
     809          10 :                                    i_val=my_type)
     810           8 :          SELECT CASE (my_type)
     811             :          CASE (do_qmmm_image_calcmatrix)
     812           8 :             qmmm_env%image_charge_pot%coeff_iterative = .FALSE.
     813             :          CASE (do_qmmm_image_iter)
     814          10 :             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          10 :                                    l_val=qmmm_env%image_charge_pot%image_restart)
     819             : 
     820             :          CALL section_vals_val_get(image_charge_section, "IMAGE_MATRIX_METHOD", &
     821          10 :                                    i_val=qmmm_env%image_charge_pot%image_matrix_method)
     822             : 
     823          10 :          IF (qmmm_env%image_charge_pot%image_matrix_method .EQ. do_eri_mme) THEN
     824           8 :             eri_mme_section => section_vals_get_subs_vals(image_charge_section, "ERI_MME")
     825           8 :             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           8 :                                        para_env=para_env)
     833             : 
     834             :          END IF
     835             :       END IF
     836             : 
     837         394 :    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         394 :    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         394 :       NULLIFY (qmmm_ff_section)
     869         394 :       qmmm_link = .FALSE.
     870         394 :       CALL section_vals_get(qmmm_section, explicit=explicit)
     871         394 :       IF (explicit) THEN
     872         394 :          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         394 :                                  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         394 :          qmmm_ff_section => section_vals_get_subs_vals(qmmm_section, "FORCEFIELD")
     880         394 :          CALL section_vals_get(qmmm_ff_section, explicit=qmmm_env%use_qmmm_ff)
     881         394 :          IF (qmmm_env%use_qmmm_ff) THEN
     882             :             CALL section_vals_val_get(qmmm_ff_section, "MULTIPLE_POTENTIAL", &
     883          20 :                                       l_val=qmmm_env%multiple_potential)
     884          20 :             CALL read_qmmm_ff_section(qmmm_ff_section, qmmm_env%inp_info)
     885             :          END IF
     886             :       END IF
     887         394 :    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         180 :    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          20 :       lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%LENNARD-JONES")
     911          20 :       wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%WILLIAMS")
     912          20 :       gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GOODWIN")
     913          20 :       gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GENPOT")
     914          20 :       CALL section_vals_get(lj_section, n_repetition=n_lj)
     915          20 :       np = n_lj
     916          20 :       IF (n_lj /= 0) THEN
     917          18 :          CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, lj_charmm=.TRUE.)
     918          18 :          CALL read_lj_section(inp_info%nonbonded, lj_section, start=0)
     919             :       END IF
     920          20 :       CALL section_vals_get(wl_section, n_repetition=n_wl)
     921          20 :       np = n_lj + n_wl
     922          20 :       IF (n_wl /= 0) THEN
     923           2 :          CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, williams=.TRUE.)
     924           2 :          CALL read_wl_section(inp_info%nonbonded, wl_section, start=n_lj)
     925             :       END IF
     926          20 :       CALL section_vals_get(gd_section, n_repetition=n_gd)
     927          20 :       np = n_lj + n_wl + n_gd
     928          20 :       IF (n_gd /= 0) THEN
     929           0 :          CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, goodwin=.TRUE.)
     930           0 :          CALL read_gd_section(inp_info%nonbonded, gd_section, start=n_lj + n_wl)
     931             :       END IF
     932          20 :       CALL section_vals_get(gp_section, n_repetition=n_gp)
     933          20 :       np = n_lj + n_wl + n_gd + n_gp
     934          20 :       IF (n_gp /= 0) THEN
     935           0 :          CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, gp=.TRUE.)
     936           0 :          CALL read_gp_section(inp_info%nonbonded, gp_section, start=n_lj + n_wl + n_gd)
     937             :       END IF
     938             :       !
     939             :       ! NONBONDED14
     940             :       !
     941          20 :       lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%LENNARD-JONES")
     942          20 :       wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%WILLIAMS")
     943          20 :       gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GOODWIN")
     944          20 :       gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GENPOT")
     945          20 :       CALL section_vals_get(lj_section, n_repetition=n_lj)
     946          20 :       np = n_lj
     947          20 :       IF (n_lj /= 0) THEN
     948           0 :          CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, lj_charmm=.TRUE.)
     949           0 :          CALL read_lj_section(inp_info%nonbonded14, lj_section, start=0)
     950             :       END IF
     951          20 :       CALL section_vals_get(wl_section, n_repetition=n_wl)
     952          20 :       np = n_lj + n_wl
     953          20 :       IF (n_wl /= 0) THEN
     954           0 :          CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, williams=.TRUE.)
     955           0 :          CALL read_wl_section(inp_info%nonbonded14, wl_section, start=n_lj)
     956             :       END IF
     957          20 :       CALL section_vals_get(gd_section, n_repetition=n_gd)
     958          20 :       np = n_lj + n_wl + n_gd
     959          20 :       IF (n_gd /= 0) THEN
     960           0 :          CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, goodwin=.TRUE.)
     961           0 :          CALL read_gd_section(inp_info%nonbonded14, gd_section, start=n_lj + n_wl)
     962             :       END IF
     963          20 :       CALL section_vals_get(gp_section, n_repetition=n_gp)
     964          20 :       np = n_lj + n_wl + n_gd + n_gp
     965          20 :       IF (n_gp /= 0) THEN
     966           0 :          CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, gp=.TRUE.)
     967           0 :          CALL read_gp_section(inp_info%nonbonded14, gp_section, start=n_lj + n_wl + n_gd)
     968             :       END IF
     969             : 
     970          20 :    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        2364 :    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         788 :       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         788 :       num_qm_atom_tot = 0
    1006         788 :       link_involv_mm = 0
    1007         788 :       nlinks = 0
    1008             :       !
    1009             :       ! QM_KINDS
    1010             :       !
    1011        1576 :       qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND")
    1012         788 :       CALL section_vals_get(qm_kinds, n_repetition=nkind)
    1013        2596 :       DO ikind = 1, nkind
    1014        1808 :          CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
    1015        6336 :          DO k = 1, n_var
    1016             :             CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
    1017        3740 :                                       i_vals=mm_indexes)
    1018        5548 :             num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
    1019             :          END DO
    1020             :       END DO
    1021             :       !
    1022             :       ! QM/MM LINKS
    1023             :       !
    1024         788 :       qmmm_links => section_vals_get_subs_vals(qmmm_section, "LINK")
    1025         788 :       CALL section_vals_get(qmmm_links, explicit=explicit)
    1026         788 :       IF (explicit) THEN
    1027         128 :          qmmm_link = .TRUE.
    1028         128 :          CALL section_vals_get(qmmm_links, n_repetition=nlinks)
    1029             :          ! Take care of the various link types
    1030         524 :          DO ikind = 1, nlinks
    1031             :             CALL section_vals_val_get(qmmm_links, "LINK_TYPE", i_rep_section=ikind, &
    1032         396 :                                       i_val=link_type)
    1033         524 :             SELECT CASE (link_type)
    1034             :             CASE (do_qmmm_link_imomm)
    1035         388 :                num_qm_atom_tot = num_qm_atom_tot + 1
    1036         388 :                link_involv_mm = link_involv_mm + 1
    1037             :             CASE (do_qmmm_link_pseudo)
    1038           8 :                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         396 :                CPABORT("")
    1043             :             END SELECT
    1044             :          END DO
    1045             :       END IF
    1046         788 :       IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) &
    1047         186 :          ALLOCATE (mm_link_scale_factor(link_involv_mm))
    1048         788 :       IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) &
    1049         186 :          ALLOCATE (fist_scale_charge_link(link_involv_mm))
    1050         788 :       IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) &
    1051         372 :          ALLOCATE (mm_link_atoms(link_involv_mm))
    1052        2364 :       IF (PRESENT(qm_atom_index)) ALLOCATE (qm_atom_index(num_qm_atom_tot))
    1053        1576 :       IF (PRESENT(qm_atom_type)) ALLOCATE (qm_atom_type(num_qm_atom_tot))
    1054        6572 :       IF (PRESENT(qm_atom_index)) qm_atom_index = 0
    1055        3680 :       IF (PRESENT(qm_atom_type)) qm_atom_type = " "
    1056         788 :       num_qm_atom_tot = 1
    1057        2596 :       DO ikind = 1, nkind
    1058        1808 :          CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
    1059        6336 :          DO k = 1, n_var
    1060             :             CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
    1061        3740 :                                       i_vals=mm_indexes)
    1062        3740 :             IF (PRESENT(qm_atom_index)) THEN
    1063       14516 :                qm_atom_index(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = mm_indexes(:)
    1064             :             END IF
    1065        3740 :             IF (PRESENT(qm_atom_type)) THEN
    1066             :                CALL section_vals_val_get(qm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
    1067        1870 :                                          c_val=qm_atom_kind)
    1068        4564 :                qm_atom_type(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = qm_atom_kind
    1069             :             END IF
    1070        5548 :             num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
    1071             :          END DO
    1072             :       END DO
    1073         982 :       IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) mm_link_scale_factor = 0.0_dp
    1074         982 :       IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) fist_scale_charge_link = 0.0_dp
    1075        1176 :       IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) mm_link_atoms = 0
    1076         788 :       IF (explicit) THEN
    1077         524 :          DO ikind = 1, nlinks
    1078         396 :             IF (PRESENT(qm_atom_type)) THEN
    1079         198 :                CALL section_vals_val_get(qmmm_links, "QM_KIND", i_rep_section=ikind, c_val=qm_link_element)
    1080         396 :                qm_atom_type(num_qm_atom_tot:num_qm_atom_tot) = TRIM(qm_link_element)//"_LINK"
    1081             :             END IF
    1082         396 :             IF (PRESENT(qm_atom_index)) THEN
    1083         396 :                CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
    1084       23508 :                CPASSERT(ALL(qm_atom_index /= mm_index))
    1085         792 :                qm_atom_index(num_qm_atom_tot:num_qm_atom_tot) = mm_index
    1086         396 :                num_qm_atom_tot = num_qm_atom_tot + 1
    1087             :             END IF
    1088         396 :             IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) THEN
    1089         388 :                CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
    1090         388 :                mm_link_atoms(ikind) = mm_index
    1091             :             END IF
    1092         396 :             IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) THEN
    1093         194 :                CALL section_vals_val_get(qmmm_links, "QMMM_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
    1094         194 :                mm_link_scale_factor(ikind) = scale_f
    1095             :             END IF
    1096         524 :             IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) THEN
    1097         194 :                CALL section_vals_val_get(qmmm_links, "FIST_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
    1098         194 :                fist_scale_charge_link(ikind) = scale_f
    1099             :             END IF
    1100             :          END DO
    1101             :       END IF
    1102         788 :       CPASSERT(num_qm_atom_tot - 1 == SIZE(qm_atom_index))
    1103             : 
    1104         788 :    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         128 :    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          64 :       INTEGER, DIMENSION(:), POINTER                     :: wrk_tmp
    1130             :       REAL(KIND=dp)                                      :: alpha, my_radius
    1131             :       TYPE(section_vals_type), POINTER                   :: qmmm_link_section
    1132             : 
    1133          64 :       NULLIFY (wrk_tmp)
    1134          64 :       n_imomm = 0
    1135          64 :       n_gho = 0
    1136          64 :       n_pseudo = 0
    1137         128 :       qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
    1138          64 :       CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
    1139          64 :       CPASSERT(nlinks /= 0)
    1140         262 :       DO ikind = 1, nlinks
    1141         198 :          CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
    1142         198 :          IF (link_type == do_qmmm_link_imomm) n_imomm = n_imomm + 1
    1143         198 :          IF (link_type == do_qmmm_link_gho) n_gho = n_gho + 1
    1144         460 :          IF (link_type == do_qmmm_link_pseudo) n_pseudo = n_pseudo + 1
    1145             :       END DO
    1146          64 :       n_tot = n_imomm + n_gho + n_pseudo
    1147          64 :       CPASSERT(n_tot /= 0)
    1148          64 :       ALLOCATE (qmmm_links)
    1149             :       NULLIFY (qmmm_links%imomm, &
    1150          64 :                qmmm_links%pseudo)
    1151             :       ! IMOMM
    1152          64 :       IF (n_imomm /= 0) THEN
    1153         186 :          ALLOCATE (qmmm_links%imomm(n_imomm))
    1154         186 :          ALLOCATE (wrk_tmp(n_imomm))
    1155         256 :          DO ikind = 1, n_imomm
    1156         194 :             NULLIFY (qmmm_links%imomm(ikind)%link)
    1157         256 :             ALLOCATE (qmmm_links%imomm(ikind)%link)
    1158             :          END DO
    1159          62 :          n_imomm = 0
    1160         256 :          DO ikind = 1, nlinks
    1161         194 :             CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
    1162         256 :             IF (link_type == do_qmmm_link_imomm) THEN
    1163         194 :                n_imomm = n_imomm + 1
    1164         194 :                CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
    1165         194 :                CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
    1166         194 :                CALL section_vals_val_get(qmmm_link_section, "ALPHA_IMOMM", i_rep_section=ikind, r_val=alpha)
    1167         194 :                CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
    1168         194 :                qmmm_links%imomm(n_imomm)%link%qm_index = qm_index
    1169         194 :                qmmm_links%imomm(n_imomm)%link%mm_index = mm_index
    1170         194 :                qmmm_links%imomm(n_imomm)%link%alpha = alpha
    1171         194 :                wrk_tmp(n_imomm) = mm_index
    1172         194 :                IF (n_rep_val == 1) THEN
    1173          64 :                   CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, r_val=my_radius)
    1174         996 :                   WHERE (mm_atom_index == mm_index) mm_el_pot_radius = my_radius
    1175         996 :                   WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
    1176             :                END IF
    1177         194 :                CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
    1178         194 :                IF (n_rep_val == 1) THEN
    1179           0 :                   CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, r_val=my_radius)
    1180           0 :                   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         256 :          DO ikind = 1, SIZE(wrk_tmp)
    1188        3514 :             IF (COUNT(wrk_tmp == wrk_tmp(ikind)) > 1) THEN
    1189           0 :                WRITE (iw, '(/A)') "In the IMOMM scheme no more than one QM atom can be bounded to the same MM atom."
    1190           0 :                WRITE (iw, '(A)') "Multiple link MM atom not allowed. Check your link sections."
    1191           0 :                CPABORT("")
    1192             :             END IF
    1193             :          END DO
    1194          62 :          DEALLOCATE (wrk_tmp)
    1195             :       END IF
    1196             :       ! PSEUDO
    1197          64 :       IF (n_pseudo /= 0) THEN
    1198           6 :          ALLOCATE (qmmm_links%pseudo(n_pseudo))
    1199           6 :          DO ikind = 1, n_pseudo
    1200           4 :             NULLIFY (qmmm_links%pseudo(ikind)%link)
    1201           6 :             ALLOCATE (qmmm_links%pseudo(ikind)%link)
    1202             :          END DO
    1203           2 :          n_pseudo = 0
    1204           6 :          DO ikind = 1, nlinks
    1205           4 :             CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
    1206           6 :             IF (link_type == do_qmmm_link_pseudo) THEN
    1207           4 :                n_pseudo = n_pseudo + 1
    1208           4 :                CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
    1209           4 :                CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
    1210           4 :                qmmm_links%pseudo(n_pseudo)%link%qm_index = qm_index
    1211           4 :                qmmm_links%pseudo(n_pseudo)%link%mm_index = mm_index
    1212             :             END IF
    1213             :          END DO
    1214             :       END IF
    1215             :       ! GHO
    1216          64 :       IF (n_gho /= 0) THEN
    1217             :          ! not yet implemented
    1218             :          ! still to define : type, implementation into QS
    1219           0 :          CPABORT("")
    1220             :       END IF
    1221          64 :    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         128 :    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          64 :       explicit = .FALSE.
    1256          64 :       move_mm_charges = .FALSE.
    1257          64 :       add_mm_charges = .FALSE.
    1258          64 :       NULLIFY (qmmm_link_section, move_section, add_section)
    1259          64 :       qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
    1260          64 :       CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
    1261          64 :       CPASSERT(nlinks /= 0)
    1262             :       icount = 0
    1263          64 :       n_move_tot = 0
    1264          64 :       n_add_tot = 0
    1265         262 :       DO ikind = 1, nlinks
    1266             :          move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
    1267         198 :                                                     i_rep_section=ikind)
    1268         198 :          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         198 :                                                    i_rep_section=ikind)
    1271         198 :          CALL section_vals_get(add_section, n_repetition=n_adds)
    1272         198 :          n_move_tot = n_move_tot + n_moves
    1273         460 :          n_add_tot = n_add_tot + n_adds
    1274             :       END DO
    1275          64 :       icount = n_move_tot + n_add_tot
    1276          64 :       IF (n_add_tot /= 0) add_mm_charges = .TRUE.
    1277          64 :       IF (n_move_tot /= 0) move_mm_charges = .TRUE.
    1278             :       !
    1279             :       ! create add_set_type
    1280             :       !
    1281          64 :       CALL create_add_set_type(added_charges, ndim=icount)
    1282             :       !
    1283             :       ! Fill in structures
    1284             :       !
    1285          64 :       icount = 0
    1286         262 :       DO ikind = 1, nlinks
    1287             :          move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
    1288         198 :                                                     i_rep_section=ikind)
    1289         198 :          CALL section_vals_get(move_section, explicit=explicit, n_repetition=n_moves)
    1290             :          !
    1291             :          ! Moving charge atoms
    1292             :          !
    1293         198 :          IF (explicit) THEN
    1294          36 :             DO i_add = 1, n_moves
    1295          26 :                icount = icount + 1
    1296          26 :                CALL section_vals_val_get(move_section, "ATOM_INDEX_1", i_val=Index1, i_rep_section=i_add)
    1297          26 :                CALL section_vals_val_get(move_section, "ATOM_INDEX_2", i_val=Index2, i_rep_section=i_add)
    1298          26 :                CALL section_vals_val_get(move_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
    1299          26 :                CALL section_vals_val_get(move_section, "RADIUS", r_val=radius, i_rep_section=i_add)
    1300          26 :                CALL section_vals_val_get(move_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
    1301          26 :                c_radius = radius
    1302          26 :                IF (n_rep_val == 1) &
    1303          26 :                   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          62 :                                      mm_atom_index=mm_atom_index, move=n_moves, Ind1=ind1)
    1309             :             END DO
    1310          10 :             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         198 :                                                    i_rep_section=ikind)
    1315         198 :          CALL section_vals_get(add_section, explicit=explicit, n_repetition=n_adds)
    1316             :          !
    1317             :          ! Adding charge atoms
    1318             :          !
    1319         460 :          IF (explicit) THEN
    1320           4 :             DO i_add = 1, n_adds
    1321           2 :                icount = icount + 1
    1322           2 :                CALL section_vals_val_get(add_section, "ATOM_INDEX_1", i_val=Index1, i_rep_section=i_add)
    1323           2 :                CALL section_vals_val_get(add_section, "ATOM_INDEX_2", i_val=Index2, i_rep_section=i_add)
    1324           2 :                CALL section_vals_val_get(add_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
    1325           2 :                CALL section_vals_val_get(add_section, "RADIUS", r_val=radius, i_rep_section=i_add)
    1326           2 :                CALL section_vals_val_get(add_section, "CHARGE", r_val=charge, i_rep_section=i_add)
    1327           2 :                CALL section_vals_val_get(add_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
    1328           2 :                c_radius = radius
    1329           2 :                IF (n_rep_val == 1) &
    1330           2 :                   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           6 :                                      mm_atom_index=mm_atom_index)
    1336             :             END DO
    1337             :          END IF
    1338             :       END DO
    1339             : 
    1340          64 :    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          28 :    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          28 :       my_move = 0
    1378          28 :       my_radius = radius
    1379          28 :       my_c_radius = c_radius
    1380          28 :       IF (PRESENT(charge)) my_charge = charge
    1381          28 :       IF (PRESENT(move)) my_move = move
    1382          28 :       i = 1
    1383          60 :       GetId: DO WHILE (i <= SIZE(mm_atom_index))
    1384          60 :          IF (Index1 == mm_atom_index(i)) EXIT GetId
    1385          60 :          i = i + 1
    1386             :       END DO GetId
    1387          28 :       IF (PRESENT(ind1)) ind1 = i
    1388          28 :       CPASSERT(i <= SIZE(mm_atom_index))
    1389          28 :       IF (.NOT. PRESENT(charge)) my_charge = mm_atom_chrg(i)/REAL(my_move, KIND=dp)
    1390          28 :       IF (my_radius == 0.0_dp) my_radius = mm_el_pot_radius(i)
    1391          28 :       IF (my_c_radius == 0.0_dp) my_c_radius = mm_el_pot_radius_corr(i)
    1392             : 
    1393          28 :       added_charges%add_env(icount)%Index1 = Index1
    1394          28 :       added_charges%add_env(icount)%Index2 = Index2
    1395          28 :       added_charges%add_env(icount)%alpha = alpha
    1396          28 :       added_charges%mm_atom_index(icount) = icount
    1397          28 :       added_charges%mm_atom_chrg(icount) = my_charge
    1398          28 :       added_charges%mm_el_pot_radius(icount) = my_radius
    1399          28 :       added_charges%mm_el_pot_radius_corr(icount) = my_c_radius
    1400          28 :    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         788 :    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         394 :       REAL(KINd=dp), DIMENSION(:), POINTER               :: vec
    1424             : 
    1425             : ! This is the vector that corrects position to apply properly the PBC
    1426             : 
    1427         394 :       tmp(1) = qm_cell_small%hmat(1, 1)
    1428         394 :       tmp(2) = qm_cell_small%hmat(2, 2)
    1429         394 :       tmp(3) = qm_cell_small%hmat(3, 3)
    1430        1576 :       CPASSERT(ALL(tmp > 0))
    1431        1576 :       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         394 :       CALL section_vals_val_get(qmmm_section, "CENTER_GRID", l_val=center_grid)
    1435         394 :       IF (center_grid) THEN
    1436          72 :          qmmm_env%utrasl = dr
    1437             :       ELSE
    1438        1504 :          qmmm_env%utrasl = 1.0_dp
    1439             :       END IF
    1440         394 :       CALL section_vals_val_get(qmmm_section, "INITIAL_TRANSLATION_VECTOR", r_vals=vec)
    1441        3152 :       qmmm_env%transl_v = vec
    1442         394 :    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          10 :    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          10 :       INTEGER, DIMENSION(:), POINTER                     :: mm_indexes
    1466             :       LOGICAL                                            :: fix_xyz, imageind_in_range
    1467          10 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind
    1468             : 
    1469          10 :       NULLIFY (mm_indexes, molecule_kind)
    1470          10 :       imageind_in_range = .FALSE.
    1471          10 :       num_image_mm_atom = 0
    1472          10 :       max_index = 0
    1473             : 
    1474             :       CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
    1475          10 :                                 n_rep_val=n_var)
    1476          20 :       DO i = 1, n_var
    1477             :          CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
    1478          10 :                                    i_rep_val=i, i_vals=mm_indexes)
    1479          20 :          num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
    1480             :       END DO
    1481             : 
    1482          30 :       ALLOCATE (qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom))
    1483             : 
    1484          30 :       qmmm_env%image_charge_pot%image_mm_list = 0
    1485          10 :       num_image_mm_atom = 1
    1486             : 
    1487          20 :       DO i = 1, n_var
    1488             :          CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
    1489          10 :                                    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          60 :                                                  + SIZE(mm_indexes) - 1) = mm_indexes(:)
    1492          20 :          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          10 :       num_image_mm_atom = num_image_mm_atom - 1
    1497             : 
    1498          10 :       max_index = SIZE(subsys_mm%particles%els)
    1499             : 
    1500          10 :       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          60 :                           .AND. (MINVAL(qmmm_env%image_charge_pot%image_mm_list) > 0)
    1503          10 :       CPASSERT(imageind_in_range)
    1504             : 
    1505          30 :       DO i = 1, num_image_mm_atom
    1506          20 :          atom_a = qmmm_env%image_charge_pot%image_mm_list(i)
    1507          60 :          IF (ANY(qm_atom_index == atom_a)) THEN
    1508           0 :             CPABORT("Image atom list must only contain MM atoms")
    1509             :          END IF
    1510          40 :          DO j = i + 1, num_image_mm_atom
    1511          10 :             atom_b = qmmm_env%image_charge_pot%image_mm_list(j)
    1512          10 :             IF (atom_a == atom_b) &
    1513          20 :                CPABORT("There are atoms doubled in image list.")
    1514             :          END DO
    1515             :       END DO
    1516             : 
    1517             :       ! check if molecules in list carry constraints
    1518          10 :       num_const_atom = 0
    1519          10 :       fix_xyz = .TRUE.
    1520          10 :       IF (ASSOCIATED(subsys_mm%molecule_kinds)) THEN
    1521          10 :          IF (ASSOCIATED(subsys_mm%molecule_kinds%els)) THEN
    1522          10 :             molecule_kind => subsys_mm%molecule_kinds%els
    1523          78 :             DO i = 1, SIZE(molecule_kind)
    1524          76 :                IF (.NOT. ASSOCIATED(molecule_kind(i)%fixd_list)) EXIT
    1525          68 :                IF (.NOT. fix_xyz) EXIT
    1526          82 :                DO j = 1, SIZE(molecule_kind(i)%fixd_list)
    1527           4 :                   IF (.NOT. fix_xyz) EXIT
    1528          80 :                   DO k = 1, num_image_mm_atom
    1529           8 :                      atom_a = qmmm_env%image_charge_pot%image_mm_list(k)
    1530          12 :                      IF (atom_a == molecule_kind(i)%fixd_list(j)%fixd) THEN
    1531           4 :                         num_const_atom = num_const_atom + 1
    1532           4 :                         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          10 :       IF (num_const_atom == num_image_mm_atom .AND. fix_xyz) THEN
    1546           2 :          qmmm_env%image_charge_pot%state_image_matrix = calc_once
    1547             :       ELSE
    1548           8 :          qmmm_env%image_charge_pot%state_image_matrix = calc_always
    1549             :       END IF
    1550             : 
    1551          10 :    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          10 :    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          10 :       logger => cp_get_default_logger()
    1571             :       iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%PROGRAM_RUN_INFO", &
    1572          10 :                                 extension=".log")
    1573          10 :       eta = qmmm_env%image_charge_pot%eta
    1574          10 :       eta_conv = cp_unit_from_cp2k(eta, "angstrom", power=-2)
    1575          10 :       V0 = qmmm_env%image_charge_pot%V0
    1576          10 :       V0_conv = cp_unit_from_cp2k(V0, "volt")
    1577             : 
    1578          10 :       IF (iw > 0) THEN
    1579           5 :          WRITE (iw, FMT="(T25,A)") "IMAGE CHARGE PARAMETERS"
    1580           5 :          WRITE (iw, FMT="(T25,A)") REPEAT("-", 23)
    1581           5 :          WRITE (iw, FMT="(/)")
    1582           5 :          WRITE (iw, FMT="(T2,A)") "INDEX OF MM ATOMS CARRYING AN IMAGE CHARGE:"
    1583           5 :          WRITE (iw, FMT="(/)")
    1584             : 
    1585          15 :          WRITE (iw, "(7X,10I6)") qmmm_env%image_charge_pot%image_mm_list
    1586           5 :          WRITE (iw, FMT="(/)")
    1587             :          WRITE (iw, "(T2,A52,T69,F12.8)") &
    1588           5 :             "WIDTH OF GAUSSIAN CHARGE DISTRIBUTION [angstrom^-2]:", eta_conv
    1589           5 :          WRITE (iw, "(T2,A26,T69,F12.8)") "EXTERNAL POTENTIAL [volt]:", V0_conv
    1590           5 :          WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 79)
    1591             :       END IF
    1592             :       CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
    1593          10 :                                         "PRINT%PROGRAM_RUN_INFO")
    1594             : 
    1595          10 :    END SUBROUTINE print_image_charge_info
    1596             : 
    1597             : END MODULE qmmm_init

Generated by: LCOV version 1.15