这篇文章大概是从高中的时候就打算写，不过考虑到长度过长，来北大之后就一直没有时间动笔。 最近刚好考完了一波期末，两学期的物理化学课程也结课，准备趁着摸鱼的这小段时间简单地介绍一下统计力学当中比较浅薄的概念。 (顺便：强力推荐化院高毅勤和刘剑两位老师合授的物理化学英文版课程， 010-01035200-1006172506-1 ，化院两年以来上的最有价值的课，虽然我是个菜逼，依然学到很多)

对于我个人来说，中学时期接触统计力学思想之前，学习经典热力学漫长又缺乏清晰内核的经历是相当痛苦的。 因此我认为从微观的角度出发，尤其是从分子动力和量子力学的角度出发，结合微观与宏观的关系来理解热力学量的成因是一条更为清楚、更加富有内涵的道路。 并且，理解微观与宏观的联系方法 (eg. Fokker-Planck equation)，对于社会统计、微观/宏观经济学、生态学等其它学科的理解也是至关重要的 (as an interesting example, see this and this 。

作为一个初级的统计力学教程兼一篇安利性质的文章，我认为不应该有过多的具体实现细节，并且希望尽量写得轻松一点。 作最简陋的打算，文章主要会把握住微观与宏观的联系这条主线，考虑相应简单地介绍数学背景和物理图像。 在此之外，如果我还有时间完善文章，作最理想的打算，则物理化学内容方面，希望尽量包括化院本科的化数、物化、中级物化课的初级内容和其它学校一些研究生课程的内容； 算法和代码方面，打算尽量讲清思路，并结合我自己的编码经验，偶尔会给出python, julia, MATLAB 或者类 C 代码的伪实现； 数学方面，希望能摆脱初级微积分的技巧性操作，尽量从抽代、实变、泛函、几何之类的角度，建立高观点来看待复变、概率、ODE、PDE等等需要用到的工具中的技巧，试图寻找比较优雅的理解问题的方法。

文章的主体结构我目前计划是这样：

目前暂时先把结构列出，把坑挖好，内容我打算慢慢施工。。 大概很长时间之后应该能完结吧。 只有章节名没有内容的，以及章节内标记[TODO]的部分都是未完成的区域。[/TODO]

有的章节后面会有“练习”和“思考问题”的部分。 “练习”是该章节非常基础的内容，应该完全掌握。 “思考问题”是基于该章节的比较困难的问题，可能需要后面的内容或者更多知识才能解决。 读者应该尽量思考，但是不一定需要解决这些问题。 (对于高中生来说全部解决是几乎不可能的，因为有时候我会放 open problem 。 当然，虽然不能解决问题，单纯地思考这些问题对于理清思路也是非常有帮助的。)

有的内容十分不容易找到中文资料。 我尽量用中文写这篇文章，但为了准确的考虑，相关名词一般用英文，第一次出现时，我可能简单翻译一下。 对于国内介绍比较少的部分，因为相关的中文学术名称还没有统一，可能只能用英文来写。

对于有的独立性较强的内容，如果篇幅很长，我会把它分离成一篇独立的文章。

因为我菜得可以，文章必定有大量谬误，请留言指正或者发邮件告诉我

主要参考：

[1] Mark E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, Oxford University Press, USA, 2010

[2] Daan Frenkel, Berend Smit, Understanding Molecular Simulation: From Algorithms to Applications, Academic Press, 2002

[3] 朗道，理论物理学教程 (中文译版)

[4] Donald A. McQuarrie, Statistical Mechanics, University Science Books, 2000

[5] David Chandler, Introduction to Modern Statistical Mechanics, Oxford University Press, 1987

[6] M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 2nd ed., 2017.

[7] David Tannor, Introduction to Quantum Mechanics: A Time-Dependent Perspective, University Science Books, 2007

数学、物理、化学和计算机科学基本背景

对于这篇文章，考虑到受众主要是高年级中学生，因此假设读者至少已经学过数学分析、高等代数(或者线性代数)等入门课程，对基本的数学分支和工具有一定了解； 物理上至少学过普通的高中物理，对牛顿定律和简单的量子论有一定了解； 化学上至少知道原子，对热力学有简单了解； 计算机科学上至少会简单使用unix系操作系统，对编程有简单概念。

在开始统计力学之前，我打算先把用到的核心知识略微回顾一遍。如果已经比较了解相关背景，可以跳过这一章。

数学复习

抽代复习

一些结构的定义

群作用

一些常用定理等

实变复习

测度论等

泛函复习

空间概念等

复变复习

复平面的拓扑

Cauchy定理

全纯/亚纯函数

留数定理

概率论复习

公理化定义

连续和离散分布

联合分布

矩

大数定律和中心极限定理

统计学复习

ODE复习

物理复习

牛顿力学复习

经典力学(Classical Mechanics)

经典力学是牛顿提出的一套适用于中等尺度世界的近似运动定律，它忽略相对论效应和量子效应。 尽管牛顿提出经典力学主要是解释天体的运动，经典力学在描述一般物体的运动上也是十分精确的。

传统的分子动力学(molecular dynamics)主要是经典力学在分子尺度的应用。 在分子尺度，量子效应不再能够忽略，但是经典力学在某些方面依然可以作为一个良好的近似。 分子动力学有很多典型应用，例如蛋白质/核酸的折叠模拟，材料科学方面表面催化/表面功能化(surface functionalization)和玻璃的机理等问题的研究，纳米技术中自组装现象的研究等等。

牛顿(Newton)运动定律

1687年左右牛顿给出三条运动定律：

1. 没有外力的情况下，物体或者以速度 \(\mathbf v\) 做匀速直线运动，或者静止不动。
2. \(\mathbf F = m\mathbf a\)
3. 力的作用是相互的，\(\mathbf{F_{AB} = – F_{BA}}\)

牛顿框架下的一些概念

牛顿的框架十分简单，不多作介绍，以下我们列出常用的定义

粒子(particle)

粒子是一个抽象概念，在经典力学下可以认为它相当于质点。在经典框架下，粒子有其质量、电荷等属性。

速度

体系中第\(i\)个粒子的速度为

\[\vec v_i(t) \triangleq \frac{d\vec r_i(t)}{dt}\]

对任何量\(q\)，简记它对时间的全导数为

\[\frac{dq}{dt} \triangleq \dot{q}\]

加速度(acceleration)

\[\vec a_i(t) \triangleq \ddot r_i(t)\]

动量(momentum)

\[\vec p_i \triangleq m\vec v_i\]

力(force)

\[\vec F_i \triangleq \frac{d\vec p_i}{dt} = \dot{\vec{p_i}}\]

路径(path)

经典力学中的路径是坐标空间中的一条连续的线，一般用来描述粒子的运动。 在这里，我们不妨认为单粒子的坐标空间是性质非常好的 \(\mathbb R^3\)，从而总是可以定义欧氏距离(Euclidean Distance)

\[|\vec r_A – \vec r_B| = \sqrt{(x_A – x_B)^2 + (y_A – y_B)^2 + (z_A – z_B)^2}\]

一个质点在经典情况下的运动，可以用Cartesian坐标

\[\vec r_i(t) = (x_i(t), y_i(t), z_i(t))\]

来完全描述。

多粒子的路径依然可以仿照以上的定义得到。 当然，目前我们使用的是一个不太清楚的定义，以后会造成一些问题。 不过，现在我们还可以暂时就用这个不清楚的定义来考虑问题，以免在事情一开始的时候就弄得过于复杂。 之后我们再在此基础上拓展。

注意以下的 Lebesgue 积分

\[\mathcal L = \int_L dl\]

并不一定总是存在，其给出某条轨迹的长度。 也就是说，这条线可能是无穷长的，也可能是有限长的。 另外，路径可能对某个变量可导，也可能不可导。 实际上，可导的情形是非常罕见的。

轨迹(trajectory)

轨迹在经典力学下指在给定的运动条件下，由运动方程生成的相空间中的一条线。 其坐标部分是一条路径。

初值(initial conditions)

我们可以确定任意一个时刻为初始时刻，或者叫零时刻。 在该时刻时 N 个粒子的 2N 个独立条件(比如它们的位置和动量)构成一组初值条件。 这一意义和常微分方程中的初值条件是相同的。

功(work)

对于位移 \(d \vec r\)，功定义为

\[dW \triangleq \vec F \cdot d\vec r\]

对其积分，可以得到沿任意路径\(L\)的功

\[W = \int_L \vec F \cdot d\vec l\]

力场(force field)

经典力学中的力场是一个高维向量，由所有粒子在所有位置和速度下的受力组成

\[\vec F(\vec r_1, \vec r_2, …, \vec r_N, \dot{\vec{r_1}}, \dot{\vec{r_2}}, …, \dot{\vec{r_N}}, t) = \begin{pmatrix}
\vec F_1(\vec r_1, \vec r_2, …, \vec r_N, \dot{\vec{r_1}}, \dot{\vec{r_2}}, …, \dot{\vec{r_N}}, t) \\
\vec F_2(\vec r_1, \vec r_2, …, \vec r_N, \dot{\vec{r_1}}, \dot{\vec{r_2}}, …, \dot{\vec{r_N}}, t) \\
… \\
\vec F_N(\vec r_1, \vec r_2, …, \vec r_N, \dot{\vec{r_1}}, \dot{\vec{r_2}}, …, \dot{\vec{r_N}}, t) \\
\end{pmatrix}\]

保守力(conservative forces)

如果某一力场 \(\vec F(\vec r_1, \vec r_2, …, \vec r_N, \dot{\vec{r_1}}, \dot{\vec{r_2}}, …, \dot{\vec{r_N}}, t)\) 是由某个标量函数 \(U(\vec r_1, \vec r_2, …, \vec r_N, \dot{\vec{r_1}}, \dot{\vec{r_2}}, …, \dot{\vec{r_N}}, t)\) 求偏微分生成的，则该力场称为保守力场，\(U(\vec r_1, \vec r_2, …, \vec r_N, \dot{\vec{r_1}}, \dot{\vec{r_2}}, …, \dot{\vec{r_N}}, t)\) 称为势场或势能函数。

常见的保守场一般都是与粒子速度无关的场，即\(\vec F(\vec r_1, \vec r_2, …, \vec r_N, t)\)

\[\vec F_i(\vec r_1, \vec r_2, …, \vec r_N, t) = -\nabla_i U(\vec r_1, \vec r_2, …, \vec r_N, t)\]

在保守力场中，由积分与路径无关的条件很容易证明沿任意路径的运动做功为

\[W = \int_A^B \vec F \cdot d\vec l = U(B) – U(A)\]

即功与路径无关。一个立即的推论是对所有环形路径\(L\)有

\[W_L = \oint_L \vec F \cdot d\vec l \equiv 0\]

练习：

假设\(U\)定义为每对粒子间距离 \(\{ r_{AB}, r_{AC}, … \}\) 的可导函数，是否存在保守力场？是否一定为保守力场？

如果\(U\)还是某些粒子速度的函数呢？

答案： 是，是，是，否

守恒(conservation)

对于任何定义在相空间的量 \(u\)，若对某个运动的轨迹 \(L\) 满足

\[u(A) = u(B), \forall A, B \in L\]

则称 \(u\) 为该运动的守恒量

动能(kinetic energy)

\[\mathcal K = \sum_i \mathcal K_i \triangleq \sum_i \frac{1}{2} m v_i^2\]

势能(potential energy)

势能又叫位能，通常定义为粒子由于其在力场中的位置所具有的能量。 实际上，在保守场中，它就是生成力场的势场 \(U(\vec p, \vec q, t)\)。

能量(energy)

能量是一个十分有用的运动积分，它是一个满足可加性的守恒量。 在经典力学中，它是动能和势能的和。

运动的解法

牛顿的运动方程实际上是一组二阶微分方程，给定粒子间的力学关系 \(\mathcal U(\vec p, \vec q, t)\) 后，所有的运动轨迹就都能确定。 若再给出独立的 2N 个条件，则求解出系统接下来的演化就只是微分方程问题。 例如初值条件或端值条件分别对应微分方程中初值问题和端值问题。

在牛顿框架下，确定一组初值之后，物体的轨迹是完全确定的。 但对于多粒子系统，其产生的微分方程复杂度过高。 对此有两种方法

• 考虑一个非常理想的简单系统。 简化的系统可以被解析地求解，用处非常有限。不过这种方法可以提供物理上的insight，可以分析不同的初值和外部条件下系统的行为。
• 采用一个较小的系统，通常有\(10^2 \sim 10^9\)个粒子， 对给定的初值和边界条件下进行数值求解(即所谓“模拟(simulations)”)。这种办法可以处理较大的体系，但对于每个不同的情况，都必须重新计算。此外，一般无法确定精确的力场，必须引入模型。

经典理论力学复习

相空间(phase space)的概念

N-particle系统的所有经典动力学(classical dynamics)可以由6N个函数描述：(或者说2dN个，d是维数。由于我们是3维，因此d=3)

\[\{\vec r_1(t), \vec r_2(t), …, \vec r_N(t), \vec p_1(t), \vec p_2(t), …, \vec p_N(t)\}\]

定义相空间向量：

\[\vec x(t) = (\vec r_1(t), \vec r_2(t), …, \vec r_N(t), \vec p_1(t), \vec p_2(t), …, \vec p_N(t))\]

则经典运动(classical motion)都可以用相空间中的一条曲线描述，参数为t。 如前所述，这就是所谓的轨迹或轨线(trajectory)。

为了对轨迹有更直观的感觉，举几个例子。 考虑一些简单的一维单粒子运动：

free particle (自由粒子)

一维自由粒子是一个最简单的散射模型，即一个沿特定方向运动的不受任何力的粒子。 按牛顿力学，它将保持匀速运动。 例如，考虑一个初速度为1，质量为1的粒子，其轨迹见 Figure [fig:phase_space_vis_1d_free_particle]

[fig:phase_space_vis_1d_free_particle]

harmonic oscillator (谐振子)

谐振子是一个最简单的束缚模型，粒子受力与其偏离平衡位置的距离成正比，这被称为 Hooke’s law (虎克定律)

\[F = -kx\]

\[m\ddot x = -kx\]

解微分方程，得

\[x(t) = x(0) \cos{\omega t} + \frac{p(0)}{m\omega} \sin{\omega t}\]

\[p(t) = p(0) \cos{\omega t} – m\omega x(0) \sin{\omega t}\]

这正是椭圆的参数方程。两式平方后线性组合，可以化为非参数的形式

\[\frac{(p(t))^2}{2m} + \frac{1}{2} m \omega^2 (x(t))^2 = C\]

\[C = \frac{(p(0))^2}{2m} + \frac{1}{2} m \omega^2 (x(0))^2\]

我们可以画出 C 不同时，相空间中的不同轨线，见 Figure [fig:phase_space_vis_1d_harmonic_oscillator]。 这些椭圆的短轴长度为 \(\sqrt{2mC}\)，长轴长度为 \(\sqrt{\frac{2C}{m\omega^2}}\)

[fig:phase_space_vis_1d_harmonic_oscillator]

separatrix

Separatrix 是动力系统中的概念，它实际上可以看成是势垒的模型。 严格地来讲，separatrix 就是微分方程中分隔两种行为模式的界面。 在这里，我们不妨考虑最简单的势垒，一个位于原点的高斯波包

\[U(x) = e^{-x^2}\]

从而力为

\[F(x) = -\frac{\partial U(x)}{\partial x} = 2 e^{-x^2} x\]

运动方程为

\[m\ddot x = F(x) = 2 e^{-x^2} x\]

设质量为1，

\[\ddot x – 2 e^{-x^2} x = 0\]

解上面的微分方程，可以简单地得到轨迹。 显然，当从负无穷向正无穷散射时，\(\dot x(0) = \sqrt 2\) 对应的曲线是 separatrix。 初速度若低于之，则将发生反射，初速度若高于之，则将发生透射。

我们画出初速度为\(\dot x(0) = 1, \sqrt 2 – \delta, \sqrt 2 + \delta, 2\)的四条轨迹，见 Figure [fig:phase_space_vis_1d_gaussian_separatrix]。

[fig:phase_space_vis_1d_gaussian_separatrix]

练习：

另一个非常简单的 separatrix 是平面上的摆，均匀重力场下单摆的运动方程为

\[\ddot \theta + {g\over \ell} \sin \theta = 0\]

\(g, \ell\)为常数，其运动自由度同样为1。 找出 separatrix，并在相空间中将其画出。 然后求解 separatrix 两侧的轨线，并各画出1条。 (提示：当使用角度 \(\theta\) 为坐标时，其对应的动量为角动量 \(m\dot \theta\)。）

答案： 参考丁同仁《常微分方程教程》或任何ODE教材。 轨线类似眼睛的形状。

相空间的可视化

相空间通常是一个非常高维的空间，要在3维空间之内画出其中的轨迹，只对1维情况下的单粒子运动是可能的(p和x两个维度，平面)。 在多粒子的情况下，几种有用的可视化方法是

[TODO]

只画出感兴趣的coordinates和对应的动量，其它维度按某个多元函数投射到单个轴上。

将某些感兴趣的坐标画在一个轴上，其对应的动量画在另一个轴上。 这是可视化化学反应的相空间的一种好办法。 具体方法见 De Leon et al., 1991

庞加莱截面(Poincare section)，（补充介绍动力系统背景)

[/TODO]

拉格朗日(Lagrangian)力学

拉格朗日力学是一套和牛顿力学等价的表述方法，也叫拉格朗日表述(Lagrangian formulations)。 一般来说拉格朗日力学要求保守力场。

首先，一个系统的动能为:

\[K(\dot{\vec {r_1}}, \dot{\vec {r_2}}, …, \dot{\vec {r_N}}) = \sum_{i=1}^N \frac{1}{2}m_i\dot{\vec {r_i}}^2\]

拉氏量(Lagrangian，或拉格朗日量)定义为动能减去势能

\[\mathcal L = K(\dot{\vec {r_1}}, \dot{\vec {r_2}}, …, \dot{\vec {r_N}}) – U({\vec {r_1}}, {\vec {r_2}}, …, {\vec {r_N}})\]

欧拉-拉格朗日方程(Eular-Lagrange equation)的形式：

\[\frac{d}{dt}(\frac{\partial\mathcal L}{\partial\dot{\vec {r_i}}})-\frac{\partial\mathcal L}{\partial{\vec {r_i}}} = 0\]

或：

\[\frac{d}{dt}(\frac{\partial\mathcal L}{\partial\dot q_i}) -\frac{\partial\mathcal L}{\partial q_i} = 0\]

将\(\mathcal L\)的形式代入，容易验证以上定义与牛顿的公式等价。 因此，用 Eular-Lagrange equation可得到系统的运动方程。

使用这套系统的好处是可以用一种系统的方法得到广义坐标下的运动方程。

广义坐标(generalized coordinates)

广义坐标一般是一套与Cartesian坐标自由度相同的坐标，因此它和Cartisian坐标一样都能完全确定系统的位置状态：

\[(q_1, q_2, …, q_{3N}) = f(\vec r_1, \vec r_2, …, \vec r_N)\]

其中

\[q_{\alpha} = f_{\alpha}(\vec r_1, \vec r_2, …, \vec r_N), \alpha = 1 \sim 3N\]

一般假设\(f\)是一一映射，即假设\(f\)存在唯一的逆\(g\)，满足

\[\vec r_{i} = g_{i}(q_1, q_2, …, q_{3N}), i = 1 \sim N\]

由链式求导法则，

\[\dot{\vec{r_i}} = \sum_{\alpha = 1}^{3N} \frac{\partial \vec r_i}{\partial q_{\alpha}} \dot{q_{\alpha}}\]

从而在广义坐标下动能可以表述为

\[\begin{aligned}
\tilde K(q_1, q_2, …, q_{3N}, \dot q_1, \dot q_2, …, \dot q_{3N}) & = \sum_{i=1}^N \frac{1}{2}m\dot{\vec {r_i}}^2 \\
& =\sum_{i=1}^N\frac{1}{2}m (\sum_{\alpha = 1}^{3N} \frac{\partial \vec r_i}{\partial q_{\alpha}} \dot{q_{\alpha}})^2 \\
& = \frac{1}{2}\sum_{\alpha = 1}^{3N}\sum_{\beta = 1}^{3N}(\sum_{i=1}^N m_i\frac{\partial \vec r_i}{\partial q_{\alpha}} \cdot\frac{\partial \vec r_i}{\partial q_{\beta}})\dot{q_{\alpha}}\dot{q_{\beta}} \\
& = \frac{1}{2}\sum_{\alpha = 1}^{3N}\sum_{\beta = 1}^{3N}G_{\alpha\beta}(q_1, q_2, …, q_{3N}) \dot{q_{\alpha}}\dot{q_{\beta}} \\
\end{aligned}\]

\(G_{\alpha\beta}(q_1, q_2, …, q_{3N})=\sum_{i=1}^N m_i \frac{\partial \vec r_i}{\partial q_{\alpha}}\cdot \frac{\partial \vec r_i}{\partial q_{\beta}}\)是一个二维张量\(\mathbf G(q_1, q_2, …, q_{3N})\)的矩阵元，该张量称为质量度量张量(mass metric tensor)，度量张量在国内物理学中一般也翻译为度规张量。 显然它是对称矩阵。 将广义坐标按编号排为列向量，则动能可以写成乘法的形式

\[\tilde K(\vec q, \vec{\dot q}) = \frac{1}{2} \vec{\dot q}^\dagger \mathbf G(q_1, q_2, …, q_{3N}) \vec{\dot q} = \frac{1}{2} \vec{\dot q}^\dagger \mathbf G \vec{\dot q}\]

从而有

\[\mathcal L = \tilde K(\vec q, \vec{\dot q}) – \tilde U(\vec q)\]

为了熟悉拉格朗日力学，我们考虑几个简单的例子。

中心势(central potential)

中心势即\(U(r) = U(|r|)\)，势能只与粒子到某个中心的距离有关。

常见的库伦势(Coulomb potential)， 重力势(Gravitational potential)， 谐振子(harmonic potential, 其三维形式为球谐(spherical harmonic)) 等均属于二次中心势。 常用的中心势还有4次势， 刚性球体势(刚球势)， 球中粒子(particle in a sphere), 球壳中粒子(particle in a spherical shell)等。 我们可以首先写出在笛卡尔坐标下某个粒子在中心势中的运动方程：

\[U(x, y, z) = f(\sqrt{x^2 + y^2 + z^2})\]

\[\begin{aligned}
\vec F_i(x_i, y_i, z_i) & = \begin{pmatrix}
-\frac{\partial}{\partial x_i}U(x_i, y_i, z_i) \\
-\frac{\partial}{\partial y_i}U(x_i, y_i, z_i) \\
-\frac{\partial}{\partial z_i}U(x_i, y_i, z_i) \\
\end{pmatrix} \\
& = \begin{pmatrix}
-\frac{\partial}{\partial x_i}f(\sqrt{x_i^2 + y_i^2 + z_i^2}) \\
-\frac{\partial}{\partial y_i}f(\sqrt{x_i^2 + y_i^2 + z_i^2}) \\
-\frac{\partial}{\partial z_i}f(\sqrt{x_i^2 + y_i^2 + z_i^2}) \\
\end{pmatrix} \\
& = \begin{pmatrix}
-f'(\sqrt{x_i^2 + y_i^2 + z_i^2})\frac{\partial}{\partial x_i}\sqrt{x_i^2 + y_i^2 + z_i^2} \\
-f'(\sqrt{x_i^2 + y_i^2 + z_i^2})\frac{\partial}{\partial y_i}\sqrt{x_i^2 + y_i^2 + z_i^2} \\
-f'(\sqrt{x_i^2 + y_i^2 + z_i^2})\frac{\partial}{\partial z_i}\sqrt{x_i^2 + y_i^2 + z_i^2} \\
\end{pmatrix} \\
& = \begin{pmatrix}
-\frac{x_if'(\sqrt{x_i^2 + y_i^2 + z_i^2})}{\sqrt{x_i^2+y_i^2+z_i^2}} \\
-\frac{y_if'(\sqrt{x_i^2 + y_i^2 + z_i^2})}{\sqrt{x_i^2+y_i^2+z_i^2}} \\
-\frac{z_if'(\sqrt{x_i^2 + y_i^2 + z_i^2})}{\sqrt{x_i^2+y_i^2+z_i^2}} \\
\end{pmatrix}
\end{aligned}\]

从而运动方程为

\[\begin{pmatrix}
-\frac{x_if'(\sqrt{x_i^2 + y_i^2 + z_i^2})}{\sqrt{x_i^2+y_i^2+z_i^2}} \\
-\frac{y_if'(\sqrt{x_i^2 + y_i^2 + z_i^2})}{\sqrt{x_i^2+y_i^2+z_i^2}} \\
-\frac{z_if'(\sqrt{x_i^2 + y_i^2 + z_i^2})}{\sqrt{x_i^2+y_i^2+z_i^2}} \\
\end{pmatrix} = \begin{pmatrix}
m\ddot x_i \\
m\ddot y_i \\
m\ddot z_i \\
\end{pmatrix}\]

这是联立的三变量二阶微分方程。 观察到在上式中，\(\sqrt{x_i^2 + y_i^2 + z_i^2}\) 多次出现，因此，采用球极座标是更加合适的坐标系统。 球极变换为

\[q_1 = r = \sqrt{x_i^2 + y_i^2 + z_i^2}\]

\[q_2 = \theta = \arctan{\sqrt{x_i^2 + y_i^2} \over z}\]

\[q_3 = \phi = \arctan{y_i \over x_i}\]

其逆变换为

\[x = r \sin\theta \cos\phi = q_1 \sin{q_2} \cos{q_3}\]

\[y = r \sin\theta \sin\phi = q_1 \sin{q_2} \sin{q_3}\]

\[z = r \cos\theta = q_1 \cos{q_2}\]

求偏导，得

\[\frac{\partial\vec r}{\partial r} = \begin{pmatrix}
\sin\theta \cos\phi \\
\sin\theta \sin\phi \\
\cos\theta \\
\end{pmatrix}\]

\[\frac{\partial\vec r}{\partial\theta} = \begin{pmatrix}
r\cos\theta \cos\phi \\
r\cos\theta \sin\phi \\
-r\sin\theta \\
\end{pmatrix}\]

\[\frac{\partial\vec r}{\partial\phi} = \begin{pmatrix}
-r\sin\theta \sin\phi \\
r\sin\theta \cos\phi \\
0 \\
\end{pmatrix}\]

回忆质量度规张量的定义，得

\[G_{\alpha\beta}(q_1, q_2, q_3) = m_i \frac{\partial \vec r_i}{\partial q_{\alpha}} \cdot \frac{\partial \vec r_i}{\partial q_{\beta}}\]

\[\begin{aligned}
\mathbf G(r, \theta, \phi) & = m_i \begin{pmatrix}
1 & 0 & 0 \\
0 & r^2 & 0 \\
0 & 0 & r^2\sin^2\theta \\
\end{pmatrix}
\end{aligned}\]

从而可以给出动能在广义坐标 \((r, \theta, \phi)\) 下的形式

\[\tilde{K} = \frac{1}{2} m_i (\dot r^2 + r^2 \dot \theta^2 + r^2\sin^2\theta \dot \phi^2)\]

势能的形式非常简单

\[\tilde U = U(r)\]

从而拉氏量在 \((r, \theta, \phi)\) 下为

\[\mathcal L = \frac{1}{2} m_i (\dot r^2 + r^2 \dot \theta^2 + r^2\sin^2\theta \dot \phi^2) – U(r)\]

由欧拉-拉格朗日方程，我们可以马上得到运动方程

\[\begin{gathered}
\frac{d}{dt}(\frac{\partial(\frac{1}{2} m_i (\dot r^2 + r^2 \dot \theta^2 + r^2\sin^2\theta \dot \phi^2) – U(r))}{\partial\dot q_i}) – \\
\frac{\partial(\frac{1}{2} m_i (\dot r^2 + r^2 \dot \theta^2 + r^2\sin^2\theta \dot \phi^2) – U(r))}{\partial q_i} = 0\end{gathered}\]

\[\left\{
\begin{array}{rl}
\frac{d}{dt}(m_i \dot r) – m_i ( r \dot \theta^2 + r \sin^2\theta \dot \phi^2) + U'(r) = 0 \\
\frac{d}{dt}(\dot{\theta } r^2 m_i) – r^2 \dot{\phi }^2 \sin (\theta ) \cos (\theta ) m_i = 0 \\
\frac{d}{dt}(r^2 \dot{\phi } \sin ^2(\theta ) m_i) =0
\end{array}\right.\]

化简，得到以下的运动方程

\[\ddot r – r \dot\theta^2 – r\sin^2\theta \dot\phi^2 + U'(r)/m = 0\]

\[\ddot\theta r^2 + 2 \dot\theta r\dot r – r^2 \dot{\phi }^2 \sin (\theta ) \cos (\theta ) = 0\]

\[r \dot\phi \left(\dot r \sin (2 \theta ) \dot \phi +r \left(\dot \theta \cos (2 \theta ) \dot\phi + \sin (2 \theta ) \ddot \phi \right)\right) = 0\]

这就是对任意方向建系所得的中心势中的运动方程。

在上面的运算中，我们注意到 \(\frac{\partial\mathcal L}{\partial\phi} \equiv 0\)，这是循环坐标(cyclic coordinate, or ignorable coordinate)的一个例子。 更具体地说，如果拉氏量对某个坐标的偏导为0，那么这个坐标就是一个循环坐标。

如果我们能找到某些守恒量，则我们可以通过转动座标系来进一步化简运动方程。 在中心势场中，轨道角动量(orbital angular momentum)是一个比较显然的守恒量，其定义为

\[\vec l = \vec r \times \vec p\]

我们转动座标系，使\(\theta=0\)方向，即原先的z轴方向指向\(\vec l\)的方向，则运动仅仅在 xy 平面内发生，有\(\theta=\pi/2, \dot\theta =0\)。 从而

\[\ddot r – r \dot\phi^2 = -U'(r)/m\] \[r^2\ddot\phi+2r\dot r \dot\phi = 0\]

注意到此时有

\[\frac{d}{dt}\left(r^2 \dot{\phi }\right) = 0\]

此即开普勒第二定律(Kepler’s second law)，或者叫“面速度守恒”(conservation of areal velocity)

二粒子系统

两个相互作用的粒子在不存在外势的情况下称为二粒子系统。 显然，我们可以立即写出其拉氏量

\[\mathcal L = \frac{1}{2}m_1 {\dot{\vec{r_1}}}^2 + \frac{1}{2}m_2 {\dot{\vec{r_2}}}^2 – U(|\vec r_1 – \vec r_2|)\]

一套更加自然的坐标是质心(center-of-mass)

\[\vec R = \frac{m_1 \vec r_1 + m_2 \vec r_2}{m_1 + m_2} = \frac{m_1}{m_1 + m_2} \vec r_1 + \frac{m_2}{m_1 + m_2} \vec r_2\]

和相对坐标(relative coordinates)

\[\vec r = \vec r_1 – \vec r_2\]

注意到这仅仅只是笛卡尔坐标的线性变换。 可以用矩阵的形式写出

\[\begin{pmatrix}
\vec R \\
\vec r
\end{pmatrix} = \begin{pmatrix}
\frac{m_1}{m_1 + m_2} & \frac{m_2}{m_1 + m_2} \\
1 & -1
\end{pmatrix} \begin{pmatrix}
\vec r_1 \\
\vec r_2
\end{pmatrix}\]

逆变换为

\[\begin{aligned}
& \begin{pmatrix}
\frac{m_1}{m_1 + m_2} & \frac{m_2}{m_1 + m_2} \\
1 & -1
\end{pmatrix} ^{-1} \\
& = \begin{pmatrix}
-\frac{1}{-\frac{m_1}{m_1+m_2}-\frac{m_2}{m_1+m_2}} & -\frac{m_2}{\left(m_1+m_2\right) \left(-\frac{m_1}{m_1+m_2}-\frac{m_2}{m_1+m_2}\right)} \\
-\frac{1}{-\frac{m_1}{m_1+m_2}-\frac{m_2}{m_1+m_2}} & \frac{m_1}{\left(m_1+m_2\right) \left(-\frac{m_1}{m_1+m_2}-\frac{m_2}{m_1+m_2}\right)} \\
\end{pmatrix} \\
& = \begin{pmatrix}
1 & \frac{m_2}{m_1+m_2} \\
1 & -\frac{m_1}{m_1+m_2} \\
\end{pmatrix}
\end{aligned}\]

因此

\[\begin{pmatrix}
\vec r_1 \\
\vec r_2
\end{pmatrix} = \begin{pmatrix}
1 & \frac{m_2}{m_1+m_2} \\
1 & -\frac{m_1}{m_1+m_2} \\
\end{pmatrix} \begin{pmatrix}
\vec R \\
\vec r
\end{pmatrix}\]

i.e.

\[\vec r_1 = \vec R + \frac{m_2}{m_1+m_2} \vec r\]

\[\vec r_2 = \vec R – \frac{m_1}{m_1+m_2} \vec r\]

代入拉氏量的表达式，可以马上得到

\[\begin{aligned}
\mathcal L & = \frac{1}{2}m_1 ({\frac{d}{dt}{(\vec R + \frac{m_2}{m_1+m_2} \vec r})})^2 + \frac{1}{2}m_2 ({\frac{d}{dt}{(\vec R – \frac{m_1}{m_1+m_2} \vec r})})^2 – U(|\vec r|) \\
& = \frac{1}{2}(m_1 + m_2) {\dot{\vec R}}^2 + \frac{1}{2}\frac{m_1 m_2}{m_1 + m_2} {\dot{\vec r}}^2 – U(|\vec r|)
\end{aligned}\]

引入总质量\(M = m_1 + m_2\)和约化质量(reduced mass)\(\mu = \frac{m_1 m_2}{m_1 + m_2}\)，化为

\[\mathcal L = \frac{1}{2} M {\dot{\vec R}}^2 + \frac{1}{2}\mu {\dot{\vec r}}^2 – U(|\vec r|)\]

同样代入欧拉-拉格朗日方程，注意到 \(\frac{\partial\mathcal L}{\partial \vec R}\)，因此 \(\vec R\) 是循环坐标，只需考虑 \(\vec r\)即可，运动方程为

\[\mu \ddot{\vec r} = -U'(|\vec r|)\frac{\vec r}{|\vec r|}\]

从上面的两个例子，我们可以看出拉格朗日力学是一种相当优雅的在坐标间变换的方法，根据对称性选择合适的坐标，能够大幅简化我们的问题。

拉格朗日力学的意义还不止如此。 首先我们注意到，在以上的介绍中，我们不加说明地直接给出了欧拉-拉格朗日方程的形式。 事实上，这一方程与一个更加普遍的原理有深刻的联系，即作用量极值原理(action extremization priciple)或作用量稳定值原理。 对经典作用量(classical action)的变分求极值或稳定值，可以得到欧拉-拉格朗日方程，即给出经典路径。 之后我们继续考虑量子力学或量子场论的路径积分表述(path integral formulation)时，还要再回到这里。

勒让德变换(Legendre transforms)

经典理论力学中，另一个等价表述是哈密顿(Hamiltonian)力学。为了阐述Hamiltonian力学，需要先介绍勒让德变换。

如果我们已知一个一元函数 \(f(x)\)，并且假设这个函数具有相当好的性质，比方说在任意的点上都可以对它求导数，即有 \(s=f'(x)\) ，如何找到一种一一映射，将 \(f(x)\) 变换为 \(s\) 的函数 \(\tilde f(s)\) ？这是勒让德变换考虑的一个问题。

考虑 \(f(x)\) 在每一点上的切线的截距 \(b(x)\) ，设\(s=f'(x)=g(x)\) ，则有 \(x=g^{-1}(s)\)，则 \(f(x)=f'(x)\cdot x + b(x)\)，因此，

\[b(x) = b(x(s)) = f(x(s)) – s \cdot x(s)\]

含有与 \(f(x)\) 相同的信息。由此，定义函数

\[\tilde f(s)  = b(x(s)) = f(x(s)) – sx(s)\]

称为 \(f(x)\) 的Legendre变换。

考虑其在多元函数上的推广，则得到广义Legendre变换 (generalized Legendre transform):

\[s_i = \frac{\partial f}{\partial x_i} = g_i(x_1, x_2, …, x_n)\]

\[\tilde f(s_1, s_2, …, s_n) = f(x_1, x_2, …, x_n) – \sum_{i=1}^n s_ix_i\]

注意：Legendre变换也可以仅对一部分变量进行变换，即可以用在变量的一个子集上。

[TODO] 实际上从更简单的角度看，勒让德变换就是底流形切丛到余切丛的微分同胚。写完微分流形的复习之后可以考虑把这部分内容整合过去。 [/TODO]

哈密顿(Hamiltonian)力学

现在我们来考虑如何得到拉格朗日力学另一个等价的表述。 上面我们介绍了对某个多元函数的勒让德变换，那么它有什么作用呢？ 回忆经典力学中动量的定义，有

\[\vec p_i = m_i \dot{\vec{r_i}}\]

拉氏量在笛卡尔系中为

\[\mathcal L = \sum_{i=1}^N \frac{1}{2}m_i{\dot{\vec{r_i}}}^2 – U(\vec r_1, \vec r_2, …, \vec r_N)\]

它是所有粒子位置和速度的多元函数。 可以注意到，对拉氏量求某个速度的偏导，可以得到

\[\frac{\partial\mathcal L}{\partial \dot{\vec{r_j}}} = \frac{\partial(\sum_{i=1}^N \frac{1}{2}m_i{\dot{\vec{r_i}}}^2 – U(\vec r_1, \vec r_2, …, \vec r_N))}{\partial \dot{\vec{r_j}}} = m_j\dot{\vec{r_j}} = \vec p_j\]

它正是该速度对应的动量。 因此，我们对 Lagrangian 做 Legendre transform，应该能获得另一个量，它是所有粒子位置和动量的多元函数

\[\begin{aligned}
& \tilde{\mathcal L}(\vec r_1, \vec r_2, …, \vec r_N, \vec p_1, \vec p_2, …, \vec p_N) \\
& = \mathcal L(\vec r_1, \vec r_2, …, \vec r_N, \dot{\vec{r_1}}(\vec p_1), \dot{\vec{r_2}}(\vec p_2), …, \dot{\vec{r_N}}(\vec p_N)) – \sum_{i=1}^N \vec p_i \cdot \dot{\vec{r_i}}(\vec p_i) \\
& = \sum_{i=1}^N \frac{1}{2}m_i{\dot{\vec{r_i}}}^2 – U(\vec r_1, \vec r_2, …, \vec r_N) – \sum_{i=1}^N \vec p_i \cdot \frac{\vec p_i}{m_i} \\
& = – \sum_{i=1}^N \frac{{\vec p_i}^2}{2m_i} – U(\vec r_1, \vec r_2, …, \vec r_N)
\end{aligned}\]

这个量的负数，即对拉格朗日量的速度做勒让德变换所得量的相反数 \(-\tilde{\mathcal L}\) ，被称为哈密顿量，用\(\mathcal H\)表示。 我们注意到，在经典力学中，哈密顿量正是系统的总能量。

\[\mathcal H(\vec r_1, \vec r_2, …, \vec r_N, \vec p_1, \vec p_2, …, \vec p_N) = \sum_{i=1}^N \frac{{\vec p_i}^2}{2m_i} + U(\vec r_1, \vec r_2, …, \vec r_N)\]

我们现在来证明当采用广义坐标时，仍然能通过上述的变换给出哈密顿量。

在广义坐标下，每个广义坐标对应的动量，即其共轭动量为

\[p_\alpha \triangleq \frac{\partial\mathcal L}{\partial\dot q_\alpha}\]

代入拉氏量的表达式，我们得到

\[p_\alpha = \sum_{\beta=1}^{3N} G_{\alpha\beta}(q_1, q_2, …, q_{3N}) \dot q_\beta\]

或者，写成乘法的形式

\[\vec p = \mathbf G(q_1, q_2, …, q_{3N}) \vec{\dot q}\]

逆变换为

\[\vec{\dot q} = \mathbf G^{-1}(q_1, q_2, …, q_{3N})\vec p\]

显然有

\[\begin{aligned}
\vec p ^\dagger \vec{\dot q} & = \vec{\dot q}^\dagger \vec p \\
2 \tilde K = \vec p ^\dagger \mathbf G^{-1} \vec p & = \vec{\dot q}^\dagger \mathbf G \vec{\dot q}
\end{aligned}\]

拉氏量在广义坐标下我们之前已经给出，写成乘法的形式为

\[\begin{aligned}
\mathcal L & = \frac{1}{2}\sum_{\alpha = 1}^{3N}\sum_{\beta = 1}^{3N}G_{\alpha\beta}(q_1, q_2, …, q_{3N}) \dot{q_{\alpha}}\dot{q_{\beta}} – \tilde U \\
& = \frac{1}{2} \vec{\dot q}^\dagger \mathbf G \vec{\dot q} – \tilde U
\end{aligned}
\label{eqn:general_coordinates_Lagrangian}\]

勒让德变换给出

\[\begin{aligned}
\tilde{\mathcal L} & = \mathcal L – \sum_{\alpha=1}^{3N} p_\alpha \dot q_\alpha \\
& = \frac{1}{2} \vec{\dot q}^\dagger \mathbf G \vec{\dot q} – \tilde U – \vec p ^\dagger \mathbf G^{-1}\vec p \\
& = \tilde K – \tilde U – 2 \tilde K \\
& = -\tilde K – \tilde U \\
& = -\tilde{\mathcal H}
\end{aligned}\]

因此广义动量下哈密顿量形式为

\[\begin{aligned}
\tilde{\mathcal H} & = \tilde K + \tilde U \\
& = \frac{1}{2} \vec{\dot q}^\dagger \mathbf G \vec{\dot q} + \tilde U \\
& = \frac{1}{2} \vec p ^\dagger \mathbf G^{-1} \vec p+ \tilde U \\
& = \frac{1}{2} \vec p ^\dagger \vec{\dot q}+ \tilde U = \frac{1}{2} \vec{\dot q}^\dagger \vec p+ \tilde U
\end{aligned}\]

由于广义坐标在接下来比笛卡尔坐标更为常用，我们接下来不再标出广义坐标下各种量的\(\sim\)号

与Lagrangian mechanics下由Eular-Lagrange equations生成运动方程类似，Hamiltonian mechanics下生成运动方程的方程是

\[\dot q_\alpha = \frac{\partial\mathcal H}{\partial p_\alpha}\]

\[\dot p_\alpha = -\frac{\partial\mathcal H}{\partial q_\alpha}\]

这两个方程被称为哈密顿运动方程(Hamilton’s equations of motion)或者哈密顿正则方程(the canonical equations of Hamilton)，后一种叫法是因为哈密顿运动方程给出的相流是一族正则变换，即时间演化是保辛结构的。我们稍后会再继续讨论这件事。

如前所述，哈密顿运动方程与拉格朗日力学、保守场下的牛顿力学是完全等价的。 要证明这一点，只要将\(\mathcal H = -\mathcal L + 2K\)代入，验证 Eular-Lagrange equations成立即可。

如果求Hamiltonian对时间的全导数，则可以发现

\[\begin{aligned}
\frac{\dd \mathcal H}{\dd t} & = \sum_{\alpha=1}^{3N} \left( \frac{\partial\mathcal H}{\partial p_\alpha}\dot p_\alpha + \frac{\partial\mathcal H}{\partial q_\alpha}\dot q_\alpha \right) \\
& = \sum_{\alpha=1}^{3N} \left( -\frac{\partial\mathcal H}{\partial p_\alpha}\frac{\partial\mathcal H}{\partial q_\alpha} + \frac{\partial\mathcal H}{\partial q_\alpha}\frac{\partial\mathcal H}{\partial p_\alpha} \right) \\
& \equiv 0
\end{aligned}\]

这说明哈密顿量是运动的守恒量。 这其实就只是能量守恒定律。

如果考虑任意相空间上不显含时的函数\(a(\vec x(t))\)随时间的演化，显然有

\[\begin{aligned}
\frac{da}{dt} & = \frac{\partial a}{\partial \vec x} \cdot \frac{d \vec x}{d t} \\
& = \sum_{\alpha=1}^{3N} \left( \frac{\partial a}{\partial p_\alpha} \dot p_\alpha + \frac{\partial a}{\partial q_\alpha} \dot q_\alpha \right) \\
& = \sum_{\alpha=1}^{3N} \left( – \frac{\partial a}{\partial p_\alpha} \frac{\partial\mathcal H}{\partial q_\alpha} + \frac{\partial a}{\partial q_\alpha} \frac{\partial\mathcal H}{\partial p_\alpha} \right) \\
\end{aligned}
\label{eqn:time_evolution_of_function_Poisson_bracket}\]

引入一个新的记号，称为泊松括号(Poisson bracket)

\[\{ a, b \} = \sum_{\alpha=1}^{3N} \left( – \frac{\partial a}{\partial p_\alpha} \frac{\partial b}{\partial q_\alpha} + \frac{\partial a}{\partial q_\alpha} \frac{\partial b}{\partial p_\alpha} \right)\]

则以上的关系可以写作

\[\dot a(t) = \{a, \mathcal H\}\]

显然，相空间上不显含时的函数 \(a(\vec x(t))\) 为守恒量的 充分必要条件是 \(\{a, \mathcal H\} \equiv 0\)。

作为一个例子，考虑动量守恒的条件。由

\[\{\mathbf P, \mathcal H\} = 0\]

知

\[\sum_i \{\vec p_i, \mathcal H\} = \{\mathbf P, \mathcal H\} = 0 = \sum_i \dot {\vec {p_i}} = \sum_i \vec F_i\]

这实际上就是系统所受外力之和为0。 (由牛顿第三定律，系统中内力之和为0)

一个立即的推论是，由正则方程，有

\[0 = \sum_i \dot {\vec {p_i}} = – \sum_i \frac{\partial \mathcal H}{\partial \vec{r_i}}\]

如果定义平移变换

\[\vec{r_i}’ = \vec r_i + \vec a, i = 1, 2, …, N\]

显然有

\[- \sum_i \frac{\partial \mathcal H}{\partial \vec{r_i}} = 0 \Leftrightarrow – \sum_i \frac{\partial \mathcal H}{\partial \vec{r_i}’} = – \sum_i \frac{\partial \mathcal H}{\partial \vec{r_i}’}\frac{\partial \vec{r_i}’}{\partial \vec{r_i}} = – \sum_i \frac{\partial \mathcal H}{\partial \vec{r_i}} = 0\]

也就是说，在平移变换下系统保持动量守恒。 以上的平移变换正是所谓平移群(translation group)的群作用。 当系统动量守恒时，在任意的平移变换下，哈密顿量保持不变。

类似于平移群与动量守恒，如果Hamiltonian在群\(\mathcal G\)群作用下守恒，则必定存在一个对应的守恒定律(conservation law)。 这被称为诺特定理(Noether’s theorem)。这是一个十分深刻的结果，在经典力学和量子力学中将用到很多次。

相空间不可压缩性(phase space incompressibility)

之前我们介绍了相空间和相空间矢量的概念。 如果我们对相空间矢量求时间导数，可以发现

\[\frac{d\vec x}{dt} = \begin{pmatrix}
\dot p_1 \\
\dot p_2 \\
… \\
\dot p_{3N} \\
\dot q_1 \\
\dot q_2 \\
… \\
\dot q_{3N}
\end{pmatrix} = \begin{pmatrix}
– \frac{\partial \mathcal H}{q_1} \\
– \frac{\partial \mathcal H}{q_2} \\
… \\
– \frac{\partial \mathcal H}{q_{3N}} \\
\frac{\partial \mathcal H}{p_1} \\
\frac{\partial \mathcal H}{p_2} \\
… \\
\frac{\partial \mathcal H}{p_{3N}} \\
\end{pmatrix} = \eta(\vec x)
\label{eqn:phase_space_velocity}\]

可以看出广义相空间速度 \(\dot{\vec x}\) 是仅依赖于 \(\vec x\) 的函数。 这与流体力学中的流场十分相似。 在流体力学中，不可压缩流体具有一种十分简单的性质：由于其不可压缩，因此流场中应该既没有源(sources)，也没有漏(sinks)。 这在数学上意味着处处散度为0

\[\nabla_{\vec x} \cdot \dot{\vec x} \equiv 0
\label{eqn:phase_space_incompressibility}\]

上式被称为相空间不可压缩性。利用正则方程很容易证明上式

\[\begin{aligned}
\nabla_{\vec x} \cdot \dot{\vec x} & = \sum_\alpha \left( \frac{\partial}{\partial p_\alpha} \dot p_\alpha + \frac{\partial}{\partial q_\alpha} \dot q_\alpha \right) \\
& = \sum_\alpha \left( – \frac{\partial}{\partial p_\alpha} \frac{\partial \mathcal H}{\partial q_\alpha} + \frac{\partial}{\partial q_\alpha} \frac{\partial \mathcal H}{\partial p_\alpha} \right) \\
& = 0
\end{aligned}\]

注意到我们以上的证明用到了正则方程，这意味着相空间不可压缩性仅仅在 Hamiltonian 力学或其等价力学中成立。 事实上，对于一般的系统，例如存在内部摩擦损耗的系统，相空间不可压缩性通常是不成立的。

练习：

现将一个理想一维谐振子放入某种油中，假设其所受阻力与速度成正比，求解其轨迹，并证明相空间不可压缩性在处处均不成立。

答案： 该体系被称为阻尼谐振子。见 Feynman Lectures on Physics 的某一卷(小时候看的，忘了是哪卷了，见谅。。应该是第一卷吧)

哈密顿运动方程的辛结构(symplectic structure)

实际上，我们可以显式地写出上一段中 \(\eta(\vec x)\) 的形式。

\[\dot{\vec x} = \mathbf M \frac{\partial \mathcal H}{\partial \vec x}\]

\[\mathbf M = \begin{pmatrix}
\mathbf 0 & \mathbf I \\
\mathbf{-I} & \mathbf 0 \\
\end{pmatrix}\]

一个动力系统如果能用上式表达，则称其具有辛结构(symplectic structure)。 考虑该系统从0时刻演化到t时刻，由于哈密顿方程对于不同初值的解唯一， \(\vec x_0\) 到 \(\vec x_t\) 是一一映射，即 \(\vec x_t = \vec x_t(\vec x_0)\)。 因此，考虑所有\(\vec x_0\) 到 \(\vec x_t\)的映射，可以将其看作相空间上的一个坐标变换。 该变换的Jacobian matrix (雅可比矩阵)是

\[\mathbf J = \frac{\partial \vec x_t}{\partial \vec x_0}\]

注意在这里我们使用了需要展开两次的 notation，具体地说，就是其矩阵元是

\[J_{kl} = \frac{\partial x_t^k}{\partial x_0^l}\]

所谓的symplectic property (辛属性)，是指满足

\[\mathbf{M=J^TMJ}\]

举例来说，之前我们计算过一维谐振子的轨迹

\[x(t) = x(0) \cos{\omega t} + \frac{p(0)}{m\omega} \sin{\omega t}\]

\[p(t) = p(0) \cos{\omega t} – m\omega x(0) \sin{\omega t}\]

为了看到这一系统具有辛属性，不妨来求其 Jacobian matrix

\[\mathbf J = \begin{pmatrix}
\frac{\partial x(t)}{\partial x(0)} & \frac{\partial x(t)}{\partial p(0)} \\
\frac{\partial p(t)}{\partial x(0)} & \frac{\partial p(t)}{\partial p(0)}
\end{pmatrix} = \begin{pmatrix}
\cos{\omega t} & \frac{\sin{\omega t}}{m\omega} \\
-m\omega \sin{\omega t} & \cos{\omega t}
\end{pmatrix}\]

容易验证

\[\begin{aligned}
\mathbf{J^TMJ} & = \begin{pmatrix}
\cos{\omega t} & -m\omega \sin{\omega t} \\
\frac{\sin{\omega t}}{m\omega} & \cos{\omega t}
\end{pmatrix} \begin{pmatrix}
0 & 1 \\
-1 & 0 \\
\end{pmatrix} \begin{pmatrix}
\cos{\omega t} & \frac{\sin{\omega t}}{m\omega} \\
-m\omega \sin{\omega t} & \cos{\omega t}
\end{pmatrix} \\
& = \begin{pmatrix}
m\omega \sin{\omega t} & \cos{\omega t} \\
-\cos{\omega t} & \frac{\sin{\omega t}}{m\omega}
\end{pmatrix} \begin{pmatrix}
\cos{\omega t} & \frac{\sin{\omega t}}{m\omega} \\
-m\omega \sin{\omega t} & \cos{\omega t}
\end{pmatrix} \\
& = \begin{pmatrix}
0 & 1 \\
-1 & 0 \\
\end{pmatrix} \\
& = \mathbf M
\end{aligned}\]

简正模分析

以上我们简单复习了一些理论力学的基本概念，现在不妨尝试一些初步的应用。 我们仍然从最简单的体系开始做起。 考虑这样一个聚合物长链最简单的模型：将所有原子或者基团视为没有体积的一些小球，相邻的小球之间有一根弹簧连接，如 Figure [fig:normal_mode_analysis_3d_harmonic_polymer_model] 所示。

[fig:normal_mode_analysis_3d_harmonic_polymer_model]

假设这些弹簧或者说键的平衡长度是 \(\{b_i\}\)，为了现在简单起见，我们再假设这些小球的质量也都相同，则该体系的 Hamiltonian 是

\[\mathcal H = \sum_{i=1}^{N} \frac{{\vec p_i}^2}{2m} + \sum_{i=1}^{N-1} (|\vec r_{i+1} – \vec r_i| – b_i)\]

即使这么简单的体系对于我们现在来说还是看起来略微有些复杂。 因此，我们暂时仅仅考虑一维时的情况，之后可以自然地将结果推广到3维。 由于所有运动被限制在一条直线上，可以将3维坐标矢量替换为1维坐标。 新的 Hamiltonian 是

\[\mathcal H = \sum_{i=1}^{N} \frac{{p_i}^2}{2m} + \frac{1}{2}m\omega^2 \sum_{i=1}^{N-1} (|x_{i+1} – x_i| – b_i)^2\]

做了如此多的限制，现在可以简单地处理这个体系了。 为了进一步简化问题，首先引入一组线性坐标变换

\[\eta_i = x_i – x_{i,0}, \\
x_{i,0} – x_{i+1,0} = b_i\]

从而在新的坐标和它的共轭动量下 Hamiltonian 可以写作

\[\mathcal H = \sum_{i=1}^{N} \frac{{p_{\eta_i}}^2}{2m} + \frac{1}{2}m\omega^2 \sum_{i=1}^{N-1} |\eta_{i+1} – \eta_i|^2\]

由正则方程，马上可以得到运动方程是

\[\dot \eta_i = \frac{p_{\eta_i}}{m}\]

\(i = 2, 3, 4 … N-1\) 时有

\[\dot p_{\eta_i} = -m\omega^2 (2\eta_i – \eta_{i+1} – \eta_{i-1})\]

在端区有

\[\dot p_{\eta_1} = -m\omega^2 (\eta_1 – \eta_2)\]

\[\dot p_{\eta_N} = -m\omega^2 (\eta_N – \eta_{N-1})\]

这一组方程是线性微分方程组，我们将其转化为线性方程组问题，非常容易可以得到它的一个基本解组。 代入

\[\eta_i = \sum_{k=1}^N C_k a_{ik} e^{i \omega_k t}\]

得

\[\begin{aligned}
& \sum_{k=1}^N \omega_k^2 C_k a_{ik} e^{i \omega_k t} \\
& = \omega^2 (2\sum_{k=1}^N C_k a_{ik} e^{i \omega_k t} – \sum_{k=1}^N C_k a_{i+1, k} e^{i \omega_k t} – \sum_{k=1}^N C_k a_{i-1, k} e^{i \omega_k t}) \\
& = \omega^2 \sum_{k=1}^N C_k e^{i \omega_k t} (2a_{ik} – a_{i+1, k} – a_{i-1, k})
\end{aligned}\]

\[\sum_{k=1}^N \omega_k^2 C_k a_{1k} e^{i \omega_k t} = \omega^2 \sum_{k=1}^N C_k e^{i \omega_k t} (a_1 – a_2)\]

\[\sum_{k=1}^N \omega_k^2 C_k a_{Nk} e^{i \omega_k t} = \omega^2 \sum_{k=1}^N C_k e^{i \omega_k t} (a_N – a_{N-1})\]

写成矩阵方程

\[\omega_k^2 \vec a_k = \omega^2 \mathbf A \vec a_k, k=1, 2, …, N\]

\[\vec a_k = (a_{1k}, a_{2k}, …, a_{Nk})^\dagger\]

\(\mathbf A\) 是如下的三对角矩阵

\[\mathbf A = \begin{pmatrix}
1 & -1 & 0 & … & 0 & 0 \\
-1 & 2 & -1 & … & 0 & 0 \\
0 & -1 & 2 & … & 0 & 0 \\
… & … & … & … & … & … \\
0 & 0 & 0 & … & 2 & -1 \\
0 & 0 & 0 & … & -1 & 1 \\
\end{pmatrix}\]

i.e.

\[\mathbf{aO} = \mathbf{Aa}\]

\[\mathbf a = (\vec a_1, \vec a_2, … , \vec a_N)\]

\[\mathbf O = diag((\omega_1^2, \omega_2^2, …, \omega_N^2)./\omega^2)\]

对角化\(\mathbf A\)，立即可以得到

\[\frac{\omega_k^2}{\omega^2} = 2 – 2\cos{\frac{(k-1)\pi}{N}}\]

所有待定系数都已经求出，到这里，我们已经解出了运动方程。 那么，这个例子和所谓的简正模分析有什么关系呢？ 我们注意到，利用对角化\(\mathbf A\)的酉矩阵\(\mathbf a\)再次引入坐标变换，

\[\vec \zeta = \mathbf a \vec \eta\]

注意到这只是简单的线性变换，计算其质量度规张量非常简单。 将计算结果代入我们之前得到的 Lagrangian 在广义坐标下的形式([eqn:general_coordinates_Lagrangian])中， 再进行 Legendre 变换即可得到体系在新坐标下的 Hamiltonian。

\[\mathcal H = \sum_{k=1}^{N} \frac{p_{\zeta_k}^2}{2m} + \sum_{k=1}^{N} \frac{1}{2}m\omega_k^2 \zeta_k^2 = \sum_{k=1}^{N} \mathcal H_k\]

现在很清楚了： 在新的坐标 \(\vec \zeta\) 下，我们解除了原本 Hamiltonian 中存在的耦合。 这意味着新坐标 \(\vec \zeta\) 的每一个分量都对应一个独立的运动模式，这就是所谓的简正模(normal mode)。 由于转变回原坐标的所有变换均为线性变换，因此，该体系的任何运动均可以用简正模的线性叠加描述。

对于一般的体系，我们同样可以进行类似以上的分析。 一般地，质量加权坐标指将粒子质量与其笛卡尔坐标耦合得到的一套广义坐标。 假设体系一个平衡构型为

\[\begin{aligned}
\vec r_0 & = (\vec r_{1,0}, \vec r_{2,0}, …, \vec r_{N,0})^\dagger \\
& = (x_{1,0}, y_{1,0}, z_{1,0}, x_{2,0}, y_{2,0}, z_{2,0}, …, x_{N,0}, y_{N,0}, z_{N,0})
\end{aligned}\]

那么相对该构型的质量加权坐标可以定义为

\[\vec q = \mathbf M^\frac{1}{2} (\vec r – \vec r_0)\]

其中 \(\mathbf M\) 是质量矩阵，定义为

\[\mathbf M = diag((m_1, m_1, m_1, m_2, m_2, m_2, …, m_N, m_N, m_N))\]

在质量加权坐标下，动能的形式可以写得比较简单。

\[\mathcal K = \frac{1}{2} \dot{\vec q}^\dagger \dot{\vec q}\]

假设平滑、允许泰勒展开且仅依赖位置的势能面，则有 \(k\) 次展开式

\[\begin{aligned}
\mathcal U(\vec r) & = \sum_{|\alpha|\leq k} \frac{D^\alpha \mathcal U(\vec r_0)}{\alpha!} (\vec r – \vec r_0)^\alpha + \sum_{|\alpha|=k} h_\alpha(\vec r)(\vec r – \vec r_0)^\alpha, \\
& \lim_{\vec r \to \vec r_0}h_\alpha(\vec r)=0.
\end{aligned}\]

在平衡位置附近展开，由虚功原理显然所有一阶项为0。 0阶项不影响运动，因此不妨首先考虑单独取出二阶项

\[\frac{1}{2} (\vec r – \vec r_0)^T (\mathbf M^T)^\frac{1}{2} \mathbf H \mathbf M^\frac{1}{2} (\vec r – \vec r_0)
= \frac{1}{2} \vec q^T \mathbf H \vec q
= \frac{1}{2} (\vec r – \vec r_0)^T \mathbf F (\vec r – \vec r_0)\]

其中，\(\mathbf H\) 称为 Hessian 矩阵，\(\mathbf F\)则是质量加权的 Hessian 矩阵

\[H_{ij} = \frac{\partial^2 \mathcal U}{\partial r_i \partial r_j}, \quad
\mathbf F = (\mathbf M^T)^\frac{1}{2} \mathbf H \mathbf M^\frac{1}{2}\]

如果我们对角化 Hessian 矩阵，则得到

\[\mathbf H = \mathbf L^T \boldsymbol \Lambda \mathbf L\]

\[\frac{1}{2} (\vec r – \vec r_0)^T (\mathbf M^T)^\frac{1}{2} \mathbf H \mathbf M^\frac{1}{2} (\vec r – \vec r_0)
= \frac{1}{2} (\vec r – \vec r_0)^T (\mathbf M^T)^\frac{1}{2} \mathbf L^T \boldsymbol \Lambda \mathbf L \mathbf M^\frac{1}{2} (\vec r – \vec r_0)\]

我们得到简正坐标

\[\vec Q = \mathbf L \mathbf M^\frac{1}{2} (\vec r – \vec r_0) = \mathbf L \vec q\]

此时势能面的二阶项表示为

\[\frac{1}{2} \vec Q^T \boldsymbol \Lambda \vec Q\]

此时，如果势能面没有更高阶项，则这里得到的简正坐标显然是一组已经解耦的坐标。

从以上的分析我们可以看出，得到简正坐标有很多种方式。 一种办法是对角化 Hessian ，然后对质量加权坐标进行变换； 另一种办法是注意到

\[\mathbf F = (\mathbf M^T)^\frac{1}{2} \mathbf H \mathbf M^\frac{1}{2} = (\mathbf M^T)^\frac{1}{2} \mathbf L^T \boldsymbol \Lambda \mathbf L \mathbf M^\frac{1}{2}\]

，因此可以对角化质量加权 Hessian 矩阵，然后对笛卡尔坐标差进行变换。

不过显然，由于我们在以上的分析中并未处理展开的二阶项以上部分，任意体系的运动并不见得总是能用有限个模式给出(见思考问题)。

为了更直观地理解简正模分析，我们来看几个例子。

首先我们来看刚刚已经解出的模型。 我们把它套到一个 N=3 的最简单情况看一看，质量不如就用1，立即可以得到简正模的频率和变换矩阵

\[\omega_k = \sqrt{\omega ^2 \left(2-2 \cos \left(\frac{\pi (k-1)}{N}\right)\right)}
= \left\{0,\omega ,\sqrt{3} \omega \right\}(k=1,2,3)\]

\[\mathbf L = \begin{pmatrix}
1 & -2 & 1 \\
-1 & 0 & 1 \\
1 & 1 & 1 \\
\end{pmatrix}, \quad
\mathbf M = \mathbf I\]

为了得到笛卡尔坐标，求变换矩阵的逆

\[\mathbf L^{-1} = \begin{pmatrix}
\frac{1}{6} & -\frac{1}{2} & \frac{1}{3} \\
-\frac{1}{3} & 0 & \frac{1}{3} \\
\frac{1}{6} & \frac{1}{2} & \frac{1}{3} \\
\end{pmatrix}\]

从而

\[\vec r = \mathbf L^{-1} \vec Q\]

\[\begin{pmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{pmatrix} = \begin{pmatrix}
\frac{1}{6} & -\frac{1}{2} & \frac{1}{3} \\
-\frac{1}{3} & 0 & \frac{1}{3} \\
\frac{1}{6} & \frac{1}{2} & \frac{1}{3} \\
\end{pmatrix} \begin{pmatrix}
A_1 \sin{(\sqrt{3} \omega t + \varphi_3)} \\
A_2 \sin{(\omega t + \varphi_2)} \\
A_3 \sin{(0t + \varphi_1)}
\end{pmatrix}\]

我们画出来看看。利用附录 [src:draw_1d_polymer] 中的程序画图，结果见图 [fig:normal_mode_analysis_1d_3_particles]

[fig:normal_mode_analysis_1d_3_particles]

不难看出，第一个简正模频率最大，这种模式被称为 不对称伸缩 (asymmetric stretch)。 第二个模式稍慢一些，叫做 对称伸缩 (symmetric stretch)。 第三个模式是任何体系都具有的一个模式，它的频率为0，即对应整个体系的平动。

完全相同的流程不难应用到更长的链上。 不过，由于高次多项式没有求根公式，想要解析求解对角化长度超过7的链的 Hessian 存在一定的困难。 不过，数值求解显然是非常简单的。 我们不妨试试画出10个粒子时的简正模，见图 [fig:normal_mode_analysis_1d_10_particles]

[fig:normal_mode_analysis_1d_10_particles]

上面的这个例子非常简单，现在我们来考虑稍微复杂一些的问题。

显然，上述例子的一个明显的问题是，它是一维的模型，而我们的现实世界是3维的。 对于3维的模型，由于自由度几乎提高了3倍，即使对于最简单的结构而言，势能面都将远远比一维的情况复杂。 比方说考虑一个水分子，它的势能将与其两个 O–H 键长，它的 H–O–H 键角都有关。 而如果我们考虑一个苯分子，则我们甚至能数出 30 个自由度。

一般地，对于非线型的三维 N 原子分子，其自由度是 \(3N – 6\)，即去掉3个平动自由度和3个转动自由度。 对于线型的分子，由于其存在一个旋转对称性，自由度将变成 \(3N – 5\)。

上面这个例子的另一个明显的问题是，它对粒子间的力学关系过于简化。 它没有考虑结构对于特定角度的偏好，此外用简单的二次势描述原子核之间的作用也并不准确。

作为我们3维情况下简正模分析的第一个例子，我们还是从最简单的入手。 不如就来考虑水分子。 为了对水分子进行简正模分析，我们有一些需要注意的地方。 首先，我们从模型进入现实，必须立即考虑的第一个问题就是单位问题。 在之前处理我们的 toy model 的时候，我们并没有考虑过任何单位的问题，遇到质量、频率等等常量时一般都随便设单位为1。 虽然处理实际问题时有时也能这样做，但是一般都会造成数据的数量级发生严重混乱。为了避免这种情况，在考虑微观问题的时候，通常会用原子单位。而处理宏观问题时，则常用 SI 单位。

水分子的势能面并不能用简单的手段得到。 由于体系比较小，我们一般想要做得精确一些。 因此，获得其势能面的一个好办法是利用量子力学来进行数值计算扫描。 具体的量子力学计算属于量子化学的内容，已经在另一篇博文里介绍，这里不再多说。

势能面算出后，我们可以立即得到水分子的平衡构型。 由于水分子的振动不会到达平面之外，我们转动整个空间，将水分子放置在 \(xz\) 平面上。

水分子的平衡构型为

H       1.8112496379040899      -0.0009820612379843       0.0000000000000000
O       0.0012030107158682       0.0015494457436835       0.0000000000000000
H      -0.4497832071840029       1.7545145430537885       0.0000000000000000

根据计算出的势能面，可以利用数值差分法计算 Hessian，见附录 [src:h2opot1]

得到

       0.5433304962802362      -0.0180087044193336       0.0000000000000000      -0.5282877659644728       0.0287668412771091       0.0000000000000000      -0.0150427303157634      -0.0107581368577755       0.0000000000000000
      -0.0180087044193336       0.0492422472250758       0.0000000000000000      -0.0339686404630287      -0.0447425612431704       0.0000000000000000       0.0519773448823623      -0.0044996859819054       0.0000000000000000
       0.0000000000000000       0.0000000000000000      -0.0000000008140854       0.0000000000000000       0.0000000000000000       0.0000000024783640       0.0000000000000000       0.0000000000000000      -0.0000000016642785
      -0.5282877659644728      -0.0339686404630287       0.0000000000000000       0.6014740736562072      -0.1112039999771998       0.0000000000000000      -0.0731863076917344       0.1451726404402285       0.0000000000000000
       0.0287668412771091      -0.0447425612431704       0.0000000000000000      -0.1112039999771998       0.5445865807590787       0.0000000000000000       0.0824371587000907      -0.4998440195159085       0.0000000000000000
       0.0000000000000000       0.0000000000000000       0.0000000024783640       0.0000000000000000       0.0000000000000000      -0.0000000049567278       0.0000000000000000       0.0000000000000000       0.0000000024783639
      -0.0150427303157634       0.0519773448823623       0.0000000000000000      -0.0731863076917344       0.0824371587000907       0.0000000000000000       0.0882290380074978      -0.1344145035824530       0.0000000000000000
      -0.0107581368577755      -0.0044996859819054       0.0000000000000000       0.1451726404402286      -0.4998440195159085       0.0000000000000000      -0.1344145035824530       0.5043437054978138       0.0000000000000000
       0.0000000000000000       0.0000000000000000      -0.0000000016642785       0.0000000000000000       0.0000000000000000       0.0000000024783639       0.0000000000000000       0.0000000000000000      -0.0000000008140853

水分子的质量非常简单。在原子单位(au)下，电子的质量为1，因此质子的质量是 1835，从而氧原子核的质量为 1835*16。质量矩阵为

\[\mathbf M_{H_2O} = diag((1835, 1835, 1835, 29360, 29360, 29360, 1835, 1835, 1835))\]

对角化 \(\mathbf M_{H_2O}^\frac{1}{2} \mathbf H_{H_2O} \mathbf M_{H_2O}^\frac{1}{2}\)，得到本征频率

\[\omega_{H_2O} = \begin{pmatrix}
0.00032296650318747373 \\
0.00030504212671449024 \\
0.00005645303081529721 \\
-2.2570086454964211\times 10^{-13} \\
-1.5184812538552072\times 10^{-12} \\
4.630303351512593\times 10^{-13} \\
-8.491102033762549\times 10^{-21} \\
2.9840575854547277\times 10^{-21} \\
9.015139585682635\times 10^{-21} \\
\end{pmatrix}\]

后6个本征值原本应为0，由于数值计算误差而残留，应该将其舍弃。 如果我们此时从原子单位换回常用的单位制，则可以得到水分子的3种振动模式的频率

\[\omega_{H_2O} = \begin{pmatrix}
3944.32527787983 \\
3833.30955982890 \\
1649.06439026534
\end{pmatrix} cm^{-1}\]

同样，可以利用变换矩阵画出水分子的模式，见图 [fig:normal_mode_analysis_H2O] 。 (懒得打开VMD了。。先放张别人的图偷个懒。 此图来自这里)

[fig:normal_mode_analysis_H2O]

练习：

积分与求和允许交换顺序的条件是什么？(提示：从 Fubini’s theorem 立即导出) 偏微分呢？

矩阵可对角化的条件是什么？为什么上面的分析中能保证 Hessian 矩阵能对角化？(提示：\(f_{ij} = f_{ji}\))

Hessian 是正定的吗？什么情况下它是正定的？

思考问题：

简正模分析时，可以解耦到有限个模式的条件是什么？

实验中观察到水分子的振动光谱峰位于 \(3755.9, 3657.1, 1594.7 cm^{-1}\)。 假如势能面是完全正确的，以上对水分子的计算准确吗？为什么会出现很大的误差？

谐振子近似不成立，该如何计算各种物理量？

作用量(action)

[TODO]

拉格朗日的作用量积分(action integral)，变分原理(variational principle)，

local functional

由拉格朗日量的变分推Eular-Lagrange方程

在光滑、性质良好的力场下，初值问题保证解的唯一性。end-point变分法得到的轨迹是否也能保证解的唯一性？（不）

其它作用量，哈密顿最小作用量原理(Hamilton’s principle of least action)

参考Goldstein, Classical Mechanics

[/TODO]

拘束系统和高斯最小拘束原理

多数时候，考虑运动系统时，其间的粒子满足一定的约束关系(Constraints)。 这些约束可能是为了简化问题(例如把振动频率很高的化学键看作是定长从而减小自由度)， 也有可能是物理上确实存在的约束(比方说我们把系统装在刚性的盒子里面，那么就会有边界约束，再比方说我们将系统和控温或者控压设施耦合，等等)

N个粒子的系统没有拘束的情况下有3N个自由度。当存在\(N_c\)个约束的时候，总自由度为\(3N – N_c\)。

holonomic constraint(完整约束) 是指仅涉及位置和时间的约束，其约束方程形式为

\[\sigma_k(q_1, q_2, …, q_{3N}, t), \quad k=1, 2, …, N_c\]

比如说，几个粒子之间存在一个刚性约束。 再比如说，粒子被限制在一个球面上运动。 这些都是典型的holonimic constraint。

而nonholonomic constraint(非完整约束) 则还涉及到速度(或者说动量)

\[\zeta (q_1, q_2, …, q_{3N}, \dot q_1, \dot q_2, …, \dot q_{3N}, t)\]

比如说，恒温，即保持系统的总动能不变，就是一个典型的nonholonimic constraint。

一种处理约束的方法是选取合适的广义坐标。 由于自由度的减少，我们总是能找到一组 \(3N – N_c\) 个广义坐标来表述我们的系统，且这 \(3N – N_c\) 个坐标之间是独立的。 确定合适的坐标之后，再用以上介绍的理论力学方法(用质量度规张量做变换，再代入EL方程或者哈密顿方程)就能方便地得到运动方程。 例如，我们之前已经介绍过中心势下合适的广义坐标，如果我们想处理被约束在球面上运动的系统，此时 \(r\) 成为非必要的坐标，因此选取坐标 \((\theta, \phi)\) 即可。

选取广义坐标的方法虽然其存在性可以保证，但对某些体系实际处理时会较为复杂。实际上，处理经典力学问题时最常用的方法是拉格朗日不定乘子法(Lagrange undetermined multipliers)。

[TODO]由变分原理推拉格朗日不定乘子法[/TODO]

注意到即使在含时的 holonomic constraints 下，Hamiltonian 依然是守恒量。

[TODO]

高斯最小拘束原理(Gauss’s principle of least constraint)

高斯运动方程(Gauss’s equation of motion)

\[m\ddot{\mathbf r} = \mathbf F – \frac{\nabla \nabla \sigma \cdot \cdot \dot{\mathbf r} \dot{\mathbf r} + \nabla \sigma \cdot \mathbf F / m}{|\nabla \sigma|^2 / m}\nabla \sigma\]

\(\dot{\mathbf r} \dot{\mathbf r}\) is the velocity–vector dyad, and \( \nabla \nabla \sigma \cdot \cdot \dot{\mathbf r} \dot{\mathbf r}\) indicates a full contraction of the two tensors \( \nabla \nabla \sigma\) and \(\dot{\mathbf r} \dot{\mathbf r}\)

高斯运动方程守恒能量和约束，并且无法从Hamiltonian或者Lagrangian得到，因此，它是一种非哈密顿动力系统(non-hamiltonian dynamical system)，见本章末。

[/TODO]

思考问题：

证明存在性：当系统处于 \(N_c\) 个约束下时，总是存在一组不超过 \(3N – N_c\) 个独立的广义坐标，可以完备地描述系统的运动。

刚体(rigid body)运动

[TODO]

刚体是指？

用处？

角速度angular velocity，转动惯量moment of inertia，角动量angular momentum，力矩torque，

欧拉角(Euler angles)

quaterions & quaternions

[/TODO]

非哈密顿系统

量子力学复习

化学复习

热力学复习

经典热力学(thermodynamics)中的概念

第一定律

第二定律

第三定律

动力学复习

经典统计力学的理论基础

从这一章开始，我们正式进入主线，也就是统计力学。 作为开头，还是从最简单的部分起步。 因此，在这一章首先考虑纯经典情况下的统计力学，即在纯粹牛顿的框架下，我们能得到何种结论。

在纯经典的框架下，得到的结果将会出现一些难以理解的矛盾。 此外，在推导的过程中我们将不得不面临许多概念和表示上的瑕疵。 这实际上向我们揭示了经典力学的不足之处，在引入量子力学之后，才能完全地解决这些矛盾。

系综(ensemble)的概念

考虑一个有多个粒子的系统。 如果说我们知道这个系统中所有粒子的总能量为 \(E\)，我们能否知道其中每个粒子的能量？ 显然，当总能量一定时，事实上有很多种将能量分配到各个粒子上的办法。 这意味着，两个系统的某些宏观属性相同，不能说明他们处在相同的状态下，即他们并不一定处在相空间的相同点上。 事实上，对于一般的多粒子系统，宏观属性相同的系统的数量是不可胜数的。

如果我们把大量满足同样的某组宏观属性约束，并且微观力学关系都相同的系统放在一起，并且它们各自是由不同的初值演化而来，从而在任意同一时刻，所有系统都具有不同的状态，那么这就构成一个所谓的 系综(ensemble)

比如说，如果将很多个初值不同的，体积(V)、粒子数(N)、粒子种类(x)和总能量(E)都相同的，粒子间无相互作用的隔离系统算作一个集合，那么就可以得到一个十分简单的系综。

当然，宏观性质并不止于这些。 同样，如果将很多个与某个热浴、或者某个恒压活塞、或者某个 particle reservoir 相连的系统算作一个集合，那么这些系统就会具有相同的温度(T)、或者压力(p)、或者渗透压(\(\Pi\))，等等宏观性质，从而我们得到各种不同的系综。

上面的这些系综在时间演化的条件下它们满足的宏观属性约束并没有任何变化，也就是说它们是处在 平衡(equilibrium) 中。 它们属于所谓的 平衡系综(equilibrium ensembles) 。 事实上，系综同样也可以由各种被外力或者外场驱动的系统定义，这类系综叫做 非平衡系综(non-equilibrium ensembles)

在经典统计力学中，宏观的可观测量可由对系综中所有成员进行平均得到。 系综中的每个成员事实上都是相空间中的一个向量或者点 \(\vec x(t)\)。 对于任意相空间上的多元函数 \(a(\vec x(t)): \mathds R^{6N} \mapsto \mathds R\)，其对应的可观测量可以形式地定义为

\[A(t) \triangleq \frac{1}{\mathcal Z} \sum_{\lambda = 1}^{\mathcal Z} a(\vec x_{\lambda}(t)) \triangleq \expval{a}\]

当然，在经典情况下，上式的求和也可能是积分。

相空间体积和刘维尔(Liouville)定律

为了继续讨论下面的问题，我们需要先研究相空间中的体积元随演化的变化的问题。

假设我们的系统从点 \(\vec x_0\) 开始演化，到t时刻到达\(\vec x_t\)。 前面我们曾说过，由于哈密顿力学下 \(\vec x_0 \leftrightarrow \vec x_t\) 是一一映射，故这一演化相当于坐标变换 \(\vec x_0 \to \vec x_t\)。 设 \(\mathbf J\) 是该变换的 Jacobian 矩阵，由于这是 \(\mathds R^{6N} \mapsto \mathds R^{6N}\) 的映射，该矩阵显然是方阵，从而我们可以求它的 determinant

\[J(\vec x_t; \vec x_0) = det \left[ \mathbf J \right]\]

如果我们考虑以上面这个点为中心的相空间中的一个体积元

\[d\vec x = dq_1 dq_2 … dq_{3N} dp_1 dp_2 … dq_{3N}\]

，则显然有

\[d\vec x_t = J(\vec x_t; \vec x_0) d\vec x_0\]

(注意不要与我们之前约定的求导记号混淆！矢量对矢量的偏导按我们的记号将展开成二阶张量，而这里的全微分符号则意味着一个体积元！)

因此，要给出 \(d\vec x_0 \to d\vec x_t\) 的关系，需要求\(det \left[ \mathbf J \right]\)。回忆线性代数的两个关系

\[det\left[ e^{\mathbf A} \right] = e^{Tr\left[ \mathbf A \right] }\]

\[d Tr\left[ \mathbf A \right] =  Tr\left[ d \mathbf A \right]\]

故

\[\begin{aligned}
det\left[ \mathbf J \right] & = e^{Tr\left[ ln \mathbf J \right]} \\
\frac{d}{dt} det\left[ \mathbf J \right] & = \frac{d}{dt} e^{Tr\left[ ln \mathbf J \right]} \\
& = e^{ Tr \left[ ln \mathbf J \right] } \frac{d}{dt} Tr \left[ ln \mathbf J \right] \\
& = det \left[ \mathbf J \right] Tr \left[ \frac{d\mathbf J}{dt} \mathbf J^{-1} \right] \\
& = det \left[ \mathbf J \right] \sum_{k, l} \frac{d\mathbf J_{kl}}{dt} \mathbf J^{-1}_{lk} \\
\end{aligned}\]

由于

\[\frac{d \mathbf J_{kl}}{dt} = \frac{d}{dt}\frac{\partial x_{tk}}{\partial x_{0l}} = \frac{\partial \dot x_{tk}}{\partial x_{0l}}\]

\[J^{-1}_{lk} = \frac{\partial x_{0l}}{\partial x_{tk}}\]

故

\[\begin{aligned}
\frac{d}{dt} det[\mathbf J] & = det[\mathbf J] \sum_{k, l} \frac{\partial \dot x_{tk}}{\partial x_{0l}} \frac{\partial x_{0l}}{\partial x_{tk}} \\
& = det[\mathbf J] \sum_{k} \frac{\partial \dot x_{tk}}{\partial x_{tk}} \\
& = det[\mathbf J] \nabla_{\vec x_t} \cdot \dot{\vec {x_t}} \\
& = 0
\end{aligned}\]

最后一个等号是根据我们之前导出的相空间不可压缩性质，即式 [eqn:phase_space_incompressibility] 。 因此，

\[\frac{d}{dt}J(\vec x_t; \vec x_0) = 0\]

又由于 \(J(\vec x_0; \vec x_0) = 1\)，

\[J(\vec x_t; \vec x_0) \equiv 1\]

\[d\vec x_t = d\vec x_0\]

此即刘维尔定理。

在我们以上的推导中，我们隐式地(implicitly)假设了相空间具有平直欧氏空间的度量(metric)。在此度量下，刘维尔定理仅仅在Hamiltonian dynamics下成立。在以后我们可以看到，对于non-Hamiltonian的系统，采用黎曼空间(Riemann space)是更加自然的选择。在此几何下，相空间的体积元与度规张量的行列式有关，并且存在广义的刘维尔定理 (Tuckermann et al., 2001)。我们将在稍后的章节继续讨论这一问题，并将其应用在非哈密顿体系的分子动力学上。

系综分布函数(distribution function)和刘维尔方程

前面我们说到，系综中含有许多系统，由于初值不同，在任意同一时刻t，各个系统在相空间中各自处于不同点。 如果我们定义 \(f(\vec x, t)\) 为系综中所有成员在相空间的分布密度，则其满足正定，归一。 因此，\(f(\vec x, t)\) 是一个概率密度 \(f(\vec x, t)\)，称为系综分布函数(ensemble distribution function)或者相空间分布函数(phase space distribution function)。

按照上面的定义，显然在某个体积元 \(d\vec x\) 中成员的数量为 \(f(\vec x, t) d\vec x\) 。 那么，在某个封闭区域 \(\Omega\) 内，成员的数量就是

\[M_\Omega = \int_\Omega f(\vec x, t) d\vec x\]

我们考虑某个轨迹上的系综分布函数 \(f(\vec x_t, t)\) ，对于确定的区域 \(\Omega\)，显然其内的成员数量随时间演化将发生变化

\[M_\Omega(t) = \int_\Omega f(\vec x_t, t) d\vec x_t\]

\[\frac{d}{dt}M_\Omega(t) = \frac{d}{dt}\int_\Omega f(\vec x_t, t) d\vec x_t\]

从另一个角度来考虑，区域 \(\Omega\) 内成员数量的变化显然也应该与分布函数的变化有关

\[\frac{d}{dt}M_\Omega(t) = \int_\Omega \frac{\partial}{\partial t}f(\vec x_t, t) d\vec x_t\]

注意由于我们不想让 \(\Omega\) 在积分时”变形”，我们不应该让 \(\vec x_t\) 演化。 因此，考虑分布函数的变化时应当只对其显式含时的部分求导，即上式中是求偏导而非全导。 (见思考问题)

我们继续考虑这个成员数量变化问题。 由于演化时，系综的成员不会凭空消失或者生成(现在可以先假设如此)，对于相空间的任意单连通光滑区域 \(\Omega\)，应在围绕该区域的超平面 \(S\) 上满足通量方程。

\[- \frac{d}{dt}M_\Omega(t) = \int_S dS \dot{\vec x}_t \cdot \hat n f(\vec x_t, t)\]

由散度定理，

\[\int_S dS \dot{\vec x}_t \cdot \hat n f(\vec x_t, t) = \int_\Omega d\vec x_t \nabla_{\vec x_t} \cdot (\dot{\vec x}_t f(\vec x_t, t))\]

\(\hat n\) 是 \(dS\) 的外向单位法向量。 联立以上几式，可得

\[\int_\Omega \frac{\partial}{\partial t}f(\vec x_t, t) d\vec x_t = – \int_\Omega d\vec x_t \nabla_{\vec x_t} \cdot (\dot{\vec x}_t f(\vec x_t, t))\]

id est

\[\int_\Omega \left( \frac{\partial}{\partial t}f(\vec x_t, t) + \nabla_{\vec x_t} \cdot (\dot{\vec x}_t f(\vec x_t, t)) \right) d\vec x_t = 0\]

由于区域\(\Omega\)的任意性，以上方程必须处处成立，即

\[\frac{\partial}{\partial t}f(\vec x_t, t) + \nabla_{\vec x_t} \cdot (\dot{\vec x}_t f(\vec x_t, t)) = 0\]

用乘法求导公式，第二项可以拆分为

\[\nabla_{\vec x_t} \cdot (\dot{\vec x}_t f(\vec x_t, t)) = \dot{\vec x}_t \cdot \nabla_{\vec x_t} f(\vec x_t, t) + f(\vec x_t, t) \nabla_{\vec x_t} \cdot \dot{\vec x}_t\]

由于相空间不可压缩性，上式第二项为0。从而

\[\frac{\partial}{\partial t}f(\vec x_t, t) + \dot{\vec x}_t \cdot \nabla_{\vec x_t} f(\vec x_t, t) = 0\]

注意到以上等式事实上就是

\[\frac{d}{dt}f(\vec x_t, t) = 0\]

此即经典力学的刘维尔方程。 我们可以看出，刘维尔方程意味着哈密顿流在相空间上保 Lebesgue 测度。

刘维尔方程也可以写成积分形式

\[f(\vec x_t, t) = f(\vec x_0, 0)\]

与我们之前得到的刘维尔定理相结合，得到

\[f(\vec x_t, t)d\vec x_t = f(\vec x_0, 0)d\vec x_0\]

其物理意义是相当明显的：一个相空间中的体积元随时间演化，其中的系综成员数不会发生变化。 这一点对于我们来说非常重要，我们稍后会看到，正是依赖于这一性质，我们才可以保证能够在任意时间点对系综求平均，从而之前对于可观测量的定义才是自洽的。

回忆之前对相空间结构的讨论。 我们曾给出，相空间矢量的广义速度仅是其自身的函数，即式 [eqn:phase_space_velocity]。 因此，我们得到刘维尔方程的另一个写法

\[\frac{\partial}{\partial t}f(\vec x_t, t) + \eta(\vec x_t, t) \cdot \nabla_{\vec x_t} f(\vec x_t, t) = 0\]

\(\eta(\vec x, t)\) 事实上是我们的速度场。 利用速度场写出以上方程的好处是，我们现在取消了对 \(t\) 的显式依赖，从而可以在上式中移除这个 dummy label (见思考问题)，得到

\[\frac{\partial}{\partial t}f(\vec x, t) + \eta(\vec x, t) \cdot \nabla_{\vec x} f(\vec x, t) = 0
\label{eqn:Liouville_equation_fixed_point}\]

进一步地，由式 [eqn:time_evolution_of_function_Poisson_bracket]，得到

\[\frac{\partial}{\partial t}f(\vec x, t) + \left\{ f(\vec x, t), \mathcal H(\vec x, t) \right\} = 0\]

这是刘维尔方程的另一个常见形式，常常被简写成

\[\frac{\partial}{\partial t}f + \left\{ f, \mathcal H \right\} = 0\]

思考问题：

对于一般的函数 \(f(\vec x, t): \mathds R^{N+1} \mapsto \mathds R\) 和 N 维实矢量 \(\vec x\) ，式

\[\frac{d}{dt}\int_\Omega f(\vec x, t) d\vec x = \int_\Omega \frac{\partial}{\partial t}f(\vec x, t) d\vec x\]

何时成立？

为什么式 [eqn:Liouville_equation_fixed_point] 中能移除 \(\vec x_t\) 对 \(t\) 的依赖？ 一般地，可以移除 dummy label 的条件是什么？

刘维尔方程的平衡解

之前我们给出了宏观可观测量的定义。 如果我们用系综分布函数写下来，应该是

\[A(t) = \int a(\vec x) f(\vec x, t) d\vec x\]

对于处在热力学平衡中的系统，我们知道其可观测量不应该随时间发生变化。 因此

\[\frac{\partial}{\partial t}f(\vec x, t) \overset{a.e.}{=\joinrel=} 0\]

由刘维尔方程，我们知道这意味着

\[\left\{ f(\vec x), \mathcal H(\vec x) \right\} \overset{a.e.}{=\joinrel=} 0\]

上式的解可以记为

\[f(\vec x) \propto \mathcal F (\mathcal H(\vec x))\]

归一化引入因子 \(\mathcal Z\)

\[\mathcal Z = \int \mathcal F (\mathcal H(\vec x)) d\vec x\]

则

\[f(\vec x) = \frac{1}{\mathcal Z} \mathcal F (\mathcal H(\vec x))\]

\(\mathcal Z\) 被称为 配分函数(partition function) 。 配分函数是统计力学中的中心概念，它是一个特定系综内所能访问到的微观态的数量的一个度量，每一个系综都有一个特定的依赖于定义这个系综的宏观物理量的配分函数，因此，各种热力学量都可以从各种配分函数或者它的偏微分等等中得到。

我们用配分函数可以写出平衡态的可观测量的计算公式

\[\begin{aligned}
A & = \int a(\vec x) f(\vec x) d\vec x \\
& = \frac{1}{\mathcal Z} \int a(\vec x) \mathcal F (\mathcal H(\vec x)) d\vec x \\
\end{aligned}\]

微正则系综和分子动力学简介

上一章我们介绍了系综的概念和一些基本的定理，我们从这一章开始讨论一些更具体的内容。

微正则系综的配分函数

现在来考虑一个最简单的系综。 为了构成这个系综，首先每个系统里要有一些粒子。 既然是最简单的系综，那么显然这些粒子应该都是相同的。 如果粒子会消失或者生成，那情况就会有点复杂，所以我们假设每个系统里粒子的数量是常数 \(N\) 。 这些粒子要运动，就要有一定的空间，因此不妨认为它们装在一个容积为 \(V\) 的容器中。 最后，如果我们想要考虑最简单的情况，那么不允许这些粒子与外界交换能量是一个不错的选择，这意味这这些系统内部的能量必定是守恒的，假设这个总能量值为 \(E\)。

这些条件放在一起，就足以构成一个最简单的系综了，这个系综就是所谓的 微正则系综(microcanonical ensemble) 。 由于微正则系综的 \(N, V, E\) 都守恒，因此也叫 NVE系综。

考虑微正则系综的系综分布函数。 从理论力学我们知道，对于哈密顿系统，其总能量正是其哈密顿量。 有

\[\mathcal H(\vec x) = E\]

我们需要选择一个合理的 \(\mathcal F\)，使得 \(\mathcal F(\mathcal H(\vec x))\) 满足以上的限制。 因此，分布函数可以写成与 Dirac \(\delta\) 相乘的形式

\[f(\vec x) = \mathcal N(\vec x) \delta(\mathcal H(\vec x) – E)\]

这样，所有位于超平面 \(\delta(\mathcal H(\vec x) – E)\) 外的点分布函数均为0，能够满足能量的约束。 那么，现在我们还有一个部分没有确定：即对于位于超平面 \(\delta(\mathcal H(\vec x) – E)\) 上的点，其分布函数如何？ 换句话说，\(\mathcal N(\vec x)\) 的形式是什么？

事实上，我们目前并没有更多的信息。 因此，这是无法决定的。 为了继续我们的讨论，我们需要引入一个假设，即所谓的 等先验概率假设(postulate of equal a priori probability)。

等先验概率假设是统计力学的一个基本假设。 我们假设，对于一个处于平衡中的隔离系统，在所有其能访问到的微观态中，它处在任何一个态下的概率都是完全相等的。 因此，如果能访问到的微观态总数为 \(\Omega(N, V, E)\) ，那么处于每个态的概率都是 \(\frac{1}{\Omega(N, V, E)}\) 。

在等概率假设下，我们可以给出分布函数的具体形式了：在等能量超平面上均匀分布，在其它位置是0，即 \(\mathcal N(\vec x) \equiv \mathcal N\) 为一归一化常数。

\[f(\vec x) = \mathcal N \delta(\mathcal H(\vec x) – E)\]

由于微观态的均匀分布，微观态总数 \(\Omega\) 显然应该与等能量超平面的面积成正比。 写成积分的形式，这个 \(6N\) 重积分是

\[\begin{aligned}
\Omega(N, V, E) & = M \int d\vec x \delta(\mathcal H(\vec x) – E) \\
& = M \int d\vec p \int_{D(V)} d\vec r \delta(\mathcal H(\vec x) – E)
\end{aligned}\]

注意我们在第二个等号中将位形空间的积分与动量空间的积分分开了。 这是由于在微正则系综中，系统是一个具有有限体积的体系，相空间是 \(\mathds R^{3N} \otimes V^{N}\)。 举例来说，如果我们的体系是放在一个长、宽、高为 \(a,b,c\) 的盒子中，那么于一顶点建系，相空间就是 \(\mathds R^{3N} \otimes [0,a]^{N} \otimes [0,b]^{N} \otimes [0,c]^{N}\) 。 上面的积分此时就应该写成

\[\begin{gathered}
\Omega(N, V, E) \\
= M \int d\vec p \int_0^a dx_1 \int_0^b dy_1 \int_0^c dz_1 \dotsi \int_0^a dx_N \int_0^b dy_N \int_0^c dz_N \delta(\mathcal H(\vec x) – E)\end{gathered}\]

上式中还有一个常数 \(M\) 没有定义。 为了获得一个 \(M\) 的取值，我们不妨借用量子力学中的不确定性关系。 由于 \(\Delta x \Delta p \sim h\) ，要数出一个区域内有多少微观态，只要数其中有多少个超立方体 \(\Delta \vec x\) 即可。 由于 \(\Delta \vec x = (\Delta x)^{3N} (\Delta p)^{3N} = h^{3N}\) ，因此，在一个厚度为 \(E_0\) 的薄壳内，每单位面积中的超立方体数量就是 \(E_0 / h^{3N}\) 。 因此，我们得到

\[M = \frac{E_0}{h^{3N}}\]

然而，以上的结果是错误的。 我们前面说到，我们所考虑的 N 个粒子是 完全相同的，这意味着如果我们交换其中的两个粒子，体系不应该发生任何变化。 既然交换前后体系没有任何区别，那么显然这两者是同一个微观态。 因此，我们在上面数微观态的数量时，多数了 \(N!\) 倍。

不幸的是，经典力学并没有办法处理这个问题。 实际上，当我们一个个写下体系中每个粒子的位置和动量矢量的时候，我们已经隐式地对它们进行了编号。 这些编号实际上相当于给每个粒子戴上了一个不同的标签，因此，它们事实上变成了可区分的粒子。

我们现在暂时放下这个问题，因为它需要加入量子力学才能完全地解决。 我们将多数的情况扣除，因此，现在得到

\[M = \frac{E_0}{N! h^{3N}}\]

结合之前对可观测量的定义，那么

\[\begin{aligned}
A = \expval{a} & = \int a(\vec x) f (\vec x) d\vec x \\
& = \frac{\int a(\vec x) \delta(\mathcal H(x) – E) d\vec x}{\int \delta(\mathcal H(x) – E) d\vec x} \\
& = \frac{M}{\Omega(N,V,E)} \int a(\vec x) \delta(\mathcal H(x) – E) d\vec x \\
\end{aligned}\]

玻尔兹曼关系

Boltzmann 和 Max Planck 给出了以下的关系

\[S(N,V,E) = k_B \ln \Omega(N,V,E)\]

\(k_B\) 称为 Boltzmann 常数，这个公式也经常被称为 Boltzmann 熵公式。 我们现在来看如何得到它。

回忆经典热力学里的一些定律。 由热力学第一定律，

\[dE = dW + dQ
= -P dV + \mu dN + T dS\]

移项可以得到

\[dS = \frac{1}{T} dE + \frac{1}{T} P dV – \frac{1}{T} \mu dN\]

同时，对于 \(S(N,V,E)\) 显然还有

\[dS = \left(\frac{\partial S}{\partial E}\right)_{N,V} dE + \left(\frac{\partial S}{\partial V}\right)_{N,E} dV + \left(\frac{\partial S}{\partial N}\right)_{V,E} dN\]

与上面比较，可得

\[\left(\frac{\partial S}{\partial E}\right)_{N,V} = \frac{1}{T}\]

\[\left(\frac{\partial S}{\partial V}\right)_{N,E} = \frac{P}{T}\]

\[\left(\frac{\partial S}{\partial N}\right)_{V,E} = -\frac{\mu}{T}\]

由可观测量的定义，

\[\frac{1}{T} = \frac{M}{\Omega(N,V,E)} \int a_{1/T}(\vec x) \delta(\mathcal H(x) – E) d\vec x = \left(\frac{\partial S}{\partial E}\right)_{N,V}
\label{eqn:observable_invT}\]

假设 \(\Omega\) 与 \(S\) 的关系为

\[S(N,V,E) = C G(\Omega(N,V,E))\]

，其中 \(C\) 为一实系数， \(G\) 为一实函数，则有

\[\left(\frac{\partial S}{\partial E}\right)_{N,V}
= C G'(\Omega(N,V,E)) \left(\frac{\partial \Omega(N,V,E)}{\partial E}\right)_{N,V}\]

由于

\[\begin{aligned}
\left(\frac{\partial \Omega(N,V,E)}{\partial E}\right)_{N,V}
& = \left(\frac{\partial (M \int d\vec x \delta(\mathcal H(\vec x) – E))}{\partial E}\right)_{N,V} \\
& = M \int \left(\frac{\partial \delta(\mathcal H(\vec x) – E)}{\partial E}\right)_{N,V} d\vec x \\
& = M \int \left(\frac{\partial \ln \delta(\mathcal H(\vec x) – E)}{\partial E}\right)_{N,V}{\delta(\mathcal H(\vec x) – E)} d\vec x
\end{aligned}
\label{eqn:Boltzmann_relation_derivation}\]

代入上式，得到

\[\left(\frac{\partial S}{\partial E}\right)_{N,V} = C G'(\Omega(N,V,E)) M \int \left(\frac{\partial \ln \delta(\mathcal H(\vec x) – E)}{\partial E}\right)_{N,V}{\delta(\mathcal H(\vec x) – E)} d\vec x\]

与式 [eqn:observable_invT] 比较，有

\[k \left(\frac{\partial \ln \delta(\mathcal H(\vec x) – E)}{\partial E}\right)_{N,V} = a_{1/T}(\vec x)\]

\[\frac{C}{k} G'(\Omega(N,V,E)) = \frac{1}{\Omega(N,V,E)}\]

将上式积分，不难看出未知函数 \(G\) 正是对数函数。 由于参数 \(C, k\) 都是任意选择的，积分，得

\[S(N,V,E) = C G(\Omega(N,V,E)) = k_B \ln \Omega(N,V,E)\]

Boltzmann 常数 \(k_B\) 正是最后的比例系数。

如果我们将 Boltzmann 关系式代入之前得到的结果，不难得到

\[\left(\frac{\partial \ln \Omega}{\partial E}\right)_{N,V} = \frac{1}{k_B T}\]

\[\left(\frac{\partial \ln \Omega}{\partial V}\right)_{N,E} = \frac{P}{k_B T}\]

\[\left(\frac{\partial \ln \Omega}{\partial N}\right)_{V,E} = -\frac{\mu}{k_B T}\]

思考问题：

回忆刘维尔方程一章中的思考问题。 对于一元实函数，积分与求导顺序可以交换的条件是什么？ 试着在可导函数和可积函数中找出一些反例，从而证明可导或者可积并非实函数积分与求导顺序可以交换的充分条件。

对于多元函数呢？ 为什么在式 [eqn:Boltzmann_relation_derivation] 中可以交换积分与求导的顺序？ Dirac \(\delta\) 函数有什么性质？

经典位力定理

我们刚刚得到了微正则系综下的配分函数。 为了加深感觉，我们来进行一点应用。 我们打算证明经典统计力学中的一个重要定理，叫做 位力定理(virial theorem) ：

\[\expval{x_i \frac{\partial \mathcal H}{\partial x_j}} = k_B T \delta_{ij}\]

其中，\(x_i\) 是第 \(i\) 个广义坐标， \(\delta_{ij}\) 是 Kronecker delta。

由平均的定义，

\[\expval{x_i \frac{\partial \mathcal H}{\partial x_j}}
= \frac{M}{\Omega(N,V,E)} \int x_i \frac{\partial \mathcal H(\vec x)}{\partial x_j} \delta(\mathcal H(\vec x) – E) d\vec x\]

引入 Heaviside \(\theta\) 函数，则有

\[\delta(\mathcal H(\vec x) – E) = \frac{\partial \theta( E – \mathcal H(\vec x) )}{\partial E}\]

从而

\[\begin{aligned}
\expval{x_i \frac{\partial \mathcal H(\vec x)}{\partial x_j}}
& = \frac{M}{\Omega(N,V,E)} \int x_i \frac{\partial \mathcal H(\vec x)}{\partial x_j} \frac{\partial \theta( E – \mathcal H(\vec x) )}{\partial E} d\vec x \\
& = \frac{M}{\Omega(N,V,E)} \frac{\partial }{\partial E} \int x_i \frac{\partial \mathcal H(\vec x)}{\partial x_j} \theta( E – \mathcal H(\vec x) ) d\vec x \\
& = \frac{M}{\Omega(N,V,E)} \frac{\partial }{\partial E} \int_{\mathcal H(\vec x) < E} x_i \frac{\partial \mathcal H(\vec x)}{\partial x_j} d\vec x \\
& = \frac{M}{\Omega(N,V,E)} \frac{\partial }{\partial E} \int_{\mathcal H(\vec x) < E} x_i \frac{\partial (\mathcal H(\vec x) – E )}{\partial x_j} d\vec x \\
& = \frac{M}{\Omega(N,V,E)} \frac{\partial }{\partial E} \int_{\mathcal H(\vec x) < E} \left( \frac{\partial x_i (\mathcal H(\vec x) – E )}{\partial x_j} – \delta_{ij} (\mathcal H(\vec x) – E ) \right)d\vec x \\
& = \frac{M \delta_{ij}}{\Omega(N,V,E)} \frac{\partial }{\partial E} \int_{\mathcal H(\vec x) < E} – (\mathcal H(\vec x) – E ) d\vec x \\
& = \frac{M \delta_{ij}}{\Omega(N,V,E)} \int_{\mathcal H(\vec x) < E} \frac{\partial }{\partial E} ( E – \mathcal H(\vec x) ) d\vec x \\
& = \frac{M \delta_{ij}}{\Omega(N,V,E)} \int_{\mathcal H(\vec x) < E} d\vec x \\
& = \frac{M \delta_{ij}}{\Omega(N,V,E)} \int \theta(E – \mathcal H(\vec x)) d\vec x \\
& = \delta_{ij} \frac{\int \theta(E – \mathcal H(\vec x)) d\vec x}{\int \delta(E – \mathcal H(\vec x)) d\vec x} \\
& = \delta_{ij} \frac{\int \theta(E – \mathcal H(\vec x)) d\vec x}{\frac{\partial}{\partial E}\int \theta(E – \mathcal H(\vec x)) d\vec x} \\
& = \delta_{ij} \left( \frac{\partial}{\partial E} \ln \int \theta(E – \mathcal H(\vec x)) d\vec x \right)^{-1} \\
& \overset{N \to \infty}{=\joinrel=} \delta_{ij} \left( \frac{\partial}{\partial E} \ln \int \delta(E – \mathcal H(\vec x)) d\vec x \right)^{-1} \\
& = \delta_{ij} \left( \frac{\partial (\ln \Omega(N,V,E) – \ln M)}{\partial E} \right)^{-1} \\
& = \delta_{ij} \left( \frac{\partial S(N,V,E)}{\partial E} \right)^{-1} \\
& = \delta_{ij} k_B T \\
\end{aligned}\]

从而位力定理得证。

在上面的运算中，我们注意到了一个新的积分

\[\int \theta(E – \mathcal H(\vec x)) d\vec x\]

显然如果我们对其乘以一个实系数 \(C\) ，则这个积分实际上是另一个系综 \((N, V, E \leq E_{max})\) 的配分函数。

\[\Sigma(N,V,E) = C \int \theta(E – \mathcal H(\vec x)) d\vec x\]

这个系综被称为 均匀系综(uniform ensemble) 。 完全类似我们处理微正则系综配分函数时的方法，不难得到

\[C(N) = \frac{1}{N! h^{3N}}\]

即其与微正则系综配分函数系数的关系为

\[M = C E_0\]

有了位力定理，在构建那些系综平均能生成宏观热力学可观测量的微观相空间函数时就方便了许多。 我们不如现在就来尝试简单使用一下位力定理，考虑选取笛卡尔坐标系，且 \(x_i = x_j = p_i\) ，那么

\[\expval{ p_i \frac{\mathcal H}{\partial p_i} }
= \expval{ \frac{p_i^2}{m_i} }
= k_B T\]

立即可以得到一个熟悉的结果

\[\expval{ T_i }
= \expval{ \frac{1}{2} m_i v_i^2 }
= \expval{ \frac{p_i^2}{2 m_i} }
= \frac{1}{2} k_B T\]

\[\sum_{i=1}^{3N} \expval{ \frac{p_i^2}{2 m_i} }
= \sum_{i=1}^{3N} \frac{1}{2} k_B T
= \frac{3}{2} N k_B T\]

热平衡的条件

微正则系综能导出的另一个结果是热平衡的条件。 我们现在来看看这一点如何实现

自由粒子和理想气体

谐振子，harmonic baths

分子动力学(molecular dynamics)

有限差分法(finite difference method, FDM)

systems subject to holonomic constraints

TODO

正则系综

TODO

等压系综

巨正则系综

蒙特卡洛方法

自由能计算

量子系综和密度矩阵

费米-狄拉克统计和玻色-爱因斯坦统计(量子理想气体)

费恩曼的路径积分方法

经典含时统计力学

量子含时统计力学

朗之万方程和广义朗之万方程

临界现象

Appendix

using Plots
using LinearAlgebra

nparticles = 3

Linv = inv([
    1 -2 1
    -1 0 1
    1 1 1
])

r0 = 1.0:1.0:nparticles

eigvals = [ sqrt(3), 1, 0 ]

q = zeros(Float64, (nparticles,nparticles))

normal_mode_indices = permutedims(repeat(1:nparticles, 1, nparticles))

pyplot(size = (600,600), legend = false)

@gif for t = 0.0:0.15:10.0
    for i = 1:nparticles
        q[i,i] = sin(eigvals[i] * t)
    end
    
    plot(
            view(normal_mode_indices, :, 1),  
            r0 + Linv * view(q, :, 1) , 
            markershape=:circle, 
            markersize=20, 
            xlims = (0, nparticles + 1), 
            ylims = (0, nparticles + 1),
        )
    for i = 2:nparticles
        plot!(
            view(normal_mode_indices, :, i),
            r0 + Linv * view(q, :, i) , 
            markershape=:circle, 
            markersize=20, 
            xlims = (0, nparticles + 1), 
            ylims = (0, nparticles + 1),
        )
    end
end

!--------------------------------------------------------------------------
!>  pes for h2o with energy and gradients for one geometry, in f90 form.
!>  modified from the f77 code "h2opes.f" in E-PAPS, which is submitted by
!>  Harry Partridge and David W. Schwenke, J. Chem. Phys., 106, 4618 (1997) 
!>          z|     /y                +----+-------------------+         
!>           | H1 /                  | i\j|   x     y     z   |       
!>           |  \/                   +----+-------------------+        
!>           |  /\                   | H1 | (1,1) (1,2) (1,3) |        
!>           | /  O --- H2           | O  | (2,1) (2,2) (2,3) |        
!>         __|/__________            | H2 | (3,1) (3,2) (3,3) |      
!>           /           x           +----+-------------------+        
!> @param [out] v: potential in au
!> @param [out] dv(3,3): gradients (dvdx, dvdy, dvdz for each atom) in a.u. 
!> @param [out] ddv(3,3,3,3): hessian matrix, each element in a.u.
!>              ddv(i,j,k,l) = d(dv(i,j), cart(k,l))
!> @param [in]  cart(3,3): Cartesian coordinates for each atom in a.u.
!> @param [in]  flag: controls calculation, 
!               0 -- v; 1 -- v,dv; 2 -- v,dv,ddv
!--------------------------------------------------------------------------
subroutine h2opot1(v,dv,ddv,cart,flag)
    implicit none
    real(8) :: v
    real(8), dimension(3,3) :: dv
    real(8), dimension(3,3,3,3) :: ddv
    real(8), dimension(3,3) :: cart
    integer :: flag
 
    ! bond coordinates: bond lengths and cosine of bond angle H1-O-H1
    real(8) :: roh1, roh2, rhh, cost
 
    ! one-body potentials 
    real(8) :: voh1, voh2, vhh
 
    ! gradients with respect to valence coordinates (r1,r2,cos(theta))
    real(8) :: dv_dr(3), ddv_drdr(3,3)
 
    ! gradients of one-body potentials: dva/dr, dvb/dr
    real(8) :: dvoh1_droh1, dvoh2_droh2, dvhh_drhh
    real(8) :: ddvoh1_droh1droh1, ddvoh2_droh2droh2, ddvhh_drhhdrhh
 
    ! partial derivatives of rhh with respect to roh1, roh2, cost
    real(8) :: drhh_dr(3), ddrhh_drdr(3,3)
 
    ! partial derivatives of vhh: dvhh/dr
    real(8) :: dvhh_dr(3), ddvhh_drdr(3,3)
 
    ! partial derivatives of vc : dvc/dr (r=r1,r2,cos(theta))
    real(8) :: dvc_dr(3), ddvc_drdr(3,3)
 
    ! partial derivatives of r1, r2, cos(theta) : dr/dc
    real(8) :: dr_dc(3,3,3), ddr_dcdc(3,3,3,3,3)
 
    ! pes parameters    
    real(8), dimension(245) :: c5z, cbasis, ccore, crest
    integer, dimension(245,3) :: idx
    real(8) :: reoh,thetae,b1,roh,alphaoh,deoh,phh1,phh2
    real(8) :: coste
 
    ! scale factors      
    real(8) :: f5z,fbasis,fcore,frest
 
    ! constant values
    real(8), parameter :: pi = 3.1415926535897932384626433832795d0
 
    ! temporary values
    integer :: i, j, k, l, ir, jr, ifirst
    real(8) :: x1, x2, x3, ex, term
    real(8) :: vcterm, dterm, ddterm, dvcterm(3), ddvcterm(3,3)
    real(8) :: exa(2), exb, dexb_dr(2), ddexb_drdr(2,2)
 
    ! fmat is an array of terms in V^{c}, see Eq.(4)
    ! dfmat is 1st derivative of fmat with respect to r1,r2,cost
    ! ddfmat is 2nd derivative of fmat with respect to r1,r2,cost
    real(8), dimension(15,3) :: fmat, dfmat, ddfmat
 
    data (idx(i,1),i=1,245)/&
       & 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2,&
       & 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,&
       & 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,&
       & 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3,&
       & 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,&
       & 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,&
       & 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6,&
       & 6, 6, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5,&
       & 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 5, 5,&
       & 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,&
       & 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6,&
       & 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 9, 9,&
       & 9, 9, 9, 9, 9/
    data (idx(i,2),i=1,245)/&
       & 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,&
       & 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,&
       & 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2,&
       & 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3,&
       & 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,&
       & 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3,&
       & 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1,&
       & 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3,&
       & 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4,&
       & 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2,&
       & 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4,&
       & 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 1, 1,&
       & 1, 1, 1, 1, 1/
    data (idx(i,3),i=1,245)/&
       & 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15, 1, 2, 3, 4, 5,&
       & 6, 7, 8, 9,10,11,12,13,14, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,&
       &12,13, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13, 1, 2, 3, 4, 5,&
       & 6, 7, 8, 9,10,11,12, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12, 1,&
       & 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,&
       &11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8,&
       & 9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 1, 2, 3, 4, 5, 6, 7, 8,&
       & 9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9,&
       & 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2,&
       & 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6,&
       & 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3,&
       & 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2,&
       & 3, 4, 5, 6, 7/
    data (c5z(i),i=1,245)/&
       & 4.2278462684916D+04, 4.5859382909906D-02, 9.4804986183058D+03,&
       & 7.5485566680955D+02, 1.9865052511496D+03, 4.3768071560862D+02,&
       & 1.4466054104131D+03, 1.3591924557890D+02,-1.4299027252645D+03,&
       & 6.6966329416373D+02, 3.8065088734195D+03,-5.0582552618154D+02,&
       &-3.2067534385604D+03, 6.9673382568135D+02, 1.6789085874578D+03,&
       &-3.5387509130093D+03,-1.2902326455736D+04,-6.4271125232353D+03,&
       &-6.9346876863641D+03,-4.9765266152649D+02,-3.4380943579627D+03,&
       & 3.9925274973255D+03,-1.2703668547457D+04,-1.5831591056092D+04,&
       & 2.9431777405339D+04, 2.5071411925779D+04,-4.8518811956397D+04,&
       &-1.4430705306580D+04, 2.5844109323395D+04,-2.3371683301770D+03,&
       & 1.2333872678202D+04, 6.6525207018832D+03,-2.0884209672231D+03,&
       &-6.3008463062877D+03, 4.2548148298119D+04, 2.1561445953347D+04,&
       &-1.5517277060400D+05, 2.9277086555691D+04, 2.6154026873478D+05,&
       &-1.3093666159230D+05,-1.6260425387088D+05, 1.2311652217133D+05,&
       &-5.1764697159603D+04, 2.5287599662992D+03, 3.0114701659513D+04,&
       &-2.0580084492150D+03, 3.3617940269402D+04, 1.3503379582016D+04,&
       &-1.0401149481887D+05,-6.3248258344140D+04, 2.4576697811922D+05,&
       & 8.9685253338525D+04,-2.3910076031416D+05,-6.5265145723160D+04,&
       & 8.9184290973880D+04,-8.0850272976101D+03,-3.1054961140464D+04,&
       &-1.3684354599285D+04, 9.3754012976495D+03,-7.4676475789329D+04,&
       &-1.8122270942076D+05, 2.6987309391410D+05, 4.0582251904706D+05,&
       &-4.7103517814752D+05,-3.6115503974010D+05, 3.2284775325099D+05,&
       & 1.3264691929787D+04, 1.8025253924335D+05,-1.2235925565102D+04,&
       &-9.1363898120735D+03,-4.1294242946858D+04,-3.4995730900098D+04,&
       & 3.1769893347165D+05, 2.8395605362570D+05,-1.0784536354219D+06,&
       &-5.9451106980882D+05, 1.5215430060937D+06, 4.5943167339298D+05,&
       &-7.9957883936866D+05,-9.2432840622294D+04, 5.5825423140341D+03,&
       & 3.0673594098716D+03, 8.7439532014842D+04, 1.9113438435651D+05,&
       &-3.4306742659939D+05,-3.0711488132651D+05, 6.2118702580693D+05,&
       &-1.5805976377422D+04,-4.2038045404190D+05, 3.4847108834282D+05,&
       &-1.3486811106770D+04, 3.1256632170871D+04, 5.3344700235019D+03,&
       & 2.6384242145376D+04, 1.2917121516510D+05,-1.3160848301195D+05,&
       &-4.5853998051192D+05, 3.5760105069089D+05, 6.4570143281747D+05,&
       &-3.6980075904167D+05,-3.2941029518332D+05,-3.5042507366553D+05,&
       & 2.1513919629391D+03, 6.3403845616538D+04, 6.2152822008047D+04,&
       &-4.8805335375295D+05,-6.3261951398766D+05, 1.8433340786742D+06,&
       & 1.4650263449690D+06,-2.9204939728308D+06,-1.1011338105757D+06,&
       & 1.7270664922758D+06, 3.4925947462024D+05,-1.9526251371308D+04,&
       &-3.2271030511683D+04,-3.7601575719875D+05, 1.8295007005531D+05,&
       & 1.5005699079799D+06,-1.2350076538617D+06,-1.8221938812193D+06,&
       & 1.5438780841786D+06,-3.2729150692367D+03, 1.0546285883943D+04,&
       &-4.7118461673723D+04,-1.1458551385925D+05, 2.7704588008958D+05,&
       & 7.4145816862032D+05,-6.6864945408289D+05,-1.6992324545166D+06,&
       & 6.7487333473248D+05, 1.4361670430046D+06,-2.0837555267331D+05,&
       & 4.7678355561019D+05,-1.5194821786066D+04,-1.1987249931134D+05,&
       & 1.3007675671713D+05, 9.6641544907323D+05,-5.3379849922258D+05,&
       &-2.4303858824867D+06, 1.5261649025605D+06, 2.0186755858342D+06,&
       &-1.6429544469130D+06,-1.7921520714752D+04, 1.4125624734639D+04,&
       &-2.5345006031695D+04, 1.7853375909076D+05,-5.4318156343922D+04,&
       &-3.6889685715963D+05, 4.2449670705837D+05, 3.5020329799394D+05,&
       & 9.3825886484788D+03,-8.0012127425648D+05, 9.8554789856472D+04,&
       & 4.9210554266522D+05,-6.4038493953446D+05,-2.8398085766046D+06,&
       & 2.1390360019254D+06, 6.3452935017176D+06,-2.3677386290925D+06,&
       &-3.9697874352050D+06,-1.9490691547041D+04, 4.4213579019433D+04,&
       & 1.6113884156437D+05,-7.1247665213713D+05,-1.1808376404616D+06,&
       & 3.0815171952564D+06, 1.3519809705593D+06,-3.4457898745450D+06,&
       & 2.0705775494050D+05,-4.3778169926622D+05, 8.7041260169714D+03,&
       & 1.8982512628535D+05,-2.9708215504578D+05,-8.8213012222074D+05,&
       & 8.6031109049755D+05, 1.0968800857081D+06,-1.0114716732602D+06,&
       & 1.9367263614108D+05, 2.8678295007137D+05,-9.4347729862989D+04,&
       & 4.4154039394108D+04, 5.3686756196439D+05, 1.7254041770855D+05,&
       &-2.5310674462399D+06,-2.0381171865455D+06, 3.3780796258176D+06,&
       & 7.8836220768478D+05,-1.5307728782887D+05,-3.7573362053757D+05,&
       & 1.0124501604626D+06, 2.0929686545723D+06,-5.7305706586465D+06,&
       &-2.6200352535413D+06, 7.1543745536691D+06,-1.9733601879064D+04,&
       & 8.5273008477607D+04, 6.1062454495045D+04,-2.2642508675984D+05,&
       & 2.4581653864150D+05,-9.0376851105383D+05,-4.4367930945690D+05,&
       & 1.5740351463593D+06, 2.4563041445249D+05,-3.4697646046367D+03,&
       &-2.1391370322552D+05, 4.2358948404842D+05, 5.6270081955003D+05,&
       &-8.5007851251980D+05,-6.1182429537130D+05, 5.6690751824341D+05,&
       &-3.5617502919487D+05,-8.1875263381402D+02,-2.4506258140060D+05,&
       & 2.5830513731509D+05, 6.0646114465433D+05,-6.9676584616955D+05,&
       & 5.1937406389690D+05, 1.7261913546007D+05,-1.7405787307472D+04,&
       &-3.8301842660567D+05, 5.4227693205154D+05, 2.5442083515211D+06,&
       &-1.1837755702370D+06,-1.9381959088092D+06,-4.0642141553575D+05,&
       & 1.1840693827934D+04,-1.5334500255967D+05, 4.9098619510989D+05,&
       & 6.1688992640977D+05, 2.2351144690009D+05,-1.8550462739570D+06,&
       & 9.6815110649918D+03,-8.1526584681055D+04,-8.0810433155289D+04,&
       & 3.4520506615177D+05, 2.5509863381419D+05,-1.3331224992157D+05,&
       &-4.3119301071653D+05,-5.9818343115856D+04, 1.7863692414573D+03,&
       & 8.9440694919836D+04,-2.5558967650731D+05,-2.2130423988459D+04,&
       & 4.4973674518316D+05,-2.2094939343618D+05/
    data (cbasis(i),i=1,245)/&
       & 6.9770019624764D-04,-2.4209870001642D+01, 1.8113927151562D+01,&
       & 3.5107416275981D+01,-5.4600021126735D+00,-4.8731149608386D+01,&
       & 3.6007189184766D+01, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       &-7.7178474355102D+01,-3.8460795013977D+01,-4.6622480912340D+01,&
       & 5.5684951167513D+01, 1.2274939911242D+02,-1.4325154752086D+02,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00,-6.0800589055949D+00,&
       & 8.6171499453475D+01,-8.4066835441327D+01,-5.8228085624620D+01,&
       & 2.0237393793875D+02, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 3.3525582670313D+02, 7.0056962392208D+01,-4.5312502936708D+01,&
       &-3.0441141194247D+02, 2.8111438108965D+02, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00,-1.2983583774779D+02, 3.9781671212935D+01,&
       &-6.6793945229609D+01,-1.9259805675433D+02, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00,-8.2855757669957D+02,-5.7003072730941D+01,&
       &-3.5604806670066D+01, 9.6277766002709D+01, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 8.8645622149112D+02,-7.6908409772041D+01,&
       & 6.8111763314154D+01, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 2.5090493428062D+02,-2.3622141780572D+02, 5.8155647658455D+02,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 2.8919570295095D+03,&
       &-1.7871014635921D+02,-1.3515667622500D+02, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00,-3.6965613754734D+03, 2.1148158286617D+02,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00,-1.4795670139431D+03,&
       & 3.6210798138768D+02, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       &-5.3552886800881D+03, 3.1006384016202D+02, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 1.6241824368764D+03, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 4.3764909606382D+03, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 1.0940849243716D+03, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 3.0743267832931D+03, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00/
    data (ccore(i),i=1,245)/&
       & 2.4332191647159D-02,-2.9749090113656D+01, 1.8638980892831D+01,&
       &-6.1272361746520D+00, 2.1567487597605D+00,-1.5552044084945D+01,&
       & 8.9752150543954D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       &-3.5693557878741D+02,-3.0398393196894D+00,-6.5936553294576D+00,&
       & 1.6056619388911D+01, 7.8061422868204D+01,-8.6270891686359D+01,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00,-3.1688002530217D+01,&
       & 3.7586725583944D+01,-3.2725765966657D+01,-5.6458213299259D+00,&
       & 2.1502613314595D+01, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 5.2789943583277D+02,-4.2461079404962D+00,-2.4937638543122D+01,&
       &-1.1963809321312D+02, 2.0240663228078D+02, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00,-6.2574211352272D+02,-6.9617539465382D+00,&
       &-5.9440243471241D+01, 1.4944220180218D+01, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00,-1.2851139918332D+03,-6.5043516710835D+00,&
       & 4.0410829440249D+01,-6.7162452402027D+01, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 1.0031942127832D+03, 7.6137226541944D+01,&
       &-2.7279242226902D+01, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       &-3.3059000871075D+01, 2.4384498749480D+01,-1.4597931874215D+02,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 1.6559579606045D+03,&
       & 1.5038996611400D+02,-7.3865347730818D+01, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00,-1.9738401290808D+03,-1.4149993809415D+02,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00,-1.2756627454888D+02,&
       & 4.1487702227579D+01, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       &-1.7406770966429D+03,-9.3812204399266D+01, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00,-1.1890301282216D+03, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 2.3723447727360D+03, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00,-1.0279968223292D+03, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 5.7153838472603D+02, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00/
    data (crest(i),i=1,245)/&
       & 0.0000000000000D+00,-4.7430930170000D+00,-1.4422132560000D+01,&
       &-1.8061146510000D+01, 7.5186735000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       &-2.7962099800000D+02, 1.7616414260000D+01,-9.9741392630000D+01,&
       & 7.1402447000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00,-7.8571336480000D+01,&
       & 5.2434353250000D+01, 7.7696745000000D+01, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 1.7799123760000D+02, 1.4564532380000D+02, 2.2347226000000D+02,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00,-4.3823284100000D+02,-7.2846553000000D+02,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00,-2.6752313750000D+02, 3.6170310000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00,&
       & 0.0000000000000D+00, 0.0000000000000D+00/
    data reoh,thetae,b1,roh,alphaoh,deoh,phh1,phh2/0.958649d0,&
       &      104.3475d0,2.0d0,0.9519607159623009d0,2.587949757553683d0,&
       &      42290.92019288289d0,16.94879431193463d0,12.66426998162947d0/
    ! @todo this is different from reference, to be checked 
    data f5z,fbasis,fcore,frest/0.99967788500000d0,&
       &      0.15860145369897d0,-1.6351695982132d0,1d0/
 
    ! for the first time, convert parameters to a.u.       
    ! and save them for future use
    save
    data ifirst/0/
    if(ifirst.eq.0)then
       ifirst=1
       ! scale parameters empirically, which leads V^{5Z} to V^{emp}
       do i=1,245
          c5z(i)=f5z*c5z(i)+fbasis*cbasis(i)+fcore*ccore(i)+frest*crest(i)
       end do
       phh1=phh1*f5z
       deoh=deoh*f5z
       !
       reoh=reoh/0.529177249d0
       b1=b1*0.529177249d0*0.529177249d0
       do i=1,245
          c5z(i)=c5z(i)*4.556335d-6 
       end do
       coste=cos(thetae*pi/1.8d2)
       phh1=phh1*exp(phh2)
       phh1=phh1*4.556335d-6
       phh2=phh2*0.529177249d0
       deoh=deoh*4.556335d-6
       roh=roh/0.529177249d0
       alphaoh=alphaoh*0.529177249d0
       c5z(1)=c5z(1)*2d0
    end if
 
    ! convert cartesian coordinates to geometry (r1, r2, cos(theta), rhh)
    call convert_cart_to_geom(roh1, roh2, cost, rhh, &
       dr_dc, ddr_dcdc, cart, flag)
 
    ! calculate V^{emp} function in the form of V^{5Z}
    ! see Eq.(4) and (10) in [JCP106,4618]
    x1=(roh1 - reoh)/reoh
    x2=(roh2 - reoh)/reoh
    x3=cost - coste
 
    ! V^{b}(rhh)
    vhh=phh1*exp(-phh2*rhh)
 
    ! V^{a}(r1)
    ex=exp(-alphaoh*(roh1-roh))
    voh1=deoh*ex*(ex-2d0)
    exa(1)=ex
 
    ! V^{a}(r2)
    ex=exp(-alphaoh*(roh2-roh))
    voh2=deoh*ex*(ex-2d0)
    exa(2)=ex
 
    ! V^{c}(r1,r2,theta)
    ! power calculations are replaced by multiplications 
    ! and then stored in arrays
    fmat(1,1)=1d0
    fmat(1,2)=1d0
    fmat(1,3)=1d0
    do j=2,15
       fmat(j,1)=fmat(j-1,1)*x1 ! ((r1-re)/re)^{j-1}
       fmat(j,2)=fmat(j-1,2)*x2 ! ((r2-re)/re)^{j-1}
       fmat(j,3)=fmat(j-1,3)*x3 ! (ct-cte)^{j-1}
    end do
    vcterm=0d0
    do j=2,245
       term=c5z(j)*(fmat(idx(j,1),1)*fmat(idx(j,2),2)  &
          &        +fmat(idx(j,2),1)*fmat(idx(j,1),2)) &
          & *fmat(idx(j,3),3)
       vcterm=vcterm+term
    end do
 
    ! final result of potential
    exb = exp(-b1*((roh1-reoh)**2+(roh2-reoh)**2))  
    v=vcterm*exb + c5z(1) + voh1 + voh2 + vhh
 
    ! --- calculate dv ---
    if (flag > 0) then 
       ! drhh/dr, r=r1,r2,cost 
       drhh_dr(1) = (roh1-roh2*cost)/rhh
       drhh_dr(2) = (roh2-roh1*cost)/rhh
       drhh_dr(3) = -roh1*roh2/rhh 
 
       ! dvhh/drhh
       dvhh_drhh = -phh2*vhh
 
       ! dvoh1/droh1
       dvoh1_droh1=2d0*deoh*alphaoh*exa(1)*(1d0-exa(1))
 
       ! dvoh2/droh2
       dvoh2_droh2=2d0*deoh*alphaoh*exa(2)*(1d0-exa(2))
 
       ! calculate dvc/dr
       dexb_dr(1)=-2d0*b1*(roh1-reoh)*exb
       dexb_dr(2)=-2d0*b1*(roh2-reoh)*exb
 
       ! 1st derivative of fmat with respect to bond coordinates
       dfmat(1,1)=0d0
       dfmat(1,2)=0d0
       dfmat(1,3)=0d0
       do j=2,15
          dfmat(j,1)=dble(j-1)/reoh*fmat(j-1,1) ! (j-1)*((r1-re)/re)^{j-2}
          dfmat(j,2)=dble(j-1)/reoh*fmat(j-1,2) ! (j-1)*((r2-re)/re)^{j-2}
          dfmat(j,3)=dble(j-1)*fmat(j-1,3)      ! (j-1)*(ct-cte)^{j-2}
       enddo
 
       ! dvc/droh1
       dterm=0d0
       do j=2,245
          term=c5z(j)*(dfmat(idx(j,1),1)*fmat(idx(j,2),2)  &
             &       + dfmat(idx(j,2),1)*fmat(idx(j,1),2)) &	
             & *fmat(idx(j,3),3)
          dterm=dterm+term
       end do
       dvcterm(1)=dterm
       dvc_dr(1) = dvcterm(1)*exb + vcterm*dexb_dr(1)
 
       ! dvc/droh2
       dterm=0d0
       do j=2,245
          term=c5z(j)*(fmat(idx(j,1),1)*dfmat(idx(j,2),2)  &
             &       + fmat(idx(j,2),1)*dfmat(idx(j,1),2)) &	
             & *fmat(idx(j,3),3)
          dterm=dterm+term
       end do
       dvcterm(2)=dterm
       dvc_dr(2) = dvcterm(2)*exb + vcterm*dexb_dr(2)
 
       ! dvc/dcost
       dterm=0d0
       do j=2,245
          term=c5z(j)*(fmat(idx(j,1),1)*fmat(idx(j,2),2)  &
             &       + fmat(idx(j,2),1)*fmat(idx(j,1),2)) &	
             & *dfmat(idx(j,3),3)
          dterm=dterm+term
       end do
       dvcterm(3)=dterm
       dvc_dr(3)=dvcterm(3)*exb
 
       ! dv/dr, r=r1,r2,cost
       dv_dr(1)= dvoh1_droh1 + dvhh_drhh*drhh_dr(1) + dvc_dr(1)
       dv_dr(2)= dvoh2_droh2 + dvhh_drhh*drhh_dr(2) + dvc_dr(2)
       dv_dr(3)= dvhh_drhh*drhh_dr(3) + dvc_dr(3) 
 
       ! finally, calculate dv using the chain rule
       do i = 1,3
          do j = 1,3
             dterm=0d0
             do ir=1,3
                dterm=dterm+dv_dr(ir)*dr_dc(ir,i,j)
             end do
             dv(i,j)=dterm 
          end do
       end do
 
       ! --- calculate ddv ---
       if (flag > 1) then  
          ! d(drhh/dr)/dr, r=r1,r2,cost
          ddrhh_drdr(1,1) = (1d0-drhh_dr(1)*drhh_dr(1))/rhh
          ddrhh_drdr(1,2) = (-cost-drhh_dr(1)*drhh_dr(2))/rhh
          ddrhh_drdr(1,3) = (-roh2-drhh_dr(1)*drhh_dr(3))/rhh
          ddrhh_drdr(2,1) = ddrhh_drdr(1,2)
          ddrhh_drdr(2,2) = (1d0-drhh_dr(2)*drhh_dr(2))/rhh
          ddrhh_drdr(2,3) = (-roh1-drhh_dr(2)*drhh_dr(3))/rhh
          ddrhh_drdr(3,1) = ddrhh_drdr(1,3)
          ddrhh_drdr(3,2) = ddrhh_drdr(2,3)
          ddrhh_drdr(3,3) = (-drhh_dr(3)*drhh_dr(3))/rhh
 
          ! 2nd derivative of Va(r1)
          ddvoh1_droh1droh1=2d0*deoh*alphaoh*alphaoh*exa(1)*(2d0*exa(1)-1d0)
 
          ! 2nd derivative of Va(r2)
          ddvoh2_droh2droh2=2d0*deoh*alphaoh*alphaoh*exa(2)*(2d0*exa(2)-1d0)
 
          ! 2nd derivatives of Vb(rhh)
          ddvhh_drhhdrhh=-phh2*dvhh_drhh
          do i = 1,3
             do j = 1,3
                ddvhh_drdr(i,j) = drhh_dr(i)*ddvhh_drhhdrhh*drhh_dr(j) &
                   + dvhh_drhh*ddrhh_drdr(i,j)
             end do
          end do
 
          ! 2nd derivative of exb with respect to bond coordinates
          ! exb = exp(-b1*((roh1-reoh)**2+(roh2-reoh)**2))  
          ddexb_drdr(1,1) = -2d0*b1*(exb + (roh1-reoh)*dexb_dr(1))
          ddexb_drdr(1,2) = -2d0*b1*(roh2-reoh)*dexb_dr(1)
          ddexb_drdr(2,1) = ddexb_drdr(1,2)
          ddexb_drdr(2,2) = -2d0*b1*(exb + (roh2-reoh)*dexb_dr(2))
 
          ! d(d(vcterm)/dr)/dr
          ddfmat(1,:)=0d0
          ddfmat(2,:)=0d0
          do j = 3,15
             ddfmat(j,1)=dble(j-1)/reoh*dfmat(j-1,1) !(j-1)/re*(j-2)*((r1-re)/re)^{j-3}
             ddfmat(j,2)=dble(j-1)/reoh*dfmat(j-1,2) !(j-1)/re*(j-2)*((r2-re)/re)^{j-3}
             ddfmat(j,3)=dble(j-1)*dfmat(j-1,3)      !(j-1)   *(j-2)*(ct-cte)^{j-3}
          end do
 
          ! --- d(dvc/dr)/dr, r=r1,r2,cost ---
          ! d(dvc/droh1)/droh1
          ddterm=0d0
          do j=2,245
             term=c5z(j)*(ddfmat(idx(j,1),1)*fmat(idx(j,2),2)  &
                &       + ddfmat(idx(j,2),1)*fmat(idx(j,1),2)) &	
                & *fmat(idx(j,3),3)
             ddterm=ddterm+term
          end do
          ddvc_drdr(1,1)=ddterm*exb+2d0*dvcterm(1)*dexb_dr(1)+vcterm*ddexb_drdr(1,1)
 
          ! d(dvc/droh1)/droh2
          ddterm=0d0
          do j=2,245
             term=c5z(j)*(dfmat(idx(j,1),1)*dfmat(idx(j,2),2)  &
                &       + dfmat(idx(j,2),1)*dfmat(idx(j,1),2)) &	
                & *fmat(idx(j,3),3)
             ddterm=ddterm+term
          end do
          ddvc_drdr(1,2)=ddterm*exb+dvcterm(1)*dexb_dr(2) &
             +dvcterm(2)*dexb_dr(1)+vcterm*ddexb_drdr(1,2)
 
          ! d(dvc/droh1)/dcost
          ddterm=0d0
          do j=2,245
             term=c5z(j)*(dfmat(idx(j,1),1)*fmat(idx(j,2),2)  &
                &       + dfmat(idx(j,2),1)*fmat(idx(j,1),2)) &	
                & *dfmat(idx(j,3),3)
             ddterm=ddterm+term
          end do
          ddvc_drdr(1,3)=ddterm*exb + dvcterm(3)*dexb_dr(1)
 
          ! d(dvc/droh2)/droh1
          ddvc_drdr(2,1)=ddvc_drdr(1,2)
 
          ! d(dvc/droh2)/droh2
          ddterm=0d0
          do j=2,245
             term=c5z(j)*(fmat(idx(j,1),1)*ddfmat(idx(j,2),2)  &
                &       + fmat(idx(j,2),1)*ddfmat(idx(j,1),2)) &	
                & *fmat(idx(j,3),3)
             ddterm=ddterm+term
          end do
          ddvc_drdr(2,2) = ddterm*exb +2d0*dvcterm(2)*dexb_dr(2)+vcterm*ddexb_drdr(2,2)
 
          ! d(dvc/droh2)/dcost
          ddterm=0d0
          do j=2,245
             term=c5z(j)*(fmat(idx(j,1),1)*dfmat(idx(j,2),2)  &
                &       + fmat(idx(j,2),1)*dfmat(idx(j,1),2)) &	
                & *dfmat(idx(j,3),3)
             ddterm=ddterm+term
          end do
          ddvc_drdr(2,3)=ddterm*exb + dvcterm(3)*dexb_dr(2)
 
          ! d(dvc/dcost)/droh1
          ddvc_drdr(3,1)=ddvc_drdr(1,3)
 
          ! d(dvc/dcost)/droh2
          ddvc_drdr(3,2)=ddvc_drdr(2,3)
 
          ! d(dvc/dcost)/dcost
          ddterm=0d0
          do j=2,245
             term=c5z(j)*(fmat(idx(j,1),1)*fmat(idx(j,2),2)  &
                &       + fmat(idx(j,2),1)*fmat(idx(j,1),2)) &	
                & *ddfmat(idx(j,3),3)
             ddterm=ddterm+term
          end do
          ddvc_drdr(3,3)=ddterm*exb
 
          ! d(dv/dr)/dr (r=r1,r2,cost)
          ddv_drdr(1,1) = ddvoh1_droh1droh1 + ddvhh_drdr(1,1) + ddvc_drdr(1,1)
          ddv_drdr(1,2) = ddvhh_drdr(1,2) + ddvc_drdr(1,2)
          ddv_drdr(1,3) = ddvhh_drdr(1,3) + ddvc_drdr(1,3)
 
          ddv_drdr(2,1) = ddv_drdr(1,2)
          ddv_drdr(2,2) = ddvoh2_droh2droh2 + ddvhh_drdr(2,2) + ddvc_drdr(2,2)
          ddv_drdr(2,3) = ddvhh_drdr(2,3) + ddvc_drdr(2,3)
 
          ddv_drdr(3,1) = ddv_drdr(1,3)
          ddv_drdr(3,2) = ddv_drdr(2,3)
          ddv_drdr(3,3) = ddvhh_drdr(3,3) + ddvc_drdr(3,3)
 
          ! finally, calculate the 2nd derivatives d(dv/dc_{i,j})/dc_{k,l} with the chain rule
          do i=1,3 ! atom H1,O,H2
             do j=1,3 ! x,y,z
                do k=1,3 ! atom H1,O,H2
                   do l=1,3 ! x,y,z
                      ddterm=0d0
                      do ir=1,3 ! dv/dr_ir
                         term=0d0
                         do jr=1,3  ! d(dv/dr_ir)/dr_jr
                            term=term + ddv_drdr(ir,jr)*dr_dc(jr,k,l)
                         end do 
                         ddterm=ddterm+ term*dr_dc(ir,i,j) + dv_dr(ir)*ddr_dcdc(ir,i,j,k,l)
                      end do 
                      ddv(i,j,k,l)=ddterm
                   end do 
                end do 
             end do 
          end do ! end of calculate ddv loop
       end if ! end of if calculate ddv 
    end if  ! end of if calculate dv
    return
 end subroutine h2opot1