src/brains/ForceFields.cpp

#include <iostream>

using namespace std;


#include <stdlib.h>

#ifdef IS_MPI
#include <mpi.h>
#endif // is_mpi


#ifdef PROFILE
#include "profiling/mdProfile.hpp"
#endif

#include "utils/simError.h"
#include "UseTheForce/ForceFields.hpp"
#include "primitives/Atom.hpp"
#include "UseTheForce/doForces_interface.h"

void ForceFields::calcRcut( void ){

#ifdef IS_MPI
  double tempBig = bigSigma;
  MPI_Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX,
                 MPI_COMM_WORLD);
#endif  //is_mpi

  //calc rCut and rList

  entry_plug->setDefaultRcut( 2.5 * bigSigma );  
    
}

void ForceFields::setRcut( double LJrcut ) {
  
#ifdef IS_MPI
  double tempBig = bigSigma;
  MPI_Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX,
                 MPI_COMM_WORLD);
#endif  //is_mpi
  
  if (LJrcut < 2.5 * bigSigma) {
    sprintf( painCave.errMsg,
             "Setting Lennard-Jones cutoff radius to %lf.\n"
             "\tThis value is smaller than %lf, which is\n"
             "\t2.5 * bigSigma, where bigSigma is the largest\n"
             "\tvalue of sigma present in the simulation.\n"
             "\tThis is potentially a problem since the LJ potential may\n"
             "\tbe appreciable at this distance.  If you don't want the\n"
             "\tsmaller cutoff, change the LJrcut variable.\n",
             LJrcut, 2.5*bigSigma);
    painCave.isFatal = 0;
    simError();
  } else {
    sprintf( painCave.errMsg,
             "Setting Lennard-Jones cutoff radius to %lf.\n"
             "\tThis value is larger than %lf, which is\n"
             "\t2.5 * bigSigma, where bigSigma is the largest\n"
             "\tvalue of sigma present in the simulation. This should\n"
             "\tnot be a problem, but could adversely effect performance.\n",
             LJrcut, 2.5*bigSigma);
    painCave.isFatal = 0;
    simError();
  }
  
  //calc rCut and rList
  
  entry_plug->setDefaultRcut( LJrcut );
}

void ForceFields::doForces( int calcPot, int calcStress ){

  int i, j, isError;
  double* frc;
  double* pos;
  double* trq;
  double* A;
  double* u_l;
  double* rc;
  double* massRatio;
  double factor;
  SimState* config;

  Molecule* myMols;
  Atom** myAtoms;
  int numAtom;
  int curIndex;
  double mtot;
  int numMol;
  int numCutoffGroups;
  CutoffGroup* myCutoffGroup;
  vector<CutoffGroup*>::iterator iterCutoff;
  double com[3];
  vector<double> rcGroup;
  
  short int passedCalcPot = (short int)calcPot;
  short int passedCalcStress = (short int)calcStress;

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

  for(i=0; i<entry_plug->n_atoms; i++){
    entry_plug->atoms[i]->zeroForces();    
  }

#ifdef PROFILE
  startProfile(pro7);
#endif
  
  for(i=0; i<entry_plug->n_mol; i++ ){
    // CalcForces in molecules takes care of mapping rigid body coordinates
    // into atomic coordinates
    entry_plug->molecules[i].calcForces();    
  }

#ifdef PROFILE
  endProfile( pro7 );
#endif

  config = entry_plug->getConfiguration();
  
  frc = config->getFrcArray();
  pos = config->getPosArray();
  trq = config->getTrqArray();
  A   = config->getAmatArray();
  u_l = config->getUlArray();

  if(entry_plug->haveCutoffGroups){
    myMols = entry_plug->molecules;
    numMol = entry_plug->n_mol;
    for(int i  = 0; i < numMol; i++){
        
      numCutoffGroups = myMols[i].getNCutoffGroups();
      for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); myCutoffGroup != NULL; 
                                                    myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
        //get center of mass of the cutoff group
        myCutoffGroup->getCOM(com); 

        rcGroup.push_back(com[0]);
        rcGroup.push_back(com[1]);
        rcGroup.push_back(com[2]);
        
      }// end for(myCutoffGroup)
      
    }//end for(int i = 0)

    rc = &rcGroup[0];
  }
  else{
    // center of mass of the group is the same as position of the atom  if cutoff group does not exist
    rc = pos;
  }

  isError = 0;
  entry_plug->lrPot = 0.0;

  for (i=0; i<9; i++) {
    entry_plug->tau[i] = 0.0;
  }


#ifdef PROFILE
  startProfile(pro8);
#endif

  doForceLoop( pos,
               rc,
               A,
               u_l,
               frc,
               trq,
               entry_plug->tau,
               &(entry_plug->lrPot), 
               &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();
  }

  // scale forces if thermodynamic integration is used
  if (entry_plug->useSolidThermInt || entry_plug->useLiquidThermInt) {
    factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
    for (i=0; i < entry_plug->n_atoms; i++) {
      for (j=0; j< 3; j++) 
        frc[3*i + j] *= factor; 
      if (entry_plug->atoms[i]->isDirectional()) {
        for (j=0; j< 3; j++) 
          trq[3*i + j] *= factor;
      }
    }
    entry_plug->vRaw = entry_plug->lrPot;
    entry_plug->lrPot *= factor;
  }

  // collect the atomic forces onto rigid bodies
  for(i=0; i<entry_plug->n_mol; i++ ){
    entry_plug->molecules[i].atoms2rigidBodies();
  }

  // do crystal restraint forces for thermodynamic integration
  if (entry_plug->useSolidThermInt){
    entry_plug->lrPot += entry_plug->restraint->Calc_Restraint_Forces(entry_plug->integrableObjects);
    entry_plug->vHarm = entry_plug->restraint->getVharm();
  }
  

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

}


void ForceFields::initFortran(int useReactionField ){
  
  int isError;
  
  isError = 0;
  initFortranFF(&useReactionField, &isError );
  
  if(isError){
    sprintf( painCave.errMsg,
             "ForceField error: There was an error initializing the forceField in fortran.\n" );
    painCave.isFatal = 1;
    simError();
  }

  
#ifdef IS_MPI
  sprintf( checkPointMsg, "ForceField successfully initialized the fortran component list.\n" );
  MPIcheckPoint();
#endif // is_mpi
  
}


