(git:374b731)
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-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! **************************************************************************************************
15 USE cp_fm_types, ONLY: cp_fm_create,&
17 USE dbcsr_api, ONLY: dbcsr_p_type
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) &
87 CALL cp_abort(__location__, &
88 "Spin constraints require a spin polarized calculation.")
89 END DO
90 IF (cdft_control%atomic_charges) THEN
91 IF (.NOT. ASSOCIATED(cdft_control%charge)) &
92 ALLOCATE (cdft_control%charge(cdft_control%natoms))
93 DO iatom = 1, cdft_control%natoms
94 CALL auxbas_pw_pool%create_pw(cdft_control%charge(iatom))
95 END DO
96 END IF
97 ! Another sanity check
98 CALL get_qs_env(qs_env, natom=natom)
99 IF (natom < cdft_control%natoms) &
100 CALL cp_abort(__location__, &
101 "The number of constraint atoms exceeds the total number of atoms.")
102 ELSE
103 DO igroup = 1, SIZE(cdft_control%group)
104 inv_vol = 1.0_dp/cdft_control%group(igroup)%weight%pw_grid%dvol
105 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 IF (cdft_control%type == outer_scf_becke_constraint) THEN
110 CALL becke_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
111 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
112 CALL hirshfeld_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
113 END IF
114 DO igroup = 1, SIZE(cdft_control%group)
115 CALL pw_scale(cdft_control%group(igroup)%weight, cdft_control%group(igroup)%weight%pw_grid%dvol)
116 END DO
117 IF (cdft_control%need_pot) cdft_control%need_pot = .false.
118 CASE DEFAULT
119 cpabort("Unknown constraint type.")
120 END SELECT
121 END IF
122
123 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, rho_ao
148
149 energy%mulliken = 0.0_dp
150
151 IF (dft_control%qs_control%mulliken_restraint) THEN
152
153 ! Test no k-points
154 cpassert(SIZE(matrix_s, 2) == 1)
155
156 CALL qs_rho_get(rho, rho_ao=rho_ao)
157
158 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 order_p=mulliken_order_p)
162 ELSE
163 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 ks_matrix=ksmat, order_p=mulliken_order_p)
167 END IF
168
169 END IF
170
171 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 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mo_derivs
193 TYPE(cp_fm_type), POINTER :: mo_coeff
194 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs, smat
195 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
196
197 NULLIFY (mo_array, mo_coeff, mo_derivs)
198
199 IF (dft_control%qs_control%s2_restraint) THEN
200 ! Test no k-points
201 cpassert(SIZE(matrix_s, 2) == 1)
202 ! adds s2_restraint energy and orbital derivatives
203 cpassert(dft_control%nspins == 2)
204 cpassert(qs_env%requires_mo_derivs)
205 ! forces are not implemented (not difficult, but ... )
206 cpassert(.NOT. calculate_forces)
207 mark_used(calculate_forces)
208 CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
209
210 ALLOCATE (fm_mo_derivs(SIZE(mo_derivs, 1))) !fm->dbcsr
211 DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
212 CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff) !fm->dbcsr
213 CALL cp_fm_create(fm_mo_derivs(i), mo_coeff%matrix_struct) !fm->dbcsr
214 CALL copy_dbcsr_to_fm(mo_derivs(i)%matrix, fm_mo_derivs(i)) !fm->dbcsr
215 END DO !fm->dbcsr
216
217 smat => matrix_s(:, 1)
218 CALL s2_restraint(mo_array, smat, fm_mo_derivs, energy%s2_restraint, &
219 dft_control%qs_control%s2_restraint_control, just_energy)
220
221 DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
222 CALL copy_fm_to_dbcsr(fm_mo_derivs(i), mo_derivs(i)%matrix) !fm->dbcsr
223 END DO !fm->dbcsr
224 DEALLOCATE (fm_mo_derivs) !fm->dbcsr
225
226 ELSE
227 energy%s2_restraint = 0.0_dp
228 END IF
229 END SUBROUTINE qs_ks_s2_restraint
230
231END 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, use_sp)
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
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.