ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/primitives/Molecule.cpp
Revision: 1908
Committed: Fri Jul 19 21:25:45 2013 UTC (11 years, 9 months ago) by gezelter
File size: 10095 byte(s)
Log Message:
Added infrastructure for region-constrained fluctuating charges.

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 gezelter 1879 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39 gezelter 1782 * [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 1908 Molecule::Molecule(int stampId, int globalIndex, const std::string& molName,
59     int region) : stampId_(stampId),
60     globalIndex_(globalIndex),
61     moleculeName_(molName),
62     region_(region),
63     constrainTotalCharge_(false) {
64 gezelter 1211 }
65    
66 gezelter 507 Molecule::~Molecule() {
67 gezelter 1211
68 tim 398 MemoryUtils::deletePointers(atoms_);
69     MemoryUtils::deletePointers(bonds_);
70     MemoryUtils::deletePointers(bends_);
71     MemoryUtils::deletePointers(torsions_);
72 gezelter 1277 MemoryUtils::deletePointers(inversions_);
73 tim 398 MemoryUtils::deletePointers(rigidBodies_);
74     MemoryUtils::deletePointers(cutoffGroups_);
75     MemoryUtils::deletePointers(constraintPairs_);
76     MemoryUtils::deletePointers(constraintElems_);
77 gezelter 1782
78 gezelter 1211 // integrableObjects_ don't own the objects
79 gezelter 246 integrableObjects_.clear();
80 gezelter 1782 fluctuatingCharges_.clear();
81 gezelter 246
82 gezelter 507 }
83 gezelter 1211
84 gezelter 507 void Molecule::addAtom(Atom* atom) {
85 gezelter 246 if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
86 gezelter 507 atoms_.push_back(atom);
87 gezelter 246 }
88 gezelter 507 }
89 gezelter 1211
90 gezelter 507 void Molecule::addBond(Bond* bond) {
91 gezelter 246 if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
92 gezelter 507 bonds_.push_back(bond);
93 gezelter 246 }
94 gezelter 507 }
95 gezelter 1211
96 gezelter 507 void Molecule::addBend(Bend* bend) {
97 gezelter 246 if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
98 gezelter 507 bends_.push_back(bend);
99 gezelter 246 }
100 gezelter 507 }
101 gezelter 1211
102 gezelter 507 void Molecule::addTorsion(Torsion* torsion) {
103 gezelter 1211 if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
104     torsions_.end()) {
105 gezelter 507 torsions_.push_back(torsion);
106 gezelter 246 }
107 gezelter 507 }
108 gezelter 1277
109     void Molecule::addInversion(Inversion* inversion) {
110     if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
111     inversions_.end()) {
112     inversions_.push_back(inversion);
113     }
114     }
115 gezelter 1211
116 gezelter 507 void Molecule::addRigidBody(RigidBody *rb) {
117 gezelter 1211 if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
118     rigidBodies_.end()) {
119 gezelter 507 rigidBodies_.push_back(rb);
120 gezelter 246 }
121 gezelter 507 }
122 gezelter 1211
123 gezelter 507 void Molecule::addCutoffGroup(CutoffGroup* cp) {
124 gezelter 1211 if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
125     cutoffGroups_.end()) {
126 gezelter 507 cutoffGroups_.push_back(cp);
127 gezelter 1211 }
128 gezelter 507 }
129 gezelter 1211
130 gezelter 507 void Molecule::addConstraintPair(ConstraintPair* cp) {
131 gezelter 1211 if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
132     constraintPairs_.end()) {
133 gezelter 507 constraintPairs_.push_back(cp);
134 gezelter 1211 }
135 gezelter 507 }
136 gezelter 1211
137 gezelter 507 void Molecule::addConstraintElem(ConstraintElem* cp) {
138 gezelter 1211 if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
139     constraintElems_.end()) {
140 gezelter 507 constraintElems_.push_back(cp);
141 gezelter 246 }
142 gezelter 507 }
143 gezelter 1211
144 gezelter 507 void Molecule::complete() {
145 gezelter 246
146     std::set<Atom*> rigidAtoms;
147 gezelter 1782 Atom* atom;
148     AtomIterator ai;
149 gezelter 246 RigidBody* rb;
150 gezelter 1782 RigidBodyIterator rbIter;
151 gezelter 2
152 gezelter 1782 // Get list of all the atoms that are part of rigid bodies
153    
154 gezelter 246 for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
155 gezelter 507 rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
156 gezelter 246 }
157 gezelter 1211
158 gezelter 1782 // add any atom that wasn't part of a rigid body to the list of integrableObjects
159    
160 gezelter 246 for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
161 gezelter 1211
162 gezelter 1782 if (rigidAtoms.find(atom) == rigidAtoms.end()) {
163 gezelter 1211
164     // If an atom does not belong to a rigid body, it is an
165     // integrable object
166    
167 gezelter 1782 integrableObjects_.push_back(atom);
168 gezelter 507 }
169 gezelter 246 }
170 gezelter 1211
171 gezelter 1782 // then add the rigid bodies themselves to the integrableObjects
172 gezelter 2
173 tim 372 for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
174     integrableObjects_.push_back(rb);
175     }
176 gezelter 1782
177     // find the atoms that are fluctuating charges and add them to the
178     // fluctuatingCharges_ vector
179    
180     for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
181     if ( atom->isFluctuatingCharge() )
182     fluctuatingCharges_.push_back( atom );
183     }
184    
185 gezelter 507 }
186 gezelter 2
187 tim 963 RealType Molecule::getMass() {
188 gezelter 246 StuntDouble* sd;
189     std::vector<StuntDouble*>::iterator i;
190 tim 963 RealType mass = 0.0;
191 gezelter 1211
192     for (sd = beginIntegrableObject(i); sd != NULL; sd =
193     nextIntegrableObject(i)){
194 gezelter 507 mass += sd->getMass();
195 gezelter 246 }
196 gezelter 1211
197     return mass;
198 gezelter 507 }
199 gezelter 2
200 gezelter 507 Vector3d Molecule::getCom() {
201 gezelter 246 StuntDouble* sd;
202     std::vector<StuntDouble*>::iterator i;
203     Vector3d com;
204 tim 963 RealType totalMass = 0;
205     RealType mass;
206 gezelter 246
207 gezelter 1211 for (sd = beginIntegrableObject(i); sd != NULL; sd =
208     nextIntegrableObject(i)){
209 gezelter 507 mass = sd->getMass();
210     totalMass += mass;
211     com += sd->getPos() * mass;
212 gezelter 2 }
213 gezelter 1211
214 gezelter 246 com /= totalMass;
215 gezelter 2
216 gezelter 246 return com;
217 gezelter 507 }
218 gezelter 2
219 gezelter 507 void Molecule::moveCom(const Vector3d& delta) {
220 gezelter 246 StuntDouble* sd;
221     std::vector<StuntDouble*>::iterator i;
222    
223 gezelter 1211 for (sd = beginIntegrableObject(i); sd != NULL; sd =
224     nextIntegrableObject(i)){
225 gezelter 507 sd->setPos(sd->getPos() + delta);
226 gezelter 1211 }
227 gezelter 507 }
228 gezelter 2
229 gezelter 507 Vector3d Molecule::getComVel() {
230 gezelter 246 StuntDouble* sd;
231     std::vector<StuntDouble*>::iterator i;
232     Vector3d velCom;
233 tim 963 RealType totalMass = 0;
234     RealType mass;
235 gezelter 246
236 gezelter 1211 for (sd = beginIntegrableObject(i); sd != NULL; sd =
237     nextIntegrableObject(i)){
238 gezelter 507 mass = sd->getMass();
239     totalMass += mass;
240     velCom += sd->getVel() * mass;
241 gezelter 246 }
242 gezelter 1211
243 gezelter 246 velCom /= totalMass;
244 gezelter 1211
245 gezelter 246 return velCom;
246 gezelter 507 }
247 gezelter 2
248 tim 963 RealType Molecule::getPotential() {
249 gezelter 2
250 gezelter 246 Bond* bond;
251     Bend* bend;
252     Torsion* torsion;
253 gezelter 1277 Inversion* inversion;
254 gezelter 246 Molecule::BondIterator bondIter;;
255     Molecule::BendIterator bendIter;
256     Molecule::TorsionIterator torsionIter;
257 gezelter 1277 Molecule::InversionIterator inversionIter;
258 gezelter 2
259 tim 963 RealType potential = 0.0;
260 gezelter 2
261 gezelter 246 for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
262 gezelter 507 potential += bond->getPotential();
263 gezelter 246 }
264 gezelter 2
265 gezelter 246 for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
266 gezelter 507 potential += bend->getPotential();
267 gezelter 2 }
268    
269 gezelter 1211 for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
270     nextTorsion(torsionIter)) {
271 gezelter 507 potential += torsion->getPotential();
272 gezelter 2 }
273 gezelter 1211
274 gezelter 1277 for (inversion = beginInversion(inversionIter); torsion != NULL;
275     inversion = nextInversion(inversionIter)) {
276     potential += inversion->getPotential();
277     }
278    
279 gezelter 246 return potential;
280 gezelter 1211
281 gezelter 507 }
282 gezelter 1211
283 cli2 1360 void Molecule::addProperty(GenericData* genData) {
284     properties_.addProperty(genData);
285     }
286    
287     void Molecule::removeProperty(const std::string& propName) {
288     properties_.removeProperty(propName);
289     }
290    
291     void Molecule::clearProperties() {
292     properties_.clearProperties();
293     }
294    
295     std::vector<std::string> Molecule::getPropertyNames() {
296     return properties_.getPropertyNames();
297     }
298    
299     std::vector<GenericData*> Molecule::getProperties() {
300     return properties_.getProperties();
301     }
302    
303     GenericData* Molecule::getPropertyByName(const std::string& propName) {
304     return properties_.getPropertyByName(propName);
305     }
306    
307 gezelter 507 std::ostream& operator <<(std::ostream& o, Molecule& mol) {
308 gezelter 246 o << std::endl;
309     o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
310     o << mol.getNAtoms() << " atoms" << std::endl;
311     o << mol.getNBonds() << " bonds" << std::endl;
312     o << mol.getNBends() << " bends" << std::endl;
313     o << mol.getNTorsions() << " torsions" << std::endl;
314 gezelter 1277 o << mol.getNInversions() << " inversions" << std::endl;
315 gezelter 246 o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
316 gezelter 1782 o << mol.getNIntegrableObjects() << " integrable objects" << std::endl;
317     o << mol.getNCutoffGroups() << " cutoff groups" << std::endl;
318     o << mol.getNConstraintPairs() << " constraint pairs" << std::endl;
319     o << mol.getNFluctuatingCharges() << " fluctuating charges" << std::endl;
320 gezelter 246 return o;
321 gezelter 507 }
322 gezelter 1211
323 gezelter 1390 }//end namespace OpenMD

Properties

Name Value
svn:keywords Author Id Revision Date