| 14 | 
  | 
 | 
| 15 | 
  | 
// some defines for ensemble and Forcefield  cases | 
| 16 | 
  | 
 | 
| 17 | 
< | 
#define NVE_ENS 0 | 
| 18 | 
< | 
#define NVT_ENS 1 | 
| 19 | 
< | 
#define NPT_ENS 2 | 
| 17 | 
> | 
#define NVE_ENS  0  | 
| 18 | 
> | 
#define NVT_ENS  1 | 
| 19 | 
> | 
#define NPTi_ENS 2 | 
| 20 | 
> | 
#define NPTf_ENS 3 | 
| 21 | 
  | 
 | 
| 22 | 
+ | 
 | 
| 23 | 
  | 
#define FF_DUFF 0 | 
| 24 | 
  | 
#define FF_LJ   1 | 
| 25 | 
  | 
 | 
| 99 | 
  | 
  n_components = the_globals->getNComponents(); | 
| 100 | 
  | 
  strcpy( force_field, the_globals->getForceField() ); | 
| 101 | 
  | 
 | 
| 102 | 
< | 
  if( !strcasecmp( force_field, "DUFF" ))    ffCase = FF_DUFF; | 
| 102 | 
> | 
  if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF; | 
| 103 | 
  | 
  else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ; | 
| 104 | 
  | 
  else{ | 
| 105 | 
  | 
    sprintf( painCave.errMsg, | 
| 112 | 
  | 
  // get the ensemble: | 
| 113 | 
  | 
  strcpy( ensemble, the_globals->getEnsemble() ); | 
| 114 | 
  | 
 | 
| 115 | 
< | 
  if( !strcasecmp( ensemble, "NVE" ))      ffCase = NVE_ENS; | 
| 116 | 
< | 
  else if( !strcasecmp( ensemble, "NVT" )) ffCase = NVT_ENS; | 
| 117 | 
< | 
  else if( !strcasecmp( ensemble, "NPT" )) ffCase = NPT_ENS; | 
| 115 | 
> | 
  if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS; | 
| 116 | 
> | 
  else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS; | 
| 117 | 
> | 
  else if( !strcasecmp( ensemble, "NPTi" ) || !strcasecmp( ensemble, "NPT") )  | 
| 118 | 
> | 
    ensembleCase = NPTi_ENS; | 
| 119 | 
> | 
  else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS; | 
| 120 | 
  | 
  else{ | 
| 121 | 
  | 
    sprintf( painCave.errMsg, | 
| 122 | 
  | 
             "SimSetup Warning. Unrecognized Ensemble -> %s, " | 
| 197 | 
  | 
    break; | 
| 198 | 
  | 
 | 
| 199 | 
  | 
  case FF_LJ: | 
| 200 | 
< | 
    the_ff = LJFF(); | 
| 200 | 
> | 
    the_ff = new LJFF(); | 
| 201 | 
  | 
    break; | 
| 202 | 
  | 
 | 
| 203 | 
  | 
  default: | 
| 483 | 
  | 
 | 
| 484 | 
  | 
  // get some of the tricky things that may still be in the globals | 
| 485 | 
  | 
 | 
| 486 | 
< | 
   | 
| 486 | 
> | 
  double boxVector[3]; | 
| 487 | 
  | 
  if( the_globals->haveBox() ){ | 
| 488 | 
< | 
    simnfo->box_x = the_globals->getBox(); | 
| 489 | 
< | 
    simnfo->box_y = the_globals->getBox(); | 
| 490 | 
< | 
    simnfo->box_z = the_globals->getBox(); | 
| 488 | 
> | 
    boxVector[0] = the_globals->getBox(); | 
| 489 | 
> | 
    boxVector[1] = the_globals->getBox(); | 
| 490 | 
> | 
    boxVector[2] = the_globals->getBox(); | 
| 491 | 
> | 
     | 
| 492 | 
> | 
    simnfo->setBox( boxVector ); | 
| 493 | 
  | 
  } | 
| 494 | 
  | 
  else if( the_globals->haveDensity() ){ | 
| 495 | 
  | 
 | 
| 496 | 
  | 
    double vol; | 
| 497 | 
  | 
    vol = (double)tot_nmol / the_globals->getDensity(); | 
| 498 | 
< | 
    simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) ); | 
| 499 | 
< | 
    simnfo->box_y = simnfo->box_x; | 
| 500 | 
< | 
    simnfo->box_z = simnfo->box_x; | 
| 498 | 
> | 
     boxVector[0] = pow( vol, ( 1.0 / 3.0 ) ); | 
| 499 | 
> | 
     boxVector[1] = boxVector[0]; | 
| 500 | 
> | 
     boxVector[2] = boxVector[0]; | 
| 501 | 
> | 
 | 
| 502 | 
> | 
    simnfo->setBox( boxVector ); | 
| 503 | 
  | 
  } | 
| 504 | 
  | 
  else{ | 
| 505 | 
  | 
    if( !the_globals->haveBoxX() ){ | 
| 508 | 
  | 
      painCave.isFatal = 1; | 
| 509 | 
  | 
      simError(); | 
| 510 | 
  | 
    } | 
| 511 | 
< | 
    simnfo->box_x = the_globals->getBoxX(); | 
| 511 | 
> | 
    boxVector[0] = the_globals->getBoxX(); | 
| 512 | 
  | 
 | 
| 513 | 
  | 
    if( !the_globals->haveBoxY() ){ | 
| 514 | 
  | 
      sprintf( painCave.errMsg, | 
| 516 | 
  | 
      painCave.isFatal = 1; | 
| 517 | 
  | 
      simError(); | 
| 518 | 
  | 
    } | 
| 519 | 
< | 
    simnfo->box_y = the_globals->getBoxY(); | 
| 519 | 
> | 
    boxVector[1] = the_globals->getBoxY(); | 
| 520 | 
  | 
 | 
| 521 | 
  | 
    if( !the_globals->haveBoxZ() ){ | 
| 522 | 
  | 
      sprintf( painCave.errMsg, | 
| 524 | 
  | 
      painCave.isFatal = 1; | 
| 525 | 
  | 
      simError(); | 
| 526 | 
  | 
    } | 
| 527 | 
< | 
    simnfo->box_z = the_globals->getBoxZ(); | 
| 527 | 
> | 
    boxVector[2] = the_globals->getBoxZ(); | 
| 528 | 
> | 
 | 
| 529 | 
> | 
    simnfo->setBox( boxVector ); | 
| 530 | 
  | 
  } | 
| 531 | 
  | 
 | 
| 532 | 
  | 
#ifdef IS_MPI | 
| 556 | 
  | 
      painCave.isFatal = 0; | 
| 557 | 
  | 
      simError(); | 
| 558 | 
  | 
      double smallest; | 
| 559 | 
< | 
      smallest = simnfo->box_x; | 
| 560 | 
< | 
      if (simnfo->box_y <= smallest) smallest = simnfo->box_y; | 
| 561 | 
< | 
      if (simnfo->box_z <= smallest) smallest = simnfo->box_z; | 
| 559 | 
> | 
      smallest = simnfo->boxLx; | 
| 560 | 
> | 
      if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy; | 
| 561 | 
> | 
      if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz; | 
| 562 | 
  | 
      simnfo->ecr = 0.5 * smallest; | 
| 563 | 
  | 
    } else { | 
| 564 | 
  | 
      simnfo->ecr        = the_globals->getECR(); | 
| 596 | 
  | 
        painCave.isFatal = 0; | 
| 597 | 
  | 
        simError(); | 
| 598 | 
  | 
        double smallest; | 
| 599 | 
< | 
        smallest = simnfo->box_x; | 
| 600 | 
< | 
        if (simnfo->box_y <= smallest) smallest = simnfo->box_y; | 
| 601 | 
< | 
        if (simnfo->box_z <= smallest) smallest = simnfo->box_z; | 
| 599 | 
> | 
        smallest = simnfo->boxLx; | 
| 600 | 
> | 
        if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy; | 
| 601 | 
> | 
        if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz; | 
| 602 | 
  | 
        simnfo->ecr = 0.5 * smallest; | 
| 603 | 
  | 
      } else { | 
| 604 | 
  | 
        simnfo->ecr        = the_globals->getECR(); | 
| 778 | 
  | 
 | 
| 779 | 
  | 
 | 
| 780 | 
  | 
  // make the integrator | 
| 781 | 
+ | 
   | 
| 782 | 
+ | 
   | 
| 783 | 
+ | 
  NVT*  myNVT = NULL; | 
| 784 | 
+ | 
  NPTi* myNPTi = NULL; | 
| 785 | 
+ | 
  NPTf* myNPTf = NULL; | 
| 786 | 
+ | 
  switch( ensembleCase ){ | 
| 787 | 
+ | 
 | 
| 788 | 
+ | 
  case NVE_ENS: | 
| 789 | 
+ | 
    new NVE( simnfo, the_ff ); | 
| 790 | 
+ | 
    break; | 
| 791 | 
  | 
 | 
| 792 | 
+ | 
  case NVT_ENS: | 
| 793 | 
+ | 
    myNVT = new NVT( simnfo, the_ff ); | 
| 794 | 
+ | 
    myNVT->setTargetTemp(the_globals->getTargetTemp()); | 
| 795 | 
  | 
 | 
| 796 | 
< | 
  if( !strcmp( ensemble, "TraPPE_Ex" ) ){ | 
| 797 | 
< | 
    new Symplectic(simnfo, the_ff, the_extendedsystem); | 
| 798 | 
< | 
  } | 
| 799 | 
< | 
  else if( !strcmp( force_field, "LJ" ) ){ | 
| 800 | 
< | 
    new Verlet( *simnfo, the_ff, the_extendedsystem ); | 
| 796 | 
> | 
    if (the_globals->haveTauThermostat())  | 
| 797 | 
> | 
      myNVT->setTauThermostat(the_globals->getTauThermostat()); | 
| 798 | 
> | 
 | 
| 799 | 
> | 
    else { | 
| 800 | 
> | 
      sprintf( painCave.errMsg, | 
| 801 | 
> | 
               "SimSetup error: If you use the NVT\n" | 
| 802 | 
> | 
               "    ensemble, you must set tauThermostat.\n"); | 
| 803 | 
> | 
      painCave.isFatal = 1; | 
| 804 | 
> | 
      simError(); | 
| 805 | 
> | 
    } | 
| 806 | 
> | 
    break; | 
| 807 | 
> | 
 | 
| 808 | 
> | 
  case NPTi_ENS: | 
| 809 | 
> | 
    myNPTi = new NPTi( simnfo, the_ff ); | 
| 810 | 
> | 
    myNPTi->setTargetTemp( the_globals->getTargetTemp()); | 
| 811 | 
> | 
 | 
| 812 | 
> | 
    if (the_globals->haveTargetPressure()) | 
| 813 | 
> | 
      myNPTi->setTargetPressure(the_globals->getTargetPressure()); | 
| 814 | 
> | 
    else { | 
| 815 | 
> | 
      sprintf( painCave.errMsg, | 
| 816 | 
> | 
               "SimSetup error: If you use a constant pressure\n" | 
| 817 | 
> | 
               "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 818 | 
> | 
      painCave.isFatal = 1; | 
| 819 | 
> | 
      simError(); | 
| 820 | 
> | 
    } | 
| 821 | 
> | 
     | 
| 822 | 
> | 
    if( the_globals->haveTauThermostat() ) | 
| 823 | 
> | 
      myNPTi->setTauThermostat( the_globals->getTauThermostat() ); | 
| 824 | 
> | 
    else{ | 
| 825 | 
> | 
      sprintf( painCave.errMsg, | 
| 826 | 
> | 
               "SimSetup error: If you use an NPT\n" | 
| 827 | 
> | 
               "    ensemble, you must set tauThermostat.\n"); | 
| 828 | 
> | 
      painCave.isFatal = 1; | 
| 829 | 
> | 
      simError(); | 
| 830 | 
> | 
    } | 
| 831 | 
> | 
 | 
| 832 | 
> | 
    if( the_globals->haveTauBarostat() ) | 
| 833 | 
> | 
      myNPTi->setTauBarostat( the_globals->getTauBarostat() ); | 
| 834 | 
> | 
    else{ | 
| 835 | 
> | 
      sprintf( painCave.errMsg, | 
| 836 | 
> | 
               "SimSetup error: If you use an NPT\n" | 
| 837 | 
> | 
               "    ensemble, you must set tauBarostat.\n"); | 
| 838 | 
> | 
      painCave.isFatal = 1; | 
| 839 | 
> | 
      simError(); | 
| 840 | 
> | 
    } | 
| 841 | 
> | 
    break; | 
| 842 | 
> | 
 | 
| 843 | 
> | 
  case NPTf_ENS: | 
| 844 | 
> | 
    myNPTf = new NPTf( simnfo, the_ff ); | 
| 845 | 
> | 
    myNPTf->setTargetTemp( the_globals->getTargetTemp()); | 
| 846 | 
> | 
 | 
| 847 | 
> | 
    if (the_globals->haveTargetPressure()) | 
| 848 | 
> | 
      myNPTf->setTargetPressure(the_globals->getTargetPressure()); | 
| 849 | 
> | 
    else { | 
| 850 | 
> | 
      sprintf( painCave.errMsg, | 
| 851 | 
> | 
               "SimSetup error: If you use a constant pressure\n" | 
| 852 | 
> | 
               "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 853 | 
> | 
      painCave.isFatal = 1; | 
| 854 | 
> | 
      simError(); | 
| 855 | 
> | 
    }     | 
| 856 | 
> | 
 | 
| 857 | 
> | 
    if( the_globals->haveTauThermostat() ) | 
| 858 | 
> | 
      myNPTf->setTauThermostat( the_globals->getTauThermostat() ); | 
| 859 | 
> | 
    else{ | 
| 860 | 
> | 
      sprintf( painCave.errMsg, | 
| 861 | 
> | 
               "SimSetup error: If you use an NPT\n" | 
| 862 | 
> | 
               "    ensemble, you must set tauThermostat.\n"); | 
| 863 | 
> | 
      painCave.isFatal = 1; | 
| 864 | 
> | 
      simError(); | 
| 865 | 
> | 
    } | 
| 866 | 
> | 
 | 
| 867 | 
> | 
    if( the_globals->haveTauBarostat() ) | 
| 868 | 
> | 
      myNPTf->setTauBarostat( the_globals->getTauBarostat() ); | 
| 869 | 
> | 
    else{ | 
| 870 | 
> | 
      sprintf( painCave.errMsg, | 
| 871 | 
> | 
               "SimSetup error: If you use an NPT\n" | 
| 872 | 
> | 
               "    ensemble, you must set tauBarostat.\n"); | 
| 873 | 
> | 
      painCave.isFatal = 1; | 
| 874 | 
> | 
      simError(); | 
| 875 | 
> | 
    } | 
| 876 | 
> | 
    break; | 
| 877 | 
> | 
     | 
| 878 | 
> | 
     | 
| 879 | 
> | 
 | 
| 880 | 
> | 
  default: | 
| 881 | 
> | 
    sprintf( painCave.errMsg, | 
| 882 | 
> | 
             "SimSetup Error. Unrecognized ensemble in case statement.\n"); | 
| 883 | 
> | 
    painCave.isFatal = 1; | 
| 884 | 
> | 
    simError(); | 
| 885 | 
  | 
  } | 
| 886 | 
+ | 
 | 
| 887 | 
  | 
 | 
| 888 | 
  | 
#ifdef IS_MPI | 
| 889 | 
  | 
  mpiSim->mpiRefresh(); | 
| 1197 | 
  | 
    have_extra =1; | 
| 1198 | 
  | 
 | 
| 1199 | 
  | 
    n_cells = (int)temp3 - 1; | 
| 1200 | 
< | 
    cellx = simnfo->box_x / temp3; | 
| 1201 | 
< | 
    celly = simnfo->box_y / temp3; | 
| 1202 | 
< | 
    cellz = simnfo->box_z / temp3; | 
| 1200 | 
> | 
    cellx = simnfo->boxLx / temp3; | 
| 1201 | 
> | 
    celly = simnfo->boxLy / temp3; | 
| 1202 | 
> | 
    cellz = simnfo->boxLz / temp3; | 
| 1203 | 
  | 
    n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells ); | 
| 1204 | 
  | 
    temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) ); | 
| 1205 | 
  | 
    n_per_extra = (int)ceil( temp1 ); | 
| 1214 | 
  | 
  } | 
| 1215 | 
  | 
  else{ | 
| 1216 | 
  | 
    n_cells = (int)temp3; | 
| 1217 | 
< | 
    cellx = simnfo->box_x / temp3; | 
| 1218 | 
< | 
    celly = simnfo->box_y / temp3; | 
| 1219 | 
< | 
    cellz = simnfo->box_z / temp3; | 
| 1217 | 
> | 
    cellx = simnfo->boxLx / temp3; | 
| 1218 | 
> | 
    celly = simnfo->boxLy / temp3; | 
| 1219 | 
> | 
    cellz = simnfo->boxLz / temp3; | 
| 1220 | 
  | 
  } | 
| 1221 | 
  | 
 | 
| 1222 | 
  | 
  current_mol = 0; |