Repository: pyscf/pyscf
Issue: shci.dryrun() incorrectly printing FCIDUMP header

Example input:

```python
from pyscf import gto, scf, lo, tools, symm                                                                                                                                                   
from pyscf.shciscf import shci                                                                                                                                                                
import numpy                                                                                                                                                                                  
                                                                                                                                                                                              
mol = gto.M(atom = 'C     0.0000    1.396792    0.0000;\\                                                                                                                                      
C     0.0000    -1.396792    0.0000;\\                                                                                                                                                         
C     1.209657    0.698396    0.0000;\\                                                                                                                                                        
C     -1.209657    -0.698396    0.0000;\\                                                                                                                                                      
C    -1.209657    0.698396    0.0000;\\                                                                                                                                                        
C     1.209657    -0.698396    0.0000;\\                                                                                                                                                       
H     0.0000    2.484212    0.0000;\\                                                                                                                                                          
H     2.151390    1.242106    0.0000;\\                                                                                                                                                        
H     -2.151390    -1.242106    0.0000;\\                                                                                                                                                      
H     -2.151390    1.242106    0.0000;\\                                                                                                                                                       
H     2.151390    -1.242106    0.0000;\\                                                                                                                                                       
H     0.0000    -2.484212    0.0000',                                                                                                                                                         
basis = "ccpvdz",                                                                                                                                                                             
verbose=5, symmetry="Cs", spin=0) # When orbitals are localized the symmetry goes down from D6h to Cs. (Cs since sigma and pi do not mix)                                                     
                                                                                                                                                                                              
mf = scf.RHF(mol)                                                                                                                                                                             
mf.scf()                                                                                                                                                                                      
                                                                                                                                                                                              
nocc = mol.nelectron//2                                                                                                                                                                       
print (6, nocc, mf.mo_coeff.shape[1])                                                                                                                                                         
                                                                                                                                                                                              
# Use Pipek-Mezey localiztion.  We should reduce the symmetry from D2h to Cs when orbs are localized.  Cs, because sigma and pi do not mix.                                                   
loc_orbco = mf.mo_coeff[:,:6]                                                                                                                                                                 
loc_orbva = lo.PM(mol, mf.mo_coeff[:,6:nocc]).kernel()                                                                                                                                        
loc_orbvi = lo.PM(mol, mf.mo_coeff[:,nocc:]).kernel()                                                                                                                                         
lmo = numpy.hstack( (loc_orbco, loc_orbva, loc_orbvi) )                                                                                                                                       
                                                                                                                                                                                              
mc = shci.SHCISCF(mf, 108, 30)                                                                                                                                                                
mc.fcisolver.sweep_epsilon = [1.e-2]                                                                                                                                                          
mc.fcisolver.sweep_iter = [0]                                                                                                                                                                 
mc.fcisolver.epsilon2 = 1.e-3                                                                                                                                                                 
                                                                                                                                                                                              
shci.dryrun(mc, mf.mo_coeff)    
```

This script stalls at run time and the header for the FCIDUMP is messed up.

Here is the FCIDUMP header:

```
 &FCI NORB=108 ,NELEC=30 ,MS2=0,
  ORBSYM=2,539767840,840969268,876033073,875967288,539778359,1887007844,1886272869,1953392942,690565683,1920098570,1529379169,538976288,538976288,113,0,0,0,0,0,0,0,0,0,0,0,0,0,786590600,11030,157882997,0,786590565,11030,786590600,11030,157882997,0,786590554,11030,0,0,2481,0,1,0,-317304672,22080,2416,0,-1,-1,327908,327680,0,0,2004117858,1600484201,2015916655,2015374385,790637618,1970217004,1867398516,539780462,1998597162,1701995880,1970426941,1663052901,1769239393,658335598,1701667187,1852402527,539764580,1701081711,1260862834,1679830055,1701869940,1852788285,1931488357,1802461813,1970426941,539777893,1852270963,1920300129,1696607333,1651471480,170483050,1836008202,1702131056,1701344288,1953063456,1936291629,1380917349,543584032,544175988,1634890337,1696625529,1701668204,1999467630,778400617,1866664458,1953853549,1948283749,1646290280,1999467625,543519593,1864389199,1752440934,
  ISYM=1,
 &END
  3.669627509713e-01    1  1  1  1
  1.015442867114e-01    2  1  2  1
  3.400483496303e-01    2  2  1  1
  3.512403472134e-01    2  2  2  2
  1.015441919753e-01    3  1  3  1
  2.718127174536e-02    3  2  3  2

```

