NoseHooverIntegrator.cpp 20.8 KB
Newer Older
1
2
3
4
5
6
7
8
/* -------------------------------------------------------------------------- *
 *                                   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.               *
 *                                                                            *
9
 * Portions copyright (c) 2019-2020 Stanford University and the Authors.      *
10
 * Authors: Andreas Krämer and Andrew C. Simmonett                            *
11
 * Contributors: Peter Eastman                                                *
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
 *                                                                            *
 * 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.                                     *
 * -------------------------------------------------------------------------- */

32
#include "openmm/NoseHooverIntegrator.h"
33
34
35
#include "openmm/Context.h"
#include "openmm/Force.h"
#include "openmm/System.h"
36
#include "openmm/NoseHooverChain.h"
37
38
39
#include "openmm/OpenMMException.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/internal/ContextImpl.h"
40
#include "openmm/internal/AssertionUtilities.h"
41
#include "openmm/kernels.h"
42
#include <iostream>
43
#include <string>
44
#include <algorithm>
45
46
47
48
49

using namespace OpenMM;
using std::string;
using std::vector;

50
NoseHooverIntegrator::NoseHooverIntegrator(double stepSize):
51
52
    forcesAreValid(false),
    hasSubsystemThermostats_(true)
53
{
54
55
    setStepSize(stepSize);
    setConstraintTolerance(1e-5);
Andy Simmonett's avatar
Andy Simmonett committed
56
    setMaximumPairDistance(0.0);
57
}
58
NoseHooverIntegrator::NoseHooverIntegrator(double temperature, double collisionFrequency, double stepSize,
59
60
61
62
                                           int chainLength, int numMTS, int numYoshidaSuzuki) : 
    forcesAreValid(false),
    hasSubsystemThermostats_(false) {

63
64
    setStepSize(stepSize);
    setConstraintTolerance(1e-5);
Andy Simmonett's avatar
Andy Simmonett committed
65
    setMaximumPairDistance(0.0);
66
    addThermostat(temperature, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
67
}
68

69
NoseHooverIntegrator::~NoseHooverIntegrator() {}
70

71
int NoseHooverIntegrator::addThermostat(double temperature, double collisionFrequency,
72
                                        int chainLength, int numMTS, int numYoshidaSuzuki) {
73
74
75
76
    if (temperature < 0)
        throw OpenMMException("NoseHooverIntegrator: temperature cannot be negative");
    if (collisionFrequency <= 0)
        throw OpenMMException("NoseHooverIntegrator: collisionFrequency must be positive");
oxdc's avatar
oxdc committed
77
    if (chainLength <= 0)
78
79
80
81
82
        throw OpenMMException("NoseHooverIntegrator: chainLength must be positive");
    if (numMTS < 0)
        throw OpenMMException("NoseHooverIntegrator: numMTS must be positive");
    if (numYoshidaSuzuki != 1 && numYoshidaSuzuki != 3 && numYoshidaSuzuki != 5 && numYoshidaSuzuki != 7)
        throw OpenMMException("NoseHooverIntegrator: numYoshidaSuzuki must be 1, 3, 5, or 7");
83
    hasSubsystemThermostats_ = false;
84
85
    return addSubsystemThermostat(std::vector<int>(), std::vector<std::pair<int, int>>(), temperature,
                                  collisionFrequency, temperature, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
86
87
}

88
int NoseHooverIntegrator::addSubsystemThermostat(const std::vector<int>& thermostatedParticles,
89
                                                 const std::vector< std::pair< int, int> > &thermostatedPairs,
90
91
                                                 double temperature, double collisionFrequency,
                                                 double relativeTemperature, double relativeCollisionFrequency,
92
                                                 int chainLength, int numMTS, int numYoshidaSuzuki) {
93
94
95
96
97
98
99
100
    if (temperature < 0)
        throw OpenMMException("NoseHooverIntegrator: temperature cannot be negative");
    if (relativeTemperature < 0)
        throw OpenMMException("NoseHooverIntegrator: relativeTemperature cannot be negative");
    if (collisionFrequency <= 0)
        throw OpenMMException("NoseHooverIntegrator: collisionFrequency must be positive");
    if (relativeCollisionFrequency <= 0)
        throw OpenMMException("NoseHooverIntegrator: relativeCollisionFrequency must be positive");
oxdc's avatar
oxdc committed
101
    if (chainLength <= 0)
102
103
104
105
106
        throw OpenMMException("NoseHooverIntegrator: chainLength must be positive");
    if (numMTS < 0)
        throw OpenMMException("NoseHooverIntegrator: numMTS must be positive");
    if (numYoshidaSuzuki != 1 && numYoshidaSuzuki != 3 && numYoshidaSuzuki != 5 && numYoshidaSuzuki != 7)
        throw OpenMMException("NoseHooverIntegrator: numYoshidaSuzuki must be 1, 3, 5, or 7");
107
    int chainID = noseHooverChains.size();
108
    // check if one thermostat already applies to all atoms or pairs
109
    if ( (chainID > 0) && (noseHooverChains[0].getThermostatedAtoms().size()+noseHooverChains[0].getThermostatedPairs().size() == 0) ) {
110
111
112
113
114
115
        throw OpenMMException(
            "Cannot add thermostat, since one of the thermostats already in the integrator applies to all particles. "
            "To manually add thermostats, use the constructor that takes only the "
            "step size as an input argument."
        );
    }
116
117
118
119
120
121
122
    int nDOF = 0; // is set in initializeThermostats()
    noseHooverChains.emplace_back(temperature, relativeTemperature,
                                 collisionFrequency, relativeCollisionFrequency,
                                 nDOF, chainLength, numMTS,
                                 numYoshidaSuzuki, chainID,
                                 thermostatedParticles, thermostatedPairs);
    return chainID;
123
124
}

125
const NoseHooverChain& NoseHooverIntegrator::getThermostat(int chainID) const {
126
127
128
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
    return noseHooverChains[chainID];
}
129

130
void NoseHooverIntegrator::initializeThermostats(const System &system) {
131

132
133
134
135
136
    allAtoms.clear();
    allPairs.clear();
    std::set<int> allAtomsSet;
    for (int atom = 0; atom < system.getNumParticles(); ++atom) allAtomsSet.insert(atom);

137
138
139
    for (auto &thermostat : noseHooverChains) {
        const auto &thermostatedParticles = thermostat.getThermostatedAtoms();
        const auto &thermostatedPairs = thermostat.getThermostatedPairs();
140

141
142
143
        for( auto& pair : thermostatedPairs ) {
            allAtomsSet.erase(pair.first);
            allAtomsSet.erase(pair.second);
144
            allPairs.emplace_back(pair.first, pair.second, thermostat.getRelativeTemperature());
145
146
        }

147
        // figure out the number of DOFs
148
        int nDOF = 3*(thermostatedParticles.size() + thermostatedPairs.size());
149
150
151
152
        for (int constraintNum = 0; constraintNum < system.getNumConstraints(); constraintNum++) {
            int particle1, particle2;
            double distance;
            system.getConstraintParameters(constraintNum, particle1, particle2, distance);
153
154
155
156
157
            bool particle1_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(),
                                                                  thermostatedParticles.end(), particle1)
                                                                    != thermostatedParticles.end())) ||
                                                      (std::find_if(thermostatedPairs.begin(),
                                                                    thermostatedPairs.end(),
158
159
                                                                    [&particle1](const std::pair<int, int>& pair){
                                                                           return pair.first == particle1 || pair.second == particle1;})
160
161
162
163
164
165
                                                                      != thermostatedPairs.end());
            bool particle2_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(),
                                                                  thermostatedParticles.end(), particle2)
                                                                    != thermostatedParticles.end())) ||
                                                      (std::find_if(thermostatedPairs.begin(),
                                                                    thermostatedPairs.end(),
166
167
                                                                    [&particle2](const std::pair<int, int>& pair){
                                                                           return pair.first == particle2 || pair.second == particle2;})
168
                                                                      != thermostatedPairs.end());
