//NAB program to provide single point 3D-RISM calculations from the commandline

//All 3D-RISM options are supported except saveprogress. 
//saveprogress will be enabled when a better format is introduced.

//Tyler Luchko 2010/12

#include "nab_netcdf.h"

//global variables to reduce the number of parameters
string name;      //program name

struct mdOptions { 
  string pdb, prmtop, rst, traj;};
struct rismOptions {
  string xvv, guv, cuv, huv, uuv, asymp, quv, chgdist, volfmt;
  int closureOrder, asympcorr;
  float buffer, solvcut;
  float grdspcX, grdspcY, grdspcZ, solvboxX, solvboxY, solvboxZ;
  int ngX, ngY, ngZ;
  float mdiis_del, mdiis_restart;
  int mdiis_nvec, mdiis_method;
  int maxstep, npropagate;
  int centering, zerofrc, apply_rism_force, polarDecomp;
  int ntwrism, verbose, progress, saveprogress;};

/////////////////
//MD options
/////////////////
struct mdOptions mdOpt;
 
/////////////////
//3D-RISM options
/////////////////
struct rismOptions rismOpt;
//no arrays in structs so 'closure' and 'tolerance' are special cases
string closure[ dynamic ];
float tolerance[ dynamic ];

///////////////////////////////////////////////////////////////////////////////
//Prints an error message to stderr
//IN:
//   mytaskid : MPI process id
//   message  : message to output to user
///////////////////////////////////////////////////////////////////////////////
int error(string message) c;
int error(string message){
  if(mytaskid==0){
    return fprintf(stderr,"ERROR:  %s\n", message);
  }
  return -1;
};

///////////////////////////////////////////////////////////////////////////////
//sets default options
//IN:
//   mdOpt : mdOptions structure
//   rismOpt : rismOptions structure
///////////////////////////////////////////////////////////////////////////////
int setDefaults(struct mdOptions mdOpt, struct rismOptions rismOpt) c;
int setDefaults(struct mdOptions mdOpt, struct rismOptions rismOpt){
  mdOpt.pdb                = NULL;
  mdOpt.prmtop             = NULL;
  mdOpt.rst                = NULL;
  mdOpt.traj               = NULL;
  rismOpt.xvv              = NULL;
  rismOpt.guv              = NULL;
  rismOpt.huv              = NULL;
  rismOpt.cuv              = NULL;
  rismOpt.uuv              = NULL;
  rismOpt.asymp            = NULL;
  rismOpt.quv              = NULL;
  rismOpt.chgdist          = NULL;
  rismOpt.volfmt           = "dx";
  closure[ 1 ]             = "kh";
  rismOpt.closureOrder     = 1;
  rismOpt.asympcorr        = 1;
  rismOpt.buffer           = 14;
  rismOpt.solvcut          = -1;
  rismOpt.grdspcX          = 0.5;
  rismOpt.grdspcY          = 0.5;
  rismOpt.grdspcZ          = 0.5;
  rismOpt.ngX              = 0;
  rismOpt.ngY              = 0;
  rismOpt.ngZ              = 0;
  rismOpt.solvboxX         = 0;
  rismOpt.solvboxY         = 0;
  rismOpt.solvboxZ         = 0;
  tolerance[ 1 ]           = 1e-5;
  rismOpt.mdiis_del        = 0.7;
  rismOpt.mdiis_restart    = 10;
  rismOpt.mdiis_nvec       = 5;
  rismOpt.mdiis_method     = 2;
  rismOpt.maxstep          = 10000;
  rismOpt.npropagate       = 5;
  rismOpt.centering        = 1;
  rismOpt.zerofrc          = 1;
  rismOpt.apply_rism_force = 0;
  rismOpt.polarDecomp      = 0;
  rismOpt.ntwrism          = 0;
  rismOpt.verbose          = 0;
  rismOpt.progress         = 1;
  rismOpt.saveprogress     = 0;
  return 0;
};

///////////////////////////////////////////////////////////////////////////////
//Prints option values
//IN:
//   mdOpt : mdOptions structure
//   rismOpt : rismOptions structure
///////////////////////////////////////////////////////////////////////////////
int printOptions(struct mdOptions mdOpt, struct rismOptions rismOpt) c;
int printOptions(struct mdOptions mdOpt, struct rismOptions rismOpt){
  setDefaults(mdOpt, rismOpt);
  
  fprintf(stderr,"\n");
  fprintf(stderr,"%-20s %s\n", "Key", "Default");
  fprintf(stderr,"%-20s %s\n", "---", "-------");

  fprintf(stderr,"%-20s %s\n", "--pdb", mdOpt.pdb);
  fprintf(stderr,"%-20s %s\n", "--prmtop", mdOpt.prmtop);
  fprintf(stderr,"%-20s %s\n", "--rst", mdOpt.rst);
  fprintf(stderr,"%-20s %s\n", "-y|--traj", mdOpt.traj);
  fprintf(stderr,"%-20s %s\n", "--xvv", rismOpt.xvv);
  fprintf(stderr,"%-20s %s\n", "--guv", rismOpt.guv);
  fprintf(stderr,"%-20s %s\n", "--huv", rismOpt.huv);
  fprintf(stderr,"%-20s %s\n", "--cuv", rismOpt.cuv);
  fprintf(stderr,"%-20s %s\n", "--uuv", rismOpt.uuv);
  fprintf(stderr,"%-20s %s\n", "--asymp", rismOpt.asymp);
  fprintf(stderr,"%-20s %s\n", "--quv", rismOpt.quv);
  fprintf(stderr,"%-20s %s\n", "--chgdist", rismOpt.chgdist);
  fprintf(stderr,"%-20s %s\n", "--volfmt", rismOpt.volfmt);
  fprintf(stderr,"%-20s %s\n", "--closure", closure[1]);
  fprintf(stderr,"%-20s %d\n", "--closureorder", rismOpt.closureOrder);
  fprintf(stderr,"%-20s %i\n", "--asympcorr", rismOpt.asympcorr);
  fprintf(stderr,"%-20s %f\n", "--buffer", rismOpt.buffer);
  fprintf(stderr,"%-20s %f\n", "--solvcut", rismOpt.solvcut);
  fprintf(stderr,"%-20s %f,%f,%f\n", "--grdspc", rismOpt.grdspcX, 
	  rismOpt.grdspcY, rismOpt.grdspcZ);
  fprintf(stderr,"%-20s %d,%d,%d\n", "--ng", rismOpt.ngX, rismOpt.ngY,
	  rismOpt.ngZ);
  fprintf(stderr,"%-20s %f,%f,%f\n", "--solvbox", rismOpt.solvboxX, 
	  rismOpt.solvboxY, rismOpt.solvboxZ);
  fprintf(stderr,"%-20s %f\n", "--tolerance", tolerance[1]);
  fprintf(stderr,"%-20s %f\n", "--mdiis_del", rismOpt.mdiis_del);
  fprintf(stderr,"%-20s %f\n", "--mdiis_restart", rismOpt.mdiis_restart);
  fprintf(stderr,"%-20s %d\n", "--mdiis_nvec", rismOpt.mdiis_nvec);
  //  fprintf(stderr,"%-20s = %d\n", "--mdiis_method", rismOpt.mdiis_method);
  fprintf(stderr,"%-20s %d\n", "--maxstep", rismOpt.maxstep);
  fprintf(stderr,"%-20s %d\n", "--npropagate", rismOpt.npropagate);
  fprintf(stderr,"%-20s %d\n", "--centering", rismOpt.centering);
  //  fprintf(stderr,"%-20s = %d\n", "--zerofrc", rismOpt.zerofrc);
  //  fprintf(stderr,"%-20s = %d\n", "--apply_rism_force", rismOpt.apply_rism_force);
  fprintf(stderr,"%-20s %d\n", "--polarDecomp", rismOpt.polarDecomp);
  fprintf(stderr,"%-20s %d\n", "--verbose", rismOpt.verbose);
  //  fprintf(stderr,"%-20s = %d\n", "--progress", rismOpt.progress);
  //  fprintf(stderr,"%-20s = %d\n", "--saveprogress", rismOpt.saveprogress);
  return 0;
};

