Examples with GeneralizedKirkwood ffx.potential.nonbonded.GeneralizedKirkwood used on opensource projects

Example 1 with GeneralizedKirkwood

use of ffx.potential.nonbonded.GeneralizedKirkwood in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method addCustomNonbondedSoftcoreForce.

/**
 * 1.) Handle interactions between non-alchemical atoms with our default
 * OpenMM NonBondedForce. Note that alchemical atoms must have eps=0 to turn
 * them off in this force.
 * <p>
 * 2.) Handle interactions between alchemical atoms and mixed non-alchemical
 * <-> alchemical interactions with an OpenMM CustomNonBondedForce.
 */
private void addCustomNonbondedSoftcoreForce() {
    VanDerWaals vdW = super.getVdwNode();
    if (vdW == null) {
        return;
    }
    /**
     * Only 6-12 LJ with arithmetic mean to define sigma and geometric mean
     * for epsilon is supported.
     */
    VanDerWaalsForm vdwForm = vdW.getVDWForm();
    if (vdwForm.vdwType != LENNARD_JONES || vdwForm.radiusRule != ARITHMETIC || vdwForm.epsilonRule != GEOMETRIC) {
        logger.info(format(" VDW Type:         %s", vdwForm.vdwType));
        logger.info(format(" VDW Radius Rule:  %s", vdwForm.radiusRule));
        logger.info(format(" VDW Epsilon Rule: %s", vdwForm.epsilonRule));
        logger.log(Level.SEVERE, String.format(" Unsuppporterd van der Waals functional form."));
        return;
    }
    // Sterics mixing rules.
    String stericsMixingRules = " epsilon = sqrt(epsilon1*epsilon2);";
    stericsMixingRules += " rmin = 0.5 * (sigma1 + sigma2) * 1.122462048309372981;";
    // Softcore Lennard-Jones, with a form equivalent to that used in FFX VanDerWaals class.
    String stericsEnergyExpression = "(vdw_lambda^beta)*epsilon*x*(x-2.0);";
    // Effective softcore distance for sterics.
    stericsEnergyExpression += " x = 1.0 / (alpha*(1.0-vdw_lambda)^2.0 + (r/rmin)^6.0);";
    // Define energy expression for sterics.
    String energyExpression = stericsEnergyExpression + stericsMixingRules;
    fixedChargeSoftcore = OpenMM_CustomNonbondedForce_create(energyExpression);
    // Get the Alpha and Beta constants from the VanDerWaals instance.
    double alpha = vdW.getAlpha();
    double beta = vdW.getBeta();
    logger.info(format(" Custom non-bonded force with alpha = %8.6f and beta = %8.6f", alpha, beta));
    OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "vdw_lambda", 1.0);
    OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "alpha", alpha);
    OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "beta", beta);
    OpenMM_CustomNonbondedForce_addPerParticleParameter(fixedChargeSoftcore, "sigma");
    OpenMM_CustomNonbondedForce_addPerParticleParameter(fixedChargeSoftcore, "epsilon");
    /**
     * Add particles.
     */
    PointerByReference alchemicalGroup = OpenMM_IntSet_create();
    PointerByReference nonAlchemicalGroup = OpenMM_IntSet_create();
    DoubleByReference charge = new DoubleByReference();
    DoubleByReference sigma = new DoubleByReference();
    DoubleByReference eps = new DoubleByReference();
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        Atom atom = atoms[i];
        OpenMM_NonbondedForce_getParticleParameters(fixedChargeNonBondedForce, i, charge, sigma, eps);
        if (atom.applyLambda()) {
            OpenMM_IntSet_insert(alchemicalGroup, i);
            logger.info(format(" Adding alchemical atom %s.", atom));
        } else {
            OpenMM_IntSet_insert(nonAlchemicalGroup, i);
        }
        PointerByReference particleParameters = OpenMM_DoubleArray_create(0);
        OpenMM_DoubleArray_append(particleParameters, sigma.getValue());
        OpenMM_DoubleArray_append(particleParameters, eps.getValue());
        OpenMM_CustomNonbondedForce_addParticle(fixedChargeSoftcore, particleParameters);
        OpenMM_DoubleArray_destroy(particleParameters);
    }
    OpenMM_CustomNonbondedForce_addInteractionGroup(fixedChargeSoftcore, alchemicalGroup, alchemicalGroup);
    OpenMM_CustomNonbondedForce_addInteractionGroup(fixedChargeSoftcore, alchemicalGroup, nonAlchemicalGroup);
    OpenMM_IntSet_destroy(alchemicalGroup);
    OpenMM_IntSet_destroy(nonAlchemicalGroup);
    Crystal crystal = super.getCrystal();
    if (crystal.aperiodic()) {
        OpenMM_CustomNonbondedForce_setNonbondedMethod(fixedChargeSoftcore, OpenMM_CustomNonbondedForce_NonbondedMethod.OpenMM_CustomNonbondedForce_NoCutoff);
    } else {
        OpenMM_CustomNonbondedForce_setNonbondedMethod(fixedChargeSoftcore, OpenMM_CustomNonbondedForce_NonbondedMethod.OpenMM_CustomNonbondedForce_CutoffPeriodic);
    }
    NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
    double off = nonbondedCutoff.off;
    double cut = nonbondedCutoff.cut;
    OpenMM_CustomNonbondedForce_setCutoffDistance(fixedChargeSoftcore, OpenMM_NmPerAngstrom * off);
    OpenMM_CustomNonbondedForce_setUseSwitchingFunction(fixedChargeSoftcore, OpenMM_True);
    OpenMM_CustomNonbondedForce_setSwitchingDistance(fixedChargeSoftcore, OpenMM_NmPerAngstrom * cut);
    if (cut == off) {
        logger.warning(" OpenMM does not properly handle cutoffs where cut == off!");
        if (cut == Double.MAX_VALUE || cut == Double.POSITIVE_INFINITY) {
            logger.info(" Detected infinite or max-value cutoff; setting cut to 1E+40 for OpenMM.");
            cut = 1E40;
        } else {
            logger.info(String.format(" Detected cut %8.4g == off %8.4g; scaling cut to 0.99 of off for OpenMM.", cut, off));
            cut *= 0.99;
        }
    }
    // Add energy parameter derivative
    OpenMM_CustomNonbondedForce_addEnergyParameterDerivative(fixedChargeSoftcore, "vdw_lambda");
    OpenMM_System_addForce(system, fixedChargeSoftcore);
    logger.log(Level.INFO, String.format(" Added fixed charge softcore sterics force."));
    GeneralizedKirkwood gk = super.getGK();
    if (gk != null) {
        logger.severe(" OpenMM alchemical methods are not supported for GB.");
        addCustomGBForce();
    }
    // Not entirely sure how to initialize this portion
    alchemicalAlchemicalStericsForce = OpenMM_CustomBondForce_create(stericsEnergyExpression);
    nonAlchemicalAlchemicalStericsForce = OpenMM_CustomBondForce_create(stericsEnergyExpression);
    // allStericsForce = (alchemicalAlchemicalStericsForce + nonAlchemicalAlchemicalStericsForce);
    // Can be reduced to two lines if I can figure out how to combine the two custom bonded sterics forces
    OpenMM_CustomBondForce_addPerBondParameter(alchemicalAlchemicalStericsForce, "rmin");
    OpenMM_CustomBondForce_addPerBondParameter(alchemicalAlchemicalStericsForce, "epsilon");
    OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "vdw_lambda", 1.0);
    OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "alpha", alpha);
    OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "beta", beta);
    OpenMM_CustomBondForce_addPerBondParameter(nonAlchemicalAlchemicalStericsForce, "rmin");
    OpenMM_CustomBondForce_addPerBondParameter(nonAlchemicalAlchemicalStericsForce, "epsilon");
    OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "vdw_lambda", 1.0);
    OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "alpha", alpha);
    OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "beta", beta);
    int range = OpenMM_NonbondedForce_getNumExceptions(fixedChargeNonBondedForce);
    IntByReference atomi = new IntByReference();
    IntByReference atomj = new IntByReference();
    int[][] torsionMask = vdW.getTorsionMask();
    for (int i = 0; i < range; i++) {
        OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, atomi, atomj, charge, sigma, eps);
        OpenMM_CustomNonbondedForce_addExclusion(fixedChargeSoftcore, atomi.getValue(), atomj.getValue());
        int[] maskI = torsionMask[atomi.getValue()];
        int jID = atomj.getValue();
        boolean epsException = false;
        for (int j = 0; j < maskI.length; j++) {
            if (maskI[j] == jID) {
                epsException = true;
                break;
            }
        }
        Atom atom1 = atoms[atomi.getValue()];
        Atom atom2 = atoms[atomj.getValue()];
        boolean bothAlchemical = false;
        boolean oneAlchemical = false;
        if (atom1.applyLambda() && atom2.applyLambda()) {
            bothAlchemical = true;
        } else if ((atom1.applyLambda() && !atom2.applyLambda()) || (!atom1.applyLambda() && atom2.applyLambda())) {
            oneAlchemical = true;
        }
        if (bothAlchemical) {
            if (epsException) {
                PointerByReference bondParameters = OpenMM_DoubleArray_create(0);
                OpenMM_DoubleArray_append(bondParameters, sigma.getValue() * 1.122462048309372981);
                OpenMM_DoubleArray_append(bondParameters, eps.getValue());
                OpenMM_CustomBondForce_addBond(alchemicalAlchemicalStericsForce, atomi.getValue(), atomj.getValue(), bondParameters);
                OpenMM_DoubleArray_destroy(bondParameters);
            }
        } else if (oneAlchemical) {
            if (epsException) {
                PointerByReference bondParameters = OpenMM_DoubleArray_create(0);
                OpenMM_DoubleArray_append(bondParameters, sigma.getValue() * 1.122462048309372981);
                OpenMM_DoubleArray_append(bondParameters, eps.getValue());
                OpenMM_CustomBondForce_addBond(nonAlchemicalAlchemicalStericsForce, atomi.getValue(), atomj.getValue(), bondParameters);
                OpenMM_DoubleArray_destroy(bondParameters);
            }
        }
    }
    /**
     *        for (int i = 0; i < range; i++){
     *            OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, atomi, atomj, charge, sigma, eps);
     *
     *            Atom atom1 = atoms[atomi.getValue()];
     *            Atom atom2 = atoms[atomj.getValue()];
     *
     *            if (atom1.applyLambda() || atom2.applyLambda()){
     *                OpenMM_NonbondedForce_setExceptionParameters(fixedChargeNonBondedForce, i, atomi.getValue(), atomj.getValue(), abs(0.0*charge.getValue()), sigma.getValue(), abs(0.0*eps.getValue()));
     *            }
     *        }
     */
    OpenMM_CustomBondForce_addEnergyParameterDerivative(alchemicalAlchemicalStericsForce, "vdw_lambda");
    OpenMM_CustomBondForce_addEnergyParameterDerivative(nonAlchemicalAlchemicalStericsForce, "vdw_lambda");
    OpenMM_System_addForce(system, alchemicalAlchemicalStericsForce);
    OpenMM_System_addForce(system, nonAlchemicalAlchemicalStericsForce);
}

