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

Comparing trunk/src/primitives/Torsion.cpp (file contents):
Revision 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43 + #include "config.h"
44 + #include <cmath>
45 +
46   #include "primitives/Torsion.hpp"
47  
48   namespace OpenMD {
# Line 47 | Line 51 | namespace OpenMD {
51                     TorsionType *tt) :
52      atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }
53  
54 <  void Torsion::calcForce(RealType& angle) {
54 >  void Torsion::calcForce(RealType& angle, bool doParticlePot) {
55  
56      Vector3d pos1 = atom1_->getPos();
57      Vector3d pos2 = atom2_->getPos();
# Line 63 | Line 67 | namespace OpenMD {
67      RealType rA = A.length();
68      Vector3d B = cross(r32, r43);
69      RealType rB = B.length();
66 –    Vector3d C = cross(r32, A);
67 –    RealType rC = C.length();
70  
71 +    /*
72 +       If either of the two cross product vectors is tiny, that means
73 +       the three atoms involved are colinear, and the torsion angle is
74 +       going to be undefined.  The easiest check for this problem is
75 +       to use the product of the two lengths.
76 +    */
77 +    if (rA * rB < OpenMD::epsilon) return;
78 +    
79      A.normalize();
80 <    B.normalize();
71 <    C.normalize();
80 >    B.normalize();  
81      
82      //  Calculate the sin and cos
83      RealType cos_phi = dot(A, B) ;
84      if (cos_phi > 1.0) cos_phi = 1.0;
85      if (cos_phi < -1.0) cos_phi = -1.0;
86 <
86 >    
87      RealType dVdcosPhi;
88      torsionType_->calcForce(cos_phi, potential_, dVdcosPhi);
89      Vector3d f1 ;
90      Vector3d f2 ;
91      Vector3d f3 ;
92 <
92 >    
93      Vector3d dcosdA = (cos_phi * A - B) /rA;
94      Vector3d dcosdB = (cos_phi * B - A) /rB;
95 <
95 >    
96      f1 = dVdcosPhi * cross(r32, dcosdA);
97      f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
98      f3 = dVdcosPhi * cross(dcosdB, r32);
# Line 92 | Line 101 | namespace OpenMD {
101      atom2_->addFrc(f2 - f1);
102      atom3_->addFrc(f3 - f2);
103      atom4_->addFrc(-f3);
104 <
105 <    atom1_->addParticlePot(potential_);
106 <    atom2_->addParticlePot(potential_);
107 <    atom3_->addParticlePot(potential_);
108 <    atom4_->addParticlePot(potential_);
109 <
110 <    angle = acos(cos_phi) /M_PI * 180.0;
111 <  }
112 <
104 >    
105 >    if (doParticlePot) {
106 >      atom1_->addParticlePot(potential_);
107 >      atom2_->addParticlePot(potential_);
108 >      atom3_->addParticlePot(potential_);
109 >      atom4_->addParticlePot(potential_);
110 >    }
111 >    
112 >    angle = acos(cos_phi) /M_PI * 180.0;    
113 >  }  
114   }

Diff Legend

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