8 |
|
#include "UseTheForce/ForceFields.hpp" |
9 |
|
#include "primitives/SRI.hpp" |
10 |
|
#include "utils/simError.h" |
11 |
+ |
#include "types/DirectionalAtomType.hpp" |
12 |
+ |
#include "UseTheForce/DarkSide/lj_interface.h" |
13 |
+ |
#include "UseTheForce/DarkSide/dipole_interface.h" |
14 |
|
#include "UseTheForce/DarkSide/sticky_interface.h" |
12 |
– |
#include "UseTheForce/DarkSide/atype_interface.h" |
13 |
– |
|
14 |
– |
//#include "UseTheForce/fortranWrappers.hpp" |
15 |
– |
|
15 |
|
|
16 |
|
#ifdef IS_MPI |
17 |
|
#include "UseTheForce/mpiForceField.h" |
109 |
|
} |
110 |
|
~LinkedAtomType(){ if( next != NULL ) delete next; } |
111 |
|
|
112 |
< |
LinkedAtomType* find(char* key){ |
112 |
> |
LinkedAtomType* find(const char* key){ |
113 |
|
if( !strcmp(name, key) ) return this; |
114 |
|
if( next != NULL ) return next->find(key); |
115 |
|
return NULL; |
207 |
|
} |
208 |
|
~LinkedBondType(){ if( next != NULL ) delete next; } |
209 |
|
|
210 |
< |
LinkedBondType* find(char* key1, char* key2){ |
210 |
> |
LinkedBondType* find(const char* key1, const char* key2){ |
211 |
|
if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this; |
212 |
|
if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this; |
213 |
|
if( next != NULL ) return next->find(key1, key2); |
277 |
|
} |
278 |
|
~LinkedBendType(){ if( next != NULL ) delete next; } |
279 |
|
|
280 |
< |
LinkedBendType* find( char* key1, char* key2, char* key3 ){ |
280 |
> |
LinkedBendType* find(const char* key1, const char* key2, const char* key3 ){ |
281 |
|
if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) |
282 |
|
&& !strcmp( nameC, key3 ) ) return this; |
283 |
|
if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 ) |
355 |
|
} |
356 |
|
~LinkedTorsionType(){ if( next != NULL ) delete next; } |
357 |
|
|
358 |
< |
LinkedTorsionType* find( char* key1, char* key2, char* key3, char* key4 ){ |
358 |
> |
LinkedTorsionType* find(const char* key1, const char* key2, const char* key3, const char* key4 ){ |
359 |
|
|
360 |
|
|
361 |
|
|
454 |
|
|
455 |
|
DUFF::DUFF(){ |
456 |
|
|
457 |
< |
char fileName[200]; |
458 |
< |
char* ffPath_env = "FORCE_PARAM_PATH"; |
460 |
< |
char* ffPath; |
461 |
< |
char temp[200]; |
457 |
> |
string fileName; |
458 |
> |
string tempString; |
459 |
|
|
460 |
|
headAtomType = NULL; |
461 |
|
currentAtomType = NULL; |
565 |
|
|
566 |
|
// generate the force file name |
567 |
|
|
568 |
< |
strcpy( fileName, "DUFF.frc" ); |
568 |
> |
fileName = "DUFF.frc"; |
569 |
|
// fprintf( stderr,"Trying to open %s\n", fileName ); |
570 |
|
|
571 |
|
// attempt to open the file in the current directory first. |
572 |
|
|
573 |
< |
frcFile = fopen( fileName, "r" ); |
573 |
> |
frcFile = fopen( fileName.c_str(), "r" ); |
574 |
|
|
575 |
|
if( frcFile == NULL ){ |
576 |
|
|
577 |
|
// next see if the force path enviorment variable is set |
578 |
< |
|
579 |
< |
ffPath = getenv( ffPath_env ); |
580 |
< |
if( ffPath == NULL ) { |
581 |
< |
STR_DEFINE(ffPath, FRC_PATH ); |
582 |
< |
} |
578 |
> |
|
579 |
> |
tempString = ffPath + "/" + fileName; |
580 |
> |
fileName = tempString; |
581 |
> |
|
582 |
> |
frcFile = fopen( fileName.c_str(), "r" ); |
583 |
|
|
587 |
– |
|
588 |
– |
strcpy( temp, ffPath ); |
589 |
– |
strcat( temp, "/" ); |
590 |
– |
strcat( temp, fileName ); |
591 |
– |
strcpy( fileName, temp ); |
592 |
– |
|
593 |
– |
frcFile = fopen( fileName, "r" ); |
594 |
– |
|
584 |
|
if( frcFile == NULL ){ |
585 |
|
|
586 |
|
sprintf( painCave.errMsg, |
588 |
|
"\t%s\n" |
589 |
|
"\tHave you tried setting the FORCE_PARAM_PATH environment " |
590 |
|
"variable?\n", |
591 |
< |
fileName ); |
591 |
> |
fileName.c_str() ); |
592 |
|
painCave.severity = OOPSE_ERROR; |
593 |
|
painCave.isFatal = 1; |
594 |
|
simError(); |
642 |
|
} |
643 |
|
|
644 |
|
|
645 |
< |
void DUFF::initForceField( int ljMixRule ){ |
645 |
> |
void DUFF::initForceField(){ |
646 |
|
|
647 |
< |
initFortran( ljMixRule, entry_plug->useReactionField ); |
647 |
> |
initFortran( entry_plug->useReactionField ); |
648 |
|
} |
649 |
|
|
661 |
– |
double DUFF::getAtomTypeMass (char* atomType) { |
650 |
|
|
663 |
– |
currentAtomType = headAtomType->find( atomType ); |
664 |
– |
if( currentAtomType == NULL ){ |
665 |
– |
sprintf( painCave.errMsg, |
666 |
– |
"AtomType error, %s not found in force file.\n", |
667 |
– |
atomType ); |
668 |
– |
painCave.isFatal = 1; |
669 |
– |
simError(); |
670 |
– |
} |
671 |
– |
|
672 |
– |
return currentAtomType->mass; |
673 |
– |
} |
674 |
– |
|
651 |
|
void DUFF::readParams( void ){ |
652 |
|
|
653 |
< |
int identNum; |
653 |
> |
int identNum, isError; |
654 |
|
|
655 |
|
atomStruct atomInfo; |
656 |
|
bondStruct bondInfo; |
657 |
|
bendStruct bendInfo; |
658 |
|
torsionStruct torsionInfo; |
659 |
+ |
|
660 |
+ |
AtomType* at; |
661 |
|
|
662 |
|
bigSigma = 0.0; |
663 |
|
|
758 |
|
} |
759 |
|
|
760 |
|
#endif // is_mpi |
761 |
< |
|
784 |
< |
|
785 |
< |
|
786 |
< |
// call new A_types in fortran |
787 |
< |
|
788 |
< |
int isError; |
789 |
< |
|
761 |
> |
|
762 |
|
// dummy variables |
763 |
< |
|
792 |
< |
int isGB = 0; |
793 |
< |
int isLJ = 1; |
794 |
< |
int isEAM =0; |
795 |
< |
int isCharge = 0; |
796 |
< |
double charge=0.0; |
797 |
< |
|
763 |
> |
|
764 |
|
currentAtomType = headAtomType->next;; |
765 |
< |
while( currentAtomType != NULL ){ |
766 |
< |
|
767 |
< |
if(currentAtomType->isDipole) entry_plug->useDipoles = 1; |
768 |
< |
if(currentAtomType->isSSD) { |
769 |
< |
entry_plug->useSticky = 1; |
770 |
< |
set_sticky_params( &(currentAtomType->w0), &(currentAtomType->v0), |
771 |
< |
&(currentAtomType->v0p), |
772 |
< |
&(currentAtomType->rl), &(currentAtomType->ru), |
773 |
< |
&(currentAtomType->rlp), &(currentAtomType->rup)); |
765 |
> |
while( currentAtomType != NULL ){ |
766 |
> |
|
767 |
> |
if( currentAtomType->name[0] != '\0' ){ |
768 |
> |
|
769 |
> |
if (currentAtomType->isSSD || currentAtomType->isDipole) |
770 |
> |
at = new DirectionalAtomType(); |
771 |
> |
else |
772 |
> |
at = new AtomType(); |
773 |
> |
|
774 |
> |
if (currentAtomType->isSSD) { |
775 |
> |
((DirectionalAtomType*)at)->setSticky(); |
776 |
> |
entry_plug->useSticky = 1; |
777 |
> |
} |
778 |
> |
|
779 |
> |
if (currentAtomType->isDipole) { |
780 |
> |
((DirectionalAtomType*)at)->setDipole(); |
781 |
> |
entry_plug->useDipoles = 1; |
782 |
> |
} |
783 |
> |
|
784 |
> |
at->setIdent(currentAtomType->ident); |
785 |
> |
at->setName(currentAtomType->name); |
786 |
> |
at->setLennardJones(); |
787 |
> |
at->complete(); |
788 |
|
} |
789 |
+ |
currentAtomType = currentAtomType->next; |
790 |
+ |
} |
791 |
+ |
|
792 |
+ |
currentAtomType = headAtomType->next;; |
793 |
+ |
while( currentAtomType != NULL ){ |
794 |
|
|
795 |
|
if( currentAtomType->name[0] != '\0' ){ |
796 |
|
isError = 0; |
797 |
< |
makeAtype( &(currentAtomType->ident), |
798 |
< |
&isLJ, |
814 |
< |
&(currentAtomType->isSSD), |
815 |
< |
&(currentAtomType->isDipole), |
816 |
< |
&isGB, |
817 |
< |
&isEAM, |
818 |
< |
&isCharge, |
819 |
< |
&(currentAtomType->epslon), |
820 |
< |
&(currentAtomType->sigma), |
821 |
< |
&charge, |
822 |
< |
&(currentAtomType->dipole), |
823 |
< |
&isError ); |
797 |
> |
newLJtype( &(currentAtomType->ident), &(currentAtomType->sigma), |
798 |
> |
&(currentAtomType->epslon), &isError); |
799 |
|
if( isError ){ |
800 |
< |
sprintf( painCave.errMsg, |
801 |
< |
"Error initializing the \"%s\" atom type in fortran\n", |
802 |
< |
currentAtomType->name ); |
803 |
< |
painCave.isFatal = 1; |
804 |
< |
simError(); |
800 |
> |
sprintf( painCave.errMsg, |
801 |
> |
"Error initializing the \"%s\" LJ type in fortran\n", |
802 |
> |
currentAtomType->name ); |
803 |
> |
painCave.isFatal = 1; |
804 |
> |
simError(); |
805 |
|
} |
806 |
+ |
|
807 |
+ |
if (currentAtomType->isDipole) { |
808 |
+ |
newDipoleType(&(currentAtomType->ident), &(currentAtomType->dipole), |
809 |
+ |
&isError); |
810 |
+ |
if( isError ){ |
811 |
+ |
sprintf( painCave.errMsg, |
812 |
+ |
"Error initializing the \"%s\" dipole type in fortran\n", |
813 |
+ |
currentAtomType->name ); |
814 |
+ |
painCave.isFatal = 1; |
815 |
+ |
simError(); |
816 |
+ |
} |
817 |
+ |
} |
818 |
+ |
|
819 |
+ |
if(currentAtomType->isSSD) { |
820 |
+ |
makeStickyType( &(currentAtomType->w0), &(currentAtomType->v0), |
821 |
+ |
&(currentAtomType->v0p), |
822 |
+ |
&(currentAtomType->rl), &(currentAtomType->ru), |
823 |
+ |
&(currentAtomType->rlp), &(currentAtomType->rup)); |
824 |
+ |
} |
825 |
+ |
|
826 |
|
} |
827 |
|
currentAtomType = currentAtomType->next; |
828 |
|
} |
1086 |
|
|
1087 |
|
#endif // is_mpi |
1088 |
|
|
1089 |
< |
entry_plug->useLJ = 1; |
1089 |
> |
entry_plug->useLennardJones = 1; |
1090 |
|
} |
1091 |
|
|
1092 |
|
|
1134 |
|
|
1135 |
|
for(int i=0; i<nAtoms; i++ ){ |
1136 |
|
|
1137 |
< |
currentAtomType = headAtomType->find( the_atoms[i]->getType() ); |
1137 |
> |
currentAtomType = headAtomType->find(the_atoms[i]->getType().c_str() ); |
1138 |
|
if( currentAtomType == NULL ){ |
1139 |
|
sprintf( painCave.errMsg, |
1140 |
|
"AtomType error, %s not found in force file.\n", |
1141 |
< |
the_atoms[i]->getType() ); |
1141 |
> |
the_atoms[i]->getType().c_str() ); |
1142 |
|
painCave.isFatal = 1; |
1143 |
|
simError(); |
1144 |
|
} |
1160 |
|
|
1161 |
|
dAtom->setJ( ji ); |
1162 |
|
|
1163 |
< |
if(!strcmp("SSD",the_atoms[i]->getType())){ |
1163 |
> |
if(!strcmp("SSD",the_atoms[i]->getType().c_str())){ |
1164 |
|
dAtom->setI( waterI ); |
1165 |
|
} |
1166 |
< |
else if(!strcmp("HEAD",the_atoms[i]->getType())){ |
1166 |
> |
else if(!strcmp("HEAD",the_atoms[i]->getType().c_str())){ |
1167 |
|
dAtom->setI( headI ); |
1168 |
|
} |
1169 |
|
else{ |
1170 |
|
sprintf(painCave.errMsg, |
1171 |
|
"AtmType error, %s does not have a moment of inertia set.\n", |
1172 |
< |
the_atoms[i]->getType() ); |
1172 |
> |
the_atoms[i]->getType().c_str() ); |
1173 |
|
painCave.isFatal = 1; |
1174 |
|
simError(); |
1175 |
|
} |
1201 |
|
void DUFF::initializeBonds( int nBonds, Bond** bondArray, |
1202 |
|
bond_pair* the_bonds ){ |
1203 |
|
int i,a,b; |
1204 |
< |
char* atomA; |
1210 |
< |
char* atomB; |
1211 |
< |
|
1204 |
> |
|
1205 |
|
Atom** the_atoms; |
1206 |
|
the_atoms = entry_plug->atoms; |
1207 |
|
|
1213 |
|
a = the_bonds[i].a; |
1214 |
|
b = the_bonds[i].b; |
1215 |
|
|
1216 |
< |
atomA = the_atoms[a]->getType(); |
1224 |
< |
atomB = the_atoms[b]->getType(); |
1225 |
< |
currentBondType = headBondType->find( atomA, atomB ); |
1216 |
> |
currentBondType = headBondType->find( the_atoms[a]->getType().c_str(), the_atoms[b]->getType().c_str() ); |
1217 |
|
if( currentBondType == NULL ){ |
1218 |
|
sprintf( painCave.errMsg, |
1219 |
|
"BondType error, %s - %s not found in force file.\n", |
1220 |
< |
atomA, atomB ); |
1220 |
> |
the_atoms[a]->getType().c_str(), the_atoms[b]->getType().c_str() ); |
1221 |
|
painCave.isFatal = 1; |
1222 |
|
simError(); |
1223 |
|
} |
1257 |
|
the_atoms = entry_plug->atoms; |
1258 |
|
|
1259 |
|
int i, a, b, c; |
1260 |
< |
char* atomA; |
1261 |
< |
char* atomB; |
1262 |
< |
char* atomC; |
1260 |
> |
string atomA; |
1261 |
> |
string atomB; |
1262 |
> |
string atomC; |
1263 |
|
|
1264 |
|
// initialize the Bends |
1265 |
|
|
1275 |
|
if( the_bends[i].isGhost ) atomC = "GHOST"; |
1276 |
|
else atomC = the_atoms[c]->getType(); |
1277 |
|
|
1278 |
< |
currentBendType = headBendType->find( atomA, atomB, atomC ); |
1278 |
> |
currentBendType = headBendType->find( atomA.c_str(), atomB.c_str(), atomC.c_str()); |
1279 |
|
if( currentBendType == NULL ){ |
1280 |
|
sprintf( painCave.errMsg, "BendType error, %s - %s - %s not found" |
1281 |
|
" in force file.\n", |
1282 |
< |
atomA, atomB, atomC ); |
1282 |
> |
atomA.c_str(), atomB.c_str(), atomC.c_str() ); |
1283 |
|
painCave.isFatal = 1; |
1284 |
|
simError(); |
1285 |
|
} |
1302 |
|
" --> central atom is not " |
1303 |
|
"correctly identified with the " |
1304 |
|
"\"ghostVectorSource = \" tag.\n", |
1305 |
< |
atomA, atomB, atomC ); |
1305 |
> |
atomA.c_str(), atomB.c_str(), atomC.c_str() ); |
1306 |
|
painCave.isFatal = 1; |
1307 |
|
simError(); |
1308 |
|
} |
1334 |
|
torsion_set* the_torsions ){ |
1335 |
|
|
1336 |
|
int i, a, b, c, d; |
1346 |
– |
char* atomA; |
1347 |
– |
char* atomB; |
1348 |
– |
char* atomC; |
1349 |
– |
char* atomD; |
1337 |
|
|
1338 |
|
CubicTorsion* cTors; |
1339 |
|
Atom** the_atoms; |
1348 |
|
c = the_torsions[i].c; |
1349 |
|
d = the_torsions[i].d; |
1350 |
|
|
1351 |
< |
atomA = the_atoms[a]->getType(); |
1352 |
< |
atomB = the_atoms[b]->getType(); |
1353 |
< |
atomC = the_atoms[c]->getType(); |
1354 |
< |
atomD = the_atoms[d]->getType(); |
1368 |
< |
currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD ); |
1351 |
> |
currentTorsionType = headTorsionType->find( the_atoms[a]->getType().c_str(), |
1352 |
> |
the_atoms[b]->getType().c_str(), |
1353 |
> |
the_atoms[c]->getType().c_str(), |
1354 |
> |
the_atoms[d]->getType().c_str() ); |
1355 |
|
if( currentTorsionType == NULL ){ |
1356 |
|
sprintf( painCave.errMsg, |
1357 |
|
"TorsionType error, %s - %s - %s - %s not found" |
1358 |
|
" in force file.\n", |
1359 |
< |
atomA, atomB, atomC, atomD ); |
1359 |
> |
the_atoms[a]->getType().c_str(), |
1360 |
> |
the_atoms[b]->getType().c_str(), |
1361 |
> |
the_atoms[c]->getType().c_str(), |
1362 |
> |
the_atoms[d]->getType().c_str() ); |
1363 |
|
painCave.isFatal = 1; |
1364 |
|
simError(); |
1365 |
|
} |