/* rastep.f -- translated by f2c (version 20061008). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib; on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */ #include "f2c.h" /* Common Block Declarations */ struct { integer assout; logical verbose; } asscom_; #define asscom_1 asscom_ struct { integer noise; } assout_; #define assout_1 assout_ /* Table of constant values */ static integer c__4 = 4; static integer c__1 = 1; static integer c__9 = 9; static integer c_n1 = -1; static integer c__5 = 5; static integer c__3 = 3; static real c_b332 = 1.f; static integer c__6 = 6; static integer c__10 = 10; static integer c__14 = 14; static integer c__0 = 0; static integer c__2 = 2; static integer c__8 = 8; static real c_b566 = .5f; static real c_b670 = 0.f; /* Main program */ int MAIN__(void) { /* Initialized data */ static char defcol[60*17] = "COLOUR###### CA ############## 0.175 0." "175 0.175 1.70" "COLOUR###### C ############## 0.175 0.1" "75 0.175 1.70" "COLOUR#######C################ 0.625 0.62" "5 0.625 1.70" "COLOUR#######N################ 0.125 0.125" " 1.000 1.60" "COLOUR#######O################ 0.750 0.050 " " 0.050 1.50" "COLOUR#######S################ 1.000 1.000 " " 0.025 1.85" "COLOUR#######H################ 1.000 1.000 " "1.000 1.20" "COLOUR#######P################ 0.050 0.750 0" ".050 1.80" "COLOUR##### CA ############### 0.175 0.175 0." "175 1.70" "COLOUR##### C ############### 0.175 0.175 0.1" "75 1.70" "COLOUR######C################# 0.625 0.625 0.62" "5 1.70" "COLOUR######N################# 0.125 0.125 1.000" " 1.60" "COLOUR######O################# 0.750 0.050 0.050 " " 1.50" "COLOUR######S################# 1.000 1.000 0.025 " "1.85" "COLOUR######H################# 1.000 1.000 1.000 1" ".20" "COLOUR######P################# 0.050 0.750 0.050 1." "80" "COLOUR######################## 1.000 0.000 1.000 2.00" ; static real problevel[50] = { .4299f,.5479f,.6334f,.7035f,.7644f,.8192f, .8694f,.9162f,.9605f,1.0026f,1.043f,1.0821f,1.12f,1.157f,1.1932f, 1.2288f,1.2638f,1.2985f,1.333f,1.3672f,1.4013f,1.4354f,1.4695f, 1.5037f,1.5382f,1.5729f,1.608f,1.6436f,1.6797f,1.7164f,1.754f, 1.7924f,1.8318f,1.8724f,1.9144f,1.958f,2.0034f,2.051f,2.1012f, 2.1544f,2.2114f,2.273f,2.3404f,2.4153f,2.5003f,2.5997f,2.7216f, 2.8829f,3.1365f,6.f }; /* Format strings */ static char fmt_800[] = "(\002******************************************" "******\002)"; static char fmt_801[] = "(\002 Spheres will bound Biso probability leve" "l\002,f5.2)"; static char fmt_802[] = "(\002 Ellipsoids will bound probability leve" "l\002,f5.2)"; static char fmt_803[] = "(\002 Corresponding critical value " " \002,f7.4)"; static char fmt_120[] = "(\0021 0 0 0 input co-ordinate + radius trans" "formation\002/\0020 1 0 0\002/\0020 0 1 0\002/4f10.3)"; static char fmt_901[] = "(a,3f8.3,6f8.4)"; static char fmt_902[] = "(a,10f8.3)"; static char fmt_903[] = "(a,3f8.3,4x,a,f8.3,4x,a,f8.3)"; static char fmt_904[] = "(a,9f7.3)"; static char fmt_905[] = "(a5,1x,a5,a1,a5,3f9.4,2x,f9.4,f9.4,f9.4)"; static char fmt_907[] = "(\002 3\002,/,3f9.3,\002 0.02\002,3f9.3,\002 0." "02\002,\002 0.5 1.0 0.3\002)"; static char fmt_151[] = "(i2,/,3(1x,f8.3),4f8.3)"; static char fmt_152[] = "(10f12.4)"; static char fmt_146[] = "(f10.3,f10.3,i10)"; static char fmt_171[] = "(\002BOUNDING_COLOR \002,3f6.3)"; static char fmt_172[] = "(\002BOUNDING_PLANE \002,i2,6f10.4)"; static char fmt_161[] = "(1x,a5,1x,a10,8x,a5,1x,a10,4x,\002Suv = \002,f8" ".4)"; static char fmt_211[] = "(\0023\002,/,3(1x,f8.3),f7.3,3(1x,f8.3),f7.3,3(" "f6.3))"; static char fmt_212[] = "(\00217\002,/,9f8.3)"; static char fmt_156[] = "(1x,a,2f8.2,f10.3)"; /* System generated locals */ integer i__1, i__2, i__3; real r__1, r__2, r__3; doublereal d__1, d__2, d__3, d__4, d__5, d__6, d__7, d__8, d__9, d__10, d__11, d__12, d__13, d__14, d__15; icilist ici__1; olist o__1; /* Builtin functions */ integer s_cmp(char *, char *, ftnlen, ftnlen), s_rsli(icilist *), do_lio( integer *, integer *, char *, ftnlen), e_rsli(void); /* Subroutine */ int s_stop(char *, ftnlen); integer f_open(olist *), s_wsle(cilist *), e_wsle(void), s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void), s_rsfe(cilist *) , e_rsfe(void), s_rsfi(icilist *), e_rsfi(void); /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen); double sqrt(doublereal); /* Local variables */ extern /* Subroutine */ int hisclean_(real *, integer *, integer *, integer *); static logical ellipses; static char atomtype[2]; static real suvlimit; static integer i__, j, k, l; static real u[16] /* was [4][4] */, x, y, z__, anis_mean__, biso_mean__; extern integer anitoquad_(real *, real *, real *, real *, real *); static integer histogram[20]; static real rg, qp[16] /* was [4][4] */, qq[16] /* was [4][4] */, rq, rr[9] /* was [3][3] */, xi, yi, zi, rv, tx, ty, tz, anis_sigma__, dx2, dy2, dz2, similarity, anisotropy, end[3], rad, rgb[15000] /* was [3][5000] */, red, det, uij[1800000] /* was [6][300000] */, vdw[5000]; extern doublereal suv_(real *, real *); static real red1, red2, ellipticity; static char card[80]; static real dmid; static integer imed, kmid; static real bmin; static integer nani, narg; static real bmax; static integer icol, iatm; static char mask[24*5000]; static logical mini; static real biso[300000]; static integer ncol, natm; static real spam[1500000] /* was [5][300000] */; static char atom[80*300000]; static integer nhyd; static real prob, xmid, temp[16] /* was [4][4] */; static doublereal xcom, ycom; static logical auto__; static integer niso; static doublereal zcom; extern /* Subroutine */ int exit_(integer *); static real umin, umax, xmax; static char test[24]; static real xmom[5], xmin, ymax, ymin, zmax, zmin, uiso, ymid, zmid, dmax__; static doublereal wsum, xsum, ysum, zsum; static integer kmax, kmin; static real blue; extern /* Subroutine */ int a2rgb_(real *, real *, real *, real *); static real dist, size; static integer imax, imin; static real vdwi; static integer jatm, jcol; static real vdwj, blue1, blue2; extern /* Subroutine */ int u2rgb_(real *, real *, real *, real *, real *, real *); static real dist2; extern doublereal tinv4_(real *, real *); extern /* Subroutine */ int tmul4_(real *, real *, real *); static real ccoef; static logical hflag; extern integer iargc_(void); static real scale; extern logical match_(char *, char *, ftnlen, ftnlen); static logical tflag; static char flags[80]; static integer fancy; static real green, anisi[300000], sum_a__, evecs[16] /* was [4][4] */, umean, sum_b__; static integer nanis; static real adist[110]; static integer hdist[110], iprob; static real total; static integer noerr; static real xspan, yspan, zspan, start[3], uprin[3], xroom, yroom, zroom, green1, green2, sum_a2__, sum_b2__, close2; static logical acflag, bcflag; extern /* Subroutine */ int jacobi_(real *, integer *, integer *, real *, real *), trnsp4_(real *, real *); static logical atflag; static real sum_ab__, radlim; extern /* Subroutine */ int getarg_(integer *, char *, ftnlen); static real eigens[4]; static char spacer[1]; static real evecit[16] /* was [4][4] */, center[3], aspect; static integer suvbad; static real usigma, radius, weight; static integer comlun; static real anisou[6]; static integer natype, hislun, nonpos, suvlun; static logical comflag; static real alonghi, quadric[10]; static char resname[3]; static logical hwhacky; static real evecinv[16] /* was [4][4] */; static integer dshells; static logical suvflag; static real pradius; static logical nohydro; static integer nerrors; /* Fortran I/O blocks */ static icilist io___31 = { 1, flags, 0, 0, 80, 1 }; static icilist io___32 = { 1, flags, 0, 0, 80, 1 }; static icilist io___33 = { 1, flags, 0, 0, 80, 1 }; static icilist io___35 = { 1, flags, 0, 0, 80, 1 }; static icilist io___43 = { 1, flags, 0, 0, 80, 1 }; static cilist io___44 = { 0, 0, 0, 0, 0 }; static cilist io___45 = { 0, 0, 0, "(/,A)", 0 }; static cilist io___46 = { 0, 0, 0, "(A)", 0 }; static cilist io___47 = { 0, 0, 0, "(A)", 0 }; static cilist io___48 = { 0, 0, 0, "(A,A)", 0 }; static cilist io___49 = { 0, 0, 0, "(A,A)", 0 }; static cilist io___52 = { 0, 0, 0, fmt_800, 0 }; static cilist io___53 = { 0, 0, 0, 0, 0 }; static cilist io___54 = { 0, 0, 0, 0, 0 }; static cilist io___55 = { 0, 0, 0, fmt_800, 0 }; static cilist io___56 = { 0, 0, 0, fmt_801, 0 }; static cilist io___57 = { 0, 0, 0, fmt_802, 0 }; static cilist io___58 = { 0, 0, 0, fmt_803, 0 }; static cilist io___59 = { 0, 0, 0, 0, 0 }; static cilist io___60 = { 0, 0, 0, 0, 0 }; static cilist io___61 = { 0, 0, 0, 0, 0 }; static cilist io___62 = { 0, 0, 0, 0, 0 }; static cilist io___65 = { 0, 6, 0, "(A,A,I5,A)", 0 }; static cilist io___66 = { 0, 6, 0, "(A)", 0 }; static cilist io___67 = { 0, 6, 0, "(A)", 0 }; static cilist io___68 = { 0, 6, 0, "(A)", 0 }; static cilist io___69 = { 0, 6, 0, "(A)", 0 }; static cilist io___70 = { 0, 6, 0, "(A)", 0 }; static cilist io___71 = { 0, 6, 0, "(A)", 0 }; static cilist io___72 = { 0, 6, 0, "(A)", 0 }; static cilist io___73 = { 0, 6, 0, "(A)", 0 }; static cilist io___74 = { 0, 6, 0, "(A)", 0 }; static cilist io___75 = { 0, 6, 0, "(A)", 0 }; static cilist io___76 = { 0, 6, 0, "(A)", 0 }; static cilist io___77 = { 0, 6, 0, "(A)", 0 }; static cilist io___86 = { 0, 5, 1, "(A80)", 0 }; static cilist io___88 = { 0, 0, 0, 0, 0 }; static icilist io___89 = { 0, card, 0, "(6X,A24,3F8.3,F6.2)", 80, 1 }; static icilist io___94 = { 1, card+28, 1, 0, 42, 1 }; static cilist io___97 = { 0, 0, 0, 0, 0 }; static cilist io___98 = { 0, 0, 0, 0, 0 }; static cilist io___99 = { 0, 0, 0, 0, 0 }; static cilist io___100 = { 0, 0, 0, 0, 0 }; static cilist io___101 = { 0, 0, 0, 0, 0 }; static cilist io___102 = { 0, 0, 0, 0, 0 }; static icilist io___113 = { 1, card, 1, "(30X,3F8.3,6X,F6.2)", 80, 1 }; static cilist io___118 = { 0, 0, 0, 0, 0 }; static cilist io___123 = { 0, 0, 0, 0, 0 }; static cilist io___124 = { 0, 0, 0, 0, 0 }; static cilist io___125 = { 0, 0, 0, 0, 0 }; static cilist io___126 = { 0, 0, 0, 0, 0 }; static cilist io___145 = { 0, 6, 0, fmt_120, 0 }; static cilist io___146 = { 0, 6, 0, "(A)", 0 }; static cilist io___147 = { 0, 6, 0, "(A)", 0 }; static cilist io___148 = { 0, 6, 0, "(A)", 0 }; static cilist io___149 = { 0, 6, 0, "(A)", 0 }; static cilist io___154 = { 0, 0, 0, "(/,A,5G13.5)", 0 }; static cilist io___155 = { 0, 0, 0, "(A)", 0 }; static cilist io___160 = { 0, 0, 0, "(3G13.6)", 0 }; static cilist io___166 = { 0, 0, 0, "(A,/,3F13.5,10X,3i2)", 0 }; static cilist io___167 = { 0, 0, 0, "(A)", 0 }; static cilist io___168 = { 0, 0, 0, "(3F13.5)", 0 }; static cilist io___171 = { 0, 6, 0, "(A)", 0 }; static cilist io___172 = { 0, 6, 0, "(A)", 0 }; static cilist io___173 = { 0, 6, 0, "(A)", 0 }; static cilist io___174 = { 0, 6, 0, "(3F13.5)", 0 }; static cilist io___175 = { 0, 6, 0, "(A)", 0 }; static cilist io___176 = { 0, 6, 0, "(A,A)", 0 }; static cilist io___177 = { 0, 6, 0, "(A,F5.2)", 0 }; static cilist io___183 = { 0, 0, 0, 0, 0 }; static cilist io___185 = { 0, 0, 0, fmt_901, 0 }; static cilist io___186 = { 0, 0, 0, fmt_902, 0 }; static cilist io___187 = { 0, 0, 0, fmt_903, 0 }; static cilist io___188 = { 0, 0, 0, fmt_904, 0 }; static cilist io___196 = { 0, 6, 0, fmt_905, 0 }; static cilist io___201 = { 0, 6, 0, fmt_907, 0 }; static cilist io___205 = { 0, 6, 0, fmt_907, 0 }; static cilist io___210 = { 0, 6, 0, fmt_151, 0 }; static cilist io___211 = { 0, 6, 0, fmt_152, 0 }; static cilist io___212 = { 0, 6, 0, "(A,/,A)", 0 }; static cilist io___213 = { 0, 6, 0, "(A)", 0 }; static cilist io___215 = { 0, 0, 0, 0, 0 }; static cilist io___216 = { 0, 0, 0, 0, 0 }; static cilist io___217 = { 0, 0, 0, 0, 0 }; static cilist io___218 = { 0, 0, 0, "(A)", 0 }; static cilist io___219 = { 0, 0, 0, "(A)", 0 }; static cilist io___220 = { 0, 0, 0, "(2F5.2,3X,F8.3,I10)", 0 }; static cilist io___228 = { 0, 0, 0, "(A)", 0 }; static cilist io___229 = { 0, 0, 0, "(A,I10)", 0 }; static cilist io___230 = { 0, 0, 0, "(A,I10)", 0 }; static cilist io___231 = { 0, 0, 0, "(A,I10)", 0 }; static cilist io___232 = { 0, 0, 0, "(A,I10)", 0 }; static cilist io___233 = { 0, 0, 0, "(A,I10)", 0 }; static cilist io___234 = { 0, 0, 0, "(A)", 0 }; static cilist io___235 = { 0, 0, 0, "(A,F7.3)", 0 }; static cilist io___236 = { 0, 0, 0, "(A)", 0 }; static cilist io___237 = { 0, 0, 0, "(A)", 0 }; static cilist io___238 = { 0, 0, 0, "(A)", 0 }; static cilist io___239 = { 0, 0, 0, "(A)", 0 }; static cilist io___240 = { 0, 0, 0, "(A,F7.3,F7.3,F7.2,I8)", 0 }; static cilist io___243 = { 0, 0, 0, "(A,4X,A,F7.3,F7.3,F7.2,I8)", 0 }; static cilist io___244 = { 0, 0, 0, "(A/A,A/A/A)", 0 }; static cilist io___246 = { 0, 0, 0, fmt_146, 0 }; static cilist io___247 = { 0, 6, 0, fmt_151, 0 }; static cilist io___248 = { 0, 0, 0, 0, 0 }; static cilist io___249 = { 0, 6, 0, "(I2,/,A,I2,/,A)", 0 }; static cilist io___250 = { 0, 6, 0, fmt_171, 0 }; static cilist io___251 = { 0, 6, 0, fmt_172, 0 }; static cilist io___252 = { 0, 6, 0, fmt_172, 0 }; static cilist io___253 = { 0, 6, 0, fmt_172, 0 }; static cilist io___254 = { 0, 6, 0, fmt_151, 0 }; static cilist io___255 = { 0, 6, 0, fmt_152, 0 }; static cilist io___256 = { 0, 6, 0, fmt_151, 0 }; static cilist io___257 = { 0, 6, 0, "(A)", 0 }; static cilist io___258 = { 0, 6, 0, "(A)", 0 }; static cilist io___259 = { 0, 0, 0, 0, 0 }; static cilist io___260 = { 0, 0, 0, "(A,A,F5.3,A)", 0 }; static cilist io___273 = { 0, 0, 0, fmt_161, 0 }; static cilist io___276 = { 0, 6, 0, fmt_211, 0 }; static cilist io___283 = { 0, 6, 0, fmt_212, 0 }; static cilist io___284 = { 0, 6, 0, fmt_211, 0 }; static cilist io___286 = { 0, 6, 0, fmt_211, 0 }; static cilist io___287 = { 0, 6, 0, fmt_211, 0 }; static cilist io___288 = { 0, 0, 0, 0, 0 }; static cilist io___289 = { 0, 0, 0, "(/)", 0 }; static cilist io___290 = { 0, 0, 0, fmt_156, 0 }; static cilist io___291 = { 0, 0, 0, fmt_156, 0 }; static cilist io___292 = { 0, 0, 0, fmt_156, 0 }; static cilist io___293 = { 0, 0, 0, fmt_156, 0 }; /* ******************************************************************************* */ /* Usage: */ /* rastep [-h] [-iso] [-Bcolor Bmin Bmax] [-prob xx] [-radius r] [-fancy[0-9]] */ /* [-tabulate histogram.file] [-by_atomtype] [-suv_check] */ /* -auto auto-orientation of viewpoint */ /* -h suppresses header records in output */ /* -iso forces isotropic B values (spheres rather than */ /* ellipsoids) even if ANISOU cards present */ /* -Bcolor Bmin Bmax */ /* color by Biso values; Bmin = dark blue, Bmax = light red */ /* -Acolor color by anisotropy; red < white (A=0.5) < green */ /* -prob xx draws ellipsoids to enclose this */ /* probability level (default = 0.50) */ /* -radius draws bonds with this radius in Angstroms */ /* (default = 0.10) */ /* -fancy[0-9] increasingly complex rendition of ellipsoids */ /* fancy0 [default] solid surface only */ /* fancy1 principle axes + transparent bounding ellipsoid */ /* fancy2 equatorial planes only */ /* fancy3 equatorial planes + transparent ellipsoid */ /* fancy4 longest principle axis only */ /* fancy5 for ORTEP lovers - one octant missing */ /* fancy6 same as fancy5, with missing octant colored grey */ /* =============================================================================== */ /* The following options are used by parvati scripts */ /* -tabulate [file]instead of creating a Raster3D input file, */ /* list all atoms with principle axes and anisotropy. */ /* Optionally write a histogram of anisotropy to speficied */ /* output file; otherwise output is to stderr */ /* output from -tabulate for this version of rastep */ /* ATOM RESNAME RESNUM EIGEN1 EIGEN2 EIGEN3 ANISOTROPY Uiso */ /* -com [file] find in shells from center of mass */ /* -nohydrogens don't plot hydrogens even if present */ /* -mini small size plot (176x208) with auto-orientation */ /* -suv_check use Suv to validate similarity of bonded ellipsoids */ /* ******************************************************************************* */ /* EAM Jul 97 - initial version */ /* EAM Dec 97 - version 2.4b release */ /* EAM Jan 98 - add tabulation option */ /* EAM May 98 - integrate with PARVATI script */ /* EAM Jun 99 - fix bug (lack of sqrt) in -iso processing, trap read errors */ /* -Acolor flag to color by anisotropy (not yet entirely satisfactory) */ /* EAM Jul 99 - V2.4l -tab output revised slightly, */ /* NPD ellipsoids colored magenta */ /* -nohydro flag to suppress drawing hydrogens */ /* -mini flag to generate smaller pictures */ /* don't draw bonds between atoms with different ALT flags */ /* EAM Aug 99 - V2.4m */ /* work harder at suppressing hydrogens */ /* add auto-orientation (NB: scaling is wrong in this case!) */ /* EAM Dec 99 - V2.5 */ /* clean up output formats a little */ /* EAM Jun 2000 - additional error reporting */ /* EAM Jul 2000 - apply Bcolor to bonds as well as atoms */ /* EAM Sep 2000 - Suv similarity test */ /* EAM Apr 2001 - V2.6 */ /* ORTEP_LIKE ellipsoids (one octant missing) */ /* error count */ /* EAM Feb 2002 - rework ANISOU and Suv processing to gain back some speed */ /* maybe I should have an array of iso/aniso flags to save */ /* time during testing? */ /* the rest of CARD() array can go too? */ /* all the tests on IF ATOM(I)(1:).eq.'ATOM' are now unneeded */ /* EAM May 2003 - expand default color table to allow for off-by-one atom names */ /* EAM Oct 2003 - trap and report 0 values for axial Uij components in input */ /* EAM May 2006 - initial arrays to 0 for gfortran */ /* I/O units for colour/co-ordinate input, specs output, user output */ /* Data structures used for auto-orientation */ /* Support for validation of similarity of bonded atoms */ /* Default to CPK colors and VDW radii */ /* Critical values for probability ellipsoids of a trivariate normal */ /* distribution. From Table 6.1 of ORTEP-III manual (Oak Ridge National */ /* Laboratory Report ORNL-6895, 1996). Tabulated below in increments of */ /* 2% in probability. Default contours enclose a probability level of */ /* 50% (critical value 1.5382). */ asscom_1.assout = 0; asscom_1.verbose = FALSE_; hflag = FALSE_; acflag = FALSE_; bcflag = FALSE_; tflag = FALSE_; atflag = FALSE_; comflag = FALSE_; suvflag = FALSE_; ellipses = TRUE_; hwhacky = FALSE_; nohydro = FALSE_; mini = FALSE_; auto__ = FALSE_; fancy = 0; prob = .5f; radius = .1f; nerrors = 0; /* Gfortran is nuts */ wsum = 0.; xsum = 0.; ysum = 0.; zsum = 0.; sum_a__ = 0.f; sum_b__ = 0.f; for (i__ = 1; i__ <= 20; ++i__) { histogram[i__ - 1] = 0; } for (i__ = 1; i__ <= 110; ++i__) { adist[i__ - 1] = 0.f; hdist[i__ - 1] = 0; } narg = iargc_(); i__ = 1; L5: getarg_(&i__, flags, (ftnlen)80); if (s_cmp(flags, "-help", (ftnlen)5, (ftnlen)5) == 0) { goto L701; } if (s_cmp(flags, "-debug", (ftnlen)6, (ftnlen)6) == 0) { asscom_1.verbose = TRUE_; } if (s_cmp(flags, "-h", (ftnlen)2, (ftnlen)2) == 0) { hflag = TRUE_; } if (s_cmp(flags, "-iso", (ftnlen)4, (ftnlen)4) == 0) { ellipses = FALSE_; } if (s_cmp(flags, "-rad", (ftnlen)4, (ftnlen)4) == 0) { ++i__; if (i__ > narg) { goto L701; } getarg_(&i__, flags, (ftnlen)80); i__1 = s_rsli(&io___31); if (i__1 != 0) { goto L701; } i__1 = do_lio(&c__4, &c__1, (char *)&radius, (ftnlen)sizeof(real)); if (i__1 != 0) { goto L701; } i__1 = e_rsli(); if (i__1 != 0) { goto L701; } if (radius < 0.f) { radius = 0.f; } } if (s_cmp(flags, "-pro", (ftnlen)4, (ftnlen)4) == 0) { ++i__; if (i__ > narg) { goto L701; } getarg_(&i__, flags, (ftnlen)80); i__1 = s_rsli(&io___32); if (i__1 != 0) { goto L701; } i__1 = do_lio(&c__4, &c__1, (char *)&prob, (ftnlen)sizeof(real)); if (i__1 != 0) { goto L701; } i__1 = e_rsli(); if (i__1 != 0) { goto L701; } if (prob <= 0.f) { s_stop("illegal probability level", (ftnlen)25); } /* If prob > 1 assume they meant it in percent */ if (prob > 1.f) { prob /= 100.f; } } if (s_cmp(flags, "-Acol", (ftnlen)5, (ftnlen)5) == 0) { acflag = TRUE_; bcflag = FALSE_; } if (s_cmp(flags, "-Bcol", (ftnlen)5, (ftnlen)5) == 0) { bcflag = TRUE_; acflag = FALSE_; ++i__; if (i__ > narg) { goto L701; } getarg_(&i__, flags, (ftnlen)80); i__1 = s_rsli(&io___33); if (i__1 != 0) { goto L701; } i__1 = do_lio(&c__4, &c__1, (char *)&bmin, (ftnlen)sizeof(real)); if (i__1 != 0) { goto L701; } i__1 = e_rsli(); if (i__1 != 0) { goto L701; } ++i__; if (i__ > narg) { goto L701; } getarg_(&i__, flags, (ftnlen)80); i__1 = s_rsli(&io___35); if (i__1 != 0) { goto L701; } i__1 = do_lio(&c__4, &c__1, (char *)&bmax, (ftnlen)sizeof(real)); if (i__1 != 0) { goto L701; } i__1 = e_rsli(); if (i__1 != 0) { goto L701; } } if (s_cmp(flags, "-fancy", (ftnlen)6, (ftnlen)6) == 0) { if (*(unsigned char *)&flags[6] == '0') { fancy = 0; } else if (*(unsigned char *)&flags[6] == '1') { fancy = 1; } else if (*(unsigned char *)&flags[6] == '2') { fancy = 2; } else if (*(unsigned char *)&flags[6] == '3') { fancy = 3; } else if (*(unsigned char *)&flags[6] == '4') { fancy = 4; } else if (*(unsigned char *)&flags[6] == '5') { fancy = 5; } else if (*(unsigned char *)&flags[6] == '6') { fancy = 6; } else { fancy = 1; } } if (s_cmp(flags, "-tab", (ftnlen)4, (ftnlen)4) == 0) { tflag = TRUE_; hflag = TRUE_; acflag = FALSE_; bcflag = FALSE_; fancy = 0; for (j = 1; j <= 300000; ++j) { anisi[j - 1] = 0.f; } hislun = 0; if (i__ >= narg) { goto L799; } i__1 = i__ + 1; getarg_(&i__1, flags, (ftnlen)80); if (*(unsigned char *)flags != '-') { hislun = 1; o__1.oerr = 0; o__1.ounit = hislun; o__1.ofnmlen = 80; o__1.ofnm = flags; o__1.orl = 0; o__1.osta = "UNKNOWN"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = 0; f_open(&o__1); } } if (s_cmp(flags, "-com", (ftnlen)4, (ftnlen)4) == 0) { comflag = TRUE_; comlun = 0; if (i__ >= narg) { goto L799; } i__1 = i__ + 1; getarg_(&i__1, flags, (ftnlen)80); if (*(unsigned char *)flags != '-') { comlun = 2; o__1.oerr = 0; o__1.ounit = comlun; o__1.ofnmlen = 80; o__1.ofnm = flags; o__1.orl = 0; o__1.osta = "UNKNOWN"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = 0; f_open(&o__1); } } if (s_cmp(flags, "-suv", (ftnlen)4, (ftnlen)4) == 0) { suvflag = TRUE_; suvlun = 0; suvlimit = .975f; if (i__ >= narg) { goto L799; } i__1 = i__ + 1; getarg_(&i__1, flags, (ftnlen)80); if (*(unsigned char *)flags != '-') { i__1 = s_rsli(&io___43); if (i__1 != 0) { goto L701; } i__1 = do_lio(&c__4, &c__1, (char *)&suvlimit, (ftnlen)sizeof( real)); if (i__1 != 0) { goto L701; } i__1 = e_rsli(); if (i__1 != 0) { goto L701; } } } if (s_cmp(flags, "-by_atom", (ftnlen)8, (ftnlen)8) == 0) { atflag = TRUE_; } if (s_cmp(flags, "-nohydro", (ftnlen)8, (ftnlen)8) == 0) { nohydro = TRUE_; } if (s_cmp(flags, "-mini", (ftnlen)8, (ftnlen)5) == 0) { mini = TRUE_; auto__ = TRUE_; } if (s_cmp(flags, "-auto", (ftnlen)8, (ftnlen)5) == 0) { auto__ = TRUE_; } ++i__; if (i__ <= narg) { goto L5; } goto L799; L701: s_wsle(&io___44); do_lio(&c__9, &c__1, "Raster3D Thermal Ellipsoid Program ", (ftnlen)35); do_lio(&c__9, &c__1, "V2.7d ", (ftnlen)8); e_wsle(); s_wsfe(&io___45); do_fio(&c__1, "syntax:", (ftnlen)7); e_wsfe(); s_wsfe(&io___46); do_fio(&c__1, "rastep\t[-h] [-iso] [-Bcolor Bmin Bmax] [-prob Plevel]", ( ftnlen)53); e_wsfe(); s_wsfe(&io___47); do_fio(&c__1, "\t[-fancy[0-6]] [-radius R] [-auto]", (ftnlen)34); e_wsfe(); s_wsfe(&io___48); do_fio(&c__1, "\t[-nohydrogens] [-suv [suv_limit]]", (ftnlen)34); e_wsfe(); s_wsfe(&io___49); do_fio(&c__1, "\t[-tabulate [tabfile]] [-by_atomtype] [-com [comfile]]", ( ftnlen)54); e_wsfe(); exit_(&c_n1); L799: /* Critical values for the radius corresponding to a sphere */ /* enclosing the requested probability level are taken from */ /* Table 6.1 of the ORTEP manual */ iprob = (prob + .01f) * 50.f; pradius = problevel[iprob - 1]; s_wsfe(&io___52); e_wsfe(); s_wsle(&io___53); do_lio(&c__9, &c__1, "Raster3D Thermal Ellipsoid Program ", (ftnlen)35); do_lio(&c__9, &c__1, "V2.7d ", (ftnlen)8); e_wsle(); s_wsle(&io___54); do_lio(&c__9, &c__1, "E A Merritt - 2 Feb 2002", (ftnlen)25); e_wsle(); s_wsfe(&io___55); e_wsfe(); if (! ellipses) { s_wsfe(&io___56); r__1 = (real) iprob / 50.f; do_fio(&c__1, (char *)&r__1, (ftnlen)sizeof(real)); e_wsfe(); } else { s_wsfe(&io___57); r__1 = (real) iprob / 50.f; do_fio(&c__1, (char *)&r__1, (ftnlen)sizeof(real)); e_wsfe(); } s_wsfe(&io___58); do_fio(&c__1, (char *)&pradius, (ftnlen)sizeof(real)); e_wsfe(); if (acflag) { s_wsle(&io___59); do_lio(&c__9, &c__1, "Atoms will be colored based on Anisotropy", ( ftnlen)41); e_wsle(); } if (bcflag) { s_wsle(&io___60); do_lio(&c__9, &c__1, "Atom colors will be assigned based on Biso", ( ftnlen)42); e_wsle(); s_wsle(&io___61); do_lio(&c__9, &c__1, " from dark blue = Bmin =", (ftnlen)27); do_lio(&c__4, &c__1, (char *)&bmin, (ftnlen)sizeof(real)); e_wsle(); s_wsle(&io___62); do_lio(&c__9, &c__1, " to light red = Bmax =", (ftnlen)27); do_lio(&c__4, &c__1, (char *)&bmax, (ftnlen)sizeof(real)); e_wsle(); umin = bmin / 78.956701824799993f; umax = bmax / 78.956701824799993f; umin = umin; umax = umax; } if (! hflag) { s_wsfe(&io___65); do_fio(&c__1, "Raster3D thermal ellipsoid program ", (ftnlen)35); do_fio(&c__1, "V2.7d ", (ftnlen)8); i__1 = (integer) (prob * 100.f + .5f); do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); do_fio(&c__1, "% probability bounds", (ftnlen)20); e_wsfe(); if (mini) { s_wsfe(&io___66); do_fio(&c__1, "22 26 tiles in x,y", (ftnlen)22); e_wsfe(); } else { s_wsfe(&io___67); do_fio(&c__1, "80 64 tiles in x,y", (ftnlen)22); e_wsfe(); } s_wsfe(&io___68); do_fio(&c__1, " 8 8 pixels (x,y) per tile", (ftnlen)31); e_wsfe(); s_wsfe(&io___69); do_fio(&c__1, "4 3x3 virtual pixels -> 2x2 pixels", (ftnlen) 42); e_wsfe(); s_wsfe(&io___70); do_fio(&c__1, "1 1 1 white background", (ftnlen)26); e_wsfe(); s_wsfe(&io___71); do_fio(&c__1, "F no, shadows are dorky", (ftnlen)31); e_wsfe(); s_wsfe(&io___72); do_fio(&c__1, "25 Phong power", (ftnlen)21); e_wsfe(); s_wsfe(&io___73); do_fio(&c__1, "0.15 secondary light contribution", (ftnlen)38); e_wsfe(); s_wsfe(&io___74); do_fio(&c__1, "0.05 ambient light contribution", (ftnlen)36); e_wsfe(); s_wsfe(&io___75); do_fio(&c__1, "0.25 specular reflection component", (ftnlen)39); e_wsfe(); s_wsfe(&io___76); do_fio(&c__1, "0.0 No perspective", (ftnlen)24); e_wsfe(); s_wsfe(&io___77); do_fio(&c__1, "1 1 1 main light source position", (ftnlen)36); e_wsfe(); } if (auto__) { aspect = .84615384615384615f; } else { aspect = 1.25f; } ncol = 0; natm = 0; nani = 0; nanis = 0; niso = 0; nhyd = 0; nonpos = 0; L10: i__1 = s_rsfe(&io___86); if (i__1 != 0) { goto L50; } i__1 = do_fio(&c__1, card, (ftnlen)80); if (i__1 != 0) { goto L50; } i__1 = e_rsfe(); if (i__1 != 0) { goto L50; } if (s_cmp(card, "COLO", (ftnlen)4, (ftnlen)4) == 0) { ++ncol; if (ncol > 5000) { s_wsle(&io___88); do_lio(&c__9, &c__1, "Colour table overflow. Increase ", (ftnlen) 33); do_lio(&c__9, &c__1, "MAXCOL and recompile.", (ftnlen)21); e_wsle(); s_stop("10", (ftnlen)2); } s_rsfi(&io___89); do_fio(&c__1, mask + (ncol - 1) * 24, (ftnlen)24); for (i__ = 1; i__ <= 3; ++i__) { do_fio(&c__1, (char *)&rgb[i__ + ncol * 3 - 4], (ftnlen)sizeof( real)); } do_fio(&c__1, (char *)&vdw[ncol - 1], (ftnlen)sizeof(real)); e_rsfi(); } else if (nohydro && s_cmp(card + 76, " H", (ftnlen)2, (ftnlen)2) == 0) { goto L10; } else if (s_cmp(card, "ANISOU", (ftnlen)6, (ftnlen)6) == 0) { ++nani; if (s_cmp(card + 12, atom + ((natm - 1) * 80 + 12), (ftnlen)15, ( ftnlen)15) != 0) { goto L14; } i__1 = s_rsli(&io___94); if (i__1 != 0) { goto L12; } for (i__ = 1; i__ <= 6; ++i__) { i__1 = do_lio(&c__4, &c__1, (char *)&uij[i__ + natm * 6 - 7], ( ftnlen)sizeof(real)); if (i__1 != 0) { goto L12; } } i__1 = e_rsli(); if (i__1 != 0) { goto L12; } for (i__ = 1; i__ <= 6; ++i__) { uij[i__ + natm * 6 - 7] *= 1e-4f; } noerr = 1; for (i__ = 1; i__ <= 3; ++i__) { if (uij[i__ + natm * 6 - 7] == 0.f) { uij[i__ + natm * 6 - 7] = 1e-4f; noerr = 0; } } if (noerr == 0) { goto L15; } } else if (s_cmp(card, "ATOM", (ftnlen)4, (ftnlen)4) == 0 || s_cmp(card, "HETA", (ftnlen)4, (ftnlen)4) == 0) { ++natm; if (natm > 300000) { s_wsle(&io___97); do_lio(&c__9, &c__1, "Atom array overflow. Increase ", (ftnlen) 31); do_lio(&c__9, &c__1, "MAXATM and recompile.", (ftnlen)21); e_wsle(); s_stop("20", (ftnlen)2); } s_copy(atom + (natm - 1) * 80, card, (ftnlen)80, (ftnlen)80); uij[natm * 6 - 6] = -1.f; } else if (s_cmp(card, "END", (ftnlen)3, (ftnlen)3) == 0) { goto L50; } goto L10; L12: s_wsle(&io___98); do_lio(&c__9, &c__1, "*** Format problem - ", (ftnlen)21); do_lio(&c__9, &c__1, card + 12, (ftnlen)58); e_wsle(); ++nerrors; goto L10; L14: s_wsle(&io___99); do_lio(&c__9, &c__1, "*** ANISOU record out of order - ", (ftnlen)33); do_lio(&c__9, &c__1, card + 12, (ftnlen)58); e_wsle(); ++nerrors; goto L10; L15: s_wsle(&io___100); do_lio(&c__9, &c__1, "*** Illegal ANISOU values - ", (ftnlen)28); do_lio(&c__9, &c__1, card + 12, (ftnlen)58); e_wsle(); ++nerrors; goto L10; /* Come here when EOF or 'END' record is reached */ L50: if (natm == 0) { s_wsle(&io___101); do_lio(&c__9, &c__1, "No atoms in input.", (ftnlen)18); e_wsle(); s_stop("30", (ftnlen)2); } /* Load default colors after any that were read in */ if (ncol <= 4983) { for (i__ = 1; i__ <= 17; ++i__) { ++ncol; ici__1.icierr = 0; ici__1.iciend = 0; ici__1.icirnum = 1; ici__1.icirlen = 60; ici__1.iciunit = defcol + (i__ - 1) * 60; ici__1.icifmt = "(6X,A24,3F8.3,F6.2)"; s_rsfi(&ici__1); do_fio(&c__1, mask + (ncol - 1) * 24, (ftnlen)24); for (j = 1; j <= 3; ++j) { do_fio(&c__1, (char *)&rgb[j + ncol * 3 - 4], (ftnlen)sizeof( real)); } do_fio(&c__1, (char *)&vdw[ncol - 1], (ftnlen)sizeof(real)); e_rsfi(); } } if (ncol == 0) { s_wsle(&io___102); do_lio(&c__9, &c__1, "No colours in input.", (ftnlen)20); e_wsle(); s_stop("40", (ftnlen)2); } xmax = -1e20f; xmin = 1e20f; ymax = -1e20f; ymin = 1e20f; zmax = -1e20f; zmin = 1e20f; i__1 = natm; for (iatm = 1; iatm <= i__1; ++iatm) { s_copy(card, atom + (iatm - 1) * 80, (ftnlen)80, (ftnlen)80); /* Do a little pre-processing to make later bookkeeping easier */ /* At least screen out non-conformant PDB files that contain */ /* something obviously not an element type in columns 77:78 */ /* Hydrogen naming conventions are totally messed up. */ if (atflag) { s_copy(resname, card + 17, (ftnlen)3, (ftnlen)3); /* if (resname.eq.'MET' .and. card(14:15).eq.'SD') then */ /* atom(iatm)(77:78) = 'SD' */ /* end if */ if (s_cmp(resname, "HOH", (ftnlen)3, (ftnlen)3) == 0 || s_cmp( resname, "H2O", (ftnlen)3, (ftnlen)3) == 0 || s_cmp( resname, "WAT", (ftnlen)3, (ftnlen)3) == 0) { s_copy(atom + ((iatm - 1) * 80 + 76), "OW", (ftnlen)2, ( ftnlen)2); } if (*(unsigned char *)&atom[(iatm - 1) * 80 + 77] >= '0' && *( unsigned char *)&atom[(iatm - 1) * 80 + 77] <= '9') { s_copy(atom + ((iatm - 1) * 80 + 76), " ", (ftnlen)2, ( ftnlen)2); } } if (nohydro && s_cmp(atom + ((iatm - 1) * 80 + 76), " ", (ftnlen)2, ( ftnlen)2) == 0) { for (i__ = 13; i__ <= 16; ++i__) { if (*(unsigned char *)&atom[(iatm - 1) * 80 + (i__ - 1)] == ' ') { goto L70; } if (*(unsigned char *)&atom[(iatm - 1) * 80 + (i__ - 1)] >= '1' && *(unsigned char *)&atom[(iatm - 1) * 80 + (i__ - 1)] <= '4') { goto L70; } if (*(unsigned char *)&atom[(iatm - 1) * 80 + (i__ - 1)] != 'H') { goto L71; } s_copy(atom + ((iatm - 1) * 80 + 76), " H", (ftnlen)2, ( ftnlen)2); L70: ; } L71: ; } s_copy(test, card + 6, (ftnlen)24, (ftnlen)24); i__2 = ncol; for (icol = 1; icol <= i__2; ++icol) { if (match_(test, mask + (icol - 1) * 24, (ftnlen)24, (ftnlen)24)) { i__3 = s_rsfi(&io___113); if (i__3 != 0) { goto L82; } i__3 = do_fio(&c__1, (char *)&x, (ftnlen)sizeof(real)); if (i__3 != 0) { goto L82; } i__3 = do_fio(&c__1, (char *)&y, (ftnlen)sizeof(real)); if (i__3 != 0) { goto L82; } i__3 = do_fio(&c__1, (char *)&z__, (ftnlen)sizeof(real)); if (i__3 != 0) { goto L82; } i__3 = do_fio(&c__1, (char *)&biso[iatm - 1], (ftnlen)sizeof( real)); if (i__3 != 0) { goto L82; } i__3 = e_rsfi(); if (i__3 != 0) { goto L82; } if (biso[iatm - 1] <= 0.f) { ++nerrors; s_wsle(&io___118); do_lio(&c__9, &c__1, "*** Illegal Biso ", (ftnlen)17); do_lio(&c__4, &c__1, (char *)&biso[iatm - 1], (ftnlen) sizeof(real)); do_lio(&c__9, &c__1, " - ", (ftnlen)3); do_lio(&c__9, &c__1, atom + ((iatm - 1) * 80 + 12), ( ftnlen)15); e_wsle(); biso[iatm - 1] = 0.f; } uiso = biso[iatm - 1] / 78.956701824799993f; spam[iatm * 5 - 5] = x; spam[iatm * 5 - 4] = y; spam[iatm * 5 - 3] = z__; spam[iatm * 5 - 2] = uiso; spam[iatm * 5 - 1] = (real) icol; rad = sqrt(uiso) * pradius; /* Computing MAX */ r__1 = xmax, r__2 = x + rad; xmax = dmax(r__1,r__2); /* Computing MIN */ r__1 = xmin, r__2 = x - rad; xmin = dmin(r__1,r__2); /* Computing MAX */ r__1 = ymax, r__2 = y + rad; ymax = dmax(r__1,r__2); /* Computing MIN */ r__1 = ymin, r__2 = y - rad; ymin = dmin(r__1,r__2); /* Computing MAX */ r__1 = zmax, r__2 = z__ + rad; zmax = dmax(r__1,r__2); /* Computing MIN */ r__1 = zmin, r__2 = z__ - rad; zmin = dmin(r__1,r__2); /* atomtype = CARD(77:78) */ /* if (atomtype.ne.' ') then */ /* weight = amass(atomtype) */ /* else */ weight = 13.4f; wsum += weight; xsum += weight * x; ysum += weight * y; zsum += weight * z__; goto L100; } /* L80: */ } s_wsle(&io___123); do_lio(&c__9, &c__1, "No colour table mask matches this atom:", ( ftnlen)39); e_wsle(); s_wsle(&io___124); do_lio(&c__9, &c__1, atom + (iatm - 1) * 80, (ftnlen)80); e_wsle(); s_stop("90", (ftnlen)2); L82: s_wsle(&io___125); do_lio(&c__9, &c__1, "Input format problem in record", (ftnlen)30); e_wsle(); s_wsle(&io___126); do_lio(&c__9, &c__1, card, (ftnlen)80); e_wsle(); s_stop("90", (ftnlen)2); L100: ; } xmid = (xmax + xmin) / 2.f; ymid = (ymax + ymin) / 2.f; zmid = (zmax + zmin) / 2.f; tx = -xmid; ty = -ymid; tz = -zmid; xcom = xsum / wsum; ycom = ysum / wsum; zcom = zsum / wsum; if (aspect >= 1.f) { /* The X direction is wider than the Y */ xroom = aspect; yroom = 1.f; zroom = 2.f; } else { xroom = 1.f; yroom = aspect; zroom = 2.f; } xspan = xmax - xmin; yspan = ymax - ymin; zspan = zmax - zmin; /* Computing MAX */ r__1 = xspan / xroom, r__2 = yspan / yroom, r__1 = max(r__1,r__2), r__2 = zspan / zroom; scale = dmax(r__1,r__2); /* Leave a little extra room as a border: */ scale /= .9f; if (mini) { /* Computing 2nd power */ r__1 = xspan; /* Computing 2nd power */ r__2 = yspan; /* Computing 2nd power */ r__3 = zspan; scale = sqrt(r__1 * r__1 + r__2 * r__2 + r__3 * r__3) * aspect; if (scale < 9.f) { scale = 9.f; } } /* These are for the center-of-mass table */ /* Computing MAX */ d__8 = (d__1 = xmax - xcom, abs(d__1)), d__9 = (d__2 = xmin - xcom, abs( d__2)); /* Computing 2nd power */ d__7 = max(d__8,d__9); /* Computing MAX */ d__11 = (d__3 = ymax - ycom, abs(d__3)), d__12 = (d__4 = ymin - ycom, abs( d__4)); /* Computing 2nd power */ d__10 = max(d__11,d__12); /* Computing MAX */ d__14 = (d__5 = zmax - zcom, abs(d__5)), d__15 = (d__6 = zmin - zcom, abs( d__6)); /* Computing 2nd power */ d__13 = max(d__14,d__15); dmax__ = d__7 * d__7 + d__10 * d__10 + d__13 * d__13; dmax__ = sqrt(dmax__); dshells = 100; if (! hflag) { s_wsfe(&io___145); do_fio(&c__1, (char *)&tx, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&ty, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&tz, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&scale, (ftnlen)sizeof(real)); e_wsfe(); s_wsfe(&io___146); do_fio(&c__1, "3 mixed object types", (ftnlen)28); e_wsfe(); s_wsfe(&io___147); do_fio(&c__1, "*", (ftnlen)1); e_wsfe(); s_wsfe(&io___148); do_fio(&c__1, "*", (ftnlen)1); e_wsfe(); s_wsfe(&io___149); do_fio(&c__1, "*", (ftnlen)1); e_wsfe(); } /* Auto-orientation */ /* 25-Aug-1999 */ /* Find Eigenvectors of moment of inertia tensor. */ /* Arrange smallest Eigenvalue along Y, largest along Z. */ /* Problems: */ /* - Could emphasize side-chain over backbone by ignoring */ /* atoms O and N, but in practice this doesn't seem to help much. */ /* - Scaling is wrong, because it was done before rotation. */ if (auto__) { i__1 = natm; for (iatm = 1; iatm <= i__1; ++iatm) { if (s_cmp(atom + (iatm - 1) * 80, "ATOM", (ftnlen)4, (ftnlen)4) != 0 && s_cmp(atom + (iatm - 1) * 80, "HETA", (ftnlen)4, ( ftnlen)4) != 0) { goto L125; } /* if (atom(iatm)(13:15).eq.' O ') goto 125 */ /* if (atom(iatm)(13:15).eq.' N ') goto 125 */ if (s_cmp(atom + ((iatm - 1) * 80 + 76), " H", (ftnlen)2, (ftnlen) 2) == 0 && nohydro) { goto L125; } x = spam[iatm * 5 - 5] - xcom; y = spam[iatm * 5 - 4] - ycom; z__ = spam[iatm * 5 - 3] - zcom; rq = x * x + y * y + z__ * z__; rv = sqrt(rq); rr[0] += x * x; rr[3] += x * y; rr[6] += x * z__; rr[1] += y * x; rr[4] += y * y; rr[7] += y * z__; rr[2] += z__ * x; rr[5] += z__ * y; rr[8] += z__ * z__; xmom[0] += 1.f; xmom[1] += rv; xmom[2] += rq; xmom[3] += rv * rq; xmom[4] += rq * rq; L125: ; } if (asscom_1.verbose) { s_wsfe(&io___154); do_fio(&c__1, " Radial moments:", (ftnlen)16); do_fio(&c__5, (char *)&xmom[0], (ftnlen)sizeof(real)); e_wsfe(); s_wsfe(&io___155); do_fio(&c__1, " Moment of inertia tensor", (ftnlen)25); e_wsfe(); } rg = sqrt(xmom[2] / xmom[0]); u[0] = xmom[2]; u[5] = xmom[2]; u[10] = xmom[2]; for (k = 1; k <= 3; ++k) { for (l = 1; l <= 3; ++l) { u[k + (l << 2) - 5] = (u[k + (l << 2) - 5] - rr[k + l * 3 - 4] ) / xmom[0]; } if (asscom_1.verbose) { s_wsfe(&io___160); for (l = 1; l <= 3; ++l) { do_fio(&c__1, (char *)&u[k + (l << 2) - 5], (ftnlen) sizeof(real)); } e_wsfe(); } } jacobi_(u, &c__3, &c__4, eigens, evecs); /* Re-order so that long axis is vertical, short axis towards viewer */ kmax = 1; if (eigens[1] > eigens[kmax - 1]) { kmax = 2; } if (eigens[2] > eigens[kmax - 1]) { kmax = 3; } kmin = 1; if (eigens[1] <= eigens[kmin - 1]) { kmin = 2; } if (eigens[2] <= eigens[kmin - 1]) { kmin = 3; } kmid = 6 - (kmin + kmax); if (asscom_1.verbose) { s_wsfe(&io___166); do_fio(&c__1, " Eigenvalues:", (ftnlen)13); do_fio(&c__1, (char *)&eigens[kmin - 1], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&eigens[kmid - 1], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&eigens[kmax - 1], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&kmin, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&kmid, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&kmax, (ftnlen)sizeof(integer)); e_wsfe(); s_wsfe(&io___167); do_fio(&c__1, " Eigenvectors:", (ftnlen)14); e_wsfe(); for (k = 1; k <= 3; ++k) { s_wsfe(&io___168); do_fio(&c__1, (char *)&evecs[k + (kmin << 2) - 5], (ftnlen) sizeof(real)); do_fio(&c__1, (char *)&evecs[k + (kmid << 2) - 5], (ftnlen) sizeof(real)); do_fio(&c__1, (char *)&evecs[k + (kmax << 2) - 5], (ftnlen) sizeof(real)); e_wsfe(); } } for (k = 1; k <= 3; ++k) { u[k - 1] = evecs[k + (kmid << 2) - 5]; u[k + 3] = evecs[k + (kmin << 2) - 5]; u[k + 7] = evecs[k + (kmax << 2) - 5]; u[k + 11] = 0.f; u[(k << 2) - 1] = 0.f; } u[15] = 1.f; /* Beware! may be left-handed at this point! */ det = tinv4_(evecinv, u); if (det < 0.f) { for (k = 1; k <= 3; ++k) { u[k + 7] = -u[k + 7]; } } /* OK, now it should be right-handed */ s_wsfe(&io___171); do_fio(&c__1, "# Auto-orientation matrix", (ftnlen)25); e_wsfe(); s_wsfe(&io___172); do_fio(&c__1, "16", (ftnlen)2); e_wsfe(); s_wsfe(&io___173); do_fio(&c__1, "ROTATION", (ftnlen)8); e_wsfe(); for (k = 1; k <= 3; ++k) { s_wsfe(&io___174); do_fio(&c__1, (char *)&u[(k << 2) - 4], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&u[(k << 2) - 3], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&u[(k << 2) - 2], (ftnlen)sizeof(real)); e_wsfe(); } s_wsfe(&io___175); do_fio(&c__1, "# End auto-orientation", (ftnlen)22); e_wsfe(); } /* Label output records */ if (! tflag) { s_wsfe(&io___176); do_fio(&c__1, "# Thermal ellipsoids from Rastep Version ", (ftnlen)41) ; do_fio(&c__1, "V2.7d ", (ftnlen)8); e_wsfe(); s_wsfe(&io___177); do_fio(&c__1, "# Probability level", (ftnlen)19); r__1 = (real) iprob / 50.f; do_fio(&c__1, (char *)&r__1, (ftnlen)sizeof(real)); e_wsfe(); } /* Write ellipsoids to input file for render */ if (fancy == 0 && ! tflag) { goto L139; } /* First, optional pass, to write fancy stuff associated with ellipsoids */ iatm = 1; L130: if (s_cmp(atom + (iatm - 1) * 80, "ATOM", (ftnlen)4, (ftnlen)4) == 0 || s_cmp(atom + (iatm - 1) * 80, "HETA", (ftnlen)4, (ftnlen)4) == 0) { x = spam[iatm * 5 - 5]; y = spam[iatm * 5 - 4]; z__ = spam[iatm * 5 - 3]; icol = spam[iatm * 5 - 1]; if (bcflag) { u2rgb_(&spam[iatm * 5 - 2], &umin, &umax, &red, &green, &blue); red *= red; green *= green; blue *= blue; } else if (acflag) { a2rgb_(&c_b332, &red, &green, &blue); red *= red; green *= green; blue *= blue; } else { red = rgb[icol * 3 - 3]; green = rgb[icol * 3 - 2]; blue = rgb[icol * 3 - 1]; } if (ellipses && uij[iatm * 6 - 6] >= 0.f) { for (i__ = 1; i__ <= 6; ++i__) { anisou[i__ - 1] = uij[i__ + iatm * 6 - 7]; } if (anitoquad_(anisou, &pradius, quadric, eigens, evecs) < 0) { s_wsle(&io___183); do_lio(&c__9, &c__1, "*** Non-positive definite ellipsoid - ", (ftnlen)38); do_lio(&c__9, &c__1, atom + ((iatm - 1) * 80 + 12), (ftnlen) 15); e_wsle(); ++nonpos; ++nerrors; biso[iatm - 1] = 0.f; goto L138; } goto L132; L132: /* Computing MAX */ r__1 = max(eigens[0],eigens[1]); radlim = pradius * dmax(r__1,eigens[2]); radlim *= 1.15f; /* Only for debugging ellipsoids */ if (asscom_1.verbose) { s_wsfe(&io___185); do_fio(&c__1, "ANISOU ", (ftnlen)7); do_fio(&c__1, (char *)&x, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&y, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&z__, (ftnlen)sizeof(real)); do_fio(&c__6, (char *)&anisou[0], (ftnlen)sizeof(real)); e_wsfe(); s_wsfe(&io___186); do_fio(&c__1, "QUADRIC", (ftnlen)7); do_fio(&c__10, (char *)&quadric[0], (ftnlen)sizeof(real)); e_wsfe(); s_wsfe(&io___187); do_fio(&c__1, "Eigenvalues", (ftnlen)11); for (i__ = 1; i__ <= 3; ++i__) { do_fio(&c__1, (char *)&eigens[i__ - 1], (ftnlen)sizeof( real)); } do_fio(&c__1, "prob", (ftnlen)4); do_fio(&c__1, (char *)&prob, (ftnlen)sizeof(real)); do_fio(&c__1, "limiting radius", (ftnlen)15); do_fio(&c__1, (char *)&radlim, (ftnlen)sizeof(real)); e_wsfe(); s_wsfe(&io___188); do_fio(&c__1, "Evecs ", (ftnlen)6); for (j = 1; j <= 3; ++j) { for (i__ = 1; i__ <= 3; ++i__) { do_fio(&c__1, (char *)&evecs[i__ + (j << 2) - 5], ( ftnlen)sizeof(real)); } } e_wsfe(); } /* Tabulate principal axes of ellipsoid for each atom */ if (tflag) { for (i__ = 1; i__ <= 3; ++i__) { /* Computing 2nd power */ r__1 = eigens[i__ - 1]; uprin[i__ - 1] = r__1 * r__1; } if (uprin[1] > uprin[0]) { umean = uprin[0]; uprin[0] = uprin[1]; uprin[1] = umean; } if (uprin[2] > uprin[0]) { umean = uprin[0]; uprin[0] = uprin[2]; uprin[2] = umean; } if (uprin[2] > uprin[1]) { umean = uprin[1]; uprin[1] = uprin[2]; uprin[2] = umean; } /* Anisotropy we define as Umin / Umax */ /* as in shelxpro output */ /* Computing MIN */ r__1 = min(uprin[0],uprin[1]); /* Computing MAX */ r__2 = max(uprin[0],uprin[1]); anisotropy = dmin(r__1,uprin[2]) / dmax(r__2,uprin[2]); /* But don't count atoms which are perfectly isotropic */ if (s_cmp(atom + ((iatm - 1) * 80 + 76), " H", (ftnlen)2, ( ftnlen)2) == 0) { ++nhyd; if (anisotropy != 1.f) { hwhacky = TRUE_; } } else if (anisotropy == 1.f) { ++niso; } else { anisi[iatm - 1] = anisotropy; sum_a__ += anisotropy; sum_b__ += biso[iatm - 1]; ++nanis; } /* Ellipticity we define as 1 / anisotropy */ if (anisotropy == 0.f) { ellipticity = 0.f; } else { ellipticity = 1.f / anisotropy; } /* Longhi et al (1997) JMB 268, 779-799. */ /* proposed another measure A = sigU / meanU */ umean = (uprin[0] + uprin[1] + uprin[2]) / 3.f; /* Computing 2nd power */ r__1 = uprin[0] - umean; /* Computing 2nd power */ r__2 = uprin[1] - umean; /* Computing 2nd power */ r__3 = uprin[2] - umean; usigma = r__1 * r__1 + r__2 * r__2 + r__3 * r__3; usigma = sqrt(usigma) / 3.f; alonghi = usigma / umean; /* Might want to check correlation with Uiso */ uiso = spam[iatm * 5 - 2]; /* Cosmetic changes to atom identifier for the sake of sorting */ /* We will force there to be exactly three entities printed. */ /* PDB format is just a mess: */ /* cols 13:16 atom */ /* col 17 alternate conf */ /* cols 18:20 residue */ /* col 22 chain */ /* cols 23:27 resnum */ for (i__ = 16; i__ >= 13; --i__) { if (*(unsigned char *)&atom[(iatm - 1) * 80 + (i__ - 1)] != ' ') { j = i__; } } for (i__ = 13; i__ <= 17; ++i__) { if (*(unsigned char *)&atom[(iatm - 1) * 80 + (i__ - 1)] != ' ') { k = i__; } } i__1 = k; for (i__ = j; i__ <= i__1; ++i__) { if (*(unsigned char *)&atom[(iatm - 1) * 80 + (i__ - 1)] == ' ') { *(unsigned char *)&atom[(iatm - 1) * 80 + (i__ - 1)] = '_'; } } /* if (ATOM(iatm)(17:17) .ne. ' ') then */ /* do i = 18,19 */ /* if (ATOM(iatm)(i:i) .eq. ' ') ATOM(iatm)(i:i) = '_' */ /* enddo */ /* endif */ *(unsigned char *)spacer = ' '; if (*(unsigned char *)&atom[(iatm - 1) * 80 + 21] != ' ') { *(unsigned char *)spacer = '_'; for (i__ = 23; i__ <= 25; ++i__) { if (*(unsigned char *)&atom[(iatm - 1) * 80 + (i__ - 1)] == ' ') { *(unsigned char *)&atom[(iatm - 1) * 80 + (i__ - 1)] = '_'; } } } s_wsfe(&io___196); do_fio(&c__1, atom + ((iatm - 1) * 80 + 12), (ftnlen)5); do_fio(&c__1, atom + ((iatm - 1) * 80 + 17), (ftnlen)5); do_fio(&c__1, spacer, (ftnlen)1); do_fio(&c__1, atom + ((iatm - 1) * 80 + 22), (ftnlen)5); do_fio(&c__1, (char *)&uprin[0], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&uprin[1], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&uprin[2], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&anisotropy, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&uiso, (ftnlen)sizeof(real)); e_wsfe(); /* Also make a histogram of anisotropies */ i__ = anisotropy * 20.f + 1; ++histogram[i__ - 1]; /* And a table of by distance from center of mass */ if (comflag && anisotropy < 1.f) { /* Computing 2nd power */ d__1 = spam[iatm * 5 - 5] - xcom; /* Computing 2nd power */ d__2 = spam[iatm * 5 - 4] - ycom; /* Computing 2nd power */ d__3 = spam[iatm * 5 - 3] - zcom; dist = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3); i__ = dist / dmax__ * (real) dshells + 1; adist[i__ - 1] += anisotropy; ++hdist[i__ - 1]; } /* End tabulation code */ } /* Skip hydrogens if requested */ if (nohydro && s_cmp(atom + ((iatm - 1) * 80 + 76), " H", (ftnlen) 2, (ftnlen)2) == 0) { goto L138; } /* Draw principal axes inside bounding ellipsoid */ if (fancy == 1) { for (i__ = 1; i__ <= 3; ++i__) { size = eigens[i__ - 1] * pradius - .02f; start[0] = x - size * evecs[(i__ << 2) - 4]; start[1] = y - size * evecs[(i__ << 2) - 3]; start[2] = z__ - size * evecs[(i__ << 2) - 2]; end[0] = x + size * evecs[(i__ << 2) - 4]; end[1] = y + size * evecs[(i__ << 2) - 3]; end[2] = z__ + size * evecs[(i__ << 2) - 2]; s_wsfe(&io___201); do_fio(&c__3, (char *)&start[0], (ftnlen)sizeof(real)); do_fio(&c__3, (char *)&end[0], (ftnlen)sizeof(real)); e_wsfe(); } } /* Draw longest principle axis only */ /* (experimental use only - not supported or documented) */ if (fancy == 4) { imax = 1; if (eigens[1] > eigens[imax - 1]) { imax = 2; } if (eigens[2] > eigens[imax - 1]) { imax = 3; } size = eigens[imax - 1] * pradius; imin = 1; if (eigens[1] < eigens[imin - 1]) { imin = 2; } if (eigens[2] < eigens[imin - 1]) { imin = 3; } imed = 6 - (imax + imin); size = size * eigens[imax - 1] / eigens[imed - 1]; start[0] = x - size * evecs[(imax << 2) - 4]; start[1] = y - size * evecs[(imax << 2) - 3]; start[2] = z__ - size * evecs[(imax << 2) - 2]; end[0] = x + size * evecs[(imax << 2) - 4]; end[1] = y + size * evecs[(imax << 2) - 3]; end[2] = z__ + size * evecs[(imax << 2) - 2]; s_wsfe(&io___205); do_fio(&c__3, (char *)&start[0], (ftnlen)sizeof(real)); do_fio(&c__3, (char *)&end[0], (ftnlen)sizeof(real)); e_wsfe(); } /* Construct 3 quadrics corresponding to the 3 orthogonal planes */ /* through the center of our ellipsoid */ if (fancy == 2 || fancy == 3) { eigens[3] = 1.f; evecs[15] = 1.f; det = tinv4_(evecinv, evecs); trnsp4_(evecit, evecinv); for (k = 1; k <= 3; ++k) { for (i__ = 1; i__ <= 4; ++i__) { for (j = 1; j <= 4; ++j) { qq[i__ + (j << 2) - 5] = 0.f; } qq[i__ + (i__ << 2) - 5] = 1.f / (eigens[i__ - 1] * eigens[i__ - 1]); } qq[k + (k << 2) - 5] = 1e3f; qq[15] = -pradius * pradius; tmul4_(temp, qq, evecinv); tmul4_(qp, evecit, temp); if (acflag) { a2rgb_(&anisotropy, &red, &green, &blue); red *= red; green *= green; blue *= blue; } s_wsfe(&io___210); do_fio(&c__1, (char *)&c__14, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&x, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&y, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&z__, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&radlim, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&red, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&green, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&blue, (ftnlen)sizeof(real)); e_wsfe(); s_wsfe(&io___211); do_fio(&c__1, (char *)&qp[0], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&qp[5], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&qp[10], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&qp[4], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&qp[9], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&qp[8], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&qp[12], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&qp[13], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&qp[14], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&qp[15], (ftnlen)sizeof(real)); e_wsfe(); } } } } L138: ++iatm; if (iatm <= natm) { goto L130; } /* Set transparency for enclosing ellipoids */ if (fancy == 1 || fancy == 3) { s_wsfe(&io___212); do_fio(&c__1, "9 Begin transparent ellipsoids", (ftnlen)30); do_fio(&c__1, "8 ", (ftnlen)2); e_wsfe(); s_wsfe(&io___213); do_fio(&c__1, " 15. 0.6 1.0 1.0 1.0 0.6 0 0 0 0", (ftnlen)39); e_wsfe(); } else if (fancy == 2 || fancy == 4) { goto L160; } L139: /* If we're just tabulating ellipticity, start making tables */ if (tflag) { total = 0.f; for (i__ = 1; i__ <= 20; ++i__) { total += histogram[i__ - 1]; } if (total == 0.f) { s_wsle(&io___215); do_lio(&c__9, &c__1, "No ANISOU records found", (ftnlen)23); e_wsle(); exit_(&c_n1); } if (nanis == 0) { s_wsle(&io___216); do_lio(&c__9, &c__1, "No anisotropic atoms found", (ftnlen)26); e_wsle(); exit_(&c_n1); } if (hwhacky) { s_wsle(&io___217); do_lio(&c__9, &c__1, "You seem to have anisotropic hydrogens", ( ftnlen)38); do_lio(&c__9, &c__1, " - is this some kind of joke?", (ftnlen)29); e_wsle(); ++nerrors; } io___218.ciunit = hislun; s_wsfe(&io___218); do_fio(&c__1, "# Anisotropy Fraction Number", (ftnlen)31); e_wsfe(); io___219.ciunit = hislun; s_wsfe(&io___219); do_fio(&c__1, "# range of atoms of atoms", (ftnlen)31); e_wsfe(); for (i__ = 1; i__ <= 20; ++i__) { io___220.ciunit = hislun; s_wsfe(&io___220); r__1 = ((real) i__ - 1.f) / 20.f; do_fio(&c__1, (char *)&r__1, (ftnlen)sizeof(real)); r__2 = (real) i__ / 20.f; do_fio(&c__1, (char *)&r__2, (ftnlen)sizeof(real)); r__3 = (real) histogram[i__ - 1] / total; do_fio(&c__1, (char *)&r__3, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&histogram[i__ - 1], (ftnlen)sizeof(integer) ); e_wsfe(); } /* Calculate mean and sigma of distribution */ sum_a2__ = 0.f; sum_b2__ = 0.f; sum_ab__ = 0.f; anis_mean__ = sum_a__ / (real) nanis; anis_sigma__ = 0.f; biso_mean__ = sum_b__ / (real) nanis; ccoef = 0.f; i__1 = natm; for (i__ = 1; i__ <= i__1; ++i__) { if (anisi[i__ - 1] != 0.f) { /* Computing 2nd power */ r__1 = anisi[i__ - 1] - anis_mean__; sum_a2__ += r__1 * r__1; /* Computing 2nd power */ r__1 = biso[i__ - 1] - biso_mean__; sum_b2__ += r__1 * r__1; sum_ab__ += (anisi[i__ - 1] - anis_mean__) * (biso[i__ - 1] - biso_mean__); } } if (nanis > 1) { anis_sigma__ = sqrt(sum_a2__ / (real) (nanis - 1)); ccoef = sum_ab__ / sqrt(sum_a2__ * sum_b2__); } io___228.ciunit = hislun; s_wsfe(&io___228); do_fio(&c__1, "#", (ftnlen)1); e_wsfe(); io___229.ciunit = hislun; s_wsfe(&io___229); do_fio(&c__1, "# number of ANISOU records:", (ftnlen)28); do_fio(&c__1, (char *)&nani, (ftnlen)sizeof(integer)); e_wsfe(); io___230.ciunit = hislun; s_wsfe(&io___230); do_fio(&c__1, "# non-isotropic atoms:", (ftnlen)28); do_fio(&c__1, (char *)&nanis, (ftnlen)sizeof(integer)); e_wsfe(); io___231.ciunit = hislun; s_wsfe(&io___231); do_fio(&c__1, "# isotropic atoms:", (ftnlen)28); do_fio(&c__1, (char *)&niso, (ftnlen)sizeof(integer)); e_wsfe(); if (nonpos > 0) { io___232.ciunit = hislun; s_wsfe(&io___232); do_fio(&c__1, "# nonpositive-definite APDs:", (ftnlen)28); do_fio(&c__1, (char *)&nonpos, (ftnlen)sizeof(integer)); e_wsfe(); } if (nhyd > 0) { io___233.ciunit = hislun; s_wsfe(&io___233); do_fio(&c__1, "# ANISOU hydrogens:", (ftnlen)28); do_fio(&c__1, (char *)&nhyd, (ftnlen)sizeof(integer)); e_wsfe(); } io___234.ciunit = hislun; s_wsfe(&io___234); do_fio(&c__1, "#", (ftnlen)1); e_wsfe(); io___235.ciunit = hislun; s_wsfe(&io___235); do_fio(&c__1, "# correlation between anisotropy and B_iso:", (ftnlen) 43); do_fio(&c__1, (char *)&ccoef, (ftnlen)sizeof(real)); e_wsfe(); io___236.ciunit = hislun; s_wsfe(&io___236); do_fio(&c__1, "#", (ftnlen)1); e_wsfe(); io___237.ciunit = hislun; s_wsfe(&io___237); do_fio(&c__1, "# Anisotropy B_iso", (ftnlen)32); e_wsfe(); io___238.ciunit = hislun; s_wsfe(&io___238); do_fio(&c__1, "# AtomType mean sigma mean number", (ftnlen)40); e_wsfe(); io___239.ciunit = hislun; s_wsfe(&io___239); do_fio(&c__1, "# --------- ----------- ------ ------", (ftnlen)40); e_wsfe(); io___240.ciunit = hislun; s_wsfe(&io___240); do_fio(&c__1, "#| Total", (ftnlen)11); do_fio(&c__1, (char *)&anis_mean__, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&anis_sigma__, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&biso_mean__, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&nanis, (ftnlen)sizeof(integer)); e_wsfe(); /* If we're tabulating ellipticity, but not by atom_type, then we're done */ if (! atflag) { goto L145; } L140: /* Find an atom type we haven't done yet, exit if none left */ /* This only works if the PDB file contains the atom type in */ /* columns 77:78 (standard since sometime in 1997, but many */ /* files do not conform to this standard) */ i__1 = natm; for (i__ = 1; i__ <= i__1; ++i__) { if (anisi[i__ - 1] != 0.f && s_cmp(atom + ((i__ - 1) * 80 + 76), " ", (ftnlen)2, (ftnlen)2) != 0) { s_copy(atomtype, atom + ((i__ - 1) * 80 + 76), (ftnlen)2, ( ftnlen)2); goto L141; } } goto L145; L141: sum_a__ = 0.f; sum_b__ = 0.f; sum_a2__ = 0.f; biso_mean__ = 0.f; natype = 0; /* Loop over atoms looking for the right ones */ i__1 = natm; for (i__ = 1; i__ <= i__1; ++i__) { if (anisi[i__ - 1] != 0.f && s_cmp(atomtype, atom + ((i__ - 1) * 80 + 76), (ftnlen)2, (ftnlen)2) == 0) { ++natype; sum_a__ += anisi[i__ - 1]; } } anis_mean__ = sum_a__ / (real) natype; anis_sigma__ = 0.f; i__1 = natm; for (i__ = 1; i__ <= i__1; ++i__) { if (anisi[i__ - 1] != 0.f && s_cmp(atomtype, atom + ((i__ - 1) * 80 + 76), (ftnlen)2, (ftnlen)2) == 0) { /* Computing 2nd power */ r__1 = anisi[i__ - 1] - anis_mean__; sum_a2__ += r__1 * r__1; biso_mean__ += biso[i__ - 1]; s_copy(atom + ((i__ - 1) * 80 + 76), " ", (ftnlen)2, (ftnlen) 2); } } if (natype > 1) { anis_sigma__ = sqrt(sum_a2__ / (real) (natype - 1)); } biso_mean__ /= (real) natype; io___243.ciunit = hislun; s_wsfe(&io___243); do_fio(&c__1, "#| ", (ftnlen)5); do_fio(&c__1, atomtype, (ftnlen)2); do_fio(&c__1, (char *)&anis_mean__, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&anis_sigma__, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&biso_mean__, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&natype, (ftnlen)sizeof(integer)); e_wsfe(); goto L140; } iatm = 1; goto L150; /* Write out plot of mean anisotropy as a function of distance from c.o.m. */ /* The hisclean routine is not strictly necessary; it collapses shells at */ /* the extreme ranges of distance that contain only a few atoms. */ L145: if (comflag) { io___244.ciunit = comlun; s_wsfe(&io___244); do_fio(&c__1, "#", (ftnlen)1); do_fio(&c__1, "# Mean anisotropy as a function of distance", (ftnlen) 43); do_fio(&c__1, " from center of mass", (ftnlen)20); do_fio(&c__1, "#", (ftnlen)1); do_fio(&c__1, "# Distance #atoms", (ftnlen)29); e_wsfe(); hisclean_(adist, hdist, &dshells, &nanis); i__1 = dshells; for (i__ = 1; i__ <= i__1; ++i__) { if (hdist[i__ - 1] > 0) { adist[i__ - 1] /= (real) hdist[i__ - 1]; dmid = dmax__ / (real) dshells * ((real) i__ - .5f); io___246.ciunit = comlun; s_wsfe(&io___246); do_fio(&c__1, (char *)&dmid, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&adist[i__ - 1], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&hdist[i__ - 1], (ftnlen)sizeof(integer) ); e_wsfe(); } } } if (suvflag) { goto L160; } exit_(&c__0); /* Second pass write a single sphere or ellipsoid for each atom */ L150: if (s_cmp(atom + (iatm - 1) * 80, "ATOM", (ftnlen)4, (ftnlen)4) == 0 || s_cmp(atom + (iatm - 1) * 80, "HETA", (ftnlen)4, (ftnlen)4) == 0) { if (nohydro && s_cmp(atom + ((iatm - 1) * 80 + 76), " H", (ftnlen)2, ( ftnlen)2) == 0) { goto L154; } x = spam[iatm * 5 - 5]; y = spam[iatm * 5 - 4]; z__ = spam[iatm * 5 - 3]; icol = spam[iatm * 5 - 1]; if (bcflag) { u2rgb_(&spam[iatm * 5 - 2], &umin, &umax, &red, &green, &blue); red *= red; green *= green; blue *= blue; } else if (acflag) { a2rgb_(&c_b332, &red, &green, &blue); red *= red; green *= green; blue *= blue; } else { red = rgb[icol * 3 - 3]; green = rgb[icol * 3 - 2]; blue = rgb[icol * 3 - 1]; } if (! ellipses) { rad = sqrt(spam[iatm * 5 - 2]) * pradius; s_wsfe(&io___247); do_fio(&c__1, (char *)&c__2, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&x, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&y, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&z__, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&rad, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&red, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&green, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&blue, (ftnlen)sizeof(real)); e_wsfe(); } else if (uij[iatm * 6 - 6] > 0.f) { for (i__ = 1; i__ <= 6; ++i__) { anisou[i__ - 1] = uij[i__ + iatm * 6 - 7]; } if (anitoquad_(anisou, &pradius, quadric, eigens, evecs) < 0) { s_wsle(&io___248); do_lio(&c__9, &c__1, "*** Non-positive definite ellipsoid - ", (ftnlen)38); do_lio(&c__9, &c__1, atom + ((iatm - 1) * 80 + 12), (ftnlen) 14); e_wsle(); ++nerrors; biso[iatm - 1] = 0.f; red = 1.f; green = 0.f; blue = 1.f; } /* Computing MAX */ r__1 = max(eigens[0],eigens[1]); radlim = pradius * dmax(r__1,eigens[2]); radlim *= 1.15f; if (acflag) { /* Computing MIN */ r__1 = min(eigens[0],eigens[1]); /* Computing MAX */ r__2 = max(eigens[0],eigens[1]); anisotropy = dmin(r__1,eigens[2]) / dmax(r__2,eigens[2]); anisotropy *= anisotropy; a2rgb_(&anisotropy, &red, &green, &blue); red *= red; green *= green; blue *= blue; } if (fancy == 5 || fancy == 6) { s_wsfe(&io___249); do_fio(&c__1, (char *)&c__8, (ftnlen)sizeof(integer)); do_fio(&c__1, " -1.0 -1.0 -1.0 -1.0 -1.0 0.0 0 0 0", ( ftnlen)38); i__1 = fancy - 1; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); do_fio(&c__1, "ORTEP_LIKE", (ftnlen)10); e_wsfe(); if (fancy == 6) { s_wsfe(&io___250); do_fio(&c__1, (char *)&c_b566, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&c_b566, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&c_b566, (ftnlen)sizeof(real)); e_wsfe(); } s_wsfe(&io___251); do_fio(&c__1, (char *)&c__0, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&x, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&y, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&z__, (ftnlen)sizeof(real)); for (i__ = 1; i__ <= 3; ++i__) { do_fio(&c__1, (char *)&evecs[i__ - 1], (ftnlen)sizeof( real)); } e_wsfe(); s_wsfe(&io___252); do_fio(&c__1, (char *)&c__0, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&x, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&y, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&z__, (ftnlen)sizeof(real)); for (i__ = 1; i__ <= 3; ++i__) { do_fio(&c__1, (char *)&evecs[i__ + 3], (ftnlen)sizeof( real)); } e_wsfe(); s_wsfe(&io___253); do_fio(&c__1, (char *)&c__0, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&x, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&y, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&z__, (ftnlen)sizeof(real)); for (i__ = 1; i__ <= 3; ++i__) { do_fio(&c__1, (char *)&evecs[i__ + 7], (ftnlen)sizeof( real)); } e_wsfe(); } s_wsfe(&io___254); do_fio(&c__1, (char *)&c__14, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&x, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&y, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&z__, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&radlim, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&red, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&green, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&blue, (ftnlen)sizeof(real)); e_wsfe(); s_wsfe(&io___255); for (i__ = 1; i__ <= 10; ++i__) { do_fio(&c__1, (char *)&quadric[i__ - 1], (ftnlen)sizeof(real)) ; } e_wsfe(); } else { rad = sqrt(spam[iatm * 5 - 2]) * pradius; s_wsfe(&io___256); do_fio(&c__1, (char *)&c__2, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&x, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&y, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&z__, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&rad, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&red, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&green, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&blue, (ftnlen)sizeof(real)); e_wsfe(); } goto L154; } L154: ++iatm; if (iatm <= natm) { goto L150; } if (fancy == 1 || fancy == 3) { s_wsfe(&io___257); do_fio(&c__1, "9 end transparent ellipsoids", (ftnlen)28); e_wsfe(); } if (fancy == 5 || fancy == 6) { s_wsfe(&io___258); do_fio(&c__1, "9 end ortep ellipsoids", (ftnlen)22); e_wsfe(); } L160: /* Write bonds to file also. Atoms are considered bonded if they lie */ /* closer to each other than 0.6 * sum of VDW radii. */ /* If two atoms of different colors are bonded, make half-bond */ /* cylinders with each color. */ if (nerrors == 0) { s_wsle(&io___259); do_lio(&c__9, &c__1, "... no errors found in input file", (ftnlen)33); e_wsle(); } if (radius == 0.f && ! suvflag) { goto L210; } if (suvflag) { io___260.ciunit = suvlun; s_wsfe(&io___260); do_fio(&c__1, "Checking for neighboring atoms with dissimilar Uij ", ( ftnlen)51); do_fio(&c__1, "(Suv < ", (ftnlen)7); do_fio(&c__1, (char *)&suvlimit, (ftnlen)sizeof(real)); do_fio(&c__1, ")...", (ftnlen)4); e_wsfe(); suvbad = 0; } i__1 = natm; for (iatm = 1; iatm <= i__1; ++iatm) { if (nohydro && s_cmp(atom + ((iatm - 1) * 80 + 76), " H", (ftnlen)2, ( ftnlen)2) == 0) { goto L202; } if (suvflag) { if (biso[iatm - 1] == 0.f) { goto L202; } if (uij[iatm * 6 - 6] < 0.f) { goto L202; } } xi = spam[iatm * 5 - 5]; yi = spam[iatm * 5 - 4]; zi = spam[iatm * 5 - 3]; icol = spam[iatm * 5 - 1]; vdwi = vdw[icol - 1]; i__2 = natm; for (jatm = iatm + 1; jatm <= i__2; ++jatm) { close2 = 4.537f; /* Computing 2nd power */ r__1 = xi - spam[jatm * 5 - 5]; dx2 = r__1 * r__1; if (dx2 > close2) { goto L201; } /* Computing 2nd power */ r__1 = yi - spam[jatm * 5 - 4]; dy2 = r__1 * r__1; /* Computing 2nd power */ r__1 = zi - spam[jatm * 5 - 3]; dz2 = r__1 * r__1; dist2 = dx2 + dy2 + dz2; /* Checking for bonded atoms with dissimilar Uij */ if (suvflag) { if (dist2 > close2) { goto L201; } if (biso[jatm - 1] == 0.f) { goto L201; } if (uij[jatm * 6 - 6] < 0.f) { goto L201; } /* DEBUG Version 2.6c used explicit test on 0.6*VdW distance */ /* DEBUG but this is very slow so I have fixed CLOSE2 at the C-S bond length */ /* DEBUG It might be worth doing a first cut using, say, 2.25A (>S-S bond) */ /* DEBUG and then only check the actual VdW distance for pairs making the cut */ /* JCOL = SPAM(5,JATM) */ /* VDWJ = VDW(JCOL) */ /* CLOSE2 = (0.6 * (VDWI + VDWJ)) **2 */ /* IF (DIST2 .GT. CLOSE2) GOTO 201 */ /* DEBUG Add this section back to see if results match V2.6c numbers */ similarity = suv_(&uij[iatm * 6 - 6], &uij[jatm * 6 - 6]); if (similarity < suvlimit) { io___273.ciunit = suvlun; s_wsfe(&io___273); do_fio(&c__1, atom + ((iatm - 1) * 80 + 12), (ftnlen)5); do_fio(&c__1, atom + ((iatm - 1) * 80 + 17), (ftnlen)10); do_fio(&c__1, atom + ((jatm - 1) * 80 + 12), (ftnlen)5); do_fio(&c__1, atom + ((jatm - 1) * 80 + 17), (ftnlen)10); do_fio(&c__1, (char *)&similarity, (ftnlen)sizeof(real)); e_wsfe(); ++suvbad; } if (tflag) { goto L201; } } /* More stringent distance test when drawing bonds */ if (nohydro && s_cmp(atom + ((jatm - 1) * 80 + 76), " H", (ftnlen) 2, (ftnlen)2) == 0) { goto L201; } jcol = spam[jatm * 5 - 1]; vdwj = vdw[jcol - 1]; /* Computing 2nd power */ r__1 = (vdwi + vdwj) * .6f; close2 = r__1 * r__1; if (dist2 > close2) { goto L201; } /* Don't draw bonds between alternate conformers of same residue */ if (*(unsigned char *)&atom[(iatm - 1) * 80 + 16] != ' ' && *( unsigned char *)&atom[(jatm - 1) * 80 + 16] != ' ' && *( unsigned char *)&atom[(iatm - 1) * 80 + 16] != *(unsigned char *)&atom[(jatm - 1) * 80 + 16]) { goto L201; } /* Atoms coloured by B value */ if (bcflag) { s_wsfe(&io___276); do_fio(&c__1, (char *)&spam[iatm * 5 - 5], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[iatm * 5 - 4], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[iatm * 5 - 3], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&radius, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&spam[jatm * 5 - 5], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[jatm * 5 - 4], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[jatm * 5 - 3], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&radius, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&c_b332, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&c_b332, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&c_b332, (ftnlen)sizeof(real)); e_wsfe(); u2rgb_(&spam[iatm * 5 - 2], &umin, &umax, &red1, &green1, & blue1); u2rgb_(&spam[jatm * 5 - 2], &umin, &umax, &red2, &green2, & blue2); s_wsfe(&io___283); do_fio(&c__1, (char *)&red1, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&green1, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&blue1, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&red2, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&green2, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&blue2, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&c_b670, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&c_b670, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&c_b670, (ftnlen)sizeof(real)); e_wsfe(); /* Same color atoms */ } else if (rgb[icol * 3 - 3] == rgb[jcol * 3 - 3] && rgb[icol * 3 - 2] == rgb[jcol * 3 - 2] && rgb[icol * 3 - 1] == rgb[ jcol * 3 - 1]) { s_wsfe(&io___284); do_fio(&c__1, (char *)&spam[iatm * 5 - 5], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[iatm * 5 - 4], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[iatm * 5 - 3], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&radius, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&spam[jatm * 5 - 5], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[jatm * 5 - 4], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[jatm * 5 - 3], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&radius, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&rgb[icol * 3 - 3], (ftnlen)sizeof(real) ); do_fio(&c__1, (char *)&rgb[icol * 3 - 2], (ftnlen)sizeof(real) ); do_fio(&c__1, (char *)&rgb[icol * 3 - 1], (ftnlen)sizeof(real) ); e_wsfe(); } else { for (k = 1; k <= 3; ++k) { center[k - 1] = (spam[k + iatm * 5 - 6] + spam[k + jatm * 5 - 6]) / 2; } s_wsfe(&io___286); do_fio(&c__1, (char *)&spam[iatm * 5 - 5], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[iatm * 5 - 4], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[iatm * 5 - 3], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&radius, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)¢er[0], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)¢er[1], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)¢er[2], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&radius, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&rgb[icol * 3 - 3], (ftnlen)sizeof(real) ); do_fio(&c__1, (char *)&rgb[icol * 3 - 2], (ftnlen)sizeof(real) ); do_fio(&c__1, (char *)&rgb[icol * 3 - 1], (ftnlen)sizeof(real) ); e_wsfe(); s_wsfe(&io___287); do_fio(&c__1, (char *)¢er[0], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)¢er[1], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)¢er[2], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&radius, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&spam[jatm * 5 - 5], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[jatm * 5 - 4], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&spam[jatm * 5 - 3], (ftnlen)sizeof( real)); do_fio(&c__1, (char *)&radius, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&rgb[jcol * 3 - 3], (ftnlen)sizeof(real) ); do_fio(&c__1, (char *)&rgb[jcol * 3 - 2], (ftnlen)sizeof(real) ); do_fio(&c__1, (char *)&rgb[jcol * 3 - 1], (ftnlen)sizeof(real) ); e_wsfe(); } L201: ; } L202: ; } L210: /* 211 FORMAT(1H3,/,11(f8.3)) */ if (suvflag) { if (suvbad == 0) { io___288.ciunit = suvlun; s_wsle(&io___288); do_lio(&c__9, &c__1, "... None!", (ftnlen)9); e_wsle(); } exit_(&c__0); } s_wsfe(&io___289); e_wsfe(); s_wsfe(&io___290); do_fio(&c__1, "X min max center-of-mass:", (ftnlen)26); do_fio(&c__1, (char *)&xmin, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&xmax, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&xcom, (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___291); do_fio(&c__1, "Y min max center-of-mass:", (ftnlen)26); do_fio(&c__1, (char *)&ymin, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&ymax, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&ycom, (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___292); do_fio(&c__1, "Z min max center-of-mass:", (ftnlen)26); do_fio(&c__1, (char *)&zmin, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&zmax, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&zcom, (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___293); do_fio(&c__1, " scale:", (ftnlen)11); do_fio(&c__1, (char *)&scale, (ftnlen)sizeof(real)); e_wsfe(); return 0; } /* MAIN__ */ logical match_(char *subj, char *mask, ftnlen subj_len, ftnlen mask_len) { /* System generated locals */ logical ret_val; /* Local variables */ static integer i__; ret_val = FALSE_; for (i__ = 1; i__ <= 24; ++i__) { if (*(unsigned char *)&subj[i__ - 1] != *(unsigned char *)&mask[i__ - 1] && *(unsigned char *)&mask[i__ - 1] != '#') { return ret_val; } /* L10: */ } ret_val = TRUE_; return ret_val; } /* match_ */ /* Dummy routine to make linker happy (called by QINP in quadric.f) */ /* Subroutine */ int transf_(real *x, real *y, real *z__) { return 0; } /* transf_ */ /* Subroutine */ int assert_(logical *logic, char *dammit, ftnlen dammit_len) { /* Builtin functions */ integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), e_wsle(void); /* Local variables */ extern /* Subroutine */ int exit_(integer *); /* Fortran I/O blocks */ static cilist io___295 = { 0, 0, 0, 0, 0 }; if (*logic) { return 0; } io___295.ciunit = assout_1.noise; s_wsle(&io___295); do_lio(&c__9, &c__1, "ERROR >>>>>> ", (ftnlen)13); do_lio(&c__9, &c__1, dammit, dammit_len); e_wsle(); /* STOP 1234 */ exit_(&c_n1); return 0; } /* assert_ */ /* CC Return RGB triple that runs from dark blue at Bmin */ /* C to light red at Bmax */ /* Subroutine */ int u2rgb_(real *uiso, real *umin, real *umax, real *r__, real *g, real *b) { static real fraction, h__, s, v; extern /* Subroutine */ int hsv2rgb_(real *, real *, real *, real *, real *, real *); fraction = (*uiso - *umin) / (*umax - *umin); if (fraction < 0.f) { fraction = 0.f; } if (fraction > 1.f) { fraction = 1.f; } h__ = (1.f - fraction) * 240.f; s = .95f; v = fraction / 4.f + .75f; hsv2rgb_(&h__, &s, &v, r__, g, b); return 0; } /* u2rgb_ */ /* CC Return RGB triple that runs from */ /* C red for A=0.0 -> white for A=0.5 -> green for A=1.0 */ /* Subroutine */ int a2rgb_(real *a, real *r__, real *g, real *b) { /* System generated locals */ real r__1; /* Local variables */ static real fraction, h__, s, v; extern /* Subroutine */ int hsv2rgb_(real *, real *, real *, real *, real *, real *); if (*a < .5f) { h__ = 0.f; } if (*a >= .5f) { h__ = 120.f; } fraction = (r__1 = (*a - .5f) * 2.f, dabs(r__1)); /* s = fraction**2 */ s = fraction; v = 1.f - s / 2.f; if (*a <= 0.f) { h__ = 300.f; v = 1.f; } hsv2rgb_(&h__, &s, &v, r__, g, b); return 0; } /* a2rgb_ */ /* CC Color format conversion from Hue/Saturation/Value to Red/Green/Blue */ /* C minimal (i.e. NO) error checking */ /* Subroutine */ int hsv2rgb_(real *h__, real *s, real *v, real *r__, real *g, real *b) { static real f; static integer i__; static real p, q, t; i__ = *h__ / 60.f; f = *h__ / 60.f - (real) i__; p = *v * (1.f - *s); q = *v * (1.f - *s * f); t = *v * (1.f - *s * (1.f - f)); if (i__ == 5) { *r__ = *v; *g = p; *b = q; } else if (i__ == 4) { *r__ = t; *g = p; *b = *v; } else if (i__ == 3) { *r__ = p; *g = q; *b = *v; } else if (i__ == 2) { *r__ = p; *g = *v; *b = t; } else if (i__ == 1) { *r__ = q; *g = *v; *b = p; } else { *r__ = *v; *g = t; *b = p; } return 0; } /* hsv2rgb_ */ /* CC This subroutine is not strictly necessary. */ /* C It smooths the histogram/curve of Anisotropy by */ /* distance from center of mass */ /* Subroutine */ int hisclean_(real *value, integer *number, integer *shells, integer *nanis) { /* System generated locals */ integer i__1; /* Local variables */ static integer minshell, maxshell, i__, min10, max10, nsum; static real vsum; /* Parameter adjustments */ --number; --value; /* Function Body */ if (*nanis < 1200) { *shells /= 2; i__1 = *shells; for (i__ = 1; i__ <= i__1; ++i__) { number[i__] = number[(i__ << 1) - 1] + number[i__ * 2]; value[i__] = value[(i__ << 1) - 1] + value[i__ * 2]; } } for (i__ = *shells; i__ >= 1; --i__) { if (number[i__] > 0) { minshell = i__; } } i__1 = *shells; for (i__ = 1; i__ <= i__1; ++i__) { if (number[i__] > 0) { maxshell = i__; } } nsum = 0; min10 = minshell; i__1 = *shells; for (i__ = minshell; i__ <= i__1; ++i__) { if (nsum + number[i__] > 10) { goto L101; } else { nsum += number[i__]; min10 = i__; } } L101: nsum = 0; max10 = maxshell; for (i__ = maxshell; i__ >= 1; --i__) { if (nsum + number[i__] > 10) { goto L102; } else { nsum += number[i__]; max10 = i__; } } L102: nsum = 0; vsum = 0.f; i__1 = min10; for (i__ = minshell; i__ <= i__1; ++i__) { nsum += number[i__]; vsum += value[i__]; number[i__] = 0; value[i__] = 0.f; } i__ = (minshell + min10) / 2; number[i__] = nsum; value[i__] = vsum; nsum = 0; vsum = 0.f; i__1 = max10; for (i__ = maxshell; i__ >= i__1; --i__) { nsum += number[i__]; vsum += value[i__]; number[i__] = 0; value[i__] = 0.f; } i__ = (maxshell + max10) / 2; number[i__] = nsum; value[i__] = vsum; return 0; } /* hisclean_ */ /* Main program alias */ int rastep_ () { MAIN__ (); return 0; }