1 subroutine readpdb(lprn,iprot,efree_temp,*)
2 C Read the PDB file and convert the peptide geometry into virtual-chain
6 include 'DIMENSIONS.ZSCOPT'
7 include 'COMMON.CONTROL'
10 include 'COMMON.CHAIN'
11 include 'COMMON.INTERACT'
12 include 'COMMON.IOUNITS'
14 include 'COMMON.NAMES'
16 integer i,j,iprot,ibeg,ishift1,ires,iires,iii,ires_old,ishift,ity,
19 double precision e1(3),e2(3),e3(3)
20 double precision dcj,efree_temp
21 character*3 seq,atom,res
22 character*80 card,prevcard
23 double precision sccor(3,20)
24 integer rescode,iterter(maxres)
29 c write (iout,*) "READPDB: UNRES_PDB",unres_pdb
34 read (ipdbin,'(a80)',end=10) card
35 do while (card(:3).eq.'END' .or. card(:3).eq.'TER')
36 read (ipdbin,'(a80)',end=10) card
39 write (iout,'(a)') card
41 c if (card(:3).eq.'END') goto 10
42 if (card(:3).eq.'END' .or.
43 & card(:3).eq.'TER'.and.prevcard(:5).eq.'REMARK') goto 10
44 if (card(:3).eq.'TER') then
46 write (iout,*) "Im here!"
47 write (iout,*) "ires_old",ires_old
49 itype(ires_old-1)=ntyp1
54 write (iout,*) "ires_old",ires_old
55 write (iout,*) "Chain ended",ires,ishift,ires_old
58 dc(j,ires)=sccor(j,iii)
61 call sccenter(ires,iii,sccor)
66 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
67 C Fish out the ATOM cards.
68 if (index(card(1:4),'ATOM').gt.0) then
69 read (card(14:16),'(a3)') atom
70 c if (atom.eq.'CA' .or. atom.eq.'CH3') then
71 read (card(23:26),*) ires
72 read (card(18:20),'(a3)') res
73 c write (iout,*) "ires",ires,ires-ishift+ishift1,
74 c & " ires_old",ires_old
75 c write (iout,*) "ishift",ishift," ishift1",ishift1
76 c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
78 iires=ires-ishift+ishift1
79 if (iires.ne.ires_old) then
80 C Calculate the CM of the preceding residue.
82 c write (iout,*) "Calculating sidechain center iii",iii
86 dc(j,ires)=sccor(j,iii)
89 call sccenter(iires-1,iii,sccor)
91 c write (iout,*) "sidechain center position calculated"
96 if (res.eq.'Cl-' .or. res.eq.'Na+') then
98 read (ipdbin,'(a80)',end=10) card
100 else if (ibeg.eq.1) then
101 c write (iout,*) "BEG ires",ires," iires",iires
104 if (res.ne.'GLY' .and. res.ne. 'ACE') then
108 ires=ires-ishift+ishift1
110 c write (iout,*) "ishift",ishift," ires",ires,
111 c & " iires",iires," ires_old",ires_old
115 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
116 ires=ires-ishift+ishift1
118 c write (iout,*) "RESSHIFT",ires,iires,ishift,ishift1
121 c write (iout,*) "ires",ires," iires",iires
123 if (res.eq.'ACE' .or. res.eq.'NHE') then
126 itype(ires)=rescode(ires,res,0)
129 ires=ires-ishift+ishift1
131 c write (iout,*) "ires_old",ires_old," ires",ires," iires",iires
133 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
136 c write (2,*) "iires",iires," res ",res," ity",ity
137 if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
138 & res.eq.'NHE'.and.atom(:2).eq.'HN') then
139 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
141 write (iout,'(2i3,2x,a,3f8.3)')
142 & iires,itype(iires),res,(c(j,ires),j=1,3)
146 sccor(j,iii)=c(j,ires)
148 c write (*,*) card(23:27),ires,itype(ires)
149 else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and.
150 & atom.ne.'N ' .and. atom.ne.'C ' .and.
151 & atom.ne.'OXT') then
153 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
156 read (ipdbin,'(a80)',end=10) card
158 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
160 if (ires.eq.0) return1
161 C Calculate the CM of the last side chain.
165 dc(j,ires)=sccor(j,iii)
168 call sccenter(iires,iii,sccor)
174 if (itype(nres).ne.10) then
178 dcj=c(j,nres-2)-c(j,nres-3)
179 c(j,nres)=c(j,nres-1)+dcj
180 c(j,2*nres)=c(j,nres)
190 c(j,2*nres)=c(j,nres)
192 if (itype(1).eq.ntyp1) then
193 nsup(iprot)=nsup(iprot)-1
196 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
197 call refsys(2,3,4,e1,e2,e3,fail)
204 c(j,1)=c(j,2)-3.8d0*e2(j)
214 C Copy the coordinates to reference coordinates
217 c cref(j,i,iprot)=c(j,i)
220 C Calculate internal coordinates.
223 & "Cartesian coordinates of the reference structure"
224 write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
225 & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
227 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
228 & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
229 & (c(j,ires+nres),j=1,3)
232 call int_from_cart(.true.,lprn)
233 call sc_loc_geom(lprn)
235 c phi_ref(i,iprot)=phi(i)
236 c theta_ref(i,iprot)=theta(i)
237 c alph_ref(i,iprot)=alph(i)
238 c omeg_ref(i,iprot)=omeg(i)
244 c---------------------------------------------------------------------------
245 subroutine int_from_cart(lside,lprn)
248 include 'DIMENSIONS.ZSCOPT'
249 include 'COMMON.LOCAL'
251 include 'COMMON.CHAIN'
252 include 'COMMON.INTERACT'
253 include 'COMMON.IOUNITS'
255 include 'COMMON.NAMES'
257 double precision dist,alpha,beta,di
258 character*3 seq,atom,res
260 double precision sccor(3,20)
265 & 'Internal coordinates of the reference structure.'
267 write (iout,'(8a)') 'Residue',' dvb',' Theta',
268 & ' Phi',' Dsc_id',' Dsc',' Alpha',
271 write (iout,'(4a)') 'Residue',' dvb',' Theta',
277 if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and.
278 & (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.6.5D0)) then
279 write (iout,'(a,i4)') 'Bad Cartesians for residue',i
282 vbld(i+1)=dist(i,i+1)
283 vbld_inv(i+1)=1.0d0/vbld(i+1)
284 if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
285 if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
287 dc(j,i)=c(j,i+1)-c(j,i)
288 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
294 c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
299 if (itype(i).ne.10) then
300 vbld_inv(i+nres)=1.0d0/di
302 dc(j,i+nres)=c(j,i+nres)-c(j,i)
303 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
306 vbld_inv(i+nres)=0.0d0
309 dc_norm(j,i+nres)=0.0d0
313 alph(i)=alpha(nres+i,i,maxres2)
314 omeg(i)=beta(nres+i,i,maxres2,i+1)
317 & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
318 & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
319 & rad2deg*alph(i),rad2deg*omeg(i)
326 write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
327 & rad2deg*theta(i),rad2deg*phi(i)
332 c-------------------------------------------------------------------------------
333 subroutine sc_loc_geom(lprn)
334 implicit real*8 (a-h,o-z)
336 include 'DIMENSIONS.ZSCOPT'
337 include 'COMMON.LOCAL'
339 include 'COMMON.CHAIN'
340 include 'COMMON.INTERACT'
341 include 'COMMON.IOUNITS'
343 include 'COMMON.NAMES'
344 include 'COMMON.CONTROL'
345 double precision x_prime(3),y_prime(3),z_prime(3)
349 dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
353 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
355 dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
359 dc_norm(j,i+nres)=0.0d0
364 costtab(i+1) =dcos(theta(i+1))
365 sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
366 cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
367 sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
368 cosfac2=0.5d0/(1.0d0+costtab(i+1))
369 cosfac=dsqrt(cosfac2)
370 sinfac2=0.5d0/(1.0d0-costtab(i+1))
371 sinfac=dsqrt(sinfac2)
373 if (it.ne.10 .and. itype(i).ne.ntyp1) then
375 C Compute the axes of tghe local cartesian coordinates system; store in
376 c x_prime, y_prime and z_prime
384 x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
385 y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
387 call vecpr(x_prime,y_prime,z_prime)
389 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
390 C to local coordinate system. Store in xx, yy, zz.
396 xx = xx + x_prime(j)*dc_norm(j,i+nres)
397 yy = yy + y_prime(j)*dc_norm(j,i+nres)
398 zz = zz + z_prime(j)*dc_norm(j,i+nres)
413 write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxtab(i),yytab(i),
419 c---------------------------------------------------------------------------
420 subroutine sccenter(ires,nscat,sccor)
423 include 'DIMENSIONS.ZSCOPT'
424 include 'COMMON.CHAIN'
425 integer ires,nscat,i,j
426 double precision sccmj
427 double precision sccor(3,20)
431 sccmj=sccmj+sccor(j,i)
433 dc(j,ires)=sccmj/nscat