OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
TorsionTypeParser.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "types/TorsionTypeParser.hpp"
46
47#include <string>
48
49#include "types/CharmmTorsionType.hpp"
50#include "types/CubicTorsionType.hpp"
51#include "types/HarmonicTorsionType.hpp"
56#include "utils/Constants.hpp"
57#include "utils/OpenMDException.hpp"
59#include "utils/StringUtils.hpp"
60
61namespace OpenMD {
62
63 TorsionTypeParser::TorsionTypeParser() : trans180_(true) {
64 stringToEnumMap_["GhostTorsion"] = ttGhostTorsion;
65 stringToEnumMap_["Cubic"] = ttCubic;
66 stringToEnumMap_["Quartic"] = ttQuartic;
67 stringToEnumMap_["Polynomial"] = ttPolynomial;
68 stringToEnumMap_["Charmm"] = ttCharmm;
69 stringToEnumMap_["Opls"] = ttOpls;
70 stringToEnumMap_["Trappe"] = ttTrappe;
71 stringToEnumMap_["Harmonic"] = ttHarmonic;
72 }
73
74 TorsionType* TorsionTypeParser::parseTypeAndPars(const std::string& type,
75 std::vector<RealType> pars) {
76 std::string line(type);
77
78 std::vector<RealType>::iterator it;
79 for (it = pars.begin(); it != pars.end(); ++it) {
80 line.append("\t");
81 line.append(std::to_string(*it));
82 }
83
84 return parseLine(line);
85 }
86
87 TorsionType* TorsionTypeParser::parseLine(const std::string& line) {
88 StringTokenizer tokenizer(line);
89 TorsionType* torsionType = NULL;
90 int nTokens = tokenizer.countTokens();
91
92 if (nTokens < 1) {
93 throw OpenMDException("TorsionTypeParser: Not enough tokens");
94 }
95
96 TorsionTypeEnum tt = getTorsionTypeEnum(tokenizer.nextToken());
97
98 nTokens -= 1;
99
100 switch (tt) {
101 case ttGhostTorsion:
102 if (nTokens < 4) {
103 throw OpenMDException("TorsionTypeParser: Not enough tokens");
104 } else {
105 RealType k3 = tokenizer.nextTokenAsDouble();
106 RealType k2 = tokenizer.nextTokenAsDouble();
107 RealType k1 = tokenizer.nextTokenAsDouble();
108 RealType k0 = tokenizer.nextTokenAsDouble();
109
110 if (trans180_)
111 torsionType = new CubicTorsionType(k3, k2, k1, k0);
112 else
113 torsionType = new CubicTorsionType(-k3, k2, -k1, k0);
114 }
115 break;
116
117 case ttCubic:
118 if (nTokens < 4) {
119 throw OpenMDException("TorsionTypeParser: Not enough tokens");
120 } else {
121 RealType k3 = tokenizer.nextTokenAsDouble();
122 RealType k2 = tokenizer.nextTokenAsDouble();
123 RealType k1 = tokenizer.nextTokenAsDouble();
124 RealType k0 = tokenizer.nextTokenAsDouble();
125
126 if (trans180_)
127 torsionType = new CubicTorsionType(k3, k2, k1, k0);
128 else
129 torsionType = new CubicTorsionType(-k3, k2, -k1, k0);
130 }
131 break;
132
133 case ttQuartic:
134 if (nTokens < 5) {
135 throw OpenMDException("TorsionTypeParser: Not enough tokens");
136 } else {
137 RealType k4 = tokenizer.nextTokenAsDouble();
138 RealType k3 = tokenizer.nextTokenAsDouble();
139 RealType k2 = tokenizer.nextTokenAsDouble();
140 RealType k1 = tokenizer.nextTokenAsDouble();
141 RealType k0 = tokenizer.nextTokenAsDouble();
142
143 if (trans180_)
144 torsionType = new QuarticTorsionType(k4, k3, k2, k1, k0);
145 else
146 torsionType = new QuarticTorsionType(k4, -k3, k2, -k1, k0);
147 }
148 break;
149
150 case ttPolynomial:
151 if (nTokens < 2 || nTokens % 2 != 0) {
152 throw OpenMDException("TorsionTypeParser: Not enough tokens");
153 } else {
154 int nPairs = nTokens / 2;
155 int power;
156 RealType coefficient;
157 torsionType = new PolynomialTorsionType();
158
160 dynamic_cast<PolynomialTorsionType*>(torsionType);
161
162 for (int i = 0; i < nPairs; ++i) {
163 power = tokenizer.nextTokenAsInt();
164 bool isOdd = power % 2 == 0 ? false : true;
165
166 coefficient = tokenizer.nextTokenAsDouble();
167
168 if (!trans180_ && isOdd) coefficient = -coefficient;
169
170 ptt->setCoefficient(power, coefficient);
171 }
172 }
173
174 break;
175
176 case ttCharmm:
177
178 if (nTokens < 3 || nTokens % 3 != 0) {
179 throw OpenMDException("TorsionTypeParser: Not enough tokens");
180 } else {
181 int nSets = nTokens / 3;
182
183 std::vector<CharmmTorsionParameter> parameters;
184 for (int i = 0; i < nSets; ++i) {
185 CharmmTorsionParameter currParam;
186 currParam.kchi = tokenizer.nextTokenAsDouble();
187 currParam.n = tokenizer.nextTokenAsInt();
188 currParam.delta = tokenizer.nextTokenAsDouble() / 180.0 *
189 Constants::PI; // convert to rad
190
191 bool isOdd = currParam.n % 2 == 0 ? false : true;
192 if (!trans180_) {
193 currParam.delta = Constants::PI - currParam.delta;
194 if (isOdd) currParam.kchi = -currParam.kchi;
195 }
196
197 parameters.push_back(currParam);
198 }
199
200 torsionType = new CharmmTorsionType(parameters);
201 }
202
203 break;
204
205 case ttOpls:
206
207 if (nTokens < 3) {
208 throw OpenMDException("TorsionTypeParser: Not enough tokens");
209 } else {
210 RealType v1 = tokenizer.nextTokenAsDouble();
211 RealType v2 = tokenizer.nextTokenAsDouble();
212 RealType v3 = tokenizer.nextTokenAsDouble();
213
214 torsionType = new OplsTorsionType(v1, v2, v3, trans180_);
215 }
216
217 break;
218
219 case ttTrappe:
220
221 if (nTokens < 4) {
222 throw OpenMDException("TorsionTypeParser: Not enough tokens");
223 } else {
224 RealType c0 = tokenizer.nextTokenAsDouble();
225 RealType c1 = tokenizer.nextTokenAsDouble();
226 RealType c2 = tokenizer.nextTokenAsDouble();
227 RealType c3 = tokenizer.nextTokenAsDouble();
228 torsionType = new TrappeTorsionType(c0, c1, c2, c3, trans180_);
229 }
230
231 break;
232
233 case ttHarmonic:
234
235 if (nTokens < 2) {
236 throw OpenMDException("TorsionTypeParser: Not enough tokens");
237 } else {
238 // Most torsions don't have specific angle information since
239 // they are cosine polynomials. This one is different,
240 // however. To match our other force field files
241 // (particularly for bends):
242 //
243 // phi0 should be read in degrees
244 // d0 should be read in kcal / mol / degrees^2
245
246 RealType degreesPerRadian = 180.0 / Constants::PI;
247
248 // convert to radians
249 RealType phi0 = tokenizer.nextTokenAsDouble() / degreesPerRadian;
250
251 // convert to kcal / mol / radians^2
252 RealType d0 = tokenizer.nextTokenAsDouble() * pow(degreesPerRadian, 2);
253
254 if (!trans180_) phi0 = Constants::PI - phi0;
255
256 torsionType = new HarmonicTorsionType(d0, phi0);
257 }
258
259 break;
260
261 case ttUnknown:
262 default:
263 throw OpenMDException("TorsionTypeParser: Unknown Torsion Type");
264 }
265 return torsionType;
266 }
267
268 TorsionTypeParser::TorsionTypeEnum TorsionTypeParser::getTorsionTypeEnum(
269 const std::string& str) {
270 std::map<std::string, TorsionTypeEnum>::iterator i;
271 i = stringToEnumMap_.find(str);
272
273 return i == stringToEnumMap_.end() ? ttUnknown : i->second;
274 }
275
276} // namespace OpenMD
"types/CharmmTorsionType.hpp" These torsion types are defined identically with functional form given ...
These torsion types are defined identically with functional form given in equation 5 in the following...
These torsion types are defined identically with functional form given in the following paper:
"types/PolynomialTorsionType.hpp"
The string tokenizer class allows an application to break a string into tokens The set of delimiters ...
"types/TrappeTorsionType.hpp" These torsion types are defined identically with functional form given ...
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.