GROMACS version: 2022
GROMACS modification: No
Here post your question
Hi All,
I am conducting a deformation MD simulation for a polymer composite using the GAFF force field with a stretching rate of 1 × 10⁻³ nm ps⁻¹. Experimentally, this material becomes stretchable in the presence of PEGs.
However, in my deformation MD simulations, Young’s modulus values for the polymer composite with and without PEGs are almost identical, which contradicts experimental results. Experimentally, Young’s modulus for the polymer composite with PEG is below 100 MPa, but my simulations yield around 6000 MPa for both cases. I also tried different stretching rates, such as 1 × 10⁻² and 1 × 10⁻⁴ nm ps⁻¹, but none of these produced any significant difference.
Can anyone suggest methods or approaches to observe a difference in Young’s modulus for the polymer composite with and without PEGs in MD simulations?
I highly appreciate your recommendations.
@MagnusL I was wondering if you could provide me with guidance on this issue. Apologies for asking so many questions about deformation simulations. I am using GROMACS 2022 for my simulations. I attempted to use the 2024 version, but the simulation failed at a stretching rate of 10-4 nm ps-1. At a rate of 10-3 nm ps-1 I observed significant fluctuations.
All transport coefficients are very potential-dependent. I am afraid I am not familar with GAFF: @Ema2022, do you expect GAFF to provide results that are in agreement with experiments in the first place?
I also tried different stretching rates, such as 1 × 10⁻² and 1 × 10⁻⁴ nm ps⁻¹, but none of these produced any significant difference.
Makes sense to me, given that the elatic behaviour doesn’t depend on the strain rate but only on the deformation. Would you expect something different? What are you relating to the uniaxial stress to compute the Young modulus? Is it the deformation velocity or the deformation?
I attempted to use the 2024 version, but the simulation failed at a stretching rate of 10-4 nm ps-1. At a rate of 10-3 nm ps-1 I observed significant fluctuations.
This is a bit odd: I would expect the simulation to fail for higher deformation velocities, not lower. This makes me wonder whether the initial condition is stable.
PS It is not a deformation rate, but a deformation velocity (this may seem nit-picking, but it makes all the difference in rheology).
Moreover, I am assuming you are deforming the simulation box in one direction: what are you doing in the other two? Are you applying another deformation to keep the volume constant? Are you applying a semiisotropic or anisotropic barostat?
Thank you for replying to my questions. The GAFF forcefield has been validated for my polymer composite, and simulations using this forcefield successfully reproduce the average lamella crystallite size, characteristic π–π stacking distance, and XRD patterns consistent with experimental results. Initial structures, including the polymer composite alone and the polymer composite with PEG, were thoroughly equilibrated using multiple annealing cycles. Equilibration was confirmed by calculating various properties to ensure system stability.
For the deformation simulations, I tested deformation velocities of 1 × 10⁻² and 1 × 10⁻⁴ nm ps⁻¹, motivated by discussions about potential artifacts arising from the stretching rate. Stress was computed as the negative of the pressure tensor component in the deformation direction
(𝜎 = −𝑃𝑧), while the engineering strain (ϵ) was defined as ϵ = 𝐿𝑡−𝐿0/𝐿0, where Lt is the length of the simulation box in the extension direction at time t, and L0 is the initial box length. I extracted the pressure tensor and box length over time using gmx_energy and calculated 𝜎 and ϵ. Young’s modulus (𝐸) was determined as the slope of the initial linear region of the stress-strain curve, specifically within the strain range of 0.3% to 3%.
Blockquote
Moreover, I am assuming you are deforming the simulation box in one direction: what are you doing in the other two?
Yes, currently, I am deforming the simulation box in the z-direction only. However, I plan to test deformations in the other two directions for both systems: the polymer composite and the polymer composite with PEG.
Blockquote
Are you applying another deformation to keep the volume constant?
I do not understand this. Could you please elaborate?
Blockquote
Are you applying a semiisotropic or anisotropic barostat?
I am currently using a semi-isotropic barostat.
One more point: I attempted z-direction deformation using a deformation velocity of 10−4 nm ps−1 with the 2024 version, this time enabling deform-init-flow (deform-init-flow = yes). So far, the simulation is running without any failures. Previously, I had this option turned off for all simulations. Could this be the reason for the observed issue, as velocities were not updated during deformation, potentially causing Young’s modulus to appear similar in both systems? However, this option is not available in GROMACS 2022.
I attempted z-direction deformation using a deformation velocity of 10−4 nm ps−1 with the 2024 version, this time enabling deform-init-flow (deform-init-flow = yes ). So far, the simulation is running without any failures.
Ok, this explains a lot: there has been a change of behaviour of the deformation code from 2023 to 2024 (this is why that option was not present in 2022). I can imagine that not initializing the flow would cause artifacts. In essence, velocities in the deformation direction are actually never updated (apart from thermalization, of course), that’s why it is important to initialize atomic velocities correctly. How are the results now?
How do we keep the volume constant in deformation simulation? @MichelePellegrino
I think that using a semiisotropic barostat is fine. Another option would be to shrink the system in the perpendicular directions while it is stretched in the other, to mantain the volume constant. I guess this could also make sense since many results in continuum mechanics are derived assuming no local volume change.
One last thing: I am no polymer expert, but I know that many polymers (especially long and massive ones) have a rather long relaxation time, i.e. it would be difficult to see them unfolding/stretching in a simulation that only lasts an handful of nanoseconds. So it could be that you are observing no difference between the systems with and without PEG simply because polymers don’t have time to change their conformation. This may imply that you need to simulate at even lower deformation speeds to see something happening…
Unfortunately, the deformation simulation failed again, unable to reach 100% strain. I feel the same perspective on polymers as you mentioned in your last comment. So, I plan to try a slower deformation over a longer time to see if it yields any differences. Once again, thank you so much for all your valuable insights.
I conducted tensile MD simulations with a deformation velocity of 1×10−6 nm ps−1 and obtained a Young’s modulus of 5 GPa for the polymer composite. Experimentally, the modulus ranges between 0.8–2.4 GPa. I believe this is a reasonable prediction, considering slightly higher values have been reported in the literature.
However, when I introduced PEG into the system, I could not achieve the experimentally reported Young’s modulus of ~100 MPa. My simulations for the polymer composite with PEG, using the same deformation velocity, resulted in a modulus of 4 GPa. At elevated temperatures, close to the glass transition temperature, the modulus decreased to 0.9 MPa for the composite with PEG (with high temp, Young’s modulus decreased).
Do you think this could be due to the PEG parameters, which I generated directly through the PolyParGen site? I would appreciate your advice on this.