src/integrators/SMLDForceManager.cpp

/*
 * Copyright (c) 2005 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.
 */

/*
Provides the force manager for Isothermal-Isobaric Langevin Dynamics where the stochastic force
is applied to the surface atoms anisotropically so as to provide a constant pressure. The 
surface atoms are determined by computing the convex hull and then triangulating that hull. The force
applied to the facets of the triangulation and mapped back onto the surface atoms. 
See: Kohanoff et.al. CHEMPHYSCHEM,2005,6,1848-1852.
*/

#include <fstream> 
#include <iostream>
#include "integrators/SMLDForceManager.hpp"
#include "math/CholeskyDecomposition.hpp"
#include "utils/OOPSEConstant.hpp"
#include "hydrodynamics/Sphere.hpp"
#include "hydrodynamics/Ellipsoid.hpp"
#if defined(HAVE_QHULL) || defined(HAVE_CGAL)
#ifdef HAVE_QHULL
#include "math/ConvexHull.hpp"
#endif

#ifdef HAVE_CGAL
#include "math/AlphaShape.hpp"
#endif
#endif
#include "utils/ElementsTable.hpp"

namespace oopse {

  SMLDForceManager::SMLDForceManager(SimInfo* info) : ForceManager(info){
    simParams = info->getSimParams();
    veloMunge = new Velocitizer(info);

    sphericalBoundaryConditions_ = false;
    if (simParams->getUseSphericalBoundaryConditions()) {
      sphericalBoundaryConditions_ = true;
      if (simParams->haveLangevinBufferRadius()) {
        langevinBufferRadius_ = simParams->getLangevinBufferRadius();
      } else {
        sprintf( painCave.errMsg,
                 "langevinBufferRadius must be specified " 
                 "when useSphericalBoundaryConditions is turned on.\n");
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();  
      }
    
      if (simParams->haveFrozenBufferRadius()) {
        frozenBufferRadius_ = simParams->getFrozenBufferRadius();
      } else {
        sprintf( painCave.errMsg,
                 "frozenBufferRadius must be specified " 
                 "when useSphericalBoundaryConditions is turned on.\n");
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();  
      }

      if (frozenBufferRadius_ < langevinBufferRadius_) {
        sprintf( painCave.errMsg,
                 "frozenBufferRadius has been set smaller than the " 
                 "langevinBufferRadius.  This is probably an error.\n");
        painCave.severity = OOPSE_WARNING;
        painCave.isFatal = 0;
        simError();  
      }
    }

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

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

    if (needHydroPropFile) {               
      if (simParams->haveHydroPropFile()) {
        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
      } else {               
        sprintf( painCave.errMsg,
                 "HydroPropFile must be set to a file name if Langevin\n"
                 "\tDynamics is specified for rigidBodies which contain more\n"
                 "\tthan one atom.  To create a HydroPropFile, run \"Hydro\".\n");
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();  
      }      

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

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

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

  std::map<std::string, HydroProp*> SMLDForceManager::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 SMLDForceManager::postCalculation(bool needStress){
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    RealType mass;
    Vector3d vel;
    Vector3d pos;
    Vector3d frc;
    Mat3x3d A;
    Mat3x3d Atrans;
    Vector3d Tb;
    Vector3d ji;
    unsigned int index = 0;
    bool doLangevinForces;
    bool freezeMolecule;
    int fdf;


    fdf = 0;

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

      doLangevinForces = true;           
      freezeMolecule = false;

      if (sphericalBoundaryConditions_) {
        
        Vector3d molPos = mol->getCom();
        RealType molRad = molPos.length();

        doLangevinForces = false;
        
        if (molRad > langevinBufferRadius_) { 
          doLangevinForces = true;
          freezeMolecule = false;
        }
        if (molRad > frozenBufferRadius_) {
          doLangevinForces = false;
          freezeMolecule = true;
        }
      }
      
      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(j)) {
          
        if (freezeMolecule) 
          fdf += integrableObject->freeze();
        
        if (doLangevinForces) {  
          vel =integrableObject->getVel(); 
          mass = integrableObject->getMass();
          if (integrableObject->isDirectional()){
            //calculate angular velocity in lab frame
            Mat3x3d I = integrableObject->getI();
            Vector3d angMom = integrableObject->getJ();
            Vector3d omega;
            
            if (integrableObject->isLinear()) {
              int linearAxis = integrableObject->linearAxis();
              int l = (linearAxis +1 )%3;
              int m = (linearAxis +2 )%3;
              omega[l] = angMom[l] /I(l, l);
              omega[m] = angMom[m] /I(m, m);
              
            } else {
              omega[0] = angMom[0] /I(0, 0);
              omega[1] = angMom[1] /I(1, 1);
              omega[2] = angMom[2] /I(2, 2);
            }

            //std::cerr << "I = " << I(0,0) << "\t" << I(1,1) << "\t" << I(2,2) << "\n\n";

            //apply friction force and torque at center of resistance
            A = integrableObject->getA();
            Atrans = A.transpose();
            //std::cerr << "A = " << integrableObject->getA() << "\n";
            //std::cerr << "Atrans = " << A.transpose() << "\n\n";
            Vector3d rcr = Atrans * hydroProps_[index]->getCOR();  
            //std::cerr << "cor = " << hydroProps_[index]->getCOR() << "\n\n\n\n";
            //std::cerr << "rcr = " << rcr << "\n\n";
            Vector3d vcdLab = vel + cross(omega, rcr);
        
            //std::cerr << "velL = " << vel << "\n\n";
            //std::cerr << "vcdL = " << vcdLab << "\n\n";
            Vector3d vcdBody = A* vcdLab;
            //std::cerr << "vcdB = " << vcdBody << "\n\n";
            Vector3d frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omega);

            //std::cerr << "xitt = " << hydroProps_[index]->getXitt() << "\n\n";
            //std::cerr << "ffB = " << frictionForceBody << "\n\n";
            Vector3d frictionForceLab = Atrans*frictionForceBody;
            //std::cerr << "ffL = " << frictionForceLab << "\n\n";
            //std::cerr << "frc = " << integrableObject->getFrc() << "\n\n"; 
            integrableObject->addFrc(frictionForceLab);
            //std::cerr << "frc = " << integrableObject->getFrc() << "\n\n"; 
            //std::cerr << "ome = " << omega << "\n\n";
            Vector3d frictionTorqueBody = - (hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omega);
            //std::cerr << "ftB = " << frictionTorqueBody << "\n\n";
            Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
            //std::cerr << "ftL = " << frictionTorqueLab << "\n\n";
            //std::cerr << "ftL2 = " << frictionTorqueLab+cross(rcr,frictionForceLab) << "\n\n";
            //std::cerr << "trq = " << integrableObject->getTrq() << "\n\n"; 
            integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
            //std::cerr << "trq = " << integrableObject->getTrq() << "\n\n"; 

            //apply random force and torque at center of resistance
            Vector3d randomForceBody;
            Vector3d randomTorqueBody;
            genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
            //std::cerr << "rfB = " << randomForceBody << "\n\n";
            //std::cerr << "rtB = " << randomTorqueBody << "\n\n";
            Vector3d randomForceLab = Atrans*randomForceBody;
            Vector3d randomTorqueLab = Atrans* randomTorqueBody;
            integrableObject->addFrc(randomForceLab);            
            //std::cerr << "rfL = " << randomForceLab << "\n\n";
            //std::cerr << "rtL = " << randomTorqueLab << "\n\n";
            //std::cerr << "rtL2 = " << randomTorqueLab + cross(rcr, randomForceLab) << "\n\n";
            integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));             
            
          } else {
            //spherical atom
            Vector3d frictionForce = -(hydroProps_[index]->getXitt() * vel);
            //std::cerr << "xitt = " << hydroProps_[index]->getXitt() << "\n\n";
            Vector3d randomForce;
            Vector3d randomTorque;
            genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
            
            integrableObject->addFrc(frictionForce+randomForce);             
          }
        }
          
