(git:34ef472)
qs_dftb_dispersion.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 dispersion in DFTB
10 !> \author JGH
11 ! **************************************************************************************************
12 MODULE qs_dftb_dispersion
13  USE atomic_kind_types, ONLY: atomic_kind_type,&
14  get_atomic_kind_set
15  USE atprop_types, ONLY: atprop_array_init,&
16  atprop_type
17  USE cp_control_types, ONLY: dft_control_type,&
18  dftb_control_type
19  USE input_constants, ONLY: dispersion_d2,&
20  dispersion_d3,&
21  dispersion_d3bj,&
22  dispersion_uff
23  USE kinds, ONLY: dp
24  USE message_passing, ONLY: mp_para_env_type
25  USE particle_types, ONLY: particle_type
26  USE qs_dftb_types, ONLY: qs_dftb_atom_type,&
27  qs_dftb_pairpot_type
28  USE qs_dftb_utils, ONLY: get_dftb_atom_param
29  USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
30  USE qs_dispersion_types, ONLY: qs_dispersion_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: get_qs_kind,&
36  qs_kind_type
37  USE qs_neighbor_list_types, ONLY: get_iterator_info,&
38  neighbor_list_iterate,&
39  neighbor_list_iterator_create,&
40  neighbor_list_iterator_p_type,&
41  neighbor_list_iterator_release,&
42  neighbor_list_set_p_type
43  USE virial_methods, ONLY: virial_pair_force
44  USE virial_types, ONLY: virial_type
45 #include "./base/base_uses.f90"
46 
47  IMPLICIT NONE
48 
49  PRIVATE
50 
51  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_dispersion'
52 
53  PUBLIC :: calculate_dftb_dispersion
54 
55 CONTAINS
56 
57 ! **************************************************************************************************
58 !> \brief ...
59 !> \param qs_env ...
60 !> \param para_env ...
61 !> \param calculate_forces ...
62 ! **************************************************************************************************
63  SUBROUTINE calculate_dftb_dispersion(qs_env, para_env, calculate_forces)
64 
65  TYPE(qs_environment_type), POINTER :: qs_env
66  TYPE(mp_para_env_type), POINTER :: para_env
67  LOGICAL, INTENT(IN) :: calculate_forces
68 
69  TYPE(dft_control_type), POINTER :: dft_control
70  TYPE(dftb_control_type), POINTER :: dftb_control
71  TYPE(qs_dispersion_type), POINTER :: dispersion_env
72  TYPE(qs_energy_type), POINTER :: energy
73 
74  CALL get_qs_env(qs_env=qs_env, &
75  energy=energy, &
76  dft_control=dft_control)
77 
78  energy%dispersion = 0._dp
79 
80  dftb_control => dft_control%qs_control%dftb_control
81  IF (dftb_control%dispersion) THEN
82  SELECT CASE (dftb_control%dispersion_type)
83  CASE (dispersion_uff)
84  CALL calculate_dispersion_uff(qs_env, para_env, calculate_forces)
85  CASE (dispersion_d3, dispersion_d3bj, dispersion_d2)
86  CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
87  CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
88  energy%dispersion, calculate_forces)
89  CASE DEFAULT
90  cpabort("")
91  END SELECT
92  END IF
93 
94  END SUBROUTINE calculate_dftb_dispersion
95 
96 ! **************************************************************************************************
97 !> \brief ...
98 !> \param qs_env ...
99 !> \param para_env ...
100 !> \param calculate_forces ...
101 ! **************************************************************************************************
102  SUBROUTINE calculate_dispersion_uff(qs_env, para_env, calculate_forces)
103 
104  TYPE(qs_environment_type), POINTER :: qs_env
105  TYPE(mp_para_env_type), POINTER :: para_env
106  LOGICAL, INTENT(IN) :: calculate_forces
107 
108  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_dispersion_uff'
109 
110  INTEGER :: atom_a, atom_b, handle, iatom, ikind, &
111  jatom, jkind, nkind
112  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
113  LOGICAL :: use_virial
114  LOGICAL, ALLOCATABLE, DIMENSION(:) :: define_kind
115  REAL(dp), ALLOCATABLE, DIMENSION(:) :: rc_kind
116  REAL(kind=dp) :: a, b, c, devdw, dij, dr, eij, evdw, fac, &
117  rc, x0ij, xij, xp
118  REAL(kind=dp), DIMENSION(3) :: fdij, rij
119  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
120  TYPE(atprop_type), POINTER :: atprop
121  TYPE(dft_control_type), POINTER :: dft_control
122  TYPE(dftb_control_type), POINTER :: dftb_control
123  TYPE(neighbor_list_iterator_p_type), &
124  DIMENSION(:), POINTER :: nl_iterator
125  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
126  POINTER :: sab_vdw
127  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
128  TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a
129  TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
130  POINTER :: dftb_potential
131  TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij
132  TYPE(qs_energy_type), POINTER :: energy
133  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
134  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
135  TYPE(virial_type), POINTER :: virial
136 
137  CALL timeset(routinen, handle)
138 
139  NULLIFY (atomic_kind_set, sab_vdw, atprop)
140 
141  CALL get_qs_env(qs_env=qs_env, &
142  energy=energy, &
143  atomic_kind_set=atomic_kind_set, &
144  qs_kind_set=qs_kind_set, &
145  virial=virial, atprop=atprop, &
146  dft_control=dft_control)
147 
148  energy%dispersion = 0._dp
149 
150  dftb_control => dft_control%qs_control%dftb_control
151 
152  IF (dftb_control%dispersion) THEN
153 
154  NULLIFY (dftb_potential)
155  CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential)
156  IF (calculate_forces) THEN
157  NULLIFY (force, particle_set)
158  CALL get_qs_env(qs_env=qs_env, &
159  particle_set=particle_set, &
160  force=force)
161  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
162  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
163  ELSE
164  use_virial = .false.
165  END IF
166  nkind = SIZE(atomic_kind_set)
167  ALLOCATE (define_kind(nkind), rc_kind(nkind))
168  DO ikind = 1, nkind
169  CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
170  CALL get_dftb_atom_param(dftb_kind_a, defined=define_kind(ikind), rcdisp=rc_kind(ikind))
171  END DO
172 
173  evdw = 0._dp
174 
175  IF (atprop%energy) THEN
176  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
177  CALL atprop_array_init(atprop%atevdw, natom=SIZE(particle_set))
178  END IF
179 
180  CALL get_qs_env(qs_env=qs_env, sab_vdw=sab_vdw)
181  CALL neighbor_list_iterator_create(nl_iterator, sab_vdw)
182  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
183  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
184  IF ((.NOT. define_kind(ikind)) .OR. (.NOT. define_kind(jkind))) cycle
185  rc = rc_kind(ikind) + rc_kind(jkind)
186  ! vdW potential
187  dr = sqrt(sum(rij(:)**2))
188  IF (dr <= rc .AND. dr > 0.001_dp) THEN
189  fac = 1._dp
190  IF (iatom == jatom) fac = 0.5_dp
191  ! retrieve information on potential
192  dftb_param_ij => dftb_potential(ikind, jkind)
193  ! vdW parameter
194  xij = dftb_param_ij%xij
195  dij = dftb_param_ij%dij
196  x0ij = dftb_param_ij%x0ij
197  a = dftb_param_ij%a
198  b = dftb_param_ij%b
199  c = dftb_param_ij%c
200  IF (dr > x0ij) THEN
201  ! This is the standard London contribution.
202  ! UFF1 - Eq. 20 (long-range)
203  xp = xij/dr
204  eij = dij*(-2._dp*xp**6 + xp**12)*fac
205  evdw = evdw + eij
206  IF (calculate_forces .AND. (dr > 0.001_dp)) THEN
207  devdw = dij*12._dp*(xp**6 - xp**12)/dr*fac
208  atom_a = atom_of_kind(iatom)
209  atom_b = atom_of_kind(jatom)
210  fdij(:) = devdw*rij(:)/dr
211  force(ikind)%dispersion(:, atom_a) = &
212  force(ikind)%dispersion(:, atom_a) - fdij(:)
213  force(jkind)%dispersion(:, atom_b) = &
214  force(jkind)%dispersion(:, atom_b) + fdij(:)
215  END IF
216  ELSE
217  ! Shorter distance.
218  ! London contribution should converge to a finite value.
219  ! Using a parabola of the form (y = A - Bx**5 -Cx**10).
220  ! Analytic parameters by forcing energy, first and second
221  ! derivatives to be continuous.
222  eij = (a - b*dr**5 - c*dr**10)*fac
223  evdw = evdw + eij
224  IF (calculate_forces .AND. (dr > 0.001_dp)) THEN
225  atom_a = atom_of_kind(iatom)
226  atom_b = atom_of_kind(jatom)
227  devdw = (-5*b*dr**4 - 10*c*dr**9)*fac
228  fdij(:) = devdw*rij(:)/dr
229  force(ikind)%dispersion(:, atom_a) = &
230  force(ikind)%dispersion(:, atom_a) - fdij(:)
231  force(jkind)%dispersion(:, atom_b) = &
232  force(jkind)%dispersion(:, atom_b) + fdij(:)
233  END IF
234  END IF
235  IF (atprop%energy) THEN
236  atprop%atevdw(iatom) = atprop%atevdw(iatom) + 0.5_dp*eij
237  atprop%atevdw(jatom) = atprop%atevdw(jatom) + 0.5_dp*eij
238  END IF
239  IF (calculate_forces .AND. (dr > 0.001_dp) .AND. use_virial) THEN
240  CALL virial_pair_force(virial%pv_virial, -1._dp, fdij, rij)
241  IF (atprop%stress) THEN
242  CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fdij, rij)
243  CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fdij, rij)
244  END IF
245  END IF
246  END IF
247  END DO
248  CALL neighbor_list_iterator_release(nl_iterator)
249 
250  ! set dispersion energy
251  CALL para_env%sum(evdw)
252  energy%dispersion = evdw
253 
254  END IF
255 
256  CALL timestop(handle)
257 
258  END SUBROUTINE calculate_dispersion_uff
259 
260 END MODULE qs_dftb_dispersion
261 
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
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.
Definition: atomic_kind_types.F:223
atprop_types
Holds information on atomic properties.
Definition: atprop_types.F:14
atprop_types::atprop_array_init
subroutine, public atprop_array_init(atarray, natom)
...
Definition: atprop_types.F:98
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::dispersion_d3
integer, parameter, public dispersion_d3
Definition: input_constants.F:1179
input_constants::dispersion_uff
integer, parameter, public dispersion_uff
Definition: input_constants.F:1179
input_constants::dispersion_d3bj
integer, parameter, public dispersion_d3bj
Definition: input_constants.F:1179
input_constants::dispersion_d2
integer, parameter, public dispersion_d2
Definition: input_constants.F:1179
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_dftb_dispersion
Calculation of dispersion in DFTB.
Definition: qs_dftb_dispersion.F:12
qs_dftb_dispersion::calculate_dftb_dispersion
subroutine, public calculate_dftb_dispersion(qs_env, para_env, calculate_forces)
...
Definition: qs_dftb_dispersion.F:64
qs_dftb_types
Definition of the DFTB parameter types.
Definition: qs_dftb_types.F:12
qs_dftb_utils
Working with the DFTB parameter types.
Definition: qs_dftb_utils.F:12
qs_dftb_utils::get_dftb_atom_param
subroutine, public get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
Definition: qs_dftb_utils.F:123
qs_dispersion_pairpot
Calculation of dispersion using pair potentials.
Definition: qs_dispersion_pairpot.F:12
qs_dispersion_pairpot::calculate_dispersion_pairpot
subroutine, public calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
...
Definition: qs_dispersion_pairpot.F:1245
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition: qs_dispersion_types.F:12
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
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.
Definition: qs_environment_types.F:502
qs_force_types
Definition: qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
virial_methods
Definition: virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition: virial_methods.F:139
virial_types
Definition: virial_types.F:13