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
File size:	20373 byte(s)
Log Message:	Fixed segfault in because of directional atom in SMIPD.
#	Content
1	/*
2	* Copyright (c) 2008 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Acknowledgement of the program authors must be made in any
10	* publication of scientific results based in part on use of the
11	* program. An acceptable form of acknowledgement is citation of
12	* the article in which the program was described (Matthew
13	* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	* Parallel Simulation Engine for Molecular Dynamics,"
16	* J. Comput. Chem. 26, pp. 252-271 (2005))
17	*
18	* 2. Redistributions of source code must retain the above copyright
19	* notice, this list of conditions and the following disclaimer.
20	*
21	* 3. Redistributions in binary form must reproduce the above copyright
22	* notice, this list of conditions and the following disclaimer in the
23	* documentation and/or other materials provided with the
24	* distribution.
25	*
26	* This software is provided "AS IS," without a warranty of any
27	* kind. All express or implied conditions, representations and
28	* warranties, including any implied warranty of merchantability,
29	* fitness for a particular purpose or non-infringement, are hereby
30	* excluded. The University of Notre Dame and its licensors shall not
31	* be liable for any damages suffered by licensee as a result of
32	* using, modifying or distributing the software or its
33	* derivatives. In no event will the University of Notre Dame or its
34	* licensors be liable for any lost revenue, profit or data, or for
35	* direct, indirect, special, consequential, incidental or punitive
36	* damages, however caused and regardless of the theory of liability,
37	* arising out of the use of or inability to use software, even if the
38	* University of Notre Dame has been advised of the possibility of
39	* such damages.
40	*/
41	#include <fstream>
42	#include <iostream>
43	#include "integrators/SMIPDForceManager.hpp"
44	#include "math/CholeskyDecomposition.hpp"
45	#include "utils/OOPSEConstant.hpp"
46	#include "hydrodynamics/Sphere.hpp"
47	#include "hydrodynamics/Ellipsoid.hpp"
48	#include "utils/ElementsTable.hpp"
49	#include "math/ConvexHull.hpp"
50	#include "math/Triangle.hpp"
51
52
53	namespace oopse {
54
55	SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
56	simParams = info->getSimParams();
57	veloMunge = new Velocitizer(info);
58
59	// Create Hull, Convex Hull for now, other options later.
60	surfaceMesh_ = new ConvexHull();
61
62
63
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	targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
84	}
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
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	// 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	"HydroPropFile must be set to a file name if SMIPDynamics\n"
140	"\tis specified for rigidBodies which contain more than one atom\n"
141	"\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	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	currShape = new Sphere(atom->getPos(), 2.0);
214	} 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	currShape = new Sphere(atom->getPos(), 2.0);
226	} 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	HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
237	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	//Build a vector of integrable objects to determine if the are surface atoms
252	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
260	}
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
296	/Compute surface Mesh/
297	surfaceMesh_->computeHull(localSites_);
298
299	/* Get area and number of surface stunt doubles and compute new variance */
300	RealType area = surfaceMesh_->getArea();
301	int nSurfaceSDs = surfaceMesh_->getNs();
302
303
304	std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
305	int nTriangles = sMesh.size();
306
307
308
309	/* Compute variance for random forces */
310
311	RealType sigma_t = sqrt(NumericConstant::PI/2.0)((targetPressure_)area/nTriangles)
312	/OOPSEConstant::energyConvert;
313
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	/* Loop over the mesh faces and apply random force to each of the faces*/
320
321
322	std::vector<Triangle>::iterator face;
323	std::vector<StuntDouble*>::iterator vertex;
324	int thisNumber = 0;
325	for (face = sMesh.begin(); face != sMesh.end(); ++face){
326
327	Triangle thisTriangle = *face;
328	std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
329
330	/* Get Random Force */
331	Vector3d unitNormal = thisTriangle.getNormal();
332	unitNormal.normalize();
333	Vector3d randomForce = -randNums[thisNumber] * unitNormal;
334	Vector3d centroid = thisTriangle.getCentroid();
335
336	Vector3d langevinForce = randomForce - gamma_t_*thisTriangle.getFacetVelocity();
337
338	for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
339	if ((*vertex) != NULL){
340	// mass = integrableObject->getMass();
341	Vector3d vertexForce = langevinForce/3.0;
342	(*vertex)->addFrc(vertexForce);
343
344	if ((*vertex)->isDirectional()){
345	Vector3d vertexPos = (*vertex)->getPos();
346	Vector3d vertexCentroidVector = vertexPos - centroid;
347	(*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
348	}
349	}
350	}
351	}
352
353	/* Now loop over all surface particles and apply the drag*/
354	/*
355	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	//apply random force and torque at center of resistance
367	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
400	//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
411	} 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
432	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
433	angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
434
435	// check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
436
437	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
447
448	} else {
449	//spherical atom
450
451	// 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	//frictionForce = -gamma_t*velStep;
476	// re-estimate velocities at full-step using friction forces:
477
478	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
479
480	// check for convergence (if the vector has converged, fdot will be 1.0):
481
482	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
493
494	}
495	*/
496	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
497	currSnapshot->setVolume(surfaceMesh_->getVolume());
498
499	ForceManager::postCalculation(needStress);
500	}
501
502	void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
503
504
505	Vector<RealType, 6> Z;
506	Vector<RealType, 6> generalForce;
507
508
509	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
562
563
564
565	}