以下程式 program geom_fo_step_structure ! This program cut the coordinate of each atom at some one step in a trajectory. character(16) :: filename,Ofilename logical :: file_exist integer :: Natom,Nstep,Nrow,step,i,j,k,l,m,n character(51) :: C(30000) character(3) :: atoms(50) 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 *, "Enter the [number of atoms] in your job" print *, " Warning : the number of atoms should be less than 50" read(*,*) Natom print *, "Enter the [maximum step No.] in your trajectory job" print *, " Warning : the product of atoms No. and steps No. should be less than 30000" read(*,*) Nstep Nrow = Natom*Nstep print *, "Enter the [atom symbols] in the label order" read(*,*) (atoms(i), i = 1, Natom) 998 print *, "Enter the [number of the step] what you want to load" read(*,*) step if (step .gt. Nstep) then print *, "The step No. what you want to load is larger than the maximum step in your job" go to 998 end if open (unit=9,status='scratch') write(9,"(i5,a4)") step,"step" rewind(9) read(9,*) Ofilename open (unit=7,file=filename,status='old',action='read') open (unit=8,file=Ofilename,status='new') do j = 1, Nrow read(7,2,advance='yes',end=3) C(j) 2 format(a51) 3 end do k = ((step-1)*Natom) + 1 l = step*Natom do m = k, l if (mod(m,Natom) .ne. 0) then n = mod(m,Natom) else n = Natom end if write(8,4,advance='yes') atoms(n),C(m) 4 format(a3,a51) end do print *, "The coordinates what you want were written in : ",Ofilename else print *, filename," doesn't exist" go to 999 end if 999 stop end -- 最後一個剪輯程式,功能是剪輯前一個程式製造出來的座標檔案,對其剪輯出所指定的步 數,並在所有座標前端加上原子符號。自動命檔名為"n"step。 行號998為防呆,避免指定步數超過最大步數。
創作者介紹
創作者 就是部落格 的頭像
mongqiu

就是部落格

mongqiu 發表在 痞客邦 留言(1) 人氣( 100 )