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