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