FORTH

Встреча с FORTH

 

Cистемы типа Matlab/Octave можно рассматривать как своего рода калькуляторы, только высокого уровня: в них есть векторы, матрицы — всё, что нужно для реализации научных вычислений. Однако, они большие и требуют мощного компьютера с операционной системой. Интересно, что давным давно, во времена, когда персональные компьютеры были немногим сильнее калькулятора, как, например, мой HP 35s существовало интересное решение.

 

Например, персональный компьютер Jupiter Ace со встроенным языком программирования FORTH. Особенность этого языка программирования в том, что ему не нужна ни операционная система, ни файлы, ровным счетом ничего, кроме чистого железа. Это позволило оснащать довольно слабые по современным меркам компьютеры даже офисным программным обеспечением, как например, компьютер Canon Cat. Набор приложений был записан в ПЗУ. В него входили: стандартный пакет офисных программ, орфографический словарь на 90 000 слов, программа связи и средства программирования на Форт и языке ассемблера.

FORTH использует технологию «шитого кода». Проще говоря, это массив вызовов подпрограмм. Указатель выполнения кода последовательно принимает адреса подпрограмм, которые исполняются. Передача параметров происходит через стек.

Шитый код используется в программировании калькулятора HP 35s

Автомат FORTH функционирует следующим образом. Из входного потока (программы на языке FORTH) читается очередное слово. Синтаксис языка крайне прост: слово — это любая последовательность символов не включающая пробел, табуляцию или перевод строки. Если принятое слово литерал (например, число), то оно кладётся на стек. Если нет, ищется в словаре. Если FORTH в состоянии интерпретации — слово исполняется, если в состоянии компиляции (определении нового слова) то адрес определения слова или подшивается к текущему определению, или немедленно исполняется (если слово имело атрибут IMMEDIATE). Если слова не было в словаре — FORTH выдает ошибку исполнения.

Контрольные слова (для определений ветвлений и циклов) имеют атрибут IMMEDIATE. FORTH позволяет управлять входным потоком слов и определять управляющие конструкции, т. е. слова с инструкциями для состояния интерпретации и компиляции.

Вот пример корректной программы на FORTH:

: НАЛИТЬ-ВОДУ КРАНЫ ОТКРЫТЬ ДО-НАПОЛНЕНИЯ КРАНЫ ЗАКРЫТЬ ; 
: ПОЛОСКАТЬ НАЛИТЬ-ВОДУ СТИРАТЬ ВЫЛИТЬ-ВОДУ ; 
: СТИРАЛЬНАЯ-МАШИНА СТИРАТЬ ВЫКРУЧИВАТЬ ПОЛОСКАТЬ ВЫКРУЧИВАТЬ ;

 Нравится? Читаем книгу Л. Броуди. Начальный курс программирования на языке ФОРТ.

Сразу отмечу, что FORTH -  не способ программирования на естественном языке, это скорее метод формализации естественного языка на строгий алгоритмический язык конечного автомата с несколькими стеками. Простота достигается за счёт применения стеков и обратной польской записи, а также за счёт отсутствия синтаксиса, типов данных и прочих абстракций программирования. Однако, если какие-то абстракции нужны - их всегда можно определить, так как вам угодно, а не создателям языка.

Nick Morgan любезно предоставил всем желающим попробовать FORTH прямо в окне браузера.

Знакомство 

Далее я перечислю места, где в настоящее время встречается FORTH. 

В космических аппаратах!

Солидный список здесь. Отдельная история - чипы RTX2000-RTX2010. Читайте.

Программирование микроконтроллеров: семейства AtmelAVR8 Atmega, Microchip 8-bit PIC18F и 16-bit PIC24, 30, 33 и Arduino UNO и многих других. Читайте тут.

Встроенный скриптовый язык:

Игровая платформа настольных игр на ПК - Zillions of Games, в ней FORTH это способ разрабатывать свои игры, см. пост.

Скриптовый язык в XYPLOT - Extensible Plotting and Data Analysis Program for 32-bit x86 GNU/Linux, в планировщике задач nnCron и других приложениях.

Подводные камни

Обратная польская нотация очень удобна в калькуляторах, однако, она не всегда подходит для решения задач. Попытки избавиться от традиционного подхода с переменными бывают трагичны. В самом деле, как отмечает в своем блоге Yossi Kreinin, для стека подходят деревья выражений




тогда как более общая структура, граф выражения уже требует операций со стеком: перестановок, копирования и т.п.
Сравните выражение в обычной и обратной польской нотации (a + b)*c+7 = a b + c * 7 + и попробуйте записать (a + b)/(a + b*c)! В итоге программы на FORTH начинают выглядеть как жонглирование элементами стека.

: siftDown                             ( a e s -- a e s)
  swap >r swap >r dup                  ( s r)
  begin                                ( s r)
    dup 2* 1+ dup r'@ <                ( s r c f)
  while                                ( s r c)
    dup 1+ dup r'@ <                   ( s r c c+1 f)
    if                                 ( s r c c+1)
      over over r@ precedes if swap then
    then drop                          ( s r c)
    over over r@ precedes              ( s r c f)
  while                                ( s r c)
    tuck r@ exchange                   ( s r)
  repeat then                          ( s r)
  drop drop r> swap r> swap            ( a e s)
;

Пример взят с http://rosettacode.org/wiki/Sorting_algorithms/Heapsort#Forth

Всё это приводит к тому, что язык Си кажется единственным, кто может с этим справиться. Одна из таких трагедий описана у Yossi Kreinin.  

Почему бы не решить проблему переменными? Они есть в FORTH. Кажущаяся проблема в том, что переменные глобальны и не связаны с единственным определяемым словом. В Си-подобных языках избегают глобальных переменных из-за того, что исчерпываются уникальные имена, сложно отслеживать значение переменной, если доступ к ней возможен из любого места программы...

Я написал кажущаяся проблема, потому что на самом деле проблемы нет вовсе. Странно, что об этой особенности FORTH не упоминается в книгах явно и с предупреждением, но глобальные переменные FORTH ведут себя не так, как в Си.

Единственная книга, в которой обсуждается способ взаимодействия определения с окружением это книга Christian Queinnec. Lisp In Small Pieces. Поведение определений в FORTH попадает в класс гиперстатического связывания. Проще показать на примере.

: foo ." foo!" ;
: bar foo ."  bar!" ;
bar foo! bar! ok.
: foo ." bizzle!" ; Redefine foo.  ok.
foo bizzle! ok.
bar foo! bar! ok.

Слово bar состоит из foo и вывода строки ."  bar!". Однако, когда термин foo оказывается переопределён, то поведение bar остаётся неизменным! Определение слова сохраняет его контекст, слова в определении имеют ровно то значение, которое было на тот момент.

Отсюда следует, что нет никакой причины беспокоиться о том, что кончатся имена переменных и так далее. Вам нужна новая переменная my_variable? Так и назовите её, ничего из того, что связано было с ней ранее не изменится. Нет смысла вводить новый синтаксис и правила обращения с локальными переменными (хотя таких попыток в FORTH было много и разных), достаточно перед тем словом, где нужна переменная определить её, а после него определить новую с тем же самым именем.

Подытожим. FORTH это

1. Шитый код
2. Обратная польская запись
3. Гиперстатическое связывание

С этих теоретических понятий уже можно собрать свой FORTH. Каноническая инструкция как это сделать написана Richard W.M. Jones, читайте по ссылке.

FORTH на вкус

Далее я расскажу о своем опыте погружения в античные времена MS DOS и существовавшей для него среды FORTH.

В dosbox'е устанавливаю F-PC 3.60 (читайте C.H. Ting. F-PC: FORTH optimized for IBM-PC // SIGFORTH Newsl. 1(1) (1989) 15-17 DOI:10.1145/382122.382928)

 

Это полноценная IDE с редактором кода

Справочником-помощью

а также отладчиком

Словом, весь джентльменский набор!

Поскольку FORTH интересовал меня с точки зрения математических вычислений, альтернатива Matlab и прочим Scientific Python, то я избавился от жонглирования стеком и режущими глаз мушками вроде f@ - чтение переменной/ f! - запись значения в переменную.

Вот как это работает.

fload sfloat.seq
fload eval.seq

macro: get-name bl word count r@ place ;
macro: $name  r@ count ;
macro: .. pad +place ;

: cell:
  here >r get-name
  " variable *" pad place $name .. 
  "  : :" .. $name .. "  *" .. $name .. "  ! ; " ..
  "  : "  .. $name .. "  *" .. $name .. "  @ ; " ..
  r> drop pad count eval ;

: float:
  here >r get-name
  " fvariable *" pad place $name ..
  "  : :" .. $name .. "  *" .. $name .. "  f! ; " ..
  "  : "  .. $name .. "  *" .. $name .. "  f@ ; " ..
  r> drop pad count eval ;

: cells:  ( n -- ) 0 do cell:  ( name ) loop ;
: floats: ( n -- ) 0 do float: ( name ) loop ;

Теперь если в каком-то определении нужны переменные, то перед ним я пишу

3 floats: x y z

Присваиваю им значения при помощи автоматически созданных слов

0.0e :x 0.0e :y 1.0e :z

и получаю их значения на стек легко и естественно

x x f* y y f*  z z f*  f+ f+ fsqrt

как если бы делал это на калькуляторе. При необходимости легко переопределить имена x, y, z - на работе предыдущих слов это никак не скажется. Гиперстатическое связывание - удобная штука!
Алгоритмы на FORTH получаются компактные, сам язык очень способствует факторизации кода. Например, сортировка:

\ Heapsort is an in-place sorting algorithm with worst case
\ and average complexity of O(n*log[n]). The basic idea is
\ to turn the array into a binary heap structure, which has
\ the property that it allows efficient retrieval and removal
\ of the maximal element. We repeatedly "remove" the maximal
\ element from the heap, thus building the sorted list from
\ back to front.

0 1 2constant endless

defer item@
defer item!
defer greater?

6 cells: arr n p q r s

: maxi :r :q :p :n :arr
  p :s
  q n < if arr q item@ arr s item@ greater? if q :s then
  then
  r n < if arr r item@ arr s item@ greater? if r :s then
  then s ;

: downheap :p :n :arr
  endless do
    arr n p p 2* 1+ p 2* 2 + maxi :q
    q p = if leave then
    arr p item@ arr q item@
    arr p item! arr q item!
    q :p
  loop ;
2 cells: arr n

: heapsort :arr
  arr @ :n
  0 n 2 - 2/ do arr n i downheap -1 +loop
  n 0 do
    arr n i - 1- item@  arr 0 item@
    arr n i - 1- item!  arr 0 item!
    arr n i - 1- 0 downheap
  loop ;

: item-get item @ ;
' item-get is item@
: item-set item ! ;
' item-set is item!
' > is greater?

Стоит обратить внимание на конструкцию

defer word
 
...
 
' method is word 

Которая позволяет отложить определение какого-либо слова в определении термина до того момента, пока оно не станет известено. А вот реализация метода нахождения корня уравнения, её вполне можно сопоставить с таковой в других языках.

\ Using Ridder's method, return the root of a function func
\ known to lie between x1 and x2.

fload vars.seq

defer f(x)
10 floats: x1 x2 xm xnew  fl fm fh fnew  s result

-1.11e30 fconstant UNLIKELY-VALUE
 1.0e-9  fconstant ACCURACY
     50   constant MAXIT
      
