30 USE dbcsr_api,
ONLY: &
31 dbcsr_add, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, &
32 dbcsr_distribution_type, dbcsr_dot, dbcsr_finalize, dbcsr_get_block_p, dbcsr_multiply, &
33 dbcsr_p_type, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_symmetric
75#include "./base/base_uses.f90"
79 INTEGER,
DIMENSION(16),
PARAMETER ::
orbptr = (/0, 1, 1, 1, &
80 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3/)
84 CHARACTER(len=*),
PARAMETER,
PRIVATE :: modulen =
'qs_dftb_matrices'
100 LOGICAL,
INTENT(IN) :: calculate_forces
102 CHARACTER(LEN=*),
PARAMETER :: routinen =
'build_dftb_matrices'
104 INTEGER :: after, atom_a, atom_b, handle, i, iatom, ic, icol, ikind, img, irow, iw, jatom, &
105 jkind, l1, l2, la, lb, llm, lmaxi, lmaxj, m, n1, n2, n_urpoly, natorb_a, natorb_b, &
106 nderivatives, ngrd, ngrdcut, nimg, nkind, spdim
107 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind
108 INTEGER,
DIMENSION(3) :: cell
109 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
110 LOGICAL :: defined, found, omit_headers, use_virial
111 REAL(kind=
dp) :: ddr, dgrd, dr, erep, erepij, f0, f1, &
112 foab, fow, s_cut, urep_cut
113 REAL(kind=
dp),
DIMENSION(0:3) :: eta_a, eta_b, skself
114 REAL(kind=
dp),
DIMENSION(10) :: urep
115 REAL(kind=
dp),
DIMENSION(2) :: surr
116 REAL(kind=
dp),
DIMENSION(3) :: drij, force_ab, force_rr, force_w, rij, &
118 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: dfblock, dsblock, fblock, fmatij, &
119 fmatji, pblock, sblock, scoeff, &
120 smatij, smatji, spxr, wblock
125 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
130 DIMENSION(:),
POINTER :: nl_iterator
136 POINTER :: dftb_potential
140 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
145 CALL timeset(routinen, handle)
152 DO l1 = 0, max(la, lb)
153 DO l2 = 0, min(l1, la, lb)
156 iptr(l1, l2, m, la, lb) = llm
163 NULLIFY (logger, virial, atprop)
166 NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
167 qs_kind_set, sab_orb, ks_env)
171 atomic_kind_set=atomic_kind_set, &
172 qs_kind_set=qs_kind_set, &
173 matrix_h_kp=matrix_h, &
174 matrix_s_kp=matrix_s, &
176 dft_control=dft_control, &
179 dftb_control => dft_control%qs_control%dftb_control
180 nimg = dft_control%nimages
182 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
184 IF (dftb_control%self_consistent .AND. calculate_forces) nderivatives = 1
185 CALL setup_matrices2(qs_env, nderivatives, nimg, matrix_s,
"OVERLAP", sab_orb)
186 CALL setup_matrices2(qs_env, 0, nimg, matrix_h,
"CORE HAMILTONIAN", sab_orb)
190 NULLIFY (dftb_potential)
191 CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential)
192 NULLIFY (particle_set)
193 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
195 IF (calculate_forces)
THEN
196 NULLIFY (rho, force, matrix_w)
199 matrix_w_kp=matrix_w, &
204 IF (
SIZE(matrix_p, 1) == 2)
THEN
206 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
207 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
208 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
209 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
213 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
216 IF (atprop%energy)
THEN
220 NULLIFY (cell_to_index)
222 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
228 nkind =
SIZE(atomic_kind_set)
233 iatom=iatom, jatom=jatom, r=rij, cell=cell)
234 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
236 defined=defined, lmax=lmaxi, skself=skself, &
237 eta=eta_a, natorb=natorb_a)
238 IF (.NOT. defined .OR. natorb_a < 1) cycle
239 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
241 defined=defined, lmax=lmaxj, eta=eta_b, natorb=natorb_b)
243 IF (.NOT. defined .OR. natorb_b < 1) cycle
246 dftb_param_ij => dftb_potential(ikind, jkind)
247 dftb_param_ji => dftb_potential(jkind, ikind)
249 ngrd = dftb_param_ij%ngrd
250 ngrdcut = dftb_param_ij%ngrdcut
251 dgrd = dftb_param_ij%dgrd
253 cpassert(dftb_param_ij%llm == dftb_param_ji%llm)
254 llm = dftb_param_ij%llm
255 fmatij => dftb_param_ij%fmat
256 smatij => dftb_param_ij%smat
257 fmatji => dftb_param_ji%fmat
258 smatji => dftb_param_ji%smat
260 n_urpoly = dftb_param_ij%n_urpoly
261 urep_cut = dftb_param_ij%urep_cut
262 urep = dftb_param_ij%urep
263 spxr => dftb_param_ij%spxr
264 scoeff => dftb_param_ij%scoeff
265 spdim = dftb_param_ij%spdim
266 s_cut = dftb_param_ij%s_cut
267 srep = dftb_param_ij%srep
268 surr = dftb_param_ij%surr
270 dr = sqrt(sum(rij(:)**2))
271 IF (nint(dr/dgrd) <= ngrdcut)
THEN
276 ic = cell_to_index(cell(1), cell(2), cell(3))
280 icol = max(iatom, jatom)
281 irow = min(iatom, jatom)
282 NULLIFY (sblock, fblock)
283 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
284 row=irow, col=icol, block=sblock, found=found)
286 CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
287 row=irow, col=icol, block=fblock, found=found)
290 IF (calculate_forces)
THEN
292 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
293 row=irow, col=icol, block=pblock, found=found)
294 cpassert(
ASSOCIATED(pblock))
296 CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
297 row=irow, col=icol, block=wblock, found=found)
298 cpassert(
ASSOCIATED(wblock))
299 IF (dftb_control%self_consistent)
THEN
301 NULLIFY (dsblocks(i)%block)
302 CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
303 row=irow, col=icol, block=dsblocks(i)%block, found=found)
309 IF (iatom == jatom .AND. dr < 0.001_dp)
THEN
312 sblock(i, i) = sblock(i, i) + 1._dp
313 fblock(i, i) = fblock(i, i) + skself(
orbptr(i))
318 llm, lmaxi, lmaxj, irow, iatom)
320 llm, lmaxi, lmaxj, irow, iatom)
321 IF (calculate_forces)
THEN
329 IF (irow == iatom) f0 = -1.0_dp
331 ALLOCATE (dfblock(n1, n2), dsblock(n1, n2))
335 dfblock = 0._dp; dsblock = 0._dp
337 drij(i) = rij(i) - ddr*f0
339 llm, lmaxi, lmaxj, irow, iatom)
341 llm, lmaxi, lmaxj, irow, iatom)
346 drij(i) = rij(i) + ddr*f0
348 llm, lmaxi, lmaxj, irow, iatom)
350 llm, lmaxi, lmaxj, irow, iatom)
352 dfblock = dfblock/(2.0_dp*ddr)
353 dsblock = dsblock/(2.0_dp*ddr)
355 foab = 2.0_dp*sum(dfblock*pblock)
356 fow = -2.0_dp*sum(dsblock*wblock)
358 force_ab(i) = force_ab(i) + foab
359 force_w(i) = force_w(i) + fow
360 IF (dftb_control%self_consistent)
THEN
361 cpassert(
ASSOCIATED(dsblocks(i + 1)%block))
362 dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock
366 IF (iatom == jatom) f0 = 0.5_dp*f0
369 IF (atprop%stress)
THEN
377 DEALLOCATE (dfblock, dsblock)
381 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp))
THEN
382 atom_a = atom_of_kind(iatom)
383 atom_b = atom_of_kind(jatom)
384 IF (irow == iatom) force_ab = -force_ab
385 IF (irow == iatom) force_w = -force_w
386 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
387 force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
388 force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - force_w(:)
389 force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + force_w(:)
395 IF ((dr <= urep_cut .OR. spdim > 0) .AND. dr > 0.001_dp)
THEN
397 CALL urep_egr(rij, dr, erepij, force_rr, &
398 n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, calculate_forces)
400 IF (atprop%energy)
THEN
401 atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*erepij
402 atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*erepij
404 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp))
THEN
405 atom_a = atom_of_kind(iatom)
406 atom_b = atom_of_kind(jatom)
407 force(ikind)%repulsive(:, atom_a) = &
408 force(ikind)%repulsive(:, atom_a) - force_rr(:)
409 force(jkind)%repulsive(:, atom_b) = &
410 force(jkind)%repulsive(:, atom_b) + force_rr(:)
413 IF (iatom == jatom) f0 = -0.5_dp
415 IF (atprop%stress)
THEN
426 DO i = 1,
SIZE(matrix_s, 1)
428 CALL dbcsr_finalize(matrix_s(i, img)%matrix)
431 DO i = 1,
SIZE(matrix_h, 1)
433 CALL dbcsr_finalize(matrix_h(i, img)%matrix)
438 CALL para_env%sum(erep)
439 energy%repulsive = erep
441 CALL section_vals_val_get(qs_env%input,
"DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
443 qs_env%input,
"DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"),
cp_p_file))
THEN
447 after = min(max(after, 1), 16)
450 output_unit=iw, omit_headers=omit_headers)
454 "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
458 qs_env%input,
"DFT%PRINT%AO_MATRICES/OVERLAP"),
cp_p_file))
THEN
462 after = min(max(after, 1), 16)
465 output_unit=iw, omit_headers=omit_headers)
468 qs_env%input,
"DFT%PRINT%AO_MATRICES/DERIVATIVES"),
cp_p_file))
THEN
469 DO i = 2,
SIZE(matrix_s, 1)
471 output_unit=iw, omit_headers=omit_headers)
477 "DFT%PRINT%AO_MATRICES/OVERLAP")
480 IF (calculate_forces)
THEN
481 IF (
SIZE(matrix_p, 1) == 2)
THEN
483 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
489 CALL timestop(handle)
501 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
503 CHARACTER(len=*),
PARAMETER :: routinen =
'build_dftb_ks_matrix'
505 INTEGER :: atom_a, handle, iatom, ikind, img, &
506 ispin, natom, nkind, nspins, &
508 REAL(kind=
dp) :: pc_ener, qmmm_el, zeff
509 REAL(kind=
dp),
DIMENSION(:),
POINTER :: mcharge, occupation_numbers
510 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: charges
513 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_p1, mo_derivs
514 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix, matrix_h, matrix_p, matrix_s
515 TYPE(dbcsr_type),
POINTER :: mo_coeff
521 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
526 CALL timeset(routinen, handle)
527 NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, &
528 ks_matrix, rho, energy)
530 cpassert(
ASSOCIATED(qs_env))
533 dft_control=dft_control, &
534 atomic_kind_set=atomic_kind_set, &
535 qs_kind_set=qs_kind_set, &
536 matrix_h_kp=matrix_h, &
539 matrix_ks_kp=ks_matrix, &
543 energy%hartree = 0.0_dp
544 energy%qmmm_el = 0.0_dp
547 nspins = dft_control%nspins
548 cpassert(
ASSOCIATED(matrix_h))
549 cpassert(
ASSOCIATED(rho))
550 cpassert(
SIZE(ks_matrix) > 0)
553 DO img = 1,
SIZE(ks_matrix, 2)
555 CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
559 IF (dft_control%qs_control%dftb_control%self_consistent)
THEN
561 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
562 matrix_s_kp=matrix_s)
564 natom =
SIZE(particle_set)
565 ALLOCATE (charges(natom, nspins))
569 ALLOCATE (mcharge(natom))
570 nkind =
SIZE(atomic_kind_set)
573 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
576 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
577 mcharge(atom_a) = zeff - sum(charges(atom_a, 1:nspins))
583 calculate_forces, just_energy)
586 calculate_forces, just_energy)
592 IF (qs_env%qmmm)
THEN
593 cpassert(
SIZE(ks_matrix, 2) == 1)
596 CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
600 CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
601 matrix_p1(ispin)%matrix, qmmm_el)
602 energy%qmmm_el = energy%qmmm_el + qmmm_el
604 pc_ener = qs_env%ks_qmmm_env%pc_ener
605 energy%qmmm_el = energy%qmmm_el + pc_ener
608 energy%total = energy%core + energy%hartree + energy%qmmm_el + energy%efield + &
609 energy%repulsive + energy%dispersion + energy%dftb3
613 IF (output_unit > 0)
THEN
614 WRITE (unit=output_unit, fmt=
"(/,(T9,A,T60,F20.10))") &
615 "Repulsive pair potential energy: ", energy%repulsive, &
616 "Zeroth order Hamiltonian energy: ", energy%core, &
617 "Charge fluctuation energy: ", energy%hartree, &
618 "London dispersion energy: ", energy%dispersion
619 IF (abs(energy%efield) > 1.e-12_dp)
THEN
620 WRITE (unit=output_unit, fmt=
"(T9,A,T60,F20.10)") &
621 "Electric field interaction energy: ", energy%efield
623 IF (dft_control%qs_control%dftb_control%dftb3_diagonal)
THEN
624 WRITE (unit=output_unit, fmt=
"(T9,A,T60,F20.10)") &
625 "DFTB3 3rd Order Energy Correction ", energy%dftb3
627 IF (qs_env%qmmm)
THEN
628 WRITE (unit=output_unit, fmt=
"(T9,A,T60,F20.10)") &
629 "QM/MM Electrostatic energy: ", energy%qmmm_el
633 "PRINT%DETAILED_ENERGY")
635 IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy)
THEN
636 cpassert(
SIZE(ks_matrix, 2) == 1)
638 TYPE(
mo_set_type),
DIMENSION(:),
POINTER :: mo_array
639 CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
640 DO ispin = 1,
SIZE(mo_derivs)
642 mo_coeff_b=mo_coeff, occupation_numbers=occupation_numbers)
643 IF (.NOT. mo_array(ispin)%use_mo_coeff_b)
THEN
646 CALL dbcsr_multiply(
'n',
'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
647 0.0_dp, mo_derivs(ispin)%matrix)
652 CALL timestop(handle)
664 TYPE(qs_environment_type),
POINTER :: qs_env
665 INTEGER,
INTENT(IN) :: nderivative
666 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
668 CHARACTER(LEN=*),
PARAMETER :: routinen =
'build_dftb_overlap'
670 INTEGER :: handle, i, iatom, icol, ikind, indder, irow, j, jatom, jkind, l1, l2, la, lb, &
671 llm, lmaxi, lmaxj, m, n1, n2, natom, natorb_a, natorb_b, ngrd, ngrdcut, nkind
672 LOGICAL :: defined, found
673 REAL(kind=dp) :: ddr, dgrd, dr, f0
674 REAL(kind=dp),
DIMENSION(0:3) :: skself
675 REAL(kind=dp),
DIMENSION(3) :: drij, rij
676 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: dsblock, dsblockm, sblock, smatij, smatji
677 REAL(kind=dp),
DIMENSION(:, :, :),
POINTER :: dsblock1
678 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
679 TYPE(block_p_type),
DIMENSION(2:10) :: dsblocks
680 TYPE(cp_logger_type),
POINTER :: logger
681 TYPE(dft_control_type),
POINTER :: dft_control
682 TYPE(dftb_control_type),
POINTER :: dftb_control
683 TYPE(neighbor_list_iterator_p_type), &
684 DIMENSION(:),
POINTER :: nl_iterator
685 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
687 TYPE(qs_dftb_atom_type),
POINTER :: dftb_kind_a, dftb_kind_b
688 TYPE(qs_dftb_pairpot_type),
DIMENSION(:, :), &
689 POINTER :: dftb_potential
690 TYPE(qs_dftb_pairpot_type),
POINTER :: dftb_param_ij, dftb_param_ji
691 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
693 CALL timeset(routinen, handle)
700 DO l1 = 0, max(la, lb)
701 DO l2 = 0, min(l1, la, lb)
704 iptr(l1, l2, m, la, lb) = llm
712 logger => cp_get_default_logger()
714 NULLIFY (atomic_kind_set, qs_kind_set, sab_orb)
716 CALL get_qs_env(qs_env=qs_env, &
717 atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
718 dft_control=dft_control)
720 dftb_control => dft_control%qs_control%dftb_control
722 NULLIFY (dftb_potential)
723 CALL get_qs_env(qs_env=qs_env, &
724 dftb_potential=dftb_potential)
726 nkind =
SIZE(atomic_kind_set)
729 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
730 CALL setup_matrices1(qs_env, nderivative, matrix_s,
'OVERLAP', sab_orb)
732 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
733 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
734 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
735 iatom=iatom, jatom=jatom, r=rij)
737 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
738 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
739 CALL get_dftb_atom_param(dftb_kind_a, &
740 defined=defined, lmax=lmaxi, skself=skself, &
743 IF (.NOT. defined .OR. natorb_a < 1) cycle
745 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
746 CALL get_dftb_atom_param(dftb_kind_b, &
747 defined=defined, lmax=lmaxj, natorb=natorb_b)
749 IF (.NOT. defined .OR. natorb_b < 1) cycle
752 dftb_param_ij => dftb_potential(ikind, jkind)
753 dftb_param_ji => dftb_potential(jkind, ikind)
755 ngrd = dftb_param_ij%ngrd
756 ngrdcut = dftb_param_ij%ngrdcut
757 dgrd = dftb_param_ij%dgrd
759 cpassert(dftb_param_ij%llm == dftb_param_ji%llm)
760 llm = dftb_param_ij%llm
761 smatij => dftb_param_ij%smat
762 smatji => dftb_param_ji%smat
764 dr = sqrt(sum(rij(:)**2))
765 IF (nint(dr/dgrd) <= ngrdcut)
THEN
767 icol = max(iatom, jatom); irow = min(iatom, jatom)
770 CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
771 row=irow, col=icol, block=sblock, found=found)
774 IF (nderivative .GT. 0)
THEN
775 DO i = 2,
SIZE(matrix_s, 1)
776 NULLIFY (dsblocks(i)%block)
777 CALL dbcsr_get_block_p(matrix=matrix_s(i)%matrix, &
778 row=irow, col=icol, block=dsblocks(i)%block, found=found)
782 IF (iatom == jatom .AND. dr < 0.001_dp)
THEN
785 sblock(i, i) = sblock(i, i) + 1._dp
789 CALL compute_block_sk(sblock, smatij, smatji, rij, ngrd, ngrdcut, dgrd, &
790 llm, lmaxi, lmaxj, irow, iatom)
792 IF (nderivative .GE. 1)
THEN
793 n1 =
SIZE(sblock, 1); n2 =
SIZE(sblock, 2)
796 ALLOCATE (dsblock1(n1, n2, 3), dsblock(n1, n2), dsblockm(n1, n2))
799 dsblock = 0._dp; dsblockm = 0.0_dp
801 f0 = 1.0_dp;
IF (irow == iatom) f0 = -1.0_dp
803 drij(i) = rij(i) - ddr*f0
804 CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
805 llm, lmaxi, lmaxj, irow, iatom)
807 drij(i) = rij(i) + ddr*f0
808 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
809 llm, lmaxi, lmaxj, irow, iatom)
811 dsblock1(:, :, i) = (dsblock + dsblockm)
812 dsblock = dsblock - dsblockm
813 dsblock = dsblock/(2.0_dp*ddr)
815 cpassert(
ASSOCIATED(dsblocks(i + 1)%block))
816 dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock
817 IF (nderivative .GT. 1)
THEN
818 indder = indder + 5 - i
819 dsblocks(indder)%block = 0.0_dp
820 dsblocks(indder)%block = dsblocks(indder)%block + &
821 (dsblock1(:, :, i) - 2.0_dp*sblock)/ddr**2
825 IF (nderivative .GT. 1)
THEN
828 dsblock = 0._dp; dsblockm = 0.0_dp
830 f0 = 1.0_dp;
IF (irow == iatom) f0 = -1.0_dp
832 drij(i) = rij(i) - ddr*f0; drij(j) = rij(j) - ddr*f0
833 CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
834 llm, lmaxi, lmaxj, irow, iatom)
836 drij(i) = rij(i) + ddr*f0; drij(j) = rij(j) + ddr*f0
837 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
838 llm, lmaxi, lmaxj, irow, iatom)
841 dsblocks(indder)%block = 0.0_dp
842 dsblocks(indder)%block = &
843 dsblocks(indder)%block + ( &
844 dsblock + dsblockm - dsblock1(:, :, i) - dsblock1(:, :, j) + 2.0_dp*sblock)/(2.0_dp*ddr**2)
849 DEALLOCATE (dsblock1, dsblock, dsblockm)
854 CALL neighbor_list_iterator_release(nl_iterator)
856 DO i = 1,
SIZE(matrix_s, 1)
857 CALL dbcsr_finalize(matrix_s(i)%matrix)
860 CALL timestop(handle)
872 SUBROUTINE setup_matrices1(qs_env, nderivative, matrices, mnames, sab_nl)
874 TYPE(qs_environment_type),
POINTER :: qs_env
875 INTEGER,
INTENT(IN) :: nderivative
876 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrices
877 CHARACTER(LEN=*) :: mnames
878 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
881 CHARACTER(1) :: symmetry_type
882 CHARACTER(LEN=default_string_length) :: matnames
883 INTEGER :: i, natom, neighbor_list_id, nkind, nmat, &
885 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf
886 INTEGER,
DIMENSION(:),
POINTER :: row_blk_sizes
887 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
888 TYPE(dbcsr_distribution_type),
POINTER :: dbcsr_dist
889 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
890 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
892 NULLIFY (particle_set, atomic_kind_set)
894 CALL get_qs_env(qs_env=qs_env, &
895 atomic_kind_set=atomic_kind_set, &
896 qs_kind_set=qs_kind_set, &
897 particle_set=particle_set, &
898 dbcsr_dist=dbcsr_dist, &
899 neighbor_list_id=neighbor_list_id)
901 nkind =
SIZE(atomic_kind_set)
902 natom =
SIZE(particle_set)
904 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
906 ALLOCATE (first_sgf(natom))
907 ALLOCATE (last_sgf(natom))
909 CALL get_particle_set(particle_set, qs_kind_set, &
910 first_sgf=first_sgf, &
914 IF (nderivative == 0) nmat = 1
915 IF (nderivative == 1) nmat = 4
916 IF (nderivative == 2) nmat = 10
919 ALLOCATE (row_blk_sizes(natom))
920 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
922 CALL dbcsr_allocate_matrix_set(matrices, nmat)
927 matnames = trim(mnames)//
" DERIVATIVE MATRIX DFTB"
928 symmetry_type = dbcsr_type_antisymmetric
929 IF (i .GT. 4) symmetry_type = dbcsr_type_symmetric
931 symmetry_type = dbcsr_type_symmetric
932 matnames = trim(mnames)//
" MATRIX DFTB"
934 ALLOCATE (matrices(i)%matrix)
935 CALL dbcsr_create(matrix=matrices(i)%matrix, &
936 name=trim(matnames), &
937 dist=dbcsr_dist, matrix_type=symmetry_type, &
938 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
939 nze=0, mutable_work=.true.)
940 CALL cp_dbcsr_alloc_block_from_nbl(matrices(i)%matrix, sab_nl)
943 DEALLOCATE (first_sgf)
944 DEALLOCATE (last_sgf)
946 DEALLOCATE (row_blk_sizes)
948 END SUBROUTINE setup_matrices1
959 SUBROUTINE setup_matrices2(qs_env, nderivative, nimg, matrices, mnames, sab_nl)
961 TYPE(qs_environment_type),
POINTER :: qs_env
962 INTEGER,
INTENT(IN) :: nderivative, nimg
963 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrices
964 CHARACTER(LEN=*) :: mnames
965 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
968 CHARACTER(1) :: symmetry_type
969 CHARACTER(LEN=default_string_length) :: matnames
970 INTEGER :: i, img, natom, neighbor_list_id, nkind, &
972 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf
973 INTEGER,
DIMENSION(:),
POINTER :: row_blk_sizes
974 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
975 TYPE(dbcsr_distribution_type),
POINTER :: dbcsr_dist
976 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
977 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
979 NULLIFY (particle_set, atomic_kind_set)
981 CALL get_qs_env(qs_env=qs_env, &
982 atomic_kind_set=atomic_kind_set, &
983 qs_kind_set=qs_kind_set, &
984 particle_set=particle_set, &
985 dbcsr_dist=dbcsr_dist, &
986 neighbor_list_id=neighbor_list_id)
988 nkind =
SIZE(atomic_kind_set)
989 natom =
SIZE(particle_set)
991 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
993 ALLOCATE (first_sgf(natom))
994 ALLOCATE (last_sgf(natom))
996 CALL get_particle_set(particle_set, qs_kind_set, &
997 first_sgf=first_sgf, &
1001 IF (nderivative == 0) nmat = 1
1002 IF (nderivative == 1) nmat = 4
1003 IF (nderivative == 2) nmat = 10
1006 ALLOCATE (row_blk_sizes(natom))
1007 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
1009 CALL dbcsr_allocate_matrix_set(matrices, nmat, nimg)
1015 matnames = trim(mnames)//
" DERIVATIVE MATRIX DFTB"
1016 symmetry_type = dbcsr_type_antisymmetric
1017 IF (i .GT. 4) symmetry_type = dbcsr_type_symmetric
1019 symmetry_type = dbcsr_type_symmetric
1020 matnames = trim(mnames)//
" MATRIX DFTB"
1022 ALLOCATE (matrices(i, img)%matrix)
1023 CALL dbcsr_create(matrix=matrices(i, img)%matrix, &
1024 name=trim(matnames), &
1025 dist=dbcsr_dist, matrix_type=symmetry_type, &
1026 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
1027 nze=0, mutable_work=.true.)
1028 CALL cp_dbcsr_alloc_block_from_nbl(matrices(i, img)%matrix, sab_nl)
1032 DEALLOCATE (first_sgf)
1033 DEALLOCATE (last_sgf)
1035 DEALLOCATE (row_blk_sizes)
1037 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...
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)
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.
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, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
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)
Get attributes of an atomic kind set.
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, 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, 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_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, 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, 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, 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_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, 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...
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.