From 414aa324633019d0b12c710bc7a7098aa9739d70 Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Mon, 8 May 2017 16:27:29 +0000 Subject: [PATCH 1/2] Add Read / Write for Nuclear Repulsion (usefull for PBC) --- ocaml/Input.ml | 2 +- ...nput_nuclei.ml => Input_nuclei_by_hand.ml} | 2 +- scripts/ezfio_interface/qp_edit_template | 12 ++-- src/Nuclei/EZFIO.cfg | 13 +++- src/Nuclei/nuclei.irp.f | 59 ++++++++++++------- 5 files changed, 59 insertions(+), 29 deletions(-) rename ocaml/{Input_nuclei.ml => Input_nuclei_by_hand.ml} (99%) diff --git a/ocaml/Input.ml b/ocaml/Input.ml index 5e012996..b1bb1060 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -7,5 +7,5 @@ include Input_bitmasks;; include Input_determinants_by_hand;; include Input_electrons;; include Input_mo_basis;; -include Input_nuclei;; +include Input_nuclei_by_hand;; include Input_auto_generated;; diff --git a/ocaml/Input_nuclei.ml b/ocaml/Input_nuclei_by_hand.ml similarity index 99% rename from ocaml/Input_nuclei.ml rename to ocaml/Input_nuclei_by_hand.ml index ca81629e..f36b6b82 100644 --- a/ocaml/Input_nuclei.ml +++ b/ocaml/Input_nuclei_by_hand.ml @@ -2,7 +2,7 @@ open Qptypes;; open Qputils;; open Core.Std;; -module Nuclei : sig +module Nuclei_by_hand : sig type t = { nucl_num : Nucl_number.t ; nucl_label : Element.t array; diff --git a/scripts/ezfio_interface/qp_edit_template b/scripts/ezfio_interface/qp_edit_template index af9b295c..9d327124 100644 --- a/scripts/ezfio_interface/qp_edit_template +++ b/scripts/ezfio_interface/qp_edit_template @@ -20,7 +20,7 @@ type keyword = | Determinants_by_hand | Electrons | Mo_basis -| Nuclei +| Nuclei_by_hand {keywords} @@ -30,7 +30,7 @@ let keyword_to_string = function | Determinants_by_hand -> "Determinants_by_hand" | Electrons -> "Electrons" | Mo_basis -> "MO basis" -| Nuclei -> "Molecule" +| Nuclei_by_hand -> "Molecule" {keywords_to_string} @@ -74,8 +74,8 @@ let get s = f Mo_basis.(read, to_rst) | Electrons -> f Electrons.(read, to_rst) - | Nuclei -> - f Nuclei.(read, to_rst) + | Nuclei_by_hand -> + f Nuclei_by_hand.(read, to_rst) | Ao_basis -> f Ao_basis.(read, to_rst) | Determinants_by_hand -> @@ -121,7 +121,7 @@ let set str s = {write} | Electrons -> write Electrons.(of_rst, write) s | Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s - | Nuclei -> write Nuclei.(of_rst, write) s + | Nuclei_by_hand -> write Nuclei_by_hand.(of_rst, write) s | Ao_basis -> () (* TODO *) | Mo_basis -> () (* TODO *) end @@ -184,7 +184,7 @@ let run check_only ?ndet ?state ezfio_filename = *) let tasks = [ - Nuclei ; + Nuclei_by_hand ; Ao_basis; Electrons ; {tasks} diff --git a/src/Nuclei/EZFIO.cfg b/src/Nuclei/EZFIO.cfg index 81c168ff..dab6b211 100644 --- a/src/Nuclei/EZFIO.cfg +++ b/src/Nuclei/EZFIO.cfg @@ -19,4 +19,15 @@ interface: ezfio, provider doc: Nuclear coordinates in the format (:, {x,y,z}) type: double precision size: (nuclei.nucl_num,3) -interface: ezfio \ No newline at end of file +interface: ezfio + +[disk_access_nuclear_repulsion] +doc: Read/Write Nuclear Repulsion from/to disk [ Write | Read | None ] +type: Disk_access +interface: ezfio,provider,ocaml +default: None + +[nuclear_repulsion] +doc: Nuclear repulsion (Computed automaticaly or Read in the EZFIO) +type:double precision +interface: ezfio diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f index 34fae989..577b8b92 100644 --- a/src/Nuclei/nuclei.irp.f +++ b/src/Nuclei/nuclei.irp.f @@ -141,34 +141,53 @@ BEGIN_PROVIDER [ double precision, positive_charge_barycentre,(3)] enddo END_PROVIDER - BEGIN_PROVIDER [ double precision, nuclear_repulsion ] +BEGIN_PROVIDER [ double precision, nuclear_repulsion ] implicit none BEGIN_DOC ! Nuclear repulsion energy END_DOC - integer :: k,l - double precision :: Z12, r2, x(3) - nuclear_repulsion = 0.d0 - do l = 1, nucl_num - do k = 1, nucl_num - if(k == l) then - cycle - endif - Z12 = nucl_charge(k)*nucl_charge(l) - x(1) = nucl_coord(k,1) - nucl_coord(l,1) - x(2) = nucl_coord(k,2) - nucl_coord(l,2) - x(3) = nucl_coord(k,3) - nucl_coord(l,3) - r2 = x(1)*x(1) + x(2)*x(2) + x(3)*x(3) - nuclear_repulsion += Z12/dsqrt(r2) - enddo - enddo - nuclear_repulsion *= 0.5d0 - + + IF (disk_access_nuclear_repulsion.EQ.'Read') THEN + print*, 'nuclear_repulsion read from disk' + LOGICAL :: has + call ezfio_has_nuclei_nuclear_repulsion(has) + if (has) then + call ezfio_get_nuclei_nuclear_repulsion(nuclear_repulsion) + else + print *, 'nuclei/nuclear_repulsion not found in EZFIO file' + stop 1 + endif + + ELSE + + integer :: k,l + double precision :: Z12, r2, x(3) + nuclear_repulsion = 0.d0 + do l = 1, nucl_num + do k = 1, nucl_num + if(k == l) then + cycle + endif + Z12 = nucl_charge(k)*nucl_charge(l) + x(1) = nucl_coord(k,1) - nucl_coord(l,1) + x(2) = nucl_coord(k,2) - nucl_coord(l,2) + x(3) = nucl_coord(k,3) - nucl_coord(l,3) + r2 = x(1)*x(1) + x(2)*x(2) + x(3)*x(3) + nuclear_repulsion += Z12/dsqrt(r2) + enddo + enddo + nuclear_repulsion *= 0.5d0 + END IF + call write_time(output_Nuclei) call write_double(output_Nuclei,nuclear_repulsion, & 'Nuclear repulsion energy') + + IF (disk_access_nuclear_repulsion.EQ.'Write') THEN + call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion) + END IF END_PROVIDER - + BEGIN_PROVIDER [ character*(128), element_name, (78)] BEGIN_DOC ! Array of the name of element, sorted by nuclear charge (integer) From 7a455d8fd611c7a24fb3b2ac315812dda4b6524e Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Mon, 8 May 2017 16:38:51 +0000 Subject: [PATCH 2/2] Add plugings to read the Molecular integrals --- plugins/read_integral/NEEDED_CHILDREN_MODULES | 1 + plugins/read_integral/README.rst | 30 +++++++++ .../read_integral/print_integrals_mo.irp.f | 61 +++++++++++++++++ plugins/read_integral/read_integrals_mo.irp.f | 66 +++++++++++++++++++ plugins/read_integral/tree_dependency.png | 0 5 files changed, 158 insertions(+) create mode 100644 plugins/read_integral/NEEDED_CHILDREN_MODULES create mode 100644 plugins/read_integral/README.rst create mode 100644 plugins/read_integral/print_integrals_mo.irp.f create mode 100644 plugins/read_integral/read_integrals_mo.irp.f create mode 100644 plugins/read_integral/tree_dependency.png diff --git a/plugins/read_integral/NEEDED_CHILDREN_MODULES b/plugins/read_integral/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..e492a3ce --- /dev/null +++ b/plugins/read_integral/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Integrals_Monoelec Integrals_Bielec diff --git a/plugins/read_integral/README.rst b/plugins/read_integral/README.rst new file mode 100644 index 00000000..02b63512 --- /dev/null +++ b/plugins/read_integral/README.rst @@ -0,0 +1,30 @@ +============= +read_integral +============= + +Warning: CAN NOT CHANGE THE NUMBER OF MO ! + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +.. image:: tree_dependency.png + +* `Integrals_Monoelec `_ +* `Integrals_Bielec `_ + +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +`print_integrals `_ + Undocumented + + +`run `_ + Undocumented + diff --git a/plugins/read_integral/print_integrals_mo.irp.f b/plugins/read_integral/print_integrals_mo.irp.f new file mode 100644 index 00000000..133e34b8 --- /dev/null +++ b/plugins/read_integral/print_integrals_mo.irp.f @@ -0,0 +1,61 @@ +program print_integrals + implicit none + + integer :: iunit + integer :: getunitandopen + + integer ::i,j,k,l + double precision :: integral + + iunit = getunitandopen('kinetic_mo','w') + do i=1,mo_tot_num + do j=1,mo_tot_num + write(iunit,*) i,j, mo_kinetic_integral(i,j) + enddo + enddo + close(iunit) + + iunit = getunitandopen('overlap_mo','w') + do i=1,mo_tot_num + do j=1,mo_tot_num + write(iunit,*) i,j, mo_overlap(i,j) + enddo + enddo + close(iunit) + + iunit = getunitandopen('nuclear_mo','w') + do i=1,mo_tot_num + do j=1,mo_tot_num + write(iunit,*) i,j, mo_nucl_elec_integral(i,j) + enddo + enddo + close(iunit) + + !iunit = getunitandopen('pseudo_mo','w') + !do i=1,mo_tot_num + ! do j=1,mo_tot_num + ! write(iunit,*) i,j, mo_pseudo_integral(i,j) + ! enddo + !enddo + !close(iunit) + + PROVIDE mo_bielec_integrals_in_map + iunit = getunitandopen('bielec_mo','w') + do l=1,mo_tot_num + do k=1,mo_tot_num + do j=l,mo_tot_num + do i=k,mo_tot_num + !if (i>=j) then + double precision :: get_mo_bielec_integral + integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + if (dabs(integral) > mo_integrals_threshold) then + write (iunit,'(4(I5,X),D22.15)') i,j,k,l, integral + endif + !end if + enddo + enddo + enddo + enddo + + close(iunit) +end diff --git a/plugins/read_integral/read_integrals_mo.irp.f b/plugins/read_integral/read_integrals_mo.irp.f new file mode 100644 index 00000000..e9868e56 --- /dev/null +++ b/plugins/read_integral/read_integrals_mo.irp.f @@ -0,0 +1,66 @@ +program print_integrals + call run +end + +subroutine run + use map_module + implicit none + + integer :: iunit + integer :: getunitandopen + + integer ::i,j,k,l + double precision :: integral + double precision, allocatable :: A(:,:) + + integer :: n_integrals + integer(key_kind), allocatable :: buffer_i(:) + real(integral_kind), allocatable :: buffer_values(:) + integer(key_kind) :: key + + call ezfio_set_mo_basis_mo_tot_num(mo_tot_num) + + allocate (A(mo_tot_num_align,mo_tot_num)) + + iunit = getunitandopen('kinetic_mo','r') + do + read (iunit,*,end=10) i,j, integral + A(i,j) = integral + enddo + 10 continue + close(iunit) + call write_one_e_integrals('mo_kinetic_integral', A, size(A,1), size(A,2)) + + + iunit = getunitandopen('nuclear_mo','r') + do + read (iunit,*,end=12) i,j, integral + A(i,j) = integral + enddo + 12 continue + close(iunit) + call write_one_e_integrals('mo_ne_integral', A, size(A,1), size(A,2)) + + call ezfio_set_integrals_monoelec_disk_access_mo_one_integrals("Read") + + allocate(buffer_i(mo_tot_num**4), buffer_values(mo_tot_num**4)) + + iunit = getunitandopen('bielec_mo','r') + n_integrals=0 + do + read (iunit,*,end=13) i,j,k,l, integral + n_integrals += 1 + call bielec_integrals_index(i, j, k, l, buffer_i(n_integrals) ) + buffer_values(n_integrals) = integral + enddo + 13 continue + close(iunit) + + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_values) + + call map_sort(mo_integrals_map) + + call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) + call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read") + +end diff --git a/plugins/read_integral/tree_dependency.png b/plugins/read_integral/tree_dependency.png new file mode 100644 index 00000000..e69de29b