Examples with HKL ffx.crystal.HKL used on opensource projects

Example 21 with HKL

use of ffx.crystal.HKL in project ffx by mjschnie.

the class CrystalStats method printScaleStats.

/**
 * print scaling and bulk solvent statistics
 */
public void printScaleStats() {
    int[] nhkl = new int[n];
    double[] scale = new double[n];
    for (HKL ih : reflectionlist.hkllist) {
        int i = ih.index();
        int b = ih.bin();
        // ignored cases
        if (Double.isNaN(fo[i][0]) || fo[i][1] <= 0.0) {
            continue;
        }
        // spline setup
        double ss = Crystal.invressq(crystal, ih);
        double fh = spline.f(ss, refinementdata.spline);
        nhkl[b]++;
        scale[b] += (fh - scale[b]) / nhkl[b];
    }
    StringBuilder sb = new StringBuilder(String.format(" Fc to Fo scale: %4.2f\n", exp(0.25 * refinementdata.model_k)));
    sb.append(" Fc to Fo spline scale: ");
    for (int i = 0; i < n; i++) {
        sb.append(String.format("%4.2f ", scale[i]));
    }
    sb.append(String.format("\n Aniso B tensor:\n"));
    sb.append(String.format("  %10.4f %10.4f %10.4f\n", refinementdata.model_b[0], refinementdata.model_b[3], refinementdata.model_b[4]));
    sb.append(String.format("  %10.4f %10.4f %10.4f\n", refinementdata.model_b[3], refinementdata.model_b[1], refinementdata.model_b[5]));
    sb.append(String.format("  %10.4f %10.4f %10.4f\n", refinementdata.model_b[4], refinementdata.model_b[5], refinementdata.model_b[2]));
    if (refinementdata.crs_fs.solventModel != SolventModel.NONE) {
        switch(refinementdata.crs_fs.solventModel) {
            case BINARY:
                sb.append(" Bulk solvent model: Binary mask\n");
                sb.append(String.format("  Probe radius: %8.3f\n  Shrink radius: %8.3f\n", refinementdata.solvent_a, refinementdata.solvent_b));
                break;
            case POLYNOMIAL:
                sb.append(" Bulk solvent model: Polynomial switch\n");
                sb.append(String.format("  a:     %8.3f\n  w:     %8.3f\n", refinementdata.solvent_a, refinementdata.solvent_b));
                break;
            case GAUSSIAN:
                sb.append(" Bulk solvent model: Gaussian\n");
                sb.append(String.format("  A: %8.3f\n  sd scale: %8.3f\n", refinementdata.solvent_a, refinementdata.solvent_b));
                break;
        }
        sb.append(String.format("  Scale: %8.3f\n  B:     %8.3f\n", refinementdata.solvent_k, refinementdata.solvent_ueq * 8.0 * Math.PI * Math.PI));
    }
    sb.append(String.format("\n -log Likelihood: %14.3f (free set: %14.3f)", refinementdata.llkr, refinementdata.llkf));
    if (print) {
        logger.info(sb.toString());
    }
}

Example 22 with HKL

use of ffx.crystal.HKL in project ffx by mjschnie.

the class DiffractionData method writeSolventMaskCNS.

/**
 * write bulk solvent mask for dataset i to a CNS map file
 *
 * @param filename output filename
 * @param i dataset to write out
 */
public void writeSolventMaskCNS(String filename, int i) {
    if (solventModel != SolventModel.NONE) {
        try {
            PrintWriter cnsfile = new PrintWriter(new BufferedWriter(new FileWriter(filename)));
            cnsfile.println(" ANOMalous=FALSE");
            cnsfile.println(" DECLare NAME=FS DOMAin=RECIprocal TYPE=COMP END");
            for (HKL ih : reflectionList[i].hkllist) {
                int j = ih.index();
                cnsfile.printf(" INDE %d %d %d FS= %.4f %.4f\n", ih.h(), ih.k(), ih.l(), refinementData[i].fsF(j), Math.toDegrees(refinementData[i].fsPhi(j)));
            }
            cnsfile.close();
        } catch (Exception e) {
            String message = "Fatal exception evaluating structure factors.\n";
            logger.log(Level.SEVERE, message, e);
            System.exit(-1);
        }
    }
}

