	program lj
c
c	molecular dynamics of Lennard-Jonesium
c
	implicit real*8 (a-h,o-z)
	parameter( npd=4000 )
	common/pos/x0(npd),y0(npd),z0(npd)
	common/vel/x1(npd),y1(npd),z1(npd)
	common/der/x2(npd),y2(npd),z2(npd),x3(npd),y3(npd),z3(npd),
     &  	x4(npd),y4(npd),z4(npd),x5(npd),y5(npd),z5(npd)
	common/for/fx(npd),fy(npd),fz(npd)
	common/dis/dax(npd),day(npd),daz(npd),x0l(npd),y0l(npd),z0l(npd)
	common/nablst/list(npd*50),nabors(npd)
	common/parm/f02,f12,f32,f42,f52
	common/propa/rdel,rmax,rdmax,rlist,fshft,eshft,eshfta,np1,ksort
	parameter( nbins=1000 )
	common/rdf/idist(nbins),voli(nbins),bin(nbins)
c
c       system size and related constants
c
        nc = 5
	np = 4*nc**3
	if(np.gt.npd) then
	  write(6,*) 'too many particles for npd - reset'
	  stop
	endif
	part = np
	np1  = np-1
c
c       operating parameters
c
	dr = 0.8
	tr = 1.2
	em = 1.0
c
c       run flags and parameters
c
	kb     = 0
	kbskp  = 500
	maxkb  = 2000
	ksort  = 10
	iwrite = 100
	kwrite = 50
	seed   = 0.3783626117
	third  = 1./3.
	pi     = 4.*datan(1.d0)
c
c       time-step and its multiples
c
	delta=0.005
	delsq=delta*delta
	deltsq=0.5*delsq/em
        fkb=iwrite*part
c
c       parameters in predictor-corrector method
c
	f02=3./16.
	f12=251./360.
	f32=11./18.
	f42=1./6.
	f52=1./60.
c
c       distances for potential cut-off, verlet list, etc.
c
	vol=part/dr
	cube=vol**third
	cube2=0.5*cube
	rc=2.5
	rlist=(rc+1.0)**2
	rdmax=cube2*cube2
	if(rlist.gt.rdmax) rlist=rdmax
	if(rc.gt.cube2) rc=cube2-0.35
	rmax=rc*rc
	vscale=third/tr
c
c       scale factors for velocities during equilibration
c
	aheat=delsq*part*3.*tr
	rtr=delta*dsqrt(3.*tr)
c
c       parameters for sampling g(r)
c
	rdel=0.025
	nr=cube2/rdel
	if(nr.gt.nbins) then
	  write(6,*) 'rdf bin array too small, nr =',nr
	  stop
	endif
	vi=2.*ksort/(dr*np*iwrite)
	do i=1,nr
	  rvol=4.*pi*rdel**3*(i**2-i+third)
	  voli(i)=vi/rvol
	  bin(i)=(i-0.5)*rdel
	enddo
c
c       shifted-force constants
c
	rrmax=1./dsqrt(rmax)
	rrmax6=rrmax**6
	eshft=rrmax6*(28.-52.*rrmax6)
	eshfta=48.*rrmax*rrmax6*(rrmax6-.5)
	fshft=eshfta
c
c       corrections for long-range interactions
c
      	rc3=rc**3
      	rc9=rc3**3
      	core=8.*pi*dr*(1./9./rc9-third/rc3)
      	corv=96.*pi*dr*(1./6./rc3-1./9./rc9)
c
c       initialize sum accumulators
c
	sumt=0.
	sume=0.
	sumv=0.
c
c       generate initial positions of atoms
c
	call fcc(cube,nc,np)
c
c       set initial velocities of atoms
c
	call intvel(aheat,rtr,part,np,seed)
c
c       assign initial accelerations based on initial positions
c
	call eval(totv,tote,cube,cube2,kb,np,nabtot)
c
c       scale accelerations and store starting positions
c
	do i=1,np
	  x2(i)=fx(i)*deltsq
	  y2(i)=fy(i)*deltsq
	  z2(i)=fz(i)*deltsq
	  x0l(i)=x0(i)
	  y0l(i)=y0(i)
	  z0l(i)=z0(i)
	enddo
c
c	output file headers
c
	write(10,*) 'Instantaneous p, T, U'
	write(10,*)
	write(11,*) 'Average p, T, U, <r^2>, E'
	write(11,*)
	write(20,*) 'Radial distribution function'
	write(20,*)
c
c       time-step loop 
c	--------------
c
	do 1000 kb=1,maxkb

	call predct(np)
	call eval(totv,tote,cube,cube2,kb,np,nabtot)
	call corr(np,deltsq,cube)
