ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/ssdFF++.cpp
Revision: 291
Committed: Wed Mar 5 20:35:54 2003 UTC (21 years, 4 months ago) by mmeineke
File size: 17469 byte(s)
Log Message:
overhauled TraPPE_Ex forcfield in C++

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