From 28afcc8c59a2271e164cceaf1c8c1ebfbca26db9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Dec 2017 09:59:46 +0100 Subject: [PATCH] Fixes #211 --- 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 - scripts/pseudo/elts_num_ele.py | 134 +++------------------- src/Nuclei/inertia.irp.f | 14 +-- src/Nuclei/nuclei.irp.f | 168 ++++------------------------ 8 files changed, 132 insertions(+), 392 deletions(-) create mode 100644 ocaml/element_create_db.ml 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/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..97e32711 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_input(k,2)-center_of_mass(2))**2 + (nucl_coord_input(k,3)-center_of_mass(3))**2) + inertia_tensor(2,2) += element_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) += element_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) -= element_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) -= element_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) -= element_mass(int(nucl_charge(k))) * ((nucl_coord_input(k,2)-center_of_mass(2)) * (nucl_coord_input(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..df5fcd00 100644 --- a/src/Nuclei/nuclei.irp.f +++ b/src/Nuclei/nuclei.irp.f @@ -260,155 +260,33 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] endif END_PROVIDER -BEGIN_PROVIDER [ character*(128), element_name, (78)] + 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 - 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' + integer :: iunit + integer, external :: getUnitAndOpen + character*(128) :: filename + 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) + print *, i, element_name(i), element_mass(i) + enddo + 10 continue + close(10) 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 -END_PROVIDER - BEGIN_PROVIDER [ double precision, center_of_mass, (3) ] implicit none BEGIN_DOC @@ -420,9 +298,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_input(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