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