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