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

2023-10-26 15:06:57 +02:00 · 2023-10-26 15:06:57 +02:00 · bae76e9243
parent ec282e5b5d 3f38c0a9da
commit bae76e9243
21 changed files with 214 additions and 28 deletions
--- a/bin/python
+++ b/bin/python
@ -0,0 +1,4 @@
+#!/bin/bash
+
+exec python3 $@
+
--- a/bin/qpsh
+++ b/bin/qpsh
@ -1,6 +1,7 @@
 #!/bin/bash

-export QP_ROOT=$(dirname "$(readlink -f "$0")")/..
+REALPATH=$( cd "$(dirname "$0")" ; pwd -P )
+export QP_ROOT=${REALPATH}/..

 bash --init-file <(cat << EOF
    [[ -f /etc/bashrc ]] && source /etc/bashrc 
--- a/bin/zcat
+++ b/bin/zcat
@ -0,0 +1,22 @@
+#!/bin/bash
+
+# On Darwin: try gzcat if available, otherwise use Python
+
+if [[ $(uname -s) = Darwin ]] ; then
+   which gzcat &> /dev/null
+   if [[ $? -eq 0 ]] ; then 
+      exec gzcat $@
+   else
+
+      exec python3 << EOF
+import sys
+import gzip
+with gzip.open("$1", "rt") as f:
+    print(f.read())
+EOF
+   fi
+else
+   command=$(which -a zcat | grep -v 'qp2/bin/' | head -1)
+   exec $command $@
+fi
+
--- a/config/gfortran.cfg
+++ b/config/gfortran.cfg
@ -10,7 +10,7 @@
 #
 #
 [COMMON]
-FC           : gfortran -g -ffree-line-length-none -I . -fPIC -march=native
+FC           : gfortran -g -ffree-line-length-none -I . -fPIC -march=native -std=legacy
 LAPACK_LIB   : -lblas -llapack
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED
--- a/config/gfortran_armpl.cfg
+++ b/config/gfortran_armpl.cfg
@ -13,7 +13,7 @@
 #
 #
 [COMMON]
-FC           : gfortran -g -ffree-line-length-none -I . -fPIC -march=native
+FC           : gfortran -g -ffree-line-length-none -I . -fPIC -march=native -std=legacy
 LAPACK_LIB   : -larmpl_lp64_mp
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED
--- a/config/gfortran_avx.cfg
+++ b/config/gfortran_avx.cfg
@ -10,7 +10,7 @@
 #
 #
 [COMMON]
-FC           : gfortran -ffree-line-length-none -I . -mavx -g -fPIC
+FC           : gfortran -ffree-line-length-none -I . -mavx -g -fPIC -std=legacy
 LAPACK_LIB   : -llapack -lblas
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DSET_NESTED
--- a/config/gfortran_debug.cfg
+++ b/config/gfortran_debug.cfg
@ -10,7 +10,7 @@
 #
 #
 [COMMON]
-FC           : gfortran -g -ffree-line-length-none -I . -fPIC
+FC           : gfortran -g -ffree-line-length-none -I . -fPIC -std=legacy
 LAPACK_LIB   : -lblas -llapack
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED
--- a/config/gfortran_macos.cfg
+++ b/config/gfortran_macos.cfg
@ -0,0 +1,62 @@
+# Common flags
+##############
+#
+# -ffree-line-length-none : Needed for IRPF90 which produces long lines
+# -lblas -llapack         : Link with libblas and liblapack libraries provided by the system
+# -I .                    : Include the curent directory (Mandatory)
+#
+# --ninja                 : Allow the utilisation of ninja. (Mandatory)
+# --align=32              : Align all provided arrays on a 32-byte boundary
+#
+#
+[COMMON]
+FC           : gfortran -ffree-line-length-none -I . -g -fPIC -std=legacy
+LAPACK_LIB   : -llapack -lblas
+IRPF90       : irpf90
+IRPF90_FLAGS : --ninja --align=32 -DSET_NESTED -DMACOS
+
+# Global options
+################
+#
+# 1 : Activate
+# 0 : Deactivate
+#
+[OPTION]
+MODE    : OPT        ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
+CACHE   : 0          ; Enable cache_compile.py
+OPENMP  : 1          ; Append OpenMP flags
+
+# Optimization flags
+####################
+#
+# -Ofast                  : Disregard strict standards compliance. Enables all -O3 optimizations.
+#                           It also enables optimizations that are not valid
+#                           for all standard-compliant programs.  It turns on
+#                           -ffast-math and the Fortran-specific
+#                           -fno-protect-parens and -fstack-arrays.
+[OPT]
+FCFLAGS : -Ofast -march=native
+
+# Profiling flags
+#################
+#
+[PROFILE]
+FC      : -p -g
+FCFLAGS : -Ofast
+
+# Debugging flags
+#################
+#
+# -fcheck=all  : Checks uninitialized variables,  array subscripts, etc...
+# -g           : Extra debugging information
+#
+[DEBUG]
+FCFLAGS : -fcheck=all -g
+
+# OpenMP flags
+#################
+#
+[OPENMP]
+FC           : -fopenmp
+IRPF90_FLAGS : --openmp
+
--- a/config/gfortran_mpi.cfg
+++ b/config/gfortran_mpi.cfg
@ -10,7 +10,7 @@
 #
 #
 [COMMON]
-FC           : mpif90 -ffree-line-length-none -I . -g -fPIC
+FC           : mpif90 -ffree-line-length-none -I . -g -fPIC -std=legacy
 LAPACK_LIB   :  -lblas -llapack
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 -DMPI -DSET_NESTED
--- a/config/gfortran_openblas.cfg
+++ b/config/gfortran_openblas.cfg
@ -10,7 +10,7 @@
 #
 #
 [COMMON]
-FC           : gfortran -g -ffree-line-length-none -I . -fPIC -march=native
+FC           : gfortran -g -ffree-line-length-none -I . -fPIC -march=native -std=legacy
 LAPACK_LIB   : -lopenblas
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED
--- a/7
+++ b/7
@ -19,7 +19,11 @@ git submodule init
 git submodule update

 # Update ARM or x86 dependencies
-ARCHITECTURE=$(uname -m)
+SYSTEM=$(uname -s)
+if [[ $SYSTEM = "Linux" ]] ; then
+  SYSTEM=""
+fi
+ARCHITECTURE=$(uname -m)$SYSTEM
 cd ${QP_ROOT}/external/qp2-dependencies
 git checkout master
 git pull
@ -275,6 +279,7 @@ EOF
              cd "\${QP_ROOT}"/external
              tar --gunzip --extract --file qp2-dependencies/zeromq-4.2.5.tar.gz
              cd zeromq-*
+              [[ "${SYSTEM}" = "Darwin" ]] && ./autogen.sh
              ./configure --prefix="\$QP_ROOT" --without-libsodium --enable-libunwind=no
              make -j 8
              make install
--- a/scripts/module/create_executables_list.sh
+++ b/scripts/module/create_executables_list.sh
@ -11,7 +11,11 @@ fi

 cd ${QP_ROOT}/data
 rm -f executables
+if [[ "$(uname -s)" = "Darwin" ]] ; then
+EXES=$(find -L ${QP_ROOT}/src -maxdepth 2 -depth -perm +111 -type f | grep -e "${QP_ROOT}/src/[^/]*/[^/]*$" |sort ) 
+else
 EXES=$(find -L ${QP_ROOT}/src -maxdepth 2 -depth -executable -type f | grep -e "${QP_ROOT}/src/[^/]*/[^/]*$" |sort ) 
+fi

 for EXE in $EXES
 do
--- a/src/ezfio_files/00.create.bats
+++ b/src/ezfio_files/00.create.bats
@ -53,7 +53,6 @@ function run {


@test "B-B" {
-  qp set_file b2_stretched.ezfio
  run b2_stretched.zmt 1 0 6-31g
 }

--- a/src/ezfio_files/ezfio.irp.f
+++ b/src/ezfio_files/ezfio.irp.f
@ -33,6 +33,8 @@ BEGIN_PROVIDER [ character*(1024), ezfio_filename ]

  call ezfio_set_file(ezfio_filename)

+IRP_IF MACOS
+IRP_ELSE
  ! Adjust out-of-memory killer flag such that the current process will be
  ! killed first by the OOM killer, allowing compute nodes to survive
  integer :: getpid
@ -40,6 +42,7 @@ BEGIN_PROVIDER [ character*(1024), ezfio_filename ]
  write(pidc,*) getpid()
  write(command,*) 'echo 15 > /proc//'//trim(adjustl(pidc))//'/oom_adj'
  call system(command)
+IRP_ENDIF

  PROVIDE file_lock

--- a/src/hartree_fock/10.hf.bats
+++ b/src/hartree_fock/10.hf.bats
@ -39,7 +39,7 @@ rm -rf $EZFIO
 qp create_ezfio -b def2-svp hcn.xyz -o $EZFIO
 qp run scf
 mv hcn_charges.xyz ${EZFIO}_point_charges.xyz
-python write_pt_charges.py ${EZFIO}
+python3 write_pt_charges.py ${EZFIO}
 qp set nuclei point_charges True
 qp run scf | tee ${EZFIO}.pt_charges.out
  energy="$(ezfio get hartree_fock energy)"
--- a/src/non_h_ints_mu/jast_deriv_utils.irp.f
+++ b/src/non_h_ints_mu/jast_deriv_utils.irp.f
@ -99,6 +99,7 @@ subroutine grad1_j12_mu(r1, r2, grad)
    stop

  endif
+  grad = -grad

  return
 end subroutine grad1_j12_mu
@ -486,6 +487,13 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
    !!!!!!!!! rho1,rho2,rho1+rho2
    call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1)
    rho_tot = rho1 + rho2
+!    if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
+    if(rho_tot.lt.1.d-10)then
+     mu_val = mu_erf 
+     mu_der = 0.d0
+     return
+    endif
+    
    if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
    inv_rho_tot = 1.d0/rho_tot
    ! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf 
@ -506,18 +514,26 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
    ! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) ] (d rho(r1) / dx1)
    !
    ! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1)
