Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief 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_operations,             ONLY: copy_dbcsr_to_fm,&
      14             :                                               copy_fm_to_dbcsr
      15             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      16             :                                               cp_fm_type
      17             :    USE dbcsr_api,                       ONLY: dbcsr_p_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       96451 :    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       96451 :       NULLIFY (dft_control)
      70       96451 :       CALL get_qs_env(qs_env, dft_control=dft_control)
      71       96451 :       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       96451 :    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       96451 :    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       96451 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ksmat, rho_ao
     148             : 
     149       96451 :       energy%mulliken = 0.0_dp
     150             : 
     151       96451 :       IF (dft_control%qs_control%mulliken_restraint) THEN
     152             : 
     153             :          ! Test no k-points
     154          30 :          CPASSERT(SIZE(matrix_s, 2) == 1)
     155             : 
     156          30 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     157             : 
     158          30 :          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          12 :                                     order_p=mulliken_order_p)
     162             :          ELSE
     163          18 :             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          18 :                                     ks_matrix=ksmat, order_p=mulliken_order_p)
     167             :          END IF
     168             : 
     169             :       END IF
     170             : 
     171       96451 :    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       96451 :    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       96451 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_mo_derivs
     193             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     194       96451 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mo_derivs, smat
     195       96451 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     196             : 
     197       96451 :       NULLIFY (mo_array, mo_coeff, mo_derivs)
     198             : 
     199       96451 :       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       96451 :          energy%s2_restraint = 0.0_dp
     228             :       END IF
     229      192902 :    END SUBROUTINE qs_ks_s2_restraint
     230             : 
     231             : END MODULE qs_ks_apply_restraints