Small cosmetic changes

implemented in ARPACK for calculating a subset of eigenvalues rather than all of
them or just the largest one.

This is invoked by setting algorithm='arnoldi' in the constructor MSSPEC.

    calc = MSSPEC(spectroscopy='EIG', algorithm='arnoldi')

This version finds NEV eigenvalues with the largest magnitudes. The value of
NEV is set in the code to be 2% of the total number. The algorithm also sets
the number of basis vectors, NCV, to be 2*NEV. Options for allowing the user
to set these values as well as other convergence criteria should be added in
later versions.

Files src/msspec/calculator.py and src/msspec/parameters.py have been modified
to accept algorithm='arnoldi'.

A new directory, src/msspec/spec/fortran/eig/ar, has been created and contains
various Fortran files implementing the Arnoldi iteration. Some files are written
in Fortran90 free-format and use some features of the Fortran 2003 standard.

Compilation requires the ARPACK library to be installed and accessible using -larpack

Compilation rules and options have been introduced/modified where necessary in the
files src/options.mk, src/msspec/spec/fortran/Makefile and
src/msspec/spec/fortran/eig_ar.mk

This version has been tested for the CaO example.

Dockerfile

@@ -1,24 +1,102 @@
 # Get the base Python image
-FROM python:latest
+FROM alpine:edge AS builder
+
+# Variables
+ARG branch="devel"
+ARG login="" password=""
+ARG folder=/opt/msspec user=msspec
 
 # Install system dependencies
-RUN apt-get update && apt-get install -y virtualenv gfortran libgtk-3-dev nano
+# tools
-
+RUN apk add bash git make gfortran python3 
-# Add a non-privileged user
+# headers
-RUN useradd -ms /bin/bash -d /opt/msspec msspec
+RUN apk add python3-dev lapack-dev musl-dev
-
+# python packages
-# Set the working directory in the container
+RUN apk add py3-virtualenv py3-pip py3-numpy-dev py3-h5py py3-lxml py3-matplotlib py3-numpy py3-pandas py3-cairo py3-scipy py3-setuptools_scm
-USER msspec
+RUN apk add --no-cache -X http://dl-cdn.alpinelinux.org/alpine/edge/testing py3-wxpython
-RUN mkdir -p /opt/msspec/code
+RUN pip install ase pint terminaltables ipython
-WORKDIR /opt/msspec/code
+# for GUI
+RUN apk add ttf-droid adwaita-icon-theme
+RUN apk add build-base
 
 # Fetch the code
-RUN git clone https://git.ipr.univ-rennes1.fr/epsi/msspec_python3.git .
+RUN mkdir -p ${folder}/code
-#COPY --chown=msspec:msspec . .
+WORKDIR ${folder}/code
+RUN git clone --branch ${branch} https://${login}:${password}@git.ipr.univ-rennes1.fr/epsi/msspec_python3.git .
 
+# Build
+RUN make pybinding NO_VENV=1 PYTHON=python3 VERBOSE=1
+RUN make -C src sdist PYTHON=python3 NO_VENV=1 VENV_PATH=${folder}/.local/src/msspec_venv && \
+    pip install src/dist/msspec*tar.gz
+
+# Add a non-privileged user
+RUN adduser -D -s /bin/bash -h ${folder} ${user}
+
+# Set the working directory in the container
+USER ${user}
+
+RUN virtualenv --system-site-packages ${folder}/.local/src/msspec_venv
+RUN make -C src frontend PYTHON=python3 NO_VENV=1 VENV_PATH=${folder}/.local/src/msspec_venv
+
+
+
+FROM alpine:edge
+# Variables
+ARG folder=/opt/msspec user=msspec
+# Install system dependencies
+RUN apk add --no-cache -X http://dl-cdn.alpinelinux.org/alpine/edge/testing \
+#  hdf5-hl cairo openblas lapack libxml2 libxslt libzlf wxwidgets-gtk3 openjpeg libimagequant \
+   nano \
+   py3-virtualenv \
+   lapack \
+   bash \
+#  git \
+#  make \
+#  gfortran \
+   python3 \
+#   ttf-droid \
+   ttf-liberation \
+   adwaita-xfce-icon-theme \
+#  python3-dev \
+#  lapack-dev \
+#  musl-dev \
+#  py3-virtualenv \
+   py3-pip \
+#  py3-numpy-dev \
+   py3-h5py \
+   py3-lxml \
+   py3-matplotlib \
+   py3-numpy \
+   py3-pandas \
+   py3-cairo \
+   py3-scipy \
+   py3-setuptools_scm \
+   py3-wxpython \
+   && pip install \
+   ase \
+   pint \
+   terminaltables \
+   ipython \
+   && pip cache purge \
+   # Add a non-privileged user
+   && adduser -D -s /bin/bash -h ${folder} ${user}
+
+# Set the working directory in the container
+USER ${user}
+WORKDIR ${folder}
 # Install msspec
-ENV PATH=/opt/msspec/.local/bin:$PATH
+#COPY --from=builder ${folder}/.local ${folder}/.local
-RUN make install VERBOSE=1
+#COPY --from=builder /usr/lib/python3.10/site-packages /usr/lib/python3.10/site-packages
+
+COPY --from=builder ${folder}/code/src/dist/msspec*tar.gz msspec.tar.gz
+RUN virtualenv --system-site-packages .local/src/msspec_venv && \
+    . .local/src/msspec_venv/bin/activate && \
+    pip install msspec.tar.gz && \
+    rm -f msspec.tar.gz && \
+    mkdir -p .local/bin
+
+COPY --from=builder ${folder}/.local/bin/msspec .local/bin/msspec
+ENV PATH=${folder}/.local/bin:$PATH
 
 # Run the msspec frontend command on startup
 ENTRYPOINT ["msspec"]

Jenkinsfile

@@ -13,9 +13,9 @@ pipeline {
         }
         stage('Syncing website...') {
             steps {
-		// echo 'Syncing website...'
+		echo 'Syncing website only in master branch, not here in devel branch...'
 		// sh 'rm -rf $HOME/www/*'
-                sh 'cp -a ./doc/build/html/* $HOME/www/'
+                // sh 'cp -a ./doc/build/html/* $HOME/www/'
             }
         }
 

Makefile

@@ -1,7 +1,7 @@
 include src/options.mk
 
 
-.PHONY: pybinding install devel venv doc clean
+.PHONY: pybinding install devel venv doc clean _attrdict
 
 
 pybinding:
@@ -10,20 +10,23 @@ pybinding:
 
 venv:
 ifeq ($(NO_VENV),0)
-	@virtualenv --python=$(PYTHON_EXE) --prompt="(msspec-$(VERSION)) " $(VENV_PATH)
+#	@virtualenv --python=$(PYTHON_EXE) --prompt="(msspec-$(VERSION)) " $(VENV_PATH)
-	$(INSIDE_VENV) \
+	$(PYTHON) -m venv $(VENV_PATH) 
-            wget https://bootstrap.pypa.io/get-pip.py && \
+	$(INSIDE_VENV) python -m ensurepip --upgrade
-            python get-pip.py                         && \
+	$(INSIDE_VENV) pip install -r src/pip.freeze
-            pip install --upgrade setuptools          && \
-            pip install -r src/pip.freeze             && \
-            rm -f get-pip.py
 endif
 
+# wget https://bootstrap.pypa.io/get-pip.py && \
+# python get-pip.py                         && \
+# rm -f get-pip.py
+# pip install --upgrade setuptools          && \
+# pip install -r src/pip.freeze             && \
+
 
 install: venv pybinding wx
 	@+$(INSIDE_VENV) $(MAKE) -C src sdist
 	@+$(INSIDE_VENV) $(MAKE) -C src frontend
-	@+$(INSIDE_VENV) pip install src/dist/msspec-$(VERSION).tar.gz
+	@+$(INSIDE_VENV) pip install --force-reinstall src/dist/msspec-$(VERSION)*.whl
 	@echo "Do not forget to check that $(INSTALL_PREFIX)/bin is set in your \$$PATH"
 
 
@@ -37,20 +40,28 @@ light: venv
 	@$(INSIDE_VENV) pip install src/
 
 
-_build_wx/wxPython.target:
+_attrdict:
+	# Check if virtualenv python version > 3.3.0
+	# If so, install the patched version of attrdict used to build the version 4.2.0 of wxPython
+	@$(INSIDE_VENV) if `python -c "import sys; exit(sys.version_info > (3,3))"`; then \
+		pip install --no-cache attrdict; \
+	else \
+		pip install thirdparty/attrdict-2.0.1.tar.gz; \
+	fi
+
+
+_build_wx/wxPython.target: _attrdict
 	@$(INSIDE_VENV) echo "Building wxPython for your `python --version 2>&1` under Linux $(DISTRO_RELEASE)..."
 	# Create a folder to build wx into
 	@mkdir -p _build_wx
-	@$(INSIDE_VENV) pip install attrdict sip
-	# TODO: attrdict is no longer compatible with collections package. The build will fail
 	# download the wheel or the source if it cannot find a wheel
-	@$(INSIDE_VENV) cd _build_wx && pip download -f https://extras.wxpython.org/wxPython4/extras/linux/gtk3/$(DISTRO_RELEASE) wxPython
+	$(INSIDE_VENV) cd _build_wx && pip download -f https://extras.wxpython.org/wxPython4/extras/linux/gtk3/$(DISTRO_RELEASE) wxPython
 	# Build the source if a tar.gz was downloaded
 	@$(INSIDE_VENV) cd _build_wx && \
           if [ -e wxPython*.tar.gz ]; then \
             tar -x --skip-old-files -vzf wxPython*.tar.gz; \
             cd `ls -d wxPython*/`; \
-            pip install requests; \
+            pip install requests sip; \
             python build.py dox etg --nodoc sip build bdist_wheel; \
             ln -sf `readlink -f dist/wxPython*.whl` ../; \
           fi;

doc/source/tutorials/RhO/RhO.py

@@ -18,7 +18,8 @@ for zi, z0 in enumerate(all_z):
     calc.set_atoms(cluster)
 
     # Compute
-    data = calc.get_theta_phi_scan(level='1s', kinetic_energy=723, data=data)
+    data = calc.get_theta_phi_scan(level='1s', kinetic_energy=723, data=data,
+                                   malloc={'NPH_M': 8000})
     dset = data[-1]
     dset.title = "{:d}) z = {:.2f} angstroms".format(zi, z0)
 

src/MANIFEST.in

@@ -1,6 +0,0 @@
-recursive-include msspec *.so
-recursive-include . SConstruct
-include setup_requirements.txt
-include requirements.txt
-include pip.freeze
-include VERSION

src/Makefile

@@ -9,16 +9,15 @@ sdist: dist/msspec-$(VERSION).tar.gz
 frontend: $(INSTALL_PREFIX)/bin/msspec
 
 
-dist/msspec-$(VERSION).tar.gz: VERSION
+dist/msspec-$(VERSION).tar.gz: msspec/VERSION
 	@echo "Creating Python source distribution..."
-	@python setup.py sdist
+	@+$(INSIDE_VENV) pip install build && python -m build
 
 
-$(INSTALL_PREFIX)/bin/msspec: msspec.sh.template VERSION
+$(INSTALL_PREFIX)/bin/msspec: msspec.sh.template msspec/VERSION
 	@echo "Installing frontend command..."
 	@mkdir -p $(dir $@)
-	@cat $< | sed -e 's/__VERSION__/$(VERSION)/' -e 's#__VENV_PATH__#$(VENV_PATH)#' > $@
+	@cat $< | sed -e 's#__VENV_PATH__#$(VENV_PATH)#' > $@
-	#@cat $< | sed 's/__VERSION__/$(VERSION)/' > $@
 	@chmod 755 $@
 
 
@@ -26,7 +25,7 @@ pybinding:
 	@echo "Building Python binding for phagen and spec..."
 	@+$(MAKE) -C msspec/phagen/fortran all
 	@+$(MAKE) -C msspec/spec/fortran all
-	@echo "$(VERSION)" > VERSION
+	@echo "$(VERSION)" > msspec/VERSION
 
 
 results: msspec/results.txt
@@ -54,7 +53,7 @@ clean::
 #	remove previous sdist
 	@rm -rf dist
 	@rm -rf *.egg*
-	@rm -f VERSION
+	@rm -f msspec/VERSION
 
 
 help:

src/msspec.sh.template

@@ -2,12 +2,11 @@
 
 SCRIPT_PATH="$0"
 SCRIPT_NAME=$(basename "$SCRIPT_PATH")
-VERSION="__VERSION__"
 VENV_PATH="__VENV_PATH__"
 
 # Check venv path
 if ! [ -d "$VENV_PATH" ]; then
-   echo "ERROR: Unable to find version $VERSION of msspec!!"
+   echo "ERROR: Unable to find msspec!!"
    exit 1
 fi
 
@@ -15,6 +14,10 @@ launch_script() {
    . "$VENV_PATH/bin/activate" && python "$@"
 }
 
+show_version () {
+    . "$VENV_PATH/bin/activate" && python -c "import msspec; print(msspec.__version__)"
+}
+
 show_help () {
    echo "Usage: 1) $SCRIPT_NAME -p [PYTHON OPTIONS] SCRIPT [ARGUMENTS...]"
    echo "       2) $SCRIPT_NAME [-l FILE | -i | -h]"
@@ -92,7 +95,7 @@ while getopts "hvil:p:eu" option; do
          ;;
       u) uninstall
          ;;
-      v) echo $VERSION
+      v) show_version
          ;;
       *|h) show_help
          ;;

src/msspec/__init__.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -16,8 +17,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/__init__.py
-# Last modified: ven. 10 avril 2020 17:22:12
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
-# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 import ase

src/msspec/calcio.py

@@ -1,5 +1,25 @@
+#!/usr/bin/env python
 # coding: utf-8
-# vim: set et ts=4 sw=4 fdm=indent mouse=a cc=+1 tw=80:
+#
+# Copyright © 2016-2020 - Rennes Physics Institute
+#
+# This file is part of msspec.
+#
+# msspec is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+# msspec is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# You should have received a copy of the GNU General Public License
+# along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Source file  : src/msspec/calcio.py
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
+
 
 """
 Module calcio
@@ -910,7 +930,7 @@ class SpecIO(object):
         if content != old_content:
             with open(filename, 'w') as fd:
                 fd.write(content)
-                LOGGER.debug(f"Writing Spec input file written in {filename}")
+                LOGGER.debug("Writing Spec input file written in {}".format(filename))
             modified = True
 
         return modified
@@ -1255,13 +1275,13 @@ class CompCurveIO(object):
         data = []
         for i in range(1, 13):
             #data.append(np.loadtxt(prefix + f'{i:02d}' + '.txt')[-1])
-            results = np.loadtxt(prefix + f'{i:02d}' + '.txt')
+            results = np.loadtxt(prefix + '{:02d}'.format(i) + '.txt')
             results = results.reshape((-1, 2))
             data.append(results[index,1])
         suffix = 'ren'
         exp = {'int': None, 'ren': None, 'chi': None, 'cdf': None}
         exp_ren = np.loadtxt(os.path.join('exp', 'div',
-                                          f'experiment_{suffix}.txt'))
+                                          'experiment_{}.txt'.format(suffix)))
         calc_ren = np.loadtxt(os.path.join('calc', 'div',
-                                           f'calculation{index:d}_{suffix}.txt'))
+                                           'calculation{:d}_{}.txt'.format(index,suffix)))
         return data, exp_ren, calc_ren

src/msspec/calculator.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -16,8 +17,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/calculator.py
-# Last modified: ven. 10 avril 2020 17:19:24
+# Last modified: Tue, 25 Oct 2022 16:21:38 +0200
-# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
+# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1666707698 +0200
 
 
 """
@@ -94,8 +95,10 @@ from msspec.parameters import TMatrixParameters
 from msspec.phagen.fortran.libphagen import main as do_phagen
 from msspec.spec.fortran import _eig_mi
 from msspec.spec.fortran import _eig_pw
+from msspec.spec.fortran import _eig_ar
 from msspec.spec.fortran import _phd_mi_noso_nosp_nosym
 from msspec.spec.fortran import _phd_se_noso_nosp_nosym
+from msspec.spec.fortran import _phd_ce_noso_nosp_nosym
 from msspec.spec.fortran import _comp_curves
 from msspec.utils import get_atom_index
 
@@ -404,6 +407,8 @@ class _MSCALCULATOR(Calculator):
                     do_spec = _phd_se_noso_nosp_nosym.run
                 elif self.global_parameters.algorithm == 'inversion':
                     do_spec = _phd_mi_noso_nosp_nosym.run
