absolute PATH in Fortran I/O

PyDuck.py

@@ -155,5 +155,5 @@ else:
 
 
 #Execute the QuAcK fortran program
-print(QuAcK_dir)
+subprocess.call([QuAcK_dir + '/bin/QuAcK', working_dir])
-subprocess.call(QuAcK_dir+'/bin/QuAcK')
+

src/QuAcK/QuAcK.f90

@@ -71,6 +71,16 @@ program QuAcK
 
   logical                       :: dotest,doRtest,doUtest,doGtest
 
+  character(len=256)            :: working_dir
+
+  ! Check if the right number of arguments is provided
+  if(command_argument_count() < 1) then
+    print *, "No working directory provided."
+    stop
+  else
+    call get_command_argument(1, working_dir)
+  endif
+
 !-------------!
 ! Hello World !
 !-------------!
@@ -95,7 +105,8 @@ program QuAcK
 ! Method selection !
 !------------------!
 
-  call read_methods(doRHF,doUHF,doGHF,doROHF,              &
+  call read_methods(working_dir,                           &
+                    doRHF,doUHF,doGHF,doROHF,              &
                     doMP2,doMP3,                           &
                     doCCD,dopCCD,doDCD,doCCSD,doCCSDT,     &
                     dodrCCD,dorCCD,docrCCD,dolCCD,         &
@@ -112,7 +123,8 @@ program QuAcK
 ! Read options for methods !
 !--------------------------!
 
-  call read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, &
+  call read_options(working_dir,                                                                &
+                    maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, &
                     reg_MP,                                                                     &
                     maxSCF_CC,thresh_CC,max_diis_CC,                                            &
                     TDA,spin_conserved,spin_flip,                                               &
@@ -133,18 +145,18 @@ program QuAcK
 ! nOrb = number of orbitals          !
 !------------------------------------!
 
-  call read_molecule(nNuc,nO,nC,nR)
+  call read_molecule(working_dir,nNuc,nO,nC,nR)
   allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
 
 ! Read geometry
 
-  call read_geometry(nNuc,ZNuc,rNuc,ENuc)
+  call read_geometry(working_dir,nNuc,ZNuc,rNuc,ENuc)
 
 !---------------------------------------!
 ! Read basis set information from PySCF !
 !---------------------------------------!
 
-  call read_basis_pyscf(nBas,nO,nV)
+  call read_basis_pyscf(working_dir,nBas,nO,nV)
 
 !--------------------------------------!
 ! Read one- and two-electron integrals !
@@ -163,8 +175,8 @@ program QuAcK
 
   call wall_time(start_int)
 
-  call read_integrals(nBas,S,T,V,Hc,ERI_AO)
+  call read_integrals(working_dir,nBas,S,T,V,Hc,ERI_AO)
-  call read_dipole_integrals(nBas,dipole_int_AO)
+  call read_dipole_integrals(working_dir,nBas,dipole_int_AO)
 
   call wall_time(end_int)
 

src/QuAcK/read_methods.f90

@@ -1,4 +1,5 @@
-subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
+subroutine read_methods(working_dir,                           &
+                        doRHF,doUHF,doGHF,doROHF,              &
                         doMP2,doMP3,                           & 
                         doCCD,dopCCD,doDCD,doCCSD,doCCSDT,     & 
                         do_drCCD,do_rCCD,do_crCCD,do_lCCD,     &
@@ -17,6 +18,10 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
 
 ! Input variables
 
+  character(len=256),intent(in) :: working_dir
+
+! Output variables
+
   logical,intent(out)           :: doRHF,doUHF,doGHF,doROHF
   logical,intent(out)           :: doMP2,doMP3
   logical,intent(out)           :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
@@ -33,12 +38,21 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
 ! Local variables
 
   character(len=1)              :: ans1,ans2,ans3,ans4,ans5,ans6
+  integer                       :: status
+  character(len=256)            :: file_path
 
-! Open file with method specification
 
-  open(unit=1,file='input/methods')
+  file_path = trim(working_dir) // '/input/methods'
+  open(unit=1, file=file_path, status='old', action='read', iostat=status)
 
-! Read mean-field methods
+    if(status /= 0) then
+
+      print *, "Error opening file: ", file_path
+      stop
+
+    else
+
+      ! Read mean-field methods
       
       doRHF  = .false.
       doUHF  = .false.
@@ -52,7 +66,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans3 == 'T') doGHF  = .true.
       if(ans4 == 'T') doROHF = .true.
       
-! Read MPn methods
+      ! Read MPn methods
       
       doMP2    = .false.
       doMP3    = .false.
@@ -62,7 +76,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans1 == 'T') doMP2    = .true.
       if(ans2 == 'T') doMP3    = .true.
       
-! Read CC methods
+      ! Read CC methods
       
       doCCD   = .false.
       dopCCD  = .false.
@@ -78,7 +92,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans4 == 'T') doCCSD  = .true.
       if(ans5 == 'T') doCCSDT = .true.
       
-! Read weird CC methods
+      ! Read weird CC methods
       
       do_drCCD = .false.
       do_rCCD  = .false.
@@ -92,7 +106,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans3 == 'T') do_crCCD = .true.
       if(ans4 == 'T') do_lCCD  = .true.
       
