#!/usr/bin/env python

from __future__ import division
from numpy import *
from sympy import * 
from scipy import * 
from scipy.integrate import * 
from scipy.optimize import * 
from scipy.interpolate import * 
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from pylab import *

import numpy
import sympy
import scipy.interpolate as interp 
import matplotlib.pyplot as plt
import re
import os
import sys
import commands
import time 
import pdb

# List of solvent type, approximate methods and type of function to be used for calculation
solvents = ['OW']; functions = ['define', 'symbolic']; methods = ['spline', 'finite']
method = methods[0] # method use to evaluate the total variance 
functype = functions[0] # Use predefined function to speed up the calculation rather than using symbolic function

Cp = 4.184; # Specific heat capacity of water under constant pressure
kB = 1.380658*6.0221367/1000.0  # Boltzmann's constant (kJ/mol/K)
temp = 298; betav = Cp/(kB*temp)
rho_OW = 0.03330001
rho_HW = 2.0*rho_OW

# Value for beta and the combined values for epsilon and sigma  
betain = betav; epsuv = 0.21146489; siguv = 3.44030500
# Initial values of the soft core potential to be optimized
Ain = 1.0; Bin = 1.0; Cin= 48.0; alphain = 2**(-Cin/6.0); alpha0 = 0

# Initial value of C and alpha for the variance optimization procedure 
DGoptv0 = [Cin, alphain]; var_deltaF_max = 1000 

# Group all input variables together and only make them a calling function
Uinall = lambda optv,epsin,sigin,lamin,rin: [epsin,sigin,Ain,Bin,optv[0],optv[1],alpha0,lamin,rin]
Ginall = lambda optv,epsin,sigin,lamin,rin: [betain,epsin,sigin,Ain,Bin,optv[0],optv[1],alpha0,lamin,rin]
DGinall = lambda optv,rhoin,epsin,sigin,lamin,rin,ilam: [rhoin,betain,epsin,sigin,Ain,Bin,optv[0],optv[1],alpha0,lamin,rin,ilam]

# Define lower and upper limit for lambda, rij and delta lambda for numerical derivative
lam_min = 0; lam_max = 1 - lam_min; dlam = 1E-3 
r_min = 0; r_min_0 = 1E-05; r_max = 10  
rcut = 1E-05; # The cutoff for values of r to prevent numerical unstability at r=0 due to LJ potential 

# Define number of data point in the chosen range of lambda and rij 
Nlam = 41; Nrij = 41; 
Nlam_spline = 21; 

# Define range of lambda and rij to be used for numerical integration 
lambda_range = linspace(lam_min,lam_max,Nlam)
rij_range = linspace(r_min,r_max,Nrij); #rij_range[0] = r_min_0  
lambdaev_spline = linspace(lam_min,lam_max,Nlam_spline) 
xigrid, yigrid = meshgrid(rij_range,lambdaev_spline) # x increment first after y increment 
xfgrid, yfgrid = meshgrid(rij_range,lambda_range) # x increment first after y increment 
umax = 100.0 # Maximum Uij value to use to cap off LJ potential
print 'rij_range = ', rij_range
print 'lambda_range = ', lambda_range

#===================================================================================================
# Functions to calculate all the properties for dilute gas limit
#===================================================================================================
# Define dHdL, ddHdL2, dHdLsq and g(r) for the general softcore case to speed up the 
# calculation rather than obtaining them from symbolic differentiation, which consumes 
# much more CPU time due to repeated calculation
def getdguljfunc(epsinv,siginv,Ainv,Binv,Cinv,alphainv,alpha1inv,laminv,rinv):

    # For any value of reffpC < reffpC_cut, cap them to reffpC_cut to prevent numerical unstability
    # This normally occurs when lambda=1, but we extended to to cases when lambda approaches 1 as well
    reffpC_cut = rcut**Cinv

    # Initialize array for Ulj
    Ulj = len(laminv)*[0]

    rijv = rinv[0]
    # Sub-components required to calculate various function in an efficient manner
    ilam = 0
    for lamv in laminv[0:,0]:
        neg6byC = -6.0/Cinv
        neg12byC = -12.0/Cinv
        epslampAx4 = 4*epsinv*lamv**Ainv
        sigpC = siginv**Cinv
        sigp6 = siginv**6.0
        sigp12 = siginv**12.0
        rpC = rijv**Cinv
        rbysig = rijv/siginv
        rbysigp2 = rbysig*rbysig
        #rbysigpC = rbysig**Cinv
        ralpha = alphainv + alpha1inv*rbysigp2
        rlamp1 = 1 - lamv
        rlampB = rlamp1**Binv
        ralphalampB = ralpha*rlampB*sigpC
        reffpC = ralphalampB + rpC
        # Replace any value of reffpC < reffpC_cut to reffpC_cut to prevent numerical unstability
        if ('%.8f' %(lamv) == '%.8f' %(1.0)):
            ir = 0
            while reffpC[ir] < reffpC_cut:
                reffpC[ir] = reffpC_cut
                ir += 1 
        reff6inv = sigp6*reffpC**neg6byC
        reff12inv = sigp12*reffpC**neg12byC

        Ulj[ilam] = epslampAx4*(reff12inv - reff6inv)
        ilam += 1

    return array(Ulj) 

