mdtools/interface_implementation/LJ_FF.cpp

#include <cstdlib>
#include <cstdio>
#include <cstring>

#include <iostream>
using namespace std;

#include "ForceFields.hpp"
#include "SRI.hpp"
#include "simError.h"


// Declare the structures that will be passed by the parser and  MPI

typedef struct{
  char name[15];
  double mass;
  double epslon;
  double sigma;
  int ident;
  int last;      //  0  -> default
                 //  1  -> in MPI: tells nodes to stop listening
} atomStruct;

int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info );

#ifdef IS_MPI

#include "mpiForceField.h"

MPI_Datatype mpiAtomStructType;

#endif


// declaration of functions needed to wrap the fortran module

extern "C" {
  void forcefactory_( char* forceName,
                      int* status,
                      void (*wrapFunction)( void (*p1)( int* ident, 
                                                        double* mass, 
                                                        double* epslon, 
                                                        double* sigma, 
                                                        int* status ),
                                            void (*p2)( int *nLocal,
                                                        int *identArray,
                                                        int *isError ),
                                            void (*p3)( double* positionArray,
                                                        double* forceArray,
                                                        double* potentialEnergy,
                                                        short int* doPotentialCalc )),
                      int forceNameLength );
}


void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon, 
                                   double* sigma, int* status ),
                        void (*p2)( int *nLocal, int *identArray, int *isError ),
                        void (*p3)( double* positionArray,double* forceArray,
                                    double* potentialEnergy,
                                    short int* doPotentialCalc ) );

void (*newLJtype)( int* ident, double* mass, double* epslon, double* sigma,
                   int* status );

void (*initLJfortran) ( int *nLocal, int *identArray, int *isError )

LJ_FF* currentLJwrap;


//****************************************************************
// begins the actual forcefield stuff.  
//****************************************************************

LJ_FF::LJ_FF(){

  char fileName[200];
  char* ffPath_env = "FORCE_PARAM_PATH";
  char* ffPath;
  char temp[200];
  char errMsg[1000];

  // do the funtion wrapping
  currentLJwrap = this;
  wrapMe();

#ifdef IS_MPI
  int i;
  
   // **********************************************************************
  // Init the atomStruct mpi type

  atomStruct atomProto; // mpiPrototype
  int atomBC[3] = {15,3,2};  // block counts
  MPI_Aint atomDspls[3];           // displacements
  MPI_Datatype atomMbrTypes[3];    // member mpi types

  MPI_Address(&atomProto.name, &atomDspls[0]);
  MPI_Address(&atomProto.mass, &atomDspls[1]);
  MPI_Address(&atomProto.ident, &atomDspls[2]);
  
  atomMbrTypes[0] = MPI_CHAR;
  atomMbrTypes[1] = MPI_DOUBLE;
  atomMbrTypes[2] = MPI_INT;
  
  for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
  
  MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
  MPI_Type_commit(&mpiAtomStructType);

  // ***********************************************************************
  
  if( worldRank == 0 ){
#endif
    
    // generate the force file name
    
    strcpy( fileName, "LJ_FF.frc" ); 
    //    fprintf( stderr,"Trying to open %s\n", fileName );
    
    // attempt to open the file in the current directory first.
    
    frcFile = fopen( fileName, "r" );
    
    if( frcFile == NULL ){
      
      // next see if the force path enviorment variable is set
      
      ffPath = getenv( ffPath_env );
      if( ffPath == NULL ) {
        sprintf( painCave.errMsg,
                 "Error opening the force field parameter file: %s\n"
                 "Have you tried setting the FORCE_PARAM_PATH environment "
                 "vairable?\n",
                 fileName );
        painCave.isFatal = 1;
        simError();
      }
      
      
      strcpy( temp, ffPath );
      strcat( temp, "/" );
      strcat( temp, fileName );
      strcpy( fileName, temp );
      
      frcFile = fopen( fileName, "r" );
      
      if( frcFile == NULL ){
        
        sprintf( painCave.errMsg,
                 "Error opening the force field parameter file: %s\n"
                 "Have you tried setting the FORCE_PARAM_PATH environment "
                 "vairable?\n",
                 fileName );
        painCave.isFatal = 1;
        simError();
      }
    }
    
#ifdef IS_MPI
  }
  
  sprintf( checkPointMsg, "LJ_FF file opened sucessfully." );
  MPIcheckPoint();
  
#endif // is_mpi
}


