NoseHooverIntegrator.cpp 15.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
/* -------------------------------------------------------------------------- *
 *                                   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.                                     *
 * -------------------------------------------------------------------------- */

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

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

49
NoseHooverIntegrator::NoseHooverIntegrator(double stepSize):
50
51
    forcesAreValid(false)
{
52
53
54
    setStepSize(stepSize);
    setConstraintTolerance(1e-5);
}
55
56
57
58
59
60
NoseHooverIntegrator::NoseHooverIntegrator(double stepSize, System &system, double temperature, int collisionFrequnency,
                                           int chainLength, int numMTS, int numYoshidaSuzuki) : forcesAreValid(false) {
    setStepSize(stepSize);
    setConstraintTolerance(1e-5);
    addThermostat(system, temperature, collisionFrequnency, chainLength, numMTS, numYoshidaSuzuki);
}
61

62
NoseHooverIntegrator::~NoseHooverIntegrator() {}
63

64
std::pair<double, double> NoseHooverIntegrator::propagateChain(std::pair<double, double> kineticEnergy, int chainID) {
65
    return nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, noseHooverChains.at(chainID), kineticEnergy, getStepSize());
66
67
}

68
int NoseHooverIntegrator::addThermostat(System& system, double temperature, double collisionFrequency,
69
                                                           int chainLength, int numMTS, int numYoshidaSuzuki) {
Andy Simmonett's avatar
Andy Simmonett committed
70
    std::vector<int> thermostatedParticles;
71
    for(int particle = 0; particle < system.getNumParticles(); ++particle) {
72
73
74
75
76
        double mass = system.getParticleMass(particle);
        if ( (mass > 0) && (mass < 0.8) ){
            std::cout << "Warning: Found particles with mass between 0.0 and 0.8 dalton. Did you mean to make a DrudeVelocityVerletIntegrator instead? "
                         "The thermostat you are about to use will not treat these particles as Drude particles!" << std::endl;
        }
77
78
79
        if(system.getParticleMass(particle) > 0) {
            thermostatedParticles.push_back(particle);
        }
80
    }
81

82
83
    return addSubsystemThermostat(system, thermostatedParticles, std::vector<std::pair<int, int>>(), temperature, temperature,
                                  collisionFrequency, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
84
85
}

86
87
88
89
90
int NoseHooverIntegrator::addSubsystemThermostat(System& system, const std::vector<int>& thermostatedParticles,
                                                 const std::vector< std::pair< int, int> > &thermostatedPairs,
                                                 double temperature, double relativeTemperature,
                                                 double collisionFrequency, double relativeCollisionFrequency,
                                                 int chainLength, int numMTS, int numYoshidaSuzuki) {
91
92
93
94
    if (context) {
        throw OpenMMException("Nose-Hoover chains cannot be added after binding this integrator to a context.");
    }
    // figure out the number of DOFs
95
    int nDOF = 3*(thermostatedParticles.size() + thermostatedPairs.size());
96
97
98
99
    for (int constraintNum = 0; constraintNum < system.getNumConstraints(); constraintNum++) {
        int particle1, particle2;
        double distance;
        system.getConstraintParameters(constraintNum, particle1, particle2, distance);
100
101
102
103
104
105
106
107
108
109
        bool particle1_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(), thermostatedParticles.end(), particle1)
                                                    != thermostatedParticles.end())) ||
                                                  (std::find_if(thermostatedPairs.begin(), thermostatedPairs.end(),
                                                   [&particle1](const std::pair<int, int>& pair){ return pair.first == particle1 || pair.second == particle1;})
                                                      != thermostatedPairs.end());
        bool particle2_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(), thermostatedParticles.end(), particle2)
                                                    != thermostatedParticles.end())) ||
                                                  (std::find_if(thermostatedPairs.begin(), thermostatedPairs.end(),
                                                   [&particle2](const std::pair<int, int>& pair){ return pair.first == particle2 || pair.second == particle2;})
                                                      != thermostatedPairs.end());
110
        if ((system.getParticleMass(particle1) > 0) && (system.getParticleMass(particle2) > 0)){
111
112
            if ((particle1_in_thermostatedParticles && !particle2_in_thermostatedParticles) ||
                 (!particle1_in_thermostatedParticles && particle2_in_thermostatedParticles)){
113
114
115
                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.");
            }
116
            if (particle1_in_thermostatedParticles && particle2_in_thermostatedParticles){
117
118
119
120
121
                nDOF -= 1;
            }
        }
    }
    int numForces = system.getNumForces();
122
    if (thermostatedPairs.size() == 0){ // remove 3 degrees of freedom from thermostats that act on absolute motions
123
        for (int forceNum = 0; forceNum < numForces; ++forceNum) {
124
            if (dynamic_cast<CMMotionRemover*>(&system.getForce(forceNum))) nDOF -= 3;
125
        }
126
127
    }

128
    // make sure that thermostats do not overlap
129
    for (auto &other_nhc: noseHooverChains) {
130
131
132
133
134
135
136
137
        for (auto &particle: thermostatedParticles){
            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());
138
139
140
141
142
143
            if (isParticleInOtherChain){
                throw OpenMMException("Found particle " + std::to_string(particle) + "in a different NoseHooverChain, "
                                      "but particles can only be thermostated by one thermostat.");
            }
        }
    }
144

145
    // create and add new chain
146
    int chainID = noseHooverChains.size();
147
148
    auto nhc = NoseHooverChain(temperature, relativeTemperature, collisionFrequency, relativeCollisionFrequency,
                               nDOF, chainLength, numMTS, numYoshidaSuzuki, chainID, thermostatedParticles, thermostatedPairs);
149
    noseHooverChains.push_back(nhc);
150
151
152
    return chainID;
}