# Function to calculate dH/dL, ddH/dL2 and (dH/dL)^2 for dilute gas limit
def getdgfunc(rhoinv,betainv,epsinv,siginv,Ainv,Binv,Cinv,alphainv,alpha1inv,laminv,rinv,ilam):

    # For any value of reffpC < reffpC_cut, cap them to reffpC_cut to prevent numerical unstability
    # This normally occurs when lambda=1, but we extended to to cases when lambda approaches 1 as well
    reffpC_cut = rcut**Cinv

    # Sub-components required to calculate various function in an efficient manner
    neg6byC = -6.0/Cinv
    neg6byCm1 = -6.0/Cinv - 1.0
    neg6byCm2 = -6.0/Cinv - 2.0
    neg12byC = -12.0/Cinv
    neg12byCm1 = -12.0/Cinv - 1.0
    neg12byCm2 = -12.0/Cinv - 2.0
    epsAx4 = 4*epsinv*Ainv
    epsAAm1x4 = epsAx4*(Ainv-1)
    epslampAx4 = 4*epsinv*laminv**Ainv
    depslampAx4_dl = epsAx4*laminv**(Ainv-1.0)
    # If A = 1, then the second derivative of lambda**A should be zero, not lambda**(A-2)
    # The if statement below ensure that the above case hold for the 2nd derivative of lambda**A
    if ('%.8f' %(Ainv) == '%.8f' %(1)):
        ddepslampAx4_dl2 = 0
    else: 
        ddepslampAx4_dl2 = epsAAm1x4*laminv**(Ainv-2.0) 
    #ddepslampAx4_dl2 = epsAAm1x4*laminv**(Ainv-2.0) 
    sigpC = siginv**Cinv
    sigp6 = siginv**6.0
    sigp12 = siginv**12.0
    rpC = rinv**Cinv
    rbysig = rinv/siginv
    rbysigp2 = rbysig*rbysig
    ralpha = alphainv + alpha1inv*rbysigp2
    rlamp1 = 1 - laminv
    rlampB = rlamp1**Binv
    ralphalampB = ralpha*rlampB*sigpC
    reffpC = ralphalampB + rpC
    if ('%.8f' %(laminv) == '%.8f' %(1)):
        ir = 0
        while reffpC[ir] < reffpC_cut:
            reffpC[ir] = reffpC_cut
            ir += 1 
    dreffpC_dl = -Binv*ralpha*sigpC*rlamp1**(Binv-1.0)
    # If B = 1, then the second derivative of (1-lambda)**B should also be zero, not (1-lambda)**(B-2)
    # The if statement below ensure that the above case hold for the 2nd derivative of (1-lambda)**B
    if ('%.8f' %(Binv) == '%.8f' %(1)):
        ddreffpC_dl2 = 0 
    else: 
        ddreffpC_dl2 = Binv*(Binv-1.0)*ralpha*sigpC*rlamp1**(Binv-2.0)
    dreffpC_dlp2 = dreffpC_dl*dreffpC_dl 
    reff6inv = sigp6*reffpC**neg6byC 
    reff6invm1 = sigp6*reffpC**neg6byCm1 
    reff6invm2 = sigp6*reffpC**neg6byCm2 
    reff12inv = sigp12*reffpC**neg12byC
    reff12invm1 = sigp12*reffpC**neg12byCm1 
    reff12invm2 = sigp12*reffpC**neg12byCm2 
    dreff6inv_dl = neg6byC*dreffpC_dl*reff6invm1
    ddreff6inv_dl2 = neg6byC*neg6byCm1*dreffpC_dlp2*reff6invm2 + neg6byC*ddreffpC_dl2*reff6invm1 
    dreff12inv_dl = neg12byC*dreffpC_dl*reff12invm1
    ddreff12inv_dl2 = neg12byC*neg12byCm1*dreffpC_dlp2*reff12invm2 + neg12byC*ddreffpC_dl2*reff12invm1 

    # All functions are define here
    Ulj = epslampAx4*(reff12inv - reff6inv)
    dUdL_lam = epslampAx4*(dreff12inv_dl - dreff6inv_dl) + depslampAx4_dl*(reff12inv - reff6inv)
    ddUdL2_lam = epslampAx4*(ddreff12inv_dl2 - ddreff6inv_dl2) + depslampAx4_dl*(dreff12inv_dl - dreff6inv_dl) + depslampAx4_dl*(dreff12inv_dl - dreff6inv_dl) + ddepslampAx4_dl2*(reff12inv - reff6inv)
    dUdLsq_lam = dUdL_lam**2.0

    nsam = len(rinv)
    gr = exp(-betainv*Ulj)

    # dHdL, ddHdL2 and dHdL^2 as a function of rij at a particular value of lambda for dilute gas limit
    dHdL_lam = dUdL_lam*gr*rhoinv*4*pi*rinv**2.0
    dHdLsq_lam = dUdLsq_lam*gr*rhoinv*4*pi*rinv**2.0
    ddHdL2_lam = ddUdL2_lam*gr*rhoinv*4*pi*rinv**2.0

    # For pure hardcore potential at lambda = 1, g(r=0) = 0. Therefore any properties that are the product of
    # g(r) would have a zero value at r=0. The if statement below prevent numerical overflow for such a case
    if ('%.8f' %(laminv) == '%.8f' %(1)):
        dHdL_lam[0:3] = 0
        dHdLsq_lam[0:3] = 0
        ddHdL2_lam[0:3] = 0

    # Determine the avergare value for dHdL, ddHdL2 and dHdL^2
    dHdL = simps(dHdL_lam,rinv)
    dHdLsq = simps(dHdLsq_lam,rinv)
    ddHdL2 = simps(ddHdL2_lam,rinv)

    return dUdL_lam, dUdLsq_lam, dHdL, dHdLsq, ddHdL2, gr 

