src/integrators/SMIPDForceManager.cpp

/*
 * Copyright (c) 2008, 2009, 2010 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. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. 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.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Vardeman & Gezelter, in progress (2009).                        
 */
#include <fstream> 
#include <iostream>
#include "integrators/SMIPDForceManager.hpp"
#include "utils/PhysicalConstants.hpp"
#include "math/ConvexHull.hpp"
#include "math/AlphaHull.hpp"
#include "math/Triangle.hpp"
#include "math/CholeskyDecomposition.hpp"

namespace OpenMD {

  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info) {
    
    simParams = info->getSimParams();
    veloMunge = new Velocitizer(info);
    
    // Create Hull, Convex Hull for now, other options later.
    
    stringToEnumMap_["ConvexHull"] = hullConvex;
    stringToEnumMap_["AlphaShape"] = hullAlphaShape;
    stringToEnumMap_["Unknown"] = hullUnknown;
    
    const std::string ht = simParams->getHULL_Method();
    
    
    
    std::map<std::string, HullTypeEnum>::iterator iter;
    iter = stringToEnumMap_.find(ht);
    hullType_ = (iter == stringToEnumMap_.end()) ? SMIPDForceManager::hullUnknown : iter->second;
    if (hullType_ == hullUnknown) {
      std::cerr << "WARNING! Hull Type Unknown!\n";
    }
    
    switch(hullType_) {
    case hullConvex :
      surfaceMesh_ = new ConvexHull();
      break;
    case hullAlphaShape :
      surfaceMesh_ = new AlphaHull(simParams->getAlpha());
      break;
    case hullUnknown :
    default :
      break;
    }
    /* Check that the simulation has target pressure and target
       temperature set */
    
    if (!simParams->haveTargetTemp()) {
      sprintf(painCave.errMsg, 
              "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "\twithout a targetTemp (K)!\n");      
      painCave.isFatal = 1;
      painCave.severity = OPENMD_ERROR;
      simError();
    } else {
      targetTemp_ = simParams->getTargetTemp();
    }
    
    if (!simParams->haveTargetPressure()) {
      sprintf(painCave.errMsg, 
              "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "\twithout a targetPressure (atm)!\n");      
      painCave.isFatal = 1;
      simError();
    } else {
      // Convert pressure from atm -> amu/(fs^2*Ang)
      targetPressure_ = simParams->getTargetPressure() / 
        PhysicalConstants::pressureConvert;
    }
   
    if (simParams->getUsePeriodicBoundaryConditions()) {
      sprintf(painCave.errMsg, 
              "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "\twith periodic boundary conditions!\n");    
      painCave.isFatal = 1;
      simError();
    } 
    
    if (!simParams->haveViscosity()) {
      sprintf(painCave.errMsg, 
              "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "\twithout a viscosity!\n");
      painCave.isFatal = 1;
      painCave.severity = OPENMD_ERROR;
      simError();
    }else{
      viscosity_ = simParams->getViscosity();
    }
    
    dt_ = simParams->getDt();
  
    variance_ = 2.0 * PhysicalConstants::kb * targetTemp_ / dt_;

    // Build a vector of integrable objects to determine if the are
    // surface atoms
    Molecule* mol;
    StuntDouble* integrableObject;
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;

    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {          
      for (integrableObject = mol->beginIntegrableObject(j); 
           integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(j)) {   
        localSites_.push_back(integrableObject);
      }
    }   
  }  
   
  void SMIPDForceManager::postCalculation(bool needStress){
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
  
    // Compute surface Mesh
    surfaceMesh_->computeHull(localSites_);

    // Get total area and number of surface stunt doubles
    RealType area = surfaceMesh_->getArea();
    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
    int nTriangles = sMesh.size();

    // Generate all of the necessary random forces
    std::vector<Vector3d>  randNums = genTriangleForces(nTriangles, variance_);

    // Loop over the mesh faces and apply external pressure to each 
    // of the faces
    std::vector<Triangle>::iterator face;
    std::vector<StuntDouble*>::iterator vertex;
    int thisFacet = 0;
    for (face = sMesh.begin(); face != sMesh.end(); ++face){
      Triangle thisTriangle = *face;
      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
      RealType thisArea = thisTriangle.getArea(); 
      Vector3d unitNormal = thisTriangle.getNormal();
      unitNormal.normalize();
      Vector3d centroid = thisTriangle.getCentroid();
      Vector3d facetVel = thisTriangle.getFacetVelocity();
      RealType thisMass = thisTriangle.getFacetMass();
      Mat3x3d hydroTensor = thisTriangle.computeHydrodynamicTensor(viscosity_);
      
      hydroTensor *= PhysicalConstants::viscoConvert;
      Mat3x3d S;
      CholeskyDecomposition(hydroTensor, S);
      
      Vector3d extPressure = -unitNormal * (targetPressure_ * thisArea) /
        PhysicalConstants::energyConvert;

      Vector3d randomForce = S * randNums[thisFacet++];
      Vector3d dragForce = -hydroTensor * facetVel;

      Vector3d langevinForce = (extPressure + randomForce + dragForce);
      
      // Apply triangle force to stuntdouble vertices
      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
        if ((*vertex) != NULL){
          Vector3d vertexForce = langevinForce / 3.0;
          (*vertex)->addFrc(vertexForce);          
        }  
      }
    } 
    
    veloMunge->removeComDrift();
    veloMunge->removeAngularDrift();
    
    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
    currSnapshot->setVolume(surfaceMesh_->getVolume());    
    ForceManager::postCalculation(needStress);   
  }
  
  
  std::vector<Vector3d> SMIPDForceManager::genTriangleForces(int nTriangles, 
                                                             RealType variance)
  {
    
    // zero fill the random vector before starting:
    std::vector<Vector3d> gaussRand;
    gaussRand.resize(nTriangles);
    std::fill(gaussRand.begin(), gaussRand.end(), V3Zero);
    
#ifdef IS_MPI
    if (worldRank == 0) {
#endif
      for (int i = 0; i < nTriangles; i++) {
        gaussRand[i][0] = randNumGen_.randNorm(0.0, variance);
        gaussRand[i][1] = randNumGen_.randNorm(0.0, variance);
        gaussRand[i][2] = randNumGen_.randNorm(0.0, variance);
      }
#ifdef IS_MPI
    }
#endif
    
    // push these out to the other processors
    
#ifdef IS_MPI
    if (worldRank == 0) {
      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
    } else {
      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
    }
#endif
    
    return gaussRand;
  }
}
Revision:	1402
Committed:	Fri Jan 8 17:15:27 2010 UTC (15 years, 3 months ago) by chuckv
File size:	8849 byte(s)
Log Message:	Added preliminary code for Alpha Hull calculation using qhull. Added preliminary support to SMIPD to support Alpha Hull. Alpha Hull does not yet add the correct things to triangle to be returned to SMPID. Preliminary changes for shadow hamiltonian integrator. Chages to md files so they will work in openMD.
#	Content
1	/*
2	* Copyright (c) 2008, 2009, 2010 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. Redistributions of source code must retain the above copyright
10	* notice, this list of conditions and the following disclaimer.
11	*
12	* 2. Redistributions in binary form must reproduce the above copyright
13	* notice, this list of conditions and the following disclaimer in the
14	* documentation and/or other materials provided with the
15	* distribution.
16	*
17	* This software is provided "AS IS," without a warranty of any
18	* kind. All express or implied conditions, representations and
19	* warranties, including any implied warranty of merchantability,
20	* fitness for a particular purpose or non-infringement, are hereby
21	* excluded. The University of Notre Dame and its licensors shall not
22	* be liable for any damages suffered by licensee as a result of
23	* using, modifying or distributing the software or its
24	* derivatives. In no event will the University of Notre Dame or its
25	* licensors be liable for any lost revenue, profit or data, or for
26	* direct, indirect, special, consequential, incidental or punitive
27	* damages, however caused and regardless of the theory of liability,
28	* arising out of the use of or inability to use software, even if the
29	* University of Notre Dame has been advised of the possibility of
30	* such damages.
31	*
32	* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33	* research, please cite the appropriate papers when you publish your
34	* work. Good starting points are:
35	*
36	* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	* [4] Vardeman & Gezelter, in progress (2009).
40	*/
41	#include <fstream>
42	#include <iostream>
43	#include "integrators/SMIPDForceManager.hpp"
44	#include "utils/PhysicalConstants.hpp"
45	#include "math/ConvexHull.hpp"
46	#include "math/AlphaHull.hpp"
47	#include "math/Triangle.hpp"
48	#include "math/CholeskyDecomposition.hpp"
49
50	namespace OpenMD {
51
52	SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info) {
53
54	simParams = info->getSimParams();
55	veloMunge = new Velocitizer(info);
56
57	// Create Hull, Convex Hull for now, other options later.
58
59	stringToEnumMap_["ConvexHull"] = hullConvex;
60	stringToEnumMap_["AlphaShape"] = hullAlphaShape;
61	stringToEnumMap_["Unknown"] = hullUnknown;
62
63	const std::string ht = simParams->getHULL_Method();
64
65
66
67	std::map<std::string, HullTypeEnum>::iterator iter;
68	iter = stringToEnumMap_.find(ht);
69	hullType_ = (iter == stringToEnumMap_.end()) ? SMIPDForceManager::hullUnknown : iter->second;
70	if (hullType_ == hullUnknown) {
71	std::cerr << "WARNING! Hull Type Unknown!\n";
72	}
73
74	switch(hullType_) {
75	case hullConvex :
76	surfaceMesh_ = new ConvexHull();
77	break;
78	case hullAlphaShape :
79	surfaceMesh_ = new AlphaHull(simParams->getAlpha());
80	break;
81	case hullUnknown :
82	default :
83	break;
84	}
85	/* Check that the simulation has target pressure and target
86	temperature set */
87
88	if (!simParams->haveTargetTemp()) {
89	sprintf(painCave.errMsg,
90	"SMIPDynamics error: You can't use the SMIPD integrator\n"
91	"\twithout a targetTemp (K)!\n");
92	painCave.isFatal = 1;
93	painCave.severity = OPENMD_ERROR;
94	simError();
95	} else {
96	targetTemp_ = simParams->getTargetTemp();
97	}
98
99	if (!simParams->haveTargetPressure()) {
100	sprintf(painCave.errMsg,
101	"SMIPDynamics error: You can't use the SMIPD integrator\n"
102	"\twithout a targetPressure (atm)!\n");
103	painCave.isFatal = 1;
104	simError();
105	} else {
106	// Convert pressure from atm -> amu/(fs^2*Ang)
107	targetPressure_ = simParams->getTargetPressure() /
108	PhysicalConstants::pressureConvert;
109	}
110
111	if (simParams->getUsePeriodicBoundaryConditions()) {
112	sprintf(painCave.errMsg,
113	"SMIPDynamics error: You can't use the SMIPD integrator\n"
114	"\twith periodic boundary conditions!\n");
115	painCave.isFatal = 1;
116	simError();
117	}
118
119	if (!simParams->haveViscosity()) {
120	sprintf(painCave.errMsg,
121	"SMIPDynamics error: You can't use the SMIPD integrator\n"
122	"\twithout a viscosity!\n");
123	painCave.isFatal = 1;
124	painCave.severity = OPENMD_ERROR;
125	simError();
126	}else{
127	viscosity_ = simParams->getViscosity();
128	}
129
130	dt_ = simParams->getDt();
131
132	variance_ = 2.0 * PhysicalConstants::kb * targetTemp_ / dt_;
133
134	// Build a vector of integrable objects to determine if the are
135	// surface atoms
136	Molecule* mol;
137	StuntDouble* integrableObject;
138	SimInfo::MoleculeIterator i;
139	Molecule::IntegrableObjectIterator j;
140
141	for (mol = info_->beginMolecule(i); mol != NULL;
142	mol = info_->nextMolecule(i)) {
143	for (integrableObject = mol->beginIntegrableObject(j);
144	integrableObject != NULL;
145	integrableObject = mol->nextIntegrableObject(j)) {
146	localSites_.push_back(integrableObject);
147	}
148	}
149	}
150
151	void SMIPDForceManager::postCalculation(bool needStress){
152	SimInfo::MoleculeIterator i;
153	Molecule::IntegrableObjectIterator j;
154	Molecule* mol;
155	StuntDouble* integrableObject;
156
157	// Compute surface Mesh
158	surfaceMesh_->computeHull(localSites_);
159
160	// Get total area and number of surface stunt doubles
161	RealType area = surfaceMesh_->getArea();
162	std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
163	int nTriangles = sMesh.size();
164
165	// Generate all of the necessary random forces
166	std::vector<Vector3d> randNums = genTriangleForces(nTriangles, variance_);
167
168	// Loop over the mesh faces and apply external pressure to each
169	// of the faces
170	std::vector<Triangle>::iterator face;
171	std::vector<StuntDouble*>::iterator vertex;
172	int thisFacet = 0;
173	for (face = sMesh.begin(); face != sMesh.end(); ++face){
174	Triangle thisTriangle = *face;
175	std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
176	RealType thisArea = thisTriangle.getArea();
177	Vector3d unitNormal = thisTriangle.getNormal();
178	unitNormal.normalize();
179	Vector3d centroid = thisTriangle.getCentroid();
180	Vector3d facetVel = thisTriangle.getFacetVelocity();
181	RealType thisMass = thisTriangle.getFacetMass();
182	Mat3x3d hydroTensor = thisTriangle.computeHydrodynamicTensor(viscosity_);
183
184	hydroTensor *= PhysicalConstants::viscoConvert;
185	Mat3x3d S;
186	CholeskyDecomposition(hydroTensor, S);
187
188	Vector3d extPressure = -unitNormal * (targetPressure_ * thisArea) /
189	PhysicalConstants::energyConvert;
190
191	Vector3d randomForce = S * randNums[thisFacet++];
192	Vector3d dragForce = -hydroTensor * facetVel;
193
194	Vector3d langevinForce = (extPressure + randomForce + dragForce);
195
196	// Apply triangle force to stuntdouble vertices
197	for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
198	if ((*vertex) != NULL){
199	Vector3d vertexForce = langevinForce / 3.0;
200	(*vertex)->addFrc(vertexForce);
201	}
202	}
203	}
204
205	veloMunge->removeComDrift();
206	veloMunge->removeAngularDrift();
207
208	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
209	currSnapshot->setVolume(surfaceMesh_->getVolume());
210	ForceManager::postCalculation(needStress);
211	}
212
213
214	std::vector<Vector3d> SMIPDForceManager::genTriangleForces(int nTriangles,
215	RealType variance)
216	{
217
218	// zero fill the random vector before starting:
219	std::vector<Vector3d> gaussRand;
220	gaussRand.resize(nTriangles);
221	std::fill(gaussRand.begin(), gaussRand.end(), V3Zero);
222
223	#ifdef IS_MPI
224	if (worldRank == 0) {
225	#endif
226	for (int i = 0; i < nTriangles; i++) {
227	gaussRand[i][0] = randNumGen_.randNorm(0.0, variance);
228	gaussRand[i][1] = randNumGen_.randNorm(0.0, variance);
229	gaussRand[i][2] = randNumGen_.randNorm(0.0, variance);
230	}
231	#ifdef IS_MPI
232	}
233	#endif
234
235	// push these out to the other processors
236
237	#ifdef IS_MPI
238	if (worldRank == 0) {
239	MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
240	} else {
241	MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
242	}
243	#endif
244
245	return gaussRand;
246	}
247	}