Example 23 with HKL

use of ffx.crystal.HKL in project ffx by mjschnie.

the class MTZWriter method write.

/**
 * <p>
 * write</p>
 */
public void write() {
    ByteOrder byteOrder = ByteOrder.nativeOrder();
    FileOutputStream fileOutputStream;
    DataOutputStream dataOutputStream;
    try {
        if (logger.isLoggable(Level.INFO)) {
            StringBuilder sb = new StringBuilder();
            sb.append(format("\n Writing MTZ HKL file: \"%s\"\n", fileName));
            logger.info(sb.toString());
        }
        fileOutputStream = new FileOutputStream(fileName);
        dataOutputStream = new DataOutputStream(fileOutputStream);
        byte[] bytes = new byte[80];
        int offset = 0;
        int writeLen = 0;
        int iMapData;
        float fMapData;
        // Header.
        StringBuilder sb = new StringBuilder();
        sb.append("MTZ ");
        dataOutputStream.writeBytes(sb.toString());
        // Header offset.
        int headerOffset = n * nCol + 21;
        ByteBuffer byteBuffer = ByteBuffer.wrap(bytes);
        byteBuffer.order(byteOrder).putInt(headerOffset);
        // 0x4441 for LE, 0x1111 for BE
        if (ByteOrder.nativeOrder().equals(ByteOrder.LITTLE_ENDIAN)) {
            iMapData = 0x4441;
        } else {
            iMapData = 0x1111;
        }
        byteBuffer.order(byteOrder).putInt(iMapData);
        dataOutputStream.write(bytes, offset, 8);
        sb = new StringBuilder();
        sb.append(" ");
        sb.setLength(68);
        dataOutputStream.writeBytes(sb.toString());
        // Data.
        Vector<String> colname = new Vector<>(nCol);
        char[] colType = new char[nCol];
        double[] res = new double[2];
        res[0] = Double.POSITIVE_INFINITY;
        res[1] = Double.NEGATIVE_INFINITY;
        float[][] colMinMax = new float[nCol][2];
        for (int i = 0; i < nCol; i++) {
            colMinMax[i][0] = Float.POSITIVE_INFINITY;
            colMinMax[i][1] = Float.NEGATIVE_INFINITY;
        }
        ReflectionSpline sigmaASpline = new ReflectionSpline(reflectionList, refinementData.sigmaa.length);
        int col = 0;
        colname.add("H");
        colType[col++] = 'H';
        colname.add("K");
        colType[col++] = 'H';
        colname.add("L");
        colType[col++] = 'H';
        writeLen += 12;
        if (mtzType != MTZType.FCONLY) {
            colname.add("FO");
            colType[col++] = 'F';
            colname.add("SIGFO");
            colType[col++] = 'Q';
            colname.add("FreeR");
            colType[col++] = 'I';
            writeLen += 12;
        }
        if (mtzType != MTZType.DATAONLY) {
            colname.add("Fs");
            colType[col++] = 'F';
            colname.add("PHIFs");
            colType[col++] = 'P';
            colname.add("Fc");
            colType[col++] = 'F';
            colname.add("PHIFc");
            colType[col++] = 'P';
            writeLen += 16;
        }
        if (mtzType == MTZType.ALL) {
            colname.add("FOM");
            colType[col++] = 'W';
            colname.add("PHIW");
            colType[col++] = 'P';
            colname.add("SigmaAs");
            colType[col++] = 'F';
            colname.add("SigmaAw");
            colType[col++] = 'Q';
            colname.add("FWT");
            colType[col++] = 'F';
            colname.add("PHWT");
            colType[col++] = 'P';
            colname.add("DELFWT");
            colType[col++] = 'F';
            colname.add("PHDELWT");
            colType[col++] = 'P';
            writeLen += 32;
        }
        for (HKL ih : reflectionList.hkllist) {
            col = 0;
            int i = ih.index();
            // Skip the 0 0 0 reflection.
            if (ih.h() == 0 && ih.k() == 0 && ih.l() == 0) {
                continue;
            }
            double ss = Crystal.invressq(crystal, ih);
            res[0] = min(ss, res[0]);
            res[1] = max(ss, res[1]);
            // HKL first (3)
            fMapData = ih.h();
            colMinMax[col][0] = min(fMapData, colMinMax[0][0]);
            colMinMax[col][1] = max(fMapData, colMinMax[0][1]);
            byteBuffer.rewind();
            byteBuffer.order(byteOrder).putFloat(fMapData);
            col++;
            fMapData = ih.k();
            colMinMax[col][0] = min(fMapData, colMinMax[1][0]);
            colMinMax[col][1] = max(fMapData, colMinMax[1][1]);
            byteBuffer.order(byteOrder).putFloat(fMapData);
            col++;
            fMapData = ih.l();
            colMinMax[col][0] = min(fMapData, colMinMax[2][0]);
            colMinMax[col][1] = max(fMapData, colMinMax[2][1]);
            byteBuffer.order(byteOrder).putFloat(fMapData);
            col++;
            if (mtzType != MTZType.FCONLY) {
                // F/sigF (2)
                fMapData = (float) refinementData.getF(i);
                if (!isNaN(fMapData)) {
                    colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                    colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                }
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                fMapData = (float) refinementData.getSigF(i);
                if (!isNaN(fMapData)) {
                    colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                    colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                }
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                // Free R (1)
                fMapData = (float) refinementData.getFreeR(i);
                if (!isNaN(fMapData)) {
                    colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                    colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                }
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
            }
            if (mtzType == MTZType.FCONLY) {
                // Fs (2)
                fMapData = (float) refinementData.fsF(i);
                colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                fMapData = (float) toDegrees(refinementData.fsPhi(i));
                colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                // Fc (unscaled!) (2)
                fMapData = (float) refinementData.fcF(i);
                if (!isNaN(fMapData)) {
                    colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                    colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                }
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                fMapData = (float) Math.toDegrees(refinementData.fcPhi(i));
                if (!isNaN(fMapData)) {
                    colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                    colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                }
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
            }
            if (mtzType == MTZType.ALL) {
                // Fs (2)
                fMapData = (float) refinementData.fsF(i);
                colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                fMapData = (float) toDegrees(refinementData.fsPhi(i));
                colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                // Fctot (2)
                fMapData = (float) refinementData.fcTotF(i);
                if (!isNaN(fMapData)) {
                    colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                    colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                }
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                fMapData = (float) toDegrees(refinementData.fcTotPhi(i));
                if (!isNaN(fMapData)) {
                    colMinMax[col][0] = Math.min(fMapData, colMinMax[col][0]);
                    colMinMax[col][1] = Math.max(fMapData, colMinMax[col][1]);
                }
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                // FOM/phase (2)
                fMapData = (float) refinementData.fomphi[i][0];
                if (!isNaN(fMapData)) {
                    colMinMax[col][0] = Math.min(fMapData, colMinMax[col][0]);
                    colMinMax[col][1] = Math.max(fMapData, colMinMax[col][1]);
                }
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                fMapData = (float) toDegrees(refinementData.fomphi[i][1]);
                colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                // Spline setup.
                double fh = spline.f(ss, refinementData.spline);
                double sa = sigmaASpline.f(ss, refinementData.sigmaa);
                double wa = sigmaASpline.f(ss, refinementData.sigmaw);
                // sigmaA/w (2)
                fMapData = (float) sa;
                if (!isNaN(fMapData)) {
                    colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                    colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                }
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                fMapData = (float) wa;
                if (!isNaN(fMapData)) {
                    colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                    colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                }
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                // Map coeffs (4).
                fMapData = (float) refinementData.FoFc2F(i);
                colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                fMapData = (float) toDegrees(refinementData.FoFc2Phi(i));
                colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                fMapData = (float) refinementData.foFc1F(i);
                colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
                fMapData = (float) toDegrees(refinementData.foFc1Phi(i));
                colMinMax[col][0] = min(fMapData, colMinMax[col][0]);
                colMinMax[col][1] = max(fMapData, colMinMax[col][1]);
                byteBuffer.order(byteOrder).putFloat(fMapData);
                col++;
            }
            dataOutputStream.write(bytes, offset, writeLen);
        }
        // Header.
        sb = new StringBuilder();
        sb.append("VERS MTZ:V1.1 ");
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        Date now = new Date();
        SimpleDateFormat sdf = new SimpleDateFormat("yyyy-MM-dd'T'HH:mm:ss ");
        sb = new StringBuilder();
        sb.append("TITLE FFX output: " + sdf.format(now));
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append(String.format("NCOL %8d %12d %8d", nCol, n, 0));
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append("SORT    0    0    0    0    0 ");
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        char cdata = spaceGroup.shortName.charAt(0);
        if (cdata == 'H') {
            cdata = 'R';
        }
        sb.append(String.format("SYMINF %3d %2d %c %5d %22s %5s", spaceGroup.getNumberOfSymOps(), spaceGroup.numPrimitiveSymEquiv, cdata, spaceGroup.number, "'" + spaceGroup.shortName + "'", spaceGroup.pointGroupName));
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        for (int i = 0; i < spaceGroup.symOps.size(); i++) {
            sb = new StringBuilder();
            sb.append("SYMM ");
            SymOp symop = spaceGroup.symOps.get(i);
            sb.append(symop.toXYZString());
            while (sb.length() < 80) {
                sb.append(" ");
            }
            dataOutputStream.writeBytes(sb.toString());
        }
        sb = new StringBuilder();
        sb.append(String.format("RESO %8.6f%13s%8.6f", res[0], " ", res[1]));
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append("VALM NAN ");
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append("NDIF        1 ");
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append("PROJECT       1 project ");
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append("CRYSTAL       1 crystal ");
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append("DATASET       1 dataset ");
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        for (int j = 0; j < nCol; j++) {
            sb = new StringBuilder();
            sb.append(String.format("COLUMN %-30s %c %17.4f %17.4f    1", colname.get(j), colType[j], colMinMax[j][0], colMinMax[j][1]));
            dataOutputStream.writeBytes(sb.toString());
        }
        sb = new StringBuilder();
        sb.append(String.format("CELL %10.4f %9.4f %9.4f %9.4f %9.4f %9.4f ", crystal.a, crystal.b, crystal.c, crystal.alpha, crystal.beta, crystal.gamma));
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append(String.format("DCELL %9d %10.4f %9.4f %9.4f %9.4f %9.4f %9.4f ", 1, crystal.a, crystal.b, crystal.c, crystal.alpha, crystal.beta, crystal.gamma));
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append("DWAVEL        1    1.00000 ");
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append("END ");
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        sb = new StringBuilder();
        sb.append("MTZENDOFHEADERS ");
        while (sb.length() < 80) {
            sb.append(" ");
        }
        dataOutputStream.writeBytes(sb.toString());
        dataOutputStream.close();
    } catch (Exception e) {
        String message = "Fatal exception evaluating structure factors.\n";
        logger.log(Level.SEVERE, message, e);
    }
}

