将Gaussian从头算动力学轨迹转换为pdb轨迹的程序GauMDconv

将Gaussian从头算动力学轨迹转换为pdb轨迹的程序GauMDconv

文/Sobereva    2012-Feb-10


1.1版下载地址(含windows下可执行文件,源代码,示例):/usr/uploads/file/20150610/20150610023408_87558.rar

 

简介:
此程序用来将Gaussian的BOMD或ADMP任务产生的轨迹转化为通用的含有多帧的pdb格式,以便被其它动力学可视化程序如VMD所进一步分析。

首先在Gaussian中执行BOMD或ADMP任务,必须产生chk文件。然后将chk文件用formchk转化为fch格式。之后启动本程序,输入fch文件名,然后输出要输出的pdb文件名即可。压缩包中的ADMP_example.fch是示例文件。如果轨迹中有多条轨迹,即ntraj设为了大于1,则本程序还会提示用户输入要读入哪条轨迹。

本程序原代码如下:

module util
contains
subroutine loclabel(fileid,label,ifound,irewind)  !此子程序用来将当前读写位置挪到fileid号文件中存在label字段的行
integer fileid,error
integer,optional :: ifound,irewind
character*80 c80
CHARACTER(LEN=*) label
if ((.not.present(irewind)).or.(present(irewind).and.irewind==1)) rewind(fileid)
do while(.true.)
 read(fileid,"(a80)",iostat=ierror) c80
 if (index(c80,label)/=0) then
  backspace(fileid)
  if (present(ifound)) ifound=1 ! Found result
  return
 end if
 if (ierror/=0) exit
end do
if (present(ifound)) ifound=0
end subroutine

end module


!!!----------- GauMDconv: Convert Gaussian AIMD trajectory from .fch to .pdb file, Version 1.0
program GauMDconv
use util
implicit real*8 (a-h,o-z)
character*2 :: name2ind_upcase(109)=(/ "H ","HE", & !Same as name2ind, but all characters are upper case, to cater to .pdb file
"LI","BE","B ","C ","N ","O ","F ","NE", & !3~10
"NA","MG","AL","SI","P ","S ","CL","AR", & !11~18
"K ","CA","SC","TI","V ","CR","MN","FE","CO","NI","CU","ZN","GA","GE","AS","SE","BR","KR", & !19~36
"RB","SR","Y ","ZR","NB","MO","TC","RU","RH","PD","AG","CD","IN","SN","SB","TE","I ","XE", & !37~54
"CS","BA","LA","CE","PR","ND","PM","SM","EU","GD","TB","DY","HO","ER","TM","YB","LU", & !55~71
"HF","TA","W ","RE","OS","IR","PT","AU","HG","TL","PB","BI","PO","AT","RN", & !72~86
"FR","RA","AC","TH","PA","U ","NP","PU","AM","CM","BK","CF","ES","FM","MD","NO","LR", & !87~103
"RF","DB","SG","BH","HS","MT" /) !104~109
character infilename*200,outfilename*200,fchtitle*70,findtraj*27
logical alive
integer,allocatable :: atmidx(:)
character,allocatable :: atmname*2(:)
real*8,allocatable :: atmx(:,:),atmy(:,:),atmz(:,:)
real*8 :: b2a=0.529177249D0
write(*,"(a)") "GauMDconv: Convert Gaussian AIMD trajectory from .fch to .pdb file, Version 1.0"
write(*,"(a)") "Programmed by Sobereva, 2012-FEB-10"
write(*,*)
write(*,*) "Input the .fch filename to be loaded"
do while(.true.)
 read(*,"(a)") infilename
 inquire(file=infilename,exist=alive)
 if (alive) exit
 write(*,*) "Cannot found this file, input again"
end do

open(10,file=infilename,status="old")
read(10,"(a)") fchtitle !标题行将直接挪到输出的pdb文件开头
call loclabel(10,'Trajectory MaxStp')
read(10,"(49x,i12)") nstep
nstep=nstep-1   !.fch文件中的MaxStp比实际的帧数大1
call loclabel(10,'Trajectory MaxJob')   !.fch文件中可能含有多条轨迹
read(10,"(49x,i12)") ntraj
itraj=1
if (ntraj>1) then
 write(*,"(a,i8,a)") "Note: There are",ntraj," trajectories, read which one?"
 read(*,*) itraj
end if
write(findtraj,"(a,i8,a)") 'Traj num',itraj,' Geometries'

call loclabel(10,'Atomic numbers')
read(10,"(49x,i12)") ncenter
allocate(atmidx(ncenter),atmname(ncenter),atmx(nstep,ncenter),atmy(nstep,ncenter),atmz(nstep,ncenter))
read(10,"(6i12)") (atmidx(iatm),iatm=1,ncenter)  !读入原子编号
atmname=name2ind_upcase(atmidx)  !获得元素符号

call loclabel(10,findtraj)  !寻找指定的那条轨迹所在位置
read(10,*)
read(10,"(5(1PE16.8))") ((atmx(ifps,iatm),atmy(ifps,iatm),atmz(ifps,iatm),iatm=1,ncenter),ifps=1,nstep)
close(10)
atmx=atmx*b2a   !.fch文件中以bohr为单位记录长度信息,转换为pdb中用的埃
atmy=atmy*b2a
atmz=atmz*b2a

write(*,*) "Input the .pdb filename to be outputted"
read(*,"(a)") outfilename
open(10,file=outfilename,status="replace")
write(10,"('TITLE     ',a)") fchtitle
do ifps=1,nstep
 write(10,"('MODEL  ',4x,i4)") ifps
 do i=1,ncenter
  write(10,"(a6,i5,1x,a4,1x,a3, 1x,a1,i4,4x,3f8.3,2f6.2,10x,a2)") &
  "HETATM",i,' '//atmname(i)//' ',"MOL",'A',1,atmx(ifps,i),atmy(ifps,i),atmz(ifps,i),1.0,0.0,atmname(i)
 end do
 write(10,"('TER')")
 write(10,"('ENDMDL')")  !每一对儿MODEL和ENDMDL中夹着一帧的信息。对于VMD等支持多帧的可视化软件,这样的pdb文件就可以作为一个轨迹来载入了
end do
write(10,"('END')")
close(10)

write(*,*) "Exporting pdb file finished! Press ENTER to exit"
pause
end program

评论已关闭