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

Comparing trunk/src/primitives/Inversion.cpp (file contents):
Revision 1290 by cli2, Wed Sep 10 19:51:45 2008 UTC vs.
Revision 1309 by gezelter, Tue Oct 21 18:23:31 2008 UTC

# Line 40 | Line 40
40   */
41  
42   #include "primitives/Inversion.hpp"
43 – #include "fstream"
43  
44   namespace oopse {
45  
# Line 57 | Line 56 | namespace oopse {
56      // is treated as atom *3* in a standard torsion form:
57  
58      Vector3d pos1 = atom2_->getPos();
59 <    Vector3d pos2 = atom1_->getPos();
60 <    Vector3d pos3 = atom4_->getPos();
61 <    Vector3d pos4 = atom3_->getPos();
59 >    Vector3d pos2 = atom3_->getPos();
60 >    Vector3d pos3 = atom1_->getPos();
61 >    Vector3d pos4 = atom4_->getPos();
62  
63 <   /*std::ofstream myfile;
64 <   myfile.open("Inversion", std::ios::app);      
65 <   myfile << atom1_->getType() << " - atom1; "
67 <              << atom2_->getType() << " - atom2; "
68 <              << atom3_->getType() << " - atom3; "
69 <              << atom4_->getType() << " - atom4; "
70 <              << std::endl;
71 < */
72 <    Vector3d r21 = pos1 - pos2;
73 <    Vector3d r32 = pos2 - pos3;
74 <    Vector3d r42 = pos2 - pos4;
63 >    Vector3d r31 = pos1 - pos3;
64 >    Vector3d r23 = pos3 - pos2;
65 >    Vector3d r43 = pos3 - pos4;
66  
67      //  Calculate the cross products and distances
68 <    Vector3d A = cross(r21, r32);
68 >    Vector3d A = cross(r31, r43);
69      RealType rA = A.length();
70 <    Vector3d B = cross(r32, r42);
70 >    Vector3d B = cross(r43, r23);
71      RealType rB = B.length();
72      //Vector3d C = cross(r23, A);
73      //RealType rC = C.length();
# Line 87 | Line 78 | namespace oopse {
78      
79      //  Calculate the sin and cos
80      RealType cos_phi = dot(A, B) ;
81 <    if (cos_phi > 1.0) {cos_phi = 1.0; std::cout << "!!!! cos_phi is bigger than 1.0"
82 <                                                 << std::endl;}
92 <    if (cos_phi < -1.0) {cos_phi = -1.0; std::cout << "!!!! cos_phi is less than -1.0"
93 <                                                   << std::endl;}
94 <    //std::cout << "We actually use this inversion!!!!" << std::endl;
81 >    if (cos_phi > 1.0) cos_phi = 1.0;
82 >    if (cos_phi < -1.0) cos_phi = -1.0;
83  
84      RealType dVdcosPhi;
97 –    //cos_phi = 2.0*cos_phi*cos_phi - 1.0;
85      inversionType_->calcForce(cos_phi, potential_, dVdcosPhi);
86      Vector3d f1 ;
87      Vector3d f2 ;
# Line 103 | Line 90 | namespace oopse {
90      Vector3d dcosdA = (cos_phi * A - B) /rA;
91      Vector3d dcosdB = (cos_phi * B - A) /rB;
92  
93 <    f1 = dVdcosPhi * cross(r32, dcosdA);
94 <    f2 = dVdcosPhi * ( cross(r42, dcosdB) - cross(r21, dcosdA));
95 <    f3 = dVdcosPhi * cross(dcosdB, r32);
93 >    f1 = dVdcosPhi * cross(r43, dcosdA);
94 >    f2 = dVdcosPhi * ( cross(r23, dcosdB) - cross(r31, dcosdA));
95 >    f3 = dVdcosPhi * cross(dcosdB, r43);
96  
97      // In OOPSE's version of an improper torsion, the central atom
98      // comes first.  However, to get the planarity in a typical cosine
# Line 122 | Line 109 | namespace oopse {
109      atom4_->addFrc(-f2);
110      atom3_->addFrc(-f3);
111  
112 +    atom1_->addParticlePot(potential_);
113 +    atom2_->addParticlePot(potential_);
114 +    atom3_->addParticlePot(potential_);
115 +    atom4_->addParticlePot(potential_);
116 +
117      angle = acos(cos_phi) /M_PI * 180.0;
118    }
119  

Diff Legend

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