source: pdbread.f@ e40e335

Last change on this file since e40e335 was e40e335, checked in by baerbaer <baerbaer@…>, 16 years ago

Initial import to BerliOS corresponding to 3.0.4

git-svn-id: svn+ssh://svn.berlios.de/svnroot/repos/smmp/trunk@1 26dc1dd8-5c4e-0410-9ffe-d298b4865968

  • Property mode set to 100644
File size: 13.3 KB
Line 
1c **************************************************************
2c
3c This file contains the subroutines: pdbread,pdbvars,atixpdb,getpar
4c
5c Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
6c Shura Hayryan, Chin-Ku
7c Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
8c Jan H. Meinke, Sandipan Mohanty
9c
10c **************************************************************
11
12 subroutine pdbread(pdbfil,ier)
13
14c ....................................................
15c PURPOSE: read protein atom coordinates from 'pdbfil'
16c (no Hydrogens, only ATOM records)
17c
18c RETURNS: 0 = no errors / 1 = error
19c
20c CALLS: iopfil,iendst
21c ......................................................
22
23 implicit real*8 (a-h,o-z)
24 implicit integer*4 (i-n)
25
26 include 'INCP.H'
27
28c -------------------------- input
29 character*(*) pdbfil
30c -------------------------- local
31 dimension cor(3)
32 character atm*4,rsn*3,rsno*3,chn,chno,
33 # rsid*5,rsido*5,line*132
34
35 natp=0
36 nchp=0
37 nrsp=0
38
39 ier=1
40
41 chno='&'
42 rsno='#&#'
43 rsido='#&#&#'
44
45 l=iendst(pdbfil)
46 if (l.gt.0) then
47 lunpdb = 99
48 else
49 write (*,'(a)')
50 # ' pdbread> empty file name to read pdb-structure'
51
52 return
53 endif
54
55 io=iopfil(lunpdb,pdbfil,'old','formatted')
56
57 if (io.le.0) then
58 write (*,'(a,/,a)')
59 # ' pdbread> ERROR opening file to read pdb-structure: ',
60 # pdbfil(1:iendst(pdbfil))
61
62 return
63 endif
64
65 1 read (lunpdb,'(a)',end=3,err=2) line
66 l=iendst(line)
67
68 if (l.lt.54.or.index(line(1:4),'ATOM').le.0) goto 1
69
70 if ( line(17:17).ne.' ' ) then
71 write (*,'(a,/,a,/,a,/,2a)')
72 # ' pdbread> found alternate atom location: ',
73 # ' !',
74 # line(:l),' in file: ',pdbfil(1:iendst(pdbfil))
75
76 close(lunpdb)
77 return
78 endif
79
80 atm=line(13:16)
81
82 if (index(atm(2:2),'H').gt.0) goto 1 ! no H
83
84 read(line,10,err=2) iat,rsn,chn,rsid,(cor(i),i=1,3)
85
86 if ((natp+1).gt.MXATP) then
87 write (*,'(a,i5,a,/,a)')
88 # ' pdbread> >MXATP (',MXATP,') ATOM lines in PDB file ',
89 # pdbfil(1:iendst(pdbfil))
90
91 close(lunpdb)
92 return
93 endif
94
95 if (chn.ne.chno) then ! new chain
96
97 if ((nchp+1).gt.MXCHP) then
98 write (*,'(a,i3,a,/,a)')
99 # ' pdbread> >MXCHP (',MXCHP,') chains in PDB file ',
100 # pdbfil(1:iendst(pdbfil))
101
102 close(lunpdb)
103 return
104 endif
105
106 if ((nrsp+1).gt.MXRSP) then
107 write (*,'(a,i3,a,/,a)')
108 # ' pdbread> >MXRSP (',MXRSP,') residues in PDB file ',
109 # pdbfil(1:iendst(pdbfil))
110
111 close(lunpdb)
112 return
113 endif
114
115 if (nchp.eq.1) then
116 nchrsp(nchp)=nrsp
117 elseif (nchp.gt.1) then
118 nchrsp(nchp)=nrsp-nchrsp(nchp)
119 endif
120
121 nchp=nchp+1
122 chno=chn
123 chnp(nchp)=chn
124
125 nchrsp(nchp)=nrsp ! -1 1st res.
126
127 if (nrsp.ge.1) then
128 nrsatp(nrsp)=natp-irsatp(nrsp)+1
129 endif
130
131 rsido=rsid
132 rsno=rsn
133
134 nrsp=nrsp+1
135 irsatp(nrsp)=natp+1
136 rsnmp(nrsp)=rsn
137 rsidp(nrsp)=rsid
138
139 elseif (rsid.ne.rsido.or.rsn.ne.rsno) then ! new residue
140
141 if ((nrsp+1).gt.MXRSP) then
142 write (*,'(a,i3,a,/,a)')
143 # ' pdbread> >MXRSP (',MXRSP,') residues in PDB file ',
144 # pdbfil(1:iendst(pdbfil))
145
146 close(lunpdb)
147 return
148 endif
149
150 nrsatp(nrsp)=natp-irsatp(nrsp)+1
151
152 rsido=rsid
153 rsno=rsn
154
155 nrsp=nrsp+1
156 irsatp(nrsp)=natp+1
157 rsnmp(nrsp)=rsn
158 rsidp(nrsp)=rsid
159
160 endif
161
162 natp=natp+1
163
164 noatp(natp)=iat
165 atnmp(natp)=atm
166
167 xatp(natp)=cor(1)
168 yatp(natp)=cor(2)
169 zatp(natp)=cor(3)
170
171 goto 1
172
173 2 write (*,'(a,/,a,/,2a)')
174 # ' pdbread> ERROR reading ATOM line ',
175 # line(:l),
176 # ' from file ',pdbfil(1:iendst(pdbfil))
177
178 close(lunpdb)
179 return
180
181 3 close(lunpdb)
182
183 if (natp.gt.0) then
184
185 if (nchp.eq.1) then
186 nchrsp(nchp)=nrsp
187 elseif (nchp.gt.1) then
188 nchrsp(nchp)=nrsp-nchrsp(nchp)
189 endif
190
191 nrsatp(nrsp)=natp-irsatp(nrsp)+1
192 ier=0
193
194 else
195
196 write (*,'(a,/,a)')
197 # ' pdbread> NO atom coordinates selected from file ',
198 # pdbfil(1:iendst(pdbfil))
199
200 endif
201
202 return
203
204 10 format(6x,i5,6x,a3,1x,a1,a5,3x,3d8.3)
205
206 end
207c **************************************************************
208
209 subroutine pdbvars()
210
211c --------------------------------------------------------------------
212c PURPOSE: sequence,indices for selected atoms (data in INCP.H)
213c & torsions from PDB to be used to build SMMP structure
214c
215c ixatp(i,)
216c = indices for SMMP atoms pointing to PDB atoms
217c (=0, if atom not selected)
218c
219c --------------------------------- ref. point & axes
220c ixrfpt(3,),rfpt(3,),xrfax(3,),yrfax(3,),zrfax(3,)
221c
222c CALLS: tolost,getmol,bldmol,addend,atixpdb,setmvs,mklist,
223c dihedr,fnd3ba,setsys,getpar,setvar,rmsdopt
224c --------------------------------------------------------------------
225
226 include 'INCL.H'
227 include 'INCP.H'
228
229 character res*3
230 dimension rm(3,3),av1(3),av2(3),h(3)
231
232
233 nml=0
234 nrs=0
235
236 do nc=1,nchp ! PDB chains
237
238c =============================== SMMP molecule
239 nml=nml+1
240 if (nml.gt.mxml) then
241 write(*,'(a,i4,2a)')' pdbvars> NUMBER of chains > '
242 # ,mxml,' in ',' ?'
243 stop
244 endif
245 ntlml=nml
246c ----------------------------- 'nmml' = ChainID
247 nmml(nml)=chnp(nc)
248
249c ======================================== get sequence
250
251 irb=nrs+1
252 ire=nrs+nchrsp(nc)
253c ----------------------------- # of 1st & last residue
254 irsml1(nml)=irb
255 irsml2(nml)=ire
256
257 do irs=irb,ire ! residues of chain 'nc'
258
259 nrs=nrs+1
260
261 if (nrs.gt.mxrs) then
262 write(*,'(a,i4,2a)') ' pdbvars> NUMBER of residues > '
263 # ,mxrs,' in ',' ?'
264 stop
265 endif
266
267 res=rsnmp(irs)
268 call tolost(res)
269
270 seq(nrs)=res
271
272 if (.not.flex.and.irs.eq.irb.and.seq(nrs)(1:3).eq.'pro')
273 # seq(nrs)='pron' ! only ECEPP/3
274
275 enddo ! residues
276
277c ======================== get initial coords. for molecule 'nml'
278c with library values for deg. of freedom
279
280 call getmol(nml) ! assemble res. data from libraries
281
282 do i=1,6 ! initialize global parameters
283 gbpr(i,nml)=zero
284 enddo
285
286 call bldmol(nml) ! co-ordinates
287 call addend(nml,'nh2 ','cooh') ! modify ends
288
289 call atixpdb(nml) ! get 'ixatp'
290
291c -------------------------- 'load' SMMP variable information
292 call setmvs(nml) ! moving sets
293 call mklist(nml) ! interaction lists
294
295c ================================= get variables for 'nml'
296
297 ii=ivrml1(nml)
298 do i=ii,ii+nvrml(nml)-1 ! SMMP torsions
299
300 isrfvr(i) = .false.
301 fxvr(i) = .false.
302 idvr(i) = i
303
304 it = ityvr(i)
305 i1 = iatvr(i)
306
307 if (it.eq.3) then ! torsion
308
309 i2=iowat(i1) ! indices of SMMP atoms
310 i3=iowat(i2)
311 i4=iowat(i3)
312
313 j1=ixatp(i1) ! inds. for corresp. PDB atoms
314 j2=ixatp(i2)
315 j3=ixatp(i3)
316 j4=ixatp(i4)
317
318 if (j1.le.0.or.j2.le.0.or.j3.le.0.or.j4.le.0) then
319
320 vlvr(i) = toat(i1) ! default value from library
321
322 else ! get value from PDB
323
324 xat(i1)=xatp(j1)
325 yat(i1)=yatp(j1)
326 zat(i1)=zatp(j1)
327
328 xat(i2)=xatp(j2)
329 yat(i2)=yatp(j2)
330 zat(i2)=zatp(j2)
331
332 xat(i3)=xatp(j3)
333 yat(i3)=yatp(j3)
334 zat(i3)=zatp(j3)
335
336 xat(i4)=xatp(j4)
337 yat(i4)=yatp(j4)
338 zat(i4)=zatp(j4)
339
340 vlvr(i) = dihedr(i1,i2,i3,i4)
341
342 isrfvr(i) = .true.
343
344 endif
345
346 elseif (it.eq.2) then ! b.angle
347 vlvr(i)=baat(i1)
348 elseif (it.eq.1) then ! b.length
349 vlvr(i)=blat(i1)
350 endif
351
352 olvlvr(i) = vlvr(i)
353
354 enddo ! SMMP vars.
355
356 nvr = ivrml1(ntlml)+nvrml(ntlml)-1
357
358c ================================= global parameters for 'nml'
359
360c +++++++++++
361 inew=0
362
363 if (inew.eq.1) then
364c ++++++++++++++++++++++++
365
366 call setvar(nml,vlvr)
367
368 nrs = irsml2(nml)-irsml1(nml)+1
369 call rmsdopt(nml,1,nrs,ixatp,xatp,yatp,zatp,0,rm,av1,av2,rmsd)
370
371c ---------------------------- retrieve ref. coords.
372c & transform acc. to opt. rmsd
373 do i=1,3
374 ii=ixrfpt(i,nml)
375
376 h(1)=xat(ii)-av1(1)
377 h(2)=yat(ii)-av1(2)
378 h(3)=zat(ii)-av1(3)
379
380 x=av2(1)
381 y=av2(2)
382 z=av2(3)
383
384 do j=1,3
385 x = x + rm(j,1) * h(j)
386 y = y + rm(j,2) * h(j)
387 z = z + rm(j,3) * h(j)
388 enddo
389
390 xat(ii)=x
391 yat(ii)=y
392 zat(ii)=z
393 enddo
394
395 call getpar(nml)
396 call bldmol(nml) ! finally build SMMP molecule
397
398c ++++++++++++++++
399 else ! old
400c ++++++++++++++++
401
402 call fnd3ba(nml,i1,i2,i3) ! three 1st bb atoms in SMMP (e.g. n,ca,c')
403
404 ixrfpt(1,nml)=i1
405 ixrfpt(2,nml)=i2
406 ixrfpt(3,nml)=i3
407
408c -------------------------------- retrieve ref. coords.
409 do i=1,3
410 ii=ixrfpt(i,nml)
411 ix=ixatp(ii)
412 if (ix.gt.0) then
413 xat(ii)=xatp(ix)
414 yat(ii)=yatp(ix)
415 zat(ii)=zatp(ix)
416 else
417 write(*,'(3a)') ' pdbvars> missing PDB atom ',nmat(ii),
418 # ' is ref. point for SMMP - cannot proceed !'
419 endif
420 enddo
421
422 call getpar(nml)
423 call setvar(nml,vlvr) ! finally build SMMP molecule
424
425 nrs = irsml2(nml)-irsml1(nml)+1
426 call rmsdopt(nml,1,nrs,ixatp,xatp,yatp,zatp,0,rm,av1,av2,rmsd)
427
428c ++++++++++
429 endif
430c ++++++++++
431
432 write(*,*) ' '
433 write(*,*) ' Initial RMSD ',rmsd
434
435 enddo ! chains(molecules)
436
437 return
438 end
439c ***************************
440 subroutine atixpdb(nml)
441
442c --------------------------------------------------------------------
443c PURPOSE: get ixatp - pointer of each SMMP atom to corresponding atom
444c of reference structure loaded in 'INCP.H'
445c (=0 if no corr. atom in ref. str.)
446c
447c CALLS: toupst
448c --------------------------------------------------------------------
449
450 include 'INCL.H'
451 include 'INCP.H'
452
453 character*4 atm
454
455
456 atm=' '
457
458 do irs=irsml1(nml),irsml2(nml) ! SMMP residues
459
460 i1=irsatp(irs)
461 i2=i1+nrsatp(irs)-1
462
463 do iat=iatrs1(irs),iatrs2(irs) ! SMMP atoms
464
465 ix=0
466
467 if (nmat(iat)(1:1).ne.'h') then ! ignore hydrogens
468
469 atm(2:4)=nmat(iat)(1:3)
470 call toupst(atm)
471
472 do i=i1,i2 ! atoms of PDB residue
473 if (index(atnmp(i),atm).gt.0) then
474 ix=i
475 goto 1
476 endif
477 enddo
478
479c write(*,'(8a)') ' pdbvars> ',atm,' not found in '
480c # ,chnp(nc),' ',rsidp(irs),' ',rsnmp(irs)
481
482 endif
483
484 1 ixatp(iat)=ix
485
486 enddo ! SMMP atoms of 'irs'
487 enddo ! residues
488
489 return
490 end
491c **************************
492 subroutine getpar(nml)
493
494 include 'INCL.H'
495
496 parameter (TOL = 1.d-12)
497
498c Obtain molecule-fixed system (J,K,L) for 1st 3 bb-atoms,
499c -> determine global parameters: shifts dX,dY,dZ
500c & angles alpha,beta,gamma [rad], put into 'gbpr'
501c
502c CALLS: none
503c
504
505 i1=ixrfpt(1,nml) ! from 'INCL.H'
506 i2=ixrfpt(2,nml)
507 i3=ixrfpt(3,nml)
508c -------------------------------------- Shifts
509 gbpr(1,nml) = xat(i1)
510 gbpr(2,nml) = yat(i1)
511 gbpr(3,nml) = zat(i1)
512
513 do i = 4,6
514 gbpr(i,nml) = 0.d0
515 enddo
516c --------------------------------- J
517 h1=xat(i2)
518 h2=yat(i2)
519 h3=zat(i2)
520
521 x1=h1-xat(i1)
522 x2=h2-yat(i1)
523 x3=h3-zat(i1)
524
525 d=sqrt(x1*x1+x2*x2+x3*x3)
526
527 x1=x1/d
528 x2=x2/d
529 x3=x3/d
530c --------------------------------- L
531 h1=xat(i3)-h1
532 h2=yat(i3)-h2
533 h3=zat(i3)-h3
534
535 z1=x2*h3-x3*h2
536 z2=x3*h1-x1*h3
537 z3=x1*h2-x2*h1
538
539 d=sqrt(z1*z1+z2*z2+z3*z3)
540
541 z1=z1/d
542 z2=z2/d
543 z3=z3/d
544
545c ---------------------------------- K
546 y1=z2*x3-z3*x2
547 y2=z3*x1-z1*x3
548 y3=z1*x2-z2*x1
549
550 if ( ( 1.d0 - abs(y3) ) .gt. TOL ) then ! ============ |beta| < PI/2
551
552c ----------------------------------------------- Y'
553 d = sqrt( y1 * y1 + y2 * y2 )
554 yp1= y1 / d
555 yp2= y2 / d
556
557 gbpr(4,nml)= atan2( -yp1, yp2 ) ! alpha
558 gbpr(5,nml)= atan2( y3, ( y1*yp1+y2*yp2 ) ) ! beta
559 gbpr(6,nml)= atan2( ( z1*yp2-z2*yp1 ),( x1*yp2-x2*yp1 ) ) ! gamma
560
561 else ! ======================= |beta| = PI/2
562
563 gbpr(4,nml) = atan2( x2, x1 ) ! alpha+gamma
564
565 if ( abs(y3) .lt. 1.d0 ) then ! beta
566 gbpr(5,nml) = asin(y3)
567 else if ( y3 .gt. 0.d0 ) then
568 gbpr(5,nml) = pi*.5d0
569 else
570 gbpr(5,nml) = -pi*.5d0
571 endif
572
573 gbpr(6,nml) = 0.0
574
575 endif
576
577 return
578 end
579
Note: See TracBrowser for help on using the repository browser.