diff --git a/ocaml/Basis.ml b/ocaml/Basis.ml index 869fb132..797d53f2 100644 --- a/ocaml/Basis.ml +++ b/ocaml/Basis.ml @@ -36,9 +36,11 @@ let read_element in_channel at_number element = -let to_string_general ~fmt ~atom_sep b = +let to_string_general ~fmt ~atom_sep ?ele_array b = let new_nucleus n = - Printf.sprintf "Atom %d" n + match ele_array with + | None -> Printf.sprintf "Atom %d" n + | Some x -> Printf.sprintf "%s" (Element.to_string x.(n-1)) in let rec do_work accu current_nucleus = function | [] -> List.rev accu @@ -56,12 +58,12 @@ let to_string_general ~fmt ~atom_sep b = do_work [new_nucleus 1] 1 b |> String.concat ~sep:"\n" -let to_string_gamess = - to_string_general ~fmt:Gto.Gamess ~atom_sep:"" +let to_string_gamess ?ele_array = + to_string_general ?ele_array ~fmt:Gto.Gamess ~atom_sep:"" -let to_string_gaussian b = +let to_string_gaussian ?ele_array b = String.concat ~sep:"\n" - [ to_string_general ~fmt:Gto.Gaussian ~atom_sep:"****" b ; "****" ] + [ to_string_general ?ele_array ~fmt:Gto.Gaussian ~atom_sep:"****" b ; "****" ] let to_string ?(fmt=Gto.Gamess) = match fmt with diff --git a/ocaml/Basis.mli b/ocaml/Basis.mli index 249c14f9..41ddc184 100644 --- a/ocaml/Basis.mli +++ b/ocaml/Basis.mli @@ -14,7 +14,7 @@ val read_element : in_channel -> Nucl_number.t -> Element.t -> (Gto.t * Nucl_number.t) list (** Convert the basis to a string *) -val to_string : ?fmt:Gto.fmt -> (Gto.t * Nucl_number.t) list -> string +val to_string : ?fmt:Gto.fmt -> ?ele_array:Element.t array -> (Gto.t * Nucl_number.t) list -> string (** Convert the basis to an MD5 hash *) val to_md5 : (Gto.t * Nucl_number.t) list -> MD5.t diff --git a/plugins/analyze_wf/NEEDED_CHILDREN_MODULES b/plugins/analyze_wf/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..aae89501 --- /dev/null +++ b/plugins/analyze_wf/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Determinants diff --git a/plugins/analyze_wf/README.rst b/plugins/analyze_wf/README.rst new file mode 100644 index 00000000..179e407d --- /dev/null +++ b/plugins/analyze_wf/README.rst @@ -0,0 +1,12 @@ +========== +analyze_wf +========== + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/analyze_wf/analyze_wf.irp.f b/plugins/analyze_wf/analyze_wf.irp.f new file mode 100644 index 00000000..6d8bffcf --- /dev/null +++ b/plugins/analyze_wf/analyze_wf.irp.f @@ -0,0 +1,70 @@ +program analyze_wf + implicit none + BEGIN_DOC +! Wave function analyzis + END_DOC + read_wf = .True. + SOFT_TOUCH read_wf + call run() +end + +subroutine run + implicit none + integer :: istate, i + integer :: class(0:mo_tot_num,5) + double precision :: occupation(mo_tot_num) + + write(*,'(A)') 'MO Occupation' + write(*,'(A)') '=============' + write(*,'(A)') '' + do istate=1,N_states + call get_occupation_from_dets(occupation,1) + write(*,'(A)') '' + write(*,'(A,I3)'), 'State ', istate + write(*,'(A)') '---------------' + write(*,'(A)') '' + write (*,'(A)') '======== ================' + class = 0 + do i=1,mo_tot_num + write (*,'(I8,X,F16.10)') i, occupation(i) + if (occupation(i) > 1.999d0) then + class(0,1) += 1 + class( class(0,1), 1) = i + else if (occupation(i) > 1.95d0) then + class(0,2) += 1 + class( class(0,2), 2) = i + else if (occupation(i) < 0.001d0) then + class(0,5) += 1 + class( class(0,5), 5) = i + else if (occupation(i) < 0.01d0) then + class(0,4) += 1 + class( class(0,4), 4) = i + else + class(0,3) += 1 + class( class(0,3), 3) = i + endif + enddo + write (*,'(A)') '======== ================' + write (*,'(A)') '' + + write (*,'(A)') 'Suggested classes' + write (*,'(A)') '-----------------' + write (*,'(A)') '' + write (*,'(A)') 'Core :' + write (*,*) (class(i,1), ',', i=1,class(0,1)) + write (*,*) '' + write (*,'(A)') 'Inactive :' + write (*,*) (class(i,2), ',', i=1,class(0,2)) + write (*,'(A)') '' + write (*,'(A)') 'Active :' + write (*,*) (class(i,3), ',', i=1,class(0,3)) + write (*,'(A)') '' + write (*,'(A)') 'Virtual :' + write (*,*) (class(i,4), ',', i=1,class(0,4)) + write (*,'(A)') '' + write (*,'(A)') 'Deleted :' + write (*,*) (class(i,5), ',', i=1,class(0,5)) + write (*,'(A)') '' + enddo + +end diff --git a/plugins/analyze_wf/occupation.irp.f b/plugins/analyze_wf/occupation.irp.f new file mode 100644 index 00000000..d426dc14 --- /dev/null +++ b/plugins/analyze_wf/occupation.irp.f @@ -0,0 +1,23 @@ +subroutine get_occupation_from_dets(occupation, istate) + implicit none + double precision, intent(out) :: occupation(mo_tot_num) + integer, intent(in) :: istate + BEGIN_DOC + ! Returns the average occupation of the MOs + END_DOC + integer :: i,j, ispin + integer :: list(N_int*bit_kind_size,2) + integer :: n_elements(2) + double precision :: c + + occupation = 0.d0 + do i=1,N_det + c = psi_coef(i,istate)*psi_coef(i,istate) + call bitstring_to_list_ab(psi_det(1,1,i), list, n_elements, N_int) + do ispin=1,2 + do j=1,n_elements(ispin) + occupation( list(j,ispin) ) += c + enddo + enddo + enddo +end diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index 118bbdf7..ed2f49bd 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -194,7 +194,6 @@ subroutine set_natural_mos double precision, allocatable :: tmp(:,:) label = "Natural" -! call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,label,-1) call mo_as_svd_vectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,mo_tot_num,label) end