Merge branch 'QuantumPackage:dev-stable' into dev-stable

.github/workflows/compilation.yml

@@ -48,6 +48,7 @@ jobs:
               ./configure -i docopt || :
               ./configure -i resultsFile || :
               ./configure -i bats || :
+              ./configure -i trexio-nohdf5 || :
               ./configure -c ./config/gfortran_debug.cfg
        - name: Compilation
          run: |

.github/workflows/configuration.yml

@@ -22,7 +22,7 @@ jobs:
        - uses: actions/checkout@v3
        - name: Install dependencies
          run: |
-           sudo apt install gfortran gcc liblapack-dev libblas-dev wget python3 make m4 pkg-config
+           sudo apt install gfortran gcc liblapack-dev libblas-dev wget python3 make m4 pkg-config hdf5
        - name: zlib
          run: |
               ./configure -i zlib || echo OK
@@ -50,6 +50,12 @@ jobs:
        - name: bats
          run: |
               ./configure -i bats || echo OK
+       - name: trexio-nohdf5
+         run: |
+              ./configure -i trexio-nohdf5 || echo OK
+       - name: trexio
+         run: |
+              ./configure -i trexio || echo OK
        - name: Final check
          run: |
               ./configure -c config/gfortran_debug.cfg

.travis.yml

@@ -1,52 +0,0 @@
-#sudo: true
-#before_script:
-#  - sudo apt-get update -q
-#  - sudo apt-get remove curl
-#  - sudo apt-get remove zlib1g-dev
-#  - sudo apt-get install autoconf
-#  - sudo rm /usr/local/bin/bats
-
-os: linux
-
-dist: bionic
-
-sudo: false
-
-compiler: gfortran
-
-addons:
-  apt:
-    packages:
-    - gfortran
-    - gcc
-    - libatlas-base-dev
-#    - liblapack-dev
-#    - libblas-dev
-    - wget
-
-env:
-  - OPAMROOT=$HOME/.opam
-
-cache:
-  directories:
-  - $HOME/.opam/
-  - $HOME/cache
-
-language: python
-python:
-    - "3.7"
-
-stages:
-  - configuration
-  - compilation
-  - testing
-
-jobs:
-  include:
-    - stage: configuration
-      script: travis/configuration.sh
-    - stage: compilation
-      script: travis/compilation.sh
-    - stage: testing
-      script: travis/testing.sh
-

RELEASE_NOTES.org

@@ -9,15 +9,23 @@
   - Configure adapted for ARM
   - Added many types of integrals
   - Accelerated four-index transformation
-
-*** TODO: take from dev
-  - [ ] Added GTOs with complex exponent
-  - Updated version of f77-zmq
   - Added transcorrelated SCF
   - Added transcorrelated CIPSI
+  - Added CCSD and CCSD(T)
+  - Added MO localization
+  - Changed coupling parameters for ROHF
+  - General Davidson algorithm
+  - Accelerated restore_symmetry
+  - Point charges in the Hamiltonian
+  - Removed cryptokit dependency in OCaml
+  - Using now standard convention in RDM
+  - Added molecular properties
+  - [ ] Added GTOs with complex exponent
+
+*** TODO: take from dev
+  - Updated version of f77-zmq
   - Started to introduce shells in AOs
   - Added ECMD UEG functional
-  - General Davidson algorithm
 
 * Version 2.2
 

bin/qp_test

@@ -46,7 +46,7 @@ def main(arguments):
                 append_bats(dirname, filenames)
     else:
         for (dirname, _, filenames) in os.walk(os.getcwd(), followlinks=False):
-            if "IRPF90_temp" not in dirname:
+            if "IRPF90_temp" not in dirname and "external" not in dirname:
                 append_bats(dirname, filenames)
     l_bats = [y for _, y in sorted(l_bats)]
 
@@ -67,6 +67,7 @@ def main(arguments):
             os.system(test+" python3 bats_to_sh.py "+bats_file+
                                     "| bash")
         else:
+#           print(" ".join(["bats", "--verbose-run", "--trace", bats_file]))
             subprocess.check_call(["bats", "--verbose-run", "--trace", bats_file], env=os.environ)
 
 

config/ifort_2021_xHost.cfg

@@ -6,7 +6,7 @@
 # --align=32                 : Align all provided arrays on a 32-byte boundary
 #
 [COMMON]
-FC           : ifort -fpic
+FC           : ifort -fpic -diag-disable 5462
 LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=64 -DINTEL 

configure

@@ -9,6 +9,8 @@ echo "QP_ROOT="$QP_ROOT
 unset CC
 unset CCXX
 
+TREXIO_VERSION=2.3.2
+
 # Force GCC instead of ICC for dependencies
 export CC=gcc
 
@@ -189,7 +191,7 @@ if [[ "${PACKAGES}.x" != ".x" ]] ; then
 fi
 
 if [[ ${PACKAGES} = all ]] ; then
-        PACKAGES="zlib ninja zeromq f77zmq gmp ocaml docopt resultsFile bats"
+        PACKAGES="zlib ninja zeromq f77zmq gmp ocaml docopt resultsFile bats trexio"
 fi
 
 
@@ -203,6 +205,33 @@ for PACKAGE in ${PACKAGES} ; do
               mv ninja "\${QP_ROOT}"/bin/
 EOF
 
+    elif [[ ${PACKAGE} = trexio-nohdf5 ]] ; then
+
+            VERSION=$TREXIO_VERSION
+            execute << EOF
+              cd "\${QP_ROOT}"/external
+              wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
+              tar -zxf trexio-${VERSION}.tar.gz
+              cd trexio-${VERSION}
+              ./configure --prefix=\${QP_ROOT} --without-hdf5
+              make -j 8 && make -j 8 check && make -j 8 install
+              cp ${QP_ROOT}/include/trexio_f.f90 ${QP_ROOT}/src/ezfio_files
+              tar -zxvf "\${QP_ROOT}"/external/qp2-dependencies/${ARCHITECTURE}/ninja.tar.gz 
+              mv ninja "\${QP_ROOT}"/bin/
+EOF
+    elif [[ ${PACKAGE} = trexio ]] ; then
+
+            VERSION=$TREXIO_VERSION
+            execute << EOF
+              cd "\${QP_ROOT}"/external
+              wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
+              tar -zxf trexio-${VERSION}.tar.gz
+              cd trexio-${VERSION}
+              ./configure --prefix=\${QP_ROOT}
+              make -j 8 && make -j 8 check && make -j 8 install
+              cp ${QP_ROOT}/include/trexio_f.f90 ${QP_ROOT}/src/ezfio_files
+EOF
+ 
 
     elif [[ ${PACKAGE} = gmp ]] ; then
 
@@ -232,7 +261,7 @@ EOF
 
             execute << EOF
               cd "\${QP_ROOT}"/external
-              tar --gunzip --extract --file qp2-dependencies/f77-zmq-4.3.2.tar.gz
+              tar --gunzip --extract --file qp2-dependencies/f77-zmq-4.3.?.tar.gz
               cd f77-zmq-*
               ./configure --prefix=\$QP_ROOT
               export ZMQ_H="\$QP_ROOT"/include/zmq.h
@@ -338,6 +367,12 @@ if [[ ${ZEROMQ} = $(not_found) ]] ; then
         fail
 fi
 
+TREXIO=$(find_lib -ltrexio)
+if [[ ${TREXIO} = $(not_found) ]] ; then
+        error "TREXIO (trexio,trexio-nohdf5) is not installed. If you don't have HDF5, use trexio-nohdf5"
+        fail
+fi
+
 F77ZMQ=$(find_lib -lzmq -lf77zmq -lpthread)
 if [[ ${F77ZMQ} = $(not_found) ]] ; then
         error "Fortran binding of ZeroMQ (f77zmq) is not installed."

data/basis/none

@@ -0,0 +1,5 @@
+$DATA
+
+HYDROGEN
+
+$END

etc/qp.rc

@@ -110,6 +110,11 @@ function qp()
       unset COMMAND
       ;;
 
+    "test")
+      shift
+      qp_test $@
+      ;;
+
     *)
       which "qp_$1" &> /dev/null
       if [[ $? -eq 0 ]] ; then
@@ -183,7 +188,7 @@ _qp_Complete()
             ;;
         esac;;
       set_file)
-        COMPREPLY=( $(compgen -W "$(for i in * ; do [[ -f ${i}/ezfio/.version ]] && echo $i ; done)" -- ${cur} ) )
+        COMPREPLY=( $(compgen -W "$(for i in */ $(find . -name ezfio | sed 's/ezfio$/.version/') ; do [[ -f $i ]] && echo ${i%/.version} ; done)" -- ${cur} ) )
         return 0
         ;;
       plugins)
@@ -215,10 +220,15 @@ _qp_Complete()
             return 0
             ;;
         esac;;
+      test)
+        COMPREPLY=( $(compgen -W "-v -a " -- $cur ) )
+        return 0
+        ;;
       *)
         COMPREPLY=( $(compgen -W 'plugins set_file \
                                    unset_file man \
                                    create_ezfio \
+                                   test \
                                    convert_output_to_ezfio \
                                    -h update' -- $cur ) )
 

external/qp2-dependencies

@@ -1 +1 @@
-Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8
+Subproject commit e0d0e02e9f5ece138d1520106954a881ab0b8db2

ocaml/Input_ao_basis.ml

@@ -247,8 +247,7 @@ end = struct
 
 
   let read () =
-    if (Ezfio.has_ao_basis_ao_basis ()) then
+      try
-      begin
         let result =
           { ao_basis        = read_ao_basis ();
             ao_num          = read_ao_num () ;
@@ -267,9 +266,8 @@ end = struct
         |> MD5.to_string
         |> Ezfio.set_ao_basis_ao_md5 ;
         Some result
-      end
+      with
-    else
+      | _ -> (Ezfio.set_ao_basis_ao_md5 "None" ; None)
-      None
   ;;
 
 

ocaml/qp_create_ezfio.ml

@@ -478,6 +478,7 @@ let run ?o b au c d m p cart xyz_file =
         let nmax =
           Nucl_number.get_max ()
         in
+
         let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function
         | [] -> accu
         | e::tail ->
@@ -520,6 +521,8 @@ let run ?o b au c d m p cart xyz_file =
       in
       let long_basis = Long_basis.of_basis basis in
       let ao_num = List.length long_basis in
+      if ao_num > 0 then
+        begin
           Ezfio.set_ao_basis_ao_num ao_num;
           Ezfio.set_ao_basis_ao_basis b;
           Ezfio.set_basis_basis b;
@@ -655,6 +658,7 @@ let run ?o b au c d m p cart xyz_file =
           match Input.Ao_basis.read () with
           | None -> failwith "Error in basis"
           | Some x -> Input.Ao_basis.write x
+     end
   in
   let () =
     try write_file () with
@@ -781,7 +785,7 @@ If a file with the same name as the basis set exists, this file will be read.  O
   run ?o:output basis au charge dummy multiplicity pseudo cart xyz_filename
   )
   with
-  | Failure txt  -> Printf.eprintf "Fatal error: %s\n%!" txt
+(*  | Failure txt  -> Printf.eprintf "Fatal error: %s\n%!" txt *)
   | Command_line.Error txt  -> Printf.eprintf "Command line error: %s\n%!" txt
 
 

scripts/Hn.py

@@ -0,0 +1,43 @@
+#!/usr/bin/env python                                                                                                     
+import sys
+from math import *
+arg = sys.argv
+#f = open('data_dft','r')
+n = int(sys.argv[1])
+r = float(sys.argv[2])
+f = open('H'+str(n)+'_'+str(r)+'.xyz','w')
+string=str(n)+"\n"
+f.write(string)
+string="\n"
+f.write(string)
+for i in range(n):
+ x = r * cos(2.* i* pi/n)
+ y = r * sin(2.* i* pi/n)
+ z = 0.
+ string="H   "+str(x)+"  "+str(y)+"  "+str(z)+"\n"
+ f.write(string)
+ 
+#lines = f.readlines()
+#cipsi_dft= []
+#
+#dissoc = []
+#dissoc.append(float(-76.0179223470363))
+#dissoc.append(float(-76.0592367866993))
+#dissoc.append(float(-76.0678739715659))
+#delta_e = []
+#
+#for line in lines:
+#  data = line.split()
+#  if(len(data)>0):
+#   dft=float(data[1])
+#   fci=float(data[2])
+#   e=fci+dft
+#   cipsi_dft.append(e)
+#
+#print(*cipsi_dft,sep="  &  ")
+#
+#for i in 0,1,2:
+# delta_e.append(1000.*(dissoc[i] - cipsi_dft[i]))
+#
+#print(*delta_e,sep="  &  ")
+#

scripts/compilation/qp_create_ninja

@@ -25,7 +25,7 @@ except ImportError:
                                       "quantum_package.rc"))
 
     print("\n".join(["", "Error:", "source %s" % f, ""]))
-    sys.exit(1)
+    raise
 
 # Compress path
 def comp_path(path):
@@ -38,7 +38,7 @@ def comp_path(path):
 
 from qp_path import QP_ROOT, QP_SRC, QP_EZFIO
 
-LIB = " -lz"
+LIB = " -lz -ltrexio"
 EZFIO_LIB = join("$QP_ROOT", "lib", "libezfio_irp.a")
 ZMQ_LIB = join("$QP_ROOT", "lib", "libf77zmq.a") + " "  + join("$QP_ROOT", "lib", "libzmq.a") + " -lstdc++ -lrt -ldl"
 ROOT_BUILD_NINJA = join("$QP_ROOT", "config", "build.ninja")

scripts/ezfio_interface/ezfio_generate_provider.py

@@ -52,7 +52,7 @@ BEGIN_PROVIDER [ %(type)s, %(name)s %(size)s ]
     %(test_null_size)s
     call ezfio_has_%(ezfio_dir)s_%(ezfio_name)s(has)
     if (has) then
-      write(6,'(A)') '.. >>>>> [ IO READ: %(name)s ] <<<<< ..'
+!      write(6,'(A)') '.. >>>>> [ IO READ: %(name)s ] <<<<< ..'
       call ezfio_get_%(ezfio_dir)s_%(ezfio_name)s(%(name)s)
     else
       print *, '%(ezfio_dir)s/%(ezfio_name)s not found in EZFIO file'
@@ -117,7 +117,7 @@ END_PROVIDER
         output = self.output
         name = self.name
         l_write = ["",
-                   "  call write_time(%(output)s)",
+                   "!  call write_time(%(output)s)",
                    ""]
 
         self.write = "\n".join(l_write) % locals()
@@ -129,7 +129,7 @@ END_PROVIDER
                 write = self.write_correspondance[self.type]
 
                 l_write = ["",
-                          "  call write_time(%(output)s)",
+                          "!  call write_time(%(output)s)",
                           "  call %(write)s(%(output)s, %(name)s, &",
                           "     '%(name)s')",
                           ""]

scripts/ezfio_interface/qp_edit_template

@@ -172,11 +172,8 @@ let run check_only ?ndet ?state ezfio_filename =
 
   (* Reorder basis set *)
   begin
-    let aos =
     match Input.Ao_basis.read() with
-      | Some x -> x
+    | Some aos -> 
-      | _ -> assert false
-    in
         let ordering = Input.Ao_basis.ordering aos in
         let test = Array.copy ordering in
         Array.sort compare test ;
@@ -191,6 +188,7 @@ let run check_only ?ndet ?state ezfio_filename =
               let new_mos = Input.Mo_basis.reorder mos ordering in
               Input.Mo_basis.write new_mos
         end
+    | _ -> ()
   end;
 
   begin

scripts/get_fci_conv.sh

@@ -0,0 +1,7 @@
+file=$1
+grep "N_det             =" $1 | cut -d "=" -f 2 > N_det_tmp
+grep "E               =" $file | cut -d "=" -f 2 > E_tmp
+grep "E+PT2           =" $file | cut -d "=" -f 2 | cut -d "+" -f 1 > E+PT2_tmp
+grep "E+rPT2          =" $file | cut -d "=" -f 2 | cut -d "+" -f 1 > E+rPT2_tmp
+paste N_det_tmp E_tmp E+PT2_tmp E+rPT2_tmp | column -s ' ' -t > $file.conv_fci
+rm N_det_tmp E_tmp E+PT2_tmp E+rPT2_tmp

scripts/get_fci_tc_conv.sh

@@ -0,0 +1,2 @@
+file=$1
+grep "Ndet,E,E+PT2,E+RPT2,|PT2|=" $file | cut -d "=" -f 2 > ${file}.conv_fci_tc

scripts/qp_exc_energy.py

@@ -0,0 +1,188 @@
+#!/usr/bin/env python
+# Computes the error on the excitation energy of a CIPSI run.
+
+def student(p,df):
+    import scipy
+    from scipy.stats import t
+    return t.ppf(p, df)
+
+
+def chi2cdf(x,k):
+    import scipy
+    import scipy.stats
+    return scipy.stats.chi2.cdf(x,k)
+
+
+def jarque_bera(data):
+
+    n    = max(len(data), 2)
+    norm = 1./ sum( [ w for (_,w) in data ] )
+
+    mu     =   sum( [ w* x        for (x,w) in data ] ) * norm
+    sigma2 =   sum( [ w*(x-mu)**2 for (x,w) in data ] ) * norm
+    if sigma2 > 0.:
+      S      = ( sum( [ w*(x-mu)**3 for (x,w) in data ] ) * norm ) / sigma2**(3./2.)
+      K      = ( sum( [ w*(x-mu)**4 for (x,w) in data ] ) * norm ) / sigma2**2
+    else:
+      S = 0.
+      K = 0.
+
+    # Value of the Jarque-Bera test
+    JB = n/6. * (S**2 + 1./4. * (K-3.)**2)
+
+    # Probability that the data comes from a Gaussian distribution
+    p = 1. - chi2cdf(JB,2)
+
+    return JB, mu, sqrt(sigma2/(n-1)), p
+
+
+
+to_eV = 27.2107362681
+import sys, os
+import scipy
+import scipy.stats
+from math import sqrt, gamma, exp
+import qp_json
+
+
+def read_data(ezfio_filename,state):
+  """ Read energies and PT2 from input file """
+  data = qp_json.load_last(ezfio_filename)
+  for method in data.keys():
+    x = data[method]
+  lines = x
+
+  print(f"State: {state}")
+
+  gs = []
+  es = []
+
+  for l in lines:
+      try:
+        pt2_0 = l['states'][0]['pt2']
+        e_0   = l['states'][0]['energy']
+        pt2_1 = l['states'][state]['pt2']
+        e_1   = l['states'][state]['energy']
+        gs.append( (e_0, pt2_0) )
+        es.append( (e_1, pt2_1) )
+      except: pass
+
+  def f(p_1, p0, p1):
+      e, pt2 = p0
+      y0, x0 = p_1
+      y1, x1 = p1
+      try:
+        alpha = (y1-y0)/(x0-x1)
+      except ZeroDivisionError:
+        alpha = 1.
+      return [e, pt2, alpha]
+
+  for l in (gs, es):
+      p_1, p0, p1 = l[0], l[0], l[1]
+      l[0] = f(p_1, p0, p1)
+
+      for i in range(1,len(l)-1):
+          p_1 = (l[i-1][0], l[i-1][1])
+          p0  = l[i]
+          p1  = l[i+1]
+          l[i] = f(p_1, p0, p1)
+
+      i = len(l)-1
+      p_1 = (l[i-1][0], l[i-1][1])
+      p0  = l[i]
+      p1  = l[-1]
+      l[i] = f(p_1, p0, p1)
+
+  return [ x+y for x,y in zip(gs,es) ]
+
+
+def compute(data):
+
+    d = []
+    for e0, p0, a0, e1, p1, a1 in data:
+        x   = (e1+p1)-(e0+p0)
+        w   = 1./sqrt(p0**2 + p1**2)
+        bias = (a1-1.)*p1 - (a0-1.)*p0
+        d.append( (x,w,bias) )
+
+    x = []
+    target = (scipy.stats.norm.cdf(1.)-0.5)*2   # = 0.6827
+
+    print("| %2s | %8s | %8s | %8s | %8s | %8s |"%( "N", "DE", "+/-", "bias", "P(G)", "J"))
+    print("|----+----------+----------+----------+----------+----------|")
+    xmax = (0.,0.,0.,0.,0.,0,0.)
+    for i in range(len(data)-1):
+        jb, mu, sigma, p = jarque_bera( [ (x,w) for (x,w,bias) in d[i:] ] )
+        bias = sum ( [ w * e for (_,w,e) in d[i:] ] ) / sum ( [ w for (_,w,_) in d[i:] ] )
+        mu = (mu+0.5*bias) * to_eV
+        sigma = sigma * to_eV
+        bias = bias * to_eV
+        n = len(data[i:])
+        beta = student(0.5*(1.+target/p) ,n)
+        err = sigma * beta + 0.5*abs(bias)
+        print("| %2d | %8.3f | %8.3f | %8.3f | %8.3f | %8.3f |"%( n, mu, err, bias, p, jb))
+        if n < 3 :
+            continue
+        y = (err, p, mu, err, jb,n,bias)
+        if p > xmax[1]: xmax = y
+        if p < 0.8:
+            continue
+        x.append(y)
+
+    x = sorted(x)
+
+    print("|----+----------+----------+----------+----------+----------|")
+    if x != []:
+      xmax = x[0]
+    _, p, mu, err, jb, n, bias = xmax
+    beta = student(0.5*(1.+target/p),n)
+    print("| %2d | %8.3f | %8.3f | %8.3f | %8.3f | %8.3f |\n"%(n, mu, err, bias, p, jb))
+
+    return mu, err, bias, p
+
+ezfio_filename = sys.argv[1]
+print(ezfio_filename)
+if len(sys.argv) > 2:
+    state = int(sys.argv[2])
+else:
+    state = 1
+data = read_data(ezfio_filename,state)
+mu, err, bias, _ = compute(data)
+print(" %s: %8.3f +/- %5.3f eV\n"%(ezfio_filename, mu, err))
+
+import numpy as np
+A = np.array( [ [ data[-1][1], 1. ],
+                [ data[-2][1], 1. ] ] )
+B = np.array( [ [ data[-1][0] ],
+                [ data[-2][0] ] ] )
+E0 = np.linalg.solve(A,B)[1]
+A = np.array( [ [ data[-1][4], 1. ],
+                [ data[-2][4], 1. ] ] )
+B = np.array( [ [ data[-1][3] ],
+                [ data[-2][3] ] ] )
+E1 = np.linalg.solve(A,B)[1]
+average_2 = (E1-E0)*to_eV
+
+A = np.array( [ [ data[-1][1], 1. ],
+                [ data[-2][1], 1. ],
+                [ data[-3][1], 1. ] ] )
+B = np.array( [ [ data[-1][0] ],
+                [ data[-2][0] ],
+                [ data[-3][0] ] ] )
+E0 = np.linalg.lstsq(A,B,rcond=None)[0][1]
+A = np.array( [ [ data[-1][4], 1. ],
+                [ data[-2][4], 1. ],
+                [ data[-3][4], 1. ] ] )
+B = np.array( [ [ data[-1][3] ],
+                [ data[-2][3] ],
+                [ data[-3][3] ] ] )
+E1 = np.linalg.lstsq(A,B,rcond=None)[0][1]
+average_3 = (E1-E0)*to_eV
+
+exc = ((data[-1][3] + data[-1][4]) - (data[-1][0] + data[-1][1])) * to_eV
+error_2 = abs(average_2 - average_3)
+error_3 = abs(average_3 - exc)
+print(" 2-3 points: %.3f +/- %.3f "% (average_3, error_2))
+print(" largest wf: %.3f +/- %.3f  "%(average_3, error_3))
+
+

scripts/qp_extract_cipsi_data.py

@@ -0,0 +1,37 @@
+#!/usr/bin/env python3
+
+import qp_json
+import sys
+
+if len(sys.argv) == 1:
+  print(f"syntax: {sys.argv[0]} EZFIO_FILE")
+
+d = qp_json.load_all(sys.argv[1])
+
+k = [ x for x in d.keys() ]
+k.sort()
+
+print("# Energy  PT2   PT2_err  rPT2  rPT2_err  exFCI\n")
+for f in k:
+  try:
+    j = d[f]["fci"]
+  except:
+     continue
+
+  print(f"# {f}")
+  for e in j:
+
+       out = f" {e['n_det']:8d}"
+
+       nstates = len(e["states"])
+       for ee in e["states"]:
+         try:
+            exc_energy = ee['ex_energy'][0]
+         except:
+            exc_energy = 0.
+         out += f" {ee['energy']:16.8f}  {ee['pt2']:e}  {ee['pt2_err']:e}  {ee['rpt2']:e}  {ee['rpt2_err']:e}  {exc_energy:16.8f}"
+       print(out)
+
+  print("\n")
+
+

