Examples with PolarizeType ffx.potential.parameters.PolarizeType used on opensource projects

Example 1 with PolarizeType

use of ffx.potential.parameters.PolarizeType in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method updateAmoebaMultipoleForce.

/**
 * Updates the Amoeba electrostatic multipolar force for change in Use
 * flags.
 *
 * @param atoms Array of all Atoms in the system
 */
private void updateAmoebaMultipoleForce(Atom[] atoms) {
    ParticleMeshEwald pme = super.getPmeNode();
    int[][] axisAtom = pme.getAxisAtoms();
    double dipoleConversion = OpenMM_NmPerAngstrom;
    double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
    double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
    double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
    double polarScale = 1.0;
    if (pme.getPolarizationType() == Polarization.NONE) {
        polarScale = 0.0;
    }
    PointerByReference dipoles = OpenMM_DoubleArray_create(3);
    PointerByReference quadrupoles = OpenMM_DoubleArray_create(9);
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        Atom atom = atoms[i];
        MultipoleType multipoleType = atom.getMultipoleType();
        PolarizeType polarType = atom.getPolarizeType();
        double useFactor = 1.0;
        if (!atoms[i].getUse() || !atoms[i].getElectrostatics()) {
            // if (!atoms[i].getUse()) {
            useFactor = 0.0;
        }
        double lambdaScale = lambda;
        if (!atom.applyLambda()) {
            lambdaScale = 1.0;
        }
        useFactor *= lambdaScale;
        /**
         * Define the frame definition.
         */
        int axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
        switch(multipoleType.frameDefinition) {
            case ZONLY:
                axisType = OpenMM_AmoebaMultipoleForce_ZOnly;
                break;
            case ZTHENX:
                axisType = OpenMM_AmoebaMultipoleForce_ZThenX;
                break;
            case BISECTOR:
                axisType = OpenMM_AmoebaMultipoleForce_Bisector;
                break;
            case ZTHENBISECTOR:
                axisType = OpenMM_AmoebaMultipoleForce_ZBisect;
                break;
            case TRISECTOR:
                axisType = OpenMM_AmoebaMultipoleForce_ThreeFold;
                break;
            default:
                break;
        }
        /**
         * Load local multipole coefficients.
         */
        for (int j = 0; j < 3; j++) {
            OpenMM_DoubleArray_set(dipoles, j, multipoleType.dipole[j] * dipoleConversion * useFactor);
        }
        int l = 0;
        for (int j = 0; j < 3; j++) {
            for (int k = 0; k < 3; k++) {
                OpenMM_DoubleArray_set(quadrupoles, l++, multipoleType.quadrupole[j][k] * quadrupoleConversion / 3.0 * useFactor);
            }
        }
        // int zaxis = 0;
        int zaxis = 1;
        // int xaxis = 0;
        int xaxis = 1;
        // int yaxis = 0;
        int yaxis = 1;
        int[] refAtoms = axisAtom[i];
        if (refAtoms != null) {
            zaxis = refAtoms[0];
            if (refAtoms.length > 1) {
                xaxis = refAtoms[1];
                if (refAtoms.length > 2) {
                    yaxis = refAtoms[2];
                }
            }
        } else {
            axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
        }
        /**
         * Add the multipole.
         */
        OpenMM_AmoebaMultipoleForce_setMultipoleParameters(amoebaMultipoleForce, i, multipoleType.charge * useFactor, dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis, polarType.thole, polarType.pdamp * dampingFactorConversion, polarType.polarizability * polarityConversion * polarScale * useFactor);
    }
    OpenMM_DoubleArray_destroy(dipoles);
    OpenMM_DoubleArray_destroy(quadrupoles);
    OpenMM_AmoebaMultipoleForce_updateParametersInContext(amoebaMultipoleForce, context);
}

Example 2 with PolarizeType

use of ffx.potential.parameters.PolarizeType in project ffx by mjschnie.

the class ParticleMeshEwaldCart method assignMultipole.

private boolean assignMultipole(int i) {
    Atom atom = atoms[i];
    AtomType atomType = atoms[i].getAtomType();
    if (atomType == null) {
        String message = " Multipoles can only be assigned to atoms that have been typed.";
        logger.severe(message);
        return false;
    }
    PolarizeType polarizeType = forceField.getPolarizeType(atomType.getKey());
    if (polarizeType != null) {
        atom.setPolarizeType(polarizeType);
    } else {
        String message = " No polarization type was found for " + atom.toString();
        logger.fine(message);
        double polarizability = 0.0;
        double thole = 0.0;
        int[] polarizationGroup = null;
        polarizeType = new PolarizeType(atomType.type, polarizability, thole, polarizationGroup);
        forceField.addForceFieldType(polarizeType);
        atom.setPolarizeType(polarizeType);
    }
    String key;
    // No reference atoms.
    key = atomType.getKey() + " 0 0";
    MultipoleType multipoleType = forceField.getMultipoleType(key);
    if (multipoleType != null) {
        atom.setMultipoleType(multipoleType);
        localMultipole[i][t000] = multipoleType.getCharge();
        localMultipole[i][t100] = multipoleType.getDipole()[0];
        localMultipole[i][t010] = multipoleType.getDipole()[1];
        localMultipole[i][t001] = multipoleType.getDipole()[2];
        localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
        localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
        localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
        localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
        localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
        localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
        axisAtom[i] = null;
        frame[i] = multipoleType.frameDefinition;
        return true;
    }
    // No bonds.
    List<Bond> bonds = atom.getBonds();
    if (bonds == null || bonds.size() < 1) {
        String message = "Multipoles can only be assigned after bonded relationships are defined.\n";
        logger.severe(message);
    }
    // 1 reference atom.
    for (Bond b : bonds) {
        Atom atom2 = b.get1_2(atom);
        key = atomType.getKey() + " " + atom2.getAtomType().getKey() + " 0";
        multipoleType = multipoleType = forceField.getMultipoleType(key);
        if (multipoleType != null) {
            int[] multipoleReferenceAtoms = new int[1];
            multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
            atom.setMultipoleType(multipoleType);
            localMultipole[i][t000] = multipoleType.getCharge();
            localMultipole[i][t100] = multipoleType.getDipole()[0];
            localMultipole[i][t010] = multipoleType.getDipole()[1];
            localMultipole[i][t001] = multipoleType.getDipole()[2];
            localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
            localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
            localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
            localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
            localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
            localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
            axisAtom[i] = multipoleReferenceAtoms;
            frame[i] = multipoleType.frameDefinition;
            return true;
        }
    }
    // 2 reference atoms.
    for (Bond b : bonds) {
        Atom atom2 = b.get1_2(atom);
        String key2 = atom2.getAtomType().getKey();
        for (Bond b2 : bonds) {
            if (b == b2) {
                continue;
            }
            Atom atom3 = b2.get1_2(atom);
            String key3 = atom3.getAtomType().getKey();
            key = atomType.getKey() + " " + key2 + " " + key3;
            multipoleType = forceField.getMultipoleType(key);
            if (multipoleType != null) {
                int[] multipoleReferenceAtoms = new int[2];
                multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
                multipoleReferenceAtoms[1] = atom3.getIndex() - 1;
                atom.setMultipoleType(multipoleType);
                atom.setMultipoleType(multipoleType);
                localMultipole[i][t000] = multipoleType.getCharge();
                localMultipole[i][t100] = multipoleType.getDipole()[0];
                localMultipole[i][t010] = multipoleType.getDipole()[1];
                localMultipole[i][t001] = multipoleType.getDipole()[2];
                localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
                localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
                localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
                localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
                localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
                localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
                axisAtom[i] = multipoleReferenceAtoms;
                frame[i] = multipoleType.frameDefinition;
                return true;
            }
        }
    }
    /**
     * 3 reference atoms.
     */
    for (Bond b : bonds) {
        Atom atom2 = b.get1_2(atom);
        String key2 = atom2.getAtomType().getKey();
        for (Bond b2 : bonds) {
            if (b == b2) {
                continue;
            }
            Atom atom3 = b2.get1_2(atom);
            String key3 = atom3.getAtomType().getKey();
            for (Bond b3 : bonds) {
                if (b == b3 || b2 == b3) {
                    continue;
                }
                Atom atom4 = b3.get1_2(atom);
                String key4 = atom4.getAtomType().getKey();
                key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
                multipoleType = forceField.getMultipoleType(key);
                if (multipoleType != null) {
                    int[] multipoleReferenceAtoms = new int[3];
                    multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
                    multipoleReferenceAtoms[1] = atom3.getIndex() - 1;
                    multipoleReferenceAtoms[2] = atom4.getIndex() - 1;
                    atom.setMultipoleType(multipoleType);
                    localMultipole[i][t000] = multipoleType.getCharge();
                    localMultipole[i][t100] = multipoleType.getDipole()[0];
                    localMultipole[i][t010] = multipoleType.getDipole()[1];
                    localMultipole[i][t001] = multipoleType.getDipole()[2];
                    localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
                    localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
                    localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
                    localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
                    localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
                    localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
                    axisAtom[i] = multipoleReferenceAtoms;
                    frame[i] = multipoleType.frameDefinition;
                    return true;
                }
            }
            List<Angle> angles = atom.getAngles();
            for (Angle angle : angles) {
                Atom atom4 = angle.get1_3(atom);
                if (atom4 != null) {
                    String key4 = atom4.getAtomType().getKey();
                    key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
                    multipoleType = forceField.getMultipoleType(key);
                    if (multipoleType != null) {
                        int[] multipoleReferenceAtoms = new int[3];
                        multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
                        multipoleReferenceAtoms[1] = atom3.getIndex() - 1;
                        multipoleReferenceAtoms[2] = atom4.getIndex() - 1;
                        atom.setMultipoleType(multipoleType);
                        localMultipole[i][t000] = multipoleType.getCharge();
                        localMultipole[i][t100] = multipoleType.getDipole()[0];
                        localMultipole[i][t010] = multipoleType.getDipole()[1];
                        localMultipole[i][t001] = multipoleType.getDipole()[2];
                        localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
                        localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
                        localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
                        localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
                        localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
                        localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
                        axisAtom[i] = multipoleReferenceAtoms;
                        frame[i] = multipoleType.frameDefinition;
                        return true;
                    }
                }
            }
        }
    }
    /**
     * Revert to a 2 reference atom definition that may include a 1-3 site.
     * For example a hydrogen on water.
     */
    for (Bond b : bonds) {
        Atom atom2 = b.get1_2(atom);
        String key2 = atom2.getAtomType().getKey();
        List<Angle> angles = atom.getAngles();
        for (Angle angle : angles) {
            Atom atom3 = angle.get1_3(atom);
            if (atom3 != null) {
                String key3 = atom3.getAtomType().getKey();
                key = atomType.getKey() + " " + key2 + " " + key3;
                multipoleType = forceField.getMultipoleType(key);
                if (multipoleType != null) {
                    int[] multipoleReferenceAtoms = new int[2];
                    multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
                    multipoleReferenceAtoms[1] = atom3.getIndex() - 1;
                    atom.setMultipoleType(multipoleType);
                    localMultipole[i][t000] = multipoleType.getCharge();
                    localMultipole[i][t100] = multipoleType.getDipole()[0];
                    localMultipole[i][t010] = multipoleType.getDipole()[1];
                    localMultipole[i][t001] = multipoleType.getDipole()[2];
                    localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
                    localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
                    localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
                    localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
                    localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
                    localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
                    axisAtom[i] = multipoleReferenceAtoms;
                    frame[i] = multipoleType.frameDefinition;
                    return true;
                }
                for (Angle angle2 : angles) {
                    Atom atom4 = angle2.get1_3(atom);
                    if (atom4 != null && atom4 != atom3) {
                        String key4 = atom4.getAtomType().getKey();
                        key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
                        multipoleType = forceField.getMultipoleType(key);
                        if (multipoleType != null) {
                            int[] multipoleReferenceAtoms = new int[3];
                            multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
                            multipoleReferenceAtoms[1] = atom3.getIndex() - 1;
                            multipoleReferenceAtoms[2] = atom4.getIndex() - 1;
                            atom.setMultipoleType(multipoleType);
                            localMultipole[i][t000] = multipoleType.getCharge();
                            localMultipole[i][t100] = multipoleType.getDipole()[0];
                            localMultipole[i][t010] = multipoleType.getDipole()[1];
                            localMultipole[i][t001] = multipoleType.getDipole()[2];
                            localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
                            localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
                            localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
                            localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
                            localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
                            localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
                            axisAtom[i] = multipoleReferenceAtoms;
                            frame[i] = multipoleType.frameDefinition;
                            return true;
                        }
                    }
                }
            }
        }
    }
    return false;
}

