comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge 8be1dfa50f into a0df33dd6a

Browse source
This commit is contained in:
Dhruv 2026-07-21 00:43:33 +02:00 • committed by GitHub
commit e284f4eda7
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
17 changed files with 4057 additions and 876 deletions
Download patch file Download diff file Expand all files Collapse all files

80
docs/methods/machine_learning/nnp.md
View file

@ -1,8 +1,84 @@
# Neural Network Potentials
Unfortunately no one has gotten around to writing this page yet :-(
CP2K supports Behler-Parrinello high-dimensional neural network potentials (HDNNPs) with
atom-centred symmetry functions (ACSF) as descriptors. NNP can drive any CP2K run that goes through
`FORCE_EVAL`, from single-point energies through geometry optimisation, molecular dynamics, biased
MD via free-energy methods, and committee-of-models extrapolation diagnostics. The same NNP path
also backs the helium-solvent interaction in path-integral runs.
In the meantime, the following links might be helpful:
The implementation reads networks trained with the [n2p2](https://github.com/CompPhysVienna/n2p2)
format (`input.nn`, `scaling.data`, `weights.<element>.data`).
## Input
The NNP method is selected with `METHOD NNP` inside `FORCE_EVAL`. The network files and one or more
model definitions are configured under `&NNP`. The complete option list is generated from the source
and is linked from the input reference manual.
A minimal committee-NNP MD input looks like
```
&FORCE_EVAL
METHOD NNP
&NNP
NNP_INPUT_FILE_NAME nnp-1/input.nn
SCALE_FILE_NAME nnp-1/scaling.data
&MODEL
WEIGHTS nnp-1/weights
&END MODEL
&MODEL
WEIGHTS nnp-2/weights
&END MODEL
! ... up to 8 typical for committee error bars ...
&END NNP
&SUBSYS
&CELL
ABC [angstrom] 12.42 12.42 12.42
&END CELL
&COORD
! ...
&END COORD
&END SUBSYS
&END FORCE_EVAL
```
Worked examples covering NVT, NPT, biased MD, and restarts live under
[`tests/NNP/regtest-1/`](https://github.com/cp2k/cp2k/tree/master/tests/NNP/regtest-1). The
path-integral helium-solute coupling is exercised by
[`tests/Pimd/regtest-2/water_in_helium_nnp.inp`](https://github.com/cp2k/cp2k/tree/master/tests/Pimd/regtest-2).
## Tuning
The default ACSF spline grid (`RAD_SPLINE_N 8192`) is calibrated for radial cutoffs around 12 bohr.
Models with substantially larger cutoffs or stricter accuracy targets may want a larger value. The
cubic-Hermite value residual scales as O(1/n^4) and the force (derivative) residual as O(1/n^3); at
the default grid both stay near machine precision for a 12 bohr cutoff (value error ~1e-14, force
error ~1e-10), with the force term the binding constraint as the cutoff grows.
The cell-list neighbour search uses a Verlet skin so the chain is only rebuilt when an atom drifts
more than `skin/2` between force evaluations. By default the skin auto-selects
`MIN(0.5 bohr, 0.1 * cutoff)`. Long stable trajectories can raise it to reduce the rebuild rate at
the cost of a larger per-atom neighbour list:
```
&NNP
! ...
VERLET_SKIN [bohr] 1.0
&END NNP
```
This is the analogue of LAMMPS's `neighbor <skin> bin` command. A negative value (the default)
selects the automatic heuristic; useful upper bounds sit near half the smallest perpendicular cell
width. See the input reference manual for the full list of NNP tuning keywords.
## Parallelism
NNP supports both MPI and OpenMP parallelism. MPI distributes the atoms across ranks, and OpenMP
parallelises the per-atom descriptor and force loops inside each rank.
## References
For background on the method and the existing implementation, the following links might be helpful
- <https://www.cp2k.org/tools:aml>
- <https://doi.org/10.1063/5.0160326>

2
src/CMakeLists.txt
View file

@ -431,10 +431,12 @@ list(
negf_subgroup_types.F
negf_vectors.F
nnp_acsf.F
nnp_cell_list.F
nnp_environment.F
nnp_environment_types.F
nnp_force.F
nnp_model.F
nnp_neighbor_interface.F
local_gemm_api.F
offload/offload_api.F
optbas_fenv_manipulation.F

49
src/input_cp2k_nnp.F
View file

@ -8,6 +8,7 @@
! **************************************************************************************************
!> \brief Creates the NNP section of the input
!> \author Christoph Schran (christoph.schran@rub.de)
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
!> \date 2020-10-10
! **************************************************************************************************
MODULE input_cp2k_nnp
@ -28,6 +29,7 @@ MODULE input_cp2k_nnp
section_release,&
section_type
USE input_val_types, ONLY: char_t,&
integer_t,&
real_t
USE kinds, ONLY: dp
#include "./base/base_uses.f90"
@ -35,7 +37,7 @@ MODULE input_cp2k_nnp
IMPLICIT NONE
PRIVATE
LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'input_cp2k_nnp'
PUBLIC :: create_nnp_section
@ -58,7 +60,7 @@ CONTAINS
CALL section_create(section, __LOCATION__, name="NNP", &
description="This section contains all information to run a "// &
"Neural Network Potential (NNP) calculation.", &
n_keywords=3, n_subsections=3, repeats=.FALSE., &
n_keywords=5, n_subsections=3, repeats=.FALSE., &
citations=[Behler2007, Behler2011, Schran2020a, Schran2020b])
NULLIFY (subsection, subsubsection, keyword)
@ -76,6 +78,49 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="RAD_SPLINE_N", &
description="Number of knots per radial group in "// &
"the cubic-Hermite spline tables that tabulate "// &
"the radial symmetry-function product "// &
"y(r) = exp(-eta*(r-rs)^2) * fcut(r). "// &
"The value error scales as O(1/n^4) and the "// &
"force (derivative) error as O(1/n^3); the "// &
"default keeps both inside the NNP regression "// &
"tolerance for a radial cutoff of ~12 bohr "// &
"(value ~1e-14, force ~1e-10). Larger cutoffs or "// &
"stricter tolerances may need a larger n, with "// &
"the force term the binding constraint. Memory "// &
"cost scales linearly in n per group.", &
repeats=.FALSE., &
n_var=1, &
type_of_var=integer_t, &
default_i_val=8192, &
usage="RAD_SPLINE_N 8192")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="VERLET_SKIN", &
description="Verlet skin distance for the NNP descriptor "// &
"neighbour cell-list. The neighbour-list cutoff is "// &
"(symmetry-function cutoff + skin); the cell-list "// &
"chain is rebuilt only when an atom drifts more "// &
"than skin/2 from its rebuild-time position, "// &
"analogous to the LAMMPS 'neighbor <skin> bin' "// &
"command. Larger skin reduces the rebuild rate "// &
"but enlarges the per-atom neighbour list. "// &
"A negative value (default) selects the "// &
"automatic heuristic MIN(0.5 bohr, 0.1*cutoff). "// &
"Useful upper bound is half the smallest "// &
"perpendicular cell width.", &
repeats=.FALSE., &
n_var=1, &
type_of_var=real_t, &
default_r_val=-1.0_dp, &
unit_str="bohr", &
usage="VERLET_SKIN [bohr] 0.5")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
! BIAS subsection
CALL section_create(subsection, __LOCATION__, name="BIAS", &
description="Section to bias the committee disagreement (sigma) by "// &

35
src/motion/helium_interactions.F
View file

@ -34,8 +34,10 @@ MODULE helium_interactions
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type
USE kinds, ONLY: dp
USE nnp_acsf, ONLY: nnp_calc_acsf
USE nnp_acsf, ONLY: nnp_calc_acsf,&
nnp_prepare_neighbor_cache
USE nnp_environment_types, ONLY: nnp_type
USE nnp_force, ONLY: nnp_scatter_dgdr_to_forces
USE nnp_model, ONLY: nnp_gradients,&
nnp_predict
USE physcon, ONLY: angstrom,&
@ -682,10 +684,8 @@ CONTAINS
OPTIONAL, POINTER :: force
INTEGER :: i, i_com, ig, ind, ind_he, j, k, m
LOGICAL :: extrapolate
REAL(KIND=dp) :: rsqr, rvect(3), threshold
REAL(KIND=dp) :: rsqr, rvect(3)
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: denergydsym
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsymdxyz
TYPE(cp_logger_type), POINTER :: logger
TYPE(nnp_type), POINTER :: nnp
TYPE(section_vals_type), POINTER :: print_section
@ -697,9 +697,6 @@ CONTAINS
helium%nnp%myforce(:, :, :) = 0.0_dp
END IF
extrapolate = .FALSE.
threshold = 0.0001d0
!fill coord array
ig = 1
DO i = 1, helium%nnp%n_ele
@ -739,6 +736,8 @@ CONTAINS
! reset flag if there's an extrapolation to report:
helium%nnp%output_expol = .FALSE.
nnp => helium%nnp
CALL nnp_prepare_neighbor_cache(nnp)
! calc atomic contribution to energy and force
!NOTE corresponds to nnp_force line with parallelization:
@ -750,15 +749,12 @@ CONTAINS
!reset input nodes and grads of ele(ind):
helium%nnp%arc(ind)%layer(1)%node(:) = 0.0_dp
nnp => helium%nnp ! work around wrong INTENT of nnp_calc_acsf
IF (PRESENT(force)) THEN
helium%nnp%arc(ind)%layer(1)%node_grad(:) = 0.0_dp
ALLOCATE (dsymdxyz(3, helium%nnp%arc(ind)%n_nodes(1), helium%nnp%num_atoms))
ALLOCATE (denergydsym(helium%nnp%arc(ind)%n_nodes(1)))
dsymdxyz(:, :, :) = 0.0_dp
CALL nnp_calc_acsf(nnp, i, dsymdxyz)
CALL nnp_calc_acsf(nnp, i, .TRUE.)
ELSE
CALL nnp_calc_acsf(nnp, i)
CALL nnp_calc_acsf(nnp, i, .FALSE.)
END IF
! input nodes filled, perform prediction:
@ -771,16 +767,12 @@ CONTAINS
IF (PRESENT(force)) THEN
denergydsym(:) = 0.0_dp
CALL nnp_gradients(helium%nnp%arc(ind), helium%nnp, i_com, denergydsym)
DO j = 1, helium%nnp%arc(ind)%n_nodes(1)
DO k = 1, helium%nnp%num_atoms
DO m = 1, 3
helium%nnp%myforce(m, k, i_com) = helium%nnp%myforce(m, k, i_com) &
- denergydsym(j)*dsymdxyz(m, j, k)
END DO
END DO
END DO
! Per-element workspace dGdr scatter; shared with the main NNP
! force path (see nnp_scatter_dgdr_to_forces in nnp_force.F).
CALL nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym, &
helium%nnp%myforce(:, :, i_com))
END IF
END DO ! end loop over committee members
@ -788,7 +780,6 @@ CONTAINS
!deallocate memory
IF (PRESENT(force)) THEN
DEALLOCATE (denergydsym)
DEALLOCATE (dsymdxyz)
END IF
END DO ! end loop over num_atoms

1993
src/nnp_acsf.F
View file

File diff suppressed because it is too large Load diff

563
src/nnp_cell_list.F Normal file
View file

@ -0,0 +1,563 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Linked-cell neighbour finder with a Verlet skin for the NNP descriptor.
!> Owns the per-nnp cell-list cache (positions, image pool, head/next
!> chain, bin geometry), rebuilds or rebins it as atoms move, and walks
!> it to fill the per-element neighbour buffers. Each &NNP force_eval
!> owns its cache, so the lifetime tracks nnp_env_release.
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
!> \author Christoph Schran (christoph.schran@rub.de)
!> \date 2026-05-21
! **************************************************************************************************
MODULE nnp_cell_list
USE cell_types, ONLY: cell_type,&
pbc,&
real_to_scaled
USE kinds, ONLY: dp
USE nnp_environment_types, ONLY: nnp_cell_list_cache_type,&
nnp_env_get,&
nnp_neighbor_type,&
nnp_type
USE nnp_neighbor_interface, ONLY: nnp_neigh_grp_grow
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
! Verlet skin: defaults to MIN(0.5 bohr, 0.1*cutoff), or the VERLET_SKIN
! keyword value (nnp%verlet_skin). The chain rebuilds once an atom drifts
! more than skin/2 from its rebuild-time position.
REAL(KIND=dp), PRIVATE, PARAMETER :: default_verlet_skin = 0.5_dp
REAL(KIND=dp), PRIVATE, PARAMETER :: rebuild_eps = 1.0E-12_dp
PUBLIC :: nnp_prepare_cell_list_cache, &
nnp_compute_neighbors_cell_list
CONTAINS
! **************************************************************************************************
!> \brief Prepare the cell-list cache for the current force evaluation: wrap
!> positions into the primary cell, then either rebuild the bin geometry
!> (cell, atom-count, cutoff or large-drift change) or rebin the chain in
!> place while the Verlet skin still holds. Call once per force evaluation,
!> before any nnp_compute_neighbors_cell_list call.
!> \param nnp NNP environment whose cell_list_cache is prepared from current positions.
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
SUBROUTINE nnp_prepare_cell_list_cache(nnp)
TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
INTEGER :: handle, i
LOGICAL :: rebuild
REAL(KIND=dp) :: dr2_max, exact_cutoff, list_cutoff, skin
REAL(KIND=dp), DIMENSION(3) :: dr
TYPE(cell_type), POINTER :: cell
CALL timeset('nnp_acsf_cell_list_prepare', handle)
NULLIFY (cell)
CALL nnp_env_get(nnp_env=nnp, cell=cell)
exact_cutoff = nnp%max_cut
IF (exact_cutoff > 0.0_dp) THEN
IF (nnp%verlet_skin >= 0.0_dp) THEN
! VERLET_SKIN keyword value (bohr).
skin = nnp%verlet_skin
ELSE
! Auto heuristic: cap at the default, scale down for small cutoffs.
skin = MIN(default_verlet_skin, 0.1_dp*exact_cutoff)
END IF
ELSE
skin = 0.0_dp
END IF
list_cutoff = exact_cutoff + skin
CALL nnp_cell_list_ensure_coord_buffers(nnp%cell_list_cache, nnp%num_atoms)
ASSOCIATE (cache => nnp%cell_list_cache)
DO i = 1, nnp%num_atoms
cache%coord_primary(:, i) = pbc(nnp%coord(:, i), cell, .TRUE.)
CALL real_to_scaled(cache%coord_scaled(:, i), &
cache%coord_primary(:, i), cell)
END DO
rebuild = .NOT. cache%initialized
IF (.NOT. rebuild) rebuild = (cache%num_atoms /= nnp%num_atoms)
IF (.NOT. rebuild) rebuild = (ABS(cache%exact_cutoff - exact_cutoff) > rebuild_eps)
IF (.NOT. rebuild) rebuild = (ABS(cache%list_cutoff - list_cutoff) > rebuild_eps)
IF (.NOT. rebuild) rebuild = (ABS(cache%verlet_skin - skin) > rebuild_eps)
IF (.NOT. rebuild) rebuild = ANY(cache%perd /= cell%perd)
IF (.NOT. rebuild) rebuild = (cache%orthorhombic .NEQV. cell%orthorhombic)
IF (.NOT. rebuild) rebuild = ANY(ABS(cache%hmat - cell%hmat) > rebuild_eps)
IF (.NOT. rebuild) rebuild = ANY(ABS(cache%h_inv - cell%h_inv) > rebuild_eps)
IF ((.NOT. rebuild) .AND. (skin > 0.0_dp)) THEN
dr2_max = 0.0_dp
DO i = 1, nnp%num_atoms
dr(:) = cache%coord_primary(:, i) - cache%ref_coord_primary(:, i)
dr(:) = pbc(dr, cell)
dr2_max = MAX(dr2_max, DOT_PRODUCT(dr, dr))
END DO
IF (dr2_max > 0.25_dp*skin**2) rebuild = .TRUE.
END IF
IF (rebuild) THEN
CALL nnp_build_cell_list_cache(cache, cell, exact_cutoff, list_cutoff, skin)
ELSE IF (cache%initialized) THEN
! On reuse the bin geometry is frozen, but the head/next chain
! is rebinned from current positions so the walk still finds every
! in-cutoff pair.
CALL nnp_cell_list_rebin_chain(cache)
END IF
END ASSOCIATE
CALL timestop(handle)
END SUBROUTINE nnp_prepare_cell_list_cache
! **************************************************************************************************
!> \brief Ensure the persistent coord buffers match num_atoms.
!> Only (re-)allocates on first call or when num_atoms changes.
!> \param cache cell-list cache whose coordinate buffers will be (re-)sized
!> \param num_atoms number of primary-cell atoms the buffers must hold
! **************************************************************************************************
SUBROUTINE nnp_cell_list_ensure_coord_buffers(cache, num_atoms)
TYPE(nnp_cell_list_cache_type), INTENT(INOUT) :: cache
INTEGER, INTENT(IN) :: num_atoms
IF (ALLOCATED(cache%coord_primary)) THEN
IF (SIZE(cache%coord_primary, 2) == num_atoms) RETURN
DEALLOCATE (cache%coord_primary)
END IF
IF (ALLOCATED(cache%coord_scaled)) DEALLOCATE (cache%coord_scaled)
IF (ALLOCATED(cache%ref_coord_primary)) DEALLOCATE (cache%ref_coord_primary)
ALLOCATE (cache%coord_primary(3, num_atoms))
ALLOCATE (cache%coord_scaled(3, num_atoms))
ALLOCATE (cache%ref_coord_primary(3, num_atoms))
! Resizing implies stale bin geometry; force the next prepare to rebuild.
cache%initialized = .FALSE.
cache%num_atoms = num_atoms
END SUBROUTINE nnp_cell_list_ensure_coord_buffers
! **************************************************************************************************
!> \brief Build the Cartesian linked-cell cache from the current positions.
!>
!> Positions must already live in cache%coord_primary / coord_scaled --
!> this routine reads them in place and resets ref_coord_primary for
!> the Verlet skin check.
!> \param cache cell-list cache to populate; coord_primary / coord_scaled must already be filled
!> \param cell simulation cell providing hmat, h_inv and per-axis periodicity
!> \param exact_cutoff largest per-pair SF cutoff used to size the exact image ring
!> \param list_cutoff cell-list bin width (exact_cutoff + Verlet skin)
!> \param skin Verlet skin allowance used by the bin-reuse displacement check
! **************************************************************************************************
SUBROUTINE nnp_build_cell_list_cache(cache, cell, exact_cutoff, list_cutoff, skin)
TYPE(nnp_cell_list_cache_type), INTENT(INOUT) :: cache
TYPE(cell_type), INTENT(IN), POINTER :: cell
REAL(KIND=dp), INTENT(IN) :: exact_cutoff, list_cutoff, skin
INTEGER :: i, img, n_images, nx, ny, nz, tx, ty, tz
INTEGER, DIMENSION(3) :: bin_index, shift
REAL(KIND=dp) :: padding, safe_cutoff
REAL(KIND=dp), DIMENSION(3) :: max_ref, min_ref, range_xyz, ref_pos
CALL nnp_cell_list_release_cell_structure(cache)
cache%exact_cutoff = exact_cutoff
cache%list_cutoff = list_cutoff
cache%verlet_skin = skin
cache%perd = cell%perd
cache%orthorhombic = cell%orthorhombic
cache%hmat = cell%hmat
cache%h_inv = cell%h_inv
cache%ref_coord_primary(:, :) = cache%coord_primary(:, :)
CALL nnp_compute_pbc_copies(cache%exact_pbc_copies, cell, exact_cutoff)
CALL nnp_compute_pbc_copies(cache%list_pbc_copies, cell, list_cutoff)
! The bin stencil clips out-of-range bins rather than wrapping, so periodic
! neighbours come only from pre-replicated images. Force at least one ring
! on each periodic axis (MAX, not +) so boxes with L > 2*cutoff keep their
! wrap-around images.
cache%image_copies = MAX(cache%list_pbc_copies, cell%perd)
! The exact image ring must fit inside the replicated pool; otherwise the
! matcher would accept a ghost that was never created and drop a real pair.
CPASSERT(ALL(cache%exact_pbc_copies <= cache%image_copies))
nx = 2*cache%image_copies(1) + 1
ny = 2*cache%image_copies(2) + 1
nz = 2*cache%image_copies(3) + 1
n_images = cache%num_atoms*nx*ny*nz
cache%n_images = n_images
ALLOCATE (cache%image_atom(n_images))
ALLOCATE (cache%image_shift(3, n_images))
ALLOCATE (cache%image_translation(3, n_images))
min_ref(:) = HUGE(1.0_dp)
max_ref(:) = -HUGE(1.0_dp)
img = 0
DO i = 1, cache%num_atoms
DO tx = -cache%image_copies(1), cache%image_copies(1)
DO ty = -cache%image_copies(2), cache%image_copies(2)
DO tz = -cache%image_copies(3), cache%image_copies(3)
img = img + 1
shift = [tx, ty, tz]
cache%image_atom(img) = i
cache%image_shift(:, img) = shift(:)
cache%image_translation(:, img) = MATMUL(cell%hmat, REAL(shift, KIND=dp))
ref_pos(:) = cache%coord_primary(:, i) + cache%image_translation(:, img)
min_ref(:) = MIN(min_ref(:), ref_pos(:))
max_ref(:) = MAX(max_ref(:), ref_pos(:))
END DO
END DO
END DO
END DO
padding = 0.5_dp*skin + rebuild_eps
cache%lower(:) = min_ref(:) - padding
cache%upper(:) = max_ref(:) + padding
range_xyz(:) = cache%upper(:) - cache%lower(:)
! Coarse grid: bin width = list_cutoff (bin_span 1, 27-bin walk). A finer
! grid changes the order in which neighbours are appended, which reorders
! the descriptor summation, so it is not used.
safe_cutoff = MAX(list_cutoff, rebuild_eps)
DO i = 1, 3
IF (range_xyz(i) > rebuild_eps) THEN
cache%nbin(i) = MAX(1, CEILING(range_xyz(i)/safe_cutoff))
cache%bin_width(i) = range_xyz(i)/REAL(cache%nbin(i), KIND=dp)
cache%bin_span(i) = MIN(cache%nbin(i) - 1, &
CEILING(list_cutoff/MAX(cache%bin_width(i), rebuild_eps)))
ELSE
cache%nbin(i) = 1
cache%bin_width(i) = 1.0_dp
cache%bin_span(i) = 0
END IF
END DO
cache%n_cells = PRODUCT(cache%nbin)
ALLOCATE (cache%head(cache%n_cells))
ALLOCATE (cache%next(n_images))
cache%head(:) = 0
cache%next(:) = 0
DO img = 1, n_images
i = cache%image_atom(img)
ref_pos(:) = cache%coord_primary(:, i) + cache%image_translation(:, img)
CALL nnp_cell_list_bin_from_position(cache, ref_pos, bin_index)
cache%next(img) = cache%head(nnp_cell_list_linear_index(cache, bin_index))
cache%head(nnp_cell_list_linear_index(cache, bin_index)) = img
END DO
cache%initialized = .TRUE.
END SUBROUTINE nnp_build_cell_list_cache
! **************************************************************************************************
!> \brief Rebin the head/next chain from current coord_primary while keeping the
!> image pool, bin geometry and ref_coord_primary frozen. Used on Verlet
!> skin reuse so the chain reflects atoms that crossed bin boundaries.
!> \param cache cell-list cache whose head/next chain is rebuilt from current coord_primary
! **************************************************************************************************
SUBROUTINE nnp_cell_list_rebin_chain(cache)
TYPE(nnp_cell_list_cache_type), INTENT(INOUT) :: cache
INTEGER :: i, img, lin_idx
INTEGER, DIMENSION(3) :: bin_index
REAL(KIND=dp), DIMENSION(3) :: ref_pos
cache%head(:) = 0
cache%next(:) = 0
DO img = 1, cache%n_images
i = cache%image_atom(img)
ref_pos(:) = cache%coord_primary(:, i) + cache%image_translation(:, img)
CALL nnp_cell_list_bin_from_position(cache, ref_pos, bin_index)
lin_idx = nnp_cell_list_linear_index(cache, bin_index)
cache%next(img) = cache%head(lin_idx)
cache%head(lin_idx) = img
END DO
END SUBROUTINE nnp_cell_list_rebin_chain
! **************************************************************************************************
!> \brief Free only the cell-structure pieces of the cache (bins and image
!> metadata). The coordinate buffers are persistent and kept, so a
!> rebuild can still read the current positions.
!> \param cache cell-list cache whose bin geometry and image metadata are freed (coord buffers kept)
! **************************************************************************************************
SUBROUTINE nnp_cell_list_release_cell_structure(cache)
TYPE(nnp_cell_list_cache_type), INTENT(INOUT) :: cache
IF (ALLOCATED(cache%image_atom)) DEALLOCATE (cache%image_atom)
IF (ALLOCATED(cache%head)) DEALLOCATE (cache%head)
IF (ALLOCATED(cache%next)) DEALLOCATE (cache%next)
IF (ALLOCATED(cache%image_shift)) DEALLOCATE (cache%image_shift)
IF (ALLOCATED(cache%image_translation)) DEALLOCATE (cache%image_translation)
cache%n_images = 0
cache%n_cells = 1
cache%exact_pbc_copies = 0
cache%list_pbc_copies = 0
cache%image_copies = 0
cache%nbin = 1
cache%bin_span = 0
cache%lower = 0.0_dp
cache%upper = 0.0_dp
cache%bin_width = 1.0_dp
END SUBROUTINE nnp_cell_list_release_cell_structure
! **************************************************************************************************
!> \brief Fill the ACSF neighbour buffers for one central atom by walking the
!> linked-cell stencil around its bin, pushing (j, dr, r) into the
!> per-species-pair slabs for entries inside each pair cutoff. The own-bin
!> self image (i == j, zero shift) is skipped.
!> \param nnp NNP environment with cell_list_cache and neighbor_interface_state ready
!> \param neighbor per-atom neighbour view filled in place (counters zeroed before this call)
!> \param i central-atom index in the global ordering
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
!> \author Christoph Schran (christoph.schran@rub.de)
! **************************************************************************************************
SUBROUTINE nnp_compute_neighbors_cell_list(nnp, neighbor, i)
TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
TYPE(nnp_neighbor_type), INTENT(INOUT) :: neighbor
INTEGER, INTENT(IN) :: i
INTEGER :: bx, by, bz, img, ind, j, neighbor_ind, &
pair_slot, s
INTEGER, DIMENSION(3) :: bin_index, current_shift
REAL(KIND=dp) :: norm
REAL(KIND=dp), DIMENSION(3) :: center, dr, image_pos
ASSOCIATE (cache => nnp%cell_list_cache, &
state => nnp%neighbor_interface_state)
ind = nnp%ele_ind(i)
center(:) = cache%coord_primary(:, i)
CALL nnp_cell_list_bin_from_position(cache, center, bin_index)
DO bx = MAX(1, bin_index(1) - cache%bin_span(1)), &
MIN(cache%nbin(1), bin_index(1) + cache%bin_span(1))
DO by = MAX(1, bin_index(2) - cache%bin_span(2)), &
MIN(cache%nbin(2), bin_index(2) + cache%bin_span(2))
DO bz = MAX(1, bin_index(3) - cache%bin_span(3)), &
MIN(cache%nbin(3), bin_index(3) + cache%bin_span(3))
img = cache%head(nnp_cell_list_linear_index(cache, [bx, by, bz]))
DO WHILE (img > 0)
j = cache%image_atom(img)
current_shift(:) = cache%image_shift(:, img)
IF (j == i .AND. ALL(current_shift == 0)) THEN
img = cache%next(img)
CYCLE
END IF
IF (.NOT. nnp_cell_list_exact_image_match(cache, i, j, current_shift)) THEN
img = cache%next(img)
CYCLE
END IF
image_pos(:) = cache%coord_primary(:, j) + cache%image_translation(:, img)
dr(:) = center(:) - image_pos(:)
norm = NORM2(dr(:))
neighbor_ind = nnp%ele_ind(j)
IF (norm < state%pair_map(ind, neighbor_ind)%max_relevant_cutoff) THEN
DO pair_slot = 1, state%pair_map(ind, neighbor_ind)%n_rad
s = state%pair_map(ind, neighbor_ind)%rad_groups(pair_slot)
IF (norm < nnp%rad(ind)%symfgrp(s)%cutoff) THEN
neighbor%n_rad(s) = neighbor%n_rad(s) + 1
IF (neighbor%n_rad(s) > neighbor%rad(s)%cap) &
CALL nnp_neigh_grp_grow(neighbor%rad(s), neighbor%n_rad(s))
neighbor%rad(s)%ind(neighbor%n_rad(s)) = j
neighbor%rad(s)%dist(1, neighbor%n_rad(s)) = dr(1)
neighbor%rad(s)%dist(2, neighbor%n_rad(s)) = dr(2)
neighbor%rad(s)%dist(3, neighbor%n_rad(s)) = dr(3)
neighbor%rad(s)%dist(4, neighbor%n_rad(s)) = norm
END IF
END DO
DO pair_slot = 1, state%pair_map(ind, neighbor_ind)%n_ang1
s = state%pair_map(ind, neighbor_ind)%ang1_groups(pair_slot)
IF (norm < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
neighbor%n_ang1(s) = neighbor%n_ang1(s) + 1
IF (neighbor%n_ang1(s) > neighbor%ang1(s)%cap) &
CALL nnp_neigh_grp_grow(neighbor%ang1(s), neighbor%n_ang1(s))
neighbor%ang1(s)%ind(neighbor%n_ang1(s)) = j
neighbor%ang1(s)%dist(1, neighbor%n_ang1(s)) = dr(1)
neighbor%ang1(s)%dist(2, neighbor%n_ang1(s)) = dr(2)
neighbor%ang1(s)%dist(3, neighbor%n_ang1(s)) = dr(3)
neighbor%ang1(s)%dist(4, neighbor%n_ang1(s)) = norm
END IF
END DO
DO pair_slot = 1, state%pair_map(ind, neighbor_ind)%n_ang2
s = state%pair_map(ind, neighbor_ind)%ang2_groups(pair_slot)
IF (norm < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
neighbor%n_ang2(s) = neighbor%n_ang2(s) + 1
IF (neighbor%n_ang2(s) > neighbor%ang2(s)%cap) &
CALL nnp_neigh_grp_grow(neighbor%ang2(s), neighbor%n_ang2(s))
neighbor%ang2(s)%ind(neighbor%n_ang2(s)) = j
neighbor%ang2(s)%dist(1, neighbor%n_ang2(s)) = dr(1)
neighbor%ang2(s)%dist(2, neighbor%n_ang2(s)) = dr(2)
neighbor%ang2(s)%dist(3, neighbor%n_ang2(s)) = dr(3)
neighbor%ang2(s)%dist(4, neighbor%n_ang2(s)) = norm
END IF
END DO
END IF
img = cache%next(img)
END DO
END DO
END DO
END DO
END ASSOCIATE
END SUBROUTINE nnp_compute_neighbors_cell_list
! **************************************************************************************************
!> \brief Check whether a ghost image corresponds to one of the exact pbc copies.
!> \param cache cell-list cache providing per-axis periodicity and scaled coordinates
!> \param i central-atom index in the primary cell
!> \param j candidate neighbour atom index in the primary cell
!> \param image_shift 3-vector of integer cell translations identifying the ghost image
!> \return .TRUE. if the (i, j, image_shift) image falls inside the exact_pbc_copies ring
! **************************************************************************************************
LOGICAL FUNCTION nnp_cell_list_exact_image_match(cache, i, j, image_shift)
TYPE(nnp_cell_list_cache_type), INTENT(IN) :: cache
INTEGER, INTENT(IN) :: i, j
INTEGER, DIMENSION(3), INTENT(IN) :: image_shift
INTEGER :: d, minimal_shift
nnp_cell_list_exact_image_match = .TRUE.
DO d = 1, 3
IF (cache%perd(d) == 1) THEN
minimal_shift = NINT(cache%coord_scaled(d, i) - cache%coord_scaled(d, j))
IF (ABS(minimal_shift - image_shift(d)) > cache%exact_pbc_copies(d)) THEN
nnp_cell_list_exact_image_match = .FALSE.
RETURN
END IF
ELSE IF (image_shift(d) /= 0) THEN
nnp_cell_list_exact_image_match = .FALSE.
RETURN
END IF
END DO
END FUNCTION nnp_cell_list_exact_image_match
! **************************************************************************************************
!> \brief Convert Cartesian position to linked-cell bin coordinates.
!> \param cache cell-list cache providing the bin origin and bin widths
!> \param pos Cartesian position (Bohr) to bin
!> \param bin_index (out) 1-based bin indices clipped to [1, nbin] on each axis
! **************************************************************************************************
SUBROUTINE nnp_cell_list_bin_from_position(cache, pos, bin_index)
TYPE(nnp_cell_list_cache_type), INTENT(IN) :: cache
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: pos
INTEGER, DIMENSION(3), INTENT(OUT) :: bin_index
INTEGER :: d
REAL(KIND=dp) :: scaled
DO d = 1, 3
scaled = (pos(d) - cache%lower(d))/MAX(cache%bin_width(d), rebuild_eps)
bin_index(d) = INT(scaled) + 1
bin_index(d) = MAX(1, MIN(cache%nbin(d), bin_index(d)))
END DO
END SUBROUTINE nnp_cell_list_bin_from_position
! **************************************************************************************************
!> \brief Flatten 3D bin indices into the head/next storage index.
!> \param cache cell-list cache providing the per-axis bin counts
!> \param bin_index 1-based 3D bin indices to flatten
!> \return linear index into cache%head
! **************************************************************************************************
INTEGER FUNCTION nnp_cell_list_linear_index(cache, bin_index)
TYPE(nnp_cell_list_cache_type), INTENT(IN) :: cache
INTEGER, DIMENSION(3), INTENT(IN) :: bin_index
nnp_cell_list_linear_index = bin_index(1) + cache%nbin(1)* &
((bin_index(2) - 1) + cache%nbin(2)*(bin_index(3) - 1))
END FUNCTION nnp_cell_list_linear_index
! **************************************************************************************************
!> \brief Number of PBC image rings along each lattice axis needed to cover a
!> sphere of radius `cutoff` from any primary-cell atom. Uses the
!> cell-vector projections onto the plane normals (correct for triclinic
!> cells) and multiplies by cell%perd so non-periodic axes return 0.
!> \param pbc_copies output: number of image rings on each of the 3 lattice axes
!> \param cell cell whose hmat/h_inv/perd/deth are read
!> \param cutoff sphere radius (Bohr) to cover
!> \author Christoph Schran (christoph.schran@rub.de)
! **************************************************************************************************
SUBROUTINE nnp_compute_pbc_copies(pbc_copies, cell, cutoff)
INTEGER, DIMENSION(3), INTENT(OUT) :: pbc_copies
TYPE(cell_type), INTENT(IN), POINTER :: cell
REAL(KIND=dp), INTENT(IN) :: cutoff
REAL(KIND=dp) :: proja, projb, projc
REAL(KIND=dp), DIMENSION(3) :: axb, axc, bxc
! A degenerate cell (deth <= 0) makes the NORM2 normalisations below
! divide by zero and the DO WHILE projection loops never terminate.
CPASSERT(cell%deth > 0.0_dp)
axb(1) = cell%hmat(2, 1)*cell%hmat(3, 2) - cell%hmat(3, 1)*cell%hmat(2, 2)
axb(2) = cell%hmat(3, 1)*cell%hmat(1, 2) - cell%hmat(1, 1)*cell%hmat(3, 2)
axb(3) = cell%hmat(1, 1)*cell%hmat(2, 2) - cell%hmat(2, 1)*cell%hmat(1, 2)
axb(:) = axb(:)/NORM2(axb(:))
axc(1) = cell%hmat(2, 1)*cell%hmat(3, 3) - cell%hmat(3, 1)*cell%hmat(2, 3)
axc(2) = cell%hmat(3, 1)*cell%hmat(1, 3) - cell%hmat(1, 1)*cell%hmat(3, 3)
axc(3) = cell%hmat(1, 1)*cell%hmat(2, 3) - cell%hmat(2, 1)*cell%hmat(1, 3)
axc(:) = axc(:)/NORM2(axc(:))
bxc(1) = cell%hmat(2, 2)*cell%hmat(3, 3) - cell%hmat(3, 2)*cell%hmat(2, 3)
bxc(2) = cell%hmat(3, 2)*cell%hmat(1, 3) - cell%hmat(1, 2)*cell%hmat(3, 3)
bxc(3) = cell%hmat(1, 2)*cell%hmat(2, 3) - cell%hmat(2, 2)*cell%hmat(1, 3)
bxc(:) = bxc(:)/NORM2(bxc(:))
proja = ABS(SUM(cell%hmat(:, 1)*bxc(:)))*0.5_dp
projb = ABS(SUM(cell%hmat(:, 2)*axc(:)))*0.5_dp
projc = ABS(SUM(cell%hmat(:, 3)*axb(:)))*0.5_dp
pbc_copies(:) = 0
DO WHILE ((pbc_copies(1) + 1)*proja <= cutoff)
pbc_copies(1) = pbc_copies(1) + 1
END DO
DO WHILE ((pbc_copies(2) + 1)*projb <= cutoff)
pbc_copies(2) = pbc_copies(2) + 1
END DO
DO WHILE ((pbc_copies(3) + 1)*projc <= cutoff)
pbc_copies(3) = pbc_copies(3) + 1
END DO
pbc_copies(:) = pbc_copies(:)*cell%perd(:)
END SUBROUTINE nnp_compute_pbc_copies
END MODULE nnp_cell_list

214
src/nnp_environment.F
View file

@ -8,6 +8,7 @@
! **************************************************************************************************
!> \brief Methods dealing with Neural Network potentials
!> \author Christoph Schran (christoph.schran@rub.de)
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
!> \date 2020-10-10
! **************************************************************************************************
MODULE nnp_environment
@ -63,7 +64,7 @@ MODULE nnp_environment
PRIVATE
LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_environment'
PUBLIC :: nnp_init
@ -231,20 +232,16 @@ CONTAINS
CHARACTER(LEN=*), INTENT(IN) :: printtag
CHARACTER(len=*), PARAMETER :: routineN = 'nnp_init_model'
INTEGER, PARAMETER :: def_str_len = 256, &
default_path_length = 256
INTEGER, PARAMETER :: def_str_len = 256
CHARACTER(len=1), ALLOCATABLE, DIMENSION(:) :: cactfnct
CHARACTER(len=2) :: ele
CHARACTER(len=def_str_len) :: dummy, line
CHARACTER(len=default_path_length) :: base_name, file_name
INTEGER :: handle, i, i_com, io, iweight, j, k, l, &
n_weight, nele, nuc_ele, symfnct_type, &
unit_nr
LOGICAL :: at_end, atom_e_found, explicit, first, &
found
CHARACTER(len=default_path_length) :: file_name
INTEGER :: handle, i, io, j, k, nele, nuc_ele, &
symfnct_type, unit_nr
LOGICAL :: atom_e_found, explicit, first, found
REAL(KIND=dp) :: energy
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: weights
REAL(KIND=dp), DIMENSION(7) :: test_array
REAL(KIND=dp), DIMENSION(:), POINTER :: work
TYPE(cp_logger_type), POINTER :: logger
@ -257,8 +254,8 @@ CONTAINS
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
unit_nr = cp_logger_get_default_unit_nr(logger)
unit_nr = cp_logger_get_default_unit_nr(logger)
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) ""
WRITE (unit_nr, *) TRIM(printtag)//"| Neural Network Potential Force Environment"
END IF
@ -272,6 +269,11 @@ CONTAINS
ALLOCATE (nnp_env%committee_stress(3, 3, nnp_env%n_committee))
CALL section_vals_val_get(nnp_env%nnp_input, "NNP_INPUT_FILE_NAME", c_val=file_name)
CALL section_vals_val_get(nnp_env%nnp_input, "RAD_SPLINE_N", i_val=nnp_env%rad_spline_n)
IF (nnp_env%rad_spline_n < 2) THEN
CPABORT("NNP| RAD_SPLINE_N must be >= 2.")
END IF
CALL section_vals_val_get(nnp_env%nnp_input, "VERLET_SKIN", r_val=nnp_env%verlet_skin)
CALL parser_create(parser, file_name, para_env=logger%para_env)
! read number of elements and cut_type and check for scale and center
@ -284,8 +286,7 @@ CONTAINS
nnp_env%normnodes = .FALSE.
nnp_env%n_hlayer = 0
IF (logger%para_env%is_source()) THEN
unit_nr = cp_logger_get_default_unit_nr(logger)
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) TRIM(printtag)//"| Reading NNP input from file: ", TRIM(file_name)
END IF
@ -306,10 +307,12 @@ CONTAINS
search_from_begin_of_file=.TRUE.)
nnp_env%scale_acsf = found
! Test if there are two keywords of this:
! parser_search_string matches substrings, so a search for
! "scale_symmetry_functions" also fires on "..._sigma". Re-search from the
! current parser position to confirm the bare keyword appears on its own line.
CALL parser_search_string(parser, "scale_symmetry_functions", .TRUE., found)
IF (found .AND. nnp_env%scale_sigma_acsf) THEN
CPWARN('Two scaling keywords in the input, we will ignore sigma scaling in this case')
CPWARN("Ignoring sigma ACSF scaling; both keywords were set.")
nnp_env%scale_sigma_acsf = .FALSE.
ELSE IF (.NOT. found .AND. nnp_env%scale_sigma_acsf) THEN
nnp_env%scale_acsf = .FALSE.
@ -634,7 +637,7 @@ CONTAINS
! read scaling information from file
IF (nnp_env%scale_acsf .OR. nnp_env%center_acsf .OR. nnp_env%scale_sigma_acsf) THEN
IF (logger%para_env%is_source()) THEN
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) TRIM(printtag)//"| Reading scaling information from file: ", TRIM(file_name)
END IF
CALL section_vals_val_get(nnp_env%nnp_input, "SCALE_FILE_NAME", &
@ -689,11 +692,59 @@ CONTAINS
END DO
END DO
CALL parser_release(parser)
! Reject degenerate scaling.data at load time. The (loc_max - loc_min)
! and sigma denominators are used unguarded inside nnp_scale_acsf;
! a zero range or zero sigma there produces silent NaN forces.
IF (nnp_env%scale_acsf) THEN
DO i = 1, nnp_env%n_ele
DO j = 1, nnp_env%n_rad(i)
IF (nnp_env%rad(i)%loc_max(j) <= nnp_env%rad(i)%loc_min(j)) THEN
WRITE (line, '(A,I0,A,I0,A,2(1X,ES13.6))') &
"scaling.data: radial sf range non-positive for element ", i, &
" sf ", j, " (loc_min, loc_max) =", &
nnp_env%rad(i)%loc_min(j), nnp_env%rad(i)%loc_max(j)
CPABORT(TRIM(line))
END IF
END DO
DO j = 1, nnp_env%n_ang(i)
IF (nnp_env%ang(i)%loc_max(j) <= nnp_env%ang(i)%loc_min(j)) THEN
WRITE (line, '(A,I0,A,I0,A,2(1X,ES13.6))') &
"scaling.data: angular sf range non-positive for element ", i, &
" sf ", j, " (loc_min, loc_max) =", &
nnp_env%ang(i)%loc_min(j), nnp_env%ang(i)%loc_max(j)
CPABORT(TRIM(line))
END IF
END DO
END DO
END IF
IF (nnp_env%scale_sigma_acsf) THEN
DO i = 1, nnp_env%n_ele
DO j = 1, nnp_env%n_rad(i)
IF (nnp_env%rad(i)%sigma(j) <= 0.0_dp) THEN
WRITE (line, '(A,I0,A,I0,A,1X,ES13.6)') &
"scaling.data: radial sf sigma non-positive for element ", i, &
" sf ", j, " sigma =", nnp_env%rad(i)%sigma(j)
CPABORT(TRIM(line))
END IF
END DO
DO j = 1, nnp_env%n_ang(i)
IF (nnp_env%ang(i)%sigma(j) <= 0.0_dp) THEN
WRITE (line, '(A,I0,A,I0,A,1X,ES13.6)') &
"scaling.data: angular sf sigma non-positive for element ", i, &
" sf ", j, " sigma =", nnp_env%ang(i)%sigma(j)
CPABORT(TRIM(line))
END IF
END DO
END DO
END IF
END IF
CALL nnp_init_acsf_groups(nnp_env)
! read weights from file
! Allocate per-layer weight tables, then load each committee member
! from disk. Allocation lives at the call site because n_committee /
! n_layer / n_nodes are owned here; the loader only fills them.
DO i = 1, nnp_env%n_ele
DO j = 2, nnp_env%n_layer
ALLOCATE (nnp_env%arc(i)%layer(j)%weights(nnp_env%arc(i)%n_nodes(j - 1), &
@ -701,14 +752,90 @@ CONTAINS
ALLOCATE (nnp_env%arc(i)%layer(j)%bweights(nnp_env%arc(i)%n_nodes(j), nnp_env%n_committee))
END DO
END DO
CALL nnp_read_committee_weights(nnp_env, model_section, printtag)
nnp_env%expol = 0
! Bias the standard deviation of committee disagreement
NULLIFY (bias_section)
explicit = .FALSE.
!HELIUM NNP does not currently define a bias term
bias_section => section_vals_get_subs_vals(nnp_env%nnp_input, "BIAS", can_return_null=.TRUE.)
IF (ASSOCIATED(bias_section)) CALL section_vals_get(bias_section, explicit=explicit)
nnp_env%bias = .FALSE.
IF (explicit) THEN
IF (nnp_env%n_committee > 1) THEN
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) "NNP| Biasing of committee disagreement enabled"
END IF
nnp_env%bias = .TRUE.
ALLOCATE (nnp_env%bias_forces(3, nnp_env%num_atoms))
ALLOCATE (nnp_env%bias_e_avrg(nnp_env%n_committee))
CALL section_vals_val_get(bias_section, "SIGMA_0", r_val=nnp_env%bias_sigma0)
CALL section_vals_val_get(bias_section, "K_B", r_val=nnp_env%bias_kb)
nnp_env%bias_e_avrg(:) = 0.0_dp
CALL section_vals_val_get(bias_section, "ALIGN_NNP_ENERGIES", explicit=explicit)
nnp_env%bias_align = explicit
IF (explicit) THEN
NULLIFY (work)
CALL section_vals_val_get(bias_section, "ALIGN_NNP_ENERGIES", r_vals=work)
IF (SIZE(work) /= nnp_env%n_committee) THEN
CPABORT("ALIGN_NNP_ENERGIES size mismatch wrt committee size.")
END IF
nnp_env%bias_e_avrg(:) = work
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) TRIM(printtag)//"| Biasing is aligned by shifting the energy prediction of the C-NNP members"
END IF
END IF
ELSE
CPWARN("NNP committee size is 1, BIAS section is ignored.")
END IF
END IF
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) TRIM(printtag)//"| NNP force environment initialized"
END IF
CALL timestop(handle)
END SUBROUTINE nnp_init_model
! **************************************************************************************************
!> \brief Read committee weights from disk into pre-allocated arc layer tables.
!>
!> One file per (committee member, element), named <base>.<nuc>.data. The file
!> is a flat list of doubles; the order on disk is layer-major over (i->j weights
!> then biases) so a single counter walks every value into its slot.
!>
!> \param nnp_env NNP environment with arc()%layer()%weights/bweights already allocated.
!> \param model_section &MODEL section block; supplies the per-member WEIGHTS base path.
!> \param printtag Log-line prefix.
! **************************************************************************************************
SUBROUTINE nnp_read_committee_weights(nnp_env, model_section, printtag)
TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp_env
TYPE(section_vals_type), POINTER :: model_section
CHARACTER(LEN=*), INTENT(IN) :: printtag
CHARACTER(len=default_path_length) :: base_name, file_name
INTEGER :: i, i_com, iweight, j, k, l, n_weight, &
unit_nr
LOGICAL :: at_end
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: weights
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_parser_type) :: parser
NULLIFY (logger)
logger => cp_get_default_logger()
unit_nr = cp_logger_get_default_unit_nr(logger)
DO i_com = 1, nnp_env%n_committee
CALL section_vals_val_get(model_section, "WEIGHTS", c_val=base_name, i_rep_section=i_com)
IF (logger%para_env%is_source()) THEN
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) TRIM(printtag)//"| Initializing weights for model: ", i_com
END IF
DO i = 1, nnp_env%n_ele
WRITE (file_name, '(A,I0.3,A)') TRIM(base_name)//".", nnp_env%nuc_ele(i), ".data"
IF (logger%para_env%is_source()) THEN
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) TRIM(printtag)//"| Reading weights from file: ", TRIM(file_name)
END IF
CALL parser_create(parser, file_name, para_env=logger%para_env)
@ -748,51 +875,6 @@ CONTAINS
END DO
END DO
!Initialize extrapolation counter
nnp_env%expol = 0
! Bias the standard deviation of committee disagreement
NULLIFY (bias_section)
explicit = .FALSE.
!HELIUM NNP does atm not allow for bias (not even defined)
bias_section => section_vals_get_subs_vals(nnp_env%nnp_input, "BIAS", can_return_null=.TRUE.)
IF (ASSOCIATED(bias_section)) CALL section_vals_get(bias_section, explicit=explicit)
nnp_env%bias = .FALSE.
IF (explicit) THEN
IF (nnp_env%n_committee > 1) THEN
IF (logger%para_env%is_source()) THEN
WRITE (unit_nr, *) "NNP| Biasing of committee disagreement enabled"
END IF
nnp_env%bias = .TRUE.
ALLOCATE (nnp_env%bias_forces(3, nnp_env%num_atoms))
ALLOCATE (nnp_env%bias_e_avrg(nnp_env%n_committee))
CALL section_vals_val_get(bias_section, "SIGMA_0", r_val=nnp_env%bias_sigma0)
CALL section_vals_val_get(bias_section, "K_B", r_val=nnp_env%bias_kb)
nnp_env%bias_e_avrg(:) = 0.0_dp
CALL section_vals_val_get(bias_section, "ALIGN_NNP_ENERGIES", explicit=explicit)
nnp_env%bias_align = explicit
IF (explicit) THEN
NULLIFY (work)
CALL section_vals_val_get(bias_section, "ALIGN_NNP_ENERGIES", r_vals=work)
IF (SIZE(work) /= nnp_env%n_committee) THEN
CPABORT("ALIGN_NNP_ENERGIES size mismatch wrt committee size.")
END IF
nnp_env%bias_e_avrg(:) = work
IF (logger%para_env%is_source()) THEN
WRITE (unit_nr, *) TRIM(printtag)//"| Biasing is aligned by shifting the energy prediction of the C-NNP members"
END IF
END IF
ELSE
CPWARN("NNP committee size is 1, BIAS section is ignored.")
END IF
END IF
IF (logger%para_env%is_source()) THEN
WRITE (unit_nr, *) TRIM(printtag)//"| NNP force environment initialized"
END IF
CALL timestop(handle)
END SUBROUTINE nnp_init_model
END SUBROUTINE nnp_read_committee_weights
END MODULE nnp_environment

