OOPSE/libmdtools/ForceFields.cpp

#include <iostream>

using namespace std;


#include <stdlib.h>

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


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

#include "simError.h"
#include "ForceFields.hpp"
#include "Atom.hpp"
#include "fortranWrappers.hpp"


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

  fortranForceLoop( 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();
  }

  for(i=0; i<entry_plug->n_mol; i++ ){
    entry_plug->molecules[i].atoms2rigidBodies();
  }


  if (entry_plug->useSolidThermInt && !entry_plug->useLiquidThermInt) {

    factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
    for (i=0; i < entry_plug->integrableObjects.size(); i++) {
      for (j=0; j< 3; j++) 
        frc[3*i + j] *= factor; 
      if (entry_plug->integrableObjects[i]->isDirectional()) {
        for (j=0; j< 3; j++) 
          trq[3*i + j] *= factor;
      }
    }
    entry_plug->vRaw = entry_plug->lrPot;
    entry_plug->lrPot *= factor;
    entry_plug->lrPot += entry_plug->restraint->Calc_Restraint_Forces(entry_plug->integrableObjects);
    entry_plug->vHarm = entry_plug->restraint->getVharm();
  }

  if (entry_plug->useLiquidThermInt) {

    factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
    for (i=0; i < entry_plug->integrableObjects.size(); i++) {
      for (j=0; j< 3; j++) 
        frc[3*i + j] *= factor; 
      if (entry_plug->integrableObjects[i]->isDirectional()) {
        for (j=0; j< 3; j++) 
          trq[3*i + j] *= factor;
      }
    }
    entry_plug->vRaw = entry_plug->lrPot;
    entry_plug->lrPot *= factor;
  }

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

}


void ForceFields::initFortran(int ljMixPolicy, int useReactionField ){
  
  int isError;
  
  isError = 0;
  initFortranFF( &ljMixPolicy, &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
  
}


void ForceFields::initRestraints(){
  int i;
  // store the initial info.
  // set the omega values to zero
  for (i=0; i<entry_plug->integrableObjects.size(); i++)
    entry_plug->integrableObjects[i]->setZangle( 0.0 );

  entry_plug->restraint->Store_Init_Info(entry_plug->integrableObjects);

}

void ForceFields::dumpzAngle(){

  // store the initial info.
  entry_plug->restraint->Write_zAngle_File(entry_plug->integrableObjects);

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