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 1114 : 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 1114 : gaussian=.FALSE.)
327 :
328 1114 : 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 1334 : 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 1334 : 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 1334 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges_mm, qpot
383 1334 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hblock, sblock
384 1334 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
385 : TYPE(dbcsr_iterator_type) :: iter
386 1334 : 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 1334 : POINTER :: sab_nl
393 1334 : 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 1334 : 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 1334 : CALL timeset(routineN, handle)
406 :
407 1334 : IF (gaussian) THEN
408 220 : nonperiodic_pot_type = pot_gauss_nonperiodic
409 220 : short_range_pot_type = pot_gauss_short_range
410 : ELSE
411 1114 : nonperiodic_pot_type = pot_tb_nonperiodic
412 1114 : 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 1334 : natom=natom)
422 1334 : dftb_control => dft_control%qs_control%dftb_control
423 1334 : xtb_control => dft_control%qs_control%xtb_control
424 :
425 1334 : IF (dft_control%qs_control%dftb) THEN
426 : do_dftb = .TRUE.
427 : do_xtb = .FALSE.
428 666 : 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 1334 : CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
436 :
437 1334 : NULLIFY (matrix_s)
438 1334 : IF (do_dftb) THEN
439 668 : CALL build_dftb_overlap(qs_env, 0, matrix_s)
440 666 : ELSE IF (do_xtb) THEN
441 666 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
442 666 : 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 4002 : ALLOCATE (qpot(natom))
446 5336 : qpot = 0.0_dp
447 1334 : pc_ener = 0.0_dp
448 :
449 : ! Create Ewald environments
450 1334 : poisson_section => section_vals_get_subs_vals(qs_env%input, "MM%POISSON")
451 24012 : ALLOCATE (ewald_env)
452 1334 : CALL ewald_env_create(ewald_env, para_env)
453 1334 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
454 1334 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
455 1334 : CALL read_ewald_section(ewald_env, ewald_section)
456 1334 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
457 1334 : ALLOCATE (ewald_pw)
458 1334 : CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=print_section)
459 :
460 1334 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, do_ipol=do_ipol)
461 1334 : IF (do_multipoles) CPABORT("No multipole force fields allowed in TB QM/MM")
462 1334 : 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 2330 : DO ipot = 1, SIZE(qmmm_env%Potentials)
469 1550 : Pot => qmmm_env%Potentials(ipot)%Pot
470 1550 : nmm = SIZE(Pot%mm_atom_index)
471 : ! get a 'clean' mm particle set
472 1550 : NULLIFY (atoms_mm)
473 1550 : CALL allocate_particle_set(atoms_mm, nmm)
474 4650 : ALLOCATE (charges_mm(nmm))
475 6200 : DO Imp = 1, nmm
476 4650 : Imm = Pot%mm_atom_index(Imp)
477 4650 : IndMM = qmmm_env%mm_atom_index(Imm)
478 37200 : atoms_mm(Imp)%r = particles_mm(IndMM)%r
479 4650 : atoms_mm(Imp)%atomic_kind => particles_mm(IndMM)%atomic_kind
480 6200 : charges_mm(Imp) = qmmm_env%mm_atom_chrg(Imm)
481 : END DO
482 1550 : IF (ewald_type == do_ewald_ewald) THEN
483 0 : CPABORT("Ewald not implemented for TB/QMMM")
484 1550 : ELSE IF (ewald_type == do_ewald_spme) THEN
485 : ! spme electrostatic potential
486 1550 : CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, particles_qm, qpot)
487 : END IF
488 1550 : CALL deallocate_particle_set(atoms_mm)
489 2330 : DEALLOCATE (charges_mm)
490 : END DO
491 780 : 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 780 : DEALLOCATE (charges_mm)
514 : END DO
515 : END IF
516 5460 : 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 780 : CALL ewald_env_get(ewald_env, rcut=rcutoff(1), alpha=alpha)
521 780 : rcutoff(2) = 0.025_dp*rcutoff(1)
522 780 : rcutoff(1) = 2.0_dp*rcutoff(1)
523 780 : nkind = SIZE(atomic_kind_set)
524 2340 : DO ikind = 1, nkind
525 1560 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
526 1560 : 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 888 : ELSE IF (do_xtb) THEN
535 888 : NULLIFY (xtb_kind)
536 888 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
537 888 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
538 888 : eta_a(0) = eta_mm
539 : END IF
540 4680 : DO i = 1, SIZE(list)
541 2340 : 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 2340 : 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 2340 : particles_qm)
550 : ! Possibly added charges
551 2340 : 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 3900 : 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 1334 : CPABORT("Unknown Ewald type!")
608 : END SELECT
609 :
610 : ! Allocate the core Hamiltonian matrix
611 1334 : CALL get_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env_loc)
612 1334 : matrix_h => ks_qmmm_env_loc%matrix_h
613 1334 : CALL dbcsr_allocate_matrix_set(matrix_h, 1)
614 1334 : ALLOCATE (matrix_h(1)%matrix)
615 : CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, &
616 1334 : name="QMMM HAMILTONIAN MATRIX")
617 1334 : CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
618 :
619 1334 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
620 5336 : DO WHILE (dbcsr_iterator_blocks_left(iter))
621 4002 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock)
622 4002 : NULLIFY (hblock)
623 : CALL dbcsr_get_block_p(matrix=matrix_h(1)%matrix, &
624 4002 : row=iatom, col=jatom, block=hblock, found=found)
625 4002 : CPASSERT(found)
626 76486 : hblock = hblock - 0.5_dp*sblock*(qpot(iatom) + qpot(jatom))
627 : END DO
628 1334 : CALL dbcsr_iterator_stop(iter)
629 :
630 1334 : ks_qmmm_env_loc%matrix_h => matrix_h
631 1334 : ks_qmmm_env_loc%pc_ener = pc_ener
632 :
633 1334 : DEALLOCATE (qpot)
634 :
635 : ! Release Ewald environment
636 1334 : CALL ewald_env_release(ewald_env)
637 1334 : DEALLOCATE (ewald_env)
638 1334 : CALL ewald_pw_release(ewald_pw)
639 1334 : DEALLOCATE (ewald_pw)
640 :
641 1334 : CALL dbcsr_deallocate_matrix_set(matrix_s)
642 :
643 1334 : CALL timestop(handle)
644 :
645 5336 : 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 1114 : 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 1114 : calc_force, Forces, Forces_added_charges, gaussian=.FALSE.)
913 :
914 1114 : 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 1334 : 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 1334 : 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 1334 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges_mm, mcharge, qpot
979 1334 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges, dsblock, Forces_MM, Forces_QM, &
980 1334 : pblock, sblock
981 1334 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
982 : TYPE(dbcsr_iterator_type) :: iter
983 1334 : 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 1334 : POINTER :: sab_nl
990 1334 : 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 1334 : 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 1334 : CALL timeset(routineN, handle)
1003 1334 : IF (gaussian) THEN
1004 220 : nonperiodic_pot_type = pot_gauss_nonperiodic
1005 220 : short_range_pot_type = pot_gauss_short_range
1006 : ELSE
1007 1114 : nonperiodic_pot_type = pot_tb_nonperiodic
1008 1114 : short_range_pot_type = pot_tb_short_range
1009 : END IF
1010 1334 : IF (calc_force) THEN
1011 38 : 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 38 : natom=number_qm_atoms)
1020 38 : dftb_control => dft_control%qs_control%dftb_control
1021 38 : xtb_control => dft_control%qs_control%xtb_control
1022 :
1023 38 : IF (dft_control%qs_control%dftb) THEN
1024 : do_dftb = .TRUE.
1025 : do_xtb = .FALSE.
1026 18 : 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 38 : 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 18 : 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 18 : basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
1040 : END IF
1041 38 : CALL qs_rho_get(rho, rho_ao=matrix_p)
1042 :
1043 38 : nspins = dft_control%nspins
1044 38 : nkind = SIZE(atomic_kind_set)
1045 : ! Mulliken charges
1046 152 : ALLOCATE (charges(number_qm_atoms, nspins))
1047 : !
1048 38 : CALL mulliken_charges(matrix_p, matrix_s(1)%matrix, para_env, charges)
1049 : !
1050 114 : ALLOCATE (mcharge(number_qm_atoms))
1051 114 : DO ikind = 1, nkind
1052 76 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
1053 76 : 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 36 : ELSE IF (do_xtb) THEN
1057 36 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1058 36 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
1059 : END IF
1060 304 : DO iatom = 1, natom
1061 114 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1062 304 : mcharge(atom_a) = zeff - SUM(charges(atom_a, 1:nspins))
1063 : END DO
1064 : END DO
1065 38 : DEALLOCATE (charges)
1066 :
1067 114 : ALLOCATE (qpot(number_qm_atoms))
1068 152 : qpot = 0.0_dp
1069 114 : ALLOCATE (Forces_QM(3, number_qm_atoms))
1070 494 : Forces_QM = 0.0_dp
1071 :
1072 : ! Create Ewald environments
1073 38 : poisson_section => section_vals_get_subs_vals(qs_env%input, "MM%POISSON")
1074 684 : ALLOCATE (ewald_env)
1075 38 : CALL ewald_env_create(ewald_env, para_env)
1076 38 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
1077 38 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
1078 38 : CALL read_ewald_section(ewald_env, ewald_section)
1079 38 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
1080 38 : ALLOCATE (ewald_pw)
1081 38 : CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=print_section)
1082 :
1083 38 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, do_ipol=do_ipol)
1084 38 : IF (do_multipoles) CPABORT("No multipole force fields allowed in DFTB QM/MM")
1085 38 : IF (do_ipol /= do_fist_pol_none) THEN
1086 0 : CPABORT("No polarizable force fields allowed in DFTB QM/MM")
1087 : END IF
1088 :
1089 0 : SELECT CASE (ewald_type)
1090 : CASE (do_ewald_pme)
1091 0 : CPABORT("PME Ewald type not implemented for DFTB/QMMM")
1092 : CASE (do_ewald_ewald, do_ewald_spme)
1093 62 : DO ipot = 1, SIZE(qmmm_env%Potentials)
1094 38 : Pot => qmmm_env%Potentials(ipot)%Pot
1095 38 : nmm = SIZE(Pot%mm_atom_index)
1096 : ! get a 'clean' mm particle set
1097 38 : NULLIFY (atoms_mm)
1098 38 : CALL allocate_particle_set(atoms_mm, nmm)
1099 114 : ALLOCATE (charges_mm(nmm))
1100 152 : DO Imp = 1, nmm
1101 114 : Imm = Pot%mm_atom_index(Imp)
1102 114 : IndMM = qmmm_env%mm_atom_index(Imm)
1103 912 : atoms_mm(Imp)%r = particles_mm(IndMM)%r
1104 114 : atoms_mm(Imp)%atomic_kind => particles_mm(IndMM)%atomic_kind
1105 152 : charges_mm(Imp) = qmmm_env%mm_atom_chrg(Imm)
1106 : END DO
1107 : ! force array for mm atoms
1108 114 : ALLOCATE (Forces_MM(3, nmm))
1109 494 : Forces_MM = 0.0_dp
1110 38 : IF (ewald_type == do_ewald_ewald) THEN
1111 0 : CPABORT("Ewald not implemented for DFTB/QMMM")
1112 38 : ELSE IF (ewald_type == do_ewald_spme) THEN
1113 : ! spme electrostatic potential
1114 : CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
1115 38 : particles_qm, qpot)
1116 : ! forces QM
1117 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
1118 38 : particles_qm, mcharge, Forces_QM)
1119 : ! forces MM
1120 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, particles_qm, mcharge, &
1121 38 : atoms_mm, charges_mm, Forces_MM)
1122 : END IF
1123 38 : CALL deallocate_particle_set(atoms_mm)
1124 38 : DEALLOCATE (charges_mm)
1125 : ! transfer MM forces
1126 950 : CALL para_env%sum(Forces_MM)
1127 152 : DO Imp = 1, nmm
1128 114 : Imm = Pot%mm_atom_index(Imp)
1129 950 : Forces(:, Imm) = Forces(:, Imm) - Forces_MM(:, Imp)
1130 : END DO
1131 62 : DEALLOCATE (Forces_MM)
1132 : END DO
1133 :
1134 24 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
1135 0 : DO ipot = 1, SIZE(qmmm_env%added_charges%Potentials)
1136 0 : Pot => qmmm_env%added_charges%Potentials(ipot)%Pot
1137 0 : nmm = SIZE(Pot%mm_atom_index)
1138 : ! get a 'clean' mm particle set
1139 0 : NULLIFY (atoms_mm)
1140 0 : CALL allocate_particle_set(atoms_mm, nmm)
1141 0 : ALLOCATE (charges_mm(nmm))
1142 0 : DO Imp = 1, nmm
1143 0 : Imm = Pot%mm_atom_index(Imp)
1144 0 : IndMM = qmmm_env%added_charges%mm_atom_index(Imm)
1145 0 : atoms_mm(Imp)%r = qmmm_env%added_charges%added_particles(IndMM)%r
1146 0 : atoms_mm(Imp)%atomic_kind => qmmm_env%added_charges%added_particles(IndMM)%atomic_kind
1147 0 : charges_mm(Imp) = qmmm_env%added_charges%mm_atom_chrg(Imm)
1148 : END DO
1149 : ! force array for mm atoms
1150 0 : ALLOCATE (Forces_MM(3, nmm))
1151 0 : Forces_MM = 0.0_dp
1152 0 : IF (ewald_type == do_ewald_ewald) THEN
1153 0 : CPABORT("Ewald not implemented for DFTB/QMMM")
1154 0 : ELSE IF (ewald_type == do_ewald_spme) THEN
1155 : ! spme electrostatic potential
1156 : CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, &
1157 0 : charges_mm, particles_qm, qpot)
1158 : ! forces QM
1159 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
1160 0 : particles_qm, mcharge, Forces_QM)
1161 : ! forces MM
1162 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, particles_qm, mcharge, &
1163 0 : atoms_mm, charges_mm, Forces_MM)
1164 : END IF
1165 0 : CALL deallocate_particle_set(atoms_mm)
1166 : ! transfer MM forces
1167 0 : CALL para_env%sum(Forces_MM)
1168 0 : DO Imp = 1, nmm
1169 0 : Imm = Pot%mm_atom_index(Imp)
1170 0 : Forces_added_charges(:, Imm) = Forces_added_charges(:, Imm) - Forces_MM(:, Imp)
1171 : END DO
1172 24 : DEALLOCATE (Forces_MM)
1173 : END DO
1174 : END IF
1175 168 : CALL para_env%sum(qpot)
1176 600 : CALL para_env%sum(Forces_QM)
1177 : ! Add the Ewald real-space term and the method-specific short-range correction
1178 : ! This is effectively using a minimum image convention!
1179 : ! Set rcutoff to values compatible with alpha Ewald
1180 24 : CALL ewald_env_get(ewald_env, rcut=rcutoff(1), alpha=alpha)
1181 24 : rcutoff(2) = 0.025_dp*rcutoff(1)
1182 24 : rcutoff(1) = 2.0_dp*rcutoff(1)
1183 24 : nkind = SIZE(atomic_kind_set)
1184 24 : iqm = 0
1185 72 : DO ikind = 1, nkind
1186 48 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1187 48 : IF (do_dftb) THEN
1188 24 : NULLIFY (dftb_kind)
1189 24 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1190 : CALL get_dftb_atom_param(dftb_kind, &
1191 24 : defined=defined, eta=eta_a, natorb=natorb)
1192 : ! use mm charge smearing for non-scc cases
1193 24 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
1194 24 : IF (.NOT. defined .OR. natorb < 1) CYCLE
1195 24 : ELSE IF (do_xtb) THEN
1196 24 : eta_a(0) = eta_mm
1197 : END IF
1198 144 : DO i = 1, SIZE(list)
1199 72 : iatom = list(i)
1200 72 : iqm = iqm + 1
1201 : CALL build_mm_pot(qpot(iatom), short_range_pot_type, eta_a(0), &
1202 : qmmm_env%Potentials, particles_mm, &
1203 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1204 72 : mm_cell, iatom, rcutoff, particles_qm)
1205 : CALL build_mm_dpot(mcharge(iatom), short_range_pot_type, eta_a(0), &
1206 : qmmm_env%Potentials, particles_mm, &
1207 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1208 : mm_cell, iatom, Forces, Forces_QM(:, iqm), &
1209 72 : rcutoff, particles_qm)
1210 : CALL build_mm_pot(qpot(iatom), pot_ewald_short_range, alpha, &
1211 : qmmm_env%Potentials, particles_mm, &
1212 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1213 72 : mm_cell, iatom, rcutoff, particles_qm)
1214 : CALL build_mm_dpot(mcharge(iatom), pot_ewald_short_range, alpha, &
1215 : qmmm_env%Potentials, particles_mm, &
1216 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1217 : mm_cell, iatom, Forces, Forces_QM(:, iqm), &
1218 72 : rcutoff, particles_qm)
1219 : ! Possibly added charges
1220 120 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
1221 : CALL build_mm_pot(qpot(iatom), short_range_pot_type, eta_a(0), &
1222 : qmmm_env%added_charges%potentials, &
1223 : qmmm_env%added_charges%added_particles, &
1224 : qmmm_env%added_charges%mm_atom_chrg, &
1225 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, rcutoff, &
1226 0 : particles_qm)
1227 : CALL build_mm_dpot( &
1228 : mcharge(iatom), short_range_pot_type, eta_a(0), &
1229 : qmmm_env%added_charges%potentials, &
1230 : qmmm_env%added_charges%added_particles, qmmm_env%added_charges%mm_atom_chrg, &
1231 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
1232 : Forces_added_charges, Forces_QM(:, iqm), &
1233 0 : rcutoff, particles_qm)
1234 : CALL build_mm_pot(qpot(iatom), pot_ewald_short_range, alpha, &
1235 : qmmm_env%added_charges%potentials, &
1236 : qmmm_env%added_charges%added_particles, &
1237 : qmmm_env%added_charges%mm_atom_chrg, &
1238 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, rcutoff, &
1239 0 : particles_qm)
1240 : CALL build_mm_dpot( &
1241 : mcharge(iatom), pot_ewald_short_range, alpha, &
1242 : qmmm_env%added_charges%potentials, &
1243 : qmmm_env%added_charges%added_particles, qmmm_env%added_charges%mm_atom_chrg, &
1244 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
1245 : Forces_added_charges, Forces_QM(:, iqm), &
1246 0 : rcutoff, particles_qm)
1247 : END IF
1248 : END DO
1249 : END DO
1250 :
1251 : CASE (do_ewald_none)
1252 : ! Directly sum the nonperiodic regularized electrostatic potential
1253 : ! calculate potential and forces from classical charges
1254 : iqm = 0
1255 42 : DO ikind = 1, nkind
1256 28 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1257 28 : IF (do_dftb) THEN
1258 16 : NULLIFY (dftb_kind)
1259 16 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1260 : CALL get_dftb_atom_param(dftb_kind, &
1261 16 : defined=defined, eta=eta_a, natorb=natorb)
1262 : ! use mm charge smearing for non-scc cases
1263 16 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
1264 16 : IF (.NOT. defined .OR. natorb < 1) CYCLE
1265 12 : ELSE IF (do_xtb) THEN
1266 12 : eta_a(0) = eta_mm
1267 : END IF
1268 84 : DO i = 1, SIZE(list)
1269 42 : iatom = list(i)
1270 42 : iqm = iqm + 1
1271 : CALL build_mm_pot(qpot(iatom), nonperiodic_pot_type, eta_a(0), &
1272 : qmmm_env%Potentials, particles_mm, &
1273 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, &
1274 42 : qmmm_env%spherical_cutoff, particles_qm)
1275 : CALL build_mm_dpot(mcharge(iatom), nonperiodic_pot_type, eta_a(0), &
1276 : qmmm_env%Potentials, particles_mm, &
1277 : qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1278 : mm_cell, iatom, Forces, Forces_QM(:, iqm), &
1279 42 : qmmm_env%spherical_cutoff, particles_qm)
1280 : ! Possibly added charges
1281 70 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
1282 : CALL build_mm_pot(qpot(iatom), nonperiodic_pot_type, eta_a(0), &
1283 : qmmm_env%added_charges%potentials, &
1284 : qmmm_env%added_charges%added_particles, &
1285 : qmmm_env%added_charges%mm_atom_chrg, &
1286 : qmmm_env%added_charges%mm_atom_index, &
1287 : mm_cell, iatom, qmmm_env%spherical_cutoff, &
1288 0 : particles_qm)
1289 : CALL build_mm_dpot(mcharge(iatom), nonperiodic_pot_type, eta_a(0), &
1290 : qmmm_env%added_charges%potentials, &
1291 : qmmm_env%added_charges%added_particles, &
1292 : qmmm_env%added_charges%mm_atom_chrg, &
1293 : qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
1294 : Forces_added_charges, &
1295 0 : Forces_QM(:, iqm), qmmm_env%spherical_cutoff, particles_qm)
1296 : END IF
1297 : END DO
1298 : END DO
1299 : CASE DEFAULT
1300 38 : CPABORT("Unknown Ewald type!")
1301 : END SELECT
1302 :
1303 : ! Transfer QM gradients to the QM particles..
1304 38 : iqm = 0
1305 114 : DO ikind = 1, nkind
1306 76 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1307 76 : IF (do_dftb) THEN
1308 40 : NULLIFY (dftb_kind)
1309 40 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1310 40 : CALL get_dftb_atom_param(dftb_kind, defined=defined, natorb=natorb)
1311 40 : IF (.NOT. defined .OR. natorb < 1) CYCLE
1312 : ELSE IF (do_xtb) THEN
1313 : !
1314 : END IF
1315 228 : DO i = 1, SIZE(list)
1316 114 : iqm = iqm + 1
1317 114 : iatom = qmmm_env%qm_atom_index(list(i))
1318 988 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
1319 : END DO
1320 : END DO
1321 :
1322 : ! derivatives from qm charges
1323 494 : Forces_QM = 0.0_dp
1324 38 : IF (SIZE(matrix_p) == 2) THEN
1325 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
1326 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1327 : END IF
1328 : !
1329 38 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
1330 152 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1331 114 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock)
1332 : !
1333 114 : IF (iatom == jatom) CYCLE
1334 : !
1335 57 : gmij = -0.5_dp*(qpot(iatom) + qpot(jatom))
1336 57 : NULLIFY (pblock)
1337 : CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
1338 57 : row=iatom, col=jatom, block=pblock, found=found)
1339 57 : CPASSERT(found)
1340 266 : DO i = 1, 3
1341 171 : NULLIFY (dsblock)
1342 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
1343 171 : row=iatom, col=jatom, block=dsblock, found=found)
1344 171 : CPASSERT(found)
1345 1485 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
1346 171 : Forces_QM(i, iatom) = Forces_QM(i, iatom) + fi
1347 456 : Forces_QM(i, jatom) = Forces_QM(i, jatom) - fi
1348 : END DO
1349 : END DO
1350 38 : CALL dbcsr_iterator_stop(iter)
1351 : !
1352 38 : IF (SIZE(matrix_p) == 2) THEN
1353 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
1354 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1355 : END IF
1356 : !
1357 : ! Transfer QM gradients to the QM particles..
1358 950 : CALL para_env%sum(Forces_QM)
1359 114 : DO ikind = 1, nkind
1360 76 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1361 228 : DO i = 1, SIZE(list)
1362 114 : iqm = list(i)
1363 114 : iatom = qmmm_env%qm_atom_index(iqm)
1364 988 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
1365 : END DO
1366 : END DO
1367 : !
1368 38 : DEALLOCATE (mcharge)
1369 : !
1370 : ! MM forces will be handled directly from the QMMM module in the same way
1371 : ! as for GPW/GAPW methods
1372 38 : DEALLOCATE (Forces_QM)
1373 38 : DEALLOCATE (qpot)
1374 :
1375 : ! Release Ewald environment
1376 38 : CALL ewald_env_release(ewald_env)
1377 38 : DEALLOCATE (ewald_env)
1378 38 : CALL ewald_pw_release(ewald_pw)
1379 38 : DEALLOCATE (ewald_pw)
1380 :
1381 152 : CALL dbcsr_deallocate_matrix_set(matrix_s)
1382 :
1383 : END IF
1384 :
1385 1334 : CALL timestop(handle)
1386 :
1387 1334 : END SUBROUTINE deriv_tb_qmmm_matrix_smeared
1388 :
1389 : ! **************************************************************************************************
1390 : !> \brief ...
1391 : !> \param qpot ...
1392 : !> \param pot_type ...
1393 : !> \param qm_alpha ...
1394 : !> \param potentials ...
1395 : !> \param particles_mm ...
1396 : !> \param mm_charges ...
1397 : !> \param mm_atom_index ...
1398 : !> \param mm_cell ...
1399 : !> \param IndQM ...
1400 : !> \param qmmm_spherical_cutoff ...
1401 : !> \param particles_qm ...
1402 : ! **************************************************************************************************
1403 7920 : SUBROUTINE build_mm_pot(qpot, pot_type, qm_alpha, potentials, &
1404 : particles_mm, mm_charges, mm_atom_index, mm_cell, IndQM, &
1405 : qmmm_spherical_cutoff, particles_qm)
1406 :
1407 : REAL(KIND=dp), INTENT(INOUT) :: qpot
1408 : INTEGER, INTENT(IN) :: pot_type
1409 : REAL(KIND=dp), INTENT(IN) :: qm_alpha
1410 : TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1411 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
1412 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges
1413 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1414 : TYPE(cell_type), POINTER :: mm_cell
1415 : INTEGER, INTENT(IN) :: IndQM
1416 : REAL(KIND=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
1417 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
1418 :
1419 : CHARACTER(len=*), PARAMETER :: routineN = 'build_mm_pot'
1420 : REAL(KIND=dp), PARAMETER :: qsmall = 1.0e-15_dp
1421 :
1422 : INTEGER :: handle, Imm, Imp, IndMM, Ipot
1423 : REAL(KIND=dp) :: dr, qeff, rt1, rt2, rt3, &
1424 : sph_chrg_factor, sr
1425 : REAL(KIND=dp), DIMENSION(3) :: r_pbc, rij
1426 : TYPE(qmmm_pot_type), POINTER :: Pot
1427 :
1428 7920 : CALL timeset(routineN, handle)
1429 : ! Loop Over MM atoms
1430 : ! Loop over Pot stores atoms with the same charge
1431 23568 : MainLoopPot: DO Ipot = 1, SIZE(Potentials)
1432 15648 : Pot => Potentials(Ipot)%Pot
1433 : ! Loop over atoms belonging to this type
1434 70512 : LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index)
1435 46944 : Imm = Pot%mm_atom_index(Imp)
1436 46944 : IndMM = mm_atom_index(Imm)
1437 187776 : r_pbc = pbc(particles_mm(IndMM)%r - particles_qm(IndQM)%r, mm_cell)
1438 46944 : rt1 = r_pbc(1)
1439 46944 : rt2 = r_pbc(2)
1440 46944 : rt3 = r_pbc(3)
1441 187776 : rij = [rt1, rt2, rt3]
1442 187776 : dr = SQRT(SUM(rij**2))
1443 46944 : qeff = mm_charges(Imm)
1444 : ! Computes the screening factor for the spherical cutoff (if defined)
1445 46944 : IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1446 28584 : CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
1447 28584 : qeff = qeff*sph_chrg_factor
1448 : END IF
1449 46944 : IF (ABS(qeff) <= qsmall) CYCLE LoopMM
1450 62592 : IF (dr > rtiny) THEN
1451 46944 : IF (pot_type == pot_tb_nonperiodic) THEN
1452 16344 : sr = gamma_rab_sr(dr, qm_alpha, eta_mm, 0.0_dp)
1453 16344 : qpot = qpot + qeff*(1.0_dp/dr - sr)
1454 : ELSE IF (pot_type == pot_tb_short_range) THEN
1455 12276 : sr = gamma_rab_sr(dr, qm_alpha, eta_mm, 0.0_dp)
1456 12276 : qpot = qpot - qeff*sr
1457 : ELSE IF (pot_type == pot_ewald_short_range) THEN
1458 14292 : sr = erfc(qm_alpha*dr)/dr
1459 14292 : qpot = qpot + qeff*sr
1460 : ELSE IF (pot_type == pot_gauss_nonperiodic) THEN
1461 2016 : sr = erf(dr/Pot%Rc)/dr
1462 2016 : qpot = qpot + qeff*sr
1463 : ELSE IF (pot_type == pot_gauss_short_range) THEN
1464 2016 : sr = erfc(dr/Pot%Rc)/dr
1465 2016 : qpot = qpot - qeff*sr
1466 : ELSE
1467 0 : CPABORT("Unknown pot_type for dr > rtiny")
1468 : END IF
1469 0 : ELSE IF (pot_type == pot_gauss_nonperiodic) THEN
1470 0 : qpot = qpot + qeff*2.0_dp/(rootpi*Pot%Rc)
1471 : END IF
1472 : END DO LoopMM
1473 : END DO MainLoopPot
1474 7920 : CALL timestop(handle)
1475 7920 : END SUBROUTINE build_mm_pot
1476 :
1477 : ! **************************************************************************************************
1478 : !> \brief ...
1479 : !> \param qcharge ...
1480 : !> \param pot_type ...
1481 : !> \param qm_alpha ...
1482 : !> \param potentials ...
1483 : !> \param particles_mm ...
1484 : !> \param mm_charges ...
1485 : !> \param mm_atom_index ...
1486 : !> \param mm_cell ...
1487 : !> \param IndQM ...
1488 : !> \param forces ...
1489 : !> \param forces_qm ...
1490 : !> \param qmmm_spherical_cutoff ...
1491 : !> \param particles_qm ...
1492 : ! **************************************************************************************************
1493 468 : SUBROUTINE build_mm_dpot(qcharge, pot_type, qm_alpha, potentials, &
1494 : particles_mm, mm_charges, mm_atom_index, mm_cell, IndQM, &
1495 234 : forces, forces_qm, qmmm_spherical_cutoff, particles_qm)
1496 :
1497 : REAL(KIND=dp), INTENT(IN) :: qcharge
1498 : INTEGER, INTENT(IN) :: pot_type
1499 : REAL(KIND=dp), INTENT(IN) :: qm_alpha
1500 : TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1501 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
1502 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges
1503 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1504 : TYPE(cell_type), POINTER :: mm_cell
1505 : INTEGER, INTENT(IN) :: IndQM
1506 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: forces
1507 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: forces_qm
1508 : REAL(KIND=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
1509 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
1510 :
1511 : CHARACTER(len=*), PARAMETER :: routineN = 'build_mm_dpot'
1512 : REAL(KIND=dp), PARAMETER :: qsmall = 1.0e-15_dp
1513 :
1514 : INTEGER :: handle, Imm, Imp, IndMM, Ipot
1515 : REAL(KIND=dp) :: dr, drm, drp, dsr, fsr, qeff, rt1, rt2, &
1516 : rt3, sph_chrg_factor
1517 : REAL(KIND=dp), DIMENSION(3) :: force_ab, r_pbc, rij
1518 : TYPE(qmmm_pot_type), POINTER :: Pot
1519 :
1520 234 : CALL timeset(routineN, handle)
1521 : ! Loop Over MM atoms
1522 : ! Loop over Pot stores atoms with the same charge
1523 606 : MainLoopPot: DO Ipot = 1, SIZE(Potentials)
1524 372 : Pot => Potentials(Ipot)%Pot
1525 : ! Loop over atoms belonging to this type
1526 1722 : LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index)
1527 1116 : Imm = Pot%mm_atom_index(Imp)
1528 1116 : IndMM = mm_atom_index(Imm)
1529 4464 : r_pbc = pbc(particles_mm(IndMM)%r - particles_qm(IndQM)%r, mm_cell)
1530 1116 : rt1 = r_pbc(1)
1531 1116 : rt2 = r_pbc(2)
1532 1116 : rt3 = r_pbc(3)
1533 4464 : rij = [rt1, rt2, rt3]
1534 4464 : dr = SQRT(SUM(rij**2))
1535 1116 : qeff = mm_charges(Imm)
1536 : ! Computes the screening factor for the spherical cutoff (if defined)
1537 : ! We neglect derivative of cutoff function for gradients!!!
1538 1116 : IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1539 684 : CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
1540 684 : qeff = qeff*sph_chrg_factor
1541 : END IF
1542 1116 : IF (ABS(qeff) <= qsmall) CYCLE LoopMM
1543 1116 : IF (dr > rtiny) THEN
1544 1116 : drp = dr + ddrmm
1545 1116 : drm = dr - ddrmm
1546 1116 : IF (pot_type == pot_tb_nonperiodic) THEN
1547 : dsr = 0.5_dp*(gamma_rab_sr(drp, qm_alpha, eta_mm, 0.0_dp) - &
1548 396 : gamma_rab_sr(drm, qm_alpha, eta_mm, 0.0_dp))/ddrmm
1549 396 : fsr = qeff*qcharge*(-1.0_dp/(dr*dr) - dsr)
1550 : ELSE IF (pot_type == pot_tb_short_range) THEN
1551 : dsr = 0.5_dp*(gamma_rab_sr(drp, qm_alpha, eta_mm, 0.0_dp) - &
1552 306 : gamma_rab_sr(drm, qm_alpha, eta_mm, 0.0_dp))/ddrmm
1553 306 : fsr = -qeff*qcharge*dsr
1554 : ELSE IF (pot_type == pot_ewald_short_range) THEN
1555 342 : dsr = 0.5_dp*(erfc(qm_alpha*drp)/drp - erfc(qm_alpha*drm)/drm)/ddrmm
1556 342 : fsr = qeff*qcharge*dsr
1557 : ELSE IF (pot_type == pot_gauss_nonperiodic) THEN
1558 : fsr = qeff*qcharge*(2.0_dp*EXP(-(dr/Pot%Rc)**2)/(rootpi*Pot%Rc*dr) - &
1559 36 : erf(dr/Pot%Rc)/dr**2)
1560 : ELSE IF (pot_type == pot_gauss_short_range) THEN
1561 : fsr = qeff*qcharge*(2.0_dp*EXP(-(dr/Pot%Rc)**2)/(rootpi*Pot%Rc*dr) + &
1562 36 : erfc(dr/Pot%Rc)/dr**2)
1563 : ELSE
1564 0 : CPABORT("Unknown pot_type for dr > rtiny")
1565 : END IF
1566 4464 : force_ab = -fsr*rij/dr
1567 : ELSE
1568 0 : force_ab = 0.0_dp
1569 : END IF
1570 : ! The array of QM forces are really the forces
1571 4464 : forces_qm(:) = forces_qm(:) - force_ab
1572 : ! The one of MM atoms are instead gradients
1573 4836 : forces(:, Imm) = forces(:, Imm) - force_ab
1574 : END DO LoopMM
1575 : END DO MainLoopPot
1576 :
1577 234 : CALL timestop(handle)
1578 :
1579 234 : END SUBROUTINE build_mm_dpot
1580 :
1581 : END MODULE qmmm_tb_methods
|
