C++ CustomNonbondedForce::addExclusion方法代码示例

void testExclusions() {
    System system;
    VerletIntegrator integrator(0.01);
    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("a*r; a=a1+a2");
    nonbonded->addPerParticleParameter("a");
    vector<double> params(1);
    vector<Vec3> positions(4);
    for (int i = 0; i < 4; i++) {
        system.addParticle(1.0);
        params[0] = i+1;
        nonbonded->addParticle(params);
        positions[i] = Vec3(i, 0, 0);
    }
    nonbonded->addExclusion(0, 1);
    nonbonded->addExclusion(1, 2);
    nonbonded->addExclusion(2, 3);
    nonbonded->addExclusion(0, 2);
    nonbonded->addExclusion(1, 3);
    system.addForce(nonbonded);
    Context context(system, integrator, platform);
    context.setPositions(positions);
    State state = context.getState(State::Forces | State::Energy);
    const vector<Vec3>& forces = state.getForces();
    ASSERT_EQUAL_VEC(Vec3(1+4, 0, 0), forces[0], TOL);
    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
    ASSERT_EQUAL_VEC(Vec3(-(1+4), 0, 0), forces[3], TOL);
    ASSERT_EQUAL_TOL((1+4)*3.0, state.getPotentialEnergy(), TOL);
}

void testInteractionGroups() {
    const int numParticles = 6;
    System system;
    VerletIntegrator integrator(0.01);
    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("v1+v2");
    nonbonded->addPerParticleParameter("v");
    vector<double> params(1, 0.001);
    for (int i = 0; i < numParticles; i++) {
        system.addParticle(1.0);
        nonbonded->addParticle(params);
        params[0] *= 10;
    }
    set<int> set1, set2, set3, set4;
    set1.insert(2);
    set2.insert(0);
    set2.insert(1);
    set2.insert(2);
    set2.insert(3);
    set2.insert(4);
    set2.insert(5);
    nonbonded->addInteractionGroup(set1, set2); // Particle 2 interacts with every other particle.
    set3.insert(0);
    set3.insert(1);
    set4.insert(4);
    set4.insert(5);
    nonbonded->addInteractionGroup(set3, set4); // Particles 0 and 1 interact with 4 and 5.
    nonbonded->addExclusion(1, 2); // Add an exclusion to make sure it gets skipped.
    system.addForce(nonbonded);
    Context context(system, integrator, platform);
    vector<Vec3> positions(numParticles);
    context.setPositions(positions);
    State state = context.getState(State::Energy);
    double expectedEnergy = 331.423; // Each digit is the number of interactions a particle particle is involved in.
    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), TOL);
}

void* CustomNonbondedForceProxy::deserialize(const SerializationNode& node) const {
    if (node.getIntProperty("version") != 1)
        throw OpenMMException("Unsupported version number");
    CustomNonbondedForce* force = NULL;
    try {
        CustomNonbondedForce* force = new CustomNonbondedForce(node.getStringProperty("energy"));
        force->setNonbondedMethod((CustomNonbondedForce::NonbondedMethod) node.getIntProperty("method"));
        force->setCutoffDistance(node.getDoubleProperty("cutoff"));
        const SerializationNode& perParticleParams = node.getChildNode("PerParticleParameters");
        for (int i = 0; i < (int) perParticleParams.getChildren().size(); i++) {
            const SerializationNode& parameter = perParticleParams.getChildren()[i];
            force->addPerParticleParameter(parameter.getStringProperty("name"));
        }
        const SerializationNode& globalParams = node.getChildNode("GlobalParameters");
        for (int i = 0; i < (int) globalParams.getChildren().size(); i++) {
            const SerializationNode& parameter = globalParams.getChildren()[i];
            force->addGlobalParameter(parameter.getStringProperty("name"), parameter.getDoubleProperty("default"));
        }
        const SerializationNode& particles = node.getChildNode("Particles");
        vector<double> params(force->getNumPerParticleParameters());
        for (int i = 0; i < (int) particles.getChildren().size(); i++) {
            const SerializationNode& particle = particles.getChildren()[i];
            for (int j = 0; j < (int) params.size(); j++) {
                stringstream key;
                key << "param";
                key << j+1;
                params[j] = particle.getDoubleProperty(key.str());
            }
            force->addParticle(params);
        }
        const SerializationNode& exclusions = node.getChildNode("Exclusions");
        for (int i = 0; i < (int) exclusions.getChildren().size(); i++) {
            const SerializationNode& exclusion = exclusions.getChildren()[i];
            force->addExclusion(exclusion.getIntProperty("p1"), exclusion.getIntProperty("p2"));
        }
        const SerializationNode& functions = node.getChildNode("Functions");
        for (int i = 0; i < (int) functions.getChildren().size(); i++) {
            const SerializationNode& function = functions.getChildren()[i];
            const SerializationNode& valuesNode = function.getChildNode("Values");
            vector<double> values;
            for (int j = 0; j < (int) valuesNode.getChildren().size(); j++)
                values.push_back(valuesNode.getChildren()[j].getDoubleProperty("v"));
            force->addFunction(function.getStringProperty("name"), values, function.getDoubleProperty("min"), function.getDoubleProperty("max"));
        }
        return force;
    }
    catch (...) {
        if (force != NULL)
            delete force;
        throw;
    }
}

