c
c Atomic Coordinates in the unit cell of a (n,m) nanotube
c
c Made by 10/31/95   by T. Takeya and R. Saito
c
c Reference: ``Physical Properties of Carbon Nanotubes''
c             Imperial College Press, 
c             by R. Saito, G. Dresselhaus, and M.S. Dresselhaus
c 
c Input: n,m  (from keyboard)
c Output: tube.xyz  (coordinates. See detail in the above book.)
c         en.xyz2   (a work file for the program A-3)
c
      implicit real*8(a-h,o-z)
c
c  acc: C-C bond distance.
c
      parameter (acc = 1.42d0)
      parameter (nk=1300)
      dimension x(nk),y(nk),z(nk)
c
      write(*,*)' Enter n,m for C_h=(n,m)'
      read(*,*) n,m
c
      call gen11(n,m,np,nq,ndr)
       write(*,*)' R = (', np, ',', nq , ')'
c     
c     a: the length of the unit vector
c     pi = 3.141592... sq3 = 1.732....
c
      sq3=sqrt(3.0d0)
      pi=4.0d0*atan(1.0d0)
      a=sqrt(3.0d0)*acc
c
c     r;|R| , c;|C_h| , t;|T|
c
      r=a*sqrt(dfloat(np*np+nq*nq+np*nq))
      c=a*sqrt(dfloat(n*n+m*m+n*m))
      t=sqrt(3.0d0)*c/dfloat(ndr)
      write(*,*) 't=',t, dfloat(n*n+m*m+n*m)
c      write(*,190)c/2.0d0/pi
c 190  format(' radius =',f10.4)
c
c     nn: the number of hexagon in the unit cell N
c
      nn=2*(n**2+m**2+m*n)/ndr
c
c     write(*,*)' N =',nn*2
c
c     rs: radius of the tube
c
      if(2*nn.gt.nk) stop 'parameter nk is too small'
c
      rs=c/(2.0d0*pi)
c
c     q1: the chiral angle for C_h
c     q2: the 'chiral' angle for R 
c     q3: the angle between C_h and R
c
      q1=atan((sq3*dfloat(m))/dfloat(2*n+m))
      q2=atan((sq3*dfloat(nq))/dfloat(2*np+nq))
      q3=q1-q2
c
c     q4: a period of an angle for the A atom
c     q5: the difference of the angle between the A and B atoms
c
      q4=2.0d0*pi/dfloat(nn)
      q5=acc*cos((pi/6.0d0)-q1)/c*2.0d0*pi
c
c     h1: 
c     h2: Delta z between the A and B atoms
c  
      h1=abs(t)/abs(sin(q3))
      h2=acc*sin((pi/6.0d0)-q1)
c
      write(*,*) 'q1: =',q1*180.0d0/pi
      write(*,*) 'q2: =',q2*180.0d0/pi
      write(*,*) 'q4: =',q4*180.0d0/pi
      write(*,*) 'q5: =',q5*180.0d0/pi
c
c    Calculate 2*nn atoms in the unit cell.
c    
c    The A atom
c     
      ii=0
      do 100 i=0,nn-1
         x1=0
         y1=0
         z1=0
          k=int(dfloat(i)*abs(r)/h1)
         x1=rs*cos(dfloat(i)*q4)
         y1=rs*sin(dfloat(i)*q4)
         z1=dfloat((dfloat(i)*abs(r)-dfloat(k)*h1))*sin(q3)
         kk2=abs(int(z1/t))+1
c
c Check the A atom is in the unit cell 0 < z1 < t
c
          if(z1.gt.t-0.02)then
            z1=z1-t*dfloat(kk2)
            endif
         if(z1.lt.-0.02) then
            z1=z1+t*dfloat(kk2)
            endif
            ii=ii+1
         x(ii)=x1
         y(ii)=y1
         z(ii)=z1
c
c     The B atoms
c
         z3=(dfloat(i)*abs(r)-dfloat(k)*h1)*sin(q3)-h2
         ii=ii+1
c
c Check the B atom is in the unit cell 0 < z3 < t
c
         if((z3.ge.-0.02).and.(z3.le.t-0.02))then
