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 {
      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();
    } 


    //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/
    */

    // Build the hydroProp map:
    std::map<std::string, HydroProp*> hydroPropMap;

    Molecule* mol;
    StuntDouble* integrableObject;
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;              
    bool needHydroPropFile = false;
    
    for (mol = info->beginMolecule(i); mol != NULL; 
         mol = info->nextMolecule(i)) {
      for (integrableObject = mol->beginIntegrableObject(j); 
           integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(j)) {
        
        if (integrableObject->isRigidBody()) {
          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
          if (rb->getNumAtoms() > 1) needHydroPropFile = true;
        }
        
      }
    }
        

    if (needHydroPropFile) {               
      if (simParams->haveHydroPropFile()) {
        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
      } else {               
        sprintf( painCave.errMsg,
                 "HydroPropFile must be set to a file name if SMIPDynamics\n"
                 "\tis specified for rigidBodies which contain more than one atom\n"
                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n"
                 "\tFor use with SMIPD, the default viscosity in Hydro should be\n"
                 "\tset to 1.0 because the friction and random forces will be\n"
                 "\tdynamically re-set assuming this is true.\n");
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();  
      }      

      for (mol = info->beginMolecule(i); mol != NULL; 
           mol = info->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); 
             integrableObject != NULL;
             integrableObject = mol->nextIntegrableObject(j)) {

          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
          if (iter != hydroPropMap.end()) {
            hydroProps_.push_back(iter->second);
          } else {
            sprintf( painCave.errMsg,
                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
            painCave.severity = OOPSE_ERROR;
            painCave.isFatal = 1;
            simError();  
          }        
        }
      }
    } else {
      
      std::map<std::string, HydroProp*> hydroPropMap;
      for (mol = info->beginMolecule(i); mol != NULL; 
           mol = info->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); 
             integrableObject != NULL;
             integrableObject = mol->nextIntegrableObject(j)) {
          Shape* currShape = NULL;

          if (integrableObject->isAtom()){
            Atom* atom = static_cast<Atom*>(integrableObject);
            AtomType* atomType = atom->getAtomType();
            if (atomType->isGayBerne()) {
              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);              
              GenericData* data = dAtomType->getPropertyByName("GayBerne");
              if (data != NULL) {
                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
                
                if (gayBerneData != NULL) {  
                  GayBerneParam gayBerneParam = gayBerneData->getData();
                  currShape = new Ellipsoid(V3Zero, 
                                            gayBerneParam.GB_l / 2.0, 
                                            gayBerneParam.GB_d / 2.0, 
                                            Mat3x3d::identity());
                } else {
                  sprintf( painCave.errMsg,
                           "Can not cast GenericData to GayBerneParam\n");
                  painCave.severity = OOPSE_ERROR;
                  painCave.isFatal = 1;
                  simError();   
                }
              } else {
                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
                painCave.severity = OOPSE_ERROR;
                painCave.isFatal = 1;
                simError();    
              }
            } else {
              if (atomType->isLennardJones()){
                GenericData* data = atomType->getPropertyByName("LennardJones");
                if (data != NULL) {
                  LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
                  if (ljData != NULL) {
                    LJParam ljParam = ljData->getData();
                    currShape = new Sphere(atom->getPos(), 2.0);
                  } else {
                    sprintf( painCave.errMsg,
                             "Can not cast GenericData to LJParam\n");
                    painCave.severity = OOPSE_ERROR;
                    painCave.isFatal = 1;
                    simError();          
                  }       
                }
              } else {
                int aNum = etab.GetAtomicNum((atom->getType()).c_str());
                if (aNum != 0) {
                  currShape = new Sphere(atom->getPos(), 2.0);
                } else {
                  sprintf( painCave.errMsg,
                           "Could not find atom type in default element.txt\n");
                  painCave.severity = OOPSE_ERROR;
                  painCave.isFatal = 1;
                  simError();          
                }
              }
            }
          }
          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
          if (iter != hydroPropMap.end()) 
            hydroProps_.push_back(iter->second);
          else {
            currHydroProp->complete();
            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
            hydroProps_.push_back(currHydroProp);
          }
        }
      }
    }

    /* 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 */
   
    RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*area/nTriangles)
       /OOPSEConstant::energyConvert;

    gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * area * area * simParams->getDt()) /
      (4.0 * nTriangles * nTriangles* OOPSEConstant::kB*simParams->getTargetTemp()) /OOPSEConstant::energyConvert;

    std::vector<RealType>  randNums = genTriangleForces(nTriangles, sigma_t);

    /* 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();
      
      /* Get Random Force */
      Vector3d unitNormal = thisTriangle.getNormal();
      unitNormal.normalize();
      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
      Vector3d centroid = thisTriangle.getCentroid();

      Vector3d langevinForce = randomForce - gamma_t_*thisTriangle.getFacetVelocity();
      
      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
        if ((*vertex) != NULL){
          // mass = integrableObject->getMass();
          Vector3d vertexForce = langevinForce/3.0;
          (*vertex)->addFrc(vertexForce);
          if (integrableObject->isDirectional()){
            Vector3d vertexPos = (*vertex)->getPos();
            Vector3d vertexCentroidVector = vertexPos - centroid;
            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
          }
        }  
      }
    } 

    /* Now loop over all surface particles and apply the drag*/
    /*
    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
      integrableObject = *vertex;
      mass = integrableObject->getMass();
      if (integrableObject->isDirectional()){
        
        // preliminaries for directional objects:
        
        A = integrableObject->getA();
        Atrans = A.transpose();
        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
        //apply random force and torque at center of resistance
        Mat3x3d I = integrableObject->getI();
        Vector3d omegaBody;
        
        // What remains contains velocity explicitly, but the velocity required
        // is at the full step: v(t + h), while we have initially the velocity
        // at the half step: v(t + h/2).  We need to iterate to converge the
        // friction force and friction torque vectors.
        
        // this is the velocity at the half-step:
        
        Vector3d vel =integrableObject->getVel();
        Vector3d angMom = integrableObject->getJ();
        
        //estimate velocity at full-step using everything but friction forces:           
        
        frc = integrableObject->getFrc();
        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
        
        Tb = integrableObject->lab2Body(integrableObject->getTrq());
        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                             
        
        Vector3d omegaLab;
        Vector3d vcdLab;
        Vector3d vcdBody;
        Vector3d frictionForceBody;
        Vector3d frictionForceLab(0.0);
        Vector3d oldFFL;  // used to test for convergence
        Vector3d frictionTorqueBody(0.0);
        Vector3d oldFTB;  // used to test for convergence
        Vector3d frictionTorqueLab;
        RealType fdot;
        RealType tdot;

        //iteration starts here:
        
        for (int k = 0; k < maxIterNum_; k++) {
          
          if (integrableObject->isLinear()) {
            int linearAxis = integrableObject->linearAxis();
            int l = (linearAxis +1 )%3;
            int m = (linearAxis +2 )%3;
            omegaBody[l] = angMomStep[l] /I(l, l);
            omegaBody[m] = angMomStep[m] /I(m, m);
            
          } else {
            omegaBody[0] = angMomStep[0] /I(0, 0);
            omegaBody[1] = angMomStep[1] /I(1, 1);
            omegaBody[2] = angMomStep[2] /I(2, 2);
          }
          
          omegaLab = Atrans * omegaBody;
          
          // apply friction force and torque at center of resistance
          
          vcdLab = velStep + cross(omegaLab, rcrLab);       
          vcdBody = A * vcdLab;
          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
          oldFFL = frictionForceLab;
          frictionForceLab = Atrans * frictionForceBody;
          oldFTB = frictionTorqueBody;
          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
          frictionTorqueLab = Atrans * frictionTorqueBody;
          
          // re-estimate velocities at full-step using friction forces:
              
          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);

          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
              
          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
          
          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
            break; // iteration ends here
        }
        
        integrableObject->addFrc(frictionForceLab);
        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));

            
      } else {
        //spherical atom

        // What remains contains velocity explicitly, but the velocity required
        // is at the full step: v(t + h), while we have initially the velocity
        // at the half step: v(t + h/2).  We need to iterate to converge the
        // friction force vector.
        
        // this is the velocity at the half-step:
        
        Vector3d vel =integrableObject->getVel();
        
        //estimate velocity at full-step using everything but friction forces:           
        
        frc = integrableObject->getFrc();
        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
        
        Vector3d frictionForce(0.0);
        Vector3d oldFF;  // used to test for convergence
        RealType fdot;
        
        //iteration starts here:
        
        for (int k = 0; k < maxIterNum_; k++) {
          
          oldFF = frictionForce;                            
          frictionForce = -hydroProps_[index]->getXitt() * velStep;
          //frictionForce = -gamma_t*velStep;
          // re-estimate velocities at full-step using friction forces:
          
          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
          
          // check for convergence (if the vector has converged, fdot will be 1.0):
          
          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
          
          if (fabs(1.0 - fdot) <= forceTolerance_)
            break; // iteration ends here
        }
        
        integrableObject->addFrc(frictionForce);
        
        
      }
  
      
  }
    */
    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));     
      }