void testParallelComputation() {
    System system;
    const int numParticles = 200;
    for (int i = 0; i < numParticles; i++)
        system.addParticle(1.0);
    CustomNonbondedForce* force = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6); sigma=0.5; eps=1");
    vector<double> params;
    for (int i = 0; i < numParticles; i++)
        force->addParticle(params);
    system.addForce(force);
    OpenMM_SFMT::SFMT sfmt;
    init_gen_rand(0, sfmt);
    vector<Vec3> positions(numParticles);
    for (int i = 0; i < numParticles; i++)
        positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt));
    for (int i = 0; i < numParticles; ++i)
        for (int j = 0; j < i; ++j) {
            Vec3 delta = positions[i]-positions[j];
            if (delta.dot(delta) < 0.1)
                force->addExclusion(i, j);
        }
    VerletIntegrator integrator1(0.01);
    Context context1(system, integrator1, platform);
    context1.setPositions(positions);
    State state1 = context1.getState(State::Forces | State::Energy);
    VerletIntegrator integrator2(0.01);
    string deviceIndex = platform.getPropertyValue(context1, CudaPlatform::CudaDeviceIndex());
    map<string, string> props;
    props[CudaPlatform::CudaDeviceIndex()] = deviceIndex+","+deviceIndex;
    Context context2(system, integrator2, platform, props);
    context2.setPositions(positions);
    State state2 = context2.getState(State::Forces | State::Energy);
    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
    for (int i = 0; i < numParticles; i++)
        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
}

void testLargeInteractionGroup() {
    const int numMolecules = 300;
    const int numParticles = numMolecules*2;
    const double boxSize = 20.0;
    
    // Create a large system.
    
    System system;
    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
    for (int i = 0; i < numParticles; i++)
        system.addParticle(1.0);
    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
    nonbonded->addPerParticleParameter("q");
    nonbonded->addPerParticleParameter("sigma");
    nonbonded->addPerParticleParameter("eps");
    vector<Vec3> positions(numParticles);
    OpenMM_SFMT::SFMT sfmt;
    init_gen_rand(0, sfmt);
    vector<double> params(3);
    for (int i = 0; i < numMolecules; i++) {
        if (i < numMolecules/2) {
            params[0] = 1.0;
            params[1] = 0.2;
            params[2] = 0.1;
            nonbonded->addParticle(params);
            params[0] = -1.0;
            params[1] = 0.1;
            nonbonded->addParticle(params);
        }
        else {
            params[0] = 1.0;
            params[1] = 0.2;
            params[2] = 0.2;
            nonbonded->addParticle(params);
            params[0] = -1.0;
            params[1] = 0.1;
            nonbonded->addParticle(params);
        }
        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
        nonbonded->addExclusion(2*i, 2*i+1);
    }
    nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
    system.addForce(nonbonded);
    
    // Compute the forces.
    
    VerletIntegrator integrator(0.01);
    Context context(system, integrator, platform);
    context.setPositions(positions);
    State state1 = context.getState(State::Forces);
    
    // Modify the force so only one particle interacts with everything else.
    
    set<int> set1, set2;
    set1.insert(151);
    for (int i = 0; i < numParticles; i++)
        set2.insert(i);
    nonbonded->addInteractionGroup(set1, set2);
    context.reinitialize();
    context.setPositions(positions);
    State state2 = context.getState(State::Forces);
    
    // The force on that one particle should be the same.
    
    ASSERT_EQUAL_VEC(state1.getForces()[151], state2.getForces()[151], 1e-4);
    
    // Modify the interaction group so it includes all interactions.  This should now reproduce the original forces
    // on all atoms.

    for (int i = 0; i < numParticles; i++)
        set1.insert(i);
    nonbonded->setInteractionGroupParameters(0, set1, set2);
    context.reinitialize();
    context.setPositions(positions);
    State state3 = context.getState(State::Forces);
    for (int i = 0; i < numParticles; i++)
        ASSERT_EQUAL_VEC(state1.getForces()[i], state3.getForces()[i], 1e-4);
}