c yes
            x2 =rs*cos(dfloat(i)*q4+q5)
            y2 =rs*sin(dfloat(i)*q4+q5)
            z2 =dfloat(dfloat(i)*abs(r)-dfloat(k)*h1)*sin(q3)-h2
            x(ii)=x2
            y(ii)=y2
            z(ii)=z2
         else
c no
            x2 =rs*cos(dfloat(i)*q4+q5)
            y2 =rs*sin(dfloat(i)*q4+q5)
            z2 =dfloat(dfloat(i)*abs(r)-dfloat(k+1)*h1)*sin(q3)-h2
            kk=abs(int(z2/t))+1
         if(z2.gt.t-0.01)then
            z2=z2-t*dfloat(kk)
            endif
         if(z2.lt.-0.01) then
            z2=z2+t*dfloat(kk)
            endif
            x(ii)=x2
            y(ii)=y2
            z(ii)=z2
             endif
 100      continue
c
c Out put to the file tube.xyz
c
          open(60,file='tube.xyz')
          write(60,*)2*nn
          write(60,*)' '
          do i=1,nn*2
             write(60,117)x(i),y(i),z(i)
             end do
 117         format('C',3f10.5)
             close(60)
c
c Out put to the file en.xyz2
c
          open(60,file='en.xyz2')
          write(60,*)2*nn
          write(60,118)t,acc
 118      format(2f25.5)
          do i=1,nn*2
             write(60,116)i,x(i),y(i),z(i)
             end do
 116         format(1i5,3f25.20)
             close(60) 
          stop
          end
c
c 10/16/95   R. Saito
c
      subroutine gen11(n,m,np,nq,ndr)
      dimension nnp(100),nnq(100)
c
      itest=1
      itest1=0
c

      nd=igcm(n,m)
      if(mod((n-m),3*nd).eq.0) then
      ndr=3*nd
      else
      ndr=nd
      endif
      if(itest.eq.1) nd=nd
      if(itest.eq.1) ndr=ndr
c
      a=sqrt(3.0)*1.421
      eps=1.0e-5
c L
      l2=n*n+m*m+n*m
      if(l2.le.0) stop 'l2.le.0'
      l=int(sqrt(dfloat(l2))+eps)
      dt=a*sqrt(dfloat(l2))/3.1415926525

c T
      nr=(2*m+n)/ndr
      ns=-(2*n+m)/ndr
      nt2=3*l2/ndr/ndr
      nt=int(sqrt(dfloat(nt2))+eps)
c N
      nn=2*l2/ndr
c R
      ichk=0
      if(nr.eq.0) then
         n60=1
         else
         n60=nr
      endif
c
      itest2=1
c
      do 60 np=-abs(n60),abs(n60)
      do 61 nq=-abs(ns),abs(ns)
      j2 = nr*nq - ns*np
      if(j2.eq.1) then
      j1 = m *np -  n*nq
      if( j1.gt.0 .and. j1.lt.nn ) then
      ichk=ichk+1
      nnp(ichk)=np
      nnq(ichk)=nq
      endif
      endif
 61   continue
 60   continue
c
      if(ichk.eq.0) then
      stop ' not found p,q strange!!'
      endif
c
      itest3=1
c
      if(ichk.ge.2) then
      if(itest3.eq.1) then
      stop ' more than 1 pair of p,q strange!!'
      endif

      if(nr.ne.0 .and. ns.ne.0) then
c
      if(itest1.eq.1) then
      do 77 i=1,ichk
      if((m*nnp(i)-n*nnq(i)).lt.nn) goto 777
 77   continue 
      endif
c
      endif
      endif
 777  continue 
c
             if(itest.eq.1) then
              np=nnp(1)
              nq=nnq(1)
             
             endif
 2    continue
      return
 1    continue
      stop

      end
c
c calculate the highest common divisor
c
      integer function igcm(ii,jj)
      i=abs(ii)
      j=abs(jj)
      if(j.gt.i) then
      iw=j
      j=i
      i=iw
      endif
      if(j.eq.0) then
      igcm=i
      return
      endif
 10   continue
      ir=mod(i,j)
      if(ir.eq.0) then
      igcm=j
      return
      else
      i=j
      j=ir
      endif
      goto 10
      end
