Merge branch 'dev' of github.com:QuantumPackage/qp2 into dev

bin/qp_gaussian

@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+#
+"""
+Runs a Quantum Package calculation using a Gaussian input file.
+
+Usage:
+      qp_gaussian INPUT
+
+"""
+
+# Requires pymatgen (https://pymatgen.org/)
+# pip install pymatgen
+
+
+import os
+import sys
+import os.path
+
+try:
+    import qp_path
+except ImportError:
+    print("source quantum_package.rc")
+
+from docopt import docopt
+import pymatgen
+from pymatgen.io.gaussian import GaussianInput
+
+
+def main(arguments):
+
+    filename = arguments["INPUT"]
+
+    with open(filename,'r') as f:
+        text = f.read()
+
+    in_file = GaussianInput.from_string(text)
+
+    d = in_file.as_dict()
+    charge = ("%d"%(d["charge"])).replace('-','m')
+    basis  = d["basis_set"]
+    mult   = d["spin_multiplicity"]
+    natoms = len(d["molecule"]["sites"])
+    with open("g09.xyz","w") as f:
+      f.write("%d\n"%natoms)
+      f.write("%s\n"%d["title"])
+      f.write("%s\n"%in_file.get_cart_coords())
+
+    if basis is None:
+       print("Basis set not found. Use '/' before basis set")
+       sys.exit(1)
+
+    command = f"rm -rf g09.ezfio"
+    os.system(command)
+
+    command = f"qp_create_ezfio -c {charge} -m {mult} g09.xyz -b {basis} -o g09.ezfio"
+    os.system(command)
+
+    command = f"rm -rf g09.xyz"
+    os.system(command)
+
+    command = f"qp_run scf g09.ezfio"
+    os.system(command)
+
+    command = f"qp_set_frozen_core g09.ezfio"
+    os.system(command)
+
+    if d["functional"] == "FCI":
+      command = f"qp_run fci g09.ezfio"
+    elif d["functional"] == "CIS":
+      command = f"qp_run cis g09.ezfio"
+    elif d["functional"] == "CISD":
+      command = f"qp_run cisd g09.ezfio"
+
+    os.system(command)
+
+    
+
+
+
+
+if __name__ == '__main__':
+    ARGUMENTS = docopt(__doc__)
+    main(ARGUMENTS)

config/gfortran.cfg

@@ -13,7 +13,7 @@
 FC           : gfortran -g -ffree-line-length-none -I . -fPIC
 LAPACK_LIB   : -lblas -llapack
 IRPF90       : irpf90
-IRPF90_FLAGS : --ninja --align=32 --assert 
+IRPF90_FLAGS : --ninja --align=32 --assert -DGNU_CHECK_OMP
 
 # Global options
 ################

configure

@@ -3,8 +3,6 @@
 # Quantum Package configuration script
 #
 
-TEMP=$(getopt -o d:c:i:h -l download:,config:,install:,help -n $0 -- "$@") || exit 1
-eval set -- "$TEMP"
 
 export QP_ROOT="$( cd "$(dirname "$0")" ; pwd -P )"
 echo "QP_ROOT="$QP_ROOT
@@ -24,15 +22,15 @@ function help()
 Quantum Package configuration script.
 
 Usage:
-  $(basename $0) -c <file>    | --config=<file>
-  $(basename $0) -h           | --help
-  $(basename $0) -i <package> | --install=<package>
+  $(basename $0) -c <file>    
+  $(basename $0) -h           
+  $(basename $0) -i <package> 
 
 Options:
-  -c, --config=<file>      Define a COMPILATION configuration file,
+  -c  <file>    Define a COMPILATION configuration file,
                 in "${QP_ROOT}/config/".
-  -h, --help               Print the HELP message
-  -i, --install=<package>  INSTALL <package>. Use at your OWN RISK:
+  -h            Print the HELP message
+  -i <package>  INSTALL <package>. Use at your OWN RISK:
                 no support will be provided for the installation of
                 dependencies.
 
@@ -68,32 +66,31 @@ function execute () {
 }
 
 PACKAGES=""
+echo $@
 
