Search in sources :

Example 1 with ComplexNumber

use of ffx.numerics.ComplexNumber in project ffx by mjschnie.

the class SigmaAMinimize method setWEstimate.

private void setWEstimate() {
    // generate initial w estimate
    ReflectionSpline spline = new ReflectionSpline(reflectionList, refinementData.nbins);
    int[] nmean = new int[refinementData.nbins];
    for (int i = 0; i < refinementData.nbins; i++) {
        nmean[i] = 0;
    }
    double mean = 0.0, tot = 0.0;
    double[][] fc = refinementData.fctot;
    double[][] fo = refinementData.fsigf;
    for (HKL ih : reflectionList.hkllist) {
        int i = ih.index();
        if (ih.allowed() == 0.0 || Double.isNaN(fc[i][0]) || Double.isNaN(fo[i][0])) {
            continue;
        }
        double s2 = Crystal.invressq(crystal, ih);
        double epsc = ih.epsilonc();
        ComplexNumber fct = new ComplexNumber(fc[i][0], fc[i][1]);
        double ecscale = spline.f(s2, refinementData.fcesq);
        double eoscale = spline.f(s2, refinementData.foesq);
        double ec = fct.times(sqrt(ecscale)).abs();
        double eo = fo[i][0] * sqrt(eoscale);
        double wi = pow(eo - ec, 2.0) / epsc;
        nmean[spline.i1()]++;
        tot++;
        x[spline.i1() + refinementData.nbins] += (wi - x[spline.i1() + refinementData.nbins]) / nmean[spline.i1()];
        mean += (wi - mean) / tot;
    }
    logger.info(String.format(" Starting mean w:    %8.3f", mean));
    logger.info(String.format(" Starting w scaling: %8.3f", 1.0 / mean));
    for (int i = 0; i < refinementData.nbins; i++) {
        x[i] -= x[i + refinementData.nbins];
        x[i] *= scaling[i];
        scaling[i + refinementData.nbins] = 1.0 / mean;
        x[i + refinementData.nbins] *= scaling[i + refinementData.nbins];
    }
    sigmaAEnergy.setScaling(scaling);
}
Also used : ReflectionSpline(ffx.crystal.ReflectionSpline) HKL(ffx.crystal.HKL) ComplexNumber(ffx.numerics.ComplexNumber)

Example 2 with ComplexNumber

use of ffx.numerics.ComplexNumber in project ffx by mjschnie.

the class CrystalStats method getRFree.

/**
 * simply return the current Rfree value
 *
 * @return rfree value as a percent
 */
public double getRFree() {
    double sum = 0.0;
    double sumfo = 0.0;
    for (HKL ih : reflectionlist.hkllist) {
        int i = ih.index();
        // ignored cases
        if (Double.isNaN(fc[i][0]) || 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);
        ComplexNumber c = new ComplexNumber(fc[i][0], fc[i][1]);
        if (refinementdata.isFreeR(i)) {
            sum += abs(abs(fo[i][0]) - fh * abs(c.abs()));
            sumfo += abs(fo[i][0]);
        }
    }
    rfree = (sum / sumfo) * 100.0;
    return rfree;
}
Also used : HKL(ffx.crystal.HKL) ComplexNumber(ffx.numerics.ComplexNumber)

Example 3 with ComplexNumber

use of ffx.numerics.ComplexNumber in project ffx by mjschnie.

the class CrystalStats method getR.

/**
 * simply return the current R value
 *
 * @return r value as a percent
 */
public double getR() {
    double numer = 0.0;
    double denom = 0.0;
    double sum = 0.0;
    double sumfo = 0.0;
    double sumall = 0.0;
    double sumfoall = 0.0;
    for (HKL ih : reflectionlist.hkllist) {
        int i = ih.index();
        // ignored cases
        if (Double.isNaN(fc[i][0]) || 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);
        ComplexNumber c = new ComplexNumber(fc[i][0], fc[i][1]);
        numer = abs(abs(fo[i][0]) - fh * abs(c.abs()));
        denom = abs(fo[i][0]);
        sumall += numer;
        sumfoall += denom;
        if (!refinementdata.isFreeR(i)) {
            sum += numer;
            sumfo += denom;
        }
    }
    rall = (sumall / sumfoall) * 100.0;
    r = (sum / sumfo) * 100.0;
    return r;
}
Also used : HKL(ffx.crystal.HKL) ComplexNumber(ffx.numerics.ComplexNumber)

Example 4 with ComplexNumber

use of ffx.numerics.ComplexNumber in project ffx by mjschnie.

the class MTZFilter method readFcs.

/**
 * Read the structure factors.
 *
 * @param mtzFile
 * @param reflectionList
 * @param fcData
 * @param properties
 * @return
 */
