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 1074 by tim, Mon Mar 1 20:01:50 2004 UTC vs.
Revision 1234 by tim, Fri Jun 4 03:15:31 2004 UTC

# Line 1 | Line 1
1   #include <math.h>
2   #include "OOPSEMinimizer.hpp"
3 + #include "ShakeMin.hpp"
4  
5   OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, ForceFields* the_ff ,
6                                                MinimizerParameterSet * param)
7                       :RealIntegrator(theInfo, the_ff), bVerbose(false), bShake(true){
8 +  dumpOut = NULL;
9 +  statOut = NULL;
10  
8 –  atoms = info->atoms;
9 –
11    tStats = new Thermo(info);
12 <  dumpOut = new DumpWriter(info);
12 <  statOut = new StatWriter(info);
12 >
13    
14    paramSet = param;
15  
# Line 18 | Line 18 | OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, Forc
18    curX = getCoor();
19    curG.resize(ndim);
20  
21 <  preMove();
21 >  //preMove();
22   }
23  
24   OOPSEMinimizer::~OOPSEMinimizer(){
25    delete tStats;
26 <  delete dumpOut;
27 <  delete statOut;
26 >  if(dumpOut)
27 >    delete dumpOut;
28 >  if(statOut)
29 >    delete statOut;
30    delete paramSet;
31   }
32  
# Line 39 | Line 41 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
41    
42    setCoor(x);
43  
44 <  if (nConstrained && bShake){
45 <    shakeStatus = shakeR();
46 <  }
44 >  //if (nConstrained && bShake){
45 >  //  shakeStatus = shakeR();
46 >  //}
47 >
48 >  shakeAlgo->doShakeR();
49        
50    calcForce(1, 1);
51    
52 <  if (nConstrained && bShake){
53 <    shakeStatus |= shakeF();
54 <  }
52 >  //if (nConstrained && bShake){
53 >  //  shakeStatus |= shakeF();
54 >  //}
55 >
56 >  shakeAlgo->doShakeF();
57    
58    x = getCoor();
59    
60  
61    index = 0;
62  
63 <  for(int i = 0; i < nAtoms; i++){
63 >  for(int i = 0; i < integrableObjects.size(); i++){
64  
65 <    if(atoms[i]->isDirectional()){
60 <      dAtom = (DirectionalAtom*) atoms[i];
61 <      dAtom->getGrad(dAtomGrad);
65 >    if (integrableObjects[i]->isDirectional()) {
66  
67 +      integrableObjects[i]->getGrad(dAtomGrad);
68 +
69        //gradient is equal to -f
70        grad[index++] = -dAtomGrad[0];
71        grad[index++] = -dAtomGrad[1];
# Line 70 | Line 76 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
76  
77      }
78      else{
79 <      atoms[i]->getFrc(force);
79 >      integrableObjects[i]->getFrc(force);
80  
81        grad[index++] = -force[0];
82        grad[index++] = -force[1];
# Line 99 | Line 105 | void OOPSEMinimizer::setCoor(vector<double>& x){
105  
106    index = 0;
107    
108 <  for(int i = 0; i < nAtoms; i++){
108 >  for(int i = 0; i < integrableObjects.size(); i++){
109      
110      position[0] = x[index++];
111      position[1] = x[index++];
112      position[2] = x[index++];
113  
114 <    atoms[i]->setPos(position);
114 >    integrableObjects[i]->setPos(position);
115  
116 <    if (atoms[i]->isDirectional()){
111 <      dAtom = (DirectionalAtom*) atoms[i];
116 >    if (integrableObjects[i]->isDirectional()){
117  
118        eulerAngle[0] = x[index++];
119        eulerAngle[1] = x[index++];
120        eulerAngle[2] = x[index++];
121  
122 <       dAtom->setEuler(eulerAngle[0], eulerAngle[1], eulerAngle[2]);
122 >      integrableObjects[i]->setEuler(eulerAngle[0],
123 >                                     eulerAngle[1],
124 >                                     eulerAngle[2]);
125  
126      }
127      
# Line 137 | Line 144 | vector<double> OOPSEMinimizer::getCoor(){
144  
145    index = 0;
146    
147 <  for(int i = 0; i < nAtoms; i++){
148 <    atoms[i]->getPos(position);
147 >  for(int i = 0; i < integrableObjects.size(); i++){
148 >    integrableObjects[i]->getPos(position);
149  
150      x[index++] = position[0];
151      x[index++] = position[1];
152      x[index++] = position[2];
153  
154 <    if (atoms[i]->isDirectional()){
148 <      dAtom = (DirectionalAtom*) atoms[i];
149 <      dAtom->getEulerAngles(eulerAngle);
154 >    if (integrableObjects[i]->isDirectional()){
155  
156 +      integrableObjects[i]->getEulerAngles(eulerAngle);
157 +      
158        x[index++] = eulerAngle[0];
159        x[index++] = eulerAngle[1];
160        x[index++] = eulerAngle[2];
# Line 160 | Line 167 | vector<double> OOPSEMinimizer::getCoor(){
167  
168   }
169  
170 + /*
171   int OOPSEMinimizer::shakeR(){
172    int i, j;
173    int done;
# Line 399 | Line 407 | int OOPSEMinimizer::shakeF(){
407      return 1;
408   }
409  
410 <    //calculate the value of object function
410 > */
411 >
412 > //calculate the value of object function
413   void OOPSEMinimizer::calcF(){
414    calcEnergyGradient(curX, curG, curF, egEvalStatus);
415   }
# Line 425 | Line 435 | void OOPSEMinimizer::calcDim(){
435  
436    ndim = 0;
437  
438 <  for(int i = 0; i < nAtoms; i++){
438 >  for(int i = 0; i < integrableObjects.size(); i++){
439      ndim += 3;
440 <    if (atoms[i]->isDirectional())
440 >    if (integrableObjects[i]->isDirectional())
441        ndim += 3;      
442    }
443   }
# Line 484 | Line 494 | void OOPSEMinimizer::printMinimizerInfo(){
494  
495   /**
496   * In thoery, we need to find the minimum along the search direction
497 < * However, function evaluation is too expensive. I
497 > * However, function evaluation is too expensive.
498   * At the very begining of the problem, we check the search direction and make sure
499   * it is a descent direction
500   * we will compare the energy of two end points,
# Line 693 | Line 703 | void OOPSEMinimizer::minimize(){
703    
704    if (bVerbose)
705      printMinimizerInfo();
706 +
707 +  dumpOut = new DumpWriter(info);
708 +  statOut = new StatWriter(info);
709    
710    init();
711  

Diff Legend

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