# Function to calculate the variance via numerical integration in rij and lambda using spline for dilute gas limit
def getdg1dsplinelinear(zoptv,escale,uinvcutoff,dudlcutoff,uljinverse,laminv,rinv):

    # For dilute gasi limit, only dH/dL and ddHdL2 for O are obtained as there is no correlation between O and H
    rho = rho_OW
    nl = len(laminv)
    nr = len(rinv)
    uljinvb = r_[uljinverse[0],zoptv] # Add Ulj at lambda = 0 in front of the array
    uljinv = r_[uljinvb,uljinverse[-1]] # Add Ulj at lambda = 1 at the end of the array

    # Initialize array for the indice of all r value
    xlj = nr*[0] 
    dxdl = nr*[0] 

    # Use 1D-spline to spline in the lambda direction for all r values 
    for ir in range(nr):
        index = index0 + ir
        tck = interp.splrep(lambdaev_spline,uljinv[index],k=1,task=0)
        # Evaluate derivatives of Ulj inverse for 1D-spline with the more crowded lambda grid
        xlj[ir] = interp.splev(laminv,tck,der=0) 
        dxdl[ir] = interp.splev(laminv,tck,der=1)
        ir += 1

    xlj = array(xlj)
    xlj = xlj.transpose()
    dxdl = array(dxdl)
    dxdl = dxdl.transpose()
    ulj = 1.0/xlj - escale - 1

    if method == 'spline':
        # Determine total variance using 1D spline in lambda ditrection 
        xlj_min = xlj[0:,0:].min()
        if xlj_min < 0:
            var_deltaF = abs(log(-xlj_min))*var_deltaF_max
        else:
            dudl = -dxdl/(xlj**2.0)
            #dudl = exp(xlj)*dxdl
            #dudl = -dxdl/(betain*xlj)
            dudlsq = dudl**2.0

            ilam = 0
            for lamv in laminv:
                ## Truncating value of the derivatives if it is too big to prevent overfloating
                #if dudl[ilam].max() > dudlcutoff:
                #    ir = 0
                #    while dudl[ilam][ir] > dudlcutoff:
                #        dudl[ilam][ir] = dudlcutoff
                #        dudlsq[ilam][ir] = dudl[ilam][ir]**2 
                #        ir += 1
                gr = exp(-betain*ulj[ilam]) 
                grw = gr*4*pi*rho*rij_range**2.0
                dHdL[ilam] = simps(dudl[ilam]*grw,rij_range)
                dHdLsq[ilam] = simps(dudlsq[ilam]*grw,rij_range)
    
                ilam += 1
                
            # Variance of the free energy using spline
            var_deltaF = simps(dHdLsq,laminv)

    else:
        # Determine total variance using 1D spline in lambda ditrection 
        xlj_min = xlj[0:,0:].min()
        if xlj_min < 0:
            var_deltaF = abs(log(-xlj_min))*var_deltaF_max
        else:
            # Determine total variance using finite difference
            # Constant required for forward, central and backward finite
            # difference derivatives
            dlamv = laminv[1] - laminv[0]
            dlamvx2h = 2*dlamv
            dlamvx12h = 12*dlamv
            dlamvxh2 = dlamv**2.0
            dlamvx12h2 = 12*dlamvxh2
            FDIL = 2 # Limit of index below which forward difference is used
            BDIL = len(laminv) - 3 # Limit of index above which backward difference is used

            # Precalculate value required for finite difference to speed up the
            # algorithsm using arraywise
            xljx1byh = xlj/dlamv
            xljx2byh = 2*xljx1byh
            xljx1by2h = xlj/dlamvx2h
            xljx3by2h = 3*xljx1by2h
            xljx4by2h = 4*xljx1by2h
            xljx1by12h = xlj/dlamvx12h
            xljx8by12h = 8*xljx1by12h
            xljx1byh2 = xlj/dlamvxh2
            xljx2byh2 = 2*xljx1byh2
            xljx4byh2 = 4*xljx1byh2
            xljx5byh2 = 5*xljx1byh2
            xljx1by12h2 = xlj/dlamvx12h2
            xljx16byh2 = 16*xljx1by12h2
            xljx30byh2 = 30*xljx1by12h2

            for ilam in range(nl):
                ilamm3 = ilam - 3
                ilamm2 = ilam - 2
                ilamm1 = ilam - 1
                ilamp1 = ilam + 1
                ilamp2 = ilam + 2
                ilamp3 = ilam + 3
                if ilam == 0:
                    ## First order correction
                    #dUdL = -uljx1byh[ilam]+uljx1byh[ilamp1]
                    #ddUdL2 = uljx1byh2[ilam]-uljx2byh2[ilamp1]+uljx1byh2[ilamp2]
                    # Second order correction
                    duljmap = -xljx3by2h[ilam]+xljx4by2h[ilamp1]-xljx1by2h[ilamp2]
                    #ddUdL2 = uljx2byh2[ilam]-uljx5byh2[ilamp1]+uljx4byh2[ilamp2]-uljx1byh2[ilamp3]
                elif ilam == 1:
                    # Second order correction
                    #dUdL = -uljx1by2h[ilamm1]+uljx1by2h[ilamp1]
                    #ddUdL2 = uljx1byh2[ilamm1]-uljx2byh2[ilam]+uljx1byh2[ilamp1]
                    ## Fourth order correction
                    duljmap = xljx1by12h[ilamm1]-xljx8by12h[ilamm1]+xljx8by12h[ilamp1]-xljx1by12h[ilamp2]
                    #ddUdL2 = -uljx1by12h2[ilamm1]+uljx16byh2[ilamm1]-uljx30byh2[ilam]+uljx16byh2[ilamp1]-uljx1by12h2[ilamp2]
                elif ilam >= BDIL:
                    ## First order correction
                    #dUdL = -uljx1byh[ilamm1]+uljx1byh[ilam]
                    #ddUdL2 = uljx1byh2[ilamm2]-uljx2byh2[ilamm1]+uljx1byh2[ilam]
                    # Second order correction
                    duljmap = xljx1by2h[ilamm2]-xljx4by2h[ilamm1]+xljx3by2h[ilam]
                    #ddUdL2 = -uljx1byh2[ilamm3]+uljx4byh2[ilamm2]-uljx5byh2[ilamm1]+uljx2byh2[ilam]
                else:
                    # Second order correction
                    #dzdl = -uljx1by2h[ilamm1]+uljx1by2h[ilamp1]
                    #ddUdL2 = uljx1byh2[ilamm1]-uljx2byh2[ilam]+uljx1byh2[ilamp1]
                    ## Fourth order correction
                    duljmap = xljx1by12h[ilamm2]-xljx8by12h[ilamm1]+xljx8by12h[ilamp1]-xljx1by12h[ilamp2]
                    #ddUdL2 = -uljx1by12h2[ilamm2]+uljx16byh2[ilamm1]-uljx30byh2[ilam]+uljx16byh2[ilamp1]-uljx1by12h2[ilamp2]
                dudl = -duljmap/(xlj[ilam]**2.0)
                #dudl = exp(xlj[ilam])*duljmap
                #dudl = -duljmap/(betain*xlj[ilam])
                dudlsq = dudl**2.0
                gr = exp(-betain*ulj[ilam]) #*4*pi*rho*rinv**2.0
                grw = gr*4*pi*rho*rinv**2.0
                dHdL[ilam] = simps(dudl*grw,rinv)
                #ddHdL2[ilam] = simps(ddUdL2*grw,rinv)
                dHdLsq[ilam] = simps(dudlsq*grw,rinv)

            # Variance of the free energy using spline
            var_deltaF = simps(dHdLsq,laminv)
            #var_deltaF = (simps(ddHdL2,lam_range) - dHdL[-1] + dHdL[0])/ betain

    print 'var(deltaF) = ', var_deltaF

    return var_deltaF

