Added Python version of mc_nvt_poly_lj. Made corresponding changes to…

… documentation. This closes #8.

mc_poly_lj_module.f90

@@ -270,7 +270,7 @@ FUNCTION potential_1 ( ri, di, i, box, j_range ) RESULT ( partial )
                       RETURN               ! Return immediately
                    END IF
 
-                   rmag     = sqrt(rab_sq)
+                   rmag     = SQRT(rab_sq)
                    sr6      = sr2**3
                    sr12     = sr6**2
                    pair%pot = 4.0*(sr12-sr6) + lambda1 + lambda2*rmag ! LJ atom-atom pair potential (force-shifted)

python_examples/CONTENTS.md

@@ -98,8 +98,7 @@ We also supply a constant-pressure version [mc_npt_hs.py](mc_npt_hs.py).
 [mc_zvt_lj.py](mc_zvt_lj.py) and  [mc_lj_module.py](mc_lj_module.py).
 
 ### 4.7 Monte Carlo program using quaternions
-We do not provide a Python version of this example.
-A Fortran version is available.
+[mc_nvt_poly_lj.py](mc_nvt_poly_lj.py) and [mc_poly_lj_module.py](mc_poly_lj_module.py).
 
 ### 4.8 Monte Carlo of hard spherocylinders
 [mc_nvt_sc.py](mc_nvt_sc.py) and [mc_sc_module.py](mc_sc_module.py).

python_examples/GUIDE.md

@@ -645,6 +645,73 @@ we make no attempt to reproduce the tests conducted with the Fortran
 versions
 for the _N=13_ chains.
 
+## Polyatomic Lennard-Jones program
+The program `mc_nvt_poly_lj.py` conducts Monte Carlo simulations of a system of rigid molecules
+composed of Lennard-Jones interaction sites.
+For simplicity the sites are taken to be identical, although the program is easily generalized.
+Molecular orientations are represented by quaternions,
+which are used to calculate the rotation matrix
+and hence the interaction site positions.
+
+We test this with the three-site model of orthoterphenyl, a fragile glassformer,
+described in the following publications amongst others.
+
+* LJ Lewis, G Wahnstrom, _Sol State Commun,_ __86,__ 295 (1993)
+* LJ Lewis, G Wahnstrom, _Phys Rev E,_ __50,__ 3865 (1994)
+* S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, _Phys Rev E,_ __65,__ 041205 (2002)
+* E La Nave, S Mossa, F Sciortino, P Tartaglia, _J Chem Phys,_ __120,__ 6128 (2004)
+
+We compare with the results of Mossa et al (2002).
+The sites are arranged at the vertices of an isosceles triangle with bond angle 75 degrees,
+LJ parameters &epsilon; = 5.276 kJ mol<sup>-1</sup>,
+&sigma;=0.483nm,
+and two equal bonds of length &sigma;.
+The program employs the usual reduced units based on &epsilon; and &sigma;
+and in these units the potential cutoff of Mossa et al (2002) is _R_<sub>c</sub>=2.612;
+the pair potential is Lennard-Jones with a shifted-force correction term, linear in _r_,
+to make the potential and its derivative vanish at _r_=_R_<sub>c</sub>.
+Apart from the temperatures
+and much longer runs (10 blocks of 20000 steps each, the same as the Fortran tests),
+default program parameters were used throughout.
+
+Tests were performed at &rho;=0.32655 which is equivalent to &rho;<sub>4</sub>=1.108g cm<sup>-3</sup>
+in Mossa et al (2002).
+Comparisons of potential energy (_PE_=_E_-3 _T_ converted to kJ/mol with a factor 5.276)
+were made with the fit given by eqn (23) of that paper.
+Note that &epsilon;/k<sub>B</sub>&asymp;635 K.
+
+_T_   | _E_        | _P_       | _T_ (K) | _PE_ (kJ/mol) | _PE_ (kJ/mol) eqn (23)
+----- | -----      | -----     | -----   | -----         | -----
+0.5   | -12.995(2) |  1.763(9) |  317    | -76.48(2)     | -76.945
+1.0   | -9.813(15) |  5.86(4)  |  635    | -67.60(8)     | -67.634
+1.5   | -6.87(1)   |  9.33(2)  |  952    | -59.99(5)     | -60.011
+2.0   | -4.10(1)   | 12.30(3)  | 1270    | -53.29(5)     | -53.265
+
+A second set of tests was performed at _T_=0.6&asymp;380K
+at the specified densities &rho;<sub>1</sub>, &hellip; &rho;<sub>5</sub> of Mossa et al (2002).
+Here the excess pressure (_P_(ex)=_P_-&rho;_T_ converted to MPa
+with a factor 77.75 based on the values of &epsilon; and &sigma;)
+is compared with the fit given by eqn (28) and the coefficients in Table III of Mossa et al (2002).
+NB the volumes to insert into the equation are those of their Table I,
+which are specific to their system size.
+In addition their eqn (29) with coefficients in Table V is a fit to their potential energy,
+which we calculate from the simulation as described above.
+
+Id    | &rho;   | _E_        | _P_     | _P_(ex) (MPa) | _P_(ex) (MPa) eqn (28) | _PE_ (kJ/mol) | _PE_ (kJ/mol) eqn (29)
+----- | -----   | -----      | -----   | -----         | -----                  | ------        | ------
+1     | 0.30533 | -11.617(7) | 0.47(2) | 22(2)         | 19.077                 | -70.79(4)     | -70.818
+2     | 0.31240 | -11.927(4) | 1.00(2) | 63(2)         | 60.143                 | -72.42(2)     | -72.289
+3     | 0.31918 | -12.16(1)  | 1.66(3) | 114(2)        | 112.798                | -73.65(5)     | -73.601
+4     | 0.32655 | -12.393(6) | 2.52(2) | 181(2)        | 177.222                | -74.88(3)     | -74.886
+5     | 0.33451 | -12.498(3) | 3.82(1) | 281(1)        | 253.510                | -75.44(2)     | -75.825
+
+In making these comparisons,
+our limited run length (10 blocks of 20000 sweeps each) should be borne in mind,
+since this system can show sluggish behaviour.
+The MD simulations of Mossa et al (2002) are reported to extend to several hundred nanoseconds
+(of order 10<sup>7</sup> MD timesteps) at the lowest temperatures.
+It should be noted that this Python code is quite slow compared to the simpler atomic LJ examples.
+
 ## DPD program
 For testing we compare with an approximate DPD equation of state for _P_.
 

