76 sap_int, calculate_forces, just_energy)
79 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix, matrix_p
82 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
84 CHARACTER(len=*),
PARAMETER :: routinen =
'build_xtb_spinpol'
86 INTEGER :: atom_a, atom_i, atom_j, handle, i, ia, iac, iatom, ib, ic, icol, ikind, iknd, &
87 irow, jatom, jkind, jknd, la, lb, na, natom, natorb, nb, nimg, nkind, nsgf, nspins
88 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
89 INTEGER,
DIMENSION(25) :: lao
90 INTEGER,
DIMENSION(3) :: cellind
91 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
92 LOGICAL :: defined, found, use_virial
93 REAL(kind=
dp) :: dr, espin, fi, fval
94 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: docg
95 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: aocg, bocg, pam, pbm, wab
96 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: wabk
97 REAL(kind=
dp),
DIMENSION(3) :: fij, rij
98 REAL(kind=
dp),
DIMENSION(3, 3) :: wall
99 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: aksb, bksb, dsblock, pamat, pbmat, sblock
100 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: dsint
106 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p_kp, matrix_s, matrix_s_kp
112 DIMENSION(:),
POINTER :: nl_iterator
117 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
121 CALL timeset(routinen, handle)
123 energy%xtb_spinpol = 0.0_dp
125 CALL get_qs_env(qs_env, dft_control=dft_control)
126 nspins = dft_control%nspins
127 nimg = dft_control%nimages
129 IF (nspins == 2)
THEN
134 qs_kind_set=qs_kind_set, &
135 particle_set=particle_set, &
136 atomic_kind_set=atomic_kind_set, &
143 atom_of_kind=atom_of_kind)
146 IF (calculate_forces)
THEN
147 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
150 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
152 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env)
155 ALLOCATE (wabk(nsgf, nsgf, nkind))
158 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
164 wabk(ia, ib, ikind) = wall(la, lb)
170 ALLOCATE (aocg(nsgf, natom), bocg(nsgf, natom))
174 matrix_s_kp => matrix_s(:, :)
175 matrix_p_kp => matrix_p(1:1, :)
176 CALL ao_charges(matrix_p_kp, matrix_s_kp, aocg, para_env)
177 matrix_p_kp => matrix_p(2:2, :)
178 CALL ao_charges(matrix_p_kp, matrix_s_kp, bocg, para_env)
180 s_matrix => matrix_s(1, 1)%matrix
181 p_matrix => matrix_p(1:1, 1)
182 CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
183 p_matrix => matrix_p(2:2, 1)
184 CALL ao_charges(p_matrix, s_matrix, bocg, para_env)
190 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
192 IF (.NOT. defined .OR. natorb < 1) cycle
193 ALLOCATE (docg(natorb), wab(natorb, natorb))
194 wab(1:natorb, 1:natorb) = wabk(1:natorb, 1:natorb, ikind)
196 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
198 docg(1:natorb) = aocg(1:natorb, atom_a) - bocg(1:natorb, atom_a)
199 espin = 0.5_dp*dot_product(docg, matmul(wab, docg))
200 energy%xtb_spinpol = energy%xtb_spinpol + espin
201 IF (atprop%energy)
THEN
202 atprop%atecoul(iatom) = atprop%atecoul(iatom) + espin
205 DEALLOCATE (docg, wab)
209 IF (calculate_forces)
THEN
211 NULLIFY (cell_to_index)
214 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
222 ikind = kind_of(irow)
223 atom_i = atom_of_kind(irow)
224 jkind = kind_of(icol)
225 atom_j = atom_of_kind(icol)
228 row=irow, col=icol, block=pamat, found=found)
231 row=irow, col=icol, block=pbmat, found=found)
239 row=irow, col=icol, block=dsblock, found=found)
242 CALL fupdate(fi, pamat, pbmat, dsblock, na, nb, &
243 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
244 aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
246 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
247 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
253 IF (use_virial .AND. 0 == 0)
THEN
254 cpassert(
ASSOCIATED(sap_int))
257 iac = ikind + nkind*(jkind - 1)
258 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist)) cycle
259 DO ia = 1, sap_int(iac)%nalist
260 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist(ia)%clist)) cycle
261 iatom = sap_int(iac)%alist(ia)%aatom
262 DO ic = 1, sap_int(iac)%alist(ia)%nclist
263 jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
264 rij = sap_int(iac)%alist(ia)%clist(ic)%rac
265 dr = sqrt(sum(rij(:)**2))
266 IF (dr > 1.e-6_dp)
THEN
267 dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
268 icol = max(iatom, jatom)
269 irow = min(iatom, jatom)
271 row=irow, col=icol, block=pamat, found=found)
274 row=irow, col=icol, block=pbmat, found=found)
276 IF (irow == iatom)
THEN
279 ALLOCATE (pam(na, nb), pbm(na, nb))
280 pam(1:na, 1:nb) = pamat(1:na, 1:nb)
281 pbm(1:na, 1:nb) = pbmat(1:na, 1:nb)
285 ALLOCATE (pam(na, nb), pbm(na, nb))
286 pam(1:na, 1:nb) = transpose(pamat(1:nb, 1:na))
287 pbm(1:na, 1:nb) = transpose(pbmat(1:nb, 1:na))
291 CALL fupdate(fi, pam, pbm, dsint(:, :, i), na, nb, &
292 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
293 aocg(1:na, iatom), aocg(1:nb, jatom), &
294 bocg(1:na, iatom), bocg(1:nb, jatom))
298 IF (iatom == jatom) fi = 0.5_dp
300 DEALLOCATE (pam, pbm)
310 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
314 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
316 dr = sqrt(sum(rij**2))
317 IF (iatom == jatom .AND. dr < 1.0e-6_dp) cycle
319 icol = max(iatom, jatom)
320 irow = min(iatom, jatom)
322 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
325 IF (irow == iatom)
THEN
328 atom_i = atom_of_kind(iatom)
329 atom_j = atom_of_kind(jatom)
334 atom_i = atom_of_kind(jatom)
335 atom_j = atom_of_kind(iatom)
340 row=irow, col=icol, block=pamat, found=found)
343 row=irow, col=icol, block=pbmat, found=found)
352 row=irow, col=icol, block=dsblock, found=found)
355 CALL fupdate(fi, pamat, pbmat, dsblock, na, nb, &
356 wabk(1:na, 1:na, iknd), wabk(1:nb, 1:nb, jknd), &
357 aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
359 force(iknd)%rho_elec(i, atom_i) = force(iknd)%rho_elec(i, atom_i) + fi
360 force(jknd)%rho_elec(i, atom_j) = force(jknd)%rho_elec(i, atom_j) - fi
365 IF (iatom == jatom) fi = 0.5_dp
376 IF (.NOT. just_energy)
THEN
378 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
387 row=irow, col=icol, block=aksb, found=found)
390 row=irow, col=icol, block=bksb, found=found)
394 ikind = kind_of(irow)
395 jkind = kind_of(icol)
397 CALL ksupdate(aksb, bksb, sblock, na, nb, fval, &
398 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
399 aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
403 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
407 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
409 icol = max(iatom, jatom)
410 irow = min(iatom, jatom)
412 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
415 ikind = kind_of(irow)
416 jkind = kind_of(icol)
419 row=irow, col=icol, block=sblock, found=found)
422 row=irow, col=icol, block=aksb, found=found)
425 row=irow, col=icol, block=bksb, found=found)
431 CALL ksupdate(aksb, bksb, sblock, na, nb, fval, &
432 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
433 aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
441 DEALLOCATE (aocg, bocg)
444 CALL timestop(handle)
457 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_matrix, matrix_p1
459 CHARACTER(len=*),
PARAMETER :: routinen =
'xtb_spinpol_hessian'
461 INTEGER :: handle, ia, ib, icol, ikind, irow, &
462 jkind, la, lb, na, natom, natorb, nb, &
463 nimg, nkind, nsgf, nspins
464 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: kind_of
465 INTEGER,
DIMENSION(25) :: lao
467 REAL(kind=
dp) :: fval
468 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: aocg1, bocg1
469 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: wabk
470 REAL(kind=
dp),
DIMENSION(3, 3) :: wall
471 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: aksb, bksb, sblock
475 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_s
479 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
482 CALL timeset(routinen, handle)
484 CALL get_qs_env(qs_env, dft_control=dft_control)
485 nspins = dft_control%nspins
486 nimg = dft_control%nimages
489 cpabort(
"No kpoints allowed in xTB response calculation")
492 IF (nspins == 2)
THEN
495 qs_kind_set=qs_kind_set, &
496 atomic_kind_set=atomic_kind_set)
500 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
502 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env)
505 ALLOCATE (wabk(nsgf, nsgf, nkind))
508 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
514 wabk(ia, ib, ikind) = wall(la, lb)
520 ALLOCATE (aocg1(nsgf, natom), bocg1(nsgf, natom))
523 s_matrix => matrix_s(1, 1)%matrix
524 p_matrix => matrix_p1(1:1)
525 CALL ao_charges(p_matrix, s_matrix, aocg1, para_env)
526 p_matrix => matrix_p1(2:2)
527 CALL ao_charges(p_matrix, s_matrix, bocg1, para_env)
535 row=irow, col=icol, block=aksb, found=found)
538 row=irow, col=icol, block=bksb, found=found)
542 ikind = kind_of(irow)
543 jkind = kind_of(icol)
545 CALL ksupdate(aksb, bksb, sblock, na, nb, fval, &
546 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
547 aocg1(1:na, irow), aocg1(1:nb, icol), &
548 bocg1(1:na, irow), bocg1(1:nb, icol))
553 DEALLOCATE (aocg1, bocg1)
556 CALL timestop(handle)
568 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_p0, matrix_p1
570 CHARACTER(len=*),
PARAMETER :: routinen =
'xtb_spinpol_hforce'
572 INTEGER :: atom_i, atom_j, handle, i, ia, ib, icol, &
573 ikind, irow, jkind, la, lb, na, natom, &
574 natorb, nb, nimg, nkind, nsgf, nspins
575 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
576 INTEGER,
DIMENSION(25) :: lao
579 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: aocg, aocg1, bocg, bocg1
580 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: wabk
581 REAL(kind=
dp),
DIMENSION(3, 3) :: wall
582 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: dsblock, p0amat, p0bmat, p1amat, p1bmat, &
586 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, p_matrix
592 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
595 CALL timeset(routinen, handle)
597 CALL get_qs_env(qs_env, dft_control=dft_control)
598 nspins = dft_control%nspins
599 nimg = dft_control%nimages
601 cpabort(
"xTB response forces for spin polarisation Hamiltonian not available")
604 IF (nspins == 2)
THEN
607 qs_kind_set=qs_kind_set, &
608 particle_set=particle_set, &
609 atomic_kind_set=atomic_kind_set)
613 atom_of_kind=atom_of_kind)
615 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
618 CALL get_qs_env(qs_env, matrix_s=matrix_s, para_env=para_env)
621 ALLOCATE (wabk(nsgf, nsgf, nkind))
624 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
630 wabk(ia, ib, ikind) = wall(la, lb)
636 s_matrix => matrix_s(1)%matrix
637 ALLOCATE (aocg(nsgf, natom), bocg(nsgf, natom))
640 p_matrix => matrix_p0(1:1)
641 CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
642 p_matrix => matrix_p0(2:2)
643 CALL ao_charges(p_matrix, s_matrix, bocg, para_env)
645 ALLOCATE (aocg1(nsgf, natom), bocg1(nsgf, natom))
648 p_matrix => matrix_p1(1:1)
649 CALL ao_charges(p_matrix, s_matrix, aocg1, para_env)
650 p_matrix => matrix_p1(2:2)
651 CALL ao_charges(p_matrix, s_matrix, bocg1, para_env)
661 ikind = kind_of(irow)
662 atom_i = atom_of_kind(irow)
663 jkind = kind_of(icol)
664 atom_j = atom_of_kind(icol)
667 row=irow, col=icol, block=p0amat, found=found)
670 row=irow, col=icol, block=p0bmat, found=found)
673 row=irow, col=icol, block=p1amat, found=found)
676 row=irow, col=icol, block=p1bmat, found=found)
684 row=irow, col=icol, block=dsblock, found=found)
688 CALL f2update(fi, p0amat, p0bmat, p1amat, p1bmat, dsblock, na, nb, &
689 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
690 aocg(1:na, irow), aocg(1:nb, icol), &
691 bocg(1:na, irow), bocg(1:nb, icol), &
692 aocg1(1:na, irow), aocg1(1:nb, icol), &
693 bocg1(1:na, irow), bocg1(1:nb, icol))
695 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
696 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
704 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, mimic, 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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
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, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, 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_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, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.