Merge pull request #11 from pfloos/work_env

Work env

PyDuck.py

@@ -1,13 +1,17 @@
 #!/usr/bin/env python
 
-import os
+import os, sys
 import argparse
 import pyscf
 from pyscf import gto
 import numpy as np
 import subprocess
 
 #Find the value of the environnement variable QUACK_ROOT. If not present we use the current repository
+if "QUACK_ROOT" not in os.environ:
+   print("Please set the QUACK_ROOT environment variable, for example:\n")
+   print("$ export QUACK_ROOT={0}".format(os.getcwd()))
+   sys.exit(1)
 QuAcK_dir=os.environ.get('QUACK_ROOT','./')
 
 #Create the argument parser object and gives a description of the script
@@ -37,7 +41,7 @@
 working_dir=args.working_dir
 
 #Read molecule
-f = open(QuAcK_dir + '/mol/' + xyz,'r')
+f = open(working_dir+'/mol/'+xyz,'r')
 lines = f.read().splitlines()
 nbAt = int(lines.pop(0))
 lines.pop(0)
@@ -70,6 +74,7 @@
 nalpha=nelec[0]
 nbeta=nelec[1]
 
+subprocess.call(['mkdir', '-p', working_dir+'/input'])
 f = open(working_dir+'/input/molecule','w')
 f.write('# nAt nEla nElb nCore nRyd\n')
 f.write(str(mol.natm)+' '+str(nalpha)+' '+str(nbeta)+' '+str(0)+' '+str(0)+'\n')
@@ -79,7 +84,8 @@
 f.close()
 
 #Compute nuclear energy and put it in a file
-subprocess.call(['rm', working_dir + '/int/ENuc.dat'])
+subprocess.call(['mkdir', '-p', working_dir+'/int'])
+subprocess.call(['rm', '-f', working_dir + '/int/ENuc.dat'])
 f = open(working_dir+'/int/ENuc.dat','w')
 f.write(str(mol.energy_nuc()))
 f.write(' ')
@@ -93,14 +99,14 @@
 x,y,z = dipole[0],dipole[1],dipole[2]
 
 norb = len(ovlp) # nBAS_AOs
-subprocess.call(['rm', working_dir + '/int/nBas.dat'])
+subprocess.call(['rm', '-f', working_dir + '/int/nBas.dat'])
 f = open(working_dir+'/int/nBas.dat','w')
 f.write(" {} ".format(str(norb)))
 f.close()
 
 
 def write_matrix_to_file(matrix,size,file,cutoff=1e-15):
-    f = open(file, 'a')
+    f = open(file, 'w')
     for i in range(size):
         for j in range(i,size):
             if abs(matrix[i][j]) > cutoff:
@@ -110,23 +116,23 @@ def write_matrix_to_file(matrix,size,file,cutoff=1e-15):
 
 #Write all 1 electron quantities in files
 #Ov,Nuc,Kin,x,y,z
-subprocess.call(['rm', working_dir + '/int/Ov.dat'])
+subprocess.call(['rm', '-f', working_dir + '/int/Ov.dat'])
 write_matrix_to_file(ovlp,norb,working_dir+'/int/Ov.dat')
-subprocess.call(['rm', working_dir + '/int/Nuc.dat'])
+subprocess.call(['rm', '-f', working_dir + '/int/Nuc.dat'])
 write_matrix_to_file(v1e,norb,working_dir+'/int/Nuc.dat')
-subprocess.call(['rm', working_dir + '/int/Kin.dat'])
+subprocess.call(['rm', '-f', working_dir + '/int/Kin.dat'])
 write_matrix_to_file(t1e,norb,working_dir+'/int/Kin.dat')
-subprocess.call(['rm', working_dir + '/int/x.dat'])
+subprocess.call(['rm', '-f', working_dir + '/int/x.dat'])
 write_matrix_to_file(x,norb,working_dir+'/int/x.dat')
-subprocess.call(['rm', working_dir + '/int/y.dat'])
+subprocess.call(['rm', '-f', working_dir + '/int/y.dat'])
 write_matrix_to_file(y,norb,working_dir+'/int/y.dat')
-subprocess.call(['rm', working_dir + '/int/z.dat'])
+subprocess.call(['rm', '-f', working_dir + '/int/z.dat'])
 write_matrix_to_file(z,norb,working_dir+'/int/z.dat')
 
 eri_ao = mol.intor('int2e')
 
 def write_tensor_to_file(tensor,size,file,cutoff=1e-15):
-    f = open(file, 'a')
+    f = open(file, 'w')
     for i in range(size):
         for j in range(i,size):
             for k in range(i,size):
@@ -140,15 +146,18 @@ def write_tensor_to_file(tensor,size,file,cutoff=1e-15):
 # Write two-electron integrals
 if print_2e:
     # (formatted)
-    subprocess.call(['rm', working_dir + '/int/ERI.dat'])
-    write_tensor_to_file(eri_ao,norb,working_dir+'/int/ERI.dat')
+    subprocess.call(['rm', '-f', working_dir + '/int/ERI.dat'])
+    write_tensor_to_file(eri_ao, norb, working_dir + '/int/ERI.dat')
 else:
     # (binary)
-    subprocess.call(['rm', working_dir + '/int/ERI.bin'])
+    subprocess.call(['rm', '-f', working_dir + '/int/ERI.bin'])
     # chem -> phys notation
     eri_ao = eri_ao.transpose(0, 2, 1, 3)
-    eri_ao.tofile('int/ERI.bin')
+    f = open(working_dir + '/int/ERI.bin', 'w')
+    eri_ao.tofile(working_dir + '/int/ERI.bin')
+    f.close()
 
 
 #Execute the QuAcK fortran program
-subprocess.call(QuAcK_dir+'/bin/QuAcK')
+subprocess.call([QuAcK_dir + '/bin/QuAcK', working_dir])
+

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_dipole_integrals(nBas,dipole_int_AO)
+  call read_integrals(working_dir,nBas,S,T,V,Hc,ERI_AO)
+  call read_dipole_integrals(working_dir,nBas,dipole_int_AO)
 
   call wall_time(end_int)