diff --git a/mc_poly_lj_module.f90 b/mc_poly_lj_module.f90 index 53217ae..f0b985f 100644 --- a/mc_poly_lj_module.f90 +++ b/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) diff --git a/python_examples/CONTENTS.md b/python_examples/CONTENTS.md index 8ac1092..a608263 100644 --- a/python_examples/CONTENTS.md +++ b/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). diff --git a/python_examples/GUIDE.md b/python_examples/GUIDE.md index 6745d04..28ce35f 100644 --- a/python_examples/GUIDE.md +++ b/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 ε = 5.276 kJ mol-1, +σ=0.483nm, +and two equal bonds of length σ. +The program employs the usual reduced units based on ε and σ +and in these units the potential cutoff of Mossa et al (2002) is _R_c=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_c. +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 ρ=0.32655 which is equivalent to ρ4=1.108g cm-3 +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 ε/kB≈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≈380K +at the specified densities ρ1, … ρ5 of Mossa et al (2002). +Here the excess pressure (_P_(ex)=_P_-ρ_T_ converted to MPa +with a factor 77.75 based on the values of ε and σ) +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 | ρ | _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 107 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_. diff --git a/python_examples/maths_module.py b/python_examples/maths_module.py index b0ad8d5..82853ed 100644 --- a/python_examples/maths_module.py +++ b/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,10 +148,9 @@ 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 @@ -152,7 +158,7 @@ def rotate_vector ( angle, axis, old ): 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 + diff --git a/python_examples/mc_nvt_poly_lj.py b/python_examples/mc_nvt_poly_lj.py new file mode 100755 index 0000000..c018cce --- /dev/null +++ b/python_examples/mc_nvt_poly_lj.py @@ -0,0 +1,203 @@ +#!/usr/bin/env python3 +# mc_nvt_poly_lj.py + +#------------------------------------------------------------------------------------------------# +# This software was written in 2016/17 # +# by Michael P. Allen / # +# and Dominic J. Tildesley ("the authors"), # +# to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), # +# published by Oxford University Press ("the publishers"). # +# # +# LICENCE # +# Creative Commons CC0 Public Domain Dedication. # +# To the extent possible under law, the authors have dedicated all copyright and related # +# and neighboring rights to this software to the PUBLIC domain worldwide. # +# This software is distributed without any warranty. # +# You should have received a copy of the CC0 Public Domain Dedication along with this software. # +# If not, see . # +# # +# DISCLAIMER # +# The authors and publishers make no warranties about the software, and disclaim liability # +# for all uses of the software, to the fullest extent permitted by applicable law. # +# The authors and publishers do not recommend use of this software for any purpose. # +# It is made freely available, solely to clarify points made in the text. When using or citing # +# the software, you should not imply endorsement by the authors or publishers. # +#------------------------------------------------------------------------------------------------# + +"""Monte Carlo, NVT ensemble, polyatomic molecule.""" + +def calc_variables ( ): + """Calculates all variables of interest. + + They are collected and returned as a list, for use in the main program. + """ + + # In this example we simulate using the shifted-force potential only + # The values of < p_sf >, < e_sf > and density should be consistent (for this potential) + # There are no long-range or delta corrections + + from averages_module import VariableType + + # Preliminary calculations + vol = box**3 # Volume + rho = n / vol # Density + + # Variables of interest, of class VariableType, containing three attributes: + # .val: the instantaneous value + # .nam: used for headings + # .method: indicating averaging method + # If not set below, .method adopts its default value of avg + # The .nam and some other attributes need only be defined once, at the start of the program, + # but for clarity and readability we assign all the values together below + + # Move acceptance ratio + m_r = VariableType ( nam = 'Move ratio', val = m_ratio, instant = False ) + + # Internal energy per molecule (shifted-force potential) + # Ideal gas contribution (assuming nonlinear molecules) plus total PE divided by N + e_sf = VariableType ( nam = 'E/N shifted force', val = 3.0*temperature + total.pot/n ) + + # Pressure (shifted-force potential) + # Ideal gas contribution plus total virial divided by V + p_sf = VariableType ( nam = 'P shifted force', val = rho*temperature + total.vir/vol ) + + # Collect together into a list for averaging + return [ m_r, e_sf, p_sf ] + +# Takes in a configuration of polyatomic molecules (positions and quaternions) +# Cubic periodic boundary conditions +# Conducts Monte Carlo at the given temperature +# Uses no special neighbour lists + +# Reads several variables and options from standard input using JSON format +# Leave input empty "{}" to accept supplied defaults + +# Positions r are divided by box length after reading in +# However, input configuration, output configuration, most calculations, and all results +# are given in simulation units defined by the model +# For example, for Lennard-Jones, sigma = 1, epsilon = 1 + +# Despite the program name, there is nothing here specific to Lennard-Jones +# The model (including the cutoff distance) is defined in mc_poly_lj_module + +import json +import sys +import numpy as np +import math +from config_io_module import read_cnf_mols, write_cnf_mols +from averages_module import run_begin, run_end, blk_begin, blk_end, blk_add +from maths_module import random_translate_vector, metropolis, random_rotate_quaternion, q_to_a +from mc_poly_lj_module import introduction, conclusion, potential, potential_1, PotentialType, na, db + +cnf_prefix = 'cnf.' +inp_tag = 'inp' +out_tag = 'out' +sav_tag = 'sav' + +print('mc_nvt_poly_lj') +print('Monte Carlo, constant-NVT ensemble, polyatomic molecule') +print('Simulation uses cut-and-shifted potential') + +# Read parameters in JSON format +try: + nml = json.load(sys.stdin) +except json.JSONDecodeError: + print('Exiting on Invalid JSON format') + sys.exit() + +# Set default values, check keys and typecheck values +defaults = {"nblock":10, "nstep":1000, "temperature":1.0, "dr_max":0.05, "de_max":0.05} +for key, val in nml.items(): + if key in defaults: + assert type(val) == type(defaults[key]), key+" has the wrong type" + else: + print('Warning', key, 'not in ',list(defaults.keys())) + +# Set parameters to input values or defaults +nblock = nml["nblock"] if "nblock" in nml else defaults["nblock"] +nstep = nml["nstep"] if "nstep" in nml else defaults["nstep"] +temperature = nml["temperature"] if "temperature" in nml else defaults["temperature"] +dr_max = nml["dr_max"] if "dr_max" in nml else defaults["dr_max"] +de_max = nml["de_max"] if "de_max" in nml else defaults["de_max"] + +introduction() +np.random.seed() + +# Write out parameters +print( "{:40}{:15d} ".format('Number of blocks', nblock) ) +print( "{:40}{:15d} ".format('Number of steps per block', nstep) ) +print( "{:40}{:15.6f}".format('Specified temperature', temperature) ) +print( "{:40}{:15.6f}".format('Maximum r displacement', dr_max) ) +print( "{:40}{:15.6f}".format('Maximum e displacement', de_max) ) + +# Read in initial configuration +n, box, r, e = read_cnf_mols ( cnf_prefix+inp_tag, quaternions=True ) +print( "{:40}{:15d} ".format('Number of particles', n) ) +print( "{:40}{:15.6f}".format('Box length', box) ) +print( "{:40}{:15.6f}".format('Density', n/box**3) ) +r = r / box # Convert positions to box units +r = r - np.rint ( r ) # Periodic boundaries + +# Calculate all bond vectors +d = np.empty ( (n,na,3), dtype=np.float_ ) +di = np.empty ( (na,3), dtype=np.float_ ) +for i, ei in enumerate(e): + ai = q_to_a ( ei ) # Rotation matrix for i + for a in range(na): + d[i,a,:] = np.dot ( db[a,:], ai ) # NB: equivalent to ai_T*db, ai_T=transpose of ai + +# Initial energy and overlap check +total = potential ( box, r, d ) +assert not total.ovr, 'Overlap in initial configuration' + +# Initialize arrays for averaging and write column headings +m_ratio = 0.0 +run_begin ( calc_variables() ) + +for blk in range(1,nblock+1): # Loop over blocks + + blk_begin() + + for stp in range(nstep): # Loop over steps + + moves = 0 + + for i in range(n): # Loop over atoms + rj = np.delete(r,i,0) # Array of all the other molecules + dj = np.delete(d,i,0) # Array of all the other molecules + partial_old = potential_1 ( r[i,:], d[i,:,:], box, rj, dj ) # Old potential, virial etc + assert not partial_old.ovr, 'Overlap in current configuration' + + ri = random_translate_vector ( dr_max/box, r[i,:] ) # Trial move to new position (in box=1 units) + ri = ri - np.rint ( ri ) # Periodic boundary correction + ei = random_rotate_quaternion ( de_max, e[i,:] ) # Trial rotation + ai = q_to_a ( ei ) # Rotation matrix for i + for a in range(na): + di[a,:] = np.dot ( db[a,:], ai ) # NB: equivalent to ai_T*db, ai_T=transpose of ai + partial_new = potential_1 ( ri, di, box, rj, dj ) # New atom potential, virial etc + + if not partial_new.ovr: # Test for non-overlapping configuration + delta = partial_new.pot - partial_old.pot # Use cut (but not shifted) potential + delta = delta / temperature + if metropolis ( delta ): # Accept Metropolis test + total = total + partial_new - partial_old # Update total values + r[i,:] = ri # Update position + e[i,:] = ei # Update quaternion + d[i,:,:] = di # Update bond vectors + moves = moves + 1 # Increment move counter + + m_ratio = moves / n + + blk_add ( calc_variables() ) + + blk_end(blk) # Output block averages + sav_tag = str(blk).zfill(3) if blk<1000 else 'sav' # Number configuration by block + write_cnf_mols ( cnf_prefix+sav_tag, n, box, r*box, e ) # Save configuration + +run_end ( calc_variables() ) + +total = potential ( box, r, d ) # Double check book-keeping +assert not total.ovr, 'Overlap in final configuration' + +write_cnf_mols ( cnf_prefix+out_tag, n, box, r*box, e ) # Save configuration +conclusion() diff --git a/python_examples/mc_poly_lj_module.py b/python_examples/mc_poly_lj_module.py new file mode 100644 index 0000000..1abdb4e --- /dev/null +++ b/python_examples/mc_poly_lj_module.py @@ -0,0 +1,215 @@ +#!/usr/bin/env python3 +# mc_poly_lj_module.py + +#------------------------------------------------------------------------------------------------# +# This software was written in 2016/17 # +# by Michael P. Allen / # +# and Dominic J. Tildesley ("the authors"), # +# to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), # +# published by Oxford University Press ("the publishers"). # +# # +# LICENCE # +# Creative Commons CC0 Public Domain Dedication. # +# To the extent possible under law, the authors have dedicated all copyright and related # +# and neighboring rights to this software to the PUBLIC domain worldwide. # +# This software is distributed without any warranty. # +# You should have received a copy of the CC0 Public Domain Dedication along with this software. # +# If not, see . # +# # +# DISCLAIMER # +# The authors and publishers make no warranties about the software, and disclaim liability # +# for all uses of the software, to the fullest extent permitted by applicable law. # +# The authors and publishers do not recommend use of this software for any purpose. # +# It is made freely available, solely to clarify points made in the text. When using or citing # +# the software, you should not imply endorsement by the authors or publishers. # +#------------------------------------------------------------------------------------------------# + +"""Routines for MC simulation, polyatomic molecule, LJ atoms.""" + +import numpy as np + +fast = True # Change this to replace NumPy potential evaluation with slower Python + +# Bond vectors in body-fixed frame +# Isosceles triangle, 3 sites, with unit bond length and bond angle alpha +# which we set to 75 degrees here +alpha = 75.0 * np.pi / 180.0 +alpha2 = alpha / 2.0 +na = 3 +db = np.array([[-np.sin(alpha2), 0.0, -np.cos(alpha2)/3.0], + [0.0, 0.0, 2.0*np.cos(alpha2)/3.0], + [np.sin(alpha2), 0.0, -np.cos(alpha2)/3.0]], dtype=np.float_) +diameter = 2.0 * np.sqrt ( np.max ( np.sum(db**2,axis=1) ) ) # Molecular diameter + +# Cutoff distance and force-shift parameters (all private) chosen as per the reference: +# S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, Phys Rev E, 65, 041205 (2002) +r_cut = 2.612 # in sigma=1 units, where r_cut = 1.2616 nm, sigma = 0.483 nm +sr_cut = 1.0/r_cut +sr_cut6 = sr_cut**6 +sr_cut12 = sr_cut6**2 +lambda1 = 4.0*(7.0*sr_cut6-13.0*sr_cut12) +lambda2 = -24.0*(sr_cut6-2.0*sr_cut12)*sr_cut + +class PotentialType: + """A composite variable for interactions.""" + + def __init__(self, pot, vir, ovr): + self.pot = pot # the potential energy cut (but not shifted) at r_cut + self.vir = vir # the virial + self.ovr = ovr # a flag indicating overlap (i.e. pot too high to use) + + def __add__(self, other): + pot = self.pot + other.pot + vir = self.vir + other.vir + ovr = self.ovr or other.ovr + return PotentialType(pot,vir,ovr) + + def __sub__(self, other): + pot = self.pot - other.pot + vir = self.vir - other.vir + ovr = self.ovr or other.ovr # This is meaningless, but inconsequential + return PotentialType(pot,vir,ovr) + +def introduction(): + """Prints out introductory statements at start of run.""" + + print('Lennard-Jones potential') + print('Cut-and-force-shifted') + print('Diameter, sigma = 1') + print('Well depth, epsilon = 1') + if fast: + print('Fast NumPy potential routine') + else: + print('Slow Python potential routine') + + print( "{:40}{:15d}".format('Number of atoms per molecule', na) ) + for i, b in enumerate(db): + print( "{}{:2d}{:15.6f}{:15.6f}{:15.6f}".format('Body-fixed atom vector',i,*b)) + print( "{:40}{:15.6f}".format('Molecular diameter', diameter) ) + + print( "{:40}{:15.6f}".format('r_cut', r_cut) ) + print( "{:40}{:15.6f}".format('Force-shift lambda1', lambda1) ) + print( "{:40}{:15.6f}".format('Force-shift lambda2', lambda2) ) + +def conclusion(): + """Prints out concluding statements at end of run.""" + + print('Program ends') + +def potential ( box, r, d ): + """Takes in box and r & d arrays, and calculates total potential etc. + + The results are returned as total, a PotentialType variable. + """ + # Actual calculation performed by function potential_1 + + n, ndim = r.shape + assert ndim==3, 'Dimension error for r' + nn, nna, ndim = d.shape + assert nna==na and ndim==3, 'Dimension error for d' + assert n==nn, 'Dimension mismatch for r and d' + + total = PotentialType ( pot=0.0, vir=0.0, ovr=False ) + + for i in range(n-1): + partial = potential_1 ( r[i,:], d[i,:,:], box, r[i+1:,:], d[i+1:,:,:] ) + if partial.ovr: + total.ovr = True + break + total = total + partial + + return total + +def potential_1 ( ri, di, box, r, d ): + """Takes in r & d of a molecule and calculates its interactions. + + Values of box and partner coordinate arrays are supplied. + The results are returned as partial, a PotentialType variable. + """ + + import numpy as np + + # partial.pot is the nonbonded cut (not shifted) potential energy of atom ri with a set of other atoms + # partial.vir is the corresponding virial of atom ri + # partial.ovr is a flag indicating overlap (potential too high) to avoid overflow + # If this is True, the values of partial.pot etc should not be used + # In general, r & d will be subsets of the complete set of simulation coordinates & bond vectors + # and none of their rows should be identical to ri, di + + # It is assumed that r has been divided by box + # Results are in LJ units where sigma = 1, epsilon = 1 + # Note that this is the force-shifted LJ potential with a linear smoothing term + # S Mossa, E La Nave, HE Stanley, C Donati, F Sciortino, P Tartaglia, Phys Rev E, 65, 041205 (2002) + + nj, ndim = r.shape + assert ndim==3, 'Dimension error for r' + nnj, nna, ndim = d.shape + assert nna==na and ndim==3, 'Dimension error for d' + assert nj==nnj, 'Dimension mismatch between r and d' + assert ri.size==3, 'Dimension error for ri' + nna, ndim = di.shape + assert nna==na and ndim==3, 'Dimension error for di' + + sr2_ovr = 1.77 # Overlap threshold (pot > 100) + rm_cut_box = ( r_cut + diameter ) / box # Molecular cutoff in box=1 units + rm_cut_box_sq = rm_cut_box**2 # squared + assert rm_cut_box<0.5, 'rm_cut/box too large' + r_cut_sq = r_cut ** 2 + + partial = PotentialType ( pot=0.0, vir=0.0, ovr=False ) # Initialize + + if fast: + rij = ri - r # Get all separation vectors from partners + rij = rij - np.rint(rij) # Periodic boundary conditions in box=1 units + rij = rij * box # Now in sigma=1 units + for a in range(na): + for b in range(na): + rab = rij + di[a,:] - d[:,b,:] # All atom-atom vectors for given a and b + rab_sq = np.sum(rab**2,axis=1) # Squared separations + in_range = rab_sq < r_cut_sq # Set flags for within cutoff + sr2 = 1.0 / rab_sq # (sigma/rab)**2 + ovr = sr2 > sr2_ovr # Set flags for any overlaps + if np.any(ovr): + partial.ovr = True + return partial + rmag = np.sqrt(rab_sq) + sr6 = sr2 ** 3 + sr12 = sr6 ** 2 + pot = np.where ( in_range, + 4.0*(sr12-sr6) + lambda1 + lambda2*rmag, 0.0 ) # LJ atom-atom pair potentials (force-shifted) + virab = np.where ( in_range, + 24.0*(2.0*sr12-sr6) - lambda2*rmag, 0.0 ) # LJ atom-atom pair virials + fab = virab * sr2 + fab = rab * fab[:,np.newaxis] # LJ atom-atom pair forces + partial = partial + PotentialType ( pot=np.sum(pot), vir=np.sum(rij*fab), ovr=False ) + + else: + for j, rj in enumerate(r): + rij = ri - rj # Separation vector + rij = rij - np.rint(rij) # Periodic boundary conditions in box=1 units + rij_sq = np.sum(rij**2) # Squared separation + if rij_sq < rm_cut_box_sq: # Check within cutoff + rij = rij * box # Now in sigma=1 units + for a in range(na): + for b in range(na): + rab = rij + di[a,:] - d[j,b,:] # Atom-atom vector, sigma=1 units + rab_sq = np.sum ( rab**2 ) # Squared atom-atom separation, sigma=1 units + + if rab_sq < r_cut_sq: # Test within potential cutoff + sr2 = 1.0 / rab_sq # (sigma/rab)**2 + ovr = sr2 > sr2_ovr # Overlap if too close + if ovr: + partial.ovr=True + return partial + rmag = np.sqrt(rab_sq) + sr6 = sr2 ** 3 + sr12 = sr6 ** 2 + pot = 4.0*(sr12-sr6) + lambda1 + lambda2*rmag # LJ atom-atom pair potential (force-shifted) + virab = 24.0*(2.0*sr12-sr6) - lambda2*rmag # LJ atom-atom pair virial + fab = rab * virab * sr2 # LJ atom-atom pair force + partial = partial + PotentialType ( pot=pot, vir=np.sum(rij*fab), ovr=ovr ) + + # Include numerical factors + partial.vir = partial.vir / 3.0 # Divide virial by 3 + + return partial