一个生成碳纳米管坐标的fortran程序

一个生成碳纳米管坐标的fortran程序

     欧阳剑   2008 年10月   ouyjian0426@163.com

 

在编写过程中,先是自己用作图软件coredraw画了一个二维蜂房结构图,在图书馆找了本一个叫成会明的人写的关于碳纳米管的书,了解了碳纳米管的一些结构参数的表示,在图上琢磨了一阵,写了这个程序,还发现没什么大错误.

!!!!!!!!!!!!以下是代码!!!!!!!!!!!!!!!!!!!!!!

!!!!运行前,创建一个名为 mn.in 的文件,用记事本打开,在里面输入m和n,用空格隔开,保存.

!!!!运行得到的xyz.dat 文件就是碳纳米管中碳原子的坐标.可以用orign软件生成图形.

!!! 目前没有去考虑c-c键长

!!! 还有就是不知道怎么去把相邻的碳原子用线连起来.

 

program ouyjian’s carbon nanotubes

implicit none

integer:: m,n,ma,na,maa,naa

real:: mb,nb,lb2,lb

real::l2,l,cosa,cosb,cosc,la2,la,c2

integer::mt,nt,dr

integer::gcd,dmn

real::lt2,lt,angc,sinc,angd,za,xa,ya,r,pi,xaa,yaa,zaa

!ange is the angle between (mb,nb)and a1

!angf is the angle between (mb,nb)and (m,n)

real::ange,cose,cosf,xb,yb,zb,anga,sinf,angg

open (unit = 15, file='mn.in')

open (unit = 5, file='mn.out')

open (unit = 9,file='xyz.dat')

Read (15,*) m,n

!l is the distance of (m,n),and it's the tube's girth

l2=m*m+n*n+m*n

l=sqrt(l2)

pi = 3.1415926535

!r is rhe radious of the tube

r=l/(2*pi)

ma = 0

na = 0

zaa=0.0

xaa=0.0

yaa=0.0

 dmn = gcd(m,n)

 if(mod((m-n),3)==0) then

  dr = 3*dmn

 else

  dr = dmn

    endif

mt = (2*n+m)/dr  

nt = -(2*m+n)/dr

lt2 = mt*mt +nt*nt + mt*nt

lt = sqrt(lt2)

cosa=(2*m+n)/(2*l)

anga=acos(cosa)

 

ma_loop:do

ma = ma + 1

  na = nt - 1

  na_loop:do

  na = na + 1

!la is the distance of (ma,na)

la2=ma*ma+na*na+ma*na

la=sqrt(la2)

!c is the distance between (m,n)and (ma,na)

c2=(m-ma)*(m-ma)+(n-na)*(n-na)+(m-ma)*(n-na)

cosb=(2*ma+na)/(2*la)

!angc is the angle between (m,n) and (ma,na)

cosc=(la2+l2-c2)/(2*la*l)

angc=acos(cosc)

sinc=sin(angc)

za=la*sinc

angd=la*cosc*2*pi/l

xa=r*sin(angd)

ya=-r*cos(angd)

  if(za<=(lt+0.001).and.la*cosc<=l.and.cosc>=0) then

    if(na<=0) then

    maa = ma

    naa = na

    xaa = xa

    yaa = ya

    zaa = za

    mb = maa - 0.333333333

    nb = naa - 0.333333333

    lb2 = mb*mb + nb*nb + mb*nb

    lb = sqrt(lb2)

    cose = (2*mb+nb)/(2*lb)

    ange = acos(cose)

    if(nb<=0)then

    cosf = cos(anga+ange)

    sinf = sin(anga+ange)

    else

      cosf=cos(anga-ange)

      sinf=sin(anga-ange)

      endif

      angg = 2*pi*lb*cosf/l

      xb = r*sin(angg)

      yb = -r*cos(angg)

      zb = lb*sinf 

    write(5,*)ma,na

    write(5,*)mb,nb

    write(9,*) xaa,'       ',yaa,'      ',zaa

    write(9,*) xb,' ',yb,' ',zb,

    else if(na>0.and.cosc>=cosa.and.cosb>=cosa) then

    maa = ma

    naa = na

    xaa = la*cosc

    xaa = xa

    yaa = ya

    zaa = za

    mb = maa - 0.333333333

    nb = naa - 0.333333333

        lb2 = mb*mb + nb*nb + mb*nb

    lb = sqrt(lb2)

    cose = (2*mb+nb)/(2*lb)

    ange = acos(cose) 

    if(nb<=0)then

    cosf = cos(anga+ange)

    sinf = sin(anga+ange)

    else

      cosf=cos(anga-ange)

      sinf=sin(anga-ange)

      endif

     

      angg = 2*pi*lb*cosf/l

      xb = r*sin(angg)

      yb = -r*cos(angg)

      zb = lb*sinf

    write(5,*)ma,na

    write(5,*)mb,nb

    write(9,*) xaa,'       ',yaa,'      ',zaa

    write(9,*) xb,' ',yb,' ',zb,

     endif

   endif

  if ( na >= n) exit na_loop

 

if ( ma >= (mt + m + 1)) exit ma_loop

  end do na_loop

end do ma_loop

close(15)

close(5)

close(9)

 

end program ouyjian’s carbon nanotubes

 

 

integer function gcd(mm,nn) result(cd)

if (mm < nn)  then

   cc = mm; mm = nn; nn = cc

end if

i = 0

do while(i <= mm)

  i = i + 1

  if(mod(mm,i) == 0.and.mod(nn,i)==0) then

    cd = i

    end if

  end do

end function gcd