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) : Integrator(info) { }

VelocityVerletIntegrator::~VelocityVerletIntegrator() { }

void VelocityVerletIntegrator::integrate() {
    double runTime = info_->run_time;

    int calcPot;
    int calcStress;


    dt_ = info_->dt;
    dt2_ = 0.5 * dt_;

    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;

    tStats = new Thermo(info_);
    statOut = new StatWriter(info_);
    dumpOut = new DumpWriter(info_);
    
    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) {
        difference = info_->getTime() + dt_ - currStatus;

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

        //notify before integrateStep
        notify()
        integrateStep(calcPot, calcStress);
        
        info_->incrTime(dt_);
        
        //notify after integratreStep        
        notify();
        
        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 IS_MPI

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

#endif // is_mpi

    }

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

    delete dumpOut;
    delete statOut;
}

void VelocityVerletIntegrator::integrateStep() { 

    moveA();
    calcForce(bool needPotential, bool needStress);
    moveB();
}


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::calcForce(bool needPotential,
                                         bool needStress) { }


} //end namespace oopse
Revision:	1774
Committed:	Tue Nov 23 23:12:23 2004 UTC (21 years, 2 months ago) by tim
File size:	4455 byte(s)
Log Message:	refactory NPT integrator
#	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	1756	VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info) { }
38	tim	1722
39	tim	1725	VelocityVerletIntegrator::~VelocityVerletIntegrator() { }
40	tim	1722
41	tim	1725	void VelocityVerletIntegrator::integrate() {
42	tim	1756	double runTime = info_->run_time;
43	tim	1722
44	tim	1756	int calcPot;
45			int calcStress;
46	tim	1722
47
48	tim	1774	dt_ = info_->dt;
49			dt2_ = 0.5 * dt_;
50	tim	1722
51	tim	1725	readyCheck();
52	tim	1722
53	tim	1725	// remove center of mass drift velocity (in case we passed in a configuration
54			// that was drifting
55			tStats->removeCOMdrift();
56	tim	1722
57	tim	1725	// initialize the forces before the first step
58	tim	1722
59	tim	1725	calcForce(1, 1);
60	tim	1722
61	tim	1725	//execute constraint algorithm to make sure at the very beginning the system is constrained
62			if (nConstrained) {
63			constrainA();
64			calcForce(1, 1);
65			constrainB();
66			}
67
68	tim	1756	if (info_->setTemp) {
69	tim	1725	thermalize();
70			}
71
72			calcPot = 0;
73			calcStress = 0;
74
75	tim	1756	tStats = new Thermo(info_);
76			statOut = new StatWriter(info_);
77			dumpOut = new DumpWriter(info_);
78
79			dumpOut->writeDump(info_->getTime());
80			statOut->writeStat(info_->getTime());
81	tim	1725
82			#ifdef IS_MPI
83
84			strcpy(checkPointMsg, "The integrator is ready to go.");
85			MPIcheckPoint();
86
87			#endif // is_mpi
88
89	tim	1756	while (info_->getTime() < runTime) {
90	tim	1774	difference = info_->getTime() + dt_ - currStatus;
91	tim	1725
92			if (difference > 0 \|\| fabs(difference) < 1e - 4) {
93			calcPot = 1;
94			calcStress = 1;
95			}
96
97	tim	1756	//notify before integrateStep
98			notify()
99	tim	1725	integrateStep(calcPot, calcStress);
100	tim	1756
101	tim	1774	info_->incrTime(dt_);
102	tim	1756
103			//notify after integratreStep
104			notify();
105
106			if (info_->setTemp) {
107			if (info_->getTime() >= currThermal) {
108	tim	1725	thermalize();
109			currThermal += thermalTime;
110			}
111			}
112
113	tim	1756	if (info_->getTime() >= currSample) {
114			dumpOut->writeDump(info_->getTime());
115	tim	1725	currSample += sampleTime;
116			}
117
118	tim	1756	if (info_->getTime() >= currStatus) {
119			statOut->writeStat(info_->getTime());
120	tim	1725	calcPot = 0;
121			calcStress = 0;
122			currStatus += statusTime;
123			}
124
125	tim	1756	if (info_->resetIntegrator) {
126			if (info_->getTime() >= currReset) {
127	tim	1725	this->resetIntegrator();
128			currReset += resetTime;
129			}
130			}
131
132			#ifdef IS_MPI
133
134			strcpy(checkPointMsg, "successfully took a time step.");
135			MPIcheckPoint();
136
137			#endif // is_mpi
138
139			}
140
141	tim	1756	dumpOut->writeFinal(info_->getTime());
142	tim	1725
143			delete dumpOut;
144			delete statOut;
145	tim	1722	}
146
147	tim	1774	void VelocityVerletIntegrator::integrateStep() {
148	tim	1725
149	tim	1774	moveA();
150			calcForce(bool needPotential, bool needStress);
151			moveB();
152			}
153	tim	1725
154	tim	1774
155	tim	1722	void VelocityVerletIntegrator::thermalize() {
156
157	tim	1756	if (!info__->have_target_temp) {
158	tim	1725	sprintf(painCave.errMsg,
159			"You can't resample the velocities without a targetTemp!\n");
160			painCave.isFatal = 1;
161			painCave.severity = OOPSE_ERROR;
162			simError();
163			return;
164			}
165
166	tim	1722	}
167
168	tim	1725	void VelocityVerletIntegrator::calcForce(bool needPotential,
169			bool needStress) { }
170	tim	1722
171
172	tim	1725	} //end namespace oopse