Physics 版 (精华区)
发信人: zjliu (秋天的萝卜), 信区: Physics
标 题: [转载] CNT的坐标生成程序
发信站: 哈工大紫丁香 (Thu Dec 4 16:05:48 2003), 站内信件
网上下的
好像是Brenner的
------------------------------------------------------------
program MkTube
implicit double precision (a-h,o-z)
parameter (MAXATM=100000)
parameter (BOND=1.418d0)
parameter (TOL=1.d-3)
character*10 filename,filename1
dimension rnuc(3,MAXATM), rout(3)
dimension avec(2),bvec(2),pvec(2),zvec(2)
dimension ur(2,2),uk(2,2)
character*80 line
c dimension iprim(2,6)
c data iprim/1,0, 0,1, 1,-1, -1,0, 0,-1, -1,1/
XMOD(A,B) = A - B*NINT(A/B)
111 PI = DACOS(-1.d0)
SCALE = DSQRT(3.d0)*BOND
ur(1,1) = 1.d0
ur(2,1) = 0.d0
ur(1,2) = 0.5d0
ur(2,2) = DSQRT(3.d0)/2.d0
do k=1,2
uk(k,1) = 2*ur(k,1) - ur(k,2)
uk(k,2) = 2*ur(k,2) - ur(k,1)
enddo
write(*,*) 'enter n1,n2'
read (*,*) n1,n2
write(*,*) 'enter number repeat units'
read (*,*) nunit
if ((n1.eq.0).and.(n2.eq.0)) then
print *,'n1,n2 both zero invalid'
stop
endif
nfile = 1
filename = 'c'
filename1 = 'c'
if (n1.lt.10) then
write (filename(nfile+1:),'(i1)') n1
write (filename1(nfile+1:),'(i1)') n1
nfile = nfile+1
else
write (filename(nfile+1:),'(i2)') n1
write (filename1(nfile+1:),'(i2)') n1
nfile = nfile+2
endif
nfile = nfile+1
filename(nfile:nfile) = 'x'
filename1(nfile:nfile) = 'x'
if (n2.lt.10) then
write (filename(nfile+1:),'(i1)') n2
write (filename1(nfile+1:),'(i1)') n2
nfile = nfile+1
else
write (filename(nfile+1:),'(i2)') n2
write (filename1(nfile+1:),'(i2)') n2
nfile = nfile+2
endif
filename(nfile+1:) = '.plu'
filename1(nfile+1:) = '.xyz'
open (unit=7,file=filename,status='unknown',
& form='formatted',carriagecontrol='list') !,access='append')
open (unit=1,file=filename1,status='unknown',
& form='formatted',carriagecontrol='list') !,access='append')
! file name determination
open (unit=8,file='innertube',status='unknown',
& form='formatted',carriagecontrol='list')
open (unit=9,file='outertube',status='unknown',
& form='formatted',carriagecontrol='list')
avec(1) = n1*ur(1,1) + n2*ur(1,2)
avec(2) = n1*ur(2,1) + n2*ur(2,2)
adist = DSQRT(avec(1)**2+avec(2)**2)
radius = adist/(2.d0*PI)
m1 = -n2
m2 = n1
bvec(1) = m1*uk(1,1) + m2*uk(1,2)
bvec(2) = m1*uk(2,1) + m2*uk(2,2)
bdist = DSQRT(bvec(1)**2+bvec(2)**2)
do k=1,2
avec(k) = avec(k)/adist
bvec(k) = bvec(k)/bdist
enddo
zmin = 2.d0
phi = 2*PI
xstep = 0.d0
imin = 0
i = 0
j1 = 0
j2 = 0
nmax = MAX(ABS(n1),ABS(n2))
do i1=-nmax,nmax
do i2=-nmax,nmax
i = i + 1
zvec(1) = i1*ur(1,1) + i2*ur(1,2)
zvec(2) = i1*ur(2,1) + i2*ur(2,2)
z = zvec(1)*bvec(1) + zvec(2)*bvec(2)
zz = zvec(1)*avec(1) + zvec(2)*avec(2)
if (z.gt.TOL) then
if ((z.lt.zmin).or.((ABS(z-zmin).lt.TOL).and.
& (ABS(zz/radius).lt.ABS(phi)))) then
imin = i
zmin = z
j1 = i1
j2 = i2
xstep = zz
phi = XMOD(zz/radius,2*PI)
endif
endif
enddo
enddo
if (imin.eq.0) then
print *,'helical twist not found'
stop
endif
c print *,'Indices of helical step operation ',j1,j2
NN = j1*n2 - j2*n1
c print *,'Test sum: ',NN
cphi1 = (FLOAT(j1*n1+j2*n2)+0.5*FLOAT(j1*n2+j2*n1))
cphi2 = FLOAT(n1*n1+n2*n2+n1*n2)
cphi = cphi1/cphi2*XMOD(DFLOAT(n2-n1),3.0d0)/3.0
cphi = cphi + FLOAT(j1-j2)/3.0
cphi = XMOD(2*cphi,2.d0)
h2 = j1*j1 + j2*j2 + j1*j2
b2 = n1*n1 + n2*n2 + n1*n2
hb = 0.5d0*(2*n1*j1+2*n2*j2+n1*j2+n2*j1)
xscale = SQRT((h2*b2-hb*hb)/b2)
c print *,'Relative phase of crossing:',cphi*PI,'(',cphi,'*PI)'
c print *,'Scale factor: ',xscale
VOL = DSQRT(3.d0)/2.d0
amotif = adist*zmin/VOL
bmotif = adist*bdist/VOL
c print *,'Number of dimers in periodic unit cell: ',bmotif
motif = NINT(amotif)
nmotif = NINT(bmotif)/motif
xmotif = bmotif/amotif
c
c print *,'Number of motifs in periodic unit cell: ',xmotif
xtwist = bdist*xstep/adist/zmin
c print *,'Number of turns per periodic repeat: ',xtwist
zunit = zmin*SCALE
n = 1
rnuc(1,n) = 0.d0
rnuc(2,n) = 0.d0
nmax = MAX(n1,n2)
uk(1,1) = 1.d0/adist
uk(2,1) = -xstep/zmin/adist
uk(1,2) = 0.d0
uk(2,2) = 1.d0/zmin
do i1=0,ABS(n1)
do i2=0,ABS(n2)
if (n.ge.motif) goto 200
zvec(1) = i1*ur(1,1) + i2*ur(1,2)
zvec(2) = i1*ur(2,1) + i2*ur(2,2)
if (zvec(1).lt.0.d0) then
zvec(1) = -zvec(1)
zvec(2) = -zvec(2)
endif
pvec(1) = avec(1)*zvec(1) + avec(2)*zvec(2)
pvec(2) = bvec(1)*zvec(1) + bvec(2)*zvec(2)
nn1 = NINT(pvec(1)*uk(1,1)+pvec(2)*uk(2,1))
nn2 = NINT(pvec(1)*uk(1,2)+pvec(2)*uk(2,2))
pvec(1) = pvec(1) - nn1*adist - nn2*xstep
pvec(2) = pvec(2) - nn2*zmin
do i=1,n
do k=1,2
zvec(k) = pvec(k) - rnuc(k,2*i-1)
enddo
aa1 = zvec(1)*uk(1,1) + zvec(2)*uk(2,1)
aa2 = zvec(1)*uk(1,2) + zvec(2)*uk(2,2)
if ((ABS(XMOD(aa1,1.d0)).lt.TOL).and.
& (ABS(XMOD(aa2,1.d0)).lt.TOL)) goto 100
enddo
n = n + 1
rnuc(1,2*n-1) = pvec(1)
rnuc(2,2*n-1) = pvec(2)
100 continue
enddo
enddo
200 xshift = avec(2)/2.d0/DSQRT(3.d0)
yshift = bvec(2)/2.d0/DSQRT(3.d0)
do i=1,motif
xx = rnuc(1,2*i-1)
zz = rnuc(2,2*i-1)
theta = (xx+xshift)/radius
rnuc(1,2*i-1) = SCALE*radius*DCOS(theta)
rnuc(2,2*i-1) = SCALE*radius*DSIN(theta)
rnuc(3,2*i-1) = SCALE*(zz+yshift)
theta = (xx-xshift)/radius
rnuc(1,2*i) = SCALE*radius*DCOS(theta)
rnuc(2,2*i) = SCALE*radius*DSIN(theta)
rnuc(3,2*i) = SCALE*(zz-yshift)
enddo
natom = 2*motif
c print *,'Number of carbons per unit cell: ',natom
c print *,'Number of carbons per unit cell: ',2*amotif
c print *,'helical axis step distance: ',zunit
c print *,'twist angle: ',phi*180.d0/PI
c print *,'cylinder radius: ',radius*SCALE
xblock = 1.23d0/zunit
c ----------------------------------------------
c write stuff for simulation
write(7,17000) n1,n2
write(7,18000) natom*nunit,0,0,0
write(7,19000) 0.0,1.0
write(7,19000) 1.0d9,1.0d9,zunit*(nunit)
c write stuff for xmol
write(1,*) natom*nunit
write(1,*) ' '
c write(7,19000) 1.0d9,1.0d9,1.0d09
nn = 0
do n=1,nunit
do i=1,natom
call screw(n,rnuc(1,i),rout,phi,zunit)
nn = nn+1
if (nn.lt.10) then
nitem = 1
elseif (nn.lt.100) then
nitem = 2
else
nitem = 3
endif
c write (7,13000) nn,rout
c write (7,20000) nn,6,rout,1
c SHIFT X,Y COORDINATE
do j=1,2
rout(j) = 10+rout(j)
end do
!rout(3)=100+rout(3)
!if (n1 == 5) then
! write inner tube coordination 'type 1'
! rout(3) = 100-6.14878+rout(3) !(100-6.14878)
! write (8,22000) 400+nn,1,12.01,rout
!
! else if (n1 == 10) then
! ! write outer tube coordination 'type 0'
! rout(3) = 100+rout(3)
! write (9,22000) nn,0,12.01,rout
! end if
!rout(3) = rout(3)
omi = 100
vx = omi*(rout(2)-10)
vy = omi*(rout(1)-10)
! Inner tubule
rout(3)=rout(3)-1.842036
!if (rout(3) .GT. -1E-10) then
!write (7,24000) nn,1,12.01,rout,vx,vy,0.0
!write (1,'(i5,3f12.6)') 6,rout
write (7,24000) 130+nn,0,12.01,rout,0.0,0.0,0.0
!else
!nn = nn-1
!end if
3500 format(i3,3f12.6)
enddo
enddo
c write other stuff
c do n=1,4
c do j=1,nn
c write(7,21000) j,0.0,0.0,0.0
c enddo
c enddo
close (7)
close (1)
stop
10500 format(a,1x,i3)
11000 format(f20.8)
12000 format(f20.12)
13000 format('C',i<nitem>,3x,3f20.12)
14000 format(q,a)
15000 format(f12.5)
16000 format(a)
17000 format('nanotube with n1,n2= ',2i3)
18000 format(4i6)
19000 format(3e20.11)
20000 format(2i5,3e20.11,i3)
21000 format(I5,3E20.11)
22000 format(2i5,4e20.11)
23000 format(i5,2a5,i5,3f8.3,3f12.6)
24000 format(2i5,4e20.11,3e20.11)
end
subroutine screw(n,in,out,twist,ucell)
implicit double precision (a-h,o-z)
integer n
double precision in(3),out(3)
if (n.eq.0) then
do i=1,3
out(i) = in(i)
enddo
else
zc = DCOS(n*twist)
zs = DSIN(n*twist)
out(1) = zc*in(1) - zs*in(2)
out(2) = zs*in(1) + zc*in(2)
out(3) = in(3) + float(n)*ucell
endif
return
end
--
╔═══════════════════╗
║★★★★★友谊第一 比赛第二★★★★★║
╚═══════════════════╝
※ 来源:．哈工大紫丁香 bbs.hit.edu.cn [FROM: 202.118.229.162]
Powered by KBS BBS 2.0 (http://dev.kcn.cn)
页面执行时间：207.316毫秒