Add Read / Write for Nuclear Repulsion (usefull for PBC)

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;

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 ->
+         | Nuclei_by_hand ->
-           f Nuclei.(read, to_rst)
+           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

@@ -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
+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

@@ -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)
+   IF (disk_access_nuclear_repulsion.EQ.'Read') THEN
-   nuclear_repulsion = 0.d0
+          print*, 'nuclear_repulsion read from disk'
-   do l = 1, nucl_num
+          LOGICAL :: has
-     do  k = 1, nucl_num
+          call ezfio_has_nuclei_nuclear_repulsion(has)
-       if(k == l) then
+          if (has) then
-         cycle
+                 call ezfio_get_nuclei_nuclear_repulsion(nuclear_repulsion)
-       endif
+          else
-       Z12 = nucl_charge(k)*nucl_charge(l)
+                 print *, 'nuclei/nuclear_repulsion not found in EZFIO file'
-       x(1) = nucl_coord(k,1) - nucl_coord(l,1)
+                 stop 1
-       x(2) = nucl_coord(k,2) - nucl_coord(l,2)
+          endif
-       x(3) = nucl_coord(k,3) - nucl_coord(l,3)
+
-       r2 = x(1)*x(1) + x(2)*x(2) + x(3)*x(3)
+   ELSE
-       nuclear_repulsion += Z12/dsqrt(r2)
+
-     enddo
+      integer                        :: k,l
-   enddo
+      double precision               :: Z12, r2, x(3)
-   nuclear_repulsion *= 0.5d0
+      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)