: fsign ( r1 -- r1 ) f0< dup not - ifloat ;
: f~ ( r1 r2 r3 -- flag ) f-rot f- fabs f> ;

: root-not-bracketed? ( fl fh -- flag )
  fdup   f0< fover  f0> and
  ( fh ) f0> ( fl ) f0< and or not ;

: ridder ( r1 r2 -- r1 ) :x2 :x1
  x1 f(x) :fl
  x2 f(x) :fh
  fl f0= if x1 exit then
  fh f0= if x2 exit then
  fl fh root-not-bracketed?
  if abort" root must be bracketed in zriddr" then
  UNLIKELY-VALUE :result
  false
  MAXIT 0
  do
    x1 x2 f+ 2.0e f/ :xm  xm f(x) :fm
    fm fdup f* fl fh f* f- fsqrt :s
    s f0= if result true leave then
    fl fh f- fsign xm x1 f- f* fm f* s f/ xm f+ :xnew
    xnew result ACCURACY f~
    if result true leave then
    xnew :result  xnew f(x) :fnew
    fnew f0= if result true leave then
    fnew fsign fm f* fm f= not
    if xm :x1  fm :fl  result :x2  fnew :fh
    else fnew fsign fl f* fl f= not
      if result :x2  fnew :fh then
      fnew fsign fh f* fh f= not
      if result :x1  fnew :fl then
    then
    x1 x2 ACCURACY f~
    if result true leave then
  loop
  if result drop
  else ." zriddr exceed maximum iterations" drop then ;

float: x
\ : func :x x fcos x x x f* f* f- ; \ 0.0-1.0 x = 0.865474
: func :x x x x f* f* 10.0e x x f* f* f- 5.0e f+ ; \ 0.6-0.8 x=0.7346

' func is f(x)
0.6e 0.8e ridder f(x) f. cr 

Не хуже чем Matlab? Учитывая, что это может работать на микроконтроллере, без колдовства с компилятором. На то можно возразить, что в Matlab'е есть средства высокого уровня, а тут где они? Продолжение следует.

Дуализм частица-волна

Луи де Бройлю были известны соотношения $E=h\nu$ (формула Планка) и $E_0=m_0c^2$ (эквивалентность массы и энергии в специальной теории относительности) где $m_0$ – масса покоя, т.е. масса в собственной системе отсчёта. Согласно теории относительности абсолютного времени нет, и у каждой точки пространства оно своё собственное – четвертая координата t. Собственное время, как рассуждал Луи де Бройль, это как личные часы. Они есть как у человека, так и у электрона. Что это значит?
Он предположил, что каждой частице свойственен внутренний периодический процесс (тик-так с частотой $\nu_0$), который служит мерой собственного времени [1]. Этот самый внутренний процесс с частотой в собственной системе отсчета позволяет связать две формулы энергии $$ E_{0}=m_{0}c^{2}=h\nu_{0} $$ Неподвижный наблюдатель воспринимает колебательный процесс частицы летящей со скоростью $v$ в точке $x$ как $$ \psi=\exp(2\pi i\nu t) $$ Переход в систему отсчёта частицы даёт в показателе экспоненты $$ \underset{\nu}{\underbrace{\nu_{0}\sqrt{1-\frac{v^{2}}{c^{2}}}}}\cdot\underset{t}{\underbrace{\frac{-t_{0}+\frac{v}{c^{2}}x}{\sqrt{1-\frac{v^{2}}{c^{2}}}}}}=\frac{v}{c^{2}}x-\nu_{0}t_{0} $$ то есть бегущую волну – волну де Бройля, которая всегда синхронизирована по фазе с внутренним процессом $\exp(2\pi i\nu_{0} t_{0})$: $$ \psi=\exp(2\pi i(\frac{v}{c^{2}}x-\nu_{0}t_{0})) $$ Подставляя формулы для энергии, получаем $$ \psi(x,t)=\exp(2\pi i(\frac{m_{0}c^{2}}{h}\frac{v}{c^{2}}\cdot x-\frac{E}{h}t_{0}))=\exp(2\pi i(\frac{m_{0}v}{h}x-\frac{E_{0}}{h}t_{0}))=\exp(\frac{i}{\hbar}(p\cdot x-E\cdot t)) $$ В таком случае частице, летящей со скоростью v соответствует длина волны $$ \lambda=\frac{h}{p} $$ где $p=mv$ – импульс частицы, а $h$ – постоянная Планка.
Наглядное представление о том, что такое частица-волна де Бройля и как волны синхронизированы с внутренним колебательным процессом было получено недавно в опытах Ива Куде. Лучше один раз увидеть!
Ив Куде исследовал поведение прыгающей капли-«блинчика» силиконового масла на поверхности жидкости в вибрационной ванне. Слой жидкости (силиконового масла) в ванне был толщиной 4 мм. Под действием вибрации на поверхности появляется рябь (волны Фарадея), если вертикальное ускорение при вибрировании больше пороговой величины (4.5g в условиях эксперимента). Однако, в подпороговом режиме вибрации с частотой 50 Гц в ванне 40x40 мм дают иную картину стоячих волн с длиной $\lambda$=6.5 мм. Капля силиконового масла размером 1 мм, аккуратно помещенная на поверхность, прыгает по этим волнам не собираясь останавливаться.

Экспериментальная установка Ива Куде. Левая верхняя часть рисунка – фотографии движения капли по поверхности. Правая верхняя часть – встреча «блинчика» с дифракционной щелью. Внизу – схема установки. 1 – «блинчик», 2 – траектория его движения, волна, 3 – затухающий всплеск от удара «блинчика» о поверхность, 4 – подтопленная стенка-препятствие, из них сделаны дифракционные щели, 5 – вибростенд.

Силиконовая капля, ведомая волнами на поверхности ведёт себя как квантовая частица. Её прыжки - внутренний колебательный процесс синхронизированный с её «волной де Бройля» соответствующей траектории движения «блинчика» 2 на рисунке.
В местах, где глубина масла меньше (подтопленная стенка 4), стоячие волны затухают и поверхность не может подкидывать «блинчик». Эти участки моделируют препятствия для прыгающей капли (дифракционные щели) не разрушая её. Так был воспроизведен знаменитый двухщелевой эксперимент: собранная статистика поведения капли ничем не отличается от таковой у электрона [2].

Статистика отклонений силиконовой капли-блинчика по данным Y. Couder and E. Fort. Single-particle diffraction and interference at a macroscopic scale. Phy. Rev. Lett., 97(15):154101, 2006.

Стоит отметить, что в отечественной науке обсуждались подобные идеи намного раньше, например, в книге томского ученого Б.Н. Родимова «Автоколебательная квантовая механика» [3].

 Литература 

1. Л. де Бройль «Волны и кванты» 93 178–180 (1967)
2. Robert Brady and Ross Anderson. Why bouncing droplets are a pretty good modelof quantum mechanics. 2014.
3. Б.Н. Родимов «Автоколебательная квантовая механика», Томск 1967.

Принцип Паули

Спин электрона

Пускай частица движется по окружности длины $2\pi r$, а через $\vec{r}$ мы обозначим позицию частицы. 

Так как $\psi(\vec{r})$ должна быть однозначной на окружности, то поворот на $2\pi$ не должен изменять функцию, т.е. $$\psi(\vec{r}) = e^{\frac{i}{\hbar}p\cdot2\pi r}=1$$ что означает, что $p\cdot r/\hbar=n$, где n – целое число. Произведение $\vec{l}=\vec{p}\times \vec{r}$ это момент импульса, что влечёт квантование $l_{z}=n\cdot\hbar$. В эксперименте Штерна-Герлаха показано, что наблюдаемый момент импульса атомов серебра (и электронов) равен 1/2. Мы с необходимостью приходим к выводу, что вращению с моментом импульса $s=1/2$ соответствуют $4\pi$ циклические движения. Такие вращения есть в 3-х мерном пространстве! Они обладают замечательным свойством, что такое вращение не перекручивает подсоединенные к телу верёвки. Трюк? В это сложно поверить не убедившись, потому вот анимация.

Принцип Паули

Перестановка двух связанных объектов в трехмерном пространстве эквивалентна вращению одного из них на $2\pi$. А так как частицы со спином $4\pi$ периодические, то это означает смену знака: первое вращение на $2\pi$ даёт -1, второе тоже -1, а $4\pi$ вращение $-1\times -1 = 1$.

  Belt trick for the exchange of two fermions from Antonio Martos de la Torre on Vimeo.

Dirac's belt trick for spin 1/2 particle from Antonio Martos de la Torre on Vimeo.

Обязательную смену знака волновой функции при перестановке двух частиц со спином S=1/2 можно трактовать как требование не перекручивать связывающую их ленту-пространство [1]. Эта лента гибкая, но до определенных пределов. Далее мы установим этот предел гибкости пространства на скручивание.

Следующий пример описывается в статье В. Вайскопфа [2]. Пусть есть два электрона с волновыми функциями $$ \psi(x_{1})=e^{\frac{i}{\hbar}p\cdot x_{1}},\quad\psi(x_{2})=e^{\frac{i}{\hbar}p\cdot x_{2}} $$ Забудем пока о электростатическом отталкивании, просто летят два электрона с импульсами $p_{1}=-p_{2}$ друг другу на встречу. Расстояние между ними равно $x=x_{2}-x_{1}$. Волновая функция системы электронов это произведение $$ \Psi(x_{1},x_{2})=e^{\frac{i}{\hbar}p\cdot x_{1}}e^{-\frac{i}{\hbar}p\cdot x_{2}} $$ однако, оно не обладает свойством антисимметрии. Легко поправить дело: $$ \Psi=e^{\frac{i}{\hbar}p\cdot(x_{2}-x_{1})}-e^{-\frac{i}{\hbar}p\cdot(x_{2}-x_{1})} $$ и это выражение можно записать как: $$ \Psi=2i\cdot\sin(\frac{p\cdot x}{\hbar}) $$ Плотность вероятности для этого состояния имеет вид $$ \rho(x)=4\cdot\sin^{2}(k\cdot x) $$ где $k = p/\hbar$ - волновое число. Если у нас возможны импульсы от максимального $4/3p_0$ и до нуля, со средним $p_0$, то плотность вероятности $\rho(x)$ представляет собой интеграл по всем волнам с различными значениями волновых чисел $k$. В итоге она имеет вид ступени:

начиная с какого-то минимального расстояния плотность вероятности встретить электроны на расстоянии менее $x_{min}$ стремится к нулю. Минимальное возможное расстояние между ними, как видно имеет порядок их средней длины волны, т.е. $$ x_{min}\approx 1/k_0 $$ Таким образом, можно считать, что силы упругости пространства не позволяют идентичным электронам подходить друг к другу ближе, чем на некий характерный размер, как будто электрон это упругий твердый шарик.

Квантовые числа

Расположим шарики-электроны плотной упаковкой, т.к. положительно заряженное ядро стягивает их к себе, а принцип Паули и кулоновское отталкивание мешают подходить им близко друг к другу. Рассмотрим образовавшуюся структуру [3].

Оболочечная структура атома: главное квантовое число, орбитальное число и магнитное число.

Мы видим, что позиция каждого шарика-электрона задаётся с помощью четырёх квантовых чисел: главного, орбитального, магнитного и спинового. Стремление к заполнению оболочки (правило 8 электронов) есть ничто иное, как попытка атома выстроить максимально плотную и симметричную структуру.

Теория двойного квартета Линнета

