Merge branch 'master' of github.com:scemama/quantum_package

README.md

@@ -41,9 +41,9 @@ Demo
 
 ### 1) Configure
 
-    $ ./configure <config_file> 
+    $ ./configure.py <config_file> 
     
-For example you can type `./configure config/gfortran.cfg`
+For example you can type `./configure.py config/gfortran.cfg`
 
 This command has two purposes :
 

configure

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 """configure
 
@@ -242,7 +242,7 @@ def checking(d_dependency):
                 version = check_output("irpf90 -v".split()).strip()
 
                 from distutils.version import LooseVersion
-                if LooseVersion(version) < LooseVersion("1.6.7"):
+                if LooseVersion(version) < LooseVersion("1.7.2"):
                     return 0
                 else:
                     return a

ocaml/myocamlbuild.ml

@@ -0,0 +1,13 @@
+open Ocamlbuild_plugin;;
+open Command;;
+
+dispatch begin function
+  | Before_rules -> 
+    begin
+    end
+  | After_rules ->
+    begin
+      flag ["ocaml";"compile";"native";"gprof"] (S [ A "-p"]);
+    end
+  | _ -> ()
+end

ocaml/qptypes_generator.ml

@@ -214,6 +214,26 @@ end = struct
     | _ -> raise (Invalid_argument (\"Wrong IO type : \"^s))
 
 end
+
+module Perturbation : sig
+  type t [@@deriving sexp]
+  val to_string : t -> string
+  val of_string : string -> t
+end = struct
+  type t = 
+  | EN
+  | Barycentric
+  [@@deriving sexp]
+
+  let to_string = function
+  | EN -> \"EN\"
+  | Barycentric -> \"Barycentric\"
+  let of_string  s = 
+    match (String.lowercase_ascii s) with
+    | \"en\"  -> EN
+    | \"barycentric\" -> Barycentric
+    | _ -> raise (Invalid_argument (\"Wrong Perturbation type : \"^s))
+end
 "
 
 

plugins/All_singles/H_apply.irp.f

@@ -1,5 +1,5 @@
 use bitmasks
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import *
 
 s = H_apply("just_1h_1p")

plugins/CID/H_apply.irp.f

@@ -1,7 +1,7 @@
 ! Generates subroutine H_apply_cid
 ! ----------------------------------
 
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import H_apply
 H = H_apply("cid",do_double_exc=True,do_mono_exc=False)
 print H

plugins/CID_selected/H_apply.irp.f

@@ -1,5 +1,5 @@
 use bitmasks
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import *
 from perturbation import perturbations
 
@@ -16,7 +16,7 @@ subroutine H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st
   integer, intent(in)            :: N_st
   double precision, intent(inout):: pt2(N_st), norm_pert(N_st), H_pert_diag(N_st)
 
-BEGIN_SHELL [  /usr/bin/env python ]
+BEGIN_SHELL [  /usr/bin/env python2 ]
 from perturbation import perturbations
 
 for perturbation in perturbations:

plugins/CIS/H_apply.irp.f

@@ -1,7 +1,7 @@
 ! Generates subroutine H_apply_cisd
 ! ----------------------------------
 
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import H_apply
 H = H_apply("cis",do_double_exc=False)
 print H

plugins/CISD/H_apply.irp.f

@@ -1,7 +1,7 @@
 ! Generates subroutine H_apply_cisd
 ! ----------------------------------
 
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import H_apply
 H = H_apply("cisd")
 print H

plugins/CISD_SC2_selected/H_apply.irp.f

@@ -1,5 +1,5 @@
 use bitmasks
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import *
 from perturbation import perturbations
 

plugins/CISD_selected/H_apply.irp.f

@@ -1,5 +1,5 @@
 use bitmasks
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import *
 from perturbation import perturbations
 
@@ -16,7 +16,7 @@ subroutine H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st
   integer, intent(in)            :: N_st
   double precision, intent(inout):: pt2(N_st), norm_pert(N_st), H_pert_diag(N_st)
 
-BEGIN_SHELL [  /usr/bin/env python ]
+BEGIN_SHELL [  /usr/bin/env python2 ]
 from perturbation import perturbations
 
 for perturbation in perturbations:

plugins/FOBOCI/H_apply.irp.f

@@ -1,5 +1,5 @@
 use bitmasks
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import *
 
 s = H_apply("just_1h_1p")

plugins/Full_CI/H_apply.irp.f

@@ -1,5 +1,5 @@
 use bitmasks
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import *
 
 s = H_apply("FCI")

plugins/Full_CI_ZMQ/energy.irp.f

@@ -12,14 +12,21 @@ BEGIN_PROVIDER [ double precision, pt2_E0_denominator, (N_states) ]
  ! E0 in the denominator of the PT2
  END_DOC
  if (initialize_pt2_E0_denominator) then
-  pt2_E0_denominator(1:N_states) = psi_energy(1:N_states)
+   if (h0_type == "EN") then
+     pt2_E0_denominator(1:N_states) = psi_energy(1:N_states)
+   else if (h0_type == "Barycentric") then
+     pt2_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
+   else
+     print *,  h0_type, ' not implemented'
+     stop
+   endif
 ! call ezfio_get_full_ci_zmq_energy(pt2_E0_denominator(1))
 ! pt2_E0_denominator(1) -= nuclear_repulsion
 ! pt2_E0_denominator(1:N_states) = HF_energy - nuclear_repulsion
-! pt2_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
   call write_double(6,pt2_E0_denominator(1)+nuclear_repulsion, 'PT2 Energy denominator')
  else
    pt2_E0_denominator = -huge(1.d0)
  endif
 END_PROVIDER
 
+

plugins/MP2/H_apply.irp.f

@@ -1,5 +1,5 @@
 use bitmasks
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import *
 from perturbation import perturbations
 

plugins/MRCC_Utils/H_apply.irp.f

@@ -1,5 +1,5 @@
 use bitmasks
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import *
 
 s = H_apply("mrcc")

plugins/MRPT_Utils/H_apply.irp.f

@@ -1,5 +1,5 @@
 use bitmasks
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import *
 
 s = H_apply("mrpt")

plugins/Perturbation/EZFIO.cfg

@@ -62,4 +62,9 @@ doc: Thresholds on selectors (fraction of the norm) for final PT2 calculation
 interface: ezfio,provider,ocaml
 default: 1.
 
+[h0_type]
+type: Perturbation
+doc: Type of zeroth-order Hamiltonian [ EN | Barycentric ]
+interface: ezfio,provider,ocaml
+default: EN
 

plugins/Perturbation/perturbation.irp.f

@@ -1,4 +1,4 @@
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 from perturbation import perturbations
 import os
 

plugins/Perturbation/perturbation.template.f

@@ -1,4 +1,4 @@
-BEGIN_SHELL [ /usr/bin/env python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 import perturbation
 END_SHELL
 

plugins/QMC/qp_convert_qmcpack_to_ezfio.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python2
 
 print "#QP -> QMCPACK"
 
@@ -333,11 +333,14 @@ if do_pseudo:
 # |_/ (/_ |_ 
 #
 
+
 psi_coef = ezfio.get_determinants_psi_coef()
 psi_det = ezfio.get_determinants_psi_det()
 bit_kind = ezfio.get_determinants_bit_kind()
 
 
+nexcitedstate = ezfio.get_determinants_n_states()
+
 print ""
 print "BEGIN_DET"
 print ""
@@ -349,7 +352,11 @@ if "QP_STATE" in os.environ:
   state = int(os.environ["QP_STATE"])-1
 else:
   state = 0
-psi_coef = psi_coef[state]
+
+psi_coef_small = psi_coef[state]
+
+
+
 
 encode = 8*bit_kind
 
@@ -359,11 +366,35 @@ def bindigits(n, bits):
 
 decode = lambda det: ''.join(bindigits(i,encode)[::-1] for i in det)[:mo_num]
 
-for coef, (det_a, det_b) in zip(psi_coef, psi_det):
+MultiDetAlpha = []
+MultiDetBeta = []
+for coef, (det_a, det_b) in zip(psi_coef_small, psi_det):
 
         print coef
-        print decode(det_a)
+        MyDetA=decode(det_a)
-        print decode(det_b)
+        MyDetB=decode(det_b)
+        print MyDetA 
+        print MyDetB
         print ''
-
+        MultiDetAlpha.append( det_a  )
+	MultiDetBeta.append( det_b  )
 print "END_DET"
+
+import h5py
+H5_qmcpack=h5py.File('MultiDet.h5','w')
+groupMultiDet=H5_qmcpack.create_group("MultiDet")
+groupMultiDet.create_dataset("NbDet",(1,),dtype="f8",data=len(psi_coef_small))
+
+groupMultiDet.create_dataset("Coeff",(len(psi_coef_small),),dtype="f8",data=psi_coef)
+groupMultiDet.create_dataset("nstate",(1,),dtype="i4",data=len(MyDetA))
+groupMultiDet.create_dataset("nexcitedstate",(1,),dtype="i4",data=nexcitedstate)
+groupMultiDet.create_dataset("Nbits",(1,),dtype="i4",data=len(det_a))
+
+print "temp=",MultiDetAlpha[0]
+mylen="S"+str(len(MyDetA))
+groupMultiDet.create_dataset("CI_Alpha",(len(psi_coef_small),len(det_a)),dtype='i8',data=MultiDetAlpha)
+
+mylen="S"+str(len(MyDetB))
+groupMultiDet.create_dataset("CI_Beta",(len(psi_coef_small),len(det_b)),dtype='i8',data=MultiDetBeta)
+
+H5_qmcpack.close()

plugins/analyze_wf/occupation.irp.f

@@ -8,13 +8,17 @@ subroutine get_occupation_from_dets(istate,occupation)
   integer                        :: i,j, ispin
   integer                        :: list(N_int*bit_kind_size,2)
   integer                        :: n_elements(2)
-  double precision               :: c
+  double precision               :: c, norm_2
   ASSERT (istate > 0) 
   ASSERT (istate <= N_states) 
   
   occupation = 0.d0
+  double precision, external :: u_dot_u
+
+  norm_2 = 1.d0/u_dot_u(psi_coef(1,istate),N_det)
+
   do i=1,N_det
-    c = psi_coef(i,istate)*psi_coef(i,istate)
+    c = psi_coef(i,istate)*psi_coef(i,istate)*norm_2
     call bitstring_to_list_ab(psi_det(1,1,i), list, n_elements, N_int)
     do ispin=1,2
       do j=1,n_elements(ispin)

plugins/pyscf/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+

plugins/pyscf/PyscfToQp.py

@@ -0,0 +1,137 @@
+
+import numpy,re,sys
+
+def pyscf2QP(cell,mf, kpts=[], int_threshold = 1E-15):
+   # The integral will be not printed in they are bellow that
+
+
+  PBC=False
+  ComputeMode= re.split('[. ]', str(mf))
+  print 'ComputeMode=',ComputeMode
+
+  for n in ComputeMode:
+     if n in ("UHF","KUHF","UKS"):
+           sys.exit('Unrestricted calculation unsupported in Quantum Package')
+     if n == "pbc":
+           PBC=True
+
+  if PBC and len(kpts) == 0:
+	sys.exit("ERROR (read!): You need to specify explicit the list of K-point (including gamma)")
+
+  print 'Performing PBC?:',PBC
+  if PBC:
+      from pyscf.pbc import ao2mo
+      from pyscf.pbc import tools
+  else:
+      from pyscf import ao2mo
+
+  natom = len(cell.atom_coords())
+  print 'n_atom',   natom
+  print 'num_elec', cell.nelectron
+  print 'nucl_num', len(cell.atom_coords())
+
+
+  print ''
+  mo_coeff = mf.mo_coeff # List of mo_coeff for each k-point
+  if not PBC:
+      nmo = mo_coeff.shape[1]
+  else:
+    nmo = mo_coeff[0].shape[1]
+  
+
+  # Wrote all the parameter need to creat a dummy EZFIO folder who will containt the integral after.
+  # More an implentation detail than a real thing
+  with open('param','w') as f:
+        f.write(' '.join(map(str,(cell.nelectron, nmo, natom))))
+  #                             _                             
+  # |\ |      _ |  _   _. ._   |_)  _  ._      |  _ o  _  ._  
+  # | \| |_| (_ | (/_ (_| |    | \ (/_ |_) |_| | _> | (_) | | 
+  #                                    |                      
+  
+  print 'mf, cell', mf.energy_nuc(), cell.energy_nuc()
+  shift = tools.pbc.madelung(cell, numpy.zeros(3))*cell.nelectron * -.5 if PBC else 0
+  e_nuc = cell.energy_nuc() + shift
+
+  print 'nucl_repul', e_nuc
+  with open('e_nuc','w') as f:
+         f.write(str(e_nuc))
+ 
+ 
+  from itertools import product
+ 
+  # ___                                              
+  #  |  ._ _|_  _   _  ._ _. |  _   |\/|  _  ._   _  
+  # _|_ | | |_ (/_ (_| | (_| | _>   |  | (_) | | (_) 
+  #                 _|                              
+  
+  if PBC:
+        h_ao = ('kinetic', mf.get_hcore(kpts=kpts) ) # Give only one k point ?
+        dummy_ao = ('nuclear', numpy.zeros( (len(kpts),nmo,nmo), dtype=numpy.float ))
+  else:
+        h_ao = ('kinetic', mf.get_hcore() )
+        dummy_ao = ('nuclear', numpy.zeros( (nmo,nmo), dtype=numpy.float ))
+    
+  def gen_mono_MO(mo_coeff,l_int,shift=0):
+         # 2Id transfortion Transformation. For now we handle only one or zero K point.
+         print 'l_int.shape=',l_int.shape
+ 
+         l_int_mo = reduce(numpy.dot, (mo_coeff.T, l_int, mo_coeff)) #This formula is only right for one kpt.
+ 
+	 print 'l_int_mo=',l_int_mo
+
+         for i,j in product(range(nmo), repeat=2):
+                 int_ = l_int_mo[i,j]
+                 yield (i+1+shift,j+1+shift, int_)
+ 
+  # Print 
+  for name, ao in (h_ao,dummy_ao):
+         with open('%s_mo' % name,'w') as f:
+		 print '%s_mo' % name
+		 if not PBC:
+                 	for mono in gen_mono_MO(mo_coeff,ao):
+                        	 f.write('%s %s %s\n'% mono)
+ 		 else:
+			for i,(m,a) in enumerate(zip(mo_coeff,ao)):
+			   for mono in gen_mono_MO(m,a,i):
+                                 f.write('%s %s %s\n'% mono)
+
+  # ___                              _    
+  #  |  ._ _|_  _   _  ._ _. |  _   |_) o 
+  # _|_ | | |_ (/_ (_| | (_| | _>   |_) | 
+  #                 _|                    
+  #
+
+  def ao2mo_amazing(mo_coeff):
+	if PBC:
+		eri_4d= mf.with_df.ao2mo(mo_coeff,compact=False)
+	else:
+		eri_4d= ao2mo.kernel(cell,mo_coeff,compact=False)
+	
+	return eri_4d.reshape((nmo,)*4)
+
+
+  def write_amazing(eri_4d, shift=0):
+
+    # HANDLE 8 FOLD by Scemama way. Maybe we can use compact=True
+    for l in range(nmo):
+      for k in range(nmo):
+        for j in range(l,nmo):
+                for i in range(max(j,k),nmo):
+                        v = eri_4d[i,k,j,l]
+                        if abs(v) > int_threshold:
+                              f.write('%s %s %s %s %s\n' % (i+1+shift,j+1+shift,k+1+shift,l+1+shift,v))
+
+
+  if PBC: 
+     eri_4d= mf.with_df.ao2mo(mo_coeff[0],compact=False)
+  else: #Molecular
+     eri_4d= ao2mo.kernel(cell,mo_coeff,compact=False)
+
+  eri_4d = eri_4d.reshape((nmo,)*4)
+
+  f = open('bielec_mo','w')
+  for i,mc in enumerate(mo_coeff):
+	eri = ao2mo_amazing(mc)
+	write_amazing(eri, nmo*i)
+
+

plugins/pyscf/README.rst

@@ -0,0 +1,21 @@
+=====
+pyscf
+=====
+Converter from Pyscf to Quatum Package for Molecules AND Solids
+Import this script in your Pyscf input. 
+
+
+Use as follow:
+```
+   from MolPyscfToQP import pyscf2QP
+   pyscf2QP(cell,mf,kpts=kpts,int_threshold = 1E-15)
+
+```
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.

plugins/pyscf/pyscf.main.irp.f

@@ -0,0 +1,38 @@
+program pyscf
+  implicit none
+  BEGIN_DOC
+! TODO
+  END_DOC
+  print *, '  _/                              '
+  print *, ' -:\_?,     _Jm####La             '
+  print *, 'J"(:" >  _]#AZ#Z#UUZ##,           '
+  print *, '_,::./   %(|i%12XmX1*1XL      _?, '
+  print *, '  \..\ _\(vmWQwodY+ia%lnL  _",/ ( '
+  print *, '   .:< ]J=mQD?WXn<uQWmmvd, -.-:=!'
+  print *, '   "{Z jC]QW|=3Zv)Bi3BmXv3  =   _7'
+  print *, '    ]h[Z6)WQ;)jZs]C;|$BZv+, : ./  '
+  print *, '    -#sJX%$Wmm#ev]hinW#Xi:` c ;   '
+  print *, '     #X#X23###1}vI$WWmX1>|,)nr"   '
+  print *, '     4XZ#Xov1v}=)vnXAX1nnv;1n"    '
+  print *, '     ]XX#ZXoovvvivnnnlvvo2*i7     '
+  print *, '     "23Z#1S2oo2XXSnnnoSo2>v"     '
+  print *, '      miX#L -~`""!!1}oSoe|i7      '
+  print *, '      4cn#m,        v221=|v[      '
+  print *, '      ]hI3Zma,;..__wXSe=+vo       '
+  print *, '      ]Zov*XSUXXZXZXSe||vo2       '
+  print *, '      ]Z#><iiii|i||||==vn2(       '
+  print *, '      ]Z#i<ii||+|=||=:{no2[       '
+  print *, '      ]ZUsiiiiivi|=||=vo22[       '
+  print *, '      ]XZvlliiIi|i=|+|vooo        '
+  print *, '      =v1llli||||=|||||lii(       '
+  print *, '      ]iillii||||||||=>=|<        '
+  print *, '      -ziiiii||||||+||==+>        '
+  print *, '       -%|+++||=|=+|=|==/         '
+  print *, '        -a>====+|====-:-          '
+  print *, '          "~,- --   /-            '
+  print *, '            -.     )>             '
+  print *, '           .~      +-             '
+  print *, '           . .... : .             '
+  print *, '            -------~              '
+  print *, ''
+end

plugins/read_integral/Gen_Ezfio_from_integral.sh

@@ -0,0 +1,17 @@
+#!/bin/bash
+
+ezfio=$1
+# Create the integral
+echo 'Create Integral'
+
+echo 'Create EZFIO'
+read nel nmo natom <<< $(cat param) 
+read e_nucl <<< $(cat e_nuc)
+./create_ezfio.py $ezfio $nel $natom $nmo $e_nucl
+#Handle the orbital consitensy check
+qp_edit -c $ezfio &> /dev/null
+cp $ezfio/{ao,mo}_basis/ao_md5 
+
+#Read the integral
+echo 'Read Integral'
+qp_run read_integrals_mo $ezfio 

plugins/read_integral/README.rst

@@ -3,6 +3,7 @@ read_integral
 =============
 
 Warning: CAN NOT CHANGE THE NUMBER OF MO !
+Scripts to read integrals and metadata and generates fake ezfio
 
 Needed Modules
 ==============

plugins/read_integral/create_ezfio.py

@@ -0,0 +1,48 @@
+#!/usr/bin/env python
+from ezfio import ezfio
+
+import sys
+filename = sys.argv[1]
+num_elec, nucl_num, mo_tot_num = map(int,sys.argv[2:5])
+
+nuclear_repulsion = float(sys.argv[5])
+ezfio.set_file(filename)
+
+#Important !
+import math
+ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.))
+ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.))
+
+
+#Important
+ezfio.set_nuclei_nucl_num(nucl_num)
+ezfio.set_nuclei_nucl_charge([0.]*nucl_num)
+ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num )
+ezfio.set_nuclei_nucl_label( ['He'] * nucl_num )
+
+ezfio.set_nuclei_disk_access_nuclear_repulsion('Read')
+ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion)
+
+# Ao num
+ao_num = mo_tot_num
+ezfio.set_ao_basis_ao_basis("Dummy one. We read MO")
+ezfio.set_ao_basis_ao_num(ao_num)
+ezfio.set_ao_basis_ao_nucl([1]*ao_num) #Maybe put a realy incorrect stuff
+
+#Just need one
+ao_prim_num_max  = 5
+
+d = [ [0] *ao_prim_num_max]*ao_num
+ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num)
+ezfio.set_ao_basis_ao_power(d)
+ezfio.set_ao_basis_ao_coef(d)
+ezfio.set_ao_basis_ao_expo(d)
+
+#Dummy one
+ao_md5 = '3b8b464dfc95f282129bde3efef3c502'
+ezfio.set_ao_basis_ao_md5(ao_md5)
+ezfio.set_mo_basis_ao_md5(ao_md5)
+
+
+ezfio.set_mo_basis_mo_tot_num(mo_tot_num)
+ezfio.set_mo_basis_mo_coef([ [0]*mo_tot_num] * ao_num)

plugins/read_integral/read_integrals_achocol.irp.f

@@ -0,0 +1,47 @@
+program read_integrals
+
+  PROVIDE ezfio_filename
+  call ezfio_set_integrals_monoelec_disk_access_ao_one_integrals("None")
+  call run
+end
+
+subroutine run
+  use map_module
+  implicit none
+  
+  integer :: iunit
+  integer :: getunitandopen
+
+  integer ::i,j,k,l
+  double precision :: integral
+  double precision, allocatable :: A(:,:)
+
+  integer             :: n_integrals 
+  integer(key_kind), allocatable   :: buffer_i(:) 
+  real(integral_kind), allocatable :: buffer_values(:)
+  integer(key_kind)  :: key
+   
+  call ezfio_set_integrals_monoelec_disk_access_ao_one_integrals("Read")
+
+  allocate(buffer_i(ao_num**4/8), buffer_values(ao_num**4/8))
+   
+  iunit = getunitandopen('bielec_ao','r')
+  n_integrals=0
+  do 
+    read (iunit,*,end=13) i,j,k,l, integral
+    n_integrals += 1
+    call bielec_integrals_index(i, j, k, l, buffer_i(n_integrals) )
+    buffer_values(n_integrals) = integral
+  enddo
+  13 continue
+  close(iunit)
+  
+  call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values)
+
+  call map_sort(ao_integrals_map)
+  call map_unique(ao_integrals_map)
+
+  call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
+  call ezfio_set_integrals_bielec_disk_access_ao_integrals('Read')
+
+end

plugins/read_integral/read_integrals_mo_chocol.irp.f

@@ -0,0 +1,86 @@
+program read_integrals
+  BEGIN_DOC
+! Reads the integrals from the following files:
+! - kinetic_mo
+! - nuclear_mo
+! - bielec_mo
+  END_DOC
+
+  integer :: iunit
+  integer :: getunitandopen
+  integer :: i,j,n
+
+  PROVIDE ezfio_filename
+  call ezfio_set_integrals_monoelec_disk_access_mo_one_integrals("None")
+
+  logical :: has
+  call ezfio_has_mo_basis_mo_tot_num(has)
+  if (.not.has) then
+
+    iunit = getunitandopen('nuclear_mo','r')
+    n=0
+    do 
+      read (iunit,*,end=12) i
+      n = max(n,i)
+    enddo
+    12 continue
+    close(iunit)
+    call ezfio_set_mo_basis_mo_tot_num(n)
+
+    call ezfio_has_ao_basis_ao_num(has)
+    mo_label = "None"
+    if (has) then
+      call huckel_guess
+    else
+      call ezfio_set_ao_basis_ao_num(n)
+    endif
+  endif
+  call run
+end
+
+subroutine run
+  use map_module
+  implicit none
+  
+  integer :: iunit
+  integer :: getunitandopen
+
+  integer ::i,j,k,l
+  double precision :: integral
+  double precision, allocatable :: A(:,:)
+
+  integer             :: n_integrals 
+  integer(key_kind), allocatable   :: buffer_i(:) 
+  real(integral_kind), allocatable :: buffer_values(:)
+  integer(key_kind)  :: key
+   
+  call ezfio_get_mo_basis_mo_tot_num(mo_tot_num)
+
+  allocate (A(mo_tot_num,mo_tot_num))
+  A = 0.d0
+  
+  iunit = getunitandopen('kinetic_mo','r')
+  do 
+    read (iunit,*,end=10) i,j, integral
+    A(i,j) = integral
+  enddo
+  10 continue
+  close(iunit)
+  call write_one_e_integrals('mo_kinetic_integral', A, size(A,1), size(A,2))
+
+
+  iunit = getunitandopen('nuclear_mo','r')
+  do 
+    read (iunit,*,end=12) i,j, integral
+    A(i,j) = integral
+  enddo
+  12 continue
+  close(iunit)
+  call write_one_e_integrals('mo_ne_integral', A, size(A,1), size(A,2))
+
+  call write_one_e_integrals('mo_pseudo_integral', mo_pseudo_integral,&
+        size(mo_pseudo_integral,1), size(mo_pseudo_integral,2))
+
+
+  call ezfio_set_integrals_monoelec_disk_access_mo_one_integrals("Read")
+end

plugins/shiftedbk/EZFIO.cfg

@@ -0,0 +1,22 @@
+[energy]
+type: double precision
+doc: Calculated energy
+interface: ezfio
+
+[thresh_dressed_ci]
+type: Threshold
+doc: Threshold on the convergence of the dressed CI energy
+interface: ezfio,provider,ocaml
+default: 1.e-5
+
+[n_it_max_dressed_ci]
+type: Strictly_positive_int
+doc: Maximum number of dressed CI iterations
+interface: ezfio,provider,ocaml
+default: 10
+
+[h0_type]
+type: Perturbation
+doc: Type of zeroth-order Hamiltonian [ EN | Barycentric ]
+interface: ezfio,provider,ocaml
+default: EN

plugins/shiftedbk/shifted_bk_routines.irp.f

@@ -82,11 +82,14 @@ BEGIN_PROVIDER [ double precision, E0_denominator, (N_states) ]
   ! E0 in the denominator of the PT2
   END_DOC
   if (initialize_E0_denominator) then
-   E0_denominator(1:N_states) = psi_energy(1:N_states)
+   if (h0_type == "EN") then
- ! call ezfio_get_full_ci_zmq_energy(pt2_E0_denominator(1))
+     E0_denominator(1:N_states) = psi_energy(1:N_states)
- ! pt2_E0_denominator(1) -= nuclear_repulsion
+   else if (h0_type == "Barycentric") then
- ! pt2_E0_denominator(1:N_states) = HF_energy - nuclear_repulsion
+     E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
- ! pt2_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
+   else
+     print *,  h0_type, ' not implemented'
+     stop
+   endif
   else
     E0_denominator = -huge(1.d0)
   endif

scripts/compilation/cache_compile.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 """
 Save the .o from a .f90
 and is the .o is asked a second time, retur it

scripts/compilation/qp_create_ninja.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 """
 Usage: qp_create_ninja.py create <config_file> (--development | --production)
@@ -25,7 +25,6 @@ except ImportError:
                                       "quantum_package.rc"))
 
     print "\n".join(["", "Error:", "source %s" % f, ""])
-    raise
     sys.exit(1)
 
 
@@ -820,13 +819,8 @@ if __name__ == "__main__":
     pickle_path = os.path.join(QP_ROOT, "config", "qp_create_ninja.pickle")
 
     if arguments["update"]:
-        try:
             with open(pickle_path, 'rb') as handle:
                 arguments = pickle.load(handle)
-        except IOError:
-            print "You need to create first my friend"
-            raise
-            sys.exit(1)
 
     elif arguments["create"]:
 
@@ -928,7 +922,6 @@ if __name__ == "__main__":
                      "- Or install a module that needs {0} with a main "]
 
             print "\n".join(l_msg).format(module.rel)
-            raise
             sys.exit(1)
 
     # ~#~#~#~#~#~#~#~#~#~#~#~ #

scripts/compilation/read_compilation_cfg.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 
 import os

scripts/entanglement.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 import sys                
 import matplotlib.pyplot as plt

scripts/ezfio_interface/ei_handler.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 """
 Welcome to the ei_handler.

