(git:ed6f26b)
Loading...
Searching...
No Matches
qmmm_tb_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief TB methods used with QMMM
10!> \author JGH
11! **************************************************************************************************
12MODULE 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 message_passing, ONLY: mp_para_env_type
40 USE mulliken, ONLY: mulliken_charges
41 USE particle_types, ONLY: allocate_particle_set,&
42 deallocate_particle_set,&
43 particle_type
44 USE pw_poisson_types, ONLY: do_ewald_ewald,&
45 do_ewald_none,&
46 do_ewald_pme,&
47 do_ewald_spme
48 USE qmmm_types_low, ONLY: qmmm_env_qm_type,&
49 qmmm_pot_p_type,&
50 qmmm_pot_type
51 USE qmmm_util, ONLY: spherical_cutoff_factor
52 USE qs_dftb_coulomb, ONLY: gamma_rab_sr
53 USE qs_dftb_matrices, ONLY: build_dftb_overlap
54 USE qs_dftb_types, ONLY: qs_dftb_atom_type
55 USE qs_dftb_utils, ONLY: get_dftb_atom_param
56 USE qs_environment_types, ONLY: get_qs_env,&
57 qs_environment_type
58 USE qs_kind_types, ONLY: get_qs_kind,&
59 qs_kind_type
60 USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
61 USE qs_ks_types, ONLY: qs_ks_env_type
62 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
63 USE qs_neighbor_lists, ONLY: build_qs_neighbor_lists
64 USE qs_overlap, ONLY: build_overlap_matrix
65 USE qs_rho_types, ONLY: qs_rho_get,&
66 qs_rho_type
67 USE spme, ONLY: spme_forces,&
68 spme_potential
69 USE xtb_types, ONLY: get_xtb_atom_param,&
70 xtb_atom_type
71#include "./base/base_uses.f90"
72
73 IMPLICIT NONE
74
75 ! small real number
76 REAL(dp), PARAMETER :: rtiny = 1.e-10_dp
77 ! eta(0) for mm atoms and non-scc qm atoms
78 REAL(dp), PARAMETER :: eta_mm = 0.47_dp
79 ! step size for qmmm finite difference
80 REAL(dp), PARAMETER :: ddrmm = 0.0001_dp
81
82 PRIVATE
83
84 CHARACTER(len=*), PARAMETER, PRIVATE :: modulen = 'qmmm_tb_methods'
85
86 PUBLIC :: build_tb_qmmm_matrix, build_tb_qmmm_matrix_zero, &
87 build_tb_qmmm_matrix_pc, deriv_tb_qmmm_matrix, &
88 deriv_tb_qmmm_matrix_pc
89
90CONTAINS
91
92! **************************************************************************************************
93!> \brief Constructs the 1-el DFTB hamiltonian
94!> \param qs_env ...
95!> \param qmmm_env ...
96!> \param particles_mm ...
97!> \param mm_cell ...
98!> \param para_env ...
99!> \author JGH 10.2014 [created]
100! **************************************************************************************************
101 SUBROUTINE build_tb_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
102
103 TYPE(qs_environment_type), POINTER :: qs_env
104 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
105 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
106 TYPE(cell_type), POINTER :: mm_cell
107 TYPE(mp_para_env_type), POINTER :: para_env
108
109 CHARACTER(len=*), PARAMETER :: routinen = 'build_tb_qmmm_matrix'
110
111 INTEGER :: handle, i, iatom, ikind, jatom, natom, &
112 natorb, nkind
113 INTEGER, DIMENSION(:), POINTER :: list
114 LOGICAL :: defined, do_dftb, do_xtb, found
115 REAL(kind=dp) :: pc_ener, zeff
116 REAL(kind=dp), DIMENSION(0:3) :: eta_a
117 REAL(kind=dp), DIMENSION(:), POINTER :: qpot
118 REAL(kind=dp), DIMENSION(:, :), POINTER :: hblock, sblock
119 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
120 TYPE(dbcsr_iterator_type) :: iter
121 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_s
122 TYPE(dft_control_type), POINTER :: dft_control
123 TYPE(dftb_control_type), POINTER :: dftb_control
124 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
125 POINTER :: sab_nl
126 TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
127 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
128 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
129 TYPE(qs_ks_env_type), POINTER :: ks_env
130 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
131 TYPE(qs_rho_type), POINTER :: rho
132 TYPE(xtb_atom_type), POINTER :: xtb_kind
133 TYPE(xtb_control_type), POINTER :: xtb_control
134
135 CALL timeset(routinen, handle)
136
137 CALL get_qs_env(qs_env=qs_env, &
138 dft_control=dft_control, &
139 atomic_kind_set=atomic_kind_set, &
140 particle_set=particles_qm, &
141 qs_kind_set=qs_kind_set, &
142 rho=rho, &
143 natom=natom)
144 dftb_control => dft_control%qs_control%dftb_control
145 xtb_control => dft_control%qs_control%xtb_control
146
147 IF (dft_control%qs_control%dftb) THEN
148 do_dftb = .true.
149 do_xtb = .false.
150 ELSEIF (dft_control%qs_control%xtb) THEN
151 do_dftb = .false.
152 do_xtb = .true.
153 ELSE
154 cpabort("TB method unknown")
155 END IF
156
157 CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
158
159 NULLIFY (matrix_s)
160 IF (do_dftb) THEN
161 CALL build_dftb_overlap(qs_env, 0, matrix_s)
162 ELSEIF (do_xtb) THEN
163 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
164 CALL build_overlap_matrix(ks_env, matrix_s, basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
165 END IF
166
167 ALLOCATE (qpot(natom))
168 qpot = 0.0_dp
169 pc_ener = 0.0_dp
170
171 nkind = SIZE(atomic_kind_set)
172 DO ikind = 1, nkind
173 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
174 IF (do_dftb) THEN
175 NULLIFY (dftb_kind)
176 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
177 CALL get_dftb_atom_param(dftb_kind, zeff=zeff, &
178 defined=defined, eta=eta_a, natorb=natorb)
179 ! use mm charge smearing for non-scc cases
180 IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
181 IF (.NOT. defined .OR. natorb < 1) cycle
182 ELSEIF (do_xtb) THEN
183 NULLIFY (xtb_kind)
184 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
185 CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
186 eta_a(0) = eta_mm
187 END IF
188 DO i = 1, SIZE(list)
189 iatom = list(i)
190 CALL build_mm_pot(qpot(iatom), 0, eta_a(0), qmmm_env%Potentials, particles_mm, &
191 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, &
192 qmmm_env%spherical_cutoff, particles_qm)
193 ! Possibly added charges
194 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
195 CALL build_mm_pot(qpot(iatom), 0, eta_a(0), qmmm_env%added_charges%potentials, &
196 qmmm_env%added_charges%added_particles, &
197 qmmm_env%added_charges%mm_atom_chrg, &
198 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
199 qmmm_env%spherical_cutoff, &
200 particles_qm)
201 END IF
202 pc_ener = pc_ener + qpot(iatom)*zeff
203 END DO
204 END DO
205
206 ! Allocate the core Hamiltonian matrix
207 CALL get_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env_loc)
208 matrix_h => ks_qmmm_env_loc%matrix_h
209 CALL dbcsr_allocate_matrix_set(matrix_h, 1)
210 ALLOCATE (matrix_h(1)%matrix)
211 CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, &
212 name="QMMM HAMILTONIAN MATRIX")
213 CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
214
215 CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
216 DO WHILE (dbcsr_iterator_blocks_left(iter))
217 CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock)
218 NULLIFY (hblock)
219 CALL dbcsr_get_block_p(matrix=matrix_h(1)%matrix, &
220 row=iatom, col=jatom, block=hblock, found=found)
221 cpassert(found)
222 hblock = hblock - 0.5_dp*sblock*(qpot(iatom) + qpot(jatom))
223 END DO
224 CALL dbcsr_iterator_stop(iter)
225
226 ks_qmmm_env_loc%matrix_h => matrix_h
227 ks_qmmm_env_loc%pc_ener = pc_ener
228
229 DEALLOCATE (qpot)
230
231 CALL dbcsr_deallocate_matrix_set(matrix_s)
232
233 CALL timestop(handle)
234
235 END SUBROUTINE build_tb_qmmm_matrix
236
237! **************************************************************************************************
238!> \brief Constructs an empty 1-el DFTB hamiltonian
239!> \param qs_env ...
240!> \param para_env ...
241!> \author JGH 10.2014 [created]
242! **************************************************************************************************
243 SUBROUTINE build_tb_qmmm_matrix_zero(qs_env, para_env)
244
245 TYPE(qs_environment_type), POINTER :: qs_env
246 TYPE(mp_para_env_type), POINTER :: para_env
247
248 CHARACTER(len=*), PARAMETER :: routinen = 'build_tb_qmmm_matrix_zero'
249
250 INTEGER :: handle
251 LOGICAL :: do_dftb, do_xtb
252 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_s
253 TYPE(dft_control_type), POINTER :: dft_control
254 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
255 POINTER :: sab_nl
256 TYPE(qs_ks_env_type), POINTER :: ks_env
257 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
258
259 CALL timeset(routinen, handle)
260
261 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
262
263 IF (dft_control%qs_control%dftb) THEN
264 do_dftb = .true.
265 do_xtb = .false.
266 ELSEIF (dft_control%qs_control%xtb) THEN
267 do_dftb = .false.
268 do_xtb = .true.
269 ELSE
270 cpabort("TB method unknown")
271 END IF
272
273 CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
274
275 NULLIFY (matrix_s)
276 IF (do_dftb) THEN
277 CALL build_dftb_overlap(qs_env, 0, matrix_s)
278 ELSEIF (do_xtb) THEN
279 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
280 CALL build_overlap_matrix(ks_env, matrix_s, basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
281 END IF
282
283 ! Allocate the core Hamiltonian matrix
284 CALL get_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env_loc)
285 matrix_h => ks_qmmm_env_loc%matrix_h
286 CALL dbcsr_allocate_matrix_set(matrix_h, 1)
287 ALLOCATE (matrix_h(1)%matrix)
288 CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, &
289 name="QMMM HAMILTONIAN MATRIX")
290 CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
291 ks_qmmm_env_loc%matrix_h => matrix_h
292 ks_qmmm_env_loc%pc_ener = 0.0_dp
293
294 CALL dbcsr_deallocate_matrix_set(matrix_s)
295
296 CALL timestop(handle)
297
298 END SUBROUTINE build_tb_qmmm_matrix_zero
299
300! **************************************************************************************************
301!> \brief Constructs the 1-el DFTB hamiltonian
302!> \param qs_env ...
303!> \param qmmm_env ...
304!> \param particles_mm ...
305!> \param mm_cell ...
306!> \param para_env ...
307!> \author JGH 10.2014 [created]
308! **************************************************************************************************
309 SUBROUTINE build_tb_qmmm_matrix_pc(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
310
311 TYPE(qs_environment_type), POINTER :: qs_env
312 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
313 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
314 TYPE(cell_type), POINTER :: mm_cell
315 TYPE(mp_para_env_type), POINTER :: para_env
316
317 CHARACTER(len=*), PARAMETER :: routinen = 'build_tb_qmmm_matrix_pc'
318
319 INTEGER :: do_ipol, ewald_type, handle, i, iatom, &
320 ikind, imm, imp, indmm, ipot, jatom, &
321 natom, natorb, nkind, nmm
322 INTEGER, DIMENSION(:), POINTER :: list
323 LOGICAL :: defined, do_dftb, do_multipoles, do_xtb, &
324 found
325 REAL(kind=dp) :: alpha, pc_ener, zeff
326 REAL(kind=dp), DIMENSION(0:3) :: eta_a
327 REAL(kind=dp), DIMENSION(2) :: rcutoff
328 REAL(kind=dp), DIMENSION(:), POINTER :: charges_mm, qpot
329 REAL(kind=dp), DIMENSION(:, :), POINTER :: hblock, sblock
330 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
331 TYPE(dbcsr_iterator_type) :: iter
332 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_s
333 TYPE(dft_control_type), POINTER :: dft_control
334 TYPE(dftb_control_type), POINTER :: dftb_control
335 TYPE(ewald_environment_type), POINTER :: ewald_env
336 TYPE(ewald_pw_type), POINTER :: ewald_pw
337 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
338 POINTER :: sab_nl
339 TYPE(particle_type), DIMENSION(:), POINTER :: atoms_mm, particles_qm
340 TYPE(qmmm_pot_type), POINTER :: pot
341 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
342 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
343 TYPE(qs_ks_env_type), POINTER :: ks_env
344 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
345 TYPE(qs_rho_type), POINTER :: rho
346 TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
347 print_section
348 TYPE(xtb_atom_type), POINTER :: xtb_kind
349 TYPE(xtb_control_type), POINTER :: xtb_control
350
351 CALL timeset(routinen, handle)
352
353 CALL get_qs_env(qs_env=qs_env, &
354 dft_control=dft_control, &
355 atomic_kind_set=atomic_kind_set, &
356 particle_set=particles_qm, &
357 qs_kind_set=qs_kind_set, &
358 rho=rho, &
359 natom=natom)
360 dftb_control => dft_control%qs_control%dftb_control
361 xtb_control => dft_control%qs_control%xtb_control
362
363 IF (dft_control%qs_control%dftb) THEN
364 do_dftb = .true.
365 do_xtb = .false.
366 ELSEIF (dft_control%qs_control%xtb) THEN
367 do_dftb = .false.
368 do_xtb = .true.
369 ELSE
370 cpabort("TB method unknown")
371 END IF
372
373 CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
374
375 NULLIFY (matrix_s)
376 IF (do_dftb) THEN
377 CALL build_dftb_overlap(qs_env, 0, matrix_s)
378 ELSEIF (do_xtb) THEN
379 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
380 CALL build_overlap_matrix(ks_env, matrix_s, basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
381 END IF
382
383 ALLOCATE (qpot(natom))
384 qpot = 0.0_dp
385 pc_ener = 0.0_dp
386
387 ! Create Ewald environments
388 poisson_section => section_vals_get_subs_vals(qs_env%input, "MM%POISSON")
389 ALLOCATE (ewald_env)
390 CALL ewald_env_create(ewald_env, para_env)
391 CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
392 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
393 CALL read_ewald_section(ewald_env, ewald_section)
394 print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
395 ALLOCATE (ewald_pw)
396 CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=print_section)
397
398 CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, do_ipol=do_ipol)
399 IF (do_multipoles) cpabort("No multipole force fields allowed in TB QM/MM")
400 IF (do_ipol /= do_fist_pol_none) cpabort("No polarizable force fields allowed in TB QM/MM")
401
402 SELECT CASE (ewald_type)
403 CASE (do_ewald_pme)
404 cpabort("PME Ewald type not implemented for TB/QMMM")
405 CASE (do_ewald_ewald, do_ewald_spme)
406 DO ipot = 1, SIZE(qmmm_env%Potentials)
407 pot => qmmm_env%Potentials(ipot)%Pot
408 nmm = SIZE(pot%mm_atom_index)
409 ! get a 'clean' mm particle set
410 NULLIFY (atoms_mm)
411 CALL allocate_particle_set(atoms_mm, nmm)
412 ALLOCATE (charges_mm(nmm))
413 DO imp = 1, nmm
414 imm = pot%mm_atom_index(imp)
415 indmm = qmmm_env%mm_atom_index(imm)
416 atoms_mm(imp)%r = particles_mm(indmm)%r
417 atoms_mm(imp)%atomic_kind => particles_mm(indmm)%atomic_kind
418 charges_mm(imp) = qmmm_env%mm_atom_chrg(imm)
419 END DO
420 IF (ewald_type == do_ewald_ewald) THEN
421 cpabort("Ewald not implemented for TB/QMMM")
422 ELSE IF (ewald_type == do_ewald_spme) THEN
423 ! spme electrostatic potential
424 CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, particles_qm, qpot)
425 END IF
426 CALL deallocate_particle_set(atoms_mm)
427 DEALLOCATE (charges_mm)
428 END DO
429 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
430 DO ipot = 1, SIZE(qmmm_env%added_charges%Potentials)
431 pot => qmmm_env%added_charges%Potentials(ipot)%Pot
432 nmm = SIZE(pot%mm_atom_index)
433 ! get a 'clean' mm particle set
434 NULLIFY (atoms_mm)
435 CALL allocate_particle_set(atoms_mm, nmm)
436 ALLOCATE (charges_mm(nmm))
437 DO imp = 1, nmm
438 imm = pot%mm_atom_index(imp)
439 indmm = qmmm_env%added_charges%mm_atom_index(imm)
440 atoms_mm(imp)%r = qmmm_env%added_charges%added_particles(indmm)%r
441 atoms_mm(imp)%atomic_kind => qmmm_env%added_charges%added_particles(indmm)%atomic_kind
442 charges_mm(imp) = qmmm_env%added_charges%mm_atom_chrg(imm)
443 END DO
444 IF (ewald_type == do_ewald_ewald) THEN
445 cpabort("Ewald not implemented for TB/QMMM")
446 ELSE IF (ewald_type == do_ewald_spme) THEN
447 ! spme electrostatic potential
448 CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, particles_qm, qpot)
449 END IF
450 CALL deallocate_particle_set(atoms_mm)
451 DEALLOCATE (charges_mm)
452 END DO
453 END IF
454 CALL para_env%sum(qpot)
455 ! Add Ewald and TB short range corrections
456 ! This is effectively using a minimum image convention!
457 ! Set rcutoff to values compatible with alpha Ewald
458 CALL ewald_env_get(ewald_env, rcut=rcutoff(1), alpha=alpha)
459 rcutoff(2) = 0.025_dp*rcutoff(1)
460 rcutoff(1) = 2.0_dp*rcutoff(1)
461 nkind = SIZE(atomic_kind_set)
462 DO ikind = 1, nkind
463 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
464 IF (do_dftb) THEN
465 NULLIFY (dftb_kind)
466 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
467 CALL get_dftb_atom_param(dftb_kind, zeff=zeff, &
468 defined=defined, eta=eta_a, natorb=natorb)
469 ! use mm charge smearing for non-scc cases
470 IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
471 IF (.NOT. defined .OR. natorb < 1) cycle
472 ELSEIF (do_xtb) THEN
473 NULLIFY (xtb_kind)
474 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
475 CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
476 eta_a(0) = eta_mm
477 END IF
478 DO i = 1, SIZE(list)
479 iatom = list(i)
480 CALL build_mm_pot(qpot(iatom), 1, eta_a(0), qmmm_env%Potentials, particles_mm, &
481 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, rcutoff, &
482 particles_qm)
483 CALL build_mm_pot(qpot(iatom), 2, alpha, qmmm_env%Potentials, particles_mm, &
484 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, rcutoff, &
485 particles_qm)
486 ! Possibly added charges
487 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
488 CALL build_mm_pot(qpot(iatom), 1, eta_a(0), qmmm_env%added_charges%potentials, &
489 qmmm_env%added_charges%added_particles, &
490 qmmm_env%added_charges%mm_atom_chrg, &
491 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, rcutoff, &
492 particles_qm)
493 CALL build_mm_pot(qpot(iatom), 2, alpha, qmmm_env%added_charges%potentials, &
494 qmmm_env%added_charges%added_particles, &
495 qmmm_env%added_charges%mm_atom_chrg, &
496 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, rcutoff, &
497 particles_qm)
498 END IF
499 pc_ener = pc_ener + qpot(iatom)*zeff
500 END DO
501 END DO
502 CASE (do_ewald_none)
503 ! Simply summing up charges with 1/R (DFTB corrected)
504 nkind = SIZE(atomic_kind_set)
505 DO ikind = 1, nkind
506 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
507 IF (do_dftb) THEN
508 NULLIFY (dftb_kind)
509 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
510 CALL get_dftb_atom_param(dftb_kind, zeff=zeff, &
511 defined=defined, eta=eta_a, natorb=natorb)
512 ! use mm charge smearing for non-scc cases
513 IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
514 IF (.NOT. defined .OR. natorb < 1) cycle
515 ELSEIF (do_xtb) THEN
516 NULLIFY (xtb_kind)
517 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
518 CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
519 eta_a(0) = eta_mm
520 END IF
521 DO i = 1, SIZE(list)
522 iatom = list(i)
523 CALL build_mm_pot(qpot(iatom), 0, eta_a(0), qmmm_env%Potentials, particles_mm, &
524 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, &
525 qmmm_env%spherical_cutoff, particles_qm)
526 ! Possibly added charges
527 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
528 CALL build_mm_pot(qpot(iatom), 0, eta_a(0), qmmm_env%added_charges%potentials, &
529 qmmm_env%added_charges%added_particles, &
530 qmmm_env%added_charges%mm_atom_chrg, &
531 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
532 qmmm_env%spherical_cutoff, &
533 particles_qm)
534 END IF
535 pc_ener = pc_ener + qpot(iatom)*zeff
536 END DO
537 END DO
538 CASE DEFAULT
539 cpabort("Unknown Ewald type!")
540 END SELECT
541
542 ! Allocate the core Hamiltonian matrix
543 CALL get_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env_loc)
544 matrix_h => ks_qmmm_env_loc%matrix_h
545 CALL dbcsr_allocate_matrix_set(matrix_h, 1)
546 ALLOCATE (matrix_h(1)%matrix)
547 CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, &
548 name="QMMM HAMILTONIAN MATRIX")
549 CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
550
551 CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
552 DO WHILE (dbcsr_iterator_blocks_left(iter))
553 CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock)
554 NULLIFY (hblock)
555 CALL dbcsr_get_block_p(matrix=matrix_h(1)%matrix, &
556 row=iatom, col=jatom, block=hblock, found=found)
557 cpassert(found)
558 hblock = hblock - 0.5_dp*sblock*(qpot(iatom) + qpot(jatom))
559 END DO
560 CALL dbcsr_iterator_stop(iter)
561
562 ks_qmmm_env_loc%matrix_h => matrix_h
563 ks_qmmm_env_loc%pc_ener = pc_ener
564
565 DEALLOCATE (qpot)
566
567 ! Release Ewald environment
568 CALL ewald_env_release(ewald_env)
569 DEALLOCATE (ewald_env)
570 CALL ewald_pw_release(ewald_pw)
571 DEALLOCATE (ewald_pw)
572
573 CALL dbcsr_deallocate_matrix_set(matrix_s)
574
575 CALL timestop(handle)
576
577 END SUBROUTINE build_tb_qmmm_matrix_pc
578
579! **************************************************************************************************
580!> \brief Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms
581!> \param qs_env ...
582!> \param qmmm_env ...
583!> \param particles_mm ...
584!> \param mm_cell ...
585!> \param para_env ...
586!> \param calc_force ...
587!> \param Forces ...
588!> \param Forces_added_charges ...
589!> \author JGH 10.2014 [created]
590! **************************************************************************************************
591 SUBROUTINE deriv_tb_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
592 calc_force, Forces, Forces_added_charges)
593
594 TYPE(qs_environment_type), POINTER :: qs_env
595 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
596 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
597 TYPE(cell_type), POINTER :: mm_cell
598 TYPE(mp_para_env_type), POINTER :: para_env
599 LOGICAL, INTENT(in), OPTIONAL :: calc_force
600 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces, forces_added_charges
601
602 CHARACTER(len=*), PARAMETER :: routinen = 'deriv_tb_qmmm_matrix'
603
604 INTEGER :: atom_a, handle, i, iatom, ikind, iqm, &
605 jatom, natom, natorb, nkind, nspins, &
606 number_qm_atoms
607 INTEGER, DIMENSION(:), POINTER :: list
608 LOGICAL :: defined, do_dftb, do_xtb, found
609 REAL(kind=dp) :: fi, gmij, zeff
610 REAL(kind=dp), DIMENSION(0:3) :: eta_a
611 REAL(kind=dp), DIMENSION(:), POINTER :: mcharge, qpot
612 REAL(kind=dp), DIMENSION(:, :), POINTER :: charges, dsblock, forces_qm, pblock, &
613 sblock
614 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
615 TYPE(dbcsr_iterator_type) :: iter
616 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s
617 TYPE(dft_control_type), POINTER :: dft_control
618 TYPE(dftb_control_type), POINTER :: dftb_control
619 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
620 POINTER :: sab_nl
621 TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
622 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
623 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
624 TYPE(qs_ks_env_type), POINTER :: ks_env
625 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
626 TYPE(qs_rho_type), POINTER :: rho
627 TYPE(xtb_atom_type), POINTER :: xtb_kind
628 TYPE(xtb_control_type), POINTER :: xtb_control
629
630 CALL timeset(routinen, handle)
631 IF (calc_force) THEN
632 NULLIFY (rho, atomic_kind_set, qs_kind_set, particles_qm)
633 CALL get_qs_env(qs_env=qs_env, &
634 rho=rho, &
635 atomic_kind_set=atomic_kind_set, &
636 qs_kind_set=qs_kind_set, &
637 ks_qmmm_env=ks_qmmm_env_loc, &
638 dft_control=dft_control, &
639 particle_set=particles_qm, &
640 natom=number_qm_atoms)
641 dftb_control => dft_control%qs_control%dftb_control
642 xtb_control => dft_control%qs_control%xtb_control
643
644 IF (dft_control%qs_control%dftb) THEN
645 do_dftb = .true.
646 do_xtb = .false.
647 ELSEIF (dft_control%qs_control%xtb) THEN
648 do_dftb = .false.
649 do_xtb = .true.
650 ELSE
651 cpabort("TB method unknown")
652 END IF
653
654 NULLIFY (matrix_s)
655 IF (do_dftb) THEN
656 CALL build_dftb_overlap(qs_env, 1, matrix_s)
657 ELSEIF (do_xtb) THEN
658 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
659 CALL build_overlap_matrix(ks_env, matrix_s, nderivative=1, &
660 basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
661 END IF
662
663 CALL qs_rho_get(rho, rho_ao=matrix_p)
664
665 nspins = dft_control%nspins
666 nkind = SIZE(atomic_kind_set)
667 ! Mulliken charges
668 ALLOCATE (charges(number_qm_atoms, nspins))
669 !
670 CALL mulliken_charges(matrix_p, matrix_s(1)%matrix, para_env, charges)
671 !
672 ALLOCATE (mcharge(number_qm_atoms))
673 DO ikind = 1, nkind
674 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
675 IF (do_dftb) THEN
676 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
677 CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
678 ELSEIF (do_xtb) THEN
679 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
680 CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
681 END IF
682 DO iatom = 1, natom
683 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
684 mcharge(atom_a) = zeff - sum(charges(atom_a, 1:nspins))
685 END DO
686 END DO
687 DEALLOCATE (charges)
688
689 ALLOCATE (qpot(number_qm_atoms))
690 qpot = 0.0_dp
691 ALLOCATE (forces_qm(3, number_qm_atoms))
692 forces_qm = 0.0_dp
693
694 ! calculate potential and forces from classical charges
695 iqm = 0
696 DO ikind = 1, nkind
697 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
698 IF (do_dftb) THEN
699 NULLIFY (dftb_kind)
700 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
701 CALL get_dftb_atom_param(dftb_kind, &
702 defined=defined, eta=eta_a, natorb=natorb)
703 ! use mm charge smearing for non-scc cases
704 IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
705 IF (.NOT. defined .OR. natorb < 1) cycle
706 ELSEIF (do_xtb) THEN
707 eta_a(0) = eta_mm
708 END IF
709 DO i = 1, SIZE(list)
710 iatom = list(i)
711 iqm = iqm + 1
712 CALL build_mm_pot(qpot(iatom), 0, eta_a(0), qmmm_env%Potentials, particles_mm, &
713 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, &
714 qmmm_env%spherical_cutoff, particles_qm)
715 CALL build_mm_dpot(mcharge(iatom), 0, eta_a(0), qmmm_env%Potentials, particles_mm, &
716 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
717 mm_cell, iatom, forces, forces_qm(:, iqm), &
718 qmmm_env%spherical_cutoff, particles_qm)
719 ! Possibly added charges
720 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
721 CALL build_mm_pot(qpot(iatom), 0, eta_a(0), qmmm_env%added_charges%potentials, &
722 qmmm_env%added_charges%added_particles, &
723 qmmm_env%added_charges%mm_atom_chrg, &
724 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
725 qmmm_env%spherical_cutoff, &
726 particles_qm)
727 CALL build_mm_dpot(mcharge(iatom), 0, eta_a(0), qmmm_env%added_charges%potentials, &
728 qmmm_env%added_charges%added_particles, &
729 qmmm_env%added_charges%mm_atom_chrg, &
730 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
731 forces_added_charges, &
732 forces_qm(:, iqm), qmmm_env%spherical_cutoff, particles_qm)
733 END IF
734 END DO
735 END DO
736
737 ! Transfer QM gradients to the QM particles..
738 iqm = 0
739 DO ikind = 1, nkind
740 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
741 IF (do_dftb) THEN
742 NULLIFY (dftb_kind)
743 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
744 CALL get_dftb_atom_param(dftb_kind, defined=defined, natorb=natorb)
745 IF (.NOT. defined .OR. natorb < 1) cycle
746 ELSEIF (do_xtb) THEN
747 ! use all kinds
748 END IF
749 DO i = 1, SIZE(list)
750 iqm = iqm + 1
751 iatom = qmmm_env%qm_atom_index(list(i))
752 particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + forces_qm(:, iqm)
753 END DO
754 END DO
755
756 ! derivatives from qm charges
757 forces_qm = 0.0_dp
758 IF (SIZE(matrix_p) == 2) THEN
759 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
760 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
761 END IF
762 !
763 CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
764 DO WHILE (dbcsr_iterator_blocks_left(iter))
765 CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock)
766 !
767 IF (iatom == jatom) cycle
768 !
769 gmij = -0.5_dp*(qpot(iatom) + qpot(jatom))
770 NULLIFY (pblock)
771 CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
772 row=iatom, col=jatom, block=pblock, found=found)
773 cpassert(found)
774 DO i = 1, 3
775 NULLIFY (dsblock)
776 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
777 row=iatom, col=jatom, block=dsblock, found=found)
778 cpassert(found)
779 fi = -2.0_dp*gmij*sum(pblock*dsblock)
780 forces_qm(i, iatom) = forces_qm(i, iatom) + fi
781 forces_qm(i, jatom) = forces_qm(i, jatom) - fi
782 END DO
783 END DO
784 CALL dbcsr_iterator_stop(iter)
785 !
786 IF (SIZE(matrix_p) == 2) THEN
787 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
788 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
789 END IF
790 !
791 ! Transfer QM gradients to the QM particles..
792 CALL para_env%sum(forces_qm)
793 DO ikind = 1, nkind
794 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
795 DO i = 1, SIZE(list)
796 iqm = list(i)
797 iatom = qmmm_env%qm_atom_index(iqm)
798 particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + forces_qm(:, iqm)
799 END DO
800 END DO
801 !
802 DEALLOCATE (mcharge)
803 !
804 ! MM forces will be handled directly from the QMMM module in the same way
805 ! as for GPW/GAPW methods
806 DEALLOCATE (forces_qm)
807 DEALLOCATE (qpot)
808
809 CALL dbcsr_deallocate_matrix_set(matrix_s)
810
811 END IF
812 CALL timestop(handle)
813
814 END SUBROUTINE deriv_tb_qmmm_matrix
815
816! **************************************************************************************************
817!> \brief Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms
818!> \param qs_env ...
819!> \param qmmm_env ...
820!> \param particles_mm ...
821!> \param mm_cell ...
822!> \param para_env ...
823!> \param calc_force ...
824!> \param Forces ...
825!> \param Forces_added_charges ...
826!> \author JGH 10.2014 [created]
827! **************************************************************************************************
828 SUBROUTINE deriv_tb_qmmm_matrix_pc(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
829 calc_force, Forces, Forces_added_charges)
830
831 TYPE(qs_environment_type), POINTER :: qs_env
832 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
833 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
834 TYPE(cell_type), POINTER :: mm_cell
835 TYPE(mp_para_env_type), POINTER :: para_env
836 LOGICAL, INTENT(in), OPTIONAL :: calc_force
837 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces, forces_added_charges
838
839 CHARACTER(len=*), PARAMETER :: routinen = 'deriv_tb_qmmm_matrix_pc'
840
841 INTEGER :: atom_a, do_ipol, ewald_type, handle, i, iatom, ikind, imm, imp, indmm, ipot, iqm, &
842 jatom, natom, natorb, nkind, nmm, nspins, number_qm_atoms
843 INTEGER, DIMENSION(:), POINTER :: list
844 LOGICAL :: defined, do_dftb, do_multipoles, do_xtb, &
845 found
846 REAL(kind=dp) :: alpha, fi, gmij, zeff
847 REAL(kind=dp), DIMENSION(0:3) :: eta_a
848 REAL(kind=dp), DIMENSION(2) :: rcutoff
849 REAL(kind=dp), DIMENSION(:), POINTER :: charges_mm, mcharge, qpot
850 REAL(kind=dp), DIMENSION(:, :), POINTER :: charges, dsblock, forces_mm, forces_qm, &
851 pblock, sblock
852 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
853 TYPE(dbcsr_iterator_type) :: iter
854 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s
855 TYPE(dft_control_type), POINTER :: dft_control
856 TYPE(dftb_control_type), POINTER :: dftb_control
857 TYPE(ewald_environment_type), POINTER :: ewald_env
858 TYPE(ewald_pw_type), POINTER :: ewald_pw
859 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
860 POINTER :: sab_nl
861 TYPE(particle_type), DIMENSION(:), POINTER :: atoms_mm, particles_qm
862 TYPE(qmmm_pot_type), POINTER :: pot
863 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
864 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
865 TYPE(qs_ks_env_type), POINTER :: ks_env
866 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
867 TYPE(qs_rho_type), POINTER :: rho
868 TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
869 print_section
870 TYPE(xtb_atom_type), POINTER :: xtb_kind
871 TYPE(xtb_control_type), POINTER :: xtb_control
872
873 CALL timeset(routinen, handle)
874 IF (calc_force) THEN
875 NULLIFY (rho, atomic_kind_set, qs_kind_set, particles_qm)
876 CALL get_qs_env(qs_env=qs_env, &
877 rho=rho, &
878 atomic_kind_set=atomic_kind_set, &
879 qs_kind_set=qs_kind_set, &
880 ks_qmmm_env=ks_qmmm_env_loc, &
881 dft_control=dft_control, &
882 particle_set=particles_qm, &
883 natom=number_qm_atoms)
884 dftb_control => dft_control%qs_control%dftb_control
885 xtb_control => dft_control%qs_control%xtb_control
886
887 IF (dft_control%qs_control%dftb) THEN
888 do_dftb = .true.
889 do_xtb = .false.
890 ELSEIF (dft_control%qs_control%xtb) THEN
891 do_dftb = .false.
892 do_xtb = .true.
893 ELSE
894 cpabort("TB method unknown")
895 END IF
896
897 NULLIFY (matrix_s)
898 IF (do_dftb) THEN
899 CALL build_dftb_overlap(qs_env, 1, matrix_s)
900 ELSEIF (do_xtb) THEN
901 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
902 CALL build_overlap_matrix(ks_env, matrix_s, nderivative=1, &
903 basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
904 END IF
905 CALL qs_rho_get(rho, rho_ao=matrix_p)
906
907 nspins = dft_control%nspins
908 nkind = SIZE(atomic_kind_set)
909 ! Mulliken charges
910 ALLOCATE (charges(number_qm_atoms, nspins))
911 !
912 CALL mulliken_charges(matrix_p, matrix_s(1)%matrix, para_env, charges)
913 !
914 ALLOCATE (mcharge(number_qm_atoms))
915 DO ikind = 1, nkind
916 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
917 IF (do_dftb) THEN
918 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
919 CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
920 ELSEIF (do_xtb) THEN
921 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
922 CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
923 END IF
924 DO iatom = 1, natom
925 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
926 mcharge(atom_a) = zeff - sum(charges(atom_a, 1:nspins))
927 END DO
928 END DO
929 DEALLOCATE (charges)
930
931 ALLOCATE (qpot(number_qm_atoms))
932 qpot = 0.0_dp
933 ALLOCATE (forces_qm(3, number_qm_atoms))
934 forces_qm = 0.0_dp
935
936 ! Create Ewald environments
937 poisson_section => section_vals_get_subs_vals(qs_env%input, "MM%POISSON")
938 ALLOCATE (ewald_env)
939 CALL ewald_env_create(ewald_env, para_env)
940 CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
941 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
942 CALL read_ewald_section(ewald_env, ewald_section)
943 print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
944 ALLOCATE (ewald_pw)
945 CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=print_section)
946
947 CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, do_ipol=do_ipol)
948 IF (do_multipoles) cpabort("No multipole force fields allowed in DFTB QM/MM")
949 IF (do_ipol /= do_fist_pol_none) cpabort("No polarizable force fields allowed in DFTB QM/MM")
950
951 SELECT CASE (ewald_type)
952 CASE (do_ewald_pme)
953 cpabort("PME Ewald type not implemented for DFTB/QMMM")
954 CASE (do_ewald_ewald, do_ewald_spme)
955 DO ipot = 1, SIZE(qmmm_env%Potentials)
956 pot => qmmm_env%Potentials(ipot)%Pot
957 nmm = SIZE(pot%mm_atom_index)
958 ! get a 'clean' mm particle set
959 NULLIFY (atoms_mm)
960 CALL allocate_particle_set(atoms_mm, nmm)
961 ALLOCATE (charges_mm(nmm))
962 DO imp = 1, nmm
963 imm = pot%mm_atom_index(imp)
964 indmm = qmmm_env%mm_atom_index(imm)
965 atoms_mm(imp)%r = particles_mm(indmm)%r
966 atoms_mm(imp)%atomic_kind => particles_mm(indmm)%atomic_kind
967 charges_mm(imp) = qmmm_env%mm_atom_chrg(imm)
968 END DO
969 ! force array for mm atoms
970 ALLOCATE (forces_mm(3, nmm))
971 forces_mm = 0.0_dp
972 IF (ewald_type == do_ewald_ewald) THEN
973 cpabort("Ewald not implemented for DFTB/QMMM")
974 ELSE IF (ewald_type == do_ewald_spme) THEN
975 ! spme electrostatic potential
976 CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
977 particles_qm, qpot)
978 ! forces QM
979 CALL spme_forces(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
980 particles_qm, mcharge, forces_qm)
981 ! forces MM
982 CALL spme_forces(ewald_env, ewald_pw, mm_cell, particles_qm, mcharge, &
983 atoms_mm, charges_mm, forces_mm)
984 END IF
985 CALL deallocate_particle_set(atoms_mm)
986 DEALLOCATE (charges_mm)
987 ! transfer MM forces
988 CALL para_env%sum(forces_mm)
989 DO imp = 1, nmm
990 imm = pot%mm_atom_index(imp)
991 forces(:, imm) = forces(:, imm) - forces_mm(:, imp)
992 END DO
993 DEALLOCATE (forces_mm)
994 END DO
995
996 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
997 DO ipot = 1, SIZE(qmmm_env%added_charges%Potentials)
998 pot => qmmm_env%added_charges%Potentials(ipot)%Pot
999 nmm = SIZE(pot%mm_atom_index)
1000 ! get a 'clean' mm particle set
1001 NULLIFY (atoms_mm)
1002 CALL allocate_particle_set(atoms_mm, nmm)
1003 ALLOCATE (charges_mm(nmm))
1004 DO imp = 1, nmm
1005 imm = pot%mm_atom_index(imp)
1006 indmm = qmmm_env%added_charges%mm_atom_index(imm)
1007 atoms_mm(imp)%r = qmmm_env%added_charges%added_particles(indmm)%r
1008 atoms_mm(imp)%atomic_kind => qmmm_env%added_charges%added_particles(indmm)%atomic_kind
1009 charges_mm(imp) = qmmm_env%added_charges%mm_atom_chrg(imm)
1010 END DO
1011 ! force array for mm atoms
1012 ALLOCATE (forces_mm(3, nmm))
1013 forces_mm = 0.0_dp
1014 IF (ewald_type == do_ewald_ewald) THEN
1015 cpabort("Ewald not implemented for DFTB/QMMM")
1016 ELSE IF (ewald_type == do_ewald_spme) THEN
1017 ! spme electrostatic potential
1018 CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, &
1019 charges_mm, particles_qm, qpot)
1020 ! forces QM
1021 CALL spme_forces(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
1022 particles_qm, mcharge, forces_qm)
1023 ! forces MM
1024 CALL spme_forces(ewald_env, ewald_pw, mm_cell, particles_qm, mcharge, &
1025 atoms_mm, charges_mm, forces_mm)
1026 END IF
1027 CALL deallocate_particle_set(atoms_mm)
1028 ! transfer MM forces
1029 CALL para_env%sum(forces_mm)
1030 DO imp = 1, nmm
1031 imm = pot%mm_atom_index(imp)
1032 forces_added_charges(:, imm) = forces_added_charges(:, imm) - forces_mm(:, imp)
1033 END DO
1034 DEALLOCATE (forces_mm)
1035 END DO
1036 END IF
1037 CALL para_env%sum(qpot)
1038 CALL para_env%sum(forces_qm)
1039 ! Add Ewald and DFTB short range corrections
1040 ! This is effectively using a minimum image convention!
1041 ! Set rcutoff to values compatible with alpha Ewald
1042 CALL ewald_env_get(ewald_env, rcut=rcutoff(1), alpha=alpha)
1043 rcutoff(2) = 0.025_dp*rcutoff(1)
1044 rcutoff(1) = 2.0_dp*rcutoff(1)
1045 nkind = SIZE(atomic_kind_set)
1046 iqm = 0
1047 DO ikind = 1, nkind
1048 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1049 IF (do_dftb) THEN
1050 NULLIFY (dftb_kind)
1051 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1052 CALL get_dftb_atom_param(dftb_kind, &
1053 defined=defined, eta=eta_a, natorb=natorb)
1054 ! use mm charge smearing for non-scc cases
1055 IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
1056 IF (.NOT. defined .OR. natorb < 1) cycle
1057 ELSEIF (do_xtb) THEN
1058 eta_a(0) = eta_mm
1059 END IF
1060 DO i = 1, SIZE(list)
1061 iatom = list(i)
1062 iqm = iqm + 1
1063 CALL build_mm_pot(qpot(iatom), 1, eta_a(0), qmmm_env%Potentials, particles_mm, &
1064 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1065 mm_cell, iatom, rcutoff, particles_qm)
1066 CALL build_mm_dpot(mcharge(iatom), 1, eta_a(0), qmmm_env%Potentials, particles_mm, &
1067 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1068 mm_cell, iatom, forces, forces_qm(:, iqm), &
1069 rcutoff, particles_qm)
1070 CALL build_mm_pot(qpot(iatom), 2, alpha, qmmm_env%Potentials, particles_mm, &
1071 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1072 mm_cell, iatom, rcutoff, particles_qm)
1073 CALL build_mm_dpot(mcharge(iatom), 2, alpha, qmmm_env%Potentials, particles_mm, &
1074 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1075 mm_cell, iatom, forces, forces_qm(:, iqm), &
1076 rcutoff, particles_qm)
1077 ! Possibly added charges
1078 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
1079 CALL build_mm_pot(qpot(iatom), 1, eta_a(0), qmmm_env%added_charges%potentials, &
1080 qmmm_env%added_charges%added_particles, &
1081 qmmm_env%added_charges%mm_atom_chrg, &
1082 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, rcutoff, &
1083 particles_qm)
1084 CALL build_mm_dpot( &
1085 mcharge(iatom), 1, eta_a(0), qmmm_env%added_charges%potentials, &
1086 qmmm_env%added_charges%added_particles, qmmm_env%added_charges%mm_atom_chrg, &
1087 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
1088 forces_added_charges, forces_qm(:, iqm), &
1089 rcutoff, particles_qm)
1090 CALL build_mm_pot(qpot(iatom), 2, alpha, qmmm_env%added_charges%potentials, &
1091 qmmm_env%added_charges%added_particles, &
1092 qmmm_env%added_charges%mm_atom_chrg, &
1093 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, rcutoff, &
1094 particles_qm)
1095 CALL build_mm_dpot( &
1096 mcharge(iatom), 2, alpha, qmmm_env%added_charges%potentials, &
1097 qmmm_env%added_charges%added_particles, qmmm_env%added_charges%mm_atom_chrg, &
1098 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
1099 forces_added_charges, forces_qm(:, iqm), &
1100 rcutoff, particles_qm)
1101 END IF
1102 END DO
1103 END DO
1104
1105 CASE (do_ewald_none)
1106 ! Simply summing up charges with 1/R (DFTB corrected)
1107 ! calculate potential and forces from classical charges
1108 iqm = 0
1109 DO ikind = 1, nkind
1110 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1111 IF (do_dftb) THEN
1112 NULLIFY (dftb_kind)
1113 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1114 CALL get_dftb_atom_param(dftb_kind, &
1115 defined=defined, eta=eta_a, natorb=natorb)
1116 ! use mm charge smearing for non-scc cases
1117 IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
1118 IF (.NOT. defined .OR. natorb < 1) cycle
1119 ELSEIF (do_xtb) THEN
1120 eta_a(0) = eta_mm
1121 END IF
1122 DO i = 1, SIZE(list)
1123 iatom = list(i)
1124 iqm = iqm + 1
1125 CALL build_mm_pot(qpot(iatom), 0, eta_a(0), qmmm_env%Potentials, particles_mm, &
1126 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, mm_cell, iatom, &
1127 qmmm_env%spherical_cutoff, particles_qm)
1128 CALL build_mm_dpot(mcharge(iatom), 0, eta_a(0), qmmm_env%Potentials, particles_mm, &
1129 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
1130 mm_cell, iatom, forces, forces_qm(:, iqm), &
1131 qmmm_env%spherical_cutoff, particles_qm)
1132 ! Possibly added charges
1133 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
1134 CALL build_mm_pot(qpot(iatom), 0, eta_a(0), qmmm_env%added_charges%potentials, &
1135 qmmm_env%added_charges%added_particles, &
1136 qmmm_env%added_charges%mm_atom_chrg, &
1137 qmmm_env%added_charges%mm_atom_index, &
1138 mm_cell, iatom, qmmm_env%spherical_cutoff, &
1139 particles_qm)
1140 CALL build_mm_dpot(mcharge(iatom), 0, eta_a(0), qmmm_env%added_charges%potentials, &
1141 qmmm_env%added_charges%added_particles, &
1142 qmmm_env%added_charges%mm_atom_chrg, &
1143 qmmm_env%added_charges%mm_atom_index, mm_cell, iatom, &
1144 forces_added_charges, &
1145 forces_qm(:, iqm), qmmm_env%spherical_cutoff, particles_qm)
1146 END IF
1147 END DO
1148 END DO
1149 CASE DEFAULT
1150 cpabort("Unknown Ewald type!")
1151 END SELECT
1152
1153 ! Transfer QM gradients to the QM particles..
1154 iqm = 0
1155 DO ikind = 1, nkind
1156 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1157 IF (do_dftb) THEN
1158 NULLIFY (dftb_kind)
1159 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1160 CALL get_dftb_atom_param(dftb_kind, defined=defined, natorb=natorb)
1161 IF (.NOT. defined .OR. natorb < 1) cycle
1162 ELSEIF (do_xtb) THEN
1163 !
1164 END IF
1165 DO i = 1, SIZE(list)
1166 iqm = iqm + 1
1167 iatom = qmmm_env%qm_atom_index(list(i))
1168 particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + forces_qm(:, iqm)
1169 END DO
1170 END DO
1171
1172 ! derivatives from qm charges
1173 forces_qm = 0.0_dp
1174 IF (SIZE(matrix_p) == 2) THEN
1175 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
1176 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1177 END IF
1178 !
1179 CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
1180 DO WHILE (dbcsr_iterator_blocks_left(iter))
1181 CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock)
1182 !
1183 IF (iatom == jatom) cycle
1184 !
1185 gmij = -0.5_dp*(qpot(iatom) + qpot(jatom))
1186 NULLIFY (pblock)
1187 CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
1188 row=iatom, col=jatom, block=pblock, found=found)
1189 cpassert(found)
1190 DO i = 1, 3
1191 NULLIFY (dsblock)
1192 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
1193 row=iatom, col=jatom, block=dsblock, found=found)
1194 cpassert(found)
1195 fi = -2.0_dp*gmij*sum(pblock*dsblock)
1196 forces_qm(i, iatom) = forces_qm(i, iatom) + fi
1197 forces_qm(i, jatom) = forces_qm(i, jatom) - fi
1198 END DO
1199 END DO
1200 CALL dbcsr_iterator_stop(iter)
1201 !
1202 IF (SIZE(matrix_p) == 2) THEN
1203 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
1204 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1205 END IF
1206 !
1207 ! Transfer QM gradients to the QM particles..
1208 CALL para_env%sum(forces_qm)
1209 DO ikind = 1, nkind
1210 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1211 DO i = 1, SIZE(list)
1212 iqm = list(i)
1213 iatom = qmmm_env%qm_atom_index(iqm)
1214 particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + forces_qm(:, iqm)
1215 END DO
1216 END DO
1217 !
1218 DEALLOCATE (mcharge)
1219 !
1220 ! MM forces will be handled directly from the QMMM module in the same way
1221 ! as for GPW/GAPW methods
1222 DEALLOCATE (forces_qm)
1223 DEALLOCATE (qpot)
1224
1225 ! Release Ewald environment
1226 CALL ewald_env_release(ewald_env)
1227 DEALLOCATE (ewald_env)
1228 CALL ewald_pw_release(ewald_pw)
1229 DEALLOCATE (ewald_pw)
1230
1231 CALL dbcsr_deallocate_matrix_set(matrix_s)
1232
1233 END IF
1234
1235 CALL timestop(handle)
1236
1237 END SUBROUTINE deriv_tb_qmmm_matrix_pc
1238
1239! **************************************************************************************************
1240!> \brief ...
1241!> \param qpot ...
1242!> \param pot_type ...
1243!> \param qm_alpha ...
1244!> \param potentials ...
1245!> \param particles_mm ...
1246!> \param mm_charges ...
1247!> \param mm_atom_index ...
1248!> \param mm_cell ...
1249!> \param IndQM ...
1250!> \param qmmm_spherical_cutoff ...
1251!> \param particles_qm ...
1252! **************************************************************************************************
1253 SUBROUTINE build_mm_pot(qpot, pot_type, qm_alpha, potentials, &
1254 particles_mm, mm_charges, mm_atom_index, mm_cell, IndQM, &
1255 qmmm_spherical_cutoff, particles_qm)
1256
1257 REAL(kind=dp), INTENT(INOUT) :: qpot
1258 INTEGER, INTENT(IN) :: pot_type
1259 REAL(kind=dp), INTENT(IN) :: qm_alpha
1260 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1261 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
1262 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1263 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1264 TYPE(cell_type), POINTER :: mm_cell
1265 INTEGER, INTENT(IN) :: indqm
1266 REAL(kind=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
1267 TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
1268
1269 CHARACTER(len=*), PARAMETER :: routinen = 'build_mm_pot'
1270 REAL(kind=dp), PARAMETER :: qsmall = 1.0e-15_dp
1271
1272 INTEGER :: handle, imm, imp, indmm, ipot
1273 REAL(kind=dp) :: dr, qeff, rt1, rt2, rt3, &
1274 sph_chrg_factor, sr
1275 REAL(kind=dp), DIMENSION(3) :: r_pbc, rij
1276 TYPE(qmmm_pot_type), POINTER :: pot
1277
1278 CALL timeset(routinen, handle)
1279 ! Loop Over MM atoms
1280 ! Loop over Pot stores atoms with the same charge
1281 mainlooppot: DO ipot = 1, SIZE(potentials)
1282 pot => potentials(ipot)%Pot
1283 ! Loop over atoms belonging to this type
1284 loopmm: DO imp = 1, SIZE(pot%mm_atom_index)
1285 imm = pot%mm_atom_index(imp)
1286 indmm = mm_atom_index(imm)
1287 r_pbc = pbc(particles_mm(indmm)%r - particles_qm(indqm)%r, mm_cell)
1288 rt1 = r_pbc(1)
1289 rt2 = r_pbc(2)
1290 rt3 = r_pbc(3)
1291 rij = (/rt1, rt2, rt3/)
1292 dr = sqrt(sum(rij**2))
1293 qeff = mm_charges(imm)
1294 ! Computes the screening factor for the spherical cutoff (if defined)
1295 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1296 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
1297 qeff = qeff*sph_chrg_factor
1298 END IF
1299 IF (abs(qeff) <= qsmall) cycle
1300 IF (dr > rtiny) THEN
1301 IF (pot_type == 0) THEN
1302 sr = gamma_rab_sr(dr, qm_alpha, eta_mm, 0.0_dp)
1303 qpot = qpot + qeff*(1.0_dp/dr - sr)
1304 ELSE IF (pot_type == 1) THEN
1305 sr = gamma_rab_sr(dr, qm_alpha, eta_mm, 0.0_dp)
1306 qpot = qpot - qeff*sr
1307 ELSE IF (pot_type == 2) THEN
1308 sr = erfc(qm_alpha*dr)/dr
1309 qpot = qpot + qeff*sr
1310 ELSE
1311 cpabort("")
1312 END IF
1313 END IF
1314 END DO loopmm
1315 END DO mainlooppot
1316 CALL timestop(handle)
1317 END SUBROUTINE build_mm_pot
1318
1319! **************************************************************************************************
1320!> \brief ...
1321!> \param qcharge ...
1322!> \param pot_type ...
1323!> \param qm_alpha ...
1324!> \param potentials ...
1325!> \param particles_mm ...
1326!> \param mm_charges ...
1327!> \param mm_atom_index ...
1328!> \param mm_cell ...
1329!> \param IndQM ...
1330!> \param forces ...
1331!> \param forces_qm ...
1332!> \param qmmm_spherical_cutoff ...
1333!> \param particles_qm ...
1334! **************************************************************************************************
1335 SUBROUTINE build_mm_dpot(qcharge, pot_type, qm_alpha, potentials, &
1336 particles_mm, mm_charges, mm_atom_index, mm_cell, IndQM, &
1337 forces, forces_qm, qmmm_spherical_cutoff, particles_qm)
1338
1339 REAL(kind=dp), INTENT(IN) :: qcharge
1340 INTEGER, INTENT(IN) :: pot_type
1341 REAL(kind=dp), INTENT(IN) :: qm_alpha
1342 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1343 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
1344 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1345 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1346 TYPE(cell_type), POINTER :: mm_cell
1347 INTEGER, INTENT(IN) :: indqm
1348 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
1349 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: forces_qm
1350 REAL(kind=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
1351 TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
1352
1353 CHARACTER(len=*), PARAMETER :: routinen = 'build_mm_dpot'
1354 REAL(kind=dp), PARAMETER :: qsmall = 1.0e-15_dp
1355
1356 INTEGER :: handle, imm, imp, indmm, ipot
1357 REAL(kind=dp) :: dr, drm, drp, dsr, fsr, qeff, rt1, rt2, &
1358 rt3, sph_chrg_factor
1359 REAL(kind=dp), DIMENSION(3) :: force_ab, r_pbc, rij
1360 TYPE(qmmm_pot_type), POINTER :: pot
1361
1362 CALL timeset(routinen, handle)
1363 ! Loop Over MM atoms
1364 ! Loop over Pot stores atoms with the same charge
1365 mainlooppot: DO ipot = 1, SIZE(potentials)
1366 pot => potentials(ipot)%Pot
1367 ! Loop over atoms belonging to this type
1368 loopmm: DO imp = 1, SIZE(pot%mm_atom_index)
1369 imm = pot%mm_atom_index(imp)
1370 indmm = mm_atom_index(imm)
1371 r_pbc = pbc(particles_mm(indmm)%r - particles_qm(indqm)%r, mm_cell)
1372 rt1 = r_pbc(1)
1373 rt2 = r_pbc(2)
1374 rt3 = r_pbc(3)
1375 rij = (/rt1, rt2, rt3/)
1376 dr = sqrt(sum(rij**2))
1377 qeff = mm_charges(imm)
1378 ! Computes the screening factor for the spherical cutoff (if defined)
1379 ! We neglect derivative of cutoff function for gradients!!!
1380 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1381 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
1382 qeff = qeff*sph_chrg_factor
1383 END IF
1384 IF (abs(qeff) <= qsmall) cycle
1385 IF (dr > rtiny) THEN
1386 drp = dr + ddrmm
1387 drm = dr - ddrmm
1388 IF (pot_type == 0) THEN
1389 dsr = 0.5_dp*(gamma_rab_sr(drp, qm_alpha, eta_mm, 0.0_dp) - &
1390 gamma_rab_sr(drm, qm_alpha, eta_mm, 0.0_dp))/ddrmm
1391 fsr = qeff*qcharge*(-1.0_dp/(dr*dr) - dsr)
1392 ELSE IF (pot_type == 1) THEN
1393 dsr = 0.5_dp*(gamma_rab_sr(drp, qm_alpha, eta_mm, 0.0_dp) - &
1394 gamma_rab_sr(drm, qm_alpha, eta_mm, 0.0_dp))/ddrmm
1395 fsr = -qeff*qcharge*dsr
1396 ELSE IF (pot_type == 2) THEN
1397 dsr = 0.5_dp*(erfc(qm_alpha*drp)/drp - erfc(qm_alpha*drm)/drm)/ddrmm
1398 fsr = qeff*qcharge*dsr
1399 ELSE
1400 cpabort("")
1401 END IF
1402 force_ab = -fsr*rij/dr
1403 ELSE
1404 force_ab = 0.0_dp
1405 END IF
1406 ! The array of QM forces are really the forces
1407 forces_qm(:) = forces_qm(:) - force_ab
1408 ! The one of MM atoms are instead gradients
1409 forces(:, imm) = forces(:, imm) - force_ab
1410 END DO loopmm
1411 END DO mainlooppot
1412
1413 CALL timestop(handle)
1414
1415 END SUBROUTINE build_mm_dpot
1416
1417END MODULE qmmm_tb_methods
1418
cell_types::pbc
Definition cell_types.F:92
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
mulliken::mulliken_charges
Definition mulliken.F:36
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:956
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:930
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1027
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:679
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:989
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
ewald_environment_types
Definition ewald_environment_types.F:14
ewald_environment_types::ewald_env_set
subroutine, public ewald_env_set(ewald_env, ewald_type, alpha, epsilon, eps_pol, gmax, ns_max, precs, o_spline, para_env, poisson_section, interaction_cutoffs, cell_hmat)
Purpose: Set the EWALD environment.
Definition ewald_environment_types.F:181
ewald_environment_types::ewald_env_create
subroutine, public ewald_env_create(ewald_env, para_env)
allocates and intitializes a ewald_env
Definition ewald_environment_types.F:222
ewald_environment_types::read_ewald_section
subroutine, public read_ewald_section(ewald_env, ewald_section)
Purpose: read the EWALD section.
Definition ewald_environment_types.F:256
ewald_environment_types::ewald_env_release
subroutine, public ewald_env_release(ewald_env)
releases the given ewald_env (see doc/ReferenceCounting.html)
Definition ewald_environment_types.F:239
ewald_environment_types::ewald_env_get
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Definition ewald_environment_types.F:124
ewald_pw_types
pw_types
Definition ewald_pw_types.F:12
ewald_pw_types::ewald_pw_release
subroutine, public ewald_pw_release(ewald_pw)
releases the memory used by the ewald_pw
Definition ewald_pw_types.F:102
ewald_pw_types::ewald_pw_create
subroutine, public ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section)
creates the structure ewald_pw_type
Definition ewald_pw_types.F:83
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_fist_pol_none
integer, parameter, public do_fist_pol_none
Definition input_constants.F:954
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mulliken
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition mulliken.F:13
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
particle_types::deallocate_particle_set
subroutine, public deallocate_particle_set(particle_set)
Deallocate a particle set.
Definition particle_types.F:87
particle_types::allocate_particle_set
subroutine, public allocate_particle_set(particle_set, nparticle)
Allocate a particle set.
Definition particle_types.F:69
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_poisson_types::do_ewald_pme
integer, parameter, public do_ewald_pme
Definition pw_poisson_types.F:75
pw_poisson_types::do_ewald_ewald
integer, parameter, public do_ewald_ewald
Definition pw_poisson_types.F:75
pw_poisson_types::do_ewald_none
integer, parameter, public do_ewald_none
Definition pw_poisson_types.F:75
pw_poisson_types::do_ewald_spme
integer, parameter, public do_ewald_spme
Definition pw_poisson_types.F:75
qmmm_tb_methods
TB methods used with QMMM.
Definition qmmm_tb_methods.F:12
qmmm_tb_methods::build_tb_qmmm_matrix_zero
subroutine, public build_tb_qmmm_matrix_zero(qs_env, para_env)
Constructs an empty 1-el DFTB hamiltonian.
Definition qmmm_tb_methods.F:244
qmmm_tb_methods::build_tb_qmmm_matrix_pc
subroutine, public build_tb_qmmm_matrix_pc(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
Constructs the 1-el DFTB hamiltonian.
Definition qmmm_tb_methods.F:310
qmmm_tb_methods::build_tb_qmmm_matrix
subroutine, public build_tb_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
Constructs the 1-el DFTB hamiltonian.
Definition qmmm_tb_methods.F:102
qmmm_tb_methods::rtiny
real(dp), parameter rtiny
Definition qmmm_tb_methods.F:76
qmmm_tb_methods::deriv_tb_qmmm_matrix
subroutine, public deriv_tb_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, forces, forces_added_charges)
Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms.
Definition qmmm_tb_methods.F:593
qmmm_tb_methods::eta_mm
real(dp), parameter eta_mm
Definition qmmm_tb_methods.F:78
qmmm_tb_methods::deriv_tb_qmmm_matrix_pc
subroutine, public deriv_tb_qmmm_matrix_pc(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, forces, forces_added_charges)
Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms.
Definition qmmm_tb_methods.F:830
qmmm_tb_methods::ddrmm
real(dp), parameter ddrmm
Definition qmmm_tb_methods.F:80
qmmm_types_low
Definition qmmm_types_low.F:13
qmmm_util
Definition qmmm_util.F:13
qmmm_util::spherical_cutoff_factor
subroutine, public spherical_cutoff_factor(spherical_cutoff, rij, factor)
Computes a spherical cutoff factor for the QMMM interactions.
Definition qmmm_util.F:615
qs_dftb_coulomb
Calculation of Coulomb contributions in DFTB.
Definition qs_dftb_coulomb.F:12
qs_dftb_coulomb::gamma_rab_sr
real(dp) function, public gamma_rab_sr(r, ga, gb, hb_para)
Computes the short-range gamma parameter from exact Coulomb interaction of normalized exp(-a*r) charg...
Definition qs_dftb_coulomb.F:550
qs_dftb_matrices
Calculation of Overlap and Hamiltonian matrices in DFTB.
Definition qs_dftb_matrices.F:12
qs_dftb_matrices::build_dftb_overlap
subroutine, public build_dftb_overlap(qs_env, nderivative, matrix_s)
...
Definition qs_dftb_matrices.F:655
qs_dftb_types
Definition of the DFTB parameter types.
Definition qs_dftb_types.F:12
qs_dftb_utils
Working with the DFTB parameter types.
Definition qs_dftb_utils.F:12
qs_dftb_utils::get_dftb_atom_param
subroutine, public get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
Definition qs_dftb_utils.F:123
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
qs_ks_qmmm_types
Definition qs_ks_qmmm_types.F:13
qs_ks_types
Definition qs_ks_types.F:15
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_lists
Generate the atomic neighbor lists.
Definition qs_neighbor_lists.F:19
qs_neighbor_lists::build_qs_neighbor_lists
subroutine, public build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_section)
Build all the required neighbor lists for Quickstep.
Definition qs_neighbor_lists.F:286
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
qs_overlap::build_overlap_matrix
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:120
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
spme
Calculate the electrostatic energy by the Smooth Particle Ewald method.
Definition spme.F:14
spme::spme_forces
subroutine, public spme_forces(ewald_env, ewald_pw, box, particle_set_a, charges_a, particle_set_b, charges_b, forces_b)
Calculate the forces on particles B for the electrostatic interaction betrween particles A and B.
Definition spme.F:568
spme::spme_potential
subroutine, public spme_potential(ewald_env, ewald_pw, box, particle_set_a, charges_a, particle_set_b, potential)
Calculate the electrostatic potential from particles A (charge A) at positions of particles B.
Definition spme.F:427
xtb_types
Definition of the xTB parameter types.
Definition xtb_types.F:20
xtb_types::get_xtb_atom_param
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, xi, kappa0, alpg, occupation, electronegativity, chmax, en, kqat2, kcn, kq)
...
Definition xtb_types.F:199
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_control_types::dftb_control_type
Definition cp_control_types.F:111
cp_control_types::xtb_control_type
Definition cp_control_types.F:138
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:186
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:56
ewald_pw_types::ewald_pw_type
Definition ewald_pw_types.F:63
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qmmm_types_low::qmmm_env_qm_type
...
Definition qmmm_types_low.F:195
qmmm_types_low::qmmm_pot_p_type
Definition qmmm_types_low.F:83
qmmm_types_low::qmmm_pot_type
Real Space Potential.
Definition qmmm_types_low.F:76
qs_dftb_types::qs_dftb_atom_type
Definition qs_dftb_types.F:27
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_ks_qmmm_types::qs_ks_qmmm_env_type
calculation environment to calculate the ks_qmmm matrix, holds the QM/MM potential and all the needed...
Definition qs_ks_qmmm_types.F:47
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:136
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
xtb_types::xtb_atom_type
Definition xtb_types.F:42