Enabling Ab-Initio Molecular Dynamics under Bias: The CP2K+SMEAGOL Interface for Integrating Density Functional Theory and Non-Equilibrium Green Functions

ملخص تقني: تمكين الديناميكا الجزيئية من المبادئ الأولية تحت تأثير الانحياز عبر واجهة CP2K+SMEAGOL

بيان المشكلة
بينما تُعد نظرية الدالة الوظيفية للكثافة (DFT) مدمجة مع دالات غرين غير المتوازنة (NEGF) نهجاً قياسياً لنمذجة النقل تحت جهود الانحياز الخارجية، فإن المنهجيات التقليدية تظل محدودة بالأنظمة الساكنة. فالنماذج الحالية تعمل عادةً ضمن إطار كانوني، حيث تحسب خصائص النقل عند مواضع ذرية ثابتة. وهذا يخلق فجوة كبيرة بين الفهم النظري والنتائج العيانية للتجارب، لا سيما في الأنظمة الكهروكيميائية (EC) حيث تكون التحولات بطبيعتها ذات طبيعة كانونية كبرى. وللوصف الدقيق للبيئة الطبيعية للتحولات الكهروكيميائية، تتطلب المحاكاة ديناميكا جزيئية (MD) تحت تأثير الانحياز ووصفاً صريحاً للحدود المفتوحة للإلكترونات التي يمكن أن تدخل وتخرج من الخلية الحسابية. علاوة على على ذلك، تفتقر تطبيقات DFT-NEGF الحالية غالباً إلى القدرة على حساب القوى المستحثة بالتيار اللازمة لتحسين الهندسة أو الديناميكا، أو تكون مقيدة بأكواد DFT محددة (مثل SIESTA) التي تفتقر إلى الميزات الموجودة في الحزم المتطورة والمُحسنة لمحاكاة الأطوار المكثفة.

المنهجية
يقدم المؤلفون واجهة جديدة بين حزمة البنية الإلكترونية الشهيرة التي تدمج بين الدوال الموجية الغاوسية وموجات المستوى CP2K وكود NEGF المسمى SMEAGOL. صُممت هذه الواجهة، التي تسمى CP2K+SMEAGL، لتكون مكتبة مستقلة يمكن ربطها مع أي برنامج DFT، مع إعادة استخدام كود SMEAGOL الأصلي قدر الإمكان.

تشمل التنفيذات التقنية الرئيسية ما يلي:

• القوى المستحثة بالتيار: تقوم الواجهة بحساب القوى غير المحافظة التي تؤثر على الأيونات تحت تأثير الانحياز. تُشتق هذه القوى من المشتق الزمني للزخم الذري، باستخدام تدرجات الهاملتوني الكامل بالنسبة للمواضع الذرية. ومن الأهمية بمكان أن التنفيذ يتضمن حد القوة الناجم عن مشتق مصفوفة الكثافة ($\Omega_{ij}$)، والذي يحسب مكون "رياح الإلكترونات" (انتقال الزخم من الإلكترونات إلى الأيونات). هذا الحد غالباً ما يتم إهماله في التنفيذات الأخرى، ولكن ثبت أنه ضروري للصحة النوعية.
• بناء مصفوفة الكثافة: لضمان إشغال الحالات المقيدة بشكل صحيح داخل نافذة الانحياز، تقوم الواجهة بتحديث حساب مصفوفة الكثافة. فبدلاً من استخدام تقييم مسار واحد أو متوسط بسيط للمساهمات من اليسار واليمين، تستخدم الواجهة مخطط المسار المزدوج الموزون (تبعاً لـ Brandbyge وآخرون) الذي يجمع التكاملات من كلا الطرفين بأوزان تحددها المساهمات غير المتوازنة.
• إعداد المحاكاة: تسمح الطريقة بإجراء ديناميكا جزيئية "بورن-أوبنهايمر" تحت تأثير الانحياز. يتم نمذجة النظام بـ "جزيء ممتد" مركزي متصل بأقطاب كهربائية شبه نهائية. يتم التعامل مع طاقات التلاشي (self-energies) للأقطاب كقيم ثابتة ومستقلة عن كثافة الشحنة في المنطقة المركزية، وهو افتراض سليم عندما تكون المنطقة المركزية كبيرة بما يكفي لحجب تغيرات الشحنة.

