Python与fortran混合使用-f2py
前言
选择python做接口解析，数值计算部分采用fortran来实现的时候， 我们可以利用python/c的接口来实现，不过这需要对f90再做个wrap，这样有点麻烦，更好的方式是通过f2py来实现。
我们的目标是：
- 掌握f2py的使用
- 复习一下简单的分子动力学内容
我们用分子动力学里的经典LJ对势，写个F90函数计算原子体系的能量和力，用f2py生成库，然后在python中调用这个库函数 原子体系的总能等于原子对势之和
简单复习一下分子动力学
原子体系的总能等于原子对势之和
\begin{equation} U = \sum_{i \lt j} u(r_{ij}) \end{equation}
\begin{equation} u(r_{ij})=4 \epsilon [(\frac{r_{ij}}{\sigma})^{-12} - (\frac{r_{ij}}{\sigma})^{-6}] \end{equation}
那么原子在x方向上的力就是：
\begin{aligned}
F_{i, x} &=-\frac{\partial U}{\partial x_{i}} \\\
&=-\frac{\partial}{\partial x_{i}} \sum_{j \neq i} u\left(r_{i j}\right) \\\
&=-\sum_{j \neq i} \frac{\partial r_{i j}}{\partial x_{i}} \frac{\partial}{\partial r_{i j}} u_{12}\left(r_{i j}\right)
\end{aligned}
由\(r_{ij}^2=(x_i - x_j)^2 + (y_i - y_j)^2 + (z_i+z_j)^2\) ， 我们有
\begin{equation} F_{i,x}= - \sum_{(j \neq i)} (\frac{x_i-x_j}{r_{ij}})\frac{\partial}{\partial r_{ij}} u(r_{ij}) \end{equation}
自然地我们得到3维情况下的力的公式，同时引入向量符号\(r_{ij} = r_j - r_i\),
\begin{aligned}
\mathbf{F}_{i} &=\sum_{j \neq i}\left(\frac{\mathbf{r}_{i j}}{r_{i j}}\right) \frac{\partial}{\partial r_{i j}} u\left(r_{i j}\right) \\\
&=\sum_{j \neq i}\left(\frac{\mathbf{r}_{i j}}{r_{i j}}\right)\left(\frac{4 \epsilon}{\sigma}\right)\left[-12\left(\frac{r_{i j}}{\sigma}\right)^{-13}+6\left(\frac{r_{i j}}{\sigma}\right)^{-7}\right] \\\
&=\sum_{j \neq i} \mathbf{r}_{i j}\left(\frac{\epsilon}{\sigma}\right)\left[-48\left(\frac{r_{i j}}{\sigma}\right)^{-14}+24\left(\frac{r_{i j}}{\sigma}\right)^{-8}\right]
\end{aligned}
考虑周期性边界条件：每个原子只考虑临近的周围镜像的原子，这个通过找2个粒子之间的最近镜像距离来确定是哪个image
\begin{equation} r_{ij}^0 = r_{ij} - L nint(r_{ij}/L) \end{equation}
L是盒子的长度，对于非正方盒子来讲 L是向量形式。
最后我们在对势函数中引入一个截断距离\(r_c\), 超过这个距离时势能为0。通常$r_c=2.5σ$或者更大值。
\begin{equation}
u(r_{i j})=\{\begin{array}{c}
4 \epsilon[(\frac{r_{i j}}{\sigma})^{-12}-(\frac{r_{i j}}{\sigma})^{-6}]-4 \epsilon[(\frac{r_{c}}{\sigma})^{-12}-(\frac{r_{c}}{\sigma})^{-6}] \quad r_{i j} \leq r_{c} \\\
0 & r_{i j}>r_{c}
\end{array}
\end{equation}
为了简化，我们采用无量纲单位来表示， 这样对势方程组的上面部分的形式就变为：
\begin{equation} u(r_{i j})=4[r_{i j}^{-12}-r_{ij}^{-6}]-4[r_{c}^{-12}-r_{c}^{-6}] \end{equation}
最终的fortran部分的实现代码，保存为ljlib.F90文件：
subroutine EnergyForces(Pos, L, rc, PEnergy, Forces, Dim, NAtom)
implicit none
integer, intent(in) :: Dim, NAtom
real(8), intent(in), dimension(0:NAtom-1, 0:Dim-1) :: Pos
real(8), intent(in) :: L, rc
real(8), intent(out) :: PEnergy
real(8), intent(inout), dimension(0:NAtom-1, 0:Dim-1) :: Forces
!f2py intent(in, out) :: Forces
real(8), dimension(Dim) :: rij, Fij, Posi
real(8) :: d2, id2, id6, id12
real(8) :: rc2, Shift
integer :: i, j
PEnergy = 0.
Forces = 0.
Shift = -4. * (rc**(-12) - rc**(-6))
rc2 = rc * rc
do i = 0, NAtom - 1
!store Pos(i,:) in a temporary array for faster access in j loop
Posi = Pos(i,:)
do j = i + 1, NAtom - 1
rij = Pos(j,:) - Posi
rij = rij - L * dnint(rij / L)
!compute only the squared distance and compare to squared cut
d2 = sum(rij * rij)
if (d2 > rc2) then
cycle
end if
id2 = 1. / d2 !inverse squared distance
id6 = id2 * id2 * id2 !inverse sixth distance
id12 = id6 * id6 !inverse twelvth distance
PEnergy = PEnergy + 4. * (id12 - id6) + Shift
Fij = rij * ((-48. * id12 + 24. * id6) * id2)
Forces(i,:) = Forces(i,:) + Fij
Forces(j,:) = Forces(j,:) - Fij
enddo
enddo
end subroutine EnergyForces
subroutine VVIntegrate(Pos, Vel, Accel, L, CutSq, dt, KEnergy, PEnergy, Dim, NAtom)
implicit none
integer, intent(in) :: Dim, NAtom
real(8), intent(in) :: L, CutSq, dt
real(8), intent(inout), dimension(0:NAtom-1, 0:Dim-1) :: Pos, Vel, Accel
!f2py intent(in, out) :: Pos, Vel, Accel
real(8), intent(out) :: KEnergy, PEnergy
external :: EnergyForces
Pos = Pos + dt * Vel + 0.5 * dt*dt * Accel
Vel = Vel + 0.5 * dt * Accel
call EnergyForces(Pos, L, CutSq, PEnergy, Accel, Dim, NAtom)
Vel = Vel + 0.5 * dt * Accel
KEnergy = 0.5 * sum(Vel*Vel)
end subroutine VVIntegrate
进入编译调试
gfortran -c liblj.F90
看看代码有无报错，恩，应该是没有的，让我们继续
准备f2py前需要的知识点
-
如果fortran代码对所有的变量都有严格的类型定义，我们很少需要使用f2py对变量进行转换， 不过对于带有intent(inout)属性的数组变量，需要添加!f2py 来告诉f2py我们想处理的数据类型;写成形式
!f2py intent(in,out) :: var -
另外必须为所有输出参数明确指定尺寸！换句话说，输出自变量不能具有假定的大小或可分配的大小
integer, allocatable, intent(out) :: array1(:) ! Not valid integer, intent(out) :: array2(:) ! Not valid integer, intent(out) :: array3(10) ! Valid integer, intent(in) :: array4(:) ! Valid -
不支持派生类型
开始运行f2py
根据上面的知识点，这里对ljlib.F90里的EnergyForces函数，添加一行
real(8), intent(inout), dimension(0:NAtom-1, 0:Dim-1) :: Forces
!f2py intent(in,out) :: Forces
对VVIntegral函数，添加一行!f2py
real(8), intent(inout), dimension(0:NAtom-1, 0:Dim-1) :: Pos, Vel, Accel
!f2py intent(in,out) :: Pos, Vel, Accel
命令行下运行
python -m numpy.f2py -c ljlib.F90 -m ljlib
在目录下会生成ljlib.cpythonxxx.so 的动态链接文件。
打开python解释器，导入这个ljlib模块
看起来python导入ljlib库完成了，f2py用起来还是很方便的吧
>>> import ljlib
>>> print(ljlib.__doc__)
This module 'ljlib' is auto-generated with f2py (version:2).
Functions:
penergy,forces = energyforces(pos,l,rc,forces,dim=shape(pos,1),natom=shape(pos,0))
pos,vel,accel,kenergy,penergy = vvintegrate(pos,vel,accel,l,cutsq,dt,dim=shape(pos,1),natom=shape(pos,0))
看看各个函数的接口
>>> print(ljlib.energyforces.__doc__)
penergy,forces = energyforces(pos,l,rc,forces,[dim,natom])
Wrapper for ``energyforces``.
Parameters
----------
pos : input rank-2 array('d') with bounds (natom,dim)
l : input float
rc : input float
forces : input rank-2 array('d') with bounds (natom,dim)
Other Parameters
----------------
dim : input int, optional
Default: shape(pos,1)
natom : input int, optional
Default: shape(pos,0)
Returns
-------
penergy : float
forces : rank-2 array('d') with bounds (natom,dim)
>>> print(ljlib.vvintegrate.__doc__)
pos,vel,accel,kenergy,penergy = vvintegrate(pos,vel,accel,l,cutsq,dt,[dim,natom])
Wrapper for ``vvintegrate``.
Parameters
----------
pos : input rank-2 array('d') with bounds (natom,dim)
vel : input rank-2 array('d') with bounds (natom,dim)
accel : input rank-2 array('d') with bounds (natom,dim)
l : input float
cutsq : input float
dt : input float
Other Parameters
----------------
dim : input int, optional
Default: shape(pos,1)
natom : input int, optional
Default: shape(pos,0)
Returns
-------
pos : rank-2 array('d') with bounds (natom,dim)
vel : rank-2 array('d') with bounds (natom,dim)
accel : rank-2 array('d') with bounds (natom,dim)
kenergy : float
penergy : float
看起来python导入ljlib库完成了，f2py用起来还是很方便的吧