Merge branch 'master' of mmka.chem.univ.gda.pl:unres
[unres.git] / source / xdrfpdb / src / xdrf2pdb-m.F
1       implicit real*8 (a-h,o-z)
2       include 'DIMENSIONS'
3       include 'COMMON.CHAIN'
4       include 'COMMON.INTERACT'
5       include 'COMMON.SBRIDGE'
6       real*4 coord(3,2*maxres)
7       real*4 prec,time,potE,uconst,t_bath,qfrag(100)
8       real*8 etot
9       character*80 arg,seqfile,pdbfile
10       character*3 sequenc(maxres)
11       character*50 tytul
12       character*8 onethree,cfreq,cntraj,citraj
13       character*3 licz
14       character*8 ucase
15       external ucase
16       logical oneletter,iblnk
17       integer rescode
18       external rescode
19       
20       ifreq=1
21       if (iargc().lt.5) then
22         print '(2a)',
23      &   "Usage: xdrf2pdb-m one/three seqfile cxfile ntraj itraj",
24      &   " [pdbfile] [freq]"
25         stop
26       endif
27       call getarg(1,onethree)
28       onethree = ucase(onethree)
29       if (onethree.eq.'ONE') then
30         oneletter = .true.
31       else if (onethree.eq.'THREE') then
32         oneletter = .false.
33       else
34         print *,"ONE or THREE must be specified"
35       endif
36       call getarg(2,seqfile)
37       open (1,file=seqfile,status='old')
38       if (oneletter) then
39         read(1,'(80a1)',end=10,err=10) (sequenc(i)(1:1),i=1,maxres)
40    10   continue
41         nres=i
42         i=0
43         do while (.not.iblnk(sequenc(i+1)(1:1)))
44 c        do while (.not.(iblnk(sequenc(i+1)(1:1)) == 0))
45           i=i+1
46         enddo 
47         nres=i
48         do i=1,nres
49           itype(i)=rescode(i,sequenc(i),1)
50         enddo
51       else
52         read(1,'(20(a3,1x))',end=11,err=11) (sequenc(i),i=1,maxres)
53    11   continue
54         nres=i
55         i=0
56         do while (.not.iblnk(sequenc(i+1)(1:1)))
57 c        do while (.not.(iblnk(sequenc(i+1)(1:1)) == 0))
58           i=i+1
59         enddo 
60         nres=i
61         do i=1,nres
62           itype(i)=rescode(i,sequenc(i),0)
63         enddo
64         print *,nres
65         print '(20(a3,1x))',(sequenc(i),i=1,nres)
66       endif
67       call getarg(3,arg)
68       iext = index(arg,'.cx') - 1
69       if (iext.lt.0) then
70         print *,"Error - not a cx file"
71         stop
72       endif
73       call getarg(4,cntraj)
74       read (cntraj,*) ntraj
75       call getarg(5,citraj)
76       read (citraj,*) itraj
77       if (iargc().gt.5) then
78         call getarg(6,pdbfile)
79       else
80         write(licz,'(bz,i3.3)') itraj
81         pdbfile=arg(:iext)//'_'//licz//'.pdb'
82       endif
83       if (iargc().gt.6) then
84         call getarg(7,cfreq)
85         read (cfreq,*) ifreq
86       endif
87 c      print *,"ifreq",ifreq," ntraj",ntraj," itraj",itraj
88       open(9,file=pdbfile)
89       nnt = 1
90       if (itype(1).eq.21) nnt = 2
91       nct=nres
92       if (itype(nres).eq.21) nct = nres-1
93       print *,"nnt",nnt," nct",nct
94       call xdrfopen(ixdrf,arg, "r", iret)
95       kk = 0
96       do while(.true.) 
97        call xdrffloat(ixdrf, time, iret)
98        print *,"time",time," iret",iret
99        if(iret.eq.0) exit
100        kk = kk + 1
101        call xdrffloat(ixdrf, potE, iret)
102        call xdrffloat(ixdrf, uconst, iret)
103        call xdrffloat(ixdrf, t_bath, iret)
104        print *,"potE",potE," uconst",uconst," t_bath",t_bath
105 c
106 c uconst_back not written in MREMD.F code
107 c#ifdef NEWUNRES
108 c       call xdrffloat(ixdrf, uconst_back, iret)
109 c       print *,"uconst_back",uconst_back
110 c#endif
111        call xdrfint(ixdrf, nss, iret) 
112        do j=1,nss
113         call xdrfint(ixdrf, ihpb(j), iret)
114         call xdrfint(ixdrf, jhpb(j), iret)
115        enddo
116        print *,"nss",nss
117        print *,"   ",(ihpb(j),"-",jhpb(j),j=1,nss)
118        call xdrfint(ixdrf, nfrag, iret)
119        call xdrfint(ixdrf, iset, iret)
120        do i=1,nfrag
121         call xdrffloat(ixdrf, qfrag(i), iret)
122        enddo
123        print *,"nfrag",nfrag,"iset",iset
124        prec=10000.0
125
126        isize=0
127        print *," call xdrf3coord"
128        call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
129        print *,"isize",isize
130
131 c       write (*,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
132 c       write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
133 c       write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
134 c       write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize)
135        if (mod(kk/ntraj,ifreq).eq.0 .and. mod(kk,ntraj).eq.itraj) then
136          if (isize .ne. nres+nct-nnt+1) then
137            print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
138          endif
139          do i=1,nres
140            do j=1,3
141              c(j,i)=coord(j,i)
142            enddo
143          enddo
144          ii = 0
145          do i=nnt,nct
146            ii = ii + 1 
147            do j=1,3
148              c(j,i+nres)=coord(j,ii+nres)
149            enddo
150          enddo
151          etot=potE
152          write (tytul,'(a,i6)') "Structure",kk
153          call pdbout(etot,tytul,9)
154        endif
155       enddo
156      
157       end