حول أهمية تفاصيل التكامل العددي لمحاكاة التدفق المتجانس On the importance of numerical integration details for homogeneous flow simulation

شارك:
المجلة: The Journal of Chemical Physics، المجلد: 164، العدد: 8
DOI: https://doi.org/10.1063/5.0315430
PMID: https://pubmed.ncbi.nlm.nih.gov/41744469
تاريخ النشر: 2026-02-26
المؤلف: Zhenyun Du وآخرون
الموضوع الرئيسي: التفاعلات الجزيئية الفيزيائية والكيميائية المتقدمة

نظرة عامة

تقدم البحث نظام تكامل عددي قابل للعكس وموفر للطاقة لمعادلات الحركة Sllod، والتي تعتبر أساسية لنمذجة التدفق المتجانس على المستوى الذري وتوقع خصائص السوائل مثل اللزوجة. يحقق النظام المقترح خطأً من رتبة $\delta t^3$، مقارنةً بمجزئات تقسيم المشغل القياسية المستخدمة في محاكاة الديناميكا الجزيئية. تم تنفيذ هذه الطريقة في برنامج LAMMPS، مما يعزز دقة المحاكاة التي تتضمن استجابات عابرة، وتدفقات مختلطة، وحالات مستقرة، خاصةً تحت معدلات قص عالية. ومن الجدير بالذكر أن الدراسة تسلط الضوء على أن الحفاظ على الطاقة أمر حاسم، حيث يمكن أن يؤدي غيابه إلى أخطاء منهجية في المتوسط الجماعي لمؤشر الضغط، مما يؤثر بشكل كبير على حسابات اللزوجة عند معدلات تدفق مرتفعة.

في الختام، يوضح المؤلفون أن نظام التكامل القابل للعكس الخاص بهم لديناميات Sllod يحافظ على الطاقة بالنسبة لمؤشرات التدفق المثلثية العامة. لا يحسن التنفيذ في LAMMPS دقة حسابات اللزوجة عند معدلات القص العالية فحسب، بل يتماشى أيضًا مع النتائج المتوقعة من نظرية الاستجابة الدقيقة، مما يعالج التناقضات المبلغ عنها سابقًا. تؤكد النتائج على أهمية الحفاظ على الطاقة في محاكاة الديناميكا الجزيئية، خاصةً في سياق تدفقات السوائل عالية المعدل.

مقدمة

تناقش مقدمة الورقة تطبيق محاكاة الديناميكا الجزيئية غير المتوازنة (NEMD)، وخاصةً باستخدام معادلات الحركة Sllod، للتحقيق في الخصائص الريولوجية للسوائل. تسهل معادلات Sllod نمذجة التدفق اللامتناهي المتجانس داخل خلية وحدة دورية لانهائية، مما يلغي الحاجة إلى نمذجة الجدران. يتم تقديم المعادلات من حيث مواقع الجسيمات، والزخم الغريب، وتدرج ملف سرعة التدفق، مع تمييز فقدان الطاقة من خلال صياغة محددة تتماشى مع التوقعات الهيدروديناميكية.

يسلط المؤلفون الضوء على الاستخدام الواسع لمعادلات Sllod في حزمة الديناميكا الجزيئية LAMMPS، مشيرين إلى أنه بينما يمكن لرموز أخرى التعامل مع تشوه خلية الوحدة الدورية، فإنها غالبًا ما تفعل ذلك بطريقة أقل منهجية. تهدف المخطوطة إلى معالجة التحديات المتعلقة بالتكامل العددي المرتبطة بمعادلات Sllod، خاصةً في الاستجابات العابرة ومعدلات التدفق العالية، وتقدم تنفيذًا جديدًا في LAMMPS يعزز الدقة. بالإضافة إلى ذلك، يستنتج المؤلفون كمية محفوظة لديناميات Sllod ذات التحكم الحراري تحت ظروف تدفق عشوائية، والتي تعمل كاختبار استقرار لمتكاملهم العددي. توسع هذه الدراسة على الدراسات السابقة من خلال النظر في مؤشرات التدفق المثلثية العشوائية، متجاوزة التركيز التقليدي على تدفق القص المستوي.

النتائج

