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

Diff Legend

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