src/types/MoleculeStamp.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.
 */
 
#include <iostream>
#include <functional>
#include "types/MoleculeStamp.hpp"
#include "utils/Tuple.hpp"

namespace oopse {
MoleculeStamp::MoleculeStamp() {
    DefineParameter(Name, "name");
    
    deprecatedKeywords_.insert("nAtoms");
    deprecatedKeywords_.insert("nBonds");
    deprecatedKeywords_.insert("nBends");
    deprecatedKeywords_.insert("nTorsions");
    deprecatedKeywords_.insert("nRigidBodies");
    deprecatedKeywords_.insert("nCutoffGroups");
    
}

MoleculeStamp::~MoleculeStamp() {

}

bool MoleculeStamp::addAtomStamp( AtomStamp* atom) {
    bool ret = addIndexSensitiveStamp(atomStamps_, atom);
    if (!ret) {
        std::cout << "multiple atoms have the same index: " << atom->getIndex() <<" in " << getName()  << " Molecule\n";
    }
    return ret;
    
}

bool MoleculeStamp::addBondStamp( BondStamp* bond) {
    bondStamps_.push_back(bond);
    return true;
}

bool MoleculeStamp::addBendStamp( BendStamp* bend) {
    bendStamps_.push_back(bend);
    return true;
}

bool MoleculeStamp::addTorsionStamp( TorsionStamp* torsion) {
    torsionStamps_.push_back(torsion);
    return true;
}

bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp* rigidbody) {
    bool ret = addIndexSensitiveStamp(rigidBodyStamps_, rigidbody);
    if (!ret) {
        std::cout << "multiple rigidbodies have the same index: " << rigidbody->getIndex() <<" in " << getName()  << " Molecule\n";
    }
    return ret;
}

bool MoleculeStamp::addCutoffGroupStamp( CutoffGroupStamp* cutoffgroup) {
    cutoffGroupStamps_.push_back(cutoffgroup);
    return true;
}

bool MoleculeStamp::addFragmentStamp( FragmentStamp* fragment) {
    return addIndexSensitiveStamp(fragmentStamps_, fragment);
}
    
void MoleculeStamp::validate() {
    DataHolder::validate();

    atom2Rigidbody.resize(getNAtoms()); 
    // negative number means atom is a free atom, does not belong to rigidbody
    //every element in atom2Rigidbody has unique negative number at the very beginning
    for(int i = 0; i < atom2Rigidbody.size(); ++i) {
        atom2Rigidbody[i] = -1 - i;
    }
    for (int i = 0; i < getNRigidBodies(); ++i) {
        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
        std::vector<int> members = rbStamp->getMembers();
        for(std::vector<int>::iterator j = members.begin(); j != members.end(); ++j) {
            atom2Rigidbody[*j] = i;                 
        }
    }

    checkAtoms();
    checkBonds();
    fillBondInfo();
    checkBends();
    checkTorsions();
    checkRigidBodies();
    checkCutoffGroups();
    checkFragments();

    int nrigidAtoms = 0;
    for (int i = 0; i < getNRigidBodies(); ++i) {
        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
        nrigidAtoms += rbStamp->getNMembers();
    }
    nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;

 }

void MoleculeStamp::checkAtoms() {
    std::vector<AtomStamp*>::iterator ai = std::find(atomStamps_.begin(), atomStamps_.end(), static_cast<AtomStamp*>(NULL));
    if (ai != atomStamps_.end()) {
        std::cout << "Error in Molecule " << getName() << ": atom[" << ai - atomStamps_.begin()<< "] is missing\n";
    }

    //make sure index is not out of range
    int natoms = getNAtoms();
    for(int i = 0; i < getNBonds(); ++i) {
        BondStamp* bondStamp = getBondStamp(i);
        if (bondStamp->getA() >=  natoms && bondStamp->getB() >= natoms) {
            std::cout << "Error in Molecule " << getName() <<  ": bond between " << bondStamp->getA() << " and " << bondStamp->getB() << " is invalid\n";
        }
    }
}

void MoleculeStamp::checkBonds() {
    //make sure bonds are unique
    std::set<std::pair<int, int> > allBonds;
    for(int i = 0; i < getNBonds(); ++i) {
        BondStamp* bondStamp= getBondStamp(i);        
        std::pair<int, int> bondPair(bondStamp->getA(), bondStamp->getB());
        //make sure bondTuple.first is always less than or equal to bondTuple.third
        if (bondPair.first > bondPair.second) {
            std::swap(bondPair.first, bondPair.second);
        }
        
        std::set<std::pair<int, int> >::iterator iter = allBonds.find(bondPair);
        if ( iter != allBonds.end()) {
            std::cout << "Error in Molecule " << getName() << ": " << "Bond appears multiple times\n";
        } else {
            allBonds.insert(bondPair);
        }
    }
    
    //make sure atoms belong to same rigidbody do not bond to each other
    for(int i = 0; i < getNBonds(); ++i) {
        BondStamp* bondStamp = getBondStamp(i);
        if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()]) {
            std::cout << "Error in Molecule " << getName() << ": "<<"bond between " << bondStamp->getA() << " and " << bondStamp->getB() << " belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
        }
    }
    
}

