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