Ab initio розрахунки динаміки гратки сегнетоелектриків Sn2P2S6

Міністерство освіти і науки, МОЛОДІ ТА СПОРТУ україни ДЕРЖАВНИЙ ВИЩИЙ НАВЧАЛЬНИЙ ЗАКЛАД

" Ужгородський національний університет"

фізичний факультет кафедрА тВЕРДОТІЛЬНОЇ ЕЛЕКТРОНІКИ

Дипломна робота з теми:

Ab initio розрахунки динаміки гратки сегнетоелектриків Sn2P2S6

УЖГОРОД — 2011

Реферат

В даному роботі було проведено ab initio розрахунки в процесі оптимізації геометрії, значення постійних гратки та координат атомів в парата сегнетофазі кристалу методом функціонала густини (DFT). Аналіз отриманих результатів свідчить про те, що в умовах оптимізації геометрії розглядувані структури зберігають топологію кристала Sn2P2S6, а розраховані значення геометричних параметрів відтворюють відповідні експериментальні результати. Обчислювали в центрі зони Брюлльена фононних мод та їх симетрії

Вступ

Кристали гексахалькогіподи-фосфати Sn (Pb)2P2S (Se)6 поєднують у собі напівпровідникові та сегнетоелектричні властивості, мають високі піроелектричні, п'єзоелектричні та електрооптичні характеристики, що робить їх перспективними для практичного використання. З метою отримання нових перспективних напівпровідникових матеріалів в даний час інтенсивно проводяться дослідження фізико-хімічної взаємодії та вирощування кристалів у псевдобінарних системах Sn2P2S6-Ме2P2S6(Me: Zn, Cd, Fe, Mn, Ni).

Задача була оптимізувати кристалу Sn2P2S6 і глибоко розуміти мікроскопічні властивості, також дослідити фононних спектрів. Велика допомога для обчислення даних задач була, ab initio методи розрахунки фізичних властивостей.

Методи квантової хімії дозволяють обчислювати рівноважні конфігурації молекул і енергії різноманітних конформацій, а також отримати строгі енергетичні характеристики та електронно-спінові стани.

У програмному коді Abinit реалізований вибір базису у вигляді приєднаних плоских хвиль. Кожен атом при цьому оточується деякої сферою (атомною сферою), всередині якої хвильова функція шукається у вигляді розкладу по базису плоских хвиль так, щоб поза сферою можливо було «зшити» цю функцію з функцією вільного електрона. У програмі Abinit, крім того, використовується метод псевдопотенціалу, що дозволяє замінити частину базисних функцій на деякий ефективний потенціал. Така заміна проводиться автоматично при розрахунку. Від користувача потрібно лише вибрати псевдопотенціал необхідного типу.

Отримавшиоптимізовануструктурукристалів Sn2P2S6, наступним кроком було обчислення фононних спектрів досліджуваних кристалів, використовуючи при цьому оптимізовані координати елементарної комірки кристалів Sn2P2S6.

Для аналізу симетрії фононних мод була використана утиліта ANADDB (Analysisof Derivative Data Base), яка дозволяє також обчислювати безліч різних фізичних властивостей, що містять похідні повної енергії (похідна Data Base-DDB).

