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

Comparing:
branches/development/src/optimization/PotentialEnergyObjectiveFunction.cpp (file contents), Revision 1741 by gezelter, Tue Jun 5 18:02:44 2012 UTC vs.
trunk/src/optimization/PotentialEnergyObjectiveFunction.cpp (file contents), Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

# Line 45 | Line 45 | namespace OpenMD{
45   namespace OpenMD{
46  
47    PotentialEnergyObjectiveFunction::PotentialEnergyObjectiveFunction(SimInfo* info, ForceManager* forceMan)
48 <    : info_(info), forceMan_(forceMan), thermo(info) {      
48 >    : info_(info), forceMan_(forceMan), thermo(info) {  
49 >    shake_ = new Shake(info_);    
50    }
51    
52  
53    
54    RealType PotentialEnergyObjectiveFunction::value(const DynamicVector<RealType>& x) {
55      setCoor(x);
56 +    shake_->constraintR();
57      forceMan_->calcForces();
58 +    shake_->constraintF();
59      return thermo.getPotential();
60    }
61    
62    void PotentialEnergyObjectiveFunction::gradient(DynamicVector<RealType>& grad, const DynamicVector<RealType>& x) {
63      
64 <    setCoor(x);        
65 <    forceMan_->calcForces();
66 <    getGrad(grad);      
64 >    setCoor(x);      
65 >    shake_->constraintR();
66 >    forceMan_->calcForces();
67 >    shake_->constraintF();
68 >    getGrad(grad);
69    }
70    
71    RealType PotentialEnergyObjectiveFunction::valueAndGradient(DynamicVector<RealType>& grad,
72                                                                const DynamicVector<RealType>& x) {
73 <  
74 <    setCoor(x);    
75 <    forceMan_->calcForces();  
73 >    setCoor(x);
74 >    shake_->constraintR();
75 >    forceMan_->calcForces();
76 >    shake_->constraintF();
77      getGrad(grad);
78      return thermo.getPotential();
79    }
# Line 78 | Line 84 | namespace OpenMD{
84      SimInfo::MoleculeIterator i;
85      Molecule::IntegrableObjectIterator  j;
86      Molecule* mol;
87 <    StuntDouble* integrableObject;    
87 >    StuntDouble* sd;    
88      int index = 0;
89      
90      for (mol = info_->beginMolecule(i); mol != NULL;
91           mol = info_->nextMolecule(i)) {
92 <      for (integrableObject = mol->beginIntegrableObject(j);
93 <           integrableObject != NULL;
94 <           integrableObject = mol->nextIntegrableObject(j)) {
92 >
93 >      for (sd = mol->beginIntegrableObject(j);  sd != NULL;
94 >           sd = mol->nextIntegrableObject(j)) {
95          
96          position[0] = x[index++];
97          position[1] = x[index++];
98          position[2] = x[index++];
99          
100 <        integrableObject->setPos(position);
100 >        sd->setPos(position);
101          
102 <        if (integrableObject->isDirectional()) {
102 >        if (sd->isDirectional()) {
103            eulerAngle[0] = x[index++];
104            eulerAngle[1] = x[index++];
105            eulerAngle[2] = x[index++];
106            
107 <          integrableObject->setEuler(eulerAngle);
108 <          }
107 >          sd->setEuler(eulerAngle);
108 >
109 >          if (sd->isRigidBody()) {
110 >            RigidBody* rb = static_cast<RigidBody*>(sd);
111 >            rb->updateAtoms();
112 >          }        
113 >
114 >        }
115        }
116      }    
117    }
# Line 108 | Line 120 | namespace OpenMD{
120      SimInfo::MoleculeIterator i;
121      Molecule::IntegrableObjectIterator  j;
122      Molecule* mol;
123 <    StuntDouble* integrableObject;    
123 >    StuntDouble* sd;    
124      std::vector<RealType> myGrad;
125      
126      int index = 0;
127      
128      for (mol = info_->beginMolecule(i); mol != NULL;
129           mol = info_->nextMolecule(i)) {
130 <      for (integrableObject = mol->beginIntegrableObject(j);
131 <           integrableObject != NULL;
132 <           integrableObject = mol->nextIntegrableObject(j)) {        
133 <        myGrad = integrableObject->getGrad();
134 <        for (size_t k = 0; k < myGrad.size(); ++k) {
130 >
131 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
132 >           sd = mol->nextIntegrableObject(j)) {        
133 >
134 >        myGrad = sd->getGrad();
135 >
136 >        for (size_t k = 0; k < myGrad.size(); ++k) {  
137            grad[index++] = myGrad[k];
138          }
139 +
140        }            
141      }        
142    }
# Line 130 | Line 145 | namespace OpenMD{
145      SimInfo::MoleculeIterator i;
146      Molecule::IntegrableObjectIterator  j;
147      Molecule* mol;
148 <    StuntDouble* integrableObject;    
148 >    StuntDouble* sd;    
149  
150      Vector3d pos;
151      Vector3d eulerAngle;
# Line 141 | Line 156 | namespace OpenMD{
156      
157      for (mol = info_->beginMolecule(i); mol != NULL;
158           mol = info_->nextMolecule(i)) {
144 –      for (integrableObject = mol->beginIntegrableObject(j);
145 –           integrableObject != NULL;
146 –           integrableObject = mol->nextIntegrableObject(j)) {        
159  
160 <        pos = integrableObject->getPos();
160 >      for (sd = mol->beginIntegrableObject(j);  sd != NULL;
161 >           sd = mol->nextIntegrableObject(j)) {        
162  
163 +        pos = sd->getPos();
164 +
165          xinit[index++] = pos[0];
166          xinit[index++] = pos[1];
167          xinit[index++] = pos[2];
168                  
169 <        if (integrableObject->isDirectional()) {
170 <          eulerAngle = integrableObject->getEuler();
169 >        if (sd->isDirectional()) {
170 >          eulerAngle = sd->getEuler();
171            xinit[index++] = eulerAngle[0];
172            xinit[index++] = eulerAngle[1];
173            xinit[index++] = eulerAngle[2];          
174          }
175 +
176        }
177      }
178      return xinit;

Diff Legend

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