-! Read CI methods
+      ! Read CI methods
       
       doCIS   = .false.
       doCIS_D = .false.
@@ -109,7 +123,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans5 == 'T') doFCI   = .true.
       if(doCIS_D)        doCIS   = .true.
       
-! Read RPA methods
+      ! Read RPA methods
       
       dophRPA  = .false.
       dophRPAx = .false.
@@ -123,7 +137,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans3 == 'T') docrRPA  = .true.
       if(ans4 == 'T') doppRPA  = .true.
       
-! Read Green's function methods
+      ! Read Green's function methods
       
       doG0F2    = .false.
       doevGF2   = .false.
@@ -141,7 +155,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans5 == 'T') doG0F3    = .true.
       if(ans6 == 'T') doevGF3   = .true.
       
-! Read GW methods
+      ! Read GW methods
       
       doG0W0    = .false.
       doevGW    = .false.
@@ -157,7 +171,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans4 == 'T') doufG0W0  = .true.
       if(ans5 == 'T') doufGW    = .true.
       
-! Read GTpp methods
+      ! Read GTpp methods
       
       doG0T0pp   = .false.
       doevGTpp   = .false.
@@ -171,7 +185,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans3 == 'T') doqsGTpp   = .true.
       if(ans4 == 'T') doufG0T0pp = .true.
       
-! Read GTeh methods
+      ! Read GTeh methods
       
       doG0T0eh   = .false.
       doevGTeh   = .false.
@@ -183,7 +197,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans2 == 'T') doevGTeh = .true.
       if(ans3 == 'T') doqsGTeh = .true.
       
-! Read test
+      ! Read test
       
       doRtest = .false.
       doUtest = .false.
@@ -195,8 +209,9 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
       if(ans2 == 'T') doUtest = .true.
       if(ans3 == 'T') doGtest = .true.
 
-! Close file 
+    endif
 
+  ! Close file
   close(unit=1)
 
 end subroutine 

src/QuAcK/read_options.f90

