Program tubes Implicit none integer I, N, L, M, NP,J,n1,l1,k,np1,v,ncar,opt double precision R1, theta,per parameter(n=80,l=80,ncar=2500) double precision x_hc(4*n*l),Z_hc(4*n*l) double precision COORD3, COORD2,COORD1, rot_x, rot_z,pi,sqt double precision eps,R(ncar,3),R3(ncar) double precision dnn integer ind(ncar) character chaine parameter (eps=1.0d-14,dnn=1.42d0) common /parametres/R1 common/angle/theta pi=4.d0*datan(1.d0) chaine=' ' k=0 sqt=dsqrt(3.d0) coord1=0.d0 coord2=0.d0 coord3=0.d0 np=1 call initialisation(L1,M,np1,OPT) write(*,'(a)') ' Computing the parameters...' call param(L1,m,per,n1) write(*,'(a12,f12.6)') ' rayon (A) :', r1 write(*,'(a12,f12.6)') ' rayon (au):', r1/0.529177d0 write(*,'(a14,f12.6)') ' periode (A) :', per write(*,'(a14,f12.6)') ' periode (au):', per/0.529177d0 write(*,'(a22,f12.6)') ' chiral angle (deg.) : ', theta*180.d0/pi write(chaine,'(I1)') OPT open(unit=20,file='CLUSTER'//chaine) if(n1.gt.ncar)then write(*,'(a,i5,a,i3,a)') ' *** ',n1,' > ',ncar,' => increase ncar .in main program ***' stop endif write(*,'(a)') ' Constructing one period of the tube...' c construction du reseau hexagonal do j=1,n do i=1,l x_hc(i+4*l*(j-1))=(dfloat(i-l/2)-1.d0)*sqt x_hc(i+l+4*l*(j-1))=sqt*(dfloat(i-l/2)+.5d0) x_hc(i+2*l+4*l*(j-1))=sqt*(dfloat(i-l/2)+.5d0) x_hc(i+3*l+4*l*(j-1))=(dfloat(i-l/2)-1.d0)*sqt z_hc(i+4*l*(j-1))=-1.d0+(dfloat(j-n/2)-1.d0)*3.d0 z_hc(i+l+4*l*(j-1))=-.5d0+(dfloat(j-n/2)-1.d0)*3.d0 z_hc(i+2*l+4*l*(j-1))=.5d0+(dfloat(j-n/2)-1.d0)*3.d0 z_hc(i+3*l+4*l*(j-1))=1.d0+(dfloat(j-n/2)-1.d0)*3.d0 enddo enddo v=0 do i=1,4*n*l COORD1=rot_x(x_hc(i)*dnn,z_hc(i)*dnn) coord3=rot_z(x_hc(i)*dnn,z_hc(i)*dnn) if(coord3.ge.-0.2d0.and. coord3.le.np*per-0.2d0 .and. &coord1.ge.0.d0.and.coord1.le.2.d0*pi*r1-eps) then call folding(COORD1,coord2) v=v+1 R(v,1)=coord1 R(v,2)=coord2 R(v,3)=coord3 k=k+1 endif enddo write(*,*) 'k,n1',k,n1 c If(k.ne.n1) stop '***Problem in the construction***' write(20,*) write(20,'(a,i2,a,i2,a)') 'tubule (',l1,',',m,') constructed ' write(20,'(a)') 'by tubule2.f, vm, lps' write(20,*) v*np1 write(*,'(a)') ' Sorting the coordinates...' do i=1,k r3(i)=r(i,3) enddo call indexx(k,R3,ind) write(*,'(a)') ' Duplicating the cell...' do i=1,np1 do j=1,v * write(20,'(a,2x,3(f16.8,2x))') ' C ', * & (R(ind(j),1)-r1)/0.529177, * & R(ind(j),2)/0.529177, * & (R(ind(j),3)+(i-1)*PER)/0.529177 write(20,'(a,2x,3(f16.8,2x))') ' C ', & (R(ind(j),1)-r1), & R(ind(j),2), & (R(ind(j),3)+(i-1)*PER) enddo enddo write(*,*) '***CONSTRUCTION SUCCESFULL***' write(*,'(a)') ' Tubule axis is the third coordinate' write(*,'(a,I4,a)') ' Output file contains ',k*np1,' atoms' close(20) end c COMPUTING THE PARAMETERS c input : c l, m : coordinate of chirtal vector c output : c a: rayon c b: chiral angle C c: length of translational vector c d: number of atoms per 1D unit cell subroutine param(l,m,c,d) implicit none double precision a,b,c integer l,m integer d double precision fac integer hcd integer r0,r1,r2 double precision radic,pi double precision dnn parameter(dnn=1.42d0) common/angle/b common/parametres/a pi=4.d0*datan(1.d0) c Computing the hcd of l and m If (M.eq.0) then hcd=l Else r0=L r1=M 5 r2=mod(r0,r1) If (r2.eq.0.) then hcd=r1 goto 18 Else r0=r1 r1=r2 Endif goto 5 18 hcd=abs(hcd) endif if(mod(l-m,3*hcd).eq.0)then fac=3.*dabs(dfloat(hcd)) else fac=dabs(dfloat(hcd)) endif c Computing rayon a a=dsqrt(3.d0)*radic(l,m)/pi/2.d0*dnn c Computing chiral angle b b=datan(dsqrt(3.d0)*dfloat(m)/(2.d0*dfloat(l)+dfloat(m))) c Computing length of translational vector c c=3.d0*radic(l,m)/fac *dnn c Computing Number of atoms in a 1d cell d d=4*(l**2+m**2+l*m)/int(fac) end c --------------------------------------------------------------- subroutine FOLDING(X,Y) double precision x double precision y double precision phi double precision R1 common /parametres/R1 phi=X/R1 X=R1*(1.d0-dcos(phi)) Y=R1*dsin(phi) return end c x FUNCTION Rotation of an angle alpha double precision function rot_x(x,y) double precision x, y, alpha common /angle/alpha rot_x=x*dcos(alpha)+y*dsin(alpha) return end c z FUNCTION Rotation of an angle alpha double precision function rot_z(x,y) double precision x, y, alpha common /angle/alpha rot_z=y*dcos(alpha)-x*dsin(alpha) return end c FUNCTION COMPUTING THE "OLD SAME RADICAL" double precision function radic(l,m) integer l,m radic=sqrt(dfloat(l)**2+dfloat(m)**2+dfloat(l)*dfloat(m)) end subroutine initialisation(L,M,np,opt) integer L,M,np,opt c character*8 chaine write(*,'(a,$)') ' Extension of file ? ' read(*,*) opt write(*,'(a,I1)') ' the output file is called CLUSTER',opt write(*,*) write(*,'(/,a,$)') ' Enter L : ' read(*,*) L write(*,'(a,$)') ' Enter M : ' read(*,*) M write(*,'(a,$)') ' Number of period(s) : ' READ(*,'(I5)') np end SUBROUTINE indexx(n,arr,indx) INTEGER n,indx(n),M,NSTACK REAL*8 arr(n) PARAMETER (M=7,NSTACK=50) INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK) REAL*8 a do 11 j=1,n indx(j)=j 11 continue jstack=0 l=1 ir=n 1 if(ir-l.lt.M)then do 13 j=l+1,ir indxt=indx(j) a=arr(indxt) do 12 i=j-1,l,-1 if(arr(indx(i)).le.a)goto 2 indx(i+1)=indx(i) 12 continue i=l-1 2 indx(i+1)=indxt 13 continue if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 itemp=indx(k) indx(k)=indx(l+1) indx(l+1)=itemp if(arr(indx(l)).gt.arr(indx(ir)))then itemp=indx(l) indx(l)=indx(ir) indx(ir)=itemp endif if(arr(indx(l+1)).gt.arr(indx(ir)))then itemp=indx(l+1) indx(l+1)=indx(ir) indx(ir)=itemp endif if(arr(indx(l)).gt.arr(indx(l+1)))then itemp=indx(l) indx(l)=indx(l+1) indx(l+1)=itemp endif i=l+1 j=ir indxt=indx(l+1) a=arr(indxt) 3 continue i=i+1 if(arr(indx(i)).lt.a)goto 3 4 continue j=j-1 if(arr(indx(j)).gt.a)goto 4 if(j.lt.i)goto 5 itemp=indx(i) indx(i)=indx(j) indx(j)=itemp goto 3 5 indx(l+1)=indx(j) indx(j)=indxt jstack=jstack+2 if(jstack.gt.NSTACK)pause 'NSTACK too small in indexx' if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 END