c
c       calculate mean square displacement & kinetic energy
c
	sumvel=0.
	tdist=0.
	do i=1,np
	  tdist=tdist+dax(i)**2+day(i)**2+daz(i)**2
	  sumvel=sumvel+x1(i)**2+y1(i)**2+z1(i)**2
	enddo
	tdist=tdist/part
	ek=sumvel/(2.*part*delsq)
c
c       instantaneous properties at kwrite intervals
c
        if(mod(kb,kwrite).eq.0) then
	  t1=sumvel/(3*delsq*part)
	  p1=t1*dr*(1.-vscale*(totv/part+corv))
	  e1=tote/part+core
	  write(10,'(i8,3e12.4)') kb,p1,t1,e1
	endif
c
c       accumulate sums for property averages
c
	sumt=sumt+sumvel
	sume=sume+tote
	sumv=sumv+totv
c
c       average property calculation at iwrite intervals
c
	if(mod(kb,iwrite).eq.0) then
	  tmp=sumt/(3.*delsq*fkb)
	  enr=sume/fkb+core
	  vir=vscale*(sumv/fkb+corv)
	  pres=tmp*dr*(1.-vir)
	  etot=enr+tmp*3./2.
	  write(11,'(i8,5e12.4)') kb,pres,tmp,enr,tdist,etot
	  write(20,*) kb
	  do i=1,nr	
	    write(20,'(f8.4,e16.7)') bin(i),idist(i)*voli(i)
	    idist(i)=0.
	  enddo
	  sumt=0.
	  sume=0.
	  sumv=0.
	endif
c
c       during start of run scale velocities for temperature
c
	if(kb.le.kbskp) call eqbrat(sumvel,aheat,np)
c
c	end time step loop
c
1000	continue

	stop
	end



	subroutine fcc(cube,nc,np)
c
c     	calculate positions of sites on face-centered cubic lattice
c
	implicit real*8 (a-h,o-z)
	parameter ( npd=4000 )
	common/pos/x0(npd),y0(npd),z0(npd)

	xl=cube/float(nc)
	y=.5*xl
	x0(1)=0.
	y0(1)=0.
	z0(1)=0.
	x0(2)=0.
	y0(2)=y
	z0(2)=y
	x0(3)=y
	y0(3)=0.
	z0(3)=y
	x0(4)=y
	y0(4)=y
	z0(4)=0.
	m=0
	do 10 i=1,nc
	  xli = xl*(i-1)
	do 10 j=1,nc
	  xlj = xl*(j-1)
	do 10 k=1,nc
	  xlk = xl*(k-1)
	  do 11 ij=1,4
	    x0(ij+m)=x0(ij)+xlk
	    y0(ij+m)=y0(ij)+xlj
	    z0(ij+m)=z0(ij)+xli
 11	  continue
	  m=m+4
 10     continue
	if (m.ne.np) then
		write(6,*) 'initialization error - wrong np'
		stop
	endif
	return
	end



	subroutine intvel(aheat,rtr,part,np,seed)
c
c       assign initial velocities to atoms
c
	implicit real*8 (a-h,o-z)
	parameter ( npd=4000 )
	common/vel/x1(npd),y1(npd),z1(npd)

	x=seed
	sumx=0.
	sumy=0.
	sumz=0.
	do i=1,np
	  call random(x,xx)
	  x=xx
	  xx=2.0*(xx-0.5)
	  call random(x,yy)
	  x=yy
	  yy=2.0*(yy-0.5)
	  call random(x,zz)
	  x=zz
	  zz=2.0*(zz-0.5)
	  xyz=1./dsqrt(xx*xx+yy*yy+zz*zz)
	  x1(i)=xx*xyz*rtr
	  y1(i)=yy*xyz*rtr
	  z1(i)=zz*xyz*rtr
	  sumx=sumx+x1(i)
	  sumy=sumy+y1(i)
	  sumz=sumz+z1(i)
	enddo
c
c       shift velocities so that total momentum=zero
c
	x=0
	do i=1,np
	  x1(i)=x1(i)-sumx/part
	  y1(i)=y1(i)-sumy/part
	  z1(i)=z1(i)-sumz/part
	  x=x+x1(i)**2+y1(i)**2+z1(i)**2
	enddo
c
c       scale velocities to desired temperature
c
	heat=dsqrt(aheat/x)
	do i=1,np
	  x1(i)=x1(i)*heat
	  y1(i)=y1(i)*heat
	  z1(i)=z1(i)*heat
	enddo

	return
	end



	subroutine predct(np)