169
170
171
172
173
174
175
176
177
            if ((system.getParticleMass(particle1) > 0) && (system.getParticleMass(particle2) > 0)){
                if ((particle1_in_thermostatedParticles && !particle2_in_thermostatedParticles) ||
                     (!particle1_in_thermostatedParticles && particle2_in_thermostatedParticles)){
                    throw OpenMMException("Cannot add only one of particles " + std::to_string(particle1) + " and " + std::to_string(particle2)
                                            + " to NoseHooverChain, because they are connected by a constraint.");
                }
                if (particle1_in_thermostatedParticles && particle2_in_thermostatedParticles){
                    nDOF -= 1;
                }
178
            }
179
180
181
182
        }

        // remove 3 degrees of freedom from thermostats that act on absolute motions
        int numForces = system.getNumForces();
183
        if (thermostatedPairs.size() == 0){
184
185
            for (int forceNum = 0; forceNum < numForces; ++forceNum) {
                if (dynamic_cast<const CMMotionRemover*>(&system.getForce(forceNum))) nDOF -= 3;
186
187
            }
        }
188

189
        // set number of DoFs for chain 
190
        thermostat.setNumDegreesOfFreedom(nDOF);
191
    }
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231


    for (int chain1 = 0; chain1 < noseHooverChains.size(); ++chain1){
        const auto& nhc = noseHooverChains[chain1];

       // make sure that thermostats do not overlap
        for (int chain2 = 0; chain2 < chain1; ++chain2){
            const auto& other_nhc = noseHooverChains[chain2];
            for (auto &particle: nhc.getThermostatedAtoms()){
                bool isParticleInOtherChain = (std::find(other_nhc.getThermostatedAtoms().begin(),
                                                         other_nhc.getThermostatedAtoms().end(),
                                                         particle) != other_nhc.getThermostatedAtoms().end()) ||
                                              (std::find_if(other_nhc.getThermostatedPairs().begin(),
                                                            other_nhc.getThermostatedPairs().end(),
                                               [&particle](const std::pair<int, int>& pair){ return pair.first == particle || pair.second == particle;})
                                                  != other_nhc.getThermostatedPairs().end());
                if (isParticleInOtherChain){
                    throw OpenMMException("Found particle " + std::to_string(particle) + "in a different NoseHooverChain, "
                                          "but particles can only be thermostated by one thermostat.");
                }
            }
            for (auto &pair: nhc.getThermostatedPairs()){
                bool isParticleInOtherChain = (std::find(other_nhc.getThermostatedAtoms().begin(),
                                                         other_nhc.getThermostatedAtoms().end(),
                                                         pair.first) != other_nhc.getThermostatedAtoms().end()) ||
                                              (std::find(other_nhc.getThermostatedAtoms().begin(),
                                                         other_nhc.getThermostatedAtoms().end(),
                                                         pair.second) != other_nhc.getThermostatedAtoms().end()) ||
                                              (std::find_if(other_nhc.getThermostatedPairs().begin(),
                                                            other_nhc.getThermostatedPairs().end(),
                                               [&pair](const std::pair<int, int>& other_pair){
                                                        return pair.first == other_pair.first || pair.first == other_pair.second ||
                                                               pair.second == other_pair.first || pair.second == other_pair.second;})
                                                  != other_nhc.getThermostatedPairs().end());
                if (isParticleInOtherChain){
                    throw OpenMMException("Found pair " + std::to_string(pair.first) + "," +
                                          std::to_string(pair.second) + " in a different NoseHooverChain, "
                                          "but particles can only be thermostated by one thermostat.");
                }
            }
232
        }
233

234
235
236
        // make sure that massless particles are not thermostated
        for(auto particle: nhc.getThermostatedAtoms()){
            double mass = system.getParticleMass(particle);
237
            if (mass == 0.0) {
238
                throw OpenMMException("Found a particle with no mass (" + std::to_string(particle) + ") in a thermostat. Massless particles cannot be thermostated.");
239
240
241
            }
        }
    }
