مبانی روش دینامیک مولکولی

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

معادلات حرکت نیوتن

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

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

 اگر جرم ذره ی i با mi و نیروی عمل کننده بر آن توسط سایر ذرات (و میدان خارجی در صورت وجود) fi لحاظ شود، حرکت این ذره طبق معادله نیوتن بصورت معادله‌ی (‏1) قابل محاسبه است. مدل برهمکنش بین‌مولکولی به یک تابع پتانسیل (r)U که شکل هندسی مولکول را توصیف می‌کند نیازمند است. r نیز نشان‌دهنده‌ی بردار موقعیت هر ذره نسبت به مبدا مختصات بوده که با مشخص شدن آن، توزیع کلی ذرات در سیستم مشخص می‌شود. معادله‌ی (‏2) برای تمامی ذرات موجود در سیستم بصورت پیوسته برای گام زمانی بسیار کوچک t∆ حل می‌شود. تحت این شرایط شبیه‌سازی برای مدت زمانی ادامه خواهد داشت که در انتهای این فاصله‌ی زمانی، انرژی پتانسیل سیستم مقدار مشخص قبلی خود را حفظ نماید (سیستم به تعادل برسد). مختصات اتم‌ها در بازه‌های مشخص زمانی بنا بر انتخاب کاربر در یک فایل خروجی ذخیره خواهد شد. مختصات اتم‌ها بعنوان تابعی از زمان، مسیر حرکت سیستم را نمایش می‌دهد. با متوسط‌گیری روی مسیر حرکت که در این فایل خروجی ذخیره شده است، بسیاری از خواص ماکروسکوپی طبق رابطه‌ی (‏3) قابل محاسبه خواهد بود.

معادلات و پارامترهای وابسته که برای تعیین انرژی پتانسیل اتم‌های سیستم به‌کار می‌روند در مجموع میدان نیرو گفته می‌شود. در طی سال‌های اخیر تعدادی میدان نیرو به منظور استفاده در شرایط کارکردی مختلف و برای نوع اتم متفاوت در حوزه‌های گسترده‌ای از علوم منتشر شده‌اند. به عنوان نمونه می‌توان به میدان نیروی چارم CHARMM، میدان نیروی OPLS-AA، میدان نیروی MMFF، میدان نیروی AMBER و … اشاره نمود. به منظور استفاده‌ی درست و قدرتمند از روش شبیه‌سازی دینامیک مولکولی، انتخاب میدان نیرویی که بصورت تجربی برای سیستم مورد نظر اعتبار سنجی شده باشد بسیار ضروری است.

بروزرسانی موقعیت اتم ها

اگر سیستم مولکولی متشکل از n ذره باشد، n دسته معادله‌ی مشابه وجود خواهد داشت. معادله‌ی (‏1)  -در بخش قبل- برای حل دستگاه n معادله‌ای با پردازشگر مناسب نبوده و بایستی آن‌را بصورت جبری بازنویسی کرد. به‌همین دلیل جمله‌ی مشتق مرتبه دوم معادله‌ی (‏1) با استفاده از بسط سری تیلور بصورت رابطه‌ی 4 بازنویسی خواهد شد.

از حاصل جمع رابطه‌های (‏4) و (‏5)، جمله‌ی (‏6) با دقت مرتبه‌ی دوم بدست می‌آید و از ترکیب آن با معادله‌ی (‏1)، می‌توان به رابطه‌ی (‏7) رسید. سرعت را می‌توان از طرح اختلاف مرکزی بصورت رابطه‌ی (8) نوشت.

اغلب پکیج‌های شبیه‌سازی از الگوریتمی بنام جهش قورباغه‌ای ورله (یا حالت‌های اصلاح شده‌ی آن) برای انتگرال‌گیری معادله‌ی حرکت (1) استفاده می‌کنند که کمی پیچیده بنظر می‌رسد. الگوریتم جهش قورباغه‌ای از موقعیت اتم r و سرعت آن V در زمان t استفاده کرده و با استفاده از تعریف نیرو -رابطه‌ی (‏2) در بخش قبل- سرعت در نیمه‌ی گام زمانی بعدی (t+∆t/2) را بدست آورده و در ادامه، موقعیت اتم در گام زمانی بعدی (t+∆t) حاصل می‌شود.

