334 adiabatic_rescale_factor, kind_set_external, &
335 rho_atom_set_external, xc_section_external, calculate_forces, &
336 composite_vxc_rho, composite_vxc_tau, composite_reference_active, &
337 direct_valence_atom_grid, atom_composite_grid)
340 LOGICAL,
INTENT(IN) :: energy_only
341 REAL(
dp),
INTENT(INOUT) :: exc1
342 REAL(
dp),
INTENT(IN),
OPTIONAL :: adiabatic_rescale_factor
344 POINTER :: kind_set_external
346 POINTER :: rho_atom_set_external
348 LOGICAL,
INTENT(IN),
OPTIONAL :: calculate_forces
350 POINTER :: composite_vxc_rho, composite_vxc_tau
351 LOGICAL,
INTENT(OUT),
OPTIONAL :: composite_reference_active
352 LOGICAL,
INTENT(IN),
OPTIONAL :: direct_valence_atom_grid, &
355 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_vxc_atom'
357 INTEGER :: adjoint_bin, adjoint_entry, adjoint_nbins, adjoint_nchannels, &
358 adjoint_tile_count(3), adjoint_tile_lower(3), adjoint_tile_upper(3), &
359 atom_composite_components, base_shift(3), bo(2), composite_descriptor_target_image, &
360 composite_image_periodicity(3), composite_local_atom, composite_local_natom, &
361 composite_nflat, composite_partition_target_image, composite_pw_nflat, composite_row, &
362 gapw_density_partition, gapw_representation, handle, ia, iat, iatom, icomponent, idir, &
363 ikind, image_i1, image_i2, image_i3, image_shell(3), image_shift(3), ir, ispin, iw, jdir, &
364 myfun, na, natom, nr, nspins
365 INTEGER :: num_pe, source_atom, target_atom, xc_deriv_method_id, xc_rho_smooth_id, zatom
366 INTEGER(KIND=int_8),
ALLOCATABLE,
DIMENSION(:) :: composite_atomic_grid_sizes, &
367 composite_local_grid_sizes
368 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: adjoint_bin_offsets, adjoint_bin_rows, &
369 composite_atom_end, composite_atom_kind, composite_atom_kind_index, composite_atom_start, &
370 composite_grid_atom, composite_local_atoms
371 INTEGER,
DIMENSION(2, 3) :: bounds
372 INTEGER,
DIMENSION(:),
POINTER :: atom_list
373 LOGICAL :: accint, atom_composite_active, atom_composite_diagnostic, &
374 atom_composite_reference, direct_valence_atom_composite, donlcc, evaluate_hard, &
375 evaluate_soft, gradient_f, image_partition_atom_composite, lsd, my_calculate_forces, &
376 native_grid_diagnostics, nlcc, one_center_kind, paw_atom, paw_pseudopotentials, &
377 requested_atom_composite_grid, rho_g_valid, skala_atom_grid, source_matrix_local, tau_f, &
378 tau_r_valid, use_atom_composite_density, use_atom_composite_gradient, &
379 use_atom_composite_tau, use_virial
380 LOGICAL,
ALLOCATABLE,
DIMENSION(:) :: composite_partition_included
381 REAL(
dp) :: agr, alpha, atom_composite_exc, atom_composite_nelec, &
382 composite_cross_cutoff_max, composite_cross_density_max, composite_cross_grad_max, &
383 composite_cross_kin_max, composite_density_max, composite_density_min, &
384 composite_grad_max, composite_kin_max, composite_kin_min, composite_tau_integral, &
385 cross_cutoff, density_cut, descriptor_window_adjoint, descriptor_window_weight, exc_h, &
386 exc_s, feature_vxc_analytic, feature_vxc_fd, feature_vxc_minus, feature_vxc_plus, &
387 feature_vxc_step, gradient_cut, local_partition_weight, my_adiabatic_rescale_factor, &
388 nlcc_density, nlcc_spin_factor
389 REAL(
dp) :: one_center_density_field_contraction, one_center_density_matrix_contraction, &
390 one_center_field_contraction, one_center_gradient_field_contraction, &
391 one_center_gradient_matrix_contraction, one_center_matrix_contraction, &
392 one_center_rho_grad_field_contraction, one_center_rho_grad_matrix_contraction, &
393 one_center_tau_field_contraction, one_center_tau_matrix_contraction, &
394 one_center_tensor_contraction, partition_adjoint, partition_scale, partition_weight, &
395 smooth_grid_contraction, smooth_input_contraction, target_partition_adjoint, tau_cut, zeff
396 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: composite_atomic_grid_weight_grad, &
397 composite_atomic_grid_weights, composite_base_grid_weights, composite_distances, &
398 composite_grid_weight_grad, composite_grid_weights, composite_partition_weights
399 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: composite_atom_coord_grad, composite_atom_coords, &
400 composite_cross_density, composite_cross_force, composite_cross_force_local, &
401 composite_cross_kin, composite_density, composite_density_grad, &
402 composite_descriptor_image_coords, composite_explicit_force, composite_grid_coord_force, &
403 composite_grid_coord_grad, composite_grid_coords, composite_kin, composite_kin_grad, &
404 composite_local_atom_coords, composite_model_atom_force, composite_moving_smooth_force, &
405 composite_nlcc_center_force, composite_nlcc_center_force_local, &
406 composite_nlcc_target_force
407 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: composite_partition_atom_coords, &
408 composite_partition_force, composite_partition_force_local, &
409 composite_partition_image_coords, composite_smooth_density_cache, &
410 composite_smooth_kin_cache, local_partition_datom
411 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :),
TARGET :: smooth_density_adjoint_storage, &
412 smooth_kin_adjoint_storage
413 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: composite_cross_grad, composite_grad, &
414 composite_grad_grad, composite_int_h, composite_int_s, composite_partition_datom, &
415 composite_partition_dstrain, composite_smooth_gradient_cache
416 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :),
TARGET :: smooth_grad_adjoint_storage
417 REAL(
dp),
DIMENSION(1, 1, 1) :: tau_d
418 REAL(
dp),
DIMENSION(1, 1, 1, 1) :: rho_d
419 REAL(
dp),
DIMENSION(2) :: composite_smooth_density_adjoint_value, &
420 composite_smooth_density_value, composite_smooth_kin_adjoint_value, &
421 composite_smooth_kin_value, cross_density, cross_density_adjoint, cross_kin, &
423 REAL(
dp),
DIMENSION(3) :: composite_point, cross_displacement, cross_spatial_derivative, &
424 fractional, image_translation, nlcc_gradient, nlcc_spatial_derivative, &
425 skala_atom_force_h, skala_atom_force_s, spatial_derivative
426 REAL(
dp),
DIMENSION(3, 1) :: local_descriptor_datom
427 REAL(
dp),
DIMENSION(3, 2) :: composite_smooth_gradient_adjoint_value, &
428 composite_smooth_gradient_value, cross_density_spatial, cross_grad, cross_grad_adjoint, &
430 REAL(
dp),
DIMENSION(3, 3) :: composite_cross_image_virial, &
431 composite_cross_image_virial_local, composite_explicit_virial, composite_feature_virial, &
432 composite_interpolation_virial, composite_partition_strain_virial, &
433 local_descriptor_dstrain, local_partition_dstrain, nlcc_hessian, skala_atom_virial, &
434 skala_atom_virial_h, skala_atom_virial_s
435 REAL(
dp),
DIMENSION(3, 3, 2) :: cross_grad_spatial
436 REAL(
dp),
DIMENSION(4) :: feature_component_analytic, &
438 REAL(
dp),
DIMENSION(:, :),
POINTER :: rho_nlcc, smooth_density_adjoint, &
439 smooth_kin_adjoint, weight_h, weight_s
440 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: composite_smooth_rho, composite_smooth_rhoa, &
441 composite_smooth_rhob, composite_smooth_tau, composite_smooth_tau_a, &
442 composite_smooth_tau_b, rho_h, rho_s, smooth_grad_adjoint, smooth_rho, smooth_rhoa, &
443 smooth_rhob, smooth_tau, smooth_tau_a, smooth_tau_b, tau_h, tau_s, vtau_h, vtau_s, vxc_h, &
445 REAL(
dp),
DIMENSION(:, :, :, :),
POINTER :: drho_h, drho_s, vxg_h, vxg_s
448 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: composite_smooth_drho, composite_smooth_drhoa, &
449 composite_smooth_drhob, smooth_drho, smooth_drhoa, smooth_drhob
456 TYPE(native_grid_interpolation_stencil_type) :: interpolation_stencil
461 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
POINTER :: smooth_rho_r, smooth_tau_r, &
462 smooth_vxc_rho, smooth_vxc_tau
464 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: my_kind_set
466 TYPE(
rho_atom_coeff),
DIMENSION(:),
POINTER :: cpc_h, cpc_s, dr_h, dr_s, int_hh, &
469 TYPE(
rho_atom_type),
DIMENSION(:),
POINTER :: my_rho_atom_set
474 TYPE(tau_basis_cache_type) :: tau_basis_cache
482 CALL timeset(routinen, handle)
485 NULLIFY (auxbas_pw_pool)
486 NULLIFY (my_kind_set)
487 NULLIFY (atomic_kind_set)
490 NULLIFY (gth_potential)
495 NULLIFY (particle_set)
499 NULLIFY (my_rho_atom_set)
501 NULLIFY (smooth_rho, smooth_rhoa, smooth_rhob, smooth_tau, smooth_tau_a, smooth_tau_b)
502 NULLIFY (composite_smooth_rho, composite_smooth_rhoa, composite_smooth_rhob, &
503 composite_smooth_tau, composite_smooth_tau_a, composite_smooth_tau_b)
504 NULLIFY (smooth_rho_g, smooth_rho_r, smooth_tau_r)
505 NULLIFY (smooth_vxc_rho, smooth_vxc_tau)
507 NULLIFY (smooth_drho(idir)%array, smooth_drhoa(idir)%array, smooth_drhob(idir)%array)
508 NULLIFY (composite_smooth_drho(idir)%array, &
509 composite_smooth_drhoa(idir)%array, &
510 composite_smooth_drhob(idir)%array)
512 NULLIFY (sgp_potential)
514 my_calculate_forces = .false.
515 IF (
PRESENT(calculate_forces)) my_calculate_forces = calculate_forces
516 IF (
PRESENT(composite_reference_active)) composite_reference_active = .false.
517 direct_valence_atom_composite = .false.
518 IF (
PRESENT(direct_valence_atom_grid))
THEN
519 direct_valence_atom_composite = direct_valence_atom_grid
521 requested_atom_composite_grid = .false.
522 IF (
PRESENT(atom_composite_grid)) requested_atom_composite_grid = atom_composite_grid
524 IF (
PRESENT(adiabatic_rescale_factor))
THEN
525 my_adiabatic_rescale_factor = adiabatic_rescale_factor
527 my_adiabatic_rescale_factor = 1.0_dp
531 dft_control=dft_control, &
534 atomic_kind_set=atomic_kind_set, &
535 qs_kind_set=my_kind_set, &
537 particle_set=particle_set, &
540 rho_atom_set=my_rho_atom_set, &
543 IF (dft_control%qs_control%gapw_xc)
THEN
544 CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct)
546 CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
549 IF (
PRESENT(kind_set_external)) my_kind_set => kind_set_external
550 IF (
PRESENT(rho_atom_set_external)) my_rho_atom_set => rho_atom_set_external
553 accint = dft_control%qs_control%gapw_control%accurate_xcint
557 IF (
PRESENT(xc_section_external)) my_xc_section => xc_section_external
564 atom_composite_diagnostic = .false.
565 atom_composite_reference = .false.
566 paw_pseudopotentials = .false.
567 native_grid_diagnostics = .false.
568 atom_composite_components = 1
569 feature_vxc_step = 3.0e-3_dp
570 IF (skala_atom_grid)
THEN
572 cpassert(
ASSOCIATED(gauxc_section))
574 i_val=gapw_representation)
576 l_val=native_grid_diagnostics)
578 IF (skala_atom_grid)
THEN
579 DO ikind = 1,
SIZE(my_kind_set)
580 NULLIFY (gth_potential, sgp_potential)
581 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
582 gth_potential=gth_potential, sgp_potential=sgp_potential)
583 paw_pseudopotentials = paw_pseudopotentials .OR. &
584 (paw_atom .AND. (
ASSOCIATED(gth_potential) .OR. &
585 ASSOCIATED(sgp_potential)))
590 "NATIVE_GRID_GAPW_ATOM_COMPOSITE_DIAGNOSTIC", &
591 l_val=atom_composite_diagnostic)
593 "NATIVE_GRID_GAPW_ATOM_COMPOSITE_REFERENCE", &
594 l_val=atom_composite_reference)
596 "NATIVE_GRID_GAPW_ATOM_COMPOSITE_COMPONENTS", &
597 i_val=atom_composite_components)
599 "NATIVE_GRID_GAPW_ATOM_COMPOSITE_FD_STEP", &
600 r_val=feature_vxc_step)
602 atom_composite_reference = atom_composite_reference .OR. &
604 paw_pseudopotentials)
605 atom_composite_reference = atom_composite_reference .OR. requested_atom_composite_grid
606 atom_composite_reference = atom_composite_reference .OR. direct_valence_atom_composite
607 atom_composite_active = atom_composite_diagnostic .OR. atom_composite_reference
608 use_atom_composite_density = atom_composite_components <= 2
609 use_atom_composite_gradient = atom_composite_components <= 2
610 use_atom_composite_tau = atom_composite_components == 1 .OR. &
611 atom_composite_components == 3
612 IF (atom_composite_active)
THEN
613 CALL ensure_native_skala_atom_grids(my_kind_set, dft_control)
615 IF (direct_valence_atom_composite)
THEN
616 use_atom_composite_density = .false.
617 use_atom_composite_gradient = .false.
618 use_atom_composite_tau = .false.
622 image_partition_atom_composite = atom_composite_active
623 composite_image_periodicity = 1
624 IF (
PRESENT(composite_reference_active)) composite_reference_active = atom_composite_reference
626 IF (skala_atom_grid)
THEN
629 use_virial =
ASSOCIATED(virial)
630 IF (use_virial) use_virial = my_calculate_forces .AND. &
631 virial%pv_calculate .AND. (.NOT. virial%pv_numer)
635 my_rho_atom_set(:)%exc_h = 0.0_dp
636 my_rho_atom_set(:)%exc_s = 0.0_dp
645 lsd = dft_control%lsd
646 nspins = dft_control%nspins
649 calc_potential=.true.)
651 gradient_f = (needs%drho .OR. needs%drho_spin) .OR. skala_atom_grid
652 tau_f = (needs%tau .OR. needs%tau_spin) .OR. skala_atom_grid
654 IF (atom_composite_active)
THEN
656 needs%rho_spin = .true.
657 needs%drho_spin = .true.
658 needs%tau_spin = .true.
665 ALLOCATE (composite_atomic_grid_sizes(
SIZE(particle_set)), &
666 composite_atom_kind(
SIZE(particle_set)), &
667 composite_atom_kind_index(
SIZE(particle_set)), &
668 composite_atom_start(
SIZE(particle_set)), &
669 composite_atom_end(
SIZE(particle_set)), &
670 composite_atom_coords(3,
SIZE(particle_set)), &
671 composite_partition_weights(
SIZE(particle_set)), &
672 composite_partition_atom_coords(3,
SIZE(particle_set)), &
673 composite_distances(
SIZE(particle_set)))
674 composite_atomic_grid_sizes = 0_int_8
675 composite_atom_kind = 0
676 composite_atom_kind_index = 0
677 composite_atom_start = 0
678 composite_atom_end = 0
679 DO iatom = 1,
SIZE(particle_set)
680 composite_atom_coords(:, iatom) = particle_set(iatom)%r
682 DO ikind = 1,
SIZE(atomic_kind_set)
683 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
684 NULLIFY (gth_potential, sgp_potential)
685 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
686 gth_potential=gth_potential, grid_atom=grid_atom, &
687 sgp_potential=sgp_potential, zatom=zatom, zeff=zeff)
689 iatom = atom_list(iat)
690 composite_atomic_grid_sizes(iatom) = int(grid_atom%nr*grid_atom%ng_sphere, kind=
int_8)
691 composite_atom_kind(iatom) = ikind
692 composite_atom_kind_index(iatom) = iat
695 IF (any(composite_atomic_grid_sizes <= 0_int_8))
THEN
696 CALL cp_abort(__location__, &
697 "The atom-composite diagnostic requires a GAPW one-center grid for every atom.")
700 composite_local_natom = 0
701 DO ikind = 1,
SIZE(atomic_kind_set)
702 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
703 NULLIFY (gth_potential, sgp_potential)
704 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
705 gth_potential=gth_potential, sgp_potential=sgp_potential, &
706 zatom=zatom, zeff=zeff)
707 bo =
get_limit(natom, para_env%num_pe, para_env%mepos)
708 composite_local_natom = composite_local_natom + max(0, bo(2) - bo(1) + 1)
710 ALLOCATE (composite_local_atoms(composite_local_natom), &
711 composite_local_grid_sizes(composite_local_natom), &
712 composite_local_atom_coords(3, composite_local_natom))
713 composite_local_atom = 0
715 DO ikind = 1,
SIZE(atomic_kind_set)
716 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
717 NULLIFY (gth_potential, sgp_potential)
718 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
719 gth_potential=gth_potential, sgp_potential=sgp_potential, &
720 zatom=zatom, zeff=zeff)
721 bo =
get_limit(natom, para_env%num_pe, para_env%mepos)
722 DO iat = bo(1), bo(2)
723 iatom = atom_list(iat)
724 composite_local_atom = composite_local_atom + 1
725 composite_local_atoms(composite_local_atom) = iatom
726 composite_local_grid_sizes(composite_local_atom) = &
727 composite_atomic_grid_sizes(iatom)
728 composite_local_atom_coords(:, composite_local_atom) = &
729 composite_atom_coords(:, iatom)
730 composite_atom_start(iatom) = composite_nflat + 1
731 composite_nflat = composite_nflat + int(composite_atomic_grid_sizes(iatom))
732 composite_atom_end(iatom) = composite_nflat
735 cpassert(composite_local_atom == composite_local_natom)
736 ALLOCATE (composite_density(composite_nflat, 2), &
737 composite_grad(composite_nflat, 3, 2), &
738 composite_kin(composite_nflat, 2), &
739 composite_grid_atom(composite_nflat), &
740 composite_smooth_density_cache(composite_nflat, 2), &
741 composite_smooth_gradient_cache(composite_nflat, 3, 2), &
742 composite_smooth_kin_cache(composite_nflat, 2), &
743 composite_grid_coords(3, composite_nflat), &
744 composite_grid_weights(composite_nflat), &
745 composite_base_grid_weights(composite_nflat), &
746 composite_atomic_grid_weights(composite_nflat))
747 composite_density = 0.0_dp
748 composite_grad = 0.0_dp
749 composite_kin = 0.0_dp
750 DO composite_local_atom = 1, composite_local_natom
751 iatom = composite_local_atoms(composite_local_atom)
752 composite_grid_atom(composite_atom_start(iatom):composite_atom_end(iatom)) = iatom
754 composite_smooth_density_cache = 0.0_dp
755 composite_smooth_gradient_cache = 0.0_dp
756 composite_smooth_kin_cache = 0.0_dp
757 composite_grid_coords = 0.0_dp
758 composite_grid_weights = 0.0_dp
759 composite_base_grid_weights = 0.0_dp
760 composite_atomic_grid_weights = 0.0_dp
762 CALL qs_rho_get(rho_struct, rho_r=smooth_rho_r, rho_g=smooth_rho_g, &
763 tau_r=smooth_tau_r, rho_g_valid=rho_g_valid, &
764 tau_r_valid=tau_r_valid)
765 cpassert(rho_g_valid)
766 cpassert(tau_r_valid)
767 cpassert(
ASSOCIATED(smooth_rho_r))
768 cpassert(
ASSOCIATED(smooth_rho_g))
769 cpassert(
ASSOCIATED(smooth_tau_r))
770 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
772 i_val=xc_deriv_method_id)
774 i_val=xc_rho_smooth_id)
779 CALL xc_rho_set_update(smooth_rho_set, smooth_rho_r, smooth_rho_g, smooth_tau_r, needs, &
780 xc_deriv_method_id, xc_rho_smooth_id, auxbas_pw_pool)
782 CALL xc_rho_set_get(smooth_rho_set, rhoa=smooth_rhoa, rhob=smooth_rhob, &
783 drhoa=smooth_drhoa, drhob=smooth_drhob, &
784 tau_a=smooth_tau_a, tau_b=smooth_tau_b)
785 CALL gather_native_grid_field(smooth_rhoa, smooth_rho_r(1)%pw_grid, para_env, &
786 composite_smooth_rhoa)
787 CALL gather_native_grid_field(smooth_rhob, smooth_rho_r(1)%pw_grid, para_env, &
788 composite_smooth_rhob)
789 CALL gather_native_grid_field(smooth_tau_a, smooth_rho_r(1)%pw_grid, para_env, &
790 composite_smooth_tau_a)
791 CALL gather_native_grid_field(smooth_tau_b, smooth_rho_r(1)%pw_grid, para_env, &
792 composite_smooth_tau_b)
794 CALL gather_native_grid_field(smooth_drhoa(idir)%array, &
795 smooth_rho_r(1)%pw_grid, para_env, &
796 composite_smooth_drhoa(idir)%array)
797 CALL gather_native_grid_field(smooth_drhob(idir)%array, &
798 smooth_rho_r(1)%pw_grid, para_env, &
799 composite_smooth_drhob(idir)%array)
802 CALL xc_rho_set_get(smooth_rho_set, rho=smooth_rho, drho=smooth_drho, &
804 CALL gather_native_grid_field(smooth_rho, smooth_rho_r(1)%pw_grid, para_env, &
805 composite_smooth_rho)
806 CALL gather_native_grid_field(smooth_tau, smooth_rho_r(1)%pw_grid, para_env, &
807 composite_smooth_tau)
809 CALL gather_native_grid_field(smooth_drho(idir)%array, &
810 smooth_rho_r(1)%pw_grid, para_env, &
811 composite_smooth_drho(idir)%array)
820 NULLIFY (rho_h, drho_h, rho_s, drho_s, weight_h, weight_s)
821 NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
822 NULLIFY (tau_h, tau_s)
823 NULLIFY (vtau_h, vtau_s)
827 DO ikind = 1,
SIZE(atomic_kind_set)
828 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
829 NULLIFY (gth_potential, sgp_potential)
830 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
831 gth_potential=gth_potential, harmonics=harmonics, &
832 grid_atom=grid_atom, sgp_potential=sgp_potential, &
833 zatom=zatom, zeff=zeff)
834 one_center_kind = .NOT. direct_valence_atom_composite
835 IF (one_center_kind)
THEN
836 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type=
"GAPW_1C")
837 one_center_kind = paw_atom
838 IF (skala_atom_grid)
THEN
839 one_center_kind = native_skala_uses_one_center_kind( &
840 paw_atom, gapw_representation, &
841 ASSOCIATED(gth_potential) .OR.
ASSOCIATED(sgp_potential), &
845 IF (.NOT. one_center_kind .AND. .NOT. atom_composite_active) cycle
848 na = grid_atom%ng_sphere
850 IF (one_center_kind)
THEN
862 weight_h => grid_atom%weight
863 alpha = dft_control%qs_control%gapw_control%aw(ikind)
864 IF (
ASSOCIATED(grid_atom%gapw_weight_s))
THEN
865 IF (grid_atom%gapw_weight_alpha /= alpha)
DEALLOCATE (grid_atom%gapw_weight_s)
867 IF (.NOT.
ASSOCIATED(grid_atom%gapw_weight_s))
THEN
868 ALLOCATE (grid_atom%gapw_weight_s(na, nr))
870 agr = 1.0_dp - exp(-alpha*grid_atom%rad2(ir))
871 grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
873 grid_atom%gapw_weight_alpha = alpha
875 weight_s => grid_atom%gapw_weight_s
877 weight_h => grid_atom%weight
878 weight_s => grid_atom%weight
885 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
887 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
893 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
894 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
895 CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
896 CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
899 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
900 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
901 CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
902 CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
906 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
907 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
908 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
909 CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
910 CALL reallocate(vtau_s, 1, na, 1, nr, 1, nspins)
917 rho_nlcc => my_kind_set(ikind)%nlcc_pot
918 IF (
ASSOCIATED(rho_nlcc)) donlcc = .true.
924 num_pe = para_env%num_pe
925 bo =
get_limit(natom, para_env%num_pe, para_env%mepos)
927 DO iat = bo(1), bo(2)
928 iatom = atom_list(iat)
930 IF (one_center_kind)
THEN
931 my_rho_atom_set(iatom)%exc_h = 0.0_dp
932 my_rho_atom_set(iatom)%exc_s = 0.0_dp
934 rho_atom => my_rho_atom_set(iatom)
938 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
940 rho_rad_s=r_s, drho_rad_h=dr_h, &
941 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
947 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
951 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
958 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
959 r_h_d, r_s_d, drho_h, drho_s)
961 CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
962 ir, rho_nlcc(:, 1), rho_h, rho_s, &
963 rho_nlcc(:, 2), drho_h, drho_s)
968 IF (atom_composite_active)
THEN
969 IF (image_partition_atom_composite)
THEN
971 composite_atom_coords, cell, iatom, composite_image_periodicity, &
972 composite_partition_image_coords, composite_partition_target_image)
974 composite_atom_coords(:, iatom:iatom), cell, 1, &
975 composite_image_periodicity, composite_descriptor_image_coords, &
976 composite_descriptor_target_image)
1001 composite_row = composite_atom_start(iatom) + (ir - 1)*na + ia - 1
1002 composite_point(1) = particle_set(iatom)%r(1) + grid_atom%rad(ir)* &
1003 grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)
1004 composite_point(2) = particle_set(iatom)%r(2) + grid_atom%rad(ir)* &
1005 grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)
1006 composite_point(3) = particle_set(iatom)%r(3) + &
1007 grid_atom%rad(ir)*grid_atom%cos_pol(ia)
1008 composite_grid_coords(:, composite_row) = composite_point
1009 IF (image_partition_atom_composite)
THEN
1011 composite_point, composite_partition_image_coords, &
1012 composite_partition_target_image, partition_weight)
1016 composite_point, composite_descriptor_image_coords, &
1017 composite_descriptor_target_image, descriptor_window_weight)
1019 descriptor_window_weight)
1022 composite_point, composite_atom_coords, cell, &
1023 composite_partition_weights, composite_partition_atom_coords, &
1024 composite_distances)
1025 partition_weight = composite_partition_weights(iatom)
1026 partition_scale = 1.0_dp
1028 composite_base_grid_weights(composite_row) = grid_atom%weight(ia, ir)
1029 composite_atomic_grid_weights(composite_row) = &
1030 composite_base_grid_weights(composite_row)*partition_scale
1031 composite_grid_weights(composite_row) = &
1032 composite_base_grid_weights(composite_row)* &
1034 CALL create_native_grid_interpolation_stencil( &
1035 interpolation_stencil, smooth_rho_r(1)%pw_grid, cell, composite_point, &
1036 image_partition_atom_composite)
1038 CALL interpolate_native_grid_fields( &
1039 composite_smooth_rhoa, composite_smooth_drhoa(1)%array, &
1040 composite_smooth_drhoa(2)%array, composite_smooth_drhoa(3)%array, &
1041 composite_smooth_tau_a, interpolation_stencil, &
1042 composite_smooth_density_value(1), &
1043 composite_smooth_gradient_value(:, 1), &
1044 composite_smooth_kin_value(1))
1045 CALL interpolate_native_grid_fields( &
1046 composite_smooth_rhob, composite_smooth_drhob(1)%array, &
1047 composite_smooth_drhob(2)%array, composite_smooth_drhob(3)%array, &
1048 composite_smooth_tau_b, interpolation_stencil, &
1049 composite_smooth_density_value(2), &
1050 composite_smooth_gradient_value(:, 2), &
1051 composite_smooth_kin_value(2))
1052 composite_smooth_density_cache(composite_row, :) = &
1053 composite_smooth_density_value
1054 composite_smooth_gradient_cache(composite_row, :, :) = &
1055 composite_smooth_gradient_value
1056 composite_smooth_kin_cache(composite_row, :) = &
1057 composite_smooth_kin_value
1059 composite_density(composite_row, ispin) = &
1060 composite_smooth_density_value(ispin)
1061 composite_grad(composite_row, :, ispin) = &
1062 composite_smooth_gradient_value(:, ispin)
1063 composite_kin(composite_row, ispin) = &
1064 composite_smooth_kin_value(ispin)
1065 IF (one_center_kind .AND. use_atom_composite_density)
THEN
1066 composite_density(composite_row, ispin) = &
1067 composite_density(composite_row, ispin) + &
1068 rho_h(ia, ir, ispin) - rho_s(ia, ir, ispin)
1070 IF (one_center_kind .AND. use_atom_composite_gradient)
THEN
1072 composite_grad(composite_row, idir, ispin) = &
1073 composite_grad(composite_row, idir, ispin) + &
1074 drho_h(idir, ia, ir, ispin) - drho_s(idir, ia, ir, ispin)
1077 IF (one_center_kind .AND. use_atom_composite_tau)
THEN
1078 composite_kin(composite_row, ispin) = &
1079 composite_kin(composite_row, ispin) + &
1080 tau_h(ia, ir, ispin) - tau_s(ia, ir, ispin)
1084 CALL interpolate_native_grid_fields( &
1085 composite_smooth_rho, composite_smooth_drho(1)%array, &
1086 composite_smooth_drho(2)%array, composite_smooth_drho(3)%array, &
1087 composite_smooth_tau, interpolation_stencil, &
1088 composite_smooth_density_value(1), &
1089 composite_smooth_gradient_value(:, 1), &
1090 composite_smooth_kin_value(1))
1091 composite_smooth_density_cache(composite_row, 1) = &
1092 composite_smooth_density_value(1)
1093 composite_smooth_gradient_cache(composite_row, :, 1) = &
1094 composite_smooth_gradient_value(:, 1)
1095 composite_smooth_kin_cache(composite_row, 1) = &
1096 composite_smooth_kin_value(1)
1097 composite_density(composite_row, :) = &
1098 0.5_dp*composite_smooth_density_value(1)
1100 composite_grad(composite_row, idir, :) = &
1101 0.5_dp*composite_smooth_gradient_value(idir, 1)
1103 composite_kin(composite_row, :) = &
1104 0.5_dp*composite_smooth_kin_value(1)
1105 IF (one_center_kind .AND. use_atom_composite_density)
THEN
1106 composite_density(composite_row, :) = &
1107 composite_density(composite_row, :) + &
1108 0.5_dp*(rho_h(ia, ir, 1) - rho_s(ia, ir, 1))
1110 IF (one_center_kind .AND. use_atom_composite_gradient)
THEN
1112 composite_grad(composite_row, idir, :) = &
1113 composite_grad(composite_row, idir, :) + &
1114 0.5_dp*(drho_h(idir, ia, ir, 1) - &
1115 drho_s(idir, ia, ir, 1))
1118 IF (one_center_kind .AND. use_atom_composite_tau)
THEN
1119 composite_kin(composite_row, :) = composite_kin(composite_row, :) + &
1120 0.5_dp*(tau_h(ia, ir, 1) - &
1124 IF (atom_composite_reference .AND. nlcc)
THEN
1125 nlcc_spin_factor = merge(1.0_dp, 0.5_dp, lsd)
1126 DO source_atom = 1,
SIZE(particle_set)
1127 NULLIFY (gth_potential, sgp_potential)
1128 CALL get_qs_kind(my_kind_set(composite_atom_kind(source_atom)), &
1129 gth_potential=gth_potential, &
1130 sgp_potential=sgp_potential)
1132 composite_point, particle_set(source_atom)%r, &
1133 gth_potential, sgp_potential, nlcc_density, &
1134 nlcc_gradient, nlcc_hessian)
1135 composite_density(composite_row, :) = &
1136 composite_density(composite_row, :) + &
1137 nlcc_spin_factor*nlcc_density
1139 composite_grad(composite_row, idir, :) = &
1140 composite_grad(composite_row, idir, :) + &
1141 nlcc_spin_factor*nlcc_gradient(idir)
1148 IF (image_partition_atom_composite)
THEN
1149 DEALLOCATE (composite_descriptor_image_coords, composite_partition_image_coords)
1151 cpassert(nr*na == composite_atom_end(iatom) - composite_atom_start(iatom) + 1)
1154 IF (.NOT. one_center_kind) cycle
1158 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
1159 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
1160 ELSE IF (gradient_f)
THEN
1161 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
1162 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
1164 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
1165 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
1169 evaluate_hard = .true.
1170 evaluate_soft = .true.
1171 skala_atom_force_h = 0.0_dp
1172 skala_atom_force_s = 0.0_dp
1173 skala_atom_virial_h = 0.0_dp
1174 skala_atom_virial_s = 0.0_dp
1175 IF (skala_atom_grid)
THEN
1176 SELECT CASE (gapw_density_partition)
1180 evaluate_soft = .false.
1182 evaluate_hard = .false.
1184 evaluate_hard = .false.
1185 evaluate_soft = .false.
1187 CALL cp_abort(__location__, &
1188 "Unknown GAUXC%NATIVE_GRID_GAPW_DENSITY_PARTITION value.")
1191 IF (atom_composite_reference)
THEN
1192 evaluate_hard = .false.
1193 evaluate_soft = .false.
1200 IF (.NOT. evaluate_hard)
THEN
1202 IF (.NOT. energy_only)
THEN
1204 IF (
ASSOCIATED(vxg_h)) vxg_h = 0.0_dp
1205 IF (
ASSOCIATED(vtau_h)) vtau_h = 0.0_dp
1207 ELSE IF (skala_atom_grid)
THEN
1209 my_xc_section, grid_atom, para_env, particle_set(iatom)%r, &
1210 rho_h, drho_h, tau_h, weight_h, lsd, nspins, na, nr, &
1211 exc_h, vxc_h, vxg_h, vtau_h, energy_only=energy_only, &
1212 atom_force=skala_atom_force_h, atom_virial=skala_atom_virial_h)
1214 CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight_h, &
1215 lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h, energy_only=energy_only, &
1216 adiabatic_rescale_factor=my_adiabatic_rescale_factor)
1218 rho_atom%exc_h = rho_atom%exc_h + exc_h
1224 IF (.NOT. evaluate_soft)
THEN
1226 IF (.NOT. energy_only)
THEN
1228 IF (
ASSOCIATED(vxg_s)) vxg_s = 0.0_dp
1229 IF (
ASSOCIATED(vtau_s)) vtau_s = 0.0_dp
1231 ELSE IF (skala_atom_grid)
THEN
1233 my_xc_section, grid_atom, para_env, particle_set(iatom)%r, &
1234 rho_s, drho_s, tau_s, weight_s, lsd, nspins, na, nr, &
1235 exc_s, vxc_s, vxg_s, vtau_s, energy_only=energy_only, &
1236 atom_force=skala_atom_force_s, atom_virial=skala_atom_virial_s)
1238 CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight_s, &
1239 lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s, energy_only=energy_only, &
1240 adiabatic_rescale_factor=my_adiabatic_rescale_factor)
1242 rho_atom%exc_s = rho_atom%exc_s + exc_s
1246 exc1 = exc1 + rho_atom%exc_h - rho_atom%exc_s
1247 IF (skala_atom_grid .AND. my_calculate_forces .AND.
ASSOCIATED(force))
THEN
1248 force(ikind)%rho_elec(:, iat) = force(ikind)%rho_elec(:, iat) + &
1249 skala_atom_force_h - skala_atom_force_s
1251 IF (skala_atom_grid .AND. use_virial)
THEN
1252 skala_atom_virial = skala_atom_virial_h - skala_atom_virial_s
1255 virial%pv_gapw(idir, jdir) = virial%pv_gapw(idir, jdir) + &
1256 skala_atom_virial(idir, jdir)
1257 virial%pv_virial(idir, jdir) = virial%pv_virial(idir, jdir) + &
1258 skala_atom_virial(idir, jdir)
1267 IF (.NOT. energy_only)
THEN
1268 NULLIFY (int_hh, int_ss)
1269 CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1270 IF (gradient_f)
THEN
1271 CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
1272 grid_atom, basis_1c, harmonics, nspins)
1275 grid_atom, basis_1c, harmonics, nspins)
1278 CALL dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, tau_basis_cache, nspins)
1281 NULLIFY (r_h, r_s, dr_h, dr_s)
1284 IF (one_center_kind)
THEN
1285 IF (tau_f)
CALL release_tau_basis_cache(tau_basis_cache)
1293 IF (atom_composite_active)
THEN
1294 ALLOCATE (composite_cross_density(composite_nflat, 2), &
1295 composite_cross_grad(composite_nflat, 3, 2), &
1296 composite_cross_kin(composite_nflat, 2))
1297 composite_cross_density = 0.0_dp
1298 composite_cross_grad = 0.0_dp
1299 composite_cross_kin = 0.0_dp
1300 composite_cross_cutoff_max = 0.0_dp
1301 IF (.NOT. direct_valence_atom_composite)
THEN
1302 DO ikind = 1,
SIZE(atomic_kind_set)
1303 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1304 NULLIFY (gth_potential, sgp_potential)
1305 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
1306 gth_potential=gth_potential, harmonics=harmonics, &
1307 grid_atom=grid_atom, sgp_potential=sgp_potential, &
1308 zatom=zatom, zeff=zeff)
1309 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type=
"GAPW_1C")
1310 IF (.NOT. native_skala_uses_one_center_kind( &
1311 paw_atom, gapw_representation, &
1312 ASSOCIATED(gth_potential) .OR.
ASSOCIATED(sgp_potential), &
1316 para_env, my_rho_atom_set, my_kind_set(ikind), atom_list, natom, nspins)
1318 na = grid_atom%ng_sphere
1319 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
1320 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
1321 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
1322 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
1323 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
1324 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
1325 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
1330 source_atom = atom_list(iat)
1331 rho_atom => my_rho_atom_set(source_atom)
1332 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1333 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s, &
1334 drho_rad_h=dr_h, drho_rad_s=dr_s, &
1335 rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1340 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
1343 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
1344 r_h_d, r_s_d, drho_h, drho_s)
1347 cross_cutoff = gapw_atom_grid_support_radius( &
1348 grid_atom, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s)
1349 IF (cross_cutoff <= 0.0_dp) cycle
1350 composite_cross_cutoff_max = max(composite_cross_cutoff_max, cross_cutoff)
1353 IF (cell%perd(idir) == 1)
THEN
1354 image_shell(idir) = ceiling( &
1355 cross_cutoff*sqrt(sum(cell%h_inv(idir, :)**2))) + 1
1369 DO composite_row = 1, composite_nflat
1370 target_atom = composite_grid_atom(composite_row)
1374 fractional(idir) = fractional(idir) + cell%h_inv(idir, jdir)* &
1375 (composite_grid_coords(jdir, composite_row) - &
1376 particle_set(source_atom)%r(jdir))
1380 base_shift(idir) = cell%perd(idir)*nint(fractional(idir))
1382 DO image_i3 = base_shift(3) - image_shell(3), &
1383 base_shift(3) + image_shell(3)
1384 DO image_i2 = base_shift(2) - image_shell(2), &
1385 base_shift(2) + image_shell(2)
1386 DO image_i1 = base_shift(1) - image_shell(1), &
1387 base_shift(1) + image_shell(1)
1388 image_shift = [image_i1, image_i2, image_i3]
1389 IF (target_atom == source_atom .AND. &
1390 all(image_shift == 0)) cycle
1391 image_translation = matmul( &
1392 cell%hmat, real(image_shift,
dp))
1393 cross_displacement = composite_grid_coords(:, composite_row) - &
1394 particle_set(source_atom)%r - &
1396 CALL interpolate_gapw_atom_grid_fields( &
1397 grid_atom, harmonics, cross_displacement, cross_cutoff, nspins, &
1398 rho_h, rho_s, drho_h, drho_s, tau_h, tau_s, &
1399 cross_density, cross_grad, cross_kin, cross_density_spatial, &
1400 cross_grad_spatial, cross_kin_spatial)
1402 IF (use_atom_composite_density)
THEN
1403 composite_cross_density(composite_row, 1:2) = &
1404 composite_cross_density(composite_row, 1:2) + &
1407 IF (use_atom_composite_gradient)
THEN
1408 composite_cross_grad(composite_row, :, 1:2) = &
1409 composite_cross_grad(composite_row, :, 1:2) + &
1412 IF (use_atom_composite_tau)
THEN
1413 composite_cross_kin(composite_row, 1:2) = &
1414 composite_cross_kin(composite_row, 1:2) + cross_kin(1:2)
1417 IF (use_atom_composite_density)
THEN
1418 composite_cross_density(composite_row, :) = &
1419 composite_cross_density(composite_row, :) + &
1420 0.5_dp*cross_density(1)
1422 IF (use_atom_composite_gradient)
THEN
1424 composite_cross_grad(composite_row, idir, :) = &
1425 composite_cross_grad(composite_row, idir, :) + &
1426 0.5_dp*cross_grad(idir, 1)
1429 IF (use_atom_composite_tau)
THEN
1430 composite_cross_kin(composite_row, :) = &
1431 composite_cross_kin(composite_row, :) + &
1442 CALL release_tau_basis_cache(tau_basis_cache)
1446 IF (native_grid_diagnostics)
THEN
1447 composite_cross_density_max = maxval(abs(composite_cross_density))
1448 composite_cross_grad_max = maxval(abs(composite_cross_grad))
1449 composite_cross_kin_max = maxval(abs(composite_cross_kin))
1450 CALL para_env%max(composite_cross_cutoff_max)
1451 CALL para_env%max(composite_cross_density_max)
1452 CALL para_env%max(composite_cross_grad_max)
1453 CALL para_env%max(composite_cross_kin_max)
1456 WRITE (unit=iw, fmt=
"(T2,A,4(1X,ES20.12))") &
1457 "SKALA_GPW| Atom-composite cross support/maxima", &
1458 composite_cross_cutoff_max, composite_cross_density_max, &
1459 composite_cross_grad_max, composite_cross_kin_max
1462 composite_density(:, :) = composite_density(:, :) + composite_cross_density(:, :)
1463 composite_grad(:, :, :) = composite_grad(:, :, :) + composite_cross_grad(:, :, :)
1464 composite_kin(:, :) = composite_kin(:, :) + composite_cross_kin(:, :)
1465 DEALLOCATE (composite_cross_density, composite_cross_grad, composite_cross_kin)
1467 cpassert(all(composite_grid_weights >= 0.0_dp))
1468 atom_composite_nelec = sum(composite_grid_weights* &
1469 (composite_density(:, 1) + composite_density(:, 2)))
1470 CALL para_env%sum(atom_composite_nelec)
1471 IF (atom_composite_reference .AND. my_calculate_forces)
THEN
1473 my_xc_section, para_env, composite_density, composite_grad, composite_kin, &
1474 composite_grid_coords, composite_grid_weights, composite_atomic_grid_weights, &
1475 composite_local_grid_sizes, composite_local_atom_coords, atom_composite_exc, &
1476 composite_density_grad, composite_grad_grad, composite_kin_grad, &
1477 composite_grid_coord_grad, composite_grid_weight_grad, &
1478 composite_atomic_grid_weight_grad, composite_atom_coord_grad)
1479 ALLOCATE (composite_cross_force(3,
SIZE(particle_set)), &
1480 composite_explicit_force(3,
SIZE(particle_set)), &
1481 composite_grid_coord_force(3,
SIZE(particle_set)), &
1482 composite_model_atom_force(3,
SIZE(particle_set)), &
1483 composite_moving_smooth_force(3,
SIZE(particle_set)), &
1484 composite_nlcc_center_force(3,
SIZE(particle_set)), &
1485 composite_nlcc_target_force(3,
SIZE(particle_set)), &
1486 composite_partition_force(3,
SIZE(particle_set)), &
1487 composite_partition_included(
SIZE(particle_set)), &
1488 composite_partition_datom(3,
SIZE(particle_set),
SIZE(particle_set)), &
1489 composite_partition_dstrain(3, 3,
SIZE(particle_set)))
1490 composite_cross_force = 0.0_dp
1491 composite_cross_image_virial = 0.0_dp
1492 composite_model_atom_force = 0.0_dp
1493 composite_grid_coord_force = 0.0_dp
1494 composite_moving_smooth_force = 0.0_dp
1495 composite_nlcc_center_force = 0.0_dp
1496 composite_nlcc_target_force = 0.0_dp
1497 composite_partition_force = 0.0_dp
1498 composite_explicit_virial = 0.0_dp
1499 composite_feature_virial = 0.0_dp
1500 composite_interpolation_virial = 0.0_dp
1501 composite_partition_strain_virial = 0.0_dp
1510 DO composite_local_atom = 1, composite_local_natom
1511 iatom = composite_local_atoms(composite_local_atom)
1512 DO composite_row = composite_atom_start(iatom), composite_atom_end(iatom)
1515 composite_smooth_gradient_value(jdir, 1) = &
1516 interpolate_native_grid( &
1517 composite_smooth_drhoa(jdir)%array, smooth_rho_r(1)%pw_grid, &
1518 cell, composite_grid_coords(:, composite_row), &
1519 image_partition_atom_composite)
1520 composite_smooth_gradient_value(jdir, 2) = &
1521 interpolate_native_grid( &
1522 composite_smooth_drhob(jdir)%array, smooth_rho_r(1)%pw_grid, &
1523 cell, composite_grid_coords(:, composite_row), &
1524 image_partition_atom_composite)
1528 composite_smooth_gradient_value(jdir, :) = 0.5_dp* &
1529 interpolate_native_grid( &
1530 composite_smooth_drho(jdir)%array, smooth_rho_r(1)%pw_grid, &
1531 cell, composite_grid_coords(:, composite_row), &
1532 image_partition_atom_composite)
1538 composite_feature_virial(jdir, idir) = &
1539 composite_feature_virial(jdir, idir) - &
1540 composite_grad_grad(composite_row, idir, ispin)* &
1541 composite_smooth_gradient_value(jdir, ispin)
1572 ALLOCATE (composite_nlcc_center_force_local(3,
SIZE(particle_set)), &
1573 composite_partition_force_local(3,
SIZE(particle_set)), &
1574 local_partition_datom(3,
SIZE(particle_set)))
1575 composite_nlcc_center_force_local = 0.0_dp
1576 composite_partition_force_local = 0.0_dp
1578 DO composite_local_atom = 1, composite_local_natom
1579 iatom = composite_local_atoms(composite_local_atom)
1580 composite_model_atom_force(:, iatom) = &
1581 composite_atom_coord_grad(:, composite_local_atom)
1582 DO composite_row = composite_atom_start(iatom), composite_atom_end(iatom)
1583 composite_grid_coord_force(:, iatom) = &
1584 composite_grid_coord_force(:, iatom) + &
1585 composite_grid_coord_grad(:, composite_row)
1587 spatial_derivative = &
1588 composite_density_grad(composite_row, 1)* &
1589 interpolate_native_grid_gradient( &
1590 composite_smooth_rhoa, smooth_rho_r(1)%pw_grid, cell, &
1591 composite_grid_coords(:, composite_row), &
1592 image_partition_atom_composite) + &
1593 composite_density_grad(composite_row, 2)* &
1594 interpolate_native_grid_gradient( &
1595 composite_smooth_rhob, smooth_rho_r(1)%pw_grid, cell, &
1596 composite_grid_coords(:, composite_row), &
1597 image_partition_atom_composite) + &
1598 composite_kin_grad(composite_row, 1)* &
1599 interpolate_native_grid_gradient( &
1600 composite_smooth_tau_a, smooth_rho_r(1)%pw_grid, cell, &
1601 composite_grid_coords(:, composite_row), &
1602 image_partition_atom_composite) + &
1603 composite_kin_grad(composite_row, 2)* &
1604 interpolate_native_grid_gradient( &
1605 composite_smooth_tau_b, smooth_rho_r(1)%pw_grid, cell, &
1606 composite_grid_coords(:, composite_row), &
1607 image_partition_atom_composite)
1609 spatial_derivative = spatial_derivative + &
1610 composite_grad_grad(composite_row, idir, 1)* &
1611 interpolate_native_grid_gradient( &
1612 composite_smooth_drhoa(idir)%array, &
1613 smooth_rho_r(1)%pw_grid, cell, &
1614 composite_grid_coords(:, composite_row), &
1615 image_partition_atom_composite) + &
1616 composite_grad_grad(composite_row, idir, 2)* &
1617 interpolate_native_grid_gradient( &
1618 composite_smooth_drhob(idir)%array, &
1619 smooth_rho_r(1)%pw_grid, cell, &
1620 composite_grid_coords(:, composite_row), &
1621 image_partition_atom_composite)
1624 spatial_derivative = 0.5_dp*sum( &
1625 composite_density_grad(composite_row, :))* &
1626 interpolate_native_grid_gradient( &
1627 composite_smooth_rho, smooth_rho_r(1)%pw_grid, cell, &
1628 composite_grid_coords(:, composite_row), &
1629 image_partition_atom_composite) + &
1630 0.5_dp*sum(composite_kin_grad(composite_row, :))* &
1631 interpolate_native_grid_gradient( &
1632 composite_smooth_tau, smooth_rho_r(1)%pw_grid, cell, &
1633 composite_grid_coords(:, composite_row), &
1634 image_partition_atom_composite)
1636 spatial_derivative = spatial_derivative + 0.5_dp*sum( &
1637 composite_grad_grad(composite_row, idir, :))* &
1638 interpolate_native_grid_gradient( &
1639 composite_smooth_drho(idir)%array, &
1640 smooth_rho_r(1)%pw_grid, &
1641 cell, composite_grid_coords(:, composite_row), &
1642 image_partition_atom_composite)
1647 composite_interpolation_virial(idir, jdir) = &
1648 composite_interpolation_virial(idir, jdir) + &
1649 spatial_derivative(idir)*( &
1650 composite_grid_coords(jdir, composite_row) - &
1651 particle_set(iatom)%r(jdir))
1654 composite_moving_smooth_force(:, iatom) = &
1655 composite_moving_smooth_force(:, iatom) + spatial_derivative
1657 nlcc_spin_factor = merge(1.0_dp, 0.5_dp, lsd)
1658 DO source_atom = 1,
SIZE(particle_set)
1659 NULLIFY (gth_potential, sgp_potential)
1660 CALL get_qs_kind(my_kind_set(composite_atom_kind(source_atom)), &
1661 gth_potential=gth_potential, &
1662 sgp_potential=sgp_potential)
1664 composite_grid_coords(:, composite_row), &
1665 particle_set(source_atom)%r, gth_potential, sgp_potential, &
1666 nlcc_density, nlcc_gradient, nlcc_hessian)
1667 nlcc_spatial_derivative = 0.0_dp
1669 nlcc_spatial_derivative = nlcc_spatial_derivative + &
1670 nlcc_spin_factor*composite_density_grad(composite_row, ispin)* &
1674 nlcc_spatial_derivative(jdir) = &
1675 nlcc_spatial_derivative(jdir) + nlcc_spin_factor* &
1676 composite_grad_grad(composite_row, idir, ispin)* &
1677 nlcc_hessian(idir, jdir)
1681 composite_nlcc_target_force(:, iatom) = &
1682 composite_nlcc_target_force(:, iatom) + nlcc_spatial_derivative
1683 composite_nlcc_center_force_local(:, source_atom) = &
1684 composite_nlcc_center_force_local(:, source_atom) - nlcc_spatial_derivative
1687 IF (image_partition_atom_composite)
THEN
1689 composite_grid_coords(:, composite_row), composite_atom_coords, cell, &
1690 iatom, local_partition_weight, local_partition_datom, &
1691 local_partition_dstrain, composite_image_periodicity)
1693 composite_grid_coords(:, composite_row), &
1694 composite_atom_coords(:, iatom:iatom), cell, 1, &
1695 descriptor_window_weight, local_descriptor_datom, &
1696 local_descriptor_dstrain, composite_image_periodicity)
1699 target_partition_adjoint = composite_base_grid_weights(composite_row)* &
1700 composite_grid_weight_grad(composite_row)
1701 descriptor_window_adjoint = composite_base_grid_weights(composite_row)* &
1702 composite_atomic_grid_weight_grad(composite_row)* &
1704 descriptor_window_weight)
1705 DO target_atom = 1,
SIZE(particle_set)
1706 composite_partition_force_local(:, target_atom) = &
1707 composite_partition_force_local(:, target_atom) + &
1708 target_partition_adjoint*local_partition_datom(:, target_atom)
1710 composite_partition_force_local(:, iatom) = &
1711 composite_partition_force_local(:, iatom) - target_partition_adjoint* &
1712 sum(local_partition_datom, dim=2)
1713 composite_partition_strain_virial = &
1714 composite_partition_strain_virial - target_partition_adjoint* &
1715 local_partition_dstrain - descriptor_window_adjoint* &
1716 local_descriptor_dstrain
1719 composite_grid_coords(:, composite_row), composite_atom_coords, cell, &
1720 composite_partition_weights, composite_partition_included, &
1721 composite_partition_datom, composite_partition_dstrain)
1722 partition_adjoint = composite_grid_weight_grad(composite_row)* &
1723 composite_atomic_grid_weights(composite_row)
1724 DO target_atom = 1,
SIZE(particle_set)
1725 composite_partition_force_local(:, target_atom) = &
1726 composite_partition_force_local(:, target_atom) + &
1727 partition_adjoint* &
1728 composite_partition_datom(:, target_atom, iatom)
1730 composite_partition_force_local(:, iatom) = &
1731 composite_partition_force_local(:, iatom) - partition_adjoint* &
1732 sum(composite_partition_datom(:, :, iatom), dim=2)
1738 composite_nlcc_center_force(:, :) = composite_nlcc_center_force(:, :) + &
1739 composite_nlcc_center_force_local
1740 composite_partition_force(:, :) = composite_partition_force(:, :) + &
1741 composite_partition_force_local
1743 DEALLOCATE (composite_nlcc_center_force_local, composite_partition_force_local, &
1744 local_partition_datom)
1749 IF (.NOT. direct_valence_atom_composite)
THEN
1750 DO ikind = 1,
SIZE(atomic_kind_set)
1751 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1752 NULLIFY (gth_potential, sgp_potential)
1753 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
1754 gth_potential=gth_potential, harmonics=harmonics, &
1755 grid_atom=grid_atom, sgp_potential=sgp_potential, &
1756 zatom=zatom, zeff=zeff)
1757 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type=
"GAPW_1C")
1758 IF (.NOT. native_skala_uses_one_center_kind( &
1759 paw_atom, gapw_representation, &
1760 ASSOCIATED(gth_potential) .OR.
ASSOCIATED(sgp_potential), &
1764 na = grid_atom%ng_sphere
1765 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
1766 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
1767 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
1768 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
1769 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
1770 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
1771 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
1774 source_atom = atom_list(iat)
1775 rho_atom => my_rho_atom_set(source_atom)
1776 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1777 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s, &
1778 drho_rad_h=dr_h, drho_rad_s=dr_s, &
1779 rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1784 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
1787 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
1788 r_h_d, r_s_d, drho_h, drho_s)
1791 cross_cutoff = gapw_atom_grid_support_radius( &
1792 grid_atom, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s)
1793 IF (cross_cutoff <= 0.0_dp) cycle
1796 IF (cell%perd(idir) == 1)
THEN
1797 image_shell(idir) = ceiling( &
1798 cross_cutoff*sqrt(sum(cell%h_inv(idir, :)**2))) + 1
1815 ALLOCATE (composite_cross_force_local(3,
SIZE(particle_set)))
1816 composite_cross_force_local = 0.0_dp
1817 composite_cross_image_virial_local = 0.0_dp
1819 DO composite_row = 1, composite_nflat
1820 target_atom = composite_grid_atom(composite_row)
1822 cross_density_adjoint = 0.0_dp
1823 cross_grad_adjoint = 0.0_dp
1824 cross_kin_adjoint = 0.0_dp
1825 IF (use_atom_composite_density)
THEN
1826 cross_density_adjoint(1:2) = &
1827 composite_density_grad(composite_row, 1:2)
1829 IF (use_atom_composite_gradient)
THEN
1830 cross_grad_adjoint(:, 1:2) = &
1831 composite_grad_grad(composite_row, :, 1:2)
1833 IF (use_atom_composite_tau)
THEN
1834 cross_kin_adjoint(1:2) = &
1835 composite_kin_grad(composite_row, 1:2)
1838 cross_density_adjoint = 0.0_dp
1839 cross_grad_adjoint = 0.0_dp
1840 cross_kin_adjoint = 0.0_dp
1841 IF (use_atom_composite_density)
THEN
1842 cross_density_adjoint(1) = 0.5_dp* &
1843 sum(composite_density_grad(composite_row, :))
1845 IF (use_atom_composite_gradient)
THEN
1847 cross_grad_adjoint(idir, 1) = 0.5_dp* &
1848 sum(composite_grad_grad(composite_row, idir, :))
1851 IF (use_atom_composite_tau)
THEN
1852 cross_kin_adjoint(1) = 0.5_dp* &
1853 sum(composite_kin_grad(composite_row, :))
1860 fractional(idir) = fractional(idir) + cell%h_inv(idir, jdir)* &
1861 (composite_grid_coords(jdir, composite_row) - &
1862 particle_set(source_atom)%r(jdir))
1866 base_shift(idir) = cell%perd(idir)*nint(fractional(idir))
1868 DO image_i3 = base_shift(3) - image_shell(3), &
1869 base_shift(3) + image_shell(3)
1870 DO image_i2 = base_shift(2) - image_shell(2), &
1871 base_shift(2) + image_shell(2)
1872 DO image_i1 = base_shift(1) - image_shell(1), &
1873 base_shift(1) + image_shell(1)
1874 image_shift = [image_i1, image_i2, image_i3]
1875 IF (target_atom == source_atom .AND. &
1876 all(image_shift == 0)) cycle
1877 image_translation = matmul( &
1878 cell%hmat, real(image_shift,
dp))
1879 cross_displacement = &
1880 composite_grid_coords(:, composite_row) - &
1881 particle_set(source_atom)%r - image_translation
1882 CALL interpolate_gapw_atom_grid_fields( &
1883 grid_atom, harmonics, cross_displacement, cross_cutoff, &
1884 nspins, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s, &
1885 cross_density, cross_grad, cross_kin, cross_density_spatial, &
1886 cross_grad_spatial, cross_kin_spatial)
1887 cross_spatial_derivative = 0.0_dp
1888 DO ispin = 1, nspins
1890 cross_spatial_derivative(idir) = &
1891 cross_spatial_derivative(idir) + &
1892 cross_density_adjoint(ispin)* &
1893 cross_density_spatial(idir, ispin) + &
1894 cross_kin_adjoint(ispin)* &
1895 cross_kin_spatial(idir, ispin)
1897 cross_spatial_derivative(idir) = &
1898 cross_spatial_derivative(idir) + &
1899 cross_grad_adjoint(jdir, ispin)* &
1900 cross_grad_spatial(jdir, idir, ispin)
1904 composite_cross_force_local(:, target_atom) = &
1905 composite_cross_force_local(:, target_atom) + &
1906 cross_spatial_derivative
1907 composite_cross_force_local(:, source_atom) = &
1908 composite_cross_force_local(:, source_atom) - &
1909 cross_spatial_derivative
1912 composite_cross_image_virial_local(idir, jdir) = &
1913 composite_cross_image_virial_local(idir, jdir) + &
1914 cross_spatial_derivative(idir)*image_translation(jdir)
1923 composite_cross_force(:, :) = &
1924 composite_cross_force(:, :) + composite_cross_force_local(:, :)
1925 composite_cross_image_virial = composite_cross_image_virial + &
1926 composite_cross_image_virial_local
1928 DEALLOCATE (composite_cross_force_local)
1931 CALL release_tau_basis_cache(tau_basis_cache)
1935 composite_explicit_force(:, :) = composite_model_atom_force(:, :) + &
1936 composite_grid_coord_force(:, :) + &
1937 composite_moving_smooth_force(:, :) + &
1938 composite_cross_force(:, :) + &
1939 composite_nlcc_center_force(:, :) + &
1940 composite_nlcc_target_force(:, :) + &
1941 composite_partition_force(:, :)
1948 DO iatom = 1,
SIZE(particle_set)
1951 composite_explicit_virial(idir, jdir) = &
1952 composite_explicit_virial(idir, jdir) - ( &
1953 composite_model_atom_force(idir, iatom) + &
1954 composite_grid_coord_force(idir, iatom) + &
1955 composite_cross_force(idir, iatom) + &
1956 composite_nlcc_center_force(idir, iatom) + &
1957 composite_nlcc_target_force(idir, iatom) + &
1958 composite_partition_force(idir, iatom))*particle_set(iatom)%r(jdir)
1962 composite_explicit_virial = composite_explicit_virial + &
1963 composite_partition_strain_virial + &
1964 composite_cross_image_virial
1965 IF (
ASSOCIATED(force))
THEN
1966 DO iatom = 1,
SIZE(particle_set)
1967 ikind = composite_atom_kind(iatom)
1968 iat = composite_atom_kind_index(iatom)
1969 cpassert(ikind > 0 .AND. iat > 0)
1970 force(ikind)%rho_elec(:, iat) = force(ikind)%rho_elec(:, iat) + &
1971 composite_explicit_force(:, iatom)
1974 IF (use_virial)
THEN
1979 virial%pv_xc = composite_feature_virial - composite_interpolation_virial
1980 virial%pv_gapw = virial%pv_gapw + composite_explicit_virial
1981 virial%pv_virial = virial%pv_virial + composite_explicit_virial
1983 IF (native_grid_diagnostics)
THEN
1984 CALL para_env%sum(composite_cross_force)
1985 CALL para_env%sum(composite_cross_image_virial)
1986 CALL para_env%sum(composite_model_atom_force)
1987 CALL para_env%sum(composite_grid_coord_force)
1988 CALL para_env%sum(composite_moving_smooth_force)
1989 CALL para_env%sum(composite_nlcc_center_force)
1990 CALL para_env%sum(composite_nlcc_target_force)
1991 CALL para_env%sum(composite_partition_force)
1992 CALL para_env%sum(composite_explicit_force)
1993 CALL para_env%sum(composite_explicit_virial)
1994 CALL para_env%sum(composite_feature_virial)
1995 CALL para_env%sum(composite_interpolation_virial)
1998 DO iatom = 1,
SIZE(particle_set)
1999 WRITE (unit=iw, fmt=
"(T2,A,1X,I0,3(1X,ES20.12))") &
2000 "SKALA_GPW| Atom-composite model-atom force", iatom, &
2001 composite_model_atom_force(:, iatom)
2002 WRITE (unit=iw, fmt=
"(T2,A,1X,I0,3(1X,ES20.12))") &
2003 "SKALA_GPW| Atom-composite grid-coordinate force", iatom, &
2004 composite_grid_coord_force(:, iatom)
2005 WRITE (unit=iw, fmt=
"(T2,A,1X,I0,3(1X,ES20.12))") &
2006 "SKALA_GPW| Atom-composite moving-smooth force", iatom, &
2007 composite_moving_smooth_force(:, iatom)
2008 WRITE (unit=iw, fmt=
"(T2,A,1X,I0,3(1X,ES20.12))") &
2009 "SKALA_GPW| Atom-composite cross-region force", iatom, &
2010 composite_cross_force(:, iatom)
2011 WRITE (unit=iw, fmt=
"(T2,A,1X,I0,3(1X,ES20.12))") &
2012 "SKALA_GPW| Atom-composite NLCC-center force", iatom, &
2013 composite_nlcc_center_force(:, iatom)
2014 WRITE (unit=iw, fmt=
"(T2,A,1X,I0,3(1X,ES20.12))") &
2015 "SKALA_GPW| Atom-composite NLCC-target force", iatom, &
2016 composite_nlcc_target_force(:, iatom)
2017 WRITE (unit=iw, fmt=
"(T2,A,1X,I0,3(1X,ES20.12))") &
2018 "SKALA_GPW| Atom-composite partition force", iatom, &
2019 composite_partition_force(:, iatom)
2020 WRITE (unit=iw, fmt=
"(T2,A,1X,I0,3(1X,ES20.12))") &
2021 "SKALA_GPW| Atom-composite explicit force", iatom, &
2022 composite_explicit_force(:, iatom)
2024 WRITE (unit=iw, fmt=
"(T2,A)") &
2025 "SKALA_GPW| Atom-composite explicit virial"
2027 WRITE (unit=iw, fmt=
"(T2,A,1X,3ES20.10)") &
2028 "SKALA_GPW|", composite_explicit_virial(idir, :)
2030 WRITE (unit=iw, fmt=
"(T2,A)") &
2031 "SKALA_GPW| Atom-composite cross-image virial"
2033 WRITE (unit=iw, fmt=
"(T2,A,1X,3ES20.10)") &
2034 "SKALA_GPW|", composite_cross_image_virial(idir, :)
2036 WRITE (unit=iw, fmt=
"(T2,A)") &
2037 "SKALA_GPW| Atom-composite feature virial"
2039 WRITE (unit=iw, fmt=
"(T2,A,1X,3ES20.10)") &
2040 "SKALA_GPW|", composite_feature_virial(idir, :)
2042 WRITE (unit=iw, fmt=
"(T2,A)") &
2043 "SKALA_GPW| Atom-composite interpolation virial"
2045 WRITE (unit=iw, fmt=
"(T2,A,1X,3ES20.10)") &
2046 "SKALA_GPW|", composite_interpolation_virial(idir, :)
2050 DEALLOCATE (composite_atom_coord_grad, composite_atomic_grid_weight_grad, &
2051 composite_cross_force, &
2052 composite_explicit_force, composite_grid_coord_force, &
2053 composite_grid_coord_grad, composite_grid_weight_grad, &
2054 composite_model_atom_force, composite_moving_smooth_force, &
2055 composite_nlcc_center_force, &
2056 composite_nlcc_target_force, &
2057 composite_partition_datom, composite_partition_dstrain, &
2058 composite_partition_force, composite_partition_included)
2061 my_xc_section, para_env, composite_density, composite_grad, composite_kin, &
2062 composite_grid_coords, composite_grid_weights, composite_atomic_grid_weights, &
2063 composite_local_grid_sizes, composite_local_atom_coords, atom_composite_exc, &
2064 composite_density_grad, composite_grad_grad, composite_kin_grad)
2066 composite_pw_nflat = product(smooth_rho_r(1)%pw_grid%npts)
2067 adjoint_nchannels = merge(2, 1, lsd)
2068 ALLOCATE (smooth_density_adjoint_storage(composite_pw_nflat, adjoint_nchannels), &
2069 smooth_grad_adjoint_storage(composite_pw_nflat, 3, adjoint_nchannels), &
2070 smooth_kin_adjoint_storage(composite_pw_nflat, adjoint_nchannels))
2071 smooth_density_adjoint_storage = 0.0_dp
2072 smooth_grad_adjoint_storage = 0.0_dp
2073 smooth_kin_adjoint_storage = 0.0_dp
2074 CALL build_native_grid_adjoint_bins( &
2075 smooth_rho_r(1)%pw_grid, cell, composite_grid_coords, &
2076 image_partition_atom_composite, adjoint_tile_count, adjoint_bin_offsets, &
2078 adjoint_nbins =
SIZE(adjoint_bin_offsets) - 1
2088 DO adjoint_bin = 1, adjoint_nbins
2089 CALL native_grid_adjoint_tile_bounds( &
2090 smooth_rho_r(1)%pw_grid, adjoint_bin, adjoint_tile_count, &
2091 adjoint_tile_lower, adjoint_tile_upper)
2092 DO adjoint_entry = adjoint_bin_offsets(adjoint_bin), &
2093 adjoint_bin_offsets(adjoint_bin + 1) - 1
2094 composite_row = adjoint_bin_rows(adjoint_entry)
2095 CALL create_native_grid_interpolation_stencil( &
2096 interpolation_stencil, smooth_rho_r(1)%pw_grid, cell, &
2097 composite_grid_coords(:, composite_row), image_partition_atom_composite)
2099 composite_smooth_density_adjoint_value = &
2100 composite_density_grad(composite_row, :)
2101 composite_smooth_gradient_adjoint_value = &
2102 composite_grad_grad(composite_row, :, :)
2103 composite_smooth_kin_adjoint_value = composite_kin_grad(composite_row, :)
2105 composite_smooth_density_adjoint_value = 0.0_dp
2106 composite_smooth_gradient_adjoint_value = 0.0_dp
2107 composite_smooth_kin_adjoint_value = 0.0_dp
2108 composite_smooth_density_adjoint_value(1) = &
2109 sum(composite_density_grad(composite_row, :))
2110 composite_smooth_gradient_adjoint_value(:, 1) = &
2111 sum(composite_grad_grad(composite_row, :, :), dim=2)
2112 composite_smooth_kin_adjoint_value(1) = &
2113 sum(composite_kin_grad(composite_row, :))
2115 CALL add_native_grid_fields_adjoint_tile( &
2116 smooth_density_adjoint_storage, smooth_grad_adjoint_storage, &
2117 smooth_kin_adjoint_storage, smooth_rho_r(1)%pw_grid, interpolation_stencil, &
2118 composite_smooth_density_adjoint_value, &
2119 composite_smooth_gradient_adjoint_value, &
2120 composite_smooth_kin_adjoint_value, adjoint_nchannels, &
2121 adjoint_tile_lower, adjoint_tile_upper)
2125 DEALLOCATE (adjoint_bin_offsets, adjoint_bin_rows)
2126 smooth_input_contraction = 0.0_dp
2133 DO composite_row = 1, composite_nflat
2134 composite_smooth_density_value = composite_smooth_density_cache(composite_row, :)
2135 composite_smooth_gradient_value = &
2136 composite_smooth_gradient_cache(composite_row, :, :)
2137 composite_smooth_kin_value = composite_smooth_kin_cache(composite_row, :)
2140 smooth_input_contraction = smooth_input_contraction + &
2141 composite_density_grad(composite_row, ispin)* &
2142 composite_smooth_density_value(ispin) + &
2143 composite_kin_grad(composite_row, ispin)* &
2144 composite_smooth_kin_value(ispin)
2146 smooth_input_contraction = smooth_input_contraction + &
2147 composite_grad_grad(composite_row, idir, ispin)* &
2148 composite_smooth_gradient_value(idir, ispin)
2151 smooth_input_contraction = smooth_input_contraction + 0.5_dp*( &
2152 composite_density_grad(composite_row, ispin)* &
2153 composite_smooth_density_value(1) + &
2154 composite_kin_grad(composite_row, ispin)* &
2155 composite_smooth_kin_value(1))
2157 smooth_input_contraction = smooth_input_contraction + 0.5_dp* &
2158 composite_grad_grad(composite_row, idir, ispin)* &
2159 composite_smooth_gradient_value(idir, 1)
2165 smooth_density_adjoint => smooth_density_adjoint_storage
2166 smooth_grad_adjoint => smooth_grad_adjoint_storage
2167 smooth_kin_adjoint => smooth_kin_adjoint_storage
2170 CALL para_env%sum(smooth_density_adjoint)
2171 CALL para_env%sum(smooth_grad_adjoint)
2172 CALL para_env%sum(smooth_kin_adjoint)
2173 CALL para_env%sum(smooth_input_contraction)
2175 smooth_vxc_rho, smooth_vxc_tau, smooth_rho_r, auxbas_pw_pool, &
2176 smooth_density_adjoint, smooth_grad_adjoint, smooth_kin_adjoint, &
2177 xc_deriv_method_id, global_grid_layout=.true.)
2178 NULLIFY (smooth_density_adjoint, smooth_grad_adjoint, smooth_kin_adjoint)
2179 DEALLOCATE (smooth_density_adjoint_storage, smooth_grad_adjoint_storage, &
2180 smooth_kin_adjoint_storage)
2181 smooth_grid_contraction = 0.0_dp
2182 DO ispin = 1, nspins
2183 smooth_grid_contraction = smooth_grid_contraction + smooth_rho_r(1)%pw_grid%dvol* &
2184 (sum(smooth_vxc_rho(ispin)%array* &
2185 smooth_rho_r(ispin)%array) + &
2186 sum(smooth_vxc_tau(ispin)%array* &
2187 smooth_tau_r(ispin)%array))
2188 IF (atom_composite_reference)
THEN
2189 cpassert(
PRESENT(composite_vxc_rho))
2190 cpassert(
PRESENT(composite_vxc_tau))
2191 cpassert(
ASSOCIATED(composite_vxc_rho))
2192 cpassert(
ASSOCIATED(composite_vxc_tau))
2193 cpassert(
SIZE(composite_vxc_rho) == nspins)
2194 cpassert(
SIZE(composite_vxc_tau) == nspins)
2195 CALL pw_axpy(smooth_vxc_rho(ispin), composite_vxc_rho(ispin), 1.0_dp)
2196 CALL pw_axpy(smooth_vxc_tau(ispin), composite_vxc_tau(ispin), 1.0_dp)
2198 CALL auxbas_pw_pool%give_back_pw(smooth_vxc_rho(ispin))
2199 CALL auxbas_pw_pool%give_back_pw(smooth_vxc_tau(ispin))
2201 CALL para_env%sum(smooth_grid_contraction)
2202 DEALLOCATE (smooth_vxc_rho, smooth_vxc_tau)
2204 one_center_field_contraction = 0.0_dp
2205 one_center_matrix_contraction = 0.0_dp
2206 one_center_density_field_contraction = 0.0_dp
2207 one_center_density_matrix_contraction = 0.0_dp
2208 one_center_gradient_field_contraction = 0.0_dp
2209 one_center_gradient_matrix_contraction = 0.0_dp
2210 one_center_rho_grad_field_contraction = 0.0_dp
2211 one_center_rho_grad_matrix_contraction = 0.0_dp
2212 one_center_tau_field_contraction = 0.0_dp
2213 one_center_tau_matrix_contraction = 0.0_dp
2214 IF (.NOT. direct_valence_atom_composite)
THEN
2215 DO ikind = 1,
SIZE(atomic_kind_set)
2216 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
2217 NULLIFY (gth_potential, sgp_potential)
2218 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
2219 gth_potential=gth_potential, harmonics=harmonics, &
2220 grid_atom=grid_atom, sgp_potential=sgp_potential, &
2221 zatom=zatom, zeff=zeff)
2222 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type=
"GAPW_1C")
2223 IF (.NOT. native_skala_uses_one_center_kind( &
2224 paw_atom, gapw_representation, &
2225 ASSOCIATED(gth_potential) .OR.
ASSOCIATED(sgp_potential), &
2229 na = grid_atom%ng_sphere
2230 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
2231 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
2232 CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
2233 CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
2234 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
2235 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
2236 CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
2237 CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
2238 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
2239 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
2240 CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
2241 CALL reallocate(vtau_s, 1, na, 1, nr, 1, nspins)
2242 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
2244 bo =
get_limit(natom, para_env%num_pe, para_env%mepos)
2246 iatom = atom_list(iat)
2247 source_matrix_local = iat >= bo(1) .AND. iat <= bo(2)
2248 rho_atom => my_rho_atom_set(iatom)
2249 NULLIFY (cpc_h, cpc_s, r_h, r_s, dr_h, dr_s, r_h_d, r_s_d, int_hh, int_ss)
2250 CALL get_rho_atom(rho_atom=rho_atom, cpc_h=cpc_h, cpc_s=cpc_s, &
2251 rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
2252 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d, &
2253 ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
2258 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
2261 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
2262 r_h_d, r_s_d, drho_h, drho_s)
2271 IF (source_matrix_local)
THEN
2272 composite_row = composite_atom_start(iatom) - 1
2275 composite_row = composite_row + 1
2277 IF (use_atom_composite_density)
THEN
2278 vxc_h(ia, ir, 1:2) = composite_density_grad(composite_row, 1:2)
2279 vxc_s(ia, ir, 1:2) = composite_density_grad(composite_row, 1:2)
2281 IF (use_atom_composite_gradient)
THEN
2283 vxg_h(idir, ia, ir, 1:2) = &
2284 composite_grad_grad(composite_row, idir, 1:2)
2285 vxg_s(idir, ia, ir, 1:2) = &
2286 composite_grad_grad(composite_row, idir, 1:2)
2289 IF (use_atom_composite_tau)
THEN
2290 vtau_h(ia, ir, 1:2) = composite_kin_grad(composite_row, 1:2)
2291 vtau_s(ia, ir, 1:2) = composite_kin_grad(composite_row, 1:2)
2294 IF (use_atom_composite_density)
THEN
2295 vxc_h(ia, ir, 1) = 0.5_dp* &
2296 sum(composite_density_grad(composite_row, :))
2297 vxc_s(ia, ir, 1) = vxc_h(ia, ir, 1)
2299 IF (use_atom_composite_gradient)
THEN
2301 vxg_h(idir, ia, ir, 1) = 0.5_dp* &
2302 sum(composite_grad_grad( &
2303 composite_row, idir, :))
2304 vxg_s(idir, ia, ir, 1) = vxg_h(idir, ia, ir, 1)
2307 IF (use_atom_composite_tau)
THEN
2308 vtau_h(ia, ir, 1) = 0.5_dp* &
2309 sum(composite_kin_grad(composite_row, :))
2310 vtau_s(ia, ir, 1) = vtau_h(ia, ir, 1)
2315 cpassert(composite_row == composite_atom_end(iatom))
2318 cross_cutoff = gapw_atom_grid_support_radius( &
2319 grid_atom, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s)
2320 IF (cross_cutoff > 0.0_dp)
THEN
2323 IF (cell%perd(idir) == 1)
THEN
2324 image_shell(idir) = ceiling( &
2325 cross_cutoff*sqrt(sum(cell%h_inv(idir, :)**2))) + 1
2328 DO composite_local_atom = 1, composite_local_natom
2329 target_atom = composite_local_atoms(composite_local_atom)
2330 DO composite_row = composite_atom_start(target_atom), &
2331 composite_atom_end(target_atom)
2332 cross_density_adjoint = 0.0_dp
2333 cross_grad_adjoint = 0.0_dp
2334 cross_kin_adjoint = 0.0_dp
2336 IF (use_atom_composite_density)
THEN
2337 cross_density_adjoint(1:2) = &
2338 composite_density_grad(composite_row, 1:2)
2340 IF (use_atom_composite_gradient)
THEN
2341 cross_grad_adjoint(:, 1:2) = &
2342 composite_grad_grad(composite_row, :, 1:2)
2344 IF (use_atom_composite_tau)
THEN
2345 cross_kin_adjoint(1:2) = &
2346 composite_kin_grad(composite_row, 1:2)
2349 IF (use_atom_composite_density)
THEN
2350 cross_density_adjoint(1) = 0.5_dp* &
2351 sum(composite_density_grad(composite_row, :))
2353 IF (use_atom_composite_gradient)
THEN
2355 cross_grad_adjoint(idir, 1) = 0.5_dp* &
2356 sum(composite_grad_grad(composite_row, idir, :))
2359 IF (use_atom_composite_tau)
THEN
2360 cross_kin_adjoint(1) = 0.5_dp* &
2361 sum(composite_kin_grad(composite_row, :))
2367 fractional(idir) = fractional(idir) + cell%h_inv(idir, jdir)* &
2368 (composite_grid_coords(jdir, composite_row) - &
2369 particle_set(iatom)%r(jdir))
2373 base_shift(idir) = cell%perd(idir)*nint(fractional(idir))
2375 DO image_i3 = base_shift(3) - image_shell(3), &
2376 base_shift(3) + image_shell(3)
2377 DO image_i2 = base_shift(2) - image_shell(2), &
2378 base_shift(2) + image_shell(2)
2379 DO image_i1 = base_shift(1) - image_shell(1), &
2380 base_shift(1) + image_shell(1)
2381 image_shift = [image_i1, image_i2, image_i3]
2382 IF (target_atom == iatom .AND. all(image_shift == 0)) cycle
2383 image_translation = matmul( &
2384 cell%hmat, real(image_shift,
dp))
2385 cross_displacement = &
2386 composite_grid_coords(:, composite_row) - &
2387 particle_set(iatom)%r - image_translation
2388 CALL add_gapw_atom_grid_interpolation_adjoint( &
2389 grid_atom, harmonics, cross_displacement, cross_cutoff, &
2390 nspins, cross_density_adjoint, cross_grad_adjoint, &
2391 cross_kin_adjoint, vxc_h, vxc_s, vxg_h, vxg_s, &
2403 CALL para_env%sum(vxc_h)
2404 CALL para_env%sum(vxc_s)
2405 CALL para_env%sum(vxg_h)
2406 CALL para_env%sum(vxg_s)
2407 CALL para_env%sum(vtau_h)
2408 CALL para_env%sum(vtau_s)
2409 IF (.NOT. source_matrix_local) cycle
2411 one_center_rho_grad_field_contraction = &
2412 one_center_rho_grad_field_contraction + sum(vxc_h*(rho_h - rho_s)) + &
2413 sum(vxg_h*(drho_h(1:3, :, :, :) - drho_s(1:3, :, :, :)))
2414 one_center_density_field_contraction = one_center_density_field_contraction + &
2415 sum(vxc_h*(rho_h - rho_s))
2416 one_center_gradient_field_contraction = one_center_gradient_field_contraction + &
2417 sum(vxg_h*(drho_h(1:3, :, :, :) - &
2418 drho_s(1:3, :, :, :)))
2419 one_center_tau_field_contraction = one_center_tau_field_contraction + &
2420 sum(vtau_h*(tau_h - tau_s))
2422 ALLOCATE (composite_int_h(
SIZE(int_hh(1)%r_coef, 1), &
2423 SIZE(int_hh(1)%r_coef, 2), nspins), &
2424 composite_int_s(
SIZE(int_ss(1)%r_coef, 1), &
2425 SIZE(int_ss(1)%r_coef, 2), nspins))
2426 DO ispin = 1, nspins
2427 composite_int_h(:, :, ispin) = int_hh(ispin)%r_coef
2428 composite_int_s(:, :, ispin) = int_ss(ispin)%r_coef
2429 int_hh(ispin)%r_coef = 0.0_dp
2430 int_ss(ispin)%r_coef = 0.0_dp
2433 grid_atom, basis_1c, harmonics, nspins)
2434 DO ispin = 1, nspins
2435 one_center_density_matrix_contraction = &
2436 one_center_density_matrix_contraction + &
2437 contract_one_center_matrix(cpc_h(ispin)%r_coef, &
2438 int_hh(ispin)%r_coef, &
2439 tau_basis_cache%n2oindex) - &
2440 contract_one_center_matrix(cpc_s(ispin)%r_coef, &
2441 int_ss(ispin)%r_coef, &
2442 tau_basis_cache%n2oindex)
2443 int_hh(ispin)%r_coef = 0.0_dp
2444 int_ss(ispin)%r_coef = 0.0_dp
2446 CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
2447 grid_atom, basis_1c, harmonics, nspins)
2448 DO ispin = 1, nspins
2449 one_center_rho_grad_matrix_contraction = &
2450 one_center_rho_grad_matrix_contraction + &
2451 contract_one_center_matrix(cpc_h(ispin)%r_coef, &
2452 int_hh(ispin)%r_coef, &
2453 tau_basis_cache%n2oindex) - &
2454 contract_one_center_matrix(cpc_s(ispin)%r_coef, &
2455 int_ss(ispin)%r_coef, &
2456 tau_basis_cache%n2oindex)
2458 CALL dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, tau_basis_cache, nspins)
2459 DO ispin = 1, nspins
2460 one_center_matrix_contraction = one_center_matrix_contraction + &
2461 contract_one_center_matrix(cpc_h(ispin)%r_coef, &
2462 int_hh(ispin)%r_coef, &
2463 tau_basis_cache%n2oindex) - &
2464 contract_one_center_matrix(cpc_s(ispin)%r_coef, &
2465 int_ss(ispin)%r_coef, &
2466 tau_basis_cache%n2oindex)
2467 IF (.NOT. atom_composite_reference)
THEN
2468 int_hh(ispin)%r_coef = composite_int_h(:, :, ispin)
2469 int_ss(ispin)%r_coef = composite_int_s(:, :, ispin)
2472 DEALLOCATE (composite_int_h, composite_int_s)
2474 CALL release_tau_basis_cache(tau_basis_cache)
2477 CALL para_env%sum(one_center_density_field_contraction)
2478 CALL para_env%sum(one_center_density_matrix_contraction)
2479 CALL para_env%sum(one_center_gradient_field_contraction)
2480 CALL para_env%sum(one_center_matrix_contraction)
2481 CALL para_env%sum(one_center_rho_grad_field_contraction)
2482 CALL para_env%sum(one_center_rho_grad_matrix_contraction)
2483 CALL para_env%sum(one_center_tau_field_contraction)
2484 one_center_field_contraction = one_center_rho_grad_field_contraction + &
2485 one_center_tau_field_contraction
2486 one_center_gradient_matrix_contraction = one_center_rho_grad_matrix_contraction - &
2487 one_center_density_matrix_contraction
2488 one_center_tau_matrix_contraction = one_center_matrix_contraction - &
2489 one_center_rho_grad_matrix_contraction
2490 feature_component_analytic(1) = sum(composite_density_grad*composite_density)
2491 feature_component_analytic(2) = sum(composite_grad_grad*composite_grad)
2492 feature_component_analytic(3) = sum(composite_kin_grad*composite_kin)
2493 feature_component_analytic(4) = sum(feature_component_analytic(1:3))
2494 CALL para_env%sum(feature_component_analytic)
2495 one_center_tensor_contraction = feature_component_analytic(4) - &
2496 smooth_input_contraction
2497 IF (atom_composite_diagnostic)
THEN
2498 DO icomponent = 1, 4
2499 SELECT CASE (icomponent)
2501 composite_density = (1.0_dp + feature_vxc_step)*composite_density
2503 composite_grad = (1.0_dp + feature_vxc_step)*composite_grad
2505 composite_kin = (1.0_dp + feature_vxc_step)*composite_kin
2507 composite_density = (1.0_dp + feature_vxc_step)*composite_density
2508 composite_grad = (1.0_dp + feature_vxc_step)*composite_grad
2509 composite_kin = (1.0_dp + feature_vxc_step)*composite_kin
2512 my_xc_section, para_env, composite_density, composite_grad, composite_kin, &
2513 composite_grid_coords, composite_grid_weights, composite_atomic_grid_weights, &
2514 composite_local_grid_sizes, composite_local_atom_coords, feature_vxc_plus)
2515 SELECT CASE (icomponent)
2517 composite_density = ((1.0_dp - feature_vxc_step)/ &
2518 (1.0_dp + feature_vxc_step))*composite_density
2520 composite_grad = ((1.0_dp - feature_vxc_step)/ &
2521 (1.0_dp + feature_vxc_step))*composite_grad
2523 composite_kin = ((1.0_dp - feature_vxc_step)/ &
2524 (1.0_dp + feature_vxc_step))*composite_kin
2526 composite_density = ((1.0_dp - feature_vxc_step)/ &
2527 (1.0_dp + feature_vxc_step))*composite_density
2528 composite_grad = ((1.0_dp - feature_vxc_step)/ &
2529 (1.0_dp + feature_vxc_step))*composite_grad
2530 composite_kin = ((1.0_dp - feature_vxc_step)/ &
2531 (1.0_dp + feature_vxc_step))*composite_kin
2534 my_xc_section, para_env, composite_density, composite_grad, composite_kin, &
2535 composite_grid_coords, composite_grid_weights, composite_atomic_grid_weights, &
2536 composite_local_grid_sizes, composite_local_atom_coords, feature_vxc_minus)
2537 SELECT CASE (icomponent)
2539 composite_density = composite_density/(1.0_dp - feature_vxc_step)
2541 composite_grad = composite_grad/(1.0_dp - feature_vxc_step)
2543 composite_kin = composite_kin/(1.0_dp - feature_vxc_step)
2545 composite_density = composite_density/(1.0_dp - feature_vxc_step)
2546 composite_grad = composite_grad/(1.0_dp - feature_vxc_step)
2547 composite_kin = composite_kin/(1.0_dp - feature_vxc_step)
2549 feature_component_fd(icomponent) = &
2550 (feature_vxc_plus - feature_vxc_minus)/(2.0_dp*feature_vxc_step)
2552 feature_vxc_analytic = feature_component_analytic(4)
2553 feature_vxc_fd = feature_component_fd(4)
2556 WRITE (unit=iw, fmt=
"(T2,A,1X,I0)") &
2557 "SKALA_GPW| Atom-composite reference components", atom_composite_components
2558 WRITE (unit=iw, fmt=
"(T2,A,1X,ES24.16)") &
2559 "SKALA_GPW| Atom-composite reference electrons", atom_composite_nelec
2560 WRITE (unit=iw, fmt=
"(T2,A,1X,ES24.16)") &
2561 "SKALA_GPW| Atom-composite reference XC energy", atom_composite_exc
2562 WRITE (unit=iw, fmt=
"(T2,A,1X,ES24.16)") &
2563 "SKALA_GPW| Atom-composite feature VXC contraction", feature_vxc_analytic
2564 WRITE (unit=iw, fmt=
"(T2,A,1X,ES24.16)") &
2565 "SKALA_GPW| Atom-composite feature finite difference", feature_vxc_fd
2566 WRITE (unit=iw, fmt=
"(T2,A,1X,ES24.16)") &
2567 "SKALA_GPW| Atom-composite feature VXC error", &
2568 feature_vxc_analytic - feature_vxc_fd
2569 WRITE (unit=iw, fmt=
"(T2,A,3(1X,ES24.16))") &
2570 "SKALA_GPW| Atom-composite smooth PW adjoint", smooth_input_contraction, &
2571 smooth_grid_contraction, smooth_grid_contraction - smooth_input_contraction
2572 WRITE (unit=iw, fmt=
"(T2,A,4(1X,ES24.16))") &
2573 "SKALA_GPW| Atom-composite one-center adjoint", one_center_tensor_contraction, &
2574 one_center_field_contraction, one_center_matrix_contraction, &
2575 one_center_matrix_contraction - one_center_tensor_contraction
2576 WRITE (unit=iw, fmt=
"(T2,A,4(1X,ES24.16))") &
2577 "SKALA_GPW| Atom-composite one-center channels", &
2578 one_center_rho_grad_field_contraction, one_center_rho_grad_matrix_contraction, &
2579 one_center_tau_field_contraction, one_center_tau_matrix_contraction
2580 WRITE (unit=iw, fmt=
"(T2,A,4(1X,ES24.16))") &
2581 "SKALA_GPW| Atom-composite one-center rho-gradient", &
2582 one_center_density_field_contraction, one_center_density_matrix_contraction, &
2583 one_center_gradient_field_contraction, one_center_gradient_matrix_contraction
2584 DO icomponent = 1, 3
2585 WRITE (unit=iw, fmt=
"(T2,A,1X,I0,3(1X,ES24.16))") &
2586 "SKALA_GPW| Atom-composite component VXC", icomponent, &
2587 feature_component_analytic(icomponent), feature_component_fd(icomponent), &
2588 feature_component_analytic(icomponent) - feature_component_fd(icomponent)
2592 IF (atom_composite_reference)
THEN
2593 exc1 = atom_composite_exc
2594 IF (native_grid_diagnostics)
THEN
2595 IF (composite_nflat > 0)
THEN
2596 composite_density_min = minval(composite_density)
2597 composite_density_max = maxval(composite_density)
2598 composite_kin_min = minval(composite_kin)
2599 composite_kin_max = maxval(composite_kin)
2600 composite_grad_max = maxval(abs(composite_grad))
2602 composite_density_min = huge(1.0_dp)
2603 composite_density_max = -huge(1.0_dp)
2604 composite_kin_min = huge(1.0_dp)
2605 composite_kin_max = -huge(1.0_dp)
2606 composite_grad_max = 0.0_dp
2608 composite_tau_integral = &
2609 sum(composite_grid_weights*sum(composite_kin, dim=2))
2610 CALL para_env%min(composite_density_min)
2611 CALL para_env%max(composite_density_max)
2612 CALL para_env%min(composite_kin_min)
2613 CALL para_env%max(composite_kin_max)
2614 CALL para_env%max(composite_grad_max)
2615 CALL para_env%sum(composite_tau_integral)
2619 WRITE (unit=iw, fmt=
"(T2,A,1X,ES24.16)") &
2620 "SKALA_GPW| Active atom-composite XC energy", atom_composite_exc
2621 IF (native_grid_diagnostics)
THEN
2622 WRITE (unit=iw, fmt=
"(T2,A,1X,ES24.16)") &
2623 "SKALA_GPW| Active atom-composite electrons", atom_composite_nelec
2624 WRITE (unit=iw, fmt=
"(T2,A,2(1X,ES24.16))") &
2625 "SKALA_GPW| Active atom-composite density range", &
2626 composite_density_min, composite_density_max
2627 WRITE (unit=iw, fmt=
"(T2,A,2(1X,ES24.16))") &
2628 "SKALA_GPW| Active atom-composite tau range", &
2629 composite_kin_min, composite_kin_max
2630 WRITE (unit=iw, fmt=
"(T2,A,1X,ES24.16)") &
2631 "SKALA_GPW| Active atom-composite tau integral", &
2632 composite_tau_integral
2633 WRITE (unit=iw, fmt=
"(T2,A,1X,ES24.16)") &
2634 "SKALA_GPW| Active atom-composite max gradient", &
2640 DEALLOCATE (composite_atomic_grid_sizes, composite_atom_kind, &
2641 composite_atom_kind_index, composite_atom_start, composite_atom_end, &
2642 composite_atom_coords, composite_local_atoms, composite_local_grid_sizes, &
2643 composite_local_atom_coords, composite_grid_atom, &
2644 composite_partition_weights, composite_partition_atom_coords, &
2645 composite_distances, composite_density, composite_grad, composite_kin, &
2646 composite_smooth_density_cache, composite_smooth_gradient_cache, &
2647 composite_smooth_kin_cache, &
2648 composite_density_grad, composite_grad_grad, composite_kin_grad, &
2649 composite_grid_coords, composite_grid_weights, &
2650 composite_base_grid_weights, composite_atomic_grid_weights)
2652 DEALLOCATE (composite_smooth_rhoa, composite_smooth_rhob, &
2653 composite_smooth_tau_a, composite_smooth_tau_b)
2655 DEALLOCATE (composite_smooth_drhoa(idir)%array, &
2656 composite_smooth_drhob(idir)%array)
2659 DEALLOCATE (composite_smooth_rho, composite_smooth_tau)
2661 DEALLOCATE (composite_smooth_drho(idir)%array)
2666 IF (.NOT. atom_composite_reference)
CALL para_env%sum(exc1)
2668 IF (
ASSOCIATED(rho_h))
DEALLOCATE (rho_h)
2669 IF (
ASSOCIATED(rho_s))
DEALLOCATE (rho_s)
2670 IF (
ASSOCIATED(vxc_h))
DEALLOCATE (vxc_h)
2671 IF (
ASSOCIATED(vxc_s))
DEALLOCATE (vxc_s)
2673 IF (gradient_f)
THEN
2674 IF (
ASSOCIATED(drho_h))
DEALLOCATE (drho_h)
2675 IF (
ASSOCIATED(drho_s))
DEALLOCATE (drho_s)
2676 IF (
ASSOCIATED(vxg_h))
DEALLOCATE (vxg_h)
2677 IF (
ASSOCIATED(vxg_s))
DEALLOCATE (vxg_s)
2681 IF (
ASSOCIATED(tau_h))
DEALLOCATE (tau_h)
2682 IF (
ASSOCIATED(tau_s))
DEALLOCATE (tau_s)
2683 IF (
ASSOCIATED(vtau_h))
DEALLOCATE (vtau_h)
2684 IF (
ASSOCIATED(vtau_s))
DEALLOCATE (vtau_s)
2689 CALL timestop(handle)