LCOV - code coverage report
Current view: top level - src - qmmm_image_charge.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 382 424 90.1 %
Date: 2024-03-29 07:50:05 Functions: 16 17 94.1 %

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

Generated by: LCOV version 1.15