(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
12MODULE 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
55CONTAINS
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
260END 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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_control_types::dftb_control_type
Definition cp_control_types.F:100
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_dftb_types::qs_dftb_atom_type
Definition qs_dftb_types.F:27
qs_dftb_types::qs_dftb_pairpot_type
Definition qs_dftb_types.F:47
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:31
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
virial_types::virial_type
Definition virial_types.F:26