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 1879 by gezelter, Sun Jun 16 15:15:42 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), 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();
75 +
76 +    if (constrainRegions_) {
77 +      SimInfo::MoleculeIterator i;
78 +      Molecule* mol;
79 +      int reg;
80 +      std::set<int> regions;      
81 +
82 +      for (mol = info_->beginMolecule(i); mol != NULL;
83 +           mol = info_->nextMolecule(i)) {
84 +        reg = mol->getRegion();
85 +        if (reg >= 0) regions.insert(reg);
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_[ (*r) ] = which;
92 +        which++;
93 +      }
94 +      regionForce_.resize( regionKeys_.size() );
95 +      regionCharges_.resize( regionKeys_.size() );
96 +    }
97 +  }
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, constrainedFrc;
71 <    
109 >    RealType totalFrc, totalMolFrc, regionFrc, constrainedFrc;
110      // accumulate the total system fluctuating charge forces
111      totalFrc = 0.0;
112 +    if (constrainRegions_) {
113 +      std::fill(regionForce_.begin(), regionForce_.end(), 0.0);
114 +      std::fill(regionCharges_.begin(), regionCharges_.end(), 0);
115 +    }
116  
117      for (mol = info_->beginMolecule(i); mol != NULL;
118           mol = info_->nextMolecule(i)) {
119  
120 +      int region = mol->getRegion();
121 +          
122        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
123             atom = mol->nextFluctuatingCharge(j)) {
80 –        totalFrc += atom->getFlucQFrc();
81 –      }
124  
125 +        RealType frc = atom->getFlucQFrc();
126 +        totalFrc += frc;
127 +        if (constrainRegions_) {
128 +          if (region >= 0) {        
129 +            regionForce_[regionKeys_[region]] += frc;
130 +            regionCharges_[regionKeys_[region]] += 1;
131 +          }
132 +        }
133 +      }
134      }
135  
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_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
150 <
150 >
151      // divide by the total number of fluctuating charges:
152      totalFrc /= info_->getNFluctuatingCharges();
153 +    
154 +    // do the same in the regions:
155 +    if (constrainRegions_) {
156 +      for (unsigned int i = 0; i < regionForce_.size(); ++i)  {
157 +        regionForce_[ i ] /= regionCharges_[ i ];
158 +      }
159 +    }
160  
161      for (mol = info_->beginMolecule(i); mol != NULL;
162           mol = info_->nextMolecule(i)) {    
163        
164 <      totalMolFrc = 0.0;
164 >      if (constrainRegions_) {
165 >        int region = mol->getRegion();
166 >        if (region >= 0)
167 >          regionFrc = regionForce_[regionKeys_[region]];
168 >        else
169 >          regionFrc = 0.0;
170 >      } else {
171 >        regionFrc = 0.0;
172 >      }
173  
174 +      totalMolFrc = 0.0;
175 +      
176        // molecular constraints can be done with a second loop.
177        if (mol->constrainTotalCharge()) {
178          for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
# Line 108 | Line 184 | namespace OpenMD {
184  
185        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
186             atom = mol->nextFluctuatingCharge(j)) {
187 <        //constrainedFrc = atom->getFlucQFrc() - totalFrc - totalMolFrc;
188 <        constrainedFrc = atom->getFlucQFrc() - totalMolFrc;
187 >        constrainedFrc = atom->getFlucQFrc() - totalFrc - totalMolFrc;
188 >
189 >        //constrainedFrc = atom->getFlucQFrc() - totalMolFrc;
190 >
191 >        if (constrainRegions_)
192 >          constrainedFrc -= regionFrc;
193 >        
194          atom->setFlucQFrc(constrainedFrc);
195        }      
196      }

Diff Legend

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