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);
        }

        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());
          }         
          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();
      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();
      restInfo_.push_back(restInfo);
    }

    return unscaledPotential_ * scalingFactor;
  }
}
Revision:	1360
Committed:	Mon Sep 7 16:31:51 2009 UTC (15 years, 7 months ago) by cli2
File size:	15741 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) 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			"Restraint warning: If you use restraints without setting\n",
82			"\tstatusTime, no restraint data will be written to the rest\n",
83			"\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
152			restraints_.push_back(rest);
153			mol->addProperty(new RestraintData("Restraint", rest));
154			restrainedMols_.push_back(mol);
155
156			} else if (myType.compare("OBJECT") == 0) {
157
158			std::string objectSelection;
159
160			if (!stamp[i]->haveObjectSelection()) {
161			sprintf(painCave.errMsg,
162			"Restraint Error: An object restraint was specified\n"
163			"\twithout providing a selection script in the\n"
164			"\tobjectSelection variable.\n");
165			painCave.isFatal = 1;
166			simError();
167			} else {
168			objectSelection = stamp[i]->getObjectSelection();
169			}
170
171			SelectionEvaluator evaluator(info);
172			SelectionManager seleMan(info);
173
174			evaluator.loadScriptString(objectSelection);
175			seleMan.setSelectionSet(evaluator.evaluate());
176			int selectionCount = seleMan.getSelectionCount();
177
178			sprintf(painCave.errMsg,
179			"Restraint Info: The specified restraint objectSelection,\n"
180			"\t\t%s\n"
181			"\twill result in %d integrable objects being\n"
182			"\trestrained.\n", objectSelection.c_str(), selectionCount);
183			painCave.isFatal = 0;
184			simError();
185
186			int selei;
187			StuntDouble* sd;
188
189			for (sd = seleMan.beginSelected(selei); sd != NULL;
190			sd = seleMan.nextSelected(selei)) {
191
192			ObjectRestraint* rest = new ObjectRestraint();
193
194			if (stamp[i]->haveDisplacementSpringConstant()) {
195			rest->setDisplacementForceConstant(stamp[i]->getDisplacementSpringConstant());
196			}
197			if (stamp[i]->haveTwistSpringConstant()) {
198			rest->setTwistForceConstant(stamp[i]->getTwistSpringConstant());
199			}
200			if (stamp[i]->haveSwingXSpringConstant()) {
201			rest->setSwingXForceConstant(stamp[i]->getSwingXSpringConstant());
202			}
203			if (stamp[i]->haveSwingYSpringConstant()) {
204			rest->setSwingYForceConstant(stamp[i]->getSwingYSpringConstant());
205			}
206			if (stamp[i]->haveRestrainedTwistAngle()) {
207			rest->setRestrainedTwistAngle(stamp[i]->getRestrainedTwistAngle());
208			}
209			if (stamp[i]->haveRestrainedSwingXAngle()) {
210			rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle());
211			}
212			if (stamp[i]->haveRestrainedSwingYAngle()) {
213			rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle());
214			}
215			restraints_.push_back(rest);
216			sd->addProperty(new RestraintData("Restraint", rest));
217			restrainedObjs_.push_back(sd);
218			}
219
220			}
221			}
222
223			// ThermodynamicIntegration subclasses RestraintForceManager, and there
224			// are times when it won't use restraints at all, so only open the
225			// restraint file if we are actually using restraints:
226
227			if (simParam->getUseRestraints()) {
228			std::string refFile = simParam->getRestraint_file();
229			RestReader* rr = new RestReader(info, refFile);
230
231			rr->readReferenceStructure();
232			}
233
234			restOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".rest";
235			restOut = new RestWriter(info_, restOutput_.c_str(), restraints_);
236
237			if(!restOut){
238			sprintf(painCave.errMsg, "Restraint error: Failed to create RestWriter\n");
239			painCave.isFatal = 1;
240			simError();
241			}
242
243			// todo: figure out the scale factor. Right now, just scale them all to 1
244			std::vector<Restraint*>::const_iterator resti;
245			for(resti=restraints_.begin(); resti != restraints_.end(); ++resti){
246			(*resti)->setScaleFactor(1.0);
247			}
248			}
249
250			RestraintForceManager::~RestraintForceManager(){
251			if (restOut)
252			delete restOut;
253			}
254
255			void RestraintForceManager::init() {
256			currRestTime_ = currSnapshot_->getTime();
257			}
258
259			void RestraintForceManager::calcForces(bool needPotential, bool needStress){
260
261			ForceManager::calcForces(needPotential, needStress);
262			RealType restPot_local, restPot;
263
264			restPot_local = doRestraints(1.0);
265
266			#ifdef IS_MPI
267			MPI::COMM_WORLD.Allreduce(&restPot_local, &restPot, 1,
268			MPI::REALTYPE, MPI::SUM);
269			#else
270			restPot = restPot_local;
271			#endif
272			currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
273			currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += restPot;
274			currSnapshot_->statData[Stats::VHARM] = restPot;
275
276			//write out forces and current positions of restrained molecules
277			if (currSnapshot_->getTime() >= currRestTime_){
278			restOut->writeRest(restInfo_);
279			currRestTime_ += restTime_;
280			}
281			}
282
283			RealType RestraintForceManager::doRestraints(RealType scalingFactor){
284			std::vector<Molecule*>::const_iterator rm;
285			GenericData* data;
286			Molecule::IntegrableObjectIterator ioi;
287			MolecularRestraint* mRest;
288			StuntDouble* sd;
289			RealType pTot;
290
291			std::vector<StuntDouble*>::const_iterator ro;
292			ObjectRestraint* oRest;
293
294			std::map<int, Restraint::RealPair> restInfo;
295
296			unscaledPotential_ = 0.0;
297
298			restInfo_.clear();
299
300			for(rm=restrainedMols_.begin(); rm != restrainedMols_.end(); ++rm){
301
302			// make sure this molecule (*rm) has a generic data for restraints:
303			data = (*rm)->getPropertyByName("Restraint");
304			if (data != NULL) {
305			// make sure we can reinterpret the generic data as restraint data:
306			RestraintData* restData= dynamic_cast<RestraintData*>(data);
307			if (restData != NULL) {
308			// make sure we can reinterpet the restraint data as a pointer to
309			// an MolecularRestraint:
310			mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
311			if (mRest == NULL) {
312			sprintf( painCave.errMsg,
313			"Can not cast RestraintData to MolecularRestraint\n");
314			painCave.severity = OOPSE_ERROR;
315			painCave.isFatal = 1;
316			simError();
317			}
318			} else {
319			sprintf( painCave.errMsg,
320			"Can not cast GenericData to RestraintData\n");
321			painCave.severity = OOPSE_ERROR;
322			painCave.isFatal = 1;
323			simError();
324			}
325			} else {
326			sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
327			painCave.severity = OOPSE_ERROR;
328			painCave.isFatal = 1;
329			simError();
330			}
331
332			// phew. At this point, we should have the pointer to the
333			// correct MolecularRestraint in the variable mRest.
334
335			Vector3d molCom = (*rm)->getCom();
336
337			std::vector<Vector3d> struc;
338			std::vector<Vector3d> forces;
339
340			for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
341			sd = (*rm)->nextIntegrableObject(ioi)) {
342			struc.push_back(sd->getPos());
343			}
344
345			mRest->setScaleFactor(scalingFactor);
346			mRest->calcForce(struc, molCom);
347			forces = mRest->getRestraintForces();
348			int index = 0;
349
350			for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
351			sd = (*rm)->nextIntegrableObject(ioi)) {
352			sd->addFrc(forces[index]);
353			struc.push_back(sd->getPos());
354			index++;
355			}
356
357			unscaledPotential_ += mRest->getUnscaledPotential();
358
359			restInfo = mRest->getRestraintInfo();
360			restInfo_.push_back(restInfo);
361			}
362
363			for(ro=restrainedObjs_.begin(); ro != restrainedObjs_.end(); ++ro){
364			// make sure this object (*ro) has a generic data for restraints:
365			data = (*ro)->getPropertyByName("Restraint");
366			if (data != NULL) {
367			// make sure we can reinterpret the generic data as restraint data:
368			RestraintData* restData= dynamic_cast<RestraintData*>(data);
369			if (restData != NULL) {
370			// make sure we can reinterpet the restraint data as a pointer to
371			// an ObjectRestraint:
372			oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
373			if (oRest == NULL) {
374			sprintf( painCave.errMsg,
375			"Can not cast RestraintData to ObjectRestraint\n");
376			painCave.severity = OOPSE_ERROR;
377			painCave.isFatal = 1;
378			simError();
379			}
380			} else {
381			sprintf( painCave.errMsg,
382			"Can not cast GenericData to RestraintData\n");
383			painCave.severity = OOPSE_ERROR;
384			painCave.isFatal = 1;
385			simError();
386			}
387			} else {
388			sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
389			painCave.severity = OOPSE_ERROR;
390			painCave.isFatal = 1;
391			simError();
392			}
393
394			// phew. At this point, we should have the pointer to the
395			// correct Object restraint in the variable oRest.
396
397			oRest->setScaleFactor(scalingFactor);
398
399			Vector3d pos = (*ro)->getPos();
400
401			if ( (*ro)->isDirectional() ) {
402
403			// directional objects may have orientational restraints as well
404			// as positional, so get the rotation matrix first:
405
406			RotMat3x3d A = (*ro)->getA();
407			oRest->calcForce(pos, A);
408			(*ro)->addFrc(oRest->getRestraintForce());
409			(*ro)->addTrq(oRest->getRestraintTorque());
410			} else {
411
412			// plain vanilla positional restraints:
413
414			oRest->calcForce(pos);
415			(*ro)->addFrc(oRest->getRestraintForce());
416			}
417
418			unscaledPotential_ += oRest->getUnscaledPotential();
419
420			restInfo = oRest->getRestraintInfo();
421			restInfo_.push_back(restInfo);
422			}
423
424			return unscaledPotential_ * scalingFactor;
425			}
426			}