From 7ebc2ac89620a8cad0254b02673e92790d3a8d70 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Dec 2017 13:48:21 +0100 Subject: [PATCH] Fixed last merge (#217) * Fixes #211 * Removed debug print * Revert input coordinates * Fix MPI --- data/list_element.txt | 174 +++++----------- ocaml/Element.ml | 1 + ocaml/Makefile | 5 +- ocaml/element_create_db.ml | 27 +++ plugins/QMC/save_for_qmcchem.irp.f | 1 - plugins/Symmetry/Symmetry.main.irp.f | 14 +- plugins/Symmetry/aos.irp.f | 38 ++-- plugins/Symmetry/find_sym.irp.f | 21 +- plugins/Symmetry/nuclei.irp.f | 92 +++++++++ scripts/pseudo/elts_num_ele.py | 134 ++---------- src/Nuclei/inertia.irp.f | 14 +- src/Nuclei/nuclei.irp.f | 299 ++++++--------------------- 12 files changed, 307 insertions(+), 513 deletions(-) create mode 100644 ocaml/element_create_db.ml create mode 100644 plugins/Symmetry/nuclei.irp.f diff --git a/data/list_element.txt b/data/list_element.txt index b2d081c1..2a9733e6 100644 --- a/data/list_element.txt +++ b/data/list_element.txt @@ -1,118 +1,56 @@ -1 H Hydrogen -2 He Helium -3 Li Lithium -4 Be Beryllium -5 B Boron -6 C Carbon -7 N Nitrogen -8 O Oxygen -9 F Fluorine -10 Ne Neon -11 Na Sodium -12 Mg Magnesium -13 Al Aluminum -14 Si Silicon -15 P Phosphorus -16 S Sulfur -17 Cl Chlorine -18 Ar Argon -19 K Potassium -20 Ca Calcium -21 Sc Scandium -22 Ti Titanium -23 V Vanadium -24 Cr Chromium -25 Mn Manganese -26 Fe Iron -27 Co Cobalt -28 Ni Nickel -29 Cu Copper -30 Zn Zinc -31 Ga Gallium -32 Ge Germanium -33 As Arsenic -34 Se Selenium -35 Br Bromine -36 Kr Krypton -37 Rb Rubidium -38 Sr Strontium -39 Y Yttrium -40 Zr Zirconium -41 Nb Niobium -42 Mo Molybdenum -43 Tc Technetium -44 Ru Ruthenium -45 Rh Rhodium -46 Pd Palladium -47 Ag Silver -48 Cd Cadmium -49 In Indium -50 Sn Tin -51 Sb Antimony -52 Te Tellurium -53 I Iodine -54 Xe Xenon -55 Cs Cesium -56 Ba Barium -57 La Lanthanum -58 Ce Cerium -59 Pr Praseodymium -60 Nd Neodymium -61 Pm Promethium -62 Sm Samarium -63 Eu Europium -64 Gd Gadolinium -65 Tb Terbium -66 Dy Dysprosium -67 Ho Holmium -68 Er Erbium -69 Tm Thulium -70 Yb Ytterbium -71 Lu Lutetium -72 Hf Hafnium -73 Ta Tantalum -74 W Tungsten -75 Re Rhenium -76 Os Osmium -77 Ir Iridium -78 Pt Platinum -79 Au Gold -80 Hg Mercury -81 Tl Thallium -82 Pb Lead -83 Bi Bismuth -84 Po Polonium -85 At Astatine -86 Rn Radon -87 Fr Francium -88 Ra Radium -89 Ac Actinium -90 Th Thorium -91 Pa Protactinium -92 U Uranium -93 Np Neptunium -94 Pu Plutonium -95 Am Americium -96 Cm Curium -97 Bk Berkelium -98 Cf Californium -99 Es Einsteinium -100 Fm Fermium -101 Md Mendelevium -102 No Nobelium -103 Lr Lawrencium -104 Rf Rutherfordium -105 Db Dubnium -106 Sg Seaborgium -107 Bh Bohrium -108 Hs Hassium -109 Mt Meitnerium -110 Ds Darmstadtium -111 Rg Roentgenium -112 Cn Copernicium -113 Uut Ununtrium -114 Fl Flerovium -115 Uup Ununpentium -116 Lv Livermorium -117 Uus Ununseptium -118 Uuo Ununoctium \ No newline at end of file + 0 X Dummy 0.000000 + 1 H Hydrogen 1.007900 + 2 He Helium 4.002600 + 3 Li Lithium 6.941000 + 4 Be Beryllium 9.012180 + 5 B Boron 10.810000 + 6 C Carbon 12.011000 + 7 N Nitrogen 14.006700 + 8 O Oxygen 15.999400 + 9 F Fluorine 18.998403 + 10 Ne Neon 20.179000 + 11 Na Sodium 22.989770 + 12 Mg Magnesium 24.305000 + 13 Al Aluminum 26.981540 + 14 Si Silicon 28.085500 + 15 P Phosphorus 30.973760 + 16 S Sulfur 32.060000 + 17 Cl Chlorine 35.453000 + 18 Ar Argon 39.948000 + 19 K Potassium 39.098300 + 20 Ca Calcium 40.080000 + 21 Sc Scandium 44.955900 + 22 Ti Titanium 47.900000 + 23 V Vanadium 50.941500 + 24 Cr Chromium 51.996000 + 25 Mn Manganese 54.938000 + 26 Fe Iron 55.933200 + 27 Co Cobalt 58.933200 + 28 Ni Nickel 58.700000 + 29 Cu Copper 63.546000 + 30 Zn Zinc 65.380000 + 31 Ga Gallium 69.720000 + 32 Ge Germanium 72.590000 + 33 As Arsenic 74.921600 + 34 Se Selenium 78.960000 + 35 Br Bromine 79.904000 + 36 Kr Krypton 83.800000 + 37 Rb Rubidium 85.467800 + 38 Sr Strontium 87.620000 + 39 Y Yttrium 88.905840 + 40 Zr Zirconium 91.224000 + 41 Nb Niobium 92.906370 + 42 Mo Molybdenum 95.950000 + 43 Tc Technetium 98.000000 + 44 Ru Ruthenium 101.070000 + 45 Rh Rhodium 102.905500 + 46 Pd Palladium 106.420000 + 47 Ag Silver 107.868200 + 48 Cd Cadmium 112.414000 + 49 In Indium 114.818000 + 50 Sn Tin 118.710000 + 51 Sb Antimony 121.760000 + 52 Te Tellurium 127.600000 + 53 I Iodine 126.904470 + 54 Xe Xenon 131.293000 + 78 Pt Platinum 195.084000 diff --git a/ocaml/Element.ml b/ocaml/Element.ml index fd08b8da..d3b68d52 100644 --- a/ocaml/Element.ml +++ b/ocaml/Element.ml @@ -499,3 +499,4 @@ let mass x = result x |> Positive_float.of_float + diff --git a/ocaml/Makefile b/ocaml/Makefile index 3534c614..b666187f 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -17,7 +17,7 @@ MLLFILES=$(wildcard *.mll) MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml Input_auto_generated.ml qp_edit.ml MLIFILES=$(wildcard *.mli) git ALL_TESTS=$(patsubst %.ml,%.byte,$(wildcard test_*.ml)) -ALL_EXE=$(patsubst %.ml,%.native,$(wildcard qp_*.ml)) qp_edit.native +ALL_EXE=$(patsubst %.ml,%.native,$(wildcard qp_*.ml)) qp_edit.native element_create_db.byte .PHONY: executables default remake_executables @@ -37,8 +37,9 @@ tests: $(ALL_TESTS) executables: $(QP_ROOT)/data/executables -$(QP_ROOT)/data/executables: remake_executables +$(QP_ROOT)/data/executables: remake_executables element_create_db.byte Qptypes.ml $(QP_ROOT)/scripts/module/create_executables_list.sh + $(QP_ROOT)/ocaml/element_create_db.byte external_libs: opam install cryptokit core diff --git a/ocaml/element_create_db.ml b/ocaml/element_create_db.ml new file mode 100644 index 00000000..03fbf3cd --- /dev/null +++ b/ocaml/element_create_db.ml @@ -0,0 +1,27 @@ +open Core +open Qptypes +open Element + +let () = + let indices = + Array.init 78 (fun i -> i) + in + Out_channel.with_file (Qpackage.root ^ "/data/list_element.txt") + ~f:(fun out_channel -> + Array.init 110 ~f:(fun i -> + let element = + try + Some (of_charge (Charge.of_int i)) + with + | _ -> None + in + match element with + | None -> "" + | Some x -> Printf.sprintf "%3d %3s %s %f\n" + i (to_string x) (to_long_string x) (Positive_float.to_float @@ mass x ) + ) + |> Array.to_list + |> String.concat ~sep:"" + |> Out_channel.output_string out_channel + ) + diff --git a/plugins/QMC/save_for_qmcchem.irp.f b/plugins/QMC/save_for_qmcchem.irp.f index 771bf618..e9fa60c4 100644 --- a/plugins/QMC/save_for_qmcchem.irp.f +++ b/plugins/QMC/save_for_qmcchem.irp.f @@ -1,7 +1,6 @@ program save_for_qmc integer :: iunit - integer, external :: get_unit_and_open logical :: exists double precision :: e_ref diff --git a/plugins/Symmetry/Symmetry.main.irp.f b/plugins/Symmetry/Symmetry.main.irp.f index 78bf126e..1d438af5 100644 --- a/plugins/Symmetry/Symmetry.main.irp.f +++ b/plugins/Symmetry/Symmetry.main.irp.f @@ -6,13 +6,13 @@ program Symmetry integer :: i, j character*8 :: sym - print *, 'Molecule is linear:', molecule_is_linear - print *, 'Has center of inversion:', molecule_has_center_of_inversion - print *, 'Has S2n improper rotation:', molecule_has_improper_rotation - print *, 'Symmetry rotation axis:', sym_rotation_axis(:) - print *, 'Group:'//point_group - print *, 'Symmetry irreps', sym_irrep(1:n_irrep) - print *, 'Symmetry operations', sym_operation(1:n_irrep) + print *, 'Molecule is linear: ', molecule_is_linear + print *, 'Has center of inversion: ', molecule_has_center_of_inversion + print *, 'Has S2n improper rotation: ', molecule_has_improper_rotation + print *, 'Symmetry rotation axis: ', sym_rotation_axis(:) + print *, 'Group: '//point_group + print *, 'Symmetry irreps : ', sym_irrep(1:n_irrep) + print *, 'Symmetry operations : ', sym_operation(1:n_irrep) print *, 'Character table' do i=1,n_irrep print *, character_table(i,:) diff --git a/plugins/Symmetry/aos.irp.f b/plugins/Symmetry/aos.irp.f index e4352948..e0aade1e 100644 --- a/plugins/Symmetry/aos.irp.f +++ b/plugins/Symmetry/aos.irp.f @@ -7,8 +7,8 @@ BEGIN_PROVIDER [ double precision, sym_box, (3,2) ] sym_box(:,:) = 0.d0 do xyz=1,3 do i=1,nucl_num - sym_box(xyz,1) = min(sym_box(xyz,1), nucl_coord(i,xyz)) - sym_box(xyz,2) = max(sym_box(xyz,2), nucl_coord(i,xyz)) + sym_box(xyz,1) = min(sym_box(xyz,1), nucl_coord_sym(i,xyz)) + sym_box(xyz,2) = max(sym_box(xyz,2), nucl_coord_sym(i,xyz)) enddo enddo sym_box(:,1) = sym_box(:,1) - 2.d0 @@ -42,24 +42,28 @@ subroutine compute_sym_ao_values(sym_points, n_sym_points, result) double precision, intent(out) :: result(n_sym_points, ao_num) integer :: i, j double precision :: x, y, z + double precision :: x2, y2, z2 result (:,:) = 0.d0 do j=1,ao_num do i=1,n_sym_points - x = sym_points(1,i) - nucl_coord_transp(1,ao_nucl(j)) - y = sym_points(2,i) - nucl_coord_transp(2,ao_nucl(j)) - z = sym_points(3,i) - nucl_coord_transp(3,ao_nucl(j)) - x = x**ao_power(j,1) - y = y**ao_power(j,2) - z = z**ao_power(j,3) -! result(i,j) = x*y*z*exp(-(x*x+y*y+z*z)) - result(i,j) = x*y*z - if (result(i,j) > 0.d0) then - result(i,j) = 1.d0 - else if (result(i,j) < 0.d0) then - result(i,j) = -1.d0 - else - result(i,j) = 0.d0 - endif + x = sym_points(1,i) - nucl_coord_sym_transp(1,ao_nucl(j)) + y = sym_points(2,i) - nucl_coord_sym_transp(2,ao_nucl(j)) + z = sym_points(3,i) - nucl_coord_sym_transp(3,ao_nucl(j)) + x2 = x*sym_molecule_rotation_inv(1,1) + y*sym_molecule_rotation_inv(2,1) + z*sym_molecule_rotation_inv(3,1) + y2 = x*sym_molecule_rotation_inv(1,2) + y*sym_molecule_rotation_inv(2,2) + z*sym_molecule_rotation_inv(3,2) + z2 = x*sym_molecule_rotation_inv(1,3) + y*sym_molecule_rotation_inv(2,3) + z*sym_molecule_rotation_inv(3,3) + x = x2**ao_power(j,1) + y = y2**ao_power(j,2) + z = z2**ao_power(j,3) + result(i,j) = x*y*z*exp(-(x*x+y*y+z*z)) +! result(i,j) = x*y*z +! if (result(i,j) > 0.d0) then +! result(i,j) = 1.d0 +! else if (result(i,j) < 0.d0) then +! result(i,j) = -1.d0 +! else +! result(i,j) = 0.d0 +! endif enddo enddo diff --git a/plugins/Symmetry/find_sym.irp.f b/plugins/Symmetry/find_sym.irp.f index 80e2fccd..95bb4e6c 100644 --- a/plugins/Symmetry/find_sym.irp.f +++ b/plugins/Symmetry/find_sym.irp.f @@ -20,7 +20,7 @@ BEGIN_PROVIDER [ logical, molecule_has_center_of_inversion ] found = .False. do j=1,nucl_num if (nucl_charge(i) /= nucl_charge(j)) cycle - point(:) = nucl_coord_transp(:,i) + nucl_coord_transp(:,j) + point(:) = nucl_coord_sym_transp(:,i) + nucl_coord_sym_transp(:,j) if (u_dot_u(point,3) < 1.d-5) then found = .True. exit @@ -52,10 +52,10 @@ BEGIN_PROVIDER [ integer, sym_rotation_axis, (3) ] sym_rotation_axis(iaxis) = iorder do i=1,nucl_num found = .False. - call sym_apply_rotation(dble(iorder),iaxis,nucl_coord_transp(1,i),point) + call sym_apply_rotation(dble(iorder),iaxis,nucl_coord_sym_transp(1,i),point) do j=1,nucl_num if (nucl_charge(i) /= nucl_charge(j)) cycle - point2(:) = nucl_coord_transp(:,j) - point(:) + point2(:) = nucl_coord_sym_transp(:,j) - point(:) if (u_dot_u(point2,3) < 1.d-5) then found = .True. exit @@ -148,10 +148,10 @@ BEGIN_PROVIDER [ logical, molecule_has_improper_rotation ] molecule_has_improper_rotation = .True. do i=1,nucl_num found = .False. - call sym_apply_improper_rotation(dble(iorder),iaxis,nucl_coord_transp(1,i),point) + call sym_apply_improper_rotation(dble(iorder),iaxis,nucl_coord_sym_transp(1,i),point) do j=1,nucl_num if (nucl_charge(i) /= nucl_charge(j)) cycle - point2(:) = nucl_coord_transp(:,j) - point(:) + point2(:) = nucl_coord_sym_transp(:,j) - point(:) if (u_dot_u(point2,3) < 1.d-5) then found = .True. exit @@ -179,7 +179,7 @@ BEGIN_PROVIDER [ logical, molecule_has_center_of_inversion ] found = .False. do j=1,nucl_num if (nucl_charge(i) /= nucl_charge(j)) cycle - point(:) = nucl_coord_transp(:,i) + nucl_coord_transp(:,j) + point(:) = nucl_coord_sym_transp(:,i) + nucl_coord_sym_transp(:,j) if (u_dot_u(point,3) < 1.d-5) then found = .True. exit @@ -208,11 +208,11 @@ BEGIN_PROVIDER [ logical, molecule_has_sigma_plane, (3) ] molecule_has_sigma_plane(iaxis) = .True. do i=1,nucl_num found = .False. - point(:) = nucl_coord_transp(:,i) + point(:) = nucl_coord_sym_transp(:,i) point(iaxis) = -point(iaxis) do j=1,nucl_num if (nucl_charge(i) /= nucl_charge(j)) cycle - point2(:) = nucl_coord_transp(:,j) - point(:) + point2(:) = nucl_coord_sym_transp(:,j) - point(:) if (u_dot_u(point2,3) < 1.d-5) then found = .True. exit @@ -233,7 +233,7 @@ BEGIN_PROVIDER [ character*16, point_group ] ! Point group of the molecule END_DOC - character*2, save :: i_to_a(24) = (/ ' 1', ' 2', ' 3', ' 4', ' 5', ' 6', ' 7', ' 8', ' 9', & + character*2, save :: i_to_a(24) = (/ '1 ', '2 ', '3 ', '4 ', '5 ', '6 ', '7 ', '8 ', '9 ', & '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '20', & '21', '22', '23', '24' /) point_group = 'C1' @@ -366,7 +366,7 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] integer :: iangle, n_sym_points double precision :: angle integer :: iop, imo, ipoint, l, i - double precision :: sym_operations_on_mos(n_irrep) + double precision :: sym_operations_on_mos(mo_tot_num) logical :: possible_irrep(n_irrep,mo_tot_num) n_sym_points = 10000 @@ -443,6 +443,7 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] sym_operations_on_mos(imo) += x enddo sym_operations_on_mos(imo) *= 1.d0/n_sym_points + print *, iop, imo, sym_operations_on_mos(imo) if (dabs(sym_operations_on_mos(imo)-1.d0) < 1.d-2) then sym_operations_on_mos(imo)=1.d0 else if (dabs(sym_operations_on_mos(imo)+1.d0) < 1.d-2) then diff --git a/plugins/Symmetry/nuclei.irp.f b/plugins/Symmetry/nuclei.irp.f new file mode 100644 index 00000000..caa4781e --- /dev/null +++ b/plugins/Symmetry/nuclei.irp.f @@ -0,0 +1,92 @@ + +BEGIN_PROVIDER [ double precision, nucl_coord_sym, (nucl_num,3) ] + implicit none + + BEGIN_DOC + ! Nuclear coordinates in standard orientation + END_DOC + + if (mpi_master) then + integer :: i + do i=1,nucl_num + nucl_coord_sym(i,1) = (nucl_coord(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,1) + & + (nucl_coord(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,1) + & + (nucl_coord(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,1) + nucl_coord_sym(i,2) = (nucl_coord(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,2) + & + (nucl_coord(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,2) + & + (nucl_coord(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,2) + nucl_coord_sym(i,3) = (nucl_coord(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,3) + & + (nucl_coord(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,3) + & + (nucl_coord(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,3) + enddo + + character*(64), parameter :: f = '(A16, 4(1X,F12.6))' + character*(64), parameter :: ft= '(A16, 4(1X,A12 ))' + double precision, parameter :: a0= 0.529177249d0 + + call write_time(output_Nuclei) + write(output_Nuclei,'(A)') '' + write(output_Nuclei,'(A)') 'Nuclear Coordinates in standard orientation (Angstroms)' + write(output_Nuclei,'(A)') '=======================================================' + write(output_Nuclei,'(A)') '' + write(output_Nuclei,ft) & + '================','============','============','============','============' + write(output_Nuclei,*) & + ' Atom Charge X Y Z ' + write(output_Nuclei,ft) & + '================','============','============','============','============' + do i=1,nucl_num + write(output_Nuclei,f) nucl_label(i), nucl_charge(i), & + nucl_coord_sym(i,1)*a0, & + nucl_coord_sym(i,2)*a0, & + nucl_coord_sym(i,3)*a0 + enddo + write(output_Nuclei,ft) & + '================','============','============','============','============' + write(output_Nuclei,'(A)') '' + + endif + + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( nucl_coord_sym, 3*nucl_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read nucl_coord_sym with MPI' + endif + IRP_ENDIF + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, nucl_coord_sym_transp, (3,nucl_num) ] + implicit none + BEGIN_DOC + ! Transposed array of nucl_coord + END_DOC + integer :: i, k + nucl_coord_sym_transp = 0.d0 + + do i=1,nucl_num + nucl_coord_sym_transp(1,i) = nucl_coord_sym(i,1) + nucl_coord_sym_transp(2,i) = nucl_coord_sym(i,2) + nucl_coord_sym_transp(3,i) = nucl_coord_sym(i,3) + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, sym_molecule_rotation, (3,3) ] + implicit none + BEGIN_DOC + ! Rotation of the molecule to go from input orientation to standard orientation + END_DOC + call find_rotation(nucl_coord, size(nucl_coord,1), nucl_coord_sym, 3, sym_molecule_rotation, 3) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, sym_molecule_rotation_inv, (3,3) ] + implicit none + BEGIN_DOC + ! Rotation of the molecule to go from standard orientation to input orientation + END_DOC + call find_rotation(nucl_coord_sym, size(nucl_coord_sym,1), nucl_coord, 3, sym_molecule_rotation_inv, 3) +END_PROVIDER + diff --git a/scripts/pseudo/elts_num_ele.py b/scripts/pseudo/elts_num_ele.py index 8f31f4f7..f0aa3179 100644 --- a/scripts/pseudo/elts_num_ele.py +++ b/scripts/pseudo/elts_num_ele.py @@ -1,119 +1,15 @@ -name_to_elec = {"X": 0, - "H": 1, - "He": 2, - "Li": 3, - "Be": 4, - "B": 5, - "C": 6, - "N": 7, - "O": 8, - "F": 9, - "Ne": 10, - "Na": 11, - "Mg": 12, - "Al": 13, - "Si": 14, - "P": 15, - "S": 16, - "Cl": 17, - "Ar": 18, - "K": 19, - "Ca": 20, - "Sc": 21, - "Ti": 22, - "V": 23, - "Cr": 24, - "Mn": 25, - "Fe": 26, - "Co": 27, - "Ni": 28, - "Cu": 29, - "Zn": 30, - "Ga": 31, - "Ge": 32, - "As": 33, - "Se": 34, - "Br": 35, - "Kr": 36, - "Rb": 37, - "Sr": 38, - "Y": 39, - "Zr": 40, - "Nb": 41, - "Mo": 42, - "Tc": 43, - "Ru": 44, - "Rh": 45, - "Pd": 46, - "Ag": 47, - "Cd": 48, - "In": 49, - "Sn": 50, - "Sb": 51, - "Te": 52, - "I": 53, - "Xe": 54, - "Cs": 55, - "Ba": 56, - "La": 57, - "Ce": 58, - "Pr": 59, - "Nd": 60, - "Pm": 61, - "Sm": 62, - "Eu": 63, - "Gd": 64, - "Tb": 65, - "Dy": 66, - "Ho": 67, - "Er": 68, - "Tm": 69, - "Yb": 70, - "Lu": 71, - "Hf": 72, - "Ta": 73, - "W": 74, - "Re": 75, - "Os": 76, - "Ir": 77, - "Pt": 78, - "Au": 79, - "Hg": 80, - "Tl": 81, - "Pb": 82, - "Bi": 83, - "Po": 84, - "At": 85, - "Rn": 86, - "Fr": 87, - "Ra": 88, - "Ac": 89, - "Th": 90, - "Pa": 91, - "U": 92, - "Np": 93, - "Pu": 94, - "Am": 95, - "Cm": 96, - "Bk": 97, - "Cf": 98, - "Es": 99, - "Fm": 100, - "Md": 101, - "No": 102, - "Lr": 103, - "Rf": 104, - "Db": 105, - "Sg": 106, - "Bh": 107, - "Hs": 108, - "Mt": 109, - "Ds": 110, - "Rg": 111, - "Cn": 112, - "Uut": 113, - "Fl": 114, - "Uup": 115, - "Lv": 116, - "Uus": 117, - "Uuo": 118} +#!/usr/bin/env python + +import os + +QP_ROOT=os.environ["QP_ROOT"] + +name_to_elec = {} +with open(QP_ROOT+"/data/list_element.txt","r") as f: + data = f.readlines() + for line in data: + b = line.split() + name_to_elec[b[1]] = int(b[0]) + +if __name__ == '__main__': + print name_to_elec diff --git a/src/Nuclei/inertia.irp.f b/src/Nuclei/inertia.irp.f index c4517cde..d2cb7cf8 100644 --- a/src/Nuclei/inertia.irp.f +++ b/src/Nuclei/inertia.irp.f @@ -6,12 +6,12 @@ BEGIN_PROVIDER [ double precision, inertia_tensor, (3,3) ] integer :: i,j,k inertia_tensor = 0.d0 do k=1,nucl_num - inertia_tensor(1,1) += mass(int(nucl_charge(k))) * ((nucl_coord_input(k,2)-center_of_mass(2))**2 + (nucl_coord_input(k,3)-center_of_mass(3))**2) - inertia_tensor(2,2) += mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1))**2 + (nucl_coord_input(k,3)-center_of_mass(3))**2) - inertia_tensor(3,3) += mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1))**2 + (nucl_coord_input(k,2)-center_of_mass(2))**2) - inertia_tensor(1,2) -= mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1)) * (nucl_coord_input(k,2)-center_of_mass(2)) ) - inertia_tensor(1,3) -= mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1)) * (nucl_coord_input(k,3)-center_of_mass(3)) ) - inertia_tensor(2,3) -= mass(int(nucl_charge(k))) * ((nucl_coord_input(k,2)-center_of_mass(2)) * (nucl_coord_input(k,3)-center_of_mass(3)) ) + inertia_tensor(1,1) += element_mass(int(nucl_charge(k))) * ((nucl_coord(k,2)-center_of_mass(2))**2 + (nucl_coord(k,3)-center_of_mass(3))**2) + inertia_tensor(2,2) += element_mass(int(nucl_charge(k))) * ((nucl_coord(k,1)-center_of_mass(1))**2 + (nucl_coord(k,3)-center_of_mass(3))**2) + inertia_tensor(3,3) += element_mass(int(nucl_charge(k))) * ((nucl_coord(k,1)-center_of_mass(1))**2 + (nucl_coord(k,2)-center_of_mass(2))**2) + inertia_tensor(1,2) -= element_mass(int(nucl_charge(k))) * ((nucl_coord(k,1)-center_of_mass(1)) * (nucl_coord(k,2)-center_of_mass(2)) ) + inertia_tensor(1,3) -= element_mass(int(nucl_charge(k))) * ((nucl_coord(k,1)-center_of_mass(1)) * (nucl_coord(k,3)-center_of_mass(3)) ) + inertia_tensor(2,3) -= element_mass(int(nucl_charge(k))) * ((nucl_coord(k,2)-center_of_mass(2)) * (nucl_coord(k,3)-center_of_mass(3)) ) enddo inertia_tensor(2,1) = inertia_tensor(1,2) inertia_tensor(3,1) = inertia_tensor(1,3) @@ -27,6 +27,6 @@ END_PROVIDER integer :: k call lapack_diagd(inertia_tensor_eigenvalues,inertia_tensor_eigenvectors,inertia_tensor,3,3) print *, 'Rotational constants (GHZ):' - print *, (1805.65468542d0/inertia_tensor_eigenvalues(k), k=3,1,-1) + print *, (1805.65468542d0/(inertia_tensor_eigenvalues(k)+1.d-32), k=3,1,-1) END_PROVIDER diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f index 43609506..910e9167 100644 --- a/src/Nuclei/nuclei.irp.f +++ b/src/Nuclei/nuclei.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [ double precision, nucl_coord_input, (nucl_num,3) ] +BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num,3) ] implicit none BEGIN_DOC @@ -8,7 +8,7 @@ BEGIN_PROVIDER [ double precision, nucl_coord_input, (nucl_num,3) ] if (mpi_master) then double precision, allocatable :: buffer(:,:) - nucl_coord_input = 0.d0 + nucl_coord = 0.d0 allocate (buffer(nucl_num,3)) buffer = 0.d0 logical :: has @@ -22,7 +22,7 @@ BEGIN_PROVIDER [ double precision, nucl_coord_input, (nucl_num,3) ] do i=1,3 do j=1,nucl_num - nucl_coord_input(j,i) = buffer(j,i) + nucl_coord(j,i) = buffer(j,i) enddo enddo deallocate(buffer) @@ -31,65 +31,6 @@ BEGIN_PROVIDER [ double precision, nucl_coord_input, (nucl_num,3) ] character*(64), parameter :: ft= '(A16, 4(1X,A12 ))' double precision, parameter :: a0= 0.529177249d0 - call write_time(output_Nuclei) - write(output_Nuclei,'(A)') '' - write(output_Nuclei,'(A)') 'Input Nuclear Coordinates (Angstroms)' - write(output_Nuclei,'(A)') '=====================================' - write(output_Nuclei,'(A)') '' - write(output_Nuclei,ft) & - '================','============','============','============','============' - write(output_Nuclei,*) & - ' Atom Charge X Y Z ' - write(output_Nuclei,ft) & - '================','============','============','============','============' - do i=1,nucl_num - write(output_Nuclei,f) nucl_label(i), nucl_charge(i), & - nucl_coord_input(i,1)*a0, & - nucl_coord_input(i,2)*a0, & - nucl_coord_input(i,3)*a0 - enddo - write(output_Nuclei,ft) & - '================','============','============','============','============' - write(output_Nuclei,'(A)') '' - - endif - - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( nucl_coord_input, 3*nucl_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read nucl_coord_input with MPI' - endif - IRP_ENDIF - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num,3) ] - implicit none - - BEGIN_DOC - ! Nuclear coordinates in standard orientation - END_DOC - - if (mpi_master) then - integer :: i - do i=1,nucl_num - nucl_coord(i,1) = (nucl_coord_input(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,1) + & - (nucl_coord_input(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,1) + & - (nucl_coord_input(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,1) - nucl_coord(i,2) = (nucl_coord_input(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,2) + & - (nucl_coord_input(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,2) + & - (nucl_coord_input(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,2) - nucl_coord(i,3) = (nucl_coord_input(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,3) + & - (nucl_coord_input(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,3) + & - (nucl_coord_input(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,3) - enddo - - character*(64), parameter :: f = '(A16, 4(1X,F12.6))' - character*(64), parameter :: ft= '(A16, 4(1X,A12 ))' - double precision, parameter :: a0= 0.529177249d0 - call write_time(output_Nuclei) write(output_Nuclei,'(A)') '' write(output_Nuclei,'(A)') 'Nuclear Coordinates (Angstroms)' @@ -207,7 +148,7 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] PROVIDE mpi_master nucl_coord nucl_charge nucl_num if (disk_access_nuclear_repulsion.EQ.'Read') then logical :: has - + if (mpi_master) then call ezfio_has_nuclei_nuclear_repulsion(has) if (has) then @@ -229,184 +170,78 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] 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 + + 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 if (mpi_master) then - call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion) + call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion) endif endif + + END_PROVIDER -BEGIN_PROVIDER [ character*(128), element_name, (78)] - BEGIN_DOC - ! Array of the name of element, sorted by nuclear charge (integer) - END_DOC - element_name(1) = 'H' - element_name(2) = 'He' - element_name(3) = 'Li' - element_name(4) = 'Be' - element_name(5) = 'B' - element_name(6) = 'C' - element_name(7) = 'N' - element_name(8) = 'O' - element_name(9) = 'F' - element_name(10) = 'Ne' - element_name(11) = 'Na' - element_name(12) = 'Mg' - element_name(13) = 'Al' - element_name(14) = 'Si' - element_name(15) = 'P' - element_name(16) = 'S' - element_name(17) = 'Cl' - element_name(18) = 'Ar' - element_name(19) = 'K' - element_name(20) = 'Ca' - element_name(21) = 'Sc' - element_name(22) = 'Ti' - element_name(23) = 'V' - element_name(24) = 'Cr' - element_name(25) = 'Mn' - element_name(26) = 'Fe' - element_name(27) = 'Co' - element_name(28) = 'Ni' - element_name(29) = 'Cu' - element_name(30) = 'Zn' - element_name(31) = 'Ga' - element_name(32) = 'Ge' - element_name(33) = 'As' - element_name(34) = 'Se' - element_name(35) = 'Br' - element_name(36) = 'Kr' - element_name(37) = 'Rb' - element_name(38) = 'Sr' - element_name(39) = 'Y' - element_name(40) = 'Zr' - element_name(41) = 'Nb' - element_name(42) = 'Mo' - element_name(43) = 'Tc' - element_name(44) = 'Ru' - element_name(45) = 'Rh' - element_name(46) = 'Pd' - element_name(47) = 'Ag' - element_name(48) = 'Cd' - element_name(49) = 'In' - element_name(50) = 'Sn' - element_name(51) = 'Sb' - element_name(52) = 'Te' - element_name(53) = 'I' - element_name(54) = 'Xe' - element_name(55) = 'Cs' - element_name(56) = 'Ba' - element_name(57) = 'La' - element_name(58) = 'Ce' - element_name(59) = 'Pr' - element_name(60) = 'Nd' - element_name(61) = 'Pm' - element_name(62) = 'Sm' - element_name(63) = 'Eu' - element_name(64) = 'Gd' - element_name(65) = 'Tb' - element_name(66) = 'Dy' - element_name(67) = 'Ho' - element_name(68) = 'Er' - element_name(69) = 'Tm' - element_name(70) = 'Yb' - element_name(71) = 'Lu' - element_name(72) = 'Hf' - element_name(73) = 'Ta' - element_name(74) = 'W' - element_name(75) = 'Re' - element_name(76) = 'Os' - element_name(77) = 'Ir' - element_name(78) = 'Pt' + BEGIN_PROVIDER [ character*(4), element_name, (0:128)] +&BEGIN_PROVIDER [ double precision, element_mass, (0:128) ] + BEGIN_DOC + ! Array of the name of element, sorted by nuclear charge (integer) + END_DOC + integer :: iunit + integer, external :: getUnitAndOpen + character*(128) :: filename + if (mpi_master) then + call getenv('QP_ROOT',filename) + filename = trim(filename)//'/data/list_element.txt' + iunit = getUnitAndOpen(filename,'r') + element_mass(:) = 0.d0 + do i=0,128 + write(element_name(i),'(I4)') i + enddo + character*(80) :: buffer, dummy + do + read(iunit,'(A80)',end=10) buffer + read(buffer,*) i ! First read i + read(buffer,*) i, element_name(i), dummy, element_mass(i) + enddo + 10 continue + close(10) + endif -END_PROVIDER - -BEGIN_PROVIDER [ double precision, mass, (0:110) ] - implicit none - BEGIN_DOC - ! Atomic masses - END_DOC - - mass( 0 ) = 0. - mass( 1 ) = 1.0079 - mass( 2 ) = 4.00260 - mass( 3 ) = 6.941 - mass( 4 ) = 9.01218 - mass( 5 ) = 10.81 - mass( 6 ) = 12.011 - mass( 7 ) = 14.0067 - mass( 8 ) = 15.9994 - mass( 9 ) = 18.998403 - mass( 10 ) = 20.179 - mass( 11 ) = 22.98977 - mass( 12 ) = 24.305 - mass( 13 ) = 26.98154 - mass( 14 ) = 28.0855 - mass( 15 ) = 30.97376 - mass( 16 ) = 32.06 - mass( 17 ) = 35.453 - mass( 18 ) = 39.948 - mass( 19 ) = 39.0983 - mass( 20 ) = 40.08 - mass( 21 ) = 44.9559 - mass( 22 ) = 47.90 - mass( 23 ) = 50.9415 - mass( 24 ) = 51.996 - mass( 25 ) = 54.9380 - mass( 26 ) = 55.9332 - mass( 27 ) = 58.9332 - mass( 28 ) = 58.70 - mass( 29 ) = 63.546 - mass( 30 ) = 65.38 - mass( 31 ) = 69.72 - mass( 32 ) = 72.59 - mass( 33 ) = 74.9216 - mass( 34 ) = 78.96 - mass( 35 ) = 79.904 - mass( 36 ) = 83.80 - mass( 37 ) = 85.4678 - mass( 38 ) = 87.62 - mass( 39 ) = 88.90584 - mass( 40 ) = 91.224 - mass( 41 ) = 92.90637 - mass( 42 ) = 95.95 - mass( 43 ) = 98. - mass( 44 ) = 101.07 - mass( 45 ) = 102.90550 - mass( 46 ) = 106.42 - mass( 47 ) = 107.8682 - mass( 48 ) = 112.414 - mass( 49 ) = 114.818 - mass( 50 ) = 118.710 - mass( 51 ) = 121.760 - mass( 52 ) = 127.60 - mass( 53 ) = 126.90447 - mass( 54 ) = 131.293 - mass( 78 ) = 195.084 + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( element_name, size(element_name)*4, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read element_name with MPI' + endif + call MPI_BCAST( element_mass, size(element_mass), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read element_name with MPI' + endif + IRP_ENDIF + END_PROVIDER BEGIN_PROVIDER [ double precision, center_of_mass, (3) ] @@ -420,9 +255,9 @@ BEGIN_PROVIDER [ double precision, center_of_mass, (3) ] s = 0.d0 do i=1,nucl_num do j=1,3 - center_of_mass(j) += nucl_coord_input(i,j)* mass(int(nucl_charge(i))) + center_of_mass(j) += nucl_coord(i,j)* element_mass(int(nucl_charge(i))) enddo - s += mass(int(nucl_charge(i))) + s += element_mass(int(nucl_charge(i))) enddo s = 1.d0/s center_of_mass(:) = center_of_mass(:)*s