68#include "./base/base_uses.f90"
74 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_ks_atom'
107 SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, &
108 kind_set_external, oce_external, sab_external, kscale, &
109 kintegral, kforce, fscale)
112 TYPE(
dbcsr_p_type),
DIMENSION(*),
INTENT(INOUT) :: ksmat, pmat
113 LOGICAL,
INTENT(IN) :: forces
114 LOGICAL,
INTENT(IN),
OPTIONAL :: tddft
116 POINTER :: rho_atom_external
118 POINTER :: kind_set_external
121 OPTIONAL,
POINTER :: sab_external
122 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: kscale, kintegral, kforce
123 REAL(kind=
dp),
DIMENSION(2),
INTENT(IN),
OPTIONAL :: fscale
125 CHARACTER(len=*),
PARAMETER :: routinen =
'update_ks_atom'
127 INTEGER :: bo(2), handle, ia_kind, iac, iat, iatom, ibc, ikind, img, ip, ispin, ja_kind, &
128 jatom, jkind, ka_kind, kac, katom, kbc, kkind, ldcpc, max_gau, max_nsgf, n_cont_a, &
129 n_cont_b, nat, natom, nimages, nkind, nl_table_num_slots, nsoctot, nspins, slot_num
130 INTEGER(KIND=int_8) :: nl_table_key
131 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind
132 INTEGER,
DIMENSION(3) :: cell_b
133 INTEGER,
DIMENSION(:),
POINTER :: atom_list, list_a, list_b
134 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
135 LOGICAL :: dista, distb, found, is_entry_null, &
136 is_task_valid, my_tddft, paw_atom, &
138 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: a_matrix, p_matrix
139 REAL(
dp),
DIMENSION(3) :: rac, rbc
140 REAL(
dp),
DIMENSION(3, 3) :: force_tmp
141 REAL(kind=
dp) :: eps_cpc, factor1, factor2
142 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: c_int_h, c_int_s, coc
143 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dcpc_h, dcpc_s
144 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: pc_h, pc_s
145 REAL(kind=
dp),
DIMENSION(2) :: force_fac
146 REAL(kind=
dp),
DIMENSION(3, 3) :: pv_virial_thread
147 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: c_coeff_hh_a, c_coeff_hh_b, &
148 c_coeff_ss_a, c_coeff_ss_b
149 TYPE(
alist_type),
POINTER :: alist_ac, alist_bc
158 DIMENSION(:),
POINTER :: nl_iterator
165 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: my_kind_set
166 TYPE(
rho_atom_coeff),
DIMENSION(:),
POINTER :: int_local_h, int_local_s
174 CALL timeset(routinen, handle)
176 NULLIFY (my_kind_set, atomic_kind_set, force, my_oce, para_env, my_rho_atom, my_sab)
177 NULLIFY (mat_h, mat_p, dft_control)
180 qs_kind_set=my_kind_set, &
181 atomic_kind_set=atomic_kind_set, &
185 rho_atom_set=my_rho_atom, &
188 dft_control=dft_control)
190 nspins = dft_control%nspins
191 nimages = dft_control%nimages
198 IF (
PRESENT(tddft)) my_tddft = tddft
200 IF (nspins == 1) factor1 = 2.0_dp
201 cpassert(nimages == 1)
203 IF (
PRESENT(kscale))
THEN
204 factor1 = factor1*kscale
205 factor2 = factor2*kscale
207 IF (
PRESENT(kintegral)) factor1 = kintegral
208 IF (
PRESENT(kforce)) factor2 = kforce
210 IF (
PRESENT(fscale)) force_fac(:) = fscale(:)
212 IF (
PRESENT(rho_atom_external)) my_rho_atom => rho_atom_external
213 IF (
PRESENT(kind_set_external)) my_kind_set => kind_set_external
214 IF (
PRESENT(oce_external)) my_oce => oce_external
215 IF (
PRESENT(sab_external)) my_sab => sab_external
218 NULLIFY (cell_to_index)
219 IF (nimages > 1)
THEN
220 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
224 eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
227 CALL get_qs_kind_set(my_kind_set, maxsgf=max_nsgf, maxgtops=max_gau, basis_type=
"GAPW_1C")
232 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
237 pv_virial_thread(:, :) = 0.0_dp
239 nkind =
SIZE(my_kind_set, 1)
241 associate(mepos => para_env%mepos, num_pe => para_env%num_pe)
244 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
245 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom)
249 DO iat = bo(1), bo(2)
250 iatom = atom_list(iat)
252 CALL prj_gather(my_rho_atom(iatom)%ga_Vlocal_gb_h(ispin)%r_coef, &
253 my_rho_atom(iatom)%int_scr_h(ispin)%r_coef, my_kind_set(ikind))
254 CALL prj_gather(my_rho_atom(iatom)%ga_Vlocal_gb_s(ispin)%r_coef, &
255 my_rho_atom(iatom)%int_scr_s(ispin)%r_coef, my_kind_set(ikind))
259 DO ip = 0, num_pe - 1
261 DO iat = bo(1), bo(2)
262 iatom = atom_list(iat)
264 CALL para_env%bcast(my_rho_atom(iatom)%int_scr_h(ispin)%r_coef, ip)
265 CALL para_env%bcast(my_rho_atom(iatom)%int_scr_s(ispin)%r_coef, ip)
273 ALLOCATE (basis_set_list(nkind))
275 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_set_a)
276 IF (
ASSOCIATED(basis_set_a))
THEN
277 basis_set_list(ikind)%gto_basis_set => basis_set_a
279 NULLIFY (basis_set_list(ikind)%gto_basis_set)
286 nl_table_num_slots = natom*natom/2
293 IF (task%iatom <= task%jatom)
THEN
294 nl_table_key = natom*task%iatom + task%jatom
296 nl_table_key = natom*task%jatom + task%iatom
331 ALLOCATE (c_int_h(max_gau*max_nsgf), c_int_s(max_gau*max_nsgf), coc(max_gau*max_gau), &
332 a_matrix(max_gau, max_gau), p_matrix(max_nsgf, max_nsgf))
334 ALLOCATE (mat_h(nspins), mat_p(nspins))
336 NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
340 ALLOCATE (dcpc_h(max_gau, max_gau), dcpc_s(max_gau, max_gau), &
341 pc_h(max_nsgf, max_gau, nspins), pc_s(max_nsgf, max_gau, nspins))
356 DO slot_num = 1, nl_table_num_slots
359 IF (.NOT. is_entry_null)
THEN
362 is_task_valid =
ASSOCIATED(task)
363 DO WHILE (is_task_valid)
369 cell_b(:) = task%cell(:)
371 basis_set_a => basis_set_list(ikind)%gto_basis_set
372 IF (.NOT.
ASSOCIATED(basis_set_a)) cycle
373 basis_set_b => basis_set_list(jkind)%gto_basis_set
374 IF (.NOT.
ASSOCIATED(basis_set_b)) cycle
376 IF (nimages > 1)
THEN
377 img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
384 NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
387 IF (iatom <= jatom)
THEN
389 row=iatom, col=jatom, &
390 block=mat_h(ispin)%array, found=found)
393 row=jatom, col=iatom, &
394 block=mat_h(ispin)%array, found=found)
400 IF (iatom <= jatom)
THEN
402 row=iatom, col=jatom, &
403 block=mat_p(ispin)%array, found=found)
406 row=jatom, col=iatom, &
407 block=mat_p(ispin)%array, found=found)
414 CALL get_qs_kind(my_kind_set(kkind), paw_atom=paw_atom)
416 IF (.NOT. paw_atom) cycle
418 iac = ikind + nkind*(kkind - 1)
419 ibc = jkind + nkind*(kkind - 1)
420 IF (.NOT.
ASSOCIATED(my_oce%intac(iac)%alist)) cycle
421 IF (.NOT.
ASSOCIATED(my_oce%intac(ibc)%alist)) cycle
423 CALL get_alist(my_oce%intac(iac), alist_ac, iatom)
424 CALL get_alist(my_oce%intac(ibc), alist_bc, jatom)
425 IF (.NOT.
ASSOCIATED(alist_ac)) cycle
426 IF (.NOT.
ASSOCIATED(alist_bc)) cycle
428 DO kac = 1, alist_ac%nclist
429 DO kbc = 1, alist_bc%nclist
430 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
432 IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0))
THEN
433 n_cont_a = alist_ac%clist(kac)%nsgf_cnt
434 n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
435 IF (n_cont_a .EQ. 0 .OR. n_cont_b .EQ. 0) cycle
437 list_a => alist_ac%clist(kac)%sgf_list
438 list_b => alist_bc%clist(kbc)%sgf_list
440 katom = alist_ac%clist(kac)%catom
442 IF (iatom == katom .AND. all(alist_ac%clist(kac)%cell == 0))
THEN
443 c_coeff_hh_a => alist_ac%clist(kac)%achint
444 c_coeff_ss_a => alist_ac%clist(kac)%acint
447 c_coeff_hh_a => alist_ac%clist(kac)%acint
448 c_coeff_ss_a => alist_ac%clist(kac)%acint
452 IF (jatom == katom .AND. all(alist_bc%clist(kbc)%cell == 0))
THEN
453 c_coeff_hh_b => alist_bc%clist(kbc)%achint
454 c_coeff_ss_b => alist_bc%clist(kbc)%acint
457 c_coeff_hh_b => alist_bc%clist(kbc)%acint
458 c_coeff_ss_b => alist_bc%clist(kbc)%acint
462 rho_at => my_rho_atom(katom)
463 nsoctot =
SIZE(c_coeff_ss_a, 2)
465 CALL get_rho_atom(rho_atom=rho_at, int_scr_h=int_local_h, int_scr_s=int_local_s)
466 CALL add_vhxca_to_ks(mat_h, c_coeff_hh_a, c_coeff_hh_b, c_coeff_ss_a, c_coeff_ss_b, &
467 int_local_h, int_local_s, nspins, iatom, jatom, nsoctot, factor1, &
468 list_a, n_cont_a, list_b, n_cont_b, c_int_h, c_int_s, a_matrix, dista, distb, coc)
472 rac = alist_ac%clist(kac)%rac
473 rbc = alist_bc%clist(kbc)%rac
475 ia_kind = atom_of_kind(iatom)
476 ja_kind = atom_of_kind(jatom)
477 ka_kind = atom_of_kind(katom)
478 rho_at => my_rho_atom(katom)
479 force_tmp(1:3, 1:3) = 0.0_dp
481 CALL get_rho_atom(rho_atom=rho_at, int_scr_h=int_local_h, int_scr_s=int_local_s)
482 IF (iatom <= jatom)
THEN
483 CALL add_vhxca_forces(mat_p, c_coeff_hh_a, c_coeff_hh_b, c_coeff_ss_a, c_coeff_ss_b, &
484 int_local_h, int_local_s, force_tmp, nspins, iatom, jatom, nsoctot, &
485 list_a, n_cont_a, list_b, n_cont_b, dcpc_h, dcpc_s, ldcpc, &
486 pc_h, pc_s, p_matrix, force_fac)
487 force_tmp = factor2*force_tmp
489 force(kkind)%vhxc_atom(1:3, ka_kind) = force(kkind)%vhxc_atom(1:3, ka_kind) + force_tmp(1:3, 3)
492 force(ikind)%vhxc_atom(1:3, ia_kind) = force(ikind)%vhxc_atom(1:3, ia_kind) + force_tmp(1:3, 1)
495 force(jkind)%vhxc_atom(1:3, ja_kind) = force(jkind)%vhxc_atom(1:3, ja_kind) + force_tmp(1:3, 2)
502 CALL add_vhxca_forces(mat_p, c_coeff_hh_b, c_coeff_hh_a, c_coeff_ss_b, c_coeff_ss_a, &
503 int_local_h, int_local_s, force_tmp, nspins, jatom, iatom, nsoctot, &
504 list_b, n_cont_b, list_a, n_cont_a, dcpc_h, dcpc_s, ldcpc, &
505 pc_h, pc_s, p_matrix, force_fac)
506 force_tmp = factor2*force_tmp
508 force(kkind)%vhxc_atom(1:3, ka_kind) = force(kkind)%vhxc_atom(1:3, ka_kind) + force_tmp(1:3, 3)
511 force(ikind)%vhxc_atom(1:3, ia_kind) = force(ikind)%vhxc_atom(1:3, ia_kind) + force_tmp(1:3, 2)
514 force(jkind)%vhxc_atom(1:3, ja_kind) = force(jkind)%vhxc_atom(1:3, ja_kind) + force_tmp(1:3, 1)
530 next_task => task%next
533 is_task_valid =
ASSOCIATED(next_task)
534 IF (is_task_valid) task => next_task
545 NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
547 DEALLOCATE (mat_h, mat_p, c_int_h, c_int_s, a_matrix, p_matrix, coc)
550 DEALLOCATE (dcpc_h, dcpc_s, pc_h, pc_s)
567 virial%pv_gapw(1:3, 1:3) = virial%pv_gapw(1:3, 1:3) + pv_virial_thread(1:3, 1:3)
568 virial%pv_virial(1:3, 1:3) = virial%pv_virial(1:3, 1:3) + pv_virial_thread(1:3, 1:3)
573 DEALLOCATE (basis_set_list)
575 CALL timestop(handle)
604 SUBROUTINE add_vhxca_to_ks(mat_h, C_hh_a, C_hh_b, C_ss_a, C_ss_b, &
605 int_local_h, int_local_s, &
606 nspins, ia, ja, nsp, factor, lista, nconta, listb, ncontb, &
607 C_int_h, C_int_s, a_matrix, dista, distb, coc)
608 TYPE(cp_2d_r_p_type),
DIMENSION(:),
POINTER :: mat_h
609 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(IN) :: c_hh_a, c_hh_b, c_ss_a, c_ss_b
610 TYPE(rho_atom_coeff),
DIMENSION(:),
POINTER :: int_local_h, int_local_s
611 INTEGER,
INTENT(IN) :: nspins, ia, ja, nsp
612 REAL(kind=dp),
INTENT(IN) :: factor
613 INTEGER,
DIMENSION(:),
INTENT(IN) :: lista
614 INTEGER,
INTENT(IN) :: nconta
615 INTEGER,
DIMENSION(:),
INTENT(IN) :: listb
616 INTEGER,
INTENT(IN) :: ncontb
617 REAL(kind=dp),
DIMENSION(:),
INTENT(OUT) :: c_int_h, c_int_s
618 REAL(dp),
DIMENSION(:, :) :: a_matrix
619 LOGICAL,
INTENT(IN) :: dista, distb
620 REAL(dp),
DIMENSION(:) :: coc
622 INTEGER :: i, ispin, j, k
623 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: h_block, int_hard, int_soft
625 NULLIFY (int_hard, int_soft)
628 h_block => mat_h(ispin)%array
630 int_hard => int_local_h(ispin)%r_coef
631 int_soft => int_local_s(ispin)%r_coef
634 IF (dista .AND. distb)
THEN
639 coc(k) = int_hard(j, i) - int_soft(j, i)
642 CALL dgemm(
'N',
'T', nsp, ncontb, nsp, 1.0_dp, coc, nsp, c_hh_b(:, :, 1),
SIZE(c_hh_b, 1), &
643 0.0_dp, c_int_h, nsp)
644 CALL dgemm(
'N',
'N', nconta, ncontb, nsp, factor, c_hh_a(:, :, 1),
SIZE(c_hh_a, 1), &
645 c_int_h, nsp, 0.0_dp, a_matrix,
SIZE(a_matrix, 1))
647 CALL dgemm(
'N',
'T', nsp, ncontb, nsp, 1.0_dp, int_hard,
SIZE(int_hard, 1), &
648 c_hh_b(:, :, 1),
SIZE(c_hh_b, 1), 0.0_dp, coc, nsp)
649 CALL dgemm(
'N',
'T', nsp, ncontb, nsp, -1.0_dp, int_soft,
SIZE(int_soft, 1), &
650 c_ss_b(:, :, 1),
SIZE(c_ss_b, 1), 1.0_dp, coc, nsp)
651 CALL dgemm(
'N',
'N', nconta, ncontb, nsp, factor, c_hh_a(:, :, 1),
SIZE(c_hh_a, 1), &
652 coc, nsp, 0.0_dp, a_matrix,
SIZE(a_matrix, 1))
654 CALL dgemm(
'N',
'N', nconta, nsp, nsp, factor, c_hh_a(:, :, 1),
SIZE(c_hh_a, 1), &
655 int_hard,
SIZE(int_hard, 1), 0.0_dp, coc, nconta)
656 CALL dgemm(
'N',
'N', nconta, nsp, nsp, -factor, c_ss_a(:, :, 1),
SIZE(c_ss_a, 1), &
657 int_soft,
SIZE(int_soft, 1), 1.0_dp, coc, nconta)
658 CALL dgemm(
'N',
'T', nconta, ncontb, nsp, 1.0_dp, coc, nconta, &
659 c_hh_b(:, :, 1),
SIZE(c_hh_b, 1), 0.0_dp, a_matrix,
SIZE(a_matrix, 1))
661 CALL dgemm(
'N',
'T', nsp, ncontb, nsp, 1.0_dp, int_hard,
SIZE(int_hard, 1), &
662 c_hh_b(:, :, 1),
SIZE(c_hh_b, 1), &
663 0.0_dp, c_int_h, nsp)
664 CALL dgemm(
'N',
'T', nsp, ncontb, nsp, 1.0_dp, int_soft,
SIZE(int_soft, 1), &
665 c_ss_b(:, :, 1),
SIZE(c_ss_b, 1), &
666 0.0_dp, c_int_s, nsp)
667 CALL dgemm(
'N',
'N', nconta, ncontb, nsp, factor, c_hh_a(:, :, 1),
SIZE(c_hh_a, 1), &
669 0.0_dp, a_matrix,
SIZE(a_matrix, 1))
670 CALL dgemm(
'N',
'N', nconta, ncontb, nsp, -factor, c_ss_a(:, :, 1),
SIZE(c_ss_a, 1), &
672 1.0_dp, a_matrix,
SIZE(a_matrix, 1))
675 CALL alist_post_align_blk(a_matrix,
SIZE(a_matrix, 1), h_block,
SIZE(h_block, 1), &
676 lista, nconta, listb, ncontb)
678 IF (dista .AND. distb)
THEN
683 coc(k) = int_hard(j, i) - int_soft(j, i)
686 CALL dgemm(
'N',
'T', nsp, nconta, nsp, 1.0_dp, coc, nsp, c_hh_a(:, :, 1),
SIZE(c_hh_a, 1), 0.0_dp, c_int_h, nsp)
687 CALL dgemm(
'N',
'N', ncontb, nconta, nsp, factor, c_hh_b(:, :, 1),
SIZE(c_hh_b, 1), &
688 c_int_h, nsp, 0.0_dp, a_matrix,
SIZE(a_matrix, 1))
690 CALL dgemm(
'N',
'N', ncontb, nsp, nsp, factor, c_hh_b(:, :, 1),
SIZE(c_hh_b, 1), &
691 int_hard,
SIZE(int_hard, 1), 0.0_dp, coc, ncontb)
692 CALL dgemm(
'N',
'N', ncontb, nsp, nsp, -factor, c_ss_b(:, :, 1),
SIZE(c_ss_b, 1), &
693 int_soft,
SIZE(int_soft, 1), 1.0_dp, coc, ncontb)
694 CALL dgemm(
'N',
'T', ncontb, nconta, nsp, 1.0_dp, coc, ncontb, &
695 c_hh_a(:, :, 1),
SIZE(c_hh_a, 1), 0.0_dp, a_matrix,
SIZE(a_matrix, 1))
697 CALL dgemm(
'N',
'T', nsp, nconta, nsp, 1.0_dp, int_hard,
SIZE(int_hard, 1), &
698 c_hh_a(:, :, 1),
SIZE(c_hh_a, 1), 0.0_dp, coc, nsp)
699 CALL dgemm(
'N',
'T', nsp, nconta, nsp, -1.0_dp, int_soft,
SIZE(int_soft, 1), &
700 c_ss_a(:, :, 1),
SIZE(c_ss_a, 1), 1.0_dp, coc, nsp)
701 CALL dgemm(
'N',
'N', ncontb, nconta, nsp, factor, c_hh_b(:, :, 1),
SIZE(c_hh_b, 1), &
702 coc, nsp, 0.0_dp, a_matrix,
SIZE(a_matrix, 1))
704 CALL dgemm(
'N',
'T', nsp, nconta, nsp, 1.0_dp, int_hard,
SIZE(int_hard, 1), &
705 c_hh_a(:, :, 1),
SIZE(c_hh_a, 1), &
706 0.0_dp, c_int_h, nsp)
707 CALL dgemm(
'N',
'T', nsp, nconta, nsp, 1.0_dp, int_soft,
SIZE(int_soft, 1), &
708 c_ss_a(:, :, 1),
SIZE(c_ss_a, 1), &
709 0.0_dp, c_int_s, nsp)
710 CALL dgemm(
'N',
'N', ncontb, nconta, nsp, factor, c_hh_b(:, :, 1),
SIZE(c_hh_b, 1), &
712 0.0_dp, a_matrix,
SIZE(a_matrix, 1))
713 CALL dgemm(
'N',
'N', ncontb, nconta, nsp, -factor, c_ss_b(:, :, 1),
SIZE(c_ss_b, 1), &
715 1.0_dp, a_matrix,
SIZE(a_matrix, 1))
718 CALL alist_post_align_blk(a_matrix,
SIZE(a_matrix, 1), h_block,
SIZE(h_block, 1), &
719 listb, ncontb, lista, nconta)
723 END SUBROUTINE add_vhxca_to_ks
751 SUBROUTINE add_vhxca_forces(mat_p, C_hh_a, C_hh_b, C_ss_a, C_ss_b, &
752 int_local_h, int_local_s, &
753 force, nspins, ia, ja, nsp, lista, nconta, listb, ncontb, &
754 dCPC_h, dCPC_s, ldCPC, PC_h, PC_s, p_matrix, force_scaling)
755 TYPE(cp_2d_r_p_type),
DIMENSION(:),
INTENT(IN), &
757 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(IN) :: c_hh_a, c_hh_b, c_ss_a, c_ss_b
758 TYPE(rho_atom_coeff),
DIMENSION(:),
POINTER :: int_local_h, int_local_s
759 REAL(dp),
DIMENSION(3, 3),
INTENT(INOUT) :: force
760 INTEGER,
INTENT(IN) :: nspins, ia, ja, nsp
761 INTEGER,
DIMENSION(:),
INTENT(IN) :: lista
762 INTEGER,
INTENT(IN) :: nconta
763 INTEGER,
DIMENSION(:),
INTENT(IN) :: listb
764 INTEGER,
INTENT(IN) :: ncontb
765 REAL(kind=dp),
DIMENSION(:, :) :: dcpc_h, dcpc_s
766 INTEGER,
INTENT(IN) :: ldcpc
767 REAL(kind=dp),
DIMENSION(:, :, :) :: pc_h, pc_s
768 REAL(kind=dp),
DIMENSION(:, :) :: p_matrix
769 REAL(kind=dp),
DIMENSION(2),
INTENT(IN) :: force_scaling
771 INTEGER :: dir, ispin
772 REAL(dp),
DIMENSION(:, :),
POINTER :: int_hard, int_soft
773 REAL(kind=dp) :: ieqj, trace
774 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: p_block
780 IF (ia == ja) ieqj = 1.0_dp
782 NULLIFY (int_hard, int_soft)
785 p_block => mat_p(ispin)%array
787 CALL alist_pre_align_blk(p_block,
SIZE(p_block, 1), p_matrix,
SIZE(p_matrix, 1), &
788 lista, nconta, listb, ncontb)
790 int_hard => int_local_h(ispin)%r_coef
791 int_soft => int_local_s(ispin)%r_coef
793 CALL dgemm(
'N',
'N', nconta, nsp, ncontb, 1.0_dp, p_matrix,
SIZE(p_matrix, 1), &
794 c_hh_b(:, :, 1),
SIZE(c_hh_b, 1), &
795 0.0_dp, pc_h(:, :, ispin),
SIZE(pc_h, 1))
796 CALL dgemm(
'N',
'N', nconta, nsp, ncontb, 1.0_dp, p_matrix(:, :),
SIZE(p_matrix, 1), &
797 c_ss_b(:, :, 1),
SIZE(c_ss_b, 1), &
798 0.0_dp, pc_s(:, :, ispin),
SIZE(pc_s, 1))
801 CALL dgemm(
'T',
'N', nsp, nsp, nconta, 1.0_dp, pc_h(:, :, ispin),
SIZE(pc_h, 1), &
802 c_hh_a(:, :, dir),
SIZE(c_hh_a, 1), &
803 0.0_dp, dcpc_h,
SIZE(dcpc_h, 1))
804 trace = trace_r_axb(dcpc_h, ldcpc, int_hard, nsp, nsp, nsp)
805 force(dir - 1, 3) = force(dir - 1, 3) + ieqj*trace*force_scaling(ispin)
806 force(dir - 1, 1) = force(dir - 1, 1) - ieqj*trace*force_scaling(ispin)
808 CALL dgemm(
'T',
'N', nsp, nsp, nconta, 1.0_dp, pc_s(:, :, ispin),
SIZE(pc_s, 1), &
809 c_ss_a(:, :, dir),
SIZE(c_ss_a, 1), &
810 0.0_dp, dcpc_s,
SIZE(dcpc_s, 1))
811 trace = trace_r_axb(dcpc_s, ldcpc, int_soft, nsp, nsp, nsp)
812 force(dir - 1, 3) = force(dir - 1, 3) - ieqj*trace*force_scaling(ispin)
813 force(dir - 1, 1) = force(dir - 1, 1) + ieqj*trace*force_scaling(ispin)
817 CALL dgemm(
'T',
'N', ncontb, nsp, nconta, 1.0_dp, p_matrix,
SIZE(p_matrix, 1), &
818 c_hh_a(:, :, 1),
SIZE(c_hh_a, 1), &
819 0.0_dp, pc_h(:, :, ispin),
SIZE(pc_h, 1))
820 CALL dgemm(
'T',
'N', ncontb, nsp, nconta, 1.0_dp, p_matrix,
SIZE(p_matrix, 1), &
821 c_ss_a(:, :, 1),
SIZE(c_ss_a, 1), &
822 0.0_dp, pc_s(:, :, ispin),
SIZE(pc_s, 1))
825 CALL dgemm(
'T',
'N', nsp, nsp, ncontb, 1.0_dp, pc_h(:, :, ispin),
SIZE(pc_h, 1), &
826 c_hh_b(:, :, dir),
SIZE(c_hh_b, 1), &
827 0.0_dp, dcpc_h,
SIZE(dcpc_h, 1))
828 trace = trace_r_axb(dcpc_h, ldcpc, int_hard, nsp, nsp, nsp)
829 force(dir - 1, 3) = force(dir - 1, 3) + ieqj*trace*force_scaling(ispin)
830 force(dir - 1, 2) = force(dir - 1, 2) - ieqj*trace*force_scaling(ispin)
832 CALL dgemm(
'T',
'N', nsp, nsp, ncontb, 1.0_dp, pc_s(:, :, ispin),
SIZE(pc_s, 1), &
833 c_ss_b(:, :, dir),
SIZE(c_ss_b, 1), &
834 0.0_dp, dcpc_s,
SIZE(dcpc_s, 1))
835 trace = trace_r_axb(dcpc_s, ldcpc, int_soft, nsp, nsp, nsp)
836 force(dir - 1, 3) = force(dir - 1, 3) - ieqj*trace*force_scaling(ispin)
837 force(dir - 1, 2) = force(dir - 1, 2) + ieqj*trace*force_scaling(ispin)
842 END SUBROUTINE add_vhxca_forces
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
All kind of helpful little routines.
pure real(dp) function, public trace_r_axb(a, lda, b, ldb, m, n)
...
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.
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Defines the basic variable types.
integer, parameter, public int_8
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.
Interface to the message passing library MPI.
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_pp, 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, 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)
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, 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, 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)
Get attributes of an atomic kind set.
routines that build the Kohn-Sham matrix contributions coming from local atomic densities
subroutine, public update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, kind_set_external, oce_external, sab_external, kscale, kintegral, kforce, fscale)
The correction to the KS matrix due to the GAPW local terms to the hartree and XC contributions is he...
Define the neighbor list data types and the corresponding functionality.
subroutine, public get_iterator_task(iterator_set, task, mepos)
Captures the current state of the iterator in a neighbor_list_task_type.
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)
...
A simple hash table of integer keys, using hash function: H(k) = (k*p) mod n + 1 where: k = key p = a...
subroutine, public nl_hash_table_is_null(hash_table, key, is_null)
Initialises a nl_hash_table object.
subroutine, public nl_hash_table_release(hash_table)
releases the hash table. Note that deallocating tasks stored in the table is the responsibility of th...
subroutine, public nl_hash_table_status(hash_table, nelements, nmax, prime)
outputs the current information about the table
recursive subroutine, public nl_hash_table_add(hash_table, key, val)
Add element to a hash table, auto resize if necessary.
subroutine, public nl_hash_table_get_from_index(hash_table, idx, val)
Retrieve value from a hash table given a specified index.
subroutine, public nl_hash_table_create(hash_table, nmax)
Creates and initialises an empty nl_hash_table object.
Routines for the construction of the coefficients for the expansion of the atomic densities rho1_hard...
subroutine, public prj_gather(ain, aout, atom)
...
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_vlocal_gb_h, ga_vlocal_gb_s, int_scr_h, int_scr_s)
...
General overlap type integrals containers.
subroutine, public alist_post_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
...
subroutine, public alist_pre_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
...
subroutine, public get_alist(sap_int, alist, atom)
...
All kind of helpful little routines.
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
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.
represent a pointer to a 2d array
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.
the object container which allows for the creation of an array of pointers to nl_hash_table objects