Line data Source code
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 TB methods used with QMMM
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE qmmm_tb_methods
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind
15 : USE cell_types, ONLY: cell_type,&
16 : pbc
17 : USE cp_control_types, ONLY: dft_control_type,&
18 : dftb_control_type,&
19 : xtb_control_type
20 : USE cp_dbcsr_api, ONLY: &
21 : dbcsr_add, dbcsr_copy, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
22 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
23 : dbcsr_p_type, dbcsr_set
24 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
25 : dbcsr_deallocate_matrix_set
26 : USE ewald_environment_types, ONLY: ewald_env_create,&
27 : ewald_env_get,&
28 : ewald_env_release,&
29 : ewald_env_set,&
30 : ewald_environment_type,&
31 : read_ewald_section
32 : USE ewald_pw_types, ONLY: ewald_pw_create,&
33 : ewald_pw_release,&
34 : ewald_pw_type
35 : USE input_constants, ONLY: do_fist_pol_none
36 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
37 : section_vals_type
38 : USE kinds, ONLY: dp
39 : USE mathconstants, ONLY: rootpi
40 : USE message_passing, ONLY: mp_para_env_type
41 : USE mulliken, ONLY: mulliken_charges
42 : USE particle_types, ONLY: allocate_particle_set,&
43 : deallocate_particle_set,&
44 : particle_type
45 : USE pw_poisson_types, ONLY: do_ewald_ewald,&
46 : do_ewald_none,&
47 : do_ewald_pme,&
48 : do_ewald_spme
49 : USE qmmm_types_low, ONLY: qmmm_env_qm_type,&
50 : qmmm_pot_p_type,&
51 : qmmm_pot_type
52 : USE qmmm_util, ONLY: spherical_cutoff_factor
53 : USE qs_dftb_coulomb, ONLY: gamma_rab_sr
54 : USE qs_dftb_matrices, ONLY: build_dftb_overlap
55 : USE qs_dftb_types, ONLY: qs_dftb_atom_type
56 : USE qs_dftb_utils, ONLY: get_dftb_atom_param
57 : USE qs_environment_types, ONLY: get_qs_env,&
58 : qs_environment_type
59 : USE qs_kind_types, ONLY: get_qs_kind,&
60 : qs_kind_type
61 : USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
62 : USE qs_ks_types, ONLY: qs_ks_env_type
63 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
64 : USE qs_neighbor_lists, ONLY: build_qs_neighbor_lists
65 : USE qs_overlap, ONLY: build_overlap_matrix
66 : USE qs_rho_types, ONLY: qs_rho_get,&
67 : qs_rho_type
68 : USE spme, ONLY: spme_forces,&
69 : spme_potential
70 : USE xtb_types, ONLY: get_xtb_atom_param,&
71 : xtb_atom_type
72 : #include "./base/base_uses.f90"
73 :
74 : IMPLICIT NONE
75 :
76 : ! small real number
77 : REAL(dp), PARAMETER :: rtiny = 1.e-10_dp
78 : ! eta(0) for mm atoms and non-scc qm atoms
79 : REAL(dp), PARAMETER :: eta_mm = 0.47_dp
80 : ! step size for qmmm finite difference
81 : REAL(dp), PARAMETER :: ddrmm = 0.0001_dp
82 : INTEGER, PARAMETER :: pot_tb_nonperiodic = 0
83 : INTEGER, PARAMETER :: pot_tb_short_range = 1
84 : INTEGER, PARAMETER :: pot_ewald_short_range = 2
85 : INTEGER, PARAMETER :: pot_gauss_nonperiodic = 3
86 : INTEGER, PARAMETER :: pot_gauss_short_range = 4
87 :
88 : PRIVATE
89 :
90 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_tb_methods'
91 :
92 : PUBLIC :: build_tb_qmmm_matrix, build_tb_qmmm_matrix_gauss, build_tb_qmmm_matrix_zero, &
93 : build_tb_qmmm_matrix_pc, deriv_tb_qmmm_matrix, deriv_tb_qmmm_matrix_gauss, &
94 : deriv_tb_qmmm_matrix_pc
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 : !> \brief Constructs the 1-el DFTB hamiltonian
100 : !> \param qs_env ...
101 : !> \param qmmm_env ...
102 : !> \param particles_mm ...
103 : !> \param mm_cell ...
104 : !> \param para_env ...
105 : !> \author JGH 10.2014 [created]
106 : ! **************************************************************************************************
107 448 : SUBROUTINE build_tb_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
108 :
109 : TYPE(qs_environment_type), POINTER :: qs_env
110 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
111 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
112 : TYPE(cell_type), POINTER :: mm_cell
113 : TYPE(mp_para_env_type), POINTER :: para_env
114 :
115 : CHARACTER(len=*), PARAMETER :: routineN = 'build_tb_qmmm_matrix'
116 :
117 : INTEGER :: handle, i, iatom, ikind, jatom, natom, &
118 : natorb, nkind
119 448 : INTEGER, DIMENSION(:), POINTER :: list
120 : LOGICAL :: defined, do_dftb, do_xtb, found
121 : REAL(KIND=dp) :: pc_ener, zeff
122 : REAL(KIND=dp), DIMENSION(0:3) :: eta_a
123 448 : REAL(KIND=dp), DIMENSION(:), POINTER :: qpot
124 448 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hblock, sblock
125 448 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
126 : TYPE(dbcsr_iterator_type) :: iter
127 448 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_s
128 : TYPE(dft_control_type), POINTER :: dft_control
129 : TYPE(dftb_control_type), POINTER :: dftb_control
130 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
131 448 : POINTER :: sab_nl
132 448 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
133 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
134 448 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
135 : TYPE(qs_ks_env_type), POINTER :: ks_env
136 : TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
137 : TYPE(qs_rho_type), POINTER :: rho
138 : TYPE(xtb_atom_type), POINTER :: xtb_kind
139 : TYPE(xtb_control_type), POINTER :: xtb_control
140 :
141 448 : CALL timeset(routineN, handle)
142 :
143 : CALL get_qs_env(qs_env=qs_env, &
144 : dft_control=dft_control, &
145 : atomic_kind_set=atomic_kind_set, &
146 : particle_set=particles_qm, &
147 : qs_kind_set=qs_kind_set, &
148 : rho=rho, &
149 448 : natom=natom)
150 448 : dftb_control => dft_control%qs_control%dftb_control
151 448 : xtb_control => dft_control%qs_control%xtb_control
152 :
153 448 : IF (dft_control%qs_control%dftb) THEN
154 : do_dftb = .TRUE.
155 : do_xtb = .FALSE.
156 224 : ELSE IF (dft_control%qs_control%xtb) THEN
157 : do_dftb = .FALSE.
158 : do_xtb = .TRUE.
159 : ELSE
160 0 : CPABORT("TB method unknown")
161 : END IF
162 :
163 448 : CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
164 :
165 448 : NULLIFY (matrix_s)
166 448 : IF (do_dftb) THEN
167 224 : CALL build_dftb_overlap(qs_env, 0, matrix_s)
168 224 : ELSE IF (do_xtb) THEN
169 224 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
170 224 : CALL build_overlap_matrix(ks_env, matrix_s, basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
171 : END IF
172 :
173 1344 : ALLOCATE (qpot(natom))
174 1792 : qpot = 0.0_dp
175 448 : pc_ener = 0.0_dp
176 :
177 448 : nkind = SIZE(atomic_kind_set)
178 1344 : DO ikind = 1, nkind
179 896 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
180 896 : IF (do_dftb) THEN
181 448 : NULLIFY (dftb_kind)
182 448 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
183 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff, &
184 448 : defined=defined, eta=eta_a, natorb=natorb)
185 : ! use mm charge smearing for non-scc cases
186 448 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
187 448 : IF (.NOT. defined .OR. natorb < 1) CYCLE
188 448 : ELSE IF (do_xtb) THEN
189 448 : NULLIFY (xtb_kind)
190 448 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
191 448 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
192 448 : eta_a(0) = eta_mm
193 : END IF
194 2688 : DO i = 1, SIZE(list)
195 1344 : iatom = list(i)
196 : CALL build_mm_pot(qpot(iatom), pot_tb_nonperiodic, eta_a(0), &
197 : qmmm_env%Potentials, particles_mm, &
198 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, &
199 1344 : qmmm_env%spherical_cutoff, particles_qm)
200 : ! Possibly added charges
201 1344 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
202 : CALL build_mm_pot(qpot(iatom), pot_tb_nonperiodic, eta_a(0), &
203 : qmmm_env%added_charges%potentials, &
204 : qmmm_env%added_charges%added_particles, &
205 : qmmm_env%added_charges%mm_atom_chrg, &
206 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
207 : qmmm_env%spherical_cutoff, &
208 0 : particles_qm)
209 : END IF
210 2240 : pc_ener = pc_ener + qpot(iatom)*zeff
211 : END DO
212 : END DO
213 :
214 : ! Allocate the core Hamiltonian matrix
215 448 : CALL get_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env_loc)
216 448 : matrix_h => ks_qmmm_env_loc%matrix_h
217 448 : CALL dbcsr_allocate_matrix_set(matrix_h, 1)
218 448 : ALLOCATE (matrix_h(1)%matrix)
219 : CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, &
220 448 : name="QMMM HAMILTONIAN MATRIX")
221 448 : CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
222 :
223 448 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
224 1792 : DO WHILE (dbcsr_iterator_blocks_left(iter))
225 1344 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock)
226 1344 : NULLIFY (hblock)
227 : CALL dbcsr_get_block_p(matrix=matrix_h(1)%matrix, &
228 1344 : row=iatom, col=jatom, block=hblock, found=found)
229 1344 : CPASSERT(found)
230 26432 : hblock = hblock - 0.5_dp*sblock*(qpot(iatom) + qpot(jatom))
231 : END DO
232 448 : CALL dbcsr_iterator_stop(iter)
233 :
234 448 : ks_qmmm_env_loc%matrix_h => matrix_h
235 448 : ks_qmmm_env_loc%pc_ener = pc_ener
236 :
237 448 : DEALLOCATE (qpot)
238 :
239 448 : CALL dbcsr_deallocate_matrix_set(matrix_s)
240 :
241 448 : CALL timestop(handle)
242 :
243 896 : END SUBROUTINE build_tb_qmmm_matrix
244 :
245 : ! **************************************************************************************************
246 : !> \brief Constructs an empty 1-el DFTB hamiltonian
247 : !> \param qs_env ...
248 : !> \param para_env ...
249 : !> \author JGH 10.2014 [created]
250 : ! **************************************************************************************************
251 8 : SUBROUTINE build_tb_qmmm_matrix_zero(qs_env, para_env)
252 :
253 : TYPE(qs_environment_type), POINTER :: qs_env
254 : TYPE(mp_para_env_type), POINTER :: para_env
255 :
256 : CHARACTER(len=*), PARAMETER :: routineN = 'build_tb_qmmm_matrix_zero'
257 :
258 : INTEGER :: handle
259 : LOGICAL :: do_dftb, do_xtb
260 8 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_s
261 : TYPE(dft_control_type), POINTER :: dft_control
262 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
263 8 : POINTER :: sab_nl
264 : TYPE(qs_ks_env_type), POINTER :: ks_env
265 : TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
266 :
267 8 : CALL timeset(routineN, handle)
268 :
269 8 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
270 :
271 8 : IF (dft_control%qs_control%dftb) THEN
272 : do_dftb = .TRUE.
273 : do_xtb = .FALSE.
274 4 : ELSE IF (dft_control%qs_control%xtb) THEN
275 : do_dftb = .FALSE.
276 : do_xtb = .TRUE.
277 : ELSE
278 0 : CPABORT("TB method unknown")
279 : END IF
280 :
281 8 : CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
282 :
283 8 : NULLIFY (matrix_s)
284 8 : IF (do_dftb) THEN
285 4 : CALL build_dftb_overlap(qs_env, 0, matrix_s)
286 4 : ELSE IF (do_xtb) THEN
287 4 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
288 4 : CALL build_overlap_matrix(ks_env, matrix_s, basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
289 : END IF
290 :
291 : ! Allocate the core Hamiltonian matrix
292 8 : CALL get_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env_loc)
293 8 : matrix_h => ks_qmmm_env_loc%matrix_h
294 8 : CALL dbcsr_allocate_matrix_set(matrix_h, 1)
295 8 : ALLOCATE (matrix_h(1)%matrix)
296 : CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, &
297 8 : name="QMMM HAMILTONIAN MATRIX")
298 8 : CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
299 8 : ks_qmmm_env_loc%matrix_h => matrix_h
300 8 : ks_qmmm_env_loc%pc_ener = 0.0_dp
301 :
302 8 : CALL dbcsr_deallocate_matrix_set(matrix_s)
303 :
304 8 : CALL timestop(handle)
305 :
306 8 : END SUBROUTINE build_tb_qmmm_matrix_zero
307 :
308 : ! **************************************************************************************************
309 : !> \brief Constructs the 1-el DFTB hamiltonian
310 : !> \param qs_env ...
311 : !> \param qmmm_env ...
312 : !> \param particles_mm ...
313 : !> \param mm_cell ...
314 : !> \param para_env ...
315 : !> \author JGH 10.2014 [created]
316 : ! **************************************************************************************************
317 1116 : SUBROUTINE build_tb_qmmm_matrix_pc(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
318 :
319 : TYPE(qs_environment_type), POINTER :: qs_env
320 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
321 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
322 : TYPE(cell_type), POINTER :: mm_cell
323 : TYPE(mp_para_env_type), POINTER :: para_env
324 :
325 : CALL build_tb_qmmm_matrix_smeared(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
326 1116 : gaussian=.FALSE.)
327 :
328 1116 : END SUBROUTINE build_tb_qmmm_matrix_pc
329 :
330 : ! **************************************************************************************************
331 : !> \brief Constructs the tight-binding QM/MM Hamiltonian for Gaussian MM charges
332 : !> \param qs_env ...
333 : !> \param qmmm_env ...
334 : !> \param particles_mm ...
335 : !> \param mm_cell ...
336 : !> \param para_env ...
337 : ! **************************************************************************************************
338 220 : SUBROUTINE build_tb_qmmm_matrix_gauss(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
339 :
340 : TYPE(qs_environment_type), POINTER :: qs_env
341 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
342 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
343 : TYPE(cell_type), POINTER :: mm_cell
344 : TYPE(mp_para_env_type), POINTER :: para_env
345 :
346 220 : IF (qmmm_env%compatibility) THEN
347 0 : CPABORT("Gaussian QM/MM coupling for tight-binding methods requires NOCOMPATIBILITY.")
348 : END IF
349 : CALL build_tb_qmmm_matrix_smeared(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
350 220 : gaussian=.TRUE.)
351 :
352 220 : END SUBROUTINE build_tb_qmmm_matrix_gauss
353 :
354 : ! **************************************************************************************************
355 : !> \brief Constructs the tight-binding QM/MM Hamiltonian for smeared electrostatic coupling
356 : !> \param qs_env ...
357 : !> \param qmmm_env ...
358 : !> \param particles_mm ...
359 : !> \param mm_cell ...
360 : !> \param para_env ...
361 : !> \param gaussian Use Gaussian MM charges instead of the tight-binding point-charge regularization
362 : ! **************************************************************************************************
363 1336 : SUBROUTINE build_tb_qmmm_matrix_smeared(qs_env, qmmm_env, particles_mm, mm_cell, para_env, gaussian)
364 :
365 : TYPE(qs_environment_type), POINTER :: qs_env
366 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
367 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
368 : TYPE(cell_type), POINTER :: mm_cell
369 : TYPE(mp_para_env_type), POINTER :: para_env
370 : LOGICAL, INTENT(IN) :: gaussian
371 :
372 : CHARACTER(len=*), PARAMETER :: routineN = 'build_tb_qmmm_matrix_smeared'
373 :
374 : INTEGER :: do_ipol, ewald_type, handle, i, iatom, ikind, imm, imp, indmm, ipot, jatom, &
375 : natom, natorb, nkind, nmm, nonperiodic_pot_type, short_range_pot_type
376 1336 : INTEGER, DIMENSION(:), POINTER :: list
377 : LOGICAL :: defined, do_dftb, do_multipoles, do_xtb, &
378 : found
379 : REAL(KIND=dp) :: alpha, pc_ener, zeff
380 : REAL(KIND=dp), DIMENSION(0:3) :: eta_a
381 : REAL(KIND=dp), DIMENSION(2) :: rcutoff
382 1336 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges_mm, qpot
383 1336 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hblock, sblock
384 1336 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
385 : TYPE(dbcsr_iterator_type) :: iter
386 1336 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_s
387 : TYPE(dft_control_type), POINTER :: dft_control
388 : TYPE(dftb_control_type), POINTER :: dftb_control
389 : TYPE(ewald_environment_type), POINTER :: ewald_env
390 : TYPE(ewald_pw_type), POINTER :: ewald_pw
391 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
392 1336 : POINTER :: sab_nl
393 1336 : TYPE(particle_type), DIMENSION(:), POINTER :: atoms_mm, particles_qm
394 : TYPE(qmmm_pot_type), POINTER :: Pot
395 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
396 1336 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
397 : TYPE(qs_ks_env_type), POINTER :: ks_env
398 : TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
399 : TYPE(qs_rho_type), POINTER :: rho
400 : TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
401 : print_section
402 : TYPE(xtb_atom_type), POINTER :: xtb_kind
403 : TYPE(xtb_control_type), POINTER :: xtb_control
404 :
405 1336 : CALL timeset(routineN, handle)
406 :
407 1336 : IF (gaussian) THEN
408 220 : nonperiodic_pot_type = pot_gauss_nonperiodic
409 220 : short_range_pot_type = pot_gauss_short_range
410 : ELSE
411 1116 : nonperiodic_pot_type = pot_tb_nonperiodic
412 1116 : short_range_pot_type = pot_tb_short_range
413 : END IF
414 :
415 : CALL get_qs_env(qs_env=qs_env, &
416 : dft_control=dft_control, &
417 : atomic_kind_set=atomic_kind_set, &
418 : particle_set=particles_qm, &
419 : qs_kind_set=qs_kind_set, &
420 : rho=rho, &
421 1336 : natom=natom)
422 1336 : dftb_control => dft_control%qs_control%dftb_control
423 1336 : xtb_control => dft_control%qs_control%xtb_control
424 :
425 1336 : IF (dft_control%qs_control%dftb) THEN
426 : do_dftb = .TRUE.
427 : do_xtb = .FALSE.
428 668 : ELSE IF (dft_control%qs_control%xtb) THEN
429 : do_dftb = .FALSE.
430 : do_xtb = .TRUE.
431 : ELSE
432 0 : CPABORT("TB method unknown")
433 : END IF
434 :
435 1336 : CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
436 :
437 1336 : NULLIFY (matrix_s)
438 1336 : IF (do_dftb) THEN
439 668 : CALL build_dftb_overlap(qs_env, 0, matrix_s)
440 668 : ELSE IF (do_xtb) THEN
441 668 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
442 668 : CALL build_overlap_matrix(ks_env, matrix_s, basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
443 : END IF
444 :
445 4008 : ALLOCATE (qpot(natom))
446 5344 : qpot = 0.0_dp
447 1336 : pc_ener = 0.0_dp
448 :
449 : ! Create Ewald environments
450 1336 : poisson_section => section_vals_get_subs_vals(qs_env%input, "MM%POISSON")
451 24048 : ALLOCATE (ewald_env)
452 1336 : CALL ewald_env_create(ewald_env, para_env)
453 1336 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
454 1336 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
455 1336 : CALL read_ewald_section(ewald_env, ewald_section)
456 1336 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
457 1336 : ALLOCATE (ewald_pw)
458 1336 : CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=print_section)
459 :
460 1336 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, do_ipol=do_ipol)
461 1336 : IF (do_multipoles) CPABORT("No multipole force fields allowed in TB QM/MM")
462 1336 : IF (do_ipol /= do_fist_pol_none) CPABORT("No polarizable force fields allowed in TB QM/MM")
463 :
464 0 : SELECT CASE (ewald_type)
465 : CASE (do_ewald_pme)
466 0 : CPABORT("PME Ewald type not implemented for TB/QMMM")
467 : CASE (do_ewald_ewald, do_ewald_spme)
468 2334 : DO ipot = 1, SIZE(qmmm_env%Potentials)
469 1552 : Pot => qmmm_env%Potentials(ipot)%Pot
470 1552 : nmm = SIZE(Pot%mm_atom_index)
471 : ! get a 'clean' mm particle set
472 1552 : NULLIFY (atoms_mm)
473 1552 : CALL allocate_particle_set(atoms_mm, nmm)
474 4656 : ALLOCATE (charges_mm(nmm))
475 6208 : DO Imp = 1, nmm
476 4656 : Imm = Pot%mm_atom_index(Imp)
477 4656 : IndMM = qmmm_env%mm_atom_index(Imm)
478 37248 : atoms_mm(Imp)%r = particles_mm(IndMM)%r
479 4656 : atoms_mm(Imp)%atomic_kind => particles_mm(IndMM)%atomic_kind
480 6208 : charges_mm(Imp) = qmmm_env%mm_atom_chrg(Imm)
481 : END DO
482 1552 : IF (ewald_type == do_ewald_ewald) THEN
483 0 : CPABORT("Ewald not implemented for TB/QMMM")
484 1552 : ELSE IF (ewald_type == do_ewald_spme) THEN
485 : ! spme electrostatic potential
486 1552 : CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, particles_qm, qpot)
487 : END IF
488 1552 : CALL deallocate_particle_set(atoms_mm)
489 2334 : DEALLOCATE (charges_mm)
490 : END DO
491 782 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
492 0 : DO ipot = 1, SIZE(qmmm_env%added_charges%Potentials)
493 0 : Pot => qmmm_env%added_charges%Potentials(ipot)%Pot
494 0 : nmm = SIZE(Pot%mm_atom_index)
495 : ! get a 'clean' mm particle set
496 0 : NULLIFY (atoms_mm)
497 0 : CALL allocate_particle_set(atoms_mm, nmm)
498 0 : ALLOCATE (charges_mm(nmm))
499 0 : DO Imp = 1, nmm
500 0 : Imm = Pot%mm_atom_index(Imp)
501 0 : IndMM = qmmm_env%added_charges%mm_atom_index(Imm)
502 0 : atoms_mm(Imp)%r = qmmm_env%added_charges%added_particles(IndMM)%r
503 0 : atoms_mm(Imp)%atomic_kind => qmmm_env%added_charges%added_particles(IndMM)%atomic_kind
504 0 : charges_mm(Imp) = qmmm_env%added_charges%mm_atom_chrg(Imm)
505 : END DO
506 0 : IF (ewald_type == do_ewald_ewald) THEN
507 0 : CPABORT("Ewald not implemented for TB/QMMM")
508 0 : ELSE IF (ewald_type == do_ewald_spme) THEN
509 : ! spme electrostatic potential
510 0 : CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, particles_qm, qpot)
511 : END IF
512 0 : CALL deallocate_particle_set(atoms_mm)
513 782 : DEALLOCATE (charges_mm)
514 : END DO
515 : END IF
516 5474 : CALL para_env%sum(qpot)
517 : ! Add the Ewald real-space term and the method-specific short-range correction
518 : ! This is effectively using a minimum image convention!
519 : ! Set rcutoff to values compatible with alpha Ewald
520 782 : CALL ewald_env_get(ewald_env, rcut=rcutoff(1), alpha=alpha)
521 782 : rcutoff(2) = 0.025_dp*rcutoff(1)
522 782 : rcutoff(1) = 2.0_dp*rcutoff(1)
523 782 : nkind = SIZE(atomic_kind_set)
524 2346 : DO ikind = 1, nkind
525 1564 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
526 1564 : IF (do_dftb) THEN
527 672 : NULLIFY (dftb_kind)
528 672 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
529 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff, &
530 672 : defined=defined, eta=eta_a, natorb=natorb)
531 : ! use mm charge smearing for non-scc cases
532 672 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
533 672 : IF (.NOT. defined .OR. natorb < 1) CYCLE
534 892 : ELSE IF (do_xtb) THEN
535 892 : NULLIFY (xtb_kind)
536 892 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
537 892 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
538 892 : eta_a(0) = eta_mm
539 : END IF
540 4692 : DO i = 1, SIZE(list)
541 2346 : iatom = list(i)
542 : CALL build_mm_pot(qpot(iatom), short_range_pot_type, eta_a(0), &
543 : qmmm_env%Potentials, particles_mm, &
544 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, rcutoff, &
545 2346 : particles_qm)
546 : CALL build_mm_pot(qpot(iatom), pot_ewald_short_range, alpha, &
547 : qmmm_env%Potentials, particles_mm, &
548 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, rcutoff, &
549 2346 : particles_qm)
550 : ! Possibly added charges
551 2346 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
552 : CALL build_mm_pot(qpot(iatom), short_range_pot_type, eta_a(0), &
553 : qmmm_env%added_charges%potentials, &
554 : qmmm_env%added_charges%added_particles, &
555 : qmmm_env%added_charges%mm_atom_chrg, &
556 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, rcutoff, &
557 0 : particles_qm)
558 : CALL build_mm_pot(qpot(iatom), pot_ewald_short_range, alpha, &
559 : qmmm_env%added_charges%potentials, &
560 : qmmm_env%added_charges%added_particles, &
561 : qmmm_env%added_charges%mm_atom_chrg, &
562 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, rcutoff, &
563 0 : particles_qm)
564 : END IF
565 3910 : pc_ener = pc_ener + qpot(iatom)*zeff
566 : END DO
567 : END DO
568 : CASE (do_ewald_none)
569 : ! Directly sum the nonperiodic regularized electrostatic potential
570 554 : nkind = SIZE(atomic_kind_set)
571 1662 : DO ikind = 1, nkind
572 1108 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
573 1108 : IF (do_dftb) THEN
574 664 : NULLIFY (dftb_kind)
575 664 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
576 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff, &
577 664 : defined=defined, eta=eta_a, natorb=natorb)
578 : ! use mm charge smearing for non-scc cases
579 664 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
580 664 : IF (.NOT. defined .OR. natorb < 1) CYCLE
581 444 : ELSE IF (do_xtb) THEN
582 444 : NULLIFY (xtb_kind)
583 444 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
584 444 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
585 444 : eta_a(0) = eta_mm
586 : END IF
587 3324 : DO i = 1, SIZE(list)
588 1662 : iatom = list(i)
589 : CALL build_mm_pot(qpot(iatom), nonperiodic_pot_type, eta_a(0), &
590 : qmmm_env%Potentials, particles_mm, &
591 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, &
592 1662 : qmmm_env%spherical_cutoff, particles_qm)
593 : ! Possibly added charges
594 1662 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
595 : CALL build_mm_pot(qpot(iatom), nonperiodic_pot_type, eta_a(0), &
596 : qmmm_env%added_charges%potentials, &
597 : qmmm_env%added_charges%added_particles, &
598 : qmmm_env%added_charges%mm_atom_chrg, &
599 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
600 : qmmm_env%spherical_cutoff, &
601 0 : particles_qm)
602 : END IF
603 2770 : pc_ener = pc_ener + qpot(iatom)*zeff
604 : END DO
605 : END DO
606 : CASE DEFAULT
607 1336 : CPABORT("Unknown Ewald type!")
608 : END SELECT
609 :
610 : ! Allocate the core Hamiltonian matrix
611 1336 : CALL get_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env_loc)
612 1336 : matrix_h => ks_qmmm_env_loc%matrix_h
613 1336 : CALL dbcsr_allocate_matrix_set(matrix_h, 1)
614 1336 : ALLOCATE (matrix_h(1)%matrix)
615 : CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, &
616 1336 : name="QMMM HAMILTONIAN MATRIX")
617 1336 : CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
618 :
619 1336 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
620 5344 : DO WHILE (dbcsr_iterator_blocks_left(iter))
621 4008 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock)
622 4008 : NULLIFY (hblock)
623 : CALL dbcsr_get_block_p(matrix=matrix_h(1)%matrix, &
624 4008 : row=iatom, col=jatom, block=hblock, found=found)
625 4008 : CPASSERT(found)
626 76624 : hblock = hblock - 0.5_dp*sblock*(qpot(iatom) + qpot(jatom))
627 : END DO
628 1336 : CALL dbcsr_iterator_stop(iter)
629 :
630 1336 : ks_qmmm_env_loc%matrix_h => matrix_h
631 1336 : ks_qmmm_env_loc%pc_ener = pc_ener
632 :
633 1336 : DEALLOCATE (qpot)
634 :
635 : ! Release Ewald environment
636 1336 : CALL ewald_env_release(ewald_env)
637 1336 : DEALLOCATE (ewald_env)
638 1336 : CALL ewald_pw_release(ewald_pw)
639 1336 : DEALLOCATE (ewald_pw)
640 :
641 1336 : CALL dbcsr_deallocate_matrix_set(matrix_s)
642 :
643 1336 : CALL timestop(handle)
644 :
645 5344 : END SUBROUTINE build_tb_qmmm_matrix_smeared
646 :
647 : ! **************************************************************************************************
648 : !> \brief Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms
649 : !> \param qs_env ...
650 : !> \param qmmm_env ...
651 : !> \param particles_mm ...
652 : !> \param mm_cell ...
653 : !> \param para_env ...
654 : !> \param calc_force ...
655 : !> \param Forces ...
656 : !> \param Forces_added_charges ...
657 : !> \author JGH 10.2014 [created]
658 : ! **************************************************************************************************
659 448 : SUBROUTINE deriv_tb_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
660 : calc_force, Forces, Forces_added_charges)
661 :
662 : TYPE(qs_environment_type), POINTER :: qs_env
663 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
664 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
665 : TYPE(cell_type), POINTER :: mm_cell
666 : TYPE(mp_para_env_type), POINTER :: para_env
667 : LOGICAL, INTENT(in), OPTIONAL :: calc_force
668 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces, Forces_added_charges
669 :
670 : CHARACTER(len=*), PARAMETER :: routineN = 'deriv_tb_qmmm_matrix'
671 :
672 : INTEGER :: atom_a, handle, i, iatom, ikind, iqm, &
673 : jatom, natom, natorb, nkind, nspins, &
674 : number_qm_atoms
675 448 : INTEGER, DIMENSION(:), POINTER :: list
676 : LOGICAL :: defined, do_dftb, do_xtb, found
677 : REAL(KIND=dp) :: fi, gmij, zeff
678 : REAL(KIND=dp), DIMENSION(0:3) :: eta_a
679 448 : REAL(KIND=dp), DIMENSION(:), POINTER :: mcharge, qpot
680 448 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges, dsblock, Forces_QM, pblock, &
681 448 : sblock
682 448 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
683 : TYPE(dbcsr_iterator_type) :: iter
684 448 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s
685 : TYPE(dft_control_type), POINTER :: dft_control
686 : TYPE(dftb_control_type), POINTER :: dftb_control
687 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
688 448 : POINTER :: sab_nl
689 448 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
690 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
691 448 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
692 : TYPE(qs_ks_env_type), POINTER :: ks_env
693 : TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
694 : TYPE(qs_rho_type), POINTER :: rho
695 : TYPE(xtb_atom_type), POINTER :: xtb_kind
696 : TYPE(xtb_control_type), POINTER :: xtb_control
697 :
698 448 : CALL timeset(routineN, handle)
699 448 : IF (calc_force) THEN
700 16 : NULLIFY (rho, atomic_kind_set, qs_kind_set, particles_qm)
701 : CALL get_qs_env(qs_env=qs_env, &
702 : rho=rho, &
703 : atomic_kind_set=atomic_kind_set, &
704 : qs_kind_set=qs_kind_set, &
705 : ks_qmmm_env=ks_qmmm_env_loc, &
706 : dft_control=dft_control, &
707 : particle_set=particles_qm, &
708 16 : natom=number_qm_atoms)
709 16 : dftb_control => dft_control%qs_control%dftb_control
710 16 : xtb_control => dft_control%qs_control%xtb_control
711 :
712 16 : IF (dft_control%qs_control%dftb) THEN
713 : do_dftb = .TRUE.
714 : do_xtb = .FALSE.
715 8 : ELSE IF (dft_control%qs_control%xtb) THEN
716 : do_dftb = .FALSE.
717 : do_xtb = .TRUE.
718 : ELSE
719 0 : CPABORT("TB method unknown")
720 : END IF
721 :
722 16 : NULLIFY (matrix_s)
723 0 : IF (do_dftb) THEN
724 8 : CALL build_dftb_overlap(qs_env, 1, matrix_s)
725 : ELSE IF (do_xtb) THEN
726 8 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
727 : CALL build_overlap_matrix(ks_env, matrix_s, nderivative=1, &
728 8 : basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
729 : END IF
730 :
731 16 : CALL qs_rho_get(rho, rho_ao=matrix_p)
732 :
733 16 : nspins = dft_control%nspins
734 16 : nkind = SIZE(atomic_kind_set)
735 : ! Mulliken charges
736 64 : ALLOCATE (charges(number_qm_atoms, nspins))
737 : !
738 16 : CALL mulliken_charges(matrix_p, matrix_s(1)%matrix, para_env, charges)
739 : !
740 48 : ALLOCATE (mcharge(number_qm_atoms))
741 48 : DO ikind = 1, nkind
742 32 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
743 32 : IF (do_dftb) THEN
744 16 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
745 16 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
746 16 : ELSE IF (do_xtb) THEN
747 16 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
748 16 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
749 : END IF
750 128 : DO iatom = 1, natom
751 48 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
752 128 : mcharge(atom_a) = zeff - SUM(charges(atom_a, 1:nspins))
753 : END DO
754 : END DO
755 16 : DEALLOCATE (charges)
756 :
757 48 : ALLOCATE (qpot(number_qm_atoms))
758 64 : qpot = 0.0_dp
759 48 : ALLOCATE (Forces_QM(3, number_qm_atoms))
760 208 : Forces_QM = 0.0_dp
761 :
762 : ! calculate potential and forces from classical charges
763 : iqm = 0
764 48 : DO ikind = 1, nkind
765 32 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
766 32 : IF (do_dftb) THEN
767 16 : NULLIFY (dftb_kind)
768 16 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
769 : CALL get_dftb_atom_param(dftb_kind, &
770 16 : defined=defined, eta=eta_a, natorb=natorb)
771 : ! use mm charge smearing for non-scc cases
772 16 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
773 16 : IF (.NOT. defined .OR. natorb < 1) CYCLE
774 16 : ELSE IF (do_xtb) THEN
775 16 : eta_a(0) = eta_mm
776 : END IF
777 96 : DO i = 1, SIZE(list)
778 48 : iatom = list(i)
779 48 : iqm = iqm + 1
780 : CALL build_mm_pot(qpot(iatom), pot_tb_nonperiodic, eta_a(0), &
781 : qmmm_env%Potentials, particles_mm, &
782 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, &
783 48 : qmmm_env%spherical_cutoff, particles_qm)
784 : CALL build_mm_dpot(mcharge(iatom), pot_tb_nonperiodic, eta_a(0), &
785 : qmmm_env%Potentials, particles_mm, &
786 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
787 : mm_cell, iatom, Forces, Forces_QM(:, iqm), &
788 48 : qmmm_env%spherical_cutoff, particles_qm)
789 : ! Possibly added charges
790 80 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
791 : CALL build_mm_pot(qpot(iatom), pot_tb_nonperiodic, eta_a(0), &
792 : qmmm_env%added_charges%potentials, &
793 : qmmm_env%added_charges%added_particles, &
794 : qmmm_env%added_charges%mm_atom_chrg, &
795 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
796 : qmmm_env%spherical_cutoff, &
797 0 : particles_qm)
798 : CALL build_mm_dpot(mcharge(iatom), pot_tb_nonperiodic, eta_a(0), &
799 : qmmm_env%added_charges%potentials, &
800 : qmmm_env%added_charges%added_particles, &
801 : qmmm_env%added_charges%mm_atom_chrg, &
802 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
803 : Forces_added_charges, &
804 0 : Forces_QM(:, iqm), qmmm_env%spherical_cutoff, particles_qm)
805 : END IF
806 : END DO
807 : END DO
808 :
809 : ! Transfer QM gradients to the QM particles..
810 : iqm = 0
811 48 : DO ikind = 1, nkind
812 32 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
813 32 : IF (do_dftb) THEN
814 16 : NULLIFY (dftb_kind)
815 16 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
816 16 : CALL get_dftb_atom_param(dftb_kind, defined=defined, natorb=natorb)
817 16 : IF (.NOT. defined .OR. natorb < 1) CYCLE
818 : ELSE IF (do_xtb) THEN
819 : ! use all kinds
820 : END IF
821 96 : DO i = 1, SIZE(list)
822 48 : iqm = iqm + 1
823 48 : iatom = qmmm_env%qm_atom_index(list(i))
824 416 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
825 : END DO
826 : END DO
827 :
828 : ! derivatives from qm charges
829 208 : Forces_QM = 0.0_dp
830 16 : IF (SIZE(matrix_p) == 2) THEN
831 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
832 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
833 : END IF
834 : !
835 16 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
836 64 : DO WHILE (dbcsr_iterator_blocks_left(iter))
837 48 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock)
838 : !
839 48 : IF (iatom == jatom) CYCLE
840 : !
841 24 : gmij = -0.5_dp*(qpot(iatom) + qpot(jatom))
842 24 : NULLIFY (pblock)
843 : CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
844 24 : row=iatom, col=jatom, block=pblock, found=found)
845 24 : CPASSERT(found)
846 112 : DO i = 1, 3
847 72 : NULLIFY (dsblock)
848 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
849 72 : row=iatom, col=jatom, block=dsblock, found=found)
850 72 : CPASSERT(found)
851 648 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
852 72 : Forces_QM(i, iatom) = Forces_QM(i, iatom) + fi
853 192 : Forces_QM(i, jatom) = Forces_QM(i, jatom) - fi
854 : END DO
855 : END DO
856 16 : CALL dbcsr_iterator_stop(iter)
857 : !
858 16 : IF (SIZE(matrix_p) == 2) THEN
859 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
860 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
861 : END IF
862 : !
863 : ! Transfer QM gradients to the QM particles..
864 400 : CALL para_env%sum(Forces_QM)
865 48 : DO ikind = 1, nkind
866 32 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
867 96 : DO i = 1, SIZE(list)
868 48 : iqm = list(i)
869 48 : iatom = qmmm_env%qm_atom_index(iqm)
870 416 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
871 : END DO
872 : END DO
873 : !
874 16 : DEALLOCATE (mcharge)
875 : !
876 : ! MM forces will be handled directly from the QMMM module in the same way
877 : ! as for GPW/GAPW methods
878 16 : DEALLOCATE (Forces_QM)
879 16 : DEALLOCATE (qpot)
880 :
881 32 : CALL dbcsr_deallocate_matrix_set(matrix_s)
882 :
883 : END IF
884 448 : CALL timestop(handle)
885 :
886 448 : END SUBROUTINE deriv_tb_qmmm_matrix
887 :
888 : ! **************************************************************************************************
889 : !> \brief Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms
890 : !> \param qs_env ...
891 : !> \param qmmm_env ...
892 : !> \param particles_mm ...
893 : !> \param mm_cell ...
894 : !> \param para_env ...
895 : !> \param calc_force ...
896 : !> \param Forces ...
897 : !> \param Forces_added_charges ...
898 : !> \author JGH 10.2014 [created]
899 : ! **************************************************************************************************
900 1116 : SUBROUTINE deriv_tb_qmmm_matrix_pc(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
901 : calc_force, Forces, Forces_added_charges)
902 :
903 : TYPE(qs_environment_type), POINTER :: qs_env
904 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
905 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
906 : TYPE(cell_type), POINTER :: mm_cell
907 : TYPE(mp_para_env_type), POINTER :: para_env
908 : LOGICAL, INTENT(IN), OPTIONAL :: calc_force
909 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces, Forces_added_charges
910 :
911 : CALL deriv_tb_qmmm_matrix_smeared(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
912 1116 : calc_force, Forces, Forces_added_charges, gaussian=.FALSE.)
913 :
914 1116 : END SUBROUTINE deriv_tb_qmmm_matrix_pc
915 :
916 : ! **************************************************************************************************
917 : !> \brief Constructs tight-binding QM/MM derivatives for Gaussian MM charges
918 : !> \param qs_env ...
919 : !> \param qmmm_env ...
920 : !> \param particles_mm ...
921 : !> \param mm_cell ...
922 : !> \param para_env ...
923 : !> \param calc_force ...
924 : !> \param Forces ...
925 : !> \param Forces_added_charges ...
926 : ! **************************************************************************************************
927 220 : SUBROUTINE deriv_tb_qmmm_matrix_gauss(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
928 : calc_force, Forces, Forces_added_charges)
929 :
930 : TYPE(qs_environment_type), POINTER :: qs_env
931 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
932 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
933 : TYPE(cell_type), POINTER :: mm_cell
934 : TYPE(mp_para_env_type), POINTER :: para_env
935 : LOGICAL, INTENT(IN), OPTIONAL :: calc_force
936 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces, Forces_added_charges
937 :
938 : CALL deriv_tb_qmmm_matrix_smeared(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
939 220 : calc_force, Forces, Forces_added_charges, gaussian=.TRUE.)
940 :
941 220 : END SUBROUTINE deriv_tb_qmmm_matrix_gauss
942 :
943 : ! **************************************************************************************************
944 : !> \brief Constructs tight-binding QM/MM derivatives for smeared electrostatic coupling
945 : !> \param qs_env ...
946 : !> \param qmmm_env ...
947 : !> \param particles_mm ...
948 : !> \param mm_cell ...
949 : !> \param para_env ...
950 : !> \param calc_force ...
951 : !> \param Forces ...
952 : !> \param Forces_added_charges ...
953 : !> \param gaussian Use Gaussian MM charges instead of the tight-binding point-charge regularization
954 : ! **************************************************************************************************
955 1336 : SUBROUTINE deriv_tb_qmmm_matrix_smeared(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
956 : calc_force, Forces, Forces_added_charges, gaussian)
957 :
958 : TYPE(qs_environment_type), POINTER :: qs_env
959 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
960 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
961 : TYPE(cell_type), POINTER :: mm_cell
962 : TYPE(mp_para_env_type), POINTER :: para_env
963 : LOGICAL, INTENT(in), OPTIONAL :: calc_force
964 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces, Forces_added_charges
965 : LOGICAL, INTENT(IN) :: gaussian
966 :
967 : CHARACTER(len=*), PARAMETER :: routineN = 'deriv_tb_qmmm_matrix_smeared'
968 :
969 : INTEGER :: atom_a, do_ipol, ewald_type, handle, i, iatom, ikind, imm, imp, indmm, ipot, iqm, &
970 : jatom, natom, natorb, nkind, nmm, nonperiodic_pot_type, nspins, number_qm_atoms, &
971 : short_range_pot_type
972 1336 : INTEGER, DIMENSION(:), POINTER :: list
973 : LOGICAL :: defined, do_dftb, do_multipoles, do_xtb, &
974 : found
975 : REAL(KIND=dp) :: alpha, fi, gmij, zeff
976 : REAL(KIND=dp), DIMENSION(0:3) :: eta_a
977 : REAL(KIND=dp), DIMENSION(2) :: rcutoff
978 1336 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges_mm, mcharge, qpot
979 1336 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges, dsblock, Forces_MM, Forces_QM, &
980 1336 : pblock, sblock
981 1336 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
982 : TYPE(dbcsr_iterator_type) :: iter
983 1336 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s
984 : TYPE(dft_control_type), POINTER :: dft_control
985 : TYPE(dftb_control_type), POINTER :: dftb_control
986 : TYPE(ewald_environment_type), POINTER :: ewald_env
987 : TYPE(ewald_pw_type), POINTER :: ewald_pw
988 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
989 1336 : POINTER :: sab_nl
990 1336 : TYPE(particle_type), DIMENSION(:), POINTER :: atoms_mm, particles_qm
991 : TYPE(qmmm_pot_type), POINTER :: Pot
992 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
993 1336 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
994 : TYPE(qs_ks_env_type), POINTER :: ks_env
995 : TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
996 : TYPE(qs_rho_type), POINTER :: rho
997 : TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
998 : print_section
999 : TYPE(xtb_atom_type), POINTER :: xtb_kind
1000 : TYPE(xtb_control_type), POINTER :: xtb_control
1001 :
1002 1336 : CALL timeset(routineN, handle)
1003 1336 : IF (gaussian) THEN
1004 220 : nonperiodic_pot_type = pot_gauss_nonperiodic
1005 220 : short_range_pot_type = pot_gauss_short_range
1006 : ELSE
1007 1116 : nonperiodic_pot_type = pot_tb_nonperiodic
1008 1116 : short_range_pot_type = pot_tb_short_range
1009 : END IF
1010 1336 : IF (calc_force) THEN
1011 40 : NULLIFY (rho, atomic_kind_set, qs_kind_set, particles_qm)
1012 : CALL get_qs_env(qs_env=qs_env, &
1013 : rho=rho, &
1014 : atomic_kind_set=atomic_kind_set, &
1015 : qs_kind_set=qs_kind_set, &
1016 : ks_qmmm_env=ks_qmmm_env_loc, &
1017 : dft_control=dft_control, &
1018 : particle_set=particles_qm, &
1019 40 : natom=number_qm_atoms)
1020 40 : dftb_control => dft_control%qs_control%dftb_control
1021 40 : xtb_control => dft_control%qs_control%xtb_control
1022 :
1023 40 : IF (dft_control%qs_control%dftb) THEN
1024 : do_dftb = .TRUE.
1025 : do_xtb = .FALSE.
1026 20 : ELSE IF (dft_control%qs_control%xtb) THEN
1027 : do_dftb = .FALSE.
1028 : do_xtb = .TRUE.
1029 : ELSE
1030 0 : CPABORT("TB method unknown")
1031 : END IF
1032 :
1033 40 : NULLIFY (matrix_s)
1034 0 : IF (do_dftb) THEN
1035 20 : CALL build_dftb_overlap(qs_env, 1, matrix_s)
1036 : ELSE IF (do_xtb) THEN
1037 20 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
1038 : CALL build_overlap_matrix(ks_env, matrix_s, nderivative=1, &
1039 20 : basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
1040 : END IF
1041 40 : CALL qs_rho_get(rho, rho_ao=matrix_p)
1042 :
1043 40 : nspins = dft_control%nspins
1044 40 : nkind = SIZE(atomic_kind_set)
1045 : ! Mulliken charges
1046 160 : ALLOCATE (charges(number_qm_atoms, nspins))
1047 : !
1048 40 : CALL mulliken_charges(matrix_p, matrix_s(1)%matrix, para_env, charges)
1049 : !
1050 120 : ALLOCATE (mcharge(number_qm_atoms))
1051 120 : DO ikind = 1, nkind
1052 80 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
1053 80 : IF (do_dftb) THEN
1054 40 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1055 40 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
1056 40 : ELSE IF (do_xtb) THEN
1057 40 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1058 40 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
1059 : END IF
1060 320 : DO iatom = 1, natom
1061 120 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1062 320 : mcharge(atom_a) = zeff - SUM(charges(atom_a, 1:nspins))
1063 : END DO
1064 : END DO
1065 40 : DEALLOCATE (charges)
1066 :
1067 120 : ALLOCATE (qpot(number_qm_atoms))
1068 160 : qpot = 0.0_dp
1069 120 : ALLOCATE (Forces_QM(3, number_qm_atoms))
1070 520 : Forces_QM = 0.0_dp
1071 :
1072 : ! Create Ewald environments
1073 40 : poisson_section => section_vals_get_subs_vals(qs_env%input, "MM%POISSON")
1074 720 : ALLOCATE (ewald_env)
1075 40 : CALL ewald_env_create(ewald_env, para_env)
1076 40 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
1077 40 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
1078 40 : CALL read_ewald_section(ewald_env, ewald_section)
1079 40 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
1080 40 : ALLOCATE (ewald_pw)
1081 40 : CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=print_section)
1082 :
1083 40 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, do_ipol=do_ipol)
1084 40 : IF (do_multipoles) CPABORT("No multipole force fields allowed in DFTB QM/MM")
1085 40 : IF (do_ipol /= do_fist_pol_none) CPABORT("No polarizable force fields allowed in DFTB QM/MM")
1086 :
1087 0 : SELECT CASE (ewald_type)
1088 : CASE (do_ewald_pme)
1089 0 : CPABORT("PME Ewald type not implemented for DFTB/QMMM")
1090 : CASE (do_ewald_ewald, do_ewald_spme)
1091 66 : DO ipot = 1, SIZE(qmmm_env%Potentials)
1092 40 : Pot => qmmm_env%Potentials(ipot)%Pot
1093 40 : nmm = SIZE(Pot%mm_atom_index)
1094 : ! get a 'clean' mm particle set
1095 40 : NULLIFY (atoms_mm)
1096 40 : CALL allocate_particle_set(atoms_mm, nmm)
1097 120 : ALLOCATE (charges_mm(nmm))
1098 160 : DO Imp = 1, nmm
1099 120 : Imm = Pot%mm_atom_index(Imp)
1100 120 : IndMM = qmmm_env%mm_atom_index(Imm)
1101 960 : atoms_mm(Imp)%r = particles_mm(IndMM)%r
1102 120 : atoms_mm(Imp)%atomic_kind => particles_mm(IndMM)%atomic_kind
1103 160 : charges_mm(Imp) = qmmm_env%mm_atom_chrg(Imm)
1104 : END DO
1105 : ! force array for mm atoms
1106 120 : ALLOCATE (Forces_MM(3, nmm))
1107 520 : Forces_MM = 0.0_dp
1108 40 : IF (ewald_type == do_ewald_ewald) THEN
1109 0 : CPABORT("Ewald not implemented for DFTB/QMMM")
1110 40 : ELSE IF (ewald_type == do_ewald_spme) THEN
1111 : ! spme electrostatic potential
1112 : CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
1113 40 : particles_qm, qpot)
1114 : ! forces QM
1115 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
1116 40 : particles_qm, mcharge, Forces_QM)
1117 : ! forces MM
1118 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, particles_qm, mcharge, &
1119 40 : atoms_mm, charges_mm, Forces_MM)
1120 : END IF
1121 40 : CALL deallocate_particle_set(atoms_mm)
1122 40 : DEALLOCATE (charges_mm)
1123 : ! transfer MM forces
1124 1000 : CALL para_env%sum(Forces_MM)
1125 160 : DO Imp = 1, nmm
1126 120 : Imm = Pot%mm_atom_index(Imp)
1127 1000 : Forces(:, Imm) = Forces(:, Imm) - Forces_MM(:, Imp)
1128 : END DO
1129 66 : DEALLOCATE (Forces_MM)
1130 : END DO
1131 :
1132 26 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
1133 0 : DO ipot = 1, SIZE(qmmm_env%added_charges%Potentials)
1134 0 : Pot => qmmm_env%added_charges%Potentials(ipot)%Pot
1135 0 : nmm = SIZE(Pot%mm_atom_index)
1136 : ! get a 'clean' mm particle set
1137 0 : NULLIFY (atoms_mm)
1138 0 : CALL allocate_particle_set(atoms_mm, nmm)
1139 0 : ALLOCATE (charges_mm(nmm))
1140 0 : DO Imp = 1, nmm
1141 0 : Imm = Pot%mm_atom_index(Imp)
1142 0 : IndMM = qmmm_env%added_charges%mm_atom_index(Imm)
1143 0 : atoms_mm(Imp)%r = qmmm_env%added_charges%added_particles(IndMM)%r
1144 0 : atoms_mm(Imp)%atomic_kind => qmmm_env%added_charges%added_particles(IndMM)%atomic_kind
1145 0 : charges_mm(Imp) = qmmm_env%added_charges%mm_atom_chrg(Imm)
1146 : END DO
1147 : ! force array for mm atoms
1148 0 : ALLOCATE (Forces_MM(3, nmm))
1149 0 : Forces_MM = 0.0_dp
1150 0 : IF (ewald_type == do_ewald_ewald) THEN
1151 0 : CPABORT("Ewald not implemented for DFTB/QMMM")
1152 0 : ELSE IF (ewald_type == do_ewald_spme) THEN
1153 : ! spme electrostatic potential
1154 : CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, &
1155 0 : charges_mm, particles_qm, qpot)
1156 : ! forces QM
1157 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
1158 0 : particles_qm, mcharge, Forces_QM)
1159 : ! forces MM
1160 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, particles_qm, mcharge, &
1161 0 : atoms_mm, charges_mm, Forces_MM)
1162 : END IF
1163 0 : CALL deallocate_particle_set(atoms_mm)
1164 : ! transfer MM forces
1165 0 : CALL para_env%sum(Forces_MM)
1166 0 : DO Imp = 1, nmm
1167 0 : Imm = Pot%mm_atom_index(Imp)
1168 0 : Forces_added_charges(:, Imm) = Forces_added_charges(:, Imm) - Forces_MM(:, Imp)
1169 : END DO
1170 26 : DEALLOCATE (Forces_MM)
1171 : END DO
1172 : END IF
1173 182 : CALL para_env%sum(qpot)
1174 650 : CALL para_env%sum(Forces_QM)
1175 : ! Add the Ewald real-space term and the method-specific short-range correction
1176 : ! This is effectively using a minimum image convention!
1177 : ! Set rcutoff to values compatible with alpha Ewald
1178 26 : CALL ewald_env_get(ewald_env, rcut=rcutoff(1), alpha=alpha)
1179 26 : rcutoff(2) = 0.025_dp*rcutoff(1)
1180 26 : rcutoff(1) = 2.0_dp*rcutoff(1)
1181 26 : nkind = SIZE(atomic_kind_set)
1182 26 : iqm = 0
1183 78 : DO ikind = 1, nkind
1184 52 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1185 52 : IF (do_dftb) THEN
1186 24 : NULLIFY (dftb_kind)
1187 24 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1188 : CALL get_dftb_atom_param(dftb_kind, &
1189 24 : defined=defined, eta=eta_a, natorb=natorb)
1190 : ! use mm charge smearing for non-scc cases
1191 24 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
1192 24 : IF (.NOT. defined .OR. natorb < 1) CYCLE
1193 28 : ELSE IF (do_xtb) THEN
1194 28 : eta_a(0) = eta_mm
1195 : END IF
1196 156 : DO i = 1, SIZE(list)
1197 78 : iatom = list(i)
1198 78 : iqm = iqm + 1
1199 : CALL build_mm_pot(qpot(iatom), short_range_pot_type, eta_a(0), &
1200 : qmmm_env%Potentials, particles_mm, &
1201 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1202 78 : mm_cell, iatom, rcutoff, particles_qm)
1203 : CALL build_mm_dpot(mcharge(iatom), short_range_pot_type, eta_a(0), &
1204 : qmmm_env%Potentials, particles_mm, &
1205 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1206 : mm_cell, iatom, Forces, Forces_QM(:, iqm), &
1207 78 : rcutoff, particles_qm)
1208 : CALL build_mm_pot(qpot(iatom), pot_ewald_short_range, alpha, &
1209 : qmmm_env%Potentials, particles_mm, &
1210 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1211 78 : mm_cell, iatom, rcutoff, particles_qm)
1212 : CALL build_mm_dpot(mcharge(iatom), pot_ewald_short_range, alpha, &
1213 : qmmm_env%Potentials, particles_mm, &
1214 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1215 : mm_cell, iatom, Forces, Forces_QM(:, iqm), &
1216 78 : rcutoff, particles_qm)
1217 : ! Possibly added charges
1218 130 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
1219 : CALL build_mm_pot(qpot(iatom), short_range_pot_type, eta_a(0), &
1220 : qmmm_env%added_charges%potentials, &
1221 : qmmm_env%added_charges%added_particles, &
1222 : qmmm_env%added_charges%mm_atom_chrg, &
1223 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, rcutoff, &
1224 0 : particles_qm)
1225 : CALL build_mm_dpot( &
1226 : mcharge(iatom), short_range_pot_type, eta_a(0), &
1227 : qmmm_env%added_charges%potentials, &
1228 : qmmm_env%added_charges%added_particles, qmmm_env%added_charges%mm_atom_chrg, &
1229 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
1230 : Forces_added_charges, Forces_QM(:, iqm), &
1231 0 : rcutoff, particles_qm)
1232 : CALL build_mm_pot(qpot(iatom), pot_ewald_short_range, alpha, &
1233 : qmmm_env%added_charges%potentials, &
1234 : qmmm_env%added_charges%added_particles, &
1235 : qmmm_env%added_charges%mm_atom_chrg, &
1236 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, rcutoff, &
1237 0 : particles_qm)
1238 : CALL build_mm_dpot( &
1239 : mcharge(iatom), pot_ewald_short_range, alpha, &
1240 : qmmm_env%added_charges%potentials, &
1241 : qmmm_env%added_charges%added_particles, qmmm_env%added_charges%mm_atom_chrg, &
1242 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
1243 : Forces_added_charges, Forces_QM(:, iqm), &
1244 0 : rcutoff, particles_qm)
1245 : END IF
1246 : END DO
1247 : END DO
1248 :
1249 : CASE (do_ewald_none)
1250 : ! Directly sum the nonperiodic regularized electrostatic potential
1251 : ! calculate potential and forces from classical charges
1252 : iqm = 0
1253 42 : DO ikind = 1, nkind
1254 28 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1255 28 : IF (do_dftb) THEN
1256 16 : NULLIFY (dftb_kind)
1257 16 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1258 : CALL get_dftb_atom_param(dftb_kind, &
1259 16 : defined=defined, eta=eta_a, natorb=natorb)
1260 : ! use mm charge smearing for non-scc cases
1261 16 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
1262 16 : IF (.NOT. defined .OR. natorb < 1) CYCLE
1263 12 : ELSE IF (do_xtb) THEN
1264 12 : eta_a(0) = eta_mm
1265 : END IF
1266 84 : DO i = 1, SIZE(list)
1267 42 : iatom = list(i)
1268 42 : iqm = iqm + 1
1269 : CALL build_mm_pot(qpot(iatom), nonperiodic_pot_type, eta_a(0), &
1270 : qmmm_env%Potentials, particles_mm, &
1271 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, &
1272 42 : qmmm_env%spherical_cutoff, particles_qm)
1273 : CALL build_mm_dpot(mcharge(iatom), nonperiodic_pot_type, eta_a(0), &
1274 : qmmm_env%Potentials, particles_mm, &
1275 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1276 : mm_cell, iatom, Forces, Forces_QM(:, iqm), &
1277 42 : qmmm_env%spherical_cutoff, particles_qm)
1278 : ! Possibly added charges
1279 70 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
1280 : CALL build_mm_pot(qpot(iatom), nonperiodic_pot_type, eta_a(0), &
1281 : qmmm_env%added_charges%potentials, &
1282 : qmmm_env%added_charges%added_particles, &
1283 : qmmm_env%added_charges%mm_atom_chrg, &
1284 : qmmm_env%added_charges%mm_atom_index, &
1285 : mm_cell, iatom, qmmm_env%spherical_cutoff, &
1286 0 : particles_qm)
1287 : CALL build_mm_dpot(mcharge(iatom), nonperiodic_pot_type, eta_a(0), &
1288 : qmmm_env%added_charges%potentials, &
1289 : qmmm_env%added_charges%added_particles, &
1290 : qmmm_env%added_charges%mm_atom_chrg, &
1291 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
1292 : Forces_added_charges, &
1293 0 : Forces_QM(:, iqm), qmmm_env%spherical_cutoff, particles_qm)
1294 : END IF
1295 : END DO
1296 : END DO
1297 : CASE DEFAULT
1298 40 : CPABORT("Unknown Ewald type!")
1299 : END SELECT
1300 :
1301 : ! Transfer QM gradients to the QM particles..
1302 40 : iqm = 0
1303 120 : DO ikind = 1, nkind
1304 80 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1305 80 : IF (do_dftb) THEN
1306 40 : NULLIFY (dftb_kind)
1307 40 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1308 40 : CALL get_dftb_atom_param(dftb_kind, defined=defined, natorb=natorb)
1309 40 : IF (.NOT. defined .OR. natorb < 1) CYCLE
1310 : ELSE IF (do_xtb) THEN
1311 : !
1312 : END IF
1313 240 : DO i = 1, SIZE(list)
1314 120 : iqm = iqm + 1
1315 120 : iatom = qmmm_env%qm_atom_index(list(i))
1316 1040 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
1317 : END DO
1318 : END DO
1319 :
1320 : ! derivatives from qm charges
1321 520 : Forces_QM = 0.0_dp
1322 40 : IF (SIZE(matrix_p) == 2) THEN
1323 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
1324 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1325 : END IF
1326 : !
1327 40 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
1328 160 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1329 120 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock)
1330 : !
1331 120 : IF (iatom == jatom) CYCLE
1332 : !
1333 60 : gmij = -0.5_dp*(qpot(iatom) + qpot(jatom))
1334 60 : NULLIFY (pblock)
1335 : CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
1336 60 : row=iatom, col=jatom, block=pblock, found=found)
1337 60 : CPASSERT(found)
1338 280 : DO i = 1, 3
1339 180 : NULLIFY (dsblock)
1340 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
1341 180 : row=iatom, col=jatom, block=dsblock, found=found)
1342 180 : CPASSERT(found)
1343 1584 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
1344 180 : Forces_QM(i, iatom) = Forces_QM(i, iatom) + fi
1345 480 : Forces_QM(i, jatom) = Forces_QM(i, jatom) - fi
1346 : END DO
1347 : END DO
1348 40 : CALL dbcsr_iterator_stop(iter)
1349 : !
1350 40 : IF (SIZE(matrix_p) == 2) THEN
1351 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
1352 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1353 : END IF
1354 : !
1355 : ! Transfer QM gradients to the QM particles..
1356 1000 : CALL para_env%sum(Forces_QM)
1357 120 : DO ikind = 1, nkind
1358 80 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1359 240 : DO i = 1, SIZE(list)
1360 120 : iqm = list(i)
1361 120 : iatom = qmmm_env%qm_atom_index(iqm)
1362 1040 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
1363 : END DO
1364 : END DO
1365 : !
1366 40 : DEALLOCATE (mcharge)
1367 : !
1368 : ! MM forces will be handled directly from the QMMM module in the same way
1369 : ! as for GPW/GAPW methods
1370 40 : DEALLOCATE (Forces_QM)
1371 40 : DEALLOCATE (qpot)
1372 :
1373 : ! Release Ewald environment
1374 40 : CALL ewald_env_release(ewald_env)
1375 40 : DEALLOCATE (ewald_env)
1376 40 : CALL ewald_pw_release(ewald_pw)
1377 40 : DEALLOCATE (ewald_pw)
1378 :
1379 160 : CALL dbcsr_deallocate_matrix_set(matrix_s)
1380 :
1381 : END IF
1382 :
1383 1336 : CALL timestop(handle)
1384 :
1385 1336 : END SUBROUTINE deriv_tb_qmmm_matrix_smeared
1386 :
1387 : ! **************************************************************************************************
1388 : !> \brief ...
1389 : !> \param qpot ...
1390 : !> \param pot_type ...
1391 : !> \param qm_alpha ...
1392 : !> \param potentials ...
1393 : !> \param particles_mm ...
1394 : !> \param mm_charges ...
1395 : !> \param mm_atom_index ...
1396 : !> \param mm_cell ...
1397 : !> \param IndQM ...
1398 : !> \param qmmm_spherical_cutoff ...
1399 : !> \param particles_qm ...
1400 : ! **************************************************************************************************
1401 7944 : SUBROUTINE build_mm_pot(qpot, pot_type, qm_alpha, potentials, &
1402 : particles_mm, mm_charges, mm_atom_index, mm_cell, IndQM, &
1403 : qmmm_spherical_cutoff, particles_qm)
1404 :
1405 : REAL(KIND=dp), INTENT(INOUT) :: qpot
1406 : INTEGER, INTENT(IN) :: pot_type
1407 : REAL(KIND=dp), INTENT(IN) :: qm_alpha
1408 : TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1409 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
1410 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges
1411 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1412 : TYPE(cell_type), POINTER :: mm_cell
1413 : INTEGER, INTENT(IN) :: IndQM
1414 : REAL(KIND=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
1415 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
1416 :
1417 : CHARACTER(len=*), PARAMETER :: routineN = 'build_mm_pot'
1418 : REAL(KIND=dp), PARAMETER :: qsmall = 1.0e-15_dp
1419 :
1420 : INTEGER :: handle, Imm, Imp, IndMM, Ipot
1421 : REAL(KIND=dp) :: dr, qeff, rt1, rt2, rt3, &
1422 : sph_chrg_factor, sr
1423 : REAL(KIND=dp), DIMENSION(3) :: r_pbc, rij
1424 : TYPE(qmmm_pot_type), POINTER :: Pot
1425 :
1426 7944 : CALL timeset(routineN, handle)
1427 : ! Loop Over MM atoms
1428 : ! Loop over Pot stores atoms with the same charge
1429 23616 : MainLoopPot: DO Ipot = 1, SIZE(Potentials)
1430 15672 : Pot => Potentials(Ipot)%Pot
1431 : ! Loop over atoms belonging to this type
1432 70632 : LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index)
1433 47016 : Imm = Pot%mm_atom_index(Imp)
1434 47016 : IndMM = mm_atom_index(Imm)
1435 188064 : r_pbc = pbc(particles_mm(IndMM)%r - particles_qm(IndQM)%r, mm_cell)
1436 47016 : rt1 = r_pbc(1)
1437 47016 : rt2 = r_pbc(2)
1438 47016 : rt3 = r_pbc(3)
1439 188064 : rij = [rt1, rt2, rt3]
1440 188064 : dr = SQRT(SUM(rij**2))
1441 47016 : qeff = mm_charges(Imm)
1442 : ! Computes the screening factor for the spherical cutoff (if defined)
1443 47016 : IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1444 28656 : CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
1445 28656 : qeff = qeff*sph_chrg_factor
1446 : END IF
1447 47016 : IF (ABS(qeff) <= qsmall) CYCLE LoopMM
1448 62688 : IF (dr > rtiny) THEN
1449 47016 : IF (pot_type == pot_tb_nonperiodic) THEN
1450 16344 : sr = gamma_rab_sr(dr, qm_alpha, eta_mm, 0.0_dp)
1451 16344 : qpot = qpot + qeff*(1.0_dp/dr - sr)
1452 : ELSE IF (pot_type == pot_tb_short_range) THEN
1453 12312 : sr = gamma_rab_sr(dr, qm_alpha, eta_mm, 0.0_dp)
1454 12312 : qpot = qpot - qeff*sr
1455 : ELSE IF (pot_type == pot_ewald_short_range) THEN
1456 14328 : sr = erfc(qm_alpha*dr)/dr
1457 14328 : qpot = qpot + qeff*sr
1458 : ELSE IF (pot_type == pot_gauss_nonperiodic) THEN
1459 2016 : sr = erf(dr/Pot%Rc)/dr
1460 2016 : qpot = qpot + qeff*sr
1461 : ELSE IF (pot_type == pot_gauss_short_range) THEN
1462 2016 : sr = erfc(dr/Pot%Rc)/dr
1463 2016 : qpot = qpot - qeff*sr
1464 : ELSE
1465 0 : CPABORT("Unknown pot_type for dr > rtiny")
1466 : END IF
1467 0 : ELSE IF (pot_type == pot_gauss_nonperiodic) THEN
1468 0 : qpot = qpot + qeff*2.0_dp/(rootpi*Pot%Rc)
1469 : END IF
1470 : END DO LoopMM
1471 : END DO MainLoopPot
1472 7944 : CALL timestop(handle)
1473 7944 : END SUBROUTINE build_mm_pot
1474 :
1475 : ! **************************************************************************************************
1476 : !> \brief ...
1477 : !> \param qcharge ...
1478 : !> \param pot_type ...
1479 : !> \param qm_alpha ...
1480 : !> \param potentials ...
1481 : !> \param particles_mm ...
1482 : !> \param mm_charges ...
1483 : !> \param mm_atom_index ...
1484 : !> \param mm_cell ...
1485 : !> \param IndQM ...
1486 : !> \param forces ...
1487 : !> \param forces_qm ...
1488 : !> \param qmmm_spherical_cutoff ...
1489 : !> \param particles_qm ...
1490 : ! **************************************************************************************************
1491 492 : SUBROUTINE build_mm_dpot(qcharge, pot_type, qm_alpha, potentials, &
1492 : particles_mm, mm_charges, mm_atom_index, mm_cell, IndQM, &
1493 246 : forces, forces_qm, qmmm_spherical_cutoff, particles_qm)
1494 :
1495 : REAL(KIND=dp), INTENT(IN) :: qcharge
1496 : INTEGER, INTENT(IN) :: pot_type
1497 : REAL(KIND=dp), INTENT(IN) :: qm_alpha
1498 : TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1499 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
1500 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges
1501 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1502 : TYPE(cell_type), POINTER :: mm_cell
1503 : INTEGER, INTENT(IN) :: IndQM
1504 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: forces
1505 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: forces_qm
1506 : REAL(KIND=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
1507 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
1508 :
1509 : CHARACTER(len=*), PARAMETER :: routineN = 'build_mm_dpot'
1510 : REAL(KIND=dp), PARAMETER :: qsmall = 1.0e-15_dp
1511 :
1512 : INTEGER :: handle, Imm, Imp, IndMM, Ipot
1513 : REAL(KIND=dp) :: dr, drm, drp, dsr, fsr, qeff, rt1, rt2, &
1514 : rt3, sph_chrg_factor
1515 : REAL(KIND=dp), DIMENSION(3) :: force_ab, r_pbc, rij
1516 : TYPE(qmmm_pot_type), POINTER :: Pot
1517 :
1518 246 : CALL timeset(routineN, handle)
1519 : ! Loop Over MM atoms
1520 : ! Loop over Pot stores atoms with the same charge
1521 630 : MainLoopPot: DO Ipot = 1, SIZE(Potentials)
1522 384 : Pot => Potentials(Ipot)%Pot
1523 : ! Loop over atoms belonging to this type
1524 1782 : LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index)
1525 1152 : Imm = Pot%mm_atom_index(Imp)
1526 1152 : IndMM = mm_atom_index(Imm)
1527 4608 : r_pbc = pbc(particles_mm(IndMM)%r - particles_qm(IndQM)%r, mm_cell)
1528 1152 : rt1 = r_pbc(1)
1529 1152 : rt2 = r_pbc(2)
1530 1152 : rt3 = r_pbc(3)
1531 4608 : rij = [rt1, rt2, rt3]
1532 4608 : dr = SQRT(SUM(rij**2))
1533 1152 : qeff = mm_charges(Imm)
1534 : ! Computes the screening factor for the spherical cutoff (if defined)
1535 : ! We neglect derivative of cutoff function for gradients!!!
1536 1152 : IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1537 720 : CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
1538 720 : qeff = qeff*sph_chrg_factor
1539 : END IF
1540 1152 : IF (ABS(qeff) <= qsmall) CYCLE LoopMM
1541 1152 : IF (dr > rtiny) THEN
1542 1152 : drp = dr + ddrmm
1543 1152 : drm = dr - ddrmm
1544 1152 : IF (pot_type == pot_tb_nonperiodic) THEN
1545 : dsr = 0.5_dp*(gamma_rab_sr(drp, qm_alpha, eta_mm, 0.0_dp) - &
1546 396 : gamma_rab_sr(drm, qm_alpha, eta_mm, 0.0_dp))/ddrmm
1547 396 : fsr = qeff*qcharge*(-1.0_dp/(dr*dr) - dsr)
1548 : ELSE IF (pot_type == pot_tb_short_range) THEN
1549 : dsr = 0.5_dp*(gamma_rab_sr(drp, qm_alpha, eta_mm, 0.0_dp) - &
1550 324 : gamma_rab_sr(drm, qm_alpha, eta_mm, 0.0_dp))/ddrmm
1551 324 : fsr = -qeff*qcharge*dsr
1552 : ELSE IF (pot_type == pot_ewald_short_range) THEN
1553 360 : dsr = 0.5_dp*(erfc(qm_alpha*drp)/drp - erfc(qm_alpha*drm)/drm)/ddrmm
1554 360 : fsr = qeff*qcharge*dsr
1555 : ELSE IF (pot_type == pot_gauss_nonperiodic) THEN
1556 : fsr = qeff*qcharge*(2.0_dp*EXP(-(dr/Pot%Rc)**2)/(rootpi*Pot%Rc*dr) - &
1557 36 : erf(dr/Pot%Rc)/dr**2)
1558 : ELSE IF (pot_type == pot_gauss_short_range) THEN
1559 : fsr = qeff*qcharge*(2.0_dp*EXP(-(dr/Pot%Rc)**2)/(rootpi*Pot%Rc*dr) + &
1560 36 : erfc(dr/Pot%Rc)/dr**2)
1561 : ELSE
1562 0 : CPABORT("Unknown pot_type for dr > rtiny")
1563 : END IF
1564 4608 : force_ab = -fsr*rij/dr
1565 : ELSE
1566 0 : force_ab = 0.0_dp
1567 : END IF
1568 : ! The array of QM forces are really the forces
1569 4608 : forces_qm(:) = forces_qm(:) - force_ab
1570 : ! The one of MM atoms are instead gradients
1571 4992 : forces(:, Imm) = forces(:, Imm) - force_ab
1572 : END DO LoopMM
1573 : END DO MainLoopPot
1574 :
1575 246 : CALL timestop(handle)
1576 :
1577 246 : END SUBROUTINE build_mm_dpot
1578 :
1579 : END MODULE qmmm_tb_methods
|
