(git:34ef472)
qs_dftb3_methods.F
Go to the documentation of this file.
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 DFTB3 Terms
10 !> \author JGH
11 ! **************************************************************************************************
13  USE atomic_kind_types, ONLY: atomic_kind_type,&
15  USE atprop_types, ONLY: atprop_type
16  USE cell_types, ONLY: cell_type
17  USE dbcsr_api, ONLY: dbcsr_add,&
18  dbcsr_get_block_p,&
19  dbcsr_iterator_blocks_left,&
20  dbcsr_iterator_next_block,&
21  dbcsr_iterator_start,&
22  dbcsr_iterator_stop,&
23  dbcsr_iterator_type,&
24  dbcsr_p_type
25  USE distribution_1d_types, ONLY: distribution_1d_type
26  USE kinds, ONLY: dp
27  USE kpoint_types, ONLY: get_kpoint_info,&
28  kpoint_type
29  USE message_passing, ONLY: mp_para_env_type
30  USE particle_types, ONLY: particle_type
31  USE qs_energy_types, ONLY: qs_energy_type
32  USE qs_environment_types, ONLY: get_qs_env,&
33  qs_environment_type
34  USE qs_force_types, ONLY: qs_force_type
35  USE qs_kind_types, ONLY: qs_kind_type
39  neighbor_list_iterator_p_type,&
41  neighbor_list_set_p_type
42  USE qs_rho_types, ONLY: qs_rho_get,&
43  qs_rho_type
44  USE sap_kind_types, ONLY: sap_int_type
46  USE virial_types, ONLY: virial_type
47 #include "./base/base_uses.f90"
48 
49  IMPLICIT NONE
50 
51  PRIVATE
52 
53  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb3_methods'
54 
55  PUBLIC :: build_dftb3_diagonal
56 
57 CONTAINS
58 
59 ! **************************************************************************************************
60 !> \brief ...
61 !> \param qs_env ...
62 !> \param ks_matrix ...
63 !> \param rho ...
64 !> \param mcharge ...
65 !> \param energy ...
66 !> \param xgamma ...
67 !> \param zeff ...
68 !> \param sap_int ...
69 !> \param calculate_forces ...
70 !> \param just_energy ...
71 ! **************************************************************************************************
72  SUBROUTINE build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeff, &
73  sap_int, calculate_forces, just_energy)
74 
75  TYPE(qs_environment_type), POINTER :: qs_env
76  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
77  TYPE(qs_rho_type), POINTER :: rho
78  REAL(dp), DIMENSION(:), INTENT(in) :: mcharge
79  TYPE(qs_energy_type), POINTER :: energy
80  REAL(dp), DIMENSION(:), INTENT(in) :: xgamma, zeff
81  TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
82  LOGICAL, INTENT(in) :: calculate_forces, just_energy
83 
84  CHARACTER(len=*), PARAMETER :: routinen = 'build_dftb3_diagonal'
85 
86  INTEGER :: atom_i, atom_j, blk, handle, i, ia, iac, &
87  iatom, ic, icol, ikind, irow, is, &
88  jatom, jkind, natom, nimg, nkind
89  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
90  INTEGER, DIMENSION(3) :: cellind
91  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
92  LOGICAL :: found, use_virial
93  REAL(kind=dp) :: dr, eb3, eloc, fi, gmij, ua, ui, uj
94  REAL(kind=dp), DIMENSION(3) :: fij, rij
95  REAL(kind=dp), DIMENSION(:, :), POINTER :: dsblock, ksblock, pblock, sblock
96  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dsint
97  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
98  TYPE(atprop_type), POINTER :: atprop
99  TYPE(cell_type), POINTER :: cell
100  TYPE(dbcsr_iterator_type) :: iter
101  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
102  TYPE(distribution_1d_type), POINTER :: local_particles
103  TYPE(kpoint_type), POINTER :: kpoints
104  TYPE(mp_para_env_type), POINTER :: para_env
105  TYPE(neighbor_list_iterator_p_type), &
106  DIMENSION(:), POINTER :: nl_iterator
107  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
108  POINTER :: n_list
109  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
110  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
111  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
112  TYPE(virial_type), POINTER :: virial
113 
114  CALL timeset(routinen, handle)
115  NULLIFY (atprop)
116 
117  ! Energy
118  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
119  qs_kind_set=qs_kind_set, atprop=atprop)
120  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
121  kind_of=kind_of, atom_of_kind=atom_of_kind)
122 
123  eb3 = 0.0_dp
124  CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
125  DO ikind = 1, SIZE(local_particles%n_el)
126  ua = xgamma(ikind)
127  DO ia = 1, local_particles%n_el(ikind)
128  iatom = local_particles%list(ikind)%array(ia)
129  eloc = -1.0_dp/6.0_dp*ua*mcharge(iatom)**3
130  eb3 = eb3 + eloc
131  IF (atprop%energy) THEN
132  ! we have to add the part not covered by 0.5*Tr(FP)
133  eloc = -0.5_dp*eloc - 0.25_dp*ua*zeff(ikind)*mcharge(iatom)**2
134  atprop%atecoul(iatom) = atprop%atecoul(iatom) + eloc
135  END IF
136  END DO
137  END DO
138  CALL get_qs_env(qs_env=qs_env, para_env=para_env)
139  CALL para_env%sum(eb3)
140  energy%dftb3 = eb3
141 
142  ! Forces and Virial
143  IF (calculate_forces) THEN
144  CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s, natom=natom, force=force, &
145  cell=cell, virial=virial, particle_set=particle_set)
146  ! virial
147  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
148  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
149  IF (SIZE(matrix_p, 1) == 2) THEN
150  DO ic = 1, SIZE(matrix_p, 2)
151  CALL dbcsr_add(matrix_p(1, ic)%matrix, matrix_p(2, ic)%matrix, &
152  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
153  END DO
154  END IF
155  !
156  nimg = SIZE(matrix_p, 2)
157  NULLIFY (cell_to_index)
158  IF (nimg > 1) THEN
159  NULLIFY (kpoints)
160  CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
161  CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
162  END IF
163  IF (nimg == 1) THEN
164  ! no k-points; all matrices have been transformed to periodic bsf
165  CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
166  DO WHILE (dbcsr_iterator_blocks_left(iter))
167  CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
168  ikind = kind_of(irow)
169  atom_i = atom_of_kind(irow)
170  ui = xgamma(ikind)
171  jkind = kind_of(icol)
172  atom_j = atom_of_kind(icol)
173  uj = xgamma(jkind)
174  !
175  gmij = -0.5_dp*(ui*mcharge(irow)**2 + uj*mcharge(icol)**2)
176  !
177  NULLIFY (pblock)
178  CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
179  row=irow, col=icol, block=pblock, found=found)
180  cpassert(found)
181  DO i = 1, 3
182  NULLIFY (dsblock)
183  CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
184  row=irow, col=icol, block=dsblock, found=found)
185  cpassert(found)
186  fi = -gmij*sum(pblock*dsblock)
187  force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
188  force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
189  fij(i) = fi
190  END DO
191  END DO
192  CALL dbcsr_iterator_stop(iter)
193  ! use dsint list
194  IF (use_virial) THEN
195  cpassert(ASSOCIATED(sap_int))
196  CALL get_qs_env(qs_env, nkind=nkind)
197  DO ikind = 1, nkind
198  DO jkind = 1, nkind
199  iac = ikind + nkind*(jkind - 1)
200  IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
201  ui = xgamma(ikind)
202  uj = xgamma(jkind)
203  DO ia = 1, sap_int(iac)%nalist
204  IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) cycle
205  iatom = sap_int(iac)%alist(ia)%aatom
206  DO ic = 1, sap_int(iac)%alist(ia)%nclist
207  jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
208  rij = sap_int(iac)%alist(ia)%clist(ic)%rac
209  dr = sqrt(sum(rij(:)**2))
210  IF (dr > 1.e-6_dp) THEN
211  dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
212  gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2)
213  icol = max(iatom, jatom)
214  irow = min(iatom, jatom)
215  NULLIFY (pblock)
216  CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
217  row=irow, col=icol, block=pblock, found=found)
218  cpassert(found)
219  DO i = 1, 3
220  IF (irow == iatom) THEN
221  fij(i) = -gmij*sum(pblock*dsint(:, :, i))
222  ELSE
223  fij(i) = -gmij*sum(transpose(pblock)*dsint(:, :, i))
224  END IF
225  END DO
226  fi = 1.0_dp
227  IF (iatom == jatom) fi = 0.5_dp
228  CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
229  IF (atprop%stress) THEN
230  CALL virial_pair_force(atprop%atstress(:, :, irow), fi*0.5_dp, fij, rij)
231  CALL virial_pair_force(atprop%atstress(:, :, icol), fi*0.5_dp, fij, rij)
232  END IF
233  END IF
234  END DO
235  END DO
236  END DO
237  END DO
238  END IF
239  ELSE
240  NULLIFY (n_list)
241  CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
242  CALL neighbor_list_iterator_create(nl_iterator, n_list)
243  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
244  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
245  iatom=iatom, jatom=jatom, r=rij, cell=cellind)
246 
247  dr = sqrt(sum(rij**2))
248  IF (iatom == jatom .AND. dr < 1.0e-6_dp) cycle
249 
250  icol = max(iatom, jatom)
251  irow = min(iatom, jatom)
252 
253  ic = cell_to_index(cellind(1), cellind(2), cellind(3))
254  cpassert(ic > 0)
255 
256  ikind = kind_of(iatom)
257  atom_i = atom_of_kind(iatom)
258  ui = xgamma(ikind)
259  jkind = kind_of(jatom)
260  atom_j = atom_of_kind(jatom)
261  uj = xgamma(jkind)
262  !
263  gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2)
264  !
265  NULLIFY (pblock)
266  CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
267  row=irow, col=icol, block=pblock, found=found)
268  cpassert(found)
269  DO i = 1, 3
270  NULLIFY (dsblock)
271  CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
272  row=irow, col=icol, block=dsblock, found=found)
273  cpassert(found)
274  IF (irow == iatom) THEN
275  fi = -gmij*sum(pblock*dsblock)
276  ELSE
277  fi = gmij*sum(pblock*dsblock)
278  END IF
279  force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
280  force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
281  fij(i) = fi
282  END DO
283  IF (use_virial) THEN
284  fi = 1.0_dp
285  IF (iatom == jatom) fi = 0.5_dp
286  CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
287  IF (atprop%stress) THEN
288  CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij)
289  CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij)
290  END IF
291  END IF
292 
293  END DO
294  CALL neighbor_list_iterator_release(nl_iterator)
295  !
296  END IF
297 
298  IF (SIZE(matrix_p, 1) == 2) THEN
299  DO ic = 1, SIZE(matrix_p, 2)
300  CALL dbcsr_add(matrix_p(1, ic)%matrix, matrix_p(2, ic)%matrix, &
301  alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
302  END DO
303  END IF
304  END IF
305 
306  ! KS matrix
307  IF (.NOT. just_energy) THEN
308  CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s, natom=natom)
309  !
310  nimg = SIZE(ks_matrix, 2)
311  NULLIFY (cell_to_index)
312  IF (nimg > 1) THEN
313  NULLIFY (kpoints)
314  CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
315  CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
316  END IF
317 
318  IF (nimg == 1) THEN
319  ! no k-points; all matrices have been transformed to periodic bsf
320  CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
321  DO WHILE (dbcsr_iterator_blocks_left(iter))
322  CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
323  ikind = kind_of(irow)
324  ui = xgamma(ikind)
325  jkind = kind_of(icol)
326  uj = xgamma(jkind)
327  gmij = -0.5_dp*(ui*mcharge(irow)**2 + uj*mcharge(icol)**2)
328  DO is = 1, SIZE(ks_matrix, 1)
329  NULLIFY (ksblock)
330  CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
331  row=irow, col=icol, block=ksblock, found=found)
332  cpassert(found)
333  ksblock = ksblock - 0.5_dp*gmij*sblock
334  END DO
335  END DO
336  CALL dbcsr_iterator_stop(iter)
337  ELSE
338  NULLIFY (n_list)
339  CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
340  CALL neighbor_list_iterator_create(nl_iterator, n_list)
341  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
342  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
343  iatom=iatom, jatom=jatom, r=rij, cell=cellind)
344 
345  icol = max(iatom, jatom)
346  irow = min(iatom, jatom)
347 
348  ic = cell_to_index(cellind(1), cellind(2), cellind(3))
349  cpassert(ic > 0)
350 
351  ikind = kind_of(iatom)
352  ui = xgamma(ikind)
353  jkind = kind_of(jatom)
354  uj = xgamma(jkind)
355  gmij = -0.5_dp*(ui*mcharge(iatom)**2 + uj*mcharge(jatom)**2)
356  !
357  NULLIFY (sblock)
358  CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
359  row=irow, col=icol, block=sblock, found=found)
360  cpassert(found)
361  DO is = 1, SIZE(ks_matrix, 1)
362  NULLIFY (ksblock)
363  CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
364  row=irow, col=icol, block=ksblock, found=found)
365  cpassert(found)
366  ksblock = ksblock - 0.5_dp*gmij*sblock
367  END DO
368 
369  END DO
370  CALL neighbor_list_iterator_release(nl_iterator)
371  !
372  END IF
373  END IF
374 
375  CALL timestop(handle)
376 
377  END SUBROUTINE build_dftb3_diagonal
378 
379 END MODULE qs_dftb3_methods
380 
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Holds information on atomic properties.
Definition: atprop_types.F:14
Handles all functions related to the CELL.
Definition: cell_types.F:15
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
Interface to the message passing library MPI.
Define the data structure for the particle information.
Calculation of DFTB3 Terms.
subroutine, public build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeff, sap_int, calculate_forces, just_energy)
...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
General overlap type integrals containers.
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.