ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/integrators/VelocityVerletIntegrator.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/integrators/VelocityVerletIntegrator.cpp (file contents):
Revision 1722 by tim, Tue Nov 9 23:11:39 2004 UTC vs.
Revision 1774 by tim, Tue Nov 23 23:12:23 2004 UTC

# Line 22 | Line 22
22   * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
23   *
24   */
25 –
26 – /**
27 –  * @file VelocityVerletIntegrator.cpp
28 –  * @author tlin
29 –  * @date 11/09/2004
30 –  * @time 16:16am
31 –  * @version 1.0
32 –  */
25  
26 + /**
27 + * @file VelocityVerletIntegrator.cpp
28 + * @author tlin
29 + * @date 11/09/2004
30 + * @time 16:16am
31 + * @version 1.0
32 + */
33 +
34   #include "integrators/VelocityVerletIntegrator.hpp"
35  
36   namespace oopse {
37 + VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info) { }
38  
39 + VelocityVerletIntegrator::~VelocityVerletIntegrator() { }
40  
41 < VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo* info){
41 > void VelocityVerletIntegrator::integrate() {
42 >    double runTime = info_->run_time;
43  
44 < }
44 >    int calcPot;
45 >    int calcStress;
46  
43 – VelocityVerletIntegrator::~VelocityVerletIntegrator() {
47  
48 < }
48 >    dt_ = info_->dt;
49 >    dt2_ = 0.5 * dt_;
50  
51 < void VelocityVerletIntegrator::integrate() {
51 >    readyCheck();
52  
53 < }
53 >    // remove center of mass drift velocity (in case we passed in a configuration
54 >    // that was drifting
55 >    tStats->removeCOMdrift();
56  
57 < void VelocityVerletIntegrator::integrateStep() {
57 >    // initialize the forces before the first step
58  
59 < }
59 >    calcForce(1, 1);
60  
61 < void VelocityVerletIntegrator::moveA() {
61 >    //execute constraint algorithm to make sure at the very beginning the system is constrained  
62 >    if (nConstrained) {
63 >        constrainA();
64 >        calcForce(1, 1);
65 >        constrainB();
66 >    }
67  
68 < }
69 <
59 < void VelocityVerletIntegrator::moveB() {
60 <
61 < }
62 <
63 < void VelocityVerletIntegrator::thermalize() {
64 <
65 < }
66 <
67 < void VelocityVerletIntegrator::velocitize() {
68 <  
69 <    Vector3d aVel;
70 <    Vector3d aJ;
71 <    Mat3x3d I;
72 <    int l, m, n; // velocity randomizer loop counts
73 <    Vector3d vdrift;
74 <    double vbar;
75 <    const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
76 <    double av2;
77 <    double kebar;
78 <    double temperature;
79 <
80 <    std::vector<Molecule*>::iterator i;
81 <    std::vector<StuntDouble*>::iterator j;
82 <    Molecule* mol;
83 <    StuntDouble* integrableObject;
84 <    gaussianSPRNG gaussStream(info_->getSeed());
85 <
86 <    if (!info->have_target_temp) {
87 <        sprintf( painCave.errMsg,
88 <             "You can't resample the velocities without a targetTemp!\n"
89 <             );
90 <        painCave.isFatal = 1;
91 <        painCave.severity = OOPSE_ERROR;
92 <        simError();
93 <        return;
68 >    if (info_->setTemp) {
69 >        thermalize();
70      }
71  
72 <    temperature   = info_->target_temp;
72 >    calcPot = 0;
73 >    calcStress = 0;
74  
75 <    kebar = kb * temperature * info_->getNdfRaw() / ( 2.0 * info_->getNdf());
76 <
75 >    tStats = new Thermo(info_);
76 >    statOut = new StatWriter(info_);
77 >    dumpOut = new DumpWriter(info_);
78      
79 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
80 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
103 <               integrableObject = mol->nextIntegrableObject(j)) {
104 <    
105 <            // uses equipartition theory to solve for vbar in angstrom/fs
79 >    dumpOut->writeDump(info_->getTime());
80 >    statOut->writeStat(info_->getTime());
81  
82 <            av2 = 2.0 * kebar / integrableObject->getMass();
108 <            vbar = sqrt( av2 );
82 > #ifdef IS_MPI
83  
84 <            // picks random velocities from a gaussian distribution
85 <            // centered on vbar
84 >    strcpy(checkPointMsg, "The integrator is ready to go.");
85 >    MPIcheckPoint();
86  
87 <            for (int k=0; k<3; k++) {
114 <                aVel[k] = vbar * gaussStream.getGaussian();
115 <            }
116 <            
117 <            integrableObject->setVel( aVel );
87 > #endif // is_mpi
88  
89 <            if(integrableObject->isDirectional()){
89 >    while (info_->getTime() < runTime) {
90 >        difference = info_->getTime() + dt_ - currStatus;
91  
92 <                I = integrableObject->getI();
92 >        if (difference > 0 || fabs(difference) < 1e - 4) {
93 >            calcPot = 1;
94 >            calcStress = 1;
95 >        }
96  
97 <                if (integrableObject->isLinear()) {
97 >        //notify before integrateStep
98 >        notify()
99 >        integrateStep(calcPot, calcStress);
100 >        
101 >        info_->incrTime(dt_);
102 >        
103 >        //notify after integratreStep        
104 >        notify();
105 >        
106 >        if (info_->setTemp) {
107 >            if (info_->getTime() >= currThermal) {
108 >                thermalize();
109 >                currThermal += thermalTime;
110 >            }
111 >        }
112  
113 <                    l= integrableObject->linearAxis();
114 <                    m = (l+1)%3;
115 <                    n = (l+2)%3;
113 >        if (info_->getTime() >= currSample) {
114 >            dumpOut->writeDump(info_->getTime());
115 >            currSample += sampleTime;
116 >        }
117  
118 <                    aJ[l] = 0.0;
119 <                    vbar = sqrt( 2.0 * kebar * I(m, m) );
120 <                    aJ[m] = vbar * gaussStream.getGaussian();
121 <                    vbar = sqrt( 2.0 * kebar * I(n, n) );
122 <                    aJ[n] = vbar * gaussStream.getGaussian();
118 >        if (info_->getTime() >= currStatus) {
119 >            statOut->writeStat(info_->getTime());
120 >            calcPot = 0;
121 >            calcStress = 0;
122 >            currStatus += statusTime;
123 >        }
124  
125 <                } else {
125 >        if (info_->resetIntegrator) {
126 >            if (info_->getTime() >= currReset) {
127 >                this->resetIntegrator();
128 >                currReset += resetTime;
129 >            }
130 >        }
131  
132 <                    for (int k = 0 ; k < 3; k++) {
138 <                        vbar = sqrt( 2.0 * kebar * I(k, k) );
139 <                        aJ[k] = vbar * gaussStream.getGaussian();
140 <                    }  
132 > #ifdef IS_MPI
133  
134 <                } // else isLinear
134 >        strcpy(checkPointMsg, "successfully took a time step.");
135 >        MPIcheckPoint();
136  
137 <                integrableObject->setJ( aJ );
137 > #endif // is_mpi
138  
139 <            }//isDirectional
147 <              
148 <       }
149 <    }//end for (mol = beginMolecule(i); ...)
139 >    }
140  
141 <    // Get the Center of Mass drift velocity.
152 <    vdrift = info_->getComVel();
153 <  
154 <    //  Corrects for the center of mass drift.
155 <    // sums all the momentum and divides by total mass.
156 <    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
157 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
158 <               integrableObject = mol->nextIntegrableObject(j)) {
141 >    dumpOut->writeFinal(info_->getTime());
142  
143 <            aVel = integrableObject->getVel();
144 <            aVel -= vdrift;
162 <            integrableObject->setVel( aVel );            
163 <        }
164 <    }  
165 <    
143 >    delete dumpOut;
144 >    delete statOut;
145   }
146  
147 < void VelocityVerletIntegrator::calcForce(int needPotential, int needStress){
147 > void VelocityVerletIntegrator::integrateStep() {
148  
149 +    moveA();
150 +    calcForce(bool needPotential, bool needStress);
151 +    moveB();
152   }
153  
172 – void VelocityVerletIntegrator::velocitize() {
154  
155 < }
155 > void VelocityVerletIntegrator::thermalize() {
156  
157 < void removeComDrift(){
157 >    if (!info__->have_target_temp) {
158 >        sprintf(painCave.errMsg,
159 >                "You can't resample the velocities without a targetTemp!\n");
160 >        painCave.isFatal = 1;
161 >        painCave.severity = OOPSE_ERROR;
162 >        simError();
163 >        return;
164 >    }
165  
166   }
167 +
168 + void VelocityVerletIntegrator::calcForce(bool needPotential,
169 +                                         bool needStress) { }
170 +
171 +
172 + } //end namespace oopse

Diff Legend

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