Revision:	221
Committed:	Tue Nov 23 22:48:31 2004 UTC (20 years, 5 months ago) by chrisfen
File size:	6148 byte(s)
Log Message:	Improvements to restraints
#	User	Rev	Content
1	gezelter	2	#include <iostream>
2
3			using namespace std;
4
5
6			#include <stdlib.h>
7
8			#ifdef IS_MPI
9			#include <mpi.h>
10			#endif // is_mpi
11
12
13			#ifdef PROFILE
14	tim	3	#include "profiling/mdProfile.hpp"
15	gezelter	2	#endif
16
17	tim	3	#include "utils/simError.h"
18			#include "UseTheForce/ForceFields.hpp"
19			#include "primitives/Atom.hpp"
20	chuckv	124	#include "UseTheForce/doForces_interface.h"
21	gezelter	2
22			void ForceFields::calcRcut( void ){
23
24			#ifdef IS_MPI
25			double tempBig = bigSigma;
26			MPI_Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX,
27			MPI_COMM_WORLD);
28			#endif //is_mpi
29
30			//calc rCut and rList
31
32			entry_plug->setDefaultRcut( 2.5 * bigSigma );
33
34			}
35
36			void ForceFields::setRcut( double LJrcut ) {
37
38			#ifdef IS_MPI
39			double tempBig = bigSigma;
40			MPI_Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX,
41			MPI_COMM_WORLD);
42			#endif //is_mpi
43
44			if (LJrcut < 2.5 * bigSigma) {
45			sprintf( painCave.errMsg,
46			"Setting Lennard-Jones cutoff radius to %lf.\n"
47			"\tThis value is smaller than %lf, which is\n"
48			"\t2.5 * bigSigma, where bigSigma is the largest\n"
49			"\tvalue of sigma present in the simulation.\n"
50			"\tThis is potentially a problem since the LJ potential may\n"
51			"\tbe appreciable at this distance. If you don't want the\n"
52			"\tsmaller cutoff, change the LJrcut variable.\n",
53			LJrcut, 2.5*bigSigma);
54			painCave.isFatal = 0;
55			simError();
56			} else {
57			sprintf( painCave.errMsg,
58			"Setting Lennard-Jones cutoff radius to %lf.\n"
59			"\tThis value is larger than %lf, which is\n"
60			"\t2.5 * bigSigma, where bigSigma is the largest\n"
61			"\tvalue of sigma present in the simulation. This should\n"
62			"\tnot be a problem, but could adversely effect performance.\n",
63			LJrcut, 2.5*bigSigma);
64			painCave.isFatal = 0;
65			simError();
66			}
67
68			//calc rCut and rList
69
70			entry_plug->setDefaultRcut( LJrcut );
71			}
72
73			void ForceFields::doForces( int calcPot, int calcStress ){
74
75			int i, j, isError;
76			double* frc;
77			double* pos;
78			double* trq;
79			double* A;
80			double* u_l;
81			double* rc;
82			double* massRatio;
83			double factor;
84			SimState* config;
85
86			Molecule* myMols;
87			Atom** myAtoms;
88			int numAtom;
89			int curIndex;
90			double mtot;
91			int numMol;
92			int numCutoffGroups;
93			CutoffGroup* myCutoffGroup;
94			vector<CutoffGroup*>::iterator iterCutoff;
95			double com[3];
96			vector<double> rcGroup;
97
98			short int passedCalcPot = (short int)calcPot;
99			short int passedCalcStress = (short int)calcStress;
100
101			// forces are zeroed here, before any are accumulated.
102			// NOTE: do not rezero the forces in Fortran.
103
104			for(i=0; i<entry_plug->n_atoms; i++){
105			entry_plug->atoms[i]->zeroForces();
106			}
107
108			#ifdef PROFILE
109			startProfile(pro7);
110			#endif
111
112			for(i=0; i<entry_plug->n_mol; i++ ){
113			// CalcForces in molecules takes care of mapping rigid body coordinates
114			// into atomic coordinates
115			entry_plug->molecules[i].calcForces();
116			}
117
118			#ifdef PROFILE
119			endProfile( pro7 );
120			#endif
121
122			config = entry_plug->getConfiguration();
123
124			frc = config->getFrcArray();
125			pos = config->getPosArray();
126			trq = config->getTrqArray();
127			A = config->getAmatArray();
128			u_l = config->getUlArray();
129
130			if(entry_plug->haveCutoffGroups){
131			myMols = entry_plug->molecules;
132			numMol = entry_plug->n_mol;
133			for(int i = 0; i < numMol; i++){
134
135			numCutoffGroups = myMols[i].getNCutoffGroups();
136			for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); myCutoffGroup != NULL;
137			myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
138			//get center of mass of the cutoff group
139			myCutoffGroup->getCOM(com);
140
141			rcGroup.push_back(com[0]);
142			rcGroup.push_back(com[1]);
143			rcGroup.push_back(com[2]);
144
145			}// end for(myCutoffGroup)
146
147			}//end for(int i = 0)
148
149			rc = &rcGroup[0];
150			}
151			else{
152			// center of mass of the group is the same as position of the atom if cutoff group does not exist
153			rc = pos;
154			}
155
156			isError = 0;
157			entry_plug->lrPot = 0.0;
158
159			for (i=0; i<9; i++) {
160			entry_plug->tau[i] = 0.0;
161			}
162
163
164			#ifdef PROFILE
165			startProfile(pro8);
166			#endif
167
168	gezelter	204	doForceLoop( pos,
169			rc,
170			A,
171			u_l,
172			frc,
173			trq,
174			entry_plug->tau,
175			&(entry_plug->lrPot),
176			&passedCalcPot,
177			&passedCalcStress,
178			&isError );
179	gezelter	2
180			#ifdef PROFILE
181			endProfile(pro8);
182			#endif
183
184
185			if( isError ){
186			sprintf( painCave.errMsg,
187			"Error returned from the fortran force calculation.\n" );
188			painCave.isFatal = 1;
189			simError();
190			}
191
192			// scale forces if thermodynamic integration is used
193			if (entry_plug->useSolidThermInt \|\| entry_plug->useLiquidThermInt) {
194			factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
195			for (i=0; i < entry_plug->n_atoms; i++) {
196			for (j=0; j< 3; j++)
197			frc[3i + j] = factor;
198			if (entry_plug->atoms[i]->isDirectional()) {
199			for (j=0; j< 3; j++)
200			trq[3i + j] = factor;
201			}
202			}
203			entry_plug->vRaw = entry_plug->lrPot;
204			entry_plug->lrPot *= factor;
205			}
206
207			// collect the atomic forces onto rigid bodies
208			for(i=0; i<entry_plug->n_mol; i++ ){
209			entry_plug->molecules[i].atoms2rigidBodies();
210			}
211
212			// do crystal restraint forces for thermodynamic integration
213			if (entry_plug->useSolidThermInt){
214			entry_plug->lrPot += entry_plug->restraint->Calc_Restraint_Forces(entry_plug->integrableObjects);
215			entry_plug->vHarm = entry_plug->restraint->getVharm();
216			}
217
218
219			#ifdef IS_MPI
220			sprintf( checkPointMsg,
221			"returned from the force calculation.\n" );
222			MPIcheckPoint();
223			#endif // is_mpi
224
225
226			}
227
228
229	gezelter	204	void ForceFields::initFortran(int useReactionField ){
230	gezelter	2
231			int isError;
232
233			isError = 0;
234	gezelter	204	initFortranFF(&useReactionField, &isError );
235	gezelter	2
236			if(isError){
237			sprintf( painCave.errMsg,
238			"ForceField error: There was an error initializing the forceField in fortran.\n" );
239			painCave.isFatal = 1;
240			simError();
241			}
242
243
244			#ifdef IS_MPI
245			sprintf( checkPointMsg, "ForceField successfully initialized the fortran component list.\n" );
246			MPIcheckPoint();
247			#endif // is_mpi
248
249			}
250
251