NoseHooverIntegrator.cpp 21.4 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-2023 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
        // set number of DoFs for chain 
182
        thermostat.setNumDegreesOfFreedom(nDOF);
183
    }
184

185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
    // 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;
        }
    }
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
232
233
234
235
236
237
238
239

    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.");
                }
            }
240
        }
241

242
243
244
        // make sure that massless particles are not thermostated
        for(auto particle: nhc.getThermostatedAtoms()){
            double mass = system.getParticleMass(particle);
245
            if (mass == 0.0) {
246
                throw OpenMMException("Found a particle with no mass (" + std::to_string(particle) + ") in a thermostat. Massless particles cannot be thermostated.");
247
248
249
            }
        }
    }
250
251

    for(const auto& atom : allAtomsSet) allAtoms.push_back(atom);
252
253
}

254
double NoseHooverIntegrator::getTemperature(int chainID) const {
255
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
256
    return noseHooverChains[chainID].getTemperature();
257
258
}

259
void NoseHooverIntegrator::setTemperature(double temperature, int chainID){
260
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
261
    noseHooverChains[chainID].setTemperature(temperature);
262
263
264

}

265
double NoseHooverIntegrator::getRelativeTemperature(int chainID) const {
266
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
267
    return noseHooverChains[chainID].getRelativeTemperature();
268
269
270
}

void NoseHooverIntegrator::setRelativeTemperature(double temperature, int chainID){
271
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
272
    noseHooverChains[chainID].setRelativeTemperature(temperature);
273
274
275
276

}

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

void NoseHooverIntegrator::setCollisionFrequency(double frequency, int chainID){
282
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
283
    noseHooverChains[chainID].setCollisionFrequency(frequency);
284
285
}

286
double NoseHooverIntegrator::getRelativeCollisionFrequency(int chainID) const {
287
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
288
    return noseHooverChains[chainID].getRelativeCollisionFrequency();
289
290
291
}

void NoseHooverIntegrator::setRelativeCollisionFrequency(double frequency, int chainID){
292
    ASSERT_VALID_INDEX(chainID, noseHooverChains);
293
    noseHooverChains[chainID].setRelativeCollisionFrequency(frequency);
294
295
296
}

double NoseHooverIntegrator::computeKineticEnergy() {
297
    forcesAreValid = false;
298
    double kE = 0.0;
299
    if(noseHooverChains.size() > 0) {
300
        for (const auto &nhc: noseHooverChains){
301
            kE += kernel.getAs<IntegrateNoseHooverStepKernel>().computeMaskedKineticEnergy(*context, nhc, true).first;
302
303
        }
    } else {
304
        kE = kernel.getAs<IntegrateNoseHooverStepKernel>().computeKineticEnergy(*context, *this);
305
306
    }
    return kE;
307
308
}

309
310
311
312
bool NoseHooverIntegrator::kineticEnergyRequiresForce() const {
    return false;
}

313
double NoseHooverIntegrator::computeHeatBathEnergy() {
314
315
    double energy = 0;
    for(auto &nhc : noseHooverChains) {
316
        if (context && (nhc.getNumDegreesOfFreedom() > 0)) {
317
            energy += kernel.getAs<IntegrateNoseHooverStepKernel>().computeHeatBathEnergy(*context, nhc);
Andy Simmonett's avatar
Andy Simmonett committed
318
        }
319
320
321
322
    }
    return energy;
}

323
void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
324
325
    if (owner != NULL && &contextRef.getOwner() != owner)
        throw OpenMMException("This Integrator is already bound to a context");
326

327
    context = &contextRef;
328
    const System& system = context->getSystem();
329
    owner = &contextRef.getOwner();
330
331
    kernel = context->getPlatform().createKernel(IntegrateNoseHooverStepKernel::Name(), contextRef);
    kernel.getAs<IntegrateNoseHooverStepKernel>().initialize(contextRef.getSystem(), *this);
332
    forcesAreValid = false;
333
334

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

353
    initializeThermostats(system);
354
355
}

356
void NoseHooverIntegrator::cleanup() {
357
    kernel = Kernel();
358
359
}

360
vector<string> NoseHooverIntegrator::getKernelNames() {
361
    std::vector<std::string> names;
362
    names.push_back(IntegrateNoseHooverStepKernel::Name());
363
364
365
    return names;
}

366
void NoseHooverIntegrator::step(int steps) {
367
368
369
    if (context == NULL)
        throw OpenMMException("This Integrator is not bound to a context!");
    for (int i = 0; i < steps; ++i) {
370
371
        if(context->updateContextState())
            forcesAreValid = false;
372
        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
373
        kernel.getAs<IntegrateNoseHooverStepKernel>().execute(*context, *this, forcesAreValid);
374
375
    }
}
376
377
378
379
380
381
382

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);
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
}

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