1 subroutine smd_pdb_natoms(filename,nt) 2 implicit none 3 character*(*) filename 4 integer nt 5c 6 character*(4) buffer 7 integer un 8c 9 un = 70 10 open(unit=un,status="old",form="formatted",file=filename) 11 12 nt = 0 13100 continue 14 read(un,'(A4)',end=200) buffer 15 if(buffer(1:4).eq."ATOM") then 16 nt = nt + 1 17 end if 18 goto 100 19200 continue 20 close(un) 21 22 end 23 24 subroutine smd_pdb_read_coords(filename,nt,c) 25 implicit none 26 character*(*) filename 27 integer nt 28 double precision c(3,nt) 29c 30 character*(180) buffer 31 character*(4) tag 32 integer i 33 integer un 34 character*(30) pname 35 36 pname = "smd_pdb_read_coords" 37c 38 un = 70 39c 40 open(unit=un,status="old",form="formatted",file=filename) 41 i = 0 42100 continue 43 read(un,'(A180)',end=200) buffer 44 if(buffer(1:4).eq."ATOM") then 45 i = i +1 46 read(buffer,*) tag,tag,tag,tag,tag, 47 > c(1,i),c(2,i),c(3,i) 48 end if 49 goto 100 50200 continue 51 close(un) 52 53 return 54 55 end 56 57 subroutine smd_pdb_read_atomres(filename,nt,ta,tr,ir) 58 implicit none 59 character*(*) filename 60 integer nt 61 character*(*) ta(nt) 62 character*(*) tr(nt) 63 integer ir(nt) 64c 65 character*(180) buffer 66 integer i 67 integer un 68 character*(30) pname 69 70 pname = "smd_pdb_read_atomres" 71c 72 un = 70 73c 74 open(unit=un,status="old",form="formatted",file=filename) 75 i = 0 76100 continue 77 read(un,'(A180)',end=200) buffer 78 if(buffer(1:4).eq."ATOM") then 79 i = i +1 80 ta(i) = buffer(13:16) 81 tr(i) = buffer(18:20) 82 read(buffer(23:26),*) ir(i) 83 end if 84 goto 100 85200 continue 86 close(un) 87 88 return 89 90 end 91 92 subroutine smd_pdb_sort_byres(nt,ta,tr,ir,c) 93 implicit none 94 integer nt 95 character*(*) ta(nt) 96 character*(*) tr(nt) 97 integer ir(nt) 98 double precision c(3,nt) 99c 100 character*(180) buffer 101 integer i 102 integer un 103 character*(30) pname 104 integer pass 105 integer sorted 106 integer itemp 107 double precision ftemp 108 character*16 stemp 109 110 pass = 1 111 sorted = 0 112 do while(sorted .eq. 0) 113 sorted = 1 114 do i = 1,nt-pass 115 if(ir(i) .gt. ir(i+1)) then 116 itemp = ir(i) 117 ir(i) = ir(i+1) 118 ir(i+1) = itemp 119c 120 stemp = ta(i) 121 ta(i) = ta(i+1) 122 ta(i+1) = stemp 123c 124 stemp = tr(i) 125 tr(i) = tr(i+1) 126 tr(i+1) = stemp 127c 128 ftemp = c(1,i) 129 c(1,i) = c(1,i+1) 130 c(1,i+1) = ftemp 131c 132 ftemp = c(2,i) 133 c(2,i) = c(2,i+1) 134 c(2,i+1) = ftemp 135c 136 ftemp = c(3,i) 137 c(3,i) = c(3,i+1) 138 c(3,i+1) = ftemp 139c 140 sorted = 0 141 endif 142 enddo 143 pass = pass +1 144 end do 145 146 return 147 148 end 149 150 subroutine smd_pdb_sort_seq_distance(nr,is,im,d) 151 implicit none 152 integer nr 153 integer is(nr) 154 integer im(nr) 155 double precision d(nr) 156c 157 integer i 158 integer pass 159 integer sorted 160 integer itemp 161 double precision ftemp 162 163 pass = 1 164 sorted = 0 165 do while(sorted .eq. 0) 166 sorted = 1 167 do i = 1,nr-pass 168 if(d(i) .gt. d(i+1)) then 169 ftemp = d(i) 170 d(i) = d(i+1) 171 d(i+1) = ftemp 172c 173 itemp = is(i) 174 is(i) = is(i+1) 175 is(i+1) = itemp 176c 177c itemp = im(i) 178c im(i) = im(i+1) 179c im(i+1) = itemp 180 181 sorted = 0 182 endif 183 enddo 184 pass = pass +1 185 end do 186 187 return 188 189 end 190 191 subroutine smd_pdb_read(filename,nt,ta,tr,ir,c) 192 implicit none 193 character*(*) filename 194 integer nt 195 character*(16) ta(nt) 196 character*(16) tr(nt) 197 integer ir(nt) 198 double precision c(3,nt) 199c 200 character*(180) buffer 201 character*(4) tag 202 integer i 203 integer un 204 character*(30) pname 205 206 pname = "smd_pdb_read" 207c 208 un = 70 209c 210 open(unit=un,status="old",form="formatted",file=filename) 211 i = 0 212100 continue 213 read(un,'(A180)',end=200) buffer 214 if(buffer(1:4).eq."ATOM") then 215 i = i +1 216 read(buffer,*) tag,tag,ta(i),tr(i),ir(i), 217 > c(1,i),c(2,i),c(3,i) 218 end if 219 goto 100 220200 continue 221 close(un) 222 223 return 224 225 end 226 227 subroutine smd_pdb_write_byseq(filename,nr,nt,im,is,ta,tr,c) 228 implicit none 229 character*(*) filename 230 integer nr 231 integer im(nr+1) 232 integer is(nr) 233 integer nt 234 character*(16) ta(nt) 235 character*(16) tr(nt) 236 double precision c(3,nt) 237c 238 character*(180) buffer 239 character*(4) tag 240 integer i 241 integer un 242 integer i0,j0,j,jb,je 243c 244 un = 70 245c 246 open(unit=un,status="unknown",form="formatted",file=filename) 247 j0=0 248 do i=1,nr 249 i0=is(i) 250 jb = im(i0) 251 je = im(i0+1)-1 252 do j=jb,je 253 j0=j0+1 254 write(un,FMT=9000)j0,ta(j),tr(j),i,c(1,j),c(2,j),c(3,j) 255 end do 256 end do 257 close(un) 258 259 return 2609000 FORMAT("ATOM",T7,I5,T13,A4,T18,A3,T23,I4,T31, 261 > F8.3,T39,F8.3,T47,F8.3,T55,F6.2) 262 263 end 264 265 subroutine smd_pdb_nres(filename,nr) 266 implicit none 267 character*(*) filename 268 integer nr 269c 270 character*(180) buffer 271 character*(4) tag 272 integer ir0,ir 273 integer un 274c 275 un = 70 276c 277 open(unit=un,status="old",form="formatted",file=filename) 278 279c reset residue arrays to be the size of number of residues only 280 nr = 0 281 ir0 = 0 282100 continue 283 read(un,'(A180)',end=200) buffer 284 if(buffer(1:4).eq."ATOM") then 285 read(buffer,*) tag,tag,tag,tag,ir 286 if(ir0.ne.ir) then 287 nr = nr + 1 288 ir0=ir 289 end if 290 end if 291 goto 100 292200 continue 293 close(un) 294 295 end 296 297 subroutine smd_pdb_nres0(nr,nt,ir) 298 implicit none 299 integer nr 300 integer nt 301 integer ir(nt) 302c 303 integer i,ir0 304 nr = 0 305 ir0 = 0 306 do i=1,nt 307 if(ir0.ne.ir(i)) then 308 nr = nr + 1 309 ir0=ir(i) 310 end if 311 end do 312 313 end 314 315 subroutine smd_pdb_cog(nt,nr,ir,c,cg) 316 implicit none 317 integer nt 318 integer nr 319 integer ir(nt) 320 double precision c(3,nt) 321 double precision cg(3,nr) 322c 323 integer i,j,ir0,nm 324 j=1 325 ir0=ir(1) 326 nm=0 327 cg=0.0d0 328 do i=1,nt 329 if(ir(i).ne.ir0) then 330 cg(:,j)=cg(:,j)/nm 331 j=j+1 332 ir0=ir(i) 333 nm=0 334 end if 335 cg(:,j)=cg(:,j)+c(:,i) 336 nm=nm+1 337 end do 338 cg(:,j)=cg(:,j)/nm 339 return 340 end 341 342 subroutine smd_pdb_cog_byname(nt,nr,tar,ta,ir,c,cg) 343C nt total number of atoms 344C nr total number of residues 345C tar atom name 346C ta atom name array 347C ir residue index array 348C c coordinate array 349C cg center of mass array 350 implicit none 351 integer nt 352 integer nr 353 character*(*) tar 354 character*(*) ta(nt) 355 integer ir(nt) 356 double precision c(3,nt) 357 double precision cg(3,nr) 358 359c 360 integer i,j,ir0,nm 361 integer s0 362 integer length 363 external length 364 365 s0 = length(tar) 366 j=1 367C first residue 368 ir0=ir(1) 369 nm=0 370 cg=0.0d0 371 do i=1,nt 372 if(ir(i).ne.ir0) then 373 if(nm.ne.0) then 374 cg(:,j)=cg(:,j)/nm 375 else 376 write(*,*) "no atoms maching pattern ",tar(1:s0) 377 stop 378 end if 379 j=j+1 380 ir0=ir(i) 381 nm=0 382 end if 383 write(*,*) "atom name",tar 384C check for name match 385 if(INDEX(ta(i),tar(1:s0)).gt.0) then 386 cg(:,j)=cg(:,j)+c(:,i) 387 nm=nm+1 388 end if 389 end do 390 cg(:,j)=cg(:,j)/nm 391 return 392 end 393 394 subroutine smd_pdb_sequence_bounds(nt,nr,ir,is,im) 395 implicit none 396 integer nt 397 integer nr 398 integer ir(nt) 399 integer is(nr) 400 integer im(nr+1) 401c 402 integer i,j,ir0,nm 403 j = 1 404 is(1) = ir(1) 405 im(1) = 1 406 do i=1,nt 407 if(ir(i).ne.is(j)) then 408 j=j+1 409 is(j) = ir(i) 410 im(j) = i 411 end if 412 end do 413 im(j+1)=i 414 415 416 return 417 end 418 419 subroutine smd_pdb_distance(nt,c,cor,cd) 420 implicit none 421 integer nt 422 double precision c(3,nt) 423 double precision cor(3) 424 double precision cd(nt) 425c 426 double precision c1(3) 427 integer i 428 429 do i=1,nt 430 c1(:)=c(:,i)-cor(:) 431 c1(:)=c1(:)*c1(:) 432 cd(i)=SUM(c1) 433 end do 434 435 cd = sqrt(cd) 436 return 437 end 438 439 subroutine smd_pdb_read_res(filename,nt,nr,tr,ir,nm) 440 implicit none 441 character*(*) filename 442 integer nt,nr,nc 443 character*(16) tr(nr) 444 integer ir(nt) 445 integer nm(nr) 446c 447 character*(30) pname 448 character*(180) buffer 449 character*(4) tag 450 character*(16) rtag,rtag0 451 integer ir0,nr0 452 integer ncenter 453 integer un 454c 455 pname = "sg_read_res" 456c 457 un = 70 458c 459 open(unit=un,status="old",form="formatted",file=filename) 460 461 ncenter = 0 462 nr0 = 0 463 rtag0 = " " 464 ir0 = 0 465100 continue 466 read(un,'(A180)',end=200) buffer 467 if(buffer(1:4).eq."ATOM") then 468 ncenter = ncenter + 1 469 read(buffer,*) tag,tag,tag,rtag,ir(ncenter) 470 if(ir0.ne.ir(ncenter)) then 471 ir0=ir(ncenter) 472 nr0 = nr0 + 1 473 tr(nr0) = rtag 474 rtag0=rtag 475 end if 476 ir(ncenter) = nr0 477 nm(nr0) = nm(nr0) + 1 478 end if 479 goto 100 480200 continue 481 482 close(un) 483 484 return 485 end 486 487 function length(string) 488 implicit none 489 integer length 490*returns length of string ignoring trailing blanks 491 character*(*) string 492 integer i 493 do i = len(string), 1, -1 494 if(string(i:i) .ne. ' ') go to 20 495 end do 49620 length = i 497 end 498 499c $Id$ 500