ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/LJ_FF.cpp
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

File Contents

# Content
1 #include <cstdlib>
2 #include <cstdio>
3 #include <cstring>
4
5 #include <iostream>
6 using namespace std;
7
8 #ifdef IS_MPI
9 #include <mpi.h>
10 #include <mpi++.h>
11 #endif //is_mpi
12
13 #include "ForceFields.hpp"
14 #include "SRI.hpp"
15 #include "simError.h"
16
17 #include <fortranWrappers.hpp>
18
19 #ifdef IS_MPI
20 #include "mpiForceField.h"
21 #endif // is_mpi
22
23
24
25 namespace {
26
27 // 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
39 int parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
40
41 #ifdef IS_MPI
42
43 MPI_Datatype mpiAtomStructType;
44
45 #endif
46 }
47
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 wrapMeFF( this );
64
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 void LJ_FF::doForces( int calcPot ){
465
466 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 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 int parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
593
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 }