در این آموزش سعی می کنیم همگرایی یه سری خواص رو بر حسب اندازه سلول واحد بررسی کنیم. سیستمی که مورد بررسی قرار میدیم مولکول آب هستش.
در انتهای مطلب فایل های ورودی لازم برای انجام این محاسبه قرار داده شده که شامل دو فایل شبه پتانسیل مربوط به اتم های اکسیژن و هیدروژن و یک فایل ورودی با پسوند fdf است.
محتوای فایل ورودی به شکل زیره:
SystemName Water molecule SystemLabel h2o NumberOfAtoms 3 NumberOfSpecies 2 %block ChemicalSpeciesLabel 1 8 O # Species index, atomic number, species label 2 1 H %endblock ChemicalSpeciesLabel AtomicCoordinatesFormat Ang %block AtomicCoordinatesAndAtomicSpecies 0.000 0.000 0.000 1 0.757 0.586 0.000 2 -0.757 0.586 0.000 2 %endblock AtomicCoordinatesAndAtomicSpecies MeshCutoff 300.0 Ry DM.Require.Energy.Convergence .true. DM.Energy.Tolerance 1.e-5 eV
در خطوط 1 تا 4 تعداد گونه ها و اتم های حاضر در سلول واحد رو وارد کرده ایم.
در خطوط 6 تا 9 انواع گونه های اتمی را وارد کرده ایم که در سیستم آب دو گونه اکسیژن و هیدروژن وجود داره.
در خطوط 11 تا 16 هم مختصات اتم ها رو وارد کرده ایم.
در واقع همانطور که در مثال قبلی هم گفته شد به این قبیل محاسبات که در اون ها از هیچ داده تجربی استفاده نمیشه محاسبات ابتدا به ساکن گفته میشه.
یه سری کمیات هم هستن که اگه تعریف نشن سیستا مقداری به صورت پیش فرض برای اون ها در نظر میگیره. مثال هایی از این کمیات و مقدار پیش فرضشون رو در زیر میبینید:
PAO.BasisSize (Basis set quality) DZP MeshCutoff (Fineness of real space integrations) 100 Ry XC.Functional (Exchange and correlation functional) LDA XC.Authors (Flavour of the exchange and correlation) CA SpinPolarized (Are we performing an spin polarized calc.) .false.
و خیلی کمایت دیگه که برای مشاهده لیست کاملشون می تونید بعد از اجرای برنامه فایل fdf.log رو ببینید.
الان زیاد نگران این ها نباشید در مثال های آینده کم کم این موارد توضیح داده میشه.
مولکول آب با شرایط مرزی دوره ای (PBC)
با اینکه سیستم ما (مولکول آب) یک سیستم غیر دوره ایه ولی سیستا از شرایط مرزی دوره ای استفاده می کنه. خوب یعنی در این شرایط محاسبات اشتباست؟
جواب این سوال «نه» هستش ما می تونیم با استفاده از یه سری تکنیک بر این مشکل غلبه کنیم. برای این کار از روش به اسم «ابر سلول» استفاده می کنیم.
استراتژی: روش ابر سلول
در این روش یک ناحیه خلأ اطراف سیستم تعریف می کنیم به شکلی که سلول های واحد عملا با هم برهمکنشی نداشته باشن. (شکل کاملا نشون میده که میخوایم چکار کنیم)
مرحله بعدی اینه که مطمئن بشیم خواص فیزیکی و شیمیایی مورد نظرمون نسبت به حجم ابرسلول همگرا شده باشن.
خوب سوالی که پیش میاد اینه که اگه ما سلول واحدی تعریف نکنیم خود سیستا چکار میکنه؟ اگه ما اینکارو نکنیم خود سیستا میاد یه سلول واحد میسازه که بردار های شبکه اون قطری هستن و اندازه اون ها به شکلیه که اولا شامل کل سیستم باشه بدون اینکه سلول های نزدیک با هم همپوشانی داشته باشن به علاوه یک لایه بافر به اندازه 10 درصد.
اولین اجرای سیستا بر روی این سیستم
دقت داشته باشید که نسبت به مثال قبلی که دو اتمی بود ما اینجا یه قدم جلوتر رفتیم و سیستمون سه اتمی شده (همانطور که میدونید توسعه دهندگان سیستا ادعا دارن که سیستا برای سیستم هایی در حد 1000 اتم به خوبی جواب میده)
خوب اول چک کنید که همه فایل های ورودی رو دارید بعد در مسیری که فایل رو گذاشتید یک ترمینال باز کنید و دستور زیر رو وارد کنید:
siesta < h2o.fdf > h2o.default.out
اینجا واژه default رو به این دلیل در نام فایل خروجی استفاده کردیم که یادمون باشه تو فایل ورودی اولیه اندازه سلول واحد داده نشده و سیستا به صورت پیش فرض اونو میسازه.
بعد از اجرای دستور فوق یک فایل خروجی خواهیم داشت.
خوب حالا اجازه بدید یه نگاهی به فایل های خروجی بیندازیم.
برای مثال فایل h2o.default.out رو نگاه کنید و به سوالات زیر جواب بدید.
- چند حلقه SCF نیاز بوده تا شرایط همگرایی حاصل بشه؟
- بعد از رسیدن به همگرایی انرژی کل سیستم جقدر بوده؟
- سلول واحدی که توسط سیستا تولید شده دارای چه ابعادیه؟
- دوقطبی الکتریکی برای این مولکول بر حسب الکترون بوهر چقدر بوده؟
حالا فایل ورودی رو تغییر بدید و یک ابرسلول به شکل زیر تعریف کنید:
حالا مقدار 8 آنگستروم رو یکی یکی تا 1 آنگستروم تغییر بدید و واسه هر کدوم یه خروجی بگیرید. می تونید این کارو بعد از تغییر ثابت شبکه و اجرای دستوری به شکل زیر انجام بدید:
siesta < h2o.fdf > h2o.your_lattice_constant.out
که به جای your_lattice_constant باید ثابت شبکه ای رو که در فایل ورودی در هر مرحله وارد می کنید قرار بدید.
خوب الان خروجی ها آماده است و تنها چیزی که نیاز دارید اینه که یه نمودار انرژی بر حسب سلول واحد رسم کنید و نتیجه رو بررسی کنید.
برای جمع آوری انرژی کل از فایل خروجی می تونید از دستور grep به شکل زیر استفاده کنید:
grep "Total =" h2o.*.out > h2o.latcon.dat
این دستور همانطور که در مثال قبلی هم گفته شد میره تمام فایل های که اول اسم اونها h2o. و اتنهای اسمشون .out هست میگرده عبارت Total = رو پیدا میکنه و خط شامل این عبارت رو در فایل h2o.latcon.dat ذخیره می کنه. بعد از این کار فایل h2o.latcon.dat رو تغییر بدید و فقط دو ستون به شکل زیر باقی بزارید:
دقت کنید که با توجه به نسخه سیستا و کامپایلری که استفاده شده ممکن این اعداد برای شما اندکی تغییر کنه.
خوب حالا داده ها رو هم دارید و می تونیم نمودار بکشیم. برای اینکار اول در ترمینال دستور gnuplot رو تایپ می کنیم و اینتر می زنیم بعد دستور زیر رو وارد می کنیم:
plot "h2o.latcon.dat" using 1:2 with lines
بعد از اجرای این دستور نموداری به شکل زیر نشون داده میشه:
بررسی نتایج
در حالت ایده آل برای یک مولکول هر ویژگی می بایست مستقل از اندازه جعبه شبیه سازی باشه. اما در مورد مولکول آب (حداقل برای انرژی) این مسئله برقرار نیست. چرا؟
مولکول آب دارای یک دوقطبیه که اندازه اون در این محاسبه به شکل زیر ماحسبه شده:
وقتی که انرژی یک سیستم غیر پریودیک با استفاده از شرایط مرزی پریودیک محاسه می کنیم، چیزی که بهش علاقه مندیم محاسبه E0 در حد L به سمت بی نهایته، که L بعد ابر سلول ما هستش. انرژی که برای یک ابر سلول محدود محاسبه میشه یعنی (E(L با E0 فرق می کنه چون چگالی بار غیر پریودیک با تصویر خودش در سلول های همسایه برهمکنش می کنه.
برای تخمینE0 از مقدار محاسبه شده (E(L نیاز داریم که رفتار مجانبی انرژی بر حسب L رو بدونیم.
از اونجایی که مولکول آب دارای دو قطبیه برهمکنش الکترواستاتیکی بین دو قطبی به صورت 1 تقسیم بر توان سوم L کاهش پیدا میکنه یعنی داریم:
که متغیر های به کار رفته در اون رو می تونید از شکل زیر متوجه بشید:
ولی این بستگی کافی نیست. چون این برهمکنش ها تغییراتی در چگالی بار غیر دوره ای ایجاد می کنه که خودش به L بستگی داره. به عبارت دیگه دو قطبی خودش هم به L وابسته است. اگه مولکول یک دو قطبی دائمی داشته باشه اونوقت دوقطبی القایی خودش تابعی از عکس توان سوم L هستش زیرا میدان تولید شده توسط یک دوقطبی به شکل عکس توان سوم L کاهش پیدا می کنه.
پس ما باید همگرایی یه کمیت دیگه رو هم بررسی کنیم و اون چیزی نیست جز دو قطبی الکتریکی!!!
برای این کار اول باید داده های مربوط به دو قطبی الکتریکی رو هم استخراج کنیم برای این کار میتونید از دستور زیر استفاده کنید:
grep "Electric dipole (a.u.)" h2o.*.out > h2o.dipole.dat
بعد از اجرای این دستور فایل h2o.dipole.dat رو ویرایش کنید و فقط ثابت شبکه مولفه های ممان دوقطبی رو نگه دارید (با توجه به جهت گیری که مولکول در این مثال داره فقط مولفه y یر صفر خواهد بود).
بعد از این کار ابتدا در ترمینال دستور gnuplot رو وارد کنید بعد از اون اینتر بزنید و دستور زیر رو وارد کنید:
plot "h2o.dipole.dat" using 1:3 with lines
بعد از اجرای دستور فوق نموداری به شکل زیر به شما نشون داده میشه:
رابطه ممان دوقطبی با اندازه L به شکل زیر هستش:
که p0 و c پارامترهایی هستن که بعد از فیت کردن می تونید بدست بیارید.
حالا اگه تو فرمول W به جای p عبارت (p(L رو قرار بدیم اونوقت سهمی که منجر به انرژی الکترواستاتیکی میشه منعکس کننده برهمکنش دوقطبی (دوقطبی القایی) هستش که از مرتبه عکس توان ششم L خواهد بود:
E0 ، a و b پارامتر هایی هستن که از فیت کردن بدست میان. چیزی که ما بهش علاقه مند هستیم مقدار E0 هستش.
حالا شما خودتون این محاسبه رو انجام بدید و یه مقداری برای E0 تخمین بزنید.
تا آموزش بعدی 🙂
[otw_shortcode_button href=”http://www.comphys.ir/wp-content/uploads/siesta/Files-h2o.tar” size=”medium” icon_position=”left” shape=”square” target=”_blank”]دانلود فایل های مورد نیاز[/otw_shortcode_button]
یک دیدگاه
Or use one of these social networks