Added putting secondary-structure information to the pdbfiles generated by xdrf2pdb...
[unres.git] / source / xdrfpdb / src-M / xdrf2pdb1.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,1000)
7       real*4 prec,potE,efree,rmsdev
8       real*8 etot
9       character*80 arg,seqfile,pdbfile
10       character*3 sequenc(maxres)
11       character*50 tytul
12       character*8 onethree,cfreq
13       character*8 ucase
14       external ucase
15       logical oneletter
16       integer rescode
17       external rescode
18       logical ilen,iblnk
19       
20       ifreq=1
21       if (iargc().lt.4) then
22         print '(a)',
23      &    "Usage: xdrf2pdb one/three seqfile cxfile conf [pdbfile]"
24         stop
25       endif
26       call getarg(1,onethree)
27       onethree = ucase(onethree)
28       if (onethree.eq.'ONE') then
29         oneletter = .true.
30       else if (onethree.eq.'THREE') then
31         oneletter = .false.
32       else
33         print *,"ONE or THREE must be specified"
34       endif
35       call getarg(2,seqfile)
36       open (1,file=seqfile,status='old')
37       if (oneletter) then
38         read(1,'(80a1)',end=10,err=10) (sequenc(i)(1:1),i=1,maxres)
39    10   continue
40         nres=i
41         i=0
42         do while (.not.iblnk(sequenc(i+1)(1:1)))
43           i=i+1
44         enddo 
45         nres=i
46         do i=1,nres
47           itype(i)=rescode(i,sequenc(i),1)
48         enddo
49       else
50         read(1,'(20(a3,1x))',end=11,err=11) (sequenc(i),i=1,maxres)
51    11   continue
52         nres=i
53         print *,"nres",nres
54         do i=1,nres
55           print *,i," ",sequenc(i)
56         enddo
57         print *
58         i=0
59         do while (.not.iblnk(sequenc(i+1)(1:1)))
60           print *,i+1," ",sequenc(i+1)," ",sequenc(i+1)(1:1)
61           i=i+1
62         enddo 
63         nres=i
64         print *,"nres",nres
65         do i=1,nres
66           itype(i)=rescode(i,sequenc(i),0)
67         enddo
68         print *,(itype(i),i=1,nres)
69       endif
70       call getarg(3,arg)
71       iext = index(arg,'.cx') - 1
72       if (iext.lt.0) then
73         print *,"Error - not a cx file"
74         stop
75       endif
76       print *,"arg ",arg
77       call getarg(4,cfreq)
78       read (cfreq,*) iconf
79       print *,"iconf",iconf
80       if (iargc().gt.4) then
81         call getarg(5,pdbfile)
82       else
83         pdbfile=arg(:iext)//'.pdb'
84       endif
85       open(9,file=pdbfile)
86       nnt = 1
87       if (itype(1).eq.21) nnt = 2
88       nct=nres
89       if (itype(nres).eq.21) nct = nres-1
90 c      if (nct.eq.nres-1) nres=nres-1
91 c      if (nnt.eq.2) nres=nres-1
92
93       print *,"nres",nres," nnt",nnt," nct",nct
94
95       call xdrfopen(ixdrf,arg, "r", iret)
96       print *,"iret",iret 
97       kk = 0
98       do while(.true.) 
99        prec=10000.0
100        isize=0
101        call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
102        call xdrfint(ixdrf, nss, iret) 
103        do j=1,nss
104         call xdrfint(ixdrf, ihpb(j), iret)
105         call xdrfint(ixdrf, jhpb(j), iret)
106        enddo
107        call xdrffloat(ixdrf, potE, iret)
108        if(iret.eq.0) exit
109        kk = kk + 1
110        call xdrffloat(ixdrf, efree, iret)
111        call xdrffloat(ixdrf, rmsdev, iret)
112        call xdrfint(ixdrf, iscor, iret)
113        if (kk.eq.iconf) then
114          print *,"pote",pote," efree",efree," rmsdev",rmsdev
115
116          print *,"isize",isize
117
118          if (isize .ne. nres+nct-nnt+1) then
119            print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
120          endif
121          do i=1,nres
122            do j=1,3
123              c(j,i)=coord(j,i)
124            enddo
125          enddo
126          ii = 0
127          do i=nnt,nct
128            ii = ii + 1 
129            do j=1,3
130              c(j,i+nres)=coord(j,ii+nres)
131            enddo
132          enddo
133          etot=potE
134          write (tytul,'(a,i6)') "Structure",kk
135          call pdbout(etot,tytul,9)
136          stop
137        endif
138       enddo
139      
140       end