| 查看: 1934 | 回复: 12 | ||||
| 【奖励】 本帖被评价9次，作者zzgyb增加金币 8.5 个 | ||||
| 当前主题已经存档。 | ||||
[资源]
产生碳纳米管坐标的Matlab程序
|
||||
|
毛主席说了，回帖给评价的同志是好同志！ 以下是一个产生任意(m,n)碳管坐标的一个Matlab程序。只需把下面的程序copy，存成creatCNT.m，即可用Matlab直接运行。比如在Matlab提示符下输入 >>creatCNT(8,0) 则产生(8,0)碳纳米管的坐标，输出文件为nanotube.xyz和nanotube.pdb格式。可以用GaussView等其他的可视化软件打开。 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [X TransVec NumAtom Diameter ChiralAngle]=createCNT(n,m) %This function creates the coordinates of a (n,m) nanotube, with n>m %Usage: % [X TransVec NumAtom Diameter ChiralAngle]=createCNT(8,0); % X -- coordinates of atoms in the nanotube unit cell % TransVec -- Translational Vector T of the nanotube % NumAtom -- Number of atoms in one unit cell of nanotube % Diamter -- diameter of the nanotube % ChiralAngle -- chiral angle of the nanotube % by bshan 2005 %change order if n n=m; m=temp; end nk=6000; acc=1.44; %acc is the C-C bond length % a: the lengh of the unit vector % pi =3.1415926 % sq3=1.732 sq3=sqrt(3.0d0); a=sq3*acc; l2 = n*n+m*m+n*m; l = sqrt(l2); dt = a*l/pi; nd = gcd(n,m); %nd is the greatest common dividor of n,m if(mod((n-m),3*nd)==0) ndr = 3*nd; else ndr = nd; end nr = (2*m+n)/ndr; %first component of T ns = -(2*n+m)/ndr; %second component of T nt2 = 3*l2/ndr/ndr; nt = floor(sqrt(nt2)); nn = 2*l2/ndr; %nn is number of haxagonals per unit cell %start searching the symmetry vector.The search range is set from %-abs(n60) to abs(n60) ichk=0; if(nr==0) n60=1; else n60=nr*4; end for i=-abs(n60):abs(n60) for j=-abs(n60):abs(n60) j2 = nr*j-ns*i; if(j2==1) j1= m*i-n*j; if((j1>0)&&(j1 nnp(ichk)=i; nnq(ichk)=j; end end end end %symmetry vector is not found. Probably the search range is not large %enough. if(ichk==0) error('not found p,q strange!!') end %more than one symmetry vector found. if(ichk>=2) error('more than 1 pair p,q strange!!') end %symmetry vector (np,nq) np = nnp(1); nq = nnq(1); %msg=sprintf('the symmetry vector is %d %d',np,nq); %disp(msg); % % r:|R| , c:|C_h|, t:|T| % lp = np*np + nq*nq + np*nq; r=a*sqrt(lp); c=a*l; t=sq3*c/ndr; % nn: the number of hexagon in the unit cell N % rs: radius of the tube if((2*nn)>nk) error('parameter nk is too small!') end rs=c/(2.0d0*pi); %msg=sprintf('radius=%f ,t=%f',rs,t); %disp(msg); % q1: the chiral angle for C_h % q2: the chiral angle for R % q3: the chiral between C_h and R q1=atan((sq3*m)/(2*n+m)); q2=atan((sq3*nq)/(2*np+nq)); q3=q1-q2; % q4: a period of an angle for the A atom % q5: the difference of the angle between the A and B atom q4=2.0*pi/nn; q5=acc*cos((pi/6.0d0)-q1)/c*2.0*pi; % h1: % h2: Delta z between the A and B atom h1=abs(t)/abs(sin(q3)); h2=acc*sin((pi/6.0)-q1); % The A atom ii=0; for i=0:nn-1 x1=0; y1=0; z1=0; k=floor(i*abs(r)/h1); x1=rs*cos(i*q4); y1=rs*sin(i*q4); z1=(i*abs(r)-k*h1)*sin(q3); kk2=abs(floor((z1+0.0001)/t)); % Check the A atom is in the unit cell 0<=z1 % boundary concition if(z1>=t-0.0001) z1=z1-t*kk2; elseif(z1<0) z1=z1+t*kk2; else end ii=ii+1; x(ii)=x1; y(ii)=y1; z(ii)=z1; % The B atom z3=(i*abs(r)-k*h1)*sin(q3)-h2; ii=ii+1; % Check the B atom is in the unit cell 0<=z3 if((z3>=0)&&(z3 y2 =rs*sin(i*q4+q5); z2 =(i*abs(r)-k*h1)*sin(q3)-h2; x(ii)=x2; y(ii)=y2; z(ii)=z2; else x2 =rs*cos(i*q4+q5); y2 =rs*sin(i*q4+q5); z2 =(i*abs(r)-(k+1)*h1)*sin(q3)-h2; kk=abs(floor(z2/t)); if(z2>=t-0.0001) z2 =z2-t*kk; elseif(z2<0) z2 =z2+t*kk; else end x(ii)=x2; y(ii)=y2; z(ii)=z2; end end %total number of atoms ntotal=2*nn; for i=1:ntotal X(i,: )=[x(i) y(i) z(i)]; end TransVec=t; NumAtom=ntotal; Diameter=rs*2; ChiralAngle=atan((sq3*n)/(2*m+n))/pi*180; %write coordinates to xyz file format fid=fopen('nanotube.xyz','w'); fprintf(fid,'%d\n',size(X,1)); fprintf(fid,'created by GUI_TB program\n'); for i=1:ntotal fprintf(fid,'C %f %f %f \n',X(i,1),X(i,2),X(i,3)); end fclose(fid); %write coordinates to pdb file format(for visualization using rasmol) fid=fopen('nanotube.pdb','w'); for i=1:ntotal fprintf(fid,'ATOM %6d C ADE A 1 %8.3f%8.3f%8.3f 1. 0.\n',i,X(i,1),X(i,2),X(i,3)); end fprintf(fid,'TER\n'); fclose(fid); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [ Last edited by csfn on 2008-12-29 at 20:11 ] |
回复此楼» 收录本帖的淘帖专辑推荐
 资源收集 | 第一性原理计算辅助工具 |
» 猜你喜欢
0703化学
已经有26人回复
-
274求调剂求调剂
已经有11人回复
-
327求调剂
已经有10人回复
-
调剂
已经有19人回复
-
275 求调剂
已经有7人回复
-
求材料调剂，一志愿郑州大学289分
已经有21人回复
-
电子信息346
已经有4人回复
-
求助071001调剂！！！
已经有7人回复
-
一志愿深大085601材料工程专业（专硕）300分可以调剂去哪
已经有12人回复
-
一志愿华东理工085601材料工程303分求调剂
已经有11人回复
» 本主题相关商家推荐: (我也要在这里推广)
2楼2007-10-16 13:14:29
3楼2007-12-04 11:06:23
4楼2007-12-25 02:24:20
5楼2007-12-25 02:24:53
6楼2007-12-25 02:30:31
7楼2007-12-25 09:26:57
8楼2007-12-27 21:04:43
9楼2008-03-03 02:53:27
