TestNoseHooverThermostat.h 12.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2019 Stanford University and the Authors.           *
 * Authors: Andreas Krämer and Andrew C. Simmonett                            *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

#include "openmm/internal/AssertionUtilities.h"
#include "openmm/NoseHooverChain.h"
#include "openmm/VelocityVerletIntegrator.h"
#include "openmm/Context.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
38
#include "openmm/CustomExternalForce.h"
39
40
41
42
#include "openmm/System.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
43
#include <iomanip>
44
#include <vector>
45
#include <algorithm>
46
47
48

using namespace OpenMM;
using namespace std;
49
#define DEBUG 0
50
51
52
53
void testHarmonicOscillator() {
    const double mass = 1.0;
    double temperature = 300;
    double frequency = 1;
54
55
    double mts = 1, ys = 1, chain_length = 3;

56
57
58
    System system;    
    system.addParticle(mass);
    vector<Vec3> positions(1);
59
    positions[0] = Vec3(0.5,0.5,0.5);
60
61
    vector<Vec3> velocities(1);
    velocities[0] = Vec3(0, 0, 0);
62
    auto harmonic_restraint = new CustomExternalForce("0.5*(x^2+y^2+z^2)");
63
64
65
66
67
68
69
70
    harmonic_restraint->addParticle(0);
    system.addForce(harmonic_restraint);
    VelocityVerletIntegrator integrator(0.001);
    integrator.addNoseHooverChainThermostat(system, temperature, frequency, chain_length, mts, ys);
    Context context(system, integrator, platform);
    context.setPositions(positions);
    context.setVelocities(velocities);

71
    double mean_temperature=0;
72
    // equilibration
73
74
    integrator.step(2000);
    for (size_t i=0; i < 2500; i++){
75
        integrator.step(10);
76
        State state = context.getState(State::Energy | State::Positions | State::Velocities);
77
78
79
        double kinetic_energy = state.getKineticEnergy();
        double temp = kinetic_energy/(0.5*3*BOLTZ);
        mean_temperature = (i*mean_temperature + temp)/(i+1);
80
81
82
83
84
85
86
87
88
89
90
        double PE = state.getPotentialEnergy();
        double time = state.getTime();
        double energy = kinetic_energy + PE + integrator.computeHeatBathEnergy();
#if DEBUG
        std::cout << " Time: " << std::setw(4) << time 
                  << " Temperature: " << std::setw(8) << std::setprecision(3) << temp
                  << " Mean Temp: " << std::setw(8) << std::setprecision(3) << mean_temperature
                  << " Kinet Energy: " << std::setw(12) << std::setprecision(8) << kinetic_energy
                  << " Poten Energy: " << std::setw(12) << std::setprecision(8) << PE
                  << " Total Energy: " << std::setw(12) << std::setprecision(8) << energy << std::endl;
#endif
91
    }
92
93
94
#if DEBUG
    std::cout << "Mean Temperature: " << mean_temperature << std::endl;;
#endif
95
96
97
    ASSERT_EQUAL_TOL(temperature, mean_temperature, 0.02);
}

