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:	1360
Committed:	Mon Sep 7 16:31:51 2009 UTC (15 years, 7 months ago) by cli2
Original Path:	*trunk/src/io/RestReader.cpp*
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	1360	/*
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	417
48	cli2	1360	namespace oopse {
49
50			void RestReader::readReferenceStructure() {
51	chrisfen	417
52	cli2	1360	// some of this comes directly from DumpReader.
53	chrisfen	417
54	cli2	1360	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	417
73	cli2	1360	needAngMom_ = false;
74	chrisfen	417
75	cli2	1360	// 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	417
79	cli2	1360	all_pos_.clear();
80			all_pos_.resize(info_->getNGlobalIntegrableObjects() );
81	chrisfen	417
82	chrisfen	996
83	cli2	1360	// 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	996
88	cli2	1360	readSet(0);
89
90			// all ObjectRestraints have been handled, now we have to worry about
91			// molecular restraints:
92	chrisfen	431
93	cli2	1360	SimInfo::MoleculeIterator i;
94			Molecule::IntegrableObjectIterator j;
95			Molecule* mol;
96			StuntDouble* sd;
97	chrisfen	996
98	cli2	1360	// 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	996
102	cli2	1360	for (mol = info_->beginMolecule(i); mol != NULL;
103			mol = info_->nextMolecule(i)) {
104	chrisfen	996
105	cli2	1360	// is this molecule restrained?
106			GenericData* data = mol->getPropertyByName("Restraint");
107	chrisfen	417
108	cli2	1360	if (data != NULL) {
109	chrisfen	996
110	cli2	1360	// make sure we can reinterpret the generic data as restraint data:
111	chrisfen	996
112	cli2	1360	RestraintData* restData= dynamic_cast<RestraintData*>(data);
113	chrisfen	423
114	cli2	1360	if (restData != NULL) {
115	chrisfen	431
116	cli2	1360	// make sure we can reinterpet the restraint data as a
117			// pointer to a MolecularRestraint:
118	chrisfen	996
119	cli2	1360	MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
120	chrisfen	996
121	cli2	1360	if (mRest != NULL) {
122	chrisfen	431
123	cli2	1360	// now we need to pack the stunt doubles for the reference
124			// structure:
125	chrisfen	990
126	cli2	1360	std::vector<Vector3d> ref;
127			int count = 0;
128			RealType mass, mTot;
129			Vector3d COM(0.0);
130	chrisfen	417
131	cli2	1360	mTot = 0.0;
132	chrisfen	417
133	cli2	1360	// 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	417	}
151
152	cli2	1360	COM /= mTot;
153	chrisfen	990
154	cli2	1360	mRest->setReferenceStructure(ref, COM);
155	chrisfen	990
156	chrisfen	417	}
157			}
158			}
159			}
160			}
161	chrisfen	431
162	cli2	1360
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	431
177	cli2	1360	int index = tokenizer.nextTokenAsInt();
178
179			StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
180	chrisfen	431
181	cli2	1360	if (integrableObject == NULL) {
182			return;
183			}
184	chrisfen	431
185	cli2	1360	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	431
210	cli2	1360	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	996
245	cli2	1360	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	417	}
263			}
264	cli2	1360
265			// keep the position in case we need it for a molecular restraint:
266	chrisfen	417
267	cli2	1360	all_pos_[index] = pos;
268	chrisfen	996
269	cli2	1360	// is this io restrained?
270			GenericData* data = integrableObject->getPropertyByName("Restraint");
271			ObjectRestraint* oRest;
272	chrisfen	996
273	cli2	1360	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	417	}
288			}
289	chrisfen	423
290	cli2	1360	}
291
292	chrisfen	996
293	chrisfen	423
294	cli2	1360	void RestReader::readFrameProperties(std::istream& inputStream) {
295			inputStream.getline(buffer, bufferSize);
296			std::string line(buffer);
297	chrisfen	423
298	cli2	1360	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	996
305	cli2	1360	// 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	417
312	cli2	1360	if(line.find("</FrameData>") != std::string::npos) {
313			break;
314	chrisfen	417	}
315
316	cli2	1360	}
317	chrisfen	985
318	chrisfen	417	}
319	chrisfen	1030
320	cli2	1360
321			}//end namespace oopse