applications/simpleBuilder/MoLocator.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 <cstdlib>
#include <cmath>

#include "utils/simError.h"
#include "applications/simpleBuilder/MoLocator.hpp"
#include "types/AtomType.hpp"
namespace oopse {
MoLocator::MoLocator( MoleculeStamp* theStamp, ForceField* theFF){

  myStamp = theStamp;
  myFF = theFF;
  nIntegrableObjects = myStamp->getNIntegrable();
  calcRef();
}

void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
    Vector3d newCoor;
    Vector3d curRefCoor;  
    RotMat3x3d rotMat = latVec2RotMat(ort);
  
    if(mol->getNIntegrableObjects() != nIntegrableObjects){
        sprintf( painCave.errMsg,
            "MoLocator error.\n"
            "  The number of integrable objects of MoleculeStamp is not the same as  that of Molecule\n");
        painCave.isFatal = 1;
        simError();
    }

    Molecule::IntegrableObjectIterator ii;
    StuntDouble* integrableObject;
    int i;
    for (integrableObject = mol->beginIntegrableObject(ii), i = 0; integrableObject != NULL; 
        integrableObject = mol->nextIntegrableObject(ii), ++i) { 

        newCoor = rotMat * refCoords[i];
        newCoor += offset;

        integrableObject->setPos( newCoor);
        integrableObject->setVel(V3Zero);

        if(integrableObject->isDirectional()){
          integrableObject->setA(rotMat * integrableObject->getA());
          integrableObject->setJ(V3Zero);  
        }        
    }
}
 
void MoLocator::calcRef( void ){
  AtomStamp* currAtomStamp;
  int nAtoms; 
  int nRigidBodies;
   std::vector<double> mass;
  Vector3d coor;
  Vector3d refMolCom;  
  int nAtomsInRb;
  double totMassInRb;
  double currAtomMass;
  double molMass;
  
  nAtoms= myStamp->getNAtoms();
  nRigidBodies = myStamp->getNRigidBodies();

  for(size_t i=0; i<nAtoms; i++){

    currAtomStamp = myStamp->getAtom(i);

    if( !currAtomStamp->havePosition() ){
      sprintf( painCave.errMsg,
                  "MoLocator error.\n"
                  "  Component %s, atom %s does not have a position specified.\n"
                  "  This means MoLocator cannot initalize it's position.\n",
                  myStamp->getID(),
                  currAtomStamp->getType() );

      painCave.isFatal = 1;
      simError();
    }

    //if atom belongs to rigidbody, just skip it
    if(myStamp->isAtomInRigidBody(i))
      continue;
    //get mass and the reference coordinate 
    else{
    
      mass.push_back(currAtomMass);
      coor.x() = currAtomStamp->getPosX();
      coor.y() = currAtomStamp->getPosY();
      coor.z() = currAtomStamp->getPosZ();
      refCoords.push_back(coor);

    }
  }

  for(int i = 0; i < nRigidBodies; i++){
    coor.x() = 0;
    coor.y() = 0;
    coor.z() = 0;
    totMassInRb = 0;

    for(int j = 0; j < nAtomsInRb; j++){

      currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
      totMassInRb +=  currAtomMass;
      
      coor.x() += currAtomStamp->getPosX() * currAtomMass;
      coor.y() += currAtomStamp->getPosY() * currAtomMass;
      coor.z() += currAtomStamp->getPosZ() * currAtomMass;
    }

    mass.push_back(totMassInRb);
    coor /= totMassInRb;
    refCoords.push_back(coor);
  }


  //calculate the reference center of mass
  molMass = 0;
  refMolCom.x() = 0;
  refMolCom.y() = 0;
  refMolCom.z() = 0;
  
  for(int i = 0; i < nIntegrableObjects; i++){
    refMolCom += refCoords[i] * mass[i];
   molMass += mass[i];
  }
  
  refMolCom /= molMass;

  //move the reference center of mass to (0,0,0) and adjust the reference coordinate 
  //of the integrabel objects
  for(int i = 0; i < nIntegrableObjects; i++)
    refCoords[i] -= refMolCom;
}


