in layer2/ObjectMap.cpp [4170:4633]
static int ObjectMapGRDStrToMap(ObjectMap * I, char *GRDStr, int bytes, int state,
int quiet)
{
/* NOTE: binary GRD reader not yet validated */
char *p;
float dens;
float *f = NULL;
int a, b, c, d, e;
float v[3], vr[3], maxd, mind;
int ok = true;
char cc[MAXLINELEN];
ObjectMapState *ms;
int got_cell = false;
int got_brick = false;
int fast_axis = 1; /* assumes fast X */
int got_ranges = false;
int normalize = false;
float mean, stdev;
double sum = 0.0;
double sumsq = 0.0;
int n_pts = 0;
int ascii = true;
int little_endian = 1, map_endian = 0;
union int_or_char_array {
int block_len;
char rev[4];
};
union int_or_char_array rev_union;
rev_union.block_len = 0;
if(state < 0)
state = I->NState;
if(I->NState <= state) {
VLACheck(I->State, ObjectMapState, state);
I->NState = state + 1;
}
ms = &I->State[state];
ObjectMapStateInit(I->Obj.G, ms);
normalize = SettingGetGlobal_b(I->Obj.G, cSetting_normalize_grd_maps);
maxd = -FLT_MAX;
mind = FLT_MAX;
p = GRDStr;
if(bytes > 4) {
if((!p[0]) || (!p[1]) || (!p[2]) || (!p[3])) {
/* if zero appears in the first four bytes, then this is a binary map */
ascii = false;
little_endian = *((char *) &little_endian);
map_endian = (*p || *(p + 1));
}
}
/* print map title */
if(ascii) {
p = ParseNCopy(cc, p, 100);
} else {
/*
according to one site on the internet...GRD binary formats look like this:
write(14) title
write(14) scale,oldmid
write(14) ivary,nbyte,intdat,extent,extent
write(14) extent,xang,yang,zang,xstart
write(14) xend,ystart,yend,zstart,zend
write(14) intx,inty,intz
do 100 i = 1,intz+1
do 100 j = 1,inty+1
100 write(14)(phimap(k,j,i),k=1,intx+1)
where: scale is the inverse grid spacing,
oldmid are x,y,z coordinates of the grid center,
intx,inty,intz = grid points per side - 1
ivary = 0
nbyte = 4
intdat = 0
xang,yang,zang = 90
extent is the absolute value of the x,y, or z coordinate of the grid
corners with the largest absolute value
xstart,xend,ystart, yend,zstart,zend are the fractional limits (-1 to 1)
of the molecule in each direction.
title is a descriptive header for the molecule
However, that doesn't make sense...in the files I see, scale and oldmid don't seem to exist!
So here's another claim which I find more credible...
character*132 toplbl !ascii header
integer*4 ivary !0 => x index varys most rapidly
integer*4 nbyte !=4, # of bytes in data
integer*4 inddat !=0, floating point data
real*4 xang,yang,zang !=90,90,90 unit cell angles
integer*4 intx,inty,intz !=igrid-1, # of intervals/grid side
real*4 extent !maximum extent of grid
real*4 xstart,xend !beginning, end of grid sides
real*4 ystart,yend !in fractional
real*4 zstart,zend !units of extent
write(14)toplbl write(14)ivary, nbyte, intdat, extent, extent, extent, xang, yang, zang, xstart, xend, ystart, yend, zstart, zend, intx, inty, intz
do k = 1,igrid
do j = 1,igrid
write(14)(phimap(i,j,k),i=1,igrid)
end do
end d
*/
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Warnings)
" ObjectMapGRD-Warning: Binary GRD reader not yet validated.\n" ENDFB(I->Obj.G);
}
if(little_endian != map_endian) {
rev_union.rev[0] = p[3];
rev_union.rev[1] = p[2];
rev_union.rev[2] = p[1];
rev_union.rev[3] = p[0];
} else {
rev_union.rev[0] = p[0]; /* gotta go char by char because of memory alignment issues ... */
rev_union.rev[1] = p[1];
rev_union.rev[2] = p[2];
rev_union.rev[3] = p[3];
}
ParseNCopy(cc, p + 4, rev_union.block_len);
/* now flip file to correct endianess */
if(little_endian != map_endian) {
char *c = p;
int cnt = bytes >> 2;
unsigned char c0, c1, c2, c3;
while(cnt) {
c0 = c[0];
c1 = c[1];
c2 = c[2];
c3 = c[3];
c[0] = c3;
c[1] = c2;
c[2] = c1;
c[3] = c0;
c += 4;
cnt--;
}
}
p += 4 + rev_union.block_len;
f = (float *) p;
}
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" ObjectMap: %s\n", cc ENDFB(I->Obj.G);
if(ascii)
p = ParseNextLine(p);
/* skip format -- we're reading float regardless... */
if(ascii)
p = ParseNextLine(p);
else {
{
int block_len_check;
int ivary;
int nbyte;
int intdat;
block_len_check = *((int *) (f++));
if(rev_union.block_len != block_len_check) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Warnings)
" ObjectMapGRD-Warning: block length not matched -- not a true GRD binary?\n"
ENDFB(I->Obj.G);
}
rev_union.block_len = *((int *) (f++));
ivary = *((int *) (f++));
nbyte = *((int *) (f++));
intdat = *((int *) (f++));
if((ivary) || (nbyte != 4) || (intdat)) {
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Warnings)
" ObjectMapGRD-Warning: funky ivary, nbyte, intdat -- not a true GRD binary?\n"
ENDFB(I->Obj.G);
}
}
}
}
/* read unit cell */
if(ok) {
if(ascii) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Dim[0]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Dim[1]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Dim[2]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Angle[0]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Angle[1]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%f", &ms->Symmetry->Crystal->Angle[2]) == 1) {
got_cell = true;
}
}
}
}
}
}
ok = got_cell;
} else {
ms->Symmetry->Crystal->Dim[0] = *(f++); /* x-extent */
ms->Symmetry->Crystal->Dim[1] = *(f++); /* y-extent */
ms->Symmetry->Crystal->Dim[2] = *(f++); /* z-extent */
ms->Symmetry->Crystal->Angle[0] = (*f++); /* xang */
ms->Symmetry->Crystal->Angle[1] = (*f++); /* yang */
ms->Symmetry->Crystal->Angle[2] = (*f++); /* zang */
got_cell = 1;
}
}
if(ascii)
p = ParseNextLine(p);
if(ascii) {
/* read brick dimensions */
if(ok) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->FDim[0]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->FDim[1]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->FDim[2]) == 1) {
got_brick = true;
ms->FDim[0]++;
ms->FDim[1]++;
ms->FDim[2]++; /* stupid 0-based counters */
}
}
}
ok = got_brick;
}
p = ParseNextLine(p);
/* read ranges */
if(ok) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &fast_axis) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Min[0]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Div[0]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Min[1]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Div[1]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Min[2]) == 1) {
p = ParseWordCopy(cc, p, 50);
if(sscanf(cc, "%d", &ms->Div[2]) == 1) {
got_ranges = true;
}
}
}
}
}
}
}
ok = got_ranges;
}
} else {
float fmin[3];
float fmax[3];
fmin[0] = *(f++);
fmax[0] = *(f++);
fmin[1] = *(f++);
fmax[1] = *(f++);
fmin[2] = *(f++);
fmax[2] = *(f++);
dump3f(fmin, "fmin");
dump3f(fmax, "fmax");
ms->FDim[0] = *((int *) f++) + 1;
ms->FDim[1] = *((int *) f++) + 1;
ms->FDim[2] = *((int *) f++) + 1;
ms->Div[0] = pymol_roundf((ms->FDim[0] - 1) / (fmax[0] - fmin[0]));
ms->Div[1] = pymol_roundf((ms->FDim[1] - 1) / (fmax[1] - fmin[1]));
ms->Div[2] = pymol_roundf((ms->FDim[2] - 1) / (fmax[2] - fmin[2]));
ms->Min[0] = pymol_roundf(ms->Div[0] * fmin[0]);
ms->Min[1] = pymol_roundf(ms->Div[1] * fmin[1]);
ms->Min[2] = pymol_roundf(ms->Div[2] * fmin[2]);
ms->Max[0] = ms->Min[0] + ms->FDim[0] - 1;
ms->Max[1] = ms->Min[1] + ms->FDim[1] - 1;
ms->Max[2] = ms->Min[2] + ms->FDim[2] - 1;
got_ranges = true;
{
int block_len_check;
block_len_check = *((int *) (f++));
if(rev_union.block_len != block_len_check) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Warnings)
" ObjectMapGRD-Warning: block length not matched -- not a true GRD binary?\n"
ENDFB(I->Obj.G);
}
}
}
if(ok) {
if(ascii) {
ms->Div[0] = ms->Div[0] - ms->Min[0];
ms->Div[1] = ms->Div[1] - ms->Min[1];
ms->Div[2] = ms->Div[2] - ms->Min[2];
ms->Max[0] = ms->Min[0] + ms->FDim[0] - 1;
ms->Max[1] = ms->Min[1] + ms->FDim[1] - 1;
ms->Max[2] = ms->Min[2] + ms->FDim[2] - 1;
}
ms->FDim[3] = 3;
if(Feedback(I->Obj.G, FB_ObjectMap, FB_Blather)) {
dump3i(ms->Div, "ms->Div");
dump3i(ms->Min, "ms->Min");
dump3i(ms->Max, "ms->Max");
dump3i(ms->FDim, "ms->FDim");
}
SymmetryUpdate(ms->Symmetry);
ms->Field = IsosurfFieldAlloc(I->Obj.G, ms->FDim);
ms->MapSource = cMapSourceGRD;
ms->Field->save_points = false;
switch (fast_axis) {
case 3: /* Fast Y - BROKEN! */
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Warnings)
" ObjectMapGRD-Warning: fast_axis %d unsupported!\n", fast_axis ENDFB(I->Obj.G);
/* intentional fall though... */
case 1: /* Fast X */
default:
for(c = 0; c < ms->FDim[2]; c++) {
v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
for(b = 0; b < ms->FDim[1]; b++) {
if(!ascii)
f++; /* skip block delimiter */
v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
for(a = 0; a < ms->FDim[0]; a++) {
v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
if(ascii) {
p = ParseNextLine(p);
p = ParseNCopy(cc, p, 24);
if(sscanf(cc, "%f", &dens) != 1) {
ok = false;
}
} else {
dens = *(f++);
}
if(ok) {
sumsq += dens * dens;
sum += dens;
n_pts++;
F3(ms->Field->data, a, b, c) = dens;
if(maxd < dens)
maxd = dens;
if(mind > dens)
mind = dens;
}
transform33f3f(ms->Symmetry->Crystal->FracToReal, v, vr);
for(e = 0; e < 3; e++) {
F4(ms->Field->points, a, b, c, e) = vr[e];
}
}
if(!ascii)
f++; /* skip fortran block delimiter */
}
}
break;
}
}
if(n_pts > 1) {
mean = (float) (sum / n_pts);
stdev = (float) sqrt1d((sumsq - (sum * sum / n_pts)) / (n_pts - 1));
if(normalize) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" ObjectMapGRDStrToMap: Normalizing: mean = %8.6f and stdev = %8.6f.\n",
mean, stdev ENDFB(I->Obj.G);
if(stdev < R_SMALL8)
stdev = 1.0F;
for(c = 0; c < ms->FDim[2]; c++)
for(b = 0; b < ms->FDim[1]; b++) {
for(a = 0; a < ms->FDim[0]; a++) {
dens = F3(ms->Field->data, a, b, c);
F3(ms->Field->data, a, b, c) = (dens - mean) / stdev;
}
}
} else {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" ObjectMapGRDStrToMap: Mean = %8.6f and stdev = %8.6f.\n",
mean, stdev ENDFB(I->Obj.G);
}
}
if(ok) {
d = 0;
for(c = 0; c < ms->FDim[2]; c += (ms->FDim[2] - 1)) {
v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
for(b = 0; b < ms->FDim[1]; b += (ms->FDim[1] - 1)) {
v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
for(a = 0; a < ms->FDim[0]; a += (ms->FDim[0] - 1)) {
v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
transform33f3f(ms->Symmetry->Crystal->FracToReal, v, vr);
copy3f(vr, ms->Corner + 3 * d);
d++;
}
}
}
v[2] = (ms->Min[2]) / ((float) ms->Div[2]);
v[1] = (ms->Min[1]) / ((float) ms->Div[1]);
v[0] = (ms->Min[0]) / ((float) ms->Div[0]);
transform33f3f(ms->Symmetry->Crystal->FracToReal, v, ms->ExtentMin);
v[2] = ((ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
v[1] = ((ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
v[0] = ((ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
transform33f3f(ms->Symmetry->Crystal->FracToReal, v, ms->ExtentMax);
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Details)
" GRDXStrToMap: Map Size %d x %d x %d\n", ms->FDim[0], ms->FDim[1], ms->FDim[2]
ENDFB(I->Obj.G);
ms->Active = true;
ObjectMapUpdateExtents(I);
if(!quiet) {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Results)
" ObjectMap: Map read. Range: %5.3f to %5.3f\n", mind, maxd ENDFB(I->Obj.G);
}
} else {
PRINTFB(I->Obj.G, FB_ObjectMap, FB_Errors)
" Error: unable to read GRD file.\n" ENDFB(I->Obj.G);
}
return (ok);
}