If the lipids were subject to image centering during the simulation, then you must account for that. Otherwise, the r^2 vs. t data will be contaminated by large jumps across the unit cell.
In order to have the box size taken from the stored trajectory frames, you must both set up CRYSTAL and use XBOX etc.; the lattice must be orthogonal (all angles 90 deg) which means only CUBIC, TETRAGONAL, or ORTHORHOMBIC lattice types. (The same restriction applies to MERGE COOR UNFOLD; only orthogonal lattice types are supported.)
Calculating diffusion accurately for just one molecule is highly problematic, regardless of the approach; standard procedure is to average over molecules and/or replicate simulations. Lipids pose bigger challenges due to the presence of short time diffusion based on "rattling in a cage", and long time diffusion based a jump model with lipids changing places. Several replicate simulations might be needed to get statistically meaningful results for a small number of lipids.
While I haven't had time to clean up them up for the Script Archive, I do have CHARMM inputs, csh batch scripts, and a Fortran90 program which do these steps:
(1) image unfold the coordinate trajectory to a scratch volume (MERGE)
(2) for each molecule, compute the r^2 vs. t difference correlation function (CORREL), in both 3D and 2D
(3) compute D from input r^2 vs t data, in 3D or 2D, with optional finite size correction (cubic only); provides standard errors based on grouping (F90 prog)
Available via email request.