src/restraints/RestraintForceManager.cpp

/*
 * Copyright (c) 2009 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 <cmath>
#include "restraints/RestraintForceManager.hpp"
#include "restraints/MolecularRestraint.hpp"
#include "restraints/ObjectRestraint.hpp"
#include "io/RestReader.hpp"
#include "utils/simError.h"
#include "utils/OOPSEConstant.hpp"
#include "utils/StringUtils.hpp"
#include "selection/SelectionEvaluator.hpp"
#include "selection/SelectionManager.hpp"
#ifdef IS_MPI
#include <mpi.h>
#endif


namespace oopse {

  RestraintForceManager::RestraintForceManager(SimInfo* info): ForceManager(info) {

    // order of affairs:
    //
    // 1) create restraints from the restraintStamps found in the MD
    // file.
    //
    // 2) Create RestraintReader to parse the input files for the ideal
    // structures.  This reader will set reference structures, and will 
    // calculate molecular centers of mass, etc.
    //
    // 3) sit around and wait for calcForces to be called.  When it comes, 
    // call the normal force manager calcForces, then loop through the 
    // restrained objects and do their restraint forces.

    currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
    Globals* simParam = info_->getSimParams();

    if (simParam->haveStatusTime()){      
      restTime_ = simParam->getStatusTime();
    } else {
      sprintf(painCave.errMsg,
              "Restraint warning: If you use restraints without setting\n"
              "\tstatusTime, no restraint data will be written to the rest\n"
              "\tfile.\n");
      painCave.isFatal = 0;
      simError();
      restTime_ = simParam->getRunTime();
    }

    int nRestraintStamps = simParam->getNRestraintStamps();
    std::vector<RestraintStamp*> stamp = simParam->getRestraintStamps();

    for (int i = 0; i < nRestraintStamps; i++){

      std::string myType = toUpperCopy(stamp[i]->getType());

      if (myType.compare("MOLECULAR")==0){

        int molIndex;
        std::vector<Vector3d> ref;
        Vector3d refCom;

        if (!stamp[i]->haveMolIndex()) {
          sprintf(painCave.errMsg,
                  "Restraint Error: A molecular restraint was specified\n"
                  "\twithout providing a value for molIndex.\n");
          painCave.isFatal = 1;
          simError();      
        } else {
          molIndex = stamp[i]->getMolIndex();
        }
        
        Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex);
        
        if (mol == NULL) {
          sprintf(painCave.errMsg,
                  "Restraint Error: A molecular restraint was specified, but\n"
                  "\tno molecule was found with global index %d.\n",
                  molIndex);
          painCave.isFatal = 1;
          simError();      
        }
        
        MolecularRestraint* rest = new MolecularRestraint();

        std::string restPre("mol_");
        std::stringstream restName;
        restName << restPre << molIndex;
        rest->setRestraintName(restName.str());
        
        if (stamp[i]->haveDisplacementSpringConstant()) {
          rest->setDisplacementForceConstant(stamp[i]->getDisplacementSpringConstant());
        }                  
        if (stamp[i]->haveTwistSpringConstant()) {
          rest->setTwistForceConstant(stamp[i]->getTwistSpringConstant());
        }        
        if (stamp[i]->haveSwingXSpringConstant()) {
          rest->setSwingXForceConstant(stamp[i]->getSwingXSpringConstant());
        }
        if (stamp[i]->haveSwingYSpringConstant()) {
          rest->setSwingYForceConstant(stamp[i]->getSwingYSpringConstant());
        }
        if (stamp[i]->haveRestrainedTwistAngle()) {
          rest->setRestrainedTwistAngle(stamp[i]->getRestrainedTwistAngle() * M_PI/180.0);
        }
        if (stamp[i]->haveRestrainedSwingYAngle()) {
          rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle() * M_PI/180.0);
        }
        if (stamp[i]->haveRestrainedSwingXAngle()) {
          rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle() * M_PI/180.0);
        }
        if (stamp[i]->havePrint()) {
          rest->setPrintRestraint(stamp[i]->getPrint());
        }

        restraints_.push_back(rest);
        mol->addProperty(new RestraintData("Restraint", rest));
        restrainedMols_.push_back(mol);
        
      } else if (myType.compare("OBJECT") == 0) {
        
        std::string objectSelection;

        if (!stamp[i]->haveObjectSelection()) {
          sprintf(painCave.errMsg,
                  "Restraint Error: An object restraint was specified\n"
                  "\twithout providing a selection script in the\n"
                  "\tobjectSelection variable.\n");
          painCave.isFatal = 1;
          simError();      
        } else {
          objectSelection = stamp[i]->getObjectSelection();
        }

        SelectionEvaluator evaluator(info);
        SelectionManager seleMan(info);
                
        evaluator.loadScriptString(objectSelection);
        seleMan.setSelectionSet(evaluator.evaluate());        
        int selectionCount = seleMan.getSelectionCount();
        
        sprintf(painCave.errMsg,
                "Restraint Info: The specified restraint objectSelection,\n"
                "\t\t%s\n"
                "\twill result in %d integrable objects being\n"
                "\trestrained.\n", objectSelection.c_str(), selectionCount);
        painCave.isFatal = 0;
        simError();                     

        int selei;
        StuntDouble* sd;

        for (sd = seleMan.beginSelected(selei); sd != NULL; 
             sd = seleMan.nextSelected(selei)) {
          
          ObjectRestraint* rest = new ObjectRestraint();
          
          if (stamp[i]->haveDisplacementSpringConstant()) {
            rest->setDisplacementForceConstant(stamp[i]->getDisplacementSpringConstant());
          }                  
          if (stamp[i]->haveTwistSpringConstant()) {
            rest->setTwistForceConstant(stamp[i]->getTwistSpringConstant());
          }        
          if (stamp[i]->haveSwingXSpringConstant()) {
            rest->setSwingXForceConstant(stamp[i]->getSwingXSpringConstant());
          }
          if (stamp[i]->haveSwingYSpringConstant()) {
            rest->setSwingYForceConstant(stamp[i]->getSwingYSpringConstant());
          }
          if (stamp[i]->haveRestrainedTwistAngle()) {
            rest->setRestrainedTwistAngle(stamp[i]->getRestrainedTwistAngle());
          }
          if (stamp[i]->haveRestrainedSwingXAngle()) {
            rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle());
          }
          if (stamp[i]->haveRestrainedSwingYAngle()) {
            rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle());
          }     
          if (stamp[i]->havePrint()) {
            rest->setPrintRestraint(stamp[i]->getPrint());
          }
    
          restraints_.push_back(rest);
          sd->addProperty(new RestraintData("Restraint", rest));
          restrainedObjs_.push_back(sd);                    
        }

      }
    }

    // ThermodynamicIntegration subclasses RestraintForceManager, and there
    // are times when it won't use restraints at all, so only open the
    // restraint file if we are actually using restraints:

    if (simParam->getUseRestraints()) {
      std::string refFile = simParam->getRestraint_file();
      RestReader* rr = new RestReader(info, refFile);

      rr->readReferenceStructure();
    }

    restOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".rest";
    restOut = new RestWriter(info_, restOutput_.c_str(), restraints_);
    
    if(!restOut){
      sprintf(painCave.errMsg, "Restraint error: Failed to create RestWriter\n");
      painCave.isFatal = 1;
      simError();
    } 

    // todo: figure out the scale factor.  Right now, just scale them all to 1
    std::vector<Restraint*>::const_iterator resti;
    for(resti=restraints_.begin(); resti != restraints_.end(); ++resti){
      (*resti)->setScaleFactor(1.0);
    }
  }

  RestraintForceManager::~RestraintForceManager(){
    if (restOut)
      delete restOut;
  }

  void RestraintForceManager::init() {
    currRestTime_ = currSnapshot_->getTime();
  }

  void RestraintForceManager::calcForces(bool needPotential, bool needStress){

    ForceManager::calcForces(needPotential, needStress);    
    RealType restPot_local, restPot;

    restPot_local = doRestraints(1.0);

#ifdef IS_MPI    
    MPI::COMM_WORLD.Allreduce(&restPot_local, &restPot, 1, 
                              MPI::REALTYPE, MPI::SUM);
#else
    restPot = restPot_local;
#endif
    currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
    currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += restPot;
    currSnapshot_->statData[Stats::VHARM] = restPot;

    //write out forces and current positions of restrained molecules    
    if (currSnapshot_->getTime() >= currRestTime_){
      restOut->writeRest(restInfo_);
      currRestTime_ += restTime_;
    }
  }

  RealType RestraintForceManager::doRestraints(RealType scalingFactor){
    std::vector<Molecule*>::const_iterator rm;
    GenericData* data;
    Molecule::IntegrableObjectIterator ioi;
    MolecularRestraint* mRest;
    StuntDouble* sd;
    RealType pTot;

    std::vector<StuntDouble*>::const_iterator ro;
    ObjectRestraint* oRest;

    std::map<int, Restraint::RealPair> restInfo;

    unscaledPotential_ = 0.0;

    restInfo_.clear();

    for(rm=restrainedMols_.begin(); rm != restrainedMols_.end(); ++rm){

      // make sure this molecule (*rm) has a generic data for restraints:
      data = (*rm)->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
          // an MolecularRestraint:
          mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
          if (mRest == NULL) {
            sprintf( painCave.errMsg,
                     "Can not cast RestraintData to MolecularRestraint\n");
            painCave.severity = OOPSE_ERROR;
            painCave.isFatal = 1;
            simError();                      
          }
        } else {
          sprintf( painCave.errMsg,
                   "Can not cast GenericData to RestraintData\n");
          painCave.severity = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();       
        }
      } else {
        sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();          
      }

      // phew.  At this point, we should have the pointer to the
      // correct MolecularRestraint in the variable mRest.

      Vector3d molCom = (*rm)->getCom();

      std::vector<Vector3d> struc;
      std::vector<Vector3d> forces;
      
      for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL; 
          sd = (*rm)->nextIntegrableObject(ioi)) {        
        struc.push_back(sd->getPos());
      }

      mRest->setScaleFactor(scalingFactor);
      mRest->calcForce(struc, molCom);
      forces = mRest->getRestraintForces();
      int index = 0;

      for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL; 
          sd = (*rm)->nextIntegrableObject(ioi)) {        
        sd->addFrc(forces[index]);
        struc.push_back(sd->getPos());
        index++;
      }
      
      unscaledPotential_ += mRest->getUnscaledPotential();

      restInfo = mRest->getRestraintInfo();

      // only collect data on restraints that we're going to print:
      if (mRest->getPrintRestraint()) 
        restInfo_.push_back(restInfo);
    }

    for(ro=restrainedObjs_.begin(); ro != restrainedObjs_.end(); ++ro){
      // make sure this object (*ro) has a generic data for restraints:
      data = (*ro)->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
          // an ObjectRestraint:
          oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
          if (oRest == NULL) {
            sprintf( painCave.errMsg,
                     "Can not cast RestraintData to ObjectRestraint\n");
            painCave.severity = OOPSE_ERROR;
            painCave.isFatal = 1;
            simError();                      
          }
        } else {
          sprintf( painCave.errMsg,
                   "Can not cast GenericData to RestraintData\n");
          painCave.severity = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();       
        }
      } else {
        sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();          
      }
      
      // phew.  At this point, we should have the pointer to the
      // correct Object restraint in the variable oRest.

      oRest->setScaleFactor(scalingFactor);
      
      Vector3d pos = (*ro)->getPos();

      if ( (*ro)->isDirectional() ) {

        // directional objects may have orientational restraints as well
        // as positional, so get the rotation matrix first:

        RotMat3x3d A = (*ro)->getA();
        oRest->calcForce(pos, A);
        (*ro)->addFrc(oRest->getRestraintForce());
        (*ro)->addTrq(oRest->getRestraintTorque());
      } else {        

        // plain vanilla positional restraints:

        oRest->calcForce(pos);
        (*ro)->addFrc(oRest->getRestraintForce());
      }

      unscaledPotential_ += oRest->getUnscaledPotential();

      restInfo = oRest->getRestraintInfo();

      // only collect data on restraints that we're going to print:
      if (oRest->getPrintRestraint()) 
        restInfo_.push_back(restInfo);
    }

    return unscaledPotential_ * scalingFactor;
  }
}
Revision:	1364
Committed:	Wed Oct 7 20:49:50 2009 UTC (15 years, 6 months ago) by cli2
File size:	16182 byte(s)
Log Message:	Bug fixes in the SelectionEvaluator and SelectionCompiler Added print option in Restraints
#	User	Rev	Content
1	cli2	1360	/*
2			* Copyright (c) 2009 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 <cmath>
43			#include "restraints/RestraintForceManager.hpp"
44			#include "restraints/MolecularRestraint.hpp"
45			#include "restraints/ObjectRestraint.hpp"
46			#include "io/RestReader.hpp"
47			#include "utils/simError.h"
48			#include "utils/OOPSEConstant.hpp"
49			#include "utils/StringUtils.hpp"
50			#include "selection/SelectionEvaluator.hpp"
51			#include "selection/SelectionManager.hpp"
52			#ifdef IS_MPI
53			#include <mpi.h>
54			#endif
55
56
57			namespace oopse {
58
59			RestraintForceManager::RestraintForceManager(SimInfo* info): ForceManager(info) {
60
61			// order of affairs:
62			//
63			// 1) create restraints from the restraintStamps found in the MD
64			// file.
65			//
66			// 2) Create RestraintReader to parse the input files for the ideal
67			// structures. This reader will set reference structures, and will
68			// calculate molecular centers of mass, etc.
69			//
70			// 3) sit around and wait for calcForces to be called. When it comes,
71			// call the normal force manager calcForces, then loop through the
72			// restrained objects and do their restraint forces.
73
74			currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
75			Globals* simParam = info_->getSimParams();
76
77			if (simParam->haveStatusTime()){
78			restTime_ = simParam->getStatusTime();
79			} else {
80			sprintf(painCave.errMsg,
81	cli2	1364	"Restraint warning: If you use restraints without setting\n"
82			"\tstatusTime, no restraint data will be written to the rest\n"
83	cli2	1360	"\tfile.\n");
84			painCave.isFatal = 0;
85			simError();
86			restTime_ = simParam->getRunTime();
87			}
88
89			int nRestraintStamps = simParam->getNRestraintStamps();
90			std::vector<RestraintStamp*> stamp = simParam->getRestraintStamps();
91
92			for (int i = 0; i < nRestraintStamps; i++){
93
94			std::string myType = toUpperCopy(stamp[i]->getType());
95
96			if (myType.compare("MOLECULAR")==0){
97
98			int molIndex;
99			std::vector<Vector3d> ref;
100			Vector3d refCom;
101
102			if (!stamp[i]->haveMolIndex()) {
103			sprintf(painCave.errMsg,
104			"Restraint Error: A molecular restraint was specified\n"
105			"\twithout providing a value for molIndex.\n");
106			painCave.isFatal = 1;
107			simError();
108			} else {
109			molIndex = stamp[i]->getMolIndex();
110			}
111
112			Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex);
113
114			if (mol == NULL) {
115			sprintf(painCave.errMsg,
116			"Restraint Error: A molecular restraint was specified, but\n"
117			"\tno molecule was found with global index %d.\n",
118			molIndex);
119			painCave.isFatal = 1;
120			simError();
121			}
122
123			MolecularRestraint* rest = new MolecularRestraint();
124
125			std::string restPre("mol_");
126			std::stringstream restName;
127			restName << restPre << molIndex;
128			rest->setRestraintName(restName.str());
129
130			if (stamp[i]->haveDisplacementSpringConstant()) {
131			rest->setDisplacementForceConstant(stamp[i]->getDisplacementSpringConstant());
132			}
133			if (stamp[i]->haveTwistSpringConstant()) {
134			rest->setTwistForceConstant(stamp[i]->getTwistSpringConstant());
135			}
136			if (stamp[i]->haveSwingXSpringConstant()) {
137			rest->setSwingXForceConstant(stamp[i]->getSwingXSpringConstant());
138			}
139			if (stamp[i]->haveSwingYSpringConstant()) {
140			rest->setSwingYForceConstant(stamp[i]->getSwingYSpringConstant());
141			}
142			if (stamp[i]->haveRestrainedTwistAngle()) {
143			rest->setRestrainedTwistAngle(stamp[i]->getRestrainedTwistAngle() * M_PI/180.0);
144			}
145			if (stamp[i]->haveRestrainedSwingYAngle()) {
146			rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle() * M_PI/180.0);
147			}
148			if (stamp[i]->haveRestrainedSwingXAngle()) {
149			rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle() * M_PI/180.0);
150			}
151	cli2	1364	if (stamp[i]->havePrint()) {
152			rest->setPrintRestraint(stamp[i]->getPrint());
153			}
154	cli2	1360
155			restraints_.push_back(rest);
156			mol->addProperty(new RestraintData("Restraint", rest));
157			restrainedMols_.push_back(mol);
158
159			} else if (myType.compare("OBJECT") == 0) {
160
161			std::string objectSelection;
162
163			if (!stamp[i]->haveObjectSelection()) {
164			sprintf(painCave.errMsg,
165			"Restraint Error: An object restraint was specified\n"
166			"\twithout providing a selection script in the\n"
167			"\tobjectSelection variable.\n");
168			painCave.isFatal = 1;
169			simError();
170			} else {
171			objectSelection = stamp[i]->getObjectSelection();
172			}
173
174			SelectionEvaluator evaluator(info);
175			SelectionManager seleMan(info);
176	cli2	1364
177	cli2	1360	evaluator.loadScriptString(objectSelection);
178			seleMan.setSelectionSet(evaluator.evaluate());
179			int selectionCount = seleMan.getSelectionCount();
180
181			sprintf(painCave.errMsg,
182			"Restraint Info: The specified restraint objectSelection,\n"
183			"\t\t%s\n"
184			"\twill result in %d integrable objects being\n"
185			"\trestrained.\n", objectSelection.c_str(), selectionCount);
186			painCave.isFatal = 0;
187			simError();
188
189			int selei;
190			StuntDouble* sd;
191
192			for (sd = seleMan.beginSelected(selei); sd != NULL;
193			sd = seleMan.nextSelected(selei)) {
194
195			ObjectRestraint* rest = new ObjectRestraint();
196
197			if (stamp[i]->haveDisplacementSpringConstant()) {
198			rest->setDisplacementForceConstant(stamp[i]->getDisplacementSpringConstant());
199			}
200			if (stamp[i]->haveTwistSpringConstant()) {
201			rest->setTwistForceConstant(stamp[i]->getTwistSpringConstant());
202			}
203			if (stamp[i]->haveSwingXSpringConstant()) {
204			rest->setSwingXForceConstant(stamp[i]->getSwingXSpringConstant());
205			}
206			if (stamp[i]->haveSwingYSpringConstant()) {
207			rest->setSwingYForceConstant(stamp[i]->getSwingYSpringConstant());
208			}
209			if (stamp[i]->haveRestrainedTwistAngle()) {
210			rest->setRestrainedTwistAngle(stamp[i]->getRestrainedTwistAngle());
211			}
212			if (stamp[i]->haveRestrainedSwingXAngle()) {
213			rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle());
214			}
215			if (stamp[i]->haveRestrainedSwingYAngle()) {
216			rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle());
217	cli2	1364	}
218			if (stamp[i]->havePrint()) {
219			rest->setPrintRestraint(stamp[i]->getPrint());
220			}
221
222	cli2	1360	restraints_.push_back(rest);
223			sd->addProperty(new RestraintData("Restraint", rest));
224			restrainedObjs_.push_back(sd);
225			}
226
227			}
228			}
229
230			// ThermodynamicIntegration subclasses RestraintForceManager, and there
231			// are times when it won't use restraints at all, so only open the
232			// restraint file if we are actually using restraints:
233
234			if (simParam->getUseRestraints()) {
235			std::string refFile = simParam->getRestraint_file();
236			RestReader* rr = new RestReader(info, refFile);
237
238			rr->readReferenceStructure();
239			}
240
241			restOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".rest";
242			restOut = new RestWriter(info_, restOutput_.c_str(), restraints_);
243
244			if(!restOut){
245			sprintf(painCave.errMsg, "Restraint error: Failed to create RestWriter\n");
246			painCave.isFatal = 1;
247			simError();
248			}
249
250			// todo: figure out the scale factor. Right now, just scale them all to 1
251			std::vector<Restraint*>::const_iterator resti;
252			for(resti=restraints_.begin(); resti != restraints_.end(); ++resti){
253			(*resti)->setScaleFactor(1.0);
254			}
255			}
256
257			RestraintForceManager::~RestraintForceManager(){
258			if (restOut)
259			delete restOut;
260			}
261
262			void RestraintForceManager::init() {
263			currRestTime_ = currSnapshot_->getTime();
264			}
265
266			void RestraintForceManager::calcForces(bool needPotential, bool needStress){
267
268			ForceManager::calcForces(needPotential, needStress);
269			RealType restPot_local, restPot;
270
271			restPot_local = doRestraints(1.0);
272
273			#ifdef IS_MPI
274			MPI::COMM_WORLD.Allreduce(&restPot_local, &restPot, 1,
275			MPI::REALTYPE, MPI::SUM);
276			#else
277			restPot = restPot_local;
278			#endif
279			currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
280			currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += restPot;
281			currSnapshot_->statData[Stats::VHARM] = restPot;
282
283			//write out forces and current positions of restrained molecules
284			if (currSnapshot_->getTime() >= currRestTime_){
285			restOut->writeRest(restInfo_);
286			currRestTime_ += restTime_;
287			}
288			}
289
290			RealType RestraintForceManager::doRestraints(RealType scalingFactor){
291			std::vector<Molecule*>::const_iterator rm;
292			GenericData* data;
293			Molecule::IntegrableObjectIterator ioi;
294			MolecularRestraint* mRest;
295			StuntDouble* sd;
296			RealType pTot;
297
298			std::vector<StuntDouble*>::const_iterator ro;
299			ObjectRestraint* oRest;
300
301			std::map<int, Restraint::RealPair> restInfo;
302
303			unscaledPotential_ = 0.0;
304
305			restInfo_.clear();
306
307			for(rm=restrainedMols_.begin(); rm != restrainedMols_.end(); ++rm){
308
309			// make sure this molecule (*rm) has a generic data for restraints:
310			data = (*rm)->getPropertyByName("Restraint");
311			if (data != NULL) {
312			// make sure we can reinterpret the generic data as restraint data:
313			RestraintData* restData= dynamic_cast<RestraintData*>(data);
314			if (restData != NULL) {
315			// make sure we can reinterpet the restraint data as a pointer to
316			// an MolecularRestraint:
317			mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
318			if (mRest == NULL) {
319			sprintf( painCave.errMsg,
320			"Can not cast RestraintData to MolecularRestraint\n");
321			painCave.severity = OOPSE_ERROR;
322			painCave.isFatal = 1;
323			simError();
324			}
325			} else {
326			sprintf( painCave.errMsg,
327			"Can not cast GenericData to RestraintData\n");
328			painCave.severity = OOPSE_ERROR;
329			painCave.isFatal = 1;
330			simError();
331			}
332			} else {
333			sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
334			painCave.severity = OOPSE_ERROR;
335			painCave.isFatal = 1;
336			simError();
337			}
338
339			// phew. At this point, we should have the pointer to the
340			// correct MolecularRestraint in the variable mRest.
341
342			Vector3d molCom = (*rm)->getCom();
343
344			std::vector<Vector3d> struc;
345			std::vector<Vector3d> forces;
346
347			for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
348			sd = (*rm)->nextIntegrableObject(ioi)) {
349			struc.push_back(sd->getPos());
350			}
351
352			mRest->setScaleFactor(scalingFactor);
353			mRest->calcForce(struc, molCom);
354			forces = mRest->getRestraintForces();
355			int index = 0;
356
357			for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
358			sd = (*rm)->nextIntegrableObject(ioi)) {
359			sd->addFrc(forces[index]);
360			struc.push_back(sd->getPos());
361			index++;
362			}
363
364			unscaledPotential_ += mRest->getUnscaledPotential();
365
366			restInfo = mRest->getRestraintInfo();
367	cli2	1364
368			// only collect data on restraints that we're going to print:
369			if (mRest->getPrintRestraint())
370			restInfo_.push_back(restInfo);
371	cli2	1360	}
372
373			for(ro=restrainedObjs_.begin(); ro != restrainedObjs_.end(); ++ro){
374			// make sure this object (*ro) has a generic data for restraints:
375			data = (*ro)->getPropertyByName("Restraint");
376			if (data != NULL) {
377			// make sure we can reinterpret the generic data as restraint data:
378			RestraintData* restData= dynamic_cast<RestraintData*>(data);
379			if (restData != NULL) {
380			// make sure we can reinterpet the restraint data as a pointer to
381			// an ObjectRestraint:
382			oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
383			if (oRest == NULL) {
384			sprintf( painCave.errMsg,
385			"Can not cast RestraintData to ObjectRestraint\n");
386			painCave.severity = OOPSE_ERROR;
387			painCave.isFatal = 1;
388			simError();
389			}
390			} else {
391			sprintf( painCave.errMsg,
392			"Can not cast GenericData to RestraintData\n");
393			painCave.severity = OOPSE_ERROR;
394			painCave.isFatal = 1;
395			simError();
396			}
397			} else {
398			sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
399			painCave.severity = OOPSE_ERROR;
400			painCave.isFatal = 1;
401			simError();
402			}
403
404			// phew. At this point, we should have the pointer to the
405			// correct Object restraint in the variable oRest.
406
407			oRest->setScaleFactor(scalingFactor);
408
409			Vector3d pos = (*ro)->getPos();
410
411			if ( (*ro)->isDirectional() ) {
412
413			// directional objects may have orientational restraints as well
414			// as positional, so get the rotation matrix first:
415
416			RotMat3x3d A = (*ro)->getA();
417			oRest->calcForce(pos, A);
418			(*ro)->addFrc(oRest->getRestraintForce());
419			(*ro)->addTrq(oRest->getRestraintTorque());
420			} else {
421
422			// plain vanilla positional restraints:
423
424			oRest->calcForce(pos);
425			(*ro)->addFrc(oRest->getRestraintForce());
426			}
427
428			unscaledPotential_ += oRest->getUnscaledPotential();
429
430			restInfo = oRest->getRestraintInfo();
431	cli2	1364
432			// only collect data on restraints that we're going to print:
433			if (oRest->getPrintRestraint())
434			restInfo_.push_back(restInfo);
435	cli2	1360	}
436
437			return unscaledPotential_ * scalingFactor;
438			}
439			}