epsi
/
msspec_python3
3
0
Fork 0
Code Issues 6 Pull Requests Releases Wiki Activity

Compare commits

base: epsi:feature/arnoldi
Branches Tags
epsi:master
epsi:devel
epsi:feature/LEED
epsi:feature/arnoldi
epsi:feature/aed
epsi:1.7.post13
epsi:1.7.post12
epsi:1.7.post11
epsi:1.7
epsi:1.6
epsi:1.5
epsi:1.4
epsi:1.3
..
compare: epsi:master
Branches Tags
epsi:devel
epsi:feature/LEED
epsi:feature/arnoldi
epsi:master
epsi:feature/aed
epsi:1.7.post13
epsi:1.7.post12
epsi:1.7.post11
epsi:1.7
epsi:1.6
epsi:1.5
epsi:1.4
epsi:1.3

No commits in common. "feature/arnoldi" and "master" have entirely different histories.
feature/ar ... master

62 changed files with 175 additions and 13025 deletions
Show Stats Expand all files Collapse all files

100
Dockerfile
View File

@ -1,102 +1,24 @@
# Get the base Python image
FROM alpine:edge AS builder
# Variables
ARG branch="devel"
ARG login="" password=""
ARG folder=/opt/msspec user=msspec
FROM python:latest
# Install system dependencies
# tools
RUN apk add bash git make gfortran python3
# headers
RUN apk add python3-dev lapack-dev musl-dev
# python packages
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
RUN apk add --no-cache -X http://dl-cdn.alpinelinux.org/alpine/edge/testing py3-wxpython
RUN pip install ase pint terminaltables ipython
# for GUI
RUN apk add ttf-droid adwaita-icon-theme
RUN apk add build-base
# Fetch the code
RUN mkdir -p ${folder}/code
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
RUN apt-get update && apt-get install -y virtualenv gfortran libgtk-3-dev nano
# Add a non-privileged user
RUN adduser -D -s /bin/bash -h ${folder} ${user}
RUN useradd -ms /bin/bash -d /opt/msspec msspec
# Set the working directory in the container
USER ${user}
USER msspec
RUN mkdir -p /opt/msspec/code
WORKDIR /opt/msspec/code
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
# Fetch the code
RUN git clone https://git.ipr.univ-rennes1.fr/epsi/msspec_python3.git .
#COPY --chown=msspec:msspec . .
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
#COPY --from=builder ${folder}/.local ${folder}/.local
#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
ENV PATH=/opt/msspec/.local/bin:$PATH
RUN make install VERBOSE=1
# Run the msspec frontend command on startup
ENTRYPOINT ["msspec"]

4
Jenkinsfile vendored
View File

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

39
Makefile
View File

@ -1,7 +1,7 @@
include src/options.mk
.PHONY: pybinding install devel venv doc clean _attrdict
.PHONY: pybinding install devel venv doc clean
pybinding:
@ -10,23 +10,20 @@ pybinding:
venv:
ifeq ($(NO_VENV),0)
# @virtualenv --python=$(PYTHON_EXE) --prompt="(msspec-$(VERSION)) " $(VENV_PATH)
$(PYTHON) -m venv $(VENV_PATH)
$(INSIDE_VENV) python -m ensurepip --upgrade
$(INSIDE_VENV) pip install -r src/pip.freeze
@virtualenv --python=$(PYTHON_EXE) --prompt="(msspec-$(VERSION)) " $(VENV_PATH)
$(INSIDE_VENV) \
wget https://bootstrap.pypa.io/get-pip.py && \
python get-pip.py && \
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 --force-reinstall src/dist/msspec-$(VERSION)*.whl
@+$(INSIDE_VENV) pip install src/dist/msspec-$(VERSION).tar.gz
@echo "Do not forget to check that $(INSTALL_PREFIX)/bin is set in your \$$PATH"
@ -40,28 +37,20 @@ light: venv
@$(INSIDE_VENV) pip install src/
_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
_build_wx/wxPython.target:
@$(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 sip; \
pip install requests; \
python build.py dox etg --nodoc sip build bdist_wheel; \
ln -sf `readlink -f dist/wxPython*.whl` ../; \
fi;

3
doc/source/tutorials/RhO/RhO.py
View File

@ -18,8 +18,7 @@ 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,
malloc={'NPH_M': 8000})
data = calc.get_theta_phi_scan(level='1s', kinetic_energy=723, data=data)
dset = data[-1]
dset.title = "{:d}) z = {:.2f} angstroms".format(zi, z0)

6
src/MANIFEST.in Normal file
View File

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

