ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/primitives/Torsion.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/primitives/Torsion.cpp (file contents):
Revision 1492 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC

# Line 1 | Line 1
1 < #include "primitives/SRI.hpp"
2 < #include "primitives/Atom.hpp"
3 < #include <math.h>
4 < #include <iostream>
5 < #include <stdlib.h>
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/Torsion.hpp"
43  
44 < void Torsion::set_atoms( Atom &a, Atom &b, Atom &c, Atom &d){
8 <  c_p_a = &a;
9 <  c_p_b = &b;
10 <  c_p_c = &c;
11 <  c_p_d = &d;
12 < }
44 > namespace oopse {
45  
46 + Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4,
47 +                 TorsionType *tt) :
48 +    atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }
49  
50 < void Torsion::calc_forces(){
51 <  
52 <  /**********************************************************************
53 <   *
54 <   * initialize vectors
20 <   *
21 <   ***********************************************************************/
22 <  
23 <  vect r_ab; /* the vector whose origin is a and end is b */
24 <  vect r_cb; /* the vector whose origin is c and end is b */
25 <  vect r_cd; /* the vector whose origin is c and end is b */
26 <  vect r_cr1; /* the cross product of r_ab and r_cb */
27 <  vect r_cr2; /* the cross product of r_cb and r_cd */
50 > void Torsion::calcForce() {
51 >    Vector3d pos1 = atom1_->getPos();
52 >    Vector3d pos2 = atom2_->getPos();
53 >    Vector3d pos3 = atom3_->getPos();
54 >    Vector3d pos4 = atom4_->getPos();
55  
56 <  double r_cr1_x2; /* the components of r_cr1 squared */
57 <  double r_cr1_y2;
58 <  double r_cr1_z2;
32 <  
33 <  double r_cr2_x2; /* the components of r_cr2 squared */
34 <  double r_cr2_y2;
35 <  double r_cr2_z2;
56 >    Vector3d r21 = pos1 - pos2;
57 >    Vector3d r32 = pos2 - pos3;
58 >    Vector3d r43 = pos3 - pos4;
59  
60 <  double r_cr1_sqr; /* the length of r_cr1 squared */
61 <  double r_cr2_sqr; /* the length of r_cr2 squared */
62 <  
63 <  double r_cr1_r_cr2; /* the length of r_cr1 * length of r_cr2 */
64 <  
65 <  double aR[3], bR[3], cR[3], dR[3];
66 <  double aF[3], bF[3], cF[3], dF[3];
60 >    //  Calculate the cross products and distances
61 >    Vector3d A = cross(r21, r32);
62 >    double rA = A.length();
63 >    Vector3d B = cross(r32, r43);
64 >    double rB = B.length();
65 >    Vector3d C = cross(r32, A);
66 >    double rC = C.length();
67  
68 <  c_p_a->getPos( aR );
69 <  c_p_b->getPos( bR );
70 <  c_p_c->getPos( cR );
71 <  c_p_d->getPos( dR );
68 >    A.normalize();
69 >    B.normalize();
70 >    C.normalize();
71 >    
72 >    //  Calculate the sin and cos
73 >    double cos_phi = dot(A, B) ;
74 >    double sin_phi = dot(C, B);
75  
76 <  r_ab.x = bR[0] - aR[0];
77 <  r_ab.y = bR[1] - aR[1];
52 <  r_ab.z = bR[2] - aR[2];
53 <  r_ab.length  = sqrt((r_ab.x * r_ab.x + r_ab.y * r_ab.y + r_ab.z * r_ab.z));
76 >    double dVdPhi;
77 >    torsionType_->calcForce(cos_phi, sin_phi, potential_, dVdPhi);
78  
79 <  r_cb.x = bR[0] - cR[0];
80 <  r_cb.y = bR[1] - cR[1];
81 <  r_cb.z = bR[2] - cR[2];
58 <  r_cb.length = sqrt((r_cb.x * r_cb.x + r_cb.y * r_cb.y + r_cb.z * r_cb.z));
79 >    Vector3d f1;
80 >    Vector3d f2;
81 >    Vector3d f3;
82  
83 <  r_cd.x = dR[0] - cR[0];
84 <  r_cd.y = dR[1] - cR[1];
85 <  r_cd.z = dR[2] - cR[2];
86 <  r_cd.length = sqrt((r_cd.x * r_cd.x + r_cd.y * r_cd.y + r_cd.z * r_cd.z));
83 >    //  Next, we want to calculate the forces.  In order
84 >    //  to do that, we first need to figure out whether the
85 >    //  sin or cos form will be more stable.  For this,
86 >    //  just look at the value of phi
87 >    //if (fabs(sin_phi) > 0.1) {
88 >        //  use the sin version to avoid 1/cos terms
89  
90 <  r_cr1.x = r_ab.y * r_cb.z - r_cb.y * r_ab.z;
91 <  r_cr1.y = r_ab.z * r_cb.x - r_cb.z * r_ab.x;
67 <  r_cr1.z = r_ab.x * r_cb.y - r_cb.x * r_ab.y;
68 <  r_cr1_x2 = r_cr1.x * r_cr1.x;
69 <  r_cr1_y2 = r_cr1.y * r_cr1.y;
70 <  r_cr1_z2 = r_cr1.z * r_cr1.z;
71 <  r_cr1_sqr = r_cr1_x2 + r_cr1_y2 + r_cr1_z2;
72 <  r_cr1.length = sqrt(r_cr1_sqr);
90 >        Vector3d dcosdA = (cos_phi * A - B) /rA;
91 >        Vector3d dcosdB = (cos_phi * B - A) /rB;
92  
93 <  r_cr2.x = r_cb.y * r_cd.z - r_cd.y * r_cb.z;
75 <  r_cr2.y = r_cb.z * r_cd.x - r_cd.z * r_cb.x;
76 <  r_cr2.z = r_cb.x * r_cd.y - r_cd.x * r_cb.y;
77 <  r_cr2_x2 = r_cr2.x * r_cr2.x;
78 <  r_cr2_y2 = r_cr2.y * r_cr2.y;
79 <  r_cr2_z2 = r_cr2.z * r_cr2.z;
80 <  r_cr2_sqr = r_cr2_x2 + r_cr2_y2 + r_cr2_z2;
81 <  r_cr2.length = sqrt(r_cr2_sqr);
93 >        double dVdcosPhi = -dVdPhi / sin_phi;
94  
95 <  r_cr1_r_cr2 = r_cr1.length * r_cr2.length;
95 >        f1 = dVdcosPhi * cross(r32, dcosdA);
96 >        f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
97 >        f3 = dVdcosPhi * cross(dcosdB, r32);
98  
99 <  /**********************************************************************
100 <   *
101 <   * dot product and angle calculations
102 <   *
89 <   ***********************************************************************/
90 <  
91 <  double cr1_dot_cr2; /* the dot product of the cr1 and cr2 vectors */
92 <  double cos_phi; /* the cosine of the torsion angle */
99 >    /** @todo fix below block, must be something wrong with the sign somewhere */
100 >    //} else {
101 >        //  This angle is closer to 0 or 180 than it is to
102 >        //  90, so use the cos version to avoid 1/sin terms
103  
104 <  cr1_dot_cr2 = r_cr1.x * r_cr2.x + r_cr1.y * r_cr2.y + r_cr1.z * r_cr2.z;
105 <  
106 <  cos_phi = cr1_dot_cr2 / r_cr1_r_cr2;
97 <  
98 <   /* adjust for the granularity of the numbers for angles near 0 or pi */
104 >        //double dVdsinPhi = dVdPhi /cos_phi;
105 >        //Vector3d dsindB = (sin_phi * B - C) /rB;
106 >        //Vector3d dsindC = (sin_phi * C - B) /rC;
107  
108 <  if(cos_phi > 1.0) cos_phi = 1.0;
101 <  if(cos_phi < -1.0) cos_phi = -1.0;
108 >        //f1.x() = dVdsinPhi*((r32.y()*r32.y() + r32.z()*r32.z())*dsindC.x() - r32.x()*r32.y()*dsindC.y() - r32.x()*r32.z()*dsindC.z());
109  
110 +        //f1.y() = dVdsinPhi*((r32.z()*r32.z() + r32.x()*r32.x())*dsindC.y() - r32.y()*r32.z()*dsindC.z() - r32.y()*r32.x()*dsindC.x());
111  
112 <  /********************************************************************
105 <   *
106 <   * This next section calculates derivatives needed for the force
107 <   * calculation
108 <   *
109 <   ********************************************************************/
112 >        //f1.z() = dVdsinPhi*((r32.x()*r32.x() + r32.y()*r32.y())*dsindC.z() - r32.z()*r32.x()*dsindC.x() - r32.z()*r32.y()*dsindC.y());
113  
114 +        //f2.x() = dVdsinPhi*(-(r32.y()*r21.y() + r32.z()*r21.z())*dsindC.x() + (2.0*r32.x()*r21.y() - r21.x()*r32.y())*dsindC.y()
115 +        //+ (2.0*r32.x()*r21.z() - r21.x()*r32.z())*dsindC.z() + dsindB.z()*r43.y() - dsindB.y()*r43.z());
116  
117 <  /* the derivatives of cos phi with respect to the x, y,
118 <     and z components of vectors cr1 and cr2. */
114 <  double d_cos_dx_cr1;
115 <  double d_cos_dy_cr1;
116 <  double d_cos_dz_cr1;
117 <  double d_cos_dx_cr2;
118 <  double d_cos_dy_cr2;
119 <  double d_cos_dz_cr2;
117 >        //f2.y() = dVdsinPhi*(-(r32.z()*r21.z() + r32.x()*r21.x())*dsindC.y() + (2.0*r32.y()*r21.z() - r21.y()*r32.z())*dsindC.z()
118 >        //+ (2.0*r32.y()*r21.x() - r21.y()*r32.x())*dsindC.x() + dsindB.x()*r43.z() - dsindB.z()*r43.x());
119  
120 <  d_cos_dx_cr1 = r_cr2.x / r_cr1_r_cr2 - (cos_phi * r_cr1.x) / r_cr1_sqr;
121 <  d_cos_dy_cr1 = r_cr2.y / r_cr1_r_cr2 - (cos_phi * r_cr1.y) / r_cr1_sqr;
123 <  d_cos_dz_cr1 = r_cr2.z / r_cr1_r_cr2 - (cos_phi * r_cr1.z) / r_cr1_sqr;
120 >        //f2.z() = dVdsinPhi*(-(r32.x()*r21.x() + r32.y()*r21.y())*dsindC.z() + (2.0*r32.z()*r21.x() - r21.z()*r32.x())*dsindC.x()
121 >        //+(2.0*r32.z()*r21.y() - r21.z()*r32.y())*dsindC.y() + dsindB.y()*r43.x() - dsindB.x()*r43.y());
122  
123 <  d_cos_dx_cr2 = r_cr1.x / r_cr1_r_cr2 - (cos_phi * r_cr2.x) / r_cr2_sqr;
126 <  d_cos_dy_cr2 = r_cr1.y / r_cr1_r_cr2 - (cos_phi * r_cr2.y) / r_cr2_sqr;
127 <  d_cos_dz_cr2 = r_cr1.z / r_cr1_r_cr2 - (cos_phi * r_cr2.z) / r_cr2_sqr;
123 >        //f3 = dVdsinPhi * cross(r32, dsindB);
124  
125 <  /***********************************************************************
130 <   *
131 <   * Calculate the actual forces and place them in the atoms.
132 <   *
133 <   ***********************************************************************/
125 >    //}
126  
127 <  double force; /*the force scaling factor */
127 >    atom1_->addFrc(f1);
128 >    atom2_->addFrc(f2 - f1);
129 >    atom3_->addFrc(f3 - f2);
130 >    atom4_->addFrc(-f3);
131 > }
132  
137 –  force = torsion_force(cos_phi);
138 –
139 –  aF[0] = force * (d_cos_dy_cr1 * r_cb.z - d_cos_dz_cr1 * r_cb.y);
140 –  aF[1] = force * (d_cos_dz_cr1 * r_cb.x - d_cos_dx_cr1 * r_cb.z);
141 –  aF[2] = force * (d_cos_dx_cr1 * r_cb.y - d_cos_dy_cr1 * r_cb.x);
142 –
143 –  bF[0] = force * (  d_cos_dy_cr1 * (r_ab.z - r_cb.z)
144 –                   - d_cos_dy_cr2 *  r_cd.z      
145 –                   + d_cos_dz_cr1 * (r_cb.y - r_ab.y)
146 –                   + d_cos_dz_cr2 *  r_cd.y);
147 –  bF[1] = force * (  d_cos_dx_cr1 * (r_cb.z - r_ab.z)
148 –                   + d_cos_dx_cr2 *  r_cd.z      
149 –                   + d_cos_dz_cr1 * (r_ab.x - r_cb.x)
150 –                   - d_cos_dz_cr2 *  r_cd.x);
151 –  bF[2] = force * (  d_cos_dx_cr1 * (r_ab.y - r_cb.y)
152 –                   - d_cos_dx_cr2 *  r_cd.y      
153 –                   + d_cos_dy_cr1 * (r_cb.x - r_ab.x)
154 –                   + d_cos_dy_cr2 *  r_cd.x);
155 –
156 –  cF[0] = force * (- d_cos_dy_cr1 *  r_ab.z
157 –                   - d_cos_dy_cr2 * (r_cb.z - r_cd.z)
158 –                   + d_cos_dz_cr1 *  r_ab.y
159 –                   - d_cos_dz_cr2 * (r_cd.y - r_cb.y));
160 –  cF[1] = force * (  d_cos_dx_cr1 *  r_ab.z
161 –                   - d_cos_dx_cr2 * (r_cd.z - r_cb.z)
162 –                   - d_cos_dz_cr1 *  r_ab.x
163 –                   - d_cos_dz_cr2 * (r_cb.x - r_cd.x));
164 –  cF[2] = force * (- d_cos_dx_cr1 *  r_ab.y
165 –                   - d_cos_dx_cr2 * (r_cb.y - r_cd.y)
166 –                   + d_cos_dy_cr1 *  r_ab.x
167 –                   - d_cos_dy_cr2 * (r_cd.x - r_cb.x));
168 –
169 –  dF[0] = force * (d_cos_dy_cr2 * r_cb.z - d_cos_dz_cr2 * r_cb.y);
170 –  dF[1] = force * (d_cos_dz_cr2 * r_cb.x - d_cos_dx_cr2 * r_cb.z);
171 –  dF[2] = force * (d_cos_dx_cr2 * r_cb.y - d_cos_dy_cr2 * r_cb.x);
172 –
173 –
174 –  c_p_a->addFrc(aF);
175 –  c_p_b->addFrc(bF);
176 –  c_p_c->addFrc(cF);
177 –  c_p_d->addFrc(dF);
133   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines