I have just started using dl_poly classic to work in Molecular Dynamics simulations. It produces a HISTORY file which records the trajectory after applying boundary conditions. Now, here I don't understand, how I am suppossed to calculate diffusivity using folded trajectories. Secondly is there a way in dl_poly to unfold the trajectories may be using some counters that they might already written (i.e. to do it without altering the source code).
However, from one of the answers based on a question here, I wrote a code to unfold the folded trajectories (am trying to verify this from a program I wrote myself for the simulation of argon liquid, not dl_poly) -
open(unit=11,file='test.out',status='unknown')
rx(:,1) = x(:,1)
do j=2,tsteps
do i=1,nm
tmp1 = x(i,j-1) - x(i,j) + boxl/2.d0
tmp1 = tmp1 - boxL*floor(tmp1/boxL)
tmp2 = boxL/2.d0 - tmp1
rx(i,j) = rx(i,j-1) + tmp2
write(11,*) i,j,x(i,j),ox(i,j),rx(i,j),ox(i,j)-rx(i,j)
end do
end do
In the above program, x(i,j)
are folded coordinates, ox(i,j)
are the unfolded ones that are there already (this is the output a simulation program I wrote) and rx(i,j)
are the coordinates I am unfolding from x(i,j)
, where i
corressponds to $i^{th}$ particle and j
corresponds to $j^{th}$ timestep.
The problem I am facing with the above program is that, in some cases, I get a difference of +boxl
or -boxl
from the actual unfolded coordinates (i.e. between, ox(i,j)
and rx(i,j)
)
This is the answer I used to write down this logic, but I might written it down wrongly as well.