24小时热门版块排行榜

返回列表

zyj8119

木虫 (著名写手)

应助: 65 (初中生)
贵宾: 0.003
金币: 915.1
散金: 1440
红花: 35
帖子: 2936
在线: 1329.4小时
虫号: 664177
注册: 2008-11-29
性别: GG
专业: 理论和计算化学

[交流] 【原创】使用MS生成mcm-41，并且插胺基的程序，与大家讨论下已有2人参与

integer natom,natom0,nho,namino
************************************************************
* Number of atoms in the original structure is 9992.
* But the parameter natom should include the number of atoms added subsequently.
* So here the value of natom is set to 15000
************************************************************
   parameter (natom=15000,natom0=9992,nho=2319,namino=435)

character a(natom)*5,fx(natom)*4,fft(natom)*5,atomname(natom)*2
integer occupation(natom),kjishu,ron,nOtemp,kstop,natom_add,Tatom
integer Ohydroxy_SN(nho),Temp_SN
double precision xc(natom),yc(natom),zc(natom)
double precision rox,roy,roz,Ohydroxy(nho,4),Otemp(nho,4)
double precision OTa(4),OTb(4),temp(3),list3(nho*3,12)
integer templist3(nho*3,3),Nra
double precision xtemp,ytemp,ztemp,xfinal,yfinal,zfinal
* integer NSiT,NSi_S,NSi_C,NCT,NC_S,NC_C,NNT,NN_S,NN_C
* character NSi_CC*1,NC_CC*1,NN_CC*1

real distance,search_step,dis1,dis2,dis3,lbond,charge(natom)

******************************************************************
*    define global variables
******************************************************************
common charge
   common xc,yc,zc


******************************************************************
*    Read the input file
******************************************************************

open(10,file='MCM41-final.car',status='old')
do 20 i=1,natom0
read(10,*)a(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i),
   &          atomname(i),charge(i)
20 continue
   close(10)

*******************************************************************
*    write the initial file to check whether the initial structure is read correctly.
*******************************************************************
open(30,file='output.car',access='append')
do 40 i=1,natom0
write(30,888)a(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i),
   &          atomname(i),charge(i)
40 continue
   close(30)

natom_add=0
Tatom=natom0+natom_add

* NSiT=0
* NSi_S=0
* NSi_C=0

* NCT=0
* NC_S=0
* NC_C=0

* NNT=0
* NN_S=0
* NN_C=0

do 45 itt=1,namino
*********************************************************
*    add Si to the chosen Oxygen atom
*********************************************************

   call HO_list(Ohydroxy,Ohydroxy_SN,kjishu)
   Nra=int(RAN2(IDUM)*kjishu)
Temp_SN=Ohydroxy_SN(Nra)
xtemp=Ohydroxy(Nra,2)
ytemp=Ohydroxy(Nra,3)
ztemp=Ohydroxy(Nra,4)
   lbond=1.90
call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN,
   &                xfinal,yfinal,zfinal)

   natom_add=natom_add+1
Tatom=natom0+natom_add

* NSiT=NSiT+1
* NSi_C=INT((NSiT+60)/99)
* NSi_S=NSiT+60-99*NSi_C
*
* if(NSi_C.eq.0)NSi_CC='T'
* if(NSi_C.eq.1)NSi_CC='U'
* if(NSi_C.eq.2)NSi_CC='V'
* if(NSi_C.eq.3)NSi_CC='W'
* if(NSi_C.eq.4)NSi_CC='X'
* if(NSi_C.eq.5)NSi_CC='Y'
* if(NSi_C.eq.6)NSi_CC='Z'

* a(Tatom)='Si'//'NSi_S'//'NSi_CC'
   a(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)

*    do 2060 i=1,nho-2
*
*    do 2065 ix=-2,kjishu+2
*       if(Ohydroxy(i,1).gt.templist(ix)-0.1.and.
*    &    Ohydroxy(i,1).lt.templist(ix)+0.1)then
* goto 2060
* endif
*2065 continue
*
*    do 2070 j=i+1,nho-1
*
*    do 2075 ix=-2,kjishu+2
*       if(Ohydroxy(j,1).gt.templist(ix)-0.1.and.
*    &       Ohydroxy(j,1).lt.templist(ix)+0.1)then
*       goto 2070
*       endif
*2075    continue
*
*       do 2080 k=j+1,nho
*
*       do 2085 ix=-2,kjishu+2
*          if(Ohydroxy(k,1).gt.templist(ix)-0.1.and.
*    &          Ohydroxy(k,1).lt.templist(ix)+0.1)then
*       goto 2080
*       endif
*2085    continue
*
*       dis1=sqrt((Ohydroxy(i,2)-Ohydroxy(j,2))**2+
*    &       (Ohydroxy(i,3)-Ohydroxy(j,3))**2+
*    &             (Ohydroxy(i,4)-Ohydroxy(j,4))**2)
*
*       dis2=sqrt((Ohydroxy(i,2)-Ohydroxy(k,2))**2+
*    &       (Ohydroxy(i,3)-Ohydroxy(k,3))**2+
*    &             (Ohydroxy(i,4)-Ohydroxy(k,4))**2)

*       dis3=sqrt((Ohydroxy(j,2)-Ohydroxy(k,2))**2+
*    &       (Ohydroxy(j,3)-Ohydroxy(k,3))**2+
*    &             (Ohydroxy(j,4)-Ohydroxy(k,4))**2)

*          if (dis1.gt.1.8.and.dis1.lt.2.6.and.
*    &       dis2.gt.1.8.and.dis1.lt.2.6.and.
*    &       dis3.gt.1.8.and.dis1.lt.2.6)then

*          kjishu=kjishu+1

*          do 2090 imm=1,4
*             list3(kjishu,imm)=Ohydroxy(i,imm)
*2090          continue
*             templist3(kjishu,1)=int(Ohydroxy(i,1))
*             charge(Ohydroxy_SN(i))=2.0

*          do 2100 imm=1,4
*             list3(kjishu,imm+4)=Ohydroxy(j,imm)
*2100          continue
*             templist3(kjishu,2)=int(Ohydroxy(j,1))
*             charge(Ohydroxy_SN(j))=2.0

*          do 2110 imm=1,4
*             list3(kjishu,imm+8)=Ohydroxy(k,imm)
*2110          continue
*             templist3(kjishu,3)=int(Ohydroxy(k,1))
*             charge(Ohydroxy_SN(k))=2.0

回复此楼

» 猜你喜欢

Cas 72-43-5需要30g，定制合成，能接单的留言已经有8人回复
求助:我三月中下旬出站，青基依托单位怎么办？已经有6人回复
北京211副教授，35岁，想重新出发，去国外做博后，怎么样？已经有8人回复
磺酰氟产物，毕不了业了！已经有5人回复
论文终于录用啦！满足毕业条件了已经有25人回复
2026年机械制造与材料应用国际会议（ICMMMA 2026）已经有3人回复
自荐读博已经有3人回复
不自信的我已经有5人回复
投稿Elsevier的杂志（返修），总是在选择OA和subscription界面被踢皮球已经有8人回复

» 本主题相关价值贴推荐，对您同样有帮助:

