Defined better PATH initialization for the code.

code/Makefile

@@ -4,11 +4,15 @@
 
 # Compilers
 
-FC = gfortran
-#FFLAGS = -fp-stack-check -g -O2 -traceback -check all
+# HPC Strasbourg:
+FC = ifort
 FFLAGS = -fstack-check -g -O2
-#LIBS = -L/b/home/quant/bsenjean/LAPACK/lapack-3.7.1 -llapacke -lrefblas
-LIBS = -L/usr/local/lib -llapack -lrefblas
+LIBS = -L/b/home/quant/bsenjean/LAPACK/lapack-3.7.1 -llapacke -lrefblas
+
+# Laptop:
+#FC = gfortran
+#FFLAGS = -fp-stack-check -g -O2 -traceback -check all
+#LIBS = -L/usr/local/lib -llapack -lrefblas
 
 MAIN1 = beta_and_derivatives.o
 OBJS1 = secant.o quadpack_double.o

code/run_PSOET.py

@@ -25,7 +25,6 @@
 from scipy.linalg import eigh
 from scipy import optimize
 
-
 def generate_Hamiltonian(L,N,t,v):
       """
       Function that generates the 
@@ -139,8 +138,8 @@ def persite_and_dblocc(occ,dimp,U,t,n_imp,beta,dbeta_dU,approx):
       persite = persite + (-4*t*np.sin(np.pi*occ[i]/2.0)/np.pi + t*dec_dt + U*dimp[i] + U*decbath_dU)/n_imp
     return persite, dblocc
 
-def write_FCIDUMP(U,h_emb,n_imp):
-   with open("FCIDUMP","w") as f:
+def write_FCIDUMP(U,h_emb,n_imp,work_directory):
+   with open(work_directory+"/FCIDUMP","w") as f:
       f.write("&FCI NORB= {}, NELEC= {}, MS2= 0,\n".format(2*n_imp,2*n_imp))
       f.write(" ORBSYM="+"1,"*(2*n_imp)+"\n")
       f.write("ISYM=1,\n")
@@ -156,8 +155,8 @@ def write_FCIDUMP(U,h_emb,n_imp):
          f.write('%15.10f %3d %3d %3d %3d\n' % (h_emb[i,i],i+1,i+1,0,0))
       f.close()
 
-def write_dmrg_conf(n_imp,m,opt):
-   with open("dmrg_"+str(opt)+".conf","w") as f:
+def write_dmrg_conf(n_imp,m,opt,work_directory):
+   with open(work_directory+"/dmrg_"+str(opt)+".conf","w") as f:
       f.write("sym c1\n")
       f.write("orbitals FCIDUMP\n")
       f.write("nelec {}\n".format(2*n_imp))
@@ -170,12 +169,12 @@ def write_dmrg_conf(n_imp,m,opt):
       f.write(str(opt))
       f.close()
 
-def read_dmrg_onepdm_twopdm(n_imp):
-    with open("spatial_onepdm.0.0.txt") as f:
+def read_dmrg_onepdm_twopdm(n_imp,work_directory):
+    with open(work_directory+"/spatial_onepdm.0.0.txt") as f:
      f.readline()
      onepdm_tmp = [[token for token in line.split()] for line in f.readlines()]
     f.close()
-    with open("spatial_twopdm.0.0.txt") as f:
+    with open(work_directory+"/spatial_twopdm.0.0.txt") as f:
      f.readline()
      twopdm_tmp = [[token for token in line.split()] for line in f.readlines()]
     f.close()
@@ -192,7 +191,7 @@ def read_dmrg_onepdm_twopdm(n_imp):
        dimp[i] = twopdm[i,i,i,i]
     return occ,dimp
 
-def self_consistent_loop(L,N,U,t,m,occ,beta,dbeta_dU,n_imp,approx,P,chem_pot,semi_opt,code_directory):
+def self_consistent_loop(L,N,U,t,m,occ,beta,dbeta_dU,n_imp,approx,P,chem_pot,semi_opt,code_directory,work_directory,dmrg_code):
    occ_exact = [N/(1.0*L)]*L
    old_occ = [0]*L
    delta_occ = 0
@@ -239,15 +238,20 @@ def self_consistent_loop(L,N,U,t,m,occ,beta,dbeta_dU,n_imp,approx,P,chem_pot,sem
    else:
       print("  {} impurities are considered. The embedded problem has {} sites with {} interacting sites.".format(n_imp,2*n_imp,n_imp))
       print("  It is solved numerically by Density Matrix Renormalization Group (https://github.com/sanshar/Block).\n") 
-      write_FCIDUMP(U,h_emb,n_imp)
-      write_dmrg_conf(n_imp,m,"onepdm")
-      write_dmrg_conf(n_imp,m,"twopdm")
+      write_FCIDUMP(U,h_emb,n_imp,work_directory)
+      write_dmrg_conf(n_imp,m,"onepdm",work_directory)
+      write_dmrg_conf(n_imp,m,"twopdm",work_directory)
       try:
-         subprocess.check_call(code_directory+"dmrg_script.sh",shell=True)
+#         subprocess.check_call(code_directory+"dmrg_script.sh",shell=True)
+          subprocess.check_call("mpirun -np 1 {0}/block.spin_adapted {1}/dmrg_onepdm.conf > {1}/dmrg1.out 2> {1}/dmrg1.err ; wait".format(dmrg_code,work_directory),shell=True)
+          subprocess.check_call("mpirun -np 1 {0}/block.spin_adapted {1}/dmrg_twopdm.conf > {1}/dmrg2.out 2> {1}/dmrg2.err ; wait".format(dmrg_code,work_directory),shell=True)
+          subprocess.check_call("cp node0/spatial_onepdm.0.0.txt {}".format(work_directory),shell=True)
+          subprocess.check_call("cp node0/spatial_twopdm.0.0.txt {}".format(work_directory),shell=True)
       except:
-         print("Could not run the dmrg_script.sh, check if it is in code_directory.")
+#         print("Could not run the dmrg_script.sh, check if it is in code_directory.")
+         print("Could not run dmrg properly.")
          sys.exit(0)
-      occ_imp,dimp = read_dmrg_onepdm_twopdm(n_imp)
+      occ_imp,dimp = read_dmrg_onepdm_twopdm(n_imp,work_directory)
 
    print("""   Step 4/5: Compute the SOET per-site energy and double occupation.
             Use the values of the occupation and the double occupation of the impurity
