1 2 PROGRAM BALLS 3* 4* Program to set up input for RENDER 5* 6* This program takes a Brookhaven PDB-format file of atomic coordinates 7* concatenated with a "colours" file that specifies the colouring and spherical 8* radius of atoms according to atom type, residue type, or any other criterion, 9* and produces a file that can be used directly as input to RENDER. 10* 11* Usage: cat colors.pdb protein.pdb | balls [-h] | render 12* 13* auxiliary files: setup.matrix optional file containing 3x3 view matrix 14* setup.angles optional file defining view matrix 15* in terms of 3 rotation angles 16* 17* I/O units for colour/co-ordinate input, specs output, user output 18* 19 INCLUDE 'VERSION.incl' 20* 21 INTEGER INPUT, OUTPUT, NOISE 22 PARAMETER (INPUT=5, OUTPUT=6, NOISE=0) 23 PARAMETER (MAXCOL=5000, MAXATM=300000) 24 REAL RGB(3,MAXCOL),RADIUS(MAXCOL) 25 REAL SPAM(7,MAXATM) 26 CHARACTER*24 MASK(MAXCOL),TEST 27 CHARACTER*80 ATOM(MAXATM),CARD 28 LOGICAL MATCH 29 logical hflag 30 character*80 flags 31c 32C Ethan Merritt Oct 1988 33C Modified to read in 3x3 view matrix (e.g. from CCP FRODO view command) 34C from file setup.matrix. Matrix is applied before 35C finding translation, center, and scale. Afterwards the input matrix 36C to RENDER is therefore the identity matrix. 37C 38 common /matrix/ matrix, postrn, coords 39 real*4 matrix(3,3), postrn(3), coords(3) 40C 41C Default to CPK colors and VDW radii 42 character*60 defcol(7) 43 data defcol / 44 & 'COLOUR#######C################ 0.625 0.625 0.625 1.70', 45 & 'COLOUR#######N################ 0.125 0.125 1.000 1.60', 46 & 'COLOUR#######O################ 0.750 0.050 0.050 1.50', 47 & 'COLOUR#######S################ 1.000 1.000 0.025 1.85', 48 & 'COLOUR#######H################ 1.000 1.000 1.000 1.20', 49 & 'COLOUR#######P################ 0.400 1.000 0.400 1.80', 50 & 'COLOUR######################## 1.000 0.000 1.000 2.00' 51 & / 52c 53 3 format(a,a) 54c 55 do i=1,3 56 do j=1,3 57 matrix(i,j)=0. 58 enddo 59 matrix(i,i)=1. 60 enddo 61 call view_matrix 62c 63c Ethan Merritt Apr 1992 64c -h option suppresses header records in output 65c 66 hflag = .false. 67 narg = iargc() 68 do i = 1, narg 69 call getarg( i, flags ) 70 if (flags(1:2) .eq. '-h') hflag = .true. 71 end do 72c 73 if (.not. hflag) then 74 WRITE(OUTPUT,'(A,A)') 'balls - Raster3D ', VERSION 75 WRITE(OUTPUT,'(A)') '80 64 tiles in x,y' 76 WRITE(OUTPUT,'(A)') ' 8 8 pixels (x,y) per tile' 77 WRITE(OUTPUT,'(A)') '4 anti-aliasing 3x3 -> 2x2 pixels' 78 WRITE(OUTPUT,'(A)') '0 0 0 black background' 79 WRITE(OUTPUT,'(A)') 'T yes, I LIKE shadows!' 80 WRITE(OUTPUT,'(A)') '25 Phong power' 81 WRITE(OUTPUT,'(A)') '0.15 secondary light contribution' 82 WRITE(OUTPUT,'(A)') '0.05 ambient light contribution' 83 WRITE(OUTPUT,'(A)') '0.25 specular reflection component' 84 WRITE(OUTPUT,'(A)') '4.0 eye position' 85 WRITE(OUTPUT,'(A)') '1 1 1 main light source position' 86 end if 87c 88 ASPECT = 1280./1024. 89 NCOL = 0 90 NATM = 0 9110 CONTINUE 92 READ(INPUT,'(A80)',END=50) CARD 93 IF (CARD(1:4).EQ.'COLO') THEN 94 NCOL = NCOL + 1 95 IF (NCOL.GT.MAXCOL) THEN 96 WRITE(NOISE,*) 'Colour table overflow. Increase ', 97 & 'MAXCOL and recompile.' 98 STOP 10 99 ENDIF 100 READ(CARD,'(6X,A24,3F8.3,F6.2)') MASK(NCOL), 101 & (RGB(I,NCOL),I=1,3),RADIUS(NCOL) 102 ELSEIF (CARD(1:4).EQ.'ATOM'.OR.CARD(1:4).EQ.'HETA') THEN 103 NATM = NATM + 1 104 IF (NATM.GT.MAXATM) THEN 105 WRITE(NOISE,*) 'Atom array overflow. Increase ', 106 & 'MAXATM and recompile.' 107 STOP 20 108 ENDIF 109 ATOM(NATM) = CARD 110 ELSEIF (CARD(1:3).EQ.'END') THEN 111 GO TO 50 112 ENDIF 113 GO TO 10 114* Come here when EOF or 'END' record is reached 11550 CONTINUE 116 IF (NATM.EQ.0) THEN 117 WRITE(NOISE,*) 'No atoms in input.' 118 STOP 30 119 ENDIF 120* Load default colors after any that were read in 121 IF (NCOL.LT.MAXCOL-8) THEN 122 DO i = 1,7 123 NCOL = NCOL + 1 124 READ(defcol(i),'(6X,A24,3F8.3,F6.2)') MASK(NCOL), 125 & (RGB(J,NCOL),J=1,3), RADIUS(NCOL) 126 ENDDO 127 ENDIF 128* 129 IF (NCOL.EQ.0) THEN 130 WRITE(NOISE,*) 'No colours in input.' 131 STOP 40 132 ENDIF 133 XMAX = -1E20 134 XMIN = 1E20 135 YMAX = -1E20 136 YMIN = 1E20 137 ZMAX = -1E20 138 ZMIN = 1E20 139 DO 100 IATM=1,NATM 140 CARD = ATOM(IATM) 141 TEST = CARD(7:30) 142 DO 80 ICOL=1,NCOL 143 IF (MATCH(TEST,MASK(ICOL))) THEN 144c READ(CARD,'(30X,3F8.3)') X,Y,Z 145c EAM Oct88 146 READ(CARD,'(30X,3F8.3)') coords 147 x = coords(1)*matrix(1,1) + coords(2)*matrix(2,1) 148 1 + coords(3)*matrix(3,1) 149 y = coords(1)*matrix(1,2) + coords(2)*matrix(2,2) 150 1 + coords(3)*matrix(3,2) 151 z = coords(1)*matrix(1,3) + coords(2)*matrix(2,3) 152 1 + coords(3)*matrix(3,3) 153c EAM Oct88 154 RAD = RADIUS(ICOL) 155 SPAM(1,IATM) = X 156 SPAM(2,IATM) = Y 157 SPAM(3,IATM) = Z 158 SPAM(4,IATM) = RAD 159 SPAM(5,IATM) = RGB(1,ICOL) 160 SPAM(6,IATM) = RGB(2,ICOL) 161 SPAM(7,IATM) = RGB(3,ICOL) 162 XMAX = MAX(XMAX,X+RAD) 163 XMIN = MIN(XMIN,X-RAD) 164 YMAX = MAX(YMAX,Y+RAD) 165 YMIN = MIN(YMIN,Y-RAD) 166 ZMAX = MAX(ZMAX,Z+RAD) 167 ZMIN = MIN(ZMIN,Z-RAD) 168 GO TO 100 169 ENDIF 17080 CONTINUE 171 WRITE(NOISE,*) 'No colour table mask matches this atom:' 172 WRITE(NOISE,*) ATOM(IATM) 173 STOP 90 174100 CONTINUE 175 XMID = (XMAX+XMIN)/2. 176 YMID = (YMAX+YMIN)/2. 177 ZMID = (ZMAX+ZMIN)/2. 178 TX = -XMID 179 TY = -YMID 180 TZ = -ZMID 181 IF (ASPECT.GE.1.) THEN 182* The X direction is wider than the Y 183 XROOM = ASPECT 184 YROOM = 1. 185 ZROOM = 2. 186 ELSE 187 XROOM = 1. 188 YROOM = ASPECT 189 ZROOM = 2. 190 ENDIF 191 XSPAN = XMAX-XMIN 192 YSPAN = YMAX-YMIN 193 ZSPAN = ZMAX-ZMIN 194 SCALE = MAX(XSPAN/XROOM,YSPAN/YROOM,ZSPAN/ZROOM) 195* Leave a little extra room as a border: 196 SCALE = SCALE / 0.90 197c 198 if (.not. hflag) then 199 WRITE(OUTPUT,120) TX,TY,TZ,SCALE 200120 FORMAT('1 0 0 0 input co-ordinate + radius transformation'/ 201 & '0 1 0 0'/ 202 & '0 0 1 0'/ 203 & 4F10.3) 204 WRITE(OUTPUT,'(A)') '3 mixed object types' 205 WRITE(OUTPUT,'(A)') '*' 206 WRITE(OUTPUT,'(A)') '*' 207 WRITE(OUTPUT,'(A)') '*' 208 end if 209c 210 DO 150 IATM=1,NATM 211 SPAM(1,IATM) = SPAM(1,IATM) + postrn(1) 212 SPAM(2,IATM) = SPAM(2,IATM) + postrn(2) 213 SPAM(3,IATM) = SPAM(3,IATM) + postrn(3) 214 WRITE(OUTPUT,140) (SPAM(I,IATM),I=1,7) 215140 FORMAT(1H2,/,7(1X,F8.3)) 216150 CONTINUE 217 write (noise,'(/)') 218 write (noise,156) 'X min max:', XMIN, XMAX 219 write (noise,156) 'Y min max:', YMIN, YMAX 220 write (noise,156) 'Z min max:', ZMIN, ZMAX 221 write (noise,156) ' scale:', SCALE 222 156 format(1x,a,3f8.2) 223 END 224 LOGICAL FUNCTION MATCH (SUBJ, MASK) 225 CHARACTER*24 SUBJ,MASK 226 MATCH = .FALSE. 227 DO 10 I = 1, 24 228 IF (SUBJ(I:I).NE.MASK(I:I) .AND. MASK(I:I).NE.'#') RETURN 22910 CONTINUE 230 MATCH = .TRUE. 231 RETURN 232 END 233 234 subroutine view_matrix 235c 236 common /matrix/ matrix, postrn, coords 237 real*4 matrix(3,3), postrn(3), coords(3) 238c 239 real*4 phiX, phiY, phiZ 240 integer noise 241 parameter (noise = 0) 242 parameter (R2D = 180./3.1415927) 243 244 postrn(1) = 0. 245 postrn(2) = 0. 246 postrn(3) = 0. 247 248 open (unit=3, file='setup.matrix', status='OLD', err=100) 249 read (3,*) ((matrix(i,j),i=1,3),j=1,3) 250 write (noise,'(1x,3f9.5)') ((matrix(i,j),i=1,3),j=1,3) 251 close (3) 252 253 det = matrix(1,1) * matrix(2,2) * matrix(3,3) 254 1 + matrix(1,2) * matrix(2,3) * matrix(3,1) 255 2 + matrix(2,1) * matrix(3,2) * matrix(1,3) 256 3 - matrix(1,3) * matrix(2,2) * matrix(3,1) 257 4 - matrix(1,2) * matrix(2,1) * matrix(3,3) 258 5 - matrix(1,1) * matrix(2,3) * matrix(3,2) 259 write (noise,'('' determinant ='',f8.3)') det 260 261 phiX = atan2( -matrix(3,2), matrix(3,3) ) 262 phiY = atan2( matrix(3,1), matrix(3,3) / cos(phiX) ) 263 phiZ = atan2( -matrix(2,1), matrix(1,1) ) 264 write (noise,3) ' View Angles from matrix',' ' 265 write (noise,2) phiZ*R2D, phiY*R2D, phiX*R2D 266 return 267 100 continue 268 269 open (unit=3, file='setup.angles', status='OLD', err=200) 270 read (3,*) phiZ, phiY, phiX 271 read (3,*,err=101,end=101) postrn(1), postrn(2), postrn(3) 272 101 continue 273 close (3) 274 write (noise,2) phiZ, phiY, phiX 275 write (noise,4) postrn(1), postrn(2), postrn(3) 276 cx = cos(phiX/R2D) 277 sx = sin(phiX/R2D) 278 cy = cos(phiY/R2D) 279 sy = sin(phiY/R2D) 280 cz = cos(phiZ/R2D) 281 sz = sin(phiZ/R2D) 282 matrix(1,1) = cz*cy 283 matrix(1,2) = sz*cx + cz*sy*sx 284 matrix(1,3) = sz*sx - cz*sy*cx 285 matrix(2,1) = -sz*cy 286 matrix(2,2) = cz*cx - sx*sy*sz 287 matrix(2,3) = cz*sx + sz*sy*cx 288 matrix(3,1) = sy 289 matrix(3,2) = -sx*cy 290 matrix(3,3) = cx*cy 291 write (noise,3) ' View Matrix from angles',' ' 292 write (noise,'(1x,3f9.5)') ((matrix(i,j),i=1,3),j=1,3) 293 return 294 200 continue 295 296 2 format(1x,' phiZ =',f8.2,' phiY =',f8.2,' phiX =',f8.2) 297 3 format(/a,a) 298 4 format(1x,' tranX =',f8.2,' tranY =',f8.2,' tranZ =',f8.2) 299 300 write (noise,*) ' No view matrix or angles provided' 301 return 302 end 303