src/brains/ForceManager.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
 * [4]  Vardeman & Gezelter, in progress (2009).
 */

/**
 * @file ForceManager.cpp
 * @author tlin
 * @date 11/09/2004
 * @time 10:39am
 * @version 1.0
 */

#include "brains/ForceManager.hpp"
#include "primitives/Molecule.hpp"
#include "UseTheForce/doForces_interface.h"
#define __OPENMD_C
#include "UseTheForce/DarkSide/fInteractionMap.h"
#include "utils/simError.h"
#include "primitives/Bond.hpp"
#include "primitives/Bend.hpp"
#include "primitives/Torsion.hpp"
#include "primitives/Inversion.hpp"
namespace OpenMD {

  void ForceManager::calcForces() {

    if (!info_->isFortranInitialized()) {
      info_->update();
    }

    preCalculation();

    calcShortRangeInteraction();

    calcLongRangeInteraction();

    postCalculation();

  }

  void ForceManager::preCalculation() {
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    Molecule::AtomIterator ai;
    Atom* atom;
    Molecule::RigidBodyIterator rbIter;
    RigidBody* rb;

    // forces are zeroed here, before any are accumulated.
    // NOTE: do not rezero the forces in Fortran.

    for (mol = info_->beginMolecule(mi); mol != NULL;
         mol = info_->nextMolecule(mi)) {
      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
  atom->zeroForcesAndTorques();
      }

      //change the positions of atoms which belong to the rigidbodies
      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
           rb = mol->nextRigidBody(rbIter)) {
  rb->zeroForcesAndTorques();
      }

    }

    // Zero out the stress tensor
    tau *= 0.0;

  }

  void ForceManager::calcShortRangeInteraction() {
    Molecule* mol;
    RigidBody* rb;
    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    Inversion* inversion;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Molecule::BondIterator bondIter;;
    Molecule::BendIterator  bendIter;
    Molecule::TorsionIterator  torsionIter;
    Molecule::InversionIterator  inversionIter;
    RealType bondPotential = 0.0;
    RealType bendPotential = 0.0;
    RealType torsionPotential = 0.0;
    RealType inversionPotential = 0.0;

    //calculate short range interactions
    for (mol = info_->beginMolecule(mi); mol != NULL;
         mol = info_->nextMolecule(mi)) {

      //change the positions of atoms which belong to the rigidbodies
      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
           rb = mol->nextRigidBody(rbIter)) {
        rb->updateAtoms();
        rb->updateAtomVel();
      }

      for (bond = mol->beginBond(bondIter); bond != NULL;
           bond = mol->nextBond(bondIter)) {
        bond->calcForce();
        bondPotential += bond->getPotential();
      }

      for (bend = mol->beginBend(bendIter); bend != NULL;
           bend = mol->nextBend(bendIter)) {

        RealType angle;
        bend->calcForce(angle);
        RealType currBendPot = bend->getPotential();

        bendPotential += bend->getPotential();
        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
        if (i == bendDataSets.end()) {
          BendDataSet dataSet;
          dataSet.prev.angle = dataSet.curr.angle = angle;
          dataSet.prev.potential = dataSet.curr.potential = currBendPot;
          dataSet.deltaV = 0.0;
          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
        }else {
          i->second.prev.angle = i->second.curr.angle;
          i->second.prev.potential = i->second.curr.potential;
          i->second.curr.angle = angle;
          i->second.curr.potential = currBendPot;
          i->second.deltaV =  fabs(i->second.curr.potential -
                                   i->second.prev.potential);
        }
      }

      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
           torsion = mol->nextTorsion(torsionIter)) {
        RealType angle;
        torsion->calcForce(angle);
        RealType currTorsionPot = torsion->getPotential();
        torsionPotential += torsion->getPotential();
        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
        if (i == torsionDataSets.end()) {
          TorsionDataSet dataSet;
          dataSet.prev.angle = dataSet.curr.angle = angle;
          dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
          dataSet.deltaV = 0.0;
          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
        }else {
          i->second.prev.angle = i->second.curr.angle;
          i->second.prev.potential = i->second.curr.potential;
          i->second.curr.angle = angle;
          i->second.curr.potential = currTorsionPot;
          i->second.deltaV =  fabs(i->second.curr.potential -
                                   i->second.prev.potential);
        }
      }

      for (inversion = mol->beginInversion(inversionIter);
     inversion != NULL;
           inversion = mol->nextInversion(inversionIter)) {
        RealType angle;
        inversion->calcForce(angle);
        RealType currInversionPot = inversion->getPotential();
        inversionPotential += inversion->getPotential();
        std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
        if (i == inversionDataSets.end()) {
          InversionDataSet dataSet;
          dataSet.prev.angle = dataSet.curr.angle = angle;
          dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
          dataSet.deltaV = 0.0;
          inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
        }else {
          i->second.prev.angle = i->second.curr.angle;
          i->second.prev.potential = i->second.curr.potential;
          i->second.curr.angle = angle;
          i->second.curr.potential = currInversionPot;
          i->second.deltaV =  fabs(i->second.curr.potential -
                                   i->second.prev.potential);
        }
      }
    }

    RealType  shortRangePotential = bondPotential + bendPotential +
      torsionPotential +  inversionPotential;
    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
    curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
    curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
    curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
    curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;

  }

  void ForceManager::calcLongRangeInteraction() {
    Snapshot* curSnapshot;
    DataStorage* config;
    RealType* frc;
    RealType* pos;
    RealType* vel;
    RealType* trq;
    RealType* A;
    RealType* electroFrame;
    RealType* rc;
    RealType* vc;
    RealType* particlePot;

    //get current snapshot from SimInfo
    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();

    //get array pointers
    config = &(curSnapshot->atomData);
    frc = config->getArrayPointer(DataStorage::dslForce);
    pos = config->getArrayPointer(DataStorage::dslPosition);
    vel = config->getArrayPointer(DataStorage::dslVelocity);
    trq = config->getArrayPointer(DataStorage::dslTorque);
    A   = config->getArrayPointer(DataStorage::dslAmat);
    electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
    particlePot = config->getArrayPointer(DataStorage::dslParticlePot);

    //calculate the center of mass of cutoff group
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    Molecule::CutoffGroupIterator ci;
    CutoffGroup* cg;
    Vector3d com;
    Vector3d comv;
    std::vector<Vector3d> rcGroup;
    std::vector<Vector3d> vcGroup;

    if(info_->getNCutoffGroups() > 0){

      for (mol = info_->beginMolecule(mi); mol != NULL;
           mol = info_->nextMolecule(mi)) {
        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
            cg = mol->nextCutoffGroup(ci)) {
              cg->getCOM(com);
              cg->getCOMV(comv);
              rcGroup.push_back(com);
              vcGroup.push_back(comv);
        }
      }// end for (mol)

      rc = rcGroup[0].getArrayPointer();
      vc = vcGroup[0].getArrayPointer();
    } else {
      // center of mass of the group is the same as position of the atom
      // if cutoff group does not exist
      rc = pos;
      vc = vel;
    }

    //initialize data before passing to fortran
    RealType longRangePotential[LR_POT_TYPES];
    RealType lrPot = 0.0;
    Vector3d totalDipole;
    Jv_ = 0.0;

    
    int isError = 0;

    for (int i=0; i<LR_POT_TYPES;i++){
      longRangePotential[i]=0.0; //Initialize array
    }

    doForceLoop(pos,
                vel,
                rc,
                vc,
                A,
                electroFrame,
                frc,
                trq,
                tau.getArrayPointer(),
                Jv_.getArrayPointer(),
                longRangePotential,
                particlePot,
                &isError );

    if( isError ){
      sprintf( painCave.errMsg,
         "Error returned from the fortran force calculation.\n" );
      painCave.isFatal = 1;
      simError();
    }
    for (int i=0; i<LR_POT_TYPES;i++){
      lrPot += longRangePotential[i]; //Quick hack
    }

    // grab the simulation box dipole moment if specified
    if (info_->getCalcBoxDipole()){
      getAccumulatedBoxDipole(totalDipole.getArrayPointer());

      curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
      curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
      curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
    }

    //store the tau and long range potential
    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
  }


  void ForceManager::postCalculation() {
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    Molecule::RigidBodyIterator rbIter;
    RigidBody* rb;
    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();

    // collect the atomic forces onto rigid bodies

    for (mol = info_->beginMolecule(mi); mol != NULL;
         mol = info_->nextMolecule(mi)) {
      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
           rb = mol->nextRigidBody(rbIter)) {
        Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
        tau += rbTau;
      }
    }

