12 |
|
#include "UseTheForce/ForceFields.hpp" |
13 |
|
#include "primitives/SRI.hpp" |
14 |
|
#include "utils/simError.h" |
15 |
< |
|
16 |
< |
|
17 |
< |
#include "UseTheForce/DarkSide/atype_interface.h" |
15 |
> |
#include "types/AtomType.hpp" |
16 |
> |
#include "types/DirectionalAtomType.hpp" |
17 |
> |
#include "UseTheForce/DarkSide/lj_interface.h" |
18 |
> |
#include "UseTheForce/DarkSide/charge_interface.h" |
19 |
> |
#include "UseTheForce/DarkSide/dipole_interface.h" |
20 |
|
#include "UseTheForce/DarkSide/sticky_interface.h" |
21 |
|
|
22 |
|
#ifdef IS_MPI |
247 |
|
|
248 |
|
WATER::WATER(){ |
249 |
|
|
250 |
< |
char fileName[200]; |
251 |
< |
char* ffPath_env = "FORCE_PARAM_PATH"; |
250 |
< |
char* ffPath; |
251 |
< |
char temp[200]; |
250 |
> |
string fileName; |
251 |
> |
string tempString; |
252 |
|
|
253 |
|
headAtomType = NULL; |
254 |
|
currentAtomType = NULL; |
308 |
|
if( worldRank == 0 ){ |
309 |
|
#endif |
310 |
|
|
311 |
< |
// generate the force file name |
312 |
< |
|
313 |
< |
strcpy( fileName, "WATER.frc" ); |
311 |
> |
// generate the force file name |
312 |
> |
|
313 |
> |
fileName = "WATER.frc"; |
314 |
> |
|
315 |
|
// fprintf( stderr,"Trying to open %s\n", fileName ); |
316 |
|
|
317 |
|
// attempt to open the file in the current directory first. |
318 |
|
|
319 |
< |
frcFile = fopen( fileName, "r" ); |
319 |
> |
frcFile = fopen( fileName.c_str(), "r" ); |
320 |
|
|
321 |
|
if( frcFile == NULL ){ |
322 |
|
|
323 |
< |
// next see if the force path enviorment variable is set |
323 |
> |
tempString = ffPath + "/" + fileName; |
324 |
> |
fileName = tempString; |
325 |
|
|
326 |
< |
ffPath = getenv( ffPath_env ); |
325 |
< |
if( ffPath == NULL ) { |
326 |
< |
STR_DEFINE(ffPath, FRC_PATH ); |
327 |
< |
} |
326 |
> |
frcFile = fopen( fileName.c_str(), "r" ); |
327 |
|
|
329 |
– |
|
330 |
– |
strcpy( temp, ffPath ); |
331 |
– |
strcat( temp, "/" ); |
332 |
– |
strcat( temp, fileName ); |
333 |
– |
strcpy( fileName, temp ); |
334 |
– |
|
335 |
– |
frcFile = fopen( fileName, "r" ); |
336 |
– |
|
328 |
|
if( frcFile == NULL ){ |
329 |
|
|
330 |
|
sprintf( painCave.errMsg, |
332 |
|
"\t%s\n" |
333 |
|
"\tHave you tried setting the FORCE_PARAM_PATH environment " |
334 |
|
"variable?\n", |
335 |
< |
fileName ); |
335 |
> |
fileName.c_str() ); |
336 |
|
painCave.severity = OOPSE_ERROR; |
337 |
|
painCave.isFatal = 1; |
338 |
|
simError(); |
390 |
|
|
391 |
|
void WATER::readParams( void ){ |
392 |
|
|
393 |
< |
int identNum; |
393 |
> |
int identNum, isError; |
394 |
|
int tempDirect0, tempDirect1; |
395 |
|
|
396 |
|
atomStruct atomInfo; |
397 |
|
directionalStruct directionalInfo; |
398 |
|
fpos_t *atomPos; |
399 |
+ |
AtomType* at; |
400 |
|
|
401 |
|
atomInfo.last = 1; // initialize last to have the last set. |
402 |
|
directionalInfo.last = 1; // if things go well, last will be set to 0 |
536 |
|
|
537 |
|
#endif // is_mpi |
538 |
|
|
547 |
– |
// call new A_types in fortran |
548 |
– |
|
549 |
– |
int isError; |
550 |
– |
|
539 |
|
// dummy variables |
552 |
– |
int isGB = 0; |
553 |
– |
int isEAM = 0; |
554 |
– |
int notDipole = 0; |
555 |
– |
int notSSD = 0; |
556 |
– |
double noDipMoment = 0.0; |
540 |
|
|
558 |
– |
|
541 |
|
currentAtomType = headAtomType->next; |
542 |
|
currentDirectionalType = headDirectionalType->next; |
543 |
|
|
544 |
|
while( currentAtomType != NULL ){ |
545 |
< |
|
546 |
< |
if( currentAtomType->isLJ ) entry_plug->useLJ = 1; |
547 |
< |
if( currentAtomType->isCharge ) entry_plug->useCharges = 1; |
548 |
< |
if( currentAtomType->isDirectional ){ |
549 |
< |
// only directional atoms can have dipoles or be sticky |
550 |
< |
if ( currentDirectionalType->isDipole ) entry_plug->useDipoles = 1; |
551 |
< |
if ( currentDirectionalType->isSticky ) { |
552 |
< |
entry_plug->useSticky = 1; |
571 |
< |
makeStickyType( &(currentDirectionalType->w0), |
572 |
< |
&(currentDirectionalType->v0), |
573 |
< |
&(currentDirectionalType->v0p), |
574 |
< |
&(currentDirectionalType->rl), |
575 |
< |
&(currentDirectionalType->ru), |
576 |
< |
&(currentDirectionalType->rlp), |
577 |
< |
&(currentDirectionalType->rup)); |
545 |
> |
if( currentAtomType->name[0] != '\0' ){ |
546 |
> |
if (currentAtomType->isDirectional) |
547 |
> |
at = new DirectionalAtomType(); |
548 |
> |
else |
549 |
> |
at = new AtomType(); |
550 |
> |
|
551 |
> |
if (currentAtomType->isLJ) { |
552 |
> |
at->setLennardJones(); |
553 |
|
} |
554 |
< |
if( currentAtomType->name[0] != '\0' ){ |
555 |
< |
isError = 0; |
556 |
< |
makeAtype( &(currentAtomType->ident), |
557 |
< |
&(currentAtomType->isLJ), |
558 |
< |
&(currentDirectionalType->isSticky), |
559 |
< |
&(currentDirectionalType->isDipole), |
560 |
< |
&isGB, |
561 |
< |
&isEAM, |
562 |
< |
&(currentAtomType->isCharge), |
588 |
< |
&(currentAtomType->epslon), |
589 |
< |
&(currentAtomType->sigma), |
590 |
< |
&(currentAtomType->charge), |
591 |
< |
&(currentDirectionalType->dipole), |
592 |
< |
&isError ); |
593 |
< |
if( isError ){ |
594 |
< |
sprintf( painCave.errMsg, |
595 |
< |
"Error initializing the \"%s\" atom type in fortran\n", |
596 |
< |
currentAtomType->name ); |
597 |
< |
painCave.isFatal = 1; |
598 |
< |
simError(); |
554 |
> |
|
555 |
> |
if (currentAtomType->isCharge) { |
556 |
> |
at->setCharge(); |
557 |
> |
} |
558 |
> |
|
559 |
> |
if (currentAtomType->isDirectional) { |
560 |
> |
if (currentDirectionalType->isDipole) { |
561 |
> |
((DirectionalAtomType*)at)->setDipole(); |
562 |
> |
entry_plug->useDipoles = 1; |
563 |
|
} |
564 |
+ |
|
565 |
+ |
if (currentDirectionalType->isSticky) { |
566 |
+ |
((DirectionalAtomType*)at)->setSticky(); |
567 |
+ |
entry_plug->useSticky = 1; |
568 |
+ |
} |
569 |
|
} |
570 |
< |
currentDirectionalType->next; |
570 |
> |
|
571 |
> |
at->setIdent(currentAtomType->ident); |
572 |
> |
at->setName(currentAtomType->name); |
573 |
> |
at->complete(); |
574 |
|
} |
575 |
+ |
currentAtomType = currentAtomType->next; |
576 |
+ |
} |
577 |
|
|
578 |
< |
else { |
579 |
< |
// use all dummy variables if this is not a directional atom |
580 |
< |
if( currentAtomType->name[0] != '\0' ){ |
581 |
< |
isError = 0; |
582 |
< |
makeAtype( &(currentAtomType->ident), |
583 |
< |
&(currentAtomType->isLJ), |
584 |
< |
¬SSD, |
585 |
< |
¬Dipole, |
586 |
< |
&isGB, |
587 |
< |
&isEAM, |
588 |
< |
&(currentAtomType->isCharge), |
589 |
< |
&(currentAtomType->epslon), |
590 |
< |
&(currentAtomType->sigma), |
591 |
< |
&(currentAtomType->charge), |
592 |
< |
&noDipMoment, |
593 |
< |
&isError ); |
578 |
> |
currentAtomType = headAtomType->next; |
579 |
> |
currentDirectionalType = headDirectionalType->next; |
580 |
> |
|
581 |
> |
while( currentAtomType != NULL ){ |
582 |
> |
|
583 |
> |
if( currentAtomType->isLJ ){ |
584 |
> |
isError = 0; |
585 |
> |
newLJtype( &(currentAtomType->ident), &(currentAtomType->sigma), |
586 |
> |
&(currentAtomType->epslon), &isError); |
587 |
> |
if( isError ){ |
588 |
> |
sprintf( painCave.errMsg, |
589 |
> |
"Error initializing the \"%s\" LJ type in fortran\n", |
590 |
> |
currentAtomType->name ); |
591 |
> |
painCave.isFatal = 1; |
592 |
> |
simError(); |
593 |
> |
} |
594 |
> |
} |
595 |
> |
|
596 |
> |
if (currentAtomType->isCharge) { |
597 |
> |
newChargeType(&(currentAtomType->ident), &(currentAtomType->charge), |
598 |
> |
&isError); |
599 |
> |
if( isError ){ |
600 |
> |
sprintf( painCave.errMsg, |
601 |
> |
"Error initializing the \"%s\" charge type in fortran\n", |
602 |
> |
currentAtomType->name ); |
603 |
> |
painCave.isFatal = 1; |
604 |
> |
simError(); |
605 |
> |
} |
606 |
> |
} |
607 |
> |
|
608 |
> |
if (currentAtomType->isDirectional){ |
609 |
> |
if (currentDirectionalType->isDipole) { |
610 |
> |
newDipoleType(&(currentAtomType->ident), |
611 |
> |
&(currentDirectionalType->dipole), |
612 |
> |
&isError); |
613 |
|
if( isError ){ |
614 |
|
sprintf( painCave.errMsg, |
615 |
< |
"Error initializing the \"%s\" atom type in fortran\n", |
616 |
< |
currentAtomType->name ); |
615 |
> |
"Error initializing the \"%s\" dipole type in fortran\n", |
616 |
> |
currentDirectionalType->name ); |
617 |
|
painCave.isFatal = 1; |
618 |
|
simError(); |
619 |
|
} |
620 |
|
} |
621 |
+ |
|
622 |
+ |
if(currentDirectionalType->isSticky) { |
623 |
+ |
makeStickyType( &(currentDirectionalType->w0), |
624 |
+ |
&(currentDirectionalType->v0), |
625 |
+ |
&(currentDirectionalType->v0p), |
626 |
+ |
&(currentDirectionalType->rl), |
627 |
+ |
&(currentDirectionalType->ru), |
628 |
+ |
&(currentDirectionalType->rlp), |
629 |
+ |
&(currentDirectionalType->rup)); |
630 |
+ |
} |
631 |
|
} |
632 |
|
currentAtomType = currentAtomType->next; |
633 |
|
} |
634 |
< |
|
634 |
> |
|
635 |
|
#ifdef IS_MPI |
636 |
|
sprintf( checkPointMsg, |
637 |
|
"WATER atom and directional structures successfully" |
638 |
|
"sent to fortran\n" ); |
639 |
|
MPIcheckPoint(); |
640 |
|
#endif // is_mpi |
641 |
< |
|
641 |
> |
|
642 |
|
} |
643 |
|
|
644 |
|
|
1077 |
|
return 1; |
1078 |
|
} |
1079 |
|
else return 0; |
1077 |
– |
} |
1078 |
– |
double WATER::getAtomTypeMass (char* atomType) { |
1079 |
– |
|
1080 |
– |
currentAtomType = headAtomType->find( atomType ); |
1081 |
– |
if( currentAtomType == NULL ){ |
1082 |
– |
sprintf( painCave.errMsg, |
1083 |
– |
"AtomType error, %s not found in force file.\n", |
1084 |
– |
atomType ); |
1085 |
– |
painCave.isFatal = 1; |
1086 |
– |
simError(); |
1087 |
– |
} |
1088 |
– |
|
1089 |
– |
return currentAtomType->mass; |
1080 |
|
} |