-while true ; do
-    case "$1" in
-        -c|--config)
-            case "$2" in
+
+while getopts "d:c:i:h" c ; do
+    case "$c" in
+        c)
+            case "$OPTARG" in
                 "") help ; break;;
-                *) if [[ -f $2 ]] ; then
-                      CONFIG="$2"
+                *) if [[ -f $OPTARG ]] ; then
+                      CONFIG="$OPTARG"
                    else
-                      error "error: configuration file $2 not found."
+                      error "error: configuration file $OPTARG not found."
                       exit 1
                    fi
-            esac
-            shift 2;;
-        -i|--install)
-            case "$2" in
+            esac;;
+        i)
+            case "$OPTARG" in
                 "") help ; break;;
-                *) PACKAGES="${PACKAGE} $2"
-            esac
-            shift 2;;
-        -h|-help|--help)
+                *) PACKAGES="${PACKAGE} $OPTARG"
+            esac;;
+        h)
             help
             exit 0;;
-        --) shift ; break ;;
         *)
-            error $(basename $0)": unknown option $1, try --help"
+            error $(basename $0)": unknown option $c, try --help"
             exit 2;;
     esac
 done

scripts/verif_omp/check_actual_setup.sh

@@ -0,0 +1,12 @@
+#!/bin/sh
+
+# go in qp2/src/fci to run check_omp_actual_setup
+# to see if we can run in parallel an omp section in another one
+echo ""
+echo "Please wait..."
+echo ""
+cd ../../src/fci
+ninja || echo "Please recompile from the root"  
+echo ""
+./check_omp_actual_setup
+cd ../../scripts/verif_omp

scripts/verif_omp/check_omp_v2.f90