scripts/script_fci_tc.sh

@@ -0,0 +1,33 @@
+ source ~/qp2/quantum_package.rc  
+ alpha=1.8
+ input=O
+ basis=cc-pvdz
+ mult=3
+ output=${input}_${basis}_al_${alpha}
+ qp create_ezfio -b ${basis} ${input}.xyz -m $mult
+ qp run scf 
+ qp set perturbation pt2_max 0.0001
+ qp set_frozen_core 
+
+########## FCI CALCULATION FOR REFERENCE 
+ qp run fci | tee ${EZFIO_FILE}.fci.out
+ qp run sort_wf 
+ mv ${EZFIO_FILE}.wf_sorted ${EZFIO_FILE}_fci.wf_sorted
+########### TC SCF CALCULATION
+ qp reset -d 
+ qp set ao_two_e_erf_ints mu_erf 0.87 
+ qp set tc_keywords j1b_type 3
+ qp set tc_keywords j1b_pen "[${alpha}]"
+ qp set tc_keywords bi_ortho True 
+ qp set tc_keywords test_cycle_tc True
+ qp set tc_keywords write_tc_integ True
+ qp set tc_keywords read_tc_integ False
+ qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
+ qp set tc_keywords write_tc_integ False
+ qp set tc_keywords read_tc_integ True
+############ TC-FCI CALCULATION
+ qp run fci_tc_bi_ortho | tee ${EZFIO_FILE}.fci_tc_bi_ortho.out
+ grep "Ndet,E,E+PT2,E+RPT2,|PT2|=" ${EZFIO_FILE}.fci_tc_bi_ortho.out | cut -d "=" -f 2 > data_al_$alpha
+ qp run sort_wf 
+ mv ${EZFIO_FILE}.wf_sorted ${EZFIO_FILE}_tc_fci.wf_sorted
+ 

scripts/utility/qp_json.py

@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+import os
+import json
+
+def fix_json(s):
+   """Properly termitates an incomplete JSON file"""
+
+   s = s.replace(' ','')
+   s = s.replace('\n','')
+   s = s.replace('\t','')
+   s = s.replace(",{}",'')
+   tmp = [ c for c in s if c in "[]{}" ]
+   tmp = "".join(tmp)
+   tmp_old = ""
+   while tmp != tmp_old:
+      tmp_old = tmp
+      tmp = tmp.replace("{}","")
+      tmp = tmp.replace("[]","")
+   while s[-1] in [ ',', '\n', ' ', '\t' ]:
+      s = s[:-1]
+   tmp = [ c for c in tmp ]
+   tmp.reverse()
+   for c in tmp:
+      if c == '[': s += "]"
+      elif c == '{': s += "}"
+   return s
+
+
+def load(filename):
+  """Loads a JSON file after calling the fix_json function."""
+  with open(filename,'r') as f:
+    data = f.read()
+  new_data = fix_json(data)
+  return json.loads(new_data)
+
+
+def load_all(ezfio_filename):
+  """Loads all JSON files of an EZFIO."""
+  d = {}
+  prefix = ezfio_filename+'/json/'
+  for filename in [ x for x in os.listdir(prefix) if x.endswith(".json")]:
+    d[filename] = load(prefix+filename)
+  return d
+
+
+def load_last(ezfio_filename):
+  """Loads last JSON file of an EZFIO."""
+  d = {}
+  prefix = ezfio_filename+'/json/'
+  l = [ x for x in os.listdir(prefix) if x.endswith(".json")]
+  l.sort()
+  filename = l[-1]
+  print(filename)
+  return load(prefix+filename)
+
+
+def fix(ezfio_filename):
+  """Fixes all JSON files in an EZFIO."""
+  d = load_all(ezfio_filename)
+  prefix = ezfio_filename+'/json/'
+  for filename in d.keys():
+        with open(prefix+filename, 'w') as json_file:
+            json.dump(d[filename], json_file)
+
+
+

src/ao_many_one_e_ints/grad2_jmu_manu.irp.f

@@ -38,7 +38,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
      !$OMP          expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,                                             &
      !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,                                           &
      !$OMP          List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b,              &
-     !$OMP          ao_overlap_abs,sq_pi_3_2)
+     !$OMP          ao_overlap_abs,sq_pi_3_2,thrsh_cycle_tc)
  !$OMP DO SCHEDULE(dynamic)
   do ipoint = 1, n_points_final_grid
     r(1) = final_grid_points(1,ipoint)
@@ -46,7 +46,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
     r(3) = final_grid_points(3,ipoint)
     do i = 1, ao_num
       do j = i, ao_num
-        if(ao_overlap_abs(j,i) .lt. 1.d-12) then
+        if(ao_overlap_abs(j,i) .lt. thrsh_cycle_tc) then
           cycle
         endif
 
@@ -58,7 +58,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
         do i_fit = 1, ng_fit_jast
           expo_fit = expo_gauss_1_erf_x_2(i_fit)
           coef_fit = -0.25d0 *  coef_gauss_1_erf_x_2(i_fit)
-          if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle
+!          if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
           int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j)
           int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss 
         enddo
@@ -81,8 +81,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
             !DIR$ FORCEINLINE
             call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
             coef_fit = -0.25d0 *  coef_gauss_1_erf_x_2(i_fit) * coef
-!           if(dabs(coef_fit*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version
+!            if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
-            if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle
 !           call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, &
 !                 expo_fit, i, j, int_fit_v, n_points_final_grid)
             int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@@ -145,14 +144,14 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
      !$OMP          expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,      &
      !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,    &
      !$OMP          List_comb_thr_b3_cent, big_array,&
-     !$OMP          ao_abs_comb_b3_j1b,ao_overlap_abs)
+     !$OMP          ao_abs_comb_b3_j1b,ao_overlap_abs,thrsh_cycle_tc)
 !
  allocate(int_fit_v(n_points_final_grid))
  !$OMP DO SCHEDULE(dynamic)
  do i = 1, ao_num
    do j = i, ao_num
 
-     if(ao_overlap_abs(j,i) .lt. 1.d-12) then
+     if(ao_overlap_abs(j,i) .lt. thrsh_cycle_tc) then
        cycle
      endif
 
@@ -161,7 +160,6 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
          coef        = List_comb_thr_b3_coef  (i_1s,j,i)
          beta        = List_comb_thr_b3_expo  (i_1s,j,i)
          int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
-!         if(dabs(coef)*dabs(int_j1b).lt.1.d-15)cycle
          B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
          B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
          B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@@ -243,7 +241,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
  !$OMP          final_grid_points, ng_fit_jast,                     &
  !$OMP          expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2,           &
  !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,sq_pi_3_2,       & 
- !$OMP          List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b)
+ !$OMP          List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b,thrsh_cycle_tc)
  !$OMP DO
   do ipoint = 1, n_points_final_grid
     r(1) = final_grid_points(1,ipoint)
@@ -260,11 +258,11 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
         ! --- --- ---
 
         int_j1b = ao_abs_comb_b3_j1b(1,j,i)
-        if(dabs(int_j1b).lt.1.d-10) cycle
+        if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
         do i_fit = 1, ng_fit_jast
           expo_fit = expo_gauss_j_mu_x_2(i_fit)
           coef_fit = coef_gauss_j_mu_x_2(i_fit)
-          if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle
+!          if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
           int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
           tmp += coef_fit * int_fit
         enddo
@@ -278,7 +276,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
           coef        = List_comb_thr_b3_coef  (i_1s,j,i)
           beta        = List_comb_thr_b3_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+!         if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@@ -288,8 +286,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
             coef_fit = coef_gauss_j_mu_x_2(i_fit)
             !DIR$ FORCEINLINE
             call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
-!            if(dabs(coef_fit*coef*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version
+!            if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
-            if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle
             int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
             tmp += coef * coef_fit * int_fit
           enddo
@@ -350,7 +347,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
  !$OMP          final_grid_points, ng_fit_jast,                     &
  !$OMP          expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf,       &
  !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,       & 
- !$OMP          List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2)
+ !$OMP          List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2,thrsh_cycle_tc)
  !$OMP DO
 
   do ipoint = 1, n_points_final_grid
@@ -369,7 +366,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
           coef        = List_comb_thr_b3_coef  (i_1s,j,i)
           beta        = List_comb_thr_b3_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+          if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@@ -392,8 +389,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
             expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist 
             coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
             sq_alpha = alpha_1s_inv * dsqrt(alpha_1s_inv)
-!            if(dabs(coef_tmp*int_j1b) .lt. 1d-10) cycle ! old version
+!            if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. thrsh_cycle_tc) cycle
-            if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. 1d-10) cycle
             
             call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
 
@@ -470,13 +466,13 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
  !$OMP          expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf,       &
  !$OMP          ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2,   &
  !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,  ao_abs_comb_b3_j1b,     &
- !$OMP          List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test)
+ !$OMP          List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test,thrsh_cycle_tc)
  !$OMP DO
   do ipoint = 1, n_points_final_grid
     do i = 1, ao_num
       do j = i, ao_num
 
-        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10) cycle
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc) cycle
 
         r(1) = final_grid_points(1,ipoint)
         r(2) = final_grid_points(2,ipoint)
@@ -489,10 +485,10 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
         ! --- --- ---
 
         int_j1b = ao_abs_comb_b3_j1b(1,j,i)
-        if(dabs(int_j1b).lt.1.d-10) cycle
+!        if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
         do i_fit = 1, ng_fit_jast
           expo_fit = expo_gauss_j_mu_1_erf(i_fit)
-          if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.1.d-15) cycle
+!          if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.thrsh_cycle_tc) cycle
           coef_fit = coef_gauss_j_mu_1_erf(i_fit)
           int_fit  = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r,  1.d+9, r)
           tmp += coef_fit * int_fit
@@ -507,7 +503,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
           coef        = List_comb_thr_b3_coef  (i_1s,j,i)
           beta        = List_comb_thr_b3_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+!          if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@@ -517,7 +513,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
           do i_fit = 1, ng_fit_jast
             expo_fit = expo_gauss_j_mu_1_erf(i_fit)
             call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
-            if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.1.d-15)cycle
+!            if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.thrsh_cycle_tc)cycle
             coef_fit = coef_gauss_j_mu_1_erf(i_fit)
 
             alpha_1s     = beta + expo_fit
@@ -527,9 +523,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
             centr_1s(3)  = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
 
             expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
-            if(expo_coef_1s .gt. 20.d0) cycle
+!            if(expo_coef_1s .gt. 20.d0) cycle
             coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
-            if(dabs(coef_tmp) .lt. 1d-08) cycle
+!            if(dabs(coef_tmp) .lt. 1d-08) cycle
 
             int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s,  1.d+9, r)
 

src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f

@@ -31,7 +31,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
  !$OMP SHARED  (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points, &
  !$OMP          List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_j1b,  &
  !$OMP          v_ij_erf_rk_cst_mu_j1b_test, mu_erf,                                   &
- !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
+ !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
  !$OMP DO
   !do ipoint = 1, 10
   do ipoint = 1, n_points_final_grid
@@ -41,7 +41,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
 
     do i = 1, ao_num
       do j = i, ao_num
-        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
 
         tmp = 0.d0
         do i_1s = 1, List_comb_thr_b2_size(j,i)
@@ -49,7 +49,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
           coef        = List_comb_thr_b2_coef  (i_1s,j,i)
           beta        = List_comb_thr_b2_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+!         if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@@ -110,7 +110,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
  !$OMP SHARED  (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points,&
  !$OMP          List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,  &
  !$OMP          x_v_ij_erf_rk_cst_mu_j1b_test, mu_erf,ao_abs_comb_b2_j1b,         &
- !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma)
+ !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,thrsh_cycle_tc)
 ! !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,expo_erfc_mu_gauss)
  !$OMP DO
   do ipoint = 1, n_points_final_grid
@@ -120,7 +120,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
 
     do i = 1, ao_num
       do j = i, ao_num
-        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
 
         tmp_x = 0.d0
         tmp_y = 0.d0
@@ -130,19 +130,11 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
           coef        = List_comb_thr_b2_coef  (i_1s,j,i)
           beta        = List_comb_thr_b2_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+ !        if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
 
-!          if(ao_prod_center(1,j,i).ne.10000.d0)then
-!           ! approximate 1 - erf(mu r12) by a gaussian * 10
-!           !DIR$ FORCEINLINE
-!           call gaussian_product(expo_erfc_mu_gauss,r,     &
-!                ao_prod_sigma(j,i),ao_prod_center(1,j,i),  & 
-!                factor_ij_1s,beta_ij,center_ij_1s)
-!           if(dabs(coef * factor_ij_1s*int_j1b*10.d0 * dsqpi_3_2 * beta_ij**(-1.5d0)).lt.1.d-10)cycle 
-!          endif
           call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints        )
           call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center,  1.d+9, r, ints_coulomb)
 
@@ -216,7 +208,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
  !$OMP          expo_gauss_j_mu_x, coef_gauss_j_mu_x,               &
  !$OMP          List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & 
  !$OMP          List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b,   &
- !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
+ !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
  !$OMP DO
   do ipoint = 1, n_points_final_grid
     r(1) = final_grid_points(1,ipoint)
@@ -225,7 +217,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
 
     do i = 1, ao_num
       do j = i, ao_num
-        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
 
         tmp = 0.d0
 
@@ -234,11 +226,11 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
         ! --- --- ---
 
         int_j1b = ao_abs_comb_b2_j1b(1,j,i)
-        if(dabs(int_j1b).lt.1.d-10) cycle
+ !      if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
         do i_fit = 1, ng_fit_jast
           expo_fit = expo_gauss_j_mu_x(i_fit)
           coef_fit = coef_gauss_j_mu_x(i_fit)
-          if(ao_overlap_abs_grid(j,i).lt.1.d-15) cycle
+ !        if(ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc) cycle
           int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
           tmp += coef_fit * int_fit
         enddo
@@ -251,7 +243,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
           coef        = List_comb_thr_b2_coef  (i_1s,j,i)
           beta        = List_comb_thr_b2_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+!          if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@@ -259,9 +251,9 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
             expo_fit = expo_gauss_j_mu_x(i_fit)
             coef_fit = coef_gauss_j_mu_x(i_fit)
             coeftot = coef * coef_fit
-            if(dabs(coeftot).lt.1.d-15)cycle
+!            if(dabs(coeftot).lt.thrsh_cycle_tc)cycle
             call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u)
-            if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle
+!            if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc)cycle
             int_fit  = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
             tmp += coef * coef_fit * int_fit
           enddo
@@ -325,7 +317,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
  !$OMP          expo_gauss_j_mu_x, coef_gauss_j_mu_x,                                  &
  !$OMP          List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size,    & 
  !$OMP          List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, &
- !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
+ !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
  !$OMP DO
   do ipoint = 1, n_points_final_grid
     r(1) = final_grid_points(1,ipoint)
@@ -334,7 +326,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
 
     do i = 1, ao_num
       do j = i, ao_num
-        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
 
         tmp = 0.d0
 
@@ -343,7 +335,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
         ! --- --- ---
 
         int_j1b = ao_abs_comb_b2_j1b(1,j,i)
-        if(dabs(int_j1b).lt.1.d-10) cycle
+!        if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
         expo_fit = expo_good_j_mu_1gauss
         int_fit  = overlap_gauss_r12_ao(r, expo_fit, i, j)
         tmp += int_fit
@@ -356,7 +348,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
           coef        = List_comb_thr_b2_coef  (i_1s,j,i)
           beta        = List_comb_thr_b2_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+!          if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@@ -364,9 +356,9 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
             expo_fit = expo_good_j_mu_1gauss
             coef_fit = 1.d0
             coeftot = coef * coef_fit
-            if(dabs(coeftot).lt.1.d-15)cycle
+            if(dabs(coeftot).lt.thrsh_cycle_tc)cycle
             call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u)
-            if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle
+            if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc)cycle
             int_fit  = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
             tmp += coef * coef_fit * int_fit
 !          enddo

src/ao_many_one_e_ints/listj1b_sorted.irp.f

@@ -3,15 +3,16 @@
 &BEGIN_PROVIDER [ integer, max_List_comb_thr_b2_size]
  implicit none
  integer :: i_1s,i,j,ipoint
- double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: coef,beta,center(3),int_j1b
  double precision :: r(3),weight,dist
- thr = 1.d-15
  List_comb_thr_b2_size = 0
+ print*,'List_all_comb_b2_size = ',List_all_comb_b2_size
+! pause
  do i = 1, ao_num
   do j = i, ao_num
    do i_1s = 1, List_all_comb_b2_size
      coef        = List_all_comb_b2_coef  (i_1s)
-     if(dabs(coef).lt.1.d-15)cycle
+     if(dabs(coef).lt.thrsh_cycle_tc)cycle
      beta        = List_all_comb_b2_expo  (i_1s)
      beta = max(beta,1.d-12)
      center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
@@ -24,7 +25,7 @@
       dist += ( center(3) - r(3) )*( center(3) - r(3) )
       int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
      enddo
-     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+     if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
       List_comb_thr_b2_size(j,i) += 1
      endif
    enddo
@@ -40,6 +41,7 @@
   list(i) = maxval(List_comb_thr_b2_size(:,i))
  enddo
  max_List_comb_thr_b2_size = maxval(list) 
+ print*,'max_List_comb_thr_b2_size = ',max_List_comb_thr_b2_size
  
 END_PROVIDER 
 
@@ -49,16 +51,15 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_j1b, ( max_List_comb_thr_b2_size ,ao_num, ao_num)]
  implicit none
  integer :: i_1s,i,j,ipoint,icount
- double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: coef,beta,center(3),int_j1b
  double precision :: r(3),weight,dist
- thr = 1.d-15
  ao_abs_comb_b2_j1b = 10000000.d0
  do i = 1, ao_num
   do j = i, ao_num
    icount = 0
    do i_1s = 1, List_all_comb_b2_size
      coef        = List_all_comb_b2_coef  (i_1s)
-     if(dabs(coef).lt.1.d-12)cycle
+     if(dabs(coef).lt.thrsh_cycle_tc)cycle
      beta        = List_all_comb_b2_expo  (i_1s)
      center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
      int_j1b = 0.d0
@@ -70,7 +71,7 @@ END_PROVIDER
       dist += ( center(3) - r(3) )*( center(3) - r(3) )
       int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
      enddo
-     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+     if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
       icount += 1
       List_comb_thr_b2_coef(icount,j,i) = coef
       List_comb_thr_b2_expo(icount,j,i) = beta
@@ -98,17 +99,17 @@ END_PROVIDER
 &BEGIN_PROVIDER [ integer, max_List_comb_thr_b3_size]
  implicit none
  integer :: i_1s,i,j,ipoint
- double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: coef,beta,center(3),int_j1b
  double precision :: r(3),weight,dist
- thr = 1.d-15
  List_comb_thr_b3_size = 0
+ print*,'List_all_comb_b3_size = ',List_all_comb_b3_size
  do i = 1, ao_num
   do j = 1, ao_num
    do i_1s = 1, List_all_comb_b3_size
      coef        = List_all_comb_b3_coef  (i_1s)
      beta        = List_all_comb_b3_expo  (i_1s)
      center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
-     if(dabs(coef).lt.thr)cycle
+     if(dabs(coef).lt.thrsh_cycle_tc)cycle
      int_j1b = 0.d0
      do ipoint = 1, n_points_extra_final_grid
       r(1:3) = final_grid_points_extra(1:3,ipoint)
@@ -118,7 +119,7 @@ END_PROVIDER
       dist += ( center(3) - r(3) )*( center(3) - r(3) )
       int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
      enddo
-     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+     if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
       List_comb_thr_b3_size(j,i) += 1
      endif
    enddo
@@ -144,9 +145,8 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, ao_abs_comb_b3_j1b, ( max_List_comb_thr_b3_size ,ao_num, ao_num)]
  implicit none
  integer :: i_1s,i,j,ipoint,icount
- double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: coef,beta,center(3),int_j1b
  double precision :: r(3),weight,dist
- thr = 1.d-15
  ao_abs_comb_b3_j1b = 10000000.d0
  do i = 1, ao_num
   do j = 1, ao_num
@@ -156,7 +156,7 @@ END_PROVIDER
      beta        = List_all_comb_b3_expo  (i_1s)
      beta = max(beta,1.d-12)
      center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
-     if(dabs(coef).lt.thr)cycle
+     if(dabs(coef).lt.thrsh_cycle_tc)cycle
      int_j1b = 0.d0
      do ipoint = 1, n_points_extra_final_grid
       r(1:3) = final_grid_points_extra(1:3,ipoint)
@@ -166,7 +166,7 @@ END_PROVIDER
       dist += ( center(3) - r(3) )*( center(3) - r(3) )
       int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
      enddo
-     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+     if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
       icount += 1
       List_comb_thr_b3_coef(icount,j,i) = coef
       List_comb_thr_b3_expo(icount,j,i) = beta
@@ -177,15 +177,5 @@ END_PROVIDER
   enddo 
  enddo
 
-! do i = 1, ao_num
-!  do j = 1, i-1
-!    do icount = 1, List_comb_thr_b3_size(j,i)
-!     List_comb_thr_b3_coef(icount,j,i) = List_comb_thr_b3_coef(icount,i,j)
-!     List_comb_thr_b3_expo(icount,j,i) = List_comb_thr_b3_expo(icount,i,j)
-!     List_comb_thr_b3_cent(1:3,icount,j,i) = List_comb_thr_b3_cent(1:3,icount,i,j)
-!    enddo
-!  enddo
-! enddo
- 
 END_PROVIDER 
 

src/ao_two_e_ints/EZFIO.cfg

@@ -11,3 +11,8 @@ interface: ezfio,provider,ocaml
 default: False
 ezfio_name: direct
 
+[do_ao_cholesky]
+type: logical
+doc: Perform Cholesky decomposition of AO integrals
+interface: ezfio,provider,ocaml
+default: False

src/ao_two_e_ints/cholesky.irp.f

@@ -0,0 +1,100 @@
+BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ]
+ implicit none
+ BEGIN_DOC
+ ! Number of Cholesky vectors in AO basis
+ END_DOC
+
+ integer :: i,j,k,l
+ double precision :: xnorm0, x, integral
+ double precision, external :: ao_two_e_integral
+
+ cholesky_ao_num_guess = 0
+ xnorm0 = 0.d0
+ x = 0.d0
+ do j=1,ao_num
+   do i=1,ao_num
+     integral = ao_two_e_integral(i,i,j,j)
+     if (integral > ao_integrals_threshold) then
+       cholesky_ao_num_guess += 1
+     else
+       x += integral
+     endif
+   enddo
+ enddo
+ print *, 'Cholesky decomposition of AO integrals'
+ print *, '--------------------------------------'
+ print *, ''
+ print *, 'Estimated Error: ', x
+ print *, 'Guess size: ', cholesky_ao_num_guess, '(', 100.d0*dble(cholesky_ao_num_guess)/dble(ao_num*ao_num), ' %)'
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, cholesky_ao_num ]
+&BEGIN_PROVIDER [ double precision, cholesky_ao, (ao_num, ao_num, cholesky_ao_num_guess) ]
+ use mmap_module
+ implicit none
+ BEGIN_DOC
+ ! Cholesky vectors in AO basis: (ik|a):
+ ! <ij|kl> = (ik|jl) = sum_a (ik|a).(a|jl)
+ END_DOC
+
+ type(c_ptr) :: ptr
+ integer :: fd, i,j,k,l, rank
+ double precision, pointer :: ao_integrals(:,:,:,:)
+ double precision, external :: ao_two_e_integral
+
+ ! Store AO integrals in a memory mapped file
+ call mmap(trim(ezfio_work_dir)//'ao_integrals', &
+   (/ int(ao_num,8), int(ao_num,8), int(ao_num,8), int(ao_num,8) /), &
+   8, fd, .False., ptr)
+ call c_f_pointer(ptr, ao_integrals, (/ao_num, ao_num, ao_num, ao_num/))
+
+ double precision :: integral
+ logical, external :: ao_two_e_integral_zero
+ !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,k,l, integral) SCHEDULE(dynamic)
+ do l=1,ao_num
+  do j=1,l
+   do k=1,ao_num
+    do i=1,k
+     if (ao_two_e_integral_zero(i,j,k,l)) cycle
+     integral = ao_two_e_integral(i,k,j,l)
+     ao_integrals(i,k,j,l) = integral
+     ao_integrals(k,i,j,l) = integral
+     ao_integrals(i,k,l,j) = integral
+     ao_integrals(k,i,l,j) = integral
+    enddo
+   enddo
+  enddo
+ enddo
+ !$OMP END PARALLEL DO
+
+ ! Call Lapack
+ cholesky_ao_num = cholesky_ao_num_guess
+ call pivoted_cholesky(ao_integrals, cholesky_ao_num, ao_integrals_threshold, ao_num*ao_num, cholesky_ao)
+ print *, 'Rank: ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
+
+ ! Remove mmap
+ double precision, external :: getUnitAndOpen
+ call munmap( &
+   (/ int(ao_num,8), int(ao_num,8), int(ao_num,8), int(ao_num,8) /), &
+   8, fd, ptr)
+ open(unit=99,file=trim(ezfio_work_dir)//'ao_integrals')
+ close(99, status='delete')
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, cholesky_ao_transp, (cholesky_ao_num, ao_num, ao_num) ]
+ implicit none
+ BEGIN_DOC
+! Transposed of the Cholesky vectors in AO basis set
+ END_DOC
+ integer :: i,j,k
+ do j=1,ao_num
+  do i=1,ao_num
+   do k=1,ao_num
+    cholesky_ao_transp(k,i,j) = cholesky_ao(i,j,k)
+   enddo
+  enddo
+ enddo
+END_PROVIDER
+

