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

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


    // 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(), ljParam.sigma/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(), etab.GetVdwRad(aNum));
                } 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);
      }
    }

    surfaceMesh_->computeHull(localSites_);
    Area0_ = surfaceMesh_->getArea();
    //variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
    
  }  

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

     /* Compute variance for random forces */
    
     variance_ = sqrt(2.0*NumericConstant::PI)*((targetPressure_/OOPSEConstant::pressureConvert)*area/nSurfaceSDs)
       /OOPSEConstant::energyConvert;
    
    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
    
    /* 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();

      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){

         // mass = integrableObject->getMass();
        Vector3d vertexForce = randomForce/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();  

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