ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/devel_omp/src/brains/ForceManager.cpp
Revision: 1598
Committed: Wed Jul 27 14:26:53 2011 UTC (13 years, 9 months ago) by mciznick
File size: 33040 byte(s)
Log Message: 
Updated ordered neighbor list generation.

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     longRangeInteractionsRapaport();
370     postCalculation();
371     }
372 chuckv 1595
373 mciznick 1598 void ForceManager::preCalculation() {
374     SimInfo::MoleculeIterator mi;
375     Molecule* mol;
376     Molecule::AtomIterator ai;
377     Atom* atom;
378     Molecule::RigidBodyIterator rbIter;
379     RigidBody* rb;
380     Molecule::CutoffGroupIterator ci;
381     CutoffGroup* cg;
382 chuckv 1595
383 mciznick 1598 // forces are zeroed here, before any are accumulated.
384 chuckv 1595
385 mciznick 1598 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
386     {
387     for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai))
388     {
389     atom->zeroForcesAndTorques();
390     }
391 chuckv 1595
392 mciznick 1598 //change the positions of atoms which belong to the rigidbodies
393     for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
394     {
395     rb->zeroForcesAndTorques();
396     }
397 chuckv 1595
398 mciznick 1598 if (info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms())
399     {
400     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
401     {
402     //calculate the center of mass of cutoff group
403     cg->updateCOM();
404     }
405     }
406     }
407 chuckv 1595
408 mciznick 1598 // Zero out the stress tensor
409     tau *= 0.0;
410 chuckv 1595
411 mciznick 1598 }
412 chuckv 1595
413 mciznick 1598 void ForceManager::shortRangeInteractions() {
414     Molecule* mol;
415     RigidBody* rb;
416     Bond* bond;
417     Bend* bend;
418     Torsion* torsion;
419     Inversion* inversion;
420     SimInfo::MoleculeIterator mi;
421     Molecule::RigidBodyIterator rbIter;
422     Molecule::BondIterator bondIter;
423     ;
424     Molecule::BendIterator bendIter;
425     Molecule::TorsionIterator torsionIter;
426     Molecule::InversionIterator inversionIter;
427     RealType bondPotential = 0.0;
428     RealType bendPotential = 0.0;
429     RealType torsionPotential = 0.0;
430     RealType inversionPotential = 0.0;
431 chuckv 1595
432 mciznick 1598 //calculate short range interactions
433     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
434     {
435 chuckv 1595
436 mciznick 1598 //change the positions of atoms which belong to the rigidbodies
437     for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
438     {
439     rb->updateAtoms();
440     }
441 chuckv 1595
442 mciznick 1598 for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter))
443     {
444     bond->calcForce();
445     bondPotential += bond->getPotential();
446     }
447 chuckv 1595
448 mciznick 1598 for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter))
449     {
450 chuckv 1595
451 mciznick 1598 RealType angle;
452     bend->calcForce(angle);
453     RealType currBendPot = bend->getPotential();
454 chuckv 1595
455 mciznick 1598 bendPotential += bend->getPotential();
456     map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
457     if (i == bendDataSets.end())
458     {
459     BendDataSet dataSet;
460     dataSet.prev.angle = dataSet.curr.angle = angle;
461     dataSet.prev.potential = dataSet.curr.potential = currBendPot;
462     dataSet.deltaV = 0.0;
463     bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
464     } else
465     {
466     i->second.prev.angle = i->second.curr.angle;
467     i->second.prev.potential = i->second.curr.potential;
468     i->second.curr.angle = angle;
469     i->second.curr.potential = currBendPot;
470     i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
471     }
472     }
473 chuckv 1595
474 mciznick 1598 for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter))
475     {
476     RealType angle;
477     torsion->calcForce(angle);
478     RealType currTorsionPot = torsion->getPotential();
479     torsionPotential += torsion->getPotential();
480     map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
481     if (i == torsionDataSets.end())
482     {
483     TorsionDataSet dataSet;
484     dataSet.prev.angle = dataSet.curr.angle = angle;
485     dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
486     dataSet.deltaV = 0.0;
487     torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
488     } else
489     {
490     i->second.prev.angle = i->second.curr.angle;
491     i->second.prev.potential = i->second.curr.potential;
492     i->second.curr.angle = angle;
493     i->second.curr.potential = currTorsionPot;
494     i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
495     }
496     }
497 chuckv 1595
498 mciznick 1598 for (inversion = mol->beginInversion(inversionIter); inversion != NULL; inversion = mol->nextInversion(
499     inversionIter))
500     {
501     RealType angle;
502     inversion->calcForce(angle);
503     RealType currInversionPot = inversion->getPotential();
504     inversionPotential += inversion->getPotential();
505     map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
506     if (i == inversionDataSets.end())
507     {
508     InversionDataSet dataSet;
509     dataSet.prev.angle = dataSet.curr.angle = angle;
510     dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
511     dataSet.deltaV = 0.0;
512     inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
513     } else
514     {
515     i->second.prev.angle = i->second.curr.angle;
516     i->second.prev.potential = i->second.curr.potential;
517     i->second.curr.angle = angle;
518     i->second.curr.potential = currInversionPot;
519     i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential);
520     }
521     }
522     }
523 chuckv 1595
524 mciznick 1598 RealType shortRangePotential = bondPotential + bendPotential + torsionPotential + inversionPotential;
525     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
526     curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
527     curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
528     curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
529     curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
530     curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
531     }
532 chuckv 1595
533 mciznick 1598 void ForceManager::longRangeInteractionsRapaport() {
534 chuckv 1595
535 mciznick 1598 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
536     DataStorage* config = &(curSnapshot->atomData);
537     DataStorage* cgConfig = &(curSnapshot->cgData);
538 chuckv 1595
539 mciznick 1598 //calculate the center of mass of cutoff group
540 chuckv 1595
541 mciznick 1598 SimInfo::MoleculeIterator mi;
542     Molecule* mol;
543     Molecule::CutoffGroupIterator ci;
544     CutoffGroup* cg;
545 chuckv 1595
546 mciznick 1598 if (info_->getNCutoffGroups() > 0)
547     {
548     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
549     {
550     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
551     {
552     // cerr << "branch1\n";
553     // cerr << "globind = " << cg->getGlobalIndex() << ":" << __LINE__ << "\n";
554     cg->updateCOM();
555    
556     // cerr << "gbI: " << cg->getGlobalIndex() << " locI: " << cg->getLocalIndex() << " x: "
557     // << cgConfig->position[cg->getLocalIndex()].x() << " y: " << cgConfig->position[cg->getLocalIndex()].y()
558     // << " z: " << cgConfig->position[cg->getLocalIndex()].z() << "\n";
559     }
560     }
561     } else
562     {
563     // center of mass of the group is the same as position of the atom
564     // if cutoff group does not exist
565     // cerr << ":" << __LINE__ << "branch2\n";
566     cgConfig->position = config->position;
567     }
568    
569     fDecomp_->zeroWorkArrays();
570     fDecomp_->distributeData();
571    
572     int atom1, atom2, topoDist;
573     CutoffGroup *cg1;
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     InteractionData idat;
585     SelfData sdat;
586     RealType mf;
587     RealType lrPot;
588     RealType vpair;
589     potVec longRangePotential(0.0);
590     potVec workPot(0.0);
591    
592     int loopStart, loopEnd;
593    
594     idat.vdwMult = &vdwMult;
595     idat.electroMult = &electroMult;
596     idat.pot = &workPot;
597     sdat.pot = fDecomp_->getEmbeddingPotential();
598     idat.vpair = &vpair;
599     idat.f1 = &f1;
600     idat.sw = &sw;
601     idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
602     idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
603    
604     loopEnd = PAIR_LOOP;
605     if (info_->requiresPrepair())
606     {
607     loopStart = PREPAIR_LOOP;
608     } else
609     {
610     loopStart = PAIR_LOOP;
611     }
612    
613     for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++)
614     {
615    
616     if (iLoop == loopStart)
617     {
618     bool update_nlist = fDecomp_->checkNeighborList();
619     if (update_nlist)
620     neighborMatW = fDecomp_->buildLayerBasedNeighborList();
621     }
622    
623     // printf("before omp loop\n");
624     //#pragma omp parallel for num_threads(3) default(none) shared(curSnapshot, idat, iLoop, sw, cerr) \
625     private(i, j, cg1, cg2, cuts, d_grp, rgrpsq, rCutSq, vij, fij, in_switching_region, rgrp, dswdr, atomListRow, atomListColumn, atom1, atom2, vpair, workPot, f1, topoDist, vdwMult, electroMult, d, r2, r, swderiv, fg, mf, dag)
626     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
627     {
628     for (cg1 = mol->beginCutoffGroup(ci); cg1 != NULL; cg1 = mol->nextCutoffGroup(ci))
629     {
630     // printf("Thread %d executes loop iteration %d\n", omp_get_thread_num(), i);
631     for (vector<CutoffGroup *>::iterator cg2 = neighborMatW[cg1->getGlobalIndex()].begin(); cg2 != neighborMatW[cg1->getGlobalIndex()].end(); ++cg2)
632     {
633    
634     cuts = fDecomp_->getGroupCutoffs(cg1->getGlobalIndex(), (*cg2)->getGlobalIndex());
635    
636     d_grp = fDecomp_->getIntergroupVector(cg1, (*cg2));
637     curSnapshot->wrapVector(d_grp);
638     rgrpsq = d_grp.lengthSquare();
639    
640     rCutSq = cuts.second;
641    
642     if (rgrpsq < rCutSq)
643     {
644     idat.rcut = &cuts.first;
645     if (iLoop == PAIR_LOOP)
646     {
647     vij = 0.0;
648     fij = V3Zero;
649     }
650    
651     in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, rgrp);
652    
653     atomListRow = fDecomp_->getAtomsInGroupRow(cg1->getGlobalIndex());
654     atomListColumn = fDecomp_->getAtomsInGroupColumn((*cg2)->getGlobalIndex());
655    
656     for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
657     {
658     atom1 = (*ia);
659    
660     for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
661     {
662     atom2 = (*jb);
663    
664     if (!fDecomp_->skipAtomPair(atom1, atom2))
665     {
666     vpair = 0.0;
667     workPot = 0.0;
668     f1 = V3Zero;
669    
670     fDecomp_->fillInteractionData(idat, atom1, atom2);
671    
672     topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
673     vdwMult = vdwScale_[topoDist];
674     electroMult = electrostaticScale_[topoDist];
675    
676     if (atomListRow.size() == 1 && atomListColumn.size() == 1)
677     {
678     idat.d = &d_grp;
679     idat.r2 = &rgrpsq;
680     // cerr << "dgrp = " << d_grp << ":" << __LINE__ << "\n";
681     } else
682     {
683     d = fDecomp_->getInteratomicVector(atom1, atom2);
684     curSnapshot->wrapVector(d);
685     r2 = d.lengthSquare();
686     // cerr << "datm = " << d << ":" << __LINE__ << "\n";
687     idat.d = &d;
688     idat.r2 = &r2;
689     }
690    
691     // cerr << "idat.d = " << *(idat.d) << ":" << __LINE__ << "\n";
692     r = sqrt(*(idat.r2));
693     idat.rij = &r;
694     // cerr << "idat.rij = " << *(idat.rij) << "\n";
695    
696     if (iLoop == PREPAIR_LOOP)
697     {
698     interactionMan_->doPrePair(idat);
699     } else
700     {
701     interactionMan_->doPair(idat);
702     fDecomp_->unpackInteractionData(idat, atom1, atom2);
703    
704     // cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << ":" << __LINE__ << "\n";
705     vij += vpair;
706     fij += f1;
707     tau -= outProduct(*(idat.d), f1);
708     }
709     }
710     }
711     }
712    
713     if (iLoop == PAIR_LOOP)
714     {
715     if (in_switching_region)
716     {
717     swderiv = vij * dswdr / rgrp;
718     fg = swderiv * d_grp;
719     fij += fg;
720    
721     if (atomListRow.size() == 1 && atomListColumn.size() == 1)
722     {
723     tau -= outProduct(*(idat.d), fg);
724     }
725    
726     for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
727     {
728     atom1 = (*ia);
729     mf = fDecomp_->getMassFactorRow(atom1);
730     // fg is the force on atom ia due to cutoff group's
731     // presence in switching region
732     fg = swderiv * d_grp * mf;
733     fDecomp_->addForceToAtomRow(atom1, fg);
734    
735     if (atomListRow.size() > 1)
736     {
737     if (info_->usesAtomicVirial())
738     {
739     // find the distance between the atom
740     // and the center of the cutoff group:
741     dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1->getGlobalIndex());
742     tau -= outProduct(dag, fg);
743     }
744     }
745     }
746     for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
747     {
748     atom2 = (*jb);
749     mf = fDecomp_->getMassFactorColumn(atom2);
750     // fg is the force on atom jb due to cutoff group's
751     // presence in switching region
752     fg = -swderiv * d_grp * mf;
753     fDecomp_->addForceToAtomColumn(atom2, fg);
754    
755     if (atomListColumn.size() > 1)
756     {
757     if (info_->usesAtomicVirial())
758     {
759     // find the distance between the atom
760     // and the center of the cutoff group:
761     dag = fDecomp_->getAtomToGroupVectorColumn(atom2, (*cg2)->getGlobalIndex());
762     tau -= outProduct(dag, fg);
763     }
764     }
765     }
766     }
767     //if (!SIM_uses_AtomicVirial) {
768     // tau -= outProduct(d_grp, fij);
769     //}
770     }
771 chuckv 1595 }
772     }
773     }
774 mciznick 1598 }// END: omp for loop
775     // printf("after omp loop\n");
776 chuckv 1595
777 mciznick 1598 if (iLoop == PREPAIR_LOOP)
778     {
779     if (info_->requiresPrepair())
780     {
781 chuckv 1595
782 mciznick 1598 fDecomp_->collectIntermediateData();
783 chuckv 1595
784 mciznick 1598 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
785     {
786     fDecomp_->fillSelfData(sdat, atom1);
787     interactionMan_->doPreForce(sdat);
788     }
789 chuckv 1595
790 mciznick 1598 fDecomp_->distributeIntermediateData();
791 chuckv 1595
792 mciznick 1598 }
793     }
794     }
795 chuckv 1595
796 mciznick 1598 fDecomp_->collectData();
797 chuckv 1595
798 mciznick 1598 if (info_->requiresSelfCorrection())
799     {
800 chuckv 1595
801 mciznick 1598 for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
802     {
803     fDecomp_->fillSelfData(sdat, atom1);
804     interactionMan_->doSelfCorrection(sdat);
805     }
806 chuckv 1595
807 mciznick 1598 }
808 chuckv 1595
809 mciznick 1598 longRangePotential = *(fDecomp_->getEmbeddingPotential()) + *(fDecomp_->getPairwisePotential());
810 chuckv 1595
811 mciznick 1598 lrPot = longRangePotential.sum();
812 chuckv 1595
813 mciznick 1598 //store the tau and long range potential
814     curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
815     curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
816     curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
817     }
818 chuckv 1595
819 mciznick 1598 void ForceManager::longRangeInteractions() {
820 chuckv 1595
821 mciznick 1598 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
822     DataStorage* config = &(curSnapshot->atomData);
823     DataStorage* cgConfig = &(curSnapshot->cgData);
824 gezelter 1581
825 mciznick 1598 //calculate the center of mass of cutoff group
826 gezelter 1545
827 mciznick 1598 SimInfo::MoleculeIterator mi;
828     Molecule* mol;
829     Molecule::CutoffGroupIterator ci;
830     CutoffGroup* cg;
831 gezelter 1581
832 mciznick 1598 if (info_->getNCutoffGroups() > 0)
833     {
834     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
835     {
836     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
837     {
838     cerr << "branch1\n";
839     cerr << "globind = " << cg->getGlobalIndex() << "\n";
840     cg->updateCOM();
841     }
842     }
843     } else
844     {
845     // center of mass of the group is the same as position of the atom
846     // if cutoff group does not exist
847     cerr << "branch2\n";
848     cgConfig->position = config->position;
849     }
850 gezelter 1581
851 mciznick 1598 fDecomp_->zeroWorkArrays();
852     fDecomp_->distributeData();
853 gezelter 1581
854 mciznick 1598 int cg1, cg2, atom1, atom2, topoDist;
855     Vector3d d_grp, dag, d;
856     RealType rgrpsq, rgrp, r2, r;
857     RealType electroMult, vdwMult;
858     RealType vij;
859     Vector3d fij, fg, f1;
860     tuple3<RealType, RealType, RealType> cuts;
861     RealType rCutSq;
862     bool in_switching_region;
863     RealType sw, dswdr, swderiv;
864     vector<int> atomListColumn, atomListRow, atomListLocal;
865     InteractionData idat;
866     SelfData sdat;
867     RealType mf;
868     RealType lrPot;
869     RealType vpair;
870     potVec longRangePotential(0.0);
871     potVec workPot(0.0);
872 gezelter 1544
873 mciznick 1598 int loopStart, loopEnd;
874 gezelter 1544
875 mciznick 1598 idat.vdwMult = &vdwMult;
876     idat.electroMult = &electroMult;
877     idat.pot = &workPot;
878     sdat.pot = fDecomp_->getEmbeddingPotential();
879     idat.vpair = &vpair;
880     idat.f1 = &f1;
881     idat.sw = &sw;
882     idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
883     idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
884 chuckv 1595
885 mciznick 1598 loopEnd = PAIR_LOOP;
886     if (info_->requiresPrepair())
887     {
888     loopStart = PREPAIR_LOOP;
889     } else
890     {
891     loopStart = PAIR_LOOP;
892     }
893 gezelter 1545
894 mciznick 1598 for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++)
895     {
896 gezelter 1545
897 mciznick 1598 if (iLoop == loopStart)
898     {
899     bool update_nlist = fDecomp_->checkNeighborList();
900     if (update_nlist)
901     neighborList = fDecomp_->buildNeighborList();
902 gezelter 1576
903 mciznick 1598 }
904 gezelter 1545
905 mciznick 1598 for (vector<pair<int, int> >::iterator it = neighborList.begin(); it != neighborList.end(); ++it)
906     {
907     cg1 = (*it).first;
908     cg2 = (*it).second;
909 gezelter 1593
910 mciznick 1598 cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
911 gezelter 1575
912 mciznick 1598 d_grp = fDecomp_->getIntergroupVector(cg1, cg2);
913     curSnapshot->wrapVector(d_grp);
914     rgrpsq = d_grp.lengthSquare();
915 gezelter 1546
916 mciznick 1598 rCutSq = cuts.second;
917 gezelter 1593
918 mciznick 1598 if (rgrpsq < rCutSq)
919     {
920     idat.rcut = &cuts.first;
921     if (iLoop == PAIR_LOOP)
922     {
923     vij = 0.0;
924     fij = V3Zero;
925     }
926 gezelter 1545
927 mciznick 1598 in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, rgrp);
928 gezelter 1545
929 mciznick 1598 atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
930     atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
931 gezelter 1545
932 mciznick 1598 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
933     {
934     atom1 = (*ia);
935 gezelter 1545
936 mciznick 1598 for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
937     {
938     atom2 = (*jb);
939 gezelter 1545
940 mciznick 1598 if (!fDecomp_->skipAtomPair(atom1, atom2))
941     {
942     vpair = 0.0;
943     workPot = 0.0;
944     f1 = V3Zero;
945 gezelter 1590
946 mciznick 1598 fDecomp_->fillInteractionData(idat, atom1, atom2);
947 gezelter 1570
948 mciznick 1598 topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
949     vdwMult = vdwScale_[topoDist];
950     electroMult = electrostaticScale_[topoDist];
951 gezelter 1590
952 mciznick 1598 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
953     {
954     idat.d = &d_grp;
955     idat.r2 = &rgrpsq;
956     cerr << "dgrp = " << d_grp << "\n";
957     } else
958     {
959     d = fDecomp_->getInteratomicVector(atom1, atom2);
960     curSnapshot->wrapVector(d);
961     r2 = d.lengthSquare();
962     cerr << "datm = " << d << "\n";
963     idat.d = &d;
964     idat.r2 = &r2;
965     }
966 gezelter 1590
967 mciznick 1598 cerr << "idat.d = " << *(idat.d) << "\n";
968     r = sqrt(*(idat.r2));
969     idat.rij = &r;
970 gezelter 1545
971 mciznick 1598 if (iLoop == PREPAIR_LOOP)
972     {
973     interactionMan_->doPrePair(idat);
974     } else
975     {
976     interactionMan_->doPair(idat);
977     fDecomp_->unpackInteractionData(idat, atom1, atom2);
978 gezelter 1545
979 mciznick 1598 cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << "\n";
980     vij += vpair;
981     fij += f1;
982     tau -= outProduct(*(idat.d), f1);
983     }
984     }
985     }
986     }
987 gezelter 1570
988 mciznick 1598 if (iLoop == PAIR_LOOP)
989     {
990     if (in_switching_region)
991     {
992     swderiv = vij * dswdr / rgrp;
993     fg = swderiv * d_grp;
994     fij += fg;
995 gezelter 1570
996 mciznick 1598 if (atomListRow.size() == 1 && atomListColumn.size() == 1)
997     {
998     tau -= outProduct(*(idat.d), fg);
999     }
1000 gezelter 1583
1001 mciznick 1598 for (vector<int>::iterator ia = atomListRow.begin(); ia != atomListRow.end(); ++ia)
1002     {
1003     atom1 = (*ia);
1004     mf = fDecomp_->getMassFactorRow(atom1);
1005     // fg is the force on atom ia due to cutoff group's
1006     // presence in switching region
1007     fg = swderiv * d_grp * mf;
1008     fDecomp_->addForceToAtomRow(atom1, fg);
1009 gezelter 1575
1010 mciznick 1598 if (atomListRow.size() > 1)
1011     {
1012     if (info_->usesAtomicVirial())
1013     {
1014     // find the distance between the atom
1015     // and the center of the cutoff group:
1016     dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
1017     tau -= outProduct(dag, fg);
1018     }
1019     }
1020     }
1021     for (vector<int>::iterator jb = atomListColumn.begin(); jb != atomListColumn.end(); ++jb)
1022     {
1023     atom2 = (*jb);
1024     mf = fDecomp_->getMassFactorColumn(atom2);
1025     // fg is the force on atom jb due to cutoff group's
1026     // presence in switching region
1027     fg = -swderiv * d_grp * mf;
1028     fDecomp_->addForceToAtomColumn(atom2, fg);
1029 gezelter 246
1030 mciznick 1598 if (atomListColumn.size() > 1)
1031     {
1032     if (info_->usesAtomicVirial())
1033     {
1034     // find the distance between the atom
1035     // and the center of the cutoff group:
1036     dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
1037     tau -= outProduct(dag, fg);
1038     }
1039     }
1040     }
1041     }
1042     //if (!SIM_uses_AtomicVirial) {
1043     // tau -= outProduct(d_grp, fij);
1044     //}
1045     }
1046     }
1047     }
1048    
1049     if (iLoop == PREPAIR_LOOP)
1050     {
1051     if (info_->requiresPrepair())
1052     {
1053    
1054     fDecomp_->collectIntermediateData();
1055    
1056     for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1057     {
1058     fDecomp_->fillSelfData(sdat, atom1);
1059     interactionMan_->doPreForce(sdat);
1060     }
1061    
1062     fDecomp_->distributeIntermediateData();
1063    
1064     }
1065     }
1066    
1067     }
1068    
1069     fDecomp_->collectData();
1070    
1071     if (info_->requiresSelfCorrection())
1072     {
1073    
1074     for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++)
1075     {
1076     fDecomp_->fillSelfData(sdat, atom1);
1077     interactionMan_->doSelfCorrection(sdat);
1078     }
1079    
1080     }
1081    
1082     longRangePotential = *(fDecomp_->getEmbeddingPotential()) + *(fDecomp_->getPairwisePotential());
1083    
1084     lrPot = longRangePotential.sum();
1085    
1086     //store the tau and long range potential
1087     curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
1088     curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
1089     curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
1090     }
1091    
1092     void ForceManager::postCalculation() {
1093     SimInfo::MoleculeIterator mi;
1094     Molecule* mol;
1095     Molecule::RigidBodyIterator rbIter;
1096     RigidBody* rb;
1097     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1098    
1099     // collect the atomic forces onto rigid bodies
1100    
1101     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
1102     {
1103     for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter))
1104     {
1105     Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
1106     tau += rbTau;
1107     }
1108     }
1109    
1110 gezelter 1126 #ifdef IS_MPI
1111 mciznick 1598 Mat3x3d tmpTau(tau);
1112     MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
1113     9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1114 gezelter 1126 #endif
1115 mciznick 1598 curSnapshot->statData.setTau(tau);
1116     }
1117 gezelter 246
1118 gezelter 1390 } //end namespace OpenMD

Properties

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