ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/visitors/AtomVisitor.cpp
Revision: 992
Committed: Wed Jun 21 18:28:27 2006 UTC (18 years, 10 months ago) by chrisfen
File size: 17257 byte(s)
Log Message:
fixed a build problem...  GB_12b_ratio no longer exists (replaced with BG_l).

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 double halfLen = gayBerneParam.GB_l/2.0;
259 c1[2] = -halfLen;
260 c2[2] = -halfLen /2;
261 c3[2] = halfLen/2;
262 c4[2] = halfLen;
263
264 }
265
266 else {
267 sprintf( painCave.errMsg,
268 "Can not cast GenericData to GayBerneParam\n");
269 painCave.severity = OOPSE_ERROR;
270 painCave.isFatal = 1;
271 simError();
272 }
273 }
274
275
276 data = datom->getPropertyByName("ATOMDATA");
277 if(data != NULL){
278 atomData = dynamic_cast<AtomData*>(data);
279 if(atomData == NULL){
280 std::cerr << "can not get Atom Data from " << datom->getType() << std::endl;
281 atomData = new AtomData;
282 haveAtomData = false;
283 } else {
284 haveAtomData = true;
285 }
286 } else {
287 atomData = new AtomData;
288 haveAtomData = false;
289 }
290
291
292 pos = datom->getPos();
293 q = datom->getQ();
294 rotMatrix = datom->getA();
295
296 // We need A^T to convert from body-fixed to space-fixed:
297 rotTrans = rotMatrix.transpose();
298
299 newVec = rotTrans * c1;
300 atomInfo = new AtomInfo;
301 atomInfo->atomTypeName = "C";
302 atomInfo->pos[0] = pos[0] + newVec[0];
303 atomInfo->pos[1] = pos[1] + newVec[1];
304 atomInfo->pos[2] = pos[2] + newVec[2];
305 atomInfo->dipole[0] = 0.0;
306 atomInfo->dipole[1] = 0.0;
307 atomInfo->dipole[2] = 0.0;
308 atomData->addAtomInfo(atomInfo);
309
310 newVec = rotTrans * c2;
311 atomInfo = new AtomInfo;
312 atomInfo->atomTypeName = "C";
313 atomInfo->pos[0] = pos[0] + newVec[0];
314 atomInfo->pos[1] = pos[1] + newVec[1];
315 atomInfo->pos[2] = pos[2] + newVec[2];
316 atomInfo->dipole[0] = 0.0;
317 atomInfo->dipole[1] = 0.0;
318 atomInfo->dipole[2] = 0.0;
319 atomData->addAtomInfo(atomInfo);
320
321 newVec = rotTrans * c3;
322 atomInfo = new AtomInfo;
323 atomInfo->atomTypeName = "C";
324 atomInfo->pos[0] = pos[0] + newVec[0];
325 atomInfo->pos[1] = pos[1] + newVec[1];
326 atomInfo->pos[2] = pos[2] + newVec[2];
327 atomInfo->dipole[0] = 0.0;
328 atomInfo->dipole[1] = 0.0;
329 atomInfo->dipole[2] = 0.0;
330 atomData->addAtomInfo(atomInfo);
331
332 newVec = rotTrans * c4;
333 atomInfo = new AtomInfo;
334 atomInfo->atomTypeName = "C";
335 atomInfo->pos[0] = pos[0] + newVec[0];
336 atomInfo->pos[1] = pos[1] + newVec[1];
337 atomInfo->pos[2] = pos[2] + newVec[2];
338 atomInfo->dipole[0] = 0.0;
339 atomInfo->dipole[1] = 0.0;
340 atomInfo->dipole[2] = 0.0;
341 atomData->addAtomInfo(atomInfo);
342
343 //add atom data into atom's property
344
345 if(!haveAtomData){
346 atomData->setID("ATOMDATA");
347 datom->addProperty(atomData);
348 }
349
350 setVisited(datom);
351
352 }
353
354 const std::string LinearAtomVisitor::toString(){
355 char buffer[65535];
356 std::string result;
357
358 sprintf(buffer ,"------------------------------------------------------------------\n");
359 result += buffer;
360
361 sprintf(buffer ,"Visitor name: %s\n", visitorName.c_str());
362 result += buffer;
363
364 sprintf(buffer , "Visitor Description: Convert linear into 4 different atoms\n");
365 result += buffer;
366
367 sprintf(buffer ,"------------------------------------------------------------------\n");
368 result += buffer;
369
370 return result;
371 }
372
373 bool GBLipidAtomVisitor::isGBLipidAtom(const std::string& atomType){
374 std::set<std::string>::iterator strIter;
375 strIter = GBLipidAtomType.find(atomType);
376
377 return strIter != GBLipidAtomType.end() ? true : false;
378 }
379
380 void GBLipidAtomVisitor::visit(DirectionalAtom* datom){
381 std::vector<AtomInfo*> atoms;
382 //we need to convert linear into 4 different atoms
383 Vector3d c1(0.0, 0.0, -6.25);
384 Vector3d c2(0.0, 0.0, -2.1);
385 Vector3d c3(0.0, 0.0, 2.1);
386 Vector3d c4(0.0, 0.0, 6.25);
387 RotMat3x3d rotMatrix;
388 RotMat3x3d rotTrans;
389 AtomInfo* atomInfo;
390 Vector3d pos;
391 Vector3d newVec;
392 Quat4d q;
393 AtomData* atomData;
394 GenericData* data;
395 bool haveAtomData;
396
397 //if atom is not GBlipid atom, just skip it
398 if(!isGBLipidAtom(datom->getType()))
399 return;
400
401 data = datom->getPropertyByName("ATOMDATA");
402 if(data != NULL){
403 atomData = dynamic_cast<AtomData*>(data);
404 if(atomData == NULL){
405 std::cerr << "can not get Atom Data from " << datom->getType() << std::endl;
406 atomData = new AtomData;
407 haveAtomData = false;
408 } else {
409 haveAtomData = true;
410 }
411 } else {
412 atomData = new AtomData;
413 haveAtomData = false;
414 }
415
416
417 pos = datom->getPos();
418 q = datom->getQ();
419 rotMatrix = datom->getA();
420
421 // We need A^T to convert from body-fixed to space-fixed:
422 rotTrans = rotMatrix.transpose();
423
424 newVec = rotTrans * c1;
425 atomInfo = new AtomInfo;
426 atomInfo->atomTypeName = "K";
427 atomInfo->pos[0] = pos[0] + newVec[0];
428 atomInfo->pos[1] = pos[1] + newVec[1];
429 atomInfo->pos[2] = pos[2] + newVec[2];
430 atomInfo->dipole[0] = 0.0;
431 atomInfo->dipole[1] = 0.0;
432 atomInfo->dipole[2] = 0.0;
433 atomData->addAtomInfo(atomInfo);
434
435 newVec = rotTrans * c2;
436 atomInfo = new AtomInfo;
437 atomInfo->atomTypeName = "K";
438 atomInfo->pos[0] = pos[0] + newVec[0];
439 atomInfo->pos[1] = pos[1] + newVec[1];
440 atomInfo->pos[2] = pos[2] + newVec[2];
441 atomInfo->dipole[0] = 0.0;
442 atomInfo->dipole[1] = 0.0;
443 atomInfo->dipole[2] = 0.0;
444 atomData->addAtomInfo(atomInfo);
445
446 newVec = rotTrans * c3;
447 atomInfo = new AtomInfo;
448 atomInfo->atomTypeName = "K";
449 atomInfo->pos[0] = pos[0] + newVec[0];
450 atomInfo->pos[1] = pos[1] + newVec[1];
451 atomInfo->pos[2] = pos[2] + newVec[2];
452 atomInfo->dipole[0] = 0.0;
453 atomInfo->dipole[1] = 0.0;
454 atomInfo->dipole[2] = 0.0;
455 atomData->addAtomInfo(atomInfo);
456
457 newVec = rotTrans * c4;
458 atomInfo = new AtomInfo;
459 atomInfo->atomTypeName = "K";
460 atomInfo->pos[0] = pos[0] + newVec[0];
461 atomInfo->pos[1] = pos[1] + newVec[1];
462 atomInfo->pos[2] = pos[2] + newVec[2];
463 atomInfo->dipole[0] = 0.0;
464 atomInfo->dipole[1] = 0.0;
465 atomInfo->dipole[2] = 0.0;
466 atomData->addAtomInfo(atomInfo);
467
468 //add atom data into atom's property
469
470 if(!haveAtomData){
471 atomData->setID("ATOMDATA");
472 datom->addProperty(atomData);
473 }
474
475 setVisited(datom);
476
477 }
478
479 const std::string GBLipidAtomVisitor::toString(){
480 char buffer[65535];
481 std::string result;
482
483 sprintf(buffer ,"------------------------------------------------------------------\n");
484 result += buffer;
485
486 sprintf(buffer ,"Visitor name: %s\n", visitorName.c_str());
487 result += buffer;
488
489 sprintf(buffer , "Visitor Description: Convert GBlipid into 4 different K atoms\n");
490 result += buffer;
491
492 sprintf(buffer ,"------------------------------------------------------------------\n");
493 result += buffer;
494
495 return result;
496 }
497
498 //----------------------------------------------------------------------------//
499
500 void DefaultAtomVisitor::visit(Atom *atom) {
501 AtomData *atomData;
502 AtomInfo *atomInfo;
503 Vector3d pos;
504
505 if (isVisited(atom))
506 return;
507
508 atomInfo = new AtomInfo;
509
510 atomData = new AtomData;
511 atomData->setID("ATOMDATA");
512
513 pos = atom->getPos();
514 atomInfo->atomTypeName = atom->getType();
515 atomInfo->pos[0] = pos[0];
516 atomInfo->pos[1] = pos[1];
517 atomInfo->pos[2] = pos[2];
518 atomInfo->dipole[0] = 0.0;
519 atomInfo->dipole[1] = 0.0;
520 atomInfo->dipole[2] = 0.0;
521
522 atomData->addAtomInfo(atomInfo);
523
524 atom->addProperty(atomData);
525
526 setVisited(atom);
527 }
528
529 void DefaultAtomVisitor::visit(DirectionalAtom *datom) {
530 AtomData *atomData;
531 AtomInfo *atomInfo;
532 Vector3d pos;
533 Vector3d u;
534
535 if (isVisited(datom))
536 return;
537
538 pos = datom->getPos();
539 if (datom->getAtomType()->isGayBerne()) {
540 u = datom->getA().transpose()*V3Z;
541 } else if (datom->getAtomType()->isMultipole()) {
542 u = datom->getElectroFrame().getColumn(2);
543 }
544 atomData = new AtomData;
545 atomData->setID("ATOMDATA");
546 atomInfo = new AtomInfo;
547
548 atomInfo->atomTypeName = datom->getType();
549 atomInfo->pos[0] = pos[0];
550 atomInfo->pos[1] = pos[1];
551 atomInfo->pos[2] = pos[2];
552 atomInfo->dipole[0] = u[0];
553 atomInfo->dipole[1] = u[1];
554 atomInfo->dipole[2] = u[2];
555
556 atomData->addAtomInfo(atomInfo);
557
558 datom->addProperty(atomData);
559
560 setVisited(datom);
561 }
562
563 const std::string DefaultAtomVisitor::toString() {
564 char buffer[65535];
565 std::string result;
566
567 sprintf(buffer,
568 "------------------------------------------------------------------\n");
569 result += buffer;
570
571 sprintf(buffer, "Visitor name: %s\n", visitorName.c_str());
572 result += buffer;
573
574 sprintf(buffer,
575 "Visitor Description: copy atom infomation into atom data\n");
576 result += buffer;
577
578 sprintf(buffer,
579 "------------------------------------------------------------------\n");
580 result += buffer;
581
582 return result;
583 }
584 } //namespace oopse

Properties

Name Value
svn:executable *