Analysis

Structure factors: S(q), X-ray and neutron

At low q the exact sum over every reciprocal-lattice vector the cell allows; above the switch-over, the Fourier transform of the partial pair distributions with a Lorch window against truncation ripples. X-ray and neutron patterns weight the partials by form factors or scattering lengths (Faber–Ziman).

Faber & Ziman 1965; Lorch 1969

When to use it

Comparing a cell with diffraction data, and neutron contrast with deuterated backbones or rings.

Source

core/src/properties.cpp

Tested by

Scattering.DeuterationPatternsSplitTheHydrogens · bench/analyze/check_analyze.py (MDAnalysis)

Departure from the reference

no q below the cell's first shell 2π/L; no instrument resolution or polarisation factors

References

1. Faber, T. E., Ziman, J. M., "A theory of the electrical properties of liquid metals III. The resistivity of binary alloys", Phil. Mag. 11, 153–173 (1965). doi:10.1080/14786436508211931
2. Lorch, E., "Neutron diffraction by germania, silica and radiation-damaged silica glasses", J. Phys. C 2, 229–237 (1969). doi:10.1088/0022-3719/2/2/305
3. Sears, V. F., "Neutron scattering lengths and cross sections", Neutron News 3, 26–37 (1992). doi:10.1080/10448639208218770
4. Prince, E. (ed.), "International Tables for Crystallography, Volume C: Mathematical, Physical and Chemical Tables", Kluwer (2004)

Free volume and pore sizes

Probe insertion on a grid: a point is free when a probe of the given radius centred there overlaps no van der Waals sphere. Bondi's fractional free volume uses the occupied van der Waals volume; each free point's pore size is the largest atom-free sphere containing it.

Bondi 1964; Gelb & Gubbins 1999

When to use it

Gas permeation, voids in glassy polymers and filled rubbers, pores of frameworks.

Source

core/src/properties.cpp · core/src/voids.cpp

Tested by

Voids.OneAtomLeavesTheFarCornerEmpty · bench/analyze/check_analyze.py

Departure from the reference

the 1.3 packing factor of Bondi is used as published

References

1. Bondi, A., "van der Waals volumes and radii", J. Phys. Chem. 68, 441–451 (1964). doi:10.1021/j100785a001
2. Gelb, L. D., Gubbins, K. E., "Pore size distributions in porous glasses: a computer simulation study", Langmuir 15, 305–308 (1999). doi:10.1021/la9808418

Mean-square displacement and diffusion

Displacements of atoms and of molecule centres over every time origin, with the drift of the whole system removed. The diffusion coefficient comes from the linear part of the molecule-centre MSD.

Einstein 1905

When to use it

Small molecules in polymers (curatives, solvents, gases). Check the log–log slope reaches 1 before trusting D.

Source

core/src/properties.cpp

Tested by

bench/analyze/check_analyze.py (MDAnalysis)

Departure from the reference

no finite-size (Yeh–Hummer) correction

References

1. Einstein, A., "Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen", Ann. Phys. 322, 549–560 (1905). doi:10.1002/andp.19053220806

Adsorption locator

Low-energy places of adsorbate molecules on a fixed substrate (a crystal surface, a filler, a polymer surface or cell) by Monte Carlo simulated annealing, as Materials Studio's Adsorption Locator does. The adsorbates, built from SMILES and added after the substrate, move as rigid bodies: translations, rotations about their centre and jumps to a random place and orientation, accepted by the Metropolis rule at a temperature falling geometrically from the start to the end value within each cycle; jumps fade out as it cools, and each cycle ends with a downhill polish of its lowest configuration with shrinking steps; later cycles restart from the best so far. The energy is the assigned force field's non-bonded interaction between molecules.

Kirkpatrick, Gelatt & Vecchi 1983; Metropolis et al. 1953

When to use it

Binding sites and energies of small molecules on surfaces and fillers: coupling agents or sizing molecules on silica, carbon or cellulose; gases, water and plasticisers on or in a polymer; comparing adsorbates by E_ads. Rigid molecules: no deformation energy; relax the result with Minimise for that.

Source

core/src/adsorption.cpp (locate_adsorption), core/src/pairmodel.cpp

Tested by

