Examples with Torsion ffx.potential.bonded.Torsion used on opensource projects

Example 11 with Torsion

use of ffx.potential.bonded.Torsion in project ffx by mjschnie.

the class RosenbluthChiAllMove method measureLysine.

public double[] measureLysine(Residue residue, boolean print) {
    if (!residue.getName().contains("LY") || (residue.getAminoAcid3() != AminoAcid3.LYS && residue.getAminoAcid3() != AminoAcid3.LYD)) {
        logger.severe("Yeah that ain't a lysine.");
    }
    double[] chi = new double[4];
    ArrayList<ROLS> torsions = residue.getTorsionList();
    Atom N = (Atom) residue.getAtomNode("N");
    Atom CA = (Atom) residue.getAtomNode("CA");
    Atom CB = (Atom) residue.getAtomNode("CB");
    Atom CD = (Atom) residue.getAtomNode("CD");
    Atom CE = (Atom) residue.getAtomNode("CE");
    Atom CG = (Atom) residue.getAtomNode("CG");
    Atom NZ = (Atom) residue.getAtomNode("NZ");
    logger.info(String.format(" Here's the atoms I found: \n%s\n%s\n%s\n%s\n%s\n%s\n%s", N, CA, CB, CD, CE, CG, NZ));
    logger.info(String.format(" Num torsions: %d", torsions.size()));
    int count = 0;
    for (ROLS rols : torsions) {
        Torsion torsion = (Torsion) rols;
        torsion.energy(false);
        logger.info(String.format(" Torsion numba %d: %s", count++, torsion));
        if (torsion.compare(N, CA, CB, CG)) {
            chi[0] = torsion.getValue();
            if (print) {
                logger.info(torsion.toString());
            }
        }
        if (torsion.compare(CA, CB, CG, CD)) {
            chi[1] = torsion.getValue();
            if (print) {
                logger.info(torsion.toString());
            }
        }
        if (torsion.compare(CB, CG, CD, CE)) {
            chi[2] = torsion.getValue();
            if (print) {
                logger.info(torsion.toString());
            }
        }
        if (torsion.compare(CG, CD, CE, NZ)) {
            chi[3] = torsion.getValue();
            if (print) {
                logger.info(torsion.toString());
            }
        }
    }
    return chi;
}

Example 12 with Torsion

use of ffx.potential.bonded.Torsion in project ffx by mjschnie.

the class RosenbluthOBMC method mcStep.

/**
 * Does all the work for a move. Moveset is a continuous 360 degree spin of
 * the chi[0] torsion. U_or in Frenkel's notation (uDep here) is the
 * associated torsion energy. Evaluation criterion: P_accept = Min( 1,
 * (Wn/Wo)*exp(-beta(U[n]-U[o]) )
 */
private boolean mcStep() {
    numMovesProposed++;
    boolean accepted;
    // Select a target residue.
    int which = ThreadLocalRandom.current().nextInt(targets.size());
    Residue target = targets.get(which);
    ResidueState origState = target.storeState();
    List<ROLS> chiList = target.getTorsionList();
    Torsion chi0 = getChiZeroTorsion(target);
    writeSnapshot("orig");
    /* Create old and new trial sets, calculate Wn and Wo, and choose a move bn.
            When doing strictly chi[0] moves, Frenkel/Smit's 'old' and 'new' configurations
            are the same state. The distinction is made here only to aid in future generalization.
         */
    List<MCMove> oldTrialSet = createTrialSet(target, origState, trialSetSize - 1);
    List<MCMove> newTrialSet = createTrialSet(target, origState, trialSetSize);
    report = new StringBuilder();
    report.append(String.format(" Rosenbluth Rotamer MC Move: %4d\n", numMovesProposed));
    report.append(String.format("    residue:   %s\n", target.toString()));
    report.append(String.format("    chi0:      %s\n", chi0.toString()));
    MCMove proposal = calculateRosenbluthFactors(target, chi0, origState, oldTrialSet, origState, newTrialSet);
    /* Calculate the independent portion of the total old-conf energy.
            Then apply the move and calculate the independent total new-conf energy.
         */
    setState(target, origState);
    writeSnapshot("uIndO");
    double uIndO = getTotalEnergy() - getTorsionEnergy(chi0);
    proposal.move();
    writeSnapshot("uIndN");
    double uIndN = getTotalEnergy() - getTorsionEnergy(chi0);
    // Apply acceptance criterion.
    double temperature = thermostat.getCurrentTemperature();
    double beta = 1.0 / (BOLTZMANN * temperature);
    double dInd = uIndN - uIndO;
    double dIndE = FastMath.exp(-beta * dInd);
    double criterion = (Wn / Wo) * FastMath.exp(-beta * (uIndN - uIndO));
    double metropolis = Math.min(1, criterion);
    double rng = ThreadLocalRandom.current().nextDouble();
    report.append(String.format("    theta:     %3.2f\n", ((RosenbluthChi0Move) proposal).theta));
    report.append(String.format("    criterion: %1.4f\n", criterion));
    report.append(String.format("       Wn/Wo:     %.2f\n", Wn / Wo));
    report.append(String.format("       uIndN,O:  %7.2f\t%7.2f\n", uIndN, uIndO));
    report.append(String.format("       dInd(E):  %7.2f\t%7.2f\n", dInd, dIndE));
    report.append(String.format("    rng:       %1.4f\n", rng));
    if (rng < metropolis) {
        numMovesAccepted++;
        report.append(String.format(" Accepted.\n"));
        accepted = true;
    } else {
        proposal.revertMove();
        report.append(String.format(" Denied.\n"));
        accepted = false;
    }
    logger.info(report.toString());
    // Cleanup.
    Wn = 0.0;
    Wo = 0.0;
    return accepted;
}

Example 13 with Torsion

use of ffx.potential.bonded.Torsion in project ffx by mjschnie.

the class ParticleMeshEwaldQI method applyMaskingRules.

/**
 * This enables PermanentRealSpaceFieldLoop and RealSpaceEnergyLoop to share
 * a code path for masking rule treatment. Defaults: d11scale = 0.0;
 * (scaleD) p12scale = p13scale = 0.0;	(scaleP) m12scale = m13scale = 0.0;
 * (scale) m14scale = 0.4 amoeba, 1/1.2 amber, 0.5 else; m15scale = 0.8
 * amoeba, 1 else.
 */
private void applyMaskingRules(boolean set, int i, double[] maskPerm, double[] maskPolD, double[] maskPolP) {
    final Atom ai = atoms[i];
    for (Atom ak : ai.get1_5s()) {
        if (maskPerm != null) {
            maskPerm[ak.getArrayIndex()] = (set) ? m15scale : 1.0;
        }
    }
    for (Torsion torsion : ai.getTorsions()) {
        Atom ak = torsion.get1_4(ai);
        if (ak != null) {
            final int index = ak.getArrayIndex();
            if (maskPerm != null) {
                maskPerm[index] = (set) ? m14scale : 1.0;
            }
            for (int member : ip11[i]) {
                if (member == index) {
                    maskPolP[index] = (set) ? intra14Scale : 1.0;
                }
            }
        }
    }
    for (Angle angle : ai.getAngles()) {
        Atom ak = angle.get1_3(ai);
        if (ak != null) {
            final int index = ak.getArrayIndex();
            if (maskPerm != null) {
                maskPerm[index] = (set) ? m13scale : 1.0;
            }
            maskPolP[index] = (set) ? p13scale : 1.0;
        }
    }
    for (Bond bond : ai.getBonds()) {
        Atom ak = bond.get1_2(ai);
        final int index = ak.getArrayIndex();
        if (maskPerm != null) {
            maskPerm[index] = (set) ? m12scale : 1.0;
        }
        maskPolP[index] = (set) ? p12scale : 1.0;
    }
    for (int index : ip11[i]) {
        maskPolD[index] = (set) ? d11scale : 1.0;
    }
}

Example 14 with Torsion

use of ffx.potential.bonded.Torsion in project ffx by mjschnie.

the class VanDerWaals method initAtomArrays.

/**
 * Allocate coordinate arrays and set up reduction indices and values.
 */
