ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/devel_omp/src/brains/ForceManager.cpp
Revision: 1608
Committed: Tue Aug 9 01:58:56 2011 UTC (13 years, 8 months ago) by mciznick
File size: 42806 byte(s)
Log Message: 
First OpenMP version.

File Contents

# User Rev Content
1 gezelter 507 /*
2 gezelter 246 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     *
4     * The University of Notre Dame grants you ("Licensee") a
5     * non-exclusive, royalty free, license to use, modify and
6     * redistribute this software in source and binary code form, provided
7     * that the following conditions are met:
8     *
9 gezelter 1390 * 1. Redistributions of source code must retain the above copyright
10 gezelter 246 * notice, this list of conditions and the following disclaimer.
11     *
12 gezelter 1390 * 2. Redistributions in binary form must reproduce the above copyright
13 gezelter 246 * notice, this list of conditions and the following disclaimer in the
14     * documentation and/or other materials provided with the
15     * distribution.
16     *
17     * This software is provided "AS IS," without a warranty of any
18     * kind. All express or implied conditions, representations and
19     * warranties, including any implied warranty of merchantability,
20     * fitness for a particular purpose or non-infringement, are hereby
21     * excluded. The University of Notre Dame and its licensors shall not
22     * be liable for any damages suffered by licensee as a result of
23     * using, modifying or distributing the software or its
24     * derivatives. In no event will the University of Notre Dame or its
25     * licensors be liable for any lost revenue, profit or data, or for
26     * direct, indirect, special, consequential, incidental or punitive
27     * damages, however caused and regardless of the theory of liability,
28     * arising out of the use of or inability to use software, even if the
29     * University of Notre Dame has been advised of the possibility of
30     * such damages.
31 gezelter 1390 *
32     * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33     * research, please cite the appropriate papers when you publish your
34     * work. Good starting points are:
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).
39     * [4] Vardeman & Gezelter, in progress (2009).
40 gezelter 246 */
41 mciznick 1598
42 gezelter 507 /**
43     * @file ForceManager.cpp
44     * @author tlin
45     * @date 11/09/2004
46     * @time 10:39am
47     * @version 1.0
48     */
49 gezelter 246
50     #include "brains/ForceManager.hpp"
51     #include "primitives/Molecule.hpp"
52 gezelter 1390 #define __OPENMD_C
53 gezelter 246 #include "utils/simError.h"
54 xsun 1215 #include "primitives/Bond.hpp"
55 tim 749 #include "primitives/Bend.hpp"
56 cli2 1275 #include "primitives/Torsion.hpp"
57     #include "primitives/Inversion.hpp"
58 gezelter 1551 #include "nonbonded/NonBondedInteraction.hpp"
59 gezelter 1549 #include "parallel/ForceMatrixDecomposition.hpp"
60 gezelter 1467
61 gezelter 1583 #include <cstdio>
62     #include <iostream>
63     #include <iomanip>
64    
65 mciznick 1598 #include <omp.h>
66    
67 gezelter 1545 using namespace std;
68 gezelter 1390 namespace OpenMD {
69 gezelter 1576
70 mciznick 1598 ForceManager::ForceManager(SimInfo * info) :
71     info_(info) {
72     forceField_ = info_->getForceField();
73     interactionMan_ = new InteractionManager();
74     fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_);
75     }
76 gezelter 1576
77 mciznick 1598 /**
78     * setupCutoffs
79     *
80     * Sets the values of cutoffRadius, switchingRadius, cutoffMethod,
81     * and cutoffPolicy
82     *
83     * cutoffRadius : realType
84     * If the cutoffRadius was explicitly set, use that value.
85     * If the cutoffRadius was not explicitly set:
86     * Are there electrostatic atoms? Use 12.0 Angstroms.
87     * No electrostatic atoms? Poll the atom types present in the
88     * simulation for suggested cutoff values (e.g. 2.5 * sigma).
89     * Use the maximum suggested value that was found.
90     *
91     * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE,
92     * or SHIFTED_POTENTIAL)
93     * If cutoffMethod was explicitly set, use that choice.
94     * If cutoffMethod was not explicitly set, use SHIFTED_FORCE
95     *
96     * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
97     * If cutoffPolicy was explicitly set, use that choice.
98     * If cutoffPolicy was not explicitly set, use TRADITIONAL
99     *
100     * switchingRadius : realType
101     * If the cutoffMethod was set to SWITCHED:
102     * If the switchingRadius was explicitly set, use that value
103     * (but do a sanity check first).
104     * If the switchingRadius was not explicitly set: use 0.85 *
105     * cutoffRadius_
106     * If the cutoffMethod was not set to SWITCHED:
107     * Set switchingRadius equal to cutoffRadius for safety.
108     */
109     void ForceManager::setupCutoffs() {
110 gezelter 1583
111 mciznick 1598 Globals* simParams_ = info_->getSimParams();
112     ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
113 gezelter 1576
114 mciznick 1598 if (simParams_->haveCutoffRadius())
115     {
116     rCut_ = simParams_->getCutoffRadius();
117     } else
118     {
119     if (info_->usesElectrostaticAtoms())
120     {
121     sprintf(painCave.errMsg, "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
122     "\tOpenMD will use a default value of 12.0 angstroms"
123     "\tfor the cutoffRadius.\n");
124     painCave.isFatal = 0;
125     painCave.severity = OPENMD_INFO;
126     simError();
127     rCut_ = 12.0;
128     } else
129     {
130     RealType thisCut;
131     set<AtomType*>::iterator i;
132     set<AtomType*> atomTypes;
133     atomTypes = info_->getSimulatedAtomTypes();
134     for (i = atomTypes.begin(); i != atomTypes.end(); ++i)
135     {
136     thisCut = interactionMan_->getSuggestedCutoffRadius((*i));
137     rCut_ = max(thisCut, rCut_);
138     }
139     sprintf(painCave.errMsg, "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
140     "\tOpenMD will use %lf angstroms.\n", rCut_);
141     painCave.isFatal = 0;
142     painCave.severity = OPENMD_INFO;
143     simError();
144     }
145     }
146 gezelter 1576
147 mciznick 1598 fDecomp_->setUserCutoff(rCut_);
148     interactionMan_->setCutoffRadius(rCut_);
149 gezelter 1576
150 mciznick 1598 map<string, CutoffMethod> stringToCutoffMethod;
151     stringToCutoffMethod["HARD"] = HARD;
152     stringToCutoffMethod["SWITCHED"] = SWITCHED;
153     stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;
154     stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
155 gezelter 1576
156 mciznick 1598 if (simParams_->haveCutoffMethod())
157     {
158     string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
159     map<string, CutoffMethod>::iterator i;
160     i = stringToCutoffMethod.find(cutMeth);
161     if (i == stringToCutoffMethod.end())
162     {
163     sprintf(painCave.errMsg, "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
164     "\tShould be one of: "
165     "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", cutMeth.c_str());
166     painCave.isFatal = 1;
167     painCave.severity = OPENMD_ERROR;
168     simError();
169     } else
170     {
171     cutoffMethod_ = i->second;
172     }
173     } else
174     {
175     sprintf(painCave.errMsg, "ForceManager::setupCutoffs: No value was set for the cutoffMethod.\n"
176     "\tOpenMD will use SHIFTED_FORCE.\n");
177     painCave.isFatal = 0;
178     painCave.severity = OPENMD_INFO;
179     simError();
180     cutoffMethod_ = SHIFTED_FORCE;
181     }
182 gezelter 1587
183 mciznick 1598 map<string, CutoffPolicy> stringToCutoffPolicy;
184     stringToCutoffPolicy["MIX"] = MIX;
185     stringToCutoffPolicy["MAX"] = MAX;
186     stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;
187 gezelter 1576
188 mciznick 1598 std::string cutPolicy;
189     if (forceFieldOptions_.haveCutoffPolicy())
190     {
191     cutPolicy = forceFieldOptions_.getCutoffPolicy();
192     } else if (simParams_->haveCutoffPolicy())
193     {
194     cutPolicy = simParams_->getCutoffPolicy();
195     }
196 gezelter 1587
197 mciznick 1598 if (!cutPolicy.empty())
198     {
199     toUpper(cutPolicy);
200     map<string, CutoffPolicy>::iterator i;
201     i = stringToCutoffPolicy.find(cutPolicy);
202 gezelter 1576
203 mciznick 1598 if (i == stringToCutoffPolicy.end())
204     {
205     sprintf(painCave.errMsg, "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
206     "\tShould be one of: "
207     "MIX, MAX, or TRADITIONAL\n", cutPolicy.c_str());
208     painCave.isFatal = 1;
209     painCave.severity = OPENMD_ERROR;
210     simError();
211     } else
212     {
213     cutoffPolicy_ = i->second;
214     }
215     } else
216     {
217     sprintf(painCave.errMsg, "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
218     "\tOpenMD will use TRADITIONAL.\n");
219     painCave.isFatal = 0;
220     painCave.severity = OPENMD_INFO;
221     simError();
222     cutoffPolicy_ = TRADITIONAL;
223     }
224 gezelter 1590
225 mciznick 1598 fDecomp_->setCutoffPolicy(cutoffPolicy_);
226 gezelter 1576
227 mciznick 1598 // create the switching function object:
228 gezelter 1576
229 mciznick 1598 switcher_ = new SwitchingFunction();
230 gezelter 1576
231 mciznick 1598 if (cutoffMethod_ == SWITCHED)
232     {
233     if (simParams_->haveSwitchingRadius())
234     {
235     rSwitch_ = simParams_->getSwitchingRadius();
236     if (rSwitch_ > rCut_)
237     {
238     sprintf(painCave.errMsg, "ForceManager::setupCutoffs: switchingRadius (%f) is larger "
239     "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
240     painCave.isFatal = 1;
241     painCave.severity = OPENMD_ERROR;
242     simError();
243     }
244     } else
245     {
246     rSwitch_ = 0.85 * rCut_;
247     sprintf(painCave.errMsg, "ForceManager::setupCutoffs: No value was set for the switchingRadius.\n"
248     "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
249     "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
250     painCave.isFatal = 0;
251     painCave.severity = OPENMD_WARNING;
252     simError();
253     }
254     } else
255     {
256     if (simParams_->haveSwitchingRadius())
257     {
258     map<string, CutoffMethod>::const_iterator it;
259     string theMeth;
260     for (it = stringToCutoffMethod.begin(); it != stringToCutoffMethod.end(); ++it)
261     {
262     if (it->second == cutoffMethod_)
263     {
264     theMeth = it->first;
265     break;
266     }
267     }
268     sprintf(painCave.errMsg, "ForceManager::setupCutoffs: the cutoffMethod (%s)\n"
269     "\tis not set to SWITCHED, so switchingRadius value\n"
270     "\twill be ignored for this simulation\n", theMeth.c_str());
271     painCave.isFatal = 0;
272     painCave.severity = OPENMD_WARNING;
273     simError();
274     }
275 gezelter 1576
276 mciznick 1598 rSwitch_ = rCut_;
277     }
278 gezelter 1576
279 mciznick 1598 // Default to cubic switching function.
280     sft_ = cubic;
281     if (simParams_->haveSwitchingFunctionType())
282     {
283     string funcType = simParams_->getSwitchingFunctionType();
284     toUpper(funcType);
285     if (funcType == "CUBIC")
286     {
287     sft_ = cubic;
288     } else
289     {
290     if (funcType == "FIFTH_ORDER_POLYNOMIAL")
291     {
292     sft_ = fifth_order_poly;
293     } else
294     {
295     // throw error
296     sprintf(painCave.errMsg,
297     "ForceManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n"
298     "\tswitchingFunctionType must be one of: "
299     "\"cubic\" or \"fifth_order_polynomial\".", funcType.c_str());
300     painCave.isFatal = 1;
301     painCave.severity = OPENMD_ERROR;
302     simError();
303     }
304     }
305     }
306     switcher_->setSwitchType(sft_);
307     switcher_->setSwitch(rSwitch_, rCut_);
308     interactionMan_->setSwitchingRadius(rSwitch_);
309     }
310 gezelter 1576
311 mciznick 1598 void ForceManager::initialize() {
312 gezelter 1576
313 mciznick 1598 if (!info_->isTopologyDone())
314     {
315 gezelter 1576
316 mciznick 1598 info_->update();
317     interactionMan_->setSimInfo(info_);
318     interactionMan_->initialize();
319 gezelter 246
320 mciznick 1598 // We want to delay the cutoffs until after the interaction
321     // manager has set up the atom-atom interactions so that we can
322     // query them for suggested cutoff values
323     setupCutoffs();
324 gezelter 246
325 mciznick 1598 info_->prepareTopology();
326     }
327 gezelter 246
328 mciznick 1598 ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
329 gezelter 246
330 mciznick 1598 // Force fields can set options on how to scale van der Waals and
331     // electrostatic interactions for atoms connected via bonds, bends
332     // and torsions in this case the topological distance between
333     // atoms is:
334     // 0 = topologically unconnected
335     // 1 = bonded together
336     // 2 = connected via a bend
337     // 3 = connected via a torsion
338 chuckv 1595
339 mciznick 1598 vdwScale_.reserve(4);
340     fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
341 chuckv 1595
342 mciznick 1598 electrostaticScale_.reserve(4);
343     fill(electrostaticScale_.begin(), electrostaticScale_.end(), 0.0);
344 chuckv 1595
345 mciznick 1598 vdwScale_[0] = 1.0;
346     vdwScale_[1] = fopts.getvdw12scale();
347     vdwScale_[2] = fopts.getvdw13scale();
348     vdwScale_[3] = fopts.getvdw14scale();
349 chuckv 1595
350 mciznick 1598 electrostaticScale_[0] = 1.0;
351     electrostaticScale_[1] = fopts.getelectrostatic12scale();
352     electrostaticScale_[2] = fopts.getelectrostatic13scale();
353     electrostaticScale_[3] = fopts.getelectrostatic14scale();
354 chuckv 1595
355 mciznick 1598 fDecomp_->distributeInitialData();
356 chuckv 1595
357 mciznick 1598 initialized_ = true;
358 chuckv 1595
359 mciznick 1598 }
360 chuckv 1595
361 mciznick 1598 void ForceManager::calcForces() {
362 chuckv 1595
363 mciznick 1598 if (!initialized_)
364     initialize();
365 chuckv 1595
366 mciznick 1598 preCalculation();
367     shortRangeInteractions();
368     // longRangeInteractions();
369 mciznick 1608 // longRangeInteractionsRapaport();
370     longRangeInteractionsParallel();
371 mciznick 1598 postCalculation();
372     }
373 chuckv 1595
374 mciznick 1598 void ForceManager::preCalculation() {
375     SimInfo::MoleculeIterator mi;
376     Molecule* mol;
377     Molecule::AtomIterator ai;
378     Atom* atom;
379     Molecule::RigidBodyIterator rbIter;
380     RigidBody* rb;
381     Molecule::CutoffGroupIterator ci;
382     CutoffGroup* cg;
383 chuckv 1595
384 mciznick 1598 // forces are zeroed here, before any are accumulated.
385 chuckv 1595
386 mciznick 1598 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
387     {
388     for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai))
389     {
390     atom->zeroForcesAndTorques();
391     }
392 chuckv 1595
393 mciznick 1598 //change the positions of atoms which belong to the rigidbodies
394     for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
395     {
396     rb->zeroForcesAndTorques();
397     }
398 chuckv 1595
399 mciznick 1598 if (info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms())
400     {
401     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
402     {
403     //calculate the center of mass of cutoff group
404     cg->updateCOM();
405     }
406     }
407     }
408 chuckv 1595
409 mciznick 1598 // Zero out the stress tensor
410     tau *= 0.0;
411 chuckv 1595
412 mciznick 1598 }
413 chuckv 1595
414 mciznick 1598 void ForceManager::shortRangeInteractions() {
415     Molecule* mol;
416     RigidBody* rb;
417     Bond* bond;
418     Bend* bend;
419     Torsion* torsion;
420     Inversion* inversion;
421     SimInfo::MoleculeIterator mi;
422     Molecule::RigidBodyIterator rbIter;
423     Molecule::BondIterator bondIter;
424     ;
425     Molecule::BendIterator bendIter;
426     Molecule::TorsionIterator torsionIter;
427     Molecule::InversionIterator inversionIter;
428     RealType bondPotential = 0.0;
429     RealType bendPotential = 0.0;
430     RealType torsionPotential = 0.0;
431     RealType inversionPotential = 0.0;
432 chuckv 1595
433 mciznick 1598 //calculate short range interactions
434     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
435     {
436 chuckv 1595
437 mciznick 1598 //change the positions of atoms which belong to the rigidbodies
438     for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
439     {
440     rb->updateAtoms();
441     }
442 chuckv 1595
443 mciznick 1598 for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter))
444     {
445     bond->calcForce();
446     bondPotential += bond->getPotential();
447     }
448 chuckv 1595
449 mciznick 1598 for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter))
450     {
451 chuckv 1595
452 mciznick 1598 RealType angle;
453     bend->calcForce(angle);
454     RealType currBendPot = bend->getPotential();
455 chuckv 1595
456 mciznick 1598 bendPotential += bend->getPotential();
457     map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
458     if (i == bendDataSets.end())
459     {
460     BendDataSet dataSet;
461     dataSet.prev.angle = dataSet.curr.angle = angle;
462     dataSet.prev.potential = dataSet.curr.potential = currBendPot;
463     dataSet.deltaV = 0.0;
464     bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
465     } else
466     {
467     i->second.prev.angle = i->second.curr.angle;
468     i->second.prev.potential = i->second.curr.potential;
469     i->second.curr.angle = angle;
470     i->second.curr.potential = currBendPot;
471     i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
472     }
473     }
474 chuckv 1595
475 mciznick 1598 for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter))
476     {
477     RealType angle;
478     torsion->calcForce(angle);
479     RealType currTorsionPot = torsion->getPotential();
480     torsionPotential += torsion->getPotential();
481     map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
482     if (i == torsionDataSets.end())
483     {
484     TorsionDataSet dataSet;
485     dataSet.prev.angle = dataSet.curr.angle = angle;
486     dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
487     dataSet.deltaV = 0.0;
488     torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
489     } else
490     {
491     i->second.prev.angle = i->second.curr.angle;
492     i->second.prev.potential = i->second.curr.potential;
493     i->second.curr.angle = angle;
494     i->second.curr.potential = currTorsionPot;
495     i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
496     }
497     }
498 chuckv 1595
499 mciznick 1598 for (inversion = mol->beginInversion(inversionIter); inversion != NULL; inversion = mol->nextInversion(
500     inversionIter))
501     {
502     RealType angle;
503     inversion->calcForce(angle);
504     RealType currInversionPot = inversion->getPotential();
505     inversionPotential += inversion->getPotential();
506     map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
507     if (i == inversionDataSets.end())
508     {
509     InversionDataSet dataSet;
510     dataSet.prev.angle = dataSet.curr.angle = angle;
511     dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
512     dataSet.deltaV = 0.0;
513     inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
514     } else
515     {
516     i->second.prev.angle = i->second.curr.angle;
517     i->second.prev.potential = i->second.curr.potential;
518     i->second.curr.angle = angle;
519     i->second.curr.potential = currInversionPot;
520     i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
521     }
522     }
523     }
524 chuckv 1595
525 mciznick 1598 RealType shortRangePotential = bondPotential + bendPotential + torsionPotential + inversionPotential;
526     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
527     curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
528     curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
529     curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
530     curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
531     curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
532     }
533 chuckv 1595
534 mciznick 1608 void ForceManager::longRangeInteractionsParallel() {
535    
536     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
537     DataStorage* config = &(curSnapshot->atomData);
538     DataStorage* cgConfig = &(curSnapshot->cgData);
539    
540     //calculate the center of mass of cutoff group
541    
542     SimInfo::MoleculeIterator mi;
543     Molecule* mol;
544     Molecule::CutoffGroupIterator ci;
545     CutoffGroup* cg;
546    
547     if (info_->getNCutoffGroups() > 0)
548     {
549     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
550     {
551     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
552     {
553     // cerr << "branch1\n";
554     // cerr << "globind = " << cg->getGlobalIndex() << ":" << __LINE__ << "\n";
555     cg->updateCOM();
556    
557     // cerr << "gbI: " << cg->getGlobalIndex() << " locI: " << cg->getLocalIndex() << " x: "
558     // << cgConfig->position[cg->getLocalIndex()].x() << " y: " << cgConfig->position[cg->getLocalIndex()].y()
559     // << " z: " << cgConfig->position[cg->getLocalIndex()].z() << "\n";
560     }
561     }
562     } else
563     {
564     // center of mass of the group is the same as position of the atom
565     // if cutoff group does not exist
566     // cerr << ":" << __LINE__ << "branch2\n";
567     cgConfig->position = config->position;
568     }
569    
570     fDecomp_->zeroWorkArrays();
571     fDecomp_->distributeData();
572    
573     int atom1, atom2, topoDist;
574     Vector3d d_grp, dag, d;
575     RealType rgrpsq, rgrp, r2, r;
576     RealType electroMult, vdwMult;
577     RealType vij;
578     Vector3d fij, fg, f1;
579     tuple3<RealType, RealType, RealType> cuts;
580     RealType rCutSq;
581     bool in_switching_region;
582     RealType sw, dswdr, swderiv;
583     vector<int> atomListColumn, atomListRow, atomListLocal;
584    
585     InteractionDataPrv idatPrv;
586    
587     SelfData sdat;
588     RealType mf;
589     RealType lrPot;
590     RealType vpair;
591     potVec longRangePotential(0.0);
592     potVec workPot(0.0);
593    
594     int loopStart, loopEnd;
595     sdat.pot = fDecomp_->getEmbeddingPotential();
596    
597     vector<CutoffGroup *> cgs;
598    
599     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
600     {
601     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
602     {
603     cgs.push_back(cg);
604     }
605     }
606    
607     loopEnd = PAIR_LOOP;
608     if (info_->requiresPrepair())
609     {
610     loopStart = PREPAIR_LOOP;
611     } else
612     {
613     loopStart = PAIR_LOOP;
614     }
615    
616     for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++)
617     {
618    
619     if (iLoop == loopStart)
620     {
621     bool update_nlist = fDecomp_->checkNeighborList();
622     if (update_nlist)
623     neighborMatW = fDecomp_->buildLayerBasedNeighborList();
624     }
625    
626     vector<CutoffGroup *>::iterator cg1;
627     vector<CutoffGroup *>::iterator cg2;
628    
629     // int nThreads = 2;
630     int chunkSize = cgs.size() / (omp_get_num_threads() * 20);
631    
632     // printf("before omp loop\n");
633     #pragma omp parallel /*num_threads(nThreads)*/ default(none) shared(curSnapshot, iLoop, cgs, chunkSize) \
634     private(cg1, cg2, cuts, d_grp, rgrpsq, rCutSq, idatPrv, vij, fij, in_switching_region, dswdr, rgrp, \
635     atomListRow, atomListColumn, atom1, atom2, topoDist, d, r2, swderiv, fg, mf, dag)
636     {
637     idatPrv.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
638     idatPrv.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
639    
640     // printf("Thread %d\n", omp_get_thread_num());
641     #pragma omp for schedule(dynamic, chunkSize)
642     for (cg1 = cgs.begin(); cg1 < cgs.end(); ++cg1)
643     {
644     for (cg2 = neighborMatW[(*cg1)->getGlobalIndex()].begin(); cg2 < neighborMatW[(*cg1)->getGlobalIndex()].end(); ++cg2)
645     {
646    
647     cuts = fDecomp_->getGroupCutoffs((*cg1)->getGlobalIndex(), (*cg2)->getGlobalIndex());
648    
649     d_grp = fDecomp_->getIntergroupVector((*cg1), (*cg2));
650     curSnapshot->wrapVector(d_grp);
651     rgrpsq = d_grp.lengthSquare();
652    
653     rCutSq = cuts.second;
654    
655     // printf("Thread %d\tcg1:%d\tcg2:%d d_grp\tx:%f\ty:%f\tz:%f\trgrpsq:%f\n", omp_get_thread_num(), (*cg1)->getGlobalIndex(), (*cg2)->getGlobalIndex(), d_grp.x(), d_grp.y(), d_grp.z(), rgrpsq);
656    
657     if (rgrpsq < rCutSq)
658     {
659     idatPrv.rcut = cuts.first;
660     if (iLoop == PAIR_LOOP)
661     {
662     vij = 0.0;
663     fij = V3Zero;
664     }
665    
666     in_switching_region = switcher_->getSwitch(rgrpsq, /*sw*/idatPrv.sw, dswdr, rgrp);
667    
668     // printf("in_switching_region:%d\trgrpsq:%f\t*idatPrv.sw:%f\tdswdr:%f\trgrp:%f\n", (in_switching_region == false ? 0 : 1), rgrpsq, idatPrv.sw, dswdr, rgrp);
669    
670     atomListRow = fDecomp_->getAtomsInGroupRow((*cg1)->getGlobalIndex());
671     atomListColumn = fDecomp_->getAtomsInGroupColumn((*cg2)->getGlobalIndex());
672    
673     for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
674     {
675     atom1 = (*ia);
676    
677     for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
678     {
679     atom2 = (*jb);
680    
681     // printf("atom1:%d atom2:%d\n", atom1, atom2);
682     if (!fDecomp_->skipAtomPair(atom1, atom2))
683     {
684     idatPrv.vpair = 0.0;
685     idatPrv.pot = 0.0;
686     idatPrv.f1 = V3Zero;
687    
688     fDecomp_->fillInteractionDataOMP(idatPrv, atom1, atom2);
689    
690     topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
691     idatPrv.vdwMult = vdwScale_[topoDist];
692     idatPrv.electroMult = electrostaticScale_[topoDist];
693    
694     // printf("topoDist:%d\tidatPrv.vdwMult:%f\tidatPrv.electroMult:%f\n", topoDist, idatPrv.vdwMult, idatPrv.electroMult);
695    
696     if (atomListRow.size() == 1 && atomListColumn.size() == 1)
697     {
698     idatPrv.d = d_grp;
699     idatPrv.r2 = rgrpsq;
700     // cerr << "dgrp = " << d_grp << ":" << __LINE__ << "\n";
701     } else
702     {
703     d = fDecomp_->getInteratomicVector(atom1, atom2);
704     curSnapshot->wrapVector(d);
705     r2 = d.lengthSquare();
706     // cerr << "datm = " << d << ":" << __LINE__ << "\n";
707     idatPrv.d = d;
708     idatPrv.r2 = r2;
709     }
710    
711     // printf("idatPrv.d x:%f\ty:%f\tz:%f\tidatPrv.r2:%f\n", (idatPrv.d).x(), (idatPrv.d).y(), (idatPrv.d).z(), idatPrv.r2);
712    
713     // cerr << "idat.d = " << *(idat.d) << ":" << __LINE__ << "\n";
714     idatPrv.rij = sqrt((idatPrv.r2));
715     // cerr << "idat.rij = " << *(idat.rij) << "\n";
716    
717     #pragma omp critical
718     {
719     interactionMan_->initializeOMP();
720     }
721    
722     if (iLoop == PREPAIR_LOOP)
723     {
724     interactionMan_->doPrePairOMP(idatPrv);
725     } else
726     {
727     interactionMan_->doPairOMP(idatPrv);
728     fDecomp_->unpackInteractionDataOMP(idatPrv, atom1, atom2);
729    
730     // cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << ":" << __LINE__ << "\n";
731     // printf("d x:%f y:%f z:%f vpair:%f f1 x:%f y:%f z:%f\n", idatPrv.d.x(), idatPrv.d.y(), idatPrv.d.z(), idatPrv.vpair, idatPrv.f1.x(), idatPrv.f1.y(), idatPrv.f1.z());
732     #pragma omp critical
733     {
734     vij += idatPrv.vpair;
735     fij += idatPrv.f1;
736     tau -= outProduct(idatPrv.d, idatPrv.f1);
737    
738     // printf("vij:%f fij x:%f y:%f z:%f\n", vij, fij.x(), fij.y(), fij.z());
739     }
740     }
741     }
742     }
743     }
744    
745     if (iLoop == PAIR_LOOP)
746     {
747     if (in_switching_region)
748     {
749     swderiv = vij * dswdr / rgrp;
750     fg = swderiv * d_grp;
751     fij += fg;
752    
753     if (atomListRow.size() == 1 && atomListColumn.size() == 1)
754     {
755     #pragma omp critical
756     {
757     tau -= outProduct(idatPrv.d, fg);
758     }
759     }
760    
761     for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
762     {
763     atom1 = (*ia);
764     mf = fDecomp_->getMassFactorRow(atom1);
765     // fg is the force on atom ia due to cutoff group's
766     // presence in switching region
767     fg = swderiv * d_grp * mf;
768     fDecomp_->addForceToAtomRowOMP(atom1, fg);
769    
770     if (atomListRow.size() > 1)
771     {
772     if (info_->usesAtomicVirial())
773     {
774     // find the distance between the atom
775     // and the center of the cutoff group:
776     dag = fDecomp_->getAtomToGroupVectorRow(atom1, (*cg1)->getGlobalIndex());
777     #pragma omp critical
778     {
779     tau -= outProduct(dag, fg);
780     }
781     }
782     }
783     }
784     for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
785     {
786     atom2 = (*jb);
787     mf = fDecomp_->getMassFactorColumn(atom2);
788     // fg is the force on atom jb due to cutoff group's
789     // presence in switching region
790     fg = -swderiv * d_grp * mf;
791     fDecomp_->addForceToAtomColumn(atom2, fg);
792    
793     if (atomListColumn.size() > 1)
794     {
795     if (info_->usesAtomicVirial())
796     {
797     // find the distance between the atom
798     // and the center of the cutoff group:
799     dag = fDecomp_->getAtomToGroupVectorColumn(atom2, (*cg2)->getGlobalIndex());
800     #pragma omp critical
801     {
802     tau -= outProduct(dag, fg);
803     }
804     }
805     }
806     }
807     }
808     //if (!SIM_uses_AtomicVirial) {
809     // tau -= outProduct(d_grp, fij);
810     //}
811     }
812     }
813     }
814     }// END: omp for loop
815     // printf("after omp loop\n");
816     }
817    
818     if (iLoop == PREPAIR_LOOP)
819     {
820     if (info_->requiresPrepair())
821     {
822    
823     fDecomp_->collectIntermediateData();
824    
825     for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
826     {
827     fDecomp_->fillSelfData(sdat, atom1);
828     interactionMan_->doPreForce(sdat);
829     }
830    
831     fDecomp_->distributeIntermediateData();
832    
833     }
834     }
835     }
836    
837     fDecomp_->collectData();
838    
839     if (info_->requiresSelfCorrection())
840     {
841    
842     for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
843     {
844     fDecomp_->fillSelfData(sdat, atom1);
845     interactionMan_->doSelfCorrection(sdat);
846     }
847    
848     }
849    
850     longRangePotential = *(fDecomp_->getEmbeddingPotential()) + *(fDecomp_->getPairwisePotential());
851    
852     lrPot = longRangePotential.sum();
853    
854     //store the tau and long range potential
855     curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
856     curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
857     curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
858     }
859    
860 mciznick 1598 void ForceManager::longRangeInteractionsRapaport() {
861 chuckv 1595
862 mciznick 1598 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
863     DataStorage* config = &(curSnapshot->atomData);
864     DataStorage* cgConfig = &(curSnapshot->cgData);
865 chuckv 1595
866 mciznick 1598 //calculate the center of mass of cutoff group
867 chuckv 1595
868 mciznick 1598 SimInfo::MoleculeIterator mi;
869     Molecule* mol;
870     Molecule::CutoffGroupIterator ci;
871     CutoffGroup* cg;
872 chuckv 1595
873 mciznick 1598 if (info_->getNCutoffGroups() > 0)
874     {
875     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
876     {
877     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
878     {
879 mciznick 1608 // cerr << "branch1\n";
880     // cerr << "globind = " << cg->getGlobalIndex() << ":" << __LINE__ << "\n";
881 mciznick 1598 cg->updateCOM();
882    
883 mciznick 1608 // cerr << "gbI: " << cg->getGlobalIndex() << " locI: " << cg->getLocalIndex() << " x: "
884     // << cgConfig->position[cg->getLocalIndex()].x() << " y: " << cgConfig->position[cg->getLocalIndex()].y()
885     // << " z: " << cgConfig->position[cg->getLocalIndex()].z() << "\n";
886 mciznick 1598 }
887     }
888     } else
889     {
890     // center of mass of the group is the same as position of the atom
891     // if cutoff group does not exist
892 mciznick 1608 // cerr << ":" << __LINE__ << "branch2\n";
893 mciznick 1598 cgConfig->position = config->position;
894     }
895    
896     fDecomp_->zeroWorkArrays();
897     fDecomp_->distributeData();
898    
899     int atom1, atom2, topoDist;
900     CutoffGroup *cg1;
901     Vector3d d_grp, dag, d;
902     RealType rgrpsq, rgrp, r2, r;
903     RealType electroMult, vdwMult;
904     RealType vij;
905     Vector3d fij, fg, f1;
906     tuple3<RealType, RealType, RealType> cuts;
907     RealType rCutSq;
908     bool in_switching_region;
909     RealType sw, dswdr, swderiv;
910     vector<int> atomListColumn, atomListRow, atomListLocal;
911     InteractionData idat;
912     SelfData sdat;
913     RealType mf;
914     RealType lrPot;
915     RealType vpair;
916     potVec longRangePotential(0.0);
917     potVec workPot(0.0);
918    
919     int loopStart, loopEnd;
920    
921     idat.vdwMult = &vdwMult;
922     idat.electroMult = &electroMult;
923     idat.pot = &workPot;
924     sdat.pot = fDecomp_->getEmbeddingPotential();
925     idat.vpair = &vpair;
926     idat.f1 = &f1;
927     idat.sw = &sw;
928     idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
929     idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
930    
931     loopEnd = PAIR_LOOP;
932     if (info_->requiresPrepair())
933     {
934     loopStart = PREPAIR_LOOP;
935     } else
936     {
937     loopStart = PAIR_LOOP;
938     }
939    
940     for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++)
941     {
942    
943     if (iLoop == loopStart)
944     {
945     bool update_nlist = fDecomp_->checkNeighborList();
946     if (update_nlist)
947     neighborMatW = fDecomp_->buildLayerBasedNeighborList();
948     }
949    
950     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
951     {
952     for (cg1 = mol->beginCutoffGroup(ci); cg1 != NULL; cg1 = mol->nextCutoffGroup(ci))
953     {
954     // printf("Thread %d executes loop iteration %d\n", omp_get_thread_num(), i);
955 mciznick 1608 for (vector<CutoffGroup *>::iterator cg2 = neighborMatW[cg1->getGlobalIndex()].begin(); cg2
956     != neighborMatW[cg1->getGlobalIndex()].end(); ++cg2)
957 mciznick 1598 {
958    
959     cuts = fDecomp_->getGroupCutoffs(cg1->getGlobalIndex(), (*cg2)->getGlobalIndex());
960    
961     d_grp = fDecomp_->getIntergroupVector(cg1, (*cg2));
962     curSnapshot->wrapVector(d_grp);
963     rgrpsq = d_grp.lengthSquare();
964    
965     rCutSq = cuts.second;
966    
967     if (rgrpsq < rCutSq)
968     {
969     idat.rcut = &cuts.first;
970     if (iLoop == PAIR_LOOP)
971     {
972     vij = 0.0;
973     fij = V3Zero;
974     }
975    
976     in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, rgrp);
977    
978     atomListRow = fDecomp_->getAtomsInGroupRow(cg1->getGlobalIndex());
979     atomListColumn = fDecomp_->getAtomsInGroupColumn((*cg2)->getGlobalIndex());
980    
981     for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
982     {
983     atom1 = (*ia);
984    
985     for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
986     {
987     atom2 = (*jb);
988    
989     if (!fDecomp_->skipAtomPair(atom1, atom2))
990     {
991     vpair = 0.0;
992     workPot = 0.0;
993     f1 = V3Zero;
994    
995     fDecomp_->fillInteractionData(idat, atom1, atom2);
996    
997     topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
998     vdwMult = vdwScale_[topoDist];
999     electroMult = electrostaticScale_[topoDist];
1000    
1001     if (atomListRow.size() == 1 && atomListColumn.size() == 1)
1002     {
1003     idat.d = &d_grp;
1004     idat.r2 = &rgrpsq;
1005 mciznick 1608 // cerr << "dgrp = " << d_grp << ":" << __LINE__ << "\n";
1006 mciznick 1598 } else
1007     {
1008     d = fDecomp_->getInteratomicVector(atom1, atom2);
1009     curSnapshot->wrapVector(d);
1010     r2 = d.lengthSquare();
1011 mciznick 1608 // cerr << "datm = " << d << ":" << __LINE__ << "\n";
1012 mciznick 1598 idat.d = &d;
1013     idat.r2 = &r2;
1014     }
1015    
1016 mciznick 1608 // cerr << "idat.d = " << *(idat.d) << ":" << __LINE__ << "\n";
1017 mciznick 1598 r = sqrt(*(idat.r2));
1018     idat.rij = &r;
1019 mciznick 1608 // cerr << "idat.rij = " << *(idat.rij) << "\n";
1020 mciznick 1598
1021     if (iLoop == PREPAIR_LOOP)
1022     {
1023     interactionMan_->doPrePair(idat);
1024     } else
1025     {
1026     interactionMan_->doPair(idat);
1027     fDecomp_->unpackInteractionData(idat, atom1, atom2);
1028    
1029 mciznick 1608 // cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << ":" << __LINE__ << "\n";
1030 mciznick 1598 vij += vpair;
1031     fij += f1;
1032     tau -= outProduct(*(idat.d), f1);
1033     }
1034     }
1035     }
1036     }
1037    
1038     if (iLoop == PAIR_LOOP)
1039     {
1040     if (in_switching_region)
1041     {
1042     swderiv = vij * dswdr / rgrp;
1043     fg = swderiv * d_grp;
1044     fij += fg;
1045    
1046     if (atomListRow.size() == 1 && atomListColumn.size() == 1)
1047     {
1048     tau -= outProduct(*(idat.d), fg);
1049     }
1050    
1051     for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
1052     {
1053     atom1 = (*ia);
1054     mf = fDecomp_->getMassFactorRow(atom1);
1055     // fg is the force on atom ia due to cutoff group's
1056     // presence in switching region
1057     fg = swderiv * d_grp * mf;
1058     fDecomp_->addForceToAtomRow(atom1, fg);
1059    
1060     if (atomListRow.size() > 1)
1061     {
1062     if (info_->usesAtomicVirial())
1063     {
1064     // find the distance between the atom
1065     // and the center of the cutoff group:
1066     dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1->getGlobalIndex());
1067     tau -= outProduct(dag, fg);
1068     }
1069     }
1070     }
1071     for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
1072     {
1073     atom2 = (*jb);
1074     mf = fDecomp_->getMassFactorColumn(atom2);
1075     // fg is the force on atom jb due to cutoff group's
1076     // presence in switching region
1077     fg = -swderiv * d_grp * mf;
1078     fDecomp_->addForceToAtomColumn(atom2, fg);
1079    
1080     if (atomListColumn.size() > 1)
1081     {
1082     if (info_->usesAtomicVirial())
1083     {
1084     // find the distance between the atom
1085     // and the center of the cutoff group:
1086     dag = fDecomp_->getAtomToGroupVectorColumn(atom2, (*cg2)->getGlobalIndex());
1087     tau -= outProduct(dag, fg);
1088     }
1089     }
1090     }
1091     }
1092     //if (!SIM_uses_AtomicVirial) {
1093     // tau -= outProduct(d_grp, fij);
1094     //}
1095     }
1096 chuckv 1595 }
1097     }
1098     }
1099 mciznick 1608 }
1100 chuckv 1595
1101 mciznick 1598 if (iLoop == PREPAIR_LOOP)
1102     {
1103     if (info_->requiresPrepair())
1104     {
1105 chuckv 1595
1106 mciznick 1598 fDecomp_->collectIntermediateData();
1107 chuckv 1595
1108 mciznick 1598 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1109     {
1110     fDecomp_->fillSelfData(sdat, atom1);
1111     interactionMan_->doPreForce(sdat);
1112     }
1113 chuckv 1595
1114 mciznick 1598 fDecomp_->distributeIntermediateData();
1115 chuckv 1595
1116 mciznick 1598 }
1117     }
1118     }
1119 chuckv 1595
1120 mciznick 1598 fDecomp_->collectData();
1121 chuckv 1595
1122 mciznick 1598 if (info_->requiresSelfCorrection())
1123     {
1124 chuckv 1595
1125 mciznick 1598 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1126     {
1127     fDecomp_->fillSelfData(sdat, atom1);
1128     interactionMan_->doSelfCorrection(sdat);
1129     }
1130 chuckv 1595
1131 mciznick 1598 }
1132 chuckv 1595
1133 mciznick 1598 longRangePotential = *(fDecomp_->getEmbeddingPotential()) + *(fDecomp_->getPairwisePotential());
1134 chuckv 1595
1135 mciznick 1598 lrPot = longRangePotential.sum();
1136 chuckv 1595
1137 mciznick 1598 //store the tau and long range potential
1138     curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
1139     curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
1140     curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
1141     }
1142 chuckv 1595
1143 mciznick 1598 void ForceManager::longRangeInteractions() {
1144 chuckv 1595
1145 mciznick 1598 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1146     DataStorage* config = &(curSnapshot->atomData);
1147     DataStorage* cgConfig = &(curSnapshot->cgData);
1148 gezelter 1581
1149 mciznick 1598 //calculate the center of mass of cutoff group
1150 gezelter 1545
1151 mciznick 1598 SimInfo::MoleculeIterator mi;
1152     Molecule* mol;
1153     Molecule::CutoffGroupIterator ci;
1154     CutoffGroup* cg;
1155 gezelter 1581
1156 mciznick 1598 if (info_->getNCutoffGroups() > 0)
1157     {
1158     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
1159     {
1160     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
1161     {
1162     cerr << "branch1\n";
1163     cerr << "globind = " << cg->getGlobalIndex() << "\n";
1164     cg->updateCOM();
1165     }
1166     }
1167     } else
1168     {
1169     // center of mass of the group is the same as position of the atom
1170     // if cutoff group does not exist
1171     cerr << "branch2\n";
1172     cgConfig->position = config->position;
1173     }
1174 gezelter 1581
1175 mciznick 1598 fDecomp_->zeroWorkArrays();
1176     fDecomp_->distributeData();
1177 gezelter 1581
1178 mciznick 1598 int cg1, cg2, atom1, atom2, topoDist;
1179     Vector3d d_grp, dag, d;
1180     RealType rgrpsq, rgrp, r2, r;
1181     RealType electroMult, vdwMult;
1182     RealType vij;
1183     Vector3d fij, fg, f1;
1184     tuple3<RealType, RealType, RealType> cuts;
1185     RealType rCutSq;
1186     bool in_switching_region;
1187     RealType sw, dswdr, swderiv;
1188     vector<int> atomListColumn, atomListRow, atomListLocal;
1189     InteractionData idat;
1190     SelfData sdat;
1191     RealType mf;
1192     RealType lrPot;
1193     RealType vpair;
1194     potVec longRangePotential(0.0);
1195     potVec workPot(0.0);
1196 gezelter 1544
1197 mciznick 1598 int loopStart, loopEnd;
1198 gezelter 1544
1199 mciznick 1598 idat.vdwMult = &vdwMult;
1200     idat.electroMult = &electroMult;
1201     idat.pot = &workPot;
1202     sdat.pot = fDecomp_->getEmbeddingPotential();
1203     idat.vpair = &vpair;
1204     idat.f1 = &f1;
1205     idat.sw = &sw;
1206     idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
1207     idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
1208 chuckv 1595
1209 mciznick 1598 loopEnd = PAIR_LOOP;
1210     if (info_->requiresPrepair())
1211     {
1212     loopStart = PREPAIR_LOOP;
1213     } else
1214     {
1215     loopStart = PAIR_LOOP;
1216     }
1217 gezelter 1545
1218 mciznick 1598 for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++)
1219     {
1220 gezelter 1545
1221 mciznick 1598 if (iLoop == loopStart)
1222     {
1223     bool update_nlist = fDecomp_->checkNeighborList();
1224     if (update_nlist)
1225     neighborList = fDecomp_->buildNeighborList();
1226 gezelter 1576
1227 mciznick 1598 }
1228 gezelter 1545
1229 mciznick 1598 for (vector<pair<int, int> >::iterator it = neighborList.begin(); it != neighborList.end(); ++it)
1230     {
1231     cg1 = (*it).first;
1232     cg2 = (*it).second;
1233 gezelter 1593
1234 mciznick 1598 cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
1235 gezelter 1575
1236 mciznick 1598 d_grp = fDecomp_->getIntergroupVector(cg1, cg2);
1237     curSnapshot->wrapVector(d_grp);
1238     rgrpsq = d_grp.lengthSquare();
1239 gezelter 1546
1240 mciznick 1598 rCutSq = cuts.second;
1241 gezelter 1593
1242 mciznick 1598 if (rgrpsq < rCutSq)
1243     {
1244     idat.rcut = &cuts.first;
1245     if (iLoop == PAIR_LOOP)
1246     {
1247     vij = 0.0;
1248     fij = V3Zero;
1249     }
1250 gezelter 1545
1251 mciznick 1598 in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, rgrp);
1252 gezelter 1545
1253 mciznick 1598 atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
1254     atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
1255 gezelter 1545
1256 mciznick 1598 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
1257     {
1258     atom1 = (*ia);
1259 gezelter 1545
1260 mciznick 1598 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
1261     {
1262     atom2 = (*jb);
1263 gezelter 1545
1264 mciznick 1598 if (!fDecomp_->skipAtomPair(atom1, atom2))
1265     {
1266     vpair = 0.0;
1267     workPot = 0.0;
1268     f1 = V3Zero;
1269 gezelter 1590
1270 mciznick 1598 fDecomp_->fillInteractionData(idat, atom1, atom2);
1271 gezelter 1570
1272 mciznick 1598 topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
1273     vdwMult = vdwScale_[topoDist];
1274     electroMult = electrostaticScale_[topoDist];
1275 gezelter 1590
1276 mciznick 1598 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
1277     {
1278     idat.d = &d_grp;
1279     idat.r2 = &rgrpsq;
1280     cerr << "dgrp = " << d_grp << "\n";
1281     } else
1282     {
1283     d = fDecomp_->getInteratomicVector(atom1, atom2);
1284     curSnapshot->wrapVector(d);
1285     r2 = d.lengthSquare();
1286     cerr << "datm = " << d << "\n";
1287     idat.d = &d;
1288     idat.r2 = &r2;
1289     }
1290 gezelter 1590
1291 mciznick 1598 cerr << "idat.d = " << *(idat.d) << "\n";
1292     r = sqrt(*(idat.r2));
1293     idat.rij = &r;
1294 gezelter 1545
1295 mciznick 1598 if (iLoop == PREPAIR_LOOP)
1296     {
1297     interactionMan_->doPrePair(idat);
1298     } else
1299     {
1300     interactionMan_->doPair(idat);
1301     fDecomp_->unpackInteractionData(idat, atom1, atom2);
1302 gezelter 1545
1303 mciznick 1598 cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << "\n";
1304     vij += vpair;
1305     fij += f1;
1306     tau -= outProduct(*(idat.d), f1);
1307     }
1308     }
1309     }
1310     }
1311 gezelter 1570
1312 mciznick 1598 if (iLoop == PAIR_LOOP)
1313     {
1314     if (in_switching_region)
1315     {
1316     swderiv = vij * dswdr / rgrp;
1317     fg = swderiv * d_grp;
1318     fij += fg;
1319 gezelter 1570
1320 mciznick 1598 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
1321     {
1322     tau -= outProduct(*(idat.d), fg);
1323     }
1324 gezelter 1583
1325 mciznick 1598 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
1326     {
1327     atom1 = (*ia);
1328     mf = fDecomp_->getMassFactorRow(atom1);
1329     // fg is the force on atom ia due to cutoff group's
1330     // presence in switching region
1331     fg = swderiv * d_grp * mf;
1332     fDecomp_->addForceToAtomRow(atom1, fg);
1333 gezelter 1575
1334 mciznick 1598 if (atomListRow.size() > 1)
1335     {
1336     if (info_->usesAtomicVirial())
1337     {
1338     // find the distance between the atom
1339     // and the center of the cutoff group:
1340     dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
1341     tau -= outProduct(dag, fg);
1342     }
1343     }
1344     }
1345     for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
1346     {
1347     atom2 = (*jb);
1348     mf = fDecomp_->getMassFactorColumn(atom2);
1349     // fg is the force on atom jb due to cutoff group's
1350     // presence in switching region
1351     fg = -swderiv * d_grp * mf;
1352     fDecomp_->addForceToAtomColumn(atom2, fg);
1353 gezelter 246
1354 mciznick 1598 if (atomListColumn.size() > 1)
1355     {
1356     if (info_->usesAtomicVirial())
1357     {
1358     // find the distance between the atom
1359     // and the center of the cutoff group:
1360     dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
1361     tau -= outProduct(dag, fg);
1362     }
1363     }
1364     }
1365     }
1366     //if (!SIM_uses_AtomicVirial) {
1367     // tau -= outProduct(d_grp, fij);
1368     //}
1369     }
1370     }
1371     }
1372    
1373     if (iLoop == PREPAIR_LOOP)
1374     {
1375     if (info_->requiresPrepair())
1376     {
1377    
1378     fDecomp_->collectIntermediateData();
1379    
1380     for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1381     {
1382     fDecomp_->fillSelfData(sdat, atom1);
1383     interactionMan_->doPreForce(sdat);
1384     }
1385    
1386     fDecomp_->distributeIntermediateData();
1387    
1388     }
1389     }
1390    
1391     }
1392    
1393     fDecomp_->collectData();
1394    
1395     if (info_->requiresSelfCorrection())
1396     {
1397    
1398     for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1399     {
1400     fDecomp_->fillSelfData(sdat, atom1);
1401     interactionMan_->doSelfCorrection(sdat);
1402     }
1403    
1404     }
1405    
1406     longRangePotential = *(fDecomp_->getEmbeddingPotential()) + *(fDecomp_->getPairwisePotential());
1407    
1408     lrPot = longRangePotential.sum();
1409    
1410     //store the tau and long range potential
1411     curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
1412     curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
1413     curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
1414     }
1415    
1416     void ForceManager::postCalculation() {
1417     SimInfo::MoleculeIterator mi;
1418     Molecule* mol;
1419     Molecule::RigidBodyIterator rbIter;
1420     RigidBody* rb;
1421     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1422    
1423     // collect the atomic forces onto rigid bodies
1424    
1425     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
1426     {
1427     for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
1428     {
1429     Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
1430     tau += rbTau;
1431     }
1432     }
1433    
1434 gezelter 1126 #ifdef IS_MPI
1435 mciznick 1598 Mat3x3d tmpTau(tau);
1436     MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
1437     9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1438 gezelter 1126 #endif
1439 mciznick 1598 curSnapshot->statData.setTau(tau);
1440     }
1441 gezelter 246
1442 gezelter 1390 } //end namespace OpenMD

Properties

Name Value
svn:executable *
svn:keywords Author Id Revision Date