Для прикладу наведено порівняння розрахованих нами фононних мод (сегнетофаза — по симетрії Ag, парафаза — по симетрії A") в центрі зони Бріллюена з існуючими теоретичними та експериментальними результатами.

РОЗДІЛ 1. СТРУКТУРА І ФІЗИЧНІ ВЛАСТИВОСТІ КРИСТАЛІВ Sn2P2S6

1.1 Кристалічна структура

Кристали тіогіподифосфату олова (Sn2P2S6) належать до сімейства халькогенідних напівпровідникових сегнетоелектриків (ХНС). Вони є іоно — ковалентними сполуками, які складаються із аніонів [P2S6]4-, що з'єднані в тривимірну гратку за допомогою катіонів Sn2+. При кімнатній температурі Sn2P2S6 володіє нецентросиметричною структурою моноклінної сингонії (Pc). Кристалографічна елементарна комірка містить дві формульні одиниці (Z=2). Аніони [P2S6]4-, що описуються деформованою точковою групою симетрії, складаються з двох деформованих тригональних пірамід PS3. Всі шість P-S зв’язків мають різну довжину в межах від 2,015 до 2,035. Піраміди з'єднані зв’язком P-P довжиною 2,202. Атоми олова розміщені між аніонами всередині неправильних поліедрів, які утворені 7 або 8 атомами сірки. Іонні зв’язки Sn-S (2,776 — 3,451) значно довші, ніж P-S і P-P зв’язки.

Проекція структури Sn2P2S6 на площину (010), що є площиною симетрії показана на рис 1.1. Параметри ґратки у сегнетоелектричній фазі: a = 9,378; b = 7,488; c = 6,513; в = 91,15?. При нагріванні кристали переходять у параелектричну фазу (Т0 = 337 К), і при T > T0 структура Sn2P2S6 описується просторовою групою P21/c. Параметри ґратки при температурі 358 К: a = 9,318; b = 7,436; c = 6,518; в = 91,15?. Із порівняння із ацентричною структурою сегнетофази слідує, що в результаті фазового переходу об'єм комірки практично не змінюється, тоді як її форма дещо деформується за рахунок видовження вздовж Z на 0,23 і вздовж Y на 0,02(осі декартової системи координат X і Y суміщені з напрямками і відповідно, а вісь Z близька до напрямку).

Рис. 1 — Елементарна комірка кристалу Sn2Р2S6 (а) та її проекція на площину (001) (б) [33]

Форми координаційних багатогранників атомів халькогену помітно не відрізняються. При T > T0 тільки атоми Sn значно зміщені в напрямку порівняно з положенням металу в центросиметричному стані. Це дозволяє вважати катіонну підгратку Sn2+ сегнетоактивною.

Механізм міграції атомів олова при переході із сегнетоелектричної в парафазу можна представити наступним чином. Всі чотири атоми Sn зміщуються приблизно на 0,26 в напрямку [100], досягаючи при цьому центросиметричних позицій. В той же час два з них (не зв’язані площиною ковзання) також зміщуються на 0,04 по [010], два інших зміщуються на таку ж відстань по. Зміна позицій атомів сірки і фосфору з наближенням до T0 є дуже незначною. Відстань від атомів олова до найближчих атомів сірки в сегнетофазі рівна 2,76, тоді як в парафазі вона складає 2,90.

Таким чином, зміщення атомів олова проходять головним чином вздовж напрямку. Ці зміщення є фізичною реалізацією параметра порядка сегнетоелектричного фазового переходу другого роду при T0 = 337 К в Sn2P2S6.

1.2 Симетрійний аналіз

У параелектричної фазі кристали системи Sn2P2S6 описуються просторовою групою P21/c[2]. Зона Бріллюєна, яка відповідає примітивній моноклінній елементарній комірці показано на рис. 1.2. Вектори оберненої гратки і лежать в площині симетрії, а перпендикулярний даній площині, тобто співпадає з віссю другого порядку. Буквами позначені високосиметричні точки, лінії і площини, значення хвильових векторів, які наведені в таблиці (1.1). В таблиці також показані елементи симетрії групи P21/c, незмінні вектори. Вони складають групу хвильового вектора.

Рис. 1.2 — Зона Брілюена, для кристалічної гратки просторової групи P21/c [5]

Приведені незвідне представлення, відповідно кожному хвильовому вектору згідно даної роботи. Таким чином, було виконане теоретично-груповий аналіз нормальних коливань кристалів для точки Г, можна найти їх розподіл по типам симетрії для любого точка зона Бріллюена.

Табл. 1.1 — Високосиметричні точки, лінії і площини в зоні Бріллюєна просторової групи P21/c

— 0.5 < м1, м2, м3? 0.5

Розглянемо результати симетрійного аналізу фононних спектрів для типу кристалів Sn2P2S6. В центросиметричній структурі аніони (P2S6)-4 знаходяться в позиціях, а катіони Sn2+ - позиціях 1[1]. Точкова група симетрії аніона (P2S6)-4 у вільному стані -. Коливальне представлення для внутрішніх мод ізольованого аніона, отримане БВ методом (рис) має вигляд:

(1.1)

Кореляційна діаграма, яка дає відношення між незвідне представленнями точкових груп, , 1 і фактор-групи парафази кристалу 2/m, статистичне розчеплення два вироджених коливань ізольованого іона (P2S6)-4 із-за коливань кристального поля. Наявність двох іонів в елементарній комірці наслідок резонансної взаємодії приводять до розчеплення кожного коливального іона на компоненти або. В кореляційну діаграму включені також зовнішні трансляційні і лібраційні коливання ізольованого іона (P2S6)-4. В кристалі вони перетворюються в типу і 3Ag+3Bg відповідно. Три степені вільності іона Sn2+ дають вклад в трансляційні коливання гратки в парафазі Sn2P2S6, розподіляючись між представленнями фактор-групу 2/m наступним чином: 3Ag+3Bg+3Au+3Bu.

Загальний вигляд незвідного представлення для оптичних нормальних коливань парафази Sn2P2S6 [8]:

(1.2)

Для центросиметричної структури парафази Sn2P2S6 виконується правило альтернативної заборони, згідно якій парні відносно центра інверсії коливання (g) активні тільки в КР, а непарні коливання (u) — тільки в ІЧ спектрах. Для акустичних коливань парафази:

(1.3)

Розрахунок БМ методом приводить до таких же результатів (табл. 1.2). Тут збережені позначення, прийняті в [212]: n1, T, T', R, n1', — відносно числа всіх нормальних коливань, трансляційно-акустичних коливань, зовнішніх коливань трансляційного типу, зовнішній коливань обертового типу, внутрішніх коливань і повне число оптичних коливань.

При переході в сегнетофазу із зміною просторової групи P21/c коливання розділяються по симетрії лише відносно площини ковзання — симетричні і асиметричні. Кореляція між типами симетрії точкових груп 2/m і m показано в табл. 1.2. Тоді незвідне представлення для оптичних коливань сегнетофази Sn2P2S6 має вид [9]:

(1.4)

Для акустичних коливань:

(1.5)

Табл. 1.2 — Аналіз нормальних коливань свобідного іона (P2S6)-4 БВ методом

Результати теоретико групового аналізу коливальних спектрів парафази Sn2P2S6 застосовані для сполук Sn2P2Sе6, Pb2P2S6, Pb2P2Se6, які володіють при кімнатній температурі центросиметричною структурою (P2/c).

1.3 Густина фононних станів і термодинамічні функції

Інформація про густину фононних станів дозволяє розраховувати температурні зміни термодинамічних функцій і відіграє важливу роль в ангармонічних процесах для конкретних систем.

Густина фононних станів Sn2P2S6 і Sn2P2Se6 визначалася із спектрів непружного розсіювання повільних нейтронів[3]. У зв’язку з тим, що елементарні комірки досліджуваних кристалів складаються із 20 атомів, то безпосередньо з експерименту визначалась функція Gs (е). Ця функція представляє собою добуток функції густини фононних станів g (е) і функції Ks (g, е), яка залежить від амплітуди когерентного розсіювання, мас і векторів поляризації окремих коливальних атомів, які входять в елементарну комірку:

(1.6)

Де rчисло атомів в елементарному комірку; aj — фактор Дебая Валлера j — атома; оjs — вектор полярізації s — моди j — атома; p — імпульс, переданий нейтроном при розсіювання. В експериментах[34] функція Gs (е) розраховується із усереднених перерізи розсіювання в області кутів 70…1200 по методу, детально описаному в роботі.

Результати дослідження динаміки гратки кристалів типу Sn2P2S6 методом непружного розсіювання повільних нейтронів представлений на рис. 1.3.

Рис. 1.3 — Густина фононних станів кристала Sn2P2S6 (на вставці залежність)

На залежності густини фононних станів від енергії, звертає на себе увагу той факт, що в інтервалі низьких частот (е < 5MeB) функція Gs (е) для Sn2P2S6 має квадратичну залежність від енергії. В області 6…7 МеВ спостерігається злам в Gs (е), що відповідає як слідує із КГ спектра, вкладу в густини фононних станів, трансляція катіонів Sn2+. В області до 20 МеВ, де моди визначаються вібраціями аніонів (P2S6)-4, виділені максимуми відсутні. Це можна пояснити значно меншим ефективним перерізом розсіювання сірки (1,2 барн.) порівнянням з оловом (4,9 барн.). При е > 20 МеВ, де появляються внутрішні деформації S — P — S коливання, а також «гнучкі «моди (накладені PS3- груп) в аніонах (P2S6)-4, спостерігається наявність піків в Gs (е) Sn2P2S6. Значення частот, при яких спостерігаються піки непогано сходяться з результатами по комбінаційному розсіюванні. Внутрішні валентні P — S моди, які знаходяться вище 70 МеВ, в даному експерименті не спостерігались враховуючи слабе заселення цих мод при кімнатній температурі[34].

На рис. 1.4 показано частотний спектр Gs (е) сегнетоелектриків Sn2P2S6. Функції Gs (е) Sn2P2S6 і Sn2P2Sе6 суттєво відрізняється одна від одної. Так в Sn2P2Sе6 відсутній злам при 6 мЕв і спостерігаються інтенсивні піки в інтервалі енергій 7…20 мЕв. При енергіях нижче 4 мЕв спектри практично ідентичні (квадратична залежність Gs (е) = bе2). У високо частотній області спостерігається перерозподіл між різними піками, що узгоджується з результатами по комбінаційному розсіювання світла[34].

Рис. 1.4 — Густина фононних станів кристала Sn2P2S6 (на вставці залежність)

Наявність спектрів густини фононних станів дозволяє розрахувати температурну поведінку термодинамічних функцій. Макроскопічні теплові властивості і елементарні коливальні стани кристала зв’язані законами статистичної фізики. В гармонічному наближенні термодинамічні функції кристалу можна отримати із статистичної суми

(1.7)

де H — Гамільтоніан в гармонічному наближенні; k — стала Больцмана; T — абсолютна температура.

Вільну енергію Гельмгольца можна представити у вигляді [13]:

(1.8)

Із (1.8) слідує, що в гармонічному наближенні термодинамічна функції (Ev=F — T (dF/dT))v, Cv=(dF/dT)v, S=-(dF/dT)v) являються адитивними функціями власних частот щj (q). Таким чином, їх можна виразити як середню від спектрального густини g (щ). Наприклад теплоємність задається виразом [14]:

