Try this:
compute a property/atom xu yu id # relevant info
#endfinder
compute molchunk all chunk/atom molecule ids once
compute molmaxid all reduce/chunk molchunk max c_a[3]
compute s_maxid all chunk/spread/atom molchunk c_molmaxid
variable is_end atom c_s_maxid==id
#parity chunking id: 1 2 3 4 5 6 ...
variable parity atom id%2 # 1 0 1 0 1 0 ...
variable oddpair atom id+v_parity # 2 2 4 4 6 6 ...
variable evnpair atom id-v_parity # 0 2 2 4 4 6 ...
compute oddchunk all chunk/atom v_oddpair ids once
compute evnchunk all chunk/atom v_evnpair ids once
# process unwrapped displacements
compute oddsum all reduce/chunk oddchunk sum c_a[1] c_a[2]
compute evnsum all reduce/chunk evnchunk sum c_a[1] c_a[2]
compute s_oddsum all chunk/spread/atom oddsum[*]
compute s_evnsum all chunk/spread/atom evnsum[*]
variable dx atom (v_parity*s_oddsum[1]+(1-v_parity)*s_evnsum[1])-xu
variable dy atom (v_parity*s_oddsum[2]+(1-v_parity)*s_evnsum[2])-yu
variable dr atom sqrt(v_dx*v_dx+v_dy*v_dy+1e-6*v_is_end) # avoid divBy0