85 charges1, mcharge1, debug_forces)
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
95 CHARACTER(len=*),
PARAMETER :: routinen =
'calc_xtb_ehess_force'
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, &
100 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 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: xgamma
108 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: gammab, gcij0, gcij1, gmcharge0, &
110 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 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: dsblock, pblock0, pblock1, sblock
117 TYPE(dbcsr_iterator_type) :: iter
118 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_s
125 DIMENSION(:),
POINTER :: nl_iterator
131 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
134 TYPE(
xtb_atom_type),
POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind
137 CALL timeset(routinen, handle)
140 IF (logger%para_env%is_source())
THEN
146 cpassert(
ASSOCIATED(matrix_p1))
149 qs_kind_set=qs_kind_set, &
150 particle_set=particle_set, &
155 dft_control=dft_control)
157 xtb_control => dft_control%qs_control%xtb_control
159 calculate_forces = .true.
160 just_energy = .false.
163 nimg = dft_control%nimages
165 cpabort(
'xTB-sTDA forces for k-points not available')
168 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
169 ALLOCATE (gchrg0(natom, 5, nmat))
171 ALLOCATE (gmcharge0(natom, nmat))
173 ALLOCATE (gchrg1(natom, 5, nmat))
175 ALLOCATE (gmcharge1(natom, nmat))
182 IF (xtb_control%old_coulomb_damping)
THEN
183 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
185 CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
190 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
191 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
193 IF (.NOT. defined .OR. natorb_a < 1) cycle
194 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
196 IF (.NOT. defined .OR. natorb_b < 1) cycle
203 ALLOCATE (gammab(ni, nj))
205 dr = sqrt(sum(rij(:)**2))
206 CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
207 gchrg0(iatom, 1:ni, 1) = gchrg0(iatom, 1:ni, 1) + matmul(gammab, charges0(jatom, 1:nj))
208 gchrg1(iatom, 1:ni, 1) = gchrg1(iatom, 1:ni, 1) + matmul(gammab, charges1(jatom, 1:nj))
209 IF (iatom /= jatom)
THEN
210 gchrg0(jatom, 1:nj, 1) = gchrg0(jatom, 1:nj, 1) + matmul(charges0(iatom, 1:ni), gammab)
211 gchrg1(jatom, 1:nj, 1) = gchrg1(jatom, 1:nj, 1) + matmul(charges1(iatom, 1:ni), gammab)
213 IF (dr > 1.e-6_dp)
THEN
214 CALL dgamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
216 gchrg0(iatom, 1:ni, i + 1) = gchrg0(iatom, 1:ni, i + 1) &
217 + matmul(gammab, charges0(jatom, 1:nj))*rij(i)/dr
218 gchrg1(iatom, 1:ni, i + 1) = gchrg1(iatom, 1:ni, i + 1) &
219 + matmul(gammab, charges1(jatom, 1:nj))*rij(i)/dr
220 IF (iatom /= jatom)
THEN
221 gchrg0(jatom, 1:nj, i + 1) = gchrg0(jatom, 1:nj, i + 1) &
222 - matmul(charges0(iatom, 1:ni), gammab)*rij(i)/dr
223 gchrg1(jatom, 1:nj, i + 1) = gchrg1(jatom, 1:nj, i + 1) &
224 - matmul(charges1(iatom, 1:ni), gammab)*rij(i)/dr
234 IF (xtb_control%coulomb_lr)
THEN
235 do_ewald = xtb_control%do_ewald
238 NULLIFY (ewald_env, ewald_pw)
240 ewald_env=ewald_env, ewald_pw=ewald_pw)
241 CALL get_cell(cell=cell, periodic=periodic, deth=deth)
242 CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
243 CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
244 CALL tb_ewald_overlap(gmcharge0, mcharge0, alpha, n_list, virial, use_virial)
245 CALL tb_ewald_overlap(gmcharge1, mcharge1, alpha, n_list, virial, use_virial)
246 SELECT CASE (ewald_type)
248 cpabort(
"Invalid Ewald type")
250 cpabort(
"Not allowed with DFTB")
252 cpabort(
"Standard Ewald not implemented in DFTB")
254 cpabort(
"PME not implemented in DFTB")
256 CALL tb_spme_zforce(ewald_env, ewald_pw, particle_set, cell, gmcharge0, mcharge0)
257 CALL tb_spme_zforce(ewald_env, ewald_pw, particle_set, cell, gmcharge1, mcharge1)
261 CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
262 DO ikind = 1,
SIZE(local_particles%n_el)
263 DO ia = 1, local_particles%n_el(ikind)
264 iatom = local_particles%list(ikind)%array(ia)
265 DO jatom = 1, iatom - 1
266 rij = particle_set(iatom)%r - particle_set(jatom)%r
268 dr = sqrt(sum(rij(:)**2))
269 IF (dr > 1.e-6_dp)
THEN
270 gmcharge0(iatom, 1) = gmcharge0(iatom, 1) + mcharge0(jatom)/dr
271 gmcharge0(jatom, 1) = gmcharge0(jatom, 1) + mcharge0(iatom)/dr
272 gmcharge1(iatom, 1) = gmcharge1(iatom, 1) + mcharge1(jatom)/dr
273 gmcharge1(jatom, 1) = gmcharge1(jatom, 1) + mcharge1(iatom)/dr
275 gmcharge0(iatom, i) = gmcharge0(iatom, i) + rij(i - 1)*mcharge0(jatom)/dr**3
276 gmcharge0(jatom, i) = gmcharge0(jatom, i) - rij(i - 1)*mcharge0(iatom)/dr**3
277 gmcharge1(iatom, i) = gmcharge1(iatom, i) + rij(i - 1)*mcharge1(jatom)/dr**3
278 gmcharge1(jatom, i) = gmcharge1(jatom, i) - rij(i - 1)*mcharge1(iatom)/dr**3
284 cpassert(.NOT. use_virial)
290 atomic_kind_set=atomic_kind_set, &
291 force=force, para_env=para_env)
292 CALL para_env%sum(gmcharge0(:, 1))
293 CALL para_env%sum(gchrg0(:, :, 1))
294 CALL para_env%sum(gmcharge1(:, 1))
295 CALL para_env%sum(gchrg1(:, :, 1))
297 IF (xtb_control%coulomb_lr)
THEN
300 gmcharge0(:, 1) = gmcharge0(:, 1) - 2._dp*alpha*
oorootpi*mcharge0(:)
301 IF (any(periodic(:) == 1))
THEN
302 gmcharge0(:, 1) = gmcharge0(:, 1) -
pi/alpha**2/deth
304 gmcharge1(:, 1) = gmcharge1(:, 1) - 2._dp*alpha*
oorootpi*mcharge1(:)
305 IF (any(periodic(:) == 1))
THEN
306 gmcharge1(:, 1) = gmcharge1(:, 1) -
pi/alpha**2/deth
313 atom_of_kind=atom_of_kind)
315 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
317 ikind = kind_of(iatom)
318 atom_i = atom_of_kind(iatom)
319 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
325 fij(i) = sum(charges0(iatom, 1:ni)*gchrg1(iatom, 1:ni, i + 1)) + &
326 sum(charges1(iatom, 1:ni)*gchrg0(iatom, 1:ni, i + 1))
328 force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
329 force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
330 force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
334 fij(i) = gmcharge1(iatom, i + 1)*mcharge0(iatom) + &
335 gmcharge0(iatom, i + 1)*mcharge1(iatom)
337 force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
338 force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
339 force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
341 IF (debug_forces)
THEN
342 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
343 CALL para_env%sum(fodeb)
344 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: P*dH[Pz] ", fodeb
347 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
349 IF (
SIZE(matrix_p0) == 2)
THEN
350 CALL dbcsr_add(matrix_p0(1)%matrix, matrix_p0(2)%matrix, &
351 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
352 CALL dbcsr_add(matrix_p1(1)%matrix, matrix_p1(2)%matrix, &
353 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
357 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
358 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
359 DO WHILE (dbcsr_iterator_blocks_left(iter))
360 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
361 ikind = kind_of(irow)
362 jkind = kind_of(icol)
365 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
366 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
372 ALLOCATE (gcij0(ni, nj))
373 ALLOCATE (gcij1(ni, nj))
378 gcij0(i, j) = 0.5_dp*(gchrg0(irow, la, 1) + gchrg0(icol, lb, 1))
379 gcij1(i, j) = 0.5_dp*(gchrg1(irow, la, 1) + gchrg1(icol, lb, 1))
382 gmij0 = 0.5_dp*(gmcharge0(irow, 1) + gmcharge0(icol, 1))
383 gmij1 = 0.5_dp*(gmcharge1(irow, 1) + gmcharge1(icol, 1))
384 atom_i = atom_of_kind(irow)
385 atom_j = atom_of_kind(icol)
387 CALL dbcsr_get_block_p(matrix=matrix_p0(1)%matrix, &
388 row=irow, col=icol, block=pblock0, found=found)
391 CALL dbcsr_get_block_p(matrix=matrix_p1(1)%matrix, &
392 row=irow, col=icol, block=pblock1, found=found)
396 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
397 row=irow, col=icol, block=dsblock, found=found)
400 fi = -2.0_dp*sum(pblock0*dsblock*gcij1) - 2.0_dp*sum(pblock1*dsblock*gcij0)
401 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
402 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
404 fi = -2.0_dp*gmij1*sum(pblock0*dsblock) - 2.0_dp*gmij0*sum(pblock1*dsblock)
405 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
406 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
408 DEALLOCATE (gcij0, gcij1)
410 CALL dbcsr_iterator_stop(iter)
411 IF (debug_forces)
THEN
412 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
413 CALL para_env%sum(fodeb)
414 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pz*H[P]*dS ", fodeb
417 IF (xtb_control%tb3_interaction)
THEN
419 ALLOCATE (xgamma(nkind))
421 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
425 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
426 CALL dftb3_diagonal_hessian_force(qs_env, mcharge0, mcharge1, &
427 matrix_p0(1)%matrix, matrix_p1(1)%matrix, xgamma)
428 IF (debug_forces)
THEN
429 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
430 CALL para_env%sum(fodeb)
431 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pz*H3[P] ", fodeb
436 IF (
SIZE(matrix_p0) == 2)
THEN
437 CALL dbcsr_add(matrix_p0(1)%matrix, matrix_p0(2)%matrix, &
438 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
439 CALL dbcsr_add(matrix_p1(1)%matrix, matrix_p1(2)%matrix, &
440 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
444 IF (qs_env%qmmm .AND. qs_env%qmmm_periodic)
THEN
445 cpabort(
"Not Available")
448 DEALLOCATE (gmcharge0, gchrg0, gmcharge1, gchrg1)
450 CALL timestop(handle)
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.
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.