src/ao_two_e_ints/map_integrals.irp.f

@@ -486,7 +486,7 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val)
   PROVIDE ao_two_e_integrals_in_map ao_integrals_map
 
   if (ao_one_e_integral_zero(j,l)) then
-    out_val = 0.d0
+    out_val(1:sze) = 0.d0
     return
   endif
 

src/bi_ort_ints/three_body_ijmk.irp.f

@@ -27,7 +27,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_direct_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_direct_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -74,7 +74,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_1_bi_ort, (mo_num, mo_num
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_cycle_1_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -121,7 +121,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_2_bi_ort, (mo_num, mo_num
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_cycle_2_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -168,7 +168,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch23_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_exch23_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -214,7 +214,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch13_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_exch13_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -261,7 +261,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch12_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_exch12_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num

src/bi_ort_ints/three_body_ijmkl.irp.f

@@ -26,7 +26,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_direct_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -75,7 +75,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -124,7 +124,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -173,7 +173,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -222,7 +222,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -271,7 +271,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_exch12_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num

src/bi_ort_ints/total_twoe_pot.irp.f

@@ -57,6 +57,7 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
 
     PROVIDE ao_tc_sym_two_e_pot_in_map
 
+    !!! TODO :: OPENMP
     do j = 1, ao_num
       do l = 1, ao_num
         do i = 1, ao_num

src/ccsd/80.ccsd_spin.bats

@@ -0,0 +1,225 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
+
+
+function run() {
+  thresh1=1e-6
+  thresh2=1e-6
+  test_exe scf || skip
+  qp set_file $1
+  qp edit --check
+  #qp run scf
+  qp set_frozen_core
+  qp set utils_cc cc_par_t true
+  qp set utils_cc cc_thresh_conv 1e-12
+  file="$(echo $1 | sed 's/.ezfio//g')"
+  qp run ccsd_spin_orb | tee $file.ccsd.out
+  energy1="$(grep 'E(CCSD)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
+  energy2="$(grep 'E(T)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
+  #rm $file.ccsd.out
+  eq $energy1 $2 $thresh1
+  eq $energy2 $3 $thresh2
+}
+
+@test "b2_stretched" {
+run  b2_stretched.ezfio -49.136487344382 -0.003497589175
+}
+
+@test "be" {
+run  be.ezfio -14.623559003577 -0.000230982022
+}
+
+@test "c2h2" {
+run  c2h2.ezfio -12.394008897618 -0.010790491561
+}
+
+@test "ch4" {
+run  ch4.ezfio -40.390721785799 -0.004476100282
+}
+
+@test "clf" {
+run  clf.ezfio -559.186562904081 -0.006577143392
+}
+
+@test "clo" {
+run  clo.ezfio -534.564874409332 -0.007584571424
+}
+
+@test "co2" {
+run  co2.ezfio -188.129602527766 -0.018040668885
+}
+
+@test "dhno" {
+run  dhno.ezfio -130.816650109473 -0.012197331453
+}
+
+@test "f2" {
+run  f2.ezfio -199.287826338097 -0.017592872692
+}
+
+@test "f" {
+run  f.ezfio -99.616644511121 -0.003624525307
+}
+
+@test "h2o2" {
+run  h2o2.ezfio -151.182552729963 -0.009511682086
+}
+
+@test "h2o" {
+run  h2o.ezfio -76.237710276526 -0.003001800577
+}
+
+@test "h2s" {
+run  h2s.ezfio -398.861214015390 -0.003300559757
+}
+
+@test "h3coh" {
+run  h3coh.ezfio -115.221296424969 -0.003566171432
+}
+
+@test "hbo" {
+run  hbo.ezfio -100.213539770415 -0.006851489212
+}
+
+@test "hcn" {
+run  hcn.ezfio -93.190247992657 -0.013418135043
+}
+
+@test "hco" {
+run  hco.ezfio -113.405413962350 -0.007973455337
+}
+
+@test "lif" {
+run  lif.ezfio -107.270402903250 -0.007742969005
+}
+
+@test "n2" {
+run  n2.ezfio -109.355358930472 -0.018477744342
+}
+
+@test "n2h4" {
+run  n2h4.ezfio -111.556885923139 -0.009048077008
+}
+
+@test "nh3" {
+run  nh3.ezfio -56.465503060954 -0.007638273755
+}
+
+@test "oh" {
+run  oh.ezfio -75.614606132774 -0.004126661739
+}
+
+@test "sih2_3b1" {
+run  sih2_3b1.ezfio -290.016780973072 -0.000497825874
+}
+
+@test "sih3" {
+run  sih3.ezfio -5.575343504534 -0.002094123268
+}
+
+@test "so" {
+run  so.ezfio -26.035945178665 -0.010594351274
+}
+
+#@test "b2_stretched" {
+#run  b2_stretched.ezfio -49.136487344382 -49.139984933557
+#}
+#
+#@test "be" {
+#run  be.ezfio -14.623559003577 -14.623789985599
+#}
+#
+#@test "c2h2" {
+#run  c2h2.ezfio -12.394008897618 -12.404799389179
+#}
+#
+#@test "ch4" {
+#run  ch4.ezfio -40.390721784961 -40.395197884406
+#}
+#
+#@test "clf" {
+#run  clf.ezfio -559.186562906072 -559.193140046904
+#}
+#
+#@test "clo" {
+#run  clo.ezfio -534.564874409333 -534.572458980757
+#}
+#
+#@test "co2" {
+#run  co2.ezfio -188.129602511724 -188.147643198675
+#}
+#
+#@test "dhno" {
+#run  dhno.ezfio -130.816650109473 -130.828847440925
+#}
+#
+#@test "f2" {
+#run  f2.ezfio -199.287826338097 -199.305419210789
+#}
+#
+#@test "f" {
+#run  f.ezfio -99.616644511120 -99.620269036428
+#}
+#
+#@test "h2o2" {
+#run  h2o2.ezfio -151.182552729963 -151.192064412049
+#}
+#
+#@test "h2o" {
+#run  h2o.ezfio -76.237710276526 -76.240712077103
+#}
+#
+#@test "h2s" {
+#run  h2s.ezfio -398.861214015416 -398.864514575146
+#}
+#
+#@test "h3coh" {
+#run  h3coh.ezfio -115.221296424969 -115.224862596401
+#}
+#
+#@test "hbo" {
+#run  hbo.ezfio -100.213539770415 -100.220391259627
+#}
+#
+#@test "hcn" {
+#run  hcn.ezfio -93.190247983000 -93.203666131216
+#}
+#
+#@test "hco" {
+#run  hco.ezfio -113.405413962350 -113.413387417687
+#}
+#
+#@test "lif" {
+#run  lif.ezfio -107.270402903211 -107.278145872216
+#}
+#
+#@test "n2" {
+#run  n2.ezfio -109.355358930472 -109.373836674814
+#}
+#
+#@test "n2h4" {
+#run  n2h4.ezfio -111.556885922642 -111.565934000556
+#}
+#
+#@test "nh3" {
+#run  nh3.ezfio -56.465503060954 -56.473141334709
+#}
+#
+#@test "oh" {
+#run  oh.ezfio -75.614606131897 -75.618732794235
+#}
+#
+#@test "sih2_3b1" {
+#run  sih2_3b1.ezfio -290.016780973071 -290.017278798946
+#}
+#
+#@test "sih3" {
+#run  sih3.ezfio -5.575343504534 -5.577437627802
+#}
+#
+#@test "so" {
+#run  so.ezfio -26.035945181998 -26.046539528491
+#}
+

src/ccsd/81.ccsd_space.bats

@@ -0,0 +1,225 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
+
+
+function run() {
+  thresh1=1e-6
+  thresh2=1e-6
+  test_exe scf || skip
+  qp set_file $1
+  qp edit --check
+  #qp run scf
+  qp set_frozen_core
+  qp set utils_cc cc_par_t true
+  qp set utils_cc cc_thresh_conv 1e-12
+  file="$(echo $1 | sed 's/.ezfio//g')"
+  qp run ccsd_space_orb | tee $file.ccsd.out
+  energy1="$(grep 'E(CCSD)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
+  energy2="$(grep 'E(T)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
+  #rm $file.ccsd.out
+  eq $energy1 $2 $thresh1
+  eq $energy2 $3 $thresh2
+}
+
+@test "b2_stretched" {
+run  b2_stretched.ezfio -49.136487344382 -0.003497589175
+}
+
+@test "be" {
+run  be.ezfio -14.623559003577 -0.000230982022
+}
+
+@test "c2h2" {
+run  c2h2.ezfio -12.394008897618 -0.010790491561
+}
+
+@test "ch4" {
+run  ch4.ezfio -40.390721785799 -0.004476100282
+}
+
+@test "clf" {
+run  clf.ezfio -559.186562904081 -0.006577143392
+}
+
+#@test "clo" {
+#run  clo.ezfio -534.564874409332 -0.007584571424
+#}
+
+@test "co2" {
+run  co2.ezfio -188.129602527766 -0.018040668885
+}
+
+#@test "dhno" {
+#run  dhno.ezfio -130.816650109473 -0.012197331453
+#}
+
+@test "f2" {
+run  f2.ezfio -199.287826338097 -0.017592872692
+}
+
+#@test "f" {
+#run  f.ezfio -99.616644511121 -0.003624525307
+#}
+
+@test "h2o2" {
+run  h2o2.ezfio -151.182552729963 -0.009511682086
+}
+
+@test "h2o" {
+run  h2o.ezfio -76.237710276526 -0.003001800577
+}
+
+@test "h2s" {
+run  h2s.ezfio -398.861214015390 -0.003300559757
+}
+
+@test "h3coh" {
+run  h3coh.ezfio -115.221296424969 -0.003566171432
+}
+
+@test "hbo" {
+run  hbo.ezfio -100.213539770415 -0.006851489212
+}
+
+@test "hcn" {
+run  hcn.ezfio -93.190247992657 -0.013418135043
+}
+
+#@test "hco" {
+#run  hco.ezfio -113.405413962350 -0.007973455337
+#}
+
+@test "lif" {
+run  lif.ezfio -107.270402903250 -0.007742969005
+}
+
+@test "n2" {
+run  n2.ezfio -109.355358930472 -0.018477744342
+}
+
+@test "n2h4" {
+run  n2h4.ezfio -111.556885923139 -0.009048077008
+}
+
+@test "nh3" {
+run  nh3.ezfio -56.465503060954 -0.007638273755
+}
+
+#@test "oh" {
+#run  oh.ezfio -75.614606132774 -0.004126661739
+#}
+
+#@test "sih2_3b1" {
+#run  sih2_3b1.ezfio -290.016780973072 -0.000497825874
+#}
+
+#@test "sih3" {
+#run  sih3.ezfio -5.575343504534 -0.002094123268
+#}
+
+#@test "so" {
+#run  so.ezfio -26.035945178665 -0.010594351274
+#}
+
+#@test "b2_stretched" {
+#run  b2_stretched.ezfio -49.136487344382 -49.139984933557
+#}
+#
+#@test "be" {
+#run  be.ezfio -14.623559003577 -14.623789985599
+#}
+#
+#@test "c2h2" {
+#run  c2h2.ezfio -12.394008897618 -12.404799389179
+#}
+#
+#@test "ch4" {
+#run  ch4.ezfio -40.390721784961 -40.395197884406
+#}
+#
+#@test "clf" {
+#run  clf.ezfio -559.186562906072 -559.193140046904
+#}
+#
+##@test "clo" {
+##run  clo.ezfio -534.564874409333 -534.572458980757
+##}
+#
+#@test "co2" {
+#run  co2.ezfio -188.129602511724 -188.147643198675
+#}
+#
+##@test "dhno" {
+##run  dhno.ezfio -130.816650109473 -130.828847440925
+##}
+#
+#@test "f2" {
+#run  f2.ezfio -199.287826338097 -199.305419210789
+#}
+#
+##@test "f" {
+##run  f.ezfio -99.616644511120 -99.620269036428
+##}
+#
+#@test "h2o2" {
+#run  h2o2.ezfio -151.182552729963 -151.192064412049
+#}
+#
+#@test "h2o" {
+#run  h2o.ezfio -76.237710276526 -76.240712077103
+#}
+#
+#@test "h2s" {
+#run  h2s.ezfio -398.861214015416 -398.864514575146
+#}
+#
+#@test "h3coh" {
+#run  h3coh.ezfio -115.221296424969 -115.224862596401
+#}
+#
+#@test "hbo" {
+#run  hbo.ezfio -100.213539770415 -100.220391259627
+#}
+#
+#@test "hcn" {
+#run  hcn.ezfio -93.190247983000 -93.203666131216
+#}
+#
+##@test "hco" {
+##run  hco.ezfio -113.405413962350 -113.413387417687
+##}
+#
+#@test "lif" {
+#run  lif.ezfio -107.270402903211 -107.278145872216
+#}
+#
+#@test "n2" {
+#run  n2.ezfio -109.355358930472 -109.373836674814
+#}
+#
+#@test "n2h4" {
+#run  n2h4.ezfio -111.556885922642 -111.565934000556
+#}
+#
+#@test "nh3" {
+#run  nh3.ezfio -56.465503060954 -56.473141334709
+#}
+#
+##@test "oh" {
+##run  oh.ezfio -75.614606131897 -75.618732794235
+##}
+#
+##@test "sih2_3b1" {
+##run  sih2_3b1.ezfio -290.016780973071 -290.017278798946
+##}
+#
+##@test "sih3" {
+##run  sih3.ezfio -5.575343504534 -5.577437627802
+##}
+#
+##@test "so" {
+##run  so.ezfio -26.035945181998 -26.046539528491
+##}
+

src/ccsd/NEED

@@ -0,0 +1,2 @@
+hartree_fock
+utils_cc

src/ccsd/README.md

@@ -0,0 +1,31 @@
+# CCSD in spin orbitals and spatial orbitals
+
+CCSD and CCSD(T) in spin orbitals for open and closed shell systems.
+CCSD and CCSD(T) in spatial orbitals for closed shell systems.
+
+## Calculations
+The program will automatically choose the version in spin or spatial orbitals
+To run the general program:  
+```
+qp run ccsd
+```
+Nevertheless, you can enforce the run in spin orbitals with
+```
+qp run ccsd_spin_orb
+```
+
+## Settings
+The settings can be changed with:
+```
+qp set utils_cc cc_#param #val
+```
+For more informations on the settings, look at the module utils_cc and its documentation. 
+
+## Org files
+The org files are stored in the directory org in order to avoid overwriting on user changes.
+The org files can be modified, to export the change to the source code, run
+```
+./TANGLE_org_mode.sh and
+mv *.irp.f ../.
+```
+

src/ccsd/ccsd.irp.f

@@ -0,0 +1,18 @@
+program ccsd
+
+  implicit none
+
+  BEGIN_DOC
+  ! CCSD program
+  END_DOC
+
+  read_wf = .True.
+  touch read_wf
+
+  if (.not. cc_ref_is_open_shell) then
+    call run_ccsd_space_orb
+  else 
+    call run_ccsd_spin_orb
+  endif
+
+end

src/ccsd/ccsd_space_orb.irp.f

@@ -0,0 +1,12 @@
+! Code
+
+program ccsd
+
+  implicit none
+
+  read_wf = .True.
+  touch read_wf
+
+  call run_ccsd_space_orb
+  
+end

src/ccsd/ccsd_spin_orb.irp.f

@@ -0,0 +1,16 @@
+! Prog
+
+program ccsd
+
+  implicit none
+
+  BEGIN_DOC
+  ! CCSD in spin orbitals
+  END_DOC
+
+  read_wf = .True.
+  touch read_wf
+
+  call run_ccsd_spin_orb
+  
+end

src/ccsd/ccsd_t_space_orb.irp.f

@@ -0,0 +1,513 @@
+! Dumb way
+
+subroutine ccsd_par_t_space(nO,nV,t1,t2,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO, nV)
+  double precision, intent(in)  :: t2(nO, nO, nV, nV)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: W(:,:,:,:,:,:)
+  double precision, allocatable :: V(:,:,:,:,:,:)
+  integer :: i,j,k,a,b,c
+
+  allocate(W(nO,nO,nO,nV,nV,nV))
+  allocate(V(nO,nO,nO,nV,nV,nV))
+
+  call form_w(nO,nV,t2,W)
+  call form_v(nO,nV,t1,W,V)
+
+  energy = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c))  !delta_ooovvv(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  energy = energy / 3d0
+
+  deallocate(V,W)
+end
+
+subroutine form_w(nO,nV,t2,W)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t2(nO, nO, nV, nV)
+  double precision, intent(out) :: W(nO, nO, nO, nV, nV, nV)
+
+  integer :: i,j,k,l,a,b,c,d
+
+  W = 0d0
+  do c = 1, nV
+    print*,'W:',c,'/',nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+
+              do d = 1, nV
+                W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+                ! chem (bd|ai)
+                ! phys <ba|di>
+                + cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+                + cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
+                + cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
+                + cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
+                + cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
+                + cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
+              enddo
+
+              do l = 1, nO
+                W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+                ! chem (ck|jl)
+                ! phys <cj|kl>
+                - cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
+                - cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
+                - cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
+                - cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
+                - cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
+                - cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
+              enddo
+
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end
+
+subroutine form_v(nO,nV,t1,w,v)
+
+implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO, nV)
+  double precision, intent(in)  :: W(nO, nO, nO, nV, nV, nV)
+  double precision, intent(out) :: V(nO, nO, nO, nV, nV, nV)
+
+  integer :: i,j,k,a,b,c
+
+  V = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+              + cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
+              + cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
+              + cc_space_v_vvoo(a,b,i,j) * t1(k,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end
+
+! Main
+
+subroutine ccsd_par_t_space_v2(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO,nV), f_o(nO), f_v(nV)
+  double precision, intent(in)  :: t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO), v_vooo(nV,nO,nO,nO)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: W(:,:,:,:,:,:)
+  double precision, allocatable :: V(:,:,:,:,:,:)
+  double precision, allocatable :: W_ijk(:,:,:), V_ijk(:,:,:)
+  double precision, allocatable :: X_vvvo(:,:,:,:), X_ovoo(:,:,:,:), X_vvoo(:,:,:,:)
+  double precision, allocatable :: T_vvoo(:,:,:,:), T_ovvo(:,:,:,:), T_vo(:,:)
+  integer                       :: i,j,k,l,a,b,c,d
+  double precision              :: e,ta,tb, delta, delta_ijk
+
+  !allocate(W(nV,nV,nV,nO,nO,nO))
+  !allocate(V(nV,nV,nV,nO,nO,nO))
+  allocate(W_ijk(nV,nV,nV), V_ijk(nV,nV,nV))
+  allocate(X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO), X_vvoo(nV,nV,nO,nO))
+  allocate(T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO), T_vo(nV,nO))
+
+  ! Temporary arrays
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,T_vvoo,T_ovvo,T_vo,X_vvvo,X_ovoo,X_vvoo, &
+  !$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
+  !$OMP PRIVATE(a,b,c,d,i,j,k,l) &
+  !$OMP DEFAULT(NONE)
+
+  !v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+  !X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j)
+
+  !$OMP DO collapse(3)
+  do i = 1, nO
+    do a = 1, nV
+      do b = 1, nV
+        do d = 1, nV
+          X_vvvo(d,b,a,i) = v_vvvo(b,a,d,i)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(3)
+  do j = 1, nO
+    do k = 1, nO
+      do c = 1, nV
+        do d = 1, nV
+          T_vvoo(d,c,k,j) = t2(k,j,c,d)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !v_vooo(c,j,k,l) * t2(i,l,a,b) &
+  !X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
+
+  !$OMP DO collapse(3)
+  do k = 1, nO
+    do j = 1, nO
+      do c = 1, nV
+        do l = 1, nO
+           X_ovoo(l,c,j,k) = v_vooo(c,j,k,l)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(3)
+  do i = 1, nO
+    do b = 1, nV
+      do a = 1, nV
+        do l = 1, nO
+          T_ovvo(l,a,b,i) = t2(i,l,a,b)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !v_vvoo(b,c,j,k) * t1(i,a) &
+  !X_vvoo(b,c,k,j) * T1_vo(a,i) &
+
+  !$OMP DO collapse(3)
+  do j = 1, nO
+    do k = 1, nO
+      do c = 1, nV
+        do b = 1, nV
+          X_vvoo(b,c,k,j) = v_vvoo(b,c,j,k)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(1)
+  do i = 1, nO
+    do a = 1, nV
+      T_vo(a,i) = t1(i,a)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call wall_time(ta)
+  energy = 0d0
+  do i = 1, nO
+    do j = 1, nO
+      do k = 1, nO
+        delta_ijk = f_o(i) + f_o(j) + f_o(k)
+        call form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_ijk)
+        call form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,W_ijk,V_ijk)
+        !$OMP PARALLEL &
+        !$OMP SHARED(energy,nV,i,j,k,W_ijk,V_ijk,f_o,f_v,delta_ijk) &
+        !$OMP PRIVATE(a,b,c,e,delta) &
+        !$OMP DEFAULT(NONE)
+        e = 0d0
+        !$OMP DO
+        do c = 1, nV
+          do b = 1, nV
+            do a = 1, nV
+              delta = 1d0 / (delta_ijk - f_v(a) - f_v(b) - f_v(c))
+              !energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c))  !delta_ooovvv(i,j,k,a,b,c)
+              e = e + (4d0 * W_ijk(a,b,c) + W_ijk(b,c,a) + W_ijk(c,a,b)) &
+                       * (V_ijk(a,b,c) - V_ijk(c,b,a)) * delta
+            enddo
+          enddo
+        enddo
+        !$OMP END DO NOWAIT
+        !$OMP CRITICAL
+        energy = energy + e
+        !$OMP END CRITICAL
+        !$OMP END PARALLEL
+      enddo
+    enddo
+    call wall_time(tb)
+    write(*,'(F12.2,A5,F12.2,A2)') dble(i)/dble(nO)*100d0, '% in ', tb - ta, ' s'
+  enddo
+
+  energy = energy / 3d0
+
+  deallocate(W_ijk,V_ijk,X_vvvo,X_ovoo,T_vvoo,T_ovvo,T_vo)
+  !deallocate(V,W)
+end
+
+! W_ijk
+
+subroutine form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV,i,j,k
+  !double precision, intent(in) :: t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO)
+  double precision, intent(in)  :: X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO)
+  double precision, intent(out) :: W(nV,nV,nV)!,nO,nO,nO)
+
+  integer :: l,a,b,c,d
+  double precision, allocatable, dimension(:,:,:) :: X, Y, Z
+
+  !W = 0d0
+  !do i = 1, nO
+  !  do j = 1, nO
+  !    do k = 1, nO
+
+  allocate(X(nV,nV,nV))
+  allocate(Y(nV,nV,nV))
+  allocate(Z(nV,nV,nV))
+
+  !$OMP PARALLEL DO
+  do b = 1, nV
+    do a = 1, nV
+      do d = 1, nV
+        Z(d,a,b) = X_vvvo(d,b,a,i)
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+  call dgemm('T','N',nV*nV,nV,nV, 1.d0, &
+     Z, nV, T_vvoo(1,1,k,j), nV, 0.d0, W, nV*nV)
+
+  !$OMP PARALLEL DO
+  do c = 1, nV
+    do a = 1, nV
+      do d = 1, nV
+        Z(d,a,c) = X_vvvo(d,c,a,i)
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+  call dgemm('T','N',nV*nV,nV,nV, 1.d0, &
+     Z, nV, T_vvoo(1,1,j,k), nV, 0.d0, Y, nV*nV)
+
+  call dgemm('T','N',nV*nV,nV,nV, 1.d0, &
+     X_vvvo(1,1,1,k), nV, T_vvoo(1,1,j,i), nV, 1.d0, Y, nV*nV)
+
+  call dgemm('T','N',nV,nV*nV,nV, 1.d0, &
+     T_vvoo(1,1,i,j), nV, X_vvvo(1,1,1,k), nV, 1.d0, W, nV)
+
+  call dgemm('T','N',nV,nV*nV,nV, 1.d0, &
+     T_vvoo(1,1,i,k), nV, X_vvvo(1,1,1,j), nV, 1.d0, Y, nV)
+
+  call dgemm('T','N',nV*nV,nV,nV, 1.d0, &
+     X_vvvo(1,1,1,j), nV, T_vvoo(1,1,k,i), nV, 1.d0, W, nV*nV)
+
+  deallocate(Z)
+
+
+  allocate(Z(nO,nV,nV))
+
+  call dgemm('T','N',nV*nV,nV,nO, -1.d0, &
+     T_ovvo(1,1,1,i), nO, X_ovoo(1,1,j,k), nO, 1.d0, W, nV*nV)
+
+  call dgemm('T','N',nV*nV,nV,nO, -1.d0, &
+     T_ovvo(1,1,1,i), nO, X_ovoo(1,1,k,j), nO, 1.d0, Y, nV*nV)
+
+  !$OMP PARALLEL DO
+  do c = 1, nV
+    do a = 1, nV
+      do l = 1, nO
+        Z(l,a,c) = T_ovvo(l,c,a,k)
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+  call dgemm('T','N',nV*nV,nV,nO, -1.d0, &
+     Z, nO, X_ovoo(1,1,i,j), nO, 1.d0, Y, nV*nV)
+
+  call dgemm('T','N',nV,nV*nV,nO, -1.d0, &
+     X_ovoo(1,1,j,i), nO, T_ovvo(1,1,1,k), nO, 1.d0, Y, nV)
+
+  call dgemm('T','N',nV,nV*nV,nO, -1.d0, &
+     X_ovoo(1,1,k,i), nO, T_ovvo(1,1,1,j), nO, 1.d0, W, nV)
+
+  !$OMP PARALLEL DO
+  do b = 1, nV
+    do a = 1, nV
+      do l = 1, nO
+        Z(l,a,b) = T_ovvo(l,b,a,j)
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+  call dgemm('T','N',nV*nV,nV,nO, -1.d0, &
+     Z, nO, X_ovoo(1,1,i,k), nO, 1.d0, W, nV*nV)
+
+  !$OMP PARALLEL DO
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        W(a,b,c) = W(a,b,c) + Y(a,c,b)
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+  deallocate(X,Y,Z)
+
+
+!  !$OMP PARALLEL &
+!  !$OMP SHARED(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W) &
+!  !$OMP PRIVATE(a,b,c,d,l) &
+!  !$OMP DEFAULT(NONE)
+!
+!  !$OMP DO collapse(2)
+!  do c = 1, nV
+!    do b = 1, nV
+!      do a = 1, nV
+!        W(a,b,c) = 0.d0
+!
+!        do d = 1, nV
+!          !W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+!          W(a,b,c) = W(a,b,c) &
+!          ! chem (bd|ai)
+!          ! phys <ba|di>
+!          !+ cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+!          !+ cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
+!          !+ cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
+!          !+ cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
+!          !+ cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
+!          !+ cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
+!         + X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j) &
+!         + X_vvvo(d,c,a,i) * T_vvoo(d,b,j,k) & ! bc kj
+!         + X_vvvo(d,a,c,k) * T_vvoo(d,b,j,i) & ! prev ac ik
+!          + X_vvvo(d,b,c,k) * T_vvoo(d,a,i,j) & ! prev ab ij
+!          + X_vvvo(d,c,b,j) * T_vvoo(d,a,i,k) & ! prev bc kj
+!          + X_vvvo(d,a,b,j) * T_vvoo(d,c,k,i) ! prev ac ik
+!        enddo
+!
+!      enddo
+!    enddo
+!  enddo
+!  !$OMP END DO nowait
+!
+!  !$OMP DO collapse(2)
+!  do c = 1, nV
+!    do b = 1, nV
+!      do a = 1, nV
+!
+!        do l = 1, nO
+!          !W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+!          W(a,b,c) = W(a,b,c) &
+!          ! chem (ck|jl)
+!          ! phys <cj|kl>
+!          !- cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
+!          !- cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
+!          !- cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
+!          !- cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
+!          !- cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
+!          !- cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
+!          - T_ovvo(l,a,b,i) * X_ovoo(l,c,j,k) &
+!          - T_ovvo(l,a,c,i) * X_ovoo(l,b,k,j) & ! bc kj
+!          - T_ovvo(l,c,a,k) * X_ovoo(l,b,i,j) & ! prev ac ik
+!          - T_ovvo(l,c,b,k) * X_ovoo(l,a,j,i) & ! prev ab ij
+!          - T_ovvo(l,b,c,j) * X_ovoo(l,a,k,i) & ! prev bc kj
+!          - T_ovvo(l,b,a,j) * X_ovoo(l,c,i,k) ! prev ac ik
+!        enddo
+!
+!      enddo
+!    enddo
+!  enddo
+!  !$OMP END DO
+!  !$OMP END PARALLEL
+
+  !    enddo
+  !  enddo
+  !enddo
+
+end
+
+! V_ijk
+
+subroutine form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,w,v)
+
+implicit none
+
+  integer, intent(in)           :: nO,nV,i,j,k
+  !double precision, intent(in)  :: t1(nO,nV)
+  double precision, intent(in)  :: T_vo(nV,nO)
+  double precision, intent(in)  :: X_vvoo(nV,nV,nO,nO)
+  double precision, intent(in)  :: W(nV,nV,nV)!,nO,nO,nO)
+  double precision, intent(out) :: V(nV,nV,nV)!,nO,nO,nO)
+
+  integer :: a,b,c
+
+  !V = 0d0
+  !do i = 1, nO
+  !  do j = 1, nO
+  !    do k = 1, nO
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,i,j,k,T_vo,X_vvoo,W,V) &
+  !$OMP PRIVATE(a,b,c) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(2)
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        !V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+        V(a,b,c) = W(a,b,c) &
+        !+ cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
+        !+ cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
+        !+ cc_space_v_vvoo(a,b,i,j) * t1(k,c)
+        + X_vvoo(b,c,k,j) * T_vo(a,i) &
+        + X_vvoo(a,c,k,i) * T_vo(b,j) &
+        + X_vvoo(a,b,j,i) * T_vo(c,k)
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  !    enddo
+  !  enddo
+  !enddo
+
+end
+

