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)( int*, int*, int* ),
                        void (*p3)( double* positionArray,double* forceArray,
                                    double* potentialEnergy,
                                    short int* doPotentialCalc ) ){
  
  
  newLJtype = p1;
  initLJfortran = p2;
  currentLJwrap->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( painCave.errMsg,
                 "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;
    
    fastForward( "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->next; //skip the first element who is a place holder.
    while( currentAtomType != NULL ){
      currentAtomType->duplicate( info );

 
      sendFrcStruct( &info, mpiAtomStructType );

      sprintf( checkPointMsg,
               "successfully sent lJ force type: \"%s\"\n",
               info.name );
      MPIcheckPoint();

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

      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
  
  double bigSigma = 0.0;
  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();

    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
  }

  
#ifdef IS_MPI
  double tempBig = bigSigma;
  MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
#endif  //is_mpi

  //calc rCut and rList

  entry_plug->rCut = 2.5 * bigSigma;
  if(entry_plug->rCut > (entry_plug->box_x / 2.0)) entry_plug->rCut = entry_plug->box_x / 2.0;
  if(entry_plug->rCut > (entry_plug->box_y / 2.0)) entry_plug->rCut = entry_plug->box_y / 2.0;
  if(entry_plug->rCut > (entry_plug->box_z / 2.0)) entry_plug->rCut = entry_plug->box_z / 2.0;

  entry_plug->rList = entry_plug->rCut + 1.0;

  // clean up the memory
  
  delete headAtomType;

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

  initFortran();
  entry_plug->refreshSim();

}

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( int calcPot ){

  int i;
  double* frc;
  double* pos;
  short int passedCalcPot = (short int)calcPot;

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

//   entry_plug->lrPot = -1;
  doLJfortran( pos, frc, &(entry_plug->lrPot), &passedCalcPot );


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