13
src/Makefile
View File

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

9
src/msspec.sh.template
View File

@ -2,11 +2,12 @@
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 msspec!!"
echo "ERROR: Unable to find version $VERSION of msspec!!"
exit 1
fi
@ -14,10 +15,6 @@ 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]"
@ -95,7 +92,7 @@ while getopts "hvil:p:eu" option; do
;;
u) uninstall
;;
v) show_version
v) echo $VERSION
;;
*|h) show_help
;;

5
src/msspec/__init__.py
View File

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

30
src/msspec/calcio.py
View File

@ -1,25 +1,5 @@
#!/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/calcio.py
# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
# vim: set et ts=4 sw=4 fdm=indent mouse=a cc=+1 tw=80:
"""
Module calcio
@ -930,7 +910,7 @@ class SpecIO(object):
if content != old_content:
with open(filename, 'w') as fd:
fd.write(content)
LOGGER.debug("Writing Spec input file written in {}".format(filename))
LOGGER.debug(f"Writing Spec input file written in {filename}")
modified = True
return modified
@ -1275,13 +1255,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 + '{:02d}'.format(i) + '.txt')
results = np.loadtxt(prefix + f'{i:02d}' + '.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',
'experiment_{}.txt'.format(suffix)))
f'experiment_{suffix}.txt'))
calc_ren = np.loadtxt(os.path.join('calc', 'div',
'calculation{:d}_{}.txt'.format(index,suffix)))
f'calculation{index:d}_{suffix}.txt'))
return data, exp_ren, calc_ren

33
src/msspec/calculator.py
View File

@ -1,5 +1,4 @@
#!/usr/bin/env python
# coding: utf-8
#
# Copyright © 2016-2020 - Rennes Physics Institute
#
@ -17,8 +16,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/calculator.py
# Last modified: Tue, 25 Oct 2022 16:21:38 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1666707698 +0200
# Last modified: ven. 10 avril 2020 17:19:24
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
"""
@ -95,10 +94,8 @@ 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
@ -407,8 +404,6 @@ 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,
@ -419,8 +414,6 @@ 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,
@ -752,7 +745,7 @@ class _PED(_MSCALCULATOR):
view = dset.add_view("E = {:.2f} eV".format(ke), title=title,
xlabel=xlabel, ylabel=ylabel,
projection='stereo', colorbar=True, autoscale=False)
projection='stereo', colorbar=True, autoscale=True)
view.select('theta', 'phi', 'cross_section')
@ -989,8 +982,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', 'arnoldi'):
LOGGER.error("Only the 'inversion', 'power' or 'arnoldi' algorithms "
if algorithm not in ('inversion', 'power'):
LOGGER.error("Only the 'inversion' or the 'power' algorithms "
"are supported in EIG spectroscopy mode")
exit(1)
self.iodata = iodata.Data('EIG Simulation')
@ -1139,7 +1132,7 @@ class RFACTOR(object):
for i in range(noif):
X, Y = args[2*i], args[2*i+1]
fname = os.path.join('calc',
'calculation{:d}.txt'.format(self.stack_count))
f'calculation{self.stack_count:d}.txt')
# And save to the working space
np.savetxt(fname, np.transpose([X, Y]))
self.stack_count += 1
@ -1147,7 +1140,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', 'calculation{:d}.txt'.format(i))
fname = os.path.join('calc', f'calculation{i:d}.txt')
self._params.calc_filename.append(fname)
# Write the input file
@ -1242,23 +1235,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(**{"calc{:d}".format(i): ycalc})
dset_rfc.add_columns(**{'variable_set{:d}'.format(i): rfc})
dset_values.add_columns(**{f"calc{i:d}": ycalc})
dset_rfc.add_columns(**{f'variable_set{i:d}': rfc})
# Plot the curves
view_values.select('x', 'yref', legend='Reference values')
title = ''
for k,v in self.best_values.items():
title += '{}={} '.format(k, v)
view_values.select('x', "calc{:d}".format(self.index),
title += f'{k}={v} '
view_values.select('x', f"calc{self.index:d}",
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', 'variable_set{:d}'.format(i),
legend="variables set #{:d}".format(i))
view_rfc.select('rfactor_number', f'variable_set{i:d}',
legend=f"variables set #{i:d}")
# Save the parameters
for p in self.get_parameters():
bundle = {'group': str(p.group),

5
src/msspec/cli.py
View File

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

5
src/msspec/create_tests_results.py
View File

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

23
src/msspec/data/__init__.py
View File

@ -1,24 +1,5 @@
#!/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/__init__.py
# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
# -*- encoding: utf-8 -*-
# vim: set fdm=indent ts=4 sw=4 sts=4 et ai tw=80 cc=+0 mouse=a nu : #
from .electron_be import electron_be

20
src/msspec/data/electron_be.py
View File

@ -1,24 +1,4 @@
#!/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

13
src/msspec/iodata.py
View File

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

15
src/msspec/iodata_gi.py
View File

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

15
src/msspec/iodata_wx.py
View File

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

26
src/msspec/looper.py
View File

@ -1,24 +1,6 @@
#!/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/looper.py
# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
# coding: utf8
# -*- encoding: future_fstrings -*-
# vim: set et sw=4 ts=4 nu tw=79 cc=+1:
from collections import OrderedDict
from functools import partial
@ -39,7 +21,7 @@ class Variable:
self.doc = doc
def __repr__(self):
return "<Variable(\'{}\')>".format(self.name)
return f"<Variable(\'{self.name}\')>"
class Sweep:
def __init__(self, key, comments="", unit=None,

5
src/msspec/misc.py
View File

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

26
src/msspec/parameters.py
View File

@ -1,5 +1,4 @@
#!/usr/bin/env python
# coding: utf-8
#
# Copyright © 2016-2020 - Rennes Physics Institute
#
@ -19,8 +18,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/parameters.py
# Last modified: Tue, 15 Feb 2022 15:37:28 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>
# Last modified: ven. 10 avril 2020 17:31:50
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
"""
@ -770,8 +769,7 @@ class GlobalParameters(BaseParameters):
Parameter('algorithm', types=str, allowed_values=('expansion',
'inversion',
'correlation',
'power',
'arnoldi'),
'power'),
default='expansion', doc="""
You can choose the algorithm used for the computation of the scattering path operator.
@ -779,7 +777,6 @@ 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
@ -841,17 +838,6 @@ 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.')
@ -2026,20 +2012,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 "
"\"{}\"".format(p.value))
f"\"{p.value}\"")
def bind_rescale(self, p):
self.compcurve_parameters.general_iscale = int(p.value)
state = "deactivated"
if p.value:
state = "activated"
LOGGER.info("Curve Comparison: Rescaling of data {}".format(state))
LOGGER.info(f"Curve Comparison: Rescaling of data {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 "
"set to \"{}\"".format(p.value))
f"set to \"{p.value}\"")