scripts/ezfio_interface/ezfio_generate_ocaml.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 """
 This program generates all the
 OCaml templates needed by qp_edit

scripts/ezfio_interface/ezfio_generate_provider.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 
 __author__ = "Applencourt PEP8"
 __date__ = "jeudi 26 mars 2015, 12:49:35 (UTC+0100)"

scripts/ezfio_interface/qp_convert_output_to_ezfio.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 """
 convert output of gamess/GAU$$IAN to ezfio
 

scripts/generate_h_apply.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 
 import os
 

scripts/module/module_handler.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 """
 Module utilitary

scripts/module/qp_module.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 """
 Usage:

scripts/module/qp_update_readme.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 """
 Updates the README.rst of a module

scripts/perturbation.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 
 import os
 

scripts/pseudo/elts_num_ele.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 
 import os
 

scripts/pseudo/put_pseudo_in_ezfio.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 """
 Create the pseudo potential for a given atom

scripts/qp_bitmasks.py

@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python2
 
 BIT_KIND_SIZE=64
 

scripts/qp_set_frozen_core.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 
 import os
 import sys

scripts/save_current_mos.sh

@@ -24,7 +24,7 @@ if [[ ! -f "${EZFIO}/mo_basis/mo_label" ]]
 then
   LABEL='no_label'
 else
