OOPSE/libmdtools/OOPSEMinimizerBase.cpp

#include "Integrator.hpp"

OOPSEMinimizerBase::OOPSEMinimizerBase(SimInfo* theInfo, ForceFields* the_ff)
                           : RealIntegrator( theInfo, the_ff ){
  atoms = info->atoms;
  tStats = new Thermo(info);
  dumpOut = new DumpWriter(info);
  statOut = new StatWriter(info);
  calcDim();

  //save the reference coordinate
  RealIntegrator::preMove();
  
}

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

/**
 *
 */


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

  setCoor(x);

  if (nConstrained){
    shakeR();
  }
      
  calcForce(1, 1);
  
  if (nConstrained){
    shakeF();
  }
  
  x = getCoor();
  

  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){

  DirectionalAtom* dAtom;
  int index;
  double position[3];
  double eulerAngle[3];


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

  x.resize(getDim());

  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(){

  DirectionalAtom* dAtom;

  dim = 0;

  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);
}


//remove constraint force along the bond direction
void OOPSEMinimizerBase::shakeF(){
  int i, j;
  int done;
  double posA[3], posB[3];
  double frcA[3], frcB[3];
  double rab[3], fpab[3];
  int a, b, ax, ay, az, bx, by, bz;
  double rma, rmb;
  double rvab;
  double gab;
  double rabsq;
  double rfab;
  int iteration;

  for (i = 0; i < nAtoms; i++){
    moving[i] = 0;
    moved[i] = 1;
  }

  done = 0;
  iteration = 0;
  while (!done && (iteration < maxIteration)){
    done = 1;

    for (i = 0; i < nConstrained; i++){
      a = constrainedA[i];
      b = constrainedB[i];

      ax = (a * 3) + 0;
      ay = (a * 3) + 1;
      az = (a * 3) + 2;

      bx = (b * 3) + 0;
      by = (b * 3) + 1;
      bz = (b * 3) + 2;

      if (moved[a] || moved[b]){

        atoms[a]->getPos(posA);
        atoms[b]->getPos(posB);

        for (j = 0; j < 3; j++)
          rab[j] = posA[j] - posB[j];

        info->wrapVector(rab);

        atoms[a]->getFrc(frcA);
        atoms[b]->getFrc(frcB);

        //rma = 1.0 / atoms[a]->getMass();
        //rmb = 1.0 / atoms[b]->getMass();
        rma = 1.0;
        rmb = 1.0;
        
        
        fpab[0] = frcA[0] * rma - frcB[0] * rmb;
        fpab[1] = frcA[1] * rma - frcB[1] * rmb;
        fpab[2] = frcA[2] * rma - frcB[2] * rmb;


          gab=fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
          
          if (gab < 1.0)
            gab = 1.0;
         
          rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
          rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];

          if (fabs(rfab) > sqrt(rabsq*gab) * 0.00001){

            gab = -rfab /(rabsq*(rma + rmb));
            
            frcA[0] = rab[0] * gab;
            frcA[1] = rab[1] * gab;
            frcA[2] = rab[2] * gab;

            atoms[a]->addFrc(frcA);
            

            frcB[0] = -rab[0] * gab;
            frcB[1] = -rab[1] * gab;
            frcB[2] = -rab[2] * gab;

            atoms[b]->addFrc(frcB);
          
            moving[a] = 1;
            moving[b] = 1;
            done = 0;
          }             
      }
    }

    for (i = 0; i < nAtoms; i++){
      moved[i] = moving[i];
      moving[i] = 0;
    }

    iteration++;
  }

}

void OOPSEMinimizerBase::shakeR(){
  int i, j;
  int done;
  double posA[3], posB[3];
  double velA[3], velB[3];
  double pab[3];
  double rab[3];
  int a, b, ax, ay, az, bx, by, bz;
  double rma, rmb;
  double dx, dy, dz;
  double rpab;
  double rabsq, pabsq, rpabsq;
  double diffsq;
  double gab;
  int iteration;

  for (i = 0; i < nAtoms; i++){
    moving[i] = 0;
    moved[i] = 1;
  }

  iteration = 0;
  done = 0;
  while (!done && (iteration < 2*maxIteration)){
    done = 1;
    for (i = 0; i < nConstrained; i++){
      a = constrainedA[i];
      b = constrainedB[i];

      ax = (a * 3) + 0;
      ay = (a * 3) + 1;
      az = (a * 3) + 2;

      bx = (b * 3) + 0;
      by = (b * 3) + 1;
      bz = (b * 3) + 2;

      if (moved[a] || moved[b]){
        atoms[a]->getPos(posA);
        atoms[b]->getPos(posB);

        for (j = 0; j < 3; j++)
          pab[j] = posA[j] - posB[j];

        //periodic boundary condition

        info->wrapVector(pab);

        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];

        rabsq = constrainedDsqr[i];
        diffsq = rabsq - pabsq;

        // the original rattle code from alan tidesley
        if (fabs(diffsq) > (tol * rabsq * 2)){
          rab[0] = oldPos[ax] - oldPos[bx];
          rab[1] = oldPos[ay] - oldPos[by];
          rab[2] = oldPos[az] - oldPos[bz];

          info->wrapVector(rab);
          
          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];

          rpabsq = rpab * rpab;


          if (rpabsq < (rabsq * -diffsq)){
#ifdef IS_MPI
            a = atoms[a]->getGlobalIndex();
            b = atoms[b]->getGlobalIndex();
#endif //is_mpi
            cerr << "Waring: constraint failure" << endl;
            gab = sqrt(rabsq/pabsq);
            rab[0] = (posA[0] - posB[0])*gab;
            rab[1]= (posA[1] - posB[1])*gab;
            rab[2] = (posA[2] - posB[2])*gab;
            
            info->wrapVector(rab);
            
            rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
            
          }

          //rma = 1.0 / atoms[a]->getMass();
          //rmb = 1.0 / atoms[b]->getMass();
          rma = 1.0;
          rmb =1.0;

          gab = diffsq / (2.0 * (rma + rmb) * rpab);

          dx = rab[0] * gab;
          dy = rab[1] * gab;
          dz = rab[2] * gab;

          posA[0] += rma * dx;
          posA[1] += rma * dy;
          posA[2] += rma * dz;

          atoms[a]->setPos(posA);

          posB[0] -= rmb * dx;
          posB[1] -= rmb * dy;
          posB[2] -= rmb * dz;

          atoms[b]->setPos(posB);

          moving[a] = 1;
          moving[b] = 1;
          done = 0;
        }
      }
    }

    for (i = 0; i < nAtoms; i++){
      moved[i] = moving[i];
      moving[i] = 0;
    }

    iteration++;
  }

  if (!done)
     cerr << "Waring: can not constraint within maxIteration" << endl;
}
Revision:	1064
Committed:	Tue Feb 24 15:44:45 2004 UTC (21 years, 2 months ago) by tim
File size:	8121 byte(s)
Log Message:	Using inherit instead of compose to implement Minimizer both versions are working
#	Content
1	#include "Integrator.hpp"
2
3	OOPSEMinimizerBase::OOPSEMinimizerBase(SimInfo* theInfo, ForceFields* the_ff)
4	: RealIntegrator( theInfo, the_ff ){
5	atoms = info->atoms;
6	tStats = new Thermo(info);
7	dumpOut = new DumpWriter(info);
8	statOut = new StatWriter(info);
9	calcDim();
10
11	//save the reference coordinate
12	RealIntegrator::preMove();
13
14	}
15
16	OOPSEMinimizerBase::~OOPSEMinimizerBase(){
17	delete tStats;
18	delete dumpOut;
19	delete statOut;
20	}
21
22	/**
23	*
24	*/
25
26
27	double OOPSEMinimizerBase::calcGradient(vector<double>& x, vector<double>& grad){
28
29	DirectionalAtom* dAtom;
30	int index;
31	double force[3];
32	double dAtomGrad[6];
33
34	setCoor(x);
35
36	if (nConstrained){
37	shakeR();
38	}
39
40	calcForce(1, 1);
41
42	if (nConstrained){
43	shakeF();
44	}
45
46	x = getCoor();
47
48
49	index = 0;
50
51	for(int i = 0; i < nAtoms; i++){
52
53	if(atoms[i]->isDirectional()){
54	dAtom = (DirectionalAtom*) atoms[i];
55	dAtom->getGrad(dAtomGrad);
56
57	//gradient = du/dx = -f
58	grad[index++] = -dAtomGrad[0];
59	grad[index++] = -dAtomGrad[1];
60	grad[index++] = -dAtomGrad[2];
61	grad[index++] = -dAtomGrad[3];
62	grad[index++] = -dAtomGrad[4];
63	grad[index++] = -dAtomGrad[5];
64
65	}
66	else{
67	atoms[i]->getFrc(force);
68
69	grad[index++] = -force[0];
70	grad[index++] = -force[1];
71	grad[index++] = -force[2];
72
73	}
74
75	}
76
77	return tStats->getPotential();
78
79	}
80
81	/**
82	*
83	*/
84
85	void OOPSEMinimizerBase::setCoor(vector<double>& x){
86
87	DirectionalAtom* dAtom;
88	int index;
89	double position[3];
90	double eulerAngle[3];
91
92
93	index = 0;
94
95	for(int i = 0; i < nAtoms; i++){
96
97	position[0] = x[index++];
98	position[1] = x[index++];
99	position[2] = x[index++];
100
101	atoms[i]->setPos(position);
102
103	if (atoms[i]->isDirectional()){
104	dAtom = (DirectionalAtom*) atoms[i];
105
106	eulerAngle[0] = x[index++];
107	eulerAngle[1] = x[index++];
108	eulerAngle[2] = x[index++];
109
110	dAtom->setEuler(eulerAngle[0], eulerAngle[1], eulerAngle[2]);
111
112	}
113
114	}
115
116	}
117
118	/**
119	*
120	*/
121	vector<double> OOPSEMinimizerBase::getCoor(){
122
123	DirectionalAtom* dAtom;
124	int index;
125	double position[3];
126	double eulerAngle[3];
127	vector<double> x;
128
129	x.resize(getDim());
130
131	index = 0;
132
133	for(int i = 0; i < nAtoms; i++){
134	atoms[i]->getPos(position);
135
136	x[index++] = position[0];
137	x[index++] = position[1];
138	x[index++] = position[2];
139
140	if (atoms[i]->isDirectional()){
141	dAtom = (DirectionalAtom*) atoms[i];
142	dAtom->getEulerAngles(eulerAngle);
143
144	x[index++] = eulerAngle[0];
145	x[index++] = eulerAngle[1];
146	x[index++] = eulerAngle[2];
147
148	}
149
150	}
151
152	return x;
153
154	}
155
156	void OOPSEMinimizerBase::calcDim(){
157
158	DirectionalAtom* dAtom;
159
160	dim = 0;
161
162	for(int i = 0; i < nAtoms; i++){
163	dim += 3;
164	if (atoms[i]->isDirectional())
165	dim += 3;
166	}
167
168	}
169
170	void OOPSEMinimizerBase::output(vector<double>& x, int iteration){
171	setCoor(x);
172	calcForce(1, 1);
173	dumpOut->writeDump(iteration);
174	statOut->writeStat(iteration);
175	}
176
177
178	//remove constraint force along the bond direction
179	void OOPSEMinimizerBase::shakeF(){
180	int i, j;
181	int done;
182	double posA[3], posB[3];
183	double frcA[3], frcB[3];
184	double rab[3], fpab[3];
185	int a, b, ax, ay, az, bx, by, bz;
186	double rma, rmb;
187	double rvab;
188	double gab;
189	double rabsq;
190	double rfab;
191	int iteration;
192
193	for (i = 0; i < nAtoms; i++){
194	moving[i] = 0;
195	moved[i] = 1;
196	}
197
198	done = 0;
199	iteration = 0;
200	while (!done && (iteration < maxIteration)){
201	done = 1;
202
203	for (i = 0; i < nConstrained; i++){
204	a = constrainedA[i];
205	b = constrainedB[i];
206
207	ax = (a * 3) + 0;
208	ay = (a * 3) + 1;
209	az = (a * 3) + 2;
210
211	bx = (b * 3) + 0;
212	by = (b * 3) + 1;
213	bz = (b * 3) + 2;
214
215	if (moved[a] \|\| moved[b]){
216
217	atoms[a]->getPos(posA);
218	atoms[b]->getPos(posB);
219
220	for (j = 0; j < 3; j++)
221	rab[j] = posA[j] - posB[j];
222
223	info->wrapVector(rab);
224
225	atoms[a]->getFrc(frcA);
226	atoms[b]->getFrc(frcB);
227
228	//rma = 1.0 / atoms[a]->getMass();
229	//rmb = 1.0 / atoms[b]->getMass();
230	rma = 1.0;
231	rmb = 1.0;
232
233
234	fpab[0] = frcA[0] * rma - frcB[0] * rmb;
235	fpab[1] = frcA[1] * rma - frcB[1] * rmb;
236	fpab[2] = frcA[2] * rma - frcB[2] * rmb;
237
238
239	gab=fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
240
241	if (gab < 1.0)
242	gab = 1.0;
243
244	rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
245	rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];
246
247	if (fabs(rfab) > sqrt(rabsqgab) 0.00001){
248
249	gab = -rfab /(rabsq*(rma + rmb));
250
251	frcA[0] = rab[0] * gab;
252	frcA[1] = rab[1] * gab;
253	frcA[2] = rab[2] * gab;
254
255	atoms[a]->addFrc(frcA);
256
257
258	frcB[0] = -rab[0] * gab;
259	frcB[1] = -rab[1] * gab;
260	frcB[2] = -rab[2] * gab;
261
262	atoms[b]->addFrc(frcB);
263
264	moving[a] = 1;
265	moving[b] = 1;
266	done = 0;
267	}
268	}
269	}
270
271	for (i = 0; i < nAtoms; i++){
272	moved[i] = moving[i];
273	moving[i] = 0;
274	}
275
276	iteration++;
277	}
278
279	}
280
281	void OOPSEMinimizerBase::shakeR(){
282	int i, j;
283	int done;
284	double posA[3], posB[3];
285	double velA[3], velB[3];
286	double pab[3];
287	double rab[3];
288	int a, b, ax, ay, az, bx, by, bz;
289	double rma, rmb;
290	double dx, dy, dz;
291	double rpab;
292	double rabsq, pabsq, rpabsq;
293	double diffsq;
294	double gab;
295	int iteration;
296
297	for (i = 0; i < nAtoms; i++){
298	moving[i] = 0;
299	moved[i] = 1;
300	}
301
302	iteration = 0;
303	done = 0;
304	while (!done && (iteration < 2*maxIteration)){
305	done = 1;
306	for (i = 0; i < nConstrained; i++){
307	a = constrainedA[i];
308	b = constrainedB[i];
309
310	ax = (a * 3) + 0;
311	ay = (a * 3) + 1;
312	az = (a * 3) + 2;
313
314	bx = (b * 3) + 0;
315	by = (b * 3) + 1;
316	bz = (b * 3) + 2;
317
318	if (moved[a] \|\| moved[b]){
319	atoms[a]->getPos(posA);
320	atoms[b]->getPos(posB);
321
322	for (j = 0; j < 3; j++)
323	pab[j] = posA[j] - posB[j];
324
325	//periodic boundary condition
326
327	info->wrapVector(pab);
328
329	pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
330
331	rabsq = constrainedDsqr[i];
332	diffsq = rabsq - pabsq;
333
334	// the original rattle code from alan tidesley
335	if (fabs(diffsq) > (tol * rabsq * 2)){
336	rab[0] = oldPos[ax] - oldPos[bx];
337	rab[1] = oldPos[ay] - oldPos[by];
338	rab[2] = oldPos[az] - oldPos[bz];
339
340	info->wrapVector(rab);
341
342	rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
343
344	rpabsq = rpab * rpab;
345
346
347	if (rpabsq < (rabsq * -diffsq)){
348	#ifdef IS_MPI
349	a = atoms[a]->getGlobalIndex();
350	b = atoms[b]->getGlobalIndex();
351	#endif //is_mpi
352	cerr << "Waring: constraint failure" << endl;
353	gab = sqrt(rabsq/pabsq);
354	rab[0] = (posA[0] - posB[0])*gab;
355	rab[1]= (posA[1] - posB[1])*gab;
356	rab[2] = (posA[2] - posB[2])*gab;
357
358	info->wrapVector(rab);
359
360	rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
361
362	}
363
364	//rma = 1.0 / atoms[a]->getMass();
365	//rmb = 1.0 / atoms[b]->getMass();
366	rma = 1.0;
367	rmb =1.0;
368
369	gab = diffsq / (2.0 * (rma + rmb) * rpab);
370
371	dx = rab[0] * gab;
372	dy = rab[1] * gab;
373	dz = rab[2] * gab;
374
375	posA[0] += rma * dx;
376	posA[1] += rma * dy;
377	posA[2] += rma * dz;
378
379	atoms[a]->setPos(posA);
380
381	posB[0] -= rmb * dx;
382	posB[1] -= rmb * dy;
383	posB[2] -= rmb * dz;
384
385	atoms[b]->setPos(posB);
386
387	moving[a] = 1;
388	moving[b] = 1;
389	done = 0;
390	}
391	}
392	}
393
394	for (i = 0; i < nAtoms; i++){
395	moved[i] = moving[i];
396	moving[i] = 0;
397	}
398
399	iteration++;
400	}
401
402	if (!done)
403	cerr << "Waring: can not constraint within maxIteration" << endl;
404	}