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

Diff Legend

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