#### NOTE
[Go to the end](#sphx-glr-download-auto-examples-force-balance-plot-force-balance-py)
to download the full example code.
# 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.
```none
Final residual (should be close to 0): 2.012106494051574e-08
Shear Modulus [MPa]: 0.04363931207507527
```
```Python
import os
from pathlib import Path
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)
def _find_fixtures_dir() -> Path:
for cur in (
[Path(__file__).resolve().parent] if "__file__" in globals() else []
) + [Path.cwd()]:
for p in [cur] + list(cur.parents):
candidate = p / "tests" / "pylimer_tools" / "fixtures"
if candidate.is_dir():
return candidate
return Path("tests/pylimer_tools/fixtures").resolve()
_fixtures_dir = _find_fixtures_dir()
universe = read_data_file(
str(_fixtures_dir / "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,
)
```
**Total running time of the script:** (0 minutes 0.169 seconds)
[`Download Jupyter notebook: plot_force_balance.ipynb`](_downloads/1060533adab89eaeb1fc311189732520/plot_force_balance.ipynb)
[`Download Python source code: plot_force_balance.py`](_downloads/482e5042e99228b5a45bc8a3084728e8/plot_force_balance.py)
[`Download zipped: plot_force_balance.zip`](_downloads/0df86cde61cca3638558dd249a5db062/plot_force_balance.zip)
[Gallery generated by Sphinx-Gallery](https://sphinx-gallery.github.io)