LJ_FF::~LJ_FF(){

#ifdef IS_MPI
  if( worldRank == 0 ){
#endif // is_mpi
    
    fclose( frcFile );
    
#ifdef IS_MPI
  }
#endif // is_mpi
}


void LJ_FF::wrapMe( void ){
  
  char* currentFF = "LJ";
  int isError = 0;
  
  forcefactory_( currentFF, &isError, LJfunctionWrapper, strlen(currentFF) );

  if( isError ){
    
    sprintf( painCave.errMsg,
             "LJ_FF error: an error was returned from fortran when the "
             "the functions were being wrapped.\n" );
    painCave.isFatal = 1;
    simError();
  }

#ifdef IS_MPI
  sprintf( checkPointMsg, "LJ_FF functions succesfully wrapped." );
  MPIcheckPoint();
#endif // is_mpi
}
  

void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon, 
                                   double* sigma, int* status ),
                        void (*p2)( void ),
                        void (*p3)( double* positionArray,double* forceArray,
                                    double* potentialEnergy,
                                    short int* doPotentialCalc ) ){
  
  
  p1 = newLJtype;
  p2 = initLJfortran;
  this->setLJfortran( p3 );
}


void LJ_FF::initializeAtoms( void ){
  
  class LinkedType {
  public:
    LinkedType(){ 
      next = NULL;
      name[0] = '\0';
    }
    ~LinkedType(){ if( next != NULL ) delete next; }

    LinkedType* find(char* key){ 
      if( !strcmp(name, key) ) return this;
      if( next != NULL ) return next->find(key);
      return NULL;
    }
    

    void add( atomStruct &info ){
    
      // check for duplicates
      
      if( !strcmp( info.name, name ) ){
        sprintf( simError.painCave,
                 "Duplicate LJ atom type \"%s\" found in "
                 "the LJ_FF param file./n",
                 name );
        painCave.isFatal = 1;
        simError();
      }
      
      if( next != NULL ) next->add(info);
      else{
        next = new LinkedType();
        strcpy(next->name, info.name);
        next->mass     = info.mass;
        next->epslon   = info.epslon;
        next->sigma    = info.sigma;
        next->ident    = info.ident;
      }
    }
    

#ifdef IS_MPI
    
    void duplicate( atomStruct &info ){
      strcpy(info.name, name);
      info.mass     = mass;
      info.epslon   = epslon;
      info.sigma    = sigma;
      info.ident    = ident;
      info.last     = 0;
    }


#endif

    char name[15];
    double mass;
    double epslon;
    double sigma;
    int ident;
    LinkedType* next;
  };
  
  LinkedType* headAtomType; 
  LinkedType* currentAtomType;
  atomStruct info;
  info.last = 1; // initialize last to have the last set. 
                 // if things go well, last will be set to 0

  int i;
  int identNum;
  
  Atom** the_atoms;
  int nAtoms;
  the_atoms = entry_plug->atoms;
  nAtoms = entry_plug->n_atoms;
  
  
#ifdef IS_MPI
  if( worldRank == 0 ){
#endif
    
    // read in the atom types.

    headAtomType = new LinkedType;
    
    fastFoward( "AtomTypes", "initializeAtoms" );

    // we are now at the AtomTypes section.
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
    
    
    // read a line, and start parseing out the atom types 

    if( eof_test == NULL ){
      sprintf( painCave.errMsg, 
               "Error in reading Atoms from force file at line %d.\n",
               lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    identNum = 1;
    // stop reading at end of file, or at next section
    while( readLine[0] != '#' && eof_test != NULL ){

      // toss comment lines
      if( readLine[0] != '!' ){
        
        // the parser returns 0 if the line was blank
        if( parseAtomLJ( readLine, lineNum, info ) ){
          info.ident = identNum;
          headAtomType->add( info );;
          identNum++;
        }
      }
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }

#ifdef IS_MPI
    
    // send out the linked list to all the other processes

    sprintf( checkPointMsg,
             "LJ_FF atom structures read successfully." );
    MPIcheckPoint();

    currentAtomType = headAtomType;
    while( currentAtomType != NULL ){
      currentAtomType->duplicate( info );
      sendFrcStruct( &info, mpiAtomStructType );
      currentAtomType = currentAtomType->next;
    }
    info.last = 1;
    sendFrcStruct( &info, mpiAtomStructType );
    
  }

  else{
    
    // listen for node 0 to send out the force params
    
    MPIcheckPoint();

    headAtomType = new LinkedType;
    recieveFrcStruct( &info, mpiAtomStructType );
    while( !info.last ){

      headAtomType->add( info );
      recieveFrcStruct( &info, mpiAtomStructType );
    }
  }
#endif // is_mpi

  // call new A_types in fortran
  
  int isError;
  currentAtomType = headAtomType;
  while( currentAtomType != NULL ){
    
    if( currentAtomType->name[0] != '\0' ){
      isError = 0;
      newLJtype( &(currentAtomType->ident),
                 &(currentAtomType->mass),
                 &(currentAtomType->epslon),
                 &(currentAtomType->sigma),
                 &isError );
      if( isError ){
        sprintf( painCave.errMsg,
                 "Error initializing the \"%s\" atom type in fortran\n",
                 currentAtomType->name );
        painCave.isFatal = 1;
        simError();
      }
    }
    currentAtomType = currentAtomType->next;
  }
      
#ifdef IS_MPI
  sprintf( checkPointMsg,
           "LJ_FF atom structures successfully sent to fortran\n" );
  MPIcheckPoint();
#endif // is_mpi

  
  // initialize the atoms
  
  Atom* thisAtom;

  for( i=0; i<nAtoms; i++ ){
    
    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
    if( currentAtomType == NULL ){
      sprintf( painCave.errMsg, 
               "AtomType error, %s not found in force file.\n",
               the_atoms[i]->getType() );
      painCave.isFatal = 1;
      simError();
    }
    
    the_atoms[i]->setMass( currentAtomType->mass );
    the_atoms[i]->setEpslon( currentAtomType->epslon );
    the_atoms[i]->setSigma( currentAtomType->sigma );
    the_atoms[i]->setIdent( currentAtomType->ident );
    the_atoms[i]->setLJ();
  }


  // clean up the memory
  
  delete headAtomType;

#ifdef IS_MPI
  sprintf( checkPointMsg, "LJ_FF atoms initialized succesfully" );
  MPIcheckPoint();
#endif // is_mpi

  initFortran();

}

void LJ_FF::initializeBonds( bond_pair* the_bonds ){
  
    if( entry_plug->n_bonds ){
      sprintf( painCave.errMsg,
               "LJ_FF does not support bonds.\n" );
      painCave.isFatal = 1;
      simError();
    }
#ifdef IS_MPI
  MPIcheckPoint();
#endif // is_mpi

}

void LJ_FF::initializeBends( bend_set* the_bends ){

    if( entry_plug->n_bends ){
      sprintf( painCave.errMsg,
               "LJ_FF does not support bends.\n" );
      painCave.isFatal = 1;
      simError();
    }
#ifdef IS_MPI
  MPIcheckPoint();
#endif // is_mpi

}

void LJ_FF::initializeTorsions( torsion_set* the_torsions ){

    if( entry_plug->n_torsions ){
      sprintf( painCave.errMsg,
               "LJ_FF does not support torsions.\n" );
      painCave.isFatal = 1;
      simError();
    }
#ifdef IS_MPI
  MPIcheckPoint();
#endif // is_mpi

}


void LJ_FF::fastForward( char* stopText, char* searchOwner ){

  int foundText = 0;
  char* the_token;

  rewind( frcFile );
  lineNum = 0;

  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  if( eof_test == NULL ){
    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
             " file is empty.\n",
             searchOwner );
    painCave.isFatal = 1;
    simError();
  }
  
  
  while( !foundText ){
    while( eof_test != NULL && readLine[0] != '#' ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }
    if( eof_test == NULL ){
      sprintf( painCave.errMsg,
               "Error fast forwarding force file for %s at "
               "line %d: file ended unexpectedly.\n",
               searchOwner,
               lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    the_token = strtok( readLine, " ,;\t#\n" );
    foundText = !strcmp( stopText, the_token );
    
    if( !foundText ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
      
      if( eof_test == NULL ){
        sprintf( painCave.errMsg, 
                 "Error fast forwarding force file for %s at "
                 "line %d: file ended unexpectedly.\n",
                 searchOwner,
                 lineNum );
        painCave.isFatal = 1;
        simError();
      } 
    }
  }  
}


int parseAtomLJ( char *lineBuffer, int lineNum,  atomStruct &info ){

  char* the_token;
  
  the_token = strtok( lineBuffer, " \n\t,;" );
  if( the_token != NULL ){
    
    strcpy( info.name, the_token );
          
    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
      sprintf( painCave.errMsg, 
               "Error parseing AtomTypes: line %d\n", lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    info.mass = atof( the_token );
    
    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
      sprintf( painCave.errMsg, 
               "Error parseing AtomTypes: line %d\n", lineNum );
      painCave.isFatal = 1;
      simError();
    }
        
    info.epslon = atof( the_token );
          
    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
      sprintf( painCave.errMsg, 
               "Error parseing AtomTypes: line %d\n", lineNum );
      painCave.isFatal = 1;
      simError();
    }
        
    info.sigma = atof( the_token );
    
    return 1;
  }
  else return 0;
}


void LJ_FF::doForces( void ){

  int i;
  double* frc;
  double* pos;
  double potE;
  short int calcPot = 0;

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

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

  frc = Atom::getFrcArray();
  pos = Atom::getPosArray();

  doLJfortran( pos, frc, potE, calcPot );
}
  
void LJ_FF::initFortran( void ){
  
  int nLocal = entry_plug->n_atoms;
  int *ident;
  int isError;
  int i;

  ident = new int[nLocal];

  for(i=0; i<nLocal; i++){
    ident[i] = entryplug->atoms[i]->getIdent();
  }

  isError = 0;
  initLJfortran( &nLocal, ident, &isError );
  
  if(isError){
    sprintf( painCave.errMsg,
             "LJ_FF error: There was an error initializing the component list in fortran.\n" );
    painCave.isFatal = 1;
    simError();
  }

  
#ifdef IS_MPI
  sprintf( checkPointMsg, "LJ_FF successfully initialized the fortran component list.\n" );
  MPIcheckPoint();
#endif // is_mpi
  
  delete[] ident;

}
  
Revision:	240
Committed:	Wed Jan 22 21:45:20 2003 UTC (21 years, 5 months ago) by chuckv
File size:	14156 byte(s)
Log Message:	Added init function to c++ force module.
#	User	Rev	Content
1	chuckv	215	#include <cstdlib>
2			#include <cstdio>
3			#include <cstring>
4
5			#include <iostream>
6			using namespace std;
7
8			#include "ForceFields.hpp"
9			#include "SRI.hpp"
10			#include "simError.h"
11
12
13	mmeineke	224	// Declare the structures that will be passed by the parser and MPI
14	chuckv	215
15			typedef struct{
16			char name[15];
17			double mass;
18			double epslon;
19			double sigma;
20	mmeineke	224	int ident;
21	chuckv	215	int last; // 0 -> default
22	mmeineke	224	// 1 -> in MPI: tells nodes to stop listening
23	chuckv	215	} atomStruct;
24	mmeineke	224
25			int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info );
26
27			#ifdef IS_MPI
28
29			#include "mpiForceField.h"
30
31	chuckv	215	MPI_Datatype mpiAtomStructType;
32
33			#endif
34
35
36	mmeineke	234	// declaration of functions needed to wrap the fortran module
37
38			extern "C" {
39			void forcefactory_( char* forceName,
40			int* status,
41	mmeineke	236	void (wrapFunction)( void (p1)( int* ident,
42			double* mass,
43			double* epslon,
44			double* sigma,
45			int* status ),
46	chuckv	240	void (p2)( int nLocal,
47			int *identArray,
48			int *isError ),
49	mmeineke	236	void (p3)( double positionArray,
50			double* forceArray,
51			double* potentialEnergy,
52			short int* doPotentialCalc )),
53			int forceNameLength );
54			}
55	mmeineke	234
56
57	mmeineke	236	void LJfunctionWrapper( void (p1)( int ident, double* mass, double* epslon,
58			double* sigma, int* status ),
59	chuckv	240	void (p2)( int nLocal, int identArray, int isError ),
60	mmeineke	236	void (p3)( double positionArray,double* forceArray,
61			double* potentialEnergy,
62			short int* doPotentialCalc ) );
63	mmeineke	234
64	mmeineke	236	void (newLJtype)( int ident, double* mass, double* epslon, double* sigma,
65			int* status );
66
67	chuckv	240	void (initLJfortran) ( int nLocal, int identArray, int isError )
68	mmeineke	236
69			LJ_FF* currentLJwrap;
70
71
72			//****************************************************************
73			// begins the actual forcefield stuff.
74			//****************************************************************
75
76	chuckv	215	LJ_FF::LJ_FF(){
77
78			char fileName[200];
79			char* ffPath_env = "FORCE_PARAM_PATH";
80			char* ffPath;
81			char temp[200];
82			char errMsg[1000];
83
84	mmeineke	236	// do the funtion wrapping
85			currentLJwrap = this;
86			wrapMe();
87
88	chuckv	215	#ifdef IS_MPI
89			int i;
90
91			// **********************************************************************
92			// Init the atomStruct mpi type
93
94			atomStruct atomProto; // mpiPrototype
95	mmeineke	224	int atomBC[3] = {15,3,2}; // block counts
96	chuckv	215	MPI_Aint atomDspls[3]; // displacements
97			MPI_Datatype atomMbrTypes[3]; // member mpi types
98
99			MPI_Address(&atomProto.name, &atomDspls[0]);
100			MPI_Address(&atomProto.mass, &atomDspls[1]);
101	mmeineke	224	MPI_Address(&atomProto.ident, &atomDspls[2]);
102	chuckv	215
103			atomMbrTypes[0] = MPI_CHAR;
104			atomMbrTypes[1] = MPI_DOUBLE;
105			atomMbrTypes[2] = MPI_INT;
106
107			for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
108
109			MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
110			MPI_Type_commit(&mpiAtomStructType);
111
112			// ***********************************************************************
113
114			if( worldRank == 0 ){
115			#endif
116
117			// generate the force file name
118
119			strcpy( fileName, "LJ_FF.frc" );
120			// fprintf( stderr,"Trying to open %s\n", fileName );
121
122			// attempt to open the file in the current directory first.
123
124			frcFile = fopen( fileName, "r" );
125
126			if( frcFile == NULL ){
127
128			// next see if the force path enviorment variable is set
129
130			ffPath = getenv( ffPath_env );
131			if( ffPath == NULL ) {
132			sprintf( painCave.errMsg,
133			"Error opening the force field parameter file: %s\n"
134			"Have you tried setting the FORCE_PARAM_PATH environment "
135			"vairable?\n",
136			fileName );
137			painCave.isFatal = 1;
138			simError();
139			}
140
141
142			strcpy( temp, ffPath );
143			strcat( temp, "/" );
144			strcat( temp, fileName );
145			strcpy( fileName, temp );
146
147			frcFile = fopen( fileName, "r" );
148
149			if( frcFile == NULL ){
150
151			sprintf( painCave.errMsg,
152			"Error opening the force field parameter file: %s\n"
153			"Have you tried setting the FORCE_PARAM_PATH environment "
154			"vairable?\n",
155			fileName );
156			painCave.isFatal = 1;
157			simError();
158			}
159			}
160
161			#ifdef IS_MPI
162			}
163
164			sprintf( checkPointMsg, "LJ_FF file opened sucessfully." );
165			MPIcheckPoint();
166
167			#endif // is_mpi
168			}
169
170
171			LJ_FF::~LJ_FF(){
172
173			#ifdef IS_MPI
174			if( worldRank == 0 ){
175			#endif // is_mpi
176
177			fclose( frcFile );
178
179			#ifdef IS_MPI
180			}
181			#endif // is_mpi
182			}
183
184	mmeineke	236
185			void LJ_FF::wrapMe( void ){
186
187			char* currentFF = "LJ";
188			int isError = 0;
189
190			forcefactory_( currentFF, &isError, LJfunctionWrapper, strlen(currentFF) );
191
192			if( isError ){
193
194			sprintf( painCave.errMsg,
195			"LJ_FF error: an error was returned from fortran when the "
196			"the functions were being wrapped.\n" );
197			painCave.isFatal = 1;
198			simError();
199			}
200
201			#ifdef IS_MPI
202			sprintf( checkPointMsg, "LJ_FF functions succesfully wrapped." );
203			MPIcheckPoint();
204			#endif // is_mpi
205			}
206
207
208			void LJfunctionWrapper( void (p1)( int ident, double* mass, double* epslon,
209			double* sigma, int* status ),
210			void (*p2)( void ),
211			void (p3)( double positionArray,double* forceArray,
212			double* potentialEnergy,
213			short int* doPotentialCalc ) ){
214
215
216			p1 = newLJtype;
217			p2 = initLJfortran;
218			this->setLJfortran( p3 );
219			}
220
221
222
223	chuckv	215	void LJ_FF::initializeAtoms( void ){
224
225			class LinkedType {
226			public:
227			LinkedType(){
228			next = NULL;
229			name[0] = '\0';
230			}
231			~LinkedType(){ if( next != NULL ) delete next; }
232
233			LinkedType* find(char* key){
234			if( !strcmp(name, key) ) return this;
235			if( next != NULL ) return next->find(key);
236			return NULL;
237			}
238
239	mmeineke	224
240	chuckv	215	void add( atomStruct &info ){
241	mmeineke	224
242			// check for duplicates
243
244			if( !strcmp( info.name, name ) ){
245			sprintf( simError.painCave,
246			"Duplicate LJ atom type \"%s\" found in "
247			"the LJ_FF param file./n",
248			name );
249			painCave.isFatal = 1;
250			simError();
251			}
252
253	chuckv	215	if( next != NULL ) next->add(info);
254			else{
255			next = new LinkedType();
256			strcpy(next->name, info.name);
257			next->mass = info.mass;
258			next->epslon = info.epslon;
259			next->sigma = info.sigma;
260	mmeineke	224	next->ident = info.ident;
261	chuckv	215	}
262			}
263
264	mmeineke	224
265			#ifdef IS_MPI
266
267	chuckv	215	void duplicate( atomStruct &info ){
268			strcpy(info.name, name);
269			info.mass = mass;
270			info.epslon = epslon;
271			info.sigma = sigma;
272	mmeineke	224	info.ident = ident;
273	chuckv	215	info.last = 0;
274			}
275
276
277			#endif
278
279			char name[15];
280			double mass;
281			double epslon;
282			double sigma;
283	mmeineke	224	int ident;
284	chuckv	215	LinkedType* next;
285			};
286
287			LinkedType* headAtomType;
288			LinkedType* currentAtomType;
289			atomStruct info;
290			info.last = 1; // initialize last to have the last set.
291			// if things go well, last will be set to 0
292
293			int i;
294	mmeineke	224	int identNum;
295	chuckv	215
296			Atom** the_atoms;
297			int nAtoms;
298			the_atoms = entry_plug->atoms;
299			nAtoms = entry_plug->n_atoms;
300
301
302			#ifdef IS_MPI
303			if( worldRank == 0 ){
304			#endif
305
306			// read in the atom types.
307	mmeineke	224
308	chuckv	215	headAtomType = new LinkedType;
309
310	mmeineke	224	fastFoward( "AtomTypes", "initializeAtoms" );
311
312	chuckv	215	// we are now at the AtomTypes section.
313
314			eof_test = fgets( readLine, sizeof(readLine), frcFile );
315			lineNum++;
316
317	mmeineke	224
318			// read a line, and start parseing out the atom types
319
320	chuckv	215	if( eof_test == NULL ){
321			sprintf( painCave.errMsg,
322			"Error in reading Atoms from force file at line %d.\n",
323			lineNum );
324			painCave.isFatal = 1;
325			simError();
326			}
327
328	mmeineke	224	identNum = 1;
329			// stop reading at end of file, or at next section
330	chuckv	215	while( readLine[0] != '#' && eof_test != NULL ){
331	mmeineke	224
332			// toss comment lines
333	chuckv	215	if( readLine[0] != '!' ){
334
335	mmeineke	224	// the parser returns 0 if the line was blank
336			if( parseAtomLJ( readLine, lineNum, info ) ){
337			info.ident = identNum;
338			headAtomType->add( info );;
339			identNum++;
340	chuckv	215	}
341			}
342			eof_test = fgets( readLine, sizeof(readLine), frcFile );
343			lineNum++;
344			}
345
346			#ifdef IS_MPI
347
348			// send out the linked list to all the other processes
349
350			sprintf( checkPointMsg,
351			"LJ_FF atom structures read successfully." );
352			MPIcheckPoint();
353
354			currentAtomType = headAtomType;
355			while( currentAtomType != NULL ){
356			currentAtomType->duplicate( info );
357			sendFrcStruct( &info, mpiAtomStructType );
358			currentAtomType = currentAtomType->next;
359			}
360			info.last = 1;
361			sendFrcStruct( &info, mpiAtomStructType );
362
363			}
364
365			else{
366
367			// listen for node 0 to send out the force params
368
369			MPIcheckPoint();
370
371			headAtomType = new LinkedType;
372			recieveFrcStruct( &info, mpiAtomStructType );
373			while( !info.last ){
374
375			headAtomType->add( info );
376			recieveFrcStruct( &info, mpiAtomStructType );
377			}
378			}
379			#endif // is_mpi
380
381	chuckv	230	// call new A_types in fortran
382	mmeineke	238
383			int isError;
384			currentAtomType = headAtomType;
385			while( currentAtomType != NULL ){
386
387	chuckv	240	if( currentAtomType->name[0] != '\0' ){
388	mmeineke	238	isError = 0;
389			newLJtype( &(currentAtomType->ident),
390			&(currentAtomType->mass),
391			&(currentAtomType->epslon),
392			&(currentAtomType->sigma),
393	chuckv	240	&isError );
394	mmeineke	238	if( isError ){
395			sprintf( painCave.errMsg,
396			"Error initializing the \"%s\" atom type in fortran\n",
397			currentAtomType->name );
398			painCave.isFatal = 1;
399			simError();
400			}
401			}
402			currentAtomType = currentAtomType->next;
403			}
404
405			#ifdef IS_MPI
406			sprintf( checkPointMsg,
407			"LJ_FF atom structures successfully sent to fortran\n" );
408			MPIcheckPoint();
409			#endif // is_mpi
410	chuckv	230
411	mmeineke	238
412	chuckv	230
413	chuckv	215	// initialize the atoms
414
415			Atom* thisAtom;
416
417			for( i=0; i<nAtoms; i++ ){
418
419			currentAtomType = headAtomType->find( the_atoms[i]->getType() );
420			if( currentAtomType == NULL ){
421			sprintf( painCave.errMsg,
422			"AtomType error, %s not found in force file.\n",
423			the_atoms[i]->getType() );
424			painCave.isFatal = 1;
425			simError();
426			}
427
428			the_atoms[i]->setMass( currentAtomType->mass );
429			the_atoms[i]->setEpslon( currentAtomType->epslon );
430			the_atoms[i]->setSigma( currentAtomType->sigma );
431	mmeineke	224	the_atoms[i]->setIdent( currentAtomType->ident );
432	chuckv	215	the_atoms[i]->setLJ();
433			}
434
435
436			// clean up the memory
437
438			delete headAtomType;
439
440			#ifdef IS_MPI
441			sprintf( checkPointMsg, "LJ_FF atoms initialized succesfully" );
442			MPIcheckPoint();
443			#endif // is_mpi
444
445	chuckv	240	initFortran();
446
447	chuckv	215	}
448
449			void LJ_FF::initializeBonds( bond_pair* the_bonds ){
450
451			if( entry_plug->n_bonds ){
452			sprintf( painCave.errMsg,
453			"LJ_FF does not support bonds.\n" );
454			painCave.isFatal = 1;
455			simError();
456			}
457			#ifdef IS_MPI
458			MPIcheckPoint();
459			#endif // is_mpi
460
461			}
462
463			void LJ_FF::initializeBends( bend_set* the_bends ){
464
465			if( entry_plug->n_bends ){
466			sprintf( painCave.errMsg,
467			"LJ_FF does not support bends.\n" );
468			painCave.isFatal = 1;
469			simError();
470			}
471			#ifdef IS_MPI
472			MPIcheckPoint();
473			#endif // is_mpi
474
475			}
476
477			void LJ_FF::initializeTorsions( torsion_set* the_torsions ){
478
479			if( entry_plug->n_torsions ){
480			sprintf( painCave.errMsg,
481			"LJ_FF does not support torsions.\n" );
482			painCave.isFatal = 1;
483			simError();
484			}
485			#ifdef IS_MPI
486			MPIcheckPoint();
487			#endif // is_mpi
488
489			}
490	mmeineke	224
491
492			void LJ_FF::fastForward( char* stopText, char* searchOwner ){
493
494			int foundText = 0;
495			char* the_token;
496
497			rewind( frcFile );
498			lineNum = 0;
499
500			eof_test = fgets( readLine, sizeof(readLine), frcFile );
501			lineNum++;
502			if( eof_test == NULL ){
503			sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
504			" file is empty.\n",
505			searchOwner );
506			painCave.isFatal = 1;
507			simError();
508			}
509
510
511			while( !foundText ){
512			while( eof_test != NULL && readLine[0] != '#' ){
513			eof_test = fgets( readLine, sizeof(readLine), frcFile );
514			lineNum++;
515			}
516			if( eof_test == NULL ){
517			sprintf( painCave.errMsg,
518			"Error fast forwarding force file for %s at "
519			"line %d: file ended unexpectedly.\n",
520			searchOwner,
521			lineNum );
522			painCave.isFatal = 1;
523			simError();
524			}
525
526			the_token = strtok( readLine, " ,;\t#\n" );
527			foundText = !strcmp( stopText, the_token );
528
529			if( !foundText ){
530			eof_test = fgets( readLine, sizeof(readLine), frcFile );
531			lineNum++;
532
533			if( eof_test == NULL ){
534			sprintf( painCave.errMsg,
535			"Error fast forwarding force file for %s at "
536			"line %d: file ended unexpectedly.\n",
537			searchOwner,
538			lineNum );
539			painCave.isFatal = 1;
540			simError();
541			}
542			}
543			}
544			}
545
546
547
548			int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info ){
549
550			char* the_token;
551
552			the_token = strtok( lineBuffer, " \n\t,;" );
553			if( the_token != NULL ){
554
555			strcpy( info.name, the_token );
556
557			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
558			sprintf( painCave.errMsg,
559			"Error parseing AtomTypes: line %d\n", lineNum );
560			painCave.isFatal = 1;
561			simError();
562			}
563
564			info.mass = atof( the_token );
565
566			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
567			sprintf( painCave.errMsg,
568			"Error parseing AtomTypes: line %d\n", lineNum );
569			painCave.isFatal = 1;
570			simError();
571			}
572
573			info.epslon = atof( the_token );
574
575			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
576			sprintf( painCave.errMsg,
577			"Error parseing AtomTypes: line %d\n", lineNum );
578			painCave.isFatal = 1;
579			simError();
580			}
581
582			info.sigma = atof( the_token );
583
584			return 1;
585			}
586			else return 0;
587			}
588	mmeineke	238
589
590			void LJ_FF::doForces( void ){
591
592			int i;
593			double* frc;
594			double* pos;
595			double potE;
596			short int calcPot = 0;
597
598			// forces are zeroed here, before any are acumulated.
599			// NOTE: do not rezero the forces in Fortran.
600
601			for(i=0; i<entry_plug->n_atoms; i++){
602			entry_plug->atoms[i]->zeroForces();
603			}
604
605			frc = Atom::getFrcArray();
606			pos = Atom::getPosArray();
607
608			doLJfortran( pos, frc, potE, calcPot );
609			}
610
611	chuckv	240	void LJ_FF::initFortran( void ){
612
613			int nLocal = entry_plug->n_atoms;
614			int *ident;
615			int isError;
616			int i;
617
618			ident = new int[nLocal];
619
620			for(i=0; i<nLocal; i++){
621			ident[i] = entryplug->atoms[i]->getIdent();
622			}
623
624			isError = 0;
625			initLJfortran( &nLocal, ident, &isError );
626
627			if(isError){
628			sprintf( painCave.errMsg,
629			"LJ_FF error: There was an error initializing the component list in fortran.\n" );
630			painCave.isFatal = 1;
631			simError();
632			}
633
634
635			#ifdef IS_MPI
636			sprintf( checkPointMsg, "LJ_FF successfully initialized the fortran component list.\n" );
637			MPIcheckPoint();
638			#endif // is_mpi
639
640			delete[] ident;
641
642			}
643