@@ -1,4 +1,5 @@
-subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, &
+subroutine read_options(working_dir,                                                                &
+                        maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, &
                         reg_MP,                                                                     &
                         maxSCF_CC,thresh_CC,max_diis_CC,                                            &
                         TDA,spin_conserved,spin_flip,                                               &
@@ -14,6 +15,10 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
 
 ! Input variables
 
+  character(len=256),intent(in) :: working_dir
+
+! Output variables
+
   integer,intent(out)           :: maxSCF_HF
   double precision,intent(out)  :: thresh_HF
   integer,intent(out)           :: max_diis_HF
@@ -70,12 +75,22 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
 ! Local variables
 
   character(len=1)              :: ans1,ans2,ans3,ans4,ans5
+  integer                       :: status
+  character(len=256)            :: file_path
 
 ! Open file with method specification
 
-  open(unit=1,file='input/options')
+  file_path = trim(working_dir) // '/input/options'
+  open(unit=1, file=file_path, status='old', action='read', iostat=status)
 
-! Read HF options
+    if(status /= 0) then
+
+      print *, "Error opening file: ", file_path
+      stop
+
+    else
+
+      ! Read HF options
     
       maxSCF_HF    = 64
       thresh_HF    = 1d-6
@@ -92,7 +107,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
       if(ans1 == 'T') dostab   = .true.
       if(ans2 == 'T') dosearch = .true.
     
-! Read MPn options
+      ! Read MPn options
     
       reg_MP = .false.
       read(1,*)
@@ -100,7 +115,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
     
       if(ans1 == 'T') reg_MP = .true.
     
-! Read CC options
+      ! Read CC options
     
       maxSCF_CC    = 64
       thresh_CC    = 1d-5
@@ -109,7 +124,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
       read(1,*)
       read(1,*) maxSCF_CC,thresh_CC,max_diis_CC
     
-! Read excited state options
+      ! Read excited state options
     
       TDA            = .false.
       spin_conserved = .false.
@@ -122,7 +137,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
       if(ans2 == 'T') spin_conserved = .true.
       if(ans3 == 'T') spin_flip      = .true.
     
-! Read GF options
+      ! Read GF options
     
       maxSCF_GF   = 64
       thresh_GF   = 1d-5
@@ -138,7 +153,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
       if(ans1 == 'T') lin_GF  = .true.
       if(ans2 == 'T') reg_GF  = .true.
     
-! Read GW options
+      ! Read GW options
     
       maxSCF_GW   = 64
       thresh_GW   = 1d-5
@@ -155,7 +170,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
       if(ans2 == 'T') TDA_W   = .true.
       if(ans3 == 'T') reg_GW  = .true.
     
-! Read GT options
+      ! Read GT options
     
       maxSCF_GT   = 64
       thresh_GT   = 1d-5
@@ -172,7 +187,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
       if(ans2 == 'T') TDA_T   = .true.
       if(ans3 == 'T') reg_GT  = .true.
     
-! Options for adiabatic connection
+      ! Options for adiabatic connection
     
       doACFDT         = .false.
       exchange_kernel = .false.
@@ -185,7 +200,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
       if(ans2 == 'T') exchange_kernel = .true.
       if(ans3 == 'T') doXBS   = .true.
     
-! Options for dynamical BSE calculations
+      ! Options for dynamical BSE calculations
     
       dophBSE  = .false.
       dophBSE2 = .false.
@@ -202,8 +217,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
       if(ans4 == 'T') dBSE     = .true.
       if(ans5 == 'F') dTDA     = .false.
 
-! Close file with options
+    endif
 
+  ! Close file with options
   close(unit=1)
 
 end subroutine 

src/utils/read_basis_pyscf.f90

@@ -1,4 +1,4 @@
-subroutine read_basis_pyscf(nBas,nO,nV)
+subroutine read_basis_pyscf(working_dir,nBas,nO,nV)
 
 ! Read basis set information from PySCF
 
@@ -8,20 +8,36 @@ subroutine read_basis_pyscf(nBas,nO,nV)
 ! Input variables
 
   integer,intent(in)            :: nO(nspin)
-
+  character(len=256),intent(in) :: working_dir
-! Local variables
 
 ! Output variables
 
   integer,intent(out)            :: nV(nspin)
   integer,intent(out)            :: nBas
 
+! Local variables
+
+  integer                       :: status
+  character(len=256)            :: file_path
+
 !------------------------------------------------------------------------
 ! Primary basis set information
 !------------------------------------------------------------------------
 
-  open(unit=3,file='int/nBas.dat')
+  file_path = trim(working_dir) // '/int/nBas.dat'
+  open(unit=3, file=file_path, status='old', action='read', iostat=status)
+
+    if(status /= 0) then
+
+      print *, "Error opening file: ", file_path
+      stop
+
+    else
+
       read(3, *) nBas
+
+    endif
+
   close(unit=3)
 
 !  write(*,'(A38)') '--------------------------------------'

src/utils/read_dipole_integrals.f90

@@ -1,4 +1,4 @@
-subroutine read_dipole_integrals(nBas,R)
+subroutine read_dipole_integrals(working_dir,nBas,R)
 
 ! Read one-electron integrals related to dipole moment from files
 
@@ -8,6 +8,7 @@ subroutine read_dipole_integrals(nBas,R)
 ! Input variables
 
   integer,intent(in)            :: nBas
+  character(len=256),intent(in) :: working_dir
 
 ! Local variables
 
@@ -19,36 +20,83 @@ subroutine read_dipole_integrals(nBas,R)
 
   double precision,intent(out)  :: R(nBas,nBas,ncart)
 
+  integer                       :: status, ios
+  character(len=256)            :: file_path
+
+
 ! Open file with integrals
 
-  open(unit=21,file='int/x.dat')
-  open(unit=22,file='int/y.dat')
-  open(unit=23,file='int/z.dat')
-
-! Read (x,y,z) integrals
-
   R(:,:,:) = 0d0
 
+  file_path = trim(working_dir) // '/int/x.dat'
+  open(unit=21, file=file_path, status='old', action='read', iostat=status)
+
+    if(status /= 0) then
+
+      print *, "Error opening file: ", file_path
+      stop
+  
+    else
+
       do
-    read(21,*,end=21) mu,nu,Dip
+        read(21,*,iostat=ios) mu,nu,Dip
+        if(ios /= 0) exit
         R(mu,nu,1) = Dip
         R(nu,mu,1) = Dip
       end do
-  21 close(unit=21)
+
+    endif
+
+  close(unit=21)
+
+  ! ---
+
+  file_path = trim(working_dir) // '/int/y.dat'
+  open(unit=22, file=file_path, status='old', action='read', iostat=status)
+
+    if(status /= 0) then
+
+      print *, "Error opening file: ", file_path
+      stop
+  
+    else
 
       do
-    read(22,*,end=22) mu,nu,Dip
+        read(22,*,iostat=ios) mu,nu,Dip
+        if(ios /= 0) exit
         R(mu,nu,2) = Dip
         R(nu,mu,2) = Dip
       end do
-  22 close(unit=22)
+
+    endif
+
+  close(unit=22)
+
+  ! ---
+
+  file_path = trim(working_dir) // '/int/z.dat'
+  open(unit=23, file=file_path, status='old', action='read', iostat=status)
+
+    if(status /= 0) then
+
+      print *, "Error opening file: ", file_path
+      stop
+  
+    else
 
       do
-    read(23,*,end=23) mu,nu,Dip
+        read(23,*,iostat=ios) mu,nu,Dip
+        if(ios /= 0) exit
         R(mu,nu,3) = Dip
         R(nu,mu,3) = Dip
       end do
-  23 close(unit=23)
+
+    endif
+
+  close(unit=23)
+
+  ! ---
+
 
 ! Print results
   if(debug) then

src/utils/read_geometry.f90

@@ -1,10 +1,14 @@
-subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc)
+subroutine read_geometry(working_dir,nNuc,ZNuc,rNuc,ENuc)
 
 ! Read molecular geometry
 
   implicit none
   include 'parameters.h'
 
