关于"Understanding Molecular Simulation: From Algorithms to Applications"分子动力学代码的问题



你需要读过这本书才能回答我的问题。

-背景- 我正在尝试学习分子动力学,我从应该是介绍性的书籍开始,但我在理解方程或其解释方面遇到了很多麻烦。我尝试了大约 6 或 7 本不同的书,从介绍章节开始,然后浏览这本书,寻找我可以复制的代码,为我提供一些方程式的背景,因为我无法理解它们。很明显,他们希望你有很多先决条件的物理知识(似乎是特定的MD知识)来理解方程如何转化为一段代码。我上过物理课程和3D物理引擎课程,但这些知识对这项工作没有帮助。 我唯一能找到的代码是在"理解分子模拟:从算法到应用">中,我有一些问题,我在书中找不到答案。如果您有任何更适合初学者的先修书籍或课程来帮助我为本书中的材料提供上下文,那就太好了 -背景结束-

这正是代码在书中的显示方式,我读到你应该在主程序结束后放置子例程,所以我就是这样做的。因为我只能找到处理非常小段代码的 fortran 教程,所以我不确定我是否应该在这里显示的所有地方都有stopendreturnimplicit none是否在它应该在的地方,以及我是否应该将每个子例程放在不同的文件中或类似的东西中。 循环与python中的循环基本相同,因此我已经能够理解一些代码,但我还没有看到更大的图景。

我面临的主要问题是,我使用的一堆变量没有在书中的任何地方声明(尽管我找到了关于它们代表其中一些变量的解释),这可能就是程序无法运行的原因。我收到此错误消息"启动失败。找不到二进制文件"。

如何声明变量delttmaxf(力)、en(能量)、ranflattice_posxvtempxmdtxrboxr2ir6iffxxvi?其中一些甚至可能不是变量,但我认为它们没有内置到 fortran 中,因为我找不到有关它们的任何信息。另外,我猜tmax将与npart(N个粒子)的数字相同,那么为什么它们会使用2个不同的变量呢?

这些变量意味着什么?(除了我能找到的fen)

每个原子的x,y,z坐标存储在哪里?我假设它们将存储在 1 个列表(或 fortran 中的任何等效列表中)中,每个列表存储 3 个值(对于 x,y,z)或 3 个列表,这些列表在每个列表中的所有 x,y,z 坐标,在init子例程中。然而,情况并非如此,因为只有xv是循环的。

en在子例程开始时立即设置为 0 时,为什么它是force()的输入?您不能只在函数中声明它而不将其作为输入吗?

program main
call init
t=0
do while (t.lt.tmax)
call force(f,en)
call integrate(f,en)
t=t+delt
call sample
enddo
stop
end
subroutine init
sumv=0
sumv2=0
do i=1, npart
x(i)=lattice_pos(i)
v(i)=(ranf()-0.5)
sumv=sumv+v(1)
sumv2=sumv2+v(i)**2
enddo
sumv=sumv/npart
sumv2=sumv2/npart
fs=sqrt(3*temp/sumv2)
do i=1,npart
v(i)=(v(i)-sumv)*fs
xm(i)=x(i)-v(i)*dt
enddo
return
subroutine force(f, en)
en=0
do i=1, npart-1
do j=i+1, npart
xr=xr-box*nint(xr/box)
r2=xr**2
if (r2.lt.rc2) then
r2i=1/r2
r6i=r2i**3
ff=48*r2i*r6i*(r6i-0.5)
f(i)=f(i)+ff*xr
f(j)=f(j)-ff*xr
en=en+4*r6i*(r6i-1)-ecut
endif
enddo
enddo
return
end
subroutine integrate(f,en)
sumv=0
sumv2=0
dp i=1,npart
xx=2*x(i)-xm(i)+delt**2*f(i)
vi=(xx-xm(i))/(2*delt)
sumv=sumv+vi
sumv2=sumv2+vi**2
xm(i)=x(i)
x(i)=xx
enddo
temp=sumv2/(3*npart)
etot=(en+0.5*sumv2)/npart
return
end
implicit none
end program main

与许多其他教科书一样,"理解分子模拟:从算法到应用"中的示例代码是伪代码。它不完整(例如,缺少声明),旨在比"仅数学"更明确地说明算法。在这本书中,这是令人困惑的,因为他们为伪代码编写了不完整的Fortran。

这仍然是一个很好的"一般问题":了解如何使用 Fortran 的老式教科书材料:-)

在 http://www.acmm.nl/molsim/frenkel_smit/index.html(阿姆斯特丹多尺度建模中心的网站),你可以找到这本书的实际Fortran代码 http://www.acmm.nl/molsim/frenkel_smit/README.html 正如伊恩·布什(Ian Bush)所写,代码是老式的,Allen&Tildesley的书有一个更新更新的代码:https://github.com/Allen-Tildesley/

至于你的问题:

  1. 你如何声明变量delt,tmax,f(force),en(energy),ranf,lattice_pos,x,v,temp,xm,dt,xr,box,r2i,r6i,ff,xx和vi?其中一些甚至可能不是变量,但我认为它们没有内置到 fortran 中,因为我找不到有关它们的任何信息。另外,我猜tmax将与npart(N个粒子)的数字相同,那么为什么它们会使用2个不同的变量呢?

    您必须阅读有关算法的章节才能理解变量。然而,它们的列表可以在本书的"xxi"页上找到。

  2. 每个原子的x,y,z坐标存储在哪里?我假设它们将存储在 1 个列表(或 fortran 中的任何等效列表中)中,每个列表存储 3 个值(对于 x,y,z)或 3 个列表,每个列表中的所有 x,y,z 坐标,在 init 子例程中。然而,情况并非如此,因为只有 x 和 v 是循环的。

    同样,伪代码有点限制。您有几个选项来声明变量:(i)作者(请参阅上面的文件链接)使用"公共块",这是一种特定于Fortran的全局变量类型。(ii) Fortran 90 和更新版本中可用的模块变量是"更好的全局变量"。(iii) 将所有数据封装在派生类型中,您可以在主程序中声明这些派生类型并传递给子例程。

    位置存储(在示例中)为数组x(npmax)y(npmax)z(npmax)npmax是粒子的最大数量。这不会显示在伪代码中。请注意,在较新版本的 Fortran(90 及更高版本)中,您可以编写向量运算,例如x = x + v*dt(其中 x 和 v 是所有粒子的位置和速度,dt 是标量时间步长)。

  3. 为什么 en 是 force() 的输入,当它在子例程开始时立即设置为 0 时?您不能只在函数中声明它而不将其作为输入吗?

    在 Fortran 中,将需要更新的变量作为参数传递给子例程是很常见的,但不是强制性的。编程风格各不相同。

最新更新