Line data Source code
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 448 : 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 448 : 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 448 : REAL(KIND=dp), DIMENSION(:), POINTER :: qpot
118 448 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hblock, sblock
119 448 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
120 : TYPE(dbcsr_iterator_type) :: iter
121 448 : 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 448 : POINTER :: sab_nl
126 448 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
127 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
128 448 : 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 448 : 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 448 : natom=natom)
144 448 : dftb_control => dft_control%qs_control%dftb_control
145 448 : xtb_control => dft_control%qs_control%xtb_control
146 :
147 448 : IF (dft_control%qs_control%dftb) THEN
148 : do_dftb = .TRUE.
149 : do_xtb = .FALSE.
150 224 : ELSEIF (dft_control%qs_control%xtb) THEN
151 : do_dftb = .FALSE.
152 : do_xtb = .TRUE.
153 : ELSE
154 0 : CPABORT("TB method unknown")
155 : END IF
156 :
157 448 : CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
158 :
159 448 : NULLIFY (matrix_s)
160 448 : IF (do_dftb) THEN
161 224 : CALL build_dftb_overlap(qs_env, 0, matrix_s)
162 224 : ELSEIF (do_xtb) THEN
163 224 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
164 224 : 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 1344 : ALLOCATE (qpot(natom))
168 1792 : qpot = 0.0_dp
169 448 : pc_ener = 0.0_dp
170 :
171 448 : nkind = SIZE(atomic_kind_set)
172 1344 : DO ikind = 1, nkind
173 896 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
174 896 : IF (do_dftb) THEN
175 448 : NULLIFY (dftb_kind)
176 448 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
177 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff, &
178 448 : defined=defined, eta=eta_a, natorb=natorb)
179 : ! use mm charge smearing for non-scc cases
180 448 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
181 448 : IF (.NOT. defined .OR. natorb < 1) CYCLE
182 448 : ELSEIF (do_xtb) THEN
183 448 : NULLIFY (xtb_kind)
184 448 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
185 448 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
186 448 : eta_a(0) = eta_mm
187 : END IF
188 2688 : DO i = 1, SIZE(list)
189 1344 : 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 1344 : qmmm_env%spherical_cutoff, particles_qm)
193 : ! Possibly added charges
194 1344 : 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 0 : particles_qm)
201 : END IF
202 2240 : pc_ener = pc_ener + qpot(iatom)*zeff
203 : END DO
204 : END DO
205 :
206 : ! Allocate the core Hamiltonian matrix
207 448 : CALL get_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env_loc)
208 448 : matrix_h => ks_qmmm_env_loc%matrix_h
209 448 : CALL dbcsr_allocate_matrix_set(matrix_h, 1)
210 448 : ALLOCATE (matrix_h(1)%matrix)
211 : CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, &
212 448 : name="QMMM HAMILTONIAN MATRIX")
213 448 : CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
214 :
215 448 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
216 1792 : DO WHILE (dbcsr_iterator_blocks_left(iter))
217 1344 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock, blk)
218 1344 : NULLIFY (hblock)
219 : CALL dbcsr_get_block_p(matrix=matrix_h(1)%matrix, &
220 1344 : row=iatom, col=jatom, block=hblock, found=found)
221 1344 : CPASSERT(found)
222 25088 : hblock = hblock - 0.5_dp*sblock*(qpot(iatom) + qpot(jatom))
223 : END DO
224 448 : CALL dbcsr_iterator_stop(iter)
225 :
226 448 : ks_qmmm_env_loc%matrix_h => matrix_h
227 448 : ks_qmmm_env_loc%pc_ener = pc_ener
228 :
229 448 : DEALLOCATE (qpot)
230 :
231 448 : CALL dbcsr_deallocate_matrix_set(matrix_s)
232 :
233 448 : CALL timestop(handle)
234 :
235 448 : 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 8 : 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 8 : 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 8 : 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 8 : CALL timeset(routineN, handle)
260 :
261 8 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
262 :
263 8 : IF (dft_control%qs_control%dftb) THEN
264 : do_dftb = .TRUE.
265 : do_xtb = .FALSE.
266 4 : ELSEIF (dft_control%qs_control%xtb) THEN
267 : do_dftb = .FALSE.
268 : do_xtb = .TRUE.
269 : ELSE
270 0 : CPABORT("TB method unknown")
271 : END IF
272 :
273 8 : CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
274 :
275 8 : NULLIFY (matrix_s)
276 8 : IF (do_dftb) THEN
277 4 : CALL build_dftb_overlap(qs_env, 0, matrix_s)
278 4 : ELSEIF (do_xtb) THEN
279 4 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
280 4 : 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 8 : CALL get_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env_loc)
285 8 : matrix_h => ks_qmmm_env_loc%matrix_h
286 8 : CALL dbcsr_allocate_matrix_set(matrix_h, 1)
287 8 : ALLOCATE (matrix_h(1)%matrix)
288 : CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, &
289 8 : name="QMMM HAMILTONIAN MATRIX")
290 8 : CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
291 8 : ks_qmmm_env_loc%matrix_h => matrix_h
292 8 : ks_qmmm_env_loc%pc_ener = 0.0_dp
293 :
294 8 : CALL dbcsr_deallocate_matrix_set(matrix_s)
295 :
296 8 : CALL timestop(handle)
297 :
298 8 : 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 1116 : 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 1116 : 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 1116 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges_mm, qpot
329 1116 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hblock, sblock
330 1116 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
331 : TYPE(dbcsr_iterator_type) :: iter
332 1116 : 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 1116 : POINTER :: sab_nl
339 1116 : 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 1116 : 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 1116 : 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 1116 : natom=natom)
360 1116 : dftb_control => dft_control%qs_control%dftb_control
361 1116 : xtb_control => dft_control%qs_control%xtb_control
362 :
363 1116 : IF (dft_control%qs_control%dftb) THEN
364 : do_dftb = .TRUE.
365 : do_xtb = .FALSE.
366 558 : ELSEIF (dft_control%qs_control%xtb) THEN
367 : do_dftb = .FALSE.
368 : do_xtb = .TRUE.
369 : ELSE
370 0 : CPABORT("TB method unknown")
371 : END IF
372 :
373 1116 : CALL build_qs_neighbor_lists(qs_env, para_env, force_env_section=qs_env%input)
374 :
375 1116 : NULLIFY (matrix_s)
376 1116 : IF (do_dftb) THEN
377 558 : CALL build_dftb_overlap(qs_env, 0, matrix_s)
378 558 : ELSEIF (do_xtb) THEN
379 558 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_nl)
380 558 : 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 3348 : ALLOCATE (qpot(natom))
384 4464 : qpot = 0.0_dp
385 1116 : pc_ener = 0.0_dp
386 :
387 : ! Create Ewald environments
388 1116 : poisson_section => section_vals_get_subs_vals(qs_env%input, "MM%POISSON")
389 14508 : ALLOCATE (ewald_env)
390 1116 : CALL ewald_env_create(ewald_env, para_env)
391 1116 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
392 1116 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
393 1116 : CALL read_ewald_section(ewald_env, ewald_section)
394 1116 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
395 1116 : ALLOCATE (ewald_pw)
396 1116 : CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=print_section)
397 :
398 1116 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, do_ipol=do_ipol)
399 1116 : IF (do_multipoles) CPABORT("No multipole force fields allowed in TB QM/MM")
400 1116 : IF (do_ipol /= do_fist_pol_none) CPABORT("No polarizable force fields allowed in TB QM/MM")
401 :
402 0 : SELECT CASE (ewald_type)
403 : CASE (do_ewald_pme)
404 0 : CPABORT("PME Ewald type not implemented for TB/QMMM")
405 : CASE (do_ewald_ewald, do_ewald_spme)
406 2004 : DO ipot = 1, SIZE(qmmm_env%Potentials)
407 1332 : Pot => qmmm_env%Potentials(ipot)%Pot
408 1332 : nmm = SIZE(Pot%mm_atom_index)
409 : ! get a 'clean' mm particle set
410 1332 : NULLIFY (atoms_mm)
411 1332 : CALL allocate_particle_set(atoms_mm, nmm)
412 3996 : ALLOCATE (charges_mm(nmm))
413 5328 : DO Imp = 1, nmm
414 3996 : Imm = Pot%mm_atom_index(Imp)
415 3996 : IndMM = qmmm_env%mm_atom_index(Imm)
416 27972 : atoms_mm(Imp)%r = particles_mm(IndMM)%r
417 3996 : atoms_mm(Imp)%atomic_kind => particles_mm(IndMM)%atomic_kind
418 5328 : charges_mm(Imp) = qmmm_env%mm_atom_chrg(Imm)
419 : END DO
420 1332 : IF (ewald_type == do_ewald_ewald) THEN
421 0 : CPABORT("Ewald not implemented for TB/QMMM")
422 1332 : ELSE IF (ewald_type == do_ewald_spme) THEN
423 : ! spme electrostatic potential
424 1332 : CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, particles_qm, qpot)
425 : END IF
426 1332 : CALL deallocate_particle_set(atoms_mm)
427 2004 : DEALLOCATE (charges_mm)
428 : END DO
429 672 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
430 0 : DO ipot = 1, SIZE(qmmm_env%added_charges%Potentials)
431 0 : Pot => qmmm_env%added_charges%Potentials(ipot)%Pot
432 0 : nmm = SIZE(Pot%mm_atom_index)
433 : ! get a 'clean' mm particle set
434 0 : NULLIFY (atoms_mm)
435 0 : CALL allocate_particle_set(atoms_mm, nmm)
436 0 : ALLOCATE (charges_mm(nmm))
437 0 : DO Imp = 1, nmm
438 0 : Imm = Pot%mm_atom_index(Imp)
439 0 : IndMM = qmmm_env%added_charges%mm_atom_index(Imm)
440 0 : atoms_mm(Imp)%r = qmmm_env%added_charges%added_particles(IndMM)%r
441 0 : atoms_mm(Imp)%atomic_kind => qmmm_env%added_charges%added_particles(IndMM)%atomic_kind
442 0 : charges_mm(Imp) = qmmm_env%added_charges%mm_atom_chrg(Imm)
443 : END DO
444 0 : IF (ewald_type == do_ewald_ewald) THEN
445 0 : CPABORT("Ewald not implemented for TB/QMMM")
446 0 : ELSE IF (ewald_type == do_ewald_spme) THEN
447 : ! spme electrostatic potential
448 0 : CALL spme_potential(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, particles_qm, qpot)
449 : END IF
450 0 : CALL deallocate_particle_set(atoms_mm)
451 672 : DEALLOCATE (charges_mm)
452 : END DO
453 : END IF
454 4704 : 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 672 : CALL ewald_env_get(ewald_env, rcut=rcutoff(1), alpha=alpha)
459 672 : rcutoff(2) = 0.025_dp*rcutoff(1)
460 672 : rcutoff(1) = 2.0_dp*rcutoff(1)
461 672 : nkind = SIZE(atomic_kind_set)
462 2016 : DO ikind = 1, nkind
463 1344 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
464 1344 : IF (do_dftb) THEN
465 672 : NULLIFY (dftb_kind)
466 672 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
467 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff, &
468 672 : defined=defined, eta=eta_a, natorb=natorb)
469 : ! use mm charge smearing for non-scc cases
470 672 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
471 672 : IF (.NOT. defined .OR. natorb < 1) CYCLE
472 672 : ELSEIF (do_xtb) THEN
473 672 : NULLIFY (xtb_kind)
474 672 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
475 672 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
476 672 : eta_a(0) = eta_mm
477 : END IF
478 4032 : DO i = 1, SIZE(list)
479 2016 : 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 2016 : 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 2016 : particles_qm)
486 : ! Possibly added charges
487 2016 : 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 0 : 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 0 : particles_qm)
498 : END IF
499 3360 : 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 444 : nkind = SIZE(atomic_kind_set)
505 1332 : DO ikind = 1, nkind
506 888 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
507 888 : IF (do_dftb) THEN
508 444 : NULLIFY (dftb_kind)
509 444 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
510 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff, &
511 444 : defined=defined, eta=eta_a, natorb=natorb)
512 : ! use mm charge smearing for non-scc cases
513 444 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
514 444 : IF (.NOT. defined .OR. natorb < 1) CYCLE
515 444 : ELSEIF (do_xtb) THEN
516 444 : NULLIFY (xtb_kind)
517 444 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
518 444 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
519 444 : eta_a(0) = eta_mm
520 : END IF
521 2664 : DO i = 1, SIZE(list)
522 1332 : 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 1332 : qmmm_env%spherical_cutoff, particles_qm)
526 : ! Possibly added charges
527 1332 : 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 0 : particles_qm)
534 : END IF
535 2220 : pc_ener = pc_ener + qpot(iatom)*zeff
536 : END DO
537 : END DO
538 : CASE DEFAULT
539 1116 : CPABORT("Unknown Ewald type!")
540 : END SELECT
541 :
542 : ! Allocate the core Hamiltonian matrix
543 1116 : CALL get_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env_loc)
544 1116 : matrix_h => ks_qmmm_env_loc%matrix_h
545 1116 : CALL dbcsr_allocate_matrix_set(matrix_h, 1)
546 1116 : ALLOCATE (matrix_h(1)%matrix)
547 : CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, &
548 1116 : name="QMMM HAMILTONIAN MATRIX")
549 1116 : CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
550 :
551 1116 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
552 4464 : DO WHILE (dbcsr_iterator_blocks_left(iter))
553 3348 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock, blk)
554 3348 : NULLIFY (hblock)
555 : CALL dbcsr_get_block_p(matrix=matrix_h(1)%matrix, &
556 3348 : row=iatom, col=jatom, block=hblock, found=found)
557 3348 : CPASSERT(found)
558 62496 : hblock = hblock - 0.5_dp*sblock*(qpot(iatom) + qpot(jatom))
559 : END DO
560 1116 : CALL dbcsr_iterator_stop(iter)
561 :
562 1116 : ks_qmmm_env_loc%matrix_h => matrix_h
563 1116 : ks_qmmm_env_loc%pc_ener = pc_ener
564 :
565 1116 : DEALLOCATE (qpot)
566 :
567 : ! Release Ewald environment
568 1116 : CALL ewald_env_release(ewald_env)
569 1116 : DEALLOCATE (ewald_env)
570 1116 : CALL ewald_pw_release(ewald_pw)
571 1116 : DEALLOCATE (ewald_pw)
572 :
573 1116 : CALL dbcsr_deallocate_matrix_set(matrix_s)
574 :
575 1116 : CALL timestop(handle)
576 :
577 2232 : 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 448 : 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 448 : 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 448 : REAL(KIND=dp), DIMENSION(:), POINTER :: mcharge, qpot
612 448 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges, dsblock, Forces_QM, pblock, &
613 448 : sblock
614 448 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
615 : TYPE(dbcsr_iterator_type) :: iter
616 448 : 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 448 : POINTER :: sab_nl
621 448 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
622 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
623 448 : 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 448 : CALL timeset(routineN, handle)
631 448 : IF (calc_force) THEN
632 16 : 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 16 : natom=number_qm_atoms)
641 16 : dftb_control => dft_control%qs_control%dftb_control
642 16 : xtb_control => dft_control%qs_control%xtb_control
643 :
644 16 : IF (dft_control%qs_control%dftb) THEN
645 : do_dftb = .TRUE.
646 : do_xtb = .FALSE.
647 8 : ELSEIF (dft_control%qs_control%xtb) THEN
648 : do_dftb = .FALSE.
649 : do_xtb = .TRUE.
650 : ELSE
651 0 : CPABORT("TB method unknown")
652 : END IF
653 :
654 16 : NULLIFY (matrix_s)
655 16 : IF (do_dftb) THEN
656 8 : CALL build_dftb_overlap(qs_env, 1, matrix_s)
657 8 : ELSEIF (do_xtb) THEN
658 8 : 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 8 : basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
661 : END IF
662 :
663 16 : CALL qs_rho_get(rho, rho_ao=matrix_p)
664 :
665 16 : nspins = dft_control%nspins
666 16 : nkind = SIZE(atomic_kind_set)
667 : ! Mulliken charges
668 64 : ALLOCATE (charges(number_qm_atoms, nspins))
669 : !
670 16 : CALL mulliken_charges(matrix_p, matrix_s(1)%matrix, para_env, charges)
671 : !
672 48 : ALLOCATE (mcharge(number_qm_atoms))
673 48 : DO ikind = 1, nkind
674 32 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
675 32 : IF (do_dftb) THEN
676 16 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
677 16 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
678 16 : ELSEIF (do_xtb) THEN
679 16 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
680 16 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
681 : END IF
682 128 : DO iatom = 1, natom
683 48 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
684 128 : mcharge(atom_a) = zeff - SUM(charges(atom_a, 1:nspins))
685 : END DO
686 : END DO
687 16 : DEALLOCATE (charges)
688 :
689 48 : ALLOCATE (qpot(number_qm_atoms))
690 64 : qpot = 0.0_dp
691 48 : ALLOCATE (Forces_QM(3, number_qm_atoms))
692 208 : Forces_QM = 0.0_dp
693 :
694 : ! calculate potential and forces from classical charges
695 : iqm = 0
696 48 : DO ikind = 1, nkind
697 32 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
698 32 : IF (do_dftb) THEN
699 16 : NULLIFY (dftb_kind)
700 16 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
701 : CALL get_dftb_atom_param(dftb_kind, &
702 16 : defined=defined, eta=eta_a, natorb=natorb)
703 : ! use mm charge smearing for non-scc cases
704 16 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
705 16 : IF (.NOT. defined .OR. natorb < 1) CYCLE
706 16 : ELSEIF (do_xtb) THEN
707 16 : eta_a(0) = eta_mm
708 : END IF
709 96 : DO i = 1, SIZE(list)
710 48 : iatom = list(i)
711 48 : 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 48 : 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 48 : qmmm_env%spherical_cutoff, particles_qm)
719 : ! Possibly added charges
720 80 : 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 0 : 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 0 : 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 48 : DO ikind = 1, nkind
740 32 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
741 32 : IF (do_dftb) THEN
742 16 : NULLIFY (dftb_kind)
743 16 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
744 16 : CALL get_dftb_atom_param(dftb_kind, defined=defined, natorb=natorb)
745 16 : IF (.NOT. defined .OR. natorb < 1) CYCLE
746 : ELSEIF (do_xtb) THEN
747 : ! use all kinds
748 : END IF
749 96 : DO i = 1, SIZE(list)
750 48 : iqm = iqm + 1
751 48 : iatom = qmmm_env%qm_atom_index(list(i))
752 368 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
753 : END DO
754 : END DO
755 :
756 : ! derivatives from qm charges
757 208 : Forces_QM = 0.0_dp
758 16 : IF (SIZE(matrix_p) == 2) THEN
759 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
760 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
761 : END IF
762 : !
763 16 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
764 64 : DO WHILE (dbcsr_iterator_blocks_left(iter))
765 48 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock, blk)
766 : !
767 48 : IF (iatom == jatom) CYCLE
768 : !
769 24 : gmij = -0.5_dp*(qpot(iatom) + qpot(jatom))
770 24 : NULLIFY (pblock)
771 : CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
772 24 : row=iatom, col=jatom, block=pblock, found=found)
773 24 : CPASSERT(found)
774 112 : DO i = 1, 3
775 72 : NULLIFY (dsblock)
776 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
777 72 : row=iatom, col=jatom, block=dsblock, found=found)
778 72 : CPASSERT(found)
779 648 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
780 72 : Forces_QM(i, iatom) = Forces_QM(i, iatom) + fi
781 192 : Forces_QM(i, jatom) = Forces_QM(i, jatom) - fi
782 : END DO
783 : END DO
784 16 : CALL dbcsr_iterator_stop(iter)
785 : !
786 16 : IF (SIZE(matrix_p) == 2) THEN
787 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
788 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
789 : END IF
790 : !
791 : ! Transfer QM gradients to the QM particles..
792 400 : CALL para_env%sum(Forces_QM)
793 48 : DO ikind = 1, nkind
794 32 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
795 96 : DO i = 1, SIZE(list)
796 48 : iqm = list(i)
797 48 : iatom = qmmm_env%qm_atom_index(iqm)
798 368 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
799 : END DO
800 : END DO
801 : !
802 16 : 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 16 : DEALLOCATE (Forces_QM)
807 16 : DEALLOCATE (qpot)
808 :
809 16 : CALL dbcsr_deallocate_matrix_set(matrix_s)
810 :
811 : END IF
812 448 : CALL timestop(handle)
813 :
814 448 : 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 1116 : 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 1116 : 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 1116 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges_mm, mcharge, qpot
850 1116 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges, dsblock, Forces_MM, Forces_QM, &
851 1116 : pblock, sblock
852 1116 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
853 : TYPE(dbcsr_iterator_type) :: iter
854 1116 : 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 1116 : POINTER :: sab_nl
861 1116 : 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 1116 : 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 1116 : CALL timeset(routineN, handle)
874 1116 : IF (calc_force) THEN
875 36 : 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 36 : natom=number_qm_atoms)
884 36 : dftb_control => dft_control%qs_control%dftb_control
885 36 : xtb_control => dft_control%qs_control%xtb_control
886 :
887 36 : IF (dft_control%qs_control%dftb) THEN
888 : do_dftb = .TRUE.
889 : do_xtb = .FALSE.
890 18 : ELSEIF (dft_control%qs_control%xtb) THEN
891 : do_dftb = .FALSE.
892 : do_xtb = .TRUE.
893 : ELSE
894 0 : CPABORT("TB method unknown")
895 : END IF
896 :
897 36 : NULLIFY (matrix_s)
898 36 : IF (do_dftb) THEN
899 18 : CALL build_dftb_overlap(qs_env, 1, matrix_s)
900 18 : ELSEIF (do_xtb) THEN
901 18 : 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 18 : basis_type_a='ORB', basis_type_b='ORB', sab_nl=sab_nl)
904 : END IF
905 36 : CALL qs_rho_get(rho, rho_ao=matrix_p)
906 :
907 36 : nspins = dft_control%nspins
908 36 : nkind = SIZE(atomic_kind_set)
909 : ! Mulliken charges
910 144 : ALLOCATE (charges(number_qm_atoms, nspins))
911 : !
912 36 : CALL mulliken_charges(matrix_p, matrix_s(1)%matrix, para_env, charges)
913 : !
914 108 : ALLOCATE (mcharge(number_qm_atoms))
915 108 : DO ikind = 1, nkind
916 72 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
917 72 : IF (do_dftb) THEN
918 36 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
919 36 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
920 36 : ELSEIF (do_xtb) THEN
921 36 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
922 36 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
923 : END IF
924 288 : DO iatom = 1, natom
925 108 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
926 288 : mcharge(atom_a) = zeff - SUM(charges(atom_a, 1:nspins))
927 : END DO
928 : END DO
929 36 : DEALLOCATE (charges)
930 :
931 108 : ALLOCATE (qpot(number_qm_atoms))
932 144 : qpot = 0.0_dp
933 108 : ALLOCATE (Forces_QM(3, number_qm_atoms))
934 468 : Forces_QM = 0.0_dp
935 :
936 : ! Create Ewald environments
937 36 : poisson_section => section_vals_get_subs_vals(qs_env%input, "MM%POISSON")
938 468 : ALLOCATE (ewald_env)
939 36 : CALL ewald_env_create(ewald_env, para_env)
940 36 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
941 36 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
942 36 : CALL read_ewald_section(ewald_env, ewald_section)
943 36 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
944 36 : ALLOCATE (ewald_pw)
945 36 : CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=print_section)
946 :
947 36 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, do_ipol=do_ipol)
948 36 : IF (do_multipoles) CPABORT("No multipole force fields allowed in DFTB QM/MM")
949 36 : IF (do_ipol /= do_fist_pol_none) CPABORT("No polarizable force fields allowed in DFTB QM/MM")
950 :
951 0 : SELECT CASE (ewald_type)
952 : CASE (do_ewald_pme)
953 0 : CPABORT("PME Ewald type not implemented for DFTB/QMMM")
954 : CASE (do_ewald_ewald, do_ewald_spme)
955 60 : DO ipot = 1, SIZE(qmmm_env%Potentials)
956 36 : Pot => qmmm_env%Potentials(ipot)%Pot
957 36 : nmm = SIZE(Pot%mm_atom_index)
958 : ! get a 'clean' mm particle set
959 36 : NULLIFY (atoms_mm)
960 36 : CALL allocate_particle_set(atoms_mm, nmm)
961 108 : ALLOCATE (charges_mm(nmm))
962 144 : DO Imp = 1, nmm
963 108 : Imm = Pot%mm_atom_index(Imp)
964 108 : IndMM = qmmm_env%mm_atom_index(Imm)
965 756 : atoms_mm(Imp)%r = particles_mm(IndMM)%r
966 108 : atoms_mm(Imp)%atomic_kind => particles_mm(IndMM)%atomic_kind
967 144 : charges_mm(Imp) = qmmm_env%mm_atom_chrg(Imm)
968 : END DO
969 : ! force array for mm atoms
970 108 : ALLOCATE (Forces_MM(3, nmm))
971 468 : Forces_MM = 0.0_dp
972 36 : IF (ewald_type == do_ewald_ewald) THEN
973 0 : CPABORT("Ewald not implemented for DFTB/QMMM")
974 36 : 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 36 : particles_qm, qpot)
978 : ! forces QM
979 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
980 36 : particles_qm, mcharge, Forces_QM)
981 : ! forces MM
982 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, particles_qm, mcharge, &
983 36 : atoms_mm, charges_mm, Forces_MM)
984 : END IF
985 36 : CALL deallocate_particle_set(atoms_mm)
986 36 : DEALLOCATE (charges_mm)
987 : ! transfer MM forces
988 900 : CALL para_env%sum(Forces_MM)
989 144 : DO Imp = 1, nmm
990 108 : Imm = Pot%mm_atom_index(Imp)
991 792 : Forces(:, Imm) = Forces(:, Imm) - Forces_MM(:, Imp)
992 : END DO
993 60 : DEALLOCATE (Forces_MM)
994 : END DO
995 :
996 24 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
997 0 : DO ipot = 1, SIZE(qmmm_env%added_charges%Potentials)
998 0 : Pot => qmmm_env%added_charges%Potentials(ipot)%Pot
999 0 : nmm = SIZE(Pot%mm_atom_index)
1000 : ! get a 'clean' mm particle set
1001 0 : NULLIFY (atoms_mm)
1002 0 : CALL allocate_particle_set(atoms_mm, nmm)
1003 0 : ALLOCATE (charges_mm(nmm))
1004 0 : DO Imp = 1, nmm
1005 0 : Imm = Pot%mm_atom_index(Imp)
1006 0 : IndMM = qmmm_env%added_charges%mm_atom_index(Imm)
1007 0 : atoms_mm(Imp)%r = qmmm_env%added_charges%added_particles(IndMM)%r
1008 0 : atoms_mm(Imp)%atomic_kind => qmmm_env%added_charges%added_particles(IndMM)%atomic_kind
1009 0 : charges_mm(Imp) = qmmm_env%added_charges%mm_atom_chrg(Imm)
1010 : END DO
1011 : ! force array for mm atoms
1012 0 : ALLOCATE (Forces_MM(3, nmm))
1013 0 : Forces_MM = 0.0_dp
1014 0 : IF (ewald_type == do_ewald_ewald) THEN
1015 0 : CPABORT("Ewald not implemented for DFTB/QMMM")
1016 0 : 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 0 : charges_mm, particles_qm, qpot)
1020 : ! forces QM
1021 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, atoms_mm, charges_mm, &
1022 0 : particles_qm, mcharge, Forces_QM)
1023 : ! forces MM
1024 : CALL spme_forces(ewald_env, ewald_pw, mm_cell, particles_qm, mcharge, &
1025 0 : atoms_mm, charges_mm, Forces_MM)
1026 : END IF
1027 0 : CALL deallocate_particle_set(atoms_mm)
1028 : ! transfer MM forces
1029 0 : CALL para_env%sum(Forces_MM)
1030 0 : DO Imp = 1, nmm
1031 0 : Imm = Pot%mm_atom_index(Imp)
1032 0 : Forces_added_charges(:, Imm) = Forces_added_charges(:, Imm) - Forces_MM(:, Imp)
1033 : END DO
1034 24 : DEALLOCATE (Forces_MM)
1035 : END DO
1036 : END IF
1037 168 : CALL para_env%sum(qpot)
1038 600 : 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 24 : CALL ewald_env_get(ewald_env, rcut=rcutoff(1), alpha=alpha)
1043 24 : rcutoff(2) = 0.025_dp*rcutoff(1)
1044 24 : rcutoff(1) = 2.0_dp*rcutoff(1)
1045 24 : nkind = SIZE(atomic_kind_set)
1046 24 : iqm = 0
1047 72 : DO ikind = 1, nkind
1048 48 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1049 48 : IF (do_dftb) THEN
1050 24 : NULLIFY (dftb_kind)
1051 24 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1052 : CALL get_dftb_atom_param(dftb_kind, &
1053 24 : defined=defined, eta=eta_a, natorb=natorb)
1054 : ! use mm charge smearing for non-scc cases
1055 24 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
1056 24 : IF (.NOT. defined .OR. natorb < 1) CYCLE
1057 24 : ELSEIF (do_xtb) THEN
1058 24 : eta_a(0) = eta_mm
1059 : END IF
1060 144 : DO i = 1, SIZE(list)
1061 72 : iatom = list(i)
1062 72 : 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 72 : 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 72 : 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 72 : 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 72 : rcutoff, particles_qm)
1077 : ! Possibly added charges
1078 120 : 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 0 : 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 0 : 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 0 : 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 0 : 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 36 : DO ikind = 1, nkind
1110 24 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1111 24 : IF (do_dftb) THEN
1112 12 : NULLIFY (dftb_kind)
1113 12 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1114 : CALL get_dftb_atom_param(dftb_kind, &
1115 12 : defined=defined, eta=eta_a, natorb=natorb)
1116 : ! use mm charge smearing for non-scc cases
1117 12 : IF (.NOT. dftb_control%self_consistent) eta_a(0) = eta_mm
1118 12 : IF (.NOT. defined .OR. natorb < 1) CYCLE
1119 12 : ELSEIF (do_xtb) THEN
1120 12 : eta_a(0) = eta_mm
1121 : END IF
1122 72 : DO i = 1, SIZE(list)
1123 36 : iatom = list(i)
1124 36 : 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 36 : 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 36 : qmmm_env%spherical_cutoff, particles_qm)
1132 : ! Possibly added charges
1133 60 : 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 0 : 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 0 : Forces_QM(:, iqm), qmmm_env%spherical_cutoff, particles_qm)
1146 : END IF
1147 : END DO
1148 : END DO
1149 : CASE DEFAULT
1150 36 : CPABORT("Unknown Ewald type!")
1151 : END SELECT
1152 :
1153 : ! Transfer QM gradients to the QM particles..
1154 36 : iqm = 0
1155 108 : DO ikind = 1, nkind
1156 72 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1157 72 : IF (do_dftb) THEN
1158 36 : NULLIFY (dftb_kind)
1159 36 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
1160 36 : CALL get_dftb_atom_param(dftb_kind, defined=defined, natorb=natorb)
1161 36 : IF (.NOT. defined .OR. natorb < 1) CYCLE
1162 : ELSEIF (do_xtb) THEN
1163 : !
1164 : END IF
1165 216 : DO i = 1, SIZE(list)
1166 108 : iqm = iqm + 1
1167 108 : iatom = qmmm_env%qm_atom_index(list(i))
1168 828 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
1169 : END DO
1170 : END DO
1171 :
1172 : ! derivatives from qm charges
1173 468 : Forces_QM = 0.0_dp
1174 36 : IF (SIZE(matrix_p) == 2) THEN
1175 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
1176 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1177 : END IF
1178 : !
1179 36 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
1180 144 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1181 108 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, sblock, blk)
1182 : !
1183 108 : IF (iatom == jatom) CYCLE
1184 : !
1185 54 : gmij = -0.5_dp*(qpot(iatom) + qpot(jatom))
1186 54 : NULLIFY (pblock)
1187 : CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
1188 54 : row=iatom, col=jatom, block=pblock, found=found)
1189 54 : CPASSERT(found)
1190 252 : DO i = 1, 3
1191 162 : NULLIFY (dsblock)
1192 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
1193 162 : row=iatom, col=jatom, block=dsblock, found=found)
1194 162 : CPASSERT(found)
1195 1458 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
1196 162 : Forces_QM(i, iatom) = Forces_QM(i, iatom) + fi
1197 432 : Forces_QM(i, jatom) = Forces_QM(i, jatom) - fi
1198 : END DO
1199 : END DO
1200 36 : CALL dbcsr_iterator_stop(iter)
1201 : !
1202 36 : IF (SIZE(matrix_p) == 2) THEN
1203 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
1204 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1205 : END IF
1206 : !
1207 : ! Transfer QM gradients to the QM particles..
1208 900 : CALL para_env%sum(Forces_QM)
1209 108 : DO ikind = 1, nkind
1210 72 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
1211 216 : DO i = 1, SIZE(list)
1212 108 : iqm = list(i)
1213 108 : iatom = qmmm_env%qm_atom_index(iqm)
1214 828 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
1215 : END DO
1216 : END DO
1217 : !
1218 36 : 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 36 : DEALLOCATE (Forces_QM)
1223 36 : DEALLOCATE (qpot)
1224 :
1225 : ! Release Ewald environment
1226 36 : CALL ewald_env_release(ewald_env)
1227 36 : DEALLOCATE (ewald_env)
1228 36 : CALL ewald_pw_release(ewald_pw)
1229 36 : DEALLOCATE (ewald_pw)
1230 :
1231 72 : CALL dbcsr_deallocate_matrix_set(matrix_s)
1232 :
1233 : END IF
1234 :
1235 1116 : CALL timestop(handle)
1236 :
1237 1116 : 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 6936 : 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 6936 : CALL timeset(routineN, handle)
1279 : ! Loop Over MM atoms
1280 : ! Loop over Pot stores atoms with the same charge
1281 20592 : MainLoopPot: DO Ipot = 1, SIZE(Potentials)
1282 13656 : Pot => Potentials(Ipot)%Pot
1283 : ! Loop over atoms belonging to this type
1284 61560 : LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index)
1285 40968 : Imm = Pot%mm_atom_index(Imp)
1286 40968 : IndMM = mm_atom_index(Imm)
1287 163872 : r_pbc = pbc(particles_mm(IndMM)%r - particles_qm(IndQM)%r, mm_cell)
1288 40968 : rt1 = r_pbc(1)
1289 40968 : rt2 = r_pbc(2)
1290 40968 : rt3 = r_pbc(3)
1291 163872 : rij = (/rt1, rt2, rt3/)
1292 163872 : dr = SQRT(SUM(rij**2))
1293 40968 : qeff = mm_charges(Imm)
1294 : ! Computes the screening factor for the spherical cutoff (if defined)
1295 40968 : IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1296 24624 : CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
1297 24624 : qeff = qeff*sph_chrg_factor
1298 : END IF
1299 40968 : IF (ABS(qeff) <= qsmall) CYCLE
1300 54624 : IF (dr > rtiny) THEN
1301 40968 : IF (pot_type == 0) THEN
1302 16344 : sr = gamma_rab_sr(dr, qm_alpha, eta_mm, 0.0_dp)
1303 16344 : qpot = qpot + qeff*(1.0_dp/dr - sr)
1304 24624 : ELSE IF (pot_type == 1) THEN
1305 12312 : sr = gamma_rab_sr(dr, qm_alpha, eta_mm, 0.0_dp)
1306 12312 : qpot = qpot - qeff*sr
1307 12312 : ELSE IF (pot_type == 2) THEN
1308 12312 : sr = erfc(qm_alpha*dr)/dr
1309 12312 : qpot = qpot + qeff*sr
1310 : ELSE
1311 0 : CPABORT("")
1312 : END IF
1313 : END IF
1314 : END DO LoopMM
1315 : END DO MainLoopPot
1316 6936 : CALL timestop(handle)
1317 6936 : 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 456 : SUBROUTINE build_mm_dpot(qcharge, pot_type, qm_alpha, potentials, &
1336 : particles_mm, mm_charges, mm_atom_index, mm_cell, IndQM, &
1337 228 : 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 228 : CALL timeset(routineN, handle)
1363 : ! Loop Over MM atoms
1364 : ! Loop over Pot stores atoms with the same charge
1365 576 : MainLoopPot: DO Ipot = 1, SIZE(Potentials)
1366 348 : Pot => Potentials(Ipot)%Pot
1367 : ! Loop over atoms belonging to this type
1368 1620 : LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index)
1369 1044 : Imm = Pot%mm_atom_index(Imp)
1370 1044 : IndMM = mm_atom_index(Imm)
1371 4176 : r_pbc = pbc(particles_mm(IndMM)%r - particles_qm(IndQM)%r, mm_cell)
1372 1044 : rt1 = r_pbc(1)
1373 1044 : rt2 = r_pbc(2)
1374 1044 : rt3 = r_pbc(3)
1375 4176 : rij = (/rt1, rt2, rt3/)
1376 4176 : dr = SQRT(SUM(rij**2))
1377 1044 : 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 1044 : IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1381 648 : CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
1382 648 : qeff = qeff*sph_chrg_factor
1383 : END IF
1384 1044 : IF (ABS(qeff) <= qsmall) CYCLE
1385 1044 : IF (dr > rtiny) THEN
1386 1044 : drp = dr + ddrmm
1387 1044 : drm = dr - ddrmm
1388 1044 : IF (pot_type == 0) THEN
1389 : dsr = 0.5_dp*(gamma_rab_sr(drp, qm_alpha, eta_mm, 0.0_dp) - &
1390 396 : gamma_rab_sr(drm, qm_alpha, eta_mm, 0.0_dp))/ddrmm
1391 396 : fsr = qeff*qcharge*(-1.0_dp/(dr*dr) - dsr)
1392 648 : ELSE IF (pot_type == 1) THEN
1393 : dsr = 0.5_dp*(gamma_rab_sr(drp, qm_alpha, eta_mm, 0.0_dp) - &
1394 324 : gamma_rab_sr(drm, qm_alpha, eta_mm, 0.0_dp))/ddrmm
1395 324 : fsr = -qeff*qcharge*dsr
1396 324 : ELSE IF (pot_type == 2) THEN
1397 324 : dsr = 0.5_dp*(erfc(qm_alpha*drp)/drp - erfc(qm_alpha*drm)/drm)/ddrmm
1398 324 : fsr = qeff*qcharge*dsr
1399 : ELSE
1400 0 : CPABORT("")
1401 : END IF
1402 4176 : force_ab = -fsr*rij/dr
1403 : ELSE
1404 0 : force_ab = 0.0_dp
1405 : END IF
1406 : ! The array of QM forces are really the forces
1407 4176 : forces_qm(:) = forces_qm(:) - force_ab
1408 : ! The one of MM atoms are instead gradients
1409 4524 : forces(:, Imm) = forces(:, Imm) - force_ab
1410 : END DO LoopMM
1411 : END DO MainLoopPot
1412 :
1413 228 : CALL timestop(handle)
1414 :
1415 228 : END SUBROUTINE build_mm_dpot
1416 :
1417 : END MODULE qmmm_tb_methods
1418 :
|