242
243

    for(const auto& atom : allAtomsSet) allAtoms.push_back(atom);
244
245
}

246
double NoseHooverIntegrator::getTemperature(int chainID) const {
247
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
248
    return noseHooverChains[chainID].getTemperature();
249
250
}

251
void NoseHooverIntegrator::setTemperature(double temperature, int chainID){
252
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
253
    noseHooverChains[chainID].setTemperature(temperature);
254
255
256

}

257
double NoseHooverIntegrator::getRelativeTemperature(int chainID) const {
258
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
259
    return noseHooverChains[chainID].getRelativeTemperature();
260
261
262
}

void NoseHooverIntegrator::setRelativeTemperature(double temperature, int chainID){
263
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
264
    noseHooverChains[chainID].setRelativeTemperature(temperature);
265
266
267
268

}

double NoseHooverIntegrator::getCollisionFrequency(int chainID) const {
269
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
270
    return noseHooverChains[chainID].getCollisionFrequency();
271
272
273
}

void NoseHooverIntegrator::setCollisionFrequency(double frequency, int chainID){
274
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
275
    noseHooverChains[chainID].setCollisionFrequency(frequency);
276
277
}

278
double NoseHooverIntegrator::getRelativeCollisionFrequency(int chainID) const {
279
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
280
    return noseHooverChains[chainID].getRelativeCollisionFrequency();
281
282
283
}

void NoseHooverIntegrator::setRelativeCollisionFrequency(double frequency, int chainID){
284
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
285
    noseHooverChains[chainID].setRelativeCollisionFrequency(frequency);
286
287
288
}

double NoseHooverIntegrator::computeKineticEnergy() {
289
    forcesAreValid = false;
290
    double kE = 0.0;
291
    if(noseHooverChains.size() > 0) {
292
        for (const auto &nhc: noseHooverChains){
293
            kE += kernel.getAs<IntegrateNoseHooverStepKernel>().computeMaskedKineticEnergy(*context, nhc, true).first;
294
295
        }
    } else {
296
        kE = kernel.getAs<IntegrateNoseHooverStepKernel>().computeKineticEnergy(*context, *this);
297
298
    }
    return kE;
299
300
}

