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 Calculation of forces for Coulomb contributions in response xTB
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE xtb_ehess_force
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind_set
15 : USE cell_types, ONLY: cell_type,&
16 : get_cell,&
17 : pbc
18 : USE cp_control_types, ONLY: dft_control_type,&
19 : xtb_control_type
20 : USE cp_dbcsr_api, ONLY: &
21 : dbcsr_add, dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
22 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_type
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_get_default_unit_nr,&
25 : cp_logger_type
26 : USE distribution_1d_types, ONLY: distribution_1d_type
27 : USE ewald_environment_types, ONLY: ewald_env_get,&
28 : ewald_environment_type
29 : USE ewald_methods_tb, ONLY: tb_ewald_overlap,&
30 : tb_spme_zforce
31 : USE ewald_pw_types, ONLY: ewald_pw_type
32 : USE kinds, ONLY: dp
33 : USE mathconstants, ONLY: oorootpi,&
34 : pi
35 : USE message_passing, ONLY: mp_para_env_type
36 : USE particle_types, ONLY: particle_type
37 : USE pw_poisson_types, ONLY: do_ewald_ewald,&
38 : do_ewald_none,&
39 : do_ewald_pme,&
40 : do_ewald_spme
41 : USE qs_energy_types, ONLY: qs_energy_type
42 : USE qs_environment_types, ONLY: get_qs_env,&
43 : qs_environment_type
44 : USE qs_force_types, ONLY: qs_force_type
45 : USE qs_kind_types, ONLY: get_qs_kind,&
46 : qs_kind_type
47 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
48 : neighbor_list_iterate,&
49 : neighbor_list_iterator_create,&
50 : neighbor_list_iterator_p_type,&
51 : neighbor_list_iterator_release,&
52 : neighbor_list_set_p_type
53 : USE qs_rho_types, ONLY: qs_rho_type
54 : USE virial_types, ONLY: virial_type
55 : USE xtb_coulomb, ONLY: dgamma_rab_sr,&
56 : gamma_rab_sr
57 : USE xtb_types, ONLY: get_xtb_atom_param,&
58 : xtb_atom_type
59 : #include "./base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 :
63 : PRIVATE
64 :
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_ehess_force'
66 :
67 : PUBLIC :: calc_xtb_ehess_force
68 :
69 : ! **************************************************************************************************
70 :
71 : CONTAINS
72 :
73 : ! **************************************************************************************************
74 : !> \brief ...
75 : !> \param qs_env ...
76 : !> \param matrix_p0 ...
77 : !> \param matrix_p1 ...
78 : !> \param charges0 ...
79 : !> \param mcharge0 ...
80 : !> \param charges1 ...
81 : !> \param mcharge1 ...
82 : !> \param debug_forces ...
83 : ! **************************************************************************************************
84 14 : SUBROUTINE calc_xtb_ehess_force(qs_env, matrix_p0, matrix_p1, charges0, mcharge0, &
85 14 : charges1, mcharge1, debug_forces)
86 :
87 : TYPE(qs_environment_type), POINTER :: qs_env
88 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p0, matrix_p1
89 : REAL(KIND=dp), DIMENSION(:, :), INTENT(in) :: charges0
90 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: mcharge0
91 : REAL(KIND=dp), DIMENSION(:, :), INTENT(in) :: charges1
92 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: mcharge1
93 : LOGICAL, INTENT(IN) :: debug_forces
94 :
95 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_xtb_ehess_force'
96 :
97 : INTEGER :: atom_i, atom_j, blk, ewald_type, handle, i, ia, iatom, icol, ikind, iounit, irow, &
98 : j, jatom, jkind, la, lb, lmaxa, lmaxb, natom, natorb_a, natorb_b, ni, nimg, nj, nkind, &
99 : nmat, za, zb
100 14 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
101 : INTEGER, DIMENSION(25) :: laoa, laob
102 : INTEGER, DIMENSION(3) :: cellind, periodic
103 : LOGICAL :: calculate_forces, defined, do_ewald, &
104 : found, just_energy, use_virial
105 : REAL(KIND=dp) :: alpha, deth, dr, etaa, etab, fi, gmij0, &
106 : gmij1, kg, rcut, rcuta, rcutb
107 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xgamma
108 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gammab, gcij0, gcij1, gmcharge0, &
109 14 : gmcharge1
110 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gchrg0, gchrg1
111 : REAL(KIND=dp), DIMENSION(3) :: fij, fodeb, rij
112 : REAL(KIND=dp), DIMENSION(5) :: kappaa, kappab
113 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, pblock0, pblock1, sblock
114 14 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
115 : TYPE(cell_type), POINTER :: cell
116 : TYPE(cp_logger_type), POINTER :: logger
117 : TYPE(dbcsr_iterator_type) :: iter
118 14 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
119 : TYPE(dft_control_type), POINTER :: dft_control
120 : TYPE(distribution_1d_type), POINTER :: local_particles
121 : TYPE(ewald_environment_type), POINTER :: ewald_env
122 : TYPE(ewald_pw_type), POINTER :: ewald_pw
123 : TYPE(mp_para_env_type), POINTER :: para_env
124 : TYPE(neighbor_list_iterator_p_type), &
125 14 : DIMENSION(:), POINTER :: nl_iterator
126 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
127 14 : POINTER :: n_list
128 14 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
129 : TYPE(qs_energy_type), POINTER :: energy
130 14 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
131 14 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
132 : TYPE(qs_rho_type), POINTER :: rho
133 : TYPE(virial_type), POINTER :: virial
134 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind
135 : TYPE(xtb_control_type), POINTER :: xtb_control
136 :
137 14 : CALL timeset(routineN, handle)
138 :
139 14 : logger => cp_get_default_logger()
140 14 : IF (logger%para_env%is_source()) THEN
141 7 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
142 : ELSE
143 : iounit = -1
144 : END IF
145 :
146 14 : CPASSERT(ASSOCIATED(matrix_p1))
147 :
148 : CALL get_qs_env(qs_env, &
149 : qs_kind_set=qs_kind_set, &
150 : particle_set=particle_set, &
151 : cell=cell, &
152 : rho=rho, &
153 : energy=energy, &
154 : virial=virial, &
155 14 : dft_control=dft_control)
156 :
157 14 : xtb_control => dft_control%qs_control%xtb_control
158 :
159 14 : calculate_forces = .TRUE.
160 14 : just_energy = .FALSE.
161 14 : use_virial = .FALSE.
162 14 : nmat = 4
163 14 : nimg = dft_control%nimages
164 14 : IF (nimg > 1) THEN
165 0 : CPABORT('xTB-sTDA forces for k-points not available')
166 : END IF
167 :
168 14 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
169 56 : ALLOCATE (gchrg0(natom, 5, nmat))
170 5030 : gchrg0 = 0._dp
171 42 : ALLOCATE (gmcharge0(natom, nmat))
172 1006 : gmcharge0 = 0._dp
173 28 : ALLOCATE (gchrg1(natom, 5, nmat))
174 5030 : gchrg1 = 0._dp
175 28 : ALLOCATE (gmcharge1(natom, nmat))
176 1006 : gmcharge1 = 0._dp
177 :
178 : ! short range contribution (gamma)
179 : ! loop over all atom pairs (sab_xtbe)
180 14 : kg = xtb_control%kg
181 14 : NULLIFY (n_list)
182 14 : IF (xtb_control%old_coulomb_damping) THEN
183 0 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
184 : ELSE
185 14 : CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
186 : END IF
187 14 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
188 24708 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
189 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
190 24694 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
191 24694 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
192 24694 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
193 24694 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
194 24694 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
195 24694 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
196 24694 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
197 : ! atomic parameters
198 24694 : CALL get_xtb_atom_param(xtb_atom_a, eta=etaa, lmax=lmaxa, kappa=kappaa, rcut=rcuta)
199 24694 : CALL get_xtb_atom_param(xtb_atom_b, eta=etab, lmax=lmaxb, kappa=kappab, rcut=rcutb)
200 : ! gamma matrix
201 24694 : ni = lmaxa + 1
202 24694 : nj = lmaxb + 1
203 98776 : ALLOCATE (gammab(ni, nj))
204 24694 : rcut = rcuta + rcutb
205 98776 : dr = SQRT(SUM(rij(:)**2))
206 24694 : CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
207 250847 : gchrg0(iatom, 1:ni, 1) = gchrg0(iatom, 1:ni, 1) + MATMUL(gammab, charges0(jatom, 1:nj))
208 250847 : gchrg1(iatom, 1:ni, 1) = gchrg1(iatom, 1:ni, 1) + MATMUL(gammab, charges1(jatom, 1:nj))
209 24694 : IF (iatom /= jatom) THEN
210 285738 : gchrg0(jatom, 1:nj, 1) = gchrg0(jatom, 1:nj, 1) + MATMUL(charges0(iatom, 1:ni), gammab)
211 285738 : gchrg1(jatom, 1:nj, 1) = gchrg1(jatom, 1:nj, 1) + MATMUL(charges1(iatom, 1:ni), gammab)
212 : END IF
213 24694 : IF (dr > 1.e-6_dp) THEN
214 24577 : CALL dgamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
215 98308 : DO i = 1, 3
216 : gchrg0(iatom, 1:ni, i + 1) = gchrg0(iatom, 1:ni, i + 1) &
217 748626 : + MATMUL(gammab, charges0(jatom, 1:nj))*rij(i)/dr
218 : gchrg1(iatom, 1:ni, i + 1) = gchrg1(iatom, 1:ni, i + 1) &
219 748626 : + MATMUL(gammab, charges1(jatom, 1:nj))*rij(i)/dr
220 98308 : IF (iatom /= jatom) THEN
221 : gchrg0(jatom, 1:nj, i + 1) = gchrg0(jatom, 1:nj, i + 1) &
222 857214 : - MATMUL(charges0(iatom, 1:ni), gammab)*rij(i)/dr
223 : gchrg1(jatom, 1:nj, i + 1) = gchrg1(jatom, 1:nj, i + 1) &
224 857214 : - MATMUL(charges1(iatom, 1:ni), gammab)*rij(i)/dr
225 : END IF
226 : END DO
227 : END IF
228 98776 : DEALLOCATE (gammab)
229 : END DO
230 14 : CALL neighbor_list_iterator_release(nl_iterator)
231 :
232 : ! 1/R contribution
233 :
234 14 : IF (xtb_control%coulomb_lr) THEN
235 14 : do_ewald = xtb_control%do_ewald
236 14 : IF (do_ewald) THEN
237 : ! Ewald sum
238 6 : NULLIFY (ewald_env, ewald_pw)
239 : CALL get_qs_env(qs_env=qs_env, &
240 6 : ewald_env=ewald_env, ewald_pw=ewald_pw)
241 6 : CALL get_cell(cell=cell, periodic=periodic, deth=deth)
242 6 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
243 6 : CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
244 6 : CALL tb_ewald_overlap(gmcharge0, mcharge0, alpha, n_list, virial, use_virial)
245 6 : CALL tb_ewald_overlap(gmcharge1, mcharge1, alpha, n_list, virial, use_virial)
246 0 : SELECT CASE (ewald_type)
247 : CASE DEFAULT
248 0 : CPABORT("Invalid Ewald type")
249 : CASE (do_ewald_none)
250 0 : CPABORT("Not allowed with DFTB")
251 : CASE (do_ewald_ewald)
252 0 : CPABORT("Standard Ewald not implemented in DFTB")
253 : CASE (do_ewald_pme)
254 0 : CPABORT("PME not implemented in DFTB")
255 : CASE (do_ewald_spme)
256 6 : CALL tb_spme_zforce(ewald_env, ewald_pw, particle_set, cell, gmcharge0, mcharge0)
257 12 : CALL tb_spme_zforce(ewald_env, ewald_pw, particle_set, cell, gmcharge1, mcharge1)
258 : END SELECT
259 : ELSE
260 : ! direct sum
261 8 : CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
262 28 : DO ikind = 1, SIZE(local_particles%n_el)
263 42 : DO ia = 1, local_particles%n_el(ikind)
264 14 : iatom = local_particles%list(ikind)%array(ia)
265 52 : DO jatom = 1, iatom - 1
266 72 : rij = particle_set(iatom)%r - particle_set(jatom)%r
267 72 : rij = pbc(rij, cell)
268 72 : dr = SQRT(SUM(rij(:)**2))
269 32 : IF (dr > 1.e-6_dp) THEN
270 18 : gmcharge0(iatom, 1) = gmcharge0(iatom, 1) + mcharge0(jatom)/dr
271 18 : gmcharge0(jatom, 1) = gmcharge0(jatom, 1) + mcharge0(iatom)/dr
272 18 : gmcharge1(iatom, 1) = gmcharge1(iatom, 1) + mcharge1(jatom)/dr
273 18 : gmcharge1(jatom, 1) = gmcharge1(jatom, 1) + mcharge1(iatom)/dr
274 72 : DO i = 2, nmat
275 54 : gmcharge0(iatom, i) = gmcharge0(iatom, i) + rij(i - 1)*mcharge0(jatom)/dr**3
276 54 : gmcharge0(jatom, i) = gmcharge0(jatom, i) - rij(i - 1)*mcharge0(iatom)/dr**3
277 54 : gmcharge1(iatom, i) = gmcharge1(iatom, i) + rij(i - 1)*mcharge1(jatom)/dr**3
278 72 : gmcharge1(jatom, i) = gmcharge1(jatom, i) - rij(i - 1)*mcharge1(iatom)/dr**3
279 : END DO
280 : END IF
281 : END DO
282 : END DO
283 : END DO
284 : CPASSERT(.NOT. use_virial)
285 : END IF
286 : END IF
287 :
288 : ! global sum of gamma*p arrays
289 : CALL get_qs_env(qs_env=qs_env, &
290 : atomic_kind_set=atomic_kind_set, &
291 14 : force=force, para_env=para_env)
292 14 : CALL para_env%sum(gmcharge0(:, 1))
293 14 : CALL para_env%sum(gchrg0(:, :, 1))
294 14 : CALL para_env%sum(gmcharge1(:, 1))
295 14 : CALL para_env%sum(gchrg1(:, :, 1))
296 :
297 14 : IF (xtb_control%coulomb_lr) THEN
298 14 : IF (do_ewald) THEN
299 : ! add self charge interaction and background charge contribution
300 212 : gmcharge0(:, 1) = gmcharge0(:, 1) - 2._dp*alpha*oorootpi*mcharge0(:)
301 12 : IF (ANY(periodic(:) == 1)) THEN
302 202 : gmcharge0(:, 1) = gmcharge0(:, 1) - pi/alpha**2/deth
303 : END IF
304 212 : gmcharge1(:, 1) = gmcharge1(:, 1) - 2._dp*alpha*oorootpi*mcharge1(:)
305 12 : IF (ANY(periodic(:) == 1)) THEN
306 202 : gmcharge1(:, 1) = gmcharge1(:, 1) - pi/alpha**2/deth
307 : END IF
308 : END IF
309 : END IF
310 :
311 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
312 : kind_of=kind_of, &
313 14 : atom_of_kind=atom_of_kind)
314 :
315 14 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
316 248 : DO iatom = 1, natom
317 234 : ikind = kind_of(iatom)
318 234 : atom_i = atom_of_kind(iatom)
319 234 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
320 234 : CALL get_xtb_atom_param(xtb_kind, lmax=ni)
321 234 : ni = ni + 1
322 : ! short range
323 234 : fij = 0.0_dp
324 936 : DO i = 1, 3
325 : fij(i) = SUM(charges0(iatom, 1:ni)*gchrg1(iatom, 1:ni, i + 1)) + &
326 2832 : SUM(charges1(iatom, 1:ni)*gchrg0(iatom, 1:ni, i + 1))
327 : END DO
328 234 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
329 234 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
330 234 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
331 : ! long range
332 234 : fij = 0.0_dp
333 936 : DO i = 1, 3
334 : fij(i) = gmcharge1(iatom, i + 1)*mcharge0(iatom) + &
335 936 : gmcharge0(iatom, i + 1)*mcharge1(iatom)
336 : END DO
337 234 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
338 234 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
339 482 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
340 : END DO
341 14 : IF (debug_forces) THEN
342 0 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
343 0 : CALL para_env%sum(fodeb)
344 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dH[Pz] ", fodeb
345 : END IF
346 :
347 14 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
348 :
349 14 : IF (SIZE(matrix_p0) == 2) THEN
350 : CALL dbcsr_add(matrix_p0(1)%matrix, matrix_p0(2)%matrix, &
351 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
352 : CALL dbcsr_add(matrix_p1(1)%matrix, matrix_p1(2)%matrix, &
353 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
354 : END IF
355 :
356 : ! no k-points; all matrices have been transformed to periodic bsf
357 14 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
358 14 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
359 4695 : DO WHILE (dbcsr_iterator_blocks_left(iter))
360 4681 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
361 4681 : ikind = kind_of(irow)
362 4681 : jkind = kind_of(icol)
363 :
364 : ! atomic parameters
365 4681 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
366 4681 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
367 4681 : CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
368 4681 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
369 :
370 4681 : ni = SIZE(sblock, 1)
371 4681 : nj = SIZE(sblock, 2)
372 18724 : ALLOCATE (gcij0(ni, nj))
373 14043 : ALLOCATE (gcij1(ni, nj))
374 19209 : DO i = 1, ni
375 52401 : DO j = 1, nj
376 33192 : la = laoa(i) + 1
377 33192 : lb = laob(j) + 1
378 33192 : gcij0(i, j) = 0.5_dp*(gchrg0(irow, la, 1) + gchrg0(icol, lb, 1))
379 47720 : gcij1(i, j) = 0.5_dp*(gchrg1(irow, la, 1) + gchrg1(icol, lb, 1))
380 : END DO
381 : END DO
382 4681 : gmij0 = 0.5_dp*(gmcharge0(irow, 1) + gmcharge0(icol, 1))
383 4681 : gmij1 = 0.5_dp*(gmcharge1(irow, 1) + gmcharge1(icol, 1))
384 4681 : atom_i = atom_of_kind(irow)
385 4681 : atom_j = atom_of_kind(icol)
386 4681 : NULLIFY (pblock0)
387 : CALL dbcsr_get_block_p(matrix=matrix_p0(1)%matrix, &
388 4681 : row=irow, col=icol, block=pblock0, found=found)
389 4681 : CPASSERT(found)
390 4681 : NULLIFY (pblock1)
391 : CALL dbcsr_get_block_p(matrix=matrix_p1(1)%matrix, &
392 4681 : row=irow, col=icol, block=pblock1, found=found)
393 4681 : CPASSERT(found)
394 18724 : DO i = 1, 3
395 14043 : NULLIFY (dsblock)
396 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
397 14043 : row=irow, col=icol, block=dsblock, found=found)
398 14043 : CPASSERT(found)
399 : ! short range
400 275571 : fi = -2.0_dp*SUM(pblock0*dsblock*gcij1) - 2.0_dp*SUM(pblock1*dsblock*gcij0)
401 14043 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
402 14043 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
403 : ! long range
404 275571 : fi = -2.0_dp*gmij1*SUM(pblock0*dsblock) - 2.0_dp*gmij0*SUM(pblock1*dsblock)
405 14043 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
406 32767 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
407 : END DO
408 23419 : DEALLOCATE (gcij0, gcij1)
409 : END DO
410 14 : CALL dbcsr_iterator_stop(iter)
411 14 : IF (debug_forces) THEN
412 0 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
413 0 : CALL para_env%sum(fodeb)
414 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*H[P]*dS ", fodeb
415 : END IF
416 :
417 14 : IF (xtb_control%tb3_interaction) THEN
418 14 : CALL get_qs_env(qs_env, nkind=nkind)
419 42 : ALLOCATE (xgamma(nkind))
420 48 : DO ikind = 1, nkind
421 34 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
422 48 : CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind))
423 : END DO
424 : ! Diagonal 3rd order correction (DFTB3)
425 14 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
426 : CALL dftb3_diagonal_hessian_force(qs_env, mcharge0, mcharge1, &
427 14 : matrix_p0(1)%matrix, matrix_p1(1)%matrix, xgamma)
428 14 : IF (debug_forces) THEN
429 0 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
430 0 : CALL para_env%sum(fodeb)
431 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*H3[P] ", fodeb
432 : END IF
433 14 : DEALLOCATE (xgamma)
434 : END IF
435 :
436 14 : IF (SIZE(matrix_p0) == 2) THEN
437 : CALL dbcsr_add(matrix_p0(1)%matrix, matrix_p0(2)%matrix, &
438 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
439 : CALL dbcsr_add(matrix_p1(1)%matrix, matrix_p1(2)%matrix, &
440 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
441 : END IF
442 :
443 : ! QMMM
444 14 : IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
445 0 : CPABORT("Not Available")
446 : END IF
447 :
448 14 : DEALLOCATE (gmcharge0, gchrg0, gmcharge1, gchrg1)
449 :
450 14 : CALL timestop(handle)
451 :
452 42 : END SUBROUTINE calc_xtb_ehess_force
453 :
454 : ! **************************************************************************************************
455 : !> \brief ...
456 : !> \param qs_env ...
457 : !> \param mcharge0 ...
458 : !> \param mcharge1 ...
459 : !> \param matrixp0 ...
460 : !> \param matrixp1 ...
461 : !> \param xgamma ...
462 : ! **************************************************************************************************
463 14 : SUBROUTINE dftb3_diagonal_hessian_force(qs_env, mcharge0, mcharge1, &
464 14 : matrixp0, matrixp1, xgamma)
465 :
466 : TYPE(qs_environment_type), POINTER :: qs_env
467 : REAL(dp), DIMENSION(:) :: mcharge0, mcharge1
468 : TYPE(dbcsr_type), POINTER :: matrixp0, matrixp1
469 : REAL(dp), DIMENSION(:) :: xgamma
470 :
471 : CHARACTER(len=*), PARAMETER :: routineN = 'dftb3_diagonal_hessian_force'
472 :
473 : INTEGER :: atom_i, atom_j, blk, handle, i, icol, &
474 : ikind, irow, jkind
475 14 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
476 : LOGICAL :: found
477 : REAL(KIND=dp) :: fi, gmijp, gmijq, ui, uj
478 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, p0block, p1block, sblock
479 14 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
480 : TYPE(dbcsr_iterator_type) :: iter
481 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
482 14 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
483 :
484 14 : CALL timeset(routineN, handle)
485 14 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
486 14 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
487 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
488 14 : kind_of=kind_of, atom_of_kind=atom_of_kind)
489 14 : CALL get_qs_env(qs_env=qs_env, force=force)
490 : ! no k-points; all matrices have been transformed to periodic bsf
491 14 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
492 4695 : DO WHILE (dbcsr_iterator_blocks_left(iter))
493 4681 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
494 4681 : ikind = kind_of(irow)
495 4681 : atom_i = atom_of_kind(irow)
496 4681 : ui = xgamma(ikind)
497 4681 : jkind = kind_of(icol)
498 4681 : atom_j = atom_of_kind(icol)
499 4681 : uj = xgamma(jkind)
500 : !
501 4681 : gmijp = ui*mcharge0(irow)*mcharge1(irow) + uj*mcharge0(icol)*mcharge1(icol)
502 4681 : gmijq = 0.5_dp*ui*mcharge0(irow)**2 + 0.5_dp*uj*mcharge0(icol)**2
503 : !
504 4681 : NULLIFY (p0block)
505 : CALL dbcsr_get_block_p(matrix=matrixp0, &
506 4681 : row=irow, col=icol, block=p0block, found=found)
507 4681 : CPASSERT(found)
508 4681 : NULLIFY (p1block)
509 : CALL dbcsr_get_block_p(matrix=matrixp1, &
510 4681 : row=irow, col=icol, block=p1block, found=found)
511 4681 : CPASSERT(found)
512 18738 : DO i = 1, 3
513 14043 : NULLIFY (dsblock)
514 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
515 14043 : row=irow, col=icol, block=dsblock, found=found)
516 14043 : CPASSERT(found)
517 275571 : fi = gmijp*SUM(p0block*dsblock) + gmijq*SUM(p1block*dsblock)
518 14043 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
519 32767 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
520 : END DO
521 : END DO
522 14 : CALL dbcsr_iterator_stop(iter)
523 :
524 14 : CALL timestop(handle)
525 :
526 28 : END SUBROUTINE dftb3_diagonal_hessian_force
527 :
528 389834 : END MODULE xtb_ehess_force
529 :
|