12
src/msspec/spec/fortran/Makefile
View File

@ -1,6 +1,6 @@
.PHONY: all phd_se phd_mi phd_ce eig_mi eig_pw eig_ar comp_curve clean
.PHONY: all phd_se phd_mi eig_mi eig_pw comp_curve clean
all: phd_se phd_mi phd_ce eig_mi eig_pw eig_ar comp_curve
all: phd_se phd_mi eig_mi eig_pw comp_curve
phd_se:
@+$(MAKE) -f phd_se_noso_nosp_nosym.mk all
@ -8,26 +8,18 @@ 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 $@

250
src/msspec/spec/fortran/eig/ar/arnoldi_mod.f90
View File

@ -1,250 +0,0 @@
!==============================================================================!
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
!==============================================================================!

1558
src/msspec/spec/fortran/eig/ar/do_main.f
View File

File diff suppressed because it is too large Load Diff

291
src/msspec/spec/fortran/eig/ar/eig_mat_ar.f90
View File

@ -1,291 +0,0 @@
!==============================================================================!
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
!==============================================================================!

103
src/msspec/spec/fortran/eig/ar/eigdif_ar.f
View File

@ -1,103 +0,0 @@
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

22
src/msspec/spec/fortran/eig/ar/main.f
View File

@ -1,22 +0,0 @@
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

4
src/msspec/spec/fortran/eig/common/eig_mat_ms.f
View File

@ -178,10 +178,6 @@ 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

37
src/msspec/spec/fortran/eig/common/save_eigenvalues.f
View File

@ -1,37 +0,0 @@
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

14
src/msspec/spec/fortran/eig_ar.mk
View File

@ -1,14 +0,0 @@
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

3
src/msspec/spec/fortran/eig_mi.mk
View File

