program mol_dyn integer i,j,k,t,tmax,p,n parameter(p=5,n=p**3,tmax=3000) real r(1:3,0:2,1:n),rij(1:3,1:n,1:n),rij2(1:n,1:n) real Fij(1:3,1:n,1:n),Fi(1:3,1:n),h2 real x(1:3,1:n),vi(1:3,1:n),v0 parameter(h2=0.001,v0=0.005) forall(k=1:3,i=1:n) r(k,0,i)=1.1*modulo((i-1)/p**(k-1),p) call random_number(x); vi=v0*x forall(k=1:3) vi(k,:)=vi(k,:)-sum(vi(k,:))/n r(:,1,:)=r(:,0,:)+sqrt(h2)*vi(:,:) Fij=0 do t=1,tmax forall(i=1:n,j=1:n,i splot for [i=7:131] 'fort.'.i with lines notitle