C++ NonbondedForce::setPMEParameters方法代码示例

void testTriclinic() {
    // Create a triclinic box containing eight particles.

    System system;
    system.setDefaultPeriodicBoxVectors(Vec3(2.5, 0, 0), Vec3(0.5, 3.0, 0), Vec3(0.7, 0.9, 3.5));
    for (int i = 0; i < 8; i++)
        system.addParticle(1.0);
    NonbondedForce* force = new NonbondedForce();
    system.addForce(force);
    force->setNonbondedMethod(NonbondedForce::PME);
    force->setCutoffDistance(1.0);
    force->setPMEParameters(3.45891, 32, 40, 48);
    for (int i = 0; i < 4; i++)
        force->addParticle(-1, 0.440104, 0.4184); // Cl parameters
    for (int i = 0; i < 4; i++)
        force->addParticle(1, 0.332840, 0.0115897); // Na parameters
    vector<Vec3> positions(8);
    positions[0] = Vec3(1.744, 2.788, 3.162);
    positions[1] = Vec3(1.048, 0.762, 2.340);
    positions[2] = Vec3(2.489, 1.570, 2.817);
    positions[3] = Vec3(1.027, 1.893, 3.271);
    positions[4] = Vec3(0.937, 0.825, 0.009);
    positions[5] = Vec3(2.290, 1.887, 3.352);
    positions[6] = Vec3(1.266, 1.111, 2.894);
    positions[7] = Vec3(0.933, 1.862, 3.490);

    // Compute the forces and energy.

    VerletIntegrator integ(0.001);
    Context context(system, integ, platform);
    context.setPositions(positions);
    State state = context.getState(State::Forces | State::Energy);

    // Compare them to values computed by Gromacs.

    double expectedEnergy = -963.370;
    vector<Vec3> expectedForce(8);
    expectedForce[0] = Vec3(4.25253e+01, -1.23503e+02, 1.22139e+02);
    expectedForce[1] = Vec3(9.74752e+01, 1.68213e+02, 1.93169e+02);
    expectedForce[2] = Vec3(-1.50348e+02, 1.29165e+02, 3.70435e+02);
    expectedForce[3] = Vec3(9.18644e+02, -3.52571e+00, -1.34772e+03);
    expectedForce[4] = Vec3(-1.61193e+02, 9.01528e+01, -7.12904e+01);
    expectedForce[5] = Vec3(2.82630e+02, 2.78029e+01, -3.72864e+02);
    expectedForce[6] = Vec3(-1.47454e+02, -2.14448e+02, -3.55789e+02);
    expectedForce[7] = Vec3(-8.82195e+02, -7.39132e+01, 1.46202e+03);
    for (int i = 0; i < 8; i++) {
        ASSERT_EQUAL_VEC(expectedForce[i], state.getForces()[i], 1e-4);
    }
    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-4);
}

void testPMEParameters() {
    // Create a cloud of random point charges.

    const int numParticles = 51;
    const double boxWidth = 4.7;
    System system;
    system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
    NonbondedForce* force = new NonbondedForce();
    system.addForce(force);
    vector<Vec3> positions(numParticles);
    OpenMM_SFMT::SFMT sfmt;
    init_gen_rand(0, sfmt);

    for (int i = 0; i < numParticles; i++) {
        system.addParticle(1.0);
        force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
        positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
    }
    force->setNonbondedMethod(NonbondedForce::PME);
    ReferencePlatform platform;
    
    // Compute the energy with an error tolerance of 1e-3.

    force->setEwaldErrorTolerance(1e-3);
    VerletIntegrator integrator1(0.01);
    Context context1(system, integrator1, platform);
    context1.setPositions(positions);
    double energy1 = context1.getState(State::Energy).getPotentialEnergy();
    
    // Try again with an error tolerance of 1e-4.

    force->setEwaldErrorTolerance(1e-4);
    VerletIntegrator integrator2(0.01);
    Context context2(system, integrator2, platform);
    context2.setPositions(positions);
    double energy2 = context2.getState(State::Energy).getPotentialEnergy();
    
    // Now explicitly set the parameters.  These should match the values that were
    // used for tolerance 1e-3.

    force->setPMEParameters(2.49291157051793, 32, 32, 32);
    VerletIntegrator integrator3(0.01);
    Context context3(system, integrator3, platform);
    context3.setPositions(positions);
    double energy3 = context3.getState(State::Energy).getPotentialEnergy();
    ASSERT_EQUAL_TOL(energy1, energy3, 1e-6);
    ASSERT(fabs((energy1-energy2)/energy1) > 1e-5);
}