# Function to calculate the variance via numerical integration in rij and lambda using spline for dilute gas limit
def getdg1dsplinecubic(zoptv,escale,uinvcutoff,dudlcutoff,uljinverse,laminv,rinv):

    # For dilute gasi limit, only dH/dL and ddHdL2 for O are obtained as there is no correlation between O and H
    rho = rho_OW
    nl = len(laminv)
    nr = len(rinv)
    uljinvb = r_[uljinverse[0],zoptv] # Add Ulj at lambda = 0 in front of the array
    uljinv = r_[uljinvb,uljinverse[-1]] # Add Ulj at lambda = 1 at the end of the array

    # Initialize array for the indice of all r value
    xlj = nr*[0] 
    dxdl = nr*[0] 

    # Use 1D-spline to spline in the lambda direction for all r values 
    for ir in range(nr):
        index = index0 + ir
        tck = interp.splrep(lambdaev_spline,uljinv[index],k=3,task=0)
        # Evaluate derivatives of Ulj inverse for 1D-spline with the more crowded lambda grid
        xlj[ir] = interp.splev(laminv,tck,der=0) 
        dxdl[ir] = interp.splev(laminv,tck,der=1)
        ir += 1

    xlj = array(xlj)
    xlj = xlj.transpose()
    dxdl = array(dxdl)
    dxdl = dxdl.transpose()
    ulj = 1.0/xlj - escale - 1

    if method == 'spline':
        # Determine total variance using 1D spline in lambda ditrection 
        xlj_min = xlj[0:,0:].min()
        if xlj_min < 0:
            var_deltaF = abs(log(-xlj_min))*var_deltaF_max
        else:
            dudl = -dxdl/(xlj**2.0)
            #dudl = exp(xlj)*dxdl
            #dudl = -dxdl/(betain*xlj)
            dudlsq = dudl**2.0

            ilam = 0
            for lamv in laminv:
                ## Truncating value of the derivatives if it is too big to prevent overfloating
                #if dudl[ilam].max() > dudlcutoff:
                #    ir = 0
                #    while dudl[ilam][ir] > dudlcutoff:
                #        dudl[ilam][ir] = dudlcutoff
                #        dudlsq[ilam][ir] = dudl[ilam][ir]**2 
                #        ir += 1
                gr = exp(-betain*ulj[ilam]) 
                grw = gr*4*pi*rho*rij_range**2.0
                dHdL[ilam] = simps(dudl[ilam]*grw,rij_range)
                dHdLsq[ilam] = simps(dudlsq[ilam]*grw,rij_range)
    
                ilam += 1
                
            # Variance of the free energy using spline
            var_deltaF = simps(dHdLsq,laminv)

    else:
        # Determine total variance using 1D spline in lambda ditrection 
        xlj_min = xlj[0:,0:].min()
        if xlj_min < 0:
            var_deltaF = abs(log(-xlj_min))*var_deltaF_max
        else:
            # Determine total variance using finite difference
            # Constant required for forward, central and backward finite
            # difference derivatives
            dlamv = laminv[1] - laminv[0]
            dlamvx2h = 2*dlamv
            dlamvx12h = 12*dlamv
            dlamvxh2 = dlamv**2.0
            dlamvx12h2 = 12*dlamvxh2
            FDIL = 2 # Limit of index below which forward difference is used
            BDIL = len(laminv) - 3 # Limit of index above which backward difference is used

            # Precalculate value required for finite difference to speed up the
            # algorithsm using arraywise
            xljx1byh = xlj/dlamv
            xljx2byh = 2*xljx1byh
            xljx1by2h = xlj/dlamvx2h
            xljx3by2h = 3*xljx1by2h
            xljx4by2h = 4*xljx1by2h
            xljx1by12h = xlj/dlamvx12h
            xljx8by12h = 8*xljx1by12h
            xljx1byh2 = xlj/dlamvxh2
            xljx2byh2 = 2*xljx1byh2
            xljx4byh2 = 4*xljx1byh2
            xljx5byh2 = 5*xljx1byh2
            xljx1by12h2 = xlj/dlamvx12h2
            xljx16byh2 = 16*xljx1by12h2
            xljx30byh2 = 30*xljx1by12h2

            for ilam in range(nl):
                ilamm3 = ilam - 3
                ilamm2 = ilam - 2
                ilamm1 = ilam - 1
                ilamp1 = ilam + 1
                ilamp2 = ilam + 2
                ilamp3 = ilam + 3
                if ilam == 0:
                    ## First order correction
                    #dUdL = -uljx1byh[ilam]+uljx1byh[ilamp1]
                    #ddUdL2 = uljx1byh2[ilam]-uljx2byh2[ilamp1]+uljx1byh2[ilamp2]
                    # Second order correction
                    duljmap = -xljx3by2h[ilam]+xljx4by2h[ilamp1]-xljx1by2h[ilamp2]
                    #ddUdL2 = uljx2byh2[ilam]-uljx5byh2[ilamp1]+uljx4byh2[ilamp2]-uljx1byh2[ilamp3]
                elif ilam == 1:
                    # Second order correction
                    #dUdL = -uljx1by2h[ilamm1]+uljx1by2h[ilamp1]
                    #ddUdL2 = uljx1byh2[ilamm1]-uljx2byh2[ilam]+uljx1byh2[ilamp1]
                    ## Fourth order correction
                    duljmap = xljx1by12h[ilamm1]-xljx8by12h[ilamm1]+xljx8by12h[ilamp1]-xljx1by12h[ilamp2]
                    #ddUdL2 = -uljx1by12h2[ilamm1]+uljx16byh2[ilamm1]-uljx30byh2[ilam]+uljx16byh2[ilamp1]-uljx1by12h2[ilamp2]
                elif ilam >= BDIL:
                    ## First order correction
                    #dUdL = -uljx1byh[ilamm1]+uljx1byh[ilam]
                    #ddUdL2 = uljx1byh2[ilamm2]-uljx2byh2[ilamm1]+uljx1byh2[ilam]
                    # Second order correction
                    duljmap = xljx1by2h[ilamm2]-xljx4by2h[ilamm1]+xljx3by2h[ilam]
                    #ddUdL2 = -uljx1byh2[ilamm3]+uljx4byh2[ilamm2]-uljx5byh2[ilamm1]+uljx2byh2[ilam]
                else:
                    # Second order correction
                    #dzdl = -uljx1by2h[ilamm1]+uljx1by2h[ilamp1]
                    #ddUdL2 = uljx1byh2[ilamm1]-uljx2byh2[ilam]+uljx1byh2[ilamp1]
                    ## Fourth order correction
                    duljmap = xljx1by12h[ilamm2]-xljx8by12h[ilamm1]+xljx8by12h[ilamp1]-xljx1by12h[ilamp2]
                    #ddUdL2 = -uljx1by12h2[ilamm2]+uljx16byh2[ilamm1]-uljx30byh2[ilam]+uljx16byh2[ilamp1]-uljx1by12h2[ilamp2]
                dudl = -duljmap/(xlj[ilam]**2.0)
                #dudl = exp(xlj[ilam])*duljmap
                #dudl = -duljmap/(betain*xlj[ilam])
                dudlsq = dudl**2.0
                gr = exp(-betain*ulj[ilam]) #*4*pi*rho*rinv**2.0
                grw = gr*4*pi*rho*rinv**2.0
                dHdL[ilam] = simps(dudl*grw,rinv)
                #ddHdL2[ilam] = simps(ddUdL2*grw,rinv)
                dHdLsq[ilam] = simps(dudlsq*grw,rinv)

            # Variance of the free energy using spline
            var_deltaF = simps(dHdLsq,laminv)
            #var_deltaF = (simps(ddHdL2,lam_range) - dHdL[-1] + dHdL[0])/ betain

    print 'var(deltaF) = ', var_deltaF

    return var_deltaF