في هذا القسم، يقدم المؤلفون نتائج التعديلات التي أجريت على برنامج محاكاة LAMMPS لضمان الحفاظ على الكمية \( H’ \) خلال محاكاة الديناميكا الجزيئية. تضمنت التعديلات الرئيسية محاذاة خطوات التكامل مع المعادلات المستنتجة، وتمكين تشوه خلية الوحدة في كل خطوة زمنية، وتصحيح حساب التدفقات المختلطة. يسمح التنفيذ المعدل بإجراء محاكاة طويلة الأمد بكفاءة وإعادة تعيين دقيقة للسرعات في كل من الإطارات المختبرية والغريبة. أجرى المؤلفون محاكاة واسعة النطاق، مولدين 40,000 حالة أولية من 50 نظامًا مستقلًا، وطبقوا معادلات الحركة Sllod لتقييم أداء تعديلاتهم.

تظهر النتائج أن كل من تنفيذات الإطار الغريب وإطار المختبر تحافظ بفعالية على \( H’ \) تحت تدفقات القص المختلفة، على عكس التنفيذ الأصلي لـ LAMMPS، الذي أظهر أخطاء منهجية كبيرة. على وجه التحديد، تم التحقق من الحفاظ على \( H’ \) من خلال مقارنات حسابات لزوجة القص، مما يكشف أن الأنظمة المعدلة أنتجت نتائج متسقة مع صياغة دالة الارتباط الزمني العابر (TTCF)، بينما أدى التنفيذ الأصلي إلى تقديرات مفرطة للزوجة بسبب عدم قدرته على الحفاظ على الطاقة بشكل صحيح. علاوة على ذلك، تشير النتائج إلى أن معدل تغير \( H’ \) يتناسب مع \( O(\delta t^2) \)، متسقًا مع التوقعات النظرية، مما يؤكد قوة نظام التكامل عبر تعقيدات التدفق المختلفة.

المناقشة

في هذا القسم، يقدم المؤلفون كمية محفوظة لمعادلات الحركة Sllod، والتي عادةً ما تكون غير هاملتونية. من خلال إدخال متغير طور إضافي، يستنتجون تعبيرًا أبسط للكمية المحفوظة التي تتضمن الطاقة الحركية الحرارية وتأخذ في الاعتبار ظروف الحدود المتغيرة. يتم التعبير عن الكمية المحفوظة، التي يشار إليها بـ \( H’ \)، على النحو التالي:

\[
H’ = \sum_i \frac{p_i \cdot p_i}{2m_i} + \Phi(q, a, b, c) + \frac{1}{2} Q p_\eta^2 + N_d k_B T \eta + \lambda,
\]

حيث تمثل \( \lambda \) مصدر طاقة يدفع التدفق. يثبت المؤلفون أنه لكي يتم الحفاظ على \( H’ \)، يجب أن يفي تطور \( \lambda \) بـ \( \dot{H}’ = 0 \)، مما يؤدي إلى مجموعة من معادلات الحركة التي تضمن أن أنظمة التكامل العددي يمكن أن تحافظ على \( H’ \) على الرغم من الطبيعة غير التبادلية لمعادلات Sllod.

تستكشف المناقشة أيضًا آثار الكمية المحفوظة في التدفقات الثابتة، مشيرةً إلى أنه بينما تتأرجح \( \langle \kappa(t) \rangle \) حول الصفر بسبب الهياكل الشبكية الدورية، تظل توزيع فضاء الطور غير متوازن. يؤكد المؤلفون على أهمية دمج معادلات الحركة بدقة وتحديث خلية الوحدة بالتزامن مع مواقع الجسيمات لتجنب الأخطاء المنهجية في حسابات الضغط. كما يسلطون الضوء على التحديات المتعلقة بتنفيذ هذه المعادلات في محاكاة الديناميكا الجزيئية على نطاق واسع، خاصةً فيما يتعلق بمعالجة ظروف الحدود وتكامل متجهات الشبكة.

Journal: The Journal of Chemical Physics, Volume: 164, Issue: 8
DOI: https://doi.org/10.1063/5.0315430
PMID: https://pubmed.ncbi.nlm.nih.gov/41744469
Publication Date: 2026-02-26
Author(s): Zhenyun Du et al.
Primary Topic: Advanced Physical and Chemical Molecular Interactions

Overview

The research presents a reversible and energy-conserving numerical integration scheme for the Sllod equations of motion, which are essential for modeling homogeneous flow at the atomic scale and predicting fluid properties like viscosity. The proposed scheme achieves an error on the order of $\delta t^3$, comparable to standard operator splitting integrators used in molecular dynamics simulations. Implemented in the LAMMPS software, this integration method enhances the accuracy of simulations involving transient responses, mixed flows, and steady states, particularly under high shear rates. Notably, the study highlights that energy conservation is critical, as its absence can lead to systematic errors in the pressure tensor’s ensemble average, significantly affecting viscosity calculations at elevated flow rates.

In conclusion, the authors demonstrate that their reversible integration scheme for Sllod dynamics is energy-conserving for general triangular flow tensors. The implementation in LAMMPS not only improves the accuracy of viscosity calculations at high shear rates but also aligns these results with those predicted by exact response theory, addressing previously reported discrepancies. The findings underscore the importance of energy conservation in molecular dynamics simulations, particularly in the context of high-rate fluid flows.

Introduction

The introduction of the paper discusses the application of Nonequilibrium Molecular Dynamics (NEMD) simulations, particularly utilizing the Sllod equations of motion, to investigate the rheological properties of fluids. The Sllod equations facilitate the modeling of homogeneous laminar flow within an infinitely periodic unit cell, eliminating the need for wall modeling. The equations are presented in terms of particle positions, peculiar momenta, and the gradient of the flow velocity profile, with energy dissipation characterized by a specific formulation that aligns with hydrodynamic expectations.

The authors highlight the widespread use of the Sllod equations in the LAMMPS molecular dynamics package, noting that while other codes can handle periodic unit cell deformation, they often do so in a less systematic manner. The manuscript aims to address numerical integration challenges associated with the Sllod equations, particularly in transient responses and high flow rates, and introduces a new implementation in LAMMPS that enhances accuracy. Additionally, the authors derive a conserved quantity for thermostatted Sllod dynamics under arbitrary flow conditions, which serves as a stability test for their numerical integrator. This work expands upon previous studies by considering arbitrary triangular flow tensors, moving beyond the traditional focus on planar shear flow.

Results

In this section, the authors present the results of modifications made to the LAMMPS simulation software to ensure the conservation of the quantity \( H’ \) during molecular dynamics simulations. Key adjustments included aligning integration steps with the derived equations, enabling unit cell deformation at each time step, and correcting the calculation of mixed flows. The modified implementation allows for efficient long-running simulations and accurate remapping of velocities in both lab and peculiar frames. The authors conducted extensive simulations, generating 40,000 initial conditions from 50 independent systems, and applied the Sllod equations of motion to assess the performance of their modifications.

The findings demonstrate that both the peculiar-frame and lab-frame implementations effectively conserve \( H’ \) under various shear flows, unlike the original LAMMPS implementation, which exhibited significant systematic errors. Specifically, the conservation of \( H’ \) was validated through comparisons of shear viscosity calculations, revealing that the modified schemes produced consistent results with the Transient Time Correlation Function (TTCF) formalism, while the original implementation led to overestimations of shear viscosity due to its inability to conserve energy properly. Furthermore, the results indicate that the rate of change of \( H’ \) scales as \( O(\delta t^2) \), consistent with theoretical expectations, confirming the robustness of the integration scheme across different flow complexities.

Discussion

In this section, the authors present a conserved quantity for the Sllod equations of motion, which are typically non-Hamiltonian. By introducing an additional phase variable, they derive a simpler expression for the conserved quantity that incorporates thermal kinetic energy and accounts for deforming boundary conditions. The conserved quantity, denoted as \( H’ \), is expressed as:

\[
H’ = \sum_i \frac{p_i \cdot p_i}{2m_i} + \Phi(q, a, b, c) + \frac{1}{2} Q p_\eta^2 + N_d k_B T \eta + \lambda,
\]

where \( \lambda \) represents an energy source driving the flow. The authors establish that for \( H’ \) to be conserved, the evolution of \( \lambda \) must satisfy \( \dot{H}’ = 0 \), leading to a set of equations of motion that ensure numerical integration schemes can preserve \( H’ \) despite the non-symplectic nature of the Sllod equations.

The discussion further explores the implications of the conserved quantity in steady-state flows, noting that while \( \langle \kappa(t) \rangle \) oscillates around zero due to periodic lattice structures, the phase space distribution remains non-equilibrium. The authors emphasize the importance of accurately integrating the equations of motion and updating the unit cell in tandem with particle positions to avoid systematic errors in pressure calculations. They also highlight the challenges of implementing these equations in large-scale molecular dynamics simulations, particularly regarding the treatment of boundary conditions and the integration of lattice vectors.