The puzzle of light-driven solids
Liquid-crystal polymer networks do something that seems almost impossible on paper: absorb a photon, convert that absorbed energy through a chain of coupled physical processes — photo-isomerization, heat generation, mechanical stress — and produce a large, reversible shape change, all within a thin film that weighs almost nothing. The applications are genuinely interesting: soft robotic grippers that operate without motors, tunable optical coatings that switch between states under light, biomedical microactuators where electrical wiring is impractical. The materials science community has been fabricating and characterizing these films for two decades and making real progress.
The simulation side is a different story. Photo-isomerization kinetics, thermal transport, and large-deformation solid mechanics each belong to different solver traditions. Most codes handle one of the three; a few handle two. When I started my PhD at IIT Madras in 2017, no off-the-shelf solver coupled all three within a single continuum framework that an experimentalist could actually hand a geometry to and get a useful answer back. So I wrote one that does.
That choice — building the simulation infrastructure first, then using it to answer scientific questions — shaped everything that followed. Seven years later, the platform had powered contributions to 13 peer-reviewed publications, been adopted by research groups in three countries working entirely from documentation, and earned the Best PhD Thesis Award from IIT Madras, selected from more than 50 candidates. The pages below describe what the platform is, what it revealed, and where it broke down in useful ways.
The platform — coupled photo-thermo-mechanical FEM
The central tool is a custom FORTRAN finite element framework embedded inside Abaqus Standard through four user-subroutine interfaces: UTEMP, UFIELD, UEXTERNALDB, and URDFIL. At roughly 2,200 lines across 15 subroutines, it implements what I think of as a solver-within-solver architecture.
Here is what that means in practice. UTEMP contains a complete 2D Bubnov-Galerkin thermal FEM solver — 4-node bilinear elements, 2×2 Gauss quadrature, Crank-Nicolson time integration. This solver runs inside an Abaqus subroutine, generating temperature fields that Abaqus itself never computes directly. Simultaneously, UFIELD solves nonlinear Corbett-Warner photo-isomerization kinetics through 21 through-thickness layers using Picard iteration, tracking how azobenzene chromophores convert between trans and cis states as light attenuates through the film following Beer-Lambert’s law. The two subroutines hand off their fields to Abaqus’s standard nonlinear mechanical solver, which computes the resulting deformation using composite S4R shell elements. The coupling is one-way sequential — photo-chemistry drives thermal, thermal drives mechanics — which is physically correct for the LCN systems I studied and avoids the convergence problems that come with full two-way coupling.
Three hand-coded linear solvers (Gauss-Seidel, Gaussian elimination, and LU decomposition with caching) are built in, so the framework has zero external numerical dependencies. That turned out to matter for adoption: a group in Eindhoven or Bangalore could run the platform without worrying about library version conflicts.
The validation covered three chromophore types, two operating environments, and light intensities from 10 to 200 mW/cm². One project I’m particularly pleased with: I used a separate 91-layer composite shell configuration to computationally explore 18 alignment-and-taper combinations for an LCN film meant to morph from flat into a cone — a shape inspired by the calla lily flower. The simulation identified the optimal 45-degree diagonal alignment configuration. The TU Eindhoven group then fabricated that specific design. Prediction first, fabrication second. That is the direction of causality I wanted.
A separate FORTRAN code of about 750 lines extended the framework to self-sustained thermal oscillations on hot plates, adding contact detection logic, temperature-dependent material properties, and a spatially varying air-temperature profile. The key finding from that work: only shapes where the center of gravity does not coincide with the midline will sustain oscillations under constant thermal stimulus. That geometric criterion gives a designer a direct rule instead of trial-and-error shape selection.
The most analytically satisfying result from the PhD came from a much simpler model than any of the above. I derived a closed-form curvature expression for bilayer actuators where the active layer covers only a fraction of the passive substrate — a problem that had been treated with numerical methods since Timoshenko’s bimetallic-strip formula in 1925 assumed full coverage. Working from strain-energy minimization with kinematically admissible displacement fields, I found that a single dimensionless parameter, m²n (combining the modulus ratio and the thickness ratio), governs three distinct design regimes. One of those regimes produces a counter-intuitive result: 40% coverage can generate higher curvature than 100% coverage. That was confirmed experimentally by collaborators in subsequent work.
Computational modeling is a recognized bottleneck in the LCE/LCN field. Most groups lack predictive multi-physics simulation capabilities, which is why an adopted-across-borders platform reads, in retrospect, as the distinctive infrastructure contribution from the PhD.
The platform’s real test was not whether it produced correct answers in my own hands. It was whether anyone else could pick it up from documentation alone and use it on their own problems. External groups did exactly that — working independently, without direct supervision, on problems I had not anticipated when I wrote the original code. Thirteen publications later, across work from Eindhoven to Bangalore, I think the answer to that test is clear.
Atomistics — when the continuum view broke down
The FEM platform answers questions about macroscopic shape change well. It does not answer questions about why the molecular-scale mechanism produces that shape change. By the time I reached the middle years of my PhD, the consensus explanation in the LCN actuation literature was what researchers called the static-geometry hypothesis: macroscopic deformation arises from the geometric difference between rod-like trans isomers and bent cis isomers accumulating in the cis state under illumination. The logic seemed reasonable — more cis isomers means more geometric mismatch, more stress, more deformation.
I was not convinced. Experimental results from the Broer group at TU Eindhoven had shown significant density reduction under dynamic illumination that the static hypothesis could not account for. The density change required continuous molecular motion, not static accumulation.
To test this, I built an all-atom molecular dynamics framework in LAMMPS using the PCFF class-2 force field for azobenzene-crosslinked liquid-crystal polymers, working with 7,278 atoms at 25 mol% crosslinking density at 300 K. The methodological core is a probabilistic dihedral switching protocol: a Python-orchestrated stochastic algorithm that modifies the C-N=N-C dihedral energy of each azobenzene mesogen at every timestep, modeling the continuous trans-cis-trans cycling that occurs under UV illumination. Earlier MD studies had modeled photo-isomerization by instantaneously switching isomer populations — a static operation. My protocol made switching dynamic and continuous.
The result was unambiguous. Dynamic trans-cis-trans cycling produces a 15.7% density reduction. Static isomerization — cis isomers accumulating and staying there — produces none. The Pearson correlation between isomerization cycling frequency and density change was 0.88. That disproves the static-geometry hypothesis computationally and reproduces the free-volume generation mechanism that experiments had been pointing toward. The work was published in the Journal of Chemical Physics in 2024, and the simulation code is publicly available.
What the atomistic work taught me, as much as any specific result, was when to trust a continuum model and when to go finer. The FEM platform tells you what a system does at the device scale given a constitutive model. The MD tells you whether the constitutive model’s underlying assumptions are right. Both are necessary, and neither substitutes for the other.
Inverse problems and design tools
Predictive FEM modeling of LCN actuation requires accurate kinetic parameters — specifically the trans-to-cis and cis-to-trans rate constants that govern how quickly chromophores respond to illumination. In solution or uncrosslinked mixtures, a single first-order kinetic term is enough. In crosslinked networks, the reduced free volume and increased steric interactions slow isomerization, and a two-term first-order model with an offset is required. That means fitting seven parameters per polymer composition from noisy UV-Vis absorbance data — a non-convex optimization problem where gradient-based methods get stuck.
I developed a genetic algorithm toolkit in MATLAB (roughly 750 lines) to solve this inverse problem systematically. The architecture uses binary encoding, roulette-wheel selection, adaptive mutation, and a two-stage hierarchical optimization: a broad global search followed by local neighborhood refinement. The result is a kinetic parameter database covering five polymer crosslinking densities, seven parameters each — the first systematic database of its kind for azobenzene LCN compositions.
This was published in Polymer in 2025 as a co-first author paper with D. Jayoti, who performed the UV-Vis spectroscopy. My contribution was the optimization framework and the database construction; hers was the experimental data the framework was fitted to. That division of labor — one person building measurement infrastructure, one building computational infrastructure, both necessary for the result — captures how most of the PhD collaborative work was structured.
The toolkit matters beyond the specific compositions we characterized. Any group modeling a new azobenzene LCN composition now has both the tool and the methodological template for extracting parameters from their own UV-Vis data, rather than guessing or borrowing parameters from a different material system.
What this work reaches into now
The platform is still in use. Collaborators at IIT Madras and at groups in three countries continue to draw on it for new problems I never modeled myself. That was the point — to build something with a longer life than the thesis that produced it.
During the PhD I mentored six M.S. and Ph.D. students on computational modeling projects. Two of those collaborations produced co-authored peer-reviewed publications: one in Soft Matter in 2021 on light-driven miniature locomotion, one in Soft Matter in 2025 on multi-wavelength LCN actuation. Setting up the computational infrastructure that made those projects tractable, and then handing it to students who could carry it forward, was as useful as any specific technical result.
The multi-scale pattern from the PhD — FEM at the continuum scale, MD at the atomistic scale, analytical mechanics at the structural scale, genetic algorithms for inverse parameter extraction — is not accidental. It reflects what the physics required at each step. When the continuum model was sufficient, I used it. When it was not, I went finer. When the problem was analytical, I derived. When it was a non-convex optimization, I reached for evolutionary methods.
That same instinct for following the physics to whatever method it demands is what connects this work to what I now do at Los Alamos: DFT and machine-learning interatomic potentials for electrochemical interfaces require the same discipline. Different domain, different tools, same architecture underneath — new domain, deep physics, reusable infrastructure, validated science.
The tools described here, including the FORTRAN FEM platform, the LAMMPS MD setup, the genetic algorithm kinetics toolkit, and the analytical bilayer curvature derivation, are described further on the /software page.