src/integrators/VelocityVerletIntegrator.cpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */
 
 /**
  * @file VelocityVerletIntegrator.cpp
  * @author tlin
  * @date 11/09/2004
  * @time 16:16am
  * @version 1.0
  */

#include "integrators/VelocityVerletIntegrator.hpp"

namespace oopse {


VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo* info){ 

}

VelocityVerletIntegrator::~VelocityVerletIntegrator() {

}

void VelocityVerletIntegrator::integrate() { 

}

void VelocityVerletIntegrator::integrateStep() { 

}

void VelocityVerletIntegrator::moveA() { 

}

void VelocityVerletIntegrator::moveB() {

}

void VelocityVerletIntegrator::thermalize() {

}

void VelocityVerletIntegrator::velocitize() { 
  
    Vector3d aVel;
    Vector3d aJ;
    Mat3x3d I;
    int l, m, n; // velocity randomizer loop counts
    Vector3d vdrift;
    double vbar;
    const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
    double av2;
    double kebar;
    double temperature;

    std::vector<Molecule*>::iterator i;
    std::vector<StuntDouble*>::iterator j;
    Molecule* mol;
    StuntDouble* integrableObject;
    gaussianSPRNG gaussStream(info_->getSeed());

    if (!info->have_target_temp) {
        sprintf( painCave.errMsg,
             "You can't resample the velocities without a targetTemp!\n"
             );
        painCave.isFatal = 1;
        painCave.severity = OOPSE_ERROR;
        simError();
        return;
    }

    temperature   = info_->target_temp;

    kebar = kb * temperature * info_->getNdfRaw() / ( 2.0 * info_->getNdf());

    
    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
               integrableObject = mol->nextIntegrableObject(j)) {
    
            // uses equipartition theory to solve for vbar in angstrom/fs

            av2 = 2.0 * kebar / integrableObject->getMass();
            vbar = sqrt( av2 );

            // picks random velocities from a gaussian distribution
            // centered on vbar

            for (int k=0; k<3; k++) {
                aVel[k] = vbar * gaussStream.getGaussian();
            }
            
            integrableObject->setVel( aVel );

            if(integrableObject->isDirectional()){

                I = integrableObject->getI();

                if (integrableObject->isLinear()) {

                    l= integrableObject->linearAxis();
                    m = (l+1)%3;
                    n = (l+2)%3;

                    aJ[l] = 0.0;
                    vbar = sqrt( 2.0 * kebar * I(m, m) );
                    aJ[m] = vbar * gaussStream.getGaussian();
                    vbar = sqrt( 2.0 * kebar * I(n, n) );
                    aJ[n] = vbar * gaussStream.getGaussian();

                } else {

                    for (int k = 0 ; k < 3; k++) {
                        vbar = sqrt( 2.0 * kebar * I(k, k) );
                        aJ[k] = vbar * gaussStream.getGaussian();
                    }   

                } // else isLinear

                integrableObject->setJ( aJ ); 

            }//isDirectional 
              
       }
    }//end for (mol = beginMolecule(i); ...)

    // Get the Center of Mass drift velocity.
    vdrift = info_->getComVel();
  
    //  Corrects for the center of mass drift.
    // sums all the momentum and divides by total mass.
    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
               integrableObject = mol->nextIntegrableObject(j)) {

            aVel = integrableObject->getVel();
            aVel -= vdrift;
            integrableObject->setVel( aVel );            
        }
    }  
    
}

void VelocityVerletIntegrator::calcForce(int needPotential, int needStress){ 

}

void VelocityVerletIntegrator::velocitize() { 

}

void removeComDrift(){ 

}
Revision:	1722
Committed:	Tue Nov 9 23:11:39 2004 UTC (19 years, 8 months ago) by tim
File size:	5004 byte(s)
Log Message:	adding ForceManager
#	User	Rev	Content
1	tim	1722	/*
2			* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3			*
4			* Contact: oopse@oopse.org
5			*
6			* This program is free software; you can redistribute it and/or
7			* modify it under the terms of the GNU Lesser General Public License
8			* as published by the Free Software Foundation; either version 2.1
9			* of the License, or (at your option) any later version.
10			* All we ask is that proper credit is given for our work, which includes
11			* - but is not limited to - adding the above copyright notice to the beginning
12			* of your source code files, and to any copyright notice that you may distribute
13			* with programs based on this work.
14			*
15			* This program is distributed in the hope that it will be useful,
16			* but WITHOUT ANY WARRANTY; without even the implied warranty of
17			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18			* GNU Lesser General Public License for more details.
19			*
20			* You should have received a copy of the GNU Lesser General Public License
21			* along with this program; if not, write to the Free Software
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			*/
33
34			#include "integrators/VelocityVerletIntegrator.hpp"
35
36			namespace oopse {
37
38
39			VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo* info){
40
41			}
42
43			VelocityVerletIntegrator::~VelocityVerletIntegrator() {
44
45			}
46
47			void VelocityVerletIntegrator::integrate() {
48
49			}
50
51			void VelocityVerletIntegrator::integrateStep() {
52
53			}
54
55			void VelocityVerletIntegrator::moveA() {
56
57			}
58
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;
94			}
95
96			temperature = info_->target_temp;
97
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
106
107			av2 = 2.0 * kebar / integrableObject->getMass();
108			vbar = sqrt( av2 );
109
110			// picks random velocities from a gaussian distribution
111			// centered on vbar
112
113			for (int k=0; k<3; k++) {
114			aVel[k] = vbar * gaussStream.getGaussian();
115			}
116
117			integrableObject->setVel( aVel );
118
119			if(integrableObject->isDirectional()){
120
121			I = integrableObject->getI();
122
123			if (integrableObject->isLinear()) {
124
125			l= integrableObject->linearAxis();
126			m = (l+1)%3;
127			n = (l+2)%3;
128
129			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();
134
135			} else {
136
137			for (int k = 0 ; k < 3; k++) {
138			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			}
165
166			}
167
168			void VelocityVerletIntegrator::calcForce(int needPotential, int needStress){
169
170			}
171
172			void VelocityVerletIntegrator::velocitize() {
173
174			}
175
176			void removeComDrift(){
177
178			}