void testSerialization() {
    // Create a Force.

    NonbondedForce force;
    force.setForceGroup(3);
    force.setNonbondedMethod(NonbondedForce::CutoffPeriodic);
    force.setSwitchingDistance(1.5);
    force.setUseSwitchingFunction(true);
    force.setCutoffDistance(2.0);
    force.setEwaldErrorTolerance(1e-3);
    force.setReactionFieldDielectric(50.0);
    force.setUseDispersionCorrection(false);
    double alpha = 0.5;
    int nx = 3, ny = 5, nz = 7;
    force.setPMEParameters(alpha, nx, ny, nz);
    double dalpha = 0.8;
    int dnx = 4, dny = 6, dnz = 7;
    force.setLJPMEParameters(dalpha, dnx, dny, dnz);
    force.addParticle(1, 0.1, 0.01);
    force.addParticle(0.5, 0.2, 0.02);
    force.addParticle(-0.5, 0.3, 0.03);
    force.addException(0, 1, 2, 0.5, 0.1);
    force.addException(1, 2, 0.2, 0.4, 0.2);
    force.addGlobalParameter("scale1", 1.0);
    force.addGlobalParameter("scale2", 2.0);
    force.addParticleParameterOffset("scale1", 2, 1.5, 2.0, 2.5);
    force.addExceptionParameterOffset("scale2", 1, -0.1, -0.2, -0.3);

    // Serialize and then deserialize it.

    stringstream buffer;
    XmlSerializer::serialize<NonbondedForce>(&force, "Force", buffer);
    NonbondedForce* copy = XmlSerializer::deserialize<NonbondedForce>(buffer);

    // Compare the two forces to see if they are identical.

    NonbondedForce& force2 = *copy;
    ASSERT_EQUAL(force.getForceGroup(), force2.getForceGroup());
    ASSERT_EQUAL(force.getNonbondedMethod(), force2.getNonbondedMethod());
    ASSERT_EQUAL(force.getSwitchingDistance(), force2.getSwitchingDistance());
    ASSERT_EQUAL(force.getUseSwitchingFunction(), force2.getUseSwitchingFunction());
    ASSERT_EQUAL(force.getCutoffDistance(), force2.getCutoffDistance());
    ASSERT_EQUAL(force.getEwaldErrorTolerance(), force2.getEwaldErrorTolerance());
    ASSERT_EQUAL(force.getReactionFieldDielectric(), force2.getReactionFieldDielectric());
    ASSERT_EQUAL(force.getUseDispersionCorrection(), force2.getUseDispersionCorrection());
    ASSERT_EQUAL(force.getNumParticles(), force2.getNumParticles());
    ASSERT_EQUAL(force.getNumExceptions(), force2.getNumExceptions());
    ASSERT_EQUAL(force.getNumGlobalParameters(), force2.getNumGlobalParameters());
    ASSERT_EQUAL(force.getNumParticleParameterOffsets(), force2.getNumParticleParameterOffsets());
    ASSERT_EQUAL(force.getNumExceptionParameterOffsets(), force2.getNumExceptionParameterOffsets());
    double alpha2;
    int nx2, ny2, nz2;
    force2.getPMEParameters(alpha2, nx2, ny2, nz2);
    ASSERT_EQUAL(alpha, alpha2);
    ASSERT_EQUAL(nx, nx2);
    ASSERT_EQUAL(ny, ny2);
    ASSERT_EQUAL(nz, nz2);    
    double dalpha2;
    int dnx2, dny2, dnz2;
    force2.getLJPMEParameters(dalpha2, dnx2, dny2, dnz2);
    ASSERT_EQUAL(dalpha, dalpha2);
    ASSERT_EQUAL(dnx, dnx2);
    ASSERT_EQUAL(dny, dny2);
    ASSERT_EQUAL(dnz, dnz2);
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
        ASSERT_EQUAL(force.getGlobalParameterName(i), force2.getGlobalParameterName(i));
        ASSERT_EQUAL(force.getGlobalParameterDefaultValue(i), force2.getGlobalParameterDefaultValue(i));
    }
    for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
        int index1, index2;
        string param1, param2;
        double charge1, sigma1, epsilon1;
        double charge2, sigma2, epsilon2;
        force.getParticleParameterOffset(i, param1, index1, charge1, sigma1, epsilon1);
        force2.getParticleParameterOffset(i, param2, index2, charge2, sigma2, epsilon2);
        ASSERT_EQUAL(index1, index1);
        ASSERT_EQUAL(param1, param2);
        ASSERT_EQUAL(charge1, charge2);
        ASSERT_EQUAL(sigma1, sigma2);
        ASSERT_EQUAL(epsilon1, epsilon2);
    }
    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
        int index1, index2;
        string param1, param2;
        double charge1, sigma1, epsilon1;
        double charge2, sigma2, epsilon2;
        force.getExceptionParameterOffset(i, param1, index1, charge1, sigma1, epsilon1);
        force2.getExceptionParameterOffset(i, param2, index2, charge2, sigma2, epsilon2);
        ASSERT_EQUAL(index1, index1);
        ASSERT_EQUAL(param1, param2);
        ASSERT_EQUAL(charge1, charge2);
        ASSERT_EQUAL(sigma1, sigma2);
        ASSERT_EQUAL(epsilon1, epsilon2);
    }
    for (int i = 0; i < force.getNumParticles(); i++) {
        double charge1, sigma1, epsilon1;
        double charge2, sigma2, epsilon2;
        force.getParticleParameters(i, charge1, sigma1, epsilon1);
        force2.getParticleParameters(i, charge2, sigma2, epsilon2);
        ASSERT_EQUAL(charge1, charge2);
//.........这里部分代码省略.........

void test_water2_dpme_energies_forces_no_exclusions() {
    const double cutoff = 7.0*OpenMM::NmPerAngstrom;
    const double dalpha = 4.0124063605;
    const int grid = 32;
    NonbondedForce* forceField = new NonbondedForce();

    vector<Vec3> positions;
    vector<double> epsvals;
    vector<double> sigvals;
    vector<pair<int, int> > bonds;
    System system;

    const int NATOMS = 6;
    double boxEdgeLength = 25*OpenMM::NmPerAngstrom;

    make_waterbox(NATOMS, boxEdgeLength, forceField,  positions, epsvals, sigvals, bonds, system, false);
    forceField->setNonbondedMethod(OpenMM::NonbondedForce::LJPME);
    forceField->setPMEParameters(0.0f, grid, grid, grid);
    forceField->setReciprocalSpaceForceGroup(1);
    forceField->setLJPMEParameters(dalpha, grid, grid, grid);
    forceField->setCutoffDistance(cutoff);
    forceField->setReactionFieldDielectric(1.0);
    system.addForce(forceField);

    // Reference calculation
    VerletIntegrator integrator(0.01);
    Platform& platform = Platform::getPlatformByName("Reference");
    Context context(system, integrator, platform);
    context.setPositions(positions);
    State state = context.getState(State::Forces | State::Energy, false, 1<<1);
    double refenergy = state.getPotentialEnergy();
    const vector<Vec3>& refforces = state.getForces();

    // Optimized CPU calculation
    CpuCalcDispersionPmeReciprocalForceKernel pme(CalcPmeReciprocalForceKernel::Name(), platform);
    IO io;
    double selfEwaldEnergy = 0;
    double dalpha6 = pow(dalpha, 6.0);
    for (int i = 0; i < NATOMS; i++) {
        io.posq.push_back((float)positions[i][0]);
        io.posq.push_back((float)positions[i][1]);
        io.posq.push_back((float)positions[i][2]);
        double c6 = 8.0f * pow(sigvals[i], 3) * epsvals[i];
        io.posq.push_back(c6);
        selfEwaldEnergy += dalpha6 * c6 * c6 / 12.0;
    }
    pme.initialize(grid, grid, grid, NATOMS, dalpha, false);
    Vec3 boxVectors[3];
    system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
    pme.beginComputation(io, boxVectors, true);
    double recenergy = pme.finishComputation(io);

    ASSERT_EQUAL_TOL(recenergy, -2.179629087, 5e-3);
    ASSERT_EQUAL_TOL(selfEwaldEnergy, 1.731404285, 1e-5);

    std::vector<Vec3> knownforces(6);
    knownforces[0] = Vec3(    -1.890360546,    -1.890723915,    -1.879662698);
    knownforces[1] = Vec3( -0.003161352455, -0.000922244929, -0.005391616425);
    knownforces[2] = Vec3( 0.0009199035545, -0.001453894176, -0.006188087146);
    knownforces[3] = Vec3(     1.887108856,     1.887241446,      1.89644647);
    knownforces[4] = Vec3( 0.0008242336483,  0.003778910089, -0.002116131106);
    knownforces[5] = Vec3(  0.004912763044,  0.002324059399, -0.002844482646);
    for (int i = 0; i < NATOMS; i++)
        ASSERT_EQUAL_VEC(refforces[i], knownforces[i], 5e-3);

    recenergy += selfEwaldEnergy;

    // See if they match.

    ASSERT_EQUAL_TOL(refenergy, recenergy, 1e-3);
    for (int i = 0; i < NATOMS; i++)
        ASSERT_EQUAL_VEC(refforces[i], Vec3(io.force[4*i], io.force[4*i+1], io.force[4*i+2]), 5e-3);

}