Physics 版 (精华区)

发信人: zjliu (秋天的萝卜), 信区: Physics
标  题: [转载]用C写的CNT坐标生成程序
发信站: 哈工大紫丁香 (Thu Dec  4 16:07:01 2003), 站内信件


这是我暑假编写的，非螺旋的两天就搞定了，但螺旋型将近花了一个星期，甚至做了纸模

型。后来又进一步改进，可以按晶胞来生成。下面的程序是按照平面石墨折叠而形成的碳

管，没有做进一步优化——后来我又编了能调整C-C键长和一个键角的CNT坐标生成程序

?
这里就不贴了。

本程序生成的文件格式是Gaussian Input
File格式*.gjc，用Chem3D可以直接读出3D结构
。


//ccnt.cpp

//************************************************************
//*  Single-Walled Carbon Nanotubes  Coordinates Generator   *
//*                     Written by                           *
//*                        aqyw                              *
//*     Department of Chemistry,Beijing Normal University    *
//*                     August,2002                          *
//*                     QQ: 11071781                         *
//************************************************************

#include<stdio.h>
#include<conio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>

//*************************global
variables*********************************
int n1,m1;  //the relative position vector of repeated point of point
(n,m)
//**************************************************************************


//======to return the reduced form of two integers n,m=========
int minDivisor(int n,int m){
  n=fabs(n);
  m=fabs(m);
  int r,l;
  if(n>m){
    r=n%m;
    l=m;
  }
  else{
    r=m%n;
    l=n;
  }
  if(r==0)
    return l;
  return minDivisor(l,r);
}
reduce_n1m1(){
  if(m1==0){
    n1=1;
    return 0;
  }
  int md=minDivisor(n1,m1);
  n1/=md;
  m1/=md;
  return 0;
}
//=============================================================