Example 2 with GeneralizedKirkwood

use of ffx.potential.nonbonded.GeneralizedKirkwood in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method updateAmoebaGeneralizedKirkwoodForce.

/**
 * Updates the AMOEBA Generalized Kirkwood force for change in Use flags.
 *
 * @param atoms Array of all Atoms in the system
 */
private void updateAmoebaGeneralizedKirkwoodForce(Atom[] atoms) {
    GeneralizedKirkwood gk = super.getGK();
    double[] overlapScale = gk.getOverlapScale();
    double[] baseRadii = gk.getBaseRadii();
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        double useFactor = 1.0;
        if (!atoms[i].getUse() || !atoms[i].getElectrostatics()) {
            // if (!atoms[i].getUse()) {
            useFactor = 0.0;
        }
        double lambdaScale = lambda;
        if (!atoms[i].applyLambda()) {
            lambdaScale = 1.0;
        }
        useFactor *= lambdaScale;
        MultipoleType multipoleType = atoms[i].getMultipoleType();
        OpenMM_AmoebaGeneralizedKirkwoodForce_setParticleParameters(amoebaGeneralizedKirkwoodForce, i, multipoleType.charge * useFactor, OpenMM_NmPerAngstrom * baseRadii[i], overlapScale[i] * useFactor);
    }
    OpenMM_AmoebaGeneralizedKirkwoodForce_updateParametersInContext(amoebaGeneralizedKirkwoodForce, context);
}

