src/primitives/Molecule.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
/**
 * @file Molecule.cpp
 * @author    tlin
 * @date  10/28/2004
 * @version 1.0
 */ 

#include <algorithm>
#include <set>

#include "primitives/Molecule.hpp"
#include "utils/MemoryUtils.hpp"
#include "utils/simError.h"

namespace oopse {
  Molecule::Molecule(int stampId, int globalIndex, const std::string& molName) 
    : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
  }
  
  Molecule::~Molecule() {
    
    MemoryUtils::deletePointers(atoms_);
    MemoryUtils::deletePointers(bonds_);
    MemoryUtils::deletePointers(bends_);
    MemoryUtils::deletePointers(torsions_);
    MemoryUtils::deletePointers(inversions_);
    MemoryUtils::deletePointers(rigidBodies_);
    MemoryUtils::deletePointers(cutoffGroups_);
    MemoryUtils::deletePointers(constraintPairs_);
    MemoryUtils::deletePointers(constraintElems_);
    // integrableObjects_ don't own the objects
    integrableObjects_.clear();
    
  }
  
  void Molecule::addAtom(Atom* atom) {
    if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
      atoms_.push_back(atom);
    }
  }
  
  void Molecule::addBond(Bond* bond) {
    if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
      bonds_.push_back(bond);
    }
  }
  
  void Molecule::addBend(Bend* bend) {
    if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
      bends_.push_back(bend);
    }
  }
  
  void Molecule::addTorsion(Torsion* torsion) {
    if (std::find(torsions_.begin(), torsions_.end(), torsion) == 
        torsions_.end()) {
      torsions_.push_back(torsion);
    }
  }

  void Molecule::addInversion(Inversion* inversion) {
    if (std::find(inversions_.begin(), inversions_.end(), inversion) == 
        inversions_.end()) {
      inversions_.push_back(inversion);
    }
  }
  
  void Molecule::addRigidBody(RigidBody *rb) {
    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == 
        rigidBodies_.end()) {
      rigidBodies_.push_back(rb);
    }
  }
  
  void Molecule::addCutoffGroup(CutoffGroup* cp) {
    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == 
        cutoffGroups_.end()) {
      cutoffGroups_.push_back(cp);
    }    
  }
  
  void Molecule::addConstraintPair(ConstraintPair* cp) {
    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == 
        constraintPairs_.end()) {
      constraintPairs_.push_back(cp);
    }    
  }
  
  void Molecule::addConstraintElem(ConstraintElem* cp) {
    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == 
        constraintElems_.end()) {
      constraintElems_.push_back(cp);
    }
  }
  
  void Molecule::complete() {
    
    std::set<Atom*> rigidAtoms;
    RigidBody* rb;
    std::vector<RigidBody*>::iterator rbIter;

    
    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
      rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
    }
    
    Atom* atom;
    AtomIterator ai;
    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
      
      if (rigidAtoms.find(*ai) == rigidAtoms.end()) {

        // If an atom does not belong to a rigid body, it is an
        // integrable object

        integrableObjects_.push_back(*ai);
      }
    }
    
    //find all free atoms (which do not belong to rigid bodies)  
    // performs the "difference" operation from set theory, the output
    // range contains a copy of every element that is contained in
    // [allAtoms.begin(), allAtoms.end()) and not contained in
    // [rigidAtoms.begin(), rigidAtoms.end()).
    //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
    //                        std::back_inserter(integrableObjects_));

    //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
    //    //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
    //    sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
    //
    //    painCave.isFatal = 1;
    //    simError();        
    //}
    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
      integrableObjects_.push_back(rb);
    } 
    //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
  }

  RealType Molecule::getMass() {
    StuntDouble* sd;
    std::vector<StuntDouble*>::iterator i;
    RealType mass = 0.0;
    
    for (sd = beginIntegrableObject(i); sd != NULL; sd = 
           nextIntegrableObject(i)){
      mass += sd->getMass();
    }
    
    return mass;    
  }

  Vector3d Molecule::getCom() {
    StuntDouble* sd;
    std::vector<StuntDouble*>::iterator i;
    Vector3d com;
    RealType totalMass = 0;
    RealType mass;
    
    for (sd = beginIntegrableObject(i); sd != NULL; sd = 
           nextIntegrableObject(i)){
      mass = sd->getMass();
      totalMass += mass;
      com += sd->getPos() * mass;    
    }
    
    com /= totalMass;

    return com;
  }

  void Molecule::moveCom(const Vector3d& delta) {
    StuntDouble* sd;
    std::vector<StuntDouble*>::iterator i;
    
    for (sd = beginIntegrableObject(i); sd != NULL; sd = 
           nextIntegrableObject(i)){
      sd->setPos(sd->getPos() + delta);
    }    
  }

  Vector3d Molecule::getComVel() {
    StuntDouble* sd;
    std::vector<StuntDouble*>::iterator i;
    Vector3d velCom;
    RealType totalMass = 0;
    RealType mass;
    
    for (sd = beginIntegrableObject(i); sd != NULL; sd = 
           nextIntegrableObject(i)){
      mass = sd->getMass();
      totalMass += mass;
      velCom += sd->getVel() * mass;    
    }
    
    velCom /= totalMass;
    
    return velCom;
  }

  RealType Molecule::getPotential() {

    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    Inversion* inversion;
    Molecule::BondIterator bondIter;;
    Molecule::BendIterator  bendIter;
    Molecule::TorsionIterator  torsionIter;
    Molecule::InversionIterator  inversionIter;

    RealType potential = 0.0;

    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
      potential += bond->getPotential();
    }

    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
      potential += bend->getPotential();
    }

    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = 
           nextTorsion(torsionIter)) {
      potential += torsion->getPotential();
    }
    
    for (inversion = beginInversion(inversionIter); torsion != NULL; 
         inversion =  nextInversion(inversionIter)) {
      potential += inversion->getPotential();
    }
    
    return potential;
    
  }
  
  void Molecule::addProperty(GenericData* genData) {
    properties_.addProperty(genData);  
  }

  void Molecule::removeProperty(const std::string& propName) {
    properties_.removeProperty(propName);  
  }

  void Molecule::clearProperties() {
    properties_.clearProperties(); 
  }

  std::vector<std::string> Molecule::getPropertyNames() {
    return properties_.getPropertyNames();  
  }
      
  std::vector<GenericData*> Molecule::getProperties() { 
    return properties_.getProperties(); 
  }

  GenericData* Molecule::getPropertyByName(const std::string& propName) {
    return properties_.getPropertyByName(propName); 
  }


  std::ostream& operator <<(std::ostream& o, Molecule& mol) {
    o << std::endl;
    o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
    o << mol.getNAtoms() << " atoms" << std::endl;
    o << mol.getNBonds() << " bonds" << std::endl;
    o << mol.getNBends() << " bends" << std::endl;
    o << mol.getNTorsions() << " torsions" << std::endl;
    o << mol.getNInversions() << " inversions" << std::endl;
    o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
    o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
    return o;
  }
  
}//end namespace oopse
Revision:	3520
Committed:	Mon Sep 7 16:31:51 2009 UTC (15 years, 7 months ago) by cli2
File size:	10102 byte(s)
Log Message:	Added new restraint infrastructure Added MolecularRestraints Added ObjectRestraints Added RestraintStamp Updated thermodynamic integration to use ObjectRestraints Added Quaternion mathematics for twist swing decompositions Significantly updated RestWriter and RestReader to use dump-like files Added selections for x, y, and z coordinates of atoms Removed monolithic Restraints class Fixed a few bugs in gradients of Euler angles in DirectionalAtom and RigidBody Added some rotational capabilities to prinicpalAxisCalculator
#	Content
1	/*
2	* 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	* 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
19	* notice, this list of conditions and the following disclaimer.
20	*
21	* 3. Redistributions in binary form must reproduce the above copyright
22	* notice, this list of conditions and the following disclaimer in the
23	* documentation and/or other materials provided with the
24	* distribution.
25	*
26	* This software is provided "AS IS," without a warranty of any
27	* kind. All express or implied conditions, representations and
28	* warranties, including any implied warranty of merchantability,
29	* fitness for a particular purpose or non-infringement, are hereby
30	* excluded. The University of Notre Dame and its licensors shall not
31	* be liable for any damages suffered by licensee as a result of
32	* using, modifying or distributing the software or its
33	* derivatives. In no event will the University of Notre Dame or its
34	* licensors be liable for any lost revenue, profit or data, or for
35	* direct, indirect, special, consequential, incidental or punitive
36	* damages, however caused and regardless of the theory of liability,
37	* arising out of the use of or inability to use software, even if the
38	* University of Notre Dame has been advised of the possibility of
39	* such damages.
40	*/
41
42	/**
43	* @file Molecule.cpp
44	* @author tlin
45	* @date 10/28/2004
46	* @version 1.0
47	*/
48
49	#include <algorithm>
50	#include <set>
51
52	#include "primitives/Molecule.hpp"
53	#include "utils/MemoryUtils.hpp"
54	#include "utils/simError.h"
55
56	namespace oopse {
57	Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
58	: stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
59	}
60
61	Molecule::~Molecule() {
62
63	MemoryUtils::deletePointers(atoms_);
64	MemoryUtils::deletePointers(bonds_);
65	MemoryUtils::deletePointers(bends_);
66	MemoryUtils::deletePointers(torsions_);
67	MemoryUtils::deletePointers(inversions_);
68	MemoryUtils::deletePointers(rigidBodies_);
69	MemoryUtils::deletePointers(cutoffGroups_);
70	MemoryUtils::deletePointers(constraintPairs_);
71	MemoryUtils::deletePointers(constraintElems_);
72	// integrableObjects_ don't own the objects
73	integrableObjects_.clear();
74
75	}
76
77	void Molecule::addAtom(Atom* atom) {
78	if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
79	atoms_.push_back(atom);
80	}
81	}
82
83	void Molecule::addBond(Bond* bond) {
84	if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
85	bonds_.push_back(bond);
86	}
87	}
88
89	void Molecule::addBend(Bend* bend) {
90	if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
91	bends_.push_back(bend);
92	}
93	}
94
95	void Molecule::addTorsion(Torsion* torsion) {
96	if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
97	torsions_.end()) {
98	torsions_.push_back(torsion);
99	}
100	}
101
102	void Molecule::addInversion(Inversion* inversion) {
103	if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
104	inversions_.end()) {
105	inversions_.push_back(inversion);
106	}
107	}
108
109	void Molecule::addRigidBody(RigidBody *rb) {
110	if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
111	rigidBodies_.end()) {
112	rigidBodies_.push_back(rb);
113	}
114	}
115
116	void Molecule::addCutoffGroup(CutoffGroup* cp) {
117	if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
118	cutoffGroups_.end()) {
119	cutoffGroups_.push_back(cp);
120	}
121	}
122
123	void Molecule::addConstraintPair(ConstraintPair* cp) {
124	if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
125	constraintPairs_.end()) {
126	constraintPairs_.push_back(cp);
127	}
128	}
129
130	void Molecule::addConstraintElem(ConstraintElem* cp) {
131	if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
132	constraintElems_.end()) {
133	constraintElems_.push_back(cp);
134	}
135	}
136
137	void Molecule::complete() {
138
139	std::set<Atom*> rigidAtoms;
140	RigidBody* rb;
141	std::vector<RigidBody*>::iterator rbIter;
142
143
144	for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
145	rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
146	}
147
148	Atom* atom;
149	AtomIterator ai;
150	for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
151
152	if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
153
154	// If an atom does not belong to a rigid body, it is an
155	// integrable object
156
157	integrableObjects_.push_back(*ai);
158	}
159	}
160
161	//find all free atoms (which do not belong to rigid bodies)
162	// performs the "difference" operation from set theory, the output
163	// range contains a copy of every element that is contained in
164	// [allAtoms.begin(), allAtoms.end()) and not contained in
165	// [rigidAtoms.begin(), rigidAtoms.end()).
166	//std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
167	// std::back_inserter(integrableObjects_));
168
169	//if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
170	// //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
171	// sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
172	//
173	// painCave.isFatal = 1;
174	// simError();
175	//}
176	for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
177	integrableObjects_.push_back(rb);
178	}
179	//integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
180	}
181
182	RealType Molecule::getMass() {
183	StuntDouble* sd;
184	std::vector<StuntDouble*>::iterator i;
185	RealType mass = 0.0;
186
187	for (sd = beginIntegrableObject(i); sd != NULL; sd =
188	nextIntegrableObject(i)){
189	mass += sd->getMass();
190	}
191
192	return mass;
193	}
194
195	Vector3d Molecule::getCom() {
196	StuntDouble* sd;
197	std::vector<StuntDouble*>::iterator i;
198	Vector3d com;
199	RealType totalMass = 0;
200	RealType mass;
201
202	for (sd = beginIntegrableObject(i); sd != NULL; sd =
203	nextIntegrableObject(i)){
204	mass = sd->getMass();
205	totalMass += mass;
206	com += sd->getPos() * mass;
207	}
208
209	com /= totalMass;
210
211	return com;
212	}
213
214	void Molecule::moveCom(const Vector3d& delta) {
215	StuntDouble* sd;
216	std::vector<StuntDouble*>::iterator i;
217
218	for (sd = beginIntegrableObject(i); sd != NULL; sd =
219	nextIntegrableObject(i)){
220	sd->setPos(sd->getPos() + delta);
221	}
222	}
223
224	Vector3d Molecule::getComVel() {
225	StuntDouble* sd;
226	std::vector<StuntDouble*>::iterator i;
227	Vector3d velCom;
228	RealType totalMass = 0;
229	RealType mass;
230
231	for (sd = beginIntegrableObject(i); sd != NULL; sd =
232	nextIntegrableObject(i)){
233	mass = sd->getMass();
234	totalMass += mass;
235	velCom += sd->getVel() * mass;
236	}
237
238	velCom /= totalMass;
239
240	return velCom;
241	}
242
243	RealType Molecule::getPotential() {
244
245	Bond* bond;
246	Bend* bend;
247	Torsion* torsion;
248	Inversion* inversion;
249	Molecule::BondIterator bondIter;;
250	Molecule::BendIterator bendIter;
251	Molecule::TorsionIterator torsionIter;
252	Molecule::InversionIterator inversionIter;
253
254	RealType potential = 0.0;
255
256	for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
257	potential += bond->getPotential();
258	}
259
260	for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
261	potential += bend->getPotential();
262	}
263
264	for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
265	nextTorsion(torsionIter)) {
266	potential += torsion->getPotential();
267	}
268
269	for (inversion = beginInversion(inversionIter); torsion != NULL;
270	inversion = nextInversion(inversionIter)) {
271	potential += inversion->getPotential();
272	}
273
274	return potential;
275
276	}
277
278	void Molecule::addProperty(GenericData* genData) {
279	properties_.addProperty(genData);
280	}
281
282	void Molecule::removeProperty(const std::string& propName) {
283	properties_.removeProperty(propName);
284	}
285
286	void Molecule::clearProperties() {
287	properties_.clearProperties();
288	}
289
290	std::vector<std::string> Molecule::getPropertyNames() {
291	return properties_.getPropertyNames();
292	}
293
294	std::vector<GenericData*> Molecule::getProperties() {
295	return properties_.getProperties();
296	}
297
298	GenericData* Molecule::getPropertyByName(const std::string& propName) {
299	return properties_.getPropertyByName(propName);
300	}
301
302
303
304
305	std::ostream& operator <<(std::ostream& o, Molecule& mol) {
306	o << std::endl;
307	o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
308	o << mol.getNAtoms() << " atoms" << std::endl;
309	o << mol.getNBonds() << " bonds" << std::endl;
310	o << mol.getNBends() << " bends" << std::endl;
311	o << mol.getNTorsions() << " torsions" << std::endl;
312	o << mol.getNInversions() << " inversions" << std::endl;
313	o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
314	o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
315	o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
316	o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
317	return o;
318	}
319
320	}//end namespace oopse