# Function to calculate the spline data point 
def getdgdataspline(uljinv,escale,laminv,rinv):

    # Initialize array for the indice of all r value
    nr = len(rinv)
    zlj = nr*[0]

    # Use 1D-spline to spline in the lambda direction for all r values 
    for ir in range(nr):
        index = index0 + ir
        tck = interp.splrep(lambdaev_spline,uljinv[index],k=3,task=0)
        # Evaluate derivatives of Ulj inverse for 1D-spline with the
        # more crowded lambda grid
        zlj[ir] = interp.splev(laminv,tck,der=0)

        ir += 1

    zlj = array(zlj)
    #ulj = exp(zlj) - escale - 1
    #ulj = -log(zlj)/betain
    ulj = 1.0/zlj-escale-1

    return ulj

# Function to evaluate all properties for the optimal pathway for dilute gas limit 
# This function evaluates the derivatives term ddLdHdL much faster than the function getdgvisual 
def getdgoptimal(optv,function,epsuv,siguv,lam_range,rij_range):

    ilam = 0
    for lamv in lam_range:
        # For dilute gasi limit, only dH/dL for O are obtained as there is no correlation between O and H 
        rho = rho_OW
        if function == 'symbolic':
            #Input for the lambdify function to calculate dHdL, ddHdL2 and dHdLsq via general symbolic function
            dginput = (optv,rho,epsuv,siguv,lamv,rij_range)
            OdHdL_dg[ilam], OdHdLsq_dg[ilam], OddHdL2_dg[ilam] = getdgdhdl(*dginput)
        elif function == 'define':
            #Input for the predefined function to calculate dHdL, ddHdL2, dHdLsq and g(r) via coding
            #DGin = DGinall(optv,rho,epsuv,siguv,1.0,rij_range,ilam)
            DGin = DGinall(optv,rho,epsuv,siguv,lamv,rij_range,ilam)
            OdHdL_dg_lam, OdHdLsq_dg_lam, OdHdL_dg[ilam], OdHdLsq_dg[ilam], OddHdL2_dg[ilam], gr_dg = getdgfunc(*DGin)        
 
        ilam += 1	

    # Constant required for forward, central and backward finite difference derivatives
    dlamv = lam_range[1] - lam_range[0]
    dlamvx2 = 2*dlamv
    dlamvx12 = 12*dlamv
    FDIL = 2 # Limit of index below which forward difference is used
    BDIL = len(lam_range) - 3 # Limit of index above which backward difference is used 
    # Evaluating ddLdHdL using all three finite difference derivatives methods with constant dlamv
    for ilam in range(len(lam_range)):
        ilamm3 = ilam - 3
        ilamm2 = ilam - 2
        ilamm1 = ilam - 1
        ilamp1 = ilam + 1
        ilamp2 = ilam + 2
        ilamp3 = ilam + 3
        if ilam <= FDIL:
            OddLdHdL_dg[ilam] = (-OdHdL_dg[ilamp2]+4*OdHdL_dg[ilamp1]-3*OdHdL_dg[ilam])/dlamvx2    
        elif ilam >= BDIL:
            OddLdHdL_dg[ilam] = (OdHdL_dg[ilamm2]-4*OdHdL_dg[ilamm1]+3*OdHdL_dg[ilam])/dlamvx2    
        else:
            OddLdHdL_dg[ilam] = (OdHdL_dg[ilamm2]-8*OdHdL_dg[ilamm1]+8*OdHdL_dg[ilamp1]-OdHdL_dg[ilamp2])/dlamvx12    

    if function == 'symbolic':
        # Use symbolic function to evaluate dHdL_lam and dHdLsq_lam
        Uin = Uinall(optv,epsuv,siguv,lam_range[-1],rij_range)
        OdHdL_dg_lam = dUdL_fun(*Uin)*gr_dg*rho*4*pi*rij_range**2.0
        OdHdLsq_dg_lam = dUdLsq_fun(*Uin)*gr_dg*rho*4*pi*rij_range**2.0

    # Combine values for O and H to get the actual total value 
    dHdL = OdHdL_dg #+ HdHdL
    dHdLsq = OdHdLsq_dg #+ HdHdLsq
    ddHdL2 = OddHdL2_dg #+ HddHdL2
    ddLdHdL = OddLdHdL_dg #+ HddLdHdL
    
    # Variance of the dHdL 
    var_dHdL = dHdLsq - dHdL**2
    var_dHdL_new = (ddHdL2 - ddLdHdL)/betain
    
    # Variance of the free energy  
    #var_deltaF = trapz(var_dHdL,lam_range) 
    #var_deltaF_new = (trapz(ddHdL2,lam_range) - dHdL[-1] + dHdL[0])/ betain
    #var_deltaF_alt = (trapz(ddHdL2-ddLdHdL,lam_range))/ betain
    var_deltaF = simps(var_dHdL,lam_range) 
    var_deltaF_new = (simps(ddHdL2,lam_range) - dHdL[-1] + dHdL[0])/ betain
    #var_deltaF_new = simps(dHdLsq,lam_range)
    var_deltaF_alt = (simps(ddHdL2-ddLdHdL,lam_range))/ betain
    print 'var(deltaF) = ' + repr(var_deltaF)
    print 'var(deltaF)_alt = ' + repr(var_deltaF_alt)
    print 'C, alpha = ', optv, ';   var(deltaF)_new = ', var_deltaF_new

    return OdHdL_dg_lam, OdHdLsq_dg_lam, dHdL, ddHdL2, ddLdHdL, dHdLsq, var_dHdL, var_dHdL_new, var_deltaF, var_deltaF_new 

