72e1a47eefd0052379c7bc9a2d991f2545834020
[unres.git] / read_ref_str.F
1       subroutine read_ref_structure(*)
2 C
3 C Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
4 C angles.
5 C
6       implicit none
7       include 'DIMENSIONS'
8       include 'DIMENSIONS.ZSCOPT'
9       include 'DIMENSIONS.FREE'
10       include 'DIMENSIONS.COMPAR'
11       include 'COMMON.IOUNITS'
12       include 'COMMON.GEO'
13       include 'COMMON.VAR'
14       include 'COMMON.INTERACT'
15       include 'COMMON.LOCAL'
16       include 'COMMON.NAMES'
17       include 'COMMON.CHAIN'
18       include 'COMMON.FFIELD'
19       include 'COMMON.SBRIDGE'
20       include 'COMMON.HEADER'
21       include 'COMMON.CONTROL'
22       include 'COMMON.CONTACTS1'
23       include 'COMMON.PEPTCONT'
24       include 'COMMON.TIME1'
25       include 'COMMON.COMPAR'
26       character*4 sequence(maxres)
27       integer rescode
28       double precision x(maxvar)
29       integer itype_pdb(maxres)
30       logical seq_comp
31       integer i,j,k,nres_pdb,iaux
32       double precision ddsc,dist
33       integer ilen
34       external ilen
35 C
36       nres0=nres
37       write (iout,*) "pdbref",pdbref
38       if (pdbref) then
39         read(inp,'(a)') pdbfile
40         write (iout,'(2a,1h.)') 'PDB data will be read from file ',
41      &    pdbfile(:ilen(pdbfile))
42         open(ipdbin,file=pdbfile,status='old',err=33)
43         goto 34 
44   33    write (iout,'(a)') 'Error opening PDB file.'
45         return1
46   34    continue
47         do i=1,nres
48           itype_pdb(i)=itype(i)
49         enddo
50         call readpdb(.true.)
51         do i=1,nres
52           iaux=itype_pdb(i)
53           itype_pdb(i)=itype(i)
54           itype(i)=iaux
55         enddo
56         close (ipdbin)
57         nres_pdb=nres
58         nres=nres0
59         nstart_seq=nnt
60         if (nsup.le.(nct-nnt+1)) then
61           do i=0,nct-nnt+1-nsup
62             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),
63      &        nsup)) then
64               do j=nnt+nsup-1,nnt,-1
65                 do k=1,3
66                   cref(k,nres+j+i)=cref(k,nres_pdb+j)
67                 enddo
68               enddo
69               do j=nnt+nsup-1,nnt,-1
70                 do k=1,3
71                   cref(k,j+i)=cref(k,j)
72                 enddo
73                 phi_ref(j+i)=phi_ref(j)
74                 theta_ref(j+i)=theta_ref(j)
75                 alph_ref(j+i)=alph_ref(j)
76                 omeg_ref(j+i)=omeg_ref(j)
77               enddo
78 #ifdef DEBUG
79               do j=nnt,nct
80                 write (iout,'(i5,3f10.5,5x,3f10.5)') 
81      &            j,(cref(k,j),k=1,3),(cref(k,j+nres),k=1,3)
82               enddo
83 #endif
84               nstart_seq=nnt+i
85               nstart_sup=nnt+i
86               goto 111
87             endif
88           enddo
89           write (iout,'(a)') 
90      &            'Error - sequences to be superposed do not match.'
91           return1
92         else
93           do i=0,nsup-(nct-nnt+1)
94             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),
95      &        nct-nnt+1)) 
96      &      then
97               nstart_sup=nstart_sup+i
98               nsup=nct-nnt+1
99               goto 111
100             endif
101           enddo 
102           write (iout,'(a)') 
103      &            'Error - sequences to be superposed do not match.'
104         endif
105   111   continue
106         write (iout,'(a,i5)') 
107      &   'Experimental structure begins at residue',nstart_seq
108       else
109         call read_angles(inp,*38)
110         goto 39
111    38   write (iout,'(a)') 'Error reading reference structure.'
112         return1
113    39   call chainbuild     
114         nstart_sup=nnt
115         nstart_seq=nnt
116         nsup=nct-nnt+1
117         do i=1,2*nres
118           do j=1,3
119             cref(j,i)=c(j,i)
120           enddo
121         enddo
122       endif
123       nend_sup=nstart_sup+nsup-1
124       do i=1,2*nres
125         do j=1,3
126           c(j,i)=cref(j,i)
127         enddo
128       enddo
129       do i=1,nres
130         do j=1,3
131           dc(j,nres+i)=cref(j,nres+i)-cref(j,i)
132         enddo
133         if (itype(i).ne.10) then
134           ddsc = dist(i,nres+i)
135           do j=1,3
136             dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
137           enddo
138         else
139           do j=1,3
140             dc_norm(j,nres+i)=0.0d0
141           enddo
142         endif
143 c        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
144 c         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
145 c         dc_norm(3,nres+i)**2
146         do j=1,3
147           dc(j,i)=c(j,i+1)-c(j,i)
148         enddo
149         ddsc = dist(i,i+1)
150         do j=1,3
151           dc_norm(j,i)=dc(j,i)/ddsc
152         enddo
153       enddo
154 c      print *,"Calling contact"
155       call contact(.true.,ncont_ref,icont_ref(1,1),
156      &  nstart_sup,nend_sup)
157 c      print *,"Calling elecont"
158       call elecont(.true.,ncont_pept_ref,
159      &   icont_pept_ref(1,1),
160      &   nstart_sup,nend_sup)
161        write (iout,'(a,i3,a,i3,a,i3,a)')
162      &    'Number of residues to be superposed:',nsup,
163      &    ' (from residue',nstart_sup,' to residue',
164      &    nend_sup,').'
165       return
166       end