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:	3520
Committed:	Mon Sep 7 16:31:51 2009 UTC (15 years, 10 months ago) by cli2
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	2204	/*
2	gezelter	1930	* 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	1960	#include "utils/StringUtils.hpp"
53	gezelter	1930
54			namespace oopse {
55	chrisfen	2101	VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info), rotAlgo(NULL) {
56	gezelter	1930	dt2 = 0.5 * dt;
57			rotAlgo = new DLM();
58			rattle = new Rattle(info);
59	chrisfen	2101	}
60
61			VelocityVerletIntegrator::~VelocityVerletIntegrator() {
62	gezelter	1930	delete rotAlgo;
63			delete rattle;
64	chrisfen	2101	}
65
66			void VelocityVerletIntegrator::initialize(){
67
68	gezelter	1930	forceMan_->init();
69	chrisfen	2101
70	chrisfen	2879	// remove center of mass drift velocity (in case we passed in a
71			// configuration that was drifting)
72	gezelter	1930	velocitizer_->removeComDrift();
73	chrisfen	2101
74	gezelter	1930	// initialize the forces before the first step
75			calcForce(true, true);
76	chrisfen	2101
77	chrisfen	2879	// execute the constraint algorithm to make sure that the system is
78			// constrained at the very beginning
79	gezelter	1930	if (info_->getNGlobalConstraints() > 0) {
80	chrisfen	2101	rattle->constraintA();
81			calcForce(true, true);
82	chrisfen	2879	rattle->constraintB();
83			//copy the current snapshot to previous snapshot
84			info_->getSnapshotManager()->advance();
85	gezelter	1930	}
86	chrisfen	2101
87	gezelter	1930	if (needVelocityScaling) {
88	chrisfen	2101	velocitizer_->velocitize(targetScalingTemp);
89	gezelter	1930	}
90
91			dumpWriter = createDumpWriter();
92	chrisfen	2101
93	gezelter	1930	statWriter = createStatWriter();
94	chrisfen	2879
95			dumpWriter->writeDumpAndEor();
96
97	gezelter	1930	//save statistics, before writeStat, we must save statistics
98			thermo.saveStat();
99			saveConservedQuantity();
100	skuang	3501	if (simParams->getUseRNEMD())
101			rnemd_->getStatus();
102
103	gezelter	1930	statWriter->writeStat(currentSnapshot_->statData);
104	chrisfen	2101
105	gezelter	1930	currSample = sampleTime + currentSnapshot_->getTime();
106	gezelter	3488	currStatus = statusTime + currentSnapshot_->getTime();
107	tim	2243	currThermal = thermalTime + currentSnapshot_->getTime();
108			if (needReset) {
109			currReset = resetTime + currentSnapshot_->getTime();
110			}
111	gezelter	3488	if (simParams->getUseRNEMD()){
112			currRNEMD = RNEMD_swapTime + currentSnapshot_->getTime();
113			}
114	gezelter	1930	needPotential = false;
115			needStress = false;
116	chrisfen	2101
117	gezelter	2204	}
118	chrisfen	2101
119	gezelter	2204	void VelocityVerletIntegrator::doIntegrate() {
120	chrisfen	2101
121
122	gezelter	2204	initialize();
123	chrisfen	2101
124	gezelter	2204	while (currentSnapshot_->getTime() < runTime) {
125	gezelter	1930
126	gezelter	2204	preStep();
127	chrisfen	2101
128	gezelter	2204	integrateStep();
129	chrisfen	2101
130	gezelter	2204	postStep();
131	chrisfen	2101
132	gezelter	2204	}
133	chrisfen	2101
134	gezelter	2204	finalize();
135	chrisfen	2101
136	gezelter	2204	}
137	gezelter	1930
138
139	gezelter	2204	void VelocityVerletIntegrator::preStep() {
140	tim	2759	RealType difference = currentSnapshot_->getTime() + dt - currStatus;
141	chrisfen	2101
142	gezelter	2204	if (difference > 0 \|\| fabs(difference) < oopse::epsilon) {
143			needPotential = true;
144			needStress = true;
145			}
146	chrisfen	2101	}
147	gezelter	1930
148	gezelter	2204	void VelocityVerletIntegrator::postStep() {
149	tim	2448
150	gezelter	2204	//save snapshot
151			info_->getSnapshotManager()->advance();
152	chrisfen	2101
153	gezelter	2204	//increase time
154			currentSnapshot_->increaseTime(dt);
155	chrisfen	2879
156	gezelter	2204	if (needVelocityScaling) {
157			if (currentSnapshot_->getTime() >= currThermal) {
158			velocitizer_->velocitize(targetScalingTemp);
159			currThermal += thermalTime;
160			}
161	chrisfen	2101	}
162	gezelter	3488	if (useRNEMD) {
163			if (currentSnapshot_->getTime() >= currRNEMD) {
164			rnemd_->doSwap();
165			currRNEMD += RNEMD_swapTime;
166			}
167			}
168
169	gezelter	2204	if (currentSnapshot_->getTime() >= currSample) {
170			dumpWriter->writeDumpAndEor();
171	gezelter	3488
172	gezelter	2204	currSample += sampleTime;
173			}
174	gezelter	3488
175	gezelter	2204	if (currentSnapshot_->getTime() >= currStatus) {
176			//save statistics, before writeStat, we must save statistics
177			thermo.saveStat();
178			saveConservedQuantity();
179	skuang	3501
180	skuang	3498	if (simParams->getUseRNEMD())
181			rnemd_->getStatus();
182	skuang	3501
183			statWriter->writeStat(currentSnapshot_->statData);
184	gezelter	3488
185	gezelter	2204	needPotential = false;
186			needStress = false;
187			currStatus += statusTime;
188			}
189	gezelter	3488
190			if (needReset && currentSnapshot_->getTime() >= currReset) {
191			resetIntegrator();
192			currReset += resetTime;
193			}
194	gezelter	2204	}
195	gezelter	1930
196
197	gezelter	2204	void VelocityVerletIntegrator::finalize() {
198			dumpWriter->writeEor();
199	chrisfen	2101
200	gezelter	2204	delete dumpWriter;
201			delete statWriter;
202	chrisfen	2101
203	gezelter	2204	dumpWriter = NULL;
204			statWriter = NULL;
205	chrisfen	2101
206	gezelter	2204	}
207	gezelter	1930
208	gezelter	2204	void VelocityVerletIntegrator::integrateStep() {
209	chrisfen	2101
210	gezelter	2204	moveA();
211			calcForce(needPotential, needStress);
212			moveB();
213			}
214	gezelter	1930
215
216	gezelter	2204	void VelocityVerletIntegrator::calcForce(bool needPotential,
217			bool needStress) {
218			forceMan_->calcForces(needPotential, needStress);
219			}
220	gezelter	1930
221	gezelter	2204	DumpWriter* VelocityVerletIntegrator::createDumpWriter() {
222			return new DumpWriter(info_);
223			}
224	gezelter	1930
225	gezelter	2204	StatWriter* VelocityVerletIntegrator::createStatWriter() {
226	tim	2380
227			std::string statFileFormatString = simParams->getStatFileFormat();
228			StatsBitSet mask = parseStatFileFormat(statFileFormatString);
229	cli2	3520
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	2204	mask.set(Stats::VRAW);
236	cli2	3520
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	2204	mask.set(Stats::VHARM);
241			}
242	tim	2238
243	chrisfen	2879	if (simParams->havePrintPressureTensor() &&
244			simParams->getPrintPressureTensor()){
245	gezelter	3448	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	2238	}
255
256	chrisfen	2917	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	3488
262	gezelter	3448	if (simParams->haveTaggedAtomPair() && simParams->havePrintTaggedPairDistance()) {
263			if (simParams->getPrintTaggedPairDistance()) {
264			mask.set(Stats::TAGGED_PAIR_DISTANCE);
265			}
266			}
267	gezelter	3488
268			if (simParams->getUseRNEMD()) {
269			mask.set(Stats::RNEMD_SWAP_TOTAL);
270			}
271	gezelter	3448
272	chrisfen	2917
273	tim	2380	return new StatWriter(info_->getStatFileName(), mask);
274	chrisfen	2101	}
275	gezelter	1930
276
277			} //end namespace oopse