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 (21 years, 9 months ago) by tim
File size:	2987 byte(s)
Log Message:	Add some lines into global.cpp to make it work with energy minimization
#	Content
1	#include "Integrator.hpp"
2
3	OOPSEMinimizerBase::OOPSEMinimizerBase(SimInfo* theInfo, ForceFields* the_ff)
4	: RealIntegrator( theInfo, the_ff ){
5	tStats = new Thermo(info);
6	dumpOut = new DumpWriter(info);
7	calcDim();
8	}
9
10	OOPSEMinimizerBase::~OOPSEMinimizerBase(){
11	delete tStats;
12	delete dumpOut;
13	}
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	vector<double> OOPSEMinimizerBase::getOptCoor(){
102	Atom** atoms;
103	DirectionalAtom* dAtom;
104	int index;
105	double position[3];
106	double eulerAngle[3];
107	vector<double> x;
108
109	x.resize(getDim());
110	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	return x;
133
134	}
135
136	void OOPSEMinimizerBase::calcDim(){
137	Atom** atoms;
138	DirectionalAtom* dAtom;
139	int dim;
140
141	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	}