در انتهای یک چرخه از این الگوریتم، سرعت در گام زمانی بعدی محاسبه می‌گردد. جایگاه محاسبه‌ی موقعیت و سرعت اتم‌ها در یک چرخه از الگوریتم در شکل ‏ نشان داده شده است. در این شکل n بیانگر گام زمانی، r ،V و a به ترتیب به موقعیت، سرعت و شتاب اتم اشاره دارد. انرژی پتانسیل و نیرو نیز با U و F نمایش داده شده است.

انرژی پتانسیل

اساسی‌ترین موضوع شبیه سازی مولکولی، بحث انرژی است. انرژی کل E مجموع انرژی جنبشی و پتانسیل است. انرژی پتانسیل مجموعه‌ی مولکول‌ها نیز به انرژی پتانسیل داخلی (ناشی از اندرکنش اتم‌های به‌هم پیوسته) و خارجی (ناشی از اندرکنش اتم‌هایی که به‌هم متصل نیستند) دسته بندی می‌شود (رابطه‌ی 12). در شبیه‌سازی دینامیک مولکولی هر مولکول (شامل چندین اتم) بصورت یکسری از نقاط (کره‌ی سخت) دارای جرم و موقعیت مشخص در نظر گرفته می‌شود که توسط فنرهایی (پیوند) بهم متصل می‌شوند. برای تعریف مقادیر طول پیوند، زاویه و پیچش و همچنین اندرکنش‌های واندرولسی، الکتروستاتیکی و … بین اتم‌ها و مولکول‌ها از عبارت میدان نیرو استفاده می‌شود. میدان نیرو شامل مجموعه‌ای از معادله‌ها و ثابت‌هایی است که برای بازتولید هندسه‌ی مولکولی و خواص موردنظر ساختارهای شیمیایی به‌کار می‌رود. نیروی ناشی از گرادیان این انرژی پتانسیل طبق رابطه‌ی (13) بدست می‌آید. انرژی ناشی از اندرکنش‌های پیوندی و غیرپیوندی در ادامه توضیح داده می‌شود.

انرژی‌های غیرپیوندی

انرژی ناشی از اندرکنش‌های غیر پیوندی بر اساس تعریف پتانسیل لنارد-جونز ( و یا باکینگهام و …) و پتانسیل کولمبیک (و یا کولمبیک اصلاح شده و …) محاسبه می‌شود. این نوع انرژی اندرکنشی بر اساس فهرست همسایگی (تعداد اتم‌های در حال اندرکنش غیرپیوندی که در شعاع مشخصی از هم قرار دارند) منظور می‌شود. در رابطه‌ی (14)، σij و ϵij پارامترهای پتانسیل لنارد جونز بوده و rij نیز فاصله‌ی بین دو اتم i و j است. معمولا پارامترهای مربوط به پتانسیل لنارد-جونز بصورت روش‌های تجربی یا نیمه تجربی تعیین می‌شود. جمله‌ی از مرتبه12 رابطه‌ی (14) براساس قانون طرد پائولی ناشی از هم‌پوشانی اوربیتال‌های الکترونی بوده و جمله‌ی از مرتبه6، بیانگر نیروی پراکندگی واندروالسی است. برای مثال، ULJ و نیروی ناشی از این پتانسیل، FLJ (منفی گرادیان ULJ) بین اتم‌های آب مدل F3C در شکل زیر نشان داده شده است. این مدل توسط لویت و همکارانش برای شبیه‌سازی سیستم‌های بیولوژیکی در سال 1996 ارائه شده است.

پارامترهای پتانسیل لنارد-جونز 𝜎𝑖𝑗 و 𝜖𝑖𝑗 بین نوع‌های مختلف اتم‌ها (atom type) با استفاده از قاعده‌ی ترکیب Lorentz-Berthelot محاسبه می‌شود (رابطه‌های 15 و 16). شکل بعد نشان دهنده‌ی مقایسه‌ی انرژی پتانسیل و نیرو برای دو نوع atom type مختلف در ارتباط با هم است (CTL3 و CTL3 بعنوان نمونه)؛

