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

Comparing trunk/src/primitives/Inversion.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2067 by gezelter, Thu Mar 5 15:35:37 2015 UTC

# Line 48 | Line 48 | namespace OpenMD {
48   namespace OpenMD {
49  
50    Inversion::Inversion(Atom *atom1, Atom *atom2, Atom *atom3,
51 <                       Atom *atom4, InversionType *it) :
52 <    atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4),
53 <    inversionType_(it) { }
51 >                       Atom *atom4, InversionType *it) :
52 >    ShortRangeInteraction(), inversionType_(it) {
53 >
54 >    atoms_.resize(4);
55 >    atoms_[0] = atom1;
56 >    atoms_[1] = atom2;
57 >    atoms_[2] = atom3;
58 >    atoms_[3] = atom4;
59 >    
60 >    inversionKey_ = inversionType_->getKey();
61 >  }
62    
63    void Inversion::calcForce(RealType& angle, bool doParticlePot) {
64      
# Line 59 | Line 67 | namespace OpenMD {
67      // version of this potential (i.e. Amber-style), the central atom
68      // is treated as atom *3* in a standard torsion form:
69  
70 <    Vector3d pos1 = atom2_->getPos();
71 <    Vector3d pos2 = atom3_->getPos();
72 <    Vector3d pos3 = atom1_->getPos();
73 <    Vector3d pos4 = atom4_->getPos();
70 >    Vector3d pos1 = atoms_[1]->getPos();
71 >    Vector3d pos2 = atoms_[2]->getPos();
72 >    Vector3d pos3 = atoms_[0]->getPos();
73 >    Vector3d pos4 = atoms_[3]->getPos();
74  
75      Vector3d r31 = pos1 - pos3;
76      Vector3d r23 = pos3 - pos2;
# Line 73 | Line 81 | namespace OpenMD {
81      RealType rA = A.length();
82      Vector3d B = cross(r43, r23);
83      RealType rB = B.length();
76 –    //Vector3d C = cross(r23, A);
77 –    //RealType rC = C.length();
84  
85      A.normalize();
86      B.normalize();
81 –    //C.normalize();
87      
88      //  Calculate the sin and cos
89      RealType cos_phi = dot(A, B) ;
# Line 86 | Line 91 | namespace OpenMD {
91      if (cos_phi < -1.0) cos_phi = -1.0;
92  
93      RealType dVdcosPhi;
94 <    inversionType_->calcForce(cos_phi, potential_, dVdcosPhi);
94 >    switch (inversionKey_) {
95 >    case itCosAngle:
96 >      inversionType_->calcForce(cos_phi, potential_, dVdcosPhi);
97 >      break;
98 >    case itAngle:
99 >      RealType phi = acos(cos_phi);
100 >      RealType dVdPhi;
101 >      inversionType_->calcForce(phi, potential_, dVdPhi);
102 >      RealType sin_phi = sqrt(1.0 - cos_phi * cos_phi);  
103 >      if (fabs(sin_phi) < 1.0E-6) {
104 >        sin_phi = 1.0E-6;
105 >      }
106 >      dVdcosPhi = dVdPhi / sin_phi;
107 >      break;
108 >    }
109 >    
110      Vector3d f1 ;
111      Vector3d f2 ;
112      Vector3d f3 ;
# Line 108 | Line 128 | namespace OpenMD {
128  
129      // Confusing enough?  Good.
130  
131 <    atom2_->addFrc(f1);
132 <    atom1_->addFrc(f2 - f1 + f3);
133 <    atom4_->addFrc(-f2);
134 <    atom3_->addFrc(-f3);
131 >    atoms_[1]->addFrc(f1);
132 >    atoms_[0]->addFrc(f2 - f1 + f3);
133 >    atoms_[3]->addFrc(-f2);
134 >    atoms_[2]->addFrc(-f3);
135  
136      if (doParticlePot) {
137 <      atom1_->addParticlePot(potential_);
138 <      atom2_->addParticlePot(potential_);
139 <      atom3_->addParticlePot(potential_);
140 <      atom4_->addParticlePot(potential_);
137 >      atoms_[0]->addParticlePot(potential_);
138 >      atoms_[1]->addParticlePot(potential_);
139 >      atoms_[2]->addParticlePot(potential_);
140 >      atoms_[3]->addParticlePot(potential_);
141      }
142      
143      angle = acos(cos_phi) /M_PI * 180.0;

Diff Legend

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