Search in sources :

Example 6 with ComplexNumber

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

the class CrystalReciprocalSpace method computeAtomicGradients.

/**
 * compute inverse FFT to determine atomic gradients
 *
 * @param hklData structure factors to apply inverse FFT
 * @param freer array of free r flags corresponding to hkldata
 * @param flag Rfree flag value
 * @param refinementMode
 * {@link RefinementMinimize.RefinementMode refinement mode}
 * @param print if true, print information on timings during the calculation
 * @see RefinementMinimize.RefinementMode
 * @see DiffractionRefinementData
 */
public void computeAtomicGradients(double[][] hklData, int[] freer, int flag, RefinementMode refinementMode, boolean print) {
    if (solvent && solventModel == SolventModel.NONE) {
        return;
    }
    // zero out the density
    for (int i = 0; i < complexFFT3DSpace; i++) {
        densityGrid[i] = 0.0;
    }
    int nfree = 0;
    StringBuilder sb = new StringBuilder();
    long symtime = -System.nanoTime();
    int nsym = crystal.spaceGroup.symOps.size();
    // int nsym = 1;
    List<SymOp> symops = crystal.spaceGroup.symOps;
    ComplexNumber c = new ComplexNumber();
    ComplexNumber cj = new ComplexNumber();
    HKL ij = new HKL();
    for (HKL ih : reflectionList.hkllist) {
        double[] fc = hklData[ih.index()];
        if (Double.isNaN(fc[0])) {
            continue;
        }
        // cross validation check!!!
        if (freer != null) {
            if (freer[ih.index()] == flag) {
                nfree++;
                continue;
            }
        }
        c.re(fc[0]);
        c.im(fc[1]);
        // scale
        c.timesIP(2.0 / fftScale);
        // apply symmetry
        for (int j = 0; j < nsym; j++) {
            cj.copy(c);
            crystal.applyTransSymRot(ih, ij, symops.get(j));
            double shift = Crystal.sym_phase_shift(ih, symops.get(j));
            int h = Crystal.mod(ij.h(), fftX);
            int k = Crystal.mod(ij.k(), fftY);
            int l = Crystal.mod(ij.l(), fftZ);
            if (h < halfFFTX + 1) {
                final int ii = iComplex3D(h, k, l, fftX, fftY);
                cj.phaseShiftIP(shift);
                densityGrid[ii] += cj.re();
                densityGrid[ii + 1] += -cj.im();
            } else {
                h = (fftX - h) % fftX;
                k = (fftY - k) % fftY;
                l = (fftZ - l) % fftZ;
                final int ii = iComplex3D(h, k, l, fftX, fftY);
                cj.phaseShiftIP(shift);
                densityGrid[ii] += cj.re();
                densityGrid[ii + 1] += cj.im();
            }
        }
    }
    symtime += System.nanoTime();
    long startTime = System.nanoTime();
    complexFFT3D.ifft(densityGrid);
    long fftTime = System.nanoTime() - startTime;
    /*
         CCP4MapWriter mapout = new CCP4MapWriter(fftX, fftY, fftZ, crystal, "/tmp/foo.map");
         mapout.write(densityGrid);
         */
    startTime = System.nanoTime();
    long permanentDensityTime = 0;
    try {
        if (solvent) {
            solventGradientRegion.setRefinementMode(refinementMode);
            parallelTeam.execute(solventGradientRegion);
        } else {
            for (int i = 0; i < nAtoms; i++) {
                optAtomicGradientWeight.set(i, 0);
            }
            atomicGradientSchedule.updateWeights(previousOptAtomicGradientWeight);
            atomicGradientRegion.setRefinementMode(refinementMode);
            parallelTeam.execute(atomicGradientRegion);
            for (int i = 0; i < nAtoms; i++) {
                previousOptAtomicGradientWeight[i] = optAtomicGradientWeight.get(i);
            }
        }
        permanentDensityTime = System.nanoTime() - startTime;
    } catch (Exception e) {
        String message = "Exception computing atomic gradients.";
        logger.log(Level.SEVERE, message, e);
    }
    if (solvent) {
        sb.append(String.format(" Solvent symmetry: %8.3f\n", symtime * toSeconds));
        sb.append(String.format(" Solvent inverse FFT: %8.3f\n", fftTime * toSeconds));
        sb.append(String.format(" Grid solvent gradients: %8.3f\n", permanentDensityTime * toSeconds));
        sb.append(String.format(" %d reflections ignored (cross validation set)\n", nfree));
    } else {
        sb.append(String.format(" Atomic symmetry: %8.3f\n", symtime * toSeconds));
        sb.append(String.format(" Atomic inverse FFT: %8.3f\n", fftTime * toSeconds));
        sb.append(String.format(" Grid atomic gradients: %8.3f\n", permanentDensityTime * toSeconds));
        sb.append(String.format(" %d reflections ignored (cross validation set)\n", nfree));
    }
    if (logger.isLoggable(Level.INFO) && print) {
        logger.info(sb.toString());
    }
}
Also used : SymOp(ffx.crystal.SymOp) HKL(ffx.crystal.HKL) ComplexNumber(ffx.numerics.ComplexNumber)