@@ -283,7 +287,7 @@ def self_consistent_loop(L,N,U,t,m,occ,beta,dbeta_dU,n_imp,approx,P,chem_pot,sem
 
 
 
-def run_psoet(L,N,U,t,n_imp,approx,code_directory,
+def run_psoet(L,N,U,t,n_imp,approx,code_directory,work_directory,dmrg_code,
               semi_opt = False,
               single_shot = False,
               chem_pot = False,
@@ -417,9 +421,9 @@ def run_psoet(L,N,U,t,n_imp,approx,code_directory,
          print("        ITERATION " + str(iteration))    
          print("#"*30+"\n")
          if chem_pot is True:
-           occ, delta_occ, persite, dblocc, persite_nexact, dblocc_nexact, mu, deltav = self_consistent_loop(L,N,U,t,m,occ,beta,dbeta_dU,n_imp,approx,P,chem_pot,semi_opt,code_directory)
+           occ, delta_occ, persite, dblocc, persite_nexact, dblocc_nexact, mu, deltav = self_consistent_loop(L,N,U,t,m,occ,beta,dbeta_dU,n_imp,approx,P,chem_pot,semi_opt,code_directory,work_directory,dmrg_code)
          else:
-           occ, delta_occ, persite, dblocc, persite_nexact, dblocc_nexact = self_consistent_loop(L,N,U,t,m,occ,beta,dbeta_dU,n_imp,approx,P,chem_pot,semi_opt,code_directory)
+           occ, delta_occ, persite, dblocc, persite_nexact, dblocc_nexact = self_consistent_loop(L,N,U,t,m,occ,beta,dbeta_dU,n_imp,approx,P,chem_pot,semi_opt,code_directory,work_directory,dmrg_code)
          f.write('%8d ' % (iteration))
          for i in range(n_imp):
             f.write('%12.8f ' % occ[i])

work_dir/Mott_transition.py

@@ -0,0 +1,26 @@
+#!/usr/bin/python 
+
+L = 400
+t = 1.0
+n_imp = 4
+approx = "iBALDA"
+description = "_singleshot"
+for U in [4.0,8.0]:
+  output = "L{}_U{}_t{}_nimp{}_{}{}_MIT_results.out".format(str(L),str(U),str(t),str(n_imp),approx,description)
+  with open("OUTPUT/"+output,"w") as f:
+     f.write("%15s %15s\n" % ("n_elec","chem_pot"))
+  for i in range(-100*int(U)-1,100*int(U)+1,1):
+    mu = i/100.0
+    energy_min = 1000
+    for N in range(2*n_imp,402,2): 
+      name = "L{}_N{}_U{}_t{}_nimp{}_{}{}.out".format(str(L),str(N),str(U),str(t),str(n_imp),approx,description)
+      with open("OUTPUT/"+name,"r") as f:
+        for i, lines in enumerate(f):
+          if i == 3:
+            energy = float(lines.split()[4+n_imp-1])
+            energy_mu_n = energy - mu*N/L
+            if energy_mu_n < energy_min:
+              energy_min = energy_mu_n
+              Nmin = N
+    with open("OUTPUT/"+output,"a") as f:
+       f.write("%15d %15.6f\n" % (Nmin, mu))

work_dir/rename_files.py

@@ -0,0 +1,15 @@
+import subprocess
+
+L=400
+t=1.0
+description="_singleshot"
+Ulist = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.5,3.0,3.5,4.0,4.5,5.0,6.0,7.0,8.0,9.0,10.0,12.0,14.0,16.0,18.0,20.0,30.0,40.0,50.0,60.0,70.0,80.0,90.0,100.0,150.0,200.0,1000.0,10000.0]
+
+for approx in ["iBALDA"]:
+ for U in Ulist:
+   for n_imp in [1,4]:
+    for N in [400]:
+     name1 = "OUTPUT/L{}_N{}_U{}_t{}_nimp{}_{}.out".format(str(L),str(N),str(U),str(t),str(n_imp),approx)
+     name2 = "OUTPUT/L{}_N{}_U{}_t{}_nimp{}_{}{}.out".format(str(L),str(N),str(U),str(t),str(n_imp),approx,description)
+     subprocess.call("mv -f " + name1 + " " + name2,shell=True)
+

work_dir/results_fctU.py

@@ -0,0 +1,33 @@
+#!/usr/bin/python 
+import os
+
+L = 400
+t = 1.0
+n_imp = 1
+approx = "2LBALDA"
+description ="_singleshot"
+Ulist = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.5,3.0,3.5,4.0,4.5,5.0,6.0,7.0,8.0,9.0,10.0,12.0,14.0,16.0,18.0,20.0,30.0,40.0,50.0,60.0,70.0,80.0,90.0,100.0,150.0,200.0,1000.0,10000.0]
+chem_pot = False
+for N in [400]:
+ output = "L{}_N{}_t{}_nimp{}_{}{}_fctU_results.out".format(str(L),str(N),str(t),str(n_imp),approx,description)
+ with open("OUTPUT/"+output,"w") as f:
+    if chem_pot is True:
+      f.write("%15s %15s %15s %15s %15s %15s %15s\n" % ("U/t","SS_energy","SS_dblocc","energy","dblocc","mu","deltav"))
+    else:
+      f.write("%15s %15s %15s %15s %15s\n" % ("U/t","SS_energy","SS_dblocc","energy","dblocc"))
+ for U in Ulist:
+    name = "L{}_N{}_U{}_t{}_nimp{}_{}{}.out".format(str(L),str(N),str(U),str(t),str(n_imp),approx,description)
+    if os.stat("OUTPUT/"+name).st_size == 0:
+       raise ValueError("File " + name + " is empty.")
+    with open("OUTPUT/"+name,"r") as f:
+       for i, lines in enumerate(f):
+         if i == 3:
+            single_shot = lines.split()
+         else:
+            pass
+       conv = os.popen("tail -n 1 %s" % "OUTPUT/"+name).read().split()
+    with open("OUTPUT/"+output,"a") as f:
+       if chem_pot is True:
+         f.write("%15f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n" % (U,float(single_shot[4+n_imp-1]),float(single_shot[6+n_imp-1]),float(conv[3+n_imp-1]),float(conv[5+n_imp-1]),float(single_shot[7+n_imp-1]),float(single_shot[8+n_imp-1])))
+       else:
+         f.write("%15f %15.8f %15.8f %15.8f %15.8f\n" % (U,float(single_shot[4+n_imp-1]),float(single_shot[6+n_imp-1]),float(conv[3+n_imp-1]),float(conv[5+n_imp-1])))

work_dir/results_fctn.py

@@ -0,0 +1,32 @@
+#!/usr/bin/python 
+import os
+
+L = 400
+t = 1.0
+n_imp = 4
+approx = "iBALDA"
+description ="_singleshot"
+chem_pot = False
+for U in [4.0,8.0]:
+ output = "L{}_U{}_t{}_nimp{}_{}{}_fctn_results.out".format(str(L),str(U),str(t),str(n_imp),approx,description)
+ with open("OUTPUT/"+output,"w") as f:
+    if chem_pot is True:
+      f.write("%15s %15s %15s %15s %15s %15s %15s\n" % ("N_elec","SS_energy","SS_dblocc","energy","dblocc","mu","deltav"))
+    else:
+      f.write("%15s %15s %15s %15s %15s\n" % ("N_elec","SS_energy","SS_dblocc","energy","dblocc"))
+ for N in  range(2*n_imp,402,2):
+    name = "L{}_N{}_U{}_t{}_nimp{}_{}{}.out".format(str(L),str(N),str(U),str(t),str(n_imp),approx,description)
+    if os.stat("OUTPUT/"+name).st_size == 0:
+       raise ValueError("File " + name + " is empty.")
+    with open("OUTPUT/"+name,"r") as f:
+       for i, lines in enumerate(f):
+         if i == 3:
+            single_shot = lines.split()
+         else:
+            pass
+       conv = os.popen("tail -n 1 %s" % "OUTPUT/"+name).read().split()
+    with open("OUTPUT/"+output,"a") as f:
+       if chem_pot is True:
+         f.write("%15d %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n" % (N,float(single_shot[4+n_imp-1]),float(single_shot[6+n_imp-1]),float(conv[3+n_imp-1]),float(conv[5+n_imp-1]),float(single_shot[7+n_imp-1]),float(single_shot[8+n_imp-1])))
+       else:
+         f.write("%15d %15.8f %15.8f %15.8f %15.8f\n" % (N,float(single_shot[4+n_imp-1]),float(single_shot[6+n_imp-1]),float(conv[3+n_imp-1]),float(conv[5+n_imp-1])))

work_dir/script.py

@@ -1,24 +1,29 @@
 #!/usr/bin/python 
 import sys
-code_directory = "/Users/bsenjean/Documents/codes/P-SOET_onHPC/P-SOET_final/code/"
+import os
+
+code_directory = os.getenv('PSOET_DIR')+"/code/"
+work_directory = os.getcwd()
+dmrg_code      = os.getenv('DMRG_code')
+
 sys.path.append(code_directory)
 import run_PSOET
 
-n_site      = 400
-n_elec      = 398
+n_site      = 12
+n_elec      = 6
 U           = 8.0
 t           = 1.0
-n_imp       = 1 # Solved analytically for n_imp = 1. DMRG impurity solver is called for n_imp > 1
+n_imp       = 2 # Solved analytically for n_imp = 1. DMRG impurity solver is called for n_imp > 1
 approx      = ["iBALDA", "2LBALDA", "SIAMBALDA"][0]
 m           = 1000 # default is 1000. Number of renormalized states in DMRG.
-semi_opt    = True # default is False. This means only dEcimp_dn(n_0) will vary, while deHxc_dn and dEHximp_dn are calculated with fixed n=N/L.
+semi_opt    = False # default is False. This means only dEcimp_dn(n_0) will vary, while deHxc_dn and dEHximp_dn are calculated with fixed n=N/L.
 chem_pot    = False # default is False. Optimized a chemical potential numerically to match the occupation between the impurity wavefunction and the KS one, like in DMET. Works for a single impurity only. 
 single_shot = False # default is False. True --> MAXITER = 1
 description = "" # default is "". Additional description added to the name of the output file. Note that semiopt, chempot and ss (single-shot) are already set by default if they are set to True.
-MAXITER     = 50 # default is 50
+MAXITER     = 2 # default is 50
 threshold   = 0.0001 # default is 0.0001
 
-run_PSOET.run_psoet(n_site,n_elec,U,t,n_imp,approx,code_directory, # this first line are mandatory arguments.
+run_PSOET.run_psoet(n_site,n_elec,U,t,n_imp,approx,code_directory,work_directory,dmrg_code, # this first line are mandatory arguments.
                     m = m,
                     semi_opt = semi_opt,
                     chem_pot = chem_pot,