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

Comparing trunk/src/flucq/FluctuatingChargeConstraints.cpp (file contents):
Revision 1908 by gezelter, Fri Jul 19 21:25:45 2013 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 UTC

# Line 40 | Line 40
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43 – #include "FluctuatingChargeConstraints.hpp"
44 – #include "primitives/Molecule.hpp"
45 –
43   #ifdef IS_MPI
44   #include <mpi.h>
45   #endif
46  
47 + #include "FluctuatingChargeConstraints.hpp"
48 + #include "primitives/Molecule.hpp"
49 +
50   namespace OpenMD {
51  
52    FluctuatingChargeConstraints::FluctuatingChargeConstraints(SimInfo* info) :
53 <    info_(info), constrainRegions_(false), hasFlucQ_(false) {
54 <    
55 <    if (info_->usesFluctuatingCharges()) {
56 <      if (info_->getNFluctuatingCharges() > 0) {
57 <        hasFlucQ_ = true;
53 >    info_(info), constrainRegions_(false), hasFlucQ_(false), initialized_(false) {
54 >  }
55 >
56 >  void FluctuatingChargeConstraints::initialize(){
57 >    if(info_->usesFluctuatingCharges()){
58 >      if(info_->getNFluctuatingCharges() > 0){
59 >        hasFlucQ_ = true;
60        }
61      }
62 +    initialized_ = true;
63    }
64  
65 +
66    void FluctuatingChargeConstraints::setConstrainRegions(bool cr) {
67      constrainRegions_ = cr;
68 +
69 +    if (!initialized_) initialize();
70 +
71      regionKeys_.clear();
72      regionForce_.clear();
73      regionCharges_.clear();
# Line 75 | Line 82 | namespace OpenMD {
82             mol = info_->nextMolecule(i)) {
83          reg = mol->getRegion();
84          if (reg >= 0) regions.insert(reg);
78 –        
85        }
86 +      // resize the keys vector to the largest found value for regions.
87 +      regionKeys_.resize( *(regions.end()) );
88 +      int which = 0;
89        for (std::set<int>::iterator r=regions.begin(); r!=regions.end(); ++r) {
90 <        regionKeys_.push_back( (*r) );
90 >        regionKeys_[ (*r) ] = which;
91 >        which++;
92        }
93        regionForce_.resize( regionKeys_.size() );
94        regionCharges_.resize( regionKeys_.size() );
# Line 87 | Line 97 | namespace OpenMD {
97  
98  
99    void FluctuatingChargeConstraints::applyConstraints() {
100 +    if (!initialized_) initialize();
101      if (!hasFlucQ_) return;
102 <    
102 >
103      SimInfo::MoleculeIterator i;
104      Molecule::FluctuatingChargeIterator  j;
105      Molecule* mol;
106      Atom* atom;
107      
108      RealType totalFrc, totalMolFrc, regionFrc, constrainedFrc;
98 –    
109      // accumulate the total system fluctuating charge forces
110      totalFrc = 0.0;
111      if (constrainRegions_) {
# Line 114 | Line 124 | namespace OpenMD {
124          RealType frc = atom->getFlucQFrc();
125          totalFrc += frc;
126          if (constrainRegions_) {
127 <          regionForce_[regionKeys_[region]] += frc;
128 <          regionCharges_[regionKeys_[region]] += 1;
127 >          if (region >= 0) {        
128 >            regionForce_[regionKeys_[region]] += frc;
129 >            regionCharges_[regionKeys_[region]] += 1;
130 >          }
131          }
132        }
133      }
# Line 123 | Line 135 | namespace OpenMD {
135   #ifdef IS_MPI
136      // in parallel, we need to add up the contributions from all
137      // processors:
138 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalFrc, 1, MPI::REALTYPE,
139 <                              MPI::SUM);
138 >    MPI_Allreduce(MPI_IN_PLACE, &totalFrc, 1, MPI_REALTYPE,
139 >                  MPI_SUM, MPI_COMM_WORLD);
140  
141      if (constrainRegions_) {
142 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &regionForce_[0],
143 <                                regionForce_.size(), MPI::REALTYPE, MPI::SUM);
144 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &regionCharges_[0],
145 <                                regionCharges_.size(), MPI::INT, MPI::SUM);
142 >      MPI_Allreduce(MPI_IN_PLACE, &regionForce_[0],
143 >                    regionForce_.size(), MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
144 >      MPI_Allreduce(MPI_IN_PLACE, &regionCharges_[0],
145 >                    regionCharges_.size(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
146      }
147  
148   #endif
# Line 140 | Line 152 | namespace OpenMD {
152      
153      // do the same in the regions:
154      if (constrainRegions_) {
155 <      for (int i = 0; i < regionForce_.size(); ++i)  {
155 >      for (unsigned int i = 0; i < regionForce_.size(); ++i)  {
156          regionForce_[ i ] /= regionCharges_[ i ];
157        }
158      }
# Line 150 | Line 162 | namespace OpenMD {
162        
163        if (constrainRegions_) {
164          int region = mol->getRegion();
165 <        regionFrc = regionForce_[regionKeys_[region]];
165 >        if (region >= 0)
166 >          regionFrc = regionForce_[regionKeys_[region]];
167 >        else
168 >          regionFrc = 0.0;
169        } else {
170          regionFrc = 0.0;
171        }
# Line 168 | Line 183 | namespace OpenMD {
183  
184        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
185             atom = mol->nextFluctuatingCharge(j)) {
186 <        //constrainedFrc = atom->getFlucQFrc() - totalFrc - totalMolFrc;
186 >        constrainedFrc = atom->getFlucQFrc() - totalFrc - totalMolFrc;
187  
188 <        constrainedFrc = atom->getFlucQFrc() - totalMolFrc;
188 >        //constrainedFrc = atom->getFlucQFrc() - totalMolFrc;
189  
190          if (constrainRegions_)
191            constrainedFrc -= regionFrc;

Diff Legend

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