转帖径向分布函数程序与简单说明 (小木虫) Posted 2020-11-02 Simyang tags: 篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了转帖径向分布函数程序与简单说明 (小木虫)相关的知识,希望对你有一定的参考价值。 径向分布函数g(r)代表了球壳内的平均数密度为离中心分子距离为r,体积为 的球壳内的瞬时分子数。具体参见李如生,《平衡和非平衡统计力学》科学出版社:1995 CODE: SUBROUTINE GR(NSWITCH) 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+1C 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,MOLSPC====================================================================C USE PBC IN X DIRECTION: SUITABLE FOR PBCX=1C 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)**2C==================================================================== 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? 以上是关于转帖径向分布函数程序与简单说明 (小木虫)的主要内容,如果未能解决你的问题,请参考以下文章 生物科研软件(从科研到发文章)(转自小木虫) neural computation这个期刊怎么样小木虫 关于towhee 回答zhang_jaj的问题-转载小木虫 学术 一个博士的经历(小木虫精华帖,留着细细体会!) [概率论与数理统计]笔记:3.3 随机向量的函数的分布与数学期望 逆分布函数法生成随机数(指数分布) R语言实现