+                elif self.global_parameters.algorithm == 'correlation':
+                    do_spec = _phd_ce_noso_nosp_nosym.run
                 else:
                     LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not "
                                  "an allowed combination.".format(self.global_parameters.spectroscopy,
@@ -414,6 +419,8 @@ class _MSCALCULATOR(Calculator):
                     do_spec = _eig_mi.run
                 elif self.global_parameters.algorithm == 'power':
                     do_spec = _eig_pw.run
+                elif self.global_parameters.algorithm == 'arnoldi':
+                    do_spec = _eig_ar.run
                 else:
                     LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not "
                                  "an allowed combination.".format(self.global_parameters.spectroscopy,
@@ -745,7 +752,7 @@ class _PED(_MSCALCULATOR):
 
                 view = dset.add_view("E = {:.2f} eV".format(ke), title=title,
                                      xlabel=xlabel, ylabel=ylabel,
-                                     projection='stereo', colorbar=True, autoscale=True)
+                                     projection='stereo', colorbar=True, autoscale=False)
                 view.select('theta', 'phi', 'cross_section')
 
 
@@ -982,8 +989,8 @@ class _EIG(_MSCALCULATOR):
         _MSCALCULATOR.__init__(self, spectroscopy='EIG', algorithm=algorithm,
                                polarization=polarization, dichroism=dichroism,
                                spinpol=spinpol, folder=folder, txt=txt)
-        if algorithm not in ('inversion', 'power'):
+        if algorithm not in ('inversion', 'power', 'arnoldi'):
-            LOGGER.error("Only the 'inversion' or the 'power' algorithms "
+            LOGGER.error("Only the 'inversion', 'power' or 'arnoldi' algorithms "
                          "are supported in EIG spectroscopy mode")
             exit(1)
         self.iodata = iodata.Data('EIG Simulation')
@@ -1132,7 +1139,7 @@ class RFACTOR(object):
         for i in range(noif):
             X, Y = args[2*i], args[2*i+1]
             fname = os.path.join('calc',
-                                 f'calculation{self.stack_count:d}.txt')
+                                 'calculation{:d}.txt'.format(self.stack_count))
             # And save to the working space
             np.savetxt(fname, np.transpose([X, Y]))
             self.stack_count += 1
@@ -1140,7 +1147,7 @@ class RFACTOR(object):
         # Update the list of input calculation files
         self._params.calc_filename = []
         for i in range(self.stack_count):
-            fname = os.path.join('calc', f'calculation{i:d}.txt')
+            fname = os.path.join('calc', 'calculation{:d}.txt'.format(i))
             self._params.calc_filename.append(fname)
 
         # Write the input file
@@ -1235,23 +1242,23 @@ class RFACTOR(object):
             dset_values.x, dset_values.yref = exp_data.T
             # Append the calculated values
             ycalc = calc_data[:,1]
-            dset_values.add_columns(**{f"calc{i:d}": ycalc})
+            dset_values.add_columns(**{"calc{:d}".format(i): ycalc})
-            dset_rfc.add_columns(**{f'variable_set{i:d}': rfc})
+            dset_rfc.add_columns(**{'variable_set{:d}'.format(i): rfc})
 
         # Plot the curves
         view_values.select('x', 'yref', legend='Reference values')
         title = ''
         for k,v in self.best_values.items():
-            title += f'{k}={v} '
+            title += '{}={} '.format(k, v)
-        view_values.select('x', f"calc{self.index:d}",
+        view_values.select('x', "calc{:d}".format(self.index),
                            legend="Best calculated values")
         view_values.set_plot_options(title=title)
 
         view_results.select('counts')
 
         for i in range(self.stack_count):
-            view_rfc.select('rfactor_number', f'variable_set{i:d}',
+            view_rfc.select('rfactor_number', 'variable_set{:d}'.format(i),
-                            legend=f"variables set #{i:d}")
+                            legend="variables set #{:d}".format(i))
         # Save the parameters
         for p in self.get_parameters():
             bundle = {'group': str(p.group),

src/msspec/cli.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -18,8 +19,8 @@
 # along with msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/cli.py
-# Last modified: jeu. 04 juin 2020 16:54:12
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
-# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 import sys

src/msspec/create_tests_results.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -18,8 +19,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/create_tests_results.py
-# Last modified: ven. 10 avril 2020 17:29:16
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
-# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 from msspec.tests import create_tests_results

src/msspec/data/__init__.py

@@ -1,5 +1,24 @@
-# -*- encoding: utf-8 -*-
+#!/usr/bin/env python
-# vim: set fdm=indent ts=4 sw=4 sts=4 et ai tw=80 cc=+0 mouse=a nu : #
+# coding: utf-8
+#
+# Copyright © 2016-2020 - Rennes Physics Institute
+#
+# This file is part of msspec.
+#
+# msspec is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+# msspec is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# You should have received a copy of the GNU General Public License
+# along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Source file  : src/msspec/data/__init__.py
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 from .electron_be import electron_be

src/msspec/data/electron_be.py

@@ -1,4 +1,24 @@
+#!/usr/bin/env python
 # coding: utf-8
+#
+# Copyright © 2016-2020 - Rennes Physics Institute
+#
+# This file is part of msspec.
+#
+# msspec is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+# msspec is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# You should have received a copy of the GNU General Public License
+# along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Source file  : src/msspec/data/electron_be.py
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 """
 Module electron_be

src/msspec/iodata.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -16,8 +17,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/iodata.py
-# Last modified: ven. 10 avril 2020 17:23:11
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
-# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 """
@@ -442,24 +443,24 @@ class DataSet(object):
             for k, v in parameters.items():
                 if k == 'Cluster':
                     continue
-                s += f"# {k}:\n"
+                s += "# {}:\n".format(k)
                 if not(isinstance(v, list)):
                     v = [v,]
                 for p in v:
-                    s += f"#       {p['name']} = {p['value']} {p['unit']}\n"
+                    s += "#       {} = {} {}\n".format(p['name'], p['value'], p['unit'])
             return s
 
         colnames = self.columns()
         with open(filename, mode) as fd:
             # write the date and time of export
             now = datetime.now()
-            fd.write(f"# Data exported on {now}\n")
+            fd.write("# Data exported on {}\n".format(now))
             fd.write(rule)
 
             # Append notes
             fd.write("# NOTES:\n")
             for line in self.notes.split('\n'):
-                fd.write(f"#    {line}\n")
+                fd.write("#    {}\n".format(line))
             fd.write(rule)
 
             # Append parameters

src/msspec/iodata_gi.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -15,9 +16,9 @@
 # You should have received a copy of the GNU General Public License
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Source file  : src/msspec/iodata.py
+# Source file  : src/msspec/iodata_gi.py
-# Last modified: ven. 10 avril 2020 17:23:11
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
-# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 """
@@ -491,24 +492,24 @@ class DataSet(object):
             for k, v in parameters.items():
                 if k == 'Cluster':
                     continue
-                s += f"# {k}:\n"
+                s += "# {}:\n".format(k)
                 if not(isinstance(v, list)):
                     v = [v,]
                 for p in v:
-                    s += f"#       {p['name']} = {p['value']} {p['unit']}\n"
+                    s += "#       {} = {} {}\n".format(p['name'], p['value'], p['unit'])
             return s
 
         colnames = self.columns()
         with open(filename, mode) as fd:
             # write the date and time of export
             now = datetime.now()
-            fd.write(f"# Data exported on {now}\n")
+            fd.write("# Data exported on {}\n".format(now))
             fd.write(rule)
 
             # Append notes
             fd.write("# NOTES:\n")
             for line in self.notes.split('\n'):
-                fd.write(f"#    {line}\n")
+                fd.write("#    {}\n".format(line))
             fd.write(rule)
 
             # Append parameters

src/msspec/iodata_wx.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -15,9 +16,9 @@
 # You should have received a copy of the GNU General Public License
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Source file  : src/msspec/iodata.py
+# Source file  : src/msspec/iodata_wx.py
-# Last modified: ven. 10 avril 2020 17:23:11
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
-# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 """
@@ -437,24 +438,24 @@ class DataSet(object):
             for k, v in parameters.items():
                 if k == 'Cluster':
                     continue
-                s += f"# {k}:\n"
+                s += "# {}:\n".format(k)
                 if not(isinstance(v, list)):
                     v = [v,]
                 for p in v:
-                    s += f"#       {p['name']} = {p['value']} {p['unit']}\n"
+                    s += "#       {} = {} {}\n".format(p['name'], p['value'], p['unit'])
             return s
 
         colnames = self.columns()
         with open(filename, mode) as fd:
             # write the date and time of export
             now = datetime.now()
-            fd.write(f"# Data exported on {now}\n")
+            fd.write("# Data exported on {}\n".format(now))
             fd.write(rule)
 
             # Append notes
             fd.write("# NOTES:\n")
             for line in self.notes.split('\n'):
-                fd.write(f"#    {line}\n")
+                fd.write("#    {}\n".format(line))
             fd.write(rule)
 
             # Append parameters

src/msspec/looper.py

@@ -1,6 +1,24 @@
-# coding: utf8
+#!/usr/bin/env python
-# -*- encoding: future_fstrings -*-
+# coding: utf-8
-# vim: set et sw=4 ts=4 nu tw=79 cc=+1:
+#
+# Copyright © 2016-2020 - Rennes Physics Institute
+#
+# This file is part of msspec.
+#
+# msspec is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+# msspec is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# You should have received a copy of the GNU General Public License
+# along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Source file  : src/msspec/looper.py
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 from collections import OrderedDict
 from functools import partial
@@ -21,7 +39,7 @@ class Variable:
         self.doc = doc
 
     def __repr__(self):
-        return f"<Variable(\'{self.name}\')>"
+        return "<Variable(\'{}\')>".format(self.name)
 
 class Sweep:
     def __init__(self, key, comments="", unit=None,

src/msspec/misc.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -18,8 +19,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/misc.py
-# Last modified: ven. 10 avril 2020 17:30:42
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
-# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 """

src/msspec/parameters.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -18,8 +19,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/parameters.py
-# Last modified: ven. 10 avril 2020 17:31:50
+# Last modified: Tue, 15 Feb 2022 15:37:28 +0100
-# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
+# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 """
@@ -769,7 +770,8 @@ class GlobalParameters(BaseParameters):
             Parameter('algorithm', types=str, allowed_values=('expansion',
                                                               'inversion',
                                                               'correlation',
-                                                              'power'),
+                                                              'power',
+                                                              'arnoldi'),
                       default='expansion', doc="""
                       You can choose the algorithm used for the computation of the scattering path operator.
 
@@ -777,6 +779,7 @@ class GlobalParameters(BaseParameters):
                         - '**expansion**', for the Rehr-Albers series expansion
                         - '**correlation**', for the correlation-expansion algorithm
                         - '**power**', for the power method approximation scheme (only for spectroscopy='EIG')
+                        - '**arnoldi**', for computing multiple eigenvalues using Arnoldi iteration (only for spectroscopy='EIG')
 
                       The series expansion algorithm is well suited for high energy since the number of terms
                       required decreases as the energy increases. The exact solution is obtained by the matrix inversion
@@ -838,6 +841,17 @@ class GlobalParameters(BaseParameters):
         self.phagen_parameters.calctype = phagen_calctype
         self.spec_parameters.calctype_spectro = spec_calctype
 
+    def bind_polarization(self, p):
+        if p.value is None:
+            ipol = 0
+        elif p.value == 'linear_qOz':
+            ipol = 1
+        elif p.value == 'linear_xOy':
+            ipol = -1
+        elif p.value == 'circular':
+            ipol = 2
+        self.spec_parameters.calctype_ipol = ipol
+
     def bind_spinpol(self, p):
         if p.value == True:
             LOGGER.error('Spin polarization is not yet enabled in the Python version.')
@@ -2012,20 +2026,20 @@ class CompCurveGeneralParameters(BaseParameters):
         value = p.allowed_values.index(p.value)
         self.compcurve_parameters.general_norm = value
         LOGGER.info("Curve Comparison: Normalization mode set to "
-                    f"\"{p.value}\"")
+                    "\"{}\"".format(p.value))
 
     def bind_rescale(self, p):
         self.compcurve_parameters.general_iscale = int(p.value)
         state = "deactivated"
         if p.value:
             state = "activated"
-        LOGGER.info(f"Curve Comparison: Rescaling of data {state}")
+        LOGGER.info("Curve Comparison: Rescaling of data {}".format(state))
 
     def bind_function(self, p):
         value = p.allowed_values.index(p.value)
         self.compcurve_parameters.general_icur = value
         LOGGER.info("Curve Comparison: Type of data used for comparison "
-                    f"set to \"{p.value}\"")
+                    "set to \"{}\"".format(p.value))
 
 
 

src/msspec/spec/fortran/Makefile

@@ -1,6 +1,6 @@
-.PHONY: all phd_se phd_mi eig_mi eig_pw comp_curve clean
+.PHONY: all phd_se phd_mi phd_ce eig_mi eig_pw eig_ar comp_curve clean
 
-all: phd_se phd_mi eig_mi eig_pw comp_curve
+all: phd_se phd_mi phd_ce eig_mi eig_pw eig_ar comp_curve
 
 phd_se:
 	@+$(MAKE) -f phd_se_noso_nosp_nosym.mk all
@@ -8,18 +8,26 @@ phd_se:
 phd_mi:
 	@+$(MAKE) -f phd_mi_noso_nosp_nosym.mk all
 
+phd_ce:
+	@+$(MAKE) -f phd_ce_noso_nosp_nosym.mk all
+
 eig_mi:
 	@+$(MAKE) -f eig_mi.mk all
 
 eig_pw:
 	@+$(MAKE) -f eig_pw.mk all
 
+eig_ar:
+	@+$(MAKE) -f eig_ar.mk all
+
 comp_curve:
 	@+$(MAKE) -f comp_curve.mk all
 
 clean::
 	@+$(MAKE) -f phd_se_noso_nosp_nosym.mk $@
 	@+$(MAKE) -f phd_mi_noso_nosp_nosym.mk $@
+	@+$(MAKE) -f phd_ce_noso_nosp_nosym.mk $@
 	@+$(MAKE) -f eig_mi.mk $@
 	@+$(MAKE) -f eig_pw.mk $@
+	@+$(MAKE) -f eig_ar.mk $@
 	@+$(MAKE) -f comp_curve.mk $@

src/msspec/spec/fortran/eig/ar/arnoldi_mod.f90

@@ -0,0 +1,250 @@
+!==============================================================================!
+ module arnoldi_mod
+!==============================================================================!
+    implicit none
+    private
+!
+    integer, parameter :: sp = kind(1.0)
+    integer, parameter :: dp = kind(1.0d0)
+    integer, parameter :: zp = kind((0.0d0,0.0d0))
+!
+    complex(zp), parameter :: one = complex(1.0d0, 0.0d0)
+    complex(zp), parameter :: zero = complex(0.0d0, 0.0d0)
+!
+! Public procedures
+!
+    public :: arnoldi_iteration 
+!
+! Private data
+!
+! ARPACK's debug common block
+!
+    integer :: logfil = 6, ndigit = -3, mcaupd = 1
+    integer :: mgetv0 = 0, msaupd = 0, msaup2 = 0, msaitr = 0, mseigt = 0,  &
+               msapps = 0, msgets = 0, mseupd = 0, mnaupd = 0, mnaup2 = 0,  &
+               mnaitr = 0, mneigh = 0, mnapps = 0, mngets = 0, mneupd = 0,  &
+               mcaup2 = 0, mcaitr = 0, mceigh = 0, mcapps = 0, mcgets = 0,  &
+               mceupd = 0
+!
+    common /debug/ logfil, ndigit, mgetv0, msaupd, msaup2, msaitr, mseigt,  &
+                   msapps, msgets, mseupd, mnaupd, mnaup2, mnaitr, mneigh,  &
+                   mnapps, mngets, mneupd, mcaupd, mcaup2, mcaitr, mceigh,  &
+                   mcapps, mcgets, mceupd
+!==============================================================================!
+ contains
+!==============================================================================!
+! Public procedures
+!==============================================================================!
+    subroutine arnoldi_iteration (ndim, a, nev, d)
+!
+! Use ARPACK routines to find a few eigenvalues (lambda) and corresponding 
+! eigenvectors (x) for the standard eigenvalue problem:
+!          
+!                        A*x = lambda*x
+! 
+! where A is a general NDIM by NDIM complex matrix
+!
+! This subroutine is based on an ARPACK test program adapted by Logan Boudet
+! as part of his M1 project in Rennes in 2022.
+! 
+    integer,     intent(in)    :: ndim
+    complex(zp), intent(in)    :: a(:,:)
+    integer,     intent(inout) :: nev           ! number of eigenvalues required
+    complex(zp), intent(out)   :: d(:)          ! vector of required eigenvalues
+!
+! Local data
+!
+    character(1), parameter :: bmat = "I"       ! standarg eigenvalue problem
+    character(2), parameter :: which = "LM"     ! find NEV eigenvalues of 
+!                                               ! largest magnitude
+!
+    integer     :: ido, ierr, info, j, lworkl, nconv, ncv 
+    integer     :: iparam(11)
+    integer     :: ipntr(14)
+    real(dp)    :: tol
+    complex(zp) :: sigma
+    logical     :: rvec
+!
+    real(dp),    allocatable :: rd(:,:)
+    real(dp),    allocatable :: rwork(:)
+    complex(zp), allocatable :: ax(:)
+    complex(zp), allocatable :: resid(:)
+    complex(zp), allocatable :: v(:,:)
+    complex(zp), allocatable :: workd(:)
+    complex(zp), allocatable :: workev(:)
+    complex(zp), allocatable :: workl(:)
+    logical,     allocatable :: select(:)
+!
+! External BLAS/LAPACK functions used
+!
+    real(dp), external :: dznrm2, dlapy2 
+!
+!
+    write(6,*)
+    write(6,*) "----------------- BEGIN OF ARNOLDI ITERATION -----------------"
+!
+! NCV is is the largest number of basis vectors that will be used in the 
+! Implicitly Restarted Arnoldi Process. Work per major iteration is 
+! proportional to NDIM*NCV*NCV.
+! 
+! Note: we must have NCV >= NEV + 2, and preferably NCV >= 2*NEV
+!
+!    ncv = max(ceiling(1.125*nev + 15), nev+3)
+    ncv = 2*nev
+!
+    iparam(11) = 0
+    ipntr(14) = 0
+!
+! stopping criteria; machine precision is used if tol <= 0
+!
+    tol = 0.0_dp
+!
+! Algorithm mode
+!
+    iparam(1) = 1           ! exact shift strategy              
+    iparam(3) = 300         ! maximum number of iterations
+    iparam(7) = 1           ! use mode1 of ZNAUPD
+!
+    lworkl = ncv*(3*ncv + 5)
+    allocate(rwork(ncv))
+    allocate(resid(ndim), v(ndim,ncv))
+    allocate(workd(3*ndim), workev(2*ncv), workl(lworkl))
+    allocate(select(ncv))
+!
+! IDO is the reverse communication parameter used to determine action to be taken
+! on return from ZNAUPD. Its initial value must be 0.
+!
+    ido = 0
+!
+! On entry, INFO == 0 instructs ZNAUPD to use a random starting vector.
+! To specify a particular starting vector, set INFO to a non-zero value.
+! The required startng vector should then be supplied in RESID.
+!    
+    info = 0
+!
+! Main reverse communication loop
+!
+    do
+        call znaupd(ido, bmat, ndim, which, nev, tol, resid, ncv, v, ndim,  &
+                    iparam, ipntr, workd, workl, lworkl, rwork, info)
+
+        if (abs(ido) /= 1) exit
+!
+! Matrix-vector multiplication y = A*x
+! Initial vector x is stored starting at  workd(ipntr(1))
+! The result y should be stored in workd(ipntr(2))
+!
+        call matvec(a, ndim, ndim, workd(ipntr(1)), workd(ipntr(2)))
+!
+    end do
+!
+    if (info < 0) then
+        write(6,*)
+        write(6,*) "Error: znaupd returned info = ", info
+        write(6,*) "Check the documentation of znaupd for more information"
+        write(6,*)
+        stop
+    end if
+!
+! No fatal errors, post-process with ZNEUPD to extract computed eigenvalues.
+! Eigenvectors may be also computed by setting RVEC = .TRUE.)
+!
+    if (info == 1) then
+        write(6,*)
+        write(6,*) "Maximum number of iterations reached."
+        write(6,*)
+    else if (info == 3) then
+        write(6,*)
+        write(6,*) "No shifts could be applied during implicit Arnoldi update"
+        write(6,*) "Try increasing NCV."
+        write(6,*)
+    end if      
+!
+    rvec = .false.
+!
+    call zneupd(rvec, 'A', select, d, v, ndim, sigma, workev, bmat, ndim,    &
+                which, nev, tol, resid, ncv, v, ndim, iparam, ipntr, workd,  &
+                workl, lworkl, rwork, ierr)
+!
+    if (ierr /= 0) then
+        write(6,*)
+        write(6,*) "Error: zneupd returned ierr = ", ierr
+        write(6,*) "Check the documentation of zneupd for more information"
+        write(6,*)
+        stop
+    end if
+!
+! Eigenvalues are returned in the one dimensional array D and if RVEC == .TRUE.
+! the corresponding eigenvectors are returned in the first NCONV == IPARAM(5) 
+! columns of the two dimensional array V 
+! 
+    nconv = iparam(5)
+!
+    if (rvec) then
+!
+! Compute the residual norm ||  A*x - lambda*x || for the NCONV accurately 
+! computed eigenvalues and eigenvectors
+!
+        allocate(rd(ncv,3), ax(ndim))
+!
+        do j = 1, nconv
+            call matvec(a, ndim, ndim, v(:,j), ax)
+            call zaxpy(ndim, -d(j), v(:,j), 1, ax, 1)
+            rd(j,1) = real(d(j), dp)
+            rd(j,2) = aimag(d(j))
+            rd(j,3) = dznrm2(ndim, ax, 1) / dlapy2(rd(j,1), rd(j,2))
+        end do
+!
+        call dmout(6, nconv, 3, rd, 2*nev, -6,  &
+                   "Ritz values (Real, Imag) and relative residuals")
+!
+        deallocate(rd, ax)
+!
+    end if
+!
+    write(6,*)
+    write(6,*) "SUMMARY"
+    write(6,*) "======="
+    write(6,*)
+    write(6,*) "Size of the matrix is ", ndim
+    write(6,*) "The number of Ritz values requested is ", nev
+    write(6,*) "The number of Arnoldi vectors generated (NCV) is ", ncv
+    write(6,*) "Portion of the spectrum: ", which
+    write(6,*) "The number of converged Ritz values is ", nconv 
+    write(6,*) "The number of implicit Arnoldi update iterations taken is ", iparam(3)
+    write(6,*) "The number of OP*x is ", iparam(9)
+    write(6,*) "The convergence criterion is ", tol
+    write(6,*) 
+!
+    nev = nconv
+!   
+    write(6,*) "------------------ END OF ARNOLDI ITERATION ------------------"
+!
+    deallocate(rwork)
+    deallocate(resid, v)
+    deallocate(workd, workev, workl)
+    deallocate(select)
+!
+    return
+    end subroutine arnoldi_iteration
+!==============================================================================!
+! Private procedures
+!==============================================================================!
+    subroutine matvec (a, n, lda, x, y)
+!
+! Compute the matrix-vector product a*x, storing the result in y
+!
+    complex(zp), intent(in)  :: a(:,:)
+    integer,     intent(in)  :: n
+    integer,     intent(in)  :: lda
+    complex(zp), intent(in)  :: x(*)
+    complex(zp), intent(out) :: y(*)
+!
+!
+    call zgemv('n', n, n, one, a, lda, x, 1, zero, y, 1)
+!
+    return
+    end subroutine matvec
+!==============================================================================!
+ end module arnoldi_mod
+!==============================================================================!

src/msspec/spec/fortran/eig/ar/eig_mat_ar.f90

@@ -0,0 +1,291 @@
+!==============================================================================!
+ subroutine eig_mat_ar (je, e_kin)
+!==============================================================================!
+! This subroutine stores the G_o T kernel matrix and computes a subset of the
+! the eigenvalues with largest magnitude using Arnoldi iteration methods from 
+! ARPACK
+!
+    use dim_mod,      only: n_gaunt, nl_m
+    use coor_mod,     only: natyp, ncorr, n_prot, sym_at
+    use outfiles_mod, only: outfile2
+    use outunits_mod, only: iuo1, iuo2
+    use trans_mod,    only: lmax, vk, tl
+!
+    use arnoldi_mod,  only: arnoldi_iteration
+!
+    implicit none
+!
+    integer, parameter :: sp = kind(1.0)
+    integer, parameter :: dp = kind(1.0d0)
+    integer, parameter :: cp = kind((0.0,0.0))
+    integer, parameter :: zp = kind((0.0d0,0.0d0))
+!
+! Subroutine arguments
+!
+    integer,  intent(in) :: je
+    real(sp), intent(in) :: e_kin
+!
+! Local data
+!
+    integer,     parameter :: ibess = 3
+    integer,     parameter :: nprint = 10
+    real(dp),    parameter :: pi = acos(-1.0_dp)
+    real(dp),    parameter :: small = 0.0001_dp
+    complex(zp), parameter :: ic = complex(0.0d0,1.0d0)
+    complex(zp), parameter :: zeroc = complex(0.0d0,0.0d0)
+    complex(zp), parameter :: four_pi_i = 4.0_dp*pi*ic
+!
+    integer     :: jatl, jlin, jtyp, jnum, lj, lmj, mj, nbtypj
+    integer     :: katl, klin, ktyp, knum, lk, lmk, mk, nbtypk
+    integer     :: j, jp, l, lin, l_max, l_min, m
+    integer     :: ndim, n_dot, n_eig, nev, nfin, nltwo, npr, n_xmax
+    real(sp)    :: eig, xj, yj, zj, xmax_l
+    real(dp)    :: attkj, xkj, ykj, zkj, rkj, zdkj, krkj
+    complex(zp) :: expkj, sum_l, tlk
+!
+    integer, save :: iout2, iout3
+!
+    real(sp),     allocatable :: w1(:), w2(:)
+    real(dp),     allocatable :: gnt(:)
+    complex(zp),  allocatable :: hl1(:), sm(:,:), ylm(:,:), w(:)
+    character(:), allocatable :: outfile, path
+!
+!
+    if (je == 1) then
+!
+! Name of second output file where eigenvalues will be written
+!
+        n_dot = index(outfile2, '.')
+        outfile = outfile2(1:n_dot)//'egv'
+        path = outfile2(1:n_dot)//'pth'
+        open(newunit=iout2, file=outfile, status='unknown')
+        open(newunit=iout3, file=path, status='unknown')
+!
+    end if
+!
+! Construction of the multiple scattering kernel matrix G_o T.
+! Elements are stored using a linear index LINJ representing (J,LJ)
+!
+! First compute Go T array dimension
+!
+    jlin = 0
+    do jtyp = 1, n_prot
+        nbtypj = natyp(jtyp)
+        lmj = lmax(jtyp,je)
+        do jnum = 1, nbtypj
+            do lj = 0, lmj
+                do mj = -lj, lj
+                    jlin = jlin + 1
+                end do
+            end do
+        end do
+    end do
+!
+    ndim = jlin
+    write(6,*) "GoT matrix SM has dimension ", ndim
+!
+    allocate(sm(ndim,ndim))
+    sm = zeroc
+!
+    nltwo = 2*nl_m
+    allocate(ylm(0:nltwo, -nltwo:nltwo))
+    ylm = zeroc
+!
+    allocate(hl1(0:nltwo))
+    hl1 = zeroc
+!
+    allocate(gnt(0:n_gaunt))
+    gnt = 0.0_dp
+!
+    jlin = 0
+    do jtyp = 1, n_prot
+        nbtypj = natyp(jtyp)
+        lmj = lmax(jtyp,je)
+        do jnum = 1, nbtypj
+            jatl = ncorr(jnum,jtyp)
+            xj = sym_at(1,jatl)
+            yj = sym_at(2,jatl)
+            zj = sym_at(3,jatl)
+            do lj = 0, lmj
+                do mj = -lj, lj
+                    jlin = jlin + 1
+!
+                    klin=0
+                    do ktyp = 1, n_prot
+                        nbtypk = natyp(ktyp)
+                        lmk = lmax(ktyp,je)
+                        do knum = 1, nbtypk
+                            katl = ncorr(knum,ktyp)
+!
+                            if (katl /= jatl) then
+                                xkj = real(sym_at(1,katl) - xj, dp)
+                                ykj = real(sym_at(2,katl) - yj, dp)
+                                zkj = real(sym_at(3,katl) - zj, dp)
+                                rkj = sqrt(xkj*xkj + ykj*ykj + zkj*zkj)
+                                krkj = real(vk(je), dp)*rkj
+                                attkj = exp(-aimag(cmplx(vk(je)))*rkj)
+                                expkj = (xkj + ic*ykj)/rkj
+                                zdkj = zkj/rkj
+                                call sph_har2(2*nl_m, zdkj, expkj, ylm, lmj+lmk)
+                                call besphe2(lmj+lmk+1, ibess, krkj, hl1)
+                            end if
+!
+                            do lk = 0,lmk
+                                l_min = abs(lk-lj)
+                                l_max = lk + lj
+                                tlk = cmplx(tl(lk,1,ktyp,je))
+                                do mk = -lk, lk
+                                    klin = klin + 1
+! 
+                                    sm(klin,jlin) = zeroc
+                                    sum_l = zeroc
+                                    if (katl /= jatl) then
+                                        call gaunt2(lk, mk, lj, mj, gnt)
+                                        do l = l_min, l_max, 2
+                                            m = mj - mk
+                                            if (abs(m) <= l) then
+                                                sum_l = sum_l + (ic**l)*hl1(l)*ylm(l,m)*gnt(l)
+                                            end if
+                                        end do
+                                        sum_l = sum_l*attkj*four_pi_i
+                                    else
+                                        sum_l = zeroc
+                                    end if
+!
+                                    sm(klin,jlin) = tlk*sum_l
+!
+                                end do
+                            end do
+                        end do
+                    end do
+                end do
+            end do
+        end do
+    end do
+!
+    deallocate(ylm, hl1, gnt)
+!
+! Compute subset of eigenvalues of SM using ARPACK
+! 
+! NEV is the number of eigenvalues required, set to 2% of NDIM
+!
+    nev = ceiling(0.02*ndim)
+    allocate(w(nev))
+!
+    call arnoldi_iteration(ndim, sm, nev, w)
+!
+    deallocate(sm)
+!
+! Save results to filestream for easy access from python
+!
+    call save_eigenvalues(w, nev, e_kin)
+!
+! Write results to OUTFILE on unit IOUT2
+!
+    write(iout2,75)
+    write(iout2,110)
+    write(iout2,80) e_kin
+    write(iout2,110)
+    write(iout2,75)
+    write(iout2,105)
+    write(iout2,75)
+!
+    allocate(w1(nev), w2(nev))
+!
+    n_eig = 0
+    xmax_l = 0.0
+    n_xmax = 0
+    do lin = 1, nev
+        eig = real(abs(w(lin)))
+        write(iout2,100) real(w(lin)), aimag(w(lin)), eig
+        if ((eig-xmax_l) > 0.0001) n_xmax = lin
+        xmax_l = max(xmax_l, eig)
+        w1(lin) = eig
+        if (eig > 1.000) then
+            n_eig = n_eig + 1
+        end if
+    end do
+!
+    write(iout2,75)
+    write(iout2,85) xmax_l
+    write(iout2,90) n_xmax
+    write(iout2,95) w(n_xmax)
+    write(iout2,75)
+!
+! Summarize results in main output file
+!
+    call ordre(nev, w1, nfin, w2)
+!
+    write(iuo1,5)
+    write(iuo1,10)
+    write(iuo1,10)
+    write(iuo1,15) w2(1)
+    write(iuo1,20) w2(nfin)
+    write(iuo1,10)
+    write(iuo1,10)
+!
+    if (n_eig >= 1) then
+        if (n_eig == 1) then
+            write(iuo1,25) ndim
+        else
+            write(iuo1,30) n_eig, ndim
+        end if
+    end if
+!
+    write(iuo1,65) n_xmax
+    write(iuo1,70) w(n_xmax)
+    write(iuo1,10)
+    write(iout3,100) real(w(n_xmax)), aimag(w(n_xmax))
+!
+    npr = nprint/5
+    write(iuo1,10)
+    write(iuo1,10)
+    write(iuo1,35) 5*npr
+    write(iuo1,10)
+    do jp = 0, npr-1
+        j = 5*jp
+        write(iuo1,40) w2(j+1), w2(j+2), w2(j+3), w2(j+4), w2(j+5)
+    enddo
+    write(iuo1,10)
+    write(iuo1,10)
+    write(iuo1,45) w2(1)
+    write(iuo2,*) e_kin, w2(1)
+    if (n_eig == 0) then
+        write(iuo1,50)
+    else
+        write(iuo1,55)
+    end if
+    write(iuo1,10)
+    write(iuo1,10)
+    write(iuo1,60)
+!
+    deallocate(w, w1, w2)
+!
+    return
+!
+  5 format(/,11X,'-----------------  EIGENVALUE ANALYSIS ','-----------------')
+ 10 format(11X,'-',54X,'-')
+ 15 format(11X,'-',14X,'MAXIMUM MODULUS : ',F9.6,13X,'-')
+ 20 format(11X,'-',14X,'MINIMUM MODULUS : ',F9.6,13X,'-')
+ 25 format(11X,'-',6X,'1 EIGENVALUE IS > 1 OF A TOTAL OF ',I8,6X,'-')
+ 30 format(11X,'-',4X,I5,' EIGENVALUES ARE > 1 OF A TOTAL OF ',I8,2X,'-')
+ 35 format(11X,'-',11X,'THE ',I3,' LARGEST EIGENVALUES ARE :',11X,'-')
+ 40 format(11X,'-',6X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,5X,'-')
+ 45 format(11X,'-',4X,'SPECTRAL RADIUS OF THE KERNEL MATRIX : ',F8.5,3X,'-')
+ 50 format(11X,'-',5X,'---> THE MULTIPLE SCATTERING SERIES ','CONVERGES',4X,'-')
+ 55 format(11X,'-',10X,'---> NO CONVERGENCE OF THE MULTIPLE',9X,'-',/,11X,'-',   &
+        18X,'SCATTERING SERIES',19X,'-')
+ 60 format(11X,'----------------------------------------','----------------',/)
+ 65 format(11X,'-',5X,'    LABEL OF LARGEST EIGENVALUE   : ',I5,8X,'-')
+ 70 format(11X,'-',5X,'    LARGEST EIGENVALUE  : ','(',F6.3,',',F6.3,')',8X,'-')
+ 75 format('              ')
+ 80 format('    KINETIC ENERGY : ',F7.2,' eV')
+ 85 format('    LARGEST MODULUS OF EIGENVALUE : ',F6.3)
+ 90 format('    LABEL OF LARGEST EIGENVALUE   : ',I5)
+ 95 format('    LARGEST EIGENVALUE            :  (',F6.3,',',F6.3,')')
+100 format(5X,F9.5,2X,F9.5,2X,F9.5)
+105 format(7X,'EIGENVALUES :',3X,'MODULUS :')
+110 format(2X,'-------------------------------')
+!==============================================================================!
+ end subroutine eig_mat_ar
+!==============================================================================!

src/msspec/spec/fortran/eig/ar/eigdif_ar.f

@@ -0,0 +1,103 @@
+C
+C=======================================================================
+C
+      SUBROUTINE EIGDIF_AR
+C
+C   This subroutine computes some of the eigenvalues of the
+C     multiple scattering matrix using Arnoldi iteration as
+C     implemented in ARPACK.
+C
+C                                      Last modified : 21 June 2023
+C
+C      INCLUDE 'spec.inc'
+      USE DIM_MOD
+      USE CONVTYP_MOD
+      USE COOR_MOD, NTCLU => NATCLU, NTP => NATYP
+      USE DEBWAL_MOD
+      USE EIGEN_MOD, NE => NE_EIG, E0 => E0_EIG, EFIN => EFIN_EIG
+      USE OUTFILES_MOD
+      USE OUTUNITS_MOD
+      USE RESEAU_MOD
+      USE TESTS_MOD
+      USE TRANS_MOD
+      USE VALIN_MOD, E1 => E0, PHLUM => PHILUM
+C
+      COMPLEX IC,ONEC
+      COMPLEX TLT(0:NT_M,4,NATM,NE_M)
+C
+C
+      DATA CONV /0.512314/
+C
+      IC=(0.,1.)
+      ONEC=(1.,0.)
+C
+      OPEN(UNIT=IUO2, FILE=OUTFILE2, STATUS='UNKNOWN')
+C
+C   Loop over the energies
+C
+      DO JE=1,NE
+        IF(NE.GT.1) THEN
+          ECIN=E0+FLOAT(JE-1)*(EFIN-E0)/FLOAT(NE-1)
+        ELSEIF(NE.EQ.1) THEN
+          ECIN=E0
+        ENDIF
+        CALL LPM(ECIN,XLPM,*6)
+        XLPM1=XLPM/A
+        IF(IPRINT.GT.0) WRITE(IUO1,56) A,XLPM1
+        IF(ITL.EQ.0) THEN
+          VK(JE)=SQRT(ECIN+VINT)*CONV*A*ONEC
+          VK2(JE)=CABS(VK(JE)*VK(JE))
+        ENDIF
+        GAMMA=1./(2.*XLPM1)
+        IF(IPOTC.EQ.0) THEN
+          VK(JE)=VK(JE)+IC*GAMMA
+        ENDIF
+        IF(I_MFP.EQ.0) THEN
+          VK(JE)=CMPLX(REAL(VK(JE)))
+          VK2(JE)=CABS(VK(JE)*VK(JE))
+        ENDIF
+        IF(I_VIB.EQ.1) THEN
+          IF(IDCM.GE.1) WRITE(IUO1,22)
+          DO JAT=1,N_PROT
+            IF(IDCM.EQ.0) THEN
+              XK2UJ2=VK2(JE)*UJ2(JAT)
+            ELSE
+              XK2UJ2=VK2(JE)*UJ_SQ(JAT)
+              WRITE(IUO1,23) JAT,UJ_SQ(JAT)*A*A
+            ENDIF
+            CALL DWSPH(JAT,JE,XK2UJ2,TLT,I_VIB)
+            DO LAT=0,LMAX(JAT,JE)
+              TL(LAT,1,JAT,JE)=TLT(LAT,1,JAT,JE)
+            ENDDO
+          ENDDO
+       ENDIF
+C
+C     Eigenvalue calculation
+C
+ckmd      IF(I_PWM.EQ.0) THEN
+ckmd        CALL EIG_MAT_MS(JE,ECIN)
+ckmd      ELSE
+ckmd        CALL SPEC_RAD_POWER(JE,ECIN)
+ckmd      ENDIF
+C
+      call eig_mat_ar(je, ecin)
+C
+C
+C     End of the loop on the energy
+C
+      ENDDO
+      GOTO 7
+C
+   6  WRITE(IUO1,55)
+C
+  22  FORMAT(16X,'INTERNAL CALCULATION OF MEAN SQUARE DISPLACEMENTS',/,
+     &25X,' BY DEBYE UNCORRELATED MODEL:',/)
+  23  FORMAT(21X,'ATOM TYPE ',I5,'  MSD = ',F8.6,' ANG**2')
+  55  FORMAT(///,12X,' <<<<<<<<<<  THIS VALUE OF ILPM IS NOT',
+     &'AVAILABLE  >>>>>>>>>>')
+  56  FORMAT(4X,'LATTICE PARAMETER A = ',F6.3,' ANGSTROEMS',4X,
+     *'MEAN FREE PATH = ',F6.3,' * A',//)
+C
+   7  RETURN
+C
+      END

src/msspec/spec/fortran/eig/ar/main.f

@@ -0,0 +1,22 @@
+      SUBROUTINE RUN(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_, 
+     &  NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_, 
+     &  NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_, 
+     &  N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_)
+
+        USE DIM_MOD
+        IMPLICIT INTEGER (A-Z)
+CF2PY   INTEGER, INTENT(IN,COPY) :: NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_
+CF2PY   INTEGER, INTENT(IN,COPY) :: NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_ 
+CF2PY   INTEGER, INTENT(IN,COPY) :: NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_
+CF2PY   INTEGER, INTENT(IN,COPY) :: N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_
+
+        CALL ALLOCATION(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_, 
+     &  NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_, 
+     &  NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_, 
+     &  N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_)
+
+        CALL DO_MAIN()
+        CALL CLOSE_ALL_FILES()
+
+      END SUBROUTINE
+

src/msspec/spec/fortran/eig/common/eig_mat_ms.f

@@ -178,6 +178,10 @@ CKMD          WRITE(IUO1,*) '                        '
 CKMD          WRITE(IUO1,*) '    --->  WORK(1),INFO =',WORK(1),INFO
 CKMD          WRITE(IUO1,*) '                        '
 CKMD      ENDIF
+C
+CKMD Save eigenvalues to unformatted stream file eigenvalues.dat
+C
+      call save_eigenvalues(w, jlin, e_kin)
 C
       N_EIG=0
 C

src/msspec/spec/fortran/eig/common/save_eigenvalues.f

@@ -0,0 +1,37 @@
+c
+c=======================================================================
+c
+      subroutine save_eigenvalues (evalues, n, ke)
+c
+      implicit none
+c
+      integer,    intent(in) :: n
+      real,       intent(in) :: ke
+      complex*16, intent(in) :: evalues(n)
+c
+c Local variables
+c
+      integer :: io
+      logical :: exists
+c
+c
+      inquire(file='eigenvalues.dat', exist=exists)
+c
+      if (exists) then
+         open(newunit=io, file='eigenvalues.dat', status='old',
+     +        form='unformatted', access='stream', action='write',
+     +        position='append')
+      else
+         open(newunit=io, file='eigenvalues.dat', status='new',
+     +        form='unformatted', access='stream', action='write')
+      end if
+c
+      write(io) ke, n, evalues(1:n)
+c
+      close(io)
+c
+      return
+      end subroutine save_eigenvalues
+c
+c=======================================================================
+c

src/msspec/spec/fortran/eig_ar.mk

@@ -0,0 +1,14 @@
+memalloc_src               := memalloc/dim_mod.f memalloc/modules.f memalloc/allocation.f 
+cluster_gen_src            := $(wildcard cluster_gen/*.f)
+common_sub_src             := $(wildcard common_sub/*.f)
+renormalization_src        := $(wildcard renormalization/*.f)
+#eig_common_src             := $(wildcard eig/common/*.f)
+eig_common_src             := $(filter-out eig/common/lapack_eig.f, $(wildcard eig/common/*.f))
+eig_ar_src                 := $(wildcard eig/ar/*.f)
+eig_ar_src_f90             := $(wildcard eig/ar/*.f90)
+
+SRCS   = $(memalloc_src) $(cluster_gen_src) $(common_sub_src) $(renormalization_src) $(eig_common_src) $(eig_ar_src_f90) $(eig_ar_src)
+MAIN_F = eig/ar/main.f
+SO     = _eig_ar.so
+
+include ../../../options.mk

src/msspec/spec/fortran/eig_mi.mk

@@ -2,7 +2,8 @@ memalloc_src               := memalloc/dim_mod.f memalloc/modules.f memalloc/all
 cluster_gen_src            := $(wildcard cluster_gen/*.f)
 common_sub_src             := $(wildcard common_sub/*.f)
 renormalization_src        := $(wildcard renormalization/*.f)
-eig_common_src             := $(wildcard eig/common/*.f)
+#eig_common_src             := $(wildcard eig/common/*.f)
+eig_common_src             := $(filter-out eig/common/lapack_eig.f, $(wildcard eig/common/*.f))
 eig_mi_src                 := $(wildcard eig/mi/*.f)
 
 SRCS   = $(memalloc_src) $(cluster_gen_src) $(common_sub_src) $(renormalization_src) $(eig_common_src) $(eig_mi_src)

src/msspec/spec/fortran/memalloc/allocation.f

@@ -25,6 +25,9 @@
         USE OUTUNITS_MOD
         USE PARCAL_MOD
         USE PARCAL_A_MOD
+        USE CORREXP_MOD
+        USE GAUNT_C_MOD
+        USE Q_ARRAY_MOD
         USE RELADS_MOD
         USE RELAX_MOD
         USE RESEAU_MOD
@@ -136,6 +139,7 @@
         CALL ALLOC_OUTUNITS()
         CALL ALLOC_PARCAL()
         CALL ALLOC_PARCAL_A()
+        CALL ALLOC_Q_ARRAY()
         CALL ALLOC_RELADS()
         CALL ALLOC_RELAX()
         CALL ALLOC_RENORM()
@@ -173,6 +177,7 @@
         CALL ALLOC_C_G()
         CALL ALLOC_C_G_A()
         CALL ALLOC_C_G_M()
+        CALL ALLOC_CORREXP()
         CALL ALLOC_DEXPFAC2()
         CALL ALLOC_DFACTSQ()
         CALL ALLOC_EIGEN()
@@ -186,6 +191,7 @@
         CALL ALLOC_SPECTRUM()
         CALL ALLOC_DIRECT()
         CALL ALLOC_DIRECT_A()
+        CALL ALLOC_GAUNT_C()
         CALL ALLOC_PATH()
         CALL ALLOC_ROT()
         CALL ALLOC_ROT_CUB()

src/msspec/spec/fortran/memalloc/dim_mod.f

@@ -34,6 +34,7 @@ C       ===============================================================
         INTEGER NCG_M
         INTEGER N_BESS, N_GAUNT
         INTEGER NLTWO
+        INTEGER NLMM
 C       ===============================================================
       CONTAINS
         SUBROUTINE INIT_DIM()
@@ -60,9 +61,10 @@ C       ===============================================================
 
 C          N_BESS=100*NL_M
 C          N_GAUNT=5*NL_M
-          N_BESS=200*NL_M
+          N_BESS=300*NL_M
           N_GAUNT=10*NL_M
 
           NLTWO=2*NL_M
+          NLMM=LINMAX*NGR_M
         END SUBROUTINE INIT_DIM
       END MODULE DIM_MOD

src/msspec/spec/fortran/memalloc/modules.f

@@ -192,6 +192,20 @@ C=======================================================================
         END SUBROUTINE ALLOC_COOR
       END MODULE COOR_MOD
 
+C=======================================================================
+      MODULE CORREXP_MOD
+        IMPLICIT NONE
+        COMPLEX*16, ALLOCATABLE, DIMENSION(:,:) :: A
+      CONTAINS
+        SUBROUTINE ALLOC_CORREXP()
+          USE DIM_MOD
+          IF (ALLOCATED(A)) THEN
+            DEALLOCATE(A)
+          ENDIF
+          ALLOCATE(A(NLMM,NLMM))
+        END SUBROUTINE ALLOC_CORREXP
+      END MODULE CORREXP_MOD
+
 C=======================================================================
       MODULE DEBWAL_MOD
         IMPLICIT NONE
@@ -417,6 +431,20 @@ C=======================================================================
         END SUBROUTINE ALLOC_PARCAL_A
       END MODULE PARCAL_A_MOD
 
+C=======================================================================
+      MODULE Q_ARRAY_MOD
+        IMPLICIT NONE
+        REAL, ALLOCATABLE, DIMENSION(:) :: Q
+      CONTAINS
+        SUBROUTINE ALLOC_Q_ARRAY()
+          USE DIM_MOD
+          IF (ALLOCATED(Q)) THEN
+            DEALLOCATE(Q)
+          ENDIF
+          ALLOCATE(Q(NGR_M))
+        END SUBROUTINE ALLOC_Q_ARRAY
+      END MODULE Q_ARRAY_MOD
+
 C=======================================================================
       MODULE RELADS_MOD
         IMPLICIT NONE
@@ -778,6 +806,20 @@ C=======================================================================
         END SUBROUTINE ALLOC_DEXPFAC
       END MODULE DEXPFAC_MOD
 
+C=======================================================================
+      MODULE GAUNT_C_MOD
+        IMPLICIT NONE
+        REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: GNT
+      CONTAINS
+        SUBROUTINE ALLOC_GAUNT_C()
+          USE DIM_MOD
+          IF (ALLOCATED(GNT)) THEN
+            DEALLOCATE(GNT)
+          ENDIF
+          ALLOCATE(GNT(0:N_GAUNT,LINMAX,LINMAX))
+        END SUBROUTINE ALLOC_GAUNT_C
+      END MODULE GAUNT_C_MOD
+
 C=======================================================================
       MODULE LOGAMAD_MOD
         IMPLICIT NONE

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym.mk

@@ -0,0 +1,11 @@
+memalloc_src               := memalloc/dim_mod.f memalloc/modules.f memalloc/allocation.f 
+cluster_gen_src            := $(wildcard cluster_gen/*.f)
+common_sub_src             := $(wildcard common_sub/*.f)
+renormalization_src        := $(wildcard renormalization/*.f)
+phd_ce_noso_nosp_nosym_src := $(filter-out phd_ce_noso_nosp_nosym/lapack_axb.f, $(wildcard phd_ce_noso_nosp_nosym/*.f))
+
+SRCS   = $(memalloc_src) $(cluster_gen_src) $(common_sub_src) $(renormalization_src) $(phd_ce_noso_nosp_nosym_src)
+MAIN_F = phd_ce_noso_nosp_nosym/main.f
+SO     = _phd_ce_noso_nosp_nosym.so
+
+include ../../../options.mk

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/cmngr.f

@@ -0,0 +1,41 @@
+C
+C======================================================================
+C
+      SUBROUTINE CMNGR(NAT,NGR,CMN)
+C
+C    input   : NAT,NGR
+C    output  : CMN
+C
+C    This subroutine calculate C(NAT-N,M-N) where,
+C                    1<=M<=NGR<=NAT,1<=N<=M
+C               C(NAT-N,M-N) is stored as CMN(N,M)
+C
+C                                              H.-F. Zhao 2007
+C
+      USE DIM_MOD
+C
+      INTEGER NAT,NGR
+C
+      REAL CMN(NGR_M,NGR_M)
+C
+      IF(NGR.GT.NAT) THEN
+        WRITE(6,*) 'NGR is larger than NAT, which is wrong'
+        STOP
+      ENDIF
+C
+      DO M=1,NGR
+        DO N=1,NGR
+          CMN(N,M)=0.
+        ENDDO
+        CMN(M,M)=1.
+      ENDDO
+C
+      DO M=1,NGR
+        DO N=M-1,1,-1
+          CMN(N,M)=CMN(N+1,M)*FLOAT(NAT-N)/FLOAT(M-N)
+        ENDDO
+      ENDDO
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/coefpq.f

@@ -0,0 +1,46 @@
+C
+C======================================================================
+C
+      SUBROUTINE COEFPQ(NAT,NGR)
+C
+C  This subroutine computes the P(n,m) and Q(n) coefficients 
+C   involved in the correlation expansion formulation
+C
+C  Reference  : equations (2.15) and (2.16) of 
+C               H. Zhao, D. Sebilleau and Z. Wu, 
+C               J. Phys.: Condens. Matter 20, 275241 (2008)
+C
+C                                                    H.-F. Zhao 2007
+C
+      USE DIM_MOD
+      USE Q_ARRAY_MOD
+C
+      INTEGER NAT,NGR
+C
+      REAL CMN(NGR_M,NGR_M),P(NGR_M,NGR_M)
+C
+C
+      IF(NGR.GT.NAT) THEN
+        WRITE(6,*) 'NGR is larger than NAT, which is wrong'
+        STOP
+      ENDIF
+C
+      CALL CMNGR(NAT,NGR,CMN)
+C
+      DO N=1,NGR
+        P(N,N)=1.
+        Q(N)=P(N,N)
+        DO M=N+1,NGR
+          P(N,M)=0.
+          DO I=N,M-1
+            P(N,M)=P(N,M)-P(N,I)*CMN(I,M)
+          ENDDO
+          Q(N)=Q(N)+P(N,M)
+C
+        ENDDO
+C
+      ENDDO
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/corexp.f

@@ -0,0 +1,47 @@
+C
+C======================================================================
+C
+      SUBROUTINE COREXP_SAVM(JE,IGR,NGR,NLM,ITYPE,IGS,TAU)
+C
+C  This subroutine call the correlation matrices calculations
+C    for a given order IGR
+C
+C                                         H.-F. Zhao : 2007
+C
+      USE DIM_MOD
+      USE COOR_MOD
+      USE Q_ARRAY_MOD
+      USE TRANS_MOD
+C
+      INTEGER NLM(NGR_M),ITYPE(NGR_M),IGS(NGR_M)
+C
+      REAL QI
+C
+      COMPLEX*16 TAU(LINMAX,LINFMAX,NATCLU_M)
+C
+C
+      DO ITYP=1,N_PROT
+        NBTYP=NATYP(ITYP)
+        NLM(IGR)=LMAX(ITYP,JE)
+        ITYPE(IGR)=ITYP
+        DO NUM=1,NBTYP
+          IGS(IGR)=NCORR(NUM,ITYP)
+C
+          IF(IGS(IGR).GT.IGS(IGR-1)) THEN
+            QI=Q(IGR)
+            CALL MPIS(IGR,NLM,ITYPE,IGS,JE,QI,TAU)
+C
+            IGR=IGR+1
+            IF(IGR.LE.NGR) THEN
+              CALL COREXP_SAVM1(JE,IGR,NGR,NLM,ITYPE,IGS,TAU)
+            ENDIF
+            IGR=IGR-1
+C
+          ENDIF
+C
+        ENDDO
+      ENDDO
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/corexp1.f

@@ -0,0 +1,19 @@
+C
+C======================================================================
+C
+      SUBROUTINE COREXP_SAVM1(JE,IGR,NGR,NLM,ITYPE,IGS,TAU)
+C
+C  This subroutine allows a recursive use of COREXP_SAVM
+C
+C                                          H.-F. Zhao : 2007
+C
+      USE DIM_MOD
+C
+      INTEGER NLM(NGR_M),ITYPE(NGR_M),IGS(NGR_M)
+      COMPLEX*16 TAU(LINMAX,LINFMAX,NATCLU_M)
+C
+      CALL COREXP_SAVM(JE,IGR,NGR,NLM,ITYPE,IGS,TAU)
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/coumat.f

@@ -0,0 +1,121 @@
+C
+C=======================================================================
+C
+      SUBROUTINE COUMAT(ITL,MI,LF,MF,DELTA,RADIAL,MATRIX)
+C
+C  This routine calculates the spin-independent PhD optical matrix
+C       elements for dipolar excitations. It is stored in
+C                          MATRIX(JDIR,JPOL)
+C
+C  Here, the conventions are :
+C
+C             IPOL=1  : linearly polarized light
+C             IPOL=2  : circularly polarized light
+C
+C             JPOL=1  : +/x polarization for circular/linear light
+C             JPOL=2  : -/y polarization for circular/linear light
+C
+C  When IDICHR=0, JDIR = 1,2 and 3 correspond respectively to the x,y
+C       and z directions for the linear polarization. But for IDICHR=1,
+C       these basis directions are those of the position of the light.
+C
+C                                     Last modified :  8 Dec 2008
+C
+      USE DIM_MOD
+C
+      USE INIT_L_MOD , L2 => NNL, L3 => LF1, L4 => LF2, L5 => ISTEP_LF
+      USE SPIN_MOD , I1 => ISPIN, N1 => NSPIN, N2 => NSPIN2, I2 => ISFLI
+     &P, I8 => IR_DIA, N3 => NSTEP
+      USE TYPCAL_MOD , I3 => IPHI, I4 => IE, I5 => ITHETA, I6 => IFTHET,
+     & I7 => IMOD, I9 => I_CP, I10 => I_EXT
+C
+      COMPLEX MATRIX(3,2),SUM_1,SUM_2,DELTA,YLM(3,-1:1),RADIAL
+      COMPLEX ONEC,IC,IL,COEF,PROD
+C
+      REAL RLM(1-NL_M:NL_M-1,1-NL_M:NL_M-1,0:NL_M-1),GNT(0:N_GAUNT)
+      REAL THETA(3),PHI(3)
+C
+      DATA PI4S3,C_LIN,SQR2 /4.188790,1.447202,1.414214/
+      DATA PIS2 /1.570796/
+C
+      ONEC=(1.,0.)
+      IC=(0.,1.)
+C
+      IF(INITL.EQ.0) GOTO 2
+C
+      M=MF-MI
+C
+      IF(MOD(LF,4).EQ.0) THEN
+        IL=ONEC
+      ELSEIF(MOD(LF,4).EQ.1) THEN
+        IL=IC
+      ELSEIF(MOD(LF,4).EQ.2) THEN
+        IL=-ONEC
+      ELSEIF(MOD(LF,4).EQ.3) THEN
+        IL=-IC
+      ENDIF
+C
+      CALL GAUNT(LI,MI,LF,MF,GNT)
+C
+      IF(ITL.EQ.0) THEN
+c        COEF=CEXP(IC*DELTA)*CONJG(IL)
+        COEF=CEXP(IC*DELTA)*IL
+      ELSE
+        IF(IDICHR.EQ.0) THEN
+c          COEF=PI4S3*CONJG(IL)
+          COEF=PI4S3*IL
+        ELSE
+c          COEF=C_LIN*CONJG(IL)
+          COEF=C_LIN*IL
+        ENDIF
+      ENDIF
+C
+      PROD=COEF*RADIAL*GNT(1)
+C
+      IF(IDICHR.EQ.0) THEN
+        YLM(1,-1)=(0.345494,0.)
+        YLM(1,0)=(0.,0.)
+        YLM(1,1)=(-0.345494,0.)
+        YLM(2,-1)=(0.,-0.345494)
+        YLM(2,0)=(0.,0.)
+        YLM(2,1)=(0.,-0.345494)
+        YLM(3,-1)=(0.,0.)
+        YLM(3,0)=(0.488602,0.)
+        YLM(3,1)=(0.,0.)
+C
+        DO JDIR=1,3
+          MATRIX(JDIR,1)=PROD*CONJG(YLM(JDIR,M))
+        ENDDO
+C
+      ELSEIF(IDICHR.GE.1) THEN
+C
+        THETA(1)=PIS2
+        PHI(1)=0.
+        THETA(2)=PIS2
+        PHI(2)=PIS2
+        THETA(3)=0.
+        PHI(3)=0.
+C
+        DO JDIR=1,3
+          CALL DJMN(THETA(JDIR),RLM,1)
+          SUM_1=RLM(-1,M,1)*PROD*CEXP((0.,-1.)*M*PHI(JDIR))
+          SUM_2=RLM(1,M,1)*PROD*CEXP((0.,-1.)*M*PHI(JDIR))
+          IF(IPOL.EQ.2) THEN
+            MATRIX(JDIR,1)=SQR2*SUM_1
+            MATRIX(JDIR,2)=SQR2*SUM_2
+          ELSEIF(ABS(IPOL).EQ.1) THEN
+            MATRIX(JDIR,1)=(SUM_2-SUM_1)
+            MATRIX(JDIR,2)=(SUM_2+SUM_1)*IC
+          ENDIF
+        ENDDO
+      ENDIF
+      GOTO 1
+C
+   2  DO JDIR=1,3
+        MATRIX(JDIR,1)=ONEC
+        MATRIX(JDIR,2)=ONEC
+      ENDDO
+C
+   1  RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/dwsph.f

@@ -0,0 +1,85 @@
+C
+C=======================================================================
+C
+      SUBROUTINE DWSPH(JTYP,JE,X,TLT,ISPEED)
+C
+C  This routine recomputes the T-matrix elements taking into account the
+C                mean square displacements.
+C
+C  When the argument X is tiny, no vibrations are taken into account
+C
+C                                      Last modified : 25 Apr 2013
+C
+      USE DIM_MOD
+C
+      USE TRANS_MOD
+C
+      DIMENSION GNT(0:N_GAUNT)
+C
+      COMPLEX TLT(0:NT_M,4,NATM,NE_M),SL1,ZEROC
+C
+      COMPLEX*16 FFL(0:2*NL_M)
+C
+      DATA PI4,EPS /12.566371,1.0E-10/
+C
+      ZEROC=(0.,0.)
+C
+      IF(X.GT.EPS) THEN
+C
+C  Standard case: vibrations
+C
+        IF(ISPEED.LT.0) THEN
+          NSUM_LB=ABS(ISPEED)
+        ENDIF
+C
+        COEF=PI4*EXP(-X)
+        NL2=2*LMAX(JTYP,JE)+2
+        IBESP=5
+        MG1=0
+        MG2=0
+C
+        CALL BESPHE(NL2,IBESP,X,FFL)
+C
+        DO L=0,LMAX(JTYP,JE)
+          XL=FLOAT(L+L+1)
+          SL1=ZEROC
+C
+          DO L1=0,LMAX(JTYP,JE)
+            XL1=FLOAT(L1+L1+1)
+            CALL GAUNT(L,MG1,L1,MG2,GNT)
+            L2MIN=ABS(L1-L)
+            IF(ISPEED.GE.0) THEN
+              L2MAX=L1+L
+            ELSEIF(ISPEED.LT.0) THEN
+              L2MAX=L2MIN+2*(NSUM_LB-1)
+            ENDIF
+            SL2=0.
+C
+            DO L2=L2MIN,L2MAX,2
+              XL2=FLOAT(L2+L2+1)
+              C=SQRT(XL1*XL2/(PI4*XL))
+              SL2=SL2+C*GNT(L2)*REAL(DREAL(FFL(L2)))
+            ENDDO
+C
+            SL1=SL1+SL2*TL(L1,1,JTYP,JE)
+          ENDDO
+C
+          TLT(L,1,JTYP,JE)=COEF*SL1
+C
+        ENDDO
+C
+      ELSE
+C
+C  Argument X tiny: no vibrations
+C
+        DO L=0,LMAX(JTYP,JE)
+C
+          TLT(L,1,JTYP,JE)=TL(L,1,JTYP,JE)
+C
+        ENDDO
+C
+      ENDIF
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/facdif.f

@@ -0,0 +1,26 @@
+C
+C=======================================================================
+C
+      SUBROUTINE FACDIF(COSTH,JAT,JE,FTHETA)
+C
+C  This routine computes the plane wave scattering factor
+C
+      USE DIM_MOD
+C
+      USE TRANS_MOD
+C
+      DIMENSION PL(0:100)
+C
+      COMPLEX FTHETA
+C
+      FTHETA=(0.,0.)
+      NL=LMAX(JAT,JE)+1
+      CALL POLLEG(NL,COSTH,PL)
+      DO 20 L=0,NL-1
+        FTHETA=FTHETA+(2*L+1)*TL(L,1,JAT,JE)*PL(L)
+  20  CONTINUE
+      FTHETA=FTHETA/VK(JE)
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/facdif1.f

@@ -0,0 +1,113 @@
+C
+C=======================================================================
+C
+      SUBROUTINE FACDIF1(VKE,RJ,RJK,THRJ,PHIRJ,BETA,GAMMA,L,M,FSPH,JAT,J
+     &E,*)
+C
+C  This routine computes a spherical wave scattering factor
+C
+C                                        Last modified : 03/04/2006
+C
+      USE DIM_MOD
+      USE APPROX_MOD
+      USE EXPFAC_MOD
+      USE TRANS_MOD
+      USE TYPCAL_MOD , I2 => IPHI, I3 => IE, I4 => ITHETA, I5 => IMOD, I
+     &6 => IPOL, I7 => I_CP, I8 => I_EXT, I9 => I_TEST
+C
+      DIMENSION PLMM(0:100,0:100)
+      DIMENSION D(1-NL_M:NL_M-1,1-NL_M:NL_M-1,0:NL_M-1)
+C
+      COMPLEX HLM(0:NO_ST_M,0:NL_M-1),HLN(0:NO_ST_M,0:NL_M-1),FSPH,RHOJ
+      COMPLEX HLM1,HLM2,HLM3,HLM4,ALMU,BLMU,SLP,SNU,SMU,VKE
+      COMPLEX RHOJK
+C
+C
+      DATA PI/3.141593/
+C
+      A=1.
+      INTER=0
+      IF(ITL.EQ.1) VKE=VK(JE)
+      RHOJ=VKE*RJ
+      RHOJK=VKE*RJK
+      HLM1=(1.,0.)
+      HLM2=(1.,0.)
+      HLM3=(1.,0.)
+      HLM4=(1.,0.)
+      IEM=1
+      CSTH=COS(BETA)
+      IF((IFTHET.EQ.0).OR.(THRJ.LT.0.0001)) THEN
+        INTER=1
+        BLMU=SQRT(4.*PI/FLOAT(2*L+1))*CEXP((0.,-1.)*M*(PHIRJ-PI))
+      ENDIF
+      CALL PLM(CSTH,PLMM,LMAX(JAT,JE))
+      IF(ISPHER.EQ.0) NO1=0
+      IF(ISPHER.EQ.1) THEN
+        IF(NO.EQ.8) THEN
+          NO1=LMAX(JAT,JE)+1
+        ELSE
+          NO1=NO
+        ENDIF
+        CALL POLHAN(ISPHER,NO1,LMAX(JAT,JE),RHOJ,HLM)
+         IF(IEM.EQ.0) THEN
+          HLM4=HLM(0,L)
+        ENDIF
+        IF(RJK.GT.0.0001) THEN
+          NDUM=0
+          CALL POLHAN(ISPHER,NDUM,LMAX(JAT,JE),RHOJK,HLN)
+        ENDIF
+        CALL DJMN(THRJ,D,L)
+        A1=ABS(D(0,M,L))
+        IF(((A1.LT.0.0001).AND.(IFTHET.EQ.1)).AND.(INTER.EQ.0)) RETURN 1
+     &
+      ENDIF
+      MUMAX=MIN0(L,NO1)
+      SMU=(0.,0.)
+      DO 10 MU=0,MUMAX
+        IF(MOD(MU,2).EQ.0) THEN
+          B=1.
+        ELSE
+          B=-1.
+          IF(SIN(BETA).LT.0.) THEN
+            A=-1.
+          ENDIF
+        ENDIF
+        IF(ISPHER.LE.1) THEN
+          ALMU=(1.,0.)
+          C=1.
+        ENDIF
+        IF(ISPHER.EQ.0) GOTO 40
+        IF(INTER.EQ.0) BLMU=CMPLX(D(M,0,L))
+        IF(MU.GT.0) THEN
+          C=B*FLOAT(L+L+1)/EXPF(MU,L)
+          ALMU=(D(M,MU,L)*CEXP((0.,-1.)*MU*GAMMA)+B*
+     *          CEXP((0.,1.)*MU*GAMMA)*D(M,-MU,L))/BLMU
+        ELSE
+          C=1.
+          ALMU=CMPLX(D(M,0,L))/BLMU
+        ENDIF
+  40    SNU=(0.,0.)
+        NU1=INT(0.5*(NO1-MU)+0.0001)
+        NUMAX=MIN0(NU1,L-MU)
+        DO 20 NU=0,NUMAX
+          SLP=(0.,0.)
+          LPMIN=MAX0(MU,NU)
+          DO 30 LP=LPMIN,LMAX(JAT,JE)
+            IF(ISPHER.EQ.1) THEN
+              HLM1=HLM(NU,LP)
+              IF(RJK.GT.0.0001) HLM3=HLN(0,LP)
+            ENDIF
+            SLP=SLP+FLOAT(2*LP+1)*TL(LP,1,JAT,JE)*HLM1*PLMM(LP,MU)*HLM3
+  30      CONTINUE
+          IF(ISPHER.EQ.1) THEN
+            HLM2=HLM(MU+NU,L)
+          ENDIF
+          SNU=SNU+SLP*HLM2
+  20    CONTINUE
+        SMU=SMU+SNU*C*ALMU*A*B
+  10  CONTINUE
+      FSPH=SMU/(VKE*HLM4)
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/gaunt_st.f

@@ -0,0 +1,126 @@
+C
+C=======================================================================
+C
+      SUBROUTINE GAUNT_ST(LMAX_T)
+C
+C   This subroutine calculates the Gaunt coefficient G(L2,L3|L1)
+C    using a downward recursion scheme due to Schulten and Gordon 
+C    for the Wigner's 3j symbols. The result is stored as GNT(L3), 
+C    making use of the selection rule M3 = M1 - M2.  
+C     
+C   Ref. :  K. Schulten and R. G. Gordon, J. Math. Phys. 16, 1961 (1975) 
+C
+C   This is the double precision version where the values are stored
+C
+C                                          Last modified : 14 May 2009
+C
+C
+      USE DIM_MOD
+      USE LOGAMAD_MOD
+      USE GAUNT_C_MOD
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+C
+      INTEGER LMAX_T
+C
+      REAL*8 F(0:N_GAUNT),G(0:N_GAUNT),A(0:N_GAUNT),A1(0:N_GAUNT)
+      REAL*8 B(0:N_GAUNT)
+C
+      DATA PI4/12.566370614359D0/
+C
+      DO L1=0,LMAX_T
+        IL1=L1*L1+L1+1
+        DO M1=-L1,L1
+          IND1=IL1+M1
+          LM1=L1+M1
+          KM1=L1-M1
+          DO L2=0,LMAX_T
+            IL2=L2*L2+L2+1
+C
+            IF(MOD(M1,2).EQ.0) THEN
+              COEF=DSQRT(DFLOAT((L1+L1+1)*(L2+L2+1))/PI4)
+            ELSE
+              COEF=-DSQRT(DFLOAT((L1+L1+1)*(L2+L2+1))/PI4)
+            ENDIF
+C
+            L12=L1+L2
+            K12=L1-L2
+            L12_1=L12+L12+1
+            L12_2=L12*L12
+            L12_21=L12*L12+L12+L12+1
+            K12_2=K12*K12
+C
+            F(L12+1)=0.D0
+            G(L12+1)=0.D0
+            A(L12+1)=0.D0
+            A1(L12+1)=0.D0
+            A1(L12)=2.D0*DSQRT(DFLOAT(L1*L2*L12_1*L12_2))
+            D1=GLD(L2+L2+1,1)-GLD(L12_1+1,1)
+            D5=0.5D0*(GLD(L1+L1+1,1)+GLD(L2+L2+1,1)-GLD(L12_1+1,1))
+            D6=GLD(L12+1,1)-GLD(L1+1,1)-GLD(L2+1,1)
+C
+            IF(MOD(K12,2).EQ.0) THEN
+              G(L12)=DEXP(D5+D6)
+            ELSE
+              G(L12)=-DEXP(D5+D6)
+            ENDIF
+C
+            DO M2=-L2,L2
+              IND2=IL2+M2
+C
+              M3=M1-M2
+              LM2=L2+M2
+              KM2=L2-M2
+C
+              DO J=1,N_GAUNT
+                GNT(J,IND2,IND1)=0.D0
+              ENDDO
+C
+              IF((ABS(M1).GT.L1).OR.(ABS(M2).GT.L2)) GOTO 10
+C
+              D2=GLD(L1+L1+1,1)-GLD(LM2+1,1)
+              D3=GLD(L12+M3+1,1)-GLD(KM2+1,1)
+              D4=GLD(L12-M3+1,1)-GLD(LM1+1,1)-GLD(KM1+1,1)
+C
+              IF(MOD(KM1-KM2,2).EQ.0) THEN
+                F(L12)=DSQRT(DEXP(D1+D2+D3+D4))
+              ELSE
+                F(L12)=-DSQRT(DEXP(D1+D2+D3+D4))
+              ENDIF
+C
+              A(L12)=2.D0*DSQRT(DFLOAT(L1*L2*L12_1*(L12_2-M3*M3)))
+              B(L12)=-DFLOAT(L12_1*((L2*L2-L1*L1-K12)*M3+L12*(L12+1)
+     1               *(M2+M1)))
+C
+              IF(ABS(M3).LE.L12) THEN 
+                GNT(L12,IND2,IND1)=COEF*F(L12)*G(L12)*
+     1                              DSQRT(DFLOAT(L12_1))
+              ENDIF
+C
+              JMIN=MAX0(ABS(K12),ABS(M3))
+C
+              DO J=L12-1,JMIN,-1
+                J1=J+1
+                J2=J+2
+                JJ=J*J
+                A1(J)=DSQRT(DFLOAT(JJ*(JJ-K12_2)*(L12_21-JJ)))
+                A(J)=DSQRT(DFLOAT((JJ-K12_2)*(L12_21-JJ)*(JJ-M3*M3)))
+                B(J)=-DFLOAT((J+J1)*(L2*(L2+1)*M3-L1*(L1+1)*M3+J*J1*
+     1                        (M2+M1)))
+                F(J)=-(DFLOAT(J1)*A(J2)*F(J2)+B(J1)*F(J1))/(DFLOAT(J2)*
+     1                 A(J1))
+                G(J)=-(DFLOAT(J1)*A1(J2)*G(J2))/(DFLOAT(J2)*A1(J1))
+C
+                IF(ABS(M3).LE.J) THEN 
+                  GNT(J,IND2,IND1)=COEF*F(J)*G(J)*DSQRT(DFLOAT(J+J1))
+                ENDIF
+              ENDDO
+C
+            ENDDO
+          ENDDO
+        ENDDO
+      ENDDO
+C
+  10  RETURN
+C
+      END
+ 

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/main.f

@@ -0,0 +1,21 @@
+      SUBROUTINE RUN(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_, 
+     &  NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_, 
+     &  NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_, 
+     &  N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_)
+
+        USE DIM_MOD
+        IMPLICIT INTEGER (A-Z)
+CF2PY   INTEGER, INTENT(IN,COPY) :: NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_
+CF2PY   INTEGER, INTENT(IN,COPY) :: NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_ 
+CF2PY   INTEGER, INTENT(IN,COPY) :: NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_
+CF2PY   INTEGER, INTENT(IN,COPY) :: N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_
+
+        CALL ALLOCATION(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_, 
+     &  NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_, 
+     &  NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_, 
+     &  N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_)
+
+        CALL MAIN_PHD_NS_CE()
+        CALL CLOSE_ALL_FILES()
+
+      END SUBROUTINE RUN

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/mpis.f

@@ -0,0 +1,280 @@
+C
+C
+C======================================================================
+C
+      SUBROUTINE MPIS(N,NLM,ITYP,IGS,JE,QI,TAU)
+C
+C
+C   This subroutine construct the correlation matrices and uses
+C    LU decomposition method to do the matrix inversion.
+C   The inverse matrix which is the contribution of a small atom group 
+C    is kept for further use.
+C
+C                                                  H. -F. Zhao : 2007
+C
+C                                    Last modified (DS) : 13 May 2009
+C
+      USE DIM_MOD
+      USE COOR_MOD
+      USE INIT_L_MOD
+      USE GAUNT_C_MOD
+      USE TRANS_MOD
+      USE CORREXP_MOD
+C
+      INTEGER NLM(NGR_M),ITYP(NGR_M),IGS(NGR_M)
+      COMPLEX*16 TAU(LINMAX,LINFMAX,NATCLU_M)
+C
+      REAL QI
+C
+      COMPLEX*16 ZEROC,ONEC,IC
+C
+      COMPLEX*16 ATTL(0:NT_M,NATM)
+      COMPLEX*16 EXPJN,ATTJN
+      COMPLEX*16 YLM(0:NLTWO,-NLTWO:NLTWO)
+      COMPLEX*16 HL1(0:NLTWO)
+      COMPLEX*16 SUM_L,SUM_L2
+      COMPLEX*16 SUM_L_A,SUM_L2_A,SUM_L_B,SUM_L2_B
+C
+      REAL*8 FOURPI
+      REAL*8 XJN,YJN,ZJN,RJN,KRJN,ZDJN
+      REAL*8 IM_VK,RE_VK
+C
+      INTEGER IPIV(NLMM),ONE_L,IN1
+C
+      COMPLEX*16 FOURPI_IC,IC_L,IC_REF,TEMP,TEMP1,TEMP2,CN1
+      COMPLEX*16 AINV(NLMM,NLMM),IN(NLMM,LINFMAX)
+C
+      DATA FOURPI /12.566370614359D0/
+C
+      ZEROC=(0.D0,0.D0)
+      ONEC=(1.D0,0.D0)
+      IC=(0.D0,1.D0)
+      IBESS=3
+      FOURPI_IC=-IC*FOURPI
+C
+      LM0=LMAX(1,JE)
+      LM0=MIN(LM0,LF2)
+      NRHS=(LM0+1)*(LM0+1)
+      INDJ=0
+C
+      NM=0
+      DO I=1,N-1
+        J=NLM(I)+1
+        NM=NM+J*J
+      ENDDO
+      L=NLM(N)
+      LNMAX=L
+      L=(L+1)*(L+1)
+      NM1=NM+1
+      NML=NM+L
+      NTYP=ITYP(N)
+C
+      DO L=0,LNMAX
+        ATTL(L,N)=DCMPLX(TL(L,1,NTYP,JE))
+      ENDDO
+      IM_VK=-DIMAG(DCMPLX(VK(JE)))
+      RE_VK=DBLE(VK(JE))
+C
+C     set up matrix blocks C((N-1)*1) and D(1*(N-1))
+C
+      I=IGS(N)
+      XN=SYM_AT(1,I)
+      YN=SYM_AT(2,I)
+      ZN=SYM_AT(3,I)
+C
+      DO J=1,N-1
+        JATL=IGS(J)
+        LJMAX=NLM(J)
+        JTYP=ITYP(J)
+        J1=J-1
+C
+        XJN=DBLE(SYM_AT(1,JATL)-XN)
+        YJN=DBLE(SYM_AT(2,JATL)-YN)
+        ZJN=DBLE(SYM_AT(3,JATL)-ZN)
+        RJN=DSQRT(XJN*XJN+YJN*YJN+ZJN*ZJN)
+        KRJN=RE_VK*RJN
+        ATTJN=FOURPI_IC*DEXP(IM_VK*RJN)
+        EXPJN=(XJN+IC*YJN)/RJN
+        ZDJN=ZJN/RJN
+        CALL SPH_HAR2(2*NL_M,ZDJN,EXPJN,YLM,LNMAX+LJMAX)
+        CALL BESPHE2(LNMAX+LJMAX+1,IBESS,KRJN,HL1)
+        DO L=0,LJMAX
+          ATTL(L,J)=ATTJN*DCMPLX(TL(L,1,JTYP,JE))
+        ENDDO
+C
+        II=NM
+        IN1=-1
+        CN1=IC
+        JJ=0
+C
+        DO LN=0,LNMAX
+          ILN=LN*LN+LN+1
+          IN1=-IN1
+          CN1=-CN1*IC
+C
+          DO MLN=-LN,LN
+            INDN=ILN+MLN
+            II=II+1
+            JJ0=J1*INDJ
+            ONE_L=-IN1
+            IC_REF=-CN1*IC
+C
+            DO LJ=0,LJMAX
+              ILJ=LJ*LJ+LJ+1
+              L_MIN=ABS(LJ-LN)
+              L_MAX=LJ+LN
+              ONE_L=-ONE_L
+              IC_REF=IC_REF*IC
+C
+C  Case MLJ equal to zero
+C 
+              JJ1=JJ0+ILJ
+              IF(LJ.GE.LN) THEN
+                IC_L=-IC_REF
+              ELSE
+                IC_L=-ONEC/IC_REF
+              ENDIF
+C
+              SUM_L=ZEROC
+              SUM_L2=ZEROC
+C
+              DO L=L_MIN,L_MAX,2
+                IC_L=-IC_L
+                IF(ABS(MLN).LE.L) THEN
+                  TEMP=IC_L*HL1(L)*GNT(L,ILJ,INDN)
+                  SUM_L=SUM_L+YLM(L,MLN)*TEMP
+                  SUM_L2=SUM_L2+DCONJG(YLM(L,MLN))*TEMP
+                ENDIF
+              ENDDO
+C
+              IF(ONE_L.EQ.-1) SUM_L2=-SUM_L2
+              A(JJ1,II)=ATTL(LJ,J)*SUM_L
+              A(II,JJ1)=ATTJN*ATTL(LN,N)*SUM_L2
+C
+C
+C  Case MLJ not equal to zero
+C 
+              DO MLJ=1,LJ
+                INDJ=ILJ+MLJ
+                INDJN=ILJ-MLJ
+                JJ1=JJ0+INDJ
+                JJ1N=JJ0+INDJN
+                MA=MLN-MLJ
+                MB=MLN+MLJ
+                IF(LJ.GE.LN) THEN
+                  IC_L=-IC_REF
+                ELSE
+                  IC_L=-ONEC/IC_REF
+                ENDIF
+C
+                SUM_L_A=ZEROC
+                SUM_L2_A=ZEROC
+                SUM_L_B=ZEROC
+                SUM_L2_B=ZEROC
+C
+                DO L=L_MIN,L_MAX,2
+                  IC_L=-IC_L
+                  IF(ABS(MA).LE.L) THEN
+                    TEMP1=IC_L*HL1(L)*GNT(L,INDJ,INDN)
+                    SUM_L_A=SUM_L_A+YLM(L,MA)*TEMP1
+                    SUM_L2_A=SUM_L2_A+DCONJG(YLM(L,MA))*TEMP1
+                  ENDIF
+                  IF(ABS(MB).LE.L) THEN
+                    TEMP2=IC_L*HL1(L)*GNT(L,INDJN,INDN)
+                    SUM_L_B=SUM_L_B+YLM(L,MB)*TEMP2
+                    SUM_L2_B=SUM_L2_B+DCONJG(YLM(L,MB))*TEMP2
+                  ENDIF
+                ENDDO
+C
+                IF(ONE_L.EQ.-1) THEN
+                   SUM_L2_A=-SUM_L2_A
+                   SUM_L2_B=-SUM_L2_B
+                ENDIF
+                A(JJ1,II)=ATTL(LJ,J)*SUM_L_A
+                A(II,JJ1)=ATTJN*ATTL(LN,N)*SUM_L2_A
+                A(JJ1N,II)=ATTL(LJ,J)*SUM_L_B
+                A(II,JJ1N)=ATTJN*ATTL(LN,N)*SUM_L2_B
+              ENDDO
+C
+C
+            ENDDO
+            JJ=JJ0+INDJ
+C
+          ENDDO
+        ENDDO
+C
+        JJ=JJ-INDN
+C
+      ENDDO
+C
+C     add B to A
+C
+      DO I=NM1,NML
+        DO J=NM1,NML
+          IF(J.EQ.I) THEN
+            A(J,I)=ONEC
+          ELSE
+            A(J,I)=ZEROC
+          ENDIF
+        ENDDO
+      ENDDO
+C
+C     construct AINV
+C
+      DO I=1,NML
+        DO J=1,NML
+          AINV(J,I)=A(J,I)
+        ENDDO
+      ENDDO
+C
+C
+C          matrix inversion(ax=b)
+C
+      CALL ZGETRF(NML,NML,AINV,NLMM,IPIV,INFO1)
+      IF(INFO1.NE.0) THEN
+        WRITE(6,*) '     --->  INFO1 =',INFO1
+      ELSE
+C
+        DO I=1,NRHS
+          DO J=1,NML
+            IF(J.EQ.I) THEN
+              IN(J,I)=(1.D0,0.D0)
+            ELSE
+              IN(J,I)=(0.D0,0.D0)
+            ENDIF
+          ENDDO
+        ENDDO
+C
+        CALL ZGETRS('N',NML,NRHS,AINV,NLMM,IPIV,IN,NLMM,INFO)
+        IF(INFO.NE.0) THEN
+          WRITE(6,*) '    --->  INFO =',INFO
+        ENDIF
+      ENDIF
+C
+C      sum of tau
+C
+      KLIN=0
+      DO K=1,N
+        KATL=IGS(K)
+        LMK=NLM(K)
+        INDKM=(LMK+1)*(LMK+1)
+C
+        DO INDJ=1,NRHS
+C
+          DO INDK=1,INDKM
+            KLIN=KLIN+1
+C
+            TAU(INDK,INDJ,KATL)=TAU(INDK,INDJ,KATL)
+     1                 +DBLE(QI)*IN(KLIN,INDJ)
+C
+          ENDDO
+          KLIN=KLIN-INDKM
+C
+        ENDDO
+        KLIN=KLIN+INDKM
+C
+      ENDDO
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/ms_cor.f

@@ -0,0 +1,165 @@
+C
+C
+C======================================================================
+C
+      SUBROUTINE MS_COR(JE,TAU)
+C
+C
+C   This subroutine calculates the scattering path operator by 
+C    the correlation expansion method.
+C
+C   The scattering path operator matrix of each small atom group 
+C    is obtained by using LU decomposition method.
+C
+C   The running time of matrix inversion subroutine used in this program
+C    scales with N^3, the memory occupied scales with N^2. We advise user to
+C    use full MS method to get the scattering path operator, i.e. directly
+C    with matrix inversion method if NGR is larger than 3. If NGR is less
+C    than 4 (i.e <=3) this subroutine will gain time.
+C
+C   This subroutine never gain memory comparing to the subrourine INV_MAT_MS
+C    as I use three large matrices stored in common, each matrix is larger or
+C    as large as the matrix used in INV_MAT_MS.
+C
+C   As I don't find a good way to solve the group problem, where all the contribution
+C    of group IGR<=NGR are collected and each small contribution has to be stored
+C    for the further larger-atom-group contribution, it's better that users change the
+C    parameter NGR_M which is set in included file 'spec.inc' to be NGR or NGR+1
+C    where NGR is the cut-off.user insterested. this subrouitne works for NGR is less
+C    than 6(<=5), if users want to calculate larger NGR, they should modify the code here
+C    to make them workable, the code is marked by 'C' in each lines (about 300 lines
+C    below here), users just release them until to the desired cut-off, the maximum is
+C    9, however, users can enlarge it if they want to. Warning ! NGR_M set in 
+C    included file should be larger than NGR and the figure listed below, don't forget
+C    to compile the code after modification.
+C
+C   Users can modify the code to make it less memory-occupied, however, no matter they
+C    do, the memories that used are more than full MS method used, so the only advantage
+C    that this code has is to gain time when NGR<=3, with command 'common' used here,
+C    the code will run faster.
+C
+C                                                      H.-F. Zhao : 2007
+C
+C                           (Photoelectron case)
+C
+C                                         Last modified : 31 Jan 2008
+C
+C
+C
+      USE DIM_MOD
+      USE COOR_MOD
+      USE INIT_L_MOD
+      USE TRANS_MOD
+      USE APPROX_MOD
+      USE CORREXP_MOD
+      USE Q_ARRAY_MOD
+C
+      COMPLEX*16 TAU1(LINMAX,LINFMAX,NATCLU_M),ONEC,ZEROC
+C
+      INTEGER NLM(NGR_M),ITYP(NGR_M),IGS(NGR_M)
+C
+      COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M),TLJ
+C
+C
+      ONEC=(1.D0,0.D0)
+      ZEROC=(0.D0,0.D0)
+C
+      LM0=LMAX(1,JE)
+      LM0=MIN(LM0,LF2)
+      NRHS=(LM0+1)*(LM0+1)
+C
+      NGR_MAX=NGR_M
+      NGR=NDIF
+C
+      IF(NGR_M.GT.NATCLU) THEN
+        WRITE(6,*) '     --->  NGR_M should be smaller than NATCLU'
+        WRITE(6,*) '     --->  it is reduced to NATCLU=',NATCLU
+        NGR_MAX=NATCLU
+      ENDIF
+C
+      IF(NGR.LT.1) THEN
+        WRITE(6,*) '     --->  NGR < 1, no expansion is done'
+        STOP
+      ELSE
+        IF(NGR.GT.NGR_MAX) THEN
+          WRITE(6,*) '     --->  NGR is too large, reduce to NGR_M=',
+     &                      NGR_MAX
+          NGR=NGR_MAX
+        ENDIF
+      ENDIF
+C
+C      Case NGR = 1
+C
+      IF(NGR.EQ.1) THEN
+        DO LJ=0,LM0
+          ILJ=LJ*LJ+LJ+1
+          TLJ=TL(LJ,1,1,JE)
+          DO MJ=-LJ,LJ
+            INDJ=ILJ+MJ
+            TAU(INDJ,INDJ,1)=TLJ
+          ENDDO
+        ENDDO
+C
+        GOTO 100
+      ENDIF
+C
+C     NGR >=2 case
+C
+C
+      DO INDJ=1,NRHS
+        TAU1(INDJ,INDJ,1)=DBLE(Q(1))*ONEC
+      ENDDO
+C
+C     Constructs the group matrix and inverses it
+C
+      IGR=1
+      LMJ=LMAX(1,JE)
+      NLM(IGR)=LMJ
+      INDJM=(LMJ+1)*(LMJ+1)
+      ITYP(IGR)=1
+      IGS(IGR)=1
+C
+      DO I=1,INDJM
+        DO J=1,INDJM
+          IF (J.EQ.I) THEN
+            A(J,I)=ONEC
+          ELSE
+            A(J,I)=ZEROC
+          ENDIF
+        ENDDO
+      ENDDO
+C
+      IGR=IGR+1
+      CALL COREXP_SAVM(JE,IGR,NGR,NLM,ITYP,IGS,TAU1)
+      IGR=IGR-1
+C
+C     TAU=TAU*tj
+C
+      DO KTYP=1,N_PROT
+        NBTYPK=NATYP(KTYP)
+        LMK=LMAX(KTYP,JE)
+        INDKM=(LMK+1)*(LMK+1)
+        DO KNUM=1,NBTYPK
+          KATL=NCORR(KNUM,KTYP)
+C
+          DO LJ=0,LM0
+            ILJ=LJ*LJ+LJ+1
+            TLJ=TL(LJ,1,1,JE)
+            DO MJ=-LJ,LJ
+              INDJ=ILJ+MJ
+C
+              DO INDK=1,INDKM
+                TAU(INDK,INDJ,KATL)=CMPLX(TAU1(INDK,INDJ,KATL))*TLJ
+              ENDDO
+C
+            ENDDO
+          ENDDO
+C
+        ENDDO
+      ENDDO
+C
+ 100  CONTINUE
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/plotfd.f

@@ -0,0 +1,106 @@
+C
+C=======================================================================
+C
+      SUBROUTINE PLOTFD(A,LMX,ITL,NL,NAT,NE)
+C
+C  This routine prepares the output for a plot of the scattering factor
+C
+      USE DIM_MOD
+C
+      USE APPROX_MOD
+      USE FDIF_MOD
+      USE INIT_L_MOD , L => LI, I2 => INITL, I3 => NNL, I4 => LF1, I5 =>
+     & LF2, I10 => ISTEP_LF
+      USE INIT_J_MOD
+      USE OUTFILES_MOD
+      USE OUTUNITS_MOD
+      USE PARCAL_MOD , N3 => NPHI, N4 => NE, N5 => NTHETA, N6 => NEPS
+      USE TYPCAL_MOD , I7 => IFTHET, I8 => IMOD, I9 => IPOL, I12 => I_CP
+     &, I13 => I_EXT, I14 => I_TEST
+      USE VALIN_MOD , U1 => THLUM, U2 => PHILUM, U3 => ELUM, N7 => NONVO
+     &L
+      USE VALFIN_MOD
+C
+C
+C
+      DIMENSION LMX(NATM,NE_M)
+C
+      COMPLEX FSPH,VKE
+C
+C
+C
+      DATA PI,CONV/3.141593,0.512314/
+C
+      OPEN(UNIT=IUO3, FILE=OUTFILE3, STATUS='UNKNOWN')
+      IF(ISPHER.EQ.0) THEN
+         L=0
+         LMAX=0
+      ELSE
+        LMAX=L
+      ENDIF
+      PHITOT=360.
+      THTOT=360.*ITHETA*(1-IPHI)+180.*ITHETA*IPHI
+      NPHI=(NFTHET+1)*IPHI+(1-IPHI)
+      NTHT=(NFTHET+1)*ITHETA*(1-IPHI)+(NFTHET/2+1)*ITHETA*IPHI+
+     *     (1-ITHETA)
+      NE=NFTHET*IE + (1-IE)
+      WRITE(IUO3,1) ISPHER,NL,NAT,L,NTHT,NPHI,NE,E0,EFIN
+      DO 10 JT=1,NTHT
+        DTHETA=THETA1+FLOAT(JT-1)*THTOT/FLOAT(MAX0(NTHT-1,1))
+        RTHETA=DTHETA*PI/180.
+        TEST=SIN(RTHETA)
+        IF(TEST.GE.0.) THEN
+          POZ=PI
+          EPS=1.
+        ELSE
+          POZ=0.
+          EPS=-1.
+        ENDIF
+        BETA=RTHETA*EPS
+        IF(ABS(TEST).LT.0.0001) THEN
+          NPHIM=1
+        ELSE
+          NPHIM=NPHI
+        ENDIF
+        DO 20 JP=1,NPHIM
+          DPHI=PHI1+FLOAT(JP-1)*PHITOT/FLOAT(MAX0(NPHI-1,1))
+          RPHI=DPHI*PI/180.
+          GAMMA=POZ-RPHI
+          DO 30 JE=1,NE
+            IF(NE.EQ.1) THEN
+              ECIN=E0
+            ELSE
+              ECIN=E0+FLOAT(JE-1)*(EFIN-E0)/FLOAT(NE-1)
+            ENDIF
+            IF(ITL.EQ.0) VKE=SQRT(ECIN-ABS(VINT))*CONV*A*(1.,0.)
+            DO 40 JAT=1,NAT
+              IF(L.GT.LMX(JAT,JE)) GOTO 90
+              DO 50 M=-LMAX,LMAX
+                CALL FACDIF1(VKE,R1,R2,THETA0,PHI0,BETA,GAMMA,L,M,FSPH,J
+     &AT,JE,*60)
+                GOTO 70
+  60            WRITE(IUO1,80)
+                STOP
+  70            REFTH=REAL(FSPH)
+                XIMFTH=AIMAG(FSPH)
+                WRITE(IUO3,5) JE,JAT,L,M,REFTH,XIMFTH,DTHETA,DPHI,ECIN
+  50          CONTINUE
+              GOTO 40
+  90          WRITE(IUO1,100) JAT
+              STOP
+  40        CONTINUE
+  30      CONTINUE
+  20    CONTINUE
+  10  CONTINUE
+      CLOSE(IUO3)
+   1  FORMAT(5X,I1,2X,I2,2X,I4,2X,I2,2X,I3,2X,I3,2X,I3,2X,F8.2,2X,F8.2)
+   5  FORMAT(1X,I3,1X,I4,1X,I2,1X,I3,1X,F6.3,1X,F6.3,1X,F6.2,1X,F6.2,1X,
+     &F8.2)
+  80  FORMAT(15X,'<<<<<  WRONG VALUE OF THETA0 : THE DENOMINATOR ','IS Z
+     &ERO  >>>>>')
+ 100  FORMAT(15X,'<<<<<  THE VALUE OF L EST IS TOO LARGE FOR ATOM',' : '
+     &,I2,'  >>>>>')
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/treat_phd.f

@@ -0,0 +1,769 @@
+C
+C=======================================================================
+C
+      SUBROUTINE TREAT_PHD(ISOM,NFICHLEC,JFICH,NP)
+C
+C  This routine sums up the calculations corresponding to different
+C        absorbers or different planes when this has to be done
+C             (parameter ISOM in the input data file).
+C
+C                                       Last modified : 24 Jan 2013
+C
+      USE DIM_MOD
+      USE OUTUNITS_MOD
+      USE TYPEXP_MOD , DUMMY => SPECTRO
+      USE VALIN_MOD
+      USE VALFIN_MOD
+C
+      PARAMETER(N_HEAD=5000,N_FILES=1000)
+C
+      CHARACTER*3 SPECTRO
+C
+      CHARACTER*13 OUTDATA
+      CHARACTER*72 HEAD(N_HEAD,N_FILES)
+C
+      REAL TAB(NDIM_M,4)
+      REAL ECIN(NE_M),DTHETA(NTH_M),DPHI(NPH_M)
+C
+C
+      DATA JVOL,JTOT/0,-1/
+C
+      REWIND IUO2
+C
+C  Reading and storing the headers:
+C
+      NHEAD=0
+      DO JLINE=1,N_HEAD
+        READ(IUO2,888) HEAD(JLINE,JFICH)
+        NHEAD=NHEAD+1
+        IF(HEAD(JLINE,JFICH)(1:6).EQ.'      ') GOTO 333
+      ENDDO
+C
+ 333  CONTINUE
+C
+      READ(IUO2,15) SPECTRO,OUTDATA
+      READ(IUO2,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE,IPH_1
+     &,I_EXT
+C
+      IF(I_EXT.EQ.2) THEN
+         IPH_1=0
+      ENDIF
+C
+      IF(ISOM.EQ.0) THEN
+C
+C........  ISOM = 0 : case of independent input files  .................
+C
+        READ(IUO2,1) NPLAN,NEMET,NTHETA,NPHI,NE
+C
+        IF(IPH_1.EQ.1) THEN
+          N_FIXED=NPHI
+          FIX0=PHI0
+          FIX1=PHI1
+          N_SCAN=NTHETA
+        ELSE
+          N_FIXED=NTHETA
+          FIX0=THETA0
+          FIX1=THETA1
+          IF(STEREO.EQ.'YES') THEN
+            NPHI=INT((PHI1-PHI0)*FLOAT(NTHETA-1)/(THETA1-THETA0)+0.0001)
+     &+1
+            IF(NTHETA*NPHI.GT.NPH_M) GOTO 37
+          ENDIF
+          N_SCAN=NPHI
+        ENDIF
+C
+        IF(I_EXT.EQ.-1) THEN
+          N_SCAN=2*N_SCAN
+        ENDIF
+C
+        IF((I_EXT.EQ.0).OR.(I_EXT.EQ.1)) THEN
+          NDP=NEMET*NTHETA*NPHI*NE
+        ELSEIF(I_EXT.EQ.-1) THEN
+          NDP=NEMET*NTHETA*NPHI*NE*2
+        ELSEIF(I_EXT.EQ.2) THEN
+          NDP=NEMET*NTHETA*NE
+          N_FIXED=NTHETA
+          N_SCAN=NPHI
+          IF((N_FIXED.GT.NTH_M).OR.(N_FIXED.GT.NPH_M)) GOTO 35
+        ENDIF
+C
+        NTT=NPLAN*NDP
+        IF(NTT.GT.NDIM_M) GOTO 5
+C
+        DO JPLAN=1,NPLAN
+         DO JEMET=1,NEMET
+          DO JE=1,NE
+C
+           DO J_FIXED=1,N_FIXED
+            IF(N_FIXED.GT.1) THEN
+              XINCRF=FLOAT(J_FIXED-1)*(FIX1-FIX0)/FLOAT(N_FIXED-1)
+            ELSEIF(N_FIXED.EQ.1) THEN
+              XINCRF=0.
+            ENDIF
+            IF(IPH_1.EQ.1) THEN
+              JPHI=J_FIXED
+            ELSE
+              THETA=THETA0+XINCRF
+              JTHETA=J_FIXED
+              IF((ABS(THETA).GT.90.).AND.(I_EXT.NE.2)) GOTO 11
+            ENDIF
+            IF(STEREO.EQ.' NO') THEN
+              N_SCAN_R=N_SCAN
+            ELSE
+              RTHETA=THETA*0.017453
+              FIX_STEP=(FIX1-FIX0)/FLOAT(N_FIXED-1)
+              N_SCAN_R=INT((PHI1-PHI0)*SIN(RTHETA)/FIX_STEP+0.0001)+1
+            ENDIF
+C
+            DO J_SCAN=1,N_SCAN_R
+             IF(IPH_1.EQ.1) THEN
+              JTHETA=J_SCAN
+             ELSE
+              JPHI=J_SCAN
+             ENDIF
+C
+             JLIN=(JPLAN-1)*NDP + (JEMET-1)*NE*N_FIXED*N_SCAN + (JE-1)*N
+     &_FIXED*N_SCAN +(JTHETA-1)*NPHI + JPHI
+C
+             IF(I_EXT.LE.0) THEN
+               IF(STEREO.EQ.' NO') THEN
+                 JPHI2=JPHI
+               ELSE
+                 JPHI2=(JTHETA-1)*NPHI+JPHI
+               ENDIF
+             ELSE
+               JPHI2=JTHETA
+             ENDIF
+C
+             READ(IUO2,2) JPL
+             IF(JPLAN.EQ.JPL) THEN
+               BACKSPACE IUO2
+               IF(IDICHR.EQ.0) THEN
+                 READ(IUO2,2) JPL,JEM,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE
+     &),TAB(JLIN,1),TAB(JLIN,2)
+                 IF(I_EXT.EQ.-1) THEN
+                   JLIN2=NTT+JLIN
+                   READ(IUO2,25) TAB(JLIN2,1),TAB(JLIN2,2)
+                 ENDIF
+               ELSE
+                 READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),DPHI(JPHI2),ECIN(J
+     &E),TAB(JLIN,1),TAB(JLIN,2),TAB(JLIN,3),TAB(JLIN,4)
+                 IF(I_EXT.EQ.-1) THEN
+                   JLIN2=NTT+JLIN
+                   READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),DPHI(JPHI2),ECIN
+     &(JE),TAB(JLIN2,1),TAB(JLIN2,2),TAB(JLIN2,3),TAB(JLIN2,4)
+                 ENDIF
+               ENDIF
+             ELSE
+               BACKSPACE IUO2
+               DO JL=JLIN,JPLAN*NDP
+                 TAB(JL,1)=0.0
+                 TAB(JL,2)=0.0
+                 TAB(JL,3)=0.0
+                 TAB(JL,4)=0.0
+               ENDDO
+               GOTO 10
+             ENDIF
+            ENDDO
+           ENDDO
+  11       CONTINUE
+          ENDDO
+         ENDDO
+  10     CONTINUE
+        ENDDO
+C
+        REWIND IUO2
+C
+C  Skipping the NHEAD lines of headers before rewriting:
+C
+        DO JLINE=1,NHEAD
+          READ(IUO2,888) HEAD(JLINE,JFICH)
+        ENDDO
+C
+        WRITE(IUO2,15) SPECTRO,OUTDATA
+        WRITE(IUO2,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE
+        WRITE(IUO2,8) NPHI,NTHETA,NE,NPLAN,ISOM
+C
+        DO JE=1,NE
+         DO JTHETA=1,NTHETA
+          IF(STEREO.EQ.' NO') THEN
+            NPHI_R=NPHI
+          ELSE
+            RTHETA=DTHETA(JTHETA)*0.017453
+            FIX_STEP=(THETA1-THETA0)/FLOAT(NTHETA-1)
+            NPHI_R=INT((PHI1-PHI0)*SIN(RTHETA)/FIX_STEP+0.0001)+1
+            NPHI=INT((PHI1-PHI0)/FIX_STEP+0.0001)+1
+          ENDIF
+          DO JPHI=1,NPHI_R
+           TOTDIF_1=0.
+           TOTDIR_1=0.
+           VOLDIF_1=0.
+           VOLDIR_1=0.
+           TOTDIF_2=0.
+           TOTDIR_2=0.
+           VOLDIF_2=0.
+           VOLDIR_2=0.
+           IF(I_EXT.EQ.-1) THEN
+             TOTDIF2_1=0.
+             TOTDIR2_1=0.
+             VOLDIF2_1=0.
+             VOLDIR2_1=0.
+             TOTDIF2_2=0.
+             TOTDIR2_2=0.
+             VOLDIF2_2=0.
+             VOLDIR2_2=0.
+           ENDIF
+C
+           DO JPLAN=1,NPLAN
+C
+            SF_1=0.
+            SR_1=0.
+            SF_2=0.
+            SR_2=0.
+            IF(I_EXT.EQ.-1) THEN
+              SF2_1=0.
+              SR2_1=0.
+              SF2_2=0.
+              SR2_2=0.
+            ENDIF
+C
+            DO JEMET=1,NEMET
+             JLIN=(JPLAN-1)*NDP + (JEMET-1)*NE*NTHETA*NPHI + (JE-1)*NTHE
+     &TA*NPHI +(JTHETA-1)*NPHI + JPHI
+             SF_1=SF_1+TAB(JLIN,2)
+             SR_1=SR_1+TAB(JLIN,1)
+             IF(I_EXT.EQ.-1) THEN
+               JLIN2=NTT+JLIN
+               SF2_1=SF2_1+TAB(JLIN2,2)
+               SR2_1=SR2_1+TAB(JLIN2,1)
+             ENDIF
+             IF(IDICHR.GE.1) THEN
+               SF_2=SF_2+TAB(JLIN,4)
+               SR_2=SR_2+TAB(JLIN,3)
+               IF(I_EXT.EQ.-1) THEN
+                 JLIN2=NTT+JLIN
+                 SF2_2=SF2_2+TAB(JLIN2,4)
+                 SR2_2=SR2_2+TAB(JLIN2,3)
+               ENDIF
+             ENDIF
+            ENDDO
+            IF(I_EXT.LE.0) THEN
+               IF(STEREO.EQ.' NO') THEN
+                 JPHI2=JPHI
+               ELSE
+                 JPHI2=(JTHETA-1)*NPHI+JPHI
+               ENDIF
+             ELSE
+              JPHI2=JTHETA
+            ENDIF
+            IF(IDICHR.EQ.0) THEN
+              WRITE(IUO2,3) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),SR
+     &_1,SF_1
+              IF(I_EXT.EQ.-1) THEN
+                WRITE(IUO2,3) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),
+     &SR2_1,SF2_1
+              ENDIF
+            ELSE
+              WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),S
+     &R_1,SF_1,SR_2,SF_2
+              IF(I_EXT.EQ.-1) THEN
+                WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE)
+     &,SR2_1,SF2_1,SR2_2,SF2_2
+              ENDIF
+            ENDIF
+            IF(JPLAN.GT.NONVOL(JFICH)) THEN
+              VOLDIF_1=VOLDIF_1+SF_1
+              VOLDIR_1=VOLDIR_1+SR_1
+              IF(I_EXT.EQ.-1) THEN
+                VOLDIF2_1=VOLDIF2_1+SF2_1
+                VOLDIR2_1=VOLDIR2_1+SR2_1
+              ENDIF
+              IF(IDICHR.GE.1) THEN
+                VOLDIF_2=VOLDIF_2+SF_2
+                VOLDIR_2=VOLDIR_1+SR_2
+                IF(I_EXT.EQ.-1) THEN
+                  VOLDIF2_2=VOLDIF2_2+SF2_2
+                  VOLDIR2_2=VOLDIR2_1+SR2_2
+                ENDIF
+              ENDIF
+            ENDIF
+            TOTDIF_1=TOTDIF_1+SF_1
+            TOTDIR_1=TOTDIR_1+SR_1
+            IF(I_EXT.EQ.-1) THEN
+              TOTDIF2_1=TOTDIF2_1+SF2_1
+              TOTDIR2_1=TOTDIR2_1+SR2_1
+            ENDIF
+            IF(IDICHR.GE.1) THEN
+              TOTDIF_2=TOTDIF_2+SF_2
+              TOTDIR_2=TOTDIR_2+SR_2
+              IF(I_EXT.EQ.-1) THEN
+                TOTDIF2_2=TOTDIF2_2+SF2_2
+                TOTDIR2_2=TOTDIR2_2+SR2_2
+              ENDIF
+            ENDIF
+           ENDDO
+           IF(IDICHR.EQ.0) THEN
+             WRITE(IUO2,3) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),VOLD
+     &IR_1,VOLDIF_1
+             IF(I_EXT.EQ.-1) THEN
+               WRITE(IUO2,3) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),VO
+     &LDIR2_1,VOLDIF2_1
+             ENDIF
+             WRITE(IUO2,3) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),TOTD
+     &IR_1,TOTDIF_1
+             IF(I_EXT.EQ.-1) THEN
+               WRITE(IUO2,3) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),TO
+     &TDIR2_1,TOTDIF2_1
+             ENDIF
+           ELSE
+             WRITE(IUO2,23) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),VOL
+     &DIR_1,VOLDIF_1,VOLDIR_2,VOLDIF_2
+            IF(I_EXT.EQ.-1) THEN
+               WRITE(IUO2,23) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),V
+     &OLDIR2_1,VOLDIF2_1,VOLDIR2_2,VOLDIF2_2
+             ENDIF
+             WRITE(IUO2,23) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),TOT
+     &DIR_1,TOTDIF_1,TOTDIR_2,TOTDIF_2
+             IF(I_EXT.EQ.-1) THEN
+               WRITE(IUO2,23) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),T
+     &OTDIR2_1,TOTDIF2_1,TOTDIR2_2,TOTDIF2_2
+             ENDIF
+           ENDIF
+          ENDDO
+         ENDDO
+        ENDDO
+C
+      ELSE
+C
+C........ ISOM not= 0 : multiple input files to be summed up  ..........
+C
+        READ(IUO2,7) NTHETA,NPHI,NE
+C
+        IF(IPH_1.EQ.1) THEN
+          N_FIXED=NPHI
+          FIX0=PHI0
+          FIX1=PHI1
+          N_SCAN=NTHETA
+        ELSE
+          N_FIXED=NTHETA
+          FIX0=THETA0
+          FIX1=THETA1
+          IF(STEREO.EQ.'YES') THEN
+            NPHI=INT((PHI1-PHI0)*FLOAT(NTHETA-1)/(THETA1-THETA0)+0.0001)
+     &+1
+            IF(NTHETA*NPHI.GT.NPH_M) GOTO 37
+          ENDIF
+          N_SCAN=NPHI
+        ENDIF
+C
+        IF(I_EXT.EQ.-1) THEN
+          N_SCAN=2*N_SCAN
+        ENDIF
+C
+        IF((I_EXT.EQ.0).OR.(I_EXT.EQ.1)) THEN
+          NDP=NTHETA*NPHI*NE
+        ELSEIF(I_EXT.EQ.-1) THEN
+          NDP=NTHETA*NPHI*NE*2
+        ELSEIF(I_EXT.EQ.2) THEN
+          NDP=NTHETA*NE
+          N_FIXED=NTHETA
+          N_SCAN=NPHI
+          IF((N_FIXED.GT.NTH_M).OR.(N_FIXED.GT.NPH_M)) GOTO 35
+        ENDIF
+C
+        NTT=NFICHLEC*NDP
+        IF(NTT.GT.NDIM_M) GOTO 5
+C
+        IF(ISOM.EQ.1) THEN
+          NPLAN=NP
+          NF=NP
+        ELSEIF(ISOM.EQ.2) THEN
+          NEMET=NFICHLEC
+          NF=NFICHLEC
+          NPLAN=1
+        ENDIF
+C
+        DO JF=1,NF
+C
+C  Reading the headers for each file:
+C
+          IF(JF.GT.1) THEN
+            DO JLINE=1,NHEAD
+              READ(IUO2,888) HEAD(JLINE,JF)
+            ENDDO
+          ENDIF
+C
+          DO JE=1,NE
+C
+           DO J_FIXED=1,N_FIXED
+            IF(N_FIXED.GT.1) THEN
+              XINCRF=FLOAT(J_FIXED-1)*(FIX1-FIX0)/FLOAT(N_FIXED-1)
+            ELSEIF(N_FIXED.EQ.1) THEN
+              XINCRF=0.
+            ENDIF
+            IF(IPH_1.EQ.1) THEN
+              JPHI=J_FIXED
+            ELSE
+              THETA=THETA0+XINCRF
+              JTHETA=J_FIXED
+              IF((ABS(THETA).GT.90.).AND.(I_EXT.NE.2)) GOTO 12
+            ENDIF
+            IF(STEREO.EQ.' NO') THEN
+              N_SCAN_R=N_SCAN
+            ELSE
+              RTHETA=THETA*0.017453
+              FIX_STEP=(FIX1-FIX0)/FLOAT(N_FIXED-1)
+              N_SCAN_R=INT((PHI1-PHI0)*SIN(RTHETA)/FIX_STEP+0.0001)+1
+            ENDIF
+C
+            DO J_SCAN=1,N_SCAN_R
+             IF(IPH_1.EQ.1) THEN
+              JTHETA=J_SCAN
+             ELSE
+              JPHI=J_SCAN
+             ENDIF
+C
+             JLIN=(JF-1)*NDP + (JE-1)*N_FIXED*N_SCAN +(JTHETA-1)*NPHI + 
+     &JPHI
+             IF(I_EXT.LE.0) THEN
+               IF(STEREO.EQ.' NO') THEN
+                 JPHI2=JPHI
+               ELSE
+                 JPHI2=(JTHETA-1)*NPHI+JPHI
+               ENDIF
+             ELSE
+               JPHI2=JTHETA
+             ENDIF
+C
+             IF(ISOM.EQ.1) THEN
+               READ(IUO2,2) JPL
+               IF(JF.EQ.JPL) THEN
+                 BACKSPACE IUO2
+                 IF(IDICHR.EQ.0) THEN
+                   READ(IUO2,2) JPL,JEM,DTHETA(JTHETA),DPHI(JPHI2),ECIN(
+     &JE),TAB(JLIN,1),TAB(JLIN,2)
+                   IF(I_EXT.EQ.-1) THEN
+                     JLIN2=NTT+JLIN
+                     READ(IUO2,25) TAB(JLIN2,1),TAB(JLIN2,2)
+                   ENDIF
+                 ELSE
+                   READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),DPHI(JPHI2),ECIN
+     &(JE),TAB(JLIN,1),TAB(JLIN,2),TAB(JLIN,3),TAB(JLIN,4)
+                   IF(I_EXT.EQ.-1) THEN
+                     JLIN2=NTT+JLIN
+                     READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),DPHI(JPHI2),EC
+     &IN(JE),TAB(JLIN2,1),TAB(JLIN2,2),TAB(JLIN2,3),TAB(JLIN2,4)
+                   ENDIF
+                 ENDIF
+               ELSE
+                 BACKSPACE IUO2
+                 DO JLINE=1,NHEAD
+                   BACKSPACE IUO2
+                 ENDDO
+                 DO JL=JLIN,JF*NDP
+                   TAB(JL,1)=0.0
+                   TAB(JL,2)=0.0
+                   TAB(JL,3)=0.0
+                   TAB(JL,4)=0.0
+                 ENDDO
+                 GOTO 13
+               ENDIF
+             ELSEIF(ISOM.EQ.2) THEN
+               IF(IDICHR.EQ.0) THEN
+                 READ(IUO2,2) JPL,JEM,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE
+     &),TAB(JLIN,1),TAB(JLIN,2)
+                 IF(I_EXT.EQ.-1) THEN
+                   JLIN2=NTT+JLIN
+                   READ(IUO2,25) TAB(JLIN2,1),TAB(JLIN2,2)
+                 ENDIF
+               ELSE
+                 READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),DPHI(JPHI2),ECIN(J
+     &E),TAB(JLIN,1),TAB(JLIN,2),TAB(JLIN,3),TAB(JLIN,4)
+                 IF(I_EXT.EQ.-1) THEN
+                   JLIN2=NTT+JLIN
+                   READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),DPHI(JPHI2),ECIN
+     &(JE),TAB(JLIN2,1),TAB(JLIN2,2),TAB(JLIN2,3),TAB(JLIN2,4)
+                   ENDIF
+               ENDIF
+             ENDIF
+            ENDDO
+ 12         CONTINUE
+           ENDDO
+          ENDDO
+ 13       CONTINUE
+        ENDDO
+C
+        REWIND IUO2
+C
+C  Writing the headers:
+C
+        DO JLINE=1,2
+          WRITE(IUO2,888) HEAD(JLINE,1)
+        ENDDO
+        DO JF=1,NFICHLEC
+          DO JLINE=3,6
+            WRITE(IUO2,888) HEAD(JLINE,JF)
+          ENDDO
+          WRITE(IUO2,888) HEAD(2,JF)
+        ENDDO
+        DO JLINE=7,NHEAD
+          WRITE(IUO2,888) HEAD(JLINE,1)
+        ENDDO
+C
+        WRITE(IUO2,15) SPECTRO,OUTDATA
+        WRITE(IUO2,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE
+        WRITE(IUO2,8) NPHI,NTHETA,NE,NPLAN,ISOM
+C
+        IF(ISOM.EQ.1) THEN
+C
+          DO JE=1,NE
+C
+           DO JTHETA=1,NTHETA
+            IF(STEREO.EQ.' NO') THEN
+              NPHI_R=NPHI
+            ELSE
+              RTHETA=DTHETA(JTHETA)*0.017453
+              FIX_STEP=(THETA1-THETA0)/FLOAT(NTHETA-1)
+              NPHI_R=INT((PHI1-PHI0)*SIN(RTHETA)/FIX_STEP+0.0001)+1
+              NPHI=INT((PHI1-PHI0)/FIX_STEP+0.0001)+1
+            ENDIF
+            DO JPHI=1,NPHI_R
+C
+             TOTDIF_1=0.
+             TOTDIR_1=0.
+             VOLDIF_1=0.
+             VOLDIR_1=0.
+             TOTDIF_2=0.
+             TOTDIR_2=0.
+             VOLDIF_2=0.
+             VOLDIR_2=0.
+             IF(I_EXT.EQ.-1) THEN
+               TOTDIF2_1=0.
+               TOTDIR2_1=0.
+               VOLDIF2_1=0.
+               VOLDIR2_1=0.
+               TOTDIF2_2=0.
+               TOTDIR2_2=0.
+               VOLDIF2_2=0.
+               VOLDIR2_2=0.
+             ENDIF
+C
+             DO JPLAN=1,NPLAN
+              JF=JPLAN
+C
+              JLIN=(JF-1)*NDP + (JE-1)*NTHETA*NPHI +(JTHETA-1)*NPHI + JP
+     &HI
+C
+              SR_1=TAB(JLIN,1)
+              SF_1=TAB(JLIN,2)
+              IF(I_EXT.EQ.-1) THEN
+                 JLIN2=NTT+JLIN
+                 SF2_1=TAB(JLIN2,2)
+                 SR2_1=TAB(JLIN2,1)
+              ENDIF
+              IF(I_EXT.LE.0) THEN
+                IF(STEREO.EQ.' NO') THEN
+                  JPHI2=JPHI
+                ELSE
+                  JPHI2=(JTHETA-1)*NPHI+JPHI
+                ENDIF
+              ELSE
+                JPHI2=JTHETA
+              ENDIF
+              IF(IDICHR.EQ.0) THEN
+                WRITE(IUO2,3) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),
+     &SR_1,SF_1
+                IF(I_EXT.EQ.-1) THEN
+                  WRITE(IUO2,3) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE
+     &),SR2_1,SF2_1
+                ENDIF
+              ELSE
+                SR_2=TAB(JLIN,3)
+                SF_2=TAB(JLIN,4)
+                IF(I_EXT.EQ.-1) THEN
+                  JLIN2=NTT+JLIN
+                  SF2_2=TAB(JLIN2,4)
+                  SR2_2=TAB(JLIN2,3)
+                ENDIF
+                WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE)
+     &,SR_1,SF_1,SR_2,SF_2
+                IF(I_EXT.EQ.-1) THEN
+                  WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(J
+     &E),SR2_1,SF2_1,SR2_2,SF2_2
+                ENDIF
+              ENDIF
+              IF(NONVOL(JPLAN).EQ.0) THEN
+                VOLDIF_1=VOLDIF_1+SF_1
+                VOLDIR_1=VOLDIR_1+SR_1
+                IF(I_EXT.EQ.-1) THEN
+                  VOLDIF2_1=VOLDIF2_1+SF2_1
+                  VOLDIR2_1=VOLDIR2_1+SR2_1
+                ENDIF
+                IF(IDICHR.GE.1) THEN
+                  VOLDIF_2=VOLDIF_2+SF_2
+                  VOLDIR_2=VOLDIR_2+SR_2
+                  IF(I_EXT.EQ.-1) THEN
+                    VOLDIF2_2=VOLDIF2_2+SF2_2
+                    VOLDIR2_2=VOLDIR2_1+SR2_2
+                  ENDIF
+                ENDIF
+              ENDIF
+              TOTDIF_1=TOTDIF_1+SF_1
+              TOTDIR_1=TOTDIR_1+SR_1
+              IF(I_EXT.EQ.-1) THEN
+                TOTDIF2_1=TOTDIF2_1+SF2_1
+                TOTDIR2_1=TOTDIR2_1+SR2_1
+              ENDIF
+              IF(IDICHR.GE.1) THEN
+                TOTDIF_2=TOTDIF_2+SF_2
+                TOTDIR_2=TOTDIR_2+SR_2
+                IF(I_EXT.EQ.-1) THEN
+                  TOTDIF2_2=TOTDIF2_2+SF2_2
+                  TOTDIR2_2=TOTDIR2_2+SR2_2
+                ENDIF
+              ENDIF
+             ENDDO
+C
+             IF(IDICHR.EQ.0) THEN
+               WRITE(IUO2,3) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),VO
+     &LDIR_1,VOLDIF_1
+               IF(I_EXT.EQ.-1) THEN
+                 WRITE(IUO2,3) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),
+     &VOLDIR2_1,VOLDIF2_1
+               ENDIF
+               WRITE(IUO2,3) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),TO
+     &TDIR_1,TOTDIF_1
+               IF(I_EXT.EQ.-1) THEN
+                 WRITE(IUO2,3) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),
+     &TOTDIR2_1,TOTDIF2_1
+               ENDIF
+             ELSE
+               WRITE(IUO2,23) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),V
+     &OLDIR_1,VOLDIF_1,VOLDIR_2,VOLDIF_2
+               IF(I_EXT.EQ.-1) THEN
+                 WRITE(IUO2,23) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE)
+     &,VOLDIR2_1,VOLDIF2_1,VOLDIR2_2,VOLDIF2_2
+               ENDIF
+               WRITE(IUO2,23) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),T
+     &OTDIR_1,TOTDIF_1,TOTDIR_2,TOTDIF_2
+               IF(I_EXT.EQ.-1) THEN
+                 WRITE(IUO2,23) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE)
+     &,TOTDIR2_1,TOTDIF2_1,TOTDIR2_2,TOTDIF2_2
+               ENDIF
+             ENDIF
+C
+            ENDDO
+           ENDDO
+          ENDDO
+        ELSEIF(ISOM.EQ.2) THEN
+          DO JE=1,NE
+C
+           DO JTHETA=1,NTHETA
+            IF(STEREO.EQ.' NO') THEN
+              NPHI_R=NPHI
+            ELSE
+              RTHETA=DTHETA(JTHETA)*0.017453
+              FIX_STEP=(THETA1-THETA0)/FLOAT(NTHETA-1)
+              NPHI_R=INT((PHI1-PHI0)*SIN(RTHETA)/FIX_STEP+0.0001)+1
+              NPHI=INT((PHI1-PHI0)/FIX_STEP+0.0001)+1
+            ENDIF
+            DO JPHI=1,NPHI_R
+C
+              SF_1=0.
+              SR_1=0.
+              SF_2=0.
+              SR_2=0.
+              IF(I_EXT.EQ.-1) THEN
+                SF2_1=0.
+                SR2_1=0.
+                SF2_2=0.
+                SR2_2=0.
+              ENDIF
+C
+              DO JEMET=1,NEMET
+               JF=JEMET
+C
+               JLIN=(JF-1)*NDP + (JE-1)*NTHETA*NPHI +(JTHETA-1)*NPHI + J
+     &PHI
+C
+               SF_1=SF_1+TAB(JLIN,2)
+               SR_1=SR_1+TAB(JLIN,1)
+               IF(I_EXT.EQ.-1) THEN
+                 JLIN2=NTT+JLIN
+                 SF2_1=SF2_1+TAB(JLIN2,2)
+                 SR2_1=SR2_1+TAB(JLIN2,1)
+               ENDIF
+               IF(IDICHR.GE.1) THEN
+                 SF_2=SF_2+TAB(JLIN,4)
+                 SR_2=SR_2+TAB(JLIN,3)
+                 IF(I_EXT.EQ.-1) THEN
+                   JLIN2=NTT+JLIN
+                   SF2_2=SF2_2+TAB(JLIN2,4)
+                   SR2_2=SR2_2+TAB(JLIN2,3)
+                 ENDIF
+               ENDIF
+              ENDDO
+              IF(I_EXT.LE.0) THEN
+                IF(STEREO.EQ.' NO') THEN
+                  JPHI2=JPHI
+                ELSE
+                  JPHI2=(JTHETA-1)*NPHI+JPHI
+                ENDIF
+              ELSE
+                JPHI2=JTHETA
+              ENDIF
+              IF(IDICHR.EQ.0) THEN
+                WRITE(IUO2,3) JPL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),SR
+     &_1,SF_1
+                IF(I_EXT.EQ.-1) THEN
+                  WRITE(IUO2,3) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE
+     &),SR2_1,SF2_1
+                ENDIF
+              ELSE
+                WRITE(IUO2,23) JPL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),S
+     &R_1,SF_1,SR_2,SF_2
+                IF(I_EXT.EQ.-1) THEN
+                  WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(J
+     &E),SR2_1,SF2_1,SR2_2,SF2_2
+                ENDIF
+              ENDIF
+            ENDDO
+           ENDDO
+          ENDDO
+        ENDIF
+      ENDIF
+C
+      GOTO 6
+C
+   5  WRITE(IUO1,4)
+      STOP
+  35  WRITE(IUO1,36) N_FIXED
+      STOP
+  37  WRITE(IUO1,38) NTHETA*NPHI
+      STOP
+C
+   1  FORMAT(2X,I3,2X,I2,2X,I4,2X,I4,2X,I4)
+   2  FORMAT(2X,I3,2X,I2,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6)
+   3  FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6)
+   4  FORMAT(//,8X,'<<<<<<<<<<  DIMENSION OF THE ARRAYS TOO SMALL ','IN 
+     &THE TREAT_PHD SUBROUTINE - INCREASE NDIM_M ','>>>>>>>>>>')
+   7  FORMAT(I4,2X,I4,2X,I4)
+   8  FORMAT(I4,2X,I4,2X,I4,2X,I3,2X,I1)
+   9  FORMAT(9(2X,I1),2X,I2)
+  15  FORMAT(2X,A3,11X,A13)
+  22  FORMAT(2X,I3,2X,I2,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6,2X,E1
+     &2.6,2X,E12.6)
+  23  FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6,2X,E12.6,2X
+     &,E12.6)
+  25  FORMAT(37X,E12.6,2X,E12.6)
+  36  FORMAT(//,4X,'<<<<<<<<<<  DIMENSION OF NTH_M OR NPH_M TOO SMALL ',
+     &'IN THE INCLUDE FILE >>>>>>>>>>',/,4X,'<<<<<<<<<<                 
+     &SHOULD BE AT LEAST ',I6,'                 >>>>>>>>>>')
+  38  FORMAT(//,8X,'<<<<<<<<<<  DIMENSION OF NPH_M TOO SMALL ','IN THE I
+     &NCLUDE FILE >>>>>>>>>>',/,8X,'<<<<<<<<<<             SHOULD BE AT 
+     &LEAST ',I6,'             >>>>>>>>>>')
+ 888  FORMAT(A72)
+C
+   6  RETURN
+C
+      END

src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/weight_sum.f

@@ -0,0 +1,335 @@
+C
+C=======================================================================
+C
+      SUBROUTINE WEIGHT_SUM(ISOM,I_EXT,I_EXT_A,JEL)
+C
+C   This subroutine performs a weighted sum of the results
+C     corresponding to different directions of the detector.
+C     The directions and weights are read from an external input file
+C
+C   JEL is the electron undetected (i.e. for which the outgoing
+C     directions are integrated over the unit sphere). It is always
+C     1 for one electron spectroscopies (PHD). For APECS, It can be
+C     1 (photoelectron) or 2 (Auger electron) or even 0 (no electron
+C     detected)
+C
+C                                         Last modified : 31 Jan 2007
+C
+      USE DIM_MOD
+      USE INFILES_MOD
+      USE INUNITS_MOD
+      USE OUTUNITS_MOD
+C
+C
+      PARAMETER(N_MAX=5810,NPM=20)
+C
+      REAL*4 W(N_MAX),W_A(N_MAX),ECIN(NE_M)
+      REAL*4 DTHETA(N_MAX),DPHI(N_MAX),DTHETAA(N_MAX),DPHIA(N_MAX)
+      REAL*4 SR_1,SF_1,SR_2,SF_2
+      REAL*4 SUMR_1(NPM,NE_M,N_MAX),SUMR_2(NPM,NE_M,N_MAX)
+      REAL*4 SUMF_1(NPM,NE_M,N_MAX),SUMF_2(NPM,NE_M,N_MAX)
+C
+      CHARACTER*3 SPECTRO,SPECTRO2
+      CHARACTER*5 LIKE
+      CHARACTER*13 OUTDATA
+C
+C
+C
+C
+      DATA JVOL,JTOT/0,-1/
+      DATA LIKE /'-like'/
+C
+      REWIND IUO2
+C
+      READ(IUO2,15) SPECTRO,OUTDATA
+      IF(SPECTRO.NE.'APC') THEN
+        READ(IUO2,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE
+        READ(IUO2,8) NPHI,NTHETA,NE,NPLAN,ISOM
+        SPECTRO2='XAS'
+      ELSE
+        READ(IUO2,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE
+        READ(IUO2,9) ISPIN_A,IDICHR_A,I_SO_A,ISFLIP_A,ICHKDIR_A,IPHI_A,I
+     &THETA_A,IE_A
+        READ(IUO2,8) NPHI,NTHETA,NE,NPLAN,ISOM
+        READ(IUO2,8) NPHI_A,NTHETA_A
+        IF(JEL.EQ.1) THEN
+          SPECTRO2='AED'
+        ELSEIF(JEL.EQ.2) THEN
+          SPECTRO2='PHD'
+        ELSEIF(JEL.EQ.0) THEN
+          SPECTRO2='XAS'
+        ENDIF
+      ENDIF
+C
+      IF(NPLAN.GT.NPM) THEN
+        WRITE(IUO1,4) NPLAN+2
+        STOP
+      ENDIF
+C
+C  Reading the number of angular points
+C
+      IF(SPECTRO.NE.'APC') THEN
+        OPEN(UNIT=IUI6, FILE=INFILE6, STATUS='OLD')
+        READ(IUI6,1) N_POINTS
+        READ(IUI6,5) I_DIM,N_DUM1,N_DUM2
+        N_POINTS_A=1
+      ELSE
+        IF(JEL.EQ.1) THEN
+          OPEN(UNIT=IUI6, FILE=INFILE6, STATUS='OLD')
+          READ(IUI6,1) N_POINTS
+          READ(IUI6,5) I_DIM,N_DUM1,N_DUM2
+          IF(I_EXT_A.EQ.0) THEN
+            N_POINTS_A=NTHETA_A*NPHI_A
+          ELSE
+            OPEN(UNIT=IUI9, FILE=INFILE9, STATUS='OLD')
+            READ(IUI9,1) N_POINTS_A
+            READ(IUI9,5) I_DIM,N_DUM1,N_DUM2
+          ENDIF
+          NTHETA0=NTHETA_A
+          NPHI0=NPHI_A
+        ELSEIF(JEL.EQ.2) THEN
+          OPEN(UNIT=IUI9, FILE=INFILE9, STATUS='OLD')
+          READ(IUI9,1) N_POINTS_A
+          READ(IUI9,5) I_DIM,N_DUM1,N_DUM2
+          IF(I_EXT.EQ.0) THEN
+            N_POINTS=NTHETA*NPHI
+          ELSE
+            OPEN(UNIT=IUI6, FILE=INFILE6, STATUS='OLD')
+            READ(IUI6,1) N_POINTS
+            READ(IUI6,5) I_DIM,N_DUM1,N_DUM2
+          ENDIF
+          NTHETA0=NTHETA
+          NPHI0=NPHI
+        ELSEIF(JEL.EQ.0) THEN
+          OPEN(UNIT=IUI6, FILE=INFILE6, STATUS='OLD')
+          OPEN(UNIT=IUI9, FILE=INFILE9, STATUS='OLD')
+          READ(IUI6,1) N_POINTS
+          READ(IUI9,1) N_POINTS_A
+          READ(IUI6,5) I_DIM,N_DUM1,N_DUM2
+          READ(IUI9,5) I_DIM,N_DUM1,N_DUM2
+        ENDIF
+      ENDIF
+C
+      IF(SPECTRO.NE.'APC') THEN
+        NANGLE=1
+      ELSE
+        IF(JEL.EQ.1) THEN
+          NANGLE=N_POINTS_A
+        ELSEIF(JEL.EQ.2) THEN
+          NANGLE=N_POINTS
+        ELSEIF(JEL.EQ.0) THEN
+          NANGLE=1
+        ENDIF
+      ENDIF
+C
+C  Initialization of the arrays
+C
+      DO JE=1,NE
+        DO JANGLE=1,NANGLE
+          DO JPLAN=1,NPLAN+2
+            SUMR_1(JPLAN,JE,JANGLE)=0.
+            SUMF_1(JPLAN,JE,JANGLE)=0.
+            IF(IDICHR.GT.0) THEN
+              SUMR_2(JPLAN,JE,JANGLE)=0.
+              SUMF_2(JPLAN,JE,JANGLE)=0.
+            ENDIF
+          ENDDO
+        ENDDO
+      ENDDO
+C
+C  Reading of the data to be angle integrated
+C
+      DO JE=1,NE
+C
+        DO JANGLE=1,N_POINTS
+          IF(I_EXT.NE.0) READ(IUI6,2) TH,PH,W(JANGLE)
+          DO JANGLE_A=1,N_POINTS_A
+            IF((I_EXT_A.NE.0).AND.(JANGLE.EQ.1)) THEN
+               READ(IUI9,2) THA,PHA,W_A(JANGLE_A)
+            ENDIF
+C
+            DO JPLAN=1,NPLAN+2
+C
+              IF(IDICHR.EQ.0) THEN
+                IF(SPECTRO.NE.'APC') THEN
+                  READ(IUO2,3) JDUM,DTHETA(JANGLE),DPHI(JANGLE),ECIN(JE)
+     &,SR_1,SF_1
+                ELSE
+                  READ(IUO2,13) JDUM,DTHETA(JANGLE),DPHI(JANGLE),ECIN(JE
+     &),DTHETAA(JANGLE_A),DPHIA(JANGLE_A),SR_1,SF_1
+                ENDIF
+              ELSE
+                IF(SPECTRO.NE.'APC') THEN
+                  READ(IUO2,23) JDUM,DTHETA(JANGLE),DPHI(JANGLE),ECIN(JE
+     &),SR_1,SF_1,SR_2,SF_2
+                ELSE
+                  READ(IUO2,24) JDUM,DTHETA(JANGLE),DPHI(JANGLE),ECIN(JE
+     &),DTHETAA(JANGLE_A),DPHIA(JANGLE_A),SR_1,SF_1,SR_2,SF_2
+                ENDIF
+              ENDIF
+C
+              IF(JEL.EQ.1) THEN
+                SUMR_1(JPLAN,JE,JANGLE_A)=SUMR_1(JPLAN,JE,JANGLE_A)+SR_1
+     &*W(JANGLE)
+                SUMF_1(JPLAN,JE,JANGLE_A)=SUMF_1(JPLAN,JE,JANGLE_A)+SF_1
+     &*W(JANGLE)
+              ELSEIF(JEL.EQ.2) THEN
+                SUMR_1(JPLAN,JE,JANGLE)=SUMR_1(JPLAN,JE,JANGLE)+SR_1*W_A
+     &(JANGLE_A)
+                SUMF_1(JPLAN,JE,JANGLE)=SUMF_1(JPLAN,JE,JANGLE)+SF_1*W_A
+     &(JANGLE_A)
+              ELSEIF(JEL.EQ.0) THEN
+                SUMR_1(JPLAN,JE,1)=SUMR_1(JPLAN,JE,1)+SR_1*W(JANGLE)*W_A
+     &(JANGLE_A)
+                SUMF_1(JPLAN,JE,1)=SUMF_1(JPLAN,JE,1)+SF_1*W(JANGLE)*W_A
+     &(JANGLE_A)
+              ENDIF
+              IF(IDICHR.GT.0) THEN
+                IF(JEL.EQ.1) THEN
+                  SUMR_2(JPLAN,JE,JANGLE_A)=SUMR_2(JPLAN,JE,JANGLE_A)+SR
+     &_2*W(JANGLE)
+                  SUMF_2(JPLAN,JE,JANGLE_A)=SUMF_2(JPLAN,JE,JANGLE_A)+SF
+     &_2*W(JANGLE)
+                ELSEIF(JEL.EQ.2) THEN
+                  SUMR_2(JPLAN,JE,JANGLE)=SUMR_2(JPLAN,JE,JANGLE)+SR_2*W
+     &_A(JANGLE_A)
+                  SUMF_2(JPLAN,JE,JANGLE)=SUMF_2(JPLAN,JE,JANGLE)+SF_2*W
+     &_A(JANGLE_A)
+                ELSEIF(JEL.EQ.0) THEN
+                  SUMR_2(JPLAN,JE,1)=SUMR_2(JPLAN,JE,1)+SR_2*W(JANGLE)*W
+     &_A(JANGLE_A)
+                  SUMF_2(JPLAN,JE,1)=SUMF_2(JPLAN,JE,1)+SF_2*W(JANGLE)*W
+     &_A(JANGLE_A)
+                ENDIF
+              ENDIF
+C
+            ENDDO
+C
+          ENDDO
+          IF(I_EXT_A.NE.0) THEN
+            REWIND IUI9
+            READ(IUI9,1) NDUM
+            READ(IUI9,1) NDUM
+          ENDIF
+        ENDDO
+C
+        IF(I_EXT.NE.0) THEN
+          REWIND IUI6
+          READ(IUI6,1) NDUM
+          READ(IUI6,1) NDUM
+        ENDIF
+      ENDDO
+C
+      CLOSE(IUI6)
+      CLOSE(IUI9)
+      REWIND IUO2
+C
+      WRITE(IUO2,16) SPECTRO2,LIKE,SPECTRO,OUTDATA
+      IF((SPECTRO.NE.'APC').OR.(JEL.EQ.0)) THEN
+        WRITE(IUO2,19) ISPIN,IDICHR,I_SO,ISFLIP
+        WRITE(IUO2,18) NE,NPLAN,ISOM
+      ELSEIF(JEL.EQ.1) THEN
+        WRITE(IUO2,20) ISPIN_A,IDICHR_A,I_SO_A,ISFLIP_A,ICHKDIR_A,IPHI_A
+     &,ITHETA_A,IE_A
+        WRITE(IUO2,21) NPHI0,NTHETA0,NE,NPLAN,ISOM
+      ELSEIF(JEL.EQ.2) THEN
+        WRITE(IUO2,20) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE
+        WRITE(IUO2,21) NPHI0,NTHETA0,NE,NPLAN,ISOM
+      ENDIF
+C
+      DO JE=1,NE
+        DO JANGLE=1,NANGLE
+          IF(SPECTRO.EQ.'APC') THEN
+            IF(JEL.EQ.1) THEN
+              THETA=DTHETAA(JANGLE)
+              PHI=DPHIA(JANGLE)
+            ELSEIF(JEL.EQ.2) THEN
+              THETA=DTHETA(JANGLE)
+              PHI=DPHI(JANGLE)
+            ENDIF
+          ENDIF
+C
+          DO JPLAN=1,NPLAN
+            IF(IDICHR.EQ.0) THEN
+              IF((SPECTRO.NE.'APC').OR.(JEL.EQ.0)) THEN
+                WRITE(IUO2,33) JPLAN,ECIN(JE),SUMR_1(JPLAN,JE,JANGLE),SU
+     &MF_1(JPLAN,JE,JANGLE)
+              ELSE
+                WRITE(IUO2,34) JPLAN,THETA,PHI,ECIN(JE),SUMR_1(JPLAN,JE,
+     &JANGLE),SUMF_1(JPLAN,JE,JANGLE)
+              ENDIF
+            ELSE
+              IF((SPECTRO.NE.'APC').OR.(JEL.EQ.0)) THEN
+                WRITE(IUO2,43) JPLAN,ECIN(JE),SUMR_1(JPLAN,JE,JANGLE),SU
+     &MF_1(JPLAN,JE,JANGLE),SUMR_2(JPLAN,JE,JANGLE),SUMF_2(JPLAN,JE,JANG
+     &LE)
+              ELSE
+                WRITE(IUO2,44) JPLAN,THETA,PHI,ECIN(JE),SUMR_1(JPLAN,JE,
+     &JANGLE),SUMF_1(JPLAN,JE,JANGLE),SUMR_2(JPLAN,JE,JANGLE),SUMF_2(JPL
+     &AN,JE,JANGLE)
+              ENDIF
+            ENDIF
+          ENDDO
+C
+          IF(IDICHR.EQ.0) THEN
+            IF((SPECTRO.NE.'APC').OR.(JEL.EQ.0)) THEN
+              WRITE(IUO2,33) JVOL,ECIN(JE),SUMR_1(NPLAN+1,JE,JANGLE),SUM
+     &F_1(NPLAN+1,JE,JANGLE)
+              WRITE(IUO2,33) JTOT,ECIN(JE),SUMR_1(NPLAN+2,JE,JANGLE),SUM
+     &F_1(NPLAN+2,JE,JANGLE)
+            ELSE
+              WRITE(IUO2,34) JVOL,THETA,PHI,ECIN(JE),SUMR_1(NPLAN+1,JE,J
+     &ANGLE),SUMF_1(NPLAN+1,JE,JANGLE)
+              WRITE(IUO2,34) JTOT,THETA,PHI,ECIN(JE),SUMR_1(NPLAN+2,JE,J
+     &ANGLE),SUMF_1(NPLAN+2,JE,JANGLE)
+            ENDIF
+          ELSE
+            IF((SPECTRO.NE.'APC').OR.(JEL.EQ.0)) THEN
+              WRITE(IUO2,43) JVOL,ECIN(JE),SUMR_1(NPLAN+1,JE,JANGLE),SUM
+     &F_1(NPLAN+1,JE,JANGLE),SUMR_2(NPLAN+1,JE,JANGLE),SUMF_2(NPLAN+1,JE
+     &,JANGLE)
+              WRITE(IUO2,43) JTOT,ECIN(JE),SUMR_1(NPLAN+2,JE,JANGLE),SUM
+     &F_1(NPLAN+2,JE,JANGLE),SUMR_2(NPLAN+2,JE,JANGLE),SUMF_2(NPLAN+2,JE
+     &,JANGLE)
+            ELSE
+              WRITE(IUO2,44) JVOL,THETA,PHI,ECIN(JE),SUMR_1(NPLAN+1,JE,J
+     &ANGLE),SUMF_1(NPLAN+1,JE,JANGLE),SUMR_2(NPLAN+1,JE,JANGLE),SUMF_2(
+     &NPLAN+1,JE,JANGLE)
+              WRITE(IUO2,44) JTOT,THETA,PHI,ECIN(JE),SUMR_1(NPLAN+2,JE,J
+     &ANGLE),SUMF_1(NPLAN+2,JE,JANGLE),SUMR_2(NPLAN+2,JE,JANGLE),SUMF_2(
+     &NPLAN+2,JE,JANGLE)
+            ENDIF
+          ENDIF
+C
+        ENDDO
+      ENDDO
+C
+   1  FORMAT(13X,I4)
+   2  FORMAT(15X,F8.3,3X,F8.3,3X,E12.6)
+   3  FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6)
+   4  FORMAT(//,8X,'<<<<<<<<<<  DIMENSION OF THE ARRAYS TOO SMALL ','IN 
+     &THE WEIGHT_SUM SUBROUTINE - INCREASE NPM TO ',I3,'>>>>>>>>>>')
+   5  FORMAT(6X,I1,1X,I3,3X,I3)
+   8  FORMAT(I4,2X,I4,2X,I4,2X,I3,2X,I1)
+   9  FORMAT(9(2X,I1),2X,I2)
+  13  FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,F6.2,2X,F6.2,2X,E12.6,2X,E
+     &12.6)
+  15  FORMAT(2X,A3,11X,A13)
+  16  FORMAT(2X,A3,A5,1X,A3,2X,A13)
+  18  FORMAT(I4,2X,I3,2X,I1)
+  19  FORMAT(4(2X,I1))
+  20   FORMAT(8(2X,I1))
+  21   FORMAT(I4,2X,I4,2X,I4,2X,I3,2X,I1)
+  23  FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6,2X,E12.6,2X
+     &,E12.6)
+  24  FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,F6.2,2X,F6.2,2X,E12.6,2X,E
+     &12.6,2X,E12.6,2X,E12.6)
+  33  FORMAT(2X,I3,2X,F8.2,2X,E12.6,2X,E12.6)
+  34  FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6)
+  43  FORMAT(2X,I3,2X,F8.2,2X,E12.6,2X,E12.6,2X,E12.6,2X,E12.6)
+  44  FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6,2X,E12.6,2X
+     &,E12.6)
+C
+      RETURN
+C
+      END

src/msspec/spec/fortran/phd_mi_noso_nosp_nosym.mk

@@ -2,7 +2,7 @@ memalloc_src               := memalloc/dim_mod.f memalloc/modules.f memalloc/all
 cluster_gen_src            := $(wildcard cluster_gen/*.f)
 common_sub_src             := $(wildcard common_sub/*.f)
 renormalization_src        := $(wildcard renormalization/*.f)
-phd_mi_noso_nosp_nosym_src := $(wildcard phd_mi_noso_nosp_nosym/*.f)
+phd_mi_noso_nosp_nosym_src := $(filter-out phd_mi_noso_nosp_nosym/lapack_axb.f, $(wildcard phd_mi_noso_nosp_nosym/*.f))
 
 SRCS   = $(memalloc_src) $(cluster_gen_src) $(common_sub_src) $(renormalization_src) $(phd_mi_noso_nosp_nosym_src)
 MAIN_F = phd_mi_noso_nosp_nosym/main.f

src/msspec/spec/fortran/phd_se_noso_nosp_nosym/pathop.f

@@ -115,7 +115,7 @@ C  Renormalization of the path
 C
       IF(I_REN.GE.1) THEN
          COEF=COEF*C_REN(JORDP)
-         write(354,*) JORDP,C_REN(JORDP)
+C         write(354,*) JORDP,C_REN(JORDP)
       ENDIF
 C
 C   Call of the subroutines used for the R-A termination matrix

src/msspec/tests.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -16,8 +17,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/tests.py
-# Last modified: ven. 10 avril 2020 17:33:28
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
-# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 import os

src/msspec/utils.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -18,8 +19,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/utils.py
-# Last modified: Thu, 06 Oct 2022 18:19:16 +0200
+# Last modified: Thu, 06 Oct 2022 18:27:24 +0200
-# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1665073156 +0200
+# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1665073644 +0200
 
 
 """
@@ -70,7 +71,7 @@ class ForeignPotential(object):
         self.phagen_data = {'types': []}
 
     def write(self, filename, prototypical_atoms):
-        LOGGER.debug(f"Writing Phagen input potential file: {filename}")
+        LOGGER.debug("Writing Phagen input potential file: {}".format(filename))
 
         def DEPRECATEDappend_atom_potential(atom):
             Z = atom.number
@@ -81,8 +82,8 @@ class ForeignPotential(object):
                     itypes.append(i)
             # Check now that we have only one type in the list
             # otherwise we do not know yet how to deal with this.
-            assert len(itypes) > 0, f"Cannot find the data for atom with Z={Z}"
+            assert len(itypes) > 0, "Cannot find the data for atom with Z={}".format(Z)
-            assert len(itypes) == 1, f"Too many datasets for atom with Z={Z}"
+            assert len(itypes) == 1, "Too many datasets for atom with Z={}".format(Z)
             # So far so good, let's write the block
             t = self.phagen_data['types'][itypes[0]]
             s = "{:<7d}{:<10d}{:1.4f}\n".format(
@@ -95,7 +96,7 @@ class ForeignPotential(object):
         def append_atom_potential(atom):
             line_fmt = "{:+1.8e} " * 4 + "\n"
             atom_type = atom.get('atom_type')
-            assert atom_type != None, f"Unable get the atom type!"
+            assert atom_type != None, "Unable get the atom type!"
             for t in self.phagen_data['types']:
                 if t['atom_type'] == atom_type:
                     s = "{:<7d}{:<10d}{:1.4f}\n".format(
@@ -138,7 +139,7 @@ class SPRKKRPotential(ForeignPotential):
         self.potfile = potfile
         self.load_sprkkr_atom_types()
         for f in exported_files:
-            LOGGER.info(f"Loading file {f}...")
+            LOGGER.info("Loading file {}...".format(f))
             # get the IT from the filename
             m=re.match('SPRKKR-IT_(?P<IT>\d+)-PHAGEN.*', os.path.basename(f))
             it = int(m.group('IT'))
@@ -192,7 +193,7 @@ class SPRKKRPotential(ForeignPotential):
             return data
 
         # load info in *.pot file
-        LOGGER.info(f"Loading SPRKKR *.pot file {self.potfile}...")
+        LOGGER.info("Loading SPRKKR *.pot file {}...".format(self.potfile))
         with open(self.potfile, 'r') as fd:
             content = fd.read()
 
@@ -233,7 +234,7 @@ class SPRKKRPotential(ForeignPotential):
                     IT = occupation['ITOQ']
                     atom = self.atoms[i]
                     atom.set('atom_type', IT)
-            LOGGER.debug(f"Site #{IQ} is type #{IT}, atom {atom}")
+            LOGGER.debug("Site #{} is type #{}, atom {}".format(IQ, IT, atom))
 
 
 
@@ -314,34 +315,13 @@ def cut_cylinder(atoms, axis="z", radius=None):
     :return: The modified atom cluster
     :rtype: ase.Atoms
     """
-    if radius is None:
+    if axis not in ('z',):
-        raise ValueError("radius not set")
+        raise ValueError("axis value != 'z' is not supported yet.")
-
+    X, Y, Z = atoms.positions.T
-    new_atoms = atoms.copy()
+    R = np.sqrt(X**2 + Y **2)
-
+    T = np.arctan2(Y, X)
-    dims = {"x": 0, "y": 1, "z": 2}
+    i = np.where(R <= radius)[0]
-    if axis in dims:
+    return atoms[i]
-        axis = dims[axis]
-    else:
-        raise ValueError("axis not valid, must be 'x','y', or 'z'")
-
-    del_list = []
-    for index, position in enumerate(new_atoms.positions):
-        # calculating the distance of the atom to the given axis
-        r = 0
-        for dim in range(3):
-            if dim != axis:
-                r = r + position[dim]**2
-        r = np.sqrt(r)
-
-        if r > radius:
-            del_list.append(index)
-
-    del_list.reverse()
-    for index in del_list:
-        del new_atoms[index]
-
-    return new_atoms
 
 
 def cut_cone(atoms, radius, z=0):
@@ -429,11 +409,15 @@ def cut_plane(atoms, x=None, y=None, z=None):
 
     dim_values = np.array(dim_values)
 
-    def constraint(coordinates):
+    X, Y, Z = atoms.positions.T
-        return np.all(np.logical_and(coordinates >= dim_values[:, 0],
+    i0 = np.where(X                     >= dim_values[0, 0])[0]
-                                     coordinates <= dim_values[:, 1]))
+    i1 = np.where(X[i0]                 <= dim_values[0, 1])[0]
+    i2 = np.where(Y[i0][i1]             >= dim_values[1, 0])[0]
+    i3 = np.where(Y[i0][i1][i2]         <= dim_values[1, 1])[0]
+    i4 = np.where(Z[i0][i1][i2][i3]     >= dim_values[2, 0])[0]
+    i5 = np.where(Z[i0][i1][i2][i3][i4] <= dim_values[2, 1])[0]
+    indices = np.arange(len(atoms))[i0][i1][i2][i3][i4][i5]
 
-    indices = np.where(list(map(constraint, atoms.positions)))[0]
     return atoms[indices]
 
 

src/msspec/version.py

@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# coding: utf-8
 #
 # Copyright © 2016-2020 - Rennes Physics Institute
 #
@@ -16,8 +17,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/version.py
-# Last modified: Thu, 06 Oct 2022 18:19:16 +0200
+# Last modified: Wed, 26 Oct 2022 17:15:24 +0200
-# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1665073156 +0200
+# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1666797324 +0200
 
 
 import os
@@ -26,28 +27,27 @@ from importlib.metadata import version
 import subprocess
 
 # find the version number
-# 1- If it fails, try to read it from the distribution file
+# 1- Try to read it from the git info
-# 2- Try to read it from the git info
+# 2- If it fails, try to read it from the VERSION file
-# 3- If it fails, try to read it from the VERSION file
+# 3- If it fails, try to read it from the distribution file
 
 PKGNAME = 'msspec'
 
 try:
-    __version__ = version(PKGNAME)
+    cmd = ["git describe|sed 's/-\([0-9]\+\)-.*/.dev\\1/g'"]
+    result = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL, shell=True)
+    __version__ = result.stdout.decode('utf-8').strip()
+    if __version__ == "":
+        raise
 except Exception as err:
-    try:
-        p = subprocess.run(["git", "describe"], capture_output=True, text=True)
-        if p.stdout not in ("", None):
-            __version__ = p.stdout.strip()
-        else:
-            raise NameError("git describe failed!")
-    except Exception as err:
     try:
         thisfile_path = os.path.abspath(__file__)
         thisfile_dir = os.path.dirname(thisfile_path)
-            versionfile = os.path.join(thisfile_dir, "../VERSION")
+        versionfile = os.path.join(thisfile_dir, "./VERSION")
         with open(versionfile, "r") as fd:
             __version__ = fd.readline().strip()
     except Exception as err:
-            print("Unable to get the version number!")
+        try:
-            __version__ = "9.9.9"
+            __version__ = version(PKGNAME)
+        except Exception as err:
+            __version__ = "0.0.0"

src/options.mk

@@ -1,6 +1,6 @@
-PYTHON   = python
+PYTHON   = python3
 PYMAJ    = 3
-PYMIN    = 6
+PYMIN    = 5
 
 FC       = gfortran
 F2PY     = f2py3 --f77exec=$(FC) --f90exec=$(FC)
@@ -31,7 +31,7 @@ IFORT_FFLAGS_DBG =
 ################################################################################
 # F2PY CONFIGURATION                                                           #
 ################################################################################
-F2PYFLAGS     = --opt=-O2
+F2PYFLAGS     = --opt=-O2 -llapack -larpack 
 F2PYFLAGS_DBG = --debug-capi --debug
 ################################################################################
 
@@ -41,7 +41,7 @@ F2PYFLAGS_DBG = --debug-capi --debug
 # /!\ DO NOT EDIT BELOW THAT LINE (unlesss you know what you're doing...)      #
 # CORE CONFIGURATION                                                           #
 ################################################################################
-VERSION:=$(shell git describe)
+VERSION:=$(shell git describe|sed 's/-\([0-9]\+\)-.*/.dev\1/g')
 VENV_PATH := $(INSTALL_PREFIX)/src/msspec_venv_$(VERSION)
 
 
@@ -103,7 +103,7 @@ endif
 
 FFLAGS = $($(PREFIX)_FFLAGS$(SUFFIX))
 
-OBJS = $(addprefix $(BUILDDIR)/, $(patsubst %.f,%.o, $(filter-out $(MAIN_F), $(SRCS))))
+OBJS = $(addprefix $(BUILDDIR)/, $(patsubst %.f90,%.o, $(patsubst %.f,%.o, $(filter-out $(MAIN_F), $(SRCS)))))
 
 .PHONY: clean obj all info
 
@@ -141,6 +141,12 @@ $(BUILDDIR)/%.o: %.f
 	$(FC) $(FFLAGS) -J $(BUILDDIR) -I $(BUILDDIR) -fPIC -o $@ -c $^ $(OUPUT_REDIRECTION)
 
 
+$(BUILDDIR)/%.o: %.f90
+	@echo "Compiling $@..."
+	mkdir -p $(basename $@)
+	$(FC) $(FFLAGS) -J $(BUILDDIR) -I $(BUILDDIR) -fPIC -o $@ -c $^ $(OUPUT_REDIRECTION)
+
+
 $(SO): $(OBJS) $(MAIN_F)
 	@echo "building Python binding $@..."      
 	mkdir -p $(BUILDDIR)

src/pip.freeze

@@ -2,7 +2,7 @@ ase
 h5py
 ipython
 lxml
-matplotlib==3.4.3
+matplotlib
 numpy
 Pint
 pandas
@@ -10,3 +10,5 @@ pycairo
 scipy
 setuptools-scm
 terminaltables
+wheel
+wxPython

src/pyproject.toml

@@ -0,0 +1,3 @@
+[build-system]
+requires = ["setuptools>=45", "setuptools_scm[toml]>=6.2"]
+build-backend = "setuptools.build_meta"

src/setup.cfg

@@ -0,0 +1,55 @@
+[metadata]
+name = msspec
+version = attr: msspec.version.__version__
+author = Didier Sébilleau, Sylvain Tricot
+author_email = sylvain.tricot@univ-rennes1.fr
+url = https://msspec.cnrs.fr
+description = A multiple scattering package for sepectroscopies using electrons to probe materials
+long_description = MsSpec is a Fortran package to compute the
+                   cross-section of several spectroscopies involving one (or more)
+                   electron(s) as the probe. This package provides a python interface to
+                   control all the steps of the calculation.
+
+                   Available spectroscopies:
+                       * Photoelectron diffraction
+                       * Auger electron diffraction
+                       * Low energy electron diffraction
+                       * X-Ray absorption spectroscopy
+                       * Auger Photoelectron coincidence spectroscopy
+                       * Computation of the spectral radius""",
+keywords = spectroscopy atom electron photon multiple scattering
+license = GPL
+classifiers =
+    Development Status :: 3 - Alpha
+    Environment :: Console
+    Intended Audience :: Science/Research
+    License :: OSI Approved :: GNU General Public License (GPL)
+    Natural Language :: English
+    Operating System :: Microsoft :: Windows :: Windows 10
+    Operating System :: POSIX :: Linux
+    Operating System :: MacOS :: MacOS X
+    Programming Language :: Fortran
+    Programming Language :: Python :: 3 :: Only
+    Topic :: Scientific/Engineering :: Physics
+
+[options]
+packages = find:
+zip_safe = False
+install_requires =
+    setuptools_scm
+    ase
+    h5py
+    ipython
+    lxml
+    matplotlib
+    numpy
+    Pint
+    pandas
+    pycairo
+    scipy
+    terminaltables
+
+[options.package_data]
+msspec.phagen = fortran/*.so
+msspec.spec = fortran/*.so
+msspec = VERSION