#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:	3465
Committed:	Tue Oct 21 16:44:00 2008 UTC (16 years, 6 months ago) by chuckv
File size:	20379 byte(s)
Log Message:	Fixed Memory Leak in ConvexHull.
#	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	targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
84	}
85
86
87	if (simParams->getUsePeriodicBoundaryConditions()) {
88	sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
89	" with periodic boundary conditions !\n");
90
91	painCave.isFatal = 1;
92	simError();
93	}
94
95
96
97
98
99	//Compute initial hull
100	/*
101	surfaceMesh_->computeHull(localSites_);
102	Area0_ = surfaceMesh_->getArea();
103	int nTriangles = surfaceMesh_->getNMeshElements();
104	// variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
105	gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
106	(4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
107	//RealType eta0 = gamma_0/
108	*/
109
110	// Build the hydroProp map:
111	std::map<std::string, HydroProp*> hydroPropMap;
112
113	Molecule* mol;
114	StuntDouble* integrableObject;
115	SimInfo::MoleculeIterator i;
116	Molecule::IntegrableObjectIterator j;
117	bool needHydroPropFile = false;
118
119	for (mol = info->beginMolecule(i); mol != NULL;
120	mol = info->nextMolecule(i)) {
121	for (integrableObject = mol->beginIntegrableObject(j);
122	integrableObject != NULL;
123	integrableObject = mol->nextIntegrableObject(j)) {
124
125	if (integrableObject->isRigidBody()) {
126	RigidBody* rb = static_cast<RigidBody*>(integrableObject);
127	if (rb->getNumAtoms() > 1) needHydroPropFile = true;
128	}
129
130	}
131	}
132
133
134	if (needHydroPropFile) {
135	if (simParams->haveHydroPropFile()) {
136	hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
137	} else {
138	sprintf( painCave.errMsg,
139	"HydroPropFile must be set to a file name if SMIPDynamics\n"
140	"\tis specified for rigidBodies which contain more than one atom\n"
141	"\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n"
142	"\tFor use with SMIPD, the default viscosity in Hydro should be\n"
143	"\tset to 1.0 because the friction and random forces will be\n"
144	"\tdynamically re-set assuming this is true.\n");
145	painCave.severity = OOPSE_ERROR;
146	painCave.isFatal = 1;
147	simError();
148	}
149
150	for (mol = info->beginMolecule(i); mol != NULL;
151	mol = info->nextMolecule(i)) {
152	for (integrableObject = mol->beginIntegrableObject(j);
153	integrableObject != NULL;
154	integrableObject = mol->nextIntegrableObject(j)) {
155
156	std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
157	if (iter != hydroPropMap.end()) {
158	hydroProps_.push_back(iter->second);
159	} else {
160	sprintf( painCave.errMsg,
161	"Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
162	painCave.severity = OOPSE_ERROR;
163	painCave.isFatal = 1;
164	simError();
165	}
166	}
167	}
168	} else {
169
170	std::map<std::string, HydroProp*> hydroPropMap;
171	for (mol = info->beginMolecule(i); mol != NULL;
172	mol = info->nextMolecule(i)) {
173	for (integrableObject = mol->beginIntegrableObject(j);
174	integrableObject != NULL;
175	integrableObject = mol->nextIntegrableObject(j)) {
176	Shape* currShape = NULL;
177
178	if (integrableObject->isAtom()){
179	Atom* atom = static_cast<Atom*>(integrableObject);
180	AtomType* atomType = atom->getAtomType();
181	if (atomType->isGayBerne()) {
182	DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
183	GenericData* data = dAtomType->getPropertyByName("GayBerne");
184	if (data != NULL) {
185	GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
186
187	if (gayBerneData != NULL) {
188	GayBerneParam gayBerneParam = gayBerneData->getData();
189	currShape = new Ellipsoid(V3Zero,
190	gayBerneParam.GB_l / 2.0,
191	gayBerneParam.GB_d / 2.0,
192	Mat3x3d::identity());
193	} else {
194	sprintf( painCave.errMsg,
195	"Can not cast GenericData to GayBerneParam\n");
196	painCave.severity = OOPSE_ERROR;
197	painCave.isFatal = 1;
198	simError();
199	}
200	} else {
201	sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
202	painCave.severity = OOPSE_ERROR;
203	painCave.isFatal = 1;
204	simError();
205	}
206	} else {
207	if (atomType->isLennardJones()){
208	GenericData* data = atomType->getPropertyByName("LennardJones");
209	if (data != NULL) {
210	LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
211	if (ljData != NULL) {
212	LJParam ljParam = ljData->getData();
213	currShape = new Sphere(atom->getPos(), 2.0);
214	} else {
215	sprintf( painCave.errMsg,
216	"Can not cast GenericData to LJParam\n");
217	painCave.severity = OOPSE_ERROR;
218	painCave.isFatal = 1;
219	simError();
220	}
221	}
222	} else {
223	int aNum = etab.GetAtomicNum((atom->getType()).c_str());
224	if (aNum != 0) {
225	currShape = new Sphere(atom->getPos(), 2.0);
226	} else {
227	sprintf( painCave.errMsg,
228	"Could not find atom type in default element.txt\n");
229	painCave.severity = OOPSE_ERROR;
230	painCave.isFatal = 1;
231	simError();
232	}
233	}
234	}
235	}
236	HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
237	std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
238	if (iter != hydroPropMap.end())
239	hydroProps_.push_back(iter->second);
240	else {
241	currHydroProp->complete();
242	hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
243	hydroProps_.push_back(currHydroProp);
244	}
245	}
246	}
247	}
248
249	/* Compute hull first time through to get the area of t=0*/
250
251	//Build a vector of integrable objects to determine if the are surface atoms
252	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
253	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
254	integrableObject = mol->nextIntegrableObject(j)) {
255	localSites_.push_back(integrableObject);
256	}
257	}
258
259
260	}
261
262	std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
263	std::map<std::string, HydroProp*> props;
264	std::ifstream ifs(filename.c_str());
265	if (ifs.is_open()) {
266
267	}
268
269	const unsigned int BufferSize = 65535;
270	char buffer[BufferSize];
271	while (ifs.getline(buffer, BufferSize)) {
272	HydroProp* currProp = new HydroProp(buffer);
273	props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
274	}
275
276	return props;
277	}
278
279	void SMIPDForceManager::postCalculation(bool needStress){
280	SimInfo::MoleculeIterator i;
281	Molecule::IntegrableObjectIterator j;
282	Molecule* mol;
283	StuntDouble* integrableObject;
284	RealType mass;
285	Vector3d pos;
286	Vector3d frc;
287	Mat3x3d A;
288	Mat3x3d Atrans;
289	Vector3d Tb;
290	Vector3d ji;
291	unsigned int index = 0;
292	int fdf;
293
294	fdf = 0;
295
296	/Compute surface Mesh/
297	surfaceMesh_->computeHull(localSites_);
298
299	/* Get area and number of surface stunt doubles and compute new variance */
300	RealType area = surfaceMesh_->getArea();
301	int nSurfaceSDs = surfaceMesh_->getNs();
302
303
304	std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
305	int nTriangles = sMesh.size();
306
307
308
309	/* Compute variance for random forces */
310
311	RealType sigma_t = sqrt(NumericConstant::PI/2.0)((targetPressure_)area/nTriangles)
312	/OOPSEConstant::energyConvert;
313
314	gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * area * area * simParams->getDt()) /
315	(4.0 * nTriangles * nTriangles* OOPSEConstant::kB*simParams->getTargetTemp()) /OOPSEConstant::energyConvert;
316
317	std::vector<RealType> randNums = genTriangleForces(nTriangles, sigma_t);
318
319	/* Loop over the mesh faces and apply random force to each of the faces*/
320
321
322	std::vector<Triangle>::iterator face;
323	std::vector<StuntDouble*>::iterator vertex;
324	int thisNumber = 0;
325	for (face = sMesh.begin(); face != sMesh.end(); ++face){
326
327	Triangle thisTriangle = *face;
328	std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
329
330	/* Get Random Force */
331	Vector3d unitNormal = thisTriangle.getNormal();
332	unitNormal.normalize();
333	Vector3d randomForce = -randNums[thisNumber] * unitNormal;
334	Vector3d centroid = thisTriangle.getCentroid();
335
336	Vector3d langevinForce = randomForce - gamma_t_*thisTriangle.getFacetVelocity();
337
338	for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
339	if ((*vertex) != NULL){
340	// mass = integrableObject->getMass();
341	Vector3d vertexForce = langevinForce/3.0;
342	(*vertex)->addFrc(vertexForce);
343	if (integrableObject->isDirectional()){
344	Vector3d vertexPos = (*vertex)->getPos();
345	Vector3d vertexCentroidVector = vertexPos - centroid;
346	(*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
347	}
348	}
349	}
350	}
351
352	/* Now loop over all surface particles and apply the drag*/
353	/*
354	std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
355	for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
356	integrableObject = *vertex;
357	mass = integrableObject->getMass();
358	if (integrableObject->isDirectional()){
359
360	// preliminaries for directional objects:
361
362	A = integrableObject->getA();
363	Atrans = A.transpose();
364	Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();
365	//apply random force and torque at center of resistance
366	Mat3x3d I = integrableObject->getI();
367	Vector3d omegaBody;
368
369	// What remains contains velocity explicitly, but the velocity required
370	// is at the full step: v(t + h), while we have initially the velocity
371	// at the half step: v(t + h/2). We need to iterate to converge the
372	// friction force and friction torque vectors.
373
374	// this is the velocity at the half-step:
375
376	Vector3d vel =integrableObject->getVel();
377	Vector3d angMom = integrableObject->getJ();
378
379	//estimate velocity at full-step using everything but friction forces:
380
381	frc = integrableObject->getFrc();
382	Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
383
384	Tb = integrableObject->lab2Body(integrableObject->getTrq());
385	Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;
386
387	Vector3d omegaLab;
388	Vector3d vcdLab;
389	Vector3d vcdBody;
390	Vector3d frictionForceBody;
391	Vector3d frictionForceLab(0.0);
392	Vector3d oldFFL; // used to test for convergence
393	Vector3d frictionTorqueBody(0.0);
394	Vector3d oldFTB; // used to test for convergence
395	Vector3d frictionTorqueLab;
396	RealType fdot;
397	RealType tdot;
398
399	//iteration starts here:
400
401	for (int k = 0; k < maxIterNum_; k++) {
402
403	if (integrableObject->isLinear()) {
404	int linearAxis = integrableObject->linearAxis();
405	int l = (linearAxis +1 )%3;
406	int m = (linearAxis +2 )%3;
407	omegaBody[l] = angMomStep[l] /I(l, l);
408	omegaBody[m] = angMomStep[m] /I(m, m);
409
410	} else {
411	omegaBody[0] = angMomStep[0] /I(0, 0);
412	omegaBody[1] = angMomStep[1] /I(1, 1);
413	omegaBody[2] = angMomStep[2] /I(2, 2);
414	}
415
416	omegaLab = Atrans * omegaBody;
417
418	// apply friction force and torque at center of resistance
419
420	vcdLab = velStep + cross(omegaLab, rcrLab);
421	vcdBody = A * vcdLab;
422	frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
423	oldFFL = frictionForceLab;
424	frictionForceLab = Atrans * frictionForceBody;
425	oldFTB = frictionTorqueBody;
426	frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
427	frictionTorqueLab = Atrans * frictionTorqueBody;
428
429	// re-estimate velocities at full-step using friction forces:
430
431	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
432	angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
433
434	// check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
435
436	fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
437	tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
438
439	if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
440	break; // iteration ends here
441	}
442
443	integrableObject->addFrc(frictionForceLab);
444	integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
445
446
447	} else {
448	//spherical atom
449
450	// What remains contains velocity explicitly, but the velocity required
451	// is at the full step: v(t + h), while we have initially the velocity
452	// at the half step: v(t + h/2). We need to iterate to converge the
453	// friction force vector.
454
455	// this is the velocity at the half-step:
456
457	Vector3d vel =integrableObject->getVel();
458
459	//estimate velocity at full-step using everything but friction forces:
460
461	frc = integrableObject->getFrc();
462	Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
463
464	Vector3d frictionForce(0.0);
465	Vector3d oldFF; // used to test for convergence
466	RealType fdot;
467
468	//iteration starts here:
469
470	for (int k = 0; k < maxIterNum_; k++) {
471
472	oldFF = frictionForce;
473	frictionForce = -hydroProps_[index]->getXitt() * velStep;
474	//frictionForce = -gamma_t*velStep;
475	// re-estimate velocities at full-step using friction forces:
476
477	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
478
479	// check for convergence (if the vector has converged, fdot will be 1.0):
480
481	fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
482
483	if (fabs(1.0 - fdot) <= forceTolerance_)
484	break; // iteration ends here
485	}
486
487	integrableObject->addFrc(frictionForce);
488
489
490	}
491
492
493	}
494	*/
495	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
496	currSnapshot->setVolume(surfaceMesh_->getVolume());
497
498	ForceManager::postCalculation(needStress);
499	}
500
501	void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
502
503
504	Vector<RealType, 6> Z;
505	Vector<RealType, 6> generalForce;
506
507
508	Z[0] = randNumGen_.randNorm(0, variance);
509	Z[1] = randNumGen_.randNorm(0, variance);
510	Z[2] = randNumGen_.randNorm(0, variance);
511	Z[3] = randNumGen_.randNorm(0, variance);
512	Z[4] = randNumGen_.randNorm(0, variance);
513	Z[5] = randNumGen_.randNorm(0, variance);
514
515	generalForce = hydroProps_[index]->getS()*Z;
516
517	force[0] = generalForce[0];
518	force[1] = generalForce[1];
519	force[2] = generalForce[2];
520	torque[0] = generalForce[3];
521	torque[1] = generalForce[4];
522	torque[2] = generalForce[5];
523
524	}
525	std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
526
527	// zero fill the random vector before starting:
528	std::vector<RealType> gaussRand;
529	gaussRand.resize(nTriangles);
530	std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
531
532
533	#ifdef IS_MPI
534	if (worldRank == 0) {
535	#endif
536	for (int i = 0; i < nTriangles; i++) {
537	gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));
538	}
539	#ifdef IS_MPI
540	}
541	#endif
542
543	// push these out to the other processors
544
545	#ifdef IS_MPI
546	if (worldRank == 0) {
547	MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
548	} else {
549	MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
550	}
551	#endif
552
553	for (int i = 0; i < nTriangles; i++) {
554	gaussRand[i] = gaussRand[i] * variance;
555	}
556
557	return gaussRand;
558	}
559
560
561
562
563
564	}