ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/integrators/SMIPDForceManager.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/integrators/SMIPDForceManager.cpp (file contents):
Revision 3515 by gezelter, Mon Jul 20 17:50:40 2009 UTC vs.
Revision 3550 by gezelter, Wed Oct 21 02:49:43 2009 UTC

# Line 1 | Line 1
1   /*
2 < * Copyright (c) 2008 The University of Notre Dame. All Rights Reserved.
2 > * Copyright (c) 2008, 2009 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
# Line 50 | Line 50 | namespace oopse {
50    SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info) {
51  
52      simParams = info->getSimParams();
53 –    thermo = new Thermo(info);
53      veloMunge = new Velocitizer(info);
54      
55      // Create Hull, Convex Hull for now, other options later.
# Line 63 | Line 62 | namespace oopse {
62      if (!simParams->haveTargetTemp()) {
63        sprintf(painCave.errMsg,
64                "SMIPDynamics error: You can't use the SMIPD integrator\n"
65 <              "   without a targetTemp!\n");      
65 >              "\twithout a targetTemp (K)!\n");      
66        painCave.isFatal = 1;
67        painCave.severity = OOPSE_ERROR;
68        simError();
# Line 74 | Line 73 | namespace oopse {
73      if (!simParams->haveTargetPressure()) {
74        sprintf(painCave.errMsg,
75                "SMIPDynamics error: You can't use the SMIPD integrator\n"
76 <              "   without a targetPressure!\n");      
76 >              "\twithout a targetPressure (atm)!\n");      
77        painCave.isFatal = 1;
78        simError();
79      } else {
# Line 86 | Line 85 | namespace oopse {
85      if (simParams->getUsePeriodicBoundaryConditions()) {
86        sprintf(painCave.errMsg,
87                "SMIPDynamics error: You can't use the SMIPD integrator\n"
88 <              "   with periodic boundary conditions!\n");    
88 >              "\twith periodic boundary conditions!\n");    
89        painCave.isFatal = 1;
90        simError();
91      }
# Line 94 | Line 93 | namespace oopse {
93      if (!simParams->haveThermalConductivity()) {
94        sprintf(painCave.errMsg,
95                "SMIPDynamics error: You can't use the SMIPD integrator\n"
96 <              "   without a thermalConductivity!\n");
96 >              "\twithout a thermalConductivity (W m^-1 K^-1)!\n");
97        painCave.isFatal = 1;
98        painCave.severity = OOPSE_ERROR;
99        simError();
100      }else{
101 <      thermalConductivity_ = simParams->getThermalConductivity();
101 >      thermalConductivity_ = simParams->getThermalConductivity() *
102 >        OOPSEConstant::thermalConductivityConvert;
103      }
104  
105      if (!simParams->haveThermalLength()) {
106        sprintf(painCave.errMsg,
107                "SMIPDynamics error: You can't use the SMIPD integrator\n"
108 <              "   without a thermalLength!\n");
108 >              "\twithout a thermalLength (Angstroms)!\n");
109        painCave.isFatal = 1;
110        painCave.severity = OOPSE_ERROR;
111        simError();
# Line 145 | Line 145 | namespace oopse {
145  
146      // Get total area and number of surface stunt doubles
147      RealType area = surfaceMesh_->getArea();
148 –    int nSurfaceSDs = surfaceMesh_->getNs();
148      std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
149      int nTriangles = sMesh.size();
150  
151      // Generate all of the necessary random forces
152      std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
153  
155 –    RealType instaTemp = thermo->getTemperature();
156 –
154      // Loop over the mesh faces and apply external pressure to each
155      // of the faces
156      std::vector<Triangle>::iterator face;
157      std::vector<StuntDouble*>::iterator vertex;
158      int thisFacet = 0;
159      for (face = sMesh.begin(); face != sMesh.end(); ++face){
163 –    
160        Triangle thisTriangle = *face;
161        std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
162        RealType thisArea = thisTriangle.getArea();
# Line 171 | Line 167 | namespace oopse {
167        RealType thisMass = thisTriangle.getFacetMass();
168  
169        // gamma is the drag coefficient normal to the face of the triangle      
170 <      RealType gamma = thermalConductivity_ * thisMass * thisArea  
171 <        / (2.0 * thermalLength_ * OOPSEConstant::kb);
172 <      
177 <      gamma *= fabs(1.0 - targetTemp_/instaTemp);      
178 <      
170 >      RealType gamma = thermalConductivity_ * thisMass * thisArea
171 >        / (2.0 * thermalLength_ * OOPSEConstant::kB);
172 >
173        RealType extPressure = - (targetPressure_ * thisArea) /
174          OOPSEConstant::energyConvert;
181 –
175        RealType randomForce = randNums[thisFacet++] * sqrt(gamma);
176        RealType dragForce = -gamma * dot(facetVel, unitNormal);
177  
178        Vector3d langevinForce = (extPressure + randomForce + dragForce) *
179          unitNormal;
180        
181 <      // Apply triangle force to stuntdouble vertices
181 >      // Apply triangle force to stuntdouble vertices
182        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
183 <        if ((*vertex) != NULL){
183 >        if ((*vertex) != NULL){
184            Vector3d vertexForce = langevinForce / 3.0;
185            (*vertex)->addFrc(vertexForce);          
193 –          if ((*vertex)->isDirectional()){          
194 –            Vector3d vertexPos = (*vertex)->getPos();
195 –            Vector3d vertexCentroidVector = vertexPos - centroid;
196 –            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
197 –          }
186          }  
187        }
188      }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines