[gmx-users] ego2xtc
Bert de Groot
bgroot at gwdg.de
Thu Nov 28 10:49:38 CET 2002
Hi,
for those that might be interested to analyse EGO trajectories with the gromacs
analysis tools (or save disk space), I attached a small tool, ego2xtc, that
converts EGO output (a series of .ego) files to a .xtc trajectory.
I attached both source and binaries for linux, sgi and dec.
cheers,
Bert
____________________________________________________________________________
Dr. Bert de Groot
Max Planck Institute for Biophysical Chemistry
Theoretical molecular biophysics group
Am Fassberg 11
37077 Goettingen, Germany
tel: +49-551-2011306, fax: +49-551-2011089
email: bgroot at gwdg.de
http://www.mpibpc.gwdg.de/abteilungen/071/bgroot
____________________________________________________________________________
-------------- next part --------------
program ego2xtc
c23456789012345678901234567890123456789012345678901234567890123456789012
implicit none
integer i,j,id1,nat,step,ret,maxat,nargs,iargc,nfiles,egonr,
& gvflen,len,dlen,skip
parameter(maxat=1000000)
character*80 egobase,segonr,snul
character*255 buf,xtcnam,egonam,dirnam
real*4 x(maxat*3),box(9),dt,time,prec,minx,miny,minz,maxx,maxy,
& maxz
real*8 xx(maxat*3)
logical gvfisf,append,bego,bnr
prec=1000.0
xtcnam='ego.xtc'
append=.TRUE.
bego=.FALSE.
bnr=.FALSE.
do i=1,9
box(i)=0.0
enddo
snul='00000000'
dirnam='./'
skip=1
c----
c---- parse command line
c----
nargs=iargc()
if (nargs.eq.0) call usage()
j=0
10 continue
j=j+1
call getarg(j,buf)
call gvstcp(buf)
if (buf(1:1).ne.'-') then
call usage()
else if (buf(1:2).eq.'-x') then
if (j.ge.nargs) call usage()
j=j+1
call getarg(j,xtcnam)
else if (buf(1:2).eq.'-e') then
if (j.ge.nargs) call usage()
j=j+1
bego=.TRUE.
call getarg(j,egonam)
else if (buf(1:2).eq.'-s') then
if (j.ge.nargs) call usage()
j=j+1
call getarg(j,buf)
read(buf,*) skip
else if (buf(1:2).eq.'-n') then
if (j.ge.nargs) call usage()
j=j+1
bnr=.TRUE.
call getarg(j,buf)
read(buf,*) nfiles
else if (buf(1:2).eq.'-p') then
if (j.ge.nargs) call usage()
j=j+1
call getarg(j,buf)
read(buf,*) prec
else if (buf(1:2).eq.'-o') then
append=.FALSE.
else if (buf(1:2).eq.'-h') then
call usage()
else
call gvsshl(buf)
write (0,'(a,a)') 'eco2xtc: unknown command line argument:',
& buf(1:gvflen(buf))
write (0,*)
call usage()
endif
if (j.lt.nargs) goto 10
if (.not.bego.or..not.bnr) call usage()
c----
c---- open xtc file
c----
call gvsshl(xtcnam)
if (gvfisf(xtcnam)) then
if (append) then
write(6,'(a,a)') 'ego2xtc: appending to ',
& xtcnam(1:gvflen(xtcnam))
call xdrfopen(id1,xtcnam,"a",ret)
else
call backitup(xtcnam)
call xdrfopen(id1,xtcnam,"w",ret)
endif
else
call xdrfopen(id1,xtcnam,"w",ret)
endif
if (ret.eq.0) then
write(0,'(a,a)') 'ego2xtc: error opening xtc file:',
& xtcnam(1:gvflen(xtcnam))
stop
endif
c----
c---- retrieve the number associated with first ego file
c----
call bbopenr(1,egonam)
close(1)
call gvsshl(egonam)
len=gvflen(egonam)
egobase=egonam(len-11:len-4)
read(egobase,*) egonr
if (len.gt.12) then
dirnam=egonam(1:len-12)
call gvsshl(dirnam)
dlen=gvflen(dirnam)
endif
c----
c---- loop over all ego files
c----
do i=1,nfiles
write(segonr,*) egonr
call gvsshl(segonr)
len=gvflen(segonr)
write(buf,*) dirnam(1:dlen),snul(1:8-len),
& segonr(1:len),'.ego'
call gvsshl(buf)
len=gvflen(buf)
call read_ego(buf,nat,step,time,xx)
do j=1,3*nat
x(j)=real(xx(j)/10.0)
enddo
minx=10000.0
miny=10000.0
minz=10000.0
maxx=-10000.0
maxy=-10000.0
maxz=-10000.0
do j=1,nat
minx=min(minx,x(3*j-2))
miny=min(miny,x(3*j-1))
minz=min(minz,x(3*j))
maxx=max(maxx,x(3*j-2))
maxy=max(maxy,x(3*j-1))
maxz=max(maxz,x(3*j))
enddo
box(1)=maxx-minx
box(5)=maxy-miny
box(9)=maxz-minz
call writextc(id1,nat,step,time,box,x,prec,ret)
egonr=egonr+skip
enddo
call xdrfclose(id1,ret)
end
subroutine usage()
implicit none
write(0,100)
write(0,110)
write(0,120)
write(0,130)
write(0,140)
write(0,150)
stop
100 format('Usage: ego2xtc -e <first_ego_file>')
110 format(' -n nr_of_ego_files')
120 format(' [-x <xtc_file>] (default ego.xtc)')
130 format(' [-s step] (default 1)')
140 format(' [-p precision] (default 1000)')
150 format(' [-o] (overwrite xtc, default: append)')
end
subroutine backitup(filnm)
C----
C---- backs up filnm fo filnm~
C----
implicit none
integer len,gvflen
character*255 syslin
character*(*) filnm
len=gvflen(filnm)
write(6,100) filnm(1:len),filnm(1:len)
if (2*len+5.gt.80) then
write(0,300)
stop
endif
write(syslin,200) filnm(1:len),filnm(1:len)
len=gvflen(syslin)
call system(syslin(1:len))
100 format('ego2xtc: backing up ',a,' to ',a,'~')
200 format('mv ',a,' ',a,'~')
300 format('ego2xtc: buffer size exceeded in subroutine backitup')
return
end
subroutine bbopenr(unit,filnm)
C----
C---- opens file in read mode in current directory
C----
implicit none
integer unit,len,gvflen
character*(*) filnm
logical gvfisf
call gvsshl(filnm)
if (gvfisf(filnm)) then
open(unit,file=filnm,status='OLD')
else
len=gvflen(filnm)
write(0,100) filnm(1:len)
stop
endif
100 format('ego2xtc: error - file not found: ',a)
return
end
subroutine read_ego(egonam,nat,step,time,x)
implicit none
integer i,nat,step,line,gvflen,len
real*8 x(*),timestep,freq
real*4 time
character*255 buf
character*(*) egonam
logical gvfisf
nat=0
line=1
len=gvflen(egonam)
call bbopenr(1,egonam)
write(6,'(a,a)') 'ego2xtc: processing ',egonam(1:len)
10 continue
read(1,100,err=60,end=70) buf
line=line+1
if (buf(1:12).ne.'[BEGINCOORD]') goto 10
read(1,*) nat,step,freq,timestep
time=step*timestep/1000.0
do i=1,nat
read(1,*,err=60,end=70) x(3*i-2),x(3*i-1),x(3*i)
line=line+1
enddo
goto 80
60 write(0,'(a,a)') 'eco2xtc: error while reading: ',egonam(1:len)
write(0,'(a,i5)') 'line:',line
stop
70 write(0,'(a,a)') 'ego2xtc: end of file while reading: ',
& egonam(1:len)
write(0,'(a,i5)') 'line:',line
stop
80 continue
c write(0,*) 'read_ego finished, nat=',nat
close(1)
100 format(a80)
end
SUBROUTINE GVSTCP (TEXT)
C///
C-------------------------------------------------------------------------------
C---- ----
C---- SHIFTS ALL CHARACTERS LEFTWARDS IN TEXT UNTIL NO BLANKS ARE LEFT. ----
C---- ----
C---- THIS FUNCTION IS WRITTEN IN FORTRAN 77. ----
C---- ----
C-------------------------------------------------------------------------------
C\\\
IMPLICIT NONE
CHARACTER*(*) TEXT
INTEGER GVFLEN, LENTXT, ICHR, ILOOP
C----
C---- SKIP EMPTY TEXT
C----
IF (GVFLEN(TEXT).LE.1) RETURN
C----
C---- SET THE START LENGTH
C----
LENTXT=GVFLEN(TEXT)
C----
C---- CHECK IF THERE ARE ANY CHARACTERS
C----
IF (LENTXT.LT.2) RETURN
C----
C---- START WITH A LEFT SHIFT
C----
CALL GVSSHL (TEXT)
C----
C---- LOOP OVER TEXT AND FIND BLANKS
C----
ICHR=1
10 LENTXT=GVFLEN(TEXT)
IF (LENTXT.LT.2) RETURN
C----
C---- LOOP OVER THE STRING TILL FIRST BLANK IS FOUND, THEN SHIFT DOWN REST
C----
DO 30 ICHR=1,LENTXT-1
IF (TEXT(ICHR:ICHR).EQ.' ') THEN
C----------
C---------- SHIFT ALL HIGHER CHARS ONE DOWN, AND REDO THE TEST
C----------
DO 20 ILOOP=ICHR,LENTXT-1
TEXT(ILOOP:ILOOP)=TEXT(ILOOP+1:ILOOP+1)
20 CONTINUE
TEXT(LENTXT:LENTXT)=' '
GOTO 10
END IF
30 CONTINUE
END
INTEGER FUNCTION GVFLEN(TEXT)
C+++
C-------------------------------------------------------------------------------
C---- ----
C---- RETURNS THE LENGTH OF TEXT ----
C---- THE LENGTH IS DEFINED AS THE SEQUENCE NUMBER IN TEXT OF THE LAST ----
C---- NON-BLANK CHARACTER. UNDEFINED CHARACTERS ARE ASSUMED TO BE BLANK. ----
C---- ----
C---- TEXT IS A CHARACTER STRING OF ANY LENGTH ----
C---- DO NOT FORGET TO DECLARE GVFLEN AS AN INTEGER IN THE CALLING ROUTINE ----
C---- ----
C-------------------------------------------------------------------------------
C===
IMPLICIT NONE
CHARACTER*(*) TEXT
INTEGER I, ICH
C----
C---- INITIALIZE
C----
GVFLEN=0
C----
C---- MEASURE THE LENGTH
C----
DO 10 I=LEN(TEXT),1,-1
ICH=ICHAR(TEXT(I:I))
IF (ICH.NE.ICHAR(' ').AND.ICH.NE.0) THEN
GVFLEN=I
GOTO 998
END IF
10 CONTINUE
998 CONTINUE
RETURN
END
SUBROUTINE GVSSHR (TEXT,LENGTH)
C+++
C-------------------------------------------------------------------------------
C---- ----
C---- A RIGHT SHIFT IS PERFORMED ON TEXT UNTIL ALL TRAILING BLANKS ARE GONE ----
C---- TEXT SHOULD BE A CHARACTER STRING OF LENGTH LENGTH. ----
C---- ----
C---- STRANGE ERRORS WILL OCCUR IF TEXT IS NOT PROPERLY DEFINED IN THE ----
C---- CALLING MODULE. CHANCES ARE THAT THESE ERRORS SHOW UP IN A TOTALLY ----
C---- DIFFERENT WAY THEN YOU WOULD EXPECT. ----
C---- ----
C-------------------------------------------------------------------------------
C===
IMPLICIT NONE
INTEGER LENGTH, I
CHARACTER*(*) TEXT
C----
C---- LOOK FOR TRIVIAL ERRORS
C----
IF (LENGTH.LT.0) THEN
WRITE(0,*) ('Length of input text is < zero in GVSSHR')
RETURN
ELSE IF (LENGTH.EQ.0) THEN
WRITE(0,*) ('Length of input text is zero in GVSSHR')
RETURN
ELSE IF (LENGTH.EQ.1) THEN
WRITE(0,*) ('Length of input text is 1 in GVSSHR')
RETURN
END IF
C----
C---- IF THE TEXT IS OK, SEE IF IT ENDS AT BLANK
C----
10 IF (TEXT(LENGTH:LENGTH).NE.' ') RETURN
C----
C---- IF THERE ARE TRAILING BLANKS, SHIFT EVERYTHING UP BY ONE, AND TEST AGAIN
C----
DO 20 I=LENGTH,2,-1
TEXT(I:I)=TEXT(I-1:I-1)
20 CONTINUE
TEXT(1:1)=' '
GOTO 10
END
SUBROUTINE GVSSHL (TEXT)
C///
C-------------------------------------------------------------------------------
C---- ----
C---- A LEFT SHIFT IS PERFORMED ON TEXT UNTIL ALL LEADING BLANCKS ARE GONE ----
C---- TEXT SHOULD BE A CHARACTER STRING OF UNLIMITTED LENGTH ----
C---- ----
C---- THIS ROUTINE IS WRITTEN IN FORTRAN 77. ----
C---- ----
C-------------------------------------------------------------------------------
C\\\
IMPLICIT NONE
CHARACTER*(*) TEXT
INTEGER GVFLEN, I, J, LENTXT
C----
C---- DETERMINE THE LENGTH OF THE TEXT
C----
LENTXT=GVFLEN(TEXT)
IF (LENTXT.GT.0) THEN
C-------
C------- IF THERE ARE LEADING BLANKS, SHIFT TO THE LEFT
C-------
DO 30 I=1,LENTXT
IF (ICHAR(TEXT(I:I)).NE.ICHAR(' ')) THEN
IF (I.EQ.1) RETURN
DO 10 J=I,LENTXT
TEXT(J-I+1:J-I+1)=TEXT(J:J)
10 CONTINUE
DO 20 J=LENTXT-I+2,LENTXT
TEXT(J:J)=' '
20 CONTINUE
RETURN
END IF
30 CONTINUE
END IF
RETURN
END
LOGICAL FUNCTION GVFISF(FILNAM)
C///
C-------------------------------------------------------------------------------
C---- ----
C---- BECOMES TRUE IF THE FILE FILNAM EXISTS. ----
C---- ----
C---- THIS FUNCTION IS WRITTEN IN FORTRAN 77. ----
C---- ----
C-------------------------------------------------------------------------------
C\\\
IMPLICIT NONE
CHARACTER*(*) FILNAM
INTEGER GVFLEN
C-------
C------- TEST FOR EMPTY FILE NAME
C-------
IF (GVFLEN(FILNAM).EQ.0) THEN
GVFISF=.FALSE.
RETURN
END IF
C-------
C------- SEE IF THE FILE EXISTS
C-------
INQUIRE (FILE=FILNAM,EXIST=GVFISF)
RETURN
END
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ego2xtc.dec
Type: application/octet-stream
Size: 73728 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20021128/19e3038b/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ego2xtc.lnx
Type: application/octet-stream
Size: 106550 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20021128/19e3038b/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ego2xtc.sgi
Type: image/x-sgi-rgba
Size: 55412 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20021128/19e3038b/attachment.bin>
More information about the gromacs.org_gmx-users
mailing list