【转帖】径向分布函数程序与简单说明 (小木虫)
径向分布函数g(r)代表了球壳内的平均数密度
为离中心分子距离为r,体积为 的球壳内的瞬时分子数。
具体参见李如生,《平衡和非平衡统计力学》科学出版社:1995
CODE:
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
PARAMETER(NM=40000,PI=3.141592653589793D0,NHIS=100)
COMMON/LCS/X0(3,-2:2*NM),X(3,-2:2*NM,5),XIN(3,-2:2*NM),
$XX0(3,-2:2*NM),XX(3,-2:2*NM,5),XXIN(3,-2:2*NM)
COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM,NC,NN,MC
COMMON/WALLS/HI(3,3),G(3,3),DH,AREA,VOLUME,SCM(3)
COMMON/PBCS/HALF,PBCX,PBCY,PBCZ
COMMON/GR_VAR/ NGR
DIMENSION H(3,3),GG(0:NHIS),R(0:NHIS)
EQUIVALENCE(X0(1,-2),H(1,1))
C *****************************************************************
C 如何确定分子数密度:DEN_IDEAL
C 取分子总数作为模拟盒中的数密度,可保证采样分子总数=总分子数?
C====================================================================
C N1=MOLSP+1
C N2=MOLSP+NC
DEN_IDEAL=MOLSP
G11=G(1,1)
G22=G(2,2)
G33=G(3,3)
G12D=G(1,2)+G(2,1)
G13D=G(1,3)+G(3,1)
G23D=G(2,3)+G(3,2)
IF(NSWITCH.EQ.0)THEN
NGR=0
DELR=HALF/NHIS
DO I=1,NHIS
GG(I)=0.D0
R(I)=0.D0
ENDDO
ELSE IF(NSWITCH.EQ.1)THEN
NGR=NGR+1
DO I=1,MOLSP-1
DO J=I+1,MOLSP
C====================================================================
C USE PBC IN X DIRECTION: SUITABLE FOR PBCX=1
C NOT GREAT PROBLEM FOR PBCX=0
C (THIS TIME USUALLY |DELTA X| < HALF)
C====================================================================
XIJ=X0(1,I)-X0(1,J)
IF(XIJ.GT.+HALF)XIJ=XIJ-PBCX
IF(XIJ.LT.-HALF)XIJ=XIJ+PBCX
YIJ=X0(2,I)-X0(2,J)
IF(YIJ.GT.+HALF)YIJ=YIJ-PBCY
IF(YIJ.LT.-HALF)YIJ=YIJ+PBCY
ZIJ=X0(3,I)-X0(3,J)
IF(ZIJ.GT.+HALF)ZIJ=ZIJ-PBCZ
IF(ZIJ.LT.-HALF)ZIJ=ZIJ+PBCZ
RSQ=XIJ*(G11*XIJ+G12D*YIJ+G13D*ZIJ)+
$ YIJ*(G22*YIJ+G23D*ZIJ)+G33*ZIJ*ZIJ
RRR=SQRT(RSQ)
RRR=RRR/H(1,1)
C====================================================================
C 以上用数组G和H的结果与下同
C RRR=SQRT(XIJ**2+YIJ**2+ZIJ**2)
C G11=H(1,1)**2
C====================================================================
IF(RRR.LT.HALF)THEN
IG=INT(RRR/DELR)
GG(IG)=GG(IG)+2
ENDIF
ENDDO
ENDDO
ELSE IF(NSWITCH.EQ.2)THEN
DO I=1,NHIS
R(I)=DELR*(I+0.5D0)
ENDDO
DO I=1,NHIS
VB=(4.D0/3.D0)*PI*(((I+1)**3-I**3)*(DELR**3))
GNID=VB*DEN_IDEAL
GG(I)=GG(I)/(NGR*MOLSP*GNID)
ENDDO
OPEN(UNIT=31,FILE="GR.DAT")
DO I=1,NHIS
WRITE(31,*)R(I),GG(I)
ENDDO
CLOSE(31)
ENDIF
RETURN
END
这样的代码看着不够明了。。。。。。
伪代码:
for (int i = 0; i < TOTN - 1; ++i)
for (int j = i + 1; j < TOTN; ++j) {
double dij = sqrt( pow(Pos[0]-Pos[j][0], 2) + pow(Pos[1]-Pos[j][1], 2) + pow(Pos[2]-Pos[j][2], 2));
int kbin = func(dij); // dij所对应的bin的序号
g(kbin) += 2;
}
// normalize
for (int k = 0; k < NBIN; ++k)
g(k) /= 4.0 * PI * r(k) * r(k) * dr * RHO; // r 为第k个bin所对应的距离值
calculate radial distribution function in molecular dynamics (转载科学网樊哲勇)
Here are the computer codes for this article:
md_rdf.cpp
find_rdf.m
test_rdf.m
Calculating radial distribution function in molecular dynamics
First I recommend a very good book on molecular dynamics (MD) simulation: the book entitled "Molecular dynamics simulation: Elementary methods" by J. M. Haile. I read this book 7 years ago when I started to learn MD simulation, and recently I enjoyed a second reading of this fantastic book. If a beginner askes me which book he/she should read about MD, I will only recommend this. This is THE BEST introductory book on MD. It tells you what is model, what is simulation, what is MD simulation, and what is the correct attitude for doing MD simulations.
In my last blog article, I have presented a Matlab code for calculating velocity autocorrelation function (VACF) and phonon density of states (PDOS) from saved velocity data. In this article, I will present a Matlab code for calculating the radial distribution function (RDF) from saved position data. The relevant definition and algorithm are nicely presented in Section 6.4 and Appendix A of Haile's book. Here I only present a C code for doing MD simulation and a Matlab code for calculating and plotting the RDF. We aim to reproduce Fig. 6.22 in Haile's book!
Step 1.
Use the C code provided above to do an MD simulation. Note that I have used a different unit systems than that used in Haile's book (he used the LJ unit system). This code only takes 14 seconds to run in my desktop. Here are my position data (there are 100 frames and each frame has 256 atoms):
r.txt
Step 2.
Write a Matlab function which can calculate the RDF for one frame of positions:
function [g] = find_rdf(r, L, pbc, Ng, rc)
% determine some parameters
N = size(r, 1); % number of particles
L_times_pbc = L .* pbc; % deal with boundary conditions
rho = N / prod(L); % global particle density
dr = rc / Ng; % bin size
% accumulate
g = zeros(Ng, 1);
for n1 = 1 : (N - 1) % sum over the atoms
for n2 = (n1 + 1) : N % skipping half of the pairs
r12 = r(n2, :) - r(n1, :); % position difference vector
r12 = r12 - round(r12 ./L ) .* L_times_pbc; % minimum image convention
d12 = sqrt(sum(r12 .* r12)); % distance
if d12 < rc % there is a cutoff
index = ceil(d12 / dr); % bin index
g(index) = g(index) + 1; % accumulate
end
end
end
% normalize
for n = 1 : Ng
g(n) = g(n) / N * 2; % 2 because half of the pairs have been skipped
dV = 4 * pi * (dr * n)^2 * dr; % volume of a spherical shell
g(n) = g(n) / dV; % now g is the local density
g(n) = g(n) / rho; % now g is the RDF
end
Step 3.
Write a Matlab script to load the position data, call the function above, and plot the results:
clear; close all;
load r.txt; % length in units of Angstrom
% parameters from MD simulation
N = 256; % number of particles
L = 5.60 * [4, 4, 4]; % box size
pbc = [1, 1, 1]; % boundary conditions
% number of bins (number of data points in the figure below)
Ng = 100;
% parameters determined automatically
rc = min(L) / 2; % the maximum radius
dr = rc / Ng; % bin size
Ns = size(r, 1) / N; % number of frames
% do the calculations
g = zeros(Ng, 1); % The RDF to be calculated
for n = 1 : Ns
r1 = r(((n - 1) * N + 1) : (n * N), :); % positions in one frame
g = g + find_rdf(r1, L, pbc, Ng, rc); % sum over frames
end
g = g / Ns; % time average in MD
% plot the data
r = (1 : Ng) * dr / 3.405;
figure;
plot(r, g, 'o-');
xlim([0, 3.5]);
ylim([0, 3.5]);
xlabel('r^{\ast}', 'fontsize', 15)
ylabel('g(r)', 'fontsize', 15)
set(gca, 'fontsize', 15);
Here is the figure I obtained:
Does it resemble Fig. 6. 22 in Haile's book?
转载于:https://www.cnblogs.com/Simulation-Campus/p/8994935.html
【转帖】径向分布函数程序与简单说明 (小木虫)相关推荐
- matlab迭代算法实例sor,SOR迭代 - 程序语言 - MATLAB/Mathematica - 小木虫论坛-学术科研互动平台...
方法一:建立了SOR.m的脚本文件,实现的是SOR迭代,程序语言如下: %SOR迭代 clear; clc; format long; i=1; n=6; H=hilb(n); X=ones(n,1) ...
- matlab寻峰代码,寻峰的函数!! - 程序语言 - MATLAB/Mathematica - 小木虫论坛-学术科研互动平台...
我这里的数据是pgm的,我将其处理成多个高斯拟合的形式,现阶段只能将其最大的那个拟合出来,其他的高斯拟合我需要找到其峰值的位置! 我把前边的语句先列举上: function [ OutArr ] = ...
- matlab模糊控制m函数,模糊控制m文件运行出错 - 程序语言 - MATLAB/Mathematica - 小木虫论坛-学术科研互动平台...
Error using parsrule (line 182) Output MF index is too high Error in readfis (line 231) out=parsrule ...
- experience(1):一个博士的经历(小木虫精华帖)
今天看到一篇一个博士生写的关于他读博士的整个过程,包含了读博的目的.在读博途中的心理变化,大论文写作建议以及对一些邮件的回复建议,其中包含大多数人会遇到的南入门和选题方面的答复还是很有借鉴意义的. [ ...
- java窗口小程序atm_简单的小程序实现ATM机操作
简单的小程序实现ATM机操作 代码如下: package Day06; import java.util.Scanner; public class TestAccount { public stat ...
- python有趣小程序代码,简单的小程序代码
谁能用python帮我写一个小程序,让用户输入任意9个数字,然后输出排序后的结果. 我只写一个函数:>>> def littleFunc(): data =[] #初始化列表 for ...
- matlab简单的程序,一段简单的matlab程序 - 程序语言 - 小木虫 - 学术 科研 互动社区...
原程序是可以运行的,为全面理解程序内容,我将分以下几个部分进行分析: 1."for x=varx"怎么理解? 请参看Matlab关于for函数的帮助文件: Syntax:for i ...
- c语言写的简单小程序,一些简单的小程序_6——C语言篇
1.写一个函数实现任意行列数的乘法表 #define _CRT_SECURE_NO_WARNINGS #include void mul(int n) { int i = 0; int j = 0; ...
- linux的命令参考手册,Linux常用命令汇总——可当作简要参考手册 - 程序语言 - 小木虫 - 学术 科研 互动社区...
基础命令 系统分区 #磁盘由盘片.机械手臂.磁头和主轴马达组成,数据写入是在盘片上面.盘片分为扇面.柱面与扇区,扇区只有512bytes大小.磁盘第一个扇区记录了"主引导分区"(可 ...
最新文章
- 第9章 数据字典(选项)管理
- java 静态内部类 弱引用_Java基础 强引用、弱引用、软引用、虚引用
- 【Clickhouse】问题记录
- 成功解决C4996: ‘fopen‘: This function or variable may be unsafe. Consider using fopen_s instead
- 【转】深入浅出理解有限状态机
- Js中Symbol对象
- python实现windows Service服务程序
- 【接口自动化测试】使用Fitness实现接口自动化测试
- 【linux】常用命令之scp命令
- 在 Mac OS X 安装gcc编辑环境,make不能用时参考
- Java 认证考试 OCAJP 经验总结
- 使用VSCode创建Java项目
- pyautogui的两天坑moveto图像识别
- python数学符号读法大全_常用数学符号读法大全
- 移动端开发的未来在哪里?-Tamic博客
- DGV:人类基因组结构变异数据库
- java中GUI中显示当前时间_javaGUI界面实现动态时间显示——Swing中的计时器Timer
- 桌面计算机 回收站图标,如何在计算机桌面上还原回收站图标?
- 一位女生程序员微自传
- 合作共赢!荷兰Swissflow成功入驻ISweek工采网!