# Commits

committed ffcf4e8

implemented substrate potentials fitting and extended the isotherm layer program to support de oliveira substrate

• Participants
• Parent commits 79fe120

# File docs/sm_final_report/sm_final_report.pdf

Binary file modified.

# File docs/sm_final_report/sm_final_report.tex

 \tableofcontents


+\begin{enumerate}
+  \item problema
+  \item teoria de oliveira
+  \item teoria ebner
+  \item approssimazioni
+  \item algoritmo e implementazione (thread pool)
+  \item risultati (plot vari)
+  \item confronto ebner - de oliveira (fitting potenziale)
+  \item conclusioni (C=2,D=2.2 artificiali per fare il salto)
+\end{enumerate}
+
+
 \section{Introduction}
 This is the first section.





-Ebner::Ebner(int N)
+Ebner::Ebner(int N, double C, double D)
 : _N(N)
 , _sigma(1.0)
 , _sigma_w(1.0)
 , _D(1.0)
 {
     _Tstar = 1.0;
-    _pot = new potentials::DeOliveiraPotentials( _N, 4.20, 7.85, 1.0, 1.0 );
+    _pot = new potentials::DeOliveiraPotentials( _N, C, D, 1.0, 1.0 );

     // 'f' as defined in ebner (12)
     _FD = _pot->f_pot(0,0);

 public:

     /// create an instance using de oliveira potentials
-    Ebner(int N);
+    Ebner(int N, double C=4.20, double D=7.85);
     /// create an instance using long distance potentials
     Ebner(int N, double sigma, double sigma_w, double alpha, double eps_k, double eps_w_k, double D, bool long_dist=true);
     virtual ~Ebner();

     dsfmt_init_gen_rand(&dsfmt,20110912);

     cmdline::parser params;
-    params.add<double>("alpha", 'a', "substrate strength", false, 1.75, cmdline::range(0.0,100.0));
+    params.add<double>("alpha", 'a', "substrate strength (for Ebner substrate formulation)", false, 1.75, cmdline::range(0.0,100.0));
     params.add<double>("tstar", 't', "normalized temperature", false, 2.40, cmdline::range(0.0,10.0));
     params.add<int>("samples", 's', "number of samples", false, 100, cmdline::range(1, 10000));
     params.add<std::string>("interaction_range", 'r', "interparticle interaction range", false, "long", cmdline::oneof<std::string>("short", "long"));
-    params.add("gzip", '\0', "gzip when transfer");
+    params.add<std::string>("substrate_potential", 'p', "substrate potential formulation", false, "ebner", cmdline::oneof<std::string>("ebner", "oliveira"));
+    params.add<double>("substr_c", 'c', "C parameter (for De Oliveira substrate formulation)", false, 4.20, cmdline::range(0.0,100.0));
+    params.add<double>("substr_d", 'd', "D parameter (for De Oliveira substrate formulation)", false, 7.85, cmdline::range(0.0,100.0));
     params.parse_check(argc, argv);

     // user parameters
     double alpha = params.get<double>("alpha");
+    double doC = params.get<double>("substr_c");
+    double doD = params.get<double>("substr_d");
     double Tstar = params.get<double>("tstar");
     int num_samples = params.get<int>("samples");
+    std::string substrate_potential = params.get<std::string>("substrate_potential");
     std::string potential_range = params.get<std::string>("interaction_range");
     bool long_range = potential_range=="long";


     const double sigma = 3.405;       // Angstrom
     const double sigma_w = 3.735;     // Angstrom
-    const double eps_k = 120.0;         // K
+    const double eps_k = 120.0;       // K
     const double eps_w_k = 160.0;     // K
-    const double D = 0.846;       // n_w * sigma_w**3
+    const double D = 0.846;           // n_w * sigma_w**3
     const int N = 20;


-    Ebner ebner(N, sigma, sigma_w, alpha, eps_k, eps_w_k, D, long_range);
-    ebner.set_Tstar(Tstar);
+    Ebner *ebner;
+    if (substrate_potential=="ebner")
+        ebner = new Ebner(N, sigma, sigma_w, alpha, eps_k, eps_w_k, D, long_range); // use substrate pot as formulated by ebner
+    else
+        ebner = new Ebner(N, doC, doD); // use substrate pot as formulated by de oliveira
+    ebner->set_Tstar(Tstar);

     const double gamma_min = -0.5;
     const double gamma_max = -0.01;
     std::vector<SampleData> samples(num_samples);
     for (unsigned int i=0; i<num_samples; i++)
     {
-        samples[i].ebner = &ebner;
+        samples[i].ebner = ebner;
         samples[i].x = x_min + (x_max-x_min)*(num_samples-1-i)/(num_samples-1);
         samples[i].gamma = log(samples[i].x)/scaling;
         samples[i].final_energy = -1.0;
     // dump events to file
     events_globals::dump_to_file("ebner_isotherm_profile.json");

+    delete ebner;
+
     // gracefully shutdown systems
     events_globals::shutdown();
     memory_globals::shutdown();



 def plot_data(data):
+    plt.subplot(1,3,1)
     plt.hold(1)
-    plt.loglog(data.x,data.y1, '-r')
-    plt.loglog(data.x,data.y2, '-b')
+    plt.plot(data.x,data.y1, '-rx')
+    plt.plot(data.x,data.y2, '-bx')
+    plt.axis([0,10,0,8])
+    plt.title('normal scale')
+    plt.subplot(1,3,2)
+    plt.hold(1)
+    plt.semilogy(data.x,data.y1, '-rx')
+    plt.semilogy(data.x,data.y2, '-bx')
+    # plt.axis([0,10,0,8])
+    plt.title('log y scale')
+    plt.subplot(1,3,3)
+    plt.hold(1)
+    plt.loglog(data.x,data.y1, '-rx')
+    plt.loglog(data.x,data.y2, '-bx')
+    plt.title('log-log scale')
     plt.show()
-


 if __name__ == '__main__':

+# plot the ising magnetization
+
+import os, sys
+import matplotlib.pyplot as plt
+
+
+class Data(object):
+    """docstring for Data"""
+    def __init__(self):
+        super(Data, self).__init__()
+        self.x = []
+        self.y1 = []
+        self.y2 = []
+    def append(self,x,y1,y2):
+        self.x.append(x)
+        self.y1.append(y1)
+        self.y2.append(y2)
+
+
+def read_file(fname):
+    f = open(fname,"r")
+    l0 = Data()
+    for line in f.readlines():
+        _x,_y1,_y2 = line.strip().split()[0:3]
+        x = float(_x)
+        y1 = float(_y1)
+        y2 = float(_y2)
+        l0.append(x,y1,y2)
+    return l0
+
+
+def plot_data(data):
+    plt.hold(1)
+    plt.plot(data.x,data.y1, '-r')
+    plt.plot(data.x,data.y2, '-b')
+    plt.axis([0,20,-10,8])
+    plt.title('normal scale')
+    plt.show()
+
+
+if __name__ == '__main__':
+    if len(sys.argv)<2:
+        print("missing filename. exit.\n")
+        sys.exit()
+
+    fname = sys.argv[1]
+    data = read_file(fname)
+    plot_data(data)

+1.000000 -5.976648 -7.850000
+2.000000 -1.160027 -0.525000
+3.000000 -0.350192 -0.155556
+4.000000 -0.148805 -0.065625
+5.000000 -0.076505 -0.033600
+6.000000 -0.044395 -0.019444
+7.000000 -0.028012 -0.012245
+8.000000 -0.018793 -0.008203
+9.000000 -0.013214 -0.005761
+10.000000 -0.009642 -0.004200
+11.000000 -0.007249 -0.003156
+12.000000 -0.005587 -0.002431
+13.000000 -0.004397 -0.001912
+14.000000 -0.003522 -0.001531
+15.000000 -0.002865 -0.001244
+16.000000 -0.002361 -0.001025
+17.000000 -0.001969 -0.000855
+18.000000 -0.001659 -0.000720
+19.000000 -0.001411 -0.000612
+20.000000 -0.001210 -0.000525
+21.000000 -0.001046 -0.000454
+22.000000 -0.000910 -0.000394
+23.000000 -0.000796 -0.000345
+24.000000 -0.000701 -0.000304
+25.000000 -0.000620 -0.000269
+26.000000 -0.000551 -0.000239
+27.000000 -0.000492 -0.000213
+28.000000 -0.000442 -0.000191
+29.000000 -0.000397 -0.000172

+1.000000 -4.357973 -4.358000
+2.000000 -0.845853 -0.845875
+3.000000 -0.255348 -0.250630
+4.000000 -0.108504 -0.105734
+5.000000 -0.055785 -0.054136
+6.000000 -0.032371 -0.031329
+7.000000 -0.020425 -0.019729
+8.000000 -0.013703 -0.013217
+9.000000 -0.009635 -0.009283
+10.000000 -0.007031 -0.006767
+11.000000 -0.005286 -0.005084
+12.000000 -0.004074 -0.003916
+13.000000 -0.003206 -0.003080
+14.000000 -0.002568 -0.002466
+15.000000 -0.002089 -0.002005
+16.000000 -0.001722 -0.001652
+17.000000 -0.001436 -0.001377
+18.000000 -0.001210 -0.001160
+19.000000 -0.001029 -0.000987
+20.000000 -0.000882 -0.000846
+21.000000 -0.000762 -0.000731
+22.000000 -0.000663 -0.000636
+23.000000 -0.000581 -0.000556
+24.000000 -0.000511 -0.000490
+25.000000 -0.000452 -0.000433
+26.000000 -0.000402 -0.000385
+27.000000 -0.000359 -0.000344
+28.000000 -0.000322 -0.000308
+29.000000 -0.000290 -0.000277

+1.000000 -4.357973 -7.850000
+2.000000 -0.845853 -0.525000
+3.000000 -0.255348 -0.155556
+4.000000 -0.108504 -0.065625
+5.000000 -0.055785 -0.033600
+6.000000 -0.032371 -0.019444
+7.000000 -0.020425 -0.012245
+8.000000 -0.013703 -0.008203
+9.000000 -0.009635 -0.005761
+10.000000 -0.007031 -0.004200
+11.000000 -0.005286 -0.003156
+12.000000 -0.004074 -0.002431
+13.000000 -0.003206 -0.001912
+14.000000 -0.002568 -0.001531
+15.000000 -0.002089 -0.001244
+16.000000 -0.001722 -0.001025
+17.000000 -0.001436 -0.000855
+18.000000 -0.001210 -0.000720
+19.000000 -0.001029 -0.000612
+20.000000 -0.000882 -0.000525
+21.000000 -0.000762 -0.000454
+22.000000 -0.000663 -0.000394
+23.000000 -0.000581 -0.000345
+24.000000 -0.000511 -0.000304
+25.000000 -0.000452 -0.000269
+26.000000 -0.000402 -0.000239
+27.000000 -0.000359 -0.000213
+28.000000 -0.000322 -0.000191
+29.000000 -0.000290 -0.000172

+1.000000 -4.980540 -4.981000
+2.000000 -0.966689 -0.966750
+3.000000 -0.291827 -0.286444
+4.000000 -0.124005 -0.120844
+5.000000 -0.063754 -0.061872
+6.000000 -0.036996 -0.035806
+7.000000 -0.023343 -0.022548
+8.000000 -0.015661 -0.015105
+9.000000 -0.011012 -0.010609
+10.000000 -0.008035 -0.007734
+11.000000 -0.006041 -0.005811
+12.000000 -0.004656 -0.004476
+13.000000 -0.003664 -0.003520
+14.000000 -0.002935 -0.002819
+15.000000 -0.002387 -0.002292
+16.000000 -0.001968 -0.001888
+17.000000 -0.001641 -0.001574
+18.000000 -0.001383 -0.001326
+19.000000 -0.001176 -0.001128
+20.000000 -0.001008 -0.000967
+21.000000 -0.000871 -0.000835
+22.000000 -0.000758 -0.000726
+23.000000 -0.000663 -0.000636
+24.000000 -0.000584 -0.000559
+25.000000 -0.000517 -0.000495
+26.000000 -0.000459 -0.000440
+27.000000 -0.000410 -0.000393
+28.000000 -0.000368 -0.000352
+29.000000 -0.000331 -0.000317

+1.000000 -4.980540 -7.850000
+2.000000 -0.966689 -0.525000
+3.000000 -0.291827 -0.155556
+4.000000 -0.124005 -0.065625
+5.000000 -0.063754 -0.033600
+6.000000 -0.036996 -0.019444
+7.000000 -0.023343 -0.012245
+8.000000 -0.015661 -0.008203
+9.000000 -0.011012 -0.005761
+10.000000 -0.008035 -0.004200
+11.000000 -0.006041 -0.003156
+12.000000 -0.004656 -0.002431
+13.000000 -0.003664 -0.001912
+14.000000 -0.002935 -0.001531
+15.000000 -0.002387 -0.001244
+16.000000 -0.001968 -0.001025
+17.000000 -0.001641 -0.000855
+18.000000 -0.001383 -0.000720
+19.000000 -0.001176 -0.000612
+20.000000 -0.001008 -0.000525
+21.000000 -0.000871 -0.000454
+22.000000 -0.000758 -0.000394
+23.000000 -0.000663 -0.000345
+24.000000 -0.000584 -0.000304
+25.000000 -0.000517 -0.000269
+26.000000 -0.000459 -0.000239
+27.000000 -0.000410 -0.000213
+28.000000 -0.000368 -0.000191
+29.000000 -0.000331 -0.000172

+1.000000 -5.976648 -5.977000
+2.000000 -1.160027 -1.160125
+3.000000 -0.350192 -0.343741
+4.000000 -0.148805 -0.145016
+5.000000 -0.076505 -0.074248
+6.000000 -0.044395 -0.042968
+7.000000 -0.028012 -0.027058
+8.000000 -0.018793 -0.018127
+9.000000 -0.013214 -0.012731
+10.000000 -0.009642 -0.009281
+11.000000 -0.007249 -0.006973
+12.000000 -0.005587 -0.005371
+13.000000 -0.004397 -0.004224
+14.000000 -0.003522 -0.003382
+15.000000 -0.002865 -0.002750
+16.000000 -0.002361 -0.002266
+17.000000 -0.001969 -0.001889
+18.000000 -0.001659 -0.001591
+19.000000 -0.001411 -0.001353
+20.000000 -0.001210 -0.001160
+21.000000 -0.001046 -0.001002
+22.000000 -0.000910 -0.000872
+23.000000 -0.000796 -0.000763
+24.000000 -0.000701 -0.000671
+25.000000 -0.000620 -0.000594
+26.000000 -0.000551 -0.000528
+27.000000 -0.000492 -0.000472
+28.000000 -0.000442 -0.000423
+29.000000 -0.000397 -0.000381

 // create substrate potetntial buffer
 void DeOliveiraPotentials::create_sub_buffer( void )
 {
-    logMessage("calculating substrate potential buffer [DE OLIVEIRA]...");
+    logMessage("calculating substrate potential buffer [DE OLIVEIRA]...%f %f",_C,_D);
     _f_sub_array.clear();
     _f_sub_array.resize(_N+1, 0.0);


+# fitting for ebner and de-oliveira substrate potentials
+import os, sys
+from math import tanh, log, exp, pi
+
+
+# ########## #
+# PARAMETERS #
+# ########## #
+k = 1.3806488e-23
+sigma = 3.405 # Angstrom
+sigma_w = 3.735 # Angstrom
+alpha = 1.75
+eps_k = 120.0   # K
+eps_w_k = 160.0 # K
+D = 0.846 # n_w * sigma_w**3
+
+A = pow(2.0,1.0/6)*sigma
+R = A/2
+
+eps = eps_k*k
+eps_w = eps_w_k*k
+# eps = 1.0
+# eps_w = 0.95
+N=30
+
+# ################ #
+# Global variables #
+# ################ #
+
+z1 = pow(2.0/5,1.0/6)*sigma_w
+
+def lennard_jones(e, s, r):
+    """calc the Lennard-Jones potential (6-12)"""
+    return 4*e*(pow(s/r,12)-pow(s/r,6))
+
+
+def Vo(r,e,C=4.20,D=7.85):
+    """calc the De Oliveira potential, representing approximate substrate potential
+    It has been normalized to set the minimum to z1, in order to make it comparable to ebner potential
+    """
+    if r<=z1:
+        return -D*e
+    else:
+        d = r/z1
+        return -C*e/(d*d*d)
+
+
+def Ve(z,e):
+    """calc the 9-3 potential, representing approximate substrate potential"""
+    return alpha*4*pi/45*e*D*(pow(sigma_w/z,9)-15.0/2*pow(sigma_w/z,3))
+
+
+
+def ebner_Vn(n):
+    """ calculates Vn/eps, where Vn is defined in ebner, (6)"""
+    e = eps_w_k/eps_k
+    zn = z1 + pow(2.0/3,0.5)*A*(n-1)
+    return Ve(zn,e)
+
+
+def oliveira_Vn(n,C=4.20,D=7.85):
+    """ calculates Vn/eps, where Vn is defined in de-oliveira"""
+    if n==1:
+        return -D
+    else:
+        return -C/(n*n*n)
+
+
+# calc the first 10 layer points
+zl = []
+for n in range(1,10):
+    zl.append(z1 + pow(2.0/3,0.5)*A*(n-1))
+
+
+alpha = 1.75
+
+
+def dump_potentials(fname,C=4.20,D=7.85):
+    ff = open(fname,'w')
+    for i in range(90,2000):
+        r = i/100.0
+        do = Vo(r,1.0,C,D)
+        eb = Ve(r,eps_w_k/eps_k)
+        ff.write("%f %f %f\n" %(r,eb,do))
+    ff.close()
+    plot_potentials(C,D)
+
+def plot_potentials(C,D):
+    import matplotlib.pyplot as plt
+    x = []
+    y1 = []
+    y2 = []
+    for i in range(9,200):
+        r = i/10.0
+        do = Vo(r,1.0,C,D)
+        eb = Ve(r,eps_w_k/eps_k)
+        x.append(r)
+        y1.append(do)
+        y2.append(eb)
+    plt.hold(1)
+    plt.plot(x,y1, '-r')
+    plt.plot(x,y2, '-b')
+    plt.axis([0,20,-8,6])
+    # plt.ylim(-10,8)
+    ax = plt.gca()
+    ax.set_xlabel("distance")
+    ax.set_ylabel('energy')
+    plt.vlines(zl,-10,10,'k','dashed')
+    plt.title('fitting with C=%f D=%f' %(C,D))
+    plt.show()
+
+def test1():
+    ff = open('pot_comp_A%.2f_orig.dat'%(alpha),'w')
+    for i in range(1,N):
+        do = oliveira_Vn(i)
+        eb = ebner_Vn(i)
+        ff.write("%f %f %f\n" %(i,eb,do))
+    ff.close()
+    dump_potentials('aaa.dat')
+
+def fit_de_oliveira():
+    """fit the C and D coeffs, matching de oliveira to ebner.
+    Since the De Oliveira potential is different for k=1 and for k>1, the
+    procedure performs 2 separate optimizations
+    """
+    # optimize D (k==1)
+    best_err = 1e8;
+    best_d = -1
+    for i in range(0,10000):
+        d = i/1000.0
+        err = abs(ebner_Vn(1)-oliveira_Vn(1,1.0,d))
+        if err<best_err:
+            best_err = err;
+            best_d = d
+    print 'final D=%f (error=%f)\n' % (best_d,best_err)
+    D = best_d
+    # optimize C (k>1)
+    best_err = 1e8;
+    best_c = -1
+    for i in range(0,10000):
+        c = i/1000.0
+        err = 0
+        for k in range(2,N):
+            err = err + abs(ebner_Vn(k)-oliveira_Vn(k,c,D))
+        if err<best_err:
+            best_err = err;
+            best_c = c
+    print 'final C=%f (error=%f)\n' % (best_c,best_err)
+    C = best_c
+    # write fitted potential
+    ff = open('pot_comp_A%.2f_fit.dat'%(alpha),'w')
+    for i in range(1,N):
+        do = oliveira_Vn(i,C,D)
+        eb = ebner_Vn(i)
+        ff.write("%f %f %f\n" %(i,eb,do))
+    ff.close()
+
+    dump_potentials('bbb.dat',C,D)
+
+
+# saved results:
+# alpha 1.75 --> C=6.767 D=4.358
+# alpha 2.0  --> C=7.734 D=4.981
+# alpha 2.4  --> C=9.281 D=5.977
+
+
+if __name__ == '__main__':
+    test1()
+    fit_de_oliveira()