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 1144 by tim, Sat May 1 18:52:38 2004 UTC

# Line 4 | Line 4 | OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, Forc
4   OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, ForceFields* the_ff ,
5                                                MinimizerParameterSet * param)
6                       :RealIntegrator(theInfo, the_ff), bVerbose(false), bShake(true){
7 +  dumpOut = NULL;
8 +  statOut = NULL;
9  
8 –  atoms = info->atoms;
9 –
10    tStats = new Thermo(info);
11 <  dumpOut = new DumpWriter(info);
12 <  statOut = new StatWriter(info);
11 >
12    
13    paramSet = param;
14  
# Line 18 | Line 17 | OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, Forc
17    curX = getCoor();
18    curG.resize(ndim);
19  
20 +  preMove();
21   }
22  
23   OOPSEMinimizer::~OOPSEMinimizer(){
24    delete tStats;
25 <  delete dumpOut;
26 <  delete statOut;
25 >  if(dumpOut)
26 >    delete dumpOut;
27 >  if(statOut)
28 >    delete statOut;
29    delete paramSet;
30   }
31  
# Line 53 | Line 55 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
55  
56    index = 0;
57  
58 <  for(int i = 0; i < nAtoms; i++){
58 >  for(int i = 0; i < integrableObjects.size(); i++){
59  
60 <    if(atoms[i]->isDirectional()){
59 <      dAtom = (DirectionalAtom*) atoms[i];
60 <      dAtom->getGrad(dAtomGrad);
60 >    if (integrableObjects[i]->isDirectional()) {
61  
62 +      integrableObjects[i]->getGrad(dAtomGrad);
63 +
64        //gradient is equal to -f
65        grad[index++] = -dAtomGrad[0];
66        grad[index++] = -dAtomGrad[1];
# Line 69 | Line 71 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
71  
72      }
73      else{
74 <      atoms[i]->getFrc(force);
74 >      integrableObjects[i]->getFrc(force);
75  
76        grad[index++] = -force[0];
77        grad[index++] = -force[1];
# Line 98 | Line 100 | void OOPSEMinimizer::setCoor(vector<double>& x){
100  
101    index = 0;
102    
103 <  for(int i = 0; i < nAtoms; i++){
103 >  for(int i = 0; i < integrableObjects.size(); i++){
104      
105      position[0] = x[index++];
106      position[1] = x[index++];
107      position[2] = x[index++];
108  
109 <    atoms[i]->setPos(position);
109 >    integrableObjects[i]->setPos(position);
110  
111 <    if (atoms[i]->isDirectional()){
110 <      dAtom = (DirectionalAtom*) atoms[i];
111 >    if (integrableObjects[i]->isDirectional()){
112  
113        eulerAngle[0] = x[index++];
114        eulerAngle[1] = x[index++];
115        eulerAngle[2] = x[index++];
116  
117 <       dAtom->setEuler(eulerAngle[0], eulerAngle[1], eulerAngle[2]);
117 >      integrableObjects[i]->setEuler(eulerAngle[0],
118 >                                     eulerAngle[1],
119 >                                     eulerAngle[2]);
120  
121      }
122      
# Line 136 | Line 139 | vector<double> OOPSEMinimizer::getCoor(){
139  
140    index = 0;
141    
142 <  for(int i = 0; i < nAtoms; i++){
143 <    atoms[i]->getPos(position);
142 >  for(int i = 0; i < integrableObjects.size(); i++){
143 >    integrableObjects[i]->getPos(position);
144  
145      x[index++] = position[0];
146      x[index++] = position[1];
147      x[index++] = position[2];
148  
149 <    if (atoms[i]->isDirectional()){
147 <      dAtom = (DirectionalAtom*) atoms[i];
148 <      dAtom->getEulerAngles(eulerAngle);
149 >    if (integrableObjects[i]->isDirectional()){
150  
151 +      integrableObjects[i]->getEulerAngles(eulerAngle);
152 +      
153        x[index++] = eulerAngle[0];
154        x[index++] = eulerAngle[1];
155        x[index++] = eulerAngle[2];
# Line 424 | Line 427 | void OOPSEMinimizer::calcDim(){
427  
428    ndim = 0;
429  
430 <  for(int i = 0; i < nAtoms; i++){
430 >  for(int i = 0; i < integrableObjects.size(); i++){
431      ndim += 3;
432 <    if (atoms[i]->isDirectional())
432 >    if (integrableObjects[i]->isDirectional())
433        ndim += 3;      
434    }
435   }
# Line 483 | Line 486 | void OOPSEMinimizer::printMinimizerInfo(){
486  
487   /**
488   * In thoery, we need to find the minimum along the search direction
489 < * However, function evaluation is too expensive. I
489 > * However, function evaluation is too expensive.
490   * At the very begining of the problem, we check the search direction and make sure
491   * it is a descent direction
492   * we will compare the energy of two end points,
# Line 517 | Line 520 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
520    double mu;
521    double eta;
522    double ftol;  
523 +  double lsTol;
524  
525    xa.resize(ndim);
526    xb.resize(ndim);
# Line 532 | Line 536 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
536    ga = curG;
537    c = a + stepSize;
538    ftol = paramSet->getFTol();
539 +  lsTol = paramSet->getLineSearchTol();
540          
541    //calculate the derivative at a = 0
542    for (size_t i = 0; i < ndim; i++)
# Line 598 | Line 603 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
603        eta = 3 *(fa -fc) /(c - a) + slopeA + slopeC;
604        mu = sqrt(eta * eta - slopeA * slopeC);      
605        b = a + (c - a) * (1 - (slopeC + mu - eta) /(slopeC - slopeA + 2 * mu));      
606 +
607 +      if (b < lsTol){
608 +        if (bVerbose)
609 +          cout << "stepSize is less than line search tolerance" << endl;
610 +        break;        
611 +      }
612      //}
613  
614      // Take a trial step to this new point - new coords in xb
# Line 684 | Line 695 | void OOPSEMinimizer::minimize(){
695    
696    if (bVerbose)
697      printMinimizerInfo();
698 +
699 +  dumpOut = new DumpWriter(info);
700 +  statOut = new StatWriter(info);
701    
702    init();
703  

Diff Legend

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