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, 7 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
#	Content
1	/*
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	VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) { }
38
39	VelocityVerletIntegrator::~VelocityVerletIntegrator() { }
40
41	void VelocityVerletIntegrator::integrate() {
42	double runTime = info->run_time;
43
44	double sampleTime = info->sampleTime;
45	double statusTime = info->statusTime;
46	double thermalTime = info->thermalTime;
47	double resetTime = info->resetTime;
48
49	double difference;
50	double currSample;
51	double currThermal;
52	double currStatus;
53	double currReset;
54
55	int calcPot,
56	calcStress;
57
58	tStats = new Thermo(info);
59	statOut = new StatWriter(info);
60	dumpOut = new DumpWriter(info);
61
62	atoms = info->atoms;
63
64	fullStep_ = info->dt;
65	halfStep_ = 0.5 * fullStep_;
66
67	readyCheck();
68
69	// remove center of mass drift velocity (in case we passed in a configuration
70	// that was drifting
71	tStats->removeCOMdrift();
72
73	// initialize the forces before the first step
74
75	calcForce(1, 1);
76
77	//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	}
177
178	void VelocityVerletIntegrator::integrateStep() { }
179
180
181	void VelocityVerletIntegrator::thermalize() {
182
183	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	}
193
194	void VelocityVerletIntegrator::velocitize(double temperature) {
195	Vector3d aVel;
196	Vector3d aJ;
197	Mat3x3d I;
198	int l;
199	int m;
200	int n; // velocity randomizer loop counts
201	Vector3d vdrift;
202	double vbar;
203	const double kb = 8.31451e - 7; // kb in amu, angstroms, fs, etc.
204	double av2;
205	double kebar;
206
207	std::vector<Molecule *>::iterator i;
208	std::vector<StuntDouble *>::iterator j;
209	Molecule * mol;
210	StuntDouble * integrableObject;
211	gaussianSPRNG gaussStream(info_->getSeed());
212
213	kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf());
214
215	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
221	// uses equipartition theory to solve for vbar in angstrom/fs
222
223	av2 = 2.0 * kebar / integrableObject->getMass();
224	vbar = sqrt(av2);
225
226	// picks random velocities from a gaussian distribution
227	// centered on vbar
228
229	for( int k = 0; k < 3; k++ ) {
230	aVel[k] = vbar * gaussStream.getGaussian();
231	}
232
233	integrableObject->setVel(aVel);
234
235	if (integrableObject->isDirectional()) {
236	I = integrableObject->getI();
237
238	if (integrableObject->isLinear()) {
239	l = integrableObject->linearAxis();
240	m = (l + 1) % 3;
241	n = (l + 2) % 3;
242
243	aJ[l] = 0.0;
244	vbar = sqrt(2.0 * kebar * I(m, m));
245	aJ[m] = vbar * gaussStream.getGaussian();
246	vbar = sqrt(2.0 * kebar * I(n, n));
247	aJ[n] = vbar * gaussStream.getGaussian();
248	} else {
249	for( int k = 0; k < 3; k++ ) {
250	vbar = sqrt(2.0 * kebar * I(k, k));
251	aJ[k] = vbar * gaussStream.getGaussian();
252	}
253	} // else isLinear
254
255	integrableObject->setJ(aJ);
256	} //isDirectional
257	}
258	} //end for (mol = beginMolecule(i); ...)
259
260	// Get the Center of Mass drift velocity.
261	vdrift = info_->getComVel();
262
263	removeComDrift(vdrift);
264
265	}
266
267	void VelocityVerletIntegrator::calcForce(bool needPotential,
268	bool needStress) { }
269
270	void VelocityVerletIntegrator::removeComDrift(const Vector3d& vdrift) {
271
272	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
288	}
289
290	} //end namespace oopse