#ifdef IS_MPI
    Mat3x3d tmpTau(tau);
    MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
                  9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
#endif
    curSnapshot->statData.setTau(tau);
    curSnapshot->statData.setJv(Jv_);
  }

} //end namespace OpenMD
Revision:	1684
Committed:	Thu Mar 1 19:04:10 2012 UTC (13 years, 4 months ago) by chuckv
File size:	12712 byte(s)
Log Message:	Fixed bug with rigidbodies not updating velocities for heat flux calculation.
#	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	chuckv	1671	*
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	chuckv	1671
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			#include "UseTheForce/doForces_interface.h"
53	gezelter	1390	#define __OPENMD_C
54	chuckv	664	#include "UseTheForce/DarkSide/fInteractionMap.h"
55	gezelter	246	#include "utils/simError.h"
56	xsun	1215	#include "primitives/Bond.hpp"
57	tim	749	#include "primitives/Bend.hpp"
58	cli2	1275	#include "primitives/Torsion.hpp"
59			#include "primitives/Inversion.hpp"
60	gezelter	1390	namespace OpenMD {
61	gezelter	246
62	gezelter	1464	void ForceManager::calcForces() {
63	chuckv	1671
64	gezelter	246	if (!info_->isFortranInitialized()) {
65	gezelter	507	info_->update();
66	gezelter	246	}
67	chuckv	1671
68	gezelter	246	preCalculation();
69	chuckv	1671
70	gezelter	246	calcShortRangeInteraction();
71
72	gezelter	1464	calcLongRangeInteraction();
73	gezelter	246
74	gezelter	1464	postCalculation();
75	chuckv	1671
76	gezelter	507	}
77	chuckv	1671
78	gezelter	507	void ForceManager::preCalculation() {
79	gezelter	246	SimInfo::MoleculeIterator mi;
80			Molecule* mol;
81			Molecule::AtomIterator ai;
82			Atom* atom;
83			Molecule::RigidBodyIterator rbIter;
84			RigidBody* rb;
85	chuckv	1671
86	gezelter	246	// forces are zeroed here, before any are accumulated.
87			// NOTE: do not rezero the forces in Fortran.
88	chuckv	1671
89			for (mol = info_->beginMolecule(mi); mol != NULL;
90	gezelter	1126	mol = info_->nextMolecule(mi)) {
91	gezelter	507	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
92	chuckv	1671	atom->zeroForcesAndTorques();
93	gezelter	507	}
94	chuckv	1671
95	gezelter	507	//change the positions of atoms which belong to the rigidbodies
96	chuckv	1671	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
97	gezelter	1126	rb = mol->nextRigidBody(rbIter)) {
98	chuckv	1671	rb->zeroForcesAndTorques();
99			}
100
101	gezelter	246	}
102	chuckv	1671
103	gezelter	1126	// Zero out the stress tensor
104			tau *= 0.0;
105	chuckv	1671
106	gezelter	507	}
107	chuckv	1671
108	gezelter	507	void ForceManager::calcShortRangeInteraction() {
109	gezelter	246	Molecule* mol;
110			RigidBody* rb;
111			Bond* bond;
112			Bend* bend;
113			Torsion* torsion;
114	cli2	1275	Inversion* inversion;
115	gezelter	246	SimInfo::MoleculeIterator mi;
116			Molecule::RigidBodyIterator rbIter;
117			Molecule::BondIterator bondIter;;
118			Molecule::BendIterator bendIter;
119			Molecule::TorsionIterator torsionIter;
120	cli2	1275	Molecule::InversionIterator inversionIter;
121	tim	963	RealType bondPotential = 0.0;
122			RealType bendPotential = 0.0;
123			RealType torsionPotential = 0.0;
124	cli2	1275	RealType inversionPotential = 0.0;
125	gezelter	246
126	chuckv	1671	//calculate short range interactions
127			for (mol = info_->beginMolecule(mi); mol != NULL;
128	gezelter	1126	mol = info_->nextMolecule(mi)) {
129	gezelter	246
130	gezelter	507	//change the positions of atoms which belong to the rigidbodies
131	chuckv	1671	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
132	gezelter	1126	rb = mol->nextRigidBody(rbIter)) {
133			rb->updateAtoms();
134	chuckv	1684	rb->updateAtomVel();
135	gezelter	507	}
136	gezelter	246
137	chuckv	1671	for (bond = mol->beginBond(bondIter); bond != NULL;
138	gezelter	1126	bond = mol->nextBond(bondIter)) {
139	tim	749	bond->calcForce();
140			bondPotential += bond->getPotential();
141	gezelter	507	}
142	gezelter	246
143	chuckv	1671	for (bend = mol->beginBend(bendIter); bend != NULL;
144	gezelter	1126	bend = mol->nextBend(bendIter)) {
145	chuckv	1671
146	gezelter	1126	RealType angle;
147			bend->calcForce(angle);
148	chuckv	1671	RealType currBendPot = bend->getPotential();
149
150	gezelter	1126	bendPotential += bend->getPotential();
151			std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
152			if (i == bendDataSets.end()) {
153			BendDataSet dataSet;
154			dataSet.prev.angle = dataSet.curr.angle = angle;
155			dataSet.prev.potential = dataSet.curr.potential = currBendPot;
156			dataSet.deltaV = 0.0;
157			bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
158			}else {
159			i->second.prev.angle = i->second.curr.angle;
160			i->second.prev.potential = i->second.curr.potential;
161			i->second.curr.angle = angle;
162			i->second.curr.potential = currBendPot;
163	chuckv	1671	i->second.deltaV = fabs(i->second.curr.potential -
164	gezelter	1126	i->second.prev.potential);
165			}
166	gezelter	507	}
167	chuckv	1671
168			for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
169	gezelter	1126	torsion = mol->nextTorsion(torsionIter)) {
170	tim	963	RealType angle;
171	gezelter	1126	torsion->calcForce(angle);
172	tim	963	RealType currTorsionPot = torsion->getPotential();
173	gezelter	1126	torsionPotential += torsion->getPotential();
174			std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
175			if (i == torsionDataSets.end()) {
176			TorsionDataSet dataSet;
177			dataSet.prev.angle = dataSet.curr.angle = angle;
178			dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
179			dataSet.deltaV = 0.0;
180			torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
181			}else {
182			i->second.prev.angle = i->second.curr.angle;
183			i->second.prev.potential = i->second.curr.potential;
184			i->second.curr.angle = angle;
185			i->second.curr.potential = currTorsionPot;
186	chuckv	1671	i->second.deltaV = fabs(i->second.curr.potential -
187	gezelter	1126	i->second.prev.potential);
188	chuckv	1671	}
189			}
190	cli2	1275
191	chuckv	1671	for (inversion = mol->beginInversion(inversionIter);
192			inversion != NULL;
193	cli2	1275	inversion = mol->nextInversion(inversionIter)) {
194			RealType angle;
195			inversion->calcForce(angle);
196			RealType currInversionPot = inversion->getPotential();
197			inversionPotential += inversion->getPotential();
198			std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
199			if (i == inversionDataSets.end()) {
200			InversionDataSet dataSet;
201			dataSet.prev.angle = dataSet.curr.angle = angle;
202			dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
203			dataSet.deltaV = 0.0;
204			inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
205			}else {
206			i->second.prev.angle = i->second.curr.angle;
207			i->second.prev.potential = i->second.curr.potential;
208			i->second.curr.angle = angle;
209			i->second.curr.potential = currInversionPot;
210	chuckv	1671	i->second.deltaV = fabs(i->second.curr.potential -
211	cli2	1275	i->second.prev.potential);
212	chuckv	1671	}
213			}
214	gezelter	246	}
215	chuckv	1671
216			RealType shortRangePotential = bondPotential + bendPotential +
217			torsionPotential + inversionPotential;
218	gezelter	246	Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
219			curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
220	tim	665	curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
221			curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
222			curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
223	cli2	1275	curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
224	chuckv	1671
225	gezelter	507	}
226	chuckv	1671
227	gezelter	1464	void ForceManager::calcLongRangeInteraction() {
228	gezelter	246	Snapshot* curSnapshot;
229			DataStorage* config;
230	tim	963	RealType* frc;
231			RealType* pos;
232	chuckv	1671	RealType* vel;
233	tim	963	RealType* trq;
234			RealType* A;
235			RealType* electroFrame;
236			RealType* rc;
237	chuckv	1682	RealType* vc;
238	chuckv	1245	RealType* particlePot;
239	chuckv	1671
240	gezelter	246	//get current snapshot from SimInfo
241			curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
242	chuckv	1671
243	gezelter	246	//get array pointers
244			config = &(curSnapshot->atomData);
245			frc = config->getArrayPointer(DataStorage::dslForce);
246			pos = config->getArrayPointer(DataStorage::dslPosition);
247	chuckv	1671	vel = config->getArrayPointer(DataStorage::dslVelocity);
248	gezelter	246	trq = config->getArrayPointer(DataStorage::dslTorque);
249			A = config->getArrayPointer(DataStorage::dslAmat);
250			electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
251	chuckv	1245	particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
252	gezelter	246
253			//calculate the center of mass of cutoff group
254			SimInfo::MoleculeIterator mi;
255			Molecule* mol;
256			Molecule::CutoffGroupIterator ci;
257			CutoffGroup* cg;
258			Vector3d com;
259	chuckv	1682	Vector3d comv;
260	gezelter	246	std::vector<Vector3d> rcGroup;
261	chuckv	1682	std::vector<Vector3d> vcGroup;
262	chuckv	1671
263	gezelter	246	if(info_->getNCutoffGroups() > 0){
264	chuckv	1671
265			for (mol = info_->beginMolecule(mi); mol != NULL;
266	gezelter	1126	mol = info_->nextMolecule(mi)) {
267	chuckv	1671	for(cg = mol->beginCutoffGroup(ci); cg != NULL;
268	gezelter	1126	cg = mol->nextCutoffGroup(ci)) {
269	chuckv	1682	cg->getCOM(com);
270			cg->getCOMV(comv);
271			rcGroup.push_back(com);
272			vcGroup.push_back(comv);
273	gezelter	246	}
274	gezelter	507	}// end for (mol)
275	chuckv	1671
276	gezelter	507	rc = rcGroup[0].getArrayPointer();
277	chuckv	1682	vc = vcGroup[0].getArrayPointer();
278	gezelter	246	} else {
279	chuckv	1671	// center of mass of the group is the same as position of the atom
280	gezelter	1126	// if cutoff group does not exist
281	gezelter	507	rc = pos;
282	chuckv	1682	vc = vel;
283	gezelter	246	}
284	chuckv	1671
285	gezelter	246	//initialize data before passing to fortran
286	tim	963	RealType longRangePotential[LR_POT_TYPES];
287			RealType lrPot = 0.0;
288	chrisfen	998	Vector3d totalDipole;
289	chuckv	1681	Jv_ = 0.0;
290	chuckv	1671
291
292
293	gezelter	246	int isError = 0;
294
295	chuckv	664	for (int i=0; i<LR_POT_TYPES;i++){
296			longRangePotential[i]=0.0; //Initialize array
297			}
298	chuckv	1671
299	xsun	1215	doForceLoop(pos,
300	chuckv	1671	vel,
301	xsun	1215	rc,
302	chuckv	1682	vc,
303	xsun	1215	A,
304			electroFrame,
305			frc,
306			trq,
307	chuckv	1671	tau.getArrayPointer(),
308			Jv_.getArrayPointer(),
309			longRangePotential,
310	chuckv	1245	particlePot,
311	xsun	1215	&isError );
312	chuckv	1671
313	gezelter	246	if( isError ){
314	gezelter	507	sprintf( painCave.errMsg,
315	chuckv	1671	"Error returned from the fortran force calculation.\n" );
316	gezelter	507	painCave.isFatal = 1;
317			simError();
318	gezelter	246	}
319	chuckv	664	for (int i=0; i<LR_POT_TYPES;i++){
320			lrPot += longRangePotential[i]; //Quick hack
321			}
322	chuckv	1671
323	chrisfen	998	// grab the simulation box dipole moment if specified
324			if (info_->getCalcBoxDipole()){
325			getAccumulatedBoxDipole(totalDipole.getArrayPointer());
326	chuckv	1671
327	chrisfen	998	curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
328			curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
329			curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
330			}
331	chuckv	1671
332			//store the tau and long range potential
333	chuckv	664	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
334	chrisfen	691	curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
335	tim	681	curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
336	gezelter	507	}
337	gezelter	246
338	chuckv	1671
339	gezelter	1464	void ForceManager::postCalculation() {
340	gezelter	246	SimInfo::MoleculeIterator mi;
341			Molecule* mol;
342			Molecule::RigidBodyIterator rbIter;
343			RigidBody* rb;
344	gezelter	1126	Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
345	chuckv	1671
346	gezelter	246	// collect the atomic forces onto rigid bodies
347	chuckv	1671
348			for (mol = info_->beginMolecule(mi); mol != NULL;
349	gezelter	1126	mol = info_->nextMolecule(mi)) {
350	chuckv	1671	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
351			rb = mol->nextRigidBody(rbIter)) {
352	gezelter	1464	Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
353			tau += rbTau;
354	gezelter	507	}
355	gezelter	1126	}
356	chuckv	1671
357	gezelter	1126	#ifdef IS_MPI
358	gezelter	1464	Mat3x3d tmpTau(tau);
359	chuckv	1671	MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
360	gezelter	1464	9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
361	gezelter	1126	#endif
362	gezelter	1464	curSnapshot->statData.setTau(tau);
363	chuckv	1671	curSnapshot->statData.setJv(Jv_);
364	gezelter	507	}
365	gezelter	246
366	gezelter	1390	} //end namespace OpenMD