///////////////////////////////////////////////////////////////////////////////
//Presents usage information for the program to stderr and exits with status 1.
//IN:
//   name : name of the program
//   mytaskid : MPI process id
///////////////////////////////////////////////////////////////////////////////
int usage() c;
int usage(){
  string empty;
  int i;
  if(mytaskid==0){
    fprintf(stderr,"USAGE: %s --pdb pdbfile --prmtop prmtopfile [--rst rstfile] [-y|--traj netCDFfile]\n",name);
    for(i=0; i< length(name); i++){
      empty = empty+" ";
    }
    fprintf(stderr,"       %s --xvv Xvv_filename [--guv Guv_rootname] [--cuv Cuv_rootname] [--huv Huv_rootname]\n",empty);
    fprintf(stderr,"       %s [--uuv Uuv_rootname] [--asymp asymp_rootname]\n",empty);
    fprintf(stderr,"       %s [--quv Quv_rootname] [--chgdist chgdist_rootname]\n",empty);
    fprintf(stderr,"       %s [--volfmt volume_format]\n",empty);
    fprintf(stderr,"       %s [--closure kh|hnc|pse(1|2|...)[ closure2[ ...]]]\n",empty);
    fprintf(stderr,"       %s [--[no]asympcorr] [--buffer distance] [--solvcut distance]\n",empty);
    fprintf(stderr,"       %s [--grdspc dx,dy,dz] [--ng nx,ny,nz] [-solvbox lx,ly,lz]\n",empty);
    fprintf(stderr,"       %s [--tolerance tol1[ tol2 [...]] [--mdiis_del step_size]\n",empty); 
    fprintf(stderr,"       %s [--mdiis_restart threshold] [--mdiis_nvec vectors]\n",empty);
    fprintf(stderr,"       %s [--maxstep 1000] [--npropagate #_old_solutions]\n",empty);
    //    fprintf(stderr,"       %s [--centering -2|1|0|1|2] [--[no]zerofrc] [--[no]apply_rism_force]\n",empty);
    fprintf(stderr,"       %s [--[no]polarDecomp] [--centering -4..4]\n",empty);
    fprintf(stderr,"       %s [--verbose 0|1|2]\n",empty);
    printOptions(mdOpt, rismOpt);
  }
  exit(1);
  return 0;
}; 

///////////////////////////////////////////////////////////////////////////////
//tests if value is a key in key/value pair.
//IN:
//   key   : the command line argument that is a key
//   value : the command line argument that should be a value
//SIDE EFFECT:
//   if 'value' starts with a '-' then an error message is printed followed
//   by calling the usage function.
///////////////////////////////////////////////////////////////////////////////
int testValue(string key, string value) c;
int testValue(string key, string value){
  if(value=~"^-" && value !~"^-[0-9]*$"){
    error(sprintf("'%s' must be followed by a value.\n",key)); 
    usage(); 
  }
  return 0;
};

///////////////////////////////////////////////////////////////////////////////
//counts the number of Comma Separated Values in a string
//IN:
//   str : string
//OUT:
//   number of comma separated values.  Includes blank values so this can be 
//   thought of as number of commas plus one.
///////////////////////////////////////////////////////////////////////////////
int numCSV(string str) c;
int numCSV(string str){
  int pos0, pos1, num;
  num = 0;
  pos0 = 0;
  pos1 = 0;
  while((pos0 = index(substr(str,pos1+1,length(str)),",")) !=0){
    num++;
    pos1 += pos0;
    if(pos1+1 > length(str)) {break;}
  }
  return ++num;
};

///////////////////////////////////////////////////////////////////////////////
//Final check to ensure options are correctly set.
//IN:
//   mdOpt : mdOptions structure
//   rismOpt : rismOptions structure
///////////////////////////////////////////////////////////////////////////////
int checkOptions(struct mdOptions mdOpt, struct rismOptions rismOpt) c;
int checkOptions(struct mdOptions mdOpt, struct rismOptions rismOpt){
  if(mdOpt.pdb==NULL){
    error("a PDB file is required\n");
    usage();
  }
  if(mdOpt.prmtop==NULL){
    error("a PRMTOP file is required\n");
    usage();
  }
  if(rismOpt.xvv==NULL){
    error("an XVV file is required\n");
    usage();
  }
  return 0;
};

///////////////////////////////////////////////////////////////////////////////
//Main
///////////////////////////////////////////////////////////////////////////////
molecule m;
float 	p_xyz[dynamic], f_xyz[dynamic], v_xyz[dynamic];
struct AmberNetcdf nc;
float 	time, dgrad, fret;
int 	i, j; 
string value, key, tmpstr;
int tempInt;
float ftmp;
file asciiTraj;
int iframe, ipos, trajDone;
int nclosure,ntolerance;
//set global variables
name = argv[1];

//closure name is set in setDefaults()
allocate closure[1];
nclosure=1;
//tolerance is set in setDefaults()
allocate tolerance[1];
ntolerance=1;
setDefaults(mdOpt, rismOpt);
////////////////////
//read command line options.  (This should be a separate function but I
//can't pass arrays of strings.)
////////////////////

i=2;
while(i<=argc){
  key = argv[i];
  value = "-";
  if(key=="--printDefaults"){
    setDefaults(mdOpt, rismOpt);
    printOptions(mdOpt,rismOpt);
    exit(0);
  }else if(key=="--pdb"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    mdOpt.pdb=value;
  }else if(key=="--rst"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    mdOpt.rst=value;
  }else if(key=="--traj" || key=="-y"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    mdOpt.traj=value;
  }else if(key=="--prmtop"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    mdOpt.prmtop=value;
  }else if(key=="--xvv"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    rismOpt.xvv=value;
  }else if(key=="--guv"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    rismOpt.guv=value;
    rismOpt.ntwrism=1;
  }else if(key=="--cuv"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    rismOpt.cuv=value;
    rismOpt.ntwrism=1;
  }else if(key=="--huv"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    rismOpt.huv=value;
    rismOpt.ntwrism=1;
  }else if(key=="--uuv"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    rismOpt.uuv=value;
    rismOpt.ntwrism=1;
  }else if(key=="--asymp"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    rismOpt.asymp=value;
    rismOpt.ntwrism=1;
  }else if(key=="--quv"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    rismOpt.quv=value;
    rismOpt.ntwrism=1;
  }else if(key=="--chgdist"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    rismOpt.chgdist=value;
    rismOpt.ntwrism=1;
  }else if(key=="--volfmt"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(value !~ "dx" && value !~ "xyzv" ){
      error("--volfmt must be one of dx or xyzv\n");
      usage();
    }
    rismOpt.volfmt=value;
  }else if(key=="--closure"){
    //find the number of closures listed
    j=i;
    while (j+1 <=argc && argv[j+1]!~"^-"){
      j++;
    }
    nclosure = j-i;
    /* allocate memory */
    deallocate(closure);
    allocate(closure[nclosure]);
    for(j=1;j<=nclosure;j++){
      i++;
      if(i<=argc) {
        value = argv[i];
      }
      testValue(key,value);
      if(value !~ "[kK][hH]" && value !~ "[hH][nN][cC]" && value !~ "[pP][sS][eE]" ){
        error("--closure must be one of kh, hnc, or pse\n");
        usage();
      }
      closure[j]=value;
    }
  }else if(key=="--closureorder"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%i",rismOpt.closureOrder) != 1 
       ||rismOpt.closureOrder<1 ){
      error("--closureOrder takes an integer > 0\n");
      usage();
    }
  }else if(key=="--asympcorr"){
    rismOpt.asympcorr=1;
  }else if(key=="--noasympcorr"){
    rismOpt.asympcorr=0;
  }else if(key=="--buffer"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%lf",rismOpt.buffer) != 1){
      error("--buffer takes a float\n");
      usage();
    }
  }else if(key=="--solvcut"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%lf",rismOpt.solvcut) != 1
       || rismOpt.solvcut <=0){
      error("--solvcut takes a float > 0\n");
      usage();
    }
  }else if(key=="--grdspc"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    tempInt=numCSV(value);
    if(tempInt ==3){
      if(sscanf(value,"%lf,%lf,%lf",
		rismOpt.grdspcX,rismOpt.grdspcY,rismOpt.grdspcZ) != 3
	 || rismOpt.grdspcX <= 0 || rismOpt.grdspcY <= 0
	 || rismOpt.grdspcZ <= 0){
	error("--grdspc takes floats > 0\n");
	usage();
      }
    }else if(tempInt ==1){
      if(sscanf(value,"%lf",rismOpt.grdspcX) != 1
	 || rismOpt.grdspcX <= 0){
	error("--grdspc takes floats\n");
	usage();
      }
      rismOpt.grdspcY = rismOpt.grdspcX;
      rismOpt.grdspcZ = rismOpt.grdspcX;
    }else{
      error("--grdspc only takes one or three values\n");
      usage();
    }
  }else if(key=="--ng"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    tempInt=numCSV(value);
    if(tempInt ==3){
      if(sscanf(value,"%i,%i,%i",
		rismOpt.ngX,rismOpt.ngY,rismOpt.ngZ) != 3 
	 || rismOpt.ngX <= 0 || rismOpt.ngY <= 0 || rismOpt.ngZ <= 0){
	error("--ng takes integers > 0\n");
	usage();
      }
    }else if(tempInt ==1){
      if(sscanf(value,"%i",rismOpt.ngX) != 1
	 || rismOpt.ngX <= 0){
	error("--ng takes integers\n");
	usage();
      }
      rismOpt.ngY = rismOpt.ngX;
      rismOpt.ngZ = rismOpt.ngX;
    }else{
      error("--ng only takes one or three values\n");
      usage();
    }
  }else if(key=="--solvbox"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    tempInt=numCSV(value);
    if(tempInt ==3){
      if(sscanf(value,"%lf,%lf,%lf",
		rismOpt.solvboxX,rismOpt.solvboxY,rismOpt.solvboxZ) != 3
	 || rismOpt.solvboxX <=0 || rismOpt.solvboxY <= 0 
	 || rismOpt.solvboxZ <= 0){
	error("--solvbox takes floats > 0\n");
	usage();
      }
    }else if(tempInt ==1){
      if(sscanf(value,"%lf",rismOpt.solvboxX) != 1
	 || rismOpt.solvboxX <=0){
	error("--solvbox takes floats > 0\n");
	usage();
      }
      rismOpt.solvboxY = rismOpt.solvboxX;
      rismOpt.solvboxZ = rismOpt.solvboxX;
    }else{
      error("--solvbox only takes one or three values\n");
      usage();
    }
  }else if(key=="--tolerance"){
    //find the number of tolerances listed
    j=i;
    while (j+1 <=argc && argv[j+1]!~"^-"){
      value = argv[j+1];
      j++;
    }
    ntolerance = j-i;
    /* allocate memory */
    deallocate(tolerance);
    allocate(tolerance[ntolerance]);
    for(j=1;j<=ntolerance;j++){
      i++;
      if(i<=argc) {
        value = argv[i];
      }
      testValue(key,value);
      if(sscanf(value,"%lf",tolerance[j]) != 1
         || tolerance[j] <=0){
        error("--tolerance takes a float > 0\n");
        usage();
      }
    }
  }else if(key=="--mdiis_del"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%lf",rismOpt.mdiis_del) != 1
       || rismOpt.mdiis_del <=0 || rismOpt.mdiis_del >2 ){
      error("--mdiis_del takes a float > 0 and <= 2\n");
      usage();
    }
  }else if(key=="--mdiis_restart"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%lf",rismOpt.mdiis_restart) != 1
       || rismOpt.mdiis_restart <=0){
      error("--mdiis_del takes a float > 0\n");
      usage();
    }
  }else if(key=="--mdiis_nvec"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%i",rismOpt.mdiis_nvec) != 1 
       ||rismOpt.mdiis_nvec<1 ){
      error("--mdiis_nvec takes an integer > 0\n");
      usage();
    }
  }else if(key=="--mdiis_method"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%i",rismOpt.mdiis_method) != 1 
       || rismOpt.mdiis_method<0 || rismOpt.mdiis_method>2 ){
      error("--mdiis_method must be 0, 1 or 2\n");
      usage();
    }
  }else if(key=="--maxstep"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%i",rismOpt.maxstep) != 1 
       ||rismOpt.maxstep<1 ){
      error("--maxstep takes an integer > 0\n");
      usage();
    }
  }else if(key=="--npropagate"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%i",rismOpt.npropagate) != 1 
       ||rismOpt.npropagate<0 ){
      error("--npropagate takes an integer >= 0\n");
      usage();
    }
  }else if(key=="--centering"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%i",rismOpt.centering) != 1 
       ||rismOpt.centering<-4 || rismOpt.centering > 4 ){
      error("--centering must be between -4 and 4\n");
      usage();
    }
