## 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
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

```