-  LABEL=$(head -1 "${EZFIO}/mo_basis/mo_label")
+  LABEL=$(head -1 "${EZFIO}/mo_basis/mo_label" | xargs) #xargs trims the result
 fi
 
 DESTINATION="save/mo_basis/${LABEL}"

scripts/utility/get_groups.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 
 import urllib

scripts/utility/is_master_repository.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 
 import subprocess
 pipe = subprocess.Popen("git config --get remote.origin.url", \

scripts/utility/qp_path.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2
 # -*- coding: utf-8 -*-
 
 import os

src/DavidsonUndressed/guess_lowest_state.irp.f

@@ -72,7 +72,7 @@ program first_guess
   call write_int(6,psi_det_size,'psi_det_size')
   TOUCH psi_det_size
 
-BEGIN_SHELL [ /usr/bin/python ]
+BEGIN_SHELL [ /usr/bin/env python2 ]
 
 template_alpha_ext = """
 do %(i2)s = %(i1)s-1,1,-1

src/Determinants/energy.irp.f

@@ -0,0 +1,31 @@
+BEGIN_PROVIDER [ double precision, diagonal_H_matrix_on_psi_det, (N_det) ]
+ implicit none
+ BEGIN_DOC
+ ! Diagonal of the Hamiltonian ordered as psi_det
+ END_DOC
+ double precision, external     :: diag_h_mat_elem
+ integer                        :: i
+ 
+ do i=1,N_det
+   diagonal_H_matrix_on_psi_det(i) = diag_h_mat_elem(psi_det(1,1,i),N_int)
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, barycentric_electronic_energy, (N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! TODO : ASCII Elephant
+ END_DOC
+ integer :: istate,i
+
+ barycentric_electronic_energy(:) = 0.d0
+
+ do istate=1,N_states
+   do i=1,N_det
+     barycentric_electronic_energy(istate) += psi_coef(i,istate)*psi_coef(i,istate)*diagonal_H_matrix_on_psi_det(i)
+   enddo
+ enddo
+
+END_PROVIDER
+

src/ZMQ/utils.irp.f

@@ -785,13 +785,14 @@ integer function zmq_abort(zmq_to_qp_run_socket)
   sze = len(trim(message))
   rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
   if (rc /= sze) then
+    print *,  'zmq_abort: rc /= sze', rc, sze
     zmq_abort = -1
     return
   endif
   
   rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
   if (trim(message(1:rc)) /= 'ok') then
-    print *,  'zmq_abort: '//trim(message(1:rc))
+    print *,  'zmq_abort: ', rc, ':', trim(message(1:rc))
     zmq_abort = -1
     return
   endif