src/ccsd/ccsd_t_space_orb_abc.irp.f

@@ -0,0 +1,252 @@
+! Main
+
+subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO,nV), f_o(nO), f_v(nV)
+  double precision, intent(in)  :: t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO), v_vooo(nV,nO,nO,nO)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: W(:,:,:,:,:,:)
+  double precision, allocatable :: V(:,:,:,:,:,:)
+  double precision, allocatable :: W_abc(:,:,:), V_abc(:,:,:)
+  double precision, allocatable :: W_cab(:,:,:), W_cba(:,:,:)
+  double precision, allocatable :: W_bca(:,:,:), V_cba(:,:,:)
+  double precision, allocatable :: X_vvvo(:,:,:,:), X_ovoo(:,:,:,:), X_vvoo(:,:,:,:)
+  double precision, allocatable :: T_vvoo(:,:,:,:), T_ovvo(:,:,:,:), T_vo(:,:)
+  integer                       :: i,j,k,l,a,b,c,d
+  double precision              :: e,ta,tb, delta, delta_abc
+
+  !allocate(W(nV,nV,nV,nO,nO,nO))
+  !allocate(V(nV,nV,nV,nO,nO,nO))
+  allocate(W_abc(nO,nO,nO), V_abc(nO,nO,nO), W_cab(nO,nO,nO))
+  allocate(W_bca(nO,nO,nO), V_cba(nO,nO,nO), W_cba(nO,nO,nO))
+  allocate(X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO), X_vvoo(nV,nV,nO,nO))
+  allocate(T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO), T_vo(nV,nO))
+
+  ! Temporary arrays
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,T_vvoo,T_ovvo,T_vo,X_vvvo,X_ovoo,X_vvoo, &
+  !$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
+  !$OMP PRIVATE(a,b,c,d,i,j,k,l) &
+  !$OMP DEFAULT(NONE)
+
+  !v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+  !X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j)
+
+  !$OMP DO collapse(3)
+  do i = 1, nO
+    do a = 1, nV
+      do b = 1, nV
+        do d = 1, nV
+          X_vvvo(d,b,a,i) = v_vvvo(b,a,d,i)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(3)
+  do j = 1, nO
+    do k = 1, nO
+      do c = 1, nV
+        do d = 1, nV
+          T_vvoo(d,c,k,j) = t2(k,j,c,d)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !v_vooo(c,j,k,l) * t2(i,l,a,b) &
+  !X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
+
+  !$OMP DO collapse(3)
+  do k = 1, nO
+    do j = 1, nO
+      do c = 1, nV
+        do l = 1, nO
+           X_ovoo(l,c,j,k) = v_vooo(c,j,k,l)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(3)
+  do i = 1, nO
+    do b = 1, nV
+      do a = 1, nV
+        do l = 1, nO
+          T_ovvo(l,a,b,i) = t2(i,l,a,b)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !v_vvoo(b,c,j,k) * t1(i,a) &
+  !X_vvoo(b,c,k,j) * T1_vo(a,i) &
+
+  !$OMP DO collapse(3)
+  do j = 1, nO
+    do k = 1, nO
+      do c = 1, nV
+        do b = 1, nV
+          X_vvoo(b,c,k,j) = v_vvoo(b,c,j,k)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(1)
+  do i = 1, nO
+    do a = 1, nV
+      T_vo(a,i) = t1(i,a)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call wall_time(ta)
+  energy = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        call form_w_abc(nO,nV,a,b,c,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_abc)
+        call form_w_abc(nO,nV,b,c,a,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_bca)
+        call form_w_abc(nO,nV,c,a,b,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_cab)
+        call form_w_abc(nO,nV,c,b,a,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_cba)
+
+        call form_v_abc(nO,nV,a,b,c,T_vo,X_vvoo,W_abc,V_abc)
+        call form_v_abc(nO,nV,c,b,a,T_vo,X_vvoo,W_cba,V_cba)
+        !$OMP PARALLEL                                               &
+            !$OMP SHARED(energy,nO,a,b,c,W_abc,W_cab,W_bca,V_abc,V_cba,f_o,f_v,delta_abc)&
+            !$OMP PRIVATE(i,j,k,e,delta)                             &
+            !$OMP DEFAULT(NONE)
+        e = 0d0
+        !$OMP DO
+        do i = 1, nO
+          do j = 1, nO
+            do k = 1, nO
+              delta = 1d0 / (f_o(i) + f_o(j) + f_o(k) - delta_abc)
+              !energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c))  !delta_ooovvv(i,j,k,a,b,c)
+              e = e + (4d0 * W_abc(i,j,k) + W_bca(i,j,k) + W_cab(i,j,k))&
+                  * (V_abc(i,j,k) - V_cba(i,j,k)) * delta
+            enddo
+          enddo
+        enddo
+        !$OMP END DO NOWAIT
+        !$OMP CRITICAL
+        energy = energy + e
+        !$OMP END CRITICAL
+        !$OMP END PARALLEL
+      enddo
+    enddo
+    call wall_time(tb)
+    write(*,'(F12.2,A5,F12.2,A2)') dble(i)/dble(nO)*100d0, '% in ', tb - ta, ' s'
+  enddo
+
+  energy = energy / 3d0
+
+  deallocate(W_abc,V_abc,W_cab,V_cba,W_bca,X_vvvo,X_ovoo,T_vvoo,T_ovvo,T_vo)
+  !deallocate(V,W)
+end
+
+
+subroutine form_w_abc(nO,nV,a,b,c,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_abc)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV,a,b,c
+  !double precision, intent(in) :: t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO)
+  double precision, intent(in)  :: X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO)
+  double precision, intent(out) :: W_abc(nO,nO,nO)
+
+  integer :: l,i,j,k,d
+
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,a,b,c,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_abc) &
+  !$OMP PRIVATE(i,j,k,d,l) &
+  !$OMP DEFAULT(NONE)
+
+  !$OMP DO collapse(3)
+  do k = 1, nO
+    do j = 1, nO
+      do i = 1, nO
+        W_abc(i,j,k) = 0.d0
+
+        do d = 1, nV
+          W_abc(i,j,k) = W_abc(i,j,k) &
+                 + X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j) &
+                 + X_vvvo(d,c,a,i) * T_vvoo(d,b,j,k) &
+                 + X_vvvo(d,a,c,k) * T_vvoo(d,b,j,i) &
+                 + X_vvvo(d,b,c,k) * T_vvoo(d,a,i,j) &
+                 + X_vvvo(d,c,b,j) * T_vvoo(d,a,i,k) &
+                 + X_vvvo(d,a,b,j) * T_vvoo(d,c,k,i)
+
+        enddo
+
+        do l = 1, nO
+          W_abc(i,j,k) = W_abc(i,j,k) &
+              - T_ovvo(l,a,b,i) * X_ovoo(l,c,j,k) &
+              - T_ovvo(l,a,c,i) * X_ovoo(l,b,k,j) & ! bc kj
+              - T_ovvo(l,c,a,k) * X_ovoo(l,b,i,j) & ! prev ac ik
+              - T_ovvo(l,c,b,k) * X_ovoo(l,a,j,i) & ! prev ab ij
+              - T_ovvo(l,b,c,j) * X_ovoo(l,a,k,i) & ! prev bc kj
+              - T_ovvo(l,b,a,j) * X_ovoo(l,c,i,k) ! prev ac ik
+        enddo
+
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+
+end
+
+
+! V_abc
+
+subroutine form_v_abc(nO,nV,a,b,c,T_vo,X_vvoo,W,V)
+
+implicit none
+
+  integer, intent(in)           :: nO,nV,a,b,c
+  !double precision, intent(in)  :: t1(nO,nV)
+  double precision, intent(in)  :: T_vo(nV,nO)
+  double precision, intent(in)  :: X_vvoo(nV,nV,nO,nO)
+  double precision, intent(in)  :: W(nO,nO,nO)
+  double precision, intent(out) :: V(nO,nO,nO)
+
+  integer :: i,j,k
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,a,b,c,T_vo,X_vvoo,W,V) &
+  !$OMP PRIVATE(i,j,k) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(2)
+  do k = 1, nO
+    do j = 1, nO
+      do i = 1, nO
+        !V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+        V(i,j,k) = W(i,j,k) &
+           + X_vvoo(b,c,k,j) * T_vo(a,i) &
+           + X_vvoo(a,c,k,i) * T_vo(b,j) &
+           + X_vvoo(a,b,j,i) * T_vo(c,k)
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+end
+

src/ccsd/ccsd_t_spin_orb.irp.f