اندرکنش کولمبیک Uc بین دو اتم باردار بصورت رابطه­‌ی (‏17) تعریف می­‌شود. qi و qj بار مربوط به هر اتم بوده و f=۱۴πεo است. همچنین εr نیز ثابت دی­ الکتریک نسبی است که از قبل مشخص است.

اندرکنش‌­های پیوندی

• واضح است که هر اتم تحت تاثیر تابع انرژی پتانسیل اتم ­های همسایه­‌ی خود قرار می­ گیرد. انرژی پیوندی بصورت حاصل جمع انرژی‌­های ناشی از تغییر طول پیوند، زاویه و چرخش در یک مولکول منظور می­‌شود (رابطه‌ی 18). در ادامه طرح شماتیک تمام چیدمان­‌های اتم­‌ها در کنار هم نشان داده شده اس

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

انرژی جنبشی

انرژی پتانسیل مجموع جملات مختلف ناشی از برهم­کنش­‌های پیوندی و غیرپیوندی است که در بخش‌ انرژی پتانسیل روش محاسبه­‌ی آن ذکر شده است. انرژی جنبشی کلی K برای n ذره در سیستم بصورت رابطه­‌ی (‏19) قابل بیان است. با تعریف دمای مطلق، می‌­توان انرژی و دما را با استفاده از رابطه‌ی 20 به‌­هم مرتبط نمود که در این رابطه، ثابت بولتزمن و Nf میزان درجه آزادی قابل شمارش است. 

هم بستگی دمایی

ترموستات در واقع الگوریتمی است که روال شبیه‌­سازی دینامیک مولکولی بر مبنای قانون نیوتن را اصلاح می­‌کند. سه ترموستات کاربردی در پکیج­‌های شبیه‌سازی مولکولی عبارتند از ترموستات برندسون (Berendsen) که روش ساده‌­ای دارد اما از دقت بالایی برخوردار نیست. ترموستات نوز-هوور (Nose-Hoover) که روش کاری آن بسیار دقیق و سریع بوده و در شبیه­‌سازی­‌هایی تحت فرآیند دما ثابت بسیار کارایی دارد. ترموستات دینامیک لانژوین (Langevin Dynamics) که از روال اتفاقی تبعیت کرده و برای کاربران پکیج نمد بسیار پرکاربرد است.

هم بستگی فشاری

باروستات الگوریتمی است که با اصلاح معادله­ ی حرکت نیوتن سعی در تولید کمیتی آماری مثل فشار دارد. هم­بسته ­سازی فشار نیز توسط باروستات­ های برندسون، دینامیک لانژوین (روش پیستون) و پارینلو-رحمان (Parrinello-Rahman) اعمال می‌­شود.

  برای توضیح بیشتر در باره این موضوع می‌­توان به منابع زیر مراجعه نمود.

[1]  Canonical sampling through velocity rescaling

[2] Thermostat algorithms for molecular dynamics simulations

[3] Molecular dynamics with coupling to an external bath     

[4] A molecular dynamics method for simulations in the canonical ensemble 

[5] Constant pressure molecular dynamics algorithms

[6] Canonical dynamics: Equilibrium phase-space distributions

[7] Strain fluctuations and elastic constants

هنگردهای ترمودینامیکی

با توجه به قانون گاز ایده ­آل، در صورت ثابت نگهداشتن حجم و دمای یک گاز، فشار آن نیز بدون تغییر باقی خواهد ماند. این موضوع به آن معنا است که متوسط حجم گاز در بازه­ ی زمانی طولانی از قانون گاز کامل تبعیت می­کند درحالی­که حجم لحظه­ ای چنین نیست. در سیستم­ های میکروسکوپی در شبیه­ سازی ­های مولکولی، حجم لحظه­ ای سیستم همواره نوساناتی با پریود زمانی کوچک دارد. بنابراین لازم است تا چندین تصویر میکروسکوپیک سیستم در بازه­ ی زمانی طولانی جمع ­آوری شده تا میانگین­ گیری معناداری از خواص ترمودینامیکی حاصل شود. در ترمودینامیک آماری، مجموع این تصاویر که متوسط­ گیری خواص را ممکن می ­سازد تحت عنوان هنگرد (Ensemble) آماری شناخته می ­شود. بنابراین بصورت خلاصه می­ توان گفت که مجموعه وضعیت­ های امکان ­پذیر برای سیستم که حالت میکروسکوپیک متفاوت اما حالت ماکروسکوپیک یا ترمودینامیکی یکسان دارند، هنگرد گفته می­ شود. همچنین بایستی اشاره کرد که احتمال وقوع هر حالت میکروسکوپی مقدار مشخصی است. شکل ‏زیر طرح شماتیک ویژگی­ های برخی از هنگردهای مرسوم در شبیه­ سازی دینامیک مولکولی را نمایش می­ دهد. در ادامه توضیح برخی از متداول ترین این هنگردها آورده شده ا

