第一种情况: 从程序自带的各种晶体及有机模型中导入体系的晶胞

1． 打开MS，由file>import>structures>metals\>pure-metals>Fe导入Fe的晶胞。

2． 由build>Surfaces>cleave Surfaces打开对话框.

    在对话框中输入要建立的晶面（hkl），选择position，其中depth控制晶面层数。

3． 进入build>Supercell，输入A 、B 、C的值，得到想要的超晶胞。

4． 到该步骤，我们已经建立了一个周期性的超晶胞。如果要做周期性计算，则应选择build>Crystals>build vaccum slab，其中真空层通常选择10埃以上。如果建立团簇模型则选择build>Symmetry>Non-periodic Structure，去掉模型的周期性，并跟据自己的实际需要删除部分原子，得到想要的团簇模型。

5． 在表面插入分子时通过菜单栏上的几个小图标添加即可。



第二种情况: 手动建模，优点是可控制晶格常数。

6． 首先从文献中查到晶体的晶格常数的实验值。

7． 打开build>Crystals>build crystals，可见到对话框。

     在对话框中选择空间群与点群，然后在Lattice Parameter中设置晶胞基矢的长度及夹角。

8． 然后打开build>Add atom，从对话框中输入坐标。这里只需输入几个有代表性的原子的坐标，不必全部输入。在坐标输入前首先在option页面中选择coordinate system，或者分数坐标或者卡迪尔坐标。

9． 以下步骤重复2-5步。

10． 需要注意的是，采取什么样的团簇并不是任意的。原因是很多模型构造出来后在优化过程中往往不收敛。要避免这个问题的办法是查阅文献，参考文献上模型进行选取，因为它们的模型通常是经过试验证实收敛的。

几点说明

1． 与高斯相比，dmol3能够计算的体系更大。如果要研究表面的吸附，而模拟表面的团簇模型又比较大，建议采用dmol3。如果计算的是局部化学反应，而体系也不是很大，则可以使用高斯。

2． 关于是否考虑周期性条件的问题

研究金属表面时，团簇计算方法在前些年由于计算量小曾经被广泛的应用过，直到现在也被很多人在使用着，主要被用来计算吸附和多个分子的共吸附等，即不考虑化学键的断裂。 近年来由于国际上计算能力的提升，人们开始考虑周期性条件，这点从JPCA，JPCB，PRL，PRB，JACS等杂志上刊出的文章里也可以看出，但是计算量要大很多。

需要注意的是，由于团簇计算方法没有考虑周期性，即在k空间里只计算了Γ点，采用该方法计算表面的化学键的断裂（即表面扩散问题等）时有可能受到质疑。

3．在研究表面时，通常把团簇固定，只优化吸附在表面的分子，这一点可以通过菜单栏上的Modify>Constraint实现。首先选定团簇中需要固定的原子，然后在下面的对话框中打勾。同时也可以在Measurement里固定部分键长和键角。

4. 关于计算参数设置

主要有几个参数需要注意

1  对于Electronic页面，需要注意的是Core treatment,对于过渡金属原子通常需要考虑相对论效应，因此一般不使用All Electron方法。其他几种方法任选。

    Basis set应为DNP，Setup下的Quality一般选fine。为了提高计算速度，一个较好的办法是先用粗糙的Basis set和Quality进行优化，然后再提高精度。

2  还有一个非常重要的选项是Electronic>More>SCF里的Use smearing。这个关键字有助于加快收敛，但是设的多大往往会产生错误的结果，它也相当于允许的误差范围。具体设置办法可参考help。

其他的关键字可酌情设置。
第三种情况，编程序输出原子坐标.XYZ文件建模。
给出一例，希望对大家有所启发！我的QQ：705357207。欢迎交流！

则产生(8,0)碳纳米管的坐标，输出文件为nanotube.xyz和nanotube.pdb格式。可以用其他的可视化软件打开。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X TransVec NumAtom Diameter ChiralAngle]=createCNT(n,m,celln)
%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

if nargin==2
celln=1;
end

%change order if n<m
if(n<m)
temp=n;
n=m;
m=temp;
end


nk=6000;
acc=1.44; %acc is the C-C bond length 1.42??
% 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<nn))
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<t
% 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<t

if((z3>=0)&&(z3<t))
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

if celln>1
xx=X;
for mnp=2:celln
X=[X;[xx(:,1:2),xx(:,3)+(mnp-1).*t]];
end
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 g_xq96 on 2010-4-22 at 22:21 ] 返回小木虫查看更多