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, 
              "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "   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 {
      // Convert pressure from atm -> amu/(fs^2*Ang)
      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();
    } 
    
    if (!simParams->haveViscosity()) {
      sprintf(painCave.errMsg, 
              "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "   without a viscosity!\n");
      painCave.isFatal = 1;
      painCave.severity = OOPSE_ERROR;
      simError();
    }else{
      viscosity_ = simParams->getViscosity();
    }
    
    dt_ = simParams->getDt();
  
    variance_ = 2.0 * OOPSEConstant::kb * targetTemp_ / dt_;

    // 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;
        }
        
      }
    }

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

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

    // 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;
    int index = 0;
  
    // Compute surface Mesh
    surfaceMesh_->computeHull(localSites_);

    // Get total area and number of surface stunt doubles
    RealType area = surfaceMesh_->getArea();
    int nSurfaceSDs = surfaceMesh_->getNs();
    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
    int nTriangles = sMesh.size();

    // Loop over the mesh faces and apply external pressure 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();
      RealType thisArea = thisTriangle.getArea(); 
      Vector3d unitNormal = thisTriangle.getNormal();
      unitNormal.normalize();
      Vector3d centroid = thisTriangle.getCentroid();
      Vector3d extPressure = - unitNormal * (targetPressure_ * thisArea)
        / OOPSEConstant::energyConvert;
      
      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
        if ((*vertex) != NULL){ 
          Vector3d vertexForce = extPressure/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 
    // and random forces
    
    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
      integrableObject = *vertex;      
      index = integrableObject->getLocalIndex();
      mass = integrableObject->getMass();
      if (integrableObject->isDirectional()){
        
        // preliminaries for directional objects:
        
        A = integrableObject->getA();
        Atrans = A.transpose();
        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  

        // apply random force and torque at center of resistance

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

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

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

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

            
      } else {
        //spherical atom

        Vector3d randomForce;
        Vector3d randomTorque;
        genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
        integrableObject->addFrc(randomForce);            
        
        // What remains contains velocity explicitly, but the velocity required
        // is at the full step: v(t + h), while we have initially the velocity
        // at the half step: v(t + h/2).  We need to iterate to converge the
        // friction force vector.
        
        // this is the velocity at the half-step:
        
        Vector3d vel =integrableObject->getVel();
        
        //estimate velocity at full-step using everything but friction forces:           
        
        frc = integrableObject->getFrc();
        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
        
        Vector3d frictionForce(0.0);
        Vector3d oldFF;  // used to test for convergence
        RealType fdot;
        
        //iteration starts here:
        
        for (int k = 0; k < maxIterNum_; k++) {
          
          oldFF = frictionForce;                            
          frictionForce = -hydroProps_[index]->getXitt() * velStep;
          //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);
                
      }       
    }
        
    veloMunge->removeComDrift();
    veloMunge->removeAngularDrift();

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