-     
    !!!!!!!!! rho1,rho2,rho1+rho2
    call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1)
    rho_tot = rho1 + rho2
+!    if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
+    if(rho_tot.lt.1.d-10)then
+     mu_val = mu_erf 
+     mu_der = 0.d0
+     return
+    endif
+    
    if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
    inv_rho_tot = 1.d0/rho_tot
-    ! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf 
-    call get_all_f_rho_simple(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2)
+    ! f(rho) = (mu_r_ct* rho)**beta_rho_power * erf(zeta_erf_mu_of_r * rho) + mu_eff * (1 - erf(zeta_erf_mu_of_r*rho))
+    call get_all_f_rho_erf(rho1,rho2,mu_r_ct,beta_rho_power,mu_erf,zeta_erf_mu_of_r,f_rho1,d_drho_f_rho1,f_rho2)
    d_dx1_f_rho1(1:3)   = d_drho_f_rho1 * grad_rho1(1:3)
    d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3)
-    mu_val = 0.5d0 * ( f_rho1 + f_rho2)
-    mu_der(1:3) = d_dx_rho_f_rho(1:3) 
+    nume   = rho1 * f_rho1 + rho2 * f_rho2
+    mu_val = nume * inv_rho_tot
+    mu_der(1:3) = inv_rho_tot*inv_rho_tot * (rho_tot * d_dx_rho_f_rho(1:3) - grad_rho1(1:3) * nume)
+     
  else
    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
    stop
