Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of the Fock matrix for SE methods
10 : !> \author JGH and TLAINO
11 : !> \par History
12 : !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : d-orbitals
13 : !> Teodoro Laino (09.2008) [tlaino] - University of Zurich : Speed-up
14 : !> Teodoro Laino (09.2008) [tlaino] - University of Zurich : Periodic SE
15 : !> Teodoro Laino (05.2009) [tlaino] - Split and module reorganization
16 : ! **************************************************************************************************
17 : MODULE se_fock_matrix
18 : USE atomic_kind_types, ONLY: atomic_kind_type
19 : USE atprop_types, ONLY: atprop_array_init,&
20 : atprop_type
21 : USE cp_control_types, ONLY: dft_control_type,&
22 : semi_empirical_control_type
23 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
24 : dbcsr_copy,&
25 : dbcsr_get_info,&
26 : dbcsr_multiply,&
27 : dbcsr_p_type,&
28 : dbcsr_type
29 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
30 : USE cp_log_handling, ONLY: cp_get_default_logger,&
31 : cp_logger_type
32 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
33 : cp_print_key_unit_nr
34 : USE input_constants, ONLY: &
35 : do_method_am1, do_method_mndo, do_method_mndod, do_method_pdg, do_method_pm3, &
36 : do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1
37 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
38 : section_vals_type
39 : USE kinds, ONLY: dp
40 : USE message_passing, ONLY: mp_para_env_type
41 : USE particle_types, ONLY: particle_type
42 : USE qs_energy_types, ONLY: qs_energy_type
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_ks_types, ONLY: qs_ks_env_type
46 : USE qs_mo_types, ONLY: get_mo_set,&
47 : mo_set_type
48 : USE qs_rho_types, ONLY: qs_rho_get,&
49 : qs_rho_type
50 : USE se_fock_matrix_coulomb, ONLY: build_fock_matrix_coul_lr_r3,&
51 : build_fock_matrix_coulomb,&
52 : build_fock_matrix_coulomb_lr
53 : USE se_fock_matrix_dbg, ONLY: dbg_energy_coulomb_lr
54 : USE se_fock_matrix_exchange, ONLY: build_fock_matrix_exchange
55 : USE semi_empirical_store_int_types, ONLY: semi_empirical_si_finalize,&
56 : semi_empirical_si_initialize,&
57 : semi_empirical_si_type
58 : #include "./base/base_uses.f90"
59 :
60 : IMPLICIT NONE
61 : PRIVATE
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_fock_matrix'
64 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
65 : LOGICAL, PARAMETER, PRIVATE :: debug_energy_coulomb_lr = .FALSE.
66 :
67 : PUBLIC :: build_se_fock_matrix
68 :
69 : CONTAINS
70 :
71 : ! **************************************************************************************************
72 : !> \brief Construction of the Fock matrix for NDDO methods
73 : !> \param qs_env ...
74 : !> \param calculate_forces ...
75 : !> \param just_energy ...
76 : !> \par History
77 : !> - Teodoro Laino [tlaino] (05.2009) - Split and module reorganization
78 : !> \author JGH
79 : ! **************************************************************************************************
80 41094 : SUBROUTINE build_se_fock_matrix(qs_env, calculate_forces, just_energy)
81 : TYPE(qs_environment_type), POINTER :: qs_env
82 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
83 :
84 : CHARACTER(len=*), PARAMETER :: routineN = 'build_se_fock_matrix'
85 :
86 : INTEGER :: handle, ispin, natom, ncol_global, &
87 : nspins, output_unit
88 : LOGICAL :: s_mstruct_changed
89 : REAL(KIND=dp) :: ecoul, qmmm_el
90 41094 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
91 41094 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
92 : TYPE(atprop_type), POINTER :: atprop
93 : TYPE(cp_logger_type), POINTER :: logger
94 41094 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_h, matrix_p, mo_derivs
95 : TYPE(dbcsr_type), POINTER :: mo_coeff
96 : TYPE(dft_control_type), POINTER :: dft_control
97 41094 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
98 : TYPE(mp_para_env_type), POINTER :: para_env
99 41094 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
100 : TYPE(qs_energy_type), POINTER :: energy
101 : TYPE(qs_ks_env_type), POINTER :: ks_env
102 : TYPE(qs_rho_type), POINTER :: rho
103 : TYPE(section_vals_type), POINTER :: scf_section
104 : TYPE(semi_empirical_control_type), POINTER :: se_control
105 : TYPE(semi_empirical_si_type), POINTER :: store_int_env
106 :
107 41094 : CALL timeset(routineN, handle)
108 41094 : NULLIFY (matrix_h, dft_control, logger, scf_section, store_int_env, se_control)
109 41094 : NULLIFY (atomic_kind_set, atprop)
110 41094 : NULLIFY (ks_env, ks_matrix, rho, energy)
111 41094 : logger => cp_get_default_logger()
112 41094 : CPASSERT(ASSOCIATED(qs_env))
113 :
114 : CALL get_qs_env(qs_env, &
115 : dft_control=dft_control, &
116 : matrix_h=matrix_h, &
117 : para_env=para_env, &
118 : se_store_int_env=store_int_env, &
119 : atprop=atprop, &
120 : atomic_kind_set=atomic_kind_set, &
121 : s_mstruct_changed=s_mstruct_changed, &
122 : ks_env=ks_env, &
123 : matrix_ks=ks_matrix, &
124 : rho=rho, &
125 41094 : energy=energy)
126 :
127 41094 : SELECT CASE (dft_control%qs_control%method_id)
128 : CASE DEFAULT
129 : ! Abort if the parameterization is an unknown one..
130 0 : CPABORT("Fock Matrix not available for the chosen parameterization! ")
131 :
132 : CASE (do_method_am1, do_method_rm1, do_method_mndo, do_method_pdg, &
133 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
134 :
135 : ! Check for properly allocation of Matrixes
136 41094 : nspins = dft_control%nspins
137 41094 : CPASSERT(((nspins >= 1) .AND. (nspins <= 2)))
138 41094 : CPASSERT(ASSOCIATED(matrix_h))
139 41094 : CPASSERT(ASSOCIATED(rho))
140 41094 : CPASSERT(SIZE(ks_matrix) > 0)
141 :
142 41094 : se_control => dft_control%qs_control%se_control
143 41094 : scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
144 41094 : CALL qs_rho_get(rho, rho_ao=matrix_p)
145 :
146 41094 : energy%qmmm_el = 0.0_dp
147 41094 : energy%total = 0.0_dp
148 :
149 84878 : DO ispin = 1, nspins
150 : ! Copy the core matrix into the fock matrix
151 84878 : CALL dbcsr_copy(ks_matrix(ispin)%matrix, matrix_h(1)%matrix)
152 : END DO
153 :
154 : ! WRITE ( *, * ) 'KS_ENV%s_mstruct_changed', ks_env%s_mstruct_changed
155 41094 : IF (atprop%energy) THEN
156 102 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
157 102 : natom = SIZE(particle_set)
158 102 : CALL atprop_array_init(atprop%atecoul, natom)
159 : END IF
160 :
161 : ! Compute Exchange and Coulomb terms
162 41094 : CALL semi_empirical_si_initialize(store_int_env, s_mstruct_changed)
163 : CALL build_fock_matrix_exchange(qs_env, ks_matrix, matrix_p, calculate_forces, &
164 41094 : store_int_env)
165 : CALL build_fock_matrix_coulomb(qs_env, ks_matrix, matrix_p, energy, calculate_forces, &
166 41094 : store_int_env)
167 :
168 : ! Debug statements for Long-Range
169 : IF (debug_energy_coulomb_lr .AND. se_control%do_ewald) THEN
170 : CALL dbg_energy_coulomb_lr(energy, ks_matrix, nspins, qs_env, matrix_p, &
171 : calculate_forces, store_int_env)
172 : END IF
173 :
174 : ! Long Range Electrostatic
175 41094 : IF (se_control%do_ewald) THEN
176 : ! Evaluate Coulomb Long-Range
177 : CALL build_fock_matrix_coulomb_lr(qs_env, ks_matrix, matrix_p, energy, calculate_forces, &
178 1792 : store_int_env)
179 :
180 : ! Possibly handle the slowly convergent term 1/R^3
181 1792 : IF (se_control%do_ewald_r3) THEN
182 : CALL build_fock_matrix_coul_lr_r3(qs_env, ks_matrix, matrix_p, energy, &
183 0 : calculate_forces)
184 : END IF
185 : END IF
186 41094 : CALL semi_empirical_si_finalize(store_int_env, s_mstruct_changed)
187 :
188 41094 : IF (atprop%energy) THEN
189 714 : atprop%atecoul = 0.5_dp*atprop%atecoul
190 : END IF
191 :
192 : ! Compute the Hartree energy
193 : ! NOTE: If we are performing SCP-NDDO, ks_matrix contains coulomb piece from SCP.
194 84878 : DO ispin = 1, nspins
195 43784 : CALL dbcsr_dot(ks_matrix(ispin)%matrix, matrix_p(ispin)%matrix, ecoul)
196 84878 : energy%hartree = energy%hartree + ecoul
197 : END DO
198 : ! WRITE ( *, * ) 'AFTER Hartree', ecoul, energy%hartree
199 :
200 : ! CALL build_fock_matrix_ph(qs_env,ks_matrix)
201 : ! QM/MM
202 41094 : IF (qs_env%qmmm) THEN
203 23944 : DO ispin = 1, nspins
204 : ! If QM/MM sumup the 1el Hamiltonian
205 : CALL dbcsr_add(ks_matrix(ispin)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
206 11972 : 1.0_dp, 1.0_dp)
207 : ! Compute QM/MM Energy
208 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
209 11972 : matrix_p(ispin)%matrix, qmmm_el)
210 23944 : energy%qmmm_el = energy%qmmm_el + qmmm_el
211 : END DO
212 : END IF
213 :
214 : ! WRITE ( *, * ) ' before TOTAL', energy%total
215 : ! Collect all the energy terms
216 41094 : energy%mulliken = 0.0_dp
217 41094 : energy%exc = 0.0_dp
218 : energy%total = energy%total + energy%core + &
219 : energy%core_overlap + &
220 : 0.5_dp*energy%hartree + &
221 : energy%qmmm_el + &
222 : energy%dispersion + &
223 41094 : energy%mulliken
224 : ! WRITE ( *, * ) ' AFTER TOTAL', energy%total
225 :
226 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", &
227 41094 : extension=".scfLog")
228 :
229 41094 : IF (output_unit > 0) THEN
230 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T60,F20.10))") &
231 31 : "Core Hamiltonian energy: ", energy%core, &
232 62 : "Two-electron integral energy: ", energy%hartree
233 31 : IF (qs_env%qmmm) THEN
234 : WRITE (UNIT=output_unit, FMT="(T3,A,T60,F20.10)") &
235 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
236 : END IF
237 : END IF
238 :
239 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
240 41094 : "PRINT%DETAILED_ENERGY")
241 :
242 : ! Here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C (plus occupation numbers)
243 82188 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
244 14946 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
245 30564 : DO ispin = 1, SIZE(mo_derivs)
246 : CALL get_mo_set(mo_set=mo_array(ispin), &
247 15618 : mo_coeff_b=mo_coeff, occupation_numbers=occupation_numbers)
248 15618 : IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
249 0 : CPABORT("")
250 : END IF
251 15618 : CALL dbcsr_get_info(mo_coeff, nfullcols_total=ncol_global)
252 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin)%matrix, mo_coeff, &
253 30564 : 0.0_dp, mo_derivs(ispin)%matrix)
254 : END DO
255 : END IF
256 :
257 : END SELECT
258 :
259 41094 : CALL timestop(handle)
260 :
261 41094 : END SUBROUTINE build_se_fock_matrix
262 :
263 : END MODULE se_fock_matrix
264 :
|
