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 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC vs.
Revision 1939 by gezelter, Thu Oct 31 18:18:57 2013 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
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), 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();
74 +
75 +    if (constrainRegions_) {
76 +      SimInfo::MoleculeIterator i;
77 +      Molecule* mol;
78 +      int reg;
79 +      std::set<int> regions;      
80 +
81 +      for (mol = info_->beginMolecule(i); mol != NULL;
82 +           mol = info_->nextMolecule(i)) {
83 +        reg = mol->getRegion();
84 +        if (reg >= 0) regions.insert(reg);
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_[ (*r) ] = which;
91 +        which++;
92 +      }
93 +      regionForce_.resize( regionKeys_.size() );
94 +      regionCharges_.resize( regionKeys_.size() );
95 +    }
96 +  }
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, constrainedFrc;
71 <    
108 >    RealType totalFrc, totalMolFrc, regionFrc, constrainedFrc;
109      // accumulate the total system fluctuating charge forces
110      totalFrc = 0.0;
111 +    if (constrainRegions_) {
112 +      std::fill(regionForce_.begin(), regionForce_.end(), 0.0);
113 +      std::fill(regionCharges_.begin(), regionCharges_.end(), 0);
114 +    }
115  
116      for (mol = info_->beginMolecule(i); mol != NULL;
117           mol = info_->nextMolecule(i)) {
118  
119 +      int region = mol->getRegion();
120 +          
121        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
122             atom = mol->nextFluctuatingCharge(j)) {
80 –        totalFrc += atom->getFlucQFrc();
81 –      }
123  
124 +        RealType frc = atom->getFlucQFrc();
125 +        totalFrc += frc;
126 +        if (constrainRegions_) {
127 +          if (region >= 0) {        
128 +            regionForce_[regionKeys_[region]] += frc;
129 +            regionCharges_[regionKeys_[region]] += 1;
130 +          }
131 +        }
132 +      }
133      }
134  
135   #ifdef IS_MPI
# Line 87 | Line 137 | namespace OpenMD {
137      // processors:
138      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalFrc, 1, MPI::REALTYPE,
139                                MPI::SUM);
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);
146 +    }
147 +
148   #endif
149 <
149 >
150      // divide by the total number of fluctuating charges:
151      totalFrc /= info_->getNFluctuatingCharges();
152 +    
153 +    // do the same in the regions:
154 +    if (constrainRegions_) {
155 +      for (unsigned int i = 0; i < regionForce_.size(); ++i)  {
156 +        regionForce_[ i ] /= regionCharges_[ i ];
157 +      }
158 +    }
159  
160      for (mol = info_->beginMolecule(i); mol != NULL;
161           mol = info_->nextMolecule(i)) {    
162        
163 <      totalMolFrc = 0.0;
163 >      if (constrainRegions_) {
164 >        int region = mol->getRegion();
165 >        if (region >= 0)
166 >          regionFrc = regionForce_[regionKeys_[region]];
167 >        else
168 >          regionFrc = 0.0;
169 >      } else {
170 >        regionFrc = 0.0;
171 >      }
172  
173 +      totalMolFrc = 0.0;
174 +      
175        // molecular constraints can be done with a second loop.
176        if (mol->constrainTotalCharge()) {
177          for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
# Line 108 | 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;
189 >
190 >        if (constrainRegions_)
191 >          constrainedFrc -= regionFrc;
192 >        
193          atom->setFlucQFrc(constrainedFrc);
194        }      
195      }

Diff Legend

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