//  }else if(key=="--zerofrc"){
//    rismOpt.zerofrc=1;
//  }else if(key=="--nozerofrc"){
//    rismOpt.zerofrc=0;
//  }else if(key=="--apply_rism_force"){
//    rismOpt.apply_rism_force=1;
//  }else if(key=="--noapply_rism_force"){
//    rismOpt.apply_rism_force=0;
  }else if(key=="--polarDecomp"){
    rismOpt.polarDecomp=1;
  }else if(key=="--nopolarDecomp"){
    rismOpt.polarDecomp=0;
  }else if(key=="--verbose"){
    i++;
    if(i<=argc) {
      value = argv[i];
    }
    testValue(key,value);
    if(sscanf(value,"%i",rismOpt.verbose) != 1 
       || rismOpt.verbose < 0 || rismOpt.verbose > 2 ){
      error("--verbose must be 0, 1 or 2\n");
      usage();
    }
  }else if(key=="--progress"){
    rismOpt.progress=1;
  }else if(key=="--noprogress"){
    rismOpt.progress=0;
  }else if(key=="--saveprogress"){
    rismOpt.saveprogress=1;
  }else if(key=="--nosaveprogress"){
    rismOpt.saveprogress=0;
  }else{
    error(sprintf("unknown option: '%s'\n", key));
    usage();
  }
  i++;
 }

