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 ((*vertex)->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:	1312
Committed:	Wed Oct 22 17:52:40 2008 UTC (16 years, 6 months ago) by chuckv
Original Path:	*trunk/src/integrators/SMIPDForceManager.cpp*
File size:	20373 byte(s)
Log Message:	Fixed segfault in because of directional atom in SMIPD.
#	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	chuckv	1308	std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
305	chuckv	1307	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	1308	std::vector<Triangle>::iterator face;
323	chuckv	1302	std::vector<StuntDouble*>::iterator vertex;
324			int thisNumber = 0;
325	chuckv	1293	for (face = sMesh.begin(); face != sMesh.end(); ++face){
326
327	chuckv	1308	Triangle thisTriangle = *face;
328			std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
329	chuckv	1302
330			/* Get Random Force */
331	chuckv	1308	Vector3d unitNormal = thisTriangle.getNormal();
332	chuckv	1302	unitNormal.normalize();
333			Vector3d randomForce = -randNums[thisNumber] * unitNormal;
334	chuckv	1308	Vector3d centroid = thisTriangle.getCentroid();
335	chuckv	1304
336	chuckv	1308	Vector3d langevinForce = randomForce - gamma_t_*thisTriangle.getFacetVelocity();
337	chuckv	1307
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	chuckv	1312
344			if ((*vertex)->isDirectional()){
345	chuckv	1307	Vector3d vertexPos = (*vertex)->getPos();
346			Vector3d vertexCentroidVector = vertexPos - centroid;
347			(*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
348			}
349			}
350	chuckv	1293	}
351	chuckv	1302	}
352	chuckv	1293
353	chuckv	1302	/* Now loop over all surface particles and apply the drag*/
354	chuckv	1307	/*
355	chuckv	1302	std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
356			for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
357			integrableObject = *vertex;
358			mass = integrableObject->getMass();
359			if (integrableObject->isDirectional()){
360
361			// preliminaries for directional objects:
362
363			A = integrableObject->getA();
364			Atrans = A.transpose();
365			Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();
366	chuckv	1306	//apply random force and torque at center of resistance
367	chuckv	1302	Mat3x3d I = integrableObject->getI();
368			Vector3d omegaBody;
369
370			// What remains contains velocity explicitly, but the velocity required
371			// is at the full step: v(t + h), while we have initially the velocity
372			// at the half step: v(t + h/2). We need to iterate to converge the
373			// friction force and friction torque vectors.
374
375			// this is the velocity at the half-step:
376
377			Vector3d vel =integrableObject->getVel();
378			Vector3d angMom = integrableObject->getJ();
379
380			//estimate velocity at full-step using everything but friction forces:
381
382			frc = integrableObject->getFrc();
383			Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
384
385			Tb = integrableObject->lab2Body(integrableObject->getTrq());
386			Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;
387
388			Vector3d omegaLab;
389			Vector3d vcdLab;
390			Vector3d vcdBody;
391			Vector3d frictionForceBody;
392			Vector3d frictionForceLab(0.0);
393			Vector3d oldFFL; // used to test for convergence
394			Vector3d frictionTorqueBody(0.0);
395			Vector3d oldFTB; // used to test for convergence
396			Vector3d frictionTorqueLab;
397			RealType fdot;
398			RealType tdot;
399	chuckv	1293
400	chuckv	1302	//iteration starts here:
401
402			for (int k = 0; k < maxIterNum_; k++) {
403
404			if (integrableObject->isLinear()) {
405			int linearAxis = integrableObject->linearAxis();
406			int l = (linearAxis +1 )%3;
407			int m = (linearAxis +2 )%3;
408			omegaBody[l] = angMomStep[l] /I(l, l);
409			omegaBody[m] = angMomStep[m] /I(m, m);
410	chuckv	1293
411	chuckv	1302	} else {
412			omegaBody[0] = angMomStep[0] /I(0, 0);
413			omegaBody[1] = angMomStep[1] /I(1, 1);
414			omegaBody[2] = angMomStep[2] /I(2, 2);
415			}
416
417			omegaLab = Atrans * omegaBody;
418
419			// apply friction force and torque at center of resistance
420
421			vcdLab = velStep + cross(omegaLab, rcrLab);
422			vcdBody = A * vcdLab;
423			frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
424			oldFFL = frictionForceLab;
425			frictionForceLab = Atrans * frictionForceBody;
426			oldFTB = frictionTorqueBody;
427			frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
428			frictionTorqueLab = Atrans * frictionTorqueBody;
429
430			// re-estimate velocities at full-step using friction forces:
431	chuckv	1293
432	chuckv	1302	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
433			angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
434	chuckv	1293
435	chuckv	1302	// check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
436	chuckv	1293
437	chuckv	1302	fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
438			tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
439
440			if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
441			break; // iteration ends here
442			}
443
444			integrableObject->addFrc(frictionForceLab);
445			integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
446	chuckv	1293
447
448	chuckv	1302	} else {
449	chuckv	1307	//spherical atom
450	chuckv	1293
451	chuckv	1302	// What remains contains velocity explicitly, but the velocity required
452			// is at the full step: v(t + h), while we have initially the velocity
453			// at the half step: v(t + h/2). We need to iterate to converge the
454			// friction force vector.
455
456			// this is the velocity at the half-step:
457
458			Vector3d vel =integrableObject->getVel();
459
460			//estimate velocity at full-step using everything but friction forces:
461
462			frc = integrableObject->getFrc();
463			Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
464
465			Vector3d frictionForce(0.0);
466			Vector3d oldFF; // used to test for convergence
467			RealType fdot;
468
469			//iteration starts here:
470
471			for (int k = 0; k < maxIterNum_; k++) {
472
473			oldFF = frictionForce;
474			frictionForce = -hydroProps_[index]->getXitt() * velStep;
475	chuckv	1307	//frictionForce = -gamma_t*velStep;
476	chuckv	1302	// re-estimate velocities at full-step using friction forces:
477	chuckv	1293
478	chuckv	1302	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
479
480			// check for convergence (if the vector has converged, fdot will be 1.0):
481	chuckv	1293
482	chuckv	1302	fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
483
484			if (fabs(1.0 - fdot) <= forceTolerance_)
485			break; // iteration ends here
486			}
487
488			integrableObject->addFrc(frictionForce);
489
490
491			}
492	chuckv	1307
493	chuckv	1302
494	chuckv	1307	}
495			*/
496	chuckv	1302	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
497			currSnapshot->setVolume(surfaceMesh_->getVolume());
498	chuckv	1293
499	chuckv	1302	ForceManager::postCalculation(needStress);
500	chuckv	1293	}
501
502	chuckv	1302	void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
503	chuckv	1293
504	chuckv	1302
505	chuckv	1293	Vector<RealType, 6> Z;
506			Vector<RealType, 6> generalForce;
507	chuckv	1302
508
509	chuckv	1293	Z[0] = randNumGen_.randNorm(0, variance);
510			Z[1] = randNumGen_.randNorm(0, variance);
511			Z[2] = randNumGen_.randNorm(0, variance);
512			Z[3] = randNumGen_.randNorm(0, variance);
513			Z[4] = randNumGen_.randNorm(0, variance);
514			Z[5] = randNumGen_.randNorm(0, variance);
515
516			generalForce = hydroProps_[index]->getS()*Z;
517
518			force[0] = generalForce[0];
519			force[1] = generalForce[1];
520			force[2] = generalForce[2];
521			torque[0] = generalForce[3];
522			torque[1] = generalForce[4];
523			torque[2] = generalForce[5];
524
525			}
526			std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
527
528			// zero fill the random vector before starting:
529			std::vector<RealType> gaussRand;
530			gaussRand.resize(nTriangles);
531			std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
532
533
534			#ifdef IS_MPI
535			if (worldRank == 0) {
536			#endif
537			for (int i = 0; i < nTriangles; i++) {
538			gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));
539			}
540			#ifdef IS_MPI
541			}
542			#endif
543
544			// push these out to the other processors
545
546			#ifdef IS_MPI
547			if (worldRank == 0) {
548			MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
549			} else {
550			MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
551			}
552			#endif
553
554			for (int i = 0; i < nTriangles; i++) {
555			gaussRand[i] = gaussRand[i] * variance;
556			}
557
558			return gaussRand;
559			}
560
561	chuckv	1306
562
563
564
565	chuckv	1293	}