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