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, isError;
  double* frc;
  double* pos;
  double* trq;
  double* A;
  double* u_l;;
  SimState* config;

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

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


#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
  
}
Revision:	1097
Committed:	Mon Apr 12 20:32:20 2004 UTC (20 years, 2 months ago) by gezelter
File size:	4100 byte(s)
Log Message:	Changes for RigidBody dynamics (Somewhat extensive)
#	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			int i, isError;
77			double* frc;
78			double* pos;
79			double* trq;
80			double* A;
81	mmeineke	443	double* u_l;;
82	mmeineke	670	SimState* config;
83	mmeineke	377
84			short int passedCalcPot = (short int)calcPot;
85			short int passedCalcStress = (short int)calcStress;
86
87	gezelter	484	// forces are zeroed here, before any are accumulated.
88	mmeineke	377	// NOTE: do not rezero the forces in Fortran.
89
90			for(i=0; i<entry_plug->n_atoms; i++){
91	gezelter	484	entry_plug->atoms[i]->zeroForces();
92	mmeineke	377	}
93
94	chuckv	892	#ifdef PROFILE
95			startProfile(pro7);
96			#endif
97
98	mmeineke	423	for(i=0; i<entry_plug->n_mol; i++ ){
99	gezelter	1097	// CalcForces in molecules takes care of mapping rigid body coordinates
100			// into atomic coordinates
101	mmeineke	428	entry_plug->molecules[i].calcForces();
102	mmeineke	389	}
103
104	chuckv	892	#ifdef PROFILE
105			endProfile( pro7 );
106			#endif
107
108	mmeineke	670	config = entry_plug->getConfiguration();
109	mmeineke	597
110	mmeineke	670	frc = config->getFrcArray();
111			pos = config->getPosArray();
112			trq = config->getTrqArray();
113			A = config->getAmatArray();
114			u_l = config->getUlArray();
115	mmeineke	561
116	mmeineke	377	isError = 0;
117			entry_plug->lrPot = 0.0;
118	mmeineke	393
119	chuckv	479	for (i=0; i<9; i++) {
120			entry_plug->tau[i] = 0.0;
121			}
122	chuckv	438
123	mmeineke	841
124	chuckv	892	#ifdef PROFILE
125			startProfile(pro8);
126			#endif
127
128	mmeineke	377	fortranForceLoop( pos,
129			A,
130			u_l,
131			frc,
132			trq,
133	chuckv	479	entry_plug->tau,
134	mmeineke	377	&(entry_plug->lrPot),
135			&passedCalcPot,
136			&passedCalcStress,
137			&isError );
138
139	chuckv	892	#ifdef PROFILE
140			endProfile(pro8);
141			#endif
142	mmeineke	841
143
144	mmeineke	377	if( isError ){
145			sprintf( painCave.errMsg,
146			"Error returned from the fortran force calculation.\n" );
147			painCave.isFatal = 1;
148			simError();
149			}
150
151	gezelter	1097	for(i=0; i<entry_plug->n_mol; i++ ){
152			entry_plug->molecules[i].atoms2rigidBodies();
153			}
154
155
156	mmeineke	377	#ifdef IS_MPI
157			sprintf( checkPointMsg,
158			"returned from the force calculation.\n" );
159			MPIcheckPoint();
160			#endif // is_mpi
161	mmeineke	597
162	mmeineke	377
163			}
164
165
166			void ForceFields::initFortran(int ljMixPolicy, int useReactionField ){
167
168			int isError;
169
170			isError = 0;
171			initFortranFF( &ljMixPolicy, &useReactionField, &isError );
172
173			if(isError){
174			sprintf( painCave.errMsg,
175			"ForceField error: There was an error initializing the forceField in fortran.\n" );
176			painCave.isFatal = 1;
177			simError();
178			}
179
180
181			#ifdef IS_MPI
182			sprintf( checkPointMsg, "ForceField successfully initialized the fortran component list.\n" );
183			MPIcheckPoint();
184			#endif // is_mpi
185
186			}