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