src/restraints/ThermoIntegrationForceManager.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 

#include "restraints/ThermoIntegrationForceManager.hpp"

#ifdef IS_MPI
#include <mpi.h>
#endif

namespace oopse {
  
  ThermoIntegrationForceManager::ThermoIntegrationForceManager(SimInfo* info): 
    RestraintForceManager(info){
    currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
    simParam = info_->getSimParams();
    
    if (simParam->haveThermodynamicIntegrationLambda()){
      tIntLambda_ = simParam->getThermodynamicIntegrationLambda();
    }
    else{
      tIntLambda_ = 1.0;
      sprintf(painCave.errMsg,
              "ThermoIntegration error: the transformation parameter\n"
              "\t(lambda) was not specified. OOPSE will use a default\n"
              "\tvalue of %f. To set lambda, use the \n"
              "\tthermodynamicIntegrationLambda variable.\n",
              tIntLambda_);
      painCave.isFatal = 0;
      simError();
    }
    
    if (simParam->haveThermodynamicIntegrationK()){
      tIntK_ = simParam->getThermodynamicIntegrationK();
    }
    else{
      tIntK_ = 1.0;
      sprintf(painCave.errMsg,
              "ThermoIntegration Warning: the tranformation parameter\n"
              "\texponent (k) was not specified. OOPSE will use a default\n"
              "\tvalue of %f. To set k, use the thermodynamicIntegrationK\n"
              "\tvariable.\n",
              tIntK_);
      painCave.isFatal = 0;
      simError();      
    }
    
    // build the scaling factor used to modulate the forces and torques
    factor_ = pow(tIntLambda_, tIntK_);
  }
  
  ThermoIntegrationForceManager::~ThermoIntegrationForceManager(){
  }
  
  void ThermoIntegrationForceManager::calcForces(bool needPotential, 
                                                 bool needStress){
    Snapshot* curSnapshot;
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    Molecule::IntegrableObjectIterator ii;
    StuntDouble* integrableObject;
    Vector3d frc;
    Vector3d trq;
    Mat3x3d tempTau;
    
    // perform the standard calcForces first
    ForceManager::calcForces(needPotential, needStress);
    
    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();

    // now scale forces and torques of all the integrableObjects
      
    for (mol = info_->beginMolecule(mi); mol != NULL; 
         mol = info_->nextMolecule(mi)) {
      for (integrableObject = mol->beginIntegrableObject(ii); 
           integrableObject != NULL; 
           integrableObject = mol->nextIntegrableObject(ii)) {
        frc = integrableObject->getFrc();
        frc *= factor_;
        integrableObject->setFrc(frc);
        
        if (integrableObject->isDirectional()){
          trq = integrableObject->getTrq();
          trq *= factor_;
          integrableObject->setTrq(trq);
        }
      }
    }
    
    // set vraw to be the unmodulated potential
    lrPot_ = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
    curSnapshot->statData[Stats::VRAW] = lrPot_;
    
    // modulate the potential and update the snapshot
    lrPot_ *= factor_;
    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
    
    // scale the pressure tensor
    tempTau = curSnapshot->statData.getTau();
    tempTau *= factor_;
    curSnapshot->statData.setTau(tempTau);

    // now, on to the applied restraining potentials (if needed):
    RealType restPot_local = 0.0;
    RealType vHarm_local = 0.0;
    
    if (simParam->getUseRestraints()) {
      // do restraints from RestraintForceManager:
      //restPot_local = doRestraints(1.0 - factor_);     
      restPot_local = doRestraints(1.0 - factor_);      
      vHarm_local = getUnscaledPotential();
    }
      
#ifdef IS_MPI
    RealType restPot;
    MPI::COMM_WORLD.Allreduce(&restPot_local, &restPot, 1, 
                              MPI::REALTYPE, MPI::SUM);
    MPI::COMM_WORLD.Allreduce(&vHarm_local, &vHarm_, 1, 
                              MPI::REALTYPE, MPI::SUM);         
    lrPot_ += restPot;
#else
    lrPot_ += restPot_local;
    vHarm_ = vHarm_local;
#endif

    // give the final values to stats
    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
    curSnapshot->statData[Stats::VHARM] = vHarm_;
  }  
}
Revision:	1360
Committed:	Mon Sep 7 16:31:51 2009 UTC (15 years, 7 months ago) by cli2
Original Path:	*trunk/src/restraints/ThermoIntegrationForceManager.cpp*
File size:	6236 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	gezelter	507	/*
2	chrisfen	417	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
41
42	cli2	1360
43	chrisfen	417	#include "restraints/ThermoIntegrationForceManager.hpp"
44
45	chrisfen	1028	#ifdef IS_MPI
46			#include <mpi.h>
47	cli2	1360	#endif
48	chrisfen	1028
49	chrisfen	417	namespace oopse {
50
51			ThermoIntegrationForceManager::ThermoIntegrationForceManager(SimInfo* info):
52	cli2	1360	RestraintForceManager(info){
53			currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
54			simParam = info_->getSimParams();
55	chrisfen	417
56	cli2	1360	if (simParam->haveThermodynamicIntegrationLambda()){
57			tIntLambda_ = simParam->getThermodynamicIntegrationLambda();
58			}
59			else{
60			tIntLambda_ = 1.0;
61			sprintf(painCave.errMsg,
62			"ThermoIntegration error: the transformation parameter\n"
63			"\t(lambda) was not specified. OOPSE will use a default\n"
64			"\tvalue of %f. To set lambda, use the \n"
65			"\tthermodynamicIntegrationLambda variable.\n",
66			tIntLambda_);
67			painCave.isFatal = 0;
68			simError();
69			}
70	chrisfen	417
71	cli2	1360	if (simParam->haveThermodynamicIntegrationK()){
72			tIntK_ = simParam->getThermodynamicIntegrationK();
73			}
74			else{
75			tIntK_ = 1.0;
76			sprintf(painCave.errMsg,
77			"ThermoIntegration Warning: the tranformation parameter\n"
78			"\texponent (k) was not specified. OOPSE will use a default\n"
79			"\tvalue of %f. To set k, use the thermodynamicIntegrationK\n"
80			"\tvariable.\n",
81			tIntK_);
82			painCave.isFatal = 0;
83			simError();
84			}
85	chrisfen	417
86	cli2	1360	// build the scaling factor used to modulate the forces and torques
87			factor_ = pow(tIntLambda_, tIntK_);
88			}
89	chrisfen	417
90			ThermoIntegrationForceManager::~ThermoIntegrationForceManager(){
91			}
92
93			void ThermoIntegrationForceManager::calcForces(bool needPotential,
94			bool needStress){
95			Snapshot* curSnapshot;
96			SimInfo::MoleculeIterator mi;
97			Molecule* mol;
98			Molecule::IntegrableObjectIterator ii;
99			StuntDouble* integrableObject;
100			Vector3d frc;
101			Vector3d trq;
102	chrisfen	423	Mat3x3d tempTau;
103	chrisfen	417
104			// perform the standard calcForces first
105			ForceManager::calcForces(needPotential, needStress);
106
107	chrisfen	419	curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
108
109	chrisfen	417	// now scale forces and torques of all the integrableObjects
110
111			for (mol = info_->beginMolecule(mi); mol != NULL;
112			mol = info_->nextMolecule(mi)) {
113			for (integrableObject = mol->beginIntegrableObject(ii);
114			integrableObject != NULL;
115			integrableObject = mol->nextIntegrableObject(ii)) {
116			frc = integrableObject->getFrc();
117			frc *= factor_;
118			integrableObject->setFrc(frc);
119
120			if (integrableObject->isDirectional()){
121			trq = integrableObject->getTrq();
122			trq *= factor_;
123			integrableObject->setTrq(trq);
124			}
125			}
126			}
127	cli2	1360
128	chrisfen	420	// set vraw to be the unmodulated potential
129			lrPot_ = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
130			curSnapshot->statData[Stats::VRAW] = lrPot_;
131	chrisfen	417
132	chrisfen	420	// modulate the potential and update the snapshot
133			lrPot_ *= factor_;
134			curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
135
136	chrisfen	423	// scale the pressure tensor
137			tempTau = curSnapshot->statData.getTau();
138			tempTau *= factor_;
139			curSnapshot->statData.setTau(tempTau);
140	cli2	1360
141			// now, on to the applied restraining potentials (if needed):
142			RealType restPot_local = 0.0;
143			RealType vHarm_local = 0.0;
144
145			if (simParam->getUseRestraints()) {
146			// do restraints from RestraintForceManager:
147			//restPot_local = doRestraints(1.0 - factor_);
148			restPot_local = doRestraints(1.0 - factor_);
149			vHarm_local = getUnscaledPotential();
150			}
151	chrisfen	417
152	cli2	1360	#ifdef IS_MPI
153			RealType restPot;
154			MPI::COMM_WORLD.Allreduce(&restPot_local, &restPot, 1,
155			MPI::REALTYPE, MPI::SUM);
156			MPI::COMM_WORLD.Allreduce(&vHarm_local, &vHarm_, 1,
157			MPI::REALTYPE, MPI::SUM);
158			lrPot_ += restPot;
159	chrisfen	1028	#else
160	cli2	1360	lrPot_ += restPot_local;
161			vHarm_ = vHarm_local;
162			#endif
163	chrisfen	1028
164	cli2	1360	// give the final values to stats
165			curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
166			curSnapshot->statData[Stats::VHARM] = vHarm_;
167			}
168	chrisfen	417	}