@ -2,8 +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)
#eig_common_src := $(wildcard eig/common/*.f)
eig_common_src := $(filter-out eig/common/lapack_eig.f, $(wildcard eig/common/*.f))
eig_common_src := $(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)

6
src/msspec/spec/fortran/memalloc/allocation.f
View File

@ -25,9 +25,6 @@
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
@ -139,7 +136,6 @@
CALL ALLOC_OUTUNITS()
CALL ALLOC_PARCAL()
CALL ALLOC_PARCAL_A()
CALL ALLOC_Q_ARRAY()
CALL ALLOC_RELADS()
CALL ALLOC_RELAX()
CALL ALLOC_RENORM()
@ -177,7 +173,6 @@
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()
@ -191,7 +186,6 @@
CALL ALLOC_SPECTRUM()
CALL ALLOC_DIRECT()
CALL ALLOC_DIRECT_A()
CALL ALLOC_GAUNT_C()
CALL ALLOC_PATH()
CALL ALLOC_ROT()
CALL ALLOC_ROT_CUB()

4
src/msspec/spec/fortran/memalloc/dim_mod.f
View File

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

42
src/msspec/spec/fortran/memalloc/modules.f
View File

@ -192,20 +192,6 @@ 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
@ -431,20 +417,6 @@ 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
@ -806,20 +778,6 @@ 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

11
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym.mk
View File

@ -1,11 +0,0 @@
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

41
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/cmngr.f
View File

@ -1,41 +0,0 @@
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

46
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/coefpq.f
View File

@ -1,46 +0,0 @@
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

47
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/corexp.f
View File

@ -1,47 +0,0 @@
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

19
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/corexp1.f
View File

@ -1,19 +0,0 @@
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

121
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/coumat.f
View File

@ -1,121 +0,0 @@
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

85
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/dwsph.f
View File

@ -1,85 +0,0 @@
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

26
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/facdif.f
View File

@ -1,26 +0,0 @@
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

113
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/facdif1.f
View File

@ -1,113 +0,0 @@
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

126
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/gaunt_st.f
View File

@ -1,126 +0,0 @@
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

5123
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/lapack_axb.f
View File

File diff suppressed because it is too large Load Diff

21
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/main.f
View File

@ -1,21 +0,0 @@
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

1700
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/main_phd_ns_ce.f
View File

File diff suppressed because it is too large Load Diff

280
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/mpis.f
View File

@ -1,280 +0,0 @@
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

165
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/ms_cor.f
View File

@ -1,165 +0,0 @@
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

1136
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/phddif_ce.f
View File

File diff suppressed because it is too large Load Diff

106
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/plotfd.f
View File

@ -1,106 +0,0 @@
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

769
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/treat_phd.f
View File

@ -1,769 +0,0 @@
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

335
src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/weight_sum.f
View File

@ -1,335 +0,0 @@
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

2
src/msspec/spec/fortran/phd_mi_noso_nosp_nosym.mk
View File

@ -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 := $(filter-out phd_mi_noso_nosp_nosym/lapack_axb.f, $(wildcard phd_mi_noso_nosp_nosym/*.f))
phd_mi_noso_nosp_nosym_src := $(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

2
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/pathop.f
View File

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

5
src/msspec/tests.py
View File

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

66
src/msspec/utils.py
View File

@ -1,5 +1,4 @@
#!/usr/bin/env python
# coding: utf-8
#
# Copyright © 2016-2020 - Rennes Physics Institute
#
@ -19,8 +18,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:27:24 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1665073644 +0200
# Last modified: Thu, 06 Oct 2022 18:19:16 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1665073156 +0200
"""
@ -71,7 +70,7 @@ class ForeignPotential(object):
self.phagen_data = {'types': []}
def write(self, filename, prototypical_atoms):
LOGGER.debug("Writing Phagen input potential file: {}".format(filename))
LOGGER.debug(f"Writing Phagen input potential file: {filename}")
def DEPRECATEDappend_atom_potential(atom):
Z = atom.number
@ -82,8 +81,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, "Cannot find the data for atom with Z={}".format(Z)
assert len(itypes) == 1, "Too many datasets for atom with Z={}".format(Z)
assert len(itypes) > 0, f"Cannot find the data for atom with Z={Z}"
assert len(itypes) == 1, f"Too many datasets for atom with Z={Z}"
# So far so good, let's write the block
t = self.phagen_data['types'][itypes[0]]
s = "{:<7d}{:<10d}{:1.4f}\n".format(
@ -96,7 +95,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, "Unable get the atom type!"
assert atom_type != None, f"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(
@ -139,7 +138,7 @@ class SPRKKRPotential(ForeignPotential):
self.potfile = potfile
self.load_sprkkr_atom_types()
for f in exported_files:
LOGGER.info("Loading file {}...".format(f))
LOGGER.info(f"Loading file {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'))
@ -193,7 +192,7 @@ class SPRKKRPotential(ForeignPotential):
return data
# load info in *.pot file
LOGGER.info("Loading SPRKKR *.pot file {}...".format(self.potfile))
LOGGER.info(f"Loading SPRKKR *.pot file {self.potfile}...")
with open(self.potfile, 'r') as fd:
content = fd.read()
@ -234,7 +233,7 @@ class SPRKKRPotential(ForeignPotential):
IT = occupation['ITOQ']
atom = self.atoms[i]
atom.set('atom_type', IT)
LOGGER.debug("Site #{} is type #{}, atom {}".format(IQ, IT, atom))
LOGGER.debug(f"Site #{IQ} is type #{IT}, atom {atom}")
@ -315,13 +314,34 @@ def cut_cylinder(atoms, axis="z", radius=None):
:return: The modified atom cluster
:rtype: ase.Atoms
"""
if axis not in ('z',):
raise ValueError("axis value != 'z' is not supported yet.")
X, Y, Z = atoms.positions.T
R = np.sqrt(X**2 + Y **2)
T = np.arctan2(Y, X)
i = np.where(R <= radius)[0]
return atoms[i]
if radius is None:
raise ValueError("radius not set")
new_atoms = atoms.copy()
dims = {"x": 0, "y": 1, "z": 2}
if axis in dims:
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):
@ -409,15 +429,11 @@ def cut_plane(atoms, x=None, y=None, z=None):
dim_values = np.array(dim_values)
X, Y, Z = atoms.positions.T
i0 = np.where(X >= dim_values[0, 0])[0]
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]
def constraint(coordinates):
return np.all(np.logical_and(coordinates >= dim_values[:, 0],
coordinates <= dim_values[:, 1]))
indices = np.where(list(map(constraint, atoms.positions)))[0]
return atoms[indices]

36
src/msspec/version.py
View File

@ -1,5 +1,4 @@
#!/usr/bin/env python
# coding: utf-8
#
# Copyright © 2016-2020 - Rennes Physics Institute
#
@ -17,8 +16,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/version.py
# Last modified: Wed, 26 Oct 2022 17:15:24 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1666797324 +0200
# Last modified: Thu, 06 Oct 2022 18:19:16 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1665073156 +0200
import os
@ -27,27 +26,28 @@ from importlib.metadata import version
import subprocess
# find the version number
# 1- 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 distribution file
# 1- If it fails, try to read it from the distribution file
# 2- Try to read it from the git info
# 3- If it fails, try to read it from the VERSION file
PKGNAME = 'msspec'
try:
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
__version__ = version(PKGNAME)
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")
with open(versionfile, "r") as fd:
__version__ = fd.readline().strip()
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:
__version__ = version(PKGNAME)
thisfile_path = os.path.abspath(__file__)
thisfile_dir = os.path.dirname(thisfile_path)
versionfile = os.path.join(thisfile_dir, "../VERSION")
with open(versionfile, "r") as fd:
__version__ = fd.readline().strip()
except Exception as err:
__version__ = "0.0.0"
print("Unable to get the version number!")
__version__ = "9.9.9"

16
src/options.mk
View File

@ -1,6 +1,6 @@
PYTHON = python3
PYTHON = python
PYMAJ = 3
PYMIN = 5
PYMIN = 6
FC = gfortran
F2PY = f2py3 --f77exec=$(FC) --f90exec=$(FC)
@ -31,7 +31,7 @@ IFORT_FFLAGS_DBG =
################################################################################
# F2PY CONFIGURATION #
################################################################################
F2PYFLAGS = --opt=-O2 -llapack -larpack
F2PYFLAGS = --opt=-O2
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|sed 's/-\([0-9]\+\)-.*/.dev\1/g')
VERSION:=$(shell git describe)
VENV_PATH := $(INSTALL_PREFIX)/src/msspec_venv_$(VERSION)
@ -103,7 +103,7 @@ endif
FFLAGS = $($(PREFIX)_FFLAGS$(SUFFIX))
OBJS = $(addprefix $(BUILDDIR)/, $(patsubst %.f90,%.o, $(patsubst %.f,%.o, $(filter-out $(MAIN_F), $(SRCS)))))
OBJS = $(addprefix $(BUILDDIR)/, $(patsubst %.f,%.o, $(filter-out $(MAIN_F), $(SRCS))))
.PHONY: clean obj all info
@ -141,12 +141,6 @@ $(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)

4
src/pip.freeze
View File

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

3
src/pyproject.toml
View File

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

55
src/setup.cfg
View File

@ -1,55 +0,0 @@
[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

BIN
thirdparty/attrdict-2.0.1.tar.gz vendored
View File

Binary file not shown.