ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/devel_omp/src/brains/ForceManager.cpp
Revision: 1614
Committed: Tue Aug 23 20:55:51 2011 UTC (13 years, 8 months ago) by mciznick
File size: 44345 byte(s)
Log Message: 
Updated scalability of OpenMP threads.

File Contents

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

Properties

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