diff --git a/.travis.yml b/.travis.yml index 5126a44c..dd28c132 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,17 +4,21 @@ # - sudo apt-get install gfortran liblapack-dev gcc # - sudo apt-get install graphviz +os: linux + dist: trusty sudo: false +compiler: gfortran + addons: apt: packages: - gfortran - gcc - - libblas-dev - - liblapack-dev +# - liblapack-dev +# - libblas-dev - graphviz cache: @@ -28,6 +32,7 @@ python: script: - ./configure ./config/travis.cfg - source ./quantum_package.rc ; qp_module.py install Full_CI Full_CI_ZMQ Hartree_Fock CAS_SD_ZMQ mrcepa0 All_singles + - source ./quantum_package.rc ; cd ~ ; install_lapack.sh ; cd $QP_ROOT - source ./quantum_package.rc ; ninja - source ./quantum_package.rc ; cd ocaml ; make ; cd - - source ./quantum_package.rc ; cd tests ; ./run_tests.sh -v diff --git a/README.md b/README.md index eebb67dd..a11c5713 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,10 @@ +## IMPORTANT + +If you have problems upgrading to the current version, consider re-installing everything from scratch including the OCaml compiler. +To do this, you will have to remove the `quantum_package` directory **and** the `$HOME/.opam` directory as well. + + + ![QP](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/qp.png) [![Build Status](https://travis-ci.org/LCPQ/quantum_package.svg?branch=master)](https://travis-ci.org/LCPQ/quantum_package) [![Gitter](https://badges.gitter.im/Join%20Chat.svg)](https://gitter.im/LCPQ/quantum_package?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) diff --git a/config/gfortran.cfg b/config/gfortran.cfg index c0aa875f..2256ccf4 100644 --- a/config/gfortran.cfg +++ b/config/gfortran.cfg @@ -11,7 +11,7 @@ # [COMMON] FC : gfortran -ffree-line-length-none -I . -LAPACK_LIB : -llapack -lblas +LAPACK_LIB : -lblas -llapack IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg index 91a12345..8a7ce9b8 100644 --- a/config/gfortran_debug.cfg +++ b/config/gfortran_debug.cfg @@ -11,7 +11,7 @@ # [COMMON] FC : gfortran -g -ffree-line-length-none -I . -LAPACK_LIB : -llapack -lblas +LAPACK_LIB : -lblas -llapack IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 --assert diff --git a/config/travis.cfg b/config/travis.cfg index c444db09..d2d75867 100644 --- a/config/travis.cfg +++ b/config/travis.cfg @@ -11,7 +11,7 @@ # [COMMON] FC : gfortran -ffree-line-length-none -I . -g -LAPACK_LIB : -llapack -lblas +LAPACK_LIB : -llapack -lrefblas -ltmglib IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 --assert diff --git a/install/scripts/install_lapack.sh b/install/scripts/install_lapack.sh new file mode 100755 index 00000000..25cbefc6 --- /dev/null +++ b/install/scripts/install_lapack.sh @@ -0,0 +1,10 @@ +#!/bin/bash -x + +git clone https://github.com/Reference-LAPACK/lapack-release.git +cd lapack-release +cp make.inc.example make.inc +make -j 8 +mv librefblas.a liblapack.a libtmglib.a $QP_ROOT/lib + + + diff --git a/install/scripts/install_ocaml.sh b/install/scripts/install_ocaml.sh index b82216d3..4d356bd3 100755 --- a/install/scripts/install_ocaml.sh +++ b/install/scripts/install_ocaml.sh @@ -5,8 +5,7 @@ QP_ROOT=$PWD cd - # Normal installation -PACKAGES="core cryptokit.1.10 ocamlfind sexplib ZMQ" -#ppx_sexp_conv +PACKAGES="core cryptokit.1.10 ocamlfind sexplib ZMQ ppx_sexp_conv ppx_deriving" # Needed for ZeroMQ export C_INCLUDE_PATH="${QP_ROOT}"/include:"${C_INCLUDE_PATH}" @@ -64,17 +63,18 @@ fi cd Downloads || exit 1 chmod +x ocaml.sh || exit 1 -echo N | ./ocaml.sh ${QP_ROOT}/bin/ || exit 1 +echo N | ./ocaml.sh ${QP_ROOT}/bin/ 4.04.2 || exit 1 ${QP_ROOT}/bin/opam config setup -a -q || exit 1 export LD_LIBRARY_PATH=${QP_ROOT}/lib:${LD_LIBRARY_PATH} export LIBRARY_PATH=${QP_ROOT}/lib:${LIBRARY_PATH} -export C_INCLUDE_PATH=${QP_ROOT}/lib:${C_INCLUDE_PATH} +export C_INCLUDE_PATH=${QP_ROOT}/lib:${QP_ROOT}/include:${C_INCLUDE_PATH} source ${HOME}/.opam/opam-init/init.sh > /dev/null 2> /dev/null || true NCPUs=$(cat /proc/cpuinfo | grep -i MHz | wc -l) +${QP_ROOT}/bin/opam install -j ${NCPUs} stdint.0.4.2 -y -q || exit 1 ${QP_ROOT}/bin/opam install -j ${NCPUs} ${PACKAGES} -y -q || exit 1 rm -f ../_build/ocaml.log diff --git a/ocaml/Address.ml b/ocaml/Address.ml index c819a463..a419806c 100644 --- a/ocaml/Address.ml +++ b/ocaml/Address.ml @@ -1,4 +1,3 @@ -open Core.Std module Tcp : sig type t @@ -8,7 +7,7 @@ module Tcp : sig end = struct type t = string let of_string x = - if not (String.is_prefix ~prefix:"tcp://" x) then + if not (String_ext.is_prefix ~prefix:"tcp://" x) then invalid_arg "Address Invalid" ; x @@ -26,7 +25,7 @@ module Ipc : sig end = struct type t = string let of_string x = - assert (String.is_prefix ~prefix:"ipc://" x); + assert (String_ext.is_prefix ~prefix:"ipc://" x); x let create name = Printf.sprintf "ipc://%s" name @@ -41,7 +40,7 @@ module Inproc : sig end = struct type t = string let of_string x = - assert (String.is_prefix ~prefix:"inproc://" x); + assert (String_ext.is_prefix ~prefix:"inproc://" x); x let create name = Printf.sprintf "inproc://%s" name diff --git a/ocaml/Atom.ml b/ocaml/Atom.ml index 72932b1f..5d385e89 100644 --- a/ocaml/Atom.ml +++ b/ocaml/Atom.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core exception AtomError of string @@ -6,7 +6,7 @@ type t = { element : Element.t ; charge : Charge.t ; coord : Point3d.t ; -} with sexp +} [@@deriving sexp] (** Read xyz coordinates of the atom *) let of_string ~units s = diff --git a/ocaml/Basis.ml b/ocaml/Basis.ml index 797d53f2..30af6577 100644 --- a/ocaml/Basis.ml +++ b/ocaml/Basis.ml @@ -1,7 +1,7 @@ -open Core.Std +open Sexplib.Std open Qptypes -type t = (Gto.t * Nucl_number.t) list with sexp +type t = (Gto.t * Nucl_number.t) list [@@deriving sexp] (** Read all the basis functions of an element *) let read in_channel at_number = @@ -16,7 +16,7 @@ let read in_channel at_number = (** Find an element in the basis set file *) let find in_channel element = - In_channel.seek in_channel 0L; + seek_in in_channel 0; let element_read = ref Element.X in while !element_read <> element do @@ -56,13 +56,13 @@ let to_string_general ~fmt ~atom_sep ?ele_array b = do_work ((Gto.to_string ~fmt g)::accu) n tail in do_work [new_nucleus 1] 1 b - |> String.concat ~sep:"\n" + |> String.concat "\n" let to_string_gamess ?ele_array = to_string_general ?ele_array ~fmt:Gto.Gamess ~atom_sep:"" let to_string_gaussian ?ele_array b = - String.concat ~sep:"\n" + String.concat "\n" [ to_string_general ?ele_array ~fmt:Gto.Gaussian ~atom_sep:"****" b ; "****" ] let to_string ?(fmt=Gto.Gamess) = diff --git a/ocaml/Bit.ml b/ocaml/Bit.ml index 28b8dac9..e5958ba6 100644 --- a/ocaml/Bit.ml +++ b/ocaml/Bit.ml @@ -1,4 +1,4 @@ -open Core.Std;; +open Core;; (* Type for bits @@ -11,7 +11,7 @@ Zero | One type t = | One | Zero -with sexp +[@@deriving sexp] let to_string = function | Zero -> "0" diff --git a/ocaml/Bit.mli b/ocaml/Bit.mli index 6dd5abdd..bc62d13c 100644 --- a/ocaml/Bit.mli +++ b/ocaml/Bit.mli @@ -1,4 +1,4 @@ -type t = One | Zero with sexp +type t = One | Zero [@@deriving sexp] (** String conversions for printing *) val to_string : t -> string diff --git a/ocaml/Bitlist.ml b/ocaml/Bitlist.ml index d7b9fc50..46d9ab9a 100644 --- a/ocaml/Bitlist.ml +++ b/ocaml/Bitlist.ml @@ -1,5 +1,5 @@ open Qptypes -open Core.Std +open Core (* Type for bits strings diff --git a/ocaml/Charge.ml b/ocaml/Charge.ml index 714a5690..8c4cc429 100644 --- a/ocaml/Charge.ml +++ b/ocaml/Charge.ml @@ -1,6 +1,6 @@ -open Core.Std +open Core -type t = float with sexp +type t = float [@@deriving sexp] let of_float x = x let of_int i = Float.of_int i diff --git a/ocaml/Charge.mli b/ocaml/Charge.mli index 07685531..9ed8e41a 100644 --- a/ocaml/Charge.mli +++ b/ocaml/Charge.mli @@ -1,4 +1,4 @@ -type t = float with sexp +type t = float [@@deriving sexp] (** Float conversion functions *) val to_float : t -> float diff --git a/ocaml/Determinant.ml b/ocaml/Determinant.ml index 3791e07e..a6d754b0 100644 --- a/ocaml/Determinant.ml +++ b/ocaml/Determinant.ml @@ -1,7 +1,7 @@ -open Core.Std;; -open Qptypes;; +open Qptypes +open Sexplib.Std -type t = int64 array with sexp +type t = int64 array [@@deriving sexp] let to_int64_array (x:t) = (x:int64 array) @@ -9,8 +9,8 @@ let to_int64_array (x:t) = (x:int64 array) let to_alpha_beta x = let x = to_int64_array x in let n_int = (Array.length x)/2 in - ( Array.init n_int ~f:(fun i -> x.(i)) , - Array.init n_int ~f:(fun i -> x.(i+n_int)) ) + ( Array.init n_int (fun i -> x.(i)) , + Array.init n_int (fun i -> x.(i+n_int)) ) let to_bitlist_couple x = @@ -28,12 +28,14 @@ let bitlist_to_string ~mo_tot_num x = let len = MO_number.to_int mo_tot_num in - List.map x ~f:(function - | Bit.Zero -> "-" - | Bit.One -> "+" - ) - |> String.concat - |> String.sub ~pos:0 ~len + let s = + List.map (function + | Bit.Zero -> "-" + | Bit.One -> "+" + ) x + |> String.concat "" + in + String.sub s 0 len @@ -77,6 +79,6 @@ let to_string ~mo_tot_num x = let (xa,xb) = to_bitlist_couple x in [ " " ; bitlist_to_string ~mo_tot_num xa ; "\n" ; " " ; bitlist_to_string ~mo_tot_num xb ] - |> String.concat + |> String.concat "" diff --git a/ocaml/Determinant.mli b/ocaml/Determinant.mli index da9fe02e..78a2f52f 100644 --- a/ocaml/Determinant.mli +++ b/ocaml/Determinant.mli @@ -5,7 +5,7 @@ * where each int64 is a list of 64 MOs. When the bit is set * to 1, the MO is occupied. *) -type t = int64 array with sexp +type t = int64 array [@@deriving sexp] (** Transform to an int64 array *) val to_int64_array : t -> int64 array diff --git a/ocaml/Element.ml b/ocaml/Element.ml index ebfd5e17..fd08b8da 100644 --- a/ocaml/Element.ml +++ b/ocaml/Element.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core open Qptypes exception ElementError of string @@ -11,7 +11,7 @@ type t = |K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr |Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe |Pt -with sexp +[@@deriving sexp] let of_string x = match (String.capitalize (String.lowercase x)) with diff --git a/ocaml/Element.mli b/ocaml/Element.mli index e78a57ad..fc6c679f 100644 --- a/ocaml/Element.mli +++ b/ocaml/Element.mli @@ -8,7 +8,7 @@ type t = |K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr |Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe |Pt -with sexp +[@@deriving sexp] (** String conversion functions *) val of_string : string -> t diff --git a/ocaml/Excitation.ml b/ocaml/Excitation.ml index d58c3093..58e18b11 100644 --- a/ocaml/Excitation.ml +++ b/ocaml/Excitation.ml @@ -1,14 +1,14 @@ -open Core.Std;; -open Qptypes;; +open Core +open Qptypes module Hole = struct - type t = MO_class.t with sexp + type t = MO_class.t [@@deriving sexp] let of_mo_class x = x let to_mo_class x = x end module Particle = struct - type t = MO_class.t with sexp + type t = MO_class.t [@@deriving sexp] let of_mo_class x = x let to_mo_class x = x end @@ -16,7 +16,7 @@ end type t = | Single of Hole.t*Particle.t | Double of Hole.t*Particle.t*Hole.t*Particle.t -with sexp;; +[@@deriving sexp] let create_single ~hole ~particle = MO_class.( @@ -29,7 +29,7 @@ let create_single ~hole ~particle = | ( _, Inactive _ ) -> failwith "Particles can not be in virtual MOs" | (h, p) -> Single ( (Hole.of_mo_class h), (Particle.of_mo_class p) ) ) -;; + let double_of_singles s1 s2 = let (h1,p1) = match s1 with @@ -40,14 +40,14 @@ let double_of_singles s1 s2 = | _ -> assert false in Double (h1,p1,h2,p2) -;; + let create_double ~hole1 ~particle1 ~hole2 ~particle2 = let s1 = create_single ~hole:hole1 ~particle:particle1 and s2 = create_single ~hole:hole2 ~particle:particle2 in double_of_singles s1 s2 -;; + let to_string = function | Single (h,p) -> @@ -68,5 +68,5 @@ let to_string = function (MO_class.to_string (Particle.to_mo_class p2)); "]"] |> String.concat ~sep:" " -;; + diff --git a/ocaml/Excitation.mli b/ocaml/Excitation.mli index 982cfd0e..ab9e083c 100644 --- a/ocaml/Excitation.mli +++ b/ocaml/Excitation.mli @@ -18,7 +18,7 @@ module Particle : type t = | Single of Hole.t * Particle.t | Double of Hole.t * Particle.t * Hole.t * Particle.t -with sexp +[@@deriving sexp] val create_single : hole:MO_class.t -> particle:MO_class.t -> t diff --git a/ocaml/GaussianPrimitive.ml b/ocaml/GaussianPrimitive.ml index cf3d7cdb..cdaabd87 100644 --- a/ocaml/GaussianPrimitive.ml +++ b/ocaml/GaussianPrimitive.ml @@ -1,10 +1,10 @@ open Qptypes -open Core.Std +open Core type t = { sym : Symmetry.t ; expo : AO_expo.t ; -} with sexp +} [@@deriving sexp] let to_string p = let { sym = s ; expo = e } = p in diff --git a/ocaml/Generic_input_of_rst.ml b/ocaml/Generic_input_of_rst.ml index 81388824..02aec5ee 100644 --- a/ocaml/Generic_input_of_rst.ml +++ b/ocaml/Generic_input_of_rst.ml @@ -1,4 +1,4 @@ -open Core.Std;; +open Core;; open Qptypes;; diff --git a/ocaml/Gto.ml b/ocaml/Gto.ml index 2c6efb88..ab265202 100644 --- a/ocaml/Gto.ml +++ b/ocaml/Gto.ml @@ -1,5 +1,5 @@ -open Core.Std open Qptypes +open Sexplib.Std exception GTO_Read_Failure of string exception End_Of_Basis @@ -11,11 +11,11 @@ type fmt = type t = { sym : Symmetry.t ; lc : ((GaussianPrimitive.t * AO_coef.t) list) -} with sexp +} [@@deriving sexp] let of_prim_coef_list pc = - let (p,c) = List.hd_exn pc in + let (p,c) = List.hd pc in let sym = p.GaussianPrimitive.sym in let rec check = function | [] -> `OK @@ -37,12 +37,12 @@ let of_prim_coef_list pc = let read_one in_channel = (* Fetch number of lines to read on first line *) let buffer = input_line in_channel in - if ( (String.strip buffer) = "" ) then + if ( (String_ext.strip buffer) = "" ) then raise End_Of_Basis; let sym_str = String.sub buffer 0 2 in let n_str = String.sub buffer 2 ((String.length buffer)-2) in - let sym = Symmetry.of_string (String.strip sym_str) in - let n = Int.of_string (String.strip n_str) in + let sym = Symmetry.of_string (String_ext.strip sym_str) in + let n = int_of_string (String_ext.strip n_str) in (* Read all the primitives *) let rec read_lines result = function | 0 -> result @@ -50,18 +50,19 @@ let read_one in_channel = begin let line_buffer = input_line in_channel in let buffer = line_buffer - |> String.split ~on:' ' - |> List.filter ~f:(fun x -> x <> "") + |> String_ext.split ~on:' ' + |> List.filter (fun x -> x <> "") in match buffer with | [ j ; expo ; coef ] -> begin - let coef = String.tr ~target:'D' ~replacement:'e' coef + let coef = + Str.global_replace (Str.regexp "D") "e" coef in let p = GaussianPrimitive.of_sym_expo sym - (AO_expo.of_float (Float.of_string expo) ) - and c = AO_coef.of_float (Float.of_string coef) in + (AO_expo.of_float (float_of_string expo) ) + and c = AO_coef.of_float (float_of_string coef) in read_lines ( (p,c)::result) (i-1) end | _ -> raise (GTO_Read_Failure line_buffer) @@ -89,7 +90,7 @@ let to_string_gamess { sym = sym ; lc = lc } = do_work (result::accu) (i+1) tail in (do_work [result] 1 lc) - |> String.concat ~sep:"\n" + |> String.concat "\n" (** Write the GTO in Gaussian format *) @@ -109,7 +110,7 @@ let to_string_gaussian { sym = sym ; lc = lc } = do_work (result::accu) (i+1) tail in (do_work [result] 1 lc) - |> String.concat ~sep:"\n" + |> String.concat "\n" (** Transform the gto to a string *) diff --git a/ocaml/Gto.mli b/ocaml/Gto.mli index 93b6c0f3..91534ebe 100644 --- a/ocaml/Gto.mli +++ b/ocaml/Gto.mli @@ -7,7 +7,7 @@ type fmt = type t = { sym : Symmetry.t ; lc : (GaussianPrimitive.t * Qptypes.AO_coef.t) list; - } with sexp + } [@@deriving sexp] (** Create from a list of GaussianPrimitive.t * Qptypes.AO_coef.t *) val of_prim_coef_list : diff --git a/ocaml/Input.ml b/ocaml/Input.ml index b1bb1060..65155f7c 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -1,6 +1,6 @@ open Qputils;; open Qptypes;; -open Core.Std;; +open Core;; include Input_ao_basis;; include Input_bitmasks;; diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index 8b0f72a2..348ddd7f 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -1,6 +1,6 @@ open Qptypes;; open Qputils;; -open Core.Std;; +open Core;; module Ao_basis : sig type t = @@ -13,7 +13,7 @@ module Ao_basis : sig ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; - } with sexp + } [@@deriving sexp] ;; val read : unit -> t option val to_string : t -> string @@ -32,7 +32,7 @@ end = struct ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; - } with sexp + } [@@deriving sexp] ;; let get_default = Qpackage.get_ezfio_default "ao_basis";; diff --git a/ocaml/Input_bi_integrals.ml b/ocaml/Input_bi_integrals.ml index c5fc4fe5..1caf38db 100644 --- a/ocaml/Input_bi_integrals.ml +++ b/ocaml/Input_bi_integrals.ml @@ -1,6 +1,6 @@ open Qptypes;; open Qputils;; -open Core.Std;; +open Core;; module Bielec_integrals : sig type t = @@ -11,7 +11,7 @@ module Bielec_integrals : sig threshold_ao : Threshold.t; threshold_mo : Threshold.t; direct : bool; - } with sexp + } [@@deriving sexp] ;; val read : unit -> t option val write : t -> unit @@ -27,7 +27,7 @@ end = struct threshold_ao : Threshold.t; threshold_mo : Threshold.t; direct : bool; - } with sexp + } [@@deriving sexp] ;; let get_default = Qpackage.get_ezfio_default "bielec_integrals";; diff --git a/ocaml/Input_bitmasks.ml b/ocaml/Input_bitmasks.ml index 0469c660..75a22a26 100644 --- a/ocaml/Input_bitmasks.ml +++ b/ocaml/Input_bitmasks.ml @@ -1,6 +1,6 @@ open Qptypes;; open Qputils;; -open Core.Std;; +open Core;; module Bitmasks : sig type t = @@ -10,7 +10,7 @@ module Bitmasks : sig generators : int64 array; n_mask_cas : Bitmask_number.t; cas : int64 array; - } with sexp + } [@@deriving sexp] ;; val read : unit -> t option val to_string : t -> string @@ -22,7 +22,7 @@ end = struct generators : int64 array; n_mask_cas : Bitmask_number.t; cas : int64 array; - } with sexp + } [@@deriving sexp] ;; let get_default = Qpackage.get_ezfio_default "bitmasks";; diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml index 6cc83745..48887ca0 100644 --- a/ocaml/Input_determinants_by_hand.ml +++ b/ocaml/Input_determinants_by_hand.ml @@ -1,6 +1,6 @@ open Qptypes;; open Qputils;; -open Core.Std;; +open Core;; module Determinants_by_hand : sig type t = @@ -11,7 +11,7 @@ module Determinants_by_hand : sig expected_s2 : Positive_float.t; psi_coef : Det_coef.t array; psi_det : Determinant.t array; - } with sexp + } [@@deriving sexp] val read : unit -> t val read_maybe : unit -> t option val write : t -> unit @@ -30,7 +30,7 @@ end = struct expected_s2 : Positive_float.t; psi_coef : Det_coef.t array; psi_det : Determinant.t array; - } with sexp + } [@@deriving sexp] ;; let get_default = Qpackage.get_ezfio_default "determinants";; diff --git a/ocaml/Input_electrons.ml b/ocaml/Input_electrons.ml index 24d0fe00..3779cfd2 100644 --- a/ocaml/Input_electrons.ml +++ b/ocaml/Input_electrons.ml @@ -1,12 +1,12 @@ open Qptypes;; open Qputils;; -open Core.Std;; +open Core;; module Electrons : sig type t = { elec_alpha_num : Elec_alpha_number.t; elec_beta_num : Elec_beta_number.t; - } with sexp + } [@@deriving sexp] ;; val read : unit -> t option val write : t -> unit @@ -18,7 +18,7 @@ end = struct type t = { elec_alpha_num : Elec_alpha_number.t; elec_beta_num : Elec_beta_number.t; - } with sexp + } [@@deriving sexp] ;; let get_default = Qpackage.get_ezfio_default "electrons";; diff --git a/ocaml/Input_mo_basis.ml b/ocaml/Input_mo_basis.ml index df47abfb..c51a7f4a 100644 --- a/ocaml/Input_mo_basis.ml +++ b/ocaml/Input_mo_basis.ml @@ -1,23 +1,29 @@ open Qptypes open Qputils -open Core.Std +open Core -type t_mo = - { mo_tot_num : MO_number.t ; - mo_label : MO_label.t; - mo_class : MO_class.t array; - mo_occ : MO_occ.t array; - mo_coef : (MO_coef.t array) array; - ao_md5 : MD5.t; - } with sexp module Mo_basis : sig - type t = t_mo + type t = + { mo_tot_num : MO_number.t ; + mo_label : MO_label.t; + mo_class : MO_class.t array; + mo_occ : MO_occ.t array; + mo_coef : (MO_coef.t array) array; + ao_md5 : MD5.t; + } [@@deriving sexp] val read : unit -> t option val to_string : t -> string val to_rst : t -> Rst_string.t end = struct - type t = t_mo + type t = + { mo_tot_num : MO_number.t ; + mo_label : MO_label.t; + mo_class : MO_class.t array; + mo_occ : MO_occ.t array; + mo_coef : (MO_coef.t array) array; + ao_md5 : MD5.t; + } [@@deriving sexp] let get_default = Qpackage.get_ezfio_default "mo_basis" let read_mo_label () = diff --git a/ocaml/Input_nuclei_by_hand.ml b/ocaml/Input_nuclei_by_hand.ml index f36b6b82..b4b8b0fe 100644 --- a/ocaml/Input_nuclei_by_hand.ml +++ b/ocaml/Input_nuclei_by_hand.ml @@ -1,6 +1,6 @@ open Qptypes;; open Qputils;; -open Core.Std;; +open Core;; module Nuclei_by_hand : sig type t = @@ -8,7 +8,7 @@ module Nuclei_by_hand : sig nucl_label : Element.t array; nucl_charge : Charge.t array; nucl_coord : Point3d.t array; - } with sexp + } [@@deriving sexp] ;; val read : unit -> t option val write : t -> unit @@ -22,7 +22,7 @@ end = struct nucl_label : Element.t array; nucl_charge : Charge.t array; nucl_coord : Point3d.t array; - } with sexp + } [@@deriving sexp] ;; let get_default = Qpackage.get_ezfio_default "nuclei";; diff --git a/ocaml/Io_ext.ml b/ocaml/Io_ext.ml new file mode 100644 index 00000000..489ed1ed --- /dev/null +++ b/ocaml/Io_ext.ml @@ -0,0 +1,24 @@ +let input_lines filename = + let in_channel = + open_in filename + in + let rec aux accu = + try + let newline = + input_line in_channel + in + aux (newline::accu) + with End_of_file -> accu + in + let result = + List.rev (aux []) + in + close_in in_channel; + result + + + +let read_all filename = + input_lines filename + |> String.concat "\n" + diff --git a/ocaml/Long_basis.ml b/ocaml/Long_basis.ml index 5c0ea6d6..06ea2ed5 100644 --- a/ocaml/Long_basis.ml +++ b/ocaml/Long_basis.ml @@ -1,7 +1,7 @@ -open Core.Std;; -open Qptypes;; +open Qptypes +open Sexplib.Std -type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list with sexp +type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp] let of_basis b = let rec do_work accu = function @@ -10,14 +10,14 @@ let of_basis b = begin let new_accu = Symmetry.Xyz.of_symmetry g.Gto.sym - |> List.rev_map ~f:(fun x-> (x,g,n)) + |> List.rev_map (fun x-> (x,g,n)) in do_work (new_accu@accu) tail end in do_work [] b |> List.rev -;; + let to_basis b = let rec do_work accu = function @@ -25,7 +25,7 @@ let to_basis b = | (s,g,n)::tail -> let first_sym = Symmetry.Xyz.of_symmetry g.Gto.sym - |> List.hd_exn + |> List.hd in let new_accu = if ( s = first_sym ) then @@ -36,19 +36,19 @@ let to_basis b = do_work new_accu tail in do_work [] b -;; + let to_string b = - let middle = List.map ~f:(fun (x,y,z) -> - "( "^((Int.to_string (Nucl_number.to_int z)))^", "^ + let middle = List.map (fun (x,y,z) -> + "( "^((string_of_int (Nucl_number.to_int z)))^", "^ (Symmetry.Xyz.to_string x)^", "^(Gto.to_string y) ^" )" ) b - |> String.concat ~sep:",\n" + |> String.concat ",\n" in "("^middle^")" -;; -include To_md5;; + +include To_md5 let to_md5 = to_md5 sexp_of_t -;; + diff --git a/ocaml/Long_basis.mli b/ocaml/Long_basis.mli index 7e69ecce..81f00539 100644 --- a/ocaml/Long_basis.mli +++ b/ocaml/Long_basis.mli @@ -5,7 +5,7 @@ open Qptypes;; * all the D orbitals are converted to xx, xy, xz, yy, yx * etc *) -type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list with sexp +type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list [@@deriving sexp] (** Transform a basis to a long basis *) val of_basis : diff --git a/ocaml/MO_class.ml b/ocaml/MO_class.ml index adf1a215..e85e9dda 100644 --- a/ocaml/MO_class.ml +++ b/ocaml/MO_class.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core open Qptypes type t = @@ -7,7 +7,7 @@ type t = | Active of MO_number.t list | Virtual of MO_number.t list | Deleted of MO_number.t list -with sexp +[@@deriving sexp] let to_string x = diff --git a/ocaml/MO_class.mli b/ocaml/MO_class.mli index 953e1afe..634fa95a 100644 --- a/ocaml/MO_class.mli +++ b/ocaml/MO_class.mli @@ -4,7 +4,7 @@ type t = | Active of Qptypes.MO_number.t list | Virtual of Qptypes.MO_number.t list | Deleted of Qptypes.MO_number.t list -with sexp +[@@deriving sexp] (** Create different excitation classes *) diff --git a/ocaml/MO_label.ml b/ocaml/MO_label.ml index 72bca28f..4e0b82f6 100644 --- a/ocaml/MO_label.ml +++ b/ocaml/MO_label.ml @@ -1,4 +1,4 @@ -open Core.Std;; +open Core;; type t = | Guess @@ -7,7 +7,7 @@ type t = | Localized | Orthonormalized | None -with sexp +[@@deriving sexp] ;; let to_string = function diff --git a/ocaml/MO_label.mli b/ocaml/MO_label.mli index d5061095..732bf1f2 100644 --- a/ocaml/MO_label.mli +++ b/ocaml/MO_label.mli @@ -5,7 +5,7 @@ type t = | Localized | Orthonormalized | None -with sexp +[@@deriving sexp] (** String representation *) val to_string : t -> string diff --git a/ocaml/Makefile b/ocaml/Makefile index 8519c973..3534c614 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -11,8 +11,8 @@ endif LIBS= PKGS= -OCAMLCFLAGS="-g -warn-error A" -OCAMLBUILD=ocamlbuild -j 0 -syntax camlp4o -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) +OCAMLCFLAGS="-g" +OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) MLLFILES=$(wildcard *.mll) MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml Input_auto_generated.ml qp_edit.ml MLIFILES=$(wildcard *.mli) git diff --git a/ocaml/Message.ml b/ocaml/Message.ml index faf5ed69..2ffd1da1 100644 --- a/ocaml/Message.ml +++ b/ocaml/Message.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core open Qptypes (** New job : Request to create a new multi-tasked job *) diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index eb4e0582..27ac4bd0 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -1,14 +1,14 @@ -open Core.Std ;; -open Qptypes ;; +open Qptypes +open Sexplib.Std -exception MultiplicityError of string;; -exception XYZError ;; +exception MultiplicityError of string +exception XYZError type t = { nuclei : Atom.t list ; elec_alpha : Elec_alpha_number.t ; elec_beta : Elec_beta_number.t ; -} with sexp +} [@@deriving sexp] let get_charge { nuclei ; elec_alpha ; elec_beta } = let result = @@ -19,7 +19,7 @@ let get_charge { nuclei ; elec_alpha ; elec_beta } = | a::rest -> (Charge.to_float a.Atom.charge) +. nucl_charge rest | [] -> 0. in - Charge.of_float (nucl_charge nuclei -. (Float.of_int result)) + Charge.of_float (nucl_charge nuclei -. (float_of_int result)) let get_multiplicity m = @@ -59,9 +59,10 @@ let name m = | a::rest -> begin let e = a.Atom.element in - match (List.Assoc.find accu e) with - | None -> build_list (List.Assoc.add accu e 1) rest - | Some i -> build_list (List.Assoc.add accu e (i+1)) rest + try + let i = List.assoc e accu in + build_list ( (e,i+1)::(List.remove_assoc e accu) ) rest + with Not_found -> build_list ( (e,1)::accu ) rest end | [] -> accu in @@ -83,7 +84,7 @@ let name m = let result = build_list [] nuclei |> build_name [c ; ", " ; mult] in - String.concat (result) + String.concat "" result let to_string_general ~f m = @@ -95,8 +96,8 @@ let to_string_general ~f m = let title = name m in - [ Int.to_string n ; title ] @ (List.map ~f nuclei) - |> String.concat ~sep:"\n" + [ string_of_int n ; title ] @ (List.map f nuclei) + |> String.concat "\n" let to_string = to_string_general ~f:(fun x -> Atom.to_string Units.Angstrom x) @@ -109,9 +110,9 @@ let of_xyz_string ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1)) ?(units=Units.Angstrom) s = - let l = String.split s ~on:'\n' - |> List.filter ~f:(fun x -> x <> "") - |> List.map ~f:(fun x -> Atom.of_string units x) + let l = String_ext.split s ~on:'\n' + |> List.filter (fun x -> x <> "") + |> List.map (fun x -> Atom.of_string units x) in let ne = ( get_charge { nuclei=l ; @@ -145,25 +146,28 @@ let of_xyz_file ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1)) ?(units=Units.Angstrom) filename = - let (x,buffer) = In_channel.read_all filename - |> String.lsplit2_exn ~on:'\n' + let lines = + match Io_ext.input_lines filename with + | natoms :: title :: rest -> + begin + try + if (int_of_string @@ String_ext.strip natoms) <= 0 then + raise XYZError + with + | _ -> raise XYZError + end; + String.concat "\n" rest + | _ -> raise XYZError in - let result = - try - (int_of_string @@ String.strip x) > 0 - with - | Failure "int_of_string" -> false - in - if not result then raise XYZError; - let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in - of_xyz_string ~charge ~multiplicity ~units buffer + of_xyz_string ~charge:charge ~multiplicity:multiplicity + ~units:units lines let of_zmt_file ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1)) ?(units=Units.Angstrom) filename = - In_channel.read_all filename + Io_ext.read_all filename |> Zmatrix.of_string |> Zmatrix.to_xyz_string |> of_xyz_string ~charge ~multiplicity ~units @@ -182,14 +186,14 @@ let of_file let distance_matrix molecule = let coord = molecule.nuclei - |> List.map ~f:(fun x -> x.Atom.coord) + |> List.map (fun x -> x.Atom.coord) |> Array.of_list in let n = Array.length coord in let result = - Array.make_matrix ~dimx:n ~dimy:n 0. + Array.make_matrix n n 0. in for i = 0 to (n-1) do @@ -203,6 +207,7 @@ let distance_matrix molecule = +open Core ;; include To_md5 let to_md5 = to_md5 sexp_of_t diff --git a/ocaml/Molecule.mli b/ocaml/Molecule.mli index f6201b18..fe80b48e 100644 --- a/ocaml/Molecule.mli +++ b/ocaml/Molecule.mli @@ -4,7 +4,7 @@ type t = { nuclei : Atom.t list; elec_alpha : Qptypes.Elec_alpha_number.t; elec_beta : Qptypes.Elec_beta_number.t; -} with sexp +} [@@deriving sexp] (** Returns the charge of the molecule *) val get_charge : t -> Charge.t diff --git a/ocaml/Multiplicity.ml b/ocaml/Multiplicity.ml index e56341c6..7570f483 100644 --- a/ocaml/Multiplicity.ml +++ b/ocaml/Multiplicity.ml @@ -1,7 +1,7 @@ -open Core.Std;; +open Core;; open Qptypes ;; -type t = Strictly_positive_int.t with sexp +type t = Strictly_positive_int.t [@@deriving sexp] let of_int = Strictly_positive_int.of_int ;; let to_int = Strictly_positive_int.to_int ;; diff --git a/ocaml/Multiplicity.mli b/ocaml/Multiplicity.mli index c6f8c6bf..0a28fb11 100644 --- a/ocaml/Multiplicity.mli +++ b/ocaml/Multiplicity.mli @@ -1,4 +1,4 @@ -type t = Qptypes.Strictly_positive_int.t with sexp +type t = Qptypes.Strictly_positive_int.t [@@deriving sexp] (** Conversion from int *) val of_int : int -> t diff --git a/ocaml/Point3d.ml b/ocaml/Point3d.ml index 5717ed39..7723ca82 100644 --- a/ocaml/Point3d.ml +++ b/ocaml/Point3d.ml @@ -1,11 +1,11 @@ -open Core.Std;; +open Core;; open Qptypes;; type t = { x : float ; y : float ; z : float ; -} with sexp +} [@@deriving sexp] let of_tuple ~units (x,y,z) = let f = match units with diff --git a/ocaml/Point3d.mli b/ocaml/Point3d.mli index 6d7428ec..476461a3 100644 --- a/ocaml/Point3d.mli +++ b/ocaml/Point3d.mli @@ -2,7 +2,7 @@ type t = { x : float; y : float; z : float; -} with sexp +} [@@deriving sexp] (** Create from a tuple of floats *) val of_tuple : units:Units.units -> float*float*float -> t diff --git a/ocaml/Primitive.mli b/ocaml/Primitive.mli index 77cb633a..f7d8809d 100644 --- a/ocaml/Primitive.mli +++ b/ocaml/Primitive.mli @@ -1,7 +1,7 @@ type t = { sym : Symmetry.t; expo : Qptypes.AO_expo.t; -} with sexp +} [@@deriving sexp] (** Conversion to string for printing *) val to_string : t -> string diff --git a/ocaml/Progress_bar.ml b/ocaml/Progress_bar.ml index 3473ac4b..ac5ae0d7 100644 --- a/ocaml/Progress_bar.ml +++ b/ocaml/Progress_bar.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core type t = { @@ -53,13 +53,13 @@ let display_tty bar = in let stop_time = let x = - Time.Span.to_float running_time + Time.Span.to_sec running_time in if (percent > 0.) then x *. 100. /. percent -. x - |> Time.Span.of_float + |> Time.Span.of_sec else - Time.Span.of_float 0. + Time.Span.of_sec 0. in Printf.eprintf "%s : [%s] %4.1f%% | %10s, ~%10s left\r%!" bar.title @@ -67,7 +67,7 @@ let display_tty bar = percent (Time.Span.to_string running_time) (stop_time |> Time.Span.to_string ); - { bar with dirty = false ; next = Time.add now (Time.Span.of_float 0.1) } + { bar with dirty = false ; next = Time.add now (Time.Span.of_sec 0.1) } let display_file bar = @@ -80,19 +80,19 @@ let display_file bar = in let stop_time = let x = - Time.Span.to_float running_time + Time.Span.to_sec running_time in if (percent > 0.) then x *. 100. /. percent -. x - |> Time.Span.of_float + |> Time.Span.of_sec else - Time.Span.of_float 0. + Time.Span.of_sec 0. in Printf.eprintf "%5.2f %% in %20s, ~%20s left\n%!" percent (Time.Span.to_string running_time) (Time.Span.to_string stop_time); - { bar with dirty = false ; next = Time.add (Time.now ()) (Time.Span.of_float 2.) } + { bar with dirty = false ; next = Time.add (Time.now ()) (Time.Span.of_sec 2.) } diff --git a/ocaml/Pseudo.ml b/ocaml/Pseudo.ml index 3791167d..8a59213c 100644 --- a/ocaml/Pseudo.ml +++ b/ocaml/Pseudo.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core open Qptypes @@ -7,7 +7,7 @@ module GaussianPrimitive_local : sig type t = { expo : AO_expo.t ; r_power : R_power.t ; - } with sexp + } [@@deriving sexp] val of_expo_r_power : AO_expo.t -> R_power.t -> t val to_string : t -> string @@ -17,7 +17,7 @@ end = struct type t = { expo : AO_expo.t ; r_power : R_power.t ; - } with sexp + } [@@deriving sexp] let of_expo_r_power dz n = { expo = dz ; r_power = n } @@ -35,7 +35,7 @@ module GaussianPrimitive_non_local : sig expo : AO_expo.t ; r_power : R_power.t ; proj : Positive_int.t - } with sexp + } [@@deriving sexp] val of_proj_expo_r_power : Positive_int.t -> AO_expo.t -> R_power.t -> t val to_string : t -> string @@ -46,7 +46,7 @@ end = struct expo : AO_expo.t ; r_power : R_power.t ; proj : Positive_int.t - } with sexp + } [@@deriving sexp] let of_proj_expo_r_power p dz n = { expo = dz ; r_power = n ; proj = p } @@ -66,7 +66,7 @@ type t = { n_elec : Positive_int.t ; local : (GaussianPrimitive_local.t * AO_coef.t ) list ; non_local : (GaussianPrimitive_non_local.t * AO_coef.t ) list -} with sexp +} [@@deriving sexp] let empty e = { element = e; diff --git a/ocaml/Qpackage.ml b/ocaml/Qpackage.ml index 8011b23b..e4759db5 100644 --- a/ocaml/Qpackage.ml +++ b/ocaml/Qpackage.ml @@ -1,4 +1,4 @@ -open Core.Std;; +open Core;; open Qptypes;; open Qputils;; diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index 36dc1102..4bde831c 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -1,4 +1,4 @@ -open Core.Std +open Sexplib (* let rec transpose = function @@ -14,12 +14,12 @@ let rec transpose = function let input_to_sexp s = let result = - String.split_lines s - |> List.filter ~f:(fun x-> - (String.strip x) <> "") - |> List.map ~f:(fun x-> - "("^(String.tr '=' ' ' x)^")") - |> String.concat + String_ext.split ~on:'\n' s + |> List.filter (fun x-> (String_ext.strip x) <> "") + |> List.map (fun x-> "("^ + (Str.global_replace (Str.regexp "=") " " x) + ^")") + |> String.concat "" in print_endline ("("^result^")"); "("^result^")" @@ -29,10 +29,10 @@ let rmdir dirname = let rec remove_one dir = Sys.chdir dir; Sys.readdir "." - |> Array.iter ~f:(fun x -> - match (Sys.is_directory x, Sys.is_file x) with - | (`Yes, _) -> remove_one x - | (_, `Yes) -> Sys.remove x + |> Array.iter (fun x -> + match (Sys.is_directory x, Sys.file_exists x) with + | (true, _) -> remove_one x + | (_, true) -> Sys.remove x | _ -> failwith ("Unable to remove file "^x^".") ); Sys.chdir ".."; diff --git a/ocaml/README_qp_edit.rst b/ocaml/README_qp_edit.rst index dfadea20..8af9cceb 100644 --- a/ocaml/README_qp_edit.rst +++ b/ocaml/README_qp_edit.rst @@ -19,7 +19,7 @@ of the block. r_y : Y_type.t ... last_r : bool - } with sexp + } [@@deriving sexp] ;; val read : unit -> t val write : t -> unit @@ -31,7 +31,7 @@ of the block. r_y : Y_type.t ... last_r : bool - } with sexp + } [@@deriving sexp] ;; let get_default = Qpackage.get_ezfio_default "new_keyword";; diff --git a/ocaml/Range.ml b/ocaml/Range.ml index 9607c6a3..91fcbcce 100644 --- a/ocaml/Range.ml +++ b/ocaml/Range.ml @@ -1,4 +1,4 @@ -open Core.Std;; +open Sexplib.Std (* A range is a string of the type: * @@ -12,14 +12,14 @@ open Core.Std;; *) -type t = int list with sexp +type t = int list [@@deriving sexp] let expand_range r = - match String.lsplit2 ~on:'-' r with + match String_ext.lsplit2 ~on:'-' r with | Some (s, f) -> begin - let start = Int.of_string s - and finish = Int.of_string f + let start = int_of_string s + and finish = int_of_string f in assert (start <= finish) ; let rec do_work = function @@ -31,9 +31,9 @@ let expand_range r = begin match r with | "" -> [] - | _ -> [Int.of_string r] + | _ -> [int_of_string r] end -;; + let of_string s = match s.[0] with @@ -43,36 +43,37 @@ let of_string s = assert (s.[0] = '[') ; assert (s.[(String.length s)-1] = ']') ; let s = String.sub s 1 ((String.length s) - 2) in - let l = String.split ~on:',' s in - let l = List.map ~f:expand_range l in - List.concat l |> List.dedup ~compare:Int.compare |> List.sort ~cmp:Int.compare -;; + let l = String_ext.split ~on:',' s in + let l = List.map expand_range l in + List.concat l + |> List.sort_uniq compare + let to_string l = let rec do_work buf symbol = function | [] -> buf | a::([] as t) -> - do_work (buf^symbol^(Int.to_string a)) "" t + do_work (buf^symbol^(string_of_int a)) "" t | a::(b::q as t) -> if (b-a = 1) then do_work buf "-" t else - do_work (buf^symbol^(Int.to_string a)^","^(Int.to_string b)) "" t + do_work (buf^symbol^(string_of_int a)^","^(string_of_int b)) "" t in let result = match l with | [] -> "[]" | h::t -> - do_work ("["^(Int.to_string h)) "" l in + do_work ("["^(string_of_int h)) "" l in (String.sub result 0 ((String.length result)))^"]" -;; + let test_module () = let s = "[72-107,36-53,126-131]" in let l = of_string s in print_string s ; print_newline () ; - List.iter ~f:(fun x -> Printf.printf "%d, " x) l ; print_newline () ; - to_string l |> print_string ; print_newline () ; -;; + List.iter (fun x -> Printf.printf "%d, " x) l ; print_newline () ; + to_string l |> print_string ; print_newline (); + diff --git a/ocaml/Range.mli b/ocaml/Range.mli index 2d56a0fa..e186ccf9 100644 --- a/ocaml/Range.mli +++ b/ocaml/Range.mli @@ -1,4 +1,4 @@ -type t = int list with sexp +type t = int list [@@deriving sexp] (** A range is a sorted list of ints in an interval. It is created using a string : diff --git a/ocaml/String_ext.ml b/ocaml/String_ext.ml new file mode 100644 index 00000000..ae8378bf --- /dev/null +++ b/ocaml/String_ext.ml @@ -0,0 +1,142 @@ +include String + +(** Split a string on a given character *) +let split ?(on=' ') str = + split_on_char on str +(* + let rec do_work ?(accu=[]) ?(left="") = function + | "" -> List.rev (left::accu) + | s -> + let new_s = + (length s) - 1 + |> sub s 1 + in + if (s.[0] = on) then + let new_accu = + left :: accu + in + do_work ~accu:new_accu new_s + else + let new_left = + concat "" [ left ; make 1 s.[0] ] + in + do_work ~accu ~left:new_left new_s + in + do_work str +*) + + +(** Strip blanks on the left of a string *) +let ltrim s = + let rec do_work s l = + match s.[0] with + | '\n' + | ' ' -> do_work (sub s 1 (l-1)) (l-1) + | _ -> s + in + let l = + length s + in + if (l > 0) then + do_work s l + else + s + +(** Strip blanks on the right of a string *) +let rtrim s = + let rec do_work s l = + let newl = + l-1 + in + match s.[newl] with + | '\n' + | ' ' -> do_work (sub s 0 (newl)) (newl) + | _ -> s + in + let l = + length s + in + if (l > 0) then + do_work s l + else + s + + +(** Strip blanks on the right and left of a string *) +let strip = String.trim + + +(** Split a string in two pieces when a character is found the 1st time from the left *) +let lsplit2_exn ?(on=' ') s = + let length = + String.length s + in + let rec do_work i = + if (i = length) then + begin + raise Not_found + end + else if (s.[i] = on) then + ( String.sub s 0 i, + String.sub s (i+1) (length-i-1) ) + else + do_work (i+1) + in + do_work 0 + + +(** Split a string in two pieces when a character is found the 1st time from the right *) +let rsplit2_exn ?(on=' ') s = + let length = + String.length s + in + let rec do_work i = + if (i = -1) then + begin + raise Not_found + end + else if (s.[i] = on) then + ( String.sub s 0 i, + String.sub s (i+1) (length-i-1) ) + else + do_work (i-1) + in + do_work length + + +let lsplit2 ?(on=' ') s = + try + Some (lsplit2_exn ~on s) + with + | Not_found -> None + + +let rsplit2 ?(on=' ') s = + try + Some (rsplit2_exn ~on s) + with + | Not_found -> None + + +let to_list s = + Array.init (String.length s) (fun i -> s.[i]) + |> Array.to_list + + +let fold ~init ~f s = + to_list s + |> List.fold_left f init + + +let is_prefix ~prefix s = + let len = + String.length prefix + in + if len > String.length s then + false + else + prefix = String.sub s 0 len + + +let of_char c = + String.make 1 c diff --git a/ocaml/Symmetry.ml b/ocaml/Symmetry.ml index 8647ae99..7b088b73 100644 --- a/ocaml/Symmetry.ml +++ b/ocaml/Symmetry.ml @@ -1,7 +1,7 @@ open Qptypes -open Core.Std +open Sexplib.Std -type t = S|P|D|F|G|H|I|J|K|L with sexp +type t = S|P|D|F|G|H|I|J|K|L [@@deriving sexp] let to_string = function | S -> "S" @@ -77,7 +77,7 @@ type st = t module Xyz = struct type t = { x: Positive_int.t ; y: Positive_int.t ; - z: Positive_int.t } with sexp + z: Positive_int.t } [@@deriving sexp] type state_type = Null | X | Y | Z (** Builds an XYZ triplet from a string. @@ -86,7 +86,7 @@ module Xyz = struct let flush state accu number = let n = if (number = "") then 1 - else (Int.of_string number) + else (int_of_string number) in match state with | X -> { x= Positive_int.(of_int ( (to_int accu.x) +n)); @@ -111,10 +111,9 @@ module Xyz = struct | 'Z'::rest | 'z'::rest -> let new_accu = flush state accu number in do_work Z new_accu "" rest - | c::rest -> do_work state accu (number^(String.of_char c)) rest + | c::rest -> do_work state accu (number^(String_ext.of_char c)) rest in - String.to_list_rev s - |> List.rev + String_ext.to_list s |> do_work Null { x=Positive_int.of_int 0 ; y=Positive_int.of_int 0 ; diff --git a/ocaml/Symmetry.mli b/ocaml/Symmetry.mli index baffeb2e..2ab63003 100644 --- a/ocaml/Symmetry.mli +++ b/ocaml/Symmetry.mli @@ -1,4 +1,4 @@ -type t = S | P | D | F | G | H | I | J | K | L with sexp +type t = S | P | D | F | G | H | I | J | K | L [@@deriving sexp] (** Creatio from strings *) val to_string : t -> string @@ -16,7 +16,7 @@ module Xyz : x : Qptypes.Positive_int.t; y : Qptypes.Positive_int.t; z : Qptypes.Positive_int.t; - } with sexp + } [@@deriving sexp] (** The string format contains the powers of x,y and z in a format like "x2z3" *) diff --git a/ocaml/TaskServer.ml b/ocaml/TaskServer.ml index 1ed403f7..45b438f2 100644 --- a/ocaml/TaskServer.ml +++ b/ocaml/TaskServer.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core open Qptypes @@ -63,7 +63,7 @@ let bind_socket ~socket_type ~socket ~port = ZMQ.Socket.bind socket @@ Printf.sprintf "tcp://*:%d" port; loop (-1) with - | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_float 1. ; loop (i-1) ) + | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_sec 1. ; loop (i-1) ) | other_exception -> raise other_exception in loop 60 diff --git a/ocaml/To_md5.ml b/ocaml/To_md5.ml index 33015d75..bc6608f9 100644 --- a/ocaml/To_md5.ml +++ b/ocaml/To_md5.ml @@ -1,5 +1,5 @@ -open Core.Std;; -open Qptypes;; +open Qptypes +open Sexplib let to_md5 sexp_of_t t = sexp_of_t t diff --git a/ocaml/_tags b/ocaml/_tags index 0935c0bb..aa36989b 100644 --- a/ocaml/_tags +++ b/ocaml/_tags @@ -1,3 +1,3 @@ -true: package(core,cryptokit,ZMQ,sexplib.syntax,str) +true: package(core,cryptokit,ZMQ,str,ppx_sexp_conv,ppx_deriving) true: thread false: profile diff --git a/ocaml/create_git_sha1.sh b/ocaml/create_git_sha1.sh index f1fb7fa6..35cbb7d5 100755 --- a/ocaml/create_git_sha1.sh +++ b/ocaml/create_git_sha1.sh @@ -4,9 +4,9 @@ SHA1=$(git log -1 | head -1 | cut -d ' ' -f 2) DATE=$(git log -1 | grep Date | cut -d ':' -f 2-) MESSAGE=$(git log -1 | tail -1 | sed 's/"/\\"/g') cat << EOF > Git.ml -open Core.Std -let sha1 = "$SHA1" |> String.strip -let date = "$DATE" |> String.strip -let message = "$MESSAGE" |> String.strip +open Core +let sha1 = "$SHA1" |> String_ext.strip +let date = "$DATE" |> String_ext.strip +let message = "$MESSAGE" |> String_ext.strip EOF diff --git a/ocaml/qp_basis_clean.ml b/ocaml/qp_basis_clean.ml index 32945ce3..2cc6d349 100644 --- a/ocaml/qp_basis_clean.ml +++ b/ocaml/qp_basis_clean.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core let filenames = let dir_name = Qpackage.root^"/data/basis/" diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index 787d4a0d..410590a2 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -1,6 +1,6 @@ open Qputils open Qptypes -open Core.Std +open Core let spec = let open Command.Spec in diff --git a/ocaml/qp_create_guess.ml b/ocaml/qp_create_guess.ml index bebfdad3..b841c350 100644 --- a/ocaml/qp_create_guess.ml +++ b/ocaml/qp_create_guess.ml @@ -1,6 +1,6 @@ open Qputils open Qptypes -open Core.Std +open Core let run ~multiplicity ezfio_file = if (not (Sys.file_exists_exn ezfio_file)) then diff --git a/ocaml/qp_find_pi_space.ml b/ocaml/qp_find_pi_space.ml new file mode 100644 index 00000000..0f5f7365 --- /dev/null +++ b/ocaml/qp_find_pi_space.ml @@ -0,0 +1,109 @@ +open Core +open Qputils +open Qptypes + +let run ?(sym="None") ezfio_filename = + Ezfio.set_file ezfio_filename ; + + let aos = + match Input.Ao_basis.read () with + | Some aos -> aos + | None -> failwith "Unable to read AOs" + and mos = + match Input.Mo_basis.read () with + | Some mos -> mos + | None -> failwith "Unable to read MOs" + in + let rec find_power symmetry accu = function + | -1 -> accu + | i -> let new_accu = + if aos.Input.Ao_basis.ao_power.(i) = symmetry then (i::accu) + else accu + in + find_power symmetry new_accu (i-1) + and n = + (AO_number.to_int aos.Input.Ao_basis.ao_num) + and m = + (MO_number.to_int mos.Input.Mo_basis.mo_tot_num) + and one = Positive_int.of_int 1 + and zero = Positive_int.of_int 0 + in + + (* Indices of the px, py and pz AOs *) + let x_indices = + find_power Symmetry.Xyz.{x=one ; y=zero ; z=zero} [] (n-1) + and y_indices = + find_power Symmetry.Xyz.{x=zero ; y=one ; z=zero} [] (n-1) + and z_indices = + find_power Symmetry.Xyz.{x=zero ; y=zero ; z=one } [] (n-1) + in + + (* Compute the relative weight of each MO on the px, py, pz spaces *) + let compute_weight mo_i list_aos = + let num = + List.fold_left ~f:(fun s i -> s +. (MO_coef.to_float @@ mos.Input.Mo_basis.mo_coef.(mo_i).(i)) ** 2.) ~init:0. list_aos + and denom = + Array.fold ~f:(fun s x -> s +. (MO_coef.to_float x) ** 2.) ~init:0. mos.Input.Mo_basis.mo_coef.(mo_i) + in + num /. denom + in + let result = + Array.init ~f:(fun mo_i -> + (mo_i+1, + compute_weight mo_i x_indices, + compute_weight mo_i y_indices, + compute_weight mo_i z_indices) ) m + |> Array.to_list + in + + let pi, sigma = + let rec aux test_xyz (accu_pi: int list) (accu_sigma: int list) = function + | [] -> (List.rev accu_pi, List.rev accu_sigma) + | (mo_i,x,y,z)::rest -> + if test_xyz (x,y,z) then + aux test_xyz (mo_i::accu_pi) accu_sigma rest + else + aux test_xyz accu_pi (mo_i::accu_sigma) rest + in + match sym with + | "x" | "X" -> aux (fun (x,y,z) -> (x>y && x>z)) [] [] result + | "y" | "Y" -> aux (fun (x,y,z) -> (y>x && y>z)) [] [] result + | "z" | "Z" -> aux (fun (x,y,z) -> (z>x && z>y)) [] [] result + | _ -> ([],[]) + in + + match sym with + | "x" | "X" | "y" | "Y" | "z" | "Z" -> + begin + Printf.printf "Pi: ["; + List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) pi; + Printf.printf "\b]\n\nSigma: ["; + List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) sigma; + Printf.printf "\b]\n" + end + | _ -> List.iter ~f:(fun (mo_i,x,y,z) -> Printf.printf "%d: (%f,%f,%f)\n" mo_i x y z) result + + + + + +let spec = + let open Command.Spec in + empty + +> flag "sym" (optional string) ~doc:"{x,y,z} Axis perpendicular to the plane" + +> anon ("ezfio_filename" %: string) + + +let command = + Command.basic + ~summary: "Quantum Package command" + ~readme:(fun () -> + "Find all the pi molecular orbitals to create a pi space. + ") + spec + (fun sym ezfio_filename () -> run ?sym ezfio_filename ) + + +let () = + Command.run command + diff --git a/ocaml/qp_overlap_of_wf.ml b/ocaml/qp_overlap_of_wf.ml index 816256fa..deea0556 100644 --- a/ocaml/qp_overlap_of_wf.ml +++ b/ocaml/qp_overlap_of_wf.ml @@ -1,6 +1,12 @@ +(** + * Computes the overlap where both Psi_0 and Psi_1 are truncated in the set + * of common determinants and normalized + *) + open Input_determinants_by_hand open Qptypes + let () = let ezfio, ezfio' = try @@ -42,16 +48,16 @@ let () = let overlap wf wf' = let result, norm, norm' = Hashtbl.fold (fun k c (accu,norm,norm') -> - let c' = - try Hashtbl.find wf' k - with Not_found -> 0. + let (c',c) = + try (Hashtbl.find wf' k, c) + with Not_found -> (0.,0.) in (accu +. c *. c' , norm +. c *. c , norm'+. c'*. c' ) ) wf (0.,0.,0.) in - result /. (norm *. norm') + result /. (sqrt (norm *. norm')) in let wf, wf' = @@ -62,5 +68,6 @@ let () = let o = overlap wf wf' in - print_float (abs_float o) + print_float (abs_float o) ; + print_newline () diff --git a/ocaml/qp_print.ml b/ocaml/qp_print.ml index ec584066..ea52bd7f 100644 --- a/ocaml/qp_print.ml +++ b/ocaml/qp_print.ml @@ -1,6 +1,6 @@ open Qputils;; open Qptypes;; -open Core.Std;; +open Core;; type input_action = | Basis diff --git a/ocaml/qp_print_basis.ml b/ocaml/qp_print_basis.ml index 52e10ce8..477ed366 100644 --- a/ocaml/qp_print_basis.ml +++ b/ocaml/qp_print_basis.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core open Qptypes let basis () = diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml index 5a656d2d..5be26539 100644 --- a/ocaml/qp_run.ml +++ b/ocaml/qp_run.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core open Qputils (* Environment variables : @@ -132,7 +132,7 @@ let run slave exe ezfio_file = Sys.remove qp_run_address_filename; let duration = Time.diff (Time.now()) time_start - |> Core.Span.to_string in + |> Time.Span.to_string in Printf.printf "Wall time : %s\n\n" duration; if (exit_code <> 0) then exit exit_code diff --git a/ocaml/qp_set_mo_class.ml b/ocaml/qp_set_mo_class.ml index 7451d87d..e82610d4 100644 --- a/ocaml/qp_set_mo_class.ml +++ b/ocaml/qp_set_mo_class.ml @@ -1,6 +1,6 @@ open Qputils open Qptypes -open Core.Std +open Core (* * Command-line arguments @@ -210,7 +210,7 @@ let get () = in let mo_tot_num = - MO_number.to_int data.Input_mo_basis.mo_tot_num + MO_number.to_int data.Input.Mo_basis.mo_tot_num in let n_int = @@ -244,7 +244,7 @@ let get () = | (MO_class.Deleted _) :: rest -> work ~del:(Printf.sprintf "%s,%d" del i) ~inact ~act ~virt ~core (i+1) rest in - work 1 (Array.to_list data.Input_mo_basis.mo_class) + work 1 (Array.to_list data.Input.Mo_basis.mo_class) diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index 06006181..600debf4 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -1,4 +1,10 @@ -open Core.Std;; +let global_replace x = + x + |> Str.global_replace (Str.regexp "Float.to_string") "string_of_float" + |> Str.global_replace (Str.regexp "Float.of_string") "float_of_string" + |> Str.global_replace (Str.regexp "Int.to_string") "string_of_int" + |> Str.global_replace (Str.regexp "Int.of_string") "int_of_string" + |> Str.global_replace (Str.regexp "String.\\(to\\|of\\)_string") "" let input_data = " * Positive_float : float @@ -118,8 +124,12 @@ let input_data = " * MD5 : string assert ((String.length x) = 32); - assert (String.fold x ~init:true ~f:(fun accu x -> - accu && (x < 'g'))); + assert ( + let a = + Array.init (String.length x) (fun i -> x.[i]) + in + Array.fold_left (fun accu x -> accu && (x < 'g')) true a + ); * Rst_string : string @@ -127,7 +137,7 @@ let input_data = " assert (x <> \"\") ; " -;; + let input_ezfio = " * MO_number : int @@ -156,25 +166,25 @@ let input_ezfio = " More than 10 billion of determinants " -;; + let untouched = " module MO_guess : sig - type t with sexp + type t [@@deriving sexp] val to_string : t -> string val of_string : string -> t end = struct type t = | Huckel | HCore - with sexp + [@@deriving sexp] let to_string = function | Huckel -> \"Huckel\" | HCore -> \"HCore\" let of_string s = - match (String.lowercase s) with + match (String.lowercase_ascii s) with | \"huckel\" -> Huckel | \"hcore\" -> HCore | _ -> raise (Invalid_argument (\"Wrong Guess type : \"^s)) @@ -182,7 +192,7 @@ end = struct end module Disk_access : sig - type t with sexp + type t [@@deriving sexp] val to_string : t -> string val of_string : string -> t end = struct @@ -190,14 +200,14 @@ end = struct | Read | Write | None - with sexp + [@@deriving sexp] let to_string = function | Read -> \"Read\" | Write -> \"Write\" | None -> \"None\" let of_string s = - match (String.lowercase s) with + match (String.lowercase_ascii s) with | \"read\" -> Read | \"write\" -> Write | \"none\" -> None @@ -206,62 +216,63 @@ end = struct end " -;; + let template = format_of_string " module %s : sig - type t with sexp + type t [@@deriving sexp] val to_%s : t -> %s val of_%s : %s %s -> t val to_string : t -> string end = struct - type t = %s with sexp + type t = %s [@@deriving sexp] let to_%s x = x let of_%s %s x = ( %s x ) let to_string x = %s.to_string x end " -;; + let parse_input input= - print_string "open Core.Std;;\nlet warning = print_string;;\n" ; + print_string "open Sexplib.Std\nlet warning = print_string\n" ; let rec parse result = function | [] -> result | ( "" , "" )::tail -> parse result tail | ( t , text )::tail -> let name,typ,params,params_val = - match String.split ~on:':' t with + match String_ext.split ~on:':' t with | [name;typ] -> (name,typ,"","") | name::typ::params::params_val -> (name,typ,params, - (String.concat params_val ~sep:":") ) + (String.concat ":" params_val) ) | _ -> assert false in - let typ = String.strip typ - and name = String.strip name in + let typ = String_ext.strip typ + and name = String_ext.strip name in let typ_cap = String.capitalize typ in let newstring = Printf.sprintf template name typ typ typ params_val typ typ - typ typ params ( String.strip text ) typ_cap + typ typ params ( String_ext.strip text ) typ_cap in List.rev (parse (newstring::result) tail ) in - String.split ~on:'*' input - |> List.map ~f:(String.lsplit2_exn ~on:'\n') + String_ext.split ~on:'*' input + |> List.map (String_ext.lsplit2_exn ~on:'\n') |> parse [] - |> String.concat + |> String.concat "" + |> global_replace |> print_string -;; + let ezfio_template = format_of_string " module %s : sig - type t with sexp + type t [@@deriving sexp] val to_%s : t -> %s val get_max : unit -> %s val of_%s : ?min:%s -> ?max:%s -> %s -> t val to_string : t -> string end = struct - type t = %s with sexp + type t = %s [@@deriving sexp] let to_string x = %s.to_string x let get_max () = if (Ezfio.has_%s ()) then @@ -287,24 +298,24 @@ end = struct end end " -;; + let parse_input_ezfio input= let parse s = match ( - String.split s ~on:'\n' - |> List.filter ~f:(fun x -> (String.strip x) <> "") + String_ext.split s ~on:'\n' + |> List.filter (fun x -> (String_ext.strip x) <> "") ) with | [] -> "" | a :: b :: c :: d :: [] -> begin - let (name,typ) = String.lsplit2_exn ~on:':' a + let (name,typ) = String_ext.lsplit2_exn ~on:':' a and ezfio_func = b - and (min, max) = String.lsplit2_exn ~on:':' c + and (min, max) = String_ext.lsplit2_exn ~on:':' c and msg = d in let (name, typ, ezfio_func, min, max, msg) = - match (List.map [ name ; typ ; ezfio_func ; min ; max ; msg ] ~f:String.strip) with + match List.map String_ext.strip [ name ; typ ; ezfio_func ; min ; max ; msg ] with | [ name ; typ ; ezfio_func ; min ; max ; msg ] -> (name, typ, ezfio_func, min, max, msg) | _ -> assert false in @@ -314,16 +325,17 @@ let parse_input_ezfio input= end | _ -> failwith "Error in input_ezfio" in - String.split ~on:'*' input - |> List.map ~f:parse - |> String.concat + String_ext.split ~on:'*' input + |> List.map parse + |> String.concat "" + |> global_replace |> print_string -;; + let () = parse_input input_data ; parse_input_ezfio input_ezfio; - print_endline untouched; + print_endline untouched diff --git a/ocaml/test_basis.ml b/ocaml/test_basis.ml index f58d30db..9007798c 100644 --- a/ocaml/test_basis.ml +++ b/ocaml/test_basis.ml @@ -1,4 +1,4 @@ -open Core.Std;; +open Core;; open Qputils;; open Qptypes;; diff --git a/ocaml/test_gto.ml b/ocaml/test_gto.ml index 423df62b..0ad6df56 100644 --- a/ocaml/test_gto.ml +++ b/ocaml/test_gto.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core open Qptypes let test_prim () = diff --git a/ocaml/test_message.ml b/ocaml/test_message.ml index 2f5592ec..6ccc381e 100644 --- a/ocaml/test_message.ml +++ b/ocaml/test_message.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core let () = Message.of_string "new_job ao_integrals tcp://127.0.0.1 inproc://ao_ints:12345" diff --git a/ocaml/test_molecule.ml b/ocaml/test_molecule.ml index 5d2935ed..e4b0f616 100644 --- a/ocaml/test_molecule.ml +++ b/ocaml/test_molecule.ml @@ -1,4 +1,4 @@ -open Core.Std ;; +open Core ;; open Qptypes ;; let test_molecule () = diff --git a/ocaml/test_pseudo.ml b/ocaml/test_pseudo.ml index 10c8aa0c..3b62680d 100644 --- a/ocaml/test_pseudo.ml +++ b/ocaml/test_pseudo.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core open Qputils open Qptypes diff --git a/ocaml/test_queuing_system.ml b/ocaml/test_queuing_system.ml index aa2fa280..728528b2 100644 --- a/ocaml/test_queuing_system.ml +++ b/ocaml/test_queuing_system.ml @@ -1,4 +1,4 @@ -open Core.Std +open Core let () = diff --git a/ocaml/test_symmetry.ml b/ocaml/test_symmetry.ml index bdc28e5c..e8323c14 100644 --- a/ocaml/test_symmetry.ml +++ b/ocaml/test_symmetry.ml @@ -1,15 +1,13 @@ -open Core.Std open Qputils open Qptypes open Symmetry let () = "SPDFGHIJKL" - |> String.to_list_rev - |> List.rev - |> List.map ~f:of_char - |> List.map ~f:Xyz.of_symmetry - |> List.iter ~f:(fun x -> List.iter x ~f:(fun y -> Xyz.to_string y |> print_endline) ; + |> String_ext.to_list + |> List.map of_char + |> List.map Xyz.of_symmetry + |> List.iter (fun x -> List.iter (fun y -> Xyz.to_string y |> print_endline) x ; print_newline ();) diff --git a/ocaml/test_task_server.ml b/ocaml/test_task_server.ml index e6a6106b..00573a9d 100644 --- a/ocaml/test_task_server.ml +++ b/ocaml/test_task_server.ml @@ -1,5 +1,3 @@ -open Core - let () = TaskServer.run 12345 diff --git a/plugins/CAS_SD_ZMQ/selection.irp.f b/plugins/CAS_SD_ZMQ/selection.irp.f index eb782cb3..b42116bd 100644 --- a/plugins/CAS_SD_ZMQ/selection.irp.f +++ b/plugins/CAS_SD_ZMQ/selection.irp.f @@ -1,28 +1,6 @@ use bitmasks -double precision function integral8(i,j,k,l) - implicit none - - integer, intent(in) :: i,j,k,l - double precision, external :: get_mo_bielec_integral - integer :: ii - ii = l-mo_integrals_cache_min - ii = ior(ii, k-mo_integrals_cache_min) - ii = ior(ii, j-mo_integrals_cache_min) - ii = ior(ii, i-mo_integrals_cache_min) - if (iand(ii, -64) /= 0) then - integral8 = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) - else - ii = l-mo_integrals_cache_min - ii = ior( ishft(ii,6), k-mo_integrals_cache_min) - ii = ior( ishft(ii,6), j-mo_integrals_cache_min) - ii = ior( ishft(ii,6), i-mo_integrals_cache_min) - integral8 = mo_integrals_cache(ii) - endif -end function - - BEGIN_PROVIDER [ integer(1), psi_phasemask, (N_int*bit_kind_size, 2, N_det)] use bitmasks implicit none @@ -287,7 +265,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2) integer :: i, j, h1, h2, p1, p2, sfix, hfix, pfix, hmob, pmob, puti double precision :: hij - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral integer, parameter :: turn3_2(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) integer, parameter :: turn2(2) = (/2,1/) @@ -300,7 +278,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) if(bannedOrb(puti)) cycle p1 = p(turn3_2(1,i), sp) p2 = p(turn3_2(2,i), sp) - hij = integral8(p1, p2, h1, h2) - integral8(p2, p1, h1, h2) + hij = mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2, p1, h1, h2) hij *= get_phase_bi(phasemask, sp, sp, h1, p1, h2, p2) vect(:, puti) += hij * coefs end do @@ -313,7 +291,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) puti = p(j, sp) if(bannedOrb(puti)) cycle pmob = p(turn2(j), sp) - hij = integral8(pfix, pmob, hfix, hmob) + hij = mo_bielec_integral(pfix, pmob, hfix, hmob) hij *= get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix) vect(:, puti) += hij * coefs end do @@ -325,7 +303,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) p2 = p(2,sfix) h1 = h(1,sfix) h2 = h(2,sfix) - hij = (integral8(p1,p2,h1,h2) - integral8(p2,p1,h1,h2)) + hij = (mo_bielec_integral(p1,p2,h1,h2) - mo_bielec_integral(p2,p1,h1,h2)) hij *= get_phase_bi(phasemask, sfix, sfix, h1, p1, h2, p2) vect(:, puti) += hij * coefs end if @@ -348,7 +326,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) logical :: ok, lbanned(mo_tot_num) integer(bit_kind) :: det(N_int, 2) double precision :: hij - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral lbanned = bannedOrb sh = 1 @@ -366,13 +344,13 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) do i=1,hole-1 if(lbanned(i)) cycle - hij = (integral8(p1, p2, i, hole) - integral8(p2, p1, i, hole)) + hij = (mo_bielec_integral(p1, p2, i, hole) - mo_bielec_integral(p2, p1, i, hole)) hij *= get_phase_bi(phasemask, sp, sp, i, p1, hole, p2) vect(:,i) += hij * coefs end do do i=hole+1,mo_tot_num if(lbanned(i)) cycle - hij = (integral8(p1, p2, hole, i) - integral8(p2, p1, hole, i)) + hij = (mo_bielec_integral(p1, p2, hole, i) - mo_bielec_integral(p2, p1, hole, i)) hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2) vect(:,i) += hij * coefs end do @@ -384,7 +362,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) p2 = p(1, sh) do i=1,mo_tot_num if(lbanned(i)) cycle - hij = integral8(p1, p2, i, hole) + hij = mo_bielec_integral(p1, p2, i, hole) hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2) vect(:,i) += hij * coefs end do @@ -787,7 +765,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral integer :: i, j, tip, ma, mi, puti, putj integer :: h1, h2, p1, p2, i1, i2 @@ -822,7 +800,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) h1 = h(1, ma) h2 = h(2, ma) - hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) + hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) if(ma == 1) then mat(:, putj, puti) += coefs * hij else @@ -841,7 +819,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) h1 = h(1,1) h2 = h(1,2) - hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) + hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) mat(:, puti, putj) += coefs * hij end do end do @@ -861,7 +839,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) i2 = turn2d(2, i, j) p1 = p(i1, ma) p2 = p(i2, ma) - hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) + hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) mat(:, puti, putj) += coefs * hij end do end do @@ -875,7 +853,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(banned(puti,putj,1)) cycle p2 = p(i, ma) - hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2) + hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2) mat(:, min(puti, putj), max(puti, putj)) += coefs * hij end do else ! tip == 4 @@ -886,7 +864,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) p2 = p(2, mi) h1 = h(1, mi) h2 = h(2, mi) - hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2) + hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2) mat(:, puti, putj) += coefs * hij end if end if @@ -905,7 +883,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) double precision, intent(in) :: coefs(N_states) double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) double precision :: hij, tmp_row(N_states, mo_tot_num), tmp_row2(N_states, mo_tot_num) - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral logical :: lbanned(mo_tot_num, 2), ok integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j, hfix, pfix, h1, h2, p1, p2, ib @@ -944,12 +922,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) tmp_row = 0d0 do putj=1, hfix-1 if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle - hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) + hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) tmp_row(1:N_states,putj) += hij * coefs(1:N_states) end do do putj=hfix+1, mo_tot_num if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle - hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) + hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) tmp_row(1:N_states,putj) += hij * coefs(1:N_states) end do @@ -969,13 +947,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) !p1 fixed putj = p1 if(.not. banned(putj,puti,bant)) then - hij = integral8(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix) + hij = mo_bielec_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix) tmp_row(:,puti) += hij * coefs end if putj = p2 if(.not. banned(putj,puti,bant)) then - hij = integral8(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix) + hij = mo_bielec_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix) tmp_row2(:,puti) += hij * coefs end if end do @@ -997,12 +975,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) tmp_row = 0d0 do putj=1,hfix-1 if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle - hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) + hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) tmp_row(:,putj) += hij * coefs end do do putj=hfix+1,mo_tot_num if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle - hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) + hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) tmp_row(:,putj) += hij * coefs end do @@ -1020,13 +998,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(lbanned(puti,ma)) cycle putj = p2 if(.not. banned(puti,putj,1)) then - hij = integral8(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1) + hij = mo_bielec_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1) tmp_row(:,puti) += hij * coefs end if putj = p1 if(.not. banned(puti,putj,1)) then - hij = integral8(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2) + hij = mo_bielec_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2) tmp_row2(:,puti) += hij * coefs end if end do @@ -1077,7 +1055,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) integer :: i, j, s, h1, h2, p1, p2, puti, putj double precision :: hij, phase - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral logical :: ok integer :: bant @@ -1096,7 +1074,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call apply_particles(mask, 1,p1,2,p2, det, ok, N_int) call i_h_j(gen, det, N_int, hij) else - hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) + hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) end if mat(:, p1, p2) += coefs(:) * hij @@ -1114,7 +1092,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int) call i_h_j(gen, det, N_int, hij) else - hij = (integral8(p1, p2, puti, putj) - integral8(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2) + hij = (mo_bielec_integral(p1, p2, puti, putj) - mo_bielec_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2) end if mat(:, puti, putj) += coefs(:) * hij end do diff --git a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f index 40d75fc4..7c30e175 100644 --- a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f +++ b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f @@ -29,7 +29,7 @@ subroutine diag_inactive_virt_and_update_mos label = "Canonical" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.false.) soft_touch mo_coef @@ -76,7 +76,7 @@ subroutine diag_inactive_virt_new_and_update_mos label = "Canonical" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.false.) soft_touch mo_coef diff --git a/plugins/FOBOCI/new_approach.irp.f b/plugins/FOBOCI/new_approach.irp.f index 8e2f2e53..6d359696 100644 --- a/plugins/FOBOCI/new_approach.irp.f +++ b/plugins/FOBOCI/new_approach.irp.f @@ -655,7 +655,7 @@ subroutine new_approach integer :: sign sign = -1 - call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,sign) + call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,sign,.true.) soft_touch mo_coef call save_mos diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f index 7d194a54..05f07b22 100644 --- a/plugins/FOBOCI/routines_foboci.irp.f +++ b/plugins/FOBOCI/routines_foboci.irp.f @@ -449,7 +449,7 @@ subroutine save_osoci_natural_mos label = "Natural" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.true.) !if(disk_access_ao_integrals == "None" .or. disk_access_ao_integrals == "Write" )then ! disk_access_ao_integrals = "Read" ! touch disk_access_ao_integrals @@ -584,7 +584,7 @@ subroutine set_osoci_natural_mos enddo label = "Natural" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.true.) soft_touch mo_coef deallocate(tmp,occ) diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f index acf19392..3e58224a 100644 --- a/plugins/Full_CI_ZMQ/selection.irp.f +++ b/plugins/Full_CI_ZMQ/selection.irp.f @@ -9,29 +9,6 @@ BEGIN_PROVIDER [ integer, fragment_count ] END_PROVIDER -double precision function integral8(i,j,k,l) - implicit none - - integer, intent(in) :: i,j,k,l - double precision, external :: get_mo_bielec_integral - integer :: ii - ii = l-mo_integrals_cache_min - ii = ior(ii, k-mo_integrals_cache_min) - ii = ior(ii, j-mo_integrals_cache_min) - ii = ior(ii, i-mo_integrals_cache_min) - if (iand(ii, -64) /= 0) then - integral8 = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) - else - ii = l-mo_integrals_cache_min - ii = ior( ishft(ii,6), k-mo_integrals_cache_min) - ii = ior( ishft(ii,6), j-mo_integrals_cache_min) - ii = ior( ishft(ii,6), i-mo_integrals_cache_min) - integral8 = mo_integrals_cache(ii) - endif -end function - - - subroutine assert(cond, msg) character(*), intent(in) :: msg logical, intent(in) :: cond @@ -135,7 +112,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2) integer :: i, j, h1, h2, p1, p2, sfix, hfix, pfix, hmob, pmob, puti double precision :: hij - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral integer, parameter :: turn3_2(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) integer, parameter :: turn2(2) = (/2,1/) @@ -148,7 +125,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) if(bannedOrb(puti)) cycle p1 = p(turn3_2(1,i), sp) p2 = p(turn3_2(2,i), sp) - hij = integral8(p1, p2, h1, h2) - integral8(p2, p1, h1, h2) + hij = mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2, p1, h1, h2) hij *= get_phase_bi(phasemask, sp, sp, h1, p1, h2, p2) vect(:, puti) += hij * coefs end do @@ -161,7 +138,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) puti = p(j, sp) if(bannedOrb(puti)) cycle pmob = p(turn2(j), sp) - hij = integral8(pfix, pmob, hfix, hmob) + hij = mo_bielec_integral(pfix, pmob, hfix, hmob) hij *= get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix) vect(:, puti) += hij * coefs end do @@ -173,7 +150,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) p2 = p(2,sfix) h1 = h(1,sfix) h2 = h(2,sfix) - hij = (integral8(p1,p2,h1,h2) - integral8(p2,p1,h1,h2)) + hij = (mo_bielec_integral(p1,p2,h1,h2) - mo_bielec_integral(p2,p1,h1,h2)) hij *= get_phase_bi(phasemask, sfix, sfix, h1, p1, h2, p2) vect(:, puti) += hij * coefs end if @@ -198,7 +175,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) logical, allocatable :: lbanned(:) integer(bit_kind) :: det(N_int, 2) double precision :: hij - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral allocate (lbanned(mo_tot_num)) lbanned = bannedOrb @@ -217,13 +194,13 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) do i=1,hole-1 if(lbanned(i)) cycle - hij = (integral8(p1, p2, i, hole) - integral8(p2, p1, i, hole)) + hij = (mo_bielec_integral(p1, p2, i, hole) - mo_bielec_integral(p2, p1, i, hole)) hij *= get_phase_bi(phasemask, sp, sp, i, p1, hole, p2) vect(1:N_states,i) += hij * coefs(1:N_states) end do do i=hole+1,mo_tot_num if(lbanned(i)) cycle - hij = (integral8(p1, p2, hole, i) - integral8(p2, p1, hole, i)) + hij = (mo_bielec_integral(p1, p2, hole, i) - mo_bielec_integral(p2, p1, hole, i)) hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2) vect(1:N_states,i) += hij * coefs(1:N_states) end do @@ -235,7 +212,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) p2 = p(1, sh) do i=1,mo_tot_num if(lbanned(i)) cycle - hij = integral8(p1, p2, i, hole) + hij = mo_bielec_integral(p1, p2, i, hole) hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2) vect(1:N_states,i) += hij * coefs(1:N_states) end do @@ -442,37 +419,82 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d fullinteresting(0) = 0 do ii=1,preinteresting(0) - i = preinteresting(ii) - mobMask(1,1) = iand(negMask(1,1), preinteresting_det(1,1,ii)) - mobMask(1,2) = iand(negMask(1,2), preinteresting_det(1,2,ii)) - nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) - do j=2,N_int - mobMask(j,1) = iand(negMask(j,1), preinteresting_det(j,1,ii)) - mobMask(j,2) = iand(negMask(j,2), preinteresting_det(j,2,ii)) - nt = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) - end do + select case (N_int) + case (1) + mobMask(1,1) = iand(negMask(1,1), preinteresting_det(1,1,ii)) + mobMask(1,2) = iand(negMask(1,2), preinteresting_det(1,2,ii)) + nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) + case (2) + mobMask(1:2,1) = iand(negMask(1:2,1), preinteresting_det(1:2,1,ii)) + mobMask(1:2,2) = iand(negMask(1:2,2), preinteresting_det(1:2,2,ii)) + nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) + & + popcnt(mobMask(2, 1)) + popcnt(mobMask(2, 2)) + case (3) + mobMask(1:3,1) = iand(negMask(1:3,1), preinteresting_det(1:3,1,ii)) + mobMask(1:3,2) = iand(negMask(1:3,2), preinteresting_det(1:3,2,ii)) + nt = 0 + do j=3,1,-1 + if (mobMask(j,1) /= 0_bit_kind) then + nt = nt+ popcnt(mobMask(j, 1)) + if (nt > 4) exit + endif + if (mobMask(j,2) /= 0_bit_kind) then + nt = nt+ popcnt(mobMask(j, 2)) + if (nt > 4) exit + endif + end do + case (4) + mobMask(1:4,1) = iand(negMask(1:4,1), preinteresting_det(1:4,1,ii)) + mobMask(1:4,2) = iand(negMask(1:4,2), preinteresting_det(1:4,2,ii)) + nt = 0 + do j=4,1,-1 + if (mobMask(j,1) /= 0_bit_kind) then + nt = nt+ popcnt(mobMask(j, 1)) + if (nt > 4) exit + endif + if (mobMask(j,2) /= 0_bit_kind) then + nt = nt+ popcnt(mobMask(j, 2)) + if (nt > 4) exit + endif + end do + case default + mobMask(1:N_int,1) = iand(negMask(1:N_int,1), preinteresting_det(1:N_int,1,ii)) + mobMask(1:N_int,2) = iand(negMask(1:N_int,2), preinteresting_det(1:N_int,2,ii)) + nt = 0 + do j=N_int,1,-1 + if (mobMask(j,1) /= 0_bit_kind) then + nt = nt+ popcnt(mobMask(j, 1)) + if (nt > 4) exit + endif + if (mobMask(j,2) /= 0_bit_kind) then + nt = nt+ popcnt(mobMask(j, 2)) + if (nt > 4) exit + endif + end do + end select - if(nt <= 4) then - interesting(0) += 1 - interesting(interesting(0)) = i + if(nt <= 4) then + i = preinteresting(ii) + interesting(0) += 1 + interesting(interesting(0)) = i minilist(1,1,interesting(0)) = preinteresting_det(1,1,ii) minilist(1,2,interesting(0)) = preinteresting_det(1,2,ii) - do j=2,N_int + do j=2,N_int minilist(j,1,interesting(0)) = preinteresting_det(j,1,ii) minilist(j,2,interesting(0)) = preinteresting_det(j,2,ii) - enddo - if(nt <= 2) then - fullinteresting(0) += 1 - fullinteresting(fullinteresting(0)) = i + enddo + if(nt <= 2) then + fullinteresting(0) += 1 + fullinteresting(fullinteresting(0)) = i fullminilist(1,1,fullinteresting(0)) = preinteresting_det(1,1,ii) fullminilist(1,2,fullinteresting(0)) = preinteresting_det(1,2,ii) - do j=2,N_int + do j=2,N_int fullminilist(j,1,fullinteresting(0)) = preinteresting_det(j,1,ii) fullminilist(j,2,fullinteresting(0)) = preinteresting_det(j,2,ii) - enddo - end if - end if - + enddo + end if + end if + end do do ii=1,prefullinteresting(0) @@ -481,12 +503,14 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i)) mobMask(1,2) = iand(negMask(1,2), psi_det_sorted(1,2,i)) nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) - do j=2,N_int + if (nt > 2) cycle + do j=N_int,2,-1 mobMask(j,1) = iand(negMask(j,1), psi_det_sorted(j,1,i)) mobMask(j,2) = iand(negMask(j,2), psi_det_sorted(j,2,i)) nt = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) + if (nt > 2) exit end do - + if(nt <= 2) then fullinteresting(0) += 1 fullinteresting(fullinteresting(0)) = i @@ -722,7 +746,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral integer :: i, j, tip, ma, mi, puti, putj integer :: h1, h2, p1, p2, i1, i2 @@ -757,7 +781,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) h1 = h(1, ma) h2 = h(2, ma) - hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) + hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) if(ma == 1) then mat(:, putj, puti) += coefs * hij else @@ -776,7 +800,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(banned(puti,putj,bant)) cycle p1 = p(turn2(i), 1) - hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) + hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) mat(:, puti, putj) += coefs * hij end do end do @@ -796,7 +820,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) i2 = turn2d(2, i, j) p1 = p(i1, ma) p2 = p(i2, ma) - hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) + hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) mat(:, puti, putj) += coefs * hij end do end do @@ -810,7 +834,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(banned(puti,putj,1)) cycle p2 = p(i, ma) - hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2) + hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2) mat(:, min(puti, putj), max(puti, putj)) += coefs * hij end do else ! tip == 4 @@ -821,7 +845,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) p2 = p(2, mi) h1 = h(1, mi) h2 = h(2, mi) - hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2) + hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2) mat(:, puti, putj) += coefs * hij end if end if @@ -841,7 +865,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp double precision :: hij, tmp_row(N_states, mo_tot_num), tmp_row2(N_states, mo_tot_num) - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral logical :: ok logical, allocatable :: lbanned(:,:) @@ -881,12 +905,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) tmp_row = 0d0 do putj=1, hfix-1 if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle - hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) + hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) tmp_row(1:N_states,putj) += hij * coefs(1:N_states) end do do putj=hfix+1, mo_tot_num if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle - hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) + hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) tmp_row(1:N_states,putj) += hij * coefs(1:N_states) end do @@ -906,13 +930,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) !p1 fixed putj = p1 if(.not. banned(putj,puti,bant)) then - hij = integral8(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix) + hij = mo_bielec_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix) tmp_row(:,puti) += hij * coefs end if putj = p2 if(.not. banned(putj,puti,bant)) then - hij = integral8(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix) + hij = mo_bielec_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix) tmp_row2(:,puti) += hij * coefs end if end do @@ -934,12 +958,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) tmp_row = 0d0 do putj=1,hfix-1 if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle - hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) + hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) tmp_row(:,putj) += hij * coefs end do do putj=hfix+1,mo_tot_num if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle - hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) + hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) tmp_row(:,putj) += hij * coefs end do @@ -957,13 +981,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(lbanned(puti,ma)) cycle putj = p2 if(.not. banned(puti,putj,1)) then - hij = integral8(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1) + hij = mo_bielec_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1) tmp_row(:,puti) += hij * coefs end if putj = p1 if(.not. banned(puti,putj,1)) then - hij = integral8(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2) + hij = mo_bielec_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2) tmp_row2(:,puti) += hij * coefs end if end do @@ -1015,7 +1039,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) integer :: i, j, s, h1, h2, p1, p2, puti, putj double precision :: hij, phase - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral logical :: ok integer :: bant @@ -1035,7 +1059,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call i_h_j(gen, det, N_int, hij) else phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) - hij = integral8(p1, p2, h1, h2) * phase + hij = mo_bielec_integral(p1, p2, h1, h2) * phase end if mat(:, p1, p2) += coefs(:) * hij end do @@ -1052,7 +1076,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int) call i_h_j(gen, det, N_int, hij) else - hij = (integral8(p1, p2, puti, putj) - integral8(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2) + hij = (mo_bielec_integral(p1, p2, puti, putj) - mo_bielec_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2) end if mat(:, puti, putj) += coefs(:) * hij end do diff --git a/plugins/Hartree_Fock/DIIS.irp.f b/plugins/Hartree_Fock/DIIS.irp.f index 614a9173..6ebb5879 100644 --- a/plugins/Hartree_Fock/DIIS.irp.f +++ b/plugins/Hartree_Fock/DIIS.irp.f @@ -94,7 +94,7 @@ END_PROVIDER do i=1,AO_num do j=1,AO_num - Xt(i,j) = X_Matrix_AO(j,i) + Xt(i,j) = S_half_inv(j,i) enddo enddo @@ -103,7 +103,7 @@ END_PROVIDER call dgemm('N','N',AO_num,AO_num,AO_num, & 1.d0, & Fock_matrix_AO,size(Fock_matrix_AO,1), & - X_Matrix_AO,size(X_Matrix_AO,1), & + S_half_inv,size(S_half_inv,1), & 0.d0, & eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1)) @@ -130,67 +130,10 @@ END_PROVIDER call dgemm('N','N',AO_num,AO_num,AO_num, & 1.d0, & - X_matrix_AO,size(X_matrix_AO,1), & + S_half_inv,size(S_half_inv,1), & scratch,size(scratch,1), & 0.d0, & eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1)) END_PROVIDER -BEGIN_PROVIDER [ double precision, X_matrix_AO, (AO_num,AO_num) ] - - BEGIN_DOC -! Matrix X = S^{-1/2} obtained by SVD - END_DOC - - implicit none - - integer :: num_linear_dependencies - integer :: LDA, LDC - double precision, allocatable :: U(:,:),Vt(:,:), D(:) - integer :: info, i, j, k - - LDA = size(AO_overlap,1) - LDC = size(X_matrix_AO,1) - - allocate( & - U(LDC,AO_num), & - Vt(LDA,AO_num), & - D(AO_num)) - - call svd( & - AO_overlap,LDA, & - U,LDC, & - D, & - Vt,LDA, & - AO_num,AO_num) - - num_linear_dependencies = 0 - do i=1,AO_num - print*,D(i) - if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then - D(i) = 0.d0 - num_linear_dependencies += 1 - else - ASSERT (D(i) > 0.d0) - D(i) = 1.d0/sqrt(D(i)) - endif - do j=1,AO_num - X_matrix_AO(j,i) = 0.d0 - enddo - enddo - write(*,*) 'linear dependencies',num_linear_dependencies -! stop - - do k=1,AO_num - if(D(k) /= 0.d0) then - do j=1,AO_num - do i=1,AO_num - X_matrix_AO(i,j) = X_matrix_AO(i,j) + U(i,k)*D(k)*Vt(k,j) - enddo - enddo - endif - enddo - - -END_PROVIDER diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg index a4b646e1..7f67d961 100644 --- a/plugins/Hartree_Fock/EZFIO.cfg +++ b/plugins/Hartree_Fock/EZFIO.cfg @@ -1,9 +1,3 @@ -[threshold_overlap_ao_eigenvalues] -type: Threshold -doc: Threshold on the magnitude of the smallest eigenvalues of the overlap matrix in the AO basis -interface: ezfio,provider,ocaml -default: 1.e-6 - [max_dim_diis] type: integer doc: Maximum size of the DIIS extrapolation procedure @@ -32,7 +26,7 @@ default: 500 type: Positive_float doc: Energy shift on the virtual MOs to improve SCF convergence interface: ezfio,provider,ocaml -default: 0.0 +default: 0.1 [scf_algorithm] type: character*(32) diff --git a/plugins/Hartree_Fock/Fock_matrix.irp.f b/plugins/Hartree_Fock/Fock_matrix.irp.f index 087933c8..0764c83f 100644 --- a/plugins/Hartree_Fock/Fock_matrix.irp.f +++ b/plugins/Hartree_Fock/Fock_matrix.irp.f @@ -263,17 +263,8 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_tot_num,mo_tot_num) BEGIN_DOC ! Fock matrix on the MO basis END_DOC - double precision, allocatable :: T(:,:) - allocate ( T(ao_num,mo_tot_num) ) - call dgemm('N','N', ao_num, mo_tot_num, ao_num, & - 1.d0, Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), & - mo_coef, size(mo_coef,1), & - 0.d0, T, size(T,1)) - call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, & - 1.d0, mo_coef,size(mo_coef,1), & - T, size(T,1), & - 0.d0, Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1)) - deallocate(T) + call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), & + Fock_matrix_mo_alpha,size(Fock_matrix_mo_alpha,1)) END_PROVIDER @@ -282,17 +273,8 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_tot_num,mo_tot_num) BEGIN_DOC ! Fock matrix on the MO basis END_DOC - double precision, allocatable :: T(:,:) - allocate ( T(ao_num,mo_tot_num) ) - call dgemm('N','N', ao_num, mo_tot_num, ao_num, & - 1.d0, Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), & - mo_coef, size(mo_coef,1), & - 0.d0, T, size(T,1)) - call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, & - 1.d0, mo_coef,size(mo_coef,1), & - T, size(T,1), & - 0.d0, Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1)) - deallocate(T) + call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), & + Fock_matrix_mo_beta,size(Fock_matrix_mo_beta,1)) END_PROVIDER BEGIN_PROVIDER [ double precision, HF_energy ] @@ -330,97 +312,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ] enddo enddo else - double precision, allocatable :: T(:,:), M(:,:) - integer :: ierr - ! F_ao = S C F_mo C^t S - allocate (T(ao_num,ao_num),M(ao_num,ao_num),stat=ierr) - if (ierr /=0 ) then - print *, irp_here, ' : allocation failed' - endif - -! ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num) -! -> M(ao_num,mo_tot_num) - call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0, & - ao_overlap, size(ao_overlap,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - -! M(ao_num,mo_tot_num) . Fock_matrix_mo (mo_tot_num,mo_tot_num) -! -> T(ao_num,mo_tot_num) - call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & - M, size(M,1), & - Fock_matrix_mo, size(Fock_matrix_mo,1), & - 0.d0, & - T, size(T,1)) - -! T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num) -! -> M(ao_num,ao_num) - call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0, & - T, size(T,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - -! M(ao_num,ao_num) . ao_overlap (ao_num,ao_num) -! -> Fock_matrix_ao(ao_num,ao_num) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & - M, size(M,1), & - ao_overlap, size(ao_overlap,1), & - 0.d0, & - Fock_matrix_ao, size(Fock_matrix_ao,1)) - - - deallocate(T) + call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), & + Fock_matrix_ao,size(Fock_matrix_ao,1)) endif END_PROVIDER -subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) - implicit none - integer, intent(in) :: LDFMO ! size(FMO,1) - integer, intent(in) :: LDFAO ! size(FAO,1) - double precision, intent(in) :: FMO(LDFMO,*) - double precision, intent(out) :: FAO(LDFAO,*) - - double precision, allocatable :: T(:,:), M(:,:) - integer :: ierr - ! F_ao = S C F_mo C^t S - allocate (T(ao_num,ao_num),M(ao_num,ao_num),stat=ierr) - if (ierr /=0 ) then - print *, irp_here, ' : allocation failed' - endif - -! ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num) -! -> M(ao_num,mo_tot_num) - call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0, & - ao_overlap, size(ao_overlap,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - -! M(ao_num,mo_tot_num) . FMO (mo_tot_num,mo_tot_num) -! -> T(ao_num,mo_tot_num) - call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & - M, size(M,1), & - FMO, size(FMO,1), & - 0.d0, & - T, size(T,1)) - -! T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num) -! -> M(ao_num,ao_num) - call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0, & - T, size(T,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - -! M(ao_num,ao_num) . ao_overlap (ao_num,ao_num) -! -> Fock_matrix_ao(ao_num,ao_num) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & - M, size(M,1), & - ao_overlap, size(ao_overlap,1), & - 0.d0, & - FAO, size(FAO,1)) - deallocate(T,M) -end diff --git a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f index 9e0d29e5..736d1a97 100644 --- a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f +++ b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f @@ -1,7 +1,7 @@ BEGIN_PROVIDER [double precision, HF_density_matrix_ao_alpha, (ao_num,ao_num) ] implicit none BEGIN_DOC - ! S^-1 x Alpha density matrix in the AO basis x S^-1 + ! S^{-1}.P_alpha.S^{-1} END_DOC call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & @@ -14,7 +14,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num,ao_num) ] implicit none BEGIN_DOC - ! S^-1 Beta density matrix in the AO basis x S^-1 + ! S^{-1}.P_beta.S^{-1} END_DOC call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & @@ -27,7 +27,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num,ao_num) ] implicit none BEGIN_DOC - ! S^-1 Density matrix in the AO basis S^-1 + ! S^{-1}.P.S^{-1} where P = C.C^t END_DOC ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1)) if (elec_alpha_num== elec_beta_num) then diff --git a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f index 5d8097d9..241721ae 100644 --- a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f +++ b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f @@ -133,7 +133,7 @@ END_DOC write(output_hartree_fock,*) if(.not.no_oa_or_av_opt)then - call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.) endif call write_double(output_hartree_fock, Energy_SCF, 'Hartree-Fock energy') diff --git a/plugins/Hartree_Fock/SCF.irp.f b/plugins/Hartree_Fock/SCF.irp.f index 5336e8e0..3d71d826 100644 --- a/plugins/Hartree_Fock/SCF.irp.f +++ b/plugins/Hartree_Fock/SCF.irp.f @@ -23,7 +23,7 @@ subroutine create_guess mo_coef = ao_ortho_lowdin_coef TOUCH mo_coef mo_label = 'Guess' - call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label) + call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.) SOFT_TOUCH mo_coef mo_label else if (mo_guess_type == "Huckel") then call huckel_guess diff --git a/plugins/Hartree_Fock/damping_SCF.irp.f b/plugins/Hartree_Fock/damping_SCF.irp.f index 2be14ed3..20a8abd7 100644 --- a/plugins/Hartree_Fock/damping_SCF.irp.f +++ b/plugins/Hartree_Fock/damping_SCF.irp.f @@ -119,7 +119,7 @@ subroutine damping_SCF write(output_hartree_fock,*) if(.not.no_oa_or_av_opt)then - call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.) endif call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy') diff --git a/plugins/Hartree_Fock/huckel.irp.f b/plugins/Hartree_Fock/huckel.irp.f index 9cb84507..c9e32ad5 100644 --- a/plugins/Hartree_Fock/huckel.irp.f +++ b/plugins/Hartree_Fock/huckel.irp.f @@ -7,25 +7,28 @@ subroutine huckel_guess double precision :: accu double precision :: c character*(64) :: label - + double precision, allocatable :: A(:,:) label = "Guess" - call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & - size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label,1) - TOUCH mo_coef - c = 0.5d0 * 1.75d0 + allocate (A(ao_num, ao_num)) + A = 0.d0 do j=1,ao_num do i=1,ao_num - Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral_diag(i) + & - ao_mono_elec_integral_diag(j)) + A(i,j) = c * ao_overlap(i,j) * (ao_mono_elec_integral_diag(i) + ao_mono_elec_integral_diag(j)) enddo - Fock_matrix_ao(j,j) = Fock_matrix_ao_alpha(j,j) + A(j,j) = ao_mono_elec_integral_diag(j) + ao_bi_elec_integral_alpha(j,j) enddo - TOUCH Fock_matrix_ao + + Fock_matrix_ao_alpha(1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num) + Fock_matrix_ao_beta (1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num) + +! TOUCH mo_coef + + TOUCH Fock_matrix_ao_alpha Fock_matrix_ao_beta mo_coef = eigenvectors_fock_matrix_mo SOFT_TOUCH mo_coef call save_mos + deallocate(A) end diff --git a/plugins/MRPT/jm_mrpt2.irp.f b/plugins/MRPT/jm_mrpt2.irp.f index 89deef2e..df9a1dbc 100644 --- a/plugins/MRPT/jm_mrpt2.irp.f +++ b/plugins/MRPT/jm_mrpt2.irp.f @@ -11,7 +11,7 @@ end subroutine routine_3 implicit none !provide fock_virt_total_spin_trace - provide delta_ij + provide delta_ij_mrpt print *, 'N_det = ', N_det print *, 'N_states = ', N_states diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 83f087bc..34d26127 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -1,5 +1,5 @@ - BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] + BEGIN_PROVIDER [ double precision, delta_ij_mrpt, (N_det,N_det,N_states) ] &BEGIN_PROVIDER [ double precision, second_order_pt_new, (N_states) ] &BEGIN_PROVIDER [ double precision, second_order_pt_new_1h, (N_states) ] &BEGIN_PROVIDER [ double precision, second_order_pt_new_1p, (N_states) ] @@ -19,7 +19,7 @@ double precision, allocatable :: delta_ij_tmp(:,:,:) - delta_ij = 0.d0 + delta_ij_mrpt = 0.d0 allocate (delta_ij_tmp(N_det,N_det,N_states)) @@ -32,7 +32,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_1h(i_state) = accu(i_state) @@ -47,7 +47,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_1p(i_state) = accu(i_state) @@ -63,7 +63,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_1h1p(i_state) = accu(i_state) @@ -79,7 +79,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_1h1p(i_state) = accu(i_state) @@ -95,7 +95,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_2h(i_state) = accu(i_state) @@ -110,7 +110,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_2p(i_state) = accu(i_state) @@ -126,7 +126,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_1h2p(i_state) = accu(i_state) @@ -142,7 +142,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_2h1p(i_state) = accu(i_state) @@ -157,7 +157,7 @@ !do i = 1, N_det ! do j = 1, N_det ! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) -! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) +! delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) ! enddo !enddo !second_order_pt_new_2h2p(i_state) = accu(i_state) @@ -168,7 +168,7 @@ call give_2h2p(contrib_2h2p) do i_state = 1, N_states do i = 1, N_det - delta_ij(i,i,i_state) += contrib_2h2p(i_state) + delta_ij_mrpt(i,i,i_state) += contrib_2h2p(i_state) enddo second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) enddo @@ -179,9 +179,9 @@ accu = 0.d0 do i_state = 1, N_states do i = 1, N_det -! write(*,'(1000(F16.10,x))')delta_ij(i,:,:) +! write(*,'(1000(F16.10,x))')delta_ij_mrpt(i,:,:) do j = i_state, N_det - accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + accu(i_state) += delta_ij_mrpt(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) enddo enddo second_order_pt_new(i_state) = accu(i_state) @@ -199,7 +199,7 @@ END_PROVIDER do i_state = 1, N_states do i = 1,N_det do j = 1,N_det - Hmatrix_dressed_pt2_new(j,i,i_state) = H_matrix_all_dets(j,i) + delta_ij(j,i,i_state) + Hmatrix_dressed_pt2_new(j,i,i_state) = H_matrix_all_dets(j,i) + delta_ij_mrpt(j,i,i_state) enddo enddo enddo @@ -214,7 +214,7 @@ END_PROVIDER do i = 1,N_det do j = i,N_det Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) = H_matrix_all_dets(j,i) & - + 0.5d0 * ( delta_ij(j,i,i_state) + delta_ij(i,j,i_state) ) + + 0.5d0 * ( delta_ij_mrpt(j,i,i_state) + delta_ij_mrpt(i,j,i_state) ) Hmatrix_dressed_pt2_new_symmetrized(i,j,i_state) = Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) enddo enddo diff --git a/plugins/analyze_wf/attachment.irp.f b/plugins/analyze_wf/attachment.irp.f new file mode 100644 index 00000000..d086aa21 --- /dev/null +++ b/plugins/analyze_wf/attachment.irp.f @@ -0,0 +1,64 @@ + BEGIN_PROVIDER [ double precision, one_body_dm_mo_detachment, (mo_tot_num,mo_tot_num,2:N_states) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_attachment, (mo_tot_num,mo_tot_num,2:N_states) ] + implicit none + BEGIN_DOC + ! Detachment and attachment density matrices in MO basis + END_DOC + integer :: i,j, k, istate + double precision :: km(mo_tot_num), kp(mo_tot_num) + + one_body_dm_mo_detachment = 0.d0 + one_body_dm_mo_attachment = 0.d0 + + do istate=2,N_states + + km(:) = 0.d0 + kp(:) = 0.d0 + do i=1,mo_tot_num + if (one_body_dm_mo_diff_eigvalues(i,istate) < 0) then + km(i) = -one_body_dm_mo_diff_eigvalues(i,istate) + else + kp(i) = one_body_dm_mo_diff_eigvalues(i,istate) + endif + enddo + + ! Attachment + do k=1,mo_tot_num + do j=1,mo_tot_num + do i=1,mo_tot_num + one_body_dm_mo_detachment(i,j,istate) = one_body_dm_mo_detachment(i,j,istate) + & + one_body_dm_mo_diff_eigvectors(i,k,istate)*km(k)* & + one_body_dm_mo_diff_eigvectors(j,k,istate) + one_body_dm_mo_attachment(i,j,istate) = one_body_dm_mo_attachment(i,j,istate) + & + one_body_dm_mo_diff_eigvectors(i,k,istate)*kp(k)* & + one_body_dm_mo_diff_eigvectors(j,k,istate) + enddo + enddo + enddo + + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, one_body_dm_ao_detachment, (ao_num,ao_num,2:N_states) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_ao_attachment, (ao_num,ao_num,2:N_states) ] + implicit none + BEGIN_DOC + ! Detachment and attachment density matrices in AO basis + END_DOC + integer :: istate + do istate=2,N_states + call mo_to_ao_no_overlap( & + one_body_dm_mo_attachment(1,1,istate), & + size(one_body_dm_mo_attachment,1), & + one_body_dm_ao_attachment(1,1,istate), & + size(one_body_dm_ao_attachment,1) ) + call mo_to_ao_no_overlap( & + one_body_dm_mo_detachment(1,1,istate), & + size(one_body_dm_mo_detachment,1), & + one_body_dm_ao_detachment(1,1,istate), & + size(one_body_dm_ao_detachment,1) ) + enddo + +END_PROVIDER + diff --git a/plugins/analyze_wf/dump_nto.irp.f b/plugins/analyze_wf/dump_nto.irp.f new file mode 100644 index 00000000..8d19c3eb --- /dev/null +++ b/plugins/analyze_wf/dump_nto.irp.f @@ -0,0 +1,19 @@ +program dump_nto + implicit none + integer :: i,j, istate + + print *, 'Phi_S' + do i=2,N_states + print *, i, Phi_S(i) + enddo + + do istate=2,N_states + do j=1,mo_tot_num + print *, 'MO: ', j, 'State:', istate, 'Eig:', one_body_dm_mo_diff_eigvalues(j,istate) + do i=1,ao_num + print *, i, transition_natorb(i,j,istate) + enddo + enddo + enddo + +end diff --git a/plugins/analyze_wf/dump_one_body_mos.irp.f b/plugins/analyze_wf/dump_one_body_mos.irp.f new file mode 100644 index 00000000..7ab841ef --- /dev/null +++ b/plugins/analyze_wf/dump_one_body_mos.irp.f @@ -0,0 +1,38 @@ +program dump_one_body_mos + implicit none + BEGIN_DOC +! Output density matrices of all the states + END_DOC + read_wf = .True. + TOUCH read_wf + call run() +end + +subroutine run + implicit none + integer :: istate + integer, parameter :: iunit = 66 + character*(64) :: filename, fmt + integer :: i,j,k + + write(fmt,'(''('',I4.4,''(X,E20.14))'')') mo_tot_num + do istate=1,N_states + write(filename,'(''state.'',I2.2)') istate + open(unit=iunit, form='formatted', file=filename) + write(iunit,*) mo_tot_num + do j=1,mo_tot_num + write(iunit,fmt) one_body_dm_mo_alpha(1:mo_tot_num,j,istate) + one_body_dm_mo_beta(1:mo_tot_num,j,istate) + enddo + enddo + + + call run2() +end +subroutine run2 + integer :: i,j, istate + print *, 'Phi_S' + do i=2,N_states + print *, i, Phi_S(i) + enddo +end + diff --git a/plugins/analyze_wf/phi_s.irp.f b/plugins/analyze_wf/phi_s.irp.f new file mode 100644 index 00000000..12bdb970 --- /dev/null +++ b/plugins/analyze_wf/phi_s.irp.f @@ -0,0 +1,128 @@ + BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvalues, (mo_tot_num, 2:N_states) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvectors, (mo_tot_num, mo_tot_num, 2:N_states) ] + implicit none + BEGIN_DOC + ! Eigenvalues and eigenvectors of one_body_dm_mo_diff + END_DOC + integer :: i,j,istate + integer :: liwork, lwork, n, info + integer, allocatable :: iwork(:) + double precision, allocatable :: work(:) + + + one_body_dm_mo_diff_eigvectors(1:mo_tot_num, 1:mo_tot_num, 2:N_states) =& + one_body_dm_mo_diff(1:mo_tot_num, 1:mo_tot_num, 2:N_states) + + n = mo_tot_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n + + allocate(work(lwork), iwork(liwork)) + + lwork=-1 + liwork=-1 + istate=2 + + call dsyevd( 'V', 'U', mo_tot_num, & + one_body_dm_mo_diff_eigvectors(1,1,istate), & + size(one_body_dm_mo_diff_eigvectors,1), & + one_body_dm_mo_diff_eigvalues(1,istate), & + work, lwork, iwork, liwork, info) + + + if (info /= 0) then + print *, irp_here//' DSYEVD failed : ', info + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(iwork,work) + + allocate(work(lwork), iwork(liwork)) + + do istate=2,N_states + call dsyevd( 'V', 'U', mo_tot_num, & + one_body_dm_mo_diff_eigvectors(1,1,istate), & + size(one_body_dm_mo_diff_eigvectors,1), & + one_body_dm_mo_diff_eigvalues(1,istate), & + work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' DSYEVD failed : ', info + stop 1 + endif + + enddo + + deallocate(iwork,work) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, transition_natorb, (ao_num,mo_tot_num,2:N_states) ] + implicit none + BEGIN_DOC + ! Natural transition molecular orbitals + END_DOC + + integer :: istate + + do istate=2,N_states + call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, & + mo_coef, size(mo_coef,1), & + one_body_dm_mo_diff_eigvectors(1,1,istate), & + size(one_body_dm_mo_diff_eigvectors,1), 0.d0, & + transition_natorb(1,1,istate), size(transition_natorb,1)) + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, phi_s, (2:N_states) ] + implicit none + BEGIN_DOC + ! + END_DOC + + integer :: i,istate + double precision, allocatable :: T(:,:), A(:,:), D(:,:) + double precision :: trace, norm + allocate(T(ao_num,ao_num), A(ao_num,ao_num), D(ao_num,ao_num)) + + do istate=2,N_states + + call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, & + S_half, size(S_half,1), & + one_body_dm_ao_attachment(1,1,istate), size(one_body_dm_ao_attachment,1), 0.d0,& + T, size(T,1)) + call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, & + T, size(T,1), & + S_half, size(S_half,1), 0.d0, & + A, size(A,1)) +! + call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, & + S_half, size(S_half,1), & + one_body_dm_ao_detachment(1,1,istate), size(one_body_dm_ao_detachment,1), 0.d0,& + T, size(T,1)) + call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, & + T, size(T,1), & + S_half, size(S_half,1), 0.d0, & + D, size(D,1)) + + trace = 0.d0 + do i=1,ao_num + trace = trace + A(i,i) + enddo + norm = 0.d0 + do i=1,ao_num + norm = norm + D(i,i) + enddo + norm = 0.5d0*(norm + trace) + + trace = 0.d0 + do i=1,mo_tot_num + trace = trace + dsqrt(A(i,i)*D(i,i)) + enddo + phi_s(istate) = trace/norm + enddo + +END_PROVIDER + diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f index f2630520..53fa74a1 100644 --- a/plugins/mrcepa0/mrcepa0_general.irp.f +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -38,12 +38,12 @@ subroutine run(N_st,energy) do while (delta_E > thresh_mrcc) iteration += 1 print *, '===============================================' - print *, 'MRCEPA0 Iteration', iteration, '/', n_it_mrcc_max + print *, 'Iteration', iteration, '/', n_it_mrcc_max print *, '===============================================' print *, '' E_old = sum(ci_energy_dressed(1:N_states)) do i=1,N_st - call write_double(6,ci_energy_dressed(i),"MRCEPA0 energy") + call write_double(6,ci_energy_dressed(i),"Energy") enddo call diagonalize_ci_dressed(lambda) E_new = sum(ci_energy_dressed(1:N_states)) @@ -61,7 +61,7 @@ subroutine run(N_st,energy) exit endif enddo - call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy") + call write_double(6,ci_energy_dressed(1),"Final energy") endif energy(1:N_st) = ci_energy_dressed(1:N_st) end diff --git a/plugins/read_integral/read_integrals_mo.irp.f b/plugins/read_integral/read_integrals_mo.irp.f index e1ff5fe8..5376b2a2 100644 --- a/plugins/read_integral/read_integrals_mo.irp.f +++ b/plugins/read_integral/read_integrals_mo.irp.f @@ -1,5 +1,10 @@ program read_integrals - + BEGIN_DOC +! Reads the integrals from the following files: +! - kinetic_mo +! - nuclear_mo +! - bielec_mo + END_DOC PROVIDE ezfio_filename call ezfio_set_integrals_monoelec_disk_access_mo_one_integrals("None") call run diff --git a/scripts/ezfio_interface/ei_handler.py b/scripts/ezfio_interface/ei_handler.py index 4d61062e..ee44a1e1 100755 --- a/scripts/ezfio_interface/ei_handler.py +++ b/scripts/ezfio_interface/ei_handler.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- """ -Welcom the ei_handler. +Welcome to the ei_handler. We will create all the ezfio related stuff from a EZFIO.cfg file. Usage: @@ -543,7 +543,7 @@ def create_ocaml_input(dict_ezfio_cfg, module_lower): template += ["open Qptypes;;", "open Qputils;;", - "open Core.Std;;", + "open Core;;", "", "module {0} : sig".format(module_lower.capitalize())] diff --git a/scripts/ezfio_interface/ezfio_generate_ocaml.py b/scripts/ezfio_interface/ezfio_generate_ocaml.py index f36441b6..244f67a3 100755 --- a/scripts/ezfio_interface/ezfio_generate_ocaml.py +++ b/scripts/ezfio_interface/ezfio_generate_ocaml.py @@ -1,10 +1,10 @@ #!/usr/bin/env python """ -This programme generate all the -ocaml template needed by qp_edit +This program generates all the +OCaml templates needed by qp_edit You can see `ezfio_generate_provider.py` -for an example of utilisation +for an example. """ import sys @@ -170,7 +170,7 @@ class EZFIO_ocaml(object): else: l_template += [" {0:<30} : {1};".format(p, t.ocaml)] - l_template += [" } with sexp", + l_template += [" } [@@deriving sexp]", ";;"] # ~#~#~#~#~#~ # @@ -352,7 +352,7 @@ class EZFIO_ocaml(object): l_template = ['open Qputils;;', 'open Qptypes;;', - 'open Core.Std;;', + 'open Core;;', ''] for m in self.l_module_lower: diff --git a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py index 6de221f3..7f4f30be 100755 --- a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py +++ b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py @@ -3,7 +3,7 @@ convert output of gamess/GAU$$IAN to ezfio Usage: - qp_convert_output_to_ezfio.py [--ezfio=] + qp_convert_output_to_ezfio.py [-o ] Option: file.out is the file to check (like gamess.out) @@ -429,8 +429,8 @@ if __name__ == '__main__': file_ = get_full_path(arguments['']) - if arguments["--ezfio"]: - ezfio_file = get_full_path(arguments["--ezfio"]) + if arguments["-o"]: + ezfio_file = get_full_path(arguments[""]) else: ezfio_file = "{0}.ezfio".format(file_) diff --git a/scripts/ezfio_interface/qp_edit_template b/scripts/ezfio_interface/qp_edit_template index 9d327124..55a35f83 100644 --- a/scripts/ezfio_interface/qp_edit_template +++ b/scripts/ezfio_interface/qp_edit_template @@ -4,7 +4,7 @@ open Qputils open Qptypes -open Core.Std +open Core (** Interactive editing of the input. diff --git a/src/AO_Basis/EZFIO.cfg b/src/AO_Basis/EZFIO.cfg index 9e548514..07712a4a 100644 --- a/src/AO_Basis/EZFIO.cfg +++ b/src/AO_Basis/EZFIO.cfg @@ -16,7 +16,7 @@ interface: ezfio, provider [ao_prim_num_max] type: integer -doc: number of primitive maximun +doc: maximum number of primitives default: =maxval(ao_basis.ao_prim_num) interface: ezfio @@ -76,3 +76,4 @@ size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio default: false + diff --git a/src/AO_Basis/ao_overlap.irp.f b/src/AO_Basis/ao_overlap.irp.f index edf48b25..83110293 100644 --- a/src/AO_Basis/ao_overlap.irp.f +++ b/src/AO_Basis/ao_overlap.irp.f @@ -129,3 +129,104 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ] !$OMP END PARALLEL DO END_PROVIDER +BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ] + implicit none + BEGIN_DOC +! S^-1 + END_DOC + call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv,size(S_inv,1)) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ] + + BEGIN_DOC +! Matrix X = S^{-1/2} obtained by SVD + END_DOC + + implicit none + + integer :: num_linear_dependencies + integer :: LDA, LDC + double precision, allocatable :: U(:,:),Vt(:,:), D(:) + integer :: info, i, j, k + double precision, parameter :: threshold_overlap_AO_eigenvalues = 1.d-6 + + LDA = size(AO_overlap,1) + LDC = size(S_half_inv,1) + + allocate( & + U(LDC,AO_num), & + Vt(LDA,AO_num), & + D(AO_num)) + + call svd( & + AO_overlap,LDA, & + U,LDC, & + D, & + Vt,LDA, & + AO_num,AO_num) + + num_linear_dependencies = 0 + do i=1,AO_num + print*,D(i) + if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then + D(i) = 0.d0 + num_linear_dependencies += 1 + else + ASSERT (D(i) > 0.d0) + D(i) = 1.d0/sqrt(D(i)) + endif + do j=1,AO_num + S_half_inv(j,i) = 0.d0 + enddo + enddo + write(*,*) 'linear dependencies',num_linear_dependencies + + do k=1,AO_num + if(D(k) /= 0.d0) then + do j=1,AO_num + do i=1,AO_num + S_half_inv(i,j) = S_half_inv(i,j) + U(i,k)*D(k)*Vt(k,j) + enddo + enddo + endif + enddo + + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num) ] + implicit none + BEGIN_DOC + ! S^{1/2} + END_DOC + + integer :: i,j,k + double precision, allocatable :: U(:,:) + double precision, allocatable :: Vt(:,:) + double precision, allocatable :: D(:) + + allocate(U(ao_num,ao_num),Vt(ao_num,ao_num),D(ao_num)) + + call svd(ao_overlap,size(ao_overlap,1),U,size(U,1),D,Vt,size(Vt,1),ao_num,ao_num) + + do i=1,ao_num + D(i) = dsqrt(D(i)) + do j=1,ao_num + S_half(j,i) = 0.d0 + enddo + enddo + + do k=1,ao_num + do j=1,ao_num + do i=1,ao_num + S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j) + enddo + enddo + enddo + + deallocate(U,Vt,D) + +END_PROVIDER + diff --git a/src/AO_Basis/aos.irp.f b/src/AO_Basis/aos.irp.f index 062ef296..871c0ee6 100644 --- a/src/AO_Basis/aos.irp.f +++ b/src/AO_Basis/aos.irp.f @@ -9,7 +9,7 @@ BEGIN_PROVIDER [ integer, ao_num_align ] ao_num_align = align_double(ao_num) END_PROVIDER - BEGIN_PROVIDER [ integer, ao_prim_num_max ] +BEGIN_PROVIDER [ integer, ao_prim_num_max ] implicit none ao_prim_num_max = 0 PROVIDE ezfio_filename diff --git a/src/Davidson/print_energy.irp.f b/src/Davidson/print_energy.irp.f new file mode 100644 index 00000000..ae6f1da2 --- /dev/null +++ b/src/Davidson/print_energy.irp.f @@ -0,0 +1,22 @@ +program print_energy + implicit none + read_wf = .true. + touch read_wf + call routine +end + +subroutine routine + implicit none + integer :: i,j + double precision :: accu,hij + + print*, 'psi_energy = ',psi_energy + nuclear_repulsion + accu = 0.d0 +! do i = 1,N_det +! do j = 1,N_det +! call i_H_j(psi_det(1,1,j),psi_det(1,1,i),N_int,hij) +! accu += psi_coef(i,1) * psi_coef(j,1) * hij +! enddo +! enddo +! print*, 'accu = ',accu + nuclear_repulsion +end diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index 9d0512f4..f4c5d866 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -2,7 +2,7 @@ type: Det_number_max doc: Max number of determinants in the wave function interface: ezfio,provider,ocaml -default: 10000 +default: 1000000 [N_det_max_property] type: Det_number_max @@ -14,7 +14,7 @@ default: 10000 type: Det_number_max doc: Maximum number of determinants diagonalized by Jacobi interface: ezfio,provider,ocaml -default: 1000 +default: 2000 [N_states] type: States_number diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index e0764d96..4fa33daa 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -15,6 +15,32 @@ enddo END_PROVIDER +BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff, (mo_tot_num,mo_tot_num,2:N_states) ] + implicit none + BEGIN_DOC + ! Difference of the one-body density matrix with respect to the ground state + END_DOC + integer :: i,j, istate + + do istate=2,N_states + do j=1,mo_tot_num + do i=1,mo_tot_num + one_body_dm_mo_diff(i,j,istate) = & + one_body_dm_mo_alpha(i,j,istate) - one_body_dm_mo_alpha(i,j,1) + & + one_body_dm_mo_beta (i,j,istate) - one_body_dm_mo_beta (i,j,1) + enddo + enddo + double precision :: trace + trace = 0.d0 + do i=1,mo_tot_num + trace += one_body_dm_mo_diff(i,i,istate) + enddo + print *, irp_here, trace + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, one_body_dm_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ] implicit none integer :: i,j,ispin,istate diff --git a/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f index f24e4974..1333753d 100644 --- a/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f +++ b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f @@ -142,8 +142,8 @@ subroutine ao_bielec_integrals_in_map_collector integer(ZMQ_PTR), external :: new_zmq_pull_socket integer(ZMQ_PTR) :: zmq_socket_pull - integer*8 :: control, accu - integer :: task_id, more, sze + integer*8 :: control, accu, sze + integer :: task_id, more zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() zmq_socket_pull = new_zmq_pull_socket() diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 190615c3..996f8464 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -327,41 +327,50 @@ subroutine insert_into_mo_integrals_map(n_integrals, & call map_update(mo_integrals_map, buffer_i, buffer_values, n_integrals, thr) end - BEGIN_PROVIDER [ integer, mo_integrals_cache_min ] -&BEGIN_PROVIDER [ integer, mo_integrals_cache_max ] + BEGIN_PROVIDER [ integer*4, mo_integrals_cache_min ] +&BEGIN_PROVIDER [ integer*4, mo_integrals_cache_max ] +&BEGIN_PROVIDER [ integer*8, mo_integrals_cache_min_8 ] +&BEGIN_PROVIDER [ integer*8, mo_integrals_cache_max_8 ] implicit none BEGIN_DOC ! Min and max values of the MOs for which the integrals are in the cache END_DOC - mo_integrals_cache_min = max(1,elec_alpha_num - 31) - mo_integrals_cache_max = min(mo_tot_num,mo_integrals_cache_min+63) + mo_integrals_cache_min_8 = max(1_8,elec_alpha_num - 63_8) + mo_integrals_cache_max_8 = min(int(mo_tot_num,8),mo_integrals_cache_min+127_8) + mo_integrals_cache_min = max(1,elec_alpha_num - 63) + mo_integrals_cache_max = min(mo_tot_num,mo_integrals_cache_min+127) END_PROVIDER -BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0:64*64*64*64) ] +BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:128_8*128_8*128_8*128_8) ] implicit none BEGIN_DOC ! Cache of MO integrals for fast access END_DOC PROVIDE mo_bielec_integrals_in_map - integer :: i,j,k,l - integer :: ii + integer*8 :: i,j,k,l + integer*4 :: i4,j4,k4,l4 + integer*8 :: ii integer(key_kind) :: idx real(integral_kind) :: integral FREE ao_integrals_cache - !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) - do l=mo_integrals_cache_min,mo_integrals_cache_max - do k=mo_integrals_cache_min,mo_integrals_cache_max - do j=mo_integrals_cache_min,mo_integrals_cache_max - do i=mo_integrals_cache_min,mo_integrals_cache_max + !$OMP PARALLEL DO PRIVATE (i,j,k,l,i4,j4,k4,l4,idx,ii,integral) + do l=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + l4 = int(l,4) + do k=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + k4 = int(k,4) + do j=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + j4 = int(j,4) + do i=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + i4 = int(i,4) !DIR$ FORCEINLINE - call bielec_integrals_index(i,j,k,l,idx) + call bielec_integrals_index(i4,j4,k4,l4,idx) !DIR$ FORCEINLINE call map_get(mo_integrals_map,idx,integral) - ii = l-mo_integrals_cache_min - ii = ior( ishft(ii,6), k-mo_integrals_cache_min) - ii = ior( ishft(ii,6), j-mo_integrals_cache_min) - ii = ior( ishft(ii,6), i-mo_integrals_cache_min) + ii = l-mo_integrals_cache_min_8 + ii = ior( ishft(ii,7), k-mo_integrals_cache_min_8) + ii = ior( ishft(ii,7), j-mo_integrals_cache_min_8) + ii = ior( ishft(ii,7), i-mo_integrals_cache_min_8) mo_integrals_cache(ii) = integral enddo enddo @@ -381,6 +390,7 @@ double precision function get_mo_bielec_integral(i,j,k,l,map) integer, intent(in) :: i,j,k,l integer(key_kind) :: idx integer :: ii + integer*8 :: ii_8 type(map_type), intent(inout) :: map real(integral_kind) :: tmp PROVIDE mo_bielec_integrals_in_map mo_integrals_cache @@ -388,18 +398,18 @@ double precision function get_mo_bielec_integral(i,j,k,l,map) ii = ior(ii, k-mo_integrals_cache_min) ii = ior(ii, j-mo_integrals_cache_min) ii = ior(ii, i-mo_integrals_cache_min) - if (iand(ii, -64) /= 0) then + if (iand(ii, -128) /= 0) then !DIR$ FORCEINLINE call bielec_integrals_index(i,j,k,l,idx) !DIR$ FORCEINLINE call map_get(map,idx,tmp) get_mo_bielec_integral = dble(tmp) else - ii = l-mo_integrals_cache_min - ii = ior( ishft(ii,6), k-mo_integrals_cache_min) - ii = ior( ishft(ii,6), j-mo_integrals_cache_min) - ii = ior( ishft(ii,6), i-mo_integrals_cache_min) - get_mo_bielec_integral = mo_integrals_cache(ii) + ii_8 = int(l,8)-mo_integrals_cache_min_8 + ii_8 = ior( ishft(ii_8,7), int(k,8)-mo_integrals_cache_min_8) + ii_8 = ior( ishft(ii_8,7), int(j,8)-mo_integrals_cache_min_8) + ii_8 = ior( ishft(ii_8,7), int(i,8)-mo_integrals_cache_min_8) + get_mo_bielec_integral = mo_integrals_cache(ii_8) endif end diff --git a/src/MOGuess/H_CORE_guess.irp.f b/src/MOGuess/H_CORE_guess.irp.f index d3e2eef9..11fef327 100644 --- a/src/MOGuess/H_CORE_guess.irp.f +++ b/src/MOGuess/H_CORE_guess.irp.f @@ -1,15 +1,9 @@ -program H_CORE_guess +program H_CORE_guess_prog BEGIN_DOC ! Produce `H_core` MO orbital ! output: mo_basis.mo_tot_num mo_basis.mo_label mo_basis.ao_md5 mo_basis.mo_coef mo_basis.mo_occ END_DOC implicit none character*(64) :: label - label = "Guess" - call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & - size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label,1) - print *, 'save mos' - call save_mos - + call hcore_guess end diff --git a/src/MOGuess/h_core_guess_routine.irp.f b/src/MOGuess/h_core_guess_routine.irp.f index 23899160..5b4ede91 100644 --- a/src/MOGuess/h_core_guess_routine.irp.f +++ b/src/MOGuess/h_core_guess_routine.irp.f @@ -7,8 +7,7 @@ subroutine hcore_guess label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label,1) - print *, 'save mos' + size(mo_mono_elec_integral,2),label,1,.false.) call save_mos SOFT_TOUCH mo_coef mo_label end diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index 8962dd00..4bcd7221 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -35,7 +35,6 @@ END_PROVIDER END_DOC integer :: i, j double precision, allocatable :: buffer(:,:) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: buffer logical :: exists PROVIDE ezfio_filename @@ -153,16 +152,17 @@ BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ] endif END_PROVIDER + subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo) implicit none BEGIN_DOC ! Transform A from the AO basis to the MO basis ! - ! C.A_ao.Ct + ! Ct.A_ao.C END_DOC integer, intent(in) :: LDA_ao,LDA_mo - double precision, intent(in) :: A_ao(LDA_ao) - double precision, intent(out) :: A_mo(LDA_mo) + double precision, intent(in) :: A_ao(LDA_ao,ao_num) + double precision, intent(out) :: A_mo(LDA_mo,mo_tot_num) double precision, allocatable :: T(:,:) allocate ( T(ao_num_align,mo_tot_num) ) @@ -189,54 +189,48 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao) ! (S.C).A_mo.(S.C)t END_DOC integer, intent(in) :: LDA_ao,LDA_mo - double precision, intent(in) :: A_mo(LDA_mo) - double precision, intent(out) :: A_ao(LDA_ao) - double precision, allocatable :: T(:,:), SC(:,:) + double precision, intent(in) :: A_mo(LDA_mo,mo_tot_num) + double precision, intent(out) :: A_ao(LDA_ao,ao_num) + double precision, allocatable :: T(:,:) - allocate ( SC(ao_num_align,mo_tot_num) ) allocate ( T(mo_tot_num_align,ao_num) ) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - - call dgemm('N','N', ao_num, mo_tot_num, ao_num, & - 1.d0, ao_overlap,size(ao_overlap,1), & - mo_coef, size(mo_coef,1), & - 0.d0, SC, ao_num_align) call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & - 1.d0, A_mo,LDA_mo, & - SC, size(SC,1), & - 0.d0, T, mo_tot_num_align) + 1.d0, A_mo,size(A_mo,1), & + S_mo_coef, size(S_mo_coef,1), & + 0.d0, T, size(T,1)) call dgemm('N','N', ao_num, ao_num, mo_tot_num, & - 1.d0, SC,size(SC,1), & - T, mo_tot_num_align, & - 0.d0, A_ao, LDA_ao) + 1.d0, S_mo_coef, size(S_mo_coef,1), & + T, size(T,1), & + 0.d0, A_ao, size(A_ao,1)) - deallocate(T,SC) + deallocate(T) end subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao) implicit none BEGIN_DOC ! Transform A from the MO basis to the S^-1 AO basis + ! Useful for density matrix END_DOC integer, intent(in) :: LDA_ao,LDA_mo - double precision, intent(in) :: A_mo(LDA_mo) - double precision, intent(out) :: A_ao(LDA_ao) + double precision, intent(in) :: A_mo(LDA_mo,mo_tot_num) + double precision, intent(out) :: A_ao(LDA_ao,ao_num) double precision, allocatable :: T(:,:) allocate ( T(mo_tot_num_align,ao_num) ) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & - 1.d0, A_mo,LDA_mo, & + 1.d0, A_mo,size(A_mo,1), & mo_coef, size(mo_coef,1), & - 0.d0, T, mo_tot_num_align) + 0.d0, T, size(T,1)) call dgemm('N','N', ao_num, ao_num, mo_tot_num, & 1.d0, mo_coef,size(mo_coef,1), & - T, mo_tot_num_align, & - 0.d0, A_ao, LDA_ao) + T, size(T,1), & + 0.d0, A_ao, size(A_ao,1)) deallocate(T) end @@ -288,18 +282,17 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) double precision, allocatable :: T(:,:) allocate ( T(ao_num_align,ao_num) ) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - call dgemm('T','N', ao_num, ao_num, ao_num, & + call dgemm('T','N', ao_num, ao_num, ao_num, & 1.d0, & - ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1), & - A_ao,LDA_ao, & - 0.d0, T, ao_num_align) + ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),& + A_ao,size(A_ao,1), & + 0.d0, T, size(T,1)) - call dgemm('N','N', ao_num, ao_num, ao_num, 1.d0, & - T, size(T,1), & + call dgemm('N','N', ao_num, ao_num, ao_num, 1.d0, & + T, size(T,1), & ao_ortho_canonical_coef_inv,size(ao_ortho_canonical_coef_inv,1),& - 0.d0, A, LDA) + 0.d0, A, size(A,1)) deallocate(T) end diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f index 750e3420..4221fce4 100644 --- a/src/MO_Basis/utils.irp.f +++ b/src/MO_Basis/utils.irp.f @@ -44,11 +44,12 @@ subroutine save_mos_truncated(n) end -subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign) +subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign,output) implicit none integer,intent(in) :: n,m, sign character*(64), intent(in) :: label double precision, intent(in) :: matrix(n,m) + logical, intent(in) :: output integer :: i,j double precision, allocatable :: mo_coef_new(:,:), R(:,:),eigvalues(:), A(:,:) @@ -76,22 +77,26 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign) mo_coef_new = mo_coef call lapack_diag(eigvalues,R,A,n,m) - write (output_mo_basis,'(A)') 'MOs are now **'//trim(label)//'**' - write (output_mo_basis,'(A)') '' - write (output_mo_basis,'(A)') 'Eigenvalues' - write (output_mo_basis,'(A)') '-----------' - write (output_mo_basis,'(A)') '' - write (output_mo_basis,'(A)') '======== ================' + if (output) then + write (output_mo_basis,'(A)') 'MOs are now **'//trim(label)//'**' + write (output_mo_basis,'(A)') '' + write (output_mo_basis,'(A)') 'Eigenvalues' + write (output_mo_basis,'(A)') '-----------' + write (output_mo_basis,'(A)') '' + write (output_mo_basis,'(A)') '======== ================' + endif if (sign == -1) then do i=1,m eigvalues(i) = -eigvalues(i) enddo endif - do i=1,m - write (output_mo_basis,'(I8,1X,F16.10)') i,eigvalues(i) - enddo - write (output_mo_basis,'(A)') '======== ================' - write (output_mo_basis,'(A)') '' + if (output) then + do i=1,m + write (output_mo_basis,'(I8,1X,F16.10)') i,eigvalues(i) + enddo + write (output_mo_basis,'(A)') '======== ================' + write (output_mo_basis,'(A)') '' + endif call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1)) deallocate(A,mo_coef_new,R,eigvalues) diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 1a3d1d36..7e81756c 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -69,7 +69,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) double precision, allocatable :: U(:,:) double precision, allocatable :: Vt(:,:) double precision, allocatable :: D(:) - double precision, allocatable :: S_half(:,:) + double precision, allocatable :: S(:,:) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j @@ -77,7 +77,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) return endif - allocate (U(ldc,n), Vt(lda,n), D(n), S_half(lda,n)) + allocate (U(ldc,n), Vt(lda,n), D(n), S(lda,n)) call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) @@ -103,13 +103,13 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(S_half,U,D,Vt,n,C,m) & + !$OMP SHARED(S,U,D,Vt,n,C,m) & !$OMP PRIVATE(i,j) !$OMP DO do j=1,n do i=1,n - S_half(i,j) = U(i,j)*D(j) + S(i,j) = U(i,j)*D(j) enddo do i=1,n U(i,j) = C(i,j) @@ -119,8 +119,8 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) !$OMP END PARALLEL - call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1)) - deallocate (U, Vt, D, S_half) + call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1)) + deallocate (U, Vt, D, S) end @@ -210,19 +210,19 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) double precision, allocatable :: U(:,:) double precision, allocatable :: Vt(:,:) double precision, allocatable :: D(:) - double precision, allocatable :: S_half(:,:) + double precision, allocatable :: S(:,:) integer :: info, i, j, k if (n < 2) then return endif - allocate(U(ldc,n),Vt(lda,n),S_half(lda,n),D(n)) + allocate(U(ldc,n),Vt(lda,n),S(lda,n),D(n)) call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(S_half,U,D,Vt,n,C,m) & + !$OMP SHARED(S,U,D,Vt,n,C,m) & !$OMP PRIVATE(i,j,k) !$OMP DO @@ -233,7 +233,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) D(i) = 1.d0/dsqrt(D(i)) endif do j=1,n - S_half(j,i) = 0.d0 + S(j,i) = 0.d0 enddo enddo !$OMP END DO @@ -243,7 +243,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !$OMP DO do j=1,n do i=1,n - S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j) + S(i,j) = S(i,j) + U(i,k)*D(k)*Vt(k,j) enddo enddo !$OMP END DO NOWAIT @@ -261,9 +261,9 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !$OMP END PARALLEL - call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1)) + call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1)) - deallocate(U,Vt,S_half,D) + deallocate(U,Vt,S,D) end diff --git a/src/Utils/constants.include.F b/src/Utils/constants.include.F index 4974fd8e..7f1e71f3 100644 --- a/src/Utils/constants.include.F +++ b/src/Utils/constants.include.F @@ -1,6 +1,6 @@ integer, parameter :: max_dim = 511 integer, parameter :: SIMD_vector = 32 -integer, parameter :: N_int_max = 16 +integer, parameter :: N_int_max = 32 double precision, parameter :: pi = dacos(-1.d0) double precision, parameter :: sqpi = dsqrt(dacos(-1.d0))