checkOptions(mdOpt,rismOpt);

////////////////////
//Run 3D-RISM
////////////////////
m = getpdb(mdOpt.pdb);
readparm(m, mdOpt.prmtop);
allocate p_xyz[3*m.natoms];
allocate f_xyz[3*m.natoms];
allocate v_xyz[3*m.natoms];
mm_options(sprintf("cut=%e",sqrt(10e37)));

mm_options(sprintf("rism=1, ntpr_rism=1, xvvfile=%s",rismOpt.xvv));
if(rismOpt.guv){
  mm_options(sprintf("guvfile=%s",rismOpt.guv));
 }
if(rismOpt.cuv){
  mm_options( sprintf("cuvfile=%s",rismOpt.cuv));
 }
if(rismOpt.huv){
  mm_options(sprintf("huvfile=%s",rismOpt.huv));
 }
if(rismOpt.uuv){
  mm_options(sprintf("uuvfile=%s",rismOpt.uuv));
 }
if(rismOpt.asymp){
  mm_options(sprintf("asympfile=%s",rismOpt.asymp));
 }
if(rismOpt.quv){
  mm_options(sprintf("quvfile=%s",rismOpt.quv));
 }
if(rismOpt.chgdist){
  mm_options(sprintf("chgdistfile=%s",rismOpt.chgdist));
 }
mm_options(sprintf("volfmt=%s",rismOpt.volfmt));
tmpstr="";
if( nclosure >= 1){
  tmpstr = sprintf("closure=%s",closure[1]);
 }
