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 1390 by gezelter, Wed Nov 25 20:02:06 2009 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 6 | Line 6
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
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   #include <iostream>
# 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 < }
51 > namespace OpenMD {
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  
59 – void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
62      Vector3d newCoor;
63      Vector3d curRefCoor;  
64      RotMat3x3d rotMat = latVec2RotMat(ort);
65 <  
65 >    
66      if(mol->getNIntegrableObjects() != nIntegrableObjects){
67 <        sprintf( painCave.errMsg,
68 <            "MoLocator error.\n"
69 <            "  The number of integrable objects of MoleculeStamp is not the same as  that of Molecule\n");
70 <        painCave.isFatal = 1;
71 <        simError();
67 >      sprintf( painCave.errMsg,
68 >               "MoLocator error.\n"
69 >               "  The number of integrable objects of MoleculeStamp is not the same as  that of Molecule\n");
70 >      painCave.isFatal = 1;
71 >      simError();
72      }
73 <
73 >    
74      Molecule::IntegrableObjectIterator ii;
75      StuntDouble* integrableObject;
76      int i;
77      for (integrableObject = mol->beginIntegrableObject(ii), i = 0; integrableObject != NULL;
78 <        integrableObject = mol->nextIntegrableObject(ii), ++i) {
79 <
80 <        newCoor = rotMat * refCoords[i];
81 <        newCoor += offset;
82 <
83 <        integrableObject->setPos( newCoor);
84 <        integrableObject->setVel(V3Zero);
85 <
86 <        if(integrableObject->isDirectional()){
87 <          integrableObject->setA(rotMat * integrableObject->getA());
88 <          integrableObject->setJ(V3Zero);  
89 <        }        
78 >         integrableObject = mol->nextIntegrableObject(ii), ++i) {
79 >      
80 >      newCoor = rotMat * refCoords[i];
81 >      newCoor += offset;
82 >    
83 >      integrableObject->setPos(newCoor);
84 >      integrableObject->setVel(V3Zero);
85 >      
86 >      if(integrableObject->isDirectional()){
87 >        integrableObject->setA(rotMat * integrableObject->getA());
88 >        integrableObject->setJ(V3Zero);  
89 >      }        
90      }
91 < }
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;
91 >  }
92    
93 <  nAtoms= myStamp->getNAtoms();
94 <  nRigidBodies = myStamp->getNRigidBodies();
95 <
96 <  for(size_t i=0; i<nAtoms; i++){
97 <
98 <    currAtomStamp = myStamp->getAtom(i);
99 <
100 <    if( !currAtomStamp->havePosition() ){
101 <      sprintf( painCave.errMsg,
102 <                  "MoLocator error.\n"
103 <                  "  Component %s, atom %s does not have a position specified.\n"
104 <                  "  This means MoLocator cannot initalize it's position.\n",
105 <                  myStamp->getID(),
106 <                  currAtomStamp->getType() );
107 <
108 <      painCave.isFatal = 1;
109 <      simError();
93 >  void MoLocator::calcRef( void ){
94 >    AtomStamp* currAtomStamp;
95 >    RigidBodyStamp* rbStamp;
96 >    int nAtoms;
97 >    int nRigidBodies;
98 >    std::vector<RealType> mass;
99 >    Vector3d coor;
100 >    Vector3d refMolCom;  
101 >    int nAtomsInRb;
102 >    RealType totMassInRb;
103 >    RealType currAtomMass;
104 >    RealType molMass;
105 >    
106 >    nAtoms= myStamp->getNAtoms();
107 >    nRigidBodies = myStamp->getNRigidBodies();
108 >    
109 >    for(size_t i=0; i<nAtoms; i++){
110 >      
111 >      currAtomStamp = myStamp->getAtomStamp(i);
112 >      
113 >      if( !currAtomStamp->havePosition() ){
114 >        sprintf( painCave.errMsg,
115 >                 "MoLocator error.\n"
116 >                 "  Component %s, atom %s does not have a position specified.\n"
117 >                 "  This means MoLocator cannot initalize it's position.\n",
118 >                 myStamp->getName().c_str(),
119 >                 currAtomStamp->getType().c_str());
120 >        
121 >        painCave.isFatal = 1;
122 >        simError();
123 >      }
124 >      
125 >      //if atom belongs to rigidbody, just skip it
126 >      if(myStamp->isAtomInRigidBody(i))
127 >        continue;
128 >      //get mass and the reference coordinate
129 >      else{
130 >        currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);  
131 >        mass.push_back(currAtomMass);
132 >        coor.x() = currAtomStamp->getPosX();
133 >        coor.y() = currAtomStamp->getPosY();
134 >        coor.z() = currAtomStamp->getPosZ();
135 >        refCoords.push_back(coor);
136 >        
137 >      }
138      }
139 <
140 <    //if atom belongs to rigidbody, just skip it
141 <    if(myStamp->isAtomInRigidBody(i))
142 <      continue;
143 <    //get mass and the reference coordinate
144 <    else{
145 <      currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);  
146 <      mass.push_back(currAtomMass);
147 <      coor.x() = currAtomStamp->getPosX();
148 <      coor.y() = currAtomStamp->getPosY();
149 <      coor.z() = currAtomStamp->getPosZ();
139 >    
140 >    for(int i = 0; i < nRigidBodies; i++){
141 >      
142 >      rbStamp = myStamp->getRigidBodyStamp(i);
143 >      nAtomsInRb = rbStamp->getNMembers();
144 >      
145 >      coor.x() = 0.0;
146 >      coor.y() = 0.0;
147 >      coor.z() = 0.0;
148 >      totMassInRb = 0.0;
149 >      
150 >      for(int j = 0; j < nAtomsInRb; j++){
151 >        
152 >        currAtomStamp = myStamp->getAtomStamp(rbStamp->getMemberAt(j));
153 >        currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
154 >        totMassInRb +=  currAtomMass;
155 >        
156 >        coor.x() += currAtomStamp->getPosX() * currAtomMass;
157 >        coor.y() += currAtomStamp->getPosY() * currAtomMass;
158 >        coor.z() += currAtomStamp->getPosZ() * currAtomMass;
159 >      }
160 >      
161 >      mass.push_back(totMassInRb);
162 >      coor /= totMassInRb;
163        refCoords.push_back(coor);
133 –
164      }
165 <  }
166 <
167 <  for(int i = 0; i < nRigidBodies; i++){
168 <    coor.x() = 0;
169 <    coor.y() = 0;
170 <    coor.z() = 0;
171 <    totMassInRb = 0;
172 <
173 <    for(int j = 0; j < nAtomsInRb; j++){
174 <
175 <      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;
165 >    
166 >    
167 >    //calculate the reference center of mass
168 >    molMass = 0;
169 >    refMolCom.x() = 0;
170 >    refMolCom.y() = 0;
171 >    refMolCom.z() = 0;
172 >    
173 >    for(int i = 0; i < nIntegrableObjects; i++){
174 >      refMolCom += refCoords[i] * mass[i];
175 >      molMass += mass[i];
176      }
177 <
178 <    mass.push_back(totMassInRb);
179 <    coor /= totMassInRb;
180 <    refCoords.push_back(coor);
177 >    
178 >    refMolCom /= molMass;
179 >    
180 >    //move the reference center of mass to (0,0,0) and adjust the reference coordinate
181 >    //of the integrabel objects
182 >    for(int i = 0; i < nIntegrableObjects; i++)
183 >      refCoords[i] -= refMolCom;
184    }
157 –
158 –
159 –  //calculate the reference center of mass
160 –  molMass = 0;
161 –  refMolCom.x() = 0;
162 –  refMolCom.y() = 0;
163 –  refMolCom.z() = 0;
185    
186 <  for(int i = 0; i < nIntegrableObjects; i++){
187 <    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) {
181 <    double mass;
186 >  RealType getAtomMass(const std::string& at, ForceField* myFF) {
187 >    RealType mass;
188      AtomType* atomType= myFF->getAtomType(at);
189      if (atomType != NULL) {
190 <        mass =     atomType->getMass();
190 >      mass =     atomType->getMass();
191      } else {
192 <        mass = 0.0;
193 <        std::cerr << "Can not find AtomType: " << at << std::endl;
192 >      mass = 0.0;
193 >      std::cerr << "Can not find AtomType: " << at << std::endl;
194      }
195      return mass;
196 < }
197 <
198 < double getMolMass(MoleculeStamp *molStamp, ForceField *myFF) {
196 >  }
197 >  
198 >  RealType getMolMass(MoleculeStamp *molStamp, ForceField *myFF) {
199      int nAtoms;
200 <    double totMass = 0;
200 >    RealType totMass = 0;
201      nAtoms = molStamp->getNAtoms();
202 <
202 >    
203      for(size_t i = 0; i < nAtoms; i++) {
204 <        AtomStamp *currAtomStamp = molStamp->getAtom(i);
205 <        totMass += getAtomMass(currAtomStamp->getType(), myFF);        
204 >      AtomStamp *currAtomStamp = molStamp->getAtomStamp(i);
205 >      totMass += getAtomMass(currAtomStamp->getType(), myFF);        
206      }
207      return totMass;
208 < }
209 < RotMat3x3d latVec2RotMat(const Vector3d& lv){
210 <
211 <    double theta =acos(lv[2]);
212 <    double phi = atan2(lv[1], lv[0]);
213 <    double psi = 0;
214 <
208 >  }
209 >  RotMat3x3d latVec2RotMat(const Vector3d& lv){
210 >    
211 >    RealType theta =acos(lv[2]);
212 >    RealType phi = atan2(lv[1], lv[0]);
213 >    RealType psi = 0;
214 >    
215      return RotMat3x3d(phi, theta, psi);
216 <
216 >    
217 >  }
218   }
212 – }
219  

Diff Legend

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