LCOV - code coverage report
Current view: top level - src - qs_integrate_potential_product.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.6 % 148 143
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            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 Build up the plane wave density by collocating the primitive Gaussian
      10              : !>      functions (pgf).
      11              : !> \par History
      12              : !>      Joost VandeVondele (02.2002)
      13              : !>            1) rewrote collocate_pgf for increased accuracy and speed
      14              : !>            2) collocate_core hack for PGI compiler
      15              : !>            3) added multiple grid feature
      16              : !>            4) new way to go over the grid
      17              : !>      Joost VandeVondele (05.2002)
      18              : !>            1) prelim. introduction of the real space grid type
      19              : !>      JGH [30.08.02] multigrid arrays independent from potential
      20              : !>      JGH [17.07.03] distributed real space code
      21              : !>      JGH [23.11.03] refactoring and new loop ordering
      22              : !>      JGH [04.12.03] OpneMP parallelization of main loops
      23              : !>      Joost VandeVondele (12.2003)
      24              : !>           1) modified to compute tau
      25              : !>      Joost removed incremental build feature
      26              : !>      Joost introduced map consistent
      27              : !>      Rewrote grid integration/collocation routines, [Joost VandeVondele,03.2007]
      28              : !>      JGH [26.06.15] modification to allow for k-points
      29              : !> \author Matthias Krack (03.04.2001)
      30              : ! **************************************************************************************************
      31              : MODULE qs_integrate_potential_product
      32              :    USE ao_util,                         ONLY: exp_radius_very_extended
      33              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      34              :                                               get_atomic_kind_set
      35              :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      36              :                                               gto_basis_set_type
      37              :    USE block_p_types,                   ONLY: block_p_type
      38              :    USE cell_types,                      ONLY: cell_type,&
      39              :                                               pbc
      40              :    USE cp_control_types,                ONLY: dft_control_type
      41              :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      42              :                                               dbcsr_finalize,&
      43              :                                               dbcsr_get_block_p,&
      44              :                                               dbcsr_p_type,&
      45              :                                               dbcsr_type
      46              :    USE cp_dbcsr_operations,             ONLY: dbcsr_deallocate_matrix_set
      47              :    USE gaussian_gridlevels,             ONLY: gridlevel_info_type
      48              :    USE grid_api,                        ONLY: grid_integrate_task_list,&
      49              :                                               integrate_pgf_product
      50              :    USE kinds,                           ONLY: default_string_length,&
      51              :                                               dp
      52              :    USE message_passing,                 ONLY: mp_comm_type
      53              :    USE orbital_pointers,                ONLY: ncoset
      54              :    USE particle_types,                  ONLY: particle_type
      55              :    USE pw_env_types,                    ONLY: pw_env_get,&
      56              :                                               pw_env_type
      57              :    USE pw_types,                        ONLY: pw_r3d_rs_type
      58              :    USE qs_environment_types,            ONLY: get_qs_env,&
      59              :                                               qs_environment_type
      60              :    USE qs_force_types,                  ONLY: qs_force_type
      61              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      62              :                                               get_qs_kind_set,&
      63              :                                               qs_kind_type
      64              :    USE realspace_grid_types,            ONLY: realspace_grid_desc_p_type,&
      65              :                                               realspace_grid_type
      66              :    USE rs_pw_interface,                 ONLY: potential_pw2rs
      67              :    USE task_list_methods,               ONLY: rs_copy_to_buffer,&
      68              :                                               rs_copy_to_matrices,&
      69              :                                               rs_distribute_matrix,&
      70              :                                               rs_gather_matrices,&
      71              :                                               rs_scatter_matrices
      72              :    USE task_list_types,                 ONLY: atom_pair_type,&
      73              :                                               task_list_type,&
      74              :                                               task_type
      75              :    USE virial_types,                    ONLY: virial_type
      76              : 
      77              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
      78              : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
      79              : !$                    omp_init_lock, omp_set_lock, &
      80              : !$                    omp_unset_lock, omp_destroy_lock
      81              : #include "./base/base_uses.f90"
      82              : 
      83              :    IMPLICIT NONE
      84              : 
      85              :    PRIVATE
      86              : 
      87              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_integrate_potential_product'
      88              : 
      89              : ! *** Public subroutines ***
      90              : ! *** Don't include this routines directly, use the interface to
      91              : ! *** qs_integrate_potential
      92              : 
      93              :    PUBLIC :: integrate_v_rspace
      94              :    PUBLIC :: integrate_v_dbasis
      95              : 
      96              : CONTAINS
      97              : 
      98              : ! **************************************************************************************************
      99              : !> \brief Integrate a potential v_rspace over the derivatives of the basis functions
     100              : !>         < da/dR | V | b > + < a | V | db/dR >
     101              : !>        Adapted from the old version of integrate_v_rspace (ED)
     102              : !> \param v_rspace ...
     103              : !> \param matrix_vhxc_dbasis ...
     104              : !> \param matrix_p ...
     105              : !> \param qs_env ...
     106              : !> \param lambda The atom index.
     107              : ! **************************************************************************************************
     108           84 :    SUBROUTINE integrate_v_dbasis(v_rspace, matrix_vhxc_dbasis, matrix_p, qs_env, lambda)
     109              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_rspace
     110              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_vhxc_dbasis
     111              :       TYPE(dbcsr_type), POINTER                          :: matrix_p
     112              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     113              :       INTEGER, INTENT(IN)                                :: lambda
     114              : 
     115              :       CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_dbasis'
     116              : 
     117              :       INTEGER :: bcol, brow, handle, i, iatom, igrid_level, ikind, ikind_old, ilevel, img, ipair, &
     118              :          ipgf, ipgf_new, iset, iset_new, iset_old, itask, ithread, jatom, jkind, jkind_old, jpgf, &
     119              :          jpgf_new, jset, jset_new, jset_old, maxco, maxpgf, maxset, maxsgf_set, na1, na2, natom, &
     120              :          nb1, nb2, ncoa, ncob, nimages, nkind, nseta, nsetb, nthread, sgfa, sgfb
     121           84 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: block_touched
     122           84 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     123           84 :                                                             npgfb, nsgfa, nsgfb
     124           84 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     125              :       LOGICAL :: atom_pair_changed, atom_pair_done, dh_duplicated, distributed_grids, found, &
     126              :          my_compute_tau, new_set_pair_coming, pab_required, scatter, use_subpatch
     127              :       REAL(KIND=dp)                                      :: eps_rho_rspace, f, prefactor, radius, &
     128              :                                                             scalef, zetp
     129              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, ra, rab, rab_inv, rb, &
     130              :                                                             rp
     131              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: my_virial_a, my_virial_b
     132           84 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     133           84 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: h_block, hab, p_block, pab, rpgfa, &
     134           84 :                                                             rpgfb, sphi_a, sphi_b, work, zeta, zetb
     135           84 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: habt, hadb, hdab, pabt, workt
     136           84 :       REAL(kind=dp), DIMENSION(:, :, :, :), POINTER      :: hadbt, hdabt
     137           84 :       TYPE(atom_pair_type), DIMENSION(:), POINTER        :: atom_pair_recv, atom_pair_send
     138           84 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     139           84 :       TYPE(block_p_type), ALLOCATABLE, DIMENSION(:)      :: vhxc_block
     140              :       TYPE(cell_type), POINTER                           :: cell
     141           84 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: deltap
     142              :       TYPE(dft_control_type), POINTER                    :: dft_control
     143              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
     144              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
     145              :       TYPE(mp_comm_type)                                 :: group
     146           84 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     147              :       TYPE(pw_env_type), POINTER                         :: pw_env
     148           84 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     149           84 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     150              :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
     151           84 :          POINTER                                         :: rs_descs
     152           84 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_rho
     153              :       TYPE(task_list_type), POINTER                      :: task_list
     154           84 :       TYPE(task_type), DIMENSION(:), POINTER             :: tasks
     155              : 
     156           84 :       CALL timeset(routineN, handle)
     157           84 :       NULLIFY (pw_env)
     158              : 
     159              :       ! get the task lists
     160           84 :       CALL get_qs_env(qs_env=qs_env, task_list=task_list)
     161           84 :       CPASSERT(ASSOCIATED(task_list))
     162              : 
     163              :       ! the information on the grids is provided through pw_env
     164              :       ! pw_env has to be the parent env for the potential grid (input)
     165              :       ! there is an option to provide an external grid
     166           84 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     167              : 
     168              :       ! *** assign from pw_env
     169           84 :       gridlevel_info => pw_env%gridlevel_info
     170              : 
     171              :       ! get all the general information on the system we are working on
     172              :       CALL get_qs_env(qs_env=qs_env, &
     173              :                       atomic_kind_set=atomic_kind_set, &
     174              :                       qs_kind_set=qs_kind_set, &
     175              :                       cell=cell, &
     176              :                       natom=natom, &
     177              :                       dft_control=dft_control, &
     178           84 :                       particle_set=particle_set)
     179              : 
     180              :       ! GAPW not implemented here
     181           84 :       CPASSERT(.NOT. dft_control%qs_control%gapw)
     182              : 
     183              :       ! *** set up the rs multi-grids
     184           84 :       CPASSERT(ASSOCIATED(pw_env))
     185           84 :       CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
     186          258 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
     187          258 :          distributed_grids = rs_rho(igrid_level)%desc%distributed
     188              :       END DO
     189              :       ! get mpi group from rs_rho
     190           84 :       group = rs_rho(1)%desc%group
     191              : 
     192              :       ! transform the potential on the rs_multigrids
     193           84 :       CALL potential_pw2rs(rs_rho, v_rspace, pw_env)
     194              : 
     195           84 :       nkind = SIZE(qs_kind_set)
     196              : 
     197              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
     198              :                            maxco=maxco, &
     199              :                            maxsgf_set=maxsgf_set, &
     200           84 :                            basis_type="ORB")
     201              : 
     202              :       ! short cuts to task list variables
     203           84 :       tasks => task_list%tasks
     204           84 :       atom_pair_send => task_list%atom_pair_send
     205           84 :       atom_pair_recv => task_list%atom_pair_recv
     206              : 
     207              :       ! needs to be consistent with rho_rspace
     208           84 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     209              : 
     210              :       !   *** Initialize working density matrix ***
     211              :       ! distributed rs grids require a matrix that will be changed
     212              :       ! whereas this is not the case for replicated grids
     213          336 :       ALLOCATE (deltap(dft_control%nimages))
     214           84 :       IF (distributed_grids) THEN
     215              :          ! this matrix has no strict sparsity pattern in parallel
     216              :          ! deltap%sparsity_id=-1
     217            0 :          CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name="DeltaP")
     218              :       ELSE
     219           84 :          deltap(1)%matrix => matrix_p
     220              :       END IF
     221           84 :       nthread = 1
     222           84 : !$    nthread = omp_get_max_threads()
     223              : 
     224              :       !   *** Allocate work storage ***
     225           84 :       NULLIFY (pabt, habt, workt)
     226          420 :       ALLOCATE (habt(maxco, maxco, 0:nthread))
     227          420 :       ALLOCATE (workt(maxco, maxsgf_set, 0:nthread))
     228          420 :       ALLOCATE (hdabt(3, maxco, maxco, 0:nthread))
     229          336 :       ALLOCATE (hadbt(3, maxco, maxco, 0:nthread))
     230          336 :       ALLOCATE (pabt(maxco, maxco, 0:nthread))
     231              : 
     232           84 :       IF (distributed_grids) THEN
     233              :          CALL rs_distribute_matrix(rs_descs, deltap, atom_pair_send, atom_pair_recv, &
     234            0 :                                    nimages, scatter=.TRUE.)
     235              :       END IF
     236              : 
     237              : !$OMP PARALLEL DEFAULT(NONE), &
     238              : !$OMP SHARED(workt,habt,hdabt,hadbt,pabt,tasks,particle_set,natom,maxset), &
     239              : !$OMP SHARED(maxpgf,matrix_vhxc_dbasis,deltap), &
     240              : !$OMP SHARED(pab_required,ncoset,rs_rho,my_compute_tau), &
     241              : !$OMP SHARED(eps_rho_rspace,force,cell), &
     242              : !$OMP SHARED(gridlevel_info,task_list,block_touched,nthread,qs_kind_set), &
     243              : !$OMP SHARED(nimages,lambda, dh_duplicated), &
     244              : !$OMP PRIVATE(ithread,work,hab,hdab,hadb,pab,iset_old,jset_old), &
     245              : !$OMP PRIVATE(ikind_old,jkind_old,iatom,jatom,iset,jset,ikind,jkind,ilevel,ipgf,jpgf), &
     246              : !$OMP PRIVATE(img,brow,bcol,orb_basis_set,first_sgfa,la_max,la_min,npgfa,nseta,nsgfa), &
     247              : !$OMP PRIVATE(rpgfa,set_radius_a,sphi_a,zeta,first_sgfb,lb_max,lb_min,npgfb), &
     248              : !$OMP PRIVATE(nsetb,nsgfb,rpgfb,set_radius_b,sphi_b,zetb,found), &
     249              : !$OMP PRIVATE(force_a,force_b,my_virial_a,my_virial_b,atom_pair_changed,h_block, vhxc_block), &
     250              : !$OMP PRIVATE(p_block,ncoa,sgfa,ncob,sgfb,rab,ra,rb,rp,zetp,f,prefactor,radius,igrid_level), &
     251              : !$OMP PRIVATE(na1,na2,nb1,nb2,use_subpatch,rab_inv,new_set_pair_coming,atom_pair_done), &
     252              : !$OMP PRIVATE(iset_new,jset_new,ipgf_new,jpgf_new,scalef), &
     253           84 : !$OMP PRIVATE(itask)
     254              : 
     255              :       IF (.NOT. ALLOCATED(vhxc_block)) ALLOCATE (vhxc_block(3))
     256              : 
     257              :       ithread = 0
     258              : !$    ithread = omp_get_thread_num()
     259              :       work => workt(:, :, ithread)
     260              :       hab => habt(:, :, ithread)
     261              :       pab => pabt(:, :, ithread)
     262              :       hdab => hdabt(:, :, :, ithread)
     263              :       hadb => hadbt(:, :, :, ithread)
     264              : 
     265              :       iset_old = -1; jset_old = -1
     266              :       ikind_old = -1; jkind_old = -1
     267              : 
     268              :       ! Here we loop over gridlevels first, finalising the matrix after each grid level is
     269              :       ! completed.  On each grid level, we loop over atom pairs, which will only access
     270              :       ! a single block of each matrix, so with OpenMP, each matrix block is only touched
     271              :       ! by a single thread for each grid level
     272              :       loop_gridlevels: DO igrid_level = 1, gridlevel_info%ngrid_levels
     273              : !$OMP BARRIER
     274              : !$OMP DO schedule (dynamic, MAX(1,task_list%npairs(igrid_level)/(nthread*50)))
     275              :          loop_pairs: DO ipair = 1, task_list%npairs(igrid_level)
     276              :          loop_tasks: DO itask = task_list%taskstart(ipair, igrid_level), task_list%taskstop(ipair, igrid_level)
     277              :             ilevel = tasks(itask)%grid_level
     278              :             img = tasks(itask)%image
     279              :             iatom = tasks(itask)%iatom
     280              :             jatom = tasks(itask)%jatom
     281              :             iset = tasks(itask)%iset
     282              :             jset = tasks(itask)%jset
     283              :             ipgf = tasks(itask)%ipgf
     284              :             jpgf = tasks(itask)%jpgf
     285              : 
     286              :             ! At the start of a block of tasks, get atom data (and kind data, if needed)
     287              :             IF (itask .EQ. task_list%taskstart(ipair, igrid_level)) THEN
     288              : 
     289              :                ikind = particle_set(iatom)%atomic_kind%kind_number
     290              :                jkind = particle_set(jatom)%atomic_kind%kind_number
     291              : 
     292              :                ra(:) = pbc(particle_set(iatom)%r, cell)
     293              : 
     294              :                IF (iatom <= jatom) THEN
     295              :                   brow = iatom
     296              :                   bcol = jatom
     297              :                ELSE
     298              :                   brow = jatom
     299              :                   bcol = iatom
     300              :                END IF
     301              : 
     302              :                IF (ikind .NE. ikind_old) THEN
     303              :                   CALL get_qs_kind(qs_kind_set(ikind), &
     304              :                                    basis_set=orb_basis_set, basis_type="ORB")
     305              : 
     306              :                   CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     307              :                                          first_sgf=first_sgfa, &
     308              :                                          lmax=la_max, &
     309              :                                          lmin=la_min, &
     310              :                                          npgf=npgfa, &
     311              :                                          nset=nseta, &
     312              :                                          nsgf_set=nsgfa, &
     313              :                                          pgf_radius=rpgfa, &
     314              :                                          set_radius=set_radius_a, &
     315              :                                          sphi=sphi_a, &
     316              :                                          zet=zeta)
     317              :                END IF
     318              : 
     319              :                IF (jkind .NE. jkind_old) THEN
     320              :                   CALL get_qs_kind(qs_kind_set(jkind), &
     321              :                                    basis_set=orb_basis_set, basis_type="ORB")
     322              :                   CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     323              :                                          first_sgf=first_sgfb, &
     324              :                                          lmax=lb_max, &
     325              :                                          lmin=lb_min, &
     326              :                                          npgf=npgfb, &
     327              :                                          nset=nsetb, &
     328              :                                          nsgf_set=nsgfb, &
     329              :                                          pgf_radius=rpgfb, &
     330              :                                          set_radius=set_radius_b, &
     331              :                                          sphi=sphi_b, &
     332              :                                          zet=zetb)
     333              : 
     334              :                END IF
     335              : 
     336              :                DO i = 1, 3
     337              :                   NULLIFY (vhxc_block(i)%block)
     338              :                   CALL dbcsr_get_block_p(matrix_vhxc_dbasis(i)%matrix, brow, bcol, vhxc_block(i)%block, found)
     339              :                   CPASSERT(found)
     340              :                END DO
     341              : 
     342              :                CALL dbcsr_get_block_p(matrix=deltap(img)%matrix, &
     343              :                                       row=brow, col=bcol, BLOCK=p_block, found=found)
     344              :                CPASSERT(found)
     345              : 
     346              :                ikind_old = ikind
     347              :                jkind_old = jkind
     348              : 
     349              :                atom_pair_changed = .TRUE.
     350              : 
     351              :             ELSE
     352              : 
     353              :                atom_pair_changed = .FALSE.
     354              : 
     355              :             END IF
     356              : 
     357              :             IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN
     358              : 
     359              :                ncoa = npgfa(iset)*ncoset(la_max(iset))
     360              :                sgfa = first_sgfa(1, iset)
     361              :                ncob = npgfb(jset)*ncoset(lb_max(jset))
     362              :                sgfb = first_sgfb(1, jset)
     363              : 
     364              :                IF (iatom <= jatom) THEN
     365              :                   work(1:ncoa, 1:nsgfb(jset)) = MATMUL(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
     366              :                                                        p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1))
     367              :                   pab(1:ncoa, 1:ncob) = MATMUL(work(1:ncoa, 1:nsgfb(jset)), TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
     368              :                ELSE
     369              :                   work(1:ncob, 1:nsgfa(iset)) = MATMUL(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1), &
     370              :                                                        p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1))
     371              :                   pab(1:ncob, 1:ncoa) = MATMUL(work(1:ncob, 1:nsgfa(iset)), TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)))
     372              :                END IF
     373              : 
     374              :                IF (iatom <= jatom) THEN
     375              :                   hab(1:ncoa, 1:ncob) = 0._dp
     376              :                   hdab(:, 1:ncoa, 1:ncob) = 0._dp
     377              :                   hadb(:, 1:ncoa, 1:ncob) = 0._dp
     378              :                ELSE
     379              :                   hab(1:ncob, 1:ncoa) = 0._dp
     380              :                   hdab(:, 1:ncob, 1:ncoa) = 0._dp
     381              :                   hadb(:, 1:ncob, 1:ncoa) = 0._dp
     382              :                END IF
     383              : 
     384              :                iset_old = iset
     385              :                jset_old = jset
     386              : 
     387              :             END IF
     388              : 
     389              :             rab = tasks(itask)%rab
     390              :             rb(:) = ra(:) + rab(:)
     391              :             zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
     392              : 
     393              :             f = zetb(jpgf, jset)/zetp
     394              :             rp(:) = ra(:) + f*rab(:)
     395              :             prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
     396              :             radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
     397              :                                               lb_min=lb_min(jset), lb_max=lb_max(jset), &
     398              :                                               ra=ra, rb=rb, rp=rp, &
     399              :                                               zetp=zetp, eps=eps_rho_rspace, &
     400              :                                               prefactor=prefactor, cutoff=1.0_dp)
     401              : 
     402              :             na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
     403              :             na2 = ipgf*ncoset(la_max(iset))
     404              :             nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
     405              :             nb2 = jpgf*ncoset(lb_max(jset))
     406              : 
     407              :             ! check whether we need to use fawzi's generalised collocation scheme
     408              :             IF (rs_rho(igrid_level)%desc%distributed) THEN
     409              :                !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks
     410              :                IF (tasks(itask)%dist_type .EQ. 2) THEN
     411              :                   use_subpatch = .TRUE.
     412              :                ELSE
     413              :                   use_subpatch = .FALSE.
     414              :                END IF
     415              :             ELSE
     416              :                use_subpatch = .FALSE.
     417              :             END IF
     418              : 
     419              :             IF (iatom <= jatom) THEN
     420              :                IF (iatom == lambda) &
     421              :                   CALL integrate_pgf_product( &
     422              :                   la_max(iset), zeta(ipgf, iset), la_min(iset), &
     423              :                   lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
     424              :                   ra, rab, rs_rho(igrid_level), &
     425              :                   hab, o1=na1 - 1, o2=nb1 - 1, &
     426              :                   radius=radius, &
     427              :                   calculate_forces=.TRUE., &
     428              :                   compute_tau=.FALSE., &
     429              :                   use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern, &
     430              :                   hdab=hdab, pab=pab)
     431              :                IF (jatom == lambda) &
     432              :                   CALL integrate_pgf_product( &
     433              :                   la_max(iset), zeta(ipgf, iset), la_min(iset), &
     434              :                   lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
     435              :                   ra, rab, rs_rho(igrid_level), &
     436              :                   hab, o1=na1 - 1, o2=nb1 - 1, &
     437              :                   radius=radius, &
     438              :                   calculate_forces=.TRUE., &
     439              :                   compute_tau=.FALSE., &
     440              :                   use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern, &
     441              :                   hadb=hadb, pab=pab)
     442              :             ELSE
     443              :                rab_inv = -rab
     444              :                IF (iatom == lambda) &
     445              :                   CALL integrate_pgf_product( &
     446              :                   lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
     447              :                   la_max(iset), zeta(ipgf, iset), la_min(iset), &
     448              :                   rb, rab_inv, rs_rho(igrid_level), &
     449              :                   hab, o1=nb1 - 1, o2=na1 - 1, &
     450              :                   radius=radius, &
     451              :                   calculate_forces=.TRUE., &
     452              :                   force_a=force_b, force_b=force_a, &
     453              :                   compute_tau=.FALSE., &
     454              :                   use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern, &
     455              :                   hadb=hadb, pab=pab)
     456              :                IF (jatom == lambda) &
     457              :                   CALL integrate_pgf_product( &
     458              :                   lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
     459              :                   la_max(iset), zeta(ipgf, iset), la_min(iset), &
     460              :                   rb, rab_inv, rs_rho(igrid_level), &
     461              :                   hab, o1=nb1 - 1, o2=na1 - 1, &
     462              :                   radius=radius, &
     463              :                   calculate_forces=.TRUE., &
     464              :                   force_a=force_b, force_b=force_a, &
     465              :                   compute_tau=.FALSE., &
     466              :                   use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern, &
     467              :                   hdab=hdab, pab=pab)
     468              :             END IF
     469              : 
     470              :             new_set_pair_coming = .FALSE.
     471              :             atom_pair_done = .FALSE.
     472              :             IF (itask < task_list%taskstop(ipair, igrid_level)) THEN
     473              :                ilevel = tasks(itask + 1)%grid_level
     474              :                img = tasks(itask + 1)%image
     475              :                iatom = tasks(itask + 1)%iatom
     476              :                jatom = tasks(itask + 1)%jatom
     477              :                iset_new = tasks(itask + 1)%iset
     478              :                jset_new = tasks(itask + 1)%jset
     479              :                ipgf_new = tasks(itask + 1)%ipgf
     480              :                jpgf_new = tasks(itask + 1)%jpgf
     481              :                IF (iset_new .NE. iset .OR. jset_new .NE. jset) THEN
     482              :                   new_set_pair_coming = .TRUE.
     483              :                END IF
     484              :             ELSE
     485              :                ! do not forget the last block
     486              :                new_set_pair_coming = .TRUE.
     487              :                atom_pair_done = .TRUE.
     488              :             END IF
     489              : 
     490              :             IF (new_set_pair_coming) THEN
     491              : 
     492              :                DO i = 1, 3
     493              :                   hdab(i, :, :) = hdab(i, :, :) + hadb(i, :, :)
     494              :                   IF (iatom <= jatom) THEN
     495              :                      work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hdab(i, 1:ncoa, 1:ncob), sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
     496              :                      vhxc_block(i)%block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
     497              :                         vhxc_block(i)%block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
     498              :                         MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
     499              :                   ELSE
     500              :                      work(1:ncob, 1:nsgfa(iset)) = MATMUL(hdab(i, 1:ncob, 1:ncoa), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
     501              :                      vhxc_block(i)%block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
     502              :                         vhxc_block(i)%block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
     503              :                         MATMUL(TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)), work(1:ncob, 1:nsgfa(iset)))
     504              :                   END IF
     505              :                END DO
     506              :             END IF  ! new_set_pair_coming
     507              : 
     508              :          END DO loop_tasks
     509              :          END DO loop_pairs
     510              : !$OMP END DO
     511              : 
     512              :          DO i = 1, 3
     513              :             CALL dbcsr_finalize(matrix_vhxc_dbasis(i)%matrix)
     514              :          END DO
     515              : 
     516              :       END DO loop_gridlevels
     517              : 
     518              : !$OMP END PARALLEL
     519              : 
     520           84 :       IF (distributed_grids) THEN
     521              :          ! Reconstruct H matrix if using distributed RS grids
     522              :          ! note send and recv direction reversed WRT collocate
     523            0 :          scatter = .FALSE.
     524              :          CALL rs_distribute_matrix(rs_descs, matrix_vhxc_dbasis, atom_pair_recv, atom_pair_send, &
     525            0 :                                    dft_control%nimages, scatter=.FALSE.)
     526              :       END IF
     527              : 
     528              :       IF (distributed_grids) THEN
     529            0 :          CALL dbcsr_deallocate_matrix_set(deltap)
     530              :       ELSE
     531          168 :          DO img = 1, dft_control%nimages
     532          168 :             NULLIFY (deltap(img)%matrix)
     533              :          END DO
     534           84 :          DEALLOCATE (deltap)
     535              :       END IF
     536              : 
     537           84 :       DEALLOCATE (pabt, habt, workt, hdabt, hadbt)
     538              : 
     539           84 :       CALL timestop(handle)
     540          252 :    END SUBROUTINE integrate_v_dbasis
     541              : 
     542              : ! **************************************************************************************************
     543              : !> \brief computes matrix elements corresponding to a given potential
     544              : !> \param v_rspace ...
     545              : !> \param hmat ...
     546              : !> \param hmat_kp ...
     547              : !> \param pmat ...
     548              : !> \param pmat_kp ...
     549              : !> \param qs_env ...
     550              : !> \param calculate_forces ...
     551              : !> \param force_adm optional scaling of force
     552              : !> \param compute_tau ...
     553              : !> \param gapw ...
     554              : !> \param basis_type ...
     555              : !> \param pw_env_external ...
     556              : !> \param task_list_external ...
     557              : !> \par History
     558              : !>      IAB (29-Apr-2010): Added OpenMP parallelisation to task loop
     559              : !>                         (c) The Numerical Algorithms Group (NAG) Ltd, 2010 on behalf of the HECToR project
     560              : !>      Some refactoring, get priorities for options correct (JGH, 04.2014)
     561              : !>      Added options to allow for k-points
     562              : !>      For a smooth transition we allow for old and new (vector) matrices (hmat, pmat) (JGH, 06.2015)
     563              : !> \note
     564              : !>     integrates a given potential (or other object on a real
     565              : !>     space grid) = v_rspace using a multi grid technique (mgrid_*)
     566              : !>     over the basis set producing a number for every element of h
     567              : !>     (should have the same sparsity structure of S)
     568              : !>     additional screening is available using the magnitude of the
     569              : !>     elements in p (? I'm not sure this is a very good idea)
     570              : !>     this argument is optional
     571              : !>     derivatives of these matrix elements with respect to the ionic
     572              : !>     coordinates can be computed as well
     573              : ! **************************************************************************************************
     574       185798 :    SUBROUTINE integrate_v_rspace(v_rspace, hmat, hmat_kp, pmat, pmat_kp, &
     575              :                                  qs_env, calculate_forces, force_adm, &
     576              :                                  compute_tau, gapw, basis_type, pw_env_external, task_list_external)
     577              : 
     578              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_rspace
     579              :       TYPE(dbcsr_p_type), INTENT(INOUT), OPTIONAL        :: hmat
     580              :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
     581              :          POINTER                                         :: hmat_kp
     582              :       TYPE(dbcsr_p_type), INTENT(IN), OPTIONAL           :: pmat
     583              :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
     584              :          POINTER                                         :: pmat_kp
     585              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     586              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     587              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: force_adm
     588              :       LOGICAL, INTENT(IN), OPTIONAL                      :: compute_tau, gapw
     589              :       CHARACTER(len=*), INTENT(IN), OPTIONAL             :: basis_type
     590              :       TYPE(pw_env_type), OPTIONAL, POINTER               :: pw_env_external
     591              :       TYPE(task_list_type), OPTIONAL, POINTER            :: task_list_external
     592              : 
     593              :       CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace'
     594              : 
     595              :       CHARACTER(len=default_string_length)               :: my_basis_type
     596              :       INTEGER                                            :: atom_a, handle, iatom, igrid_level, &
     597              :                                                             ikind, img, maxco, maxsgf_set, natoms, &
     598              :                                                             nimages, nkind
     599       185798 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     600              :       LOGICAL                                            :: calculate_virial, distributed_grids, &
     601              :                                                             do_kp, my_compute_tau, my_gapw, &
     602              :                                                             pab_required
     603              :       REAL(KIND=dp)                                      :: admm_scal_fac
     604       185798 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: forces_array
     605              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: virial_matrix
     606       185798 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     607              :       TYPE(cell_type), POINTER                           :: cell
     608       185798 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: deltap, dhmat
     609              :       TYPE(dft_control_type), POINTER                    :: dft_control
     610              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
     611              :       TYPE(mp_comm_type)                                 :: group
     612       185798 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     613              :       TYPE(pw_env_type), POINTER                         :: pw_env
     614       185798 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     615       185798 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     616       185798 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_v
     617              :       TYPE(task_list_type), POINTER                      :: task_list, task_list_soft
     618              :       TYPE(virial_type), POINTER                         :: virial
     619              : 
     620       185798 :       CALL timeset(routineN, handle)
     621              : 
     622              :       ! we test here if the provided operator matrices are consistent
     623       185798 :       CPASSERT(PRESENT(hmat) .OR. PRESENT(hmat_kp))
     624       185798 :       do_kp = .FALSE.
     625       185798 :       IF (PRESENT(hmat_kp)) do_kp = .TRUE.
     626       185798 :       IF (PRESENT(pmat)) THEN
     627        32016 :          CPASSERT(PRESENT(hmat))
     628       153782 :       ELSE IF (PRESENT(pmat_kp)) THEN
     629       124928 :          CPASSERT(PRESENT(hmat_kp))
     630              :       END IF
     631              : 
     632       185798 :       NULLIFY (pw_env)
     633              : 
     634              :       ! this routine works in two modes:
     635              :       ! normal mode : <a| V | b>
     636              :       ! tau mode    : < nabla a| V | nabla b>
     637       185798 :       my_compute_tau = .FALSE.
     638       185798 :       IF (PRESENT(compute_tau)) my_compute_tau = compute_tau
     639              : 
     640              :       ! this sets the basis set to be used. GAPW(==soft basis) overwrites basis_type
     641              :       ! default is "ORB"
     642       185798 :       my_gapw = .FALSE.
     643       185798 :       IF (PRESENT(gapw)) my_gapw = gapw
     644       185798 :       IF (PRESENT(basis_type)) THEN
     645        13064 :          my_basis_type = basis_type
     646              :       ELSE
     647       172734 :          my_basis_type = "ORB"
     648              :       END IF
     649              : 
     650              :       ! get the task lists
     651              :       ! task lists have to be in sync with basis sets
     652              :       ! there is an option to provide the task list from outside (not through qs_env)
     653              :       ! outside option has highest priority
     654              :       CALL get_qs_env(qs_env=qs_env, &
     655              :                       task_list=task_list, &
     656       185798 :                       task_list_soft=task_list_soft)
     657       185798 :       IF (.NOT. my_basis_type == "ORB") THEN
     658        12464 :          CPASSERT(PRESENT(task_list_external))
     659              :       END IF
     660       185798 :       IF (my_gapw) task_list => task_list_soft
     661       185798 :       IF (PRESENT(task_list_external)) task_list => task_list_external
     662       185798 :       CPASSERT(ASSOCIATED(task_list))
     663              : 
     664              :       ! the information on the grids is provided through pw_env
     665              :       ! pw_env has to be the parent env for the potential grid (input)
     666              :       ! there is an option to provide an external grid
     667       185798 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     668       185798 :       IF (PRESENT(pw_env_external)) pw_env => pw_env_external
     669              : 
     670              :       ! get all the general information on the system we are working on
     671              :       CALL get_qs_env(qs_env=qs_env, &
     672              :                       atomic_kind_set=atomic_kind_set, &
     673              :                       qs_kind_set=qs_kind_set, &
     674              :                       cell=cell, &
     675              :                       natom=natoms, &
     676              :                       dft_control=dft_control, &
     677              :                       particle_set=particle_set, &
     678              :                       force=force, &
     679       185798 :                       virial=virial)
     680              : 
     681       185798 :       admm_scal_fac = 1.0_dp
     682       185798 :       IF (PRESENT(force_adm)) admm_scal_fac = force_adm
     683              : 
     684       185798 :       CPASSERT(ASSOCIATED(pw_env))
     685       185798 :       CALL pw_env_get(pw_env, rs_grids=rs_v)
     686              : 
     687              :       ! get mpi group from rs_v
     688       185798 :       group = rs_v(1)%desc%group
     689              : 
     690              :       ! assign from pw_env
     691       185798 :       gridlevel_info => pw_env%gridlevel_info
     692              : 
     693              :       ! transform the potential on the rs_multigrids
     694       185798 :       CALL potential_pw2rs(rs_v, v_rspace, pw_env)
     695              : 
     696       185798 :       nimages = dft_control%nimages
     697       185798 :       IF (nimages > 1) THEN
     698         1512 :          CPASSERT(do_kp)
     699              :       END IF
     700       185798 :       nkind = SIZE(qs_kind_set)
     701       185798 :       calculate_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. calculate_forces
     702       185798 :       pab_required = (PRESENT(pmat) .OR. PRESENT(pmat_kp)) .AND. calculate_forces
     703              : 
     704              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
     705              :                            maxco=maxco, &
     706              :                            maxsgf_set=maxsgf_set, &
     707       185798 :                            basis_type=my_basis_type)
     708              : 
     709       185798 :       distributed_grids = .FALSE.
     710       921044 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
     711       921044 :          IF (rs_v(igrid_level)%desc%distributed) THEN
     712          230 :             distributed_grids = .TRUE.
     713              :          END IF
     714              :       END DO
     715              : 
     716       557394 :       ALLOCATE (forces_array(3, natoms))
     717              : 
     718       185798 :       IF (pab_required) THEN
     719              :          ! initialize the working pmat structures
     720       106882 :          ALLOCATE (deltap(nimages))
     721        23301 :          IF (do_kp) THEN
     722        26844 :             DO img = 1, nimages
     723        26844 :                deltap(img)%matrix => pmat_kp(img)%matrix
     724              :             END DO
     725              :          ELSE
     726        16718 :             deltap(1)%matrix => pmat%matrix
     727              :          END IF
     728              : 
     729              :          ! Distribute matrix blocks.
     730        23301 :          IF (distributed_grids) THEN
     731           36 :             CALL rs_scatter_matrices(deltap, task_list%pab_buffer, task_list, group)
     732              :          ELSE
     733        23265 :             CALL rs_copy_to_buffer(deltap, task_list%pab_buffer, task_list)
     734              :          END IF
     735        23301 :          DEALLOCATE (deltap)
     736              :       END IF
     737              : 
     738              :       ! Map all tasks from the grids
     739              :       CALL grid_integrate_task_list(task_list=task_list%grid_task_list, &
     740              :                                     compute_tau=my_compute_tau, &
     741              :                                     calculate_forces=calculate_forces, &
     742              :                                     calculate_virial=calculate_virial, &
     743              :                                     pab_blocks=task_list%pab_buffer, &
     744              :                                     rs_grids=rs_v, &
     745              :                                     hab_blocks=task_list%hab_buffer, &
     746              :                                     forces=forces_array, &
     747       185798 :                                     virial=virial_matrix)
     748              : 
     749       185798 :       IF (calculate_forces) THEN
     750        23301 :          CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
     751              : !$OMP PARALLEL DO DEFAULT(NONE)  PRIVATE(atom_a, ikind) &
     752        23301 : !$OMP             SHARED(natoms, force, forces_array, atom_of_kind, kind_of, admm_scal_fac)
     753              :          DO iatom = 1, natoms
     754              :             atom_a = atom_of_kind(iatom)
     755              :             ikind = kind_of(iatom)
     756              :             force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + admm_scal_fac*forces_array(:, iatom)
     757              :          END DO
     758              : !$OMP END PARALLEL DO
     759        23301 :          DEALLOCATE (atom_of_kind, kind_of)
     760              :       END IF
     761              : 
     762       185798 :       IF (calculate_virial) THEN
     763        46917 :          virial%pv_virial = virial%pv_virial + admm_scal_fac*virial_matrix
     764              :       END IF
     765              : 
     766              :       ! Gather all matrix images into a single array.
     767       859922 :       ALLOCATE (dhmat(nimages))
     768       185798 :       IF (PRESENT(hmat_kp)) THEN
     769       124928 :          CPASSERT(.NOT. PRESENT(hmat))
     770       366586 :          DO img = 1, nimages
     771       366586 :             dhmat(img)%matrix => hmat_kp(img)%matrix
     772              :          END DO
     773              :       ELSE
     774        60870 :          CPASSERT(PRESENT(hmat) .AND. nimages == 1)
     775        60870 :          dhmat(1)%matrix => hmat%matrix
     776              :       END IF
     777              : 
     778              :       ! Distribute matrix blocks.
     779       185798 :       IF (distributed_grids) THEN
     780          218 :          CALL rs_gather_matrices(task_list%hab_buffer, dhmat, task_list, group)
     781              :       ELSE
     782       185580 :          CALL rs_copy_to_matrices(task_list%hab_buffer, dhmat, task_list)
     783              :       END IF
     784       185798 :       DEALLOCATE (dhmat)
     785              : 
     786       185798 :       CALL timestop(handle)
     787              : 
     788       557394 :    END SUBROUTINE integrate_v_rspace
     789              : 
     790              : END MODULE qs_integrate_potential_product
        

Generated by: LCOV version 2.0-1