c
c       use fifth-order taylor series to integrate
c
	implicit real*8 (a-h,o-z)
	parameter ( npd=4000 )
	common/pos/x0(npd),y0(npd),z0(npd)
	common/vel/x1(npd),y1(npd),z1(npd)
	common/der/x2(npd),y2(npd),z2(npd),x3(npd),y3(npd),z3(npd)
     &	     ,x4(npd),y4(npd),z4(npd),x5(npd),y5(npd),z5(npd)
	common/for/fx(npd),fy(npd),fz(npd)
 
	do i=1,np
	  x0(i)=x0(i)+x1(i)+x2(i)+x3(i)+x4(i)+x5(i)
	  y0(i)=y0(i)+y1(i)+y2(i)+y3(i)+y4(i)+y5(i)
	  z0(i)=z0(i)+z1(i)+z2(i)+z3(i)+z4(i)+z5(i)
	  x1(i)=x1(i)+2.*x2(i)+3.*x3(i)+4.*x4(i)+5.*x5(i)
	  y1(i)=y1(i)+2.*y2(i)+3.*y3(i)+4.*y4(i)+5.*y5(i)
	  z1(i)=z1(i)+2.*z2(i)+3.*z3(i)+4.*z4(i)+5.*z5(i)
	  x2(i)=x2(i)+3.*x3(i)+6.*x4(i)+10.*x5(i)
	  y2(i)=y2(i)+3.*y3(i)+6.*y4(i)+10.*y5(i)
	  z2(i)=z2(i)+3.*z3(i)+6.*z4(i)+10.*z5(i)
	  x3(i)=x3(i)+4.*x4(i)+10.*x5(i)
	  y3(i)=y3(i)+4.*y4(i)+10.*y5(i)
	  z3(i)=z3(i)+4.*z4(i)+10.*z5(i)
	  x4(i)=x4(i)+5.*x5(i)
	  y4(i)=y4(i)+5.*y5(i)
	  z4(i)=z4(i)+5.*z5(i)
c
c       initialize force arrays
c
	  fx(i)=0.
	  fy(i)=0.
	  fz(i)=0.
	enddo

	return
	end



	subroutine eval(totv,tote,cube,cube2,kb,np,nabtot)
c
c       evaluate forces on atoms using pairwise additive lennard-jones
c	  (6-12) intermolecular potential
c
	implicit real*8 (a-h,o-z)
	parameter ( npd=4000 )
	common/pos/x0(npd),y0(npd),z0(npd)
	common/for/fx(npd),fy(npd),fz(npd)
	common/nablst/list(npd*50),nabors(npd)
	common/propa/rdel,rmax,rdmax,rlist,fshft,eshft,eshfta,np1,ksort
	parameter( nbins=1000 )
	common/rdf/idist(nbins),voli(nbins),bin(nbins)
 
	k=0
	tote=0.
	totv=0.
c
c       set flag for list up-date
c
	lchk=2
	if(mod(kb,ksort).eq.0) lchk=1
c
c       outer loop over atoms
c
	do 500 i=1,np1
	if(lchk.eq.2) then
c
c       use neighbor-list
c
	  jbegin=nabors(i)
	  jend=nabors(i+1)-1
	  if(jbegin.gt.jend) go to 500
	else
c
c       loop over all atoms at ksort intervals
c
          nabors(i)=k+1
	  jbegin=i+1
	  jend=np
	endif
c
c       store position of atom i
c
        xi=x0(i)
	yi=y0(i)
	zi=z0(i)
c
c       inner loops over atoms
c
	do 495 jx= jbegin,jend
	  j=jx
	  if(lchk.eq.2)j=list(jx)
c
c       distance between atom i and atom j
c
	  x = xi - x0(j)
	  y = yi - y0(j)
	  z = zi - z0(j)
c
c       find image of atom j closest to atom i
c
	  xr=abs(x)
	  yr=abs(y)
	  zr=abs(z)
	  if(xr.gt.cube2) x=(xr-cube)*sign(1.d0,x)
	  if(yr.gt.cube2) y=(yr-cube)*sign(1.d0,y)
	  if(zr.gt.cube2) z=(zr-cube)*sign(1.d0,z)
c
c       calculate r(ij)
c
	  rsq=x*x+y*y+z*z
	  rctr=dsqrt(rsq)
c
c       increment counter for rdf & update list at intervals
c
	if(lchk.eq.1) then
	  if ( rsq .gt. rdmax ) go to 495
	  ij = rctr/rdel + 1.
	  idist(ij)=idist(ij)+1
	  if ( rsq .gt. rlist ) go to 495
	    k=k+1
	    list(k)=j
	endif
c
c       evaluate forces on atoms
c
        if(rsq.gt.rmax)goto 495
	  rsi=1./rsq
	  r6=rsi**3
	  rpl = 48.*r6*(r6-.5)-rctr*fshft
	  rp=rpl*rsi
	  repx=rp*x
	  fx(i)=fx(i)+repx
	  fx(j)=fx(j)-repx
	  repy=rp*y
	  fy(i)=fy(i)+repy
	  fy(j)=fy(j)-repy
	  repz=rp*z
	  fz(i)=fz(i)+repz
	  fz(j)=fz(j)-repz
