src/integrators/VelocityVerletIntegrator.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
/**
 * @file VelocityVerletIntegrator.cpp
 * @author tlin
 * @date 11/09/2004
 * @time 16:16am
 * @version 1.0
 */

#include "integrators/VelocityVerletIntegrator.hpp"
#include "integrators/DLM.hpp"
#include "utils/StringUtils.hpp"

namespace oopse {
  VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info), rotAlgo(NULL) { 
    dt2 = 0.5 * dt;
    rotAlgo = new DLM();
    rattle = new Rattle(info);
  }
  
  VelocityVerletIntegrator::~VelocityVerletIntegrator() { 
    delete rotAlgo;
    delete rattle;
  }
  
  void VelocityVerletIntegrator::initialize(){
    
    forceMan_->init();
    
    // remove center of mass drift velocity (in case we passed in a
    // configuration that was drifting)
    velocitizer_->removeComDrift();
    
    // initialize the forces before the first step
    calcForce(true, true);
    
    // execute the constraint algorithm to make sure that the system is
    // constrained at the very beginning  
    if (info_->getNGlobalConstraints() > 0) {
      rattle->constraintA();
      calcForce(true, true);
      rattle->constraintB();      
      //copy the current snapshot to previous snapshot
      info_->getSnapshotManager()->advance();
    }
    
    if (needVelocityScaling) {
      velocitizer_->velocitize(targetScalingTemp);
    }
    
    dumpWriter = createDumpWriter();
    
    statWriter = createStatWriter();
 
    dumpWriter->writeDumpAndEor();

    //save statistics, before writeStat,  we must save statistics
    thermo.saveStat();
    saveConservedQuantity();
    if (simParams->getUseRNEMD())
      rnemd_->getStatus();

    statWriter->writeStat(currentSnapshot_->statData);
    
    currSample = sampleTime + currentSnapshot_->getTime();
    currStatus =  statusTime + currentSnapshot_->getTime();
    currThermal = thermalTime + currentSnapshot_->getTime();
    if (needReset) {
      currReset = resetTime + currentSnapshot_->getTime();
    }
    if (simParams->getUseRNEMD()){
      currRNEMD = RNEMD_swapTime + currentSnapshot_->getTime();
    }
    needPotential = false;
    needStress = false;       
    
  }

  void VelocityVerletIntegrator::doIntegrate() {
  
  
    initialize();
  
    while (currentSnapshot_->getTime() < runTime) {
    
      preStep();
    
      integrateStep();
    
      postStep();
    
    }
  
    finalize();
  
  }


  void VelocityVerletIntegrator::preStep() {
    RealType difference = currentSnapshot_->getTime() + dt - currStatus;
  
    if (difference > 0 || fabs(difference) < oopse::epsilon) {
      needPotential = true;
      needStress = true;   
    }
  }

  void VelocityVerletIntegrator::postStep() {

    //save snapshot
    info_->getSnapshotManager()->advance();
  
    //increase time
    currentSnapshot_->increaseTime(dt);        
   
    if (needVelocityScaling) {
      if (currentSnapshot_->getTime() >= currThermal) {
        velocitizer_->velocitize(targetScalingTemp);
        currThermal += thermalTime;
      }
    }
    if (useRNEMD) {
      if (currentSnapshot_->getTime() >= currRNEMD) {
        rnemd_->doSwap();
        currRNEMD += RNEMD_swapTime;
      }
    }
    
    if (currentSnapshot_->getTime() >= currSample) {
      dumpWriter->writeDumpAndEor();
      
      currSample += sampleTime;
    }
    
    if (currentSnapshot_->getTime() >= currStatus) {
      //save statistics, before writeStat,  we must save statistics
      thermo.saveStat();
      saveConservedQuantity();

      if (simParams->getUseRNEMD())
        rnemd_->getStatus();

      statWriter->writeStat(currentSnapshot_->statData);
      
      needPotential = false;
      needStress = false;
      currStatus += statusTime;
    }
    
    if (needReset && currentSnapshot_->getTime() >= currReset) {    
      resetIntegrator();
      currReset += resetTime;
    }        
  }


  void VelocityVerletIntegrator::finalize() {
    dumpWriter->writeEor();
  
    delete dumpWriter;
    delete statWriter;
  
    dumpWriter = NULL;
    statWriter = NULL;
  
  }

  void VelocityVerletIntegrator::integrateStep() { 
  
    moveA();
    calcForce(needPotential, needStress);
    moveB();
  }


  void VelocityVerletIntegrator::calcForce(bool needPotential,
                                           bool needStress) { 
    forceMan_->calcForces(needPotential, needStress);
  }

  DumpWriter* VelocityVerletIntegrator::createDumpWriter() {
    return new DumpWriter(info_);
  }

  StatWriter* VelocityVerletIntegrator::createStatWriter() {

    std::string statFileFormatString = simParams->getStatFileFormat();
    StatsBitSet mask = parseStatFileFormat(statFileFormatString);
   
    // if we're doing a thermodynamic integration, we'll want the raw
    // potential as well as the full potential:

 
    if (simParams->getUseThermodynamicIntegration()) 
      mask.set(Stats::VRAW);

    // if we've got restraints turned on, we'll also want a report of the
    // total harmonic restraints
    if (simParams->getUseRestraints()){
      mask.set(Stats::VHARM);
    }

    if (simParams->havePrintPressureTensor() && 
        simParams->getPrintPressureTensor()){
      mask.set(Stats::PRESSURE_TENSOR_XX);
      mask.set(Stats::PRESSURE_TENSOR_XY);
      mask.set(Stats::PRESSURE_TENSOR_XZ);
      mask.set(Stats::PRESSURE_TENSOR_YX);
      mask.set(Stats::PRESSURE_TENSOR_YY);
      mask.set(Stats::PRESSURE_TENSOR_YZ);
      mask.set(Stats::PRESSURE_TENSOR_ZX);
      mask.set(Stats::PRESSURE_TENSOR_ZY);
      mask.set(Stats::PRESSURE_TENSOR_ZZ);
    }
    
    if (simParams->getAccumulateBoxDipole()) {
      mask.set(Stats::BOX_DIPOLE_X);
      mask.set(Stats::BOX_DIPOLE_Y);
      mask.set(Stats::BOX_DIPOLE_Z);
    }
   
    if (simParams->haveTaggedAtomPair() && simParams->havePrintTaggedPairDistance()) {
      if (simParams->getPrintTaggedPairDistance()) {
        mask.set(Stats::TAGGED_PAIR_DISTANCE);
      }
    }

    if (simParams->getUseRNEMD()) {
      mask.set(Stats::RNEMD_SWAP_TOTAL);
    }
      

     return new StatWriter(info_->getStatFileName(), mask);
  }


} //end namespace oopse
Revision:	1360
Committed:	Mon Sep 7 16:31:51 2009 UTC (15 years, 7 months ago) by cli2
Original Path:	*trunk/src/integrators/VelocityVerletIntegrator.cpp*
File size:	8053 byte(s)
Log Message:	Added new restraint infrastructure Added MolecularRestraints Added ObjectRestraints Added RestraintStamp Updated thermodynamic integration to use ObjectRestraints Added Quaternion mathematics for twist swing decompositions Significantly updated RestWriter and RestReader to use dump-like files Added selections for x, y, and z coordinates of atoms Removed monolithic Restraints class Fixed a few bugs in gradients of Euler angles in DirectionalAtom and RigidBody Added some rotational capabilities to prinicpalAxisCalculator
#	User	Rev	Content
1	gezelter	507	/*
2	gezelter	246	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
41
42			/**
43			* @file VelocityVerletIntegrator.cpp
44			* @author tlin
45			* @date 11/09/2004
46			* @time 16:16am
47			* @version 1.0
48			*/
49
50			#include "integrators/VelocityVerletIntegrator.hpp"
51			#include "integrators/DLM.hpp"
52	tim	276	#include "utils/StringUtils.hpp"
53	gezelter	246
54			namespace oopse {
55	chrisfen	417	VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info), rotAlgo(NULL) {
56	gezelter	246	dt2 = 0.5 * dt;
57			rotAlgo = new DLM();
58			rattle = new Rattle(info);
59	chrisfen	417	}
60
61			VelocityVerletIntegrator::~VelocityVerletIntegrator() {
62	gezelter	246	delete rotAlgo;
63			delete rattle;
64	chrisfen	417	}
65
66			void VelocityVerletIntegrator::initialize(){
67
68	gezelter	246	forceMan_->init();
69	chrisfen	417
70	chrisfen	993	// remove center of mass drift velocity (in case we passed in a
71			// configuration that was drifting)
72	gezelter	246	velocitizer_->removeComDrift();
73	chrisfen	417
74	gezelter	246	// initialize the forces before the first step
75			calcForce(true, true);
76	chrisfen	417
77	chrisfen	993	// execute the constraint algorithm to make sure that the system is
78			// constrained at the very beginning
79	gezelter	246	if (info_->getNGlobalConstraints() > 0) {
80	chrisfen	417	rattle->constraintA();
81			calcForce(true, true);
82	chrisfen	993	rattle->constraintB();
83			//copy the current snapshot to previous snapshot
84			info_->getSnapshotManager()->advance();
85	gezelter	246	}
86	chrisfen	417
87	gezelter	246	if (needVelocityScaling) {
88	chrisfen	417	velocitizer_->velocitize(targetScalingTemp);
89	gezelter	246	}
90
91			dumpWriter = createDumpWriter();
92	chrisfen	417
93	gezelter	246	statWriter = createStatWriter();
94	chrisfen	993
95			dumpWriter->writeDumpAndEor();
96
97	gezelter	246	//save statistics, before writeStat, we must save statistics
98			thermo.saveStat();
99			saveConservedQuantity();
100	skuang	1341	if (simParams->getUseRNEMD())
101			rnemd_->getStatus();
102
103	gezelter	246	statWriter->writeStat(currentSnapshot_->statData);
104	chrisfen	417
105	gezelter	246	currSample = sampleTime + currentSnapshot_->getTime();
106	gezelter	1329	currStatus = statusTime + currentSnapshot_->getTime();
107	tim	546	currThermal = thermalTime + currentSnapshot_->getTime();
108			if (needReset) {
109			currReset = resetTime + currentSnapshot_->getTime();
110			}
111	gezelter	1329	if (simParams->getUseRNEMD()){
112			currRNEMD = RNEMD_swapTime + currentSnapshot_->getTime();
113			}
114	gezelter	246	needPotential = false;
115			needStress = false;
116	chrisfen	417
117	gezelter	507	}
118	chrisfen	417
119	gezelter	507	void VelocityVerletIntegrator::doIntegrate() {
120	chrisfen	417
121
122	gezelter	507	initialize();
123	chrisfen	417
124	gezelter	507	while (currentSnapshot_->getTime() < runTime) {
125	gezelter	246
126	gezelter	507	preStep();
127	chrisfen	417
128	gezelter	507	integrateStep();
129	chrisfen	417
130	gezelter	507	postStep();
131	chrisfen	417
132	gezelter	507	}
133	chrisfen	417
134	gezelter	507	finalize();
135	chrisfen	417
136	gezelter	507	}
137	gezelter	246
138
139	gezelter	507	void VelocityVerletIntegrator::preStep() {
140	tim	963	RealType difference = currentSnapshot_->getTime() + dt - currStatus;
141	chrisfen	417
142	gezelter	507	if (difference > 0 \|\| fabs(difference) < oopse::epsilon) {
143			needPotential = true;
144			needStress = true;
145			}
146	chrisfen	417	}
147	gezelter	246
148	gezelter	507	void VelocityVerletIntegrator::postStep() {
149	tim	749
150	gezelter	507	//save snapshot
151			info_->getSnapshotManager()->advance();
152	chrisfen	417
153	gezelter	507	//increase time
154			currentSnapshot_->increaseTime(dt);
155	chrisfen	993
156	gezelter	507	if (needVelocityScaling) {
157			if (currentSnapshot_->getTime() >= currThermal) {
158			velocitizer_->velocitize(targetScalingTemp);
159			currThermal += thermalTime;
160			}
161	chrisfen	417	}
162	gezelter	1329	if (useRNEMD) {
163			if (currentSnapshot_->getTime() >= currRNEMD) {
164			rnemd_->doSwap();
165			currRNEMD += RNEMD_swapTime;
166			}
167			}
168
169	gezelter	507	if (currentSnapshot_->getTime() >= currSample) {
170			dumpWriter->writeDumpAndEor();
171	gezelter	1329
172	gezelter	507	currSample += sampleTime;
173			}
174	gezelter	1329
175	gezelter	507	if (currentSnapshot_->getTime() >= currStatus) {
176			//save statistics, before writeStat, we must save statistics
177			thermo.saveStat();
178			saveConservedQuantity();
179	skuang	1341
180	skuang	1338	if (simParams->getUseRNEMD())
181			rnemd_->getStatus();
182	skuang	1341
183			statWriter->writeStat(currentSnapshot_->statData);
184	gezelter	1329
185	gezelter	507	needPotential = false;
186			needStress = false;
187			currStatus += statusTime;
188			}
189	gezelter	1329
190			if (needReset && currentSnapshot_->getTime() >= currReset) {
191			resetIntegrator();
192			currReset += resetTime;
193			}
194	gezelter	507	}
195	gezelter	246
196
197	gezelter	507	void VelocityVerletIntegrator::finalize() {
198			dumpWriter->writeEor();
199	chrisfen	417
200	gezelter	507	delete dumpWriter;
201			delete statWriter;
202	chrisfen	417
203	gezelter	507	dumpWriter = NULL;
204			statWriter = NULL;
205	chrisfen	417
206	gezelter	507	}
207	gezelter	246
208	gezelter	507	void VelocityVerletIntegrator::integrateStep() {
209	chrisfen	417
210	gezelter	507	moveA();
211			calcForce(needPotential, needStress);
212			moveB();
213			}
214	gezelter	246
215
216	gezelter	507	void VelocityVerletIntegrator::calcForce(bool needPotential,
217			bool needStress) {
218			forceMan_->calcForces(needPotential, needStress);
219			}
220	gezelter	246
221	gezelter	507	DumpWriter* VelocityVerletIntegrator::createDumpWriter() {
222			return new DumpWriter(info_);
223			}
224	gezelter	246
225	gezelter	507	StatWriter* VelocityVerletIntegrator::createStatWriter() {
226	tim	681
227			std::string statFileFormatString = simParams->getStatFileFormat();
228			StatsBitSet mask = parseStatFileFormat(statFileFormatString);
229	cli2	1360
230			// if we're doing a thermodynamic integration, we'll want the raw
231			// potential as well as the full potential:
232
233
234			if (simParams->getUseThermodynamicIntegration())
235	gezelter	507	mask.set(Stats::VRAW);
236	cli2	1360
237			// if we've got restraints turned on, we'll also want a report of the
238			// total harmonic restraints
239			if (simParams->getUseRestraints()){
240	gezelter	507	mask.set(Stats::VHARM);
241			}
242	tim	541
243	chrisfen	993	if (simParams->havePrintPressureTensor() &&
244			simParams->getPrintPressureTensor()){
245	gezelter	1291	mask.set(Stats::PRESSURE_TENSOR_XX);
246			mask.set(Stats::PRESSURE_TENSOR_XY);
247			mask.set(Stats::PRESSURE_TENSOR_XZ);
248			mask.set(Stats::PRESSURE_TENSOR_YX);
249			mask.set(Stats::PRESSURE_TENSOR_YY);
250			mask.set(Stats::PRESSURE_TENSOR_YZ);
251			mask.set(Stats::PRESSURE_TENSOR_ZX);
252			mask.set(Stats::PRESSURE_TENSOR_ZY);
253			mask.set(Stats::PRESSURE_TENSOR_ZZ);
254	tim	541	}
255
256	chrisfen	998	if (simParams->getAccumulateBoxDipole()) {
257			mask.set(Stats::BOX_DIPOLE_X);
258			mask.set(Stats::BOX_DIPOLE_Y);
259			mask.set(Stats::BOX_DIPOLE_Z);
260			}
261	gezelter	1329
262	gezelter	1291	if (simParams->haveTaggedAtomPair() && simParams->havePrintTaggedPairDistance()) {
263			if (simParams->getPrintTaggedPairDistance()) {
264			mask.set(Stats::TAGGED_PAIR_DISTANCE);
265			}
266			}
267	gezelter	1329
268			if (simParams->getUseRNEMD()) {
269			mask.set(Stats::RNEMD_SWAP_TOTAL);
270			}
271	gezelter	1291
272	chrisfen	998
273	tim	681	return new StatWriter(info_->getStatFileName(), mask);
274	chrisfen	417	}
275	gezelter	246
276
277			} //end namespace oopse