OOPSE/libmdtools/OOPSEMinimizerBase.cpp

#include "Integrator.hpp"

OOPSEMinimizerBase::OOPSEMinimizerBase(SimInfo* theInfo, ForceFields* the_ff)
                           : RealIntegrator( theInfo, the_ff ){
  tStats = new Thermo(info);
  dumpOut = new DumpWriter(info);
  calcDim();
}

OOPSEMinimizerBase::~OOPSEMinimizerBase(){
  delete tStats;
  delete dumpOut;
}

/**
 *
 */


double OOPSEMinimizerBase::calcGradient(vector<double>& x, vector<double>& grad){
  Atom** atoms;  
  DirectionalAtom* dAtom;
  int index;
  double force[3];
  double dAtomGrad[6];

  setOptCoor(x);

  atoms = this->atoms;
  index = 0;

  for(int i = 0; i < nAtoms; i++){

    if(atoms[i]->isDirectional()){
      dAtom = (DirectionalAtom*) atoms[i];
      dAtom->getGrad(dAtomGrad);

      grad[index++] = dAtomGrad[0];
      grad[index++] = dAtomGrad[1];
      grad[index++] = dAtomGrad[2];
      grad[index++] = dAtomGrad[3];
      grad[index++] = dAtomGrad[4];
      grad[index++] = dAtomGrad[5];

    }
    else{
      atoms[i]->getFrc(force);

      grad[index++] = force[0];
      grad[index++] = force[1];
      grad[index++] = force[2];

    }
    
  }
 
  return tStats->getPotential();
  
}

/**
 *
 */

void OOPSEMinimizerBase::setOptCoor(vector<double>& x){
  Atom** atoms;  
  DirectionalAtom* dAtom;
  int index;
  double position[3];
  double eulerAngle[3];

  atoms = this->atoms;
  index = 0;
  
  for(int i = 0; i < nAtoms; i++){
    
    position[0] = x[index++];
    position[1] = x[index++];
    position[2] = x[index++];

    atoms[i]->setPos(position);

    if (atoms[i]->isDirectional()){
      dAtom = (DirectionalAtom*) atoms[i];
 
      eulerAngle[0] = x[index++];
      eulerAngle[1] = x[index++];
      eulerAngle[2] = x[index++];

       dAtom->setEuler(eulerAngle[0], eulerAngle[1], eulerAngle[2]);

    }
    
  }
  
}

/**
 *
 */
vector<double> OOPSEMinimizerBase::getOptCoor(){
  Atom** atoms;  
  DirectionalAtom* dAtom;
  int index;
  double position[3];
  double eulerAngle[3];
  vector<double> x;

  x.resize(getDim());
  atoms = this->atoms;
  index = 0;
  
  for(int i = 0; i < nAtoms; i++){
    atoms[i]->getPos(position);

    x[index++] = position[0];
    x[index++] = position[1];
    x[index++] = position[2];

    if (atoms[i]->isDirectional()){
      dAtom = (DirectionalAtom*) atoms[i];
      dAtom->getEulerAngles(eulerAngle);

      x[index++] = eulerAngle[0];
      x[index++] = eulerAngle[1];
      x[index++] = eulerAngle[2];
      
    }
    
  }

  return x;

}
         
void OOPSEMinimizerBase::calcDim(){
  Atom** atoms;  
  DirectionalAtom* dAtom;
  int dim;

  dim = 0;

  atoms = this->atoms;
  
  for(int i = 0; i < nAtoms; i++){
    dim += 3;
    if (atoms[i]->isDirectional())
      dim += 3;      
  }

}

void OOPSEMinimizerBase::output(vector<double>& x, int iteration){
  setOptCoor(x);
  dumpOut->writeDump(iteration);
}
Revision:	1031
Committed:	Fri Feb 6 18:58:06 2004 UTC (20 years, 5 months ago) by tim
File size:	2987 byte(s)
Log Message:	Add some lines into global.cpp to make it work with energy minimization
#	User	Rev	Content
1	tim	1011	#include "Integrator.hpp"
2
3			OOPSEMinimizerBase::OOPSEMinimizerBase(SimInfo* theInfo, ForceFields* the_ff)
4			: RealIntegrator( theInfo, the_ff ){
5			tStats = new Thermo(info);
6	tim	1031	dumpOut = new DumpWriter(info);
7			calcDim();
8	tim	1011	}
9
10			OOPSEMinimizerBase::~OOPSEMinimizerBase(){
11			delete tStats;
12	tim	1031	delete dumpOut;
13	tim	1011	}
14
15			/**
16			*
17			*/
18
19
20			double OOPSEMinimizerBase::calcGradient(vector<double>& x, vector<double>& grad){
21			Atom** atoms;
22			DirectionalAtom* dAtom;
23			int index;
24			double force[3];
25			double dAtomGrad[6];
26
27			setOptCoor(x);
28
29			atoms = this->atoms;
30			index = 0;
31
32			for(int i = 0; i < nAtoms; i++){
33
34			if(atoms[i]->isDirectional()){
35			dAtom = (DirectionalAtom*) atoms[i];
36			dAtom->getGrad(dAtomGrad);
37
38			grad[index++] = dAtomGrad[0];
39			grad[index++] = dAtomGrad[1];
40			grad[index++] = dAtomGrad[2];
41			grad[index++] = dAtomGrad[3];
42			grad[index++] = dAtomGrad[4];
43			grad[index++] = dAtomGrad[5];
44
45			}
46			else{
47			atoms[i]->getFrc(force);
48
49			grad[index++] = force[0];
50			grad[index++] = force[1];
51			grad[index++] = force[2];
52
53			}
54
55			}
56
57			return tStats->getPotential();
58
59			}
60
61			/**
62			*
63			*/
64
65			void OOPSEMinimizerBase::setOptCoor(vector<double>& x){
66			Atom** atoms;
67			DirectionalAtom* dAtom;
68			int index;
69			double position[3];
70			double eulerAngle[3];
71
72			atoms = this->atoms;
73			index = 0;
74
75			for(int i = 0; i < nAtoms; i++){
76
77			position[0] = x[index++];
78			position[1] = x[index++];
79			position[2] = x[index++];
80
81			atoms[i]->setPos(position);
82
83			if (atoms[i]->isDirectional()){
84			dAtom = (DirectionalAtom*) atoms[i];
85
86			eulerAngle[0] = x[index++];
87			eulerAngle[1] = x[index++];
88			eulerAngle[2] = x[index++];
89
90			dAtom->setEuler(eulerAngle[0], eulerAngle[1], eulerAngle[2]);
91
92			}
93
94			}
95
96			}
97
98			/**
99			*
100			*/
101	tim	1031	vector<double> OOPSEMinimizerBase::getOptCoor(){
102	tim	1011	Atom** atoms;
103			DirectionalAtom* dAtom;
104			int index;
105			double position[3];
106			double eulerAngle[3];
107	tim	1031	vector<double> x;
108	tim	1011
109	tim	1031	x.resize(getDim());
110	tim	1011	atoms = this->atoms;
111			index = 0;
112
113			for(int i = 0; i < nAtoms; i++){
114			atoms[i]->getPos(position);
115
116			x[index++] = position[0];
117			x[index++] = position[1];
118			x[index++] = position[2];
119
120			if (atoms[i]->isDirectional()){
121			dAtom = (DirectionalAtom*) atoms[i];
122			dAtom->getEulerAngles(eulerAngle);
123
124			x[index++] = eulerAngle[0];
125			x[index++] = eulerAngle[1];
126			x[index++] = eulerAngle[2];
127
128			}
129
130			}
131
132	tim	1031	return x;
133
134	tim	1011	}
135
136	tim	1031	void OOPSEMinimizerBase::calcDim(){
137			Atom** atoms;
138			DirectionalAtom* dAtom;
139			int dim;
140	tim	1011
141	tim	1031	dim = 0;
142
143			atoms = this->atoms;
144
145			for(int i = 0; i < nAtoms; i++){
146			dim += 3;
147			if (atoms[i]->isDirectional())
148			dim += 3;
149			}
150
151			}
152
153			void OOPSEMinimizerBase::output(vector<double>& x, int iteration){
154			setOptCoor(x);
155			dumpOut->writeDump(iteration);
156			}