diff --git a/src/CalcTB_laguerre.jl b/src/CalcTB_laguerre.jl index c8662e9c..10ec6301 100644 --- a/src/CalcTB_laguerre.jl +++ b/src/CalcTB_laguerre.jl @@ -528,7 +528,7 @@ function make_coefs(at_list, dim; datH=missing, datS=missing, cutoff=18.01, min_ for a in at_arr #println("a $a") - println(typeof(a)) +# println(typeof(a)) if !haskey(lim, a) lim[a] = 0.02 end @@ -568,9 +568,9 @@ function make_coefs(at_list, dim; datH=missing, datS=missing, cutoff=18.01, min_ end end - for x in [dim, datH, datS, totH, totS, data_info, inds_int, at_list, orbs, cutoff, min_dist, dist_frontier2, version, lim, repval] - println(typeof(x)) - end +# for x in [dim, datH, datS, totH, totS, data_info, inds_int, at_list, orbs, cutoff, min_dist, dist_frontier2, version, lim, repval] +# println(typeof(x)) +# end return coefs(dim, datH, datS, totH, totS, data_info, inds_int, at_list, orbs, cutoff, min_dist, dist_frontier2, version, lim, repval) @@ -3938,12 +3938,14 @@ function calc_tb_prepare_fast(reference_tbc::tb_crys; use_threebody=false, use_t cut = cutoff_fn(dist, cutoff_onX - cutoff_length, cutoff_onX) end - ad = 2.0*dist - expa=exp.(-0.5*ad) + if dist > 1e-5 + ad = 2.0*dist + expa=exp.(-0.5*ad) - rho[a1, 1] += (1.0 * expa) * cut - rho[a1, 2] += (1.0 .- ad) * expa * cut - + rho[a1, 1] += (1.0 * expa) * cut + rho[a1, 2] += (1.0 .- ad) * expa * cut + # println("ADD RHO FIT $a1 dist $dist 1 $((1.0 * expa) * cut) 2 $( (1.0 .- ad) * expa * cut) cut $cut") + end for o1 = orb2ind[a1] a1a,t1a,s1 = ind2orb[o1] @@ -4006,12 +4008,12 @@ function calc_tb_prepare_fast(reference_tbc::tb_crys; use_threebody=false, use_t for o = orb2ind[a] aa,t,s = ind2orb[o] ind = ind_conversion[(o,o,c_zero)] - println("t $t ", typeof(t)) +# println("t $t ", typeof(t)) at_set = Set((t,)) - println(at_set) - println("at_set ", at_set) - println(keys(eam_arrays)) - println("RHO $(rho[a,:]) vals $([rho[a,1]^2, rho[a,2]^2, rho[a,1]*rho[a,2]])") +# println(at_set) +# println("at_set ", at_set) +# println(keys(eam_arrays)) +# println("RHO $(rho[a,:]) vals $([rho[a,1]^2, rho[a,2]^2, rho[a,1]*rho[a,2]])") eam_arrays[at_set][1][ind,:] += [rho[a,1]^2, rho[a,2]^2, rho[a,1]*rho[a,2]] end end @@ -7892,6 +7894,8 @@ function calc_tb_LV(crys::crystal, database=missing; reference_tbc=missing, verb rho_th[a1, 1, id] += lag_arr[1] * cut_on rho_th[a1, 2, id] += lag_arr[2] * cut_on + +# println("add rho $a1 dist $dist_a 1 $(lag_arr[1] * cut_on) 2 $(lag_arr[2] * cut_on)") core!(cham, a1, a2, t1, t2, norb, orbs_arr, DAT_IND_ARR, lag_arr, DAT_ARR, cut_a, H, S, lmn_arr, sym_arr, sym_arrS) core_onsite!(c_zero, a1, a2, t1, t2, norb, orbs_arr, DAT_IND_ARR_O, lag_arr, DAT_ARR, cut_on, H, lmn_arr, sym_arr, sym_arrS) @@ -7906,11 +7910,12 @@ function calc_tb_LV(crys::crystal, database=missing; reference_tbc=missing, verb if use_eam rho = sum(rho_th, dims=3) +# println("rho $rho") for a = 1:crys.nat t = crys.stypes[a] d = database[(:eam, t)].datH temp = d[1]*rho[a,1]^2 + d[2]*rho[a,2]^2 + d[3]*rho[a,1]*rho[a,2] - println("a $a d $d rho[a,:] $(rho[a,:]) temp $temp vals $([rho[a,1]^2, rho[a,2]^2, rho[a,1]*rho[a,2]]) ") +# println("a $a d $d rho[a,:] $(rho[a,:]) temp $temp vals $([rho[a,1]^2, rho[a,2]^2, rho[a,1]*rho[a,2]]) ") for ox = 1:norb[a] oxx = orbs_arr[a,ox,1] H[ oxx, oxx, c_zero] += temp diff --git a/test/test_laguerre_eam.jl b/test/test_laguerre_eam.jl index 4589f5b8..0b74c2c4 100644 --- a/test/test_laguerre_eam.jl +++ b/test/test_laguerre_eam.jl @@ -21,31 +21,33 @@ function testEAM() tbc_list = [] @testset "testing laguerre fitting fake example" begin + atom = :H + @suppress begin - c1 = makecrys([7.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0], [:H], units="Bohr"); - c2 = makecrys([8.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0], [:H], units="Bohr"); - c3 = makecrys([8.5 0 0; 0 8.5 0; 0 0 7.0], [0 0 0], [:H], units="Bohr"); - c4 = makecrys([6.0 0 0; 0 6.0 0; 0 0 6.0], [0 0 0], [:H], units="Bohr"); + c1 = makecrys([7.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0], [atom], units="Bohr"); + c2 = makecrys([8.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0], [atom], units="Bohr"); + c3 = makecrys([8.5 0 0; 0 8.5 0; 0 0 7.0], [0 0 0], [atom], units="Bohr"); + c4 = makecrys([6.0 0 0; 0 6.0 0; 0 0 6.0], [0 0 0], [atom], units="Bohr"); - c5 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.49 0 0 ], [:H, :H], units="Bohr"); - c6 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.45 0 0 ], [:H, :H], units="Bohr"); - c7 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.40 0 0 ], [:H, :H], units="Bohr"); - c8 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.35 0 0 ], [:H, :H], units="Bohr"); - c9 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.47 0 0 ], [:H, :H], units="Bohr"); + c5 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.49 0 0 ], [atom, atom], units="Bohr"); + c6 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.45 0 0 ], [atom, atom], units="Bohr"); + c7 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.40 0 0 ], [atom, atom], units="Bohr"); + c8 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.35 0 0 ], [atom, atom], units="Bohr"); + c9 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.47 0 0 ], [atom, atom], units="Bohr"); - c4a = makecrys([5.0 0 0; 0 5.0 0; 0 0 5.0], [0 0 0], [:H], units="Bohr"); - c4b = makecrys([4.0 0 0; 0 4.0 0; 0 0 4.0], [0 0 0], [:H], units="Bohr"); - c4c = makecrys([3.0 0 0; 0 3.0 0; 0 0 3.0], [0 0 0], [:H], units="Bohr"); + c4a = makecrys([5.0 0 0; 0 5.0 0; 0 0 5.0], [0 0 0], [atom], units="Bohr"); + c4b = makecrys([4.0 0 0; 0 4.0 0; 0 0 4.0], [0 0 0], [atom], units="Bohr"); + c4c = makecrys([3.0 0 0; 0 3.0 0; 0 0 3.0], [0 0 0], [atom], units="Bohr"); begin # if true database = Dict() - database[(:H, :H)] = ThreeBodyTB.CalcTB.make_coefs(Set([:H, :H]), 2) - database[(:eam, :H)] = ThreeBodyTB.CalcTB.make_coefs(Set([:H, :H]), 0) - database[(:H, :H)].datH .= 0.0 + database[(atom, atom)] = ThreeBodyTB.CalcTB.make_coefs(Set([atom, atom]), 2) + database[(:eam, atom)] = ThreeBodyTB.CalcTB.make_coefs(Set([atom, atom]), 0) +# database[(atom, atom)].datH .= 0.0 tbc1 = ThreeBodyTB.CalcTB.calc_tb_LV(c1, database, use_threebody=false, use_threebody_onsite=false, use_eam=true); tbc2 = ThreeBodyTB.CalcTB.calc_tb_LV(c2, database, use_threebody=false, use_threebody_onsite=false, use_eam=true); @@ -62,16 +64,18 @@ function testEAM() tbc_list = [tbc1, tbc2, tbc3, tbc4, tbc5, tbc6, tbc7, tbc8, tbc9, tbc4a, tbc4b, tbc4c] - newdatabase = ThreeBodyTB.FitTB.do_fitting(tbc_list, fit_threebody=false, fit_threebody_onsite=false,do_plot=false) + newdatabase = ThreeBodyTB.FitTB.do_fitting(tbc_list, fit_threebody=false, fit_threebody_onsite=false,do_plot=false, fit_eam=true) # newdatabase = ThreeBodyTB.FitTB.do_fitting_recursive(tbc_list, fit_threebody=false, do_plot=false) # println("newdartabase") # println(newdatabase[("Li", "Li")]) - println(newdatabase[(:H, :H)].datH) - println(database[(:H, :H)].datH) - println(sum(abs.(newdatabase[(:H, :H)].datH .- database[(:H, :H)].datH))) - @test sum(abs.(newdatabase[(:H, :H)].datH .- database[(:H, :H)].datH)) ≤ 1.5e-3 +# println(newdatabase[(atom, atom)].datH) +# println(database[(atom, atom)].datH) +# println(sum(abs.(newdatabase[(atom, atom)].datH .- database[(atom, atom)].datH))) + @test sum(abs.(newdatabase[(atom, atom)].datH .- database[(atom, atom)].datH)) ≤ 1.5e-3 + + @test sum(abs.(newdatabase[(:eam, atom)].datH .- database[(:eam, atom)].datH)) ≤ 1.5e-3 end