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

Example 1 with SymOp

use of ffx.crystal.SymOp 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());
    }
}

Example 2 with SymOp

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

the class RotamerOptimization method evaluateDistance.

/**
 * Evaluates the pairwise distance between two residues' rotamers under any
 * symmetry operator; does "lazy loading" for the distance matrix.
 *
 * @param i Residue i
 * @param ri Rotamer for i
 * @param j Residue j
 * @param rj Rotamer for j
 * @return Shortest distance
 */
private double evaluateDistance(int i, int ri, int j, int rj) {
    Residue resi = allResiduesArray[i];
    Rotamer[] rotamersI = resi.getRotamers(library);
    Rotamer roti = rotamersI[ri];
    double[][] xi;
    if (roti.equals(resi.getRotamer())) {
        xi = resi.storeCoordinateArray();
    } else {
        ResidueState origI = resi.storeState();
        RotamerLibrary.applyRotamer(resi, roti);
        xi = resi.storeCoordinateArray();
        resi.revertState(origI);
    }
    Residue resj = allResiduesArray[j];
    Rotamer[] rotamersJ = resj.getRotamers(library);
    Rotamer rotj = rotamersJ[rj];
    double[][] xj;
    if (rotj.equals(resj.getRotamer())) {
        xj = resj.storeCoordinateArray();
    } else {
        ResidueState origJ = resj.storeState();
        RotamerLibrary.applyRotamer(resj, rotj);
        xj = resj.storeCoordinateArray();
        resj.revertState(origJ);
    }
    Crystal crystal = molecularAssembly.getCrystal();
    int nSymm = crystal.spaceGroup.getNumberOfSymOps();
    double minDist = Double.MAX_VALUE;
    for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
        SymOp symOp = crystal.spaceGroup.getSymOp(iSymOp);
        double dist = interResidueDistance(xi, xj, symOp);
        minDist = dist < minDist ? dist : minDist;
    }
    return minDist;
}

Example 3 with SymOp

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

the class NCSRestraint method residual.

public double residual(boolean gradient, boolean print) {
    if (nAtoms % nSymm != 0) {
        return 0;
    }
    if (lambdaTerm) {
        dEdL = 0.0;
        d2EdL2 = 0.0;
        fill(lambdaGradient, 0.0);
    }
    double residual = 0.0;
    int nAsymmAtoms = nAtoms / nSymm;
    double fx2 = forceConstant * 2.0;
    for (int i = 1; i < nSymm; i++) {
        SymOp symOp = spaceGroup.getSymOp(i);
        ncsCrystal.getTransformationOperator(symOp, transOp);
        int offset = nAsymmAtoms * i;
        for (int j = 0; j < nAsymmAtoms; j++) {
            // Reference atom of the asymmetric unit
            Atom atom1 = atoms[j];
            atom1.getXYZ(a1);
            if (atom1.isHydrogen()) {
                continue;
            }
            /**
             * Apply the symmetry operator to superpose the reference atom
             * with its symmetry mate.
             */
            ncsCrystal.applySymOp(a1, a1, symOp);
            // Symmetry mate atom.
            Atom atom2 = atoms[offset + j];
            atom2.getXYZ(a2);
            // Compute their separation.
            dx[0] = a1[0] - a2[0];
            dx[1] = a1[1] - a2[1];
            dx[2] = a1[2] - a2[2];
            // Apply the minimum image convention.
            double r2 = rsq(dx);
            // double r2 = ncsCrystal.image(dx);
            // logger.info(String.format(" %d %16.8f", j, Math.sqrt(r2)));
            residual += r2;
            if (gradient || lambdaTerm) {
                final double dedx = dx[0] * fx2;
                final double dedy = dx[1] * fx2;
                final double dedz = dx[2] * fx2;
                atom2.addToXYZGradient(-lambdaPow * dedx, -lambdaPow * dedy, -lambdaPow * dedz);
                /**
                 * Rotate the force on the reference atom back into the
                 * asymmetric unit.
                 */
                final double dedx1 = dedx * transOp[0][0] + dedy * transOp[1][0] + dedz * transOp[2][0];
                final double dedy1 = dedx * transOp[0][1] + dedy * transOp[1][1] + dedz * transOp[2][1];
                final double dedz1 = dedx * transOp[0][2] + dedy * transOp[1][2] + dedz * transOp[2][2];
                atom1.addToXYZGradient(lambdaPow * dedx1, lambdaPow * dedy1, lambdaPow * dedz1);
                if (lambdaTerm) {
                    int j3 = j * 3;
                    lambdaGradient[j3] += dLambdaPow * dedx1;
                    lambdaGradient[j3 + 1] += dLambdaPow * dedy1;
                    lambdaGradient[j3 + 2] += dLambdaPow * dedz1;
                    int oj3 = (offset + j) * 3;
                    lambdaGradient[oj3] -= dLambdaPow * dedx;
                    lambdaGradient[oj3 + 1] -= dLambdaPow * dedy;
                    lambdaGradient[oj3 + 2] -= dLambdaPow * dedz;
                }
            }
        }
    }
    if (lambdaTerm) {
        dEdL = dLambdaPow * forceConstant * residual;
        d2EdL2 = d2LambdaPow * forceConstant * residual;
    }
    return forceConstant * residual * lambdaPow;
}