Changed files:
pyscf/dmrgscf/dmrgci.py, pyscf/shciscf/shci.py, pyscf/shciscf/test/test_shci.py

diff --git a/pyscf/dmrgscf/dmrgci.py b/pyscf/dmrgscf/dmrgci.py
index 85bc62ed0..699dbae2e 100644
--- a/pyscf/dmrgscf/dmrgci.py
+++ b/pyscf/dmrgscf/dmrgci.py
@@ -443,7 +443,7 @@ class DMRGCI(lib.StreamObject):
 
         return onepdm, twopdm, threepdm
 
-    def make_rdm3(self, state, norb, nelec, dt=numpy.dtype('Float64'), filetype = "binary", link_index=None, **kwargs):
+    def make_rdm3(self, state, norb, nelec, dt=numpy.float64, filetype = "binary", link_index=None, **kwargs):
         import os
 
         if self.has_threepdm == False:
@@ -488,7 +488,7 @@ class DMRGCI(lib.StreamObject):
             fname = os.path.join(self.scratchDirectory,"node0", "spatial_threepdm.%d.%d.bin" %(state, state))
             fnameout = os.path.join(self.scratchDirectory,"node0", "spatial_threepdm.%d.%d.bin.unpack" %(state, state))
             libunpack.unpackE3(ctypes.c_char_p(fname.encode()), ctypes.c_char_p(fnameout.encode()), ctypes.c_int(norb))
-            E3 = numpy.fromfile(fnameout, dtype=numpy.dtype('Float64'))
+            E3 = numpy.fromfile(fnameout, dtype=numpy.float64)
             E3 = numpy.reshape(E3, (norb, norb, norb, norb, norb, norb), order='F')
           else:
             print('Reading binary 3RDM from BLOCK')
@@ -518,7 +518,7 @@ class DMRGCI(lib.StreamObject):
         print('')
         return E3
 
-    def make_rdm4(self, state, norb, nelec, dt=numpy.dtype('Float64'), filetype = "binary", link_index=None, **kwargs):
+    def make_rdm4(self, state, norb, nelec, dt=numpy.float64, filetype = "binary", link_index=None, **kwargs):
         import os
 
         if self.has_fourpdm == False:
@@ -567,7 +567,7 @@ class DMRGCI(lib.StreamObject):
             fname = os.path.join(self.scratchDirectory,"node0", "spatial_fourpdm.%d.%d.bin" %(state, state))
             fnameout = os.path.join(self.scratchDirectory,"node0", "spatial_fourpdm.%d.%d.bin.unpack" %(state, state))
             libunpack.unpackE4(ctypes.c_char_p(fname.encode()), ctypes.c_char_p(fnameout.encode()), ctypes.c_int(norb))
-            E4 = numpy.fromfile(fnameout, dtype=numpy.dtype('Float64'))
+            E4 = numpy.fromfile(fnameout, dtype=numpy.float64)
             E4 = numpy.reshape(E4, (norb, norb, norb, norb, norb, norb, norb, norb), order='F')
           else:
             print('Reading binary 4RDM from BLOCK')
diff --git a/pyscf/shciscf/shci.py b/pyscf/shciscf/shci.py
index f1e54feac..7bb062ae0 100644
--- a/pyscf/shciscf/shci.py
+++ b/pyscf/shciscf/shci.py
@@ -464,7 +464,7 @@ class SHCI(pyscf.lib.StreamObject):
                   state,
                   norb,
                   nelec,
-                  dt=numpy.dtype('Float64'),
+                  dt=numpy.float64,
                   filetype="binary",
                   link_index=None,
                   bypass=False,
@@ -552,7 +552,7 @@ class SHCI(pyscf.lib.StreamObject):
                   state,
                   norb,
                   nelec,
-                  dt=numpy.dtype('Float64'),
+                  dt=numpy.float64,
                   filetype="binary",
                   link_index=None,
                   bypass=False,
