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

    if (!simParams->haveViscosity()) {
      sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
      painCave.isFatal = 1;
      painCave.severity = OOPSE_ERROR;
      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 */
   

    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/3.14;

      RealType f_normal = simParams->getViscosity()*hydroLength*1.439326479e4;
      RealType extPressure = -(targetPressure_ * thisArea)/OOPSEConstant::energyConvert;
      RealType randomForce = randNums[thisNumber] * sqrt(2.0 * f_normal * OOPSEConstant::kb*targetTemp_/simParams->getDt());
      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 " << randomForce << " dragForce " << dragForce <<  " hydro  " << hydroLength << std::endl;


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

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