Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 121829 : 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 121829 : NULLIFY (dft_control)
70 121829 : CALL get_qs_env(qs_env, dft_control=dft_control)
71 121829 : IF (dft_control%qs_control%cdft) THEN
72 3410 : cdft_control => dft_control%qs_control%cdft_control
73 : ! Test no k-points
74 3410 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
75 3410 : IF (do_kpoints) CPABORT("CDFT constraints with k-points not supported.")
76 :
77 6820 : SELECT CASE (cdft_control%type)
78 : CASE (outer_scf_becke_constraint, outer_scf_hirshfeld_constraint)
79 3410 : IF (cdft_control%need_pot) THEN
80 : ! First SCF iteraration => allocate storage
81 476 : DO igroup = 1, SIZE(cdft_control%group)
82 252 : ALLOCATE (cdft_control%group(igroup)%weight)
83 252 : 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 476 : .AND. dft_control%nspins == 1) THEN
87 : CALL cp_abort(__LOCATION__, &
88 0 : "Spin constraints require a spin polarized calculation.")
89 : END IF
90 : END DO
91 224 : IF (cdft_control%atomic_charges) THEN
92 110 : IF (.NOT. ASSOCIATED(cdft_control%charge)) THEN
93 40 : ALLOCATE (cdft_control%charge(cdft_control%natoms))
94 : END IF
95 334 : DO iatom = 1, cdft_control%natoms
96 334 : CALL auxbas_pw_pool%create_pw(cdft_control%charge(iatom))
97 : END DO
98 : END IF
99 : ! Another sanity check
100 224 : CALL get_qs_env(qs_env, natom=natom)
101 224 : IF (natom < cdft_control%natoms) THEN
102 : CALL cp_abort(__LOCATION__, &
103 0 : "The number of constraint atoms exceeds the total number of atoms.")
104 : END IF
105 : ELSE
106 7320 : DO igroup = 1, SIZE(cdft_control%group)
107 4134 : inv_vol = 1.0_dp/cdft_control%group(igroup)%weight%pw_grid%dvol
108 7320 : CALL pw_scale(cdft_control%group(igroup)%weight, inv_vol)
109 : END DO
110 : END IF
111 : ! Build/Integrate CDFT constraints with selected population analysis method
112 3410 : IF (cdft_control%type == outer_scf_becke_constraint) THEN
113 3324 : CALL becke_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
114 86 : ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
115 86 : CALL hirshfeld_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
116 : END IF
117 7796 : DO igroup = 1, SIZE(cdft_control%group)
118 7796 : CALL pw_scale(cdft_control%group(igroup)%weight, cdft_control%group(igroup)%weight%pw_grid%dvol)
119 : END DO
120 3410 : IF (cdft_control%need_pot) cdft_control%need_pot = .FALSE.
121 : CASE DEFAULT
122 3410 : CPABORT("Unknown constraint type.")
123 : END SELECT
124 : END IF
125 :
126 121829 : END SUBROUTINE qs_ks_cdft_constraint
127 :
128 : ! **************************************************************************************************
129 : !> \brief ...
130 : !> \param energy ...
131 : !> \param dft_control ...
132 : !> \param just_energy ...
133 : !> \param para_env ...
134 : !> \param ks_matrix ...
135 : !> \param matrix_s ...
136 : !> \param rho ...
137 : !> \param mulliken_order_p ...
138 : ! **************************************************************************************************
139 121829 : SUBROUTINE qs_ks_mulliken_restraint(energy, dft_control, just_energy, para_env, &
140 : ks_matrix, matrix_s, rho, mulliken_order_p)
141 :
142 : TYPE(qs_energy_type), POINTER :: energy
143 : TYPE(dft_control_type), POINTER :: dft_control
144 : LOGICAL, INTENT(in) :: just_energy
145 : TYPE(mp_para_env_type), POINTER :: para_env
146 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_s
147 : TYPE(qs_rho_type), POINTER :: rho
148 : REAL(KIND=dp) :: mulliken_order_p
149 :
150 121829 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, rho_ao
151 :
152 121829 : energy%mulliken = 0.0_dp
153 :
154 121829 : IF (dft_control%qs_control%mulliken_restraint) THEN
155 :
156 : ! Test no k-points
157 48 : CPASSERT(SIZE(matrix_s, 2) == 1)
158 :
159 48 : CALL qs_rho_get(rho, rho_ao=rho_ao)
160 :
161 48 : IF (just_energy) THEN
162 : CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
163 : para_env, matrix_s(1, 1)%matrix, rho_ao, energy=energy%mulliken, &
164 18 : order_p=mulliken_order_p)
165 : ELSE
166 30 : ksmat => ks_matrix(:, 1)
167 : CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
168 : para_env, matrix_s(1, 1)%matrix, rho_ao, energy=energy%mulliken, &
169 30 : ks_matrix=ksmat, order_p=mulliken_order_p)
170 : END IF
171 :
172 : END IF
173 :
174 121829 : END SUBROUTINE qs_ks_mulliken_restraint
175 :
176 : ! **************************************************************************************************
177 : !> \brief ...
178 : !> \param dft_control ...
179 : !> \param qs_env ...
180 : !> \param matrix_s ...
181 : !> \param energy ...
182 : !> \param calculate_forces ...
183 : !> \param just_energy ...
184 : ! **************************************************************************************************
185 121829 : SUBROUTINE qs_ks_s2_restraint(dft_control, qs_env, matrix_s, &
186 : energy, calculate_forces, just_energy)
187 :
188 : TYPE(dft_control_type), POINTER :: dft_control
189 : TYPE(qs_environment_type), POINTER :: qs_env
190 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
191 : TYPE(qs_energy_type), POINTER :: energy
192 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
193 :
194 : INTEGER :: i
195 121829 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mo_derivs
196 : TYPE(cp_fm_type), POINTER :: mo_coeff
197 121829 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs, smat
198 121829 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
199 :
200 121829 : NULLIFY (mo_array, mo_coeff, mo_derivs)
201 :
202 121829 : IF (dft_control%qs_control%s2_restraint) THEN
203 : ! Test no k-points
204 0 : CPASSERT(SIZE(matrix_s, 2) == 1)
205 : ! adds s2_restraint energy and orbital derivatives
206 0 : CPASSERT(dft_control%nspins == 2)
207 0 : CPASSERT(qs_env%requires_mo_derivs)
208 : ! forces are not implemented (not difficult, but ... )
209 0 : CPASSERT(.NOT. calculate_forces)
210 : MARK_USED(calculate_forces)
211 0 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
212 :
213 0 : ALLOCATE (fm_mo_derivs(SIZE(mo_derivs, 1))) !fm->dbcsr
214 0 : DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
215 0 : CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff) !fm->dbcsr
216 0 : CALL cp_fm_create(fm_mo_derivs(i), mo_coeff%matrix_struct) !fm->dbcsr
217 0 : CALL copy_dbcsr_to_fm(mo_derivs(i)%matrix, fm_mo_derivs(i)) !fm->dbcsr
218 : END DO !fm->dbcsr
219 :
220 0 : smat => matrix_s(:, 1)
221 : CALL s2_restraint(mo_array, smat, fm_mo_derivs, energy%s2_restraint, &
222 0 : dft_control%qs_control%s2_restraint_control, just_energy)
223 :
224 0 : DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
225 0 : CALL copy_fm_to_dbcsr(fm_mo_derivs(i), mo_derivs(i)%matrix) !fm->dbcsr
226 : END DO !fm->dbcsr
227 0 : DEALLOCATE (fm_mo_derivs) !fm->dbcsr
228 :
229 : ELSE
230 121829 : energy%s2_restraint = 0.0_dp
231 : END IF
232 243658 : END SUBROUTINE qs_ks_s2_restraint
233 :
234 : END MODULE qs_ks_apply_restraints
|
