program n_particles_LJ_PBC parameter(n=3) real x(n),y(n),x_old(n),y_old(n),vx(n),vy(n) real force_x(n),force_y(n),force_x_old(n),force_y_old(n) integer pgopen do i=1,n write(*,*)'What is initial vx of particle',i,'of',n,'?' read(*,*) vx(i) write(*,*)'What is initial vy of particle',i,'of',n,'?' read(*,*) vy(i) end do x_range = 10.0 y_range = 10.0 C In this code, n must be >= 2. do i=1,n x(i) = 5.0*(i-1)/(n-1) y(i) = 5.0*(i-1)/(n-1) end do if( pgopen('/xwin') <= 0 ) stop call pgenv(-0.15,x_range+0.15,-0.15,y_range+0.15,1,-1) call pgsfs(2) dt = 0.0001 dt2 = dt*dt do i_step =1, 10000000 c find forces from pairs do i=1,n force_x(i)=0.0 force_y(i)=0.0 end do do i=1,n-1 do j=i+1,n c call lj_force(x(i),y(i),x(j),y(j),force_x_tmp,force_y_tmp) call PBC_lj_force(x(i),y(i),x(j),y(j),force_x_tmp, & force_y_tmp,x_range,y_range) force_x(i)=force_x(i)+force_x_tmp force_y(i)=force_y(i)+force_y_tmp force_x(j)=force_x(j)-force_x_tmp force_y(j)=force_y(j)-force_y_tmp end do end do c do i=1,n c call wall_reflect(x(i),y(i),vx(i),vy(i),x_range,y_range) call PBC(x(i),y(i),x_range,y_range) x(i) = x(i) + vx(i) * dt + 0.5 * force_x(i) * dt2 y(i) = y(i) + vy(i) * dt + 0.5 * force_y(i) * dt2 vx(i) = vx(i) + (force_x(i)+force_x_old(i))/2 * dt vy(i) = vy(i) + (force_y(i)+force_y_old(i))/2 * dt call pgbbuf call pgsci(0) call pgcirc(x_old(i),y_old(i),0.1) call pgsci(1) call pgcirc(x(i),y(i),0.1) call pgebuf x_old(i) = x(i) y_old(i) = y(i) force_x_old(i) = force_x(i) force_y_old(i) = force_y(i) enddo enddo call pgclos end subroutine PBC_distance(x1,y1,x2,y2,dis,x2mx1,y2my1,x_range, &y_range) dis_min=sqrt( (0.5*x_range)**2 + (0.5*y_range)**2 ) do i=-1,1 do j=-1,1 x2mx1_temp=i*x_range+x2-x1 y2my1_temp=j*y_range+y2-y1 r=sqrt( (x2mx1_temp)**2 + (y2my1_temp)**2 ) if(r < dis_min)then dis_min = r x2mx1=x2mx1_temp y2my1=y2my1_temp dis=r endif end do end do end subroutine subroutine PBC_lj_force(x1,y1,x2,y2,force_x,force_y,x_range, &y_range) sigma = 1.0 epsilon = 1.0 c r = sqrt( (x2-x1)**2 + (y2-y1)**2 ) call PBC_distance(x1,y1,x2,y2,dis,x2mx1,y2my1,x_range,y_range) r=dis s_o_r = sigma / r force = (-1)*24*epsilon*(2*s_o_r**12-s_o_r**6)/r c force_x = force*(x2-x1)/r c force_y = force*(y2-y1)/r force_x = force* x2mx1 /r force_y = force* y2my1 /r end subroutine subroutine lj_force(x1,y1,x2,y2,force_x,force_y) sigma = 1.0 epsilon = 1.0 r = sqrt( (x2-x1)**2 + (y2-y1)**2 ) s_o_r = sigma / r force = (-1)*24*epsilon*(2*s_o_r**12-s_o_r**6)/r force_x = force*(x2-x1)/r force_y = force*(y2-y1)/r end subroutine subroutine wall_reflect(x,y,vx,vy,x_range,y_range) if(x>x_range) vx = (-1)*vx if(x<0.0) vx = (-1)*vx if(y>y_range) vy = (-1)*vy if(y<0.0) vy = (-1)*vy end subroutine subroutine PBC(x,y,x_range,y_range) if(x>x_range) x = x - x_range if(x<0.0) x = x + x_range if(y>y_range) y = y - y_range if(y<0.0) y = y + y_range end subroutine