#===================================================================================================
# Performing optimization via cubic spline 
#===================================================================================================
# Define the array for dHdl and dHdL^2 
dHdL_dg = numpy.zeros(len(lambda_range),float)
dHdLsq_dg = numpy.zeros(len(dHdL_dg),float)
ddHdL2_dg = numpy.zeros(len(dHdL_dg),float)
ddLdHdL_dg = numpy.zeros(len(dHdL_dg),float)
OdHdL_dg = numpy.zeros(len(lambda_range),float)
OdHdLsq_dg = numpy.zeros(len(dHdL_dg),float)
OddHdL2_dg = numpy.zeros(len(dHdL_dg),float)
OddLdHdL_dg = numpy.zeros(len(dHdL_dg),float)
dHdL = numpy.zeros(len(lambda_range),float)
dHdLsq = numpy.zeros(len(dHdL),float)
ddHdL2 = numpy.zeros(len(dHdL),float)

bounds = ([1, 200],[1E-5,100])
uscale = epsuv # scaling valing to be added the potential for inverse mapping  
uinvcut = 1E-08 # Cutoff for the inverse mapping of the LJ potential 
dudlcut = 1E+50 #Cutoff for the derivatives for the inverse mapping of the LJ potential 
DGin_full = (epsuv,siguv,Ain,Bin,Cin,alphain,alpha0,yigrid,xigrid)
uljfull = getdguljfunc(*DGin_full)        
uljmap = 1.0/(uljfull+uscale+1) # inverse mapping of Ulj for spline
uljmap0 = uljmap[1:-1].flatten(0) # Turn array into Rank-1 array
uljmapall = uljmap.flatten(0) # 1D-array of inverse mapping of Ulj for spline
rigrid = xigrid.flatten(0) # Change rgrid into a 1D or rank 1 array
ligrid = yigrid.flatten(0) # Change lgrid into a 1D or rank 1 array
rfgrid = xfgrid.flatten(0) # Change rgrid into a 1D or rank 1 array
lfgrid = yfgrid.flatten(0) # Change lgrid into a 1D or rank 1 array
index0 = arange(0,len(rigrid),len(uljmap[0])) # index where all points have the same lambda value

#===================================================================================================
# Evaluate and write the initial potential to be optimized to a file
#===================================================================================================

# Spline the intial Ulj data to obtain more intermediate data points
uinitlist = getdgdataspline(uljmapall,uscale,lambda_range,rij_range)
uinit = uinitlist.flatten(1)

# Evaluate the initial starting potential
DGopt = (DGoptv0,functype,epsuv,siguv,lambda_range,rij_range)
dUdL_lam_dg, dUdLsq_lam_dg, dHdL_dg, ddHdL2_dg, ddLdHdL_dg, dHdLsq_dg, var_dHdL_dg, var_dHdL_new_dg, var_deltaF_dg, var_deltaF_new_dg = getdgoptimal(*DGopt)

# Write value of the starting potential to be optimized using finite difference with dilute gas approximation 
dhdl_dg_initial='initial-DG_spline-cubic_finite-CD_map-A%i-B%i-C%i-alpha%.4f-%i_lambda-%i_rij.dg'%(Ain,Bin,Cin,alphain,Nlam_spline,Nrij) 
datafile = open(dhdl_dg_initial, 'w')
datafile.write("%16s %16s %16s %16s %16s %16s %16s %16s %16s %16s\n" %('Nlam','Nrij','epsilon','sigma','Ain','Bin','Cin', 'alphain','alpha1in', 'var(deltaF)'))
datafile.write("%16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n"  %(Nlam,Nrij,epsuv,siguv,Ain,Bin,Cin,alphain,alpha0,var_deltaF_new_dg))
datafile.write("%16s %16s %16s\n" %('lambda', 'rij', 'Ulj'))
for k in range(len(rfgrid)):
    datafile.write("%16.8e %16.8e %16.8e\n"  %(lfgrid[k], rfgrid[k], uinit[k]))