struct BendLessThan : public std::binary_function<IntTuple4, IntTuple4, bool> {
    bool operator()(IntTuple3 b1, IntTuple3 b2) {
        return b1.first < b2.first
             || (!(b2.first < b1.first) && b1.second < b2.second)
             || (!(b2.first < b1.first) && !(b2.second < b2.second) && b1.third < b2.third);
    }
};

void MoleculeStamp::checkBends() {
    for(int i = 0; i < getNBends(); ++i) {
        BendStamp* bendStamp = getBendStamp(i);
        std::vector<int> bendAtoms =  bendStamp->getMembers();
        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
        if (j != bendAtoms.end()) {
            std::cout << "Error in Molecule " << getName();
        }

        if (bendAtoms.size() == 2 ) {
            if (!bendStamp->haveGhostVectorSource()) {
                std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing\n";
            }else{
                int ghostIndex = bendStamp->getGhostVectorSource();
                if (ghostIndex < getNAtoms()) {
                    if (std::find(bendAtoms.begin(), bendAtoms.end(), ghostIndex) == bendAtoms.end()) {
                      std::cout <<  "Error in Molecule " << getName() << ": ghostVectorSouce "<< ghostIndex<<"is invalid\n"; 
                    }
                    if (!getAtomStamp(ghostIndex)->haveOrientation()) {
                        std::cout <<  "Error in Molecule " << getName() << ": ghost atom must be a directioanl atom\n"; 
                    }
                }else {
                    std::cout << "Error in Molecule " << getName() <<  ": ghostVectorsource " << ghostIndex<< "  is invalid\n";
                }
            }
        } else if (bendAtoms.size() == 3 && bendStamp->haveGhostVectorSource()) {
            std::cout <<  "Error in Molecule " << getName() << ": normal bend should not have ghostVectorSouce\n"; 
        }
    }

    for(int i = 0; i < getNBends(); ++i) {
        BendStamp* bendStamp = getBendStamp(i);
        std::vector<int> bendAtoms =  bendStamp->getMembers();
        std::vector<int> rigidSet(getNRigidBodies(), 0);
        std::vector<int>::iterator j;
        for( j = bendAtoms.begin(); j != bendAtoms.end(); ++j) {
            int rigidbodyIndex = atom2Rigidbody[*j];
            if (rigidbodyIndex >= 0) {
                ++rigidSet[rigidbodyIndex];
                if (rigidSet[rigidbodyIndex] > 1) {
                    std::cout << "Error in Molecule " << getName() << ": ";
                    //std::cout << "atoms of bend " <<  << "belong to same rigidbody " << rigidbodyIndex << "\n";                    
                }
            }
        }
    } 
    
    
    std::set<IntTuple3, BendLessThan> allBends;
    std::set<IntTuple3, BendLessThan>::iterator iter;
    for(int i = 0; i < getNBends(); ++i) {
        BendStamp* bendStamp= getBendStamp(i);
        std::vector<int> bend = bendStamp->getMembers();
        if (bend.size() == 2) {
        // in case we have two ghost bend. For example, 
        // bend {
        // members (0, 1);
        //   ghostVectorSource = 0;
        // }
        // and
        // bend {
        //   members (0, 1);
        // ghostVectorSource = 0;
        // }
        // In order to distinguish them. we expand them to Tuple3.
        // the first one is expanded to (0, 0, 1) while the second one is expaned to (0, 1, 1)
             int ghostIndex = bendStamp->getGhostVectorSource();
             std::vector<int>::iterator j = std::find(bend.begin(), bend.end(), ghostIndex);
             if (j != bend.end()) {
                bend.insert(j, ghostIndex);
             }
        }
        
        IntTuple3 bendTuple(bend[0], bend[1], bend[2]);
        //make sure bendTuple.first is always less than or equal to bendTuple.third
        if (bendTuple.first > bendTuple.third) {
            std::swap(bendTuple.first, bendTuple.third);
        }
        
        iter = allBends.find(bendTuple);
        if ( iter != allBends.end()) {
            std::cout << "Error in Molecule " << getName() << ": " << "Bend appears multiple times\n";
        } else {
            allBends.insert(bendTuple);
        }
    }

    for (int i = 0; i < getNBonds(); ++i) {
        BondStamp* bondStamp = getBondStamp(i);
        int a = bondStamp->getA();
        int b = bondStamp->getB();

        AtomStamp* atomA = getAtomStamp(a);
        AtomStamp* atomB = getAtomStamp(b);

        //find bend c--a--b
        AtomStamp::AtomIter ai;
        for(int c= atomA->getFirstBonedAtom(ai);c != -1;c = atomA->getNextBonedAtom(ai))
        {
            if(b == c)
                continue;          
            
            IntTuple3 newBend(c, a, b);
            if (newBend.first > newBend.third) {
                std::swap(newBend.first, newBend.third);
            }

            if (allBends.find(newBend) == allBends.end() ) {                
                allBends.insert(newBend);
                BendStamp * newBendStamp = new BendStamp();
                newBendStamp->setMembers(newBend);
                addBendStamp(newBendStamp);
            }
        }        

        //find bend a--b--c
        for(int c= atomB->getFirstBonedAtom(ai);c != -1;c = atomB->getNextBonedAtom(ai))
        {
            if(a == c)
                continue;          

            IntTuple3 newBend( a, b, c);
            if (newBend.first > newBend.third) {
                std::swap(newBend.first, newBend.third);
            }            
            if (allBends.find(newBend) == allBends.end() ) {                
                allBends.insert(newBend);
                BendStamp * newBendStamp = new BendStamp();
                newBendStamp->setMembers(newBend);
                addBendStamp(newBendStamp);
            }
        }        
    }

}

struct TorsionLessThan : public std::binary_function<IntTuple4, IntTuple4, bool> {
    bool operator()(IntTuple4 t1, IntTuple4 t2) {

        return t1.first < t2.first
             || (!(t2.first < t1.first) && t1.second < t2.second)
             || (!(t2.first < t1.first) && !(t2.second < t2.second) && t1.third < t2.third)
             ||(!(t2.first < t1.first) && !(t2.second < t2.second) && !(t2.third < t1.third) && t1.fourth < t2.fourth);
    }


};


void MoleculeStamp::checkTorsions() {
    for(int i = 0; i < getNTorsions(); ++i) {
        TorsionStamp* torsionStamp = getTorsionStamp(i);
        std::vector<int> torsionAtoms =  torsionStamp->getMembers();
        std::vector<int> rigidSet(getNRigidBodies(), 0);
        std::vector<int>::iterator j;
        for( j = torsionAtoms.begin(); j != torsionAtoms.end(); ++j) {
            int rigidbodyIndex = atom2Rigidbody[*j];
            if (rigidbodyIndex >= 0) {
                ++rigidSet[rigidbodyIndex];
                if (rigidSet[rigidbodyIndex] > 1) {
                    std::cout << "Error in Molecule " << getName() << ": ";
                    //std::cout << "atoms of torsion " <<  << "belong to same rigidbody " << rigidbodyIndex << "\n";                    
                }
            }
        }
    }     

    std::set<IntTuple4, TorsionLessThan> allTorsions;
    std::set<IntTuple4, TorsionLessThan>::iterator iter;
     for(int i = 0; i < getNTorsions(); ++i) {
         TorsionStamp* torsionStamp= getTorsionStamp(i);
         std::vector<int> torsion = torsionStamp->getMembers();
         if (torsion.size() == 3) {
             int ghostIndex = torsionStamp->getGhostVectorSource();
             std::vector<int>::iterator j = std::find(torsion.begin(), torsion.end(), ghostIndex);
             if (j != torsion.end()) {
                torsion.insert(j, ghostIndex);
             }
         }

        IntTuple4 torsionTuple(torsion[0], torsion[1], torsion[2], torsion[3]);
        if (torsionTuple.first > torsionTuple.fourth) {
            std::swap(torsionTuple.first, torsionTuple.fourth);
            std::swap(torsionTuple.second, torsionTuple.third);                    
        }                

         iter = allTorsions.find(torsionTuple);
         if ( iter == allTorsions.end()) {
            allTorsions.insert(torsionTuple);
         } else {
            std::cout << "Error in Molecule " << getName() << ": " << "Torsion appears multiple times\n";
         }
     }

    for (int i = 0; i < getNBonds(); ++i) {
        BondStamp* bondStamp = getBondStamp(i);
        int b = bondStamp->getA();
        int c = bondStamp->getB();

        AtomStamp* atomB = getAtomStamp(b);
        AtomStamp* atomC = getAtomStamp(c);

        AtomStamp::AtomIter ai2;
        AtomStamp::AtomIter ai3;

        for(int a = atomB->getFirstBonedAtom(ai2);a != -1;a = atomB->getNextBonedAtom(ai2))
        {
            if(a == c)
                continue;

            for(int d = atomC->getFirstBonedAtom(ai3);d != -1;d = atomC->getNextBonedAtom(ai3))
            {
                if(d == b)
                    continue;
                
                IntTuple4 newTorsion(a, b, c, d);
                //make sure the first element is always less than or equal to the fourth element in IntTuple4
                if (newTorsion.first > newTorsion.fourth) {
                    std::swap(newTorsion.first, newTorsion.fourth);
                    std::swap(newTorsion.second, newTorsion.third);                    
                }                
                if (allTorsions.find(newTorsion) == allTorsions.end() ) {                
                    allTorsions.insert(newTorsion);
                    TorsionStamp * newTorsionStamp = new TorsionStamp();
                    newTorsionStamp->setMembers(newTorsion);
                    addTorsionStamp(newTorsionStamp);                    
                }            
            }
        }    
    }

}

void MoleculeStamp::checkRigidBodies() {
     std::vector<RigidBodyStamp*>::iterator ri = std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(), static_cast<RigidBodyStamp*>(NULL));
     if (ri != rigidBodyStamps_.end()) {
         std::cout << "Error in Molecule " << getName() << ":rigidBody[" <<  ri - rigidBodyStamps_.begin()<< "] is missing\n";
     }

    for (int i = 0; i < getNRigidBodies(); ++i) {
        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
        std::vector<int> rigidAtoms =  rbStamp ->getMembers();
        std::vector<int>::iterator j =std::find_if(rigidAtoms.begin(), rigidAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
        if (j != rigidAtoms.end()) {
            std::cout << "Error in Molecule " << getName();
        }
        
    }    
}

void MoleculeStamp::checkCutoffGroups() {

    for(int i = 0; i < getNCutoffGroups(); ++i) {
        CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
        std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
        std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
        if (j != cutoffGroupAtoms.end()) {
            std::cout << "Error in Molecule " << getName() << ": cutoffGroup" << " is out of range\n"; 
        }
    }    
}

void MoleculeStamp::checkFragments() {

    std::vector<FragmentStamp*>::iterator fi = std::find(fragmentStamps_.begin(), fragmentStamps_.end(), static_cast<FragmentStamp*>(NULL));
    if (fi != fragmentStamps_.end()) {
        std::cout << "Error in Molecule " << getName() << ":fragment[" <<  fi - fragmentStamps_.begin()<< "] is missing\n";
    }
    
}

void MoleculeStamp::fillBondInfo() {

    for (int i = 0; i < getNBonds(); ++i) {
        BondStamp* bondStamp = getBondStamp(i);
        int a = bondStamp->getA();
        int b = bondStamp->getB();
        AtomStamp* atomA = getAtomStamp(a);
        AtomStamp* atomB = getAtomStamp(b);
        atomA->addBond(i);
        atomA->addBondedAtom(b);
        atomB->addBond(i);        
        atomB->addBondedAtom(a);

    }
}


//Function Name: isBondInSameRigidBody
//Return true is both atoms of the bond belong to the same rigid body, otherwise return false
bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
  int rbA;
  int rbB;
  int consAtomA;
  int consAtomB;

  if (!isAtomInRigidBody(bond->getA(),rbA, consAtomA))
    return false;

  if(!isAtomInRigidBody(bond->getB(),rbB, consAtomB) )
    return false;

  if(rbB == rbA)
    return true;
  else
    return false;
}

// Function Name: isAtomInRigidBody 
//return false if atom does not belong to a rigid body, otherwise return true 
bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
  int whichRigidBody;
  int consAtomIndex;

  return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
   
}

// Function Name: isAtomInRigidBody 
//return false if atom does not belong to a rigid body otherwise return true and set whichRigidBody 
//and consAtomIndex
//atomIndex : the index of atom in component
//whichRigidBody: the index of rigidbody in component
//consAtomIndex:  the position of joint atom apears in  rigidbody's definition
bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
  RigidBodyStamp* rbStamp;
  int numRb;
  int numAtom;

  whichRigidBody = -1;
  consAtomIndex = -1;

  numRb = this->getNRigidBodies();
  
  for(int i = 0 ; i < numRb; i++){
    rbStamp = this->getRigidBodyStamp(i);
    numAtom = rbStamp->getNMembers();
    for(int j = 0; j < numAtom; j++)
      if (rbStamp->getMemberAt(j) == atomIndex){
        whichRigidBody = i;
        consAtomIndex = j;
        return true;
      }
  }

  return false;
   
}

