Force Balance

This example demonstrates how to use the MEHP implementation with slip-springs to homogenize a polymer network and predict its equilibrium shear modulus.

Final residual (should be close to 0): 3.560906924081682e-08
Shear Modulus [MPa]:  0.054569399499852754

import os

import numpy as np

from pylimer_tools.io.bead_spring_parameter_provider import (
    ParameterType,
    Parameters,
    get_parameters_for_polymer,
)
from pylimer_tools.io.read_lammps_output_file import read_data_file
from pylimer_tools_cpp import MEHPForceBalance, Universe

# Get parameters for conversion factors
params = get_parameters_for_polymer(
    "PDMS", parameter_type=ParameterType.GAUSSIAN)
assert isinstance(params, Parameters)

# Load your network (replace with your file)
universe = read_data_file(
    os.path.join(
        os.getcwd(),
        "../..",
        "tests/pylimer_tools/fixtures/structure/network_100_a_46.structure.out",
    )
)
assert isinstance(universe, Universe)

# 2. MEHPForceBalance: Hookean springs, allows slip-links
mehp_fb = MEHPForceBalance.construct_with_random_sliplinks(
    universe,
    nr_of_sliplinks_to_sample=int(
        params.get_entanglement_density() * universe.get_volume()
    ),
    upper_sampling_cutoff=params.get_sampling_cutoff(),
)
mehp_fb.run_force_relaxation(max_nr_of_steps=1000)

print(
    "Final residual (should be close to 0):", mehp_fb.get_displacement_residual_norm()
)

# apply conversion factors, measure resulting shear modulus
r02_slope = params.get("R02")
r02_slope_magnitude = r02_slope.to(params.get("distance_units") ** 2).magnitude
kbt = params.get("T") * params.get("kb")
gamma_conversion_factor = (
    (kbt / ((params.get("distance_units")) ** 3)).to("MPa").magnitude
)

shear_modulus = (
    gamma_conversion_factor
    * np.sum(mehp_fb.get_gamma_factors(r02_slope_magnitude))
    / universe.get_volume()
)
print(
    "Shear Modulus [MPa]: ",
    shear_modulus,
)