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);
  statOut = new StatWriter(info);
  calcDim();
}

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

/**
 *
 */


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

  setCoor(x);
  calcForce(1, 1);
  
  atoms = info->atoms;
  index = 0;

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

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

      //gradient = du/dx = -f
      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::setCoor(vector<double>& x){
  Atom** atoms;  
  DirectionalAtom* dAtom;
  int index;
  double position[3];
  double eulerAngle[3];

  atoms = info->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::getCoor(){
  Atom** atoms;  
  DirectionalAtom* dAtom;
  int index;
  double position[3];
  double eulerAngle[3];
  vector<double> x;

  x.resize(getDim());
  atoms = info->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;

  dim = 0;

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

}

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