Jastrow optimization script

bin/jast_opt.py

@@ -0,0 +1,124 @@
+#!/usr/bin/env python3
+
+import scipy as sp
+import scipy.optimize
+import numpy as np
+import sys
+import os
+import time
+import subprocess
+
+sys.path.insert(0,"/home/scemama/qmcchem/EZFIO/Python/")
+
+from ezfio import ezfio
+
+# PARAMETERS
+thresh = 1.e-3
+block_time = 30
+
+def main():
+    filename = sys.argv[1]
+    ezfio.set_file(filename)
+
+    def make_atom_map():
+        labels = {}
+        dimension = 0
+        for i,k in enumerate(ezfio.nuclei_nucl_label):
+            if k in labels:
+                labels[k].append(i)
+            else:
+                labels[k] = [dimension, i]
+                dimension += 1
+        atom_map = [[] for i in range(dimension)]
+        for atom in labels.keys():
+            l = labels[atom]
+            atom_map[l[0]] = l[1:]
+        return atom_map
+    atom_map = make_atom_map()
+
+    def get_params_pen():
+        d = ezfio.jastrow_jast_pen
+        return np.array([d[m[0]] for m in atom_map])
+
+    def get_energy():
+        buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', filename],
+                                         encoding='UTF-8')
+        if buffer.strip() != "":
+            buffer = buffer.splitlines()[-1]
+            _, energy, error = [float(x) for x in buffer.split()]
+            return energy, error
+        else:
+            return None, None
+
+    def set_params_pen(x):
+        x = np.abs(x)
+        y=list(ezfio.jastrow_jast_pen)
+        for i,m in enumerate(atom_map):
+            for j in m:
+                y[j] = x[i]
+        ezfio.set_jastrow_jast_pen(y)
+
+    def run_qmc():
+        subprocess.check_output(['qmcchem', 'run', filename])
+
+    def stop_qmc():
+        subprocess.check_output(['qmcchem', 'stop', filename])
+
+    def set_vmc_params():
+        subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+                                 '-m', 'VMC',
+                                 '-l', str(block_time),
+                                 '--time-step=0.3',
+                                 '--stop-time=3600',
+                                 '--norm=1.e-5',
+                                 '-w', '10',
+                                 filename])
+
+    memo_energy = {}
+    def f(x):
+        print ("x = %s"%str(x))
+        h = str(x)
+        if h in memo_energy:
+            return memo_energy[h]
+        set_params_pen(x)
+        set_vmc_params()
+        pid = os.fork()
+        if pid == 0:
+            run_qmc()
+            os._exit(os.EX_OK)
+        else:
+            import atexit
+            atexit.register(stop_qmc)
+
+            err = thresh+1.
+            time.sleep(block_time/2)
+            while err > thresh:
+                time.sleep(block_time)
+                energy, err = get_energy()
+                if energy is not None:
+                    print(" %f  %f"%(energy, err))
+                else:
+                    err = thresh+1.
+
+                # Check if PID is still running
+                try: os.kill(pid,0)
+                except OSError: err = 0.
+            stop_qmc()
+        os.wait()
+        memo_energy[h] = energy
+        return energy
+
+
+    def run():
+        x = get_params_pen()
+        opt = sp.optimize.minimize(f,x,method="Powell",
+                             options= {'disp':True, 'ftol':thresh,'xtol':0.02})
+        print("x = "+str(opt))
+        set_params_pen(opt['x'])
+    run()
+
+if __name__ == '__main__':
+    main()
+
+
+

ocaml/Qmcchem_dataserver.ml

@@ -65,7 +65,7 @@ let run ?(daemon=true) ezfio_filename =
 
 
   (** Maximum size of the blocks file before compressing *)
-  let max_file_size = ref ( 64 * 1024 )
+  let max_file_size = ref ( 64/64 * 1024 )
   in
 
 

ocaml/Random_variable.ml

@@ -678,7 +678,6 @@ let compress_files () =
          | Not_found -> true
        )
     |> List.rev_map (fun x -> dir_name^x)
-    |> List.rev
   in
 
   let out_channel_dir = 

src/JASTROW/jastrow_param.irp.f

@@ -40,7 +40,7 @@ BEGIN_PROVIDER [ real, jast_a_up_up ]
 ! a_{up up} parameters of the Jastrow
   END_DOC
   include '../types.F'
-  jast_a_up_up = 0.25
+  jast_a_up_up = 0.5
   call get_jastrow_jast_a_up_up(jast_a_up_up)
 
 END_PROVIDER
@@ -62,7 +62,7 @@ BEGIN_PROVIDER [ real, jast_b_up_up ]
 ! b_{up up} parameters of the Jastrow
   END_DOC
   include '../types.F'
-  jast_b_up_up = 5.
+  jast_b_up_up = 1.
   call get_jastrow_jast_b_up_up(jast_b_up_up)
 
 END_PROVIDER
@@ -73,7 +73,7 @@ BEGIN_PROVIDER [ real, jast_b_up_dn ]
 ! b_{up dn} parameters of the Jastrow
   END_DOC
   include '../types.F'
-  jast_b_up_dn = 5.
+  jast_b_up_dn = 1.
   call get_jastrow_jast_b_up_dn(jast_b_up_dn)
 
 END_PROVIDER
@@ -106,7 +106,7 @@ BEGIN_PROVIDER [ real, jast_eeN_e_b, (nucl_num) ]
 ! b parameters of the electron-electron-Nucleus component of the Jastrow
   END_DOC
   include '../types.F'
-  jast_eeN_e_b(:) = 3.
+  jast_eeN_e_b(:) = 1.
   call get_jastrow_jast_eeN_e_b(jast_eeN_e_b)
 
 END_PROVIDER

src/ZMQ/sockets.irp.f

@@ -6,7 +6,7 @@ use f77_zmq
 BEGIN_PROVIDER [ character*(48), dataserver_address ]
   implicit none
   BEGIN_DOC  
-! Adderss of the data server
+! Address of the data server
   END_DOC
   dataserver_address = trim(http_server)
   integer :: i