ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/types/MoleculeStamp.cpp
Revision: 2204
Committed: Fri Apr 15 22:04:00 2005 UTC (19 years, 2 months ago) by gezelter
File size: 16010 byte(s)
Log Message:
xemacs has been drafted to perform our indentation services

File Contents

# 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 #include <stdlib.h>
43 #include <stdio.h>
44 #include <string.h>
45 #include <iostream>
46
47 #include "types/MoleculeStamp.hpp"
48
49 char MoleculeStamp::errMsg[500];
50
51 MoleculeStamp::MoleculeStamp(){
52
53 n_atoms = 0;
54 n_bonds = 0;
55 n_bends = 0;
56 n_torsions = 0;
57 n_rigidbodies = 0;
58 n_cutoffgroups = 0;
59 n_integrable = 0;
60
61 unhandled = NULL;
62 atoms = NULL;
63 bonds = NULL;
64 bends = NULL;
65 torsions = NULL;
66 rigidBodies = NULL;
67 cutoffGroups = NULL;
68
69 have_name = 0;
70 have_atoms = 0;
71 have_bonds = 0;
72 have_bends = 0;
73 have_torsions = 0;
74 have_rigidbodies = 0;
75 have_cutoffgroups = 0;
76
77 }
78
79 MoleculeStamp::~MoleculeStamp(){
80 int i;
81
82 if( unhandled != NULL) delete unhandled;
83
84 if( rigidBodies != NULL ) {
85 for( i=0; i<n_rigidbodies; i++ ) delete rigidBodies[i];
86 }
87 delete[] rigidBodies;
88
89 if( cutoffGroups != NULL ) {
90 for( i=0; i<n_cutoffgroups; i++ ) delete cutoffGroups[i];
91 }
92 delete[] cutoffGroups;
93
94 if( atoms != NULL ){
95 for( i=0; i<n_atoms; i++ ) delete atoms[i];
96 }
97 delete[] atoms;
98
99 if( bonds != NULL ){
100 for( i=0; i<n_bonds; i++ ) delete bonds[i];
101 }
102 delete[] bonds;
103
104 if( bends != NULL ){
105 for( i=0; i<n_bends; i++ ) delete bends[i];
106 }
107 delete[] bends;
108
109 if( torsions != NULL ){
110 for( i=0; i<n_torsions; i++ ) delete torsions[i];
111 }
112 delete[] torsions;
113
114
115
116
117 }
118
119 char* MoleculeStamp::assignString( char* lhs, char* rhs ){
120
121 if( !strcmp( lhs, "name" ) ){
122 strcpy( name, rhs );
123 have_name = 1;
124 }
125 else{
126 if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
127 else unhandled->add( lhs, rhs );
128 have_extras = 1;
129 }
130 return NULL;
131 }
132
133 char* MoleculeStamp::assignDouble( char* lhs, double rhs ){
134 int i;
135
136 if( !strcmp( lhs, "nAtoms" ) ){
137 n_atoms = (int)rhs;
138
139 if( have_atoms ){
140 sprintf( errMsg,
141 "MoleculeStamp error, n_atoms already declared"
142 " for molecule: %s\n",
143 name);
144 return strdup( errMsg );
145 }
146 have_atoms = 1;
147 atoms = new AtomStamp*[n_atoms];
148 for( i=0; i<n_atoms; i++ ) atoms[i] = NULL;
149 }
150
151 else if( !strcmp( lhs, "nBonds" ) ){
152 n_bonds = (int)rhs;
153
154 if( have_bonds ){
155 sprintf( errMsg,
156 "MoleculeStamp error, n_bonds already declared for"
157 " molecule: %s\n",
158 name);
159 return strdup( errMsg );
160 }
161 have_bonds = 1;
162 bonds = new BondStamp*[n_bonds];
163 for( i=0; i<n_bonds; i++ ) bonds[i] = NULL;
164 }
165
166 else if( !strcmp( lhs, "nBends" ) ){
167 n_bends = (int)rhs;
168
169 if( have_bends ){
170 sprintf( errMsg,
171 "MoleculeStamp error, n_bends already declared for"
172 " molecule: %s\n",
173 name);
174 return strdup( errMsg );
175 }
176 have_bends = 1;
177 bends = new BendStamp*[n_bends];
178 for( i=0; i<n_bends; i++ ) bends[i] = NULL;
179 }
180
181 else if( !strcmp( lhs, "nTorsions" ) ){
182 n_torsions = (int)rhs;
183
184 if( have_torsions ){
185 sprintf( errMsg,
186 "MoleculeStamp error, n_torsions already declared for"
187 " molecule: %s\n",
188 name );
189 return strdup( errMsg );
190 }
191 have_torsions = 1;
192 torsions = new TorsionStamp*[n_torsions];
193 for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
194 }
195
196 else if( !strcmp( lhs, "nRigidBodies" ) ){
197 n_rigidbodies = (int)rhs;
198
199 if( have_rigidbodies ){
200 sprintf( errMsg,
201 "MoleculeStamp error, n_rigidbodies already declared for"
202 " molecule: %s\n",
203 name );
204 return strdup( errMsg );
205 }
206 have_rigidbodies = 1;
207 rigidBodies = new RigidBodyStamp*[n_rigidbodies];
208 for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
209 }
210
211 else if( !strcmp( lhs, "nCutoffGroups" ) ){
212 n_cutoffgroups = (int)rhs;
213
214 if( have_cutoffgroups ){
215 sprintf( errMsg,
216 "MoleculeStamp error, n_cutoffgroups already declared for"
217 " molecule: %s\n",
218 name );
219 return strdup( errMsg );
220 }
221 have_cutoffgroups = 1;
222 cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
223 for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
224 }
225
226 else{
227 if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
228 else unhandled->add( lhs, rhs );
229 have_extras = 1;
230 }
231 return NULL;
232 }
233
234 char* MoleculeStamp::assignInt( char* lhs, int rhs ){
235 int i;
236
237 if( !strcmp( lhs, "nAtoms" ) ){
238 n_atoms = rhs;
239
240 if( have_atoms ){
241 sprintf( errMsg,
242 "MoleculeStamp error, n_atoms already declared for"
243 " molecule: %s\n",
244 name);
245 return strdup( errMsg );
246 }
247 have_atoms = 1;
248 atoms = new AtomStamp*[n_atoms];
249 for( i=0; i<n_atoms; i++ ) atoms[i] = NULL;
250 }
251
252 else if( !strcmp( lhs, "nBonds" ) ){
253 n_bonds = rhs;
254
255 if( have_bonds ){
256 sprintf( errMsg,
257 "MoleculeStamp error, n_bonds already declared for"
258 " molecule: %s\n",
259 name);
260 return strdup( errMsg );
261 }
262 have_bonds = 1;
263 bonds = new BondStamp*[n_bonds];
264 for( i=0; i<n_bonds; i++ ) bonds[i] = NULL;
265 }
266
267 else if( !strcmp( lhs, "nBends" ) ){
268 n_bends = rhs;
269
270 if( have_bends ){
271 sprintf( errMsg,
272 "MoleculeStamp error, n_bends already declared for"
273 " molecule: %s\n",
274 name );
275 return strdup( errMsg );
276 }
277 have_bends = 1;
278 bends = new BendStamp*[n_bends];
279 for( i=0; i<n_bends; i++ ) bends[i] = NULL;
280 }
281
282 else if( !strcmp( lhs, "nTorsions" ) ){
283 n_torsions = rhs;
284
285 if( have_torsions ){
286 sprintf( errMsg,
287 "MoleculeStamp error, n_torsions already declared for"
288 " molecule: %s\n",
289 name);
290 return strdup( errMsg );
291 }
292 have_torsions = 1;
293 torsions = new TorsionStamp*[n_torsions];
294 for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
295 }
296
297 else if( !strcmp( lhs, "nRigidBodies" ) ){
298 n_rigidbodies = rhs;
299
300 if( have_rigidbodies ){
301 sprintf( errMsg,
302 "MoleculeStamp error, n_rigidbodies already declared for"
303 " molecule: %s\n",
304 name);
305 return strdup( errMsg );
306 }
307 have_rigidbodies = 1;
308 rigidBodies = new RigidBodyStamp*[n_rigidbodies];
309 for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
310 }
311 else if( !strcmp( lhs, "nCutoffGroups" ) ){
312 n_cutoffgroups = rhs;
313
314 if( have_cutoffgroups ){
315 sprintf( errMsg,
316 "MoleculeStamp error, n_cutoffgroups already declared for"
317 " molecule: %s\n",
318 name);
319 return strdup( errMsg );
320 }
321 have_cutoffgroups = 1;
322 cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
323 for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
324 }
325 else{
326 if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
327 else unhandled->add( lhs, rhs );
328 have_extras = 1;
329 }
330 return NULL;
331 }
332
333 char* MoleculeStamp::addAtom( AtomStamp* the_atom, int atomIndex ){
334
335 if( have_atoms && atomIndex < n_atoms ) atoms[atomIndex] = the_atom;
336 else {
337 if( have_atoms ){
338 sprintf( errMsg, "MoleculeStamp error, %d out of nAtoms range",
339 atomIndex );
340 return strdup( errMsg );
341 }
342 else return strdup("MoleculeStamp error, nAtoms not given before"
343 " first atom declaration." );
344 }
345
346 return NULL;
347 }
348
349 char* MoleculeStamp::addRigidBody( RigidBodyStamp* the_rigidbody,
350 int rigidBodyIndex ){
351
352 if( have_rigidbodies && rigidBodyIndex < n_rigidbodies )
353 rigidBodies[rigidBodyIndex] = the_rigidbody;
354 else {
355 if( have_rigidbodies ){
356 sprintf( errMsg, "MoleculeStamp error, %d out of nRigidBodies range",
357 rigidBodyIndex );
358 return strdup( errMsg );
359 }
360 else return strdup("MoleculeStamp error, nRigidBodies not given before"
361 " first rigidBody declaration." );
362 }
363
364 return NULL;
365 }
366
367 char* MoleculeStamp::addCutoffGroup( CutoffGroupStamp* the_cutoffgroup,
368 int cutoffGroupIndex ){
369
370 if( have_cutoffgroups && cutoffGroupIndex < n_cutoffgroups )
371 cutoffGroups[cutoffGroupIndex] = the_cutoffgroup;
372 else {
373 if( have_cutoffgroups ){
374 sprintf( errMsg, "MoleculeStamp error, %d out of nCutoffGroups range",
375 cutoffGroupIndex );
376 return strdup( errMsg );
377 }
378 else return strdup("MoleculeStamp error, nCutoffGroups not given before"
379 " first CutoffGroup declaration." );
380 }
381
382 return NULL;
383 }
384
385 char* MoleculeStamp::addBond( BondStamp* the_bond, int bondIndex ){
386
387
388 if( have_bonds && bondIndex < n_bonds ) bonds[bondIndex] = the_bond;
389 else{
390 if( have_bonds ){
391 sprintf( errMsg, "MoleculeStamp error, %d out of nBonds range",
392 bondIndex );
393 return strdup( errMsg );
394 }
395 else return strdup("MoleculeStamp error, nBonds not given before"
396 "first bond declaration." );
397 }
398
399 return NULL;
400 }
401
402 char* MoleculeStamp::addBend( BendStamp* the_bend, int bendIndex ){
403
404
405 if( have_bends && bendIndex < n_bends ) bends[bendIndex] = the_bend;
406 else{
407 if( have_bends ){
408 sprintf( errMsg, "MoleculeStamp error, %d out of nBends range",
409 bendIndex );
410 return strdup( errMsg );
411 }
412 else return strdup("MoleculeStamp error, nBends not given before"
413 "first bend declaration." );
414 }
415
416 return NULL;
417 }
418
419 char* MoleculeStamp::addTorsion( TorsionStamp* the_torsion, int torsionIndex ){
420
421
422 if( have_torsions && torsionIndex < n_torsions )
423 torsions[torsionIndex] = the_torsion;
424 else{
425 if( have_torsions ){
426 sprintf( errMsg, "MoleculeStamp error, %d out of nTorsions range",
427 torsionIndex );
428 return strdup( errMsg );
429 }
430 else return strdup("MoleculeStamp error, nTorsions not given before"
431 "first torsion declaration." );
432 }
433
434 return NULL;
435 }
436
437 char* MoleculeStamp::checkMe( void ){
438
439 int i;
440 short int no_atom, no_rigidbody, no_cutoffgroup;
441
442 if( !have_name ) return strdup( "MoleculeStamp error. Molecule's name"
443 " was not given.\n" );
444
445 if( !have_atoms ){
446 return strdup( "MoleculeStamp error. Molecule contains no atoms." );
447 }
448
449 no_rigidbody = 0;
450 for( i=0; i<n_rigidbodies; i++ ){
451 if( rigidBodies[i] == NULL ) no_rigidbody = 1;
452 }
453
454 if( no_rigidbody ){
455 sprintf( errMsg,
456 "MoleculeStamp error. Not all of the RigidBodies were"
457 " declared in molecule \"%s\".\n", name );
458 return strdup( errMsg );
459 }
460
461 no_cutoffgroup = 0;
462 for( i=0; i<n_cutoffgroups; i++ ){
463 if( cutoffGroups[i] == NULL ) no_cutoffgroup = 1;
464 }
465
466 if( no_cutoffgroup ){
467 sprintf( errMsg,
468 "MoleculeStamp error. Not all of the CutoffGroups were"
469 " declared in molecule \"%s\".\n", name );
470 return strdup( errMsg );
471 }
472
473 no_atom = 0;
474 for( i=0; i<n_atoms; i++ ){
475 if( atoms[i] == NULL ) no_atom = 1;
476 }
477
478 if( no_atom ){
479 sprintf( errMsg,
480 "MoleculeStamp error. Not all of the atoms were"
481 " declared in molecule \"%s\".\n", name );
482 return strdup( errMsg );
483 }
484
485 n_integrable = n_atoms;
486 for (i = 0; i < n_rigidbodies; i++)
487 n_integrable = n_integrable - rigidBodies[i]->getNMembers() + 1; //rigidbody is an integrable object
488
489 if (n_integrable <= 0 || n_integrable > n_atoms) {
490 sprintf( errMsg,
491 "MoleculeStamp error. n_integrable is either <= 0 or"
492 " greater than n_atoms in molecule \"%s\".\n", name );
493 return strdup( errMsg );
494 }
495
496 return NULL;
497 }
498
499
500 //Function Name: isBondInSameRigidBody
501 //Return true is both atoms of the bond belong to the same rigid body, otherwise return false
502 bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
503 int rbA;
504 int rbB;
505 int consAtomA;
506 int consAtomB;
507
508 if (!isAtomInRigidBody(bond->getA(),rbA, consAtomA))
509 return false;
510
511 if(!isAtomInRigidBody(bond->getB(),rbB, consAtomB) )
512 return false;
513
514 if(rbB == rbA)
515 return true;
516 else
517 return false;
518 }
519
520 // Function Name: isAtomInRigidBody
521 //return false if atom does not belong to a rigid body, otherwise return true
522 bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
523 int whichRigidBody;
524 int consAtomIndex;
525
526 return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
527
528 }
529
530 // Function Name: isAtomInRigidBody
531 //return false if atom does not belong to a rigid body otherwise return true and set whichRigidBody
532 //and consAtomIndex
533 //atomIndex : the index of atom in component
534 //whichRigidBody: the index of rigidbody in component
535 //consAtomIndex: the position of joint atom apears in rigidbody's definition
536 bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
537 RigidBodyStamp* rbStamp;
538 int numRb;
539 int numAtom;
540
541 whichRigidBody = -1;
542 consAtomIndex = -1;
543
544 numRb = this->getNRigidBodies();
545
546 for(int i = 0 ; i < numRb; i++){
547 rbStamp = this->getRigidBody(i);
548 numAtom = rbStamp->getNMembers();
549 for(int j = 0; j < numAtom; j++)
550 if (rbStamp->getMember(j) == atomIndex){
551 whichRigidBody = i;
552 consAtomIndex = j;
553 return true;
554 }
555 }
556
557 return false;
558
559 }
560
561 //return the position of joint atom apears in rigidbody's definition
562 //for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
563 //actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
564 std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
565 RigidBodyStamp* rbStamp1;
566 RigidBodyStamp* rbStamp2;
567 int natomInRb1;
568 int natomInRb2;
569 int atomIndex1;
570 int atomIndex2;
571 std::vector<std::pair<int, int> > jointAtomIndexPair;
572
573 rbStamp1 = this->getRigidBody(rb1);
574 natomInRb1 =rbStamp1->getNMembers();
575
576 rbStamp2 = this->getRigidBody(rb2);
577 natomInRb2 =rbStamp2->getNMembers();
578
579 for(int i = 0; i < natomInRb1; i++){
580 atomIndex1 = rbStamp1->getMember(i);
581
582 for(int j= 0; j < natomInRb1; j++){
583 atomIndex2 = rbStamp2->getMember(j);
584
585 if(atomIndex1 == atomIndex2){
586 jointAtomIndexPair.push_back(std::make_pair(i, j));
587 break;
588 }
589
590 }//end for(j =0)
591
592 }//end for (i = 0)
593
594 return jointAtomIndexPair;
595 }