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.
#	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	/* 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	"HydroPropFile must be set to a file name if SMIPDynamics\n"
125	"\tis specified for rigidBodies which contain more than one atom\n"
126	"\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	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	HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
222	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	//variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
247
248	}
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
284	/Compute surface Mesh/
285	surfaceMesh_->computeHull(localSites_);
286
287	/* Get area and number of surface stunt doubles and compute new variance */
288	RealType area = surfaceMesh_->getArea();
289	int nSurfaceSDs = surfaceMesh_->getNs();
290
291	/* Compute variance for random forces */
292
293	variance_ = sqrt(2.0NumericConstant::PI)((targetPressure_/OOPSEConstant::pressureConvert)*area/nSurfaceSDs)
294	/OOPSEConstant::energyConvert;
295
296	std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
297	std::vector<RealType> randNums = genTriangleForces(sMesh.size(),variance_);
298
299	/* Loop over the mesh faces and apply random force to each of the faces*/
300
301
302	std::vector<Triangle*>::iterator face;
303	std::vector<StuntDouble*>::iterator vertex;
304	int thisNumber = 0;
305	for (face = sMesh.begin(); face != sMesh.end(); ++face){
306
307	Triangle* thisTriangle = *face;
308	std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
309
310	/* Get Random Force */
311	Vector3d unitNormal = thisTriangle->getNormal();
312	unitNormal.normalize();
313	Vector3d randomForce = -randNums[thisNumber] * unitNormal;
314	Vector3d centroid = thisTriangle->getCentroid();
315
316	for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
317
318	// mass = integrableObject->getMass();
319	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	}
328	}
329
330	/* Now loop over all surface particles and apply the drag*/
331
332	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
344
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
378	//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
389	} 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
410	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
411	angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
412
413	// check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
414
415	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
425
426	} else {
427	//spherical atom
428
429
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
457	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
458
459	// check for convergence (if the vector has converged, fdot will be 1.0):
460
461	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
477	ForceManager::postCalculation(needStress);
478	}
479
480	void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
481
482
483	Vector<RealType, 6> Z;
484	Vector<RealType, 6> generalForce;
485
486
487	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	}