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

Diff Legend

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