#### NOTE
[Go to the end](#sphx-glr-download-auto-examples-force-balance-plot-force-balance-2-py)
to download the full example code.
# Force Balance 2
MEHPForceBalance2 is a faster implementation of the MEHP force balance method,
which uses static links to model entanglements instead of slip-links.
See Bernhard and Gusev [[BG25](../../acknowledgements.html.md#id12)]
```none
Final residual (should be close to 0): 1.7153887116022568e-13
Shear Modulus [MPa]: 0.11323753391944498
```
```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 MEHPForceBalance2, 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 = MEHPForceBalance2(
universe,
nr_of_entanglements_to_sample=int(
params.get_entanglement_density() * universe.get_volume()
),
upper_sampling_cutoff=params.get_sampling_cutoff(),
entanglements_as_springs=False,
)
mehp_fb.run_force_relaxation()
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.213 seconds)
[`Download Jupyter notebook: plot_force_balance_2.ipynb`](_downloads/01077d6c84a59ddcb66597da8c21b737/plot_force_balance_2.ipynb)
[`Download Python source code: plot_force_balance_2.py`](_downloads/9419b20ac815b8fc487781636094dad8/plot_force_balance_2.py)
[`Download zipped: plot_force_balance_2.zip`](_downloads/d15e58ba856b3d8f5c29b12e83f74b6a/plot_force_balance_2.zip)
[Gallery generated by Sphinx-Gallery](https://sphinx-gallery.github.io)