Reftraj: add option to wrap coordinates before calculation (#5479)

This commit is contained in:
HE Zilong 2026-06-29 04:20:33 +08:00 committed by GitHub
parent ed4c696b7e
commit 31876b56ae
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
3 changed files with 54 additions and 3 deletions

View file

@ -48,7 +48,10 @@ MODULE input_cp2k_md
USE kinds, ONLY: dp
USE reftraj_types, ONLY: REFTRAJ_EVAL_ENERGY,&
REFTRAJ_EVAL_ENERGY_FORCES,&
REFTRAJ_EVAL_NONE
REFTRAJ_EVAL_NONE,&
REFTRAJ_WRAP_CENTRAL,&
REFTRAJ_WRAP_NONE,&
REFTRAJ_WRAP_POSITIVE
USE string_utilities, ONLY: newline,&
s2a
#include "../base/base_uses.f90"
@ -606,6 +609,22 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="WRAP", &
description="Whether and how atoms would be wrapped to inside the "// &
"simulation box according to periodic boundary condition before any "// &
"subsequent analysis. Useful when there is need to convert the raw, "// &
"unwrapped coordinates produced from `MOTION%PRINT%TRAJECTORY` for "// &
"post-processing. This keyword will be ignored when requesting MSD "// &
"calculation because wrapping coordinates makes displacement discontinuous.", &
default_i_val=REFTRAJ_WRAP_NONE, &
enum_i_vals=[REFTRAJ_WRAP_NONE, REFTRAJ_WRAP_POSITIVE, REFTRAJ_WRAP_CENTRAL], &
enum_c_vals=s2a("NONE", "POSITIVE", "CENTRAL"), &
enum_desc=s2a("Do not wrap atoms at all and use coordinates exactly as input", &
"Wrap atoms so that fractional coordinates are within [0, 1]", &
"Wrap atoms so that fractional coordinates are within [-1/2, +1/2]"))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL create_msd_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)

View file

@ -32,7 +32,8 @@ MODULE integrator
USE cell_methods, ONLY: init_cell,&
read_xyz_comment
USE cell_types, ONLY: cell_type,&
parse_cell_line
parse_cell_line,&
pbc
USE constraint, ONLY: rattle_control,&
shake_control,&
shake_roll_control,&
@ -88,6 +89,9 @@ MODULE integrator
USE qs_environment_types, ONLY: get_qs_env
USE reftraj_types, ONLY: REFTRAJ_EVAL_ENERGY_FORCES,&
REFTRAJ_EVAL_NONE,&
REFTRAJ_WRAP_CENTRAL,&
REFTRAJ_WRAP_NONE,&
REFTRAJ_WRAP_POSITIVE,&
reftraj_type
USE reftraj_util, ONLY: compute_msd_reftraj
USE rt_propagation_methods, ONLY: propagation_step
@ -1688,6 +1692,25 @@ CONTAINS
CALL init_cell(cell)
END IF
! Wrap coordinates if requested
SELECT CASE (reftraj_env%info%wrap)
CASE (REFTRAJ_WRAP_NONE)
! Do Nothing
CASE (REFTRAJ_WRAP_POSITIVE)
! Wrap to positive range
DO i = 1, nparticle
particle_set(i)%r(1:3) = pbc(particle_set(i)%r(1:3), cell, positive_range=.TRUE.)
END DO
CASE (REFTRAJ_WRAP_CENTRAL)
! Wrap to halfway, i.e. origin is at the center
DO i = 1, nparticle
particle_set(i)%r(1:3) = pbc(particle_set(i)%r(1:3), cell)
END DO
CASE DEFAULT
! Should not reach here
CPABORT("Option invalid or unavailable for reftraj_env%info%wrap")
END SELECT
![ADAPT] update input structure with new coordinates, make new labels
CALL qmmmx_update_force_env(force_env, force_env%root_section)
! no pointers into force_env%subsys to update

View file

@ -33,9 +33,12 @@ MODULE reftraj_types
INTEGER, PARAMETER, PUBLIC :: REFTRAJ_EVAL_NONE = 101
INTEGER, PARAMETER, PUBLIC :: REFTRAJ_EVAL_ENERGY = 102
INTEGER, PARAMETER, PUBLIC :: REFTRAJ_EVAL_ENERGY_FORCES = 103
INTEGER, PARAMETER, PUBLIC :: REFTRAJ_WRAP_NONE = 201
INTEGER, PARAMETER, PUBLIC :: REFTRAJ_WRAP_POSITIVE = 202
INTEGER, PARAMETER, PUBLIC :: REFTRAJ_WRAP_CENTRAL = 203
! **************************************************************************************************
!> \brief parameters related to the analysis of previously generated trajecorties
!> \brief parameters related to the analysis of previously generated trajectories
!> \author MI
! **************************************************************************************************
TYPE reftraj_info_type
@ -43,6 +46,7 @@ MODULE reftraj_types
INTEGER :: last_snapshot = 0
INTEGER :: stride = 0
INTEGER :: eval = REFTRAJ_EVAL_NONE
INTEGER :: wrap = REFTRAJ_WRAP_NONE
LOGICAL :: variable_volume = .FALSE.
LOGICAL :: msd = .FALSE.
TYPE(cp_parser_type), POINTER :: traj_parser => NULL()
@ -117,9 +121,14 @@ CONTAINS
CALL section_vals_val_get(reftraj_section, "LAST_SNAPSHOT", i_val=reftraj%info%last_snapshot)
CALL section_vals_val_get(reftraj_section, "STRIDE", i_val=reftraj%info%stride)
CALL section_vals_val_get(reftraj_section, "EVAL", i_val=reftraj%info%eval)
CALL section_vals_val_get(reftraj_section, "WRAP", i_val=reftraj%info%wrap)
CALL section_vals_val_get(reftraj_section, "MSD%_SECTION_PARAMETERS_", &
l_val=reftraj%info%msd)
IF (reftraj%info%msd .AND. (reftraj%info%wrap /= REFTRAJ_WRAP_NONE)) THEN
CPWARN("MSD calculation is incompatible with wrapping coordinates")
reftraj%info%wrap = REFTRAJ_WRAP_NONE
END IF
END SUBROUTINE create_reftraj