+! Input variables
+
+  character(len=256),intent(in) :: working_dir
+
 ! Ouput variables
 
   integer,intent(in)            :: nNuc
@@ -20,12 +24,23 @@ subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc)
 
   double precision,intent(out)  :: ZNuc(nNuc),rNuc(nNuc,ncart),ENuc
 
+  integer                       :: status
+  character(len=256)            :: file_path
+
 ! Open file with geometry specification
 
-  open(unit=10,file='input/molecule')
+  file_path = trim(working_dir) // '/input/molecule'
-  open(unit=11,file='input/molecule.xyz')
+  open(unit=10, file=file_path, status='old', action='read', iostat=status)
 
-! Read geometry and create xyz file for integrals
+    if(status /= 0) then
+
+      print *, "Error opening file: ", file_path
+      stop
+
+    else
+
+      ! Read geometry and create xyz file for integrals
+      open(unit=11,file=trim(working_dir) // '/input/molecule.xyz')
       
       read(10,*) 
       read(10,*) 
@@ -40,23 +55,39 @@ subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc)
         ZNuc(i) = dble(element_number(El))
       end do
 
-! Compute nuclear repulsion energy
+      close(unit=11)
+
+    endif
+
+  close(unit=10)
+
+  ! ---
+
+  file_path = trim(working_dir) // '/int/ENuc.dat'
+  open(unit=3, file=file_path, status='old', action='read', iostat=status)
+
+    if(status /= 0) then
+
+      print *, "Error opening file: ", file_path
+      stop
+
+    else
 
