LCOV - code coverage report
Current view: top level - src - qs_integrate_potential_product.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 144 149 96.6 %
Date: 2024-03-28 07:31:50 Functions: 2 2 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief 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_operations,             ONLY: dbcsr_deallocate_matrix_set
      42             :    USE dbcsr_api,                       ONLY: dbcsr_copy,&
      43             :                                               dbcsr_finalize,&
      44             :                                               dbcsr_get_block_p,&
      45             :                                               dbcsr_p_type,&
      46             :                                               dbcsr_type
      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          42 :    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          42 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: block_touched
     122          42 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     123          42 :                                                             npgfb, nsgfa, nsgfb
     124          42 :       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          42 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     133          42 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: h_block, hab, p_block, pab, rpgfa, &
     134          42 :                                                             rpgfb, sphi_a, sphi_b, work, zeta, zetb
     135          42 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: habt, hadb, hdab, pabt, workt
     136          42 :       REAL(kind=dp), DIMENSION(:, :, :, :), POINTER      :: hadbt, hdabt
     137          42 :       TYPE(atom_pair_type), DIMENSION(:), POINTER        :: atom_pair_recv, atom_pair_send
     138          42 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     139          42 :       TYPE(block_p_type), ALLOCATABLE, DIMENSION(:)      :: vhxc_block
     140             :       TYPE(cell_type), POINTER                           :: cell
     141          42 :       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          42 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     147             :       TYPE(pw_env_type), POINTER                         :: pw_env
     148          42 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     149          42 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     150             :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
     151          42 :          POINTER                                         :: rs_descs
     152          42 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_rho
     153             :       TYPE(task_list_type), POINTER                      :: task_list
     154          42 :       TYPE(task_type), DIMENSION(:), POINTER             :: tasks
     155             : 
     156          42 :       CALL timeset(routineN, handle)
     157          42 :       NULLIFY (pw_env)
     158             : 
     159             :       ! get the task lists
     160          42 :       CALL get_qs_env(qs_env=qs_env, task_list=task_list)
     161          42 :       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          42 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     167             : 
     168             :       ! *** assign from pw_env
     169          42 :       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          42 :                       particle_set=particle_set)
     179             : 
     180             :       ! GAPW not implemented here
     181          42 :       CPASSERT(.NOT. dft_control%qs_control%gapw)
     182             : 
     183             :       ! *** set up the rs multi-grids
     184          42 :       CPASSERT(ASSOCIATED(pw_env))
     185          42 :       CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
     186         138 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
     187         138 :          distributed_grids = rs_rho(igrid_level)%desc%distributed
     188             :       END DO
     189             :       ! get mpi group from rs_rho
     190          42 :       group = rs_rho(1)%desc%group
     191             : 
     192             :       ! transform the potential on the rs_multigrids
     193          42 :       CALL potential_pw2rs(rs_rho, v_rspace, pw_env)
     194             : 
     195          42 :       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          42 :                            basis_type="ORB")
     201             : 
     202             :       ! short cuts to task list variables
     203          42 :       tasks => task_list%tasks
     204          42 :       atom_pair_send => task_list%atom_pair_send
     205          42 :       atom_pair_recv => task_list%atom_pair_recv
     206             : 
     207             :       ! needs to be consistent with rho_rspace
     208          42 :       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         168 :       ALLOCATE (deltap(dft_control%nimages))
     214          42 :       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          42 :          deltap(1)%matrix => matrix_p
     220             :       END IF
     221          42 :       nthread = 1
     222          42 : !$    nthread = omp_get_max_threads()
     223             : 
     224             :       !   *** Allocate work storage ***
     225          42 :       NULLIFY (pabt, habt, workt)
     226         210 :       ALLOCATE (habt(maxco, maxco, 0:nthread))
     227         210 :       ALLOCATE (workt(maxco, maxsgf_set, 0:nthread))
     228         210 :       ALLOCATE (hdabt(3, maxco, maxco, 0:nthread))
     229         210 :       ALLOCATE (hadbt(3, maxco, maxco, 0:nthread))
     230         210 :       ALLOCATE (pabt(maxco, maxco, 0:nthread))
     231             : 
     232          42 :       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          42 : !$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          42 :       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          42 :       IF (distributed_grids) THEN
     529           0 :          CALL dbcsr_deallocate_matrix_set(deltap)
     530             :       ELSE
     531          84 :          DO img = 1, dft_control%nimages
     532          84 :             NULLIFY (deltap(img)%matrix)
     533             :          END DO
     534          42 :          DEALLOCATE (deltap)
     535             :       END IF
     536             : 
     537          42 :       DEALLOCATE (pabt, habt, workt, hdabt, hadbt)
     538             : 
     539          42 :       CALL timestop(handle)
     540         126 :    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      178779 :    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      178779 :       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      178779 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: forces_array
     605             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: virial_matrix
     606      178779 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     607             :       TYPE(cell_type), POINTER                           :: cell
     608      178779 :       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      178779 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     613             :       TYPE(pw_env_type), POINTER                         :: pw_env
     614      178779 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     615      178779 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     616      178779 :       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      178779 :       CALL timeset(routineN, handle)
     621             : 
     622             :       ! we test here if the provided operator matrices are consistent
     623      178779 :       CPASSERT(PRESENT(hmat) .OR. PRESENT(hmat_kp))
     624      178779 :       do_kp = .FALSE.
     625      178779 :       IF (PRESENT(hmat_kp)) do_kp = .TRUE.
     626      178779 :       IF (PRESENT(pmat)) THEN
     627       31604 :          CPASSERT(PRESENT(hmat))
     628      147175 :       ELSE IF (PRESENT(pmat_kp)) THEN
     629      118781 :          CPASSERT(PRESENT(hmat_kp))
     630             :       END IF
     631             : 
     632      178779 :       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      178779 :       my_compute_tau = .FALSE.
     638      178779 :       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      178779 :       my_gapw = .FALSE.
     643      178779 :       IF (PRESENT(gapw)) my_gapw = gapw
     644      178779 :       IF (PRESENT(basis_type)) THEN
     645       12484 :          my_basis_type = basis_type
     646             :       ELSE
     647      166295 :          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      178779 :                       task_list_soft=task_list_soft)
     657      178779 :       IF (.NOT. my_basis_type == "ORB") THEN
     658       11884 :          CPASSERT(PRESENT(task_list_external))
     659             :       END IF
     660      178779 :       IF (my_gapw) task_list => task_list_soft
     661      178779 :       IF (PRESENT(task_list_external)) task_list => task_list_external
     662      178779 :       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      178779 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     668      178779 :       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      178779 :                       virial=virial)
     680             : 
     681      178779 :       admm_scal_fac = 1.0_dp
     682      178779 :       IF (PRESENT(force_adm)) admm_scal_fac = force_adm
     683             : 
     684      178779 :       CPASSERT(ASSOCIATED(pw_env))
     685      178779 :       CALL pw_env_get(pw_env, rs_grids=rs_v)
     686             : 
     687             :       ! get mpi group from rs_v
     688      178779 :       group = rs_v(1)%desc%group
     689             : 
     690             :       ! assign from pw_env
     691      178779 :       gridlevel_info => pw_env%gridlevel_info
     692             : 
     693             :       ! transform the potential on the rs_multigrids
     694      178779 :       CALL potential_pw2rs(rs_v, v_rspace, pw_env)
     695             : 
     696      178779 :       nimages = dft_control%nimages
     697      178779 :       IF (nimages > 1) THEN
     698        1318 :          CPASSERT(do_kp)
     699             :       END IF
     700      178779 :       nkind = SIZE(qs_kind_set)
     701      178779 :       calculate_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. calculate_forces
     702      178779 :       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      178779 :                            basis_type=my_basis_type)
     708             : 
     709      178779 :       distributed_grids = .FALSE.
     710      888179 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
     711      888179 :          IF (rs_v(igrid_level)%desc%distributed) THEN
     712         230 :             distributed_grids = .TRUE.
     713             :          END IF
     714             :       END DO
     715             : 
     716      536337 :       ALLOCATE (forces_array(3, natoms))
     717             : 
     718      178779 :       IF (pab_required) THEN
     719             :          ! initialize the working pmat structures
     720      103014 :          ALLOCATE (deltap(nimages))
     721       22521 :          IF (do_kp) THEN
     722       25692 :             DO img = 1, nimages
     723       25692 :                deltap(img)%matrix => pmat_kp(img)%matrix
     724             :             END DO
     725             :          ELSE
     726       16140 :             deltap(1)%matrix => pmat%matrix
     727             :          END IF
     728             : 
     729             :          ! Distribute matrix blocks.
     730       22521 :          IF (distributed_grids) THEN
     731          36 :             CALL rs_scatter_matrices(deltap, task_list%pab_buffer, task_list, group)
     732             :          ELSE
     733       22485 :             CALL rs_copy_to_buffer(deltap, task_list%pab_buffer, task_list)
     734             :          END IF
     735       22521 :          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      178779 :                                     virial=virial_matrix)
     748             : 
     749      178779 :       IF (calculate_forces) THEN
     750       22521 :          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       22521 : !$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       22521 :          DEALLOCATE (atom_of_kind, kind_of)
     760             :       END IF
     761             : 
     762      178779 :       IF (calculate_virial) THEN
     763       43901 :          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      821406 :       ALLOCATE (dhmat(nimages))
     768      178779 :       IF (PRESENT(hmat_kp)) THEN
     769      118781 :          CPASSERT(.NOT. PRESENT(hmat))
     770      343852 :          DO img = 1, nimages
     771      343852 :             dhmat(img)%matrix => hmat_kp(img)%matrix
     772             :          END DO
     773             :       ELSE
     774       59998 :          CPASSERT(PRESENT(hmat) .AND. nimages == 1)
     775       59998 :          dhmat(1)%matrix => hmat%matrix
     776             :       END IF
     777             : 
     778             :       ! Distribute matrix blocks.
     779      178779 :       IF (distributed_grids) THEN
     780         218 :          CALL rs_gather_matrices(task_list%hab_buffer, dhmat, task_list, group)
     781             :       ELSE
     782      178561 :          CALL rs_copy_to_matrices(task_list%hab_buffer, dhmat, task_list)
     783             :       END IF
     784      178779 :       DEALLOCATE (dhmat)
     785             : 
     786      178779 :       CALL timestop(handle)
     787             : 
     788      536337 :    END SUBROUTINE integrate_v_rspace
     789             : 
     790             : END MODULE qs_integrate_potential_product

Generated by: LCOV version 1.15