src/integrators/SMIPDForceManager.cpp

/*
 * Copyright (c) 2008 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.
 */
#include <fstream> 
#include <iostream>
#include "integrators/SMIPDForceManager.hpp"
#include "math/CholeskyDecomposition.hpp"
#include "utils/OOPSEConstant.hpp"
#include "hydrodynamics/Sphere.hpp"
#include "hydrodynamics/Ellipsoid.hpp"
#include "utils/ElementsTable.hpp"
#include "math/ConvexHull.hpp"
#include "math/Triangle.hpp"


namespace oopse {
  
  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
    simParams = info->getSimParams();
    veloMunge = new Velocitizer(info);
    
    // Create Hull, Convex Hull for now, other options later.
    surfaceMesh_ = new ConvexHull();

    
    /* Check that the simulation has target pressure and target
       temperature set*/
    
    if (!simParams->haveTargetTemp()) {
      sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
      painCave.isFatal = 1;
      painCave.severity = OOPSE_ERROR;
      simError();
    } else {
      targetTemp_ = simParams->getTargetTemp();
    }
    
    if (!simParams->haveTargetPressure()) {
      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "   without a targetPressure!\n");
      
      painCave.isFatal = 1;
      simError();
    } else {
      /* Convert pressure from atm -> amu/(fs^2*Ang)*/
      targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
    }

   
    if (simParams->getUsePeriodicBoundaryConditions()) {
      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "   with periodic boundary conditions !\n");
      
      painCave.isFatal = 1;
      simError();
    } 

    if (!simParams->haveViscosity()) {
      sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
      painCave.isFatal = 1;
      painCave.severity = OOPSE_ERROR;
      simError();
    }else{
      viscosity_ = simParams->getViscosity();
    }

    dt_ = simParams->getDt();

    
    //Compute initial hull
    /*
    surfaceMesh_->computeHull(localSites_);
    Area0_ = surfaceMesh_->getArea();
    int nTriangles = surfaceMesh_->getNMeshElements();
    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
    //RealType eta0 = gamma_0/
    */


    Molecule* mol;
    StuntDouble* integrableObject;
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;              

    
    /* Compute hull first time through to get the area of t=0*/

    //Build a vector of integrable objects to determine if the are surface atoms 
    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);
      }
    }


  }  

  std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
    std::map<std::string, HydroProp*> props;
    std::ifstream ifs(filename.c_str());
    if (ifs.is_open()) {
      
    }
    
    const unsigned int BufferSize = 65535;
    char buffer[BufferSize];   
    while (ifs.getline(buffer, BufferSize)) {
      HydroProp* currProp = new HydroProp(buffer);
      props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
    }

    return props;
  }
   
  void SMIPDForceManager::postCalculation(bool needStress){
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    RealType mass;
    Vector3d pos;
    Vector3d frc;
    Mat3x3d A;
    Mat3x3d Atrans;
    Vector3d Tb;
    Vector3d ji;
    unsigned int index = 0;
    int fdf;
   
    fdf = 0;
  
    /*Compute surface Mesh*/
    surfaceMesh_->computeHull(localSites_);

    /* Get area and number of surface stunt doubles and compute new variance */
     RealType area = surfaceMesh_->getArea();
     int nSurfaceSDs = surfaceMesh_->getNs();

    
    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
    int nTriangles = sMesh.size();


     /* Compute variance for random forces */
    std::vector<RealType>  randNums = genTriangleForces(nTriangles, 1.0);

    /* Loop over the mesh faces and apply random force to each of the faces*/    
    std::vector<Triangle>::iterator face;
    std::vector<StuntDouble*>::iterator vertex;
    int thisNumber = 0;
    for (face = sMesh.begin(); face != sMesh.end(); ++face){
     
      Triangle thisTriangle = *face;
      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
      RealType thisArea = thisTriangle.getArea(); 
      // RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*thisArea) /OOPSEConstant::energyConvert;
      // gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * thisArea * thisArea * simParams->getDt()) /(4.0 * OOPSEConstant::kB*simParams->getTargetTemp());

      /* Get Random Force */
      Vector3d unitNormal = thisTriangle.getNormal();
      unitNormal.normalize();
      //Vector3d randomForce = -randNums[thisNumber] * sigma_t * unitNormal;
      Vector3d centroid = thisTriangle.getCentroid();
      Vector3d facetVel = thisTriangle.getFacetVelocity();
      RealType hydroLength = thisTriangle.getIncircleRadius()*2.0/NumericConstant::PI;

      RealType f_normal = viscosity_*hydroLength*OOPSEConstant::viscoConvert;
      RealType extPressure = -(targetPressure_ * thisArea)/OOPSEConstant::energyConvert;
      RealType randomForce = randNums[thisNumber++] * sqrt(2.0 * f_normal * OOPSEConstant::kb*targetTemp_/dt_);

      RealType dragForce = -f_normal * dot(facetVel, unitNormal);


      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
      
      //      Vector3d dragForce = - gamma_t_ * dot(facetVel, unitNormal) * unitNormal / OOPSEConstant::energyConvert;
      
      // std::cout << " " << randomForce << " " << f_normal <<   std::endl;

      /* Apply triangle force to stuntdouble vertices */
      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
        if ((*vertex) != NULL){
          // mass = integrableObject->getMass();
          Vector3d vertexForce = langevinForce/3.0;
          (*vertex)->addFrc(vertexForce);

          if ((*vertex)->isDirectional()){

            Vector3d vertexPos = (*vertex)->getPos();
            Vector3d vertexCentroidVector = vertexPos - centroid;
            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
          }
        }  
      }
    } 

    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
    currSnapshot->setVolume(surfaceMesh_->getVolume());    
    ForceManager::postCalculation(needStress);   
  }

  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {

    
    Vector<RealType, 6> Z;
    Vector<RealType, 6> generalForce;
    

    Z[0] = randNumGen_.randNorm(0, variance);
    Z[1] = randNumGen_.randNorm(0, variance);
    Z[2] = randNumGen_.randNorm(0, variance);
    Z[3] = randNumGen_.randNorm(0, variance);
    Z[4] = randNumGen_.randNorm(0, variance);
    Z[5] = randNumGen_.randNorm(0, variance);
     
    generalForce = hydroProps_[index]->getS()*Z;
    
    force[0] = generalForce[0];
    force[1] = generalForce[1];
    force[2] = generalForce[2];
    torque[0] = generalForce[3];
    torque[1] = generalForce[4];
    torque[2] = generalForce[5];
    
} 
  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {

    // zero fill the random vector before starting:
    std::vector<RealType> gaussRand;
    gaussRand.resize(nTriangles);
    std::fill(gaussRand.begin(), gaussRand.end(), 0.0);

    
#ifdef IS_MPI
    if (worldRank == 0) {
#endif
      for (int i = 0; i < nTriangles; i++) {
        //gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));     
        gaussRand[i] = randNumGen_.randNorm(0.0, 1.0);     
      }
#ifdef IS_MPI
    }
#endif

    // push these out to the other processors

#ifdef IS_MPI
    if (worldRank == 0) {
      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
    } else {
      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
    }
#endif
   
    for (int i = 0; i < nTriangles; i++) {
      gaussRand[i] = gaussRand[i] * variance;
    }
   
    return gaussRand;
  }


}
Revision:	3480
Committed:	Tue Nov 25 16:32:47 2008 UTC (16 years, 5 months ago) by chuckv
File size:	10696 byte(s)
Log Message:	Code clean-up. Removed old langevin code.
#	Content
1	/*
2	* Copyright (c) 2008 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	#include <fstream>
42	#include <iostream>
43	#include "integrators/SMIPDForceManager.hpp"
44	#include "math/CholeskyDecomposition.hpp"
45	#include "utils/OOPSEConstant.hpp"
46	#include "hydrodynamics/Sphere.hpp"
47	#include "hydrodynamics/Ellipsoid.hpp"
48	#include "utils/ElementsTable.hpp"
49	#include "math/ConvexHull.hpp"
50	#include "math/Triangle.hpp"
51
52
53	namespace oopse {
54
55	SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
56	simParams = info->getSimParams();
57	veloMunge = new Velocitizer(info);
58
59	// Create Hull, Convex Hull for now, other options later.
60	surfaceMesh_ = new ConvexHull();
61
62
63
64	/* Check that the simulation has target pressure and target
65	temperature set*/
66
67	if (!simParams->haveTargetTemp()) {
68	sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
69	painCave.isFatal = 1;
70	painCave.severity = OOPSE_ERROR;
71	simError();
72	} else {
73	targetTemp_ = simParams->getTargetTemp();
74	}
75
76	if (!simParams->haveTargetPressure()) {
77	sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
78	" without a targetPressure!\n");
79
80	painCave.isFatal = 1;
81	simError();
82	} else {
83	/* Convert pressure from atm -> amu/(fs^2Ang)/
84	targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
85	}
86
87
88	if (simParams->getUsePeriodicBoundaryConditions()) {
89	sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
90	" with periodic boundary conditions !\n");
91
92	painCave.isFatal = 1;
93	simError();
94	}
95
96	if (!simParams->haveViscosity()) {
97	sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
98	painCave.isFatal = 1;
99	painCave.severity = OOPSE_ERROR;
100	simError();
101	}else{
102	viscosity_ = simParams->getViscosity();
103	}
104
105	dt_ = simParams->getDt();
106
107
108
109	//Compute initial hull
110	/*
111	surfaceMesh_->computeHull(localSites_);
112	Area0_ = surfaceMesh_->getArea();
113	int nTriangles = surfaceMesh_->getNMeshElements();
114	// variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
115	gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
116	(4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
117	//RealType eta0 = gamma_0/
118	*/
119
120
121	Molecule* mol;
122	StuntDouble* integrableObject;
123	SimInfo::MoleculeIterator i;
124	Molecule::IntegrableObjectIterator j;
125
126
127
128	/* Compute hull first time through to get the area of t=0*/
129
130	//Build a vector of integrable objects to determine if the are surface atoms
131	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
132	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
133	integrableObject = mol->nextIntegrableObject(j)) {
134	localSites_.push_back(integrableObject);
135	}
136	}
137
138
139	}
140
141	std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
142	std::map<std::string, HydroProp*> props;
143	std::ifstream ifs(filename.c_str());
144	if (ifs.is_open()) {
145
146	}
147
148	const unsigned int BufferSize = 65535;
149	char buffer[BufferSize];
150	while (ifs.getline(buffer, BufferSize)) {
151	HydroProp* currProp = new HydroProp(buffer);
152	props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
153	}
154
155	return props;
156	}
157
158	void SMIPDForceManager::postCalculation(bool needStress){
159	SimInfo::MoleculeIterator i;
160	Molecule::IntegrableObjectIterator j;
161	Molecule* mol;
162	StuntDouble* integrableObject;
163	RealType mass;
164	Vector3d pos;
165	Vector3d frc;
166	Mat3x3d A;
167	Mat3x3d Atrans;
168	Vector3d Tb;
169	Vector3d ji;
170	unsigned int index = 0;
171	int fdf;
172
173	fdf = 0;
174
175	/Compute surface Mesh/
176	surfaceMesh_->computeHull(localSites_);
177
178	/* Get area and number of surface stunt doubles and compute new variance */
179	RealType area = surfaceMesh_->getArea();
180	int nSurfaceSDs = surfaceMesh_->getNs();
181
182
183	std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
184	int nTriangles = sMesh.size();
185
186
187
188	/* Compute variance for random forces */
189	std::vector<RealType> randNums = genTriangleForces(nTriangles, 1.0);
190
191	/* Loop over the mesh faces and apply random force to each of the faces*/
192	std::vector<Triangle>::iterator face;
193	std::vector<StuntDouble*>::iterator vertex;
194	int thisNumber = 0;
195	for (face = sMesh.begin(); face != sMesh.end(); ++face){
196
197	Triangle thisTriangle = *face;
198	std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
199	RealType thisArea = thisTriangle.getArea();
200	// RealType sigma_t = sqrt(NumericConstant::PI/2.0)((targetPressure_)thisArea) /OOPSEConstant::energyConvert;
201	// gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * thisArea * thisArea * simParams->getDt()) /(4.0 * OOPSEConstant::kB*simParams->getTargetTemp());
202
203	/* Get Random Force */
204	Vector3d unitNormal = thisTriangle.getNormal();
205	unitNormal.normalize();
206	//Vector3d randomForce = -randNums[thisNumber] * sigma_t * unitNormal;
207	Vector3d centroid = thisTriangle.getCentroid();
208	Vector3d facetVel = thisTriangle.getFacetVelocity();
209	RealType hydroLength = thisTriangle.getIncircleRadius()*2.0/NumericConstant::PI;
210
211	RealType f_normal = viscosity_hydroLengthOOPSEConstant::viscoConvert;
212	RealType extPressure = -(targetPressure_ * thisArea)/OOPSEConstant::energyConvert;
213	RealType randomForce = randNums[thisNumber++] * sqrt(2.0 * f_normal * OOPSEConstant::kb*targetTemp_/dt_);
214
215	RealType dragForce = -f_normal * dot(facetVel, unitNormal);
216
217
218	Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
219
220	// Vector3d dragForce = - gamma_t_ * dot(facetVel, unitNormal) * unitNormal / OOPSEConstant::energyConvert;
221
222	// std::cout << " " << randomForce << " " << f_normal << std::endl;
223
224	/* Apply triangle force to stuntdouble vertices */
225	for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
226	if ((*vertex) != NULL){
227	// mass = integrableObject->getMass();
228	Vector3d vertexForce = langevinForce/3.0;
229	(*vertex)->addFrc(vertexForce);
230
231	if ((*vertex)->isDirectional()){
232
233	Vector3d vertexPos = (*vertex)->getPos();
234	Vector3d vertexCentroidVector = vertexPos - centroid;
235	(*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
236	}
237	}
238	}
239	}
240
241	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
242	currSnapshot->setVolume(surfaceMesh_->getVolume());
243	ForceManager::postCalculation(needStress);
244	}
245
246	void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
247
248
249	Vector<RealType, 6> Z;
250	Vector<RealType, 6> generalForce;
251
252
253	Z[0] = randNumGen_.randNorm(0, variance);
254	Z[1] = randNumGen_.randNorm(0, variance);
255	Z[2] = randNumGen_.randNorm(0, variance);
256	Z[3] = randNumGen_.randNorm(0, variance);
257	Z[4] = randNumGen_.randNorm(0, variance);
258	Z[5] = randNumGen_.randNorm(0, variance);
259
260	generalForce = hydroProps_[index]->getS()*Z;
261
262	force[0] = generalForce[0];
263	force[1] = generalForce[1];
264	force[2] = generalForce[2];
265	torque[0] = generalForce[3];
266	torque[1] = generalForce[4];
267	torque[2] = generalForce[5];
268
269	}
270	std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
271
272	// zero fill the random vector before starting:
273	std::vector<RealType> gaussRand;
274	gaussRand.resize(nTriangles);
275	std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
276
277
278
279	#ifdef IS_MPI
280	if (worldRank == 0) {
281	#endif
282	for (int i = 0; i < nTriangles; i++) {
283	//gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));
284	gaussRand[i] = randNumGen_.randNorm(0.0, 1.0);
285	}
286	#ifdef IS_MPI
287	}
288	#endif
289
290	// push these out to the other processors
291
292	#ifdef IS_MPI
293	if (worldRank == 0) {
294	MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
295	} else {
296	MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
297	}
298	#endif
299
300	for (int i = 0; i < nTriangles; i++) {
301	gaussRand[i] = gaussRand[i] * variance;
302	}
303
304	return gaussRand;
305	}
306
307
308	}