(git:98357aa)
Loading...
Searching...
No Matches
qs_ks_apply_restraints.F
Go to the documentation of this file.
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! **************************************************************************************************
13 USE cp_dbcsr_api, ONLY: dbcsr_p_type
16 USE cp_fm_types, ONLY: cp_fm_create,&
21 USE kinds, ONLY: dp
24 USE pw_methods, ONLY: pw_scale
32 USE qs_mo_types, ONLY: get_mo_set,&
34 USE qs_rho_types, ONLY: qs_rho_get,&
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
47 PUBLIC :: qs_ks_cdft_constraint
48
49CONTAINS
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 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 NULLIFY (dft_control)
70 CALL get_qs_env(qs_env, dft_control=dft_control)
71 IF (dft_control%qs_control%cdft) THEN
72 cdft_control => dft_control%qs_control%cdft_control
73 ! Test no k-points
74 CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
75 IF (do_kpoints) cpabort("CDFT constraints with k-points not supported.")
76
77 SELECT CASE (cdft_control%type)
79 IF (cdft_control%need_pot) THEN
80 ! First SCF iteraration => allocate storage
81 DO igroup = 1, SIZE(cdft_control%group)
82 ALLOCATE (cdft_control%group(igroup)%weight)
83 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 .AND. dft_control%nspins == 1) THEN
87 CALL cp_abort(__location__, &
88 "Spin constraints require a spin polarized calculation.")
89 END IF
90 END DO
91 IF (cdft_control%atomic_charges) THEN
92 IF (.NOT. ASSOCIATED(cdft_control%charge)) THEN
93 ALLOCATE (cdft_control%charge(cdft_control%natoms))
94 END IF
95 DO iatom = 1, cdft_control%natoms
96 CALL auxbas_pw_pool%create_pw(cdft_control%charge(iatom))
97 END DO
98 END IF
99 ! Another sanity check
100 CALL get_qs_env(qs_env, natom=natom)
101 IF (natom < cdft_control%natoms) THEN
102 CALL cp_abort(__location__, &
103 "The number of constraint atoms exceeds the total number of atoms.")
104 END IF
105 ELSE
106 DO igroup = 1, SIZE(cdft_control%group)
107 inv_vol = 1.0_dp/cdft_control%group(igroup)%weight%pw_grid%dvol
108 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 IF (cdft_control%type == outer_scf_becke_constraint) THEN
113 CALL becke_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
114 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
115 CALL hirshfeld_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
116 END IF
117 DO igroup = 1, SIZE(cdft_control%group)
118 CALL pw_scale(cdft_control%group(igroup)%weight, cdft_control%group(igroup)%weight%pw_grid%dvol)
119 END DO
120 IF (cdft_control%need_pot) cdft_control%need_pot = .false.
121 CASE DEFAULT
122 cpabort("Unknown constraint type.")
123 END SELECT
124 END IF
125
126 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, rho_ao
151
152 energy%mulliken = 0.0_dp
153
154 IF (dft_control%qs_control%mulliken_restraint) THEN
155
156 ! Test no k-points
157 cpassert(SIZE(matrix_s, 2) == 1)
158
159 CALL qs_rho_get(rho, rho_ao=rho_ao)
160
161 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 order_p=mulliken_order_p)
165 ELSE
166 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 ks_matrix=ksmat, order_p=mulliken_order_p)
170 END IF
171
172 END IF
173
174 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 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mo_derivs
196 TYPE(cp_fm_type), POINTER :: mo_coeff
197 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs, smat
198 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
199
200 NULLIFY (mo_array, mo_coeff, mo_derivs)
201
202 IF (dft_control%qs_control%s2_restraint) THEN
203 ! Test no k-points
204 cpassert(SIZE(matrix_s, 2) == 1)
205 ! adds s2_restraint energy and orbital derivatives
206 cpassert(dft_control%nspins == 2)
207 cpassert(qs_env%requires_mo_derivs)
208 ! forces are not implemented (not difficult, but ... )
209 cpassert(.NOT. calculate_forces)
210 mark_used(calculate_forces)
211 CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
212
213 ALLOCATE (fm_mo_derivs(SIZE(mo_derivs, 1))) !fm->dbcsr
214 DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
215 CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff) !fm->dbcsr
216 CALL cp_fm_create(fm_mo_derivs(i), mo_coeff%matrix_struct) !fm->dbcsr
217 CALL copy_dbcsr_to_fm(mo_derivs(i)%matrix, fm_mo_derivs(i)) !fm->dbcsr
218 END DO !fm->dbcsr
219
220 smat => matrix_s(:, 1)
221 CALL s2_restraint(mo_array, smat, fm_mo_derivs, energy%s2_restraint, &
222 dft_control%qs_control%s2_restraint_control, just_energy)
223
224 DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
225 CALL copy_fm_to_dbcsr(fm_mo_derivs(i), mo_derivs(i)%matrix) !fm->dbcsr
226 END DO !fm->dbcsr
227 DEALLOCATE (fm_mo_derivs) !fm->dbcsr
228
229 ELSE
230 energy%s2_restraint = 0.0_dp
231 END IF
232 END SUBROUTINE qs_ks_s2_restraint
233
234END MODULE qs_ks_apply_restraints
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public cdft_charge_constraint
integer, parameter, public outer_scf_becke_constraint
integer, parameter, public outer_scf_hirshfeld_constraint
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Interface to the message passing library MPI.
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition mulliken.F:13
subroutine, public mulliken_restraint(mulliken_restraint_control, para_env, s_matrix, p_matrix, energy, order_p, ks_matrix, w_matrix)
computes the energy and density matrix derivate of a constraint on the mulliken charges
Definition mulliken.F:73
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Subroutines for building CDFT constraints.
subroutine, public becke_constraint(qs_env, calc_pot, calculate_forces)
Driver routine for calculating a Becke constraint.
subroutine, public hirshfeld_constraint(qs_env, calc_pot, calculate_forces)
Driver routine for calculating a Hirshfeld constraint.
Defines CDFT control structures.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Set of routines to apply restraints to the KS hamiltonian.
subroutine, public qs_ks_s2_restraint(dft_control, qs_env, matrix_s, energy, calculate_forces, just_energy)
...
subroutine, public qs_ks_mulliken_restraint(energy, dft_control, just_energy, para_env, ks_matrix, matrix_s, rho, mulliken_order_p)
...
subroutine, public qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
Apply a CDFT constraint.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Methods related to (\cal S)^2 (i.e. spin)
subroutine, public s2_restraint(mos, matrix_s, mo_derivs, energy, s2_restraint_control, just_energy)
restrains/constrains the value of s2 in a calculation
represent a full matrix
stores all the informations relevant to an mpi environment
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
keeps the density in various representations, keeping track of which ones are valid.