24 dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_symmetric
79#include "./base/base_uses.f90"
83 INTEGER,
DIMENSION(16),
PARAMETER ::
orbptr = [0, 1, 1, 1, &
84 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3]
88 CHARACTER(len=*),
PARAMETER,
PRIVATE :: modulen =
'qs_dftb_matrices'
89 REAL(kind=
dp),
PARAMETER,
PRIVATE :: dftb_fd_deriv_step = 1.0e-3_dp
105 LOGICAL,
INTENT(IN) :: calculate_forces
107 CHARACTER(LEN=*),
PARAMETER :: routinen =
'build_dftb_matrices'
109 INTEGER :: after, atom_a, atom_b, handle, i, iatom, ic, icol, ikind, img, irow, iw, jatom, &
110 jkind, l1, l2, la, lb, llm, lmaxi, lmaxj, m, n1, n2, n_urpoly, natorb_a, natorb_b, &
111 nderivatives, ngrd, ngrdcut, nimg, nkind, spdim
112 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind
113 INTEGER,
DIMENSION(3) :: cell
114 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
115 LOGICAL :: defined, found, omit_headers, use_virial
116 REAL(kind=
dp) :: ddr, dgrd, dr, erep, erepij, f0, foab, &
118 REAL(kind=
dp),
DIMENSION(0:3) :: eta_a, eta_b, skself
119 REAL(kind=
dp),
DIMENSION(10) :: urep
120 REAL(kind=
dp),
DIMENSION(2) :: surr
121 REAL(kind=
dp),
DIMENSION(3) :: drij, force_ab, force_rr, force_w, rij, &
123 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: dfblock, dsblock, fblock, fmatij, &
124 fmatji, pblock, sblock, scoeff, &
125 smatij, smatji, spxr, wblock
130 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
135 DIMENSION(:),
POINTER :: nl_iterator
141 POINTER :: dftb_potential
145 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
150 CALL timeset(routinen, handle)
157 DO l1 = 0, max(la, lb)
158 DO l2 = 0, min(l1, la, lb)
161 iptr(l1, l2, m, la, lb) = llm
168 NULLIFY (logger, virial, atprop)
171 NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
172 qs_kind_set, sab_orb, ks_env)
176 atomic_kind_set=atomic_kind_set, &
177 qs_kind_set=qs_kind_set, &
178 matrix_h_kp=matrix_h, &
179 matrix_s_kp=matrix_s, &
181 dft_control=dft_control, &
184 dftb_control => dft_control%qs_control%dftb_control
185 nimg = dft_control%nimages
187 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
189 IF (dftb_control%self_consistent .AND. calculate_forces) nderivatives = 1
190 CALL setup_matrices2(qs_env, nderivatives, nimg, matrix_s,
"OVERLAP", sab_orb)
191 CALL setup_matrices2(qs_env, 0, nimg, matrix_h,
"CORE HAMILTONIAN", sab_orb)
195 NULLIFY (dftb_potential)
196 CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential)
197 NULLIFY (particle_set)
198 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
200 IF (calculate_forces)
THEN
201 NULLIFY (rho, force, matrix_w)
204 matrix_w_kp=matrix_w, &
209 IF (
SIZE(matrix_p, 1) == 2)
THEN
211 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
212 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
213 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
214 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
218 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
221 IF (atprop%energy)
THEN
225 NULLIFY (cell_to_index)
227 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
233 nkind =
SIZE(atomic_kind_set)
238 iatom=iatom, jatom=jatom, r=rij, cell=cell)
239 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
241 defined=defined, lmax=lmaxi, skself=skself, &
242 eta=eta_a, natorb=natorb_a)
243 IF (.NOT. defined .OR. natorb_a < 1) cycle
244 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
246 defined=defined, lmax=lmaxj, eta=eta_b, natorb=natorb_b)
248 IF (.NOT. defined .OR. natorb_b < 1) cycle
251 dftb_param_ij => dftb_potential(ikind, jkind)
252 dftb_param_ji => dftb_potential(jkind, ikind)
254 ngrd = dftb_param_ij%ngrd
255 ngrdcut = dftb_param_ij%ngrdcut
256 dgrd = dftb_param_ij%dgrd
257 ddr = dgrd*dftb_fd_deriv_step
258 cpassert(dftb_param_ij%llm == dftb_param_ji%llm)
259 llm = dftb_param_ij%llm
260 fmatij => dftb_param_ij%fmat
261 smatij => dftb_param_ij%smat
262 fmatji => dftb_param_ji%fmat
263 smatji => dftb_param_ji%smat
265 n_urpoly = dftb_param_ij%n_urpoly
266 urep_cut = dftb_param_ij%urep_cut
267 urep = dftb_param_ij%urep
268 spxr => dftb_param_ij%spxr
269 scoeff => dftb_param_ij%scoeff
270 spdim = dftb_param_ij%spdim
271 s_cut = dftb_param_ij%s_cut
272 srep = dftb_param_ij%srep
273 surr = dftb_param_ij%surr
275 dr = sqrt(sum(rij(:)**2))
276 IF (nint(dr/dgrd) <= ngrdcut)
THEN
281 ic = cell_to_index(cell(1), cell(2), cell(3))
285 icol = max(iatom, jatom)
286 irow = min(iatom, jatom)
287 NULLIFY (sblock, fblock)
289 row=irow, col=icol, block=sblock, found=found)
292 row=irow, col=icol, block=fblock, found=found)
295 IF (calculate_forces)
THEN
298 row=irow, col=icol, block=pblock, found=found)
299 cpassert(
ASSOCIATED(pblock))
302 row=irow, col=icol, block=wblock, found=found)
303 cpassert(
ASSOCIATED(wblock))
304 IF (dftb_control%self_consistent)
THEN
306 NULLIFY (dsblocks(i)%block)
308 row=irow, col=icol, block=dsblocks(i)%block, found=found)
314 IF (iatom == jatom .AND. dr < 0.001_dp)
THEN
317 sblock(i, i) = sblock(i, i) + 1._dp
318 fblock(i, i) = fblock(i, i) + skself(
orbptr(i))
323 llm, lmaxi, lmaxj, irow, iatom)
325 llm, lmaxi, lmaxj, irow, iatom)
326 IF (calculate_forces)
THEN
334 IF (irow == iatom) f0 = -1.0_dp
336 ALLOCATE (dfblock(n1, n2), dsblock(n1, n2))
340 dfblock = 0._dp; dsblock = 0._dp
342 drij(i) = rij(i) - ddr*f0
344 llm, lmaxi, lmaxj, irow, iatom)
346 llm, lmaxi, lmaxj, irow, iatom)
351 drij(i) = rij(i) + ddr*f0
353 llm, lmaxi, lmaxj, irow, iatom)
355 llm, lmaxi, lmaxj, irow, iatom)
357 dfblock = dfblock/(2.0_dp*ddr)
358 dsblock = dsblock/(2.0_dp*ddr)
360 foab = 2.0_dp*sum(dfblock*pblock)
361 fow = -2.0_dp*sum(dsblock*wblock)
363 force_ab(i) = force_ab(i) + foab
364 force_w(i) = force_w(i) + fow
365 IF (dftb_control%self_consistent)
THEN
366 cpassert(
ASSOCIATED(dsblocks(i + 1)%block))
367 dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock
371 IF (iatom == jatom) f0 = 0.5_dp*f0
375 DEALLOCATE (dfblock, dsblock)
379 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp))
THEN
380 atom_a = atom_of_kind(iatom)
381 atom_b = atom_of_kind(jatom)
382 IF (irow == iatom) force_ab = -force_ab
383 IF (irow == iatom) force_w = -force_w
384 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
385 force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
386 force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - force_w(:)
387 force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + force_w(:)
393 IF ((dr <= urep_cut .OR. spdim > 0) .AND. dr > 0.001_dp)
THEN
395 CALL urep_egr(rij, dr, erepij, force_rr, &
396 n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, calculate_forces)
398 IF (atprop%energy)
THEN
399 atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*erepij
400 atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*erepij
402 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp))
THEN
403 atom_a = atom_of_kind(iatom)
404 atom_b = atom_of_kind(jatom)
405 force(ikind)%repulsive(:, atom_a) = &
406 force(ikind)%repulsive(:, atom_a) - force_rr(:)
407 force(jkind)%repulsive(:, atom_b) = &
408 force(jkind)%repulsive(:, atom_b) + force_rr(:)
411 IF (iatom == jatom) f0 = -0.5_dp
420 DO i = 1,
SIZE(matrix_s, 1)
425 DO i = 1,
SIZE(matrix_h, 1)
432 CALL para_env%sum(erep)
433 energy%repulsive = erep
435 CALL section_vals_val_get(qs_env%input,
"DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
437 qs_env%input,
"DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"),
cp_p_file))
THEN
441 after = min(max(after, 1), 16)
444 output_unit=iw, omit_headers=omit_headers)
448 "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
452 qs_env%input,
"DFT%PRINT%AO_MATRICES/OVERLAP"),
cp_p_file))
THEN
456 after = min(max(after, 1), 16)
459 output_unit=iw, omit_headers=omit_headers)
462 qs_env%input,
"DFT%PRINT%AO_MATRICES/DERIVATIVES"),
cp_p_file))
THEN
463 DO i = 2,
SIZE(matrix_s, 1)
465 output_unit=iw, omit_headers=omit_headers)
471 "DFT%PRINT%AO_MATRICES/OVERLAP")
474 IF (calculate_forces)
THEN
475 IF (
SIZE(matrix_p, 1) == 2)
THEN
477 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
479 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, alpha_scalar=1.0_dp, &
485 CALL timestop(handle)
497 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
499 CHARACTER(len=*),
PARAMETER :: routinen =
'build_dftb_ks_matrix'
501 INTEGER :: atom_a, handle, iatom, ikind, img, &
502 ispin, natom, nkind, nspins, &
504 REAL(kind=
dp) :: pc_ener, qmmm_el, zeff
505 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: mix_charge
506 REAL(kind=
dp),
DIMENSION(:),
POINTER :: mcharge, occupation_numbers
507 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: charges
510 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_p1, mo_derivs
511 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix, matrix_h, matrix_p, matrix_s
518 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
524 CALL timeset(routinen, handle)
525 NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
526 ks_matrix, rho, energy, scf_env)
528 cpassert(
ASSOCIATED(qs_env))
531 dft_control=dft_control, &
532 atomic_kind_set=atomic_kind_set, &
533 qs_kind_set=qs_kind_set, &
534 matrix_h_kp=matrix_h, &
537 matrix_ks_kp=ks_matrix, &
541 energy%hartree = 0.0_dp
542 energy%qmmm_el = 0.0_dp
545 nspins = dft_control%nspins
546 cpassert(
ASSOCIATED(matrix_h))
547 cpassert(
ASSOCIATED(rho))
548 cpassert(
SIZE(ks_matrix) > 0)
551 DO img = 1,
SIZE(ks_matrix, 2)
553 CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
557 IF (dft_control%qs_control%dftb_control%self_consistent)
THEN
559 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
560 matrix_s_kp=matrix_s)
562 natom =
SIZE(particle_set)
563 ALLOCATE (charges(natom, nspins))
567 ALLOCATE (mcharge(natom))
568 nkind =
SIZE(atomic_kind_set)
571 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
574 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
575 mcharge(atom_a) = zeff - sum(charges(atom_a, 1:nspins))
580 IF ((.NOT. dft_control%qs_control%do_ls_scf) .AND. &
582 CALL get_qs_env(qs_env=qs_env, scf_env=scf_env)
583 ALLOCATE (mix_charge(
SIZE(mcharge), 1))
584 mix_charge(:, 1) = mcharge(:)
585 CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, &
586 mix_charge, para_env, scf_env%iter_count, &
587 scc_mixer=dft_control%qs_control%dftb_control%tblite_scc_mixer, &
588 tblite_mixer_iterations= &
589 dft_control%qs_control%dftb_control%tblite_mixer_iterations, &
590 tblite_mixer_damping=dft_control%qs_control%dftb_control%tblite_mixer_damping, &
591 tblite_mixer_memory=dft_control%qs_control%dftb_control%tblite_mixer_memory, &
592 tblite_mixer_omega0=dft_control%qs_control%dftb_control%tblite_mixer_omega0, &
593 tblite_mixer_min_weight= &
594 dft_control%qs_control%dftb_control%tblite_mixer_min_weight, &
595 tblite_mixer_max_weight= &
596 dft_control%qs_control%dftb_control%tblite_mixer_max_weight, &
597 tblite_mixer_weight_factor= &
598 dft_control%qs_control%dftb_control%tblite_mixer_weight_factor)
599 mcharge(:) = mix_charge(:, 1)
600 DEALLOCATE (mix_charge)
604 calculate_forces, just_energy)
607 calculate_forces, just_energy)
613 IF (qs_env%qmmm)
THEN
614 cpassert(
SIZE(ks_matrix, 2) == 1)
617 CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
621 CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
622 matrix_p1(ispin)%matrix, qmmm_el)
623 energy%qmmm_el = energy%qmmm_el + qmmm_el
625 pc_ener = qs_env%ks_qmmm_env%pc_ener
626 energy%qmmm_el = energy%qmmm_el + pc_ener
629 energy%total = energy%core + energy%hartree + energy%qmmm_el + energy%efield + &
630 energy%repulsive + energy%dispersion + energy%dftb3
632 IF (dft_control%qs_control%dftb_control%self_consistent)
THEN
635 IF (output_unit > 0)
THEN
636 WRITE (unit=output_unit, fmt=
"(/,(T9,A,T60,F20.10))") &
637 "Repulsive pair potential energy: ", energy%repulsive, &
638 "Zeroth order Hamiltonian energy: ", energy%core, &
639 "Charge fluctuation energy: ", energy%hartree, &
640 "London dispersion energy: ", energy%dispersion
641 IF (abs(energy%efield) > 1.e-12_dp)
THEN
642 WRITE (unit=output_unit, fmt=
"(T9,A,T60,F20.10)") &
643 "Electric field interaction energy: ", energy%efield
645 IF (dft_control%qs_control%dftb_control%dftb3_diagonal)
THEN
646 WRITE (unit=output_unit, fmt=
"(T9,A,T60,F20.10)") &
647 "DFTB3 3rd Order Energy Correction ", energy%dftb3
649 IF (qs_env%qmmm)
THEN
650 WRITE (unit=output_unit, fmt=
"(T9,A,T60,F20.10)") &
651 "QM/MM Electrostatic energy: ", energy%qmmm_el
655 "PRINT%DETAILED_ENERGY")
658 IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy)
THEN
659 cpassert(
SIZE(ks_matrix, 2) == 1)
661 TYPE(
mo_set_type),
DIMENSION(:),
POINTER :: mo_array
662 CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
663 DO ispin = 1,
SIZE(mo_derivs)
665 mo_coeff_b=mo_coeff, occupation_numbers=occupation_numbers)
666 IF (.NOT. mo_array(ispin)%use_mo_coeff_b)
THEN
669 CALL dbcsr_multiply(
'n',
'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
670 0.0_dp, mo_derivs(ispin)%matrix)
675 CALL timestop(handle)
687 TYPE(qs_environment_type),
POINTER :: qs_env
688 INTEGER,
INTENT(IN) :: nderivative
689 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
691 CHARACTER(LEN=*),
PARAMETER :: routinen =
'build_dftb_overlap'
693 INTEGER :: handle, i, iatom, icol, ikind, indder, irow, j, jatom, jkind, l1, l2, la, lb, &
694 llm, lmaxi, lmaxj, m, n1, n2, natom, natorb_a, natorb_b, ngrd, ngrdcut, nkind
695 LOGICAL :: defined, found
696 REAL(kind=dp) :: ddr, dgrd, dr, f0
697 REAL(kind=dp),
DIMENSION(0:3) :: skself
698 REAL(kind=dp),
DIMENSION(3) :: drij, rij
699 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: dsblock, dsblockm, sblock, smatij, smatji
700 REAL(kind=dp),
DIMENSION(:, :, :),
POINTER :: dsblock1
701 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
702 TYPE(block_p_type),
DIMENSION(2:10) :: dsblocks
703 TYPE(cp_logger_type),
POINTER :: logger
704 TYPE(dft_control_type),
POINTER :: dft_control
705 TYPE(dftb_control_type),
POINTER :: dftb_control
706 TYPE(neighbor_list_iterator_p_type), &
707 DIMENSION(:),
POINTER :: nl_iterator
708 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
710 TYPE(qs_dftb_atom_type),
POINTER :: dftb_kind_a, dftb_kind_b
711 TYPE(qs_dftb_pairpot_type),
DIMENSION(:, :), &
712 POINTER :: dftb_potential
713 TYPE(qs_dftb_pairpot_type),
POINTER :: dftb_param_ij, dftb_param_ji
714 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
716 CALL timeset(routinen, handle)
723 DO l1 = 0, max(la, lb)
724 DO l2 = 0, min(l1, la, lb)
727 iptr(l1, l2, m, la, lb) = llm
735 logger => cp_get_default_logger()
737 NULLIFY (atomic_kind_set, qs_kind_set, sab_orb)
739 CALL get_qs_env(qs_env=qs_env, &
740 atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
741 dft_control=dft_control)
743 dftb_control => dft_control%qs_control%dftb_control
745 NULLIFY (dftb_potential)
746 CALL get_qs_env(qs_env=qs_env, &
747 dftb_potential=dftb_potential)
749 nkind =
SIZE(atomic_kind_set)
752 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
753 CALL setup_matrices1(qs_env, nderivative, matrix_s,
'OVERLAP', sab_orb)
755 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
756 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
757 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
758 iatom=iatom, jatom=jatom, r=rij)
760 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
761 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
762 CALL get_dftb_atom_param(dftb_kind_a, &
763 defined=defined, lmax=lmaxi, skself=skself, &
766 IF (.NOT. defined .OR. natorb_a < 1) cycle
768 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
769 CALL get_dftb_atom_param(dftb_kind_b, &
770 defined=defined, lmax=lmaxj, natorb=natorb_b)
772 IF (.NOT. defined .OR. natorb_b < 1) cycle
775 dftb_param_ij => dftb_potential(ikind, jkind)
776 dftb_param_ji => dftb_potential(jkind, ikind)
778 ngrd = dftb_param_ij%ngrd
779 ngrdcut = dftb_param_ij%ngrdcut
780 dgrd = dftb_param_ij%dgrd
781 ddr = dgrd*dftb_fd_deriv_step
782 cpassert(dftb_param_ij%llm == dftb_param_ji%llm)
783 llm = dftb_param_ij%llm
784 smatij => dftb_param_ij%smat
785 smatji => dftb_param_ji%smat
787 dr = sqrt(sum(rij(:)**2))
788 IF (nint(dr/dgrd) <= ngrdcut)
THEN
790 icol = max(iatom, jatom); irow = min(iatom, jatom)
793 CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
794 row=irow, col=icol, block=sblock, found=found)
797 IF (nderivative > 0)
THEN
798 DO i = 2,
SIZE(matrix_s, 1)
799 NULLIFY (dsblocks(i)%block)
800 CALL dbcsr_get_block_p(matrix=matrix_s(i)%matrix, &
801 row=irow, col=icol, block=dsblocks(i)%block, found=found)
805 IF (iatom == jatom .AND. dr < 0.001_dp)
THEN
808 sblock(i, i) = sblock(i, i) + 1._dp
812 CALL compute_block_sk(sblock, smatij, smatji, rij, ngrd, ngrdcut, dgrd, &
813 llm, lmaxi, lmaxj, irow, iatom)
815 IF (nderivative >= 1)
THEN
816 n1 =
SIZE(sblock, 1); n2 =
SIZE(sblock, 2)
819 ALLOCATE (dsblock1(n1, n2, 3), dsblock(n1, n2), dsblockm(n1, n2))
822 dsblock = 0._dp; dsblockm = 0.0_dp
824 f0 = 1.0_dp;
IF (irow == iatom) f0 = -1.0_dp
826 drij(i) = rij(i) - ddr*f0
827 CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
828 llm, lmaxi, lmaxj, irow, iatom)
830 drij(i) = rij(i) + ddr*f0
831 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
832 llm, lmaxi, lmaxj, irow, iatom)
834 dsblock1(:, :, i) = (dsblock + dsblockm)
835 dsblock = dsblock - dsblockm
836 dsblock = dsblock/(2.0_dp*ddr)
838 cpassert(
ASSOCIATED(dsblocks(i + 1)%block))
839 dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock
840 IF (nderivative > 1)
THEN
841 indder = indder + 5 - i
842 dsblocks(indder)%block = 0.0_dp
843 dsblocks(indder)%block = dsblocks(indder)%block + &
844 (dsblock1(:, :, i) - 2.0_dp*sblock)/ddr**2
848 IF (nderivative > 1)
THEN
851 dsblock = 0._dp; dsblockm = 0.0_dp
853 f0 = 1.0_dp;
IF (irow == iatom) f0 = -1.0_dp
855 drij(i) = rij(i) - ddr*f0; drij(j) = rij(j) - ddr*f0
856 CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
857 llm, lmaxi, lmaxj, irow, iatom)
859 drij(i) = rij(i) + ddr*f0; drij(j) = rij(j) + ddr*f0
860 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
861 llm, lmaxi, lmaxj, irow, iatom)
864 dsblocks(indder)%block = 0.0_dp
865 dsblocks(indder)%block = &
866 dsblocks(indder)%block + ( &
867 dsblock + dsblockm - dsblock1(:, :, i) - dsblock1(:, :, j) + 2.0_dp*sblock)/(2.0_dp*ddr**2)
872 DEALLOCATE (dsblock1, dsblock, dsblockm)
877 CALL neighbor_list_iterator_release(nl_iterator)
879 DO i = 1,
SIZE(matrix_s, 1)
880 CALL dbcsr_finalize(matrix_s(i)%matrix)
883 CALL timestop(handle)
895 SUBROUTINE setup_matrices1(qs_env, nderivative, matrices, mnames, sab_nl)
897 TYPE(qs_environment_type),
POINTER :: qs_env
898 INTEGER,
INTENT(IN) :: nderivative
899 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrices
900 CHARACTER(LEN=*) :: mnames
901 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
904 CHARACTER(1) :: symmetry_type
905 CHARACTER(LEN=default_string_length) :: matnames
906 INTEGER :: i, natom, neighbor_list_id, nkind, nmat, &
908 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf
909 INTEGER,
DIMENSION(:),
POINTER :: row_blk_sizes
910 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
911 TYPE(dbcsr_distribution_type),
POINTER :: dbcsr_dist
912 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
913 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
915 NULLIFY (particle_set, atomic_kind_set)
917 CALL get_qs_env(qs_env=qs_env, &
918 atomic_kind_set=atomic_kind_set, &
919 qs_kind_set=qs_kind_set, &
920 particle_set=particle_set, &
921 dbcsr_dist=dbcsr_dist, &
922 neighbor_list_id=neighbor_list_id)
924 nkind =
SIZE(atomic_kind_set)
925 natom =
SIZE(particle_set)
927 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
929 ALLOCATE (first_sgf(natom))
930 ALLOCATE (last_sgf(natom))
932 CALL get_particle_set(particle_set, qs_kind_set, &
933 first_sgf=first_sgf, &
937 IF (nderivative == 0) nmat = 1
938 IF (nderivative == 1) nmat = 4
939 IF (nderivative == 2) nmat = 10
942 ALLOCATE (row_blk_sizes(natom))
943 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
945 CALL dbcsr_allocate_matrix_set(matrices, nmat)
950 matnames = trim(mnames)//
" DERIVATIVE MATRIX DFTB"
951 symmetry_type = dbcsr_type_antisymmetric
952 IF (i > 4) symmetry_type = dbcsr_type_symmetric
954 symmetry_type = dbcsr_type_symmetric
955 matnames = trim(mnames)//
" MATRIX DFTB"
957 ALLOCATE (matrices(i)%matrix)
958 CALL dbcsr_create(matrix=matrices(i)%matrix, &
959 name=trim(matnames), &
960 dist=dbcsr_dist, matrix_type=symmetry_type, &
961 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
963 CALL cp_dbcsr_alloc_block_from_nbl(matrices(i)%matrix, sab_nl)
966 DEALLOCATE (first_sgf)
967 DEALLOCATE (last_sgf)
969 DEALLOCATE (row_blk_sizes)
971 END SUBROUTINE setup_matrices1
982 SUBROUTINE setup_matrices2(qs_env, nderivative, nimg, matrices, mnames, sab_nl)
984 TYPE(qs_environment_type),
POINTER :: qs_env
985 INTEGER,
INTENT(IN) :: nderivative, nimg
986 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrices
987 CHARACTER(LEN=*) :: mnames
988 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
991 CHARACTER(1) :: symmetry_type
992 CHARACTER(LEN=default_string_length) :: matnames
993 INTEGER :: i, img, natom, neighbor_list_id, nkind, &
995 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf
996 INTEGER,
DIMENSION(:),
POINTER :: row_blk_sizes
997 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
998 TYPE(dbcsr_distribution_type),
POINTER :: dbcsr_dist
999 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
1000 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1002 NULLIFY (particle_set, atomic_kind_set)
1004 CALL get_qs_env(qs_env=qs_env, &
1005 atomic_kind_set=atomic_kind_set, &
1006 qs_kind_set=qs_kind_set, &
1007 particle_set=particle_set, &
1008 dbcsr_dist=dbcsr_dist, &
1009 neighbor_list_id=neighbor_list_id)
1011 nkind =
SIZE(atomic_kind_set)
1012 natom =
SIZE(particle_set)
1014 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
1016 ALLOCATE (first_sgf(natom))
1017 ALLOCATE (last_sgf(natom))
1019 CALL get_particle_set(particle_set, qs_kind_set, &
1020 first_sgf=first_sgf, &
1024 IF (nderivative == 0) nmat = 1
1025 IF (nderivative == 1) nmat = 4
1026 IF (nderivative == 2) nmat = 10
1029 ALLOCATE (row_blk_sizes(natom))
1030 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
1032 CALL dbcsr_allocate_matrix_set(matrices, nmat, nimg)
1038 matnames = trim(mnames)//
" DERIVATIVE MATRIX DFTB"
1039 symmetry_type = dbcsr_type_antisymmetric
1040 IF (i > 4) symmetry_type = dbcsr_type_symmetric
1042 symmetry_type = dbcsr_type_symmetric
1043 matnames = trim(mnames)//
" MATRIX DFTB"
1045 ALLOCATE (matrices(i, img)%matrix)
1046 CALL dbcsr_create(matrix=matrices(i, img)%matrix, &
1047 name=trim(matnames), &
1048 dist=dbcsr_dist, matrix_type=symmetry_type, &
1049 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
1050 mutable_work=.true.)
1051 CALL cp_dbcsr_alloc_block_from_nbl(matrices(i, img)%matrix, sab_nl)
1055 DEALLOCATE (first_sgf)
1056 DEALLOCATE (last_sgf)
1058 DEALLOCATE (row_blk_sizes)
1060 END SUBROUTINE setup_matrices2
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.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Holds information on atomic properties.
subroutine, public atprop_array_init(atarray, natom)
...
collect pointers to a block of reals
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_finalize(matrix)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Routines that link DBCSR and CP2K concepts together.
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Calculation of electric field contributions in TB.
subroutine, public efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
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, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Retrieve information from a kpoint environment.
Interface to the message passing library MPI.
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
subroutine, public charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count, scc_mixer, tblite_mixer_iterations, tblite_mixer_damping, tblite_mixer_memory, tblite_mixer_omega0, tblite_mixer_min_weight, tblite_mixer_max_weight, tblite_mixer_weight_factor)
Driver for TB SCC variable mixing, calls the requested method.
Calculation of Coulomb contributions in DFTB.
subroutine, public build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
Calculation of Overlap and Hamiltonian matrices in DFTB.
subroutine, public build_dftb_ks_matrix(qs_env, calculate_forces, just_energy)
...
integer, dimension(16), parameter orbptr
subroutine, public build_dftb_matrices(qs_env, para_env, calculate_forces)
...
subroutine, public build_dftb_overlap(qs_env, nderivative, matrix_s)
...
Definition of the DFTB parameter types.
Working with the DFTB parameter types.
subroutine, public urep_egr(rv, r, erep, derep, n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, dograd)
...
subroutine, public compute_block_sk(block, smatij, smatji, rij, ngrd, ngrdcut, dgrd, llm, lmaxi, lmaxj, irow, iatom)
...
integer, dimension(0:3, 0:3, 0:3, 0:3, 0:3), public iptr
subroutine, public get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
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.
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, 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.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition and initialisation of the mo data type.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
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...
module that contains the definitions of the scf types
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Provides all information about an atomic kind.
type for the atomic properties
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.