(1.9)

1.4 Дисперсія фононних спектрів

Опис коливального спектра кристала в гармонічному наближенні характеризується сукупністю законів дисперсії нормальних коливань, де щ — частота; - хвильовий вектор; j — номер дисперсного вітка. За допомогою методу теоретичного-групового аналізу можна встановити загальну характеристику дисперсних віток, яка випливає із симетрії елементарної комірки — розподілення по незвідним представленням, наявність вироджень в симетричних точках зони Бріллюена. Розрахунки фононних дисперсій дають можливість отримати форми коливань гратки при різних. Вони представляють особливий інтерес для вивчення динаміки сегнетоелектричних ФП, в яких основну роль відіграє м’яка оптична мода. Вид дисперсії м’якої вітки, а також других коливань, дозволяє дослідити характер взаємодій між ними. 34]

Задача по розрахунку частот і форм коливання кристалічної гратки в гармонічному наближенні може бути сформовано наступним чином [16, 17]. Розглянемо кристал із n атомів в елементарній комірці. Нехай вектор описує рівноважне положення k атома в l комірка. Відхилення частинок від положення рівноваги позначаємо uб (lk), де б = X, Y, Z, k=1,2…, n. Величину uб (lk), можна розглядати як узагальнені координати в так званому механічному просторі коливань. Розкладемо потенціальну енергію коливань по степені зміщення, вважаючи їх малими. Перший член розкладання описується потенціальною енергією в положенні рівноваги (всі uб (lk)=0) і може бути вибраний рівним нулю. Другий член містить першу похідну, також рівний нулю в положенні рівноваги. Для потенціальної енергії U отримає:

(1.10)

де. (1.11)

В такому наближенні в рівнянні руху фігурують тільки лінійні комбінації зміщення. Вони задаються виразами:

(1.12)

де mk — маса k-го атома. Коефіціент визначає силу, діючу в напрямку б на атом k в комірку l, коли атом в комірку зміщуються на одиницю довжини в напрямку в, і не залежить від вибору систем відліку.

Порівнюємо результати розрахунків коливального спектра кристала Sn2P2S6 з експериментальними даними. Згідно теоретично-групового аналізу в сегнетофазі для зовнішніх повносиметричних коливань A' очікуються, крім 2 акустичних, 7 трансляційних і 3 лібраційних оптичних мод, а для асиметричних A'' - 1 акустична, 8 трансляційних і три лібраційних оптичних мод. Нижче T0 аніони (P2S6)-4 деформуються, що приводить до зміщення трансляційних і лібраційних мод.

Відмитимо, що мінімальне значення хвильового вектора, для яких виконувалися розрахунки, складали qmin=0,001qmax. Це обмеження викликано необхідністю обліку ефектів запізнення для хвиль більших довжини. В подальшому, говорячи про розрахунку в центрі зони Бріллюена, будемо розуміти qmin. Аналіз в центрі зони Бріллюена власних векторів підтверджує їх розділення по типам симетрії (рис. 1.5). Однак більшість коливань, представляє собою змішаний характер руху. Здебільшого трансляційними являються моди при 67, 70 і 138 cm-1, які містять 92, 98і і 96% трансляцій відповідно. В моди при 42, 52 і 105 cm-1 вклад лібрації складає 75, 81 і 54%. При віддаленні від центру зони Бріллюена характер власних векторів змінюється. При цьому розділення по типам симетрії можливий тільки для? Y.

З рис 1.5 видно, що вітки симетрії A'' однакової частоти при наближенні до центру зони Бріллюена по напрямках qx і qz, де вони поперечні. В той же час частоти даних мод при qy >0 розміщенні вище, так як коливання A'', що поширюються перпендикулярно в площині симетрії, являються повздовжніми.

Верхня границя зовнішніх коливань кристалу Sn2P2S6, визначена методом співставлення коливальних спектрів його аналогів і вільного іона (P2S6)-4 [3], а також по зміненню нахилу частотної залежності параметрів мод Грюнайзена [12], які розміщені в околі 150cm-1. Отримані в розрахунках частоти верхньої оптичної моди при q0 рівна 137 cm-1.

Нахил акустичних віток повинен співпадати з швидкостями пружних хвиль. Наприклад, для повздовжніх акустичних хвиль по осі Y виміряна швидкість при 295 K складає 2,9*105 cm/c [15], відповідає пунктирна лінія на рис. 1.5. Це значення задовільно узгоджується з розрахунком 2,5*105 cm/c. Зауважимо, що для всіх напрямків характер акустичних віток задовільняє симетрійного вимогам. Наприклад, верхні акустичні вітки завжди являється повздовжніми або квазіповздовжніми.

Вони добре узгоджуються з експериментальними даними для Sn2P2S6. На рис. 1.5 точками показано значення частот, отримані із спектрів КГ при 295 К.

Динамічні розрахунки виконані з використанням структурних даних, результаті з спектроскопічних ультразвукових досліджень параелектричної фази кристалів Sn2P2S6 і Sn2P2Se6.

Кінцевий набір параметрів вибрано по найкращому узгодженню результатів з експериментальними значеннями частот коливання гратки в центрі зони Бріллюєна і швидкості звука. Використані значення Z i R рівні 1,2; 0.6; -0.6 i -0.6 для зарядів і 2,45, 1,08, 1,40 і 1,55 для ефективних радіусів Sn, P, S i Se відповідно. Отримані при таким наборів дисперсійні криві представлені на рис. 1.6. Відмітимо, що значення параметрів для параелектричної фази Sn2P2S6 близькі до розрахунків для динаміки сегнетофазі цього кристала.

