NoseHooverIntegrator.cpp 21 KB
Newer Older
1
2
3
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
Evan Pretti's avatar
Evan Pretti committed
4
5
 * This is part of the OpenMM molecular simulation toolkit.                   *
 * See https://openmm.org/development.                                        *
6
 *                                                                            *
7
 * Portions copyright (c) 2019-2024 Stanford University and the Authors.      *
8
 * Authors: Andreas Krämer and Andrew C. Simmonett                            *
9
 * Contributors: Peter Eastman                                                *
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 *                                                                            *
 * 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.                                     *
 * -------------------------------------------------------------------------- */

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

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

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

61
NoseHooverIntegrator::~NoseHooverIntegrator() {}
62

63
int NoseHooverIntegrator::addThermostat(double temperature, double collisionFrequency,
64
                                        int chainLength, int numMTS, int numYoshidaSuzuki) {
65
66
67
68
    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
69
    if (chainLength <= 0)
70
71
72
73
74
        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");
75
    hasSubsystemThermostats_ = false;
76
77
    return addSubsystemThermostat(std::vector<int>(), std::vector<std::pair<int, int>>(), temperature,
                                  collisionFrequency, temperature, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
78
79
}

80
int NoseHooverIntegrator::addSubsystemThermostat(const std::vector<int>& thermostatedParticles,
81
                                                 const std::vector< std::pair< int, int> > &thermostatedPairs,
82
83
                                                 double temperature, double collisionFrequency,
                                                 double relativeTemperature, double relativeCollisionFrequency,
84
                                                 int chainLength, int numMTS, int numYoshidaSuzuki) {
85
86
87
88
89
90
91
92
    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
93
    if (chainLength <= 0)
94
95
96
97
98
        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");
99
    int chainID = noseHooverChains.size();
100
    // check if one thermostat already applies to all atoms or pairs
101
    if ( (chainID > 0) && (noseHooverChains[0].getThermostatedAtoms().size()+noseHooverChains[0].getThermostatedPairs().size() == 0) ) {
102
103
104
105
106
107
        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."
        );
    }
108
109
110
111
112
113
114
    int nDOF = 0; // is set in initializeThermostats()
    noseHooverChains.emplace_back(temperature, relativeTemperature,
                                 collisionFrequency, relativeCollisionFrequency,
                                 nDOF, chainLength, numMTS,
                                 numYoshidaSuzuki, chainID,
                                 thermostatedParticles, thermostatedPairs);
    return chainID;
115
116
}

117
const NoseHooverChain& NoseHooverIntegrator::getThermostat(int chainID) const {
118
119
120
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
    return noseHooverChains[chainID];
}
121

122
void NoseHooverIntegrator::initializeThermostats(const System &system) {
123

124
125
126
127
128
    allAtoms.clear();
    allPairs.clear();
    std::set<int> allAtomsSet;
    for (int atom = 0; atom < system.getNumParticles(); ++atom) allAtomsSet.insert(atom);

129
130
131
    for (auto &thermostat : noseHooverChains) {
        const auto &thermostatedParticles = thermostat.getThermostatedAtoms();
        const auto &thermostatedPairs = thermostat.getThermostatedPairs();
132

133
134
135
        for( auto& pair : thermostatedPairs ) {
            allAtomsSet.erase(pair.first);
            allAtomsSet.erase(pair.second);
136
            allPairs.emplace_back(pair.first, pair.second, thermostat.getRelativeTemperature());
137
138
        }

139
        // figure out the number of DOFs
140
        int nDOF = 3*(thermostatedParticles.size() + thermostatedPairs.size());
141
142
143
144
        for (int constraintNum = 0; constraintNum < system.getNumConstraints(); constraintNum++) {
            int particle1, particle2;
            double distance;
            system.getConstraintParameters(constraintNum, particle1, particle2, distance);
145
146
147
148
149
            bool particle1_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(),
                                                                  thermostatedParticles.end(), particle1)
                                                                    != thermostatedParticles.end())) ||
                                                      (std::find_if(thermostatedPairs.begin(),
                                                                    thermostatedPairs.end(),
150
151
                                                                    [&particle1](const std::pair<int, int>& pair){
                                                                           return pair.first == particle1 || pair.second == particle1;})
152
153
154
155
156
157
                                                                      != thermostatedPairs.end());
            bool particle2_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(),
                                                                  thermostatedParticles.end(), particle2)
                                                                    != thermostatedParticles.end())) ||
                                                      (std::find_if(thermostatedPairs.begin(),
                                                                    thermostatedPairs.end(),
158
159
                                                                    [&particle2](const std::pair<int, int>& pair){
                                                                           return pair.first == particle2 || pair.second == particle2;})
160
                                                                      != thermostatedPairs.end());
161
162
163
164
165
166
167
168
169
            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;
                }
170
            }
171
172
        }

173
        // set number of DoFs for chain 
174
        thermostat.setNumDegreesOfFreedom(nDOF);
175
    }
176

177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
    // If there is a CMMotionRemover, remove 3 degrees of freedom from the largest thermostat.
    for (int force = 0; force < system.getNumForces(); ++force) {
        if (dynamic_cast<const CMMotionRemover*>(&system.getForce(force))) {
            int largest = 0;
            int largestSize = noseHooverChains[0].getThermostatedAtoms().size() + noseHooverChains[0].getThermostatedPairs().size();
            for (int chain = 1; chain < noseHooverChains.size(); chain++) {
                int size = noseHooverChains[chain].getThermostatedAtoms().size() + noseHooverChains[chain].getThermostatedPairs().size();
                if (size > largestSize) {
                    largest = chain;
                    largestSize = size;
                }
            }
            noseHooverChains[largest].setNumDegreesOfFreedom(noseHooverChains[largest].getNumDegreesOfFreedom()-3);
            break;
        }
    }
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
    double kE = 0.0;
290
    if(noseHooverChains.size() > 0) {
291
        for (const auto &nhc: noseHooverChains){
292
            kE += kernel.getAs<IntegrateNoseHooverStepKernel>().computeMaskedKineticEnergy(*context, nhc, true).first;
293
294
        }
    } else {
295
        kE = kernel.getAs<IntegrateNoseHooverStepKernel>().computeKineticEnergy(*context, *this);
296
297
    }
    return kE;
298
299
}

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

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

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

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

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

343
    initializeThermostats(system);
344
345
}

346
void NoseHooverIntegrator::cleanup() {
347
    kernel = Kernel();
348
349
}

350
vector<string> NoseHooverIntegrator::getKernelNames() {
Peter Eastman's avatar
Peter Eastman committed
351
    return {IntegrateNoseHooverStepKernel::Name()};
352
353
}

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

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);
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
}

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);
}