src/brains/ForceManager.cpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */
 
 /**
  * @file ForceManager.cpp
  * @author tlin
  * @date 11/09/2004
  * @time 10:39am
  * @version 1.0
  */

#include "brains/ForceManager.hpp"

namespace oopse {

void ForceManager::calcForces(bool needPotential, bool needStress) {

    Snapshot* curSnapshot;
    DataStorage* config;
    double* frc;
    double* pos;
    double* trq;
    double* A;
    double* u_l;
    double* rc;

    std::vector<Molecule*>::iterator mi;
    Molecule* mol;
    std::vector<Atom*>::iterator ai;
    Atom* atom;
    std::vector<Atom*>::iterator ci;
    CutoffGroup* cg;
    Vector3d com;
    std::vector<Vector3d> rcGroup;

    double longRangePotential;
    double tau[9];
    short int passedCalcPot;
    short int passedCalcStress;
    int isError;     

    // 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->zeroForces();
        }
    }// end for (mol)

#ifdef PROFILE
    startProfile(pro7);
#endif
    //calculate short range interactions    
    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
        mol->calcForces();
    }// end for (mol)

#ifdef PROFILE
    endProfile( pro7 );
#endif

    //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);
    trq = config->getArrayPointer(DataStorage::dslTorque);
    A   = config->getArrayPointer(DataStorage::dslAmat);
    u_l = config->getArrayPointer(DataStorage::dslUnitVector);

    //calculate the center of mass of cutoff group
    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)
            rcGroup.push_back(com);
        }
    }// end for (mol)
       
        rc = rcGroup[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;
    }
  

    //initialize data before passing to fortran
    longRangePotential = 0.0;
    for (int i=0; i<9; i++) {
        tau[i] = 0.0;
    }    
    passedCalcPot = (short int)needPotential;
    passedCalcStress = (short int)needStress;
    isError = 0;

#ifdef PROFILE
    startProfile(pro8);
#endif

    doForceLoop( pos,
            rc,
            A,
            u_l,
            frc,
            trq,
            tau,
            &longRangePotential, 
            &passedCalcPot,
            &passedCalcStress,
            &isError );

#ifdef PROFILE
    endProfile(pro8);
#endif

    if( isError ){
        sprintf( painCave.errMsg,
             "Error returned from the fortran force calculation.\n" );
        painCave.isFatal = 1;
        simError();
    }

    //store the tau and long range potential    
    curSnapshot->statData[Stats::LONGRANGEPOTENTIAL] = longRangePotential;
    curSnapshot->statData[Stats::TAUXX] = tau[0];
    curSnapshot->statData[Stats::TAUXY] = tau[1];
    curSnapshot->statData[Stats::TAUXZ] = tau[2];
    curSnapshot->statData[Stats::TAUYX] = tau[3];
    curSnapshot->statData[Stats::TAUYY] = tau[4];
    curSnapshot->statData[Stats::TAUYZ] = tau[5];
    curSnapshot->statData[Stats::TAUZX] = tau[6];
    curSnapshot->statData[Stats::TAUZY] = tau[7];
    curSnapshot->statData[Stats::TAUZZ] = tau[8];

    // collect the atomic forces onto rigid bodies
    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
        mol->atoms2rigidBodies();
    }

#ifdef IS_MPI
    sprintf( checkPointMsg,
       "returned from the force calculation.\n" );
    MPIcheckPoint();
#endif // is_mpi

}


} //end namespace oopse
Revision:	1722
Committed:	Tue Nov 9 23:11:39 2004 UTC (19 years, 9 months ago) by tim
File size:	5226 byte(s)
Log Message:	adding ForceManager
#	User	Rev	Content
1	tim	1722	/*
2			* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3			*
4			* Contact: oopse@oopse.org
5			*
6			* This program is free software; you can redistribute it and/or
7			* modify it under the terms of the GNU Lesser General Public License
8			* as published by the Free Software Foundation; either version 2.1
9			* of the License, or (at your option) any later version.
10			* All we ask is that proper credit is given for our work, which includes
11			* - but is not limited to - adding the above copyright notice to the beginning
12			* of your source code files, and to any copyright notice that you may distribute
13			* with programs based on this work.
14			*
15			* This program is distributed in the hope that it will be useful,
16			* but WITHOUT ANY WARRANTY; without even the implied warranty of
17			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18			* GNU Lesser General Public License for more details.
19			*
20			* You should have received a copy of the GNU Lesser General Public License
21			* along with this program; if not, write to the Free Software
22			* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23			*
24			*/
25
26			/**
27			* @file ForceManager.cpp
28			* @author tlin
29			* @date 11/09/2004
30			* @time 10:39am
31			* @version 1.0
32			*/
33
34			#include "brains/ForceManager.hpp"
35
36			namespace oopse {
37
38			void ForceManager::calcForces(bool needPotential, bool needStress) {
39
40			Snapshot* curSnapshot;
41			DataStorage* config;
42			double* frc;
43			double* pos;
44			double* trq;
45			double* A;
46			double* u_l;
47			double* rc;
48
49			std::vector<Molecule*>::iterator mi;
50			Molecule* mol;
51			std::vector<Atom*>::iterator ai;
52			Atom* atom;
53			std::vector<Atom*>::iterator ci;
54			CutoffGroup* cg;
55			Vector3d com;
56			std::vector<Vector3d> rcGroup;
57
58			double longRangePotential;
59			double tau[9];
60			short int passedCalcPot;
61			short int passedCalcStress;
62			int isError;
63
64			// forces are zeroed here, before any are accumulated.
65			// NOTE: do not rezero the forces in Fortran.
66			for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
67			for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
68			atom->zeroForces();
69			}
70			}// end for (mol)
71
72			#ifdef PROFILE
73			startProfile(pro7);
74			#endif
75			//calculate short range interactions
76			for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
77			mol->calcForces();
78			}// end for (mol)
79
80			#ifdef PROFILE
81			endProfile( pro7 );
82			#endif
83
84			//get current snapshot from SimInfo
85			curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot()
86
87			//get array pointers
88			config = &(curSnapshot->atomData);
89			frc = config->getArrayPointer(DataStorage::dslForce);
90			pos = config->getArrayPointer(DataStorage::dslPosition);
91			trq = config->getArrayPointer(DataStorage::dslTorque);
92			A = config->getArrayPointer(DataStorage::dslAmat);
93			u_l = config->getArrayPointer(DataStorage::dslUnitVector);
94
95			//calculate the center of mass of cutoff group
96			if(info_->getNCutoffGroups() > 0){
97
98			for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
99			for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
100			cg->getCOM(com)
101			rcGroup.push_back(com);
102			}
103			}// end for (mol)
104
105			rc = rcGroup[0]->getArrayPointer();
106			}
107			else{
108			// center of mass of the group is the same as position of the atom if cutoff group does not exist
109			rc = pos;
110			}
111
112
113			//initialize data before passing to fortran
114			longRangePotential = 0.0;
115			for (int i=0; i<9; i++) {
116			tau[i] = 0.0;
117			}
118			passedCalcPot = (short int)needPotential;
119			passedCalcStress = (short int)needStress;
120			isError = 0;
121
122			#ifdef PROFILE
123			startProfile(pro8);
124			#endif
125
126			doForceLoop( pos,
127			rc,
128			A,
129			u_l,
130			frc,
131			trq,
132			tau,
133			&longRangePotential,
134			&passedCalcPot,
135			&passedCalcStress,
136			&isError );
137
138			#ifdef PROFILE
139			endProfile(pro8);
140			#endif
141
142			if( isError ){
143			sprintf( painCave.errMsg,
144			"Error returned from the fortran force calculation.\n" );
145			painCave.isFatal = 1;
146			simError();
147			}
148
149			//store the tau and long range potential
150			curSnapshot->statData[Stats::LONGRANGEPOTENTIAL] = longRangePotential;
151			curSnapshot->statData[Stats::TAUXX] = tau[0];
152			curSnapshot->statData[Stats::TAUXY] = tau[1];
153			curSnapshot->statData[Stats::TAUXZ] = tau[2];
154			curSnapshot->statData[Stats::TAUYX] = tau[3];
155			curSnapshot->statData[Stats::TAUYY] = tau[4];
156			curSnapshot->statData[Stats::TAUYZ] = tau[5];
157			curSnapshot->statData[Stats::TAUZX] = tau[6];
158			curSnapshot->statData[Stats::TAUZY] = tau[7];
159			curSnapshot->statData[Stats::TAUZZ] = tau[8];
160
161			// collect the atomic forces onto rigid bodies
162			for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
163			mol->atoms2rigidBodies();
164			}
165
166			#ifdef IS_MPI
167			sprintf( checkPointMsg,
168			"returned from the force calculation.\n" );
169			MPIcheckPoint();
170			#endif // is_mpi
171
172			}
173
174
175			} //end namespace oopse