LCOV - code coverage report
Current view: top level - src - cp_ddapc_types.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.3 % 164 158
Test Date: 2025-07-25 12:55:17 Functions: 71.4 % 7 5

            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 contains information regarding the decoupling/recoupling method of Bloechl
      10              : !> \author Teodoro Laino
      11              : ! **************************************************************************************************
      12              : MODULE cp_ddapc_types
      13              :    USE cell_methods,                    ONLY: read_cell
      14              :    USE cell_types,                      ONLY: cell_release,&
      15              :                                               cell_type
      16              :    USE cp_ddapc_methods,                ONLY: ddapc_eval_AmI,&
      17              :                                               ddapc_eval_gfunc,&
      18              :                                               ewald_ddapc_pot,&
      19              :                                               solvation_ddapc_pot
      20              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      21              :                                               cp_logger_get_default_io_unit,&
      22              :                                               cp_logger_type
      23              :    USE cp_output_handling,              ONLY: cp_printkey_is_on
      24              :    USE ewald_spline_util,               ONLY: Setup_Ewald_Spline
      25              :    USE input_section_types,             ONLY: section_vals_get,&
      26              :                                               section_vals_get_subs_vals,&
      27              :                                               section_vals_type,&
      28              :                                               section_vals_val_get
      29              :    USE kinds,                           ONLY: dp
      30              :    USE mathconstants,                   ONLY: pi
      31              :    USE message_passing,                 ONLY: mp_para_env_type
      32              :    USE particle_types,                  ONLY: particle_type
      33              :    USE pw_grid_types,                   ONLY: pw_grid_type
      34              :    USE pw_grids,                        ONLY: pw_grid_release
      35              :    USE pw_poisson_types,                ONLY: pw_poisson_multipole
      36              :    USE pw_pool_types,                   ONLY: pw_pool_release,&
      37              :                                               pw_pool_type
      38              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      39              :                                               pw_r3d_rs_type
      40              : #include "./base/base_uses.f90"
      41              : 
      42              :    IMPLICIT NONE
      43              :    PRIVATE
      44              :    LOGICAL, PRIVATE, PARAMETER          :: debug_this_module = .TRUE.
      45              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc_types'
      46              : 
      47              :    PUBLIC :: cp_ddapc_type, cp_ddapc_create, cp_ddapc_release
      48              :    PUBLIC :: cp_ddapc_ewald_type, cp_ddapc_ewald_create, cp_ddapc_ewald_release
      49              : 
      50              : ! **************************************************************************************************
      51              : !> \author Teodoro Laino
      52              : ! **************************************************************************************************
      53              :    TYPE cp_ddapc_type
      54              :       REAL(KIND=dp) :: c0 = 0.0_dp
      55              :       REAL(KIND=dp), DIMENSION(:, :), POINTER     :: AmI => NULL()
      56              :       REAL(KIND=dp), DIMENSION(:, :), POINTER     :: Md => NULL() ! decoupling
      57              :       REAL(KIND=dp), DIMENSION(:, :), POINTER     :: Mr => NULL() ! recoupling
      58              :       REAL(KIND=dp), DIMENSION(:, :), POINTER     :: Mt => NULL() ! decoupling+recoupling
      59              :       REAL(KIND=dp), DIMENSION(:, :), POINTER     :: Ms => NULL() ! solvation
      60              :       REAL(KIND=dp), POINTER, DIMENSION(:, :)  :: gfunc => NULL()
      61              :       REAL(KIND=dp), POINTER, DIMENSION(:)    :: w => NULL()
      62              :    END TYPE cp_ddapc_type
      63              : 
      64              : ! **************************************************************************************************
      65              :    TYPE cp_ddapc_ewald_type
      66              :       LOGICAL                                    :: do_decoupling = .FALSE.
      67              :       LOGICAL                                    :: do_qmmm_periodic_decpl = .FALSE.
      68              :       LOGICAL                                    :: do_solvation = .FALSE.
      69              :       LOGICAL                                    :: do_property = .FALSE.
      70              :       LOGICAL                                    :: do_restraint = .FALSE.
      71              :       TYPE(section_vals_type), POINTER :: ewald_section => NULL()
      72              :       TYPE(pw_pool_type), POINTER :: pw_pool_qm => NULL(), pw_pool_mm => NULL()
      73              :       TYPE(pw_grid_type), POINTER :: pw_grid_qm => NULL(), pw_grid_mm => NULL()
      74              :       TYPE(pw_r3d_rs_type), POINTER :: coeff_qm => NULL(), coeff_mm => NULL()
      75              :    END TYPE cp_ddapc_ewald_type
      76              : 
      77              : CONTAINS
      78              : 
      79              : ! **************************************************************************************************
      80              : !> \brief ...
      81              : !> \param cp_para_env ...
      82              : !> \param cp_ddapc_env ...
      83              : !> \param cp_ddapc_ewald ...
      84              : !> \param particle_set ...
      85              : !> \param radii ...
      86              : !> \param cell ...
      87              : !> \param super_cell ...
      88              : !> \param rho_tot_g ...
      89              : !> \param gcut ...
      90              : !> \param iw2 ...
      91              : !> \param Vol ...
      92              : !> \param force_env_section ...
      93              : !> \author Tedoro Laino
      94              : !> \note NB receive cp_para_env to pass down to parallelized ewald_ddapc_pot()
      95              : ! **************************************************************************************************
      96          248 :    SUBROUTINE cp_ddapc_create(cp_para_env, cp_ddapc_env, cp_ddapc_ewald, &
      97              :                               particle_set, radii, cell, super_cell, rho_tot_g, gcut, iw2, Vol, &
      98              :                               force_env_section)
      99              :       TYPE(mp_para_env_type), POINTER                    :: cp_para_env
     100              :       TYPE(cp_ddapc_type), INTENT(OUT)                   :: cp_ddapc_env
     101              :       TYPE(cp_ddapc_ewald_type), POINTER                 :: cp_ddapc_ewald
     102              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     103              :       REAL(kind=dp), DIMENSION(:), POINTER               :: radii
     104              :       TYPE(cell_type), POINTER                           :: cell, super_cell
     105              :       TYPE(pw_c1d_gs_type), INTENT(IN)                   :: rho_tot_g
     106              :       REAL(KIND=dp), INTENT(IN)                          :: gcut
     107              :       INTEGER, INTENT(IN)                                :: iw2
     108              :       REAL(KIND=dp), INTENT(IN)                          :: Vol
     109              :       TYPE(section_vals_type), POINTER                   :: force_env_section
     110              : 
     111              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cp_ddapc_create'
     112              : 
     113              :       INTEGER                                            :: handle
     114              :       TYPE(section_vals_type), POINTER                   :: param_section, solvation_section
     115              : 
     116          248 :       CALL timeset(routineN, handle)
     117              :       NULLIFY (cp_ddapc_env%AmI, &
     118          248 :                cp_ddapc_env%Md, &
     119          248 :                cp_ddapc_env%Mt, &
     120          248 :                cp_ddapc_env%Mr, &
     121          248 :                cp_ddapc_env%Ms, &
     122          248 :                cp_ddapc_env%gfunc, &
     123          248 :                cp_ddapc_env%w)
     124              :       ! Evaluates gfunc and AmI
     125          248 :       CALL ddapc_eval_gfunc(cp_ddapc_env%gfunc, cp_ddapc_env%w, gcut, rho_tot_g, radii)
     126              :       CALL ddapc_eval_AmI(cp_ddapc_env%AmI, &
     127              :                           cp_ddapc_env%c0, &
     128              :                           cp_ddapc_env%gfunc, &
     129              :                           cp_ddapc_env%w, &
     130              :                           particle_set, &
     131              :                           gcut, &
     132              :                           rho_tot_g, &
     133              :                           radii, &
     134              :                           iw2, &
     135          248 :                           Vol)
     136          248 :       IF (cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. &
     137              :           cp_ddapc_ewald%do_decoupling) THEN
     138              :          !
     139              :          ! Evaluate the matrix for the Classical contribution to the coupling/decoupling scheme
     140              :          !
     141          114 :          param_section => cp_ddapc_ewald%ewald_section
     142              :          !NB parallelized ewald_ddapc_pot() needs cp_para_env
     143              :          CALL ewald_ddapc_pot(cp_para_env, cp_ddapc_ewald%coeff_qm, &
     144              :                               1.0_dp, &
     145              :                               cell, &
     146              :                               param_section, &
     147              :                               particle_set, &
     148              :                               cp_ddapc_env%Md, &
     149          114 :                               radii)
     150          114 :          IF (cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. cp_ddapc_ewald%do_decoupling) THEN
     151          456 :             ALLOCATE (cp_ddapc_env%Mt(SIZE(cp_ddapc_env%Md, 1), SIZE(cp_ddapc_env%Md, 2)))
     152          114 :             IF (cp_ddapc_ewald%do_decoupling) THEN
     153              :                ! Just decoupling
     154         4900 :                cp_ddapc_env%Mt = cp_ddapc_env%Md
     155              :             ELSE
     156              :                ! QMMM periodic calculation
     157              :                !NB parallelized ewald_ddapc_pot() needs cp_para_env
     158              :                CALL ewald_ddapc_pot(cp_para_env, cp_ddapc_ewald%coeff_mm, -1.0_dp, super_cell, param_section, &
     159           62 :                                     particle_set, cp_ddapc_env%Mr, radii)
     160        41762 :                cp_ddapc_env%Mt = cp_ddapc_env%Md + cp_ddapc_env%Mr
     161              :             END IF
     162              :          END IF
     163              :       END IF
     164          248 :       IF (cp_ddapc_ewald%do_solvation) THEN
     165              :          ! Spherical Solvation model
     166           26 :          solvation_section => section_vals_get_subs_vals(force_env_section, "DFT%SCRF")
     167              :          CALL solvation_ddapc_pot(solvation_section, &
     168           26 :                                   particle_set, cp_ddapc_env%Ms, radii)
     169              :       END IF
     170          248 :       CALL timestop(handle)
     171          248 :    END SUBROUTINE cp_ddapc_create
     172              : 
     173              : ! **************************************************************************************************
     174              : !> \brief ...
     175              : !> \param cp_ddapc_env ...
     176              : !> \par History
     177              : !>      none
     178              : !> \author Teodoro Laino - [tlaino]
     179              : ! **************************************************************************************************
     180          248 :    SUBROUTINE cp_ddapc_release(cp_ddapc_env)
     181              :       TYPE(cp_ddapc_type), INTENT(INOUT)                 :: cp_ddapc_env
     182              : 
     183          248 :       IF (ASSOCIATED(cp_ddapc_env%AmI)) THEN
     184          248 :          DEALLOCATE (cp_ddapc_env%AmI)
     185              :       END IF
     186          248 :       IF (ASSOCIATED(cp_ddapc_env%Mt)) THEN
     187          114 :          DEALLOCATE (cp_ddapc_env%Mt)
     188              :       END IF
     189          248 :       IF (ASSOCIATED(cp_ddapc_env%Md)) THEN
     190          114 :          DEALLOCATE (cp_ddapc_env%Md)
     191              :       END IF
     192          248 :       IF (ASSOCIATED(cp_ddapc_env%Mr)) THEN
     193           62 :          DEALLOCATE (cp_ddapc_env%Mr)
     194              :       END IF
     195          248 :       IF (ASSOCIATED(cp_ddapc_env%Ms)) THEN
     196           26 :          DEALLOCATE (cp_ddapc_env%Ms)
     197              :       END IF
     198          248 :       IF (ASSOCIATED(cp_ddapc_env%gfunc)) THEN
     199          248 :          DEALLOCATE (cp_ddapc_env%gfunc)
     200              :       END IF
     201          248 :       IF (ASSOCIATED(cp_ddapc_env%w)) THEN
     202          248 :          DEALLOCATE (cp_ddapc_env%w)
     203              :       END IF
     204              : 
     205          248 :    END SUBROUTINE cp_ddapc_release
     206              : 
     207              : ! **************************************************************************************************
     208              : !> \brief ...
     209              : !> \param cp_ddapc_ewald ...
     210              : !> \param qmmm_decoupl ...
     211              : !> \param qm_cell ...
     212              : !> \param force_env_section ...
     213              : !> \param subsys_section ...
     214              : !> \param para_env ...
     215              : !> \par History
     216              : !>      none
     217              : !> \author Teodoro Laino - [tlaino]
     218              : ! **************************************************************************************************
     219        37020 :    SUBROUTINE cp_ddapc_ewald_create(cp_ddapc_ewald, qmmm_decoupl, qm_cell, &
     220              :                                     force_env_section, subsys_section, para_env)
     221              :       TYPE(cp_ddapc_ewald_type), POINTER                 :: cp_ddapc_ewald
     222              :       LOGICAL, INTENT(IN)                                :: qmmm_decoupl
     223              :       TYPE(cell_type), POINTER                           :: qm_cell
     224              :       TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
     225              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     226              : 
     227              :       INTEGER                                            :: my_val, npts(3), output_unit
     228         7404 :       INTEGER, DIMENSION(:), POINTER                     :: ngrids
     229              :       LOGICAL                                            :: analyt, decoupling, &
     230              :                                                             do_qmmm_periodic_decpl, do_restraint, &
     231              :                                                             do_restraintB, do_solvation
     232              :       REAL(KIND=dp)                                      :: hmat(3, 3)
     233         7404 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: gx, gy, gz, LG
     234              :       TYPE(cell_type), POINTER                           :: dummy_cell, mm_cell
     235              :       TYPE(cp_logger_type), POINTER                      :: logger
     236              :       TYPE(section_vals_type), POINTER :: cell_section, grid_print_section, multipole_section, &
     237              :          poisson_section, printC_section, qmmm_per_section, restraint_section, restraint_sectionB, &
     238              :          solvation_section
     239              : 
     240        14808 :       logger => cp_get_default_logger()
     241         7404 :       output_unit = cp_logger_get_default_io_unit(logger)
     242         7404 :       CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald))
     243         7404 :       ALLOCATE (cp_ddapc_ewald)
     244              :       NULLIFY (cp_ddapc_ewald%pw_grid_mm, &
     245              :                cp_ddapc_ewald%pw_grid_qm, &
     246              :                cp_ddapc_ewald%ewald_section, &
     247              :                cp_ddapc_ewald%pw_pool_mm, &
     248              :                cp_ddapc_ewald%pw_pool_qm, &
     249              :                cp_ddapc_ewald%coeff_mm, &
     250              :                cp_ddapc_ewald%coeff_qm)
     251         7404 :       NULLIFY (multipole_section)
     252              : 
     253         7404 :       poisson_section => section_vals_get_subs_vals(force_env_section, "DFT%POISSON")
     254         7404 :       solvation_section => section_vals_get_subs_vals(force_env_section, "DFT%SCRF")
     255         7404 :       qmmm_per_section => section_vals_get_subs_vals(force_env_section, "QMMM%PERIODIC")
     256         7404 :       printC_section => section_vals_get_subs_vals(force_env_section, "PROPERTIES%FIT_CHARGE")
     257         7404 :       restraint_section => section_vals_get_subs_vals(force_env_section, "DFT%QS%DDAPC_RESTRAINT")
     258              :       restraint_sectionB => section_vals_get_subs_vals(force_env_section, &
     259         7404 :                                                        "PROPERTIES%ET_COUPLING%DDAPC_RESTRAINT_A")
     260         7404 :       CALL section_vals_get(solvation_section, explicit=do_solvation)
     261         7404 :       CALL section_vals_get(poisson_section, explicit=decoupling)
     262         7404 :       CALL section_vals_get(restraint_section, explicit=do_restraint)
     263         7404 :       CALL section_vals_get(restraint_sectionB, explicit=do_restraintB)
     264         7404 :       do_qmmm_periodic_decpl = qmmm_decoupl
     265         7404 :       cp_ddapc_ewald%do_solvation = do_solvation
     266         7404 :       cp_ddapc_ewald%do_qmmm_periodic_decpl = do_qmmm_periodic_decpl
     267         7404 :       cp_ddapc_ewald%do_property = cp_printkey_is_on(logger%iter_info, printC_section)
     268         7404 :       cp_ddapc_ewald%do_restraint = do_restraint .OR. do_restraintB
     269              :       ! Determining the tasks and further check
     270         7404 :       IF (do_qmmm_periodic_decpl .AND. decoupling) THEN
     271              :          ! Check than an additional POISSON section has not been defined. In case write a warning
     272            0 :          IF (output_unit > 0) &
     273              :             WRITE (output_unit, '(T2,"WARNING",A)') &
     274            0 :             "A calculation with the QMMM periodic model has been requested.", &
     275            0 :             "The explicit POISSON section in DFT section will be IGNORED.", &
     276            0 :             "QM Electrostatic controlled only by the PERIODIC section in QMMM section"
     277            0 :          decoupling = .FALSE.
     278              :       END IF
     279         7404 :       IF (decoupling) THEN
     280              :          ! Simple decoupling technique
     281         2038 :          CALL section_vals_val_get(poisson_section, "POISSON_SOLVER", i_val=my_val)
     282           20 :          SELECT CASE (my_val)
     283              :          CASE (pw_poisson_multipole)
     284           20 :             multipole_section => section_vals_get_subs_vals(poisson_section, "MULTIPOLE")
     285              :          CASE DEFAULT
     286         2038 :             decoupling = .FALSE.
     287              :          END SELECT
     288              :       END IF
     289         7404 :       cp_ddapc_ewald%do_decoupling = decoupling
     290         7404 :       IF (cp_ddapc_ewald%do_qmmm_periodic_decpl) THEN
     291              :          ! QMMM periodic
     292           38 :          multipole_section => section_vals_get_subs_vals(qmmm_per_section, "MULTIPOLE")
     293              :       END IF
     294         7404 :       cp_ddapc_ewald%ewald_section => multipole_section
     295         7404 :       IF (cp_ddapc_ewald%do_decoupling .OR. cp_ddapc_ewald%do_qmmm_periodic_decpl) THEN
     296              :          ! Do we do the calculation analytically or interpolating the g-space factor?
     297           58 :          CALL section_vals_val_get(multipole_section, "ANALYTICAL_GTERM", l_val=analyt)
     298           58 :          IF (.NOT. analyt) THEN
     299           22 :             CALL section_vals_val_get(multipole_section, "ngrids", i_vals=ngrids)
     300           88 :             npts = ngrids
     301              : 
     302           22 :             NULLIFY (LG, gx, gy, gz)
     303          286 :             hmat = qm_cell%hmat
     304           22 :             CALL eval_lg(multipole_section, hmat, qm_cell%deth, LG, gx, gy, gz)
     305           22 :             grid_print_section => section_vals_get_subs_vals(force_env_section, "PRINT%GRID_INFORMATION")
     306              :             CALL Setup_Ewald_Spline(pw_grid=cp_ddapc_ewald%pw_grid_qm, pw_pool=cp_ddapc_ewald%pw_pool_qm, &
     307              :                                     coeff=cp_ddapc_ewald%coeff_qm, LG=LG, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, &
     308              :                                     param_section=multipole_section, tag="ddapc", &
     309           22 :                                     print_section=grid_print_section)
     310           22 :             DEALLOCATE (LG)
     311           22 :             DEALLOCATE (gx)
     312           22 :             DEALLOCATE (gy)
     313           22 :             DEALLOCATE (gz)
     314           22 :             IF (cp_ddapc_ewald%do_qmmm_periodic_decpl) THEN
     315           16 :                NULLIFY (mm_cell, dummy_cell)
     316           16 :                cell_section => section_vals_get_subs_vals(subsys_section, "CELL")
     317           16 :                CALL read_cell(mm_cell, dummy_cell, cell_section=cell_section, para_env=para_env)
     318          208 :                hmat = mm_cell%hmat
     319           16 :                CALL eval_lg(multipole_section, hmat, mm_cell%deth, LG, gx, gy, gz)
     320           16 :                grid_print_section => section_vals_get_subs_vals(force_env_section, "PRINT%GRID_INFORMATION")
     321              :                CALL Setup_Ewald_Spline(pw_grid=cp_ddapc_ewald%pw_grid_mm, pw_pool=cp_ddapc_ewald%pw_pool_mm, &
     322              :                                        coeff=cp_ddapc_ewald%coeff_mm, LG=LG, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, &
     323              :                                        param_section=multipole_section, tag="ddapc", &
     324           16 :                                        print_section=grid_print_section)
     325           16 :                DEALLOCATE (LG)
     326           16 :                DEALLOCATE (gx)
     327           16 :                DEALLOCATE (gy)
     328           16 :                DEALLOCATE (gz)
     329           16 :                CALL cell_release(dummy_cell)
     330           16 :                CALL cell_release(mm_cell)
     331              :             END IF
     332              :          END IF
     333              :       END IF
     334         7404 :    END SUBROUTINE cp_ddapc_ewald_create
     335              : 
     336              : ! **************************************************************************************************
     337              : !> \brief ...
     338              : !> \param multipole_section ...
     339              : !> \param hmat ...
     340              : !> \param deth ...
     341              : !> \param LG ...
     342              : !> \param gx ...
     343              : !> \param gy ...
     344              : !> \param gz ...
     345              : !> \par History
     346              : !>      none
     347              : !> \author Teodoro Laino - [tlaino]
     348              : ! **************************************************************************************************
     349           38 :    SUBROUTINE eval_lg(multipole_section, hmat, deth, LG, gx, gy, gz)
     350              :       TYPE(section_vals_type), POINTER                   :: multipole_section
     351              :       REAL(KIND=dp), INTENT(IN)                          :: hmat(3, 3), deth
     352              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: LG, gx, gy, gz
     353              : 
     354              :       INTEGER                                            :: i, k1, k2, k3, n_rep, ndim, nmax1, &
     355              :                                                             nmax2, nmax3
     356              :       REAL(KIND=dp)                                      :: alpha, eps, fac, fs, fvec(3), galpha, &
     357              :                                                             gsq, gsqi, rcut, tol, tol1
     358              : 
     359           38 :       rcut = MIN(hmat(1, 1), hmat(2, 2), hmat(3, 3))/2.0_dp
     360           38 :       CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep)
     361           38 :       IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut)
     362           38 :       CALL section_vals_val_get(multipole_section, "EWALD_PRECISION", r_val=eps)
     363           38 :       eps = MIN(ABS(eps), 0.5_dp)
     364           38 :       tol = SQRT(ABS(LOG(eps*rcut)))
     365           38 :       alpha = SQRT(ABS(LOG(eps*rcut*tol)))/rcut
     366           38 :       galpha = 1.0_dp/(4.0_dp*alpha*alpha)
     367           38 :       tol1 = SQRT(-LOG(eps*rcut*(2.0_dp*tol*alpha)**2))
     368           38 :       nmax1 = NINT(0.25_dp + hmat(1, 1)*alpha*tol1/pi)
     369           38 :       nmax2 = NINT(0.25_dp + hmat(2, 2)*alpha*tol1/pi)
     370           38 :       nmax3 = NINT(0.25_dp + hmat(3, 3)*alpha*tol1/pi)
     371           38 :       fac = 1.e0_dp/deth
     372          152 :       fvec = 2.0_dp*pi/(/hmat(1, 1), hmat(2, 2), hmat(3, 3)/)
     373           38 :       ndim = (nmax1 + 1)*(2*nmax2 + 1)*(2*nmax3 + 1) - 1
     374          114 :       ALLOCATE (LG(ndim))
     375           76 :       ALLOCATE (gx(ndim))
     376           76 :       ALLOCATE (gy(ndim))
     377           76 :       ALLOCATE (gz(ndim))
     378              : 
     379           38 :       i = 0
     380          216 :       DO k1 = 0, nmax1
     381         2234 :          DO k2 = -nmax2, nmax2
     382        30350 :             DO k3 = -nmax3, nmax3
     383        28154 :                IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) CYCLE
     384        28116 :                i = i + 1
     385        28116 :                fs = 2.0_dp; IF (k1 == 0) fs = 1.0_dp
     386        28116 :                gx(i) = fvec(1)*REAL(k1, KIND=dp)
     387        28116 :                gy(i) = fvec(2)*REAL(k2, KIND=dp)
     388        28116 :                gz(i) = fvec(3)*REAL(k3, KIND=dp)
     389        28116 :                gsq = gx(i)*gx(i) + gy(i)*gy(i) + gz(i)*gz(i)
     390        28116 :                gsqi = fs/gsq
     391        30172 :                LG(i) = fac*gsqi*EXP(-galpha*gsq)
     392              :             END DO
     393              :          END DO
     394              :       END DO
     395              : 
     396           76 :    END SUBROUTINE eval_lg
     397              : 
     398              : ! **************************************************************************************************
     399              : !> \brief ...
     400              : !> \param cp_ddapc_ewald ...
     401              : !> \par History
     402              : !>      none
     403              : !> \author Teodoro Laino - [tlaino]
     404              : ! **************************************************************************************************
     405         7413 :    SUBROUTINE cp_ddapc_ewald_release(cp_ddapc_ewald)
     406              :       TYPE(cp_ddapc_ewald_type), POINTER                 :: cp_ddapc_ewald
     407              : 
     408         7413 :       IF (ASSOCIATED(cp_ddapc_ewald)) THEN
     409         7404 :          IF (ASSOCIATED(cp_ddapc_ewald%coeff_qm)) THEN
     410           22 :             CALL cp_ddapc_ewald%pw_pool_qm%give_back_pw(cp_ddapc_ewald%coeff_qm)
     411           22 :             DEALLOCATE (cp_ddapc_ewald%coeff_qm)
     412              :          END IF
     413         7404 :          IF (ASSOCIATED(cp_ddapc_ewald%coeff_mm)) THEN
     414           16 :             CALL cp_ddapc_ewald%pw_pool_mm%give_back_pw(cp_ddapc_ewald%coeff_mm)
     415           16 :             DEALLOCATE (cp_ddapc_ewald%coeff_mm)
     416              :          END IF
     417         7404 :          IF (ASSOCIATED(cp_ddapc_ewald%pw_pool_qm)) THEN
     418           22 :             CALL pw_pool_release(cp_ddapc_ewald%pw_pool_qm)
     419           22 :             CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald%pw_pool_qm))
     420              :          END IF
     421         7404 :          IF (ASSOCIATED(cp_ddapc_ewald%pw_pool_mm)) THEN
     422           16 :             CALL pw_pool_release(cp_ddapc_ewald%pw_pool_mm)
     423           16 :             CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald%pw_pool_mm))
     424              :          END IF
     425         7404 :          IF (ASSOCIATED(cp_ddapc_ewald%pw_grid_qm)) THEN
     426           22 :             CALL pw_grid_release(cp_ddapc_ewald%pw_grid_qm)
     427           22 :             CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald%pw_grid_qm))
     428              :          END IF
     429         7404 :          IF (ASSOCIATED(cp_ddapc_ewald%pw_grid_mm)) THEN
     430           16 :             CALL pw_grid_release(cp_ddapc_ewald%pw_grid_mm)
     431           16 :             CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald%pw_grid_mm))
     432              :          END IF
     433         7404 :          DEALLOCATE (cp_ddapc_ewald)
     434              :       END IF
     435              : 
     436         7413 :    END SUBROUTINE cp_ddapc_ewald_release
     437              : 
     438            0 : END MODULE cp_ddapc_types
        

Generated by: LCOV version 2.0-1