        ++index;
    
      }
    }    

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

    ForceManager::postCalculation(needStress);   
  }

void SMLDForceManager::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:	1210
Committed:	Wed Jan 23 03:45:33 2008 UTC (17 years, 3 months ago) by gezelter
File size:	17654 byte(s)
Log Message:	Removed older version of openbabel from our code. We now have a configure check to see if openbabel is installed and then we link to the stuff we need. Conversion to OOPSE's md format is handled by only one application (atom2md), so most of the work went on there. ElementsTable still needs some work to function in parallel.
#	User	Rev	Content
1	chuckv	1145	/*
2			* Copyright (c) 2005 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
42			/*
43			Provides the force manager for Isothermal-Isobaric Langevin Dynamics where the stochastic force
44			is applied to the surface atoms anisotropically so as to provide a constant pressure. The
45			surface atoms are determined by computing the convex hull and then triangulating that hull. The force
46			applied to the facets of the triangulation and mapped back onto the surface atoms.
47			See: Kohanoff et.al. CHEMPHYSCHEM,2005,6,1848-1852.
48			*/
49
50			#include <fstream>
51			#include <iostream>
52			#include "integrators/SMLDForceManager.hpp"
53			#include "math/CholeskyDecomposition.hpp"
54			#include "utils/OOPSEConstant.hpp"
55			#include "hydrodynamics/Sphere.hpp"
56			#include "hydrodynamics/Ellipsoid.hpp"
57	cpuglis	1204	#if defined(HAVE_QHULL) \|\| defined(HAVE_CGAL)
58			#ifdef HAVE_QHULL
59	chuckv	1145	#include "math/ConvexHull.hpp"
60	cpuglis	1204	#endif
61
62			#ifdef HAVE_CGAL
63			#include "math/AlphaShape.hpp"
64			#endif
65			#endif
66	gezelter	1210	#include "utils/ElementsTable.hpp"
67	chuckv	1145
68			namespace oopse {
69
70			SMLDForceManager::SMLDForceManager(SimInfo* info) : ForceManager(info){
71			simParams = info->getSimParams();
72			veloMunge = new Velocitizer(info);
73
74			sphericalBoundaryConditions_ = false;
75			if (simParams->getUseSphericalBoundaryConditions()) {
76			sphericalBoundaryConditions_ = true;
77			if (simParams->haveLangevinBufferRadius()) {
78			langevinBufferRadius_ = simParams->getLangevinBufferRadius();
79			} else {
80			sprintf( painCave.errMsg,
81			"langevinBufferRadius must be specified "
82			"when useSphericalBoundaryConditions is turned on.\n");
83			painCave.severity = OOPSE_ERROR;
84			painCave.isFatal = 1;
85			simError();
86			}
87
88			if (simParams->haveFrozenBufferRadius()) {
89			frozenBufferRadius_ = simParams->getFrozenBufferRadius();
90			} else {
91			sprintf( painCave.errMsg,
92			"frozenBufferRadius must be specified "
93			"when useSphericalBoundaryConditions is turned on.\n");
94			painCave.severity = OOPSE_ERROR;
95			painCave.isFatal = 1;
96			simError();
97			}
98
99			if (frozenBufferRadius_ < langevinBufferRadius_) {
100			sprintf( painCave.errMsg,
101			"frozenBufferRadius has been set smaller than the "
102			"langevinBufferRadius. This is probably an error.\n");
103			painCave.severity = OOPSE_WARNING;
104			painCave.isFatal = 0;
105			simError();
106			}
107			}
108
109			// Build the hydroProp map:
110			std::map<std::string, HydroProp*> hydroPropMap;
111
112			Molecule* mol;
113			StuntDouble* integrableObject;
114			SimInfo::MoleculeIterator i;
115			Molecule::IntegrableObjectIterator j;
116			bool needHydroPropFile = false;
117
118			for (mol = info->beginMolecule(i); mol != NULL;
119			mol = info->nextMolecule(i)) {
120			for (integrableObject = mol->beginIntegrableObject(j);
121			integrableObject != NULL;
122			integrableObject = mol->nextIntegrableObject(j)) {
123
124			if (integrableObject->isRigidBody()) {
125			RigidBody* rb = static_cast<RigidBody*>(integrableObject);
126			if (rb->getNumAtoms() > 1) needHydroPropFile = true;
127			}
128
129			}
130			}
131
132
133			if (needHydroPropFile) {
134			if (simParams->haveHydroPropFile()) {
135			hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
136			} else {
137			sprintf( painCave.errMsg,
138			"HydroPropFile must be set to a file name if Langevin\n"
139			"\tDynamics is specified for rigidBodies which contain more\n"
140			"\tthan one atom. To create a HydroPropFile, run \"Hydro\".\n");
141			painCave.severity = OOPSE_ERROR;
142			painCave.isFatal = 1;
143			simError();
144			}
145
146			for (mol = info->beginMolecule(i); mol != NULL;
147			mol = info->nextMolecule(i)) {
148			for (integrableObject = mol->beginIntegrableObject(j);
149			integrableObject != NULL;
150			integrableObject = mol->nextIntegrableObject(j)) {
151
152			std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
153			if (iter != hydroPropMap.end()) {
154			hydroProps_.push_back(iter->second);
155			} else {
156			sprintf( painCave.errMsg,
157			"Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
158			painCave.severity = OOPSE_ERROR;
159			painCave.isFatal = 1;
160			simError();
161			}
162			}
163			}
164			} else {
165
166			std::map<std::string, HydroProp*> hydroPropMap;
167			for (mol = info->beginMolecule(i); mol != NULL;
168			mol = info->nextMolecule(i)) {
169			for (integrableObject = mol->beginIntegrableObject(j);
170			integrableObject != NULL;
171			integrableObject = mol->nextIntegrableObject(j)) {
172			Shape* currShape = NULL;
173	gezelter	1210
174			if (integrableObject->isAtom()){
175			Atom* atom = static_cast<Atom*>(integrableObject);
176			AtomType* atomType = atom->getAtomType();
177	chuckv	1145	if (atomType->isGayBerne()) {
178	gezelter	1210	DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
179	chuckv	1145	GenericData* data = dAtomType->getPropertyByName("GayBerne");
180			if (data != NULL) {
181			GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
182
183			if (gayBerneData != NULL) {
184			GayBerneParam gayBerneParam = gayBerneData->getData();
185			currShape = new Ellipsoid(V3Zero,
186	gezelter	1210	gayBerneParam.GB_l / 2.0,
187	chuckv	1145	gayBerneParam.GB_d / 2.0,
188			Mat3x3d::identity());
189			} else {
190			sprintf( painCave.errMsg,
191			"Can not cast GenericData to GayBerneParam\n");
192			painCave.severity = OOPSE_ERROR;
193			painCave.isFatal = 1;
194			simError();
195			}
196			} else {
197			sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
198			painCave.severity = OOPSE_ERROR;
199			painCave.isFatal = 1;
200			simError();
201			}
202	gezelter	1210	} else {
203			if (atomType->isLennardJones()){
204			GenericData* data = atomType->getPropertyByName("LennardJones");
205			if (data != NULL) {
206			LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
207			if (ljData != NULL) {
208			LJParam ljParam = ljData->getData();
209			currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
210			} else {
211			sprintf( painCave.errMsg,
212			"Can not cast GenericData to LJParam\n");
213			painCave.severity = OOPSE_ERROR;
214			painCave.isFatal = 1;
215			simError();
216			}
217			}
218			} else {
219			int obanum = etab.GetAtomicNum((atom->getType()).c_str());
220			if (obanum != 0) {
221			currShape = new Sphere(atom->getPos(), etab.GetVdwRad(obanum));
222	chuckv	1145	} else {
223			sprintf( painCave.errMsg,
224	gezelter	1210	"Could not find atom type in default element.txt\n");
225	chuckv	1145	painCave.severity = OOPSE_ERROR;
226			painCave.isFatal = 1;
227			simError();
228	gezelter	1210	}
229	chuckv	1145	}
230			}
231			}
232			HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
233			std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
234			if (iter != hydroPropMap.end())
235			hydroProps_.push_back(iter->second);
236			else {
237			currHydroProp->complete();
238			hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
239			hydroProps_.push_back(currHydroProp);
240			}
241			}
242			}
243			}
244			variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
245			}
246
247			std::map<std::string, HydroProp*> SMLDForceManager::parseFrictionFile(const std::string& filename) {
248			std::map<std::string, HydroProp*> props;
249			std::ifstream ifs(filename.c_str());
250			if (ifs.is_open()) {
251
252			}
253
254			const unsigned int BufferSize = 65535;
255			char buffer[BufferSize];
256			while (ifs.getline(buffer, BufferSize)) {
257			HydroProp* currProp = new HydroProp(buffer);
258			props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
259			}
260
261			return props;
262			}
263
264			void SMLDForceManager::postCalculation(bool needStress){
265			SimInfo::MoleculeIterator i;
266			Molecule::IntegrableObjectIterator j;
267			Molecule* mol;
268			StuntDouble* integrableObject;
269	gezelter	1210	RealType mass;
270	chuckv	1145	Vector3d vel;
271			Vector3d pos;
272			Vector3d frc;
273			Mat3x3d A;
274			Mat3x3d Atrans;
275			Vector3d Tb;
276			Vector3d ji;
277			unsigned int index = 0;
278			bool doLangevinForces;
279			bool freezeMolecule;
280			int fdf;
281
282	gezelter	1210
283
284	chuckv	1145	fdf = 0;
285
286			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
287
288			doLangevinForces = true;
289			freezeMolecule = false;
290
291			if (sphericalBoundaryConditions_) {
292
293			Vector3d molPos = mol->getCom();
294			RealType molRad = molPos.length();
295
296			doLangevinForces = false;
297
298			if (molRad > langevinBufferRadius_) {
299			doLangevinForces = true;
300			freezeMolecule = false;
301			}
302			if (molRad > frozenBufferRadius_) {
303			doLangevinForces = false;
304			freezeMolecule = true;
305			}
306			}
307
308			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
309			integrableObject = mol->nextIntegrableObject(j)) {
310
311			if (freezeMolecule)
312			fdf += integrableObject->freeze();
313
314			if (doLangevinForces) {
315			vel =integrableObject->getVel();
316	gezelter	1210	mass = integrableObject->getMass();
317	chuckv	1145	if (integrableObject->isDirectional()){
318			//calculate angular velocity in lab frame
319			Mat3x3d I = integrableObject->getI();
320			Vector3d angMom = integrableObject->getJ();
321			Vector3d omega;
322
323			if (integrableObject->isLinear()) {
324			int linearAxis = integrableObject->linearAxis();
325			int l = (linearAxis +1 )%3;
326			int m = (linearAxis +2 )%3;
327			omega[l] = angMom[l] /I(l, l);
328			omega[m] = angMom[m] /I(m, m);
329
330			} else {
331			omega[0] = angMom[0] /I(0, 0);
332			omega[1] = angMom[1] /I(1, 1);
333			omega[2] = angMom[2] /I(2, 2);
334			}
335	gezelter	1210
336			//std::cerr << "I = " << I(0,0) << "\t" << I(1,1) << "\t" << I(2,2) << "\n\n";
337
338	chuckv	1145	//apply friction force and torque at center of resistance
339			A = integrableObject->getA();
340			Atrans = A.transpose();
341	gezelter	1210	//std::cerr << "A = " << integrableObject->getA() << "\n";
342			//std::cerr << "Atrans = " << A.transpose() << "\n\n";
343	chuckv	1145	Vector3d rcr = Atrans * hydroProps_[index]->getCOR();
344	gezelter	1210	//std::cerr << "cor = " << hydroProps_[index]->getCOR() << "\n\n\n\n";
345			//std::cerr << "rcr = " << rcr << "\n\n";
346	chuckv	1145	Vector3d vcdLab = vel + cross(omega, rcr);
347	gezelter	1210
348			//std::cerr << "velL = " << vel << "\n\n";
349			//std::cerr << "vcdL = " << vcdLab << "\n\n";
350	chuckv	1145	Vector3d vcdBody = A* vcdLab;
351	gezelter	1210	//std::cerr << "vcdB = " << vcdBody << "\n\n";
352	chuckv	1145	Vector3d frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omega);
353	gezelter	1210
354			//std::cerr << "xitt = " << hydroProps_[index]->getXitt() << "\n\n";
355			//std::cerr << "ffB = " << frictionForceBody << "\n\n";
356	chuckv	1145	Vector3d frictionForceLab = Atrans*frictionForceBody;
357	gezelter	1210	//std::cerr << "ffL = " << frictionForceLab << "\n\n";
358			//std::cerr << "frc = " << integrableObject->getFrc() << "\n\n";
359	chuckv	1145	integrableObject->addFrc(frictionForceLab);
360	gezelter	1210	//std::cerr << "frc = " << integrableObject->getFrc() << "\n\n";
361			//std::cerr << "ome = " << omega << "\n\n";
362	chuckv	1145	Vector3d frictionTorqueBody = - (hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omega);
363	gezelter	1210	//std::cerr << "ftB = " << frictionTorqueBody << "\n\n";
364	chuckv	1145	Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
365	gezelter	1210	//std::cerr << "ftL = " << frictionTorqueLab << "\n\n";
366			//std::cerr << "ftL2 = " << frictionTorqueLab+cross(rcr,frictionForceLab) << "\n\n";
367			//std::cerr << "trq = " << integrableObject->getTrq() << "\n\n";
368	chuckv	1145	integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
369	gezelter	1210	//std::cerr << "trq = " << integrableObject->getTrq() << "\n\n";
370
371	chuckv	1145	//apply random force and torque at center of resistance
372			Vector3d randomForceBody;
373			Vector3d randomTorqueBody;
374			genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
375	gezelter	1210	//std::cerr << "rfB = " << randomForceBody << "\n\n";
376			//std::cerr << "rtB = " << randomTorqueBody << "\n\n";
377	chuckv	1145	Vector3d randomForceLab = Atrans*randomForceBody;
378			Vector3d randomTorqueLab = Atrans* randomTorqueBody;
379			integrableObject->addFrc(randomForceLab);
380	gezelter	1210	//std::cerr << "rfL = " << randomForceLab << "\n\n";
381			//std::cerr << "rtL = " << randomTorqueLab << "\n\n";
382			//std::cerr << "rtL2 = " << randomTorqueLab + cross(rcr, randomForceLab) << "\n\n";
383	chuckv	1145	integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));
384
385			} else {
386			//spherical atom
387			Vector3d frictionForce = -(hydroProps_[index]->getXitt() * vel);
388	gezelter	1210	//std::cerr << "xitt = " << hydroProps_[index]->getXitt() << "\n\n";
389	chuckv	1145	Vector3d randomForce;
390			Vector3d randomTorque;
391			genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
392
393			integrableObject->addFrc(frictionForce+randomForce);
394			}
395			}
396
397			++index;
398
399			}
400			}
401
402			info_->setFdf(fdf);
403	gezelter	1210	veloMunge->removeComDrift();
404	chuckv	1145	// Remove angular drift if we are not using periodic boundary conditions.
405			if(!simParams->getUsePeriodicBoundaryConditions())
406			veloMunge->removeAngularDrift();
407
408			ForceManager::postCalculation(needStress);
409			}
410
411			void SMLDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
412
413
414			Vector<RealType, 6> Z;
415			Vector<RealType, 6> generalForce;
416
417			Z[0] = randNumGen_.randNorm(0, variance);
418			Z[1] = randNumGen_.randNorm(0, variance);
419			Z[2] = randNumGen_.randNorm(0, variance);
420			Z[3] = randNumGen_.randNorm(0, variance);
421			Z[4] = randNumGen_.randNorm(0, variance);
422			Z[5] = randNumGen_.randNorm(0, variance);
423
424
425			generalForce = hydroProps_[index]->getS()*Z;
426
427			force[0] = generalForce[0];
428			force[1] = generalForce[1];
429			force[2] = generalForce[2];
430			torque[0] = generalForce[3];
431			torque[1] = generalForce[4];
432			torque[2] = generalForce[5];
433
434	gezelter	1210	}
435	chuckv	1145
436			}