@ -676,8 +692,17 @@ subroutine get_all_f_rho_simple(rho1,rho2,alpha,mu0,beta,f_rho1,d_drho_f_rho1,f_
 double precision, intent(in) :: rho1,rho2,alpha,mu0,beta
 double precision, intent(out):: f_rho1,d_drho_f_rho1,f_rho2
 double precision :: tmp
- call f_mu_and_deriv_mu_simple(rho1,alpha,mu0,beta,f_rho1,d_drho_f_rho1)
- call f_mu_and_deriv_mu_simple(rho2,alpha,mu0,beta,f_rho2,tmp)
+ if(rho1.lt.1.d-10)then
+  f_rho1 = 0.d0
+  d_drho_f_rho1 = 0.d0
+ else
+  call f_mu_and_deriv_mu_simple(rho1,alpha,mu0,beta,f_rho1,d_drho_f_rho1)
+ endif
+ if(rho2.lt.1.d-10)then
+  f_rho2 = 0.d0
+ else
+  call f_mu_and_deriv_mu_simple(rho2,alpha,mu0,beta,f_rho2,tmp)
+ endif
 end

 subroutine f_mu_and_deriv_mu_simple(rho,alpha,mu0,beta,f_mu,d_drho_f_mu)
@ -691,10 +716,53 @@ subroutine f_mu_and_deriv_mu_simple(rho,alpha,mu0,beta,f_mu,d_drho_f_mu)
 END_DOC
 double precision, intent(in)  :: rho,alpha,mu0,beta
 double precision, intent(out) :: f_mu,d_drho_f_mu
- f_mu = alpha * (rho)**beta + mu0 
- d_drho_f_mu = alpha * beta * rho**(beta-1.d0) 
+ f_mu = alpha**beta * (rho)**beta + mu0 
+ d_drho_f_mu = alpha**beta * beta * rho**(beta-1.d0) 

 end

 ! ---

+subroutine f_mu_and_deriv_mu_erf(rho,alpha,zeta,mu0,beta,f_mu,d_drho_f_mu)
+ implicit none
+  include 'constants.include.F'
+ BEGIN_DOC
+! function giving mu as a function of rho
+!
+! f_mu = (alpha * rho)**zeta * erf(beta * rho) + mu0 * (1 - erf(beta*rho))
+!
+! and its derivative with respect to rho d_drho_f_mu
+!
+! d_drho_f_mu = 2 beta/sqrt(pi) * exp(-(beta*rho)**2) * ( (alpha*rho)**zeta - mu0) 
+!               + alpha * zeta * (alpha *rho)**(zeta-1)  * erf(beta*rho)
+ END_DOC
+ double precision, intent(in)  :: rho,alpha,mu0,beta,zeta
+ double precision, intent(out) :: f_mu,d_drho_f_mu
+ f_mu = (alpha * rho)**zeta * derf(beta * rho) + mu0 * (1.d0 - derf(beta*rho))
+ d_drho_f_mu = 2.d0  * beta * inv_sq_pi * dexp(-(beta*rho)**2) * ( (alpha*rho)**zeta - mu0) & 
+             + alpha * zeta * (alpha *rho)**(zeta-1)  * derf(beta*rho)
+
+end
+
+
+subroutine get_all_f_rho_erf(rho1,rho2,alpha,zeta,mu0,beta,f_rho1,d_drho_f_rho1,f_rho2)
+ implicit none
+ BEGIN_DOC
+! returns the values f(mu(r1)), f(mu(r2)) and d/drho(1) f(mu(r1))
+! with f_mu = (alpha * rho)**zeta * erf(beta * rho) + mu0 * (1 - erf(beta*rho))
+ END_DOC
+ double precision, intent(in) :: rho1,rho2,alpha,mu0,beta,zeta
+ double precision, intent(out):: f_rho1,d_drho_f_rho1,f_rho2
+ double precision :: tmp
+ if(rho1.lt.1.d-10)then
+  f_rho1 = mu_erf
+  d_drho_f_rho1 = 0.d0
+ else
+  call f_mu_and_deriv_mu_erf(rho1,alpha,zeta,mu0,beta,f_rho1,d_drho_f_rho1)
+ endif
+ if(rho2.lt.1.d-10)then
+  f_rho2 = mu_erf
+ else
+  call f_mu_and_deriv_mu_erf(rho2,alpha,zeta,mu0,beta,f_rho2,tmp)
+ endif
+end
--- a/src/non_h_ints_mu/plot_mu_of_r.irp.f
+++ b/src/non_h_ints_mu/plot_mu_of_r.irp.f
@ -13,9 +13,9 @@ subroutine routine_print
 integer :: i_unit_output,getUnitAndOpen
 output=trim(ezfio_filename)//'.mu_of_r'
 i_unit_output = getUnitAndOpen(output,'w')
- integer :: ipoint,nx
- double precision :: xmax,xmin,r(3),dx
- double precision :: mu_val, mu_der(3),dm_a,dm_b,grad
+ integer :: ipoint,nx,i
+ double precision :: xmax,xmin,r(3),dx,sigma
+ double precision :: mu_val, mu_der(3),dm_a,dm_b,grad,grad_dm_a(3), grad_dm_b(3)
 xmax =  5.D0
 xmin = -5.D0
 nx = 10000
@ -24,10 +24,15 @@ subroutine routine_print
 r(1) = xmin
 do ipoint = 1, nx
  call mu_r_val_and_grad(r, r, mu_val, mu_der)
-  call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
+  call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b)
+  sigma = 0.d0
+  do i = 1,3
+   sigma += grad_dm_a(i)**2
+  enddo
+  sigma=dsqrt(sigma)
  grad = mu_der(1)**2 + mu_der(2)**2 + mu_der(3)**2 
  grad = dsqrt(grad)