Example 24 with HKL

use of ffx.crystal.HKL in project ffx by mjschnie.

the class MTZFilter method readFile.

/**
 * {@inheritDoc}
 */
@Override
public boolean readFile(File mtzFile, ReflectionList reflectionList, DiffractionRefinementData refinementData, CompositeConfiguration properties) {
    int nRead, nIgnore, nRes, nFriedel, nCut;
    ByteOrder byteOrder = ByteOrder.nativeOrder();
    FileInputStream fileInputStream;
    DataInputStream dataInputStream;
    boolean transpose = false;
    StringBuilder sb = new StringBuilder();
    try {
        fileInputStream = new FileInputStream(mtzFile);
        dataInputStream = new DataInputStream(fileInputStream);
        byte[] headerOffset = new byte[4];
        byte[] bytes = new byte[80];
        int offset = 0;
        // Eat "MTZ" title.
        dataInputStream.read(bytes, offset, 4);
        String mtzstr = null;
        // Header offset.
        dataInputStream.read(headerOffset, offset, 4);
        // Machine stamp.
        dataInputStream.read(bytes, offset, 4);
        ByteBuffer byteBuffer = ByteBuffer.wrap(bytes);
        int stamp = byteBuffer.order(ByteOrder.BIG_ENDIAN).getInt();
        String stampString = Integer.toHexString(stamp);
        switch(stampString.charAt(0)) {
            case '1':
            case '3':
                if (byteOrder.equals(ByteOrder.LITTLE_ENDIAN)) {
                    byteOrder = ByteOrder.BIG_ENDIAN;
                }
                break;
            case '4':
                if (byteOrder.equals(ByteOrder.BIG_ENDIAN)) {
                    byteOrder = ByteOrder.LITTLE_ENDIAN;
                }
                break;
        }
        byteBuffer = ByteBuffer.wrap(headerOffset);
        int headerOffsetI = byteBuffer.order(byteOrder).getInt();
        // skip to header and parse
        dataInputStream.skipBytes((headerOffsetI - 4) * 4);
        for (Boolean parsing = true; parsing; dataInputStream.read(bytes, offset, 80)) {
            mtzstr = new String(bytes);
            parsing = parseHeader(mtzstr);
        }
        // column identifiers
        foString = sigFoString = rFreeString = null;
        if (properties != null) {
            foString = properties.getString("fostring", null);
            sigFoString = properties.getString("sigfostring", null);
            rFreeString = properties.getString("rfreestring", null);
        }
        h = k = l = fo = sigFo = rFree = -1;
        fPlus = sigFPlus = fMinus = sigFMinus = rFreePlus = rFreeMinus = -1;
        boolean print = true;
        parseColumns(print);
        if (h < 0 || k < 0 || l < 0) {
            String message = "Fatal error in MTZ file - no H K L indexes?\n";
            logger.log(Level.SEVERE, message);
            return false;
        }
        // Reopen to start at beginning.
        fileInputStream = new FileInputStream(mtzFile);
        dataInputStream = new DataInputStream(fileInputStream);
        // Skip initial header.
        dataInputStream.skipBytes(80);
        // Check if HKLs need to be transposed or not.
        float[] data = new float[nColumns];
        HKL mate = new HKL();
        int nPosIgnore = 0;
        int nTransIgnore = 0;
        int nZero = 0;
        int none = 0;
        for (int i = 0; i < nReflections; i++) {
            for (int j = 0; j < nColumns; j++) {
                dataInputStream.read(bytes, offset, 4);
                byteBuffer = ByteBuffer.wrap(bytes);
                data[j] = byteBuffer.order(byteOrder).getFloat();
            }
            int ih = (int) data[h];
            int ik = (int) data[k];
            int il = (int) data[l];
            boolean friedel = reflectionList.findSymHKL(ih, ik, il, mate, false);
            HKL hklpos = reflectionList.getHKL(mate);
            if (hklpos == null) {
                nPosIgnore++;
            }
            friedel = reflectionList.findSymHKL(ih, ik, il, mate, true);
            HKL hkltrans = reflectionList.getHKL(mate);
            if (hkltrans == null) {
                nTransIgnore++;
            }
            if (rFree > 0) {
                if (((int) data[rFree]) == 0) {
                    nZero++;
                } else if (((int) data[rFree]) == 1) {
                    none++;
                }
            }
            if (rFreePlus > 0) {
                if (((int) data[rFreePlus]) == 0) {
                    nZero++;
                } else if (((int) data[rFreePlus]) == 1) {
                    none++;
                }
            }
            if (rFreeMinus > 0) {
                if (((int) data[rFreeMinus]) == 0) {
                    nZero++;
                } else if (((int) data[rFreeMinus]) == 1) {
                    none++;
                }
            }
        }
        if (nPosIgnore > nTransIgnore) {
            transpose = true;
        }
        if (none > (nZero * 2) && refinementData.rfreeflag < 0) {
            refinementData.setFreeRFlag(0);
            sb.append(format(" Setting R free flag to %d based on MTZ file data.\n", refinementData.rfreeflag));
        } else if (nZero > (none * 2) && refinementData.rfreeflag < 0) {
            refinementData.setFreeRFlag(1);
            sb.append(format(" Setting R free flag to %d based on MTZ file data.\n", refinementData.rfreeflag));
        } else if (refinementData.rfreeflag < 0) {
            refinementData.setFreeRFlag(0);
            sb.append(format(" Setting R free flag to MTZ default: %d\n", refinementData.rfreeflag));
        }
        // reopen to start at beginning
        fileInputStream = new FileInputStream(mtzFile);
        dataInputStream = new DataInputStream(fileInputStream);
        // skip initial header
        dataInputStream.skipBytes(80);
        // read in data
        double[][] anofSigF = new double[refinementData.n][4];
        for (int i = 0; i < refinementData.n; i++) {
            anofSigF[i][0] = anofSigF[i][1] = anofSigF[i][2] = anofSigF[i][3] = Double.NaN;
        }
        nRead = nIgnore = nRes = nFriedel = nCut = 0;
        for (int i = 0; i < nReflections; i++) {
            for (int j = 0; j < nColumns; j++) {
                dataInputStream.read(bytes, offset, 4);
                byteBuffer = ByteBuffer.wrap(bytes);
                data[j] = byteBuffer.order(byteOrder).getFloat();
            }
            int ih = (int) data[h];
            int ik = (int) data[k];
            int il = (int) data[l];
            boolean friedel = reflectionList.findSymHKL(ih, ik, il, mate, transpose);
            HKL hkl = reflectionList.getHKL(mate);
            if (hkl != null) {
                if (fo > 0 && sigFo > 0) {
                    if (refinementData.fsigfcutoff > 0.0) {
                        if ((data[fo] / data[sigFo]) < refinementData.fsigfcutoff) {
                            nCut++;
                            continue;
                        }
                    }
                    if (friedel) {
                        anofSigF[hkl.index()][2] = data[fo];
                        anofSigF[hkl.index()][3] = data[sigFo];
                        nFriedel++;
                    } else {
                        anofSigF[hkl.index()][0] = data[fo];
                        anofSigF[hkl.index()][1] = data[sigFo];
                    }
                } else {
                    if (fPlus > 0 && sigFPlus > 0) {
                        if (refinementData.fsigfcutoff > 0.0) {
                            if ((data[fPlus] / data[sigFPlus]) < refinementData.fsigfcutoff) {
                                nCut++;
                                continue;
                            }
                        }
                        anofSigF[hkl.index()][0] = data[fPlus];
                        anofSigF[hkl.index()][1] = data[sigFPlus];
                    }
                    if (fMinus > 0 && sigFMinus > 0) {
                        if (refinementData.fsigfcutoff > 0.0) {
                            if ((data[fMinus] / data[sigFMinus]) < refinementData.fsigfcutoff) {
                                nCut++;
                                continue;
                            }
                        }
                        anofSigF[hkl.index()][2] = data[fMinus];
                        anofSigF[hkl.index()][3] = data[sigFMinus];
                    }
                }
                if (rFree > 0) {
                    refinementData.setFreeR(hkl.index(), (int) data[rFree]);
                } else {
                    if (rFreePlus > 0 && rFreeMinus > 0) {
                        // not sure what the correct thing to do here is?
                        refinementData.setFreeR(hkl.index(), (int) data[rFreePlus]);
                    } else if (rFreePlus > 0) {
                        refinementData.setFreeR(hkl.index(), (int) data[rFreePlus]);
                    } else if (rFreeMinus > 0) {
                        refinementData.setFreeR(hkl.index(), (int) data[rFreeMinus]);
                    }
                }
                nRead++;
            } else {
                HKL tmp = new HKL(ih, ik, il);
                if (!reflectionList.resolution.inInverseResSqRange(Crystal.invressq(reflectionList.crystal, tmp))) {
                    nRes++;
                } else {
                    nIgnore++;
                }
            }
        }
        // Set up fsigf from F+ and F-.
        refinementData.generate_fsigf_from_anofsigf(anofSigF);
        // Log results.
        if (logger.isLoggable(Level.INFO)) {
            sb.append(format(" MTZ file type (machine stamp): %s\n", stampString));
            sb.append(format(" HKL data is %s\n", transpose ? "transposed" : "not transposed"));
            sb.append(format(" HKL read in:                             %d\n", nRead));
            sb.append(format(" HKL read as friedel mates:               %d\n", nFriedel));
            sb.append(format(" HKL NOT read in (too high resolution):   %d\n", nRes));
            sb.append(format(" HKL NOT read in (not in internal list?): %d\n", nIgnore));
            sb.append(format(" HKL NOT read in (F/sigF cutoff):         %d\n", nCut));
            sb.append(format(" HKL in internal list:                    %d\n", reflectionList.hkllist.size()));
            logger.info(sb.toString());
        }
        if (rFree < 0 && rFreePlus < 0 && rFreeMinus < 0) {
            refinementData.generateRFree();
        }
    } catch (EOFException e) {
        String message = " MTZ end of file reached.";
        logger.log(Level.WARNING, message, e);
        return false;
    } catch (IOException e) {
        String message = " MTZ IO Exception.";
        logger.log(Level.WARNING, message, e);
        return false;
    }
    return true;
}

