20 USE dbcsr_api,
ONLY: dbcsr_get_block_p,&
21 dbcsr_iterator_blocks_left,&
22 dbcsr_iterator_next_block,&
23 dbcsr_iterator_start,&
29 ewald_environment_type
49 neighbor_list_iterator_p_type,&
51 neighbor_list_set_p_type
56 #include "./base/base_uses.f90"
62 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'xtb_ehess'
78 TYPE(qs_environment_type),
POINTER :: qs_env
79 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_matrix
80 REAL(
dp),
DIMENSION(:, :) :: charges1
81 REAL(
dp),
DIMENSION(:) :: mcharge1, mcharge
83 CHARACTER(len=*),
PARAMETER :: routinen =
'xtb_coulomb_hessian'
85 INTEGER :: blk, ewald_type, handle, i, ia, iatom, icol, ikind, irow, is, j, jatom, jkind, &
86 la, lb, lmaxa, lmaxb, natom, natorb_a, natorb_b, ni, nj, nkind, nmat, za, zb
87 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: kind_of
88 INTEGER,
DIMENSION(25) :: laoa, laob
89 INTEGER,
DIMENSION(3) :: cellind, periodic
90 LOGICAL :: defined, do_ewald, found
91 REAL(kind=
dp) :: alpha, deth, dr, etaa, etab, gmij, kg, &
93 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: xgamma
94 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: gammab, gcij, gmcharge
95 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: gchrg
96 REAL(kind=
dp),
DIMENSION(3) :: rij
97 REAL(kind=
dp),
DIMENSION(5) :: kappaa, kappab
98 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: ksblock, sblock
99 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
100 TYPE(cell_type),
POINTER :: cell
101 TYPE(dbcsr_iterator_type) :: iter
102 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_s
103 TYPE(dft_control_type),
POINTER :: dft_control
104 TYPE(distribution_1d_type),
POINTER :: local_particles
105 TYPE(ewald_environment_type),
POINTER :: ewald_env
106 TYPE(ewald_pw_type),
POINTER :: ewald_pw
107 TYPE(mp_para_env_type),
POINTER :: para_env
108 TYPE(neighbor_list_iterator_p_type), &
109 DIMENSION(:),
POINTER :: nl_iterator
110 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
112 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
113 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
114 TYPE(virial_type),
POINTER :: virial
115 TYPE(xtb_atom_type),
POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind
116 TYPE(xtb_control_type),
POINTER :: xtb_control
118 CALL timeset(routinen, handle)
121 matrix_s_kp=matrix_s, &
122 qs_kind_set=qs_kind_set, &
123 particle_set=particle_set, &
125 dft_control=dft_control)
127 xtb_control => dft_control%qs_control%xtb_control
129 IF (dft_control%nimages /= 1)
THEN
130 cpabort(
"No kpoints allowed in xTB response calculation")
133 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
135 ALLOCATE (gchrg(natom, 5, nmat))
137 ALLOCATE (gmcharge(natom, nmat))
144 IF (xtb_control%old_coulomb_damping)
THEN
145 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
147 CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
152 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
153 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
155 IF (.NOT. defined .OR. natorb_a < 1) cycle
156 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
158 IF (.NOT. defined .OR. natorb_b < 1) cycle
165 ALLOCATE (gammab(ni, nj))
167 dr = sqrt(sum(rij(:)**2))
168 CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
169 gchrg(iatom, 1:ni, 1) = gchrg(iatom, 1:ni, 1) + matmul(gammab, charges1(jatom, 1:nj))
170 IF (iatom /= jatom)
THEN
171 gchrg(jatom, 1:nj, 1) = gchrg(jatom, 1:nj, 1) + matmul(charges1(iatom, 1:ni), gammab)
179 IF (xtb_control%coulomb_lr)
THEN
180 do_ewald = xtb_control%do_ewald
183 NULLIFY (ewald_env, ewald_pw)
186 ewald_env=ewald_env, ewald_pw=ewald_pw)
187 CALL get_cell(cell=cell, periodic=periodic, deth=deth)
188 CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
189 CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
191 SELECT CASE (ewald_type)
193 cpabort(
"Invalid Ewald type")
195 cpabort(
"Not allowed with DFTB")
197 cpabort(
"Standard Ewald not implemented in DFTB")
199 cpabort(
"PME not implemented in DFTB")
202 gmcharge, mcharge1, .false., virial, .false.)
207 local_particles=local_particles)
208 DO ikind = 1,
SIZE(local_particles%n_el)
209 DO ia = 1, local_particles%n_el(ikind)
210 iatom = local_particles%list(ikind)%array(ia)
211 DO jatom = 1, iatom - 1
212 rij = particle_set(iatom)%r - particle_set(jatom)%r
214 dr = sqrt(sum(rij(:)**2))
215 IF (dr > 1.e-6_dp)
THEN
216 gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge1(jatom)/dr
217 gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge1(iatom)/dr
226 CALL get_qs_env(qs_env=qs_env, para_env=para_env)
227 CALL para_env%sum(gmcharge(:, 1))
228 CALL para_env%sum(gchrg(:, :, 1))
230 IF (xtb_control%coulomb_lr)
THEN
233 gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*
oorootpi*mcharge1(:)
234 IF (any(periodic(:) == 1))
THEN
235 gmcharge(:, 1) = gmcharge(:, 1) -
pi/alpha**2/deth
240 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
244 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
245 DO WHILE (dbcsr_iterator_blocks_left(iter))
246 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
247 ikind = kind_of(irow)
248 jkind = kind_of(icol)
251 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
252 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
258 ALLOCATE (gcij(ni, nj))
263 gcij(i, j) = gchrg(irow, la, 1) + gchrg(icol, lb, 1)
266 gmij = gmcharge(irow, 1) + gmcharge(icol, 1)
267 DO is = 1,
SIZE(ks_matrix)
269 CALL dbcsr_get_block_p(matrix=ks_matrix(is)%matrix, &
270 row=irow, col=icol, block=ksblock, found=found)
272 ksblock = ksblock - gcij*sblock
273 ksblock = ksblock - gmij*sblock
277 CALL dbcsr_iterator_stop(iter)
279 IF (xtb_control%tb3_interaction)
THEN
281 ALLOCATE (xgamma(nkind))
283 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
287 CALL dftb3_diagonal_hessian(qs_env, ks_matrix, mcharge, mcharge1, xgamma)
291 IF (qs_env%qmmm .AND. qs_env%qmmm_periodic)
THEN
292 cpabort(
"QMMM not available in xTB response calculations")
295 DEALLOCATE (gmcharge, gchrg)
297 CALL timestop(handle)
309 SUBROUTINE dftb3_diagonal_hessian(qs_env, ks_matrix, mcharge, mcharge1, xgamma)
311 TYPE(qs_environment_type),
POINTER :: qs_env
312 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_matrix
313 REAL(
dp),
DIMENSION(:) :: mcharge, mcharge1, xgamma
315 CHARACTER(len=*),
PARAMETER :: routinen =
'dftb3_diagonal_hessian'
317 INTEGER :: blk, handle, icol, ikind, irow, is, jkind
318 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: kind_of
320 REAL(kind=
dp) :: gmij, ui, uj
321 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: ksblock, sblock
322 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
323 TYPE(dbcsr_iterator_type) :: iter
324 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_s
325 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
327 CALL timeset(routinen, handle)
329 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
330 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
333 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
334 DO WHILE (dbcsr_iterator_blocks_left(iter))
335 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
336 ikind = kind_of(irow)
338 jkind = kind_of(icol)
340 gmij = ui*mcharge(irow)*mcharge1(irow) + uj*mcharge(icol)*mcharge1(icol)
341 DO is = 1,
SIZE(ks_matrix)
343 CALL dbcsr_get_block_p(matrix=ks_matrix(is)%matrix, &
344 row=irow, col=icol, block=ksblock, found=found)
346 ksblock = ksblock + gmij*sblock
349 CALL dbcsr_iterator_stop(iter)
351 CALL timestop(handle)
353 END SUBROUTINE dftb3_diagonal_hessian
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
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.
Handles all functions related to the CELL.
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Calculation of Ewald contributions in DFTB.
subroutine, public tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial, atprop)
...
subroutine, public tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, gmcharge, mcharge, calculate_forces, virial, use_virial, atprop)
...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public oorootpi
real(kind=dp), parameter, public pi
Interface to the message passing library MPI.
Define the data structure for the particle information.
functions related to the poisson solver on regular grids
integer, parameter, public do_ewald_pme
integer, parameter, public do_ewald_ewald
integer, parameter, public do_ewald_none
integer, parameter, public do_ewald_spme
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.
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.
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)
...
Calculation of Coulomb contributions in xTB.
subroutine, public gamma_rab_sr(gmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
Computes the short-range gamma parameter from Nataga-Mishimoto-Ohno-Klopman formula for xTB WARNING: ...
Calculation of Coulomb Hessian contributions in xTB.
subroutine, public xtb_coulomb_hessian(qs_env, ks_matrix, charges1, mcharge1, mcharge)
...
Definition of the xTB parameter types.
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, occupation, electronegativity, chmax)
...