@@ -0,0 +1,376 @@
+! v1
+
+subroutine ccsd_par_t_spin(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,v_vvvo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: f_o(nO), f_v(nV), f_ov(nO,nV)
+  double precision, intent(in)  :: v_ooov(nO,nO,nO,nV)
+  double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO), v_vvvo(nV,nV,nV,nO)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: t3(:,:,:,:,:,:), s(:,:)
+  double precision              :: e_t, e_st, e_dt, delta_abc, delta
+  integer                       :: i,j,k,l,m,a,b,c,d,e
+
+  allocate(t3(nO,nO,nO,nV,nV,nV), s(nO,nV))
+
+  t3 = 0d0
+  
+  ! T3
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              do e = 1, nV
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+                  + t2(j,k,a,e) * v_vvvo(b,c,e,i) &
+                  - t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
+                  - t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
+                  - t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
+                  - t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
+                  + t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
+                  + t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
+                  + t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
+                  + t2(j,i,c,e) * v_vvvo(b,a,e,k)   ! + P(ik) P(ac)
+              enddo
+              do m = 1, nO
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+                  + t2(m,i,b,c) * v_ooov(j,k,m,a) &
+                  - t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
+                  - t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
+                  - t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
+                  - t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
+                  + t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
+                  + t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
+                  + t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
+                  + t2(m,k,b,a) * v_ooov(j,i,m,c)   ! + P(ik) P(ac)
+              enddo
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) * (1d0 / delta)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  
+
+  ! E_T
+  e_t = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              e_t = e_t + t3(i,j,k,a,b,c) * delta * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  e_t = e_t / 36d0
+
+  ! E_ST
+  s = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  e_st = 0d0
+  do a = 1, nV
+    do i = 1, nO
+      e_st = e_st + s(i,a) * t1(i,a)
+    enddo
+  enddo
+  e_st = e_st * 0.25d0
+
+  ! E_DT
+  e_dt = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              e_dt = e_dt + t2(i,j,a,b) * f_ov(k,c) * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  e_dt = e_dt * 0.25d0
+
+  ! (T)
+  !print*,e_t,e_st,e_dt
+  energy = e_t + e_st + e_dt
+
+  deallocate(t3,s)
+  
+end
+
+! v2
+
+subroutine ccsd_par_t_spin_v2(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: f_o(nO), f_v(nV), f_ov(nO,nV)
+  double precision, intent(in)  :: v_ooov(nO,nO,nO,nV)
+  double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: t3_bc(:,:,:,:), s(:,:), e_t(:), e_dt(:)
+  double precision, allocatable :: A_vovv(:,:,:,:), v_vvvo(:,:,:,:)
+  double precision, allocatable :: T_voov(:,:,:,:), B_ooov(:,:,:,:)
+  double precision              :: e_st, delta_abc, delta, ta, tb
+  integer                       :: i,j,k,l,m,a,b,c,d,e
+
+  allocate(t3_bc(nO,nO,nO,nV), s(nO,nV), e_t(nV), e_dt(nV))
+  allocate(A_vovv(nV,nO,nV,nV),v_vvvo(nV,nV,nV,nO),T_voov(nV,nO,nO,nV),B_ooov(nO,nO,nO,nV))
+
+  call gen_v_spin(cc_nV_m,cc_nV_m,cc_nV_m,cc_nO_m, &
+       cc_nV_S,cc_nV_S,cc_nV_S,cc_nO_S, &
+       cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin, &
+       nV,nV,nV,nO, v_vvvo)
+  
+  ! Init 
+  s = 0d0
+  e_t = 0d0
+  e_st = 0d0
+  e_dt = 0d0
+
+  call wall_time(ta)
+  !$OMP PARALLEL &
+  !$OMP PRIVATE(i,j,k,m,a,b,c,e) &
+  !$OMP SHARED(A_vovv,ta,tb,t3_bc,s,e_t,e_st,e_dt,t2,v_vvvo,v_ooov, &
+  !$OMP v_vvoo,f_o,f_v,f_ov,delta,delta_abc,nO,nV,T_voov,B_ooov) &
+  !$OMP DEFAULT(NONE)
+
+  !$OMP DO collapse(3)
+  do c = 1, nV
+    do b = 1, nV
+      do i = 1, nO
+        do e = 1, nV
+          A_vovv(e,i,b,c) = v_vvvo(b,c,e,i)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$omp do collapse(3)
+  do a = 1, nV
+    do k = 1, nO
+      do j = 1, nO
+        do e = 1, nV
+          T_voov(e,j,k,a) = t2(j,k,a,e)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$omp end do nowait
+
+  !$omp do collapse(3)
+  do a = 1, nV
+    do k = 1, nO
+      do j = 1, nO
+        do m = 1, nO
+          B_ooov(m,j,k,a) = v_ooov(j,k,m,a)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$omp end do
+
+  do c = 1, nV
+    do b = 1, nV
+
+      ! T3(:,:,:,:,b,c)
+      ! Init
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              t3_bc(i,j,k,a) = 0d0
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              do e = 1, nV
+                t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
+                   !+ t2(j,k,a,e) * v_vvvo(b,c,e,i) &
+                   !- t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
+                   !- t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
+                   !- t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
+                   !- t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
+                   !+ t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
+                   !+ t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
+                   !+ t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
+                   !+ t2(j,i,c,e) * v_vvvo(b,a,e,k)   ! + P(ik) P(ac)
+                   + T_voov(e,j,k,a) * A_vovv(e,i,b,c) &
+                   - T_voov(e,i,k,a) * A_vovv(e,j,b,c) & ! - P(ij)
+                   - T_voov(e,j,i,a) * A_vovv(e,k,b,c) & ! - P(ik)
+                   - T_voov(e,j,k,b) * A_vovv(e,i,a,c) & ! - P(ab)
+                   - T_voov(e,j,k,c) * A_vovv(e,i,b,a) & ! - P(ac)
+                   + T_voov(e,i,k,b) * A_vovv(e,j,a,c) & ! + P(ij) P(ab)
+                   + T_voov(e,i,k,c) * A_vovv(e,j,b,a) & ! + P(ij) P(ac)
+                   + T_voov(e,j,i,b) * A_vovv(e,k,a,c) & ! + P(ik) P(ab)
+                   + T_voov(e,j,i,c) * A_vovv(e,k,b,a)   ! + P(ik) P(ac)
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              do m = 1, nO
+                t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
+                   !+ t2(m,i,b,c) * v_ooov(j,k,m,a) &
+                   !- t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
+                   !- t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
+                   !- t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
+                   !- t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
+                   !+ t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
+                   !+ t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
+                   !+ t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
+                   !+ t2(m,k,b,a) * v_ooov(j,i,m,c)   ! + P(ik) P(ac)
+                   + t2(m,i,b,c) * B_ooov(m,j,k,a) &
+                   - t2(m,j,b,c) * B_ooov(m,i,k,a) & ! - P(ij)
+                   - t2(m,k,b,c) * B_ooov(m,j,i,a) & ! - P(ik)
+                   - t2(m,i,a,c) * B_ooov(m,j,k,b) & ! - P(ab)
+                   - t2(m,i,b,a) * B_ooov(m,j,k,c) & ! - P(ac)
+                   + t2(m,j,a,c) * B_ooov(m,i,k,b) & ! + P(ij) P(ab)
+                   + t2(m,j,b,a) * B_ooov(m,i,k,c) & ! + P(ij) P(ac)
+                   + t2(m,k,a,c) * B_ooov(m,j,i,b) & ! + P(ik) P(ab)
+                   + t2(m,k,b,a) * B_ooov(m,j,i,c)   ! + P(ik) P(ac)
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+
+      !$OMP DO
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+               delta = f_o(i) + f_o(j) + f_o(k) - delta_abc 
+               t3_bc(i,j,k,a) = t3_bc(i,j,k,a) * (1d0 / delta)
+             enddo
+           enddo
+         enddo
+       enddo
+      !$OMP END DO
+
+      ! E_T
+      !$OMP DO
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              e_t(a) = e_t(a) + t3_bc(i,j,k,a) * delta * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO nowait
+
+      ! E_ST
+      !$OMP DO
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO nowait
+
+      ! E_DT
+      !$OMP DO
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              e_dt(a) = e_dt(a) + t2(i,j,a,b) * f_ov(k,c) * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+    enddo
+    !$OMP MASTER
+    call wall_time(tb)
+    write(*,'(A1,F6.2,A5,F10.2,A2)') ' ', dble(c)/dble(nV)*100d0, '% in ', tb-ta, ' s'
+    !$OMP END MASTER
+  enddo
+  !$OMP END PARALLEL
+
+  do a = 2, nV
+    e_t(1) = e_t(1) + e_t(a)
+  enddo
+  
+  do a = 2, nV
+    e_dt(1) = e_dt(1) + e_dt(a)
+  enddo
+
+  e_t = e_t / 36d0
+  
+  do a = 1, nV
+    do i = 1, nO
+      e_st = e_st + s(i,a) * t1(i,a)
+    enddo
+  enddo
+  e_st = e_st * 0.25d0
+
+  e_dt = e_dt * 0.25d0
+
+  ! (T)
+  !print*,e_t(1),e_st,e_dt(1)
+  energy = e_t(1) + e_st + e_dt(1)
+
+  deallocate(t3_bc,s)
+  
+end

src/ccsd/org/ccsd_t_space_orb.org

@@ -0,0 +1,428 @@
+Ref:
+Integral-Direct and Parallel Implementation of the CCSD(T) Method:
+Algorithmic Developments and Large-Scale Applications
+László Gyevi-Nagy, Mihály Kállay, and Péter R. Nagy
+J. Chem. Theory Comput. 2020, 16, 1, 366–384
+https://doi.org/10.1021/acs.jctc.9b00957
+
+* Dumb way
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine ccsd_par_t_space(nO,nV,t1,t2,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO, nV)
+  double precision, intent(in)  :: t2(nO, nO, nV, nV)
+  double precision, intent(out) :: energy
+  
+  double precision, allocatable :: W(:,:,:,:,:,:)
+  double precision, allocatable :: V(:,:,:,:,:,:)
+  integer :: i,j,k,a,b,c
+  
+  allocate(W(nO,nO,nO,nV,nV,nV))
+  allocate(V(nO,nO,nO,nV,nV,nV))
+
+  call form_w(nO,nV,t2,W)   
+  call form_v(nO,nV,t1,W,V)
+
+  energy = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c))  !delta_ooovvv(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  
+  energy = energy / 3d0
+  
+  deallocate(V,W)
+end
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine form_w(nO,nV,t2,W)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t2(nO, nO, nV, nV)
+  double precision, intent(out) :: W(nO, nO, nO, nV, nV, nV)
+  
+  integer :: i,j,k,l,a,b,c,d
+
+  W = 0d0
+  do c = 1, nV
+    print*,'W:',c,'/',nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+
+              do d = 1, nV
+                W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+                ! chem (bd|ai)
+                ! phys <ba|di>
+                + cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+                + cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
+                + cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
+                + cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
+                + cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
+                + cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
+              enddo
+
+              do l = 1, nO
+                W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+                ! chem (ck|jl)
+                ! phys <cj|kl>
+                - cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
+                - cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
+                - cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
+                - cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
+                - cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
+                - cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
+              enddo
+
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine form_v(nO,nV,t1,w,v)
+
+implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO, nV)
+  double precision, intent(in)  :: W(nO, nO, nO, nV, nV, nV)
+  double precision, intent(out) :: V(nO, nO, nO, nV, nV, nV)
+
+  integer :: i,j,k,a,b,c
+
+  V = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+              + cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
+              + cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
+              + cc_space_v_vvoo(a,b,i,j) * t1(k,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end
+#+END_SRC
+
+* Better way
+** Main
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine ccsd_par_t_space_v2(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO,nV), f_o(nO), f_v(nV)
+  double precision, intent(in)  :: t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO), v_vooo(nV,nO,nO,nO)
+  double precision, intent(out) :: energy
+  
+  double precision, allocatable :: W(:,:,:,:,:,:)
+  double precision, allocatable :: V(:,:,:,:,:,:)
+  double precision, allocatable :: W_ijk(:,:,:), V_ijk(:,:,:)
+  double precision, allocatable :: X_vvvo(:,:,:,:), X_ovoo(:,:,:,:), X_vvoo(:,:,:,:)
+  double precision, allocatable :: T_vvoo(:,:,:,:), T_ovvo(:,:,:,:), T_vo(:,:)
+  integer                       :: i,j,k,l,a,b,c,d
+  double precision              :: e,ta,tb, delta, delta_ijk
+ 
+  !allocate(W(nV,nV,nV,nO,nO,nO))
+  !allocate(V(nV,nV,nV,nO,nO,nO))
+  allocate(W_ijk(nV,nV,nV), V_ijk(nV,nV,nV))
+  allocate(X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO), X_vvoo(nV,nV,nO,nO))
+  allocate(T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO), T_vo(nV,nO))
+
+  ! Temporary arrays
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,T_vvoo,T_ovvo,T_vo,X_vvvo,X_ovoo,X_vvoo, &
+  !$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
+  !$OMP PRIVATE(a,b,c,d,i,j,k,l) &
+  !$OMP DEFAULT(NONE)
+  
+  !v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+  !X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j)
+  
+  !$OMP DO collapse(3)
+  do i = 1, nO
+    do a = 1, nV
+      do b = 1, nV
+        do d = 1, nV
+          X_vvvo(d,b,a,i) = v_vvvo(b,a,d,i)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(3)
+  do j = 1, nO
+    do k = 1, nO
+      do c = 1, nV
+        do d = 1, nV
+          T_vvoo(d,c,k,j) = t2(k,j,c,d)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+ 
+  !v_vooo(c,j,k,l) * t2(i,l,a,b) &
+  !X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
+
+  !$OMP DO collapse(3)
+  do k = 1, nO
+    do j = 1, nO
+      do c = 1, nV
+        do l = 1, nO
+           X_ovoo(l,c,j,k) = v_vooo(c,j,k,l)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(3)
+  do i = 1, nO
+    do b = 1, nV
+      do a = 1, nV
+        do l = 1, nO
+          T_ovvo(l,a,b,i) = t2(i,l,a,b)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+                     
+  !v_vvoo(b,c,j,k) * t1(i,a) &
+  !X_vvoo(b,c,k,j) * T1_vo(a,i) &
+  
+  !$OMP DO collapse(3)
+  do j = 1, nO
+    do k = 1, nO
+      do c = 1, nV
+        do b = 1, nV
+          X_vvoo(b,c,k,j) = v_vvoo(b,c,j,k)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(1)
+  do i = 1, nO
+    do a = 1, nV
+      T_vo(a,i) = t1(i,a)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call wall_time(ta)
+  energy = 0d0
+  do i = 1, nO
+    do j = 1, nO
+      do k = 1, nO
+        delta_ijk = f_o(i) + f_o(j) + f_o(k)
+        call form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_ijk)
+        call form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,W_ijk,V_ijk)
+        !$OMP PARALLEL &
+        !$OMP SHARED(energy,nV,i,j,k,W_ijk,V_ijk,f_o,f_v,delta_ijk) &
+        !$OMP PRIVATE(a,b,c,e,delta) &
+        !$OMP DEFAULT(NONE)
+        e = 0d0
+        !$OMP DO
+        do c = 1, nV
+          do b = 1, nV
+            do a = 1, nV
+              delta = 1d0 / (delta_ijk - f_v(a) - f_v(b) - f_v(c))
+              !energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c))  !delta_ooovvv(i,j,k,a,b,c)
+              e = e + (4d0 * W_ijk(a,b,c) + W_ijk(b,c,a) + W_ijk(c,a,b)) &
+                       * (V_ijk(a,b,c) - V_ijk(c,b,a)) * delta
+            enddo
+          enddo
+        enddo
+        !$OMP END DO
+        !$OMP CRITICAL
+        energy = energy + e
+        !$OMP END CRITICAL
+        !$OMP END PARALLEL
+      enddo
+    enddo
+    call wall_time(tb)
+    write(*,'(F12.2,A5,F12.2,A2)') dble(i)/dble(nO)*100d0, '% in ', tb - ta, ' s'
+  enddo
+  
+  energy = energy / 3d0
+
+  deallocate(W_ijk,V_ijk,X_vvvo,X_ovoo,T_vvoo,T_ovvo,T_vo)
+  !deallocate(V,W)
+end
+#+END_SRC
+
+** W_ijk
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV,i,j,k
+  !double precision, intent(in) :: t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO)
+  double precision, intent(in)  :: X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO)
+  double precision, intent(out) :: W(nV,nV,nV)!,nO,nO,nO)
+  
+  integer :: l,a,b,c,d
+
+  !W = 0d0
+  !do i = 1, nO
+  !  do j = 1, nO
+  !    do k = 1, nO
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W) &
+  !$OMP PRIVATE(a,b,c,d,l) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(2)
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        W(a,b,c) = 0d0
+
+        do d = 1, nV
+          !W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+          W(a,b,c) = W(a,b,c) &
+          ! chem (bd|ai)
+          ! phys <ba|di>
+          !+ cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+          !+ cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
+          !+ cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
+          !+ cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
+          !+ cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
+          !+ cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
+          + X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j) &
+          + X_vvvo(d,c,a,i) * T_vvoo(d,b,j,k) & ! bc kj
+          + X_vvvo(d,a,c,k) * T_vvoo(d,b,j,i) & ! prev ac ik
+          + X_vvvo(d,b,c,k) * T_vvoo(d,a,i,j) & ! prev ab ij
+          + X_vvvo(d,c,b,j) * T_vvoo(d,a,i,k) & ! prev bc kj
+          + X_vvvo(d,a,b,j) * T_vvoo(d,c,k,i) ! prev ac ik
+        enddo
+        
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(2)
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+         
+        do l = 1, nO
+          !W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+          W(a,b,c) = W(a,b,c) &
+          ! chem (ck|jl)
+          ! phys <cj|kl>
+          !- cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
+          !- cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
+          !- cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
+          !- cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
+          !- cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
+          !- cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
+          - X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
+          - X_ovoo(l,b,k,j) * T_ovvo(l,a,c,i) & ! bc kj
+          - X_ovoo(l,b,i,j) * T_ovvo(l,c,a,k) & ! prev ac ik
+          - X_ovoo(l,a,j,i) * T_ovvo(l,c,b,k) & ! prev ab ij
+          - X_ovoo(l,a,k,i) * T_ovvo(l,b,c,j) & ! prev bc kj
+          - X_ovoo(l,c,i,k) * T_ovvo(l,b,a,j) ! prev ac ik
+        enddo
+
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  
+  !    enddo
+  !  enddo
+  !enddo
+
+end
+#+END_SRC
+
+** V_ijk
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,w,v)
+
+implicit none
+
+  integer, intent(in)           :: nO,nV,i,j,k
+  !double precision, intent(in)  :: t1(nO,nV)
+  double precision, intent(in)  :: T_vo(nV,nO)
+  double precision, intent(in)  :: X_vvoo(nV,nV,nO,nO)
+  double precision, intent(in)  :: W(nV,nV,nV)!,nO,nO,nO)
+  double precision, intent(out) :: V(nV,nV,nV)!,nO,nO,nO)
+
+  integer :: a,b,c
+
+  !V = 0d0
+  !do i = 1, nO
+  !  do j = 1, nO
+  !    do k = 1, nO
+  
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,i,j,k,T_vo,X_vvoo,W,V) &
+  !$OMP PRIVATE(a,b,c) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(2)
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        !V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+        V(a,b,c) = W(a,b,c) &
+        !+ cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
+        !+ cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
+        !+ cc_space_v_vvoo(a,b,i,j) * t1(k,c)
+        + X_vvoo(b,c,k,j) * T_vo(a,i) &
+        + X_vvoo(a,c,k,i) * T_vo(b,j) &
+        + X_vvoo(a,b,j,i) * T_vo(c,k)
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  
+  !    enddo
+  !  enddo
+  !enddo
+
+end
+#+END_SRC

src/ccsd/org/ccsd_t_spin_orb.org

@@ -0,0 +1,385 @@
+* CCSD(T) spin orb
+
+Ref:
+John D. Watts, Jürgen Gauss, and Rodney J. Bartlett
+J. Chem. Phys. 98, 8718 (1993)
+http://dx.doi.org/10.1063/1.464480
+
+** v1
+#+begin_src f90 :comments org :tangle ccsd_t_spin_orb.irp.f
+subroutine ccsd_par_t_spin(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,v_vvvo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: f_o(nO), f_v(nV), f_ov(nO,nV)
+  double precision, intent(in)  :: v_ooov(nO,nO,nO,nV)
+  double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO), v_vvvo(nV,nV,nV,nO)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: t3(:,:,:,:,:,:), s(:,:)
+  double precision              :: e_t, e_st, e_dt, delta_abc, delta
+  integer                       :: i,j,k,l,m,a,b,c,d,e
+
+  allocate(t3(nO,nO,nO,nV,nV,nV), s(nO,nV))
+
+  t3 = 0d0
+  
+  ! T3
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              do e = 1, nV
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+                  + t2(j,k,a,e) * v_vvvo(b,c,e,i) &
+                  - t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
+                  - t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
+                  - t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
+                  - t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
+                  + t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
+                  + t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
+                  + t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
+                  + t2(j,i,c,e) * v_vvvo(b,a,e,k)   ! + P(ik) P(ac)
+              enddo
+              do m = 1, nO
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+                  + t2(m,i,b,c) * v_ooov(j,k,m,a) &
+                  - t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
+                  - t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
+                  - t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
+                  - t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
+                  + t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
+                  + t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
+                  + t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
+                  + t2(m,k,b,a) * v_ooov(j,i,m,c)   ! + P(ik) P(ac)
+              enddo
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) * (1d0 / delta)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  
+
+  ! E_T
+  e_t = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              e_t = e_t + t3(i,j,k,a,b,c) * delta * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  e_t = e_t / 36d0
+
+  ! E_ST
+  s = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  e_st = 0d0
+  do a = 1, nV
+    do i = 1, nO
+      e_st = e_st + s(i,a) * t1(i,a)
+    enddo
+  enddo
+  e_st = e_st * 0.25d0
+
+  ! E_DT
+  e_dt = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              e_dt = e_dt + t2(i,j,a,b) * f_ov(k,c) * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  e_dt = e_dt * 0.25d0
+
+  ! (T)
+  !print*,e_t,e_st,e_dt
+  energy = e_t + e_st + e_dt
+
+  deallocate(t3,s)
+  
+end
+#+end_src
+
+** v2
+#+begin_src f90 :comments org :tangle ccsd_t_spin_orb.irp.f
+subroutine ccsd_par_t_spin_v2(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: f_o(nO), f_v(nV), f_ov(nO,nV)
+  double precision, intent(in)  :: v_ooov(nO,nO,nO,nV)
+  double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: t3_bc(:,:,:,:), s(:,:), e_t(:), e_dt(:)
+  double precision, allocatable :: A_vovv(:,:,:,:), v_vvvo(:,:,:,:)
+  double precision, allocatable :: T_voov(:,:,:,:), B_ooov(:,:,:,:)
+  double precision              :: e_st, delta_abc, delta, ta, tb
+  integer                       :: i,j,k,l,m,a,b,c,d,e
+
+  allocate(t3_bc(nO,nO,nO,nV), s(nO,nV), e_t(nV), e_dt(nV))
+  allocate(A_vovv(nV,nO,nV,nV),v_vvvo(nV,nV,nV,nO),T_voov(nV,nO,nO,nV),B_ooov(nO,nO,nO,nV))
+
+  call gen_v_spin(cc_nV_m,cc_nV_m,cc_nV_m,cc_nO_m, &
+       cc_nV_S,cc_nV_S,cc_nV_S,cc_nO_S, &
+       cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin, &
+       nV,nV,nV,nO, v_vvvo)
+  
+  ! Init 
+  s = 0d0
+  e_t = 0d0
+  e_st = 0d0
+  e_dt = 0d0
+
+  call wall_time(ta)
+  !$OMP PARALLEL &
+  !$OMP PRIVATE(i,j,k,m,a,b,c,e) &
+  !$OMP SHARED(A_vovv,ta,tb,t3_bc,s,e_t,e_st,e_dt,t2,v_vvvo,v_ooov, &
+  !$OMP v_vvoo,f_o,f_v,f_ov,delta,delta_abc,nO,nV,T_voov,B_ooov) &
+  !$OMP DEFAULT(NONE)
+
+  !$OMP DO collapse(3)
+  do c = 1, nV
+    do b = 1, nV
+      do i = 1, nO
+        do e = 1, nV
+          A_vovv(e,i,b,c) = v_vvvo(b,c,e,i)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$omp do collapse(3)
+  do a = 1, nV
+    do k = 1, nO
+      do j = 1, nO
+        do e = 1, nV
+          T_voov(e,j,k,a) = t2(j,k,a,e)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$omp end do nowait
+
+  !$omp do collapse(3)
+  do a = 1, nV
+    do k = 1, nO
+      do j = 1, nO
+        do m = 1, nO
+          B_ooov(m,j,k,a) = v_ooov(j,k,m,a)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$omp end do
+
+  do c = 1, nV
+    do b = 1, nV
+
+      ! T3(:,:,:,:,b,c)
+      ! Init
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              t3_bc(i,j,k,a) = 0d0
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              do e = 1, nV
+                t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
+                   !+ t2(j,k,a,e) * v_vvvo(b,c,e,i) &
+                   !- t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
+                   !- t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
+                   !- t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
+                   !- t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
+                   !+ t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
+                   !+ t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
+                   !+ t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
+                   !+ t2(j,i,c,e) * v_vvvo(b,a,e,k)   ! + P(ik) P(ac)
+                   + T_voov(e,j,k,a) * A_vovv(e,i,b,c) &
+                   - T_voov(e,i,k,a) * A_vovv(e,j,b,c) & ! - P(ij)
+                   - T_voov(e,j,i,a) * A_vovv(e,k,b,c) & ! - P(ik)
+                   - T_voov(e,j,k,b) * A_vovv(e,i,a,c) & ! - P(ab)
+                   - T_voov(e,j,k,c) * A_vovv(e,i,b,a) & ! - P(ac)
+                   + T_voov(e,i,k,b) * A_vovv(e,j,a,c) & ! + P(ij) P(ab)
+                   + T_voov(e,i,k,c) * A_vovv(e,j,b,a) & ! + P(ij) P(ac)
+                   + T_voov(e,j,i,b) * A_vovv(e,k,a,c) & ! + P(ik) P(ab)
+                   + T_voov(e,j,i,c) * A_vovv(e,k,b,a)   ! + P(ik) P(ac)
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              do m = 1, nO
+                t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
+                   !+ t2(m,i,b,c) * v_ooov(j,k,m,a) &
+                   !- t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
+                   !- t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
+                   !- t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
+                   !- t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
+                   !+ t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
+                   !+ t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
+                   !+ t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
+                   !+ t2(m,k,b,a) * v_ooov(j,i,m,c)   ! + P(ik) P(ac)
+                   + t2(m,i,b,c) * B_ooov(m,j,k,a) &
+                   - t2(m,j,b,c) * B_ooov(m,i,k,a) & ! - P(ij)
+                   - t2(m,k,b,c) * B_ooov(m,j,i,a) & ! - P(ik)
+                   - t2(m,i,a,c) * B_ooov(m,j,k,b) & ! - P(ab)
+                   - t2(m,i,b,a) * B_ooov(m,j,k,c) & ! - P(ac)
+                   + t2(m,j,a,c) * B_ooov(m,i,k,b) & ! + P(ij) P(ab)
+                   + t2(m,j,b,a) * B_ooov(m,i,k,c) & ! + P(ij) P(ac)
+                   + t2(m,k,a,c) * B_ooov(m,j,i,b) & ! + P(ik) P(ab)
+                   + t2(m,k,b,a) * B_ooov(m,j,i,c)   ! + P(ik) P(ac)
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+
+      !$OMP DO
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+               delta = f_o(i) + f_o(j) + f_o(k) - delta_abc 
+               t3_bc(i,j,k,a) = t3_bc(i,j,k,a) * (1d0 / delta)
+             enddo
+           enddo
+         enddo
+       enddo
+      !$OMP END DO
+
+      ! E_T
+      !$OMP DO
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              e_t(a) = e_t(a) + t3_bc(i,j,k,a) * delta * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO nowait
+
+      ! E_ST
+      !$OMP DO
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO nowait
+
+      ! E_DT
+      !$OMP DO
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              e_dt(a) = e_dt(a) + t2(i,j,a,b) * f_ov(k,c) * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+    enddo
+    !$OMP MASTER
+    call wall_time(tb)
+    write(*,'(A1,F6.2,A5,F10.2,A2)') ' ', dble(c)/dble(nV)*100d0, '% in ', tb-ta, ' s'
+    !$OMP END MASTER
+  enddo
+  !$OMP END PARALLEL
+
+  do a = 2, nV
+    e_t(1) = e_t(1) + e_t(a)
+  enddo
+  
+  do a = 2, nV
+    e_dt(1) = e_dt(1) + e_dt(a)
+  enddo
+
+  e_t = e_t / 36d0
+  
+  do a = 1, nV
+    do i = 1, nO
+      e_st = e_st + s(i,a) * t1(i,a)
+    enddo
+  enddo
+  e_st = e_st * 0.25d0
+
+  e_dt = e_dt * 0.25d0
+
+  ! (T)
+  !print*,e_t(1),e_st,e_dt(1)
+  energy = e_t(1) + e_st + e_dt(1)
+
+  deallocate(t3_bc,s)
+  
+end
+#+end_src

src/cipsi/NEED

@@ -1,5 +1,7 @@
+json
 perturbation
 zmq
 mpi
 iterations
 csf
+mol_properties

src/cipsi/cipsi.irp.f

@@ -16,7 +16,6 @@ subroutine run_cipsi
   double precision, external :: memory_of_double
   PROVIDE H_apply_buffer_allocated
 
-  N_iter = 1
   threshold_generators = 1.d0
   SOFT_TOUCH threshold_generators
 
@@ -76,7 +75,6 @@ subroutine run_cipsi
         )
       write(*,'(A)')  '--------------------------------------------------------------------------------'
 
-
     to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
     to_select = max(N_states_diag, to_select)
     if (do_pt2) then
@@ -106,9 +104,10 @@ subroutine run_cipsi
 
     call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
 
-    call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
+    call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
     call print_extrapolated_energy()
-    N_iter += 1
+    call print_mol_properties()
+    call write_cipsi_json(pt2_data,pt2_data_err)
 
     if (qp_stop()) exit
 
@@ -129,13 +128,13 @@ subroutine run_cipsi
     if (qp_stop()) exit
   enddo
 
-  if (.not.qp_stop()) then
+  ! If stopped because N_det > N_det_max, do an extra iteration to compute the PT2
-    if (N_det < N_det_max) then
+  if ((.not.qp_stop()).and.                                          &
-        call diagonalize_CI
+        (N_det > N_det_max) .and.                                    &
-        call save_wavefunction
+        (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and.    &
-        call save_energy(psi_energy_with_nucl_rep, zeros)
+        (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and.&
-    endif
+        (correlation_energy_ratio <= correlation_energy_ratio_max)   &
-
+        ) then
     if (do_pt2) then
       call pt2_dealloc(pt2_data)
       call pt2_dealloc(pt2_data_err)
@@ -154,10 +153,13 @@ subroutine run_cipsi
     call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
     call print_summary(psi_energy_with_nucl_rep(1:N_states), &
       pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2)
-    call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
+    call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
     call print_extrapolated_energy()
+    call print_mol_properties()
+    call write_cipsi_json(pt2_data,pt2_data_err)
   endif
   call pt2_dealloc(pt2_data)
   call pt2_dealloc(pt2_data_err)
 
 end
+

src/cipsi/environment.irp.f

@@ -7,7 +7,9 @@ BEGIN_PROVIDER [ integer, nthreads_pt2 ]
  character*(32) :: env
  call getenv('QP_NTHREADS_PT2',env)
  if (trim(env) /= '') then
+   call lock_io()
    read(env,*) nthreads_pt2
+   call unlock_io()
    call write_int(6,nthreads_pt2,'Target number of threads for PT2')
  endif
 END_PROVIDER

src/cipsi/selection.irp.f

@@ -312,9 +312,6 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
   end do
   deallocate(indices)
 
-!  !$OMP CRITICAL
-!  print *,  'Step1: ', i_generator, preinteresting(0)
-!  !$OMP END CRITICAL
 
   allocate(banned(mo_num, mo_num,2), bannedOrb(mo_num, 2))
   allocate (mat(N_states, mo_num, mo_num))
@@ -466,17 +463,8 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
           fullinteresting(sze+1) = i
         end if
       end do
-
       allocate (fullminilist (N_int, 2, fullinteresting(0)), &
                     minilist (N_int, 2,     interesting(0)) )
-!      if(pert_2rdm)then
-!        allocate(coef_fullminilist_rev(N_states,fullinteresting(0)))
-!        do i=1,fullinteresting(0)
-!          do j = 1, N_states
-!            coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j)
-!          enddo
-!        enddo
-!      endif
 
       do i=1,fullinteresting(0)
         fullminilist(:,:,i) = psi_det_sorted(:,:,fullinteresting(i))
@@ -524,33 +512,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
             call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
             if(fullMatch) cycle
-! !$OMP CRITICAL
-!  print *,  'Step3: ', i_generator, h1, interesting(0)
-! !$OMP END CRITICAL
 
             call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
-
-!            if(.not.pert_2rdm)then
             call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
-!            else
-!             call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0))
-!            endif
           end if
         enddo
         if(s1 /= s2) monoBdo = .false.
       enddo
       deallocate(fullminilist,minilist)