В 1961 году Линнет выдвинул интересную модификацию правила октета Льюиса [4, 5]. Он предположил, что ключевым принципом построения оболочки должно быть максимальное отталкивание электронов одного спина. Учитывая, что они же стягиваются ядром, получается плотная упаковка - тетраэдр. Устойчивой оболочкой Линнет считал ориентацию двух тетраэдров, обеспечивающую их максимальное отталкивание, т.е. 4+4=8. Пока мы не рассматриваем спин, его правило не отличается от правила октета, однако, оно приводит к геометрической трактовке связи. Например, отличие однократной, двойной и тройной связей выглядит так:

Интересно, что предсказываемые соотношения для длин связей находятся в прекрасном согласии с наблюдаемой геометрией молекул. Более того, его принцип позволяет объяснить электронную структуру молекулы кислорода, для которой основное состояние - триплет, два неспаренных электрона. Правило октета в данном случае бессильно.

Электронная структура кислорода. Черные точки - электроны со спином вверх, белые - со спином вниз.

Долгое время было загадкой, почему две молекулы NO (свободные радикалы) не образуют устойчивый димер, в отличие от CN, которые легко димеризуются в дициан $(CN)_2$. С позиции теории двойного квартета, структура O=N-N=O потребовала бы пространственного совмещения тетраэдров электронов разного спина, что невыгодно с точки зрения кулоновского отталкивания. Напротив, дициан позволяет минимизировать отталкивание электронов.
Принцип плотной упаковки электронов-шариков описывает все типы химических связей: ковалентные кратные связи, ионные и металлическую связь. Последняя представляет собой частный случай ионной связи. В роли анионов выступают просто электроны. Расположение вершин тетраэдров (часто они искаженные) можно найти с помощью метода локализованных орбиталей Бойза [6].

Силовое поле электронов

Как известно, принцип плотнейшей упаковки хорошо соблюдается для высоко симметричных кристаллов. Атом, как многоэлектронная система с выделенным центром, как мы уже видели, хорошо описывается плотной упаковкой. Интересно, что если рассчитать энергию связи электронов в атоме исходя не из сложного уравнения Шредингера, а из существенно более простой задачи оптимального расположения шариков в некоем силовом поле (электростатика + отталкивание Паули), то результат оказывается в количественном согласии с квантовой теорией [7].

Полная энергия связи электронов с ядром у атомов: сплошные круги - метод силового поля, пустые окружности - квантовомеханический расчёт методом Хартри-Фока.

Потенциал Паули в работе [7] выбран в форме $$ V(r,p)=\frac{\xi^{2}}{4\alpha r^{2}\mu}\cdot\exp\left\{ \alpha\left[1-\left(\frac{r\cdot p}{\xi}\right)^{4}\right]\right\} $$ где $r$ - межэлектронное расстояние (для электронов одного спина!), $p$ - межэлектронный импульс, $\mu=0.5$ - приведённая масса двух электронов, а $\xi=2.767$ и $\alpha=5.0$ - параметры потенциала, подобранные для воспроизведения решения атома водорода и электронного газа.
Сведение квантовомеханической задачи к молекулярной механике фермионов открывает в перспективе легкий и очень быстрый способ моделирования свойств материалов и химических реакций на обычном персональном компьютере. Однако, простые аппроксимации потенциала Паули хорошо работают только для атомов, с более сложными структурами еще предстоит исследовательская работа. Обнадеживающие результаты и выявленные проблемы найдены в работе Julius T. Su [8].

Литература

1. Hartung, R.W. Pauli principle in Euclidean geometry. American Journal of Physics, 47(10) (1979) 900–910.
2. Вайскопф, В. Современная физика в элементарном изложении. УФН 103(1) (1971) 155-179.
3. Stevens, P.S. A Geometric Analogue of the Electron Cloud. Proceedings of the National Academy of Sciences 56(3) (1966) 789-793.
4. Дей К., Селбин Д. Теоретическая неорганическая химия. М.: "Химия" 1976. Стр. 197.
5. Luder, W.F. The electron repulsion theory of the chemical bond. I. New models of atomic structure. Journal of Chemical Education, 44(4) (1967) 206.
6. Duke, B.J. Linnett's double quartet theory and localised orbitals. Journal of Molecular Structure: THEOCHEM, 152(3-4) (1987) 319-330; Foster, J.M. and Boys, S. Canonical configurational interaction procedure. Reviews of Modern Physics, 32(2) (1960) 300.
7. Lawrence, W. and Cohen, J.S. Fermion molecular dynamics in atomic, molecular, and optical physics, Contemporary Physics, 39(3) (1998) 163-175.
8. Julius Su and William A. Goddard III (advisor) An electron force field for simulating large scale excited electron dynamics. (2008).

Термодинамические потенциалы и преобразование Лежандра