python_examples/maths_module.py

@@ -52,11 +52,9 @@ def random_perpendicular_vector ( old ):
     # Note that we do not require the reference vector to be of unit length
     # However we do require its length to be greater than a small tolerance!
 
-    tol = 1.e-6
-
     assert old.size==3, 'Error in old vector dimension'
     norm = np.sum ( old**2 ) # Old squared length
-    assert norm>tol, "{}{:15.3e}{:15.3e}".format('old normalization error', norm, tol)
+    assert not np.isclose(norm,0.0,atol=1.e-6), 'old too small {} {} {}'.format(*old)
     n = old / np.sqrt(norm) # Normalized old vector
 
     while True: # Loop until generated vector is not too small
@@ -90,6 +88,19 @@ def random_quaternion():
     f = np.sqrt ( ( 1.0 - norm1 ) / norm2 )
     return np.array ( ( zeta[0], zeta[1], beta[0]*f, beta[1]*f ), dtype=np.float_ ) # Random quaternion
 
+def random_rotate_quaternion ( angle_max, old ):
+    """Returns a unit quaternion rotated by a maximum angle (in radians) relative to the old quaternion."""
+    import numpy as np
+
+    # Note that the reference quaternion should be normalized and we test for this
+    assert old.size==4, 'Error in old quaternion dimension'
+    assert np.isclose(np.sum(old**2),1.0), 'old normalization error {} {} {} {}'.format(*old)
+
+    axis = random_vector()                             # Choose random unit vector
+    angle = ( 2.0*np.random.rand() - 1.0 ) * angle_max # Uniform random angle in desired range
+    e = rotate_quaternion ( angle, axis, old )         # General rotation function
+    return e
+
 def random_translate_vector ( dr_max, old ):
     """Returns a vector translated by a random amount."""
 
@@ -105,17 +116,13 @@ def random_rotate_vector ( angle_max, old ):
     """Returns a vector rotated by a small amount relative to the old one."""
 
     import numpy as np
-    
+
     # A small randomly chosen vector is added to the old one, and the result renormalized
     # Provided angle_max is << 1, it is approximately the maximum rotation angle (in radians)
     # The magnitude of the rotation is not uniformly sampled, but this should not matter
 
     # Note that the old vector should be normalized and we test for this
-
-    tol = 1.e-6
-
-    norm = np.sum ( old**2 ) # Old squared length
-    assert np.fabs ( norm - 1.0 ) < tol, "{}{:20.8e}{:20.8e}".format('old normalization error', norm, tol)
+    assert np.isclose(np.sum(old**2),1.0), 'old normalization error {} {} {}'.format(*old)
 
     # Choose new orientation by adding random small vector
     e    = old + angle_max * random_vector ()
