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 1819 by tim, Wed Dec 1 22:45:49 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 <
38 <
39 < VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo* info){
40 <
37 > VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info) {
38 >    dt2 = 0.5 * dt;
39   }
40  
41 < VelocityVerletIntegrator::~VelocityVerletIntegrator() {
41 > VelocityVerletIntegrator::~VelocityVerletIntegrator() { }
42  
43 < }
43 > void VelocityVerletIntegrator::initialize(){
44  
45 < void VelocityVerletIntegrator::integrate() {
45 >    currSample = 0.0;
46 >    currStatus = 0.0;
47 >    currThermal = 0.0;
48 >    needPotential = false;
49 >    needStress = false;      
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->writeFrame();
73 >    statWriter->writeStat(currentSnapshot_.statData);
74  
75 < void VelocityVerletIntegrator::moveB() {
60 <
75 >    
76   }
77 +        
78 + void VelocityVerletIntegrator::doIntegrate() {
79  
63 – void VelocityVerletIntegrator::thermalize() {
80  
81 < }
81 >    initialize();
82  
83 < void VelocityVerletIntegrator::velocitize() {
84 <  
85 <    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;
83 >    while (currentSnapshot_->getTime() < runTime) {
84 >        
85 >        preStep();
86  
87 <    std::vector<Molecule*>::iterator i;
81 <    std::vector<StuntDouble*>::iterator j;
82 <    Molecule* mol;
83 <    StuntDouble* integrableObject;
84 <    gaussianSPRNG gaussStream(info_->getSeed());
87 >        integrateStep();
88  
89 <    if (!info->have_target_temp) {
90 <        sprintf( painCave.errMsg,
91 <             "You can't resample the velocities without a targetTemp!\n"
92 <             );
90 <        painCave.isFatal = 1;
91 <        painCave.severity = OOPSE_ERROR;
92 <        simError();
93 <        return;
89 >        info_->incrTime(dt);
90 >        
91 >        postStep();
92 >
93      }
94  
95 <    temperature   = info_->target_temp;
95 >    finalize();
96 >    
97 > }
98  
98 –    kebar = kb * temperature * info_->getNdfRaw() / ( 2.0 * info_->getNdf());
99  
100 <    
101 <    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
100 > void VelocityVerletIntegrator::preStep() {
101 >        double difference = currentSnapshot_->getTime() + dt - currStatus;
102  
103 <            av2 = 2.0 * kebar / integrableObject->getMass();
104 <            vbar = sqrt( av2 );
103 >        if (difference > 0 || fabs(difference) < oopse::epsilon) {
104 >            needPotential = true;
105 >            needStress = true;  
106 >        }
107  
108 <            // picks random velocities from a gaussian distribution
111 <            // centered on vbar
108 > }
109  
110 <            for (int k=0; k<3; k++) {
111 <                aVel[k] = vbar * gaussStream.getGaussian();
110 > void VelocityVerletIntegrator::postStep() {
111 >        
112 >        
113 >        if (needVelocityScaling) {
114 >            if (currentSnapshot_->getTime() >= currThermal) {
115 >                velocitizer_->velocitize(targetScalingTemp);
116 >                currThermal += thermalTime;
117              }
118 <            
117 <            integrableObject->setVel( aVel );
118 >        }
119  
120 <            if(integrableObject->isDirectional()){
120 >        if (currentSnapshot_->getTime() >= currSample) {
121 >            dumpWriter->writeFrame(currentSnapshot_->getTime());
122 >            currSample += sampleTime;
123 >        }
124  
125 <                I = integrableObject->getI();
125 >        if (currentSnapshot_->getTime() >= currStatus) {
126 >            statWriter->writeStat(currentSnapshot_.statData);
127 >            needPotential = false;
128 >            needStress = false;
129 >            currStatus += statusTime;
130 >        }
131 > }
132  
123 –                if (integrableObject->isLinear()) {
133  
134 <                    l= integrableObject->linearAxis();
126 <                    m = (l+1)%3;
127 <                    n = (l+2)%3;
134 > void VelocityVerletIntegrator::finalize() {
135  
136 <                    aJ[l] = 0.0;
130 <                    vbar = sqrt( 2.0 * kebar * I(m, m) );
131 <                    aJ[m] = vbar * gaussStream.getGaussian();
132 <                    vbar = sqrt( 2.0 * kebar * I(n, n) );
133 <                    aJ[n] = vbar * gaussStream.getGaussian();
136 >    dumpWriter->writeFrame();
137  
138 <                } else {
138 >    delete dumpWriter;
139 >    delete statWriter;
140  
141 <                    for (int k = 0 ; k < 3; k++) {
142 <                        vbar = sqrt( 2.0 * kebar * I(k, k) );
139 <                        aJ[k] = vbar * gaussStream.getGaussian();
140 <                    }  
141 <
142 <                } // else isLinear
143 <
144 <                integrableObject->setJ( aJ );
145 <
146 <            }//isDirectional
147 <              
148 <       }
149 <    }//end for (mol = beginMolecule(i); ...)
150 <
151 <    // 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)) {
159 <
160 <            aVel = integrableObject->getVel();
161 <            aVel -= vdrift;
162 <            integrableObject->setVel( aVel );            
163 <        }
164 <    }  
141 >    dumpWriter = NULL;
142 >    statWriter = NULL;
143      
144   }
145  
146 < void VelocityVerletIntegrator::calcForce(int needPotential, int needStress){
146 > void VelocityVerletIntegrator::integrateStep() {
147  
148 +    moveA();
149 +    calcForce(needPotential, needStress);
150 +    moveB();
151   }
152  
172 – void VelocityVerletIntegrator::velocitize() {
153  
154 + void VelocityVerletIntegrator::calcForce(bool needPotential,
155 +                                         bool needStress) {
156 +    forceMan_->calcForces(needPotential, needStress)
157   }
158  
176 – void removeComDrift(){
159  
160 < }
160 > } //end namespace oopse

Diff Legend

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