c
c       accummulate energy and virial
c
	  tote = tote + 4.*r6*(r6-1)+eshft+rctr*eshfta
	  totv=totv-rpl
 495    continue
 500    continue
	if(lchk.ne.1) return
	nabors(np)=k+1
	nabtot=k

	return
	end



	subroutine corr(np,deltsq,cube)
c
c       correct predicted positions and their derivatives
c
	implicit real*8 (a-h,o-z)
	parameter ( npd=4000 )
	common/pos/x0(npd),y0(npd),z0(npd)
	common/vel/x1(npd),y1(npd),z1(npd)
	common/der/x2(npd),y2(npd),z2(npd),x3(npd),y3(npd),z3(npd)
     &	    ,x4(npd),y4(npd),z4(npd),x5(npd),y5(npd),z5(npd)
	common/for/fx(npd),fy(npd),fz(npd)
	common/dis/dax(npd),day(npd),daz(npd),x0l(npd),y0l(npd),z0l(npd)
	common/parm/f02,f12,f32,f42,f52
c
	do i=1,np
	    xerr=x2(i)-deltsq*fx(i)
	    yerr=y2(i)-deltsq*fy(i)
	    zerr=z2(i)-deltsq*fz(i)
	    x0(i)=x0(i)-xerr*f02
	    x1(i)=x1(i)-xerr*f12
	    x2(i)=x2(i)-xerr
	    x3(i)=x3(i)-xerr*f32
	    x4(i)=x4(i)-xerr*f42
	    x5(i)=x5(i)-xerr*f52
	    y0(i)=y0(i)-yerr*f02
	    y1(i)=y1(i)-yerr*f12
	    y2(i)=y2(i)-yerr
	    y3(i)=y3(i)-yerr*f32
	    y4(i)=y4(i)-yerr*f42
	    y5(i)=y5(i)-yerr*f52
	    z0(i)=z0(i)-zerr*f02
	    z1(i)=z1(i)-zerr*f12
	    z2(i)=z2(i)-zerr
	    z3(i)=z3(i)-zerr*f32
	    z4(i)=z4(i)-zerr*f42
	    z5(i)=z5(i)-zerr*f52
c
c     displacements
c
	    dax(i)=dax(i)-x0(i)+x0l(i)
	    day(i)=day(i)-y0(i)+y0l(i)
	    daz(i)=daz(i)-z0(i)+z0l(i)
c
c     check for conservation of particles in primary cell
c
	    xx=x0(i)
	    yy=y0(i)
	    zz=z0(i)
	    if(xx.le.0.)   xx=xx+cube
	    if(xx.ge.cube) xx=xx-cube
	    if(yy.le.0.)   yy=yy+cube
	    if(yy.ge.cube) yy=yy-cube
	    if(zz.le.0.)   zz=zz+cube
	    if(zz.ge.cube) zz=zz-cube
c
c     store new positions
c
	    x0(i)=xx
	    y0(i)=yy
	    z0(i)=zz
	    x0l(i)=xx
	    y0l(i)=yy
	    z0l(i)=zz
	enddo

	return
	end



	subroutine eqbrat(sumvel,aheat,np)
c
c       scale velocities during initial time-steps
c
	implicit real*8 (a-h,o-z)
	parameter ( npd=4000 )
	common/vel/x1(npd),y1(npd),z1(npd)
c
	heat=dsqrt(aheat/sumvel)
	do i=1,np
	  x1(i)=x1(i)*heat
	  y1(i)=y1(i)*heat
	  z1(i)=z1(i)*heat
	enddo

	return
	end



	subroutine random(x,xr)
c
c       subroutine for pseudo-random numbers
c
	implicit real*8 (a-h,o-z)
	data k,j,m,rm/5701,3612,566927,566927./
	ix=int(x*rm)
	kr=j*ix+k
	irand=mod(kr,m)
	xr=(float(irand)+0.5)/rm
	return
	end



	block data
	implicit real*8 (a-h,o-z)
	parameter ( npd=4000, np6 = npd*6, np9 = npd*9 )
	common/der/x2(npd),y2(npd),z2(npd),x3(npd),y3(npd),z3(npd)
     a	    ,x4(npd),y4(npd),z4(npd),x5(npd),y5(npd),z5(npd)
	common/for/fx(npd),fy(npd),fz(npd)
	common/dis/dax(npd),day(npd),daz(npd),x0l(npd),y0l(npd),z0l(npd)
	parameter( nbins=1000 )
	common/rdf/idist(nbins),voli(nbins),bin(nbins)
	data x3,y3,z3,x4,y4,z4,x5,y5,z5/np9*0.0/
	data dax,day,daz,fx,fy,fz/np6*0.0/
	data idist/1000*0/
	end

