24小时热门版块排行榜

返回列表

【奖励】本帖被评价9次，作者zzgyb增加金币 8.5 个

当前主题已经存档。

zzgyb

荣誉版主 (文坛精英)

[资源] 产生碳纳米管坐标的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 if(n temp=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             ichk=ichk+1;
            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 %  If it is outside the unit cell, translate it back using periodic
%  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   x2 =rs*cos(i*q4+q5);
  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 ]

回复此楼

» 收录本帖的淘帖专辑推荐

资源收集	第一性原理计算辅助工具

» 猜你喜欢

拟解决的关键科学问题还要不要写已经有8人回复
26申博已经有3人回复
存款400万可以在学校里躺平吗已经有22人回复
最失望的一年已经有4人回复
国自然申请面上模板最新2026版出了吗？已经有19人回复
请教限项目规定已经有3人回复
基金委咋了？2026年的指南还没有出来？已经有10人回复
基金申报已经有6人回复
推荐一本书已经有13人回复
疑惑？已经有5人回复

» 本主题相关商家推荐: (我也要在这里推广)

1楼 2007-10-16 12:22:39

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

sui2066

木虫 (职业作家)

应助: 5 (幼儿园)
金币: 6188.4
帖子: 3600
在线: 809小时
虫号: 382693

★★★★★ 五星级,优秀推荐

wo lai zuo sofa

回复此楼

2楼2007-10-16 13:14:29

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

yoyojn

木虫 (著名写手)

应助: 3 (幼儿园)
金币: 3918.4
帖子: 2550
在线: 425.2小时
虫号: 313258

运行结果

??? function [X TransVec NumAtom Diameter ChiralAngle]=createCNT(n,m)
|
Error: Function definitions are not permitted at the prompt or in scripts.