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, 6 months ago) by mmeineke
File size: 17469 byte(s)
Log Message:
overhauled TraPPE_Ex forcfield in C++

File Contents

# User Rev Content
1 mmeineke 287 #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 mmeineke 291 SSD_FF* currentSSDwrap;
78 mmeineke 287
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     }