Algol 68 Genie

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
Deze website is gearchiveerd door de KB, nationale bibliotheek.

This website is archived by the National Library of the Netherlands.