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
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
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Acknowledgement of the program authors must be made in any
10	* publication of scientific results based in part on use of the
11	* program. An acceptable form of acknowledgement is citation of
12	* the article in which the program was described (Matthew
13	* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	* Parallel Simulation Engine for Molecular Dynamics,"
16	* J. Comput. Chem. 26, pp. 252-271 (2005))
17	*
18	* 2. Redistributions of source code must retain the above copyright
19	* notice, this list of conditions and the following disclaimer.
20	*
21	* 3. Redistributions in binary form must reproduce the above copyright
22	* notice, this list of conditions and the following disclaimer in the
23	* documentation and/or other materials provided with the
24	* distribution.
25	*
26	* This software is provided "AS IS," without a warranty of any
27	* kind. All express or implied conditions, representations and
28	* warranties, including any implied warranty of merchantability,
29	* fitness for a particular purpose or non-infringement, are hereby
30	* excluded. The University of Notre Dame and its licensors shall not
31	* be liable for any damages suffered by licensee as a result of
32	* using, modifying or distributing the software or its
33	* derivatives. In no event will the University of Notre Dame or its
34	* licensors be liable for any lost revenue, profit or data, or for
35	* direct, indirect, special, consequential, incidental or punitive
36	* damages, however caused and regardless of the theory of liability,
37	* arising out of the use of or inability to use software, even if the
38	* University of Notre Dame has been advised of the possibility of
39	* such damages.
40	*/
41
42
43	#include "restraints/ThermoIntegrationForceManager.hpp"
44
45	#ifdef IS_MPI
46	#include <mpi.h>
47	#endif
48
49	namespace oopse {
50
51	ThermoIntegrationForceManager::ThermoIntegrationForceManager(SimInfo* info):
52	RestraintForceManager(info){
53	currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
54	simParam = info_->getSimParams();
55
56	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
71	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
86	// build the scaling factor used to modulate the forces and torques
87	factor_ = pow(tIntLambda_, tIntK_);
88	}
89
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	Mat3x3d tempTau;
103
104	// perform the standard calcForces first
105	ForceManager::calcForces(needPotential, needStress);
106
107	curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
108
109	// 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
128	// set vraw to be the unmodulated potential
129	lrPot_ = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
130	curSnapshot->statData[Stats::VRAW] = lrPot_;
131
132	// modulate the potential and update the snapshot
133	lrPot_ *= factor_;
134	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
135
136	// scale the pressure tensor
137	tempTau = curSnapshot->statData.getTau();
138	tempTau *= factor_;
139	curSnapshot->statData.setTau(tempTau);
140
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
152	#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	#else
160	lrPot_ += restPot_local;
161	vHarm_ = vHarm_local;
162	#endif
163
164	// give the final values to stats
165	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
166	curSnapshot->statData[Stats::VHARM] = vHarm_;
167	}
168	}