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
Line 
1! **************************************************************
2!
3! This file contains the subroutines: pdbread,pdbvars,atixpdb,getpar
4!
5! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
6! Shura Hayryan, Chin-Ku
7! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
8! Jan H. Meinke, Sandipan Mohanty
9!
10! **************************************************************
11
12 subroutine pdbread(pdbfil,ier)
13
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! ......................................................
22
23 include 'INCP.H'
24
25 double precision cor
26
27 integer ier, l, iendst, lunpdb, io, iopfil, iat, i
28
29! -------------------------- input
30 character*(*) pdbfil
31! -------------------------- local
32 dimension cor(3)
33 character atm*4,rsn*3,rsno*3,chn,chno,
34 & rsid*5,rsido*5,line*132, logString*255
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
50 write (logString, '(a)')
51 & ' pdbread> empty file name to read pdb-structure'
52
53 return
54 endif
55
56 io=iopfil(lunpdb,pdbfil,'old','formatted')
57
58 if (io.le.0) then
59 write (logString, '(a,/,a)')
60 & ' pdbread> ERROR opening file to read pdb-structure: ',
61 & pdbfil(1:iendst(pdbfil))
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
72 write (logString, '(a,/,a,/,a,/,2a)')
73 & ' pdbread> found alternate atom location: ',
74 & ' !',
75 & line(:l),' in file: ',pdbfil(1:iendst(pdbfil))
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
88 write (logString, '(a,i5,a,/,a)')
89 & ' pdbread> >MXATP (',MXATP,') ATOM lines in PDB file ',
90 & pdbfil(1:iendst(pdbfil))
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
99 write (logString, '(a,i3,a,/,a)')
100 & ' pdbread> >MXCHP (',MXCHP,') chains in PDB file ',
101 & pdbfil(1:iendst(pdbfil))
102
103 close(lunpdb)
104 return
105 endif
106
107 if ((nrsp+1).gt.MXRSP) then
108 write (logString, '(a,i3,a,/,a)')
109 & ' pdbread> >MXRSP (',MXRSP,') residues in PDB file ',
110 & pdbfil(1:iendst(pdbfil))
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
143 write (logString, '(a,i3,a,/,a)')
144 & ' pdbread> >MXRSP (',MXRSP,') residues in PDB file ',
145 & pdbfil(1:iendst(pdbfil))
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
174 2 write (logString, '(a,/,a,/,2a)')
175 & ' pdbread> ERROR reading ATOM line ',
176 & line(:l),
177 & ' from file ',pdbfil(1:iendst(pdbfil))
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
197 write (logString, '(a,/,a)')
198 & ' pdbread> NO atom coordinates selected from file ',
199 & pdbfil(1:iendst(pdbfil))
200
201 endif
202
203 return
204
205 10 format(6x,i5,6x,a3,1x,a1,a5,3x,3d8.3)
206
207 end
208! **************************************************************
209
210 subroutine pdbvars()
211
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! --------------------------------------------------------------------
226
227 include 'INCL.H'
228 include 'INCP.H'
229
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
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
244! =============================== SMMP molecule
245 nml=nml+1
246 if (nml.gt.mxml) then
247 write (logString, '(a,i4,2a)')' pdbvars> NUMBER of chains > '
248 & ,mxml,' in ',' ?'
249 stop
250 endif
251 ntlml=nml
252! ----------------------------- 'nmml' = ChainID
253 nmml(nml)=chnp(nc)
254
255! ======================================== get sequence
256
257 irb=nrs+1
258 ire=nrs+nchrsp(nc)
259! ----------------------------- # of 1st & last residue
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
268 write (logString, '(a,i4,2a)')
269 & ' pdbvars> NUMBER of residues > ', mxrs, ' in ',' ?'
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')
279 & seq(nrs)='pron' ! only ECEPP/3
280
281 enddo ! residues
282
283! ======================== get initial coords. for molecule 'nml'
284! with library values for deg. of freedom
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
297! -------------------------- 'load' SMMP variable information
298 call setmvs(nml) ! moving sets
299 call mklist(nml) ! interaction lists
300
301! ================================= get variables for 'nml'
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
364! ================================= global parameters for 'nml'
365
366! +++++++++++
367 inew=0
368
369 if (inew.eq.1) then
370! ++++++++++++++++++++++++
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
377! ---------------------------- retrieve ref. coords.
378! & transform acc. to opt. rmsd
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
404! ++++++++++++++++
405 else ! old
406! ++++++++++++++++
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
414! -------------------------------- retrieve ref. coords.
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
423 write (logString, '(3a)') ' pdbvars> missing PDB atom ',
424 & nmat(ii), ' is ref. point for SMMP - cannot proceed !'
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
434! ++++++++++
435 endif
436! ++++++++++
437
438 write (logString, *) ' '
439 write (logString, *) ' Initial RMSD ',rmsd
440
441 enddo ! chains(molecules)
442
443 return
444 end
445! ***************************
446 subroutine atixpdb(nml)
447
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! --------------------------------------------------------------------
455
456 include 'INCL.H'
457 include 'INCP.H'
458
459 integer irs, nml, i1, i2, iat, ix, i
460
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
486
487! write (logString, '(8a)') ' pdbvars> ',atm,' not found in '
488! # ,chnp(nc),' ',rsidp(irs),' ',rsnmp(irs)
489
490 endif
491
492 1 ixatp(iat)=ix
493
494 enddo ! SMMP atoms of 'irs'
495 enddo ! residues
496
497 return
498 end
499! **************************
500 subroutine getpar(nml)
501
502 include 'INCL.H'
503
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
509 parameter (TOL = 1.d-12)
510
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'
514!
515! CALLS: none
516!
517
518 i1=ixrfpt(1,nml) ! from 'INCL.H'
519 i2=ixrfpt(2,nml)
520 i3=ixrfpt(3,nml)
521! -------------------------------------- Shifts
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
529! --------------------------------- J
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
543! --------------------------------- L
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
558! ---------------------------------- K
559 y1=z2*x3-z3*x2
560 y2=z3*x1-z1*x3
561 y3=z1*x2-z2*x1
562
563 if ( ( 1.d0 - abs(y3) ) .gt. TOL ) then ! ============ |beta| < PI/2
564
565! ----------------------------------------------- Y'
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.