LCOV - code coverage report
Current view: top level - src - qmmm_image_charge.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 90.0 % 418 376
Test Date: 2025-07-25 12:55:17 Functions: 94.1 % 17 16

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Routines for image charge calculation within QM/MM
      10              : !> \par History
      11              : !>      12.2011 created
      12              : !> \author Dorothea Golze
      13              : ! **************************************************************************************************
      14              : MODULE qmmm_image_charge
      15              :    USE ao_util,                         ONLY: exp_radius_very_extended
      16              :    USE cell_types,                      ONLY: cell_type,&
      17              :                                               pbc
      18              :    USE cp_control_types,                ONLY: dft_control_type
      19              :    USE cp_eri_mme_interface,            ONLY: cp_eri_mme_param,&
      20              :                                               cp_eri_mme_update_local_counts
      21              :    USE cp_files,                        ONLY: close_file,&
      22              :                                               open_file
      23              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      24              :                                               cp_logger_type
      25              :    USE cp_output_handling,              ONLY: cp_p_file,&
      26              :                                               cp_print_key_finished_output,&
      27              :                                               cp_print_key_generate_filename,&
      28              :                                               cp_print_key_should_output,&
      29              :                                               cp_print_key_unit_nr
      30              :    USE eri_mme_integrate,               ONLY: eri_mme_2c_integrate
      31              :    USE input_constants,                 ONLY: calc_always,&
      32              :                                               calc_once,&
      33              :                                               calc_once_done,&
      34              :                                               do_eri_gpw,&
      35              :                                               do_eri_mme
      36              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      37              :                                               section_vals_type,&
      38              :                                               section_vals_val_get
      39              :    USE kinds,                           ONLY: default_path_length,&
      40              :                                               dp
      41              :    USE mathconstants,                   ONLY: pi
      42              :    USE mathlib,                         ONLY: invmat_symm
      43              :    USE memory_utilities,                ONLY: reallocate
      44              :    USE message_passing,                 ONLY: mp_para_env_type
      45              :    USE pw_env_types,                    ONLY: pw_env_get,&
      46              :                                               pw_env_type
      47              :    USE pw_methods,                      ONLY: pw_axpy,&
      48              :                                               pw_integral_ab,&
      49              :                                               pw_scale,&
      50              :                                               pw_transfer,&
      51              :                                               pw_zero
      52              :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
      53              :    USE pw_poisson_types,                ONLY: pw_poisson_type
      54              :    USE pw_pool_types,                   ONLY: pw_pool_type
      55              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      56              :                                               pw_r3d_rs_type
      57              :    USE qmmm_types_low,                  ONLY: qmmm_env_qm_type
      58              :    USE qs_collocate_density,            ONLY: calculate_rho_metal,&
      59              :                                               calculate_rho_single_gaussian
      60              :    USE qs_energy_types,                 ONLY: qs_energy_type
      61              :    USE qs_environment_types,            ONLY: get_qs_env,&
      62              :                                               qs_environment_type
      63              :    USE qs_integrate_potential,          ONLY: integrate_pgf_product
      64              :    USE realspace_grid_types,            ONLY: realspace_grid_desc_type,&
      65              :                                               realspace_grid_type,&
      66              :                                               transfer_pw2rs
      67              :    USE util,                            ONLY: get_limit
      68              :    USE virial_types,                    ONLY: virial_type
      69              : #include "./base/base_uses.f90"
      70              : 
      71              :    IMPLICIT NONE
      72              :    PRIVATE
      73              : 
      74              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      75              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_image_charge'
      76              : 
      77              :    PUBLIC :: calculate_image_pot, &
      78              :              integrate_potential_devga_rspace, &
      79              :              conditional_calc_image_matrix, &
      80              :              add_image_pot_to_hartree_pot, &
      81              :              print_image_coefficients
      82              : 
      83              : !***
      84              : CONTAINS
      85              : ! **************************************************************************************************
      86              : !> \brief determines coefficients by solving image_matrix*coeff=-pot_const by
      87              : !>        Gaussian elimination or in an iterative fashion and calculates
      88              : !>        image/metal potential with these coefficients
      89              : !> \param v_hartree_rspace Hartree potential in real space
      90              : !> \param rho_hartree_gspace Kohn Sham density in reciprocal space
      91              : !> \param energy structure where energies are stored
      92              : !> \param qmmm_env qmmm environment
      93              : !> \param qs_env qs environment
      94              : ! **************************************************************************************************
      95           60 :    SUBROUTINE calculate_image_pot(v_hartree_rspace, rho_hartree_gspace, energy, &
      96              :                                   qmmm_env, qs_env)
      97              : 
      98              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_hartree_rspace
      99              :       TYPE(pw_c1d_gs_type), INTENT(IN)                   :: rho_hartree_gspace
     100              :       TYPE(qs_energy_type), POINTER                      :: energy
     101              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     102              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     103              : 
     104              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_image_pot'
     105              : 
     106              :       INTEGER                                            :: handle
     107              : 
     108           60 :       CALL timeset(routineN, handle)
     109              : 
     110           60 :       IF (qmmm_env%image_charge_pot%coeff_iterative) THEN
     111              :          !calculate preconditioner matrix for CG if necessary
     112           12 :          IF (qs_env%calc_image_preconditioner) THEN
     113            2 :             IF (qmmm_env%image_charge_pot%image_restart) THEN
     114              :                CALL restart_image_matrix(image_matrix=qs_env%image_matrix, &
     115            0 :                                          qs_env=qs_env, qmmm_env=qmmm_env)
     116              :             ELSE
     117              :                CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
     118            2 :                                            qs_env=qs_env, qmmm_env=qmmm_env)
     119              :             END IF
     120              :          END IF
     121              :          CALL calc_image_coeff_iterative(v_hartree_rspace=v_hartree_rspace, &
     122              :                                          coeff=qs_env%image_coeff, qmmm_env=qmmm_env, &
     123           12 :                                          qs_env=qs_env)
     124              : 
     125              :       ELSE
     126              :          CALL calc_image_coeff_gaussalgorithm(v_hartree_rspace=v_hartree_rspace, &
     127              :                                               coeff=qs_env%image_coeff, qmmm_env=qmmm_env, &
     128           48 :                                               qs_env=qs_env)
     129              :       END IF
     130              : 
     131              :       ! calculate the image/metal potential with the optimized coefficients
     132           60 :       ALLOCATE (qs_env%ks_qmmm_env%v_metal_rspace)
     133              :       CALL calculate_potential_metal(v_metal_rspace= &
     134              :                                      qs_env%ks_qmmm_env%v_metal_rspace, coeff=qs_env%image_coeff, &
     135              :                                      rho_hartree_gspace=rho_hartree_gspace, &
     136           60 :                                      energy=energy, qs_env=qs_env)
     137              : 
     138           60 :       CALL timestop(handle)
     139              : 
     140           60 :    END SUBROUTINE calculate_image_pot
     141              : 
     142              : ! **************************************************************************************************
     143              : !> \brief determines coefficients by solving the linear set of equations
     144              : !>        image_matrix*coeff=-pot_const using a Gaussian elimination scheme
     145              : !> \param v_hartree_rspace Hartree potential in real space
     146              : !> \param coeff expansion coefficients of the image charge density, i.e.
     147              : !>        rho_metal=sum_a c_a*g_a
     148              : !> \param qmmm_env qmmm environment
     149              : !> \param qs_env qs environment
     150              : ! **************************************************************************************************
     151           48 :    SUBROUTINE calc_image_coeff_gaussalgorithm(v_hartree_rspace, coeff, qmmm_env, &
     152              :                                               qs_env)
     153              : 
     154              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_hartree_rspace
     155              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: coeff
     156              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     157              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     158              : 
     159              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_image_coeff_gaussalgorithm'
     160              : 
     161              :       INTEGER                                            :: handle, info, natom
     162              :       REAL(KIND=dp)                                      :: eta, V0
     163              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pot_const
     164              : 
     165           48 :       CALL timeset(routineN, handle)
     166              : 
     167              :       NULLIFY (pot_const)
     168              : 
     169              :       !minus sign V0: account for the fact that v_hartree has the opposite sign
     170           48 :       V0 = -qmmm_env%image_charge_pot%V0
     171           48 :       eta = qmmm_env%image_charge_pot%eta
     172           48 :       natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
     173              : 
     174          144 :       ALLOCATE (pot_const(natom))
     175           48 :       IF (.NOT. ASSOCIATED(coeff)) THEN
     176           16 :          ALLOCATE (coeff(natom))
     177              :       END IF
     178          144 :       coeff = 0.0_dp
     179              : 
     180              :       CALL integrate_potential_ga_rspace(v_hartree_rspace, qmmm_env, qs_env, &
     181           48 :                                          pot_const)
     182              :       !add integral V0*ga(r)
     183          144 :       pot_const(:) = -pot_const(:) + V0*SQRT((pi/eta)**3)
     184              : 
     185              :       !solve linear system of equations T*coeff=-pot_const
     186              :       !LU factorization of T by DGETRF done in calculate_image_matrix
     187              :       CALL dgetrs('N', natom, 1, qs_env%image_matrix, natom, qs_env%ipiv, &
     188           48 :                   pot_const, natom, info)
     189           48 :       CPASSERT(info == 0)
     190              : 
     191          240 :       coeff = pot_const
     192              : 
     193           48 :       DEALLOCATE (pot_const)
     194              : 
     195           48 :       CALL timestop(handle)
     196              : 
     197           48 :    END SUBROUTINE calc_image_coeff_gaussalgorithm
     198              : 
     199              : ! **************************************************************************************************
     200              : !> \brief determines image coefficients iteratively
     201              : !> \param v_hartree_rspace Hartree potential in real space
     202              : !> \param coeff expansion coefficients of the image charge density, i.e.
     203              : !>        rho_metal=sum_a c_a*g_a
     204              : !> \param qmmm_env qmmm environment
     205              : !> \param qs_env qs environment
     206              : ! **************************************************************************************************
     207           12 :    SUBROUTINE calc_image_coeff_iterative(v_hartree_rspace, coeff, qmmm_env, &
     208              :                                          qs_env)
     209              : 
     210              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_hartree_rspace
     211              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: coeff
     212              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     213              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     214              : 
     215              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_image_coeff_iterative'
     216              : 
     217              :       INTEGER                                            :: handle, iter_steps, natom, output_unit
     218              :       REAL(KIND=dp)                                      :: alpha, eta, rsnew, rsold, V0
     219           12 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: Ad, d, pot_const, r, vmetal_const, z
     220              :       TYPE(cp_logger_type), POINTER                      :: logger
     221              :       TYPE(pw_r3d_rs_type)                               :: auxpot_Ad_rspace, v_metal_rspace_guess
     222              :       TYPE(section_vals_type), POINTER                   :: input
     223              : 
     224           12 :       CALL timeset(routineN, handle)
     225              : 
     226           12 :       NULLIFY (pot_const, vmetal_const, logger, input)
     227           12 :       logger => cp_get_default_logger()
     228              : 
     229              :       !minus sign V0: account for the fact that v_hartree has the opposite sign
     230           12 :       V0 = -qmmm_env%image_charge_pot%V0
     231           12 :       eta = qmmm_env%image_charge_pot%eta
     232           12 :       natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
     233              : 
     234           36 :       ALLOCATE (pot_const(natom))
     235           24 :       ALLOCATE (vmetal_const(natom))
     236           24 :       ALLOCATE (r(natom))
     237           24 :       ALLOCATE (d(natom))
     238           24 :       ALLOCATE (z(natom))
     239           24 :       ALLOCATE (Ad(natom))
     240           12 :       IF (.NOT. ASSOCIATED(coeff)) THEN
     241            4 :          ALLOCATE (coeff(natom))
     242              :       END IF
     243              : 
     244              :       CALL integrate_potential_ga_rspace(v_hartree_rspace, qmmm_env, qs_env, &
     245           12 :                                          pot_const)
     246              : 
     247              :       !add integral V0*ga(r)
     248           36 :       pot_const(:) = -pot_const(:) + V0*SQRT((pi/eta)**3)
     249              : 
     250              :       !initial guess for coeff
     251           36 :       coeff = 1.0_dp
     252           36 :       d = 0.0_dp
     253           36 :       z = 0.0_dp
     254           36 :       r = 0.0_dp
     255           12 :       rsold = 0.0_dp
     256           12 :       rsnew = 0.0_dp
     257           12 :       iter_steps = 0
     258              : 
     259              :       !calculate first guess of image/metal potential
     260              :       CALL calculate_potential_metal(v_metal_rspace=v_metal_rspace_guess, &
     261           12 :                                      coeff=coeff, qs_env=qs_env)
     262              :       CALL integrate_potential_ga_rspace(potential=v_metal_rspace_guess, &
     263           12 :                                          qmmm_env=qmmm_env, qs_env=qs_env, int_res=vmetal_const)
     264              : 
     265              :       ! modify coefficients iteratively
     266           60 :       r = pot_const - vmetal_const
     267          276 :       z = MATMUL(qs_env%image_matrix, r)
     268           60 :       d = z
     269           36 :       rsold = DOT_PRODUCT(r, z)
     270              : 
     271            6 :       DO
     272              :          !calculate A*d
     273           54 :          Ad = 0.0_dp
     274              :          CALL calculate_potential_metal(v_metal_rspace= &
     275           18 :                                         auxpot_Ad_rspace, coeff=d, qs_env=qs_env)
     276              :          CALL integrate_potential_ga_rspace(potential= &
     277              :                                             auxpot_Ad_rspace, qmmm_env=qmmm_env, &
     278           18 :                                             qs_env=qs_env, int_res=Ad)
     279              : 
     280           54 :          alpha = rsold/DOT_PRODUCT(d, Ad)
     281           90 :          coeff = coeff + alpha*d
     282              : 
     283           90 :          r = r - alpha*Ad
     284          414 :          z = MATMUL(qs_env%image_matrix, r)
     285           54 :          rsnew = DOT_PRODUCT(r, z)
     286           18 :          iter_steps = iter_steps + 1
     287              :          ! SQRT(rsnew) < 1.0E-08
     288           18 :          IF (rsnew < 1.0E-16) THEN
     289           12 :             CALL auxpot_Ad_rspace%release()
     290              :             EXIT
     291              :          END IF
     292           30 :          d = z + rsnew/rsold*d
     293            6 :          rsold = rsnew
     294           24 :          CALL auxpot_Ad_rspace%release()
     295              :       END DO
     296              : 
     297              :       ! print iteration info
     298              :       CALL get_qs_env(qs_env=qs_env, &
     299           12 :                       input=input)
     300              :       output_unit = cp_print_key_unit_nr(logger, input, &
     301              :                                          "QMMM%PRINT%PROGRAM_RUN_INFO", &
     302           12 :                                          extension=".qmmmLog")
     303           12 :       IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A,T74,I7)") &
     304            6 :          "Number of iteration steps for determination of image coefficients:", iter_steps
     305              :       CALL cp_print_key_finished_output(output_unit, logger, input, &
     306           12 :                                         "QMMM%PRINT%PROGRAM_RUN_INFO")
     307              : 
     308           12 :       IF (iter_steps .LT. 25) THEN
     309           12 :          qs_env%calc_image_preconditioner = .FALSE.
     310              :       ELSE
     311            0 :          qs_env%calc_image_preconditioner = .TRUE.
     312              :       END IF
     313              : 
     314           12 :       CALL v_metal_rspace_guess%release()
     315           12 :       DEALLOCATE (pot_const)
     316           12 :       DEALLOCATE (vmetal_const)
     317           12 :       DEALLOCATE (r)
     318           12 :       DEALLOCATE (d, z)
     319           12 :       DEALLOCATE (Ad)
     320              : 
     321           12 :       CALL timestop(handle)
     322              : 
     323           24 :    END SUBROUTINE calc_image_coeff_iterative
     324              : 
     325              : ! ****************************************************************************
     326              : !> \brief calculates the integral V(r)*ga(r)
     327              : !> \param potential any potential
     328              : !> \param qmmm_env qmmm environment
     329              : !> \param qs_env qs environment
     330              : !> \param int_res result of the integration
     331              : !> \param atom_num atom index, needed when calculating image_matrix
     332              : !> \param atom_num_ref index of reference atom, needed when calculating
     333              : !>        image_matrix
     334              : ! **************************************************************************************************
     335           98 :    SUBROUTINE integrate_potential_ga_rspace(potential, qmmm_env, qs_env, int_res, &
     336              :                                             atom_num, atom_num_ref)
     337              : 
     338              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: potential
     339              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     340              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     341              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: int_res
     342              :       INTEGER, INTENT(IN), OPTIONAL                      :: atom_num, atom_num_ref
     343              : 
     344              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_potential_ga_rspace'
     345              : 
     346              :       INTEGER                                            :: atom_a, atom_b, atom_ref, handle, iatom, &
     347              :                                                             j, k, natom, npme
     348           98 :       INTEGER, DIMENSION(:), POINTER                     :: cores
     349              :       REAL(KIND=dp)                                      :: eps_rho_rspace, radius
     350              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     351           98 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hab
     352              :       TYPE(cell_type), POINTER                           :: cell
     353              :       TYPE(dft_control_type), POINTER                    :: dft_control
     354              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     355              :       TYPE(pw_env_type), POINTER                         :: pw_env
     356              :       TYPE(realspace_grid_desc_type), POINTER            :: auxbas_rs_desc
     357              :       TYPE(realspace_grid_type), POINTER                 :: rs_v
     358              : 
     359           98 :       CALL timeset(routineN, handle)
     360              : 
     361           98 :       NULLIFY (cores, hab, cell, auxbas_rs_desc, pw_env, para_env, &
     362           98 :                dft_control, rs_v)
     363           98 :       ALLOCATE (hab(1, 1))
     364              : 
     365           98 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     366              :       CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, &
     367           98 :                       auxbas_rs_grid=rs_v)
     368           98 :       CALL transfer_pw2rs(rs_v, potential)
     369              : 
     370              :       CALL get_qs_env(qs_env=qs_env, &
     371              :                       cell=cell, &
     372              :                       dft_control=dft_control, &
     373           98 :                       para_env=para_env, pw_env=pw_env)
     374              : 
     375           98 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     376              : 
     377           98 :       natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
     378           98 :       k = 1
     379           98 :       IF (PRESENT(atom_num)) k = atom_num
     380              : 
     381           98 :       CALL reallocate(cores, 1, natom - k + 1)
     382          294 :       int_res = 0.0_dp
     383           98 :       npme = 0
     384          290 :       cores = 0
     385              : 
     386          290 :       DO iatom = k, natom
     387          290 :          IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
     388              :             ! replicated realspace grid, split the atoms up between procs
     389          192 :             IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
     390           96 :                npme = npme + 1
     391           96 :                cores(npme) = iatom
     392              :             END IF
     393              :          ELSE
     394            0 :             npme = npme + 1
     395            0 :             cores(npme) = iatom
     396              :          END IF
     397              :       END DO
     398              : 
     399          194 :       DO j = 1, npme
     400              : 
     401           96 :          iatom = cores(j)
     402           96 :          atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
     403              : 
     404           96 :          IF (PRESENT(atom_num) .AND. PRESENT(atom_num_ref)) THEN
     405              :             ! shift the function since potential only calculate for ref atom
     406            6 :             atom_b = qmmm_env%image_charge_pot%image_mm_list(k)
     407            6 :             atom_ref = qmmm_env%image_charge_pot%image_mm_list(atom_num_ref)
     408              :             ra(:) = pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell) &
     409              :                     - pbc(qmmm_env%image_charge_pot%particles_all(atom_b)%r, cell) &
     410           24 :                     + pbc(qmmm_env%image_charge_pot%particles_all(atom_ref)%r, cell)
     411              : 
     412              :          ELSE
     413           90 :             ra(:) = pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell)
     414              :          END IF
     415              : 
     416           96 :          hab(1, 1) = 0.0_dp
     417              : 
     418              :          radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
     419              :                                            ra=ra, rb=ra, rp=ra, &
     420              :                                            zetp=qmmm_env%image_charge_pot%eta, eps=eps_rho_rspace, &
     421           96 :                                            prefactor=1.0_dp, cutoff=1.0_dp)
     422              : 
     423              :          CALL integrate_pgf_product(0, qmmm_env%image_charge_pot%eta, 0, &
     424              :                                     0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
     425              :                                     rs_v, hab, o1=0, o2=0, &
     426              :                                     radius=radius, calculate_forces=.FALSE., &
     427           96 :                                     use_subpatch=.TRUE., subpatch_pattern=0)
     428              : 
     429          194 :          int_res(iatom) = hab(1, 1)
     430              : 
     431              :       END DO
     432              : 
     433          490 :       CALL para_env%sum(int_res)
     434              : 
     435           98 :       DEALLOCATE (hab, cores)
     436              : 
     437           98 :       CALL timestop(handle)
     438              : 
     439           98 :    END SUBROUTINE integrate_potential_ga_rspace
     440              : 
     441              : ! **************************************************************************************************
     442              : !> \brief calculates the image forces on the MM atoms
     443              : !> \param potential any potential, in this case: Hartree potential
     444              : !> \param coeff expansion coefficients of the image charge density, i.e.
     445              : !>        rho_metal=sum_a c_a*g_a
     446              : !> \param forces structure storing the force contribution of the image charges
     447              : !>        for the metal (MM) atoms
     448              : !> \param qmmm_env qmmm environment
     449              : !> \param qs_env qs environment
     450              : ! **************************************************************************************************
     451           20 :    SUBROUTINE integrate_potential_devga_rspace(potential, coeff, forces, qmmm_env, &
     452              :                                                qs_env)
     453              : 
     454              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: potential
     455              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: coeff
     456              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: forces
     457              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     458              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     459              : 
     460              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_potential_devga_rspace'
     461              : 
     462              :       INTEGER                                            :: atom_a, handle, iatom, j, natom, npme
     463           20 :       INTEGER, DIMENSION(:), POINTER                     :: cores
     464              :       LOGICAL                                            :: use_virial
     465              :       REAL(KIND=dp)                                      :: eps_rho_rspace, radius
     466              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, ra
     467           20 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hab, pab
     468              :       TYPE(cell_type), POINTER                           :: cell
     469              :       TYPE(dft_control_type), POINTER                    :: dft_control
     470              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     471              :       TYPE(pw_env_type), POINTER                         :: pw_env
     472              :       TYPE(realspace_grid_desc_type), POINTER            :: auxbas_rs_desc
     473              :       TYPE(realspace_grid_type), POINTER                 :: rs_v
     474              :       TYPE(virial_type), POINTER                         :: virial
     475              : 
     476           20 :       CALL timeset(routineN, handle)
     477              : 
     478           20 :       NULLIFY (cores, hab, pab, cell, auxbas_rs_desc, pw_env, para_env, &
     479           20 :                dft_control, rs_v, virial)
     480           20 :       use_virial = .FALSE.
     481              : 
     482           20 :       ALLOCATE (hab(1, 1))
     483           20 :       ALLOCATE (pab(1, 1))
     484              : 
     485           20 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     486              :       CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, &
     487           20 :                       auxbas_rs_grid=rs_v)
     488           20 :       CALL transfer_pw2rs(rs_v, potential)
     489              : 
     490              :       CALL get_qs_env(qs_env=qs_env, &
     491              :                       cell=cell, &
     492              :                       dft_control=dft_control, &
     493              :                       para_env=para_env, pw_env=pw_env, &
     494           20 :                       virial=virial)
     495              : 
     496           20 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     497              : 
     498              :       IF (use_virial) THEN
     499            0 :          CPABORT("Virial not implemented for image charge method")
     500              :       END IF
     501              : 
     502           20 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     503              : 
     504           20 :       natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
     505              : 
     506           20 :       IF (.NOT. ASSOCIATED(forces)) THEN
     507           30 :          ALLOCATE (forces(3, natom))
     508              :       END IF
     509              : 
     510          180 :       forces(:, :) = 0.0_dp
     511              : 
     512           20 :       CALL reallocate(cores, 1, natom)
     513           20 :       npme = 0
     514           60 :       cores = 0
     515              : 
     516           60 :       DO iatom = 1, natom
     517           60 :          IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
     518              :             ! replicated realspace grid, split the atoms up between procs
     519           40 :             IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
     520           20 :                npme = npme + 1
     521           20 :                cores(npme) = iatom
     522              :             END IF
     523              :          ELSE
     524            0 :             npme = npme + 1
     525            0 :             cores(npme) = iatom
     526              :          END IF
     527              :       END DO
     528              : 
     529           40 :       DO j = 1, npme
     530              : 
     531           20 :          iatom = cores(j)
     532           20 :          atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
     533           20 :          ra(:) = pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell)
     534           20 :          hab(1, 1) = 0.0_dp
     535           20 :          pab(1, 1) = 1.0_dp
     536           20 :          force_a(:) = 0.0_dp
     537           20 :          force_b(:) = 0.0_dp
     538              : 
     539              :          radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
     540              :                                            ra=ra, rb=ra, rp=ra, &
     541              :                                            zetp=qmmm_env%image_charge_pot%eta, eps=eps_rho_rspace, &
     542              :                                            pab=pab, o1=0, o2=0, &  ! without map_consistent
     543           20 :                                            prefactor=1.0_dp, cutoff=1.0_dp)
     544              : 
     545              :          CALL integrate_pgf_product(0, qmmm_env%image_charge_pot%eta, 0, &
     546              :                                     0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
     547              :                                     rs_v, hab, pab, o1=0, o2=0, &
     548              :                                     radius=radius, calculate_forces=.TRUE., &
     549              :                                     force_a=force_a, force_b=force_b, use_subpatch=.TRUE., &
     550           20 :                                     subpatch_pattern=0)
     551              : 
     552           80 :          force_a(:) = coeff(iatom)*force_a(:)
     553          100 :          forces(:, iatom) = force_a(:)
     554              : 
     555              :       END DO
     556              : 
     557          340 :       CALL para_env%sum(forces)
     558              : 
     559           20 :       DEALLOCATE (hab, pab, cores)
     560              : 
     561              :       ! print info on gradients if wanted
     562           20 :       CALL print_gradients_image_atoms(forces, qs_env)
     563              : 
     564           20 :       CALL timestop(handle)
     565              : 
     566           20 :    END SUBROUTINE integrate_potential_devga_rspace
     567              : 
     568              : !****************************************************************************
     569              : !> \brief calculate image matrix T depending on constraints on image atoms
     570              : !>        in case coefficients are estimated not iteratively
     571              : !> \param qs_env qs environment
     572              : !> \param qmmm_env qmmm environment
     573              : ! **************************************************************************************************
     574           20 :    SUBROUTINE conditional_calc_image_matrix(qs_env, qmmm_env)
     575              : 
     576              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     577              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     578              : 
     579           20 :       IF (.NOT. qmmm_env%image_charge_pot%coeff_iterative) THEN
     580           28 :          SELECT CASE (qmmm_env%image_charge_pot%state_image_matrix)
     581              :          CASE (calc_always)
     582              :             CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
     583           12 :                                         ipiv=qs_env%ipiv, qs_env=qs_env, qmmm_env=qmmm_env)
     584              :          CASE (calc_once)
     585              :             !if all image atoms are fully constrained, calculate image matrix
     586              :             !only for the first MD or GEO_OPT step
     587              :             CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
     588            2 :                                         ipiv=qs_env%ipiv, qs_env=qs_env, qmmm_env=qmmm_env)
     589            2 :             qmmm_env%image_charge_pot%state_image_matrix = calc_once_done
     590            2 :             IF (qmmm_env%center_qm_subsys0) &
     591              :                CALL cp_warn(__LOCATION__, &
     592              :                             "The image atoms are fully "// &
     593              :                             "constrained and the image matrix is only calculated once. "// &
     594            0 :                             "To be safe, set CENTER to NEVER ")
     595              :          CASE (calc_once_done)
     596              :             ! do nothing image matrix is stored
     597              :          CASE DEFAULT
     598           16 :             CPABORT("No initialization for image charges available?")
     599              :          END SELECT
     600              :       END IF
     601              : 
     602           20 :    END SUBROUTINE conditional_calc_image_matrix
     603              : 
     604              : !****************************************************************************
     605              : !> \brief calculate image matrix T
     606              : !> \param image_matrix matrix T
     607              : !> \param ipiv pivoting prior to DGETRS (for Gaussian elimination)
     608              : !> \param qs_env qs environment
     609              : !> \param qmmm_env qmmm environment
     610              : ! **************************************************************************************************
     611           16 :    SUBROUTINE calculate_image_matrix(image_matrix, ipiv, qs_env, qmmm_env)
     612              : 
     613              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: image_matrix
     614              :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: ipiv
     615              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     616              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     617              : 
     618              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_image_matrix'
     619              : 
     620              :       INTEGER                                            :: handle, natom, output_unit, stat
     621              :       TYPE(cp_logger_type), POINTER                      :: logger
     622              :       TYPE(section_vals_type), POINTER                   :: input
     623              : 
     624           16 :       CALL timeset(routineN, handle)
     625           16 :       NULLIFY (input, logger)
     626              : 
     627           16 :       logger => cp_get_default_logger()
     628              : 
     629           16 :       natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
     630              : 
     631           16 :       IF (.NOT. ASSOCIATED(image_matrix)) THEN
     632           40 :          ALLOCATE (image_matrix(natom, natom))
     633              :       END IF
     634           16 :       IF (PRESENT(ipiv)) THEN
     635           14 :          IF (.NOT. ASSOCIATED(ipiv)) THEN
     636           24 :             ALLOCATE (ipiv(natom))
     637              :          END IF
     638           42 :          ipiv = 0
     639              :       END IF
     640              : 
     641           16 :       CALL get_qs_env(qs_env, input=input)
     642              :       !print info
     643              :       output_unit = cp_print_key_unit_nr(logger, input, &
     644              :                                          "QMMM%PRINT%PROGRAM_RUN_INFO", &
     645           16 :                                          extension=".qmmmLog")
     646           16 :       IF (qmmm_env%image_charge_pot%coeff_iterative) THEN
     647            2 :          IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A)") &
     648            1 :             "Calculating image matrix"
     649              :       ELSE
     650           14 :          IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T2,A)") &
     651            7 :             "Calculating image matrix"
     652              :       END IF
     653              :       CALL cp_print_key_finished_output(output_unit, logger, input, &
     654           16 :                                         "QMMM%PRINT%PROGRAM_RUN_INFO")
     655              : 
     656              :       ! Calculate image matrix using either GPW or MME method
     657           20 :       SELECT CASE (qmmm_env%image_charge_pot%image_matrix_method)
     658              :       CASE (do_eri_gpw)
     659            4 :          CALL calculate_image_matrix_gpw(image_matrix, qs_env, qmmm_env)
     660              :       CASE (do_eri_mme)
     661           12 :          CALL calculate_image_matrix_mme(image_matrix, qs_env, qmmm_env)
     662              :       CASE DEFAULT
     663           16 :          CPABORT("Unknown method for calculating image matrix")
     664              :       END SELECT
     665              : 
     666           16 :       IF (qmmm_env%image_charge_pot%coeff_iterative) THEN
     667              :          !inversion --> preconditioner matrix for CG
     668            2 :          CALL invmat_symm(qs_env%image_matrix)
     669            2 :          CALL write_image_matrix(qs_env%image_matrix, qs_env)
     670              :       ELSE
     671              :          !pivoting prior to DGETRS (Gaussian elimination)
     672           14 :          IF (PRESENT(ipiv)) THEN
     673           14 :             CALL dgetrf(natom, natom, image_matrix, natom, ipiv, stat)
     674           14 :             CPASSERT(stat == 0)
     675              :          END IF
     676              :       END IF
     677              : 
     678           16 :       CALL timestop(handle)
     679              : 
     680           16 :    END SUBROUTINE calculate_image_matrix
     681              : 
     682              : ! **************************************************************************************************
     683              : !> \brief calculate image matrix T using GPW method
     684              : !> \param image_matrix matrix T
     685              : !> \param qs_env qs environment
     686              : !> \param qmmm_env qmmm environment
     687              : ! **************************************************************************************************
     688            4 :    SUBROUTINE calculate_image_matrix_gpw(image_matrix, qs_env, qmmm_env)
     689              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: image_matrix
     690              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     691              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     692              : 
     693              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_image_matrix_gpw'
     694              : 
     695              :       INTEGER                                            :: handle, iatom, iatom_ref, natom
     696              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: int_res
     697              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     698              :       TYPE(pw_c1d_gs_type)                               :: rho_gb, vb_gspace
     699              :       TYPE(pw_env_type), POINTER                         :: pw_env
     700              :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     701              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     702              :       TYPE(pw_r3d_rs_type)                               :: vb_rspace
     703              : 
     704            4 :       CALL timeset(routineN, handle)
     705            4 :       NULLIFY (pw_env, auxbas_pw_pool, poisson_env, para_env, int_res)
     706              : 
     707            4 :       natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
     708           12 :       ALLOCATE (int_res(natom))
     709              : 
     710           28 :       image_matrix = 0.0_dp
     711              : 
     712            4 :       CALL get_qs_env(qs_env, pw_env=pw_env, para_env=para_env)
     713              : 
     714              :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     715            4 :                       poisson_env=poisson_env)
     716            4 :       CALL auxbas_pw_pool%create_pw(rho_gb)
     717            4 :       CALL auxbas_pw_pool%create_pw(vb_gspace)
     718            4 :       CALL auxbas_pw_pool%create_pw(vb_rspace)
     719              : 
     720              :       ! calculate vb only once for one reference atom
     721            4 :       iatom_ref = 1 !
     722              :       !collocate gaussian of reference MM atom on grid
     723            4 :       CALL pw_zero(rho_gb)
     724            4 :       CALL calculate_rho_single_gaussian(rho_gb, qs_env, iatom_ref)
     725              :       !calculate potential vb like hartree potential
     726            4 :       CALL pw_zero(vb_gspace)
     727            4 :       CALL pw_poisson_solve(poisson_env, rho_gb, vhartree=vb_gspace)
     728            4 :       CALL pw_zero(vb_rspace)
     729            4 :       CALL pw_transfer(vb_gspace, vb_rspace)
     730            4 :       CALL pw_scale(vb_rspace, vb_rspace%pw_grid%dvol)
     731              : 
     732           12 :       DO iatom = 1, natom
     733              :          !calculate integral vb_rspace*ga
     734           24 :          int_res = 0.0_dp
     735              :          CALL integrate_potential_ga_rspace(vb_rspace, qs_env%qmmm_env_qm, &
     736              :                                             qs_env, int_res, atom_num=iatom, &
     737            8 :                                             atom_num_ref=iatom_ref)
     738           40 :          image_matrix(iatom, iatom:natom) = int_res(iatom:natom)
     739           28 :          image_matrix(iatom + 1:natom, iatom) = int_res(iatom + 1:natom)
     740              :       END DO
     741              : 
     742            4 :       CALL vb_gspace%release()
     743            4 :       CALL vb_rspace%release()
     744            4 :       CALL rho_gb%release()
     745              : 
     746            4 :       DEALLOCATE (int_res)
     747              : 
     748            4 :       CALL timestop(handle)
     749            4 :    END SUBROUTINE calculate_image_matrix_gpw
     750              : 
     751              : ! **************************************************************************************************
     752              : !> \brief calculate image matrix T using MME (MiniMax-Ewald) method
     753              : !> \param image_matrix matrix T
     754              : !> \param qs_env qs environment
     755              : !> \param qmmm_env qmmm environment
     756              : ! **************************************************************************************************
     757           12 :    SUBROUTINE calculate_image_matrix_mme(image_matrix, qs_env, qmmm_env)
     758              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: image_matrix
     759              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     760              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     761              : 
     762              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_image_matrix_mme'
     763              : 
     764              :       INTEGER                                            :: atom_a, handle, iatom, natom
     765              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: zeta
     766              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: ra
     767              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     768              : 
     769           12 :       CALL timeset(routineN, handle)
     770           12 :       NULLIFY (para_env)
     771              : 
     772           12 :       natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
     773           60 :       ALLOCATE (zeta(natom), ra(3, natom))
     774              : 
     775           36 :       zeta(:) = qmmm_env%image_charge_pot%eta
     776              : 
     777           36 :       DO iatom = 1, natom
     778           24 :          atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
     779          108 :          ra(:, iatom) = qmmm_env%image_charge_pot%particles_all(atom_a)%r(:)
     780              :       END DO
     781              : 
     782           12 :       CALL get_qs_env(qs_env, para_env=para_env)
     783              : 
     784              :       CALL integrate_s_mme(qmmm_env%image_charge_pot%eri_mme_param, &
     785           12 :                            zeta, zeta, ra, ra, image_matrix, para_env)
     786              : 
     787           12 :       CALL timestop(handle)
     788           24 :    END SUBROUTINE calculate_image_matrix_mme
     789              : 
     790              : ! **************************************************************************************************
     791              : !> \brief high-level integration routine for 2c integrals over s-type functions.
     792              : !>        Parallelization over pairs of functions.
     793              : !> \param param ...
     794              : !> \param zeta ...
     795              : !> \param zetb ...
     796              : !> \param ra ...
     797              : !> \param rb ...
     798              : !> \param hab ...
     799              : !> \param para_env ...
     800              : ! **************************************************************************************************
     801           12 :    SUBROUTINE integrate_s_mme(param, zeta, zetb, ra, rb, hab, para_env)
     802              :       TYPE(cp_eri_mme_param), INTENT(INOUT)              :: param
     803              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zeta, zetb
     804              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: ra, rb
     805              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: hab
     806              :       TYPE(mp_para_env_type), INTENT(IN), POINTER        :: para_env
     807              : 
     808              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'integrate_s_mme'
     809              : 
     810              :       INTEGER                                            :: G_count, handle, ipgf, ipgf_prod, jpgf, &
     811              :                                                             npgf_prod, npgfa, npgfb, R_count
     812              :       INTEGER, DIMENSION(2)                              :: limits
     813              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
     814              : 
     815           12 :       CALL timeset(routineN, handle)
     816           12 :       G_count = 0; R_count = 0
     817              : 
     818           84 :       hab(:, :) = 0.0_dp
     819              : 
     820           12 :       npgfa = SIZE(zeta)
     821           12 :       npgfb = SIZE(zetb)
     822           12 :       npgf_prod = npgfa*npgfb ! total number of integrals
     823              : 
     824           12 :       limits = get_limit(npgf_prod, para_env%num_pe, para_env%mepos)
     825              : 
     826           36 :       DO ipgf_prod = limits(1), limits(2)
     827           24 :          ipgf = (ipgf_prod - 1)/npgfb + 1
     828           24 :          jpgf = MOD(ipgf_prod - 1, npgfb) + 1
     829           96 :          rab(:) = ra(:, ipgf) - rb(:, jpgf)
     830              :          CALL eri_mme_2c_integrate(param%par, 0, 0, 0, 0, zeta(ipgf), &
     831           36 :                                    zetb(jpgf), rab, hab, ipgf - 1, jpgf - 1, G_count=G_count, R_count=R_count)
     832              :       END DO
     833              : 
     834           12 :       CALL cp_eri_mme_update_local_counts(param, para_env, G_count_2c=G_count, R_count_2c=R_count)
     835          156 :       CALL para_env%sum(hab)
     836           12 :       CALL timestop(handle)
     837              : 
     838           12 :    END SUBROUTINE integrate_s_mme
     839              : 
     840              : ! **************************************************************************************************
     841              : !> \brief calculates potential of the metal (image potential) given a set of
     842              : !>        coefficients coeff
     843              : !> \param v_metal_rspace potential generated by rho_metal in real space
     844              : !> \param coeff expansion coefficients of the image charge density, i.e.
     845              : !>        rho_metal=sum_a c_a*g_a
     846              : !> \param rho_hartree_gspace Kohn Sham density in reciprocal space
     847              : !> \param energy structure where energies are stored
     848              : !> \param qs_env qs environment
     849              : ! **************************************************************************************************
     850          180 :    SUBROUTINE calculate_potential_metal(v_metal_rspace, coeff, rho_hartree_gspace, energy, &
     851              :                                         qs_env)
     852              : 
     853              :       TYPE(pw_r3d_rs_type), INTENT(OUT)                  :: v_metal_rspace
     854              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: coeff
     855              :       TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL         :: rho_hartree_gspace
     856              :       TYPE(qs_energy_type), OPTIONAL, POINTER            :: energy
     857              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     858              : 
     859              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_potential_metal'
     860              : 
     861              :       INTEGER                                            :: handle
     862              :       REAL(KIND=dp)                                      :: en_external, en_vmetal_rhohartree, &
     863              :                                                             total_rho_metal
     864              :       TYPE(pw_c1d_gs_type)                               :: rho_metal, v_metal_gspace
     865              :       TYPE(pw_env_type), POINTER                         :: pw_env
     866              :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     867              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     868              : 
     869           90 :       CALL timeset(routineN, handle)
     870              : 
     871           90 :       NULLIFY (pw_env, auxbas_pw_pool, poisson_env)
     872           90 :       en_vmetal_rhohartree = 0.0_dp
     873           90 :       en_external = 0.0_dp
     874              : 
     875           90 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     876              :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     877           90 :                       poisson_env=poisson_env)
     878              : 
     879           90 :       CALL auxbas_pw_pool%create_pw(rho_metal)
     880              : 
     881           90 :       CALL auxbas_pw_pool%create_pw(v_metal_gspace)
     882              : 
     883           90 :       CALL auxbas_pw_pool%create_pw(v_metal_rspace)
     884              : 
     885           90 :       CALL pw_zero(rho_metal)
     886              :       CALL calculate_rho_metal(rho_metal, coeff, total_rho_metal=total_rho_metal, &
     887           90 :                                qs_env=qs_env)
     888              : 
     889           90 :       CALL pw_zero(v_metal_gspace)
     890              :       CALL pw_poisson_solve(poisson_env, rho_metal, &
     891           90 :                             vhartree=v_metal_gspace)
     892              : 
     893           90 :       IF (PRESENT(rho_hartree_gspace)) THEN
     894              :          en_vmetal_rhohartree = 0.5_dp*pw_integral_ab(v_metal_gspace, &
     895           60 :                                                       rho_hartree_gspace)
     896           60 :          en_external = qs_env%qmmm_env_qm%image_charge_pot%V0*total_rho_metal
     897           60 :          energy%image_charge = en_vmetal_rhohartree - 0.5_dp*en_external
     898              :          CALL print_image_energy_terms(en_vmetal_rhohartree, en_external, &
     899           60 :                                        total_rho_metal, qs_env)
     900              :       END IF
     901              : 
     902           90 :       CALL pw_zero(v_metal_rspace)
     903           90 :       CALL pw_transfer(v_metal_gspace, v_metal_rspace)
     904           90 :       CALL pw_scale(v_metal_rspace, v_metal_rspace%pw_grid%dvol)
     905           90 :       CALL v_metal_gspace%release()
     906           90 :       CALL rho_metal%release()
     907              : 
     908           90 :       CALL timestop(handle)
     909              : 
     910           90 :    END SUBROUTINE calculate_potential_metal
     911              : 
     912              : ! ****************************************************************************
     913              : !> \brief Add potential of metal (image charge pot) to Hartree Potential
     914              : !> \param v_hartree Hartree potential (in real space)
     915              : !> \param v_metal potential generated by rho_metal (in real space)
     916              : !> \param qs_env qs environment
     917              : ! **************************************************************************************************
     918           60 :    SUBROUTINE add_image_pot_to_hartree_pot(v_hartree, v_metal, qs_env)
     919              : 
     920              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: v_hartree
     921              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_metal
     922              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     923              : 
     924              :       CHARACTER(len=*), PARAMETER :: routineN = 'add_image_pot_to_hartree_pot'
     925              : 
     926              :       INTEGER                                            :: handle, output_unit
     927              :       TYPE(cp_logger_type), POINTER                      :: logger
     928              :       TYPE(section_vals_type), POINTER                   :: input
     929              : 
     930           60 :       CALL timeset(routineN, handle)
     931              : 
     932           60 :       NULLIFY (input, logger)
     933           60 :       logger => cp_get_default_logger()
     934              : 
     935              :       !add image charge potential
     936           60 :       CALL pw_axpy(v_metal, v_hartree)
     937              : 
     938              :       ! print info
     939              :       CALL get_qs_env(qs_env=qs_env, &
     940           60 :                       input=input)
     941              :       output_unit = cp_print_key_unit_nr(logger, input, &
     942              :                                          "QMMM%PRINT%PROGRAM_RUN_INFO", &
     943           60 :                                          extension=".qmmmLog")
     944           60 :       IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A)") &
     945           30 :          "Adding image charge potential to the Hartree potential."
     946              :       CALL cp_print_key_finished_output(output_unit, logger, input, &
     947           60 :                                         "QMMM%PRINT%PROGRAM_RUN_INFO")
     948              : 
     949           60 :       CALL timestop(handle)
     950              : 
     951           60 :    END SUBROUTINE add_image_pot_to_hartree_pot
     952              : 
     953              : !****************************************************************************
     954              : !> \brief writes image matrix T to file when used as preconditioner for
     955              : !>        calculating image coefficients iteratively
     956              : !> \param image_matrix matrix T
     957              : !> \param qs_env qs environment
     958              : ! **************************************************************************************************
     959            2 :    SUBROUTINE write_image_matrix(image_matrix, qs_env)
     960              : 
     961              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: image_matrix
     962              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     963              : 
     964              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'write_image_matrix'
     965              : 
     966              :       CHARACTER(LEN=default_path_length)                 :: filename
     967              :       INTEGER                                            :: handle, rst_unit
     968              :       TYPE(cp_logger_type), POINTER                      :: logger
     969              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     970              :       TYPE(section_vals_type), POINTER                   :: print_key, qmmm_section
     971              : 
     972            2 :       CALL timeset(routineN, handle)
     973              : 
     974            2 :       NULLIFY (qmmm_section, print_key, logger, para_env)
     975            2 :       logger => cp_get_default_logger()
     976            2 :       rst_unit = -1
     977              : 
     978              :       CALL get_qs_env(qs_env=qs_env, para_env=para_env, &
     979            2 :                       input=qmmm_section)
     980              : 
     981              :       print_key => section_vals_get_subs_vals(qmmm_section, &
     982            2 :                                               "QMMM%PRINT%IMAGE_CHARGE_RESTART")
     983              : 
     984            2 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     985              :                                            qmmm_section, "QMMM%PRINT%IMAGE_CHARGE_RESTART"), &
     986              :                 cp_p_file)) THEN
     987              : 
     988              :          rst_unit = cp_print_key_unit_nr(logger, qmmm_section, &
     989              :                                          "QMMM%PRINT%IMAGE_CHARGE_RESTART", &
     990              :                                          extension=".Image", &
     991              :                                          file_status="REPLACE", &
     992              :                                          file_action="WRITE", &
     993            2 :                                          file_form="UNFORMATTED")
     994              : 
     995            2 :          IF (rst_unit > 0) filename = cp_print_key_generate_filename(logger, &
     996              :                                                                      print_key, extension=".IMAGE", &
     997            1 :                                                                      my_local=.FALSE.)
     998              : 
     999            2 :          IF (rst_unit > 0) THEN
    1000            7 :             WRITE (rst_unit) image_matrix
    1001              :          END IF
    1002              : 
    1003              :          CALL cp_print_key_finished_output(rst_unit, logger, qmmm_section, &
    1004            2 :                                            "QMMM%PRINT%IMAGE_CHARGE_RESTART")
    1005              :       END IF
    1006              : 
    1007            2 :       CALL timestop(handle)
    1008              : 
    1009            2 :    END SUBROUTINE write_image_matrix
    1010              : 
    1011              : !****************************************************************************
    1012              : !> \brief restarts image matrix T when used as preconditioner for calculating
    1013              : !>        image coefficients iteratively
    1014              : !> \param image_matrix matrix T
    1015              : !> \param qs_env qs environment
    1016              : !> \param qmmm_env qmmm environment
    1017              : ! **************************************************************************************************
    1018            0 :    SUBROUTINE restart_image_matrix(image_matrix, qs_env, qmmm_env)
    1019              : 
    1020              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: image_matrix
    1021              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1022              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
    1023              : 
    1024              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'restart_image_matrix'
    1025              : 
    1026              :       CHARACTER(LEN=default_path_length)                 :: image_filename
    1027              :       INTEGER                                            :: handle, natom, output_unit, rst_unit
    1028              :       LOGICAL                                            :: exist
    1029              :       TYPE(cp_logger_type), POINTER                      :: logger
    1030              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1031              :       TYPE(section_vals_type), POINTER                   :: qmmm_section
    1032              : 
    1033            0 :       CALL timeset(routineN, handle)
    1034              : 
    1035            0 :       NULLIFY (qmmm_section, logger, para_env)
    1036            0 :       logger => cp_get_default_logger()
    1037            0 :       exist = .FALSE.
    1038            0 :       rst_unit = -1
    1039              : 
    1040            0 :       natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
    1041              : 
    1042            0 :       IF (.NOT. ASSOCIATED(image_matrix)) THEN
    1043            0 :          ALLOCATE (image_matrix(natom, natom))
    1044              :       END IF
    1045              : 
    1046            0 :       image_matrix = 0.0_dp
    1047              : 
    1048              :       CALL get_qs_env(qs_env=qs_env, para_env=para_env, &
    1049            0 :                       input=qmmm_section)
    1050              : 
    1051              :       CALL section_vals_val_get(qmmm_section, "QMMM%IMAGE_CHARGE%IMAGE_RESTART_FILE_NAME", &
    1052            0 :                                 c_val=image_filename)
    1053              : 
    1054            0 :       INQUIRE (FILE=image_filename, exist=exist)
    1055              : 
    1056            0 :       IF (exist) THEN
    1057            0 :          IF (para_env%is_source()) THEN
    1058              :             CALL open_file(file_name=image_filename, &
    1059              :                            file_status="OLD", &
    1060              :                            file_form="UNFORMATTED", &
    1061              :                            file_action="READ", &
    1062            0 :                            unit_number=rst_unit)
    1063              : 
    1064            0 :             READ (rst_unit) qs_env%image_matrix
    1065              :          END IF
    1066              : 
    1067            0 :          CALL para_env%bcast(qs_env%image_matrix)
    1068              : 
    1069            0 :          IF (para_env%is_source()) CALL close_file(unit_number=rst_unit)
    1070              : 
    1071              :          output_unit = cp_print_key_unit_nr(logger, qmmm_section, &
    1072              :                                             "QMMM%PRINT%PROGRAM_RUN_INFO", &
    1073            0 :                                             extension=".qmmmLog")
    1074            0 :          IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A)") &
    1075            0 :             "Restarted image matrix"
    1076              :       ELSE
    1077            0 :          CPABORT("Restart file for image matrix not found")
    1078              :       END IF
    1079              : 
    1080            0 :       qmmm_env%image_charge_pot%image_restart = .FALSE.
    1081              : 
    1082            0 :       CALL timestop(handle)
    1083              : 
    1084            0 :    END SUBROUTINE restart_image_matrix
    1085              : 
    1086              : ! ****************************************************************************
    1087              : !> \brief Print info on image gradients on image MM atoms
    1088              : !> \param forces structure storing the force contribution of the image charges
    1089              : !>        for the metal (MM) atoms (actually these are only the gradients)
    1090              : !> \param qs_env qs environment
    1091              : ! **************************************************************************************************
    1092           20 :    SUBROUTINE print_gradients_image_atoms(forces, qs_env)
    1093              : 
    1094              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: forces
    1095              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1096              : 
    1097              :       INTEGER                                            :: atom_a, iatom, natom, output_unit
    1098              :       REAL(KIND=dp), DIMENSION(3)                        :: sum_gradients
    1099              :       TYPE(cp_logger_type), POINTER                      :: logger
    1100              :       TYPE(section_vals_type), POINTER                   :: input
    1101              : 
    1102           20 :       NULLIFY (input, logger)
    1103           20 :       logger => cp_get_default_logger()
    1104              : 
    1105           20 :       sum_gradients = 0.0_dp
    1106           20 :       natom = SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
    1107              : 
    1108           60 :       DO iatom = 1, natom
    1109          180 :          sum_gradients(:) = sum_gradients(:) + forces(:, iatom)
    1110              :       END DO
    1111              : 
    1112           20 :       CALL get_qs_env(qs_env=qs_env, input=input)
    1113              : 
    1114              :       output_unit = cp_print_key_unit_nr(logger, input, &
    1115           20 :                                          "QMMM%PRINT%DERIVATIVES", extension=".Log")
    1116           20 :       IF (output_unit > 0) THEN
    1117              :          WRITE (unit=output_unit, fmt="(/1X,A)") &
    1118            0 :             "Image gradients [a.u.] on MM image charge atoms after QMMM calculation: "
    1119              :          WRITE (unit=output_unit, fmt="(T4,A4,T27,A1,T50,A1,T74,A1)") &
    1120            0 :             "Atom", "X", "Y", "Z"
    1121            0 :          DO iatom = 1, natom
    1122            0 :             atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
    1123              :             WRITE (unit=output_unit, fmt="(T2,I6,T22,ES12.5,T45,ES12.5,T69,ES12.5)") &
    1124            0 :                atom_a, forces(:, iatom)
    1125              :          END DO
    1126              : 
    1127            0 :          WRITE (unit=output_unit, fmt="(T2,A)") REPEAT("-", 79)
    1128              :          WRITE (unit=output_unit, fmt="(T2,A,T22,ES12.5,T45,ES12.5,T69,ES12.5)") &
    1129            0 :             "sum gradients:", sum_gradients
    1130            0 :          WRITE (unit=output_unit, fmt="(/)")
    1131              :       END IF
    1132              : 
    1133              :       CALL cp_print_key_finished_output(output_unit, logger, input, &
    1134           20 :                                         "QMMM%PRINT%DERIVATIVES")
    1135              : 
    1136           20 :    END SUBROUTINE print_gradients_image_atoms
    1137              : 
    1138              : ! ****************************************************************************
    1139              : !> \brief Print image coefficients
    1140              : !> \param image_coeff expansion coefficients of the image charge density
    1141              : !> \param qs_env qs environment
    1142              : ! **************************************************************************************************
    1143           10 :    SUBROUTINE print_image_coefficients(image_coeff, qs_env)
    1144              : 
    1145              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: image_coeff
    1146              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1147              : 
    1148              :       INTEGER                                            :: atom_a, iatom, natom, output_unit
    1149              :       REAL(KIND=dp)                                      :: normalize_factor, sum_coeff
    1150              :       TYPE(cp_logger_type), POINTER                      :: logger
    1151              :       TYPE(section_vals_type), POINTER                   :: input
    1152              : 
    1153           10 :       NULLIFY (input, logger)
    1154           10 :       logger => cp_get_default_logger()
    1155              : 
    1156           10 :       sum_coeff = 0.0_dp
    1157           10 :       natom = SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
    1158           10 :       normalize_factor = SQRT((qs_env%qmmm_env_qm%image_charge_pot%eta/pi)**3)
    1159              : 
    1160           30 :       DO iatom = 1, natom
    1161           30 :          sum_coeff = sum_coeff + image_coeff(iatom)
    1162              :       END DO
    1163              : 
    1164           10 :       CALL get_qs_env(qs_env=qs_env, input=input)
    1165              : 
    1166              :       output_unit = cp_print_key_unit_nr(logger, input, &
    1167           10 :                                          "QMMM%PRINT%IMAGE_CHARGE_INFO", extension=".Log")
    1168           10 :       IF (output_unit > 0) THEN
    1169            2 :          WRITE (unit=output_unit, fmt="(/)")
    1170              :          WRITE (unit=output_unit, fmt="(T2,A)") &
    1171            2 :             "Image charges [a.u.] after QMMM calculation: "
    1172            2 :          WRITE (unit=output_unit, fmt="(T4,A4,T67,A)") "Atom", "Image charge"
    1173            2 :          WRITE (unit=output_unit, fmt="(T4,A,T67,A)") REPEAT("-", 4), REPEAT("-", 12)
    1174              : 
    1175            6 :          DO iatom = 1, natom
    1176            4 :             atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
    1177              :             !opposite sign for image_coeff; during the calculation they have
    1178              :             !the 'wrong' sign to ensure consistency with v_hartree which has
    1179              :             !the opposite sign
    1180              :             WRITE (unit=output_unit, fmt="(T2,I6,T65,ES16.9)") &
    1181            6 :                atom_a, -image_coeff(iatom)/normalize_factor
    1182              :          END DO
    1183              : 
    1184            2 :          WRITE (unit=output_unit, fmt="(T2,A)") REPEAT("-", 79)
    1185              :          WRITE (unit=output_unit, fmt="(T2,A,T65,ES16.9)") &
    1186            2 :             "sum image charges:", -sum_coeff/normalize_factor
    1187              :       END IF
    1188              : 
    1189              :       CALL cp_print_key_finished_output(output_unit, logger, input, &
    1190           10 :                                         "QMMM%PRINT%IMAGE_CHARGE_INFO")
    1191              : 
    1192           10 :    END SUBROUTINE print_image_coefficients
    1193              : 
    1194              : ! ****************************************************************************
    1195              : !> \brief Print detailed image charge energies
    1196              : !> \param en_vmetal_rhohartree energy contribution of the image charges
    1197              : !>        without external potential, i.e. 0.5*integral(v_metal*rho_hartree)
    1198              : !> \param en_external additional energy contribution of the image charges due
    1199              : !>        to an external potential, i.e. V0*total_rho_metal
    1200              : !> \param total_rho_metal total induced image charge density
    1201              : !> \param qs_env qs environment
    1202              : ! **************************************************************************************************
    1203           60 :    SUBROUTINE print_image_energy_terms(en_vmetal_rhohartree, en_external, &
    1204              :                                        total_rho_metal, qs_env)
    1205              : 
    1206              :       REAL(KIND=dp), INTENT(IN)                          :: en_vmetal_rhohartree, en_external, &
    1207              :                                                             total_rho_metal
    1208              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1209              : 
    1210              :       INTEGER                                            :: output_unit
    1211              :       TYPE(cp_logger_type), POINTER                      :: logger
    1212              :       TYPE(section_vals_type), POINTER                   :: input
    1213              : 
    1214           60 :       NULLIFY (input, logger)
    1215           60 :       logger => cp_get_default_logger()
    1216              : 
    1217           60 :       CALL get_qs_env(qs_env=qs_env, input=input)
    1218              : 
    1219              :       output_unit = cp_print_key_unit_nr(logger, input, &
    1220           60 :                                          "QMMM%PRINT%IMAGE_CHARGE_INFO", extension=".Log")
    1221              : 
    1222           60 :       IF (output_unit > 0) THEN
    1223              :          WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
    1224            6 :             "Total induced charge density [a.u.]:", total_rho_metal
    1225            6 :          WRITE (unit=output_unit, fmt="(T3,A)") "Image energy terms: "
    1226              :          WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
    1227            6 :             "Coulomb energy of QM and image charge density [a.u.]:", en_vmetal_rhohartree
    1228              :          WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
    1229            6 :             "External potential energy term [a.u.]:", -0.5_dp*en_external
    1230              :          WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
    1231            6 :             "Total image charge energy [a.u.]:", en_vmetal_rhohartree - 0.5_dp*en_external
    1232              :       END IF
    1233              : 
    1234              :       CALL cp_print_key_finished_output(output_unit, logger, input, &
    1235           60 :                                         "QMMM%PRINT%IMAGE_CHARGE_INFO")
    1236              : 
    1237           60 :    END SUBROUTINE print_image_energy_terms
    1238              : 
    1239           30 : END MODULE qmmm_image_charge
        

Generated by: LCOV version 2.0-1