我看還是貼在這邊比較方便,可以寫中文註解。

而且隨時想要就複製下來。

以下程式
--
program collecting_energy
! This program collectes the energies at each step in a trajectory.
character(16) :: filename,Ofilename
logical :: file_exist
integer :: Nstep,Nrow,i
character(16) :: tim(6000),ke(6000),pe(6000),C2(6000),C4(6000)
print *, "Please enter a [filename] for loading your data"
print *, " Warning : the length of the filenam must be less than 16"
read(*,*) filename
inquire (file=filename,exist=file_exist)
if (file_exist .eq. .true.) then
print *, "Please enter a [filename for output] file"
read(*,*) Ofilename
print *, "Enter the [maximum step No.] in your trajectory job"
print *, " Warning : the maximum step should be less than
2000"
read(*,*) Nstep
Nrow = Nstep*3
open (unit=7,file=filename,status='old',action='read')
open (unit=8,file=Ofilename,status='new')
do i = 1, Nrow
read(7,2,advance='yes',end=4)
tim(i),C2(i),ke(i),C4(i),pe(i)
if (mod(i,3) .eq. 1) then
write(8,3) tim(i),ke(i),pe(i)
end if
2 format (5a16)
3 format (3a16)
end do
else
print *, filename," doesn't exist"
go to 999
end if
4 print *, "The energies were written in : ",Ofilename
999 stop
end
--
這是一個可以剪輯從G03的fchk檔抓下來的BOMD在每一步的時間和動能位能的小執行檔。

最後修正版。

自動剪出軌跡每步的時間,動能及位能。

執行時會要求你輸入你從fchk剪輯出來能量部分的檔案的檔名,如果你輸入錯檔名會自動

告訴你檔案不存在。

當你放在執行檔同資料夾的能量列檔案存在時,會自動剪輯出時間、動能、位能。

read(7,2,advance='yes',end=4)該段語法由end修正回圈結束,同樣防呆以免誤植過大之

最大步數。

mongqiu 發表在 痞客邦 PIXNET 留言(0) 人氣()