27 USE dbcsr_api,
ONLY: dbcsr_add,&
29 dbcsr_iterator_blocks_left,&
30 dbcsr_iterator_next_block,&
31 dbcsr_iterator_start,&
37 ewald_environment_type
66 neighbor_list_iterator_p_type,&
68 neighbor_list_set_p_type
78 #include "./base/base_uses.f90"
84 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'xtb_coulomb'
102 calculate_forces, just_energy)
104 TYPE(qs_environment_type),
POINTER :: qs_env
105 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
106 TYPE(qs_rho_type),
POINTER :: rho
107 REAL(
dp),
DIMENSION(:, :),
INTENT(in) :: charges
108 REAL(
dp),
DIMENSION(:),
INTENT(in) :: mcharge
109 TYPE(qs_energy_type),
POINTER :: energy
110 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
112 CHARACTER(len=*),
PARAMETER :: routinen =
'build_xtb_coulomb'
114 INTEGER :: atom_i, atom_j, blk, ewald_type, handle, i, ia, iac, iatom, ic, icol, ikind, img, &
115 irow, is, j, jatom, jkind, la, lb, lmaxa, lmaxb, natom, natorb_a, natorb_b, ni, nimg, nj, &
117 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
118 INTEGER,
DIMENSION(25) :: laoa, laob
119 INTEGER,
DIMENSION(3) :: cellind, periodic
120 INTEGER,
DIMENSION(5) :: occ
121 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
122 LOGICAL :: defined, do_ewald, do_gamma_stress, &
124 REAL(kind=
dp) :: alpha, deth, dr, ecsr, etaa, etab, f1, &
125 f2, fi, gmij, kg, rcut, rcuta, rcutb, &
127 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: xgamma, zeffk
128 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: gammab, gcij, gmcharge
129 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: gchrg
130 REAL(kind=
dp),
DIMENSION(25) :: gcint
131 REAL(kind=
dp),
DIMENSION(3) :: fij, rij
132 REAL(kind=
dp),
DIMENSION(5) :: kappaa, kappab
133 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: dsblock, ksblock, pblock, sblock
134 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: dsint
135 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
136 TYPE(atprop_type),
POINTER :: atprop
137 TYPE(cell_type),
POINTER :: cell
138 TYPE(dbcsr_iterator_type) :: iter
139 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s
140 TYPE(dft_control_type),
POINTER :: dft_control
141 TYPE(distribution_1d_type),
POINTER :: local_particles
142 TYPE(ewald_environment_type),
POINTER :: ewald_env
143 TYPE(ewald_pw_type),
POINTER :: ewald_pw
144 TYPE(kpoint_type),
POINTER :: kpoints
145 TYPE(mp_para_env_type),
POINTER :: para_env
146 TYPE(neighbor_list_iterator_p_type), &
147 DIMENSION(:),
POINTER :: nl_iterator
148 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
150 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
151 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
152 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
153 TYPE(sap_int_type),
DIMENSION(:),
POINTER :: sap_int
154 TYPE(virial_type),
POINTER :: virial
155 TYPE(xtb_atom_type),
POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind
156 TYPE(xtb_control_type),
POINTER :: xtb_control
158 CALL timeset(routinen, handle)
160 NULLIFY (matrix_p, matrix_s, virial, atprop, dft_control)
163 qs_kind_set=qs_kind_set, &
164 particle_set=particle_set, &
168 dft_control=dft_control)
170 xtb_control => dft_control%qs_control%xtb_control
173 IF (calculate_forces)
THEN
174 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
177 do_gamma_stress = .false.
178 IF (.NOT. just_energy .AND. use_virial)
THEN
179 IF (dft_control%nimages == 1) do_gamma_stress = .true.
182 IF (atprop%energy)
THEN
183 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
184 natom =
SIZE(particle_set)
188 IF (calculate_forces)
THEN
194 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
195 ALLOCATE (gchrg(natom, 5, nmat))
197 ALLOCATE (gmcharge(natom, nmat))
204 IF (xtb_control%old_coulomb_damping)
THEN
205 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
207 CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
212 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
213 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
215 IF (.NOT. defined .OR. natorb_a < 1) cycle
216 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
218 IF (.NOT. defined .OR. natorb_b < 1) cycle
225 ALLOCATE (gammab(ni, nj))
227 dr = sqrt(sum(rij(:)**2))
228 CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
229 gchrg(iatom, 1:ni, 1) = gchrg(iatom, 1:ni, 1) + matmul(gammab, charges(jatom, 1:nj))
230 IF (iatom /= jatom)
THEN
231 gchrg(jatom, 1:nj, 1) = gchrg(jatom, 1:nj, 1) + matmul(charges(iatom, 1:ni), gammab)
233 IF (calculate_forces)
THEN
234 IF (dr > 1.e-6_dp)
THEN
235 CALL dgamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
237 gchrg(iatom, 1:ni, i + 1) = gchrg(iatom, 1:ni, i + 1) &
238 + matmul(gammab, charges(jatom, 1:nj))*rij(i)/dr
239 IF (iatom /= jatom)
THEN
240 gchrg(jatom, 1:nj, i + 1) = gchrg(jatom, 1:nj, i + 1) &
241 - matmul(charges(iatom, 1:ni), gammab)*rij(i)/dr
245 gcint(1:ni) = matmul(gammab, charges(jatom, 1:nj))
247 fij(i) = -sum(charges(iatom, 1:ni)*gcint(1:ni))*rij(i)/dr
250 IF (iatom == jatom) fi = 0.5_dp
252 IF (atprop%stress)
THEN
265 IF (xtb_control%coulomb_lr)
THEN
266 do_ewald = xtb_control%do_ewald
269 NULLIFY (ewald_env, ewald_pw)
271 ewald_env=ewald_env, ewald_pw=ewald_pw)
272 CALL get_cell(cell=cell, periodic=periodic, deth=deth)
273 CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
274 CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
276 use_virial, atprop=atprop)
277 SELECT CASE (ewald_type)
279 cpabort(
"Invalid Ewald type")
281 cpabort(
"Not allowed with DFTB")
283 cpabort(
"Standard Ewald not implemented in DFTB")
285 cpabort(
"PME not implemented in DFTB")
288 gmcharge, mcharge, calculate_forces, virial, &
289 use_virial, atprop=atprop)
294 local_particles=local_particles)
295 DO ikind = 1,
SIZE(local_particles%n_el)
296 DO ia = 1, local_particles%n_el(ikind)
297 iatom = local_particles%list(ikind)%array(ia)
298 DO jatom = 1, iatom - 1
299 rij = particle_set(iatom)%r - particle_set(jatom)%r
301 dr = sqrt(sum(rij(:)**2))
302 IF (dr > 1.e-6_dp)
THEN
303 gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)/dr
304 gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)/dr
306 gmcharge(iatom, i) = gmcharge(iatom, i) + rij(i - 1)*mcharge(jatom)/dr**3
307 gmcharge(jatom, i) = gmcharge(jatom, i) - rij(i - 1)*mcharge(iatom)/dr**3
313 cpassert(.NOT. use_virial)
319 atomic_kind_set=atomic_kind_set, &
320 force=force, para_env=para_env)
321 CALL para_env%sum(gmcharge(:, 1))
322 CALL para_env%sum(gchrg(:, :, 1))
324 IF (xtb_control%coulomb_lr)
THEN
327 gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*
oorootpi*mcharge(:)
328 IF (any(periodic(:) == 1))
THEN
329 gmcharge(:, 1) = gmcharge(:, 1) -
pi/alpha**2/deth
337 atom_of_kind=atom_of_kind)
340 ikind = kind_of(iatom)
341 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
344 ecsr = ecsr + sum(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, 1))
347 energy%hartree = energy%hartree + 0.5_dp*ecsr
348 energy%hartree = energy%hartree + 0.5_dp*sum(mcharge(:)*gmcharge(:, 1))
350 IF (atprop%energy)
THEN
351 CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
352 DO ikind = 1,
SIZE(local_particles%n_el)
353 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
356 zeff = sum(real(occ, kind=
dp))
357 DO ia = 1, local_particles%n_el(ikind)
358 iatom = local_particles%list(ikind)%array(ia)
359 atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
360 0.5_dp*sum(real(occ(1:ni), kind=
dp)*gchrg(iatom, 1:ni, 1))
361 atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
362 0.5_dp*zeff*gmcharge(iatom, 1)
367 IF (calculate_forces)
THEN
369 ikind = kind_of(iatom)
370 atom_i = atom_of_kind(iatom)
371 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
376 fij(i) = sum(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, i + 1))
378 force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
379 force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
380 force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
383 fij(i) = gmcharge(iatom, i + 1)*mcharge(iatom)
385 force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
386 force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
387 force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
391 IF (.NOT. just_energy)
THEN
392 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
395 nimg = dft_control%nimages
396 NULLIFY (cell_to_index)
399 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
403 IF (calculate_forces .AND.
SIZE(matrix_p, 1) == 2)
THEN
405 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
406 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
411 IF (do_gamma_stress)
THEN
418 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
419 DO WHILE (dbcsr_iterator_blocks_left(iter))
420 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
421 ikind = kind_of(irow)
422 jkind = kind_of(icol)
425 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
426 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
432 ALLOCATE (gcij(ni, nj))
437 gcij(i, j) = 0.5_dp*(gchrg(irow, la, 1) + gchrg(icol, lb, 1))
440 gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1))
441 DO is = 1,
SIZE(ks_matrix, 1)
443 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
444 row=irow, col=icol, block=ksblock, found=found)
446 ksblock = ksblock - gcij*sblock
447 ksblock = ksblock - gmij*sblock
449 IF (calculate_forces)
THEN
450 atom_i = atom_of_kind(irow)
451 atom_j = atom_of_kind(icol)
453 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
454 row=irow, col=icol, block=pblock, found=found)
458 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
459 row=irow, col=icol, block=dsblock, found=found)
463 fi = -2.0_dp*sum(pblock*dsblock*gcij)
464 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
465 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
468 fi = -2.0_dp*gmij*sum(pblock*dsblock)
469 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
470 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
476 CALL dbcsr_iterator_stop(iter)
478 IF (do_gamma_stress)
THEN
481 iac = ikind + nkind*(jkind - 1)
482 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist)) cycle
484 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
485 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
488 DO ia = 1, sap_int(iac)%nalist
489 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist(ia)%clist)) cycle
490 iatom = sap_int(iac)%alist(ia)%aatom
491 DO ic = 1, sap_int(iac)%alist(ia)%nclist
492 jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
493 rij = sap_int(iac)%alist(ia)%clist(ic)%rac
494 dr = sqrt(sum(rij(:)**2))
495 IF (dr > 1.e-6_dp)
THEN
496 dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
497 ALLOCATE (gcij(ni, nj))
502 gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
505 gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
506 icol = max(iatom, jatom)
507 irow = min(iatom, jatom)
509 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
510 row=irow, col=icol, block=pblock, found=found)
515 IF (irow == iatom)
THEN
516 f1 = -2.0_dp*sum(pblock*dsint(:, :, i)*gcij)
517 f2 = -2.0_dp*gmij*sum(pblock*dsint(:, :, i))
519 f1 = -2.0_dp*sum(transpose(pblock)*dsint(:, :, i)*gcij)
520 f2 = -2.0_dp*gmij*sum(transpose(pblock)*dsint(:, :, i))
526 IF (iatom == jatom) fi = 0.5_dp
528 IF (atprop%stress)
THEN
540 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
544 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
546 icol = max(iatom, jatom)
547 irow = min(iatom, jatom)
549 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
553 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
554 row=irow, col=icol, block=sblock, found=found)
558 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
559 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
565 ALLOCATE (gcij(ni, nj))
568 IF (irow == iatom)
THEN
571 gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
575 gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
579 gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
580 DO is = 1,
SIZE(ks_matrix, 1)
582 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
583 row=irow, col=icol, block=ksblock, found=found)
585 ksblock = ksblock - gcij*sblock
586 ksblock = ksblock - gmij*sblock
589 IF (calculate_forces)
THEN
590 atom_i = atom_of_kind(iatom)
591 atom_j = atom_of_kind(jatom)
592 IF (irow /= iatom)
THEN
597 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
598 row=irow, col=icol, block=pblock, found=found)
602 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
603 row=irow, col=icol, block=dsblock, found=found)
607 fi = -2.0_dp*sum(pblock*dsblock*gcij)
608 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
609 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
612 fi = -2.0_dp*gmij*sum(pblock*dsblock)
613 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
614 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
618 dr = sqrt(sum(rij(:)**2))
619 IF (dr > 1.e-6_dp)
THEN
621 IF (iatom == jatom) fi = 0.5_dp
623 IF (atprop%stress)
THEN
636 IF (calculate_forces .AND.
SIZE(matrix_p, 1) == 2)
THEN
638 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
639 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
644 IF (xtb_control%tb3_interaction)
THEN
646 ALLOCATE (zeffk(nkind), xgamma(nkind))
648 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
653 sap_int, calculate_forces, just_energy)
654 DEALLOCATE (zeffk, xgamma)
658 IF (qs_env%qmmm .AND. qs_env%qmmm_periodic)
THEN
660 calculate_forces, just_energy)
663 IF (do_gamma_stress)
THEN
667 CALL timestop(handle)
692 SUBROUTINE gamma_rab_sr(gmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
693 REAL(
dp),
DIMENSION(:, :),
INTENT(INOUT) :: gmat
694 REAL(
dp),
INTENT(IN) :: rab
695 INTEGER,
INTENT(IN) :: nla
696 REAL(
dp),
DIMENSION(:),
INTENT(IN) :: kappaa
697 REAL(
dp),
INTENT(IN) :: etaa
698 INTEGER,
INTENT(IN) :: nlb
699 REAL(
dp),
DIMENSION(:),
INTENT(IN) :: kappab
700 REAL(
dp),
INTENT(IN) :: etab, kg, rcut
702 REAL(kind=
dp),
PARAMETER :: rsmooth = 1.0_dp
705 REAL(kind=
dp) :: fcut, r, rk, x
706 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: eta
708 ALLOCATE (eta(nla, nlb))
713 eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
714 eta(i, j) = 2._dp/eta(i, j)
719 IF (rab < 1.e-6_dp)
THEN
721 gmat(:, :) = eta(:, :)
722 ELSEIF (rab > rcut)
THEN
727 IF (rab < rcut - rsmooth)
THEN
730 r = rab - (rcut - rsmooth)
732 fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
734 gmat(:, :) = fcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg) - fcut/rab
762 SUBROUTINE dgamma_rab_sr(dgmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
763 REAL(
dp),
DIMENSION(:, :),
INTENT(INOUT) :: dgmat
764 REAL(
dp),
INTENT(IN) :: rab
765 INTEGER,
INTENT(IN) :: nla
766 REAL(
dp),
DIMENSION(:),
INTENT(IN) :: kappaa
767 REAL(
dp),
INTENT(IN) :: etaa
768 INTEGER,
INTENT(IN) :: nlb
769 REAL(
dp),
DIMENSION(:),
INTENT(IN) :: kappab
770 REAL(
dp),
INTENT(IN) :: etab, kg, rcut
772 REAL(kind=
dp),
PARAMETER :: rsmooth = 1.0_dp
775 REAL(kind=
dp) :: dfcut, fcut, r, rk, x
776 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: eta
778 ALLOCATE (eta(nla, nlb))
782 eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
783 eta(i, j) = 2._dp/eta(i, j)
787 IF (rab < 1.e-6)
THEN
790 ELSEIF (rab > rcut)
THEN
795 IF (rab < rcut - rsmooth)
THEN
799 r = rab - (rcut - rsmooth)
801 fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
802 dfcut = -30._dp*x**4 + 60._dp*x**3 - 30._dp*x**2
803 dfcut = dfcut/rsmooth
805 dgmat(:, :) = dfcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg)
806 dgmat(:, :) = dgmat(:, :) - dfcut/rab + fcut/rab**2
807 dgmat(:, :) = dgmat(:, :) - fcut/(rk + eta(:, :))*(1._dp/(rk + eta(:, :)))**(1._dp/kg)*rk/rab
821 TYPE(qs_environment_type),
POINTER :: qs_env
822 TYPE(sap_int_type),
DIMENSION(:),
POINTER :: sap_int
824 CHARACTER(LEN=*),
PARAMETER :: routinen =
'xtb_dsint_list'
826 INTEGER :: handle, i, iac, iatom, ikind, ilist, iset, jatom, jkind, jneighbor, jset, ldsab, &
827 n1, n2, natorb_a, natorb_b, ncoa, ncob, nkind, nlist, nneighbor, nseta, nsetb, sgfa, sgfb
828 INTEGER,
DIMENSION(3) :: cell
829 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
831 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa, first_sgfb
834 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: owork
835 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: oint, sint
836 REAL(kind=
dp),
DIMENSION(3) :: rij
837 REAL(kind=
dp),
DIMENSION(:),
POINTER :: set_radius_a, set_radius_b
838 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
839 TYPE(clist_type),
POINTER :: clist
840 TYPE(dft_control_type),
POINTER :: dft_control
841 TYPE(gto_basis_set_p_type),
DIMENSION(:),
POINTER :: basis_set_list
842 TYPE(gto_basis_set_type),
POINTER :: basis_set_a, basis_set_b
843 TYPE(neighbor_list_iterator_p_type), &
844 DIMENSION(:),
POINTER :: nl_iterator
845 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
847 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
848 TYPE(xtb_atom_type),
POINTER :: xtb_atom_a, xtb_atom_b
850 CALL timeset(routinen, handle)
853 cpassert(.NOT.
ASSOCIATED(sap_int))
854 ALLOCATE (sap_int(nkind*nkind))
855 DO i = 1, nkind*nkind
856 NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
857 sap_int(i)%nalist = 0
861 qs_kind_set=qs_kind_set, &
862 dft_control=dft_control, &
866 ALLOCATE (basis_set_list(nkind))
873 jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, &
874 inode=jneighbor, cell=cell, r=rij)
875 iac = ikind + nkind*(jkind - 1)
877 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
879 IF (.NOT. defined .OR. natorb_a < 1) cycle
880 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
882 IF (.NOT. defined .OR. natorb_b < 1) cycle
884 dr = sqrt(sum(rij(:)**2))
887 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist))
THEN
888 sap_int(iac)%a_kind = ikind
889 sap_int(iac)%p_kind = jkind
890 sap_int(iac)%nalist = nlist
891 ALLOCATE (sap_int(iac)%alist(nlist))
893 NULLIFY (sap_int(iac)%alist(i)%clist)
894 sap_int(iac)%alist(i)%aatom = 0
895 sap_int(iac)%alist(i)%nclist = 0
898 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist(ilist)%clist))
THEN
899 sap_int(iac)%alist(ilist)%aatom = iatom
900 sap_int(iac)%alist(ilist)%nclist = nneighbor
901 ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
903 sap_int(iac)%alist(ilist)%clist(i)%catom = 0
906 clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
910 ALLOCATE (clist%acint(natorb_a, natorb_b, 3))
911 NULLIFY (clist%achint)
914 NULLIFY (clist%sgf_list)
917 basis_set_a => basis_set_list(ikind)%gto_basis_set
918 IF (.NOT.
ASSOCIATED(basis_set_a)) cycle
919 basis_set_b => basis_set_list(jkind)%gto_basis_set
920 IF (.NOT.
ASSOCIATED(basis_set_b)) cycle
922 first_sgfa => basis_set_a%first_sgf
923 la_max => basis_set_a%lmax
924 la_min => basis_set_a%lmin
925 npgfa => basis_set_a%npgf
926 nseta = basis_set_a%nset
927 nsgfa => basis_set_a%nsgf_set
928 rpgfa => basis_set_a%pgf_radius
929 set_radius_a => basis_set_a%set_radius
930 scon_a => basis_set_a%scon
931 zeta => basis_set_a%zet
933 first_sgfb => basis_set_b%first_sgf
934 lb_max => basis_set_b%lmax
935 lb_min => basis_set_b%lmin
936 npgfb => basis_set_b%npgf
937 nsetb = basis_set_b%nset
938 nsgfb => basis_set_b%nsgf_set
939 rpgfb => basis_set_b%pgf_radius
940 set_radius_b => basis_set_b%set_radius
941 scon_b => basis_set_b%scon
942 zetb => basis_set_b%zet
944 ldsab = get_memory_usage(qs_kind_set,
"ORB",
"ORB")
945 ALLOCATE (oint(ldsab, ldsab, 4), owork(ldsab, ldsab))
946 ALLOCATE (sint(natorb_a, natorb_b, 4))
950 ncoa = npgfa(iset)*
ncoset(la_max(iset))
951 n1 = npgfa(iset)*(
ncoset(la_max(iset)) -
ncoset(la_min(iset) - 1))
952 sgfa = first_sgfa(1, iset)
954 IF (set_radius_a(iset) + set_radius_b(jset) < dr) cycle
955 ncob = npgfb(jset)*
ncoset(lb_max(jset))
956 n2 = npgfb(jset)*(
ncoset(lb_max(jset)) -
ncoset(lb_min(jset) - 1))
957 sgfb = first_sgfb(1, jset)
958 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
959 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
960 rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
963 CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
964 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
965 CALL block_add(
"IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), &
966 sgfa, sgfb, trans=.false.)
971 clist%acint(1:natorb_a, 1:natorb_b, 1:3) = sint(1:natorb_a, 1:natorb_b, 2:4)
973 DEALLOCATE (oint, owork, sint)
978 DEALLOCATE (basis_set_list)
980 CALL timestop(handle)
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
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.
Holds information on atomic properties.
subroutine, public atprop_array_init(atarray, natom)
...
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
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
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.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
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
Calculation of QMMM Coulomb contributions in TB.
subroutine, public build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
Calculation of DFTB3 Terms.
subroutine, public build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeff, sap_int, calculate_forces, just_energy)
...
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.
Some utility functions for the calculation of integrals.
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
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)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
General overlap type integrals containers.
subroutine, public release_sap_int(sap_int)
...
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
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: ...
subroutine, public dgamma_rab_sr(dgmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
Computes the derivative of the short-range gamma parameter from Nataga-Mishimoto-Ohno-Klopman formula...
subroutine, public build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, calculate_forces, just_energy)
...
subroutine, public xtb_dsint_list(qs_env, sap_int)
...
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)
...