Merge pull request #192 from TApplencourt/master

Add stuff useful for PBC support

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;;

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;

plugins/read_integral/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Integrals_Monoelec Integrals_Bielec

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 <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec>`_
+* `Integrals_Bielec <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Bielec>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`print_integrals <http://github.com/LCPQ/quantum_package/tree/master/plugins/read_integral/read_integrals_mo.irp.f#L1>`_
+  Undocumented
+
+
+`run <http://github.com/LCPQ/quantum_package/tree/master/plugins/read_integral/read_integrals_mo.irp.f#L8>`_
+  Undocumented
+

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

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

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}

src/Nuclei/EZFIO.cfg

@@ -20,3 +20,14 @@ doc:  Nuclear coordinates in the format (:, {x,y,z})
 type: double precision 
 size: (nuclei.nucl_num,3)
 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

src/Nuclei/nuclei.irp.f

@@ -146,6 +146,20 @@ END_PROVIDER
    BEGIN_DOC
    ! Nuclear repulsion energy
    END_DOC
+
+   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
@@ -163,10 +177,15 @@ END_PROVIDER
          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)]