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 : CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
183 14 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
184 24708 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
185 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
186 24694 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
187 24694 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
188 24694 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
189 24694 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
190 24694 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
191 24694 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
192 24694 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
193 : ! atomic parameters
194 24694 : CALL get_xtb_atom_param(xtb_atom_a, eta=etaa, lmax=lmaxa, kappa=kappaa, rcut=rcuta)
195 24694 : CALL get_xtb_atom_param(xtb_atom_b, eta=etab, lmax=lmaxb, kappa=kappab, rcut=rcutb)
196 : ! gamma matrix
197 24694 : ni = lmaxa + 1
198 24694 : nj = lmaxb + 1
199 98776 : ALLOCATE (gammab(ni, nj))
200 24694 : rcut = rcuta + rcutb
201 98776 : dr = SQRT(SUM(rij(:)**2))
202 24694 : CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
203 250847 : gchrg0(iatom, 1:ni, 1) = gchrg0(iatom, 1:ni, 1) + MATMUL(gammab, charges0(jatom, 1:nj))
204 250847 : gchrg1(iatom, 1:ni, 1) = gchrg1(iatom, 1:ni, 1) + MATMUL(gammab, charges1(jatom, 1:nj))
205 24694 : IF (iatom /= jatom) THEN
206 285738 : gchrg0(jatom, 1:nj, 1) = gchrg0(jatom, 1:nj, 1) + MATMUL(charges0(iatom, 1:ni), gammab)
207 285738 : gchrg1(jatom, 1:nj, 1) = gchrg1(jatom, 1:nj, 1) + MATMUL(charges1(iatom, 1:ni), gammab)
208 : END IF
209 24694 : IF (dr > 1.e-6_dp) THEN
210 24577 : CALL dgamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
211 98308 : DO i = 1, 3
212 : gchrg0(iatom, 1:ni, i + 1) = gchrg0(iatom, 1:ni, i + 1) &
213 748626 : + MATMUL(gammab, charges0(jatom, 1:nj))*rij(i)/dr
214 : gchrg1(iatom, 1:ni, i + 1) = gchrg1(iatom, 1:ni, i + 1) &
215 748626 : + MATMUL(gammab, charges1(jatom, 1:nj))*rij(i)/dr
216 98308 : IF (iatom /= jatom) THEN
217 : gchrg0(jatom, 1:nj, i + 1) = gchrg0(jatom, 1:nj, i + 1) &
218 857214 : - MATMUL(charges0(iatom, 1:ni), gammab)*rij(i)/dr
219 : gchrg1(jatom, 1:nj, i + 1) = gchrg1(jatom, 1:nj, i + 1) &
220 857214 : - MATMUL(charges1(iatom, 1:ni), gammab)*rij(i)/dr
221 : END IF
222 : END DO
223 : END IF
224 98776 : DEALLOCATE (gammab)
225 : END DO
226 14 : CALL neighbor_list_iterator_release(nl_iterator)
227 :
228 : ! 1/R contribution
229 :
230 14 : IF (xtb_control%coulomb_lr) THEN
231 14 : do_ewald = xtb_control%do_ewald
232 14 : IF (do_ewald) THEN
233 : ! Ewald sum
234 6 : NULLIFY (ewald_env, ewald_pw)
235 : CALL get_qs_env(qs_env=qs_env, &
236 6 : ewald_env=ewald_env, ewald_pw=ewald_pw)
237 6 : CALL get_cell(cell=cell, periodic=periodic, deth=deth)
238 6 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
239 6 : CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
240 6 : CALL tb_ewald_overlap(gmcharge0, mcharge0, alpha, n_list, virial, use_virial)
241 6 : CALL tb_ewald_overlap(gmcharge1, mcharge1, alpha, n_list, virial, use_virial)
242 0 : SELECT CASE (ewald_type)
243 : CASE DEFAULT
244 0 : CPABORT("Invalid Ewald type")
245 : CASE (do_ewald_none)
246 0 : CPABORT("Not allowed with DFTB")
247 : CASE (do_ewald_ewald)
248 0 : CPABORT("Standard Ewald not implemented in DFTB")
249 : CASE (do_ewald_pme)
250 0 : CPABORT("PME not implemented in DFTB")
251 : CASE (do_ewald_spme)
252 6 : CALL tb_spme_zforce(ewald_env, ewald_pw, particle_set, cell, gmcharge0, mcharge0)
253 12 : CALL tb_spme_zforce(ewald_env, ewald_pw, particle_set, cell, gmcharge1, mcharge1)
254 : END SELECT
255 : ELSE
256 : ! direct sum
257 8 : CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
258 28 : DO ikind = 1, SIZE(local_particles%n_el)
259 42 : DO ia = 1, local_particles%n_el(ikind)
260 14 : iatom = local_particles%list(ikind)%array(ia)
261 52 : DO jatom = 1, iatom - 1
262 72 : rij = particle_set(iatom)%r - particle_set(jatom)%r
263 72 : rij = pbc(rij, cell)
264 72 : dr = SQRT(SUM(rij(:)**2))
265 32 : IF (dr > 1.e-6_dp) THEN
266 18 : gmcharge0(iatom, 1) = gmcharge0(iatom, 1) + mcharge0(jatom)/dr
267 18 : gmcharge0(jatom, 1) = gmcharge0(jatom, 1) + mcharge0(iatom)/dr
268 18 : gmcharge1(iatom, 1) = gmcharge1(iatom, 1) + mcharge1(jatom)/dr
269 18 : gmcharge1(jatom, 1) = gmcharge1(jatom, 1) + mcharge1(iatom)/dr
270 72 : DO i = 2, nmat
271 54 : gmcharge0(iatom, i) = gmcharge0(iatom, i) + rij(i - 1)*mcharge0(jatom)/dr**3
272 54 : gmcharge0(jatom, i) = gmcharge0(jatom, i) - rij(i - 1)*mcharge0(iatom)/dr**3
273 54 : gmcharge1(iatom, i) = gmcharge1(iatom, i) + rij(i - 1)*mcharge1(jatom)/dr**3
274 72 : gmcharge1(jatom, i) = gmcharge1(jatom, i) - rij(i - 1)*mcharge1(iatom)/dr**3
275 : END DO
276 : END IF
277 : END DO
278 : END DO
279 : END DO
280 : CPASSERT(.NOT. use_virial)
281 : END IF
282 : END IF
283 :
284 : ! global sum of gamma*p arrays
285 : CALL get_qs_env(qs_env=qs_env, &
286 : atomic_kind_set=atomic_kind_set, &
287 14 : force=force, para_env=para_env)
288 14 : CALL para_env%sum(gmcharge0(:, 1))
289 14 : CALL para_env%sum(gchrg0(:, :, 1))
290 14 : CALL para_env%sum(gmcharge1(:, 1))
291 14 : CALL para_env%sum(gchrg1(:, :, 1))
292 :
293 14 : IF (xtb_control%coulomb_lr) THEN
294 14 : IF (do_ewald) THEN
295 : ! add self charge interaction and background charge contribution
296 212 : gmcharge0(:, 1) = gmcharge0(:, 1) - 2._dp*alpha*oorootpi*mcharge0(:)
297 12 : IF (ANY(periodic(:) == 1)) THEN
298 202 : gmcharge0(:, 1) = gmcharge0(:, 1) - pi/alpha**2/deth
299 : END IF
300 212 : gmcharge1(:, 1) = gmcharge1(:, 1) - 2._dp*alpha*oorootpi*mcharge1(:)
301 12 : IF (ANY(periodic(:) == 1)) THEN
302 202 : gmcharge1(:, 1) = gmcharge1(:, 1) - pi/alpha**2/deth
303 : END IF
304 : END IF
305 : END IF
306 :
307 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
308 : kind_of=kind_of, &
309 14 : atom_of_kind=atom_of_kind)
310 :
311 14 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
312 248 : DO iatom = 1, natom
313 234 : ikind = kind_of(iatom)
314 234 : atom_i = atom_of_kind(iatom)
315 234 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
316 234 : CALL get_xtb_atom_param(xtb_kind, lmax=ni)
317 234 : ni = ni + 1
318 : ! short range
319 234 : fij = 0.0_dp
320 936 : DO i = 1, 3
321 : fij(i) = SUM(charges0(iatom, 1:ni)*gchrg1(iatom, 1:ni, i + 1)) + &
322 2832 : SUM(charges1(iatom, 1:ni)*gchrg0(iatom, 1:ni, i + 1))
323 : END DO
324 234 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
325 234 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
326 234 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
327 : ! long range
328 234 : fij = 0.0_dp
329 936 : DO i = 1, 3
330 : fij(i) = gmcharge1(iatom, i + 1)*mcharge0(iatom) + &
331 936 : gmcharge0(iatom, i + 1)*mcharge1(iatom)
332 : END DO
333 234 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
334 234 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
335 482 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
336 : END DO
337 14 : IF (debug_forces) THEN
338 0 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
339 0 : CALL para_env%sum(fodeb)
340 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dH[Pz] ", fodeb
341 : END IF
342 :
343 14 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
344 :
345 14 : IF (SIZE(matrix_p0) == 2) THEN
346 : CALL dbcsr_add(matrix_p0(1)%matrix, matrix_p0(2)%matrix, &
347 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
348 : CALL dbcsr_add(matrix_p1(1)%matrix, matrix_p1(2)%matrix, &
349 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
350 : END IF
351 :
352 : ! no k-points; all matrices have been transformed to periodic bsf
353 14 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
354 14 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
355 4695 : DO WHILE (dbcsr_iterator_blocks_left(iter))
356 4681 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
357 4681 : ikind = kind_of(irow)
358 4681 : jkind = kind_of(icol)
359 :
360 : ! atomic parameters
361 4681 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
362 4681 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
363 4681 : CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
364 4681 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
365 :
366 4681 : ni = SIZE(sblock, 1)
367 4681 : nj = SIZE(sblock, 2)
368 18724 : ALLOCATE (gcij0(ni, nj))
369 14043 : ALLOCATE (gcij1(ni, nj))
370 19209 : DO i = 1, ni
371 52401 : DO j = 1, nj
372 33192 : la = laoa(i) + 1
373 33192 : lb = laob(j) + 1
374 33192 : gcij0(i, j) = 0.5_dp*(gchrg0(irow, la, 1) + gchrg0(icol, lb, 1))
375 47720 : gcij1(i, j) = 0.5_dp*(gchrg1(irow, la, 1) + gchrg1(icol, lb, 1))
376 : END DO
377 : END DO
378 4681 : gmij0 = 0.5_dp*(gmcharge0(irow, 1) + gmcharge0(icol, 1))
379 4681 : gmij1 = 0.5_dp*(gmcharge1(irow, 1) + gmcharge1(icol, 1))
380 4681 : atom_i = atom_of_kind(irow)
381 4681 : atom_j = atom_of_kind(icol)
382 4681 : NULLIFY (pblock0)
383 : CALL dbcsr_get_block_p(matrix=matrix_p0(1)%matrix, &
384 4681 : row=irow, col=icol, block=pblock0, found=found)
385 4681 : CPASSERT(found)
386 4681 : NULLIFY (pblock1)
387 : CALL dbcsr_get_block_p(matrix=matrix_p1(1)%matrix, &
388 4681 : row=irow, col=icol, block=pblock1, found=found)
389 4681 : CPASSERT(found)
390 18724 : DO i = 1, 3
391 14043 : NULLIFY (dsblock)
392 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
393 14043 : row=irow, col=icol, block=dsblock, found=found)
394 14043 : CPASSERT(found)
395 : ! short range
396 275571 : fi = -2.0_dp*SUM(pblock0*dsblock*gcij1) - 2.0_dp*SUM(pblock1*dsblock*gcij0)
397 14043 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
398 14043 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
399 : ! long range
400 275571 : fi = -2.0_dp*gmij1*SUM(pblock0*dsblock) - 2.0_dp*gmij0*SUM(pblock1*dsblock)
401 14043 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
402 32767 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
403 : END DO
404 23419 : DEALLOCATE (gcij0, gcij1)
405 : END DO
406 14 : CALL dbcsr_iterator_stop(iter)
407 14 : IF (debug_forces) THEN
408 0 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
409 0 : CALL para_env%sum(fodeb)
410 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*H[P]*dS ", fodeb
411 : END IF
412 :
413 14 : IF (xtb_control%tb3_interaction) THEN
414 14 : CALL get_qs_env(qs_env, nkind=nkind)
415 42 : ALLOCATE (xgamma(nkind))
416 48 : DO ikind = 1, nkind
417 34 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
418 48 : CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind))
419 : END DO
420 : ! Diagonal 3rd order correction (DFTB3)
421 14 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
422 : CALL dftb3_diagonal_hessian_force(qs_env, mcharge0, mcharge1, &
423 14 : matrix_p0(1)%matrix, matrix_p1(1)%matrix, xgamma)
424 14 : IF (debug_forces) THEN
425 0 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
426 0 : CALL para_env%sum(fodeb)
427 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*H3[P] ", fodeb
428 : END IF
429 14 : DEALLOCATE (xgamma)
430 : END IF
431 :
432 14 : IF (SIZE(matrix_p0) == 2) THEN
433 : CALL dbcsr_add(matrix_p0(1)%matrix, matrix_p0(2)%matrix, &
434 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
435 : CALL dbcsr_add(matrix_p1(1)%matrix, matrix_p1(2)%matrix, &
436 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
437 : END IF
438 :
439 : ! QMMM
440 14 : IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
441 0 : CPABORT("Not Available")
442 : END IF
443 :
444 14 : DEALLOCATE (gmcharge0, gchrg0, gmcharge1, gchrg1)
445 :
446 14 : CALL timestop(handle)
447 :
448 42 : END SUBROUTINE calc_xtb_ehess_force
449 :
450 : ! **************************************************************************************************
451 : !> \brief ...
452 : !> \param qs_env ...
453 : !> \param mcharge0 ...
454 : !> \param mcharge1 ...
455 : !> \param matrixp0 ...
456 : !> \param matrixp1 ...
457 : !> \param xgamma ...
458 : ! **************************************************************************************************
459 14 : SUBROUTINE dftb3_diagonal_hessian_force(qs_env, mcharge0, mcharge1, &
460 14 : matrixp0, matrixp1, xgamma)
461 :
462 : TYPE(qs_environment_type), POINTER :: qs_env
463 : REAL(dp), DIMENSION(:) :: mcharge0, mcharge1
464 : TYPE(dbcsr_type), POINTER :: matrixp0, matrixp1
465 : REAL(dp), DIMENSION(:) :: xgamma
466 :
467 : CHARACTER(len=*), PARAMETER :: routineN = 'dftb3_diagonal_hessian_force'
468 :
469 : INTEGER :: atom_i, atom_j, blk, handle, i, icol, &
470 : ikind, irow, jkind
471 14 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
472 : LOGICAL :: found
473 : REAL(KIND=dp) :: fi, gmijp, gmijq, ui, uj
474 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, p0block, p1block, sblock
475 14 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
476 : TYPE(dbcsr_iterator_type) :: iter
477 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
478 14 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
479 :
480 14 : CALL timeset(routineN, handle)
481 14 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
482 14 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
483 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
484 14 : kind_of=kind_of, atom_of_kind=atom_of_kind)
485 14 : CALL get_qs_env(qs_env=qs_env, force=force)
486 : ! no k-points; all matrices have been transformed to periodic bsf
487 14 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
488 4695 : DO WHILE (dbcsr_iterator_blocks_left(iter))
489 4681 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
490 4681 : ikind = kind_of(irow)
491 4681 : atom_i = atom_of_kind(irow)
492 4681 : ui = xgamma(ikind)
493 4681 : jkind = kind_of(icol)
494 4681 : atom_j = atom_of_kind(icol)
495 4681 : uj = xgamma(jkind)
496 : !
497 4681 : gmijp = ui*mcharge0(irow)*mcharge1(irow) + uj*mcharge0(icol)*mcharge1(icol)
498 4681 : gmijq = 0.5_dp*ui*mcharge0(irow)**2 + 0.5_dp*uj*mcharge0(icol)**2
499 : !
500 4681 : NULLIFY (p0block)
501 : CALL dbcsr_get_block_p(matrix=matrixp0, &
502 4681 : row=irow, col=icol, block=p0block, found=found)
503 4681 : CPASSERT(found)
504 4681 : NULLIFY (p1block)
505 : CALL dbcsr_get_block_p(matrix=matrixp1, &
506 4681 : row=irow, col=icol, block=p1block, found=found)
507 4681 : CPASSERT(found)
508 18738 : DO i = 1, 3
509 14043 : NULLIFY (dsblock)
510 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
511 14043 : row=irow, col=icol, block=dsblock, found=found)
512 14043 : CPASSERT(found)
513 275571 : fi = gmijp*SUM(p0block*dsblock) + gmijq*SUM(p1block*dsblock)
514 14043 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
515 32767 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
516 : END DO
517 : END DO
518 14 : CALL dbcsr_iterator_stop(iter)
519 :
520 14 : CALL timestop(handle)
521 :
522 28 : END SUBROUTINE dftb3_diagonal_hessian_force
523 :
524 389834 : END MODULE xtb_ehess_force
525 :
|
