20 gto_basis_set_p_type,&
54 USE dbcsr_api,
ONLY: &
55 dbcsr_add_block_node, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, &
56 dbcsr_deallocate_matrix, dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_block_p, &
57 dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry
83 pw_integrate_function,&
110 neighbor_list_iterator_p_type,&
112 neighbor_list_set_p_type
119 realspace_grid_desc_type,&
120 realspace_grid_type,&
132 #include "./base/base_uses.f90"
141 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_linres_current'
145 REAL(dp),
POINTER,
DIMENSION(:, :) :: r
147 REAL(dp),
DIMENSION(3, 3, 3),
PARAMETER :: Levi_Civita = reshape((/ &
148 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
149 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
150 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), (/3, 3, 3/))
181 TYPE(current_env_type) :: current_env
182 TYPE(qs_environment_type),
POINTER :: qs_env
183 INTEGER,
INTENT(IN) :: ib
185 CHARACTER(LEN=*),
PARAMETER :: routinen =
'current_build_current'
187 CHARACTER(LEN=default_path_length) :: ext, filename, my_pos
188 INTEGER :: handle, idir, iib, iiib, ispin, istate, &
189 j, jstate, nao, natom, nmo, nspins, &
190 nstates(2), output_unit, unit_nr
191 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf
192 INTEGER,
DIMENSION(:),
POINTER :: row_blk_sizes
193 LOGICAL :: append_cube, gapw, mpi_io
194 REAL(dp) :: dk(3), jrho_tot_g(3, 3), &
195 jrho_tot_r(3, 3), maxocc, scale_fac
196 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :) :: ddk
197 REAL(dp),
EXTERNAL :: ddot
198 TYPE(cell_type),
POINTER :: cell
199 TYPE(cp_2d_i_p_type),
DIMENSION(:),
POINTER :: center_list
200 TYPE(cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: p_psi1, psi1
201 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: psi0_order
202 TYPE(cp_fm_type),
DIMENSION(:, :),
POINTER :: psi1_d, psi1_p, psi1_rxp
203 TYPE(cp_fm_type),
POINTER :: mo_coeff
204 TYPE(cp_logger_type),
POINTER :: logger
205 TYPE(dbcsr_distribution_type),
POINTER :: dbcsr_dist
206 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: density_matrix0, density_matrix_a, &
207 density_matrix_ii, density_matrix_iii
208 TYPE(dft_control_type),
POINTER :: dft_control
209 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
210 TYPE(mp_para_env_type),
POINTER :: para_env
211 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
213 TYPE(particle_list_type),
POINTER :: particles
214 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
215 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: jrho1_g
216 TYPE(pw_env_type),
POINTER :: pw_env
217 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
218 TYPE(pw_r3d_rs_type) :: wf_r
219 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: jrho1_r
220 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
221 TYPE(qs_matrix_pools_type),
POINTER :: mpools
222 TYPE(qs_subsys_type),
POINTER :: subsys
223 TYPE(realspace_grid_desc_type),
POINTER :: auxbas_rs_desc
224 TYPE(section_vals_type),
POINTER :: current_section
226 CALL timeset(routinen, handle)
228 NULLIFY (logger, current_section, density_matrix0, density_matrix_a, &
229 density_matrix_ii, density_matrix_iii, cell, dft_control, mos, &
230 particle_set, pw_env, auxbas_rs_desc, auxbas_pw_pool, &
231 para_env, center_list, mo_coeff, jrho1_r, jrho1_g, &
232 psi1_p, psi1_d, psi1_rxp, sab_all, qs_kind_set)
239 center_list=center_list, &
243 psi0_order=psi0_order, &
250 dft_control=dft_control, &
257 particle_set=particle_set, &
258 qs_kind_set=qs_kind_set, &
259 dbcsr_dist=dbcsr_dist)
263 gapw = dft_control%qs_control%gapw
264 nspins = dft_control%nspins
265 natom =
SIZE(particle_set, 1)
268 ALLOCATE (psi1(nspins), p_psi1(nspins))
270 CALL cp_fm_create(psi1(ispin), psi0_order(ispin)%matrix_struct)
271 CALL cp_fm_create(p_psi1(ispin), psi0_order(ispin)%matrix_struct)
283 ALLOCATE (first_sgf(natom))
284 ALLOCATE (last_sgf(natom))
286 first_sgf=first_sgf, &
288 ALLOCATE (row_blk_sizes(natom))
289 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
290 DEALLOCATE (first_sgf)
291 DEALLOCATE (last_sgf)
297 ALLOCATE (density_matrix0(ispin)%matrix)
298 CALL dbcsr_create(matrix=density_matrix0(ispin)%matrix, &
299 name=
"density_matrix0", &
300 dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
301 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
302 nze=0, mutable_work=.true.)
306 ALLOCATE (density_matrix_a(ispin)%matrix)
307 CALL dbcsr_copy(density_matrix_a(ispin)%matrix, density_matrix0(ispin)%matrix, &
308 name=
"density_matrix_a")
311 ALLOCATE (density_matrix_ii(ispin)%matrix)
312 CALL dbcsr_copy(density_matrix_ii(ispin)%matrix, density_matrix0(ispin)%matrix, &
313 name=
"density_matrix_ii")
316 ALLOCATE (density_matrix_iii(ispin)%matrix)
317 CALL dbcsr_copy(density_matrix_iii(ispin)%matrix, density_matrix0(ispin)%matrix, &
318 name=
"density_matrix_iii")
321 DEALLOCATE (row_blk_sizes)
334 mo_coeff => psi0_order(ispin)
337 CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc)
341 CALL dbcsr_set(density_matrix0(ispin)%matrix, 0.0_dp)
343 matrix_v=mo_coeff, matrix_g=mo_coeff, &
344 ncol=nmo, alpha=maxocc)
347 ALLOCATE (ddk(3, nmo))
352 associate(psi_a_ib => psi1(ispin), psi_buf => p_psi1(ispin))
358 DO istate = 1, current_env%nbr_center(ispin)
359 dk(1:3) = current_env%centers_set(ispin)%array(1:3, istate)
363 DO j = center_list(ispin)%array(1, istate), center_list(ispin)%array(1, istate + 1) - 1
364 jstate = center_list(ispin)%array(2, j)
365 CALL cp_fm_to_fm(psi1_p(ispin, iib), psi_a_ib, 1, jstate, jstate)
366 CALL cp_fm_to_fm(psi1_p(ispin, iiib), psi_buf, 1, jstate, jstate)
367 ddk(:, jstate) = dk(1:3)
380 IF (current_env%full)
THEN
389 CALL dbcsr_set(density_matrix_a(ispin)%matrix, 0.0_dp)
391 matrix_v=mo_coeff, matrix_g=psi_a_ib, &
392 ncol=nmo, alpha=maxocc)
396 CALL dbcsr_set(density_matrix_iii(ispin)%matrix, 0.0_dp)
398 matrix_v=mo_coeff, matrix_g=psi1_p(ispin, iiib), &
399 ncol=nmo, alpha=maxocc)
402 CALL dbcsr_set(density_matrix_ii(ispin)%matrix, 0.0_dp)
404 matrix_v=mo_coeff, matrix_g=psi1_p(ispin, iib), &
405 ncol=nmo, alpha=maxocc)
412 CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_r=jrho1_r)
413 CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_g=jrho1_g)
414 associate(jrho_rspace => jrho1_r(ispin), jrho_gspace => jrho1_g(ispin))
415 CALL pw_zero(jrho_rspace)
416 CALL pw_zero(jrho_gspace)
418 density_matrix_a(ispin)%matrix, &
419 density_matrix_ii(ispin)%matrix, &
420 density_matrix_iii(ispin)%matrix, &
421 ib, idir, jrho_rspace, jrho_gspace, qs_env, &
424 scale_fac = cell%deth/
twopi
425 CALL pw_scale(jrho_rspace, scale_fac)
426 CALL pw_scale(jrho_gspace, scale_fac)
428 jrho_tot_g(idir, ib) = pw_integrate_function(jrho_gspace, isign=-1)
429 jrho_tot_r(idir, ib) = pw_integrate_function(jrho_rspace, isign=-1)
432 IF (output_unit > 0)
THEN
433 WRITE (output_unit,
'(T2,2(A,E24.16))')
'Integrated j_'&
434 &//achar(idir + 119)//achar(ib + 119)//
'(r): G-space=', &
435 jrho_tot_g(idir, ib),
' R-space=', jrho_tot_r(idir, ib)
451 CALL calculate_jrho_atom_coeff(qs_env, current_env, &
455 density_matrix_iii, &
462 CALL calculate_jrho_atom_rad(qs_env, current_env, idir)
465 CALL calculate_jrho_atom(current_env, qs_env, ib, idir)
470 IF (btest(cp_print_key_should_output(logger%iter_info, current_section,&
471 &
"PRINT%CURRENT_CUBES"), cp_p_file))
THEN
472 append_cube = section_get_lval(current_section,
"PRINT%CURRENT_CUBES%APPEND")
474 IF (append_cube)
THEN
478 CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
479 auxbas_pw_pool=auxbas_pw_pool)
481 CALL auxbas_pw_pool%create_pw(wf_r)
485 CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_r=jrho1_r)
487 CALL pw_axpy(jrho1_r(ispin), wf_r, 1.0_dp)
494 cpabort(
"GAPW needs to be finalized")
498 WRITE (ext,
'(a2,I1,a2,I1,a5)')
"iB", ib,
"_d", idir,
".cube"
499 WRITE (ext,
'(a2,a1,a2,a1,a5)')
"iB", achar(ib + 119),
"_d", achar(idir + 119),
".cube"
500 unit_nr = cp_print_key_unit_nr(logger, current_section,
"PRINT%CURRENT_CUBES", &
501 extension=trim(ext), middle_name=trim(filename), &
502 log_filename=.false., file_position=my_pos, &
505 CALL cp_pw_to_cube(wf_r, unit_nr,
"RESPONSE CURRENT DENSITY ", &
506 particles=particles, &
507 stride=section_get_ivals(current_section,
"PRINT%CURRENT_CUBES%STRIDE"), &
509 CALL cp_print_key_finished_output(unit_nr, logger, current_section,&
510 &
"PRINT%CURRENT_CUBES", mpi_io=mpi_io)
513 CALL auxbas_pw_pool%give_back_pw(wf_r)
517 IF (output_unit > 0)
THEN
518 WRITE (output_unit,
'(T2,A,E24.16)')
'CheckSum R-integrated j=', &
519 sqrt(ddot(9, jrho_tot_r(1, 1), 1, jrho_tot_r(1, 1), 1))
524 CALL cp_fm_release(psi1)
525 CALL cp_fm_release(p_psi1)
527 CALL dbcsr_deallocate_matrix_set(density_matrix0)
528 CALL dbcsr_deallocate_matrix_set(density_matrix_a)
529 CALL dbcsr_deallocate_matrix_set(density_matrix_ii)
530 CALL dbcsr_deallocate_matrix_set(density_matrix_iii)
533 CALL timestop(handle)
577 current_rs, current_gs, qs_env, current_env, soft_valid, retain_rsgrid)
579 TYPE(dbcsr_type),
POINTER :: mat_d0, mat_jp, mat_jp_rii, mat_jp_riii
580 INTEGER,
INTENT(IN) :: ib, idir
581 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: current_rs
582 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: current_gs
583 TYPE(qs_environment_type),
POINTER :: qs_env
584 TYPE(current_env_type) :: current_env
585 LOGICAL,
INTENT(IN),
OPTIONAL :: soft_valid, retain_rsgrid
587 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_jrho_resp'
588 INTEGER,
PARAMETER :: max_tasks = 2000
590 CHARACTER(LEN=default_string_length) :: basis_type
591 INTEGER :: adbmdab_func, bcol, brow, cindex, curr_tasks, handle, i, iatom, iatom_old, idir2, &
592 igrid_level, iib, iiib, ikind, ikind_old, ipgf, iset, iset_old, itask, ithread, jatom, &
593 jatom_old, jkind, jkind_old, jpgf, jset, jset_old, maxco, maxpgf, maxset, maxsgf, &
594 maxsgf_set, na1, na2, natom, nb1, nb2, ncoa, ncob, nimages, nkind, nseta, nsetb, ntasks, &
596 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
598 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa, first_sgfb
599 LOGICAL :: atom_pair_changed, den_found, &
600 den_found_a, distributed_rs_grids, &
601 do_igaim, my_retain_rsgrid, my_soft
602 REAL(dp),
DIMENSION(:, :, :),
POINTER :: my_current, my_gauge, my_rho
603 REAL(kind=dp) :: eps_rho_rspace, f, kind_radius_a, &
604 kind_radius_b, lxo2, lyo2, lzo2, &
605 prefactor, radius, scale, scale2, zetp
606 REAL(kind=dp),
DIMENSION(3) :: ra, rab, rb, rp
607 REAL(kind=dp),
DIMENSION(:),
POINTER :: set_radius_a, set_radius_b
608 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: jp_block_a, jp_block_b, jp_block_c, jp_block_d, &
609 jpab_a, jpab_b, jpab_c, jpab_d, jpblock_a, jpblock_b, jpblock_c, jpblock_d, rpgfa, rpgfb, &
610 sphi_a, sphi_b, work, zeta, zetb
611 REAL(kind=dp),
DIMENSION(:, :, :),
POINTER :: jpabt_a, jpabt_b, jpabt_c, jpabt_d, workt
612 TYPE(atom_pair_type),
DIMENSION(:),
POINTER :: atom_pair_recv, atom_pair_send
613 TYPE(cell_type),
POINTER :: cell
614 TYPE(cube_info_type),
DIMENSION(:),
POINTER :: cube_info
615 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: deltajp_a, deltajp_b, deltajp_c, &
617 TYPE(dbcsr_type),
POINTER :: mat_a, mat_b, mat_c, mat_d
618 TYPE(dft_control_type),
POINTER :: dft_control
619 TYPE(gridlevel_info_type),
POINTER :: gridlevel_info
620 TYPE(gto_basis_set_p_type),
DIMENSION(:),
POINTER :: basis_set_list
621 TYPE(gto_basis_set_type),
POINTER :: basis_set_a, basis_set_b, orb_basis_set
622 TYPE(mp_para_env_type),
POINTER :: para_env
623 TYPE(neighbor_list_iterator_p_type), &
624 DIMENSION(:),
POINTER :: nl_iterator
625 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
627 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
628 TYPE(pw_env_type),
POINTER :: pw_env
629 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
630 TYPE(qs_kind_type),
POINTER :: qs_kind
631 TYPE(realspace_grid_desc_p_type),
DIMENSION(:), &
633 TYPE(realspace_grid_type),
DIMENSION(:),
POINTER :: rs_current, rs_rho
634 TYPE(realspace_grid_type),
DIMENSION(:, :), &
636 TYPE(task_type),
DIMENSION(:),
POINTER :: tasks
638 NULLIFY (qs_kind, cell, dft_control, orb_basis_set, rs_rho, &
639 qs_kind_set, sab_orb, particle_set, rs_current, pw_env, &
640 rs_descs, para_env, set_radius_a, set_radius_b, la_max, la_min, &
641 lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, rpgfa, rpgfb, &
642 sphi_a, sphi_b, zeta, zetb, first_sgfa, first_sgfb, tasks, &
643 workt, mat_a, mat_b, mat_c, mat_d, rs_gauge)
644 NULLIFY (deltajp_a, deltajp_b, deltajp_c, deltajp_d)
645 NULLIFY (jp_block_a, jp_block_b, jp_block_c, jp_block_d)
646 NULLIFY (jpblock_a, jpblock_b, jpblock_c, jpblock_d)
647 NULLIFY (jpabt_a, jpabt_b, jpabt_c, jpabt_d)
648 NULLIFY (atom_pair_send, atom_pair_recv)
650 CALL timeset(routinen, handle)
655 do_igaim = current_env%gauge .EQ. current_gauge_atom
660 IF (do_igaim) mat_d => mat_d0
662 my_retain_rsgrid = .false.
663 IF (
PRESENT(retain_rsgrid)) my_retain_rsgrid = retain_rsgrid
665 CALL get_qs_env(qs_env=qs_env, &
666 qs_kind_set=qs_kind_set, &
668 dft_control=dft_control, &
669 particle_set=particle_set, &
674 IF (do_igaim)
CALL get_current_env(current_env=current_env, rs_gauge=rs_gauge)
677 CALL set_vecp(ib, iib, iiib)
682 IF (idir .NE. ib)
THEN
683 CALL set_vecp_rev(idir, ib, idir2)
684 scale2 = fac_vecp(idir, ib, idir2)
688 gridlevel_info => pw_env%gridlevel_info
689 cube_info => pw_env%cube_info
692 cpassert(
ASSOCIATED(sab_orb))
694 cpassert(
ASSOCIATED(pw_env))
695 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
697 distributed_rs_grids = .false.
698 DO igrid_level = 1, gridlevel_info%ngrid_levels
699 IF (.NOT. all(rs_descs(igrid_level)%rs_desc%perd == 1))
THEN
700 distributed_rs_grids = .true.
703 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
707 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
710 maxsgf_set=maxsgf_set)
712 lxo2 = sqrt(sum(cell%hmat(:, 1)**2))/2.0_dp
713 lyo2 = sqrt(sum(cell%hmat(:, 2)**2))/2.0_dp
714 lzo2 = sqrt(sum(cell%hmat(:, 3)**2))/2.0_dp
717 IF (
PRESENT(soft_valid)) my_soft = soft_valid
719 basis_type =
"ORB_SOFT"
724 nkind =
SIZE(qs_kind_set)
726 CALL reallocate(jpabt_a, 1, maxco, 1, maxco, 0, nthread - 1)
727 CALL reallocate(jpabt_b, 1, maxco, 1, maxco, 0, nthread - 1)
728 CALL reallocate(jpabt_c, 1, maxco, 1, maxco, 0, nthread - 1)
729 CALL reallocate(jpabt_d, 1, maxco, 1, maxco, 0, nthread - 1)
730 CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
731 CALL reallocate_tasks(tasks, max_tasks)
734 curr_tasks =
SIZE(tasks)
737 natom =
SIZE(particle_set)
742 nimages = dft_control%nimages
743 cpassert(nimages == 1)
747 qs_kind => qs_kind_set(ikind)
749 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_set)
751 IF (.NOT.
ASSOCIATED(orb_basis_set)) cycle
753 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, npgf=npgfa, nset=nseta)
754 maxset = max(nseta, maxset)
755 maxpgf = max(maxval(npgfa), maxpgf)
762 ALLOCATE (deltajp_a(1), deltajp_b(1), deltajp_c(1), deltajp_d(1))
763 IF (distributed_rs_grids)
THEN
764 ALLOCATE (deltajp_a(1)%matrix, deltajp_b(1)%matrix, deltajp_c(1)%matrix)
766 ALLOCATE (deltajp_d(1)%matrix)
769 CALL dbcsr_create(deltajp_a(1)%matrix, template=mat_a, name=
'deltajp_a')
770 CALL dbcsr_create(deltajp_b(1)%matrix, template=mat_a, name=
'deltajp_b')
771 CALL dbcsr_create(deltajp_c(1)%matrix, template=mat_a, name=
'deltajp_c')
772 IF (do_igaim)
CALL dbcsr_create(deltajp_d(1)%matrix, template=mat_a, name=
'deltajp_d')
774 deltajp_a(1)%matrix => mat_a
775 deltajp_b(1)%matrix => mat_b
776 deltajp_c(1)%matrix => mat_c
777 IF (do_igaim) deltajp_d(1)%matrix => mat_d
780 ALLOCATE (basis_set_list(nkind))
782 qs_kind => qs_kind_set(ikind)
783 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
784 IF (
ASSOCIATED(basis_set_a))
THEN
785 basis_set_list(ikind)%gto_basis_set => basis_set_a
787 NULLIFY (basis_set_list(ikind)%gto_basis_set)
790 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
791 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
792 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rab)
793 basis_set_a => basis_set_list(ikind)%gto_basis_set
794 IF (.NOT.
ASSOCIATED(basis_set_a)) cycle
795 basis_set_b => basis_set_list(jkind)%gto_basis_set
796 IF (.NOT.
ASSOCIATED(basis_set_b)) cycle
797 ra(:) =
pbc(particle_set(iatom)%r, cell)
799 first_sgfa => basis_set_a%first_sgf
800 la_max => basis_set_a%lmax
801 la_min => basis_set_a%lmin
802 npgfa => basis_set_a%npgf
803 nseta = basis_set_a%nset
804 nsgfa => basis_set_a%nsgf_set
805 rpgfa => basis_set_a%pgf_radius
806 set_radius_a => basis_set_a%set_radius
807 kind_radius_a = basis_set_a%kind_radius
808 sphi_a => basis_set_a%sphi
809 zeta => basis_set_a%zet
811 first_sgfb => basis_set_b%first_sgf
812 lb_max => basis_set_b%lmax
813 lb_min => basis_set_b%lmin
814 npgfb => basis_set_b%npgf
815 nsetb = basis_set_b%nset
816 nsgfb => basis_set_b%nsgf_set
817 rpgfb => basis_set_b%pgf_radius
818 set_radius_b => basis_set_b%set_radius
819 kind_radius_b = basis_set_b%kind_radius
820 sphi_b => basis_set_b%sphi
821 zetb => basis_set_b%zet
823 IF (abs(rab(1)) > lxo2 .OR. abs(rab(2)) > lyo2 .OR. abs(rab(3)) > lzo2)
THEN
830 CALL dbcsr_get_block_p(matrix=mat_a, row=brow, col=bcol, &
831 block=jp_block_a, found=den_found_a)
832 CALL dbcsr_get_block_p(matrix=mat_b, row=brow, col=bcol, &
833 block=jp_block_b, found=den_found)
834 CALL dbcsr_get_block_p(matrix=mat_c, row=brow, col=bcol, &
835 block=jp_block_c, found=den_found)
836 IF (do_igaim)
CALL dbcsr_get_block_p(matrix=mat_d, row=brow, col=bcol, &
837 block=jp_block_d, found=den_found)
839 IF (.NOT.
ASSOCIATED(jp_block_a)) cycle
841 IF (distributed_rs_grids)
THEN
842 NULLIFY (jpblock_a, jpblock_b, jpblock_c, jpblock_d)
843 CALL dbcsr_add_block_node(deltajp_a(1)%matrix, brow, bcol, jpblock_a)
844 jpblock_a = jp_block_a
845 CALL dbcsr_add_block_node(deltajp_b(1)%matrix, brow, bcol, jpblock_b)
846 jpblock_b = jp_block_b
847 CALL dbcsr_add_block_node(deltajp_c(1)%matrix, brow, bcol, jpblock_c)
848 jpblock_c = jp_block_c
850 CALL dbcsr_add_block_node(deltajp_d(1)%matrix, brow, bcol, jpblock_d)
851 jpblock_d = jp_block_d
854 jpblock_a => jp_block_a
855 jpblock_b => jp_block_b
856 jpblock_c => jp_block_c
857 IF (do_igaim) jpblock_d => jp_block_d
860 CALL task_list_inner_loop(tasks, ntasks, curr_tasks, rs_descs, &
861 dft_control, cube_info, gridlevel_info, cindex, &
862 iatom, jatom, rpgfa, rpgfb, zeta, zetb, kind_radius_b, &
863 set_radius_a, set_radius_b, ra, rab, &
864 la_max, la_min, lb_max, lb_min, npgfa, npgfb, nseta, nsetb)
867 CALL neighbor_list_iterator_release(nl_iterator)
869 DEALLOCATE (basis_set_list)
871 IF (distributed_rs_grids)
THEN
872 CALL dbcsr_finalize(deltajp_a(1)%matrix)
873 CALL dbcsr_finalize(deltajp_b(1)%matrix)
874 CALL dbcsr_finalize(deltajp_c(1)%matrix)
875 IF (do_igaim)
CALL dbcsr_finalize(deltajp_d(1)%matrix)
879 CALL distribute_tasks(rs_descs=rs_descs, ntasks=ntasks, natoms=natom, tasks=tasks, &
880 atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
881 symmetric=.false., reorder_rs_grid_ranks=.true., &
882 skip_load_balance_distributed=.false.)
884 ALLOCATE (rs_current(gridlevel_info%ngrid_levels))
886 DO igrid_level = 1, gridlevel_info%ngrid_levels
889 IF (rs_descs(igrid_level)%rs_desc%distributed .AND. .NOT. my_retain_rsgrid)
THEN
890 CALL rs_grid_release(rs_rho(igrid_level))
891 CALL rs_grid_create(rs_rho(igrid_level), rs_descs(igrid_level)%rs_desc)
893 CALL rs_grid_zero(rs_rho(igrid_level))
894 CALL rs_grid_create(rs_current(igrid_level), rs_descs(igrid_level)%rs_desc)
895 CALL rs_grid_zero(rs_current(igrid_level))
900 IF (.NOT. current_env%gauge_init .AND. do_igaim)
THEN
901 CALL current_set_gauge(current_env, qs_env)
902 current_env%gauge_init = .true.
907 DO igrid_level = 1, gridlevel_info%ngrid_levels
908 my_rho => rs_rho(igrid_level)%r
909 my_current => rs_current(igrid_level)%r
910 IF (lbound(my_rho, 3) .NE. lbound(my_current, 3) .OR. &
911 lbound(my_rho, 2) .NE. lbound(my_current, 2) .OR. &
912 lbound(my_rho, 1) .NE. lbound(my_current, 1) .OR. &
913 ubound(my_rho, 3) .NE. ubound(my_current, 3) .OR. &
914 ubound(my_rho, 2) .NE. ubound(my_current, 2) .OR. &
915 ubound(my_rho, 1) .NE. ubound(my_current, 1))
THEN
916 WRITE (*, *)
'LBOUND(my_rho,3),LBOUND(my_current,3)', lbound(my_rho, 3), lbound(my_current, 3)
917 WRITE (*, *)
'LBOUND(my_rho,2),LBOUND(my_current,2)', lbound(my_rho, 2), lbound(my_current, 2)
918 WRITE (*, *)
'LBOUND(my_rho,1),LBOUND(my_current,1)', lbound(my_rho, 1), lbound(my_current, 1)
919 WRITE (*, *)
'UBOUND(my_rho,3),UBOUND(my_current,3)', ubound(my_rho, 3), ubound(my_current, 3)
920 WRITE (*, *)
'UBOUND(my_rho,2),UBOUND(my_current,2)', ubound(my_rho, 2), ubound(my_current, 2)
921 WRITE (*, *)
'UBOUND(my_rho,1),UBOUND(my_current,1)', ubound(my_rho, 1), ubound(my_current, 1)
925 my_gauge => rs_gauge(1, igrid_level)%r
926 IF (lbound(my_rho, 3) .NE. lbound(my_gauge, 3) .OR. &
927 lbound(my_rho, 2) .NE. lbound(my_gauge, 2) .OR. &
928 lbound(my_rho, 1) .NE. lbound(my_gauge, 1) .OR. &
929 ubound(my_rho, 3) .NE. ubound(my_gauge, 3) .OR. &
930 ubound(my_rho, 2) .NE. ubound(my_gauge, 2) .OR. &
931 ubound(my_rho, 1) .NE. ubound(my_gauge, 1))
THEN
932 WRITE (*, *)
'LBOUND(my_rho,3),LBOUND(my_gauge,3)', lbound(my_rho, 3), lbound(my_gauge, 3)
933 WRITE (*, *)
'LBOUND(my_rho,2),LBOUND(my_gauge,2)', lbound(my_rho, 2), lbound(my_gauge, 2)
934 WRITE (*, *)
'LBOUND(my_rho,1),LBOUND(my_gauge,1)', lbound(my_rho, 1), lbound(my_gauge, 1)
935 WRITE (*, *)
'UBOUND(my_rho,3),UbOUND(my_gauge,3)', ubound(my_rho, 3), ubound(my_gauge, 3)
936 WRITE (*, *)
'UBOUND(my_rho,2),UBOUND(my_gauge,2)', ubound(my_rho, 2), ubound(my_gauge, 2)
937 WRITE (*, *)
'UBOUND(my_rho,1),UBOUND(my_gauge,1)', ubound(my_rho, 1), ubound(my_gauge, 1)
945 IF (distributed_rs_grids)
THEN
946 CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_a, &
947 atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
948 nimages=nimages, scatter=.true.)
949 CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_b, &
950 atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
951 nimages=nimages, scatter=.true.)
952 CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_c, &
953 atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
954 nimages=nimages, scatter=.true.)
955 IF (do_igaim)
CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_d, &
956 atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
957 nimages=nimages, scatter=.true.)
961 jpab_a => jpabt_a(:, :, ithread)
962 jpab_b => jpabt_b(:, :, ithread)
963 jpab_c => jpabt_c(:, :, ithread)
964 IF (do_igaim) jpab_d => jpabt_d(:, :, ithread)
965 work => workt(:, :, ithread)
967 iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
968 ikind_old = -1; jkind_old = -1
970 loop_tasks:
DO itask = 1, ntasks
971 igrid_level = tasks(itask)%grid_level
972 cindex = tasks(itask)%image
973 iatom = tasks(itask)%iatom
974 jatom = tasks(itask)%jatom
975 iset = tasks(itask)%iset
976 jset = tasks(itask)%jset
977 ipgf = tasks(itask)%ipgf
978 jpgf = tasks(itask)%jpgf
981 cpassert(tasks(itask)%dist_type .NE. 2)
983 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old)
THEN
985 ikind = particle_set(iatom)%atomic_kind%kind_number
986 jkind = particle_set(jatom)%atomic_kind%kind_number
988 IF (iatom .NE. iatom_old) ra(:) =
pbc(particle_set(iatom)%r, cell)
993 IF (ikind .NE. ikind_old)
THEN
994 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
995 basis_type=basis_type)
997 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
998 first_sgf=first_sgfa, &
1005 set_radius=set_radius_a, &
1010 IF (jkind .NE. jkind_old)
THEN
1012 CALL get_qs_kind(qs_kind_set(jkind), &
1013 basis_set=orb_basis_set, basis_type=basis_type)
1015 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1016 first_sgf=first_sgfb, &
1017 kind_radius=kind_radius_b, &
1024 set_radius=set_radius_b, &
1030 CALL dbcsr_get_block_p(matrix=deltajp_a(1)%matrix, row=brow, col=bcol, &
1031 block=jp_block_a, found=den_found)
1032 CALL dbcsr_get_block_p(matrix=deltajp_b(1)%matrix, row=brow, col=bcol, &
1033 block=jp_block_b, found=den_found)
1034 CALL dbcsr_get_block_p(matrix=deltajp_c(1)%matrix, row=brow, col=bcol, &
1035 block=jp_block_c, found=den_found)
1036 IF (do_igaim)
CALL dbcsr_get_block_p(matrix=deltajp_d(1)%matrix, row=brow, col=bcol, &
1037 block=jp_block_d, found=den_found)
1039 IF (.NOT.
ASSOCIATED(jp_block_a)) &
1040 cpabort(
"p_block not associated in deltap")
1047 atom_pair_changed = .true.
1051 atom_pair_changed = .false.
1055 IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset)
THEN
1057 ncoa = npgfa(iset)*
ncoset(la_max(iset))
1058 sgfa = first_sgfa(1, iset)
1059 ncob = npgfb(jset)*
ncoset(lb_max(jset))
1060 sgfb = first_sgfb(1, jset)
1063 CALL dgemm(
"N",
"N", ncoa, nsgfb(jset), nsgfa(iset), &
1064 1.0_dp, sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
1065 jp_block_a(sgfa, sgfb),
SIZE(jp_block_a, 1), &
1066 0.0_dp, work(1, 1), maxco)
1067 CALL dgemm(
"N",
"T", ncoa, ncob, nsgfb(jset), &
1068 1.0_dp, work(1, 1), maxco, &
1069 sphi_b(1, sgfb),
SIZE(sphi_b, 1), &
1070 0.0_dp, jpab_a(1, 1), maxco)
1072 CALL dgemm(
"N",
"N", ncoa, nsgfb(jset), nsgfa(iset), &
1073 1.0_dp, sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
1074 jp_block_b(sgfa, sgfb),
SIZE(jp_block_b, 1), &
1075 0.0_dp, work(1, 1), maxco)
1076 CALL dgemm(
"N",
"T", ncoa, ncob, nsgfb(jset), &
1077 1.0_dp, work(1, 1), maxco, &
1078 sphi_b(1, sgfb),
SIZE(sphi_b, 1), &
1079 0.0_dp, jpab_b(1, 1), maxco)
1081 CALL dgemm(
"N",
"N", ncoa, nsgfb(jset), nsgfa(iset), &
1082 1.0_dp, sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
1083 jp_block_c(sgfa, sgfb),
SIZE(jp_block_c, 1), &
1084 0.0_dp, work(1, 1), maxco)
1085 CALL dgemm(
"N",
"T", ncoa, ncob, nsgfb(jset), &
1086 1.0_dp, work(1, 1), maxco, &
1087 sphi_b(1, sgfb),
SIZE(sphi_b, 1), &
1088 0.0_dp, jpab_c(1, 1), maxco)
1091 CALL dgemm(
"N",
"N", ncoa, nsgfb(jset), nsgfa(iset), &
1092 1.0_dp, sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
1093 jp_block_d(sgfa, sgfb),
SIZE(jp_block_d, 1), &
1094 0.0_dp, work(1, 1), maxco)
1095 CALL dgemm(
"N",
"T", ncoa, ncob, nsgfb(jset), &
1096 1.0_dp, work(1, 1), maxco, &
1097 sphi_b(1, sgfb),
SIZE(sphi_b, 1), &
1098 0.0_dp, jpab_d(1, 1), maxco)
1108 adbmdab_func = grid_func_adbmdab_x
1110 adbmdab_func = grid_func_adbmdab_y
1112 adbmdab_func = grid_func_adbmdab_z
1114 cpabort(
"invalid idir")
1117 rab(:) = tasks(itask)%rab
1118 rb(:) = ra(:) + rab(:)
1119 zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
1120 f = zetb(jpgf, jset)/zetp
1121 prefactor = exp(-zeta(ipgf, iset)*f*dot_product(rab, rab))
1122 rp(:) = ra(:) + f*rab(:)
1124 na1 = (ipgf - 1)*
ncoset(la_max(iset)) + 1
1125 na2 = ipgf*
ncoset(la_max(iset))
1126 nb1 = (jpgf - 1)*
ncoset(lb_max(jset)) + 1
1127 nb2 = jpgf*
ncoset(lb_max(jset))
1136 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1137 lb_min=lb_min(jset), lb_max=lb_max(jset), &
1138 ra=ra, rb=rb, rp=rp, zetp=zetp, eps=eps_rho_rspace, &
1139 prefactor=prefactor, cutoff=1.0_dp)
1141 CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1142 la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1143 ra, rab, scale, jpab_a, na1 - 1, nb1 - 1, &
1144 rs_current(igrid_level), &
1145 radius=radius, ga_gb_function=adbmdab_func)
1149 IF (scale2 .NE. 0.0_dp)
THEN
1150 CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1151 la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1152 ra, rab, scale2, jpab_d, na1 - 1, nb1 - 1, &
1153 rs_rho(igrid_level), &
1154 radius=radius, ga_gb_function=grid_func_ab)
1160 current_env%rs_buf(igrid_level)%r(:, :, :) = 0.0_dp
1161 CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1162 la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1163 ra, rab, scale, jpab_b, na1 - 1, nb1 - 1, &
1165 ga_gb_function=adbmdab_func, &
1166 rsgrid=current_env%rs_buf(igrid_level))
1167 CALL collocate_gauge_ortho(rsgrid=current_env%rs_buf(igrid_level), &
1168 rsbuf=rs_current(igrid_level), &
1169 rsgauge=rs_gauge(iiib, igrid_level), &
1170 cube_info=cube_info(igrid_level), radius=radius, &
1171 zeta=zeta(ipgf, iset), zetb=zetb(jpgf, jset), &
1172 ra=ra, rab=rab, ir=iiib)
1178 current_env%rs_buf(igrid_level)%r(:, :, :) = 0.0_dp
1179 CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1180 la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1181 ra, rab, scale, jpab_c, na1 - 1, nb1 - 1, &
1183 ga_gb_function=adbmdab_func, &
1184 rsgrid=current_env%rs_buf(igrid_level))
1185 CALL collocate_gauge_ortho(rsgrid=current_env%rs_buf(igrid_level), &
1186 rsbuf=rs_current(igrid_level), &
1187 rsgauge=rs_gauge(iib, igrid_level), &
1188 cube_info=cube_info(igrid_level), radius=radius, &
1189 zeta=zeta(ipgf, iset), zetb=zetb(jpgf, jset), &
1190 ra=ra, rab=rab, ir=iib)
1196 CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1197 la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1198 ra, rab, scale, jpab_b, na1 - 1, nb1 - 1, &
1199 rs_current(igrid_level), &
1201 ga_gb_function=encode_ardbmdarb_func(idir=idir, ir=iiib))
1206 CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1207 la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1208 ra, rab, scale, jpab_c, na1 - 1, nb1 - 1, &
1209 rs_current(igrid_level), &
1211 ga_gb_function=encode_ardbmdarb_func(idir=idir, ir=iib))
1218 DO igrid_level = 1, gridlevel_info%ngrid_levels
1219 CALL rs_grid_mult_and_add(rs_current(igrid_level), rs_rho(igrid_level), &
1220 rs_gauge(idir2, igrid_level), 1.0_dp)
1225 IF (distributed_rs_grids)
THEN
1226 CALL dbcsr_deallocate_matrix(deltajp_a(1)%matrix)
1227 CALL dbcsr_deallocate_matrix(deltajp_b(1)%matrix)
1228 CALL dbcsr_deallocate_matrix(deltajp_c(1)%matrix)
1229 IF (do_igaim)
CALL dbcsr_deallocate_matrix(deltajp_d(1)%matrix)
1231 DEALLOCATE (deltajp_a, deltajp_b, deltajp_c, deltajp_d)
1233 DEALLOCATE (jpabt_a, jpabt_b, jpabt_c, jpabt_d, workt, tasks)
1235 IF (
ASSOCIATED(atom_pair_send))
DEALLOCATE (atom_pair_send)
1236 IF (
ASSOCIATED(atom_pair_recv))
DEALLOCATE (atom_pair_recv)
1238 CALL density_rs2pw(pw_env, rs_current, current_rs, current_gs)
1239 DO i = 1,
SIZE(rs_current)
1240 CALL rs_grid_release(rs_current(i))
1243 DO i = 1,
SIZE(rs_rho)
1244 IF (rs_descs(i)%rs_desc%distributed .AND. .NOT. my_retain_rsgrid)
THEN
1245 CALL rs_grid_release(rs_rho(i))
1250 DEALLOCATE (rs_current)
1252 CALL timestop(handle)
1262 FUNCTION encode_ardbmdarb_func(idir, ir)
RESULT(func)
1263 INTEGER,
INTENT(IN) :: idir, ir
1266 cpassert(1 <= idir .AND. idir <= 3 .AND. 1 <= ir .AND. ir <= 3)
1267 SELECT CASE (10*idir + ir)
1269 func = grid_func_ardbmdarb_xx
1271 func = grid_func_ardbmdarb_xy
1273 func = grid_func_ardbmdarb_xz
1275 func = grid_func_ardbmdarb_yx
1277 func = grid_func_ardbmdarb_yy
1279 func = grid_func_ardbmdarb_yz
1281 func = grid_func_ardbmdarb_zx
1283 func = grid_func_ardbmdarb_zy
1285 func = grid_func_ardbmdarb_zz
1287 cpabort(
"invalid idir or iiiB")
1289 END FUNCTION encode_ardbmdarb_func
1304 SUBROUTINE collocate_gauge_ortho(rsgrid, rsbuf, rsgauge, cube_info, radius, ra, rab, zeta, zetb, ir)
1305 TYPE(realspace_grid_type) :: rsgrid, rsbuf, rsgauge
1306 TYPE(cube_info_type),
INTENT(IN) :: cube_info
1307 REAL(kind=dp),
INTENT(IN) :: radius
1308 REAL(kind=dp),
DIMENSION(3),
INTENT(IN) :: ra, rab
1309 REAL(kind=dp),
INTENT(IN) :: zeta, zetb
1310 INTEGER,
INTENT(IN) :: ir
1312 INTEGER :: cmax, i, ig, igmax, igmin, j, j2, jg, &
1313 jg2, jgmin, k, k2, kg, kg2, kgmin, &
1314 length, offset, sci, start
1315 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: map
1316 INTEGER,
DIMENSION(3) :: cubecenter, lb_cube, ng, ub_cube
1317 INTEGER,
DIMENSION(:),
POINTER :: sphere_bounds
1318 REAL(kind=dp) :: f, point(3, 4), res(4), x, y, y2, z, z2, &
1320 REAL(kind=dp),
DIMENSION(3) :: dr, rap, rb, rbp, roffset, rp
1321 REAL(kind=dp),
DIMENSION(:, :, :),
POINTER :: gauge, grid, grid_buf
1323 cpassert(rsgrid%desc%orthorhombic)
1324 NULLIFY (sphere_bounds)
1326 grid => rsgrid%r(:, :, :)
1327 grid_buf => rsbuf%r(:, :, :)
1328 gauge => rsgauge%r(:, :, :)
1334 rbp(:) = rap(:) - rab(:)
1335 rp(:) = ra(:) + rap(:)
1336 rb(:) = ra(:) + rab(:)
1339 ng(:) = rsgrid%desc%npts(:)
1340 dr(1) = rsgrid%desc%dh(1, 1)
1341 dr(2) = rsgrid%desc%dh(2, 2)
1342 dr(3) = rsgrid%desc%dh(3, 3)
1345 CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds)
1346 cmax = maxval(ub_cube)
1355 ALLOCATE (map(-cmax:cmax, 3))
1357 CALL compute_cube_center(cubecenter, rsgrid%desc, zeta, zetb, ra, rab)
1358 roffset(:) = rp(:) - real(cubecenter(:), dp)*dr(:)
1362 IF (rsgrid%desc%perd(i) == 1)
THEN
1365 offset =
modulo(cubecenter(i) + start, ng(i)) + 1 - start
1366 length = min(ub_cube(i), ng(i) - offset) - start
1367 DO ig = start, start + length
1368 map(ig, i) = ig + offset
1370 IF (start + length .GE. ub_cube(i))
EXIT
1371 start = start + length + 1
1375 offset =
modulo(cubecenter(i) + lb_cube(i) + rsgrid%desc%lb(i) - rsgrid%lb_local(i), ng(i)) + 1 - lb_cube(i)
1377 IF (ub_cube(i) + offset > ubound(grid, i) .OR. lb_cube(i) + offset < lbound(grid, i))
THEN
1380 DO ig = lb_cube(i), ub_cube(i)
1381 map(ig, i) = ig + offset
1388 kgmin = sphere_bounds(sci)
1394 jgmin = sphere_bounds(sci)
1396 z = (real(kg, dp) + real(cubecenter(3), dp))*dr(3)
1397 z2 = (real(kg2, dp) + real(cubecenter(3), dp))*dr(3)
1402 igmin = sphere_bounds(sci)
1405 y = (real(jg, dp) + real(cubecenter(2), dp))*dr(2)
1406 y2 = (real(jg2, dp) + real(cubecenter(2), dp))*dr(2)
1407 DO ig = igmin, igmax
1409 x = (real(ig, dp) + real(cubecenter(1), dp))*dr(1)
1410 point(1, 1) = x; point(2, 1) = y; point(3, 1) = z
1411 point(1, 2) = x; point(2, 2) = y2; point(3, 2) = z
1412 point(1, 3) = x; point(2, 3) = y; point(3, 3) = z2
1413 point(1, 4) = x; point(2, 4) = y2; point(3, 4) = z2
1415 res(1) = (point(ir, 1) - rb(ir)) - gauge(i, j, k)
1416 res(2) = (point(ir, 2) - rb(ir)) - gauge(i, j2, k)
1417 res(3) = (point(ir, 3) - rb(ir)) - gauge(i, j, k2)
1418 res(4) = (point(ir, 4) - rb(ir)) - gauge(i, j2, k2)
1420 grid_buf(i, j, k) = grid_buf(i, j, k) + grid(i, j, k)*res(1)
1421 grid_buf(i, j2, k) = grid_buf(i, j2, k) + grid(i, j2, k)*res(2)
1422 grid_buf(i, j, k2) = grid_buf(i, j, k2) + grid(i, j, k2)*res(3)
1423 grid_buf(i, j2, k2) = grid_buf(i, j2, k2) + grid(i, j2, k2)*res(4)
1427 END SUBROUTINE collocate_gauge_ortho
1434 SUBROUTINE current_set_gauge(current_env, qs_env)
1436 TYPE(current_env_type) :: current_env
1437 TYPE(qs_environment_type),
POINTER :: qs_env
1439 CHARACTER(LEN=*),
PARAMETER :: routinen =
'current_set_gauge'
1443 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :) :: box_data
1444 INTEGER :: handle, igrid_level, nbox(3), gauge
1445 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: box_ptr
1446 TYPE(realspace_grid_desc_p_type),
DIMENSION(:), &
1448 TYPE(pw_env_type),
POINTER :: pw_env
1449 TYPE(realspace_grid_type),
DIMENSION(:, :),
POINTER :: rs_gauge
1451 TYPE(box_type),
DIMENSION(:, :, :),
POINTER :: box
1452 LOGICAL :: use_old_gauge_atom
1454 NULLIFY (rs_gauge, box)
1456 CALL timeset(routinen, handle)
1458 CALL get_current_env(current_env=current_env, &
1459 use_old_gauge_atom=use_old_gauge_atom, &
1460 rs_gauge=rs_gauge, &
1463 IF (gauge .EQ. current_gauge_atom)
THEN
1464 CALL get_qs_env(qs_env=qs_env, &
1466 CALL pw_env_get(pw_env=pw_env, &
1470 IF (use_old_gauge_atom)
THEN
1471 CALL box_atoms(qs_env)
1473 CALL box_atoms_new(current_env, qs_env, box)
1477 DO igrid_level = pw_env%gridlevel_info%ngrid_levels, 1, -1
1480 CALL rs_grid_create(rs_gauge(idir, igrid_level), rs_descs(igrid_level)%rs_desc)
1483 IF (use_old_gauge_atom)
THEN
1484 CALL collocate_gauge(current_env, qs_env, &
1485 rs_gauge(1, igrid_level), &
1486 rs_gauge(2, igrid_level), &
1487 rs_gauge(3, igrid_level))
1489 CALL collocate_gauge_new(current_env, qs_env, &
1490 rs_gauge(1, igrid_level), &
1491 rs_gauge(2, igrid_level), &
1492 rs_gauge(3, igrid_level), &
1498 ALLOCATE (current_env%rs_buf(pw_env%gridlevel_info%ngrid_levels))
1499 DO igrid_level = 1, pw_env%gridlevel_info%ngrid_levels
1500 CALL rs_grid_create(current_env%rs_buf(igrid_level), rs_descs(igrid_level)%rs_desc)
1503 DEALLOCATE (box_ptr, box_data)
1504 CALL deallocate_box(box)
1507 CALL timestop(handle)
1515 SUBROUTINE box_atoms(qs_env)
1516 TYPE(qs_environment_type),
POINTER :: qs_env
1518 REAL(kind=dp),
PARAMETER :: box_size_guess = 5.0_dp
1520 INTEGER :: i, iatom, ibox, ii, jbox, kbox, natms
1521 REAL(dp) :: offset(3)
1522 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :) :: ratom
1523 TYPE(cell_type),
POINTER :: cell
1524 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
1525 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1527 CALL get_qs_env(qs_env=qs_env, &
1528 qs_kind_set=qs_kind_set, &
1530 particle_set=particle_set)
1532 natms =
SIZE(particle_set, 1)
1533 ALLOCATE (ratom(3, natms))
1536 nbox(1) = ceiling(cell%hmat(1, 1)/box_size_guess)
1537 nbox(2) = ceiling(cell%hmat(2, 2)/box_size_guess)
1538 nbox(3) = ceiling(cell%hmat(3, 3)/box_size_guess)
1540 dbox(1) = cell%hmat(1, 1)/real(nbox(1), dp)
1541 dbox(2) = cell%hmat(2, 2)/real(nbox(2), dp)
1542 dbox(3) = cell%hmat(3, 3)/real(nbox(3), dp)
1544 ALLOCATE (box_ptr(0:nbox(1), 0:nbox(2) - 1, 0:nbox(3) - 1), box_data(3, natms))
1545 box_data(:, :) = huge(0.0_dp)
1546 box_ptr(:, :, :) = huge(0)
1548 offset(1) = cell%hmat(1, 1)*0.5_dp
1549 offset(2) = cell%hmat(2, 2)*0.5_dp
1550 offset(3) = cell%hmat(3, 3)*0.5_dp
1552 ratom(:, iatom) =
pbc(particle_set(iatom)%r(:), cell) + offset(:)
1556 DO kbox = 0, nbox(3) - 1
1557 DO jbox = 0, nbox(2) - 1
1558 box_ptr(0, jbox, kbox) = i
1559 DO ibox = 0, nbox(1) - 1
1562 IF (int(ratom(1, iatom)/dbox(1)) .EQ. ibox .AND. &
1563 int(ratom(2, iatom)/dbox(2)) .EQ. jbox .AND. &
1564 int(ratom(3, iatom)/dbox(3)) .EQ. kbox)
THEN
1565 box_data(:, i) = ratom(:, iatom) - offset(:)
1570 box_ptr(ibox + 1, jbox, kbox) = box_ptr(ibox, jbox, kbox) + ii
1576 DO kbox = 0, nbox(3) - 1
1577 DO jbox = 0, nbox(2) - 1
1578 DO ibox = 0, nbox(1) - 1
1579 WRITE (*, *)
'box=', ibox, jbox, kbox
1580 WRITE (*, *)
'nbr atom=', box_ptr(ibox + 1, jbox, kbox) - box_ptr(ibox, jbox, kbox)
1581 DO iatom = box_ptr(ibox, jbox, kbox), box_ptr(ibox + 1, jbox, kbox) - 1
1582 WRITE (*, *)
'iatom=', iatom
1583 WRITE (*,
'(A,3E14.6)')
'coor=', box_data(:, iatom)
1590 END SUBROUTINE box_atoms
1600 SUBROUTINE collocate_gauge(current_env, qs_env, rs_grid_x, rs_grid_y, rs_grid_z)
1602 TYPE(current_env_type) :: current_env
1603 TYPE(qs_environment_type),
POINTER :: qs_env
1604 TYPE(realspace_grid_type),
INTENT(IN) :: rs_grid_x, rs_grid_y, rs_grid_z
1606 INTEGER :: i, iatom, ibeg, ibox, iend,
imax,
imin, &
1607 j, jatom, jbox, jmax, jmin, k, kbox, &
1608 kmax, kmin, lb(3), lb_local(3), natms, &
1610 REAL(kind=dp) :: ab, buf_tmp, dist, dr(3), &
1611 gauge_atom_radius, offset(3), pa, pb, &
1612 point(3), pra(3), r(3), res(3), summe, &
1614 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: buf, nrm_atms_pnt
1615 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: atms_pnt, ratom
1616 REAL(kind=dp),
DIMENSION(:, :, :),
POINTER :: grid_x, grid_y, grid_z
1617 TYPE(cell_type),
POINTER :: cell
1618 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
1619 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1623 CALL get_current_env(current_env=current_env, &
1624 gauge_atom_radius=gauge_atom_radius)
1626 CALL get_qs_env(qs_env=qs_env, &
1627 qs_kind_set=qs_kind_set, &
1629 particle_set=particle_set)
1631 natms =
SIZE(particle_set, 1)
1632 dr(1) = rs_grid_x%desc%dh(1, 1)
1633 dr(2) = rs_grid_x%desc%dh(2, 2)
1634 dr(3) = rs_grid_x%desc%dh(3, 3)
1635 lb(:) = rs_grid_x%desc%lb(:)
1636 lb_local(:) = rs_grid_x%lb_local(:)
1637 grid_x => rs_grid_x%r(:, :, :)
1638 grid_y => rs_grid_y%r(:, :, :)
1639 grid_z => rs_grid_z%r(:, :, :)
1640 ng(:) = ubound(grid_x)
1641 offset(1) = cell%hmat(1, 1)*0.5_dp
1642 offset(2) = cell%hmat(2, 2)*0.5_dp
1643 offset(3) = cell%hmat(3, 3)*0.5_dp
1644 ALLOCATE (buf(natms), ratom(3, natms), atms_pnt(3, natms), nrm_atms_pnt(natms))
1651 point(3) = real(k - 1 + lb_local(3) - lb(3), dp)*dr(3)
1652 point(2) = real(j - 1 + lb_local(2) - lb(2), dp)*dr(2)
1653 point(1) = real(i - 1 + lb_local(1) - lb(1), dp)*dr(1)
1654 point =
pbc(point, cell)
1658 kmin = int((point(3) + offset(3) - gauge_atom_radius)/dbox(3))
1659 kmax = int((point(3) + offset(3) + gauge_atom_radius)/dbox(3))
1660 IF (kmax - kmin + 1 .GT. nbox(3))
THEN
1664 DO kbox = kmin, kmax
1665 jmin = int((point(2) + offset(2) - gauge_atom_radius)/dbox(2))
1666 jmax = int((point(2) + offset(2) + gauge_atom_radius)/dbox(2))
1667 IF (jmax - jmin + 1 .GT. nbox(2))
THEN
1671 DO jbox = jmin, jmax
1672 imin = int((point(1) + offset(1) - gauge_atom_radius)/dbox(1))
1673 imax = int((point(1) + offset(1) + gauge_atom_radius)/dbox(1))
1674 IF (
imax -
imin + 1 .GT. nbox(1))
THEN
1680 iend = box_ptr(
modulo(ibox, nbox(1)) + 1,
modulo(jbox, nbox(2)),
modulo(kbox, nbox(3))) - 1
1681 DO iatom = ibeg, iend
1682 r(:) =
pbc(box_data(:, iatom) - point(:), cell) + point(:)
1683 dist = (r(1) - point(1))**2 + (r(2) - point(2))**2 + (r(3) - point(3))**2
1684 IF (dist .LT. gauge_atom_radius**2)
THEN
1685 natms_local = natms_local + 1
1686 ratom(:, natms_local) = r(:)
1692 atms_pnt(1, natms_local) = x
1693 atms_pnt(2, natms_local) = y
1694 atms_pnt(3, natms_local) = z
1695 nrm_atms_pnt(natms_local) = sqrt(x*x + y*y + z*z)
1702 IF (natms_local .GT. 0)
THEN
1705 DO iatom = 1, natms_local
1707 pra(1) = atms_pnt(1, iatom)
1708 pra(2) = atms_pnt(2, iatom)
1709 pra(3) = atms_pnt(3, iatom)
1710 pa = nrm_atms_pnt(iatom)
1711 DO jatom = 1, natms_local
1712 IF (iatom .EQ. jatom) cycle
1713 pb = nrm_atms_pnt(jatom)
1714 x = pra(1) - atms_pnt(1, jatom)
1715 y = pra(2) - atms_pnt(2, jatom)
1716 z = pra(3) - atms_pnt(3, jatom)
1717 ab = sqrt(x*x + y*y + z*z)
1720 tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
1721 tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
1722 tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
1723 buf_tmp = buf_tmp*0.5_dp*(1.0_dp - tmp)
1725 buf(iatom) = buf_tmp
1731 DO iatom = 1, natms_local
1732 res(1) = res(1) + ratom(1, iatom)*buf(iatom)
1733 res(2) = res(2) + ratom(2, iatom)*buf(iatom)
1734 res(3) = res(3) + ratom(3, iatom)*buf(iatom)
1735 summe = summe + buf(iatom)
1737 res(1) = res(1)/summe
1738 res(2) = res(2)/summe
1739 res(3) = res(3)/summe
1740 grid_x(i, j, k) = point(1) - res(1)
1741 grid_y(i, j, k) = point(2) - res(2)
1742 grid_z(i, j, k) = point(3) - res(3)
1744 grid_x(i, j, k) = 0.0_dp
1745 grid_y(i, j, k) = 0.0_dp
1746 grid_z(i, j, k) = 0.0_dp
1752 DEALLOCATE (buf, ratom, atms_pnt, nrm_atms_pnt)
1754 END SUBROUTINE collocate_gauge
1762 SUBROUTINE box_atoms_new(current_env, qs_env, box)
1763 TYPE(current_env_type) :: current_env
1764 TYPE(qs_environment_type),
POINTER :: qs_env
1765 TYPE(box_type),
DIMENSION(:, :, :),
POINTER :: box
1767 CHARACTER(LEN=*),
PARAMETER :: routinen =
'box_atoms_new'
1769 INTEGER :: handle, i, iatom, ibeg, ibox, iend, &
1770 ifind, ii,
imax,
imin, j, jatom, jbox, &
1771 jmax, jmin, k, kbox, kmax, kmin, &
1773 REAL(dp) :: gauge_atom_radius, offset(3), scale
1774 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :) :: ratom
1775 REAL(dp),
DIMENSION(:, :),
POINTER :: r_ptr
1776 REAL(kind=dp) :: box_center(3), box_center_wrap(3), &
1777 box_size_guess, r(3)
1778 TYPE(cell_type),
POINTER :: cell
1779 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
1780 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1782 CALL timeset(routinen, handle)
1784 CALL get_qs_env(qs_env=qs_env, &
1785 qs_kind_set=qs_kind_set, &
1787 particle_set=particle_set)
1789 CALL get_current_env(current_env=current_env, &
1790 gauge_atom_radius=gauge_atom_radius)
1794 box_size_guess = gauge_atom_radius/scale
1796 natms =
SIZE(particle_set, 1)
1797 ALLOCATE (ratom(3, natms))
1801 nbox(1) = ceiling(cell%hmat(1, 1)/box_size_guess)
1802 nbox(2) = ceiling(cell%hmat(2, 2)/box_size_guess)
1803 nbox(3) = ceiling(cell%hmat(3, 3)/box_size_guess)
1804 dbox(1) = cell%hmat(1, 1)/real(nbox(1), dp)
1805 dbox(2) = cell%hmat(2, 2)/real(nbox(2), dp)
1806 dbox(3) = cell%hmat(3, 3)/real(nbox(3), dp)
1807 ALLOCATE (box_ptr(0:nbox(1), 0:nbox(2) - 1, 0:nbox(3) - 1), box_data(3, natms))
1808 box_data(:, :) = huge(0.0_dp)
1809 box_ptr(:, :, :) = huge(0)
1811 offset(1) = cell%hmat(1, 1)*0.5_dp
1812 offset(2) = cell%hmat(2, 2)*0.5_dp
1813 offset(3) = cell%hmat(3, 3)*0.5_dp
1815 ratom(:, iatom) =
pbc(particle_set(iatom)%r(:), cell)
1819 DO kbox = 0, nbox(3) - 1
1820 DO jbox = 0, nbox(2) - 1
1821 box_ptr(0, jbox, kbox) = i
1822 DO ibox = 0, nbox(1) - 1
1825 IF (
modulo(floor(ratom(1, iatom)/dbox(1)), nbox(1)) .EQ. ibox .AND. &
1826 modulo(floor(ratom(2, iatom)/dbox(2)), nbox(2)) .EQ. jbox .AND. &
1827 modulo(floor(ratom(3, iatom)/dbox(3)), nbox(3)) .EQ. kbox)
THEN
1828 box_data(:, i) = ratom(:, iatom)
1833 box_ptr(ibox + 1, jbox, kbox) = box_ptr(ibox, jbox, kbox) + ii
1839 DO kbox = 0, nbox(3) - 1
1840 DO jbox = 0, nbox(2) - 1
1841 DO ibox = 0, nbox(1) - 1
1842 IF (box_ptr(ibox + 1, jbox, kbox) - box_ptr(ibox, jbox, kbox) .GT. 0)
THEN
1843 WRITE (*, *)
'box=', ibox, jbox, kbox
1844 WRITE (*, *)
'nbr atom=', box_ptr(ibox + 1, jbox, kbox) - box_ptr(ibox, jbox, kbox)
1845 DO iatom = box_ptr(ibox, jbox, kbox), box_ptr(ibox + 1, jbox, kbox) - 1
1846 WRITE (*,
'(A,I3,3E14.6)')
'coor=', iatom, box_data(:, iatom)
1855 ALLOCATE (box(0:nbox(1) - 1, 0:nbox(2) - 1, 0:nbox(3) - 1))
1858 DO k = 0, nbox(3) - 1
1859 DO j = 0, nbox(2) - 1
1860 DO i = 0, nbox(1) - 1
1862 box_center(1) = (real(i, dp) + 0.5_dp)*dbox(1)
1863 box_center(2) = (real(j, dp) + 0.5_dp)*dbox(2)
1864 box_center(3) = (real(k, dp) + 0.5_dp)*dbox(3)
1865 box_center_wrap =
pbc(box_center, cell)
1869 kmin = floor((box_center(3) - gauge_atom_radius)/dbox(3))
1870 kmax = floor((box_center(3) + gauge_atom_radius)/dbox(3))
1871 IF (kmax - kmin + 1 .GT. nbox(3))
THEN
1875 DO kbox = kmin, kmax
1876 jmin = floor((box_center(2) - gauge_atom_radius)/dbox(2))
1877 jmax = floor((box_center(2) + gauge_atom_radius)/dbox(2))
1878 IF (jmax - jmin + 1 .GT. nbox(2))
THEN
1882 DO jbox = jmin, jmax
1883 imin = floor((box_center(1) - gauge_atom_radius)/dbox(1))
1884 imax = floor((box_center(1) + gauge_atom_radius)/dbox(1))
1885 IF (
imax -
imin + 1 .GT. nbox(1))
THEN
1891 iend = box_ptr(
modulo(ibox, nbox(1)) + 1,
modulo(jbox, nbox(2)),
modulo(kbox, nbox(3))) - 1
1892 DO iatom = ibeg, iend
1893 r =
pbc(box_center_wrap(:) - box_data(:, iatom), cell)
1894 IF (abs(r(1)) .LE. (scale + 0.5_dp)*dbox(1) .AND. &
1895 abs(r(2)) .LE. (scale + 0.5_dp)*dbox(2) .AND. &
1896 abs(r(3)) .LE. (scale + 0.5_dp)*dbox(3))
THEN
1897 natms_local = natms_local + 1
1898 ratom(:, natms_local) = box_data(:, iatom)
1906 box(i, j, k)%n = natms_local
1907 NULLIFY (box(i, j, k)%r)
1908 IF (natms_local .GT. 0)
THEN
1909 ALLOCATE (box(i, j, k)%r(3, natms_local))
1910 r_ptr => box(i, j, k)%r
1911 CALL dcopy(3*natms_local, ratom(1, 1), 1, r_ptr(1, 1), 1)
1918 DO k = 0, nbox(3) - 1
1919 DO j = 0, nbox(2) - 1
1920 DO i = 0, nbox(1) - 1
1921 IF (box(i, j, k)%n .GT. 0)
THEN
1923 WRITE (*, *)
'box=', i, j, k
1924 box_center(1) = (real(i, dp) + 0.5_dp)*dbox(1)
1925 box_center(2) = (real(j, dp) + 0.5_dp)*dbox(2)
1926 box_center(3) = (real(k, dp) + 0.5_dp)*dbox(3)
1927 box_center =
pbc(box_center, cell)
1928 WRITE (*,
'(A,3E14.6)')
'box_center=', box_center
1929 WRITE (*, *)
'nbr atom=', box(i, j, k)%n
1930 r_ptr => box(i, j, k)%r
1931 DO iatom = 1, box(i, j, k)%n
1932 WRITE (*,
'(A,I3,3E14.6)')
'coor=', iatom, r_ptr(:, iatom)
1933 r(:) =
pbc(box_center(:) - r_ptr(:, iatom), cell)
1934 IF (abs(r(1)) .GT. (scale + 0.5_dp)*dbox(1) .OR. &
1935 abs(r(2)) .GT. (scale + 0.5_dp)*dbox(2) .OR. &
1936 abs(r(3)) .GT. (scale + 0.5_dp)*dbox(3))
THEN
1937 WRITE (*, *)
'error too many atoms'
1938 WRITE (*, *)
'dist=', abs(r(:))
1939 WRITE (*, *)
'large_dist=', (scale + 0.5_dp)*dbox
1950 DO k = 0, nbox(3) - 1
1951 DO j = 0, nbox(2) - 1
1952 DO i = 0, nbox(1) - 1
1953 box_center(1) = (real(i, dp) + 0.5_dp)*dbox(1)
1954 box_center(2) = (real(j, dp) + 0.5_dp)*dbox(2)
1955 box_center(3) = (real(k, dp) + 0.5_dp)*dbox(3)
1956 box_center =
pbc(box_center, cell)
1957 r_ptr => box(i, j, k)%r
1959 r(:) =
pbc(box_center(:) - ratom(:, iatom), cell)
1961 DO jatom = 1, box(i, j, k)%n
1962 IF (sum(abs(ratom(:, iatom) - r_ptr(:, jatom))) .LT. 1e-10_dp) ifind = 1
1965 IF (ifind .EQ. 0)
THEN
1967 IF (dot_product(r, r) .LT. (gauge_atom_radius*gauge_atom_radius))
THEN
1968 WRITE (*, *)
'error atom too close'
1969 WRITE (*, *)
'iatom', iatom
1970 WRITE (*, *)
'box_center', box_center
1971 WRITE (*, *)
'ratom', ratom(:, iatom)
1972 WRITE (*, *)
'gauge_atom_radius', gauge_atom_radius
1984 CALL timestop(handle)
1986 END SUBROUTINE box_atoms_new
1997 SUBROUTINE collocate_gauge_new(current_env, qs_env, rs_grid_x, rs_grid_y, rs_grid_z, box)
1999 TYPE(current_env_type) :: current_env
2000 TYPE(qs_environment_type),
POINTER :: qs_env
2001 TYPE(realspace_grid_type),
INTENT(IN) :: rs_grid_x, rs_grid_y, rs_grid_z
2002 TYPE(box_type),
DIMENSION(:, :, :),
POINTER :: box
2004 CHARACTER(LEN=*),
PARAMETER :: routinen =
'collocate_gauge_new'
2006 INTEGER :: delta_lb(3), handle, i, iatom, ib, ibe, ibox, ibs, ie, is, j, jatom, jb, jbe, &
2007 jbox, jbs, je, js, k, kb, kbe, kbox, kbs, ke, ks, lb(3), lb_local(3), natms, &
2008 natms_local0, natms_local1, ng(3)
2009 REAL(dp),
DIMENSION(:, :),
POINTER :: r_ptr
2010 REAL(kind=dp) :: ab, box_center(3), buf_tmp, dist, dr(3), &
2011 gauge_atom_radius, offset(3), pa, pb, &
2012 point(3), pra(3), r(3), res(3), summe, &
2014 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: buf, nrm_atms_pnt
2015 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: atms_pnt, ratom
2016 REAL(kind=dp),
DIMENSION(:, :, :),
POINTER :: grid_x, grid_y, grid_z
2017 TYPE(cell_type),
POINTER :: cell
2018 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
2019 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2021 CALL timeset(routinen, handle)
2024 CALL get_current_env(current_env=current_env, &
2025 gauge_atom_radius=gauge_atom_radius)
2027 CALL get_qs_env(qs_env=qs_env, &
2028 qs_kind_set=qs_kind_set, &
2030 particle_set=particle_set)
2032 natms =
SIZE(particle_set, 1)
2033 dr(1) = rs_grid_x%desc%dh(1, 1)
2034 dr(2) = rs_grid_x%desc%dh(2, 2)
2035 dr(3) = rs_grid_x%desc%dh(3, 3)
2036 lb(:) = rs_grid_x%desc%lb(:)
2037 lb_local(:) = rs_grid_x%lb_local(:)
2038 grid_x => rs_grid_x%r(:, :, :)
2039 grid_y => rs_grid_y%r(:, :, :)
2040 grid_z => rs_grid_z%r(:, :, :)
2041 ng(:) = ubound(grid_x)
2042 delta_lb(:) = lb_local(:) - lb(:)
2043 offset(1) = cell%hmat(1, 1)*0.5_dp
2044 offset(2) = cell%hmat(2, 2)*0.5_dp
2045 offset(3) = cell%hmat(3, 3)*0.5_dp
2046 ALLOCATE (buf(natms), ratom(3, natms), atms_pnt(3, natms), nrm_atms_pnt(natms))
2049 ibs = floor(real(delta_lb(1), dp)*dr(1)/dbox(1))
2050 ibe = floor(real(ng(1) - 1 + delta_lb(1), dp)*dr(1)/dbox(1))
2051 jbs = floor(real(delta_lb(2), dp)*dr(2)/dbox(2))
2052 jbe = floor(real(ng(2) - 1 + delta_lb(2), dp)*dr(2)/dbox(2))
2053 kbs = floor(real(delta_lb(3), dp)*dr(3)/dbox(3))
2054 kbe = floor(real(ng(3) - 1 + delta_lb(3), dp)*dr(3)/dbox(3))
2060 ibox =
modulo(ib, nbox(1))
2061 jbox =
modulo(jb, nbox(2))
2062 kbox =
modulo(kb, nbox(3))
2064 is = max(ceiling(real(ib, dp)*dbox(1)/dr(1)), delta_lb(1)) - delta_lb(1) + 1
2065 ie = min(floor(real(ib + 1, dp)*dbox(1)/dr(1)), ng(1) - 1 + delta_lb(1)) - delta_lb(1) + 1
2066 js = max(ceiling(real(jb, dp)*dbox(2)/dr(2)), delta_lb(2)) - delta_lb(2) + 1
2067 je = min(floor(real(jb + 1, dp)*dbox(2)/dr(2)), ng(2) - 1 + delta_lb(2)) - delta_lb(2) + 1
2068 ks = max(ceiling(real(kb, dp)*dbox(3)/dr(3)), delta_lb(3)) - delta_lb(3) + 1
2069 ke = min(floor(real(kb + 1, dp)*dbox(3)/dr(3)), ng(3) - 1 + delta_lb(3)) - delta_lb(3) + 1
2073 IF (real(ks - 1 + delta_lb(3), dp)*dr(3) .LT. real(kb, dp)*dbox(3) .OR. &
2074 REAL(ke - 1 + delta_lb(3), dp)*dr(3) .GT.
REAL(kb + 1, dp)*dbox(3)) then
2075 WRITE (*, *)
'box_k', real(kb, dp)*dbox(3), real(kb + 1, dp)*dbox(3)
2076 WRITE (*, *)
'point_k', real(ks - 1 + delta_lb(3), dp)*dr(3), real(ke - 1 + delta_lb(3), dp)*dr(3)
2077 WRITE (*, *)
'ibox', ibox,
'jbox', jbox,
'kbox', kbox
2078 WRITE (*, *)
'is,ie', is, ie,
' js,je', js, je,
' ks,ke', ks, ke
2079 WRITE (*, *)
'ibs,ibe', ibs, ibe,
' jbs,jbe', jbs, jbe,
' kbs,kbe', kbs, kbe
2080 cpabort(
"we stop_k")
2082 IF (real(js - 1 + delta_lb(2), dp)*dr(2) .LT. real(jb, dp)*dbox(2) .OR. &
2083 REAL(je - 1 + delta_lb(2), dp)*dr(2) .GT.
REAL(jb + 1, dp)*dbox(2)) then
2084 WRITE (*, *)
'box_j', real(jb, dp)*dbox(2), real(jb + 1, dp)*dbox(2)
2085 WRITE (*, *)
'point_j', real(js - 1 + delta_lb(2), dp)*dr(2), real(je - 1 + delta_lb(2), dp)*dr(2)
2086 WRITE (*, *)
'is,ie', is, ie,
' js,je', js, je,
' ks,ke', ks, ke
2087 WRITE (*, *)
'ibs,ibe', ibs, ibe,
' jbs,jbe', jbs, jbe,
' kbs,kbe', kbs, kbe
2088 cpabort(
"we stop_j")
2090 IF (real(is - 1 + delta_lb(1), dp)*dr(1) .LT. real(ib, dp)*dbox(1) .OR. &
2091 REAL(ie - 1 + delta_lb(1), dp)*dr(1) .GT.
REAL(ib + 1, dp)*dbox(1)) then
2092 WRITE (*, *)
'box_i', real(ib, dp)*dbox(1), real(ib + 1, dp)*dbox(1)
2093 WRITE (*, *)
'point_i', real(is - 1 + delta_lb(1), dp)*dr(1), real(ie - 1 + delta_lb(1), dp)*dr(1)
2094 WRITE (*, *)
'is,ie', is, ie,
' js,je', js, je,
' ks,ke', ks, ke
2095 WRITE (*, *)
'ibs,ibe', ibs, ibe,
' jbs,jbe', jbs, jbe,
' kbs,kbe', kbs, kbe
2096 cpabort(
"we stop_i")
2101 box_center(1) = (real(ibox, dp) + 0.5_dp)*dbox(1)
2102 box_center(2) = (real(jbox, dp) + 0.5_dp)*dbox(2)
2103 box_center(3) = (real(kbox, dp) + 0.5_dp)*dbox(3)
2106 natms_local0 = box(ibox, jbox, kbox)%n
2107 r_ptr => box(ibox, jbox, kbox)%r
2110 IF (natms_local0 .GT. 0)
THEN
2116 point(1) = real(i - 1 + delta_lb(1), dp)*dr(1)
2117 point(2) = real(j - 1 + delta_lb(2), dp)*dr(2)
2118 point(3) = real(k - 1 + delta_lb(3), dp)*dr(3)
2119 point =
pbc(point, cell)
2123 DO iatom = 1, natms_local0
2124 r(:) =
pbc(r_ptr(:, iatom) - point(:), cell) + point(:)
2125 dist = (r(1) - point(1))**2 + (r(2) - point(2))**2 + (r(3) - point(3))**2
2126 IF (dist .LT. gauge_atom_radius**2)
THEN
2127 natms_local1 = natms_local1 + 1
2128 ratom(:, natms_local1) = r(:)
2134 atms_pnt(1, natms_local1) = x
2135 atms_pnt(2, natms_local1) = y
2136 atms_pnt(3, natms_local1) = z
2137 nrm_atms_pnt(natms_local1) = sqrt(x*x + y*y + z*z)
2142 IF (natms_local1 .GT. 0)
THEN
2145 DO iatom = 1, natms_local1
2147 pra(1) = atms_pnt(1, iatom)
2148 pra(2) = atms_pnt(2, iatom)
2149 pra(3) = atms_pnt(3, iatom)
2150 pa = nrm_atms_pnt(iatom)
2151 DO jatom = 1, natms_local1
2152 IF (iatom .EQ. jatom) cycle
2153 pb = nrm_atms_pnt(jatom)
2154 x = pra(1) - atms_pnt(1, jatom)
2155 y = pra(2) - atms_pnt(2, jatom)
2156 z = pra(3) - atms_pnt(3, jatom)
2157 ab = sqrt(x*x + y*y + z*z)
2160 tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
2161 tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
2162 tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
2163 buf_tmp = buf_tmp*0.5_dp*(1.0_dp - tmp)
2165 buf(iatom) = buf_tmp
2171 DO iatom = 1, natms_local1
2172 res(1) = res(1) + ratom(1, iatom)*buf(iatom)
2173 res(2) = res(2) + ratom(2, iatom)*buf(iatom)
2174 res(3) = res(3) + ratom(3, iatom)*buf(iatom)
2175 summe = summe + buf(iatom)
2177 res(1) = res(1)/summe
2178 res(2) = res(2)/summe
2179 res(3) = res(3)/summe
2180 grid_x(i, j, k) = point(1) - res(1)
2181 grid_y(i, j, k) = point(2) - res(2)
2182 grid_z(i, j, k) = point(3) - res(3)
2184 grid_x(i, j, k) = 0.0_dp
2185 grid_y(i, j, k) = 0.0_dp
2186 grid_z(i, j, k) = 0.0_dp
2198 grid_x(i, j, k) = 0.0_dp
2199 grid_y(i, j, k) = 0.0_dp
2200 grid_z(i, j, k) = 0.0_dp
2211 DEALLOCATE (buf, ratom, atms_pnt, nrm_atms_pnt)
2213 CALL timestop(handle)
2215 END SUBROUTINE collocate_gauge_new
2221 SUBROUTINE deallocate_box(box)
2222 TYPE(box_type),
DIMENSION(:, :, :),
POINTER :: box
2226 IF (
ASSOCIATED(box))
THEN
2227 DO k = lbound(box, 3), ubound(box, 3)
2228 DO j = lbound(box, 2), ubound(box, 2)
2229 DO i = lbound(box, 1), ubound(box, 1)
2230 IF (
ASSOCIATED(box(i, j, k)%r))
THEN
2231 DEALLOCATE (box(i, j, k)%r)
2238 END SUBROUTINE deallocate_box
2239 END SUBROUTINE current_set_gauge
2249 TYPE(current_env_type) :: current_env
2250 TYPE(qs_environment_type),
POINTER :: qs_env
2251 INTEGER,
INTENT(IN) :: ib
2253 IF (current_env%full)
THEN
2254 CALL current_build_chi_many_centers(current_env, qs_env, ib)
2255 ELSE IF (current_env%nbr_center(1) > 1)
THEN
2256 CALL current_build_chi_many_centers(current_env, qs_env, ib)
2258 CALL current_build_chi_one_center(current_env, qs_env, ib)
2269 SUBROUTINE current_build_chi_many_centers(current_env, qs_env, iB)
2271 TYPE(current_env_type) :: current_env
2272 TYPE(qs_environment_type),
POINTER :: qs_env
2273 INTEGER,
INTENT(IN) :: ib
2275 CHARACTER(LEN=*),
PARAMETER :: routinen =
'current_build_chi_many_centers'
2277 INTEGER :: handle, icenter, idir, idir2, ii, iib, iii, iiib, ispin, istate, j, jstate, &
2278 max_states, nao, natom, nbr_center(2), nmo, nspins, nstate_loc, nstates(2), output_unit
2279 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf
2280 INTEGER,
DIMENSION(:),
POINTER :: row_blk_sizes
2281 LOGICAL :: chi_pbc, gapw
2282 REAL(dp) :: chi(3), chi_tmp, contrib, contrib2, &
2283 dk(3), int_current(3), &
2284 int_current_tmp, maxocc
2285 TYPE(cell_type),
POINTER :: cell
2286 TYPE(cp_2d_i_p_type),
DIMENSION(:),
POINTER :: center_list
2287 TYPE(cp_2d_r_p_type),
DIMENSION(:),
POINTER :: centers_set
2288 TYPE(cp_fm_struct_type),
POINTER :: tmp_fm_struct
2289 TYPE(cp_fm_type) :: psi0, psi_d, psi_p1, psi_p2, psi_rxp
2290 TYPE(cp_fm_type),
DIMENSION(3) :: p_rxp, r_p1, r_p2
2291 TYPE(cp_fm_type),
DIMENSION(9, 3) :: rr_p1, rr_p2, rr_rxp
2292 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: psi0_order
2293 TYPE(cp_fm_type),
DIMENSION(:, :),
POINTER :: psi1_d, psi1_p, psi1_rxp
2294 TYPE(cp_fm_type),
POINTER :: mo_coeff
2295 TYPE(cp_logger_type),
POINTER :: logger
2296 TYPE(dbcsr_distribution_type),
POINTER :: dbcsr_dist
2297 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: op_mom_ao, op_p_ao
2298 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: op_mom_der_ao
2299 TYPE(dft_control_type),
POINTER :: dft_control
2300 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
2301 TYPE(mp_para_env_type),
POINTER :: para_env
2302 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
2303 POINTER :: sab_all, sab_orb
2304 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
2305 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2309 CALL timeset(routinen, handle)
2311 NULLIFY (dft_control, mos, para_env, mo_coeff, op_mom_ao, &
2312 op_mom_der_ao, center_list, centers_set, &
2313 op_p_ao, psi1_p, psi1_rxp, psi1_d, &
2314 cell, particle_set, qs_kind_set)
2316 logger => cp_get_default_logger()
2317 output_unit = cp_logger_get_default_io_unit(logger)
2319 CALL get_qs_env(qs_env=qs_env, &
2320 dft_control=dft_control, &
2322 para_env=para_env, &
2324 dbcsr_dist=dbcsr_dist, &
2325 particle_set=particle_set, &
2326 qs_kind_set=qs_kind_set, &
2330 nspins = dft_control%nspins
2331 gapw = dft_control%qs_control%gapw
2333 CALL get_current_env(current_env=current_env, &
2336 nbr_center=nbr_center, &
2337 center_list=center_list, &
2338 centers_set=centers_set, &
2340 psi1_rxp=psi1_rxp, &
2343 psi0_order=psi0_order)
2347 DO ispin = 1, nspins
2348 DO icenter = 1, nbr_center(ispin)
2349 max_states = max(max_states, center_list(ispin)%array(1, icenter + 1)&
2350 & - center_list(ispin)%array(1, icenter))
2356 CALL dbcsr_allocate_matrix_set(op_mom_ao, 9)
2357 CALL dbcsr_allocate_matrix_set(op_mom_der_ao, 9, 3)
2360 natom =
SIZE(particle_set, 1)
2361 ALLOCATE (first_sgf(natom))
2362 ALLOCATE (last_sgf(natom))
2363 CALL get_particle_set(particle_set, qs_kind_set, &
2364 first_sgf=first_sgf, &
2366 ALLOCATE (row_blk_sizes(natom))
2367 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
2368 DEALLOCATE (first_sgf)
2369 DEALLOCATE (last_sgf)
2372 ALLOCATE (op_mom_ao(1)%matrix)
2373 CALL dbcsr_create(matrix=op_mom_ao(1)%matrix, &
2375 dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
2376 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2377 nze=0, mutable_work=.true.)
2378 CALL cp_dbcsr_alloc_block_from_nbl(op_mom_ao(1)%matrix, sab_all)
2381 ALLOCATE (op_mom_der_ao(1, idir2)%matrix)
2382 CALL dbcsr_copy(op_mom_der_ao(1, idir2)%matrix, op_mom_ao(1)%matrix, &
2383 "op_mom_der_ao"//
"-"//trim(adjustl(cp_to_string(idir2))))
2386 DO idir = 2,
SIZE(op_mom_ao, 1)
2387 ALLOCATE (op_mom_ao(idir)%matrix)
2388 CALL dbcsr_copy(op_mom_ao(idir)%matrix, op_mom_ao(1)%matrix, &
2389 "op_mom_ao"//
"-"//trim(adjustl(cp_to_string(idir))))
2391 ALLOCATE (op_mom_der_ao(idir, idir2)%matrix)
2392 CALL dbcsr_copy(op_mom_der_ao(idir, idir2)%matrix, op_mom_ao(1)%matrix, &
2393 "op_mom_der_ao"//
"-"//trim(adjustl(cp_to_string(idir*idir2))))
2397 CALL dbcsr_allocate_matrix_set(op_p_ao, 3)
2398 ALLOCATE (op_p_ao(1)%matrix)
2399 CALL dbcsr_create(matrix=op_p_ao(1)%matrix, &
2401 dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
2402 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2403 nze=0, mutable_work=.true.)
2404 CALL cp_dbcsr_alloc_block_from_nbl(op_p_ao(1)%matrix, sab_orb)
2407 ALLOCATE (op_p_ao(idir)%matrix)
2408 CALL dbcsr_copy(op_p_ao(idir)%matrix, op_p_ao(1)%matrix, &
2409 "op_p_ao"//
"-"//trim(adjustl(cp_to_string(idir))))
2413 DEALLOCATE (row_blk_sizes)
2417 mo_coeff => psi0_order(1)
2418 NULLIFY (tmp_fm_struct)
2419 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
2420 ncol_global=max_states, para_env=para_env, &
2421 context=mo_coeff%matrix_struct%context)
2422 CALL cp_fm_create(psi0, tmp_fm_struct)
2423 CALL cp_fm_create(psi_d, tmp_fm_struct)
2424 CALL cp_fm_create(psi_rxp, tmp_fm_struct)
2425 CALL cp_fm_create(psi_p1, tmp_fm_struct)
2426 CALL cp_fm_create(psi_p2, tmp_fm_struct)
2427 CALL cp_fm_struct_release(tmp_fm_struct)
2429 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
2430 ncol_global=max_states, para_env=para_env, &
2431 context=mo_coeff%matrix_struct%context)
2433 CALL cp_fm_create(p_rxp(idir), tmp_fm_struct)
2434 CALL cp_fm_create(r_p1(idir), tmp_fm_struct)
2435 CALL cp_fm_create(r_p2(idir), tmp_fm_struct)
2437 CALL cp_fm_create(rr_rxp(idir2, idir), tmp_fm_struct)
2438 CALL cp_fm_create(rr_p1(idir2, idir), tmp_fm_struct)
2439 CALL cp_fm_create(rr_p2(idir2, idir), tmp_fm_struct)
2442 CALL cp_fm_struct_release(tmp_fm_struct)
2447 CALL build_lin_mom_matrix(qs_env, op_p_ao)
2452 CALL set_vecp(ib, iib, iiib)
2453 DO ispin = 1, nspins
2456 nmo = nstates(ispin)
2457 mo_coeff => psi0_order(ispin)
2458 CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc)
2462 int_current = 0.0_dp
2465 DO icenter = 1, nbr_center(ispin)
2468 dk(1:3) = centers_set(ispin)%array(1:3, icenter)
2473 CALL dbcsr_set(op_mom_ao(idir)%matrix, 0.0_dp)
2475 CALL dbcsr_set(op_mom_der_ao(idir, idir2)%matrix, 0.0_dp)
2478 CALL rrc_xyz_der_ao(op_mom_ao, op_mom_der_ao, qs_env, dk, order=2, &
2479 minimum_image=.false., soft=gapw)
2482 CALL cp_fm_set_all(psi0, 0.0_dp)
2483 CALL cp_fm_set_all(psi_rxp, 0.0_dp)
2484 CALL cp_fm_set_all(psi_d, 0.0_dp)
2485 CALL cp_fm_set_all(psi_p1, 0.0_dp)
2486 CALL cp_fm_set_all(psi_p2, 0.0_dp)
2487 nstate_loc = center_list(ispin)%array(1, icenter + 1) - center_list(ispin)%array(1, icenter)
2489 DO j = center_list(ispin)%array(1, icenter), center_list(ispin)%array(1, icenter + 1) - 1
2490 istate = center_list(ispin)%array(2, j)
2493 CALL cp_fm_to_fm(mo_coeff, psi0, 1, istate, jstate)
2495 CALL cp_fm_to_fm(psi1_rxp(ispin, ib), psi_rxp, 1, istate, jstate)
2496 IF (current_env%full)
CALL cp_fm_to_fm(psi1_d(ispin, ib), psi_d, 1, istate, jstate)
2499 CALL cp_fm_to_fm(psi1_p(ispin, iib), psi_p1, 1, istate, jstate)
2500 CALL cp_fm_to_fm(psi1_p(ispin, iiib), psi_p2, 1, istate, jstate)
2506 IF (current_env%full)
CALL cp_fm_scale_and_add(1.0_dp, psi_rxp, -1.0_dp, psi_d)
2509 CALL set_vecp(idir, ii, iii)
2510 CALL cp_dbcsr_sm_fm_multiply(op_p_ao(idir)%matrix, psi_rxp, &
2511 p_rxp(idir), ncol=nstate_loc, alpha=1.e0_dp)
2512 IF (iiib .EQ. iii .OR. iiib .EQ. ii)
THEN
2513 CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix, psi_p1, &
2514 r_p1(idir), ncol=nstate_loc, alpha=1.e0_dp)
2516 IF (iib .EQ. iii .OR. iib .EQ. ii)
THEN
2517 CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix, psi_p2, &
2518 r_p2(idir), ncol=nstate_loc, alpha=1.e0_dp)
2521 IF (idir2 .EQ. ii .OR. idir2 .EQ. iii)
THEN
2522 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2, idir)%matrix, psi_rxp, &
2523 rr_rxp(idir2, idir), ncol=nstate_loc, alpha=1.e0_dp)
2526 IF (idir2 .EQ. ind_m2(ii, iiib) .OR. idir2 .EQ. ind_m2(iii, iiib) .OR. idir2 .EQ. iiib)
THEN
2527 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2, idir)%matrix, psi_p1, &
2528 rr_p1(idir2, idir), ncol=nstate_loc, alpha=1.e0_dp)
2531 IF (idir2 .EQ. ind_m2(ii, iib) .OR. idir2 .EQ. ind_m2(iii, iib) .OR. idir2 .EQ. iib)
THEN
2532 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2, idir)%matrix, psi_p2, &
2533 rr_p2(idir2, idir), ncol=nstate_loc, alpha=1.e0_dp)
2546 int_current_tmp = 0.0_dp
2549 CALL set_vecp(idir, ii, iii)
2554 CALL cp_fm_trace(psi0, rr_rxp(ii, iii), contrib)
2555 chi_tmp = chi_tmp + 2.0_dp*contrib
2558 CALL cp_fm_trace(psi0, rr_rxp(iii, ii), contrib)
2559 chi_tmp = chi_tmp - 2.0_dp*contrib
2564 CALL cp_fm_trace(psi0, p_rxp(iii), contrib)
2565 IF (.NOT. chi_pbc) chi_tmp = chi_tmp + 2.0_dp*dk(ii)*contrib
2566 int_current_tmp = int_current_tmp + 2.0_dp*contrib
2569 CALL cp_fm_trace(psi0, p_rxp(ii), contrib2)
2570 IF (.NOT. chi_pbc) chi_tmp = chi_tmp - 2.0_dp*dk(iii)*contrib2
2576 idir2 = ind_m2(ii, iib)
2577 CALL cp_fm_trace(psi0, rr_p2(idir2, iii), contrib)
2578 chi_tmp = chi_tmp - 2.0_dp*contrib
2580 IF (iib == iii)
THEN
2581 CALL cp_fm_trace(psi0, r_p2(ii), contrib2)
2582 chi_tmp = chi_tmp - contrib2
2586 idir2 = ind_m2(iii, iib)
2587 CALL cp_fm_trace(psi0, rr_p2(idir2, ii), contrib)
2588 chi_tmp = chi_tmp + 2.0_dp*contrib
2591 CALL cp_fm_trace(psi0, r_p2(iii), contrib2)
2592 chi_tmp = chi_tmp + contrib2
2600 CALL cp_fm_trace(psi0, rr_p2(iib, iii), contrib)
2601 IF (.NOT. chi_pbc) chi_tmp = chi_tmp - 2.0_dp*dk(ii)*contrib
2602 int_current_tmp = int_current_tmp - 2.0_dp*contrib
2605 CALL cp_fm_trace(psi0, rr_p2(iib, ii), contrib2)
2606 IF (.NOT. chi_pbc) chi_tmp = chi_tmp + 2.0_dp*dk(iii)*contrib2
2612 idir2 = ind_m2(ii, iiib)
2613 CALL cp_fm_trace(psi0, rr_p1(idir2, iii), contrib)
2614 chi_tmp = chi_tmp + 2.0_dp*contrib
2616 IF (iiib == iii)
THEN
2617 CALL cp_fm_trace(psi0, r_p1(ii), contrib2)
2618 chi_tmp = chi_tmp + contrib2
2622 idir2 = ind_m2(iii, iiib)
2623 CALL cp_fm_trace(psi0, rr_p1(idir2, ii), contrib)
2624 chi_tmp = chi_tmp - 2.0_dp*contrib
2626 IF (iiib == ii)
THEN
2627 CALL cp_fm_trace(psi0, r_p1(iii), contrib2)
2628 chi_tmp = chi_tmp - contrib2
2635 CALL cp_fm_trace(psi0, rr_p1(iiib, iii), contrib)
2636 IF (.NOT. chi_pbc) chi_tmp = chi_tmp + 2.0_dp*dk(ii)*contrib
2637 int_current_tmp = int_current_tmp + 2.0_dp*contrib
2640 CALL cp_fm_trace(psi0, rr_p1(iiib, ii), contrib2)
2641 IF (.NOT. chi_pbc) chi_tmp = chi_tmp - 2.0_dp*dk(iii)*contrib2
2644 chi(idir) = chi(idir) + maxocc*chi_tmp
2645 int_current(iii) = int_current(iii) + int_current_tmp
2651 current_env%chi_tensor(idir, ib, ispin) = current_env%chi_tensor(idir, ib, ispin) + &
2653 IF (output_unit > 0)
THEN
2664 CALL dbcsr_deallocate_matrix_set(op_mom_ao)
2665 CALL dbcsr_deallocate_matrix_set(op_mom_der_ao)
2666 CALL dbcsr_deallocate_matrix_set(op_p_ao)
2668 CALL cp_fm_release(psi0)
2669 CALL cp_fm_release(psi_rxp)
2670 CALL cp_fm_release(psi_d)
2671 CALL cp_fm_release(psi_p1)
2672 CALL cp_fm_release(psi_p2)
2674 CALL cp_fm_release(p_rxp(idir))
2675 CALL cp_fm_release(r_p1(idir))
2676 CALL cp_fm_release(r_p2(idir))
2678 CALL cp_fm_release(rr_rxp(idir2, idir))
2679 CALL cp_fm_release(rr_p1(idir2, idir))
2680 CALL cp_fm_release(rr_p2(idir2, idir))
2684 CALL timestop(handle)
2686 END SUBROUTINE current_build_chi_many_centers
2694 SUBROUTINE current_build_chi_one_center(current_env, qs_env, iB)
2696 TYPE(current_env_type) :: current_env
2697 TYPE(qs_environment_type),
POINTER :: qs_env
2698 INTEGER,
INTENT(IN) :: ib
2700 CHARACTER(LEN=*),
PARAMETER :: routinen =
'current_build_chi_one_center'
2702 INTEGER :: handle, idir, idir2, iib, iiib, ispin, jdir, jjdir, kdir, max_states, nao, natom, &
2703 nbr_center(2), nmo, nspins, nstates(2), output_unit
2704 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf
2705 INTEGER,
DIMENSION(:),
POINTER :: row_blk_sizes
2706 LOGICAL :: chi_pbc, gapw
2707 REAL(dp) :: chi(3), contrib, dk(3), int_current(3), &
2709 TYPE(cell_type),
POINTER :: cell
2710 TYPE(cp_2d_i_p_type),
DIMENSION(:),
POINTER :: center_list
2711 TYPE(cp_2d_r_p_type),
DIMENSION(:),
POINTER :: centers_set
2712 TYPE(cp_fm_type) :: buf
2713 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: psi0_order
2714 TYPE(cp_fm_type),
DIMENSION(:, :),
POINTER :: psi1_p, psi1_rxp
2715 TYPE(cp_fm_type),
POINTER :: mo_coeff
2716 TYPE(cp_logger_type),
POINTER :: logger
2717 TYPE(dbcsr_distribution_type),
POINTER :: dbcsr_dist
2718 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: op_mom_ao, op_p_ao
2719 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: op_mom_der_ao
2720 TYPE(dft_control_type),
POINTER :: dft_control
2721 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
2722 TYPE(mp_para_env_type),
POINTER :: para_env
2723 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
2724 POINTER :: sab_all, sab_orb
2725 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
2726 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2730 CALL timeset(routinen, handle)
2732 NULLIFY (dft_control, mos, para_env, mo_coeff, op_mom_ao, &
2733 op_mom_der_ao, center_list, centers_set, &
2734 op_p_ao, psi1_p, psi1_rxp, cell, psi0_order)
2736 logger => cp_get_default_logger()
2737 output_unit = cp_logger_get_default_io_unit(logger)
2739 CALL get_qs_env(qs_env=qs_env, &
2740 dft_control=dft_control, &
2742 para_env=para_env, &
2744 particle_set=particle_set, &
2745 qs_kind_set=qs_kind_set, &
2748 dbcsr_dist=dbcsr_dist)
2750 nspins = dft_control%nspins
2751 gapw = dft_control%qs_control%gapw
2753 CALL get_current_env(current_env=current_env, &
2756 nbr_center=nbr_center, &
2757 center_list=center_list, &
2758 centers_set=centers_set, &
2760 psi1_rxp=psi1_rxp, &
2762 psi0_order=psi0_order)
2764 max_states = maxval(nstates(1:nspins))
2768 CALL dbcsr_allocate_matrix_set(op_mom_ao, 9)
2769 CALL dbcsr_allocate_matrix_set(op_mom_der_ao, 9, 3)
2772 natom =
SIZE(particle_set, 1)
2773 ALLOCATE (first_sgf(natom))
2774 ALLOCATE (last_sgf(natom))
2775 CALL get_particle_set(particle_set, qs_kind_set, &
2776 first_sgf=first_sgf, &
2778 ALLOCATE (row_blk_sizes(natom))
2779 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
2780 DEALLOCATE (first_sgf)
2781 DEALLOCATE (last_sgf)
2784 ALLOCATE (op_mom_ao(1)%matrix)
2785 CALL dbcsr_create(matrix=op_mom_ao(1)%matrix, &
2787 dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
2788 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2789 nze=0, mutable_work=.true.)
2790 CALL cp_dbcsr_alloc_block_from_nbl(op_mom_ao(1)%matrix, sab_all)
2793 ALLOCATE (op_mom_der_ao(1, idir2)%matrix)
2794 CALL dbcsr_copy(op_mom_der_ao(1, idir2)%matrix, op_mom_ao(1)%matrix, &
2795 "op_mom_der_ao"//
"-"//trim(adjustl(cp_to_string(idir2))))
2798 DO idir = 2,
SIZE(op_mom_ao, 1)
2799 ALLOCATE (op_mom_ao(idir)%matrix)
2800 CALL dbcsr_copy(op_mom_ao(idir)%matrix, op_mom_ao(1)%matrix, &
2801 "op_mom_ao"//
"-"//trim(adjustl(cp_to_string(idir))))
2803 ALLOCATE (op_mom_der_ao(idir, idir2)%matrix)
2804 CALL dbcsr_copy(op_mom_der_ao(idir, idir2)%matrix, op_mom_ao(1)%matrix, &
2805 "op_mom_der_ao"//
"-"//trim(adjustl(cp_to_string(idir*idir2))))
2809 CALL dbcsr_allocate_matrix_set(op_p_ao, 3)
2810 ALLOCATE (op_p_ao(1)%matrix)
2811 CALL dbcsr_create(matrix=op_p_ao(1)%matrix, &
2813 dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
2814 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2815 nze=0, mutable_work=.true.)
2816 CALL cp_dbcsr_alloc_block_from_nbl(op_p_ao(1)%matrix, sab_orb)
2819 ALLOCATE (op_p_ao(idir)%matrix)
2820 CALL dbcsr_copy(op_p_ao(idir)%matrix, op_p_ao(1)%matrix, &
2821 "op_p_ao"//
"-"//trim(adjustl(cp_to_string(idir))))
2825 DEALLOCATE (row_blk_sizes)
2828 CALL build_lin_mom_matrix(qs_env, op_p_ao)
2833 CALL set_vecp(ib, iib, iiib)
2834 DO ispin = 1, nspins
2836 cpassert(nbr_center(ispin) == 1)
2839 nmo = nstates(ispin)
2840 mo_coeff => psi0_order(ispin)
2841 CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc)
2844 CALL cp_fm_create(buf, mo_coeff%matrix_struct)
2848 int_current = 0.0_dp
2852 dk(1:3) = centers_set(ispin)%array(1:3, 1)
2857 CALL dbcsr_set(op_mom_ao(idir)%matrix, 0.0_dp)
2859 CALL dbcsr_set(op_mom_der_ao(idir, idir2)%matrix, 0.0_dp)
2862 CALL rrc_xyz_der_ao(op_mom_ao, op_mom_der_ao, qs_env, dk, order=2, &
2863 minimum_image=.false., soft=gapw)
2877 IF (.NOT. chi_pbc)
THEN
2878 CALL cp_dbcsr_sm_fm_multiply(op_p_ao(idir)%matrix, mo_coeff, &
2879 buf, ncol=nmo, alpha=1.e0_dp)
2882 IF (levi_civita(kdir, jdir, idir) .EQ. 0.0_dp) cycle
2883 CALL cp_fm_trace(buf, psi1_rxp(ispin, ib), contrib)
2884 chi(kdir) = chi(kdir) - levi_civita(kdir, jdir, idir)*2.0_dp*dk(jdir)*contrib
2899 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(jdir, idir)%matrix, mo_coeff, &
2900 buf, ncol=nmo, alpha=1.e0_dp)
2902 IF (levi_civita(kdir, jdir, idir) .EQ. 0.0_dp) cycle
2903 CALL cp_fm_trace(buf, psi1_rxp(ispin, ib), contrib)
2904 chi(kdir) = chi(kdir) - levi_civita(kdir, jdir, idir)*2.0_dp*contrib
2907 IF (.NOT. chi_pbc)
THEN
2908 IF (jdir .EQ. iib)
THEN
2911 IF (levi_civita(kdir, jjdir, idir) .EQ. 0.0_dp) cycle
2912 CALL cp_fm_trace(buf, psi1_p(ispin, iiib), contrib)
2913 chi(kdir) = chi(kdir) + levi_civita(kdir, jjdir, idir)*2.0_dp*dk(jjdir)*contrib
2918 IF (jdir .EQ. iiib)
THEN
2921 IF (levi_civita(kdir, jjdir, idir) .EQ. 0.0_dp) cycle
2922 CALL cp_fm_trace(buf, psi1_p(ispin, iib), contrib)
2923 chi(kdir) = chi(kdir) - levi_civita(kdir, jjdir, idir)*2.0_dp*dk(jjdir)*contrib
2939 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(ind_m2(jdir, iib), idir)%matrix, mo_coeff, &
2940 buf, ncol=nmo, alpha=1.e0_dp)
2942 IF (levi_civita(kdir, jdir, idir) .EQ. 0.0_dp) cycle
2943 CALL cp_fm_trace(buf, psi1_p(ispin, iiib), contrib)
2944 chi(kdir) = chi(kdir) + levi_civita(kdir, jdir, idir)*2.0_dp*contrib
2947 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(ind_m2(jdir, iiib), idir)%matrix, mo_coeff, &
2948 buf, ncol=nmo, alpha=1.e0_dp)
2950 IF (levi_civita(kdir, jdir, idir) .EQ. 0.0_dp) cycle
2951 CALL cp_fm_trace(buf, psi1_p(ispin, iib), contrib)
2952 chi(kdir) = chi(kdir) - levi_civita(kdir, jdir, idir)*2.0_dp*contrib
2963 CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix, mo_coeff, &
2964 buf, ncol=nmo, alpha=1.e0_dp)
2967 IF (levi_civita(kdir, idir, jdir) .EQ. 0.0_dp) cycle
2968 IF (iib == jdir)
THEN
2969 CALL cp_fm_trace(buf, psi1_p(ispin, iiib), contrib)
2970 chi(kdir) = chi(kdir) + levi_civita(kdir, idir, jdir)*contrib
2977 IF (levi_civita(kdir, idir, jdir) .EQ. 0.0_dp) cycle
2978 IF (iiib == jdir)
THEN
2979 CALL cp_fm_trace(buf, psi1_p(ispin, iib), contrib)
2980 chi(kdir) = chi(kdir) - levi_civita(kdir, idir, jdir)*contrib
2992 current_env%chi_tensor(idir, ib, ispin) = current_env%chi_tensor(idir, ib, ispin) + &
2994 IF (output_unit > 0)
THEN
3000 CALL cp_fm_release(buf)
3004 CALL dbcsr_deallocate_matrix_set(op_mom_ao)
3005 CALL dbcsr_deallocate_matrix_set(op_mom_der_ao)
3006 CALL dbcsr_deallocate_matrix_set(op_p_ao)
3008 CALL timestop(handle)
3010 END SUBROUTINE current_build_chi_one_center
static int imax(int x, int y)
Returns the larger of two given integer (missing from the C standard)
static int imin(int x, int y)
Returns the smaller of two given integer (missing from the C standard)
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
static GRID_HOST_DEVICE int ncoset(const int l)
Number of Cartesian orbitals up to given angular momentum quantum.
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
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.
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Handles all functions related to the CELL.
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...
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_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
basic linear algebra operations for full matrices
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
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...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
for a given dr()/dh(r) this will provide the bounds to be used if one wants to go over a sphere-subre...
subroutine, public compute_cube_center(cube_center, rs_desc, zeta, zetb, ra, rab)
unifies the computation of the cube center, so that differences in implementation,...
subroutine, public return_cube(info, radius, lb_cube, ub_cube, sphere_bounds)
...
Fortran API for the grid package, which is written in C.
integer, parameter, public grid_func_adbmdab_z
integer, parameter, public grid_func_adbmdab_y
integer, parameter, public grid_func_ardbmdarb_yx
integer, parameter, public grid_func_ardbmdarb_zz
integer, parameter, public grid_func_ardbmdarb_xy
integer, parameter, public grid_func_ardbmdarb_zx
integer, parameter, public grid_func_ardbmdarb_xx
integer, parameter, public grid_func_ardbmdarb_yz
integer, parameter, public grid_func_ab
integer, parameter, public grid_func_ardbmdarb_yy
integer, parameter, public grid_func_adbmdab_x
subroutine, public collocate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, scale, pab, o1, o2, rsgrid, ga_gb_function, radius, use_subpatch, subpatch_pattern)
low level collocation of primitive gaussian functions
integer, parameter, public grid_func_ardbmdarb_zy
integer, parameter, public grid_func_ardbmdarb_xz
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Definition of mathematical constants and functions.
real(kind=dp), parameter, public twopi
Utility routines for the memory handling.
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
represent a simple array based list of the given type
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.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
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.
given the response wavefunctions obtained by the application of the (rxp), p, and ((dk-dl)xp) operato...
subroutine, public calculate_jrho_atom(current_env, qs_env, iB, idir)
...
subroutine, public calculate_jrho_atom_rad(qs_env, current_env, idir)
...
subroutine, public calculate_jrho_atom_coeff(qs_env, current_env, mat_d0, mat_jp, mat_jp_rii, mat_jp_riii, iB, idir)
Calculate the expansion coefficients for the atomic terms of the current densitiy in GAPW.
given the response wavefunctions obtained by the application of the (rxp), p, and ((dk-dl)xp) operato...
subroutine, public current_build_current(current_env, qs_env, iB)
First calculate the density matrixes, for each component of the current they are 3 because of the r d...
subroutine, public current_build_chi(current_env, qs_env, iB)
...
subroutine, public calculate_jrho_resp(mat_d0, mat_jp, mat_jp_rii, mat_jp_riii, iB, idir, current_rs, current_gs, qs_env, current_env, soft_valid, retain_rsgrid)
Calculation of the idir component of the response current density in the presence of a constant magne...
Calculate the operators p rxp and D needed in the optimization of the different contribution of the f...
integer function, public ind_m2(ii, iii)
...
subroutine, public set_vecp(i1, i2, i3)
...
real(dp) function, public fac_vecp(a, b, c)
...
subroutine, public fm_scale_by_pbc_ac(matrix, ra, rc, cell, ixyz)
scale a matrix as a_ij = a_ij * pbc(rc(:,j),ra(:,i))(ixyz)
subroutine, public set_vecp_rev(i1, i2, i3)
...
Type definitiona for linear response calculations.
subroutine, public get_current_env(current_env, simple_done, simple_converged, full_done, nao, nstates, gauge, list_cubes, statetrueindex, gauge_name, basisfun_center, nbr_center, center_list, centers_set, psi1_p, psi1_rxp, psi1_D, p_psi0, rxp_psi0, jrho1_atom_set, jrho1_set, chi_tensor, chi_tensor_loc, gauge_atom_radius, rs_gauge, use_old_gauge_atom, chi_pbc, psi0_order)
...
wrapper for the pools of matrixes
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)
...
subroutine, public rrc_xyz_der_ao(op, op_der, qs_env, rc, order, minimum_image, soft)
Calculation of the multipole operators integrals and of its derivatives of the type [\mu | op | d(\nu...
subroutine, public build_lin_mom_matrix(qs_env, matrix)
Calculation of the linear momentum matrix <mu|∂|nu> over Cartesian Gaussian functions.
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...
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
subroutine, public rs_grid_create(rs, desc)
...
subroutine, public rs_grid_mult_and_add(rs1, rs2, rs3, scalar)
rs1(i) = rs1(i) + rs2(i)*rs3(i)
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Transfers densities from PW to RS grids and potentials from PW to RS.
subroutine, public density_rs2pw(pw_env, rs_rho, rho, rho_gspace)
given partial densities on the realspace multigrids, computes the full density on the plane wave grid...
generate the tasks lists used by collocate and integrate routines
subroutine, public rs_distribute_matrix(rs_descs, pmats, atom_pair_send, atom_pair_recv, nimages, scatter, hmats)
redistributes the matrix so that it can be used in realspace operations i.e. according to the task li...
subroutine, public distribute_tasks(rs_descs, ntasks, natoms, tasks, atom_pair_send, atom_pair_recv, symmetric, reorder_rs_grid_ranks, skip_load_balance_distributed)
Assembles tasks to be performed on local grid.
subroutine, public task_list_inner_loop(tasks, ntasks, curr_tasks, rs_descs, dft_control, cube_info, gridlevel_info, cindex, iatom, jatom, rpgfa, rpgfb, zeta, zetb, kind_radius_b, set_radius_a, set_radius_b, ra, rab, la_max, la_min, lb_max, lb_min, npgfa, npgfb, nseta, nsetb)
...
subroutine, public reallocate_tasks(tasks, new_size)
Grow an array of tasks while preserving the existing entries.