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() {
    double runTime = info->run_time;

    double sampleTime = info->sampleTime;
    double statusTime = info->statusTime;
    double thermalTime = info->thermalTime;
    double resetTime = info->resetTime;

    double difference;
    double currSample;
    double currThermal;
    double currStatus;
    double currReset;

    int calcPot,
    calcStress;

    tStats = new Thermo(info);
    statOut = new StatWriter(info);
    dumpOut = new DumpWriter(info);

    atoms = info->atoms;

    fullStep_ = info->dt;
    halfStep_ = 0.5 * fullStep_;

    readyCheck();

    // remove center of mass drift velocity (in case we passed in a configuration
    // that was drifting
    tStats->removeCOMdrift();

    // initialize the forces before the first step

    calcForce(1, 1);

    //execute constraint algorithm to make sure at the very beginning the system is constrained  
    if (nConstrained) {
        constrainA();
        calcForce(1, 1);
        constrainB();
    }

    if (info->setTemp) {
        thermalize();
    }

    calcPot = 0;
    calcStress = 0;
    currSample = sampleTime + info->getTime();
    currThermal = thermalTime + info->getTime();
    currStatus = statusTime + info->getTime();
    currReset = resetTime + info->getTime();

    dumpOut->writeDump(info->getTime());
    statOut->writeStat(info->getTime());

#ifdef IS_MPI

    strcpy(checkPointMsg, "The integrator is ready to go.");
    MPIcheckPoint();

#endif // is_mpi

    while (info->getTime() < runTime && !stopIntegrator()) {
        difference = info->getTime() + fullStep_ - currStatus;

        if (difference > 0 || fabs(difference) < 1e - 4) {
            calcPot = 1;
            calcStress = 1;
        }

#ifdef PROFILE

        startProfile(pro1);

#endif

        integrateStep(calcPot, calcStress);

#ifdef PROFILE

        endProfile(pro1);

        startProfile(pro2);

#endif // profile

        info->incrTime(fullStep_);

        if (info->setTemp) {
            if (info->getTime() >= currThermal) {
                thermalize();
                currThermal += thermalTime;
            }
        }

        if (info->getTime() >= currSample) {
            dumpOut->writeDump(info->getTime());
            currSample += sampleTime;
        }

        if (info->getTime() >= currStatus) {
            statOut->writeStat(info->getTime());
            calcPot = 0;
            calcStress = 0;
            currStatus += statusTime;
        }

        if (info->resetIntegrator) {
            if (info->getTime() >= currReset) {
                this->resetIntegrator();
                currReset += resetTime;
            }
        }

#ifdef PROFILE

        endProfile(pro2);

#endif //profile

#ifdef IS_MPI

        strcpy(checkPointMsg, "successfully took a time step.");
        MPIcheckPoint();

#endif // is_mpi

    }

    dumpOut->writeFinal(info->getTime());

    delete dumpOut;
    delete statOut;
}

void VelocityVerletIntegrator::integrateStep() { }


void VelocityVerletIntegrator::thermalize() {

    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;
    }

}

void VelocityVerletIntegrator::velocitize(double temperature) {
    Vector3d aVel;
    Vector3d aJ;
    Mat3x3d I;
    int l;
    int m;
    int 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;

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

    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();

    removeComDrift(vdrift);

}

void VelocityVerletIntegrator::calcForce(bool needPotential,
                                         bool needStress) { }

