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

namespace oopse {

template<class ContainerType>
std::string containerToString(ContainerType& cont) {
    std::ostringstream oss;
    oss << "(";
    typename ContainerType::iterator i = cont.begin();
    if (i != cont.end()) {
        oss << *i;
        ++i;
    }
    for (; i != cont.end();++i) {
        oss << ", ";
        oss << *i;
    }
    oss << ")";
    return oss.str();
}

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<< "Error in Molecule " << getName()  << ": multiple atoms have the same indices"<< atom->getIndex() <<"\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<< "Error in Molecule " << getName()  << ": multiple rigidbodies have the same indices: " << rigidbody->getIndex() <<"\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";
    }

}

void MoleculeStamp::checkBonds() {
    //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(" << bondStamp->getA() << ", " << bondStamp->getB() << ") is invalid\n";
        }
    }
    
    //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(" <<iter->first << ", "<< iter->second << ")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(" << bondStamp->getA() << ", " << 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() << " : atoms of bend" << containerToString(bendAtoms) << "have invalid indices\n";
        }

        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() << ": bend" << containerToString(bendAtoms) << " 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 < getNBends(); ++i) {
        TorsionStamp* torsionStamp = getTorsionStamp(i);
        std::vector<int> torsionAtoms =  torsionStamp ->getMembers();
        std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
        if (j != torsionAtoms.end()) {
            std::cout << "Error in Molecule " << getName() << ": atoms of torsion" << containerToString(torsionAtoms) << " have invalid indices\n"; 
        }
    }
    
    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() << ": torsion" << containerToString(torsionAtoms) << "is invalid\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){
  return atom2Rigidbody[atomIndex] >=0 ;
   
}

// 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){

  
  whichRigidBody = -1;
  consAtomIndex = -1;

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