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. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. 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.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
 
#ifdef IS_MPI
#include <mpi.h>
#endif

#include "config.h"
#include <cmath>

#include "restraints/RestraintForceManager.hpp"
#include "restraints/MolecularRestraint.hpp"
#include "restraints/ObjectRestraint.hpp"
#include "io/RestReader.hpp"
#include "utils/simError.h"
#include "utils/PhysicalConstants.hpp"
#include "utils/StringUtils.hpp"
#include "selection/SelectionEvaluator.hpp"
#include "selection/SelectionManager.hpp"

namespace OpenMD {

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

    std::vector<int> stuntDoubleIndex;

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

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

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

        int molIndex;
        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();
        }
        
        if (molIndex < 0) {
          sprintf(painCave.errMsg,
                  "Restraint Error: A molecular restraint was specified\n"
                  "\twith a molIndex that was less than 0\n");
          painCave.isFatal = 1;
          simError();      
        }
        if (molIndex >= info_->getNGlobalMolecules()) {
          sprintf(painCave.errMsg,
                  "Restraint Error: A molecular restraint was specified with\n"
                  "\ta molIndex that was greater than the total number of molecules\n");
          painCave.isFatal = 1;
          simError();      
        }
       
        Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex);
        
        if (mol == NULL) {
#ifdef IS_MPI
          // getMoleculeByGlobalIndex returns a NULL in parallel if
          // this proc doesn't have the molecule.  Do a quick check to
          // make sure another processor is supposed to have it.
         
          int myrank;
          MPI_Comm_rank( MPI_COMM_WORLD, &myrank);

          if (info_->getMolToProc(molIndex) == myrank) {
         
            // If we were supposed to have it but got a null, then freak out.
#endif

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

#ifdef IS_MPI
          }
#endif
        }
        

#ifdef IS_MPI
        // only handle this molecular restraint if this processor owns the
        // molecule
        int myrank;
        MPI_Comm_rank( MPI_COMM_WORLD, &myrank);
        if (info_->getMolToProc(molIndex) == myrank) {

#endif

        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);
#ifdef IS_MPI
        }
#endif        
      } 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();

#ifdef IS_MPI    
        MPI_Allreduce(MPI_IN_PLACE, &selectionCount, 1, MPI_INT, MPI_SUM, 
                      MPI_COMM_WORLD);
#endif
                
        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.severity = OPENMD_INFO;
        painCave.isFatal = 0;
        simError();                     

        int selei;
        StuntDouble* sd;

        for (sd = seleMan.beginSelected(selei); sd != NULL; 
             sd = seleMan.nextSelected(selei)) {
          stuntDoubleIndex.push_back(sd->getGlobalIntegrableObjectIndex());

          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, stuntDoubleIndex);
      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(){

    ForceManager::calcForces();    
    RealType restPot(0.0);

    restPot = doRestraints(1.0);

#ifdef IS_MPI    
    MPI_Allreduce(MPI_IN_PLACE, &restPot, 1, MPI_REALTYPE, MPI_SUM, 
                  MPI_COMM_WORLD);
#endif

    currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
    currSnapshot_->setRestraintPotential(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;

    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 = OPENMD_ERROR;
            painCave.isFatal = 1;
            simError();                      
          }
        } else {
          sprintf( painCave.errMsg,
                   "Can not cast GenericData to RestraintData\n");
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError();       
        }
      } else {
        sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
        painCave.severity = OPENMD_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();      

      // only collect data on restraints that we're going to print:
      if (mRest->getPrintRestraint()) 
        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 = OPENMD_ERROR;
            painCave.isFatal = 1;
            simError();                      
          }
        } else {
          sprintf( painCave.errMsg,
                   "Can not cast GenericData to RestraintData\n");
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError();       
        }
      } else {
        sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
        painCave.severity = OPENMD_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();

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

    return unscaledPotential_ * scalingFactor;
  }
}
Revision:	2024
Committed:	Thu Oct 16 19:13:51 2014 UTC (10 years, 6 months ago) by gezelter
File size:	17771 byte(s)
Log Message:	Added Radial and Z-projected velocity autocorrelation functions Started to add SequentialProps program Mucking about with angular restraint potentials
#	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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	cli2	1360	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	cli2	1360	* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31	gezelter	1390	*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	gezelter	1879	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	gezelter	1782	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	cli2	1360	*/
42
43	gezelter	1938	#ifdef IS_MPI
44			#include <mpi.h>
45			#endif
46
47	gezelter	1782	#include "config.h"
48	cli2	1360	#include <cmath>
49	gezelter	1782
50	cli2	1360	#include "restraints/RestraintForceManager.hpp"
51			#include "restraints/MolecularRestraint.hpp"
52			#include "restraints/ObjectRestraint.hpp"
53			#include "io/RestReader.hpp"
54			#include "utils/simError.h"
55	gezelter	1390	#include "utils/PhysicalConstants.hpp"
56	cli2	1360	#include "utils/StringUtils.hpp"
57			#include "selection/SelectionEvaluator.hpp"
58			#include "selection/SelectionManager.hpp"
59
60	gezelter	1390	namespace OpenMD {
61	cli2	1360
62			RestraintForceManager::RestraintForceManager(SimInfo* info): ForceManager(info) {
63
64			// order of affairs:
65			//
66			// 1) create restraints from the restraintStamps found in the MD
67			// file.
68			//
69			// 2) Create RestraintReader to parse the input files for the ideal
70			// structures. This reader will set reference structures, and will
71			// calculate molecular centers of mass, etc.
72			//
73			// 3) sit around and wait for calcForces to be called. When it comes,
74			// call the normal force manager calcForces, then loop through the
75			// restrained objects and do their restraint forces.
76
77			currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
78			Globals* simParam = info_->getSimParams();
79
80			if (simParam->haveStatusTime()){
81			restTime_ = simParam->getStatusTime();
82			} else {
83			sprintf(painCave.errMsg,
84	cli2	1364	"Restraint warning: If you use restraints without setting\n"
85			"\tstatusTime, no restraint data will be written to the rest\n"
86	cli2	1360	"\tfile.\n");
87			painCave.isFatal = 0;
88			simError();
89			restTime_ = simParam->getRunTime();
90			}
91
92			int nRestraintStamps = simParam->getNRestraintStamps();
93			std::vector<RestraintStamp*> stamp = simParam->getRestraintStamps();
94
95	cli2	1407	std::vector<int> stuntDoubleIndex;
96
97	cli2	1360	for (int i = 0; i < nRestraintStamps; i++){
98
99			std::string myType = toUpperCopy(stamp[i]->getType());
100
101			if (myType.compare("MOLECULAR")==0){
102
103			int molIndex;
104			Vector3d refCom;
105
106			if (!stamp[i]->haveMolIndex()) {
107			sprintf(painCave.errMsg,
108			"Restraint Error: A molecular restraint was specified\n"
109			"\twithout providing a value for molIndex.\n");
110			painCave.isFatal = 1;
111			simError();
112			} else {
113			molIndex = stamp[i]->getMolIndex();
114			}
115
116	gezelter	1366	if (molIndex < 0) {
117	cli2	1360	sprintf(painCave.errMsg,
118	gezelter	1366	"Restraint Error: A molecular restraint was specified\n"
119			"\twith a molIndex that was less than 0\n");
120	cli2	1360	painCave.isFatal = 1;
121			simError();
122			}
123	gezelter	1366	if (molIndex >= info_->getNGlobalMolecules()) {
124			sprintf(painCave.errMsg,
125			"Restraint Error: A molecular restraint was specified with\n"
126			"\ta molIndex that was greater than the total number of molecules\n");
127			painCave.isFatal = 1;
128			simError();
129			}
130
131			Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex);
132	cli2	1360
133	gezelter	1366	if (mol == NULL) {
134			#ifdef IS_MPI
135			// getMoleculeByGlobalIndex returns a NULL in parallel if
136			// this proc doesn't have the molecule. Do a quick check to
137			// make sure another processor is supposed to have it.
138
139	gezelter	1969	int myrank;
140			MPI_Comm_rank( MPI_COMM_WORLD, &myrank);
141
142	gezelter	1366	if (info_->getMolToProc(molIndex) == myrank) {
143	cli2	1407
144	gezelter	1366	// If we were supposed to have it but got a null, then freak out.
145			#endif
146
147			sprintf(painCave.errMsg,
148			"Restraint Error: A molecular restraint was specified, but\n"
149			"\tno molecule was found with global index %d.\n",
150			molIndex);
151			painCave.isFatal = 1;
152			simError();
153
154			#ifdef IS_MPI
155			}
156			#endif
157			}
158
159	gezelter	1398
160			#ifdef IS_MPI
161			// only handle this molecular restraint if this processor owns the
162			// molecule
163	gezelter	1969	int myrank;
164			MPI_Comm_rank( MPI_COMM_WORLD, &myrank);
165	gezelter	1398	if (info_->getMolToProc(molIndex) == myrank) {
166	cli2	1407
167	gezelter	1398	#endif
168
169	cli2	1360	MolecularRestraint* rest = new MolecularRestraint();
170
171			std::string restPre("mol_");
172			std::stringstream restName;
173			restName << restPre << molIndex;
174			rest->setRestraintName(restName.str());
175
176			if (stamp[i]->haveDisplacementSpringConstant()) {
177			rest->setDisplacementForceConstant(stamp[i]->getDisplacementSpringConstant());
178			}
179			if (stamp[i]->haveTwistSpringConstant()) {
180			rest->setTwistForceConstant(stamp[i]->getTwistSpringConstant());
181			}
182			if (stamp[i]->haveSwingXSpringConstant()) {
183			rest->setSwingXForceConstant(stamp[i]->getSwingXSpringConstant());
184			}
185			if (stamp[i]->haveSwingYSpringConstant()) {
186			rest->setSwingYForceConstant(stamp[i]->getSwingYSpringConstant());
187			}
188			if (stamp[i]->haveRestrainedTwistAngle()) {
189			rest->setRestrainedTwistAngle(stamp[i]->getRestrainedTwistAngle() * M_PI/180.0);
190			}
191			if (stamp[i]->haveRestrainedSwingYAngle()) {
192			rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle() * M_PI/180.0);
193			}
194			if (stamp[i]->haveRestrainedSwingXAngle()) {
195			rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle() * M_PI/180.0);
196			}
197	cli2	1364	if (stamp[i]->havePrint()) {
198			rest->setPrintRestraint(stamp[i]->getPrint());
199			}
200	cli2	1360
201			restraints_.push_back(rest);
202			mol->addProperty(new RestraintData("Restraint", rest));
203			restrainedMols_.push_back(mol);
204	gezelter	1398	#ifdef IS_MPI
205			}
206			#endif
207	cli2	1360	} else if (myType.compare("OBJECT") == 0) {
208
209			std::string objectSelection;
210
211			if (!stamp[i]->haveObjectSelection()) {
212			sprintf(painCave.errMsg,
213			"Restraint Error: An object restraint was specified\n"
214			"\twithout providing a selection script in the\n"
215			"\tobjectSelection variable.\n");
216			painCave.isFatal = 1;
217			simError();
218			} else {
219			objectSelection = stamp[i]->getObjectSelection();
220			}
221
222			SelectionEvaluator evaluator(info);
223			SelectionManager seleMan(info);
224	cli2	1364
225	cli2	1360	evaluator.loadScriptString(objectSelection);
226			seleMan.setSelectionSet(evaluator.evaluate());
227			int selectionCount = seleMan.getSelectionCount();
228	gezelter	2020
229			#ifdef IS_MPI
230			MPI_Allreduce(MPI_IN_PLACE, &selectionCount, 1, MPI_INT, MPI_SUM,
231			MPI_COMM_WORLD);
232			#endif
233
234	cli2	1360	sprintf(painCave.errMsg,
235			"Restraint Info: The specified restraint objectSelection,\n"
236			"\t\t%s\n"
237			"\twill result in %d integrable objects being\n"
238			"\trestrained.\n", objectSelection.c_str(), selectionCount);
239	gezelter	1973	painCave.severity = OPENMD_INFO;
240	cli2	1360	painCave.isFatal = 0;
241			simError();
242
243			int selei;
244			StuntDouble* sd;
245
246			for (sd = seleMan.beginSelected(selei); sd != NULL;
247			sd = seleMan.nextSelected(selei)) {
248	cli2	1407	stuntDoubleIndex.push_back(sd->getGlobalIntegrableObjectIndex());
249
250	cli2	1360	ObjectRestraint* rest = new ObjectRestraint();
251
252			if (stamp[i]->haveDisplacementSpringConstant()) {
253			rest->setDisplacementForceConstant(stamp[i]->getDisplacementSpringConstant());
254			}
255			if (stamp[i]->haveTwistSpringConstant()) {
256			rest->setTwistForceConstant(stamp[i]->getTwistSpringConstant());
257			}
258			if (stamp[i]->haveSwingXSpringConstant()) {
259			rest->setSwingXForceConstant(stamp[i]->getSwingXSpringConstant());
260			}
261			if (stamp[i]->haveSwingYSpringConstant()) {
262			rest->setSwingYForceConstant(stamp[i]->getSwingYSpringConstant());
263			}
264			if (stamp[i]->haveRestrainedTwistAngle()) {
265			rest->setRestrainedTwistAngle(stamp[i]->getRestrainedTwistAngle());
266			}
267			if (stamp[i]->haveRestrainedSwingXAngle()) {
268			rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle());
269			}
270			if (stamp[i]->haveRestrainedSwingYAngle()) {
271			rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle());
272	cli2	1364	}
273			if (stamp[i]->havePrint()) {
274			rest->setPrintRestraint(stamp[i]->getPrint());
275			}
276
277	cli2	1360	restraints_.push_back(rest);
278			sd->addProperty(new RestraintData("Restraint", rest));
279			restrainedObjs_.push_back(sd);
280			}
281
282			}
283			}
284
285			// ThermodynamicIntegration subclasses RestraintForceManager, and there
286			// are times when it won't use restraints at all, so only open the
287			// restraint file if we are actually using restraints:
288
289	cli2	1407	if (simParam->getUseRestraints()) {
290	cli2	1360	std::string refFile = simParam->getRestraint_file();
291	cli2	1407	RestReader* rr = new RestReader(info, refFile, stuntDoubleIndex);
292	cli2	1360	rr->readReferenceStructure();
293			}
294
295			restOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".rest";
296			restOut = new RestWriter(info_, restOutput_.c_str(), restraints_);
297			if(!restOut){
298			sprintf(painCave.errMsg, "Restraint error: Failed to create RestWriter\n");
299			painCave.isFatal = 1;
300			simError();
301			}
302
303			// todo: figure out the scale factor. Right now, just scale them all to 1
304			std::vector<Restraint*>::const_iterator resti;
305			for(resti=restraints_.begin(); resti != restraints_.end(); ++resti){
306			(*resti)->setScaleFactor(1.0);
307			}
308			}
309
310			RestraintForceManager::~RestraintForceManager(){
311			if (restOut)
312			delete restOut;
313			}
314
315			void RestraintForceManager::init() {
316			currRestTime_ = currSnapshot_->getTime();
317			}
318
319	gezelter	1464	void RestraintForceManager::calcForces(){
320	cli2	1360
321	gezelter	1464	ForceManager::calcForces();
322	gezelter	2020	RealType restPot(0.0);
323	cli2	1360
324	gezelter	2020	restPot = doRestraints(1.0);
325	cli2	1360
326			#ifdef IS_MPI
327	gezelter	2020	MPI_Allreduce(MPI_IN_PLACE, &restPot, 1, MPI_REALTYPE, MPI_SUM,
328			MPI_COMM_WORLD);
329	cli2	1360	#endif
330	gezelter	2020
331	cli2	1360	currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
332	gezelter	1782	currSnapshot_->setRestraintPotential(restPot);
333	cli2	1360
334			//write out forces and current positions of restrained molecules
335			if (currSnapshot_->getTime() >= currRestTime_){
336			restOut->writeRest(restInfo_);
337			currRestTime_ += restTime_;
338			}
339			}
340
341			RealType RestraintForceManager::doRestraints(RealType scalingFactor){
342			std::vector<Molecule*>::const_iterator rm;
343			GenericData* data;
344			Molecule::IntegrableObjectIterator ioi;
345			MolecularRestraint* mRest;
346			StuntDouble* sd;
347
348			std::vector<StuntDouble*>::const_iterator ro;
349			ObjectRestraint* oRest;
350
351			std::map<int, Restraint::RealPair> restInfo;
352
353			unscaledPotential_ = 0.0;
354
355			restInfo_.clear();
356
357			for(rm=restrainedMols_.begin(); rm != restrainedMols_.end(); ++rm){
358
359			// make sure this molecule (*rm) has a generic data for restraints:
360			data = (*rm)->getPropertyByName("Restraint");
361			if (data != NULL) {
362			// make sure we can reinterpret the generic data as restraint data:
363			RestraintData* restData= dynamic_cast<RestraintData*>(data);
364			if (restData != NULL) {
365			// make sure we can reinterpet the restraint data as a pointer to
366			// an MolecularRestraint:
367			mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
368			if (mRest == NULL) {
369			sprintf( painCave.errMsg,
370			"Can not cast RestraintData to MolecularRestraint\n");
371	gezelter	1390	painCave.severity = OPENMD_ERROR;
372	cli2	1360	painCave.isFatal = 1;
373			simError();
374			}
375			} else {
376			sprintf( painCave.errMsg,
377			"Can not cast GenericData to RestraintData\n");
378	gezelter	1390	painCave.severity = OPENMD_ERROR;
379	cli2	1360	painCave.isFatal = 1;
380			simError();
381			}
382			} else {
383			sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
384	gezelter	1390	painCave.severity = OPENMD_ERROR;
385	cli2	1360	painCave.isFatal = 1;
386			simError();
387			}
388
389			// phew. At this point, we should have the pointer to the
390			// correct MolecularRestraint in the variable mRest.
391
392			Vector3d molCom = (*rm)->getCom();
393
394			std::vector<Vector3d> struc;
395			std::vector<Vector3d> forces;
396	gezelter	1973
397	cli2	1360
398			for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
399	cli2	1407	sd = (*rm)->nextIntegrableObject(ioi)) {
400	cli2	1360	struc.push_back(sd->getPos());
401			}
402
403			mRest->setScaleFactor(scalingFactor);
404			mRest->calcForce(struc, molCom);
405			forces = mRest->getRestraintForces();
406			int index = 0;
407
408			for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
409			sd = (*rm)->nextIntegrableObject(ioi)) {
410			sd->addFrc(forces[index]);
411			struc.push_back(sd->getPos());
412			index++;
413			}
414
415	gezelter	2020	unscaledPotential_ += mRest->getUnscaledPotential();
416	cli2	1360
417	cli2	1364	// only collect data on restraints that we're going to print:
418			if (mRest->getPrintRestraint())
419	gezelter	2020	restInfo = mRest->getRestraintInfo();
420	cli2	1364	restInfo_.push_back(restInfo);
421	cli2	1360	}
422
423			for(ro=restrainedObjs_.begin(); ro != restrainedObjs_.end(); ++ro){
424			// make sure this object (*ro) has a generic data for restraints:
425			data = (*ro)->getPropertyByName("Restraint");
426			if (data != NULL) {
427			// make sure we can reinterpret the generic data as restraint data:
428			RestraintData* restData= dynamic_cast<RestraintData*>(data);
429			if (restData != NULL) {
430			// make sure we can reinterpet the restraint data as a pointer to
431			// an ObjectRestraint:
432			oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
433			if (oRest == NULL) {
434			sprintf( painCave.errMsg,
435			"Can not cast RestraintData to ObjectRestraint\n");
436	gezelter	1390	painCave.severity = OPENMD_ERROR;
437	cli2	1360	painCave.isFatal = 1;
438			simError();
439			}
440			} else {
441			sprintf( painCave.errMsg,
442			"Can not cast GenericData to RestraintData\n");
443	gezelter	1390	painCave.severity = OPENMD_ERROR;
444	cli2	1360	painCave.isFatal = 1;
445			simError();
446			}
447			} else {
448			sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
449	gezelter	1390	painCave.severity = OPENMD_ERROR;
450	cli2	1360	painCave.isFatal = 1;
451			simError();
452			}
453
454			// phew. At this point, we should have the pointer to the
455			// correct Object restraint in the variable oRest.
456			oRest->setScaleFactor(scalingFactor);
457
458			Vector3d pos = (*ro)->getPos();
459
460			if ( (*ro)->isDirectional() ) {
461
462			// directional objects may have orientational restraints as well
463			// as positional, so get the rotation matrix first:
464
465			RotMat3x3d A = (*ro)->getA();
466			oRest->calcForce(pos, A);
467			(*ro)->addFrc(oRest->getRestraintForce());
468			(*ro)->addTrq(oRest->getRestraintTorque());
469			} else {
470
471			// plain vanilla positional restraints:
472
473			oRest->calcForce(pos);
474			(*ro)->addFrc(oRest->getRestraintForce());
475			}
476
477			unscaledPotential_ += oRest->getUnscaledPotential();
478
479	cli2	1364	// only collect data on restraints that we're going to print:
480			if (oRest->getPrintRestraint())
481	gezelter	2020	restInfo = oRest->getRestraintInfo();
482	cli2	1364	restInfo_.push_back(restInfo);
483	cli2	1360	}
484
485			return unscaledPotential_ * scalingFactor;
486			}
487			}