301
302
303
304
bool NoseHooverIntegrator::kineticEnergyRequiresForce() const {
    return false;
}

305
double NoseHooverIntegrator::computeHeatBathEnergy() {
306
307
    double energy = 0;
    for(auto &nhc : noseHooverChains) {
308
        if (context && (nhc.getNumDegreesOfFreedom() > 0)) {
309
            energy += kernel.getAs<IntegrateNoseHooverStepKernel>().computeHeatBathEnergy(*context, nhc);
Andy Simmonett's avatar
Andy Simmonett committed
310
        }
311
312
313
314
    }
    return energy;
}

315
void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
316
317
    if (owner != NULL && &contextRef.getOwner() != owner)
        throw OpenMMException("This Integrator is already bound to a context");
318

319
    context = &contextRef;
320
    const System& system = context->getSystem();
321
    owner = &contextRef.getOwner();
322
323
    kernel = context->getPlatform().createKernel(IntegrateNoseHooverStepKernel::Name(), contextRef);
    kernel.getAs<IntegrateNoseHooverStepKernel>().initialize(contextRef.getSystem(), *this);
324
    forcesAreValid = false;
325
326

    // check for drude particles and build the Nose-Hoover Chains
327
    for (auto& thermostat: noseHooverChains){
328
        // if there are no thermostated particles or pairs in the lists this is a regular thermostat for the whole (non-Drude) system
329
330
        if ( (thermostat.getThermostatedAtoms().size() == 0) && (thermostat.getThermostatedPairs().size() == 0) ){
            std::vector<int> thermostatedParticles;
331
332
333
            for(int particle = 0; particle < system.getNumParticles(); ++particle) {
                double mass = system.getParticleMass(particle);
                if ( (mass > 0) && (mass < 0.8) ){
334
                    std::cout << "Warning: Found particles with mass between 0.0 and 0.8 dalton. Did you mean to make a DrudeNoseHooverIntegrator instead? "
335
336
337
                                 "The thermostat you are about to use will not treat these particles as Drude particles!" << std::endl;
                }
                if(system.getParticleMass(particle) > 0) {
338
                   thermostatedParticles.push_back(particle);
339
340
                }
            }
341
            thermostat.setThermostatedAtoms(thermostatedParticles);
342
343
344
        }
    }

345
    initializeThermostats(system);
346
347
}

348
void NoseHooverIntegrator::cleanup() {
349
    kernel = Kernel();
350
351
}

352
vector<string> NoseHooverIntegrator::getKernelNames() {
353
    std::vector<std::string> names;
354
    names.push_back(IntegrateNoseHooverStepKernel::Name());
355
356
357
    return names;
}

358
void NoseHooverIntegrator::step(int steps) {
359
360
361
    if (context == NULL)
        throw OpenMMException("This Integrator is not bound to a context!");
    for (int i = 0; i < steps; ++i) {
362
363
        if(context->updateContextState())
            forcesAreValid = false;
364
        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
365
        kernel.getAs<IntegrateNoseHooverStepKernel>().execute(*context, *this, forcesAreValid);
366
367
    }
}
368
369
370
371
372
373
374

void NoseHooverIntegrator::createCheckpoint(std::ostream& stream) const {
    kernel.getAs<IntegrateNoseHooverStepKernel>().createCheckpoint(*context, stream);
}

void NoseHooverIntegrator::loadCheckpoint(std::istream& stream) {
    kernel.getAs<IntegrateNoseHooverStepKernel>().loadCheckpoint(*context, stream);
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
}

void NoseHooverIntegrator::serializeParameters(SerializationNode& node) const {
    node.setIntProperty("version", 1);
    vector<vector<double> > positions, velocities;
    kernel.getAs<IntegrateNoseHooverStepKernel>().getChainStates(*context, positions, velocities);
    for (int i = 0; i < positions.size(); i++) {
        SerializationNode& chain = node.createChildNode("Chain");
        for (int j = 0; j < positions[i].size(); j++)
            chain.createChildNode("Bead").setDoubleProperty("position", positions[i][j]).setDoubleProperty("velocity", velocities[i][j]);
    }
}

void NoseHooverIntegrator::deserializeParameters(const SerializationNode& node) {
    if (node.getIntProperty("version") != 1)
        throw OpenMMException("Unsupported version number");
    int numChains = node.getChildren().size();
    vector<vector<double> > positions(numChains), velocities(numChains);
    for (int i = 0; i < numChains; i++) {
        auto& chain = node.getChildren()[i];
        for (auto& bead : chain.getChildren()) {
            positions[i].push_back(bead.getDoubleProperty("position"));
            velocities[i].push_back(bead.getDoubleProperty("velocity"));
        }
    }
    kernel.getAs<IntegrateNoseHooverStepKernel>().setChainStates(*context, positions, velocities);
}