public boolean readFcs(File mtzFile, ReflectionList reflectionList, DiffractionRefinementData fcData, CompositeConfiguration properties) {
    int nRead, nIgnore, nRes, nFriedel, nCut;
    ByteOrder byteOrder = ByteOrder.nativeOrder();
    FileInputStream fileInputStream;
    DataInputStream dataInputStream;
    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 mtzString = 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)) {
            mtzString = new String(bytes);
            parsing = parseHeader(mtzString);
        }
        // Column identifiers.
        fc = phiC = fs = phiS = -1;
        boolean print = true;
        parseFcColumns(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);
        float[] data = new float[nColumns];
        HKL mate = new HKL();
        // Read in data.
        ComplexNumber complexNumber = new ComplexNumber();
        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, false);
            HKL hkl = reflectionList.getHKL(mate);
            if (hkl != null) {
                if (fc > 0 && phiC > 0) {
                    complexNumber.re(data[fc] * cos(toRadians(data[phiC])));
                    complexNumber.im(data[fc] * sin(toRadians(data[phiC])));
                    fcData.setFc(hkl.index(), complexNumber);
                }
                if (fs > 0 && phiS > 0) {
                    complexNumber.re(data[fs] * cos(toRadians(data[phiS])));
                    complexNumber.im(data[fs] * sin(toRadians(data[phiS])));
                    fcData.setFs(hkl.index(), complexNumber);
                }
                nRead++;
            } else {
                HKL tmp = new HKL(ih, ik, il);
                if (!reflectionList.resolution.inInverseResSqRange(Crystal.invressq(reflectionList.crystal, tmp))) {
                    nRes++;
                } else {
                    nIgnore++;
                }
            }
        }
        if (logger.isLoggable(Level.INFO)) {
            sb.append(format(" MTZ file type (machine stamp): %s\n", stampString));
            sb.append(format(" Fc HKL read in:                             %d\n", nRead));
            sb.append(format(" Fc HKL read as friedel mates:               %d\n", nFriedel));
            sb.append(format(" Fc HKL NOT read in (too high resolution):   %d\n", nRes));
            sb.append(format(" Fc HKL NOT read in (not in internal list?): %d\n", nIgnore));
            sb.append(format(" Fc 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());
        }
    } 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;
}
Also used : HKL(ffx.crystal.HKL) ByteOrder(java.nio.ByteOrder) ComplexNumber(ffx.numerics.ComplexNumber) IOException(java.io.IOException) DataInputStream(java.io.DataInputStream) ByteBuffer(java.nio.ByteBuffer) FileInputStream(java.io.FileInputStream) EOFException(java.io.EOFException)

Example 5 with ComplexNumber

use of ffx.numerics.ComplexNumber in project ffx by mjschnie.

the class CrystalReciprocalSpaceTest method test1N7SPermanent.

/**
 * Test of permanent method, of class CrystalReciprocalSpace.
 */
@Test
public void test1N7SPermanent() {
    String filename = "ffx/xray/structures/1N7S.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());
    PotentialsUtils potutil = new PotentialsUtils();
    MolecularAssembly mola = potutil.open(structure);
    CompositeConfiguration properties = mola.getProperties();
    Crystal crystal = new Crystal(39.767, 51.750, 132.938, 90.00, 90.00, 90.00, "P212121");
    Resolution resolution = new Resolution(1.45);
    ReflectionList reflectionList = new ReflectionList(crystal, resolution);
    DiffractionRefinementData refinementData = new DiffractionRefinementData(properties, reflectionList);
    mola.finalize(true, mola.getForceField());
    ForceFieldEnergy energy = mola.getPotentialEnergy();
    List<Atom> atomList = mola.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(-828.584, -922.704);
    HKL hkl = reflectionList.getHKL(1, 1, 4);
    ComplexNumber a = refinementData.getFc(hkl.index());
    System.out.println("1 1 4: " + a.toString() + " | " + b.toString() + " | " + a.divides(b).toString());
    assertEquals("1 1 4 reflection should be correct", -753.4722104328416, a.re(), 0.0001);
    assertEquals("1 1 4 reflection should be correct", -1012.1341308707799, a.im(), 0.0001);
    b.re(-70.4582);
    b.im(-486.142);
    hkl = reflectionList.getHKL(2, 1, 10);
    a = refinementData.getFc(hkl.index());
    System.out.println("2 1 10: " + a.toString() + " | " + b.toString() + " | " + a.divides(b).toString());
    assertEquals("2 1 10 reflection should be correct", -69.39660884054359, a.re(), 0.0001);
    assertEquals("2 1 10 reflection should be correct", -412.0147625765328, a.im(), 0.0001);
}
Also used : ParallelTeam(edu.rit.pj.ParallelTeam) CompositeConfiguration(org.apache.commons.configuration.CompositeConfiguration) HKL(ffx.crystal.HKL) ForceFieldEnergy(ffx.potential.ForceFieldEnergy) ReflectionList(ffx.crystal.ReflectionList) ComplexNumber(ffx.numerics.ComplexNumber) Atom(ffx.potential.bonded.Atom) MolecularAssembly(ffx.potential.MolecularAssembly) File(java.io.File) PotentialsUtils(ffx.potential.utils.PotentialsUtils) Crystal(ffx.crystal.Crystal) Resolution(ffx.crystal.Resolution) Test(org.junit.Test)

Aggregations

HKL (ffx.crystal.HKL)8 ComplexNumber (ffx.numerics.ComplexNumber)8 ParallelTeam (edu.rit.pj.ParallelTeam)2 Crystal (ffx.crystal.Crystal)2 ReflectionList (ffx.crystal.ReflectionList)2 ReflectionSpline (ffx.crystal.ReflectionSpline)2 Resolution (ffx.crystal.Resolution)2 ForceFieldEnergy (ffx.potential.ForceFieldEnergy)2 MolecularAssembly (ffx.potential.MolecularAssembly)2 Atom (ffx.potential.bonded.Atom)2 File (java.io.File)2 CompositeConfiguration (org.apache.commons.configuration.CompositeConfiguration)2 Test (org.junit.Test)2 SymOp (ffx.crystal.SymOp)1 ForceField (ffx.potential.parameters.ForceField)1 ForceFieldFilter (ffx.potential.parsers.ForceFieldFilter)1 PDBFilter (ffx.potential.parsers.PDBFilter)1 PotentialsUtils (ffx.potential.utils.PotentialsUtils)1 DataInputStream (java.io.DataInputStream)1 EOFException (java.io.EOFException)1