source: pdbread.f

Last change on this file was 38d77eb, checked in by baerbaer <baerbaer@…>, 14 years ago

Redirected standard out to logString.

SMMP produced a lot of log messages. This became an issue when run in massively
parallel environments. I replaced all writes to standard out to a write to logString.
The next step is to pass this string to a function that writes the messages to a log
file according to their importance and the chosen log level.

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

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