Свойства термодинамической системы называются параметрами состояния. Параметр состояния – некая величина, характеристика системы, изменение которой определяется только начальным и конечным состоянием системы, т.е. $$ \Delta V=V_{2}-V_{1} $$ не зависит от характера процесса изменения его состояния (путь от $V_{1}$ к $V_{2}$). Параметры делятся на две группы: 1) экстенсивные – зависящие от количества вещества и обладающие свойством аддитивности, т.е. $$ V=V_{A}+V_{B} $$ для двух систем A, B (самый простой пример – объём V) и 2) интенсивные – не зависящие от количеств вещества. Интенсивные параметры обладают следующим важным свойством. Они определяются контактным равновесием.
Пусть система поделена две части A и B, объёмов $V_{A}$ и $V_{B}$ находящиеся в контакте. Условием равновесия двух объёмов $V_{A}$ и $V_{B}$ является равенство давлений. Это несложно себе представить: подвижный поршень разделяет две камеры A и B.
Исследовать систему, прибегая к контактным измерениям, например, приложив термометр, значительно проще, чем при помощи измерения и контроля экстенсивных параметров. Так, не существует прибора, при помощи которого можно непосредственно измерить энтропию и нет приспособления для удерживания её постоянной. В этом главная проблема при использовании фундаментального уравнения $$ dU=TdS-PdV $$ Следовательно, оно должно быть преобразовано таким образом, чтобы работать не с функцией состояния экстенсивных параметров $U(S,\,V)$, а с другой функцией, также сохраняющей полное описание системы, но имеющей зависимость от интенсивных параметров.
Сделаем это с помощью преобразования Лежандра.
Пусть есть функция двух независимых переменных $f(x,y)$. Тогда её полный дифференциал $$ df=\underset{u}{\underbrace{\left(\frac{\partial f}{\partial x}\right)_{y}}}dx+\underset{w}{\underbrace{\left(\frac{\partial f}{\partial y}\right)_{x}}}dy $$ где пары ($u=\left(\frac{\partial f}{\partial x}\right)_{y}$, $x$) и ($w=\left(\frac{\partial f}{\partial y}\right)_{x}$, $y$) называются сопряженные пары переменных. Запишем дифференциал произведения $$ d(w\cdot y)=ydw+wdy $$ и вычтем его из исходного выражения для $df(x,y)$ $$ d(f-w\cdot y)=udx+wdy-ydw-wdy=udx-ydw $$ Полученное выражение, в свою очередь, является полным дифференциалом новой функции $g(x,w)$, описывающий ту же самую систему, что и исходная функция$f(x,y)$. $$ g=f-w\cdot y $$ Переход от $f(x,y)$ к $g(x,\left(\frac{\partial f}{\partial y}\right)_{x})$ называется преобразованием Лежандра функции $f(x,y)$. Легко убедиться в том, что оно обратимо.
Воспользуемся преобразованием Лежандра для выведения новой функции состояния вместо $U(S,V)$. На первом шаге из $dU=TdS-PdV$ получаем $$ H(S,P)=U+P\cdot V $$ $$ dH=TdS+VdP $$ что есть ничто иное, как энтальпия. Наконец, подвергнув преобразованию Лежандра $H(S,P)$ найдём свободную энергию Гиббса $$ G(P,T)=H-T\cdot S $$ Таким способом из фундаментального уравнения термодинамики получаются все термодинамические потенциалы.

Взаимосвязь термодинамики и статистической физики

Пусть есть два тела, большое – назовем его резервуар, настолько большое, что оно совсем не остынет если к нему присоединить совсем маленькое тело - систему. Такая пара называется в статистической механике канонический ансамбль. Пусть это маленькое тело состоит из двух независимых подсистем A и B. Для независимых систем энергии складываются $$ E=E_{A}+E_{B} $$ а вероятности реализации той или иной конфигурации умножаются $$ P(E)=P(E_{A})\cdot P(E_{B}) $$ где $E_{A}$, $E_{B}$ – энергии, а $P(E_{A})$ и $P(E_{B})$ – вероятности нахождения подсистем A в конфигурации с энергией $E_{A}$ и $E_{B}$, соответственно. Связать суммирование энергий и умножение вероятностей можно единственным способом: $$ P(E)=C\cdot e^{-\beta E} $$ Эта функция определяет вероятность найти систему в одном из состояний с энергией равной E. Важно следующее. Это не вероятность того, что система имеет энергию E, так как может оказаться несколько конфигураций с такой энергией. Стремление системы к минимальной энергии реализуется как увеличение вероятности конфигурации с наименьшей энергией, то есть $\beta\geq0$. Показатель экспоненты должен быть безразмерный, для этого определим температуру как $$ \beta=1/k_{B}T $$ Сумма вероятностей всех возможных состояний должна быть равна единице, следовательно, нужен нормировочный фактор $С=1/Z$, который мы положим $$ Z(\beta)=\sum_{E}e^{-\beta E} $$ и назовём статсуммой (Z от нем. Zustandssumme – сумма по состояниям). Обозначим внутреннюю энергию U как среднюю энергию системы: $$ U=\sum_{n}E_{n}\cdot P_{n}(E_{n}) $$ и найдём её полный дифференциал $$ dU=\sum_{n}E_{n}dP_{n}+\sum_{n}P_{n}dE_{n} $$ Свяжем это уравнение с классической термодинамикой. $$ dU=TdS-pdV $$ Это фундаментальное уравнение термодинамики. Согласно второму началу термодинамики $\delta Q=TdS$, где $\delta Q$ - тепло, T - температура, а S - энтропия. $$ TdS=\sum_{n}E_{n}dP_{n} $$ Тепловая энергия $\delta Q$ связана с изменениями вероятностей реализаций конфигураций системы за счёт подведённого тепла. При нагреве вероятность системе оказаться в состоянии с большей энергией возрастает. Вторая часть фундаментального уравнения термодинамики описывает совершенную системой работу: $$ -pdV=\sum_{n}P_{n}dE_{n} $$ и как видно, она связана с изменениями энергий самих конфигураций. Эти изменения энергий ушли на совершение работы.
Для определения энтропии по Больцману возьмём логарифм от $P(E_n)=1/Z\cdot e^{-\beta E_n}$ и выразим $E_n$: $$ E_{n}=-\frac{1}{\beta}(\ln P_{n}+\ln Z) $$ Затем, используя найденное соотношение $-pdV=\sum_{n}P_{n}dE_{n}$ запишем $$ TdS=\sum_{n}E_{n}dP_{n}=-\frac{1}{\beta}\left[\sum_{n}\ln P_{n}dP_{n}+\ln Z\cdot\sum_{n}dP_{n})\right] $$ так как $\sum_{n}P_{n}\equiv1$, то $\sum_{n}dP_{n}=0$ и член с $\ln Z$ можно убрать $$ TdS=\sum_{n}E_{n}dP_{n}=-\frac{1}{\beta}\sum_{n}\ln P_{n}dP_{n} $$ Заметим, что $$ d\sum_{n}P_{n}\ln P_{n}=\sum_{n}\ln P_{n}dP_{n}+\sum_{n}P_{n}\cdot\frac{dP_{n}}{P_{n}}=\sum_{n}\ln P_{n}dP_{n} $$ что приводит к $$ TdS=\sum_{n}E_{n}dP_{n}=-\frac{1}{\beta}d\left(\sum_{n}P_{n}\ln P_{n}\right) $$ откуда окончательно находим формулу Больцмана $$ S=-k_{B}\sum_{n}P_{n}\ln P_{n} $$ где $k_B$ - постоянная Больцмана. Имея на руках выражения для энергии U и энтропии S, можно определить остальные термодинамические величины. Для средней энергии и энтропии после несложных манипуляций с алгеброй находим $$ U=\sum_{n}E_{n}P_{n}=\frac{1}{Z}\sum_{n}E_{n}e^{-\beta E_{n}}=k_{B}T^{2}\cdot\frac{\partial\ln Z}{\partial T} $$ $$ S=\sum_{n}P_{n}\ln P_{n}=\frac{\partial}{\partial T}\left(k_{B}T\ln Z\right) $$ Последние выражения используются при расчёте термодинамических функций из первых принципов, решением уравнения Шредингера для интересующей системы $$ H\Psi_{n}=E_{n}\Psi_{n} $$