Example 4 with SymOp

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

the class NeighborList method getNeighborIndices.

/**
 * Returns the indices of atoms neighboring the passed set of atom indices.
 * Returned Set is exclusive of the passed-in indices.
 *
 * @param atomIndices Atom indices
 * @param maxDist Maximum distance to consider
 * @return Set of neighboring atoms indices
 */
public Set<Integer> getNeighborIndices(List<Integer> atomIndices, double maxDist) {
    double md2 = maxDist * maxDist;
    Set<Integer> neighbors = new HashSet<>();
    for (Integer intI : atomIndices) {
        int i = intI.intValue();
        double[] xyzI = new double[3];
        atoms[i].getXYZ(xyzI);
        for (int iSymm = 0; iSymm < nSymm; iSymm++) {
            SymOp symOp = crystal.spaceGroup.getSymOp(iSymm);
            int[] listi = lists[iSymm][i];
            int nListI = listi.length;
            for (int j = 0; j < nListI; j++) {
                int indexJ = listi[j];
                if (atomIndices.contains(indexJ) || neighbors.contains(indexJ)) {
                    continue;
                }
                double[] xyzJ = new double[3];
                atoms[indexJ].getXYZ(xyzJ);
                crystal.applySymOp(xyzJ, xyzJ, symOp);
                for (int k = 0; k < 3; k++) {
                    xyzJ[k] -= xyzI[k];
                }
                if (crystal.image(xyzJ) <= md2) {
                    neighbors.add(indexJ);
                }
            }
        }
    }
    /*Set<Integer> neighbors = atomIndices.parallelStream().
                mapToInt(Integer::intValue).
                flatMap(i -> {
                    double[] xyzI = new double[3];
                    atoms[i].getXYZ(xyzI);
                    return IntStream.range(0, nSymm).flatMap(iSymm -> {
                        SymOp symOp = crystal.spaceGroup.getSymOp(iSymm);
                        int[] listi = lists[iSymm][i];
                        return IntStream.range(0, listi.length).
                            map(j -> listi[j]).
                            filter(indexJ -> {
                                if (atomIndices.contains(indexJ)) {
                                    return false;
                                }
                                double[] xyzJ = new double[3];
                                atoms[indexJ].getXYZ(xyzJ);
                                crystal.applySymOp(xyzJ, xyzJ, symOp);
                                for (int k = 0; k < 3; k++) {
                                    xyzJ[k] -= xyzI[k];
                                }
                                return crystal.image(xyzJ) <= md2;
                            });
                    }); // Filtering on distinct elements here requires boxing & unboxing.
                }).boxed().collect(Collectors.toSet());*/
    return neighbors;
}

Example 5 with SymOp

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

the class RotamerOptimization method distanceMatrix.

/**
 * Calculates a residue-residue distance matrix.
 * <p>
 * Residue-residue distance is defined as the shortest atom-atom distance in
 * any possible rotamer-rotamer pair if the residues are neighbors (central
 * atom-central atom distances are within a cutoff). Otherewise, distances
 * are set to a default of Double.MAX_VALUE.
 * </p>
 * <p>
 * The intent of using a neighbor list is to avoid tediously searching
 * rotamer- rotamer pairs when two residues are so far apart we will never
 * need the exact distance. We use the distance matrix for adding residues
 * to the sliding window and determining whether to set 3-body energy to
 * 0.0.
 * </p>
 * <p>
 * If the central atoms are too distant from each other, we can safely
 * assume no atoms will ever be close enough for addition to sliding window
 * or to cause explicit calculation of 3-body energy.
 * </p>
 */