Example 3 with GeneralizedKirkwood

use of ffx.potential.nonbonded.GeneralizedKirkwood in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method addCustomGBForce.

private void addCustomGBForce() {
    GeneralizedKirkwood gk = super.getGK();
    if (gk == null) {
        return;
    }
    customGBForce = OpenMM_CustomGBForce_create();
    OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "q");
    OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "radius");
    OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "scale");
    OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "solventDielectric", 78.3);
    OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "soluteDielectric", 1.0);
    // Factor of 0.1 for Ang to nm.
    OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "dOffset", gk.getDielecOffset() * OpenMM_NmPerAngstrom);
    OpenMM_CustomGBForce_addComputedValue(customGBForce, "I", // "step(r+sr2-or1)*0.5*((1/L^3-1/U^3)/3+(1/U^4-1/L^4)/8*(r-sr2*sr2/r)+0.25*(1/U^2-1/L^2)/r+C);"
    "0.5*((1/L^3-1/U^3)/3.0+(1/U^4-1/L^4)/8.0*(r-sr2*sr2/r)+0.25*(1/U^2-1/L^2)/r+C);" + "U=r+sr2;" + // + "C=2*(1/or1-1/L)*step(sr2-r-or1);"
    "C=2/3*(1/or1^3-1/L^3)*step(sr2-r-or1);" + // + "D=step(r-sr2)*(r-sr2) + (1-step(r-sr2))*(sr2-r);"
    "L = step(sr2 - r1r)*sr2mr + (1 - step(sr2 - r1r))*L;" + "sr2mr = sr2 - r;" + "r1r = radius1 + r;" + "L = step(r1sr2 - r)*radius1 + (1 - step(r1sr2 - r))*L;" + "r1sr2 = radius1 + sr2;" + "L = r - sr2;" + "sr2 = scale2 * radius2;" + "or1 = radius1; or2 = radius2", OpenMM_CustomGBForce_ParticlePairNoExclusions);
    OpenMM_CustomGBForce_addComputedValue(customGBForce, "B", // "psi=I*or; or=radius-0.009"
    "step(BB-radius)*BB + (1 - step(BB-radius))*radius;" + "BB = 1 / ( (3.0*III)^(1.0/3.0) );" + "III = step(II)*II + (1 - step(II))*1.0e-9/3.0;" + "II = maxI - I;" + "maxI = 1/(3.0*radius^3)", OpenMM_CustomGBForce_SingleParticle);
    double sTens = gk.getSurfaceTension();
    logger.info(String.format(" FFX surface tension: %9.5g kcal/mol/Ang^2", sTens));
    sTens *= OpenMM_KJPerKcal;
    // 100 square Angstroms per square nanometer.
    sTens *= 100.0;
    logger.info(String.format(" OpenMM surface tension: %9.5g kJ/mol/nm^2", sTens));
    String surfaceTension = Double.toString(sTens);
    OpenMM_CustomGBForce_addEnergyTerm(customGBForce, surfaceTension + "*(radius+0.14+dOffset)^2*((radius+dOffset)/B)^6/6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B", OpenMM_CustomGBForce_SingleParticle);
    /**
     * Particle pair term is the generalized Born cross term.
     */
    OpenMM_CustomGBForce_addEnergyTerm(customGBForce, "-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;" + "f=sqrt(r^2+B1*B2*exp(-r^2/(2.455*B1*B2)))", OpenMM_CustomGBForce_ParticlePair);
    double[] baseRadii = gk.getBaseRadii();
    double[] overlapScale = gk.getOverlapScale();
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    PointerByReference doubleArray = OpenMM_DoubleArray_create(0);
    for (int i = 0; i < nAtoms; i++) {
        MultipoleType multipoleType = atoms[i].getMultipoleType();
        OpenMM_DoubleArray_append(doubleArray, multipoleType.charge);
        OpenMM_DoubleArray_append(doubleArray, OpenMM_NmPerAngstrom * baseRadii[i]);
        OpenMM_DoubleArray_append(doubleArray, overlapScale[i]);
        OpenMM_CustomGBForce_addParticle(customGBForce, doubleArray);
        OpenMM_DoubleArray_resize(doubleArray, 0);
    }
    OpenMM_DoubleArray_destroy(doubleArray);
    double cut = gk.getCutoff();
    OpenMM_CustomGBForce_setCutoffDistance(customGBForce, cut);
    OpenMM_Force_setForceGroup(customGBForce, 1);
    OpenMM_System_addForce(system, customGBForce);
    logger.log(Level.INFO, " Added generalized Born force");
}

Example 4 with GeneralizedKirkwood

use of ffx.potential.nonbonded.GeneralizedKirkwood in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method addGKForce.

private void addGKForce() {
    GeneralizedKirkwood gk = super.getGK();
    amoebaGeneralizedKirkwoodForce = OpenMM_AmoebaGeneralizedKirkwoodForce_create();
    OpenMM_AmoebaGeneralizedKirkwoodForce_setSolventDielectric(amoebaGeneralizedKirkwoodForce, 78.3);
    OpenMM_AmoebaGeneralizedKirkwoodForce_setSoluteDielectric(amoebaGeneralizedKirkwoodForce, 1.0);
    double[] overlapScale = gk.getOverlapScale();
    double[] baseRadii = gk.getBaseRadii();
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        MultipoleType multipoleType = atoms[i].getMultipoleType();
        OpenMM_AmoebaGeneralizedKirkwoodForce_addParticle(amoebaGeneralizedKirkwoodForce, multipoleType.charge, OpenMM_NmPerAngstrom * baseRadii[i], overlapScale[i]);
    }
    OpenMM_AmoebaGeneralizedKirkwoodForce_setProbeRadius(amoebaGeneralizedKirkwoodForce, 1.4 * OpenMM_NmPerAngstrom);
    NonPolar nonpolar = gk.getNonPolarModel();
    switch(nonpolar) {
        case BORN_SOLV:
        case BORN_CAV_DISP:
        default:
            // Configure a Born Radii based surface area term.
            double surfaceTension = gk.getSurfaceTension() * OpenMM_KJPerKcal * OpenMM_AngstromsPerNm * OpenMM_AngstromsPerNm;
            OpenMM_AmoebaGeneralizedKirkwoodForce_setIncludeCavityTerm(amoebaGeneralizedKirkwoodForce, OpenMM_True);
            OpenMM_AmoebaGeneralizedKirkwoodForce_setSurfaceAreaFactor(amoebaGeneralizedKirkwoodForce, -surfaceTension);
            break;
        case CAV:
        case CAV_DISP:
        case HYDROPHOBIC_PMF:
        case NONE:
            // This NonPolar model does not use a Born Radii based surface area term.
            OpenMM_AmoebaGeneralizedKirkwoodForce_setIncludeCavityTerm(amoebaGeneralizedKirkwoodForce, OpenMM_False);
            break;
    }
    OpenMM_System_addForce(system, amoebaGeneralizedKirkwoodForce);
    switch(nonpolar) {
        case CAV_DISP:
        case BORN_CAV_DISP:
            addWCAForce();
            break;
        case CAV:
        case HYDROPHOBIC_PMF:
        case BORN_SOLV:
        case NONE:
        default:
    }
    logger.log(Level.INFO, " Added generalized Kirkwood force.");
}

Example 5 with GeneralizedKirkwood

use of ffx.potential.nonbonded.GeneralizedKirkwood 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();
    }
}