98
void testDimerBox(bool constrain=true) {
99
    // Check conservation of system + bath energy for a harmonic oscillator
100
    System system;
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
    int numMolecules = 20;
    double boxLength = 2; // nm
    Vec3 a(boxLength, 0.0, 0.0);
    Vec3 b(0.0, boxLength, 0.0);
    Vec3 c(0.0, 0.0, boxLength);
    double mass = 20;
    double bondForceConstant = 30000; //0.001;
    double bondLength = 0.1;
    double bondLengthSquared = bondLength * bondLength;
    int numDOF = 0;
    NonbondedForce* forceField = new NonbondedForce();
    HarmonicBondForce* bondForce = new HarmonicBondForce();
    for(int molecule = 0; molecule < numMolecules; ++molecule) {
        int particle1 = system.addParticle(mass);
        int particle2 = system.addParticle(mass);
        forceField->addParticle(0.0, 0.1, 1.0);
        forceField->addParticle(0.0, 0.1, 1.0);
        forceField->addException(particle1, particle2, 0, 0, 0);
        bondForce->addBond(particle1, particle2, bondLength, bondForceConstant);
        numDOF += 6;
        if (constrain) {
            system.addConstraint(particle1, particle2, bondLength);
            numDOF -= 1;
        }
    }
    forceField->setCutoffDistance(.99*boxLength/2);
    forceField->setSwitchingDistance(.88*boxLength/2);
    forceField->setUseSwitchingFunction(true);
    forceField->setUseDispersionCorrection(false);
    forceField->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
    system.addForce(forceField);
    system.addForce(bondForce);
    system.setDefaultPeriodicBoxVectors(a, b, c);
    std::vector<Vec3> positions(numMolecules*2);
    OpenMM_SFMT::SFMT sfmt;
    init_gen_rand(0, sfmt);
    for (int i = 0; i < numMolecules; i++) {
        while (true) {
            Vec3 pos = Vec3(boxLength*genrand_real2(sfmt), boxLength*genrand_real2(sfmt), boxLength*genrand_real2(sfmt));
            Vec3 pos1 = pos + Vec3(0,0, bondLength/2);
            Vec3 pos2 = pos + Vec3(0,0,-bondLength/2);
            double minDist = 2*boxLength;
            for (int j = 0; j < i; j++) {
                Vec3 delta = pos1-positions[j];
                minDist = std::min(minDist, sqrt(delta.dot(delta)));
                delta = pos2-positions[j];
                minDist = std::min(minDist, sqrt(delta.dot(delta)));
            }
            if (minDist > 0.15) {
                positions[2*i+0] = pos1;
                positions[2*i+1] = pos2;
                break;
            }
        }
    }
 
157
158
    VelocityVerletIntegrator integrator(0.001);
    double temperature = 300; // kelvin
159
    double collisionFrequency = 25; // 1/ps
160
161
162
163
164
    int numMTS = 3;
    int numYS = 3;
    int chainLength = 5;
    integrator.addNoseHooverChainThermostat(system, temperature, collisionFrequency,
                                            chainLength, numMTS, numYS);
165
166
    Context context(system, integrator, platform);
    context.setPositions(positions);
167
168
    context.setVelocitiesToTemperature(temperature);
    integrator.step(1500);
169

170
    int nSteps = 2000;
171
    double mean_temp = 0.0;
172
173
174
175
176
177
178
179
180
181
182
183
    std::vector<double> energies(nSteps);
    for (int i = 0; i < nSteps; ++i) {
        integrator.step(1);
        State state = context.getState(State::Energy | (constrain ? State::Positions : 0));
        if (constrain) {
            auto positions = state.getPositions();
            for(int i = 0; i < numMolecules; ++i) {
                Vec3 delta = positions[2*i+1] - positions[2*i];
                double dR2 = delta.dot(delta);
                ASSERT_EQUAL_TOL(bondLengthSquared, dR2, 1e-4);
            }
        }
184
185
        double KE = state.getKineticEnergy();
        double PE = state.getPotentialEnergy();
186
        double time = state.getTime();
187
188
        double instantaneous_temperature = 2 * KE / (BOLTZ * numDOF);
        mean_temp = (i*mean_temp + instantaneous_temperature)/(i+1);
189
        double energy = KE + PE + integrator.computeHeatBathEnergy();
190
191
192
193
194
195
196
#if DEBUG
        if(i % 10 == 0) std::cout << " Time: " << std::setw(4) << time 
                                  << " Temperature: " << std::setw(8) << std::setprecision(3) << temperature
                                  << " Mean Temperature: " << std::setw(8) << std::setprecision(3) << mean_temp
                                  << " Total Energy: " << std::setw(12) << std::setprecision(8) << energy << std::endl;
#endif
        energies[i] = energy;
197
    }
198
199
200
201
202
203
204
205
206
207
    double sum = std::accumulate(energies.begin(), energies.end(), 0.0);
    double mean = sum / energies.size();
    double sq_sum = std::inner_product(energies.begin(), energies.end(), energies.begin(), 0.0);
    double std = std::sqrt(sq_sum / energies.size() - mean * mean);
    double relative_std = std / mean;
    //std::cout << "Relative STD of energy " << relative_std << std::endl;
    // Check mean temperature
    ASSERT_EQUAL_TOL(temperature, temperature, 1e-2);
    // Check fluctuation of conserved (total bath + system) energy
    ASSERT_EQUAL_TOL(relative_std, 0, 1e-4);
208
209
210
}

