LCOV - code coverage report
Current view: top level - src - qs_cdft_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 91.5 % 846 774
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

            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 Subroutines for building CDFT constraints
      10              : !> \par   History
      11              : !>                 separated from et_coupling [03.2017]
      12              : !> \author Nico Holmberg [03.2017]
      13              : ! **************************************************************************************************
      14              : MODULE qs_cdft_methods
      15              :    USE ao_util,                         ONLY: exp_radius_very_extended
      16              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      17              :                                               get_atomic_kind,&
      18              :                                               get_atomic_kind_set
      19              :    USE cell_types,                      ONLY: cell_type,&
      20              :                                               pbc
      21              :    USE cp_control_types,                ONLY: dft_control_type
      22              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23              :                                               cp_logger_type
      24              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      25              :                                               cp_print_key_unit_nr
      26              :    USE cp_realspace_grid_cube,          ONLY: cp_cube_to_pw
      27              :    USE grid_api,                        ONLY: GRID_FUNC_AB,&
      28              :                                               collocate_pgf_product
      29              :    USE hirshfeld_types,                 ONLY: hirshfeld_type
      30              :    USE input_constants,                 ONLY: cdft_alpha_constraint,&
      31              :                                               cdft_beta_constraint,&
      32              :                                               cdft_charge_constraint,&
      33              :                                               cdft_magnetization_constraint,&
      34              :                                               outer_scf_becke_constraint,&
      35              :                                               outer_scf_hirshfeld_constraint
      36              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      37              :                                               section_vals_type
      38              :    USE kahan_sum,                       ONLY: accurate_dot_product
      39              :    USE kinds,                           ONLY: dp
      40              :    USE message_passing,                 ONLY: mp_para_env_type
      41              :    USE particle_types,                  ONLY: particle_type
      42              :    USE pw_env_types,                    ONLY: pw_env_get,&
      43              :                                               pw_env_type
      44              :    USE pw_methods,                      ONLY: pw_axpy,&
      45              :                                               pw_copy,&
      46              :                                               pw_integral_ab,&
      47              :                                               pw_integrate_function,&
      48              :                                               pw_set,&
      49              :                                               pw_zero
      50              :    USE pw_pool_types,                   ONLY: pw_pool_type
      51              :    USE pw_types,                        ONLY: pw_r3d_rs_type
      52              :    USE qs_cdft_types,                   ONLY: becke_constraint_type,&
      53              :                                               cdft_control_type,&
      54              :                                               cdft_group_type,&
      55              :                                               hirshfeld_constraint_type
      56              :    USE qs_cdft_utils,                   ONLY: becke_constraint_init,&
      57              :                                               cdft_constraint_print,&
      58              :                                               cdft_print_hirshfeld_density,&
      59              :                                               hfun_scale,&
      60              :                                               hirshfeld_constraint_init
      61              :    USE qs_energy_types,                 ONLY: qs_energy_type
      62              :    USE qs_environment_types,            ONLY: get_qs_env,&
      63              :                                               qs_environment_type
      64              :    USE qs_force_types,                  ONLY: qs_force_type
      65              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      66              :                                               qs_kind_type
      67              :    USE qs_rho0_types,                   ONLY: get_rho0_mpole,&
      68              :                                               mpole_rho_atom,&
      69              :                                               rho0_mpole_type
      70              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      71              :                                               qs_rho_type
      72              :    USE qs_subsys_types,                 ONLY: qs_subsys_get,&
      73              :                                               qs_subsys_type
      74              :    USE realspace_grid_types,            ONLY: realspace_grid_desc_type,&
      75              :                                               realspace_grid_type,&
      76              :                                               rs_grid_create,&
      77              :                                               rs_grid_release,&
      78              :                                               rs_grid_zero,&
      79              :                                               transfer_rs2pw
      80              : #include "./base/base_uses.f90"
      81              : 
      82              :    IMPLICIT NONE
      83              : 
      84              :    PRIVATE
      85              : 
      86              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cdft_methods'
      87              :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      88              : 
      89              : ! *** Public subroutines ***
      90              : 
      91              :    PUBLIC :: becke_constraint, hirshfeld_constraint
      92              : 
      93              : CONTAINS
      94              : 
      95              : ! **************************************************************************************************
      96              : !> \brief Driver routine for calculating a Becke constraint
      97              : !> \param qs_env the qs_env where to build the constraint
      98              : !> \param calc_pot if the potential needs to be recalculated or just integrated
      99              : !> \param calculate_forces logical if potential has to be calculated or only_energy
     100              : !> \par   History
     101              : !>        Created 01.2007 [fschiff]
     102              : !>        Extended functionality 12/15-12/16 [Nico Holmberg]
     103              : ! **************************************************************************************************
     104         2948 :    SUBROUTINE becke_constraint(qs_env, calc_pot, calculate_forces)
     105              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     106              :       LOGICAL                                            :: calc_pot, calculate_forces
     107              : 
     108              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'becke_constraint'
     109              : 
     110              :       INTEGER                                            :: handle
     111              :       TYPE(cdft_control_type), POINTER                   :: cdft_control
     112              :       TYPE(dft_control_type), POINTER                    :: dft_control
     113              : 
     114         2948 :       CALL timeset(routineN, handle)
     115         2948 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     116         2948 :       cdft_control => dft_control%qs_control%cdft_control
     117         2948 :       IF (dft_control%qs_control%cdft .AND. cdft_control%type == outer_scf_becke_constraint) THEN
     118         2948 :          IF (calc_pot) THEN
     119              :             ! Initialize the Becke constraint environment
     120          190 :             CALL becke_constraint_init(qs_env)
     121              :             ! Calculate the Becke weight function and possibly the gradients
     122          190 :             CALL becke_constraint_low(qs_env)
     123              :          END IF
     124              :          ! Integrate the Becke constraint
     125         2948 :          CALL cdft_constraint_integrate(qs_env)
     126              :          ! Calculate forces
     127         2948 :          IF (calculate_forces) CALL cdft_constraint_force(qs_env)
     128              :       END IF
     129         2948 :       CALL timestop(handle)
     130              : 
     131         2948 :    END SUBROUTINE becke_constraint
     132              : 
     133              : ! **************************************************************************************************
     134              : !> \brief Low level routine to build a Becke weight function and its gradients
     135              : !> \param qs_env the qs_env where to build the constraint
     136              : !> \param just_gradients optional logical which determines if only the gradients should be calculated
     137              : !> \par   History
     138              : !>        Created 03.2017 [Nico Holmberg]
     139              : ! **************************************************************************************************
     140          200 :    SUBROUTINE becke_constraint_low(qs_env, just_gradients)
     141              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     142              :       LOGICAL, OPTIONAL                                  :: just_gradients
     143              : 
     144              :       CHARACTER(len=*), PARAMETER :: routineN = 'becke_constraint_low'
     145              : 
     146              :       INTEGER                                            :: handle, i, iatom, igroup, ind(3), ip, j, &
     147              :                                                             jatom, jp, k, natom, np(3), nskipped
     148          200 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: catom
     149              :       INTEGER, DIMENSION(2, 3)                           :: bo, bo_conf
     150              :       LOGICAL                                            :: in_memory, my_just_gradients
     151          200 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: is_constraint, skip_me
     152          200 :       LOGICAL, ALLOCATABLE, DIMENSION(:, :)              :: atom_in_group
     153              :       REAL(kind=dp)                                      :: dist1, dist2, dmyexp, dvol, eps_cavity, &
     154              :                                                             my1, my1_homo, myexp, sum_cell_f_all, &
     155              :                                                             th, tmp_const
     156          200 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: cell_functions, ds_dR_i, ds_dR_j, &
     157          200 :                                                             sum_cell_f_group
     158          200 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: d_sum_Pm_dR, dP_i_dRi
     159          200 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: dP_i_dRj
     160              :       REAL(kind=dp), DIMENSION(3)                        :: cell_v, dist_vec, dmy_dR_i, dmy_dR_j, &
     161              :                                                             dr, dr1_r2, dr_i_dR, dr_ij_dR, &
     162              :                                                             dr_j_dR, grid_p, r, r1, shift
     163          200 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: cutoffs
     164              :       TYPE(becke_constraint_type), POINTER               :: becke_control
     165              :       TYPE(cdft_control_type), POINTER                   :: cdft_control
     166          200 :       TYPE(cdft_group_type), DIMENSION(:), POINTER       :: group
     167              :       TYPE(cell_type), POINTER                           :: cell
     168              :       TYPE(dft_control_type), POINTER                    :: dft_control
     169          200 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     170          200 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: charge
     171              : 
     172          200 :       NULLIFY (cutoffs, cell, dft_control, particle_set, group, charge, cdft_control)
     173          200 :       CALL timeset(routineN, handle)
     174              :       ! Get simulation environment
     175              :       CALL get_qs_env(qs_env, &
     176              :                       cell=cell, &
     177              :                       particle_set=particle_set, &
     178              :                       natom=natom, &
     179          200 :                       dft_control=dft_control)
     180          200 :       cdft_control => dft_control%qs_control%cdft_control
     181          200 :       becke_control => cdft_control%becke_control
     182          200 :       group => cdft_control%group
     183          200 :       cutoffs => becke_control%cutoffs
     184          200 :       IF (cdft_control%atomic_charges) &
     185          106 :          charge => cdft_control%charge
     186          200 :       in_memory = .FALSE.
     187          200 :       IF (cdft_control%save_pot) THEN
     188           82 :          in_memory = becke_control%in_memory
     189              :       END IF
     190          200 :       eps_cavity = becke_control%eps_cavity
     191              :       ! Decide if only gradients need to be calculated
     192          200 :       my_just_gradients = .FALSE.
     193          200 :       IF (PRESENT(just_gradients)) my_just_gradients = just_gradients
     194           10 :       IF (my_just_gradients) THEN
     195           10 :          in_memory = .TRUE.
     196              :          !  Pairwise distances need to be recalculated
     197           10 :          IF (becke_control%vector_buffer%store_vectors) THEN
     198           30 :             ALLOCATE (becke_control%vector_buffer%distances(natom))
     199           30 :             ALLOCATE (becke_control%vector_buffer%distance_vecs(3, natom))
     200           40 :             IF (in_memory) ALLOCATE (becke_control%vector_buffer%pair_dist_vecs(3, natom, natom))
     201           20 :             ALLOCATE (becke_control%vector_buffer%position_vecs(3, natom))
     202              :          END IF
     203           40 :          ALLOCATE (becke_control%vector_buffer%R12(natom, natom))
     204           40 :          DO i = 1, 3
     205           40 :             cell_v(i) = cell%hmat(i, i)
     206              :          END DO
     207           20 :          DO iatom = 1, natom - 1
     208           30 :             DO jatom = iatom + 1, natom
     209           40 :                r = particle_set(iatom)%r
     210           40 :                r1 = particle_set(jatom)%r
     211           40 :                DO i = 1, 3
     212           30 :                   r(i) = MODULO(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
     213           40 :                   r1(i) = MODULO(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
     214              :                END DO
     215           40 :                dist_vec = (r - r1) - ANINT((r - r1)/cell_v)*cell_v
     216           10 :                IF (becke_control%vector_buffer%store_vectors) THEN
     217           40 :                   becke_control%vector_buffer%position_vecs(:, iatom) = r(:)
     218           40 :                   IF (iatom == 1 .AND. jatom == natom) becke_control%vector_buffer%position_vecs(:, jatom) = r1(:)
     219              :                   IF (in_memory) THEN
     220           40 :                      becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom) = dist_vec(:)
     221           40 :                      becke_control%vector_buffer%pair_dist_vecs(:, jatom, iatom) = -dist_vec(:)
     222              :                   END IF
     223              :                END IF
     224           40 :                becke_control%vector_buffer%R12(iatom, jatom) = SQRT(DOT_PRODUCT(dist_vec, dist_vec))
     225           20 :                becke_control%vector_buffer%R12(jatom, iatom) = becke_control%vector_buffer%R12(iatom, jatom)
     226              :             END DO
     227              :          END DO
     228              :       END IF
     229          600 :       ALLOCATE (catom(cdft_control%natoms))
     230              :       IF (cdft_control%save_pot .OR. &
     231          200 :           becke_control%cavity_confine .OR. &
     232              :           becke_control%should_skip) THEN
     233          576 :          ALLOCATE (is_constraint(natom))
     234          616 :          is_constraint = .FALSE.
     235              :       END IF
     236              :       ! This boolean is needed to prevent calculation of atom pairs ji when the pair ij has
     237              :       ! already been calculated (data for pair ji is set using symmetry)
     238              :       ! With gradient precomputation, symmetry exploited for both weight function and gradients
     239          600 :       ALLOCATE (skip_me(natom))
     240          596 :       DO i = 1, cdft_control%natoms
     241          396 :          catom(i) = cdft_control%atoms(i)
     242              :          ! Notice that here is_constraint=.TRUE. also for dummy atoms to properly compute their Becke charges
     243              :          ! A subsequent check (atom_in_group) ensures that the gradients of these dummy atoms are correct
     244              :          IF (cdft_control%save_pot .OR. &
     245          396 :              becke_control%cavity_confine .OR. &
     246              :              becke_control%should_skip) &
     247          578 :             is_constraint(catom(i)) = .TRUE.
     248              :       END DO
     249         2000 :       bo = group(1)%weight%pw_grid%bounds_local
     250          800 :       dvol = group(1)%weight%pw_grid%dvol
     251          800 :       dr = group(1)%weight%pw_grid%dr
     252          800 :       np = group(1)%weight%pw_grid%npts
     253          800 :       shift = -REAL(MODULO(np, 2), dp)*dr/2.0_dp
     254          800 :       DO i = 1, 3
     255          800 :          cell_v(i) = cell%hmat(i, i)
     256              :       END DO
     257              :       ! If requested, allocate storage for gradients
     258          200 :       IF (in_memory) THEN
     259           72 :          bo_conf = bo
     260              :          ! With confinement active, we dont need to store gradients outside
     261              :          ! the confinement bounds since they vanish for all particles
     262           72 :          IF (becke_control%cavity_confine) THEN
     263           64 :             bo_conf(1, 3) = becke_control%confine_bounds(1)
     264           64 :             bo_conf(2, 3) = becke_control%confine_bounds(2)
     265              :          END IF
     266          288 :          ALLOCATE (atom_in_group(SIZE(group), natom))
     267          392 :          atom_in_group = .FALSE.
     268          160 :          DO igroup = 1, SIZE(group)
     269              :             ALLOCATE (group(igroup)%gradients(3*natom, bo_conf(1, 1):bo_conf(2, 1), &
     270              :                                               bo_conf(1, 2):bo_conf(2, 2), &
     271          528 :                                               bo_conf(1, 3):bo_conf(2, 3)))
     272     23089200 :             group(igroup)%gradients = 0.0_dp
     273          264 :             ALLOCATE (group(igroup)%d_sum_const_dR(3, natom))
     274          792 :             group(igroup)%d_sum_const_dR = 0.0_dp
     275          336 :             DO ip = 1, SIZE(group(igroup)%atoms)
     276          264 :                atom_in_group(igroup, group(igroup)%atoms(ip)) = .TRUE.
     277              :             END DO
     278              :          END DO
     279              :       END IF
     280              :       ! Allocate remaining work
     281          600 :       ALLOCATE (sum_cell_f_group(SIZE(group)))
     282          600 :       ALLOCATE (cell_functions(natom))
     283          200 :       IF (in_memory) THEN
     284           72 :          ALLOCATE (ds_dR_j(3))
     285           72 :          ALLOCATE (ds_dR_i(3))
     286          216 :          ALLOCATE (d_sum_Pm_dR(3, natom))
     287          288 :          ALLOCATE (dP_i_dRj(3, natom, natom))
     288          144 :          ALLOCATE (dP_i_dRi(3, natom))
     289          200 :          th = 1.0e-8_dp
     290              :       END IF
     291              :       ! Build constraint
     292         4208 :       DO k = bo(1, 1), bo(2, 1)
     293       170960 :          DO j = bo(1, 2), bo(2, 2)
     294      7373448 :             DO i = bo(1, 3), bo(2, 3)
     295              :                ! If the grid point is too far from all constraint atoms and cavity confinement is active,
     296              :                ! we can skip this grid point as it does not contribute to the weight or gradients
     297      7202688 :                IF (becke_control%cavity_confine) THEN
     298      6424576 :                   IF (becke_control%cavity%array(k, j, i) < eps_cavity) CYCLE
     299              :                END IF
     300     19695324 :                ind = (/k, j, i/)
     301      4923831 :                grid_p(1) = k*dr(1) + shift(1)
     302      4923831 :                grid_p(2) = j*dr(2) + shift(2)
     303      4923831 :                grid_p(3) = i*dr(3) + shift(3)
     304      4923831 :                nskipped = 0
     305     15531637 :                cell_functions = 1.0_dp
     306     15531637 :                skip_me = .FALSE.
     307     15531637 :                IF (becke_control%vector_buffer%store_vectors) becke_control%vector_buffer%distances = 0.0_dp
     308      4923831 :                IF (in_memory) THEN
     309     15260751 :                   d_sum_Pm_dR = 0.0_dp
     310      3790808 :                   DO igroup = 1, SIZE(group)
     311     20552160 :                      group(igroup)%d_sum_const_dR = 0.0_dp
     312              :                   END DO
     313     15260751 :                   dP_i_dRi = 0.0_dp
     314              :                END IF
     315              :                ! Iterate over all atoms in the system
     316     12863480 :                DO iatom = 1, natom
     317     10232484 :                   IF (skip_me(iatom)) THEN
     318       393721 :                      cell_functions(iatom) = 0.0_dp
     319       393721 :                      IF (becke_control%should_skip) THEN
     320       252029 :                         IF (is_constraint(iatom)) nskipped = nskipped + 1
     321       252029 :                         IF (nskipped == cdft_control%natoms) THEN
     322            0 :                            IF (in_memory) THEN
     323            0 :                               IF (becke_control%cavity_confine) THEN
     324            0 :                                  becke_control%cavity%array(k, j, i) = 0.0_dp
     325              :                               END IF
     326              :                            END IF
     327              :                            EXIT
     328              :                         END IF
     329              :                      END IF
     330              :                      CYCLE
     331              :                   END IF
     332      9838763 :                   IF (becke_control%vector_buffer%store_vectors) THEN
     333      9838763 :                      IF (becke_control%vector_buffer%distances(iatom) .EQ. 0.0_dp) THEN
     334     35989816 :                         r = becke_control%vector_buffer%position_vecs(:, iatom)
     335     35989816 :                         dist_vec = (r - grid_p) - ANINT((r - grid_p)/cell_v)*cell_v
     336     35989816 :                         dist1 = SQRT(DOT_PRODUCT(dist_vec, dist_vec))
     337     35989816 :                         becke_control%vector_buffer%distance_vecs(:, iatom) = dist_vec
     338      8997454 :                         becke_control%vector_buffer%distances(iatom) = dist1
     339              :                      ELSE
     340      3365236 :                         dist_vec = becke_control%vector_buffer%distance_vecs(:, iatom)
     341              :                         dist1 = becke_control%vector_buffer%distances(iatom)
     342              :                      END IF
     343              :                   ELSE
     344            0 :                      r = particle_set(iatom)%r
     345            0 :                      DO ip = 1, 3
     346            0 :                         r(ip) = MODULO(r(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
     347              :                      END DO
     348            0 :                      dist_vec = (r - grid_p) - ANINT((r - grid_p)/cell_v)*cell_v
     349            0 :                      dist1 = SQRT(DOT_PRODUCT(dist_vec, dist_vec))
     350              :                   END IF
     351     12469759 :                   IF (dist1 .LE. cutoffs(iatom)) THEN
     352      2173196 :                      IF (in_memory) THEN
     353       761650 :                         IF (dist1 .LE. th) dist1 = th
     354      3046600 :                         dr_i_dR(:) = dist_vec(:)/dist1
     355              :                      END IF
     356      6980971 :                      DO jatom = 1, natom
     357      6980971 :                         IF (jatom .NE. iatom) THEN
     358              :                            ! Using pairwise symmetry, execute block only for such j<i
     359              :                            ! that have previously not been looped over
     360              :                            ! Note that if skip_me(jatom) = .TRUE., this means that the outer
     361              :                            ! loop over iatom skipped this index when iatom=jatom, but we still
     362              :                            ! need to compute the pair for iatom>jatom
     363      2634579 :                            IF (jatom < iatom) THEN
     364      1286576 :                               IF (.NOT. skip_me(jatom)) CYCLE
     365              :                            END IF
     366      1723481 :                            IF (becke_control%vector_buffer%store_vectors) THEN
     367      1723481 :                               IF (becke_control%vector_buffer%distances(jatom) .EQ. 0.0_dp) THEN
     368      4940120 :                                  r1 = becke_control%vector_buffer%position_vecs(:, jatom)
     369      4940120 :                                  dist_vec = (r1 - grid_p) - ANINT((r1 - grid_p)/cell_v)*cell_v
     370      4940120 :                                  dist2 = SQRT(DOT_PRODUCT(dist_vec, dist_vec))
     371      4940120 :                                  becke_control%vector_buffer%distance_vecs(:, jatom) = dist_vec
     372      1235030 :                                  becke_control%vector_buffer%distances(jatom) = dist2
     373              :                               ELSE
     374      1953804 :                                  dist_vec = becke_control%vector_buffer%distance_vecs(:, jatom)
     375              :                                  dist2 = becke_control%vector_buffer%distances(jatom)
     376              :                               END IF
     377              :                            ELSE
     378            0 :                               r1 = particle_set(jatom)%r
     379            0 :                               DO ip = 1, 3
     380            0 :                                  r1(ip) = MODULO(r1(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
     381              :                               END DO
     382            0 :                               dist_vec = (r1 - grid_p) - ANINT((r1 - grid_p)/cell_v)*cell_v
     383            0 :                               dist2 = SQRT(DOT_PRODUCT(dist_vec, dist_vec))
     384              :                            END IF
     385      1723481 :                            IF (in_memory) THEN
     386       484606 :                               IF (becke_control%vector_buffer%store_vectors) THEN
     387      1938424 :                                  dr1_r2 = becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom)
     388              :                               ELSE
     389            0 :                                  dr1_r2 = (r - r1) - ANINT((r - r1)/cell_v)*cell_v
     390              :                               END IF
     391       484606 :                               IF (dist2 .LE. th) dist2 = th
     392       484606 :                               tmp_const = (becke_control%vector_buffer%R12(iatom, jatom)**3)
     393      1938424 :                               dr_ij_dR(:) = dr1_r2(:)/tmp_const
     394              :                               !derivative w.r.t. Rj
     395      1938424 :                               dr_j_dR = dist_vec(:)/dist2
     396      1938424 :                              dmy_dR_j(:) = -(dr_j_dR(:)/becke_control%vector_buffer%R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dR(:))
     397              :                               !derivative w.r.t. Ri
     398      1938424 :                               dmy_dR_i(:) = dr_i_dR(:)/becke_control%vector_buffer%R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dR(:)
     399              :                            END IF
     400              :                            ! myij
     401      1723481 :                            my1 = (dist1 - dist2)/becke_control%vector_buffer%R12(iatom, jatom)
     402      1723481 :                            IF (becke_control%adjust) THEN
     403      1111478 :                               my1_homo = my1 ! Homonuclear quantity needed for gradient
     404      1111478 :                               my1 = my1 + becke_control%aij(iatom, jatom)*(1.0_dp - my1**2)
     405              :                            END IF
     406              :                            ! f(myij)
     407      1723481 :                            myexp = 1.5_dp*my1 - 0.5_dp*my1**3
     408      1723481 :                            IF (in_memory) THEN
     409       484606 :                               dmyexp = 1.5_dp - 1.5_dp*my1**2
     410              :                               tmp_const = (1.5_dp**2)*dmyexp*(1 - myexp**2)* &
     411       484606 :                                           (1.0_dp - ((1.5_dp*myexp - 0.5_dp*(myexp**3))**2))
     412              :                               ! d s(myij)/d R_i
     413      1938424 :                               ds_dR_i(:) = -0.5_dp*tmp_const*dmy_dR_i(:)
     414              :                               ! d s(myij)/d R_j
     415      1938424 :                               ds_dR_j(:) = -0.5_dp*tmp_const*dmy_dR_j(:)
     416       484606 :                               IF (becke_control%adjust) THEN
     417              :                                  tmp_const = 1.0_dp - 2.0_dp*my1_homo* &
     418       268771 :                                              becke_control%aij(iatom, jatom)
     419      1075084 :                                  ds_dR_i(:) = ds_dR_i(:)*tmp_const
     420              :                                  ! tmp_const is same for both since aij=-aji and myij=-myji
     421      1075084 :                                  ds_dR_j(:) = ds_dR_j(:)*tmp_const
     422              :                               END IF
     423              :                            END IF
     424              :                            ! s(myij) = f[f(f{myij})]
     425      1723481 :                            myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
     426      1723481 :                            myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
     427      1723481 :                            tmp_const = 0.5_dp*(1.0_dp - myexp)
     428      1723481 :                            cell_functions(iatom) = cell_functions(iatom)*tmp_const
     429      1723481 :                            IF (in_memory) THEN
     430       484606 :                               IF (ABS(tmp_const) .LE. th) tmp_const = tmp_const + th
     431              :                               ! P_i independent part of dP_i/dR_i
     432      1938424 :                               dP_i_dRi(:, iatom) = dP_i_dRi(:, iatom) + ds_dR_i(:)/tmp_const
     433              :                               ! P_i independent part of dP_i/dR_j
     434      1938424 :                               dP_i_dRj(:, iatom, jatom) = ds_dR_j(:)/tmp_const
     435              :                            END IF
     436              : 
     437      1723481 :                            IF (dist2 .LE. cutoffs(jatom)) THEN
     438       911098 :                               tmp_const = 0.5_dp*(1.0_dp + myexp) ! s(myji)
     439       911098 :                               cell_functions(jatom) = cell_functions(jatom)*tmp_const
     440       911098 :                               IF (in_memory) THEN
     441       277044 :                                  IF (ABS(tmp_const) .LE. th) tmp_const = tmp_const + th
     442              :                                  ! P_j independent part of dP_j/dR_i
     443              :                                  ! d s(myji)/d R_i = -d s(myij)/d R_i
     444      1108176 :                                  dP_i_dRj(:, jatom, iatom) = -ds_dR_i(:)/tmp_const
     445              :                                  ! P_j independent part of dP_j/dR_j
     446              :                                  ! d s(myji)/d R_j = -d s(myij)/d R_j
     447      1108176 :                                  dP_i_dRi(:, jatom) = dP_i_dRi(:, jatom) - ds_dR_j(:)/tmp_const
     448              :                               END IF
     449              :                            ELSE
     450       812383 :                               skip_me(jatom) = .TRUE.
     451              :                            END IF
     452              :                         END IF
     453              :                      END DO ! jatom
     454      2173196 :                      IF (in_memory) THEN
     455              :                         ! Final value of dP_i_dRi
     456      3046600 :                         dP_i_dRi(:, iatom) = cell_functions(iatom)*dP_i_dRi(:, iatom)
     457              :                         ! Update relevant sums with value
     458      3046600 :                         d_sum_Pm_dR(:, iatom) = d_sum_Pm_dR(:, iatom) + dP_i_dRi(:, iatom)
     459       761650 :                         IF (is_constraint(iatom)) THEN
     460      1682312 :                            DO igroup = 1, SIZE(group)
     461       920662 :                               IF (.NOT. atom_in_group(igroup, iatom)) CYCLE
     462      1380999 :                               DO jp = 1, SIZE(group(igroup)%atoms)
     463      1380999 :                                  IF (iatom == group(igroup)%atoms(jp)) THEN
     464              :                                     ip = jp
     465              :                                     EXIT
     466              :                                  END IF
     467              :                               END DO
     468              :                               group(igroup)%d_sum_const_dR(1:3, iatom) = group(igroup)%d_sum_const_dR(1:3, iatom) + &
     469      4444298 :                                                                          group(igroup)%coeff(ip)*dP_i_dRi(:, iatom)
     470              :                            END DO
     471              :                         END IF
     472      2284950 :                         DO jatom = 1, natom
     473      2284950 :                            IF (jatom .NE. iatom) THEN
     474              :                               ! Final value of dP_i_dRj
     475      3046600 :                               dP_i_dRj(:, iatom, jatom) = cell_functions(iatom)*dP_i_dRj(:, iatom, jatom)
     476              :                               ! Update where needed
     477      3046600 :                               d_sum_Pm_dR(:, jatom) = d_sum_Pm_dR(:, jatom) + dP_i_dRj(:, iatom, jatom)
     478       761650 :                               IF (is_constraint(iatom)) THEN
     479      1682312 :                                  DO igroup = 1, SIZE(group)
     480       920662 :                                     IF (.NOT. atom_in_group(igroup, iatom)) CYCLE
     481       920662 :                                     ip = -1
     482      1380999 :                                     DO jp = 1, SIZE(group(igroup)%atoms)
     483      1380999 :                                        IF (iatom == group(igroup)%atoms(jp)) THEN
     484              :                                           ip = jp
     485              :                                           EXIT
     486              :                                        END IF
     487              :                                     END DO
     488              :                                     group(igroup)%d_sum_const_dR(1:3, jatom) = group(igroup)%d_sum_const_dR(1:3, jatom) + &
     489              :                                                                                group(igroup)%coeff(ip)* &
     490      4444298 :                                                                                dP_i_dRj(:, iatom, jatom)
     491              :                                  END DO
     492              :                               END IF
     493              :                            END IF
     494              :                         END DO
     495              :                      END IF
     496              :                   ELSE
     497      7665567 :                      cell_functions(iatom) = 0.0_dp
     498      7665567 :                      skip_me(iatom) = .TRUE.
     499      7665567 :                      IF (becke_control%should_skip) THEN
     500      4629324 :                         IF (is_constraint(iatom)) nskipped = nskipped + 1
     501      4629324 :                         IF (nskipped == cdft_control%natoms) THEN
     502      2292835 :                            IF (in_memory) THEN
     503       897142 :                               IF (becke_control%cavity_confine) THEN
     504       897142 :                                  becke_control%cavity%array(k, j, i) = 0.0_dp
     505              :                               END IF
     506              :                            END IF
     507              :                            EXIT
     508              :                         END IF
     509              :                      END IF
     510              :                   END IF
     511              :                END DO !iatom
     512      4923831 :                IF (nskipped == cdft_control%natoms) CYCLE
     513              :                ! Sum up cell functions
     514      5442400 :                sum_cell_f_group = 0.0_dp
     515      5442400 :                DO igroup = 1, SIZE(group)
     516     11107611 :                   DO ip = 1, SIZE(group(igroup)%atoms)
     517              :                      sum_cell_f_group(igroup) = sum_cell_f_group(igroup) + group(igroup)%coeff(ip)* &
     518      8476615 :                                                 cell_functions(group(igroup)%atoms(ip))
     519              :                   END DO
     520              :                END DO
     521      2630996 :                sum_cell_f_all = 0.0_dp
     522      8435725 :                DO ip = 1, natom
     523      8435725 :                   sum_cell_f_all = sum_cell_f_all + cell_functions(ip)
     524              :                END DO
     525              :                ! Gradients at (k,j,i)
     526      2630996 :                IF (in_memory .AND. ABS(sum_cell_f_all) .GT. 0.0_dp) THEN
     527      1072432 :                   DO igroup = 1, SIZE(group)
     528      2248084 :                      DO iatom = 1, natom
     529              :                         group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) = &
     530              :                            group(igroup)%d_sum_const_dR(1:3, iatom)/sum_cell_f_all - sum_cell_f_group(igroup)* &
     531      5290434 :                            d_sum_Pm_dR(1:3, iatom)/(sum_cell_f_all**2)
     532              :                      END DO
     533              :                   END DO
     534              :                END IF
     535              :                ! Weight function(s) at (k,j,i)
     536      2797748 :                IF (.NOT. my_just_gradients .AND. ABS(sum_cell_f_all) .GT. 0.000001) THEN
     537      2744726 :                   DO igroup = 1, SIZE(group)
     538      2744726 :                      group(igroup)%weight%array(k, j, i) = sum_cell_f_group(igroup)/sum_cell_f_all
     539              :                   END DO
     540      1282159 :                   IF (cdft_control%atomic_charges) THEN
     541      2164389 :                      DO iatom = 1, cdft_control%natoms
     542      2164389 :                         charge(iatom)%array(k, j, i) = cell_functions(catom(iatom))/sum_cell_f_all
     543              :                      END DO
     544              :                   END IF
     545              :                END IF
     546              :             END DO
     547              :          END DO
     548              :       END DO
     549              :       ! Release storage
     550          200 :       IF (in_memory) THEN
     551           72 :          DEALLOCATE (ds_dR_j)
     552           72 :          DEALLOCATE (ds_dR_i)
     553           72 :          DEALLOCATE (d_sum_Pm_dR)
     554           72 :          DEALLOCATE (dP_i_dRj)
     555           72 :          DEALLOCATE (dP_i_dRi)
     556          160 :          DO igroup = 1, SIZE(group)
     557          160 :             DEALLOCATE (group(igroup)%d_sum_const_dR)
     558              :          END DO
     559           72 :          DEALLOCATE (atom_in_group)
     560           72 :          IF (becke_control%vector_buffer%store_vectors) THEN
     561           72 :             DEALLOCATE (becke_control%vector_buffer%pair_dist_vecs)
     562              :          END IF
     563              :       END IF
     564          200 :       NULLIFY (cutoffs)
     565          200 :       IF (ALLOCATED(is_constraint)) &
     566          192 :          DEALLOCATE (is_constraint)
     567          200 :       DEALLOCATE (catom)
     568          200 :       DEALLOCATE (cell_functions)
     569          200 :       DEALLOCATE (skip_me)
     570          200 :       DEALLOCATE (sum_cell_f_group)
     571          200 :       DEALLOCATE (becke_control%vector_buffer%R12)
     572          200 :       IF (becke_control%vector_buffer%store_vectors) THEN
     573          200 :          DEALLOCATE (becke_control%vector_buffer%distances)
     574          200 :          DEALLOCATE (becke_control%vector_buffer%distance_vecs)
     575          200 :          DEALLOCATE (becke_control%vector_buffer%position_vecs)
     576              :       END IF
     577          200 :       CALL timestop(handle)
     578              : 
     579          400 :    END SUBROUTINE becke_constraint_low
     580              : 
     581              : ! **************************************************************************************************
     582              : !> \brief Driver routine for calculating a Hirshfeld constraint
     583              : !> \param qs_env ...
     584              : !> \param calc_pot ...
     585              : !> \param calculate_forces ...
     586              : ! **************************************************************************************************
     587           86 :    SUBROUTINE hirshfeld_constraint(qs_env, calc_pot, calculate_forces)
     588              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     589              :       LOGICAL                                            :: calc_pot, calculate_forces
     590              : 
     591              :       CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint'
     592              : 
     593              :       INTEGER                                            :: handle
     594              :       TYPE(cdft_control_type), POINTER                   :: cdft_control
     595              :       TYPE(dft_control_type), POINTER                    :: dft_control
     596              : 
     597           86 :       CALL timeset(routineN, handle)
     598           86 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     599           86 :       cdft_control => dft_control%qs_control%cdft_control
     600           86 :       IF (dft_control%qs_control%cdft .AND. cdft_control%type == outer_scf_hirshfeld_constraint) THEN
     601           86 :          IF (calc_pot) THEN
     602              :             ! Initialize the Hirshfeld constraint environment
     603           22 :             CALL hirshfeld_constraint_init(qs_env)
     604              :             ! Calculate the Hirshfeld weight function and possibly the gradients
     605           22 :             CALL hirshfeld_constraint_low(qs_env)
     606              :          END IF
     607              :          ! Integrate the Hirshfeld constraint
     608           86 :          CALL cdft_constraint_integrate(qs_env)
     609              :          ! Calculate forces
     610           86 :          IF (calculate_forces) CALL cdft_constraint_force(qs_env)
     611              :       END IF
     612           86 :       CALL timestop(handle)
     613              : 
     614           86 :    END SUBROUTINE hirshfeld_constraint
     615              : 
     616              : ! **************************************************************************************************
     617              : !> \brief Calculates Hirshfeld constraints
     618              : !> \param qs_env ...
     619              : !> \param just_gradients ...
     620              : ! **************************************************************************************************
     621           24 :    SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients)
     622              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     623              :       LOGICAL, OPTIONAL                                  :: just_gradients
     624              : 
     625              :       CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint_low'
     626              : 
     627              :       INTEGER :: atom_a, atoms_memory, atoms_memory_num, handle, i, iatom, iex, igroup, ikind, &
     628              :          ithread, j, k, natom, npme, nthread, num_atoms, num_species, numexp, subpatch_pattern
     629           24 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: num_species_small
     630              :       INTEGER, DIMENSION(2, 3)                           :: bo
     631              :       INTEGER, DIMENSION(3)                              :: lb_pw, lb_rs, npts, ub_pw, ub_rs
     632           24 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores
     633              :       LOGICAL                                            :: my_just_gradients
     634           24 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: compute_charge, is_constraint
     635              :       REAL(kind=dp)                                      :: alpha, coef, eps_rho_rspace, exp_eval, &
     636              :                                                             prefactor, radius
     637           24 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: coefficients
     638              :       REAL(kind=dp), DIMENSION(3)                        :: dr_pw, dr_rs, origin, r2, r_pbc, ra
     639           24 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab
     640           24 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     641              :       TYPE(cdft_control_type), POINTER                   :: cdft_control
     642              :       TYPE(cell_type), POINTER                           :: cell
     643              :       TYPE(dft_control_type), POINTER                    :: dft_control
     644              :       TYPE(hirshfeld_constraint_type), POINTER           :: hirshfeld_control
     645              :       TYPE(hirshfeld_type), POINTER                      :: hirshfeld_env
     646              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     647           24 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     648              :       TYPE(pw_env_type), POINTER                         :: pw_env
     649              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     650           24 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:)    :: pw_single_dr
     651           24 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
     652              :       TYPE(qs_rho_type), POINTER                         :: rho
     653              :       TYPE(realspace_grid_desc_type), POINTER            :: auxbas_rs_desc
     654          936 :       TYPE(realspace_grid_type)                          :: rs_rho_all, rs_rho_constr
     655              :       TYPE(realspace_grid_type), ALLOCATABLE, &
     656           24 :          DIMENSION(:)                                    :: rs_single, rs_single_charge, rs_single_dr
     657              : 
     658           24 :       NULLIFY (atom_list, atomic_kind_set, dft_control, &
     659           24 :                hirshfeld_env, particle_set, pw_env, auxbas_pw_pool, para_env, &
     660           24 :                auxbas_rs_desc, cdft_control, pab, &
     661           24 :                hirshfeld_control, cell, rho_r, rho)
     662              : 
     663           24 :       CALL timeset(routineN, handle)
     664              :       CALL get_qs_env(qs_env, &
     665              :                       atomic_kind_set=atomic_kind_set, &
     666              :                       particle_set=particle_set, &
     667              :                       natom=natom, &
     668              :                       cell=cell, &
     669              :                       rho=rho, &
     670              :                       dft_control=dft_control, &
     671              :                       para_env=para_env, &
     672           24 :                       pw_env=pw_env)
     673           24 :       CALL qs_rho_get(rho, rho_r=rho_r)
     674              : 
     675           24 :       num_atoms = natom
     676              : 
     677           24 :       cdft_control => dft_control%qs_control%cdft_control
     678           24 :       hirshfeld_control => cdft_control%hirshfeld_control
     679           24 :       hirshfeld_env => hirshfeld_control%hirshfeld_env
     680              : 
     681              :       ! Check if only gradient should be calculated, if gradients should be precomputed
     682           24 :       my_just_gradients = .FALSE.
     683           24 :       IF (PRESENT(just_gradients)) my_just_gradients = just_gradients
     684            2 :       IF (my_just_gradients) THEN
     685            2 :          cdft_control%in_memory = .TRUE.
     686            2 :          hirshfeld_control%print_density = .FALSE.
     687              :       END IF
     688              : 
     689           72 :       ALLOCATE (coefficients(natom))
     690           72 :       ALLOCATE (is_constraint(natom))
     691              : 
     692           24 :       subpatch_pattern = 0
     693           24 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     694           24 :       radius = 100.0_dp
     695              : 
     696           24 :       dr_pw(1) = rho_r(1)%pw_grid%dr(1)
     697           24 :       dr_pw(2) = rho_r(1)%pw_grid%dr(2)
     698           24 :       dr_pw(3) = rho_r(1)%pw_grid%dr(3)
     699           96 :       lb_pw(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
     700           96 :       ub_pw(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
     701           96 :       npts = rho_r(1)%pw_grid%npts
     702           24 :       origin(1) = (dr_pw(1)*npts(1))*0.5_dp
     703           24 :       origin(2) = (dr_pw(2)*npts(2))*0.5_dp
     704           24 :       origin(3) = (dr_pw(3)*npts(3))*0.5_dp
     705              : 
     706              :       CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
     707           24 :                       auxbas_pw_pool=auxbas_pw_pool)
     708           24 :       CALL rs_grid_create(rs_rho_all, auxbas_rs_desc)
     709           24 :       CALL rs_grid_zero(rs_rho_all)
     710              : 
     711           24 :       dr_rs(1) = rs_rho_all%desc%dh(1, 1)
     712           24 :       dr_rs(2) = rs_rho_all%desc%dh(2, 2)
     713           24 :       dr_rs(3) = rs_rho_all%desc%dh(3, 3)
     714           24 :       lb_rs(1) = LBOUND(rs_rho_all%r(:, :, :), 1)
     715           24 :       lb_rs(2) = LBOUND(rs_rho_all%r(:, :, :), 2)
     716           24 :       lb_rs(3) = LBOUND(rs_rho_all%r(:, :, :), 3)
     717           24 :       ub_rs(1) = UBOUND(rs_rho_all%r(:, :, :), 1)
     718           24 :       ub_rs(2) = UBOUND(rs_rho_all%r(:, :, :), 2)
     719           24 :       ub_rs(3) = UBOUND(rs_rho_all%r(:, :, :), 3)
     720              : 
     721              :       ! For each CDFT group
     722           48 :       DO igroup = 1, SIZE(cdft_control%group)
     723              : 
     724           24 :          IF (igroup == 2 .AND. .NOT. cdft_control%in_memory) THEN
     725            0 :             CALL rs_grid_zero(rs_rho_all)
     726              :          END IF
     727          240 :          bo = cdft_control%group(igroup)%weight%pw_grid%bounds_local
     728              : 
     729              :          ! Coefficients
     730           72 :          coefficients(:) = 0.0_dp
     731           72 :          is_constraint = .FALSE.
     732           70 :          DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
     733           46 :             coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
     734           70 :             is_constraint(cdft_control%group(igroup)%atoms(i)) = .TRUE.
     735              :          END DO
     736              : 
     737              :          ! rs_rho_constr: Sum of isolated Gaussian densities over constraint atoms in this constraint group
     738           24 :          CALL rs_grid_create(rs_rho_constr, auxbas_rs_desc)
     739           24 :          CALL rs_grid_zero(rs_rho_constr)
     740              : 
     741              :          ! rs_single: Gaussian density over single atoms when required
     742           24 :          IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
     743            0 :             ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic(cdft_control%natoms))
     744            0 :             ALLOCATE (rs_single(cdft_control%natoms))
     745            0 :             DO i = 1, cdft_control%natoms
     746            0 :                CALL rs_grid_create(rs_single(i), auxbas_rs_desc)
     747            0 :                CALL rs_grid_zero(rs_single(i))
     748              :             END DO
     749              :          END IF
     750              : 
     751              :          ! Setup pw
     752           24 :          CALL pw_zero(cdft_control%group(igroup)%weight)
     753              : 
     754           24 :          CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
     755           24 :          CALL pw_set(cdft_control%group(igroup)%hw_rho_total_constraint, 1.0_dp)
     756              : 
     757           24 :          IF (igroup == 1) THEN
     758           24 :             CALL auxbas_pw_pool%create_pw(cdft_control%hw_rho_total)
     759           24 :             CALL pw_set(cdft_control%hw_rho_total, 1.0_dp)
     760              : 
     761           24 :             IF (hirshfeld_control%print_density) THEN
     762            0 :                DO iatom = 1, cdft_control%natoms
     763            0 :                   CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic(iatom))
     764            0 :                   CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic(iatom), 1.0_dp)
     765              :                END DO
     766              :             END IF
     767              :          END IF
     768              : 
     769           24 :          IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
     770           40 :             ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge(cdft_control%natoms))
     771          192 :             ALLOCATE (rs_single_charge(cdft_control%natoms))
     772           24 :             ALLOCATE (compute_charge(natom))
     773           24 :             compute_charge = .FALSE.
     774              : 
     775           24 :             DO i = 1, cdft_control%natoms
     776           16 :                CALL rs_grid_create(rs_single_charge(i), auxbas_rs_desc)
     777           16 :                CALL rs_grid_zero(rs_single_charge(i))
     778           24 :                compute_charge(cdft_control%atoms(i)) = .TRUE.
     779              :             END DO
     780              : 
     781           24 :             DO iatom = 1, cdft_control%natoms
     782           16 :                CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom))
     783           24 :                CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom), 1.0_dp)
     784              :             END DO
     785              :          END IF
     786              : 
     787           24 :          ALLOCATE (pab(1, 1))
     788           24 :          nthread = 1
     789           24 :          ithread = 0
     790              : 
     791           72 :          DO ikind = 1, SIZE(atomic_kind_set)
     792           48 :             numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
     793           48 :             IF (numexp <= 0) CYCLE
     794           48 :             CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
     795          144 :             ALLOCATE (cores(num_species))
     796              : 
     797          180 :             DO iex = 1, numexp
     798          132 :                alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
     799          132 :                coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
     800          132 :                npme = 0
     801          264 :                cores = 0
     802          264 :                DO iatom = 1, num_species
     803          132 :                   atom_a = atom_list(iatom)
     804          132 :                   ra(:) = pbc(particle_set(atom_a)%r, cell)
     805          264 :                   IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
     806          128 :                      IF (MODULO(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
     807           64 :                         npme = npme + 1
     808           64 :                         cores(npme) = iatom
     809              :                      END IF
     810              :                   ELSE
     811            4 :                      npme = npme + 1
     812            4 :                      cores(npme) = iatom
     813              :                   END IF
     814              :                END DO
     815          248 :                DO j = 1, npme
     816           68 :                   iatom = cores(j)
     817           68 :                   atom_a = atom_list(iatom)
     818           68 :                   pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
     819           68 :                   ra(:) = pbc(particle_set(atom_a)%r, cell)
     820              : 
     821           68 :                   IF (hirshfeld_control%use_atomic_cutoff) THEN
     822              :                      radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
     823              :                                                        ra=ra, rb=ra, rp=ra, &
     824              :                                                        zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
     825              :                                                        pab=pab, o1=0, o2=0, &  ! without map_consistent
     826           68 :                                                        prefactor=1.0_dp, cutoff=0.0_dp)
     827              :                   END IF
     828              : 
     829           68 :                   IF (igroup == 1) THEN
     830              :                      CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
     831              :                                                 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
     832              :                                                 rs_rho_all, radius=radius, &
     833              :                                                 ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., &
     834           68 :                                                 subpatch_pattern=subpatch_pattern)
     835              :                   END IF
     836              : 
     837           68 :                   IF (is_constraint(atom_a)) THEN
     838              :                      CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
     839              :                                                 (/0.0_dp, 0.0_dp, 0.0_dp/), coefficients(atom_a), &
     840              :                                                 pab, 0, 0, rs_rho_constr, &
     841              :                                                 radius=radius, &
     842              :                                                 ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., &
     843           67 :                                                 subpatch_pattern=subpatch_pattern)
     844              :                   END IF
     845              : 
     846           68 :                   IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
     847            0 :                      IF (is_constraint(atom_a)) THEN
     848            0 :                      DO iatom = 1, cdft_control%natoms
     849            0 :                         IF (atom_a == cdft_control%atoms(iatom)) EXIT
     850              :                      END DO
     851            0 :                      CPASSERT(iatom <= cdft_control%natoms)
     852              :                      CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
     853              :                                                 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
     854              :                                                 rs_single(iatom), radius=radius, &
     855              :                                                 ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., &
     856            0 :                                                 subpatch_pattern=subpatch_pattern)
     857              :                      END IF
     858              :                   END IF
     859              : 
     860          200 :                   IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
     861           22 :                      IF (compute_charge(atom_a)) THEN
     862           33 :                         DO iatom = 1, cdft_control%natoms
     863           33 :                            IF (atom_a == cdft_control%atoms(iatom)) EXIT
     864              :                         END DO
     865           22 :                         CPASSERT(iatom <= cdft_control%natoms)
     866              :                         CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
     867              :                                                    (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
     868              :                                                    rs_single_charge(iatom), radius=radius, &
     869              :                                                    ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., &
     870           22 :                                                    subpatch_pattern=subpatch_pattern)
     871              :                      END IF
     872              :                   END IF
     873              : 
     874              :                END DO
     875              :             END DO
     876          120 :             DEALLOCATE (cores)
     877              :          END DO
     878           24 :          DEALLOCATE (pab)
     879              : 
     880           24 :          IF (igroup == 1) THEN
     881           24 :             CALL transfer_rs2pw(rs_rho_all, cdft_control%hw_rho_total)
     882              :          END IF
     883              : 
     884           24 :          CALL transfer_rs2pw(rs_rho_constr, cdft_control%group(igroup)%hw_rho_total_constraint)
     885           24 :          CALL rs_grid_release(rs_rho_constr)
     886              : 
     887              :          ! Calculate weight function
     888              :          CALL hfun_scale(cdft_control%group(igroup)%weight%array, &
     889              :                          cdft_control%group(igroup)%hw_rho_total_constraint%array, &
     890              :                          cdft_control%hw_rho_total%array, divide=.TRUE., &
     891           24 :                          small=hirshfeld_control%eps_cutoff)
     892              : 
     893              :          ! Calculate charges
     894           24 :          IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
     895           24 :             DO i = 1, cdft_control%natoms
     896           16 :                CALL transfer_rs2pw(rs_single_charge(i), cdft_control%group(igroup)%hw_rho_atomic_charge(i))
     897              :                CALL hfun_scale(cdft_control%charge(i)%array, &
     898              :                                cdft_control%group(igroup)%hw_rho_atomic_charge(i)%array, &
     899              :                                cdft_control%hw_rho_total%array, divide=.TRUE., &
     900           24 :                                small=hirshfeld_control%eps_cutoff)
     901              :             END DO
     902              :          END IF
     903              : 
     904              :          ! Print atomic densities if requested
     905           48 :          IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
     906            0 :             DO i = 1, cdft_control%natoms
     907            0 :                CALL transfer_rs2pw(rs_single(i), cdft_control%group(igroup)%hw_rho_atomic(i))
     908              :             END DO
     909            0 :             CALL cdft_print_hirshfeld_density(qs_env)
     910              :          END IF
     911              : 
     912              :       END DO
     913              : 
     914           48 :       DO igroup = 1, SIZE(cdft_control%group)
     915              : 
     916           24 :          CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
     917              : 
     918           24 :          IF (.NOT. cdft_control%in_memory .AND. igroup == 1) THEN
     919           22 :             CALL auxbas_pw_pool%give_back_pw(cdft_control%hw_rho_total)
     920              :          END IF
     921              : 
     922           24 :          IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
     923            0 :             DO i = 1, cdft_control%natoms
     924            0 :                CALL rs_grid_release(rs_single(i))
     925            0 :                CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic(i))
     926              :             END DO
     927            0 :             DEALLOCATE (rs_single)
     928            0 :             DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic)
     929              :          END IF
     930              : 
     931           48 :          IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
     932           24 :             DO i = 1, cdft_control%natoms
     933           16 :                CALL rs_grid_release(rs_single_charge(i))
     934           24 :                CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(i))
     935              :             END DO
     936           24 :             DEALLOCATE (rs_single_charge)
     937            8 :             DEALLOCATE (compute_charge)
     938            8 :             DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge)
     939              :          END IF
     940              : 
     941              :       END DO
     942              : 
     943           24 :       IF (cdft_control%in_memory) THEN
     944            4 :          DO igroup = 1, SIZE(cdft_control%group)
     945              :             ALLOCATE (cdft_control%group(igroup)%gradients_x(1*natom, lb_pw(1):ub_pw(1), &
     946           12 :                                                              lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
     947       195284 :             cdft_control%group(igroup)%gradients_x(:, :, :, :) = 0.0_dp
     948              :          END DO
     949              :       END IF
     950              : 
     951           24 :       IF (cdft_control%in_memory) THEN
     952            4 :          DO igroup = 1, SIZE(cdft_control%group)
     953              : 
     954            2 :             ALLOCATE (pab(1, 1))
     955            2 :             nthread = 1
     956            2 :             ithread = 0
     957            2 :             atoms_memory = hirshfeld_control%atoms_memory
     958              : 
     959            6 :             DO ikind = 1, SIZE(atomic_kind_set)
     960            4 :                numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
     961            4 :                IF (numexp <= 0) CYCLE
     962            4 :                CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
     963              : 
     964           16 :                ALLOCATE (pw_single_dr(num_species))
     965           96 :                ALLOCATE (rs_single_dr(num_species))
     966              : 
     967            8 :                DO i = 1, num_species
     968            4 :                   CALL auxbas_pw_pool%create_pw(pw_single_dr(i))
     969            8 :                   CALL pw_zero(pw_single_dr(i))
     970              :                END DO
     971              : 
     972           16 :                atoms_memory_num = SIZE((/(j, j=1, num_species, atoms_memory)/))
     973              : 
     974              :                ! Can't store all pw grids, therefore split into groups of size atom_memory
     975              :                ! Ideally this code should be re-written to be more memory efficient
     976            4 :                IF (num_species > atoms_memory) THEN
     977            0 :                   ALLOCATE (num_species_small(atoms_memory_num + 1))
     978            0 :                   num_species_small(1:atoms_memory_num) = (/(j, j=1, num_species, atoms_memory)/)
     979            0 :                   num_species_small(atoms_memory_num + 1) = num_species
     980              :                ELSE
     981            4 :                   ALLOCATE (num_species_small(2))
     982           12 :                   num_species_small(:) = (/1, num_species/)
     983              :                END IF
     984              : 
     985            8 :                DO k = 1, SIZE(num_species_small) - 1
     986            4 :                   IF (num_species > atoms_memory) THEN
     987            0 :                      ALLOCATE (cores(num_species_small(k + 1) - (num_species_small(k) - 1)))
     988              :                   ELSE
     989           12 :                      ALLOCATE (cores(num_species))
     990              :                   END IF
     991              : 
     992            8 :                   DO i = num_species_small(k), num_species_small(k + 1)
     993            4 :                      CALL rs_grid_create(rs_single_dr(i), auxbas_rs_desc)
     994            8 :                      CALL rs_grid_zero(rs_single_dr(i))
     995              :                   END DO
     996           36 :                   DO iex = 1, numexp
     997              : 
     998           32 :                      alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
     999           32 :                      coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
    1000           32 :                      prefactor = 2.0_dp*alpha
    1001           32 :                      npme = 0
    1002           64 :                      cores = 0
    1003              : 
    1004           64 :                      DO iatom = 1, SIZE(cores)
    1005           32 :                         atom_a = atom_list(iatom + (num_species_small(k) - 1))
    1006           32 :                         ra(:) = pbc(particle_set(atom_a)%r, cell)
    1007              : 
    1008           64 :                         IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
    1009           32 :                            IF (MODULO(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
    1010           16 :                               npme = npme + 1
    1011           16 :                               cores(npme) = iatom
    1012              :                            END IF
    1013              :                         ELSE
    1014            0 :                            npme = npme + 1
    1015            0 :                            cores(npme) = iatom
    1016              :                         END IF
    1017              :                      END DO
    1018           52 :                      DO j = 1, npme
    1019           16 :                         iatom = cores(j)
    1020           16 :                         atom_a = atom_list(iatom + (num_species_small(k) - 1))
    1021           16 :                         pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
    1022           16 :                         ra(:) = pbc(particle_set(atom_a)%r, cell)
    1023              :                         subpatch_pattern = 0
    1024              : 
    1025              :                         ! Calculate cutoff
    1026           16 :                         IF (hirshfeld_control%use_atomic_cutoff) THEN
    1027              :                            radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
    1028              :                                                              ra=ra, rb=ra, rp=ra, &
    1029              :                                                              zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
    1030              :                                                              pab=pab, o1=0, o2=0, &  ! without map_consistent
    1031           16 :                                                              prefactor=1.0_dp, cutoff=0.0_dp)
    1032              :                         END IF
    1033              : 
    1034              :                         CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
    1035              :                                                    (/0.0_dp, 0.0_dp, 0.0_dp/), prefactor, &
    1036              :                                                    pab, 0, 0, rs_single_dr(iatom + (num_species_small(k) - 1)), &
    1037              :                                                    radius=radius, &
    1038              :                                                    ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., &
    1039           48 :                                                    subpatch_pattern=subpatch_pattern)
    1040              : 
    1041              :                      END DO
    1042              :                   END DO
    1043              : 
    1044            8 :                   DO iatom = num_species_small(k), num_species_small(k + 1)
    1045            4 :                      CALL transfer_rs2pw(rs_single_dr(iatom), pw_single_dr(iatom))
    1046            8 :                      CALL rs_grid_release(rs_single_dr(iatom))
    1047              :                   END DO
    1048              : 
    1049            8 :                   DEALLOCATE (cores)
    1050              :                END DO
    1051              : 
    1052            8 :                DO iatom = 1, num_species
    1053            4 :                   atom_a = atom_list(iatom)
    1054       134564 :                   cdft_control%group(igroup)%gradients_x(atom_a, :, :, :) = pw_single_dr(iatom)%array(:, :, :)
    1055            8 :                   CALL auxbas_pw_pool%give_back_pw(pw_single_dr(iatom))
    1056              :                END DO
    1057              : 
    1058            8 :                DEALLOCATE (rs_single_dr)
    1059            4 :                DEALLOCATE (num_species_small)
    1060           10 :                DEALLOCATE (pw_single_dr)
    1061              :             END DO
    1062            4 :             DEALLOCATE (pab)
    1063              :          END DO
    1064              :       END IF
    1065              : 
    1066           24 :       IF (cdft_control%in_memory) THEN
    1067            4 :          DO igroup = 1, SIZE(cdft_control%group)
    1068              :             ALLOCATE (cdft_control%group(igroup)%gradients_y(1*num_atoms, lb_pw(1):ub_pw(1), &
    1069           12 :                                                              lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
    1070              :             ALLOCATE (cdft_control%group(igroup)%gradients_z(1*num_atoms, lb_pw(1):ub_pw(1), &
    1071           10 :                                                              lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
    1072       195282 :             cdft_control%group(igroup)%gradients_y(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
    1073       195284 :             cdft_control%group(igroup)%gradients_z(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
    1074              :          END DO
    1075              :       END IF
    1076              : 
    1077              :       ! Calculate gradient if requested
    1078           24 :       IF (cdft_control%in_memory) THEN
    1079              : 
    1080            4 :          DO igroup = 1, SIZE(cdft_control%group)
    1081              : 
    1082              :             ! Coefficients
    1083            6 :             coefficients(:) = 0.0_dp
    1084            6 :             is_constraint = .FALSE.
    1085            6 :             DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
    1086            4 :                coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
    1087            6 :                is_constraint(cdft_control%group(igroup)%atoms(i)) = .TRUE.
    1088              :             END DO
    1089              : 
    1090           84 :             DO k = lb_pw(3), ub_pw(3)
    1091         3282 :                DO j = lb_pw(2), ub_pw(2)
    1092        67280 :                   DO i = lb_pw(1), ub_pw(1)
    1093       195200 :                   DO iatom = 1, natom
    1094              : 
    1095       512000 :                      ra(:) = particle_set(iatom)%r
    1096              : 
    1097       192000 :                      IF (cdft_control%hw_rho_total%array(i, j, k) > hirshfeld_control%eps_cutoff) THEN
    1098              : 
    1099              :                         exp_eval = (coefficients(iatom) - &
    1100              :                                     cdft_control%group(igroup)%weight%array(i, j, k))/ &
    1101       128000 :                                    cdft_control%hw_rho_total%array(i, j, k)
    1102              : 
    1103       512000 :                         r2 = (/i*dr_pw(1), j*dr_pw(2), k*dr_pw(3)/) + origin
    1104       128000 :                         r_pbc = pbc(ra, r2, cell)
    1105              : 
    1106              :                         ! Store gradient d/dR_x w, including term: (r_x - R_x)
    1107              :                         cdft_control%group(igroup)%gradients_x(iatom, i, j, k) = &
    1108              :                            cdft_control%group(igroup)%gradients_x(iatom, i, j, k)* &
    1109       128000 :                            r_pbc(1)*exp_eval
    1110              : 
    1111              :                         ! Store gradient d/dR_y w, including term: (r_y - R_y)
    1112              :                         cdft_control%group(igroup)%gradients_y(iatom, i, j, k) = &
    1113              :                            cdft_control%group(igroup)%gradients_y(iatom, i, j, k)* &
    1114       128000 :                            r_pbc(2)*exp_eval
    1115              : 
    1116              :                         ! Store gradient d/dR_z w, including term:(r_z - R_z)
    1117              :                         cdft_control%group(igroup)%gradients_z(iatom, i, j, k) = &
    1118              :                            cdft_control%group(igroup)%gradients_z(iatom, i, j, k)* &
    1119       128000 :                            r_pbc(3)*exp_eval
    1120              : 
    1121              :                      END IF
    1122              :                   END DO
    1123              :                   END DO
    1124              :                END DO
    1125              :             END DO
    1126              :          END DO
    1127            2 :          CALL auxbas_pw_pool%give_back_pw(cdft_control%hw_rho_total)
    1128              :       END IF
    1129              : 
    1130           24 :       CALL rs_grid_release(rs_rho_all)
    1131              : 
    1132           24 :       IF (ALLOCATED(coefficients)) DEALLOCATE (coefficients)
    1133           24 :       IF (ALLOCATED(is_constraint)) DEALLOCATE (is_constraint)
    1134              : 
    1135           24 :       CALL timestop(handle)
    1136              : 
    1137           72 :    END SUBROUTINE hirshfeld_constraint_low
    1138              : 
    1139              : ! **************************************************************************************************
    1140              : !> \brief Calculates the value of a CDFT constraint by integrating the product of the CDFT
    1141              : !>        weight function and the realspace electron density
    1142              : !> \param qs_env ...
    1143              : ! **************************************************************************************************
    1144         3034 :    SUBROUTINE cdft_constraint_integrate(qs_env)
    1145              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1146              : 
    1147              :       CHARACTER(len=*), PARAMETER :: routineN = 'cdft_constraint_integrate'
    1148              : 
    1149              :       INTEGER                                            :: handle, i, iatom, igroup, ikind, ivar, &
    1150              :                                                             iw, jatom, natom, nvar
    1151              :       LOGICAL                                            :: is_becke, paw_atom
    1152              :       REAL(kind=dp)                                      :: dvol, eps_cavity, sign
    1153         3034 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: dE, strength, target_val
    1154         3034 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: electronic_charge, gapw_offset
    1155              :       TYPE(becke_constraint_type), POINTER               :: becke_control
    1156              :       TYPE(cdft_control_type), POINTER                   :: cdft_control
    1157         3034 :       TYPE(cdft_group_type), DIMENSION(:), POINTER       :: group
    1158              :       TYPE(cp_logger_type), POINTER                      :: logger
    1159              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1160              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1161         3034 :       TYPE(mpole_rho_atom), DIMENSION(:), POINTER        :: mp_rho
    1162         3034 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1163         3034 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: charge, rho_r
    1164              :       TYPE(qs_energy_type), POINTER                      :: energy
    1165         3034 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1166              :       TYPE(qs_rho_type), POINTER                         :: rho
    1167              :       TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
    1168              :       TYPE(section_vals_type), POINTER                   :: cdft_constraint_section
    1169              : 
    1170         3034 :       NULLIFY (para_env, dft_control, particle_set, rho_r, energy, rho, &
    1171         3034 :                logger, cdft_constraint_section, qs_kind_set, mp_rho, &
    1172         3034 :                rho0_mpole, group, charge)
    1173         3034 :       CALL timeset(routineN, handle)
    1174         3034 :       logger => cp_get_default_logger()
    1175              :       CALL get_qs_env(qs_env, &
    1176              :                       particle_set=particle_set, &
    1177              :                       rho=rho, &
    1178              :                       natom=natom, &
    1179              :                       dft_control=dft_control, &
    1180              :                       para_env=para_env, &
    1181         3034 :                       qs_kind_set=qs_kind_set)
    1182         3034 :       CALL qs_rho_get(rho, rho_r=rho_r)
    1183         3034 :       CPASSERT(ASSOCIATED(qs_kind_set))
    1184         3034 :       cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
    1185         3034 :       iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog")
    1186         3034 :       cdft_control => dft_control%qs_control%cdft_control
    1187         3034 :       is_becke = (cdft_control%type == outer_scf_becke_constraint)
    1188         3034 :       becke_control => cdft_control%becke_control
    1189         3034 :       IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) &
    1190            0 :          CPABORT("Becke control has not been allocated.")
    1191         3034 :       group => cdft_control%group
    1192              :       ! Initialize
    1193         3034 :       nvar = SIZE(cdft_control%target)
    1194         9102 :       ALLOCATE (strength(nvar))
    1195         6068 :       ALLOCATE (target_val(nvar))
    1196         6068 :       ALLOCATE (dE(nvar))
    1197         7044 :       strength(:) = cdft_control%strength(:)
    1198         7044 :       target_val(:) = cdft_control%target(:)
    1199         7044 :       sign = 1.0_dp
    1200         7044 :       dE = 0.0_dp
    1201         3034 :       dvol = group(1)%weight%pw_grid%dvol
    1202         3034 :       IF (cdft_control%atomic_charges) THEN
    1203         1598 :          charge => cdft_control%charge
    1204         6392 :          ALLOCATE (electronic_charge(cdft_control%natoms, dft_control%nspins))
    1205        11018 :          electronic_charge = 0.0_dp
    1206              :       END IF
    1207              :       ! Calculate value of constraint i.e. int ( rho(r) w(r) dr)
    1208         9102 :       DO i = 1, dft_control%nspins
    1209        14088 :          DO igroup = 1, SIZE(group)
    1210         8020 :             SELECT CASE (group(igroup)%constraint_type)
    1211              :             CASE (cdft_charge_constraint)
    1212           16 :                sign = 1.0_dp
    1213              :             CASE (cdft_magnetization_constraint)
    1214           16 :                IF (i == 1) THEN
    1215              :                   sign = 1.0_dp
    1216              :                ELSE
    1217            8 :                   sign = -1.0_dp
    1218              :                END IF
    1219              :             CASE (cdft_alpha_constraint)
    1220         1944 :                sign = 1.0_dp
    1221         1944 :                IF (i == 2) CYCLE
    1222              :             CASE (cdft_beta_constraint)
    1223         1944 :                sign = 1.0_dp
    1224         1944 :                IF (i == 1) CYCLE
    1225              :             CASE DEFAULT
    1226         8020 :                CPABORT("Unknown constraint type.")
    1227              :             END SELECT
    1228        12144 :             IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
    1229              :                ! With external control, we can use cavity_mat as a mask to kahan sum
    1230          180 :                eps_cavity = becke_control%eps_cavity
    1231          180 :                IF (igroup /= 1) &
    1232              :                   CALL cp_abort(__LOCATION__, &
    1233            0 :                                 "Multiple constraints not yet supported by parallel mixed calculations.")
    1234              :                dE(igroup) = dE(igroup) + sign*accurate_dot_product(group(igroup)%weight%array, rho_r(i)%array, &
    1235          180 :                                                                    becke_control%cavity_mat, eps_cavity)*dvol
    1236              :             ELSE
    1237         5896 :                dE(igroup) = dE(igroup) + sign*pw_integral_ab(group(igroup)%weight, rho_r(i), local_only=.TRUE.)
    1238              :             END IF
    1239              :          END DO
    1240         9102 :          IF (cdft_control%atomic_charges) THEN
    1241         9420 :             DO iatom = 1, cdft_control%natoms
    1242         9420 :                electronic_charge(iatom, i) = pw_integral_ab(charge(iatom), rho_r(i), local_only=.TRUE.)
    1243              :             END DO
    1244              :          END IF
    1245              :       END DO
    1246         3034 :       CALL get_qs_env(qs_env, energy=energy)
    1247         3034 :       CALL para_env%sum(dE)
    1248         3034 :       IF (cdft_control%atomic_charges) THEN
    1249         1598 :          CALL para_env%sum(electronic_charge)
    1250              :       END IF
    1251              :       ! Use fragment densities as reference value (= Becke deformation density)
    1252         3034 :       IF (cdft_control%fragment_density .AND. .NOT. cdft_control%fragments_integrated) THEN
    1253           10 :          CALL prepare_fragment_constraint(qs_env)
    1254              :       END IF
    1255         3034 :       IF (dft_control%qs_control%gapw) THEN
    1256              :          ! GAPW: add core charges (rho_hard - rho_soft)
    1257           46 :          IF (cdft_control%fragment_density) &
    1258              :             CALL cp_abort(__LOCATION__, &
    1259            0 :                           "Fragment constraints not yet compatible with GAPW.")
    1260          184 :          ALLOCATE (gapw_offset(nvar, dft_control%nspins))
    1261          230 :          gapw_offset = 0.0_dp
    1262           46 :          CALL get_qs_env(qs_env, rho0_mpole=rho0_mpole)
    1263           46 :          CALL get_rho0_mpole(rho0_mpole, mp_rho=mp_rho)
    1264          138 :          DO i = 1, dft_control%nspins
    1265          230 :             DO igroup = 1, SIZE(group)
    1266          460 :                DO iatom = 1, SIZE(group(igroup)%atoms)
    1267          276 :                   SELECT CASE (group(igroup)%constraint_type)
    1268              :                   CASE (cdft_charge_constraint)
    1269            0 :                      sign = 1.0_dp
    1270              :                   CASE (cdft_magnetization_constraint)
    1271            0 :                      IF (i == 1) THEN
    1272              :                         sign = 1.0_dp
    1273              :                      ELSE
    1274            0 :                         sign = -1.0_dp
    1275              :                      END IF
    1276              :                   CASE (cdft_alpha_constraint)
    1277            0 :                      sign = 1.0_dp
    1278            0 :                      IF (i == 2) CYCLE
    1279              :                   CASE (cdft_beta_constraint)
    1280            0 :                      sign = 1.0_dp
    1281            0 :                      IF (i == 1) CYCLE
    1282              :                   CASE DEFAULT
    1283          276 :                      CPABORT("Unknown constraint type.")
    1284              :                   END SELECT
    1285          276 :                   jatom = group(igroup)%atoms(iatom)
    1286          276 :                   CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind)
    1287          276 :                   CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
    1288          368 :                   IF (paw_atom) THEN
    1289          276 :                      gapw_offset(igroup, i) = gapw_offset(igroup, i) + sign*group(igroup)%coeff(iatom)*mp_rho(jatom)%q0(i)
    1290              :                   END IF
    1291              :                END DO
    1292              :             END DO
    1293              :          END DO
    1294           46 :          IF (cdft_control%atomic_charges) THEN
    1295          184 :             DO iatom = 1, cdft_control%natoms
    1296          138 :                jatom = cdft_control%atoms(iatom)
    1297          138 :                CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind)
    1298          138 :                CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
    1299          184 :                IF (paw_atom) THEN
    1300          414 :                   DO i = 1, dft_control%nspins
    1301          414 :                      electronic_charge(iatom, i) = electronic_charge(iatom, i) + mp_rho(jatom)%q0(i)
    1302              :                   END DO
    1303              :                END IF
    1304              :             END DO
    1305              :          END IF
    1306          138 :          DO i = 1, dft_control%nspins
    1307          230 :             DO ivar = 1, nvar
    1308          184 :                dE(ivar) = dE(ivar) + gapw_offset(ivar, i)
    1309              :             END DO
    1310              :          END DO
    1311           46 :          DEALLOCATE (gapw_offset)
    1312              :       END IF
    1313              :       ! Update constraint value and energy
    1314         7044 :       cdft_control%value(:) = dE(:)
    1315         3034 :       energy%cdft = 0.0_dp
    1316         7044 :       DO ivar = 1, nvar
    1317         7044 :          energy%cdft = energy%cdft + (dE(ivar) - target_val(ivar))*strength(ivar)
    1318              :       END DO
    1319              :       ! Print constraint info and atomic CDFT charges
    1320         3034 :       CALL cdft_constraint_print(qs_env, electronic_charge)
    1321              :       ! Deallocate tmp storage
    1322         3034 :       DEALLOCATE (dE, strength, target_val)
    1323         3034 :       IF (cdft_control%atomic_charges) DEALLOCATE (electronic_charge)
    1324         3034 :       CALL cp_print_key_finished_output(iw, logger, cdft_constraint_section, "PROGRAM_RUN_INFO")
    1325         3034 :       CALL timestop(handle)
    1326              : 
    1327         6068 :    END SUBROUTINE cdft_constraint_integrate
    1328              : 
    1329              : ! **************************************************************************************************
    1330              : !> \brief Calculates atomic forces due to a CDFT constraint (Becke or Hirshfeld)
    1331              : !> \param qs_env ...
    1332              : ! **************************************************************************************************
    1333          118 :    SUBROUTINE cdft_constraint_force(qs_env)
    1334              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1335              : 
    1336              :       CHARACTER(len=*), PARAMETER :: routineN = 'cdft_constraint_force'
    1337              : 
    1338              :       INTEGER                                            :: handle, i, iatom, igroup, ikind, ispin, &
    1339              :                                                             j, k, natom, nvar
    1340          118 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
    1341              :       INTEGER, DIMENSION(2, 3)                           :: bo
    1342              :       INTEGER, DIMENSION(3)                              :: lb, ub
    1343              :       REAL(kind=dp)                                      :: dvol, eps_cavity, sign
    1344          118 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: strength
    1345          118 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: cutoffs
    1346          118 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1347              :       TYPE(becke_constraint_type), POINTER               :: becke_control
    1348              :       TYPE(cdft_control_type), POINTER                   :: cdft_control
    1349          118 :       TYPE(cdft_group_type), DIMENSION(:), POINTER       :: group
    1350              :       TYPE(cell_type), POINTER                           :: cell
    1351              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1352              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1353          118 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1354          118 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
    1355          118 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
    1356              :       TYPE(qs_rho_type), POINTER                         :: rho
    1357              : 
    1358          118 :       CALL timeset(routineN, handle)
    1359          118 :       NULLIFY (atomic_kind_set, cell, para_env, dft_control, particle_set, &
    1360          118 :                rho, rho_r, force, cutoffs, becke_control, group)
    1361              : 
    1362              :       CALL get_qs_env(qs_env, &
    1363              :                       atomic_kind_set=atomic_kind_set, &
    1364              :                       natom=natom, &
    1365              :                       particle_set=particle_set, &
    1366              :                       cell=cell, &
    1367              :                       rho=rho, &
    1368              :                       force=force, &
    1369              :                       dft_control=dft_control, &
    1370          118 :                       para_env=para_env)
    1371          118 :       CALL qs_rho_get(rho, rho_r=rho_r)
    1372              : 
    1373          118 :       cdft_control => dft_control%qs_control%cdft_control
    1374          118 :       becke_control => cdft_control%becke_control
    1375          118 :       group => cdft_control%group
    1376          118 :       nvar = SIZE(cdft_control%target)
    1377          354 :       ALLOCATE (strength(nvar))
    1378          252 :       strength(:) = cdft_control%strength(:)
    1379          118 :       cutoffs => cdft_control%becke_control%cutoffs
    1380          118 :       eps_cavity = cdft_control%becke_control%eps_cavity
    1381              : 
    1382              :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
    1383              :                                atom_of_kind=atom_of_kind, &
    1384          118 :                                kind_of=kind_of)
    1385          252 :       DO igroup = 1, SIZE(cdft_control%group)
    1386          402 :          ALLOCATE (cdft_control%group(igroup)%integrated(3, natom))
    1387         1324 :          cdft_control%group(igroup)%integrated = 0.0_dp
    1388              :       END DO
    1389              : 
    1390          472 :       lb(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
    1391          472 :       ub(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
    1392         1180 :       bo = cdft_control%group(1)%weight%pw_grid%bounds_local
    1393          118 :       dvol = cdft_control%group(1)%weight%pw_grid%dvol
    1394          118 :       sign = 1.0_dp
    1395              : 
    1396          118 :       IF (cdft_control%type == outer_scf_becke_constraint) THEN
    1397          116 :          IF (.NOT. cdft_control%becke_control%in_memory) THEN
    1398           10 :             CALL becke_constraint_low(qs_env, just_gradients=.TRUE.)
    1399              :          END IF
    1400              : 
    1401            2 :       ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
    1402            2 :          IF (.NOT. cdft_control%in_memory) THEN
    1403            2 :             CALL hirshfeld_constraint_low(qs_env, just_gradients=.TRUE.)
    1404              :          END IF
    1405              :       END IF
    1406              : 
    1407              :       ! If no Becke Gaussian confinement
    1408          118 :       IF (.NOT. ASSOCIATED(becke_control%cavity_mat)) THEN
    1409              :          ! No external control
    1410         2106 :          DO k = bo(1, 1), bo(2, 1)
    1411        93674 :             DO j = bo(1, 2), bo(2, 2)
    1412      4518732 :                DO i = bo(1, 3), bo(2, 3)
    1413              :                   ! First check if this grid point should be skipped
    1414      4425152 :                   IF (cdft_control%becke_control%cavity_confine) THEN
    1415      4273152 :                      IF (cdft_control%becke_control%cavity%array(k, j, i) < eps_cavity) CYCLE
    1416              :                   END IF
    1417              : 
    1418      2713222 :                   DO igroup = 1, SIZE(cdft_control%group)
    1419      8512463 :                      DO iatom = 1, natom
    1420      9537059 :                         DO ispin = 1, dft_control%nspins
    1421              : 
    1422      5449748 :                            SELECT CASE (cdft_control%group(igroup)%constraint_type)
    1423              :                            CASE (cdft_charge_constraint)
    1424            0 :                               sign = 1.0_dp
    1425              :                            CASE (cdft_magnetization_constraint)
    1426            0 :                               IF (ispin == 1) THEN
    1427              :                                  sign = 1.0_dp
    1428              :                               ELSE
    1429            0 :                                  sign = -1.0_dp
    1430              :                               END IF
    1431              :                            CASE (cdft_alpha_constraint)
    1432       412880 :                               sign = 1.0_dp
    1433       412880 :                               IF (ispin == 2) CYCLE
    1434              :                            CASE (cdft_beta_constraint)
    1435       412880 :                               sign = 1.0_dp
    1436       412880 :                               IF (ispin == 1) CYCLE
    1437              :                            CASE DEFAULT
    1438      5449748 :                               CPABORT("Unknown constraint type.")
    1439              :                            END SELECT
    1440              : 
    1441      7761742 :                            IF (cdft_control%type == outer_scf_becke_constraint) THEN
    1442              : 
    1443              :                               cdft_control%group(igroup)%integrated(:, iatom) = &
    1444              :                                  cdft_control%group(igroup)%integrated(:, iatom) + sign* &
    1445              :                                  cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
    1446              :                                  *rho_r(ispin)%array(k, j, i) &
    1447     19123472 :                                  *dvol
    1448              : 
    1449       256000 :                            ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
    1450              : 
    1451              :                               cdft_control%group(igroup)%integrated(1, iatom) = &
    1452              :                                  cdft_control%group(igroup)%integrated(1, iatom) + sign* &
    1453              :                                  cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
    1454              :                                  *rho_r(ispin)%array(k, j, i) &
    1455       256000 :                                  *dvol
    1456              : 
    1457              :                               cdft_control%group(igroup)%integrated(2, iatom) = &
    1458              :                                  cdft_control%group(igroup)%integrated(2, iatom) + sign* &
    1459              :                                  cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
    1460              :                                  *rho_r(ispin)%array(k, j, i) &
    1461       256000 :                                  *dvol
    1462              : 
    1463              :                               cdft_control%group(igroup)%integrated(3, iatom) = &
    1464              :                                  cdft_control%group(igroup)%integrated(3, iatom) + sign* &
    1465              :                                  cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
    1466              :                                  *rho_r(ispin)%array(k, j, i) &
    1467       256000 :                                  *dvol
    1468              : 
    1469              :                            END IF
    1470              : 
    1471              :                         END DO
    1472              :                      END DO
    1473              :                   END DO
    1474              :                END DO
    1475              :             END DO
    1476              :          END DO
    1477              : 
    1478              :          ! If Becke Gaussian confinement
    1479              :       ELSE
    1480         1224 :          DO k = LBOUND(cdft_control%becke_control%cavity_mat, 1), UBOUND(cdft_control%becke_control%cavity_mat, 1)
    1481        61848 :             DO j = LBOUND(cdft_control%becke_control%cavity_mat, 2), UBOUND(cdft_control%becke_control%cavity_mat, 2)
    1482      2969728 :                DO i = LBOUND(cdft_control%becke_control%cavity_mat, 3), UBOUND(cdft_control%becke_control%cavity_mat, 3)
    1483              : 
    1484              :                   ! First check if this grid point should be skipped
    1485      2793472 :                   IF (cdft_control%becke_control%cavity_mat(k, j, i) < eps_cavity) CYCLE
    1486              : 
    1487      1747632 :                   DO igroup = 1, SIZE(group)
    1488      5327368 :                      DO iatom = 1, natom
    1489      5912424 :                         DO ispin = 1, dft_control%nspins
    1490      3378528 :                            SELECT CASE (group(igroup)%constraint_type)
    1491              :                            CASE (cdft_charge_constraint)
    1492            0 :                               sign = 1.0_dp
    1493              :                            CASE (cdft_magnetization_constraint)
    1494            0 :                               IF (ispin == 1) THEN
    1495              :                                  sign = 1.0_dp
    1496              :                               ELSE
    1497            0 :                                  sign = -1.0_dp
    1498              :                               END IF
    1499              :                            CASE (cdft_alpha_constraint)
    1500            0 :                               sign = 1.0_dp
    1501            0 :                               IF (ispin == 2) CYCLE
    1502              :                            CASE (cdft_beta_constraint)
    1503            0 :                               sign = 1.0_dp
    1504            0 :                               IF (ispin == 1) CYCLE
    1505              :                            CASE DEFAULT
    1506      3378528 :                               CPABORT("Unknown constraint type.")
    1507              :                            END SELECT
    1508              : 
    1509              :                            ! Integrate gradient of weight function
    1510      5067792 :                            IF (cdft_control%type == outer_scf_becke_constraint) THEN
    1511              : 
    1512              :                               cdft_control%group(igroup)%integrated(:, iatom) = &
    1513              :                                  cdft_control%group(igroup)%integrated(:, iatom) + sign* &
    1514              :                                  cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
    1515              :                                  *rho_r(ispin)%array(k, j, i) &
    1516     13514112 :                                  *dvol
    1517              : 
    1518            0 :                            ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
    1519              : 
    1520              :                               cdft_control%group(igroup)%integrated(1, iatom) = &
    1521              :                                  cdft_control%group(igroup)%integrated(1, iatom) + sign* &
    1522              :                                  cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
    1523              :                                  *rho_r(ispin)%array(k, j, i) &
    1524            0 :                                  *dvol
    1525              : 
    1526              :                               cdft_control%group(igroup)%integrated(2, iatom) = &
    1527              :                                  cdft_control%group(igroup)%integrated(2, iatom) + sign* &
    1528              :                                  cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
    1529              :                                  *rho_r(ispin)%array(k, j, i) &
    1530            0 :                                  *dvol
    1531              : 
    1532              :                               cdft_control%group(igroup)%integrated(3, iatom) = &
    1533              :                                  cdft_control%group(igroup)%integrated(3, iatom) + sign* &
    1534              :                                  cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
    1535              :                                  *rho_r(ispin)%array(k, j, i) &
    1536            0 :                                  *dvol
    1537              : 
    1538              :                            END IF
    1539              : 
    1540              :                         END DO
    1541              :                      END DO
    1542              :                   END DO
    1543              :                END DO
    1544              :             END DO
    1545              :          END DO
    1546              :       END IF
    1547              : 
    1548          118 :       IF (.NOT. cdft_control%transfer_pot) THEN
    1549           98 :          IF (cdft_control%type == outer_scf_becke_constraint) THEN
    1550          208 :             DO igroup = 1, SIZE(group)
    1551          208 :                DEALLOCATE (cdft_control%group(igroup)%gradients)
    1552              :             END DO
    1553            2 :          ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
    1554            4 :             DO igroup = 1, SIZE(group)
    1555            2 :                DEALLOCATE (cdft_control%group(igroup)%gradients_x)
    1556            2 :                DEALLOCATE (cdft_control%group(igroup)%gradients_y)
    1557            4 :                DEALLOCATE (cdft_control%group(igroup)%gradients_z)
    1558              :             END DO
    1559              :          END IF
    1560              :       END IF
    1561              : 
    1562          252 :       DO igroup = 1, SIZE(group)
    1563         2396 :          CALL para_env%sum(group(igroup)%integrated)
    1564              :       END DO
    1565              : 
    1566              :       ! Update force only on master process. Otherwise force due to constraint becomes multiplied
    1567              :       ! by the number of processes when the final force%rho_elec is constructed in qs_force
    1568              :       ! by mp_summing [the final integrated(:,:) is distributed on all processors]
    1569          118 :       IF (para_env%is_source()) THEN
    1570          150 :          DO igroup = 1, SIZE(group)
    1571          308 :             DO iatom = 1, natom
    1572          158 :                ikind = kind_of(iatom)
    1573          158 :                i = atom_of_kind(iatom)
    1574         1185 :                force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + group(igroup)%integrated(:, iatom)*strength(igroup)
    1575              :             END DO
    1576              :          END DO
    1577              :       END IF
    1578              : 
    1579          118 :       DEALLOCATE (strength)
    1580          252 :       DO igroup = 1, SIZE(group)
    1581          252 :          DEALLOCATE (group(igroup)%integrated)
    1582              :       END DO
    1583          118 :       NULLIFY (group)
    1584              : 
    1585          118 :       CALL timestop(handle)
    1586              : 
    1587          236 :    END SUBROUTINE cdft_constraint_force
    1588              : 
    1589              : ! **************************************************************************************************
    1590              : !> \brief Prepare CDFT fragment constraints. Fragment densities are read from cube files, multiplied
    1591              : !>        by the CDFT weight functions and integrated over the realspace grid.
    1592              : !> \param qs_env ...
    1593              : ! **************************************************************************************************
    1594           10 :    SUBROUTINE prepare_fragment_constraint(qs_env)
    1595              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1596              : 
    1597              :       CHARACTER(len=*), PARAMETER :: routineN = 'prepare_fragment_constraint'
    1598              : 
    1599              :       INTEGER                                            :: handle, i, iatom, igroup, natom, &
    1600              :                                                             nelectron_total, nfrag_spins
    1601              :       LOGICAL                                            :: is_becke, needs_spin_density
    1602              :       REAL(kind=dp)                                      :: dvol, multiplier(2), nelectron_frag
    1603              :       TYPE(becke_constraint_type), POINTER               :: becke_control
    1604              :       TYPE(cdft_control_type), POINTER                   :: cdft_control
    1605           10 :       TYPE(cdft_group_type), DIMENSION(:), POINTER       :: group
    1606              :       TYPE(cp_logger_type), POINTER                      :: logger
    1607              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1608              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1609              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1610              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1611           10 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:)    :: rho_frag
    1612              :       TYPE(qs_subsys_type), POINTER                      :: subsys
    1613              : 
    1614           10 :       NULLIFY (para_env, dft_control, logger, subsys, pw_env, auxbas_pw_pool, group)
    1615           10 :       CALL timeset(routineN, handle)
    1616           10 :       logger => cp_get_default_logger()
    1617              :       CALL get_qs_env(qs_env, &
    1618              :                       natom=natom, &
    1619              :                       dft_control=dft_control, &
    1620           10 :                       para_env=para_env)
    1621              : 
    1622           10 :       cdft_control => dft_control%qs_control%cdft_control
    1623           10 :       is_becke = (cdft_control%type == outer_scf_becke_constraint)
    1624           10 :       becke_control => cdft_control%becke_control
    1625           10 :       IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) &
    1626            0 :          CPABORT("Becke control has not been allocated.")
    1627           10 :       group => cdft_control%group
    1628           10 :       dvol = group(1)%weight%pw_grid%dvol
    1629              :       ! Fragment densities are meaningful only for some calculation types
    1630           10 :       IF (.NOT. qs_env%single_point_run) &
    1631              :          CALL cp_abort(__LOCATION__, &
    1632              :                        "CDFT fragment constraints are only compatible with single "// &
    1633            0 :                        "point calculations (run_type ENERGY or ENERGY_FORCE).")
    1634           10 :       IF (dft_control%qs_control%gapw) &
    1635              :          CALL cp_abort(__LOCATION__, &
    1636            0 :                        "CDFT fragment constraint not compatible with GAPW.")
    1637           30 :       needs_spin_density = .FALSE.
    1638           30 :       multiplier = 1.0_dp
    1639           10 :       nfrag_spins = 1
    1640           22 :       DO igroup = 1, SIZE(group)
    1641           10 :          SELECT CASE (group(igroup)%constraint_type)
    1642              :          CASE (cdft_charge_constraint)
    1643              :             ! Do nothing
    1644              :          CASE (cdft_magnetization_constraint)
    1645            6 :             needs_spin_density = .TRUE.
    1646              :          CASE (cdft_alpha_constraint, cdft_beta_constraint)
    1647              :             CALL cp_abort(__LOCATION__, &
    1648              :                           "CDFT fragment constraint not yet compatible with "// &
    1649            0 :                           "spin specific constraints.")
    1650              :          CASE DEFAULT
    1651           12 :             CPABORT("Unknown constraint type.")
    1652              :          END SELECT
    1653              :       END DO
    1654           10 :       IF (needs_spin_density) THEN
    1655           12 :          nfrag_spins = 2
    1656           12 :          DO i = 1, 2
    1657           12 :             IF (cdft_control%flip_fragment(i)) multiplier(i) = -1.0_dp
    1658              :          END DO
    1659              :       END IF
    1660              :       ! Read fragment reference densities
    1661           68 :       ALLOCATE (cdft_control%fragments(nfrag_spins, 2))
    1662           34 :       ALLOCATE (rho_frag(nfrag_spins))
    1663           10 :       CALL get_qs_env(qs_env, pw_env=pw_env)
    1664           10 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
    1665              :       ! Total density (rho_alpha + rho_beta)
    1666           10 :       CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 1))
    1667              :       CALL cp_cube_to_pw(cdft_control%fragments(1, 1), &
    1668           10 :                          cdft_control%fragment_a_fname, 1.0_dp)
    1669           10 :       CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 2))
    1670              :       CALL cp_cube_to_pw(cdft_control%fragments(1, 2), &
    1671           10 :                          cdft_control%fragment_b_fname, 1.0_dp)
    1672              :       ! Spin difference density (rho_alpha - rho_beta) if needed
    1673           10 :       IF (needs_spin_density) THEN
    1674            4 :          CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 1))
    1675              :          CALL cp_cube_to_pw(cdft_control%fragments(2, 1), &
    1676            4 :                             cdft_control%fragment_a_spin_fname, multiplier(1))
    1677            4 :          CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 2))
    1678              :          CALL cp_cube_to_pw(cdft_control%fragments(2, 2), &
    1679            4 :                             cdft_control%fragment_b_spin_fname, multiplier(2))
    1680              :       END IF
    1681              :       ! Sum up fragments
    1682           24 :       DO i = 1, nfrag_spins
    1683           14 :          CALL auxbas_pw_pool%create_pw(rho_frag(i))
    1684           14 :          CALL pw_copy(cdft_control%fragments(i, 1), rho_frag(i))
    1685           14 :          CALL pw_axpy(cdft_control%fragments(i, 2), rho_frag(i), 1.0_dp)
    1686           14 :          CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 1))
    1687           24 :          CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 2))
    1688              :       END DO
    1689           10 :       DEALLOCATE (cdft_control%fragments)
    1690              :       ! Check that the number of electrons is consistent
    1691           10 :       CALL get_qs_env(qs_env, subsys=subsys)
    1692           10 :       CALL qs_subsys_get(subsys, nelectron_total=nelectron_total)
    1693           10 :       nelectron_frag = pw_integrate_function(rho_frag(1))
    1694           10 :       IF (NINT(nelectron_frag) /= nelectron_total) &
    1695              :          CALL cp_abort(__LOCATION__, &
    1696              :                        "The number of electrons in the reference and interacting "// &
    1697            0 :                        "configurations does not match. Check your fragment cube files.")
    1698              :       ! Update constraint target value i.e. perform integration w_i*rho_frag_{tot/spin}*dr
    1699           22 :       cdft_control%target = 0.0_dp
    1700           22 :       DO igroup = 1, SIZE(group)
    1701           12 :          IF (group(igroup)%constraint_type == cdft_charge_constraint) THEN
    1702              :             i = 1
    1703              :          ELSE
    1704            6 :             i = 2
    1705              :          END IF
    1706           22 :          IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
    1707              :             cdft_control%target(igroup) = cdft_control%target(igroup) + &
    1708              :                                           accurate_dot_product(group(igroup)%weight%array, rho_frag(i)%array, &
    1709            0 :                                                                becke_control%cavity_mat, becke_control%eps_cavity)*dvol
    1710              :          ELSE
    1711              :             cdft_control%target(igroup) = cdft_control%target(igroup) + &
    1712           12 :                                           pw_integral_ab(group(igroup)%weight, rho_frag(i), local_only=.TRUE.)
    1713              :          END IF
    1714              :       END DO
    1715           34 :       CALL para_env%sum(cdft_control%target)
    1716              :       ! Calculate reference atomic charges int( w_i * rho_frag * dr )
    1717           10 :       IF (cdft_control%atomic_charges) THEN
    1718           40 :          ALLOCATE (cdft_control%charges_fragment(cdft_control%natoms, nfrag_spins))
    1719           24 :          DO i = 1, nfrag_spins
    1720           44 :             DO iatom = 1, cdft_control%natoms
    1721              :                cdft_control%charges_fragment(iatom, i) = &
    1722           34 :                   pw_integral_ab(cdft_control%charge(iatom), rho_frag(i), local_only=.TRUE.)
    1723              :             END DO
    1724              :          END DO
    1725           78 :          CALL para_env%sum(cdft_control%charges_fragment)
    1726              :       END IF
    1727           24 :       DO i = 1, nfrag_spins
    1728           24 :          CALL auxbas_pw_pool%give_back_pw(rho_frag(i))
    1729              :       END DO
    1730           10 :       DEALLOCATE (rho_frag)
    1731           10 :       cdft_control%fragments_integrated = .TRUE.
    1732              : 
    1733           10 :       CALL timestop(handle)
    1734              : 
    1735           10 :    END SUBROUTINE prepare_fragment_constraint
    1736              : 
    1737              : END MODULE qs_cdft_methods
        

Generated by: LCOV version 2.0-1