@@ -1002,7 +1002,7 @@ def writeSHCIConfFile(SHCI, nelec, Restart):
     f.write('nroots %r\\n' % SHCI.nroots)
 
     # Variational Keyword Section
-    f.write('#variational\\n')
+    f.write('\\n#variational\\n')
     if (not Restart):
         schedStr = make_sched(SHCI)
         f.write(schedStr)
@@ -1022,7 +1022,7 @@ def writeSHCIConfFile(SHCI, nelec, Restart):
         f.write('fullrestart\\n')
 
     # Perturbative Keyword Section
-    f.write('#pt\\n')
+    f.write('\\n#pt\\n')
     if (SHCI.stochastic == False):
         f.write('deterministic \\n')
     else:
@@ -1033,8 +1033,8 @@ def writeSHCIConfFile(SHCI, nelec, Restart):
     f.write('sampleN %i\\n' % SHCI.sampleN)
 
     # Miscellaneous Keywords
-    f.write('#misc\\n')
-    f.write('noio \\n')
+    f.write('\\n#misc\\n')
+    f.write('io \\n')
     if (SHCI.scratchDirectory != ""):
         if not os.path.exists(SHCI.scratchDirectory):
             os.makedirs(SHCI.scratchDirectory)
@@ -1480,12 +1480,41 @@ def writeSOCIntegrals(mc,
 
 
 def dryrun(mc, mo_coeff=None):
-    '''Generate FCIDUMP and SHCI config file'''
+    '''
+    Generate FCIDUMP and SHCI input.dat file for a give multiconfigurational 
+    object. This method DOES NOT run a single iteration using Dice it just 
+    generates all the inputs you'd need to do so.
+
+    Args:
+        mc: a CASSCF/CASCI object (or RHF object)
+
+    Kwargs: 
+        mo_coeff: `np.ndarray` of the molecular orbital coefficients. 
+            Default is None.
+
+    Examples:
+
+    >>> from pyscf import gto, scf, mcscf
+    >>> from pyscf.shciscf import shci
+    >>> mol = gto.M(atom="H 0 0 0; H 0 0 2; O 0 1 1;", symmetry=True)
+    >>> mf = scf.RHF(mol).run()
+    >>> mc = mcscf.CASCI(mf, 3, 4)
+    >>> mc.fcisolver = shci.SHCI(mol)
+    >>> shci.dryrun(mc)
+
+    '''
     if mo_coeff is None:
         mo_coeff = mc.mo_coeff
-    #bak, mc.fcisolver.onlywriteIntegral = mc.fcisolver.onlywriteIntegral, True
-    mc.casci(mo_coeff)
-    #mc.fcisolver.onlywriteIntegral = bak
+    
+    # Set orbsym b/c it's needed in `writeIntegralFile()`
+    if mc.fcisolver.orbsym is not []:
+        mc.fcisolver.orbsym = getattr(mo_coeff, "orbsym", [])
+
+    #mc.kernel(mo_coeff) # Works, but runs full CASCI/CASSCF
+    h1e, ecore = mc.get_h1eff(mo_coeff)
+    h2e = mc.get_h2eff(mo_coeff)
+    writeIntegralFile(mc.fcisolver, h1e, h2e, mc.ncas, mc.nelecas, ecore)
+    writeSHCIConfFile(mc.fcisolver, mc.nelecas, False)
 
 
 #mc is the CASSCF object
diff --git a/pyscf/shciscf/test/test_shci.py b/pyscf/shciscf/test/test_shci.py
index 620d7abcc..a1f337ce7 100644
--- a/pyscf/shciscf/test/test_shci.py
+++ b/pyscf/shciscf/test/test_shci.py
@@ -176,7 +176,6 @@ class KnownValues(unittest.TestCase):
         # Number of orbital and electrons
         ncas = 8
         nelecas = 12
-        dimer_atom = 'O'
 
         mc = mcscf.CASCI(mf, ncas, nelecas)
         e_casscf = mc.kernel()[0]
@@ -199,8 +198,7 @@ class KnownValues(unittest.TestCase):
         # Number of orbital and electrons
         ncas = 8
         nelecas = 12
-        dimer_atom = 'O'
-
+        
         mc = mcscf.CASSCF(mf, ncas, nelecas)
         e_casscf = mc.kernel()[0]