private void distanceMatrix() {
    distanceMatrix = new double[numResidues - 1][][][];
    long numDistances = 0L;
    for (int i = 0; i < (numResidues - 1); i++) {
        Residue residuei = allResiduesArray[i];
        int lengthRi;
        try {
            if (checkIfForced(residuei)) {
                lengthRi = 1;
            } else {
                lengthRi = residuei.getRotamers(library).length;
            }
        } catch (IndexOutOfBoundsException ex) {
            if (useForcedResidues) {
                logger.warning(ex.toString());
            } else {
                logIfMaster(format(" Non-forced Residue i %s has null rotamers.", residuei.toFormattedString(false, true)), Level.WARNING);
            }
            continue;
        }
        distanceMatrix[i] = new double[lengthRi][][];
        for (int ri = 0; ri < lengthRi; ri++) {
            distanceMatrix[i][ri] = new double[numResidues][];
            for (int j = (i + 1); j < numResidues; j++) {
                Residue residuej = allResiduesArray[j];
                int lengthRj;
                try {
                    if (checkIfForced(residuej)) {
                        lengthRj = 1;
                    } else {
                        lengthRj = residuej.getRotamers(library).length;
                    }
                } catch (IndexOutOfBoundsException ex) {
                    if (useForcedResidues) {
                        logger.warning(ex.toString());
                    } else {
                        logIfMaster(format(" Residue j %s has null rotamers.", residuej.toFormattedString(false, true)));
                    }
                    continue;
                }
                distanceMatrix[i][ri][j] = new double[lengthRj];
                numDistances += lengthRj;
                if (!lazyMatrix) {
                    fill(distanceMatrix[i][ri][j], Double.MAX_VALUE);
                } else {
                    fill(distanceMatrix[i][ri][j], -1.0);
                }
            }
        }
    }
    logger.info(format(" Number of pairwise distances: %d", numDistances));
    if (!lazyMatrix) {
        ResidueState[] orig = ResidueState.storeAllCoordinates(allResiduesList);
        int nMultiRes = 0;
        /**
         * Build a list that contains one atom from each Residues: CA from
         * amino acids, C1 from nucleic acids, or the first atom otherwise.
         */
        Atom[] atoms = new Atom[numResidues];
        for (int i = 0; i < numResidues; i++) {
            Residue residuei = allResiduesArray[i];
            atoms[i] = residuei.getReferenceAtom();
            if (residuei instanceof MultiResidue) {
                ++nMultiRes;
            }
        }
        /**
         * Use of the pre-existing ParallelTeam causes a conflict when
         * MultiResidues must re-init the force field. Temporary solution
         * for sequence optimization: if > 1 residue optimized, run on only
         * one thread.
         */
        int nThreads = 1;
        if (molecularAssembly.getPotentialEnergy().getParallelTeam() != null) {
            nThreads = (nMultiRes > 1) ? 1 : molecularAssembly.getPotentialEnergy().getParallelTeam().getThreadCount();
        } else {
            // Suggested: nThreads = (nMultiRes > 1) ? 1 : ParallelTeam.getDefaultThreadCount();
            nThreads = 16;
        }
        ParallelTeam parallelTeam = new ParallelTeam(nThreads);
        Crystal crystal = molecularAssembly.getCrystal();
        int nSymm = crystal.spaceGroup.getNumberOfSymOps();
        logger.info("\n Computing Residue Distance Matrix");
        double nlistCutoff = Math.max(Math.max(distance, twoBodyCutoffDist), threeBodyCutoffDist);
        /**
         * I think this originated from the fact that side-chain
         * (and later nucleic acid) atoms could be fairly distant
         * from the reference atom.
         */
        double magicNumberBufferOfUnknownOrigin = 25.0;
        nlistCutoff += magicNumberBufferOfUnknownOrigin;
        NeighborList neighborList = new NeighborList(null, crystal, atoms, nlistCutoff, 0.0, parallelTeam);
        // Expand coordinates
        double[][] xyz = new double[nSymm][3 * numResidues];
        double[] in = new double[3];
        double[] out = new double[3];
        for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
            SymOp symOp = crystal.spaceGroup.getSymOp(iSymOp);
            for (int i = 0; i < numResidues; i++) {
                int i3 = i * 3;
                int iX = i3 + 0;
                int iY = i3 + 1;
                int iZ = i3 + 2;
                Atom atom = atoms[i];
                in[0] = atom.getX();
                in[1] = atom.getY();
                in[2] = atom.getZ();
                crystal.applySymOp(in, out, symOp);
                xyz[iSymOp][iX] = out[0];
                xyz[iSymOp][iY] = out[1];
                xyz[iSymOp][iZ] = out[2];
            }
        }
        // Build the residue neighbor-list.
        int[][][] lists = new int[nSymm][numResidues][];
        boolean[] use = new boolean[numResidues];
        fill(use, true);
        boolean forceRebuild = true;
        boolean printLists = false;
        long neighborTime = -System.nanoTime();
        neighborList.buildList(xyz, lists, use, forceRebuild, printLists);
        neighborTime += System.nanoTime();
        logger.info(format(" Built residue neighbor list:           %8.3f sec", neighborTime * 1.0e-9));
        DistanceRegion distanceRegion = new DistanceRegion(parallelTeam.getThreadCount(), numResidues, crystal, lists, neighborList.getPairwiseSchedule());
        long parallelTime = -System.nanoTime();
        try {
            parallelTeam.execute(distanceRegion);
        } catch (Exception e) {
            String message = " Exception compting residue distance matrix.";
            logger.log(Level.SEVERE, message, e);
        }
        parallelTime += System.nanoTime();
        logger.info(format(" Pairwise distance matrix:              %8.3f sec\n", parallelTime * 1.0e-9));
        ResidueState.revertAllCoordinates(allResiduesList, orig);
        try {
            parallelTeam.shutdown();
        } catch (Exception ex) {
            logger.warning(format(" Exception shutting down parallel team for the distance matrix: %s", ex.toString()));
        }
    }
}