-  ENuc = 0
-  open(unit=3,file='int/ENuc.dat')
       read(3,*) ENuc
+
+    endif
+
   close(unit=3)
       
-  ! do i=1,nNuc-1
+  ! Compute nuclear repulsion energy
+  !ENuc = 0.d0
+  !do i=1,nNuc-1
   !  do j=i+1,nNuc
   !    RAB = (rNuc(i,1)-rNuc(j,1))**2 + (rNuc(i,2)-rNuc(j,2))**2 + (rNuc(i,3)-rNuc(j,3))**2
   !    ENuc = ENuc + ZNuc(i)*ZNuc(j)/(AntoBo*sqrt(RAB))
   !  end do
-  ! end do
+  !end do
 
-! Close file with geometry specification
-  close(unit=10)
-  close(unit=11)
 
 ! Print geometry
   write(*,'(A28)') '------------------'

src/utils/read_integrals.f90

@@ -1,4 +1,4 @@
-subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
+subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G)
 
 ! Read one- and two-electron integrals from files
 
@@ -8,6 +8,7 @@ subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
 ! Input variables
 
   integer,intent(in)            :: nBas_AOs
+  character(len=256),intent(in) :: working_dir
 
 ! Local variables
 
@@ -24,6 +25,9 @@ subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
   double precision,intent(out)  :: Hc(nBas_AOs,nBas_AOs)
   double precision,intent(out)  :: G(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
 
+  integer                       :: status, ios
+  character(len=256)            :: file_path
+
 ! Open file with integrals
 
   debug = .false.
@@ -32,46 +36,67 @@ subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
 
   print*, 'Scaling integrals by ',lambda
 
-  open(unit=8 ,file='int/Ov.dat')
-  open(unit=9 ,file='int/Kin.dat')
-  open(unit=10,file='int/Nuc.dat')
 
-  open(unit=21,file='int/x.dat')
+  ! ---
-  open(unit=22,file='int/y.dat')
-  open(unit=23,file='int/z.dat')
-
-! Read overlap integrals
 
+  ! Read overlap integrals
+  file_path = trim(working_dir) // '/int/Ov.dat'
+  open(unit=8, file=file_path, status='old', action='read', iostat=status)
+    if(status /= 0) then
+      print *, "Error opening file: ", file_path
+      stop
+    else
       S(:,:) = 0d0
       do 
-    read(8,*,end=8) mu,nu,Ov
+        read(8,*,iostat=ios) mu,nu,Ov
+        if(ios /= 0) exit
         S(mu,nu) = Ov
         S(nu,mu) = Ov
       end do
-  8 close(unit=8)
+    endif
+  close(unit=8)
 
-! Read kinetic integrals
+  ! ---
 
+  ! Read kinetic integrals
+  file_path = trim(working_dir) // '/int/Kin.dat'
+  open(unit=9, file=file_path, status='old', action='read', iostat=status)
+    if(status /= 0) then
+      print *, "Error opening file: ", file_path
+      stop
+    else
       T(:,:) = 0d0
       do 
-    read(9,*,end=9) mu,nu,Kin
+        read(9,*,iostat=ios) mu,nu,Kin
+        if(ios /= 0) exit
         T(mu,nu) = Kin
         T(nu,mu) = Kin
       end do
-  9 close(unit=9)
+    endif
+  close(unit=9)
 
-! Read nuclear integrals
+  ! ---
 
+  ! Read nuclear integrals
+  file_path = trim(working_dir) // '/int/Nuc.dat'
+  open(unit=10, file=file_path, status='old', action='read', iostat=status)
+    if(status /= 0) then
+      print *, "Error opening file: ", file_path
+      stop
+    else
       V(:,:) = 0d0
       do 
-    read(10,*,end=10) mu,nu,Nuc
+        read(10,*,iostat=ios) mu,nu,Nuc
+        if(ios /= 0) exit
         V(mu,nu) = Nuc
         V(nu,mu) = Nuc
       end do
-  10 close(unit=10)
+    endif
+  close(unit=10)
 
-! Define core Hamiltonian
+  ! ---
 
+  ! Define core Hamiltonian
   Hc(:,:) = T(:,:) + V(:,:)
 
 ! Read 2e-integrals
@@ -94,9 +119,16 @@ subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
 !  11 close(unit=11)
 
   ! binary file
-  open(unit=11, file='int/ERI.bin', form='unformatted', access='stream')
+  file_path = trim(working_dir) // '/int/ERI.bin'
+  open(unit=11, file=file_path, status='old', action='read', form='unformatted', access='stream', iostat=status)
+    if(status /= 0) then
+      print *, "Error opening file: ", file_path
+      stop
+    else
       read(11) G
-  close(11)
+    endif
+  close(unit=11)
+
 
 
 ! Print results

src/utils/read_molecule.f90

@@ -1,4 +1,4 @@
-subroutine read_molecule(nNuc,nO,nC,nR)
+subroutine read_molecule(working_dir,nNuc,nO,nC,nR)
 
 ! Read number of atoms and number of electrons 
 
@@ -11,6 +11,10 @@ subroutine read_molecule(nNuc,nO,nC,nR)
   integer                       :: nCore
   integer                       :: nRyd
 
+! Input variables
+
+  character(len=256),intent(in) :: working_dir
+
 ! Output variables
 
   integer,intent(out)           :: nNuc
@@ -18,11 +22,22 @@ subroutine read_molecule(nNuc,nO,nC,nR)
   integer,intent(out)           :: nC(nspin)
   integer,intent(out)           :: nR(nspin)
 
+  integer                       :: status
+  character(len=256)            :: file_path
+
 ! Open file with geometry specification
 
-  open(unit=1,file='input/molecule')
+  file_path = trim(working_dir) // '/input/molecule'
+  open(unit=1, file=file_path, status='old', action='read', iostat=status)
 
-! Read number of atoms and number of electrons
+    if(status /= 0) then
+
+      print *, "Error opening file: ", file_path
+      stop
+
+    else
+
+      ! Read number of atoms and number of electrons
     
       read(1,*) 
       read(1,*) nNuc,nO(1),nO(2),nCore,nRyd
@@ -37,7 +52,7 @@ subroutine read_molecule(nNuc,nO,nC,nR)
       nC(:) = nCore/2
       nR(:) = nRyd/2
     
-! Print results
+      ! Print results
     
       write(*,'(A28)') '----------------------'
       write(*,'(A28,1X,I16)') 'Number of atoms',nNuc
@@ -55,8 +70,9 @@ subroutine read_molecule(nNuc,nO,nC,nR)
       write(*,'(A28)') '----------------------'
       write(*,*)
 
-! Close file with geometry specification
+    endif
 
+  ! Close file with geometry specification
   close(unit=1)
 
 end subroutine