	program lj_vv_2
c
c	lj_vv.f modified for a Langevin thermostat - changes in CAPS
c
c	molecular dynamics of Lennard-Jonesium
c
	parameter( npd=4000 )
	common/pos/x0(npd),y0(npd),z0(npd)
	common/vel/x1(npd),y1(npd),z1(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)
	COMMON/THERM/GAMMA,CONST
c
c       system size and related constants
c
        nc = 6
	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.0
	em = 1.0
c
c       run flags and parameters
c
	kb     = 0
	maxkb  = 10000
	kbskp  = 1000
	ksort  = 10
	iwrite = 500
	kwrite = 100
	seed   = 0.3783626117
	third  = 1./3.
	pi     = 4.*atan(1.0)
c
c       time-step and its multiples
c
	delta=0.005
	delt1=0.5*delta/em
	delt2=delt1*delta
	fkb=iwrite*part
C
C	THERMOSTAT PARAMETERS
C
	gamma = 0.125
	const = sqrt(6.*tr*gamma/delta)
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
	aheat=part*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./sqrt(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,part,np,seed)
c
c       assign initial accelerations based on initial positions
c
	call eval(totv,tote,cube,cube2,kb,np,nabtot)
c
c       store starting positions
c
	do i=1,np
	  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
c
c	first velocity-verlet step
c
	do i=1,np
		x0(i)=x0(i)+x1(i)*delta+fx(i)*delt2
		y0(i)=y0(i)+y1(i)*delta+fy(i)*delt2
		z0(i)=z0(i)+z1(i)*delta+fz(i)*delt2
		x1(i)=x1(i)+fx(i)*delt1
		y1(i)=y1(i)+fy(i)*delt1
		z1(i)=z1(i)+fz(i)*delt1
		fx(i)=0.
		fy(i)=0.
		fz(i)=0.
	enddo
c
c	force evaluation
c
	call eval(totv,tote,cube,cube2,kb,np,nabtot)
C
C	THERMOSTAT
C
	IF(KB.GT.KBSKP) CALL LANGEVIN(NP,SEED)
c
c       second velocity-verlet step
c
        do i=1,np
		x1(i)=x1(i)+fx(i)*delt1
		y1(i)=y1(i)+fy(i)*delt1
		z1(i)=z1(i)+fz(i)*delt1
        enddo
c
c	periodic BC and displacements
c
	call correct(np,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)
c
c       instantaneous properties at kwrite intervals
c
        if(mod(kb,kwrite).eq.0) then
	  t1=sumvel/(3*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.*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       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
	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,part,np,seed)
c
c       assign initial velocities to atoms
c
	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)
	  x1(i)=xx
	  call random(x,yy)
	  x=yy
	  yy=2.0*(yy-0.5)
	  y1(i)=yy
	  call random(x,zz)
	  x=zz
	  zz=2.0*(zz-0.5)
	  z1(i)=zz
	  sumx=sumx+x1(i)
	  sumy=sumy+y1(i)
	  sumz=sumz+z1(i)
	enddo
c
	SEED=X
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=sqrt(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 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
	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.ge.cube2) x=(xr-cube)*sign(1.0,x)
	  if(yr.ge.cube2) y=(yr-cube)*sign(1.0,y)
	  if(zr.ge.cube2) z=(zr-cube)*sign(1.0,z)
c
c       calculate r(ij)
c
	  rsq=x*x+y*y+z*z
	  rctr=sqrt(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 correct(np,cube)
c
c       corrects for periodic boundary conditions
c
	parameter ( npd=4000 )
	common/pos/x0(npd),y0(npd),z0(npd)
	common/dis/dax(npd),day(npd),daz(npd),x0l(npd),y0l(npd),z0l(npd)
c
	do i=1,np
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 LANGEVIN(NP,SEED)
	PARAMETER ( NPD=4000 )
        COMMON/VEL/X1(NPD),Y1(NPD),Z1(NPD)
        COMMON/FOR/FX(NPD),FY(NPD),FZ(NPD)
        COMMON/THERM/GAMMA,CONSt

	X=SEED
	DO I=1,NP
          CALL RANDOM(X,XX)
          X=XX
	  RF=2.0*(X-0.5)*CONST
          FX(I)=FX(I)-GAMMA*X1(I)+RF
	  CALL RANDOM(X,XX)
	  X=XX
          RF=2.0*(X-0.5)*CONST
          FY(I)=FY(I)-GAMMA*Y1(I)+RF
	  CALL RANDOM(X,XX)
	  X=XX
          RF=2.0*(X-0.5)*CONST
          FZ(I)=FZ(I)-GAMMA*Z1(I)+RF
	ENDDO
	SEED=X

	RETURN
	END


	subroutine eqbrat(sumvel,aheat,np)
c
c       scale velocities during initial time-steps
c
	parameter ( npd=4000 )
	common/vel/x1(npd),y1(npd),z1(npd)
c
	heat=sqrt(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
	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
	parameter ( npd=4000, np6 = npd*6 )
	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 dax,day,daz,fx,fy,fz/np6*0.0/
	data idist/1000*0/
	end