Рис. 1.5 — Дисперсійні криві сегнетофази кристала Sn2P2S6; здовж qx i qz вітки A' зображений суцільними лініями, А" - пунктирними лініями; точками показані частоти, спостерігаючи в спектрах КГ; світлі - А', темні А"; хвильові вектори приведені в частках векторів оберненої гратки Рис. 1.6 — Дисперсія фононів параелектричної фази кристалу Sn2P2S6 в моделі жорстких іонів; дисперсійні віток симетричних мод (А для qy і А' для qz) зображені суцільними лініями не симетричних мод (В для qy і А" для qz) — штриг-пунктирними; нахил штрихових ліній відповідає швиткостям ультразвукових хвиль визначених експериментально Згідно результатів досліджень комбінаційних спектрів розсіювання кристалів Sn2P2S6(Se)6 і твердих розчинів на їх основі, частоти внутрішніх коливань знаходяться в інтервалі 30…140 см-1 і 20…120 см-1 для Sn2P2S6(Se)6. Отримані розрахунки значення частот при q? 0 (рис. 1.6) добре співпадають з експериментальними даними. Розподіл зовнішніх оптичних коливань по типам симетрії точкових груп параелектричної фази наступний: 4Bu+6Ag+6Bg+5Au. Аналіз в центрі зони Бріллюена власних векторів розрахованих коливань підтверджує такі їх симетрійні розподілення. При віддаленні від центру зони Бріллюена в площині симетрії XZ коливання розділяються на 10A'+11A". В напрямку q¦Y дисперсійні вітки мають симетрію 11A+10B.

Спостерігається не співпадіння частот при наближенні до центру ЗБ по різним напрямках. Воно обумовлене горизонтально-поперечним (LO — TO) розщепленням коливань. Із рис. 1.6 видно, що неполярні вітки Ag i Bg такого розчеплення не знайдені. Моди Bu i Au поляризовані в напрямках X, Z i Y. аналіз частот і поляризація цих коливань по різним направленням поблизу центра ЗБ засвідчує про виконані умови щLO? щTO. Вироджені вітки на границі ЗБ в напрямках qz i qy узгоджуються з результатами симетрійного аналізу.

Результати вимірювань швидкостей горизонтальних і поперечних пружних хвиль в напрямку X, Y i Z для кристалів Sn2P2S6 і Sn2P2Se6 задовільно узгоджується з розрахунками (пунктирні лінії на рис. 1.6).

Розділ 2. Квантово-хімічні методи

2.1 Основні квантово-хімічні методи

Квантово-хімічні методи базуються на розв’язанні рівняння Шредінгера (2.1) різноманітними методами.

(2.1)

Першим швидким числовим методом був метод Хартрі-Фока, де хвильова функція системи електронів шукалася як лінійна комбінація атомних хвильових функцій. Таке наближення носить назву Методу Молекулярних Орбіт як Лінійних Комбінацій Атомних Орбіт (ММО-ЛКАО). Особливість методу Хартрі-Фока полягає в тому, що міжелектронні взаємодії замінюються взаємодією кожного електрону з усередненим полем, яке утворюється всіма електронами[38].

Із назви методу видно, що молекулярна орбіта апроксимується лінійною комбінацією атомних орбіт. Базис атомних орбіт задається користувачем і називається базисом розрахунку. Як правило, для апроксимації атомних орбіт вибирається в якості базисної функції лінійна комбінація гауссіан. Точність розрахунку в значній мірі визначається величиною базису.

Обмежений метод Хартрі-Фока придатний для систем із замкнутими електронними оболонками, де на кожну орбіту знаходяться по два спарених електрона. Для розімкнутих систем застосовується необмежений метод Хартрі-Фока, який використовує вдвічі менше процесорного часу, ніж обмежений.

Метод Хартрі-Фока має два суттєві недоліки. По-перше, відсутнє врахування енергії кореляції або енергії Ван-дер-Вальсовських взаємодій. За рахунок усереднення поля всіх електронів, кореляція руху електронів, яка може існувати в системі усереднюється і не враховується. По-друге, в такому розрахунку не враховуються взаємодії збуджених електронів з незбудженими, обчислення енергії збуджених станів не вірні. Із-за цього метод Хартрі-Фока неможна використовувати для розрахунку спектрів.

Покращити розрахунок в межах методу Хартрі-Фока можна, скориставшись теорією збурень Мюлера-Плессета, яка дозволяє враховувати додаткові міжелектронні взаємодії, які не входять в рівняння ХФ як збурення.

Існують ще два потужні методи розрахунку: розрахунок за допомогою Теорії Функціоналу Густини (DFT) і розрахунок методом Багато-Конфігураційної Взаємодії (MCSCF).

Теорія функціоналу густини дозволяє обчислити енергію системи для заданого емпіричного функціоналу. Метод дуже зручний у використанні, тобто не потрібно задавати ніяких спеціальних параметрів розрахунку крім власне функціоналу густини. Розрахунок за допомогою DFT дозволяє коректно врахувати енергію кореляції.

MCSCF, напроти, розрахунок дуже складний як для коректного налаштування, так і для процесора: MCSCF — найдовший вид розрахунку. Його доцільно використовувати якщо потрібна дуже висока точність обчислення і якщо вам необхідно розраховувати спектр молекули. Метод базується на представленні скінченої хвильової функції лінійною комбінацію основного і збудженого стану. Можлива комбінація MCSCF і методу теорії збурень.

Іноді, навпаки потрібне спрощення розрахунку. Існують ряд напівемпіричнх методів для розрахунку типу Хартрі-Фока: PM3, AM1. Такий розрахунок займає дуже мало процесорного часу і добре підходить для оцінки розрахунків.

Існують додаткові поправки, які покращують розрахунок. До таких поправок відноситься врахування оточення (розчинника). Він може проводитись в рамках Моделі Поляризованого Континууму (РСМ) або в моделі СамоУзгодженого Реакційного Поля Кірквуда-Озагера (SCRF). Обидва метода конструюють порожнину, куди поміщається система і обчислюються наведені на неперервному діелектрику досліджуваної системи. SCRF будує сферичну область, а РСМ будує область із складною поверхнею, яка утворюється при суміщенні сфер із центром на кожному атомі системи. РСМ на сьогодні являється найкращим методом врахування оточення, однак, якщо ваша система має велику кількість атомів, то необхідно самостійно задавати порожнину, що може виявитись незручним. Крім цього в програмному пакеті GAMESS РСМ не підтримується деякими видами розрахунків, наприклад UHF, ROHF.

Існує численна кількість програм, які реалізують квантово-хімічні методи. Найвідоміші: Gaussian, US-GAMMES, UK-GAMMES, PC-GAMMES, Abinit. 23].

2.2 Теорія функціоналу густини

Найважливіша характеристика електронного основного стану — це його енергія Е3. За допомогою методів хвильових функцій Е можна знайти або безпосередньо шляхом наближеного розв’язання рівняння Шредінгера або використовуючи принцип Релея-Рітца:

(2.2)

де — пробна хвильова функція основного стану. Цей принцип можна сформулювати на мові електронної густини. При фіксованій густині можна визначити локальний мінімум енергії таким чином:

. (2.3)

Тут T і U — кінетична и потенціальна енергія відповідно, і

. (2.4)

кристал фононний стан густина Цей функціонал не вимагає точного знання потенціалу. Далі можна мінімізувати енергію:

. (2.5)

У разі невиродженого основного стану мінімум досягається тоді, коли є густиною основного стану, а для виродженого основного стану — коли є одна з густин основного стану (будь-яка). Отже, складна задача пошуку мінімуму (2.1) по ЗN-мірним пробним функцій (N — число частинок в системі) зведена до достатньо тривіальної за своїм виглядом задачі знаходження мінімуму по 3-мірним пробним функцій. Насправді визначення (2.3) для приводить нас до мінімізації по ЗN-мірним пробним хвильовим функціям. Подальший розвиток цього формалізму в дусі наближення Хартрі призводить до узгоджених рівнянь, які називаються рівняннями Кона-Шема:

(2.6)

де — так званий обмінно-кореляційний потенціал, і - одночасткова енергія і одночасткова хвильова функція (для j-ої частинки), а густина визначається наступним чином:

(2.7)

(2.8)

Потрібно, однак, відзначити, що до цих пір існує проблеми вибору обмінно-кореляційного потенціалу і вибору базисних хвильових функцій. Ніяких особливих обмежень на форму рівняння (2.5) не ставляться, тому широке поширення одержав спосіб вибору цього локального потенціалу у вигляді потенціалу для однорідного електронного газу — так зване наближення локальної густини (local density approximation, LDA):

(2.9)

де — обмінно-кореляційна енергія однорідного електронного газу, а — радіус сфери, що припадає на один електрон, визначається як .

Це наближення, широко поширене при дослідженні різних об'єктів, тим не менше, погано працює для ізоляторів і деяких інших систем і в даний час використовується, в основному, як початкове наближення. Зробимо деякі зауваження:

1. Наближення локальної густини з високою точністю описує типові метали, а також з достатньо хорошою точністю описує системи з високою електронною густиною, наприклад, перехідні метали.

2. Розв’язання рівнянь Кона-Шема в наближенні локальної густини не набагато складніше, ніж розв’язання рівнянь Хартрі, і набагато простіше, ніж розв’язання рівнянь Хартрі - Фока. При цьому точність розрахунку енергії обміну становить порядку 10%, в той же час звичайно менша по величині кореляційна.

3. З практики розрахунків відомо, що ПЛП дає енергії іонізації атомів, енергії дисоціації молекул і енергії зв’язку твердих тіл з непоганою точністю, звично 10 — 20%. При цьому довжини зв’язків і, отже, геометрична будова молекул і твердих тіл виходять в ПЛП, як правило, з дуже хорошою точністю ~1 — 2%.

4. Наближення локальної густини, як і наближення густини (узагальнення ПЛП для систем з неспареними спінами) спіну, може виявитися непридатним в деяких випадках, наприклад, для систем з важкими ферміонами, коли кореляційні ефекти електрон-електронного взаємодії (для d-, f-електронів) такі сильні, що ці системи втрачають всяку схожість з невзаємодіючим електронним газом.

Рис. 2.1 — Ілюстрація ідеї наближення локальної густини. Розподіл електронного газу в деякому нескінченно малому об'ємі можна вважати однорідним, таким чином, його внесок в обмінно-кореляційну енергію по величині рівний внеску однорідного електронного газу такого ж об'єму. Вісь абсцис пропорційна радіусу сфери, що доводиться на один електрон. По осі ординат відкладена обмінно-кореляційна енергія однорідного електронного газу енергія істотно завищується, як правило, приблизно в 2 рази. В більшості випадків обидві помилки частково скорочуються Алгоритм розрахунку за допомогою теорії функціонала локальної густини, що включає самоузгоджене визначення електронних потенціалів і густини наведено на наступній блок-схемі (2.2).

Вхідні дані: Як вхідні дані задаються координати атомів, заряди ядер, повне число електронів. У разі використання методу псевдопотенціалів

Рис. 2.2 — Блок-схема алгоритму розрахунку за допомогою теорії функціоналу густини.

задається також форма псевдопотенціалу для всіх різних типів атомів в системі. Якщо як базисні функції використовуються плоскі хвилі, необхідно перевіряти результати розрахунку, оскільки вони залежать від багатьох параметрів: енергії обрізання, числа к-точок при інтегруванні по зоні Бріллюена та ін.