-  write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a+dm_b,grad
+  write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a+dm_b,grad,sigma/dm_a
  r(1) += dx
 enddo
 end
--- a/src/nuclei/write_pt_charges.py
+++ b/src/nuclei/write_pt_charges.py
@ -52,7 +52,7 @@ fncharges.write(" "+str(n_charges)+'\n')
 fncharges.close()
 mv_in_ezfio(EZFIO,tmp)

-# Write the file containing the charges and set in EZFIO folder 
+# Write the file containing the charges and set in EZFIO folder
 tmp="pts_charge_z"
 fcharges = open(tmp,'w')
 fcharges.write(" 1\n")
--- a/src/tc_keywords/EZFIO.cfg
+++ b/src/tc_keywords/EZFIO.cfg
@ -158,13 +158,19 @@ default: 0
 type: double precision
 doc: a parameter used to define mu(r)
 interface: ezfio, provider, ocaml
-default: 6.203504908994001e-1
+default: 1.5

 [beta_rho_power]
 type: double precision
 doc: a parameter used to define mu(r)
 interface: ezfio, provider, ocaml
-default: 0.5
+default: 0.33333
+
+[zeta_erf_mu_of_r]
+type: double precision
+doc: a parameter used to define mu(r)
+interface: ezfio, provider, ocaml
+default: 1.

 [thr_degen_tc]
 type: Threshold