@@ -0,0 +1,175 @@
+program check_omp_v2
+ 
+  use omp_lib
+
+  implicit none
+
+  integer :: accu, accu2
+  integer :: s, n_setting
+  logical :: verbose, test_versions
+  logical, allocatable :: is_working(:)
+
+  verbose = .False.
+  test_versions = .True.
+  n_setting = 4
+
+  allocate(is_working(n_setting))
+
+  is_working = .False.
+
+  ! set the number of threads
+  call omp_set_num_threads(2) 
+
+  do s = 1, n_setting
+
+    accu = 0
+    accu2 = 0
+
+    call omp_set_max_active_levels(1)
+    call omp_set_nested(.False.)
+
+    if (s==1) then
+      !call set_multiple_levels_omp()
+      cycle
+    elseif (s==2) then
+      call omp_set_max_active_levels(5)
+    elseif (s==3) then
+      call omp_set_nested(.True.)
+    else
+      call omp_set_nested(.True.)
+      call omp_set_max_active_levels(5)
+    endif
+
+    ! Level 1
+    !$OMP PARALLEL 
+      if (verbose) then
+        print*,'Num threads level 1:',omp_get_num_threads()
+      endif
+
+      ! Level 2
+      !$OMP PARALLEL 
+        if (verbose) then
+          print*,'Num threads level 2:',omp_get_num_threads()
+        endif
+
+        ! Level 3
+        !$OMP PARALLEL 
+          if (verbose) then
+            print*,'Num threads level 3:',omp_get_num_threads()
+          endif
+
+          call check_omp_in_subroutine(accu2)
+     
+          ! Level 4
+          !$OMP PARALLEL 
+
+            if (verbose) then
+              print*,'Num threads level 4:',omp_get_num_threads()
+            endif
+
+           !$OMP ATOMIC    
+           accu = accu + 1
+           !$OMP END ATOMIC 
+       
+          !$OMP END PARALLEL
+     
+     
+        !$OMP END PARALLEL
+
+
+      !$OMP END PARALLEL
+
+
+    !$OMP END PARALLEL
+
+  if (verbose) then
+    print*,'Setting:',s,'accu=',accu
+    print*,'Setting:',s,'accu2=',accu2
+  endif
+
+  if (accu == 16 .and. accu2 == 16) then
+    is_working(s) = .True.
+  endif
+
+  enddo
+
+  if (verbose) then
+    if (is_working(2)) then
+      print*,'The parallelization works on 4 levels with:'
+      print*,'call omp_set_max_active_levels(5)'
+      print*,''
+      print*,'Please use the irpf90 flags -DSET_MAX_ACT in qp2/config/${compiler_name}.cfg'
+    elseif (is_working(3)) then
+      print*,'The parallelization works on 4 levels with:'
+      print*,'call omp_set_nested(.True.)'
+      print*,''
+      print*,'Please use the irpf90 flag -DSET_NESTED in qp2/config/${compiler_name}.cfg'
+    elseif (is_working(4)) then
+      print*,'The parallelization works on 4 levels with:'
+      print*,'call omp_set_nested(.True.)'
+      print*,'+'
+      print*,'call omp_set_max_active_levels(5)'
+      print*,''
+      print*,'Please use the irpf90 flags -DSET_NESTED -DSET_MAX_ACT in qp2/config/${compiler_name}.cfg'
+    else
+      print*,'The parallelization on multiple levels does not work with:'
+      print*,'call omp_set_max_active_levels(5)'
+      print*,'or'
+      print*,'call omp_set_nested(.True.)'
+      print*,'or'
+      print*,'call omp_set_nested(.True.)'
+      print*,'+'
+      print*,'call omp_set_max_active_levels(5)'
+      print*,''
+      print*,'Try an other compiler and good luck...'
+    endif
+
+   ! if (is_working(1)) then
+   !   print*,''
+   !   print*,'=========================================================='
+   !   print*,'Your actual set up works for parallelization with 4 levels'
+   !   print*,'=========================================================='
+   !   print*,''
+   ! else
+   !   print*,''
+   !   print*,'==================================================================='
+   !   print*,'Your actual set up does not work for parallelization with 4 levels'
+   !   print*,'Please look at the previous messages to understand the requirements'
+   !   print*,'==================================================================='
+   !   print*,''
+   ! endif
+  endif
+
+  ! List of working flags
+  if (test_versions) then
+    print*,'Tests:',is_working(2:4)
+  endif
+
+  ! IRPF90_FLAGS
+  if (is_working(2)) then
+    print*,'-DSET_MAX_ACT'
+  elseif (is_working(3)) then
+    print*,'-DSET_NESTED'
+  elseif (is_working(4)) then
+    print*,'-DSET_MAX_ACT -DSET_NESTED'
+  else
+    print*,'ERROR'
+  endif
+
+end
+
+subroutine check_omp_in_subroutine(accu2)
+
+  implicit none
+ 
+  integer, intent(inout) :: accu2
+
+  !$OMP PARALLEL
+
+  !$OMP ATOMIC
+  accu2 = accu2 + 1
+  !$OMP END ATOMIC
+
+  !$OMP END PARALLEL
+
+end

scripts/verif_omp/check_required_setup.sh

@@ -0,0 +1,19 @@
+#!/bin/sh
+
+# take one argument which is the compiler used
+# return the required IRPF90_FLAGS for the $1 compiler
+
+if [ -z "$1" ]
+then
+    echo "Give the compiler in argument"
+else
+
+$1 --version > /dev/null \
+&& $1 -O0 -fopenmp check_omp_v2.f90 \
+&& ./a.out | tail -n 1
+
+
+# if there is an error or if the compiler is not found
+$1 --version > /dev/null || echo 'compiler not found'
+
+fi

scripts/verif_omp/study_omp.sh

@@ -0,0 +1,30 @@
+#!/bin/sh
+
+# list of compilers
+list_comp="ifort gfortran-7 gfortran-8 gfortran-9"
+
+# file to store the results
+FILE=results.dat
+
+touch $FILE
+rm $FILE
+
+# Comments
+echo "1: omp_set_max_active_levels(5)" >> $FILE
+echo "2: omp_set_nested(.True.)" >> $FILE
+echo "3: 1 + 2" >> $FILE
+echo "" >> $FILE
+echo "1 2 3" >> $FILE
+
+# loop on the comp
+for comp in $list_comp
+do
+	$comp --version > /dev/null \
+        && $comp -O0 -fopenmp check_omp_v2.f90 \
+	&& echo $(./a.out | grep "Tests:" | cut -d ":" -f2- ) $(echo " : ") $($comp --version | head -n 1) >> $FILE
+
+done
+
+# Display
+cat $FILE
+

