xdrf added to cluster
[unres.git] / source / cluster / wham / src-M / readpdb.f
1       subroutine readpdb
2 C Read the PDB file and convert the peptide geometry into virtual-chain 
3 C geometry.
4       implicit real*8 (a-h,o-z)
5       include 'DIMENSIONS'
6       include 'COMMON.LOCAL'
7       include 'COMMON.VAR'
8       include 'COMMON.CHAIN'
9       include 'COMMON.INTERACT'
10       include 'COMMON.IOUNITS'
11       include 'COMMON.GEO'
12       include 'COMMON.NAMES'
13       character*3 seq,atom,res
14       character*80 card
15       dimension sccor(3,20)
16       integer rescode
17       call permut(symetr)
18       ibeg=1
19       do
20         read (ipdbin,'(a80)',end=10) card
21         if (card(:3).eq.'END') then
22           goto 10
23         else if (card(:3).eq.'TER') then
24 C End current chain
25           ires_old=ires+1 
26           itype(ires_old)=21
27           ibeg=2
28 c          write (iout,*) "Chain ended",ires,ishift,ires_old
29           call sccenter(ires,iii,sccor)
30         endif
31 C Fish out the ATOM cards.
32         if (index(card(1:4),'ATOM').gt.0) then  
33           read (card(14:16),'(a3)') atom
34           if (atom.eq.'CA' .or. atom.eq.'CH3') then
35 C Calculate the CM of the preceding residue.
36             if (ibeg.eq.0) then
37               call sccenter(ires,iii,sccor)
38             endif
39 C Start new residue.
40 c            write (iout,'(a80)') card
41             read (card(24:26),*) ires
42             read (card(18:20),'(a3)') res
43             if (ibeg.eq.1) then
44               ishift=ires-1
45               if (res.ne.'GLY' .and. res.ne. 'ACE') then
46                 ishift=ishift-1
47                 itype(1)=21
48               endif
49 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
50               ibeg=0          
51             else if (ibeg.eq.2) then
52 c Start a new chain
53               ishift=-ires_old+ires-1
54 c              write (iout,*) "New chain started",ires,ishift
55               ibeg=0
56             endif
57             ires=ires-ishift
58 c            write (2,*) "ires",ires," ishift",ishift
59             if (res.eq.'ACE') then
60               ity=10
61             else
62               itype(ires)=rescode(ires,res,0)
63             endif
64             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
65             write (iout,'(2i3,2x,a,3f8.3)') 
66      &       ires,itype(ires),res,(c(j,ires),j=1,3)
67             iii=1
68             do j=1,3
69               sccor(j,iii)=c(j,ires)
70             enddo
71           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and. 
72      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
73             iii=iii+1
74             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
75           endif
76         endif
77       enddo
78    10 write (iout,'(a,i5)') ' Nres: ',ires
79 C Calculate dummy residue coordinates inside the "chain" of a multichain
80 C system
81       nres=ires
82       do i=2,nres-1
83 c        write (iout,*) i,itype(i)
84         if (itype(i).eq.21) then
85 c          write (iout,*) "dummy",i,itype(i)
86           do j=1,3
87             c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
88 c            c(j,i)=(c(j,i-1)+c(j,i+1))/2
89             dc(j,i)=c(j,i)
90           enddo
91         endif
92       enddo
93 C Calculate the CM of the last side chain.
94       call sccenter(ires,iii,sccor)
95       nsup=nres
96       nstart_sup=1
97       if (itype(nres).ne.10) then
98         nres=nres+1
99         itype(nres)=21
100         do j=1,3
101           dcj=c(j,nres-2)-c(j,nres-3)
102           c(j,nres)=c(j,nres-1)+dcj
103           c(j,2*nres)=c(j,nres)
104         enddo
105       endif
106       do i=2,nres-1
107         do j=1,3
108           c(j,i+nres)=dc(j,i)
109         enddo
110       enddo
111       do j=1,3
112         c(j,nres+1)=c(j,1)
113         c(j,2*nres)=c(j,nres)
114       enddo
115       if (itype(1).eq.21) then
116         nsup=nsup-1
117         nstart_sup=2
118         do j=1,3
119           dcj=c(j,4)-c(j,3)
120           c(j,1)=c(j,2)-dcj
121           c(j,nres+1)=c(j,1)
122         enddo
123       endif
124 C Calculate internal coordinates.
125       do ires=1,nres
126         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') 
127      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
128      &    (c(j,nres+ires),j=1,3)
129       enddo
130       call int_from_cart(.true.,.false.)
131       do i=1,nres-1
132         do j=1,3
133           dc(j,i)=c(j,i+1)-c(j,i)
134           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
135         enddo
136       enddo
137       do i=2,nres-1
138         do j=1,3
139           dc(j,i+nres)=c(j,i+nres)-c(j,i)
140           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
141         enddo
142 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
143 c     &   vbld_inv(i+nres)
144       enddo
145 c      call chainbuild
146 C Copy the coordinates to reference coordinates
147       do i=1,2*nres
148         do j=1,3
149           cref(j,i)=c(j,i)
150         enddo
151       enddo
152
153       ishift_pdb=ishift
154       return
155       end
156 c---------------------------------------------------------------------------
157       subroutine int_from_cart(lside,lprn)
158       implicit real*8 (a-h,o-z)
159       include 'DIMENSIONS'
160       include 'COMMON.LOCAL'
161       include 'COMMON.VAR'
162       include 'COMMON.CHAIN'
163       include 'COMMON.INTERACT'
164       include 'COMMON.IOUNITS'
165       include 'COMMON.GEO'
166       include 'COMMON.NAMES'
167       character*3 seq,atom,res
168       character*80 card
169       dimension sccor(3,20)
170       integer rescode
171       logical lside,lprn
172       if (lprn) then 
173         write (iout,'(/a)') 
174      &  'Internal coordinates calculated from crystal structure.'
175         if (lside) then 
176           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
177      & '       Phi','    Dsc_id','       Dsc','     Alpha',
178      & '     Omega'
179         else 
180           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
181      & '       Phi'
182         endif
183       endif
184       call flush(iout)
185       do i=nnt+1,nct
186         iti=itype(i)
187 c        write (iout,*) i,dist(i,i-1)
188         if (dist(i,i-1).lt.2.0D0 .or. dist(i,i-1).gt.5.0D0) then
189           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
190           stop
191         endif
192         theta(i+1)=alpha(i-1,i,i+1)
193         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
194       enddo
195       if (lside) then
196         do i=2,nres-1
197           do j=1,3
198             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
199           enddo
200           iti=itype(i)
201           di=dist(i,nres+i)
202           if (iti.ne.10) then
203             alph(i)=alpha(nres+i,i,maxres2)
204             omeg(i)=beta(nres+i,i,maxres2,i+1)
205           endif
206           if (lprn)
207      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
208      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,rad2deg*alph(i),
209      &    rad2deg*omeg(i)
210         enddo
211       else if (lprn) then
212         do i=2,nres
213           iti=itype(i)
214           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
215      &    rad2deg*theta(i),rad2deg*phi(i)
216         enddo
217       endif
218       return
219       end
220 c---------------------------------------------------------------------------
221       subroutine sccenter(ires,nscat,sccor)
222       implicit real*8 (a-h,o-z)
223       include 'DIMENSIONS'
224       include 'COMMON.CHAIN'
225       dimension sccor(3,20)
226       do j=1,3
227         sccmj=0.0D0
228         do i=1,nscat
229           sccmj=sccmj+sccor(j,i) 
230         enddo
231         dc(j,ires)=sccmj/nscat
232       enddo
233       return
234       end