From 701544903135005f3730752f24732947f678c25c Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 23 Mar 2015 17:18:48 +0100
Subject: [PATCH] Beginning work on CASINO interface

---
 src/Dets/save_for_casino.irp.f | 183 +++++++++++++++++++++++++++++++++
 1 file changed, 183 insertions(+)
 create mode 100644 src/Dets/save_for_casino.irp.f

diff --git a/src/Dets/save_for_casino.irp.f b/src/Dets/save_for_casino.irp.f
new file mode 100644
index 00000000..0c18d1e2
--- /dev/null
+++ b/src/Dets/save_for_casino.irp.f
@@ -0,0 +1,183 @@
+subroutine save_casino
+ use bitmasks
+ implicit none
+ character*(128) :: message
+ integer                        :: getUnitAndOpen, iunit
+ integer, allocatable           :: itmp(:)
+ real, allocatable              :: rtmp(:)
+ PROVIDE ezfio_filename
+
+ iunit = getUnitAndOpen('gwfn.data','w')
+ print *, 'Title?'
+ read(*,*) message
+ write(iunit,'(A)') trim(message)
+ write(iunit,'(A)') ''
+ write(iunit,'(A)') 'BASIC_INFO'
+ write(iunit,'(A)') '----------'
+ write(iunit,'(A)') 'Generated by:'
+ write(iunit,'(A)') 'Quantum package'
+ write(iunit,'(A)') 'Method:'
+ print *, 'Method?'
+ read(*,*) message
+ write(iunit,'(A)') trim(message)
+ write(iunit,'(A)') 'DFT Functional:'
+ write(iunit,'(A)') 'none'
+ write(iunit,'(A)') 'Periodicity:'
+ write(iunit,'(A)') '0'
+ write(iunit,'(A)') 'Spin unrestricted:'
+ write(iunit,'(A)') '.false.'
+ write(iunit,'(A)') 'nuclear-nuclear repulsion energy (au/atom):'
+ write(iunit,*) nuclear_repulsion
+ write(iunit,'(A)') 'Number of electrons per primitive cell:'
+ write(iunit,*) elec_num
+ write(iunit,*) ''
+
+
+ write(iunit,*) 'GEOMETRY'
+ write(iunit,'(A)') '--------'
+ write(iunit,'(A)') 'Number of atoms:'
+ write(iunit,*) nucl_num
+ write(iunit,'(A)') 'Atomic positions (au):'
+ integer                        :: i
+ do i=1,nucl_num
+   write(iunit,'(3(1PE20.13))') nucl_coord(i,1:3)
+ enddo
+ write(iunit,'(A)') 'Atomic numbers for each atom:'
+ ! Add 200 if pseudopotential
+ allocate(itmp(nucl_num))
+ do i=1,nucl_num
+   itmp(i) = int(nucl_charge(i))
+ enddo
+ write(iunit,'(8(I10))') itmp(1:nucl_num)
+ deallocate(itmp)
+ write(iunit,'(A)') 'Valence charges for each atom:'
+ write(iunit,'(4(1PE20.13))') nucl_charge(1:nucl_num)
+ write(iunit,'(A)') ''
+
+
+ write(iunit,'(A)') 'BASIS SET'
+ write(iunit,'(A)') '---------'
+ write(iunit,'(A)') 'Number of Gaussian centres'
+ write(iunit,*) nucl_num
+ write(iunit,'(A)') 'Number of shells per primitive cell'
+ integer :: icount
+ icount = 0
+ do i=1,ao_num
+  if (ao_l(i) == ao_power(i,1)) then
+    icount += 1
+  endif
+ enddo
+ write(iunit,*) icount
+ write(iunit,'(A)') 'Number of basis functions (''AO'') per primitive cell'
+ icount = 0
+ do i=1,ao_num
+  if (ao_l(i) == ao_power(i,1)) then
+    icount += 2*ao_l(i)+1
+  endif
+ enddo
+ write(iunit,*) icount
+ write(iunit,'(A)') 'Number of Gaussian primitives per primitive cell'
+ allocate(itmp(ao_num))
+ integer :: l
+ l=0
+ do i=1,ao_num
+   if (ao_l(i) == ao_power(i,1)) then
+     l += 1
+     itmp(l) = ao_prim_num(i)
+   endif
+ enddo
+ write(iunit,'(8(I10))') sum(itmp(1:l))
+ write(iunit,'(A)') 'Highest shell angular momentum (s/p/d/f... 1/2/3/4...)'
+ write(iunit,*) maxval(ao_l(1:ao_num))+1
+ write(iunit,'(A)') 'Code for shell types (s/sp/p/d/f... 1/2/3/4/5...)'
+ l=0
+ do i=1,ao_num
+   if (ao_l(i) == ao_power(i,1)) then
+     l += 1
+     if (ao_l(i) > 0) then
+       itmp(l) = ao_l(i)+2
+     else
+       itmp(l) = ao_l(i)+1
+     endif
+   endif
+ enddo
+ write(iunit,'(8(I10))') itmp(1:l)
+ write(iunit,'(A)') 'Number of primitive Gaussians in each shell'
+ l=0
+ do i=1,ao_num
+   if (ao_l(i) == ao_power(i,1)) then
+     l += 1
+     itmp(l) = ao_prim_num(i)
+   endif
+ enddo
+ write(iunit,'(8(I10))') itmp(1:l)
+ deallocate(itmp)
+ write(iunit,'(A)') 'Sequence number of first shell on each centre'
+ allocate(itmp(nucl_num))
+ l=0
+ icount = 1
+ itmp(icount) = 1
+ do i=1,ao_num
+  if (ao_l(i) == ao_power(i,1)) then
+    l = l+1
+    if (ao_nucl(i) == icount) then
+      continue
+    else if (ao_nucl(i) == icount+1) then
+      icount += 1
+      itmp(icount) = l
+    else
+      print *, 'Problem in order of centers of basis functions'
+      stop 1
+    endif
+  endif
+ enddo
+ ! Check
+ if (icount /= nucl_num) then
+   print *,  'Error :'
+   print *,  ' icount :', icount
+   print *,  ' nucl_num:', nucl_num
+   stop 2
+ endif
+ write(iunit,'(8(I10))') itmp(1:nucl_num)
+ deallocate(itmp)
+ write(iunit,'(A)') 'Exponents of Gaussian primitives'
+ allocate(rtmp(ao_num))
+ l=0
+ do i=1,ao_num
+   if (ao_l(i) == ao_power(i,1)) then
+     do j=1,ao_prim_num(i)
+       l+=1
+       rtmp(l) = ao_expo(i,ao_prim_num(i)-j+1)
+     enddo
+   endif
+ enddo
+ write(iunit,'(4(1PE20.13))') rtmp(1:l)
+ write(iunit,'(A)') 'Normalized contraction coefficients'
+ l=0
+ integer :: j
+ do i=1,ao_num
+   if (ao_l(i) == ao_power(i,1)) then
+     do j=1,ao_prim_num(i)
+       l+=1
+       rtmp(l) = ao_coef(i,ao_prim_num(i)-j+1)
+     enddo
+   endif
+ enddo
+ write(iunit,'(4(1PE20.13))') rtmp(1:l)
+ deallocate(rtmp)
+ write(iunit,'(A)') 'Position of each shell (au)'
+ l=0
+ do i=1,ao_num
+   if (ao_l(i) == ao_power(i,1)) then
+     write(iunit,'(3(1PE20.13))') nucl_coord( ao_nucl(i), 1:3 )
+   endif
+ enddo
+ write(iunit,'(A)') 
+ 
+
+ close(iunit) 
+end
+
+program prog_save_casino
+  call save_casino
+end