Handle periodic image coordinates for native-grid SKALA (#5321)

Co-authored-by: Thomas D. Kuehne <tkuehne@cp2k.org>
This commit is contained in:
Dynamics of Condensed Matter 2026-05-31 13:29:29 +02:00 committed by GitHub
parent 1446772304
commit 59397cdcc8
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
7 changed files with 371 additions and 14 deletions

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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},