private void initAtomArrays() {
    if (esvTerm) {
        atoms = esvSystem.getExtendedAtoms();
        nAtoms = atoms.length;
    }
    if (atomClass == null || nAtoms > atomClass.length || lambdaTerm || esvTerm) {
        atomClass = new int[nAtoms];
        coordinates = new double[nAtoms * 3];
        reduced = new double[nSymm][nAtoms * 3];
        reducedXYZ = reduced[0];
        reductionIndex = new int[nAtoms];
        reductionValue = new double[nAtoms];
        bondMask = new int[nAtoms][];
        angleMask = new int[nAtoms][];
        if (vdwForm.vdwType == VanDerWaalsForm.VDW_TYPE.LENNARD_JONES) {
            torsionMask = new int[nAtoms][];
        } else {
            torsionMask = null;
        }
        use = new boolean[nAtoms];
        isSoft = new boolean[nAtoms];
        softCore = new boolean[2][nAtoms];
        lambdaGradX = null;
        lambdaGradY = null;
        lambdaGradZ = null;
        switch(atomicDoubleArrayImpl) {
            case MULTI:
                gradX = new MultiDoubleArray(threadCount, nAtoms);
                gradY = new MultiDoubleArray(threadCount, nAtoms);
                gradZ = new MultiDoubleArray(threadCount, nAtoms);
                if (lambdaTerm) {
                    lambdaGradX = new MultiDoubleArray(threadCount, nAtoms);
                    lambdaGradY = new MultiDoubleArray(threadCount, nAtoms);
                    lambdaGradZ = new MultiDoubleArray(threadCount, nAtoms);
                }
                break;
            case PJ:
                gradX = new PJDoubleArray(threadCount, nAtoms);
                gradY = new PJDoubleArray(threadCount, nAtoms);
                gradZ = new PJDoubleArray(threadCount, nAtoms);
                if (lambdaTerm) {
                    lambdaGradX = new PJDoubleArray(threadCount, nAtoms);
                    lambdaGradY = new PJDoubleArray(threadCount, nAtoms);
                    lambdaGradZ = new PJDoubleArray(threadCount, nAtoms);
                }
                break;
            case ADDER:
            default:
                gradX = new AdderDoubleArray(threadCount, nAtoms);
                gradY = new AdderDoubleArray(threadCount, nAtoms);
                gradZ = new AdderDoubleArray(threadCount, nAtoms);
                if (lambdaTerm) {
                    lambdaGradX = new AdderDoubleArray(threadCount, nAtoms);
                    lambdaGradY = new AdderDoubleArray(threadCount, nAtoms);
                    lambdaGradZ = new AdderDoubleArray(threadCount, nAtoms);
                }
                break;
        }
    }
    /**
     * Initialize all atoms to be used in the energy.
     */
    fill(use, true);
    fill(isSoft, false);
    fill(softCore[HARD], false);
    fill(softCore[SOFT], false);
    softCoreInit = false;
    // Needs initialized regardless of esvTerm.
    esvAtoms = new boolean[nAtoms];
    esvLambda = new double[nAtoms];
    atomEsvID = new int[nAtoms];
    fill(esvAtoms, false);
    fill(esvLambda, 1.0);
    fill(atomEsvID, -1);
    if (esvTerm) {
        updateEsvLambda();
    }
    lambdaFactors = new LambdaFactors[threadCount];
    for (int i = 0; i < threadCount; i++) {
        if (esvTerm) {
            lambdaFactors[i] = new LambdaFactorsESV();
        } else if (lambdaTerm) {
            lambdaFactors[i] = new LambdaFactorsOSRW();
        } else {
            lambdaFactors[i] = new LambdaFactors();
        }
    }
    for (int i = 0; i < nAtoms; i++) {
        Atom ai = atoms[i];
        assert (i == ai.getXyzIndex() - 1);
        double[] xyz = ai.getXYZ(null);
        int i3 = i * 3;
        coordinates[i3 + XX] = xyz[XX];
        coordinates[i3 + YY] = xyz[YY];
        coordinates[i3 + ZZ] = xyz[ZZ];
        AtomType atomType = ai.getAtomType();
        if (atomType == null) {
            logger.severe(ai.toString());
            // Severe no longer guarantees program crash.
            continue;
        }
        String vdwIndex = forceField.getString(ForceField.ForceFieldString.VDWINDEX, "Class");
        if (vdwIndex.equalsIgnoreCase("Type")) {
            atomClass[i] = atomType.type;
        } else {
            atomClass[i] = atomType.atomClass;
        }
        VDWType type = forceField.getVDWType(Integer.toString(atomClass[i]));
        if (type == null) {
            logger.info(" No VdW type for atom class " + atomClass[i]);
            logger.severe(" No VdW type for atom " + ai.toString());
            return;
        }
        ai.setVDWType(type);
        ArrayList<Bond> bonds = ai.getBonds();
        int numBonds = bonds.size();
        if (type.reductionFactor > 0.0 && numBonds == 1) {
            Bond bond = bonds.get(0);
            Atom heavyAtom = bond.get1_2(ai);
            // Atom indexes start at 1
            reductionIndex[i] = heavyAtom.getIndex() - 1;
            reductionValue[i] = type.reductionFactor;
        } else {
            reductionIndex[i] = i;
            reductionValue[i] = 0.0;
        }
        bondMask[i] = new int[numBonds];
        for (int j = 0; j < numBonds; j++) {
            Bond bond = bonds.get(j);
            bondMask[i][j] = bond.get1_2(ai).getIndex() - 1;
        }
        ArrayList<Angle> angles = ai.getAngles();
        int numAngles = 0;
        for (Angle angle : angles) {
            Atom ak = angle.get1_3(ai);
            if (ak != null) {
                numAngles++;
            }
        }
        angleMask[i] = new int[numAngles];
        int j = 0;
        for (Angle angle : angles) {
            Atom ak = angle.get1_3(ai);
            if (ak != null) {
                angleMask[i][j++] = ak.getIndex() - 1;
            }
        }
        if (vdwForm.scale14 != 1.0) {
            ArrayList<Torsion> torsions = ai.getTorsions();
            int numTorsions = 0;
            for (Torsion torsion : torsions) {
                Atom ak = torsion.get1_4(ai);
                if (ak != null) {
                    numTorsions++;
                }
            }
            torsionMask[i] = new int[numTorsions];
            j = 0;
            for (Torsion torsion : torsions) {
                Atom ak = torsion.get1_4(ai);
                if (ak != null) {
                    torsionMask[i][j++] = ak.getIndex() - 1;
                }
            }
        }
    }
}