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 Langevin Dynamics\n"
                 "\tis specified for rigidBodies which contain more than one atom\n"
                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\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();
     //RealType nSurfaceSDs = surfaceMesh_->getNs();

//    std::cerr << "Surface Area is: " << area << " nSurfaceSDs is: " << nSurfaceSDs << std::endl;

    /* Compute variance for random forces */
    
//    variance_ = sqrt(2.0*NumericConstant::PI)*(targetPressure_*area/nSurfaceSDs);
    
//    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;
/*    
    for (face = sMesh.begin(); face != sMesh.end(); ++face){
     
      Triangle* thisTriangle = *face;
      std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();

      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
            std::cout << (*vertex)->getPos() << std::endl;
         // mass = integrableObject->getMass();

         //      integrableObject->addFrc(randomForce);           
      }


    }
  */ 
    /*
    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {

      
      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(j)) {
          
          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);

          
          }
          
        ++index;
    
      }
    }    
    */
    // info_->setFdf(fdf);
    // veloMunge->removeComDrift();
    // Remove angular drift if we are not using periodic boundary conditions.
    //if(!simParams->getUsePeriodicBoundaryConditions()) 
    //  veloMunge->removeAngularDrift();

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