diff --git a/src/motion/input_cp2k_md.F b/src/motion/input_cp2k_md.F
index daae38b8a7..c4396f20be 100644
--- a/src/motion/input_cp2k_md.F
+++ b/src/motion/input_cp2k_md.F
@@ -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)
diff --git a/src/motion/integrator.F b/src/motion/integrator.F
index b0f72df75a..96de427a30 100644
--- a/src/motion/integrator.F
+++ b/src/motion/integrator.F
@@ -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
diff --git a/src/motion/reftraj_types.F b/src/motion/reftraj_types.F
index d9b5f736c5..4f9b9ce4a5 100644
--- a/src/motion/reftraj_types.F
+++ b/src/motion/reftraj_types.F
@@ -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