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 1907 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 1908 by gezelter, Fri Jul 19 21:25:45 2013 UTC

# Line 50 | Line 50 | namespace OpenMD {
50   namespace OpenMD {
51  
52    FluctuatingChargeConstraints::FluctuatingChargeConstraints(SimInfo* info) :
53 <    info_(info), hasFlucQ_(false) {
53 >    info_(info), constrainRegions_(false), hasFlucQ_(false) {
54      
55      if (info_->usesFluctuatingCharges()) {
56        if (info_->getNFluctuatingCharges() > 0) {
# Line 59 | Line 59 | namespace OpenMD {
59      }
60    }
61  
62 +  void FluctuatingChargeConstraints::setConstrainRegions(bool cr) {
63 +    constrainRegions_ = cr;
64 +    regionKeys_.clear();
65 +    regionForce_.clear();
66 +    regionCharges_.clear();
67 +
68 +    if (constrainRegions_) {
69 +      SimInfo::MoleculeIterator i;
70 +      Molecule* mol;
71 +      int reg;
72 +      std::set<int> regions;      
73 +
74 +      for (mol = info_->beginMolecule(i); mol != NULL;
75 +           mol = info_->nextMolecule(i)) {
76 +        reg = mol->getRegion();
77 +        if (reg >= 0) regions.insert(reg);
78 +        
79 +      }
80 +      for (std::set<int>::iterator r=regions.begin(); r!=regions.end(); ++r) {
81 +        regionKeys_.push_back( (*r) );
82 +      }
83 +      regionForce_.resize( regionKeys_.size() );
84 +      regionCharges_.resize( regionKeys_.size() );
85 +    }
86 +  }
87 +
88 +
89    void FluctuatingChargeConstraints::applyConstraints() {
90      if (!hasFlucQ_) return;
91      
# Line 67 | Line 94 | namespace OpenMD {
94      Molecule* mol;
95      Atom* atom;
96      
97 <    RealType totalFrc, totalMolFrc, constrainedFrc;
97 >    RealType totalFrc, totalMolFrc, regionFrc, constrainedFrc;
98      
99      // accumulate the total system fluctuating charge forces
100      totalFrc = 0.0;
101 +    if (constrainRegions_) {
102 +      std::fill(regionForce_.begin(), regionForce_.end(), 0.0);
103 +      std::fill(regionCharges_.begin(), regionCharges_.end(), 0);
104 +    }
105  
106      for (mol = info_->beginMolecule(i); mol != NULL;
107           mol = info_->nextMolecule(i)) {
108  
109 +      int region = mol->getRegion();
110 +          
111        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
112             atom = mol->nextFluctuatingCharge(j)) {
80 –        totalFrc += atom->getFlucQFrc();
81 –      }
113  
114 +        RealType frc = atom->getFlucQFrc();
115 +        totalFrc += frc;
116 +        if (constrainRegions_) {
117 +          regionForce_[regionKeys_[region]] += frc;
118 +          regionCharges_[regionKeys_[region]] += 1;
119 +        }
120 +      }
121      }
122  
123   #ifdef IS_MPI
# Line 87 | Line 125 | namespace OpenMD {
125      // processors:
126      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalFrc, 1, MPI::REALTYPE,
127                                MPI::SUM);
128 +
129 +    if (constrainRegions_) {
130 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &regionForce_[0],
131 +                                regionForce_.size(), MPI::REALTYPE, MPI::SUM);
132 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &regionCharges_[0],
133 +                                regionCharges_.size(), MPI::INT, MPI::SUM);
134 +    }
135 +
136   #endif
137 <
137 >
138      // divide by the total number of fluctuating charges:
139      totalFrc /= info_->getNFluctuatingCharges();
140 +    
141 +    // do the same in the regions:
142 +    if (constrainRegions_) {
143 +      for (int i = 0; i < regionForce_.size(); ++i)  {
144 +        regionForce_[ i ] /= regionCharges_[ i ];
145 +      }
146 +    }
147  
148      for (mol = info_->beginMolecule(i); mol != NULL;
149           mol = info_->nextMolecule(i)) {    
150        
151 <      totalMolFrc = 0.0;
151 >      if (constrainRegions_) {
152 >        int region = mol->getRegion();
153 >        regionFrc = regionForce_[regionKeys_[region]];
154 >      } else {
155 >        regionFrc = 0.0;
156 >      }
157  
158 +      totalMolFrc = 0.0;
159 +      
160        // molecular constraints can be done with a second loop.
161        if (mol->constrainTotalCharge()) {
162          for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
# Line 109 | Line 169 | namespace OpenMD {
169        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
170             atom = mol->nextFluctuatingCharge(j)) {
171          //constrainedFrc = atom->getFlucQFrc() - totalFrc - totalMolFrc;
172 +
173          constrainedFrc = atom->getFlucQFrc() - totalMolFrc;
174 +
175 +        if (constrainRegions_)
176 +          constrainedFrc -= regionFrc;
177 +        
178          atom->setFlucQFrc(constrainedFrc);
179        }      
180      }

Diff Legend

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