ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/OOPSEMinimizer.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/OOPSEMinimizer.cpp (file contents):
Revision 1064 by tim, Tue Feb 24 15:44:45 2004 UTC vs.
Revision 1108 by tim, Wed Apr 14 15:37:41 2004 UTC

# Line 5 | Line 5 | OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, Forc
5                                                MinimizerParameterSet * param)
6                       :RealIntegrator(theInfo, the_ff), bVerbose(false), bShake(true){
7  
8 –  atoms = info->atoms;
9 –
8    tStats = new Thermo(info);
9    dumpOut = new DumpWriter(info);
10    statOut = new StatWriter(info);
# Line 18 | Line 16 | OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, Forc
16    curX = getCoor();
17    curG.resize(ndim);
18  
19 +  preMove();
20   }
21  
22   OOPSEMinimizer::~OOPSEMinimizer(){
# Line 53 | Line 52 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
52  
53    index = 0;
54  
55 <  for(int i = 0; i < nAtoms; i++){
55 >  for(int i = 0; i < integrableObjects.size(); i++){
56  
57 <    if(atoms[i]->isDirectional()){
59 <      dAtom = (DirectionalAtom*) atoms[i];
60 <      dAtom->getGrad(dAtomGrad);
57 >    if (integrableObjects[i]->isDirectional()) {
58  
59 +      integrableObjects[i]->getGrad(dAtomGrad);
60 +
61        //gradient is equal to -f
62        grad[index++] = -dAtomGrad[0];
63        grad[index++] = -dAtomGrad[1];
# Line 69 | Line 68 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
68  
69      }
70      else{
71 <      atoms[i]->getFrc(force);
71 >      integrableObjects[i]->getFrc(force);
72  
73        grad[index++] = -force[0];
74        grad[index++] = -force[1];
# Line 98 | Line 97 | void OOPSEMinimizer::setCoor(vector<double>& x){
97  
98    index = 0;
99    
100 <  for(int i = 0; i < nAtoms; i++){
100 >  for(int i = 0; i < integrableObjects.size(); i++){
101      
102      position[0] = x[index++];
103      position[1] = x[index++];
104      position[2] = x[index++];
105  
106 <    atoms[i]->setPos(position);
106 >    integrableObjects[i]->setPos(position);
107  
108 <    if (atoms[i]->isDirectional()){
110 <      dAtom = (DirectionalAtom*) atoms[i];
108 >    if (integrableObjects[i]->isDirectional()){
109  
110        eulerAngle[0] = x[index++];
111        eulerAngle[1] = x[index++];
112        eulerAngle[2] = x[index++];
113  
114 <       dAtom->setEuler(eulerAngle[0], eulerAngle[1], eulerAngle[2]);
114 >      integrableObjects[i]->setEuler(eulerAngle[0],
115 >                                     eulerAngle[1],
116 >                                     eulerAngle[2]);
117  
118      }
119      
# Line 136 | Line 136 | vector<double> OOPSEMinimizer::getCoor(){
136  
137    index = 0;
138    
139 <  for(int i = 0; i < nAtoms; i++){
140 <    atoms[i]->getPos(position);
139 >  for(int i = 0; i < integrableObjects.size(); i++){
140 >    integrableObjects[i]->getPos(position);
141  
142      x[index++] = position[0];
143      x[index++] = position[1];
144      x[index++] = position[2];
145  
146 <    if (atoms[i]->isDirectional()){
147 <      dAtom = (DirectionalAtom*) atoms[i];
148 <      dAtom->getEulerAngles(eulerAngle);
146 >    if (integrableObjects[i]->isDirectional()){
147  
148 +      integrableObjects[i]->getEulerAngles(eulerAngle);
149 +      
150        x[index++] = eulerAngle[0];
151        x[index++] = eulerAngle[1];
152        x[index++] = eulerAngle[2];
# Line 424 | Line 424 | void OOPSEMinimizer::calcDim(){
424  
425    ndim = 0;
426  
427 <  for(int i = 0; i < nAtoms; i++){
427 >  for(int i = 0; i < integrableObjects.size(); i++){
428      ndim += 3;
429 <    if (atoms[i]->isDirectional())
429 >    if (integrableObjects[i]->isDirectional())
430        ndim += 3;      
431    }
432   }
# Line 483 | Line 483 | void OOPSEMinimizer::printMinimizerInfo(){
483  
484   /**
485   * In thoery, we need to find the minimum along the search direction
486 < * However, function evaluation is too expensive. I
486 > * However, function evaluation is too expensive.
487   * At the very begining of the problem, we check the search direction and make sure
488   * it is a descent direction
489   * we will compare the energy of two end points,
# Line 517 | Line 517 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
517    double mu;
518    double eta;
519    double ftol;  
520 +  double lsTol;
521  
522    xa.resize(ndim);
523    xb.resize(ndim);
# Line 532 | Line 533 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
533    ga = curG;
534    c = a + stepSize;
535    ftol = paramSet->getFTol();
536 +  lsTol = paramSet->getLineSearchTol();
537          
538    //calculate the derivative at a = 0
539    for (size_t i = 0; i < ndim; i++)
# Line 598 | Line 600 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
600        eta = 3 *(fa -fc) /(c - a) + slopeA + slopeC;
601        mu = sqrt(eta * eta - slopeA * slopeC);      
602        b = a + (c - a) * (1 - (slopeC + mu - eta) /(slopeC - slopeA + 2 * mu));      
603 +
604 +      if (b < lsTol){
605 +        if (bVerbose)
606 +          cout << "stepSize is less than line search tolerance" << endl;
607 +        break;        
608 +      }
609      //}
610  
611      // Take a trial step to this new point - new coords in xb

Diff Legend

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