scripts/verif_omp/update_comp.sh

@@ -0,0 +1,49 @@
+#!/bin/bash
+
+# Compiler
+COMP=$1
+
+# Path to file.cfg
+config_PATH="../../config/"
+END="*.cfg"
+CONFIG="/config/"
+
+#LIST=${config_PATH}${COMP}${END} # without ${QP_ROOT}
+LIST=${QP_ROOT}${CONFIG}${COMP}${END}
+
+if [ -z "$1" ]
+then
+    echo "Give the compiler in argument"
+else
+
+    # List of the config files for the compiler
+    #list_files=$(ls ../../config/$comp*.cfg) #does not give the right list
+    list_files=${LIST}
+    echo "Files that will be modified:"
+    echo $list_files
+   
+    # Flags that must be added 
+    FLAGS=$(./check_required_setup.sh $COMP)
+
+    # Add the flags
+    for file in $list_files
+    do
+        echo $file
+        BASE="IRPF90_FLAGS : --ninja"
+        ACTUAL=$(grep "$BASE" $file)
+        
+        # To have only one time each flag	
+	grep " -DSET_MAX_ACT" $file && ${ACTUAL/" -DSET_MAX"/""}
+	grep " -DSET_NESTED" $file && ${ACTUAL/" -DSET_NESTED"/""}
+        SPACE=" "
+
+	NEW=${ACTUAL}${SPACE}${FLAGS}
+
+	# Debug
+        #echo ${NEW}
+
+        sed "s/${ACTUAL}/${NEW}/" $file
+        # -i # to change the files
+    done
+
+fi

src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f

@@ -28,6 +28,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)]
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
+  use omp_lib
   implicit none
   BEGIN_DOC
   ! Local pseudo-potential
@@ -42,7 +43,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
 
   double precision               :: wall_1, wall_2, wall_0
   integer                        :: thread_num
-  integer                        :: omp_get_thread_num
   double precision               :: c
   double precision               :: Z
 
@@ -158,6 +158,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
 
 
  BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_non_local, (ao_num,ao_num)]
+  use omp_lib
   implicit none
   BEGIN_DOC
   ! Non-local pseudo-potential
@@ -169,7 +170,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
   integer                        :: power_A(3),power_B(3)
   integer                        :: i,j,k,l,m
   double precision               :: Vloc, Vpseudo
-  integer                        :: omp_get_thread_num
 
   double precision               :: wall_1, wall_2, wall_0
   integer                        :: thread_num

src/cis/EZFIO.cfg

@@ -5,4 +5,3 @@ interface: ezfio
 size: (determinants.n_states)
 
 
-

src/cis/cis.irp.f

@@ -79,6 +79,6 @@ subroutine run
   call ezfio_set_cis_energy(CI_energy)
   psi_coef = ci_eigenvectors
   SOFT_TOUCH psi_coef
-  call save_wavefunction_truncated(1.d-12)
+  call save_wavefunction_truncated(save_threshold)
 
 end

src/cisd/cisd.irp.f

@@ -63,7 +63,7 @@ subroutine run
   endif
   psi_coef = ci_eigenvectors
   SOFT_TOUCH psi_coef
-  call save_wavefunction
+  call save_wavefunction_truncated(save_threshold)
   call ezfio_set_cisd_energy(CI_energy)
 
   do i = 1,N_states

src/csf/sigma_vector.irp.f

@@ -59,6 +59,11 @@
       endif
     endif
     ncfg = ncfgpersomo - ncfgprev
+    if(i .EQ. 0 .OR. i .EQ. 1) then
+      dimcsfpercfg = 1
+    elseif( i .EQ. 3) then
+      dimcsfpercfg = 2
+    else
       if(iand(MS,1) .EQ. 0) then
         !dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1))))
         binom1 = dexp(logabsgamma(1.0d0*(i+1))                            &
@@ -78,6 +83,7 @@
                     - logabsgamma(1.0d0*(i-(((i+3)/2)+1)+1)));
         dimcsfpercfg = max(1,nint(binom1 - binom2))
       endif
