src/io/RestReader.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 "io/DumpReader.hpp" 
#include "io/RestReader.hpp"
#include "primitives/Molecule.hpp"
#include "restraints/ObjectRestraint.hpp"
#include "restraints/MolecularRestraint.hpp"

namespace oopse { 
      
  void RestReader::readReferenceStructure() {

    // some of this comes directly from DumpReader.

    if (!isScanned_) 
      scanFile(); 
         
    int storageLayout = info_->getSnapshotManager()->getStorageLayout(); 
     
    if (storageLayout & DataStorage::dslPosition) { 
      needPos_ = true; 
    } else {
      needPos_ = false; 
    } 
     
    needVel_ = false;
     
    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) { 
      needQuaternion_ = true; 
    } else { 
      needQuaternion_ = false; 
    } 

    needAngMom_ = false;

    // We need temporary storage to keep track of all StuntDouble positions
    // in case some of the restraints are molecular (i.e. if they use
    // multiple SD positions to determine restrained orientations or positions:

    all_pos_.clear();
    all_pos_.resize(info_->getNGlobalIntegrableObjects() );


    // Restraint files are just standard dump files, but with the reference
    // structure stored in the first frame (frame 0).
    // RestReader overloads readSet and explicitly handles all of the
    // ObjectRestraints in that method:

    readSet(0); 
    
    // all ObjectRestraints have been handled, now we have to worry about
    // molecular restraints:

    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator j;
    Molecule* mol;
    StuntDouble* sd;

    // no need to worry about parallel molecules, as molecules are not
    // split across processor boundaries.  Just loop over all molecules
    // we know about:

    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {          

      // is this molecule restrained?    
      GenericData* data = mol->getPropertyByName("Restraint");
      
      if (data != NULL) {

        // make sure we can reinterpret the generic data as restraint data:

        RestraintData* restData= dynamic_cast<RestraintData*>(data);        

        if (restData != NULL) {

          // make sure we can reinterpet the restraint data as a
          // pointer to a MolecularRestraint:

          MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());

          if (mRest != NULL) {          

            // now we need to pack the stunt doubles for the reference
            // structure:

            std::vector<Vector3d> ref;
            int count = 0;
            RealType mass, mTot;
            Vector3d COM(0.0);
            
            mTot = 0.0;
            
            // loop over the stunt doubles in this molecule in the order we
            // will be looping them in the restraint code:
            
            for (sd = mol->beginIntegrableObject(j); sd != NULL; 
                 sd = mol->nextIntegrableObject(j)) {
              
              // push back the reference positions of the stunt
              // doubles from the *globally* sorted array of
              // positions:
              
              ref.push_back( all_pos_[sd->getGlobalIndex()] );
              
              mass = sd->getMass();
              
              COM = COM + mass * ref[count];
              mTot = mTot + mass;
              count = count + 1;
            }
            
            COM /= mTot;

            mRest->setReferenceStructure(ref, COM);         

          }
        }
      }
    }
  }

   
  void RestReader::parseDumpLine(const std::string& line) {        
    StringTokenizer tokenizer(line); 
    int nTokens; 
     
    nTokens = tokenizer.countTokens(); 
     
    if (nTokens < 2) { 
      sprintf(painCave.errMsg, 
              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); 
      painCave.isFatal = 1; 
      simError(); 
    } 

    int index = tokenizer.nextTokenAsInt();
 
    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);

    if (integrableObject == NULL) {
      return;
    }   

    std::string type = tokenizer.nextToken(); 
    int size = type.size();
    Vector3d pos;
    Vector3d vel;
    Quat4d q;
    Vector3d ji;
    Vector3d force;
    Vector3d torque;
          
    for(int i = 0; i < size; ++i) {
      switch(type[i]) {
        
        case 'p': {
            pos[0] = tokenizer.nextTokenAsDouble(); 
            pos[1] = tokenizer.nextTokenAsDouble(); 
            pos[2] = tokenizer.nextTokenAsDouble(); 
            break;
        }
        case 'v' : {
            vel[0] = tokenizer.nextTokenAsDouble(); 
            vel[1] = tokenizer.nextTokenAsDouble(); 
            vel[2] = tokenizer.nextTokenAsDouble(); 
            break;
        }

        case 'q' : {
           if (integrableObject->isDirectional()) { 
              
             q[0] = tokenizer.nextTokenAsDouble(); 
             q[1] = tokenizer.nextTokenAsDouble(); 
             q[2] = tokenizer.nextTokenAsDouble(); 
             q[3] = tokenizer.nextTokenAsDouble(); 
              
             RealType qlen = q.length(); 
             if (qlen < oopse::epsilon) { //check quaternion is not equal to 0 
                
               sprintf(painCave.errMsg, 
                       "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n"); 
               painCave.isFatal = 1; 
               simError(); 
                
             }  
              
             q.normalize(); 
           }            
           break;
        }  
        case 'j' : {
          if (integrableObject->isDirectional()) {
             ji[0] = tokenizer.nextTokenAsDouble(); 
             ji[1] = tokenizer.nextTokenAsDouble(); 
             ji[2] = tokenizer.nextTokenAsDouble(); 
          }
          break;
        }  
        case 'f': {
          force[0] = tokenizer.nextTokenAsDouble(); 
          force[1] = tokenizer.nextTokenAsDouble(); 
          force[2] = tokenizer.nextTokenAsDouble();           

          break;
        }
        case 't' : {
           torque[0] = tokenizer.nextTokenAsDouble(); 
           torque[1] = tokenizer.nextTokenAsDouble(); 
           torque[2] = tokenizer.nextTokenAsDouble();           

           break;
        }
        default: {
               sprintf(painCave.errMsg, 
                       "DumpReader Error: %s is an unrecognized type\n", type.c_str()); 
               painCave.isFatal = 1; 
               simError(); 
          break;   
        }
          
      }
    }
     
    // keep the position in case we need it for a molecular restraint:
    
    all_pos_[index] = pos;

    // is this io restrained?    
    GenericData* data = integrableObject->getPropertyByName("Restraint");
    ObjectRestraint* oRest;

    if (data != NULL) {
      // make sure we can reinterpret the generic data as restraint data:
      RestraintData* restData= dynamic_cast<RestraintData*>(data);        
      if (restData != NULL) {
        // make sure we can reinterpet the restraint data as a pointer to
        // an ObjectRestraint:
        oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
        if (oRest != NULL) {          
          if (integrableObject->isDirectional()) {
            oRest->setReferenceStructure(pos, q.toRotationMatrix3());
          } else {                           
            oRest->setReferenceStructure(pos);
          }
        }
      }
    }

  }
   

  void RestReader::readFrameProperties(std::istream& inputStream) {
    inputStream.getline(buffer, bufferSize);
    std::string line(buffer);

    if (line.find("<FrameData>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "DumpReader Error: Missing <FrameData>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    }

    // restraints don't care about frame data (unless we need to wrap
    // coordinates, but we'll worry about that later), so 
    // we'll just scan ahead until the end of the frame data:

    while(inputStream.getline(buffer, bufferSize)) {
      line = buffer;
      
      if(line.find("</FrameData>") != std::string::npos) {
        break;
      }
      
    }

  }

   
}//end namespace oopse 
Revision:	3520
Committed:	Mon Sep 7 16:31:51 2009 UTC (16 years, 2 months ago) by cli2
File size:	10381 byte(s)
Log Message:	Added new restraint infrastructure Added MolecularRestraints Added ObjectRestraints Added RestraintStamp Updated thermodynamic integration to use ObjectRestraints Added Quaternion mathematics for twist swing decompositions Significantly updated RestWriter and RestReader to use dump-like files Added selections for x, y, and z coordinates of atoms Removed monolithic Restraints class Fixed a few bugs in gradients of Euler angles in DirectionalAtom and RigidBody Added some rotational capabilities to prinicpalAxisCalculator
#	User	Rev	Content
1	cli2	3520	/*
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 "io/DumpReader.hpp"
43			#include "io/RestReader.hpp"
44			#include "primitives/Molecule.hpp"
45			#include "restraints/ObjectRestraint.hpp"
46			#include "restraints/MolecularRestraint.hpp"
47	chrisfen	2101
48	cli2	3520	namespace oopse {
49
50			void RestReader::readReferenceStructure() {
51	chrisfen	2101
52	cli2	3520	// some of this comes directly from DumpReader.
53	chrisfen	2101
54	cli2	3520	if (!isScanned_)
55			scanFile();
56
57			int storageLayout = info_->getSnapshotManager()->getStorageLayout();
58
59			if (storageLayout & DataStorage::dslPosition) {
60			needPos_ = true;
61			} else {
62			needPos_ = false;
63			}
64
65			needVel_ = false;
66
67			if (storageLayout & DataStorage::dslAmat \|\| storageLayout & DataStorage::dslElectroFrame) {
68			needQuaternion_ = true;
69			} else {
70			needQuaternion_ = false;
71			}
72	chrisfen	2101
73	cli2	3520	needAngMom_ = false;
74	chrisfen	2101
75	cli2	3520	// We need temporary storage to keep track of all StuntDouble positions
76			// in case some of the restraints are molecular (i.e. if they use
77			// multiple SD positions to determine restrained orientations or positions:
78	chrisfen	2101
79	cli2	3520	all_pos_.clear();
80			all_pos_.resize(info_->getNGlobalIntegrableObjects() );
81	chrisfen	2101
82	chrisfen	2903
83	cli2	3520	// Restraint files are just standard dump files, but with the reference
84			// structure stored in the first frame (frame 0).
85			// RestReader overloads readSet and explicitly handles all of the
86			// ObjectRestraints in that method:
87	chrisfen	2903
88	cli2	3520	readSet(0);
89
90			// all ObjectRestraints have been handled, now we have to worry about
91			// molecular restraints:
92	chrisfen	2115
93	cli2	3520	SimInfo::MoleculeIterator i;
94			Molecule::IntegrableObjectIterator j;
95			Molecule* mol;
96			StuntDouble* sd;
97	chrisfen	2903
98	cli2	3520	// no need to worry about parallel molecules, as molecules are not
99			// split across processor boundaries. Just loop over all molecules
100			// we know about:
101	chrisfen	2903
102	cli2	3520	for (mol = info_->beginMolecule(i); mol != NULL;
103			mol = info_->nextMolecule(i)) {
104	chrisfen	2903
105	cli2	3520	// is this molecule restrained?
106			GenericData* data = mol->getPropertyByName("Restraint");
107	chrisfen	2101
108	cli2	3520	if (data != NULL) {
109	chrisfen	2903
110	cli2	3520	// make sure we can reinterpret the generic data as restraint data:
111	chrisfen	2903
112	cli2	3520	RestraintData* restData= dynamic_cast<RestraintData*>(data);
113	chrisfen	2107
114	cli2	3520	if (restData != NULL) {
115	chrisfen	2115
116	cli2	3520	// make sure we can reinterpet the restraint data as a
117			// pointer to a MolecularRestraint:
118	chrisfen	2903
119	cli2	3520	MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
120	chrisfen	2903
121	cli2	3520	if (mRest != NULL) {
122	chrisfen	2115
123	cli2	3520	// now we need to pack the stunt doubles for the reference
124			// structure:
125	chrisfen	2868
126	cli2	3520	std::vector<Vector3d> ref;
127			int count = 0;
128			RealType mass, mTot;
129			Vector3d COM(0.0);
130	chrisfen	2101
131	cli2	3520	mTot = 0.0;
132	chrisfen	2101
133	cli2	3520	// loop over the stunt doubles in this molecule in the order we
134			// will be looping them in the restraint code:
135
136			for (sd = mol->beginIntegrableObject(j); sd != NULL;
137			sd = mol->nextIntegrableObject(j)) {
138
139			// push back the reference positions of the stunt
140			// doubles from the globally sorted array of
141			// positions:
142
143			ref.push_back( all_pos_[sd->getGlobalIndex()] );
144
145			mass = sd->getMass();
146
147			COM = COM + mass * ref[count];
148			mTot = mTot + mass;
149			count = count + 1;
150	chrisfen	2101	}
151
152	cli2	3520	COM /= mTot;
153	chrisfen	2868
154	cli2	3520	mRest->setReferenceStructure(ref, COM);
155	chrisfen	2868
156	chrisfen	2101	}
157			}
158			}
159			}
160			}
161	chrisfen	2115
162	cli2	3520
163
164			void RestReader::parseDumpLine(const std::string& line) {
165			StringTokenizer tokenizer(line);
166			int nTokens;
167
168			nTokens = tokenizer.countTokens();
169
170			if (nTokens < 2) {
171			sprintf(painCave.errMsg,
172			"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
173			painCave.isFatal = 1;
174			simError();
175			}
176	chrisfen	2115
177	cli2	3520	int index = tokenizer.nextTokenAsInt();
178
179			StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
180	chrisfen	2115
181	cli2	3520	if (integrableObject == NULL) {
182			return;
183			}
184	chrisfen	2115
185	cli2	3520	std::string type = tokenizer.nextToken();
186			int size = type.size();
187			Vector3d pos;
188			Vector3d vel;
189			Quat4d q;
190			Vector3d ji;
191			Vector3d force;
192			Vector3d torque;
193
194			for(int i = 0; i < size; ++i) {
195			switch(type[i]) {
196
197			case 'p': {
198			pos[0] = tokenizer.nextTokenAsDouble();
199			pos[1] = tokenizer.nextTokenAsDouble();
200			pos[2] = tokenizer.nextTokenAsDouble();
201			break;
202			}
203			case 'v' : {
204			vel[0] = tokenizer.nextTokenAsDouble();
205			vel[1] = tokenizer.nextTokenAsDouble();
206			vel[2] = tokenizer.nextTokenAsDouble();
207			break;
208			}
209	chrisfen	2115
210	cli2	3520	case 'q' : {
211			if (integrableObject->isDirectional()) {
212
213			q[0] = tokenizer.nextTokenAsDouble();
214			q[1] = tokenizer.nextTokenAsDouble();
215			q[2] = tokenizer.nextTokenAsDouble();
216			q[3] = tokenizer.nextTokenAsDouble();
217
218			RealType qlen = q.length();
219			if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
220
221			sprintf(painCave.errMsg,
222			"DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
223			painCave.isFatal = 1;
224			simError();
225
226			}
227
228			q.normalize();
229			}
230			break;
231			}
232			case 'j' : {
233			if (integrableObject->isDirectional()) {
234			ji[0] = tokenizer.nextTokenAsDouble();
235			ji[1] = tokenizer.nextTokenAsDouble();
236			ji[2] = tokenizer.nextTokenAsDouble();
237			}
238			break;
239			}
240			case 'f': {
241			force[0] = tokenizer.nextTokenAsDouble();
242			force[1] = tokenizer.nextTokenAsDouble();
243			force[2] = tokenizer.nextTokenAsDouble();
244	chrisfen	2903
245	cli2	3520	break;
246			}
247			case 't' : {
248			torque[0] = tokenizer.nextTokenAsDouble();
249			torque[1] = tokenizer.nextTokenAsDouble();
250			torque[2] = tokenizer.nextTokenAsDouble();
251
252			break;
253			}
254			default: {
255			sprintf(painCave.errMsg,
256			"DumpReader Error: %s is an unrecognized type\n", type.c_str());
257			painCave.isFatal = 1;
258			simError();
259			break;
260			}
261
262	chrisfen	2101	}
263			}
264	cli2	3520
265			// keep the position in case we need it for a molecular restraint:
266	chrisfen	2101
267	cli2	3520	all_pos_[index] = pos;
268	chrisfen	2903
269	cli2	3520	// is this io restrained?
270			GenericData* data = integrableObject->getPropertyByName("Restraint");
271			ObjectRestraint* oRest;
272	chrisfen	2903
273	cli2	3520	if (data != NULL) {
274			// make sure we can reinterpret the generic data as restraint data:
275			RestraintData* restData= dynamic_cast<RestraintData*>(data);
276			if (restData != NULL) {
277			// make sure we can reinterpet the restraint data as a pointer to
278			// an ObjectRestraint:
279			oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
280			if (oRest != NULL) {
281			if (integrableObject->isDirectional()) {
282			oRest->setReferenceStructure(pos, q.toRotationMatrix3());
283			} else {
284			oRest->setReferenceStructure(pos);
285			}
286			}
287	chrisfen	2101	}
288			}
289	chrisfen	2107
290	cli2	3520	}
291
292	chrisfen	2903
293	chrisfen	2107
294	cli2	3520	void RestReader::readFrameProperties(std::istream& inputStream) {
295			inputStream.getline(buffer, bufferSize);
296			std::string line(buffer);
297	chrisfen	2107
298	cli2	3520	if (line.find("<FrameData>") == std::string::npos) {
299			sprintf(painCave.errMsg,
300			"DumpReader Error: Missing <FrameData>\n");
301			painCave.isFatal = 1;
302			simError();
303			}
304	chrisfen	2903
305	cli2	3520	// restraints don't care about frame data (unless we need to wrap
306			// coordinates, but we'll worry about that later), so
307			// we'll just scan ahead until the end of the frame data:
308
309			while(inputStream.getline(buffer, bufferSize)) {
310			line = buffer;
311	chrisfen	2101
312	cli2	3520	if(line.find("</FrameData>") != std::string::npos) {
313			break;
314	chrisfen	2101	}
315
316	cli2	3520	}
317	chrisfen	2814
318	chrisfen	2101	}
319	chrisfen	2992
320	cli2	3520
321			}//end namespace oopse