ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new-templateless/OOPSE/libmdtools/Integrator.cpp
(Generate patch)

Comparing:
trunk/OOPSE/libmdtools/Integrator.cpp (file contents), Revision 784 by mmeineke, Wed Sep 24 19:34:39 2003 UTC vs.
branches/new-templateless/OOPSE/libmdtools/Integrator.cpp (file contents), Revision 852 by mmeineke, Thu Nov 6 18:20:47 2003 UTC

# Line 1 | Line 1
1   #include <iostream>
2 < #include <cstdlib>
3 < #include <cmath>
2 > #include <stdlib.h>
3 > #include <math.h>
4  
5   #ifdef IS_MPI
6   #include "mpiSimulation.hpp"
# Line 11 | Line 11
11   #include "simError.h"
12  
13  
14 < template<typename T> Integrator<T>::Integrator(SimInfo* theInfo,
15 <                                               ForceFields* the_ff){
14 > Integrator::Integrator(SimInfo* theInfo, ForceFields* the_ff){
15 >  
16    info = theInfo;
17    myFF = the_ff;
18    isFirst = 1;
# Line 25 | Line 25 | template<typename T> Integrator<T>::Integrator(SimInfo
25    if (info->the_integrator != NULL){
26      delete info->the_integrator;
27    }
28 <  
28 >
29    nAtoms = info->n_atoms;
30  
31    // check for constraints
# Line 42 | Line 42 | template<typename T> Integrator<T>::Integrator(SimInfo
42    checkConstraints();
43   }
44  
45 < template<typename T> Integrator<T>::~Integrator(){
45 > Integrator::~Integrator(){
46    if (nConstrained){
47      delete[] constrainedA;
48      delete[] constrainedB;
# Line 53 | Line 53 | template<typename T> Integrator<T>::~Integrator(){
53    }
54   }
55  
56 < template<typename T> void Integrator<T>::checkConstraints(void){
56 > void Integrator::checkConstraints(void){
57    isConstrained = 0;
58  
59    Constraint* temp_con;
# Line 145 | Line 145 | template<typename T> void Integrator<T>::checkConstrai
145   }
146  
147  
148 < template<typename T> void Integrator<T>::integrate(void){
149 <  int i, j;                         // loop counters
148 > void Integrator::integrate(void){
149  
150    double runTime = info->run_time;
151    double sampleTime = info->sampleTime;
# Line 154 | Line 153 | template<typename T> void Integrator<T>::integrate(voi
153    double thermalTime = info->thermalTime;
154    double resetTime = info->resetTime;
155  
156 +  
157  
158    double currSample;
159    double currThermal;
160    double currStatus;
161    double currReset;
162 <  
162 >
163    int calcPot, calcStress;
164 –  int isError;
164  
165    tStats = new Thermo(info);
166    statOut = new StatWriter(info);
167    dumpOut = new DumpWriter(info);
168  
169    atoms = info->atoms;
171 –  DirectionalAtom* dAtom;
170  
171    dt = info->dt;
172    dt2 = 0.5 * dt;
# Line 182 | Line 180 | template<typename T> void Integrator<T>::integrate(voi
180    if (nConstrained){
181      preMove();
182      constrainA();
183 <    calcForce(1, 1);    
183 >    calcForce(1, 1);
184      constrainB();
185    }
186    
# Line 201 | Line 199 | template<typename T> void Integrator<T>::integrate(voi
199    statOut->writeStat(info->getTime());
200  
201  
204 –
202   #ifdef IS_MPI
203    strcpy(checkPointMsg, "The integrator is ready to go.");
204    MPIcheckPoint();
# Line 230 | Line 227 | template<typename T> void Integrator<T>::integrate(voi
227      }
228  
229      if (info->getTime() >= currStatus){
230 <      statOut->writeStat(info->getTime());
231 <      calcPot = 0;
230 >      statOut->writeStat(info->getTime());
231 >      calcPot = 0;
232        calcStress = 0;
233        currStatus += statusTime;
234 <    }
234 >    }
235  
236      if (info->resetIntegrator){
237        if (info->getTime() >= currReset){
# Line 249 | Line 246 | template<typename T> void Integrator<T>::integrate(voi
246   #endif // is_mpi
247    }
248  
252 –  dumpOut->writeFinal(info->getTime());
249  
250 +  // write the last frame
251 +  dumpOut->writeDump(info->getTime());
252 +
253    delete dumpOut;
254    delete statOut;
255   }
256  
257 < template<typename T> void Integrator<T>::integrateStep(int calcPot,
258 <                                                       int calcStress){
257 > void Integrator::integrateStep(int calcPot,
258 >                               int calcStress){
259    // Position full step, and velocity half step
260    preMove();
261  
# Line 294 | Line 293 | template<typename T> void Integrator<T>::integrateStep
293   }
294  
295  
296 < template<typename T> void Integrator<T>::moveA(void){
296 > void Integrator::moveA(void){
297    int i, j;
298    DirectionalAtom* dAtom;
299    double Tb[3], ji[3];
# Line 345 | Line 344 | template<typename T> void Integrator<T>::moveA(void){
344   }
345  
346  
347 < template<typename T> void Integrator<T>::moveB(void){
347 > void Integrator::moveB(void){
348    int i, j;
349    DirectionalAtom* dAtom;
350    double Tb[3], ji[3];
# Line 367 | Line 366 | template<typename T> void Integrator<T>::moveB(void){
366      if (atoms[i]->isDirectional()){
367        dAtom = (DirectionalAtom *) atoms[i];
368  
369 <      // get and convert the torque to body frame      
369 >      // get and convert the torque to body frame
370  
371        dAtom->getTrq(Tb);
372        dAtom->lab2Body(Tb);
# Line 389 | Line 388 | template<typename T> void Integrator<T>::moveB(void){
388    }
389   }
390  
391 < template<typename T> void Integrator<T>::preMove(void){
391 > void Integrator::preMove(void){
392    int i, j;
393    double pos[3];
394  
# Line 404 | Line 403 | template<typename T> void Integrator<T>::preMove(void)
403    }
404   }
405  
406 < template<typename T> void Integrator<T>::constrainA(){
407 <  int i, j, k;
406 > void Integrator::constrainA(){
407 >  int i, j;
408    int done;
409    double posA[3], posB[3];
410    double velA[3], velB[3];
# Line 548 | Line 547 | template<typename T> void Integrator<T>::constrainA(){
547  
548   }
549  
550 < template<typename T> void Integrator<T>::constrainB(void){
551 <  int i, j, k;
550 > void Integrator::constrainB(void){
551 >  int i, j;
552    int done;
553    double posA[3], posB[3];
554    double velA[3], velB[3];
# Line 558 | Line 557 | template<typename T> void Integrator<T>::constrainB(vo
557    int a, b, ax, ay, az, bx, by, bz;
558    double rma, rmb;
559    double dx, dy, dz;
560 <  double rabsq, pabsq, rvab;
562 <  double diffsq;
560 >  double rvab;
561    double gab;
562    int iteration;
563  
# Line 649 | Line 647 | template<typename T> void Integrator<T>::constrainB(vo
647    }
648   }
649  
650 < template<typename T> void Integrator<T>::rotationPropagation
650 > void Integrator::rotationPropagation
651   ( DirectionalAtom* dAtom, double ji[3] ){
652  
653    double angle;
# Line 660 | Line 658 | template<typename T> void Integrator<T>::rotationPropa
658  
659    dAtom->getA(A);
660    dAtom->getI(I);
661 <  
662 <  // rotate about the x-axis      
661 >
662 >  // rotate about the x-axis
663    angle = dt2 * ji[0] / I[0][0];
664 <  this->rotate( 1, 2, angle, ji, A );
665 <  
664 >  this->rotate( 1, 2, angle, ji, A );
665 >
666    // rotate about the y-axis
667    angle = dt2 * ji[1] / I[1][1];
668    this->rotate( 2, 0, angle, ji, A );
669 <  
669 >
670    // rotate about the z-axis
671    angle = dt * ji[2] / I[2][2];
672    this->rotate( 0, 1, angle, ji, A);
673 <  
673 >
674    // rotate about the y-axis
675    angle = dt2 * ji[1] / I[1][1];
676    this->rotate( 2, 0, angle, ji, A );
677 <  
677 >
678    // rotate about the x-axis
679    angle = dt2 * ji[0] / I[0][0];
680    this->rotate( 1, 2, angle, ji, A );
681 <  
682 <  dAtom->setA( A  );    
681 >
682 >  dAtom->setA( A  );
683   }
684  
685 < template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
686 <                                                double angle, double ji[3],
687 <                                                double A[3][3]){
685 > void Integrator::rotate(int axes1, int axes2,
686 >                        double angle, double ji[3],
687 >                        double A[3][3]){
688    int i, j, k;
689    double sinAngle;
690    double cosAngle;
# Line 749 | Line 747 | template<typename T> void Integrator<T>::rotate(int ax
747      }
748    }
749  
750 <  // rotate the Rotation matrix acording to:
750 >  // rotate the Rotation matrix acording to:
751    //            A[][] = A[][] * transpose(rot[][])
752  
753  
# Line 767 | Line 765 | template<typename T> void Integrator<T>::rotate(int ax
765    }
766   }
767  
768 < template<typename T> void Integrator<T>::calcForce(int calcPot, int calcStress){
768 > void Integrator::calcForce(int calcPot, int calcStress){
769    myFF->doForces(calcPot, calcStress);
770   }
771  
772 < template<typename T> void Integrator<T>::thermalize(){
772 > void Integrator::thermalize(){
773    tStats->velocitize();
774   }
775  
776 < template<typename T> double Integrator<T>::getConservedQuantity(void){
776 > double Integrator::getConservedQuantity(void){
777    return tStats->getTotalE();
778   }
779 +
780 + char* Integrator::getAdditionalParameters(void){
781 +  //By default, return a null string
782 +  //The reason we use string instead of char* is that if we use char*, we will
783 +  //return a pointer point to local variable which might cause problem
784 +  addParams[0] = '\0';
785 +  return addParams;
786 + }

Diff Legend

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