О теореме Ирншоу

В предыдущей заметке я привёл оценки размеров молекул и энергий химических взаимодействий. Они по своей природе электростатические. Однако, свести все взаимодействия атомов к притяжению-отталкиванию статических зарядов δ+ и δ- не получится. Теорема Ирншоу утверждает, что всякая неподвижная конфигурация точечных зарядов неустойчива в пространстве.

Доказательство теоремы Ирншоу опирается на теорему Гаусса. Согласно которой поток вектора напряжённости электрического поля $\vec{\Phi}{}_{E}$ через замкнутую поверхность S пропорционален заключённому внутри этой поверхности электрическому заряду Q: $$ \vec{\Phi}{}_{E}=\oint_{S}\vec{E}\cdot d\vec{S}=\frac{Q}{\epsilon_{0}} $$ Другими словами, число силовых линий поля, входящих в поверхность S огораживающую заряд Q, минус число силовых линий выходящих из неё, пропорциональны величине заключенного внутри заряда. Когда внутри ничего нет поток равен нулю. Сколько силовых линий поля зашло, столько и вышло.

Слева – иллюстрация действия теоремы Гаусса. Суммарный поток вектора напряженности (число силовых линий) равен нулю, так как внутри объёма, ограниченного поверхностью S пусто. Справа – противоречие теореме Гаусса, вызванное допущением, будто бы существует устойчивое равновесие для системы точечных зарядов.

Теорема Ирншоу доказывается от противного. Допустим, что какая-то система неподвижных точечных зарядов находится в устойчивом равновесии. Под устойчивостью следует понимать не то, что внешние силы, действующие на заряд равны нулю, а то, что любое малое отклонение вызывает силы, возвращающие заряды на исходные положения.

В этой системе зарядов силы скомпенсированы, но равновесие не устойчиво.

Возьмем любой заряд q этой системы, находящийся в равновесии в точке A, схема на рисунке справа. Если заряд q сместится в близкую точку A', то должна возникнуть сила, направленная к исходной точке A, стремящаяся вернуть заряд на своё место. Пусть $\vec{E}$ – электрическое поле, создаваемое всеми остальными зарядами, кроме q. В точке A' оно должно быть направлено к A, каково бы не было направление смещения AA'. Поле $\vec{E}$ действует на заряд q возвращая его в равновесное положение. Заряд q ограничен замкнутой поверхностью S, так, чтобы все прочие заряды не попадали внутрь. Так как поле $\vec{E}$ на поверхности S всегда направлено к A, то поток вектора $\vec{E}$ через поверхность S ненулевой. Это противоречит теореме Гаусса. Поток должен быть равен нулю, так как это поле создаётся всеми зарядами, расположенными вне S. Противоречие. Такое поле не может существовать, значит, не создаётся устойчивое равновесие.
Однако, как мы знаем, атомы и молекулы, составленные из практически точечных зарядов устойчивы. Речь идёт о неустойчивости системы точечных зарядов. В модели атома Томсона атом представляет собой «булку с изюмом», точечные электроны находятся в статическом равновесии внутри облака положительного заряда ядра. В результате опытов Резерфорда по рассеянию альфа-частиц обнаружено, что ядра скорее точки, как и электроны, а не облака. На основании этого опыта модель Томсона была отвергнута..
Следовательно, атомные системы представляют собой не закрепленные в пространстве заряды, а движущиеся заряды, то есть ядра и электроны в динамике. Законы систем движущихся зарядов описываются в электродинамике.

Статистические справочные таблицы на калькуляторе

Задача. Есть калькулятор, но нет под рукой статистических таблиц. Например, нужны таблицы критических точек распределения Стьюдента для вычисления доверительного интервала. Взять компьютер с Excel? Не спортивно.

Большая точность не нужна, можно воспользоваться приближенными формулами. Идея приведённых ниже формул состоит в том, что преобразованием аргумента все распределения можно так или иначе свести к нормальному. Аппроксимации должны обеспечивать как вычисление кумулятивной функции распределения, так и расчет обратной к ней функции.