//return the position of joint atom apears in  rigidbody's definition
//for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
//actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
  RigidBodyStamp* rbStamp1;
  RigidBodyStamp* rbStamp2;
  int natomInRb1;
  int natomInRb2;
  int atomIndex1;
  int atomIndex2;
  std::vector<std::pair<int, int> > jointAtomIndexPair;
  
  rbStamp1 = this->getRigidBodyStamp(rb1);
  natomInRb1 =rbStamp1->getNMembers();

  rbStamp2 = this->getRigidBodyStamp(rb2);
  natomInRb2 =rbStamp2->getNMembers();

  for(int i = 0; i < natomInRb1; i++){
    atomIndex1 = rbStamp1->getMemberAt(i);
      
    for(int j= 0; j < natomInRb1; j++){
      atomIndex2 = rbStamp2->getMemberAt(j);

      if(atomIndex1 == atomIndex2){
        jointAtomIndexPair.push_back(std::make_pair(i, j));
        break;
      }
      
    }

  }

  return jointAtomIndexPair;
}

}
Revision:	2483
Committed:	Mon Dec 5 18:23:30 2005 UTC (18 years, 7 months ago) by tim
File size:	21544 byte(s)
Log Message:	adding more checking for topology
#	User	Rev	Content
1	gezelter	2204	/*
2	gezelter	1930	* 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	gezelter	1490	#include <iostream>
43	tim	2483	#include <functional>
44	tim	1492	#include "types/MoleculeStamp.hpp"
45	tim	2483	#include "utils/Tuple.hpp"
46	gezelter	1490
47	tim	2469	namespace oopse {
48			MoleculeStamp::MoleculeStamp() {
49			DefineParameter(Name, "name");
50
51			deprecatedKeywords_.insert("nAtoms");
52			deprecatedKeywords_.insert("nBonds");
53			deprecatedKeywords_.insert("nBends");
54			deprecatedKeywords_.insert("nTorsions");
55			deprecatedKeywords_.insert("nRigidBodies");
56			deprecatedKeywords_.insert("nCutoffGroups");
57
58			}
59	gezelter	1490
60	tim	2469	MoleculeStamp::~MoleculeStamp() {
61	gezelter	1490
62	tim	2469	}
63	gezelter	1490
64	tim	2469	bool MoleculeStamp::addAtomStamp( AtomStamp* atom) {
65			bool ret = addIndexSensitiveStamp(atomStamps_, atom);
66			if (!ret) {
67			std::cout << "multiple atoms have the same index: " << atom->getIndex() <<" in " << getName() << " Molecule\n";
68			}
69			return ret;
70
71			}
72	gezelter	1490
73	tim	2469	bool MoleculeStamp::addBondStamp( BondStamp* bond) {
74			bondStamps_.push_back(bond);
75			return true;
76	gezelter	1490	}
77
78	tim	2469	bool MoleculeStamp::addBendStamp( BendStamp* bend) {
79			bendStamps_.push_back(bend);
80			return true;
81			}
82	gezelter	1490
83	tim	2469	bool MoleculeStamp::addTorsionStamp( TorsionStamp* torsion) {
84			torsionStamps_.push_back(torsion);
85			return true;
86			}
87	gezelter	1490
88	tim	2469	bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp* rigidbody) {
89			bool ret = addIndexSensitiveStamp(rigidBodyStamps_, rigidbody);
90			if (!ret) {
91			std::cout << "multiple rigidbodies have the same index: " << rigidbody->getIndex() <<" in " << getName() << " Molecule\n";
92			}
93			return ret;
94			}
95	gezelter	1490
96	tim	2469	bool MoleculeStamp::addCutoffGroupStamp( CutoffGroupStamp* cutoffgroup) {
97			cutoffGroupStamps_.push_back(cutoffgroup);
98			return true;
99	gezelter	1490	}
100
101	tim	2469	bool MoleculeStamp::addFragmentStamp( FragmentStamp* fragment) {
102			return addIndexSensitiveStamp(fragmentStamps_, fragment);
103	gezelter	1490	}
104
105	tim	2469	void MoleculeStamp::validate() {
106			DataHolder::validate();
107	gezelter	1490
108	tim	2483	atom2Rigidbody.resize(getNAtoms());
109			// negative number means atom is a free atom, does not belong to rigidbody
110			//every element in atom2Rigidbody has unique negative number at the very beginning
111			for(int i = 0; i < atom2Rigidbody.size(); ++i) {
112			atom2Rigidbody[i] = -1 - i;
113			}
114			for (int i = 0; i < getNRigidBodies(); ++i) {
115			RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
116			std::vector<int> members = rbStamp->getMembers();
117			for(std::vector<int>::iterator j = members.begin(); j != members.end(); ++j) {
118			atom2Rigidbody[*j] = i;
119			}
120			}
121
122			checkAtoms();
123			checkBonds();
124			fillBondInfo();
125			checkBends();
126			checkTorsions();
127			checkRigidBodies();
128			checkCutoffGroups();
129			checkFragments();
130
131			int nrigidAtoms = 0;
132			for (int i = 0; i < getNRigidBodies(); ++i) {
133			RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
134			nrigidAtoms += rbStamp->getNMembers();
135			}
136			nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
137
138			}
139
140			void MoleculeStamp::checkAtoms() {
141	tim	2469	std::vector<AtomStamp>::iterator ai = std::find(atomStamps_.begin(), atomStamps_.end(), static_cast<AtomStamp>(NULL));
142			if (ai != atomStamps_.end()) {
143			std::cout << "Error in Molecule " << getName() << ": atom[" << ai - atomStamps_.begin()<< "] is missing\n";
144	gezelter	1490	}
145
146	tim	2469	//make sure index is not out of range
147			int natoms = getNAtoms();
148			for(int i = 0; i < getNBonds(); ++i) {
149			BondStamp* bondStamp = getBondStamp(i);
150			if (bondStamp->getA() >= natoms && bondStamp->getB() >= natoms) {
151			std::cout << "Error in Molecule " << getName() << ": bond between " << bondStamp->getA() << " and " << bondStamp->getB() << " is invalid\n";
152			}
153	gezelter	1490	}
154	tim	2483	}
155
156			void MoleculeStamp::checkBonds() {
157			//make sure bonds are unique
158			std::set<std::pair<int, int> > allBonds;
159			for(int i = 0; i < getNBonds(); ++i) {
160			BondStamp* bondStamp= getBondStamp(i);
161			std::pair<int, int> bondPair(bondStamp->getA(), bondStamp->getB());
162			//make sure bondTuple.first is always less than or equal to bondTuple.third
163			if (bondPair.first > bondPair.second) {
164			std::swap(bondPair.first, bondPair.second);
165			}
166
167			std::set<std::pair<int, int> >::iterator iter = allBonds.find(bondPair);
168			if ( iter != allBonds.end()) {
169			std::cout << "Error in Molecule " << getName() << ": " << "Bond appears multiple times\n";
170			} else {
171			allBonds.insert(bondPair);
172			}
173			}
174
175			//make sure atoms belong to same rigidbody do not bond to each other
176			for(int i = 0; i < getNBonds(); ++i) {
177			BondStamp* bondStamp = getBondStamp(i);
178			if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()]) {
179			std::cout << "Error in Molecule " << getName() << ": "<<"bond between " << bondStamp->getA() << " and " << bondStamp->getB() << " belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
180			}
181			}
182
183			}
184
185			struct BendLessThan : public std::binary_function<IntTuple4, IntTuple4, bool> {
186			bool operator()(IntTuple3 b1, IntTuple3 b2) {
187			return b1.first < b2.first
188			\|\| (!(b2.first < b1.first) && b1.second < b2.second)
189			\|\| (!(b2.first < b1.first) && !(b2.second < b2.second) && b1.third < b2.third);
190			}
191			};
192
193			void MoleculeStamp::checkBends() {
194	tim	2469	for(int i = 0; i < getNBends(); ++i) {
195			BendStamp* bendStamp = getBendStamp(i);
196			std::vector<int> bendAtoms = bendStamp->getMembers();
197	tim	2483	std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
198	tim	2469	if (j != bendAtoms.end()) {
199			std::cout << "Error in Molecule " << getName();
200			}
201	gezelter	1490
202	tim	2483	if (bendAtoms.size() == 2 ) {
203			if (!bendStamp->haveGhostVectorSource()) {
204			std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing\n";
205			}else{
206			int ghostIndex = bendStamp->getGhostVectorSource();
207			if (ghostIndex < getNAtoms()) {
208			if (std::find(bendAtoms.begin(), bendAtoms.end(), ghostIndex) == bendAtoms.end()) {
209			std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce "<< ghostIndex<<"is invalid\n";
210			}
211			if (!getAtomStamp(ghostIndex)->haveOrientation()) {
212			std::cout << "Error in Molecule " << getName() << ": ghost atom must be a directioanl atom\n";
213			}
214			}else {
215			std::cout << "Error in Molecule " << getName() << ": ghostVectorsource " << ghostIndex<< " is invalid\n";
216			}
217			}
218			} else if (bendAtoms.size() == 3 && bendStamp->haveGhostVectorSource()) {
219			std::cout << "Error in Molecule " << getName() << ": normal bend should not have ghostVectorSouce\n";
220	tim	2469	}
221	gezelter	1490	}
222
223	tim	2469	for(int i = 0; i < getNBends(); ++i) {
224			BendStamp* bendStamp = getBendStamp(i);
225			std::vector<int> bendAtoms = bendStamp->getMembers();
226			std::vector<int> rigidSet(getNRigidBodies(), 0);
227			std::vector<int>::iterator j;
228			for( j = bendAtoms.begin(); j != bendAtoms.end(); ++j) {
229			int rigidbodyIndex = atom2Rigidbody[*j];
230			if (rigidbodyIndex >= 0) {
231			++rigidSet[rigidbodyIndex];
232			if (rigidSet[rigidbodyIndex] > 1) {
233			std::cout << "Error in Molecule " << getName() << ": ";
234			//std::cout << "atoms of bend " << << "belong to same rigidbody " << rigidbodyIndex << "\n";
235			}
236			}
237			}
238	tim	2483	}
239
240
241			std::set<IntTuple3, BendLessThan> allBends;
242			std::set<IntTuple3, BendLessThan>::iterator iter;
243			for(int i = 0; i < getNBends(); ++i) {
244			BendStamp* bendStamp= getBendStamp(i);
245			std::vector<int> bend = bendStamp->getMembers();
246			if (bend.size() == 2) {
247			// in case we have two ghost bend. For example,
248			// bend {
249			// members (0, 1);
250			// ghostVectorSource = 0;
251			// }
252			// and
253			// bend {
254			// members (0, 1);
255			// ghostVectorSource = 0;
256			// }
257			// In order to distinguish them. we expand them to Tuple3.
258			// the first one is expanded to (0, 0, 1) while the second one is expaned to (0, 1, 1)
259			int ghostIndex = bendStamp->getGhostVectorSource();
260			std::vector<int>::iterator j = std::find(bend.begin(), bend.end(), ghostIndex);
261			if (j != bend.end()) {
262			bend.insert(j, ghostIndex);
263			}
264			}
265
266			IntTuple3 bendTuple(bend[0], bend[1], bend[2]);
267			//make sure bendTuple.first is always less than or equal to bendTuple.third
268			if (bendTuple.first > bendTuple.third) {
269			std::swap(bendTuple.first, bendTuple.third);
270			}
271
272			iter = allBends.find(bendTuple);
273			if ( iter != allBends.end()) {
274			std::cout << "Error in Molecule " << getName() << ": " << "Bend appears multiple times\n";
275			} else {
276			allBends.insert(bendTuple);
277			}
278			}
279
280			for (int i = 0; i < getNBonds(); ++i) {
281			BondStamp* bondStamp = getBondStamp(i);
282			int a = bondStamp->getA();
283			int b = bondStamp->getB();
284
285			AtomStamp* atomA = getAtomStamp(a);
286			AtomStamp* atomB = getAtomStamp(b);
287
288			//find bend c--a--b
289			AtomStamp::AtomIter ai;
290			for(int c= atomA->getFirstBonedAtom(ai);c != -1;c = atomA->getNextBonedAtom(ai))
291			{
292			if(b == c)
293			continue;
294
295			IntTuple3 newBend(c, a, b);
296			if (newBend.first > newBend.third) {
297			std::swap(newBend.first, newBend.third);
298			}
299
300			if (allBends.find(newBend) == allBends.end() ) {
301			allBends.insert(newBend);
302			BendStamp * newBendStamp = new BendStamp();
303			newBendStamp->setMembers(newBend);
304			addBendStamp(newBendStamp);
305			}
306			}
307
308			//find bend a--b--c
309			for(int c= atomB->getFirstBonedAtom(ai);c != -1;c = atomB->getNextBonedAtom(ai))
310			{
311			if(a == c)
312			continue;
313
314			IntTuple3 newBend( a, b, c);
315			if (newBend.first > newBend.third) {
316			std::swap(newBend.first, newBend.third);
317			}
318			if (allBends.find(newBend) == allBends.end() ) {
319			allBends.insert(newBend);
320			BendStamp * newBendStamp = new BendStamp();
321			newBendStamp->setMembers(newBend);
322			addBendStamp(newBendStamp);
323			}
324			}
325			}
326
327			}
328
329			struct TorsionLessThan : public std::binary_function<IntTuple4, IntTuple4, bool> {
330			bool operator()(IntTuple4 t1, IntTuple4 t2) {
331
332			return t1.first < t2.first
333			\|\| (!(t2.first < t1.first) && t1.second < t2.second)
334			\|\| (!(t2.first < t1.first) && !(t2.second < t2.second) && t1.third < t2.third)
335			\|\|(!(t2.first < t1.first) && !(t2.second < t2.second) && !(t2.third < t1.third) && t1.fourth < t2.fourth);
336			}
337
338
339
340			};
341
342
343			void MoleculeStamp::checkTorsions() {
344	tim	2469	for(int i = 0; i < getNTorsions(); ++i) {
345			TorsionStamp* torsionStamp = getTorsionStamp(i);
346			std::vector<int> torsionAtoms = torsionStamp->getMembers();
347			std::vector<int> rigidSet(getNRigidBodies(), 0);
348			std::vector<int>::iterator j;
349			for( j = torsionAtoms.begin(); j != torsionAtoms.end(); ++j) {
350			int rigidbodyIndex = atom2Rigidbody[*j];
351			if (rigidbodyIndex >= 0) {
352			++rigidSet[rigidbodyIndex];
353			if (rigidSet[rigidbodyIndex] > 1) {
354			std::cout << "Error in Molecule " << getName() << ": ";
355			//std::cout << "atoms of torsion " << << "belong to same rigidbody " << rigidbodyIndex << "\n";
356			}
357			}
358			}
359	tim	2483	}
360	gezelter	1490
361	tim	2483	std::set<IntTuple4, TorsionLessThan> allTorsions;
362			std::set<IntTuple4, TorsionLessThan>::iterator iter;
363			for(int i = 0; i < getNTorsions(); ++i) {
364			TorsionStamp* torsionStamp= getTorsionStamp(i);
365			std::vector<int> torsion = torsionStamp->getMembers();
366			if (torsion.size() == 3) {
367			int ghostIndex = torsionStamp->getGhostVectorSource();
368			std::vector<int>::iterator j = std::find(torsion.begin(), torsion.end(), ghostIndex);
369			if (j != torsion.end()) {
370			torsion.insert(j, ghostIndex);
371			}
372			}
373	gezelter	1490
374	tim	2483	IntTuple4 torsionTuple(torsion[0], torsion[1], torsion[2], torsion[3]);
375			if (torsionTuple.first > torsionTuple.fourth) {
376			std::swap(torsionTuple.first, torsionTuple.fourth);
377			std::swap(torsionTuple.second, torsionTuple.third);
378			}
379	gezelter	1490
380	tim	2483	iter = allTorsions.find(torsionTuple);
381			if ( iter == allTorsions.end()) {
382			allTorsions.insert(torsionTuple);
383			} else {
384			std::cout << "Error in Molecule " << getName() << ": " << "Torsion appears multiple times\n";
385			}
386			}
387
388			for (int i = 0; i < getNBonds(); ++i) {
389			BondStamp* bondStamp = getBondStamp(i);
390			int b = bondStamp->getA();
391			int c = bondStamp->getB();
392
393			AtomStamp* atomB = getAtomStamp(b);
394			AtomStamp* atomC = getAtomStamp(c);
395
396			AtomStamp::AtomIter ai2;
397			AtomStamp::AtomIter ai3;
398
399			for(int a = atomB->getFirstBonedAtom(ai2);a != -1;a = atomB->getNextBonedAtom(ai2))
400			{
401			if(a == c)
402			continue;
403
404			for(int d = atomC->getFirstBonedAtom(ai3);d != -1;d = atomC->getNextBonedAtom(ai3))
405			{
406			if(d == b)
407			continue;
408
409			IntTuple4 newTorsion(a, b, c, d);
410			//make sure the first element is always less than or equal to the fourth element in IntTuple4
411			if (newTorsion.first > newTorsion.fourth) {
412			std::swap(newTorsion.first, newTorsion.fourth);
413			std::swap(newTorsion.second, newTorsion.third);
414			}
415			if (allTorsions.find(newTorsion) == allTorsions.end() ) {
416			allTorsions.insert(newTorsion);
417			TorsionStamp * newTorsionStamp = new TorsionStamp();
418			newTorsionStamp->setMembers(newTorsion);
419			addTorsionStamp(newTorsionStamp);
420			}
421			}
422			}
423			}
424
425			}
426
427			void MoleculeStamp::checkRigidBodies() {
428			std::vector<RigidBodyStamp>::iterator ri = std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(), static_cast<RigidBodyStamp>(NULL));
429			if (ri != rigidBodyStamps_.end()) {
430			std::cout << "Error in Molecule " << getName() << ":rigidBody[" << ri - rigidBodyStamps_.begin()<< "] is missing\n";
431			}
432
433	tim	2469	for (int i = 0; i < getNRigidBodies(); ++i) {
434			RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
435	tim	2483	std::vector<int> rigidAtoms = rbStamp ->getMembers();
436			std::vector<int>::iterator j =std::find_if(rigidAtoms.begin(), rigidAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
437			if (j != rigidAtoms.end()) {
438			std::cout << "Error in Molecule " << getName();
439			}
440
441			}
442			}
443	gezelter	1490
444	tim	2483	void MoleculeStamp::checkCutoffGroups() {
445	gezelter	1490
446	tim	2483	for(int i = 0; i < getNCutoffGroups(); ++i) {
447			CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
448			std::vector<int> cutoffGroupAtoms = cutoffGroupStamp ->getMembers();
449			std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
450			if (j != cutoffGroupAtoms.end()) {
451			std::cout << "Error in Molecule " << getName() << ": cutoffGroup" << " is out of range\n";
452			}
453			}
454	gezelter	1490	}
455
456	tim	2483	void MoleculeStamp::checkFragments() {
457	gezelter	1490
458	tim	2483	std::vector<FragmentStamp>::iterator fi = std::find(fragmentStamps_.begin(), fragmentStamps_.end(), static_cast<FragmentStamp>(NULL));
459			if (fi != fragmentStamps_.end()) {
460			std::cout << "Error in Molecule " << getName() << ":fragment[" << fi - fragmentStamps_.begin()<< "] is missing\n";
461			}
462
463	gezelter	1490	}
464
465	tim	2483	void MoleculeStamp::fillBondInfo() {
466	gezelter	1490
467	tim	2483	for (int i = 0; i < getNBonds(); ++i) {
468			BondStamp* bondStamp = getBondStamp(i);
469			int a = bondStamp->getA();
470			int b = bondStamp->getB();
471			AtomStamp* atomA = getAtomStamp(a);
472			AtomStamp* atomB = getAtomStamp(b);
473			atomA->addBond(i);
474			atomA->addBondedAtom(b);
475			atomB->addBond(i);
476			atomB->addBondedAtom(a);
477
478			}
479	gezelter	1490	}
480
481	tim	2483
482
483	gezelter	1490	//Function Name: isBondInSameRigidBody
484			//Return true is both atoms of the bond belong to the same rigid body, otherwise return false
485			bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
486			int rbA;
487			int rbB;
488			int consAtomA;
489			int consAtomB;
490
491			if (!isAtomInRigidBody(bond->getA(),rbA, consAtomA))
492			return false;
493
494			if(!isAtomInRigidBody(bond->getB(),rbB, consAtomB) )
495			return false;
496
497			if(rbB == rbA)
498			return true;
499			else
500			return false;
501			}
502
503			// Function Name: isAtomInRigidBody
504			//return false if atom does not belong to a rigid body, otherwise return true
505			bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
506			int whichRigidBody;
507			int consAtomIndex;
508
509			return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
510
511			}
512
513			// Function Name: isAtomInRigidBody
514			//return false if atom does not belong to a rigid body otherwise return true and set whichRigidBody
515			//and consAtomIndex
516			//atomIndex : the index of atom in component
517			//whichRigidBody: the index of rigidbody in component
518			//consAtomIndex: the position of joint atom apears in rigidbody's definition
519			bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
520			RigidBodyStamp* rbStamp;
521			int numRb;
522			int numAtom;
523
524			whichRigidBody = -1;
525			consAtomIndex = -1;
526
527			numRb = this->getNRigidBodies();
528
529			for(int i = 0 ; i < numRb; i++){
530	tim	2469	rbStamp = this->getRigidBodyStamp(i);
531	gezelter	1490	numAtom = rbStamp->getNMembers();
532			for(int j = 0; j < numAtom; j++)
533	tim	2469	if (rbStamp->getMemberAt(j) == atomIndex){
534	gezelter	1490	whichRigidBody = i;
535			consAtomIndex = j;
536			return true;
537			}
538			}
539
540			return false;
541
542			}
543
544			//return the position of joint atom apears in rigidbody's definition
545			//for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
546			//actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
547	gezelter	2204	std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
548	gezelter	1490	RigidBodyStamp* rbStamp1;
549			RigidBodyStamp* rbStamp2;
550			int natomInRb1;
551			int natomInRb2;
552			int atomIndex1;
553			int atomIndex2;
554	gezelter	2204	std::vector<std::pair<int, int> > jointAtomIndexPair;
555	gezelter	1490
556	tim	2469	rbStamp1 = this->getRigidBodyStamp(rb1);
557	gezelter	1490	natomInRb1 =rbStamp1->getNMembers();
558
559	tim	2469	rbStamp2 = this->getRigidBodyStamp(rb2);
560	gezelter	1490	natomInRb2 =rbStamp2->getNMembers();
561
562			for(int i = 0; i < natomInRb1; i++){
563	tim	2469	atomIndex1 = rbStamp1->getMemberAt(i);
564	gezelter	1490
565			for(int j= 0; j < natomInRb1; j++){
566	tim	2469	atomIndex2 = rbStamp2->getMemberAt(j);
567	gezelter	1490
568			if(atomIndex1 == atomIndex2){
569	gezelter	1930	jointAtomIndexPair.push_back(std::make_pair(i, j));
570	gezelter	1490	break;
571			}
572
573	tim	2469	}
574	gezelter	1490
575	tim	2469	}
576	gezelter	1490
577			return jointAtomIndexPair;
578			}
579	tim	2469
580			}