ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/utils/MoLocator.cpp
(Generate patch)

Comparing trunk/src/utils/MoLocator.cpp (file contents):
Revision 490 by tim, Wed Apr 13 18:41:17 2005 UTC vs.
Revision 770 by tim, Fri Dec 2 15:38:03 2005 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 47 | Line 47
47   #include "utils/simError.h"
48   #include "utils/MoLocator.hpp"
49   #include "types/AtomType.hpp"
50 – namespace oopse {
51 – MoLocator::MoLocator( MoleculeStamp* theStamp, ForceField* theFF){
50  
51 <  myStamp = theStamp;
52 <  myFF = theFF;
53 <  nIntegrableObjects = myStamp->getNIntegrable();
54 <  calcRef();
55 < }
56 <
57 < void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
51 > namespace oopse {
52 >  MoLocator::MoLocator( MoleculeStamp* theStamp, ForceField* theFF){
53 >    
54 >    myStamp = theStamp;
55 >    myFF = theFF;
56 >    nIntegrableObjects = myStamp->getNIntegrable();
57 >    calcRef();
58 >  }
59 >  
60 >  void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
61      Vector3d newCoor;
62      Vector3d curRefCoor;  
63      RotMat3x3d rotMat = latVec2RotMat(ort);
64 <  
64 >    
65      if(mol->getNIntegrableObjects() != nIntegrableObjects){
66 <        sprintf( painCave.errMsg,
67 <            "MoLocator error.\n"
68 <            "  The number of integrable objects of MoleculeStamp is not the same as  that of Molecule\n");
69 <        painCave.isFatal = 1;
70 <        simError();
66 >      sprintf( painCave.errMsg,
67 >               "MoLocator error.\n"
68 >               "  The number of integrable objects of MoleculeStamp is not the same as  that of Molecule\n");
69 >      painCave.isFatal = 1;
70 >      simError();
71      }
72 <
72 >    
73      Molecule::IntegrableObjectIterator ii;
74      StuntDouble* integrableObject;
75      int i;
76      for (integrableObject = mol->beginIntegrableObject(ii), i = 0; integrableObject != NULL;
77 <        integrableObject = mol->nextIntegrableObject(ii), ++i) {
78 <
79 <        newCoor = rotMat * refCoords[i];
80 <        newCoor += offset;
81 <
82 <        integrableObject->setPos( newCoor);
83 <        integrableObject->setVel(V3Zero);
84 <
85 <        if(integrableObject->isDirectional()){
86 <          integrableObject->setA(rotMat * integrableObject->getA());
87 <          integrableObject->setJ(V3Zero);  
88 <        }        
77 >         integrableObject = mol->nextIntegrableObject(ii), ++i) {
78 >      
79 >      newCoor = rotMat * refCoords[i];
80 >      newCoor += offset;
81 >      
82 >      integrableObject->setPos( newCoor);
83 >      integrableObject->setVel(V3Zero);
84 >      
85 >      if(integrableObject->isDirectional()){
86 >        integrableObject->setA(rotMat * integrableObject->getA());
87 >        integrableObject->setJ(V3Zero);  
88 >      }        
89      }
90 < }
90 <
91 < void MoLocator::calcRef( void ){
92 <  AtomStamp* currAtomStamp;
93 <  int nAtoms;
94 <  int nRigidBodies;
95 <   std::vector<double> mass;
96 <  Vector3d coor;
97 <  Vector3d refMolCom;  
98 <  int nAtomsInRb;
99 <  double totMassInRb;
100 <  double currAtomMass;
101 <  double molMass;
90 >  }
91    
92 <  nAtoms= myStamp->getNAtoms();
93 <  nRigidBodies = myStamp->getNRigidBodies();
94 <
95 <  for(size_t i=0; i<nAtoms; i++){
96 <
97 <    currAtomStamp = myStamp->getAtom(i);
98 <
99 <    if( !currAtomStamp->havePosition() ){
100 <      sprintf( painCave.errMsg,
101 <                  "MoLocator error.\n"
102 <                  "  Component %s, atom %s does not have a position specified.\n"
103 <                  "  This means MoLocator cannot initalize it's position.\n",
104 <                  myStamp->getID(),
105 <                  currAtomStamp->getType() );
106 <
107 <      painCave.isFatal = 1;
108 <      simError();
92 >  void MoLocator::calcRef( void ){
93 >    AtomStamp* currAtomStamp;
94 >    RigidBodyStamp* rbStamp;
95 >    int nAtoms;
96 >    int nRigidBodies;
97 >    std::vector<double> mass;
98 >    Vector3d coor;
99 >    Vector3d refMolCom;  
100 >    int nAtomsInRb;
101 >    double totMassInRb;
102 >    double currAtomMass;
103 >    double molMass;
104 >    
105 >    nAtoms= myStamp->getNAtoms();
106 >    nRigidBodies = myStamp->getNRigidBodies();
107 >    
108 >    for(size_t i=0; i<nAtoms; i++){
109 >      
110 >      currAtomStamp = myStamp->getAtomStamp(i);
111 >      
112 >      if( !currAtomStamp->havePosition() ){
113 >        sprintf( painCave.errMsg,
114 >                 "MoLocator error.\n"
115 >                 "  Component %s, atom %s does not have a position specified.\n"
116 >                 "  This means MoLocator cannot initalize it's position.\n",
117 >                 myStamp->getName().c_str(),
118 >                 currAtomStamp->getType().c_str());
119 >        
120 >        painCave.isFatal = 1;
121 >        simError();
122 >      }
123 >      
124 >      //if atom belongs to rigidbody, just skip it
125 >      if(myStamp->isAtomInRigidBody(i))
126 >        continue;
127 >      //get mass and the reference coordinate
128 >      else{
129 >        currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);  
130 >        mass.push_back(currAtomMass);
131 >        coor.x() = currAtomStamp->getPosX();
132 >        coor.y() = currAtomStamp->getPosY();
133 >        coor.z() = currAtomStamp->getPosZ();
134 >        refCoords.push_back(coor);
135 >        
136 >      }
137      }
138 <
139 <    //if atom belongs to rigidbody, just skip it
140 <    if(myStamp->isAtomInRigidBody(i))
141 <      continue;
142 <    //get mass and the reference coordinate
143 <    else{
144 <      currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);  
145 <      mass.push_back(currAtomMass);
146 <      coor.x() = currAtomStamp->getPosX();
147 <      coor.y() = currAtomStamp->getPosY();
148 <      coor.z() = currAtomStamp->getPosZ();
138 >    
139 >    for(int i = 0; i < nRigidBodies; i++){
140 >      
141 >      rbStamp = myStamp->getRigidBodyStamp(i);
142 >      nAtomsInRb = rbStamp->getNMembers();
143 >      
144 >      coor.x() = 0.0;
145 >      coor.y() = 0.0;
146 >      coor.z() = 0.0;
147 >      totMassInRb = 0.0;
148 >      
149 >      for(int j = 0; j < nAtomsInRb; j++){
150 >        
151 >        currAtomStamp = myStamp->getAtomStamp(rbStamp->getMemberAt(j));
152 >        currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
153 >        totMassInRb +=  currAtomMass;
154 >        
155 >        coor.x() += currAtomStamp->getPosX() * currAtomMass;
156 >        coor.y() += currAtomStamp->getPosY() * currAtomMass;
157 >        coor.z() += currAtomStamp->getPosZ() * currAtomMass;
158 >      }
159 >      
160 >      mass.push_back(totMassInRb);
161 >      coor /= totMassInRb;
162        refCoords.push_back(coor);
133 –
163      }
164 <  }
165 <
166 <  for(int i = 0; i < nRigidBodies; i++){
167 <    coor.x() = 0;
168 <    coor.y() = 0;
169 <    coor.z() = 0;
170 <    totMassInRb = 0;
171 <
172 <    for(int j = 0; j < nAtomsInRb; j++){
173 <
174 <      currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
146 <      totMassInRb +=  currAtomMass;
147 <      
148 <      coor.x() += currAtomStamp->getPosX() * currAtomMass;
149 <      coor.y() += currAtomStamp->getPosY() * currAtomMass;
150 <      coor.z() += currAtomStamp->getPosZ() * currAtomMass;
164 >    
165 >    
166 >    //calculate the reference center of mass
167 >    molMass = 0;
168 >    refMolCom.x() = 0;
169 >    refMolCom.y() = 0;
170 >    refMolCom.z() = 0;
171 >    
172 >    for(int i = 0; i < nIntegrableObjects; i++){
173 >      refMolCom += refCoords[i] * mass[i];
174 >      molMass += mass[i];
175      }
176 <
177 <    mass.push_back(totMassInRb);
178 <    coor /= totMassInRb;
179 <    refCoords.push_back(coor);
176 >    
177 >    refMolCom /= molMass;
178 >    
179 >    //move the reference center of mass to (0,0,0) and adjust the reference coordinate
180 >    //of the integrabel objects
181 >    for(int i = 0; i < nIntegrableObjects; i++)
182 >      refCoords[i] -= refMolCom;
183    }
157 –
158 –
159 –  //calculate the reference center of mass
160 –  molMass = 0;
161 –  refMolCom.x() = 0;
162 –  refMolCom.y() = 0;
163 –  refMolCom.z() = 0;
184    
185 <  for(int i = 0; i < nIntegrableObjects; i++){
166 <    refMolCom += refCoords[i] * mass[i];
167 <   molMass += mass[i];
168 <  }
169 <  
170 <  refMolCom /= molMass;
171 <
172 <  //move the reference center of mass to (0,0,0) and adjust the reference coordinate
173 <  //of the integrabel objects
174 <  for(int i = 0; i < nIntegrableObjects; i++)
175 <    refCoords[i] -= refMolCom;
176 < }
177 <
178 <
179 <
180 < double getAtomMass(const std::string& at, ForceField* myFF) {
185 >  double getAtomMass(const std::string& at, ForceField* myFF) {
186      double mass;
187      AtomType* atomType= myFF->getAtomType(at);
188      if (atomType != NULL) {
189 <        mass =     atomType->getMass();
189 >      mass =     atomType->getMass();
190      } else {
191 <        mass = 0.0;
192 <        std::cerr << "Can not find AtomType: " << at << std::endl;
191 >      mass = 0.0;
192 >      std::cerr << "Can not find AtomType: " << at << std::endl;
193      }
194      return mass;
195 < }
196 <
197 < double getMolMass(MoleculeStamp *molStamp, ForceField *myFF) {
195 >  }
196 >  
197 >  double getMolMass(MoleculeStamp *molStamp, ForceField *myFF) {
198      int nAtoms;
199      double totMass = 0;
200      nAtoms = molStamp->getNAtoms();
201 <
201 >    
202      for(size_t i = 0; i < nAtoms; i++) {
203 <        AtomStamp *currAtomStamp = molStamp->getAtom(i);
204 <        totMass += getAtomMass(currAtomStamp->getType(), myFF);        
203 >      AtomStamp *currAtomStamp = molStamp->getAtomStamp(i);
204 >      totMass += getAtomMass(currAtomStamp->getType(), myFF);        
205      }
206      return totMass;
207 < }
208 < RotMat3x3d latVec2RotMat(const Vector3d& lv){
209 <
207 >  }
208 >  RotMat3x3d latVec2RotMat(const Vector3d& lv){
209 >    
210      double theta =acos(lv[2]);
211      double phi = atan2(lv[1], lv[0]);
212      double psi = 0;
213 <
213 >    
214      return RotMat3x3d(phi, theta, psi);
215 <
215 >    
216 >  }
217   }
212 – }
218  

Diff Legend

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