شبیه سازی دینامیک مولکولی – به زبان ساده

دینامیک مولکولی (MD)، روشی برای آنالیز حرکت فیزیکی اتم‌ها و مولکول‌ها است. به اتم‌ها و مولکول‌ها در مدت زمان مشخصی فرصت داده می‌شود تا با یکدیگر برهم‌کنش انجام دهند تا از طریق آن، دیدگاهی در خصوص تکامل دینامیکی سیستم بدست آید. خط سیر اتم‌ها و مولکول‌ها را به کمک حل عددی معادله حرکت نیوتون تعیین می‌کنند. از آن‌جایی که یک سیستم مولکولی به طور معمول شامل تعداد اجزای بسیار زیادی است، به طور تحلیلی نمی‌توان خواص چنین سیستم‌های پیچیده‌ای را تعیین کرد. به همین دلیل، از روش‌های عددی کمک می‌گیرند. در این آموزش، موارد مقدماتی را جهت آشنایی با دینامیک مولکولی به زبان ساده بررسی می‌کنیم.

مقدمه

اساس نظری برای دینامیک مولکولی، شامل بررسی دیدگاه‌های بسیاری از دانشمندان در خصوص مکانیک تحلیلی از جمله اویلر،‌ همیلتون، لاگرانژ و نیوتون است. ساده‌ترین شکل دینامیک مولکولی، چیزی بیشتر از قانون دوم نیوتون را در بر دارد. درحالیکه شبیه‌سازی دینامیک مولکولی به شدت به کامپیوترها وابسته است اما نگاهی هم به توسعه دو موضوع فیزیک یعنی نسبیت و مکانیک کوانتوم دارد.

فیلم آموزش ترمودینامیک مهندسی شیمی – جامع و کاربردی در فرادرس

کلیک کنید

نسبیت خاص، انتقال اطلاعات را در سرعت‌های بالاتر از نور شرح می‌دهد. مکانیک کوانتوم نیز، اصل عدم قطعیت را شامل می‌شود. در این بین، شبیه‌سازی دینامیک مولکولی، نیازمند اطلاعات کاملی در خصوص مکان و تکانه در هر لحظه است.

برهم‌کنش‌ها و معادله حرکت

ابتدایی‌ترین مدل میکروسکوپی برای مواد با سه حالت جامد، مایع و گاز، بر اساس برهم‌کنش ذرات کروی با یکدیگر بنا شده است. این برهم‌کنش‌ها در ساده‌ترین حالت خود، بین جفت‌های اتمی اتفاق می‌افتند که در نهایت، دو ویژگی اساسی را در نیروهای بین‌اتمی بوجود می‌آورند. اولین ویژگی، مقاومت در برابر متراکم شدن و دومین ویژگی، تشکیل پیوند بین اتم‌ها در حالت‌های مایع و جامد است. توابع پتانسیلی که این ویژگی‌ها را توصیف کنند، شکل‌های مختلفی خواهند داشت و اگر به دقت انتخاب شوند، مدل صحیحی از مواد را بدست می‌دهند.

از جمله بهترین مدل‌ها،‌ مدلی است که برای آرگون مایع موسوم به «لنارد جونز» (LJ) پیشنهاد شده است. انرژی پتانسیل در این رابطه برای یک جفت اتم $$i$$ و $$j$$، که در $$r_i$$ و $$r_j$$ قرار دارند به صورت زیر تعریف می‌شود:

