PROGRAM md
!
!Purpose:
!The main program to calculate the thermal conductivity of solid argon in a NVT ensemble.
!
!Record of revisions:
! Date Programmer Description of change
! ==== ========== ======================
! 11/24/2004 Xuan Zheng Force calculation.
! 11/27/2004 Xuan Zheng Modify the initialization subroutine to calculate the
! old position using the randomly generated velocity.
! Add the integrate subroutine.
! 12/04/2004 Xuan Zheng Add common blocks.
! 12/07/2004 Xuan Zheng Fix the error caused by the drift of the mass center
! z direction.
IMPLICIT none
!
!Declare variables.
!qid: id of the run.
!fpar: file name of the paprameter file.
!finitial: file name of the initial conditin output.
!lpx: length of the qid.
!natoms: number of atoms in the box.
!rho: density (in reduced units).
!box: size of the box.
!hbox: half of the box.
!temp: inintial temperature.
!nstep: total number of steps.
!vx(natoms): velocity of atom i in the x direction.
!vy(natoms): velocity of atom i in the y direction.
!vz(natoms): velocity of atom i in the z direction.
!fx(natoms), fy(natoms), fz(natoms): force on every atom along x, y, z directions.
!rcut: curoff radius.
!upot: potential energy.
!press: pressure of the system.
!x(natoms), y(natoms), z(natoms): positions of the atoms.
!xo(natoms), yo(natoms), zo(natoms): old positons of the atoms.
!xn(natoms), yn(natoms), zn(natoms): new positions of the atoms.
!step: the number of the current step.
!nwarm: number of warm up steps.
!ukin: kinetic energy of the system.
!kin(10000): output kinetic energy.
!write_scaler: intervals between write out scalers.
COMMON /block1/ x, y, z
COMMON /block2/ xo, yo, zo
COMMON /block3/ vx, vy, vz, ukin
COMMON /block4/ fx, fy, fz, upot
COMMON /block5/ box, tstep
COMMON /block6/ rcut, point, list
COMMON /block7/ pot, j2x, j2y, j2z
COMMON /block8/ qstep, nquench, quench_times, quench_interval
COMMON /block9/ temp
CHARACTER qid*8, fpar*14, finitial*14, flambda*14, fkin*14,
& fpos*14, fcur*14
INTEGER natoms, lpx, i, j, nstep, step, write_scaler, nscaler,
& quench_interval, quench_times, nquench, qstep, maxnab
PARAMETER (natoms=108)
PARAMETER (maxnab=natoms*150)
INTEGER point(natoms), list(maxnab), mavg, nequ
PARAMETER (mavg=10000, nequ=100000)
DOUBLE PRECISION tstep, rho, box, hbox, temp,
& x(natoms), y(natoms), z(natoms),
& vx(natoms), vy(natoms), vz(natoms),
& fx(natoms), fy(natoms), fz(natoms),
& xo(natoms), yo(natoms), zo(natoms),
& xn(natoms), yn(natoms), zn(natoms),
& rcut, upot, press, ukin,
& kin(100000), pot(natoms), j2x, j2y, j2z, jx, jy, jz,
& jxi(5000000), jyi(5000000), jzi(5000000),
& sigma, eps, kb, mass, realt,
& jijj, volume, temp2, lambda,
& heatcorr(mavg), lambdacoeff, unitconv, potential(100000)
PARAMETER (sigma=3.405D-10)
PARAMETER (eps=1.65398276D-21)
PARAMETER (kb=1.38062D-23, mass=6.67545D-26)
!Read parameters from qid.par.
!
! Format of the qid.par file:
! rcut
! rho
! temp
! nstep
! tstep
! quench_interval
! quench_times
! write_scaler
WRITE (*,*) 'Input the run id:'
READ (*,'(a)') qid
lpx=index(qid,' ')-1
fpar=qid(1:lpx)//'.par'
finitial=qid(1:lpx)//'.ini'
flambda=qid(1:lpx)//'.lbd'
fkin=qid(1:lpx)//'.erg'
fpos=qid(1:lpx)//'.pos'
fcur=qid(1:lpx)//'.cur'
OPEN (UNIT=9,FILE=fpar,STATUS='OLD')
READ (9,*) rcut
WRITE (*,*) rcut
READ (9,*) rho
WRITE (*,*) rho
READ (9,*) temp
WRITE (*,*) temp
READ (9,*) nstep
WRITE (*,*) nstep
READ (9,*) tstep
WRITE (*,*) tstep
READ (9,*) quench_interval
WRITE (*,*) quench_interval
READ (9,*) quench_times
WRITE(*,*) quench_times
READ(9,*) write_scaler
WRITE(*,*) write_scaler
! Calculate the box size and half of the box.
box=(DBLE(natoms)/rho)**(1.D00/3.D00)
realt=DSQRT(mass*sigma*sigma/eps)
! write(*,*) 'realt=', realt
! write(*,*) 'box', box
hbox=0.5D00*box
nquench=0
step=0
qstep=0
CALL initial
WRITE (*,*) 'Initialization done.'
! write (*,*) 'temp after initial',temp
CALL neighborlist
WRITE (*,*) 'Neighbor list generated.'
! write (*,*) 'temp after neighbor',temp
! call force
! write(*,*) temp
! call integrate
! write(*,*) temp
WRITE(*,*) 'MD loop.....'
!MD loop.
DO i=1, nstep
step=step+1
qstep=qstep+1
! write(*,*) pot(1)
CALL force
CALL heatcurrent(jx, jy, jz)
! write(*,*) pot(1)
jxi(i)=jx
jyi(i)=jy
jzi(i)=jz
! write(*,*) temp
CALL integrate
! write(*,*) temp
nscaler=(step/write_scaler)+1
kin(nscaler)=ukin
! write(*,*) upot
potential(nscaler)=upot
ENDDO
!calculate the thermal conductivity.
! DO i=1, natoms
! write(*,*) pot(i)
! ENDDO
jijj=0.0D0
WRITE(*,*) 'MD loop done.'
DO i=1, mavg
heatcorr(i)=0.0D00
DO j=nequ, nstep-i
heatcorr(i)=heatcorr(i)+jxi(i+j)*jxi(j)
& +jyi(i+j)*jyi(j)+jzi(i+j)*jzi(j)
ENDDO
heatcorr(i)=heatcorr(i)/(nstep-i-nequ)
jijj=jijj+heatcorr(i)
ENDDO
WRITE (*,*) 'Heat correlation function generated.'
! write(*,*) jijj
unitconv=eps*eps*sigma*sigma/realt/realt
jijj=jijj*unitconv
volume=box*box*box*sigma*sigma*sigma
temp2=temp*temp*eps*eps/kb/kb
! write(*,*) 'temp2', temp2
lambdacoeff=tstep*realt/volume/temp2/kb/3
lambda=lambdacoeff*jijj
write(*,*) 'thermal conductivity', lambda
DO i=1, mavg
heatcorr(i)=unitconv*lambdacoeff*heatcorr(i)
ENDDO
OPEN (UNIT=11,FILE=flambda,STATUS='unknown')
WRITE(11,*) 'temperature=', temp, temp*119.8, 'K'
WRITE(11,*) 'thermal conductivity=', lambda, 'W m-1 K-1'
DO i=1, mavg
WRITE (11, *) i, heatcorr(i)
ENDDO
! potential=0.0D0
OPEN (UNIT=10,FILE=fkin,STATUS='unknown')
DO i=1, nscaler
! DO j=1, natoms
! potential=pot(j)+potential
! ENDDO
! potential=0.5D0*potential
WRITE (10,*) i, kin(i), potential(i)
! potential=0.0D0
ENDDO
OPEN (UNIT=12,FILE=fpos,STATUS='unknown')
WRITE(12,*) natoms
WRITE(12,*) box
DO i=1, natoms
WRITE (12, *) x(i), y(i), z(i)
ENDDO
OPEN (UNIT=13,FILE=fcur,STATUS='unknown')
DO i=nequ, nequ+1000
WRITE (13, *) i, jxi(i)
ENDDO
END
! End of the main program.
! Begin of the subroutine to initialize the system.
SUBROUTINE initial
!
!Purpose:
!To distribute the atoms in a cubic lattice, and assign every atom a velocity.
!The velocity is normalized to the desired temperature.
!Calculate the previous position using the generated velocity.
!
IMPLICIT none
COMMON /block1/ x, y, z
COMMON /block2/ xo, yo, zo
COMMON /block3/ vx, vy, vz, ukin
COMMON /block5/ box, tstep
COMMON /block9/ temp
INTEGER natoms
PARAMETER (natoms=108)
DOUBLE PRECISION box, tstep
DOUBLE PRECISION x(natoms), y(natoms), z(natoms),
& vx(natoms), vy(natoms), vz(natoms), ukin,
& xo(natoms), yo(natoms), zo(natoms)
INTEGER n, iref, ix, iy, iz, m, i
DOUBLE PRECISION del, v2, vx0, vy0, vz0, vx0t, vy0t, vz0t,
& temp, rangauss, f
del=(box**3)**(1.D00/3.D00)
! write (*,*) box, del
! WRITE(*,*) natoms
n=DNINT((DBLE(natoms)*0.25D00)**(1.D00/3.D00))
! write (*,*) n
del=0.5D00*del/DBLE(n)
! write(*,*) 'del=', del
!sublattice 1.
x(1)=0.0D00
y(1)=0.0D00
z(1)=0.0D00
!sublattice 2.
x(2)=del
y(2)=del
z(2)=0.0D0
!sublattice 3.
x(3)=0.0D00
y(3)=del
z(3)=del
!sublattice 4.
x(4)=del
y(4)=0.0D00
z(4)=del
m=0
! Put the atoms into a cubic lattice.
DO iz=1, n
DO iy=1, n
DO ix=1, n
DO iref=1,4
x(iref+m)=x(iref)+2.0D00*del*DBLE(ix-1)
y(iref+m)=y(iref)+2.0D00*del*DBLE(iy-1)
z(iref+m)=z(iref)+2.0D00*del*DBLE(iz-1)
ENDDO
m=m+4
END DO
END DO
END DO
! Set up initial velocities of the atoms.
vx0=0.D0
vy0=0.D0
vz0=0.D0
v2=0.0D00
DO i=1, natoms
vx(i)=rangauss()
vy(i)=rangauss()
vz(i)=rangauss()
vx0=vx0+vx(i)
vy0=vy0+vy(i)
vz0=vz0+vz(i)
v2=v2+vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)
ENDDO
! Set center of mass movement to zero. And scale to the desired temperature.
vx0=vx0/DBLE(natoms)
vy0=vy0/DBLE(natoms)
vz0=vz0/DBLE(natoms)
vx0t=0.D0
vy0t=0.D0
vz0t=0.D0
f = DSQRT(3.0D0*DBLE(natoms)*temp/v2)
! write(*,*) 'temp=', temp
v2=0.D0
DO i=1, natoms
vx(i)=(vx(i)-vx0)*f
vy(i)=(vy(i)-vy0)*f
vz(i)=(vz(i)-vz0)*f
vx0t=vx0t+vx(i)
vy0t=vy0t+vy(i)
vz0t=vz0t+vz(i)
v2=v2+vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)
ENDDO
v2=v2/3.0D0/DBLE(natoms)
vx0t=vx0t/natoms
vy0t=vy0t/natoms
vz0t=vz0t/natoms
! write(*,*) 'inittemp', v2
! write(*,*) 'velocity of center', vx0t, vy0t, vz0t
!Calcuate the previous position using the generated velocity.
DO i=1, natoms
xo(i)=x(i)-tstep*vx(i)
yo(i)=y(i)-tstep*vy(i)
zo(i)=z(i)-tstep*vz(i)
! write(*,*) xo(i), yo(i), zo(i)
! write(*,*) x(i), y(i), z(i)
ENDDO
END
!End of initialization.
!
SUBROUTINE neighborlist
!purpose:
!subroutine to calculate the verlet neighbor list. Because it is solid,
!the neighbor list do not need update.
!
IMPLiCIT none
COMMON /block1/ x, y, z
COMMON /block5/ box, tstep
COMMON /block6/ rcut, point, list
INTEGER natoms, maxnab
PARAMETER (natoms=108, maxnab=natoms*150)
INTEGER point(natoms), list(maxnab)
DOUBLE PRECISION x(natoms), y(natoms), z(natoms), box, tstep, rcut
DOUBLE PRECISION xi, yi, zi, dx, dy, dz, d2, boxi, rlist2
INTEGER i, nlist,j
boxi=1/box
rlist2=(rcut+0.1D0)*(rcut+0.1D0)
nlist=0
WRITE(*,*) 'Generating neighbor list.......'
DO i=1, natoms-1
point(i)=nlist+1
xi=x(i)
yi=y(i)
zi=z(i)
! WRITE(*,*) 'neighbor list'
DO j=i+1, natoms
dx=xi-x(j)
dy=yi-y(j)
dz=zi-z(j)
! write (*,*) dx
dx=dx-DNINT(dx*boxi)*box
dy=dy-DNINT(dy*boxi)*box
dz=dz-DNINT(dz*boxi)*box
d2=dx*dx+dy*dy+dz*dz
! write(*,*) d2
IF (d2 .LT. rlist2) THEN
nlist=nlist+1
! write(*,*) nlist
list(nlist)=j
ENDIF
ENDDO
ENDDO
! WRITE (*,*) 'Neighbor list done'
point(natoms)=nlist+1
END
!
!
SUBROUTINE force
!
!Purpose:
!To calculate the force, the potential energy, and the heat current.
!
!
IMPLICIT none
COMMON /block1/ x, y, z
COMMON /block3/ vx, vy, vz, ukin
COMMON /block4/ fx, fy, fz, upot
COMMON /block5/ box, tstep
COMMON /block6/ rcut, point, list
COMMON /block7/ pot, j2x, j2y, j2z
!Declare variables:
!fx(i)
!fy(i)
!fz(i)
!upot: potential energy of the system.
!press: pressure of the system.
!dx, dy, dz: distance of two atoms in the x, y, z direction, respectively.
!d2: dx^2+dy^2+dz^2.
!d2i: inverse of d2.
!d6i=d2i^3
!rcutoff: cutoff radius
!rcutoff2: square of the cutoff radius.
!ecut: potential energy at the cutoff radius.
!
INTEGER natoms, i, j, jbeg, jend, jnab, maxnab
PARAMETER (natoms=108, maxnab=natoms*150)
INTEGER point(natoms), list(maxnab)
DOUBLE PRECISION x(natoms), y(natoms), z(natoms), xi, yi, zi
DOUBLE PRECISION fx(natoms), fy(natoms), fz(natoms)
DOUBLE PRECISION vx(natoms), vy(natoms), vz(natoms)
DOUBLE PRECISION upot, virial, press, rcut, rcut2, rcut2i, rcut6i
DOUBLE PRECISION dx, dy, dz, d2, d2i, d6i, d12i, ff, ecut, tstep
DOUBLE PRECISION box, hbox, boxi, uij, vij, ukin
DOUBLE PRECISION pot(natoms), j2x, j2y, j2z, fvx, fvy, fvz
rcut2=rcut*rcut
rcut2i=1.0D0/rcut2
rcut6i=rcut2i*rcut2i*rcut2i
ecut=4.0D00*rcut6i*(rcut6i-1.0D00)
boxi=1/box
!Set forces, potential energy and pressure to zero.
DO i=1, natoms
fx(i)=0.0D0
fy(i)=0.0D0
fz(i)=0.0D0
pot(i)=0.0D00
END DO
upot=0.0D00
virial=0.0D00
press=0.0D0
j2x=0.0D0
j2y=0.0D0
j2z=0.0D0
DO i=1, natoms-1
jbeg=point(i)
jend=point(i+1)-1
IF (jbeg .LE. jend) THEN
xi=x(i)
yi=y(i)
zi=z(i)
DO jnab=jbeg, jend
j=list(jnab)
! write(*,*) j
dx=xi-x(j)
dy=yi-y(j)
dz=zi-z(j)
dx=dx-DNINT(dx*boxi)*box
dy=dy-DNINT(dy*boxi)*box
dz=dz-DNINT(dz*boxi)*box
d2=dx*dx+dy*dy+dz*dz
IF (d2 .LT. rcut2) THEN
d2i=1.0D0/d2
d6i=d2i*d2i*d2i
d12i=d6i*d6i
uij=4.0D00*(d12i-d6i)
upot=upot+uij-ecut
vij=uij+4.0D0*d12i
virial=virial+vij
ff=6.0D00*vij*d2i
pot(i)=pot(i)+uij-ecut
pot(j)=pot(j)+uij-ecut
fx(i)=fx(i)+ff*dx
fy(i)=fy(i)+ff*dy
fz(i)=fz(i)+ff*dz
fx(j)=fx(j)-ff*dx
fy(j)=fy(j)-ff*dy
fz(j)=fz(j)-ff*dz
!second term of j(t).
fvx=dx*ff*(vx(i)+vx(j))
fvy=dy*ff*(vy(i)+vy(j))
fvz=dz*ff*(vz(i)+vz(j))
j2x=j2x+dx*(fvx+fvy+fvz)
j2y=j2y+dy*(fvx+fvy+fvz)
j2z=j2z+dz*(fvx+fvy+fvz)
ENDIF
ENDDO
ENDIF
ENDDO
END
SUBROUTINE heatcurrent(jx, jy, jz)
!
!purpose:
!to calculate the heat current.
!
IMPLICIT none
COMMON /block3/ vx, vy, vz, ukin
COMMON /block4/ fx, fy, fz, upot
COMMON /block5/ box, tstep
COMMON /block7/ pot, j2x, j2y, j2z
INTEGER natoms, i
PARAMETER (natoms=108)
DOUBLE PRECISION vx(natoms), vy(natoms), vz(natoms), ukin,
& fx(natoms), fy(natoms), fz(natoms), tstep, pot(natoms),
& j2x, j2y, j2z, jx, jy, jz, e(natoms), box, upot
jx=0.0D0
jy=0.0D0
jz=0.0D0
DO i=1, natoms
e(i)=0.0D0
ENDDO
DO i=1, natoms
e(i)=vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)
e(i)=e(i)+pot(i)
e(i)=0.5D0*e(i)
! write(*,*) pot(i)
jx=jx+vx(i)*e(i)
jy=jy+vy(i)*e(i)
jz=jz+vz(i)*e(i)
ENDDO
jx=jx+0.5D0*j2x
jy=jy+0.5D0*j2y
jz=jz+0.5D0*j2z
END
!
SUBROUTINE integrate
!
!Purpose:
!To integrate the equation of motion using Verlet algorithm.
!
IMPLICIT none
COMMON /block1/ x, y, z
COMMON /block2/ xo, yo, zo
COMMON /block3/ vx, vy, vz, ukin
COMMON /block4/ fx, fy, fz, upot
COMMON /block5/ box, tstep
COMMON /block8/ qstep, nquench, quench_times, quench_interval
COMMON /block9/ temp
INTEGER i, step, natoms, nwarm, nquench, qstep, quench_times
PARAMETER (natoms=108)
INTEGER quench_interval
DOUBLE PRECISION x(natoms), y(natoms), z(natoms)
DOUBLE PRECISION xo(natoms), yo(natoms), zo(natoms)
DOUBLE PRECISION fx(natoms), fy(natoms), fz(natoms)
DOUBLE PRECISION xn(natoms), yn(natoms), zn(natoms)
DOUBLE PRECISION vx(natoms), vy(natoms), vz(natoms), upot
DOUBLE PRECISION ukin, tstep, tstep2, i2tstep, scale, box, temp
ukin=0.0D0
tstep2=tstep*tstep
i2tstep=0.5D0/tstep
!Integrate the equation of motion using Verlet algorithm.
DO i=1, natoms
xn(i)=2.0D0*x(i)-xo(i)+fx(i)*tstep2
yn(i)=2.0D0*y(i)-yo(i)+fy(i)*tstep2
zn(i)=2.0D0*z(i)-zo(i)+fz(i)*tstep2
vx(i)=(xn(i)-xo(i))*i2tstep
vy(i)=(yn(i)-yo(i))*i2tstep
vz(i)=(zn(i)-zo(i))*i2tstep
!write(*,*) vx(i), vy(i), vz(i)
ukin=ukin+vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)
! write (*,*) 'ukin=', ukin
ENDDO
! write(*,*) zn(5), zo(5), i2tstep
!For the warm up steps, use the velocity scaling to get the exact temperature.
IF(qstep.EQ. quench_interval .AND. nquench .LE. quench_times) THEN
! write(*,*) ukin
scale=DSQRT(temp*3.0D0*DBLE(natoms)/ukin)
qstep=0
nquench=nquench+1
ELSE
scale=1.0D0
ENDIF
! write (*,*) temp, scale
ukin=0.0D0
DO i=1, natoms
!Scale the velocity of the atoms to the initial temperature.
vx(i)=scale*vx(i)
vy(i)=scale*vy(i)
vz(i)=scale*vz(i)
xn(i)=xo(i)+2.0D0*vx(i)*tstep
yn(i)=yo(i)+2.0D0*vy(i)*tstep
zn(i)=zo(i)+2.0D0*vz(i)*tstep
ukin=ukin+vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)
xo(i)=x(i)
yo(i)=y(i)
zo(i)=z(i)
x(i)=xn(i)
y(i)=yn(i)
z(i)=zn(i)
!Put the atoms back into the box. If the new positons are moved, the old positions must be moved
!too.
IF (x(i) .GT. box) THEN
x(i)=x(i)-box
xo(i)=xo(i)-box
ELSEIF (x(i) .LT. 0.0D0) THEN
x(i)=x(i)+box
xo(i)=xo(i)+box
ENDIF
IF (y(i) .GT. box) THEN
y(i)=y(i)-box
yo(i)=yo(i)-box
ELSEIF (y(i) .LT. 0.0D0) THEN
y(i)=y(i)+box
yo(i)=yo(i)+box
ENDIF
IF (z(i) .GT. box) THEN
z(i)=z(i)-box
zo(i)=zo(i)-box
ELSEIF (z(i) .LT. 0.0D0) THEN
z(i)=z(i)+box
zo(i)=zo(i)+box
ENDIF
ENDDO
END
FUNCTION ran_uniform()
!
!Purpose:
!To generate a serial of uniformly distributed random numbers with the seed idum.
!
IMPLICIT none
INTEGER idum
DOUBLE PRECISION ran_uniform,ran2
SAVE idum
DATA idum/-10/
ran_uniform = ran2(idum)
RETURN
END
FUNCTION ran2(idum)
IMPLICIT none
INTEGER idum,im1,im2,imm1,ia1,ia2,
& iq1,iq2,ir1,ir2,ntab,ndiv
DOUBLE PRECISION ran2,am,eps,rnmx
PARAMETER (im1=2147483563,im2=2147483399,
& am=1.0d0/im1,imm1=im1-1,ia1=40014,
& ia2=40692,iq1=53668,iq2=52774,ir1=12211,
& ir2=3791,ntab=32,ndiv=1+imm1/ntab,
& eps=1.2d-7,rnmx=1.0d0-eps)
INTEGER idum2,j,k,iv(Ntab),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/ntab*0/, iy/0/
IF (idum.LE.0) THEN
idum=MAX(-idum,1)
idum2=idum
DO j=ntab+8,1,-1
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
IF (idum.LT.0) idum=idum+im1
IF (j.LE.ntab) iv(j)=idum
ENDDO
iy=iv(1)
ENDIF
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
IF (idum.LT.0) idum=idum+im1
k=idum2/iq2
idum2=ia2*(idum2-k*iq2)-k*ir2
IF (idum2.LT.0) idum2=idum2+im2
j=1+iy/ndiv
iy=iv(j)-idum2
iv(j)=idum
IF(iy.LT.1)iy=iy+imm1
ran2=Min(am*iy,rnmx)
RETURN
END
FUNCTION rangauss()
!
!Purpose:
!To generate random numbers from a gaussian distribution.
!
DOUBLE PRECISION ran_uniform, rangauss, v1, v2, rsq
100 v1=2.0D0*ran_uniform()-1.0D0
v2=2.0D0*ran_uniform()-1.0D0
rsq=v1*v1+v2*v2
IF (rsq .GE. 1.0D0 .OR. rsq .LE. 0.0D0) GOTO 100
rangauss=v1*DSQRT(-2.0D0*DLOG(rsq)/rsq)
RETURN
END