333 |
|
|
334 |
|
// get some of the tricky things that may still be in the globals |
335 |
|
|
336 |
– |
if( simnfo->n_dipoles ){ |
337 |
– |
|
338 |
– |
if( !the_globals->haveRRF() ){ |
339 |
– |
sprintf( painCave.errMsg, |
340 |
– |
"SimSetup Error, system has dipoles, but no rRF was set.\n"); |
341 |
– |
painCave.isFatal = 1; |
342 |
– |
simError(); |
343 |
– |
} |
344 |
– |
if( !the_globals->haveDielectric() ){ |
345 |
– |
sprintf( painCave.errMsg, |
346 |
– |
"SimSetup Error, system has dipoles, but no" |
347 |
– |
" dielectric was set.\n" ); |
348 |
– |
painCave.isFatal = 1; |
349 |
– |
simError(); |
350 |
– |
} |
351 |
– |
|
352 |
– |
simnfo->rRF = the_globals->getRRF(); |
353 |
– |
simnfo->dielectric = the_globals->getDielectric(); |
354 |
– |
} |
355 |
– |
|
356 |
– |
#ifdef IS_MPI |
357 |
– |
strcpy( checkPointMsg, "rRf and dielectric check out" ); |
358 |
– |
MPIcheckPoint(); |
359 |
– |
#endif // is_mpi |
336 |
|
|
337 |
|
if( the_globals->haveBox() ){ |
338 |
|
simnfo->box_x = the_globals->getBox(); |
402 |
|
} |
403 |
|
|
404 |
|
|
405 |
< |
|
405 |
> |
if (the_globals->getUseRF() ) { |
406 |
> |
simnfo->useReactionField = 1; |
407 |
> |
|
408 |
> |
if( !the_globals->haveECR() ){ |
409 |
> |
sprintf( painCave.errMsg, |
410 |
> |
"SimSetup Warning: using default value of 1/2 the smallest " |
411 |
> |
"box length for the electrostaticCutoffRadius.\n" |
412 |
> |
"I hope you have a very fast processor!\n"); |
413 |
> |
painCave.isFatal = 0; |
414 |
> |
simError(); |
415 |
> |
double smallest; |
416 |
> |
smallest = simnfo->box_x; |
417 |
> |
if (simnfo->box_y <= smallest) smallest = simnfo->box_y; |
418 |
> |
if (simnfo->box_z <= smallest) smallest = simnfo->box_z; |
419 |
> |
simnfo->ecr = 0.5 * smallest; |
420 |
> |
} else { |
421 |
> |
simnfo->ecr = the_globals->getECR(); |
422 |
> |
} |
423 |
|
|
424 |
+ |
if( !the_globals->haveEST() ){ |
425 |
+ |
sprintf( painCave.errMsg, |
426 |
+ |
"SimSetup Warning: using default value of 0.05 * the " |
427 |
+ |
"electrostaticCutoffRadius for the electrostaticSkinThickness\n" |
428 |
+ |
); |
429 |
+ |
painCave.isFatal = 0; |
430 |
+ |
simError(); |
431 |
+ |
simnfo->est = 0.05 * simnfo->ecr; |
432 |
+ |
} else { |
433 |
+ |
simnfo->est = the_globals->getEST(); |
434 |
+ |
} |
435 |
+ |
|
436 |
+ |
if(!the_globals->haveDielectric() ){ |
437 |
+ |
sprintf( painCave.errMsg, |
438 |
+ |
"SimSetup Error: You are trying to use Reaction Field without" |
439 |
+ |
"setting a dielectric constant!\n" |
440 |
+ |
); |
441 |
+ |
painCave.isFatal = 1; |
442 |
+ |
simError(); |
443 |
+ |
} |
444 |
+ |
simnfo->dielectric = the_globals->getDielectric(); |
445 |
+ |
} else { |
446 |
+ |
if (simnfo->n_dipoles) { |
447 |
+ |
|
448 |
+ |
if( !the_globals->haveECR() ){ |
449 |
+ |
sprintf( painCave.errMsg, |
450 |
+ |
"SimSetup Warning: using default value of 1/2 the smallest" |
451 |
+ |
"box length for the electrostaticCutoffRadius.\n" |
452 |
+ |
"I hope you have a very fast processor!\n"); |
453 |
+ |
painCave.isFatal = 0; |
454 |
+ |
simError(); |
455 |
+ |
double smallest; |
456 |
+ |
smallest = simnfo->box_x; |
457 |
+ |
if (simnfo->box_y <= smallest) smallest = simnfo->box_y; |
458 |
+ |
if (simnfo->box_z <= smallest) smallest = simnfo->box_z; |
459 |
+ |
simnfo->ecr = 0.5 * smallest; |
460 |
+ |
} else { |
461 |
+ |
simnfo->ecr = the_globals->getECR(); |
462 |
+ |
} |
463 |
+ |
|
464 |
+ |
if( !the_globals->haveEST() ){ |
465 |
+ |
sprintf( painCave.errMsg, |
466 |
+ |
"SimSetup Warning: using default value of 5% of the" |
467 |
+ |
"electrostaticCutoffRadius for the " |
468 |
+ |
"electrostaticSkinThickness\n" |
469 |
+ |
); |
470 |
+ |
painCave.isFatal = 0; |
471 |
+ |
simError(); |
472 |
+ |
simnfo->est = 0.05 * simnfo->ecr; |
473 |
+ |
} else { |
474 |
+ |
simnfo->est = the_globals->getEST(); |
475 |
+ |
} |
476 |
+ |
} |
477 |
+ |
} |
478 |
|
|
479 |
+ |
#ifdef IS_MPI |
480 |
+ |
strcpy( checkPointMsg, "electrostatic parameters check out" ); |
481 |
+ |
MPIcheckPoint(); |
482 |
+ |
#endif // is_mpi |
483 |
|
|
484 |
|
if( the_globals->haveInitialConfig() ){ |
485 |
|
|
671 |
|
#endif // is_mpi |
672 |
|
} |
673 |
|
|
674 |
+ |
|
675 |
+ |
void SimSetup::makeMolecules( void ){ |
676 |
+ |
|
677 |
+ |
int i, j, exI, exJ, tempEx, stampID, atomOffset; |
678 |
+ |
molInit info; |
679 |
+ |
DirectionalAtom* dAtom; |
680 |
+ |
AtomStamp* currentAtom; |
681 |
+ |
BondStamp* currentBond; |
682 |
+ |
BendStamp* currentBend; |
683 |
+ |
TorsionStamp* currentTorsion; |
684 |
+ |
|
685 |
+ |
//init the forceField paramters |
686 |
+ |
|
687 |
+ |
the_ff->readParams(); |
688 |
+ |
|
689 |
+ |
|
690 |
+ |
// init the molecules |
691 |
+ |
|
692 |
+ |
atomOffset = 0; |
693 |
+ |
for(i=0; i<simnfo->n_mol; i++){ |
694 |
+ |
|
695 |
+ |
stampID = the_molecules[i].getStampID(); |
696 |
+ |
|
697 |
+ |
info.nAtoms = comp_stamps[stampID]->getNAtoms(); |
698 |
+ |
info.nBonds = comp_stamps[stampID]->getNBonds(); |
699 |
+ |
info.nBends = comp_stamps[stampID]->getNBends(); |
700 |
+ |
info.nTorsions = comp_stamps[stampID]->getNTorsions(); |
701 |
+ |
|
702 |
+ |
info.myAtoms = &the_atoms[atomOffset]; |
703 |
+ |
info.myBonds = new Bond*[info.nBonds]; |
704 |
+ |
info.myBends = new Bend*[info.nBends]; |
705 |
+ |
info.myTorsions = new Torsions*[info.nTorsions]; |
706 |
+ |
|
707 |
+ |
theBonds = new bond_pair[info.nBonds]; |
708 |
+ |
theBends = new bend_set[info.nBends]; |
709 |
+ |
theTorsions = new torsion_set[info.nTorsions]; |
710 |
+ |
|
711 |
+ |
// make the Atoms |
712 |
+ |
|
713 |
+ |
for(j=0; j<info.nAtoms; j++){ |
714 |
+ |
|
715 |
+ |
currentAtom = theComponents[stampID]->getAtom( j ); |
716 |
+ |
if( currentAtom->haveOrientation() ){ |
717 |
+ |
|
718 |
+ |
dAtom = new DirectionalAtom(j + atomOffset); |
719 |
+ |
simnfo->n_oriented++; |
720 |
+ |
info.myAtoms[j] = dAtom; |
721 |
+ |
|
722 |
+ |
ux = currentAtom->getOrntX(); |
723 |
+ |
uy = currentAtom->getOrntY(); |
724 |
+ |
uz = currentAtom->getOrntZ(); |
725 |
+ |
|
726 |
+ |
uSqr = (ux * ux) + (uy * uy) + (uz * uz); |
727 |
+ |
|
728 |
+ |
u = sqrt( uSqr ); |
729 |
+ |
ux = ux / u; |
730 |
+ |
uy = uy / u; |
731 |
+ |
uz = uz / u; |
732 |
+ |
|
733 |
+ |
dAtom->setSUx( ux ); |
734 |
+ |
dAtom->setSUy( uy ); |
735 |
+ |
dAtom->setSUz( uz ); |
736 |
+ |
} |
737 |
+ |
else{ |
738 |
+ |
info.myAtoms[j] = new GeneralAtom(j + atomOffset); |
739 |
+ |
} |
740 |
+ |
info.myAtoms[j]->setType( currentAtom->getType() ); |
741 |
+ |
|
742 |
+ |
#ifdef IS_MPI |
743 |
+ |
|
744 |
+ |
info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] ); |
745 |
+ |
|
746 |
+ |
#endif // is_mpi |
747 |
+ |
} |
748 |
+ |
|
749 |
+ |
// make the bonds |
750 |
+ |
for(j=0; j<nBonds; j++){ |
751 |
+ |
|
752 |
+ |
currentBond = comp_stamps[stampID]->getBond( j ); |
753 |
+ |
theBonds[j].a = currentBond->getA() + atomOffset; |
754 |
+ |
theBonds[j].b = currentBond->getB() + atomOffset; |
755 |
+ |
|
756 |
+ |
exI = theBonds[i].a; |
757 |
+ |
exJ = theBonds[i].b; |
758 |
+ |
|
759 |
+ |
// exclude_I must always be the smaller of the pair |
760 |
+ |
if( exI > exJ ){ |
761 |
+ |
tempEx = exI; |
762 |
+ |
exI = exJ; |
763 |
+ |
exJ = tempEx; |
764 |
+ |
} |
765 |
+ |
#ifdef IS_MPI |
766 |
+ |
|
767 |
+ |
the_excludes[index*2] = |
768 |
+ |
the_atoms[exI]->getGlobalIndex() + 1; |
769 |
+ |
the_excludes[index*2 + 1] = |
770 |
+ |
the_atoms[exJ]->getGlobalIndex() + 1; |
771 |
+ |
|
772 |
+ |
#else // isn't MPI |
773 |
+ |
|
774 |
+ |
the_excludes[index*2] = exI + 1; |
775 |
+ |
the_excludes[index*2 + 1] = exJ + 1; |
776 |
+ |
// fortran index from 1 (hence the +1 in the indexing) |
777 |
+ |
|
778 |
+ |
#endif //is_mpi |
779 |
+ |
|
780 |
+ |
} |
781 |
+ |
|
782 |
+ |
|
783 |
+ |
|
784 |
+ |
|
785 |
+ |
|
786 |
+ |
|
787 |
+ |
|
788 |
+ |
|
789 |
+ |
|
790 |
+ |
|
791 |
+ |
|
792 |
+ |
|
793 |
+ |
|
794 |
+ |
|
795 |
+ |
|
796 |
|
void SimSetup::makeAtoms( void ){ |
797 |
|
|
798 |
|
int i, j, k, index; |