datafile.close()

#===================================================================================================
# Perform the 1D-spline optimization here 
#===================================================================================================
# Set array of bounds for value of the optimized function
uljmap_min= uljmapall.min()
uljmap_max = uljmapall.max()
bound = [uljmapall.min(), uljmapall.max()]
bounds = []
for ilam in range(len(uljmap0)):
    bounds.append(bound)
# Convert list of bounds into tuple to be used for the optimization
bounds = tuple(bounds)

print 'start DG optimize', time.ctime()
# First perform the optimization using linear interpolation scheme as it
# is very fast
zopt, var_deltaF, DG_dict = fmin_l_bfgs_b(getdg1dsplinelinear, uljmap0, args=(uscale,uinvcut,dudlcut,uljmap,lambda_range,rij_range), approx_grad=1, bounds=bounds, m=10, factr=1e03, pgtol=1e-04, epsilon=1e-04, maxfun=1e+04)
# Then refine the results using cubic interpolation scheme
print 'refine 1D-spline from linear to cubic' 
uinvopt, var_deltaF, DG_dict = fmin_l_bfgs_b(getdg1dsplinecubic, zopt, args=(uscale,uinvcut,dudlcut,uljmap,lambda_range,rij_range), approx_grad=1, bounds=bounds, m=10, factr=1e03, pgtol=1e-05, epsilon=1e-05, maxfun=1e+04)
uinvoptv = r_[uljmap[0],uinvopt]
uinvoptv = r_[uinvoptv,uljmap[-1]]
print 'optimized var(deltaF) = ', var_deltaF
print 'finish DG optimize', time.ctime()

#===================================================================================================
# Done with 1D-spline optimization then write results to file
#===================================================================================================
# Spline the intial Ulj data to obtain more intermediate data points
uoptvlist = getdgdataspline(uinvoptv,uscale,lambda_range,rij_range)
uoptv = uoptvlist.flatten(1)

# Data for optimal B-spline using dilute gas approximation 
dhdl_dg_optimal='optimal-DG_spline-cubic_finite-CD_map-A%i-B%i-C%i-alpha%.4f-%i_lambda-%i_rij.dg'%(Ain,Bin,Cin,alphain,Nlam_spline,Nrij) 
datafile = open(dhdl_dg_optimal, 'w')
datafile.write("%16s %16s %16s %16s %16s %16s %16s %16s %16s %16s\n" %('Nlam','Nrij','epsilon','sigma','Ain','Bin','Cin', 'alphain','alpha1in', 'var(deltaF)'))
datafile.write("%16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n"  %(Nlam,Nrij,epsuv,siguv,Ain,Bin,Cin,alphain,alpha0,var_deltaF))
datafile.write("%16s %16s %16s\n" %('lambda', 'rij', 'Ulj'))
for k in range(len(rfgrid)):
    datafile.write("%16.8e %16.8e %16.8e\n"  %(lfgrid[k], rfgrid[k], uoptv[k]))
datafile.close()