void testCoulombLennardJones() {
    const int numMolecules = 300;
    const int numParticles = numMolecules*2;
    const double boxSize = 20.0;

    // Create two systems: one with a NonbondedForce, and one using a CustomNonbondedForce to implement the same interaction.

    System standardSystem;
    System customSystem;
    for (int i = 0; i < numParticles; i++) {
        standardSystem.addParticle(1.0);
        customSystem.addParticle(1.0);
    }
    NonbondedForce* standardNonbonded = new NonbondedForce();
    CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
    customNonbonded->addPerParticleParameter("q");
    customNonbonded->addPerParticleParameter("sigma");
    customNonbonded->addPerParticleParameter("eps");
    vector<Vec3> positions(numParticles);
    vector<Vec3> velocities(numParticles);
    OpenMM_SFMT::SFMT sfmt;
    init_gen_rand(0, sfmt);

    vector<double> params(3);
    for (int i = 0; i < numMolecules; i++) {
        if (i < numMolecules/2) {
            standardNonbonded->addParticle(1.0, 0.2, 0.1);
            params[0] = 1.0;
            params[1] = 0.2;
            params[2] = 0.1;
            customNonbonded->addParticle(params);
            standardNonbonded->addParticle(-1.0, 0.1, 0.1);
            params[0] = -1.0;
            params[1] = 0.1;
            customNonbonded->addParticle(params);
        }
        else {
            standardNonbonded->addParticle(1.0, 0.2, 0.2);
            params[0] = 1.0;
            params[1] = 0.2;
            params[2] = 0.2;
            customNonbonded->addParticle(params);
            standardNonbonded->addParticle(-1.0, 0.1, 0.2);
            params[0] = -1.0;
            params[1] = 0.1;
            customNonbonded->addParticle(params);
        }
        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
        velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
        velocities[2*i+1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
        standardNonbonded->addException(2*i, 2*i+1, 0.0, 1.0, 0.0);
        customNonbonded->addExclusion(2*i, 2*i+1);
    }
    standardNonbonded->setNonbondedMethod(NonbondedForce::NoCutoff);
    customNonbonded->setNonbondedMethod(CustomNonbondedForce::NoCutoff);
    standardSystem.addForce(standardNonbonded);
    customSystem.addForce(customNonbonded);
    VerletIntegrator integrator1(0.01);
    VerletIntegrator integrator2(0.01);
    Context context1(standardSystem, integrator1, platform);
    Context context2(customSystem, integrator2, platform);
    context1.setPositions(positions);
    context2.setPositions(positions);
    context1.setVelocities(velocities);
    context2.setVelocities(velocities);
    State state1 = context1.getState(State::Forces | State::Energy);
    State state2 = context2.getState(State::Forces | State::Energy);
    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
    for (int i = 0; i < numParticles; i++) {
        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
    }
}

void* CustomNonbondedForceProxy::deserialize(const SerializationNode& node) const {
    int version = node.getIntProperty("version");
    if (version < 1 || version > 2)
        throw OpenMMException("Unsupported version number");
    CustomNonbondedForce* force = NULL;
    try {
        CustomNonbondedForce* force = new CustomNonbondedForce(node.getStringProperty("energy"));
        force->setForceGroup(node.getIntProperty("forceGroup", 0));
        force->setNonbondedMethod((CustomNonbondedForce::NonbondedMethod) node.getIntProperty("method"));
        force->setCutoffDistance(node.getDoubleProperty("cutoff"));
        force->setUseSwitchingFunction(node.getBoolProperty("useSwitchingFunction", false));
        force->setSwitchingDistance(node.getDoubleProperty("switchingDistance", -1.0));
        force->setUseLongRangeCorrection(node.getBoolProperty("useLongRangeCorrection", false));
        const SerializationNode& perParticleParams = node.getChildNode("PerParticleParameters");
        for (int i = 0; i < (int) perParticleParams.getChildren().size(); i++) {
            const SerializationNode& parameter = perParticleParams.getChildren()[i];
            force->addPerParticleParameter(parameter.getStringProperty("name"));
        }
        const SerializationNode& globalParams = node.getChildNode("GlobalParameters");
        for (int i = 0; i < (int) globalParams.getChildren().size(); i++) {
            const SerializationNode& parameter = globalParams.getChildren()[i];
            force->addGlobalParameter(parameter.getStringProperty("name"), parameter.getDoubleProperty("default"));
        }
        if (version > 1) {
            const SerializationNode& energyDerivs = node.getChildNode("EnergyParameterDerivatives");
            for (int i = 0; i < (int) energyDerivs.getChildren().size(); i++) {
                const SerializationNode& parameter = energyDerivs.getChildren()[i];
                force->addEnergyParameterDerivative(parameter.getStringProperty("name"));
            }
        }
        const SerializationNode& particles = node.getChildNode("Particles");
        vector<double> params(force->getNumPerParticleParameters());
        for (int i = 0; i < (int) particles.getChildren().size(); i++) {
            const SerializationNode& particle = particles.getChildren()[i];
            for (int j = 0; j < (int) params.size(); j++) {
                stringstream key;
                key << "param";
                key << j+1;
                params[j] = particle.getDoubleProperty(key.str());
            }
            force->addParticle(params);
        }
        const SerializationNode& exclusions = node.getChildNode("Exclusions");
        for (int i = 0; i < (int) exclusions.getChildren().size(); i++) {
            const SerializationNode& exclusion = exclusions.getChildren()[i];
            force->addExclusion(exclusion.getIntProperty("p1"), exclusion.getIntProperty("p2"));
        }
        const SerializationNode& functions = node.getChildNode("Functions");
        for (int i = 0; i < (int) functions.getChildren().size(); i++) {
            const SerializationNode& function = functions.getChildren()[i];
            if (function.hasProperty("type")) {
                force->addTabulatedFunction(function.getStringProperty("name"), function.decodeObject<TabulatedFunction>());
            }
            else {
                // This is an old file created before TabulatedFunction existed.

                const SerializationNode& valuesNode = function.getChildNode("Values");
                vector<double> values;
                for (int j = 0; j < (int) valuesNode.getChildren().size(); j++)
                    values.push_back(valuesNode.getChildren()[j].getDoubleProperty("v"));
                force->addTabulatedFunction(function.getStringProperty("name"), new Continuous1DFunction(values, function.getDoubleProperty("min"), function.getDoubleProperty("max")));
            }
        }
        bool hasInteractionGroups = false; // Older files will be missing this block.
        for (int i = 0; i < (int) node.getChildren().size(); i++) {
            if (node.getChildren()[i].getName() == "InteractionGroups")
                hasInteractionGroups = true;
        }
        if (hasInteractionGroups) {
            const SerializationNode& interactionGroups = node.getChildNode("InteractionGroups");
            for (int i = 0; i < (int) interactionGroups.getChildren().size(); i++) {
                const SerializationNode& interactionGroup = interactionGroups.getChildren()[i];
                // Get set 1.
                const SerializationNode& set1node = interactionGroup.getChildNode("Set1");
                std::set<int> set1;
                for (int j = 0; j < (int) set1node.getChildren().size(); j++)
                    set1.insert(set1node.getChildren()[j].getIntProperty("index"));
                // Get set 2.
                const SerializationNode& set2node = interactionGroup.getChildNode("Set2");
                std::set<int> set2;
                for (int j = 0; j < (int) set2node.getChildren().size(); j++)
                    set2.insert(set2node.getChildren()[j].getIntProperty("index"));
                force->addInteractionGroup(set1, set2);
            }
        }
        return force;
    }
    catch (...) {
        if (force != NULL)
            delete force;
        throw;
    }
}