| 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 | 
> | 
#define NPTim_ENS 4 | 
| 22 | 
> | 
#define NPTfm_ENS 5 | 
| 23 | 
  | 
 | 
| 24 | 
+ | 
 | 
| 25 | 
  | 
#define FF_DUFF 0 | 
| 26 | 
  | 
#define FF_LJ   1 | 
| 27 | 
  | 
 | 
| 101 | 
  | 
  n_components = the_globals->getNComponents(); | 
| 102 | 
  | 
  strcpy( force_field, the_globals->getForceField() ); | 
| 103 | 
  | 
 | 
| 104 | 
< | 
  if( !strcasecmp( force_field, "DUFF" ))    ffCase = FF_DUFF; | 
| 104 | 
> | 
  if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF; | 
| 105 | 
  | 
  else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ; | 
| 106 | 
  | 
  else{ | 
| 107 | 
  | 
    sprintf( painCave.errMsg, | 
| 114 | 
  | 
  // get the ensemble: | 
| 115 | 
  | 
  strcpy( ensemble, the_globals->getEnsemble() ); | 
| 116 | 
  | 
 | 
| 117 | 
< | 
  if( !strcasecmp( ensemble, "NVE" ))      ffCase = NVE_ENS; | 
| 118 | 
< | 
  else if( !strcasecmp( ensemble, "NVT" )) ffCase = NVT_ENS; | 
| 119 | 
< | 
  else if( !strcasecmp( ensemble, "NPT" )) ffCase = NPT_ENS; | 
| 117 | 
> | 
  if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS; | 
| 118 | 
> | 
  else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS; | 
| 119 | 
> | 
  else if( !strcasecmp( ensemble, "NPTi" ) || !strcasecmp( ensemble, "NPT") )  | 
| 120 | 
> | 
    ensembleCase = NPTi_ENS; | 
| 121 | 
> | 
  else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS; | 
| 122 | 
> | 
  else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS; | 
| 123 | 
> | 
  else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS; | 
| 124 | 
  | 
  else{ | 
| 125 | 
  | 
    sprintf( painCave.errMsg, | 
| 126 | 
  | 
             "SimSetup Warning. Unrecognized Ensemble -> %s, " | 
| 201 | 
  | 
    break; | 
| 202 | 
  | 
 | 
| 203 | 
  | 
  case FF_LJ: | 
| 204 | 
< | 
    the_ff = LJFF(); | 
| 204 | 
> | 
    the_ff = new LJFF(); | 
| 205 | 
  | 
    break; | 
| 206 | 
  | 
 | 
| 207 | 
  | 
  default: | 
| 487 | 
  | 
 | 
| 488 | 
  | 
  // get some of the tricky things that may still be in the globals | 
| 489 | 
  | 
 | 
| 490 | 
< | 
   | 
| 490 | 
> | 
  double boxVector[3]; | 
| 491 | 
  | 
  if( the_globals->haveBox() ){ | 
| 492 | 
< | 
    simnfo->box_x = the_globals->getBox(); | 
| 493 | 
< | 
    simnfo->box_y = the_globals->getBox(); | 
| 494 | 
< | 
    simnfo->box_z = the_globals->getBox(); | 
| 492 | 
> | 
    boxVector[0] = the_globals->getBox(); | 
| 493 | 
> | 
    boxVector[1] = the_globals->getBox(); | 
| 494 | 
> | 
    boxVector[2] = the_globals->getBox(); | 
| 495 | 
> | 
     | 
| 496 | 
> | 
    simnfo->setBox( boxVector ); | 
| 497 | 
  | 
  } | 
| 498 | 
  | 
  else if( the_globals->haveDensity() ){ | 
| 499 | 
  | 
 | 
| 500 | 
  | 
    double vol; | 
| 501 | 
  | 
    vol = (double)tot_nmol / the_globals->getDensity(); | 
| 502 | 
< | 
    simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) ); | 
| 503 | 
< | 
    simnfo->box_y = simnfo->box_x; | 
| 504 | 
< | 
    simnfo->box_z = simnfo->box_x; | 
| 502 | 
> | 
     boxVector[0] = pow( vol, ( 1.0 / 3.0 ) ); | 
| 503 | 
> | 
     boxVector[1] = boxVector[0]; | 
| 504 | 
> | 
     boxVector[2] = boxVector[0]; | 
| 505 | 
> | 
 | 
| 506 | 
> | 
    simnfo->setBox( boxVector ); | 
| 507 | 
  | 
  } | 
| 508 | 
  | 
  else{ | 
| 509 | 
  | 
    if( !the_globals->haveBoxX() ){ | 
| 512 | 
  | 
      painCave.isFatal = 1; | 
| 513 | 
  | 
      simError(); | 
| 514 | 
  | 
    } | 
| 515 | 
< | 
    simnfo->box_x = the_globals->getBoxX(); | 
| 515 | 
> | 
    boxVector[0] = the_globals->getBoxX(); | 
| 516 | 
  | 
 | 
| 517 | 
  | 
    if( !the_globals->haveBoxY() ){ | 
| 518 | 
  | 
      sprintf( painCave.errMsg, | 
| 520 | 
  | 
      painCave.isFatal = 1; | 
| 521 | 
  | 
      simError(); | 
| 522 | 
  | 
    } | 
| 523 | 
< | 
    simnfo->box_y = the_globals->getBoxY(); | 
| 523 | 
> | 
    boxVector[1] = the_globals->getBoxY(); | 
| 524 | 
  | 
 | 
| 525 | 
  | 
    if( !the_globals->haveBoxZ() ){ | 
| 526 | 
  | 
      sprintf( painCave.errMsg, | 
| 528 | 
  | 
      painCave.isFatal = 1; | 
| 529 | 
  | 
      simError(); | 
| 530 | 
  | 
    } | 
| 531 | 
< | 
    simnfo->box_z = the_globals->getBoxZ(); | 
| 531 | 
> | 
    boxVector[2] = the_globals->getBoxZ(); | 
| 532 | 
> | 
 | 
| 533 | 
> | 
    simnfo->setBox( boxVector ); | 
| 534 | 
  | 
  } | 
| 535 | 
  | 
 | 
| 536 | 
  | 
#ifdef IS_MPI | 
| 560 | 
  | 
      painCave.isFatal = 0; | 
| 561 | 
  | 
      simError(); | 
| 562 | 
  | 
      double smallest; | 
| 563 | 
< | 
      smallest = simnfo->box_x; | 
| 564 | 
< | 
      if (simnfo->box_y <= smallest) smallest = simnfo->box_y; | 
| 565 | 
< | 
      if (simnfo->box_z <= smallest) smallest = simnfo->box_z; | 
| 563 | 
> | 
      smallest = simnfo->boxLx; | 
| 564 | 
> | 
      if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy; | 
| 565 | 
> | 
      if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz; | 
| 566 | 
  | 
      simnfo->ecr = 0.5 * smallest; | 
| 567 | 
  | 
    } else { | 
| 568 | 
  | 
      simnfo->ecr        = the_globals->getECR(); | 
| 600 | 
  | 
        painCave.isFatal = 0; | 
| 601 | 
  | 
        simError(); | 
| 602 | 
  | 
        double smallest; | 
| 603 | 
< | 
        smallest = simnfo->box_x; | 
| 604 | 
< | 
        if (simnfo->box_y <= smallest) smallest = simnfo->box_y; | 
| 605 | 
< | 
        if (simnfo->box_z <= smallest) smallest = simnfo->box_z; | 
| 603 | 
> | 
        smallest = simnfo->boxLx; | 
| 604 | 
> | 
        if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy; | 
| 605 | 
> | 
        if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz; | 
| 606 | 
  | 
        simnfo->ecr = 0.5 * smallest; | 
| 607 | 
  | 
      } else { | 
| 608 | 
  | 
        simnfo->ecr        = the_globals->getECR(); | 
| 782 | 
  | 
 | 
| 783 | 
  | 
 | 
| 784 | 
  | 
  // make the integrator | 
| 785 | 
+ | 
   | 
| 786 | 
+ | 
   | 
| 787 | 
+ | 
  NVT*  myNVT = NULL; | 
| 788 | 
+ | 
  NPTi* myNPTi = NULL; | 
| 789 | 
+ | 
  NPTf* myNPTf = NULL; | 
| 790 | 
+ | 
  NPTim* myNPTim = NULL; | 
| 791 | 
+ | 
  NPTfm* myNPTfm = NULL; | 
| 792 | 
  | 
 | 
| 793 | 
  | 
  switch( ensembleCase ){ | 
| 794 | 
  | 
 | 
| 796 | 
  | 
    new NVE( simnfo, the_ff ); | 
| 797 | 
  | 
    break; | 
| 798 | 
  | 
 | 
| 799 | 
+ | 
  case NVT_ENS: | 
| 800 | 
+ | 
    myNVT = new NVT( simnfo, the_ff ); | 
| 801 | 
+ | 
    myNVT->setTargetTemp(the_globals->getTargetTemp()); | 
| 802 | 
+ | 
 | 
| 803 | 
+ | 
    if (the_globals->haveTauThermostat())  | 
| 804 | 
+ | 
      myNVT->setTauThermostat(the_globals->getTauThermostat()); | 
| 805 | 
+ | 
 | 
| 806 | 
+ | 
    else { | 
| 807 | 
+ | 
      sprintf( painCave.errMsg, | 
| 808 | 
+ | 
               "SimSetup error: If you use the NVT\n" | 
| 809 | 
+ | 
               "    ensemble, you must set tauThermostat.\n"); | 
| 810 | 
+ | 
      painCave.isFatal = 1; | 
| 811 | 
+ | 
      simError(); | 
| 812 | 
+ | 
    } | 
| 813 | 
+ | 
    break; | 
| 814 | 
+ | 
 | 
| 815 | 
+ | 
  case NPTi_ENS: | 
| 816 | 
+ | 
    myNPTi = new NPTi( simnfo, the_ff ); | 
| 817 | 
+ | 
    myNPTi->setTargetTemp( the_globals->getTargetTemp()); | 
| 818 | 
+ | 
 | 
| 819 | 
+ | 
    if (the_globals->haveTargetPressure()) | 
| 820 | 
+ | 
      myNPTi->setTargetPressure(the_globals->getTargetPressure()); | 
| 821 | 
+ | 
    else { | 
| 822 | 
+ | 
      sprintf( painCave.errMsg, | 
| 823 | 
+ | 
               "SimSetup error: If you use a constant pressure\n" | 
| 824 | 
+ | 
               "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 825 | 
+ | 
      painCave.isFatal = 1; | 
| 826 | 
+ | 
      simError(); | 
| 827 | 
+ | 
    } | 
| 828 | 
+ | 
     | 
| 829 | 
+ | 
    if( the_globals->haveTauThermostat() ) | 
| 830 | 
+ | 
      myNPTi->setTauThermostat( the_globals->getTauThermostat() ); | 
| 831 | 
+ | 
    else{ | 
| 832 | 
+ | 
      sprintf( painCave.errMsg, | 
| 833 | 
+ | 
               "SimSetup error: If you use an NPT\n" | 
| 834 | 
+ | 
               "    ensemble, you must set tauThermostat.\n"); | 
| 835 | 
+ | 
      painCave.isFatal = 1; | 
| 836 | 
+ | 
      simError(); | 
| 837 | 
+ | 
    } | 
| 838 | 
+ | 
 | 
| 839 | 
+ | 
    if( the_globals->haveTauBarostat() ) | 
| 840 | 
+ | 
      myNPTi->setTauBarostat( the_globals->getTauBarostat() ); | 
| 841 | 
+ | 
    else{ | 
| 842 | 
+ | 
      sprintf( painCave.errMsg, | 
| 843 | 
+ | 
               "SimSetup error: If you use an NPT\n" | 
| 844 | 
+ | 
               "    ensemble, you must set tauBarostat.\n"); | 
| 845 | 
+ | 
      painCave.isFatal = 1; | 
| 846 | 
+ | 
      simError(); | 
| 847 | 
+ | 
    } | 
| 848 | 
+ | 
    break; | 
| 849 | 
+ | 
 | 
| 850 | 
+ | 
  case NPTf_ENS: | 
| 851 | 
+ | 
    myNPTf = new NPTf( simnfo, the_ff ); | 
| 852 | 
+ | 
    myNPTf->setTargetTemp( the_globals->getTargetTemp()); | 
| 853 | 
+ | 
 | 
| 854 | 
+ | 
    if (the_globals->haveTargetPressure()) | 
| 855 | 
+ | 
      myNPTf->setTargetPressure(the_globals->getTargetPressure()); | 
| 856 | 
+ | 
    else { | 
| 857 | 
+ | 
      sprintf( painCave.errMsg, | 
| 858 | 
+ | 
               "SimSetup error: If you use a constant pressure\n" | 
| 859 | 
+ | 
               "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 860 | 
+ | 
      painCave.isFatal = 1; | 
| 861 | 
+ | 
      simError(); | 
| 862 | 
+ | 
    }     | 
| 863 | 
+ | 
 | 
| 864 | 
+ | 
    if( the_globals->haveTauThermostat() ) | 
| 865 | 
+ | 
      myNPTf->setTauThermostat( the_globals->getTauThermostat() ); | 
| 866 | 
+ | 
    else{ | 
| 867 | 
+ | 
      sprintf( painCave.errMsg, | 
| 868 | 
+ | 
               "SimSetup error: If you use an NPT\n" | 
| 869 | 
+ | 
               "    ensemble, you must set tauThermostat.\n"); | 
| 870 | 
+ | 
      painCave.isFatal = 1; | 
| 871 | 
+ | 
      simError(); | 
| 872 | 
+ | 
    } | 
| 873 | 
+ | 
 | 
| 874 | 
+ | 
    if( the_globals->haveTauBarostat() ) | 
| 875 | 
+ | 
      myNPTf->setTauBarostat( the_globals->getTauBarostat() ); | 
| 876 | 
+ | 
    else{ | 
| 877 | 
+ | 
      sprintf( painCave.errMsg, | 
| 878 | 
+ | 
               "SimSetup error: If you use an NPT\n" | 
| 879 | 
+ | 
               "    ensemble, you must set tauBarostat.\n"); | 
| 880 | 
+ | 
      painCave.isFatal = 1; | 
| 881 | 
+ | 
      simError(); | 
| 882 | 
+ | 
    } | 
| 883 | 
+ | 
    break; | 
| 884 | 
+ | 
     | 
| 885 | 
+ | 
  case NPTim_ENS: | 
| 886 | 
+ | 
    myNPTim = new NPTim( simnfo, the_ff ); | 
| 887 | 
+ | 
    myNPTim->setTargetTemp( the_globals->getTargetTemp()); | 
| 888 | 
+ | 
 | 
| 889 | 
+ | 
    if (the_globals->haveTargetPressure()) | 
| 890 | 
+ | 
      myNPTim->setTargetPressure(the_globals->getTargetPressure()); | 
| 891 | 
+ | 
    else { | 
| 892 | 
+ | 
      sprintf( painCave.errMsg, | 
| 893 | 
+ | 
               "SimSetup error: If you use a constant pressure\n" | 
| 894 | 
+ | 
               "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 895 | 
+ | 
      painCave.isFatal = 1; | 
| 896 | 
+ | 
      simError(); | 
| 897 | 
+ | 
    } | 
| 898 | 
+ | 
     | 
| 899 | 
+ | 
    if( the_globals->haveTauThermostat() ) | 
| 900 | 
+ | 
      myNPTim->setTauThermostat( the_globals->getTauThermostat() ); | 
| 901 | 
+ | 
    else{ | 
| 902 | 
+ | 
      sprintf( painCave.errMsg, | 
| 903 | 
+ | 
               "SimSetup error: If you use an NPT\n" | 
| 904 | 
+ | 
               "    ensemble, you must set tauThermostat.\n"); | 
| 905 | 
+ | 
      painCave.isFatal = 1; | 
| 906 | 
+ | 
      simError(); | 
| 907 | 
+ | 
    } | 
| 908 | 
+ | 
 | 
| 909 | 
+ | 
    if( the_globals->haveTauBarostat() ) | 
| 910 | 
+ | 
      myNPTim->setTauBarostat( the_globals->getTauBarostat() ); | 
| 911 | 
+ | 
    else{ | 
| 912 | 
+ | 
      sprintf( painCave.errMsg, | 
| 913 | 
+ | 
               "SimSetup error: If you use an NPT\n" | 
| 914 | 
+ | 
               "    ensemble, you must set tauBarostat.\n"); | 
| 915 | 
+ | 
      painCave.isFatal = 1; | 
| 916 | 
+ | 
      simError(); | 
| 917 | 
+ | 
    } | 
| 918 | 
+ | 
    break; | 
| 919 | 
+ | 
 | 
| 920 | 
+ | 
  case NPTfm_ENS: | 
| 921 | 
+ | 
    myNPTfm = new NPTfm( simnfo, the_ff ); | 
| 922 | 
+ | 
    myNPTfm->setTargetTemp( the_globals->getTargetTemp()); | 
| 923 | 
+ | 
 | 
| 924 | 
+ | 
    if (the_globals->haveTargetPressure()) | 
| 925 | 
+ | 
      myNPTfm->setTargetPressure(the_globals->getTargetPressure()); | 
| 926 | 
+ | 
    else { | 
| 927 | 
+ | 
      sprintf( painCave.errMsg, | 
| 928 | 
+ | 
               "SimSetup error: If you use a constant pressure\n" | 
| 929 | 
+ | 
               "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 930 | 
+ | 
      painCave.isFatal = 1; | 
| 931 | 
+ | 
      simError(); | 
| 932 | 
+ | 
    } | 
| 933 | 
+ | 
     | 
| 934 | 
+ | 
    if( the_globals->haveTauThermostat() ) | 
| 935 | 
+ | 
      myNPTfm->setTauThermostat( the_globals->getTauThermostat() ); | 
| 936 | 
+ | 
    else{ | 
| 937 | 
+ | 
      sprintf( painCave.errMsg, | 
| 938 | 
+ | 
               "SimSetup error: If you use an NPT\n" | 
| 939 | 
+ | 
               "    ensemble, you must set tauThermostat.\n"); | 
| 940 | 
+ | 
      painCave.isFatal = 1; | 
| 941 | 
+ | 
      simError(); | 
| 942 | 
+ | 
    } | 
| 943 | 
+ | 
 | 
| 944 | 
+ | 
    if( the_globals->haveTauBarostat() ) | 
| 945 | 
+ | 
      myNPTfm->setTauBarostat( the_globals->getTauBarostat() ); | 
| 946 | 
+ | 
    else{ | 
| 947 | 
+ | 
      sprintf( painCave.errMsg, | 
| 948 | 
+ | 
               "SimSetup error: If you use an NPT\n" | 
| 949 | 
+ | 
               "    ensemble, you must set tauBarostat.\n"); | 
| 950 | 
+ | 
      painCave.isFatal = 1; | 
| 951 | 
+ | 
      simError(); | 
| 952 | 
+ | 
    } | 
| 953 | 
+ | 
    break; | 
| 954 | 
+ | 
 | 
| 955 | 
  | 
  default: | 
| 956 | 
  | 
    sprintf( painCave.errMsg, | 
| 957 | 
  | 
             "SimSetup Error. Unrecognized ensemble in case statement.\n"); | 
| 1272 | 
  | 
    have_extra =1; | 
| 1273 | 
  | 
 | 
| 1274 | 
  | 
    n_cells = (int)temp3 - 1; | 
| 1275 | 
< | 
    cellx = simnfo->box_x / temp3; | 
| 1276 | 
< | 
    celly = simnfo->box_y / temp3; | 
| 1277 | 
< | 
    cellz = simnfo->box_z / temp3; | 
| 1275 | 
> | 
    cellx = simnfo->boxLx / temp3; | 
| 1276 | 
> | 
    celly = simnfo->boxLy / temp3; | 
| 1277 | 
> | 
    cellz = simnfo->boxLz / temp3; | 
| 1278 | 
  | 
    n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells ); | 
| 1279 | 
  | 
    temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) ); | 
| 1280 | 
  | 
    n_per_extra = (int)ceil( temp1 ); | 
| 1289 | 
  | 
  } | 
| 1290 | 
  | 
  else{ | 
| 1291 | 
  | 
    n_cells = (int)temp3; | 
| 1292 | 
< | 
    cellx = simnfo->box_x / temp3; | 
| 1293 | 
< | 
    celly = simnfo->box_y / temp3; | 
| 1294 | 
< | 
    cellz = simnfo->box_z / temp3; | 
| 1292 | 
> | 
    cellx = simnfo->boxLx / temp3; | 
| 1293 | 
> | 
    celly = simnfo->boxLy / temp3; | 
| 1294 | 
> | 
    cellz = simnfo->boxLz / temp3; | 
| 1295 | 
  | 
  } | 
| 1296 | 
  | 
 | 
| 1297 | 
  | 
  current_mol = 0; |