国家自然基金pdf格式的生成已经有12人回复
经典网格生成书籍【转载】已经有20人回复
求助：讨论不同压力下蒸汽消耗的关系已经有9人回复
ms5.5 在win8 总是出问题已经有13人回复
格式试剂怎么做GC或GC-MS？已经有10人回复
MCM-41表面接枝已经有7人回复
单氰胺在酸性条件下（浓盐酸）会不会生成氢氰酸已经有20人回复
不知道怎么回答审稿人的问题，求助大家,(关于分子筛MCM-41的）已经有10人回复
【求助成功】MS模拟MCM-41过程中forcite filed计算的问题已经有12人回复

好好学习，天天向上。

1楼 2010-05-31 22:19:03

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

zyj8119

木虫 (著名写手)

应助: 65 (初中生)
贵宾: 0.003
金币: 915.1
散金: 1440
红花: 35
帖子: 2936
在线: 1329.4小时
虫号: 664177
注册: 2008-11-29
性别: GG
专业: 理论和计算化学

*          goto 2061
*          endif

*2080    continue
*2070    continue
*2060  continue

*    open(2120,file='list3.car',access='append')
* do 2130 i=1,kjishu
*    write(2120,777)int(list3(i,1)),list3(i,2),list3(i,3),list3(i,4),
*    &             templist3(i,1)
*    write(2120,777)int(list3(i,5)),list3(i,6),list3(i,7),list3(i,8),
*    &             templist3(i,2)
*    write(2120,777)int(list3(i,9)),list3(i,10),list3(i,11),list3(i,12)
*    &             ,templist3(i,3)
* write(2120,*)''
*2130  continue
*777 Format(I4,3X,F12.9,2X,F13.9,3X,F12.9,3X,I4)
*    close(2120)

***************************************************
*    renew hydroxy oxygen list
***************************************************
*    do 2135 i=1,nho
* Ohydroxy_SN(i)=0
*    do 2136 j=1,4
*       Ohydroxy(i,j)=0
*2136    continue
*2135  continue

*    kjishu=0
*    do 2140 i=1,natom0
*    if (charge(i).eq.1.0) then
*    kjishu=kjishu+1
*    Ohydroxy_SN(kjishu)=i
*       Ohydroxy(kjishu,1)=kjishu
*    Ohydroxy(kjishu,2)=xc(i)
*    Ohydroxy(kjishu,3)=yc(i)
*    Ohydroxy(kjishu,4)=zc(i)
*    endif
*2140  continue

*    open(2150,file='hydrogen oxygen-1.car',access='append')
*    do 2160 i=1,nho
* write(2150,*)(Ohydroxy(i,j),j=1,4),Ohydroxy_SN(i)
*2160  continue
*    close(2150)

*    stop
************************************
* choose a initial oxygen atom randomly
************************************

*    ron=int(ran2(idum)*nho)
*    rox=Ohydroxy(nho,2)
*    roy=Ohydroxy(nho,3)
*    roz=Ohydroxy(nho,4)

*    nOtemp=0
*    do 60 i=1,nho
*    distance=sqrt((rox-Ohydroxy(i,2))**2+(roy-Ohydroxy(i,3))**2+
*    &    (roz-Ohydroxy(i,4))**2)
*       if(distance.gt.1.0.and.distance.lt.2.6)then
*          nOtemp=nOtemp+1
*       Otemp(nOtemp,1)=nOtemp
*       Otemp(nOtemp,2)=Ohydroxy(i,2)
*    Otemp(nOtemp,3)=Ohydroxy(i,3)
*    Otemp(nOtemp,4)=Ohydroxy(i,4)
*       endif
*60 continue

*    kstop=0
*    do 70 i=1,nOtemp-1
*    do 80 j=i+1,nOtemp
*    distance=sqrt((Ohydroxy(i,2)-Ohydroxy(j,2))**2+
*    &                (Ohydroxy(i,3)-Ohydroxy(j,3))**2+
*    &                (Ohydroxy(i,4)-Ohydroxy(j,4))**2)
*    if(distance.gt.1.0.and.distance.lt.2.6)then

*************************************************
*    Mark of interrupt service routine 2
*************************************************
*    kstop=1

*    OTa(1)=i
*    OTa(2)=Otemp(i,2)
*    OTa(3)=Otemp(i,3)
*    OTa(4)=Otemp(i,4)

*    OTb(1)=j
*    OTb(2)=Otemp(j,2)
*    OTb(3)=Otemp(j,3)
*    OTb(4)=Otemp(j,4)

*    goto 90
*       endif
*80    continue
*70 continue

********************************************
*    warning 1
********************************************
*    if(kstop.eq.0)then
* write(*,*)'Warning',kstop,': no tri-grafting any more'
* stop
* endif
********************************************
*    warning 1
********************************************

*    kstop=0
*90 search_step=0.01
*    do 100 i=-300,300
*       do 110 j=-300,300
*       do 120 k=-300,300
*       temp(1)=rox+search_step*i
*       temp(2)=roy+search_step*j
*       temp(3)=roz+search_step*k

*       dis1=sqrt((temp(1)-rox)**2+
*    &             (temp(2)-roy)**2+
*    &             (temp(3)-roz)**2)

*       dis2=sqrt((temp(1)-OTa(2))**2+
*    &             (temp(2)-OTa(3))**2+
*    &             (temp(3)-OTa(4))**2)

*       dis3=sqrt((temp(1)-OTb(2))**2+
*    &             (temp(2)-OTb(3))**2+
*    &             (temp(3)-OTb(4))**2)

*          if (dis1.gt.1.8.and.dis1.lt.2.0.and.
*    &       dis2.gt.1.8.and.dis1.lt.2.0.and.
*    &       dis3.gt.1.8.and.dis1.lt.2.0)then

*             do 130 im=1,natom+natom_add
*    distance=sqrt((temp(1)-xc(im))**2+
*    &                      (temp(2)-yc(im))**2+
*    &                      (temp(3)-zc(im))**2)
*                if(distance.le.3.0)then
* goto 120
* endif
*130          continue
*************************************************
*    Mark of interrupt service routine 2
*************************************************
*          kstop=2

*          natom_add=natom_add+1
*       Tatom=natom0+natom_add
*       a(Tatom)='Si'
*       xc(Tatom)=temp(1)
*       yc(Tatom)=temp(2)
*       zc(Tatom)=temp(3)
*       fx(Tatom)='XXXX'
*       occupation(Tatom)=1
*       fft(Tatom)='Si3'
*          atomname(Tatom)='Si'
*       charge(Tatom)=3.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)
      Temp_SN=Tatom
      call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN,
   &                xfinal,yfinal,zfinal)

         natom_add=natom_add+1
      Tatom=natom0+natom_add

