mdtools/libmdCode/LJ_FF.cpp

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

#include <iostream>
using namespace std;

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

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