Hello everyone,
I am looking for a critical assessment of the methodology used in my undergraduate research project. I am particularly interested in whether the coupling between Monte Carlo radiation transport and classical molecular dynamics is physically justified, and what should be improved before extending the work into a publication.
The system is hydrated Na-montmorillonite, intended as a simplified model of a clay material relevant to radioactive-waste barrier systems.
Current workflow
My workflow is approximately:
1. Build the montmorillonite structure
I constructed a Na-montmorillonite model with two layers of interlayer water. The model includes isomorphic substitutions in the tetrahedral and octahedral sheets and compensating interlayer cations.
The production model used a (6x6x3) supercell with dimensions of approximately (31x54x81) Å.
For the molecular dynamics part I use ClayFF-type parameters for the clay and SPC/E-based parameters for water.
2. Simulate γ irradiation
I consider 1.3 MeV γ-rays.
The first step is a Geant4 calculation of the energy deposition in the model. I considered the Compton interaction and the corresponding maximum kinetic energy of the recoil electron, which is about 1.086 MeV for a 1.3 MeV photon.
Because my objective was to model the local consequences of radiation energy deposition rather than the macroscopic photon transport itself, I subsequently used electrons in the corresponding energy range to represent the secondary electrons produced in the material.
The absorbed dose in the model was approximately 0.3 MGy.
3. Introduce radiation-chemical products
The deposited energy was used to estimate the number of ionization/radiolysis events. For the particular simulation I obtained up to approximately 18 radiation-chemical species in one supercell.
These species were then explicitly introduced into the atomistic model according to the assumed early-time products of water radiolysis.
4. Molecular dynamics
The radiation-modified structure is then relaxed in LAMMPS.
I treat this relaxation as a model of the post-radiation structural relaxation of the system. After relaxation, I calculate diffusion of interlayer water and exchange cations.
The main transport quantity is the self-diffusion coefficient obtained from the Einstein relation:
I also considered the Green-Kubo/VACF formulation as an independent way of evaluating diffusion.
In addition, I analyze RDF, MSD, SASA and the interlayer spacing.
What I would particularly like feedback on
There are several methodological points where I am not confident that my approximations are sufficiently rigorous:
- Photon → electron treatment
Is replacing the 1.3 MeV γ field by electrons in the 0–1.086 MeV range a defensible approximation for this purpose?
I suspect that the important issue is not simply the energy range, but the actual secondary-electron energy spectrum and its weighting. Should I explicitly sample the Compton electron spectrum generated by Geant4 instead of using a simple energy range?
- Conversion from deposited energy to radiolysis products
Is it acceptable to estimate the number of radiolysis events from the deposited energy / mean energy per ion pair and then explicitly introduce the corresponding chemical species into a classical MD model?
Or would this require a more explicit radiation-chemistry model (for example Geant4-DNA chemistry, another radiolysis model, or experimental G-values)?
- Force-field consistency
After irradiation, the chemical composition of the system is no longer identical to the initial montmorillonite + water system.
I therefore introduce parameters for OH, HO₂, H₂O₂ and related species and perform classical MD with ClayFF/SPC/E-derived parameters.
Is this a reasonable approximation for studying the subsequent structural relaxation and diffusion, or does the presence of newly formed chemical species make a reactive force field essential?
- Interpretation of the resulting diffusion coefficients
My main hypothesis is that radiation-induced chemistry modifies the interlayer environment and therefore changes the diffusion of water and exchange cations.
Would you consider the comparison of MSD-derived diffusion coefficients between irradiated and non-irradiated systems sufficient evidence for such a conclusion, assuming that the trajectories are long enough and statistically converged?
- Validation
What would you consider the minimum set of validation calculations for such a workflow?
For example:
-
convergence with respect to supercell size;
-
multiple independent irradiation realizations;
-
multiple independent MD trajectories;
-
dependence on absorbed dose;
-
comparison with experimental diffusion/interlayer-spacing data;
-
comparison with an alternative force field;
-
independent evaluation of diffusion using Green-Kubo/VACF.
I am especially interested in criticism of the physical validity of the multiscale coupling itself, rather than only technical LAMMPS issues.
The project is currently a proof-of-concept study, so I am fully aware that several assumptions are simplified. My goal now is to determine which assumptions are acceptable for a qualitative study and which ones must be changed before moving toward a quantitative publication.
Any criticism of the methodology is very welcome.
Thank you.
P.S. Regarding the force field choice: to the best of my knowledge, there is currently no verified ReaxFF parameterization optimized for modeling complex clay-water-radiolysis products simultaneously. To bypass this limitation in our classical MD setup, I utilized SASA (Solvent-Accessible Surface Area) analysis. By evaluating SASA, I established a geometric proxy for the probability of radical-pore and radical-ion interactions, treating the spatial availability of the interlayer environment as a key factor in the post-irradiation behavior.