النتائج الرئيسية والتحقق من الصحة
قام المؤلفون بالتحقق من صحة واجهة CP2K+SMEAGOL مقابل الحسابات المرجعية التي أجريت باستخدام SIESTA+SMEAGOL عبر عدة أنظمة:

1. قوى جهد الصفر (سلك ذهب لانهائي): نجحت الواجهة في إعادة إنتاج القوى الذرية عند جهد الصفر، مطابقةً حسابات CP2K وحدها. كما أظهرت أن إهمال حد مصفوفة كثافة الطاقة ($\Omega$) يؤدي إلى قوى غير صحيحة نوعياً.
2. المكثف ذو الصفائح المتوازية: تعيد الواجهة إنتاج هبوط الجهد المتوقع عبر المكثف، حيث تظهر جهداً ثابتاً في الأقطاب وهبوطاً خطياً في الفراغ. وتؤكد أن مخطط المسار المزدوج الموزون ضروري لإنتاج كثافات شحنة متماثلة وجهود هارتري صحيحة عند الانحياز العالي، بينما تعطي تقييمات المسار الواحد نتائج غير متماثلة وخاطئة.
3. تحسين هندسة وصلة Au-H$_2$-Au: تحت تأثير الانحياز المطبق (حتى 1.5 فولت)، تتنبأ الواجهة بشكل صحيح باستطالة رابطة H-H والإزاحات الذرية، بما يتوافق مع نتائج SIESTA+SMEAGOL السابقة. وتُعزى هذه التغييرات إلى التيار الكهربائي وليس المجال الكهربائي.
4. الديناميكا الجزيئية لسلك ذهب مذاب: أجرى المؤلفون أول عمليات محاكاة للديناميكا الجزيئية بنظام DFT-NEGF لنظام طور مكثف واسع النطاق (سلك ذهب أحادي الذرة مذاب في 166 جزيء ماء).
   • الاستقرار: كانت الديناميكا مستقرة مع حد أدنى من الانجراف الطاقي طويل الأمد، رغم طبيعة النظام المفتوح.
   • هجرة الإلكترونات: في مسار زمني مدته 200 فمتوثانية عند انحياز 1 فولت، لاحظت المحاكاة تأثيرات هجرة الإلكترونات، حيث انتقل الحد الأدنى لكثافة الإلكترونات المحلية من ذرة السلك المركزية إلى الذرة اليمنى، وهو ما يتوافق مع حجب الشحنة الموجبة للقطب الأيمن.
   • النفاذية: أظهرت أطياف النفاذية للسلك المذاب تشابهاً نوعياً مع سلك مثالي لانهائي، مع وجود اختلافات تُعزى إلى قنوات نفق إضافية فتحتها جزيئات الماء.

الأداء
تتمثل العقبة الحسابية الرئيسية في تقييم الكثافة الإلكترونية كدالة تكامل لدالة غرين، مما يتطلب 96 نقطة طاقة (64 على المحور الحقيقي، و32 على المحور المركب) لكل خطوة SCF، بدلاً من عملية قطر واحدة للهاملتوني.

• التوازي: يستخدم الكود MPI للحسابات في الفضاء الحقيقي ونقاط k، وOpenMP لحلقات DO وإجراءات LAPACK. يعمل الأداء بشكل جيد حتى 32 خيط OpenMP.
• التكلفة: بالنسبة لنظام سلك الذهب المذاب، كانت عملية حساب CP2K+SMEAGOL أبطأ بمقدار 3.2 مرة لكل خطوة SCF مقارنة بحساب CP2K القياسي، رغم استخدام 8 أضعاف عدد النوى (2048 مقابل 256). ويشير المؤلفون إلى أن هذا يمثل زيادة في التكلفة بمقدار رتبة مقدارية، مما يجعل مسارات الديناميكا الجزيئية الطويلة أمراً صعباً دون تسريع إضافي (مثل ScaLAPACK أو حقول القوة القائمة على التعلم الآلي).

الأهمية والادعاءات
يزعم البحث أن CP2K+SMEAGOL هي الأكثر غنى بالميزات والاكتمال بين تطبيقات DFT-NEGF المتاحة لـ CP2K حتى الآن، متفوقة على تطبيق CP2K الأصلي المقتصر على نقل نقطة $\Gamma$ بدون قوى.

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