for(i=2;i<=nclosure;i++){
  tmpstr = sprintf("%s,%s",tmpstr,closure[i]);
}
mm_options(sprintf("%s, closureOrder=%i",tmpstr, rismOpt.closureOrder));
mm_options(sprintf("buffer=%g, solvcut=%g",rismOpt.buffer,rismOpt.solvcut));
mm_options(sprintf("grdspcx=%g, grdspcy=%g, grdspcz=%g",rismOpt.grdspcX,rismOpt.grdspcY,rismOpt.grdspcZ));
mm_options(sprintf("solvboxx=%g, solvboxy=%g, solvboxz=%g",rismOpt.solvboxX,rismOpt.solvboxY,rismOpt.solvboxZ));
mm_options(sprintf("ngx=%d, ngy=%d, ngz=%d",rismOpt.ngX,rismOpt.ngY,rismOpt.ngZ));
tmpstr="";
if( ntolerance >= 1){
  tmpstr = sprintf("tolerance=%g",tolerance[1]);
 }
for(i=2;i<=ntolerance;i++){
  tmpstr = sprintf("%s,%g",tmpstr,tolerance[i]);
}
mm_options(sprintf("%s, mdiis_del=%g",tmpstr,rismOpt.mdiis_del));
mm_options(sprintf("mdiis_restart=%g",rismOpt.mdiis_restart));
mm_options(sprintf("mdiis_nvec=%d, mdiis_method=%d",rismOpt.mdiis_nvec,rismOpt.mdiis_method));
mm_options(sprintf("maxstep=%d, npropagate=%d",rismOpt.maxstep,rismOpt.npropagate));
mm_options(sprintf("centering=%d, zerofrc=%d, apply_rism_force=%d",rismOpt.centering,rismOpt.zerofrc,rismOpt.apply_rism_force));
mm_options(sprintf("polarDecomp=%d",rismOpt.polarDecomp));
mm_options(sprintf("ntwrism=%d, verbose=%d, progress=%d",rismOpt.ntwrism,rismOpt.verbose,rismOpt.progress));
mm_options(sprintf("asympCorr=%d",rismOpt.asympcorr));
mm_options(sprintf("saveprogress=%d",rismOpt.saveprogress));
mme_init( m, NULL, "::Z", p_xyz, NULL);

if(mdOpt.traj){
  //try NetCDF first
    if(netcdfLoad( nc, mdOpt.traj)==0){
      netcdfLoad( nc, mdOpt.traj );
      if(mytaskid==0) printf("\nProcessing NetCDF trajectory: %s\n",mdOpt.traj);
      while ( netcdfGetNextFrame(nc,p_xyz,NULL,NULL) ) {  
        if(mytaskid==0) printf("\nFrame: %d of %d\n",nc.currentFrame, nc.ncframe);
        mme(p_xyz,f_xyz, nc.currentFrame );
      }
      netcdfClose(nc);
    }else{ 
      //Assume ASCII
      if(mytaskid==0){
        if(numtasks > 1){
          error("ASCII trajectories not supported for more than one process");
          exit(1);
        }
        printf("\nProcessing ASCII trajectory: %s\n",mdOpt.traj);
      }
      trajDone = 0;
      asciiTraj = fopen( mdOpt.traj, "r" );
      if(asciiTraj==NULL){
        error(sprintf("Failed to open '%s'",mdOpt.traj));
        exit(1);
      }
      getline(asciiTraj);
      for(iframe=1;;iframe++){
          for(ipos=1; ipos<=3*m.natoms; ipos++){
            if(fscanf(asciiTraj,"%lf",p_xyz[ipos])<1){
              trajDone=1;
              break;
            }
          }          
          if(trajDone) break;
        if(mytaskid==0){
          printf("\nFrame: %d\n",iframe);
        }
        mme(p_xyz,f_xyz, iframe );
      }
    }

 }else{
  if(mdOpt.rst){
    if(mytaskid==0) printf("\nProcessing restart file: %s\n",mdOpt.rst);
    getxv(mdOpt.rst,m.natoms,time,p_xyz,v_xyz);  
  }else{
    if(mytaskid==0) printf("\nProcessing PDB file: %s\n",mdOpt.pdb);
    setxyz_from_mol( m, NULL, p_xyz );
  }
   mme(p_xyz,f_xyz, 1 );
  
 }

if(mytaskid==0){ printf("\n3D-RISM processing complete.\n");}

mme_rism_max_memory();
mme_timer();