ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/restraints/ThermoIntegrationForceManager.cpp
(Generate patch)

Comparing trunk/src/restraints/ThermoIntegrationForceManager.cpp (file contents):
Revision 419 by chrisfen, Thu Mar 10 16:15:13 2005 UTC vs.
Revision 1028 by chrisfen, Thu Aug 31 20:12:12 2006 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 42 | Line 42
42   #include <cmath>
43   #include "restraints/ThermoIntegrationForceManager.hpp"
44   #include "integrators/Integrator.hpp"
45 + #include "math/SquareMatrix3.hpp"
46   #include "primitives/Molecule.hpp"
47   #include "utils/simError.h"
48   #include "utils/OOPSEConstant.hpp"
49   #include "utils/StringUtils.hpp"
50  
51 + #ifdef IS_MPI
52 + #include <mpi.h>
53 + #define TAKE_THIS_TAG_REAL 2
54 + #endif //is_mpi
55 +
56   namespace oopse {
57    
58    ThermoIntegrationForceManager::ThermoIntegrationForceManager(SimInfo* info):
59 <  ForceManager(info){
60 <    currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
61 <    simParam = info_->getSimParams();
59 >    ForceManager(info){
60 >      currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
61 >      simParam = info_->getSimParams();
62      
63 <    if (simParam->haveThermIntLambda()){
64 <      tIntLambda_ = simParam->getThermIntLambda();
65 <    }
66 <    else{
67 <      tIntLambda_ = 1.0;
68 <      sprintf(painCave.errMsg,
69 <              "ThermoIntegration error: the transformation parameter (lambda) was\n"
70 <              "\tnot specified. OOPSE will use a default value of %f. To set\n"
71 <              "\tlambda, use the thermodynamicIntegrationLambda variable.\n",
72 <              tIntLambda_);
73 <      painCave.isFatal = 0;
74 <      simError();
75 <    }
63 >      if (simParam->haveThermodynamicIntegrationLambda()){
64 >        tIntLambda_ = simParam->getThermodynamicIntegrationLambda();
65 >      }
66 >      else{
67 >        tIntLambda_ = 1.0;
68 >        sprintf(painCave.errMsg,
69 >                "ThermoIntegration error: the transformation parameter\n"
70 >                "\t(lambda) was not specified. OOPSE will use a default\n"
71 >                "\tvalue of %f. To set lambda, use the \n"
72 >                "\tthermodynamicIntegrationLambda variable.\n",
73 >                tIntLambda_);
74 >        painCave.isFatal = 0;
75 >        simError();
76 >      }
77      
78 <    if (simParam->haveThermIntK()){
79 <      tIntK_ = simParam->getThermIntK();
80 <    }
81 <    else{
82 <      tIntK_ = 1.0;
83 <      sprintf(painCave.errMsg,
84 <              "ThermoIntegration Warning: the tranformation parameter exponent\n"
85 <              "\t(k) was not specified. OOPSE will use a default value of %f.\n"
86 <              "\tTo set k, use the thermodynamicIntegrationK variable.\n",
87 <              tIntK_);
88 <      painCave.isFatal = 0;
89 <      simError();      
90 <    }
78 >      if (simParam->haveThermodynamicIntegrationK()){
79 >        tIntK_ = simParam->getThermodynamicIntegrationK();
80 >      }
81 >      else{
82 >        tIntK_ = 1.0;
83 >        sprintf(painCave.errMsg,
84 >                "ThermoIntegration Warning: the tranformation parameter\n"
85 >                "\texponent (k) was not specified. OOPSE will use a default\n"
86 >                "\tvalue of %f. To set k, use the thermodynamicIntegrationK\n"
87 >                "\tvariable.\n",
88 >                tIntK_);
89 >        painCave.isFatal = 0;
90 >        simError();      
91 >      }
92      
93 <    if (simParam->getUseSolidThermInt()) {
94 <      // build a restraint object
95 <      restraint_ =  new Restraints(info_, tIntLambda_, tIntK_);
93 >      if (simParam->getUseSolidThermInt()) {
94 >        // build a restraint object
95 >        restraint_ =  new Restraints(info_, tIntLambda_, tIntK_);
96        
97 <    }
97 >      }
98      
99 <    // build the scaling factor used to modulate the forces and torques
100 <    factor_ = pow(tIntLambda_, tIntK_);
99 >      // build the scaling factor used to modulate the forces and torques
100 >      factor_ = pow(tIntLambda_, tIntK_);
101  
102 <    printf("%f is the factor\n",factor_);
95 <    
96 <  }
102 >    }
103    
104    ThermoIntegrationForceManager::~ThermoIntegrationForceManager(){
105    }
# Line 107 | Line 113 | namespace oopse {
113      StuntDouble* integrableObject;
114      Vector3d frc;
115      Vector3d trq;
116 +    Mat3x3d tempTau;
117      
118      // perform the standard calcForces first
119      ForceManager::calcForces(needPotential, needStress);
# Line 130 | Line 137 | namespace oopse {
137            integrableObject->setTrq(trq);
138          }
139        }
133 –      
134 –      // set vraw to be the unmodulated potential
135 –      lrPot_ = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
136 –      curSnapshot->statData[Stats::VRAW] = lrPot_;
137 –      
138 –      // modulate the potential and update the snapshot
139 –      lrPot_ *= factor_;
140 –      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
140      }
141 +  
142 +    // set vraw to be the unmodulated potential
143 +    lrPot_ = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
144 +    curSnapshot->statData[Stats::VRAW] = lrPot_;
145      
146 <    // do crystal restraint forces for thermodynamic integration
146 >    // modulate the potential and update the snapshot
147 >    lrPot_ *= factor_;
148 >    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
149 >    
150 >    // scale the pressure tensor
151 >    tempTau = curSnapshot->statData.getTau();
152 >    tempTau *= factor_;
153 >    curSnapshot->statData.setTau(tempTau);
154 > #ifndef IS_MPI
155 >    // do the single processor crystal restraint forces for
156 >    // thermodynamic integration
157      if (simParam->getUseSolidThermInt()) {
158        
159        lrPot_ += restraint_->Calc_Restraint_Forces();
# Line 149 | Line 162 | namespace oopse {
162        vHarm_ = restraint_->getVharm();
163        curSnapshot->statData[Stats::VHARM] = vHarm_;
164      }
165 <    
165 > #else
166 >    double tempLRPot = 0.0;
167 >    double tempVHarm = 0.0;
168 >    MPI_Status ierr;
169 >    int nproc;
170 >    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
171 >
172 >    // do the MPI crystal restraint forces for each processor
173 >    if (simParam->getUseSolidThermInt()) {
174 >      tempLRPot = restraint_->Calc_Restraint_Forces();
175 >      tempVHarm = restraint_->getVharm();
176 >    }
177 >
178 >    // master receives and accumulates the restraint info
179 >    if (worldRank == 0) {
180 >      for(int i = 0 ; i < nproc; ++i) {
181 >        if (i == worldRank) {
182 >          lrPot_ += tempLRPot;
183 >          vHarm_ += tempVHarm;
184 >        } else {
185 >          MPI_Recv(&tempLRPot, 1, MPI_REALTYPE, i,
186 >                   TAKE_THIS_TAG_REAL, MPI_COMM_WORLD, &ierr);
187 >          MPI_Recv(&tempVHarm, 1, MPI_REALTYPE, i,
188 >                   TAKE_THIS_TAG_REAL, MPI_COMM_WORLD, &ierr);
189 >          lrPot_ += tempLRPot;
190 >          vHarm_ += tempVHarm;
191 >        }
192 >      }
193 >
194 >      // give the final values to stats
195 >      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
196 >      curSnapshot->statData[Stats::VHARM] = vHarm_;
197 >
198 >    } else {
199 >      // pack up and send the appropriate info to the master node
200 >      for(int j = 1; j < nproc; ++j) {
201 >        if (worldRank == j) {
202 >
203 >          MPI_Send(&tempLRPot, 1, MPI_REALTYPE, 0,
204 >                   TAKE_THIS_TAG_REAL, MPI_COMM_WORLD);
205 >          MPI_Send(&tempVHarm, 1, MPI_REALTYPE, 0,
206 >                   TAKE_THIS_TAG_REAL, MPI_COMM_WORLD);
207 >        }
208 >      }
209 >    }
210 > #endif //is_mpi
211    }
212    
213   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines