From 5035ffad287b5dd9205c4b54429ee0468efbee0f Mon Sep 17 00:00:00 2001 From: gorges <58849467+gorges97@users.noreply.github.com> Date: Fri, 11 Nov 2022 16:54:45 +0100 Subject: [PATCH 01/11] "-p" now plots only the highest peak of an isotope pattern --- src/isotopes.f90 | 70 +++++++++++++----------------------------------- 1 file changed, 18 insertions(+), 52 deletions(-) diff --git a/src/isotopes.f90 b/src/isotopes.f90 index 72aad2b..87c26c1 100644 --- a/src/isotopes.f90 +++ b/src/isotopes.f90 @@ -35,6 +35,8 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & real(wp) :: list_masses(nrnd) real(wp),allocatable :: isotope_masses(:) + real(wp) :: mipmax, iipmax ! mass and intensity of highest isotope peak + integer :: indexipmax ! index of peak with highest intensity real(wp) :: current_mass real(wp), allocatable :: exact_intensity(:) @@ -468,57 +470,8 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & enddo - !Calculate NO isotope pattern - if ( no_isotopes ) then - xmass=0 - do i = 1, ntot - iti = iat_save(i) - p1 = 0.0_wp - do iso = 1, niso(iti) - if ( prob(iti,iso) > p1 ) then - x = massiso(iti,iso) - p1 = prob(iti,iso) - endif - enddo - xmass = xmass + x - enddo - - current_mass = xmass / real(abs(z_chrg),wp) - - there = .true. - if ( current_mass > mzmin ) then - !> loop over all entries in the list -noiso: do - loop = loop + 1 - !>> write if there is no entry - if ( list_masses(loop) == 0.0_wp ) then - there = .false. - exit noiso - !>> true if already in list, end - elseif ( abs(list_masses(loop) - current_mass) < 1.0d-10 ) then - there = .true. - store_int(loop) = store_int(loop) + 1 - exit noiso - !>> false if not in list, store - elseif ( abs(list_masses(loop) - current_mass) > 1.0d-10 ) then - there = .false. - cycle noiso - endif - enddo noiso - loop = 0 - - !> if it is not in the list, add it - if ( .not. there ) then - index_mass = index_mass + 1 - list_masses(index_mass) = current_mass - store_int(index_mass) = store_int(index_mass) + 1 - !write(*,*) list_masses(index_mass) - endif - - endif - - ! Calculate isotope pattern - else + + ! Calculate isotope pattern ! loop over random number runs do n = 1, nrnd @@ -577,7 +530,6 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & enddo - endif allocate(exact_intensity(index_mass)) @@ -592,6 +544,20 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & !write(*,*) exact_intensity(loop) enddo + + + +! if no isotopes we take here the peak of the isotope pattern with the highest intensity +if ( no_isotopes ) then +iipmax = maxval(exact_intensity) +indexipmax = maxloc(exact_intensity, dim = 1) +mipmax = isotope_masses(indexipmax) +deallocate(exact_intensity,isotope_masses) +index_mass = 1 +allocate(exact_intensity(1),isotope_masses(1)) +exact_intensity(1)=iipmax +isotope_masses(1)=mipmax +endif end subroutine isotope From 78428c6b1b62c16aa04e23adf937904803f33c9f Mon Sep 17 00:00:00 2001 From: gorges <58849467+gorges97@users.noreply.github.com> Date: Mon, 12 Dec 2022 11:24:51 +0100 Subject: [PATCH 02/11] fixed formatting --- src/isotopes.f90 | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/isotopes.f90 b/src/isotopes.f90 index 87c26c1..f3894f4 100644 --- a/src/isotopes.f90 +++ b/src/isotopes.f90 @@ -547,17 +547,17 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & -! if no isotopes we take here the peak of the isotope pattern with the highest intensity -if ( no_isotopes ) then -iipmax = maxval(exact_intensity) -indexipmax = maxloc(exact_intensity, dim = 1) -mipmax = isotope_masses(indexipmax) -deallocate(exact_intensity,isotope_masses) -index_mass = 1 -allocate(exact_intensity(1),isotope_masses(1)) -exact_intensity(1)=iipmax -isotope_masses(1)=mipmax -endif + ! if no isotopes we take here the peak of the isotope pattern with the highest intensity + if ( no_isotopes ) then + iipmax = maxval(exact_intensity) + indexipmax = maxloc(exact_intensity, dim = 1) + mipmax = isotope_masses(indexipmax) + deallocate(exact_intensity,isotope_masses) + index_mass = 1 + allocate(exact_intensity(1),isotope_masses(1)) + exact_intensity(1)=iipmax + isotope_masses(1)=mipmax + endif end subroutine isotope From 9eb257f935ee2fc86b9b15d7cca76f2151847cd9 Mon Sep 17 00:00:00 2001 From: gorges <58849467+gorges97@users.noreply.github.com> Date: Mon, 25 Sep 2023 19:14:37 +0200 Subject: [PATCH 03/11] added stable isotopes for elements up to Bismuth --- src/isotopes.f90 | 475 ++++++++++++++++++++++++++++++++++++++++++++++- src/main.F90 | 10 +- 2 files changed, 469 insertions(+), 16 deletions(-) diff --git a/src/isotopes.f90 b/src/isotopes.f90 index f3894f4..692c25e 100644 --- a/src/isotopes.f90 +++ b/src/isotopes.f90 @@ -50,35 +50,56 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso=0 - ! 1 H (Hydrogen) + ! 1 H (Hydrogen) niso(1) = 2 prob(1,1) = 0.0115_wp prob(1,2) = 99.9885_wp massiso(1,1) = 2.014101_wp massiso(1,2) = 1.007825_wp - ! 5 B (Boron) + ! 2 He (Helium) + niso(2) = 2 + prob(2,1) = 0.000134_wp + prob(2,2) = 99.99866_wp + massiso(2,1) = 3.016029309_wp + massiso(2,2) = 4.002603249_wp + + ! 3 Li (Lithium) + + niso(3) = 2 + prob(3,1) = 7.59_wp + prob(3,2) = 92.41_wp + massiso(3,1) = 6.0151223_wp + massiso(3,2) = 7.0160041_wp + + ! 4 Be (Beryllium) + + niso(4) = 1 + prob(4,1) = 100.0_wp + massiso(4,1) = 9.0121822_wp + + ! 5 B (Boron) niso(5) = 2 prob(5,1) = 19.9_wp prob(5,2) = 80.1_wp massiso(5,1) = 10.0129370_wp massiso(5,2) = 11.0093055_wp - ! 6 C (Carbon) + ! 6 C (Carbon) niso(6) = 2 prob(6,1) = 1.07_wp prob(6,2) = 98.93_wp massiso(6,1) = 13.003354_wp massiso(6,2) = 12.00_wp - ! 7 N (Nitrogen) + ! 7 N (Nitrogen) niso(7) = 2 prob(7,1) = 0.368_wp prob(7,2) = 99.632_wp massiso(7,1) = 15.000108_wp massiso(7,2) = 14.003074_wp - ! 8 O (Oxygen) + ! 8 O (Oxygen) niso(8) = 3 prob(8,1) = 0.038_wp prob(8,2) = 0.205_wp @@ -87,11 +108,34 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(8,2) = 17.999159_wp massiso(8,3) = 15.994914_wp - ! 9 F (Fluorine) + ! 9 F (Fluorine) niso(9) = 1 prob(9,1) = 100.0_wp massiso(9,1) = 18.998403_wp - + + ! 10 Ne (Neon) + niso(10) = 3 + prob(10,1) = 90.48_wp + prob(10,2) = 0.27_wp + prob(10,3) = 9.25_wp + massiso(10,1) = 19.992440176_wp + massiso(10,2) = 20.99384674_wp + massiso(10,3) = 21.99138550_wp + ! 11 Na (Sodium) + niso(11) = 1 + prob(11,1) = 100.0_wp + massiso(11,1) = 22.98976966_wp + + ! 12 Mg (Magnesium) + niso(12) = 3 + prob(12,1) = 78.99_wp + prob(12,2) = 10.00_wp + prob(12,3) = 11.01_wp + massiso(12,1) = 23.98504187_wp + massiso(12,2) = 24.98583700_wp + massiso(12,3) = 25.98259300_wp + + ! 13 Al (Aluminium) niso(13) = 1 prob(13,1) = 100.0_wp @@ -129,6 +173,24 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(17,1) = 34.968852_wp massiso(17,2) = 36.965902_wp + ! 18 Ar (Argon) + niso(18) = 3 + prob(18,1) = 0.3365_wp + prob(18,2) = 0.0632_wp + prob(18,3) = 99.6003_wp + massiso(18,1) = 35.9675462_wp + massiso(18,2) = 37.9627322_wp + massiso(18,3) = 39.96238312_wp + + ! 19 K (Potassium) + niso(19) = 3 + prob(19,1) = 93.2581_wp + prob(19,2) = 0.0117_wp + prob(19,3) = 6.7302_wp + massiso(19,1) = 38.9637069_wp + massiso(19,2) = 39.96399867_wp + massiso(19,3) = 40.96182597_wp + ! 20 Ca (Calcium) niso(20) = 5 prob(20,1) = 96.94_wp @@ -140,7 +202,12 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(20,2) = 41.958618_wp massiso(20,3) = 42.958766_wp massiso(20,4) = 43.955482_wp - massiso(20,5) = 45.953688_wp + massiso(20,5) = 45.953688_wp + + ! 21 Sc (Scandium) + niso(21) = 1 + prob(21,1) = 100.0_wp + massiso(21,1) = 44.9559102_wp ! 22 Ti (Titanium) niso(22) = 5 @@ -154,6 +221,13 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(22,3) = 47.947946_wp massiso(22,4) = 48.947870_wp massiso(22,5) = 49.944791_wp + + ! 23 V (Vanadium) + niso(23) = 2 + prob(23,1) = 0.250_wp + prob(23,2) = 99.750_wp + massiso(23,1) = 49.9471627_wp + massiso(23,2) = 50.943963_wp ! 24 Cr (Chromium) niso(24) = 4 @@ -219,6 +293,14 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(30,3) = 66.927127_wp massiso(30,4) = 67.924884_wp massiso(30,5) = 69.925319_wp + + ! 31 Ga (Gallium) + + niso(31) = 2 + prob(31,1) = 60.108_wp + prob(31,2) = 39.892_wp + massiso(31,1) = 68.925581_wp + massiso(31,2) = 70.9247073_wp ! 32 Ge (Germanium) niso(32) = 5 @@ -259,6 +341,46 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & prob(35,2) = 49.31_wp massiso(35,1) = 78.91833_wp massiso(35,2) = 80.91629_wp + + ! 36 Kr (Krypton) + niso(36) = 6 + prob(36,1) = 0.355_wp + prob(36,2) = 2.286_wp + prob(36,3) = 11.593_wp + prob(36,4) = 11.500_wp + prob(36,5) = 56.987_wp + prob(36,6) = 17.279_wp + massiso(36,1) = 77.920388_wp + massiso(36,2) = 79.916379_wp + massiso(36,3) = 81.9134850_wp + massiso(36,4) = 82.914137_wp + massiso(36,5) = 83.911508_wp + massiso(36,6) = 85.910615_wp + + ! 37 Rb (Rubidium) + niso(37) = 2 + prob(37,1) = 72.17_wp + prob(37,2) = 27.83_wp + massiso(37,1) = 84.9117924_wp + massiso(37,2) = 86.9091858_wp + + ! 38 Sr (Strontium) + niso(38) = 4 + prob(38,1) = 0.56_wp + prob(38,2) = 9.86_wp + prob(38,3) = 7.00_wp + prob(38,4) = 82.58_wp + massiso(38,1) = 83.913426_wp + massiso(38,2) = 85.9092647_wp + massiso(38,3) = 86.908816_wp + massiso(38,4) = 87.9056167_wp + + ! 39 Y (Yttrium) + + niso(39) = 1 + prob(39,1) = 100.0_wp + massiso(39,1) = 88.9058485_wp + ! 40 Zr (Zirconium) niso(40) = 5 @@ -272,6 +394,35 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(40,3) = 91.905040_wp massiso(40,4) = 93.906315_wp massiso(40,5) = 95.908273_wp + + + ! 41 Nb (Niobium) + niso(41) = 1 + prob(41,1) = 100.0_wp + massiso(41,1) = 92.9063762_wp + + ! 42 Mo (Molybdenum) + niso(42) = 7 + prob(42,1) = 14.77_wp + prob(42,2) = 9.23_wp + prob(42,3) = 15.90_wp + prob(42,4) = 16.68_wp + prob(42,5) = 9.56_wp + prob(42,6) = 24.19_wp + prob(42,7) = 9.67_wp + massiso(42,1) = 91.906810_wp + massiso(42,2) = 93.9050876_wp + massiso(42,3) = 94.9058406_wp + massiso(42,4) = 95.9046780_wp + massiso(42,5) = 96.9060201_wp + massiso(42,6) = 97.9054069_wp + massiso(42,7) = 99.907476_wp + + + ! 43 Tc (Technetium) no stable isotopes + + + ! 44 Ru (Ruthenium) niso(44) = 7 @@ -289,6 +440,12 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(44,5) = 100.905582_wp massiso(44,6) = 101.906323_wp massiso(44,7) = 103.905433_wp + + ! 45 Rh (Rhodium) + + niso(45) = 1 + prob(45,1) = 100.0_wp + massiso(45,1) = 102.905504_wp ! 46 Pd (Palladium) niso(46) = 6 @@ -304,6 +461,43 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(46,4) = 105.903486_wp massiso(46,5) = 107.903892_wp massiso(46,6) = 109.905153_wp + + + ! 47 Ag (Silver) + + niso(47) = 2 + prob(47,1) = 51.839_wp + prob(47,2) = 48.161_wp + massiso(47,1) = 106.905093_wp + massiso(47,2) = 108.904756_wp + + ! 48 Cd (Cadmium) + + niso(48) = 8 + prob(48,1) = 1.25_wp + prob(48,2) = 0.89_wp + prob(48,3) = 12.49_wp + prob(48,4) = 12.80_wp + prob(48,5) = 24.13_wp + prob(48,6) = 12.22_wp + prob(48,7) = 28.73_wp + prob(48,8) = 7.49_wp + massiso(48,1) = 105.906458_wp + massiso(48,2) = 107.904183_wp + massiso(48,3) = 109.903006_wp + massiso(48,4) = 110.904182_wp + massiso(48,5) = 111.9027577_wp + massiso(48,6) = 112.9044014_wp + massiso(48,7) = 113.9033586_wp + massiso(48,8) = 115.904756_wp + + + ! 49 In (Indium) + niso(49) = 2 + prob(49,1) = 4.29_wp + prob(49,2) = 95.71_wp + massiso(49,1) = 112.904062_wp + massiso(49,2) = 114.903879_wp ! 50 Sn (Tin) niso(50) = 10 @@ -359,6 +553,50 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & prob(53,1) = 100.00_wp massiso(53,1) = 126.904473_wp + + ! 54 Xe (Xenon) + niso(54) = 9 + prob(54,1) = 0.0952_wp + prob(54,2) = 0.0890_wp + prob(54,3) = 1.9102_wp + prob(54,4) = 26.4006_wp + prob(54,5) = 4.0710_wp + prob(54,6) = 21.2324_wp + prob(54,7) = 26.9086_wp + prob(54,8) = 10.4357_wp + prob(54,9) = 8.8573_wp + massiso(54,1) = 123.9058954_wp + massiso(54,2) = 125.904268_wp + massiso(54,3) = 127.9035305_wp + massiso(54,4) = 128.9047799_wp + massiso(54,5) = 129.9035089_wp + massiso(54,6) = 130.9050828_wp + massiso(54,7) = 131.9041546_wp + massiso(54,8) = 133.9053945_wp + massiso(54,9) = 135.907220_wp + + ! 55 Cs (Caesium) + niso(55) = 1 + prob(55,1) = 100.0_wp + massiso(55,1) = 132.905447_wp + + ! 56 Ba (Barium) + niso(56) = 7 + prob(56,1) = 0.106_wp + prob(56,2) = 0.101_wp + prob(56,3) = 2.417_wp + prob(56,4) = 6.592_wp + prob(56,5) = 7.854_wp + prob(56,6) = 11.232_wp + prob(56,7) = 71.698_wp + massiso(56,1) = 129.906311_wp + massiso(56,2) = 131.905056_wp + massiso(56,3) = 133.904504_wp + massiso(56,4) = 134.905684_wp + massiso(56,5) = 135.904571_wp + massiso(56,6) = 136.905822_wp + massiso(56,7) = 137.905242_wp + ! 57 La (Lanthanum) niso(57) = 2 prob(57,1) = 0.09_wp @@ -366,6 +604,178 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(57,1) = 137.907112_wp massiso(57,2) = 138.909477_wp + + ! 58 Ce (Cerium) + niso(58) = 4 + prob(58,1) = 0.185_wp + prob(58,2) = 0.251_wp + prob(58,3) = 88.450_wp + prob(58,4) = 11.114_wp + massiso(58,1) = 135.907140_wp + massiso(58,2) = 137.905986_wp + massiso(58,3) = 139.905435_wp + massiso(58,4) = 141.909241_wp + + ! 59 Pr (Praseodymium) + niso(59) = 1 + prob(59,1) = 100.0_wp + massiso(59,1) = 140.907648_wp + + ! 60 Nd (Neodymium) + + niso(60) = 7 + prob(60,1) = 27.2_wp + prob(60,2) = 12.2_wp + prob(60,3) = 23.8_wp + prob(60,4) = 8.3_wp + prob(60,5) = 17.2_wp + prob(60,6) = 5.7_wp + prob(60,7) = 5.6_wp + massiso(60,1) = 141.907719_wp + massiso(60,2) = 142.909810_wp + massiso(60,3) = 143.910083_wp + massiso(60,4) = 144.912569_wp + massiso(60,5) = 145.913113_wp + massiso(60,6) = 147.916889_wp + massiso(60,7) = 149.920887_wp + + ! 61 Pm (Promethium) no stable isotopes + + ! 62 Sm (Samarium) + niso(62) = 7 + prob(62,1) = 3.07_wp + prob(62,2) = 14.99_wp + prob(62,3) = 11.24_wp + prob(62,4) = 13.82_wp + prob(62,5) = 7.38_wp + prob(62,6) = 26.75_wp + prob(62,7) = 22.75_wp + massiso(62,1) = 143.911996_wp + massiso(62,2) = 146.914894_wp + massiso(62,3) = 147.914818_wp + massiso(62,4) = 148.917180_wp + massiso(62,5) = 149.917272_wp + massiso(62,6) = 151.919729_wp + massiso(62,7) = 153.922206_wp + + ! 63 Eu (Europium) + niso(63) = 2 + prob(63,1) = 47.81_wp + prob(63,2) = 52.19_wp + massiso(63,1) = 150.919846_wp + massiso(63,2) = 152.921227_wp + + ! 64 Gd (Gadolinium) + niso(64) = 7 + prob(64,1) = 0.20_wp + prob(64,2) = 2.18_wp + prob(64,3) = 14.80_wp + prob(64,4) = 20.47_wp + prob(64,5) = 15.65_wp + prob(64,6) = 24.84_wp + prob(64,7) = 21.86_wp + massiso(64,1) = 151.919789_wp + massiso(64,2) = 153.920862_wp + massiso(64,3) = 154.922619_wp + massiso(64,4) = 155.922120_wp + massiso(64,5) = 156.923957_wp + massiso(64,6) = 157.924101_wp + massiso(64,7) = 159.927051_wp + + ! 65 Tb (Terbium) + niso(65) = 1 + prob(65,1) = 100.0_wp + massiso(65,1) = 158.925343_wp + + + ! 66 Dy (Dysprosium) + + niso(66) = 7 + prob(66,1) = 0.056_wp + prob(66,2) = 0.095_wp + prob(66,3) = 2.34_wp + prob(66,4) = 18.889_wp + prob(66,5) = 25.475_wp + prob(66,6) = 24.896_wp + prob(66,7) = 28.260_wp + massiso(66,1) = 155.924278_wp + massiso(66,2) = 157.924405_wp + massiso(66,3) = 159.925194_wp + massiso(66,4) = 160.926930_wp + massiso(66,5) = 161.926795_wp + massiso(66,6) = 162.928728_wp + massiso(66,7) = 163.929171_wp + + ! 67 Ho (Holmium) + niso(67) = 1 + prob(67,1) = 100.0_wp + massiso(67,1) = 164.930319_wp + + ! 68 Er (Erbium) + niso(68) = 6 + prob(68,1) = 0.139_wp + prob(68,2) = 1.601_wp + prob(68,3) = 33.503_wp + prob(68,4) = 22.869_wp + prob(68,5) = 26.978_wp + prob(68,6) = 14.910_wp + massiso(68,1) = 161.928775_wp + massiso(68,2) = 163.929197_wp + massiso(68,3) = 165.930290_wp + massiso(68,4) = 166.932046_wp + massiso(68,5) = 167.932368_wp + massiso(68,6) = 169.935461_wp + ! 69 Tm (Thulium) + niso(69) = 1 + prob(69,1) = 100.0_wp + massiso(69,1) = 168.934211_wp + + ! 70 Yb (Ytterbium) + niso(70) = 7 + prob(70,1) = 0.13_wp + prob(70,2) = 3.04_wp + prob(70,3) = 14.28_wp + prob(70,4) = 21.83_wp + prob(70,5) = 16.13_wp + prob(70,6) = 31.83_wp + prob(70,7) = 12.76_wp + massiso(70,1) = 167.933895_wp + massiso(70,2) = 169.934759_wp + massiso(70,3) = 170.936323_wp + massiso(70,4) = 171.936378_wp + massiso(70,5) = 172.938207_wp + massiso(70,6) = 173.938858_wp + massiso(70,7) = 175.942569_wp + + ! 71 Lu (Lutetium) + niso(71) = 2 + prob(71,1) = 97.41_wp + prob(71,2) = 2.59_wp + massiso(71,1) = 174.9407682_wp + massiso(71,2) = 175.9426827_wp + + ! 72 Hf (Hafnium) + niso(72) = 6 + prob(72,1) = 0.16_wp + prob(72,2) = 5.26_wp + prob(72,3) = 18.60_wp + prob(72,4) = 27.28_wp + prob(72,5) = 13.62_wp + prob(72,6) = 35.08_wp + massiso(72,1) = 173.940042_wp + massiso(72,2) = 175.941403_wp + massiso(72,3) = 176.943204_wp + massiso(72,4) = 177.9436981_wp + massiso(72,5) = 178.9458154_wp + massiso(72,6) = 179.9465488_wp + + ! 73 Ta (Tantalum) + niso(73) = 2 + prob(73,1) = 0.012_wp + prob(73,2) = 99.988_wp + massiso(73,1) = 179.947466_wp + massiso(73,2) = 180.947996_wp + ! 74 W (Tungsten) niso(74) = 5 prob(74,1) = 0.12_wp @@ -378,6 +788,37 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(74,3) = 182.950223_wp massiso(74,4) = 183.950931_wp massiso(74,5) = 185.954364_wp + + ! 75 Re (Rhenium) + niso(75) = 2 + prob(75,1) = 37.40_wp + prob(75,2) = 62.60_wp + massiso(75,1) = 184.952955_wp + massiso(75,2) = 186.955750_wp + + ! 76 Os (Osmium) + niso(76) = 7 + prob(76,1) = 0.02_wp + prob(76,2) = 1.59_wp + prob(76,3) = 1.96_wp + prob(76,4) = 13.24_wp + prob(76,5) = 16.15_wp + prob(76,6) = 26.26_wp + prob(76,7) = 40.78_wp + massiso(76,1) = 183.952491_wp + massiso(76,2) = 185.953838_wp + massiso(76,3) = 186.9557476_wp + massiso(76,4) = 187.9558357_wp + massiso(76,5) = 188.958145_wp + massiso(76,6) = 189.958445_wp + massiso(76,7) = 191.961479_wp + + ! 77 Ir (Iridium) + niso(77) = 2 + prob(77,1) = 37.3_wp + prob(77,2) = 62.7_wp + massiso(77,1) = 190.960591_wp + massiso(77,2) = 192.96293_wp ! 78 Pt (Platinum) niso(78) = 5 @@ -391,6 +832,11 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & massiso(78,3) = 194.964791_wp massiso(78,4) = 195.964951_wp massiso(78,5) = 197.967893_wp + + ! 79 Au (Gold) + niso(79) = 1 + prob(79,1) = 100.0_wp + massiso(79,1) = 196.966551_wp ! 80 Hg (Mercury) niso(80) = 7 @@ -431,6 +877,12 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & niso(83) = 1 prob(83,1) = 100.0_wp massiso(83,1) = 208.980398_wp + + ! 84 Po (Polonium) no stable isotopes + + ! 85 At (Astatine) no stable isotopes + + ! 86 Rn (Radon) no stable isotopes ! postprocessing starts here @@ -449,7 +901,7 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & do j = 1, niso(i) sum_prob = sum_prob + prob(i,j) enddo - if ( niso(i) > 0 .and. abs(sum_prob - 1.0_wp) > 0.01_wp ) stop 'internal isotope error 1' + if ( niso(i) > 0 .and. abs(sum_prob - 1.0_wp) > 0.01_wp ) stop 'internal isotope error 1, wrong isotopic abundances in code' enddo endif @@ -466,7 +918,10 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & !> if number of isotopes == 0 <- often errors if wrong .res file is read (is (hopefully) fixed) do i=1,ntot - if(niso(iat_save(i)) == 0) stop 'internal isotope error 2' + if(niso(iat_save(i)) == 0) then + write(*,*) 'internal isotope error 2. no isotopes implemented for element', iat_save(i) + stop + end if enddo diff --git a/src/main.F90 b/src/main.F90 index 6d4a5c4..a094a10 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -1114,10 +1114,6 @@ program plotms ! exp_entries, exp_mass, exp_int,score,tmax) call match(added_thr_masses, added_thr_intensities, new_length_thr, & new_length_exp, added_exp_masses, added_exp_intensities,score,tmax) - write(*,*) "added_thr_masses, added_thr_intensities, new_length_thr, & - new_length_exp, added_exp_masses, added_exp_intensities,score,tmax" - write(*,*) added_thr_masses, added_thr_intensities, new_length_thr, & - new_length_exp, added_exp_masses, added_exp_intensities,score,tmax write(*,*) write(*,*)"!!!!!!!!!!!!!!!!!!!!!!! " write(*,*)" Matching score: " @@ -1187,7 +1183,8 @@ subroutine match(added_masses, added_intensities, new_length, & allocate(w_exp(exp_entries)) norm = 0.0_wp do i = 1, exp_entries - w_exp(i) = exp_mass(i)**2 * exp_int(i)**0.6_wp + ! w_exp(i) = exp_mass(i)**2 * exp_int(i)**0.6_wp + w_exp(i) = exp_mass(i)**3 * exp_int(i)**0.6_wp norm = norm + w_exp(i)**2 enddo @@ -1201,7 +1198,8 @@ subroutine match(added_masses, added_intensities, new_length, & allocate(w_thr(new_length)) norm = 0.0_wp do j = 1, new_length - w_thr(j) = (added_masses(j))**2 * ( added_intensities(j) )**0.6_wp + !w_thr(j) = (added_masses(j))**2 * ( added_intensities(j) )**0.6_wp + w_thr(j) = (added_masses(j))**3 * ( added_intensities(j) )**0.6_wp norm = norm + w_thr(j)**2 enddo From c12b38a315c7a8c96121d6771cde5549b59202f6 Mon Sep 17 00:00:00 2001 From: gorges <58849467+gorges97@users.noreply.github.com> Date: Tue, 26 Sep 2023 13:31:16 +0200 Subject: [PATCH 04/11] upgraded mctc-lib to pass test suite --- subprojects/json-fortran-8.2.5.wrap | 2 ++ subprojects/mctc-lib.wrap | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) create mode 100644 subprojects/json-fortran-8.2.5.wrap diff --git a/subprojects/json-fortran-8.2.5.wrap b/subprojects/json-fortran-8.2.5.wrap new file mode 100644 index 0000000..f3675a6 --- /dev/null +++ b/subprojects/json-fortran-8.2.5.wrap @@ -0,0 +1,2 @@ +[wrap-redirect] +filename = mctc-lib/subprojects/json-fortran-8.2.5.wrap diff --git a/subprojects/mctc-lib.wrap b/subprojects/mctc-lib.wrap index 6b71458..f3f6e67 100644 --- a/subprojects/mctc-lib.wrap +++ b/subprojects/mctc-lib.wrap @@ -1,4 +1,4 @@ [wrap-git] directory = mctc-lib url = https://github.com/grimme-lab/mctc-lib -revision = v0.2.3 +revision = v0.3.1 From 421d777beb9c7383409d8252bb0ba5e61336ced0 Mon Sep 17 00:00:00 2001 From: gorges <58849467+gorges97@users.noreply.github.com> Date: Tue, 26 Sep 2023 14:29:33 +0200 Subject: [PATCH 05/11] deleted subprojects wrap because of compiler errors --- subprojects/json-fortran-8.2.5.wrap | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 subprojects/json-fortran-8.2.5.wrap diff --git a/subprojects/json-fortran-8.2.5.wrap b/subprojects/json-fortran-8.2.5.wrap deleted file mode 100644 index f3675a6..0000000 --- a/subprojects/json-fortran-8.2.5.wrap +++ /dev/null @@ -1,2 +0,0 @@ -[wrap-redirect] -filename = mctc-lib/subprojects/json-fortran-8.2.5.wrap From 4a9d8b6847b5a473783a1a73949322bf122f0f7e Mon Sep 17 00:00:00 2001 From: gorges <58849467+gorges97@users.noreply.github.com> Date: Tue, 26 Sep 2023 15:22:02 +0200 Subject: [PATCH 06/11] update to version 6.2.0 --- src/main.F90 | 34 ++++++++++++++++++++++++++++++++++ src/version.f90 | 6 +++--- 2 files changed, 37 insertions(+), 3 deletions(-) diff --git a/src/main.F90 b/src/main.F90 index a094a10..ae9cf1a 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -236,6 +236,14 @@ program plotms call readl(arg,xx,nn) total_charge = real(xx(1),wp) + ! print version of program + case("--version") + call print_version + stop + + ! print help + case("-h","--help") + call print_help case default error stop ' -- Unrecognized Keyword -- ' @@ -1331,6 +1339,32 @@ subroutine calcfr(pp,pair,sum4) end subroutine calcfr +subroutine print_help() +implicit none + + write(*,*) + write(*,'(1x, ''usage :'')') + write(*,'(1x, ''plotms [options]'')') + write(*,*) + write(*,'(1x, ''plotms needs a qcxms.res qcxms_cid.res tmpqcxms.res or tmpqcxms_cid.res file from a qcxms calculation in the current directory.'')') + write(*,*) + write(*,'(/,1x,''command line arguments'')') + write(*,'(5x,''-a: count charges from -1000 (cthr) to cthr2'')') + write(*,'(5x,''-v: print spectra "mass % intensity counts Int. exptl" to stdout;'')') + write(*,'(5x,'' with "Int. exptl" (experimental) taken from exp.dat but not all'')') + write(*,'(5x,'' exp peaks are exported if no theoretical counterpart exists'')') + write(*,'(5x,''-f : give name of res file'')') + write(*,'(5x,''-t: couting ions with charge x to y (give the value, e.g. "-t 1 2" for charge 1 to 2)'')') + write(*,'(5x,''-w [real]: broadening the charges by an SD '')') + write(*,'(5x,''-s: Take only secondary and tertiary fragmentations (give the value, e.g. "-s 2" for secondary)'')') + write(*,'(5x,''-m [int] set minimum value for m/z, so 100% value will be calc. for higher values (CID)'')') + write(*,'(5x,''-p calculate NO isotope pattern, highest peak of isotope pattern is used'')') + write(*,'(5x,''-i [real]: set minimum intensity that is considered in rel. intensity percent'')') + write(*,'(5x,''-e: provide exp. input file (in .csv format)'')') + write(*,'(5x,''--version: print version of this program'')') + stop ' [-h] displayed. exit.' +end subroutine print_help + ! Get the absolute path to file ! Note, this function replaces shell tilde ~/ with the explicit home dir string. function fullpath(fname) diff --git a/src/version.f90 b/src/version.f90 index 5993a91..207bec5 100644 --- a/src/version.f90 +++ b/src/version.f90 @@ -10,13 +10,13 @@ subroutine print_version write(*,'(6x,''* P l o t M S *'')') write(*,'(6x,''* * * * * * * * * * * * * * * * * *'')') write(*,'(6x,''* QCxMS spectra plotting tool *'')') - write(*,'(6x,''* - v. 6.1.0 - *'')') - write(*,'(6x,''* 25. Jul 2022 *'')') + write(*,'(6x,''* - v. 6.2.0 - *'')') + write(*,'(6x,''* 26. Sep 2023 *'')') write(*,'(6x,''* * * * *'')') write(*,'(6x,''* S. Grimme *'')') write(*,'(6x,''* J. Koopman *'')') write(*,'(6x,''* * * * * * * * * * * * * * * * * *'')') - write(*,'(6x,'' Contributor: C.Bauer, T.Kind '')') + write(*,'(6x,'' Contributor: C.Bauer, J. Gorges, T.Kind '')') write(*,*) end subroutine print_version From 965ae9e164af29d7f4687402f28ffd6df1e6a07e Mon Sep 17 00:00:00 2001 From: Johannes Gorges <58849467+gorges97@users.noreply.github.com> Date: Tue, 20 May 2025 10:30:30 +0200 Subject: [PATCH 07/11] resolved bug when using -p and -m flag --- src/entries.f90 | 7 +- src/entries.i90 | 95 ++ src/isotopes.f90 | 9 +- src/isotopes.i90 | 1022 ++++++++++++++++++++ src/main.F90 | 18 +- src/main.i90 | 1392 +++++++++++++++++++++++++++ subprojects/json-fortran-8.2.5.wrap | 2 + 7 files changed, 2536 insertions(+), 9 deletions(-) create mode 100644 src/entries.i90 create mode 100644 src/isotopes.i90 create mode 100644 src/main.i90 create mode 100644 subprojects/json-fortran-8.2.5.wrap diff --git a/src/entries.f90 b/src/entries.f90 index 7dd1927..e050c1e 100644 --- a/src/entries.f90 +++ b/src/entries.f90 @@ -5,7 +5,7 @@ module count_entries contains subroutine check_entries(index_mass, isotope_masses, exact_intensity, & - list_masses, intensity, count_mass, chrg) + list_masses, intensity, count_mass, chrg, int_masses) integer :: loop, loop2, index_mass integer :: count_mass @@ -31,6 +31,7 @@ subroutine check_entries(index_mass, isotope_masses, exact_intensity, & !real (wp) :: save_list(1000) logical :: there = .true. + logical :: int_masses xmass = 0 loop = 0 @@ -47,11 +48,13 @@ subroutine check_entries(index_mass, isotope_masses, exact_intensity, & !> loop over the list from the isotope subroutine outer: do loop2 = 1, index_mass loop = 0 - + if (int_masses) isotope_masses(loop2) = int(isotope_masses(loop2)) !> loop over all entries of new (check) list inner: do loop = loop + 1 + if (int_masses) list_masses(loop) = int(list_masses(loop)) mass_diff = abs(list_masses(loop) - isotope_masses(loop2)) + if ( list_masses(loop) == 0.0_wp ) then there = .false. diff --git a/src/entries.i90 b/src/entries.i90 new file mode 100644 index 0000000..d3585b6 --- /dev/null +++ b/src/entries.i90 @@ -0,0 +1,95 @@ +# 1 "/tmp1/coding/PlotMS/src/entries.f90" +module count_entries + use xtb_mctc_accuracy, only: wp + implicit none + + contains + + subroutine check_entries(index_mass, isotope_masses, exact_intensity, & + list_masses, intensity, count_mass, chrg, int_masses) + + integer :: loop, loop2, index_mass + integer :: count_mass + integer :: sum_index, i + integer :: check + +!real(wp) :: intensity(index_mass) + real(wp) :: xmass + real(wp) :: isotope_masses(index_mass) + real(wp) :: exact_intensity(index_mass) + real(wp) :: chrg + + real(wp) :: list_masses(10000) + real(wp) :: intensity(10000) + + real(wp) :: mass_diff + +!real(wp), allocatable :: list_masses(:) +!real(wp), allocatable :: save_list(:) +!real(wp), allocatable :: save_int(:) +!real(wp), allocatable :: intensity(:) +!real (wp) :: save_list(1000) +!real (wp) :: save_list(1000) + + logical :: there = .true. + logical :: int_masses + + xmass = 0 + loop = 0 + loop2 = 0 + check = 0 +!count_mass = count_mass + index_mass + if(count_mass == 0)then + sum_index = index_mass + else + sum_index = count_mass + index_mass + endif + + +!> loop over the list from the isotope subroutine +outer: do loop2 = 1, index_mass + loop = 0 + if (int_masses) isotope_masses(loop2) = int(isotope_masses(loop2)) +!> loop over all entries of new (check) list +inner: do + loop = loop + 1 + if (int_masses) list_masses(loop) = int(list_masses(loop)) + mass_diff = abs(list_masses(loop) - isotope_masses(loop2)) + + + if ( list_masses(loop) == 0.0_wp ) then + there = .false. +!write(*,*) 'NULL' + exit inner + +!elseif ( list_masses(loop) == isotope_masses(loop2) ) then + elseif ( mass_diff < 1.0d0-10 .or. mass_diff == 0.0_wp) then + there = .true. + intensity(loop) = intensity(loop) + 1 * exact_intensity(loop2) & + * abs(chrg) + if (loop2 < index_mass ) cycle outer + if (loop2 == index_mass ) exit inner + + +!>> false if not in list, store +!elseif ( list_masses(loop) /= isotope_masses(loop2) ) then + elseif ( mass_diff > 1.0d0-10 ) then + there = .false. + if (loop == sum_index ) exit inner + endif + enddo inner + + + if ( .not. there ) then + count_mass = count_mass + 1 + list_masses(count_mass) = isotope_masses(loop2) + intensity(count_mass) = intensity(count_mass) + 1 * exact_intensity(loop2) & + * abs(chrg) + endif + + enddo outer + + + end subroutine check_entries + +end module count_entries diff --git a/src/isotopes.f90 b/src/isotopes.f90 index 692c25e..877d5ee 100644 --- a/src/isotopes.f90 +++ b/src/isotopes.f90 @@ -26,7 +26,7 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & integer :: loop, index_mass ! integer :: tmp_intensity integer :: store_int(1000) - integer :: mzmin + real(wp) :: mzmin integer :: z_chrg real(wp) :: rnd(nrnd,maxatm) @@ -984,6 +984,11 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & endif + + enddo + write(*,*) 'number of peaks', index_mass + do i = 1, index_mass + write(*,*) list_masses(i), store_int(i) enddo @@ -1003,7 +1008,7 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & ! if no isotopes we take here the peak of the isotope pattern with the highest intensity - if ( no_isotopes ) then + if ( no_isotopes .and. index_mass .gt. 0 ) then iipmax = maxval(exact_intensity) indexipmax = maxloc(exact_intensity, dim = 1) mipmax = isotope_masses(indexipmax) diff --git a/src/isotopes.i90 b/src/isotopes.i90 new file mode 100644 index 0000000..44f3faa --- /dev/null +++ b/src/isotopes.i90 @@ -0,0 +1,1022 @@ +# 1 "/tmp1/coding/PlotMS/src/isotopes.f90" +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +! compute the isotopic mass distribution for +! given chemical formula (ntot atoms with OZ iat_save()) +! returns nsig masses in mass() with probability +! this is a quick and dirty MC algorithm and its run-time depends +! critically on the number of trials nrnd + +! TK See: ATOMIC WEIGHTS OF THE ELEMENTS: REVIEW 2000 (IUPAC Technical Report) +! TK Pure Appl. Chem., Vol. 75, No. 6, pp. 683–800, 2003. + +module isotope_pattern + use xtb_mctc_accuracy, only: wp + implicit none + + contains + +subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & + no_isotopes, index_mass, exact_intensity, isotope_masses, z_chrg) + + integer :: ntot,iat_save(*),maxatm + integer :: niso(200) + integer :: n,i,j,iso,iti + integer :: counter + integer,parameter :: nrnd = 50000 + integer :: loop, index_mass +! integer :: tmp_intensity + integer :: store_int(1000) + integer :: mzmin + integer :: z_chrg + + real(wp) :: rnd(nrnd,maxatm) + real(wp) :: r,sum_prob + real(wp) :: prob(200,10),massiso(200,10),p1,p2,x,xmass + real(wp) :: list_masses(nrnd) + real(wp),allocatable :: isotope_masses(:) + + real(wp) :: mipmax, iipmax ! mass and intensity of highest isotope peak + integer :: indexipmax ! index of peak with highest intensity + real(wp) :: current_mass + + real(wp), allocatable :: exact_intensity(:) +!real(wp) :: exact_intensity(1000) + + logical :: no_isotopes + logical :: there + + niso=0 + prob=0 + massiso=0 + + +! 1 H (Hydrogen) + niso(1) = 2 + prob(1,1) = 0.0115_wp + prob(1,2) = 99.9885_wp + massiso(1,1) = 2.014101_wp + massiso(1,2) = 1.007825_wp + +! 2 He (Helium) + niso(2) = 2 + prob(2,1) = 0.000134_wp + prob(2,2) = 99.99866_wp + massiso(2,1) = 3.016029309_wp + massiso(2,2) = 4.002603249_wp + +! 3 Li (Lithium) + + niso(3) = 2 + prob(3,1) = 7.59_wp + prob(3,2) = 92.41_wp + massiso(3,1) = 6.0151223_wp + massiso(3,2) = 7.0160041_wp + +! 4 Be (Beryllium) + + niso(4) = 1 + prob(4,1) = 100.0_wp + massiso(4,1) = 9.0121822_wp + +! 5 B (Boron) + niso(5) = 2 + prob(5,1) = 19.9_wp + prob(5,2) = 80.1_wp + massiso(5,1) = 10.0129370_wp + massiso(5,2) = 11.0093055_wp + +! 6 C (Carbon) + niso(6) = 2 + prob(6,1) = 1.07_wp + prob(6,2) = 98.93_wp + massiso(6,1) = 13.003354_wp + massiso(6,2) = 12.00_wp + +! 7 N (Nitrogen) + niso(7) = 2 + prob(7,1) = 0.368_wp + prob(7,2) = 99.632_wp + massiso(7,1) = 15.000108_wp + massiso(7,2) = 14.003074_wp + +! 8 O (Oxygen) + niso(8) = 3 + prob(8,1) = 0.038_wp + prob(8,2) = 0.205_wp + prob(8,3) = 99.757_wp + massiso(8,1) = 16.999131_wp + massiso(8,2) = 17.999159_wp + massiso(8,3) = 15.994914_wp + +! 9 F (Fluorine) + niso(9) = 1 + prob(9,1) = 100.0_wp + massiso(9,1) = 18.998403_wp + +! 10 Ne (Neon) + niso(10) = 3 + prob(10,1) = 90.48_wp + prob(10,2) = 0.27_wp + prob(10,3) = 9.25_wp + massiso(10,1) = 19.992440176_wp + massiso(10,2) = 20.99384674_wp + massiso(10,3) = 21.99138550_wp +! 11 Na (Sodium) + niso(11) = 1 + prob(11,1) = 100.0_wp + massiso(11,1) = 22.98976966_wp + +! 12 Mg (Magnesium) + niso(12) = 3 + prob(12,1) = 78.99_wp + prob(12,2) = 10.00_wp + prob(12,3) = 11.01_wp + massiso(12,1) = 23.98504187_wp + massiso(12,2) = 24.98583700_wp + massiso(12,3) = 25.98259300_wp + + +! 13 Al (Aluminium) + niso(13) = 1 + prob(13,1) = 100.0_wp + massiso(13,1) = 26.981538_wp + +! 14 Si (Silicon) + niso(14) = 3 + prob(14,1) = 92.223_wp + prob(14,2) = 4.685_wp + prob(14,3) = 3.092_wp + massiso(14,1) = 27.976926_wp + massiso(14,2) = 28.976494_wp + massiso(14,3) = 29.973770_wp + +! 15 P (Phosphorus) + niso(15) = 1 + prob(15,1) = 100.0_wp + massiso(15,1) = 30.973761_wp + +! 16 S (Sulfur) + niso(16) = 4 + prob(16,1) = 0.02_wp + prob(16,2) = 0.76_wp + prob(16,3) = 4.29_wp + prob(16,4) = 94.93_wp + massiso(16,1) = 35.967080_wp + massiso(16,2) = 32.971458_wp + massiso(16,3) = 33.967867_wp + massiso(16,4) = 31.972071_wp + +! 17 Cl (Chlorine) + niso(17) = 2 + prob(17,1) = 75.76_wp + prob(17,2) = 24.24_wp + massiso(17,1) = 34.968852_wp + massiso(17,2) = 36.965902_wp + +! 18 Ar (Argon) + niso(18) = 3 + prob(18,1) = 0.3365_wp + prob(18,2) = 0.0632_wp + prob(18,3) = 99.6003_wp + massiso(18,1) = 35.9675462_wp + massiso(18,2) = 37.9627322_wp + massiso(18,3) = 39.96238312_wp + +! 19 K (Potassium) + niso(19) = 3 + prob(19,1) = 93.2581_wp + prob(19,2) = 0.0117_wp + prob(19,3) = 6.7302_wp + massiso(19,1) = 38.9637069_wp + massiso(19,2) = 39.96399867_wp + massiso(19,3) = 40.96182597_wp + +! 20 Ca (Calcium) + niso(20) = 5 + prob(20,1) = 96.94_wp + prob(20,2) = 0.65_wp + prob(20,3) = 0.14_wp + prob(20,4) = 2.09_wp + prob(20,5) = 0.18_wp + massiso(20,1) = 39.962591_wp + massiso(20,2) = 41.958618_wp + massiso(20,3) = 42.958766_wp + massiso(20,4) = 43.955482_wp + massiso(20,5) = 45.953688_wp + +! 21 Sc (Scandium) + niso(21) = 1 + prob(21,1) = 100.0_wp + massiso(21,1) = 44.9559102_wp + +! 22 Ti (Titanium) + niso(22) = 5 + prob(22,1) = 8.0_wp + prob(22,2) = 7.3_wp + prob(22,3) = 73.8_wp + prob(22,4) = 5.5_wp + prob(22,5) = 5.4_wp + massiso(22,1) = 45.952632_wp + massiso(22,2) = 46.951763_wp + massiso(22,3) = 47.947946_wp + massiso(22,4) = 48.947870_wp + massiso(22,5) = 49.944791_wp + +! 23 V (Vanadium) + niso(23) = 2 + prob(23,1) = 0.250_wp + prob(23,2) = 99.750_wp + massiso(23,1) = 49.9471627_wp + massiso(23,2) = 50.943963_wp + +! 24 Cr (Chromium) + niso(24) = 4 + prob(24,1) = 4.35_wp + prob(24,2) = 83.79_wp + prob(24,3) = 9.50_wp + prob(24,4) = 2.37_wp + massiso(24,1) = 49.946044_wp + massiso(24,2) = 51.940507_wp + massiso(24,3) = 52.940649_wp + massiso(24,4) = 53.938880_wp + +! 25 Mn (Manganese) + niso(25) = 1 + prob(25,1) = 100.0_wp + massiso(25,1) = 54.938045_wp + +! 26 Fe (Iron) + niso(26) = 4 + prob(26,1) = 5.845_wp + prob(26,2) = 91.754_wp + prob(26,3) = 2.119_wp + prob(26,4) = 0.2_wp + massiso(26,1) = 53.939_wp + massiso(26,2) = 55.934_wp + massiso(26,3) = 56.935_wp + massiso(26,4) = 57.933_wp + +! 27 Co (Cobalt) + niso(27) = 1 + prob(27,1) = 100.000_wp + massiso(27,1) = 58.933195_wp + +! 28 Ni (Nickel) + niso(28) = 5 + prob(28,1) = 68.08_wp + prob(28,2) = 26.22_wp + prob(28,3) = 1.14_wp + prob(28,4) = 3.63_wp + prob(28,5) = 0.93_wp + massiso(28,1) = 57.935343_wp + massiso(28,2) = 59.930786_wp + massiso(28,3) = 60.931056_wp + massiso(28,4) = 61.928345_wp + massiso(28,5) = 63.927966_wp + +! 29 Cu (Copper) + niso(29) = 2 + prob(29,1) = 69.15_wp + prob(29,2) = 30.85_wp + massiso(29,1) = 62.929597_wp + massiso(29,2) = 64.927789_wp + +! 30 Zn (Zinc) + niso(30) = 5 + prob(30,1) = 48.6_wp + prob(30,2) = 27.9_wp + prob(30,3) = 4.1_wp + prob(30,4) = 18.8_wp + prob(30,5) = 0.6_wp + massiso(30,1) = 63.929142_wp + massiso(30,2) = 65.926033_wp + massiso(30,3) = 66.927127_wp + massiso(30,4) = 67.924884_wp + massiso(30,5) = 69.925319_wp + +! 31 Ga (Gallium) + + niso(31) = 2 + prob(31,1) = 60.108_wp + prob(31,2) = 39.892_wp + massiso(31,1) = 68.925581_wp + massiso(31,2) = 70.9247073_wp + +! 32 Ge (Germanium) + niso(32) = 5 + prob(32,1) = 20.38_wp + prob(32,2) = 27.31_wp + prob(32,3) = 7.76_wp + prob(32,4) = 36.72_wp + prob(32,5) = 7.83_wp + massiso(32,1) = 69.924247_wp + massiso(32,2) = 71.922076_wp + massiso(32,3) = 72.923459_wp + massiso(32,4) = 73.921178_wp + massiso(32,5) = 75.921402_wp + +! 33 As (Arsenic) + niso(33) = 1 + prob(33,1) = 100.0_wp + massiso(33,1) = 74.921596_wp + +! 34 Se (Selenium) + niso(34) = 6 + prob(34,1) = 0.89_wp + prob(34,2) = 9.37_wp + prob(34,3) = 7.63_wp + prob(34,4) = 23.77_wp + prob(34,5) = 49.61_wp + prob(34,6) = 8.73_wp + massiso(34,1) = 73.922476_wp + massiso(34,2) = 75.919213_wp + massiso(34,3) = 76.919914_wp + massiso(34,4) = 77.917309_wp + massiso(34,5) = 79.916521_wp + massiso(34,6) = 81.916699_wp + +! 35 Br (Bromine) + niso(35) = 2 + prob(35,1) = 50.69_wp + prob(35,2) = 49.31_wp + massiso(35,1) = 78.91833_wp + massiso(35,2) = 80.91629_wp + +! 36 Kr (Krypton) + niso(36) = 6 + prob(36,1) = 0.355_wp + prob(36,2) = 2.286_wp + prob(36,3) = 11.593_wp + prob(36,4) = 11.500_wp + prob(36,5) = 56.987_wp + prob(36,6) = 17.279_wp + massiso(36,1) = 77.920388_wp + massiso(36,2) = 79.916379_wp + massiso(36,3) = 81.9134850_wp + massiso(36,4) = 82.914137_wp + massiso(36,5) = 83.911508_wp + massiso(36,6) = 85.910615_wp + +! 37 Rb (Rubidium) + niso(37) = 2 + prob(37,1) = 72.17_wp + prob(37,2) = 27.83_wp + massiso(37,1) = 84.9117924_wp + massiso(37,2) = 86.9091858_wp + +! 38 Sr (Strontium) + niso(38) = 4 + prob(38,1) = 0.56_wp + prob(38,2) = 9.86_wp + prob(38,3) = 7.00_wp + prob(38,4) = 82.58_wp + massiso(38,1) = 83.913426_wp + massiso(38,2) = 85.9092647_wp + massiso(38,3) = 86.908816_wp + massiso(38,4) = 87.9056167_wp + +! 39 Y (Yttrium) + + niso(39) = 1 + prob(39,1) = 100.0_wp + massiso(39,1) = 88.9058485_wp + + +! 40 Zr (Zirconium) + niso(40) = 5 + prob(40,1) = 51.45_wp + prob(40,2) = 11.22_wp + prob(40,3) = 17.15_wp + prob(40,4) = 17.38_wp + prob(40,5) = 2.80_wp + massiso(40,1) = 89.904704_wp + massiso(40,2) = 90.905646_wp + massiso(40,3) = 91.905040_wp + massiso(40,4) = 93.906315_wp + massiso(40,5) = 95.908273_wp + + +! 41 Nb (Niobium) + niso(41) = 1 + prob(41,1) = 100.0_wp + massiso(41,1) = 92.9063762_wp + +! 42 Mo (Molybdenum) + niso(42) = 7 + prob(42,1) = 14.77_wp + prob(42,2) = 9.23_wp + prob(42,3) = 15.90_wp + prob(42,4) = 16.68_wp + prob(42,5) = 9.56_wp + prob(42,6) = 24.19_wp + prob(42,7) = 9.67_wp + massiso(42,1) = 91.906810_wp + massiso(42,2) = 93.9050876_wp + massiso(42,3) = 94.9058406_wp + massiso(42,4) = 95.9046780_wp + massiso(42,5) = 96.9060201_wp + massiso(42,6) = 97.9054069_wp + massiso(42,7) = 99.907476_wp + + +! 43 Tc (Technetium) no stable isotopes + + + + +! 44 Ru (Ruthenium) + niso(44) = 7 + prob(44,1) = 5.554_wp + prob(44,2) = 1.873_wp + prob(44,3) = 12.761_wp + prob(44,4) = 12.607_wp + prob(44,5) = 17.062_wp + prob(44,6) = 31.551_wp + prob(44,7) = 18.623_wp + massiso(44,1) = 95.907598_wp + massiso(44,2) = 97.905287_wp + massiso(44,3) = 98.905939_wp + massiso(44,4) = 99.904219_wp + massiso(44,5) = 100.905582_wp + massiso(44,6) = 101.906323_wp + massiso(44,7) = 103.905433_wp + +! 45 Rh (Rhodium) + + niso(45) = 1 + prob(45,1) = 100.0_wp + massiso(45,1) = 102.905504_wp + +! 46 Pd (Palladium) + niso(46) = 6 + prob(46,1) = 1.02_wp + prob(46,2) = 11.15_wp + prob(46,3) = 22.34_wp + prob(46,4) = 27.33_wp + prob(46,5) = 26.47_wp + prob(46,6) = 11.73_wp + massiso(46,1) = 101.905609_wp + massiso(46,2) = 103.904036_wp + massiso(46,3) = 104.905085_wp + massiso(46,4) = 105.903486_wp + massiso(46,5) = 107.903892_wp + massiso(46,6) = 109.905153_wp + + +! 47 Ag (Silver) + + niso(47) = 2 + prob(47,1) = 51.839_wp + prob(47,2) = 48.161_wp + massiso(47,1) = 106.905093_wp + massiso(47,2) = 108.904756_wp + +! 48 Cd (Cadmium) + + niso(48) = 8 + prob(48,1) = 1.25_wp + prob(48,2) = 0.89_wp + prob(48,3) = 12.49_wp + prob(48,4) = 12.80_wp + prob(48,5) = 24.13_wp + prob(48,6) = 12.22_wp + prob(48,7) = 28.73_wp + prob(48,8) = 7.49_wp + massiso(48,1) = 105.906458_wp + massiso(48,2) = 107.904183_wp + massiso(48,3) = 109.903006_wp + massiso(48,4) = 110.904182_wp + massiso(48,5) = 111.9027577_wp + massiso(48,6) = 112.9044014_wp + massiso(48,7) = 113.9033586_wp + massiso(48,8) = 115.904756_wp + + +! 49 In (Indium) + niso(49) = 2 + prob(49,1) = 4.29_wp + prob(49,2) = 95.71_wp + massiso(49,1) = 112.904062_wp + massiso(49,2) = 114.903879_wp + +! 50 Sn (Tin) + niso(50) = 10 + prob(50,1) = 0.97_wp + prob(50,2) = 0.66_wp + prob(50,3) = 0.34_wp + prob(50,4) = 14.54_wp + prob(50,5) = 7.68_wp + prob(50,6) = 24.22_wp + prob(50,7) = 8.59_wp + prob(50,8) = 32.58_wp + prob(50,9) = 4.63_wp + prob(50,10) = 5.79_wp + massiso(50,1) = 111.904848_wp + massiso(50,2) = 113.902779_wp + massiso(50,3) = 114.903342_wp + massiso(50,4) = 115.901741_wp + massiso(50,5) = 116.902952_wp + massiso(50,6) = 117.901603_wp + massiso(50,7) = 118.903308_wp + massiso(50,8) = 119.902194_wp + massiso(50,9) = 121.903439_wp + massiso(50,10) = 123.905274_wp + +! 51 Sb (Antimony) + niso(51) = 2 + prob(51,1) = 57.21_wp + prob(51,2) = 42.79_wp + massiso(51,1) = 120.903816_wp + massiso(51,2) = 122.904214_wp + +! 52 Te (Tellurium) + niso(52) = 8 + prob(52,1) = 0.09_wp + prob(52,2) = 2.55_wp + prob(52,3) = 0.89_wp + prob(52,4) = 4.74_wp + prob(52,5) = 7.07_wp + prob(52,6) = 18.84_wp + prob(52,7) = 31.75_wp + prob(52,8) = 34.09_wp + massiso(52,1) = 119.90402_wp + massiso(52,2) = 121.903044_wp + massiso(52,3) = 122.904270_wp + massiso(52,4) = 123.902817_wp + massiso(52,5) = 124.904431_wp + massiso(52,6) = 125.903312_wp + massiso(52,7) = 127.904463_wp + massiso(52,8) = 129.906224_wp + +! 53 I (Iodine) + niso(53) = 1 + prob(53,1) = 100.00_wp + massiso(53,1) = 126.904473_wp + + +! 54 Xe (Xenon) + niso(54) = 9 + prob(54,1) = 0.0952_wp + prob(54,2) = 0.0890_wp + prob(54,3) = 1.9102_wp + prob(54,4) = 26.4006_wp + prob(54,5) = 4.0710_wp + prob(54,6) = 21.2324_wp + prob(54,7) = 26.9086_wp + prob(54,8) = 10.4357_wp + prob(54,9) = 8.8573_wp + massiso(54,1) = 123.9058954_wp + massiso(54,2) = 125.904268_wp + massiso(54,3) = 127.9035305_wp + massiso(54,4) = 128.9047799_wp + massiso(54,5) = 129.9035089_wp + massiso(54,6) = 130.9050828_wp + massiso(54,7) = 131.9041546_wp + massiso(54,8) = 133.9053945_wp + massiso(54,9) = 135.907220_wp + +! 55 Cs (Caesium) + niso(55) = 1 + prob(55,1) = 100.0_wp + massiso(55,1) = 132.905447_wp + +! 56 Ba (Barium) + niso(56) = 7 + prob(56,1) = 0.106_wp + prob(56,2) = 0.101_wp + prob(56,3) = 2.417_wp + prob(56,4) = 6.592_wp + prob(56,5) = 7.854_wp + prob(56,6) = 11.232_wp + prob(56,7) = 71.698_wp + massiso(56,1) = 129.906311_wp + massiso(56,2) = 131.905056_wp + massiso(56,3) = 133.904504_wp + massiso(56,4) = 134.905684_wp + massiso(56,5) = 135.904571_wp + massiso(56,6) = 136.905822_wp + massiso(56,7) = 137.905242_wp + +! 57 La (Lanthanum) + niso(57) = 2 + prob(57,1) = 0.09_wp + prob(57,2) = 99.91_wp + massiso(57,1) = 137.907112_wp + massiso(57,2) = 138.909477_wp + + +! 58 Ce (Cerium) + niso(58) = 4 + prob(58,1) = 0.185_wp + prob(58,2) = 0.251_wp + prob(58,3) = 88.450_wp + prob(58,4) = 11.114_wp + massiso(58,1) = 135.907140_wp + massiso(58,2) = 137.905986_wp + massiso(58,3) = 139.905435_wp + massiso(58,4) = 141.909241_wp + +! 59 Pr (Praseodymium) + niso(59) = 1 + prob(59,1) = 100.0_wp + massiso(59,1) = 140.907648_wp + +! 60 Nd (Neodymium) + + niso(60) = 7 + prob(60,1) = 27.2_wp + prob(60,2) = 12.2_wp + prob(60,3) = 23.8_wp + prob(60,4) = 8.3_wp + prob(60,5) = 17.2_wp + prob(60,6) = 5.7_wp + prob(60,7) = 5.6_wp + massiso(60,1) = 141.907719_wp + massiso(60,2) = 142.909810_wp + massiso(60,3) = 143.910083_wp + massiso(60,4) = 144.912569_wp + massiso(60,5) = 145.913113_wp + massiso(60,6) = 147.916889_wp + massiso(60,7) = 149.920887_wp + +! 61 Pm (Promethium) no stable isotopes + +! 62 Sm (Samarium) + niso(62) = 7 + prob(62,1) = 3.07_wp + prob(62,2) = 14.99_wp + prob(62,3) = 11.24_wp + prob(62,4) = 13.82_wp + prob(62,5) = 7.38_wp + prob(62,6) = 26.75_wp + prob(62,7) = 22.75_wp + massiso(62,1) = 143.911996_wp + massiso(62,2) = 146.914894_wp + massiso(62,3) = 147.914818_wp + massiso(62,4) = 148.917180_wp + massiso(62,5) = 149.917272_wp + massiso(62,6) = 151.919729_wp + massiso(62,7) = 153.922206_wp + +! 63 Eu (Europium) + niso(63) = 2 + prob(63,1) = 47.81_wp + prob(63,2) = 52.19_wp + massiso(63,1) = 150.919846_wp + massiso(63,2) = 152.921227_wp + +! 64 Gd (Gadolinium) + niso(64) = 7 + prob(64,1) = 0.20_wp + prob(64,2) = 2.18_wp + prob(64,3) = 14.80_wp + prob(64,4) = 20.47_wp + prob(64,5) = 15.65_wp + prob(64,6) = 24.84_wp + prob(64,7) = 21.86_wp + massiso(64,1) = 151.919789_wp + massiso(64,2) = 153.920862_wp + massiso(64,3) = 154.922619_wp + massiso(64,4) = 155.922120_wp + massiso(64,5) = 156.923957_wp + massiso(64,6) = 157.924101_wp + massiso(64,7) = 159.927051_wp + +! 65 Tb (Terbium) + niso(65) = 1 + prob(65,1) = 100.0_wp + massiso(65,1) = 158.925343_wp + + +! 66 Dy (Dysprosium) + + niso(66) = 7 + prob(66,1) = 0.056_wp + prob(66,2) = 0.095_wp + prob(66,3) = 2.34_wp + prob(66,4) = 18.889_wp + prob(66,5) = 25.475_wp + prob(66,6) = 24.896_wp + prob(66,7) = 28.260_wp + massiso(66,1) = 155.924278_wp + massiso(66,2) = 157.924405_wp + massiso(66,3) = 159.925194_wp + massiso(66,4) = 160.926930_wp + massiso(66,5) = 161.926795_wp + massiso(66,6) = 162.928728_wp + massiso(66,7) = 163.929171_wp + +! 67 Ho (Holmium) + niso(67) = 1 + prob(67,1) = 100.0_wp + massiso(67,1) = 164.930319_wp + +! 68 Er (Erbium) + niso(68) = 6 + prob(68,1) = 0.139_wp + prob(68,2) = 1.601_wp + prob(68,3) = 33.503_wp + prob(68,4) = 22.869_wp + prob(68,5) = 26.978_wp + prob(68,6) = 14.910_wp + massiso(68,1) = 161.928775_wp + massiso(68,2) = 163.929197_wp + massiso(68,3) = 165.930290_wp + massiso(68,4) = 166.932046_wp + massiso(68,5) = 167.932368_wp + massiso(68,6) = 169.935461_wp +! 69 Tm (Thulium) + niso(69) = 1 + prob(69,1) = 100.0_wp + massiso(69,1) = 168.934211_wp + +! 70 Yb (Ytterbium) + niso(70) = 7 + prob(70,1) = 0.13_wp + prob(70,2) = 3.04_wp + prob(70,3) = 14.28_wp + prob(70,4) = 21.83_wp + prob(70,5) = 16.13_wp + prob(70,6) = 31.83_wp + prob(70,7) = 12.76_wp + massiso(70,1) = 167.933895_wp + massiso(70,2) = 169.934759_wp + massiso(70,3) = 170.936323_wp + massiso(70,4) = 171.936378_wp + massiso(70,5) = 172.938207_wp + massiso(70,6) = 173.938858_wp + massiso(70,7) = 175.942569_wp + +! 71 Lu (Lutetium) + niso(71) = 2 + prob(71,1) = 97.41_wp + prob(71,2) = 2.59_wp + massiso(71,1) = 174.9407682_wp + massiso(71,2) = 175.9426827_wp + +! 72 Hf (Hafnium) + niso(72) = 6 + prob(72,1) = 0.16_wp + prob(72,2) = 5.26_wp + prob(72,3) = 18.60_wp + prob(72,4) = 27.28_wp + prob(72,5) = 13.62_wp + prob(72,6) = 35.08_wp + massiso(72,1) = 173.940042_wp + massiso(72,2) = 175.941403_wp + massiso(72,3) = 176.943204_wp + massiso(72,4) = 177.9436981_wp + massiso(72,5) = 178.9458154_wp + massiso(72,6) = 179.9465488_wp + +! 73 Ta (Tantalum) + niso(73) = 2 + prob(73,1) = 0.012_wp + prob(73,2) = 99.988_wp + massiso(73,1) = 179.947466_wp + massiso(73,2) = 180.947996_wp + +! 74 W (Tungsten) + niso(74) = 5 + prob(74,1) = 0.12_wp + prob(74,2) = 26.50_wp + prob(74,3) = 14.31_wp + prob(74,4) = 30.64_wp + prob(74,5) = 28.43_wp + massiso(74,1) = 179.946704_wp + massiso(74,2) = 181.948204_wp + massiso(74,3) = 182.950223_wp + massiso(74,4) = 183.950931_wp + massiso(74,5) = 185.954364_wp + +! 75 Re (Rhenium) + niso(75) = 2 + prob(75,1) = 37.40_wp + prob(75,2) = 62.60_wp + massiso(75,1) = 184.952955_wp + massiso(75,2) = 186.955750_wp + +! 76 Os (Osmium) + niso(76) = 7 + prob(76,1) = 0.02_wp + prob(76,2) = 1.59_wp + prob(76,3) = 1.96_wp + prob(76,4) = 13.24_wp + prob(76,5) = 16.15_wp + prob(76,6) = 26.26_wp + prob(76,7) = 40.78_wp + massiso(76,1) = 183.952491_wp + massiso(76,2) = 185.953838_wp + massiso(76,3) = 186.9557476_wp + massiso(76,4) = 187.9558357_wp + massiso(76,5) = 188.958145_wp + massiso(76,6) = 189.958445_wp + massiso(76,7) = 191.961479_wp + +! 77 Ir (Iridium) + niso(77) = 2 + prob(77,1) = 37.3_wp + prob(77,2) = 62.7_wp + massiso(77,1) = 190.960591_wp + massiso(77,2) = 192.96293_wp + +! 78 Pt (Platinum) + niso(78) = 5 + prob(78,1) = 0.782_wp + prob(78,2) = 32.97_wp + prob(78,3) = 33.83_wp + prob(78,4) = 25.24_wp + prob(78,5) = 7.16_wp + massiso(78,1) = 191.961038_wp + massiso(78,2) = 193.962680_wp + massiso(78,3) = 194.964791_wp + massiso(78,4) = 195.964951_wp + massiso(78,5) = 197.967893_wp + +! 79 Au (Gold) + niso(79) = 1 + prob(79,1) = 100.0_wp + massiso(79,1) = 196.966551_wp + +! 80 Hg (Mercury) + niso(80) = 7 + prob(80,1) = 0.15_wp + prob(80,2) = 10.0_wp + prob(80,3) = 16.87_wp + prob(80,4) = 23.10_wp + prob(80,5) = 13.19_wp + prob(80,6) = 29.86_wp + prob(80,7) = 6.87_wp + massiso(80,1) = 195.965833_wp + massiso(80,2) = 197.966769_wp + massiso(80,3) = 198.968279_wp + massiso(80,4) = 199.968326_wp + massiso(80,5) = 200.970302_wp + massiso(80,6) = 201.970643_wp + massiso(80,7) = 203.973494_wp + +! 81 Tl (Thallium) + niso(81) = 2 + prob(81,1) = 29.52_wp + prob(81,2) = 70.48_wp + massiso(81,1) = 202.972344_wp + massiso(81,2) = 204.974427_wp + +! 82 Pb (Lead) + niso(82) = 4 + prob(82,1) = 1.3_wp + prob(82,2) = 24.1_wp + prob(82,3) = 22.1_wp + prob(82,4) = 52.4_wp + massiso(82,1) = 203.973043_wp + massiso(82,2) = 205.974465_wp + massiso(82,3) = 206.975897_wp + massiso(82,4) = 207.976652_wp + +! 83 Bi (Bismuth) + niso(83) = 1 + prob(83,1) = 100.0_wp + massiso(83,1) = 208.980398_wp + +! 84 Po (Polonium) no stable isotopes + +! 85 At (Astatine) no stable isotopes + +! 86 Rn (Radon) no stable isotopes + + +! postprocessing starts here + + prob = prob * 0.01_wp + list_masses = 0.0_wp + index_mass = 0 + loop = 0 + store_int = 0.0_wp + +! Has to be done only once to check for correct isotope entries +!> loop over all mases to element 83 (Bismuth) + if (counter == 1) then + do i = 1, 83 + sum_prob = 0 + do j = 1, niso(i) + sum_prob = sum_prob + prob(i,j) + enddo + if ( niso(i) > 0 .and. abs(sum_prob - 1.0_wp) > 0.01_wp ) stop 'internal isotope error 1, wrong isotopic abundances in code' + enddo + endif + +!> this checks if the fragment has atms larger than 100 (no clue why this is checked) +!write(*,*) 'ntot', ntot +!do i = 1, ntot +! write(*,*) i, iat_save(i) +! if ( iat_save(i) > 100 ) then +! niso (iat_save(i)) = 1 +! prob (iat_save(i),1) = 1.0_wp +! massiso(iat_save(i),1) = iat_save(i) - 100.0_wp +! endif +!enddo + +!> if number of isotopes == 0 <- often errors if wrong .res file is read (is (hopefully) fixed) + do i=1,ntot + if(niso(iat_save(i)) == 0) then + write(*,*) 'internal isotope error 2. no isotopes implemented for element', iat_save(i) + stop + end if + enddo + + + +! Calculate isotope pattern + +! loop over random number runs + do n = 1, nrnd + xmass = 0 + do i = 1, ntot + iti = iat_save(i) + r = rnd(n,i) + p1 = 0.0_wp + p2 = prob(iti,1) + do iso = 1, niso(iti) + if ( r >= p1 .and. r <= p2 ) then + x = massiso(iti,iso) + exit + endif + p1 = p2 + p2 = p2 + prob(iti,iso+1) + enddo + xmass = xmass + x + enddo + + current_mass = xmass / float(abs(z_chrg)) + + there = .true. +!> only values larger than user input + if ( current_mass > mzmin ) then + +!> loop over all entries in the list +inner: do + loop = loop + 1 + +!>> write if there is no entry + if ( list_masses(loop) == 0.0_wp ) then + there = .false. + exit inner +!>> true if already in list, end + elseif ( abs(list_masses(loop) - current_mass) < 1.0d-10 ) then + there = .true. + store_int(loop) = store_int(loop) + 1 + exit inner +!>> false if not in list, store + elseif ( abs(list_masses(loop) - current_mass) > 1.0d-10 ) then + there = .false. + cycle inner + endif + enddo inner + loop = 0 + +!> if it is not in the list, add it + if ( .not. there ) then + index_mass = index_mass + 1 + list_masses(index_mass) = current_mass + store_int(index_mass) = store_int(index_mass) + 1 +!write(*,*) list_masses(index_mass) + endif + endif + + + enddo + + + allocate(exact_intensity(index_mass)) + allocate(isotope_masses(index_mass)) + + + do loop = 1, index_mass + isotope_masses(loop) = list_masses(loop) +!exact_intensity(loop) = real(store_int(loop),wp) / sum(real(store_int,wp)) + exact_intensity(loop) = float(store_int(loop)) / sum(float(store_int)) +!write(*,*) isotope_masses(loop), store_int(loop) +!write(*,*) exact_intensity(loop) + enddo + + + + +! if no isotopes we take here the peak of the isotope pattern with the highest intensity + if ( no_isotopes ) then + iipmax = maxval(exact_intensity) + indexipmax = maxloc(exact_intensity, dim = 1) + mipmax = isotope_masses(indexipmax) + deallocate(exact_intensity,isotope_masses) + index_mass = 1 + allocate(exact_intensity(1),isotope_masses(1)) + exact_intensity(1)=iipmax + isotope_masses(1)=mipmax + endif + + +end subroutine isotope + + +end module isotope_pattern diff --git a/src/main.F90 b/src/main.F90 index ae9cf1a..f5667ac 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -51,7 +51,8 @@ program plotms integer :: iat (10000) integer :: nat (10000) integer :: iat_save (10000) - integer :: isec,jcoll,jsec,ial,jal(0:10),kal(0:30),nf,irun,mzmin + integer :: isec,jcoll,jsec,ial,jal(0:10),kal(0:30),nf,irun + real(wp) :: mzmin integer :: maxrun ! TK number of peaks in in-silico spectra, needed for JCAMP-DX export integer :: io, io_spec, io_raw, io_mass, io_csv, io_exp, io_jcamp @@ -88,6 +89,8 @@ program plotms logical :: expdat, nistdat logical :: ex,ex1,ex2,ex3,ex4 logical :: noIso, Plasma + logical :: int_masses + ! fname= or result file, xname contains the mass.agr plot file character(len=80) :: xname @@ -121,10 +124,11 @@ program plotms iarg = 0 norm = 1.0 nagrfile = 410 - mzmin = 10 + mzmin = 10.0_wp spec = 0 noIso = .false. Plasma = .false. + int_masses = .false. cthr = 1.d-3 chrg_wdth = 0.0_wp @@ -181,6 +185,8 @@ program plotms verbose = .true. case("-p","--noisotope") noIso = .true. + case("--int_masses") + int_masses = .true. case("-t","--unity") iarg = iarg + 1 @@ -219,7 +225,7 @@ program plotms iarg = iarg + 1 call get_argument(iarg, arg) call readl(arg,xx,nn) - mzmin=int(xx(1)) + mzmin=xx(1) ! set minimum intensity that is considered for the output ! default: >1%. For all (>0%): set = 0 @@ -341,7 +347,7 @@ program plotms ! -m if(mzmin > 10)then - write(*,'('' - Only m/z values greater than '',i4,'' are considered'')') mzmin + write(*,'('' - Only m/z values greater than '',f7.4,'' are considered'')') mzmin endif ! -i @@ -353,6 +359,7 @@ program plotms if ( noIso ) then write(*,'('' - No isotope pattern calculated '')') endif + ! -t (choose if unity intensities or normal) !if(cthr >= 0)then @@ -561,7 +568,7 @@ program plotms !> sort the single fragment intensities over the entire list of frags if (index_mass > 0) then call check_entries( index_mass, isotope_masses, exact_intensity, & - list_masses, intensity, count_mass, chrg) + list_masses, intensity, count_mass, chrg,int_masses) endif deallocate(exact_intensity) @@ -1360,6 +1367,7 @@ subroutine print_help() write(*,'(5x,''-m [int] set minimum value for m/z, so 100% value will be calc. for higher values (CID)'')') write(*,'(5x,''-p calculate NO isotope pattern, highest peak of isotope pattern is used'')') write(*,'(5x,''-i [real]: set minimum intensity that is considered in rel. intensity percent'')') + write(*,'(5x,''--int_masses : round m/z values to integers'')') write(*,'(5x,''-e: provide exp. input file (in .csv format)'')') write(*,'(5x,''--version: print version of this program'')') stop ' [-h] displayed. exit.' diff --git a/src/main.i90 b/src/main.i90 new file mode 100644 index 0000000..20fc1f8 --- /dev/null +++ b/src/main.i90 @@ -0,0 +1,1392 @@ +# 1 "/tmp1/coding/PlotMS/src/main.F90" +! ============================================================================ +! Name : plotms.f90 +!! Authors : S. Grimme +!! Co-authors : J.Koopman, C.Bauer +!! Contributions : T. Kind (FiehnLab 2013) +!! Version : 5.1 (Aug 06 2021) +!! Copyright : S. Grimme +!! Description : plot mass spectra from QCxMS +!! ============================================================================ +! +!!==================================================================== +!! Original PlotMS Program for extraction of mass spectra +!! from quantum mechanical simulations used within QCxMS +!! +!! ********************************************* +!! * S. Grimme * +! * Mulliken Center for Theoretical Chemistry * +! * Universitaet Bonn * +! * 2008-22 * +! ********************************************* +! +! Please cite as: +! S.Grimme, Angew.Chem.Int.Ed. 52 (2013) 6306-6312. +! +! Involves: +! - Matching score -- C.Bauer +! - JCAMP-DX Files -- T.Kind +! +!==================================================================== + +! 5.0: Changing name from QCEIMS to QCxMS +! 5.1: Changing from pseudo-random distr. to BoxMuller distribution + + +program plotms + use count_entries, only: check_entries + use isotope_pattern, only: isotope + use version, only: print_version + use qcxms_boxmuller, only: vary_energies + use readcommon + use mctc_env, only : error_type, get_argument + use xtb_mctc_accuracy, only: wp + use xtb_mctc_convert + use xtb_mctc_symbols, only: toSymbol + implicit none + + integer :: n,i,j,k,kk,kkk,kkkk,nn,ntot + integer :: atm_types + integer :: maxatm,imax,imin,nagrfile + integer :: spec + integer :: iat (10000) + integer :: nat (10000) + integer :: iat_save (10000) + integer :: isec,jcoll,jsec,ial,jal(0:10),kal(0:30),nf,irun,mzmin + integer :: maxrun +! TK number of peaks in in-silico spectra, needed for JCAMP-DX export + integer :: io, io_spec, io_raw, io_mass, io_csv, io_exp, io_jcamp + integer :: counter + integer :: z_chrg + integer :: index_mass, count_mass, sum_index + integer :: sorted_index + integer :: list_length + integer :: exp_entries + integer :: iarg, narg + + real(wp) :: xx(100),tmax,r,rms,norm,cthr,cthr2 + real(wp) :: chrg,chrg2,checksum,score +!real(wp) :: dum + real(wp) :: chrg_wdth + real(wp) :: total_charge + real(wp) :: lowest + real(wp) :: min_intensity + real(wp), allocatable :: exp_mass (:) + real(wp), allocatable :: exp_int (:) + real(wp), allocatable :: reduced_intensities(:), reduced_masses(:) + real(wp), allocatable :: sorted_masses(:), sorted_intensities(:) + real(wp), allocatable :: exact_intensity(:) + real(wp), allocatable :: isotope_masses(:) + real(wp), allocatable :: rnd(:,:) + real(wp), allocatable :: checksum2(:) +!real(wp), allocatable :: list_masses(:) +!real(wp), allocatable :: intensity(:) + real(wp) :: list_masses(10000) + real(wp) :: intensity(10000) + real(wp) :: mmax + + logical :: sel,verbose,small + logical :: expdat, nistdat + logical :: ex,ex1,ex2,ex3,ex4 + logical :: noIso, Plasma + logical :: int_masses + + +! fname= or result file, xname contains the mass.agr plot file + character(len=80) :: xname + character(len=:), allocatable :: fname,fname1,fname2,fname3,fname4 + character(len=4096), external :: fullpath + character(len=:), allocatable :: exp_dat + character(len=:), allocatable :: arg + +!!!!!!!!!!!!!!!!!!!!! + integer :: sm + integer, allocatable :: int_exp(:), int_thr(:) + integer :: store_check_list + integer :: new_length_thr, new_length_exp + real(wp) :: store_intensities + real(wp), allocatable :: added_thr_intensities(:), added_exp_intensities(:) + real(wp), allocatable :: added_thr_masses(:), added_exp_masses(:) + real(wp) :: store_thr_int(1000) + real(wp) :: store_exp_int(1000) + real(wp) :: store_thr_mass(1000) + real(wp) :: store_exp_mass(1000) +!real(wp) :: store_check_list +!integer, allocatable :: added_masses(:) +!integer, allocatable :: exp_mass (:) +!integer :: store_mass(1000) +!!!!!!!!!!!!!!!!!!!!! + + expdat = .false. + verbose = .false. + small = .false. + isec = 0 + iarg = 0 + norm = 1.0 + nagrfile = 410 + mzmin = 10 + spec = 0 + noIso = .false. + Plasma = .false. + int_masses = .false. + + cthr = 1.d-3 + chrg_wdth = 0.0_wp + total_charge = 1.0_wp + cthr2 = 0.01 ! only used for information, not in program + min_intensity = 0.0_wp + + z_chrg = 1 + + +! edit this path name to some standard xmgrace plot file +! JK changed to universal location/file + +!xname='~/.mass_raw.agr' +!xname='~/.mass_jay.agr' + fname='' + exp_dat="exp.dat" + +! start loop processing aruments + +! comand line parameters +!> -a count charges from -1000 (cthr) to cthr2 +!> -v print spectra "mass % intensity counts Int. exptl" to stdout; +! with "Int. exptl" (experimental) taken from exp.dat but not all +! exp peaks are exported if no theoretical counterpart exists +!> -f filename or -f +!> -t couting ions with charge x to y (give the value, e.g. "-t 1 2" +! for charge 1 to 2) +!> -w broadening the charges by an SD +!> -s Take only secondary and tertiary fragmentations (give the value, +! e.g. "-s 2" for secondary) +! -m set minimum value for m/z, so 100% value will be calc. for higher values (CID) +! -p calculate NO isotope pattern +! -i set minimum intensity that is considered in rel. intensity percent +! -e provide exp. input file (in .csv format) + + narg = command_argument_count() + + do while(iarg < narg) + iarg = iarg + 1 + call get_argument(iarg, arg) + + select case(arg) +!if(index(arg(i),'-a') /= 0) cthr = -1000.0_wp +!if(index(arg(i),'-v') /= 0) verbose = .true. +!if(index(arg(i),'-f') /= 0) fname = arg(i+1) +!if(index(arg(i),'-p') /= 0) noIso = .true. +!if(index(arg(i),'-e') /= 0) exp_dat = arg(i+1) + + + case("-a") + cthr = -1000.0_wp + case("-v","--verbose") + verbose = .true. + case("-p","--noisotope") + noIso = .true. + case("--int_masses") + int_masses = .true. + + case("-t","--unity") + iarg = iarg + 1 + call get_argument(iarg, arg) + call readl(arg,xx,nn) + cthr = real(xx(1),wp) + +!> -f filename or -f + case("-f","--file") + iarg = iarg + 1 + call get_argument(iarg, arg) + fname = arg + +! -e provide exp. input file (in .csv format) + case("-e","--exp") + iarg = iarg + 1 + call get_argument(iarg, arg) + exp_dat = arg + +! set width distr. of boltz sampling + case("-w","--width") + iarg = iarg + 1 + call get_argument(iarg, arg) + call readl(arg,xx,nn) + chrg_wdth = real(xx(1),wp) + +! set number of cascading runs + case("-s","--cascades") + iarg = iarg + 1 + call get_argument(iarg, arg) + call readl(arg,xx,nn) + isec = int(xx(1)) + +! set minimum mass + case("-m","--min") + iarg = iarg + 1 + call get_argument(iarg, arg) + call readl(arg,xx,nn) + mzmin=int(xx(1)) + +! set minimum intensity that is considered for the output +! default: >1%. For all (>0%): set = 0 + case("-i","--mzmin") + iarg = iarg + 1 + call get_argument(iarg, arg) + call readl(arg,xx,nn) + min_intensity = real(xx(1),wp) + +! get initial charge of system <- not used anymore + case("-c","--charge") + iarg = iarg + 1 + call get_argument(iarg, arg) + call readl(arg,xx,nn) + total_charge = real(xx(1),wp) + +! print version of program + case("--version") + call print_version + stop + +! print help + case("-h","--help") + call print_help + + case default + error stop ' -- Unrecognized Keyword -- ' + end select + enddo + + + xname = '~/.mass_raw.agr' + +! fname contains the results from each calculation or the temporary result tmpqcxms.res +! xname contains the xmgrace plot file + + if ( Plasma ) then + write(*,*) 'Print neutral PLASMA spectrum' + fname = 'qcxms_neutrals.res' + spec = 1 + + elseif ( fname == '' ) then + fname1 = 'qcxms.res' + fname2 = 'qcxms_cid.res' + fname3 = 'tmpqcxms.res' + fname4 = 'tmpqcxms_cid.res' + +! check if either EI or CID res file exist + inquire ( file = fname1, exist = ex1 ) + inquire ( file = fname2, exist = ex2 ) + inquire ( file = fname3, exist = ex3 ) + inquire ( file = fname4, exist = ex4 ) + + if ( ex1 .and. ex2 ) stop 'cannot process both EI and CID .res files!' + if ( ex3 .and. ex4 ) stop 'cannot process both EI and CID tmp.res files!' + +! none exist + if ( .not. ex1 .and. .not. ex2 .and. .not. ex3 .and. .not. ex4 ) stop '! res file does not exist !' + +! if EI exists + if ( ex1 .or. ex3 ) then + spec = 1 + if ( ex1 ) then + fname = 'qcxms.res' + else + fname = 'tmpqcxms.res' + endif + endif + +! if CID exists + if ( ex2 .or. ex4 ) then + spec = 2 + if ( ex2 ) then + fname = 'qcxms_cid.res' + else + fname = 'tmpqcxms_cid.res' + endif + endif +! for -f option + elseif (index(fname,'cid') /= 0) then + spec = 2 ! CID + else + spec = 1 ! EI + endif + + call print_version + + write(*,'(6x,''-> Reading .res file: '',2x,(a))')trim(fname) + +!> write out if compared + if ( exp_dat == 'exp.dat') then + inquire(file='exp.dat',exist=expdat) + elseif (exp_dat /= 'exp.dat') then + inquire(file=exp_dat,exist=expdat) + endif + if (expdat) write(*,'(6x,''-> Experimental file: '',2x, (a))') exp_dat + +!> give info about xmgrace file (if verbose) + if (verbose) write(*,'(6x '' -> xmgrace file body '',(a) )')trim(xname) + write(*,*) + +! ------------------------------------------------------------------------------------------------------! +! execute arguments + + if (narg > 0 ) then + write(*,'(50(''!''))') + write(*,'('' The following settings are being used : '')') + write(*,*) + +! -w + if(chrg_wdth > 0)then + write(*,'('' - Broadening of the charges by an SD, wdth :'',f4.1)')chrg_wdth*100. + endif + +! -s + if(isec /= 0)then + write(*,'('' - Taking only secondary, tertiary ... fragmentations '')') isec + endif + +! -m + if(mzmin > 10)then + write(*,'('' - Only m/z values greater than '',i4,'' are considered'')') mzmin + endif + +! -i + if ( min_intensity > 0.0_wp ) then + write(*,'('' - Minimum signal intensity: '',f7.4,'' % (relative)'')') min_intensity + endif + +! -p + if ( noIso ) then + write(*,'('' - No isotope pattern calculated '')') + endif + + +! -t (choose if unity intensities or normal) +!if(cthr >= 0)then +! write(*,'('' - couting ions with charge from '',f4.1,'' to '',f4.1)') cthr,cthr2 + total_charge +!else +! write(*,'('' - counting all fragments with unity charge (frag. overview)'')') +!endif + write(*,*) + write(*,'(50(''!''))') + endif + +! ------------------------------------------------------------------------------------------------------! + maxrun = 0 + i=1 + maxatm=0 + counter = 0 +! ------------------------------------------------------------------------------------------------------! + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! read .res file for the first time: +! charge(chrg2), trajectory(irun),number of collision event(icoll), numbers of (secondary-)fragments (jsec,nf), +! no. of different atom types in the fragment (atm_types), atomic number (iat), +! amount of this atom typ (nat) <== from kk = 1 to k +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + open ( file=fname, newunit=io_spec, status='old' ) + + do + + counter = counter + 1 +!write(*,*) 'COUNTER', counter + + if (spec == 1) then !EI + read(io_spec,*,iostat=io) chrg2, z_chrg, irun, jsec, nf, atm_types, (iat(kk),nat(kk), kk = 1, atm_types) + if(io<0) exit !EOF + if(io>0) stop ' -- Fail in read -- ' !fail + + elseif(spec == 2) then !CID + read(io_spec,*,iostat=io) chrg2, z_chrg, irun, jcoll, jsec, nf, atm_types, (iat(kk), nat(kk), kk = 1, atm_types) + if(io<0) exit !EOF + if(io>0) stop ' -- Fail in read -- ' !fail + + else + write(*,*) 'S T O P - Something wrong in plotms - no spec' + stop + endif + + + if (isec > 0 .and. isec /= jsec) cycle !check tertiary,etc fragmentation (isec) + + if ( z_chrg > 0 ) then + if (chrg2 > cthr) then !default: chrg2 > 0.001 + ntot = sum(nat(1:atm_types)) + if ( ntot > maxatm ) maxatm = ntot !get highest number of atoms in fragment + if ( z_chrg > int(total_charge) ) total_charge = real( z_chrg, wp ) + endif + elseif ( z_chrg < 0 ) then + if (chrg2 < -1.0_wp*cthr) then !default: chrg2 > 0.001 + ntot = sum(nat(1:atm_types)) + if ( ntot > maxatm ) maxatm = ntot !get highest number of atoms in fragment + if ( z_chrg < int(total_charge) ) total_charge = real( z_chrg, wp ) + endif + endif + + +!if (chrg2 > cthr) then !default: chrg2 > 0.001 +! ntot = sum(nat(1:atm_types)) +! if ( ntot > maxatm ) maxatm = ntot !get highest number of atoms in fragment +!endif + + i = i + 1 !count number of single fragments with charge > chrt + if ( irun > maxrun ) maxrun = irun !save highest run number + + enddo +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!write(*,*) 'The maximum charge of the system is ', total_charge + +!> allocate the variables + allocate (checksum2(maxrun)) + + n = i - 1 ! n = no. of fragments with amount maxatm of atoms + + write(*,*) + write(*,'(''--------------------------------------- '')') + write(*,'(i4, '' fragments with '', i3 ,'' atoms max.'')') n, maxatm + write(*,'(''--------------------------------------- '')') + write(*,*) + +!close file + close(io_spec) + +! initialize the random number array (efficiency) + allocate(rnd(50000,maxatm)) + do i = 1, 50000 + do j = 1, maxatm + call random_number(r) + rnd(i,j) = r + enddo + enddo + + +! read it again + write(*,*) 'Computing ...' + write(*,*) + +! contains qcxms.res as standard option + open ( file=fname, newunit=io_spec, status='old') + imin=100000 + i = 1 + counter = 0 + ial = 0 + jal = 0 + kal = 0 + checksum = 0.0_wp + checksum2 = 0.0_wp + index_mass= 0 + intensity = 0 + count_mass = 0 + list_masses = 0.0_wp + sum_index = 0 + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!> start the loop, processing the qcxms.res file, line-by-line +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + do + + counter = counter + 1 + +!> read the fragment entry from qcxms.res line by line + if (spec == 1) then !EI + read(io_spec,*,iostat=io) chrg2, z_chrg, irun, jsec, nf, atm_types, (iat(kk),nat(kk),kk=1,atm_types) + if(io<0) exit !EOF + if(io>0) stop !fail + + elseif(spec == 2) then !CID + read(io_spec,*,iostat=io) chrg2, z_chrg, irun, jcoll, jsec, nf, atm_types, (iat(kk),nat(kk),kk=1,atm_types) + if(io<0) exit !EOF + if(io>0) stop !fail + + endif + + + sel = .false. + chrg = chrg2 + +!> sum up the charges + checksum2(irun) = checksum2(irun) + chrg2 + +!> manual setting (see above) + if ( cthr < 0 ) chrg = total_charge + +!> count the moment where fragments were created (first/second MD or collision...) +!if ( (z_chrg > 0 .and. chrg > cthr) .or. (z_chrg < 0 .and. chrg < -1* cthr) ) then ! default: yes + if ( abs(chrg) > cthr ) then ! default: yes + sel = .true. + ial = ial + 1 !count total amount of fragmentations + jal(jsec) = jal(jsec) + 1 !count amount of first, sec., tert., ... fragmentations + +!> if not fragmented (i.e. smaller/larger than total_charge) count it to kal(1) +!if ( abs(chrg) < abs(total_charge) ) then +! kal(jcoll) = kal(jcoll) + 1 !count the collisions +!else +! kal(1) = kal(jcoll) + 1 !no fragmentation +!endif + + +!> count the collision at which most framentation occured + if ( spec == 2 ) then + if ( chrg == total_charge ) then + kal(1) = kal(1) + 1 !count the collisions + else + kal(jcoll) = kal(jcoll) + 1 !no fragmentation + endif + endif + + endif +!---- + if ( isec > 0 .and. isec /= jsec ) sel = .false. !Only if asked, not default +!---- + +!> get the specifics of the current fragment entry + if ( sel ) then ! default +! sum all atoms of all types (total atoms) + ntot = sum(nat(1:atm_types)) +! all types of atoms + kkkk = 0 + do kk = 1, atm_types +! all atoms of this type + do kkk = 1, nat(kk) + kkkk = kkkk + 1 + iat_save(kkkk) = iat(kk) + enddo + enddo + + checksum = checksum + chrg !Check total charge + +!> compute isotope patterns via monte carlo and save intensites +! of all possible combinations + call isotope (counter, mzmin, ntot, iat_save, maxatm, rnd, & + noIso, index_mass, exact_intensity, isotope_masses, z_chrg) + +!>> distribute charge (if set as input) + if(chrg_wdth > 1.0d-6) chrg = vary_energies(chrg,chrg_wdth) + +!> sort the single fragment intensities over the entire list of frags + if (index_mass > 0) then + call check_entries( index_mass, isotope_masses, exact_intensity, & + list_masses, intensity, count_mass, chrg,int_masses) + endif + + deallocate(exact_intensity) + deallocate(isotope_masses) + + endif + + enddo +!> end the loop +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + list_length=0 + +!if (z_chrg == 1) extra_mass = z_chrg * amutoau !init for electron mass + +!> find the highest intensity + tmax = maxval(intensity) + +!> normalize to 100% + intensity(:) = intensity(:)/tmax * 100.0_wp + + allocate(reduced_masses(count_mass), reduced_intensities(count_mass)) + + reduced_masses = 0.0_wp + reduced_intensities = 0.0_wp + +!> sort out everything we don't like + do i = 1, count_mass + if ( (intensity(i)) >= min_intensity ) then + list_length = list_length + 1 + reduced_intensities(list_length) = intensity(i) + reduced_masses(list_length) = list_masses(i) + endif + enddo + +!> sort the list with increasing masses + allocate(sorted_masses(list_length), sorted_intensities(list_length)) + lowest = 0.0_wp + +!>> find the lowest entry and save it in new lists, repeat with all values +! larger than the last value (lowest) + do i = 1, list_length + sorted_index = minloc(reduced_masses,1, mask = reduced_masses > lowest) + lowest = reduced_masses(sorted_index) + sorted_masses(i) = reduced_masses(sorted_index) + sorted_intensities(i) = reduced_intensities(sorted_index) +!write(*,*) i, sorted_index, sorted_masses(i), sorted_intensities(i) + enddo + + +!> write out the sum of charges + write(*,*) n,' (charged) fragments done.' + write(*,*) + write(*,*) 'checksum of charge :', checksum + write(*,*) + + k=0 +!> count absolute values larger than 1e-6 + do i = 1, maxrun + if ( abs(checksum2(i) ) > 1.0d-6) k = k + 1 + if ( abs(checksum2(i) ) > 1.0d-6 .and. abs(checksum2(i) - total_charge ) > 1.0d-3)then + write(*,*) 'checksum error for trj', i,' chrg=',checksum2(i) + endif + enddo + + write(*,'(''* * * * * * * * * * * * * * * * * * * * '')') + write(*,*) k,' successfull runs!! ' + write(*,'(''* * * * * * * * * * * * * * * * * * * * '')') + + write(*,*) + write(*,*) 'ion sources:' + write(*,'(''% inital run'',f6.1)') 100.0_wp * float(jal(1)) / float(ial) + + do j=2,8 + write(*,'(''% '',i1,''th'',f6.1)') j ,100.0_wp * float(jal(j))/float(ial) + enddo + + if ( spec == 2 ) then + write(*,*) + write(*,'(''% collisions '',f6.1)') 100.0_wp * float(kal(1)) / float(ial) + + do j=2,16 + write(*,'(''% '',i2,''th'',f6.1)') j ,100.0_wp * float(kal(j))/float(ial) + enddo + endif + + write(*,*) + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! read the provided experimental file (exp.dat) ! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + imax=0 + imin=0 + + inquire(file='nist.dat',exist=nistdat) + + + if(nistdat) then + write(*,*) 'Reading nist.dat ...' + open( file = 'nist.dat', newunit = io_exp, status = 'old') + +!> read the experimental (NIST/JCAMP) file +rd: do + read(io_exp,'(a)',iostat=iocheck)line + if (iocheck < 0) exit rd + +!> read the maximum intensity of exp. + if(index(line,'##MAXY=') /= 0)then + line(7:7)=' ' + call readl(line,xx,nn) + norm=xx(nn)/100.0_wp + endif + +!> allocate the number of entries + if(index(line,'##NPOINTS=') /= 0)then + line(10:10)=' ' + call readl(line,xx,nn) + exp_entries = nint(xx(nn)) +!write(*,*) 'EXP ENTRIES', exp_entries + allocate (exp_mass(exp_entries), & + & exp_int(exp_entries)) + endif + + if(index(line,'##PEAK TABLE') /= 0)then + kk=0 + do + read(io_exp,'(a)',iostat=iocheck)line + if (iocheck < 0) exit rd + +! TK JCAMP DX for MS data has "##END=" + if(index(line,'##END') /= 0) cycle + + do k=1,80 + if(line(k:k) == ',')line(k:k)=' ' + enddo + + call readl(line,xx,nn) + + do k = 1, nn/2 + kk = kk + 1 + exp_mass(kk) = xx(2 * k - 1) + exp_int (kk) = xx(2 * k) +!if(exp_mass(kk) > imax) imax = exp_mass(kk) +!if(exp_int (kk) < imin) imin = exp_int (kk) + enddo + enddo + endif + enddo rd + + imax = maxval(exp_mass) +!imin = minval(exp_mass) + + mmax = maxval(exp_int) + +! bring exp to the right order of magnitude + do i = 1, exp_entries + exp_int(i) = exp_int(i) / norm + enddo + endif + +!!!!!!!!!!!!!!!!!!!!! +!> if it is a csv file + + + + if (expdat) then + write(*,*) 'Reading ', exp_dat,' ...' + open( file = exp_dat, newunit = io_exp, status = 'old') + + exp_entries = 0 + +!>> count the entries + do + read(io_exp,'(a)',iostat=iocheck)!line + if (iocheck < 0) exit + + exp_entries = exp_entries + 1 + enddo + + rewind(io_exp) + +!>> allocate memory + allocate (exp_mass(exp_entries), & + & exp_int(exp_entries)) + + kk = 0 + iocheck=0 + +!>> count the intensities and masses + do kk = 1, exp_entries + read(io_exp,'(a)',iostat=iocheck)line + if (iocheck < 0) exit + + do k=1,80 + if(line(k:k) == ',') then + line(k:k)=' ' + endif + enddo + + call readl(line,xx,nn) + + exp_mass(kk) = xx(1) + exp_int (kk) = xx(2) + enddo + close(io_exp) + + imax = maxval(exp_mass) +!imin = minval(exp_mass) + + mmax = maxval(exp_int) + + do kk = 1, exp_entries + exp_int (kk) = (exp_int(kk) / mmax) *100.0_wp + enddo + + + endif + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!> get the minimum m/z value that is to be plotted + imin = mzmin + +!> get the maximum m/z value that is to be plotted + if (imax < nint(maxval(sorted_masses))) imax = nint(maxval(sorted_masses)) + +! counts of structures in the 100% peak + write(*,*) 'Theoretical counts in 100 % signal:', idint(tmax) + +!----------------------------------------------------------------------------- +!> verbose information -> not yet available anymore with acc masses +!if(verbose)then +! write(*,*)'mass % intensity counts Int. exptl' +! do i = 1, list_length !count_mass +! if (sorted_masses(i) > mzmin) then +! dum = 0 +! do j = 1, exp_entries +! !if(int(iexp(1,j)) == i) dum = iexp(2,j) +! if(int(iexp(1,j)) == i) dum = iexp(2,j) +! enddo + +! write(*,'(i8,F8.2,i8,F8.2)') & +! !i, 100.0_wp *sorted_masses(i)/tmax, idint(tmass(i)), dum/norm +! !i, 100.0_wp *sorted_intensities(i)/tmax, sorted_masses(i), dum/norm +! i, sorted_intensities(i), sorted_masses(i), dum/norm +! endif +! enddo +!endif +!----------------------------------------------------------------------------- + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! P R I N T +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + write(*,*) + write(*,'('' Writing ...'')') +!> print out comma-separated-value (csv) file + + write(*,'('' - result.csv file'')') + open( file = 'result.csv', newunit = io_csv, status='replace') + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!> print out JCAMP file + + write(*,'('' - result.jdx file'')') + open( file = 'result.jdx', newunit = io_jcamp, status='replace') + +! TK minimum JCAMP-DX header + write(io_jcamp,"(A)")'##TITLE=Theoretical in-silico spectrum (QCxMS)' + write(io_jcamp,"(A)")'##JCAMP-DX=4.24' + write(io_jcamp,"(A)")'##DATA TYPE=MASS SPECTRUM' + + write(io_jcamp,"(A)")'##XUNITS=M/Z' + write(io_jcamp,"(A)")'##YUNITS=RELATIVE INTENSITY' + +!> number of in-silico spectra + write(io_jcamp,'(A, I0)')'##NPOINTS=' , list_length + write(io_jcamp,"(A)")'##PEAK TABLE=(XY..XY) 1' + +! >>>> Write out the results into the files <<<< ! + do i = 1, list_length + write(io_csv, '(f12.6, 1x, a1, 3x, f26.18)') & + sorted_masses(i),',', sorted_intensities(i) + write(io_jcamp,'(f12.6,4x, f26.18)') & + sorted_masses(i), sorted_intensities(i) + enddo + + write(io_jcamp,"(A)")'##END=' + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!> print out XMGRACE file + +!> when the template exists + xname = trim(fullpath(xname)) + inquire(file=xname, exist=ex) + if(ex)then + open(file=xname, newunit=io_raw) + +! Original file called .mass_raw.agr + if( small ) open( file='small_mass.agr', newunit=io_mass) + if( .not. small ) open( file='mass.agr', newunit=io_mass) + + write(*,'('' - mass.agr file'')') + do i = 1, nagrfile + read (io_raw,'(a)') line + if( .not. small )then + if(index(line,'world 10') /= 0)then + write(line,'(''@ world '',i3,'', -105,'',i4,'', 100'')') imin, imax + 5 + endif + if(index(line,'xaxis tick major') /= 0)then + line='@ xaxis tick major 20' + if(imax > 200) line='@ xaxis tick major 50' + endif + + if(index(line,'@ s1 symbol size') /= 0)then + line='@ s1 symbol size 0.160000' + if(imax > 200) line='@ s1 symbol size 0.100000' + endif + + if(index(line,'@ s2 symbol size') /= 0)then + line='@ s2 symbol size 0.160000' + if(imax > 200) line='@ s2 symbol size 0.100000' + endif +!kleines Format + else + if(index(line,'world 10') /= 0)then + write(line,'(''@ world '',i3,'', 0,'',i3,'', 100'')') imin, imax + 5 + endif +!if(index(line,'xaxis tick major') /= 0)then +!line='@ xaxis tick major 10' +!if(imax > 200) line='@ xaxis tick major 50' +!endif + + if(index(line,'@ s1 symbol size') /= 0)then + line='@ s1 symbol size 0.320000' +!if(imax > 200) line='@ s1 symbol size 0.100000' + endif + endif + write(io_mass,'(a)')line + enddo + +!> Write out the results into the mass.agr file ! + do i = 1, list_length ! count_mass + write(io_mass,*) sorted_masses(i), sorted_intensities(i) + enddo + + write(io_mass,*)'&' +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TK establish again if expdat (experimental JCAMPDX) exists +! TK exp_mass(i) - contains the exp masses and exp_int(i) contains the abundance +! TK exp_entries contains number of experimental spectra + + if(nistdat .or. expdat)then + + write(io_mass,*)'@target G0.S2' + write(io_mass,*)'@type bar' + do i=1,exp_entries + write(io_mass,*) exp_mass(i),-exp_int(i) + enddo + write(io_mass,*)'&' + + endif + + endif !ex + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +! compute deviation exp-theor. + if(expdat .or. nistdat)then + + rms = 0 + nn = 0 + kkk = 0 + kkkk = 0 + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!> this is for integer masses +!> calculate integer list of theor. masses +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + store_thr_int = 0 + + allocate(int_thr(list_length)) + +!> transfer to integer + do j = 1, list_length + int_thr(j) = nint(sorted_masses(j)) + enddo + + + new_length_thr = 0 + sm = 0 + + do i = 1, list_length + store_intensities = 0.0_wp + + store_check_list = int_thr(i) + + if(sm /= store_check_list) new_length_thr = new_length_thr + 1 + + do j = 1, list_length + if (store_check_list == int_thr(j))then + store_intensities = store_intensities + sorted_intensities(j) + endif + + enddo + + sm = store_check_list + + store_thr_int (new_length_thr) = store_intensities + store_thr_mass(new_length_thr) = store_check_list + +!write(*,*) store_check_list, new_length_thr, added_intensities(new_length_thr) + + enddo + + allocate(added_thr_intensities(new_length_thr), & + added_thr_masses(new_length_thr)) + + do i = 1, new_length_thr + added_thr_masses(i) = store_thr_mass(i) + added_thr_intensities(i) = (store_thr_int(i)/maxval(store_thr_int)*100.0_wp) +! write(*,*) added_masses(i), added_intensities(i) + enddo + + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!> for experimental masses in integer +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + store_exp_int = 0 + allocate(int_exp(exp_entries)) + + do j = 1, exp_entries + int_exp(j) = nint(exp_mass(j)) + enddo + + new_length_exp = 0 + sm = 0 + + do i = 1, exp_entries + store_intensities = 0.0_wp + + store_check_list = int_exp(i) + + if(sm /= store_check_list) new_length_exp = new_length_exp + 1 + + do j = 1, exp_entries + if (store_check_list == int_exp(j))then + store_intensities = store_intensities + exp_int(j) + endif + + enddo + + sm = store_check_list + + store_exp_int(new_length_exp) = store_intensities + store_exp_mass(new_length_exp) = store_check_list + +!write(*,*) store_check_list, new_length_exp, added_intensities(new_length_exp) + + enddo + + allocate(added_exp_intensities(new_length_exp), & + added_exp_masses(new_length_exp)) + + do i = 1, new_length_exp + added_exp_masses(i) = store_exp_mass(i) + added_exp_intensities(i) = (store_exp_int(i)/maxval(store_exp_int)*100.0_wp) +! write(*,*) added_masses(i), added_intensities(i) + enddo + + +ol: do i = 1, new_length_exp +il: do j = 1, new_length_thr + if ( added_exp_intensities(i) > 5.0_wp )then +!if ( nint(exp_mass(i)) == added_masses(j) ) then + if ( added_exp_masses(i) == added_thr_masses(j) ) then + r = added_thr_intensities(j) - added_exp_intensities(i) + kkk = kkk + 1 + if ( added_thr_intensities(j) > 2.5_wp) kkkk = kkkk + 1 + rms = rms + abs(r) + endif + endif + enddo il + enddo ol + +! new_length = 0 +! sm = 0 +! do i = 1, list_length +! +! store_check_list = sorted_masses(i) +! +! store_intensities = 0.0_wp +! +! do i = 1, list_length +! if (sorted_masses(i) - store_check_list <= 1.0d-2)then +! store_intensities = store_intensities + sorted_intensities(j) +! store_mass = store_intensities + sorted_intensities(j) +!! store_intensities = sorted_masses(i) - sorted_masses(i-1) +!! write(*,*) i, sorted_masses(i-1), sorted_masses(i), store_intensities +! endif +! enddo +! +! if(sm /= store_check_list) new_length = new_length + 1 +! sm = store_check_list +! +! stop +! +! +! +!ol: do i = 1, exp_entries +!il: do j = 1, list_length +! if ( exp_int(i) > 5.0_wp )then +! !if ( nint(exp_mass(i)) == added_masses(j) ) then +! if ( exp_mass(i) == sorted_masses(j) ) then +! r = added_intensities(j) - exp_int(i) +! kkk = kkk + 1 +! if ( added_intensities(j) > 2.5_wp) kkkk = kkkk + 1 +! rms = rms + abs(r) +! endif +! endif +! enddo il +! enddo ol + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + write(*,*)'MAD(exptl./theor.) = ', rms / kkk + write(*,*)'# exptl. > 5 % = ', kkk + write(*,*)'% correctly found = ', 100 * float(kkkk)/float(kkk) + + + close(io_spec) + close(io_raw) + close(io_mass) + close(io_csv) + close(io_jcamp) + +! compute mass spectral match score + score = 0.0_wp + +!write(*,*) +!write(*,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' +!write(*,*) ' - NIST mass comparison not supported in v6.0 - ' +!write(*,*) ' - no automatic score comparison -' +!write(*,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' +!stop ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' + + +!call match_accurate(sorted_masses, sorted_intensities, list_length, & +! exp_entries, exp_mass, exp_int,score,tmax) +!call match(added_masses, added_intensities, new_length, & +! exp_entries, exp_mass, exp_int,score,tmax) + call match(added_thr_masses, added_thr_intensities, new_length_thr, & + new_length_exp, added_exp_masses, added_exp_intensities,score,tmax) + write(*,*) + write(*,*)"!!!!!!!!!!!!!!!!!!!!!!! " + write(*,*)" Matching score: " + write(*,'(6F10.3)') score + write(*,*)"!!!!!!!!!!!!!!!!!!!!!!! " + write(*,*) + write(*,*)"Composite match score, see " + write(*,*)"Stein, S. E.; Scott, D. R. J. Am. Soc. Mass Spectrom. 1994, 5, 859-866" +!write(*,*)"For our implementation, see " +!write(*,*)"Bauer, C. A.; Grimme,S. J. Phys. Chem. A 2014, 118, 11479-11484" + + + endif + +end program plotms + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +subroutine match(added_masses, added_intensities, new_length, & + exp_entries, exp_mass, exp_int,score,tmax) + use xtb_mctc_accuracy, only: wp + implicit none + + integer :: i,j + integer :: pp !peak pair number pp + integer :: numcalc ! number of calculated peaks + integer :: new_length, exp_entries +! integer :: entry, entry2 + integer :: cnt + +!real(wp) :: tmass(10000) +!real(wp) :: iexp(2,10000) + real(wp) :: exp_mass(exp_entries) + real(wp) :: exp_int(exp_entries) + real(wp) :: score + real(wp) :: tmax + real(wp),allocatable :: w_exp(:) !weighted vectors + real(wp),allocatable :: w_thr(:) !weighted vectors +!real(wp) :: w(2,10000) !weighted vectors +!real(wp), allocatable :: w_exp(:) !weighted vectors +!real(wp), allocatable :: w_thr(:) !weighted vectors + real(wp) :: sum1,sum2,sum3,sum4,dot,fr +!real(wp) :: norm(2) + real(wp) :: norm +!real(wp) :: p(2,10000) ! peak pair matrix + real(wp), allocatable :: p(:,:) ! peak pair matrix + real(wp) :: added_masses(new_length) + real(wp) :: added_intensities(new_length) +!!!!!!!!!!!!!!! +!integer :: exp_mass(exp_entries) +!integer :: added_masses(new_length) +!!!!!!!!!!!!!!! + + +! numcalc = new_length +! entry = 0 + +!allocate( w_exp(exp_entries), & +! w_thr(new_length)) + +!> weighted spectral vectors, m**3, int**0.6 scaling +!w = 0.0_wp + + pp = 0 +!> get weighting of experimental spectrum + allocate(w_exp(exp_entries)) + norm = 0.0_wp + do i = 1, exp_entries +! w_exp(i) = exp_mass(i)**2 * exp_int(i)**0.6_wp + w_exp(i) = exp_mass(i)**3 * exp_int(i)**0.6_wp + norm = norm + w_exp(i)**2 + enddo + + norm = sqrt(norm) + + do i = 1, exp_entries + w_exp(i) = w_exp(i) / norm + enddo + +!> get weighting of calculated spectrum + allocate(w_thr(new_length)) + norm = 0.0_wp + do j = 1, new_length +!w_thr(j) = (added_masses(j))**2 * ( added_intensities(j) )**0.6_wp + w_thr(j) = (added_masses(j))**3 * ( added_intensities(j) )**0.6_wp + norm = norm + w_thr(j)**2 + enddo + + norm = sqrt(norm) + + do j = 1, new_length + w_thr(j) = w_thr(j) / norm + enddo + + + dot = 0.0_wp + sum1 = 0.0_wp + sum2 = 0.0_wp + sum3 = 0.0_wp + + + do i = 1, exp_entries + do j = 1, new_length + if ( exp_mass(i) == added_masses(j)) then + pp = pp + 1 + sum1 = sum1 + w_exp(i) * w_thr(j) + sum2 = sum2 + (w_exp(i))**2 + sum3 = sum3 + (w_thr(j))**2 + endif + enddo + enddo + + + dot = sum1**2 / (sum2 * sum3) +!write(*,*) 'DOT', dot + + deallocate(w_exp, w_thr) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! masses and intensity scalement for the second term +! m**0 +! int**1 + allocate(w_exp(exp_entries), w_thr(new_length)) + + norm=0.0_wp + + do i=1, exp_entries +!norm = norm + w_exp(i)**2 ! + exp_int(i) !w(1,i)**2 + norm = norm + exp_int(i)**2 !w(1,i)**2 + enddo + + norm = sqrt(norm) + + do i=1, exp_entries + w_exp(i) = exp_int(i) / norm + enddo + +!!!!!!!!!!!!!!!!!!!! + norm=0.0_wp + do j=1, new_length + norm = norm + added_intensities(j)**2 !w(2,i)**2 + enddo + + norm = sqrt(norm) + + do j=1, new_length + w_thr(j) = added_intensities(j) / norm + enddo + +!pp = 0 +!p = 0.0_wp + sum4 = 0.0_wp + fr = 0.0_wp + cnt = 0 + + allocate (p(2,pp)) + + do i = 1, exp_entries + do j = 1, new_length + if ( exp_mass(i) == added_masses(j)) then + cnt = cnt + 1 + p(1,cnt) = w_exp(i) + p(2,cnt) = w_thr(j) + endif + enddo + enddo + +! ardous loop + call calcfr(pp,p,sum4) + if (pp == new_length) sum4 = sum4 + 1.0_wp +!write(*,*) 'SUM4', sum4 +!write(*,*) 'numcalc', numcalc +!write(*,*) 'numcalc', new_length + + fr = sum4 / float(pp) + score = (new_length * dot + pp * fr) / (new_length + pp) + write(*,*) "Fr is:", fr + write(*,*) "Dot is:", dot + write(*,*) "new_length", new_length + write(*,*) "PP is:", pp + return +end subroutine match + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! computes the second expression in the mass spec match score +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +subroutine calcfr(pp,pair,sum4) + use xtb_mctc_accuracy, only: wp + implicit none + + integer :: i + integer :: pp + + real(wp) :: pair(2,pp) + real(wp) :: sum4 + + sum4 = 0.0_wp + + do i = 2, pp + if (abs((pair(1,i) * pair(2,i-1)) / (pair(1,i-1) * pair(2,i))) < 1.0_wp)then + + sum4 = sum4 + (pair(1,i) * pair(2,i-1)) / (pair(1,i-1) * pair(2,i)) + + elseif (abs((pair(1,i) * pair(2,i-1)) / (pair(1,i-1) * pair(2,i))) > 1.0d0)then + + sum4 = sum4 + 1 / ((pair(1,i) * pair(2,i-1)) / (pair(1,i-1) * pair(2,i))) + + elseif (abs((pair(1,i) * pair(2,i-1)) / (pair(1,i-1) * pair(2,i))) == 1.0d0)then + sum4 = sum4 + 1.0_wp + endif + enddo + +end subroutine calcfr + + +subroutine print_help() +implicit none + + write(*,*) + write(*,'(1x, ''usage :'')') + write(*,'(1x, ''plotms [options]'')') + write(*,*) + write(*,'(1x, ''plotms needs a qcxms.res qcxms_cid.res tmpqcxms.res or tmpqcxms_cid.res file from a qcxms calculation in the current directory.'')') + write(*,*) + write(*,'(/,1x,''command line arguments'')') + write(*,'(5x,''-a: count charges from -1000 (cthr) to cthr2'')') + write(*,'(5x,''-v: print spectra "mass % intensity counts Int. exptl" to stdout;'')') + write(*,'(5x,'' with "Int. exptl" (experimental) taken from exp.dat but not all'')') + write(*,'(5x,'' exp peaks are exported if no theoretical counterpart exists'')') + write(*,'(5x,''-f : give name of res file'')') + write(*,'(5x,''-t: couting ions with charge x to y (give the value, e.g. "-t 1 2" for charge 1 to 2)'')') + write(*,'(5x,''-w [real]: broadening the charges by an SD '')') + write(*,'(5x,''-s: Take only secondary and tertiary fragmentations (give the value, e.g. "-s 2" for secondary)'')') + write(*,'(5x,''-m [int] set minimum value for m/z, so 100% value will be calc. for higher values (CID)'')') + write(*,'(5x,''-p calculate NO isotope pattern, highest peak of isotope pattern is used'')') + write(*,'(5x,''-i [real]: set minimum intensity that is considered in rel. intensity percent'')') + write(*,'(5x,''--int_masses : round m/z values to integers'')') + write(*,'(5x,''-e: provide exp. input file (in .csv format)'')') + write(*,'(5x,''--version: print version of this program'')') + stop ' [-h] displayed. exit.' +end subroutine print_help + +! Get the absolute path to file +! Note, this function replaces shell tilde ~/ with the explicit home dir string. +function fullpath(fname) + implicit none + character(len=*), intent(in) :: fname + character(len=4096) :: absdir, fullpath + if(fname(1:1) .eq. '/') then + fullpath = fname + else if(fname(1:2) .eq. '~/') then + call getenv('HOME', absdir) + fullpath = absdir(:len_trim(absdir)) // fname(2:len_trim(fname)) + else + call getcwd(absdir) + fullpath = absdir(:len_trim(absdir)) // '/' // fname(:len_trim(fname)) + endif +end function fullpath + diff --git a/subprojects/json-fortran-8.2.5.wrap b/subprojects/json-fortran-8.2.5.wrap new file mode 100644 index 0000000..f3675a6 --- /dev/null +++ b/subprojects/json-fortran-8.2.5.wrap @@ -0,0 +1,2 @@ +[wrap-redirect] +filename = mctc-lib/subprojects/json-fortran-8.2.5.wrap From 77f762a5d75fc0afc8dc9ca15e48da4bec323344 Mon Sep 17 00:00:00 2001 From: Johannes Gorges <58849467+gorges97@users.noreply.github.com> Date: Fri, 20 Jun 2025 11:20:22 +0200 Subject: [PATCH 08/11] removed uneccesarry output --- src/isotopes.f90 | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/isotopes.f90 b/src/isotopes.f90 index 877d5ee..fccc15f 100644 --- a/src/isotopes.f90 +++ b/src/isotopes.f90 @@ -986,12 +986,7 @@ subroutine isotope(counter, mzmin, ntot, iat_save, maxatm, rnd, & enddo - write(*,*) 'number of peaks', index_mass - do i = 1, index_mass - write(*,*) list_masses(i), store_int(i) - enddo - - + allocate(exact_intensity(index_mass)) allocate(isotope_masses(index_mass)) From b1b6f1cd5b14e950123d0cc95857dad5757e916e Mon Sep 17 00:00:00 2001 From: Johannes Gorges <58849467+gorges97@users.noreply.github.com> Date: Thu, 10 Jul 2025 11:07:26 +0200 Subject: [PATCH 09/11] Update build.yml --- .github/workflows/build.yml | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 4c19df2..05431f0 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -51,8 +51,14 @@ jobs: PYTHON_V: 3.8 steps: - - name: Checkout code + - name: Checkout code (with submodules) uses: actions/checkout@v4 + with: + submodules: recursive + - name: Prepare Meson subprojects + run: | + meson subprojects download + - uses: actions/setup-python@v4 with: From 3d54d70fa6935fdb83f3eed4c87e402bf5aee107 Mon Sep 17 00:00:00 2001 From: Johannes Gorges <58849467+gorges97@users.noreply.github.com> Date: Thu, 10 Jul 2025 11:11:11 +0200 Subject: [PATCH 10/11] removed json.wrap for subproject --- subprojects/json-fortran-8.2.5.wrap | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 subprojects/json-fortran-8.2.5.wrap diff --git a/subprojects/json-fortran-8.2.5.wrap b/subprojects/json-fortran-8.2.5.wrap deleted file mode 100644 index f3675a6..0000000 --- a/subprojects/json-fortran-8.2.5.wrap +++ /dev/null @@ -1,2 +0,0 @@ -[wrap-redirect] -filename = mctc-lib/subprojects/json-fortran-8.2.5.wrap From 3fa7fc178a50e611d63d41c6f218678fde25ee76 Mon Sep 17 00:00:00 2001 From: Johannes Gorges <58849467+gorges97@users.noreply.github.com> Date: Thu, 10 Jul 2025 11:12:47 +0200 Subject: [PATCH 11/11] reworked build file --- .github/workflows/build.yml | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 05431f0..74463f6 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -51,13 +51,8 @@ jobs: PYTHON_V: 3.8 steps: - - name: Checkout code (with submodules) + - name: Checkout code uses: actions/checkout@v4 - with: - submodules: recursive - - name: Prepare Meson subprojects - run: | - meson subprojects download - uses: actions/setup-python@v4