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_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 96949 : 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 96949 : NULLIFY (dft_control)
70 96949 : CALL get_qs_env(qs_env, dft_control=dft_control)
71 96949 : 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 96949 : 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 96949 : 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 96949 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, rho_ao
148 :
149 96949 : energy%mulliken = 0.0_dp
150 :
151 96949 : 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 96949 : 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 96949 : 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 96949 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mo_derivs
193 : TYPE(cp_fm_type), POINTER :: mo_coeff
194 96949 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs, smat
195 96949 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
196 :
197 96949 : NULLIFY (mo_array, mo_coeff, mo_derivs)
198 :
199 96949 : 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 96949 : energy%s2_restraint = 0.0_dp
228 : END IF
229 193898 : END SUBROUTINE qs_ks_s2_restraint
230 :
231 : END MODULE qs_ks_apply_restraints
|