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(simParams->getViscosity(),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 */
    
     RealType TD_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(),TD_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();  
        //apply random force and torque at center of resistance

        Vector3d randomForceBody;
        Vector3d randomTorqueBody;
        genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
        Vector3d randomForceLab = Atrans * randomForceBody;
        Vector3d randomTorqueLab = Atrans * randomTorqueBody;
        integrableObject->addFrc(randomForceLab);            
        integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));             
        

        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

        Vector3d randomForce;
        Vector3d randomTorque;
        genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
        integrableObject->addFrc(randomForce);   
         
        // 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:	1306
Committed:	Thu Oct 16 18:25:36 2008 UTC (16 years, 6 months ago) by chuckv
File size:	20257 byte(s)
Log Message:	Returned old LD forces to SMIPDForceManager.
#	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	1306	HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),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	1306	variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
247	chuckv	1301
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	1306	RealType TD_variance = sqrt(2.0NumericConstant::PI)((targetPressure_/OOPSEConstant::pressureConvert)*area/nSurfaceSDs)
294	chuckv	1304	/OOPSEConstant::energyConvert;
295	chuckv	1301
296	chuckv	1302	std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
297	chuckv	1306	std::vector<RealType> randNums = genTriangleForces(sMesh.size(),TD_variance);
298	chuckv	1302
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	1306	//apply random force and torque at center of resistance
344	chuckv	1293
345	chuckv	1306	Vector3d randomForceBody;
346			Vector3d randomTorqueBody;
347			genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
348			Vector3d randomForceLab = Atrans * randomForceBody;
349			Vector3d randomTorqueLab = Atrans * randomTorqueBody;
350			integrableObject->addFrc(randomForceLab);
351			integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));
352
353
354	chuckv	1302
355			Mat3x3d I = integrableObject->getI();
356			Vector3d omegaBody;
357
358			// What remains contains velocity explicitly, but the velocity required
359			// is at the full step: v(t + h), while we have initially the velocity
360			// at the half step: v(t + h/2). We need to iterate to converge the
361			// friction force and friction torque vectors.
362
363			// this is the velocity at the half-step:
364
365			Vector3d vel =integrableObject->getVel();
366			Vector3d angMom = integrableObject->getJ();
367
368			//estimate velocity at full-step using everything but friction forces:
369
370			frc = integrableObject->getFrc();
371			Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
372
373			Tb = integrableObject->lab2Body(integrableObject->getTrq());
374			Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;
375
376			Vector3d omegaLab;
377			Vector3d vcdLab;
378			Vector3d vcdBody;
379			Vector3d frictionForceBody;
380			Vector3d frictionForceLab(0.0);
381			Vector3d oldFFL; // used to test for convergence
382			Vector3d frictionTorqueBody(0.0);
383			Vector3d oldFTB; // used to test for convergence
384			Vector3d frictionTorqueLab;
385			RealType fdot;
386			RealType tdot;
387	chuckv	1293
388	chuckv	1302	//iteration starts here:
389
390			for (int k = 0; k < maxIterNum_; k++) {
391
392			if (integrableObject->isLinear()) {
393			int linearAxis = integrableObject->linearAxis();
394			int l = (linearAxis +1 )%3;
395			int m = (linearAxis +2 )%3;
396			omegaBody[l] = angMomStep[l] /I(l, l);
397			omegaBody[m] = angMomStep[m] /I(m, m);
398	chuckv	1293
399	chuckv	1302	} else {
400			omegaBody[0] = angMomStep[0] /I(0, 0);
401			omegaBody[1] = angMomStep[1] /I(1, 1);
402			omegaBody[2] = angMomStep[2] /I(2, 2);
403			}
404
405			omegaLab = Atrans * omegaBody;
406
407			// apply friction force and torque at center of resistance
408
409			vcdLab = velStep + cross(omegaLab, rcrLab);
410			vcdBody = A * vcdLab;
411			frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
412			oldFFL = frictionForceLab;
413			frictionForceLab = Atrans * frictionForceBody;
414			oldFTB = frictionTorqueBody;
415			frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
416			frictionTorqueLab = Atrans * frictionTorqueBody;
417
418			// re-estimate velocities at full-step using friction forces:
419	chuckv	1293
420	chuckv	1302	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
421			angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
422	chuckv	1293
423	chuckv	1302	// check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
424	chuckv	1293
425	chuckv	1302	fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
426			tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
427
428			if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
429			break; // iteration ends here
430			}
431
432			integrableObject->addFrc(frictionForceLab);
433			integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
434	chuckv	1293
435
436	chuckv	1302	} else {
437			//spherical atom
438	chuckv	1293
439	chuckv	1306	Vector3d randomForce;
440			Vector3d randomTorque;
441			genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
442			integrableObject->addFrc(randomForce);
443
444	chuckv	1302	// What remains contains velocity explicitly, but the velocity required
445			// is at the full step: v(t + h), while we have initially the velocity
446			// at the half step: v(t + h/2). We need to iterate to converge the
447			// friction force vector.
448
449			// this is the velocity at the half-step:
450
451			Vector3d vel =integrableObject->getVel();
452
453			//estimate velocity at full-step using everything but friction forces:
454
455			frc = integrableObject->getFrc();
456			Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
457
458			Vector3d frictionForce(0.0);
459			Vector3d oldFF; // used to test for convergence
460			RealType fdot;
461
462			//iteration starts here:
463
464			for (int k = 0; k < maxIterNum_; k++) {
465
466			oldFF = frictionForce;
467			frictionForce = -hydroProps_[index]->getXitt() * velStep;
468
469			// re-estimate velocities at full-step using friction forces:
470	chuckv	1293
471	chuckv	1302	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
472
473			// check for convergence (if the vector has converged, fdot will be 1.0):
474	chuckv	1293
475	chuckv	1302	fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
476
477			if (fabs(1.0 - fdot) <= forceTolerance_)
478			break; // iteration ends here
479			}
480
481			integrableObject->addFrc(frictionForce);
482
483
484			}
485
486
487			}
488			Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
489			currSnapshot->setVolume(surfaceMesh_->getVolume());
490	chuckv	1293
491	chuckv	1302	ForceManager::postCalculation(needStress);
492	chuckv	1293	}
493
494	chuckv	1302	void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
495	chuckv	1293
496	chuckv	1302
497	chuckv	1293	Vector<RealType, 6> Z;
498			Vector<RealType, 6> generalForce;
499	chuckv	1302
500
501	chuckv	1293	Z[0] = randNumGen_.randNorm(0, variance);
502			Z[1] = randNumGen_.randNorm(0, variance);
503			Z[2] = randNumGen_.randNorm(0, variance);
504			Z[3] = randNumGen_.randNorm(0, variance);
505			Z[4] = randNumGen_.randNorm(0, variance);
506			Z[5] = randNumGen_.randNorm(0, variance);
507
508			generalForce = hydroProps_[index]->getS()*Z;
509
510			force[0] = generalForce[0];
511			force[1] = generalForce[1];
512			force[2] = generalForce[2];
513			torque[0] = generalForce[3];
514			torque[1] = generalForce[4];
515			torque[2] = generalForce[5];
516
517			}
518			std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
519
520			// zero fill the random vector before starting:
521			std::vector<RealType> gaussRand;
522			gaussRand.resize(nTriangles);
523			std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
524
525
526			#ifdef IS_MPI
527			if (worldRank == 0) {
528			#endif
529			for (int i = 0; i < nTriangles; i++) {
530			gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));
531			}
532			#ifdef IS_MPI
533			}
534			#endif
535
536			// push these out to the other processors
537
538			#ifdef IS_MPI
539			if (worldRank == 0) {
540			MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
541			} else {
542			MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
543			}
544			#endif
545
546			for (int i = 0; i < nTriangles; i++) {
547			gaussRand[i] = gaussRand[i] * variance;
548			}
549
550			return gaussRand;
551			}
552
553	chuckv	1306
554
555
556
557	chuckv	1293	}