Physical Adsorption Project
Physical Adsorption Project
Click to ftp this document in Text Format
Current Free Energy Fortran Program
PROGRAM JPFEN1
c latest version: 22 Sept 1996, programmers JP and JAV
c For neon/graphite W/ Data output of Energy (F), Entropy (S),
c Specific Heat (C), Pressure (P). Starting program used for Xe and Kr
c in the late 70' and early '80's by J.A. Venables and coworkers
c (J. Physique 38 (1977) C4-105; Surf. Sci. 105 (1981) 536).
c
DIMENSION EM(36,45,6),T(36),X1(36),ZFE(3,36),TH(3),
+ EP(3,36),EMU(36)
Character FNAME*20
REAL*8 H,K,M,LAM
COMMON/SQRTS/SQR2,SQR3,SQR7,SQR12,SQR28,SQI12,SQI28,PI
COMMON/R1/R1,CONS,N1
c
Write(*,40)
40 Format(' Output file:')
Call namget(fname)
Open(unit=7,file=fname,status='new')
WRITE(*,50)
Write(7,50)
50 FORMAT(1x,'NEON/GRAPHITE LENNARD JONES +CELL MODEL')
c
PI=3.1415926536
SQR2=SQRT(2.)
SQR3=SQRT(3.)
SQR7=SQRT(7.)
SQR12=SQRT(12.)
SQR28=SQRT(28.)
SQI12=1./SQR12
SQI28=1./SQR28
c
c Atomic and Potential constants
c
N1=100 !grid division for cell model
H=6.626176E-4 !Planck (exponents to avoid overflow)
K=1.380662E+07 !Boltzmann
M=20.183*1.6726E+03 !Neon mass
CONS=0.0 !3-body forces not used in this version
EAV=-496.0 !Average Ne/graphite potential: ref G.L.Price,
R1=3.116 !Ne radius (Angstroms)
TH(1)=400.0 !Surface Sci. 46 (1974) 697-702, table 1
TH(2)=60.0 !Einstein temperature actually used for ZFE
TH(3)=80.0
A1 =3.09 !Lattice parmeters
DA=0.10
T1 = 12 !Temperatures
DT = 0.25
WRITE(*,850)EAV,TH(1),TH(2),TH(3)
WRITE(7,850)EAV,TH(1),TH(2),TH(3)
850 FORMAT(1x,'ADSORPTION ENERGY=',1F6.1,' THETA Z VALUES:',3F6.1)
c
c Temperature loop
c
DO 1000 IT=1,35
T(IT)=T1+(IT-(1))*DT
LAM=H*1.0E+10/SQRT(2.0*PI*M*K*T(IT)) !Thermal de Broglie lambda
RLAM=LAM*LAM
c
c EMU is MU_o, standard free energy of 3D gas (T)
c
EMU(IT)=-T(IT)*(ALOG((T(IT)*K)/(RLAM*LAM)))
c
DO 1500 IK=1,3 !To test effects of Einstein frequency on ZFE
QE=EXP(-TH(IK)/(2*T(IT)))/(1-EXP(-TH(IK)/T(IT)))
VZA=LAM*QE
ZFE(IK,IT)=-T(IT)*ALOG(QE) !Z-free energy, perp to substrate
1500 CONTINUE
c
c Lattice Parameter (a) Loop, 3 values
c
DO 2000 IM=1,3
X1(IM)=A1+(IM-1)*DA
CALL LATDYN(X1(IM),V6,VH2,T(IT),RLAM) !Other argts in common/R1/
EM(IM,IT,1)=V6
EM(IM,IT,2)=VH2
EM(IM,IT,3)=ZFE(IM,IT)
c Free energy (F) as function of lattice parameter and T
EM(IM,IT,4)=V6+VH2+ZFE(2,IT)+EAV !Use TH(2) in ZFE
IF(IT.EQ.1)GOTO 2000
c Entropy (S) by differentiation of F (check sign?)
EM(IM,IT,5)=-(EM(IM,IT,4)-EM(IM,IT-1,4))/(T(IT)-T(IT-1))
IF(IT.EQ.2)GOTO 2000
c Specific Heat (C) by differentiation of S (check sign?)
EM(IM,IT,6)=T(IT-1)*(EM(IM,IT,5)-EM(IM,IT-1,5))/(T(IT)-T(IT-1))
2000 CONTINUE
c
c Pressure from Free Energy F(a) and dF/da. EP is a calculation of
c presssure in Torr, from calculation of Ln(p) in Pascal. Notes of
c June 20, 1996. Uses 3 a-values to fit F data to a parabola locally.
c JP still checking whether values are reasonable.
c
EDF1=(X1(1)/(2.*DA))*(EM(2,IT,4)-EM(1,IT,4))
EP(1,IT)=(DEXP((EM(1,IT,4)-EDF1-EMU(IT))/T(IT)))*(1/133.32)
c
EDF2=(X1(2)/(4.*DA))*(EM(3,IT,4)-EM(1,IT,4))
EP(2,IT)=(DEXP((EM(2,IT,4)-EDF2-EMU(IT))/T(IT)))*(1/133.32)
c
EDF3=(X1(3)/(2.*DA))*(EM(3,IT,4)-EM(2,IT,4))
EP(3,IT)=(DEXP((EM(3,IT,4)-EDF3-EMU(IT))/T(IT)))*(1/133.32)
c
1000 CONTINUE
c
c Write Output Tables
c
DO 4000 II=0,1
IF(II.EQ.1) GOTO 4000
IF(II.EQ.0) WRITE(*,600)
IF(II.EQ.0) WRITE(7,600)
600 FORMAT(///10X,'TABLE OF FREE ENERGY : LATTICE PARAMETER VS T'//)
IF(II.EQ.0) WRITE(*,650)
IF(II.EQ.0) WRITE(7,650)
650 FORMAT(1x,'T',10x,'F1',10x,'F2',10x,'F3',
+ 20x,'S1',10x,'S2',10x,'S3',10x,'C1',10x,
+ 'C2',10x,'C3',10X,'P1',10X,'P2',10X,'P3',10x,'EMU')
DO 3500 IT=1,35
WRITE (*,700)T(IT),(EM(IM,IT,4),IM=1,3),(EM(IM,IT,5),IM=1,3),
+ (EM(IM,IT,6),IM=1,3),(EP(IM,IT),IM=1,3),(EMU(IT))
WRITE (7,700)T(IT),(EM(IM,IT,4),IM=1,3),(EM(IM,IT,5),IM=1,3),
+ (EM(IM,IT,6),IM=1,3),(EP(IM,IT),IM=1,3),(EMU(IT))
700 FORMAT(1X,1F10.3,3F12.4,1x,3F12.4,1x,3F12.4,
+ 1x,3E12.4,1x,1F12.4)
3500 CONTINUE
4000 CONTINUE
STOP
END
c
SUBROUTINE LATDYN(X1,V6,VH2,T,RLAM)
c This s/r does the 2D potential sums (up to ninth neighbors)
c and computes the cell model in 2D parallel to the substrate.
c
DIMENSION X(11),Z(11),RD(4),Y(11)
COMMON/SQRTS/SQR2,SQR3,SQR7,SQR12,SQR28,SQI12,SQI28,PI
COMMON/R1/R1,CONS,N1
DATA Z/2.,4.,4.,2.,2.,4.,4.,4.,4.,2.,4./
c
B=0.0
V4=0.0
V6=0.0
VH2=0.0
X(1)=X1 !X(i) are distances
X(2)=X1
X(3)=X1*SQR3
X(4)=X(3)
X(5)=2.0*X1
X(6)=X(5)
X(7)=X1*SQR7
X(8)=X(7)
X(9)=X(7)
X(10)=0.0
X(11)=0.0
Y(1) =(Z(1)/T)*P0BWL(X(1)) !Y(i) are potential energies
Y(2) =2.0*Y(1) !Z(i) are multiplicities associated
Y(4)=(Z(4)/T)*P0BWL(X(4)) !with the individual X(i) in the
Y(3)=2.0*Y(4) !hexagonal adsorbed layer
Y(5)=(Z(5)/T)*P0BWL(X(5))
Y(6)=2.0*Y(5)
Y(7)=(Z(7)/T)*P0BWL(X(7))
Y(8)=Y(7)
Y(9)=Y(7)
Y(10)=0.0
Y(11)=0.0
c
DO 19 I=1,9
V6=V6+Y(I)*T/2
19 CONTINUE
V6=V6+CONS/(X(1)**9) !Add effect of 3-body forces if present
c
c 2D cell model for vibrations parallel to the substrate
c
DO 21 N=1,N1
RN=FLOAT(N)
R3=(RN-0.5)*R1/FLOAT(N1) !displacement from cell center
RD(1)=R3
RD(2)=R3/2.0
RD(3)=R3*SQR3/2.0
RD(4)=0.0
B=0.0
c
DO 23 I=1,3
V2=X(I)+RD(I)
V3=X(I)-RD(I)
IF(R1/V3-1000.0) 26,26,25
26 B=B+Y(I)-(Z(I)/(2.0*T))*(P0BWL(V2)+P0BWL(V3))
IF(ABS(B)-80.) 23,23,25
23 CONTINUE
c
V4=V4+2.0*PI*EXP(B)*R3*R1/FLOAT(N1)
21 CONTINUE
25 VH2=-T*ALOG(V4/RLAM) !Cell model formula
EN=V6+VH2 !Lateral contribution to free energy
RETURN
END
c
FUNCTION P0BWL(RR)
c Lennard-Jones (12-6) Neon Potential
COMMON/R1/R1,CONS,N1
c R1=3.116 !Ne radius, all neighbor potl, in main program
EPS=36.09 !Depth of Ne-Ne potential well (K), note units
R=R1/RR !Ref: G.K. Horton, in Rare Gas Solids,
R3=R*R*R !vol I, Chapter 1, Table IV, page 88.
R6=R3*R3
P0BWL=EPS*R6*(R6-2.0) !12-6 formula
RETURN
END
c
subroutine namget(fname)
character fname*20
read(*,100) fname
100 format (A20)
return
end