diff --git a/src/skala_gpw_features.F b/src/skala_gpw_features.F index 23b87ce893..dd065741e6 100644 --- a/src/skala_gpw_features.F +++ b/src/skala_gpw_features.F @@ -64,7 +64,7 @@ MODULE skala_gpw_features CONTAINS ! ************************************************************************************************** -!> \brief Build a flat SKALA molecular feature dictionary from a local GPW grid. +!> \brief Build a flat SKALA feature dictionary from a local GPW grid. !> \param features ... !> \param rho_set ... !> \param rho_r ... @@ -97,6 +97,7 @@ CONTAINS LOGICAL :: my_requires_coordinate_grad, & my_requires_grad REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: global_feature, local_feature + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atom_coords_pbc REAL(KIND=dp), DIMENSION(3) :: grid_point, owner_coord REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho, rhoa, rhob, tau_a, tau_b, tau_total TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drhoa, drhob @@ -123,12 +124,16 @@ CONTAINS CALL skala_gpw_feature_release(features) ALLOCATE (local_owner(nflat_local), local_feature(nreal_per_point*nflat_local), & - feature_counts(nproc), feature_displs(nproc), real_counts(nproc), real_displs(nproc)) + feature_counts(nproc), feature_displs(nproc), real_counts(nproc), & + real_displs(nproc), atom_coords_pbc(3, natom)) ALLOCATE (features%feature_index(bo(1, 1):bo(2, 1), & bo(1, 2):bo(2, 2), & bo(1, 3):bo(2, 3))) features%feature_index = 0 local_feature = 0.0_dp + DO iatom = 1, natom + atom_coords_pbc(:, iatom) = pbc(particle_set(iatom)%r, cell, positive_range=.TRUE.) + END DO IF (nspins == 1) THEN CALL xc_rho_set_get(rho_set, rho=rho, drho=drho, tau=tau_total) @@ -144,13 +149,13 @@ CONTAINS local_row = local_row + 1 real_base = nreal_per_point*(local_row - 1) grid_point = grid_coordinate(pw_grid, [i, j, k]) - owner = nearest_atom(grid_point, particle_set, cell) + owner = nearest_atom(grid_point, atom_coords_pbc, cell) local_owner(local_row) = owner features%feature_index(i, j, k) = local_row - owner_coord = pbc(particle_set(owner)%r, cell, positive_range=.TRUE.) - local_feature(real_base + 11:real_base + 13) = owner_coord + & - pbc(owner_coord, grid_point, cell) + owner_coord = atom_coords_pbc(:, owner) + local_feature(real_base + 11:real_base + 13) = & + nearest_image_coordinate(owner_coord, grid_point, cell) local_feature(real_base + 14) = pw_grid%dvol IF (PRESENT(weights)) THEN IF (ASSOCIATED(weights)) THEN @@ -227,7 +232,7 @@ CONTAINS features%atomic_grid_weights = pw_grid%dvol DO iatom = 1, natom - features%coarse_0_atomic_coords(:, iatom) = pbc(particle_set(iatom)%r, cell, positive_range=.TRUE.) + features%coarse_0_atomic_coords(:, iatom) = atom_coords_pbc(:, iatom) END DO DO ipt = 1, nflat @@ -269,8 +274,9 @@ CONTAINS CALL add_feature_tensors(features, my_requires_grad, my_requires_coordinate_grad) features%active = .TRUE. - DEALLOCATE (atom_offset, atom_position, feature_counts, feature_displs, global_feature, & - global_owner, local_feature, local_owner, local_to_global, real_counts, real_displs) + DEALLOCATE (atom_coords_pbc, atom_offset, atom_position, feature_counts, feature_displs, & + global_feature, global_owner, local_feature, local_owner, local_to_global, & + real_counts, real_displs) END SUBROUTINE skala_gpw_feature_build @@ -374,13 +380,13 @@ CONTAINS ! ************************************************************************************************** !> \brief Assign a grid point to the nearest periodic atom. !> \param grid_point ... -!> \param particle_set ... +!> \param atom_coords ... !> \param cell ... !> \return ... ! ************************************************************************************************** - FUNCTION nearest_atom(grid_point, particle_set, cell) RESULT(owner) + FUNCTION nearest_atom(grid_point, atom_coords, cell) RESULT(owner) REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: grid_point - TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: atom_coords TYPE(cell_type), POINTER :: cell INTEGER :: owner @@ -390,8 +396,8 @@ CONTAINS owner = 1 best_r2 = HUGE(1.0_dp) - DO iatom = 1, SIZE(particle_set) - rij = pbc(grid_point, particle_set(iatom)%r, cell) + DO iatom = 1, SIZE(atom_coords, 2) + rij = pbc(grid_point, atom_coords(:, iatom), cell) r2 = SUM(rij**2) IF (r2 < best_r2) THEN best_r2 = r2 @@ -401,4 +407,20 @@ CONTAINS END FUNCTION nearest_atom +! ************************************************************************************************** +!> \brief Return the nearest periodic image of an atom center to a GPW grid point. +!> \param atom_coord ... +!> \param grid_point ... +!> \param cell ... +!> \return ... +! ************************************************************************************************** + FUNCTION nearest_image_coordinate(atom_coord, grid_point, cell) RESULT(coord) + REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: atom_coord, grid_point + TYPE(cell_type), POINTER :: cell + REAL(KIND=dp), DIMENSION(3) :: coord + + coord = atom_coord + pbc(atom_coord, grid_point, cell) + + END FUNCTION nearest_image_coordinate + END MODULE skala_gpw_features diff --git a/tests/QS/regtest-gauxc/AR4_NATIVE_SKALA_GPW_IMAGE_COORDS.inp b/tests/QS/regtest-gauxc/AR4_NATIVE_SKALA_GPW_IMAGE_COORDS.inp new file mode 100644 index 0000000000..4981399bca --- /dev/null +++ b/tests/QS/regtest-gauxc/AR4_NATIVE_SKALA_GPW_IMAGE_COORDS.inp @@ -0,0 +1,61 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT_NAME AR4_NATIVE_SKALA_GPW_IMAGE_COORDS + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + MULTIPLICITY 1 + POTENTIAL_FILE_NAME GTH_POTENTIALS + UKS FALSE + &MGRID + CUTOFF 150 + REL_CUTOFF 30 + &END MGRID + &POISSON + PERIODIC XYZ + &END POISSON + &QS + EPS_DEFAULT 1.0E-8 + EXTRAPOLATION USE_PREV_WF + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-5 + IGNORE_CONVERGENCE_FAILURE T + MAX_SCF 1 + SCF_GUESS ATOMIC + &DIAGONALIZATION + &END DIAGONALIZATION + &END SCF + &XC + &XC_FUNCTIONAL + &GAUXC + FUNCTIONAL PBE + MODEL SKALA + NATIVE_GRID T + NATIVE_GRID_DIAGNOSTICS T + &END GAUXC + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.30 5.30 5.30 + PERIODIC XYZ + &END CELL + &COORD + Ar 5.400000 5.400000 5.400000 + Ar -5.200000 2.750000 8.050000 + Ar 8.050000 -5.200000 2.750000 + Ar 2.750000 8.050000 -5.200000 + &END COORD + &KIND Ar + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q8 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-gauxc/H2O_NATIVE_SKALA_GPW_IMAGE_COORDS.inp b/tests/QS/regtest-gauxc/H2O_NATIVE_SKALA_GPW_IMAGE_COORDS.inp new file mode 100644 index 0000000000..10747dfd4f --- /dev/null +++ b/tests/QS/regtest-gauxc/H2O_NATIVE_SKALA_GPW_IMAGE_COORDS.inp @@ -0,0 +1,64 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT_NAME H2O_NATIVE_SKALA_GPW_IMAGE_COORDS + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + MULTIPLICITY 1 + POTENTIAL_FILE_NAME GTH_POTENTIALS + UKS FALSE + &MGRID + CUTOFF 150 + REL_CUTOFF 30 + &END MGRID + &POISSON + PERIODIC XYZ + &END POISSON + &QS + EPS_DEFAULT 1.0E-8 + EXTRAPOLATION USE_PREV_WF + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-5 + IGNORE_CONVERGENCE_FAILURE T + MAX_SCF 1 + SCF_GUESS ATOMIC + &DIAGONALIZATION + &END DIAGONALIZATION + &END SCF + &XC + &XC_FUNCTIONAL + &GAUXC + FUNCTIONAL PBE + MODEL SKALA + NATIVE_GRID T + NATIVE_GRID_DIAGNOSTICS T + &END GAUXC + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + PERIODIC XYZ + &END CELL + &COORD + O 9.000000 -3.000000 3.000000 + H 9.000000 -2.241398 3.504284 + H 9.000000 -3.758602 3.504284 + &END COORD + &KIND O + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q6 + &END KIND + &KIND H + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-gauxc/H2P_NATIVE_SKALA_GPW_UKS_PBC_FORCE_DEBUG.inp b/tests/QS/regtest-gauxc/H2P_NATIVE_SKALA_GPW_UKS_PBC_FORCE_DEBUG.inp new file mode 100644 index 0000000000..3c2ed09a33 --- /dev/null +++ b/tests/QS/regtest-gauxc/H2P_NATIVE_SKALA_GPW_UKS_PBC_FORCE_DEBUG.inp @@ -0,0 +1,68 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT_NAME H2P_NATIVE_SKALA_GPW_UKS_PBC_FORCE_DEBUG + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + CHECK_ATOM_FORCE 2 z + DEBUG_FORCES T + DEBUG_STRESS_TENSOR F + DX 1.0E-4 + EPS_NO_ERROR_CHECK 5.0E-5 + STOP_ON_MISMATCH T +&END DEBUG + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + CHARGE 1 + MULTIPLICITY 2 + POTENTIAL_FILE_NAME GTH_POTENTIALS + UKS TRUE + &MGRID + CUTOFF 40 + REL_CUTOFF 10 + &END MGRID + &POISSON + PERIODIC XYZ + &END POISSON + &QS + EPS_DEFAULT 1.0E-9 + EXTRAPOLATION USE_PREV_WF + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-5 + MAX_SCF 20 + SCF_GUESS ATOMIC + &DIAGONALIZATION + &END DIAGONALIZATION + &END SCF + &XC + &XC_FUNCTIONAL + &GAUXC + FUNCTIONAL PBE + MODEL SKALA + NATIVE_GRID T + NATIVE_GRID_DIAGNOSTICS F + &END GAUXC + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC XYZ + &END CELL + &COORD + H 2.0 2.0 3.7 + H 2.0 2.0 0.5 + &END COORD + &KIND H + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-gauxc/H2_NATIVE_SKALA_GPW_PBC_FORCE_DEBUG.inp b/tests/QS/regtest-gauxc/H2_NATIVE_SKALA_GPW_PBC_FORCE_DEBUG.inp new file mode 100644 index 0000000000..ac4cda95ee --- /dev/null +++ b/tests/QS/regtest-gauxc/H2_NATIVE_SKALA_GPW_PBC_FORCE_DEBUG.inp @@ -0,0 +1,67 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT_NAME H2_NATIVE_SKALA_GPW_PBC_FORCE_DEBUG + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + CHECK_ATOM_FORCE 2 z + DEBUG_FORCES T + DEBUG_STRESS_TENSOR F + DX 1.0E-4 + EPS_NO_ERROR_CHECK 5.0E-5 + STOP_ON_MISMATCH T +&END DEBUG + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + MULTIPLICITY 1 + POTENTIAL_FILE_NAME GTH_POTENTIALS + UKS FALSE + &MGRID + CUTOFF 40 + REL_CUTOFF 10 + &END MGRID + &POISSON + PERIODIC XYZ + &END POISSON + &QS + EPS_DEFAULT 1.0E-9 + EXTRAPOLATION USE_PREV_WF + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-5 + MAX_SCF 20 + SCF_GUESS ATOMIC + &DIAGONALIZATION + &END DIAGONALIZATION + &END SCF + &XC + &XC_FUNCTIONAL + &GAUXC + FUNCTIONAL PBE + MODEL SKALA + NATIVE_GRID T + NATIVE_GRID_DIAGNOSTICS F + &END GAUXC + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC XYZ + &END CELL + &COORD + H 2.0 2.0 3.7 + H 2.0 2.0 0.5 + &END COORD + &KIND H + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-gauxc/OH_NATIVE_SKALA_GPW_UKS_IMAGE_COORDS.inp b/tests/QS/regtest-gauxc/OH_NATIVE_SKALA_GPW_UKS_IMAGE_COORDS.inp new file mode 100644 index 0000000000..9b9663375f --- /dev/null +++ b/tests/QS/regtest-gauxc/OH_NATIVE_SKALA_GPW_UKS_IMAGE_COORDS.inp @@ -0,0 +1,63 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT_NAME OH_NATIVE_SKALA_GPW_UKS_IMAGE_COORDS + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + MULTIPLICITY 2 + POTENTIAL_FILE_NAME GTH_POTENTIALS + UKS TRUE + &MGRID + CUTOFF 150 + REL_CUTOFF 30 + &END MGRID + &POISSON + PERIODIC XYZ + &END POISSON + &QS + EPS_DEFAULT 1.0E-8 + EXTRAPOLATION USE_PREV_WF + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-5 + IGNORE_CONVERGENCE_FAILURE T + MAX_SCF 1 + SCF_GUESS ATOMIC + &DIAGONALIZATION + &END DIAGONALIZATION + &END SCF + &XC + &XC_FUNCTIONAL + &GAUXC + FUNCTIONAL PBE + MODEL SKALA + NATIVE_GRID T + NATIVE_GRID_DIAGNOSTICS T + &END GAUXC + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + PERIODIC XYZ + &END CELL + &COORD + O 9.0 -3.0 3.0 + H 9.0 -3.0 3.97 + &END COORD + &KIND O + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q6 + &END KIND + &KIND H + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-gauxc/TEST_FILES.toml b/tests/QS/regtest-gauxc/TEST_FILES.toml index eb46c61a7d..1ca6761ca5 100644 --- a/tests/QS/regtest-gauxc/TEST_FILES.toml +++ b/tests/QS/regtest-gauxc/TEST_FILES.toml @@ -5,14 +5,22 @@ "H2_SKALA_ENERGY_CHUNKED.inp" = [{matcher="E_total", tol=1e-8, ref=-0.979366068078563}] "H2_NATIVE_SKALA_GPW.inp" = [{matcher="E_total", tol=1e-8, ref=-0.976732566963415}] "H2_NATIVE_SKALA_GPW_FORCE.inp" = [{matcher="M072", tol=1e-5, ref=3.03111884E-04}] +"H2_NATIVE_SKALA_GPW_PBC_FORCE_DEBUG.inp" = [{matcher="DEBUG_force_sum", tol=5e-5, ref=0.0}] +"H2P_NATIVE_SKALA_GPW_UKS_PBC_FORCE_DEBUG.inp" = [{matcher="DEBUG_force_sum", tol=5e-5, ref=0.0}] "H2O_NATIVE_GPW_PBE_REFERENCE.inp" = [{matcher="E_total", tol=1e-8, ref=-17.200873708850686}] "H2O_NATIVE_SKALA_GPW.inp" = [{matcher="E_total", tol=5e-6, ref=-17.067481570036010}, {matcher="SKALA_GPW_feature_electrons", tol=1e-8, ref=8.00008704097}, {matcher="SKALA_GPW_feature_spin_moment", tol=1e-12, ref=0.0}] +"H2O_NATIVE_SKALA_GPW_IMAGE_COORDS.inp" = [{matcher="E_total", tol=5e-6, ref=-17.067481570036010}, + {matcher="SKALA_GPW_feature_electrons", tol=1e-8, ref=8.00008704097}, + {matcher="SKALA_GPW_feature_spin_moment", tol=1e-12, ref=0.0}] "OH_NATIVE_GPW_PBE_UKS_REFERENCE.inp" = [{matcher="E_total", tol=5e-6, ref=-16.524536848018680}] "OH_NATIVE_SKALA_GPW_UKS.inp" = [{matcher="E_total", tol=5e-6, ref=-16.408200534919416}, {matcher="SKALA_GPW_feature_electrons", tol=1e-8, ref=7.00008702590}, {matcher="SKALA_GPW_feature_spin_moment", tol=1e-8, ref=1.00001243227}] +"OH_NATIVE_SKALA_GPW_UKS_IMAGE_COORDS.inp" = [{matcher="E_total", tol=5e-6, ref=-16.408200534919416}, + {matcher="SKALA_GPW_feature_electrons", tol=1e-8, ref=7.00008702590}, + {matcher="SKALA_GPW_feature_spin_moment", tol=1e-8, ref=1.00001243227}] "AR4_NATIVE_GPW_PBE_REFERENCE.inp" = [{matcher="E_total", tol=1e-8, ref=-84.254262954861559}] "AR4_NATIVE_GPW_PBE_WRAPPED_REFERENCE.inp" = [{matcher="E_total", tol=1e-8, ref=-84.254262954861559}] "NE4_NATIVE_GPW_PBE_REFERENCE.inp" = [{matcher="E_total", tol=1e-8, ref=-139.361323676360371}] @@ -24,6 +32,10 @@ {matcher="SKALA_GPW_feature_electrons", tol=1e-8, ref=31.9999999584}, {matcher="SKALA_GPW_feature_spin_moment", tol=1e-12, ref=0.0}, {matcher="SKALA_GPW_feature_weight_sum", tol=1e-8, ref=1004.67180773}] +"AR4_NATIVE_SKALA_GPW_IMAGE_COORDS.inp" = [{matcher="E_total", tol=5e-6, ref=-84.344774381927250}, + {matcher="SKALA_GPW_feature_electrons", tol=1e-8, ref=31.9999999584}, + {matcher="SKALA_GPW_feature_spin_moment", tol=1e-12, ref=0.0}, + {matcher="SKALA_GPW_feature_weight_sum", tol=1e-8, ref=1004.67180773}] "AR4_NATIVE_SKALA_GPW_WRAPPED.inp" = [{matcher="E_total", tol=5e-6, ref=-84.344779494109700}, {matcher="SKALA_GPW_feature_electrons", tol=1e-8, ref=31.9999999584}, {matcher="SKALA_GPW_feature_spin_moment", tol=1e-12, ref=0.0},