void main(int paraN, char *paraStr[]){
  int n,m,cellN; //(n,m) and number of unit cells
  double bondL; //bond length(angstrom)
  FILE *fp;  //file ponter
  char fileNm[20]; //file name

  if(paraN<5){
    printf("Usage:   cocnt n m l(number of unit cells) destination\n");
    exit(0);
  }
  printf("\n--Cell-split SWCNT Atom Gaussian Input File (.GJC) Generator");
  printf("\n--(C)Written by Yang Wang, Sep.17th.\n\n");

  n=atoi(paraStr[1]);
  m=atoi(paraStr[2]);
  cellN=atoi(paraStr[3]);
  bondL=1.420;
  strcpy(fileNm,paraStr[4]);

  double a=bondL*sqrt(3),b=bondL;  //two vectors in a hexagon
  double phi=atan(sqrt(3)*m/(2*n+m));  //angle of the rolling vector
  double c=a*sqrt(n*n+m*m+n*m);  //perimeter of the tube
  double r=c/(2*M_PI); //radius of the tube
  char atmSbl='C'; //atom symbol
  int ttlAtmN; //totle atom number
  int ucAtmN;  //totle atom number in a unit cell

  //====================to determine a unit cell=====================
  n1=n+2*m;
  m1=n-m;
  reduce_n1m1();  //to determine n1,m1

  double ucL=a*sqrt(n1*n1+m1*m1+n1*m1);  //length of unit cell
  double pTtlAtmN=ucL/(3*sqrt(3)*b*b)*c*4*cellN;  //proximate total atom
numbe
r
  //======================================================================
  double atmCord[3][1900];  //array to store atom's 3D cartisian
coordinates
  double theta[1900]; //polar angle of each atom

  //first atom's coordinates(atom A1,in first spiral)
  atmCord[0][1]=0;  //z
  atmCord[1][1]=r;  //x
  atmCord[2][1]=0;  //y
  theta[1]=0;  //theta

  //second atom(atom B1,in first spiral)
  atmCord[0][2]=atmCord[0][1]-b*sin(M_PI/6-phi);
  theta[2]=b*cos(M_PI/6-phi)/r+theta[1];
  atmCord[1][2]=r*cos(theta[2]);
  atmCord[2][2]=r*sin(theta[2]);

  //third atom(atom A2,in first spiral)
  atmCord[0][3]=atmCord[0][2]+b*sin(M_PI/6+phi);
  theta[3]=b*cos(M_PI/6+phi)/r+theta[2];
  atmCord[1][3]=r*cos(theta[3]);
  atmCord[2][3]=r*sin(theta[3]);

  double Sz=atmCord[0][3]-atmCord[0][1];  //z-component of position vector
S b
etween A1 and A2
  double Stheta=theta[3]-theta[1];  //polar angle of an S vector

  double Wz=-a*sin(M_PI/3-phi)-atmCord[0][1];  //z-component of position
vecto
r W between A1 atoms in two adjacent spirals
  double Wtheta=a*cos(M_PI/3-phi)/r-theta[1];  //polar angle of an W
vector

  int atmInd=2;  //atom index number
  int i,j,k,l;  //loop index
  int j0=2;  //the atom index number for initiation
  int j1=(m==0)?2*n:1900;//nPerS*tbN;  //the maxium number of a whole
spiral c
ircle in a primitive cell
  atmCord[1][2]=r*cos(theta[2]);
  atmCord[2][2]=r*sin(theta[2]);

  //third atom(atom A2,in first spiral)
  atmCord[0][3]=atmCord[0][2]+b*sin(M_PI/6+phi);
  theta[3]=b*cos(M_PI/6+phi)/r+theta[2];
  atmCord[1][3]=r*cos(theta[3]);
  atmCord[2][3]=r*sin(theta[3]);

  double Sz=atmCord[0][3]-atmCord[0][1];  //z-component of position vector
S b
etween A1 and A2
  double Stheta=theta[3]-theta[1];  //polar angle of an S vector

  double Wz=-a*sin(M_PI/3-phi)-atmCord[0][1];  //z-component of position
vecto
r W between A1 atoms in two adjacent spirals
  double Wtheta=a*cos(M_PI/3-phi)/r-theta[1];  //polar angle of an W
vector

  int atmInd=2;  //atom index number
  int i,j,k,l;  //loop index
  int j0=2;  //the atom index number for initiation
  int j1=(m==0)?2*n:1900;//nPerS*tbN;  //the maxium number of a whole
spiral c
ircle in a primitive cell
  int k1=(m==0)?2:m;  //the number of spirals
  if(pTtlAtmN>1900){
    printf("Warning:There are too many atoms(%.0f).The maxium number of
atoms
is 1900!\n",pTtlAtmN);
    getch();
    printf("End of this task!");
    exit(0);
  }
  for(k=0;k<k1;k++){
    //=======for the second or later spirals=================
    if(k!=0){
      double da=(m==0)?0:k*sin(M_PI/3-phi)/sin(phi)+1;  //offset from
original
 A1 atom to the modified A1 atom which will

make the bottom of the tube more flat
      double fda=floor(da);
      da=(da==fda)?fda-1:fda;

      atmInd++;
      //the first atom A1 im current spiral
      atmCord[0][atmInd]=atmCord[0][1]+k*Wz;
      theta[atmInd]=theta[1]+k*Wtheta;
      atmCord[1][atmInd]=r*cos(theta[atmInd]);
      atmCord[2][atmInd]=r*sin(theta[atmInd]);

      //to determine the nearest atom above the A1 atom in former spiral
      for(l=atmInd+1;l<=atmInd+da;l++){
    atmCord[0][l]= atmCord[0][l-1]+Sz;
    theta[l]=theta[l-1]+Stheta;
    atmCord[1][l]=r*cos(theta[l]);
    atmCord[2][l]=r*sin(theta[l]);
      }
      //the first atom A1 im current spiral
      atmCord[0][atmInd]=atmCord[0][l-1];
      atmCord[1][atmInd]=atmCord[1][l-1];
      atmCord[2][atmInd]=atmCord[2][l-1];
      theta[atmInd]=theta[l-1];

      atmInd++;
      //the second atom B1 im current spiral
      atmCord[0][atmInd]=atmCord[0][atmInd-1]-b*sin(M_PI/6-phi);
      theta[atmInd]=b*cos(M_PI/6-phi)/r+theta[atmInd-1];
      atmCord[1][atmInd]=r*cos(theta[atmInd]);
      atmCord[2][atmInd]=r*sin(theta[atmInd]);

    }
    //=========================================================

    //determine the coordinates in the same spiral by the second former
atom's

    for(j=j0;j<j1;j++){
      atmInd++;
      atmCord[0][atmInd]=atmCord[0][atmInd-2]+Sz;
      if(atmCord[0][atmInd]>=ucL||fabs(atmCord[0][atmInd]-ucL)<=1E-6){
    atmInd--;
    break;
      }
      theta[atmInd]=theta[atmInd-2]+Stheta;
      atmCord[1][atmInd]=r*cos(theta[atmInd]);
      atmCord[2][atmInd]=r*sin(theta[atmInd]);
    }
  }
  ucAtmN=atmInd;


  //=================create other cells==========================
  for(i=1;i<cellN;i++){
    for(j=0;j<ucAtmN;j++){
      atmInd++;
      //translate former cell to create a new cell
      atmCord[0][atmInd]=atmCord[0][atmInd-ucAtmN]+ucL;
      atmCord[1][atmInd]=atmCord[1][atmInd-ucAtmN];
      atmCord[2][atmInd]=atmCord[2][atmInd-ucAtmN];
    }
  }
  //=============================================================
  ttlAtmN=atmInd;
  //produce a Guassian input file (.GJC)
  fp=fopen(strcat(fileNm,".GJC"),"w");
  fprintf(fp,"  0  1");
  for(i=1;i<=ttlAtmN;i++){
    fprintf(fp,"\n%c
0%12f%12f%12f",atmSbl,atmCord[1][i],atmCord[2][i],atmCord
[0][i]);
  }
  fclose(fp);
  printf("Some parameters:\n");
  printf("\nOriginal radius: %f\n",r);
  printf("\nUnit cell length=%f\nNumber of cells:%i\n",ucL,cellN);
  printf("\nTotle atom number: %i\nNumber of Atoms in a unit cell:
%i\n",ttlAt
mN,ucAtmN);
  printf("\nThe file %s has been created successfully!\n",fileNm);
}


--
╔═══════════════════╗
║★★★★★友谊第一  比赛第二★★★★★║
╚═══════════════════╝

※ 来源:．哈工大紫丁香 bbs.hit.edu.cn [FROM: 202.118.229.162]