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