383
src/nnp_environment_types.F
View file

@ -8,6 +8,7 @@
! **************************************************************************************************
!> \brief Data types for neural network potentials
!> \author Christoph Schran (christoph.schran@rub.de)
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
!> \date 2020-10-10
! **************************************************************************************************
MODULE nnp_environment_types
@ -45,20 +46,29 @@ MODULE nnp_environment_types
PRIVATE
LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_environment_types'
!> derived data types
PUBLIC :: nnp_type
PUBLIC :: nnp_arc_type
PUBLIC :: nnp_neighbor_type
PUBLIC :: nnp_neigh_grp_type
PUBLIC :: nnp_symfgrp_type
PUBLIC :: nnp_acsf_rad_type
PUBLIC :: nnp_acsf_ang_type
PUBLIC :: nnp_cell_list_cache_type
PUBLIC :: nnp_neighbor_pair_map_type
PUBLIC :: nnp_dGdr_grp_type
PUBLIC :: nnp_neighbor_workspace_type
PUBLIC :: nnp_neighbor_interface_state_type
! Public subroutines ***
PUBLIC :: nnp_env_release, &
nnp_env_set, &
nnp_env_get
nnp_env_get, &
nnp_cell_list_cache_release, &
nnp_neighbor_interface_state_release
INTEGER, PARAMETER, PUBLIC :: &
nnp_cut_cos = 1, &
@ -105,9 +115,11 @@ MODULE nnp_environment_types
INTEGER :: n_committee = -1
INTEGER :: n_hlayer = -1
INTEGER :: n_layer = -1
INTEGER :: rad_spline_n = 8192 ! knots/radial group (RAD_SPLINE_N)
REAL(KIND=dp) :: verlet_skin = -1.0_dp ! cell-list skin, bohr (VERLET_SKIN; <0 = auto)
INTEGER, ALLOCATABLE, DIMENSION(:) :: n_hnodes
INTEGER, ALLOCATABLE, DIMENSION(:) :: actfnct
INTEGER :: expol = -1 ! extrapolation coutner
INTEGER :: expol = -1 ! extrapolation counter
LOGICAL :: output_expol = .FALSE. ! output extrapolation
! structures for calculation
INTEGER :: num_atoms = -1
@ -136,6 +148,12 @@ MODULE nnp_environment_types
REAL(KIND=dp) :: bias_sigma = -1.0_dp
REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: bias_forces
REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: bias_e_avrg
! Per-nnp persistent neighbour-finder state, owned by nnp_type so committee
! and MIX force_evals carry independent caches whose lifetime tracks
! nnp_env_release. Allocated by nnp_prepare_neighbor_cache, freed in
! nnp_env_release; operated on by nnp_cell_list and nnp_neighbor_interface.
TYPE(nnp_cell_list_cache_type), ALLOCATABLE :: cell_list_cache
TYPE(nnp_neighbor_interface_state_type), ALLOCATABLE :: neighbor_interface_state
END TYPE nnp_type
! **************************************************************************************************
@ -153,6 +171,26 @@ MODULE nnp_environment_types
INTEGER, DIMENSION(:), ALLOCATABLE :: ele_ind
CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: ele
REAL(KIND=dp) :: cutoff = -1.0_dp
! Packed per-member parameters: contiguous arrays sized n_symf so the inner
! SF loop streams memory instead of indexing ang(ind)%{eta,zeta,lam,prefzeta}
! through symf(sf). pack_izeta/pack_use_int_zeta skip NINT in the inner loop.
REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: pack_eta
REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: pack_zeta
REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: pack_lam
REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: pack_prefzeta
INTEGER, DIMENSION(:), ALLOCATABLE :: pack_izeta
LOGICAL, DIMENSION(:), ALLOCATABLE :: pack_use_int_zeta
! Group-shared Hermite cubic spline tables (radial groups only). All SFs in
! a group share grp%cutoff, hence one uniform grid (spline_n nodes, spacing
! spline_dx). y(r) and y'(r) are packed sf-first so nnp_calc_rad streams them
! after computing the interpolation parameters once per call.
LOGICAL :: spline_built = .FALSE.
INTEGER :: spline_n = 0
REAL(KIND=dp) :: spline_dx = 1.0_dp
REAL(KIND=dp) :: spline_dx_inv = 1.0_dp
REAL(KIND=dp) :: spline_x_max = 0.0_dp
REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: spline_y ! (n_symf, n_grid)
REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: spline_dy ! (n_symf, n_grid)
END TYPE nnp_symfgrp_type
! **************************************************************************************************
@ -221,25 +259,163 @@ MODULE nnp_environment_types
END TYPE nnp_acsf_ang_type
! **************************************************************************************************
!> \brief Contains neighbors list of an atom
!> \param dist - distance vectors + norm DIM(4,nat)
!> \param n - number of neighbors
!> \brief Per-SF-group dense neighbour container. cap is the allocated slab size;
!> the live count is nnp_neighbor_type%n_rad/n_ang1/n_ang2 for the matching
!> group. ind(j) is the j-th neighbour's atom index; dist(1:3, j) is its
!> displacement and dist(4, j) its norm. Sized lazily by nnp_neigh_grp_grow.
! **************************************************************************************************
TYPE nnp_neigh_grp_type
INTEGER :: cap = 0
INTEGER, ALLOCATABLE, DIMENSION(:) :: ind ! (cap)
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dist ! (4, cap)
END TYPE nnp_neigh_grp_type
! **************************************************************************************************
!> \brief Contains neighbors list of an atom (per-group dense layout).
!> \param n_rad/n_ang1/n_ang2 - live neighbor counts per symfgrp
!> \param rad/ang1/ang2 - per-group dense (ind, dist) containers
!> \author Christoph Schran (christoph.schran@rub.de)
!> \date 2020-10-10
! **************************************************************************************************
TYPE nnp_neighbor_type
INTEGER, DIMENSION(3) :: pbc_copies = -1
INTEGER, DIMENSION(:), ALLOCATABLE :: n_rad
INTEGER, DIMENSION(:), ALLOCATABLE :: n_ang1
INTEGER, DIMENSION(:), ALLOCATABLE :: n_ang2
INTEGER, DIMENSION(:, :), ALLOCATABLE :: ind_rad
INTEGER, DIMENSION(:, :), ALLOCATABLE :: ind_ang1
INTEGER, DIMENSION(:, :), ALLOCATABLE :: ind_ang2
REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dist_rad
REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dist_ang1
REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dist_ang2
INTEGER, DIMENSION(3) :: pbc_copies = -1
INTEGER, DIMENSION(:), ALLOCATABLE :: n_rad
INTEGER, DIMENSION(:), ALLOCATABLE :: n_ang1
INTEGER, DIMENSION(:), ALLOCATABLE :: n_ang2
TYPE(nnp_neigh_grp_type), ALLOCATABLE, DIMENSION(:) :: rad
TYPE(nnp_neigh_grp_type), ALLOCATABLE, DIMENSION(:) :: ang1
TYPE(nnp_neigh_grp_type), ALLOCATABLE, DIMENSION(:) :: ang2
END TYPE nnp_neighbor_type
! **************************************************************************************************
!> \brief Linked-cell + Verlet-skin cache for the NNP descriptor neighbour walk.
!> Persisted across MD steps so the image pool, bin geometry, and
!> reference positions survive between force evaluations and the
!> head/next chain can be reused as long as no atom has drifted more
!> than skin/2. Operated on by nnp_cell_list.F.
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
TYPE nnp_cell_list_cache_type
LOGICAL :: initialized = .FALSE.
INTEGER :: num_atoms = -1
INTEGER :: n_images = 0
INTEGER :: n_cells = 1
INTEGER, DIMENSION(3) :: exact_pbc_copies = 0
INTEGER, DIMENSION(3) :: list_pbc_copies = 0
INTEGER, DIMENSION(3) :: image_copies = 0
INTEGER, DIMENSION(3) :: nbin = 1
INTEGER, DIMENSION(3) :: bin_span = 0
INTEGER, DIMENSION(3) :: perd = 0
LOGICAL :: orthorhombic = .FALSE.
REAL(KIND=dp) :: exact_cutoff = -1.0_dp
REAL(KIND=dp) :: list_cutoff = -1.0_dp
REAL(KIND=dp) :: verlet_skin = 0.0_dp
REAL(KIND=dp), DIMENSION(3) :: lower = 0.0_dp
REAL(KIND=dp), DIMENSION(3) :: upper = 0.0_dp
REAL(KIND=dp), DIMENSION(3) :: bin_width = 1.0_dp
REAL(KIND=dp), DIMENSION(3, 3) :: hmat = 0.0_dp
REAL(KIND=dp), DIMENSION(3, 3) :: h_inv = 0.0_dp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord_primary
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord_scaled
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ref_coord_primary
INTEGER, ALLOCATABLE, DIMENSION(:) :: image_atom
INTEGER, ALLOCATABLE, DIMENSION(:) :: head
INTEGER, ALLOCATABLE, DIMENSION(:) :: next
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: image_shift
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: image_translation
END TYPE nnp_cell_list_cache_type
! **************************************************************************************************
!> \brief Species-pair routing table. For one (central element, neighbour element)
!> pair, lists which radial / angular symmetry-function groups need that
!> neighbour and the worst-case relevant cutoff.
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
TYPE nnp_neighbor_pair_map_type
INTEGER :: n_rad = 0
INTEGER :: n_ang1 = 0
INTEGER :: n_ang2 = 0
REAL(KIND=dp) :: max_relevant_cutoff = 0.0_dp
INTEGER, ALLOCATABLE, DIMENSION(:) :: rad_groups
INTEGER, ALLOCATABLE, DIMENSION(:) :: ang1_groups
INTEGER, ALLOCATABLE, DIMENSION(:) :: ang2_groups
END TYPE nnp_neighbor_pair_map_type
! **************************************************************************************************
!> \brief Per-(element, SF-group) dG_k/dr buffer, sized to the group's n_symf (no
!> max_*_symf padding) and the observed peak neighbour count for the group.
!> Grown lazily (1.5x) so it settles at the true peak.
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
TYPE nnp_dGdr_grp_type
INTEGER :: cap = 0
INTEGER :: n_symf = 0
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: data ! (3, n_symf, cap)
END TYPE nnp_dGdr_grp_type
! **************************************************************************************************
!> \brief Reusable per-element scratch / persistent caches for the ACSF
!> descriptor and force assembly. See nnp_neighbor_interface for the
!> routines that populate and grow these buffers.
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
TYPE nnp_neighbor_workspace_type
INTEGER :: max_rad_symf = 0
INTEGER :: max_ang_symf = 0
INTEGER :: n_input_nodes = 0
! High-water mark of the per-element angular-cache slab, tracked so the 1D
! cutoff caches settle at the true per-element peak.
INTEGER :: cache_cap = 0
! There is a bug with some older compilers preventing requiring an explicit
! initialization of allocatable components (see mp2_types.F).
#if defined(FTN_NO_DEFAULT_INIT)
TYPE(nnp_neighbor_type) :: neighbor = nnp_neighbor_type( &
pbc_copies=-1, n_rad=NULL(), &
n_ang1=NULL(), n_ang2=NULL(), &
rad=NULL(), ang1=NULL(), ang2=NULL())
#else
TYPE(nnp_neighbor_type) :: neighbor = nnp_neighbor_type()
#endif
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radial_sym
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: radial_force
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: angular_sym
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: angular_force
! Per-element persistent cutoff caches reused across this element's central
! atoms. Sized lazily by nnp_workspace_grow_caches to MAX(n_ang1)/MAX(n_ang2).
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fc_cache1
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dfc_cache1
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fc_cache2
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dfc_cache2
! Per-neighbour dG_k/dr storage, replacing the global
! dsymdxyz(3, n_input_nodes, num_atoms) slab. The per-group buffers are
! dense in (3, n_symf_for_group, n_neighbors_in_group), which shrinks the
! working set and keeps each group walk a contiguous stride.
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: self_dGdr
TYPE(nnp_dGdr_grp_type), ALLOCATABLE, DIMENSION(:) :: dGdr_rad
TYPE(nnp_dGdr_grp_type), ALLOCATABLE, DIMENSION(:) :: dGdr_ang_jj
TYPE(nnp_dGdr_grp_type), ALLOCATABLE, DIMENSION(:) :: dGdr_ang_kk
END TYPE nnp_neighbor_workspace_type
! **************************************************************************************************
!> \brief Persistent neighbour-interface state. Owned by nnp_type so that the
!> per-nnp pair-routing tables and per-element workspaces survive
!> across MD steps without leaking between independent &NNP
!> force_evals. Operated on by nnp_neighbor_interface.F.
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
TYPE nnp_neighbor_interface_state_type
LOGICAL :: initialized = .FALSE.
INTEGER :: n_ele = 0
INTEGER, ALLOCATABLE, DIMENSION(:) :: n_rad
INTEGER, ALLOCATABLE, DIMENSION(:) :: n_ang
INTEGER, ALLOCATABLE, DIMENSION(:) :: n_radgrp
INTEGER, ALLOCATABLE, DIMENSION(:) :: n_anggrp
TYPE(nnp_neighbor_pair_map_type), ALLOCATABLE, &
DIMENSION(:, :) :: pair_map
TYPE(nnp_neighbor_workspace_type), ALLOCATABLE, &
DIMENSION(:) :: workspace
END TYPE nnp_neighbor_interface_state_type
! **************************************************************************************************
!> \brief Data type for artificial neural networks
!> \author Christoph Schran (christoph.schran@rub.de)
@ -280,9 +456,15 @@ CONTAINS
IF (ASSOCIATED(nnp_env%rad)) THEN
DO i = 1, nnp_env%n_ele
DO j = 1, nnp_env%rad(i)%n_symfgrp
DEALLOCATE (nnp_env%rad(i)%symfgrp(j)%symf, &
nnp_env%rad(i)%symfgrp(j)%ele, &
nnp_env%rad(i)%symfgrp(j)%ele_ind)
IF (ALLOCATED(nnp_env%rad(i)%symfgrp(j)%symf)) &
DEALLOCATE (nnp_env%rad(i)%symfgrp(j)%symf, &
nnp_env%rad(i)%symfgrp(j)%ele, &
nnp_env%rad(i)%symfgrp(j)%ele_ind)
IF (ALLOCATED(nnp_env%rad(i)%symfgrp(j)%spline_y)) &
DEALLOCATE (nnp_env%rad(i)%symfgrp(j)%spline_y)
IF (ALLOCATED(nnp_env%rad(i)%symfgrp(j)%spline_dy)) &
DEALLOCATE (nnp_env%rad(i)%symfgrp(j)%spline_dy)
nnp_env%rad(i)%symfgrp(j)%spline_built = .FALSE.
END DO
DEALLOCATE (nnp_env%rad(i)%y, &
nnp_env%rad(i)%funccut, &
@ -302,9 +484,22 @@ CONTAINS
IF (ASSOCIATED(nnp_env%ang)) THEN
DO i = 1, nnp_env%n_ele
DO j = 1, nnp_env%ang(i)%n_symfgrp
DEALLOCATE (nnp_env%ang(i)%symfgrp(j)%symf, &
nnp_env%ang(i)%symfgrp(j)%ele, &
nnp_env%ang(i)%symfgrp(j)%ele_ind)
IF (ALLOCATED(nnp_env%ang(i)%symfgrp(j)%symf)) &
DEALLOCATE (nnp_env%ang(i)%symfgrp(j)%symf, &
nnp_env%ang(i)%symfgrp(j)%ele, &
nnp_env%ang(i)%symfgrp(j)%ele_ind)
IF (ALLOCATED(nnp_env%ang(i)%symfgrp(j)%pack_eta)) &
DEALLOCATE (nnp_env%ang(i)%symfgrp(j)%pack_eta)
IF (ALLOCATED(nnp_env%ang(i)%symfgrp(j)%pack_zeta)) &
DEALLOCATE (nnp_env%ang(i)%symfgrp(j)%pack_zeta)
IF (ALLOCATED(nnp_env%ang(i)%symfgrp(j)%pack_lam)) &
DEALLOCATE (nnp_env%ang(i)%symfgrp(j)%pack_lam)
IF (ALLOCATED(nnp_env%ang(i)%symfgrp(j)%pack_prefzeta)) &
DEALLOCATE (nnp_env%ang(i)%symfgrp(j)%pack_prefzeta)
IF (ALLOCATED(nnp_env%ang(i)%symfgrp(j)%pack_izeta)) &
DEALLOCATE (nnp_env%ang(i)%symfgrp(j)%pack_izeta)
IF (ALLOCATED(nnp_env%ang(i)%symfgrp(j)%pack_use_int_zeta)) &
DEALLOCATE (nnp_env%ang(i)%symfgrp(j)%pack_use_int_zeta)
END DO
DEALLOCATE (nnp_env%ang(i)%y, &
nnp_env%ang(i)%funccut, &
@ -368,11 +563,16 @@ CONTAINS
IF (ALLOCATED(nnp_env%committee_forces)) DEALLOCATE (nnp_env%committee_forces)
IF (ALLOCATED(nnp_env%committee_stress)) DEALLOCATE (nnp_env%committee_stress)
IF (ALLOCATED(nnp_env%atoms)) DEALLOCATE (nnp_env%atoms)
IF (ALLOCATED(nnp_env%nnp_forces)) DEALLOCATE (nnp_env%nnp_forces)
IF (ASSOCIATED(nnp_env%subsys)) THEN
CALL cp_subsys_release(nnp_env%subsys)
IF (ALLOCATED(nnp_env%cell_list_cache)) THEN
CALL nnp_cell_list_cache_release(nnp_env%cell_list_cache)
DEALLOCATE (nnp_env%cell_list_cache)
END IF
IF (ALLOCATED(nnp_env%neighbor_interface_state)) THEN
CALL nnp_neighbor_interface_state_release(nnp_env%neighbor_interface_state)
DEALLOCATE (nnp_env%neighbor_interface_state)
END IF
IF (ASSOCIATED(nnp_env%subsys)) THEN
CALL cp_subsys_release(nnp_env%subsys)
END IF
@ -385,6 +585,139 @@ CONTAINS
END SUBROUTINE nnp_env_release
! **************************************************************************************************
!> \brief Free the allocatable parts of a cell-list cache. Co-located with the
!> type definition to avoid creating a USE-cycle through nnp_cell_list.
!> \param cache cell-list cache whose allocatable image-pool and bin arrays are freed
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
SUBROUTINE nnp_cell_list_cache_release(cache)
TYPE(nnp_cell_list_cache_type), INTENT(INOUT) :: cache
IF (ALLOCATED(cache%coord_primary)) DEALLOCATE (cache%coord_primary)
IF (ALLOCATED(cache%coord_scaled)) DEALLOCATE (cache%coord_scaled)
IF (ALLOCATED(cache%ref_coord_primary)) DEALLOCATE (cache%ref_coord_primary)
IF (ALLOCATED(cache%image_atom)) DEALLOCATE (cache%image_atom)
IF (ALLOCATED(cache%head)) DEALLOCATE (cache%head)
IF (ALLOCATED(cache%next)) DEALLOCATE (cache%next)
IF (ALLOCATED(cache%image_shift)) DEALLOCATE (cache%image_shift)
IF (ALLOCATED(cache%image_translation)) DEALLOCATE (cache%image_translation)
cache%initialized = .FALSE.
cache%num_atoms = -1
cache%n_images = 0
END SUBROUTINE nnp_cell_list_cache_release
! **************************************************************************************************
!> \brief Free the allocatable parts of a neighbour-interface state. Co-located
!> with the type definition to avoid creating a USE-cycle through
!> nnp_neighbor_interface.
!> \param state neighbour-interface state whose pair maps and per-element workspaces are freed
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
SUBROUTINE nnp_neighbor_interface_state_release(state)
TYPE(nnp_neighbor_interface_state_type), &
INTENT(INOUT) :: state
INTEGER :: i, j
IF (ALLOCATED(state%pair_map)) THEN
DO i = 1, SIZE(state%pair_map, 1)
DO j = 1, SIZE(state%pair_map, 2)
IF (ALLOCATED(state%pair_map(i, j)%rad_groups)) DEALLOCATE (state%pair_map(i, j)%rad_groups)
IF (ALLOCATED(state%pair_map(i, j)%ang1_groups)) DEALLOCATE (state%pair_map(i, j)%ang1_groups)
IF (ALLOCATED(state%pair_map(i, j)%ang2_groups)) DEALLOCATE (state%pair_map(i, j)%ang2_groups)
END DO
END DO
DEALLOCATE (state%pair_map)
END IF
IF (ALLOCATED(state%workspace)) THEN
DO i = 1, SIZE(state%workspace)
CALL nnp_workspace_release(state%workspace(i))
END DO
DEALLOCATE (state%workspace)
END IF
IF (ALLOCATED(state%n_rad)) DEALLOCATE (state%n_rad)
IF (ALLOCATED(state%n_ang)) DEALLOCATE (state%n_ang)
IF (ALLOCATED(state%n_radgrp)) DEALLOCATE (state%n_radgrp)
IF (ALLOCATED(state%n_anggrp)) DEALLOCATE (state%n_anggrp)
state%initialized = .FALSE.
state%n_ele = 0
END SUBROUTINE nnp_neighbor_interface_state_release
! **************************************************************************************************
!> \brief Free the allocatable parts of one per-element workspace.
!> \param workspace per-element workspace whose scratch, cutoff caches and dGdr buffers are freed
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
SUBROUTINE nnp_workspace_release(workspace)
TYPE(nnp_neighbor_workspace_type), INTENT(INOUT) :: workspace
INTEGER :: s
IF (ALLOCATED(workspace%neighbor%rad)) THEN
DO s = 1, SIZE(workspace%neighbor%rad)
IF (ALLOCATED(workspace%neighbor%rad(s)%ind)) DEALLOCATE (workspace%neighbor%rad(s)%ind)
IF (ALLOCATED(workspace%neighbor%rad(s)%dist)) DEALLOCATE (workspace%neighbor%rad(s)%dist)
END DO
DEALLOCATE (workspace%neighbor%rad)
END IF
IF (ALLOCATED(workspace%neighbor%ang1)) THEN
DO s = 1, SIZE(workspace%neighbor%ang1)
IF (ALLOCATED(workspace%neighbor%ang1(s)%ind)) DEALLOCATE (workspace%neighbor%ang1(s)%ind)
IF (ALLOCATED(workspace%neighbor%ang1(s)%dist)) DEALLOCATE (workspace%neighbor%ang1(s)%dist)
END DO
DEALLOCATE (workspace%neighbor%ang1)
END IF
IF (ALLOCATED(workspace%neighbor%ang2)) THEN
DO s = 1, SIZE(workspace%neighbor%ang2)
IF (ALLOCATED(workspace%neighbor%ang2(s)%ind)) DEALLOCATE (workspace%neighbor%ang2(s)%ind)
IF (ALLOCATED(workspace%neighbor%ang2(s)%dist)) DEALLOCATE (workspace%neighbor%ang2(s)%dist)
END DO
DEALLOCATE (workspace%neighbor%ang2)
END IF
IF (ALLOCATED(workspace%neighbor%n_rad)) DEALLOCATE (workspace%neighbor%n_rad)
IF (ALLOCATED(workspace%neighbor%n_ang1)) DEALLOCATE (workspace%neighbor%n_ang1)
IF (ALLOCATED(workspace%neighbor%n_ang2)) DEALLOCATE (workspace%neighbor%n_ang2)
workspace%neighbor%pbc_copies = -1
IF (ALLOCATED(workspace%radial_sym)) DEALLOCATE (workspace%radial_sym)
IF (ALLOCATED(workspace%radial_force)) DEALLOCATE (workspace%radial_force)
IF (ALLOCATED(workspace%angular_sym)) DEALLOCATE (workspace%angular_sym)
IF (ALLOCATED(workspace%angular_force)) DEALLOCATE (workspace%angular_force)
IF (ALLOCATED(workspace%fc_cache1)) DEALLOCATE (workspace%fc_cache1)
IF (ALLOCATED(workspace%dfc_cache1)) DEALLOCATE (workspace%dfc_cache1)
IF (ALLOCATED(workspace%fc_cache2)) DEALLOCATE (workspace%fc_cache2)
IF (ALLOCATED(workspace%dfc_cache2)) DEALLOCATE (workspace%dfc_cache2)
IF (ALLOCATED(workspace%self_dGdr)) DEALLOCATE (workspace%self_dGdr)
IF (ALLOCATED(workspace%dGdr_rad)) THEN
DO s = 1, SIZE(workspace%dGdr_rad)
IF (ALLOCATED(workspace%dGdr_rad(s)%data)) DEALLOCATE (workspace%dGdr_rad(s)%data)
END DO
DEALLOCATE (workspace%dGdr_rad)
END IF
IF (ALLOCATED(workspace%dGdr_ang_jj)) THEN
DO s = 1, SIZE(workspace%dGdr_ang_jj)
IF (ALLOCATED(workspace%dGdr_ang_jj(s)%data)) DEALLOCATE (workspace%dGdr_ang_jj(s)%data)
END DO
DEALLOCATE (workspace%dGdr_ang_jj)
END IF
IF (ALLOCATED(workspace%dGdr_ang_kk)) THEN
DO s = 1, SIZE(workspace%dGdr_ang_kk)
IF (ALLOCATED(workspace%dGdr_ang_kk(s)%data)) DEALLOCATE (workspace%dGdr_ang_kk(s)%data)
END DO
DEALLOCATE (workspace%dGdr_ang_kk)
END IF
workspace%max_rad_symf = 0
workspace%max_ang_symf = 0
workspace%n_input_nodes = 0
workspace%cache_cap = 0
END SUBROUTINE nnp_workspace_release
! **************************************************************************************************
!> \brief Returns various attributes of the nnp environment
!> \param nnp_env ...

219
src/nnp_force.F
View file

@ -8,6 +8,7 @@
! **************************************************************************************************
!> \brief Methods dealing with Neural Network potentials
!> \author Christoph Schran (christoph.schran@rub.de)
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
!> \date 2020-10-10
! **************************************************************************************************
MODULE nnp_force
@ -29,7 +30,8 @@ MODULE nnp_force
USE kinds, ONLY: default_path_length,&
default_string_length,&
dp
USE nnp_acsf, ONLY: nnp_calc_acsf
USE nnp_acsf, ONLY: nnp_calc_acsf,&
nnp_prepare_neighbor_cache
USE nnp_environment_types, ONLY: nnp_env_get,&
nnp_type
USE nnp_model, ONLY: nnp_gradients,&
@ -44,10 +46,10 @@ MODULE nnp_force
PRIVATE
LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_force'
PUBLIC :: nnp_calc_energy_force
PUBLIC :: nnp_calc_energy_force, nnp_scatter_dgdr_to_forces
CONTAINS
@ -64,12 +66,13 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'nnp_calc_energy_force'
INTEGER :: handle, i, i_com, ig, ind, istart, j, k, &
m, mecalc
INTEGER :: handle, handle_loop, i, i_com, ig, ind, &
istart, j, k, m, max_input_nodes, &
mecalc, n_input_nodes
INTEGER, ALLOCATABLE, DIMENSION(:) :: allcalc
LOGICAL :: calc_stress
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: denergydsym
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsymdxyz, stress
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: stress
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_subsys_type), POINTER :: subsys
@ -85,6 +88,7 @@ CONTAINS
logger => cp_get_default_logger()
CPASSERT(ASSOCIATED(nnp))
CPASSERT(nnp%n_committee >= 1)
CALL nnp_env_get(nnp_env=nnp, particle_set=particle_set, &
subsys=subsys, local_particles=local_particles, &
atomic_kind_set=atomic_kind_set)
@ -93,9 +97,9 @@ CONTAINS
virial=virial)
calc_stress = virial%pv_availability .AND. (.NOT. virial%pv_numer)
IF (calc_stress .AND. .NOT. calc_forces) CPABORT('Stress cannot be calculated without forces')
IF (calc_stress .AND. .NOT. calc_forces) &
CPABORT('Stress cannot be calculated without forces')
! Initialize energy and gradient:
nnp%atomic_energy(:, :) = 0.0_dp
IF (calc_forces) nnp%myforce(:, :, :) = 0.0_dp
IF (calc_forces) nnp%committee_forces(:, :, :) = 0.0_dp
@ -134,11 +138,32 @@ CONTAINS
! reset extrapolation status
nnp%output_expol = .FALSE.
CALL nnp_prepare_neighbor_cache(nnp)
max_input_nodes = 0
IF (calc_forces) THEN
DO i = 1, nnp%n_ele
max_input_nodes = MAX(max_input_nodes, nnp%arc(i)%n_nodes(1))
END DO
! dGdr lives in per-element workspace; no global (3, n_sf, num_atoms) slab.
IF (calc_stress) THEN
ALLOCATE (stress(3, 3, max_input_nodes))
stress(:, :, :) = 0.0_dp
END IF
ALLOCATE (denergydsym(max_input_nodes))
denergydsym(:) = 0.0_dp
END IF
! Per-rank timer around the per-atom work. CP2K aggregates the label as
! avg/max across ranks, so max - avg measures the load imbalance.
CALL timeset('nnp_per_rank_atom_loop', handle_loop)
! calc atomic contribution to energy and force
DO i = istart, istart + mecalc - 1
! determine index of atom type and offset
ind = nnp%ele_ind(i)
n_input_nodes = nnp%arc(ind)%n_nodes(1)
! reset input nodes of ele(ind):
nnp%arc(ind)%layer(1)%node(:) = 0.0_dp
@ -147,17 +172,14 @@ CONTAINS
IF (calc_forces) THEN
!reset input grads of ele(ind):
nnp%arc(ind)%layer(1)%node_grad(:) = 0.0_dp
ALLOCATE (dsymdxyz(3, nnp%arc(ind)%n_nodes(1), nnp%num_atoms))
dsymdxyz(:, :, :) = 0.0_dp
IF (calc_stress) THEN
ALLOCATE (stress(3, 3, nnp%arc(ind)%n_nodes(1)))
stress(:, :, :) = 0.0_dp
CALL nnp_calc_acsf(nnp, i, dsymdxyz, stress)
stress(:, :, 1:n_input_nodes) = 0.0_dp
CALL nnp_calc_acsf(nnp, i, .TRUE., stress(:, :, 1:n_input_nodes))
ELSE
CALL nnp_calc_acsf(nnp, i, dsymdxyz)
CALL nnp_calc_acsf(nnp, i, .TRUE.)
END IF
ELSE
CALL nnp_calc_acsf(nnp, i)
CALL nnp_calc_acsf(nnp, i, .FALSE.)
END IF
DO i_com = 1, nnp%n_committee
@ -168,34 +190,35 @@ CONTAINS
! predict forces
IF (calc_forces) THEN
ALLOCATE (denergydsym(nnp%arc(ind)%n_nodes(1)))
denergydsym(:) = 0.0_dp
CALL nnp_gradients(nnp%arc(ind), nnp, i_com, denergydsym)
DO j = 1, nnp%arc(ind)%n_nodes(1)
DO k = 1, nnp%num_atoms
DO m = 1, 3
nnp%myforce(m, k, i_com) = nnp%myforce(m, k, i_com) - denergydsym(j)*dsymdxyz(m, j, k)
END DO
END DO
IF (calc_stress) THEN
denergydsym(1:n_input_nodes) = 0.0_dp
CALL nnp_gradients(nnp%arc(ind), nnp, i_com, denergydsym(1:n_input_nodes))
! Force scatter over the per-element workspace dGdr arrays; see
! nnp_scatter_dgdr_to_forces below, also reused by the helium-NNP
! coupling in src/motion/helium_interactions.F.
CALL nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym(1:n_input_nodes), &
nnp%myforce(:, :, i_com))
IF (calc_stress) THEN
DO j = 1, n_input_nodes
nnp%committee_stress(:, :, i_com) = nnp%committee_stress(:, :, i_com) - &
denergydsym(j)*stress(:, :, j)
END IF
END DO
DEALLOCATE (denergydsym)
END DO
END IF
END IF
END DO
!deallocate memory
IF (calc_forces) THEN
DEALLOCATE (dsymdxyz)
IF (calc_stress) THEN
DEALLOCATE (stress)
END IF
END IF
END DO ! loop over num_atoms
CALL timestop(handle_loop)
IF (calc_forces) THEN
DEALLOCATE (denergydsym)
IF (calc_stress) THEN
DEALLOCATE (stress)
END IF
END IF
! calculate energy:
CALL logger%para_env%sum(nnp%atomic_energy(:, :))
nnp%committee_energy(:) = SUM(nnp%atomic_energy, 1)
@ -335,7 +358,7 @@ CONTAINS
print_key => section_vals_get_subs_vals(print_section, "FORCES")
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
IF (unit_nr > 0) CALL nnp_print_forces(nnp, print_key)
CALL nnp_print_forces(nnp, print_key)
END IF
print_key => section_vals_get_subs_vals(print_section, "FORCES_SIGMA")
@ -426,7 +449,6 @@ CONTAINS
NULLIFY (logger)
logger => cp_get_default_logger()
! TODO use write_particle_coordinates from particle_methods.F
DO i = 1, nnp%n_committee
WRITE (fmt_string, *) i
WRITE (middle_name, "(A,A)") "nnp-forces-", ADJUSTL(TRIM(fmt_string))
@ -612,7 +634,6 @@ CONTAINS
CHARACTER(len=default_path_length) :: fmt_string
INTEGER :: i
! TODO use write_particle_coordinates from particle_methods.F
WRITE (unit_nr, *) nnp%num_atoms
WRITE (unit_nr, "(A,F20.9)") "NNP bias forces [a.u.] for bias energy [a.u]=", nnp%bias_energy
@ -671,4 +692,124 @@ CONTAINS
END SUBROUTINE nnp_print_sumforces
! **************************************************************************************************
!> \brief Scatter the per-neighbour dG/dr arrays held in the nnp_neighbor_workspace
!> into a per-atom Cartesian force array. Three contributions per central
!> atom i:
!> 1) self: self_dGdr(:, m) -> force_xyz(:, i)
!> 2) radial: walk neighbor%rad(s)%ind(j), read dGdr_rad(s)%data(:, sf, j)
!> 3) angular: walk neighbor%ang1(s)%ind / ang2(s)%ind, read dGdr_ang_jj/kk
!> No global (n_sf, num_atoms) slab is touched. Reused by the helium-NNP
!> coupling in helium_interactions.
!>
!> Precondition: nnp_calc_acsf(nnp, i, calc_forces=.TRUE.[, stress]) must
!> have run for this same atom i immediately before; the per-element
!> workspace it fills is overwritten on every ACSF call.
!> \param nnp ...
!> \param ind central-atom element index
!> \param i central atom index (absolute)
!> \param denergydsym dE/dG_k for k = 1..n_input_nodes(ind)
!> \param force_xyz (3, num_atoms) destination -- accumulated, NOT overwritten
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
SUBROUTINE nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym, force_xyz)
TYPE(nnp_type), INTENT(IN) :: nnp
INTEGER, INTENT(IN) :: ind, i
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: denergydsym
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force_xyz
INTEGER :: j, k, k_atom, m, n_ang1_s, n_ang2_s, &
n_input_nodes, n_symf_s, off, s, sf
LOGICAL :: homo_grp
REAL(KIND=dp) :: de
ASSOCIATE (workspace => nnp%neighbor_interface_state%workspace(ind), &
neighbor => nnp%neighbor_interface_state%workspace(ind)%neighbor, &
self_dGdr => nnp%neighbor_interface_state%workspace(ind)%self_dGdr)
n_input_nodes = workspace%n_input_nodes
! Self contribution
DO j = 1, n_input_nodes
de = denergydsym(j)
force_xyz(1, i) = force_xyz(1, i) - de*self_dGdr(1, j)
force_xyz(2, i) = force_xyz(2, i) - de*self_dGdr(2, j)
force_xyz(3, i) = force_xyz(3, i) - de*self_dGdr(3, j)
END DO
! Radial neighbours
DO s = 1, nnp%rad(ind)%n_symfgrp
n_symf_s = nnp%rad(ind)%symfgrp(s)%n_symf
ASSOCIATE (rad_buf => workspace%dGdr_rad(s)%data)
DO j = 1, neighbor%n_rad(s)
k_atom = neighbor%rad(s)%ind(j)
DO sf = 1, n_symf_s
m = nnp%rad(ind)%symfgrp(s)%symf(sf)
de = denergydsym(m)
force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*rad_buf(1, sf, j)
force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*rad_buf(2, sf, j)
force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*rad_buf(3, sf, j)
END DO
END DO
END ASSOCIATE
END DO
! Angular neighbours
off = nnp%n_rad(ind)
DO s = 1, nnp%ang(ind)%n_symfgrp
n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
n_ang1_s = neighbor%n_ang1(s)
homo_grp = (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2))
IF (homo_grp) THEN
n_ang2_s = n_ang1_s ! kk buffer is also indexed in ang1(s)%ind
ELSE
n_ang2_s = neighbor%n_ang2(s)
END IF
ASSOCIATE (jj_buf => workspace%dGdr_ang_jj(s)%data, &
kk_buf => workspace%dGdr_ang_kk(s)%data)
! jj-side
DO j = 1, n_ang1_s
k_atom = neighbor%ang1(s)%ind(j)
DO sf = 1, n_symf_s
m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
de = denergydsym(m)
force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*jj_buf(1, sf, j)
force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*jj_buf(2, sf, j)
force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*jj_buf(3, sf, j)
END DO
END DO
! kk-side: homo reads ang1(s)%ind, hetero reads ang2(s)%ind
IF (homo_grp) THEN
DO k = 1, n_ang2_s
k_atom = neighbor%ang1(s)%ind(k)
DO sf = 1, n_symf_s
m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
de = denergydsym(m)
force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*kk_buf(1, sf, k)
force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*kk_buf(2, sf, k)
force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*kk_buf(3, sf, k)
END DO
END DO
ELSE
DO k = 1, n_ang2_s
k_atom = neighbor%ang2(s)%ind(k)
DO sf = 1, n_symf_s
m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
de = denergydsym(m)
force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*kk_buf(1, sf, k)
force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*kk_buf(2, sf, k)
force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*kk_buf(3, sf, k)
END DO
END DO
END IF
END ASSOCIATE
END DO
END ASSOCIATE
END SUBROUTINE nnp_scatter_dgdr_to_forces
END MODULE nnp_force

4
src/nnp_model.F
View file

@ -28,7 +28,7 @@ MODULE nnp_model
PRIVATE
LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_model'
PUBLIC :: nnp_write_arc, &
@ -169,7 +169,7 @@ CONTAINS
TYPE(nnp_arc_type), INTENT(INOUT) :: arc
TYPE(nnp_type), INTENT(INOUT) :: nnp
INTEGER, INTENT(IN) :: i_com
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: denergydsym
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: denergydsym
CHARACTER(len=*), PARAMETER :: routineN = 'nnp_gradients'

465
src/nnp_neighbor_interface.F Normal file
View file

@ -0,0 +1,465 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Per-nnp persistent neighbour-interface state for the NNP hot path.
!> Separates neighbour bookkeeping from the ACSF and network loops:
!> species-pair routing is precomputed once and per-element work
!> buffers are reused. State lives on the parent nnp_type, so each
!> &NNP force_eval keeps its own copy and the lifetime tracks
!> nnp_env_release.
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
!> \author Christoph Schran (christoph.schran@rub.de)
!> \date 2026-05-21
! **************************************************************************************************
MODULE nnp_neighbor_interface
USE kinds, ONLY: dp
USE nnp_environment_types, ONLY: nnp_dGdr_grp_type,&
nnp_neigh_grp_type,&
nnp_neighbor_interface_state_release,&
nnp_neighbor_workspace_type,&
nnp_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
PUBLIC :: nnp_neighbor_interface_prepare, &
nnp_neighbor_interface_reset_neighbor, &
nnp_grp_grow_dGdr, &
nnp_neigh_grp_grow, &
nnp_workspace_grow_caches
CONTAINS
! **************************************************************************************************
!> \brief Ensure pair-routing metadata and reusable workspaces are ready for the current NNP model.
!> \param nnp NNP environment whose neighbor_interface_state will be (re)built if needed.
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
SUBROUTINE nnp_neighbor_interface_prepare(nnp)
TYPE(nnp_type), INTENT(INOUT) :: nnp
CHARACTER(len=*), PARAMETER :: routineN = 'nnp_neighbor_interface_prepare'
INTEGER :: handle
LOGICAL :: rebuild
CALL timeset(routineN, handle)
rebuild = .NOT. nnp%neighbor_interface_state%initialized
IF (.NOT. rebuild) CALL nnp_neighbor_interface_needs_rebuild(nnp, rebuild)
IF (rebuild) THEN
CALL nnp_neighbor_interface_state_release(nnp%neighbor_interface_state)
CALL nnp_neighbor_interface_build_pair_maps(nnp)
END IF
CALL timestop(handle)
END SUBROUTINE nnp_neighbor_interface_prepare
! **************************************************************************************************
!> \brief Reset per-group neighbour counters for one central element before refilling
!> the reusable buffers. Zeroes n_rad/n_ang1/n_ang2; the (ind, dist) slabs
!> are kept allocated and overwritten in place by the next push.
!> \param nnp NNP environment whose neighbor_interface_state will have its counters cleared.
!> \param ind central-element index (1..nnp%n_ele) whose workspace is reset.
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
SUBROUTINE nnp_neighbor_interface_reset_neighbor(nnp, ind)
TYPE(nnp_type), INTENT(INOUT) :: nnp
INTEGER, INTENT(IN) :: ind
nnp%neighbor_interface_state%workspace(ind)%neighbor%n_rad(:) = 0
nnp%neighbor_interface_state%workspace(ind)%neighbor%n_ang1(:) = 0
nnp%neighbor_interface_state%workspace(ind)%neighbor%n_ang2(:) = 0
END SUBROUTINE nnp_neighbor_interface_reset_neighbor
! **************************************************************************************************
!> \brief Check whether the persistent state matches the current NNP model. Assumes
!> the per-group routing (ele_ind, cutoff) is fixed after nnp_init_acsf_groups,
!> so only element and per-element group counts are compared.
!> \param nnp NNP environment whose persistent neighbour-interface state will be inspected
!> \param rebuild (out) .TRUE. if element count or per-element group sizes have changed
! **************************************************************************************************
SUBROUTINE nnp_neighbor_interface_needs_rebuild(nnp, rebuild)
TYPE(nnp_type), INTENT(IN) :: nnp
LOGICAL, INTENT(OUT) :: rebuild
INTEGER :: i
rebuild = .FALSE.
IF (nnp%neighbor_interface_state%n_ele /= nnp%n_ele) THEN
rebuild = .TRUE.
RETURN
END IF
IF (.NOT. ALLOCATED(nnp%neighbor_interface_state%n_rad)) THEN
rebuild = .TRUE.
RETURN
END IF
DO i = 1, nnp%n_ele
IF (nnp%neighbor_interface_state%n_rad(i) /= nnp%n_rad(i) .OR. &
nnp%neighbor_interface_state%n_ang(i) /= nnp%n_ang(i) .OR. &
nnp%neighbor_interface_state%n_radgrp(i) /= nnp%rad(i)%n_symfgrp .OR. &
nnp%neighbor_interface_state%n_anggrp(i) /= nnp%ang(i)%n_symfgrp) THEN
rebuild = .TRUE.
RETURN
END IF
END DO
END SUBROUTINE nnp_neighbor_interface_needs_rebuild
! **************************************************************************************************
!> \brief Build species-pair routing tables and initialize reusable workspaces.
!> \param nnp NNP environment; pair_map and workspace arrays are (re-)allocated on its neighbor_interface_state
! **************************************************************************************************
SUBROUTINE nnp_neighbor_interface_build_pair_maps(nnp)
TYPE(nnp_type), INTENT(INOUT) :: nnp
INTEGER :: i
ASSOCIATE (state => nnp%neighbor_interface_state)
state%n_ele = nnp%n_ele
ALLOCATE (state%n_rad(nnp%n_ele))
ALLOCATE (state%n_ang(nnp%n_ele))
ALLOCATE (state%n_radgrp(nnp%n_ele))
ALLOCATE (state%n_anggrp(nnp%n_ele))
ALLOCATE (state%pair_map(nnp%n_ele, nnp%n_ele))
ALLOCATE (state%workspace(nnp%n_ele))
DO i = 1, nnp%n_ele
state%n_rad(i) = nnp%n_rad(i)
state%n_ang(i) = nnp%n_ang(i)
state%n_radgrp(i) = nnp%rad(i)%n_symfgrp
state%n_anggrp(i) = nnp%ang(i)%n_symfgrp
END DO
DO i = 1, nnp%n_ele
CALL nnp_neighbor_interface_build_pair_map_for_element(nnp, i)
CALL nnp_neighbor_interface_init_workspace_metadata(nnp, i)
END DO
state%initialized = .TRUE.
END ASSOCIATE
END SUBROUTINE nnp_neighbor_interface_build_pair_maps
! **************************************************************************************************
!> \brief Build all species-pair routes for one central element.
!> \param nnp NNP environment providing the radial/angular SF groups
!> \param ind central-element index whose pair-map row is built
! **************************************************************************************************
SUBROUTINE nnp_neighbor_interface_build_pair_map_for_element(nnp, ind)
TYPE(nnp_type), INTENT(INOUT) :: nnp
INTEGER, INTENT(IN) :: ind
INTEGER :: idx, neighbor_ind, s
DO neighbor_ind = 1, nnp%n_ele
ASSOCIATE (pair_map => nnp%neighbor_interface_state%pair_map(ind, neighbor_ind))
pair_map%n_rad = 0
pair_map%n_ang1 = 0
pair_map%n_ang2 = 0
pair_map%max_relevant_cutoff = 0.0_dp
DO s = 1, nnp%rad(ind)%n_symfgrp
IF (nnp%rad(ind)%symfgrp(s)%ele_ind(1) == neighbor_ind) pair_map%n_rad = pair_map%n_rad + 1
END DO
DO s = 1, nnp%ang(ind)%n_symfgrp
IF (nnp%ang(ind)%symfgrp(s)%ele_ind(1) == neighbor_ind) pair_map%n_ang1 = pair_map%n_ang1 + 1
IF (nnp%ang(ind)%symfgrp(s)%ele_ind(2) == neighbor_ind) pair_map%n_ang2 = pair_map%n_ang2 + 1
END DO
ALLOCATE (pair_map%rad_groups(MAX(1, pair_map%n_rad)))
ALLOCATE (pair_map%ang1_groups(MAX(1, pair_map%n_ang1)))
ALLOCATE (pair_map%ang2_groups(MAX(1, pair_map%n_ang2)))
idx = 0
DO s = 1, nnp%rad(ind)%n_symfgrp
IF (nnp%rad(ind)%symfgrp(s)%ele_ind(1) == neighbor_ind) THEN
idx = idx + 1
pair_map%rad_groups(idx) = s
pair_map%max_relevant_cutoff = MAX(pair_map%max_relevant_cutoff, nnp%rad(ind)%symfgrp(s)%cutoff)
END IF
END DO
idx = 0
DO s = 1, nnp%ang(ind)%n_symfgrp
IF (nnp%ang(ind)%symfgrp(s)%ele_ind(1) == neighbor_ind) THEN
idx = idx + 1
pair_map%ang1_groups(idx) = s
pair_map%max_relevant_cutoff = MAX(pair_map%max_relevant_cutoff, nnp%ang(ind)%symfgrp(s)%cutoff)
END IF
END DO
idx = 0
DO s = 1, nnp%ang(ind)%n_symfgrp
IF (nnp%ang(ind)%symfgrp(s)%ele_ind(2) == neighbor_ind) THEN
idx = idx + 1
pair_map%ang2_groups(idx) = s
pair_map%max_relevant_cutoff = MAX(pair_map%max_relevant_cutoff, nnp%ang(ind)%symfgrp(s)%cutoff)
END IF
END DO
END ASSOCIATE
END DO
END SUBROUTINE nnp_neighbor_interface_build_pair_map_for_element
! **************************************************************************************************
!> \brief Cache max scratch sizes for one central element.
!> \param nnp NNP environment providing per-element SF group definitions
!> \param ind central-element index whose scratch workspace metadata is cached
! **************************************************************************************************
SUBROUTINE nnp_neighbor_interface_init_workspace_metadata(nnp, ind)
TYPE(nnp_type), INTENT(INOUT) :: nnp
INTEGER, INTENT(IN) :: ind
INTEGER :: s
ASSOCIATE (workspace => nnp%neighbor_interface_state%workspace(ind))
workspace%max_rad_symf = 0
workspace%max_ang_symf = 0
workspace%n_input_nodes = nnp%n_rad(ind) + nnp%n_ang(ind)
DO s = 1, nnp%rad(ind)%n_symfgrp
workspace%max_rad_symf = MAX(workspace%max_rad_symf, nnp%rad(ind)%symfgrp(s)%n_symf)
END DO
DO s = 1, nnp%ang(ind)%n_symfgrp
workspace%max_ang_symf = MAX(workspace%max_ang_symf, nnp%ang(ind)%symfgrp(s)%n_symf)
END DO
! Per-SF scratch reused inside nnp_calc_rad / nnp_calc_ang, sized by max_*_symf.
IF (ALLOCATED(workspace%radial_sym)) DEALLOCATE (workspace%radial_sym)
IF (ALLOCATED(workspace%radial_force)) DEALLOCATE (workspace%radial_force)
IF (ALLOCATED(workspace%angular_sym)) DEALLOCATE (workspace%angular_sym)
IF (ALLOCATED(workspace%angular_force)) DEALLOCATE (workspace%angular_force)
ALLOCATE (workspace%radial_sym(MAX(1, workspace%max_rad_symf)))
ALLOCATE (workspace%radial_force(3, MAX(1, workspace%max_rad_symf)))
ALLOCATE (workspace%angular_sym(MAX(1, workspace%max_ang_symf)))
ALLOCATE (workspace%angular_force(3, 3, MAX(1, workspace%max_ang_symf)))
! 1D angular cutoff caches; start small and grow lazily to the peak
! per-element angular neighbour count (nnp_workspace_grow_caches).
IF (ALLOCATED(workspace%fc_cache1)) DEALLOCATE (workspace%fc_cache1)
IF (ALLOCATED(workspace%dfc_cache1)) DEALLOCATE (workspace%dfc_cache1)
IF (ALLOCATED(workspace%fc_cache2)) DEALLOCATE (workspace%fc_cache2)
IF (ALLOCATED(workspace%dfc_cache2)) DEALLOCATE (workspace%dfc_cache2)
workspace%cache_cap = 8
ALLOCATE (workspace%fc_cache1(workspace%cache_cap))
ALLOCATE (workspace%dfc_cache1(workspace%cache_cap))
ALLOCATE (workspace%fc_cache2(workspace%cache_cap))
ALLOCATE (workspace%dfc_cache2(workspace%cache_cap))
! self_dGdr sized once here; independent of the candidate pool.
IF (ALLOCATED(workspace%self_dGdr)) DEALLOCATE (workspace%self_dGdr)
ALLOCATE (workspace%self_dGdr(3, MAX(1, workspace%n_input_nodes)))
! Per-group neighbour counters and dense (ind, dist) containers, sized to
! n_symfgrp here; the per-group slabs grow lazily via nnp_neigh_grp_grow.
CALL nnp_release_neighbor_local(workspace)
ALLOCATE (workspace%neighbor%n_rad(MAX(1, nnp%rad(ind)%n_symfgrp)))
ALLOCATE (workspace%neighbor%n_ang1(MAX(1, nnp%ang(ind)%n_symfgrp)))
ALLOCATE (workspace%neighbor%n_ang2(MAX(1, nnp%ang(ind)%n_symfgrp)))
ALLOCATE (workspace%neighbor%rad(MAX(1, nnp%rad(ind)%n_symfgrp)))
ALLOCATE (workspace%neighbor%ang1(MAX(1, nnp%ang(ind)%n_symfgrp)))
ALLOCATE (workspace%neighbor%ang2(MAX(1, nnp%ang(ind)%n_symfgrp)))
workspace%neighbor%n_rad(:) = 0
workspace%neighbor%n_ang1(:) = 0
workspace%neighbor%n_ang2(:) = 0
workspace%neighbor%pbc_copies = 0
! Per-group dG/dr buffers: container arrays sized to n_symfgrp and each
! entry's n_symf recorded. The %data slab stays unallocated until
! nnp_grp_grow_dGdr sees a real neighbour count.
CALL nnp_release_dGdr_grp_array(workspace%dGdr_rad)
CALL nnp_release_dGdr_grp_array(workspace%dGdr_ang_jj)
CALL nnp_release_dGdr_grp_array(workspace%dGdr_ang_kk)
ALLOCATE (workspace%dGdr_rad(MAX(1, nnp%rad(ind)%n_symfgrp)))
ALLOCATE (workspace%dGdr_ang_jj(MAX(1, nnp%ang(ind)%n_symfgrp)))
ALLOCATE (workspace%dGdr_ang_kk(MAX(1, nnp%ang(ind)%n_symfgrp)))
DO s = 1, nnp%rad(ind)%n_symfgrp
workspace%dGdr_rad(s)%n_symf = nnp%rad(ind)%symfgrp(s)%n_symf
workspace%dGdr_rad(s)%cap = 0
END DO
DO s = 1, nnp%ang(ind)%n_symfgrp
workspace%dGdr_ang_jj(s)%n_symf = nnp%ang(ind)%symfgrp(s)%n_symf
workspace%dGdr_ang_jj(s)%cap = 0
workspace%dGdr_ang_kk(s)%n_symf = nnp%ang(ind)%symfgrp(s)%n_symf
workspace%dGdr_ang_kk(s)%cap = 0
END DO
END ASSOCIATE
END SUBROUTINE nnp_neighbor_interface_init_workspace_metadata
! **************************************************************************************************
!> \brief Release the workspace's neighbor%(rad,ang1,ang2) per-group slabs before
!> re-allocation, keeping the persistent caches that survive a rebuild.
!> \param workspace per-element scratch workspace whose neighbor%(rad,ang1,ang2) slabs will be released
! **************************************************************************************************
SUBROUTINE nnp_release_neighbor_local(workspace)
TYPE(nnp_neighbor_workspace_type), INTENT(INOUT) :: workspace
INTEGER :: s
IF (ALLOCATED(workspace%neighbor%rad)) THEN
DO s = 1, SIZE(workspace%neighbor%rad)
IF (ALLOCATED(workspace%neighbor%rad(s)%ind)) DEALLOCATE (workspace%neighbor%rad(s)%ind)
IF (ALLOCATED(workspace%neighbor%rad(s)%dist)) DEALLOCATE (workspace%neighbor%rad(s)%dist)
END DO
DEALLOCATE (workspace%neighbor%rad)
END IF
IF (ALLOCATED(workspace%neighbor%ang1)) THEN
DO s = 1, SIZE(workspace%neighbor%ang1)
IF (ALLOCATED(workspace%neighbor%ang1(s)%ind)) DEALLOCATE (workspace%neighbor%ang1(s)%ind)
IF (ALLOCATED(workspace%neighbor%ang1(s)%dist)) DEALLOCATE (workspace%neighbor%ang1(s)%dist)
END DO
DEALLOCATE (workspace%neighbor%ang1)
END IF
IF (ALLOCATED(workspace%neighbor%ang2)) THEN
DO s = 1, SIZE(workspace%neighbor%ang2)
IF (ALLOCATED(workspace%neighbor%ang2(s)%ind)) DEALLOCATE (workspace%neighbor%ang2(s)%ind)
IF (ALLOCATED(workspace%neighbor%ang2(s)%dist)) DEALLOCATE (workspace%neighbor%ang2(s)%dist)
END DO
DEALLOCATE (workspace%neighbor%ang2)
END IF
IF (ALLOCATED(workspace%neighbor%n_rad)) DEALLOCATE (workspace%neighbor%n_rad)
IF (ALLOCATED(workspace%neighbor%n_ang1)) DEALLOCATE (workspace%neighbor%n_ang1)
IF (ALLOCATED(workspace%neighbor%n_ang2)) DEALLOCATE (workspace%neighbor%n_ang2)
workspace%neighbor%pbc_copies = -1
END SUBROUTINE nnp_release_neighbor_local
! **************************************************************************************************
!> \brief Release every per-group dG/dr buffer in a container array, then deallocate the container.
!> \param grps per-group dG/dr container array to release (each entry's %data is freed first)
! **************************************************************************************************
SUBROUTINE nnp_release_dGdr_grp_array(grps)
TYPE(nnp_dGdr_grp_type), ALLOCATABLE, &
INTENT(INOUT) :: grps(:)
INTEGER :: s
IF (.NOT. ALLOCATED(grps)) RETURN
DO s = 1, SIZE(grps)
IF (ALLOCATED(grps(s)%data)) DEALLOCATE (grps(s)%data)
grps(s)%cap = 0
grps(s)%n_symf = 0
END DO
DEALLOCATE (grps)
END SUBROUTINE nnp_release_dGdr_grp_array
! **************************************************************************************************
!> \brief Ensure a per-group dG/dr buffer holds n_needed neighbours, growing by
!> 1.5x (with an additive floor) so reallocation amortizes to O(1) over a
!> trajectory. grp%data is allocated on return (cap >= 1), so callers can
!> ASSOCIATE-bind it even for a zero-trip loop.
!> \param grp per-element dGdr group whose %data slab is reallocated if undersized
!> \param n_needed minimum required capacity (third dim of grp%data)
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
SUBROUTINE nnp_grp_grow_dGdr(grp, n_needed)
TYPE(nnp_dGdr_grp_type), INTENT(INOUT) :: grp
INTEGER, INTENT(IN) :: n_needed
INTEGER :: new_cap
IF (ALLOCATED(grp%data) .AND. grp%cap >= n_needed) RETURN
IF (ALLOCATED(grp%data)) DEALLOCATE (grp%data)
new_cap = MAX(MAX(1, n_needed), INT(grp%cap*1.5_dp) + 8)
ALLOCATE (grp%data(3, MAX(1, grp%n_symf), new_cap))
grp%cap = new_cap
END SUBROUTINE nnp_grp_grow_dGdr
! **************************************************************************************************
!> \brief Ensure a per-group (ind, dist) neighbour buffer holds n_needed entries.
!> Called from inside the linked-cell push loop, so existing entries are
!> preserved on grow via MOVE_ALLOC. 1.5x growth amortizes reallocation to
!> O(1) over a trajectory.
!> \param grp per-species-pair neighbour group whose (ind, dist) buffers grow
!> \param n_needed minimum required entry count
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
SUBROUTINE nnp_neigh_grp_grow(grp, n_needed)
TYPE(nnp_neigh_grp_type), INTENT(INOUT) :: grp
INTEGER, INTENT(IN) :: n_needed
INTEGER :: n_old, new_cap
INTEGER, ALLOCATABLE :: new_ind(:)
REAL(KIND=dp), ALLOCATABLE :: new_dist(:, :)
IF (ALLOCATED(grp%dist) .AND. grp%cap >= n_needed) RETURN
new_cap = MAX(MAX(8, n_needed), INT(grp%cap*1.5_dp) + 8)
IF (ALLOCATED(grp%dist)) THEN
n_old = grp%cap
ALLOCATE (new_dist(4, new_cap))
ALLOCATE (new_ind(new_cap))
IF (n_old > 0) THEN
new_dist(:, 1:n_old) = grp%dist(:, 1:n_old)
new_ind(1:n_old) = grp%ind(1:n_old)
END IF
CALL MOVE_ALLOC(new_dist, grp%dist)
CALL MOVE_ALLOC(new_ind, grp%ind)
ELSE
ALLOCATE (grp%dist(4, new_cap))
ALLOCATE (grp%ind(new_cap))
END IF
grp%cap = new_cap
END SUBROUTINE nnp_neigh_grp_grow
! **************************************************************************************************
!> \brief Ensure the four per-element angular cutoff caches hold at least n_needed
!> entries. These 1D scratch arrays are reused across angular groups within
!> one nnp_calc_acsf call, so the peak is MAX_s(n_ang1) and MAX_s(n_ang2).
!> 1.5x lazy growth.
!> \param workspace per-element workspace whose fc_cache1/dfc_cache1/fc_cache2/dfc_cache2 grow
!> \param n_needed minimum required cache length
!> \author Dhruv Sharma (ds2173@cam.ac.uk)
! **************************************************************************************************
SUBROUTINE nnp_workspace_grow_caches(workspace, n_needed)
TYPE(nnp_neighbor_workspace_type), INTENT(INOUT) :: workspace
INTEGER, INTENT(IN) :: n_needed
INTEGER :: new_cap
IF (workspace%cache_cap >= n_needed) RETURN
new_cap = MAX(MAX(8, n_needed), INT(workspace%cache_cap*1.5_dp) + 8)
! Contents not preserved: nnp_fill_fc_dfc_cache overwrites fully before any read.
IF (ALLOCATED(workspace%fc_cache1)) DEALLOCATE (workspace%fc_cache1)
IF (ALLOCATED(workspace%dfc_cache1)) DEALLOCATE (workspace%dfc_cache1)
IF (ALLOCATED(workspace%fc_cache2)) DEALLOCATE (workspace%fc_cache2)
IF (ALLOCATED(workspace%dfc_cache2)) DEALLOCATE (workspace%dfc_cache2)
ALLOCATE (workspace%fc_cache1(new_cap))
ALLOCATE (workspace%dfc_cache1(new_cap))
ALLOCATE (workspace%fc_cache2(new_cap))
ALLOCATE (workspace%dfc_cache2(new_cap))
workspace%cache_cap = new_cap
END SUBROUTINE nnp_workspace_grow_caches
END MODULE nnp_neighbor_interface

37
src/start/input_cp2k_motion.F
View file

@ -2334,7 +2334,7 @@ CONTAINS
CALL section_create(subsection, __LOCATION__, name="NNP", &
description="This section contains all information to run an helium-solute "// &
"interaction Neural Network Potential (NNP) calculation.", &
n_keywords=2, n_subsections=3, repeats=.FALSE.)
n_keywords=4, n_subsections=3, repeats=.FALSE.)
CALL keyword_create(keyword, __LOCATION__, name="NNP_INPUT_FILE_NAME", &
description="File containing the input information for the setup "// &
@ -2350,6 +2350,41 @@ CONTAINS
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="RAD_SPLINE_N", &
description="Number of knots per radial group in the cubic-Hermite "// &
"spline tables that tabulate the radial symmetry-function product "// &
"y(r) = exp(-eta*(r-rs)^2) * fcut(r). Cubic-Hermite per-evaluation "// &
"error scales as O(1/n^4); the default keeps the spline residual "// &
"inside the NNP regression tolerance for a radial cutoff of "// &
"~12 bohr. Models with larger cutoffs or stricter tolerances may "// &
"need a larger n. Memory cost scales linearly in n per radial group.", &
repeats=.FALSE., &
n_var=1, &
type_of_var=integer_t, &
default_i_val=8192, &
usage="RAD_SPLINE_N 8192")
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="VERLET_SKIN", &
description="Verlet skin distance for the NNP descriptor neighbour "// &
"cell-list. The neighbour-list cutoff is (symmetry-function cutoff "// &
"+ skin); the cell-list chain is rebuilt only when an atom drifts "// &
"more than skin/2 from its rebuild-time position, analogous to the "// &
"LAMMPS 'neighbor <skin> bin' command. Larger skin reduces the "// &
"rebuild rate but enlarges the per-atom neighbour list. A negative "// &
"value (default) selects the automatic heuristic "// &
"MIN(0.5 bohr, 0.1*cutoff). Useful upper bound is half the smallest "// &
"perpendicular cell width.", &
repeats=.FALSE., &
n_var=1, &
type_of_var=real_t, &
default_r_val=-1.0_dp, &
unit_str="bohr", &
usage="VERLET_SKIN [bohr] 0.5")
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
NULLIFY (subsubsection)
CALL section_create(subsubsection, __LOCATION__, name="SR_CUTOFF", &
description="Section for failsafe short range cutoffs for the NNPs, "// &

4
tests/HUGE_TESTS_SUPPRESSIONS
View file

@ -167,4 +167,8 @@ QS/regtest-smeagol-2/au-h2-au-1x1x1-1V-forces.inp
# Tests produced 11.69 MiB of output.
QS/regtest-openpmd/H2O-geo-ot-lumo-all-openpmd.inp
# Test produced 9.87 MiB of output (per-step &FORCES print x 8 committee
# models x finite-difference stress perturbations).
NNP/regtest-1/H2O-64_C-NNP_MD-NpT-numeric.inp
#EOF

851
tests/NNP/regtest-1/H2O-256_C-NNP_MD-bitexact.inp Normal file
View file

@ -0,0 +1,851 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT H2O-256_C-NNP_MD-bitexact
RUN_TYPE MD
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 5
TEMPERATURE 300
TIMESTEP 0.5
&PRINT
FORCE_LAST
&ENERGY
&EACH
MD 50
&END EACH
&END ENERGY
&END PRINT
&THERMOSTAT
REGION GLOBAL
TYPE CSVR
&CSVR
TIMECON 30.00
&END CSVR
&END THERMOSTAT
&END MD
&PRINT
&TRAJECTORY
&EACH
MD 500
&END EACH
&END TRAJECTORY
&END PRINT
&END MOTION
&FORCE_EVAL
METHOD NNP
&NNP
NNP_INPUT_FILE_NAME NNP/bulkH2O-jcp2020-cnnp/nnp-1/input.nn
SCALE_FILE_NAME NNP/bulkH2O-jcp2020-cnnp/nnp-1/scaling.data
&MODEL
WEIGHTS NNP/bulkH2O-jcp2020-cnnp/nnp-1/weights
&END MODEL
&MODEL
WEIGHTS NNP/bulkH2O-jcp2020-cnnp/nnp-2/weights
&END MODEL
&MODEL
WEIGHTS NNP/bulkH2O-jcp2020-cnnp/nnp-3/weights
&END MODEL
&MODEL
WEIGHTS NNP/bulkH2O-jcp2020-cnnp/nnp-4/weights
&END MODEL
&MODEL
WEIGHTS NNP/bulkH2O-jcp2020-cnnp/nnp-5/weights
&END MODEL
&MODEL
WEIGHTS NNP/bulkH2O-jcp2020-cnnp/nnp-6/weights
&END MODEL
&MODEL
WEIGHTS NNP/bulkH2O-jcp2020-cnnp/nnp-7/weights
&END MODEL
&MODEL
WEIGHTS NNP/bulkH2O-jcp2020-cnnp/nnp-8/weights
&END MODEL
&PRINT
&ENERGIES
&EACH
MD 100
&END EACH
&END ENERGIES
&END PRINT
&END NNP
&SUBSYS
&CELL
ABC [angstrom] 24.84000000 24.84000000 12.42000000
PERIODIC XYZ
&END CELL
&COORD
O 17.85470000 14.89360000 4.30373000
H 18.81710000 14.81910000 4.30668000
H 17.60460000 14.06810000 4.71635000
O 5.08473000 -2.09187000 -12.52020000
H 5.79642000 -1.44510000 -12.70560000
H 5.42957000 -2.91898000 -12.19440000
O 4.02996000 7.60544000 -3.68524000
H 4.27003000 6.94768000 -4.41446000
H 4.77458000 8.28889000 -3.67208000
O -12.22650000 18.89220000 12.25860000
H -13.08640000 19.19600000 12.55890000
H -12.45140000 18.59980000 11.38190000
O 1.36461000 4.68750000 16.07470000
H 0.77220000 5.33188000 15.70890000
H 0.92914100 4.08655000 16.73020000
O -3.42189000 3.99641000 12.57220000
H -2.47934000 4.05205000 12.38670000
H -3.57781000 4.18579000 13.49940000
O -12.45810000 5.26919000 -2.78020000
H -11.67410000 4.62996000 -2.73820000
H -13.15450000 4.83976000 -3.28648000
O -12.87460000 14.22990000 10.97500000
H -12.40850000 13.48810000 11.29410000
H -12.17020000 14.92310000 10.79950000
O 16.21660000 10.92160000 2.09580000
H 15.77660000 10.15310000 2.51719000
H 16.70080000 10.54170000 1.31827000
O 10.17360000 12.92940000 3.09205000
H 10.77360000 12.63090000 3.73607000
H 9.56179000 12.13220000 2.94022000
O 7.77605000 4.98249000 -2.29855000
H 8.07228000 4.60108000 -1.47153000
H 6.99714000 5.40667000 -2.04022000
O 14.58280000 2.43851000 2.33359000
H 13.85980000 2.38788000 1.70989000
H 14.40950000 3.36726000 2.62982000
O -6.66475000 9.71604000 21.73940000
H -6.28708000 10.17630000 20.91490000
H -5.93185000 9.37252000 22.22920000
O 1.00730000 7.62893000 -3.71784000
H 1.95122000 7.70838000 -3.49284000
H 0.78209600 6.70589000 -3.39957000
O -2.00359000 8.81074000 4.64761000
H -2.33580000 9.40933000 3.96349000
H -1.02960000 8.88476000 4.52207000
O 12.25120000 10.23950000 9.10079000
H 12.66980000 10.95650000 8.67039000
H 12.71420000 9.38279000 8.86118000
O -0.24441300 15.84350000 13.25460000
H 0.00470418 16.62720000 12.75340000
H -0.42289200 15.22620000 12.52560000
O 17.25350000 -6.21439000 18.95750000
H 17.26240000 -7.18011000 18.92590000
H 16.44650000 -5.96877000 18.52680000
O 1.57232000 8.55013000 3.48535000
H 1.78649000 8.35399000 4.39869000
H 0.91633100 7.79040000 3.31060000
O 9.73687000 10.53440000 -12.33500000
H 9.33966000 11.43670000 -12.37350000
H 10.71770000 10.62350000 -12.51260000
O 3.19186000 -1.70399000 6.69143000
H 2.64086000 -0.91080900 6.70096000
H 3.29213000 -1.99342000 7.64120000
O 5.01286000 3.41556000 7.01155000
H 5.25540000 3.30163000 7.95343000
H 5.82406000 3.47012000 6.52344000
O 16.21560000 -0.86663700 -2.85105000
H 17.05640000 -1.46010000 -2.85949000
H 15.82330000 -0.94576600 -1.95729000
O 18.39240000 15.45230000 9.37872000
H 17.70210000 15.50250000 10.02480000
H 19.03780000 16.17220000 9.59111000
O 12.43740000 -8.71873000 -6.55194000
H 13.24560000 -9.18214000 -6.25474000
H 11.74430000 -9.43587000 -6.45253000
O -2.12768000 -4.80632000 13.63960000
H -1.96159000 -3.92263000 13.27480000
H -3.14023000 -4.92484000 13.58540000
O 5.47695000 -12.78860000 18.12950000
H 6.04948000 -12.92690000 18.88220000
H 4.74944000 -13.43130000 18.27050000
O 9.06380000 13.52840000 12.92680000
H 9.07296000 14.42340000 12.64730000
H 9.45712000 13.49050000 13.82760000
O 6.54756000 -11.95970000 12.30680000
H 5.90365000 -11.48740000 12.84550000
H 7.46789000 -11.66100000 12.52900000
O -3.25432000 -2.19516000 15.18330000
H -3.21119000 -2.38945000 14.23370000
H -4.18029000 -2.00474000 15.49650000
O 19.16120000 20.64270000 4.17300000
H 18.25990000 20.42210000 4.00929000
H 19.17920000 21.59520000 4.03867000
O 16.62900000 4.75612000 3.61679000
H 17.41650000 4.29018000 3.87743000
H 15.90400000 4.43119000 4.10392000
O -14.46860000 1.26932000 -6.11559000
H -13.90410000 0.58683400 -6.47353000
H -14.78000000 1.12280000 -5.19548000
O 2.42216000 12.91590000 4.10819000
H 1.94154000 13.51420000 3.57260000
H 2.79204000 12.34870000 3.45220000
O 17.30610000 14.63850000 1.43734000
H 16.39600000 14.75480000 1.71543000
H 17.68200000 14.67460000 2.27931000
O 3.42317000 13.87760000 8.37823000
H 4.23289000 13.97000000 7.90067000
H 3.52615000 13.14510000 8.97284000
O -4.91217000 23.61640000 20.12810000
H -4.02584000 23.47060000 19.75830000
H -4.83605000 24.14340000 20.97620000
O 20.30930000 8.75645000 -1.55507000
H 21.03230000 9.38853000 -1.42815000
H 20.44040000 8.12820000 -2.30558000
O 17.63270000 6.21611000 -1.45303000
H 17.18510000 6.65247000 -2.18930000
H 17.59600000 6.94078000 -0.77261100
O -11.85650000 -3.09341000 6.38427000
H -11.91870000 -3.99992000 6.73656000
H -10.96070000 -2.82121000 6.56532000
O 14.34220000 -5.83286000 5.75163000
H 14.20720000 -6.54550000 5.13562000
H 14.05200000 -6.07263000 6.61672000
O -1.09658000 6.20827000 5.58702000
H -2.04987000 5.94104000 5.55373000
H -0.68250800 5.41005000 5.89003000
O 8.13776000 2.00921000 16.81480000
H 8.76908000 1.47311000 16.26590000
H 8.54559000 2.38894000 17.63890000
O 0.05855300 -0.59401100 17.45460000
H 0.96214000 -0.31451700 17.27990000
H 0.08964280 -1.51268000 17.79400000
O 27.67350000 0.60729400 12.06030000
H 28.34720000 0.76104100 12.69580000
H 27.19490000 -0.15785700 12.45860000
O 9.89541000 -7.95756000 8.17742000
H 9.62897000 -8.79966000 8.58524000
H 9.28724000 -7.36014000 8.69880000
O -4.41370000 20.18930000 8.05879000
H -3.76101000 20.81950000 7.70885000
H -4.49934000 19.50700000 7.38700000
O 16.99780000 3.63080000 11.63730000
H 17.31810000 3.30370000 12.51290000
H 17.08550000 4.57727000 11.68570000
O 6.48319000 10.92640000 3.40103000
H 5.63763000 10.99260000 2.85932000
H 6.32943000 11.66920000 4.05745000
O 4.49883000 7.78536000 13.14360000
H 3.57897000 7.92015000 12.89780000
H 4.40785000 7.58602000 14.13040000
O 22.38010000 -2.58249000 19.81360000
H 22.50240000 -2.69404000 18.86030000
H 23.27000000 -2.40926000 20.13910000
O 7.60549000 20.11670000 1.52318000
H 7.40576000 20.44790000 2.44617000
H 6.89393000 19.54060000 1.16986000
O -9.80851000 15.18180000 5.82547000
H -8.98420000 15.49160000 6.23087000
H -9.53181000 14.61140000 5.10912000
O 1.89878000 3.39080000 9.95468000
H 2.53193000 3.64424000 10.67490000
H 2.50766000 2.74971000 9.48478000
O 9.54794000 13.95950000 -3.42910000
H 8.80101000 13.54530000 -2.93926000
H 10.27920000 13.98280000 -2.79475000
O -3.92553000 -7.74473000 -9.61729000
H -3.86420000 -6.99741000 -8.98787000
H -3.75817000 -8.54726000 -9.04343000
O 0.44092400 10.71270000 12.05860000
H 0.83931100 9.90666000 12.42840000
H 0.35887600 10.57580000 11.11070000
O 12.12770000 6.70656000 3.12351000
H 11.86650000 6.60036000 4.03971000
H 11.35860000 6.98032000 2.60943000
O 26.48550000 8.45521000 -11.83130000
H 26.44880000 8.37541000 -10.86690000
H 25.91510000 7.70284000 -12.22220000
O 8.76592000 6.61206000 4.83144000
H 9.47545000 7.33383000 4.67087000
H 7.97602000 7.11755000 4.74199000
O -4.62340000 16.47700000 6.19871000
H -3.86477000 16.34810000 6.84392000
H -4.32266000 17.35310000 5.96346000
O 13.40050000 0.18383700 -4.66240000
H 14.24900000 0.62818200 -4.64394000
H 12.99810000 0.46925800 -5.47766000
O 4.34821000 7.30338000 3.29043000
H 3.53537000 7.62296000 3.70824000
H 4.26206000 6.32561000 3.39956000
O 19.59060000 -11.77550000 22.11300000
H 19.37940000 -11.65870000 23.07720000
H 19.24050000 -10.92170000 21.74620000
O 17.85470000 27.31360000 4.30373000
H 18.81710000 27.23910000 4.30668000
H 17.60460000 26.48810000 4.71635000
O 5.08473000 10.32813000 -12.52020000
H 5.79642000 10.97490000 -12.70560000
H 5.42957000 9.50102000 -12.19440000
O 4.02996000 20.02544000 -3.68524000
H 4.27003000 19.36768000 -4.41446000
H 4.77458000 20.70889000 -3.67208000
O -12.22650000 31.31220000 12.25860000
H -13.08640000 31.61600000 12.55890000
H -12.45140000 31.01980000 11.38190000
O 1.36461000 17.10750000 16.07470000
H 0.77220000 17.75188000 15.70890000
H 0.92914100 16.50655000 16.73020000
O -3.42189000 16.41641000 12.57220000
H -2.47934000 16.47205000 12.38670000
H -3.57781000 16.60579000 13.49940000
O -12.45810000 17.68919000 -2.78020000
H -11.67410000 17.04996000 -2.73820000
H -13.15450000 17.25976000 -3.28648000
O -12.87460000 26.64990000 10.97500000
H -12.40850000 25.90810000 11.29410000
H -12.17020000 27.34310000 10.79950000
O 16.21660000 23.34160000 2.09580000
H 15.77660000 22.57310000 2.51719000
H 16.70080000 22.96170000 1.31827000
O 10.17360000 25.34940000 3.09205000
H 10.77360000 25.05090000 3.73607000
H 9.56179000 24.55220000 2.94022000
O 7.77605000 17.40249000 -2.29855000
H 8.07228000 17.02108000 -1.47153000
H 6.99714000 17.82667000 -2.04022000
O 14.58280000 14.85851000 2.33359000
H 13.85980000 14.80788000 1.70989000
H 14.40950000 15.78726000 2.62982000
O -6.66475000 22.13604000 21.73940000
H -6.28708000 22.59630000 20.91490000
H -5.93185000 21.79252000 22.22920000
O 1.00730000 20.04893000 -3.71784000
H 1.95122000 20.12838000 -3.49284000
H 0.78209600 19.12589000 -3.39957000
O -2.00359000 21.23074000 4.64761000
H -2.33580000 21.82933000 3.96349000
H -1.02960000 21.30476000 4.52207000
O 12.25120000 22.65950000 9.10079000
H 12.66980000 23.37650000 8.67039000
H 12.71420000 21.80279000 8.86118000
O -0.24441300 28.26350000 13.25460000
H 0.00470418 29.04720000 12.75340000
H -0.42289200 27.64620000 12.52560000
O 17.25350000 6.20561000 18.95750000
H 17.26240000 5.23989000 18.92590000
H 16.44650000 6.45123000 18.52680000
O 1.57232000 20.97013000 3.48535000
H 1.78649000 20.77399000 4.39869000
H 0.91633100 20.21040000 3.31060000
O 9.73687000 22.95440000 -12.33500000
H 9.33966000 23.85670000 -12.37350000
H 10.71770000 23.04350000 -12.51260000
O 3.19186000 10.71601000 6.69143000
H 2.64086000 11.50919100 6.70096000
H 3.29213000 10.42658000 7.64120000
O 5.01286000 15.83556000 7.01155000
H 5.25540000 15.72163000 7.95343000
H 5.82406000 15.89012000 6.52344000
O 16.21560000 11.55336300 -2.85105000
H 17.05640000 10.95990000 -2.85949000
H 15.82330000 11.47423400 -1.95729000
O 18.39240000 27.87230000 9.37872000
H 17.70210000 27.92250000 10.02480000
H 19.03780000 28.59220000 9.59111000
O 12.43740000 3.70127000 -6.55194000
H 13.24560000 3.23786000 -6.25474000
H 11.74430000 2.98413000 -6.45253000
O -2.12768000 7.61368000 13.63960000
H -1.96159000 8.49737000 13.27480000
H -3.14023000 7.49516000 13.58540000
O 5.47695000 -0.36860000 18.12950000
H 6.04948000 -0.50690000 18.88220000
H 4.74944000 -1.01130000 18.27050000
O 9.06380000 25.94840000 12.92680000
H 9.07296000 26.84340000 12.64730000
H 9.45712000 25.91050000 13.82760000
O 6.54756000 0.46030000 12.30680000
H 5.90365000 0.93260000 12.84550000
H 7.46789000 0.75900000 12.52900000
O -3.25432000 10.22484000 15.18330000
H -3.21119000 10.03055000 14.23370000
H -4.18029000 10.41526000 15.49650000
O 19.16120000 33.06270000 4.17300000
H 18.25990000 32.84210000 4.00929000
H 19.17920000 34.01520000 4.03867000
O 16.62900000 17.17612000 3.61679000
H 17.41650000 16.71018000 3.87743000
H 15.90400000 16.85119000 4.10392000
O -14.46860000 13.68932000 -6.11559000
H -13.90410000 13.00683400 -6.47353000
H -14.78000000 13.54280000 -5.19548000
O 2.42216000 25.33590000 4.10819000
H 1.94154000 25.93420000 3.57260000
H 2.79204000 24.76870000 3.45220000
O 17.30610000 27.05850000 1.43734000
H 16.39600000 27.17480000 1.71543000
H 17.68200000 27.09460000 2.27931000
O 3.42317000 26.29760000 8.37823000
H 4.23289000 26.39000000 7.90067000
H 3.52615000 25.56510000 8.97284000
O -4.91217000 36.03640000 20.12810000
H -4.02584000 35.89060000 19.75830000
H -4.83605000 36.56340000 20.97620000
O 20.30930000 21.17645000 -1.55507000
H 21.03230000 21.80853000 -1.42815000
H 20.44040000 20.54820000 -2.30558000
O 17.63270000 18.63611000 -1.45303000
H 17.18510000 19.07247000 -2.18930000
H 17.59600000 19.36078000 -0.77261100
O -11.85650000 9.32659000 6.38427000
H -11.91870000 8.42008000 6.73656000
H -10.96070000 9.59879000 6.56532000
O 14.34220000 6.58714000 5.75163000
H 14.20720000 5.87450000 5.13562000
H 14.05200000 6.34737000 6.61672000
O -1.09658000 18.62827000 5.58702000
H -2.04987000 18.36104000 5.55373000
H -0.68250800 17.83005000 5.89003000
O 8.13776000 14.42921000 16.81480000
H 8.76908000 13.89311000 16.26590000
H 8.54559000 14.80894000 17.63890000
O 0.05855300 11.82598900 17.45460000
H 0.96214000 12.10548300 17.27990000
H 0.08964280 10.90732000 17.79400000
O 27.67350000 13.02729400 12.06030000
H 28.34720000 13.18104100 12.69580000
H 27.19490000 12.26214300 12.45860000
O 9.89541000 4.46244000 8.17742000
H 9.62897000 3.62034000 8.58524000
H 9.28724000 5.05986000 8.69880000
O -4.41370000 32.60930000 8.05879000
H -3.76101000 33.23950000 7.70885000
H -4.49934000 31.92700000 7.38700000
O 16.99780000 16.05080000 11.63730000
H 17.31810000 15.72370000 12.51290000
H 17.08550000 16.99727000 11.68570000
O 6.48319000 23.34640000 3.40103000
H 5.63763000 23.41260000 2.85932000
H 6.32943000 24.08920000 4.05745000
O 4.49883000 20.20536000 13.14360000
H 3.57897000 20.34015000 12.89780000
H 4.40785000 20.00602000 14.13040000
O 22.38010000 9.83751000 19.81360000
H 22.50240000 9.72596000 18.86030000
H 23.27000000 10.01074000 20.13910000
O 7.60549000 32.53670000 1.52318000
H 7.40576000 32.86790000 2.44617000
H 6.89393000 31.96060000 1.16986000
O -9.80851000 27.60180000 5.82547000
H -8.98420000 27.91160000 6.23087000
H -9.53181000 27.03140000 5.10912000
O 1.89878000 15.81080000 9.95468000
H 2.53193000 16.06424000 10.67490000
H 2.50766000 15.16971000 9.48478000
O 9.54794000 26.37950000 -3.42910000
H 8.80101000 25.96530000 -2.93926000
H 10.27920000 26.40280000 -2.79475000
O -3.92553000 4.67527000 -9.61729000
H -3.86420000 5.42259000 -8.98787000
H -3.75817000 3.87274000 -9.04343000
O 0.44092400 23.13270000 12.05860000
H 0.83931100 22.32666000 12.42840000
H 0.35887600 22.99580000 11.11070000
O 12.12770000 19.12656000 3.12351000
H 11.86650000 19.02036000 4.03971000
H 11.35860000 19.40032000 2.60943000
O 26.48550000 20.87521000 -11.83130000
H 26.44880000 20.79541000 -10.86690000
H 25.91510000 20.12284000 -12.22220000
O 8.76592000 19.03206000 4.83144000
H 9.47545000 19.75383000 4.67087000
H 7.97602000 19.53755000 4.74199000
O -4.62340000 28.89700000 6.19871000
H -3.86477000 28.76810000 6.84392000
H -4.32266000 29.77310000 5.96346000
O 13.40050000 12.60383700 -4.66240000
H 14.24900000 13.04818200 -4.64394000
H 12.99810000 12.88925800 -5.47766000
O 4.34821000 19.72338000 3.29043000
H 3.53537000 20.04296000 3.70824000
H 4.26206000 18.74561000 3.39956000
O 19.59060000 0.64450000 22.11300000
H 19.37940000 0.76130000 23.07720000
H 19.24050000 1.49830000 21.74620000
O 30.27470000 14.89360000 4.30373000
H 31.23710000 14.81910000 4.30668000
H 30.02460000 14.06810000 4.71635000
O 17.50473000 -2.09187000 -12.52020000
H 18.21642000 -1.44510000 -12.70560000
H 17.84957000 -2.91898000 -12.19440000
O 16.44996000 7.60544000 -3.68524000
H 16.69003000 6.94768000 -4.41446000
H 17.19458000 8.28889000 -3.67208000
O 0.19350000 18.89220000 12.25860000
H -0.66640000 19.19600000 12.55890000
H -0.03140000 18.59980000 11.38190000
O 13.78461000 4.68750000 16.07470000
H 13.19220000 5.33188000 15.70890000
H 13.34914100 4.08655000 16.73020000
O 8.99811000 3.99641000 12.57220000
H 9.94066000 4.05205000 12.38670000
H 8.84219000 4.18579000 13.49940000
O -0.03810000 5.26919000 -2.78020000
H 0.74590000 4.62996000 -2.73820000
H -0.73450000 4.83976000 -3.28648000
O -0.45460000 14.22990000 10.97500000
H 0.01150000 13.48810000 11.29410000
H 0.24980000 14.92310000 10.79950000
O 28.63660000 10.92160000 2.09580000
H 28.19660000 10.15310000 2.51719000
H 29.12080000 10.54170000 1.31827000
O 22.59360000 12.92940000 3.09205000
H 23.19360000 12.63090000 3.73607000
H 21.98179000 12.13220000 2.94022000
O 20.19605000 4.98249000 -2.29855000
H 20.49228000 4.60108000 -1.47153000
H 19.41714000 5.40667000 -2.04022000
O 27.00280000 2.43851000 2.33359000
H 26.27980000 2.38788000 1.70989000
H 26.82950000 3.36726000 2.62982000
O 5.75525000 9.71604000 21.73940000
H 6.13292000 10.17630000 20.91490000
H 6.48815000 9.37252000 22.22920000
O 13.42730000 7.62893000 -3.71784000
H 14.37122000 7.70838000 -3.49284000
H 13.20209600 6.70589000 -3.39957000
O 10.41641000 8.81074000 4.64761000
H 10.08420000 9.40933000 3.96349000
H 11.39040000 8.88476000 4.52207000
O 24.67120000 10.23950000 9.10079000
H 25.08980000 10.95650000 8.67039000
H 25.13420000 9.38279000 8.86118000
O 12.17558700 15.84350000 13.25460000
H 12.42470418 16.62720000 12.75340000
H 11.99710800 15.22620000 12.52560000
O 29.67350000 -6.21439000 18.95750000
H 29.68240000 -7.18011000 18.92590000
H 28.86650000 -5.96877000 18.52680000
O 13.99232000 8.55013000 3.48535000
H 14.20649000 8.35399000 4.39869000
H 13.33633100 7.79040000 3.31060000
O 22.15687000 10.53440000 -12.33500000
H 21.75966000 11.43670000 -12.37350000
H 23.13770000 10.62350000 -12.51260000
O 15.61186000 -1.70399000 6.69143000
H 15.06086000 -0.91080900 6.70096000
H 15.71213000 -1.99342000 7.64120000
O 17.43286000 3.41556000 7.01155000
H 17.67540000 3.30163000 7.95343000
H 18.24406000 3.47012000 6.52344000
O 28.63560000 -0.86663700 -2.85105000
H 29.47640000 -1.46010000 -2.85949000
H 28.24330000 -0.94576600 -1.95729000
O 30.81240000 15.45230000 9.37872000
H 30.12210000 15.50250000 10.02480000
H 31.45780000 16.17220000 9.59111000
O 24.85740000 -8.71873000 -6.55194000
H 25.66560000 -9.18214000 -6.25474000
H 24.16430000 -9.43587000 -6.45253000
O 10.29232000 -4.80632000 13.63960000
H 10.45841000 -3.92263000 13.27480000
H 9.27977000 -4.92484000 13.58540000
O 17.89695000 -12.78860000 18.12950000
H 18.46948000 -12.92690000 18.88220000
H 17.16944000 -13.43130000 18.27050000
O 21.48380000 13.52840000 12.92680000
H 21.49296000 14.42340000 12.64730000
H 21.87712000 13.49050000 13.82760000
O 18.96756000 -11.95970000 12.30680000
H 18.32365000 -11.48740000 12.84550000
H 19.88789000 -11.66100000 12.52900000
O 9.16568000 -2.19516000 15.18330000
H 9.20881000 -2.38945000 14.23370000
H 8.23971000 -2.00474000 15.49650000
O 31.58120000 20.64270000 4.17300000
H 30.67990000 20.42210000 4.00929000
H 31.59920000 21.59520000 4.03867000
O 29.04900000 4.75612000 3.61679000
H 29.83650000 4.29018000 3.87743000
H 28.32400000 4.43119000 4.10392000
O -2.04860000 1.26932000 -6.11559000
H -1.48410000 0.58683400 -6.47353000
H -2.36000000 1.12280000 -5.19548000
O 14.84216000 12.91590000 4.10819000
H 14.36154000 13.51420000 3.57260000
H 15.21204000 12.34870000 3.45220000
O 29.72610000 14.63850000 1.43734000
H 28.81600000 14.75480000 1.71543000
H 30.10200000 14.67460000 2.27931000
O 15.84317000 13.87760000 8.37823000
H 16.65289000 13.97000000 7.90067000
H 15.94615000 13.14510000 8.97284000
O 7.50783000 23.61640000 20.12810000
H 8.39416000 23.47060000 19.75830000
H 7.58395000 24.14340000 20.97620000
O 32.72930000 8.75645000 -1.55507000
H 33.45230000 9.38853000 -1.42815000
H 32.86040000 8.12820000 -2.30558000
O 30.05270000 6.21611000 -1.45303000
H 29.60510000 6.65247000 -2.18930000
H 30.01600000 6.94078000 -0.77261100
O 0.56350000 -3.09341000 6.38427000
H 0.50130000 -3.99992000 6.73656000
H 1.45930000 -2.82121000 6.56532000
O 26.76220000 -5.83286000 5.75163000
H 26.62720000 -6.54550000 5.13562000
H 26.47200000 -6.07263000 6.61672000
O 11.32342000 6.20827000 5.58702000
H 10.37013000 5.94104000 5.55373000
H 11.73749200 5.41005000 5.89003000
O 20.55776000 2.00921000 16.81480000
H 21.18908000 1.47311000 16.26590000
H 20.96559000 2.38894000 17.63890000
O 12.47855300 -0.59401100 17.45460000
H 13.38214000 -0.31451700 17.27990000
H 12.50964280 -1.51268000 17.79400000
O 40.09350000 0.60729400 12.06030000
H 40.76720000 0.76104100 12.69580000
H 39.61490000 -0.15785700 12.45860000
O 22.31541000 -7.95756000 8.17742000
H 22.04897000 -8.79966000 8.58524000
H 21.70724000 -7.36014000 8.69880000
O 8.00630000 20.18930000 8.05879000
H 8.65899000 20.81950000 7.70885000
H 7.92066000 19.50700000 7.38700000
O 29.41780000 3.63080000 11.63730000
H 29.73810000 3.30370000 12.51290000
H 29.50550000 4.57727000 11.68570000
O 18.90319000 10.92640000 3.40103000
H 18.05763000 10.99260000 2.85932000
H 18.74943000 11.66920000 4.05745000
O 16.91883000 7.78536000 13.14360000
H 15.99897000 7.92015000 12.89780000
H 16.82785000 7.58602000 14.13040000
O 34.80010000 -2.58249000 19.81360000
H 34.92240000 -2.69404000 18.86030000
H 35.69000000 -2.40926000 20.13910000
O 20.02549000 20.11670000 1.52318000
H 19.82576000 20.44790000 2.44617000
H 19.31393000 19.54060000 1.16986000
O 2.61149000 15.18180000 5.82547000
H 3.43580000 15.49160000 6.23087000
H 2.88819000 14.61140000 5.10912000
O 14.31878000 3.39080000 9.95468000
H 14.95193000 3.64424000 10.67490000
H 14.92766000 2.74971000 9.48478000
O 21.96794000 13.95950000 -3.42910000
H 21.22101000 13.54530000 -2.93926000
H 22.69920000 13.98280000 -2.79475000
O 8.49447000 -7.74473000 -9.61729000
H 8.55580000 -6.99741000 -8.98787000
H 8.66183000 -8.54726000 -9.04343000
O 12.86092400 10.71270000 12.05860000
H 13.25931100 9.90666000 12.42840000
H 12.77887600 10.57580000 11.11070000
O 24.54770000 6.70656000 3.12351000
H 24.28650000 6.60036000 4.03971000
H 23.77860000 6.98032000 2.60943000
O 38.90550000 8.45521000 -11.83130000
H 38.86880000 8.37541000 -10.86690000
H 38.33510000 7.70284000 -12.22220000
O 21.18592000 6.61206000 4.83144000
H 21.89545000 7.33383000 4.67087000
H 20.39602000 7.11755000 4.74199000
O 7.79660000 16.47700000 6.19871000
H 8.55523000 16.34810000 6.84392000
H 8.09734000 17.35310000 5.96346000
O 25.82050000 0.18383700 -4.66240000
H 26.66900000 0.62818200 -4.64394000
H 25.41810000 0.46925800 -5.47766000
O 16.76821000 7.30338000 3.29043000
H 15.95537000 7.62296000 3.70824000
H 16.68206000 6.32561000 3.39956000
O 32.01060000 -11.77550000 22.11300000
H 31.79940000 -11.65870000 23.07720000
H 31.66050000 -10.92170000 21.74620000
O 30.27470000 27.31360000 4.30373000
H 31.23710000 27.23910000 4.30668000
H 30.02460000 26.48810000 4.71635000
O 17.50473000 10.32813000 -12.52020000
H 18.21642000 10.97490000 -12.70560000
H 17.84957000 9.50102000 -12.19440000
O 16.44996000 20.02544000 -3.68524000
H 16.69003000 19.36768000 -4.41446000
H 17.19458000 20.70889000 -3.67208000
O 0.19350000 31.31220000 12.25860000
H -0.66640000 31.61600000 12.55890000
H -0.03140000 31.01980000 11.38190000
O 13.78461000 17.10750000 16.07470000
H 13.19220000 17.75188000 15.70890000
H 13.34914100 16.50655000 16.73020000
O 8.99811000 16.41641000 12.57220000
H 9.94066000 16.47205000 12.38670000
H 8.84219000 16.60579000 13.49940000
O -0.03810000 17.68919000 -2.78020000
H 0.74590000 17.04996000 -2.73820000
H -0.73450000 17.25976000 -3.28648000
O -0.45460000 26.64990000 10.97500000
H 0.01150000 25.90810000 11.29410000
H 0.24980000 27.34310000 10.79950000
O 28.63660000 23.34160000 2.09580000
H 28.19660000 22.57310000 2.51719000
H 29.12080000 22.96170000 1.31827000
O 22.59360000 25.34940000 3.09205000
H 23.19360000 25.05090000 3.73607000
H 21.98179000 24.55220000 2.94022000
O 20.19605000 17.40249000 -2.29855000
H 20.49228000 17.02108000 -1.47153000
H 19.41714000 17.82667000 -2.04022000
O 27.00280000 14.85851000 2.33359000
H 26.27980000 14.80788000 1.70989000
H 26.82950000 15.78726000 2.62982000
O 5.75525000 22.13604000 21.73940000
H 6.13292000 22.59630000 20.91490000
H 6.48815000 21.79252000 22.22920000
O 13.42730000 20.04893000 -3.71784000
H 14.37122000 20.12838000 -3.49284000
H 13.20209600 19.12589000 -3.39957000
O 10.41641000 21.23074000 4.64761000
H 10.08420000 21.82933000 3.96349000
H 11.39040000 21.30476000 4.52207000
O 24.67120000 22.65950000 9.10079000
H 25.08980000 23.37650000 8.67039000
H 25.13420000 21.80279000 8.86118000
O 12.17558700 28.26350000 13.25460000
H 12.42470418 29.04720000 12.75340000
H 11.99710800 27.64620000 12.52560000
O 29.67350000 6.20561000 18.95750000
H 29.68240000 5.23989000 18.92590000
H 28.86650000 6.45123000 18.52680000
O 13.99232000 20.97013000 3.48535000
H 14.20649000 20.77399000 4.39869000
H 13.33633100 20.21040000 3.31060000
O 22.15687000 22.95440000 -12.33500000
H 21.75966000 23.85670000 -12.37350000
H 23.13770000 23.04350000 -12.51260000
O 15.61186000 10.71601000 6.69143000
H 15.06086000 11.50919100 6.70096000
H 15.71213000 10.42658000 7.64120000
O 17.43286000 15.83556000 7.01155000
H 17.67540000 15.72163000 7.95343000
H 18.24406000 15.89012000 6.52344000
O 28.63560000 11.55336300 -2.85105000
H 29.47640000 10.95990000 -2.85949000
H 28.24330000 11.47423400 -1.95729000
O 30.81240000 27.87230000 9.37872000
H 30.12210000 27.92250000 10.02480000
H 31.45780000 28.59220000 9.59111000
O 24.85740000 3.70127000 -6.55194000
H 25.66560000 3.23786000 -6.25474000
H 24.16430000 2.98413000 -6.45253000
O 10.29232000 7.61368000 13.63960000
H 10.45841000 8.49737000 13.27480000
H 9.27977000 7.49516000 13.58540000
O 17.89695000 -0.36860000 18.12950000
H 18.46948000 -0.50690000 18.88220000
H 17.16944000 -1.01130000 18.27050000
O 21.48380000 25.94840000 12.92680000
H 21.49296000 26.84340000 12.64730000
H 21.87712000 25.91050000 13.82760000
O 18.96756000 0.46030000 12.30680000
H 18.32365000 0.93260000 12.84550000
H 19.88789000 0.75900000 12.52900000
O 9.16568000 10.22484000 15.18330000
H 9.20881000 10.03055000 14.23370000
H 8.23971000 10.41526000 15.49650000
O 31.58120000 33.06270000 4.17300000
H 30.67990000 32.84210000 4.00929000
H 31.59920000 34.01520000 4.03867000
O 29.04900000 17.17612000 3.61679000
H 29.83650000 16.71018000 3.87743000
H 28.32400000 16.85119000 4.10392000
O -2.04860000 13.68932000 -6.11559000
H -1.48410000 13.00683400 -6.47353000
H -2.36000000 13.54280000 -5.19548000
O 14.84216000 25.33590000 4.10819000
H 14.36154000 25.93420000 3.57260000
H 15.21204000 24.76870000 3.45220000
O 29.72610000 27.05850000 1.43734000
H 28.81600000 27.17480000 1.71543000
H 30.10200000 27.09460000 2.27931000
O 15.84317000 26.29760000 8.37823000
H 16.65289000 26.39000000 7.90067000
H 15.94615000 25.56510000 8.97284000
O 7.50783000 36.03640000 20.12810000
H 8.39416000 35.89060000 19.75830000
H 7.58395000 36.56340000 20.97620000
O 32.72930000 21.17645000 -1.55507000
H 33.45230000 21.80853000 -1.42815000
H 32.86040000 20.54820000 -2.30558000
O 30.05270000 18.63611000 -1.45303000
H 29.60510000 19.07247000 -2.18930000
H 30.01600000 19.36078000 -0.77261100
O 0.56350000 9.32659000 6.38427000
H 0.50130000 8.42008000 6.73656000
H 1.45930000 9.59879000 6.56532000
O 26.76220000 6.58714000 5.75163000
H 26.62720000 5.87450000 5.13562000
H 26.47200000 6.34737000 6.61672000
O 11.32342000 18.62827000 5.58702000
H 10.37013000 18.36104000 5.55373000
H 11.73749200 17.83005000 5.89003000
O 20.55776000 14.42921000 16.81480000
H 21.18908000 13.89311000 16.26590000
H 20.96559000 14.80894000 17.63890000
O 12.47855300 11.82598900 17.45460000
H 13.38214000 12.10548300 17.27990000
H 12.50964280 10.90732000 17.79400000
O 40.09350000 13.02729400 12.06030000
H 40.76720000 13.18104100 12.69580000
H 39.61490000 12.26214300 12.45860000
O 22.31541000 4.46244000 8.17742000
H 22.04897000 3.62034000 8.58524000
H 21.70724000 5.05986000 8.69880000
O 8.00630000 32.60930000 8.05879000
H 8.65899000 33.23950000 7.70885000
H 7.92066000 31.92700000 7.38700000
O 29.41780000 16.05080000 11.63730000
H 29.73810000 15.72370000 12.51290000
H 29.50550000 16.99727000 11.68570000
O 18.90319000 23.34640000 3.40103000
H 18.05763000 23.41260000 2.85932000
H 18.74943000 24.08920000 4.05745000
O 16.91883000 20.20536000 13.14360000
H 15.99897000 20.34015000 12.89780000
H 16.82785000 20.00602000 14.13040000
O 34.80010000 9.83751000 19.81360000
H 34.92240000 9.72596000 18.86030000
H 35.69000000 10.01074000 20.13910000
O 20.02549000 32.53670000 1.52318000
H 19.82576000 32.86790000 2.44617000
H 19.31393000 31.96060000 1.16986000
O 2.61149000 27.60180000 5.82547000
H 3.43580000 27.91160000 6.23087000
H 2.88819000 27.03140000 5.10912000
O 14.31878000 15.81080000 9.95468000
H 14.95193000 16.06424000 10.67490000
H 14.92766000 15.16971000 9.48478000
O 21.96794000 26.37950000 -3.42910000
H 21.22101000 25.96530000 -2.93926000
H 22.69920000 26.40280000 -2.79475000
O 8.49447000 4.67527000 -9.61729000
H 8.55580000 5.42259000 -8.98787000
H 8.66183000 3.87274000 -9.04343000
O 12.86092400 23.13270000 12.05860000
H 13.25931100 22.32666000 12.42840000
H 12.77887600 22.99580000 11.11070000
O 24.54770000 19.12656000 3.12351000
H 24.28650000 19.02036000 4.03971000
H 23.77860000 19.40032000 2.60943000
O 38.90550000 20.87521000 -11.83130000
H 38.86880000 20.79541000 -10.86690000
H 38.33510000 20.12284000 -12.22220000
O 21.18592000 19.03206000 4.83144000
H 21.89545000 19.75383000 4.67087000
H 20.39602000 19.53755000 4.74199000
O 7.79660000 28.89700000 6.19871000
H 8.55523000 28.76810000 6.84392000
H 8.09734000 29.77310000 5.96346000
O 25.82050000 12.60383700 -4.66240000
H 26.66900000 13.04818200 -4.64394000
H 25.41810000 12.88925800 -5.47766000
O 16.76821000 19.72338000 3.29043000
H 15.95537000 20.04296000 3.70824000
H 16.68206000 18.74561000 3.39956000
O 32.01060000 0.64450000 22.11300000
H 31.79940000 0.76130000 23.07720000
H 31.66050000 1.49830000 21.74620000
&END COORD
&END SUBSYS
&END FORCE_EVAL

14
tests/NNP/regtest-1/TEST_FILES.toml
View file

@ -8,4 +8,18 @@
"H2O-64_C-NNP_MD-NpT-numeric.inp" = [{matcher="M002", tol=1.0E-14, ref=0.375124841615E+04}]
"H2O-64_C-NNP_biased_MD.inp" = [{matcher="M002", tol=1.0E-14, ref=0.375124232248E+04}]
"H2O-64_C-NNP_biased_MD_restart.inp" = [{matcher="M002", tol=1.0E-14, ref=0.375124124406E+04}]
# H2O-256_C-NNP_MD-bitexact: non-cubic 256-water NVT (24.84 x 24.84
# x 12.42 A, two axes with L > 2*r_cut). Locks the step-zero
# energy, final-step energy and conserved quantity at 1e-14
# relative. Runtime ~3 s on one rank. The bit-exact conserved
# quantity check may need relaxing in future (Lyapunov-divergent FP
# across compilers), but without it integrator and neighbour-cache
# bugs that let energy drift and blow up the temperature can pass
# silently.
"H2O-256_C-NNP_MD-bitexact.inp" = [
{matcher="M_INIT_ENERGY", tol=1.0E-14, ref=15004.992479346592518},
{matcher="M002", tol=1.0E-14, ref=0.150049870566E+05},
{matcher="M_CONS_QTY", tol=1.0E-14, ref=0.150060861414E+05},
]
#EOF

11
tests/matcher_classes.py
View file

@ -25,11 +25,17 @@ class Matcher(Protocol):
# ======================================================================================
class GenericMatcher(Matcher):
def __init__(
self, pattern: str, col: int, regex: bool = False, abs_value: bool = False
self,
pattern: str,
col: int,
regex: bool = False,
abs_value: bool = False,
first: bool = False,
):
self.pattern = pattern
self.regex_mode = regex
self.abs_value = abs_value
self.first = first
if not regex:
for c in r"[]()|+*?":
pattern = pattern.replace(c, f"\\{c}")
@ -41,7 +47,8 @@ class GenericMatcher(Matcher):
assert isinstance(tol, float) or isinstance(ref, int)
assert isinstance(ref, float) or isinstance(ref, int)
# grep result
for line in reversed(output.split("\n")):
lines = output.split("\n")
for line in lines if self.first else reversed(lines):
match = self.regex.search(line)
if match:
if self.regex_mode and match.groups():

9
tests/matchers.py
View file

@ -307,4 +307,13 @@ registry["E_RIRS_LUMO"] = GenericMatcher(r"G0W0 conduction band minimum", col=6)
# Floquet Calculations
registry["Quasienergy"] = GenericMatcher(r" 4", col=2)
registry["Floquet_DOS"] = GenericMatcher(r"-1.690", col=2)
# NNP MD matchers. M_INIT_ENERGY passes first=True because ENERGY|Total
# FORCE_EVAL is printed once per MD step, and the default (last-line)
# strategy would return step N rather than the step-0 initial value.
# M_CONS_QTY reads the final-step MD|Conserved quantity.
registry["M_INIT_ENERGY"] = GenericMatcher(
r"ENERGY| Total FORCE_EVAL", col=9, first=True
)
registry["M_CONS_QTY"] = GenericMatcher(r"MD| Conserved quantity", col=5)
# EOF