import cv2

import numpy as np particle_id = 1

import [Link] as plt for cnt in contours:
# -------------------------------- area = [Link](cnt)
# 1. Load grayscale image # Filter noise
img = [Link]("[Link]", if area < 5:
cv2.IMREAD_GRAYSCALE)
continue
if img is None:
M = [Link](cnt)
raise ValueError("Image not found")
if M["m00"] == 0:
# --------------------------------
continue
# 2. Threshold (black dots -> white)
cx = M["m10"] / M["m00"]
binary = [Link](
cy = M["m01"] / M["m00"]
img, 0, 255,
[Link]((particle_id, cx, cy))
cv2.THRESH_BINARY_INV +
cv2.THRESH_OTSU) particle_id += 1

# -------------------------------- # --------------------------------

# 3. Remove small noise # 6. Print particle positions

# -------------------------------- # --------------------------------

kernel = [Link]((3, 3), np.uint8) print("Particle ID | x (pixels) | y

(pixels)")
binary = [Link](binary,
cv2.MORPH_OPEN, kernel) print("--------------------------------------------"
)
# --------------------------------
for pid, x, y in positions:
# 4. Find contours
print(f"{pid:^11} | {x:12.2f} | {y:12.2f}")
# --------------------------------
# --------------------------------
contours, _ = [Link](binary,
# 7. Visualization
cv2.RETR_EXTERNAL,
# --------------------------------
cv2.CHAIN_APPROX_SIMPLE)
img_color = [Link](img,
# -------------------------------- cv2.COLOR_GRAY2BGR)
# 5. Detect particle positions for pid, cx, cy in positions:
# -------------------------------- [Link](img_color, (int(cx), int(cy)), 4,
(0, 0, 255), -1)
positions = [] # (particle_id, x, y)
 [Link](img_color, str(pid),  inverted threshold
(THRESH_BINARY_INV), so:
(int(cx)+5, int(cy)-5),

cv2.FONT_HERSHEY_SIMPLEX,

{
binary (x , y )=
0 if I (x , y)>T
255 if I (x , y)≤T

0.4, (255, 0, 0), 1)  This makes black dots (low

intensity) → white spots in the
[Link](img_color) binary image.
[Link]("Detected Particle Positions")
 THRESH_BINARY_INV +
[Link]("off") THRESH_OTSU:
o THRESH_BINARY_INV:
[Link]()
inverts the mapping.
o THRESH_OTSU: tells
 cv2 → OpenCV: library for image
OpenCV to automatically
processing (reading images, thresholding,
pick the best threshold T by
contours, etc.)
minimizing intra-class
variance.
 numpy (as np) → Handles arrays,
matrices, numerical operations.
Otsu’s method (idea only):
It assumes the histogram has two peaks
 [Link] (as plt) →
(background and object) and chooses T that
For plotting and showing images. best separates them statistically.
[Link]("[Link]",
o particles are black,
cv2.IMREAD_GRAYSCALE): background is white.
o Contours are found on white
 Reads the image from file. regions → so we want
 IMREAD_GRAYSCALE converts it particles to become white.
to a single-channel grayscale o Inversion flips black→white
image, where each pixel is a number and white→black, which is
from 0 to 255. what we need.
o 0 = black
o 255 = white Kernel : This creates a 3×3 matrix of ones:
o values in between = gray.

[ ]
 if the original color image had 1 1 1
channels ( R , G , B ) , grayscale kernel= 1 1 1
typically uses a weighted sum like: 1 1 1

I (x , y )=0.299 R+0.587 G+ 0.114 BGoal:

 It defines a small neighborhood
Convert the grayscale image into a binary around each pixel.
image:
 Opening = erosion followed by
 Background → black dilation
 Particles → white
 m pq=∑ ∑ x y f ( x , y )
p q
1. Erosion: shrinks white regions.
Any small isolated white specks get x y

removed.
2. Dilation: grows white regions back, where f (x , y )is 1 inside the particle, 0
but noise that disappeared doesn’t outside (for binary).
come back.
Important ones here:
Intuitively:
 m00: area (sum of all ones)
 Noise → small → destroyed by  m10: ∑ x
erosion → gone.  m01: ∑ y
 Real particles → larger → slightly
shrunk, then restored in dilation. Then the centroid (center of mass) is:

A contour is essentially a curve that traces m10 m01

the boundary of a connected white region. xc= , y c=
m00 m00
Each white blob (particle) in the binary
Concept Summary
image will ideally give you one contour.
1. Read image → grayscale matrix.
Parameters:
2. Threshold → binary image (0 or
255).
 binary → the cleaned binary
3. Morphological opening → clean
image.
noise.
 cv2.RETR_EXTERNAL:
4. Find contours → each contour ≈
o Only retrieves the outermost
one particle.
contours.
5. Compute centroid ( x c , y c )using
o If a shape has a hole, we
image moments:
ignore the inner contour.
o For isolated dots, outer
m10 m01
contour is enough. xc= , y c=
 cv2.CHAIN_APPROX_SIMPLE: m00 m00
o Compresses
horizontal/vertical/diagonal 6. Store & print particle positions.
segments and saves memory. 7. Visualize positions on image.
o Instead of storing every pixel
on the boundary, it stores
only key points.

 [Link](cnt) returns a
dictionary of spatial moments, e.g.: how to do amplitude comparision and
o m00, m10, m01, m20, m02, wavelength comparision
etc.
Here is the clearest, step-by-step
Mathematically, for a binary region: procedure to do Amplitude Comparison
and Wavelength Comparison between two
(or more) stripe images using ImageJ +
Origin.
This is exactly the workflow used in surface Step 3 — Generate the intensity plot
patterning / wrinkling / optics papers.
Analyze → Plot Profile

A window will appear.

⭐ OVERVIEW OF WHAT YOU MUST
DO Click:

For each image, extract:  List → save as

Image1_Profile.txt
1. Wavelength (λ)  Repeat for Image 2
– From FFT or peak–peak distance
2. Amplitude (A)
– From sinusoidal fitting
🔶 PART 2 — FIND WAVELENGTH (λ)
After you have A₁, A₂ and λ₁, λ₂, you
create two bar/line plots: OPTION A (easy): Direct measurement
in ImageJ
 Amplitude Comparison Plot
 Wavelength Comparison Plot 1. From the line profile → peaks are
visible
2. Measure distance between peaks:
o Analyze → Measure
🔶 PART 1 — EXTRACT VALUES 3. Average multiple peak distances:
FROM EACH IMAGE (ImageJ)
1
λ= ∑(x i +1−x i)
Step 1 — Set the scale N

1. Draw a line on the scale bar

2. Analyze → Set Scale OPTION B (better): FFT method
3. Enter known length (e.g., 20 µm)
1. Save the profile as TXT from ImageJ
2. Import into Origin (File →
Import → Single ASCII)
Step 2 — Draw a profile line 3. Run:
o Analysis → Signal
The line MUST be perpendicular to the Processing → FFT
stripes. 4. Locate dominant frequency f

Use: 1
λ=
f
 Straight Line Tool
 Make it long enough to cross many
periods (15–25 stripes) 🔶 PART 3 — FIND AMPLITUDE (A)
BEST METHOD: Sinusoidal fit in Origin Image Amplitude (A) Wavelength (λ)
Image 1 A₁ λ₁
Import profile (Intensity vs Distance) into Image 2 A₂ λ₂
Origin.
Example (from my automatic calculation
Step 1 — Open nonlinear fitting earlier):
Analysis → Fitting → Nonlinear
Image A (intensity) λ (µm)
Curve Fit → Open Dialog
1 77.98 53.03
Step 2 — Choose model 2 76.26 5.73

Go to:
🔶 PART 5 — PLOT AMPLITUDE
Function → Built-in → Sinusoidal → Sine COMPARISON

Which fits: METHOD 1: Bar Graph (most common)

I (x)= y 0+ A sin ⁡(2 πfx+ϕ ) 1. Highlight Image and Amplitude

columns
2. Plot → Column/Bar →
Step 3 — Initial guesses
Column
 y 0 ≈ mean intensity
Set:
 A ≈(I max −I min )/2
 f ≈ 1/ λ  Title: Amplitude Comparison
 ϕ=0  Y-axis label: Amplitude (Intensity
Units)
Step 4 — Fit

Click Fit
METHOD 2: Line/Scatter Plot
You will get:
1. Highlight both columns
 Amplitude A 2. Plot → Basic 2D → Scatter
 Frequency f → wavelength λ 3. Connect points with a line
 Fit errors (± value)

Repeat for Image 2.

🔶 PART 6 — PLOT WAVELENGTH
COMPARISON

🔶 PART 4 — BUILD THE Exactly same steps:

COMPARISON TABLE
1. Highlight Image and Wavelength
Create a worksheet in Origin like this:
 2. Plot → Column/Bar → ✅ STEP 5 — Format the plot
Column
Double-click plot to open Plot Details:
Title: Wavelength Comparison
Y-axis label: Wavelength (µm)  Change bar color (Fill tab)
 Add title: right-click → Add Text
 Edit axes labels in Axis Dialog

⭐ FINAL OUTPUT YOU SHOULD Set y-axis label:

HAVE
Amplitude (intensity units)
You will produce two graphs:

1. Amplitude Comparison Graph

 X-axis: Image 1, Image 2

 Y-axis: amplitude A
 (Optional) error bars from fitting

2. Wavelength Comparison Graph

 X-axis: Image 1, Image 2

 Y-axis: wavelength λ (in µm)

These two graphs are what you will include

in

If you have amplitude errors (standard error

or standard deviation), create a 3rd column:

Image Amplitude Error

Image 1 77.5 1.0
Image 2 75.8 1.0

Then:

1. Click on the plotted bars to select

the dataset
2. Right-click → Add Error Bars Option 1: Integrated Multi-
3. Choose Y Error → Select the Error
column Scale Approach
Now your bar graph will match the style in "Multi-Scale Polymer Engineering: From
publications. Thin Film Instabilities to Nanocomposite
Dynamics and Electrospun
Architectures"
Strengths: "Mechanical Behavior and Dynamic
Response of Polymer Systems Under
 Emphasizes integration across scales Confinement and External Fields"
 Shows breadth of expertise
 Action-oriented ("Engineering") Strengths:
 Clear methodology scope
 Very broad, includes all three
Option 2: Fundamental projects
 Emphasizes fundamental
Physics Emphasis mechanisms
 Suitable for interdisciplinary journals
"Interfacial Mechanics and Dissipative  Leaves room for future work
Dynamics in Polymer Systems: Buckling
Instabilities, Nanoparticle Dispersions, Option 5: Processing-
and Electrokinetic Processing"
Performance Connection
Strengths:
"Controlling Polymer Morphology and
 Emphasizes fundamental physics Properties: From Processing Parameters
 Technical and rigorous language in Electrospinning to Interfacial
 Shows experimental + simulation Engineering in Composites"
approach
 Direct connection to conferences Strengths:
(APS, Soft Matter)
 Emphasizes control/design aspect
Option 3: Structure-  Processing + characterization focus
 Very practical and applications-
Function Focus oriented
 Good for industry connections
"Structure-Property Relationships in
Stressed Polymers: Wrinkling
Phenomena, Viscoelastic Composites, and
Fiber Morphology Control" COMPREHENSIVE
Strengths: PRESENTATION
 Material science terminology
 Directly addresses properties
 Clear experimental deliverables
 Highly citable structure OUTLINE WITH
CONTENT
Option 4: Broader Systems
Approach
Full Structure  Title: "What We're
Investigating"
(Recommended 20-25  3-4 fundamental questions,
e.g.,:
slides for defense)  Q1: "How do stressed
films transition from
SECTION 1: smooth to wrinkled?"
 Q2: "How do grafted
MOTIVATION & polymer chains
CONTEXT (3-4 slides) dissipate energy in
composites?"
Slide 1: Title Slide  Q3: "What controls
fiber morphology in
 Your name, advisor, electrospinning?"
institution, date
 Visual: Department logo or
striking polymer image
SECTION 2A: PROJECT
Slide 2: The Challenge
1 - FILM BUCKLING (3-4
 Title: "Why Polymer Physics slides)
Matters"
 Content: 3-4 key societal Slide 5: Film Buckling—Motivation
challenges:
 Need for multi-  Title: "Stress-Induced
functional materials Instabilities in Thin Films"
 Bridging experiment  Visual: AFM images showing
and simulation smooth → wrinkled transition
 Scaling from lab to  Text:
applications  Wrinkling = common
 Visual: Industrial applications failure mechanism
(filters, composites, coatings)  Controllable pattern
formation opportunity
Slide 3: Your PhD Scope  Relevant to: devices,
sensors, templates
 Title: "Three Connected
Research Threads" Slide 6: Experimental Approach
 Show framework diagram
[see Chart #41 I generated]  Title: "From Stress
 Text: Brief 2-3 line Application to Pattern
description of each project Imaging"
 Visual: Three distinct  Flow: Sample → Apply
icons/images for each project Uniaxial Stress → AFM Scan
→ Image Analysis
Slide 4: Research Questions
  Characterization tools: AFM,  Visual: Schematic of
SEM, Optical microscope nanoparticle (core) + grafted
 Visual: Schematic of stress polymer chains
apparatus + sample holder  Key: Polymer chains act as
mechanical dampers
Slide 7: Key Results  Characterization methods:
DLS, Rheology, TEM, MD
 Title: "Wavelength Scaling
with Film Thickness" Slide 10: Rheological Response
 Visual: 2-3 AFM images
showing wrinkle patterns at  Title: "Viscoelastic Behavior:

Key finding: λ ∝ h^α

different stresses Storage vs Loss Modulus"
  Visual: G' and G'' vs
(wavelength vs thickness frequency plot
power law)  Text: PGNPs show enhanced
 Add quantitative plot of damping compared to
wavelength vs stress unfilled polymers
 Key metric: Loss angle (tan
Slide 8: Simulation & Theory δ) increases with PGNP
loading
 Title: "Molecular Dynamics
Insights" Slide 11: Molecular Dynamics
 LAMMPS simulation Analysis
methodology (coarse-grained
polymer chains)  Title: "Chain Motion & Energy
 Show: Stress-strain curves Dissipation at Nanoscale"
from MD vs experiment  Visual: Molecular trajectory
 Key insight: Viscoelastic snapshots or mean-squared
effects dominate at certain displacement
strain rates  Text: Grafted chains exhibit
two relaxation timescales
 Finding: Energy dissipation
proportional to chain mobility
SECTION 2B: PROJECT Slide 12: Application Impact
2 - POLYMER-GRAFTED
NANOPARTICLES (3-4  Title: "Toward Smart
Damping Materials"
slides)  2-3 applications: Vibration
isolators, shock absorbers,
Slide 9: PGNPs—System Overview aerospace materials
 Visual: Photos or schematics
 Title: "Dissipation in Polymer- of end-use devices
Grafted Nanoparticle
Composites"
SECTION 2C: PROJECT Slide 16: Applications & Future
Work
3 - ELECTROSPINNING
(3-4 slides)  Title: "From Lab Fibers to
Functional Materials"
 Applications: Air filters (μm
Slide 13: Electrospinning—
range), tissue scaffolds, drug
Phenomena
delivery vehicles
 Future directions: Composite
 Title: "From Jet to Fiber: fibers, core-sheath
Morphology Control" structures, hybrid assemblies
 Visual: Electrospinning setup
diagram (needle, voltage,
collector)
 Text: High electric field →
polymer jet → ultrathin fiber SECTION 3:
 Morphologies: Continuous
fibers, beads, beads-on-
INTEGRATION &
fibers IMPACT (2-3 slides)
Slide 14: Parameter Optimization Slide 17: Connecting the Three
Projects
 Title: "Controlling Fiber
Morphology"  Title: "Unified Perspective:
 Visual: Matrix showing effect Polymers Under Stress"
of voltage, feed rate,  Text:
concentration on structure  All three projects
 Key result: Transition map involve polymer
(beads → beads-on-fibers → deformation/processing
continuous fibers)  Combine experiment,
 Characterization: SEM simulation,
images + fiber diameter characterization
distributions  Bridge macroscopic
behavior ↔ molecular
Slide 15: Structural mechanisms
Characterization  Visual: Flow diagram showing
connections (interfacial
 Title: "Fiber Properties and mechanics, polymer
Analysis" dynamics, processing
 Visual: SEM images at effects)
different magnifications
 Metrics: Average fiber Slide 18: Broader Contributions
diameter, size distribution,
porosity  Title: "Significance to
 Compare: Expected vs Polymer Science"
observed (from theory)  3 contributions:
 1. New understanding of Slide 21: Thank You & Questions
stress-induced pattern
formation  Large, centered "Thank You!"
2. Molecular mechanisms  Your contact info, lab
of dissipation in affiliation
composites  Invitation to ask questions
3. Process-structure-
property relations in
electrospinning
Visual: Icons or small images

representing each
PRACTICAL
contribution PRESENTATION
Slide 19: Future Directions & Open
DELIVERY TIPS
Questions
Narrative Structure
 Title: "What's Next?" (Swath & Dive Pattern)
 4-5 possibilities:
 Can buckling be used
1. SWATH (10 min / ~10 slides):
for self-assembly?
 High-level overview of
 Can we engineer
all three projects
PGNPs for
 Why they matter, how
temperature-
they're connected
responsive damping?
 What methods you
 Multi-material
used
electrospun fibers?
 Preliminary big-picture
 Machine learning to
findings
optimize
electrospinning 2. DIVE (15 min / ~8-10 slides):
parameters?  Select ONE project
 In-situ characterization (typically the most
during stress developed)
application?  Go into detailed
methodology
 Show raw data,
Slide 20: Conclusion
analysis steps, derived
 Title: "Summary: From insights
 Trace logic from
Understanding to Innovation"
 Key message: Polymers are experiment →
versatile; understanding observation →
mechanics at all scales interpretation
enables smart design 3. RECAP & OUTLOOK (3 min /
 Visual: Composite image of ~2 slides):
three projects  Synthesize all three
projects
  Future directions, next one answers (Story-
broader impact driven approach)
 Closing statement 5. Visuals > Text: Use AFM/SEM
images, molecular
Delivery Strategy snapshots, plots rather than
bullet-heavy text
 Practice timing: 1 minute per
slide average
 Pause points: After major
results, before difficult RECOMMENDED
concepts SECTION HEADERS
 Audience engagement: Ask
rhetorical questions ("What FOR YOUR SPECIFIC
happens if we increase RESEARCH
voltage?")
 Backup slides: 2-3 extra
detailed slides for tough For Project 1 (Film
questions (experimental Buckling):
data, supplementary
analysis)  "Surface Instabilities as a
Design Tool"
 "Elasticity Meets
Viscoelasticity: Wrinkling in
REFERENCES FOR Polymer Films"
"From Microscopic Stress to
DESIGNING YOUR 
Macroscopic Patterns"
SLIDES
For Project 2 (PGNPs):
From the research I found:
 "Molecular Engineering of
1. Data-to-Ink Ratio: Maximize Damping Materials"
information density;  "Soft Matter Mechanics:
minimize non-data ink Polymer Chains as Energy
(decorations, gridlines) Dissipators"
2. Color Science: Dark theme  "From Nanoparticles to
(slate/teal) is professional; Functional Composites"
ensure 4.5:1 contrast ratio
for text readability
3. Typography: Use sans-serif For Project 3
(Arial, Helvetica) for (Electrospinning):
presentations; 18pt minimum
for body text
 "Electrokinetics Meets
4. Narrative Flow: Each slide Polymer Physics"
should raise a question the
  "Process Control in Nanofiber
Fabrication"
 "Engineering Fiber This integrated presentation
Architecture Through Field framework, combined with the
Control" SOTA slides above, will position
your research for a compelling
defense and future publications.
Focus on narrative coherence,
FINAL CHECKLIST visual clarity, and connecting the
three projects through shared
BEFORE DEFENSE physics principles.

✅ Content:

 Each slide title is a statement

(not a question) 1. [Link]
 All figures have captions and 711819211/Polymer-
unit labels Nanocomposites-Thesis-PDF
 Equations appear only when 2. [Link]
discussed articlehtml/2025/sc/d4sc04600e
 Results shown before 3. https://
conclusions [Link]/
doi/full/10.1002/aesr.202400077
4. [Link]
✅ Design:
literature-selections/nanocomposites-
materials/dissertation/
 Consistent color scheme
5. [Link]
throughout
10.63823/20250105
 Font sizes ≥18pt for body, 6. [Link]
≥24pt for titles polymers-2026
 Dark background (reduces 7. [Link]
eye strain during long full/10.1080/17458080600670015
presentations) 8. [Link]
 No more than 4-5 bullet 9. https://
points per slide [Link]
10. [Link]
✅ Delivery: science/article/abs/pii/
S1359029403000049
 Practiced out loud (20-25 11. [Link]
slides in 45-50 minutes) 3523cb76-af24-43ce-89ff-
 Transitions are smooth and c43a4587495d
logical 12. https://
 Backup slides prepared for [Link]/
likely questions abstract/
 Timer set for your defense 13. [Link]
duration articles/PMC6589095/
14. [Link] 32. [Link]
album_invisible/2024/tae-yeon- v=TF2D6zlme2c
kong-ph-d-disseration-defense 33. [Link]
15. [Link] tag/data-analysis/
members-activities/ICPP- 34. https://
10_Program_web.pdf [Link]
16. https:// m/best-poster-presentation
[Link]/ 35. [Link]
doi/full/10.1002/adma.202500162 v=CM1gcTkYje0
17. [Link] 36. [Link]
polymera2026 slideshow/data-visualizationpptx-
18. [Link] 252223193/252223193
9/1/011319/2835574/ 37. [Link]
Electrospinning-research-and- content/uploads/2016/03/
products-The-road-and Good_Presentation_Guidelines_2pag
19. [Link] [Link]
jacsau.1c00297 38. [Link]
20. [Link] g=883178
articles/PMC12295322/ 39. [Link]
21. [Link] 40. [Link]
designing-poster-presentation uploads/2025/03/2025-Guidelines-
22. [Link] for-Successful-Scientific-
2024/12/03/my-tips-for-scientific- Presentations-_Updated2025.pdf
conference-presentations/
23. [Link]
[Link]/
eposters
24. https://
[Link]/
2024/speaker-guidelines
25. [Link]
presentation/
26. [Link]
dive/
27. [Link]
v=lIv2vDE6ENo
Overview of the
28. [Link] Simulation Approach
creating-a-scientific-presentation
29. [Link] The paper uses a bead-spring
slideshow/prepare-your-phd-defense- polymer model simulated with the
presentation/7445444 ESPResSo++ MD package to
30. [Link] investigate equilibrated polymer
visualization-slide-types melts. The simulations start from
31. [Link] previously equilibrated high
2024/[Link] molecular weight polymer melts
using a hierarchical backmapping 2. Lennard-Jones (LJ) Potential
procedure. The excluded volume interactions
between all monomers (both intra-
Model Details and and intermolecular):

[( ) ( ) ]
12 6
σ σ
System Setup U LJ ( r)=4 ϵ
r
−
r
This is

typically truncated at a cutoff

Polymer Model distance rc = 2.5σ.

The polymers are modeled as 3. Bond Bending Potential

semiflexible bead-spring chains Controls chain stiffness through
with specific parameters: three-body interactions:
U bend =k θ [1−cos ⁡θ] where θ is the
 Chain lengths: N = 500, angle between consecutive bond
1000, and 2000 monomers vectors, and kθ = 1.5ε provides
 Number of chains: nc = 1000 semiflexibility.
chains
 Monomer density: ρ = Step-by-Step
0.85σ⁻³ (where σ is the LJ
length unit) Reproduction
 Bending stiffness: kθ = 1.5 Protocol
(primary parameter used)
 Root-mean-square bond
length: ℓb ≈ 0.964σ Step 1: Initial
 Flory characteristic ratio: C∞ Configuration
≈ 2.88
Generation
Force Field What's happening: You need to
Components create the initial polymer melt
structure. Since generating
The total potential energy consists equilibrated high molecular weight
of three terms: polymer melts is extremely
challenging, the paper uses pre-
1. Bond Stretching Potential (FENE) equilibrated configurations
The finitely extensible nonlinear obtained through hierarchical
elastic (FENE) potential prevents backmapping.
bond crossing and maintains chain
connectivity: Implementation:

( )
2
−1 2 r
U FENE= k 0 R 0 ln ⁡ 1 − 2 where  Start with a soft-sphere
2 R0
coarse-grained model at low
typical values are k₀ = 30ε/σ² and resolution
R₀ = 1.5σ.
  Gradually increase resolution  Run NPT (constant pressure-
through sequential temperature) MD simulation
backmapping  Temperature: T = 1.0 (in
 Final step: Apply full MD reduced units, kBT/ε = 1.0)
simulation with the bead-  Pressure: Start at high
spring model pressure (e.g., P = 200 atm)
 Alternative: Build random to speed equilibration
polymer chains in a  Time step: δt = 0.01τ (where
simulation box and use a τ is the LJ time unit: τ =
"push-off" procedure σ√(m/ε))
followed by extensive  Duration: Run until density
equilibration stabilizes (~10⁵-10⁶ time
steps)
Use Langevin thermostat or
Step 2: Energy 
Nose-Hoover thermostat
Minimization
What's happening: Remove any
Step 4: Equilibration
overlaps or high-energy at Target Density
configurations from the initial
structure. What's happening: Equilibrate at
the final target density ρ =
Implementation: 0.85σ⁻³.

 Use steepest descent or Implementation:

conjugate gradient
minimization  Switch to NVT (constant
 Minimize until forces are volume-temperature)
below a threshold (e.g., ensemble or continue NPT at
Fmax < 0.01ε/σ) lower pressure
 This prevents numerical  Temperature: T = 1.0
instabilities during MD  Continue until all chain
conformations are relaxed
Check equilibration criteria:
Step 3: Initial 

end distance ⟨R²e⟩

 Mean square end-to-
Equilibration at High
Pressure reaches steady value

gyration ⟨R²g⟩
 Mean square radius of

What's happening: Equilibrate the

system at the target temperature stabilizes
 Potential energy
while allowing the density to
adjust. converges
 Required time: At least
Implementation: several reptation times τd ~
τ₀N³·⁴/Ne
Characteristic Time Scales: Step 6: Production
 τ₀ ≈ 2.89τ (elementary time) Run for Dynamic
 τe ≈ 1.98×10³τ Properties
(entanglement time)
 τR ≈ 6.44×10⁵τ for N=500
What's happening: Monitor time
(Rouse time)
evolution to calculate dynamic
 τd ≈ 2.97×10⁷τ for N=500
properties.
(disentanglement time)
Implementation:
Step 5: Production
Run for Static  Continue NVT MD at T = 1.0
 Record monomer positions
Properties frequently (every 10-100
time steps)
What's happening: Collect  Track positions up to time t
equilibrium configurations to ~ 10⁷τ
calculate static properties.  Calculate mean square
displacements:
Implementation:  g₁(t): MSD of inner
monomers
 Run NVT MD at T = 1.0  g₂(t): MSD relative to
 Save configurations every center of mass
~10³-10⁴ time steps  g₃(t): MSD of center of
 Duration: Multiple τd to mass
ensure independent samples
 Calculate: Expected Scaling Behaviors:

end distance: ⟨R²e⟩

 Mean square end-to-
 t < τ₀: g₁(t) ~ t¹ (ballistic)

gyration: ⟨R²g⟩
 Mean square radius of  τ₀ < t < τe: g₁(t) ~ t^(1/2)
(Rouse)

distance: ⟨R²(s)⟩ vs
 Mean square internal  τe < t < τR: g₁(t) ~ t^(1/4)
(constrained
separation s Rouse/reptation)

⟨cos θ(s)⟩
 Bond-bond correlation:  τR < t < τd: g₁(t) ~ t^(1/2)
(reptation)
 Structure factor: Sc(q)  t > τd: g₁(t) ~ t¹ (diffusion)

Expected Results: Step 7: Primitive

⟨R²e⟩/⟨R²g⟩ ≈ 6 (ideal chain

Path Analysis (PPA)
⟨R²e⟩ ∝ N and ⟨R²g⟩ ∝ N
behavior)
What's happening: Identify

(Gaussian scaling) entanglements by finding the
shortest path connecting chain
ends while avoiding chain  Stress relaxation modulus:
crossings. G(t) = (V/kBT)⟨SAF(t)⟩

Implementation: Method 2 - Step Strain:

 Fix all chain endpoints in  Apply small uniaxial

space elongation (λ ≈ 1.2)
 Turn off intrachain excluded  Monitor normal stress decay
volume and bending σnorm(t)
interactions  Calculate: G(t) =
 Keep interchain interactions σnorm(t)/(λ² - 1/λ)
active
 Minimize system energy by Expected Results:
"cooling" to T = 0
 The resulting contracted  t < τe: G(t) ~ t^(-1/2)
paths are the primitive paths (Rouse regime)
 Calculate entanglement  τe < t << τd: G(t) ≈ G⁰N =
length: Ne,PPA = (4/5)(ρkBT/Ne) (plateau
ℓK^(pp)/ℓb^(pp) modulus)
 Entanglement length from
Physical Meaning: The primitive plateau: Ne ≈ 28±2 for kθ =
path represents the average tube 1.5
confinement experienced by each
chain. The entanglement length
Ne,PPA ≈ 28 monomers for kθ =
Key Physical
1.5. Insights
Step 8: Stress Equilibration Challenge: The
longest relaxation time τd scales
Relaxation and as N³·⁴, making equilibration of
Viscoelasticity long chains computationally
expensive. For N = 2000, you need
What's happening: Calculate the simulations extending to ~10⁸τ
stress relaxation modulus G(t) to time units.
characterize viscoelastic
properties. Ideal Chain Behavior: Semiflexible
chains (kθ = 1.5) show near-ideal
Method 1 - Green-Kubo Relation: Gaussian statistics in the melt due
to screening of excluded volume
 Calculate off-diagonal stress interactions, with minimal
tensor components σαβ(t) deviations from ideality compared
 Compute stress to fully flexible chains.

SAFαβ(t) = ⟨σ̄αβ(t)σ̄αβ(0)⟩
autocorrelation function:
Entanglement Effects: The
entanglement length Ne ≈ 28
defines the crossover from Rouse positions = [start_pos]
to reptation dynamics, with chains
confined to tube-like regions for N for i in range(n_monomers - 1):
>> Ne.
# Random direction
This comprehensive protocol
allows you to reproduce the static theta = [Link](0, [Link])
and dynamic properties reported in
the paper, providing deep insights phi = [Link](0, 2*[Link])
into polymer melt physics across
multiple time and length scales.

# Python script to create initial polymer melt dx = bond_length * [Link](theta) *

configuration [Link](phi)

import numpy as np dy = bond_length * [Link](theta) *

[Link](phi)

dz = bond_length * [Link](theta)
# System parameters

n_chains = 1000
new_pos = positions[-1] +
n_monomers = 500 # Can be 500, 1000, or [Link]([dx, dy, dz])
2000
[Link](new_pos)
density = 0.85

sigma = 1.0
return [Link](positions)

# Calculate box size

# Generate all chains
n_total = n_chains * n_monomers
all_positions = []
volume = n_total / density
all_molecules = []
box_length = volume**(1/3)
all_bonds = []

bond_id = 1
# Create random walk polymer chains
atom_id = 1
def create_polymer_chain(start_pos,
n_monomers, bond_length=0.97):
for chain_id in range(n_chains):

# Random starting position # Write LAMMPS data file

start_pos = [Link](0, with open('polymer_melt.data', 'w') as f:

box_length, 3)
[Link]('# Polymer melt system\n\n')

[Link](f'{n_total} atoms\n')
# Create chain
[Link](f'{len(all_bonds)} bonds\n')
chain_pos =
create_polymer_chain(start_pos, [Link](f'{n_total - n_chains} angles\n\n')
n_monomers)

[Link]('1 atom types\n')

# Wrap into box
[Link]('1 bond types\n')
chain_pos = chain_pos % box_length
[Link]('1 angle types\n\n')

# Store positions and molecule IDs

[Link](f'0.0 {box_length} xlo xhi\n')
for pos in chain_pos:
[Link](f'0.0 {box_length} ylo yhi\n')
all_positions.append(pos)
[Link](f'0.0 {box_length} zlo zhi\n\n')
all_molecules.append(chain_id + 1)

[Link]('Masses\n\n')
# Create bonds (connect consecutive
monomers) [Link]('1 1.0\n\n')

for i in range(n_monomers - 1):

all_bonds.append([bond_id, 1, atom_id [Link]('Atoms\n\n')

+ i, atom_id + i + 1])
for i, (pos, mol) in
bond_id += 1 enumerate(zip(all_positions,
all_molecules)):

[Link](f'{i+1} {mol} 1 {pos[0]:.6f}

atom_id += n_monomers {pos[1]:.6f} {pos[2]:.6f}\n')
 #
===============================
[Link]('\nBonds\n\n') =============

for bond in all_bonds: # INITIALIZATION

[Link](f'{bond[0]} {bond[1]} #
{bond[2]} {bond[3]}\n') ===============================
=============

units lj # Use reduced LJ units

# Create angle list
atom_style molecular # Atoms in
[Link]('\nAngles\n\n') molecules with bonds

angle_id = 1 boundary ppp # Periodic

boundaries in all directions
atom_id = 1
neighbor 0.3 bin # Skin distance
for chain_id in range(n_chains): for neighbor list

for i in range(n_monomers - 2): neigh_modify every 1 delay 0 check yes

[Link](f'{angle_id} 1 {atom_id+i}
{atom_id+i+1} {atom_id+i+2}\n')
#
angle_id += 1 ===============================
=============
atom_id += n_monomers
# READ CONFIGURATION

#
print(f'Created polymer melt with ===============================
{n_chains} chains of {n_monomers} =============
monomers')
read_data polymer_melt.data
print(f'Box length: {box_length:.2f} sigma')

# Polymer melt simulation - reproducing

Hsu & Kremer JCP 2016 #
===============================
# Semiflexible bead-spring model with =============
FENE bonds
# FORCE FIELD PARAMETERS
# #
=============================== ===============================
============= =============

# Lennard-Jones potential (all monomers # Compute properties

interact)
compute myTemp all temp
pair_style lj/cut 2.5 # Cutoff at 2.5
sigma compute myPE all pe

pair_coeff 1 1 1.0 1.0 2.5 # epsilon=1.0, compute myKE all ke

sigma=1.0, cutoff=2.5

# Compute chain properties

# FENE bond potential
compute gyration all gyration/molecule
bond_style fene
compute rg_avg all reduce ave
bond_coeff 1 30.0 1.5 1.0 1.0 c_gyration

# Parameters: K=30.0, R0=1.5, epsilon=1.0,

sigma=1.0
# Compute end-to-end distance

compute mol_chunk all chunk/atom

# Bond angle potential (bending stiffness) molecule

angle_style cosine compute end2end all gyration/chunk

mol_chunk
angle_coeff 1 1.5

# k_theta = 1.5 (semiflexible chains)

# Compute MSD

compute msd all msd com yes

#
===============================
=============
# Thermo output
# SETTINGS AND COMPUTES
thermo_style custom step temp press pe ke
etotal density c_rg_avg
thermo 1000 dump 1 dump_equil1 all custom
10000 [Link] id mol type x y z

#
=============================== run 100000
=============
undump 1
# ENERGY MINIMIZATION
unfix 1
#
=============================== unfix 2
=============

minimize 1.0e-4 1.0e-6 1000 10000

#
reset_timestep 0 ===============================
=============

# EQUILIBRATION STAGE 2: Target

# Temperature NPT
===============================
============= #
===============================
# EQUILIBRATION STAGE 1: High =============
Temperature NVT
fix 1 all npt temp 1.0 1.0 1.0 iso 0.0
# 0.0 10.0
===============================
============= timestep 0.01

velocity all create 2.0 482937 mom yes

rot yes dist gaussian
dump 2 dump_equil2 all custom
10000 [Link] id mol type x y z

fix 1 all nve/limit 0.05 run 500000

fix 2 all langevin 2.0 1.0 1.0 498437 undump 2

zero yes
unfix 1

timestep 0.005
#
===============================
============= # Detailed thermo output

# EQUILIBRATION STAGE 3: Long NVT thermo_style custom step temp press pe ke

Equilibration etotal density c_rg_avg

# thermo 5000
===============================
=============

fix 1 all nvt temp 1.0 1.0 1.0 # Run production

timestep 0.01 run 10000000

dump 3 dump_equil3 all custom write_data final_config.data

50000 [Link] id mol type x y z

run 5000000
print "Production run completed!"
undump 3

#
===============================
=============
# Dynamic properties analysis - continuation
# PRODUCTION RUN: Static Properties from equilibrated state

# # Compute MSDs and stress autocorrelation

===============================
=============

reset_timestep 0 read_data final_config.data

# Output detailed trajectories # Force field (same as before)

dump 4 dump_prod all custom 10000 pair_style lj/cut 2.5

[Link] id mol type x y z
pair_coeff 1 1 1.0 1.0 2.5
dump_modify 4 sort id
 # High-frequency output for dynamics

bond_style fene thermo 1000

bond_coeff 1 30.0 1.5 1.0 1.0 thermo_style custom step temp

c_msd_all[4] c_press_components[*]

angle_style cosine
dump 1 dump_dyn all custom 100
angle_coeff 1 1.5 [Link] id mol type x y z vx vy
vz

dump_modify 1 sort id
# Set velocities

velocity all create 1.0 328493 mom yes

rot yes run 10000000

# Compute MSDs write_data final_dynamic.data

compute msd_all all msd com yes

compute msd_inner all msd com no # Primitive path analysis to find

entanglements

# Read equilibrated configuration

# Compute stress for autocorrelation

compute stress all stress/atom NULL

read_data final_config.data
compute press_components all reduce
ave c_stress[1] c_stress[2] c_stress[3]

# Turn off intrachain LJ and angles

# Output for correlation analysis pair_style lj/cut 2.5

fix 1 all nvt temp 1.0 1.0 1.0 pair_coeff 1 1 1.0 1.0 2.5

timestep 0.01

bond_style fene
bond_coeff 1 30.0 1.5 1.0 1.0

print "Primitive path analysis

completed!"
# Disable angles (set k=0 for contraction)

angle_style cosine

angle_coeff 1 0.0

# Freeze chain ends

group ends id <= 1000:500:500000 #

Every 500th atom is chain end

fix 1 ends setforce 0.0 0.0 0.0

# Energy minimization to contract chains

minimize 1.0e-8 1.0e-10 100000

1000000

# Compute primitive path length

compute gyration all gyration/molecule

compute rg_pp all reduce ave

c_gyration

thermo_style custom step pe c_rg_pp

thermo 100

# Output contracted configuration

write_data primitive_path.data