Example 7 with ComplexNumber

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

the class CrystalStats method printRStats.

/**
 * print R factors and associated statistics in a binned fashion
 */
public void printRStats() {
    double[][] res = new double[n][2];
    double[] nhkl = new double[n + 1];
    double[][] rb = new double[n + 1][2];
    double[][] sumfo = new double[n + 1][2];
    double[][] s = new double[n + 1][4];
    double numer = 0.0;
    double denom = 0.0;
    double sumall = 0.0;
    double sumfoall = 0.0;
    ReflectionSpline sigmaaspline = new ReflectionSpline(reflectionlist, refinementdata.sigmaa.length);
    for (int i = 0; i < n; i++) {
        res[i][0] = Double.NEGATIVE_INFINITY;
        res[i][1] = Double.POSITIVE_INFINITY;
    }
    for (HKL ih : reflectionlist.hkllist) {
        int i = ih.index();
        int b = ih.bin();
        // 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);
        double sa = sigmaaspline.f(ss, refinementdata.sigmaa);
        double wa = sigmaaspline.f(ss, refinementdata.sigmaw);
        double eoscale = sigmaaspline.f(ss, refinementdata.foesq);
        // determine res limits of each bin
        double rs = Crystal.res(crystal, ih);
        if (rs > res[b][0]) {
            res[b][0] = rs;
        }
        if (rs < res[b][1]) {
            res[b][1] = rs;
        }
        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]);
        if (refinementdata.isFreeR(i)) {
            rb[b][1] += numer;
            sumfo[b][1] += denom;
            rb[n][1] += numer;
            sumfo[n][1] += denom;
            sumall += numer;
            sumfoall += denom;
        } else {
            rb[b][0] += numer;
            sumfo[b][0] += denom;
            rb[n][0] += numer;
            sumfo[n][0] += denom;
            sumall += numer;
            sumfoall += denom;
        }
        nhkl[b]++;
        nhkl[n]++;
        s[b][0] += (sa - s[b][0]) / nhkl[b];
        s[b][1] += (wa - s[b][1]) / nhkl[b];
        s[b][2] += ((wa / Math.sqrt(eoscale)) - s[b][2]) / nhkl[b];
        s[b][3] += (fomphi[i][0] - s[b][3]) / nhkl[b];
        s[n][0] += (sa - s[n][0]) / nhkl[n];
        s[n][1] += (wa - s[n][1]) / nhkl[n];
        s[n][2] += ((wa / Math.sqrt(eoscale)) - s[n][2]) / nhkl[n];
        s[n][3] += (fomphi[i][0] - s[n][3]) / nhkl[n];
    }
    StringBuilder sb = new StringBuilder(String.format("\n %15s | %7s | %7s | %7s | %7s | %7s | %7s\n", "Res. Range", "  R", "Rfree", "s", "w(E)", "w(F)", "FOM"));
    for (int i = 0; i < n; i++) {
        sb.append(String.format(" %7.3f %7.3f | ", res[i][0], res[i][1]));
        sb.append(String.format("%7.2f | %7.2f | %7.4f | %7.4f | %7.2f | %7.4f\n", (rb[i][0] / sumfo[i][0]) * 100.0, (rb[i][1] / sumfo[i][1]) * 100.0, s[i][0], s[i][1], s[i][2], s[i][3]));
    }
    sb.append(String.format(" %7.3f %7.3f | ", res[0][0], res[n - 1][1]));
    sb.append(String.format("%7.2f | %7.2f | %7.4f | %7.4f | %7.2f | %7.4f\n", (rb[n][0] / sumfo[n][0]) * 100.0, (rb[n][1] / sumfo[n][1]) * 100.0, s[n][0], s[n][1], s[n][2], s[n][3]));
    sb.append(" s and w are analagous to D and sum_wc");
    reslow = res[0][0];
    reshigh = res[n - 1][1];
    highreslow = res[n - 1][0];
    highreshigh = res[n - 1][1];
    r = (rb[n][0] / sumfo[n][0]) * 100.0;
    rfree = (rb[n][1] / sumfo[n][1]) * 100.0;
    rall = (sumall / sumfoall) * 100.0;
    highr = (rb[n - 1][0] / sumfo[n - 1][0]) * 100.0;
    highrfree = (rb[n - 1][1] / sumfo[n - 1][1]) * 100.0;
    if (print) {
        logger.info(sb.toString());
    }
}
Also used : ReflectionSpline(ffx.crystal.ReflectionSpline) HKL(ffx.crystal.HKL) ComplexNumber(ffx.numerics.ComplexNumber)

Example 8 with ComplexNumber

use of ffx.numerics.ComplexNumber 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);
}
Also used : ParallelTeam(edu.rit.pj.ParallelTeam) CompositeConfiguration(org.apache.commons.configuration.CompositeConfiguration) ForceFieldFilter(ffx.potential.parsers.ForceFieldFilter) 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) ForceField(ffx.potential.parameters.ForceField) File(java.io.File) PDBFilter(ffx.potential.parsers.PDBFilter) 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