@@ -141,18 +148,17 @@ def metropolis ( delta ):
 
 def rotate_vector ( angle, axis, old ):
 
-    """Returns a vector rotated from the old one by specified angle about specified axis."""
+    """Returns a vector rotated from the old one by angle about axis."""
 
     import numpy as np
-    from math import isclose
 
     # Note that the axis vector should be normalized and we test for this
     # In general, the old vector need not be normalized, and the same goes for the result
     # although quite often in our applications they will be
 
     assert old.size == 3, 'Incorrect size of old'
     assert axis.size == 3, 'Incorrect size of axis'
-    assert isclose(np.sum(axis**2),1.0), 'Non-unit vector {} {} {}'.format(*axis)
+    assert np.isclose(np.sum(axis**2),1.0), 'Non-unit vector {} {} {}'.format(*axis)
 
     c    = np.cos ( angle )
     s    = np.sin ( angle )
@@ -163,6 +169,38 @@ def rotate_vector ( angle, axis, old ):
 
     return e
 
+def rotate_quaternion ( angle, axis, old ):
+    """Returns a quaternion rotated by angle about axis relative to old quaternion."""
+
+    import numpy as np
+
+    # Note that the axis vector should be normalized and we test for this
+    # In general, the old quaternion need not be normalized, and the same goes for the result
+    # although in our applications we only ever use unit quaternions (to represent orientations)
+    assert old.size==4, 'Error in old quaternion dimension'
+    assert axis.size==3, 'Error in axis dimension'
+    assert np.isclose (np.sum(axis**2),1.0), 'axis normalization error {} {} {}'.format(*axis)
+
+    # Standard formula for rotation quaternion, using half angles
+    rot = np.sin(0.5*angle) * axis
+    rot = np.array([np.cos(0.5*angle),rot[0],rot[1],rot[2]],dtype=np.float_)
+
+    e = quatmul ( rot, old ) # Apply rotation to old quaternion
+    return e
+
+def quatmul ( a, b ):
+    """Returns quaternion product of two supplied quaternions."""
+
+    import numpy as np
+
+    assert a.size==4, 'Error in a dimension'
+    assert b.size==4, 'Error in b dimension'
+
+    return np.array ( [ a[0]*b[0] - a[1]*b[1] - a[2]*b[2] - a[3]*b[3],
+                        a[1]*b[0] + a[0]*b[1] - a[3]*b[2] + a[2]*b[3],
+                        a[2]*b[0] + a[3]*b[1] + a[0]*b[2] - a[1]*b[3],
+                        a[3]*b[0] - a[2]*b[1] + a[1]*b[2] + a[0]*b[3] ], dtype=np.float_ )
+
 def nematic_order ( e ):
     """Returns a nematic orientational order parameter."""
 
@@ -187,3 +225,25 @@ def nematic_order ( e ):
 
     evals = np.linalg.eigvalsh(q)
     return evals[2]
+
+def q_to_a ( q ):
+    """Returns a 3x3 rotation matrix calculated from supplied quaternion."""
+
+    import numpy as np
+
+    # The rows of the rotation matrix correspond to unit vectors of the molecule in the space-fixed frame
+    # The third row  a(3,:) is "the" axis of the molecule, for uniaxial molecules
+    # Use a to convert space-fixed to body-fixed axes thus: db = np.dot(a,ds)
+    # Use transpose of a to convert body-fixed to space-fixed axes thus: ds = np.dot(db,a)
+
+    # The supplied quaternion should be normalized and we check for this
+    assert np.isclose(np.sum(q**2),1.0), 'quaternion normalization error {} {} {} {}'.format(*q)
+
+    # Write out row by row, for clarity
+    a = np.empty( (3,3), dtype=np.float_ )
+    a[0,:] = [ q[0]**2+q[1]**2-q[2]**2-q[3]**2,   2*(q[1]*q[2]+q[0]*q[3]),       2*(q[1]*q[3]-q[0]*q[2])     ]
+    a[1,:] = [     2*(q[1]*q[2]-q[0]*q[3]),   q[0]**2-q[1]**2+q[2]**2-q[3]**2,   2*(q[2]*q[3]+q[0]*q[1])     ]
+    a[2,:] = [     2*(q[1]*q[3]+q[0]*q[2]),       2*(q[2]*q[3]-q[0]*q[1]),   q[0]**2-q[1]**2-q[2]**2+q[3]**2 ]
+
+    return a
+