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

Helium-Solute NNP for interaction [cschran/lduran/hforbert] (#3043)

Browse source 
Helium-Solute NNP for interaction [cschran/lduran/hforbert]
This commit is contained in:
hforbert 2023-10-17 16:52:59 +02:00 • committed by GitHub
parent 54ae6c8529
commit 715739a36e
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
21 changed files with 4626 additions and 110 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/CMakeLists.txt
View file

@ -1036,6 +1036,7 @@ list(
motion/helium_interactions.F
motion/helium_io.F
motion/helium_methods.F
motion/helium_nnp.F
motion/helium_sampling.F
motion/helium_types.F
motion/helium_worm.F

4
src/f77_interface.F
View file

@ -109,8 +109,7 @@ MODULE f77_interface
mp_perf_env_type,&
rm_mp_perf_env
USE nnp_environment, ONLY: nnp_init
USE nnp_environment_types, ONLY: nnp_env_create,&
nnp_type
USE nnp_environment_types, ONLY: nnp_type
USE offload_api, ONLY: offload_get_chosen_device,&
offload_get_device_count,&
offload_set_chosen_device
@ -851,7 +850,6 @@ CONTAINS
CASE (do_nnp)
ALLOCATE (nnp_env)
CALL nnp_env_create(nnp_env)
CALL nnp_init(nnp_env, root_section, my_para_env, force_env_section=force_env_section, &
subsys_section=subsys_section, use_motion_section=use_motion_section)
CALL force_env_create(my_force_env, root_section, nnp_env=nnp_env, para_env=my_para_env, &

196
src/motion/helium_interactions.F
View file

@ -12,10 +12,13 @@
! **************************************************************************************************
MODULE helium_interactions
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE helium_common, ONLY: helium_eval_chain,&
helium_eval_expansion,&
helium_pbc,&
helium_spline
USE helium_nnp, ONLY: helium_nnp_print
USE helium_types, ONLY: e_id_interact,&
e_id_kinetic,&
e_id_potential,&
@ -26,8 +29,15 @@ MODULE helium_interactions
helium_solvent_type
USE input_constants, ONLY: helium_sampling_worm,&
helium_solute_intpot_mwater,&
helium_solute_intpot_nnp,&
helium_solute_intpot_none
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_environment_types, ONLY: nnp_type
USE nnp_model, ONLY: nnp_gradients,&
nnp_predict
USE physcon, ONLY: angstrom,&
kelvin
USE pint_types, ONLY: pint_env_type
@ -334,6 +344,7 @@ CONTAINS
!> \param force calculated force (if requested)
!> \par History
!> 2019-09 Added multiple-time striding in imag. time [cschran]
!> 2023-07-23 Modified to work with NNP solute-solvent interactions [lduran]
!> \author Lukasz Walewski
! **************************************************************************************************
SUBROUTINE helium_bead_solute_e_f(pint_env, helium, helium_part_index, &
@ -386,6 +397,26 @@ CONTAINS
)
END IF
CASE (helium_solute_intpot_nnp)
IF (PRESENT(force)) THEN
force(:, :) = 0.0_dp
my_force => force(qi, :)
CALL helium_intpot_nnp( &
pint_env%x(qi, :), &
helium, &
helium_r, &
energy, &
my_force &
)
ELSE
CALL helium_intpot_nnp( &
pint_env%x(qi, :), &
helium, &
helium_r, &
energy &
)
END IF
CASE (helium_solute_intpot_none)
energy = 0.0_dp
IF (PRESENT(force)) THEN
@ -629,6 +660,171 @@ CONTAINS
END SUBROUTINE helium_intpot_model_water
! ***************************************************************************
!> \brief Calculate helium-solute interaction energy and forces between one
!> helium bead and the corresponding solute time slice using NNP.
!> \param solute_x solute positions ARR(3*NATOMS)
!> to global atom indices
!> \param helium only needed for helium_pbc call at the moment
!> \param helium_x helium bead position ARR(3)
!> \param energy calculated interaction energy
!> \param force (optional) calculated force
!> \date 2023-02-22
!> \author Laura Duran
! **************************************************************************************************
SUBROUTINE helium_intpot_nnp(solute_x, helium, helium_x, energy, force)
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: solute_x
TYPE(helium_solvent_type), INTENT(IN) :: helium
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: helium_x
REAL(KIND=dp), INTENT(OUT) :: energy
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
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), 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
NULLIFY (logger)
logger => cp_get_default_logger()
IF (PRESENT(force)) THEN
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
IF (helium%nnp%ele(i) == 'He') THEN
ind_he = ig
DO m = 1, 3
helium%nnp%coord(m, ig) = helium_x(m)
END DO
ig = ig + 1
END IF
DO j = 1, helium%solute_atoms
IF (helium%nnp%ele(i) == helium%solute_element(j)) THEN
DO m = 1, 3
helium%nnp%coord(m, ig) = solute_x(3*(j - 1) + m)
END DO
ig = ig + 1
END IF
END DO
END DO
! check for hard core condition
IF (ASSOCIATED(helium%nnp_sr_cut)) THEN
DO i = 1, helium%nnp%num_atoms
IF (i == ind_he) CYCLE
rvect(:) = helium%nnp%coord(:, i) - helium%nnp%coord(:, ind_he)
CALL helium_pbc(helium, rvect)
rsqr = rvect(1)*rvect(1) + rvect(2)*rvect(2) + rvect(3)*rvect(3)
IF (rsqr < helium%nnp_sr_cut(helium%nnp%ele_ind(i))) THEN
energy = 0.3_dp + 1.0_dp/rsqr
IF (PRESENT(force)) THEN
force = 0.0_dp
END IF
RETURN
END IF
END DO
END IF
! reset flag if there's an extrapolation to report:
helium%nnp%output_expol = .FALSE.
! calc atomic contribution to energy and force
!NOTE corresponds to nnp_force line with parallelization:
!DO i = istart, istart + mecalc - 1
DO i = 1, helium%nnp%num_atoms
!determine index of atom type
ind = helium%nnp%ele_ind(i)
!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)
ELSE
CALL nnp_calc_acsf(nnp, i)
END IF
! input nodes filled, perform prediction:
DO i_com = 1, helium%nnp%n_committee !loop over committee members
! Predict energy
CALL nnp_predict(helium%nnp%arc(ind), helium%nnp, i_com)
helium%nnp%atomic_energy(i, i_com) = helium%nnp%arc(ind)%layer(helium%nnp%n_layer)%node(1) ! + helium%nnp%atom_energies(ind)
!Gradients
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
END IF
END DO ! end loop over committee members
!deallocate memory
IF (PRESENT(force)) THEN
DEALLOCATE (denergydsym)
DEALLOCATE (dsymdxyz)
END IF
END DO ! end loop over num_atoms
! calculate energy:
helium%nnp%committee_energy(:) = SUM(helium%nnp%atomic_energy, 1)
energy = SUM(helium%nnp%committee_energy)/REAL(helium%nnp%n_committee, dp)
helium%nnp%nnp_potential_energy = energy
IF (PRESENT(force)) THEN
! bring myforce to force array
DO j = 1, helium%nnp%num_atoms
DO k = 1, 3
helium%nnp%committee_forces(k, j, :) = helium%nnp%myforce(k, j, :)
END DO
END DO
helium%nnp%nnp_forces(:, :) = SUM(helium%nnp%committee_forces, DIM=3)/REAL(helium%nnp%n_committee, dp)
! project out helium force entry
ig = 1
DO j = 1, helium%nnp%num_atoms
IF (j == ind_he) CYCLE
DO k = 1, 3
force(3*(helium%nnp%sort(ig) - 1) + k) = helium%nnp%nnp_forces(k, j)
END DO
ig = ig + 1
END DO
END IF
! print properties if requested
print_section => section_vals_get_subs_vals(helium%nnp%nnp_input, "PRINT")
CALL helium_nnp_print(helium%nnp, print_section, ind_he)
RETURN
END SUBROUTINE helium_intpot_nnp
! ***************************************************************************
!> \brief Helium-helium pair interaction potential.
!> \param r ...

6
src/motion/helium_io.F
View file

@ -41,7 +41,7 @@ MODULE helium_io
USE input_constants, ONLY: &
fmt_id_pdb, fmt_id_xyz, helium_cell_shape_cube, helium_cell_shape_octahedron, &
helium_sampling_ceperley, helium_sampling_worm, helium_solute_intpot_mwater, &
helium_solute_intpot_none, perm_cycle, perm_plain
helium_solute_intpot_nnp, helium_solute_intpot_none, perm_cycle, perm_plain
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
@ -179,6 +179,8 @@ CONTAINS
!> \brief Write helium parameters to the output unit
!> \param helium ...
!> \date 2009-06-03
!> \par History
!> 2023-07-23 Modified to work with NNP solute-solvent interactions [lduran]
!> \author Lukasz Walewski
! **************************************************************************************************
SUBROUTINE helium_write_setup(helium)
@ -364,6 +366,8 @@ CONTAINS
WRITE (stmp1, *) "NONE"
CASE (helium_solute_intpot_mwater)
WRITE (stmp1, *) "MWATER"
CASE (helium_solute_intpot_nnp)
WRITE (stmp1, *) "NNP"
CASE DEFAULT
WRITE (stmp1, *) "***UNKNOWN***"
END SELECT

88
src/motion/helium_methods.F
View file

@ -19,8 +19,12 @@ MODULE helium_methods
get_cell
USE cp_files, ONLY: close_file,&
open_file
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_log_handling, ONLY: cp_add_default_logger,&
cp_get_default_logger,&
cp_logger_create,&
cp_logger_release,&
cp_logger_type,&
cp_rm_default_logger
USE cp_output_handling, ONLY: cp_printkey_is_on
USE cp_subsys_types, ONLY: cp_subsys_get,&
cp_subsys_type
@ -33,6 +37,7 @@ MODULE helium_methods
USE helium_interactions, ONLY: helium_vij
USE helium_io, ONLY: helium_write_line,&
helium_write_setup
USE helium_nnp, ONLY: helium_init_nnp
USE helium_sampling, ONLY: helium_sample
USE helium_types, ONLY: helium_solvent_p_type,&
helium_solvent_type,&
@ -46,6 +51,7 @@ MODULE helium_methods
helium_cell_shape_octahedron,&
helium_sampling_ceperley,&
helium_sampling_worm,&
helium_solute_intpot_nnp,&
helium_solute_intpot_none
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
@ -58,7 +64,8 @@ MODULE helium_methods
max_line_length
USE mathconstants, ONLY: pi,&
twopi
USE message_passing, ONLY: mp_comm_type
USE message_passing, ONLY: mp_para_env_type
USE nnp_environment_types, ONLY: nnp_env_release
USE parallel_rng_types, ONLY: GAUSSIAN,&
UNIFORM,&
rng_stream_p_type,&
@ -99,6 +106,7 @@ CONTAINS
!> \param solute ...
!> \par History
!> 2016-07-14 Modified to work with independent helium_env [cschran]
!> 2023-07-23 Modified to work with NNP solute-solvent interactions [lduran]
!> \author hforbert
! **************************************************************************************************
SUBROUTINE helium_create(helium_env, input, solute)
@ -118,9 +126,9 @@ CONTAINS
REAL(KIND=dp) :: cgeof, dx, he_mass, mHe, rtmp, T, tau, &
tcheck, x1, x_spline
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pot_transfer
TYPE(cp_logger_type), POINTER :: logger
TYPE(mp_comm_type) :: new_comm
TYPE(section_vals_type), POINTER :: helium_section, input_worm
TYPE(cp_logger_type), POINTER :: logger, tmplogger
TYPE(mp_para_env_type), POINTER :: new_comm
TYPE(section_vals_type), POINTER :: helium_section, input_worm, nnp_section
CALL timeset(routineN, handle)
@ -163,17 +171,20 @@ CONTAINS
ELSE
color_sub = 1
END IF
ALLOCATE (new_comm)
CALL new_comm%from_split(logger%para_env, color_sub)
! release new_comm for processors without helium_env
IF (mepos == 0) THEN
CALL new_comm%free()
DEALLOCATE (new_comm)
NULLIFY (new_comm)
END IF
NULLIFY (helium_env)
IF (mepos .GT. 0) THEN
ALLOCATE (helium_env(mepos))
DO k = 1, mepos
helium_env(k)%comm = new_comm
helium_env(k)%comm => new_comm
NULLIFY (helium_env(k)%env_all)
helium_env(k)%env_all => env_all
ALLOCATE (helium_env(k)%helium)
@ -338,24 +349,24 @@ CONTAINS
SELECT CASE (helium_env(k)%helium%cell_shape)
CASE (helium_cell_shape_octahedron)
helium_env(k)%helium%cell_m(1, 1) = helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(2, 1) = 0.0_dp
helium_env(k)%helium%cell_m(3, 1) = 0.0_dp
helium_env(k)%helium%cell_m(1, 2) = 0.0_dp
helium_env(k)%helium%cell_m(2, 2) = helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(3, 2) = 0.0_dp
helium_env(k)%helium%cell_m(1, 3) = helium_env(k)%helium%cell_size/2.0_dp
helium_env(k)%helium%cell_m(2, 3) = helium_env(k)%helium%cell_size/2.0_dp
helium_env(k)%helium%cell_m(3, 3) = helium_env(k)%helium%cell_size/2.0_dp
helium_env(k)%helium%cell_m_inv(1, 1) = 1.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(2, 1) = 0.0_dp
helium_env(k)%helium%cell_m_inv(3, 1) = 0.0_dp
helium_env(k)%helium%cell_m_inv(1, 2) = 0.0_dp
helium_env(k)%helium%cell_m_inv(2, 2) = 1.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(3, 2) = 0.0_dp
helium_env(k)%helium%cell_m_inv(1, 3) = -1.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(2, 3) = -1.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(3, 3) = 2.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(1, 1) = -0.5_dp*helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(2, 1) = 0.5_dp*helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(3, 1) = 0.5_dp*helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(1, 2) = 0.5_dp*helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(2, 2) = -0.5_dp*helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(3, 2) = 0.5_dp*helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(1, 3) = 0.5_dp*helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(2, 3) = 0.5_dp*helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m(3, 3) = -0.5_dp*helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(1, 1) = 0.0_dp
helium_env(k)%helium%cell_m_inv(2, 1) = 1.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(3, 1) = 1.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(1, 2) = 1.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(2, 2) = 0.0_dp
helium_env(k)%helium%cell_m_inv(3, 2) = 1.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(1, 3) = 1.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(2, 3) = 1.0_dp/helium_env(k)%helium%cell_size
helium_env(k)%helium%cell_m_inv(3, 3) = 0.0_dp
CASE (helium_cell_shape_cube)
helium_env(k)%helium%cell_m(1, 1) = helium_env(k)%helium%cell_size
@ -760,6 +771,21 @@ CONTAINS
! set the interaction potential type
CALL section_vals_val_get(helium_section, "SOLUTE_INTERACTION", &
i_val=helium_env(k)%helium%solute_interaction)
IF (helium_env(k)%helium%solute_interaction .EQ. helium_solute_intpot_nnp) THEN
IF (k == 1) THEN
NULLIFY (nnp_section)
nnp_section => section_vals_get_subs_vals(helium_section, "NNP")
CALL section_vals_get(nnp_section, explicit=explicit)
msg_str = "NNP section not explicitly stated. Using default file names."
IF (.NOT. explicit) CPWARN(msg_str)
END IF
ALLOCATE (helium_env(k)%helium%nnp)
CALL cp_logger_create(tmplogger, para_env=helium_env(k)%comm, template_logger=logger)
CALL cp_add_default_logger(tmplogger)
CALL helium_init_nnp(helium_env(k)%helium, helium_env(k)%helium%nnp, nnp_section)
CALL cp_rm_default_logger()
CALL cp_logger_release(tmplogger)
END IF
IF (helium_env(k)%helium%solute_interaction .EQ. helium_solute_intpot_none) THEN
WRITE (msg_str, '(A,I0,A)') &
"Solute found but no helium-solute interaction selected "// &
@ -835,6 +861,7 @@ CONTAINS
DO k = 1, SIZE(helium_env)
IF (k .EQ. 1) THEN
CALL helium_env(k)%comm%free()
DEALLOCATE (helium_env(k)%comm)
DEALLOCATE (helium_env(k)%env_all)
END IF
NULLIFY (helium_env(k)%env_all)
@ -990,6 +1017,17 @@ CONTAINS
NULLIFY (helium_env(k)%helium%ename)
END IF
! NNP interaction
IF (ASSOCIATED(helium_env(k)%helium%nnp)) THEN
CALL nnp_env_release(helium_env(k)%helium%nnp)
DEALLOCATE (helium_env(k)%helium%nnp)
NULLIFY (helium_env(k)%helium%nnp)
END IF
IF (ASSOCIATED(helium_env(k)%helium%nnp_sr_cut)) THEN
DEALLOCATE (helium_env(k)%helium%nnp_sr_cut)
NULLIFY (helium_env(k)%helium%nnp_sr_cut)
END IF
DEALLOCATE (helium_env(k)%helium)
END DO

334
src/motion/helium_nnp.F Normal file
View file

@ -0,0 +1,334 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2023 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Methods dealing with Neural Network interaction potential
!> \author Laura Duran
!> \date 2023-02-17
! **************************************************************************************************
MODULE helium_nnp
USE bibliography, ONLY: Behler2007,&
Behler2011,&
Schran2020a,&
Schran2020b,&
cite_reference
USE cell_methods, ONLY: cell_create
USE cell_types, ONLY: cell_release,&
cell_type,&
pbc
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_output_handling, ONLY: cp_p_file,&
cp_print_key_finished_output,&
cp_print_key_should_output,&
cp_print_key_unit_nr
USE cp_units, ONLY: cp_unit_from_cp2k
USE helium_types, ONLY: helium_solvent_type
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: default_path_length,&
default_string_length,&
dp
USE nnp_environment, ONLY: nnp_init_model
USE nnp_environment_types, ONLY: nnp_env_get,&
nnp_env_set,&
nnp_type
USE periodic_table, ONLY: get_ptable_info
USE physcon, ONLY: angstrom
#include "../base/base_uses.f90"
IMPLICIT NONE
PRIVATE
LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_nnp'
PUBLIC :: helium_init_nnp, &
helium_nnp_print
CONTAINS
! ***************************************************************************
!> \brief Read and initialize all the information for neural network potentials
!> \param helium ...
!> \param nnp ...
!> \param input ...
!> \date 2023-02-21
!> \author lduran
! **************************************************************************************************
SUBROUTINE helium_init_nnp(helium, nnp, input)
TYPE(helium_solvent_type), INTENT(INOUT) :: helium
TYPE(nnp_type), POINTER :: nnp
TYPE(section_vals_type), POINTER :: input
CHARACTER(len=default_path_length) :: msg_str
CHARACTER(len=default_string_length) :: elem
INTEGER :: i, ig, is, j
INTEGER, DIMENSION(3) :: periodicity
LOGICAL :: found
TYPE(cell_type), POINTER :: he_cell
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: sr_cutoff_section
CALL cite_reference(Behler2007)
CALL cite_reference(Behler2011)
CALL cite_reference(Schran2020a)
CALL cite_reference(Schran2020b)
NULLIFY (logger)
logger => cp_get_default_logger()
CALL nnp_env_set(nnp_env=nnp, nnp_input=input)
nnp%num_atoms = helium%solute_atoms + 1
CALL nnp_init_model(nnp, "HELIUM NNP")
periodicity = 0
IF (helium%periodic) periodicity = 1
NULLIFY (he_cell)
CALL cell_create(he_cell, hmat=helium%cell_m, &
periodic=periodicity, tag="HELIUM NNP")
CALL nnp_env_set(nnp, cell=he_cell)
CALL cell_release(he_cell)
! Set up arrays for calculation:
ALLOCATE (nnp%ele_ind(nnp%num_atoms))
ALLOCATE (nnp%nuc_atoms(nnp%num_atoms))
ALLOCATE (nnp%coord(3, nnp%num_atoms))
ALLOCATE (nnp%nnp_forces(3, nnp%num_atoms))
ALLOCATE (nnp%atoms(nnp%num_atoms))
ALLOCATE (nnp%sort(nnp%num_atoms - 1))
!fill arrays, assume that order will not change during simulation
ig = 1
is = 1
DO i = 1, nnp%n_ele
IF (nnp%ele(i) == 'He') THEN
nnp%atoms(ig) = 'He'
CALL get_ptable_info(nnp%atoms(ig), number=nnp%nuc_atoms(ig))
nnp%ele_ind(ig) = i
ig = ig + 1
END IF
DO j = 1, helium%solute_atoms
IF (nnp%ele(i) == helium%solute_element(j)) THEN
nnp%atoms(ig) = nnp%ele(i)
CALL get_ptable_info(nnp%atoms(ig), number=nnp%nuc_atoms(ig))
nnp%ele_ind(ig) = i
nnp%sort(is) = j
ig = ig + 1
is = is + 1
END IF
END DO
END DO
ALLOCATE (helium%nnp_sr_cut(nnp%n_ele))
helium%nnp_sr_cut = 0.0_dp
sr_cutoff_section => section_vals_get_subs_vals(nnp%nnp_input, "SR_CUTOFF")
CALL section_vals_get(sr_cutoff_section, n_repetition=is)
DO i = 1, is
CALL section_vals_val_get(sr_cutoff_section, "ELEMENT", c_val=elem, i_rep_section=i)
found = .FALSE.
DO ig = 1, nnp%n_ele
IF (TRIM(nnp%ele(ig)) == TRIM(elem)) THEN
found = .TRUE.
CALL section_vals_val_get(sr_cutoff_section, "RADIUS", r_val=helium%nnp_sr_cut(ig), &
i_rep_section=i)
END IF
END DO
IF (.NOT. found) THEN
msg_str = "SR_CUTOFF for element "//TRIM(elem)//" defined but not found in NNP"
CPWARN(msg_str)
END IF
END DO
helium%nnp_sr_cut(:) = helium%nnp_sr_cut(:)**2
RETURN
END SUBROUTINE helium_init_nnp
! **************************************************************************************************
!> \brief Print properties according to the requests in input file
!> \param nnp ...
!> \param print_section ...
!> \param ind_he ...
!> \date 2023-07-31
!> \author Laura Duran
! **************************************************************************************************
SUBROUTINE helium_nnp_print(nnp, print_section, ind_he)
TYPE(nnp_type), INTENT(INOUT) :: nnp
TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
INTEGER, INTENT(IN) :: ind_he
INTEGER :: unit_nr
LOGICAL :: file_is_new
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: print_key
NULLIFY (logger, print_key)
logger => cp_get_default_logger()
print_key => section_vals_get_subs_vals(print_section, "ENERGIES")
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
middle_name="helium-nnp-energies", is_new_file=file_is_new)
IF (unit_nr > 0) CALL helium_nnp_print_energies(nnp, unit_nr, file_is_new)
CALL cp_print_key_finished_output(unit_nr, logger, print_key)
END IF
print_key => section_vals_get_subs_vals(print_section, "FORCES_SIGMA")
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
middle_name="helium-nnp-forces-std", is_new_file=file_is_new)
IF (unit_nr > 0) CALL helium_nnp_print_force_sigma(nnp, unit_nr, file_is_new)
CALL cp_print_key_finished_output(unit_nr, logger, print_key)
END IF
CALL logger%para_env%sum(nnp%output_expol)
IF (nnp%output_expol) THEN
print_key => section_vals_get_subs_vals(print_section, "EXTRAPOLATION")
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
middle_name="-NNP-He-extrapolation")
IF (unit_nr > 0) CALL helium_nnp_print_expol(nnp, unit_nr, ind_he)
CALL cp_print_key_finished_output(unit_nr, logger, print_key)
END IF
END IF
END SUBROUTINE helium_nnp_print
! **************************************************************************************************
!> \brief Print NNP energies and standard deviation sigma
!> \param nnp ...
!> \param unit_nr ...
!> \param file_is_new ...
!> \date 2023-07-31
!> \author Laura Duran
! **************************************************************************************************
SUBROUTINE helium_nnp_print_energies(nnp, unit_nr, file_is_new)
TYPE(nnp_type), INTENT(INOUT) :: nnp
INTEGER, INTENT(IN) :: unit_nr
LOGICAL, INTENT(IN) :: file_is_new
CHARACTER(len=default_path_length) :: fmt_string
INTEGER :: i
REAL(KIND=dp) :: std
IF (file_is_new) THEN
WRITE (unit_nr, "(A1,1X,A20)", ADVANCE='no') "#", "NNP Average [a.u.],"
WRITE (unit_nr, "(A20)", ADVANCE='no') "NNP sigma [a.u.]"
DO i = 1, nnp%n_committee
WRITE (unit_nr, "(A17,I3)", ADVANCE='no') "NNP", i
END DO
WRITE (unit_nr, "(A)") ""
END IF
fmt_string = "(2X,2(F20.9))"
WRITE (fmt_string, "(A,I3,A)") "(2X", nnp%n_committee + 2, "(F20.9))"
std = SUM((SUM(nnp%atomic_energy, 1) - nnp%nnp_potential_energy)**2)
std = std/REAL(nnp%n_committee, dp)
std = SQRT(std)
WRITE (unit_nr, fmt_string) nnp%nnp_potential_energy, std, SUM(nnp%atomic_energy, 1)
END SUBROUTINE helium_nnp_print_energies
! **************************************************************************************************
!> \brief Print standard deviation sigma of NNP forces
!> \param nnp ...
!> \param unit_nr ...
!> \param file_is_new ...
!> \date 2023-07-31
!> \author Laura Duran
! **************************************************************************************************
SUBROUTINE helium_nnp_print_force_sigma(nnp, unit_nr, file_is_new)
TYPE(nnp_type), INTENT(INOUT) :: nnp
INTEGER, INTENT(IN) :: unit_nr
LOGICAL, INTENT(IN) :: file_is_new
INTEGER :: i, ig, j
REAL(KIND=dp), DIMENSION(3) :: var
IF (unit_nr > 0) THEN
IF (file_is_new) THEN
WRITE (unit_nr, "(A,1X,A)") "# NNP sigma of forces [a.u.] x, y, z coordinates"
END IF
ig = 1
DO i = 1, nnp%num_atoms
IF (nnp%ele(i) == 'He') THEN
var = 0.0_dp
DO j = 1, nnp%n_committee
var = var + (nnp%committee_forces(:, i, j) - nnp%nnp_forces(:, i))**2
END DO
WRITE (unit_nr, "(A4,1X,3E20.10)") nnp%ele(i), var
END IF
ig = ig + 1
END DO
END IF
END SUBROUTINE helium_nnp_print_force_sigma
! **************************************************************************************************
!> \brief Print structures with extrapolation warning
!> \param nnp ...
!> \param unit_nr ...
!> \param ind_he ...
!> \date 2023-10-11
!> \author Harald Forbert (harald.forbert@rub.de)
! **************************************************************************************************
SUBROUTINE helium_nnp_print_expol(nnp, unit_nr, ind_he)
TYPE(nnp_type), INTENT(INOUT) :: nnp
INTEGER, INTENT(IN) :: unit_nr, ind_he
CHARACTER(len=default_path_length) :: fmt_string
INTEGER :: i
REAL(KIND=dp) :: mass, unit_conv
REAL(KIND=dp), DIMENSION(3) :: com
TYPE(cell_type), POINTER :: cell
NULLIFY (cell)
CALL nnp_env_get(nnp_env=nnp, cell=cell)
nnp%expol = nnp%expol + 1
WRITE (unit_nr, *) nnp%num_atoms
WRITE (unit_nr, "(A,1X,I6)") "HELIUM-NNP extrapolation point N =", nnp%expol
! move to COM of solute and wrap the box
! coord not needed afterwards, therefore manipulation ok
com = 0.0_dp
mass = 0.0_dp
DO i = 1, nnp%num_atoms
IF (i == ind_he) CYCLE
CALL get_ptable_info(nnp%atoms(i), amass=unit_conv)
com(:) = com(:) + nnp%coord(:, i)*unit_conv
mass = mass + unit_conv
END DO
com(:) = com(:)/mass
DO i = 1, nnp%num_atoms
nnp%coord(:, i) = nnp%coord(:, i) - com(:)
nnp%coord(:, i) = pbc(nnp%coord(:, i), cell)
END DO
unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM("angstrom"))
fmt_string = "(A4,1X,3F20.10)"
DO i = 1, nnp%num_atoms
WRITE (unit_nr, fmt_string) &
nnp%atoms(i), &
nnp%coord(1, i)*unit_conv, &
nnp%coord(2, i)*unit_conv, &
nnp%coord(3, i)*unit_conv
END DO
END SUBROUTINE helium_nnp_print_expol
END MODULE helium_nnp

