src/integrators/LDForceManager.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.
 */
#include <fstream> 
#include <iostream>
#include "integrators/LDForceManager.hpp"
#include "math/CholeskyDecomposition.hpp"
#include "utils/OOPSEConstant.hpp"
#include "hydrodynamics/Sphere.hpp"
#include "hydrodynamics/Ellipsoid.hpp"
#include "utils/ElementsTable.hpp"

namespace oopse {

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