2016-02-03 00:37:03 +01:00
open Qputils
open Qptypes
2017-08-18 18:28:33 +02:00
open Core
2014-08-27 16:38:13 +02:00
2017-04-04 17:35:29 +02:00
type element =
| Element of Element . t
| Int_elem of ( Nucl_number . t * Element . t )
2015-11-17 22:25:26 +01:00
2016-01-25 15:44:15 +01:00
(* * Handle dummy atoms placed on bonds *)
2015-11-17 22:25:26 +01:00
let dummy_centers ~ threshold ~ molecule ~ nuclei =
let d =
Molecule . distance_matrix molecule
in
let n =
Array . length d
in
let nuclei =
Array . of_list nuclei
in
let rec aux accu = function
| ( - 1 , _ ) -> accu
| ( i , - 1 ) -> aux accu ( i - 1 , i - 1 )
| ( i , j ) when ( i > j ) ->
2019-01-14 15:20:51 +01:00
let new_accu =
2015-11-17 22:25:26 +01:00
let x , y =
2019-01-14 15:20:51 +01:00
Element . covalent_radius ( nuclei . ( i ) ) . Atom . element | > Positive_float . to_float ,
2015-11-17 22:25:26 +01:00
Element . covalent_radius ( nuclei . ( j ) ) . Atom . element | > Positive_float . to_float
2019-01-14 15:20:51 +01:00
in
2015-11-17 22:25:26 +01:00
let r =
( x + . y ) * . threshold
in
if d . ( i ) . ( j ) < r then
( i , x , j , y , d . ( i ) . ( j ) ) :: accu
else
accu
in aux new_accu ( i , j - 1 )
| ( i , j ) when ( i = j ) -> aux accu ( i , j - 1 )
| _ -> assert false
in
aux [] ( n - 1 , n - 1 )
| > List . map ~ f : ( fun ( i , x , j , y , r ) ->
2019-01-14 15:20:51 +01:00
let f =
2015-11-17 22:25:26 +01:00
x /. ( x + . y )
in
2019-01-14 15:20:51 +01:00
let u =
2015-11-17 22:25:26 +01:00
Point3d . of_tuple ~ units : Units . Bohr
2019-01-14 15:20:51 +01:00
( nuclei . ( i ) . Atom . coord . Point3d . x + .
( nuclei . ( j ) . Atom . coord . Point3d . x -. nuclei . ( i ) . Atom . coord . Point3d . x ) * . f ,
nuclei . ( i ) . Atom . coord . Point3d . y + .
( nuclei . ( j ) . Atom . coord . Point3d . y -. nuclei . ( i ) . Atom . coord . Point3d . y ) * . f ,
nuclei . ( i ) . Atom . coord . Point3d . z + .
( nuclei . ( j ) . Atom . coord . Point3d . z -. nuclei . ( i ) . Atom . coord . Point3d . z ) * . f )
2015-11-17 22:25:26 +01:00
in
Atom . { element = Element . X ; charge = Charge . of_int 0 ; coord = u }
)
2019-01-14 15:20:51 +01:00
2016-01-25 15:44:15 +01:00
(* * Returns the list of available basis sets *)
2015-11-25 23:46:09 +01:00
let list_basis () =
2019-01-14 15:20:51 +01:00
let basis_list =
" python2 " ^ Qpackage . root ^ " /external/emsl/EMSL_api.py list_basis "
| > Unix . open_process_in
2015-11-25 23:46:09 +01:00
| > In_channel . input_lines
2019-01-14 15:20:51 +01:00
| > List . map ~ f : ( fun x ->
2015-11-25 23:46:09 +01:00
match String . split x ~ on : '\'' with
| [] -> " "
2019-01-14 15:20:51 +01:00
| a :: []
2015-11-25 23:46:09 +01:00
| _ :: a :: _ -> String . strip a
)
2019-01-14 15:20:51 +01:00
in
2018-10-17 19:27:58 +02:00
List . sort basis_list ~ compare : String . ascending
2015-11-25 23:46:09 +01:00
2015-11-17 22:25:26 +01:00
2016-01-25 15:44:15 +01:00
(* * Run the program *)
2017-11-20 12:25:12 +01:00
let run ? o b au c d m p cart xyz_file =
2014-08-27 16:38:13 +02:00
(* Read molecule *)
let molecule =
2017-11-20 12:25:12 +01:00
if au then
( Molecule . of_file xyz_file ~ charge : ( Charge . of_int c )
~ multiplicity : ( Multiplicity . of_int m ) ~ units : Units . Bohr )
else
( Molecule . of_file xyz_file ~ charge : ( Charge . of_int c )
~ multiplicity : ( Multiplicity . of_int m ) )
2014-08-27 16:38:13 +02:00
in
2015-11-17 22:25:26 +01:00
let dummy =
dummy_centers ~ threshold : d ~ molecule ~ nuclei : molecule . Molecule . nuclei
in
let nuclei =
molecule . Molecule . nuclei @ dummy
in
2014-11-16 18:05:04 +01:00
2016-01-25 15:44:15 +01:00
(* * * * * * * * * *
Basis set
* * * * * * * * * * )
2015-11-24 17:21:38 +01:00
let basis_table =
Hashtbl . Poly . create ()
in
2014-11-16 18:05:04 +01:00
(* Open basis set channels *)
let basis_channel element =
2015-11-17 22:25:26 +01:00
let key =
2017-04-04 17:35:29 +02:00
match element with
| Element e -> Element . to_string e
| Int_elem ( i , e ) -> Printf . sprintf " %d,%s " ( Nucl_number . to_int i ) ( Element . to_string e )
2015-11-17 22:25:26 +01:00
in
2014-11-16 18:05:04 +01:00
match Hashtbl . find basis_table key with
2019-01-14 15:20:51 +01:00
| Some in_channel ->
2014-11-16 18:05:04 +01:00
in_channel
2018-10-17 19:27:58 +02:00
| None -> raise Caml . Not_found
2014-11-16 18:05:04 +01:00
in
2019-01-14 15:20:51 +01:00
2015-01-16 00:48:09 +01:00
let temp_filename =
Filename . temp_file " qp_create_ " " .basis "
in
2019-01-14 15:20:51 +01:00
let () =
2015-11-22 15:18:51 +01:00
Sys . remove temp_filename
in
2015-11-24 17:21:38 +01:00
let fetch_channel basis =
let command =
2019-01-14 15:20:51 +01:00
Qpackage . root ^ " /scripts/get_basis \" " ^ temp_filename
2015-11-24 17:21:38 +01:00
^ " . " ^ basis ^ " \" \" " ^ basis ^ " \" "
in
2019-01-14 15:20:51 +01:00
let long_basis =
2016-01-25 23:38:07 +01:00
Qpackage . root ^ " /data/basis/ " ^ basis
in
2019-01-14 15:20:51 +01:00
match
2016-01-25 23:38:07 +01:00
Sys . is_file basis ,
Sys . is_file long_basis
with
| ` Yes , _ -> In_channel . create basis
| ` No , ` Yes -> In_channel . create long_basis
2019-01-14 15:20:51 +01:00
| _ ->
2015-11-24 17:21:38 +01:00
begin
2019-01-14 15:20:51 +01:00
let filename =
2015-11-24 17:21:38 +01:00
Unix . open_process_in command
| > In_channel . input_all
| > String . strip
in
let new_channel =
In_channel . create filename
in
Unix . unlink filename ;
new_channel
end
in
2014-11-16 18:05:04 +01:00
let rec build_basis = function
| [] -> ()
2019-01-14 15:20:51 +01:00
| elem_and_basis_name :: rest ->
2014-11-16 18:05:04 +01:00
begin
match ( String . lsplit2 ~ on : ':' elem_and_basis_name ) with
| None -> (* Principal basis *)
begin
2019-01-14 15:20:51 +01:00
let basis =
elem_and_basis_name
2014-11-16 18:05:04 +01:00
in
2015-01-12 16:58:22 +01:00
let new_channel =
2015-11-24 17:21:38 +01:00
fetch_channel basis
2015-01-12 16:58:22 +01:00
in
List . iter nuclei ~ f : ( fun elem ->
2019-01-14 15:20:51 +01:00
let key =
2015-11-17 22:25:26 +01:00
Element . to_string elem . Atom . element
2014-11-16 18:05:04 +01:00
in
match Hashtbl . add basis_table ~ key : key ~ data : new_channel with
| ` Ok -> ()
| ` Duplicate -> ()
)
end
| Some ( key , basis ) -> (* Aux basis *)
begin
2019-01-14 15:20:51 +01:00
let elem =
2017-04-04 17:35:29 +02:00
try
Element ( Element . of_string key )
with Element . ElementError _ ->
2019-01-14 15:20:51 +01:00
let result =
2017-04-04 17:35:29 +02:00
match ( String . split ~ on : ',' key ) with
| i :: k :: [] -> ( Nucl_number . of_int @@ int_of_string i , Element . of_string k )
| _ -> failwith " Expected format is int,Element:basis "
in Int_elem result
2015-11-24 17:21:38 +01:00
and basis =
String . lowercase basis
2014-11-16 18:05:04 +01:00
in
2019-01-14 15:20:51 +01:00
let key =
2017-04-04 17:35:29 +02:00
match elem with
| Element e -> Element . to_string e
| Int_elem ( i , e ) -> Printf . sprintf " %d,%s " ( Nucl_number . to_int i ) ( Element . to_string e )
2014-11-16 18:05:04 +01:00
in
2015-11-24 17:21:38 +01:00
let new_channel =
fetch_channel basis
2015-01-12 16:58:22 +01:00
in
begin
match Hashtbl . add basis_table ~ key : key ~ data : new_channel with
| ` Ok -> ()
2017-04-04 17:35:29 +02:00
| ` Duplicate ->
2019-01-14 15:20:51 +01:00
let e =
2017-04-04 17:35:29 +02:00
match elem with
| Element e -> e
| Int_elem ( _ , e ) -> e
in
failwith ( " Duplicate definition of basis for " ^ ( Element . to_long_string e ) )
2015-01-12 16:58:22 +01:00
end
end
2014-11-16 18:05:04 +01:00
end ;
build_basis rest
in
2015-01-12 16:58:22 +01:00
String . split ~ on : '|' b
2019-01-14 15:20:51 +01:00
| > List . rev_map ~ f : String . strip
2014-11-16 18:05:04 +01:00
| > build_basis ;
2014-08-27 16:38:13 +02:00
2016-01-25 15:44:15 +01:00
2019-01-14 15:20:51 +01:00
(* * * * * * * * * * * * * * *
Pseudopotential
2016-01-25 15:44:15 +01:00
* * * * * * * * * * * * * * * )
let pseudo_table =
Hashtbl . Poly . create ()
in
(* Open pseudo channels *)
let pseudo_channel element =
let key =
Element . to_string element
in
2019-01-14 15:20:51 +01:00
Hashtbl . find pseudo_table key
2016-01-25 15:44:15 +01:00
in
let temp_filename =
Filename . temp_file " qp_create_ " " .pseudo "
in
let () =
Sys . remove temp_filename
in
let fetch_channel pseudo =
2016-01-25 23:38:07 +01:00
let long_pseudo =
Qpackage . root ^ " /data/pseudo/ " ^ pseudo
in
2019-01-14 15:20:51 +01:00
match
2016-01-25 23:38:07 +01:00
Sys . is_file pseudo ,
Sys . is_file long_pseudo
with
| ` Yes , _ -> In_channel . create pseudo
| ` No , ` Yes -> In_channel . create long_pseudo
2016-01-25 15:44:15 +01:00
| _ -> failwith ( " Pseudo file " ^ pseudo ^ " not found. " )
in
let rec build_pseudo = function
| [] -> ()
| elem_and_pseudo_name :: rest ->
begin
match ( String . lsplit2 ~ on : ':' elem_and_pseudo_name ) with
| None -> (* Principal pseudo *)
begin
let pseudo =
elem_and_pseudo_name
in
let new_channel =
fetch_channel pseudo
in
List . iter nuclei ~ f : ( fun elem ->
let key =
Element . to_string elem . Atom . element
in
match Hashtbl . add pseudo_table ~ key : key ~ data : new_channel with
| ` Ok -> ()
| ` Duplicate -> ()
)
end
| Some ( key , pseudo ) -> (* Aux pseudo *)
begin
let elem =
Element . of_string key
and pseudo =
String . lowercase pseudo
in
let key =
Element . to_string elem
in
let new_channel =
fetch_channel pseudo
in
begin
match Hashtbl . add pseudo_table ~ key : key ~ data : new_channel with
| ` Ok -> ()
| ` Duplicate -> failwith ( " Duplicate definition of pseudo for " ^ ( Element . to_long_string elem ) )
end
end
end ;
build_pseudo rest
in
2019-01-14 15:20:51 +01:00
let () =
2016-01-25 15:44:15 +01:00
match p with
| None -> ()
| Some p ->
String . split ~ on : '|' p
| > List . rev_map ~ f : String . strip
| > build_pseudo
in
2014-08-27 16:38:13 +02:00
(* Build EZFIO File name *)
let ezfio_file =
match o with
| Some x -> x
| None ->
begin
match String . rsplit2 ~ on : '.' xyz_file with
2019-01-14 15:20:51 +01:00
| Some ( x , " xyz " )
2016-09-08 22:37:05 +02:00
| Some ( x , " zmt " ) -> x ^ " .ezfio "
2014-08-27 16:38:13 +02:00
| _ -> xyz_file ^ " .ezfio "
end
in
if Sys . file_exists_exn ezfio_file then
failwith ( ezfio_file ^ " already exists " ) ;
2019-01-14 15:20:51 +01:00
let write_file () =
2016-02-03 00:37:03 +01:00
(* Create EZFIO *)
Ezfio . set_file ezfio_file ;
(* Write Pseudo *)
let pseudo =
List . map nuclei ~ f : ( fun x ->
match pseudo_channel x . Atom . element with
| Some channel -> Pseudo . read_element channel x . Atom . element
| None -> Pseudo . empty x . Atom . element
2019-01-14 15:20:51 +01:00
)
2016-02-03 00:37:03 +01:00
in
2016-01-25 15:44:15 +01:00
2019-01-14 15:20:51 +01:00
let molecule =
2016-02-03 00:37:03 +01:00
let n_elec_to_remove =
List . fold pseudo ~ init : 0 ~ f : ( fun accu x ->
accu + ( Positive_int . to_int x . Pseudo . n_elec ) )
2016-01-25 15:44:15 +01:00
in
2019-01-14 15:20:51 +01:00
{ Molecule . elec_alpha =
2016-02-03 00:37:03 +01:00
( Elec_alpha_number . to_int molecule . Molecule . elec_alpha )
- n_elec_to_remove / 2
| > Elec_alpha_number . of_int ;
2019-01-14 15:20:51 +01:00
Molecule . elec_beta =
2016-02-03 00:37:03 +01:00
( Elec_beta_number . to_int molecule . Molecule . elec_beta )
- ( n_elec_to_remove - n_elec_to_remove / 2 )
| > Elec_beta_number . of_int ;
2019-01-14 15:20:51 +01:00
Molecule . nuclei =
2016-02-03 00:37:03 +01:00
let charges =
List . map pseudo ~ f : ( fun x -> Positive_int . to_int x . Pseudo . n_elec
| > Float . of_int )
| > Array . of_list
in
List . mapi molecule . Molecule . nuclei ~ f : ( fun i x ->
{ x with Atom . charge = ( Charge . to_float x . Atom . charge ) -. charges . ( i )
| > Charge . of_float }
)
}
2016-01-25 15:44:15 +01:00
in
2016-02-03 00:37:03 +01:00
let nuclei =
molecule . Molecule . nuclei @ dummy
2016-01-25 15:44:15 +01:00
in
2019-01-14 15:20:51 +01:00
2016-02-03 00:37:03 +01:00
(* Write Electrons *)
Ezfio . set_electrons_elec_alpha_num ( Elec_alpha_number . to_int
molecule . Molecule . elec_alpha ) ;
Ezfio . set_electrons_elec_beta_num ( Elec_beta_number . to_int
molecule . Molecule . elec_beta ) ;
(* Write Nuclei *)
let labels =
List . map ~ f : ( fun x -> Element . to_string x . Atom . element ) nuclei
2019-01-14 15:20:51 +01:00
and charges =
List . map ~ f : ( fun x -> Atom . ( Charge . to_float x . charge ) ) nuclei
and coords =
2016-02-03 00:37:03 +01:00
( List . map ~ f : ( fun x -> x . Atom . coord . Point3d . x ) nuclei ) @
( List . map ~ f : ( fun x -> x . Atom . coord . Point3d . y ) nuclei ) @
( List . map ~ f : ( fun x -> x . Atom . coord . Point3d . z ) nuclei ) in
let nucl_num = ( List . length labels ) in
Ezfio . set_nuclei_nucl_num nucl_num ;
2019-01-14 15:20:51 +01:00
Ezfio . set_nuclei_nucl_label ( Ezfio . ezfio_array_of_list
2016-02-03 00:37:03 +01:00
~ rank : 1 ~ dim : [| nucl_num |] ~ data : labels ) ;
2019-01-14 15:20:51 +01:00
Ezfio . set_nuclei_nucl_charge ( Ezfio . ezfio_array_of_list
2016-02-03 00:37:03 +01:00
~ rank : 1 ~ dim : [| nucl_num |] ~ data : charges ) ;
2019-01-14 15:20:51 +01:00
Ezfio . set_nuclei_nucl_coord ( Ezfio . ezfio_array_of_list
2016-02-03 00:37:03 +01:00
~ rank : 2 ~ dim : [| nucl_num ; 3 |] ~ data : coords ) ;
(* Write pseudopotential *)
2019-01-14 15:20:51 +01:00
let () =
2016-02-03 00:37:03 +01:00
match p with
| None -> Ezfio . set_pseudo_do_pseudo false
| _ -> Ezfio . set_pseudo_do_pseudo true
in
2019-01-14 15:20:51 +01:00
let klocmax =
List . fold pseudo ~ init : 0 ~ f : ( fun accu x ->
let x =
2016-02-03 00:37:03 +01:00
List . length x . Pseudo . local
in
if ( x > accu ) then x
else accu
)
2019-01-14 15:20:51 +01:00
and lmax =
List . fold pseudo ~ init : 0 ~ f : ( fun accu x ->
let x =
2016-02-03 00:37:03 +01:00
List . fold x . Pseudo . non_local ~ init : 0 ~ f : ( fun accu ( x , _ ) ->
2019-01-14 15:20:51 +01:00
let x =
2017-05-31 23:47:18 +02:00
Positive_int . to_int x . Pseudo . GaussianPrimitive_non_local . proj
2016-02-03 00:37:03 +01:00
in
if ( x > accu ) then x
else accu
)
2016-01-25 15:44:15 +01:00
in
if ( x > accu ) then x
else accu
)
in
2016-02-03 00:37:03 +01:00
2019-01-14 15:20:51 +01:00
let kmax =
2016-02-03 00:37:03 +01:00
Array . init ( lmax + 1 ) ~ f : ( fun i ->
2019-01-14 15:20:51 +01:00
List . map pseudo ~ f : ( fun x ->
List . filter x . Pseudo . non_local ~ f : ( fun ( y , _ ) ->
2017-05-31 23:47:18 +02:00
( Positive_int . to_int y . Pseudo . GaussianPrimitive_non_local . proj ) = i )
2016-02-03 00:37:03 +01:00
| > List . length )
| > List . fold ~ init : 0 ~ f : ( fun accu x ->
if accu > x then accu else x )
)
| > Array . fold ~ init : 0 ~ f : ( fun accu i ->
if i > accu then i else accu )
2016-01-25 22:45:23 +01:00
in
2019-01-14 15:20:51 +01:00
let () =
2016-02-03 00:37:03 +01:00
Ezfio . set_pseudo_pseudo_klocmax klocmax ;
Ezfio . set_pseudo_pseudo_kmax kmax ;
Ezfio . set_pseudo_pseudo_lmax lmax ;
let tmp_array_v_k , tmp_array_dz_k , tmp_array_n_k =
Array . make_matrix ~ dimx : klocmax ~ dimy : nucl_num 0 . ,
Array . make_matrix ~ dimx : klocmax ~ dimy : nucl_num 0 . ,
Array . make_matrix ~ dimx : klocmax ~ dimy : nucl_num 0
2016-01-25 22:45:23 +01:00
in
2019-01-14 15:20:51 +01:00
List . iteri pseudo ~ f : ( fun j x ->
2016-02-03 00:37:03 +01:00
List . iteri x . Pseudo . local ~ f : ( fun i ( y , c ) ->
tmp_array_v_k . ( i ) . ( j ) <- AO_coef . to_float c ;
let y , z =
2019-01-14 15:20:51 +01:00
AO_expo . to_float y . Pseudo . GaussianPrimitive_local . expo ,
R_power . to_int y . Pseudo . GaussianPrimitive_local . r_power
2016-02-03 00:37:03 +01:00
in
tmp_array_dz_k . ( i ) . ( j ) <- y ;
tmp_array_n_k . ( i ) . ( j ) <- z ;
)
) ;
2019-01-14 15:20:51 +01:00
let concat_2d tmp_array =
let data =
2016-02-03 00:37:03 +01:00
Array . map tmp_array ~ f : Array . to_list
| > Array . to_list
| > List . concat
in
Ezfio . ezfio_array_of_list ~ rank : 2 ~ dim : [| nucl_num ; klocmax |] ~ data
2019-01-14 15:20:51 +01:00
in
concat_2d tmp_array_v_k
2016-02-03 00:37:03 +01:00
| > Ezfio . set_pseudo_pseudo_v_k ;
2019-01-14 15:20:51 +01:00
concat_2d tmp_array_dz_k
2016-02-03 00:37:03 +01:00
| > Ezfio . set_pseudo_pseudo_dz_k ;
2019-01-14 15:20:51 +01:00
concat_2d tmp_array_n_k
2016-02-03 00:37:03 +01:00
| > Ezfio . set_pseudo_pseudo_n_k ;
2019-01-14 15:20:51 +01:00
let tmp_array_v_kl , tmp_array_dz_kl , tmp_array_n_kl =
2016-02-03 00:37:03 +01:00
Array . init ( lmax + 1 ) ~ f : ( fun _ ->
( Array . make_matrix ~ dimx : kmax ~ dimy : nucl_num 0 . ) ) ,
Array . init ( lmax + 1 ) ~ f : ( fun _ ->
( Array . make_matrix ~ dimx : kmax ~ dimy : nucl_num 0 . ) ) ,
Array . init ( lmax + 1 ) ~ f : ( fun _ ->
( Array . make_matrix ~ dimx : kmax ~ dimy : nucl_num 0 ) )
in
2019-01-14 15:20:51 +01:00
List . iteri pseudo ~ f : ( fun j x ->
2016-02-03 00:37:03 +01:00
let last_idx =
Array . create ~ len : ( lmax + 1 ) 0
in
2019-01-14 15:20:51 +01:00
List . iter x . Pseudo . non_local ~ f : ( fun ( y , c ) ->
2016-02-03 00:37:03 +01:00
let k , y , z =
2017-05-31 23:47:18 +02:00
Positive_int . to_int y . Pseudo . GaussianPrimitive_non_local . proj ,
AO_expo . to_float y . Pseudo . GaussianPrimitive_non_local . expo ,
R_power . to_int y . Pseudo . GaussianPrimitive_non_local . r_power
2016-02-03 00:37:03 +01:00
in
2019-01-14 15:20:51 +01:00
let i =
last_idx . ( k )
2016-02-03 00:37:03 +01:00
in
tmp_array_v_kl . ( k ) . ( i ) . ( j ) <- AO_coef . to_float c ;
tmp_array_dz_kl . ( k ) . ( i ) . ( j ) <- y ;
tmp_array_n_kl . ( k ) . ( i ) . ( j ) <- z ;
last_idx . ( k ) <- i + 1 ;
)
2019-01-14 15:20:51 +01:00
) ;
let concat_3d tmp_array =
let data =
Array . map tmp_array ~ f : ( fun x ->
2016-02-03 00:37:03 +01:00
Array . map x ~ f : Array . to_list
| > Array . to_list
| > List . concat )
| > Array . to_list
| > List . concat
in
Ezfio . ezfio_array_of_list ~ rank : 3 ~ dim : [| nucl_num ; kmax ; lmax + 1 |] ~ data
in
concat_3d tmp_array_v_kl
| > Ezfio . set_pseudo_pseudo_v_kl ;
concat_3d tmp_array_dz_kl
| > Ezfio . set_pseudo_pseudo_dz_kl ;
concat_3d tmp_array_n_kl
| > Ezfio . set_pseudo_pseudo_n_kl ;
2016-01-25 22:45:23 +01:00
in
2019-01-14 15:20:51 +01:00
2016-01-25 15:44:15 +01:00
2016-02-03 00:37:03 +01:00
(* Write Basis set *)
let basis =
2014-11-15 11:01:30 +01:00
2016-02-03 00:37:03 +01:00
let nmax =
Nucl_number . get_max ()
in
let rec do_work ( accu : ( Atom . t * Nucl_number . t ) list ) ( n : int ) = function
| [] -> accu
| e :: tail ->
let new_accu =
( e , ( Nucl_number . of_int ~ max : nmax n ) ) :: accu
in
do_work new_accu ( n + 1 ) tail
in
let result = do_work [] 1 nuclei
| > List . rev
| > List . map ~ f : ( fun ( x , i ) ->
2019-01-14 15:20:51 +01:00
try
2016-02-03 00:37:03 +01:00
let e =
match x . Atom . element with
| Element . X -> Element . H
| e -> e
in
2017-04-04 17:35:29 +02:00
let key =
Int_elem ( i , x . Atom . element )
in
try
Basis . read_element ( basis_channel key ) i e
2018-10-17 19:27:58 +02:00
with Caml . Not_found ->
2017-04-04 17:35:29 +02:00
let key =
Element x . Atom . element
in
try
Basis . read_element ( basis_channel key ) i e
2018-10-17 19:27:58 +02:00
with Caml . Not_found ->
2017-04-04 17:35:29 +02:00
failwith ( Printf . sprintf " Basis not found for atom %d (%s) " ( Nucl_number . to_int i )
( Element . to_string x . Atom . element ) )
2016-02-03 00:37:03 +01:00
with
| End_of_file -> failwith
( " Element " ^ ( Element . to_string x . Atom . element ) ^ " not found in basis set. " )
2019-01-14 15:20:51 +01:00
)
2016-02-03 00:37:03 +01:00
| > List . concat
in
(* close all in_channels *)
result
2015-11-22 15:18:51 +01:00
in
2016-02-03 00:37:03 +01:00
let long_basis = Long_basis . of_basis basis in
let ao_num = List . length long_basis in
Ezfio . set_ao_basis_ao_num ao_num ;
Ezfio . set_ao_basis_ao_basis b ;
2019-01-14 15:20:51 +01:00
let ao_prim_num = List . map long_basis ~ f : ( fun ( _ , g , _ ) -> List . length g . Gto . lc )
2016-02-03 00:37:03 +01:00
and ao_nucl = List . map long_basis ~ f : ( fun ( _ , _ , n ) -> Nucl_number . to_int n )
2019-01-14 15:20:51 +01:00
and ao_power =
2016-02-03 00:37:03 +01:00
let l = List . map long_basis ~ f : ( fun ( x , _ , _ ) -> x ) in
( List . map l ~ f : ( fun t -> Positive_int . to_int Symmetry . Xyz . ( t . x ) ) ) @
( List . map l ~ f : ( fun t -> Positive_int . to_int Symmetry . Xyz . ( t . y ) ) ) @
2019-01-14 15:20:51 +01:00
( List . map l ~ f : ( fun t -> Positive_int . to_int Symmetry . Xyz . ( t . z ) ) )
2014-09-18 17:01:43 +02:00
in
2016-02-03 00:37:03 +01:00
let ao_prim_num_max = List . fold ~ init : 0 ~ f : ( fun s x ->
if x > s then x
else s ) ao_prim_num
2014-09-18 17:01:43 +02:00
in
2019-01-14 15:20:51 +01:00
let gtos =
2016-02-03 00:37:03 +01:00
List . map long_basis ~ f : ( fun ( _ , x , _ ) -> x )
2014-09-18 17:01:43 +02:00
in
2019-01-14 15:20:51 +01:00
let create_expo_coef ec =
let coefs =
2016-02-03 00:37:03 +01:00
begin match ec with
| ` Coefs -> List . map gtos ~ f : ( fun x ->
List . map x . Gto . lc ~ f : ( fun ( _ , coef ) -> AO_coef . to_float coef ) )
| ` Expos -> List . map gtos ~ f : ( fun x ->
List . map x . Gto . lc ~ f : ( fun ( prim , _ ) -> AO_expo . to_float
2017-05-31 23:47:18 +02:00
prim . GaussianPrimitive . expo ) )
2016-02-03 00:37:03 +01:00
end
in
let rec get_n n accu = function
2019-01-14 15:20:51 +01:00
| [] -> List . rev accu
2016-02-03 00:37:03 +01:00
| h :: tail ->
let y =
begin match List . nth h n with
| Some x -> x
| None -> 0 .
end
2019-01-14 15:20:51 +01:00
in
2016-02-03 00:37:03 +01:00
get_n n ( y :: accu ) tail
in
let rec build accu = function
| n when n = ao_prim_num_max -> accu
| n -> build ( accu @ ( get_n n [] coefs ) ) ( n + 1 )
in
build [] 0
in
let ao_coef = create_expo_coef ` Coefs
and ao_expo = create_expo_coef ` Expos
in
let () =
Ezfio . set_ao_basis_ao_prim_num ( Ezfio . ezfio_array_of_list
~ rank : 1 ~ dim : [| ao_num |] ~ data : ao_prim_num ) ;
Ezfio . set_ao_basis_ao_nucl ( Ezfio . ezfio_array_of_list
~ rank : 1 ~ dim : [| ao_num |] ~ data : ao_nucl ) ;
Ezfio . set_ao_basis_ao_power ( Ezfio . ezfio_array_of_list
~ rank : 2 ~ dim : [| ao_num ; 3 |] ~ data : ao_power ) ;
Ezfio . set_ao_basis_ao_coef ( Ezfio . ezfio_array_of_list
~ rank : 2 ~ dim : [| ao_num ; ao_prim_num_max |] ~ data : ao_coef ) ;
Ezfio . set_ao_basis_ao_expo ( Ezfio . ezfio_array_of_list
~ rank : 2 ~ dim : [| ao_num ; ao_prim_num_max |] ~ data : ao_expo ) ;
Ezfio . set_ao_basis_ao_cartesian ( cart ) ;
in
match Input . Ao_basis . read () with
| None -> failwith " Error in basis "
| Some x -> Input . Ao_basis . write x
2014-09-17 23:47:13 +02:00
in
2016-01-25 15:44:15 +01:00
let () =
2016-02-03 00:37:03 +01:00
try write_file () with
2019-01-14 15:20:51 +01:00
| ex ->
2016-02-03 00:37:03 +01:00
begin
begin
match Sys . is_directory ezfio_file with
| ` Yes -> rmdir ezfio_file
| _ -> ()
2019-01-14 15:20:51 +01:00
end ;
2016-02-03 00:37:03 +01:00
raise ex ;
end
2019-01-11 19:10:12 +01:00
in
ignore @@ Sys . command ( " qp_edit -c " ^ ezfio_file ) ;
print_endline ezfio_file
2019-01-14 15:20:51 +01:00
2015-11-17 22:25:26 +01:00
2014-08-27 16:38:13 +02:00
2015-11-25 23:46:09 +01:00
2019-01-11 19:10:12 +01:00
let () =
2015-11-25 23:46:09 +01:00
2019-01-14 19:52:41 +01:00
let open Command_line in
begin
" Creates an EZFIO directory from a standard xyz file or from a z-matrix file in Gaussian format. The basis set is defined as a single string if all the atoms are taken from the same basis set, otherwise specific elements can be defined as follows:
2015-11-25 23:46:09 +01:00
2019-01-14 19:52:41 +01:00
- b \ " cc-pcvdz | H:cc-pvdz | C:6-31g \"
- b \ " cc-pvtz | 1,H:sto-3g | 3,H:6-31g \"
2015-11-25 23:46:09 +01:00
2019-01-14 19:52:41 +01:00
If a file with the same name as the basis set exists , this file will be read . Otherwise , the basis set is obtained from the database .
" |> set_description_doc ;
set_header_doc ( Sys . argv . ( 0 ) ^ " - Quantum Package command " ) ;
2014-11-16 18:05:04 +01:00
2019-01-14 19:52:41 +01:00
[ { opt = Optional ; short = 'o' ; long = " output " ;
arg = With_arg " EZFIO_DIR " ;
doc = " Name of the created EZFIO directory. " } ;
2014-11-16 18:05:04 +01:00
2019-01-14 19:52:41 +01:00
{ opt = Mandatory ; short = 'b' ; long = " basis " ;
arg = With_arg " <string> " ;
doc = " Name of basis set. If <string>=show, the list of all basis sets is displayed. " } ;
2015-11-24 17:31:22 +01:00
2019-01-14 19:52:41 +01:00
{ opt = Optional ; short = 'a' ; long = " au " ;
arg = Without_arg ;
doc = " Input geometry is in atomic units. " } ;
2019-01-11 19:10:12 +01:00
2019-01-14 19:52:41 +01:00
{ opt = Optional ; short = 'c' ; long = " charge " ;
arg = With_arg " <int> " ;
doc = " Total charge of the molecule. Default is 0. " } ;
{ opt = Optional ; short = 'd' ; long = " dummy " ;
arg = With_arg " <float> " ;
doc = " Add dummy atoms. x * (covalent radii of the atoms). " } ;
{ opt = Optional ; short = 'm' ; long = " multiplicity " ;
arg = With_arg " <int> " ;
doc = " Spin multiplicity (2S+1) of the molecule. Default is 1. " } ;
{ opt = Optional ; short = 'p' ; long = " pseudo " ;
arg = With_arg " <string> " ;
doc = " Name of the pseudopotential. " } ;
{ opt = Optional ; short = 'x' ; long = " cartesian " ;
arg = Without_arg ;
doc = " Compute AOs in the Cartesian basis set (6d, 10f, ...). " } ;
anonymous " FILE " Mandatory " Input file in xyz format or z-matrix. " ;
]
| > set_specs
end ;
2019-01-14 15:20:51 +01:00
2019-01-11 19:10:12 +01:00
(* Handle options *)
let output =
2019-01-14 15:20:51 +01:00
Command_line . get " output "
2019-01-11 19:10:12 +01:00
in
2015-11-17 22:25:26 +01:00
2019-01-11 19:10:12 +01:00
let basis =
match Command_line . get " basis " with
2019-01-14 19:52:41 +01:00
| None -> assert false
2019-01-11 19:10:12 +01:00
| Some x -> x
in
2014-08-27 16:38:13 +02:00
2019-01-11 19:10:12 +01:00
let au =
Command_line . get_bool " au "
in
2015-11-17 22:25:26 +01:00
2019-01-11 19:10:12 +01:00
let charge =
match Command_line . get " charge " with
| None -> 0
| Some x -> int_of_string x
in
let dummy =
match Command_line . get " dummy " with
| None -> 0 .
| Some x -> float_of_string x
in
let multiplicity =
match Command_line . get " multiplicity " with
| None -> 1
| Some n -> int_of_string n
in
let pseudo =
2019-01-14 15:20:51 +01:00
Command_line . get " pseudo "
2019-01-11 19:10:12 +01:00
in
let cart =
2019-01-14 15:20:51 +01:00
Command_line . get_bool " cartesian "
2019-01-11 19:10:12 +01:00
in
2019-01-14 19:52:41 +01:00
if basis = " show " then
begin
list_basis ()
| > List . iter ~ f : print_endline ;
exit 0
end ;
2019-01-11 19:10:12 +01:00
let xyz_filename =
match Command_line . anon_args () with
| [ x ] -> x
| _ -> ( Command_line . help () ; failwith " input file is missing " )
in
2014-08-27 16:38:13 +02:00
2019-01-11 19:10:12 +01:00
run ? o : output basis au charge dummy multiplicity pseudo cart xyz_filename
2014-08-27 16:38:13 +02:00