20
src/motion/helium_sampling.F
View file

@ -13,8 +13,12 @@
MODULE helium_sampling
USE cp_external_control, ONLY: external_control
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_log_handling, ONLY: cp_add_default_logger,&
cp_get_default_logger,&
cp_logger_create,&
cp_logger_release,&
cp_logger_type,&
cp_rm_default_logger
USE cp_output_handling, ONLY: cp_add_iter_level,&
cp_iterate,&
cp_rm_iter_level
@ -188,7 +192,12 @@ CONTAINS
TYPE(cp_logger_type), POINTER :: logger
NULLIFY (logger)
logger => cp_get_default_logger()
IF (SIZE(helium_env) < 1) THEN
logger => cp_get_default_logger()
ELSE
CALL cp_logger_create(logger, para_env=helium_env(1)%comm, template_logger=cp_get_default_logger())
CALL cp_add_default_logger(logger)
END IF
DO k = 1, SIZE(helium_env)
@ -390,6 +399,11 @@ CONTAINS
END IF
END IF
IF (SIZE(helium_env) > 0) THEN
CALL cp_rm_default_logger()
CALL cp_logger_release(logger)
END IF
RETURN
END SUBROUTINE helium_sample

8
src/motion/helium_types.F
View file

@ -21,7 +21,8 @@ MODULE helium_types
USE kinds, ONLY: default_string_length,&
dp,&
int_8
USE message_passing, ONLY: mp_comm_type
USE message_passing, ONLY: mp_para_env_type
USE nnp_environment_types, ONLY: nnp_type
USE parallel_rng_types, ONLY: rng_stream_type
USE splines_types, ONLY: spline_data_type
#include "../base/base_uses.f90"
@ -274,6 +275,9 @@ MODULE helium_types
LOGICAL :: interaction_pot_scan = .FALSE.!< whether to perform solute-helium interaction scan
TYPE(nnp_type), POINTER :: nnp => NULL() !< neural network potential
REAL(KIND=dp), DIMENSION(:), POINTER :: nnp_sr_cut => NULL() !< hard core cutoff in addition to the nnp
! temporary arrays for optimization
!
INTEGER, DIMENSION(:), POINTER :: itmp_atoms_1d => NULL()!< DIM(atoms) - same as permutation
@ -295,7 +299,7 @@ MODULE helium_types
! ***************************************************************************
TYPE helium_solvent_p_type
TYPE(helium_solvent_type), POINTER :: helium => NULL()
TYPE(mp_comm_type) :: comm = mp_comm_type()
TYPE(mp_para_env_type), POINTER :: comm => NULL()
INTEGER, DIMENSION(:), POINTER :: env_all => NULL()
END TYPE helium_solvent_p_type

