ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/visitors/AtomVisitor.cpp
(Generate patch)

Comparing trunk/src/visitors/AtomVisitor.cpp (file contents):
Revision 954 by tim, Wed May 10 01:44:48 2006 UTC vs.
Revision 1219 by xsun, Wed Jan 23 21:23:32 2008 UTC

# Line 82 | Line 82 | namespace oopse {
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
85 >    //we need to convert SSD into 4 different atoms
86 >    //one oxygen atom, two hydrogen atoms and one pseudo atom which is the center of
87 >    //the mass 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);
# Line 205 | Line 205 | namespace oopse {
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 +
218 +  bool TREDAtomVisitor::isTREDAtom(const std::string&atomType) {
219 +    std::set<std::string>::iterator strIter;
220 +    strIter = tredAtomType.find(atomType);
221 +    return strIter != tredAtomType.end() ? true : false;
222 +  }
223 +
224 +  void TREDAtomVisitor::visit(DirectionalAtom *datom) {
225 +    std::vector<AtomInfo*>atoms;
226 +
227 +    // we need to convert a TRED into 4 different atoms:
228 +    // one oxygen atom, two hydrogen atoms, and one atom which is the center of
229 +    // the mass of the water with a dipole moment
230 +    Vector3d h1(0.0, -0.75695, 0.5206);
231 +    Vector3d h2(0.0, 0.75695, 0.5206);
232 +    Vector3d ox(0.0, 0.0, -0.0654);
233 +    Vector3d u(0, 0, 1);
234 +    RotMat3x3d   rotMatrix;
235 +    RotMat3x3d   rotTrans;
236 +    AtomInfo *   atomInfo;
237 +    Vector3d     pos;
238 +    Vector3d     newVec;
239 +    Quat4d       q;
240 +    AtomData *   atomData;
241 +    GenericData *data;
242 +    bool         haveAtomData;
243 +
244 +    // if the atom is not a TRED atom, skip it
245 +    if (!isTREDAtom(datom->getType()))
246 +      return;
247 +
248 +    data = datom->getPropertyByName("ATOMDATA");
249 +
250 +    if (data != NULL) {
251 +      atomData = dynamic_cast<AtomData *>(data);
252 +
253 +      if (atomData == NULL) {
254 +        std::cerr << "can not get Atom Data from " << datom->getType() << std::endl;
255 +        atomData = new AtomData;
256 +        haveAtomData = false;
257 +      } else
258 +        haveAtomData = true;
259 +    } else {
260 +      atomData = new AtomData;
261 +      haveAtomData = false;
262 +    }
263 +
264 +    pos = datom->getPos();
265 +    q = datom->getQ();
266 +    rotMatrix = datom->getA();
267 +
268 +    // We need A^T to convert from body-fixed to space-fixed:
269 +    // transposeMat3(rotMatrix, rotTrans);
270 +    rotTrans = rotMatrix.transpose();
271 +
272 +    // center of mass of the water molecule
273 +    // matVecMul3(rotTrans, u, newVec);
274 +    newVec = rotTrans * u;
275 +
276 +    atomInfo = new AtomInfo;
277 +    atomInfo->atomTypeName = "TRED";
278 +    atomInfo->pos[0] = pos[0];
279 +    atomInfo->pos[1] = pos[1];
280 +    atomInfo->pos[2] = pos[2];
281 +    atomInfo->dipole[0] = newVec[0];
282 +    atomInfo->dipole[1] = newVec[1];
283 +    atomInfo->dipole[2] = newVec[2];
284 +
285 +    atomData->addAtomInfo(atomInfo);
286 +
287 +    // oxygen
288 +    // matVecMul3(rotTrans, ox, newVec);
289 +    newVec = rotTrans * ox;
290 +
291 +    atomInfo = new AtomInfo;
292 +    atomInfo->atomTypeName = "O";
293 +    atomInfo->pos[0] = pos[0] + newVec[0];
294 +    atomInfo->pos[1] = pos[1] + newVec[1];
295 +    atomInfo->pos[2] = pos[2] + newVec[2];
296 +    atomInfo->dipole[0] = 0.0;
297 +    atomInfo->dipole[1] = 0.0;
298 +    atomInfo->dipole[2] = 0.0;
299 +    atomData->addAtomInfo(atomInfo);
300 +
301 +    // hydrogen1
302 +    // matVecMul3(rotTrans, h1, newVec);
303 +    newVec = rotTrans * h1;
304 +    atomInfo = new AtomInfo;
305 +    atomInfo->atomTypeName = "H";
306 +    atomInfo->pos[0] = pos[0] + newVec[0];
307 +    atomInfo->pos[1] = pos[1] + newVec[1];
308 +    atomInfo->pos[2] = pos[2] + newVec[2];
309 +    atomInfo->dipole[0] = 0.0;
310 +    atomInfo->dipole[1] = 0.0;
311 +    atomInfo->dipole[2] = 0.0;
312 +    atomData->addAtomInfo(atomInfo);
313 +
314 +    // hydrogen2
315 +    // matVecMul3(rotTrans, h2, newVec);
316 +    newVec = rotTrans * h2;
317 +    atomInfo = new AtomInfo;
318 +    atomInfo->atomTypeName = "H";
319 +    atomInfo->pos[0] = pos[0] + newVec[0];
320 +    atomInfo->pos[1] = pos[1] + newVec[1];
321 +    atomInfo->pos[2] = pos[2] + newVec[2];
322 +    atomInfo->dipole[0] = 0.0;
323 +    atomInfo->dipole[1] = 0.0;
324 +    atomInfo->dipole[2] = 0.0;
325 +    atomData->addAtomInfo(atomInfo);
326 +
327 +    // add atom data into atom's property
328 +
329 +    if (!haveAtomData) {
330 +      atomData->setID("ATOMDATA");
331 +      datom->addProperty(atomData);
332 +    }
333 +
334 +    setVisited(datom);
335 +  }
336 +
337 +  const std::string TREDAtomVisitor::toString() {
338 +    char   buffer[65535];
339 +    std::string result;
340 +
341 +    sprintf(buffer,
342 +            "------------------------------------------------------------------\n");
343      result += buffer;
344  
345 +    sprintf(buffer, "Visitor name: %s\n", visitorName.c_str());
346 +    result += buffer;
347 +
348      sprintf(buffer,
349 +            "Visitor Description: Convert the TRED atom into 4 different atoms\n");
350 +    result += buffer;
351 +
352 +    sprintf(buffer,
353              "------------------------------------------------------------------\n");
354      result += buffer;
355  
356      return result;
357    }
358  
359 +
360    bool LinearAtomVisitor::isLinearAtom(const std::string& atomType){
361      std::set<std::string>::iterator strIter;
362      strIter = linearAtomType.find(atomType);
# Line 241 | Line 384 | namespace oopse {
384      AtomData* atomData;
385      GenericData* data;
386      bool haveAtomData;
387 <
388 <    //if atom is not SSD atom, just skip it
389 <    if(!isLinearAtom(datom->getType()))
387 >    AtomType* atomType;
388 >    //if atom is not linear atom, just skip it
389 >    if(!isLinearAtom(datom->getType()) || !datom->getAtomType()->isGayBerne())
390        return;
391  
392 +    //setup GayBerne type in fortran side
393 +    data = datom->getAtomType()->getPropertyByName("GayBerne");
394 +    if (data != NULL) {
395 +       GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
396 +
397 +       if (gayBerneData != NULL) {
398 +           GayBerneParam gayBerneParam = gayBerneData->getData();
399 +
400 +                          // double halfLen = gayBerneParam.GB_sigma * gayBerneParam.GB_l2b_ratio/2.0;
401 +                          double halfLen = gayBerneParam.GB_l/2.0;
402 +                          c1[2] = -halfLen;
403 +              c2[2] = -halfLen /2;
404 +              c3[2] = halfLen/2;
405 +              c4[2] = halfLen;
406 +                
407 +            }
408 +            
409 +              else {
410 +                    sprintf( painCave.errMsg,
411 +                           "Can not cast GenericData to GayBerneParam\n");
412 +                    painCave.severity = OOPSE_ERROR;
413 +                    painCave.isFatal = 1;
414 +                    simError();          
415 +        }            
416 +    }
417 +
418 +
419      data = datom->getPropertyByName("ATOMDATA");
420      if(data != NULL){
421        atomData = dynamic_cast<AtomData*>(data);  
# Line 353 | Line 523 | namespace oopse {
523    void GBLipidAtomVisitor::visit(DirectionalAtom* datom){
524      std::vector<AtomInfo*> atoms;
525      //we need to convert linear into 4 different atoms
526 <    Vector3d c1(0.0, 0.0, -6.25);
527 <    Vector3d c2(0.0, 0.0, -2.1);
528 <    Vector3d c3(0.0, 0.0,  2.1);
529 <    Vector3d c4(0.0, 0.0,  6.25);
526 >    //Vector3d c1(0.0, 0.0, -6.25);
527 >    //Vector3d c2(0.0, 0.0, -2.1);
528 >    //Vector3d c3(0.0, 0.0,  2.1);
529 >    //Vector3d c4(0.0, 0.0,  6.25);
530 >    Vector3d c1(0.0, 0.0, 0.0);
531 >    Vector3d c2(0.0, 0.0, 1.0);
532      RotMat3x3d rotMatrix;
533      RotMat3x3d rotTrans;
534      AtomInfo* atomInfo;
535      Vector3d pos;
536      Vector3d newVec;
537 +    Vector3d dVec;
538      Quat4d q;
539      AtomData* atomData;
540      GenericData* data;
# Line 395 | Line 568 | namespace oopse {
568      rotTrans = rotMatrix.transpose();
569  
570      newVec = rotTrans * c1;
571 +    dVec = rotTrans * c2;
572      atomInfo = new AtomInfo;
573 +    atomInfo->atomTypeName = "GB";
574 +    atomInfo->pos[0] = pos[0] + newVec[0];
575 +    atomInfo->pos[1] = pos[1] + newVec[1];
576 +    atomInfo->pos[2] = pos[2] + newVec[2];
577 +    atomInfo->dipole[0] = dVec[0];
578 +    atomInfo->dipole[1] = dVec[1];
579 +    atomInfo->dipole[2] = dVec[2];
580 +    atomInfo->hasVector = true;
581 +    atomInfo->charge = 3.0;
582 +    atomInfo->hasCharge = true;
583 +    atomData->addAtomInfo(atomInfo);
584 +
585 +    //newVec = rotTrans * c2;
586 +    //atomInfo = new AtomInfo;
587 +    //atomInfo->atomTypeName = "K";
588 +    //atomInfo->pos[0] = pos[0] + newVec[0];
589 +    //atomInfo->pos[1] = pos[1] + newVec[1];
590 +    //atomInfo->pos[2] = pos[2] + newVec[2];
591 +    //atomInfo->dipole[0] = 0.0;
592 +    //atomInfo->dipole[1] = 0.0;
593 +    //atomInfo->dipole[2] = 0.0;
594 +    //atomData->addAtomInfo(atomInfo);
595 +
596 +    //newVec = rotTrans * c3;
597 +    //atomInfo = new AtomInfo;
598 +    //atomInfo->atomTypeName = "K";
599 +    //atomInfo->pos[0] = pos[0] + newVec[0];
600 +    //atomInfo->pos[1] = pos[1] + newVec[1];
601 +    //atomInfo->pos[2] = pos[2] + newVec[2];
602 +    //atomInfo->dipole[0] = 0.0;
603 +    //atomInfo->dipole[1] = 0.0;
604 +    //atomInfo->dipole[2] = 0.0;
605 +    //atomData->addAtomInfo(atomInfo);
606 +
607 +    //newVec = rotTrans * c4;
608 +    //atomInfo = new AtomInfo;
609 +    //atomInfo->atomTypeName = "K";
610 +    //atomInfo->pos[0] = pos[0] + newVec[0];
611 +    //atomInfo->pos[1] = pos[1] + newVec[1];
612 +    //atomInfo->pos[2] = pos[2] + newVec[2];
613 +    //atomInfo->dipole[0] = 0.0;
614 +    //atomInfo->dipole[1] = 0.0;
615 +    //atomInfo->dipole[2] = 0.0;
616 +    //atomData->addAtomInfo(atomInfo);
617 +
618 +    //add atom data into atom's property
619 +
620 +    if(!haveAtomData){
621 +      atomData->setID("ATOMDATA");
622 +      datom->addProperty(atomData);
623 +    }
624 +
625 +    setVisited(datom);
626 +
627 +  }
628 +
629 +  const std::string GBLipidAtomVisitor::toString(){
630 +    char buffer[65535];
631 +    std::string result;
632 +  
633 +    sprintf(buffer ,"------------------------------------------------------------------\n");
634 +    result += buffer;
635 +
636 +    sprintf(buffer ,"Visitor name: %s\n", visitorName.c_str());
637 +    result += buffer;
638 +
639 +    sprintf(buffer , "Visitor Description: Convert GBlipid into 4 different K atoms\n");
640 +    result += buffer;
641 +
642 +    sprintf(buffer ,"------------------------------------------------------------------\n");
643 +    result += buffer;
644 +
645 +    return result;
646 +  }
647 +
648 +  bool Ring5gbAtomVisitor::isRing5gbAtom(const std::string& atomType){
649 +    std::set<std::string>::iterator strIter;
650 +    strIter = Ring5gbAtomType.find(atomType);
651 +
652 +    return strIter != Ring5gbAtomType.end() ? true : false;
653 +  }
654 +
655 +  void Ring5gbAtomVisitor::visit(DirectionalAtom* datom){
656 +    std::vector<AtomInfo*> atoms;
657 +    //we need to convert linear into 4 different atoms
658 +    Vector3d c1(0.0, 0.0, -5.5);
659 +    Vector3d c2(0.0, 0.0, -1.8);
660 +    Vector3d c3(0.0, 0.0,  1.8);
661 +    Vector3d c4(0.0, 0.0,  5.5);
662 +    RotMat3x3d rotMatrix;
663 +    RotMat3x3d rotTrans;
664 +    AtomInfo* atomInfo;
665 +    Vector3d pos;
666 +    Vector3d newVec;
667 +    Vector3d dVec;
668 +    Quat4d q;
669 +    AtomData* atomData;
670 +    GenericData* data;
671 +    bool haveAtomData;
672 +
673 +    //if atom is not Ring5GB atom, just skip it
674 +    if(!isRing5gbAtom(datom->getType()))
675 +      return;
676 +
677 +    data = datom->getPropertyByName("ATOMDATA");
678 +    if(data != NULL){
679 +      atomData = dynamic_cast<AtomData*>(data);  
680 +      if(atomData == NULL){
681 +        std::cerr << "can not get Atom Data from " << datom->getType() << std::endl;
682 +        atomData = new AtomData;
683 +        haveAtomData = false;      
684 +      } else {
685 +        haveAtomData = true;
686 +      }
687 +    } else {
688 +      atomData = new AtomData;
689 +      haveAtomData = false;
690 +    }
691 +  
692 +  
693 +    pos = datom->getPos();
694 +    q = datom->getQ();
695 +    rotMatrix = datom->getA();
696 +
697 +    // We need A^T to convert from body-fixed to space-fixed:  
698 +    rotTrans = rotMatrix.transpose();
699 +
700 +    newVec = rotTrans * c1;
701 +    atomInfo = new AtomInfo;
702      atomInfo->atomTypeName = "K";
703      atomInfo->pos[0] = pos[0] + newVec[0];
704      atomInfo->pos[1] = pos[1] + newVec[1];
# Line 449 | Line 752 | namespace oopse {
752  
753    }
754  
755 <  const std::string GBLipidAtomVisitor::toString(){
755 >  const std::string Ring5gbAtomVisitor::toString(){
756      char buffer[65535];
757      std::string result;
758    
# Line 459 | Line 762 | namespace oopse {
762      sprintf(buffer ,"Visitor name: %s\n", visitorName.c_str());
763      result += buffer;
764  
765 <    sprintf(buffer , "Visitor Description: Convert GBlipid into 4 different K atoms\n");
765 >    sprintf(buffer , "Visitor Description: Convert Ring5GB into 4 different K atoms\n");
766 >    result += buffer;
767 >
768 >    sprintf(buffer ,"------------------------------------------------------------------\n");
769      result += buffer;
770  
771 +    return result;
772 +  }
773 +
774 +  bool HeadAtomVisitor::isHeadAtom(const std::string& atomType){
775 +    std::set<std::string>::iterator strIter;
776 +    strIter = HeadAtomType.find(atomType);
777 +
778 +    return strIter != HeadAtomType.end() ? true : false;
779 +  }
780 +
781 +  void HeadAtomVisitor::visit(DirectionalAtom* datom){
782 +    std::vector<AtomInfo*> atoms;
783 +    //we need to convert linear into 2 different atoms
784 +    Vector3d c1(0.0, 0.0, -1.5);
785 +    Vector3d c2(0.0, 0.0, 1.5);
786 +    RotMat3x3d rotMatrix;
787 +    RotMat3x3d rotTrans;
788 +    AtomInfo* atomInfo;
789 +    Vector3d pos;
790 +    Vector3d newVec;
791 +    Vector3d dVec;
792 +    Quat4d q;
793 +    AtomData* atomData;
794 +    GenericData* data;
795 +    bool haveAtomData;
796 +
797 +    //if atom is not Head atom, just skip it
798 +    if(!isHeadAtom(datom->getType()))
799 +      return;
800 +
801 +    data = datom->getPropertyByName("ATOMDATA");
802 +    if(data != NULL){
803 +      atomData = dynamic_cast<AtomData*>(data);  
804 +      if(atomData == NULL){
805 +        std::cerr << "can not get Atom Data from " << datom->getType() << std::endl;
806 +        atomData = new AtomData;
807 +        haveAtomData = false;      
808 +      } else {
809 +        haveAtomData = true;
810 +      }
811 +    } else {
812 +      atomData = new AtomData;
813 +      haveAtomData = false;
814 +    }
815 +  
816 +  
817 +    pos = datom->getPos();
818 +    q = datom->getQ();
819 +    rotMatrix = datom->getA();
820 +
821 +    // We need A^T to convert from body-fixed to space-fixed:  
822 +    rotTrans = rotMatrix.transpose();
823 +
824 +    newVec = rotTrans * c1;
825 +    atomInfo = new AtomInfo;
826 +    atomInfo->atomTypeName = "C";
827 +    atomInfo->pos[0] = pos[0] + newVec[0];
828 +    atomInfo->pos[1] = pos[1] + newVec[1];
829 +    atomInfo->pos[2] = pos[2] + newVec[2];
830 +    atomInfo->dipole[0] = 0.0;
831 +    atomInfo->dipole[1] = 0.0;
832 +    atomInfo->dipole[2] = 0.0;
833 +    atomData->addAtomInfo(atomInfo);
834 +
835 +    newVec = rotTrans * c2;
836 +    atomInfo = new AtomInfo;
837 +    atomInfo->atomTypeName = "O";
838 +    atomInfo->pos[0] = pos[0] + newVec[0];
839 +    atomInfo->pos[1] = pos[1] + newVec[1];
840 +    atomInfo->pos[2] = pos[2] + newVec[2];
841 +    atomInfo->dipole[0] = 0.0;
842 +    atomInfo->dipole[1] = 0.0;
843 +    atomInfo->dipole[2] = 0.0;
844 +    atomData->addAtomInfo(atomInfo);
845 +
846 +    //add atom data into atom's property
847 +
848 +    if(!haveAtomData){
849 +      atomData->setID("ATOMDATA");
850 +      datom->addProperty(atomData);
851 +    }
852 +
853 +    setVisited(datom);
854 +
855 +  }
856 +
857 +  const std::string HeadAtomVisitor::toString(){
858 +    char buffer[65535];
859 +    std::string result;
860 +  
861      sprintf(buffer ,"------------------------------------------------------------------\n");
862      result += buffer;
863  
864 +    sprintf(buffer ,"Visitor name: %s\n", visitorName.c_str());
865 +    result += buffer;
866 +
867 +    sprintf(buffer , "Visitor Description: Convert HEAD into C atom and O atom\n");
868 +    result += buffer;
869 +
870 +    sprintf(buffer ,"------------------------------------------------------------------\n");
871 +    result += buffer;
872 +
873      return result;
874    }
875  
876 +
877    //----------------------------------------------------------------------------//
878  
879    void DefaultAtomVisitor::visit(Atom *atom) {
# Line 509 | Line 915 | namespace oopse {
915        return;
916  
917      pos = datom->getPos();
918 <    u = datom->getElectroFrame().getColumn(2);
919 <
918 >    if (datom->getAtomType()->isGayBerne()) {
919 >        u = datom->getA().transpose()*V3Z;        
920 >    } else if (datom->getAtomType()->isMultipole()) {
921 >        u = datom->getElectroFrame().getColumn(2);
922 >    }
923      atomData = new AtomData;
924      atomData->setID("ATOMDATA");
925      atomInfo = new AtomInfo;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines