ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/primitives/Molecule.cpp
Revision: 1718
Committed: Thu May 24 01:29:59 2012 UTC (12 years, 11 months ago) by gezelter
File size: 9913 byte(s)
Log Message:
Adding fluctuating charges, still a work in progress

File Contents

# User Rev Content
1 gezelter 507 /*
2 gezelter 246 * 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 gezelter 1390 * 1. Redistributions of source code must retain the above copyright
10 gezelter 246 * notice, this list of conditions and the following disclaimer.
11     *
12 gezelter 1390 * 2. Redistributions in binary form must reproduce the above copyright
13 gezelter 246 * notice, this list of conditions and the following disclaimer in the
14     * documentation and/or other materials provided with the
15     * distribution.
16     *
17     * This software is provided "AS IS," without a warranty of any
18     * kind. All express or implied conditions, representations and
19     * warranties, including any implied warranty of merchantability,
20     * fitness for a particular purpose or non-infringement, are hereby
21     * excluded. The University of Notre Dame and its licensors shall not
22     * be liable for any damages suffered by licensee as a result of
23     * using, modifying or distributing the software or its
24     * derivatives. In no event will the University of Notre Dame or its
25     * licensors be liable for any lost revenue, profit or data, or for
26     * direct, indirect, special, consequential, incidental or punitive
27     * damages, however caused and regardless of the theory of liability,
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 gezelter 1390 *
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 gezelter 1665 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40     * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 gezelter 246 */
42    
43     /**
44     * @file Molecule.cpp
45     * @author tlin
46     * @date 10/28/2004
47     * @version 1.0
48     */
49 gezelter 2
50 gezelter 246 #include <algorithm>
51     #include <set>
52 gezelter 2
53 tim 3 #include "primitives/Molecule.hpp"
54 gezelter 246 #include "utils/MemoryUtils.hpp"
55 tim 3 #include "utils/simError.h"
56 gezelter 2
57 gezelter 1390 namespace OpenMD {
58 gezelter 507 Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
59 gezelter 1718 : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName),
60     constrainTotalCharge_(false) {
61 gezelter 1211 }
62    
63 gezelter 507 Molecule::~Molecule() {
64 gezelter 1211
65 tim 398 MemoryUtils::deletePointers(atoms_);
66     MemoryUtils::deletePointers(bonds_);
67     MemoryUtils::deletePointers(bends_);
68     MemoryUtils::deletePointers(torsions_);
69 gezelter 1277 MemoryUtils::deletePointers(inversions_);
70 tim 398 MemoryUtils::deletePointers(rigidBodies_);
71     MemoryUtils::deletePointers(cutoffGroups_);
72     MemoryUtils::deletePointers(constraintPairs_);
73     MemoryUtils::deletePointers(constraintElems_);
74 gezelter 1715
75 gezelter 1211 // integrableObjects_ don't own the objects
76 gezelter 246 integrableObjects_.clear();
77 gezelter 1715 fluctuatingCharges_.clear();
78 gezelter 246
79 gezelter 507 }
80 gezelter 1211
81 gezelter 507 void Molecule::addAtom(Atom* atom) {
82 gezelter 246 if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
83 gezelter 507 atoms_.push_back(atom);
84 gezelter 246 }
85 gezelter 507 }
86 gezelter 1211
87 gezelter 507 void Molecule::addBond(Bond* bond) {
88 gezelter 246 if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
89 gezelter 507 bonds_.push_back(bond);
90 gezelter 246 }
91 gezelter 507 }
92 gezelter 1211
93 gezelter 507 void Molecule::addBend(Bend* bend) {
94 gezelter 246 if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
95 gezelter 507 bends_.push_back(bend);
96 gezelter 246 }
97 gezelter 507 }
98 gezelter 1211
99 gezelter 507 void Molecule::addTorsion(Torsion* torsion) {
100 gezelter 1211 if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
101     torsions_.end()) {
102 gezelter 507 torsions_.push_back(torsion);
103 gezelter 246 }
104 gezelter 507 }
105 gezelter 1277
106     void Molecule::addInversion(Inversion* inversion) {
107     if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
108     inversions_.end()) {
109     inversions_.push_back(inversion);
110     }
111     }
112 gezelter 1211
113 gezelter 507 void Molecule::addRigidBody(RigidBody *rb) {
114 gezelter 1211 if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
115     rigidBodies_.end()) {
116 gezelter 507 rigidBodies_.push_back(rb);
117 gezelter 246 }
118 gezelter 507 }
119 gezelter 1211
120 gezelter 507 void Molecule::addCutoffGroup(CutoffGroup* cp) {
121 gezelter 1211 if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
122     cutoffGroups_.end()) {
123 gezelter 507 cutoffGroups_.push_back(cp);
124 gezelter 1211 }
125 gezelter 507 }
126 gezelter 1211
127 gezelter 507 void Molecule::addConstraintPair(ConstraintPair* cp) {
128 gezelter 1211 if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
129     constraintPairs_.end()) {
130 gezelter 507 constraintPairs_.push_back(cp);
131 gezelter 1211 }
132 gezelter 507 }
133 gezelter 1211
134 gezelter 507 void Molecule::addConstraintElem(ConstraintElem* cp) {
135 gezelter 1211 if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
136     constraintElems_.end()) {
137 gezelter 507 constraintElems_.push_back(cp);
138 gezelter 246 }
139 gezelter 507 }
140 gezelter 1211
141 gezelter 507 void Molecule::complete() {
142 gezelter 246
143     std::set<Atom*> rigidAtoms;
144 gezelter 1715 Atom* atom;
145     AtomIterator ai;
146 gezelter 246 RigidBody* rb;
147 gezelter 1715 RigidBodyIterator rbIter;
148 gezelter 2
149 gezelter 1715 // Get list of all the atoms that are part of rigid bodies
150    
151 gezelter 246 for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
152 gezelter 507 rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
153 gezelter 246 }
154 gezelter 1211
155 gezelter 1715 // add any atom that wasn't part of a rigid body to the list of integrableObjects
156    
157 gezelter 246 for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
158 gezelter 1211
159 gezelter 1715 if (rigidAtoms.find(atom) == rigidAtoms.end()) {
160 gezelter 1211
161     // If an atom does not belong to a rigid body, it is an
162     // integrable object
163    
164 gezelter 1715 integrableObjects_.push_back(atom);
165 gezelter 507 }
166 gezelter 246 }
167 gezelter 1211
168 gezelter 1715 // then add the rigid bodies themselves to the integrableObjects
169 gezelter 2
170 tim 372 for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
171     integrableObjects_.push_back(rb);
172     }
173 gezelter 1715
174     // find the atoms that are fluctuating charges and add them to the
175     // fluctuatingCharges_ vector
176    
177     for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
178     if ( atom->isFluctuatingCharge() )
179     fluctuatingCharges_.push_back( atom );
180     }
181    
182 gezelter 507 }
183 gezelter 2
184 tim 963 RealType Molecule::getMass() {
185 gezelter 246 StuntDouble* sd;
186     std::vector<StuntDouble*>::iterator i;
187 tim 963 RealType mass = 0.0;
188 gezelter 1211
189     for (sd = beginIntegrableObject(i); sd != NULL; sd =
190     nextIntegrableObject(i)){
191 gezelter 507 mass += sd->getMass();
192 gezelter 246 }
193 gezelter 1211
194     return mass;
195 gezelter 507 }
196 gezelter 2
197 gezelter 507 Vector3d Molecule::getCom() {
198 gezelter 246 StuntDouble* sd;
199     std::vector<StuntDouble*>::iterator i;
200     Vector3d com;
201 tim 963 RealType totalMass = 0;
202     RealType mass;
203 gezelter 246
204 gezelter 1211 for (sd = beginIntegrableObject(i); sd != NULL; sd =
205     nextIntegrableObject(i)){
206 gezelter 507 mass = sd->getMass();
207     totalMass += mass;
208     com += sd->getPos() * mass;
209 gezelter 2 }
210 gezelter 1211
211 gezelter 246 com /= totalMass;
212 gezelter 2
213 gezelter 246 return com;
214 gezelter 507 }
215 gezelter 2
216 gezelter 507 void Molecule::moveCom(const Vector3d& delta) {
217 gezelter 246 StuntDouble* sd;
218     std::vector<StuntDouble*>::iterator i;
219    
220 gezelter 1211 for (sd = beginIntegrableObject(i); sd != NULL; sd =
221     nextIntegrableObject(i)){
222 gezelter 507 sd->setPos(sd->getPos() + delta);
223 gezelter 1211 }
224 gezelter 507 }
225 gezelter 2
226 gezelter 507 Vector3d Molecule::getComVel() {
227 gezelter 246 StuntDouble* sd;
228     std::vector<StuntDouble*>::iterator i;
229     Vector3d velCom;
230 tim 963 RealType totalMass = 0;
231     RealType mass;
232 gezelter 246
233 gezelter 1211 for (sd = beginIntegrableObject(i); sd != NULL; sd =
234     nextIntegrableObject(i)){
235 gezelter 507 mass = sd->getMass();
236     totalMass += mass;
237     velCom += sd->getVel() * mass;
238 gezelter 246 }
239 gezelter 1211
240 gezelter 246 velCom /= totalMass;
241 gezelter 1211
242 gezelter 246 return velCom;
243 gezelter 507 }
244 gezelter 2
245 tim 963 RealType Molecule::getPotential() {
246 gezelter 2
247 gezelter 246 Bond* bond;
248     Bend* bend;
249     Torsion* torsion;
250 gezelter 1277 Inversion* inversion;
251 gezelter 246 Molecule::BondIterator bondIter;;
252     Molecule::BendIterator bendIter;
253     Molecule::TorsionIterator torsionIter;
254 gezelter 1277 Molecule::InversionIterator inversionIter;
255 gezelter 2
256 tim 963 RealType potential = 0.0;
257 gezelter 2
258 gezelter 246 for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
259 gezelter 507 potential += bond->getPotential();
260 gezelter 246 }
261 gezelter 2
262 gezelter 246 for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
263 gezelter 507 potential += bend->getPotential();
264 gezelter 2 }
265    
266 gezelter 1211 for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
267     nextTorsion(torsionIter)) {
268 gezelter 507 potential += torsion->getPotential();
269 gezelter 2 }
270 gezelter 1211
271 gezelter 1277 for (inversion = beginInversion(inversionIter); torsion != NULL;
272     inversion = nextInversion(inversionIter)) {
273     potential += inversion->getPotential();
274     }
275    
276 gezelter 246 return potential;
277 gezelter 1211
278 gezelter 507 }
279 gezelter 1211
280 cli2 1360 void Molecule::addProperty(GenericData* genData) {
281     properties_.addProperty(genData);
282     }
283    
284     void Molecule::removeProperty(const std::string& propName) {
285     properties_.removeProperty(propName);
286     }
287    
288     void Molecule::clearProperties() {
289     properties_.clearProperties();
290     }
291    
292     std::vector<std::string> Molecule::getPropertyNames() {
293     return properties_.getPropertyNames();
294     }
295    
296     std::vector<GenericData*> Molecule::getProperties() {
297     return properties_.getProperties();
298     }
299    
300     GenericData* Molecule::getPropertyByName(const std::string& propName) {
301     return properties_.getPropertyByName(propName);
302     }
303    
304 gezelter 507 std::ostream& operator <<(std::ostream& o, Molecule& mol) {
305 gezelter 246 o << std::endl;
306     o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
307     o << mol.getNAtoms() << " atoms" << std::endl;
308     o << mol.getNBonds() << " bonds" << std::endl;
309     o << mol.getNBends() << " bends" << std::endl;
310     o << mol.getNTorsions() << " torsions" << std::endl;
311 gezelter 1277 o << mol.getNInversions() << " inversions" << std::endl;
312 gezelter 246 o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
313 gezelter 1715 o << mol.getNIntegrableObjects() << " integrable objects" << std::endl;
314     o << mol.getNCutoffGroups() << " cutoff groups" << std::endl;
315     o << mol.getNConstraintPairs() << " constraint pairs" << std::endl;
316     o << mol.getNFluctuatingCharges() << " fluctuating charges" << std::endl;
317 gezelter 246 return o;
318 gezelter 507 }
319 gezelter 1211
320 gezelter 1390 }//end namespace OpenMD

Properties

Name Value
svn:keywords Author Id Revision Date