10
src/motion/pint_methods.F
View file

@ -623,10 +623,12 @@ CONTAINS
pint_env%beadwise_constraints = .FALSE.
CALL section_vals_val_get(constraint_section, "PIMD_BEADWISE_CONSTRAINT", &
l_val=pint_env%beadwise_constraints)
IF (pint_env%beadwise_constraints) THEN
CALL pint_write_line("Using beadwise constraints")
ELSE
CALL pint_write_line("Using centroid constraints")
IF (pint_env%simpar%constraint) THEN
IF (pint_env%beadwise_constraints) THEN
CALL pint_write_line("Using beadwise constraints")
ELSE
CALL pint_write_line("Using centroid constraints")
END IF
END IF
IF (explicit) THEN

12
src/nnp_acsf.F
View file

@ -1159,7 +1159,7 @@ CONTAINS
rad = 0
DO s = 1, nnp%n_rad(i)
IF (nnp%rad(i)%ele(s) == nnp%rad(i)%symfgrp(j)%ele(1)) THEN
IF (ABS(nnp%rad(i)%funccut(s) - nnp%rad(i)%symfgrp(j)%cutoff) < 1.0e-7_dp) THEN
IF (ABS(nnp%rad(i)%funccut(s) - nnp%rad(i)%symfgrp(j)%cutoff) <= 1.0e-5_dp) THEN
rad = rad + 1
nnp%rad(i)%symfgrp(j)%symf(rad) = s
END IF
@ -1174,7 +1174,7 @@ CONTAINS
nnp%ang(i)%ele2(s) == nnp%ang(i)%symfgrp(j)%ele(2)) .OR. &
(nnp%ang(i)%ele1(s) == nnp%ang(i)%symfgrp(j)%ele(2) .AND. &
nnp%ang(i)%ele2(s) == nnp%ang(i)%symfgrp(j)%ele(1))) THEN
IF (ABS(nnp%ang(i)%funccut(s) - nnp%ang(i)%symfgrp(j)%cutoff) < 1.0e-7_dp) THEN
IF (ABS(nnp%ang(i)%funccut(s) - nnp%ang(i)%symfgrp(j)%cutoff) <= 1.0e-5_dp) THEN
ang = ang + 1
nnp%ang(i)%symfgrp(j)%symf(ang) = s
END IF
@ -1189,12 +1189,14 @@ CONTAINS
!> \brief Write symmetry function information
!> \param nnp ...
!> \param para_env ...
!> \param printtag ...
!> \date 2020-10-10
!> \author Christoph Schran (christoph.schran@rub.de)
! **************************************************************************************************
SUBROUTINE nnp_write_acsf(nnp, para_env)
SUBROUTINE nnp_write_acsf(nnp, para_env, printtag)
TYPE(nnp_type), INTENT(INOUT) :: nnp
TYPE(mp_para_env_type), POINTER :: para_env
CHARACTER(LEN=*), INTENT(IN) :: printtag
CHARACTER(len=default_string_length) :: my_label
INTEGER :: i, j, unit_nr
@ -1203,7 +1205,7 @@ CONTAINS
NULLIFY (logger)
logger => cp_get_default_logger()
my_label = "NNP| "
my_label = TRIM(printtag)//"| "
IF (para_env%is_source()) THEN
unit_nr = cp_logger_get_default_unit_nr(logger)
WRITE (unit_nr, '(1X,A,1X,10(I2,1X))') TRIM(my_label)//" Activation functions:", nnp%actfnct(:)
@ -1225,6 +1227,8 @@ CONTAINS
END DO
END IF
RETURN
END SUBROUTINE nnp_write_acsf
! **************************************************************************************************

68
src/nnp_environment.F
View file

@ -72,6 +72,7 @@ MODULE nnp_environment
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_environment'
PUBLIC :: nnp_init
PUBLIC :: nnp_init_model
CONTAINS
@ -229,7 +230,7 @@ CONTAINS
CALL distribution_1d_release(local_particles)
CALL distribution_1d_release(local_molecules)
CALL nnp_init_model(nnp_env=nnp_env)
CALL nnp_init_model(nnp_env=nnp_env, printtag="NNP")
CALL timestop(handle)
@ -238,11 +239,13 @@ CONTAINS
! **************************************************************************************************
!> \brief Initialize the Neural Network Potential
!> \param nnp_env ...
!> \param printtag ...
!> \date 2020-10-10
!> \author Christoph Schran (christoph.schran@rub.de)
! **************************************************************************************************
SUBROUTINE nnp_init_model(nnp_env)
SUBROUTINE nnp_init_model(nnp_env, printtag)
TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp_env
CHARACTER(LEN=*), INTENT(IN) :: printtag
CHARACTER(len=*), PARAMETER :: routineN = 'nnp_init_model'
INTEGER, PARAMETER :: def_str_len = 256, &
@ -274,7 +277,7 @@ CONTAINS
IF (logger%para_env%is_source()) THEN
unit_nr = cp_logger_get_default_unit_nr(logger)
WRITE (unit_nr, *) ""
WRITE (unit_nr, *) "NNP| Neural Network Potential Force Environment"
WRITE (unit_nr, *) TRIM(printtag)//"| Neural Network Potential Force Environment"
END IF
model_section => section_vals_get_subs_vals(nnp_env%nnp_input, "MODEL")
@ -300,7 +303,7 @@ CONTAINS
IF (logger%para_env%is_source()) THEN
unit_nr = cp_logger_get_default_unit_nr(logger)
WRITE (unit_nr, *) "NNP| Reading NNP input from file: ", TRIM(file_name)
WRITE (unit_nr, *) TRIM(printtag)//"| Reading NNP input from file: ", TRIM(file_name)
END IF
CALL parser_search_string(parser, "number_of_elements", .TRUE., found, line, &
@ -308,7 +311,8 @@ CONTAINS
IF (found) THEN
READ (line, *) dummy, nnp_env%n_ele
ELSE
CPABORT("NNP| number of elements missing in NNP_INPUT_FILE")
CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| number of elements missing in NNP_INPUT_FILE")
END IF
CALL parser_search_string(parser, "scale_symmetry_functions_sigma", .TRUE., found, &
@ -353,7 +357,8 @@ CONTAINS
IF (found) THEN
READ (line, *) dummy, nnp_env%cut_type
ELSE
CPABORT("NNP| no cutoff type specified in NNP_INPUT_FILE")
CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| no cutoff type specified in NNP_INPUT_FILE")
END IF
CALL parser_search_string(parser, "global_hidden_layers_short", .TRUE., found, line, &
@ -361,7 +366,8 @@ CONTAINS
IF (found) THEN
READ (line, *) dummy, nnp_env%n_hlayer
ELSE
CPABORT("NNP| number of hidden layers missing in NNP_INPUT_FILE")
CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| number of hidden layers missing in NNP_INPUT_FILE")
END IF
nnp_env%n_layer = nnp_env%n_hlayer + 2
@ -395,7 +401,8 @@ CONTAINS
EXIT
END IF
ELSE
CPABORT("NNP| elements not specified in NNP_INPUT_FILE")
CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| elements not specified in NNP_INPUT_FILE")
END IF
END DO
@ -417,7 +424,8 @@ CONTAINS
END DO
IF (i == nele) EXIT
ELSE
CPABORT("NNP| atom energies are not specified")
CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| atom energies are not specified")
END IF
END DO
END IF
@ -427,7 +435,8 @@ CONTAINS
IF (found) THEN
READ (line, *) dummy, nnp_env%n_hnodes(:)
ELSE
CPABORT("NNP| global_nodes_short not specified in NNP_INPUT_FILE")
CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"NNP| global_nodes_short not specified in NNP_INPUT_FILE")
END IF
CALL parser_search_string(parser, "global_activation_short", .TRUE., found, line, &
@ -435,7 +444,8 @@ CONTAINS
IF (found) THEN
READ (line, *) dummy, cactfnct(:)
ELSE
CPABORT("NNP| global_activation_short not specified in NNP_INPUT_FILE")
CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| global_activation_short not specified in NNP_INPUT_FILE")
END IF
DO i = 1, nnp_env%n_hlayer + 1
@ -459,7 +469,8 @@ CONTAINS
CASE ("h")
nnp_env%actfnct(i) = nnp_actfnct_quad
CASE DEFAULT
CPABORT("NNP| Activation function unkown")
CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| Activation function unkown")
END SELECT
END DO
@ -483,13 +494,15 @@ CONTAINS
ELSE IF (symfnct_type .EQ. 3) THEN
nnp_env%n_ang(i) = nnp_env%n_ang(i) + 1
ELSE
CPABORT("NNP| Symmetry function type not supported")
CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| Symmetry function type not supported")
END IF
END IF
END DO
first = .FALSE.
ELSE
IF (first) CPABORT("NNP| no symfunction_short specified in NNP_INPUT_FILE")
IF (first) CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| no symfunction_short specified in NNP_INPUT_FILE")
! no additional symfnct found
EXIT
END IF
@ -586,13 +599,15 @@ CONTAINS
nnp_env%max_cut = nnp_env%ang(i)%funccut(nnp_env%n_ang(i))
END IF
ELSE
CPABORT("NNP| Symmetry function type not supported")
CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| Symmetry function type not supported")
END IF
END IF
END DO
first = .FALSE.
ELSE
IF (first) CPABORT("NNP| no symfunction_short specified in NNP_INPUT_FILE")
IF (first) CALL cp_abort(__LOCATION__, TRIM(printtag)// &
"| no symfunction_short specified in NNP_INPUT_FILE")
! no additional symfnct found
EXIT
END IF
@ -621,13 +636,13 @@ CONTAINS
! sort symmetry functions and output information
CALL nnp_sort_acsf(nnp_env)
CALL nnp_write_acsf(nnp_env, logger%para_env)
CALL nnp_write_arc(nnp_env, logger%para_env)
CALL nnp_write_acsf(nnp_env, logger%para_env, printtag)
CALL nnp_write_arc(nnp_env, logger%para_env, printtag)
! 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
WRITE (unit_nr, *) "NNP| Reading scaling information from file: ", TRIM(file_name)
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", &
c_val=file_name)
@ -696,12 +711,12 @@ CONTAINS
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
WRITE (unit_nr, *) "NNP| Initializing weights for model: ", i_com
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
WRITE (unit_nr, *) "NNP| Reading weights from file: ", TRIM(file_name)
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)
n_weight = 0
@ -744,8 +759,11 @@ CONTAINS
nnp_env%expol = 0
! Bias the standard deviation of committee disagreement
bias_section => section_vals_get_subs_vals(nnp_env%nnp_input, "BIAS")
CALL section_vals_get(bias_section, explicit=explicit)
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
@ -768,7 +786,7 @@ CONTAINS
END IF
nnp_env%bias_e_avrg(:) = work
IF (logger%para_env%is_source()) THEN
WRITE (unit_nr, *) "NNP| Biasing is aligned by shifting the energy prediction of the C-NNP members"
WRITE (unit_nr, *) TRIM(printtag)//"| Biasing is aligned by shifting the energy prediction of the C-NNP members"
END IF
END IF
ELSE
@ -777,7 +795,7 @@ CONTAINS
END IF
IF (logger%para_env%is_source()) THEN
WRITE (unit_nr, *) "NNP| NNP force environment initialized"
WRITE (unit_nr, *) TRIM(printtag)//"| NNP force environment initialized"
END IF
CALL timestop(handle)

64
src/nnp_environment_types.F
View file

@ -58,8 +58,7 @@ MODULE nnp_environment_types
! Public subroutines ***
PUBLIC :: nnp_env_release, &
nnp_env_set, &
nnp_env_get, &
nnp_env_create
nnp_env_get
INTEGER, PARAMETER, PUBLIC :: &
nnp_cut_cos = 1, &
@ -82,8 +81,8 @@ MODULE nnp_environment_types
!> \date 2020-10-10
! **************************************************************************************************
TYPE nnp_type
TYPE(nnp_acsf_rad_type), DIMENSION(:), POINTER :: rad ! DIM(n_ele)
TYPE(nnp_acsf_ang_type), DIMENSION(:), POINTER :: ang ! DIM(n_ele)
TYPE(nnp_acsf_rad_type), DIMENSION(:), POINTER :: rad => NULL() ! DIM(n_ele)
TYPE(nnp_acsf_ang_type), DIMENSION(:), POINTER :: ang => NULL() ! DIM(n_ele)
INTEGER, DIMENSION(:), ALLOCATABLE :: n_rad ! # radial symfnct for this element
INTEGER, DIMENSION(:), ALLOCATABLE :: n_ang ! # angular symfnct for this element
INTEGER :: n_ele ! # elements
@ -102,7 +101,7 @@ MODULE nnp_environment_types
REAL(KIND=dp) :: scmin !scale
REAL(KIND=dp) :: max_cut !largest cutoff
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: atom_energies !DIM(n_ele)
TYPE(nnp_arc_type), POINTER, DIMENSION(:) :: arc ! DIM(n_ele)
TYPE(nnp_arc_type), POINTER, DIMENSION(:) :: arc => NULL() ! DIM(n_ele)
INTEGER :: n_committee
INTEGER :: n_hlayer
INTEGER :: n_layer
@ -122,11 +121,12 @@ MODULE nnp_environment_types
ALLOCATABLE, DIMENSION(:) :: atoms
REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: nnp_forces
REAL(KIND=dp) :: nnp_potential_energy
TYPE(cp_subsys_type), POINTER :: subsys
TYPE(section_vals_type), POINTER :: nnp_input, &
force_env_input
TYPE(cell_type), POINTER :: cell_ref
LOGICAL :: use_ref_cell
TYPE(cp_subsys_type), POINTER :: subsys => NULL()
TYPE(section_vals_type), POINTER :: nnp_input => NULL()
TYPE(section_vals_type), POINTER :: force_env_input => NULL()
TYPE(cell_type), POINTER :: cell => NULL()
TYPE(cell_type), POINTER :: cell_ref => NULL()
LOGICAL :: use_ref_cell = .FALSE.
! bias
LOGICAL :: bias
LOGICAL :: bias_align
@ -265,28 +265,6 @@ MODULE nnp_environment_types
CONTAINS
! **************************************************************************************************
!> \brief Create data structure that holds all the information for neural
!> network potentials
!> \param nnp_env ...
!> \date 2020-10-10
!> \author Christoph Schran (christoph.schran@rub.de)
! **************************************************************************************************
SUBROUTINE nnp_env_create(nnp_env)
TYPE(nnp_type), INTENT(OUT) :: nnp_env
NULLIFY (nnp_env%rad)
NULLIFY (nnp_env%ang)
NULLIFY (nnp_env%arc)
NULLIFY (nnp_env%subsys)
NULLIFY (nnp_env%nnp_input)
NULLIFY (nnp_env%force_env_input)
NULLIFY (nnp_env%cell_ref)
nnp_env%use_ref_cell = .FALSE.
END SUBROUTINE nnp_env_create
! **************************************************************************************************
!> \brief Release data structure that holds all the information for neural
!> network potentials
@ -398,6 +376,9 @@ CONTAINS
IF (ASSOCIATED(nnp_env%subsys)) THEN
CALL cp_subsys_release(nnp_env%subsys)
END IF
IF (ASSOCIATED(nnp_env%cell)) THEN
CALL cell_release(nnp_env%cell)
END IF
IF (ASSOCIATED(nnp_env%cell_ref)) THEN
CALL cell_release(nnp_env%cell_ref)
END IF
@ -466,6 +447,10 @@ CONTAINS
END IF
IF (PRESENT(nnp_forces)) nnp_forces = nnp_env%nnp_forces
! note cell will be overwritten if subsys is associated
! helium_env uses nnp without subsys
IF (PRESENT(cell)) cell => nnp_env%cell
IF (PRESENT(subsys)) subsys => nnp_env%subsys
IF (ASSOCIATED(nnp_env%subsys)) THEN
CALL cp_subsys_get(nnp_env%subsys, &
@ -503,6 +488,7 @@ CONTAINS
!> \param local_molecules All molecules on this particular node
!> \param nnp_input ...
!> \param force_env_input Pointer to the force_env input section
!> \param cell ...
!> \param cell_ref The reference simulation cell
!> \param use_ref_cell Logical which indicates if reference
!> simulation cell is used
@ -515,7 +501,7 @@ CONTAINS
SUBROUTINE nnp_env_set(nnp_env, nnp_forces, subsys, &
atomic_kind_set, particle_set, local_particles, &
molecule_kind_set, molecule_set, local_molecules, &
nnp_input, force_env_input, cell_ref, &
nnp_input, force_env_input, cell, cell_ref, &
use_ref_cell, nnp_potential_energy)
TYPE(nnp_type), INTENT(INOUT) :: nnp_env
@ -532,7 +518,7 @@ CONTAINS
POINTER :: molecule_set
TYPE(distribution_1d_type), OPTIONAL, POINTER :: local_molecules
TYPE(section_vals_type), OPTIONAL, POINTER :: nnp_input, force_env_input
TYPE(cell_type), OPTIONAL, POINTER :: cell_ref
TYPE(cell_type), OPTIONAL, POINTER :: cell, cell_ref
LOGICAL, INTENT(IN), OPTIONAL :: use_ref_cell
REAL(KIND=dp), INTENT(IN), OPTIONAL :: nnp_potential_energy
@ -554,6 +540,16 @@ CONTAINS
END IF
nnp_env%subsys => subsys
END IF
IF (PRESENT(cell)) THEN
IF (ASSOCIATED(cell)) THEN
CALL cell_retain(cell)
CALL cell_release(nnp_env%cell)
nnp_env%cell => cell
END IF
IF (ASSOCIATED(nnp_env%subsys)) THEN
CALL cp_subsys_set(nnp_env%subsys, cell=cell)
END IF
END IF
IF (PRESENT(atomic_kind_set)) THEN
CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
CALL cp_subsys_set(nnp_env%subsys, atomic_kinds=atomic_kinds)

8
src/nnp_model.F
View file

@ -41,10 +41,12 @@ CONTAINS
!> \brief Write neural network architecture information
!> \param nnp ...
!> \param para_env ...
!> \param printtag ...
! **************************************************************************************************
SUBROUTINE nnp_write_arc(nnp, para_env)
SUBROUTINE nnp_write_arc(nnp, para_env, printtag)
TYPE(nnp_type), INTENT(IN) :: nnp
TYPE(mp_para_env_type), INTENT(IN) :: para_env
CHARACTER(LEN=*), INTENT(IN) :: printtag
CHARACTER(len=default_string_length) :: my_label
INTEGER :: i, j, unit_nr
@ -53,7 +55,7 @@ CONTAINS
NULLIFY (logger)
logger => cp_get_default_logger()
my_label = "NNP| "
my_label = TRIM(printtag)//"| "
IF (para_env%is_source()) THEN
unit_nr = cp_logger_get_default_unit_nr(logger)
DO i = 1, nnp%n_ele
@ -66,6 +68,8 @@ CONTAINS
END DO
END IF
RETURN
END SUBROUTINE nnp_write_arc
! **************************************************************************************************

112
src/start/input_cp2k_motion.F
View file

@ -36,7 +36,7 @@ MODULE input_cp2k_motion
helium_cell_shape_octahedron, helium_forces_average, helium_forces_last, &
helium_mdist_exponential, helium_mdist_gaussian, helium_mdist_linear, &
helium_mdist_quadratic, helium_mdist_singlev, helium_mdist_uniform, &
helium_sampling_ceperley, helium_sampling_worm, helium_solute_intpot_mwater, &
helium_sampling_ceperley, helium_sampling_worm, helium_solute_intpot_mwater, helium_solute_intpot_nnp, &
helium_solute_intpot_none, integrate_exact, integrate_numeric, ls_2pnt, ls_3pnt, ls_fit, &
ls_gold, ls_none, matrix_init_cholesky, matrix_init_diagonal, numerical, perm_cycle, &
perm_plain, propagator_pimd, propagator_rpmd, propagator_cmd, transformation_normal, transformation_stage
@ -61,7 +61,8 @@ MODULE input_cp2k_motion
section_type
USE input_val_types, ONLY: integer_t, &
logical_t, &
real_t
real_t, &
char_t
USE kinds, ONLY: dp
USE string_utilities, ONLY: s2a
#include "../base/base_uses.f90"
@ -2108,7 +2109,7 @@ CONTAINS
CALL section_create(section, __LOCATION__, name="HELIUM", &
description="The section that controls optional helium solvent"// &
" environment (highly experimental, not for general use yet)", &
n_keywords=31, n_subsections=9, repeats=.FALSE.)
n_keywords=31, n_subsections=11, repeats=.FALSE.)
NULLIFY (keyword)
CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
@ -2156,15 +2157,17 @@ CONTAINS
CALL keyword_create(keyword, __LOCATION__, name="SOLUTE_INTERACTION", &
description="Interaction potential between helium and the solute", &
usage="SOLUTE_INTERACTION (NONE | MWATER)", &
usage="SOLUTE_INTERACTION (NONE | MWATER | NNP)", &
default_i_val=helium_solute_intpot_none, &
enum_c_vals=s2a("NONE", "MWATER"), &
enum_c_vals=s2a("NONE", "MWATER", "NNP"), &
enum_i_vals=(/ &
helium_solute_intpot_none, &
helium_solute_intpot_mwater/), &
helium_solute_intpot_mwater, &
helium_solute_intpot_nnp/), &
enum_desc=s2a( &
"No interaction with solute", &
"Test interaction with wrong Water"))
"Test interaction with wrong Water", &
"Interaction with NNP"))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
@ -2229,6 +2232,101 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
! Helium-solute interaction NNP
NULLIFY (subsection)
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.)
CALL keyword_create(keyword, __LOCATION__, name="NNP_INPUT_FILE_NAME", &
description="File containing the input information for the setup "// &
"of the NNP (n2p2/RuNNer format). ", &
repeats=.FALSE., default_lc_val="input.nn")
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="SCALE_FILE_NAME", &
description="File containing the scaling information for the symmetry "// &
"functions of the NNP. ", &
repeats=.FALSE., default_lc_val="scaling.data")
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, "// &
"if the distance between solvent and specified solute element becomes "// &
"smaller than the given cutoff, an artifical repulsive potential is "// &
"introduced. Note this is only meant to prevent such configurations, "// &
"not to physically sample them.", &
n_keywords=2, n_subsections=0, repeats=.TRUE.)
CALL keyword_create(keyword, __LOCATION__, name="ELEMENT", &
description="Solute element for which the short range cutoff is in effect", &
repeats=.FALSE., default_c_val="none")
CALL section_add_keyword(subsubsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="RADIUS", &
description="Short range cutoff in Angstrom, below this cutoff, the energy "// &
"is replaced by a sizable positive value plus a 1/r**2 term to guide particles "// &
"away from each other.", &
default_r_val=cp_unit_to_cp2k(0.0_dp, "angstrom"), &
repeats=.FALSE., type_of_var=real_t, unit_str="angstrom")
CALL section_add_keyword(subsubsection, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(subsection, subsubsection)
CALL section_release(subsubsection)
NULLIFY (subsubsection)
CALL section_create(subsubsection, __LOCATION__, name="MODEL", &
description="Section for a single NNP model. If this section is repeated, "// &
"a committee model (C-NNP)is used where the NNP members share the same "// &
"symmetry functions. ", &
n_keywords=1, n_subsections=0, repeats=.TRUE.)
CALL keyword_create(keyword, __LOCATION__, name="WEIGHTS", &
description="File containing the weights for the artificial neural "// &
"networks of the NNP. The specified name is extended by .XXX.data ", &
repeats=.FALSE., default_lc_val="weights")
CALL section_add_keyword(subsubsection, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(subsection, subsubsection)
CALL section_release(subsubsection)
! Create the PRINT subsection
NULLIFY (subsubsection)
CALL section_create(subsubsection, __LOCATION__, name="PRINT", &
description="Section of possible print options in NNP code.", &
n_keywords=0, n_subsections=3, repeats=.FALSE.)
NULLIFY (print_key, keyword)
CALL cp_print_key_section_create(print_key, __LOCATION__, "ENERGIES", &
description="Controls the printing of the NNP energies.", &
print_level=medium_print_level, common_iter_levels=1)
CALL section_add_subsection(subsubsection, print_key)
CALL section_release(print_key)
CALL cp_print_key_section_create(print_key, __LOCATION__, "FORCES_SIGMA", &
description="Controls the printing of the STD per atom of the NNP forces.", &
print_level=medium_print_level, common_iter_levels=1)
CALL section_add_subsection(subsubsection, print_key)
CALL section_release(print_key)
CALL cp_print_key_section_create(print_key, __LOCATION__, "EXTRAPOLATION", &
description="If activated, output structures with extrapolation "// &
"warning in xyz-format", &
print_level=medium_print_level, common_iter_levels=1)
CALL section_add_subsection(subsubsection, print_key)
CALL section_release(print_key)
CALL section_add_subsection(subsection, subsubsection)
CALL section_release(subsubsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection) ! release NNP subsection
! Ceperley's sampling algorithm
NULLIFY (subsection)
CALL section_create(subsection, __LOCATION__, name="CEPERLEY", &
description="Enables sampling with Ceperley's algorithm", &

1
tests/Pimd/regtest-2/TEST_FILES
View file

@ -16,4 +16,5 @@ he32_only_worm.inp 40 1e-11 3.284
he32_only_worm_restart.inp 40 1e-11 3.1056893819001927E-004
water-in-helium-worm.inp 9 1e-11 -4.2435303435780906E-004
water-in-helium-striding.inp 9 1e-11 9.2149634888323367E-002
water_in_helium_nnp.inp 9 1e-11 0.188542848154504
#EOF

88
tests/Pimd/regtest-2/input.nn Normal file
View file

@ -0,0 +1,88 @@
number_of_elements 3
elements O H He
global_hidden_layers_short 2
global_nodes_short 25 25
global_activation_short t t l
cutoff_type 1
# radial H H
symfunction_short H 2 H 0.00000 0.0 6.000 ! eta rshift funccutoff
symfunction_short H 2 H 0.01284 0.0 6.000 ! eta rshift funccutoff
symfunction_short H 2 H 0.04711 0.0 6.000 ! eta rshift funccutoff
symfunction_short H 2 H 0.11421 0.0 6.000 ! eta rshift funccutoff
symfunction_short H 2 H 0.38554 0.0 6.000 ! eta rshift funccutoff
symfunction_short H 2 H 1.14214 0.0 6.000 ! eta rshift funccutoff
# radial H O / O H
symfunction_short H 2 O 0.00000 0.0 4.000 ! eta rshift funccutoff
symfunction_short H 2 O 0.01284 0.0 4.000 ! eta rshift funccutoff
symfunction_short H 2 O 0.04711 0.0 4.000 ! eta rshift funccutoff
symfunction_short H 2 O 0.11421 0.0 4.000 ! eta rshift funccutoff
symfunction_short H 2 O 0.38554 0.0 4.000 ! eta rshift funccutoff
symfunction_short H 2 O 1.14214 0.0 4.000 ! eta rshift funccutoff
symfunction_short O 2 H 0.00000 0.0 4.000 ! eta rshift funccutoff
symfunction_short O 2 H 0.01284 0.0 4.000 ! eta rshift funccutoff
symfunction_short O 2 H 0.04711 0.0 4.000 ! eta rshift funccutoff
symfunction_short O 2 H 0.11421 0.0 4.000 ! eta rshift funccutoff
symfunction_short O 2 H 0.38554 0.0 4.000 ! eta rshift funccutoff
symfunction_short O 2 H 1.14214 0.0 4.000 ! eta rshift funccutoff
#
# radial H He / He H
symfunction_short H 2 He 0.00000 0.0 22.000 ! eta rshift funccutoff
symfunction_short H 2 He 0.01284 0.0 22.000 ! eta rshift funccutoff
symfunction_short H 2 He 0.04711 0.0 22.000 ! eta rshift funccutoff
symfunction_short H 2 He 0.11421 0.0 22.000 ! eta rshift funccutoff
symfunction_short H 2 He 0.38554 0.0 22.000 ! eta rshift funccutoff
#
symfunction_short He 2 H 0.00000 0.0 22.000 ! eta rshift funccutoff
symfunction_short He 2 H 0.01284 0.0 22.000 ! eta rshift funccutoff
symfunction_short He 2 H 0.04711 0.0 22.000 ! eta rshift funccutoff
symfunction_short He 2 H 0.11421 0.0 22.000 ! eta rshift funccutoff
symfunction_short He 2 H 0.38554 0.0 22.000 ! eta rshift funccutoff
#
# radial O He / He O
symfunction_short O 2 He 0.00000 0.0 22.000 ! eta rshift funccutoff
symfunction_short O 2 He 0.01284 0.0 22.000 ! eta rshift funccutoff
symfunction_short O 2 He 0.04711 0.0 22.000 ! eta rshift funccutoff
symfunction_short O 2 He 0.11421 0.0 22.000 ! eta rshift funccutoff
symfunction_short O 2 He 0.38554 0.0 22.000 ! eta rshift funccutoff
#
symfunction_short He 2 O 0.00000 0.0 22.000 ! eta rshift funccutoff
symfunction_short He 2 O 0.01284 0.0 22.000 ! eta rshift funccutoff
symfunction_short He 2 O 0.04711 0.0 22.000 ! eta rshift funccutoff
symfunction_short He 2 O 0.11421 0.0 22.000 ! eta rshift funccutoff
symfunction_short He 2 O 0.38554 0.0 22.000 ! eta rshift funccutoff
#
# angular
symfunction_short H 3 H He 0.0 1.0 2.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short H 3 H He 0.0 1.0 16.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short H 3 H He 0.0 -1.0 2.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short H 3 H He 0.0 -1.0 16.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short O 3 H H 0.0 1.0 2.0 6.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short O 3 H H 0.0 -1.0 2.0 6.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short O 3 H H 0.0 -1.0 16.0 6.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short O 3 He H 0.0 1.0 2.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short O 3 He H 0.0 1.0 16.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short O 3 He H 0.0 -1.0 2.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short O 3 He H 0.0 -1.0 16.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short H 3 O H 0.0 1.0 2.0 6.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short H 3 O H 0.0 1.0 16.0 6.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short H 3 O H 0.0 -1.0 2.0 6.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short He 3 O H 0.0 1.0 2.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short He 3 O H 0.0 1.0 16.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short He 3 O H 0.0 -1.0 2.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short He 3 H H 0.0 1.0 2.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short He 3 H H 0.0 1.0 16.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
symfunction_short He 3 H H 0.0 -1.0 2.0 22.00000 ! central_atom type neighbor_atom1 neighbor_atom2 eta lambda zeta funccutoff
scale_symmetry_functions
center_symmetry_functions
#normalize_nodes

59
tests/Pimd/regtest-2/scaling.data Normal file
View file

@ -0,0 +1,59 @@
1 1 0.056817790 0.750692450 0.535422570
1 2 0.049035400 0.733801881 0.511770326
1 3 0.033094917 0.690558470 0.453923969
1 4 0.015326365 0.613129380 0.359844582
1 5 0.000681635 0.379056130 0.145144367
1 6 0.000000116 0.099158979 0.014140815
1 7 0.303742678 1.374101777 0.950667957
1 8 0.235810207 1.287115579 0.845399854
1 9 0.120225731 1.081696249 0.620860609
1 10 0.032412417 0.771739614 0.345340187
1 11 0.000162353 0.213805507 0.039411695
1 12 0.000000000 0.011145065 0.000284075
1 13 0.000000000 0.971823088 0.674142370
1 14 0.000000000 0.904632499 0.448357163
1 15 0.000000000 0.747181452 0.200857852
1 16 0.000000000 0.513834679 0.067521608
1 17 0.000000000 0.113060914 0.004983123
1 18 0.000173036 0.089741938 0.014016891
1 19 0.159049241 1.528033231 0.924189394
1 20 0.021928251 0.741271301 0.264708443
1 21 0.000000000 2.966883486 0.538884027
1 22 0.000000000 2.907583844 0.942197930
1 23 0.000000000 2.010270751 0.098215870
1 24 0.000000000 1.984445110 0.269188161
2 1 0.000000000 2.840136700 2.022427109
2 2 0.000000000 0.925422032 0.681974012
2 3 0.000000000 2.479841745 1.345071490
2 4 0.000000000 0.763221458 0.458412591
2 5 0.000000000 1.737259739 0.602573557
2 6 0.000000000 0.456334883 0.193352815
2 7 0.000000000 0.979460492 0.202564825
2 8 0.000000000 0.166699990 0.045058226
2 9 0.000000000 0.179451002 0.014949368
2 10 0.000000000 0.002840834 0.000295471
2 11 0.000000000 0.269412607 0.019135491
2 12 0.000000000 0.039142871 0.003080234
2 13 0.000000000 4.553843569 2.800523813
2 14 0.000000000 4.827097878 3.161841792
2 15 0.000000000 3.450986149 1.434508332
2 16 0.000000000 4.185789633 2.325371247
3 1 0.904224152 2.066881809 1.606267709
3 2 0.833469552 2.007545008 1.535310978
3 3 0.670584553 1.857396706 1.361771907
3 4 0.438119010 1.595201594 1.079533747
3 5 0.078461817 0.862496097 0.435433102
3 6 0.000656033 0.155828033 0.042422445
3 7 0.000000000 0.925422032 0.681974012
3 8 0.000000000 0.763221458 0.458412591
3 9 0.000000000 0.456334883 0.193352815
3 10 0.000000000 0.166699990 0.045058226
3 11 0.000000000 0.002840834 0.000295471
3 12 0.226700653 1.450755585 0.737233133
3 13 0.026515912 0.630416258 0.242989934
3 14 0.000024620 0.427149376 0.024691249
3 15 0.000000000 2.773800494 0.993854746
3 16 0.000000000 3.196696257 1.300572740
3 17 0.000000000 1.553136834 0.176102454
3 18 0.000000000 1.767843771 0.364970085
-0.0022195411 0.0037903600

104
tests/Pimd/regtest-2/water_in_helium_nnp.inp Normal file
View file

@ -0,0 +1,104 @@
&GLOBAL
PROJECT_NAME pwater-in-helium
RUN_TYPE PINT
PRINT_LEVEL low
&END GLOBAL
&MOTION
&PINT
P 16
PROC_PER_REPLICA 1
NUM_STEPS 5
DT 0.5
NRESPA 1
PROPAGATOR RPMD
HARM_INT EXACT
TEMP 7.5
T_TOL 0.0
TRANSFORMATION NORMAL
&PILE
&END PILE
&HELIUM
SOLUTE_INTERACTION NNP
POTENTIAL_FILE_NAME lj-test-pot.dat
COORD_INIT_TEMP -1
GET_FORCES LAST
NATOMS 20
NBEADS 16
NUM_ENV 1
INOROT 80
IROT 1000
PERIODIC F
CELL_SHAPE OCTAHEDRON
SAMPLING_METHOD WORM
DROPLET_RADIUS 20
&WORM
CENTROID_DRMAX 0.2
STAGING_L 5
ALLOW_OPEN T
OPEN_CLOSE_SCALE 100.0
SHOW_STATISTICS T
CENTROID_WEIGHT 5
STAGING_WEIGHT 10
OPEN_CLOSE_WEIGHT 20
HEAD_TAIL_WEIGHT 10
CRAWL_WEIGHT 10
CRAWL_REPETITION 2
SWAP_WEIGHT 100
&END WORM
&NNP
NNP_INPUT_FILE_NAME ./input.nn
SCALE_FILE_NAME ./scaling.data
&MODEL
WEIGHTS ./weights
&END MODEL
&SR_CUTOFF
ELEMENT H
RADIUS 1.25
&END SR_CUTOFF
&SR_CUTOFF
ELEMENT O
RADIUS 2.05
&END SR_CUTOFF
! &PRINT
! &EXTRAPOLATION
! &END EXTRAPOLATION
! &END PRINT
&END NNP
&END HELIUM
&END PINT
&END MOTION
&FORCE_EVAL
METHOD FIST
&MM
&POISSON
&EWALD
GMAX 25
&END EWALD
&END POISSON
&FORCEFIELD
PARM_FILE_NAME ../../Fist/sample_pot/water.pot
PARMTYPE CHM
&CHARGE
ATOM OT
CHARGE -0.8476
&END CHARGE
&CHARGE
ATOM HT
CHARGE 0.4238
&END CHARGE
&END FORCEFIELD
&END MM
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
PERIODIC NONE
&END CELL
&TOPOLOGY
COORD_FILE_NAME ../../Fist/sample_pdb/water_1.pdb
COORD_FILE_FORMAT PDB
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

1301
tests/Pimd/regtest-2/weights.001.data Normal file
View file

File diff suppressed because it is too large Load diff

1101
tests/Pimd/regtest-2/weights.002.data Normal file
View file

File diff suppressed because it is too large Load diff

1151
tests/Pimd/regtest-2/weights.008.data Normal file
View file

File diff suppressed because it is too large Load diff