OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
MolecularRestraint.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 "restraints/MolecularRestraint.hpp"
46
47#include <utility>
48
49#include "math/SVD.hpp"
51
52// using namespace JAMA;
53
54namespace OpenMD {
55
56 void MolecularRestraint::calcForce(std::vector<Vector3d> struc,
57 Vector3d molCom) {
58 assert(struc.size() == ref_.size());
59
60 std::vector<Vector3d>::iterator it;
61
62 // clear out initial values:
63 pot_ = 0.0;
64 for (it = forces_.begin(); it != forces_.end(); ++it)
65 (*it) = 0.0;
66
67 if (restType_ & rtDisplacement) {
68 Vector3d del = molCom - refCom_;
69
70 RealType r = del.length();
71 RealType p = 0.5 * kDisp_ * r * r;
72
73 pot_ += p;
74
75 if (printRest_) restInfo_[rtDisplacement] = std::make_pair(r, p);
76
77 for (it = forces_.begin(); it != forces_.end(); ++it)
78 (*it) += -kDisp_ * del * scaleFactor_;
79 }
80
81 if (restType_ & rtAbsoluteZ) {
82 RealType r = molCom(2) - posZ0_;
83 RealType p = 0.5 * kAbs_ * r * r;
84 Vector3d frc = Vector3d(0.0, 0.0, -kAbs_ * r);
85
86 pot_ += p;
87
88 if (printRest_) restInfo_[rtAbsoluteZ] = std::make_pair(r, p);
89
90 for (it = forces_.begin(); it != forces_.end(); ++it)
91 (*it) += frc * scaleFactor_;
92 }
93
94 // rtDisplacement is 1, rtAbsolute is 2, so anything higher than 3
95 // that requires orientations:
96 if (restType_ > 3) {
97 Vector3d tBody(0.0);
98
99 Mat3x3d R(0.0);
100
101 for (unsigned int n = 0; n < struc.size(); n++) {
102 /*
103 * First migrate the center of mass:
104 */
105 struc[n] -= molCom;
106
107 /*
108 * correlation matrix R:
109 * R(i,j) = sum(over n): y(n,i) * x(n,j)
110 * where x(n) and y(n) are two vector sets
111 */
112
113 R += outProduct(struc[n], ref_[n]);
114 }
115
116 // SVD class uses dynamic matrices, so we must wrap the correlation
117 // matrix before calling SVD and then unwrap the results into Mat3x3d
118 // and Vector3d before we use them.
119
120 DynamicRectMatrix<RealType> Rtmp(3, 3, 0.0);
121 DynamicRectMatrix<RealType> vtmp(3, 3);
122 DynamicVector<RealType> stmp(3);
123 DynamicRectMatrix<RealType> wtmp(3, 3);
124
125 Rtmp.setSubMatrix(0, 0, R);
126
127 // Heavy lifting goes here:
128
129 JAMA::SVD<RealType> svd(Rtmp);
130
131 svd.getU(vtmp);
132 svd.getSingularValues(stmp);
133 svd.getV(wtmp);
134
135 Mat3x3d v;
136 Vector3d s;
137 Mat3x3d w_tr;
138
139 vtmp.getSubMatrix(0, 0, v);
140 stmp.getSubVector(0, s);
141 wtmp.getSubMatrix(0, 0, w_tr);
142
143 bool is_reflection = (v.determinant() * w_tr.determinant()) < 0.0;
144
145 if (is_reflection) {
146 v(2, 0) = -v(2, 0);
147 v(2, 1) = -v(2, 1);
148 v(2, 2) = -v(2, 2);
149 }
150
151 RotMat3x3d Atrans = v * w_tr.transpose();
152 RotMat3x3d A = Atrans.transpose();
153
154 Quat4d quat = A.toQuaternion();
155
156 RealType twistAngle;
157 RealType swingX, swingY;
158
159 quat.toTwistSwing(twistAngle, swingX, swingY);
160
161 RealType dVdtwist, dVdswingX, dVdswingY;
162 RealType dTwist, dSwingX, dSwingY;
163 RealType p;
164
165 if (restType_ & rtTwist) {
166 dTwist = twistAngle - twist0_;
167 /// dVdtwist = kTwist_ * sin(dTwist) ;
168 /// p = kTwist_ * (1.0 - cos(dTwist) ) ;
169 dVdtwist = kTwist_ * dTwist;
170 p = 0.5 * kTwist_ * dTwist * dTwist;
171 pot_ += p;
172 tBody -= dVdtwist * V3Z;
173 if (printRest_) restInfo_[rtTwist] = std::make_pair(twistAngle, p);
174 }
175
176 if (restType_ & rtSwingX) {
177 dSwingX = swingX - swingX0_;
178 /// dVdswingX = kSwingX_ * 2.0 * sin(2.0 * dSwingX);
179 /// p = kSwingX_ * (1.0 - cos(2.0 * dSwingX));
180 dVdswingX = kSwingX_ * dSwingX;
181 p = 0.5 * kSwingX_ * dSwingX * dSwingX;
182 pot_ += p;
183 tBody -= dVdswingX * V3X;
184 if (printRest_) restInfo_[rtSwingX] = std::make_pair(swingX, p);
185 }
186 if (restType_ & rtSwingY) {
187 dSwingY = swingY - swingY0_;
188 /// dVdswingY = kSwingY_ * 2.0 * sin(2.0 * dSwingY);
189 /// p = kSwingY_ * (1.0 - cos(2.0 * dSwingY));
190 dVdswingY = kSwingY_ * dSwingY;
191 p = 0.5 * kSwingX_ * dSwingY * dSwingY;
192 pot_ += p;
193 tBody -= dVdswingY * V3Y;
194 if (printRest_) restInfo_[rtSwingY] = std::make_pair(swingY, p);
195 }
196
197 RealType t2 = dot(tBody, tBody);
198
199 Vector3d rLab, rBody, txr, fBody, fLab;
200
201 for (unsigned int i = 0; i < struc.size(); i++) {
202 rLab = struc[i];
203 rBody = A * rLab;
204
205 txr = cross(tBody, rBody);
206 fBody = txr * t2;
207 fLab = Atrans * fBody;
208 fLab *= scaleFactor_;
209
210 forces_[i] += fLab;
211 }
212
213 // test the force vectors and see if it is the right orientation
214 // std::cout << struc.size() << std::endl << std::endl;
215 // for (int i = 0; i != struc.size(); ++i){
216 // std::cout << "H\t" << struc[i].x() << "\t" << struc[i].y() <<
217 // "\t" << struc[i].z() << "\t"; std::cout << forces_[i].x() <<
218 // "\t"
219 // << forces_[i].y()
220 // << "\t" << forces_[i].z() << std::endl;
221 // }
222 }
223 }
224} // namespace OpenMD
Singular Value Decomposition.
Definition SVD.hpp:33
Quaternion< Real > toQuaternion()
Returns the quaternion from this rotation matrix.
Real length()
Returns the length of this vector.
Definition Vector.hpp:393
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Vector3< Real > cross(const Vector3< Real > &v1, const Vector3< Real > &v2)
Returns the cross product of two Vectors.
Definition Vector3.hpp:136
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.