-!      if(pert_2rdm)then
-!        deallocate(coef_fullminilist_rev)
-!      endif
     enddo
   enddo
   deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)
   deallocate(banned, bannedOrb,mat)
 end subroutine
-
-
-
 subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
   use bitmasks
   use selection_types
@@ -606,18 +580,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 ! to a determinant of the future. In that case, the determinant will be
 ! detected as already generated when generating in the future with a
 ! double excitation.
-!
-!      if (.not.do_singles) then
-!        if ((h1 == p1) .or. (h2 == p2)) then
-!          cycle
-!        endif
-!      endif
-!
-!      if (.not.do_doubles) then
-!        if ((h1 /= p1).and.(h2 /= p2)) then
-!          cycle
-!        endif
-!      endif
 ! -----
 
       if(bannedOrb(p2, s2)) cycle
@@ -974,13 +936,10 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
 
         call get_mask_phase(psi_det_sorted(1,1,interesting(i)), phasemask,N_int)
         if(nt == 4) then
-!          call get_d2_reference(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
           call get_d2(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
         else if(nt == 3) then
-!          call get_d1_reference(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
           call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
         else
-!          call get_d0_reference(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
           call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
         end if
     else if(nt == 4) then
@@ -1540,8 +1499,6 @@ subroutine past_d2(banned, p, sp)
   end if
 end
 
-
-
 subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
   use bitmasks
   implicit none

src/cipsi/stochastic_cipsi.irp.f

@@ -15,7 +15,6 @@ subroutine run_stochastic_cipsi
   double precision, external :: memory_of_double
   PROVIDE H_apply_buffer_allocated distributed_davidson mo_two_e_integrals_in_map
 
-  N_iter = 1
   threshold_generators = 1.d0
   SOFT_TOUCH threshold_generators
 
@@ -96,9 +95,10 @@ subroutine run_stochastic_cipsi
 
     call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
 
-    call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
+    call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
     call print_extrapolated_energy()
-    N_iter += 1
+    call print_mol_properties()
+    call write_cipsi_json(pt2_data,pt2_data_err)
 
     if (qp_stop()) exit
 
@@ -118,13 +118,13 @@ subroutine run_stochastic_cipsi
     if (qp_stop()) exit
   enddo
 
-  if (.not.qp_stop()) then
+  ! If stopped because N_det > N_det_max, do an extra iteration to compute the PT2
-    if (N_det < N_det_max) then
+  if ((.not.qp_stop()).and.                                          &
-        call diagonalize_CI
+        (N_det > N_det_max) .and.                                    &
-        call save_wavefunction
+        (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and.    &
-        call save_energy(psi_energy_with_nucl_rep, zeros)
+        (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and.&
-    endif
+        (correlation_energy_ratio <= correlation_energy_ratio_max)   &
-
+        ) then
     call pt2_dealloc(pt2_data)
     call pt2_dealloc(pt2_data_err)
     call pt2_alloc(pt2_data, N_states)
@@ -134,8 +134,10 @@ subroutine run_stochastic_cipsi
     call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
     call print_summary(psi_energy_with_nucl_rep, &
        pt2_data , pt2_data_err, N_det, N_configuration, N_states, psi_s2)
-    call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
+    call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
     call print_extrapolated_energy()
+    call print_mol_properties()
+    call write_cipsi_json(pt2_data,pt2_data_err)
   endif
   call pt2_dealloc(pt2_data)
   call pt2_dealloc(pt2_data_err)

src/cipsi/write_cipsi_json.irp.f

@@ -0,0 +1,53 @@
+subroutine write_cipsi_json(pt2_data, pt2_data_err)
+    use selection_types
+    implicit none
+    BEGIN_DOC
+!  Writes JSON data for CIPSI runs
+    END_DOC
+    type(pt2_type), intent(in) :: pt2_data, pt2_data_err
+    integer :: i,j,k
+
+    call lock_io
+    character*(64), allocatable    :: fmtk(:)
+    integer :: N_states_p, N_iter_p
+    N_states_p = min(N_states,N_det)
+    N_iter_p = min(N_iter,8)
+    allocate(fmtk(0:N_iter_p))
+    fmtk(:) = '(''   '',E22.15,'','')'
+    fmtk(N_iter_p) = '(''   '',E22.15)'
+
+    write(json_unit, json_dict_uopen_fmt)
+    write(json_unit, json_int_fmt) 'n_det', N_det
+    if (s2_eig) then
+      write(json_unit, json_int_fmt) 'n_cfg', N_configuration
+      if (only_expected_s2) then
+        write(json_unit, json_int_fmt) 'n_csf', N_csf
+      endif
+    endif
+    write(json_unit, json_array_open_fmt) 'states'
+    do k=1,N_states_p
+      write(json_unit, json_dict_uopen_fmt)
+      write(json_unit, json_real_fmt) 'energy', psi_energy_with_nucl_rep(k)
+      write(json_unit, json_real_fmt) 's2', psi_s2(k)
+      write(json_unit, json_real_fmt) 'pt2', pt2_data % pt2(k)
+      write(json_unit, json_real_fmt) 'pt2_err', pt2_data_err % pt2(k)
+      write(json_unit, json_real_fmt) 'rpt2', pt2_data % rpt2(k)
+      write(json_unit, json_real_fmt) 'rpt2_err', pt2_data_err % rpt2(k)
+      write(json_unit, json_real_fmt) 'variance', pt2_data % variance(k)
+      write(json_unit, json_real_fmt) 'variance_err', pt2_data_err % variance(k)
+      write(json_unit, json_array_open_fmt) 'ex_energy'
+      do i=2,N_iter_p
+          write(json_unit, fmtk(i)) extrapolated_energy(i,k)
+      enddo
+      write(json_unit, json_array_close_fmtx)
+      if (k < N_states_p) then
+        write(json_unit, json_dict_close_fmt)
+      else
+        write(json_unit, json_dict_close_fmtx)
+      endif
+    enddo
+    write(json_unit, json_array_close_fmtx)
+    write(json_unit, json_dict_close_fmt)
+    deallocate(fmtk)
+    call unlock_io
+end

src/cipsi_tc_bi_ortho/NEED

@@ -1,6 +1,7 @@
+json
 mpi
 perturbation
 zmq
-iterations_tc
+iterations
 csf
 tc_bi_ortho

src/cipsi_tc_bi_ortho/cipsi.irp.f

@@ -64,7 +64,7 @@ subroutine run_cipsi
   endif
 
   if (N_det > N_det_max) then
-    psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det)
+    psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det)
     psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states)
     N_det = N_det_max
     soft_touch N_det psi_det psi_coef

src/cipsi_tc_bi_ortho/pt2.irp.f

@@ -52,7 +52,7 @@ subroutine pt2_tc_bi_ortho
 !  call routine_save_right
 
   if (N_det > N_det_max) then
-    psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det)
+    psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det)
     psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states)
     N_det = N_det_max
     soft_touch N_det psi_det psi_coef

src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f

@@ -134,7 +134,6 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
   PROVIDE psi_det_hii selection_weight pseudo_sym
   PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max
   PROVIDE excitation_beta_max  excitation_alpha_max excitation_max
-  PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
 
   if (h0_type == 'CFG') then
     PROVIDE psi_configuration_hii det_to_configuration

src/cipsi_tc_bi_ortho/selection.irp.f

@@ -181,7 +181,6 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
   PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order
   PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
   PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp_tc
-  PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
 
   PROVIDE banned_excitation
 
@@ -464,15 +463,15 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
 
       do i = 1, fullinteresting(0)
         do k = 1, N_int
-          fullminilist(k,1,i) = psi_det_sorted_tc(k,1,fullinteresting(i))
+          fullminilist(k,1,i) = psi_selectors(k,1,fullinteresting(i))
-          fullminilist(k,2,i) = psi_det_sorted_tc(k,2,fullinteresting(i))
+          fullminilist(k,2,i) = psi_selectors(k,2,fullinteresting(i))
         enddo
       enddo
 
       do i = 1, interesting(0)
         do k = 1, N_int
-          minilist(k,1,i) = psi_det_sorted_tc(k,1,interesting(i))
+          minilist(k,1,i) = psi_selectors(k,1,interesting(i))
-          minilist(k,2,i) = psi_det_sorted_tc(k,2,interesting(i))
+          minilist(k,2,i) = psi_selectors(k,2,interesting(i))
         enddo
       enddo
 
@@ -616,7 +615,6 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
 
 
   PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc
-  PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
 
 
   mat   = 0d0
@@ -628,7 +626,10 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
     negMask(i,2) = not(mask(i,2))
   end do
 
+!  print*,'in selection '
   do i = 1, N_sel
+!    call debug_det(det(1,1,i),N_int)
+!    print*,i,dabs(psi_selectors_coef_transp_tc(1,2,i) * psi_selectors_coef_transp_tc(1,1,i))
     if(interesting(i) < 0) then
       stop 'prefetch interesting(i) and det(i)'
     endif
@@ -674,9 +675,6 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
           perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
           perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
         end do
-!        call get_d3_h  ( det(1,1,i), bannedOrb, banned, mat         , mask, p, sp, psi_selectors_coef_transp_tc (1, interesting(i)) )
-!        call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_r, mat_l, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) &
-!                       , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) )
 
         call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
         call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int)
@@ -916,8 +914,29 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
            psi_h_alpha = mat_l(istate, p1, p2)
            alpha_h_psi = mat_r(istate, p1, p2)
           endif
+          val = 4.d0 * psi_h_alpha * alpha_h_psi 
+          tmp = dsqrt(delta_E * delta_E + val)
+!          if (delta_E < 0.d0) then
+!              tmp = -tmp
+!          endif
+          e_pert(istate) = 0.25 * val / delta_E
+!          e_pert(istate) = 0.5d0 * (tmp - delta_E)
+          if(dsqrt(dabs(tmp)).gt.1.d-4.and.dabs(alpha_h_psi).gt.1.d-4)then
+           coef(istate)   = e_pert(istate) / psi_h_alpha 
+          else
            coef(istate)   = alpha_h_psi / delta_E 
-          e_pert(istate) = coef(istate) * psi_h_alpha
+          endif
+
+          if(selection_tc == 1)then
+           if(e_pert(istate).lt.0.d0)then 
+            e_pert(istate)=0.d0
+           else 
+            e_pert(istate)=-e_pert(istate)
+           endif
+          else if(selection_tc == -1)then
+           if(e_pert(istate).gt.0.d0)e_pert(istate)=0.d0
+          endif
+          
 !         if(selection_tc     ==  1 )then
 !          if(e_pert(istate).lt.0.d0)then
 !           e_pert(istate) = 0.d0
@@ -930,7 +949,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         enddo
 
 
-      do_diag = sum(dabs(coef)) > 0.001d0 .and. N_states > 1
 
       do istate = 1, N_states
 

src/cipsi_tc_bi_ortho/selection_buffer.irp.f

@@ -125,7 +125,11 @@ subroutine merge_selection_buffers(b1, b2)
   enddo
   b2%det => detmp
   b2%val => val
+!  if(selection_tc == 1)then
+!   b2%mini = max(b2%mini,b2%val(b2%N))
+!  else
    b2%mini = min(b2%mini,b2%val(b2%N))
+!  endif
   b2%cur = nmwen
 end
 
@@ -157,7 +161,11 @@ subroutine sort_selection_buffer(b)
   end do
   deallocate(b%det,iorder)
   b%det => detmp
+!  if(selection_tc == 1)then
+!   b%mini = max(b%mini,b%val(b%N))
+!  else
    b%mini = min(b%mini,b%val(b%N))
+!  endif
   b%cur = nmwen
 end subroutine
 

src/cipsi_tc_bi_ortho/slave_cipsi.irp.f

@@ -17,7 +17,6 @@ end
 
 subroutine provide_everything
   PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag
-  PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
 
   PROVIDE pt2_e0_denominator mo_num N_int ci_energy mpi_master zmq_state zmq_context
   PROVIDE psi_det psi_coef threshold_generators state_average_weight
@@ -312,7 +311,6 @@ subroutine run_slave_main
           PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order
           PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
           PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted_tc
-          PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
 
           PROVIDE psi_det_hii selection_weight pseudo_sym pt2_min_parallel_tasks
 

src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f

@@ -54,7 +54,7 @@ subroutine run_stochastic_cipsi
 
 
 !  if (N_det > N_det_max) then
-!    psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det)
+!    psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det)
 !    psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states)
 !    N_det = N_det_max
 !    soft_touch N_det psi_det psi_coef
@@ -78,6 +78,8 @@ subroutine run_stochastic_cipsi
         (N_det < N_det_max) .and.                                    &
         (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max)          &
         )
+      print*,'maxval(abs(pt2_data % pt2(1:N_states)))',maxval(abs(pt2_data % pt2(1:N_states)))
+      print*,pt2_max
       write(*,'(A)')  '--------------------------------------------------------------------------------'
 
 
@@ -92,7 +94,15 @@ subroutine run_stochastic_cipsi
     call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection
 !    stop
 
-    N_iter += 1
+    call print_summary(psi_energy_with_nucl_rep, &
+       pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2)
+
+    call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
+
+    call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
+    call print_extrapolated_energy()
+!    call print_mol_properties()
+    call write_cipsi_json(pt2_data,pt2_data_err)
 
     if (qp_stop()) exit
 
@@ -106,6 +116,7 @@ subroutine run_stochastic_cipsi
     ept2(N_iter-1) = E_tc + nuclear_repulsion + (pt2_data % pt2(1))/norm
     pt1(N_iter-1) = dsqrt(pt2_data % overlap(1,1))
     call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
+!    stop
     if (qp_stop()) exit
   enddo
 !  print*,'data to extrapolate '

src/cipsi_tc_bi_ortho/write_cipsi_json.irp.f

@@ -0,0 +1,53 @@
+subroutine write_cipsi_json(pt2_data, pt2_data_err)
+    use selection_types
+    implicit none
+    BEGIN_DOC
+!  Writes JSON data for CIPSI runs
+    END_DOC
+    type(pt2_type), intent(in) :: pt2_data, pt2_data_err
+    integer :: i,j,k
+
+    call lock_io
+    character*(64), allocatable    :: fmtk(:)
+    integer :: N_states_p, N_iter_p
+    N_states_p = min(N_states,N_det)
+    N_iter_p = min(N_iter,8)
+    allocate(fmtk(0:N_iter_p))
+    fmtk(:) = '(''   '',E22.15,'','')'
+    fmtk(N_iter_p) = '(''   '',E22.15)'
+
+    write(json_unit, json_dict_uopen_fmt)
+    write(json_unit, json_int_fmt) 'n_det', N_det
+    if (s2_eig) then
+      write(json_unit, json_int_fmt) 'n_cfg', N_configuration
+      if (only_expected_s2) then
+        write(json_unit, json_int_fmt) 'n_csf', N_csf
+      endif
+    endif
+    write(json_unit, json_array_open_fmt) 'states'
+    do k=1,N_states_p
+      write(json_unit, json_dict_uopen_fmt)
+      write(json_unit, json_real_fmt) 'energy', psi_energy_with_nucl_rep(k)
+      write(json_unit, json_real_fmt) 's2', psi_s2(k)
+      write(json_unit, json_real_fmt) 'pt2', pt2_data % pt2(k)
+      write(json_unit, json_real_fmt) 'pt2_err', pt2_data_err % pt2(k)
+      write(json_unit, json_real_fmt) 'rpt2', pt2_data % rpt2(k)
+      write(json_unit, json_real_fmt) 'rpt2_err', pt2_data_err % rpt2(k)
+      write(json_unit, json_real_fmt) 'variance', pt2_data % variance(k)
+      write(json_unit, json_real_fmt) 'variance_err', pt2_data_err % variance(k)
+      write(json_unit, json_array_open_fmt) 'ex_energy'
+      do i=2,N_iter_p
+          write(json_unit, fmtk(i)) extrapolated_energy(i,k)
+      enddo
+      write(json_unit, json_array_close_fmtx)
+      if (k < N_states_p) then
+        write(json_unit, json_dict_close_fmt)
+      else
+        write(json_unit, json_dict_close_fmtx)
+      endif
+    enddo
+    write(json_unit, json_array_close_fmtx)
+    write(json_unit, json_dict_close_fmt)
+    deallocate(fmtk)
+    call unlock_io
+end

src/cisd/30.cisd.bats