ست.

هنگرد میکرو کانونیکال

در هنگرد میکروکانونیکال (Microcanonical ensemble یا NVE)، سیستم از تغییرات مول N (تعداد ذرات)، حجم V و انرژی E جدا شده است. این هنگرد که پایه ­ای ترین هنگرد در دینامیک مولکولی است معادل یک فرآیند بدون تبادل حرارت است. به بیان دیگر در این هنگرد انرژی جنبشی و پتاسیل با هم مبادله می­ شوند ولی انرژی کل E ثابت باقی می­ ماند.

هنگرد کانونیکال

زمانی که خواص ترمودینامیکی یک سیستم در دما و حجم ثابت محاسبه شوند، هنگرد کانونیکال (NVT) استفاده شده است. تعداد ذرات N و حجم V برای کل سیستم ثابت در نظر گرفته می­ شود اما انرژی کل E متغیری در حال نوسان برای این هنگرد است. این هنگرد گاهی دینامیک مولکولی دما ثابت نامیده می شود (Constant Temprature Molecular Dynamics). در NVT انرژی گرماگیری و گرمادهی (Exothermic) فرآیند با یک ترموستات مبادله می­ گردد. چندین نمونه الگوریتم ترموستات که انرژی را از مرزهای دینامیک مولکولی به روش­ های واقعی و غیرواقعی اضافه یا کم می­ کند هنگرد کانونیکال را تقریب می­ زند. دمای ثابت سیستم به این معنا نیست که هر کدام از حالت­ های سیستم در دمای T است بلکه انرژی تمامی حالت­ های سیستم توزیع مشخصی دارد. 

هنگرد NPT

وقتی سیستمی در فشار و دمای مشخصی شبیه­ سازی می­شود، هنگرد فشار-دما ثابت، یا به اصطلاح NPT برآن اعمال شده است. در این هنگرد، مقدار انرژی و حجم نوسان دارند. ضمنا علاوه بر ترموستات، یک باروستات هم برای کنترل فشار مورد نیاز است. این نوع هنگرد بیشتر با شرایط آزمایشگاهی با یک مخزن باز به سمت دما و فشار محیط متناظر است. در شبیه­ سازی غشاهای بیولوژیکی، کنترل فشار ایزوتروپ مناسب نیست.

هنگرد NPnAT نوعی از هنگرد NPT بوده که در شبیه­ سازی سیستم­ های سطح مشترکی بسیار مفید است . شکل زیر طرحوره­ ای از سیستم سطح مشترکی با دو مایع امتزاج­ ناپذیر را نشان می­ دهد. تفاوت بین NPT و NPnAT در نحوه ­ی کنترل فشار است. در هنگرد NPT فشار بصورت ایزوتروپیک در تمامی جهات کنترل می­شود. به این معنا که تمامی ضلع­ های جعبه با حفظ نسبت یکسان طی شبیه ­سازی دینامیک مولکولی در گذر زمان بروزرسانی می­ شوند. این در حالی است که در هنگرد NPnAT فشار فقط در یک جهت کنترل شده و طول ضلع جعبه نیز تنها در این جهت دچار تغییر می ­شود. بایستی توجه داشت که سطح مماس با فشار کنترلی ثابت باقی می ­ماند. این ویژگی زمانی بسیار کارا است که سطح مشترک مورد بررسی، سطح ثابت در این هنگرد در نظر گرفته شود. در این حالت فشار توده­ ای در هر فاز (مایعات دو طرف سطح مشترک) با فشار نهایی مشخصی کنترل می­ ش