Adsorption.FindsTheScanMinimum (the lowest energy of an atom over a lattice equals a fine grid scan's, −0.8980 kcal/mol at the hollow site); Adsorption.RigidMoleculesAndDeDn

Departure from the reference

damped shifted force electrostatics, not Ewald; rigid adsorbates (one conformer each)

References

1. Kirkpatrick, S., Gelatt, C. D., Vecchi, M. P., "Optimization by simulated annealing", Science 220, 671–680 (1983). doi:10.1126/science.220.4598.671

Sorption: Widom insertion and grand-canonical Monte Carlo

How much of a gas or small molecule a fixed host takes up — a polymer or rubber cell, a filler, a porous crystal — as Materials Studio's Sorption computes it. The sorbate is a rigid molecule typed with the host's force field. Widom test-particle insertion averages the Boltzmann factor of one molecule placed at random places and orientations: the excess chemical potential, the Henry constant and the solubility coefficient follow. Grand-canonical Monte Carlo at each pressure (insertions, deletions, translations and rotations; an ideal-gas reservoir, fugacity = pressure) gives the loading, its error from ten blocks, and the isosteric heat from the fluctuations of energy and number. A gas mixture is several rigid sorbates at once, each at its share of the pressure (f_i = y_i p): insertions and deletions pick a gas at random, and each gas gets its loading, its Widom terms alone, its isosteric heat from the multicomponent fluctuations, and the adsorption selectivity over the first gas. A rigid model (TraPPE's CO₂, N₂, O₂: bond and angle constants 0) takes its force field's own bond lengths and angles before the run (TraPPE CO₂: C=O 1.16 Å, linear); with a united-atom host the sorbate keeps its types' own charges (TraPPE CO₂ +0.70/−0.35 e). A sorbate gets its force field's own preparation (united-atom sites; a model's massless charge sites such as TraPPE N₂'s centre, added as virtual sites).

Widom 1963; Frenkel & Smit 2002

When to use it

Gas solubility in rubbers and polymers (CO₂, CH₄, N₂, O₂, water) for barrier and permeability estimates (P = D·S with D from Diffusion), uptake of fillers and porous crystals. The host is rigid: swelling and the host's relaxation around the gas are not included; equilibrated, full-density cells are needed for polymer solubilities. Mixtures: CO₂/N₂ or CO₂/CH₄ selectivity of a membrane polymer, water against a gas in a rubber, competing curative fragments.

Source

core/src/sorption.cpp (sorption)

Tested by

Sorption.WidomEqualsTheGridIntegral (⟨e^(−βU)⟩ 0.03985 ± 0.00055 against the exact grid average 0.03995); Sorption.GcmcIdealGasAndHenryLimit (⟨N⟩ = βpV for an ideal gas; the Henry regime equals Widom's); Sorption.BinaryMixtureIdealAndHenry (an ideal binary gas takes y_i βpV of each, 6.02 + 13.91 for 6 + 14; in a host at low loading each gas follows its own Henry law and the selectivity 0.886 equals the Widom ratio 0.879); Sorption.TrappeCo2VirialAndGcmc (TraPPE CO₂: B₂(300 K) −111 cm³/mol against −121 measured; GCMC of the pure gas at 10 bar 56.6 ± 0.5 molecules against 57.1 from the model's virial expansion); bench/ff/check_trappe_co2.py (CO₂ σ, ε, charges as gmso's TraPPE file; 100 CO₂ CAPS against LAMMPS −31.3923 / −31.3921 kcal/mol); Sorption.TrappeN2AgainstNistBenchmark (TraPPE N₂ against NIST's benchmark simulations of the model, 110 K, 1.0 mol/L, Z = 0.8780 ± 0.0007: B₂ −125.7 against B₂ + B₃ρ = −122.0 cm³/mol; GCMC at that state's fugacity 0.984 ± 0.010 mol/L); TraPPE CO₂ against NIST's benchmark simulations of the model at 300 K: B₂ −111.3 against (Z − 1)/ρ = −108.5 cm³/mol at 0.1 mol/L, GCMC at the fugacity of the 0.5 mol/L state 0.493 ± 0.005 mol/L

Departure from the reference

ideal-gas fugacity (no equation of state); rigid host and sorbate; damped shifted force electrostatics

References

1. Frenkel, D., Smit, B., "Understanding Molecular Simulation: From Algorithms to Applications", Academic Press (2002)
2. Potoff, J. J., Siepmann, J. I., "Vapor-liquid equilibria of mixtures containing alkanes, carbon dioxide, and nitrogen", AIChE J. 47, 1676–1682 (2001). doi:10.1002/aic.690470719
3. Widom, B., "Some topics in the theory of fluids", J. Chem. Phys. 39, 2808–2812 (1963). doi:10.1063/1.1734110

Shear viscosity by Green–Kubo

The shear viscosity is the time integral of the pressure tensor's autocorrelation in equilibrium. CAPS runs NVT dynamics of a copy of the frame (Nosé–Hoover, weakly coupled so the dynamics are nearly Newtonian), samples the full pressure tensor (kinetic and virial) every few femtoseconds, and averages the autocorrelation of all five independent components of its traceless symmetric part. The running integral should level off; its mean over the last third of the window is the value, and blocks of the run give the error. The result feeds the Yeh–Hummer finite-size correction of diffusion coefficients.

Daivis & Evans 1994

When to use it

Solvents and small molecules converge in 100–200 ps; polymer melts relax over nanoseconds to microseconds, and the report says when the integral has not levelled off (then the value is at best a lower bound).

Source

core/src/mechanics.cpp (green_kubo_viscosity, viscosity_green_kubo)

Tested by

Mechanics.GreenKuboViscosity (an Ornstein–Uhlenbeck stress with a known integral, within 6 %)

Departure from the reference

one long run split in blocks rather than independent replicas; no time decomposition of the integral

References

1. Daivis, P. J., Evans, D. J., "Comparison of constant pressure and constant volume nonequilibrium simulations of sheared model decane", J. Chem. Phys. 100, 541–547 (1994). doi:10.1063/1.466970
2. Martyna, G. J., Klein, M. L., Tuckerman, M., "Nosé–Hoover chains: the canonical ensemble via continuous dynamics", J. Chem. Phys. 97, 2635–2643 (1992). doi:10.1063/1.463940

Elastic constants by static strain

The minimised cell is strained by ±ε along each of the six components, re-minimised with the pair list frozen, and the stress read from the virial; C_ij are the central differences. Isotropic moduli are Hill averages of the Voigt and Reuss bounds.

Theodorou & Suter 1986

When to use it

Glassy polymers and crystals at 0 K. Rubbers above Tg need the fluctuation method or tensile runs.

Source

core/src/mechanics.cpp

Tested by

Mechanics.StaticConstantsOfBravaisLattice · bench/mechanics/check_born_lammps.py

Departure from the reference

none

References

1. Theodorou, D. N., Suter, U. W., "Atomistic modeling of mechanical properties of polymeric glasses", Macromolecules 19, 139–154 (1986). doi:10.1021/ma00155a022

Elastic constants from stress fluctuations

At finite temperature the elastic constants are the Born term (second strain derivatives of the energy) minus the stress fluctuations, plus the kinetic term, averaged over an NVT run.

Lutsko 1989; Clavier et al. 2017

When to use it

Finite-temperature constants. For bonded polymers the fluctuation term converges slowly: errors stay at the GPa level over 100 ps, which CAPS reports.

Source

core/src/mechanics.cpp

Tested by

Mechanics.FluctuationConstantsOfColdCrystal · Mechanics.BornMatrixMatchesLatticeSum

Departure from the reference

none

References

1. Lutsko, J. F., "Generalized expressions for the calculation of elastic constants by computer simulation", J. Appl. Phys. 65, 2991–2997 (1989). doi:10.1063/1.342716
2. Clavier, G., Desbiens, N., Bourasseau, E., Lachet, V., Brusselle-Dupend, N., Rousseau, B., "Computation of elastic constants of solids using molecular simulation: comparison of constant volume and constant pressure ensemble methods", Mol. Simul. 43, 1413–1422 (2017). doi:10.1080/08927022.2017.1313418

χ from pair contacts

The Flory–Huggins χ(T) of two components from the energies of their molecules in contact. Each component is one rigid molecule (a solvent, or a repeat unit with its ends capped by hydrogen). Random relative orientations and directions are drawn; the second molecule slides in until the van der Waals surfaces touch, and the pair energy is recorded. The energies are Boltzmann-averaged at each temperature, the coordination numbers come from packing as many molecules as fit around one, and χ(T) is fitted to A + B/T.

Fan, Olafson, Blanco & Hsu 1992

When to use it

A predictive screen of solvents and blend partners when no measured solubility parameters exist. In CAPS's check (bench/chi/check_contacts.py) it agrees with known behaviour for 11 of 16 polymer–solvent pairs (Hildebrand: 12 of 16) and fails where single strong contacts dominate: THF, toluene with natural rubber, and polystyrene/PVME (called immiscible; it is miscible). Self-mixing controls resolve χ to about ±0.15.

Source

core/src/chipair.cpp

Tested by

ChiPair.ArgonPairsTouchAtTheBondiSumAndPackUnderTwelve · ChiPair.SelfMixingControlIsZeroWithinItsErrorAndWaterIsNot · bench/chi/check_contacts.py

Departure from the reference

rigid molecules with one conformer each; one molecule per lattice site of either kind, with no correction for their different sizes; coordination numbers by random sequential packing

References

1. Fan, C. F., Olafson, B. D., Blanco, M., Hsu, S. L., "Application of molecular simulation to derive phase diagrams of binary mixtures", Macromolecules 25, 3667–3676 (1992)
2. Blanco, M., "Molecular silverware. I. General solutions to excluded volume constrained problems", J. Comput. Chem. 12, 237–247 (1991)

Polyhedral template matching

Each particle's nearest neighbours, with the particle itself, are compared with ideal FCC, HCP, BCC, icosahedral and simple-cubic neighbourhoods. Both point sets are centred and scaled to unit mean distance; neighbours are matched to template points from anchor pairs whose angles agree, and the optimal rotation (Horn's quaternion) gives the RMSD. The template with the least RMSD below the cutoff is the particle's structure; its rotation is the lattice orientation, and the deformation left after it gives a shear strain.

Larsen, Schmidt & Schiøtz 2016

When to use it

Crystalline fillers and metal particles in a composite: grains, surfaces, stacking faults and local strain; the RMSD table shows where to set the cutoff.

Source

core/src/pipeline.cpp (step_ptm) · core/src/superpose.cpp

Tested by

Structure.PolyhedralTemplateMatching

Departure from the reference

correspondences from anchor pairs and nearest-point assignment, not the paper's graph-canonical search; the shear strain is relative to the particle's own lattice scale

References

1. Larsen, P. M., Schmidt, S., Schiøtz, J., "Robust structural identification via polyhedral template matching", Modelling Simul. Mater. Sci. Eng. 24, 055007 (2016). doi:10.1088/0965-0393/24/5/055007

Voronoi and radical tessellation

Each atom's cell is built on its own: a box around it is clipped by the bisecting plane of every neighbour, nearest first, until no farther neighbour can reach the cell (periodic images included). The radical (power) variant moves each plane so that atoms with larger van der Waals radii get larger cells. The cells tile the box exactly; their volumes, faces and the Voronoi index ⟨n3 n4 n5 n6⟩ (faces with 3, 4, 5, 6 edges) follow.

cell-by-cell construction as in Rycroft 2009; radical planes Gellatly & Finney 1982

When to use it

Free volume per atom in a polymer (radical), local packing in amorphous metals and fillers; the grid methods remain for a quick look.

Source

core/src/voronoi.cpp · core/src/pipeline.cpp (step_voronoi)

Tested by

Structure.ExactVoronoiCells

Departure from the reference

none known; cells of atoms without a periodic cell around them may be open (reported)

References

1. Rycroft, C. H., "VORO++: A three-dimensional Voronoi cell library in C++", Chaos 19, 041111 (2009). doi:10.1063/1.3215722
2. Gellatly, B. J., Finney, J. L., "Characterisation of models of multicomponent amorphous metals: the radical alternative to the Voronoi polyhedron", J. Non-Cryst. Solids 50, 313–329 (1982). doi:10.1016/0022-3093(82)90093-X

Wigner–Seitz defect analysis

Every particle is assigned to the reference site nearest to it — the site whose Wigner–Seitz cell holds it. A site without a particle is a vacancy; each particle beyond the first on a site is an interstitial. The reference is a frame of the trajectory (the undamaged lattice) or a file.

When to use it

Damage in crystalline fillers or substrates: after deformation, irradiation or cutting.

Source

core/src/pipeline.cpp (step_wigner_seitz)

Tested by

Structure.WignerSeitzDefects

Departure from the reference

none known

Backbone conformation: torsions, angles, bonds

The dihedral angle of every four consecutive backbone atoms of every chain and frame, as a distribution with the trans and gauche fractions; the backbone bond angles and lengths with it.

IUPAC sign convention (cis 0°, trans ±180°); the trans/gauche boundaries of the rotational isomeric state model

When to use it

Chain stiffness and packing: a melt's t/g ratio against the isolated chain's, the effect of a filler surface or of strain on conformations.

Source

core/src/properties.cpp

Tested by

tests/test_core.cpp Analyze.BackboneConformation (an all-trans zigzag and one gauche torsion)

Departure from the reference

every backbone torsion counts alike (ester and amide torsions with C–C ones)

Velocity autocorrelation and vibrational density of states

The velocity autocorrelation over every frame as a time origin; its integral gives the Green–Kubo diffusion coefficient and its cosine transform (mass-weighted, Hann window) the vibrational density of states.

Dickey & Paskin 1969; the spectrum reaches the Nyquist wavenumber 1/(2cΔt)

When to use it

Vibrational spectra of a cell (needs frames a few fs apart: C–H stretches at 3000 cm⁻¹ need ≤ 5 fs); D of small molecules as a check on the MSD.

Source

core/src/properties.cpp

Tested by

tests/test_core.cpp Analyze.VacfSpectrumPeak (a 1000 cm⁻¹ oscillation peaks at 1000 cm⁻¹)

Departure from the reference

no quantum correction of the classical spectrum

References

1. Dickey, J. M., Paskin, A., "Computer simulation of the lattice dynamics of solids", Phys. Rev. 188, 1407–1418 (1969). doi:10.1103/PhysRev.188.1407

Hydrogen bonds

Donor–hydrogen···acceptor contacts between N, O and F by distance and angle, counted per frame; the lifetime from the intermittent correlation of each bond's existence.

Luzar & Chandler 1996

When to use it

Epoxidised and hydroxylated rubbers, maleic acid and PBS ends, water: how many hydrogen bonds hold the network and how long they last.

Source

core/src/properties.cpp

Tested by

tests/test_core.cpp Analyze.HydrogenBonds (two waters, one bond; turned away, none)

Departure from the reference

the lifetime is a lower bound when C(t) has not decayed over the trajectory

References

1. Luzar, A., Chandler, D., "Hydrogen-bond kinetics in liquid water", Nature 379, 55–57 (1996). doi:10.1038/379055a0

van Hove self-correlation and the non-Gaussian parameter

How far atoms move in a time t: the distribution of their displacements over every time origin, and how far it is from a Gaussian.

α₂ = 0 for Gaussian (free, Fickian) displacements; a maximum marks heterogeneous dynamics near the glass transition

When to use it

Polymer melts and rubbers near Tg; diffusing small molecules in a matrix (hopping shows as a second peak of 4πr²G_s).

Source

core/src/properties.cpp

Tested by

tests/test_core.cpp Analyze.VanHoveAndP2 (Gaussian steps: α₂ ≈ 0; half the atoms frozen: α₂ > 0.4)

Departure from the reference

Orientational correlation P₂(r)

Whether backbone segments near each other point the same way: the second Legendre polynomial of the angle between two chords, against the distance between their centres.

1 parallel, 0 random, −½ perpendicular

When to use it

Local alignment in a stretched rubber, chain order next to a fibre or filler surface, the onset of crystallisation (PBS, PE).

Source

core/src/properties.cpp

Tested by

tests/test_core.cpp Analyze.VanHoveAndP2 (parallel all-trans chains: P₂ > 0.95)

Departure from the reference

at most 4000 chords a frame (every n-th)

Static dielectric constant from dipole fluctuations

The fluctuations of the cell's total dipole moment over an equilibrium trajectory give its static dielectric constant (conducting surroundings, as Ewald sums with tin-foil boundaries assume).

Neumann 1983; molecules made whole across the cell walls before M is summed

When to use it

Insulating rubbers and fillers for cables, dielectric elastomers: ε of the matrix and of composites (with long NVT or NPT runs; the running ε shows convergence).

Source

core/src/properties.cpp

Tested by

tests/test_core.cpp Analyze.DielectricFromDipoleFluctuations (a flipping ±1 e pair in a 10 Å cell at 300 K: 3.332)

Departure from the reference

no electronic polarisability: ε∞ = 1 (add it for a comparison with experiment)

References

1. Neumann, M., "Dipole moment fluctuation formulas in computer simulations of polar systems", Molecular Physics 50, 841–858 (1983). doi:10.1080/00268978300102721

Heat capacity, compressibility and expansion from fluctuations

An equilibrium NPT run gives Cp, the isothermal compressibility κ_T (and the bulk modulus B = 1/κ_T) and the thermal expansion α_P from the fluctuations of the enthalpy and the volume; an NVT run gives Cv. One run, no finite differences.

the frames carry no velocities: the kinetic part is exact in classical statistics, ⟨δK²⟩ = 3N(k_BT)²/2, uncorrelated with the configuration

When to use it

Rubber compounds and glassy polymers: κ_T and B for the bulk response, α_P below and above Tg (compare with the slope of V(T) from the Tg scan), Cp for thermal modelling.

Source

core/src/properties.cpp

Tested by

tests/test_core.cpp Analyze.ResponseFunctionsFromFluctuations (argon beyond the cut-off: U = 0, the volume's variance alone; NVT gives 3R/2M)

Departure from the reference

classical heat capacity: every vibration holds k_BT, so Cp exceeds experiment where high-frequency modes are frozen out; fluctuations converge slowly (the ± is the spread of the blocks)

References

1. Allen, M. P., Tildesley, D. J., "Computer Simulation of Liquids", Oxford University Press (2017). doi:10.1093/oso/9780198803195.001.0001

Compliance, directional moduli and sound speeds

From the elastic constants: the compliance S = C⁻¹, Young's moduli and Poisson's ratios along the cell axes, the shear moduli, how far from isotropic the solid is, and its longitudinal, transverse and Debye-mean sound speeds.

K and G are the Hill averages; ν_ij is the strain along j over the strain along i when stretched along i; v_m = [(2/v_T³ + 1/v_L³)/3]^(−1/3)

When to use it

Oriented or filled rubbers and fibre composites: E along and across the fibres, the anisotropy of a stretched network; sound speeds for acoustic and damping comparisons (ultrasonic measurements).

Source

core/src/mech_props.cpp

Tested by

tests/test_mechanics.cpp Mechanics.ComplianceAndSoundSpeeds (an isotropic tensor: E_x, ν_xy, G, A^U = 0 and both sound speeds as the formulas)

Departure from the reference

the speeds are of the isotropic average; along a direction of an anisotropic solid they follow from the Christoffel equation

References

1. Ranganathan, S. I., Ostoja-Starzewski, M., "Universal elastic anisotropy index", Physical Review Letters 101, 055504 (2008). doi:10.1103/PhysRevLett.101.055504
2. Anderson, O. L., "A simplified method for calculating the Debye temperature from elastic constants", Journal of Physics and Chemistry of Solids 24, 909–917 (1963). doi:10.1016/0022-3697(63)90067-2

Normal modes

The harmonic vibrations of a structure at a minimum of its force field: wavenumbers, the IR spectrum of the fixed charges, the vibrational density of states and the quantum harmonic zero-point energy, entropy and heat capacity. A mode plays in the view from its card.

H by central differences of the forces (±0.005 Å, the pair list frozen); the rigid translations and rotations projected out (Miller, Handy & Adams 1980); λ < 0 shown as a negative (imaginary) wavenumber

When to use it

Checking a force field's bonded terms against IR/Raman bands (C=C of a diene, C=O of an ester, S–S of a crosslink), the vibrational entropy of a conformer, and whether a minimised structure is really a minimum (no imaginary modes).

Source

core/src/normal_modes.cpp

Tested by

tests/test_normal_modes.cpp (the eigen-solver on a random matrix; a harmonic diatomic gives √(2k/μ)/2πc with five rigid motions projected; water has three real modes)

Departure from the reference

harmonic, fixed charges: no anharmonicity, no polarisation in the intensities; ZPE and S from the force field's frequencies

References

1. Miller, W. H., Handy, N. C., Adams, J. E., "Reaction path Hamiltonian for polyatomic molecules", The Journal of Chemical Physics 72, 99–112 (1980). doi:10.1063/1.438959

Shear viscosity η(γ̇) by NEMD

Planar Couette flow imposed on the structure: the x velocity grows along y at the shear rate, the cell tilts with the flow, and the steady-state shear stress gives the viscosity at that rate. A sweep of rates shows shear thinning and the first normal-stress difference of polymer melts.

SLLOD (Evans & Morriss 1984) with Lees–Edwards boundaries by the cell's tilt (b_x grows at γ̇ b_y, flipped by a at ±a_x/2); the thermostat acts on the peculiar momenta p; the first quarter of each run left out

When to use it

Melt viscosity and shear thinning of rubber compounds and oils at processing-like (but far higher) rates; the power-law index n for flow models; N₁ > 0 shows the chains orienting in the flow. For long runs on a cluster, Export writes the same flow for LAMMPS (protocol Shear viscosity, NEMD): the box made triclinic, fix nvt/sllod with fix deform xy erate remap v, and η = −P_xy/γ̇ in mPa·s averaged in blocks into viscosity.dat; one file per rate.

Source

core/src/mechanics.cpp

Tested by

tests/test_dynamics.cpp Dynamics.SllodShearViscosityOfTheLjLiquid and tests/test_mechanics.cpp Mechanics.NemdShearThinningOfTheLjLiquid (the Lennard-Jones triple point: η* ≈ 2.1 at γ̇* = 1, lower at 2, ⟨T⟩ held)

Departure from the reference

simulated rates are 10⁹–10¹¹ s⁻¹; viscous heat leaves T above the target by about the heating rate × τ (τ 20 fs by default; the mean T is reported per rate)

References

1. Evans, D. J., Morriss, G. P., "Nonlinear-response theory for steady planar Couette flow", Physical Review A 30, 1528–1530 (1984). doi:10.1103/PhysRevA.30.1528

Conformer search

The low-energy shapes of a molecule under its force field: many starts minimised — random staggered torsions, or snapshots of a 1000 K run quenched — and the minima grouped when their heavy atoms superpose, with energies and Boltzmann populations. They open as frames to step through.

rotatable bonds: acyclic single bonds with a heavy atom beyond each end, set to 60°, 180° or 300° (± 15° noise); superposition by Horn's quaternion method

When to use it

Crosslinkers, curatives, plasticisers and monomers before packing or reacting: the lowest conformer as the building block, the spread of shapes, the anti/gauche balance of a chain segment; the anneal method also shows the distribution of minima.

Source

core/src/conformers.cpp

Tested by

tests/test_conformers.cpp (butane: anti lowest, gauche a fraction of a kcal/mol above, populations sum to 1; benzene one conformer; the anneal search minimises every snapshot)

Departure from the reference

vacuum, energies only (no vibrational entropy, no solvent); rings keep their shape; mirror images are separate only when the heavy atoms do not superpose

References

1. Horn, B. K. P., "Closed-form solution of absolute orientation using unit quaternions", Journal of the Optical Society of America A 4, 629–642 (1987). doi:10.1364/JOSAA.4.000629

Creep under constant stress

A constant true stress along one axis, the other faces at the ambient pressure: the strain against time, the creep compliance J(t) = ε(t)/σ, the creep rate at the end and the lateral contraction.

Berendsen per axis: dε_k = −(β/τ_p)(P_target − P_kk) dt/3 every ten steps; the barostat's τ_p sets how fast the cell answers, so the first few τ_p are its response

When to use it

Rubber and fibre–matrix composites under a dead load: the elastic part, the early creep and how much the lateral faces contract; compare compounds at one stress and temperature. Export writes it for LAMMPS too (protocol Creep at constant stress): fix npt with the stressed axis at −σ and the others at P, axes uncoupled, strain against time in creep.dat.

Source

core/src/mechanics.cpp

Tested by

tests/test_mechanics.cpp Mechanics.StressControlAndCreep (an argon crystal: P_zz held at −50 MPa, the lateral axes at 1 atm, strain and lateral contraction positive)

Departure from the reference

nanosecond times only: the long-time creep of a polymer is far beyond MD; the barostat (not the inertia of the sample) sets the first response

Sliding friction between walls

One wall molecule slid along x at a set speed over the film between it and a held wall, at the gap the structure has: the force the film puts on the moving wall gives the friction force, the shear stress at the wall and the friction coefficient at that load.

F the summed force of the other atoms on the moving wall's atoms, averaged over the steady part (block errors); A the cell's xy area; the moving wall's height fixed (constant gap)

When to use it

Rubber against a fibre or filler surface, a lubricant or processing aid between surfaces: the stress to slide at a given gap, compared across surface treatments (silane, sizing) at one speed and temperature.

Source

core/src/mechanics.cpp

Tested by

tests/test_mechanics.cpp Mechanics.SlidingFrictionOfALjFilm (an argon liquid between argon walls at γ̇* ≈ 0.5: 47 MPa against the no-slip estimate η γ̇ ≈ 50 MPa, the temperature held)

Departure from the reference

MD speeds (m/s and up) far exceed a tribometer's; the thermostat acts on the flow direction too; constant gap, not constant load

Solvation free energy (thermodynamic integration)

The free energy of moving one molecule from the gas into the rest of the frame: its interactions with everything else are switched off in windows on a copy, first its charges and then its Lennard-Jones with soft-core pairs, and the mean derivative in each window is integrated over λ.

λ = 0 the solute decoupled (gas), λ = 1 coupled (solvated), so the integrals give G(solvated) − G(gas); trapezoid rule over the windows; solute–rest pairs only (its internal interactions stay on in both states); α = 0.5; each soft-core pair shifted to zero at the cut-off; pairwise (DSF) electrostatics, no tail correction; errors from block averages per window

When to use it

Solubility of a plasticiser, curative, antioxidant or water in a rubber matrix compared across matrices (the difference of ΔG between two matrices is the transfer free energy and partition coefficient, ln K = −ΔΔG/RT); hydration free energies of small molecules to check a force field.

Source

core/src/free_energy.cpp

Tested by

tests/test_dynamics.cpp Dynamics.SolvationTiMatchesWidomInsertion (an LJ particle in an LJ liquid: TI against Widom insertion, within the errors); Dynamics.AlchemicalPairsAndTheirDerivatives (λ = 1 equals the plain pairs; ∂U/∂λ and the forces against finite differences)

Departure from the reference

no long-range (PME) electrostatics or tail correction during the windows; NVT at the frame's volume (no pressure–volume term, small for a dense phase); short windows give large errors for big solutes

References

1. Kirkwood, J. G., "Statistical mechanics of fluid mixtures", J. Chem. Phys. 3, 300–313 (1935). doi:10.1063/1.1749657
2. Beutler, T. C., Mark, A. E., van Schaik, R. C., Gerber, P. R., van Gunsteren, W. F., "Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations", Chem. Phys. Lett. 222, 529–539 (1994). doi:10.1016/0009-2614(94)00397-1

Fitting a curve

Any curve in Analyze fitted over a chosen x range: a line or quadratic (solved directly), a power law, an exponential decay, a stretched exponential (KWW), an Arrhenius law, or Gaussian peaks on a baseline (Levenberg–Marquardt). Each parameter with its standard error, R², and the fitted line drawn over the curve.

J the Jacobian at the minimum (central differences); log–log curves are fitted on the logarithms of both axes

When to use it

The power-law exponent of an MSD, a KWW relaxation time and β of a stress or bond autocorrelation, an activation energy from rates against temperature, peak positions and widths of a scattering curve.

Source

studio/CapsStudio/ViewModels/CurveFit.cs

Tested by

the Studio self-test fits back an exponential decay (a, τ, c to 10⁻⁴) and two Gaussians (centres to 10⁻³)

Departure from the reference

the errors assume independent points with equal variance: points of a time correlation are not independent, so treat the ± as a lower bound

Crystallinity index from a diffraction curve

A simulated (or imported) WAXS curve over q fitted with Gaussian peaks on a baseline: peaks at least as wide as the halo width are the amorphous halo, narrower ones crystalline reflections, and the crystallinity index is the crystalline share of the peak area.

σ ≥ the halo width (0.12 Å⁻¹ by default) → amorphous; a finite cell broadens crystalline peaks by about 2π/L, so a small cell needs a wider threshold

When to use it

Strain-induced crystallisation of natural rubber, semicrystalline fibres (PET, PA, cellulose) in a composite: the crystalline share of a simulated cell, or of a measured pattern for comparison.

Source

studio/CapsStudio/ViewModels/AnalyzeViewModel.cs

Tested by

the Studio self-test: a narrow and a wide Gaussian of known areas give X_c to 0.1 %

Departure from the reference

no Lorentz or polarisation corrections and no diffuse disorder term (Ruland's k): an index for comparison, not an absolute crystallinity

References

1. Ruland, W., "X-ray determination of crystallinity and diffuse disorder scattering", Acta Crystallographica 14, 1180–1185 (1961). doi:10.1107/S0365110X61003429
2. Hermans, P. H., Weidinger, A., "Quantitative X-ray investigations on the crystallinity of cellulose fibers. A background analysis", Journal of Applied Physics 19, 491–506 (1948). doi:10.1063/1.1698162

Where a sorbate sits: density maps

During each GCMC pressure's production steps the sorbed molecules' centres are counted on a grid over the cell; the map shows, projected onto a cell face, where in the host the gas or small molecule actually resides — free-volume pockets, filler surfaces, the interphase.

about 2,000 samples per pressure; the map summed over the cell is the mean loading

When to use it

Gas barrier of rubbers and composites: whether CO₂, O₂ or water gathers at a filler surface or in the matrix's free volume, and how that changes with pressure.

Source

core/src/sorption.cpp

Tested by

tests/test_sorption.cpp Sorption.DensityMapOfAnIdealGas (the map sums to the loading and stays even for an ideal gas)

Departure from the reference

the host is fixed (no swelling); centres only, not the molecules' orientation