void testNHCPropagation() {
211
212
213
214
215
216
    // test if the velocity scale factor goes to one for a single particle
    // with no forces in the system
    for (int numMTS = 1; numMTS < 5; numMTS++){
        for (int numYS = 1; numYS <=5; numYS+=2){
            for (int chainLength = 2; chainLength < 6; chainLength += 2){
            double temperature = 300; // kelvin
217
            double collisionFrequency = 10; // 1/ps
218
219
220
221
222
223
224
225
226
            VelocityVerletIntegrator integrator(0.001);
            // make system
            System system;
            double mass = 1;
            system.addParticle(mass);
            int chainID = integrator.addNoseHooverChainThermostat(system, temperature, collisionFrequency,
                                                                  chainLength, numMTS, numYS);
            Context context(system, integrator, platform);
            // propagate the chain
227
228
229
230
            double velocity = 1;
            double temp, scale, kineticEnergy;
            double mean_temp=0, mean_scale=0;
            for(int i = 0; i < 10000; ++i) {
231
232
233
                kineticEnergy = 3 * 0.5 * mass * velocity * velocity;
                scale = integrator.propagateChain(kineticEnergy, chainID);
                velocity *= scale;
234
235
236
                temp = 2* kineticEnergy / BOLTZ / 3;
                mean_temp = (i*mean_temp + temp)/(i+1);
                mean_scale = (i*mean_scale + scale)/(i+1);
237
            }
238
239
240
            //std::cout << mean_scale <<  " " << mean_temp << std::endl;
            ASSERT_EQUAL_TOL(1, mean_scale,  1e-2);
            ASSERT_EQUAL_TOL(temperature, mean_temp,  0.25);
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
            }
        }
    }
}

void testPropagateChainConsistentWithPythonReference() {
    VelocityVerletIntegrator integrator(0.001);
    System system;
    double mass = 1;
    system.addParticle(mass);
    double kineticEnergy = 1e6;
    double temperature=300, collisionFrequency=1, chainLength=3, numMTS=3, numYS=3;
    int chainID = integrator.addNoseHooverChainThermostat(system, temperature, collisionFrequency,
                                                                  chainLength, numMTS, numYS);
    Context context(system, integrator, platform);
    double scale = integrator.propagateChain(kineticEnergy, chainID);
257
258
259
#if DEBUG
    std::cout << std::setw(12) << std::setprecision(10) << scale << std::endl;
#endif
260
    ASSERT_EQUAL_TOL(0.9674732261005896, scale, 1e-5)
261
262
263
264
265
266
}

void runPlatformTests();

int main(int argc, char* argv[]) {
    try {
267
        initializeTests(argc, argv);
268
        testHarmonicOscillator();
269
270
271
272
        bool constrain;
        constrain = false; testDimerBox(constrain);
        constrain = true; testDimerBox(constrain);
        runPlatformTests();
273
274
275
276
277
278
279
280
    }
    catch(const exception& e) {
        cout << "exception: " << e.what() << endl;
        return 1;
    }
    cout << "Done" << endl;
    return 0;
}