1 PROGRAM RASTEP 2******************************************************************************** 3* 4* Usage: 5* rastep [-h] [-iso] [-Bcolor Bmin Bmax] [-prob xx] [-radius r] [-fancy[0-9]] 6* [-tabulate histogram.file] [-by_atomtype] [-suv_check] [-cn_check] 7* 8* 9* -auto auto-orientation of viewpoint 10* -h suppresses header records in output 11* -iso forces isotropic B values (spheres rather than 12* ellipsoids) even if ANISOU cards present 13* -aniso converts isotropic B values into Uij terms so that they 14* can be included in the analysis 15* -Bcolor Bmin Bmax 16* color by Biso values; Bmin = dark blue, Bmax = light red 17* -Acolor color by anisotropy; red < white (A=0.5) < green 18* -prob xx draws ellipsoids to enclose this 19* probability level (default = 0.50) 20* -radius draws bonds with this radius in Angstroms 21* (default = 0.10) 22* -fancy[0-9] increasingly complex rendition of ellipsoids 23* fancy0 [default] solid surface only 24* fancy1 principle axes + transparent bounding ellipsoid 25* fancy2 equatorial planes only 26* fancy3 equatorial planes + transparent ellipsoid 27* fancy4 longest principle axis only 28* fancy5 for ORTEP lovers - one octant missing 29* fancy6 same as fancy5, with missing octant colored grey 30* 31*=============================================================================== 32* The following options are used by parvati scripts 33* 34* -tabulate [file]instead of creating a Raster3D input file, 35* list all atoms with principle axes and anisotropy. 36* Optionally write a histogram of anisotropy to speficied 37* output file; otherwise output is to stderr 38* 39* output from -tabulate for this version of rastep 40* ATOM RESNAME RESNUM EIGEN1 EIGEN2 EIGEN3 ANISOTROPY Uiso 41* 42* -com [file] find <anisotropy> in shells from center of mass 43* -nohydrogens don't plot hydrogens even if present 44* -mini small size plot (176x208) with auto-orientation 45* -suv_check use Suv to validate similarity of bonded ellipsoids 46* -cn_check use CCuij to check similarity across C-N bonds 47* 48******************************************************************************** 49* 50* EAM Jul 97 - initial version 51* EAM Dec 97 - version 2.4b release 52* EAM Jan 98 - add tabulation option 53* EAM May 98 - integrate with PARVATI script 54* EAM Jun 99 - fix bug (lack of sqrt) in -iso processing, trap read errors 55* -Acolor flag to color by anisotropy (not yet entirely satisfactory) 56* EAM Jul 99 - V2.4l -tab output revised slightly, 57* NPD ellipsoids colored magenta 58* -nohydro flag to suppress drawing hydrogens 59* -mini flag to generate smaller pictures 60* don't draw bonds between atoms with different ALT flags 61* EAM Aug 99 - V2.4m 62* work harder at suppressing hydrogens 63* add auto-orientation (NB: scaling is wrong in this case!) 64* EAM Dec 99 - V2.5 65* clean up output formats a little 66* EAM Jun 2000 - additional error reporting 67* EAM Jul 2000 - apply Bcolor to bonds as well as atoms 68* EAM Sep 2000 - Suv similarity test 69* EAM Apr 2001 - V2.6 70* ORTEP_LIKE ellipsoids (one octant missing) 71* error count 72* EAM Feb 2002 - rework ANISOU and Suv processing to gain back some speed 73* maybe I should have an array of iso/aniso flags to save 74* time during testing? 75* the rest of CARD() array can go too? 76* all the tests on IF ATOM(I)(1:).eq.'ATOM' are now unneeded 77* EAM May 2003 - expand default color table to allow for off-by-one atom names 78* EAM Oct 2003 - trap and report 0 values for axial Uij components in input 79* EAM May 2006 - initial arrays to 0 for gfortran 80* EAM Feb 2008 - 2.7s initialize even more arrays for gfortran 81* EAM Aug 2009 - slightly revised output format for -tabulate 82* EAM Sep 2009 - Explicitly test CCuij for C-N bonds 83* EAM Oct 2009 - report CCuij = if C or N is NPD 84* -aniso flag 85* EAM Mar 2010 - CCuij for phosphate backbone also 86* 87* I/O units for colour/co-ordinate input, specs output, user output 88* 89 INCLUDE 'VERSION.incl' 90* 91 INTEGER INPUT, OUTPUT, NOISE 92 PARAMETER (INPUT=5, OUTPUT=6, NOISE=0) 93 PARAMETER (MAXCOL=5000, MAXATM=300000) 94 REAL RGB(3,MAXCOL), VDW(MAXCOL) 95 REAL SPAM(5,MAXATM) 96 REAL UIJ(6,MAXATM) 97 real center(3) 98 CHARACTER*24 MASK(MAXCOL),TEST 99 CHARACTER*80 ATOM(MAXATM),CARD 100 character*3 resname 101 character*1 spacer 102 LOGICAL MATCH 103 logical hflag, ellipses, bcflag, tflag, atflag, comflag 104 logical acflag, nohydro, mini, auto, aniflag 105 integer fancy 106 character*132 flags 107c 108c Data structures used for auto-orientation 109 real Rr(3,3), U(4,4), Xmom(5) 110c 111 COMMON /ASSCOM/ assout, verbose 112 integer assout 113 logical verbose 114c 115 real quadric(10), anisou(6) 116 real eigens(4), evecs(4,4), evecinv(4,4), evecit(4,4) 117 real qq(4,4), qp(4,4), temp(4,4) 118c 119 external anitoquad 120 integer anitoquad 121c 122 real problevel(50) 123c 124 real start(3),end(3) 125 real MARGIN 126 parameter (MARGIN = 1.15) 127c 128 real Uprin(3), Umean, Usigma, anisotropy, ellipticity 129 integer histogram(20), hislun 130 real anisi(MAXATM), sum_a, sum_a2, anis_mean, anis_sigma 131 real Biso(MAXATM), sum_b, sum_b2, sum_ab 132 integer nanis, niso, nhyd, nonpos 133 logical hwhacky 134 integer nerrors 135c 136 character*2 atomtype 137 integer natype 138c 139 real*8 wsum, xsum, ysum, zsum 140 real*8 xcom, ycom, zcom 141 real adist(110) 142 integer hdist(110), comlun, dshells 143c 144c Support for validation of similarity of bonded atoms 145 logical suvflag, cnflag 146 integer suvlun, suvbad, cnlun, cnbad, ccuijpair 147 character*1 prevchain 148 character*10 name_i, name_j 149 real anisov(6) 150 real cn_min_ccuij 151c 152c Default to CPK colors and VDW radii 153 integer DEFCOLS 154 parameter (DEFCOLS = 17) 155 character*60 defcol(DEFCOLS) 156 data defcol / 157 & 'COLOUR###### CA ############## 0.175 0.175 0.175 1.70', 158 & 'COLOUR###### C ############## 0.175 0.175 0.175 1.70', 159 & 'COLOUR#######C################ 0.625 0.625 0.625 1.70', 160 & 'COLOUR#######N################ 0.125 0.125 1.000 1.60', 161 & 'COLOUR#######O################ 0.750 0.050 0.050 1.50', 162 & 'COLOUR#######S################ 1.000 1.000 0.025 1.85', 163 & 'COLOUR#######H################ 1.000 1.000 1.000 1.20', 164 & 'COLOUR#######P################ 0.050 0.750 0.050 1.80', 165 & 'COLOUR##### CA ############### 0.175 0.175 0.175 1.70', 166 & 'COLOUR##### C ############### 0.175 0.175 0.175 1.70', 167 & 'COLOUR######C################# 0.625 0.625 0.625 1.70', 168 & 'COLOUR######N################# 0.125 0.125 1.000 1.60', 169 & 'COLOUR######O################# 0.750 0.050 0.050 1.50', 170 & 'COLOUR######S################# 1.000 1.000 0.025 1.85', 171 & 'COLOUR######H################# 1.000 1.000 1.000 1.20', 172 & 'COLOUR######P################# 0.050 0.750 0.050 1.80', 173 & 'COLOUR######################## 1.000 0.000 1.000 2.00' 174 & / 175c 176c Critical values for probability ellipsoids of a trivariate normal 177c distribution. From Table 6.1 of ORTEP-III manual (Oak Ridge National 178c Laboratory Report ORNL-6895, 1996). Tabulated below in increments of 179c 2% in probability. Default contours enclose a probability level of 180c 50% (critical value 1.5382). 181c 182 data problevel / 0.4299, 0.5479, 0.6334, 0.7035, 0.7644, 183 & 0.8192, 0.8694, 0.9162, 0.9605, 1.0026, 184 & 1.0430, 1.0821, 1.1200, 1.1570, 1.1932, 185 & 1.2288, 1.2638, 1.2985, 1.3330, 1.3672, 186 & 1.4013, 1.4354, 1.4695, 1.5037, 1.5382, 187 & 1.5729, 1.6080, 1.6436, 1.6797, 1.7164, 188 & 1.7540, 1.7924, 1.8318, 1.8724, 1.9144, 189 & 1.9580, 2.0034, 2.0510, 2.1012, 2.1544, 190 & 2.2114, 2.2730, 2.3404, 2.4153, 2.5003, 191 & 2.5997, 2.7216, 2.8829, 3.1365, 6.0000 / 192c 193c 194 assout = noise 195 verbose = .false. 196 hflag = .false. 197 acflag = .false. 198 bcflag = .false. 199 tflag = .false. 200 atflag = .false. 201 comflag = .false. 202 suvflag = .false. 203 cnflag = .false. 204 ellipses = .true. 205 hwhacky = .false. 206 nohydro = .false. 207 mini = .false. 208 auto = .false. 209 aniflag = .false. 210 fancy = 0 211 prob = 0.50 212 radius = 0.10 213 nerrors = 0 214c 215 prevchain = '@' 216c 217c Gfortran is nuts 218 wsum = 0 219 xsum = 0 220 ysum = 0 221 zsum = 0 222 sum_a = 0 223 sum_b = 0 224 do i=1,20 225 histogram(i) = 0 226 enddo 227 do i=1,110 228 adist(i) = 0 229 hdist(i) = 0 230 enddo 231 do i = 1,3 232 do j = 1,3 233 Rr(i,j) = 0 234 enddo 235 enddo 236 do i = 1,4 237 eigens(i) = 0 238 do j = 1,4 239 evecs(i,j) = 0 240 evecinv(i,j) = 0 241 evecit(i,j) = 0 242 qq(i,j) = 0 243 qp(i,j) = 0 244 temp(i,j) = 0 245 U(i,j) = 0 246 enddo 247 enddo 248 do i = 1,5 249 Xmom(i) = 0 250 enddo 251c 252 narg = iargc() 253 i = 1 254 5 continue 255 call getarg( i, flags ) 256 if (flags(1:5) .eq. '-help') goto 701 257 if (flags(1:6) .eq. '-debug') verbose = .true. 258 if (flags(1:2) .eq. '-h') hflag = .true. 259 if (flags(1:4) .eq. '-iso') ellipses = .false. 260 if (flags(1:6) .eq. '-aniso') aniflag = .true. 261 if (flags(1:4) .eq. '-rad') then 262 i = i + 1 263 if (i.gt.narg) goto 701 264 call getarg( i, flags ) 265 read (flags,*,err=701) radius 266 if (radius.lt.0) radius = 0.0 267 end if 268 if (flags(1:4) .eq. '-pro') then 269 i = i + 1 270 if (i.gt.narg) goto 701 271 call getarg( i, flags ) 272 read (flags,*,err=701) prob 273 if (prob.le.0.) stop 'illegal probability level' 274* If prob > 1 assume they meant it in percent 275 if (prob.gt.1.) prob = prob / 100. 276 end if 277 if (flags(1:5) .eq. '-Acol') then 278 acflag = .true. 279 bcflag = .false. 280 end if 281 if (flags(1:5) .eq. '-Bcol') then 282 bcflag = .true. 283 acflag = .false. 284 i = i + 1 285 if (i.gt.narg) goto 701 286 call getarg( i, flags ) 287 read (flags,*,err=701) Bmin 288 i = i + 1 289 if (i.gt.narg) goto 701 290 call getarg( i, flags ) 291 read (flags,*,err=701) Bmax 292 endif 293 if (flags(1:6) .eq. '-fancy') then 294 if (flags(7:7).eq.'0') then 295 fancy = 0 296 else if (flags(7:7).eq.'1') then 297 fancy = 1 298 else if (flags(7:7).eq.'2') then 299 fancy = 2 300 else if (flags(7:7).eq.'3') then 301 fancy = 3 302 else if (flags(7:7).eq.'4') then 303 fancy = 4 304 else if (flags(7:7).eq.'5') then 305 fancy = 5 306 else if (flags(7:7).eq.'6') then 307 fancy = 6 308 else 309 fancy = 1 310 endif 311 endif 312 if (flags(1:4) .eq. '-tab') then 313 tflag = .true. 314 hflag = .true. 315 acflag = .false. 316 bcflag = .false. 317 fancy = 0 318 do j=1,MAXATM 319 anisi(j) = 0.0 320 end do 321 hislun = NOISE 322 if (i.ge.narg) goto 799 323 call getarg(i+1,flags) 324 if (flags(1:1) .ne. '-') then 325 hislun = 1 326 open(unit=hislun,file=flags,status='UNKNOWN' 327 & ) 328 i = i + 1 329 endif 330 endif 331 if (flags(1:4) .eq. '-com') then 332 comflag = .true. 333 comlun = NOISE 334 if (i.ge.narg) goto 799 335 call getarg(i+1,flags) 336 if (flags(1:1) .ne. '-') then 337 comlun = 2 338 open(unit=comlun,file=flags,status='UNKNOWN' 339 & ) 340 i = i + 1 341 endif 342 endif 343 if (flags(1:4) .eq. '-suv') then 344 suvflag = .true. 345 suvlun = NOISE 346 suvlimit = 0.975 347 if (i.ge.narg) goto 799 348 call getarg(i+1,flags) 349 if (flags(1:1) .ne. '-') then 350 read (flags,*,err=701) suvlimit 351 i = i + 1 352 endif 353 endif 354 if (flags(1:9) .eq. '-cn_check') then 355 cnflag = .true. 356 cnlun = NOISE 357 cn_min_ccuij = 0.95 358 if (i.ge.narg) goto 799 359 call getarg(i+1,flags) 360 if (flags(1:1) .ne. '-') then 361 cnlun = 3 362 open(unit=cnlun,file=flags,status='UNKNOWN') 363 i = i + 1 364 endif 365 endif 366 if (flags(1:8) .eq. '-by_atom') then 367 atflag = .true. 368 endif 369 if (flags(1:8) .eq. '-nohydro') then 370 nohydro = .true. 371 endif 372 if (flags(1:8) .eq. '-mini') then 373 mini = .true. 374 auto = .true. 375 endif 376 if (flags(1:8) .eq. '-auto') then 377 auto = .true. 378 endif 379c 380 i = i + 1 381 if (i.le.narg) goto 5 382 goto 799 383 701 continue 384 write (noise,*) 'Raster3D Thermal Ellipsoid Program ', 385 & VERSION 386 write (noise,'(/,A)') 'syntax:' 387 write (noise,'(A)') 388 & 'rastep [-h] [-iso] [-Bcolor Bmin Bmax] [-prob Plevel]' 389 write (noise,'(A)') 390 & ' [-fancy[0-6]] [-radius R] [-auto]' 391 write (noise,'(A,A)') 392 & ' [-nohydrogens] [-suv [suv_limit]] [-cn_check [cnfile]]' 393 write (noise,'(A,A)') 394 & ' [-tabulate [tabfile]] [-by_atomtype] [-com [comfile]]' 395 call exit(-1) 396 799 continue 397 398c 399c Critical values for the radius corresponding to a sphere 400c enclosing the requested probability level are taken from 401c Table 6.1 of the ORTEP manual 402 iprob = (prob+0.01)*50. 403 pradius = problevel(iprob) 404 405c 406 write (noise,800) 407 write (noise,*) 'Raster3D Thermal Ellipsoid Program ', 408 & VERSION 409 write (noise,*) 'E A Merritt - 11 Mar 2010' 410 write (noise,800) 411 800 format('************************************************') 412c 413 if (.not.ellipses) then 414 write (noise,801) float(iprob)/50. 415 801 format(' Spheres will bound Biso probability level', f5.2) 416 else 417 write (noise,802) float(iprob)/50. 418 802 format(' Ellipsoids will bound probability level', f5.2) 419 endif 420 write (noise,803) pradius 421 803 format(' Corresponding critical value ', f7.4) 422 if (aniflag) write (noise,*) 423 & ' Isotropic atoms will be included in the analysis' 424c 425 if (acflag) then 426 write (noise,*) 'Atoms will be colored based on Anisotropy' 427 endif 428c 429 if (bcflag) then 430 write (noise,*) 'Atom colors will be assigned based on Biso' 431 write (noise,*) ' from dark blue = Bmin =', Bmin 432 write (noise,*) ' to light red = Bmax =', Bmax 433 Umin = Bmin / (8. * 3.14159*3.14159) 434 Umax = Bmax / (8. * 3.14159*3.14159) 435 Umin = Umin 436 Umax = Umax 437 endif 438c 439 if (.not. hflag) then 440 WRITE(OUTPUT,'(A,A,I5,A)') 441 & 'Raster3D thermal ellipsoid program ',VERSION, 442 & INT(prob*100.+0.5), '% probability bounds' 443 if (mini) then 444 WRITE(OUTPUT,'(A)') '22 26 tiles in x,y' 445 else 446 WRITE(OUTPUT,'(A)') '80 64 tiles in x,y' 447 endif 448 WRITE(OUTPUT,'(A)') ' 8 8 pixels (x,y) per tile' 449 WRITE(OUTPUT,'(A)') '4 3x3 virtual pixels -> 2x2 pixels' 450 WRITE(OUTPUT,'(A)') '1 1 1 white background' 451 WRITE(OUTPUT,'(A)') 'F no, shadows are dorky' 452 WRITE(OUTPUT,'(A)') '25 Phong power' 453 WRITE(OUTPUT,'(A)') '0.15 secondary light contribution' 454 WRITE(OUTPUT,'(A)') '0.05 ambient light contribution' 455 WRITE(OUTPUT,'(A)') '0.25 specular reflection component' 456 WRITE(OUTPUT,'(A)') '0.0 No perspective' 457 WRITE(OUTPUT,'(A)') '1 1 1 main light source position' 458 end if 459 460c 461 if (auto) then 462 ASPECT = 22./26. 463 else 464 ASPECT = 80./64. 465 endif 466c 467 NCOL = 0 468 NATM = 0 469 NANI = 0 470 nanis = 0 471 niso = 0 472 nhyd = 0 473 nonpos = 0 47410 CONTINUE 475 READ(INPUT,'(A80)',END=50) CARD 476 IF (CARD(1:4).EQ.'COLO') THEN 477 NCOL = NCOL + 1 478 IF (NCOL.GT.MAXCOL) THEN 479 WRITE(NOISE,*) 'Colour table overflow. Increase ', 480 & 'MAXCOL and recompile.' 481 STOP 10 482 ENDIF 483 READ(CARD,'(6X,A24,3F8.3,F6.2)') MASK(NCOL), 484 & (RGB(I,NCOL),I=1,3), VDW(NCOL) 485 ELSEIF (nohydro .AND. CARD(77:78).EQ.' H') THEN 486 goto 10 487 ELSEIF (CARD(1:6).EQ.'ANISOU') THEN 488 nani = nani + 1 489 if (card(13:27).ne.atom(natm)(13:27)) goto 14 490 read (card(29:70),*,err=12,end=12) (uij(i,natm),i=1,6) 491 do i=1,6 492 uij(i,natm) = uij(i,natm) * 0.0001 493 enddo 494 noerr = 1 495 do i=1,3 496 if (uij(i,natm).eq.0.0) then 497 uij(i,natm) = 0.0001 498 noerr = 0 499 endif 500 enddo 501 if (noerr.eq.0) goto 15 502 ELSEIF (CARD(1:4).EQ.'ATOM'.OR.CARD(1:4).EQ.'HETA') THEN 503 NATM = NATM + 1 504 IF (NATM.GT.MAXATM) THEN 505 WRITE(NOISE,*) 'Atom array overflow. Increase ', 506 & 'MAXATM and recompile.' 507 STOP 20 508 ENDIF 509 ATOM(NATM) = CARD 510 uij(1,natm) = -1.0 511 ELSEIF (CARD(1:3).EQ.'END') THEN 512 GO TO 50 513 ENDIF 514 GO TO 10 51512 write(noise,*) '*** Format problem - ', card(13:70) 516 nerrors = nerrors + 1 517 goto 10 51814 write(noise,*) '*** ANISOU record out of order - ', card(13:70) 519 nerrors = nerrors + 1 520 goto 10 52115 write(noise,*) '*** Illegal ANISOU values - ', card(13:70) 522 nerrors = nerrors + 1 523 goto 10 524* Come here when EOF or 'END' record is reached 52550 CONTINUE 526 IF (NATM.EQ.0) THEN 527 WRITE(NOISE,*) 'No atoms in input.' 528 STOP 30 529 ENDIF 530* Load default colors after any that were read in 531 IF (NCOL.LE.MAXCOL-DEFCOLS) THEN 532 DO i = 1,DEFCOLS 533 NCOL = NCOL + 1 534 READ(defcol(i),'(6X,A24,3F8.3,F6.2)') MASK(NCOL), 535 & (RGB(J,NCOL),J=1,3), VDW(NCOL) 536 ENDDO 537 ENDIF 538* 539 IF (NCOL.EQ.0) THEN 540 WRITE(NOISE,*) 'No colours in input.' 541 STOP 40 542 ENDIF 543 XMAX = -1E20 544 XMIN = 1E20 545 YMAX = -1E20 546 YMIN = 1E20 547 ZMAX = -1E20 548 ZMIN = 1E20 549 DO 100 IATM=1,NATM 550 CARD = ATOM(IATM) 551c 552c Do a little pre-processing to make later bookkeeping easier 553c At least screen out non-conformant PDB files that contain 554c something obviously not an element type in columns 77:78 555c Hydrogen naming conventions are totally messed up. 556c 557 if (atflag) then 558 resname = card(18:20) 559c if (resname.eq.'MET' .and. card(14:15).eq.'SD') then 560c atom(iatm)(77:78) = 'SD' 561c end if 562 if ( resname.eq.'HOH' .or. resname.eq.'H2O' 563 & .or. resname.eq.'WAT') then 564 atom(iatm)(77:78) = 'OW' 565 end if 566 if (atom(iatm)(78:78).ge.'0' .and. atom(iatm)(78:78).le.'9') 567 & atom(iatm)(77:78) = ' ' 568 end if 569 if (nohydro .and. atom(iatm)(77:78).eq.' ') then 570 do 70 i = 13,16 571 if (atom(iatm)(i:i).eq.' ') goto 70 572 if (atom(iatm)(i:i).ge.'1' 573 & .and. atom(iatm)(i:i).le.'4') goto 70 574 if (atom(iatm)(i:i).ne.'H') goto 71 575 atom(iatm)(77:78) = ' H' 576 70 continue 577 71 continue 578 end if 579 580c 581 TEST = CARD(7:30) 582 DO 80 ICOL=1,NCOL 583 IF (MATCH(TEST,MASK(ICOL))) THEN 584 READ(CARD,'(30X,3F8.3,6X,F6.2)',end=82,err=82) 585 & X,Y,Z, Biso(IATM) 586 IF (Biso(IATM).LE.0.0) THEN 587 nerrors = nerrors + 1 588 write(noise,*) '*** Illegal Biso ',Biso(IATM),' - ', 589 & atom(iatm)(13:27) 590 Biso(IATM) = 0.0 591 ENDIF 592 Uiso = Biso(IATM) / (8. * 3.14159*3.14159) 593 SPAM(1,IATM) = X 594 SPAM(2,IATM) = Y 595 SPAM(3,IATM) = Z 596 SPAM(4,IATM) = Uiso 597 SPAM(5,IATM) = ICOL 598 RAD = sqrt(Uiso) * PRADIUS 599 XMAX = MAX(XMAX,X+RAD) 600 XMIN = MIN(XMIN,X-RAD) 601 YMAX = MAX(YMAX,Y+RAD) 602 YMIN = MIN(YMIN,Y-RAD) 603 ZMAX = MAX(ZMAX,Z+RAD) 604 ZMIN = MIN(ZMIN,Z-RAD) 605c atomtype = CARD(77:78) 606c if (atomtype.ne.' ') then 607c weight = amass(atomtype) 608c else 609 weight = 13.4 610 wsum = wsum + weight 611 xsum = xsum + weight * X 612 ysum = ysum + weight * Y 613 zsum = zsum + weight * Z 614 if (uij(1,iatm).lt.0 .and. aniflag) then 615 uij(1,iatm) = Uiso 616 uij(2,iatm) = Uiso 617 uij(3,iatm) = Uiso 618 uij(4,iatm) = 0.0 619 uij(5,iatm) = 0.0 620 uij(6,iatm) = 0.0 621c write(0,*) 'Setting Uiso for atom ',card 622 endif 623 GO TO 100 624 ENDIF 62580 CONTINUE 626 WRITE(NOISE,*) 'No colour table mask matches this atom:' 627 WRITE(NOISE,*) ATOM(IATM) 628 STOP 90 62982 continue 630 write(noise,*) 'Input format problem in record' 631 write(noise,*) CARD 632 STOP 90 633100 CONTINUE 634 XMID = (XMAX+XMIN)/2. 635 YMID = (YMAX+YMIN)/2. 636 ZMID = (ZMAX+ZMIN)/2. 637 TX = -XMID 638 TY = -YMID 639 TZ = -ZMID 640 xcom = xsum / wsum 641 ycom = ysum / wsum 642 zcom = zsum / wsum 643 IF (ASPECT.GE.1.) THEN 644* The X direction is wider than the Y 645 XROOM = ASPECT 646 YROOM = 1. 647 ZROOM = 2. 648 ELSE 649 XROOM = 1. 650 YROOM = ASPECT 651 ZROOM = 2. 652 ENDIF 653 XSPAN = XMAX-XMIN 654 YSPAN = YMAX-YMIN 655 ZSPAN = ZMAX-ZMIN 656 SCALE = MAX(XSPAN/XROOM,YSPAN/YROOM,ZSPAN/ZROOM) 657* Leave a little extra room as a border: 658 SCALE = SCALE / 0.90 659 if (mini) then 660 scale = sqrt(xspan**2+yspan**2+zspan**2) * aspect 661 if (scale .lt. 9.) scale = 9. 662 end if 663* 664* These are for the center-of-mass table 665 DMAX = MAX( ABS(XMAX-XCOM), ABS(XMIN-XCOM) )**2 666 & + MAX( ABS(YMAX-YCOM), ABS(YMIN-YCOM) )**2 667 & + MAX( ABS(ZMAX-ZCOM), ABS(ZMIN-ZCOM) )**2 668 DMAX = SQRT(DMAX) 669 dshells = 100 670* 671 if (.not. hflag) then 672 WRITE(OUTPUT,120) TX,TY,TZ,SCALE 673120 FORMAT('1 0 0 0 input co-ordinate + radius transformation'/ 674 & '0 1 0 0'/ 675 & '0 0 1 0'/ 676 & 4F10.3) 677 WRITE(OUTPUT,'(A)') '3 mixed object types' 678 WRITE(OUTPUT,'(A)') '*' 679 WRITE(OUTPUT,'(A)') '*' 680 WRITE(OUTPUT,'(A)') '*' 681 end if 682c 683c Auto-orientation 684c 25-Aug-1999 685c Find Eigenvectors of moment of inertia tensor. 686c Arrange smallest Eigenvalue along Y, largest along Z. 687c Problems: 688c - Could emphasize side-chain over backbone by ignoring 689c atoms O and N, but in practice this doesn't seem to help much. 690c - Scaling is wrong, because it was done before rotation. 691c 692 if (auto) then 693 do 125 iatm = 1, natm 694 if (atom(iatm)(1:4).ne.'ATOM' 695 & .and. atom(iatm)(1:4).ne.'HETA') goto 125 696C if (atom(iatm)(13:15).eq.' O ') goto 125 697C if (atom(iatm)(13:15).eq.' N ') goto 125 698 if (atom(iatm)(77:78).eq.' H' .and. nohydro) goto 125 699 x = spam(1,iatm) - xcom 700 y = spam(2,iatm) - ycom 701 z = spam(3,iatm) - zcom 702 Rq = (x*x + y*y + z*z) 703 Rv = sqrt(Rq) 704 Rr(1,1) = Rr(1,1) + x*x 705 Rr(1,2) = Rr(1,2) + x*y 706 Rr(1,3) = Rr(1,3) + x*z 707 Rr(2,1) = Rr(2,1) + y*x 708 Rr(2,2) = Rr(2,2) + y*y 709 Rr(2,3) = Rr(2,3) + y*z 710 Rr(3,1) = Rr(3,1) + z*x 711 Rr(3,2) = Rr(3,2) + z*y 712 Rr(3,3) = Rr(3,3) + z*z 713 Xmom(1) = Xmom(1) + 1. 714 Xmom(2) = Xmom(2) + Rv 715 Xmom(3) = Xmom(3) + Rq 716 Xmom(4) = Xmom(4) + Rv*Rq 717 Xmom(5) = Xmom(5) + Rq*Rq 718 125 continue 719 if (verbose) then 720 write (NOISE,'(/,A,5G13.5)') ' Radial moments:',Xmom 721 write (NOISE,'(A)') ' Moment of inertia tensor' 722 end if 723 Rg = sqrt(Xmom(3)/Xmom(1)) 724 U(1,1) = Xmom(3) 725 U(2,2) = Xmom(3) 726 U(3,3) = Xmom(3) 727 do k = 1,3 728 do l = 1,3 729 U(k,l) = (U(k,l) - Rr(k,l)) / Xmom(1) 730 end do 731 if (verbose) write (NOISE,'(3G13.6)') (U(k,l),l=1,3) 732 end do 733 call jacobi( U, 3, 4, Eigens, Evecs ) 734c Re-order so that long axis is vertical, short axis towards viewer 735 kmax = 1 736 if (Eigens(2).gt.Eigens(kmax)) kmax = 2 737 if (Eigens(3).gt.Eigens(kmax)) kmax = 3 738 kmin = 1 739 if (Eigens(2).le.Eigens(kmin)) kmin = 2 740 if (Eigens(3).le.Eigens(kmin)) kmin = 3 741 kmid = 6 - (kmin + kmax) 742 if (verbose) then 743 write (NOISE,'(A,/,3F13.5,10X,3i2)') ' Eigenvalues:', 744 & Eigens(kmin),Eigens(kmid),Eigens(kmax) 745 & ,kmin,kmid,kmax 746 write (NOISE,'(A)') ' Eigenvectors:' 747 do k = 1,3 748 write (NOISE,'(3F13.5)') 749 & Evecs(k,kmin),Evecs(k,kmid),Evecs(k,kmax) 750 enddo 751 end if 752 do k = 1,3 753 U(k,1) = Evecs(k,kmid) 754 U(k,2) = Evecs(k,kmin) 755 U(k,3) = Evecs(k,kmax) 756 U(k,4) = 0.0 757 U(4,k) = 0.0 758 enddo 759 U(4,4) = 1.0 760c Beware! may be left-handed at this point! 761 det = tinv4( evecinv, U ) 762 if (det .lt. 0.0) then 763 do k = 1,3 764 U(k,3) = -U(k,3) 765 enddo 766 endif 767c OK, now it should be right-handed 768 WRITE(OUTPUT,'(A)') '# Auto-orientation matrix' 769 WRITE(OUTPUT,'(A)') '16' 770 WRITE(OUTPUT,'(A)') 'ROTATION' 771 do k = 1,3 772 write (OUTPUT,'(3F13.5)') U(1,k),U(2,k),U(3,k) 773 enddo 774 WRITE(OUTPUT,'(A)') '# End auto-orientation' 775 end if 776c 777c Label output records 778c 779 if (.not. tflag) then 780 WRITE(OUTPUT,'(A,A)') 781 & '# Thermal ellipsoids from Rastep Version ',VERSION 782 WRITE(OUTPUT,'(A,F5.2)') '# Probability level',float(iprob)/50. 783 end if 784c 785c Write ellipsoids to input file for render 786c 787 IF (fancy.eq.0 .and. .not.tflag) GOTO 139 788c 789c First, optional pass, to write fancy stuff associated with ellipsoids 790c 791 IATM = 1 792 130 CONTINUE 793 IF (ATOM(IATM)(1:4).EQ.'ATOM' .OR. 794 & ATOM(IATM)(1:4).EQ.'HETA') THEN 795 X = SPAM(1,IATM) 796 Y = SPAM(2,IATM) 797 Z = SPAM(3,IATM) 798 ICOL = SPAM(5,IATM) 799 if (bcflag) then 800 call U2RGB( SPAM(4,IATM), Umin, Umax, RED, GREEN, BLUE ) 801 RED = RED*RED 802 GREEN = GREEN*GREEN 803 BLUE = BLUE*BLUE 804 else if (acflag) then 805 call A2RGB( 1.0, RED, GREEN, BLUE ) 806 RED = RED*RED 807 GREEN = GREEN*GREEN 808 BLUE = BLUE*BLUE 809 else 810 RED = RGB(1,ICOL) 811 GREEN = RGB(2,ICOL) 812 BLUE = RGB(3,ICOL) 813 endif 814 IF (ellipses .and. uij(1,iatm).ge.0.) THEN 815 do i=1,6 816 anisou(i) = uij(i,iatm) 817 enddo 818 if (anitoquad(anisou,pradius,quadric,eigens,evecs).lt.0)then 819 write(noise,*) '*** Non-positive definite ellipsoid - ', 820 & atom(iatm)(13:27) 821 nonpos = nonpos + 1 822 nerrors = nerrors + 1 823 Biso(iatm) = 0.0 824 goto 138 825 endif 826 goto 132 827 132 continue 828 radlim = pradius * max( eigens(1),eigens(2),eigens(3) ) 829 radlim = radlim * MARGIN 830c 831c Only for debugging ellipsoids 832 if (verbose) then 833 write (noise,901) 'ANISOU ',X,Y,Z,ANISOU 834 write (noise,902) 'QUADRIC',QUADRIC 835 write (noise,903) 'Eigenvalues', (EIGENS(i),i=1,3), 836 & 'prob', prob,'limiting radius', radlim 837 write (noise,904) 'Evecs ',((evecs(i,j),i=1,3),j=1,3) 838 endif 839901 format(a,3f8.3,6f8.4) 840902 format(a,10f8.3) 841903 format(a,3f8.3,4x,a,f8.3,4x,a,f8.3) 842904 format(a,9f7.3) 843c 844c Tabulate principal axes of ellipsoid for each atom 845 if (tflag) then 846 do i=1,3 847 Uprin(i) = eigens(i)**2 848 enddo 849 if (Uprin(2).gt.Uprin(1)) then 850 Umean = Uprin(1) 851 Uprin(1) = Uprin(2) 852 Uprin(2) = Umean 853 endif 854 if (Uprin(3).gt.Uprin(1)) then 855 Umean = Uprin(1) 856 Uprin(1) = Uprin(3) 857 Uprin(3) = Umean 858 endif 859 if (Uprin(3).gt.Uprin(2)) then 860 Umean = Uprin(2) 861 Uprin(2) = Uprin(3) 862 Uprin(3) = Umean 863 endif 864c 865c Anisotropy we define as Umin / Umax 866c as in shelxpro output 867 anisotropy = min(Uprin(1),Uprin(2),Uprin(3)) 868 & / max(Uprin(1),Uprin(2),Uprin(3)) 869c 870c But don't count atoms which are perfectly isotropic 871 if (atom(iatm)(77:78).eq.' H') then 872 nhyd = nhyd + 1 873 if (anisotropy .ne. 1.0) hwhacky = .true. 874 else if (anisotropy .eq. 1.0) then 875 niso = niso + 1 876 else 877 anisi(iatm) = anisotropy 878 sum_A = sum_A + anisotropy 879 sum_B = sum_B + Biso(iatm) 880 nanis = nanis + 1 881 end if 882c 883c Ellipticity we define as 1 / anisotropy 884 if (anisotropy.eq.0) then 885 ellipticity = 0 886 else 887 ellipticity = 1. / anisotropy 888 end if 889c 890c Longhi et al (1997) JMB 268, 779-799. 891c proposed another measure A = sigU / meanU 892 Umean = (Uprin(1) + Uprin(2) + Uprin(3)) / 3.0 893 Usigma = (Uprin(1)-Umean)**2 894 & + (Uprin(2)-Umean)**2 + (Uprin(3)-Umean)**2 895 Usigma = sqrt( Usigma ) / 3.0 896 alonghi = Usigma / Umean 897c 898c Might want to check correlation with Uiso 899 Uiso = SPAM(4,iatm) 900c 901c Cosmetic changes to atom identifier for the sake of sorting 902c We will force there to be exactly three entities printed. 903c PDB format is just a mess: 904c cols 13:16 atom 905c col 17 alternate conf 906c cols 18:20 residue 907c col 22 chain 908c cols 23:27 resnum 909c 910 do i = 16, 13, -1 911 if (ATOM(iatm)(i:i) .ne. ' ') j = i 912 enddo 913 do i = 13, 17 914 if (ATOM(iatm)(i:i) .ne. ' ') k = i 915 enddo 916 do i = j, k 917 if (ATOM(iatm)(i:i) .eq. ' ') ATOM(iatm)(i:i) = '_' 918 enddo 919c if (ATOM(iatm)(17:17) .ne. ' ') then 920c do i = 18,19 921c if (ATOM(iatm)(i:i) .eq. ' ') ATOM(iatm)(i:i) = '_' 922c enddo 923c endif 924 spacer = ' ' 925 if (ATOM(iatm)(22:22) .ne. ' ') then 926 spacer = '_' 927 do i = 23,25 928 if (ATOM(iatm)(i:i) .eq. ' ') ATOM(iatm)(i:i) = '_' 929 enddo 930 endif 931 write (output,905) ATOM(iatm)(13:17), 932 & ATOM(iatm)(18:22),spacer,ATOM(iatm)(23:27), 933 & Uprin(1),Uprin(2),Uprin(3),anisotropy,Uiso 934905 format(A5,1X,A5,A1,A5,3F9.4,2X,F9.4,F9.4,F9.4) 935c 936c Also make a histogram of anisotropies 937 i = anisotropy * 20. + 1 938 histogram(i) = histogram(i) + 1 939c 940c And a table of <anis> by distance from center of mass 941 if (comflag .and. anisotropy.lt.1.0) then 942 dist = sqrt( (SPAM(1,iatm)-XCOM)**2 943 & + (SPAM(2,iatm)-YCOM)**2 944 & + (SPAM(3,iatm)-ZCOM)**2 ) 945 i = (dist/dmax) * float(dshells) + 1 946 adist(i) = adist(i) + anisotropy 947 hdist(i) = hdist(i) + 1 948 endif 949c 950c End tabulation code 951 endif 952 953c 954c Skip hydrogens if requested 955 if (nohydro .and. ATOM(IATM)(77:78).eq.' H') goto 138 956 957c 958c Draw principal axes inside bounding ellipsoid 959 if (fancy.eq.1) then 960 do i=1,3 961 size = eigens(i) * pradius - 0.02 962 start(1) = x - size*evecs(1,i) 963 start(2) = y - size*evecs(2,i) 964 start(3) = z - size*evecs(3,i) 965 end(1) = x + size*evecs(1,i) 966 end(2) = y + size*evecs(2,i) 967 end(3) = z + size*evecs(3,i) 968 write (output,907) start, end 969907 format(' 3',/, 970 & 3f9.3,' 0.02',3f9.3,' 0.02',' 0.5 1.0 0.3') 971 enddo 972 endif 973c 974c Draw longest principle axis only 975c (experimental use only - not supported or documented) 976 if (fancy.eq.4) then 977 imax = 1 978 if (eigens(2).gt.eigens(imax)) imax = 2 979 if (eigens(3).gt.eigens(imax)) imax = 3 980 size = eigens(imax) * pradius 981 imin = 1 982 if (eigens(2).lt.eigens(imin)) imin = 2 983 if (eigens(3).lt.eigens(imin)) imin = 3 984 imed = 6 - (imax+imin) 985 size = size * eigens(imax)/eigens(imed) 986 start(1) = x - size*evecs(1,imax) 987 start(2) = y - size*evecs(2,imax) 988 start(3) = z - size*evecs(3,imax) 989 end(1) = x + size*evecs(1,imax) 990 end(2) = y + size*evecs(2,imax) 991 end(3) = z + size*evecs(3,imax) 992 write (output,907) start, end 993 endif 994c 995c Construct 3 quadrics corresponding to the 3 orthogonal planes 996c through the center of our ellipsoid 997 if (fancy.eq.2 .or. fancy.eq.3) then 998 eigens(4) = 1.0 999 evecs(4,4)= 1.0 1000 det = tinv4( evecinv, evecs ) 1001 call trnsp4( evecit, evecinv ) 1002 do k = 1,3 1003 do i = 1,4 1004 do j = 1,4 1005 QQ(i,j) = 0.0 1006 enddo 1007 QQ(i,i) = 1. / (eigens(i)*eigens(i)) 1008 enddo 1009 QQ(k,k) = 1000. 1010 QQ(4,4) = -pradius*pradius 1011 call tmul4( TEMP, QQ, evecinv ) 1012 call tmul4( QP, evecit, TEMP ) 1013 if (acflag) then 1014 call A2RGB( anisotropy, red, green, blue ) 1015 red = red*red 1016 green = green*green 1017 blue = blue*blue 1018 endif 1019 write (output,151) 14, X,Y,Z, radlim, red, green, blue 1020 write (output,152) QP(1,1),QP(2,2),QP(3,3),QP(1,2),QP(2,3), 1021 & QP(1,3),QP(1,4),QP(2,4),QP(3,4),QP(4,4) 1022 enddo 1023 endif 1024c 1025 endif 1026 ENDIF 1027 138 continue 1028 IATM = IATM + 1 1029 IF (IATM.LE.NATM) GOTO 130 1030c 1031c Set transparency for enclosing ellipoids 1032 if (fancy.eq.1 .or. fancy.eq.3) then 1033 write (output,'(A,/,A)') '9 Begin transparent ellipsoids','8 ' 1034 write (output,'(A)') ' 15. 0.6 1.0 1.0 1.0 0.6 0 0 0 0' 1035 else if (fancy.eq.2 .or. fancy.eq.4) then 1036 goto 160 1037 endif 1038c 1039 139 CONTINUE 1040c 1041c If we're just tabulating ellipticity, start making tables 1042c 1043 if (tflag) then 1044 total = 0 1045 do i = 1,20 1046 total = total + histogram(i) 1047 end do 1048 if (total.eq.0) then 1049 write (noise,*) 'No ANISOU records found' 1050 call exit(-1) 1051 end if 1052 if (nanis.eq.0) then 1053 write (noise,*) 'No anisotropic atoms found' 1054 call exit(-1) 1055 end if 1056 if (hwhacky) then 1057 write(noise,*) 1058 & 'You seem to have anisotropic hydrogens', 1059 & ' - is this some kind of joke?' 1060 nerrors = nerrors + 1 1061 end if 1062 1063 write (hislun,'(A)') '# Anisotropy Fraction Number' 1064 write (hislun,'(A)') '# range of atoms of atoms' 1065 do i = 1,20 1066 write (hislun,'(2F5.2,3X,F8.3,I10)') 1067 & (float(i)-1.)/20., float(i)/20., 1068 & float(histogram(i))/total, histogram(i) 1069 end do 1070c Calculate mean and sigma of distribution 1071 sum_a2 = 0.0 1072 sum_b2 = 0.0 1073 sum_ab = 0.0 1074 anis_mean = sum_A / float(nanis) 1075 anis_sigma = 0.0 1076 biso_mean = sum_B / float(nanis) 1077 ccoef = 0.0 1078 do i = 1,NATM 1079 if (anisi(i).ne.0) then 1080 sum_a2 = sum_a2 + (anisi(i) - anis_mean)**2 1081 sum_b2 = sum_b2 + (Biso(i) - biso_mean)**2 1082 sum_ab = sum_ab 1083 & + (anisi(i)-anis_mean) * (Biso(i)-biso_mean) 1084 end if 1085 end do 1086 if (nanis.gt.1) then 1087 anis_sigma = sqrt( sum_a2 / float(nanis-1) ) 1088 ccoef = sum_ab / sqrt(sum_a2 * sum_b2) 1089 end if 1090 1091 write (hislun,'(A)') '#' 1092 write (hislun,'(A,I10)')'# number of ANISOU records:',nani 1093 write (hislun,'(A,I10)')'# non-isotropic atoms:',nanis 1094 write (hislun,'(A,I10)')'# isotropic atoms:',niso 1095 if (nonpos.gt.0) 1096 & write (hislun,'(A,I10)')'# nonpositive-definite APDs:',nonpos 1097 if (nhyd.gt.0) 1098 & write (hislun,'(A,I10)')'# ANISOU hydrogens:',nhyd 1099 write (hislun,'(A)') '#' 1100 write (hislun,'(A,F7.3)') 1101 & '# correlation of anisotropy with B_iso:', ccoef 1102 write (hislun,'(A)') '#' 1103 write (hislun,'(A)') '# Anisotropy B_iso' 1104 write (hislun,'(A)') '# AtomType mean sigma mean number' 1105 write (hislun,'(A)') '# --------- ----------- ------ ------' 1106 write (hislun,'(A,F7.3,F7.3,F7.2,I8)') 1107 & '#| Total',anis_mean,anis_sigma,biso_mean,nanis 1108 1109c 1110c If we're tabulating ellipticity, but not by atom_type, then we're done 1111c 1112 if (.not. atflag) goto 145 1113 140 continue 1114c 1115C Find an atom type we haven't done yet, exit if none left 1116c This only works if the PDB file contains the atom type in 1117c columns 77:78 (standard since sometime in 1997, but many 1118c files do not conform to this standard) 1119c 1120 do i = 1, NATM 1121 if (anisi(i).ne.0.0 .and. atom(i)(77:78).ne.' ') then 1122 atomtype = atom(i)(77:78) 1123 goto 141 1124 end if 1125 end do 1126 goto 145 1127 141 continue 1128 1129 sum_A = 0.0 1130 sum_B = 0.0 1131 sum_a2 = 0.0 1132 biso_mean = 0.0 1133 natype = 0 1134c 1135c Loop over atoms looking for the right ones 1136c 1137 do i = 1, NATM 1138 if (anisi(i).ne.0.0 .and. atomtype.eq.atom(i)(77:78)) then 1139 natype = natype + 1 1140 sum_A = sum_A + anisi(i) 1141 end if 1142 end do 1143 anis_mean = sum_A / float(natype) 1144 anis_sigma = 0.0 1145 do i = 1, NATM 1146 if (anisi(i).ne.0.0 .and. atomtype.eq.atom(i)(77:78)) then 1147 sum_a2 = sum_a2 + (anisi(i)-anis_mean)**2 1148 biso_mean = biso_mean + Biso(i) 1149 atom(i)(77:78) = ' ' 1150 end if 1151 end do 1152 if (natype.gt.1) anis_sigma = sqrt( sum_a2 / float(natype-1) ) 1153 biso_mean = biso_mean / float(natype) 1154 1155 write (hislun,'(A,4X,A,F7.3,F7.3,F7.2,I8)') 1156 & '#| ',atomtype,anis_mean,anis_sigma,biso_mean,natype 1157 goto 140 1158 end if 1159 IATM = 1 1160 goto 150 1161c 1162c Write out plot of mean anisotropy as a function of distance from c.o.m. 1163c The hisclean routine is not strictly necessary; it collapses shells at 1164c the extreme ranges of distance that contain only a few atoms. 1165c 1166 145 continue 1167 if (comflag) then 1168 write (comlun,'(A/A,A/A/A)') '#', 1169 & '# Mean anisotropy as a function of distance', 1170 & ' from center of mass', 1171 & '#', 1172 & '# Distance <anis> #atoms' 1173 call hisclean( adist, hdist, dshells, nanis ) 1174 do i = 1, dshells 1175 if (hdist(i).gt.0) then 1176 adist(i) = adist(i) / float(hdist(i)) 1177 dmid = (dmax/float(dshells))*(float(i)-0.5) 1178 write (comlun,146) dmid, adist(i), hdist(i) 1179 endif 1180 enddo 1181 146 format(F10.3,F10.3,I10) 1182 endif 1183 if (suvflag.or.cnflag) goto 160 1184 call exit(0) 1185 1186c 1187c Second pass write a single sphere or ellipsoid for each atom 1188c 1189 150 CONTINUE 1190 IF (ATOM(IATM)(1:4).EQ.'ATOM' .OR. 1191 & ATOM(IATM)(1:4).EQ.'HETA') THEN 1192 if (nohydro .and. ATOM(IATM)(77:78).eq.' H') goto 154 1193 X = SPAM(1,IATM) 1194 Y = SPAM(2,IATM) 1195 Z = SPAM(3,IATM) 1196 ICOL = SPAM(5,IATM) 1197 if (bcflag) then 1198 call U2RGB( SPAM(4,IATM), Umin, Umax, RED, GREEN, BLUE ) 1199 RED = RED*RED 1200 GREEN = GREEN*GREEN 1201 BLUE = BLUE*BLUE 1202 else if (acflag) then 1203 call A2RGB( 1.0, RED, GREEN, BLUE ) 1204 RED = RED*RED 1205 GREEN = GREEN*GREEN 1206 BLUE = BLUE*BLUE 1207 else 1208 RED = RGB(1,ICOL) 1209 GREEN = RGB(2,ICOL) 1210 BLUE = RGB(3,ICOL) 1211 endif 1212 IF (.not. ellipses) THEN 1213 RAD = sqrt(SPAM(4,IATM)) * PRADIUS 1214 WRITE(OUTPUT,151) 2, X,Y,Z,RAD,RED,GREEN,BLUE 1215 ELSE IF (uij(1,iatm).gt.0.) THEN 1216 do i=1,6 1217 anisou(i) = uij(i,iatm) 1218 enddo 1219 if (anitoquad(anisou,pradius,quadric,eigens,evecs).lt.0)then 1220 write(noise,*) '*** Non-positive definite ellipsoid - ', 1221 & atom(iatm)(13:26) 1222 nerrors = nerrors + 1 1223 Biso(iatm) = 0.0 1224 red = 1.0 1225 green = 0.0 1226 blue = 1.0 1227 endif 1228 radlim = pradius * max( eigens(1),eigens(2),eigens(3) ) 1229 radlim = radlim * MARGIN 1230 if (acflag) then 1231 anisotropy = min( eigens(1),eigens(2),eigens(3) ) 1232 & / max( eigens(1),eigens(2),eigens(3) ) 1233 anisotropy = anisotropy * anisotropy 1234 call A2RGB( anisotropy, red, green, blue ) 1235 red = red*red 1236 green = green*green 1237 blue = blue*blue 1238 endif 1239 if (fancy.eq.5 .or. fancy.eq.6) then 1240 write (output,'(I2,/,A,I2,/,A)') 8, 1241 & ' -1.0 -1.0 -1.0 -1.0 -1.0 0.0 0 0 0', 1242 & fancy-1, 'ORTEP_LIKE' 1243 if (fancy.eq.6) write (output,171) 0.5, 0.5, 0.5 1244 write (output,172) 0, x,y,z, (evecs(i,1),i=1,3) 1245 write (output,172) 0, x,y,z, (evecs(i,2),i=1,3) 1246 write (output,172) 0, x,y,z, (evecs(i,3),i=1,3) 1247 171 format('BOUNDING_COLOR ',3F6.3) 1248 172 format('BOUNDING_PLANE ',I2,6F10.4) 1249 endif 1250 write (output,151) 14, x,y,z,radlim,red,green,blue 1251 write (output,152) (quadric(i),i=1,10) 1252 ELSE 1253 RAD = sqrt(SPAM(4,IATM)) * PRADIUS 1254 WRITE(OUTPUT,151) 2, X,Y,Z,RAD,RED,GREEN,BLUE 1255 151 FORMAT(I2,/,3(1X,F8.3),4F8.3) 1256 152 FORMAT(10F12.4) 1257 ENDIF 1258 goto 154 1259 ENDIF 1260 154 continue 1261 IATM = IATM + 1 1262 IF (IATM.LE.NATM) GOTO 150 1263 IF (fancy.eq.1 .or. fancy.eq.3) then 1264 write (output,'(A)') '9 end transparent ellipsoids' 1265 endif 1266 IF (fancy.eq.5 .or. fancy.eq.6) then 1267 write (output,'(A)') '9 end ortep ellipsoids' 1268 endif 1269 160 continue 1270c 1271c Write bonds to file also. Atoms are considered bonded if they lie 1272c closer to each other than 0.6 * sum of VDW radii. 1273C If two atoms of different colors are bonded, make half-bond 1274C cylinders with each color. 1275C 1276 if (nerrors.eq.0) write(noise,*) 1277 & '... no errors found in input file' 1278 if (radius.eq.0.0 .and. .not.suvflag .and. .not.cnflag) goto 210 1279 if (suvflag) then 1280 write (suvlun,'(A,A,F5.3,A)') 1281 & 'Checking for neighboring atoms with dissimilar Uij ', 1282 & '(Suv < ',suvlimit,')...' 1283 suvbad = 0 1284 endif 1285 if (cnflag) then 1286 write (NOISE,'(A,F5.3)') 1287 & 'Checking for C-N linkages with CCuij < ', cn_min_ccuij 1288 cnbad = 0 1289 endif 1290c 1291 DO 202 IATM=1,NATM 1292 IF (nohydro .AND. ATOM(IATM)(77:78).EQ.' H') GOTO 202 1293 if (suvflag .or. cnflag) then 1294 if (Biso(IATM).eq.0.0) goto 202 1295 if (uij(1,IATM).lt.0.) goto 202 1296 endif 1297 XI = SPAM(1,IATM) 1298 YI = SPAM(2,IATM) 1299 ZI = SPAM(3,IATM) 1300 ICOL = SPAM(5,IATM) 1301 VDWI = VDW(ICOL) 1302 DO 201 JATM=IATM+1,NATM 1303 CLOSE2 = 4.537 1304 DX2 = (XI - SPAM(1,JATM))**2 1305 if (dx2 .gt. close2) goto 201 1306 DY2 = (YI - SPAM(2,JATM))**2 1307 if (dy2 .gt. close2) goto 201 1308 DZ2 = (ZI - SPAM(3,JATM))**2 1309 DIST2 = DX2 + DY2 + DZ2 1310c 1311c Checking for bonded atoms with dissimilar Uij 1312 if (suvflag .or. cnflag) then 1313 IF (DIST2 .GT. CLOSE2) GOTO 201 1314 IF (Biso(JATM).eq.0.0) goto 180 1315 if (uij(1,jatm).lt.0.) goto 180 1316CDEBUG Version 2.6c used explicit test on 0.6*VdW distance 1317CDEBUG but this is very slow so I have fixed CLOSE2 at the C-S bond length 1318CDEBUG It might be worth doing a first cut using, say, 2.25A (>S-S bond) 1319CDEBUG and then only check the actual VdW distance for pairs making the cut 1320C JCOL = SPAM(5,JATM) 1321C VDWJ = VDW(JCOL) 1322C CLOSE2 = (0.6 * (VDWI + VDWJ)) **2 1323C IF (DIST2 .GT. CLOSE2) GOTO 201 1324CDEBUG Add this section back to see if results match V2.6c numbers 1325 similarity = Suv( uij(1,iatm), uij(1,jatm) ) 1326 if (suvflag .and. (similarity .lt. suvlimit)) then 1327 write (suvlun,161) 1328 & ATOM(IATM)(13:17),ATOM(IATM)(18:27), 1329 & ATOM(JATM)(13:17),ATOM(JATM)(18:27), 1330 & similarity 1331 suvbad = suvbad + 1 1332 161 format(1X,A5,1X,A10,8X,A5,1X,A10,4X,'Suv = ',F8.4) 1333 endif 1334c 1335c Check in particular for C-N links that may be TLS group boundaries 1336 180 continue 1337 if (cnflag) then 1338 ccuijpair = 0 1339 if ((ATOM(IATM)(13:15) .eq. ' C ') .and. (ATOM(JATM)(13:15) .eq. ' N ')) then 1340 ccuijpair = 1 1341 endif 1342 if ((ATOM(IATM)(13:15) .eq. ' O3') .and. (ATOM(JATM)(13:15) .eq. ' P ')) then 1343 ccuijpair = 2 1344 endif 1345 1346 if (ccuijpair.ne.0) then 1347 cc = CCuv( uij(1,iatm), uij(1,jatm) ) 1348 if (ATOM(IATM)(22:22).ne.prevchain) then 1349 prevchain = ATOM(IATM)(22:22) 1350 write (cnlun,'(1x)') 1351 write (cnlun,'(1x)') 1352 endif 1353 name_i = ATOM(IATM)(18:27) 1354 name_j = ATOM(JATM)(18:27) 1355 if (name_i(5:5) .eq. ' ') name_i(5:5) = 'X' 1356 if (name_j(5:5) .eq. ' ') name_j(5:5) = 'X' 1357 do ichar = 6,9 1358 if (name_i(ichar:ichar).eq.'_') name_i(ichar:ichar) = ' ' 1359 if (name_j(ichar:ichar).eq.'_') name_j(ichar:ichar) = ' ' 1360 enddo 1361 write (cnlun,182) name_i, name_j, cc 1362 if (cc .lt. cn_min_ccuij) then 1363 if (ccuijpair.eq.1) write (NOISE,183) name_i, name_j, cc 1364 if (ccuijpair.eq.2) write (NOISE,184) name_i, name_j, cc 1365 cnbad = cnbad + 1 1366 endif 1367 182 format(1X,A10,8X,A10,4X,'CCuij = ',F8.4) 1368 183 format(1X,A15,8X,A15,4X,'C-N CCuij = ',F8.4) 1369 184 format(1X,A15,8X,A15,4X,'O3''-P CCuij = ',F8.4) 1370 endif 1371 endif 1372 1373 if (tflag) goto 201 1374 endif 1375 1376c More stringent distance test when drawing bonds 1377 IF (nohydro .AND. ATOM(JATM)(77:78).EQ.' H') GOTO 201 1378 JCOL = SPAM(5,JATM) 1379 VDWJ = VDW(JCOL) 1380 CLOSE2 = (0.6 * (VDWI + VDWJ)) **2 1381 IF (DIST2 .GT. CLOSE2) GOTO 201 1382 1383c Don't draw bonds between alternate conformers of same residue 1384 IF (ATOM(IATM)(17:17).ne.' '.and.ATOM(JATM)(17:17).ne.' ' 1385 & .and. ATOM(IATM)(17:17).ne.ATOM(JATM)(17:17)) goto 201 1386c 1387c Atoms coloured by B value 1388 if (bcflag) then 1389 write(output,211) 1390 1 SPAM(1,IATM),SPAM(2,IATM),SPAM(3,IATM),radius, 1391 2 SPAM(1,JATM),SPAM(2,JATM),SPAM(3,JATM),radius, 1392 3 1.0, 1.0, 1.0 1393 call U2RGB( SPAM(4,IATM), Umin, Umax, RED1, GREEN1, BLUE1 ) 1394 call U2RGB( SPAM(4,JATM), Umin, Umax, RED2, GREEN2, BLUE2 ) 1395 write(output,212) 1396 1 RED1,GREEN1,BLUE1, RED2,GREEN2,BLUE2, 0., 0., 0. 1397c Same color atoms 1398 elseif (RGB(1,ICOL) .EQ. RGB(1,JCOL) .AND. 1399 1 RGB(2,ICOL) .EQ. RGB(2,JCOL) .AND. 1400 2 RGB(3,ICOL) .EQ. RGB(3,JCOL)) THEN 1401 WRITE(OUTPUT,211) 1402 1 SPAM(1,IATM),SPAM(2,IATM),SPAM(3,IATM),radius, 1403 2 SPAM(1,JATM),SPAM(2,JATM),SPAM(3,JATM),radius, 1404 3 RGB(1,ICOL),RGB(2,ICOL),RGB(3,ICOL) 1405 ELSE 1406 DO K=1,3 1407 center(K) = (SPAM(K,IATM)+SPAM(K,JATM))/2 1408 ENDDO 1409 WRITE(OUTPUT,211) 1410 1 SPAM(1,IATM),SPAM(2,IATM),SPAM(3,IATM),radius, 1411 2 center(1),center(2),center(3),radius, 1412 3 RGB(1,ICOL),RGB(2,ICOL),RGB(3,ICOL) 1413 WRITE(OUTPUT,211) 1414 1 center(1),center(2),center(3),radius, 1415 2 SPAM(1,JATM),SPAM(2,JATM),SPAM(3,JATM),radius, 1416 3 RGB(1,JCOL),RGB(2,JCOL),RGB(3,JCOL) 1417 ENDIF 1418 201 CONTINUE 1419 202 CONTINUE 1420 210 CONTINUE 1421 1422C211 FORMAT(1H3,/,11(f8.3)) 1423211 FORMAT(1H3,/,3(1x,f8.3),f7.3,3(1x,f8.3),f7.3,3(f6.3)) 1424212 FORMAT(2H17,/,9f8.3) 1425c 1426 if (suvflag) then 1427 if (suvbad.eq.0) write (suvlun,*) '... No Suv outliers' 1428 endif 1429 if (cnflag) then 1430 if (cnbad.eq.0) write (suvlun,*) '... No CCuij outliers' 1431 endif 1432 if (suvflag .or. cnflag) call exit(0) 1433c 1434 write (noise,'(/)') 1435 write (noise,156) 'X min max center-of-mass:', XMIN, XMAX, xcom 1436 write (noise,156) 'Y min max center-of-mass:', YMIN, YMAX, ycom 1437 write (noise,156) 'Z min max center-of-mass:', ZMIN, ZMAX, zcom 1438 write (noise,156) ' scale:', SCALE 1439 156 format(1x,a,2f8.2,f10.3) 1440 END 1441 1442 1443 LOGICAL FUNCTION MATCH (SUBJ, MASK) 1444 CHARACTER*24 SUBJ,MASK 1445 MATCH = .FALSE. 1446 DO 10 I = 1, 24 1447 IF (SUBJ(I:I).NE.MASK(I:I) .AND. MASK(I:I).NE.'#') RETURN 144810 CONTINUE 1449 MATCH = .TRUE. 1450 RETURN 1451 END 1452 1453C Dummy routine to make linker happy (called by QINP in quadric.f) 1454 SUBROUTINE TRANSF (X,Y,Z) 1455 RETURN 1456 END 1457 1458 SUBROUTINE ASSERT (LOGIC, DAMMIT) 1459 LOGICAL LOGIC 1460 CHARACTER*(*) DAMMIT 1461 COMMON /ASSOUT/ NOISE 1462 IF (LOGIC) RETURN 1463 WRITE (NOISE,*) 'ERROR >>>>>> ',DAMMIT 1464C STOP 1234 1465 CALL EXIT(-1) 1466 END 1467 1468 1469CCC Return RGB triple that runs from dark blue at Bmin 1470CC to light red at Bmax 1471C 1472 subroutine U2rgb( Uiso, Umin, Umax, r, g, b ) 1473 real Uiso, Umin, Umax, r, g, b 1474c 1475 real fraction, h, s, v 1476c 1477 fraction = (Uiso-Umin) / (Umax-Umin) 1478 if (fraction.lt.0.) fraction = 0. 1479 if (fraction.gt.1.) fraction = 1. 1480 h = 240. * (1.-fraction) 1481 s = 0.95 1482 v = 0.75 + fraction/4. 1483 call hsv2rgb( h, s, v, r, g, b ) 1484 return 1485 end 1486 1487 1488CCC Return RGB triple that runs from 1489CC red for A=0.0 -> white for A=0.5 -> green for A=1.0 1490C 1491 subroutine A2rgb( A, r, g, b ) 1492 real A, r, g, b 1493c 1494 real fraction, h, s, v 1495c 1496 if (A .lt. 0.5) h = 0.0 1497 if (A .ge. 0.5) h = 120. 1498 fraction = abs( (A-0.5)*2.0 ) 1499c s = fraction**2 1500 s = fraction 1501 v = 1.0 - s/2.0 1502 if (A .le. 0.0) then 1503 h = 300. 1504 v = 1.0 1505 endif 1506 call hsv2rgb( h, s, v, r, g, b ) 1507 return 1508 end 1509 1510 1511CCC Color format conversion from Hue/Saturation/Value to Red/Green/Blue 1512CC minimal (i.e. NO) error checking 1513C 1514 subroutine hsv2rgb( h, s, v, r, g, b ) 1515 real h, s, v, r, g, b 1516c 1517 real f, p, q, t 1518 integer i 1519c 1520 i = h /60. 1521 f = h/60. - float(i) 1522 p = v * (1. - s) 1523 q = v * (1. - s*f) 1524 t = v * (1. - s*(1. - f)) 1525 if (i.eq.5) then 1526 r = v 1527 g = p 1528 b = q 1529 else if (i.eq.4) then 1530 r = t 1531 g = p 1532 b = v 1533 else if (i.eq.3) then 1534 r = p 1535 g = q 1536 b = v 1537 else if (i.eq.2) then 1538 r = p 1539 g = v 1540 b = t 1541 else if (i.eq.1) then 1542 r = q 1543 g = v 1544 b = p 1545 else 1546 r = v 1547 g = t 1548 b = p 1549 endif 1550 return 1551 end 1552 1553CCC This subroutine is not strictly necessary. 1554CC It smooths the histogram/curve of Anisotropy by 1555C distance from center of mass 1556 1557 subroutine hisclean( value, number, shells, nanis ) 1558c 1559 integer shells 1560 real value(shells) 1561 integer number(shells) 1562 integer nanis 1563c 1564 integer minshell, maxshell, min10, max10 1565 integer nsum 1566 real vsum 1567c 1568 if (nanis .lt. 1200) then 1569 shells = shells / 2 1570 do i = 1, shells 1571 number(i) = number(2*i-1) + number(2*i) 1572 value(i) = value(2*i-1) + value(2*i) 1573 enddo 1574 endif 1575c 1576 do i = shells, 1, -1 1577 if (number(i).gt.0) minshell = i 1578 enddo 1579 do i = 1, shells 1580 if (number(i).gt.0) maxshell = i 1581 enddo 1582 nsum = 0 1583 min10 = minshell 1584 do i = minshell, shells 1585 if (nsum + number(i) .gt. 10) then 1586 goto 101 1587 else 1588 nsum = nsum + number(i) 1589 min10 = i 1590 endif 1591 enddo 1592101 continue 1593 nsum = 0 1594 max10 = maxshell 1595 do i = maxshell, 1, -1 1596 if (nsum + number(i) .gt. 10) then 1597 goto 102 1598 else 1599 nsum = nsum + number(i) 1600 max10 = i 1601 endif 1602 enddo 1603102 continue 1604c 1605 nsum = 0 1606 vsum = 0.0 1607 do i = minshell, min10 1608 nsum = nsum + number(i) 1609 vsum = vsum + value(i) 1610 number(i) = 0 1611 value(i) = 0.0 1612 enddo 1613 i = (minshell + min10) / 2 1614 number(i) = nsum 1615 value(i) = vsum 1616c 1617 nsum = 0 1618 vsum = 0.0 1619 do i = maxshell, max10, -1 1620 nsum = nsum + number(i) 1621 vsum = vsum + value(i) 1622 number(i) = 0 1623 value(i) = 0.0 1624 enddo 1625 i = (maxshell + max10) / 2 1626 number(i) = nsum 1627 value(i) = vsum 1628c 1629 return 1630 end 1631 1632