*       NCT=NCT+1
*       NC_C=INT(NCT/99)
*       NC_S=NCT-99*NC_C

*       if(NC_C.eq.0)NC_CC=''
*       if(NC_C.eq.1)NC_CC='A'
*       if(NC_C.eq.2)NC_CC='B'
*       if(NC_C.eq.3)NC_CC='C'
*       if(NC_C.eq.4)NC_CC='D'
*       if(NC_C.eq.5)NC_CC='E'
*       if(NC_C.eq.6)NC_CC='F'

*       a(Tatom)='C'//'NC_S'//'NC_CC'
         a(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)
      Temp_SN=Tatom
      call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN,
   &                xfinal,yfinal,zfinal)

         natom_add=natom_add+1
      Tatom=natom0+natom_add

*       NCT=NCT+1
*       NC_C=INT(NCT/99)
*       NC_S=NCT-99*NC_C

*       if(NC_C.eq.0)NC_CC=''
*       if(NC_C.eq.1)NC_CC='A'
*       if(NC_C.eq.2)NC_CC='B'
*       if(NC_C.eq.3)NC_CC='C'
*       if(NC_C.eq.4)NC_CC='D'
*       if(NC_C.eq.5)NC_CC='E'
*       if(NC_C.eq.6)NC_CC='F'

*       a(Tatom)='C'//'NC_S'//'NC_CC'
         a(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)
      Temp_SN=Tatom
      call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN,
   &                xfinal,yfinal,zfinal)

         natom_add=natom_add+1
      Tatom=natom0+natom_add

*       NCT=NCT+1
*       NC_C=INT(NCT/99)
*       NC_S=NCT-99*NC_C

