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 QMMM Hamiltonian integral matrix <a|\sum_i q_i|b> for
10 : !> semi-empirical methods
11 : !> \author Teodoro Laino 04.2007 [created]
12 : ! **************************************************************************************************
13 : MODULE qmmm_se_energy
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind
16 : USE cell_types, ONLY: cell_type,&
17 : pbc
18 : USE cp_control_types, ONLY: dft_control_type
19 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
20 : dbcsr_get_block_p,&
21 : dbcsr_p_type,&
22 : dbcsr_set
23 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
24 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
25 : USE cp_log_handling, ONLY: cp_get_default_logger,&
26 : cp_logger_type
27 : USE cp_output_handling, ONLY: cp_p_file,&
28 : cp_print_key_finished_output,&
29 : cp_print_key_should_output,&
30 : cp_print_key_unit_nr
31 : USE input_constants, ONLY: &
32 : do_method_am1, do_method_mndo, do_method_mndod, do_method_pchg, do_method_pdg, &
33 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1, &
34 : do_qmmm_coulomb, do_qmmm_gauss, do_qmmm_none, do_qmmm_pcharge, do_qmmm_swave
35 : USE kinds, ONLY: dp
36 : USE message_passing, ONLY: mp_para_env_type
37 : USE multipole_types, ONLY: do_multipole_none
38 : USE particle_types, ONLY: particle_type
39 : USE qmmm_types_low, ONLY: qmmm_env_qm_type,&
40 : qmmm_pot_p_type,&
41 : qmmm_pot_type
42 : USE qmmm_util, ONLY: spherical_cutoff_factor
43 : USE qs_energy_types, ONLY: qs_energy_type
44 : USE qs_environment_types, ONLY: get_qs_env,&
45 : qs_environment_type
46 : USE qs_kind_types, ONLY: get_qs_kind,&
47 : qs_kind_type
48 : USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
49 : USE qs_ks_types, ONLY: qs_ks_env_type,&
50 : set_ks_env
51 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
52 : USE qs_neighbor_lists, ONLY: build_qs_neighbor_lists
53 : USE qs_overlap, ONLY: build_overlap_matrix
54 : USE semi_empirical_int_arrays, ONLY: se_orbital_pointer
55 : USE semi_empirical_integrals, ONLY: corecore,&
56 : rotnuc
57 : USE semi_empirical_types, ONLY: get_se_param,&
58 : se_int_control_type,&
59 : se_taper_type,&
60 : semi_empirical_create,&
61 : semi_empirical_release,&
62 : semi_empirical_type,&
63 : setup_se_int_control_type
64 : USE semi_empirical_utils, ONLY: get_se_type,&
65 : se_param_set_default
66 : #include "./base/base_uses.f90"
67 :
68 : IMPLICIT NONE
69 :
70 : PRIVATE
71 :
72 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_se_energy'
73 :
74 : PUBLIC :: build_se_qmmm_matrix
75 :
76 : CONTAINS
77 :
78 : ! **************************************************************************************************
79 : !> \brief Constructs the 1-el semi-empirical hamiltonian
80 : !> \param qs_env ...
81 : !> \param qmmm_env ...
82 : !> \param particles_mm ...
83 : !> \param mm_cell ...
84 : !> \param para_env ...
85 : !> \author Teodoro Laino 04.2007 [created]
86 : ! **************************************************************************************************
87 1446 : SUBROUTINE build_se_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
88 :
89 : TYPE(qs_environment_type), POINTER :: qs_env
90 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
91 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
92 : TYPE(cell_type), POINTER :: mm_cell
93 : TYPE(mp_para_env_type), POINTER :: para_env
94 :
95 : CHARACTER(len=*), PARAMETER :: routineN = 'build_se_qmmm_matrix'
96 :
97 : INTEGER :: handle, i, iatom, ikind, itype, iw, &
98 : natom, natorb_a, nkind
99 1446 : INTEGER, DIMENSION(:), POINTER :: list
100 : LOGICAL :: anag, defined, found
101 : REAL(KIND=dp) :: enuclear
102 1446 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block_a
103 1446 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
104 : TYPE(cp_logger_type), POINTER :: logger
105 1446 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
106 : TYPE(dft_control_type), POINTER :: dft_control
107 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
108 1446 : POINTER :: sab_orb
109 1446 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
110 : TYPE(qs_energy_type), POINTER :: energy
111 1446 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
112 : TYPE(qs_ks_env_type), POINTER :: ks_env
113 : TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
114 : TYPE(se_int_control_type) :: se_int_control
115 : TYPE(se_taper_type), POINTER :: se_taper
116 : TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm
117 :
118 1446 : CALL timeset(routineN, handle)
119 1446 : NULLIFY (logger)
120 1446 : logger => cp_get_default_logger()
121 :
122 1446 : NULLIFY (matrix_s, atomic_kind_set, qs_kind_set, energy)
123 1446 : NULLIFY (se_kind_a, se_kind_mm, se_taper, particles_qm, ks_env, sab_orb)
124 1446 : CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
125 : CALL get_qs_env(qs_env, &
126 : ks_env=ks_env, &
127 : matrix_s=matrix_s, &
128 : energy=energy, &
129 1446 : sab_orb=sab_orb)
130 :
131 : CALL build_overlap_matrix(ks_env, matrix_s=matrix_s, &
132 : matrix_name="OVERLAP", &
133 : basis_type_a="ORB", &
134 : basis_type_b="ORB", &
135 1446 : sab_nl=sab_orb)
136 :
137 1446 : CALL set_ks_env(ks_env, matrix_s=matrix_s)
138 : CALL get_qs_env(qs_env=qs_env, &
139 : se_taper=se_taper, &
140 : atomic_kind_set=atomic_kind_set, &
141 : qs_kind_set=qs_kind_set, &
142 : ks_qmmm_env=ks_qmmm_env_loc, &
143 : dft_control=dft_control, &
144 1446 : particle_set=particles_qm)
145 :
146 1446 : SELECT CASE (dft_control%qs_control%method_id)
147 : CASE (do_method_am1, do_method_rm1, do_method_mndo, do_method_pdg, &
148 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
149 : ! Go on with the calculation..
150 : CASE DEFAULT
151 : ! Otherwise stop..
152 1446 : CPABORT("Method not available")
153 : END SELECT
154 1446 : anag = dft_control%qs_control%se_control%analytical_gradients
155 : ! Setup type for SE integral control
156 : CALL setup_se_int_control_type( &
157 : se_int_control, shortrange=.FALSE., do_ewald_r3=.FALSE., &
158 : do_ewald_gks=.FALSE., integral_screening=dft_control%qs_control%se_control%integral_screening, &
159 1446 : max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.)
160 :
161 : ! Allocate the core Hamiltonian matrix
162 1446 : CALL dbcsr_allocate_matrix_set(ks_qmmm_env_loc%matrix_h, 1)
163 1446 : ALLOCATE (ks_qmmm_env_loc%matrix_h(1)%matrix)
164 :
165 : CALL dbcsr_copy(ks_qmmm_env_loc%matrix_h(1)%matrix, matrix_s(1)%matrix, &
166 1446 : name="QMMM HAMILTONIAN MATRIX")
167 1446 : CALL dbcsr_set(ks_qmmm_env_loc%matrix_h(1)%matrix, 0.0_dp)
168 :
169 2382 : SELECT CASE (qmmm_env%qmmm_coupl_type)
170 : CASE (do_qmmm_coulomb, do_qmmm_gauss, do_qmmm_swave, do_qmmm_pcharge)
171 : ! Create a fake semi-empirical type to handle the classical atom
172 936 : CALL semi_empirical_create(se_kind_mm)
173 936 : CALL se_param_set_default(se_kind_mm, 0, do_method_pchg)
174 936 : itype = get_se_type(se_kind_mm%typ)
175 936 : nkind = SIZE(atomic_kind_set)
176 936 : enuclear = 0.0_dp
177 2860 : Kinds: DO ikind = 1, nkind
178 1924 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=list)
179 1924 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
180 : CALL get_se_param(se_kind_a, &
181 : defined=defined, &
182 1924 : natorb=natorb_a)
183 1924 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
184 11280 : Atoms: DO i = 1, SIZE(list)
185 4572 : iatom = list(i)
186 : ! Give back block
187 4572 : NULLIFY (h_block_a)
188 : CALL dbcsr_get_block_p(matrix=ks_qmmm_env_loc%matrix_h(1)%matrix, &
189 4572 : row=iatom, col=iatom, BLOCK=h_block_a, found=found)
190 :
191 6496 : IF (ASSOCIATED(h_block_a)) THEN
192 22410 : h_block_a = 0.0_dp
193 : ! Real QM/MM computation
194 : CALL build_se_qmmm_matrix_low(h_block_a, &
195 : se_kind_a, &
196 : se_kind_mm, &
197 : qmmm_env%Potentials, &
198 : particles_mm, &
199 : qmmm_env%mm_atom_chrg, &
200 : qmmm_env%mm_atom_index, &
201 : mm_cell, &
202 : iatom, &
203 : enuclear, &
204 : itype, &
205 : se_taper, &
206 : se_int_control, &
207 : anag, &
208 : qmmm_env%spherical_cutoff, &
209 2286 : particles_qm)
210 : ! Possibly added charges
211 2286 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
212 : CALL build_se_qmmm_matrix_low(h_block_a, &
213 : se_kind_a, &
214 : se_kind_mm, &
215 : qmmm_env%added_charges%potentials, &
216 : qmmm_env%added_charges%added_particles, &
217 : qmmm_env%added_charges%mm_atom_chrg, &
218 : qmmm_env%added_charges%mm_atom_index, &
219 : mm_cell, &
220 : iatom, &
221 : enuclear, &
222 : itype, &
223 : se_taper, &
224 : se_int_control, &
225 : anag, &
226 : qmmm_env%spherical_cutoff, &
227 0 : particles_qm)
228 : END IF
229 : END IF
230 : END DO Atoms
231 : END DO Kinds
232 936 : CALL para_env%sum(enuclear)
233 936 : energy%qmmm_nu = enuclear
234 936 : CALL semi_empirical_release(se_kind_mm)
235 : CASE (do_qmmm_none)
236 : ! Zero Matrix
237 1446 : CALL dbcsr_set(ks_qmmm_env_loc%matrix_h(1)%matrix, 0.0_dp)
238 : END SELECT
239 1446 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
240 : qs_env%input, "QMMM%PRINT%QMMM_MATRIX"), cp_p_file)) THEN
241 : iw = cp_print_key_unit_nr(logger, qs_env%input, "QMMM%PRINT%QMMM_MATRIX", &
242 196 : extension=".Log")
243 : CALL cp_dbcsr_write_sparse_matrix(ks_qmmm_env_loc%matrix_h(1)%matrix, 4, 6, qs_env, para_env, &
244 196 : scale=1.0_dp, output_unit=iw)
245 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
246 196 : "QMMM%PRINT%QMMM_MATRIX")
247 : END IF
248 :
249 1446 : CALL timestop(handle)
250 :
251 1446 : END SUBROUTINE build_se_qmmm_matrix
252 :
253 : ! **************************************************************************************************
254 : !> \brief Low Level : Constructs the 1-el semi-empirical hamiltonian block
255 : !> \param h_block_a ...
256 : !> \param se_kind_a ...
257 : !> \param se_kind_mm ...
258 : !> \param potentials ...
259 : !> \param particles_mm ...
260 : !> \param mm_charges ...
261 : !> \param mm_atom_index ...
262 : !> \param mm_cell ...
263 : !> \param IndQM ...
264 : !> \param enuclear ...
265 : !> \param itype ...
266 : !> \param se_taper ...
267 : !> \param se_int_control ...
268 : !> \param anag ...
269 : !> \param qmmm_spherical_cutoff ...
270 : !> \param particles_qm ...
271 : !> \author Teodoro Laino 04.2007 [created]
272 : ! **************************************************************************************************
273 2286 : SUBROUTINE build_se_qmmm_matrix_low(h_block_a, se_kind_a, se_kind_mm, potentials, &
274 : particles_mm, mm_charges, mm_atom_index, &
275 : mm_cell, IndQM, enuclear, itype, se_taper, se_int_control, anag, &
276 : qmmm_spherical_cutoff, particles_qm)
277 :
278 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block_a
279 : TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm
280 : TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
281 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
282 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges
283 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
284 : TYPE(cell_type), POINTER :: mm_cell
285 : INTEGER, INTENT(IN) :: IndQM
286 : REAL(KIND=dp), INTENT(INOUT) :: enuclear
287 : INTEGER, INTENT(IN) :: itype
288 : TYPE(se_taper_type), POINTER :: se_taper
289 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
290 : LOGICAL, INTENT(IN) :: anag
291 : REAL(KIND=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
292 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
293 :
294 : CHARACTER(len=*), PARAMETER :: routineN = 'build_se_qmmm_matrix_low'
295 :
296 : INTEGER :: handle, i1, i1L, i2, Imm, Imp, IndMM, &
297 : Ipot, j1, j1L
298 : REAL(KIND=dp) :: enuc, rt1, rt2, rt3, sph_chrg_factor
299 : REAL(KIND=dp), DIMENSION(3) :: r_pbc, rij
300 : REAL(KIND=dp), DIMENSION(45) :: e1b
301 : TYPE(qmmm_pot_type), POINTER :: Pot
302 :
303 2286 : CALL timeset(routineN, handle)
304 : ! Loop Over MM atoms
305 : ! Loop over Pot stores atoms with the same charge
306 6246 : MainLoopPot: DO Ipot = 1, SIZE(Potentials)
307 3960 : Pot => Potentials(Ipot)%Pot
308 : ! Loop over atoms belonging to this type
309 8437473 : LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index)
310 8431227 : Imm = Pot%mm_atom_index(Imp)
311 8431227 : IndMM = mm_atom_index(Imm)
312 33724908 : r_pbc = pbc(particles_mm(IndMM)%r - particles_qm(IndQM)%r, mm_cell)
313 8431227 : rt1 = r_pbc(1)
314 8431227 : rt2 = r_pbc(2)
315 8431227 : rt3 = r_pbc(3)
316 33724908 : rij = (/rt1, rt2, rt3/)
317 8431227 : se_kind_mm%zeff = mm_charges(Imm)
318 : ! Computes the screening factor for the spherical cutoff (if defined)
319 8431227 : IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
320 910272 : CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
321 910272 : se_kind_mm%zeff = se_kind_mm%zeff*sph_chrg_factor
322 : END IF
323 8431227 : IF (ABS(se_kind_mm%zeff) <= EPSILON(0.0_dp)) CYCLE
324 : CALL rotnuc(se_kind_a, se_kind_mm, rij, itype=itype, e1b=e1b, anag=anag, &
325 7672448 : se_int_control=se_int_control, se_taper=se_taper)
326 : CALL corecore(se_kind_a, se_kind_mm, rij, itype=itype, enuc=enuc, anag=anag, &
327 7672448 : se_int_control=se_int_control, se_taper=se_taper)
328 7672448 : enuclear = enuclear + enuc
329 : ! Contribution to the iatom block
330 : ! Computation of the QMMM core matrix
331 7672448 : i2 = 0
332 25592743 : DO i1L = 1, se_kind_a%natorb
333 17916335 : i1 = se_orbital_pointer(i1L)
334 38404109 : DO j1L = 1, i1L - 1
335 20487774 : j1 = se_orbital_pointer(j1L)
336 20487774 : i2 = i2 + 1
337 20487774 : h_block_a(i1, j1) = h_block_a(i1, j1) + e1b(i2)
338 38404109 : h_block_a(j1, i1) = h_block_a(i1, j1)
339 : END DO
340 17916335 : j1 = se_orbital_pointer(j1L)
341 17916335 : i2 = i2 + 1
342 26347562 : h_block_a(i1, j1) = h_block_a(i1, j1) + e1b(i2)
343 : END DO
344 : END DO LoopMM
345 : END DO MainLoopPot
346 2286 : CALL timestop(handle)
347 2286 : END SUBROUTINE build_se_qmmm_matrix_low
348 :
349 : END MODULE qmmm_se_energy
|
