(git:e7e05ae)
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-2024 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_operations, ONLY: dbcsr_allocate_matrix_set,&
21  dbcsr_deallocate_matrix_set
22  USE dbcsr_api, ONLY: &
23  dbcsr_add, dbcsr_copy, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
24  dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
25  dbcsr_p_type, dbcsr_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 
90 CONTAINS
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 :: blk, handle, i, iatom, ikind, jatom, &
112  natom, 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, blk)
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 :: blk, do_ipol, ewald_type, handle, i, &
320  iatom, ikind, imm, imp, indmm, ipot, &
321  jatom, 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, blk)
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, blk, handle, i, iatom, ikind, &
605  iqm, jatom, natom, natorb, nkind, &
606  nspins, 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, blk)
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, blk, do_ipol, ewald_type, handle, i, iatom, ikind, imm, imp, indmm, ipot, &
842  iqm, 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, blk)
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 
1417 END MODULE qmmm_tb_methods
1418 
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
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_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:180
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:221
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:255
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:238
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:123
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:103
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:84
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:941
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:724
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::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::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::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::eta_mm
real(dp), parameter eta_mm
Definition: qmmm_tb_methods.F:78
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:564
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:663
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
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, 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_r3d_rs_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_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
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:287
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:608
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:467
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, occupation, electronegativity, chmax)
...
Definition: xtb_types.F:175