--- a/src/tools/print_sorted_wf_coef.irp.f
+++ b/src/tools/print_sorted_wf_coef.irp.f
@ -13,7 +13,7 @@ subroutine routine
 output=trim(ezfio_filename)//'.wf_sorted'
 i_unit_output = getUnitAndOpen(output,'w')
 do i= 1, N_det
-  write(i_unit_output,*)i,dabs(psi_coef_sorted(i,1))
+  write(i_unit_output,*)i,dabs(psi_coef_sorted(i,1)),dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1))
 enddo

 end
--- a/src/utils/memory.irp.f
+++ b/src/utils/memory.irp.f
@ -33,6 +33,8 @@ subroutine resident_memory(value)
  call usleep(10)

  value = 0.d0
+IRP_IF MACOS
+IRP_ELSE
  iunit = getUnitAndOpen('/proc/self/status','r')
  do
    read(iunit,*,err=10,end=20) key, value
@ -43,6 +45,7 @@ subroutine resident_memory(value)
  end do
  20 continue
  close(iunit)
+IRP_ENDIF
  value = value / (1024.d0*1024.d0)
  call unlock_io()
 end function
@ -58,6 +61,9 @@ subroutine total_memory(value)
  double precision, intent(out) :: value

  call lock_io()
+  value = 0.d0
+IRP_IF MACOS
+IRP_ELSE
  iunit = getUnitAndOpen('/proc/self/status','r')
  do
    read(iunit,*,err=10,end=20) key, value
@ -68,6 +74,7 @@ subroutine total_memory(value)
  end do
  20 continue
  close(iunit)
+IRP_ENDIF
  value = value / (1024.d0*1024.d0)
  call unlock_io()
 end function