# | Line 1 | Line 1 | |
---|---|---|
1 | < | #include "primitives/SRI.hpp" |
2 | < | #include "primitives/Atom.hpp" |
1 | > | /* |
2 | > | * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 | > | * |
4 | > | * The University of Notre Dame grants you ("Licensee") a |
5 | > | * non-exclusive, royalty free, license to use, modify and |
6 | > | * redistribute this software in source and binary code form, provided |
7 | > | * that the following conditions are met: |
8 | > | * |
9 | > | * 1. Acknowledgement of the program authors must be made in any |
10 | > | * publication of scientific results based in part on use of the |
11 | > | * program. An acceptable form of acknowledgement is citation of |
12 | > | * the article in which the program was described (Matthew |
13 | > | * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
14 | > | * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
15 | > | * Parallel Simulation Engine for Molecular Dynamics," |
16 | > | * J. Comput. Chem. 26, pp. 252-271 (2005)) |
17 | > | * |
18 | > | * 2. Redistributions of source code must retain the above copyright |
19 | > | * notice, this list of conditions and the following disclaimer. |
20 | > | * |
21 | > | * 3. Redistributions in binary form must reproduce the above copyright |
22 | > | * notice, this list of conditions and the following disclaimer in the |
23 | > | * documentation and/or other materials provided with the |
24 | > | * distribution. |
25 | > | * |
26 | > | * This software is provided "AS IS," without a warranty of any |
27 | > | * kind. All express or implied conditions, representations and |
28 | > | * warranties, including any implied warranty of merchantability, |
29 | > | * fitness for a particular purpose or non-infringement, are hereby |
30 | > | * excluded. The University of Notre Dame and its licensors shall not |
31 | > | * be liable for any damages suffered by licensee as a result of |
32 | > | * using, modifying or distributing the software or its |
33 | > | * derivatives. In no event will the University of Notre Dame or its |
34 | > | * licensors be liable for any lost revenue, profit or data, or for |
35 | > | * direct, indirect, special, consequential, incidental or punitive |
36 | > | * damages, however caused and regardless of the theory of liability, |
37 | > | * arising out of the use of or inability to use software, even if the |
38 | > | * University of Notre Dame has been advised of the possibility of |
39 | > | * such damages. |
40 | > | */ |
41 | > | |
42 | > | #include "primitives/Bend.hpp" |
43 | ||
44 | < | #include <math.h> |
5 | < | #include <iostream> |
6 | < | #include <stdlib.h> |
44 | > | namespace oopse { |
45 | ||
46 | < | void Bend::set_atoms( Atom &a, Atom &b, Atom &c){ |
46 | > | /**@todo still a lot left to improve*/ |
47 | > | void Bend::calcForce(double& angle) { |
48 | > | Vector3d pos1 = atom1_->getPos(); |
49 | > | Vector3d pos2 = atom2_->getPos(); |
50 | > | Vector3d pos3 = atom3_->getPos(); |
51 | ||
52 | < | c_p_a = &a; |
53 | < | c_p_b = &b; |
12 | < | c_p_c = &c; |
52 | > | Vector3d r21 = pos1 - pos2; |
53 | > | double d21 = r21.length(); |
54 | ||
55 | < | c_potential_E = 0.0; |
15 | < | } |
55 | > | double d21inv = 1.0 / d21; |
56 | ||
57 | + | Vector3d r23 = pos3 - pos2; |
58 | + | double d23 = r23.length(); |
59 | ||
60 | < | void Bend::calc_forces(){ |
19 | < | |
20 | < | double dx,dy,dz,gx,gy,gz,dx2,dy2,dz2,gx2,gy2,gz2; |
21 | < | double rij2, rkj2, riji2, rkji2, dot, denom, cosang, angl; |
22 | < | |
23 | < | double sina2, sinai; |
60 | > | double d23inv = 1.0 / d23; |
61 | ||
62 | < | double comf2, comf3, comf4; |
26 | < | double dcsidx, dcsidy, dcsidz, dcskdx, dcskdy, dcskdz; |
27 | < | // double dcsjdx, dcsjdy, dcsjdz; |
28 | < | double dadxi, dadyi, dadzi; |
29 | < | double dadxk, dadyk, dadzk;//, dadxj, dadyj, dadzj; |
30 | < | double daxi, dayi, dazi, daxk, dayk, dazk, daxj, dayj, dazj; |
31 | < | |
32 | < | double aR[3], bR[3], cR[3]; |
33 | < | double aF[3], bF[3], cF[3]; |
62 | > | double cosTheta = dot(r21, r23) / (d21 * d23); |
63 | ||
64 | < | c_p_a->getPos( aR ); |
65 | < | c_p_b->getPos( bR ); |
66 | < | c_p_c->getPos( cR ); |
67 | < | |
64 | > | //check roundoff |
65 | > | if (cosTheta > 1.0) { |
66 | > | cosTheta = 1.0; |
67 | > | } else if (cosTheta < -1.0) { |
68 | > | cosTheta = -1.0; |
69 | > | } |
70 | ||
71 | < | dx = aR[0] - bR[0]; |
41 | < | dy = aR[1] - bR[1]; |
42 | < | dz = aR[2] - bR[2]; |
43 | < | |
44 | < | gx = cR[0] - bR[0]; |
45 | < | gy = cR[1] - bR[1]; |
46 | < | gz = cR[2] - bR[2]; |
47 | < | |
48 | < | dx2 = dx * dx; |
49 | < | dy2 = dy * dy; |
50 | < | dz2 = dz * dz; |
71 | > | double theta = acos(cosTheta); |
72 | ||
73 | < | gx2 = gx * gx; |
53 | < | gy2 = gy * gy; |
54 | < | gz2 = gz * gz; |
55 | < | |
56 | < | rij2 = dx2 + dy2 + dz2; |
57 | < | rkj2 = gx2 + gy2 + gz2; |
58 | < | |
59 | < | riji2 = 1.0 / rij2; |
60 | < | rkji2 = 1.0 / rkj2; |
73 | > | double dVdTheta; |
74 | ||
75 | < | dot = dx * gx + dy * gy + dz * gz; |
76 | < | denom = sqrt((riji2 * rkji2)); |
77 | < | cosang = dot * denom; |
75 | > | bendType_->calcForce(theta, potential_, dVdTheta); |
76 | > | //std::cout << atom1_->getType() << "\t" << atom2_->getType() << "\t" << atom3_->getType() << "\t"; |
77 | > | //std::cout << "theta = " << theta/M_PI * 180.0 <<", potential = " << potential_ << std::endl; |
78 | ||
79 | < | if(cosang > 1.0)cosang = 1.0; |
67 | < | if(cosang < -1.0) cosang = -1.0; |
79 | > | double sinTheta = sqrt(1.0 - cosTheta * cosTheta); |
80 | ||
81 | < | angl = acos(cosang); |
82 | < | angl = angl * 180.0 / M_PI; |
81 | > | if (fabs(sinTheta) < 1.0E-6) { |
82 | > | sinTheta = 1.0E-6; |
83 | > | } |
84 | ||
85 | < | sina2 = 1.0 - cosang*cosang; |
86 | < | if(fabs(sina2) < 1.0E-12 ) sina2 = 1.0E-12; |
74 | < | sinai = 1.0 / sqrt(sina2); |
85 | > | double commonFactor1 = dVdTheta / sinTheta * d21inv; |
86 | > | double commonFactor2 = dVdTheta / sinTheta * d23inv; |
87 | ||
88 | < | comf2 = cosang * riji2; |
89 | < | comf3 = cosang * rkji2; |
78 | < | comf4 = bend_force(angl); |
88 | > | Vector3d force1 = commonFactor1 * (r23 * d23inv - r21*d21inv*cosTheta); |
89 | > | Vector3d force3 = commonFactor2 * (r21 * d21inv - r23*d23inv*cosTheta); |
90 | ||
91 | + | //total force in current bend is zero |
92 | + | Vector3d force2 = force1 + force3; |
93 | + | force2 *= -1.0; |
94 | ||
95 | < | dcsidx = gx*denom - comf2*dx; |
96 | < | dcsidy = gy*denom - comf2*dy; |
97 | < | dcsidz = gz*denom - comf2*dz; |
84 | < | |
85 | < | dcskdx = dx*denom - comf3*gx; |
86 | < | dcskdy = dy*denom - comf3*gy; |
87 | < | dcskdz = dz*denom - comf3*gz; |
88 | < | |
89 | < | // dcsjdx = -dcsidx - dcskdx; |
90 | < | // dcsjdy = -dcsidy - dcskdy; |
91 | < | // dcsjdz = -dcsidz - dcskdz; |
95 | > | atom1_->addFrc(force1); |
96 | > | atom2_->addFrc(force2); |
97 | > | atom3_->addFrc(force3); |
98 | ||
99 | < | dadxi = -sinai*dcsidx; |
100 | < | dadyi = -sinai*dcsidy; |
95 | < | dadzi = -sinai*dcsidz; |
99 | > | angle = theta /M_PI * 180.0; |
100 | > | } |
101 | ||
102 | < | dadxk = -sinai*dcskdx; |
98 | < | dadyk = -sinai*dcskdy; |
99 | < | dadzk = -sinai*dcskdz; |
100 | < | |
101 | < | // dadxj = -dadxi - dadxk; |
102 | < | // dadyj = -dadyi - dadyk; |
103 | < | // dadzj = -dadzi - dadzk; |
104 | < | |
105 | < | daxi = comf4*dadxi; |
106 | < | dayi = comf4*dadyi; |
107 | < | dazi = comf4*dadzi; |
108 | < | |
109 | < | daxk = comf4*dadxk; |
110 | < | dayk = comf4*dadyk; |
111 | < | dazk = comf4*dadzk; |
112 | < | |
113 | < | daxj = -daxi - daxk; |
114 | < | dayj = -dayi - dayk; |
115 | < | dazj = -dazi - dazk; |
116 | < | |
117 | < | aF[0] = daxi; |
118 | < | aF[1] = dayi; |
119 | < | aF[2] = dazi; |
120 | < | |
121 | < | bF[0] = daxj; |
122 | < | bF[1] = dayj; |
123 | < | bF[2] = dazj; |
124 | < | |
125 | < | cF[0] = daxk; |
126 | < | cF[1] = dayk; |
127 | < | cF[2] = dazk; |
128 | < | |
129 | < | c_p_a->addFrc(aF); |
130 | < | c_p_b->addFrc(bF); |
131 | < | c_p_c->addFrc(cF); |
132 | < | |
133 | < | return; |
134 | < | } |
102 | > | } //end namespace oopse |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |