Manual
User Manual:
Open the PDF directly: View PDF 
.
Page Count: 13
Page | 1
INTPOL 使用手册
三原子通用插值程序 (多体展开/样条/再生核/混合方法)
版本 1.1，2015 年12 月19 日。
1
简介
INTPOL 用于从 ab initio 离散原子坐标和能量数据，通过多种插值方法，得到精确的连续的势
能面。INTPOL 可生成势能面的 Fortran 语言源程序，并且可提供势能面的解析一阶和二阶导数。
支持多种原始数据的坐标：内坐标、内坐标和角度坐标的组合、雅克比坐标。可以采用可选的多
体展开方法降低误差。可以对每个坐标指定不同的插值方法。
对于生成的势能面，可以进一步计算插值在原始数据点的均方根误差。可以采用二次最速下
降法搜索势能面中的最小值点和过渡态。对给定点可以计算简谐振动频率和零点能。
下面是使用时需要输入的指令概述。详细的解释请参考下一节。
第一次使用：
unzip INTPOL.1.1.zip
g++ -O3 intpol.cpp -o intpol.x
生成新的势能面程序并计算系数：
cp template.f90 intpol.x <工作目录>
cd <工作目录>
./intpol.x <输入文件>
cd <输出目录>
./compile.sh
./prep.x
检查插值误差：
./test.x
进行几何优化 (寻找最小值和过渡态) 和频率计算：
./run.x
Page | 2
2
详细运行步骤
解压源程序压缩包。
unzip INTPOL.1.1.zip
有两个主要文件：
intpol.cpp 是驱动文件，用于读取输入文件和生成预编译头文件。
template.f90 是主程序模板文件，包含插值程序相关的所有 fortran 代码片段。注意不能直接
编译该文件，因为该文件中很多代码片段是相冲突的。
下面是运行方法：
2.1
编译
intpol.cpp
编译 intpol.cpp 得到可执行文件 intpol.x。注意该程序和具体的需要插值的体系无关，所以只
需要编译一次。以后无论更改输入文件还是更改体系，都无需重新编译。输入如下命令：
g++ -O3 intpol.cpp -o intpol.x
2.2
复制
intpol.x
和
template.f90
到工作目录
工作目录是自己选定的任意目录。下面以附带的用于示例的 HClLi 体系为例，工作目录是
examples/hclli 目录。因此需要将 intpol.x 和 template.f90 复制到 examples/hclli ：
cp template.f90 intpol.x examples/hclli
下面的操作都在工作目录进行，所以转到工作目录：
cd examples/hclli
2.3
准备输入文件及相关
ab initio
数据
输入文件描述采用何种插值方法及各种参数细节。输入文件的具体语法将在第 2节进行解
释。这里 examples/hclli 目录中已经准备好一个输入文件 input-r2s.txt，这个输入文件指定的是二
维再生核+一维样条的混合方法。
这里已经准备好 HClLi 体系的势能离散点数据和三个二体势数据，共 4个文件，存放在
orig_data 文件夹内。当用户自己准备这些离散点数据时，文件名和路径都可以任意指定，只要在
input-r2s.txt 输入文件中写明这些数据文件所在位置的相对路径即可 (相对工作目录的路径)。
Page | 3
2.4
读取输入文件
接下来读取并自动分析输入文件：
./intpol.x input-r2s.txt
这一步会生成许多输出文件，它们会被存放在一个指定的输出目录 hclli-r2s。该目录的名字是
在input-r2s.txt 中给定的。下面进入到这个输出目录 hclli-r2s：
cd hclli-r2s
2.5
编译
在hclli-r2s 会有一个自动生成的脚本 compile.sh 用于编译。运行该脚本得到插值程序和可执
行程序：
./compile.sh
2.6
检查结果
编译完会得到 prep.f90，test.f90，calc.f90，run.f90，prep.x，test.x，run.x 几个文件。
prep.f90 这是用于计算系数的 fortran 源程序。无论采用何种插值方法，在进行计算之前，都
需要先计算系数。计算系数将产生系数文件 *-prep.txt。该文件可以重复利用，只要不再更改势能
面参数，就不用重新计算系数。
test.f90 这是用于检查势能面插值效果 (计算误差) 的fortran 源程序。
calc.f90 这是势能面求值程序，可以和别的外部程序 (比如动力学程序) 一起编译。
run.f90 这是几何优化和频率计算程序。只有输入文件中要求计算频率或进行几何优化时，才
会产生该文件。该程序会将优化得到的坐标和频率输出到屏幕上。
prep.x 和 test.x 分别是编译 prep.f90 和 test.f90 得到的，可以先执行 ./prep.x 计算系数，然后
执行 ./test.x 计算势能面误差。如果不关心误差也可以不执行 ./test.x。
run.x 是编译 run.f90 得到的。如果想直接从屏幕上读取几何优化和频率计算结果，则执行
./run.x
如果想保存该结果供以后引用，则可以执行
./run.x > run.out
Page | 4
除此之外还有一些 txt 文件会被程序引用。每个程序并不需要所有这些 txt 文件。下面*表示
特定的体系名称。
*-pes-main.txt 是经过预处理的主要的ab initio 数据。
*-pes-mbrab.txt，*-pes-mbrbc.txt，*-pes-mbrac.txt 是经过预处理的二体势 ab initio 数据。
*-pes-param.txt 是包含距离型再生核插值指数和再生核正规化参数的参数文件。如果指定了
几何优化或者频率计算，该文件还会包含几何优化的起始点和其他参数。该文件可能包含几个
namelist。第一个 namelist 是rkhs，包含再生核相关的参数。
在rkhs 这个 namelist 里面，第一个变量 r_exp 是再生核插值指数，第二个变量 r_alpha 是正
规化参数。每个变量第一个数对应多体展开的再生核。第二三四个数分别对应第一二三个坐标的
再生核。如果没有使用多体展开或对应坐标没有使用再生核，对应的参数会是 0并且不起作用。
可以修改该文件而不需重新编译程序。但是修改该文件以后，一般需要重新进行系数计算。
如果指定了几何优化或者频率计算，会有 namelist point1, point2, …, opt1, opt2 等等。
在point 系列的namelist 里面，只有一个变量 x，有三个分量，表示三个坐标。这三个坐标
指定了要进行频率计算或者几何优化的初始点。如果想在其他的点进行几何优化或者频率计算，
只需修改该参数文件而无需重新编译程序。point 系列的 namelist 出现顺序和次数是和输入文件
中的 Point 指令出现的顺序和次数一一对应的。
在opt 系列的 namelist 里面，包含了每次进行几何优化的参数。iter 指定进行几何优化的最
大迭代步数。root 指定优化到最小值点 (1) 或者优化到过渡态 (2)。step_len 指定每次迭代的最大
步长。step_conv 和 f_conv 分别是判断迭代收敛的步长和力 (即能量导数) 阈值。只有这两个条件
同时满足才会判定收敛。可以直接修改这些参数然后重新执行程序而无需重新编译。
*-pes-prep.txt 是计算系数得到的系数文件。
*-pes-manybody.txt 是包含二体势 ab initio 数据和对应系数的系数文件。
prep.f90 在执行时需要*-pes-main.txt，*-pes-param.txt。如果使用了多体展开，还需要 *-pes-
mbrab.txt，*-pes-mbrbc.txt，*-pes-mbrac.txt。该程序执行完会产生*-pes-prep.txt。如果使用了多
体展开，还会产生*-pes-manybody.txt。
calc.f90/test.f90/run.f90 在执行时需要*-pes-main.txt，*-pes-param.txt，*-pes-prep.txt。如果
使用了多体展开，还需要*-pes-manybody.txt。
在 calc.f90 中，调用势能面子程序的方法是：
call intpol_main(0) 初始化。
call intpol_main(1, x, y, z, v) 计算坐标 (x, y, z) 的势能值 v。
Page | 5
call intpol_main(3, x, y, z, 0d0, dv) 计算坐标 (x, ,y, z) 的势能对坐标一阶导数 dv(1:3)。各分量分
别为 dx, dy, dz。
call intpol_main(6, x, y, z, 0d0, dv) 计算坐标 (x, y, z) 的势能对坐标二阶导数 dv(1:6)。各分量分
别为 dxdx，dxdy，dxdz，dydy，dydz，dzdz。
有时在别的程序直接写 call intpol_main 时会出现一些问题。这时可以在对应的 subroutine 的
implicit none 语句后面加上下面的 interface。之后便可以成功调用 intpol_main。
interface
subroutine intpol_main(do_type, x, y, z, v, dv)
implicit none
integer :: do_type
real(8), optional :: x, y, z, v, dv(0:do_type - 1)
end subroutine intpol_main
end interface
3
输入语法
输入文件中可以使用感叹号 ! 加入注释语句。参数赋值的格式为：
参数名 = 参数值;
当要求多个参数值时，各参数值可以用空格或 Tab 隔开。每条赋值句用分号;结尾，多个赋值
可以写在一行上。除了 Point, Frequencies 这些对应几何优化和频率计算的指令以外，参数之间没
有顺序要求。
有些参数需要放在对应的组中。一个组可以包含一些特定的参数。组的格式为：
{ 组名: 参数名 = 参数值; 参数名 = 参数值; … }
不需要放在组中的参数称为全局参数。
允许的组名为：ManyBody, Spline, RKHS, Optimization。分别用于指定多体展开，样条插值和
再生核插值的细节参数。如果没有使用某些插值方法，则不用在输入文件中包括这些组。如果对
多个插值步骤使用了某一方法，可能需要使用对应的组多次。除了 Optimization 属于对应几何优
化的指令外，其他组之间没有顺序要求。
组名和参数名不区分大小写。下面依次介绍所有的参数。
Page | 6
3.1
全局参数
3.1.1 ManyBody
可能的值：True 或 False。表示是否使用多体展开。
如果使用多体展开，则需提供独立的二体数据。虽然在旧版程序中二体数据可以自动从三体
数据得到。二体数据要尽量光滑。二体数据也可以单独计算二体势得到，但是要注意这时需要给
每个二体数据单独指定 VCut 和 VMin 参数 (在ManyBody 组中)。
多体展开的原理是减少三体项的变化幅度，从而减小整体误差。多体展开是否可以改进样条
方法的插值效果尚不明确。由于再生核对函数渐近区特性有限制，当插值方法包括再生核部分
时，强烈建议启用多体展开。
3.1.2 AtomName
值：三个字符串，用空格隔开。依次对应 A、B、C原子的元素名称或自定义名称。
如果进行频率计算或几何优化，这些名称将会被用到频率计算和几何优化的输出中。
3.1.3 AtomMass
值：三个实数，可以以 amu 为单位。依次对应 A、B、C原子的质量。
原子质量用于振动频率计算，对势能面插值不影响。为了得到准确的振动频率，必须使用
amu (相对原子质量) 单位。如果不进行频率计算，则可以任意选择单位。如果使用雅克比坐标，
将会使用该参数计算双原子质心，用以实现坐标变换。
3.1.4 Coordinate
值：三个特定格式字符串，见下面的举例。表示三个插值坐标的含义。
如果是距离坐标，用 R开头；如果是角度坐标，用 X开头。后面使用 ABCM 代表原子。M用
于雅克比坐标所需的中点位置。例如
RAB RBC RAC 表示选取了三个内坐标
RAB RBC XABC 表示选取了两个距离坐标及其夹角
RAM RBC XAMC 表示选取 A+BC 构形的雅克比坐标
XBMA RBM RAC 表示选取 B+AC 构形的雅克比坐标
距离坐标的单位可以是任意的，但是需要统一单位，通常可以取原子单位。不同距离单位对
再生核插值会影响结果。角度单位必须是度(0~360)，通常考虑到对称性只需取遍 0~180。如果要
进行频率计算，要求距离单位必须为原子单位。
Page | 7
3.1.5 Method
值：下面各种给定的值之一。表示插值方法。这里给出对每个坐标依次采用哪种方法进行插
值。只需给出方法名称。每种方法的具体参数在后面的小节中指定。通常，S表示样条，R表示
再生核。
推荐的方法：SSS、R2S。可以得到准确的一阶和二阶导数。
可能的取值有下面几种组合：
SSS 三维立方样条 (即纯三维样条，属于分三步的分步插值。部分区域对称性可能校差，准确
性较高。步数：3)
R3 对三个坐标整体使用三维再生核 (即纯再生核。整体对称性好，系数计算很慢且需要大量
内存，准确性可能不高，可用于完全非网格的点分布，但并不建议使用高度非网格分布，因其会
降低势能面求值速度。步数：1)
RRR 对三个坐标分别使用三个一维再生核 (即分步插值的再生核。对称性可能校差，速度
快，要求网格点分布。导数会很不准确，可能导致几何优化难以收敛。步数：3)
R2R 对前两个坐标整体使用二维再生核，对最后一个坐标使用一维再生核 (与纯再生核相比，
降低了整体对称性，提高了系数计算速度，提高了准确性。前两个坐标可以使用非网格点分布。
导数不准确。步数：2)
RR2 对第一个坐标使用一维再生核，对后两个坐标整体使用二维再生核 (求值速度会非常慢，
导数不准确。步数：2)
R2S 对前两个坐标整体使用二维再生核，对最后一个坐标使用一维样条 (即传统的混合方法，
属于分两步的分步插值。整体对称性与纯三维样条相比更高，计算速度快，准确性较高。可以给
出准确的解析一阶和二阶导数。步数：2)
SR2 对第一个坐标使用一维样条，对后两个坐标整体使用二维再生核 (步数：2)
SRR 对三个坐标依次使用样条、再生核、再生核 (分三步的混合分步插值。步数：3)
RSR 对三个坐标依次使用再生核、样条、再生核 (分三步的混合分步插值。步数：3)
RRS 对三个坐标依次使用再生核、再生核、样条 (分三步的混合分步插值。步数：3)
RSS 对三个坐标依次使用再生核、样条、样条 (分三步的混合分步插值。步数：3)
SRS 对三个坐标依次使用样条、再生核、样条 (分三步的混合分步插值。步数：3)
SSR 对三个坐标依次使用样条、样条、再生核 (分三步的混合分步插值。步数：3)
Page | 8
对上述所有方法，即便说明需要网格分布，实际上对于非网格分布算法也能运行，但是视具
体情况可能一定程度上会降低插值效果 (当然如果非网格部分的点选取合适也可以提升差插值效
果)。可以使用某种非网格分布，然后计算插值误差，同时检查增加的点是否能提升频率的准确
性，就可以知道该种分布是否可以改善插值效果。
当需要求一阶导数或二阶导数时，对多步方法，若后续步骤为再生核，可能导数的效果不
好，因为再生核的形式可能不适合对导数函数进行插值。但是该问题不影响势能值的效果。当分
步插值中，第二或第三步为再生核时，求值速度会显著变慢。
3.1.6 MainFile
值：文件相对路径。指定总能量文件名。
该文件包含总能量，如果使用了多体展开，也不需要用户手动从中减去二体势。
3.1.7 MainColumn
值：四个整数。分别为第一、二、三个坐标和能量值所在列的序号。列的序号从 1开始。
注意，尽管在不同的插值方法中可以对每一坐标的处理方式做调整，但是对于分步插值而
言，无论选择何种插值方法，排在后面的坐标一定会后处理。而不同的处理顺序也会对插值效果
造成影响。所以，有时可以通过更改此参数，调整数据列中三个坐标的给出顺序，来改善插值效
果。当修改此顺序时，注意 参数“Coordinate”也应做相应调整，以保持对应。
3.1.8 OutputDir
值：文件相对路径。指定输出文件夹名称。
经过处理的数据文件、参数文件、程序生成脚本、头文件等会输出到此文件夹。程序模板文
件会复制到此文件夹。
3.1.9 VCut
值：实数。能量截断值。能量 (减去 VMin 之前的) 大于该值的点的能量将被设置为该值。
3.1.10 VMin
值：实数。能量都会被减去该值。该值只是一个能量参考点的更改，不影响插值效果。可以
设为 0或能量最小值或者任何自定的参考点。更改参考点后，插值程序输出的数值会产生平移，
即都以该值为参考点。
3.2
组：
ManyBody
当 ManyBody 参数的值为 True 时应提供该组。如果不使用多体展开，则可以不使用此组。
3.2.1 Method
值：R 或 S。指定二体势插值方法。S表示样条，R表示再生核。
Page | 9
3.2.2 E1
值：实数。一体势能量 (减去全局参数 VMin 之前的)。该能量一般是三个原子相距无穷远时的
能量值。
3.2.3 VCut
值：三个实数。分别表示 rab、rbc、rac 的二体势能量截断值 (减去各自的VMin 之前的)。如
果二体势原始数据的能量参考点和总能量参考点是一样的，那么该值都设为和全局参数 VCut 一样
即可。
3.2.4 VMin
值：三个实数。rab、rbc、rac 的二体势能量都会被减去该值。如果二体势原始数据的能量参
考点和总能量参考点是一样的，那么该值都设为和全局参数 VMin 一样即可。
3.2.5 RABFile
值：文件相对路径。指定 rab 二体势能量文件名。
3.2.6 RBCFile
值：文件相对路径。指定 rbc 二体势能量文件名。
3.2.7 RACFile
值：文件相对路径。指定 rac 二体势能量文件名。
3.2.8 RABColumn
值：两个整数。分别为 RABFile 文件中 rab 坐标和能量值所在列的序号。列的序号从 1开始。
3.2.9 RBCColumn
值：两个整数。分别为 RBCFile 文件中 rbc 坐标和能量值所在列的序号。列的序号从 1开始。
3.2.10 RACColumn
值：两个整数。分别为 RACFile 文件中 rac 坐标和能量值所在列的序号。列的序号从 1开始。
3.3
组：
Spline
该组指定样条插值的具体参数。视具体情况，可以有多个。
3.3.1 Dim
值：一个整数，0, 1, 2, 或3。
0 表示该组描述多体展开的插值参数。1 2 3 分别表示该组描述第一二三个坐标的插值参数。
3.3.2 Boundary
值：natural 或clamp。表示样条插值的边界条件。
natural 自然样条边界条件。要求边界处二阶导数为零。当该坐标为距离坐标时建议选用。
Page | 10
clamp 钳制边界条件。要求边界处一阶导数为零。当该坐标为角度坐标时建议选用。
3.3.3 Fast
值：True 或False。是否使用快速样条。
只有该小节对应维度为 2或3时，才能使用快速样条。第一个坐标的求值理论上不能被加
速。快速样条会造成较大的浮点误差，导致与 ab initio 固定小数位数的随机偏离，该偏离一般对
计算影响不大。如果不使用快速样条，将降低求值速度。也可以部分使用快速样条。
注意：快速样条在目前版本的程序中尚未实现，所以该值必须设为 False。
3.4
组：
RKHS
该组指定再生核插值的具体参数。视具体情况，可以有多个。
3.4.1 Dim
值：一个整数，0, 1, 2, 或3。
0 表示该组描述多体展开的插值参数。1 2 3 分别表示该组描述第一二三个坐标的插值参数。
3.4.2 Type
值：distancelike 或anglelike。表示再生核类型。
distancelike 距离型再生核。当该坐标为距离坐标时建议选用。
anglelike 角度型再生核。当该坐标为角度坐标时建议选用。
3.4.3 Exp
值：实数。表示距离型再生核的指数。
该指数越大，表示该坐标对应的两个原子之间相互作用随距离增大衰减越快。只对距离型再
生核需要指定，通常取为 6。对角度型再生核应设为 0，表示此参数不适用。
3.4.4 Regularization
值：实数。表示正规化参数。
用于降低求解系数的病态性。该系数越大，方程病态性越低，拟合出的函数越光滑，但是与
ab initio 数据的偏离越高。该系数越小，拟合出的函数越不光滑，与 ab initio 的偏离越小。当数据
点不多时，该参数取为 0也是可以的。对二体项的插值建议取为 0。其他情况通常取为 1E-14 左
右。当采用二维或二维以上的再生核时，对应的各个维度的正规化参数必须一样。
Page | 11
3.5
频率和几何优化
频率计算和几何优化通过指令 Frequencies/Freq, Point 和 Optimization/Opt 来实现。其中
Frequencies 和 Point 是参数名，Optimization 是组名，这意味着 Frequencies 无需放在大括号中，
而 Optmization 必须放在大括号中并且加冒号。
通常，频率需要在几何优化后的点 (最小值或过渡态) 进行计算。在其他一般点进行频率计算
没有意义，因为不满足简谐近似的条件 (一阶导数为零)。所以一般的执行顺序是先指定一个初始
点，然后优化到该点附近的过渡态或者最小值，然后在优化的点计算频率。对应的输入指令如
下：
Point = 1.0 1.0 90.0;
{ Optimization: Type = minimum; }
Frequencies;
其中 Point 指令指定了一个初始点，Optimization: Type = minimum 指令将寻找该点附近的最
小值，最后 Frequencies 指令计算该最小值点的频率。不带参数值 Point 指令将计算当前点的能量
和导数，比如
Point = 1.0 1.0 90.0;
{ Optimization: Type = minimum; }
Point;
的意思是，指定初始点为 (1, 1, 90)，然后优化到附近的最小值点，最后计算最小值点的能量
和导数。最后一步计算出的导数理论上应该非常接近零。
除此之外，也可进行多个不同点的几何优化和频率计算，只要按顺序执行这些指令：
Point = 1.0 1.0 90.0;
{ Optimization: Type = minimum; }
Frequencies;
Point = 2.0 2.0 30.0;
{ Optimization: Type = transition; }
Frequencies;
下面是每个指令的详细解释。
3.5.1 Point
值：没有值或给出三个坐标。若给出坐标，三个坐标依次对应 Coordinate 参数指定的含义。
用于指定几何优化或者频率计算的点，并在该点进行单点能计算和导数计算 (基于构造的势能
面)。如果不给出坐标，则在当前点进行单点能和导数计算。在任何 Frequencies 或者
Page | 12
Optimization 指令之前，必须有至少一个带坐标值的Point 指令。因为无法在不给定点的情况下计
算频率或进行几何优化。
3.5.2 Frequencies/Freq
值：无值。用于在当前点计算振动频率和零点能。结果在执行产生的 run.x 程序时输出。
3.6
组：
Optimization/Opt
该组指定几何优化的具体参数。每个参数都有对应的默认值。如果采取默认值，则可以省略
该参数。几何优化的起始点通过在该组之前的 Point 指令指定。若要省略所有参数进行几何优
化，可以在输入文件中写：
{ Optimization: }
下面给出该组中可以包含的参数及其默认值。这些参数在编译后都可以在产生的参数文件中
更改。几何优化可能因为不收敛而失败。如果几何优化失败，将不会执行该指令后的所有指令而
直接退出。几何优化会转换到内坐标中迭代，完成后后转换回原坐标。
3.6.1 Type
值：transition 或 minimum。默认值：minimum。
minimum 指定优化到最小值点。transition 指定优化到过渡态。
3.6.2 Iter
值：整数。默认值：100。
指定几何优化的最大迭代次数。如果步长过小通常需要较大的迭代次数。如果经过最大迭代
次数而不能收敛，将会判定几何优化失败。
3.6.3 StepLength
值：实数。默认值：0.1。
指定几何优化的最大步长。步长过大会降低精度。步长过小可能导致收敛慢。
3.6.4 StepConv
值：实数。默认值：1e-8。指定几何优化的步长收敛阈值。小于此步长并且满足其他收敛条
件，将停止迭代。
3.6.5 ForceConv
值：实数。默认值：1e-8。指定几何优化的导数收敛阈值。导数矢量的长度小于此值并且满
足其他收敛条件，将停止迭代。
Page | 13
4
参考样例
程序附带了 NH2 体系和 HClLi 体系的混合(R2S)插值、三维样条插值、多体展开的三维样条插
值和三维再生核插值的输入文件样例。其中 NH2 体系是非网格点的情形。通过阅读这些样例的输
入文件，可以了解各种插值方法的参数指定细节。按照第一节所示的命令分别编译这些输入文
件，并检查结果，可以检查各种插值方法的适用性。
5
引用
程序作者：Huanchen Zhai
Email：stczhc@gmail.com
混合方法引用文献：Huanchen Zhai, and Shi Ying Lin. “A Fast Hybrid Method for Constructing
Multidimensional Potential Energy Surfaces From ab initio Calculations: A New Global Analytic PES of
NH2 System”, Chemical Physics, 2015, 455 (7), pp 57-64.
再生核方法引用文献：T.S. Ho, H. Rabitz, J. Chem. Phys. 104 (1996) 2584.
样条方法引用文献：W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, NUMERICAL
RECIPES IN FORTRAN 77, Cambridge University Press, Cambridge, 1992.
几何优化引用文献：(a) Jun-Qiang Sun and Klaus Ruedenberg, J. Chem. Phys. 99 (1993) 5257. (b)
Jun-Qiang Sun, Klaus Ruedenberg and Gregory J. Atchity, J. Chem. Phys. 99 (1993) 5276. (c) Jun-Qiang
Sun and Klaus Ruedenberg, J. Chem. Phys. 101 (1994) 2157.