@@ -10,6 +10,7 @@ function run() {
   qp set determinants n_states  2
   qp set davidson threshold_davidson 1.e-12
   qp set davidson n_states_diag 24
+  qp run cis
   qp run cisd 
   energy1="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 1)"
   energy2="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 2)"
@@ -20,26 +21,31 @@ function run() {
 
 @test "B-B" { # 
   qp set_file b2_stretched.ezfio
+  qp set_frozen_core
   run  -49.120607088648597 -49.055152453388231
 }
 
 @test "SiH2_3B1" { # 1.53842s 3.53856s
   qp set_file sih2_3b1.ezfio
+  qp set_frozen_core
   run -290.015949171697 -289.805036176618
 }
 
 @test "HBO" { # 4.42968s  19.6099s
   qp set_file hbo.ezfio
+  qp set_frozen_core
   run -100.2019254455993 -99.79484127741013 
 }
 
 @test "HCO" { # 6.6077s 28.6801s
   qp set_file hco.ezfio
+  qp set_frozen_core
   run -113.39088802205114 -113.22204293108558
 }
 
 @test "H2O" { # 7.0651s 30.6642s
   qp set_file h2o.ezfio
+  qp set_frozen_core
   run -76.22975602077072 -75.80609108747208 
 }
 
@@ -50,6 +56,7 @@ function run() {
 @test "H2S" { # 7.42152s 32.5461s
   [[ -n $TRAVIS ]] && skip
   qp set_file h2s.ezfio
+  qp set_frozen_core
   run -398.853701416768 -398.519020035337
 }
 
@@ -70,6 +77,7 @@ function run() {
 @test "OH" { # 18.2159s 1.28453m
   [[ -n $TRAVIS ]] && skip
   qp set_file oh.ezfio
+  qp set_frozen_core
   run -75.6087472926588 -75.5370393736601
 }
 
@@ -83,6 +91,7 @@ function run() {
 @test "SiH3" { # 20.2202s 1.38648m
   [[ -n $TRAVIS ]] && skip
   qp set_file sih3.ezfio
+  qp set_frozen_core
   run -5.57096611856522  -5.30950347928823
 }
 
@@ -103,6 +112,7 @@ function run() {
 @test "H3COH" { # 24.7248s 1.85043m
   [[ -n $TRAVIS ]] && skip
   qp set_file h3coh.ezfio
+  qp set_frozen_core
   run -115.204958752377 -114.755913828245
 }
 
@@ -117,6 +127,7 @@ function run() {
 @test "ClF" { # 30.3225s
   [[ -n $TRAVIS ]] && skip
   qp set_file clf.ezfio
+  qp set_frozen_core
   run -559.162476603880  -558.792395927088
 }
 
@@ -130,6 +141,7 @@ function run() {
 @test "ClO" { # 37.6949s
   [[ -n $TRAVIS ]] && skip
   qp set_file clo.ezfio
+  qp set_frozen_core
   run -534.5404021326773 -534.3818725793897 
 }
 
@@ -150,6 +162,7 @@ function run() {
 @test "SO" { # 51.2476s
   [[ -n $TRAVIS ]] && skip
   qp set_file so.ezfio
+  qp set_frozen_core
   run -26.0131812819785 -25.7053111980226
 }
 

src/cisd/cisd.irp.f

@@ -69,7 +69,9 @@ subroutine run
   do i = 1,N_states
     k = maxloc(dabs(psi_coef_sorted(1:N_det,i)),dim=1)
     delta_E  = CI_electronic_energy(i) - diag_h_mat_elem(psi_det_sorted(1,1,k),N_int)
+    if (elec_alpha_num + elec_beta_num >= 4) then
       cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2)
+    endif
   enddo
   print *,  'N_det = ', N_det
   print*,''
@@ -78,13 +80,16 @@ subroutine run
   do i = 1,N_states
     print *,  i, CI_energy(i)
   enddo
+  if (elec_alpha_num + elec_beta_num >= 4) then
     print*,''
     print*,'******************************'
     print *,  'CISD+Q Energies'
     do i = 1,N_states
       print *,  i, cisdq(i)
     enddo
+  endif
   if (N_states > 1) then
+    if (elec_alpha_num + elec_beta_num >= 4) then
       print*,''
       print*,'******************************'
       print*,'Excitation energies (au)    (CISD+Q)'
@@ -95,9 +100,23 @@ subroutine run
       print*,'******************************'
       print*,'Excitation energies (eV)    (CISD+Q)'
       do i = 2, N_states
-      print*, i ,(CI_energy(i) - CI_energy(1))/0.0367502d0, &
+        print*, i ,(CI_energy(i) - CI_energy(1)) * ha_to_ev, &
-        (cisdq(i) - cisdq(1)) / 0.0367502d0
+          (cisdq(i) - cisdq(1)) * ha_to_ev
       enddo
+    else
+      print*,''
+      print*,'******************************'
+      print*,'Excitation energies (au)    (CISD)'
+      do i = 2, N_states
+        print*, i ,CI_energy(i) - CI_energy(1)
+      enddo
+      print*,''
+      print*,'******************************'
+      print*,'Excitation energies (eV)    (CISD)'
+      do i = 2, N_states
+        print*, i ,(CI_energy(i) - CI_energy(1)) * ha_to_ev
+      enddo
+    endif
   endif
 
 end

src/davidson/davidson_parallel.irp.f

@@ -150,7 +150,9 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, N_st, sze,
       exit
     endif
     if(task_id == 0) exit
+    call lock_io()
     read (msg,*) imin, imax, ishift, istep
+    call unlock_io()
     integer :: k
     do k=imin,imax
       v_t(:,k) = 0.d0
@@ -555,7 +557,9 @@ BEGIN_PROVIDER [ integer, nthreads_davidson ]
  character*(32) :: env
  call getenv('QP_NTHREADS_DAVIDSON',env)
  if (trim(env) /= '') then
+   call lock_io()
    read(env,*) nthreads_davidson
+   call unlock_io()
    call write_int(6,nthreads_davidson,'Target number of threads for <Psi|H|Psi>')
  endif
 END_PROVIDER

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -466,6 +466,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
        double precision, allocatable :: work(:)
 
        y = h
+!       y = h_p   ! Doesn't work for non-singlets
        lwork = -1
        allocate(work(1))
        call dsygv(1,'V','U',shift2,y,size(y,1), &

src/davidson/input.irp.f

@@ -30,10 +30,5 @@ BEGIN_PROVIDER [ integer, n_states_diag  ]
     endif
   IRP_ENDIF
 
-  call write_time(6)
-  if (mpi_master) then
-    write(6, *) 'Read  n_states_diag'
-  endif
-
 END_PROVIDER
 

src/determinants/h_apply.irp.f

@@ -69,9 +69,15 @@ subroutine resize_H_apply_buffer(new_size,iproc)
   END_DOC
   PROVIDE H_apply_buffer_allocated
 
+
   ASSERT (new_size > 0)
   ASSERT (iproc >= 0)
   ASSERT (iproc < nproc)
+  if (N_det < 0) call abort() !irp_here//': N_det < 0')
+  if (N_int <= 0) call abort() !irp_here//': N_int <= 0')
+  if (new_size <= 0) call abort() !irp_here//': new_size <= 0')
+  if (iproc < 0) call abort() !irp_here//': iproc < 0')
+  if (iproc >= nproc) call abort() !irp_here//': iproc >= nproc')
 
   allocate ( buffer_det(N_int,2,new_size),                           &
       buffer_coef(new_size,N_states),                                &
@@ -126,8 +132,10 @@ subroutine copy_H_apply_buffer_to_wf
 
 
   ASSERT (N_int > 0)
-  ASSERT (N_det > 0)
+  ASSERT (N_det >= 0)
 
+  N_det_old = N_det
+  if (N_det > 0) then
      allocate ( buffer_det(N_int,2,N_det), buffer_coef(N_det,N_states) )
 
      ! Backup determinants
@@ -151,6 +159,7 @@ subroutine copy_H_apply_buffer_to_wf
        enddo
        ASSERT ( j == N_det_old )
      enddo
+  endif
 
   ! Update N_det
   N_det = N_det_old
@@ -164,6 +173,7 @@ subroutine copy_H_apply_buffer_to_wf
     TOUCH psi_det_size
   endif
 
+  if (N_det_old > 0) then
     ! Restore backup in resized array
     do i=1,N_det_old
       psi_det(:,:,i) = buffer_det(:,:,i)
@@ -175,6 +185,7 @@ subroutine copy_H_apply_buffer_to_wf
         psi_coef(i,k) = buffer_coef(i,k)
       enddo
     enddo
+  endif
 
   ! Copy new buffers
 
@@ -339,3 +350,33 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
   call omp_unset_lock(H_apply_buffer_lock(1,iproc))
 end
 
+
+subroutine replace_wf(N_det_new, LDA, psi_coef_new, psi_det_new)
+  use omp_lib
+  implicit none
+  BEGIN_DOC
+! Replaces the wave function.
+! After calling this subroutine, N_det, psi_det and psi_coef need to be touched
+  END_DOC
+  integer, intent(in)            :: N_det_new, LDA
+  double precision, intent(in)   :: psi_coef_new(LDA,N_states)
+  integer(bit_kind), intent(in)  :: psi_det_new(N_int,2,N_det_new)
+
+  integer                        :: i,j
+
+  PROVIDE H_apply_buffer_allocated
+
+  if (N_det_new <= 0) call abort() !irp_here//': N_det_new <= 0')
+  if (N_int <= 0) call abort() !irp_here//': N_int <= 0')
+  if (LDA < N_det_new) call abort() !irp_here//': LDA < N_det_new')
+
+  do j=0,nproc-1
+    H_apply_buffer(j)%N_det = 0
+  enddo
+  N_det = 0
+  SOFT_TOUCH N_det
+  call fill_H_apply_buffer_no_selection(N_det_new,psi_det_new,N_int,0)
+  call copy_h_apply_buffer_to_wf
+  psi_coef(1:N_det_new,1:N_states) = psi_coef_new(1:N_det_new,1:N_states)
+
+end

src/determinants/tr_density_matrix.irp.f

@@ -0,0 +1,313 @@
+BEGIN_PROVIDER [double precision, one_e_tr_dm_mo, (mo_num, mo_num, N_states, N_states)]
+
+  implicit none
+
+  BEGIN_DOC
+  ! One body transition density matrix for all pairs of states n and m, < Psi^n | a_i^\dagger a_a | Psi^m >
+  END_DOC
+
+  integer                        :: j,k,l,m,k_a,k_b,n
+  integer                        :: occ(N_int*bit_kind_size,2)
+  double precision               :: ck, cl, ckl
+  double precision               :: phase
+  integer                        :: h1,h2,p1,p2,s1,s2, degree
+  integer(bit_kind)              :: tmp_det(N_int,2), tmp_det2(N_int)
+  integer                        :: exc(0:2,2),n_occ(2)
+  double precision, allocatable  :: tmp_a(:,:,:,:), tmp_b(:,:,:,:)
+  integer                        :: krow, kcol, lrow, lcol
+
+  PROVIDE psi_det
+
+  one_e_tr_dm_mo = 0d0
+ 
+  !$OMP PARALLEL DEFAULT(NONE)                                      &
+      !$OMP PRIVATE(j,k,k_a,k_b,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc,&
+      !$OMP  tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
+      !$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,  &
+      !$OMP  elec_beta_num,one_e_tr_dm_mo,N_det,&
+      !$OMP  mo_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns,&
+      !$OMP  psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns,&
+      !$OMP  psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique,&
+      !$OMP  psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values,&
+      !$OMP  N_det_alpha_unique,N_det_beta_unique,irp_here)
+  allocate(tmp_a(mo_num,mo_num,N_states,N_states), tmp_b(mo_num,mo_num,N_states,N_states) )
+  tmp_a = 0.d0
+  !$OMP DO SCHEDULE(dynamic,64)
+  do k_a=1,N_det
+    krow = psi_bilinear_matrix_rows(k_a)
+    ASSERT (krow <= N_det_alpha_unique)
+
+    kcol = psi_bilinear_matrix_columns(k_a)
+    ASSERT (kcol <= N_det_beta_unique)
+
+    tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
+    tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol)
+
+    ! Diagonal part
+    ! -------------
+
+    call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
+    do m=1,N_states
+      do n = 1, N_states
+        ck = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(k_a,n)
+        do l=1,elec_alpha_num
+          j = occ(l,1)
+          tmp_a(j,j,m,n) += ck
+        enddo
+      enddo
+    enddo
+
+    if (k_a == N_det) cycle
+    l = k_a+1
+    lrow = psi_bilinear_matrix_rows(l)
+    lcol = psi_bilinear_matrix_columns(l)
+    ! Fix beta determinant, loop over alphas
+    do while ( lcol == kcol )
+      tmp_det2(:) = psi_det_alpha_unique(:, lrow)
+      call get_excitation_degree_spin(tmp_det(1,1),tmp_det2,degree,N_int)
+      if (degree == 1) then
+        exc = 0
+        call get_single_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int)
+        call decode_exc_spin(exc,h1,p1,h2,p2)
+        do m=1,N_states
+          do n = 1, N_states
+            ckl = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(l,n) * phase
+            tmp_a(h1,p1,m,n) += ckl
+            ckl = psi_bilinear_matrix_values(k_a,n)*psi_bilinear_matrix_values(l,m) * phase
+            tmp_a(p1,h1,m,n) += ckl
+          enddo
+        enddo
+      endif
+      l = l+1
+      if (l>N_det) exit
+      lrow = psi_bilinear_matrix_rows(l)
+      lcol = psi_bilinear_matrix_columns(l)
+    enddo
+
+  enddo
+  !$OMP END DO NOWAIT
+
+  !$OMP CRITICAL
+  one_e_tr_dm_mo(:,:,:,:) = one_e_tr_dm_mo(:,:,:,:) + tmp_a(:,:,:,:)
+  !$OMP END CRITICAL
+  deallocate(tmp_a)
+  !$OMP BARRIER
+
+  tmp_b = 0.d0
+  !$OMP DO SCHEDULE(dynamic,64)
+  do k_b=1,N_det
+    krow = psi_bilinear_matrix_transp_rows(k_b)
+    ASSERT (krow <= N_det_alpha_unique)
+
+    kcol = psi_bilinear_matrix_transp_columns(k_b)
+    ASSERT (kcol <= N_det_beta_unique)
+
+    tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
+    tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol)
+
+    ! Diagonal part
+    ! -------------
+
+    call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
+    do m=1,N_states
+      do n = 1, N_states
+        ck = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(k_b,n)
+        do l=1,elec_beta_num
+          j = occ(l,2)
+          tmp_b(j,j,m,n) += ck
+        enddo
+      enddo
+    enddo
+
+    if (k_b == N_det) cycle
+    l = k_b+1
+    lrow = psi_bilinear_matrix_transp_rows(l)
+    lcol = psi_bilinear_matrix_transp_columns(l)
+    ! Fix beta determinant, loop over alphas
+    do while ( lrow == krow )
+      tmp_det2(:) = psi_det_beta_unique(:, lcol)
+      call get_excitation_degree_spin(tmp_det(1,2),tmp_det2,degree,N_int)
+      if (degree == 1) then
+        exc = 0
+        call get_single_excitation_spin(tmp_det(1,2),tmp_det2,exc,phase,N_int)
+        call decode_exc_spin(exc,h1,p1,h2,p2)
+        do m=1,N_states
+          do n = 1, N_states
+            ckl = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(l,n) * phase
+            tmp_b(h1,p1,m,n) += ckl
+            ckl = psi_bilinear_matrix_transp_values(k_b,n)*psi_bilinear_matrix_transp_values(l,m) * phase
+            tmp_b(p1,h1,m,n) += ckl
+          enddo
+        enddo
+      endif
+      l = l+1
+      if (l>N_det) exit
+      lrow = psi_bilinear_matrix_transp_rows(l)
+      lcol = psi_bilinear_matrix_transp_columns(l)
+    enddo
+
+  enddo
+  !$OMP END DO NOWAIT
+  !$OMP CRITICAL
+  one_e_tr_dm_mo(:,:,:,:)  = one_e_tr_dm_mo(:,:,:,:)  + tmp_b(:,:,:,:)
+  !$OMP END CRITICAL
+
+  deallocate(tmp_b)
+  !$OMP END PARALLEL
+
+END_PROVIDER
+ BEGIN_PROVIDER [ double precision, one_e_tr_dm_mo_alpha, (mo_num,mo_num,N_states,N_states) ]
+&BEGIN_PROVIDER [ double precision, one_e_tr_dm_mo_beta, (mo_num,mo_num,N_states,N_states) ]
+  implicit none
+  BEGIN_DOC
+  ! $\alpha$ and $\beta$ one-body transition density matrices for all pairs of states
+  END_DOC
+
+  integer                        :: j,k,l,m,n,k_a,k_b
+  integer                        :: occ(N_int*bit_kind_size,2)
+  double precision               :: ck, cl, ckl
+  double precision               :: phase
+  integer                        :: h1,h2,p1,p2,s1,s2, degree
+  integer(bit_kind)              :: tmp_det(N_int,2), tmp_det2(N_int)
+  integer                        :: exc(0:2,2),n_occ(2)
+  double precision, allocatable  :: tmp_a(:,:,:,:), tmp_b(:,:,:,:)
+  integer                        :: krow, kcol, lrow, lcol
+
+  PROVIDE psi_det
+
+  one_e_tr_dm_mo_alpha = 0.d0
+  one_e_tr_dm_mo_beta  = 0.d0
+  !$OMP PARALLEL DEFAULT(NONE)                                      &
+      !$OMP PRIVATE(j,k,k_a,k_b,l,m,n,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc,&
+      !$OMP  tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
+      !$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,  &
+      !$OMP  elec_beta_num,one_e_tr_dm_mo_alpha,one_e_tr_dm_mo_beta,N_det,&
+      !$OMP  mo_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns,&
+      !$OMP  psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns,&
+      !$OMP  psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique,&
+      !$OMP  psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values,&
+      !$OMP  N_det_alpha_unique,N_det_beta_unique,irp_here)
+  allocate(tmp_a(mo_num,mo_num,N_states,N_states), tmp_b(mo_num,mo_num,N_states,N_states) )
+  tmp_a = 0.d0
+  !$OMP DO SCHEDULE(dynamic,64)
+  do k_a=1,N_det
+    krow = psi_bilinear_matrix_rows(k_a)
+    ASSERT (krow <= N_det_alpha_unique)
+
+    kcol = psi_bilinear_matrix_columns(k_a)
+    ASSERT (kcol <= N_det_beta_unique)
+
+    tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
+    tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol)
+
+    ! Diagonal part
+    ! -------------
+
+    call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
+    do m=1,N_states
+      do n = 1, N_states
+        ck = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(k_a,n)
+        do l=1,elec_alpha_num
+          j = occ(l,1)
+          tmp_a(j,j,m,n) += ck
+        enddo
+      enddo
+    enddo
+
+    if (k_a == N_det) cycle
+    l = k_a+1
+    lrow = psi_bilinear_matrix_rows(l)
+    lcol = psi_bilinear_matrix_columns(l)
+    ! Fix beta determinant, loop over alphas
+    do while ( lcol == kcol )
+      tmp_det2(:) = psi_det_alpha_unique(:, lrow)
+      call get_excitation_degree_spin(tmp_det(1,1),tmp_det2,degree,N_int)
+      if (degree == 1) then
+        exc = 0
+        call get_single_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int)
+        call decode_exc_spin(exc,h1,p1,h2,p2)
+        do m=1,N_states
+          do n = 1, N_states
+            ckl = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(l,n) * phase
+            tmp_a(h1,p1,m,n) += ckl
+            tmp_a(p1,h1,m,n) += ckl
+          enddo
+        enddo
+      endif
+      l = l+1
+      if (l>N_det) exit
+      lrow = psi_bilinear_matrix_rows(l)
+      lcol = psi_bilinear_matrix_columns(l)
+    enddo
+
+  enddo
+  !$OMP END DO NOWAIT
+
+  !$OMP CRITICAL
+  one_e_tr_dm_mo_alpha(:,:,:,:) = one_e_tr_dm_mo_alpha(:,:,:,:) + tmp_a(:,:,:,:)
+  !$OMP END CRITICAL
+  deallocate(tmp_a)
+
+  tmp_b = 0.d0
+  !$OMP DO SCHEDULE(dynamic,64)
+  do k_b=1,N_det
+    krow = psi_bilinear_matrix_transp_rows(k_b)
+    ASSERT (krow <= N_det_alpha_unique)
+
+    kcol = psi_bilinear_matrix_transp_columns(k_b)
+    ASSERT (kcol <= N_det_beta_unique)
+
+    tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
+    tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol)
+
+    ! Diagonal part
+    ! -------------
+
+    call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
+    do m=1,N_states
+      do n = 1, N_states
+        ck = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(k_b,n)
+        do l=1,elec_beta_num
+          j = occ(l,2)
+          tmp_b(j,j,m,n) += ck
+        enddo
+      enddo
+    enddo
+
+    if (k_b == N_det) cycle
+    l = k_b+1
+    lrow = psi_bilinear_matrix_transp_rows(l)
+    lcol = psi_bilinear_matrix_transp_columns(l)
+    ! Fix beta determinant, loop over alphas
+    do while ( lrow == krow )
+      tmp_det2(:) = psi_det_beta_unique(:, lcol)
+      call get_excitation_degree_spin(tmp_det(1,2),tmp_det2,degree,N_int)
+      if (degree == 1) then
+        exc = 0
+        call get_single_excitation_spin(tmp_det(1,2),tmp_det2,exc,phase,N_int)
+        call decode_exc_spin(exc,h1,p1,h2,p2)
+        do m=1,N_states
+          do n = 1, N_states
+            ckl = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(l,n) * phase
+            tmp_b(h1,p1,m,n) += ckl
+            tmp_b(p1,h1,m,n) += ckl
+          enddo
+        enddo
+      endif
+      l = l+1
+      if (l>N_det) exit
+      lrow = psi_bilinear_matrix_transp_rows(l)
+      lcol = psi_bilinear_matrix_transp_columns(l)
+    enddo
+
+  enddo
+  !$OMP END DO NOWAIT
+  !$OMP CRITICAL
+  one_e_tr_dm_mo_beta(:,:,:,:)  = one_e_tr_dm_mo_beta(:,:,:,:)  + tmp_b(:,:,:,:)
+  !$OMP END CRITICAL
+
+  deallocate(tmp_b)
+  !$OMP END PARALLEL
+
+END_PROVIDER
+

src/ezfio_files/lock.irp.f

@@ -9,4 +9,21 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), file_lock ]
   call omp_init_lock(file_lock)
 END_PROVIDER
 
+! These functions need to be called because internal read and write are not thread safe.
+subroutine lock_io()
+  implicit none
+  BEGIN_DOC
+! Needs to be called because before doing I/O because internal read and write
+! are not thread safe.
+  END_DOC
+  call omp_set_lock(file_lock)
+end subroutine lock_io()
 
+subroutine unlock_io()
+  implicit none
+  BEGIN_DOC
+! Needs to be called because afterdoing I/O because internal read and write
+! are not thread safe.
+  END_DOC
+  call omp_unset_lock(file_lock)
+end subroutine lock_io()

src/fci/fci.irp.f

@@ -39,12 +39,19 @@ program fci
   if (.not.is_zmq_slave) then
     PROVIDE psi_det psi_coef mo_two_e_integrals_in_map
 
+    write(json_unit,json_array_open_fmt) 'fci'
+
     if (do_pt2) then
       call run_stochastic_cipsi
     else
       call run_cipsi
     endif
 
+    write(json_unit,json_dict_uopen_fmt)
+    write(json_unit,json_dict_close_fmtx)
+    write(json_unit,json_array_close_fmtx)
+    call json_close
+
   else
     PROVIDE mo_two_e_integrals_in_map pt2_min_parallel_tasks
 

src/fci_tc_bi/NEED

@@ -1,3 +1,4 @@
+json
 tc_bi_ortho
 davidson_undressed
 cipsi_tc_bi_ortho

src/fci_tc_bi/diagonalize_ci.irp.f

@@ -49,9 +49,8 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
       psi_coef(i,j) = dabs(psi_l_coef_bi_ortho(i,j) * psi_r_coef_bi_ortho(i,j))   
     enddo
   enddo
-  SOFT_TOUCH  eigval_left_tc_bi_orth  eigval_right_tc_bi_orth  leigvec_tc_bi_orth  reigvec_tc_bi_orth norm_ground_left_right_bi_orth psi_coef  psi_l_coef_bi_ortho psi_r_coef_bi_ortho
+  SOFT_TOUCH eigval_left_tc_bi_orth  eigval_right_tc_bi_orth leigvec_tc_bi_orth reigvec_tc_bi_orth norm_ground_left_right_bi_orth 
-
+  SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho psi_coef  
-
 
   call save_tc_bi_ortho_wavefunction
 end

src/fci_tc_bi/fci_tc_bi_ortho.irp.f

@@ -62,6 +62,7 @@ subroutine run_cipsi_tc
       endif
     endif
     ! ---
+    write(json_unit,json_array_open_fmt) 'fci_tc'
 
     if (do_pt2) then
       call run_stochastic_cipsi
@@ -69,6 +70,11 @@ subroutine run_cipsi_tc
       call run_cipsi
     endif
 
+    write(json_unit,json_dict_uopen_fmt)
+    write(json_unit,json_dict_close_fmtx)
+    write(json_unit,json_array_close_fmtx)
+    call json_close
+
   else
     PROVIDE mo_bi_ortho_tc_one_e mo_bi_ortho_tc_two_e pt2_min_parallel_tasks
     if(elec_alpha_num+elec_beta_num.ge.3)then

src/fci_tc_bi/generators.irp.f

@@ -46,6 +46,11 @@ END_PROVIDER
  psi_det_sorted_tc_gen =  psi_det_sorted_tc
  psi_coef_sorted_tc_gen = psi_coef_sorted_tc
  psi_det_sorted_tc_gen_order = psi_det_sorted_tc_order
+ integer :: i
+! do i = 1,N_det
+!  print*,'i = ',i
+!  call debug_det(psi_det_sorted_tc(1,1,i),N_int)
+! enddo
 END_PROVIDER
 
 

src/fci_tc_bi/save_energy.irp.f

@@ -4,6 +4,6 @@ subroutine save_energy(E,pt2)
 ! Saves the energy in |EZFIO|.
   END_DOC
   double precision, intent(in) :: E(N_states), pt2(N_states)
-  call ezfio_set_fci_tc_energy(E(1:N_states))
+  call ezfio_set_fci_tc_bi_energy(E(1:N_states))
-  call ezfio_set_fci_tc_energy_pt2(E(1:N_states)+pt2(1:N_states))
+  call ezfio_set_fci_tc_bi_energy_pt2(E(1:N_states)+pt2(1:N_states))
 end

src/fci_tc_bi/selectors.irp.f

@@ -18,15 +18,6 @@ BEGIN_PROVIDER [ integer, N_det_selectors]
   double precision               :: norm, norm_max
   call write_time(6)
   N_det_selectors = N_det
-  norm = 1.d0
-  do i=1,N_det
-    norm = norm - psi_average_norm_contrib_tc(i)
-    if (norm - 1.d-10 < 1.d0 - threshold_selectors) then
-      N_det_selectors = i
-      exit
-    endif
-  enddo
-  N_det_selectors = max(N_det_selectors,N_det_generators)
   call write_int(6,N_det_selectors,'Number of selectors')
 END_PROVIDER
 
@@ -50,8 +41,6 @@ END_PROVIDER
       psi_selectors_coef(i,k)      = psi_coef_sorted_tc_gen(i,k)
       psi_selectors_coef_tc(i,1,k) = psi_l_coef_sorted_bi_ortho(i,k)
       psi_selectors_coef_tc(i,2,k) = psi_r_coef_sorted_bi_ortho(i,k)
-!      psi_selectors_coef_tc(i,1,k) = 1.d0
-!      psi_selectors_coef_tc(i,2,k) = 1.d0
     enddo
   enddo
 
@@ -74,25 +63,6 @@ END_PROVIDER
   enddo
 END_PROVIDER
 
- BEGIN_PROVIDER [ double precision, psi_selectors_rcoef_bi_orth_transp, (N_states, psi_det_size) ]
-&BEGIN_PROVIDER [ double precision, psi_selectors_lcoef_bi_orth_transp, (N_states, psi_det_size) ]
-
-  implicit none
-  integer :: i, k
-
-  psi_selectors_rcoef_bi_orth_transp = 0.d0
-  psi_selectors_lcoef_bi_orth_transp = 0.d0
-
-  print*,'N_det,N_det_selectors',N_det,N_det_selectors
-  do i = 1, N_det_selectors
-    do k = 1, N_states
-      psi_selectors_rcoef_bi_orth_transp(k,i) = psi_r_coef_sorted_bi_ortho(i,k)
-      psi_selectors_lcoef_bi_orth_transp(k,i) = psi_l_coef_sorted_bi_ortho(i,k)
-    enddo
-  enddo
-
-END_PROVIDER
-
 BEGIN_PROVIDER [ integer, psi_selectors_size ]
  implicit none
  psi_selectors_size = psi_det_size

src/hartree_fock/NEED

@@ -1,3 +1,4 @@
 ao_one_e_ints
 ao_two_e_ints
 scf_utils
+json

src/hartree_fock/fock_matrix_hf.irp.f

@@ -15,16 +15,22 @@
  double precision, allocatable  :: ao_two_e_integral_alpha_tmp(:,:)
  double precision, allocatable  :: ao_two_e_integral_beta_tmp(:,:)
 
+ if (do_ao_cholesky) then   ! Use Cholesky-decomposed integrals
+   ao_two_e_integral_alpha(:,:) = ao_two_e_integral_alpha_chol(:,:)
+   ao_two_e_integral_beta (:,:) = ao_two_e_integral_beta_chol (:,:)
+
+ else ! Use integrals in AO basis set
+
    ao_two_e_integral_alpha = 0.d0
    ao_two_e_integral_beta  = 0.d0
    if (do_direct_integrals) then
 
      !$OMP PARALLEL DEFAULT(NONE)                                    &
-       !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,keys,values,p,q,r,s,i0,j0,k0,l0, &
+         !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,keys,values,p,q,r,s,i0,j0,k0,l0,&
-       !$OMP ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, c0, c1, c2, &
+         !$OMP ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, c0, c1, c2,&
-       !$OMP local_threshold)&
+         !$OMP local_threshold)                                      &
          !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
-       !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz, &
+         !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz,&
          !$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta)
 
      allocate(keys(1), values(1))
@@ -105,7 +111,7 @@
      double precision, allocatable  :: values(:)
 
      !$OMP PARALLEL DEFAULT(NONE)                                    &
-       !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, &
+         !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max,&
          !$OMP  n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp)&
          !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
          !$OMP  ao_integrals_map, ao_two_e_integral_alpha, ao_two_e_integral_beta)
@@ -150,6 +156,81 @@
      !$OMP END PARALLEL
 
    endif
