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:	1210
Committed:	Wed Jan 23 03:45:33 2008 UTC (17 years, 3 months ago) by gezelter
Original Path:	*trunk/src/integrators/LDForceManager.cpp*
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	895	/*
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	1120	#include <iostream>
43	tim	895	#include "integrators/LDForceManager.hpp"
44			#include "math/CholeskyDecomposition.hpp"
45	tim	904	#include "utils/OOPSEConstant.hpp"
46	gezelter	956	#include "hydrodynamics/Sphere.hpp"
47			#include "hydrodynamics/Ellipsoid.hpp"
48	gezelter	1210	#include "utils/ElementsTable.hpp"
49	gezelter	956
50	tim	895	namespace oopse {
51
52			LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info){
53	gezelter	983	simParams = info->getSimParams();
54			veloMunge = new Velocitizer(info);
55
56	gezelter	945	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	895
81	gezelter	945	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	956
91			// Build the hydroProp map:
92	gezelter	981	std::map<std::string, HydroProp*> hydroPropMap;
93	gezelter	956
94	tim	895	Molecule* mol;
95			StuntDouble* integrableObject;
96	gezelter	956	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	945	integrableObject = mol->nextIntegrableObject(j)) {
105	gezelter	956
106			if (integrableObject->isRigidBody()) {
107			RigidBody* rb = static_cast<RigidBody*>(integrableObject);
108			if (rb->getNumAtoms() > 1) needHydroPropFile = true;
109	gezelter	945	}
110
111			}
112	tim	895	}
113	gezelter	956
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	971
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	981	std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
135	tim	971	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	956	}
146			} else {
147	gezelter	981
148			std::map<std::string, HydroProp*> hydroPropMap;
149	gezelter	956	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	1185
156			if (integrableObject->isAtom()){
157			Atom* atom = static_cast<Atom*>(integrableObject);
158			AtomType* atomType = atom->getAtomType();
159	gezelter	956	if (atomType->isGayBerne()) {
160	xsun	1185	DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
161	gezelter	956	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	1185	gayBerneParam.GB_l / 2.0,
169	gezelter	981	gayBerneParam.GB_d / 2.0,
170	gezelter	956	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	1185	} 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	956	} else {
205			sprintf( painCave.errMsg,
206	xsun	1185	"Could not find atom type in default element.txt\n");
207	gezelter	956	painCave.severity = OOPSE_ERROR;
208			painCave.isFatal = 1;
209			simError();
210	xsun	1185	}
211	gezelter	956	}
212			}
213			}
214	gezelter	981	HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
215			std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
216	gezelter	956	if (iter != hydroPropMap.end())
217			hydroProps_.push_back(iter->second);
218			else {
219	gezelter	981	currHydroProp->complete();
220			hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
221			hydroProps_.push_back(currHydroProp);
222	gezelter	956	}
223			}
224			}
225			}
226	tim	904	variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
227	gezelter	981	}
228	gezelter	956
229	gezelter	981	std::map<std::string, HydroProp*> LDForceManager::parseFrictionFile(const std::string& filename) {
230			std::map<std::string, HydroProp*> props;
231	tim	895	std::ifstream ifs(filename.c_str());
232			if (ifs.is_open()) {
233	gezelter	945
234	tim	895	}
235	gezelter	945
236	tim	895	const unsigned int BufferSize = 65535;
237			char buffer[BufferSize];
238			while (ifs.getline(buffer, BufferSize)) {
239	gezelter	981	HydroProp* currProp = new HydroProp(buffer);
240			props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
241	tim	895	}
242	gezelter	981
243	tim	895	return props;
244			}
245	gezelter	981
246	gezelter	1126	void LDForceManager::postCalculation(bool needStress){
247	tim	895	SimInfo::MoleculeIterator i;
248			Molecule::IntegrableObjectIterator j;
249			Molecule* mol;
250			StuntDouble* integrableObject;
251	xsun	1185	RealType mass;
252	tim	895	Vector3d vel;
253			Vector3d pos;
254			Vector3d frc;
255			Mat3x3d A;
256	tim	904	Mat3x3d Atrans;
257	tim	895	Vector3d Tb;
258			Vector3d ji;
259			unsigned int index = 0;
260	gezelter	945	bool doLangevinForces;
261			bool freezeMolecule;
262			int fdf;
263	gezelter	983
264	xsun	1185
265
266	chuckv	1120	fdf = 0;
267	gezelter	983
268	tim	895	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
269	gezelter	970
270			doLangevinForces = true;
271			freezeMolecule = false;
272
273	gezelter	945	if (sphericalBoundaryConditions_) {
274
275			Vector3d molPos = mol->getCom();
276	tim	963	RealType molRad = molPos.length();
277	chuckv	1120
278	gezelter	945	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	956	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
291			integrableObject = mol->nextIntegrableObject(j)) {
292	gezelter	945
293	gezelter	956	if (freezeMolecule)
294			fdf += integrableObject->freeze();
295
296	chuckv	1120	if (doLangevinForces) {
297	tim	895	vel =integrableObject->getVel();
298	xsun	1185	mass = integrableObject->getMass();
299	tim	895	if (integrableObject->isDirectional()){
300	gezelter	945	//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	1185
318			//std::cerr << "I = " << I(0,0) << "\t" << I(1,1) << "\t" << I(2,2) << "\n\n";
319
320	gezelter	945	//apply friction force and torque at center of resistance
321			A = integrableObject->getA();
322			Atrans = A.transpose();
323	xsun	1185	//std::cerr << "A = " << integrableObject->getA() << "\n";
324			//std::cerr << "Atrans = " << A.transpose() << "\n\n";
325	gezelter	981	Vector3d rcr = Atrans * hydroProps_[index]->getCOR();
326	xsun	1185	//std::cerr << "cor = " << hydroProps_[index]->getCOR() << "\n\n\n\n";
327			//std::cerr << "rcr = " << rcr << "\n\n";
328	gezelter	945	Vector3d vcdLab = vel + cross(omega, rcr);
329	xsun	1185
330			//std::cerr << "velL = " << vel << "\n\n";
331			//std::cerr << "vcdL = " << vcdLab << "\n\n";
332	gezelter	945	Vector3d vcdBody = A* vcdLab;
333	xsun	1185	//std::cerr << "vcdB = " << vcdBody << "\n\n";
334	gezelter	981	Vector3d frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omega);
335	xsun	1185
336			//std::cerr << "xitt = " << hydroProps_[index]->getXitt() << "\n\n";
337			//std::cerr << "ffB = " << frictionForceBody << "\n\n";
338	gezelter	945	Vector3d frictionForceLab = Atrans*frictionForceBody;
339	xsun	1185	//std::cerr << "ffL = " << frictionForceLab << "\n\n";
340			//std::cerr << "frc = " << integrableObject->getFrc() << "\n\n";
341	gezelter	945	integrableObject->addFrc(frictionForceLab);
342	xsun	1185	//std::cerr << "frc = " << integrableObject->getFrc() << "\n\n";
343			//std::cerr << "ome = " << omega << "\n\n";
344	gezelter	981	Vector3d frictionTorqueBody = - (hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omega);
345	xsun	1185	//std::cerr << "ftB = " << frictionTorqueBody << "\n\n";
346	gezelter	945	Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
347	xsun	1185	//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	945	integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
351	xsun	1185	//std::cerr << "trq = " << integrableObject->getTrq() << "\n\n";
352
353	gezelter	945	//apply random force and torque at center of resistance
354			Vector3d randomForceBody;
355			Vector3d randomTorqueBody;
356			genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
357	xsun	1185	//std::cerr << "rfB = " << randomForceBody << "\n\n";
358			//std::cerr << "rtB = " << randomTorqueBody << "\n\n";
359	gezelter	945	Vector3d randomForceLab = Atrans*randomForceBody;
360			Vector3d randomTorqueLab = Atrans* randomTorqueBody;
361			integrableObject->addFrc(randomForceLab);
362	xsun	1185	//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	945	integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));
366
367	tim	895	} else {
368	gezelter	945	//spherical atom
369	gezelter	981	Vector3d frictionForce = -(hydroProps_[index]->getXitt() * vel);
370	xsun	1185	//std::cerr << "xitt = " << hydroProps_[index]->getXitt() << "\n\n";
371	gezelter	945	Vector3d randomForce;
372			Vector3d randomTorque;
373			genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
374
375			integrableObject->addFrc(frictionForce+randomForce);
376	tim	895	}
377	gezelter	956	}
378	gezelter	945
379	gezelter	956	++index;
380	tim	895
381			}
382	gezelter	956	}
383	chuckv	1120
384	gezelter	945	info_->setFdf(fdf);
385	gezelter	983	veloMunge->removeComDrift();
386			// Remove angular drift if we are not using periodic boundary conditions.
387			if(!simParams->getUsePeriodicBoundaryConditions())
388			veloMunge->removeAngularDrift();
389
390	gezelter	1126	ForceManager::postCalculation(needStress);
391	tim	895	}
392
393	tim	963	void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
394	tim	904
395	tim	906
396	tim	963	Vector<RealType, 6> Z;
397			Vector<RealType, 6> generalForce;
398	tim	904
399	tim	895	Z[0] = randNumGen_.randNorm(0, variance);
400			Z[1] = randNumGen_.randNorm(0, variance);
401			Z[2] = randNumGen_.randNorm(0, variance);
402	xsun	1185	//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	895	Z[3] = randNumGen_.randNorm(0, variance);
406			Z[4] = randNumGen_.randNorm(0, variance);
407			Z[5] = randNumGen_.randNorm(0, variance);
408	tim	904
409
410	gezelter	981	generalForce = hydroProps_[index]->getS()*Z;
411	tim	904
412	tim	895	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	1185	}
420	tim	895
421			}