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