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 1901 by tim, Tue Jan 4 22:18:36 2005 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 + #include "integrators/DLM.hpp"
36  
37   namespace oopse {
38 <
39 <
40 < VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo* info){
41 <
38 > VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info), rotAlgo(NULL) {
39 >    dt2 = 0.5 * dt;
40 >    rotAlgo = new DLM();
41 >    rattle = new Rattle(info);
42   }
43  
44 < VelocityVerletIntegrator::~VelocityVerletIntegrator() {
45 <
44 > VelocityVerletIntegrator::~VelocityVerletIntegrator() {
45 >    delete rotAlgo;
46 >    delete rattle;
47   }
48  
49 < void VelocityVerletIntegrator::integrate() {
49 > void VelocityVerletIntegrator::initialize(){
50  
51 < }
51 >    // remove center of mass drift velocity (in case we passed in a configuration
52 >    // that was drifting
53 >    velocitizer_->removeComDrift();
54  
55 < void VelocityVerletIntegrator::integrateStep() {
55 >    // initialize the forces before the first step
56 >    calcForce(true, true);
57  
58 < }
58 >    //execute constraint algorithm to make sure at the very beginning the system is constrained  
59 >    //if (nConstrained) {
60 >    //    constrainA();
61 >    //    calcForce(true, true);
62 >    //    constrainB();
63 >    //}
64  
65 < void VelocityVerletIntegrator::moveA() {
65 >    if (needVelocityScaling) {
66 >        velocitizer_->velocitize(targetScalingTemp);
67 >    }
68 >    
69 >    dumpWriter = createDumpWriter();
70 >    statWriter = createStatWriter();
71  
72 < }
72 >    dumpWriter->writeDump();
73  
74 < void VelocityVerletIntegrator::moveB() {
74 >    //save statistics, before writeStat,  we must save statistics
75 >    thermo.saveStat();
76 >    saveConservedQuantity();
77 >    statWriter->writeStat(currentSnapshot_->statData);
78  
79 +    currSample = sampleTime + currentSnapshot_->getTime();
80 +    currStatus =  statusTime + currentSnapshot_->getTime();;
81 +    currThermal = thermalTime +  + currentSnapshot_->getTime();
82 +    needPotential = false;
83 +    needStress = false;      
84 +  
85   }
86 +        
87 + void VelocityVerletIntegrator::doIntegrate() {
88  
63 – void VelocityVerletIntegrator::thermalize() {
89  
90 < }
90 >    initialize();
91  
92 < void VelocityVerletIntegrator::velocitize() {
93 <  
94 <    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;
92 >    while (currentSnapshot_->getTime() < runTime) {
93 >        
94 >        preStep();
95  
96 <    std::vector<Molecule*>::iterator i;
97 <    std::vector<StuntDouble*>::iterator j;
98 <    Molecule* mol;
83 <    StuntDouble* integrableObject;
84 <    gaussianSPRNG gaussStream(info_->getSeed());
96 >        integrateStep();
97 >        
98 >        postStep();
99  
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;
100      }
101  
102 <    temperature   = info_->target_temp;
97 <
98 <    kebar = kb * temperature * info_->getNdfRaw() / ( 2.0 * info_->getNdf());
99 <
102 >    finalize();
103      
104 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
102 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
103 <               integrableObject = mol->nextIntegrableObject(j)) {
104 <    
105 <            // uses equipartition theory to solve for vbar in angstrom/fs
104 > }
105  
107 –            av2 = 2.0 * kebar / integrableObject->getMass();
108 –            vbar = sqrt( av2 );
106  
107 <            // picks random velocities from a gaussian distribution
108 <            // centered on vbar
107 > void VelocityVerletIntegrator::preStep() {
108 >        double difference = currentSnapshot_->getTime() + dt - currStatus;
109  
110 <            for (int k=0; k<3; k++) {
111 <                aVel[k] = vbar * gaussStream.getGaussian();
112 <            }
113 <            
117 <            integrableObject->setVel( aVel );
110 >        if (difference > 0 || fabs(difference) < oopse::epsilon) {
111 >            needPotential = true;
112 >            needStress = true;  
113 >        }
114  
115 <            if(integrableObject->isDirectional()){
115 > }
116  
117 <                I = integrableObject->getI();
117 > void VelocityVerletIntegrator::postStep() {
118  
119 <                if (integrableObject->isLinear()) {
119 >        //save snapshot
120 >        info_->getSnapshotManager()->advance();
121  
122 <                    l= integrableObject->linearAxis();
123 <                    m = (l+1)%3;
124 <                    n = (l+2)%3;
122 >        //increase time
123 >        currentSnapshot_->increaseTime(dt);        
124 >        
125 >        if (needVelocityScaling) {
126 >            if (currentSnapshot_->getTime() >= currThermal) {
127 >                velocitizer_->velocitize(targetScalingTemp);
128 >                currThermal += thermalTime;
129 >            }
130 >        }
131  
132 <                    aJ[l] = 0.0;
133 <                    vbar = sqrt( 2.0 * kebar * I(m, m) );
134 <                    aJ[m] = vbar * gaussStream.getGaussian();
135 <                    vbar = sqrt( 2.0 * kebar * I(n, n) );
133 <                    aJ[n] = vbar * gaussStream.getGaussian();
132 >        if (currentSnapshot_->getTime() >= currSample) {
133 >            dumpWriter->writeDump();
134 >            currSample += sampleTime;
135 >        }
136  
137 <                } else {
137 >        if (currentSnapshot_->getTime() >= currStatus) {
138 >            //save statistics, before writeStat,  we must save statistics
139 >            thermo.saveStat();
140 >            saveConservedQuantity();
141 >            statWriter->writeStat(currentSnapshot_->statData);
142  
143 <                    for (int k = 0 ; k < 3; k++) {
144 <                        vbar = sqrt( 2.0 * kebar * I(k, k) );
145 <                        aJ[k] = vbar * gaussStream.getGaussian();
146 <                    }  
143 >            needPotential = false;
144 >            needStress = false;
145 >            currStatus += statusTime;
146 >        }
147  
148 <                } // else isLinear
148 >        
149 > }
150  
144 –                integrableObject->setJ( aJ );
151  
152 <            }//isDirectional
147 <              
148 <       }
149 <    }//end for (mol = beginMolecule(i); ...)
152 > void VelocityVerletIntegrator::finalize() {
153  
154 <    // 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)) {
154 >    dumpWriter->writeDump();
155  
156 <            aVel = integrableObject->getVel();
157 <            aVel -= vdrift;
158 <            integrableObject->setVel( aVel );            
159 <        }
160 <    }  
156 >    delete dumpWriter;
157 >    delete statWriter;
158 >
159 >    dumpWriter = NULL;
160 >    statWriter = NULL;
161      
162   }
163  
164 < void VelocityVerletIntegrator::calcForce(int needPotential, int needStress){
164 > void VelocityVerletIntegrator::integrateStep() {
165  
166 +    moveA();
167 +    calcForce(needPotential, needStress);
168 +    moveB();
169   }
170  
172 – void VelocityVerletIntegrator::velocitize() {
171  
172 + void VelocityVerletIntegrator::calcForce(bool needPotential,
173 +                                         bool needStress) {
174 +    forceMan_->calcForces(needPotential, needStress);
175   }
176  
177 < void removeComDrift(){
177 > DumpWriter* VelocityVerletIntegrator::createDumpWriter() {
178 >    return new DumpWriter(info_, info_->getDumpFileName());
179 > }
180  
181 + StatWriter* VelocityVerletIntegrator::createStatWriter() {
182 +    return new StatWriter(info_->getStatFileName());
183   }
184 +
185 +
186 + } //end namespace oopse

Diff Legend

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