double getAtomMass(const std::string& at, ForceField* myFF) {
    double mass;
    AtomType* atomType= myFF->getAtomType(at);
    if (atomType != NULL) {
        mass =     atomType->getMass();
    } else {
        mass = 0.0;
        std::cerr << "Can not find AtomType: " << at << std::endl;
    }
    return mass;
}

double getMolMass(MoleculeStamp *molStamp, ForceField *myFF) {
    int nAtoms;
    double totMass = 0;
    nAtoms = molStamp->getNAtoms();

    for(size_t i = 0; i < nAtoms; i++) {
        AtomStamp *currAtomStamp = molStamp->getAtom(i);
        totMass += getAtomMass(currAtomStamp->getType(), myFF);         
    }
    return totMass;
}
RotMat3x3d latVec2RotMat(const Vector3d& lv){

    double theta =acos(lv[2]);
    double phi = atan2(lv[1], lv[0]);
    double psi = 0;

    return RotMat3x3d(phi, theta, psi);

}
}

Revision:	1930
Committed:	Wed Jan 12 22:41:40 2005 UTC (19 years, 7 months ago) by gezelter
File size:	6544 byte(s)
Log Message:	merging new_design branch into OOPSE-2.0
#	User	Rev	Content
1	gezelter	1930	/*
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	tim	1501	#include <iostream>
43
44			#include <cstdlib>
45			#include <cmath>
46
47			#include "utils/simError.h"
48			#include "applications/simpleBuilder/MoLocator.hpp"
49	gezelter	1930	#include "types/AtomType.hpp"
50			namespace oopse {
51			MoLocator::MoLocator( MoleculeStamp* theStamp, ForceField* theFF){
52	tim	1501
53			myStamp = theStamp;
54			myFF = theFF;
55			nIntegrableObjects = myStamp->getNIntegrable();
56	gezelter	1930	calcRef();
57	tim	1501	}
58
59			void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
60	gezelter	1930	Vector3d newCoor;
61			Vector3d curRefCoor;
62			RotMat3x3d rotMat = latVec2RotMat(ort);
63	tim	1501
64	gezelter	1930	if(mol->getNIntegrableObjects() != nIntegrableObjects){
65			sprintf( painCave.errMsg,
66			"MoLocator error.\n"
67			" The number of integrable objects of MoleculeStamp is not the same as that of Molecule\n");
68			painCave.isFatal = 1;
69			simError();
70			}
71	tim	1501
72	gezelter	1930	Molecule::IntegrableObjectIterator ii;
73			StuntDouble* integrableObject;
74			int i;
75			for (integrableObject = mol->beginIntegrableObject(ii), i = 0; integrableObject != NULL;
76			integrableObject = mol->nextIntegrableObject(ii), ++i) {
77	tim	1501
78	gezelter	1930	newCoor = rotMat * refCoords[i];
79			newCoor += offset;
80	tim	1501
81	gezelter	1930	integrableObject->setPos( newCoor);
82			integrableObject->setVel(V3Zero);
83	tim	1501
84	gezelter	1930	if(integrableObject->isDirectional()){
85			integrableObject->setA(rotMat * integrableObject->getA());
86			integrableObject->setJ(V3Zero);
87			}
88	tim	1501	}
89			}
90
91	gezelter	1930	void MoLocator::calcRef( void ){
92	tim	1501	AtomStamp* currAtomStamp;
93			int nAtoms;
94			int nRigidBodies;
95	gezelter	1930	std::vector<double> mass;
96	tim	1501	Vector3d coor;
97			Vector3d refMolCom;
98			int nAtomsInRb;
99			double totMassInRb;
100			double currAtomMass;
101			double molMass;
102
103			nAtoms= myStamp->getNAtoms();
104			nRigidBodies = myStamp->getNRigidBodies();
105
106			for(size_t i=0; i<nAtoms; i++){
107
108			currAtomStamp = myStamp->getAtom(i);
109
110			if( !currAtomStamp->havePosition() ){
111			sprintf( painCave.errMsg,
112			"MoLocator error.\n"
113			" Component %s, atom %s does not have a position specified.\n"
114			" This means MoLocator cannot initalize it's position.\n",
115			myStamp->getID(),
116			currAtomStamp->getType() );
117
118			painCave.isFatal = 1;
119			simError();
120			}
121
122			//if atom belongs to rigidbody, just skip it
123			if(myStamp->isAtomInRigidBody(i))
124			continue;
125			//get mass and the reference coordinate
126			else{
127	gezelter	1930
128	tim	1501	mass.push_back(currAtomMass);
129	tim	1681	coor.x() = currAtomStamp->getPosX();
130			coor.y() = currAtomStamp->getPosY();
131			coor.z() = currAtomStamp->getPosZ();
132	tim	1501	refCoords.push_back(coor);
133
134			}
135			}
136
137			for(int i = 0; i < nRigidBodies; i++){
138	tim	1681	coor.x() = 0;
139			coor.y() = 0;
140			coor.z() = 0;
141	tim	1501	totMassInRb = 0;
142
143			for(int j = 0; j < nAtomsInRb; j++){
144
145	gezelter	1930	currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
146	tim	1501	totMassInRb += currAtomMass;
147
148	tim	1681	coor.x() += currAtomStamp->getPosX() * currAtomMass;
149			coor.y() += currAtomStamp->getPosY() * currAtomMass;
150			coor.z() += currAtomStamp->getPosZ() * currAtomMass;
151	tim	1501	}
152
153			mass.push_back(totMassInRb);
154			coor /= totMassInRb;
155			refCoords.push_back(coor);
156			}
157
158
159			//calculate the reference center of mass
160			molMass = 0;
161	tim	1681	refMolCom.x() = 0;
162			refMolCom.y() = 0;
163			refMolCom.z() = 0;
164	tim	1501
165			for(int i = 0; i < nIntegrableObjects; i++){
166			refMolCom += refCoords[i] * mass[i];
167			molMass += mass[i];
168			}
169
170			refMolCom /= molMass;
171
172			//move the reference center of mass to (0,0,0) and adjust the reference coordinate
173			//of the integrabel objects
174			for(int i = 0; i < nIntegrableObjects; i++)
175			refCoords[i] -= refMolCom;
176			}
177
178
179
180	gezelter	1930	double getAtomMass(const std::string& at, ForceField* myFF) {
181			double mass;
182			AtomType* atomType= myFF->getAtomType(at);
183			if (atomType != NULL) {
184			mass = atomType->getMass();
185			} else {
186			mass = 0.0;
187			std::cerr << "Can not find AtomType: " << at << std::endl;
188			}
189			return mass;
190			}
191	tim	1501
192	gezelter	1930	double getMolMass(MoleculeStamp molStamp, ForceField myFF) {
193			int nAtoms;
194			double totMass = 0;
195			nAtoms = molStamp->getNAtoms();
196
197			for(size_t i = 0; i < nAtoms; i++) {
198			AtomStamp *currAtomStamp = molStamp->getAtom(i);
199			totMass += getAtomMass(currAtomStamp->getType(), myFF);
200			}
201			return totMass;
202	tim	1501	}
203	gezelter	1930	RotMat3x3d latVec2RotMat(const Vector3d& lv){
204	tim	1501
205	gezelter	1930	double theta =acos(lv[2]);
206			double phi = atan2(lv[1], lv[0]);
207			double psi = 0;
208
209			return RotMat3x3d(phi, theta, psi);
210
211			}
212			}
213