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"


#ifdef IS_MPI

#include "mpiForceField.h"


// Declare the structures that will be passed by MPI

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

#endif


LJ_FF::LJ_FF(){

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

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

  atomStruct atomProto; // mpiPrototype
  int atomBC[3] = {15,3,1};  // 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.last, &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;
    }
    
#ifdef IS_MPI
    void add( atomStruct &info ){
      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;
      }
    }
    
    void duplicate( atomStruct &info ){
      strcpy(info.name, name);
      info.mass     = mass;
      info.epslon   = epslon;
      info.sigma    = sigma;
      info.last     = 0;
    }


#endif

    char name[15];
    double mass;
    double epslon;
    double sigma;
    LinkedType* next;
  };
  
  LinkedType* headAtomType; 
  LinkedType* currentAtomType;
  LinkedType* tempAtomType;

#ifdef IS_MPI
  atomStruct info;
  info.last = 1; // initialize last to have the last set. 
                 // if things go well, last will be set to 0
#endif
  

  char readLine[500];
  char* the_token;
  char* eof_test;
  int foundAtom = 0;
  int lineNum = 0;
  int i;

  
  //////////////////////////////////////////////////
  // a quick water fix 

  double waterI[3][3]; 
  waterI[0][0] = 1.76958347772500;
  waterI[0][1] = 0.0;
  waterI[0][2] = 0.0;

  waterI[1][0] = 0.0;
  waterI[1][1] = 0.614537057924513;
  waterI[1][2] = 0.0;

  waterI[2][0] = 0.0;
  waterI[2][1] = 0.0;
  waterI[2][2] = 1.15504641980049;


  double headI[3][3]; 
  headI[0][0] = 1125;
  headI[0][1] = 0.0;
  headI[0][2] = 0.0;

  headI[1][0] = 0.0;
  headI[1][1] = 1125;
  headI[1][2] = 0.0;

  headI[2][0] = 0.0;
  headI[2][1] = 0.0;
  headI[2][2] = 250;

  
  //////////////////////////////////////////////////

  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.
    
    rewind( frcFile );
    headAtomType = new LinkedType;
    currentAtomType = headAtomType;
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
    if( eof_test == NULL ){
      sprintf( painCave.errMsg, "Error in reading Atoms from force file.\n" );
      painCave.isFatal = 1;
      simError();
    }
    
    
    while( !foundAtom ){
      while( eof_test != NULL && readLine[0] != '#' ){
        eof_test = fgets( readLine, sizeof(readLine), frcFile );
        lineNum++;
      }
      if( eof_test == NULL ){
        sprintf( painCave.errMsg, 
                 "Error in reading Atoms from force file at line %d.\n",
                 lineNum );
        painCave.isFatal = 1;
        simError();
      }
      
      the_token = strtok( readLine, " ,;\t#\n" );
      foundAtom = !strcmp( "AtomTypes", the_token );
      
      if( !foundAtom ){
        eof_test = fgets( readLine, sizeof(readLine), frcFile );
        lineNum++;
        
        if( eof_test == NULL ){
          sprintf( painCave.errMsg, 
                   "Error in reading Atoms from force file at line %d.\n",
                   lineNum );
          painCave.isFatal = 1;
          simError();
        } 
      }
    }
    
    // we are now at the AtomTypes section.
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
    
    if( eof_test == NULL ){
      sprintf( painCave.errMsg, 
               "Error in reading Atoms from force file at line %d.\n",
               lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    while( readLine[0] != '#' && eof_test != NULL ){
      
      if( readLine[0] != '!' ){
        
        the_token = strtok( readLine, " \n\t,;" );
        if( the_token != NULL ){
          
          strcpy( currentAtomType->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();
          }
          
          sscanf( the_token, "%lf", &currentAtomType->mass );
          
          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            sprintf( painCave.errMsg, 
                     "Error parseing AtomTypes: line %d\n", lineNum );
            painCave.isFatal = 1;
            simError();
          }
          
          sscanf( the_token, "%lf", &currentAtomType->epslon );
          
          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            sprintf( painCave.errMsg, 
                     "Error parseing AtomTypes: line %d\n", lineNum );
            painCave.isFatal = 1;
            simError();
          }
          
          sscanf( the_token, "%lf", &currentAtomType->sigma );
        }
      }
      
      tempAtomType = new LinkedType;
      currentAtomType->next = tempAtomType;
      currentAtomType = tempAtomType;
      
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }

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

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

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

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

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

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

  
  // initialize the atoms
  
  Atom* thisAtom;

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


  // clean up the memory
  
  delete headAtomType;

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

}

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

}
Revision:	215
Committed:	Thu Dec 19 21:59:51 2002 UTC (22 years, 10 months ago) by chuckv
File size:	9955 byte(s)
Log Message:	+ added lennard-jones force module and corresponding class. + created forceFactory directory.
#	User	Rev	Content
1	chuckv	215	#include <cstdlib>
2			#include <cstdio>
3			#include <cstring>
4
5			#include <iostream>
6			using namespace std;
7
8			#include "ForceFields.hpp"
9			#include "SRI.hpp"
10			#include "simError.h"
11
12
13			#ifdef IS_MPI
14
15			#include "mpiForceField.h"
16
17
18			// Declare the structures that will be passed by MPI
19
20			typedef struct{
21			char name[15];
22			double mass;
23			double epslon;
24			double sigma;
25			int last; // 0 -> default
26			// 1 -> tells nodes to stop listening
27			} atomStruct;
28			MPI_Datatype mpiAtomStructType;
29
30			#endif
31
32
33
34			LJ_FF::LJ_FF(){
35
36			char fileName[200];
37			char* ffPath_env = "FORCE_PARAM_PATH";
38			char* ffPath;
39			char temp[200];
40			char errMsg[1000];
41
42			#ifdef IS_MPI
43			int i;
44
45			// **********************************************************************
46			// Init the atomStruct mpi type
47
48			atomStruct atomProto; // mpiPrototype
49			int atomBC[3] = {15,3,1}; // block counts
50			MPI_Aint atomDspls[3]; // displacements
51			MPI_Datatype atomMbrTypes[3]; // member mpi types
52
53			MPI_Address(&atomProto.name, &atomDspls[0]);
54			MPI_Address(&atomProto.mass, &atomDspls[1]);
55			MPI_Address(&atomProto.last, &atomDspls[2]);
56
57			atomMbrTypes[0] = MPI_CHAR;
58			atomMbrTypes[1] = MPI_DOUBLE;
59			atomMbrTypes[2] = MPI_INT;
60
61			for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
62
63			MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
64			MPI_Type_commit(&mpiAtomStructType);
65
66			// ***********************************************************************
67
68			if( worldRank == 0 ){
69			#endif
70
71			// generate the force file name
72
73			strcpy( fileName, "LJ_FF.frc" );
74			// fprintf( stderr,"Trying to open %s\n", fileName );
75
76			// attempt to open the file in the current directory first.
77
78			frcFile = fopen( fileName, "r" );
79
80			if( frcFile == NULL ){
81
82			// next see if the force path enviorment variable is set
83
84			ffPath = getenv( ffPath_env );
85			if( ffPath == NULL ) {
86			sprintf( painCave.errMsg,
87			"Error opening the force field parameter file: %s\n"
88			"Have you tried setting the FORCE_PARAM_PATH environment "
89			"vairable?\n",
90			fileName );
91			painCave.isFatal = 1;
92			simError();
93			}
94
95
96			strcpy( temp, ffPath );
97			strcat( temp, "/" );
98			strcat( temp, fileName );
99			strcpy( fileName, temp );
100
101			frcFile = fopen( fileName, "r" );
102
103			if( frcFile == NULL ){
104
105			sprintf( painCave.errMsg,
106			"Error opening the force field parameter file: %s\n"
107			"Have you tried setting the FORCE_PARAM_PATH environment "
108			"vairable?\n",
109			fileName );
110			painCave.isFatal = 1;
111			simError();
112			}
113			}
114
115			#ifdef IS_MPI
116			}
117
118			sprintf( checkPointMsg, "LJ_FF file opened sucessfully." );
119			MPIcheckPoint();
120
121			#endif // is_mpi
122			}
123
124
125			LJ_FF::~LJ_FF(){
126
127			#ifdef IS_MPI
128			if( worldRank == 0 ){
129			#endif // is_mpi
130
131			fclose( frcFile );
132
133			#ifdef IS_MPI
134			}
135			#endif // is_mpi
136			}
137
138			void LJ_FF::initializeAtoms( void ){
139
140			class LinkedType {
141			public:
142			LinkedType(){
143			next = NULL;
144			name[0] = '\0';
145			}
146			~LinkedType(){ if( next != NULL ) delete next; }
147
148			LinkedType* find(char* key){
149			if( !strcmp(name, key) ) return this;
150			if( next != NULL ) return next->find(key);
151			return NULL;
152			}
153
154			#ifdef IS_MPI
155			void add( atomStruct &info ){
156			if( next != NULL ) next->add(info);
157			else{
158			next = new LinkedType();
159			strcpy(next->name, info.name);
160			next->mass = info.mass;
161			next->epslon = info.epslon;
162			next->sigma = info.sigma;
163			}
164			}
165
166			void duplicate( atomStruct &info ){
167			strcpy(info.name, name);
168			info.mass = mass;
169			info.epslon = epslon;
170			info.sigma = sigma;
171			info.last = 0;
172			}
173
174
175			#endif
176
177			char name[15];
178			double mass;
179			double epslon;
180			double sigma;
181			LinkedType* next;
182			};
183
184			LinkedType* headAtomType;
185			LinkedType* currentAtomType;
186			LinkedType* tempAtomType;
187
188			#ifdef IS_MPI
189			atomStruct info;
190			info.last = 1; // initialize last to have the last set.
191			// if things go well, last will be set to 0
192			#endif
193
194
195			char readLine[500];
196			char* the_token;
197			char* eof_test;
198			int foundAtom = 0;
199			int lineNum = 0;
200			int i;
201
202
203			//////////////////////////////////////////////////
204			// a quick water fix
205
206			double waterI[3][3];
207			waterI[0][0] = 1.76958347772500;
208			waterI[0][1] = 0.0;
209			waterI[0][2] = 0.0;
210
211			waterI[1][0] = 0.0;
212			waterI[1][1] = 0.614537057924513;
213			waterI[1][2] = 0.0;
214
215			waterI[2][0] = 0.0;
216			waterI[2][1] = 0.0;
217			waterI[2][2] = 1.15504641980049;
218
219
220			double headI[3][3];
221			headI[0][0] = 1125;
222			headI[0][1] = 0.0;
223			headI[0][2] = 0.0;
224
225			headI[1][0] = 0.0;
226			headI[1][1] = 1125;
227			headI[1][2] = 0.0;
228
229			headI[2][0] = 0.0;
230			headI[2][1] = 0.0;
231			headI[2][2] = 250;
232
233
234
235			//////////////////////////////////////////////////
236
237			Atom** the_atoms;
238			int nAtoms;
239			the_atoms = entry_plug->atoms;
240			nAtoms = entry_plug->n_atoms;
241
242
243			#ifdef IS_MPI
244			if( worldRank == 0 ){
245			#endif
246
247			// read in the atom types.
248
249			rewind( frcFile );
250			headAtomType = new LinkedType;
251			currentAtomType = headAtomType;
252
253			eof_test = fgets( readLine, sizeof(readLine), frcFile );
254			lineNum++;
255			if( eof_test == NULL ){
256			sprintf( painCave.errMsg, "Error in reading Atoms from force file.\n" );
257			painCave.isFatal = 1;
258			simError();
259			}
260
261
262			while( !foundAtom ){
263			while( eof_test != NULL && readLine[0] != '#' ){
264			eof_test = fgets( readLine, sizeof(readLine), frcFile );
265			lineNum++;
266			}
267			if( eof_test == NULL ){
268			sprintf( painCave.errMsg,
269			"Error in reading Atoms from force file at line %d.\n",
270			lineNum );
271			painCave.isFatal = 1;
272			simError();
273			}
274
275			the_token = strtok( readLine, " ,;\t#\n" );
276			foundAtom = !strcmp( "AtomTypes", the_token );
277
278			if( !foundAtom ){
279			eof_test = fgets( readLine, sizeof(readLine), frcFile );
280			lineNum++;
281
282			if( eof_test == NULL ){
283			sprintf( painCave.errMsg,
284			"Error in reading Atoms from force file at line %d.\n",
285			lineNum );
286			painCave.isFatal = 1;
287			simError();
288			}
289			}
290			}
291
292			// we are now at the AtomTypes section.
293
294			eof_test = fgets( readLine, sizeof(readLine), frcFile );
295			lineNum++;
296
297			if( eof_test == NULL ){
298			sprintf( painCave.errMsg,
299			"Error in reading Atoms from force file at line %d.\n",
300			lineNum );
301			painCave.isFatal = 1;
302			simError();
303			}
304
305			while( readLine[0] != '#' && eof_test != NULL ){
306
307			if( readLine[0] != '!' ){
308
309			the_token = strtok( readLine, " \n\t,;" );
310			if( the_token != NULL ){
311
312			strcpy( currentAtomType->name, the_token );
313
314			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
315			sprintf( painCave.errMsg,
316			"Error parseing AtomTypes: line %d\n", lineNum );
317			painCave.isFatal = 1;
318			simError();
319			}
320
321			sscanf( the_token, "%lf", &currentAtomType->mass );
322
323			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
324			sprintf( painCave.errMsg,
325			"Error parseing AtomTypes: line %d\n", lineNum );
326			painCave.isFatal = 1;
327			simError();
328			}
329
330			sscanf( the_token, "%lf", &currentAtomType->epslon );
331
332			if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
333			sprintf( painCave.errMsg,
334			"Error parseing AtomTypes: line %d\n", lineNum );
335			painCave.isFatal = 1;
336			simError();
337			}
338
339			sscanf( the_token, "%lf", &currentAtomType->sigma );
340			}
341			}
342
343			tempAtomType = new LinkedType;
344			currentAtomType->next = tempAtomType;
345			currentAtomType = tempAtomType;
346
347			eof_test = fgets( readLine, sizeof(readLine), frcFile );
348			lineNum++;
349			}
350
351			#ifdef IS_MPI
352
353			// send out the linked list to all the other processes
354
355			sprintf( checkPointMsg,
356			"LJ_FF atom structures read successfully." );
357			MPIcheckPoint();
358
359			currentAtomType = headAtomType;
360			while( currentAtomType != NULL ){
361			currentAtomType->duplicate( info );
362			sendFrcStruct( &info, mpiAtomStructType );
363			currentAtomType = currentAtomType->next;
364			}
365			info.last = 1;
366			sendFrcStruct( &info, mpiAtomStructType );
367
368			}
369
370			else{
371
372			// listen for node 0 to send out the force params
373
374			MPIcheckPoint();
375
376			headAtomType = new LinkedType;
377			recieveFrcStruct( &info, mpiAtomStructType );
378			while( !info.last ){
379
380			headAtomType->add( info );
381			recieveFrcStruct( &info, mpiAtomStructType );
382			}
383			}
384			#endif // is_mpi
385
386
387			// initialize the atoms
388
389			Atom* thisAtom;
390
391			for( i=0; i<nAtoms; i++ ){
392
393			currentAtomType = headAtomType->find( the_atoms[i]->getType() );
394			if( currentAtomType == NULL ){
395			sprintf( painCave.errMsg,
396			"AtomType error, %s not found in force file.\n",
397			the_atoms[i]->getType() );
398			painCave.isFatal = 1;
399			simError();
400			}
401
402			the_atoms[i]->setMass( currentAtomType->mass );
403			the_atoms[i]->setEpslon( currentAtomType->epslon );
404			the_atoms[i]->setSigma( currentAtomType->sigma );
405			the_atoms[i]->setLJ();
406			}
407
408
409			// clean up the memory
410
411			delete headAtomType;
412
413			#ifdef IS_MPI
414			sprintf( checkPointMsg, "LJ_FF atoms initialized succesfully" );
415			MPIcheckPoint();
416			#endif // is_mpi
417
418			}
419
420			void LJ_FF::initializeBonds( bond_pair* the_bonds ){
421
422			if( entry_plug->n_bonds ){
423			sprintf( painCave.errMsg,
424			"LJ_FF does not support bonds.\n" );
425			painCave.isFatal = 1;
426			simError();
427			}
428			#ifdef IS_MPI
429			MPIcheckPoint();
430			#endif // is_mpi
431
432			}
433
434			void LJ_FF::initializeBends( bend_set* the_bends ){
435
436			if( entry_plug->n_bends ){
437			sprintf( painCave.errMsg,
438			"LJ_FF does not support bends.\n" );
439			painCave.isFatal = 1;
440			simError();
441			}
442			#ifdef IS_MPI
443			MPIcheckPoint();
444			#endif // is_mpi
445
446			}
447
448			void LJ_FF::initializeTorsions( torsion_set* the_torsions ){
449
450			if( entry_plug->n_torsions ){
451			sprintf( painCave.errMsg,
452			"LJ_FF does not support torsions.\n" );
453			painCave.isFatal = 1;
454			simError();
455			}
456			#ifdef IS_MPI
457			MPIcheckPoint();
458			#endif // is_mpi
459
460			}