[請益] 程式問題請教

看板Fortran作者mr24 (=大頭=)時間12年前 (2012/04/25 18:23)推噓2(2推 0噓 9→)

留言11則, 2人參與討論串1/1

先PO一下我原做這程式的動機我是有一個檔案他會記錄一個碳球由60顆碳所組成的空間座標檔他會隨著我跑模擬後而計錄每一步階的原子位置變化而我的目的是去計算這顆碳球的表面積，他是由20個六角形跟12個五角形所組成，我已經將角形是由哪幾顆原子所組成的紀錄好，只差帶入公式去計算，而我在這是採用海龍公式。這邊開始PO我的主要程式部分，我只先PO其中一小塊面積的計算部分順便說明參數定義 iTimeStep 表我總共跑了多少步階 iBins 表我有多少顆原子 n 是為了先給我的原子檔案上編號 x,y,z我是替每一顆原子所上的編號 open(33,file='input.dat',status='unknown') open(66,file='out.dat',status='unknown') open(67,file='area.dat',status='unknown') n=0 iTimeStep=0 321 format(1x, i3,1x, a1, 3(1x,e15.8)) 322 format(1x,i6,1x,e15.8) do while (1 .eq. 1) read(33,*,end=25) iBins read(33,*) iTimeStep=iTimeStep+1 write (66,*)iTimeStep write (66,*)iBins do i = 1, iBins read(33,*) c,x(i),y(i),z(i) n=n+1 write(66,321)n,c,x(i),y(i),z(i) !!aa1 if(i.eq.1)then!!!!!!!!!!!!!!!!!!!!!!!!!!!A x1=x(i) y1=y(i) z1=z(i) endif if(i.eq.2)then x2=x(i) y2=y(i) z2=z(i) endif if(i.eq.3)then x3=x(i) y3=y(i) z3=z(i) endif if(i.eq.4)then x4=x(i) y4=y(i) z4=z(i) endif if(i.eq.5)then x5=x(i) y5=y(i) z5=z(i) endif 裡面的 if 是希望能夠讀原子編號，i=1時讀編號一的原子座標，依此類推... 當讀完時候開始做裡面積算 cx=(x1+x2+x3+x4+x5)/5 cy=(y1+y2+y3+y4+y5)/5 cz=(z1+z2+z3+z4+z5)/5 !calculate the dis from center to each point lc1=(((cx-x1)**2)+((cy-y1)**2)+((cz-z1)**2))**0.5 lc2=(((cx-x2)**2)+((cy-y2)**2)+((cz-z2)**2))**0.5 lc3=(((cx-x3)**2)+((cy-y3)**2)+((cz-z3)**2))**0.5 lc4=(((cx-x4)**2)+((cy-y4)**2)+((cz-z4)**2))**0.5 lc5=(((cx-x5)**2)+((cy-y5)**2)+((cz-z5)**2))**0.5 !calculate the dis from point to point l12=(((x1-x2)**2)+((y1-y2)**2)+((z1-z2)**2))**0.5 l23=(((x2-x3)**2)+((y2-y3)**2)+((z2-z3)**2))**0.5 l34=(((x3-x4)**2)+((y3-y4)**2)+((z3-z4)**2))**0.5 l45=(((x4-x5)**2)+((y4-y5)**2)+((z4-z5)**2))**0.5 l51=(((x5-x1)**2)+((y5-y1)**2)+((z5-z1)**2))**0.5 s1=(lc1+lc2+l12)/2 s2=(lc2+lc3+l23)/2 s3=(lc3+lc4+l34)/2 s4=(lc4+lc5+l45)/2 s5=(lc5+lc1+l51)/2 !!!!!!!!!!!計算每一三角形面積 a1=(s1*(s1-lc1)*(s1-lc2)*(s1-l12))**0.5 a2=(s2*(s2-lc2)*(s2-lc3)*(s2-l23))**0.5 a3=(s3*(s3-lc3)*(s3-lc4)*(s3-l34))**0.5 a4=(s4*(s4-lc4)*(s4-lc5)*(s4-l45))**0.5 a5=(s5*(s5-lc5)*(s5-lc1)*(s5-l51))**0.5 !!!!!五邊形面積 aa1=a1+a2+a3+a4+a5 enddo t5a=aa1 totala=t5a write(67,322)iTimeStep, totala n=0 enddo 25 continue end =================================================================== 我每一個原子都給他不同變數，但計算部分例如cx,cy,cz,lc1,l23,s2,a1 等共同部分則沒有個別設立目前單獨計算一塊面積時跟我用excel(一開始不會寫只好用excel去計算)算時幾乎一樣，小數三位後才不同，但我開始計算多塊加總後卻差別越來越大有懷疑是否之後的計算面積是否真的有代入新座標，還是仍用第一組座標還有不知我共同計算部分的參數是否也該個別獨立宣告，也懷疑是否是這邊共同計算部分的問題。內容有點多，希望板上的高手能幫忙指點一下，最近很煩惱 >"< 萬分感謝! -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 140.116.75.123