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"

#include <fortranWrappers.hpp>

#ifdef IS_MPI
#include "mpiForceField.h"
#endif // is_mpi


namespace {

  // 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 parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
  
#ifdef IS_MPI
  
  MPI_Datatype mpiAtomStructType;
  
#endif
}


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

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

  int i, isError;
  double* frc;
  double* pos;
  double* tau;
  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();
  tau = entry_plug->tau;

  isError = 0;
  doLJfortran( pos, frc, &(entry_plug->lrPot), tau, &passedCalcPot, &isError );


  if( isError ){
    sprintf( painCave.errMsg,
             "Error returned from the fortran force calculation.\n" );
    painCave.isFatal = 1;
    simError();
  }

#ifdef IS_MPI
  sprintf( checkPointMsg,
           "returned from the force calculation.\n" );
  MPIcheckPoint();
#endif // is_mpi

}
  
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;

}

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