ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/visitors/AtomVisitor.cpp
Revision: 989
Committed: Sat Jun 17 17:02:33 2006 UTC (18 years, 10 months ago) by tim
File size: 17212 byte(s)
Log Message:
output the principle axis of GB

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 <cstring>
43 #include "visitors/AtomVisitor.hpp"
44 #include "primitives/DirectionalAtom.hpp"
45 #include "primitives/RigidBody.hpp"
46
47 namespace oopse {
48 void BaseAtomVisitor::visit(RigidBody *rb) {
49 //vector<Atom*> myAtoms;
50 //vector<Atom*>::iterator atomIter;
51
52 //myAtoms = rb->getAtoms();
53
54 //for(atomIter = myAtoms.begin(); atomIter != myAtoms.end(); ++atomIter)
55 // (*atomIter)->accept(this);
56 }
57
58 void BaseAtomVisitor::setVisited(Atom *atom) {
59 GenericData *data;
60 data = atom->getPropertyByName("VISITED");
61
62 //if visited property is not existed, add it as new property
63 if (data == NULL) {
64 data = new GenericData();
65 data->setID("VISITED");
66 atom->addProperty(data);
67 }
68 }
69
70 bool BaseAtomVisitor::isVisited(Atom *atom) {
71 GenericData *data;
72 data = atom->getPropertyByName("VISITED");
73 return data == NULL ? false : true;
74 }
75
76 bool SSDAtomVisitor::isSSDAtom(const std::string&atomType) {
77 std::set<std::string>::iterator strIter;
78 strIter = ssdAtomType.find(atomType);
79 return strIter != ssdAtomType.end() ? true : false;
80 }
81
82 void SSDAtomVisitor::visit(DirectionalAtom *datom) {
83 std::vector<AtomInfo*>atoms;
84
85 //we need to convert SSD into 4 differnet atoms
86 //one oxygen atom, two hydrogen atoms and one pseudo atom which is the center of the mass
87 //of the water with a dipole moment
88 Vector3d h1(0.0, -0.75695, 0.5206);
89 Vector3d h2(0.0, 0.75695, 0.5206);
90 Vector3d ox(0.0, 0.0, -0.0654);
91 Vector3d u(0, 0, 1);
92 RotMat3x3d rotMatrix;
93 RotMat3x3d rotTrans;
94 AtomInfo * atomInfo;
95 Vector3d pos;
96 Vector3d newVec;
97 Quat4d q;
98 AtomData * atomData;
99 GenericData *data;
100 bool haveAtomData;
101
102 //if atom is not SSD atom, just skip it
103 if (!isSSDAtom(datom->getType()))
104 return;
105
106 data = datom->getPropertyByName("ATOMDATA");
107
108 if (data != NULL) {
109 atomData = dynamic_cast<AtomData *>(data);
110
111 if (atomData == NULL) {
112 std::cerr << "can not get Atom Data from " << datom->getType() << std::endl;
113 atomData = new AtomData;
114 haveAtomData = false;
115 } else
116 haveAtomData = true;
117 } else {
118 atomData = new AtomData;
119 haveAtomData = false;
120 }
121
122 pos = datom->getPos();
123 q = datom->getQ();
124 rotMatrix = datom->getA();
125
126 // We need A^T to convert from body-fixed to space-fixed:
127 //transposeMat3(rotMatrix, rotTrans);
128 rotTrans = rotMatrix.transpose();
129
130 //center of mass of the water molecule
131 //matVecMul3(rotTrans, u, newVec);
132 newVec = rotTrans * u;
133
134 atomInfo = new AtomInfo;
135 atomInfo->atomTypeName = "X";
136 atomInfo->pos[0] = pos[0];
137 atomInfo->pos[1] = pos[1];
138 atomInfo->pos[2] = pos[2];
139 atomInfo->dipole[0] = newVec[0];
140 atomInfo->dipole[1] = newVec[1];
141 atomInfo->dipole[2] = newVec[2];
142
143 atomData->addAtomInfo(atomInfo);
144
145 //oxygen
146 //matVecMul3(rotTrans, ox, newVec);
147 newVec = rotTrans * ox;
148
149 atomInfo = new AtomInfo;
150 atomInfo->atomTypeName = "O";
151 atomInfo->pos[0] = pos[0] + newVec[0];
152 atomInfo->pos[1] = pos[1] + newVec[1];
153 atomInfo->pos[2] = pos[2] + newVec[2];
154 atomInfo->dipole[0] = 0.0;
155 atomInfo->dipole[1] = 0.0;
156 atomInfo->dipole[2] = 0.0;
157 atomData->addAtomInfo(atomInfo);
158
159 //hydrogen1
160 //matVecMul3(rotTrans, h1, newVec);
161 newVec = rotTrans * h1;
162 atomInfo = new AtomInfo;
163 atomInfo->atomTypeName = "H";
164 atomInfo->pos[0] = pos[0] + newVec[0];
165 atomInfo->pos[1] = pos[1] + newVec[1];
166 atomInfo->pos[2] = pos[2] + newVec[2];
167 atomInfo->dipole[0] = 0.0;
168 atomInfo->dipole[1] = 0.0;
169 atomInfo->dipole[2] = 0.0;
170 atomData->addAtomInfo(atomInfo);
171
172 //hydrogen2
173 //matVecMul3(rotTrans, h2, newVec);
174 newVec = rotTrans * h2;
175 atomInfo = new AtomInfo;
176 atomInfo->atomTypeName = "H";
177 atomInfo->pos[0] = pos[0] + newVec[0];
178 atomInfo->pos[1] = pos[1] + newVec[1];
179 atomInfo->pos[2] = pos[2] + newVec[2];
180 atomInfo->dipole[0] = 0.0;
181 atomInfo->dipole[1] = 0.0;
182 atomInfo->dipole[2] = 0.0;
183 atomData->addAtomInfo(atomInfo);
184
185 //add atom data into atom's property
186
187 if (!haveAtomData) {
188 atomData->setID("ATOMDATA");
189 datom->addProperty(atomData);
190 }
191
192 setVisited(datom);
193 }
194
195 const std::string SSDAtomVisitor::toString() {
196 char buffer[65535];
197 std::string result;
198
199 sprintf(buffer,
200 "------------------------------------------------------------------\n");
201 result += buffer;
202
203 sprintf(buffer, "Visitor name: %s\n", visitorName.c_str());
204 result += buffer;
205
206 sprintf(buffer,
207 "Visitor Description: Convert SSD into 4 different atoms\n");
208 result += buffer;
209
210 sprintf(buffer,
211 "------------------------------------------------------------------\n");
212 result += buffer;
213
214 return result;
215 }
216
217 bool LinearAtomVisitor::isLinearAtom(const std::string& atomType){
218 std::set<std::string>::iterator strIter;
219 strIter = linearAtomType.find(atomType);
220
221 return strIter != linearAtomType.end() ? true : false;
222 }
223
224 void LinearAtomVisitor::addGayBerneAtomType(const std::string& atomType){
225 linearAtomType.insert(atomType);
226 }
227
228 void LinearAtomVisitor::visit(DirectionalAtom* datom){
229 std::vector<AtomInfo*> atoms;
230 //we need to convert linear into 4 different atoms
231 Vector3d c1(0.0, 0.0, -1.8);
232 Vector3d c2(0.0, 0.0, -0.6);
233 Vector3d c3(0.0, 0.0, 0.6);
234 Vector3d c4(0.0, 0.0, 1.8);
235 RotMat3x3d rotMatrix;
236 RotMat3x3d rotTrans;
237 AtomInfo* atomInfo;
238 Vector3d pos;
239 Vector3d newVec;
240 Quat4d q;
241 AtomData* atomData;
242 GenericData* data;
243 bool haveAtomData;
244 AtomType* atomType;
245 //if atom is not SSD atom, just skip it
246 if(!isLinearAtom(datom->getType()) || !datom->getAtomType()->isGayBerne())
247 return;
248
249 //setup GayBerne type in fortran side
250 data = datom->getAtomType()->getPropertyByName("GayBerne");
251 if (data != NULL) {
252 GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
253
254 if (gayBerneData != NULL) {
255 GayBerneParam gayBerneParam = gayBerneData->getData();
256
257 double halfLen = gayBerneParam.GB_sigma * gayBerneParam.GB_l2b_ratio/2.0;
258 c1[2] = -halfLen;
259 c2[2] = -halfLen /2;
260 c3[2] = halfLen/2;
261 c4[2] = halfLen;
262
263 }
264
265 else {
266 sprintf( painCave.errMsg,
267 "Can not cast GenericData to GayBerneParam\n");
268 painCave.severity = OOPSE_ERROR;
269 painCave.isFatal = 1;
270 simError();
271 }
272 }
273
274
275 data = datom->getPropertyByName("ATOMDATA");
276 if(data != NULL){
277 atomData = dynamic_cast<AtomData*>(data);
278 if(atomData == NULL){
279 std::cerr << "can not get Atom Data from " << datom->getType() << std::endl;
280 atomData = new AtomData;
281 haveAtomData = false;
282 } else {
283 haveAtomData = true;
284 }
285 } else {
286 atomData = new AtomData;
287 haveAtomData = false;
288 }
289
290
291 pos = datom->getPos();
292 q = datom->getQ();
293 rotMatrix = datom->getA();
294
295 // We need A^T to convert from body-fixed to space-fixed:
296 rotTrans = rotMatrix.transpose();
297
298 newVec = rotTrans * c1;
299 atomInfo = new AtomInfo;
300 atomInfo->atomTypeName = "C";
301 atomInfo->pos[0] = pos[0] + newVec[0];
302 atomInfo->pos[1] = pos[1] + newVec[1];
303 atomInfo->pos[2] = pos[2] + newVec[2];
304 atomInfo->dipole[0] = 0.0;
305 atomInfo->dipole[1] = 0.0;
306 atomInfo->dipole[2] = 0.0;
307 atomData->addAtomInfo(atomInfo);
308
309 newVec = rotTrans * c2;
310 atomInfo = new AtomInfo;
311 atomInfo->atomTypeName = "C";
312 atomInfo->pos[0] = pos[0] + newVec[0];
313 atomInfo->pos[1] = pos[1] + newVec[1];
314 atomInfo->pos[2] = pos[2] + newVec[2];
315 atomInfo->dipole[0] = 0.0;
316 atomInfo->dipole[1] = 0.0;
317 atomInfo->dipole[2] = 0.0;
318 atomData->addAtomInfo(atomInfo);
319
320 newVec = rotTrans * c3;
321 atomInfo = new AtomInfo;
322 atomInfo->atomTypeName = "C";
323 atomInfo->pos[0] = pos[0] + newVec[0];
324 atomInfo->pos[1] = pos[1] + newVec[1];
325 atomInfo->pos[2] = pos[2] + newVec[2];
326 atomInfo->dipole[0] = 0.0;
327 atomInfo->dipole[1] = 0.0;
328 atomInfo->dipole[2] = 0.0;
329 atomData->addAtomInfo(atomInfo);
330
331 newVec = rotTrans * c4;
332 atomInfo = new AtomInfo;
333 atomInfo->atomTypeName = "C";
334 atomInfo->pos[0] = pos[0] + newVec[0];
335 atomInfo->pos[1] = pos[1] + newVec[1];
336 atomInfo->pos[2] = pos[2] + newVec[2];
337 atomInfo->dipole[0] = 0.0;
338 atomInfo->dipole[1] = 0.0;
339 atomInfo->dipole[2] = 0.0;
340 atomData->addAtomInfo(atomInfo);
341
342 //add atom data into atom's property
343
344 if(!haveAtomData){
345 atomData->setID("ATOMDATA");
346 datom->addProperty(atomData);
347 }
348
349 setVisited(datom);
350
351 }
352
353 const std::string LinearAtomVisitor::toString(){
354 char buffer[65535];
355 std::string result;
356
357 sprintf(buffer ,"------------------------------------------------------------------\n");
358 result += buffer;
359
360 sprintf(buffer ,"Visitor name: %s\n", visitorName.c_str());
361 result += buffer;
362
363 sprintf(buffer , "Visitor Description: Convert linear into 4 different atoms\n");
364 result += buffer;
365
366 sprintf(buffer ,"------------------------------------------------------------------\n");
367 result += buffer;
368
369 return result;
370 }
371
372 bool GBLipidAtomVisitor::isGBLipidAtom(const std::string& atomType){
373 std::set<std::string>::iterator strIter;
374 strIter = GBLipidAtomType.find(atomType);
375
376 return strIter != GBLipidAtomType.end() ? true : false;
377 }
378
379 void GBLipidAtomVisitor::visit(DirectionalAtom* datom){
380 std::vector<AtomInfo*> atoms;
381 //we need to convert linear into 4 different atoms
382 Vector3d c1(0.0, 0.0, -6.25);
383 Vector3d c2(0.0, 0.0, -2.1);
384 Vector3d c3(0.0, 0.0, 2.1);
385 Vector3d c4(0.0, 0.0, 6.25);
386 RotMat3x3d rotMatrix;
387 RotMat3x3d rotTrans;
388 AtomInfo* atomInfo;
389 Vector3d pos;
390 Vector3d newVec;
391 Quat4d q;
392 AtomData* atomData;
393 GenericData* data;
394 bool haveAtomData;
395
396 //if atom is not GBlipid atom, just skip it
397 if(!isGBLipidAtom(datom->getType()))
398 return;
399
400 data = datom->getPropertyByName("ATOMDATA");
401 if(data != NULL){
402 atomData = dynamic_cast<AtomData*>(data);
403 if(atomData == NULL){
404 std::cerr << "can not get Atom Data from " << datom->getType() << std::endl;
405 atomData = new AtomData;
406 haveAtomData = false;
407 } else {
408 haveAtomData = true;
409 }
410 } else {
411 atomData = new AtomData;
412 haveAtomData = false;
413 }
414
415
416 pos = datom->getPos();
417 q = datom->getQ();
418 rotMatrix = datom->getA();
419
420 // We need A^T to convert from body-fixed to space-fixed:
421 rotTrans = rotMatrix.transpose();
422
423 newVec = rotTrans * c1;
424 atomInfo = new AtomInfo;
425 atomInfo->atomTypeName = "K";
426 atomInfo->pos[0] = pos[0] + newVec[0];
427 atomInfo->pos[1] = pos[1] + newVec[1];
428 atomInfo->pos[2] = pos[2] + newVec[2];
429 atomInfo->dipole[0] = 0.0;
430 atomInfo->dipole[1] = 0.0;
431 atomInfo->dipole[2] = 0.0;
432 atomData->addAtomInfo(atomInfo);
433
434 newVec = rotTrans * c2;
435 atomInfo = new AtomInfo;
436 atomInfo->atomTypeName = "K";
437 atomInfo->pos[0] = pos[0] + newVec[0];
438 atomInfo->pos[1] = pos[1] + newVec[1];
439 atomInfo->pos[2] = pos[2] + newVec[2];
440 atomInfo->dipole[0] = 0.0;
441 atomInfo->dipole[1] = 0.0;
442 atomInfo->dipole[2] = 0.0;
443 atomData->addAtomInfo(atomInfo);
444
445 newVec = rotTrans * c3;
446 atomInfo = new AtomInfo;
447 atomInfo->atomTypeName = "K";
448 atomInfo->pos[0] = pos[0] + newVec[0];
449 atomInfo->pos[1] = pos[1] + newVec[1];
450 atomInfo->pos[2] = pos[2] + newVec[2];
451 atomInfo->dipole[0] = 0.0;
452 atomInfo->dipole[1] = 0.0;
453 atomInfo->dipole[2] = 0.0;
454 atomData->addAtomInfo(atomInfo);
455
456 newVec = rotTrans * c4;
457 atomInfo = new AtomInfo;
458 atomInfo->atomTypeName = "K";
459 atomInfo->pos[0] = pos[0] + newVec[0];
460 atomInfo->pos[1] = pos[1] + newVec[1];
461 atomInfo->pos[2] = pos[2] + newVec[2];
462 atomInfo->dipole[0] = 0.0;
463 atomInfo->dipole[1] = 0.0;
464 atomInfo->dipole[2] = 0.0;
465 atomData->addAtomInfo(atomInfo);
466
467 //add atom data into atom's property
468
469 if(!haveAtomData){
470 atomData->setID("ATOMDATA");
471 datom->addProperty(atomData);
472 }
473
474 setVisited(datom);
475
476 }
477
478 const std::string GBLipidAtomVisitor::toString(){
479 char buffer[65535];
480 std::string result;
481
482 sprintf(buffer ,"------------------------------------------------------------------\n");
483 result += buffer;
484
485 sprintf(buffer ,"Visitor name: %s\n", visitorName.c_str());
486 result += buffer;
487
488 sprintf(buffer , "Visitor Description: Convert GBlipid into 4 different K atoms\n");
489 result += buffer;
490
491 sprintf(buffer ,"------------------------------------------------------------------\n");
492 result += buffer;
493
494 return result;
495 }
496
497 //----------------------------------------------------------------------------//
498
499 void DefaultAtomVisitor::visit(Atom *atom) {
500 AtomData *atomData;
501 AtomInfo *atomInfo;
502 Vector3d pos;
503
504 if (isVisited(atom))
505 return;
506
507 atomInfo = new AtomInfo;
508
509 atomData = new AtomData;
510 atomData->setID("ATOMDATA");
511
512 pos = atom->getPos();
513 atomInfo->atomTypeName = atom->getType();
514 atomInfo->pos[0] = pos[0];
515 atomInfo->pos[1] = pos[1];
516 atomInfo->pos[2] = pos[2];
517 atomInfo->dipole[0] = 0.0;
518 atomInfo->dipole[1] = 0.0;
519 atomInfo->dipole[2] = 0.0;
520
521 atomData->addAtomInfo(atomInfo);
522
523 atom->addProperty(atomData);
524
525 setVisited(atom);
526 }
527
528 void DefaultAtomVisitor::visit(DirectionalAtom *datom) {
529 AtomData *atomData;
530 AtomInfo *atomInfo;
531 Vector3d pos;
532 Vector3d u;
533
534 if (isVisited(datom))
535 return;
536
537 pos = datom->getPos();
538 if (datom->getAtomType()->isGayBerne()) {
539 u = datom->getA().transpose()*V3Z;
540 } else if (datom->getAtomType()->isMultipole()) {
541 u = datom->getElectroFrame().getColumn(2);
542 }
543 atomData = new AtomData;
544 atomData->setID("ATOMDATA");
545 atomInfo = new AtomInfo;
546
547 atomInfo->atomTypeName = datom->getType();
548 atomInfo->pos[0] = pos[0];
549 atomInfo->pos[1] = pos[1];
550 atomInfo->pos[2] = pos[2];
551 atomInfo->dipole[0] = u[0];
552 atomInfo->dipole[1] = u[1];
553 atomInfo->dipole[2] = u[2];
554
555 atomData->addAtomInfo(atomInfo);
556
557 datom->addProperty(atomData);
558
559 setVisited(datom);
560 }
561
562 const std::string DefaultAtomVisitor::toString() {
563 char buffer[65535];
564 std::string result;
565
566 sprintf(buffer,
567 "------------------------------------------------------------------\n");
568 result += buffer;
569
570 sprintf(buffer, "Visitor name: %s\n", visitorName.c_str());
571 result += buffer;
572
573 sprintf(buffer,
574 "Visitor Description: copy atom infomation into atom data\n");
575 result += buffer;
576
577 sprintf(buffer,
578 "------------------------------------------------------------------\n");
579 result += buffer;
580
581 return result;
582 }
583 } //namespace oopse

Properties

Name Value
svn:executable *