Проблема вибору базисних функцій, які використовуються при розв’язанні рівнянь Кона-Шема полягає в тому, що їх число повинне бути, з одного боку, досить великим, щоб якомога ближче описувати точний базис (який в загальному випадку нескінченний). З іншого боку, число базисних функцій має бути, по можливості, найменшим, щоб забезпечити стабільність розрахунку і збільшити його швидкість. Зазвичай для вирішення тих чи інших задач використовуються готові набори базисних функцій, властивості яких відомі (http://www.wfu.edu/natalie/papers/pwpaw/ periodictable / periodictable. html).

У програмному коді Abinit реалізований вибір базису у вигляді приєднаних плоских хвиль. Кожен атом при цьому оточується деякої сферою (атомною сферою), всередині якої хвильова функція шукається у вигляді розкладу по базису плоских хвиль так, щоб поза сферою можливо було «зшити» цю функцію з функцією вільного електрона. У програмі Abinit, крім того, використовується метод псевдопотенціалу, що дозволяє замінити частину базисних функцій на деякий ефективний потенціал. Така заміна проводиться автоматично при розрахунку. Від користувача потрібно лише вибрати псевдопотенціал необхідного типу.

2.3 Наближення теорії псевдопотенціалів

Стандартні алгоритми діагоналізації дозволяють розрахувати фіксовану кількість власних векторів гамільтоніану за час, що пропорційний N3, де N — кількість базисних функцій гамільтоніан. Очевидно, треба прагнути до того, щоб власні вектори вимагали для свого розкладу мінімальну кількість базисних векторів або, застосовано до нашого випадку розрахунок хвильових функцій електронів, — членів розкладу хвильових функцій електронів по базисним функціям. Для зменшення цієї кількості членів розкладу прагнуть до того, щоб використовувані базисні функції були максимально схожими на власні хвильові функції електронів в даній молекулярній або кристалічній структурі. Одним із можливих наборів базисних функцій є набор плоских хвиль, головними перевагами якого є властивість повноти та ортонормованості. На жаль, його пряме застосування для опису хвильових функцій електронів обмежено тим, що хвильові функції поблизу ядер атомів сильно осцилюють, тому для їх розкладу по плоским хвилям необхідно використати величезну кількість членів розкладу.

Суть теорії псевдопотенціалу полягає в тому, щоб шляхом деякої процедури зменшити степінь осциляцій розрахованих валентних хвильових функцій поблизу ядра атома. При цьому враховуються тільки валентні електрони, оскільки відомо, що більшість фізичних властивостей системи залежить саме від поведінки валентних електронів. Внутрішні електрони вважаються незмінними, тобто поведінка їх хвильових функцій не змінюється при зміні зовнішнього хімічного оточення атома. Ці електрони будуть приводити тільки до зміни ефективного заряду іона. Ця процедура є еквівалентною заміні сильного електрон-іонного потенціалу більш слабким псевдо потенціалом, який визначає всі явно виражені властивості валентних електронів, в тому числі й релятивіських. Таким чином, досліджувану систему замінюють системою, що складається із псевдовалентних електронів і псевдоіонів. Властивості псевдопотенціолу такі, що його потенціал поза деяким радіусом обрізки rc співпадає з потенціалом модельованого іону, а всередині цієї сфери набагато слабший. В цьому випадку рівняння Шредінгера розв’язується набагато легше всередині сфери радіуса rc, оскільки шукана хвильова функція розкладається по меншій кількості базисних функцій. Для того, щоб сконструювати псевдопотенціал вводять два припущення:

всі електрони розбивають на дві групи: інертні електрони остову та валентні електрони, при цьому стани остову в системі залишаються такими ж як у вільному атомі. Таким чином, теорія псевдопотенціалу трактує речовину як суму іонних остовів, що складається із ядра атому та інертних електронів серцевини, та валентних електронів, що знаходяться між атомами;

замінюють багато електронну задачу задачею про самоузгоджене поле, де ефекти взаємодії між електронами враховуються через самоузгоджений потенціал, який враховує також обмінні властивості та електронну кореляцію. Найбільш поширені наступні псевдопотенціали:

Псевдопотенціал Керкера. У 1980 році ним була запропонована схема побудови псевдопотенціалу, що зберігає нормування. Він запропонував безвузлову псевдохвильову функцію експоненціального вигляду з аргументом у вигляді полінома четвертого ступеня, що має коректну поведінку при великих r. А саме, у даному псевдопотенціалі електронна хвильова функція усередині радіусу обрізки rC замінюється аналітичною функцією наступного вигляду:

(2.10)

Функція f® вибирається так, щоб підсумковий потенціал був гладким, що не має особливостей. Одна з можливих форм цієї функції наступна:

(2.11)

де (2.12)

Щоб уникнути сингулярності в псевдопотенціалі, в (2.39) відсутній лінійний член. Для визначення коефіцієнтів в поліномі p® використовуються умови збереження нормування.

Псевдопотенціал Хаманна, Шлютера і Чанга. У цій схемі спочатку розраховується повний електронний потенціал для ізольованого атома. Цей потенціал множиться на гладку короткодіючу хвильову функцію, що зменшує вплив сильнодіючої та сингулярної частин потенціалу, яка пов’язана з притяганням. Псевдохвильова функція має вигляд:

(2.13)

де функція рівна нулю при і поводиться пропорційно xl+1 при малих x. Вибрані асимптоти f (x) і g (x) гарантують, що потенціал буде скінченим при r=0. д1 — найменший корінь квадратного рівняння, тоді умова нормування псевдохвильової функції має вид:

(2.14)

Кінцевий псевдопотенціал відтворюється з псевдохвильової функції через процедуру інвертування радіального рівняння КонаШема.

Псевдопотенціал Вандербільта. Цей псевдопотенціал був запропонований Вандербільтом і незалежно від нього Блохом у 1990 р. Важливою його особливістю є послаблення умови збереження нормування. Псевдохвильова функція при r < rC представляється настільки згладженою, наскільки це взагалі можливо, тим самим значно зменшуючи необхідну енергію відриву. Нестача заряду, що виходить в результаті, компенсується введенням приєднаних атомно-центрових зарядів. Приєднані заряди визначаються як різниця електронної густини між повною і псевдохвильовою хвильовою функціями. Радіус обрізки псевдопотенціалу не залежить від максимуму хвильової функції, а вибирається на половині відстані між найближчими сусідами. Проте для приєднаних зарядів вибирається маленький радіус обрізки для точного відновлення розподілу моментів і заряду повної хвильової функції.

Псевдопотенціал Вандербільта має широку область застосувань. Головним чином він використовується у випадку перехідних металів, дозволяючи проводити розрахунки з високими швидкістю і точністю. Проте у методу є недоліки, пов’язані, головним чином, з складною схемою конструювання псевдопотенціалу, оскільки необхідно враховувати достатньо багато параметрів.

Псевдопотенціал Труллєра-Мартінса. Відмінність схеми конструювання цього псевдопотенціалу від схеми Керкера полягає у тому, що поліном в експоненті для схеми Труллєра-Мартінса є наступним:

(2.15)

Коефіцієнти підібрані таким чином, щоб потенціал підпорядковувався критеріям збереження нормування. Крім того, вимагається, щоб псевдохвильова функція мала чотири безперервні похідні в інтервалі r

РОЗДІЛ 3. РЕЗУЛЬТАТИ AB INITIO РОЗРАХУНКІВ КРИСТАЛІВ Sn2P2S6

3.1 Рівноважна геометрична структура

3.1.1 Методика розрахунку геометричної структури

Розрахунки просторової структури досліджуваних кристалів Sn2P2S6 проведено в рамках теорії функціоналу густини з використанням програмного пакету ABINIT [23, 24]. Ця програма заснована на базисі плоских хвиль та псевдопотенціалів і використовує ефективний алгоритм швидкого перетворення Фур'є (для перетворення хвильових функцій між реальним та оберненим простором) [25], зонний метод спряженого градієнту та алгоритм спряженого градієнту потенціала, що дозволяє визначити самоузгоджений потенціал [27], повну енергію та оптимізувати геометричну конфігурацію системи. Для опису електронів атомного остову використано псевдопотенціали, розроблені Трулльє і Мартінсом. Необхідний базисний набір плоских хвиль задавався параметром обмеження максимальної кінетичної енергії - 15 Ry. Самоузгодження по матриці густини проводилося до досягнення величини 10−10 а.о.е. для різниці повних енергій елементарної комірки між двома послідовними ітераціями. Інтегрування по зоні Бріллюена проводилося по спеціальним точкам типу Монкхорста-Пака. В обидвох структурних модифікаціях для кристалу Sn2P2S6 вико-ристовувався набір спеціальних точок 4?4?4, тобто нескінченний кристал замінювався циклічною системою із 64 примітивних комірок (такому набору відповідає 16 точок в неприведеній частині зони Бріллюена). Оптимізація атомних координат і параметрів елементарної комірки була проведена з використанням методу BFGS і продовжувалася до тих пір, поки сили, що діють на атоми не стали меншими за 2,5· 10−3 еВ/?, а механічні напруги в комірці меншими за 1,5· 10−3 ГПа.

3.1.2 Рівноважна геометрична структура

Як відомо [31], елементарна комірка кристалу Sn2P2S6 в парафазі містить дві формульні одиниці. При пониженні температури кристал зазнає фазового переходу 2 роду із вихідної параелектричної фази (P21/n) в сегнетоелектричну (Pn). Виходячи з цього структурну оптимізацію геометрії проведено із дотриманням симетрії. Результати пошуку рівноважної гео-метрії у вигляді залежності енергії елементарної комірки від її об'єму зображено на рис. 2. Як видно з рисунку в результаті проведених розрахунків отримано дві стабільні конфігурації: сегнетофазу, повна енергія елементарної

Рис. 3.1 — Залежність енергії елементарної комірки від її об'єму в процесі оптимізації геометрії [33]

комірки якої Е=-182,5369Ха (-4967,082еВ) та парафазу Е=-182,5365Ха (-4967,073еВ). Як і передбачалося, менш симетрична низькотемпературна кон-фігурація кристалу Sn2P2S6 (сегнетофаза) є більш стійкою. Порівняння осно-вних розрахованих кристалографічних параметрів гратки з відповідними експериментальними [32, 33] наведено в таблицях 3.1 і 3.2. Відхилення оптимізованих значень постійних гратки для кожної із кристалічних модифікацій від експериментальних становить 3−6%. В таблиці 2 наведено координати тільки тих симетрично нееквівалентних атомів елементарної комірки, які займають позиції Уайкова з вільними параметрами. Координати інших атомів однозначно визначаються заданням просторової групи симетрії та координат атомів зайнятих позицій Уайкова.

Табл. 3.1 — Розраховані та експериментальні [29, 30] (наведено в дужках) значення відносних координат симетрійно нееквіваленних атомів Співставлення розрахованих та експериментальних значень [31, 32] дозволяє стверджувати, що відносні координати атомів в сегнетота парафазах кристалу Sn2P2S6 відтворюються з хорошою точністю (відхилення становить 2−5%).

3.2 Фононні спектри кристалів Sn2P2S6

3.2.1 Методика розрахунку фононних спектрів

Отримавши оптимізовану структуру кристалів Sn2P2S6 наступним кроком було обчислення фононних спектрів досліджуваних кристалів, вико-ристовуючи при цьому оптимізовані координати елементарної комірки кристалів Sn2P2S6. Крім оптимізованих координат та енергії обрізки, який для обчислення фононних спектрів був рівний 35 Ry точність обчислення фононних частот із перших принципів контролюється таким параметром, як число спеціальних точок, що використовуються для підсумовування по ЗБ. В наших розрахунках використовувалася k-сітка 4?4?4 (із зсувом від початку координат), що відповідає 11 точках у непроводимій зоні ЗБ. Метод Фур'є-інтерполяції дозволив за допомогою розрахунку динамічних матриць у деякому наборі точок ЗБ отримати силові сталі в реальному просторі і потім використовувати їх для розрахунку фононних мод в довільній точці ЗБ. Якість отриманої таким чином дисперсійної залежності обумовлена тільки числом точок так званої q-сітки. В наших розрахунках, так само як і при інтегруванні по ЗБ, використовувалася високоефективна q-сітка 4?4?4 (без зсуву від початку координат).

Як показали обчислення з енергією обрізки при 50 Ry, описані вище параметри забезпечують дуже хорошу збіжність фононних частот в точці Г.

3.2.2 Розрахунок симетричних елементів

У вихідному log файлі, який формується під час обчислення фононних спектрів досліджуваних кристалів наявна інформація про симетричні елемен-ти кристалів Sn2P2S6. Зокрема маємо такі розрахунки:

Для Sn2P2S6 (cегнетофаза):

symspgr: the symmetry operation no. 1 is the identity symspgr: the symmetry operation no. 2 is an inversion symaxes: the symmetry operation no. 3 is a 2₁-axis symplanes: the symmetry operation no. 4 is an n plane symspgr: spgroup= 14 P2₁/c.

Для Sn2P2S6 (парафаза):

symspgr: the symmetry operation no. 1 is the identity

symplanes: the symmetry operation no. 2 is an n plane

symspgr: spgroup= 7 Pc (=Cs²)

3.2.3 Аналіз симетрії фононних мод

Для аналізу симетрії фононних мод була використана утиліта ANADDB (Analysis of Derivative DataBase), яка дозволяє також обчислювати безліч різних фізичних властивостей, що містять похідні повної енергії (похідна DataBase — DDB). Вихідний файл містить наступну інформацію:

Symmetry characters of vibration mode # 1

degenerate with vibration modes # 2 to 3

3.0 -3.0 -1.0 1.0 => Au+2Bu

Symmetry characters of vibration mode # 4

1.0 -1.0 -1.0 1.0 => Bu

Symmetry characters of vibration mode # 5

1.0 1.0 1.0 1.0 => Ag

…

Symmetry characters of vibration mode # 60

1.0 -1.0 1.0 -1.0 => Au

В результаті маємо для акустичних і оптичних мод в сегнетофазі неприводиме представлення в точці Г:

Аналогічно для парафази маємо такий результат:

Отримані результати співпадають з теоретико-груповим аналізом, проведеним в [x].

3.2.4 Фононні моди в цетрі ЗБ

В табл. 3.2 і 3.3 для прикладу наведено порівняння розрахованих нами фононних мод (сегнетофаза — по симетрії Ag, парафаза — по симетрії A") в центрі зони Бріллюєна з існуючими теоретичними та експериментальними результатами. Порівняння теоретичних і експериментальних частот показує в цілому задовільне узгодження між трьома наборами даних. Зокрема слід відмітити, що наші розрахунки краще відтворюють експериментальні результати ніж відповідні розрахунки з. Теоретичні частоти найбільш близькі до експериментальних в нижній і середній частині спектру, тоді як для коливальних мод у верхній частині спектру спостерігається заниження до 20 cm-1 в порівнянні з відповідними теоретичними і експериментальним результатами.

Табл. 3.2 — Фононні спектри (Сегнетофаза — симетрія Ag)

Табл. 3.3 — Фононні спектри (Парафаза — симетрія A")

Тому питання про ступінь близькості отриманих нами результатів до експериментальним вимагає більш ретельного аналізу та порівняння з можливо великою її кількістю. Тим не менш, систематичне завищення або заниження фононних частот при розрахунках в рамках теорії функціонала щільності раніше вже зазначалося і являється характерною особливістю розрахунків із перших принципів. Як показано в [24], причиною завищення фононних мод є використання наближення локальної густини (узагальненого градієнтного наближення) при описанні обмінно-кореляційної взаємодії.

Висновки

Проаналізовано літературу по темі динаміки гратки сегнетоелектриків Sn2P2S6, а також освоєно методику квантовохімічних розрахунків, що реалізована в програмному пакеті ABINIT.

Розраховані в процесі оптимізації геометрії значення постійних гратки та координат атомів в парата сегнетофазі кристалу практично відтворюють експериментальні результати. На основі результатів аналізу енергетичної стабільності структурних модифікацій кристалу Sn2P2S6, встановлено, що ацентрична конфігурація (сегнетоелектрична фаза) досліджуваної сполуки є більш стабільною.

Обчислені в центрі зони Брюлльена фононних мод та їх симетрії. Деякі завищення частот фононних мод зумовлене використанням наближення локальної густини (узагальненого градієнтного наближення) при описані обмінно-кореляційної взаємодії.

Література

1. Dittmar G., Schafer H. Die Struktur des Di — Zinn — Hexathiohypodiphosphats Sn2P2S6 //Zs. Naturforsch. 1974. Vol. 29 b. N 5−6. S. 312−317.

2. Carpentier C.D., Nitsche R. Vapour growth and crystal data of the thio (seleno) — hypodiphosphates Sn2P2S6, Sn2P2Se6, Pb2P2S6, Pb2P2Se6 and their mixed crystals //Mat. Res. Bull. 1974. Vol. 9. № 4. P. 401−410.

3. Ризак И. М. Динамики решетки и термодинамика сегнетоэлектриков системы Sn (Pb)2P2S (Se)6. Автореф. Дис. … канд. физ.-мат. наук. Ужгород. 1992.

4. Ризак В. М., Высочанский Ю. М., Грабар А. А., Сливка В. Ю. Несоразмерная фаза сегнетоэлектриков системы Sn2P2S6 — Sn2P2Se6 в модели жестких ионов //Физика твердого тела. 1989. Т. 31. № 7. C. 154−160.

5. Высочанский Ю. М., Фурцев В. Г., Хома М. М. Расщепление сегнетоэлектрического фазового перехода в поле лазерного излучения и ого самофокусировка //Журн. эксперимен. и теорет. физики 1985. T. 89. № 9. C.939−945.

6. Ковалев О. В. Неприводимые представления пространственных групп. К., 1961.

7. Багавантам С., Венкатарайду Т. Теория групп и ее применение к физическим проблемам. М., 1959.

8. Fateley W.G., Mc Devitt N.T., Bentiey F.F. Infared and Raman selection rules for lattice vibrations: The correlation method //Appl. Specrosc. 1971. Vol. 25. № 2. P. 155−173.

9. Высочанский Ю. М., Сливка В. Ю., Ворошилов Ю. В. и др. Поляризационные спектры комбинационного расеяния сегнетополупроводника //Физика Твердого тела. 1979. Т. 21. № 1. С. 211−215.

10. Burger V.H., Falius H. Schwingungsspectren und Kraftkonstanten von PS4−3 P2S6−4 und isoelektrichen Molekeln SiCl4 und SiCl6 //Zs. anorg. allg Che 1968. Vol. 363. S. 24−32.

11. Genzel L., Martin T.P., Perry C.H. Model for long wavelength opticalphonon modes of mixed crystal //Phys. Status Solidi (b). 1974. Vol. 62. № 1. P. 83−92.

12. Грабар А. А., Высочанский Ю. М., Мельник Н. Н., Субботин С. И., Панфилов В. В., Сливка В. Ю. Влияние гидростатического давления на колебательные спектры сегнетоэлектрика Sn2P2S6 //Физика твердого тела. 1984. Т. 26.№ 1. С. 65−68.

13. Киттель Ч. Квантовая теория твердых тел. М., 1967. 491 с.

14. Рейсленд Дж. Физика фононов. М., 1975. 365 с.

15. Valevichius V., Samulionis V., Skritskij. Orientatinal dependence of ultrasonic velocity near the phase transition in Sn2P2S6 single crystals //Ferroelectrics. 1988. Vol. 79. № 1. P.225.

16. Борн М. Кунь Х. Динамическая теория кристаллических решеток. М., 1958.

17. Venkataraman G., Sahni V. C. External vibrations in complex crystals //Rev. Mod. Phys. 1970. Vol. 42. N 4. P. 409−470.

18. Высочанский Ю. М., Герзанич Е. И., Грабар А. А., Майор М. М., Перечинский С. И., Сало Л. А., Сливка А. Г., Сливка В. Ю., Тягур Ю. И., Хома М. М. Фазовые диаграммы смешанных кристаллов типа Sn2P2S6: феноменологический анализ и динамика решетки // Кристаллография. — 1985. Т.30, № 6. С.1116−1120.

19. Різак І.М., Чобаль О.І., Різак В.М. Abinitio розрахунки електронних і коли-вальних спектрів сегнетоелектричнихкристалів Sn2P2S6 // Науковий ВісникУжгородського університету. Серія Фі-зика. — 2006. — Випуск 19. — С. 26−31.

20. Майор М. М., Бовтун В. П., Поплавко Ю. М. и др. Диэлектрические свой-ства кристаллов Sn2P2S6// ФТТ.- 1984.-Т.26. № 3. С. 659−664.

21. Hohenberg P., Kohn W. Inhomogeneouselectron gas // Phys. Rev. — 1964.-V.136. P.864−871.

22. Gonze X.; Beuken J.-M.; Caracas R.; DetrauxF.; Fuchs M.; Rignanese G.-M.;Sindic L.; Verstraete M.; Zerah G.;Jollet F.; Torrent M.; Roy A.; Mikami M.;Ghosez Ph.; Raty J.-Y.; Allan D.C. Firstprinciplecomputation of material properties: the ABINIT software project// Comput. Materials Science. -2002. — V.25. -No.3. -P. 478−492.