10 |
|
|
11 |
|
#ifdef IS_MPI |
12 |
|
#include "mpiBASS.h" |
13 |
+ |
#include "mpiSimulation.hpp" |
14 |
|
#include "bassDiag.hpp" |
15 |
|
#endif |
16 |
|
|
74 |
|
|
75 |
|
MakeStamps *the_stamps; |
76 |
|
Globals* the_globals; |
77 |
< |
int i; |
77 |
> |
int i, j; |
78 |
|
|
79 |
|
// get the stamps and globals; |
80 |
|
the_stamps = stamps; |
93 |
|
if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF(); |
94 |
|
else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF(); |
95 |
|
else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF(); |
96 |
+ |
else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF(); |
97 |
|
else{ |
98 |
|
sprintf( painCave.errMsg, |
99 |
|
"SimSetup Error. Unrecognized force field -> %s\n", |
189 |
|
if( currentStamp == NULL ){ |
190 |
|
sprintf( painCave.errMsg, |
191 |
|
"SimSetup error: Component \"%s\" was not found in the " |
192 |
< |
"list of declared molecules\n" |
192 |
> |
"list of declared molecules\n", |
193 |
|
id ); |
194 |
|
painCave.isFatal = 1; |
195 |
|
simError(); |
215 |
|
tot_bends = 0; |
216 |
|
tot_torsions = 0; |
217 |
|
for( i=0; i<n_components; i++ ){ |
218 |
< |
|
219 |
< |
tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); |
220 |
< |
tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); |
221 |
< |
tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); |
218 |
> |
|
219 |
> |
tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); |
220 |
> |
tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); |
221 |
> |
tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); |
222 |
|
tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); |
223 |
|
} |
224 |
|
|
236 |
|
|
237 |
|
// divide the molecules among processors here. |
238 |
|
|
239 |
< |
new mpiSimulation( simnfo ); |
239 |
> |
mpiSim = new mpiSimulation( simnfo ); |
240 |
|
|
241 |
< |
simnfo->mpiSim->divideLabor( n_components, comp_stamps, components_nmol ); |
241 |
> |
|
242 |
> |
|
243 |
> |
globalIndex = mpiSim->divideLabor(); |
244 |
> |
|
245 |
> |
|
246 |
> |
|
247 |
> |
// set up the local variables |
248 |
> |
|
249 |
> |
int localMol, allMol; |
250 |
> |
int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; |
251 |
> |
|
252 |
> |
allMol = 0; |
253 |
> |
localMol = 0; |
254 |
> |
local_atoms = 0; |
255 |
> |
local_bonds = 0; |
256 |
> |
local_bends = 0; |
257 |
> |
local_torsions = 0; |
258 |
> |
for( i=0; i<n_components; i++ ){ |
259 |
> |
|
260 |
> |
for( j=0; j<components_nmol[i]; j++ ){ |
261 |
> |
|
262 |
> |
if( mpiSim->getMyMolStart() <= allMol && |
263 |
> |
allMol <= mpiSim->getMyMolEnd() ){ |
264 |
> |
|
265 |
> |
local_atoms += comp_stamps[i]->getNAtoms(); |
266 |
> |
local_bonds += comp_stamps[i]->getNBonds(); |
267 |
> |
local_bends += comp_stamps[i]->getNBends(); |
268 |
> |
local_torsions += comp_stamps[i]->getNTorsions(); |
269 |
> |
localMol++; |
270 |
> |
} |
271 |
> |
allMol++; |
272 |
> |
} |
273 |
> |
} |
274 |
> |
local_SRI = local_bonds + local_bends + local_torsions; |
275 |
> |
|
276 |
|
|
277 |
+ |
simnfo->n_atoms = mpiSim->getMyNlocal(); |
278 |
+ |
|
279 |
+ |
if( local_atoms != simnfo->n_atoms ){ |
280 |
+ |
sprintf( painCave.errMsg, |
281 |
+ |
"SimSetup error: mpiSim's localAtom (%d) and SimSetup's" |
282 |
+ |
" localAtom (%d) are note equal.\n", |
283 |
+ |
simnfo->n_atoms, |
284 |
+ |
local_atoms ); |
285 |
+ |
painCave.isFatal = 1; |
286 |
+ |
simError(); |
287 |
+ |
} |
288 |
+ |
|
289 |
+ |
simnfo->n_bonds = local_bonds; |
290 |
+ |
simnfo->n_bends = local_bends; |
291 |
+ |
simnfo->n_torsions = local_torsions; |
292 |
+ |
simnfo->n_SRI = local_SRI; |
293 |
+ |
simnfo->n_mol = localMol; |
294 |
+ |
|
295 |
+ |
strcpy( checkPointMsg, "Passed nlocal consistency check." ); |
296 |
+ |
MPIcheckPoint(); |
297 |
+ |
|
298 |
+ |
|
299 |
|
#endif // is_mpi |
300 |
|
|
301 |
|
|
319 |
|
simnfo->excludes = the_excludes; |
320 |
|
|
321 |
|
|
264 |
– |
// initialize the arrays |
265 |
– |
|
266 |
– |
the_ff->setSimInfo( simnfo ); |
267 |
– |
|
268 |
– |
makeAtoms(); |
269 |
– |
|
270 |
– |
if( tot_bonds ){ |
271 |
– |
makeBonds(); |
272 |
– |
} |
273 |
– |
|
274 |
– |
if( tot_bends ){ |
275 |
– |
makeBends(); |
276 |
– |
} |
277 |
– |
|
278 |
– |
if( tot_torsions ){ |
279 |
– |
makeTorsions(); |
280 |
– |
} |
281 |
– |
|
282 |
– |
|
322 |
|
// get some of the tricky things that may still be in the globals |
323 |
|
|
324 |
|
if( simnfo->n_dipoles ){ |
391 |
|
#endif // is_mpi |
392 |
|
|
393 |
|
|
394 |
+ |
// initialize the arrays |
395 |
|
|
396 |
< |
// if( the_globals->haveInitialConfig() ){ |
357 |
< |
// InitializeFromFile* fileInit; |
358 |
< |
// fileInit = new InitializeFromFile( the_globals->getInitialConfig() ); |
396 |
> |
the_ff->setSimInfo( simnfo ); |
397 |
|
|
398 |
< |
// fileInit->read_xyz( simnfo ); // default velocities on |
398 |
> |
makeAtoms(); |
399 |
|
|
400 |
< |
// delete fileInit; |
401 |
< |
// } |
402 |
< |
// else{ |
400 |
> |
if( tot_bonds ){ |
401 |
> |
makeBonds(); |
402 |
> |
} |
403 |
|
|
404 |
+ |
if( tot_bends ){ |
405 |
+ |
makeBends(); |
406 |
+ |
} |
407 |
+ |
|
408 |
+ |
if( tot_torsions ){ |
409 |
+ |
makeTorsions(); |
410 |
+ |
} |
411 |
+ |
|
412 |
+ |
|
413 |
+ |
|
414 |
+ |
|
415 |
+ |
|
416 |
+ |
|
417 |
+ |
if( the_globals->haveInitialConfig() ){ |
418 |
+ |
|
419 |
+ |
InitializeFromFile* fileInit; |
420 |
+ |
#ifdef IS_MPI // is_mpi |
421 |
+ |
if( worldRank == 0 ){ |
422 |
+ |
#endif //is_mpi |
423 |
+ |
fileInit = new InitializeFromFile( the_globals->getInitialConfig() ); |
424 |
+ |
#ifdef IS_MPI |
425 |
+ |
}else fileInit = new InitializeFromFile( NULL ); |
426 |
+ |
#endif |
427 |
+ |
fileInit->read_xyz( simnfo ); // default velocities on |
428 |
+ |
|
429 |
+ |
delete fileInit; |
430 |
+ |
} |
431 |
+ |
else{ |
432 |
+ |
|
433 |
|
#ifdef IS_MPI |
434 |
|
|
435 |
|
// no init from bass |
443 |
|
|
444 |
|
initFromBass(); |
445 |
|
|
446 |
< |
#endif // is_mpi |
447 |
< |
|
446 |
> |
|
447 |
> |
#endif |
448 |
> |
} |
449 |
> |
|
450 |
|
#ifdef IS_MPI |
451 |
|
strcpy( checkPointMsg, "Successfully read in the initial configuration" ); |
452 |
|
MPIcheckPoint(); |
456 |
|
|
457 |
|
|
458 |
|
|
459 |
< |
// } |
459 |
> |
|
460 |
|
|
461 |
|
#ifdef IS_MPI |
462 |
|
if( worldRank == 0 ){ |
567 |
|
if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet(); |
568 |
|
|
569 |
|
|
570 |
< |
// make the longe range forces and the integrator |
570 |
> |
// // make the longe range forces and the integrator |
571 |
|
|
572 |
< |
new AllLong( simnfo ); |
572 |
> |
// new AllLong( simnfo ); |
573 |
|
|
574 |
< |
if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo ); |
574 |
> |
if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff ); |
575 |
|
if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo ); |
576 |
|
if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo ); |
577 |
+ |
if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff ); |
578 |
+ |
|
579 |
|
} |
580 |
|
|
581 |
|
void SimSetup::makeAtoms( void ){ |
584 |
|
double ux, uy, uz, uSqr, u; |
585 |
|
AtomStamp* current_atom; |
586 |
|
DirectionalAtom* dAtom; |
587 |
< |
int molIndex, molStart, molEnd, nMemb; |
587 |
> |
int molIndex, molStart, molEnd, nMemb, lMolIndex; |
588 |
|
|
589 |
< |
|
589 |
> |
lMolIndex = 0; |
590 |
|
molIndex = 0; |
591 |
|
index = 0; |
592 |
|
for( i=0; i<n_components; i++ ){ |
593 |
|
|
594 |
|
for( j=0; j<components_nmol[i]; j++ ){ |
595 |
|
|
596 |
< |
molStart = index; |
597 |
< |
nMemb = comp_stamps[i]->getNAtoms(); |
598 |
< |
for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){ |
596 |
> |
#ifdef IS_MPI |
597 |
> |
if( mpiSim->getMyMolStart() <= molIndex && |
598 |
> |
molIndex <= mpiSim->getMyMolEnd() ){ |
599 |
> |
#endif // is_mpi |
600 |
|
|
601 |
< |
current_atom = comp_stamps[i]->getAtom( k ); |
602 |
< |
if( current_atom->haveOrientation() ){ |
603 |
< |
|
604 |
< |
dAtom = new DirectionalAtom(index); |
605 |
< |
simnfo->n_oriented++; |
606 |
< |
the_atoms[index] = dAtom; |
607 |
< |
|
608 |
< |
ux = current_atom->getOrntX(); |
609 |
< |
uy = current_atom->getOrntY(); |
610 |
< |
uz = current_atom->getOrntZ(); |
611 |
< |
|
612 |
< |
uSqr = (ux * ux) + (uy * uy) + (uz * uz); |
613 |
< |
|
614 |
< |
u = sqrt( uSqr ); |
615 |
< |
ux = ux / u; |
616 |
< |
uy = uy / u; |
617 |
< |
uz = uz / u; |
618 |
< |
|
619 |
< |
dAtom->setSUx( ux ); |
620 |
< |
dAtom->setSUy( uy ); |
621 |
< |
dAtom->setSUz( uz ); |
601 |
> |
molStart = index; |
602 |
> |
nMemb = comp_stamps[i]->getNAtoms(); |
603 |
> |
for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){ |
604 |
> |
|
605 |
> |
current_atom = comp_stamps[i]->getAtom( k ); |
606 |
> |
if( current_atom->haveOrientation() ){ |
607 |
> |
|
608 |
> |
dAtom = new DirectionalAtom(index); |
609 |
> |
simnfo->n_oriented++; |
610 |
> |
the_atoms[index] = dAtom; |
611 |
> |
|
612 |
> |
ux = current_atom->getOrntX(); |
613 |
> |
uy = current_atom->getOrntY(); |
614 |
> |
uz = current_atom->getOrntZ(); |
615 |
> |
|
616 |
> |
uSqr = (ux * ux) + (uy * uy) + (uz * uz); |
617 |
> |
|
618 |
> |
u = sqrt( uSqr ); |
619 |
> |
ux = ux / u; |
620 |
> |
uy = uy / u; |
621 |
> |
uz = uz / u; |
622 |
> |
|
623 |
> |
dAtom->setSUx( ux ); |
624 |
> |
dAtom->setSUy( uy ); |
625 |
> |
dAtom->setSUz( uz ); |
626 |
> |
} |
627 |
> |
else{ |
628 |
> |
the_atoms[index] = new GeneralAtom(index); |
629 |
> |
} |
630 |
> |
the_atoms[index]->setType( current_atom->getType() ); |
631 |
> |
the_atoms[index]->setIndex( index ); |
632 |
> |
|
633 |
> |
// increment the index and repeat; |
634 |
> |
index++; |
635 |
|
} |
636 |
< |
else{ |
637 |
< |
the_atoms[index] = new GeneralAtom(index); |
638 |
< |
} |
639 |
< |
the_atoms[index]->setType( current_atom->getType() ); |
640 |
< |
the_atoms[index]->setIndex( index ); |
636 |
> |
|
637 |
> |
molEnd = index -1; |
638 |
> |
the_molecules[lMolIndex].setNMembers( nMemb ); |
639 |
> |
the_molecules[lMolIndex].setStartAtom( molStart ); |
640 |
> |
the_molecules[lMolIndex].setEndAtom( molEnd ); |
641 |
> |
the_molecules[lMolIndex].setStampID( i ); |
642 |
> |
lMolIndex++; |
643 |
|
|
644 |
< |
// increment the index and repeat; |
558 |
< |
index++; |
644 |
> |
#ifdef IS_MPI |
645 |
|
} |
646 |
< |
|
647 |
< |
molEnd = index -1; |
562 |
< |
the_molecules[molIndex].setNMembers( nMemb ); |
563 |
< |
the_molecules[molIndex].setStartAtom( molStart ); |
564 |
< |
the_molecules[molIndex].setEndAtom( molEnd ); |
646 |
> |
#endif //is_mpi |
647 |
> |
|
648 |
|
molIndex++; |
566 |
– |
|
649 |
|
} |
650 |
|
} |
651 |
|
|
652 |
+ |
#ifdef IS_MPI |
653 |
+ |
for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] ); |
654 |
+ |
|
655 |
+ |
delete[] globalIndex; |
656 |
+ |
|
657 |
+ |
mpiSim->mpiRefresh(); |
658 |
+ |
#endif //IS_MPI |
659 |
+ |
|
660 |
|
the_ff->initializeAtoms(); |
661 |
|
} |
662 |
|
|
663 |
|
void SimSetup::makeBonds( void ){ |
664 |
|
|
665 |
< |
int i, j, k, index, offset; |
665 |
> |
int i, j, k, index, offset, molIndex; |
666 |
|
bond_pair* the_bonds; |
667 |
|
BondStamp* current_bond; |
668 |
|
|
669 |
|
the_bonds = new bond_pair[tot_bonds]; |
670 |
|
index = 0; |
671 |
|
offset = 0; |
672 |
+ |
molIndex = 0; |
673 |
+ |
|
674 |
|
for( i=0; i<n_components; i++ ){ |
675 |
|
|
676 |
|
for( j=0; j<components_nmol[i]; j++ ){ |
677 |
|
|
678 |
< |
for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){ |
679 |
< |
|
680 |
< |
current_bond = comp_stamps[i]->getBond( k ); |
681 |
< |
the_bonds[index].a = current_bond->getA() + offset; |
682 |
< |
the_bonds[index].b = current_bond->getB() + offset; |
683 |
< |
|
684 |
< |
the_excludes[index].i = the_bonds[index].a; |
685 |
< |
the_excludes[index].j = the_bonds[index].b; |
686 |
< |
|
687 |
< |
// increment the index and repeat; |
688 |
< |
index++; |
678 |
> |
#ifdef IS_MPI |
679 |
> |
if( mpiSim->getMyMolStart() <= molIndex && |
680 |
> |
molIndex <= mpiSim->getMyMolEnd() ){ |
681 |
> |
#endif // is_mpi |
682 |
> |
|
683 |
> |
for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){ |
684 |
> |
|
685 |
> |
current_bond = comp_stamps[i]->getBond( k ); |
686 |
> |
the_bonds[index].a = current_bond->getA() + offset; |
687 |
> |
the_bonds[index].b = current_bond->getB() + offset; |
688 |
> |
|
689 |
> |
the_excludes[index].i = the_bonds[index].a; |
690 |
> |
the_excludes[index].j = the_bonds[index].b; |
691 |
> |
|
692 |
> |
// increment the index and repeat; |
693 |
> |
index++; |
694 |
> |
} |
695 |
> |
offset += comp_stamps[i]->getNAtoms(); |
696 |
> |
|
697 |
> |
#ifdef IS_MPI |
698 |
|
} |
699 |
< |
offset += comp_stamps[i]->getNAtoms(); |
700 |
< |
} |
699 |
> |
#endif //is_mpi |
700 |
> |
|
701 |
> |
molIndex++; |
702 |
> |
} |
703 |
|
} |
704 |
|
|
705 |
|
the_ff->initializeBonds( the_bonds ); |
707 |
|
|
708 |
|
void SimSetup::makeBends( void ){ |
709 |
|
|
710 |
< |
int i, j, k, index, offset; |
710 |
> |
int i, j, k, index, offset, molIndex; |
711 |
|
bend_set* the_bends; |
712 |
|
BendStamp* current_bend; |
713 |
|
|
714 |
|
the_bends = new bend_set[tot_bends]; |
715 |
|
index = 0; |
716 |
|
offset = 0; |
717 |
+ |
molIndex = 0; |
718 |
|
for( i=0; i<n_components; i++ ){ |
719 |
|
|
720 |
|
for( j=0; j<components_nmol[i]; j++ ){ |
721 |
|
|
722 |
< |
for( k=0; k<comp_stamps[i]->getNBends(); k++ ){ |
722 |
> |
#ifdef IS_MPI |
723 |
> |
if( mpiSim->getMyMolStart() <= molIndex && |
724 |
> |
molIndex <= mpiSim->getMyMolEnd() ){ |
725 |
> |
#endif // is_mpi |
726 |
|
|
727 |
< |
current_bend = comp_stamps[i]->getBend( k ); |
728 |
< |
the_bends[index].a = current_bend->getA() + offset; |
729 |
< |
the_bends[index].b = current_bend->getB() + offset; |
730 |
< |
the_bends[index].c = current_bend->getC() + offset; |
731 |
< |
|
732 |
< |
the_excludes[index + tot_bonds].i = the_bends[index].a; |
733 |
< |
the_excludes[index + tot_bonds].j = the_bends[index].c; |
734 |
< |
|
735 |
< |
// increment the index and repeat; |
736 |
< |
index++; |
727 |
> |
for( k=0; k<comp_stamps[i]->getNBends(); k++ ){ |
728 |
> |
|
729 |
> |
current_bend = comp_stamps[i]->getBend( k ); |
730 |
> |
the_bends[index].a = current_bend->getA() + offset; |
731 |
> |
the_bends[index].b = current_bend->getB() + offset; |
732 |
> |
the_bends[index].c = current_bend->getC() + offset; |
733 |
> |
|
734 |
> |
the_excludes[index + tot_bonds].i = the_bends[index].a; |
735 |
> |
the_excludes[index + tot_bonds].j = the_bends[index].c; |
736 |
> |
|
737 |
> |
// increment the index and repeat; |
738 |
> |
index++; |
739 |
> |
} |
740 |
> |
offset += comp_stamps[i]->getNAtoms(); |
741 |
> |
|
742 |
> |
#ifdef IS_MPI |
743 |
|
} |
744 |
< |
offset += comp_stamps[i]->getNAtoms(); |
744 |
> |
#endif //is_mpi |
745 |
> |
|
746 |
> |
molIndex++; |
747 |
|
} |
748 |
|
} |
749 |
|
|
752 |
|
|
753 |
|
void SimSetup::makeTorsions( void ){ |
754 |
|
|
755 |
< |
int i, j, k, index, offset; |
755 |
> |
int i, j, k, index, offset, molIndex; |
756 |
|
torsion_set* the_torsions; |
757 |
|
TorsionStamp* current_torsion; |
758 |
|
|
759 |
|
the_torsions = new torsion_set[tot_torsions]; |
760 |
|
index = 0; |
761 |
|
offset = 0; |
762 |
+ |
molIndex = 0; |
763 |
|
for( i=0; i<n_components; i++ ){ |
764 |
|
|
765 |
|
for( j=0; j<components_nmol[i]; j++ ){ |
766 |
|
|
767 |
+ |
#ifdef IS_MPI |
768 |
+ |
if( mpiSim->getMyMolStart() <= molIndex && |
769 |
+ |
molIndex <= mpiSim->getMyMolEnd() ){ |
770 |
+ |
#endif // is_mpi |
771 |
+ |
|
772 |
|
for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){ |
773 |
|
|
774 |
|
current_torsion = comp_stamps[i]->getTorsion( k ); |
784 |
|
index++; |
785 |
|
} |
786 |
|
offset += comp_stamps[i]->getNAtoms(); |
787 |
+ |
|
788 |
+ |
#ifdef IS_MPI |
789 |
+ |
} |
790 |
+ |
#endif //is_mpi |
791 |
+ |
|
792 |
+ |
molIndex++; |
793 |
|
} |
794 |
|
} |
795 |
|
|