#!/usr/local/bin/bash
#
# Multiwfn start-up script, written by Xing (stecue@gmail.com)
# Copyright (C) 2017 Multiwfn developers and gMultiwfn project
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#

MW_PROG="/usr/local/lib/gMultiwfn/Multiwfn"

export Multiwfnpath="$HOME/.config/gMultiwfn"
if [ ! -d "$Multiwfnpath" ]
then
    mkdir -p "$Multiwfnpath"
fi
fsettings="${Multiwfnpath}/settings.ini"
if [ ! -s "${fsettings}" ]
then
    echo "No settings.ini found! Creating a new config file with default optons in"
    echo "          ${Multiwfnpath}"
    echo -n "..."
cat > "${fsettings}" <<EOF
// This file store parameters used by Multiwfn and can be modified by user. If this file is missing, Multiwfn will use default parameters.
// Texts followed by "//" are comment, those after the first semicolon are allowed options or valid data range, the content after the second semicolon is default value.
// Before and after each value a space is needed, don't insert any space between variable name and equal sign.


//Below are the parameters that can affect calculation results
  iuserfunc= 0 // Determine which user defined function to be used, see the end of section 2.7 of the manual for detail. You can also use option 2 in main function 1000 to set this parameter; (-inf,inf) integer; 0
  refxyz= 0,0,0 // The X Y and Z coordinates of reference point for calculating exchange-correlation density, correlation hole, correlation factor and source function, the unit is Bohr. You can also set the value in main function 1000 (not documented); (-inf,inf) float; 0,0,0
  iDFTxcsel= 84 // Used to select DFT exchange-correlation functional, see Section 2.7 of the manual for explanation; [0,inf) integer; 84
  paircorrtype= 3 // Determine which type of correlation effect will be taken into consideration in calculation of exchange-correlation density, correlation hole and correlation factor, 1=Only exchange part, 2=Only Coulomb part, 3=Both exchange and Coulomb parts; 1/2/3; 3
  pairfunctype= 1 // Determine which function and which spin will be calculated by real space function 17, see corresponding description in Section 2.6 for detail. =1/2: correlation hole for alpha/beta electrons. =4/5: correlation factor for alpha/beta electrons. =7/8: Exchange-correlation density for alpha/beta electrons. =10/11/12: Pair density for alpha/beta/total pair; 1/2/4/5/7/8/10/11/12; 1
  iautointgrid= 1 // 0/1= Disable/enable Multiwfn to change radpot and sphpot value in certain function adaptively to get good balence between accuracy and cost; 0/1; 1
  radpot= 75 // Number of points for Becke numerical quadrature in radial; [1,inf) integer; 75
  sphpot= 434 // Number of points for Becke numerical quadrature on spherical surface; 110/170/230/266/302/434/590/770/974/1454/1730/2030/2354/2702/3074/3470/3890/4334/4802/5294/5810; 434
  radcut= 10 // Unit is bohr. For saving time in Becke's numerical quadrature, discard all the points if the distance between which and reference atom is longer than radcut. If radcut=0, no points will be discarded. Generally default value is safe; [0,inf) float; 10
  expcutoff= -40 // When evaluating exp(x) term, if x is more negative than expcutoff, then this evaluation will be skipped for saving time. You can set it to 1 to disable this threatment; (-inf,inf) float; -40
  espprecutoff= 0 // The cutoff of "prefactor" for electrostatic potential evaluation, increase this parameter will speed up ESP evaluation, but decrease calculation accuracy; [0,inf) float; 0
  RDG_maxrho= 0.05 // When calculate reduced density gradient, if electron density >= RDG_maxrho, the function value will be set to 100. If RDG_maxrho= 0.0, this option will be disabled; [0,inf) float; 0.05
  RDGprodens_maxrho= 0.1 // The same as RDG_maxrho, but for the one with promolecule approximation; [0,inf) float; 0.1
  ELF_addminimal= 1 // If add a minial value(1D-5) to D(r) for fixing the drawback in the definition of Becke's ELF; 0/1=No/Yes; 1
  ELFLOL_type= 0 // =0: Use ELF/LOL defined by Becke, =1 Use ELF/LOL defined by Tsirelson, =2 Use special ELF/LOL definition (see manual) =3: Only get D/D0 term of ELF; 0/1/2; 0
  ELFLOL_cut= 0.0 // If this value is set to x, then ELF(or LOL) will be zero where ELF(or LOL) is less than x; [0,1] float; 0.0
  iALIEdecomp= 0 // If =1, output contribution of each occupied orbital to total average local ionization energy in main function 1; 0/1; 0
  srcfuncmode= 1 // The mode of source function, see description of source function in Multiwfn manual; 1/2; 1
  atomdenscut= 1 // During the RDG analyses with promolecular approximation, if the distance between current point and an atom is longer than cutoff value, the atom will be ignored to reduce time consuming. =1/2/3/4: if larger than 2.5/2.2/1.8/1.5 times of atomic vdW radius then the atom will be ignored, so 1 and 4 are the most accurate and inaccurate one, respectively. =0 Don't employ the cutoff treatment; 0/1/2/3; 1
  aug1D= 1.5 // Default augment size(in bohr) for calculating function value in a line; (0,inf) float; 1.5
  aug2D= 4.5 // Default augment size(in bohr) for calculating function value in a plane; (0,inf) float; 4.5
  aug3D= 6.0 // Default augment size(in bohr) for calculating function value in a spatial scope; (0,inf) float; 6.0
  num1Dpoints= 3000 // The number of points for plotting function in a line, DOS and spectrum, if the property selected is total ESP, the value will be divided by 6 automatically to reduce computational time; [1,inf) integer; 3000
  nprevorbgrid= 120000 // The number of points to be calculated for previewing orbitals in main function 0; [1,inf) integer; 120000
  bndordthres= 0.05 // Threshold for outputting bond order; [0,inf) float; 0.05
  compthres= 0.5 // Threshold for printing orbital composition in main function 8, unit is %; [0,inf) float; 0.5
  compthresCDA= 1 // Similar to compthres, but this parameter applies to CDA module; [0,inf) float; 1
  ispheratm= 1 // If do sphericalization for free atom wavefunctions, 0/1=no/yes; 0/1; 1
  laplfac= 1 // Multiply electron density Laplacian with this value; (-inf,inf) float; 1.0
  ipolarpara= 0 // If equals 1, spin polarization parameter function will be used instead of spin density; 0/1; 0
  ishowchgtrans= 0 // =1: Show charge transfer between atoms in atomic dipole moment correction (ADC) process or in CM5 calculation, =0: don't show; 0/1; 0
  SpherIVgroup= 0 // Determine how to generate and sphericalize wavefunction of atoms in IV main group, =0: Use sp3 configuration =1: Use s2p2 configuration then rotate and duplicate occupied orbital, then equipartition occupation number. If you use "atomwfn" folder to store calculated atomic wavefunctions, don't forget to recalculate them!; 0/1; 0
  MCvolmethod= 2 // Definition of molecular volume, the value will be evaluated by Monte Carlo method. =1: Superposition of vdW sphere of atoms, =2: Encompassed by certain isosurface of density; 1/2; 2
  readEDF= 1 // If read EDF field in the wfx file produced by Gaussian to represent inner core density. 0=Don't read, 1=Read; 0/1; 1
  isupplyEDF= 2 // Control how to supply EDF field for the input file containing GTF information, of course this option is only relevant when pseudopotential basis set was involved (For .wfx file, this option is only meaningful when EDF field is not presented, or it is presented but readEDF is set to zero). =0: Don't supply EDF field, =1: Read EDF field from atomic .wfx files produced by Gaussian, =2: Supply EDF field from built-in EDF library (edflib.f90); integer; 2
  idelvirorb= 1 // =1: At some stages automatically delete virtual orbitals higher than LUMO+10 to speed up calculation of real space functions, =0: don't delete any orbital automatically; 0/1; 1
  ifchprog= 1 // =1: The .fch is produced by Gaussian, =2: The .fch is produced by Q-Chem; 1/2; 1
  ishowptESP= 1 // =1: Also show electrostatic potential when showing real space functions at a given point (or at a critical point), =0 Don't show it for saving computational time; 0/1; 1
  imolsurparmode= 1 // =1: Use Voronoi-like partition of molecular surface, =2: Partition the molecular surface based on square root of Bondi atomic vdW radius; 1/2; 1
  steric_addminimal= 0.0001 // Add which value to the electron density terms in the expression of steric potential, force and charge; (-inf,inf) float; 0.0001
  steric_potcutrho= 0 // If electron density at a point is lower than this value, then steric potential will be set to a constant "steric_potcons", and steric charge will be 0. This value has different meaning for damped form of Steric potential/force, see manual; [0,inf) float; 0.0
  steric_potcons= 0 // See above; (-inf,inf) float; 0.0
  NICSnptlim= 8000 // In option 4 of main function 200, the maximal number of atoms+Bq in each Gaussian NMR task; [1,inf); 8000
  iplaneextdata= 0 // In plane map plotting module (main function 4), if directly loading data from external file instead of calculating them by Multiwfn. 0=No, 1=Yes; 0/1; 0
  igenP= 1 // =1: Generate or load density matrix when loading .fch/.molden/.gms files. =0: Don't do this for saving memory, but some functions will not work; 0/1; 1
  igenDbas= 0 // =1: Generate dipole moment integral matrix between basis functions when loading .fch/.molden/.gms files. =0: Don't generate it for saving time; 0/1; 0
  igenMagbas= 0 // =1: Generate magnetic moment integral matrix between basis functions when loading .fch/.molden/.gms files. =0: Don't generate it for saving time; 0/1; 0
  iloadasCart= 0 // Convert spherical harmonic type of basis functions (if any) to Cartesian type when loading .fch or .molden file; 0/1; 0


//Below are the parameters involved in plotting
  symbolsize= 8 // Symbol size for drawing scatter graph; [1,inf) integer; 8
  pleatmlabsize= 50 // Size of atomic labels when plotting plane graph; [1,inf) integer; 50
  disshowlabel= 0.5 // When drawing color-filled/contour/gradient/vector field map, if the distance between an atom nucleus/critical point and the plane of interest is less than this value (bohr), the atom/critical point will be labelled on the graph; [0,inf) float; 0.5
  iatom_on_contour_far= 0 // When showing atom labels in plane map, if also show the labels of the atoms that beyond "disshowlabel" as light face type; 0/1=Yes/No; 0/1; 0
  iatmlabtype= 1 // When showing atom labels in plane map, =1 only plot element symbol, =2 only plot atom index, =3 plot both; 1/2/3; 1
  iatmlabtype3D= 3 // When showing atom labels in 3D map, =1 only plot element symbol, =2 only plot atom index, =3 plot both; 1/2/3; 3
  graphformat= png // Format for outputting graph file; ps/eps/pdf/wmf/gif/tiff/bmp; png
  graph1Dsize= 1600,1000 // Width and height of outputted graph file for line plotting; [1,inf) integer; 1600,1000
  graph2Dsize= 1500,1200 // Width and height of outputted graph file for 2D plotting. If you want to adjust this value, you'd better ensure that the ratio of of the two values is 5:4, e.g. 3000:2400 and 1000:800; [1,inf) integer; 1500,1200
  graph3Dsize= 1400,1400 // Width and height of outputted graph file for 3D plotting; [1,inf) integer; 1400,1400
  numdigxyz= 2,2,3 // The number of digits after the decimal point of the labels on X, Y and Z axes of plane plots; [0,inf) integer; 2,2,3
  numdiglinexy= 3,3 // Similar to numdigxyz, but for line plots; [0,inf) integer; 3,3
  numdigctr= 3 // The number of digits after the decimal point of label on contour lines; [0,inf) integer; 3
  fillcoloritpxy= 5,5 // The number of interpolation steps between grids in X and Y directions when drawing filled color map, increase it will make the graph prettier but plotting speed will be slow down; [1,inf) integer; 5,5
  inowhiteblack= 0 // If set to 1, when plotting color-filled map, the regions where the value is higher than and lower than the color scale will not be colored as white and black, respectively; 0/1; 0
  isurfstyle= 2 // Use which style to plot interbasin surface in 3D graph. =1 use a bunch of paths, =2 use solid surface consists of triangle tiles, =3 use solid surface consists of cylinders; 1/2/3; 2
  bondRGB= 0.1,1.0,0.1 //The red, green and blue components of the color of bonds. 0,0,0 corresponds to black, 1,1,1 corresponds to white; [0,1] float; 0.1,1.0,0.1
  atmlabRGB= 0.0,0.0,0.0 //The red, green and blue components of the atomic labels in 3D plot. 0,0,0 corresponds to black, 1,1,1 corresponds to white; [0,1] float; 0.0,0.0,0.0
  CP_RGB= 0.72,0,0.72, 1,0.5,0, 1,1,0, 0,1,0 // RGB color for (3,-3), (3,-1), (3,+1) and (3,+3) type critical points; [0,1] float; 0.72,0,0.72, 1,0.5,0, 1,1,0, 0,1,0
  atmcolorfile= none // The path of the file defining the color of atom spheres in 3D plots. "examples\element_color.txt" is a template file, in which red, green and blue components (between 0.0~1.0) are defined for every element, you can alter color for specific element by changing corresponding values. If write "none", default colors will be used.


//Below are the parameters about system environment
  nthreads= 0 //How many threads are used in parallel mode, if equals 1, parallel mode will be disabled; [1,inf) integer; 4
  ompstacksize= 80000000 //For parallel implementation in windows version, this variable sets stacksize for each thread, unit is in byte. If available memory is insufficient, this value should be decreased. If your task size is huge and crashes during calculation, try to increase ompstacksize. ompstacksize doesn't affect Linux or Mac OS X version, for which you should modify KMP_STACKSIZE environment variable instead; [0,inf) integer; 80000000
  gaupath= "/sob/g16/g16" //The path of Gaussian executable file (not the GUI version such as g09w.exe), quotation marks are needed! Must be less than 80 characters.
  isilent= 0 // If pop up graph or GUI interface immediately after calculation of line, plane... is finished or changing plotting parameter. This option is useful for batch processing, so that user needn't to manually close graph or GUI interface that popped up; 0/1; 0
  isys= 2 // Used to determine current operation system, 1=windows 2=Linux 3=MacOS; 1/2/3; 1
  imodlayout= 0 // =1: For Windows version, use different layout in main function 0 to avoid truncation of orbital list; 0/1; 0
  outmedinfo= 0 // Output some intermediate informations to file, only for debugging; 0/1; 0
  iwfntmptype= 1 // When atomic wavefunction files are needed to be generated, =1 delete the old "wfntmp" folder (if presents) and then use the newly built "wfntmp" folder, =2 Build and use "wfntmp" folder with different suffix, e.g. wfntmp0001, wfntmp0002... so that you can simultaneously run multiple instances without conflict; 1/2; 1
  ispecial= 0 // Sob only knows. 1=For RCY, 2=For shubin's 2nd project, Hirshfeld density as reference, or obtain Renyi entropy; 0/1/2; 0

//The last opened file name (There must be a space line in present file)
lastfile= D:\raman\MPBA_DA.log
EOF
    echo "Done!"
    echo "Starting Multiwfn with OMP_NUM_THREADS enabled..."
    if [ x$OMP_NUM_THREADS = 'x' ]
    then
        nCPUs=`sysctl -n hw.ncpu`
        echo "Warning! OMP_NUM_THREADS not set. As many as $nCPUs threads might be used."
    fi
    echo ""
    ${MW_PROG} "$@"
else
    nCPUs=`grep "nthreads" "${fsettings}"|sed -e 's/.*nthreads[ =]*\([0-9]*\).*/\1/g'`
    if [ x$nCPUs = 'x' ]
    then
        echo "Please check your settings.ini file. Wrong nthreads setting."
        exit
    fi
    if [ $nCPUs -eq 0 ]
    then
        echo "Starting Multiwfn with OMP_NUM_THREADS enabled..."
        if [ x$OMP_NUM_THREADS = 'x' ]
        then
            nCPUs=`sysctl -n hw.ncpu`
            echo "Warning! OMP_NUM_THREADS not set. As many as $nCPUs threads might be used."
        fi
    fi
    echo ""
    ${MW_PROG} "$@"
fi