*       if(NC_C.eq.0)NC_CC=''
*       if(NC_C.eq.1)NC_CC='A'
*       if(NC_C.eq.2)NC_CC='B'
*       if(NC_C.eq.3)NC_CC='C'
*       if(NC_C.eq.4)NC_CC='D'
*       if(NC_C.eq.5)NC_CC='E'
*       if(NC_C.eq.6)NC_CC='F'

*       a(Tatom)='C'//'NC_S'//'NC_CC'
         a(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)

赞一下

回复此楼

好好学习，天天向上。

2楼2010-05-31 22:19:24

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

zyj8119

木虫 (著名写手)

应助: 65 (初中生)
贵宾: 0.003
金币: 915.1
散金: 1440
红花: 35
帖子: 2936
在线: 1329.4小时
虫号: 664177
注册: 2008-11-29
性别: GG
专业: 理论和计算化学

******************************************************************
*       add N atom on the chosen C3 atom
******************************************************************

      lbond=1.53
      xtemp=xc(Tatom)
      ytemp=yc(Tatom)
      ztemp=zc(Tatom)
      Temp_SN=Tatom
      call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN,
   &                xfinal,yfinal,zfinal)

         natom_add=natom_add+1
      Tatom=natom0+natom_add

*       NNT=NNT+1
*       NN_C=INT(NNT/99)
*       NN_S=NNT-99*NN_C

*       if(NN_C.eq.0)NN_CC=''
*       if(NN_C.eq.1)NN_CC='A'
*       if(NN_C.eq.2)NN_CC='B'
*       if(NN_C.eq.3)NN_CC='C'
*       if(NN_C.eq.4)NN_CC='D'
*       if(NN_C.eq.5)NN_CC='E'
*       if(NN_C.eq.6)NN_CC='F'

*       a(Tatom)='N'//'NN_S'//'NN_CC'
         a(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

****************************************************
*    rewrite the atom serial.
****************************************************
   open(46,file='Satom.car',status='old')
do 47 i=natom0+1,Tatom
read(46,1888)a(i)
47 continue
   close(46)
1888  format(A5)

*       endif
*
*120       continue
*110    continue
*100 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

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
*    To make a list of hydroxy oxygen.
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
   subroutine HO_list(Ohydroxy,Ohydroxy_SN,kjishu)

integer nho,natom,natom0
parameter (natom=15000,natom0=9992,nho=2319)
integer kjishu,Ohydroxy_SN(nho)
double precision Ohydroxy(nho,4)
*******************************************************************
*    declare the global varialbe
*******************************************************************
real charge(natom)
double precision xc(natom),yc(natom),zc(natom)
common charge
common xc,yc,zc

   kjishu=0
   do 50 i=1,natom0
   if (charge(i).eq.1.0) then
   kjishu=kjishu+1
   Ohydroxy_SN(kjishu)=i
      Ohydroxy(kjishu,1)=kjishu
   Ohydroxy(kjishu,2)=xc(i)
   Ohydroxy(kjishu,3)=yc(i)
   Ohydroxy(kjishu,4)=zc(i)
   endif
50 continue
******************************************************************
*    write out the hydroxy oxygen list
******************************************************************
   open(51,file='HO_list.car',access='append')
   do 52 i=1,kjishu
write(51,999)int(Ohydroxy(i,1)),Ohydroxy(i,2),
   &          Ohydroxy(i,3),Ohydroxy(i,4),Ohydroxy_SN(i)
52 continue
   write(51,*)'***********************************************'
write(51,*)''
   close(51)
999 Format(I4,3X,F12.9,2X,F13.9,3X,F12.9,3X,I4)
   endsubroutine
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
*    add new atom
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
subroutine addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN,
   &                xfinal,yfinal,zfinal)

integer natom
parameter(natom=15000)

real search_step,lbond,minlbond,maxlbond
double precision xtemp,ytemp,ztemp,targetx,targety,targetz
double precision dis1,dis2,xc(natom),yc(natom),zc(natom)
double precision xfinal,yfinal,zfinal
double precision center(14,2),xcenter,ycenter,lc,lct
integer ia,ib,ic,ja,jb,jc
integer Tatom,Temp_SN

