LCOV - code coverage report
Current view: top level - src - qs_ks_apply_restraints.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 76.5 % 68 52
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 3 3

            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 Set of routines to apply restraints to the KS hamiltonian
      10              : ! **************************************************************************************************
      11              : MODULE qs_ks_apply_restraints
      12              :    USE cp_control_types,                ONLY: dft_control_type
      13              :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type
      14              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      15              :                                               copy_fm_to_dbcsr
      16              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      17              :                                               cp_fm_type
      18              :    USE input_constants,                 ONLY: cdft_charge_constraint,&
      19              :                                               outer_scf_becke_constraint,&
      20              :                                               outer_scf_hirshfeld_constraint
      21              :    USE kinds,                           ONLY: dp
      22              :    USE message_passing,                 ONLY: mp_para_env_type
      23              :    USE mulliken,                        ONLY: mulliken_restraint
      24              :    USE pw_methods,                      ONLY: pw_scale
      25              :    USE pw_pool_types,                   ONLY: pw_pool_type
      26              :    USE qs_cdft_methods,                 ONLY: becke_constraint,&
      27              :                                               hirshfeld_constraint
      28              :    USE qs_cdft_types,                   ONLY: cdft_control_type
      29              :    USE qs_energy_types,                 ONLY: qs_energy_type
      30              :    USE qs_environment_types,            ONLY: get_qs_env,&
      31              :                                               qs_environment_type
      32              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      33              :                                               mo_set_type
      34              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      35              :                                               qs_rho_type
      36              :    USE s_square_methods,                ONLY: s2_restraint
      37              : #include "./base/base_uses.f90"
      38              : 
      39              :    IMPLICIT NONE
      40              : 
      41              :    PRIVATE
      42              : 
      43              :    LOGICAL, PARAMETER :: debug_this_module = .TRUE.
      44              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_apply_restraints'
      45              : 
      46              :    PUBLIC :: qs_ks_mulliken_restraint, qs_ks_s2_restraint
      47              :    PUBLIC :: qs_ks_cdft_constraint
      48              : 
      49              : CONTAINS
      50              : 
      51              : ! **************************************************************************************************
      52              : !> \brief Apply a CDFT constraint
      53              : !> \param qs_env the qs_env where to apply the constraint
      54              : !> \param auxbas_pw_pool the pool that owns the real space grid where the CDFT potential is defined
      55              : !> \param calculate_forces if forces should be calculated
      56              : !> \param cdft_control the CDFT control type
      57              : ! **************************************************************************************************
      58       103397 :    SUBROUTINE qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
      59              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      60              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
      61              :       LOGICAL, INTENT(in)                                :: calculate_forces
      62              :       TYPE(cdft_control_type), POINTER                   :: cdft_control
      63              : 
      64              :       INTEGER                                            :: iatom, igroup, natom
      65              :       LOGICAL                                            :: do_kpoints
      66              :       REAL(KIND=dp)                                      :: inv_vol
      67              :       TYPE(dft_control_type), POINTER                    :: dft_control
      68              : 
      69       103397 :       NULLIFY (dft_control)
      70       103397 :       CALL get_qs_env(qs_env, dft_control=dft_control)
      71       103397 :       IF (dft_control%qs_control%cdft) THEN
      72         3034 :          cdft_control => dft_control%qs_control%cdft_control
      73              :          ! Test no k-points
      74         3034 :          CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
      75         3034 :          IF (do_kpoints) CPABORT("CDFT constraints with k-points not supported.")
      76              : 
      77         6068 :          SELECT CASE (cdft_control%type)
      78              :          CASE (outer_scf_becke_constraint, outer_scf_hirshfeld_constraint)
      79         3034 :             IF (cdft_control%need_pot) THEN
      80              :                ! First SCF iteraration => allocate storage
      81          452 :                DO igroup = 1, SIZE(cdft_control%group)
      82          240 :                   ALLOCATE (cdft_control%group(igroup)%weight)
      83          240 :                   CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%weight)
      84              :                   ! Sanity check
      85              :                   IF (cdft_control%group(igroup)%constraint_type /= cdft_charge_constraint &
      86          240 :                       .AND. dft_control%nspins == 1) &
      87              :                      CALL cp_abort(__LOCATION__, &
      88          212 :                                    "Spin constraints require a spin polarized calculation.")
      89              :                END DO
      90          212 :                IF (cdft_control%atomic_charges) THEN
      91          110 :                   IF (.NOT. ASSOCIATED(cdft_control%charge)) &
      92           40 :                      ALLOCATE (cdft_control%charge(cdft_control%natoms))
      93          334 :                   DO iatom = 1, cdft_control%natoms
      94          334 :                      CALL auxbas_pw_pool%create_pw(cdft_control%charge(iatom))
      95              :                   END DO
      96              :                END IF
      97              :                ! Another sanity check
      98          212 :                CALL get_qs_env(qs_env, natom=natom)
      99          212 :                IF (natom < cdft_control%natoms) &
     100              :                   CALL cp_abort(__LOCATION__, &
     101            0 :                                 "The number of constraint atoms exceeds the total number of atoms.")
     102              :             ELSE
     103         6592 :                DO igroup = 1, SIZE(cdft_control%group)
     104         3770 :                   inv_vol = 1.0_dp/cdft_control%group(igroup)%weight%pw_grid%dvol
     105         6592 :                   CALL pw_scale(cdft_control%group(igroup)%weight, inv_vol)
     106              :                END DO
     107              :             END IF
     108              :             ! Build/Integrate CDFT constraints with selected population analysis method
     109         3034 :             IF (cdft_control%type == outer_scf_becke_constraint) THEN
     110         2948 :                CALL becke_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
     111           86 :             ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
     112           86 :                CALL hirshfeld_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
     113              :             END IF
     114         7044 :             DO igroup = 1, SIZE(cdft_control%group)
     115         7044 :                CALL pw_scale(cdft_control%group(igroup)%weight, cdft_control%group(igroup)%weight%pw_grid%dvol)
     116              :             END DO
     117         3034 :             IF (cdft_control%need_pot) cdft_control%need_pot = .FALSE.
     118              :          CASE DEFAULT
     119         3034 :             CPABORT("Unknown constraint type.")
     120              :          END SELECT
     121              :       END IF
     122              : 
     123       103397 :    END SUBROUTINE qs_ks_cdft_constraint
     124              : 
     125              : ! **************************************************************************************************
     126              : !> \brief ...
     127              : !> \param energy ...
     128              : !> \param dft_control ...
     129              : !> \param just_energy ...
     130              : !> \param para_env ...
     131              : !> \param ks_matrix ...
     132              : !> \param matrix_s ...
     133              : !> \param rho ...
     134              : !> \param mulliken_order_p ...
     135              : ! **************************************************************************************************
     136       103397 :    SUBROUTINE qs_ks_mulliken_restraint(energy, dft_control, just_energy, para_env, &
     137              :                                        ks_matrix, matrix_s, rho, mulliken_order_p)
     138              : 
     139              :       TYPE(qs_energy_type), POINTER                      :: energy
     140              :       TYPE(dft_control_type), POINTER                    :: dft_control
     141              :       LOGICAL, INTENT(in)                                :: just_energy
     142              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     143              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: ks_matrix, matrix_s
     144              :       TYPE(qs_rho_type), POINTER                         :: rho
     145              :       REAL(KIND=dp)                                      :: mulliken_order_p
     146              : 
     147       103397 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ksmat, rho_ao
     148              : 
     149       103397 :       energy%mulliken = 0.0_dp
     150              : 
     151       103397 :       IF (dft_control%qs_control%mulliken_restraint) THEN
     152              : 
     153              :          ! Test no k-points
     154           48 :          CPASSERT(SIZE(matrix_s, 2) == 1)
     155              : 
     156           48 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     157              : 
     158           48 :          IF (just_energy) THEN
     159              :             CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
     160              :                                     para_env, matrix_s(1, 1)%matrix, rho_ao, energy=energy%mulliken, &
     161           18 :                                     order_p=mulliken_order_p)
     162              :          ELSE
     163           30 :             ksmat => ks_matrix(:, 1)
     164              :             CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
     165              :                                     para_env, matrix_s(1, 1)%matrix, rho_ao, energy=energy%mulliken, &
     166           30 :                                     ks_matrix=ksmat, order_p=mulliken_order_p)
     167              :          END IF
     168              : 
     169              :       END IF
     170              : 
     171       103397 :    END SUBROUTINE qs_ks_mulliken_restraint
     172              : 
     173              : ! **************************************************************************************************
     174              : !> \brief ...
     175              : !> \param dft_control ...
     176              : !> \param qs_env ...
     177              : !> \param matrix_s ...
     178              : !> \param energy ...
     179              : !> \param calculate_forces ...
     180              : !> \param just_energy ...
     181              : ! **************************************************************************************************
     182       103397 :    SUBROUTINE qs_ks_s2_restraint(dft_control, qs_env, matrix_s, &
     183              :                                  energy, calculate_forces, just_energy)
     184              : 
     185              :       TYPE(dft_control_type), POINTER                    :: dft_control
     186              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     187              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s
     188              :       TYPE(qs_energy_type), POINTER                      :: energy
     189              :       LOGICAL, INTENT(in)                                :: calculate_forces, just_energy
     190              : 
     191              :       INTEGER                                            :: i
     192       103397 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_mo_derivs
     193              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     194       103397 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mo_derivs, smat
     195       103397 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     196              : 
     197       103397 :       NULLIFY (mo_array, mo_coeff, mo_derivs)
     198              : 
     199       103397 :       IF (dft_control%qs_control%s2_restraint) THEN
     200              :          ! Test no k-points
     201            0 :          CPASSERT(SIZE(matrix_s, 2) == 1)
     202              :          ! adds s2_restraint energy and orbital derivatives
     203            0 :          CPASSERT(dft_control%nspins == 2)
     204            0 :          CPASSERT(qs_env%requires_mo_derivs)
     205              :          ! forces are not implemented (not difficult, but ... )
     206            0 :          CPASSERT(.NOT. calculate_forces)
     207              :          MARK_USED(calculate_forces)
     208            0 :          CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
     209              : 
     210            0 :          ALLOCATE (fm_mo_derivs(SIZE(mo_derivs, 1))) !fm->dbcsr
     211            0 :          DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
     212            0 :             CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff) !fm->dbcsr
     213            0 :             CALL cp_fm_create(fm_mo_derivs(i), mo_coeff%matrix_struct) !fm->dbcsr
     214            0 :             CALL copy_dbcsr_to_fm(mo_derivs(i)%matrix, fm_mo_derivs(i)) !fm->dbcsr
     215              :          END DO !fm->dbcsr
     216              : 
     217            0 :          smat => matrix_s(:, 1)
     218              :          CALL s2_restraint(mo_array, smat, fm_mo_derivs, energy%s2_restraint, &
     219            0 :                            dft_control%qs_control%s2_restraint_control, just_energy)
     220              : 
     221            0 :          DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
     222            0 :             CALL copy_fm_to_dbcsr(fm_mo_derivs(i), mo_derivs(i)%matrix) !fm->dbcsr
     223              :          END DO !fm->dbcsr
     224            0 :          DEALLOCATE (fm_mo_derivs) !fm->dbcsr
     225              : 
     226              :       ELSE
     227       103397 :          energy%s2_restraint = 0.0_dp
     228              :       END IF
     229       206794 :    END SUBROUTINE qs_ks_s2_restraint
     230              : 
     231              : END MODULE qs_ks_apply_restraints
        

Generated by: LCOV version 2.0-1