CODE: PROGRAM MAIN
PARAMETER(num=14797,blx=61.634002,bly=62.406401,blz=63.429759,
* namino=200)
CHARACTER name(num)*9,fft(num)*5,atomname(num)*4,fx(num)*4
INTEGER occupation(num),Tatom,itt,i
real xc(num),yc(num),zc(num),charge(num)
real xtemp,ytemp,ztemp,xfinal,yfinal,zfinal
OPEN(10,file='2.car',status='old')
do 20 i=1,num
read(10,*)name(i),xc(i),yc(i),zc(i),fx(i),occupation(i),
* fft(i),atomname(i),charge(num)
20 continue
close(10)
do 30 i=1,num
xc(i)=0
yc(i)=0
zc(i)=0
30 continue
natom_add=0
Tatom=num+natom_add
*********************************************************
* add Si to the chosen Oxygen atom
*********************************************************
do 45 itt=1,namino
xtemp=blx*RAN2(IDUM)
ytemp=bly*RAN2(IDUM)
ztemp=blz*RAN2(IDUM)
lbond=1.90
call addatom(xtemp,ytemp,ztemp,lbond,Tatom,xfinal,yfinal,zfinal)
natom_add=natom_add+1
Tatom=Tatom+natom_add
name(Tatom)='Si'
xc(Tatom)=xfinal
yc(Tatom)=yfinal
zc(Tatom)=zfinal
fx(Tatom)='XXXX'
occupation(Tatom)=1
fft(Tatom)='Si3'
atomname(Tatom)='Si'
charge(Tatom)=5.000
open(140,file='amino.car',access='append')
write(140,888)a(Tatom),xc(Tatom),
& yc(Tatom),zc(Tatom),
& fx(Tatom),occupation(Tatom),
& fft(Tatom),atomname(Tatom),
& charge(Tatom)
close(140)
******************************************************************
* add C1 atom on the chosen Si atom
******************************************************************
lbond=1.93
xtemp=xc(Tatom)
ytemp=yc(Tatom)
ztemp=zc(Tatom)
call addatom(xtemp,ytemp,ztemp,lbond,Tatom,xfinal,yfinal,
* zfinal)
natom_add=natom_add+1
Tatom=Tatom+natom_add
name(Tatom)='C'
xc(Tatom)=xfinal
yc(Tatom)=yfinal
zc(Tatom)=zfinal
fx(Tatom)='XXXX'
occupation(Tatom)=1
fft(Tatom)='C_3'
atomname(Tatom)='C'
charge(Tatom)=6.000
open(150,file='amino.car',access='append')
write(150,888)a(Tatom),xc(Tatom),
& yc(Tatom),zc(Tatom),
& fx(Tatom),occupation(Tatom),
& fft(Tatom),atomname(Tatom),
& charge(Tatom)
close(150)
******************************************************************
* add C2 atom on the chosen C1 atom
******************************************************************
lbond=1.50
xtemp=xc(Tatom)
ytemp=yc(Tatom)
ztemp=zc(Tatom)
call addatom(xtemp,ytemp,ztemp,lbond,Tatom,
& xfinal,yfinal,zfinal)
natom_add=natom_add+1
Tatom=Tatom+natom_add
name(Tatom)='C'
xc(Tatom)=xfinal
yc(Tatom)=yfinal
zc(Tatom)=zfinal
fx(Tatom)='XXXX'
occupation(Tatom)=1
fft(Tatom)='C_3'
atomname(Tatom)='C'
charge(Tatom)=7.000
open(160,file='amino.car',access='append')
write(160,888)a(Tatom),xc(Tatom),
& yc(Tatom),zc(Tatom),
& fx(Tatom),occupation(Tatom),
& fft(Tatom),atomname(Tatom),
& charge(Tatom)
close(160)
******************************************************************
* add C3 atom on the chosen C2 atom
******************************************************************
lbond=1.50
xtemp=xc(Tatom)
ytemp=yc(Tatom)
ztemp=zc(Tatom)
call addatom(xtemp,ytemp,ztemp,lbond,Tatom,
& xfinal,yfinal,zfinal)
natom_add=natom_add+1
Tatom=Tatom+natom_add
name(Tatom)='C'
xc(Tatom)=xfinal
yc(Tatom)=yfinal
zc(Tatom)=zfinal
fx(Tatom)='XXXX'
occupation(Tatom)=1
fft(Tatom)='C_3'
atomname(Tatom)='C'
charge(Tatom)=8.000
open(170,file='amino.car',access='append')
write(170,888)a(Tatom),xc(Tatom),
& yc(Tatom),zc(Tatom),
& fx(Tatom),occupation(Tatom),
& fft(Tatom),atomname(Tatom),
& charge(Tatom)
close(170)
******************************************************************
* add N atom on the chosen C3 atom
******************************************************************
lbond=1.53
xtemp=xc(Tatom)
ytemp=yc(Tatom)
ztemp=zc(Tatom)
call addatom(xtemp,ytemp,ztemp,lbond,Tatom,
& xfinal,yfinal,zfinal)
natom_add=natom_add+1
Tatom=Tatom+natom_add
name(Tatom)='N'
xc(Tatom)=xfinal
yc(Tatom)=yfinal
zc(Tatom)=zfinal
fx(Tatom)='XXXX'
occupation(Tatom)=1
fft(Tatom)='N_3'
atomname(Tatom)='N'
charge(Tatom)=9.000
open(180,file='amino.car',access='append')
write(180,888)a(Tatom),xc(Tatom),
& yc(Tatom),zc(Tatom),
& fx(Tatom),occupation(Tatom),
& fft(Tatom),atomname(Tatom),
& charge(Tatom)
close(180)
45 continue
*********************************************************************************
* write the output file
*********************************************************************************
open(190,file='output-final.car',access='append')
do 200 i=1,Tatom
write(190,888)a(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i),
& atomname(i),charge(i)
200 continue
close(190)
888 format(A5,3X,F12.9,2X,F13.9,3X,F12.9,1X,A4,1X,I1,6X,A5,3X,A2,F7.3)
*********************************************************************************
*********************************************************************************
end
**********************************add atom************************************
SUBROUTINE addatom(xtemp,ytemp,ztemp,lbond,Tatom,xfinal,yfinal,
* zfinal)
REAL center(14,2)
INTEGER Tatom
open(35,file='center.txt',status='old')
do 36 i=1,14
read(35,*)center(i,1),center(i,2)
36 continue
close(35)
lc=1000
xcenter=0
ycenter=0
do 40 j=1,14
lct=sqrt((xtemp-center(i,1))**2+(ytemp-center(i,2)**2))
if(lct.lt.lc)then
lc=lct
*******把整个体分解成无数的平行的平面,确定整个平面的每个分割区域的中心*******************
xcenter=center(i,1)
ycenter=center(i,2)
end if
40 continue
if(ytemp.le.ycenter)then
ja=20
jb=0
jc=-1
jd=j
else
ja=-20
jb=0
jc=1
jd=j
endif
*****************向负方向搜索**********************************
if(xtemp.le.xcenter-5)then
ia=20
ib=0
ic=-1
minlbond=lbond-0.1
maxlbond=lbond+0.1
search_step=0.1
do 1000 i=ia,ib,ic
do 1010 j=ja,jb,jc
do 1020 k=-20,20
targetx=xtemp+search_step*i
targety=ytemp+search_step*j
targetz=ztemp+search_step*k
dis1=sqrt((targetx-xtemp)**2+
& (targety-ytemp)**2+
& (targetz-ztemp)**2)
if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then
do 1030 im=1,Tatom
dis2=sqrt((targetx-xc(im))**2+
& (targety-yc(im))**2+
& (targetz-zc(im))**2)
if(dis2.lt.3.0)then
goto 1020
end if
1030 continue
goto 1040
end if
1020 continue
1010 continue
1000 continue
1040 xfinal=targetx
yfinal=targety
zfinal=targetz
end if
*****************向正方向搜索**********************************
if(xtemp.ge.xcenter+5)then
ia=20
ib=-20
ic=-1
minlbond=lbond-0.1
maxlbond=lbond+0.1
search_step=0.1
do 1050 i=ia,ib,ic
do 1060 j=ja,jb,jc
do 1070 k=-20,20
targetx=xtemp+search_step*i
targety=ytemp+search_step*j
targetz=ztemp+search_step*k
dis1=sqrt((targetx-xtemp)**2+
& (targety-ytemp)**2+
& (targetz-ztemp)**2)
if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then
do 1080 im=1,Tatom
dis2=sqrt((targetx-xc(im))**2+
& (targety-yc(im))**2+
& (targetz-zc(im))**2)
if(dis2.lt.3.0)then
goto 1070
endif
1080 continue
end if
goto 1090
1070 continue
1060 continue
1050 continue
1090 xfinal=targetx
yfinal=targety
zfinal=targetz
end if
end SUBROUTINE
******************************************************************************
FUNCTION RAN2(IDUM)
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
*IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
*NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if (idum.le.0) then
idum=max(-idum,1)
idum2=idum
do 11 j=NTAB+8,1,-1
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
if (j.le.NTAB) iv(j)=idum
11 continue
iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if (idum2.lt.0) idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
END