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->useThermInt) {
    myStunts = entry_plug->integrableObjects;

    factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
    for (i=0; i < myStunts.size(); i++) {
      for (j=0; j< 3; j++) 
        frc[3*i + j] *= factor; 
      if (myStunts[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(myStunts);
    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 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(){

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

}

void ForceFields::dumpzAngle(){

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

}
Revision:	1180
Committed:	Thu May 20 20:24:07 2004 UTC (20 years, 1 month ago) by chrisfen
File size:	6040 byte(s)
Log Message:	Several additions... Restraints.cpp and .hpp were included for restraining particles in thermodynamic integration. By including these, changes were made in Integrator, SimInfo, ForceFeilds, SimSetup, StatWriter, and possibly some other files. Two bass keywords were also added for performing thermodynamic integration: a lambda value one and a k power one.
#	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	1180	if (entry_plug->useThermInt) {
202			myStunts = entry_plug->integrableObjects;
203
204			factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
205			for (i=0; i < myStunts.size(); i++) {
206			for (j=0; j< 3; j++)
207			frc[3i + j] = factor;
208			if (myStunts[i]->isDirectional()) {
209			for (j=0; j< 3; j++)
210			trq[3i + j] = factor;
211			}
212			}
213			entry_plug->vRaw = entry_plug->lrPot;
214			entry_plug->lrPot *= factor;
215			entry_plug->lrPot += entry_plug->restraint->Calc_Restraint_Forces(myStunts);
216			entry_plug->vHarm = entry_plug->restraint->getVharm();
217			}
218
219	mmeineke	377	#ifdef IS_MPI
220			sprintf( checkPointMsg,
221			"returned from the force calculation.\n" );
222			MPIcheckPoint();
223			#endif // is_mpi
224	mmeineke	597
225	mmeineke	377
226			}
227
228
229			void ForceFields::initFortran(int ljMixPolicy, int useReactionField ){
230
231			int isError;
232
233			isError = 0;
234			initFortranFF( &ljMixPolicy, &useReactionField, &isError );
235
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	chrisfen	1180
251
252			void ForceFields::initRestraints(){
253
254			// store the initial info.
255			entry_plug->restraint->Store_Init_Info();
256
257			}
258
259			void ForceFields::dumpzAngle(){
260
261			// store the initial info.
262			entry_plug->restraint->Write_zAngle_File();
263
264			}