$$\begin{equation}<br /> u\left( r _ { i j }\right) = \left\{\begin{array}{ll}<br /> {4 \epsilon\left[\left(\frac{\sigma}{r _ { i j}}\right) ^ {1 2 }-\left(\frac{\sigma}{r _ { i j}}\right) ^ { 6 }\right]} & {r_{i j}<r_{c}} \\<br /> {0} & {r _ {i j} \geq r_{c}}<br /> \end{array}\right.<br /> \end{equation}$$

در این رابطه:

• $$r _ { i j } = r_ i - r_ j$$
• $$r _ { i j } = |r _ { i j }|$$
• $$\epsilon$$: بیانگر قدرت برهم‌کنش
• $$\sigma$$: مقیاس طول

این برهم‌کنش‌ها در بازه‌های بسیار نزدیک، به صورت دافعه‌ هستند و بعد از آن به صورت جاذبه عمل می‌کنند که این جاذبه در نهایت، در $$r_c$$ از بین می‌رود. این برهم‌کنش‌ها شامل جفت‌های اتمی می‌شوند که هر جفت را به طور جداگانه مورد بررسی قرار می‌دهند به گونه‌ای که اتم‌ها موجود در همسایگی، هیچ اثری بر نیروهای بین آن‌ها ندارند. این رابطه را به صورت ساده شده زیر نیز نشان می‌دهند:

$$\begin{equation} u\left( r _ {i j }\right)=\left\{\begin{array}{ll} {4 \epsilon\left[\left(\frac{\sigma}{r _ {i j}}\right)^{1 2}-\left(\frac{\sigma}{r_{i j}}\right)^{6}\right]+\epsilon} & { r _ { i j} < r _ {c } = 2 ^ { 1 / 6} \sigma} \\ {0} & {r _ {i j} \geq r _ { c}} \end{array}\right. \end{equation}$$

مدلی که از طریق این رابطه بدست می‌آید، قدری پیچیده‌تر از تجمع ذرات توپی‌شکل و نرم و برخورد آن‌ها است. نیروی متناظر با $$u ( r )$$ عبارتست از:

$$\begin{equation} f = - \nabla u ( r) \end{equation}$$

بنابراین،‌ نیرویی که اتم $$j$$ بر اتم $$i$$ وارد می‌کند برابر است با:

$$\begin{equation} f _ { i j } = \left(\frac{ 4 8 \epsilon}{\sigma ^ { 2 }}\right)\left[\left(\frac{\sigma}{r _ {i j}}\right) ^ {1 4 } - \frac{ 1}{ 2}\left(\frac{\sigma}{ r _ { i j}}\right) ^ {8}\right] r _ { i j} \end{equation}$$

مقدار عبارت بالا برای مقادیر $$r _ {ij} < r _ c$$ و در غیر از این حالت، برابر با صفر است. هرقدر میزان $$r$$ به $$r_c$$ نزدیک‌تر (بزرگتر) شود،‌ نیرو به صفر نزدیک می‌شود. بنابراین هیج ناپیوستگی در $$r_c$$ برای نیرو و پتانسیل دیده نمی‌شود. معادله حرکت در قانون دوم نیوتون،‌ به صورت رابطه زیر تعریف می‌شود:

$$\begin{equation} m \ddot{r}_{i}=f_{i}=\sum_{j=1 \atop j \neq i}^{N_{m}} f_{i j} \end{equation}$$

در این رابطه، $$m$$ جرم اتمی است. قانون سوم نیوتون نیز بیان می‌کند که $$f _ {ij} = - f _ {ji}$$، بنابراین، هر جفت‌ اتمی را تنها یک بار مورد بررسی قرار می‌دهیم. مقدار کار انجام شده نیز با $$N _ m ^ 2$$ متناسب است، در نتیجه، برای مدل‌هایی که در آن‌ها $$r_c$$ در مقایسه با اندازه ظرف، کوچک باشند، بهتر است که جفت‌های اتمی را برای $$r {i j } \leq r _ c$$ مشخص کنیم تا مقادیر محاسبات کاهش یابند.

واحدهای بدون بعد در دینامیک مولکولی

در ادامه قصد داریم تا برخی واحدهای بدون بعد در دینامیک مولکولی را بیان کنیم. دلایل مختلفی برای این کار وجود دارند که از آن‌جمله می‌توان به کار با مقادیر عددی اشاره کرد که به جای اینکه اعداد بسیار کوچکی باشند، مقادیری نزدیک به واحد دارند. از جمله مزایای دیگر واحدهای بدون بعد، این است که معادله حرکت نیز ساده خواهد شد. با بدون بعد کردن ابعاد و ساده شدن روابط، در حقیقت، یک مدل ساده می‌تواند مسائل مختلفی را توصیف کند و زمانی که خصوصیات، به صورت بدون بعد اندازه‌گیری شدند، به راحتی می‌توان آن‌ها را به مقیاس‌های مختلف تبدیل کرد.

در مطالعات دینامیک مولکولی با استفاده از مدل LJ، مناسب‌ترین واحدهای بدون بعد با انتخاب $$\sigma$$، $$m$$ و $$\epsilon$$ به ترتیب به صورت واحدهای طول، جرم و انرژی و به شکل زیر، تعریف می‌شوند:

• طول: $$r \rightarrow r \sigma$$
• انرژی: $$e \rightarrow e \epsilon$$
• زمان:‍ $$\begin{equation} t \rightarrow t \sqrt{m \sigma^{2} / \epsilon} \end{equation}$$

در نتیجه، شکل اصلی معادله حرکت به صورت واحدهای دینامیک مولکولی، به شکل زیر خواهد بود:

$$\begin{equation} \ddot{r}_{i}=48 \sum_{j(\neq i)}\left(r_{i j}^{-14}-\frac{1}{2} r_{i j}^{-8}\right) r_{i j} \end{equation}$$

 انرژی‌های جنبشی و پتانسیل نیز به ازای هر اتم، به شکل زیر تعریف می‌شوند:

$$\begin{equation}<br /> \begin{array}{l}<br /> {E_{K}=\frac{1}{2 N_{m}} \sum_{i=1}^{N_{m}} v_{i}^{2}} \\<br /> {E_{U}=\frac{4}{N_{m}} \sum_{1 \leq i<j \leq N_{m}}\left(r_{i j}^{-12}-r_{i j}^{-6}\right)}<br /> \end{array}<br /> \end{equation}$$

در رابطه بالا، $$v_i$$ را سرعت در نظر می‌گیرند. نمودار انرژی‌های برهم‌کنش LJ و «کره‌های نرم»‌ (Soft Sphere) در واحدهای دینامیک مولکولی بر حسب فاصله را در تصویر زیر مشاهده می‌کنید.

به همین صورت، دما در یک سیستم دو یا سه بعدی نیز با رابطه زیر تعریف می‌شود:

$$\begin{equation} T=\frac{1}{d N_{m}} \sum_{i} v_{i}^{2} \end{equation}$$

اگر مدل بمنظور شبیه‌سازی دینامیک مولکولی آرگون مایع در نظر گرفته شده باشد، ارتباط بین واحدهای بدون بعد دینامیک مولکولی و داده‌های فیزیکی واقعی به صورت زیر خواهد بود:

طول به صورت $$\sigma=3.4 \mathrm{A}$$ تعریف می‌شود.

واحدهای انرژی به صورت $$\epsilon / k_{B}=120 \mathrm{K}$$ خواهند بود به گونه‌ای که: $$\epsilon = 120 \times 1.3806 \times 10 ^ {-16} erg /atom$$

اگر $$m$$ را به صورت جرم اتم آرگون و برابر با مقدار $$m=39.95 \times 1.6747 \times 10^{-24} \mathrm{g}$$ در نظر بگیریم، واحد زمانی دینامیک مولکولی برابر خواهد بود با $$2.161 \times 10 ^ {-12} s$$. در نتیجه، اگر میزان «گام زمانی» (Timestep) را برابر با $$\Delta t = 0.005$$ در نظر بگیریم، این میزان متناظر با مقدار $$10 ^ {-14} s$$ خواهد بود.

در نهایت، اگر تعداد $$N_m$$ اتم، فضایی مکعبی با ضلع L را اشغال کنند، چگالی معمول در مایع و برابر با $$0.942 g/cm^ 3$$، بیان می‌کند که $$L=4.142 N_{m}^{1 / 3} \mathrm{A}$$

شرایط مرزی

سیستم‌های «محدود» (Finite) و «نامحدود» (Infinite) با یکدیگر بسیار متفاوت هستند و میزان بزرگی یک سیستم کوچک، برای بیان رفتار یک سیستم نامحدود، پاسخی یکتا ندارد. شبیه‌سازی برای محفظه‌ای صورت می‌گیرد که دیواره‌های آن در حکم مرزی سخت هستند که اتم‌ها برای خروج از محدوده شبیه‌سازی، به آن‌ها برخورد می‌کنند. در سیستم‌هایی با اندازه ماکروسکوپی، تنها کسر کوچکی از اتم‌ها به میزان کافی به دیواره‌ها نزدیک هستند تا انحراف از فضای قالب را تجربه کنند.

به طور مثال، سیستمی سه بعدی با چگالی مایع و $$N _ m = 10 ^ {21}$$ را در نظر بگیرید. از آن‌جایی که تعداد اتم‌های نزدیک دیواره، از مرتبه $$N _ m ^ {2/3}$$ است، مقدار آن‌ها برابر با $$10 ^ {14}$$ خواهد بود. اما به طور معمول، برای مقادیری از $$N_ m = 1000$$، در حدود 500 اتم به دیواره‌ها نزدیک هستند. به همین دلیل، تنها زمانی این شبیه‌سازی موفق خواهد بود که بخواهیم رفتار اتم‌ها را در نزدیکی دیواره بررسی کنیم.

می‌توانیم یک سیستم محدود اما فارغ از دیواره فیزیکی را با استفاده از «شرایط مرزی تناوبی» (Periodic Boundary Condition) ایجاد کنیم که آن را در تصویر زیر مشاهده می‌کنید. این شرایط مرزی تناوبی را می‌توان به صورت آرایه‌ای از نمونه‌های مشابه محدوده شبیه‌سازی معرفی کرد. این شرایط مرزی تناوبی، دو تاثیر دارد: اول این‌که اتمی که ناحیه شبیه‌سازی را ترک کند، به سرعت به ناحیه دیگر وارد می‌شود. دومین مورد این است که اتمی که در فاصله $$r_c$$ از مرز قرار دارد، با اتم‌ها در ناحیه مجاور برهم‌کنش انجام می‌دهد که به این حالت «اثر کمربندی» (Wraparound Effect) می‌گویند.

اثر «Wraparound» را باید در معادلات حرکت و محاسبات برهم‌کنش‌ها لحاظ کرد. بعد از انجام هر «گام محاسباتی» (Integration Step) باید مختصات به طور مجدد ارزیابی شوند و اگر اتمی به خارج از محدوده حرکت کرده باشد، مختصات آن به طور مجدد تنظیم شوند تا دوباره در محدوده قرار بگیرد. به طور مثال، اگر مختصات $$x$$ را طوری تنظیم کرده باشیم که بین $$- L _ x /2$$ و $$L _ x /2$$ قرار بگیرد - که در این جا $$L_x$$، اندازه محدوده در جهت $$x$$ است - تنظیمات را می‌توانیم به صورت زیر اعمال کنیم. در صورتیکه محدوده مورد نظر به صورت مستطیلی دو بعدی یا منشور در سه بعد باشد،‌ کار با مرزهای تناوبی بسیار ساده خواهد بود:

• اگر $$\begin{equation} r _ {i x} \geq L _ { x } / 2 \end{equation}$$، آن را با $$r _ {i x } - L_ x$$ جایگزین کنید.
• در غیر اینصورت، اگر $$\begin{equation} r _ {i x} < -L _ { x } / 2 \end{equation}$$، آن را با $$\begin{equation} r _ {i x} + L _ { x } \end{equation}$$ جایگزین کنید.

برای لحاظ کردن برهم‌کنش‌ها در محاسبات نیز، به طور مشابه، به صورت زیر عمل می‌کنیم:

• اگر $$\begin {equation} r _ { i j x } \geq L _ { x } / 2 \end {equation}$$، آن را با $$\begin {equation} r _ { i j x }- L _ { x } \end {equation}$$ جایگزین کنید.
• در غیر اینصورت،‌ اگر $$\begin{equation} r _ {i j x} < -L _ { x } / 2 \end{equation}$$،  آن‌را با $$\begin{equation} r _ {i j x} + L _ { x } \end{equation}$$ جایگزین کنید.

شرایط اولیه