src/integrators/SMIPDForceManager.cpp

/*
 * Copyright (c) 2008 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
#include <fstream> 
#include <iostream>
#include "integrators/SMIPDForceManager.hpp"
#include "math/CholeskyDecomposition.hpp"
#include "utils/OOPSEConstant.hpp"
#include "hydrodynamics/Sphere.hpp"
#include "hydrodynamics/Ellipsoid.hpp"
#include "utils/ElementsTable.hpp"
#include "math/ConvexHull.hpp"
#include "math/Triangle.hpp"


namespace oopse {

  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
    simParams = info->getSimParams();
    veloMunge = new Velocitizer(info);
    
    // Create Hull, Convex Hull for now, other options later.
    surfaceMesh_ = new ConvexHull();

    
    /* Check that the simulation has target pressure and target
       temperature set*/

    if (!simParams->haveTargetTemp()) {
      sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
      painCave.isFatal = 1;
      painCave.severity = OOPSE_ERROR;
      simError();
    } else {
      targetTemp_ = simParams->getTargetTemp();
    }

    if (!simParams->haveTargetPressure()) {
      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "   without a targetPressure!\n");
      
      painCave.isFatal = 1;
      simError();
    } else {
      targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
    }

   
    if (simParams->getUsePeriodicBoundaryConditions()) {
      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "   with periodic boundary conditions !\n");
      
      painCave.isFatal = 1;
      simError();
    } 


    //Compute initial hull
    /*
    surfaceMesh_->computeHull(localSites_);
    Area0_ = surfaceMesh_->getArea();
    int nTriangles = surfaceMesh_->getNMeshElements();
    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
    //RealType eta0 = gamma_0/
    */

    // Build the hydroProp map:
    std::map<std::string, HydroProp*> hydroPropMap;

    Molecule* mol;
    StuntDouble* integrableObject;
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;              
    bool needHydroPropFile = false;
    
    for (mol = info->beginMolecule(i); mol != NULL; 
         mol = info->nextMolecule(i)) {
      for (integrableObject = mol->beginIntegrableObject(j); 
           integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(j)) {
        
        if (integrableObject->isRigidBody()) {
          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
          if (rb->getNumAtoms() > 1) needHydroPropFile = true;
        }
        
      }
    }
        

    if (needHydroPropFile) {               
      if (simParams->haveHydroPropFile()) {
        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
      } else {               
        sprintf( painCave.errMsg,
                 "HydroPropFile must be set to a file name if SMIPDynamics\n"
                 "\tis specified for rigidBodies which contain more than one atom\n"
                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n"
                 "\tFor use with SMIPD, the default viscosity in Hydro should be\n"
                 "\tset to 1.0 because the friction and random forces will be\n"
                 "\tdynamically re-set assuming this is true.\n");
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();  
      }      

      for (mol = info->beginMolecule(i); mol != NULL; 
           mol = info->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); 
             integrableObject != NULL;
             integrableObject = mol->nextIntegrableObject(j)) {

          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
          if (iter != hydroPropMap.end()) {
            hydroProps_.push_back(iter->second);
          } else {
            sprintf( painCave.errMsg,
                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
            painCave.severity = OOPSE_ERROR;
            painCave.isFatal = 1;
            simError();  
          }        
        }
      }
    } else {
      
      std::map<std::string, HydroProp*> hydroPropMap;
      for (mol = info->beginMolecule(i); mol != NULL; 
           mol = info->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); 
             integrableObject != NULL;
             integrableObject = mol->nextIntegrableObject(j)) {
          Shape* currShape = NULL;

          if (integrableObject->isAtom()){
            Atom* atom = static_cast<Atom*>(integrableObject);
            AtomType* atomType = atom->getAtomType();
            if (atomType->isGayBerne()) {
              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);              
              GenericData* data = dAtomType->getPropertyByName("GayBerne");
              if (data != NULL) {
                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
                
                if (gayBerneData != NULL) {  
                  GayBerneParam gayBerneParam = gayBerneData->getData();
                  currShape = new Ellipsoid(V3Zero, 
                                            gayBerneParam.GB_l / 2.0, 
                                            gayBerneParam.GB_d / 2.0, 
                                            Mat3x3d::identity());
                } else {
                  sprintf( painCave.errMsg,
                           "Can not cast GenericData to GayBerneParam\n");
                  painCave.severity = OOPSE_ERROR;
                  painCave.isFatal = 1;
                  simError();   
                }
              } else {
                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
                painCave.severity = OOPSE_ERROR;
                painCave.isFatal = 1;
                simError();    
              }
            } else {
              if (atomType->isLennardJones()){
                GenericData* data = atomType->getPropertyByName("LennardJones");
                if (data != NULL) {
                  LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
                  if (ljData != NULL) {
                    LJParam ljParam = ljData->getData();
                    currShape = new Sphere(atom->getPos(), 2.0);
                  } else {
                    sprintf( painCave.errMsg,
                             "Can not cast GenericData to LJParam\n");
                    painCave.severity = OOPSE_ERROR;
                    painCave.isFatal = 1;
                    simError();          
                  }       
                }
              } else {
                int aNum = etab.GetAtomicNum((atom->getType()).c_str());
                if (aNum != 0) {
                  currShape = new Sphere(atom->getPos(), 2.0);
                } else {
                  sprintf( painCave.errMsg,
                           "Could not find atom type in default element.txt\n");
                  painCave.severity = OOPSE_ERROR;
                  painCave.isFatal = 1;
                  simError();          
                }
              }
            }
          }
          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
          if (iter != hydroPropMap.end()) 
            hydroProps_.push_back(iter->second);
          else {
            currHydroProp->complete();
            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
            hydroProps_.push_back(currHydroProp);
          }
        }
      }
    }

    /* Compute hull first time through to get the area of t=0*/

    //Build a vector of integrable objects to determine if the are surface atoms 
    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(j)) {   
        localSites_.push_back(integrableObject);
      }
    }


  }  

  std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
    std::map<std::string, HydroProp*> props;
    std::ifstream ifs(filename.c_str());
    if (ifs.is_open()) {
      
    }
    
    const unsigned int BufferSize = 65535;
    char buffer[BufferSize];   
    while (ifs.getline(buffer, BufferSize)) {
      HydroProp* currProp = new HydroProp(buffer);
      props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
    }

    return props;
  }
   
  void SMIPDForceManager::postCalculation(bool needStress){
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    RealType mass;
    Vector3d pos;
    Vector3d frc;
    Mat3x3d A;
    Mat3x3d Atrans;
    Vector3d Tb;
    Vector3d ji;
    unsigned int index = 0;
    int fdf;
   
    fdf = 0;
  
    /*Compute surface Mesh*/
    surfaceMesh_->computeHull(localSites_);

    /* Get area and number of surface stunt doubles and compute new variance */
     RealType area = surfaceMesh_->getArea();
     int nSurfaceSDs = surfaceMesh_->getNs();

    
    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
    int nTriangles = sMesh.size();


     /* Compute variance for random forces */
   
    RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*area/nTriangles)
       /OOPSEConstant::energyConvert;

    gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * area * area * simParams->getDt()) /
      (4.0 * nTriangles * nTriangles* OOPSEConstant::kB*simParams->getTargetTemp()) /OOPSEConstant::energyConvert;

    std::vector<RealType>  randNums = genTriangleForces(nTriangles, sigma_t);

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

      Vector3d langevinForce = randomForce - gamma_t_*thisTriangle->getFacetVelocity();
      
      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
        if ((*vertex) != NULL){
          // mass = integrableObject->getMass();
          Vector3d vertexForce = langevinForce/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
        Mat3x3d I = integrableObject->getI();
        Vector3d omegaBody;
        
        // What remains contains velocity explicitly, but the velocity required
        // is at the full step: v(t + h), while we have initially the velocity
        // at the half step: v(t + h/2).  We need to iterate to converge the
        // friction force and friction torque vectors.
        
        // this is the velocity at the half-step:
        
        Vector3d vel =integrableObject->getVel();
        Vector3d angMom = integrableObject->getJ();
        
        //estimate velocity at full-step using everything but friction forces:           
        
        frc = integrableObject->getFrc();
        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
        
        Tb = integrableObject->lab2Body(integrableObject->getTrq());
        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                             
        
        Vector3d omegaLab;
        Vector3d vcdLab;
        Vector3d vcdBody;
        Vector3d frictionForceBody;
        Vector3d frictionForceLab(0.0);
        Vector3d oldFFL;  // used to test for convergence
        Vector3d frictionTorqueBody(0.0);
        Vector3d oldFTB;  // used to test for convergence
        Vector3d frictionTorqueLab;
        RealType fdot;
        RealType tdot;

        //iteration starts here:
        
        for (int k = 0; k < maxIterNum_; k++) {
          
          if (integrableObject->isLinear()) {
            int linearAxis = integrableObject->linearAxis();
            int l = (linearAxis +1 )%3;
            int m = (linearAxis +2 )%3;
            omegaBody[l] = angMomStep[l] /I(l, l);
            omegaBody[m] = angMomStep[m] /I(m, m);
            
          } else {
            omegaBody[0] = angMomStep[0] /I(0, 0);
            omegaBody[1] = angMomStep[1] /I(1, 1);
            omegaBody[2] = angMomStep[2] /I(2, 2);
          }
          
          omegaLab = Atrans * omegaBody;
          
          // apply friction force and torque at center of resistance
          
          vcdLab = velStep + cross(omegaLab, rcrLab);       
          vcdBody = A * vcdLab;
          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
          oldFFL = frictionForceLab;
          frictionForceLab = Atrans * frictionForceBody;
          oldFTB = frictionTorqueBody;
          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
          frictionTorqueLab = Atrans * frictionTorqueBody;
          
          // re-estimate velocities at full-step using friction forces:
              
          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);

          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
              
          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
          
          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
            break; // iteration ends here
        }
        
        integrableObject->addFrc(frictionForceLab);
        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));

            
      } else {
        //spherical atom

        // What remains contains velocity explicitly, but the velocity required
        // is at the full step: v(t + h), while we have initially the velocity
        // at the half step: v(t + h/2).  We need to iterate to converge the
        // friction force vector.
        
        // this is the velocity at the half-step:
        
        Vector3d vel =integrableObject->getVel();
        
        //estimate velocity at full-step using everything but friction forces:           
        
        frc = integrableObject->getFrc();
        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
        
        Vector3d frictionForce(0.0);
        Vector3d oldFF;  // used to test for convergence
        RealType fdot;
        
        //iteration starts here:
        
        for (int k = 0; k < maxIterNum_; k++) {
          
          oldFF = frictionForce;                            
          frictionForce = -hydroProps_[index]->getXitt() * velStep;
          //frictionForce = -gamma_t*velStep;
          // re-estimate velocities at full-step using friction forces:
          
          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
          
          // check for convergence (if the vector has converged, fdot will be 1.0):
          
          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
          
          if (fabs(1.0 - fdot) <= forceTolerance_)
            break; // iteration ends here
        }
        
        integrableObject->addFrc(frictionForce);
        
        
      }
  
      
  }
    */
    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
    currSnapshot->setVolume(surfaceMesh_->getVolume());
    
    ForceManager::postCalculation(needStress);   
  }

  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {

    
    Vector<RealType, 6> Z;
    Vector<RealType, 6> generalForce;
    

    Z[0] = randNumGen_.randNorm(0, variance);
    Z[1] = randNumGen_.randNorm(0, variance);
    Z[2] = randNumGen_.randNorm(0, variance);
    Z[3] = randNumGen_.randNorm(0, variance);
    Z[4] = randNumGen_.randNorm(0, variance);
    Z[5] = randNumGen_.randNorm(0, variance);
     
    generalForce = hydroProps_[index]->getS()*Z;
    
    force[0] = generalForce[0];
    force[1] = generalForce[1];
    force[2] = generalForce[2];
    torque[0] = generalForce[3];
    torque[1] = generalForce[4];
    torque[2] = generalForce[5];
    
} 
  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {

    // zero fill the random vector before starting:
    std::vector<RealType> gaussRand;
    gaussRand.resize(nTriangles);
    std::fill(gaussRand.begin(), gaussRand.end(), 0.0);


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