ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/LJ_FF.cpp
Revision: 348
Committed: Fri Mar 14 21:33:10 2003 UTC (21 years, 5 months ago) by mmeineke
File size: 12342 byte(s)
Log Message:
mostly compiles. interface twixt c and fortran is broken. (c needs to be brought up to date with fortran.)

File Contents

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