Example 3 with PolarizeType

use of ffx.potential.parameters.PolarizeType in project ffx by mjschnie.

the class ParticleMeshEwaldCart method initAtomArrays.

private void initAtomArrays() {
    if (localMultipole == null || localMultipole.length < nAtoms) {
        localMultipole = new double[nAtoms][10];
        frame = new MultipoleType.MultipoleFrameDefinition[nAtoms];
        axisAtom = new int[nAtoms][];
        cartMultipolePhi = new double[nAtoms][tensorCount];
        directDipole = new double[nAtoms][3];
        directDipoleCR = new double[nAtoms][3];
        cartesianDipolePhi = new double[nAtoms][tensorCount];
        cartesianDipolePhiCR = new double[nAtoms][tensorCount];
        ip11 = new int[nAtoms][];
        ip12 = new int[nAtoms][];
        ip13 = new int[nAtoms][];
        thole = new double[nAtoms];
        ipdamp = new double[nAtoms];
        polarizability = new double[nAtoms];
        realSpaceSchedule = new PairwiseSchedule(maxThreads, nAtoms, realSpaceRanges);
        if (scfAlgorithm == SCFAlgorithm.CG) {
            rsd = new double[3][nAtoms];
            rsdCR = new double[3][nAtoms];
            rsdPre = new double[3][nAtoms];
            rsdPreCR = new double[3][nAtoms];
            conj = new double[3][nAtoms];
            conjCR = new double[3][nAtoms];
            vec = new double[3][nAtoms];
            vecCR = new double[3][nAtoms];
        }
        if (scfPredictor != SCFPredictor.NONE) {
            if (lambdaTerm) {
                predictorInducedDipole = new double[3][predictorOrder][nAtoms][3];
                predictorInducedDipoleCR = new double[3][predictorOrder][nAtoms][3];
            } else {
                predictorInducedDipole = new double[1][predictorOrder][nAtoms][3];
                predictorInducedDipoleCR = new double[1][predictorOrder][nAtoms][3];
            }
        }
        /**
         * Initialize per-thread memory for collecting the gradient, torque,
         * field and chain-rule field.
         */
        grad = new double[maxThreads][3][nAtoms];
        torque = new double[maxThreads][3][nAtoms];
        field = new double[maxThreads][3][nAtoms];
        fieldCR = new double[maxThreads][3][nAtoms];
        if (lambdaTerm) {
            lambdaGrad = new double[maxThreads][3][nAtoms];
            lambdaTorque = new double[maxThreads][3][nAtoms];
        }
        isSoft = new boolean[nAtoms];
        use = new boolean[nAtoms];
        coordinates = new double[nSymm][3][nAtoms];
        globalMultipole = new double[nSymm][nAtoms][10];
        inducedDipole = new double[nSymm][nAtoms][3];
        inducedDipoleCR = new double[nSymm][nAtoms][3];
        /**
         * The size of reduced neighbor list depends on the size of the real
         * space cutoff.
         */
        realSpaceLists = new int[nSymm][nAtoms][];
        realSpaceCounts = new int[nSymm][nAtoms];
        preconditionerLists = new int[nSymm][nAtoms][preconditionerListSize];
        preconditionerCounts = new int[nSymm][nAtoms];
    }
    /**
     * Initialize the soft core lambda mask to false for all atoms.
     */
    fill(isSoft, false);
    /**
     * Initialize the use mask to true for all atoms.
     */
    fill(use, true);
    /**
     * Assign multipole parameters.
     */
    assignMultipoles();
    /**
     * Assign polarization groups.
     */
    assignPolarizationGroups();
    /**
     * Fill the thole, inverse polarization damping and polarizability
     * arrays.
     */
    for (Atom ai : atoms) {
        PolarizeType polarizeType = ai.getPolarizeType();
        int index = ai.getIndex() - 1;
        thole[index] = polarizeType.thole;
        ipdamp[index] = polarizeType.pdamp;
        if (!(ipdamp[index] > 0.0)) {
            ipdamp[index] = Double.POSITIVE_INFINITY;
        } else {
            ipdamp[index] = 1.0 / ipdamp[index];
        }
        polarizability[index] = polarizeType.polarizability;
    }
}

Example 4 with PolarizeType

use of ffx.potential.parameters.PolarizeType in project ffx by mjschnie.

the class ParticleMeshEwaldQI method initAtomArrays.

private void initAtomArrays() {
    if (localMultipole == null || localMultipole.length < nAtoms || lambdaTerm || esvTerm) {
        localMultipole = new double[nAtoms][10];
        frame = new MultipoleType.MultipoleFrameDefinition[nAtoms];
        axisAtom = new int[nAtoms][];
        cartMultipolePhi = new double[nAtoms][tensorCount];
        directDipole = new double[nAtoms][3];
        directDipoleCR = new double[nAtoms][3];
        cartesianDipolePhi = new double[nAtoms][tensorCount];
        cartesianDipolePhiCR = new double[nAtoms][tensorCount];
        ip11 = new int[nAtoms][];
        ip12 = new int[nAtoms][];
        ip13 = new int[nAtoms][];
        thole = new double[nAtoms];
        ipdamp = new double[nAtoms];
        polarizability = new double[nAtoms];
        realSpaceSchedule = new PairwiseSchedule(maxThreads, nAtoms, realSpaceRanges);
        if (scfAlgorithm == SCFAlgorithm.CG) {
            rsd = new double[3][nAtoms];
            rsdCR = new double[3][nAtoms];
            rsdPre = new double[3][nAtoms];
            rsdPreCR = new double[3][nAtoms];
            conj = new double[3][nAtoms];
            conjCR = new double[3][nAtoms];
            vec = new double[3][nAtoms];
            vecCR = new double[3][nAtoms];
        }
        /**
         * Initialize per-thread memory for collecting the gradient, torque,
         * field and chain-rule field.
         */
        grad = new double[maxThreads][3][nAtoms];
        torque = new double[maxThreads][3][nAtoms];
        field = new double[maxThreads][3][nAtoms];
        fieldCR = new double[maxThreads][3][nAtoms];
        if (lambdaTerm) {
            lambdaGrad = new double[maxThreads][3][nAtoms];
            lambdaTorque = new double[maxThreads][3][nAtoms];
        }
        isSoft = new boolean[nAtoms];
        fill(isSoft, false);
        use = new boolean[nAtoms];
        fill(use, true);
        coordinates = new double[nSymm][3][nAtoms];
        globalMultipole = new double[nSymm][nAtoms][10];
        inducedDipole = new double[nSymm][nAtoms][3];
        inducedDipoleCR = new double[nSymm][nAtoms][3];
        if (scfPredictor != null) {
            if (esvTerm) {
                throw new UnsupportedOperationException();
            }
            scfPredictor.setInducedDipoleReferences(inducedDipole, inducedDipoleCR, lambdaTerm);
        }
        /* ESV flag array, initialized regardless of esvTerm. */
        // True for other ESV residue atoms.
        esvAtomsScaled = new boolean[nAtoms];
        esvAtomsScaledAlpha = new boolean[nAtoms];
        fill(esvAtomsScaled, false);
        fill(esvAtomsScaledAlpha, false);
        lambdaFactors = new LambdaFactors[maxThreads];
        for (int i = 0; i < maxThreads; i++) {
            if (esvTerm) {
                // Invoked every time through inner loops.
                lambdaFactors[i] = new LambdaFactorsESV();
            } else if (lambdaTerm) {
                // Invoked on calls to setLambda().
                lambdaFactors[i] = new LambdaFactorsOSRW();
            } else {
                // Invoked never; inoperative defaults.
                lambdaFactors[i] = LambdaDefaults;
            }
        }
        /**
         * The size of reduced neighbor list depends on the size of the real
         * space cutoff.
         */
        realSpaceLists = new int[nSymm][nAtoms][];
        realSpaceCounts = new int[nSymm][nAtoms];
        preconditionerLists = new int[nSymm][nAtoms][preconditionerListSize];
        preconditionerCounts = new int[nSymm][nAtoms];
    }
    /**
     * Assign multipole parameters and polarization groups.
     */
    for (int i = 0; i < nAtoms; i++) {
        MultipoleType.multipoleTypeFactory(atoms[i], forceField);
        localMultipole[i] = atoms[i].getMultipoleType().getMultipole();
        axisAtom[i] = atoms[i].getAxisAtomIndices();
        frame[i] = atoms[i].getMultipoleType().frameDefinition;
    }
    /**
     * Assign polarization groups.
     */
    assignPolarizationGroups();
    /**
     * Fill the thole, inverse polarization damping and polarizability
     * arrays.
     */
    for (Atom ai : atoms) {
        PolarizeType polarizeType = ai.getPolarizeType();
        int index = ai.getIndex() - 1;
        thole[index] = polarizeType.thole;
        ipdamp[index] = polarizeType.pdamp;
        if (!(ipdamp[index] > 0.0)) {
            ipdamp[index] = Double.POSITIVE_INFINITY;
        } else {
            ipdamp[index] = 1.0 / ipdamp[index];
        }
        polarizability[index] = polarizeType.polarizability;
    }
    if (esvTerm) {
        updateEsvLambda();
    }
}