void VelocityVerletIntegrator::removeComDrift(const Vector3d& vdrift) {

    std::vector<Molecule *>::iterator i;
    std::vector<StuntDouble *>::iterator j;
    Molecule * mol;
    StuntDouble * integrableObject;
    
    //  Corrects for the center of mass drift.
    // sums all the momentum and divides by total mass.
    for( mol = info_->beginMolecule(i); mol != NULL;
        mol = info_->nextMolecule(i) ) {
        for( integrableObject = mol->beginIntegrableObject(j);
            integrableObject != NULL;
            integrableObject = mol->nextIntegrableObject(j) ) {
            integrableObject->setVel(integrableObject->getVel() - vdrift);
        }
    }

}

} //end namespace oopse
Revision:	1725
Committed:	Wed Nov 10 22:01:06 2004 UTC (19 years, 8 months ago) by tim
File size:	7969 byte(s)
Log Message:	another painful day (1) SimCreator, SimInfo, mpiSimulation (2) DumpReader, DumpWriter (InitializeFrom File will be removed) (3) ForceField (at least LJ) and BondType, BendType, TorsionType (4)Integrator (5)oopse.cpp (6)visitors & Dump2XYZ (7)SimpleBuilder (8)Constraint & ZConstraint
#	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	tim	1725	/**
27			* @file VelocityVerletIntegrator.cpp
28			* @author tlin
29			* @date 11/09/2004
30			* @time 16:16am
31			* @version 1.0
32			*/
33
34	tim	1722	#include "integrators/VelocityVerletIntegrator.hpp"
35
36			namespace oopse {
37	tim	1725	VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) { }
38	tim	1722
39	tim	1725	VelocityVerletIntegrator::~VelocityVerletIntegrator() { }
40	tim	1722
41	tim	1725	void VelocityVerletIntegrator::integrate() {
42			double runTime = info->run_time;
43	tim	1722
44	tim	1725	double sampleTime = info->sampleTime;
45			double statusTime = info->statusTime;
46			double thermalTime = info->thermalTime;
47			double resetTime = info->resetTime;
48	tim	1722
49	tim	1725	double difference;
50			double currSample;
51			double currThermal;
52			double currStatus;
53			double currReset;
54	tim	1722
55	tim	1725	int calcPot,
56			calcStress;
57	tim	1722
58	tim	1725	tStats = new Thermo(info);
59			statOut = new StatWriter(info);
60			dumpOut = new DumpWriter(info);
61	tim	1722
62	tim	1725	atoms = info->atoms;
63	tim	1722
64	tim	1725	fullStep_ = info->dt;
65			halfStep_ = 0.5 * fullStep_;
66	tim	1722
67	tim	1725	readyCheck();
68	tim	1722
69	tim	1725	// remove center of mass drift velocity (in case we passed in a configuration
70			// that was drifting
71			tStats->removeCOMdrift();
72	tim	1722
73	tim	1725	// initialize the forces before the first step
74	tim	1722
75	tim	1725	calcForce(1, 1);
76	tim	1722
77	tim	1725	//execute constraint algorithm to make sure at the very beginning the system is constrained
78			if (nConstrained) {
79			constrainA();
80			calcForce(1, 1);
81			constrainB();
82			}
83
84			if (info->setTemp) {
85			thermalize();
86			}
87
88			calcPot = 0;
89			calcStress = 0;
90			currSample = sampleTime + info->getTime();
91			currThermal = thermalTime + info->getTime();
92			currStatus = statusTime + info->getTime();
93			currReset = resetTime + info->getTime();
94
95			dumpOut->writeDump(info->getTime());
96			statOut->writeStat(info->getTime());
97
98			#ifdef IS_MPI
99
100			strcpy(checkPointMsg, "The integrator is ready to go.");
101			MPIcheckPoint();
102
103			#endif // is_mpi
104
105			while (info->getTime() < runTime && !stopIntegrator()) {
106			difference = info->getTime() + fullStep_ - currStatus;
107
108			if (difference > 0 \|\| fabs(difference) < 1e - 4) {
109			calcPot = 1;
110			calcStress = 1;
111			}
112
113			#ifdef PROFILE
114
115			startProfile(pro1);
116
117			#endif
118
119			integrateStep(calcPot, calcStress);
120
121			#ifdef PROFILE
122
123			endProfile(pro1);
124
125			startProfile(pro2);
126
127			#endif // profile
128
129			info->incrTime(fullStep_);
130
131			if (info->setTemp) {
132			if (info->getTime() >= currThermal) {
133			thermalize();
134			currThermal += thermalTime;
135			}
136			}
137
138			if (info->getTime() >= currSample) {
139			dumpOut->writeDump(info->getTime());
140			currSample += sampleTime;
141			}
142
143			if (info->getTime() >= currStatus) {
144			statOut->writeStat(info->getTime());
145			calcPot = 0;
146			calcStress = 0;
147			currStatus += statusTime;
148			}
149
150			if (info->resetIntegrator) {
151			if (info->getTime() >= currReset) {
152			this->resetIntegrator();
153			currReset += resetTime;
154			}
155			}
156
157			#ifdef PROFILE
158
159			endProfile(pro2);
160
161			#endif //profile
162
163			#ifdef IS_MPI
164
165			strcpy(checkPointMsg, "successfully took a time step.");
166			MPIcheckPoint();
167
168			#endif // is_mpi
169
170			}
171
172			dumpOut->writeFinal(info->getTime());
173
174			delete dumpOut;
175			delete statOut;
176	tim	1722	}
177
178	tim	1725	void VelocityVerletIntegrator::integrateStep() { }
179
180
181	tim	1722	void VelocityVerletIntegrator::thermalize() {
182
183	tim	1725	if (!info_->have_target_temp) {
184			sprintf(painCave.errMsg,
185			"You can't resample the velocities without a targetTemp!\n");
186			painCave.isFatal = 1;
187			painCave.severity = OOPSE_ERROR;
188			simError();
189			return;
190			}
191
192	tim	1722	}
193
194	tim	1725	void VelocityVerletIntegrator::velocitize(double temperature) {
195	tim	1722	Vector3d aVel;
196			Vector3d aJ;
197			Mat3x3d I;
198	tim	1725	int l;
199			int m;
200			int n; // velocity randomizer loop counts
201	tim	1722	Vector3d vdrift;
202			double vbar;
203	tim	1725	const double kb = 8.31451e - 7; // kb in amu, angstroms, fs, etc.
204	tim	1722	double av2;
205			double kebar;
206
207	tim	1725	std::vector<Molecule *>::iterator i;
208			std::vector<StuntDouble *>::iterator j;
209			Molecule * mol;
210			StuntDouble * integrableObject;
211	tim	1722	gaussianSPRNG gaussStream(info_->getSeed());
212
213	tim	1725	kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf());
214	tim	1722
215	tim	1725	for( mol = info_->beginMolecule(i); mol != NULL;
216			mol = info_->nextMolecule(i) ) {
217			for( integrableObject = mol->beginIntegrableObject(j);
218			integrableObject != NULL;
219			integrableObject = mol->nextIntegrableObject(j) ) {
220	tim	1722
221			// uses equipartition theory to solve for vbar in angstrom/fs
222
223			av2 = 2.0 * kebar / integrableObject->getMass();
224	tim	1725	vbar = sqrt(av2);
225	tim	1722
226			// picks random velocities from a gaussian distribution
227			// centered on vbar
228
229	tim	1725	for( int k = 0; k < 3; k++ ) {
230	tim	1722	aVel[k] = vbar * gaussStream.getGaussian();
231			}
232
233	tim	1725	integrableObject->setVel(aVel);
234	tim	1722
235	tim	1725	if (integrableObject->isDirectional()) {
236	tim	1722	I = integrableObject->getI();
237
238			if (integrableObject->isLinear()) {
239	tim	1725	l = integrableObject->linearAxis();
240			m = (l + 1) % 3;
241			n = (l + 2) % 3;
242	tim	1722
243			aJ[l] = 0.0;
244	tim	1725	vbar = sqrt(2.0 * kebar * I(m, m));
245	tim	1722	aJ[m] = vbar * gaussStream.getGaussian();
246	tim	1725	vbar = sqrt(2.0 * kebar * I(n, n));
247	tim	1722	aJ[n] = vbar * gaussStream.getGaussian();
248			} else {
249	tim	1725	for( int k = 0; k < 3; k++ ) {
250			vbar = sqrt(2.0 * kebar * I(k, k));
251	tim	1722	aJ[k] = vbar * gaussStream.getGaussian();
252	tim	1725	}
253	tim	1722	} // else isLinear
254
255	tim	1725	integrableObject->setJ(aJ);
256			} //isDirectional
257			}
258			} //end for (mol = beginMolecule(i); ...)
259	tim	1722
260			// Get the Center of Mass drift velocity.
261			vdrift = info_->getComVel();
262
263	tim	1725	removeComDrift(vdrift);
264
265	tim	1722	}
266
267	tim	1725	void VelocityVerletIntegrator::calcForce(bool needPotential,
268			bool needStress) { }
269	tim	1722
270	tim	1725	void VelocityVerletIntegrator::removeComDrift(const Vector3d& vdrift) {
271	tim	1722
272	tim	1725	std::vector<Molecule *>::iterator i;
273			std::vector<StuntDouble *>::iterator j;
274			Molecule * mol;
275			StuntDouble * integrableObject;
276
277			// Corrects for the center of mass drift.
278			// sums all the momentum and divides by total mass.
279			for( mol = info_->beginMolecule(i); mol != NULL;
280			mol = info_->nextMolecule(i) ) {
281			for( integrableObject = mol->beginIntegrableObject(j);
282			integrableObject != NULL;
283			integrableObject = mol->nextIntegrableObject(j) ) {
284			integrableObject->setVel(integrableObject->getVel() - vdrift);
285			}
286			}
287	tim	1722
288			}
289
290	tim	1725	} //end namespace oopse