+    endif
     n_CSF += ncfg * dimcsfpercfg
     if(cfg_seniority_index(i+2) > ncfgprev) then
       ncfgprev = cfg_seniority_index(i+2)

src/davidson/davidson_parallel.irp.f

@@ -508,7 +508,8 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
   endif
 
 
-  call omp_set_max_active_levels(5)
+  !call omp_set_max_active_levels(5)
+  call set_multiple_levels_omp()
 
   !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread)
   ithread = omp_get_thread_num()

src/davidson/davidson_parallel_csf.irp.f

@@ -464,7 +464,9 @@ subroutine H_u_0_nstates_zmq(v_0,u_0,N_st,sze)
     print *,  irp_here, ': Failed in zmq_set_running'
   endif
 
-  call omp_set_max_active_levels(4)
+  !call omp_set_max_active_levels(4)
+  call set_multiple_levels_omp()
+
   !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread)
   ithread = omp_get_thread_num()
   if (ithread == 0 ) then

src/davidson/davidson_parallel_nos2.irp.f

@@ -464,7 +464,9 @@ subroutine H_u_0_nstates_zmq(v_0,u_0,N_st,sze)
     print *,  irp_here, ': Failed in zmq_set_running'
   endif
 
-  call omp_set_max_active_levels(4)
+  !call omp_set_max_active_levels(4)
+  call set_multiple_levels_omp()
+
   !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread)
   ithread = omp_get_thread_num()
   if (ithread == 0 ) then

src/determinants/EZFIO.cfg

@@ -136,3 +136,9 @@ doc: If |true|, discard any Slater determinants with an interaction smaller than
 interface: ezfio,provider,ocaml
 default: False
 
+[save_threshold]
+type: Threshold
+doc: Cut-off to apply to the CI coefficients when the wave function is stored
+interface: ezfio,provider,ocaml
+default: 0.
+

src/fci/check_omp_actual_setup.irp.f

@@ -0,0 +1,174 @@
+program check_omp_actual_setup
+ 
+  use omp_lib
+
+  implicit none
+
+  integer :: accu, accu2
+  integer :: s, n_setting
+  logical :: verbose, test_versions
+  logical, allocatable :: is_working(:)
+
+  verbose = .True.
+  test_versions = .False.
+  n_setting = 4
+
+  allocate(is_working(n_setting))
+
+  is_working = .False.
+
+  ! set the number of threads
+  call omp_set_num_threads(2) 
+
+  do s = 1, n_setting
+
+    accu = 0
+    accu2 = 0
+
+    call omp_set_max_active_levels(1)
+    call omp_set_nested(.False.)
+
+    if (s==1) then
+      call set_multiple_levels_omp()
+    elseif (s==2) then
+      call omp_set_max_active_levels(5)
+    elseif (s==3) then
+      call omp_set_nested(.True.)
+    else
+      call omp_set_nested(.True.)
+      call omp_set_max_active_levels(5)
+    endif
+
+    ! Level 1
+    !$OMP PARALLEL 
+      if (verbose) then
+        print*,'Num threads level 1:',omp_get_num_threads()
+      endif
+
+      ! Level 2
+      !$OMP PARALLEL 
+        if (verbose) then
+          print*,'Num threads level 2:',omp_get_num_threads()
+        endif
+
+        ! Level 3
+        !$OMP PARALLEL 
+          if (verbose) then
+            print*,'Num threads level 3:',omp_get_num_threads()
+          endif
+
+          call check_omp_in_subroutine(accu2)
+     
+          ! Level 4
+          !$OMP PARALLEL 
+
+            if (verbose) then
+              print*,'Num threads level 4:',omp_get_num_threads()
+            endif
+
+           !$OMP ATOMIC    
+           accu = accu + 1
+           !$OMP END ATOMIC 
+       
+          !$OMP END PARALLEL
+     
+     
+        !$OMP END PARALLEL
+
+
+      !$OMP END PARALLEL
+
+
+    !$OMP END PARALLEL
+
+  if (verbose) then
+    print*,'Setting:',s,'accu=',accu
+    print*,'Setting:',s,'accu2=',accu2
+  endif
+
+  if (accu == 16 .and. accu2 == 16) then
+    is_working(s) = .True.
+  endif
+
+  enddo
+
+  if (verbose) then
+    if (is_working(2)) then
+      print*,'The parallelization works on 4 levels with:'
+      print*,'call omp_set_max_active_levels(5)'
+      print*,''
+      print*,'Please use the irpf90 flags -DSET_MAX_ACT in qp2/config/${compiler_name}.cfg'
+    elseif (is_working(3)) then
+      print*,'The parallelization works on 4 levels with:'
+      print*,'call omp_set_nested(.True.)'
+      print*,''
+      print*,'Please use the irpf90 flag -DSET_NESTED in qp2/config/${compiler_name}.cfg'
+    elseif (is_working(4)) then
+      print*,'The parallelization works on 4 levels with:'
+      print*,'call omp_set_nested(.True.)'
+      print*,'+'
+      print*,'call omp_set_max_active_levels(5)'
+      print*,''
+      print*,'Please use the irpf90 flags -DSET_NESTED -DSET_MAX_ACT in qp2/config/${compiler_name}.cfg'
+    else
+      print*,'The parallelization on multiple levels does not work with:'
+      print*,'call omp_set_max_active_levels(5)'
+      print*,'or'
+      print*,'call omp_set_nested(.True.)'
+      print*,'or'
+      print*,'call omp_set_nested(.True.)'
+      print*,'+'
+      print*,'call omp_set_max_active_levels(5)'
+      print*,''
+      print*,'Try an other compiler and good luck...'
+    endif
+
+    if (is_working(1)) then
+      print*,''
+      print*,'=========================================================='
+      print*,'Your actual set up works for parallelization with 4 levels'
+      print*,'=========================================================='
+      print*,''
+    else
+      print*,''
+      print*,'==================================================================='
+      print*,'Your actual set up does not work for parallelization with 4 levels'
+      print*,'Please look at the previous messages to understand the requirements'
+      print*,'==================================================================='
+      print*,''
+    endif
+  endif
+
+  ! List of working flags
+  if (test_versions) then
+    print*,is_working(2:4)
+  endif
+
+  ! IRPF90_FLAGS
+  if (is_working(2)) then
+    print*,'-DSET_MAX_ACT'
+  elseif (is_working(3)) then
+    print*,'-DSET_NESTED'
+  elseif (is_working(4)) then
+    print*,'-DSET_MAX_ACT -DSET_NESTED'
+  else
+    print*,'ERROR'
+  endif
+
+end
+
+subroutine check_omp_in_subroutine(accu2)
+
+  implicit none
+ 
+  integer, intent(inout) :: accu2
+
+  !$OMP PARALLEL
+
+  !$OMP ATOMIC
+  accu2 = accu2 + 1
+  !$OMP END ATOMIC
+
+  !$OMP END PARALLEL
+
+end

src/utils/set_multiple_levels_omp.irp.f

@@ -0,0 +1,16 @@
+subroutine set_multiple_levels_omp()
+
+! Doc : idk
+
+  implicit none
+
+  IRP_IF SET_MAX_ACT
+    !print*,'SET_MAX_ACT: True, call omp_set_max_active_levels(5)'
+    call omp_set_max_active_levels(5)
+  IRP_ENDIF
+  IRP_IF SET_NESTED
+    !print*,'SET_NESTED: True, call omp_set_nested(.True.)'
+    call omp_set_nested(.True.)
+  IRP_ENDIF
+
+end

src/utils/util.irp.f

@@ -300,12 +300,12 @@ subroutine wall_time(t)
 end
 
 BEGIN_PROVIDER [ integer, nproc ]
+  use omp_lib
   implicit none
   BEGIN_DOC
   ! Number of current OpenMP threads
   END_DOC
 
-  integer                        :: omp_get_num_threads
   nproc = 1
   !$OMP PARALLEL
   !$OMP MASTER