From ef3a52f54d5ef00d4abbc310b7e5beec1851b7e2 Mon Sep 17 00:00:00 2001
From: eginer <giner.emmanuel@gmail.com>
Date: Thu, 23 Feb 2023 16:38:40 +0100
Subject: [PATCH] added write_pt_charges.py

---
 src/hartree_fock/10.hf.bats    | 20 +++++----
 src/nuclei/write_pt_charges.py | 80 ++++++++++++++++++++++++++++++++++
 2 files changed, 91 insertions(+), 9 deletions(-)
 create mode 100644 src/nuclei/write_pt_charges.py

diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats
index 20b59603..df566032 100644
--- a/src/hartree_fock/10.hf.bats
+++ b/src/hartree_fock/10.hf.bats
@@ -34,23 +34,28 @@ cat > hcn_charges.xyz << EOF
 0.5   -2.0   0.0   0.0
 EOF
 
-rm -rf hcn.ezfio
-qp create_ezfio -b def2-svp hcn.xyz
+EZFIO=hcn_pt_charges
+rm -rf $EZFIO
+qp create_ezfio -b def2-svp hcn.xyz -o $EZFIO
 qp run scf
-mv hcn_charges.xyz hcn.ezfio_point_charges.xyz
-python write_pt_charges.py hcn.ezfio
+mv hcn_charges.xyz ${EZFIO}_point_charges.xyz
+python write_pt_charges.py ${EZFIO}
 qp set nuclei point_charges True
-qp run scf | tee hcn.ezfio.pt_charges.out
+qp run scf | tee ${EZFIO}.pt_charges.out
   energy="$(ezfio get hartree_fock energy)"
-rm -rf hcn.ezfio
 good=-92.76613324421798
   eq $energy $good $thresh
+rm -rf $EZFIO
 }
 
 @test "point charges" { 
  run_pt_charges
 }
 
+@test "HCN" { # 7.792500 8.51926s
+  run hcn.ezfio -92.88717500035233
+}
+
 @test "B-B" { # 3s
   run b2_stretched.ezfio -48.9950585434279
 }
@@ -124,9 +129,6 @@ good=-92.76613324421798
   run ch4.ezfio -40.19961807784367
 }
 
-@test "HCN" { # 7.792500 8.51926s
-  run hcn.ezfio -92.88717500035233
-}
 
 @test "N2" { # 8.648100 13.754s
   run n2.ezfio  -108.9834897852979
diff --git a/src/nuclei/write_pt_charges.py b/src/nuclei/write_pt_charges.py
new file mode 100644
index 00000000..6dbcd5b8
--- /dev/null
+++ b/src/nuclei/write_pt_charges.py
@@ -0,0 +1,80 @@
+#!/usr/bin/env python  
+import os
+import sys
+
+# First argument is the EZFIO file
+# It reads a file EZFIO_point_charges.xyz written in this way:
+# charge x y z (Angstrom)
+# for all charges
+
+
+def zip_in_ezfio(ezfio,tmp):
+  tmpzip=tmp+".gz"
+  cmdzip="gzip -c "+tmp+" > "+tmpzip
+  os.system(cmdzip)
+  os.system("rm "+tmp)
+  cmdmv="mv "+tmpzip+" "+EZFIO+"/nuclei/"+tmpzip
+  os.system(cmdmv)
+
+def mv_in_ezfio(ezfio,tmp):
+  cmdmv="mv "+tmp+" "+EZFIO+"/nuclei/"+tmp
+  os.system(cmdmv)
+
+
+# Getting the EZFIO
+EZFIO=sys.argv[1]
+EZFIO=EZFIO.replace("/", "")
+print(EZFIO)
+
+# Reading the point charges and convert the Angstrom geometry in Bohr for QP
+f = open(EZFIO+'_point_charges.xyz','r')
+lines = f.readlines()
+convert_angs_to_bohr=1.8897259885789233
+
+n_charges=0
+coord_x=[]
+coord_y=[]
+coord_z=[]
+charges=[]
+for line in lines:
+  data = line.split()
+  if(len(data)>0):
+   n_charges += 1
+   charges.append(str(data[0]))
+   coord_x.append(str(convert_angs_to_bohr*float(data[1])))
+   coord_y.append(str(convert_angs_to_bohr*float(data[2])))
+   coord_z.append(str(convert_angs_to_bohr*float(data[3])))
+
+# Write the file containing the number of charges and set in EZFIO folder
+tmp="n_pts_charge"
+fncharges = open(tmp,'w')
+fncharges.write(" "+str(n_charges)+'\n')
+fncharges.close()
+mv_in_ezfio(EZFIO,tmp)
+
+# Write the file containing the charges and set in EZFIO folder 
+tmp="pts_charge_z"
+fcharges = open(tmp,'w')
+fcharges.write(" 1\n")
+fcharges.write(" "+str(n_charges)+'\n')
+for i in range(n_charges):
+ fcharges.write(charges[i]+'\n')
+fcharges.close()
+zip_in_ezfio(EZFIO,tmp)
+
+# Write the file containing the charge coordinates and set in EZFIO folder
+tmp="pts_charge_coord"
+fcoord = open(tmp,'w')
+fcoord.write("  2\n")
+fcoord.write("                   "+str(n_charges)+'                    3\n')
+#fcoord.write(" "+'   3 '+str(n_charges)+' \n')
+for i in range(n_charges):
+ fcoord.write('   '+coord_x[i]+'\n')
+for i in range(n_charges):
+ fcoord.write('   '+coord_y[i]+'\n')
+for i in range(n_charges):
+ fcoord.write('   '+coord_z[i]+'\n')
+fcoord.close()
+zip_in_ezfio(EZFIO,tmp)
+
+