From 28afcc8c59a2271e164cceaf1c8c1ebfbca26db9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Dec 2017 09:59:46 +0100 Subject: [PATCH 1/9] 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 From ac37a499d32732f79e76a8a5f7ac8fafece12dbf Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Dec 2017 10:07:13 +0100 Subject: [PATCH 2/9] Removed debug print --- src/Nuclei/nuclei.irp.f | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f index df5fcd00..0d14ae7e 100644 --- a/src/Nuclei/nuclei.irp.f +++ b/src/Nuclei/nuclei.irp.f @@ -280,7 +280,6 @@ END_PROVIDER 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) From 17e0518410142a3621854b45d70863359a1662eb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Dec 2017 11:49:48 +0100 Subject: [PATCH 3/9] Revert input coordinates --- 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 +++++++++++ src/Nuclei/inertia.irp.f | 12 +- src/Nuclei/nuclei.irp.f | 223 ++++++++++++--------------- 6 files changed, 235 insertions(+), 165 deletions(-) create mode 100644 plugins/Symmetry/nuclei.irp.f 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/src/Nuclei/inertia.irp.f b/src/Nuclei/inertia.irp.f index 97e32711..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) += 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)) ) + 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) diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f index 0d14ae7e..7e4de39f 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)' @@ -205,35 +146,36 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] END_DOC 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 - call ezfio_get_nuclei_nuclear_repulsion(nuclear_repulsion) - else - print *, 'nuclei/nuclear_repulsion not found in EZFIO file' - stop 1 + if (mpi_master) then + if (disk_access_nuclear_repulsion.EQ.'Read') then + logical :: has + + if (mpi_master) then + 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 + print*, 'Read nuclear_repulsion' endif - print*, 'Read nuclear_repulsion' - endif - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( nuclear_repulsion, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read nuclear_repulsion with MPI' - endif - IRP_ENDIF - - - else - - integer :: k,l - double precision :: Z12, r2, x(3) - nuclear_repulsion = 0.d0 - do l = 1, nucl_num + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( nuclear_repulsion, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read nuclear_repulsion with MPI' + endif + IRP_ENDIF + + + 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 @@ -245,45 +187,76 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] 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) + 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) + endif endif - endif + + endif + + 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 [ character*(4), element_name, (0:128)] + 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 - 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) + 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 + 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) ] @@ -297,7 +270,7 @@ 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)* element_mass(int(nucl_charge(i))) + center_of_mass(j) += nucl_coord(i,j)* element_mass(int(nucl_charge(i))) enddo s += element_mass(int(nucl_charge(i))) enddo From d031148512c5b4d848896ce67273df7f680b0515 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Dec 2017 11:54:20 +0100 Subject: [PATCH 4/9] Fix MPI --- src/Nuclei/nuclei.irp.f | 119 ++++++++++++++++++---------------------- 1 file changed, 52 insertions(+), 67 deletions(-) diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f index 7e4de39f..910e9167 100644 --- a/src/Nuclei/nuclei.irp.f +++ b/src/Nuclei/nuclei.irp.f @@ -146,76 +146,61 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] END_DOC PROVIDE mpi_master nucl_coord nucl_charge nucl_num - if (mpi_master) then - if (disk_access_nuclear_repulsion.EQ.'Read') then - logical :: has - - if (mpi_master) then - 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 - print*, 'Read nuclear_repulsion' - endif - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( nuclear_repulsion, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read nuclear_repulsion with MPI' - endif - IRP_ENDIF - - - 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 - end if + if (disk_access_nuclear_repulsion.EQ.'Read') then + logical :: has - 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) + if (mpi_master) then + 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 + print*, 'Read nuclear_repulsion' endif - - endif - - 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 - + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( nuclear_repulsion, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read nuclear_repulsion with MPI' + endif + IRP_ENDIF + + + 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 + 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) + endif + endif + + END_PROVIDER BEGIN_PROVIDER [ character*(4), element_name, (0:128)] From f838765834cb456fcf34735b2cb8495cd72d944e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Dec 2017 18:26:24 +0100 Subject: [PATCH 5/9] Multi-state modification for qmcpack --- plugins/QMC/README.rst | 8 ++++++++ plugins/QMC/qp_convert_qmcpack_to_ezfio.py | 5 ++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/plugins/QMC/README.rst b/plugins/QMC/README.rst index 7e942878..129ca95e 100644 --- a/plugins/QMC/README.rst +++ b/plugins/QMC/README.rst @@ -2,6 +2,14 @@ QmcChem Module ============== +For multi-state calculations, to extract state 2 use: + +`` +QP_STATE=2 qp_run save_for_qmcpack x.ezfio +`` +(state 1 is the ground state). + + Documentation ============= diff --git a/plugins/QMC/qp_convert_qmcpack_to_ezfio.py b/plugins/QMC/qp_convert_qmcpack_to_ezfio.py index 8e8cea7b..acf48708 100755 --- a/plugins/QMC/qp_convert_qmcpack_to_ezfio.py +++ b/plugins/QMC/qp_convert_qmcpack_to_ezfio.py @@ -345,7 +345,10 @@ print "mo_num", mo_num print "det_num", n_det print "" -state = 0 +if "QP_STATE" in os.environ: + state = int(os.environ["QP_STATE"])-1 +else: + state = 0 psi_coef = psi_coef[state] encode = 8*bit_kind From c4a2b1f9a0e5819c5451beb0ca9faf73d1e9524f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Dec 2017 18:28:53 +0100 Subject: [PATCH 6/9] Removed useless file --- data/list_element.txt | 56 -------------------------------- plugins/Symmetry/aos.irp.f | 57 +++++++++++++++++++++------------ plugins/Symmetry/find_sym.irp.f | 25 +++++++++------ 3 files changed, 52 insertions(+), 86 deletions(-) delete mode 100644 data/list_element.txt diff --git a/data/list_element.txt b/data/list_element.txt deleted file mode 100644 index 2a9733e6..00000000 --- a/data/list_element.txt +++ /dev/null @@ -1,56 +0,0 @@ - 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/plugins/Symmetry/aos.irp.f b/plugins/Symmetry/aos.irp.f index e0aade1e..a20c1a08 100644 --- a/plugins/Symmetry/aos.irp.f +++ b/plugins/Symmetry/aos.irp.f @@ -24,9 +24,13 @@ subroutine generate_sym_coord(n_sym_points,result) END_DOC integer :: i, xyz - call random_number(result) + do i=1,n_sym_points + call random_number(result(1,i)) + call random_number(result(2,i)) + call random_number(result(3,i)) + enddo do xyz=1,3 - result(xyz,:) = sym_box(xyz,1) + result(xyz,:) * (sym_box(xyz,2)-sym_box(xyz,1)) + result(xyz,1:n_sym_points) = sym_box(xyz,1) + result(xyz,:) * (sym_box(xyz,2)-sym_box(xyz,1)) enddo end @@ -43,27 +47,35 @@ subroutine compute_sym_ao_values(sym_points, n_sym_points, result) integer :: i, j double precision :: x, y, z double precision :: x2, y2, z2 + integer :: k result (:,:) = 0.d0 +print *, sym_molecule_rotation_inv +print *, '' +print *, sym_molecule_rotation +stop do j=1,ao_num do i=1,n_sym_points - 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 + x2 = sym_points(1,i) + y2 = sym_points(2,i) + z2 = sym_points(3,i) + x = x2*sym_molecule_rotation_inv(1,1) + y2*sym_molecule_rotation_inv(2,1) + z2*sym_molecule_rotation_inv(3,1) + y = x2*sym_molecule_rotation_inv(1,2) + y2*sym_molecule_rotation_inv(2,2) + z2*sym_molecule_rotation_inv(3,2) + z = x2*sym_molecule_rotation_inv(1,3) + y2*sym_molecule_rotation_inv(2,3) + z2*sym_molecule_rotation_inv(3,3) + x = x - nucl_coord_transp(1,ao_nucl(j)) + y = y - nucl_coord_transp(2,ao_nucl(j)) + z = z - nucl_coord_transp(3,ao_nucl(j)) + x2 = x*x + y*y + z*z + result(i,j) = 0.d0 + do k=1,ao_prim_num(j) + result(i,j) += ao_coef_normalized_ordered_transp(k,j)*exp(-ao_expo_ordered_transp(k,j)*x2) + enddo +print *, real(x), ao_power(j,1), real(y), ao_power(j,2), real(z), ao_power(j,3) + x = x**ao_power(j,1) + y = y**ao_power(j,2) + z = z**ao_power(j,3) +print *, result(i,j) + result(i,j) = x*y*z*result(i,j) + print *, result(i,j) enddo enddo @@ -81,6 +93,11 @@ subroutine compute_sym_mo_values(sym_points, n_sym_points, result) double precision, allocatable :: tmp(:,:) allocate(tmp(n_sym_points,ao_num)) call compute_sym_ao_values(sym_points,n_sym_points,tmp) +integer :: i +do i=1,ao_num + print *, tmp(:,i) +enddo +stop call dgemm('N','N',n_sym_points,mo_tot_num,ao_num, & 1.d0, tmp,size(tmp,1), mo_coef, size(mo_coef,1), & 0.d0, result,size(result,1)) diff --git a/plugins/Symmetry/find_sym.irp.f b/plugins/Symmetry/find_sym.irp.f index 95bb4e6c..c7e2ba37 100644 --- a/plugins/Symmetry/find_sym.irp.f +++ b/plugins/Symmetry/find_sym.irp.f @@ -369,7 +369,7 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] double precision :: sym_operations_on_mos(mo_tot_num) logical :: possible_irrep(n_irrep,mo_tot_num) - n_sym_points = 10000 + n_sym_points = 10 allocate(val(n_sym_points,mo_tot_num,2), sym_points(3,n_sym_points), ref_points(3,n_sym_points)) call generate_sym_coord(n_sym_points,ref_points) @@ -402,10 +402,12 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] call sym_apply_diagonal_reflexion(molecule_principal_axis,ref_points(1,ipoint),sym_points(1,ipoint)) enddo else if (sym_operation(iop) == 'C2''') then + angle = 2.d0 do ipoint=1,n_sym_points call sym_apply_rotation(angle,molecule_secondary_axis,ref_points(1,ipoint),sym_points(1,ipoint)) enddo else if (sym_operation(iop) == 'C2"') then + angle = 2.d0 do ipoint=1,n_sym_points call sym_apply_rotation(angle,molecule_ternary_axis,ref_points(1,ipoint),sym_points(1,ipoint)) enddo @@ -414,9 +416,12 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] if (sym_operation(iop)(l:l) == '^') exit enddo read(sym_operation(iop)(2:l-1), *) iangle - l=1 - read(sym_operation(iop)(l+1:), *, err=10, end=10) l - 10 continue + if (l == len(sym_operation(iop))+1) then + l=1 + else + read(sym_operation(iop)(l+1:), *, err=10, end=10) l + 10 continue + endif angle = dble(iangle)/(dble(l)) if (sym_operation(iop)(1:1) == 'C') then do ipoint=1,n_sym_points @@ -432,23 +437,23 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] call compute_sym_mo_values(sym_points,n_sym_points,val(1,1,1)) print *, sym_operation(iop) + double precision :: icount do imo=1,mo_tot_num sym_operations_on_mos(imo) = 0.d0 + icount = 0 do ipoint=1,n_sym_points double precision :: x + if (dabs(val(ipoint,imo,1)) < 1.d-5) cycle + icount += 1.d0 x = val(ipoint,imo,1)/val(ipoint,imo,2) if (dabs(x) > 1.d0) then x = 1.d0/x endif sym_operations_on_mos(imo) += x enddo - sym_operations_on_mos(imo) *= 1.d0/n_sym_points + sym_operations_on_mos(imo) *= 1.d0/icount 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 - sym_operations_on_mos(imo)=-1.d0 - endif + print *, val(:,imo,1) do i=1,n_irrep if (character_table(i,iop) /= sym_operations_on_mos(imo)) then possible_irrep(i,imo) = .False. From e8df53857ae14ade933a698a1a5779feb3ffceae Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Dec 2017 22:19:29 +0100 Subject: [PATCH 7/9] Fixed QMCPack pseudo bug --- plugins/QMC/qp_convert_qmcpack_to_ezfio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plugins/QMC/qp_convert_qmcpack_to_ezfio.py b/plugins/QMC/qp_convert_qmcpack_to_ezfio.py index acf48708..94f0c347 100755 --- a/plugins/QMC/qp_convert_qmcpack_to_ezfio.py +++ b/plugins/QMC/qp_convert_qmcpack_to_ezfio.py @@ -26,7 +26,7 @@ if do_pseudo: l_element_raw = data_raw.split("\n") l_element = [element_raw.split() for element_raw in l_element_raw] - d_z = dict((abr, z) for (z, abr, ele) in l_element) + d_z = dict((abr, z) for (z, abr, ele, _) in filter(lambda x: x != [], l_element) ) else: print "do_pseudo False" From 4a04e769c5ac2ea4bd0f99b6b84838574e75fc91 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 20 Dec 2017 18:02:55 +0100 Subject: [PATCH 8/9] pin old core library --- install/scripts/install_ocaml.sh | 2 +- plugins/Symmetry/Symmetry.main.irp.f | 2 +- plugins/Symmetry/aos.irp.f | 28 +++------- plugins/Symmetry/find_sym.irp.f | 27 ++++++---- plugins/Symmetry/nuclei.irp.f | 78 +++++++++++++++++++--------- plugins/Symmetry/sym_operation.irp.f | 25 ++++----- src/Nuclei/inertia.irp.f | 5 +- 7 files changed, 91 insertions(+), 76 deletions(-) diff --git a/install/scripts/install_ocaml.sh b/install/scripts/install_ocaml.sh index 88d38845..f322bd0b 100755 --- a/install/scripts/install_ocaml.sh +++ b/install/scripts/install_ocaml.sh @@ -5,7 +5,7 @@ QP_ROOT=$PWD cd - # Normal installation -PACKAGES="core cryptokit.1.10 ocamlfind sexplib ZMQ ppx_sexp_conv ppx_deriving" +PACKAGES="core.v0.9.1 cryptokit.1.10 ocamlfind sexplib.v0.9.1 ZMQ ppx_sexp_conv ppx_deriving" # Needed for ZeroMQ export C_INCLUDE_PATH="${QP_ROOT}"/include:"${C_INCLUDE_PATH}" diff --git a/plugins/Symmetry/Symmetry.main.irp.f b/plugins/Symmetry/Symmetry.main.irp.f index 1d438af5..9540295f 100644 --- a/plugins/Symmetry/Symmetry.main.irp.f +++ b/plugins/Symmetry/Symmetry.main.irp.f @@ -15,7 +15,7 @@ program Symmetry print *, 'Symmetry operations : ', sym_operation(1:n_irrep) print *, 'Character table' do i=1,n_irrep - print *, character_table(i,:) + print *, i, real(character_table(i,:)) enddo PROVIDE mo_sym end diff --git a/plugins/Symmetry/aos.irp.f b/plugins/Symmetry/aos.irp.f index 5515e1b7..1ed567bc 100644 --- a/plugins/Symmetry/aos.irp.f +++ b/plugins/Symmetry/aos.irp.f @@ -45,38 +45,27 @@ subroutine compute_sym_ao_values(sym_points, n_sym_points, result) double precision, intent(in) :: sym_points(3,n_sym_points) double precision, intent(out) :: result(n_sym_points, ao_num) integer :: i, j + double precision :: point(3) double precision :: x, y, z double precision :: x2, y2, z2 integer :: k - + result (:,:) = 0.d0 -print *, sym_molecule_rotation_inv -print *, '' -print *, sym_molecule_rotation -stop do j=1,ao_num do i=1,n_sym_points - x2 = sym_points(1,i) - y2 = sym_points(2,i) - z2 = sym_points(3,i) - x = x2*sym_molecule_rotation_inv(1,1) + y2*sym_molecule_rotation_inv(2,1) + z2*sym_molecule_rotation_inv(3,1) - y = x2*sym_molecule_rotation_inv(1,2) + y2*sym_molecule_rotation_inv(2,2) + z2*sym_molecule_rotation_inv(3,2) - z = x2*sym_molecule_rotation_inv(1,3) + y2*sym_molecule_rotation_inv(2,3) + z2*sym_molecule_rotation_inv(3,3) - x = x - nucl_coord_transp(1,ao_nucl(j)) - y = y - nucl_coord_transp(2,ao_nucl(j)) - z = z - nucl_coord_transp(3,ao_nucl(j)) + call point_to_input_orientation(sym_points(:,i), point) + x = point(1) - nucl_coord_transp(1,ao_nucl(j)) + y = point(2) - nucl_coord_transp(2,ao_nucl(j)) + z = point(3) - nucl_coord_transp(3,ao_nucl(j)) x2 = x*x + y*y + z*z result(i,j) = 0.d0 do k=1,ao_prim_num(j) result(i,j) += ao_coef_normalized_ordered_transp(k,j)*exp(-ao_expo_ordered_transp(k,j)*x2) enddo -print *, real(x), ao_power(j,1), real(y), ao_power(j,2), real(z), ao_power(j,3) x = x**ao_power(j,1) y = y**ao_power(j,2) z = z**ao_power(j,3) -print *, result(i,j) result(i,j) = x*y*z*result(i,j) - print *, result(i,j) enddo enddo @@ -94,11 +83,6 @@ subroutine compute_sym_mo_values(sym_points, n_sym_points, result) double precision, allocatable :: tmp(:,:) allocate(tmp(n_sym_points,ao_num)) call compute_sym_ao_values(sym_points,n_sym_points,tmp) -integer :: i -do i=1,ao_num - print *, tmp(:,i) -enddo -stop call dgemm('N','N',n_sym_points,mo_tot_num,ao_num, & 1.d0, tmp,size(tmp,1), mo_coef, size(mo_coef,1), & 0.d0, result,size(result,1)) diff --git a/plugins/Symmetry/find_sym.irp.f b/plugins/Symmetry/find_sym.irp.f index c7e2ba37..38c0a6b7 100644 --- a/plugins/Symmetry/find_sym.irp.f +++ b/plugins/Symmetry/find_sym.irp.f @@ -84,7 +84,7 @@ END_PROVIDER END_DOC molecule_principal_axis = maxloc(sym_rotation_axis,1) if (molecule_principal_axis == 1) then - if (sym_rotation_axis(2) >= sym_rotation_axis(3)) then + if (sym_rotation_axis(2) > sym_rotation_axis(3)) then molecule_secondary_axis = 2 molecule_ternary_axis = 3 else @@ -92,7 +92,7 @@ END_PROVIDER molecule_ternary_axis = 2 endif else if (molecule_principal_axis == 2) then - if (sym_rotation_axis(1) >= sym_rotation_axis(3)) then + if (sym_rotation_axis(1) > sym_rotation_axis(3)) then molecule_secondary_axis = 1 molecule_ternary_axis = 3 else @@ -100,7 +100,7 @@ END_PROVIDER molecule_ternary_axis = 1 endif else if (molecule_principal_axis == 3) then - if (sym_rotation_axis(1) >= sym_rotation_axis(2)) then + if (sym_rotation_axis(1) > sym_rotation_axis(2)) then molecule_secondary_axis = 1 molecule_ternary_axis = 2 else @@ -375,11 +375,11 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] call generate_sym_coord(n_sym_points,ref_points) call compute_sym_mo_values(ref_points,n_sym_points,val(1,1,2)) + possible_irrep = .True. do iop=1,n_irrep if (sym_operation(iop) == 'E') then cycle endif - possible_irrep = .True. if (sym_operation(iop) == 'i') then do ipoint=1,n_sym_points @@ -395,11 +395,12 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] enddo else if (sym_operation(iop) == 'sv') then do ipoint=1,n_sym_points - call sym_apply_reflexion(molecule_secondary_axis,ref_points(1,ipoint),sym_points(1,ipoint)) + call sym_apply_reflexion(molecule_ternary_axis,ref_points(1,ipoint),sym_points(1,ipoint)) enddo else if (sym_operation(iop) == 'sd') then + angle = dble(maxval(sym_rotation_axis)) do ipoint=1,n_sym_points - call sym_apply_diagonal_reflexion(molecule_principal_axis,ref_points(1,ipoint),sym_points(1,ipoint)) + call sym_apply_diagonal_reflexion(angle,molecule_principal_axis,ref_points(1,ipoint),sym_points(1,ipoint)) enddo else if (sym_operation(iop) == 'C2''') then angle = 2.d0 @@ -452,17 +453,23 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] sym_operations_on_mos(imo) += x enddo sym_operations_on_mos(imo) *= 1.d0/icount - print *, iop, imo, sym_operations_on_mos(imo) - print *, val(:,imo,1) + 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 + sym_operations_on_mos(imo) = -1.d0 + else if (dabs(sym_operations_on_mos(imo)) < 1.d-2) then + sym_operations_on_mos(imo) = 0.d0 + endif + print *, imo, sym_operations_on_mos(imo) do i=1,n_irrep - if (character_table(i,iop) /= sym_operations_on_mos(imo)) then + if (dabs(character_table(i,iop) - sym_operations_on_mos(imo)) > 1.d-2) then possible_irrep(i,imo) = .False. endif enddo enddo enddo do imo=1,mo_tot_num - print *, imo + print *, 'MO ', imo do i=1,n_irrep if (possible_irrep(i,imo)) then print *, sym_irrep(i) diff --git a/plugins/Symmetry/nuclei.irp.f b/plugins/Symmetry/nuclei.irp.f index caa4781e..d680393d 100644 --- a/plugins/Symmetry/nuclei.irp.f +++ b/plugins/Symmetry/nuclei.irp.f @@ -1,3 +1,56 @@ +subroutine point_to_standard_orientation(point_in,point_out) + implicit none + double precision, intent(in) :: point_in(3) + double precision, intent(out) :: point_out(3) + BEGIN_DOC + ! Returns the coordinates of a point in the standard orientation + END_DOC + double precision :: point_tmp(3) + + point_tmp(1) = point_in(1) - center_of_mass(1) + point_tmp(2) = point_in(2) - center_of_mass(2) + point_tmp(3) = point_in(3) - center_of_mass(3) + + point_out(1) = point_tmp(1)*inertia_tensor_eigenvectors(1,1) + & + point_tmp(2)*inertia_tensor_eigenvectors(2,1) + & + point_tmp(3)*inertia_tensor_eigenvectors(3,1) + + point_out(2) = point_tmp(1)*inertia_tensor_eigenvectors(1,2) + & + point_tmp(2)*inertia_tensor_eigenvectors(2,2) + & + point_tmp(3)*inertia_tensor_eigenvectors(3,2) + + point_out(3) = point_tmp(1)*inertia_tensor_eigenvectors(1,3) + & + point_tmp(2)*inertia_tensor_eigenvectors(2,3) + & + point_tmp(3)*inertia_tensor_eigenvectors(3,3) + +end + +subroutine point_to_input_orientation(point_in,point_out) + implicit none + double precision, intent(in) :: point_in(3) + double precision, intent(out) :: point_out(3) + BEGIN_DOC + ! Returns the coordinates of a point in the input orientation + END_DOC + double precision :: point_tmp(3) + + point_tmp(1) = point_in(1)*inertia_tensor_eigenvectors(1,1) + & + point_in(2)*inertia_tensor_eigenvectors(1,2) + & + point_in(3)*inertia_tensor_eigenvectors(1,3) + + point_tmp(2) = point_in(1)*inertia_tensor_eigenvectors(2,1) + & + point_in(2)*inertia_tensor_eigenvectors(2,2) + & + point_in(3)*inertia_tensor_eigenvectors(2,3) + + point_tmp(3) = point_in(1)*inertia_tensor_eigenvectors(3,1) + & + point_in(2)*inertia_tensor_eigenvectors(3,2) + & + point_in(3)*inertia_tensor_eigenvectors(3,3) + + point_out(1) = point_tmp(1) + center_of_mass(1) + point_out(2) = point_tmp(2) + center_of_mass(2) + point_out(3) = point_tmp(3) + center_of_mass(3) + +end BEGIN_PROVIDER [ double precision, nucl_coord_sym, (nucl_num,3) ] implicit none @@ -9,15 +62,7 @@ BEGIN_PROVIDER [ double precision, nucl_coord_sym, (nucl_num,3) ] 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) + call point_to_standard_orientation(nucl_coord(i,:), nucl_coord_sym(i,:)) enddo character*(64), parameter :: f = '(A16, 4(1X,F12.6))' @@ -74,19 +119,4 @@ BEGIN_PROVIDER [ double precision, nucl_coord_sym_transp, (3,nucl_num) ] 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/plugins/Symmetry/sym_operation.irp.f b/plugins/Symmetry/sym_operation.irp.f index 93e580e2..cfc86621 100644 --- a/plugins/Symmetry/sym_operation.irp.f +++ b/plugins/Symmetry/sym_operation.irp.f @@ -15,25 +15,18 @@ subroutine sym_apply_reflexion(iaxis,point_in,point_out) end -subroutine sym_apply_diagonal_reflexion(iaxis,point_in,point_out) +subroutine sym_apply_diagonal_reflexion(angle,iaxis,point_in,point_out) implicit none integer, intent(in) :: iaxis - double precision, intent(in) :: point_in(3) + double precision, intent(in) :: point_in(3), angle double precision, intent(out) :: point_out(3) - if (iaxis == 1) then - point_out(1) = point_in(1) - point_out(2) = point_in(3) - point_out(3) = point_in(2) - else if (iaxis == 2) then - point_out(1) = point_in(3) - point_out(2) = point_in(2) - point_out(3) = point_in(1) - else if (iaxis == 3) then - point_out(1) = point_in(2) - point_out(2) = point_in(1) - point_out(3) = point_in(3) - endif - + double precision :: point_tmp1(3), point_tmp2(3) + integer :: iaxis2 + iaxis2 = mod(iaxis,3)+1 + iaxis2 = mod(iaxis2,3)+1 + call sym_apply_rotation(-angle,iaxis,point_in,point_tmp1) + call sym_apply_reflexion(iaxis2,point_tmp1,point_tmp2) + call sym_apply_rotation(angle,iaxis,point_tmp2,point_out) end subroutine sym_apply_rotation(angle,iaxis,point_in,point_out) diff --git a/src/Nuclei/inertia.irp.f b/src/Nuclei/inertia.irp.f index d2cb7cf8..c79fa243 100644 --- a/src/Nuclei/inertia.irp.f +++ b/src/Nuclei/inertia.irp.f @@ -25,8 +25,9 @@ END_PROVIDER ! Eigenvectors/eigenvalues of the inertia_tensor. Used to find normal orientation. END_DOC integer :: k - call lapack_diagd(inertia_tensor_eigenvalues,inertia_tensor_eigenvectors,inertia_tensor,3,3) + call lapack_diagd(inertia_tensor_eigenvalues,inertia_tensor_eigenvectors,-inertia_tensor,3,3) + inertia_tensor_eigenvalues = -inertia_tensor_eigenvalues print *, 'Rotational constants (GHZ):' - print *, (1805.65468542d0/(inertia_tensor_eigenvalues(k)+1.d-32), k=3,1,-1) + print *, (1805.65468542d0/(inertia_tensor_eigenvalues(k)+1.d-32), k=1,3) END_PROVIDER