+ endif
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, ao_two_e_integral_alpha_chol, (ao_num, ao_num) ]
+&BEGIN_PROVIDER [ double precision, ao_two_e_integral_beta_chol ,  (ao_num, ao_num) ]
+ use map_module
+ implicit none
+ BEGIN_DOC
+ ! Alpha and Beta Fock matrices in AO basis set
+ END_DOC
+
+ integer :: m,n,l,s,j
+ double precision :: integral
+ double precision, allocatable :: X(:), X2(:,:,:,:), X3(:,:,:,:)
+
+ allocate (X(cholesky_ao_num))
+
+
+ ! X(j) = \sum_{mn} SCF_density_matrix_ao(m,n) * cholesky_ao(m,n,j)
+ call dgemm('T','N',cholesky_ao_num,1,ao_num*ao_num,1.d0,            &
+     cholesky_ao, ao_num*ao_num,                                     &
+     SCF_density_matrix_ao, ao_num*ao_num,0.d0,                      &
+     X, cholesky_ao_num)
+!
+
+ ! ao_two_e_integral_alpha(m,n) = \sum_{j} cholesky_ao(m,n,j) * X(j)
+ call dgemm('N','N',ao_num*ao_num,1,cholesky_ao_num, 1.d0,           &
+     cholesky_ao, ao_num*ao_num,                                     &
+     X, cholesky_ao_num, 0.d0,                                       &
+     ao_two_e_integral_alpha_chol, ao_num*ao_num)
+
+ deallocate(X)
+
+ ao_two_e_integral_beta_chol = ao_two_e_integral_alpha_chol
+
+
+ allocate(X2(ao_num,ao_num,cholesky_ao_num,2))
+
+! ao_two_e_integral_alpha_chol (l,s) -= cholesky_ao(l,m,j) * SCF_density_matrix_ao_beta (m,n) * cholesky_ao(n,s,j)
+
+ call dgemm('N','N',ao_num,ao_num*cholesky_ao_num,ao_num, 1.d0,      &
+     SCF_density_matrix_ao_alpha, ao_num,                            &
+     cholesky_ao, ao_num, 0.d0,                                      &
+     X2(1,1,1,1), ao_num)
+
+ call dgemm('N','N',ao_num,ao_num*cholesky_ao_num,ao_num, 1.d0,      &
+     SCF_density_matrix_ao_beta, ao_num,                             &
+     cholesky_ao, ao_num, 0.d0,                                      &
+     X2(1,1,1,2), ao_num)
+
+ allocate(X3(ao_num,cholesky_ao_num,ao_num,2))
+
+ do s=1,ao_num
+  do j=1,cholesky_ao_num
+   do m=1,ao_num
+    X3(m,j,s,1) = X2(m,s,j,1)
+    X3(m,j,s,2) = X2(m,s,j,2)
+   enddo
+  enddo
+ enddo
+
+ deallocate(X2)
+
+ call dgemm('N','N',ao_num,ao_num,ao_num*cholesky_ao_num, -1.d0,     &
+     cholesky_ao, ao_num,                                            &
+     X3(1,1,1,1), ao_num*cholesky_ao_num, 1.d0,                      &
+     ao_two_e_integral_alpha_chol, ao_num)
+
+ call dgemm('N','N',ao_num,ao_num,ao_num*cholesky_ao_num, -1.d0,     &
+     cholesky_ao, ao_num,                                            &
+     X3(1,1,1,2), ao_num*cholesky_ao_num, 1.d0,                      &
+     ao_two_e_integral_beta_chol, ao_num)
+
+ deallocate(X3)
 
 END_PROVIDER
 

src/hartree_fock/scf.irp.f

@@ -80,9 +80,14 @@ subroutine run
 
   mo_label = 'Orthonormalized'
 
-  call Roothaan_Hall_SCF
+  write(json_unit,json_array_open_fmt) 'scf'
-  call ezfio_set_hartree_fock_energy(SCF_energy)
 
+  call Roothaan_Hall_SCF
+
+  write(json_unit,json_array_close_fmtx)
+  call json_close
+
+  call ezfio_set_hartree_fock_energy(SCF_energy)
 end
 
 

src/iterations/EZFIO.cfg

@@ -1,24 +0,0 @@
-[n_iter]
-interface: ezfio
-doc: Number of saved iterations
-type:integer
-default: 1
-
-[n_det_iterations]
-interface: ezfio, provider
-doc: Number of determinants at each iteration
-type: integer
-size: (100)
- 
-[energy_iterations]
-interface: ezfio, provider
-doc: The variational energy at each iteration
-type: double precision 
-size: (determinants.n_states,100)
-
-[pt2_iterations]
-interface: ezfio, provider
-doc: The |PT2| correction at each iteration
-type: double precision 
-size: (determinants.n_states,100)
-

src/iterations/io.irp.f

@@ -1,37 +0,0 @@
-BEGIN_PROVIDER [ integer, n_iter  ]
-  implicit none
-  BEGIN_DOC
-! number of iterations
-  END_DOC
-
-  logical                        :: has
-  PROVIDE ezfio_filename
-  if (mpi_master) then
-
-      double precision :: zeros(N_states,100)
-      integer :: izeros(100)
-      zeros = 0.d0
-      izeros = 0
-      call ezfio_set_iterations_n_iter(0)
-      call ezfio_set_iterations_energy_iterations(zeros)
-      call ezfio_set_iterations_pt2_iterations(zeros)
-      call ezfio_set_iterations_n_det_iterations(izeros)
-      n_iter = 1
-  endif
-  IRP_IF MPI_DEBUG
-    print *,  irp_here, mpi_rank
-    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
-  IRP_ENDIF
-  IRP_IF MPI
-    include 'mpif.h'
-    integer :: ierr
-    call MPI_BCAST( n_iter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
-    if (ierr /= MPI_SUCCESS) then
-      stop 'Unable to read n_iter with MPI'
-    endif
-  IRP_ENDIF
-
-  call write_time(6)
-
-END_PROVIDER
-

src/iterations/iterations.irp.f

@@ -1,42 +1,65 @@
-BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter,N_states) ]
+BEGIN_PROVIDER [ integer, N_iter  ]
   implicit none
   BEGIN_DOC
- ! Extrapolated energy, using E_var = f(PT2) where PT2=0
+! Number of CIPSI iterations
   END_DOC
- integer :: i
+
- do i=1,min(N_states,N_det)
+  N_iter = 0
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, N_iter_max ]
+ implicit none
+ BEGIN_DOC
+ ! Max number of iterations for extrapolations
+ END_DOC
+ N_iter_max = 8
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, energy_iterations , (n_states,N_iter_max) ]
+&BEGIN_PROVIDER [ double precision, pt2_iterations , (n_states,N_iter_max) ]
+&BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter_max,N_states) ]
+  implicit none
+  BEGIN_DOC
+! The energy at each iteration for the extrapolations
+  END_DOC
+
+   energy_iterations = 0.d0
+   pt2_iterations = 0.d0
+   extrapolated_energy = 0.d0
+END_PROVIDER
+
+subroutine increment_n_iter(e, pt2_data)
+  use selection_types
+  implicit none
+  BEGIN_DOC
+! Does what is necessary to increment n_iter
+  END_DOC
+  double precision, intent(in) :: e(*)
+  type(pt2_type), intent(in)   :: pt2_data
+  integer :: k, i
+
+  if (N_det < N_states) return
+
+  if (N_iter < N_iter_max) then
+    N_iter += 1
+  else
+    do k=2,N_iter
+      energy_iterations(1:N_states,k-1) = energy_iterations(1:N_states,k)
+      pt2_iterations(1:N_states,k-1) = pt2_iterations(1:N_states,k)
+    enddo
+  endif
+  energy_iterations(1:N_states,N_iter) = e(1:N_states)
+  pt2_iterations(1:N_states,N_iter) = pt2_data % rpt2(1:N_states)
+
+  if (N_iter < 2) then
+    extrapolated_energy(1,:) = energy_iterations(:,1) + pt2_iterations(:,1)
+    extrapolated_energy(2,:) = energy_iterations(:,2) + pt2_iterations(:,2)
+  else
+    do i=1,N_states
       call extrapolate_data(N_iter,                               &
           energy_iterations(i,1:N_iter),                          &
              pt2_iterations(i,1:N_iter),                          &
           extrapolated_energy(1:N_iter,i))
     enddo
-END_PROVIDER
-
-
-subroutine save_iterations(e_, pt2_,n_)
-  implicit none
-  BEGIN_DOC
-! Update the energy in the EZFIO file.
-  END_DOC
-  integer, intent(in) :: n_
-  double precision, intent(in) :: e_(N_states), pt2_(N_states)
-  integer :: i
-
-  if (N_iter == 101) then
-    do i=2,N_iter-1
-      energy_iterations(1:N_states,N_iter-1) = energy_iterations(1:N_states,N_iter)
-      pt2_iterations(1:N_states,N_iter-1) = pt2_iterations(1:N_states,N_iter) 
-    enddo
-    N_iter = N_iter-1
-    TOUCH N_iter
   endif
-
-  energy_iterations(1:N_states,N_iter) = e_(1:N_states)
-     pt2_iterations(1:N_states,N_iter) = pt2_(1:N_states)
-  n_det_iterations(N_iter) = n_
-  call ezfio_set_iterations_N_iter(N_iter)
-  call ezfio_set_iterations_energy_iterations(energy_iterations)
-  call ezfio_set_iterations_pt2_iterations(pt2_iterations)
-  call ezfio_set_iterations_n_det_iterations(n_det_iterations)
 end
-

src/iterations/print_extrapolation.irp.f

@@ -5,10 +5,14 @@ subroutine print_extrapolated_energy
   END_DOC
 
   integer :: i,k
+  integer :: N_states_p, N_iter_p
 
   if (N_iter< 2) then
     return
   endif
+  N_states_p = min(N_states,N_det)
+  N_iter_p = min(N_iter, 8)
+
   write(*,'(A)') ''
   write(*,'(A)') 'Extrapolated energies'
   write(*,'(A)') '------------------------'
@@ -20,20 +24,20 @@ subroutine print_extrapolated_energy
   write(*,*)  '=========== ', '==================='
   write(*,*)  'minimum PT2 ', 'Extrapolated energy'
   write(*,*)  '=========== ', '==================='
-  do k=2,min(N_iter,8)
+  do k=2,N_iter_p
-    write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,N_iter+1-k), extrapolated_energy(k,1)
+    write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,N_iter_p+1-k), extrapolated_energy(k,1)
   enddo
   write(*,*)  '=========== ', '==================='
 
-  do i=2, min(N_states,N_det)
+  do i=2, N_states_p
     print *,  ''
     print *,  'State ', i
     print *,  ''
     write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
     write(*,*)  'minimum PT2 ', 'Extrapolated energy ', '  Excitation (a.u)  ', '  Excitation (eV)  '
     write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
-    do k=2,min(N_iter,8)
+    do k=2,N_iter_p
-      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), &
+      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,k), extrapolated_energy(k,i), &
           extrapolated_energy(k,i) - extrapolated_energy(k,1), &
           (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
     enddo

src/iterations_tc/EZFIO.cfg

@@ -1,24 +0,0 @@
-[n_iter]
-interface: ezfio
-doc: Number of saved iterations
-type:integer
-default: 1
-
-[n_det_iterations]
-interface: ezfio, provider
-doc: Number of determinants at each iteration
-type: integer
-size: (100)
- 
-[energy_iterations]
-interface: ezfio, provider
-doc: The variational energy at each iteration
-type: double precision 
-size: (determinants.n_states,100)
-
-[pt2_iterations]
-interface: ezfio, provider
-doc: The |PT2| correction at each iteration
-type: double precision 
-size: (determinants.n_states,100)
-

src/iterations_tc/io.irp.f

@@ -1,37 +0,0 @@
-BEGIN_PROVIDER [ integer, n_iter  ]
-  implicit none
-  BEGIN_DOC
-! number of iterations
-  END_DOC
-
-  logical                        :: has
-  PROVIDE ezfio_filename
-  if (mpi_master) then
-
-      double precision :: zeros(N_states,100)
-      integer :: izeros(100)
-      zeros = 0.d0
-      izeros = 0
-      call ezfio_set_iterations_n_iter(0)
-      call ezfio_set_iterations_energy_iterations(zeros)
-      call ezfio_set_iterations_pt2_iterations(zeros)
-      call ezfio_set_iterations_n_det_iterations(izeros)
-      n_iter = 1
-  endif
-  IRP_IF MPI_DEBUG
-    print *,  irp_here, mpi_rank
-    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
-  IRP_ENDIF
-  IRP_IF MPI
-    include 'mpif.h'
-    integer :: ierr
-    call MPI_BCAST( n_iter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
-    if (ierr /= MPI_SUCCESS) then
-      stop 'Unable to read n_iter with MPI'
-    endif
-  IRP_ENDIF
-
-  call write_time(6)
-
-END_PROVIDER
-

src/iterations_tc/iterations.irp.f

@@ -1,43 +0,0 @@
-BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter,N_states) ]
- implicit none
- BEGIN_DOC
- ! Extrapolated energy, using E_var = f(PT2) where PT2=0
- END_DOC
-! integer :: i
- extrapolated_energy = 0.D0
-END_PROVIDER 
-
- subroutine get_extrapolated_energy(Niter,ept2,pt1,extrap_energy)
- implicit none
- integer, intent(in)  :: Niter
- double precision, intent(in) :: ept2(Niter),pt1(Niter),extrap_energy(Niter)
- call extrapolate_data(Niter,ept2,pt1,extrap_energy)
- end
-
-subroutine save_iterations(e_, pt2_,n_)
-  implicit none
-  BEGIN_DOC
-! Update the energy in the EZFIO file.
-  END_DOC
-  integer, intent(in) :: n_
-  double precision, intent(in) :: e_(N_states), pt2_(N_states)
-  integer :: i
-
-  if (N_iter == 101) then
-    do i=2,N_iter-1
-      energy_iterations(1:N_states,N_iter-1) = energy_iterations(1:N_states,N_iter)
-      pt2_iterations(1:N_states,N_iter-1) = pt2_iterations(1:N_states,N_iter) 
-    enddo
-    N_iter = N_iter-1
-    TOUCH N_iter
-  endif
-
-  energy_iterations(1:N_states,N_iter) = e_(1:N_states)
-     pt2_iterations(1:N_states,N_iter) = pt2_(1:N_states)
-  n_det_iterations(N_iter) = n_
-  call ezfio_set_iterations_N_iter(N_iter)
-  call ezfio_set_iterations_energy_iterations(energy_iterations)
-  call ezfio_set_iterations_pt2_iterations(pt2_iterations)
-  call ezfio_set_iterations_n_det_iterations(n_det_iterations)
-end
-

src/iterations_tc/print_extrapolation.irp.f

@@ -1,46 +0,0 @@
-subroutine print_extrapolated_energy
-  implicit none
-  BEGIN_DOC
-! Print the extrapolated energy in the output
-  END_DOC
-
-  integer :: i,k
-
-  if (N_iter< 2) then
-    return
-  endif
-  write(*,'(A)') ''
-  write(*,'(A)') 'Extrapolated energies'
-  write(*,'(A)') '------------------------'
-  write(*,'(A)') ''
-
-  print *,  ''
-  print *,  'State ', 1
-  print *,  ''
-  write(*,*)  '=========== ', '==================='
-  write(*,*)  'minimum PT2 ', 'Extrapolated energy'
-  write(*,*)  '=========== ', '==================='
-  do k=2,min(N_iter,8)
-    write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,N_iter+1-k), extrapolated_energy(k,1)
-  enddo
-  write(*,*)  '=========== ', '==================='
-
-  do i=2, min(N_states,N_det)
-    print *,  ''
-    print *,  'State ', i
-    print *,  ''
-    write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
-    write(*,*)  'minimum PT2 ', 'Extrapolated energy ', '  Excitation (a.u)  ', '  Excitation (eV)  '
-    write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
-    do k=2,min(N_iter,8)
-      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), &
-          extrapolated_energy(k,i) - extrapolated_energy(k,1), &
-          (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
-    enddo
-    write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
-  enddo
-
-  print *,  ''
-
-end subroutine
-

src/iterations_tc/print_summary.irp.f

@@ -1,104 +0,0 @@
-subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s2_)
-  use selection_types
-  implicit none
-  BEGIN_DOC
-! Print the extrapolated energy in the output
-  END_DOC
-
-  integer, intent(in)            :: n_det_, n_configuration_, n_st
-  double precision, intent(in)   :: e_(n_st), s2_(n_st)
-  type(pt2_type)  , intent(in)   :: pt2_data, pt2_data_err
-  integer                        :: i, k
-  integer                        :: N_states_p
-  character*(9)                  :: pt2_string
-  character*(512)                :: fmt
-
-  if (do_pt2) then
-    pt2_string = '        '
-  else
-    pt2_string = '(approx)'
-  endif
-
-  N_states_p = min(N_det_,n_st)
-
-  print *, ''
-  print '(A,I12)',  'Summary at N_det = ', N_det_
-  print '(A)',      '-----------------------------------'
-  print *, ''
-
-  write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
-  write(*,fmt)
-  write(fmt,*) '(13X,', N_states_p, '(6X,A7,1X,I6,10X))'
-  write(*,fmt) ('State',k, k=1,N_states_p)
-  write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
-  write(*,fmt)
-  write(fmt,*) '(A13,', N_states_p, '(1X,F14.8,15X))'
-  write(*,fmt) '# E          ', e_(1:N_states_p)
-  if (N_states_p > 1) then
-    write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1)
-    write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0
-  endif
-  write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))'
-  write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p)
-  write(*,fmt) '# rPT2'//pt2_string, (pt2_data % rpt2(k), pt2_data_err % rpt2(k), k=1,N_states_p)
-  write(*,'(A)') '#'
-  write(*,fmt) '# E+PT2      ', (e_(k)+pt2_data % pt2(k),pt2_data_err % pt2(k), k=1,N_states_p)
-  write(*,fmt) '# E+rPT2     ', (e_(k)+pt2_data % rpt2(k),pt2_data_err % rpt2(k), k=1,N_states_p)
-  if (N_states_p > 1) then
-    write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), &
-      dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p)
-    write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, &
-      dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p)
-  endif
-  write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
-  write(*,fmt)
-  print *,  ''
-
-  print *,  'N_det             = ', N_det_
-  print *,  'N_states          = ', n_st
-  if (s2_eig) then
-    print *,  'N_cfg             = ', N_configuration_
-    if (only_expected_s2) then
-      print *,  'N_csf             = ', N_csf
-    endif
-  endif
-  print *,  ''
-
-  do k=1, N_states_p
-    print*,'* State ',k
-    print *,  '< S^2 >         = ', s2_(k)
-    print *,  'E               = ', e_(k)
-    print *,  'Variance        = ', pt2_data % variance(k), ' +/- ', pt2_data_err % variance(k)
-    print *,  'PT norm         = ', dsqrt(pt2_data % overlap(k,k)), ' +/- ', 0.5d0*dsqrt(pt2_data % overlap(k,k)) * pt2_data_err % overlap(k,k) / (pt2_data % overlap(k,k))
-    print *,  'PT2             = ', pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k)
-    print *,  'rPT2            = ', pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k)
-    print *,  'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k)
-    print *,  'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k)
-    print *,  ''
-  enddo
-
-  print *,  '-----'
-  if(n_st.gt.1)then
-    print *, 'Variational Energy difference (au | eV)'
-    do i=2, N_states_p
-      print*,'Delta E = ', (e_(i) - e_(1)), &
-        (e_(i) - e_(1)) * 27.211396641308d0
-    enddo
-    print *,  '-----'
-    print*, 'Variational + perturbative Energy difference (au | eV)'
-    do i=2, N_states_p
-      print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), &
-        (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0
-    enddo
-    print *,  '-----'
-    print*, 'Variational + renormalized perturbative Energy difference (au | eV)'
-    do i=2, N_states_p
-      print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), &
-        (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0
-    enddo
-  endif
-
-!  call print_energy_components()
-
-end subroutine
-

src/json/EZFIO.cfg

@@ -0,0 +1,5 @@
+[empty]
+type: logical
+doc: Needed to create the json directory
+interface: ezfio
+

src/json/NEED

@@ -0,0 +1 @@
+ezfio_files

src/json/README.rst

@@ -0,0 +1,5 @@
+====
+json
+====
+
+JSON files to simplify getting output information from QP.

src/json/json.irp.f

@@ -0,0 +1,45 @@
+BEGIN_PROVIDER [ character*(128), json_filename ]
+ implicit none
+ BEGIN_DOC
+ ! Fortran unit of the JSON file
+ END_DOC
+ integer, external            :: getUnitAndOpen
+ integer                      :: counter
+ character*(128)              :: prefix
+ logical                      :: exists
+
+ prefix = trim(ezfio_filename)//'/json/'
+
+ call lock_io
+ exists = .True.
+ counter = 0
+ do while (exists)
+   counter += 1
+   write(json_filename, '(A,I5.5,A)') trim(prefix), counter, '.json'
+   INQUIRE(FILE=trim(json_filename), EXIST=exists)
+ enddo
+ call unlock_io
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, json_unit]
+ implicit none
+ BEGIN_DOC
+ ! Unit file for JSON output
+ END_DOC
+ integer, external :: getUnitAndOpen
+ call ezfio_set_json_empty(.False.)
+ call lock_io
+ json_unit = getUnitAndOpen(json_filename, 'w')
+ write(json_unit, '(A)') '{'
+ call unlock_io
+END_PROVIDER
+
+subroutine json_close
+ call lock_io
+ write(json_unit, '(A)') '}'
+ close(json_unit)
+ call unlock_io
+ FREE json_unit
+end
+

src/json/json_formats.irp.f

@@ -0,0 +1,46 @@
+ BEGIN_PROVIDER [ character*(64), json_int_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_int_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_real_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_real_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_str_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_str_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_true_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_true_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_false_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_false_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_array_open_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_array_uopen_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_array_close_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_array_close_uopen_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_array_close_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_dict_open_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_dict_uopen_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_dict_close_uopen_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_dict_close_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_dict_close_fmtx ]
+ implicit none
+ BEGIN_DOC
+ ! Formats for JSON output.
+ ! x: used to mark the last write (no comma)
+ END_DOC
+ json_int_fmt    = '(''   "'',A,''": '',I10,'','')'
+ json_int_fmtx   = '(''   "'',A,''": '',I10)'
+ json_real_fmt   = '(''   "'',A,''": '',E22.15,'','')'
+ json_real_fmtx  = '(''   "'',A,''": '',E22.15)'
+ json_str_fmt    = '(''   "'',A,''": "'',A,''",'')'
+ json_str_fmtx   = '(''   "'',A,''": "'',A,''"'')'
+ json_true_fmt   = '(''   "'',A,''": true,'')'
+ json_true_fmtx  = '(''   "'',A,''": true'')'
+ json_false_fmt  = '(''   "'',A,''": false,'')'
+ json_false_fmtx = '(''   "'',A,''": false'')'
+ json_array_open_fmt   = '(''   "'',A,''": ['')'
+ json_array_uopen_fmt  = '(''   ['')'
+ json_array_close_fmt  = '(''   ],'')'
+ json_array_close_uopen_fmt  = '(''   ], ['')'
+ json_array_close_fmtx = '(''   ]'')'
+ json_dict_open_fmt    = '(''   "'',A,''": {'')'
+ json_dict_uopen_fmt   = '(''   {'')'
+ json_dict_close_fmt   = '(''   },'')'
+ json_dict_close_uopen_fmt   = '(''   }, {'')'
+ json_dict_close_fmtx  = '(''   }'')'
+END_PROVIDER

src/kohn_sham/ks_scf.irp.f

@@ -90,7 +90,11 @@ subroutine run
 
 ! Choose SCF algorithm
 
+  write(json_unit,*) '"scf" : ['
   call Roothaan_Hall_SCF
+  write(json_unit,*) ']'
+
+  call json_close
 
 end
 

src/kohn_sham/print_mos.irp.f

@@ -0,0 +1,23 @@
+program print_mos
+ implicit none
+ integer :: i,nx
+ double precision :: r(3), xmax, dx, accu
+ double precision, allocatable :: mos_array(:)
+ double precision:: alpha,envelop
+ allocate(mos_array(mo_num))
+ xmax = 5.d0
+ nx = 1000
+ dx=xmax/dble(nx)
+ r = 0.d0
+ alpha = 0.5d0
+ do i = 1, nx
+  call give_all_mos_at_r(r,mos_array)
+  accu = mos_array(3)**2+mos_array(4)**2+mos_array(5)**2
+  accu = dsqrt(accu)
+  envelop = (1.d0 - dexp(-alpha * r(3)**2))
+  write(33,'(100(F16.10,X))')r(3), mos_array(1), mos_array(2), accu, envelop
+  r(3) += dx
+ enddo
+
+end
+

src/kohn_sham_rs/rs_ks_scf.irp.f

@@ -93,7 +93,10 @@ subroutine run
 
   level_shift += 1.d0
   touch level_shift
+  write(json_unit,*) '"scf" : ['
   call Roothaan_Hall_SCF
+  write(json_unit,*) ']'
+  call json_close
   call ezfio_set_kohn_sham_rs_energy(SCF_energy)
 
  write(*, '(A22,X,F16.10)') 'one_e_energy = ',one_e_energy

src/mo_basis/mo_class.irp.f

@@ -1,8 +1,3 @@
-! DO NOT MODIFY BY HAND
-! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
-! from file /home/eginer/programs/qp2/src/mo_basis/EZFIO.cfg
-
-
 BEGIN_PROVIDER [ character*(32), mo_class , (mo_num) ]
   implicit none
   BEGIN_DOC
@@ -35,6 +30,4 @@ BEGIN_PROVIDER [ character*(32), mo_class , (mo_num) ]
     endif
   IRP_ENDIF
 
-  call write_time(6)
-
 END_PROVIDER