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, 9 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.
#	User	Rev	Content
1	chuckv	1293	/*
2	chuckv	1402	* Copyright (c) 2008, 2009, 2010 The University of Notre Dame. All Rights Reserved.
3	chuckv	1293	*
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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	chuckv	1293	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	chuckv	1293	* 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	gezelter	1390	*
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	chuckv	1293	*/
41			#include <fstream>
42			#include <iostream>
43			#include "integrators/SMIPDForceManager.hpp"
44	gezelter	1390	#include "utils/PhysicalConstants.hpp"
45	chuckv	1293	#include "math/ConvexHull.hpp"
46	chuckv	1402	#include "math/AlphaHull.hpp"
47	chuckv	1293	#include "math/Triangle.hpp"
48	gezelter	1398	#include "math/CholeskyDecomposition.hpp"
49	chuckv	1293
50	gezelter	1390	namespace OpenMD {
51	gezelter	1324
52	gezelter	1344	SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info) {
53	chuckv	1402
54	chuckv	1293	simParams = info->getSimParams();
55	chuckv	1325	veloMunge = new Velocitizer(info);
56	chuckv	1293
57			// Create Hull, Convex Hull for now, other options later.
58	gezelter	1324
59	chuckv	1402	stringToEnumMap_["ConvexHull"] = hullConvex;
60			stringToEnumMap_["AlphaShape"] = hullAlphaShape;
61			stringToEnumMap_["Unknown"] = hullUnknown;
62	chuckv	1293
63	chuckv	1402	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	chuckv	1293	/* Check that the simulation has target pressure and target
86	chuckv	1325	temperature set */
87	chuckv	1323
88	chuckv	1293	if (!simParams->haveTargetTemp()) {
89	gezelter	1324	sprintf(painCave.errMsg,
90			"SMIPDynamics error: You can't use the SMIPD integrator\n"
91	gezelter	1357	"\twithout a targetTemp (K)!\n");
92	chuckv	1293	painCave.isFatal = 1;
93	gezelter	1390	painCave.severity = OPENMD_ERROR;
94	chuckv	1293	simError();
95			} else {
96			targetTemp_ = simParams->getTargetTemp();
97			}
98	chuckv	1323
99	chuckv	1293	if (!simParams->haveTargetPressure()) {
100	gezelter	1324	sprintf(painCave.errMsg,
101			"SMIPDynamics error: You can't use the SMIPD integrator\n"
102	gezelter	1357	"\twithout a targetPressure (atm)!\n");
103	chuckv	1293	painCave.isFatal = 1;
104			simError();
105			} else {
106	gezelter	1324	// Convert pressure from atm -> amu/(fs^2*Ang)
107			targetPressure_ = simParams->getTargetPressure() /
108	gezelter	1390	PhysicalConstants::pressureConvert;
109	chuckv	1293	}
110
111			if (simParams->getUsePeriodicBoundaryConditions()) {
112	gezelter	1324	sprintf(painCave.errMsg,
113			"SMIPDynamics error: You can't use the SMIPD integrator\n"
114	gezelter	1357	"\twith periodic boundary conditions!\n");
115	chuckv	1293	painCave.isFatal = 1;
116			simError();
117			}
118	gezelter	1324
119	gezelter	1398	if (!simParams->haveViscosity()) {
120	gezelter	1324	sprintf(painCave.errMsg,
121			"SMIPDynamics error: You can't use the SMIPD integrator\n"
122	gezelter	1398	"\twithout a viscosity!\n");
123	chuckv	1316	painCave.isFatal = 1;
124	gezelter	1390	painCave.severity = OPENMD_ERROR;
125	chuckv	1316	simError();
126	chuckv	1323	}else{
127	gezelter	1398	viscosity_ = simParams->getViscosity();
128	chuckv	1316	}
129	gezelter	1324
130	chuckv	1323	dt_ = simParams->getDt();
131	gezelter	1324
132	gezelter	1390	variance_ = 2.0 * PhysicalConstants::kb * targetTemp_ / dt_;
133	chuckv	1307
134	gezelter	1344	// Build a vector of integrable objects to determine if the are
135			// surface atoms
136	chuckv	1293	Molecule* mol;
137			StuntDouble* integrableObject;
138			SimInfo::MoleculeIterator i;
139	gezelter	1344	Molecule::IntegrableObjectIterator j;
140	chuckv	1323
141	gezelter	1324	for (mol = info_->beginMolecule(i); mol != NULL;
142			mol = info_->nextMolecule(i)) {
143			for (integrableObject = mol->beginIntegrableObject(j);
144			integrableObject != NULL;
145	chuckv	1293	integrableObject = mol->nextIntegrableObject(j)) {
146			localSites_.push_back(integrableObject);
147			}
148	gezelter	1324	}
149	chuckv	1293	}
150	chuckv	1325
151	chuckv	1293	void SMIPDForceManager::postCalculation(bool needStress){
152			SimInfo::MoleculeIterator i;
153			Molecule::IntegrableObjectIterator j;
154			Molecule* mol;
155			StuntDouble* integrableObject;
156	chuckv	1325
157	gezelter	1324	// Compute surface Mesh
158	chuckv	1293	surfaceMesh_->computeHull(localSites_);
159	chuckv	1325
160	gezelter	1324	// Get total area and number of surface stunt doubles
161			RealType area = surfaceMesh_->getArea();
162	chuckv	1308	std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
163	chuckv	1307	int nTriangles = sMesh.size();
164
165	gezelter	1344	// Generate all of the necessary random forces
166	gezelter	1398	std::vector<Vector3d> randNums = genTriangleForces(nTriangles, variance_);
167	gezelter	1344
168	chuckv	1325	// Loop over the mesh faces and apply external pressure to each
169			// of the faces
170	chuckv	1308	std::vector<Triangle>::iterator face;
171	chuckv	1302	std::vector<StuntDouble*>::iterator vertex;
172	gezelter	1344	int thisFacet = 0;
173	chuckv	1293	for (face = sMesh.begin(); face != sMesh.end(); ++face){
174	chuckv	1308	Triangle thisTriangle = *face;
175			std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
176	chuckv	1316	RealType thisArea = thisTriangle.getArea();
177	chuckv	1308	Vector3d unitNormal = thisTriangle.getNormal();
178	chuckv	1302	unitNormal.normalize();
179	chuckv	1308	Vector3d centroid = thisTriangle.getCentroid();
180	gezelter	1344	Vector3d facetVel = thisTriangle.getFacetVelocity();
181	gezelter	1355	RealType thisMass = thisTriangle.getFacetMass();
182	gezelter	1398	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	gezelter	1344
191	gezelter	1398	Vector3d randomForce = S * randNums[thisFacet++];
192			Vector3d dragForce = -hydroTensor * facetVel;
193	gezelter	1377
194	gezelter	1398	Vector3d langevinForce = (extPressure + randomForce + dragForce);
195	gezelter	1355
196	chuckv	1375	// Apply triangle force to stuntdouble vertices
197	chuckv	1325	for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
198	chuckv	1375	if ((*vertex) != NULL){
199	gezelter	1344	Vector3d vertexForce = langevinForce / 3.0;
200			(*vertex)->addFrc(vertexForce);
201	chuckv	1307	}
202	chuckv	1293	}
203	chuckv	1302	}
204	gezelter	1355
205	chuckv	1325	veloMunge->removeComDrift();
206			veloMunge->removeAngularDrift();
207	gezelter	1355
208	chuckv	1302	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
209	chuckv	1316	currSnapshot->setVolume(surfaceMesh_->getVolume());
210	chuckv	1302	ForceManager::postCalculation(needStress);
211	chuckv	1293	}
212	gezelter	1344
213	gezelter	1355
214	gezelter	1398	std::vector<Vector3d> SMIPDForceManager::genTriangleForces(int nTriangles,
215	gezelter	1344	RealType variance)
216			{
217	chuckv	1325
218	gezelter	1344	// zero fill the random vector before starting:
219	gezelter	1398	std::vector<Vector3d> gaussRand;
220	gezelter	1344	gaussRand.resize(nTriangles);
221	gezelter	1398	std::fill(gaussRand.begin(), gaussRand.end(), V3Zero);
222	chuckv	1325
223	gezelter	1344	#ifdef IS_MPI
224			if (worldRank == 0) {
225			#endif
226			for (int i = 0; i < nTriangles; i++) {
227	gezelter	1398	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	gezelter	1344	}
231			#ifdef IS_MPI
232			}
233			#endif
234	chuckv	1325
235	gezelter	1344	// push these out to the other processors
236
237			#ifdef IS_MPI
238			if (worldRank == 0) {
239	gezelter	1398	MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
240	gezelter	1344	} else {
241	gezelter	1398	MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
242	gezelter	1344	}
243			#endif
244
245			return gaussRand;
246			}
247	chuckv	1293	}