Начнём с нормального распределения. $$\Phi(z)=P=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{z}{\sqrt{2}}\right)\right]$$ $$z=\Phi^{-1}(P)=\sqrt{2}\cdot\mathrm{erf}^{-1}(2P-1)$$
Для него требуется вычислить функцию \(\mathrm{erf}(x)\) и обратную к ней. Я воспользовался приближением [1]: $$ \mathrm{erf}(x)=\mathrm{sign}(x)\cdot\sqrt{1-\exp\left(-x^{2}\cdot\frac{\frac{4}{\pi}+ax^{2}}{1+ax^{2}}\right)} $$ $$ \mathrm{erf}^{-1}(x)=\mathrm{sign}(x)\cdot\sqrt{-t_{2}+\sqrt{t_{2}^{2}-\frac{1}{a}\cdot\ln t_{1}}} $$ где \(t_1\) и \(t_2\) вспомогательные переменные: $$ t_{1}=1-x^{2},\:t_{2}=\frac{2}{\pi a}+\frac{\ln t_{1}}{2} $$ а константа \(a=0.147\). Ниже дан код на языке Octave.

function y = erfa(x)
  a  = 0.147;
  x2 = x**2; t = x2*(4/pi + a*x2)/(1 + a*x2);
  y  = sign(x)*sqrt(1 - exp(-t));
endfunction
function y = erfinva(x)
  a  = 0.147; 
  t1 = 1 - x**2; t2 = 2/pi/a + log(t1)/2;
  y  = sign(x)*sqrt(-t2 + sqrt(t2**2 - log(t1)/a));
endfunction

function y = normcdfa(x)
  y = 1/2*(1 + erfa(x/sqrt(2)));
endfunction
function y = norminva(x)
  y = sqrt(2)*erfinva(2*x - 1);
endfunction

Теперь, когда есть функции нормального распределения, приведём аргумент и вычислим t-распределение Стьюдента [2]: $$ F_{t}(x,n)=\Phi\left(\sqrt{\frac{1}{t_{1}}\cdot\ln(1+\frac{x^{2}}{n})}\right) $$ $$ t=F_{t}^{-1}(P,n)=\sqrt{n\cdot\exp\left(\Phi^{-1}(P)^{2}\cdot t_{1}\right)-n} $$ где вспомогательная переменная \(t_1\) есть $$ t_{1}=\frac{n-1.5}{(n-1)^{2}} $$

function y = tcdfa(x,n)
  t1 = (n - 1.5)/(n - 1)**2;
 y = normcdfa(sqrt(1/t1*log(1 + x**2/n)));
endfunction
function y = tinva(x,n)
  t1 = (n - 1.5)/(n - 1)**2;
  y  = sqrt(n*exp(t1*norminva(x)**2) - n);
endfunction

Идея приближенного вычисления распределения \(\chi^{2}\) наглядно представлена формулами [3]: $$ \sigma^{2}=\frac{2}{9n},\:\mu=1-\sigma^{2} $$ $$ F_{\chi^{2}}(x,n)=\Phi\left(\frac{\left(\frac{x}{n}\right)^{1/3}-\mu}{\sigma}\right) $$ $$ \chi^{2}=F_{\chi^{2}}^{-1}(P,n)=n\cdot\left(\Phi^{-1}(P)\cdot\sigma+\mu\right)^{3} $$

function y = chi2cdfa(x,n)
  s2 = 2/9/n; mu = 1 - s2;
  y  = normcdfa(((x/n)**(1/3) - mu)/sqrt(s2));
endfunction
function y = chi2inva(x,n)
 s2 = 2/9/n; mu = 1 - s2;
  y = n*(norminva(x)*sqrt(s2) + mu)**3;
endfunction

Распределение Фишера (для \(n/k\geq3\) и \(n\geq3\)) находится в два шага. Сначала аргумент преобразуется к вычислению распределения Фишера через распределение \(\chi^{2}\) [4], а его мы уже знаем, как вычислить. $$ \sigma^{2}=\frac{2}{9n},\:\mu=1-\sigma^{2} $$ $$ \lambda=\frac{2n+k\cdot x/3+(k-2)}{2n+4k\cdot x/3} $$ $$ F_{f}(x;k,n)=\Phi\left(\frac{\left(\lambda\cdot x\right)^{1/3}-\mu}{\sigma}\right) $$ Найдём обратную функцию, решив квадратное уравнение. $$ q=\left(\Phi^{-1}(P)\cdot\sigma+\mu\right)^{3} $$ $$ b=2n+k-2-4/3\cdot kq $$ $$ D=b^{2}+8/3\cdot knq $$ $$ x=F_{f}^{-1}(P;k,n)=\frac{-b+\sqrt{D}}{2k/3} $$

function y = fcdfa(x,k,n)
  mu = 1-2/9/k; s = sqrt(2/9/k);
  lambda = (2*n + k*x/3 + k-2)/(2*n + 4*k*x/3);
  normcdfa(((lambda*x)**(1/3)-mu)/s)
endfunction
function y = finva(x,k,n)
  mu = 1-2/9/k; s = sqrt(2/9/k);
  q = (norminva(x)*s + mu)**3;
  b = 2*n + k-2 -4/3*k*q;
  d = b**2 + 8/3*k*n*q;
  y = (sqrt(d) - b)/(2*k/3);
endfunction

Список литературы

1. Sergei Winitzki. A handy approximation for the error function and its inverse. February 6, 2008.
2.  Gleason J.R. A note on a proposed Student t approximation // Computational statistics & data analysis. – 2000. – Vol. 34. – №. 1. – Pp. 63-66.
3.  Wilson E.B., Hilferty M.M. The distribution of chi-square // Proceedings of the National Academy of Sciences. – 1931. – Vol. 17. – №. 12. – Pp. 684-688.
4.  Li B. and Martin E.B. An approximation to the F-distribution using the chi-square distribution. Computational statistics & data analysis. – 2002. Vol. 40. – №. 1. pp. 21-26.