OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
SwitchingFunction.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 "nonbonded/SwitchingFunction.hpp"
46
47#include <cmath>
48#include <cstdio>
49#include <iostream>
50
51#include "utils/simError.h"
52
53using namespace std;
54namespace OpenMD {
55
56 SwitchingFunction::SwitchingFunction() :
57 functionType_(cubic), haveSpline_(false), isCubic_(true), np_(150) {
58 switchSpline_ = std::make_shared<CubicSpline>();
59 }
60
61 void SwitchingFunction::setSwitchType(SwitchingFunctionType sft) {
62 if ((sft == fifth_order_poly) || (sft == cubic)) {
63 if (haveSpline_) {
64 switchSpline_ = std::make_shared<CubicSpline>();
65 setSwitch(rin_, rout_);
66 }
67 } else {
68 snprintf(
69 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
70 "SwitchingFunction::setSwitchType was given unknown function type\n");
71 painCave.severity = OPENMD_ERROR;
72 painCave.isFatal = 1;
73 simError();
74 }
75 }
76
77 void SwitchingFunction::setSwitch(RealType rinner, RealType router) {
78 if (router < rinner) {
79 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
80 "SwitchingFunction::setSwitch was given rinner (%lf) which was\n"
81 "\tlarger than router (%lf).\n",
82 rinner, router);
83 painCave.severity = OPENMD_ERROR;
84 painCave.isFatal = 1;
85 simError();
86 }
87 if (router < 0.0) {
88 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
89 "SwitchingFunction::setSwitch was given router (%lf) which was\n"
90 "\tless than zero.\n",
91 router);
92 painCave.severity = OPENMD_ERROR;
93 painCave.isFatal = 1;
94 simError();
95 }
96 if (rinner < 0.0) {
97 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
98 "SwitchingFunction::setSwitch was given rinner (%lf) which was\n"
99 "\tless than zero.\n",
100 router);
101 painCave.severity = OPENMD_ERROR;
102 painCave.isFatal = 1;
103 simError();
104 }
105
106 rin_ = rinner;
107 rout_ = router;
108 rin2_ = rin_ * rin_;
109 rout2_ = rout_ * rout_;
110
111 if ((router - rinner) < 1.0e-8) {
112 // no reason to set up spline if the switching region is tiny
113 return;
114 }
115
116 // setup r -> sw lookup spline
117 if (functionType_ == fifth_order_poly) {
118 isCubic_ = false;
119 RealType c0 = 1.0;
120 RealType c3 = -10.0;
121 RealType c4 = 15.0;
122 RealType c5 = -6.0;
123
124 RealType dx, r, yval, rval, rval2, rval3, rval4, rval5;
125 RealType rvaldi, rvaldi2, rvaldi3, rvaldi4, rvaldi5;
126
127 dx = (rout_ - rin_) / RealType(np_ - 1);
128
129 for (int i = 0; i < np_; i++) {
130 r = rin_ + RealType(i) * dx;
131 rval = (r - rin_);
132 rval2 = rval * rval;
133 rval3 = rval2 * rval;
134 rval4 = rval2 * rval2;
135 rval5 = rval3 * rval2;
136 rvaldi = 1.0 / (rout_ - rin_);
137 rvaldi2 = rvaldi * rvaldi;
138 rvaldi3 = rvaldi2 * rvaldi;
139 rvaldi4 = rvaldi2 * rvaldi2;
140 rvaldi5 = rvaldi3 * rvaldi2;
141 yval = c0 + c3 * rval3 * rvaldi3 + c4 * rval4 * rvaldi4 +
142 c5 * rval5 * rvaldi5;
143 switchSpline_->addPoint(r, yval);
144 }
145 } else {
146 // cubic splines only need 2 points to do a cubic switching function...
147 isCubic_ = true;
148 switchSpline_->addPoint(rin_, 1.0);
149 switchSpline_->addPoint(rout_, 0.0);
150 }
151 haveSpline_ = true;
152 return;
153 }
154
155 bool SwitchingFunction::getSwitch(const RealType& r2, RealType& sw,
156 RealType& dswdr, RealType& r) {
157 sw = 1.0;
158 dswdr = 0.0;
159
160 bool in_switching_region = false;
161
162 if (r2 > rin2_) {
163 if (r2 > rout2_) {
164 sw = 0.0;
165 } else {
166 in_switching_region = true;
167 r = sqrt(r2);
168 switchSpline_->getValueAndDerivativeAt(r, sw, dswdr);
169 }
170 }
171 return in_switching_region;
172 }
173} // namespace OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.