*****************************************************************
*    difine global variable
*****************************************************************
   common xc,yc,zc
*****************************************************************
*****************************************************************

   open(3000,file='center.txt',status='old')
do 3010 i=1,14
read(3000,*)(center(i,j),j=1,2)
3010  continue
   close(3000)

*    open(3011,file='center-1.txt',access='append')
* do 3012 i=1,14
*    write(3011,*)(center(i,j),j=1,2)
*3012  continue
*    close(3011)

   lc=1000
xcenter=0
ycenter=0
do 3020 i=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)
endif
3020  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
         if(im.ne.Temp_SN)then
            dis2=sqrt((targetx-xc(im))**2+
   &             (targety-yc(im))**2+
   &             (targetz-zc(im))**2)
            if(dis2.lt.3.0)then
            goto 1020
            endif
         endif

*          if(im.ne.Temp_SN.and.dis2.lt.3.0)then
*          goto 1020
*          endif
1030       continue
      goto 1040
      endif
1020    continue
1010 continue
1000  continue

1040  xfinal=targetx
   yfinal=targety
zfinal=targetz

赞一下

回复此楼

好好学习，天天向上。

3楼2010-05-31 22:19:49

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

zyj8119

木虫 (著名写手)

应助: 65 (初中生)
贵宾: 0.003
金币: 915.1
散金: 1440
红花: 35
帖子: 2936
在线: 1329.4小时
虫号: 664177
注册: 2008-11-29
性别: GG
专业: 理论和计算化学

elseif(xtemp.gt.xcenter+5)then
   ia=-20
   ib=0
   ic=1

   minlbond=lbond-0.1
maxlbond=lbond+0.1

search_step=0.1
do 1001 i=ia,ib,ic
   do 1011 j=ja,jb,jc
      do 1021 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 1031 im=1,Tatom
         if(im.ne.Temp_SN)then
            dis2=sqrt((targetx-xc(im))**2+
   &             (targety-yc(im))**2+
   &             (targetz-zc(im))**2)
            if(dis2.lt.3.0)then
            goto 1021
            endif
         endif

*       if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then
*          do 1031 im=1,Tatom
*       dis2=sqrt((targetx-xc(im))**2+
*    &             (targety-yc(im))**2+
*    &             (targetz-zc(im))**2)
*          if(im.ne.Temp_SN.and.dis2.lt.3.0)then
*          goto 1021
*          endif
1031       continue
      goto 1041
      endif
1021    continue
1011 continue
1001  continue

1041  xfinal=targetx
   yfinal=targety
zfinal=targetz

   else
      ia=20
   ib=-20
   ic=-1

   minlbond=lbond-0.1
maxlbond=lbond+0.1

search_step=0.1
do 1002 j=ja,jb,jc
   do 1012 i=ia,ib,ic
      do 1022 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 1032 im=1,Tatom
         if(im.ne.Temp_SN)then
            dis2=sqrt((targetx-xc(im))**2+
   &             (targety-yc(im))**2+
   &             (targetz-zc(im))**2)
            if(dis2.lt.3.0)then
            goto 1022
            endif
         endif

*       if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then
*          do 1032 im=1,Tatom
*       dis2=sqrt((targetx-xc(im))**2+
*    &             (targety-yc(im))**2+
*    &             (targetz-zc(im))**2)
*          if(im.ne.Temp_SN.and.dis2.lt.3.0)then
*          goto 1022
*          endif

1032       continue
      goto 1042
      endif
1022    continue
1012 continue
1002  continue

1042  xfinal=targetx
   yfinal=targety
zfinal=targetz

endif

endsubroutine
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
*    Random Number Generator
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

                  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
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

赞一下

回复此楼

好好学习，天天向上。

4楼2010-05-31 22:20:06

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