Example 25 with HKL

use of ffx.crystal.HKL in project ffx by mjschnie.

the class CrystalReciprocalSpaceTest method test1NSFPermanent.

@Test
public void test1NSFPermanent() {
    String filename = "ffx/xray/structures/1NSF.pdb";
    int index = filename.lastIndexOf(".");
    String name = filename.substring(0, index);
    // load the structure
    ClassLoader cl = this.getClass().getClassLoader();
    File structure = new File(cl.getResource(filename).getPath());
    // load any properties associated with it
    CompositeConfiguration properties = Keyword.loadProperties(structure);
    Crystal crystal = new Crystal(115.996, 115.996, 44.13, 90.0, 90.0, 120.0, "P6");
    Resolution resolution = new Resolution(1.89631);
    ReflectionList reflectionList = new ReflectionList(crystal, resolution);
    DiffractionRefinementData refinementData = new DiffractionRefinementData(properties, reflectionList);
    ForceFieldFilter forceFieldFilter = new ForceFieldFilter(properties);
    ForceField forceField = forceFieldFilter.parse();
    // associate molecular assembly with the structure, set up forcefield
    MolecularAssembly molecularAssembly = new MolecularAssembly(name);
    molecularAssembly.setFile(structure);
    molecularAssembly.setForceField(forceField);
    PDBFilter pdbFile = new PDBFilter(structure, molecularAssembly, forceField, properties);
    pdbFile.readFile();
    pdbFile.applyAtomProperties();
    molecularAssembly.finalize(true, forceField);
    ForceFieldEnergy energy = ForceFieldEnergy.energyFactory(molecularAssembly, pdbFile.getCoordRestraints());
    List<Atom> atomList = molecularAssembly.getAtomList();
    Atom[] atomArray = atomList.toArray(new Atom[atomList.size()]);
    // set up FFT and run it
    ParallelTeam parallelTeam = new ParallelTeam();
    CrystalReciprocalSpace crs = new CrystalReciprocalSpace(reflectionList, atomArray, parallelTeam, parallelTeam);
    crs.computeAtomicDensity(refinementData.fc);
    // tests
    ComplexNumber b = new ComplexNumber(-496.999, 431.817);
    HKL hkl = reflectionList.getHKL(1, 9, 4);
    ComplexNumber a = refinementData.getFc(hkl.index());
    System.out.println("1 9 4: " + a.toString() + " | " + b.toString() + " | " + a.divides(b).toString());
    assertEquals("1 9 4 reflection should be correct", -493.7799429881329, a.re(), 0.0001);
    assertEquals("1 9 4 reflection should be correct", 460.7022632345927, a.im(), 0.0001);
    b.re(-129.767);
    b.im(-76.9812);
    hkl = reflectionList.getHKL(5, 26, 8);
    a = refinementData.getFc(hkl.index());
    System.out.println("5 26 8: " + a.toString() + " | " + b.toString() + " | " + a.divides(b).toString());
    assertEquals("5 26 8 reflection should be correct", -123.05535567943377, a.re(), 0.0001);
    assertEquals("5 26 8 reflection should be correct", -74.59007322382718, a.im(), 0.0001);
}