153
double NoseHooverIntegrator::getTemperature(int chainID) const {
154
155
156
157
158
    if (chainID >= noseHooverChains.size()) {
        throw OpenMMException("Cannot get temperature for chainID " + std::to_string(chainID)
                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
        );
    }
159
    return noseHooverChains.at(chainID).getDefaultTemperature();
160
161
}

162
void NoseHooverIntegrator::setTemperature(double temperature, int chainID){
163
164
165
166
167
    if (chainID >= noseHooverChains.size()) {
        throw OpenMMException("Cannot set temperature for chainID " + std::to_string(chainID)
                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
        );
    }
168
    noseHooverChains.at(chainID).setDefaultTemperature(temperature);
169
170
171

}

172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
double NoseHooverIntegrator::getRelativeTemperature(int chainID) const {
    if (chainID >= noseHooverChains.size()) {
        throw OpenMMException("Cannot get relative temperature for chainID " + std::to_string(chainID)
                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
        );
    }
    return noseHooverChains.at(chainID).getDefaultRelativeTemperature();
}

void NoseHooverIntegrator::setRelativeTemperature(double temperature, int chainID){
    if (chainID >= noseHooverChains.size()) {
        throw OpenMMException("Cannot set relative temperature for chainID " + std::to_string(chainID)
                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
        );
    }
    noseHooverChains.at(chainID).setDefaultRelativeTemperature(temperature);

}

double NoseHooverIntegrator::getCollisionFrequency(int chainID) const {
192
193
194
195
196
    if (chainID >= noseHooverChains.size()) {
        throw OpenMMException("Cannot get collision frequency for chainID " + std::to_string(chainID)
                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
        );
    }
197
    return noseHooverChains.at(chainID).getDefaultCollisionFrequency();
198
199
200
}

void NoseHooverIntegrator::setCollisionFrequency(double frequency, int chainID){
201
202
203
204
205
    if (chainID >= noseHooverChains.size()) {
        throw OpenMMException("Cannot set collision frequency for chainID " + std::to_string(chainID)
                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
        );
    }
206
    noseHooverChains.at(chainID).setDefaultCollisionFrequency(frequency);
207
208
}

209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
double NoseHooverIntegrator::getRelativeCollisionFrequency(int chainID) const {
    if (chainID >= noseHooverChains.size()) {
        throw OpenMMException("Cannot get relative collision frequency for chainID " + std::to_string(chainID)
                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
        );
    }
    return noseHooverChains.at(chainID).getDefaultRelativeCollisionFrequency();
}

void NoseHooverIntegrator::setRelativeCollisionFrequency(double frequency, int chainID){
    if (chainID >= noseHooverChains.size()) {
        throw OpenMMException("Cannot set relative collision frequency for chainID " + std::to_string(chainID)
                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
        );
    }
    noseHooverChains.at(chainID).setDefaultRelativeCollisionFrequency(frequency);
}

double NoseHooverIntegrator::computeKineticEnergy() {
228
229
230
    double kE = 0.0;
    if(noseHooverChains.size()) {
        for (const auto &nhc: noseHooverChains){
231
            kE += nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, true).first;
232
233
234
235
236
        }
    } else {
        kE = vvKernel.getAs<IntegrateVelocityVerletStepKernel>().computeKineticEnergy(*context, *this);
    }
    return kE;
237
238
}

239
double NoseHooverIntegrator::computeHeatBathEnergy() {
240
241
    double energy = 0;
    for(auto &nhc : noseHooverChains) {
242
        energy += nhcKernel.getAs<NoseHooverChainKernel>().computeHeatBathEnergy(*context, nhc);
243
244
245
246
    }
    return energy;
}

247
void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
248
249
250
251
252
253
    if (owner != NULL && &contextRef.getOwner() != owner)
        throw OpenMMException("This Integrator is already bound to a context");
    context = &contextRef;
    owner = &contextRef.getOwner();
    vvKernel = context->getPlatform().createKernel(IntegrateVelocityVerletStepKernel::Name(), contextRef);
    vvKernel.getAs<IntegrateVelocityVerletStepKernel>().initialize(contextRef.getSystem(), *this);
254
255
    nhcKernel = context->getPlatform().createKernel(NoseHooverChainKernel::Name(), contextRef);
    nhcKernel.getAs<NoseHooverChainKernel>().initialize();
256
    forcesAreValid = false;
257
258
}

259
void NoseHooverIntegrator::cleanup() {
260
    vvKernel = Kernel();
261
    nhcKernel = Kernel();
262
263
}

264
vector<string> NoseHooverIntegrator::getKernelNames() {
265
    std::vector<std::string> names;
266
    names.push_back(NoseHooverChainKernel::Name());
267
268
269
270
    names.push_back(IntegrateVelocityVerletStepKernel::Name());
    return names;
}

271
void NoseHooverIntegrator::step(int steps) {
272
273
    if (context == NULL)
        throw OpenMMException("This Integrator is not bound to a context!");
274
    std::pair<double, double> scale, kineticEnergy;
275
276
    for (int i = 0; i < steps; ++i) {
        context->updateContextState();
277
        for(auto &nhc : noseHooverChains) {
278
            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, false);
279
280
            scale = nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, nhc, kineticEnergy, getStepSize());
            nhcKernel.getAs<NoseHooverChainKernel>().scaleVelocities(*context, nhc, scale);
281
        }
282
        vvKernel.getAs<IntegrateVelocityVerletStepKernel>().execute(*context, *this, forcesAreValid);
283
        for(auto &nhc : noseHooverChains) {
284
            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, false);
285
286
            scale = nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, nhc, kineticEnergy, getStepSize());
            nhcKernel.getAs<NoseHooverChainKernel>().scaleVelocities(*context, nhc, scale);
287
        }
288
289
    }
}