wca.a68
1 COMMENT
2
3 @section Synopsis
4
5 Compute equilibrium Lennard Jones thermodynamic data.
6
7 More accurate schemes exist but this suffices to calibrate simulation software.
8 Literature:
9 J.D. Weeks, D. Chandler and H.C. Andersen, J. Chem. Phys. 54(12) 5237
10 L. Verlet and J.J. Weis, Phys. Rev. A 5(2) 939
11 R.J. Baxter, J. Chem. Phys. 52 4559
12 J.W. Perram, Mol. Phys. 30 1505
13
14 COMMENT
15
16 INT resolt = 512, # Channels for rdf #
17 REAL r zero = 2^(1/6); # Radius where f=0 #
18 REAL r max = r zero; # Maximum radius for rdf #
19
20 OP R = (INT n) REAL: (n - 0.5) / resolt * r max,
21 N = (REAL r) INT: ENTIER (1 + resolt * r / r max);
22
23 OP U = (REAL r) REAL: (REAL r6 = r^-6; 4 * r6 * (r6 - 1)),
24 UR = (REAL r) REAL: (r < r zero | U r + 1 | 0),
25 F = (REAL r) REAL: (REAL r6 = r^-6; 24 * r6 / r * (2 * r6 - 1)),
26 G = (REAL r) REAL: rdf[(N r < resolt | N r | resolt)];
27
28 PROC simpson = (REAL a, b, PROC (REAL) REAL f, REAL epsilon) REAL:
29 BEGIN
30 # Depending on error term, integrate recursively #
31 PROC adaptive = (REAL a, m, b, # and values # f a, f m, f b) REAL:
32 IF REAL l = (a + m) / 2, r = (m + b) / 2;
33 REAL f l = f(l), f r = f(r);
34 (b - a) *
35 ABS (f a - 4 * f l + 6 * f m - 4 * f r + f b) / 180 > epsilon
36 THEN adaptive (a, l, m, f a, f l, f m) +
37 adaptive (m, r, b, f m, f r, f b)
38 ELSE (b - a) * (f a + 4 * f l + 2 * f m + 4 * f r + f b) / 12
39 FI;
40 REAL m = (a + b) / 2;
41 adaptive (a, m, b, f(a), f(m), f(b))
42 END;
43
44 REAL ref rho, ref d, [resolt] REAL rdf;
45
46 # Read parameters #
47
48 write(("rho = "));
49 REAL rho := read real;
50 write(("kT = "));
51 REAL beta := 1 / read real;
52
53 # Hard sphere radius is first order Baxter term #
54
55 ref d := simpson(0.0, r max,
56 (REAL r) REAL: 1 - (r < 0.7 | 0 | exp(-beta * UR r)),
57 small real);
58 write(("d = ", fixed(ref d, 0, 6), newline));
59 ref rho := rho * ref d^3;
60 INT chns = 610, INT len = 100; REAL delta = 1 / len;
61 [1 .. 150] REAL t, [1 .. 2100] REAL table;
62 REAL eta, cx, a, a1, a2, a3, b1, b2, b3, b, c1, c2, c;
63
64 # Hard sphere RDF by Wiener-Hopf factorisation #
65
66 eta := pi / 6 * ref rho;
67 cx := 1 - eta;
68 a1 := 1 + 2 * eta;
69 a2 := eta * a1;
70 a3 := a1 / cx^2 - 3 * a2 * delta / cx^3;
71 a := 0.5 * delta * a3;
72 b1 := -3 * eta;
73 b2 := cx * eta - a2;
74 b3 := b1 - 3 * delta * b2 / cx + delta * a1;
75 b := 0.5 * delta * b3 / cx^2;
76 c1 := -b1;
77 c2 := 0.5 * delta * c1 / cx - 1;
78 c := 0.5 * delta * c2 / cx;
79 FOR m TO len
80 DO REAL r = (m - 0.5) * delta;
81 t[m] := c + b * r + a * r^2
82 OD;
83
84 # Fill table #
85
86 FOR n TO chns
87 DO REAL sum := 0;
88 FOR m TO len
89 DO sum +:= (n > m | table[n - m] | 0.5 - n + m - len) * t[m]
90 OD;
91 table[n] := 2 * pi * ref rho * sum
92 OD;
93
94 # Construct RDF #
95
96 FOR m TO UPB rdf
97 DO IF REAL r = R m / ref d; REAL z := r - 1.0;
98 z < 0
99 THEN rdf[m] := 0
100 ELIF z >= 30
101 THEN rdf[m] := 1
102 ELSE z := z / delta + 0.5;
103 INT n = (ENTIER z = 0 | 1 | ENTIER z);
104 REAL y = z - n;
105 REAL h low = table[n] / (len + n - 0.5),
106 h hig = table[n + 1] / (len + n + 0.5);
107 rdf[m] := 1 + h low * (1 - y) + h hig * y
108 FI
109 OD;
110
111 # Extend to the core with Percus-Yevick DCF #
112
113 a1 := -((1 + 2 * eta)^2) / (1 - eta)^4;
114 a2 := (1 + 0.5 * eta)^2 / (1 - eta)^4;
115 FOR n TO UPB rdf
116 WHILE REAL r := R n / ref d;
117 r < 1
118 DO rdf[n] := -(a1 + 6 * eta * a2 * r + 0.5 * a1 * eta * r^3)
119 OD;
120
121 # Construct LJ rdf #
122
123 FOR n TO UPB rdf
124 DO REAL r = R n;
125 rdf[n] *:= (r < 0.7 | 0 | exp (-beta * (r < r max | UR r | 0)))
126 OD;
127
128 # Thermodynamic data #
129
130 write(("N = ", fixed(rho * simpson(0.7, r max,
131 (REAL r) REAL: 4 * pi * r^2 * G r, 1e-3), 0, 2), new line));
132 write(("E/N = ", fixed(rho / 2 * beta * simpson(0.7, r max,
133 (REAL r) REAL: 4 * pi * r^2 * G r * U r, 1e-3), 0, 2), new line));
134 write(("p = ", fixed(rho^2 / 6 * beta * simpson(0.7, r max,
135 (REAL r) REAL: 4 * pi * r^2 * G r * r * F r, 1e-3), 0, 2), new line)); # Virial equation #
136
137 open(standout, "rdf.dat", standout channel);
138 FOR i TO resolt
139 DO write((fixed (R i, -10, 6), ",",
140 fixed(G R i, -10, 6), new line))
141 OD