Example 5 with PolarizeType

use of ffx.potential.parameters.PolarizeType in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method addAmoebaMultipoleForce.

private void addAmoebaMultipoleForce() {
    ParticleMeshEwald pme = super.getPmeNode();
    if (pme == null) {
        return;
    }
    int[][] axisAtom = pme.getAxisAtoms();
    double dipoleConversion = OpenMM_NmPerAngstrom;
    double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
    double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
    double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
    amoebaMultipoleForce = OpenMM_AmoebaMultipoleForce_create();
    double polarScale = 1.0;
    if (pme.getPolarizationType() != Polarization.MUTUAL) {
        OpenMM_AmoebaMultipoleForce_setPolarizationType(amoebaMultipoleForce, OpenMM_AmoebaMultipoleForce_Direct);
        if (pme.getPolarizationType() == Polarization.NONE) {
            polarScale = 0.0;
        }
    } else {
        ForceField forceField = molecularAssembly.getForceField();
        String algorithm = forceField.getString(ForceField.ForceFieldString.SCF_ALGORITHM, "CG");
        ParticleMeshEwald.SCFAlgorithm scfAlgorithm;
        try {
            algorithm = algorithm.replaceAll("-", "_").toUpperCase();
            scfAlgorithm = ParticleMeshEwald.SCFAlgorithm.valueOf(algorithm);
        } catch (Exception e) {
            scfAlgorithm = ParticleMeshEwald.SCFAlgorithm.CG;
        }
        switch(scfAlgorithm) {
            case EPT:
                logger.info(" Using extrapolated perturbation theory approximation instead of full SCF calculations. Not supported in FFX reference implementation.");
                OpenMM_AmoebaMultipoleForce_setPolarizationType(amoebaMultipoleForce, OpenMM_AmoebaMultipoleForce_Extrapolated);
                PointerByReference exptCoefficients = OpenMM_DoubleArray_create(4);
                OpenMM_DoubleArray_set(exptCoefficients, 0, -0.154);
                OpenMM_DoubleArray_set(exptCoefficients, 1, 0.017);
                OpenMM_DoubleArray_set(exptCoefficients, 2, 0.657);
                OpenMM_DoubleArray_set(exptCoefficients, 3, 0.475);
                OpenMM_AmoebaMultipoleForce_setExtrapolationCoefficients(amoebaMultipoleForce, exptCoefficients);
                OpenMM_DoubleArray_destroy(exptCoefficients);
                break;
            case CG:
            case SOR:
            default:
                OpenMM_AmoebaMultipoleForce_setPolarizationType(amoebaMultipoleForce, OpenMM_AmoebaMultipoleForce_Mutual);
                break;
        }
    }
    PointerByReference dipoles = OpenMM_DoubleArray_create(3);
    PointerByReference quadrupoles = OpenMM_DoubleArray_create(9);
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        Atom atom = atoms[i];
        MultipoleType multipoleType = atom.getMultipoleType();
        PolarizeType polarType = atom.getPolarizeType();
        /**
         * Define the frame definition.
         */
        int axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
        switch(multipoleType.frameDefinition) {
            case ZONLY:
                axisType = OpenMM_AmoebaMultipoleForce_ZOnly;
                break;
            case ZTHENX:
                axisType = OpenMM_AmoebaMultipoleForce_ZThenX;
                break;
            case BISECTOR:
                axisType = OpenMM_AmoebaMultipoleForce_Bisector;
                break;
            case ZTHENBISECTOR:
                axisType = OpenMM_AmoebaMultipoleForce_ZBisect;
                break;
            case TRISECTOR:
                axisType = OpenMM_AmoebaMultipoleForce_ThreeFold;
                break;
            default:
                break;
        }
        double useFactor = 1.0;
        if (!atoms[i].getUse() || !atoms[i].getElectrostatics()) {
            // if (!atoms[i].getUse()) {
            useFactor = 0.0;
        }
        // Should be 1.0 at this point.
        double lambdaScale = lambda;
        if (!atom.applyLambda()) {
            lambdaScale = 1.0;
        }
        useFactor *= lambdaScale;
        /**
         * Load local multipole coefficients.
         */
        for (int j = 0; j < 3; j++) {
            OpenMM_DoubleArray_set(dipoles, j, multipoleType.dipole[j] * dipoleConversion * useFactor);
        }
        int l = 0;
        for (int j = 0; j < 3; j++) {
            for (int k = 0; k < 3; k++) {
                OpenMM_DoubleArray_set(quadrupoles, l++, multipoleType.quadrupole[j][k] * quadrupoleConversion * useFactor / 3.0);
            }
        }
        // int zaxis = 0;
        int zaxis = 1;
        // int xaxis = 0;
        int xaxis = 1;
        // int yaxis = 0;
        int yaxis = 1;
        int[] refAtoms = axisAtom[i];
        if (refAtoms != null) {
            zaxis = refAtoms[0];
            if (refAtoms.length > 1) {
                xaxis = refAtoms[1];
                if (refAtoms.length > 2) {
                    yaxis = refAtoms[2];
                }
            }
        } else {
            axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
            logger.info(String.format(" Atom type %s", atom.getAtomType().toString()));
        }
        double charge = multipoleType.charge * useFactor;
        /**
         * Add the multipole.
         */
        OpenMM_AmoebaMultipoleForce_addMultipole(amoebaMultipoleForce, charge, dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis, polarType.thole, polarType.pdamp * dampingFactorConversion, polarType.polarizability * polarityConversion * polarScale);
    }
    OpenMM_DoubleArray_destroy(dipoles);
    OpenMM_DoubleArray_destroy(quadrupoles);
    Crystal crystal = super.getCrystal();
    if (!crystal.aperiodic()) {
        OpenMM_AmoebaMultipoleForce_setNonbondedMethod(amoebaMultipoleForce, OpenMM_AmoebaMultipoleForce_PME);
        OpenMM_AmoebaMultipoleForce_setCutoffDistance(amoebaMultipoleForce, pme.getEwaldCutoff() * OpenMM_NmPerAngstrom);
        OpenMM_AmoebaMultipoleForce_setAEwald(amoebaMultipoleForce, pme.getEwaldCoefficient() / OpenMM_NmPerAngstrom);
        double ewaldTolerance = 1.0e-04;
        OpenMM_AmoebaMultipoleForce_setEwaldErrorTolerance(amoebaMultipoleForce, ewaldTolerance);
        PointerByReference gridDimensions = OpenMM_IntArray_create(3);
        ReciprocalSpace recip = pme.getReciprocalSpace();
        OpenMM_IntArray_set(gridDimensions, 0, recip.getXDim());
        OpenMM_IntArray_set(gridDimensions, 1, recip.getYDim());
        OpenMM_IntArray_set(gridDimensions, 2, recip.getZDim());
        OpenMM_AmoebaMultipoleForce_setPmeGridDimensions(amoebaMultipoleForce, gridDimensions);
        OpenMM_IntArray_destroy(gridDimensions);
    } else {
        OpenMM_AmoebaMultipoleForce_setNonbondedMethod(amoebaMultipoleForce, OpenMM_AmoebaMultipoleForce_NoCutoff);
    }
    OpenMM_AmoebaMultipoleForce_setMutualInducedMaxIterations(amoebaMultipoleForce, 500);
    OpenMM_AmoebaMultipoleForce_setMutualInducedTargetEpsilon(amoebaMultipoleForce, pme.getPolarEps());
    int[][] ip11 = pme.getPolarization11();
    int[][] ip12 = pme.getPolarization12();
    int[][] ip13 = pme.getPolarization13();
    ArrayList<Integer> list12 = new ArrayList<>();
    ArrayList<Integer> list13 = new ArrayList<>();
    ArrayList<Integer> list14 = new ArrayList<>();
    PointerByReference covalentMap = OpenMM_IntArray_create(0);
    for (int i = 0; i < nAtoms; i++) {
        Atom ai = atoms[i];
        list12.clear();
        list13.clear();
        list14.clear();
        for (Bond bond : ai.getBonds()) {
            int index = bond.get1_2(ai).getIndex() - 1;
            OpenMM_IntArray_append(covalentMap, index);
            list12.add(index);
        }
        OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i, OpenMM_AmoebaMultipoleForce_Covalent12, covalentMap);
        OpenMM_IntArray_resize(covalentMap, 0);
        for (Angle angle : ai.getAngles()) {
            Atom ak = angle.get1_3(ai);
            if (ak != null) {
                int index = ak.getIndex() - 1;
                if (!list12.contains(index)) {
                    list13.add(index);
                    OpenMM_IntArray_append(covalentMap, index);
                }
            }
        }
        OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i, OpenMM_AmoebaMultipoleForce_Covalent13, covalentMap);
        OpenMM_IntArray_resize(covalentMap, 0);
        for (Torsion torsion : ai.getTorsions()) {
            Atom ak = torsion.get1_4(ai);
            if (ak != null) {
                int index = ak.getIndex() - 1;
                if (!list12.contains(index) && !list13.contains(index)) {
                    list14.add(index);
                    OpenMM_IntArray_append(covalentMap, index);
                }
            }
        }
        OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i, OpenMM_AmoebaMultipoleForce_Covalent14, covalentMap);
        OpenMM_IntArray_resize(covalentMap, 0);
        for (Atom ak : ai.get1_5s()) {
            int index = ak.getIndex() - 1;
            if (!list12.contains(index) && !list13.contains(index) && !list14.contains(index)) {
                OpenMM_IntArray_append(covalentMap, index);
            }
        }
        OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i, OpenMM_AmoebaMultipoleForce_Covalent15, covalentMap);
        OpenMM_IntArray_resize(covalentMap, 0);
        for (int j = 0; j < ip11[i].length; j++) {
            OpenMM_IntArray_append(covalentMap, ip11[i][j]);
        }
        OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i, OpenMM_AmoebaMultipoleForce_PolarizationCovalent11, covalentMap);
        OpenMM_IntArray_resize(covalentMap, 0);
    // for (int j = 0; j < ip12[i].length; j++) {
    // OpenMM_IntArray_append(covalentMap, ip12[i][j]);
    // }
    // OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
    // OpenMM_AmoebaMultipoleForce_PolarizationCovalent12, covalentMap);
    // OpenMM_IntArray_resize(covalentMap, 0);
    // 
    // for (int j = 0; j < ip13[i].length; j++) {
    // OpenMM_IntArray_append(covalentMap, ip13[i][j]);
    // }
    // OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
    // OpenMM_AmoebaMultipoleForce_PolarizationCovalent13, covalentMap);
    // OpenMM_IntArray_resize(covalentMap, 0);
    // 
    // OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
    // OpenMM_AmoebaMultipoleForce_PolarizationCovalent14, covalentMap);
    }
    OpenMM_IntArray_destroy(covalentMap);
    OpenMM_System_addForce(system, amoebaMultipoleForce);
    OpenMM_Force_setForceGroup(amoebaMultipoleForce, 1);
    logger.log(Level.INFO, " Added polarizable multipole force.");
    GeneralizedKirkwood gk = super.getGK();
    if (gk != null) {
        addGKForce();
    }
}