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

Compare commits

base: epsi:devel
Branches Tags
epsi:master
epsi:devel
epsi:feature/LEED
epsi:feature/arnoldi
epsi:feature/aed
epsi:1.7.post14
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:master
epsi:feature/LEED
epsi:feature/arnoldi
epsi:feature/aed
epsi:1.7.post14
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. "devel" and "master" have entirely different histories.
devel ... master

25 changed files with 398 additions and 25856 deletions
Show Stats Expand all files Collapse all files

5
src/msspec/__init__.py
View File

@ -17,8 +17,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/__init__.py
# Last modified: Wed, 18 Jun 2025 13:49:16 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
import ase
@ -38,6 +38,5 @@ def init_msspec():
ase.atom.names['RA_cut_off'] = ('RA_cuts_off', 1)
ase.atom.names['atom_type'] = ('atom_types', None)
ase.atoms.Atoms.absorber = None
ase.atoms.Atoms.emitter = property(lambda self: self.absorber, lambda self,i: setattr(self, "absorber", i))
init_msspec()

16
src/msspec/calculator.py
View File

@ -17,8 +17,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/calculator.py
# Last modified: Mon, 23 Jun 2025 13:58:23 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
# Last modified: Tue, 25 Oct 2022 16:21:38 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1666707698 +0200
"""
@ -113,14 +113,13 @@ class _MSCALCULATOR(Calculator):
def __init__(self, spectroscopy='PED', algorithm='expansion',
polarization=None, dichroism=None, spinpol=False,
folder='./calc', txt='-', **kwargs):
self.txt = txt
stdout = sys.stdout
#if isinstance(txt, str) and txt != '-':
# stdout = open(txt, 'w')
if isinstance(txt, str) and txt != '-':
stdout = open(txt, 'w')
#elif isinstance(txt, buffer):
# stdout = txt
#elif txt == None:
# stdout = open('/dev/null', 'a')
elif txt == None:
stdout = open('/dev/null', 'a')
#set_log_output(stdout)
########################################################################
LOGGER.debug('Initialization of %s', self.__class__.__name__)
@ -173,7 +172,6 @@ class _MSCALCULATOR(Calculator):
self.global_parameters, self.phagen_parameters, self.spec_parameters)
# initialize all parameters with defaults values
self.spec_parameters.output_log = txt
LOGGER.info("Set default values =========================================")
for p in (list(self.global_parameters) +
list(self.muffintin_parameters) +
@ -381,7 +379,7 @@ class _MSCALCULATOR(Calculator):
'LI_M' : get_li(self.spec_parameters.extra_level) + 1,
'NEMET_M' : 1,
'NO_ST_M' : self.spec_parameters.calc_no,
'NDIF_M' : 18,
'NDIF_M' : 10,
'NSO_M' : 2,
'NTEMP_M' : 1,
'NODES_EX_M' : 3,

243
src/msspec/iodata.py
View File

@ -17,7 +17,7 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/iodata.py
# Last modified: Wed, 18 Jun 2025 13:30:09 +0200
# Last modified: Wed, 26 Feb 2025 11:10:17 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
@ -70,7 +70,8 @@ Here is an example of how to store values in a Data object:
import os
import sys
from looseversion import LooseVersion
from distutils.version import LooseVersion
from distutils.version import StrictVersion
from io import StringIO
from datetime import datetime
@ -79,9 +80,10 @@ from ase.io.extxyz import read_xyz, write_xyz
import h5py
import numpy as np
from lxml import etree
#from matplotlib.backends.backend_wxagg import FigureCanvasWx as FigureCanvas
from matplotlib.backends.backend_agg import FigureCanvasAgg
#from matplotlib.backends.backend_cairo import FigureCanvasCairo as FigureCanvasAgg
from matplotlib.figure import Figure
from matplotlib import pyplot as plt
from terminaltables import AsciiTable
import msspec
@ -334,15 +336,6 @@ class DataSet(object):
except:
pass
def get_views(self):
"""Returns all the defined views in the dataset.
:return: A list of view
:rtype: List of :py:class:`iodata._DataSetView`
"""
return self._views
@property
def views(self):
"""Returns all the defined views in the dataset.
@ -372,13 +365,7 @@ class DataSet(object):
mydset.add_parameter(name='Spectrometer', group='misc', value='Omicron', unit='')
"""
group = kwargs.get('group')
name = kwargs.get('name')
r = self.get_parameter(group=group, name=name)
if r:
r.update(**kwargs)
else:
self._parameters.append(kwargs)
self._parameters.append(kwargs)
def parameters(self):
"""
@ -411,39 +398,30 @@ class DataSet(object):
p.append(_)
return p[0] if len(p) == 1 else p
def set_cluster(self, cluster):
clusbuf = StringIO()
cluster.info['absorber'] = cluster.absorber
write_xyz(clusbuf, cluster)
self.add_parameter(group='Cluster', name='cluster', value=clusbuf.getvalue(), hidden="True")
def get_cluster(self):
"""Get all the atoms in the cluster.
:return: The cluster
:rtype: :py:class:`ase.Atoms`
"""
try:
p = self.get_parameter(group='Cluster', name='cluster')['value']
s = StringIO()
s.write(self.get_parameter(group='Cluster', name='cluster')['value'])
s.seek(0)
#return ase.io.read(s, format='xyz')
cluster = list(read_xyz(s))[-1]
return cluster
except:
return None
p = self.get_parameter(group='Cluster', name='cluster')['value']
s = StringIO()
s.write(self.get_parameter(group='Cluster', name='cluster')['value'])
s.seek(0)
#return ase.io.read(s, format='xyz')
cluster = list(read_xyz(s))[-1]
return cluster
def select(self, *args, **kwargs):
condition = kwargs.get('where', 'True')
indices = []
def export_views(self, folder, dpi=100):
for view in self.get_views():
def export_views(self, folder):
for view in self.views():
f = view.get_figure()
fname = os.path.join(folder, view.title) + '.png'
f.savefig(fname, dpi=dpi)
f.savefig(fname)
def export(self, filename="", mode="w"):
@ -510,7 +488,9 @@ class DataSet(object):
if isinstance(value, t):
fmt = f
break
#fd.write(' ')
fd.write(fmt.format(value))
#fd.write(str(value) + ', ')
fd.write('\n')
def __getitem__(self, itemspec):
@ -555,6 +535,7 @@ class DataSet(object):
def __len__(self):
try:
#length = len(self._col_arrays[0])
length = 0
for array in self._col_arrays:
length = max(length, len(array))
@ -690,7 +671,6 @@ class Data(object):
return
else:
data_grp = fd.create_group('DATA')
data_grp.attrs['dset_names'] = titles
meta_grp = fd.create_group('MsSpec viewer metainfo')
data_grp.attrs['title'] = self.title
@ -701,7 +681,6 @@ class Data(object):
continue
grp = data_grp.create_group(dset.title)
grp.attrs['notes'] = dset.notes
grp.attrs['col_names'] = dset.columns()
for col_name in dset.columns():
data = dset[col_name]
grp.create_dataset(col_name, data=data)
@ -712,7 +691,7 @@ class Data(object):
# xmlize views
for dset in self._datasets:
views_node = etree.SubElement(root, 'views', dataset=dset.title)
for view in dset.get_views():
for view in dset.views():
view_el = etree.fromstring(view.to_xml())
views_node.append(view_el)
@ -733,7 +712,7 @@ class Data(object):
self._dirty = False
LOGGER.info('Data saved in {}'.format(os.path.abspath(filename)))
def export(self, folder, overwrite=False, dpi=150):
def export(self, folder, overwrite=False):
os.makedirs(folder, exist_ok=overwrite)
for dset in self._datasets:
dset_name = dset.title.replace(' ', '_')
@ -741,7 +720,7 @@ class Data(object):
os.makedirs(p, exist_ok=overwrite)
fname = os.path.join(p, dset_name) + '.txt'
dset.export(fname)
dset.export_views(p, dpi=dpi)
dset.export_views(p)
@staticmethod
def load(filename):
@ -758,20 +737,12 @@ class Data(object):
views = {}
output.title = fd['DATA'].attrs['title']
try:
dset_names = fd['DATA'].attrs['dset_names']
except:
dset_names = [_ for _ in fd['DATA']]
for dset_name in dset_names:
for dset_name in fd['DATA'] :
parameters[dset_name] = []
views[dset_name] = []
dset = output.add_dset(dset_name)
dset.notes = fd['DATA'][dset_name].attrs['notes']
try:
col_names = fd['DATA'][dset_name].attrs['col_names']
except:
col_names = [_ for _ in fd['DATA'][dset_name]]
for h5dset in col_names:
for h5dset in fd['DATA'][dset_name]:
dset.add_columns(**{h5dset: fd['DATA'][dset_name][h5dset][...]})
try:
@ -899,19 +870,10 @@ class _DataSetView(object):
data.append(values)
return data
def plot(self):
f = self.get_figure(backend='plt')
return f, f.get_axes()[0]
def get_figure(self, backend=None):
def get_figure(self):
opts = self._plotopts
if backend is None:
figure = Figure()
else:
figure = plt.figure(num="[{}][{}]".format(self.dataset.title, self.title))
figure = Figure()
axes = None
proj = opts['projection']
scale = opts['scale']
@ -947,6 +909,7 @@ class _DataSetView(object):
axes.set_yticklabels(theta_ticks)
cbar = figure.colorbar(im)
#im.set_clim(0, 0.0275)
elif proj == 'polar':
values[0] = np.radians(values[0])
@ -996,6 +959,7 @@ class _DataSetView(object):
root = etree.Element('view', name=self.title)
for key, value in list(plotopts.items()):
root.attrib[key] = str(value)
#root.attrib['dataset_name'] = self.dataset.title
for tags, cond, legend in zip(self._selection_tags,
self._selection_conditions,
@ -1013,12 +977,21 @@ class _DataSetView(object):
def from_xml(self, xmlstr):
root = etree.fromstring(xmlstr)
self.title = root.attrib['name']
#self._plotopts['title'] = root.attrib['title']
#self._plotopts['xlabel'] = root.attrib['xlabel']
# self._plotopts['ylabel'] = root.attrib['ylabel']
# self._plotopts['grid'] = bool(root.attrib['grid'])
# self._plotopts['colorbar'] = bool(root.attrib['colorbar'])
# self._plotopts['projection'] = root.attrib['projection']
# self._plotopts['marker'] = root.attrib['marker']
for key in list(self._plotopts.keys()):
try:
self._plotopts[key] = eval(root.attrib.get(key))
except:
self._plotopts[key] = root.attrib.get(key)
legends = []
conditions = []
tags = []
@ -1078,6 +1051,8 @@ if has_gui:
if is_hidden == "True":
continue
group = datatree.get(p['group'], [])
#strval = str(p['value'] * p['unit'] if p['unit'] else p['value'])
#group.append("{:s} = {:s}".format(p['name'], strval))
group.append("{} = {} {}".format(p['name'], p['value'], p['unit']))
datatree[p['group']] = group
@ -1103,7 +1078,7 @@ if has_gui:
self._filename = None
self._current_dset = None
wx.Frame.__init__(self, None, title="", size=(800, 600))
wx.Frame.__init__(self, None, title="", size=(640, 480))
self.Bind(wx.EVT_CLOSE, self.on_close)
@ -1119,6 +1094,7 @@ if has_gui:
# Add the notebook to hold all graphs
self.notebooks = {}
sizer = wx.BoxSizer(wx.VERTICAL)
#sizer.Add(self.notebook)
self.SetSizer(sizer)
self.Bind(wx.EVT_NOTEBOOK_PAGE_CHANGED, self.on_page_changed)
@ -1135,8 +1111,9 @@ if has_gui:
for dset in self.data:
nb = wx.Notebook(self, -1)
self.notebooks[dset.title] = nb
#self.GetSizer().Add(nb, 1, wx.ALL|wx.EXPAND)
self.GetSizer().Add(nb, proportion=1, flag=wx.ALL|wx.EXPAND)
for view in dset.get_views():
for view in dset.views():
self.create_page(nb, view)
self.create_menu()
@ -1270,6 +1247,13 @@ if has_gui:
cluster_viewer = ClusterViewer(win, size=wx.Size(480, 340))
dset = self.data[self._current_dset]
#s = StringIO()
#s.write(dset.get_parameter(group='Cluster', name='cluster')['value'])
#_s = dset.get_parameter(group='Cluster', name='cluster')['value']
#print(_s)
# rewind to the begining of the string
#s.seek(0)
#atoms = ase.io.read(s, format='xyz')
atoms = dset.get_cluster()
cluster_viewer.set_atoms(atoms, rescale=True, center=True)
cluster_viewer.rotate_atoms(0., 180.)
@ -1307,10 +1291,6 @@ if has_gui:
self.Layout()
self.update_statusbar()
self._current_dset = name
has_cluster = True if self.data[self._current_dset].get_cluster() is not None else False
menu_item = self.GetMenuBar().FindItemById(302)
menu_item.Enable(has_cluster)
def create_page(self, nb, view):
# Get the matplotlib figure
@ -1344,6 +1324,95 @@ if has_gui:
nb.AddPage(p, view.title)
canvas.draw()
def OLDcreate_page(self, nb, view):
opts = view._plotopts
p = wx.Panel(nb, -1)
figure = Figure()
axes = None
proj = opts['projection']
scale = opts['scale']
if proj == 'rectilinear':
axes = figure.add_subplot(111, projection='rectilinear')
elif proj in ('polar', 'ortho', 'stereo'):
axes = figure.add_subplot(111, projection='polar')
canvas = FigureCanvas(p, -1, figure)
sizer = wx.BoxSizer(wx.VERTICAL)
toolbar = NavigationToolbar2WxAgg(canvas)
toolbar.Realize()
sizer.Add(toolbar, 0, wx.ALL|wx.EXPAND)
toolbar.update()
sizer.Add(canvas, 5, wx.ALL|wx.EXPAND)
p.SetSizer(sizer)
p.Fit()
p.Show()
for values, label in zip(view.get_data(), opts['legend']):
# if we have only one column to plot, select a bar graph
if np.shape(values)[0] == 1:
xvalues = list(range(len(values[0])))
axes.bar(xvalues, values[0], label=label,
picker=5)
axes.set_xticks(xvalues)
else:
if proj in ('ortho', 'stereo'):
theta, phi, Xsec = cols2matrix(*values)
theta_ticks = np.arange(0, 91, 15)
if proj == 'ortho':
R = np.sin(np.radians(theta))
R_ticks = np.sin(np.radians(theta_ticks))
elif proj == 'stereo':
R = 2 * np.tan(np.radians(theta/2.))
R_ticks = 2 * np.tan(np.radians(theta_ticks/2.))
#R = np.tan(np.radians(theta/2.))
X, Y = np.meshgrid(np.radians(phi), R)
im = axes.pcolormesh(X, Y, Xsec)
axes.set_yticks(R_ticks)
axes.set_yticklabels(theta_ticks)
figure.colorbar(im)
elif proj == 'polar':
values[0] = np.radians(values[0])
axes.plot(*values, label=label, picker=5,
marker=opts['marker'])
else:
if scale == 'semilogx':
pltcmd = axes.semilogx
elif scale == 'semilogy':
pltcmd = axes.semilogy
elif scale == 'log':
pltcmd = axes.loglog
else:
pltcmd = axes.plot
pltcmd(*values, label=label, picker=5,
marker=opts['marker'])
axes.grid(opts['grid'])
axes.set_title(opts['title'])
axes.set_xlabel(opts['xlabel'])
axes.set_ylabel(opts['ylabel'])
axes.set_xlim(*opts['xlim'])
axes.set_ylim(*opts['ylim'])
if label:
axes.legend()
axes.autoscale(enable=opts['autoscale'])
# MPL events
figure.canvas.mpl_connect('motion_notify_event', self.on_mpl_motion)
figure.canvas.mpl_connect('pick_event', self.on_mpl_pick)
nb.AddPage(p, view.title)
def update_statusbar(self):
sb = self.GetStatusBar()
menu_id = self.GetMenuBar().FindMenu('Datasets')
@ -1377,7 +1446,41 @@ if has_gui:
if __name__ == "__main__":
if False:
data = Data('all my data')
dset = data.add_dset('Dataset 0')
X = np.arange(0, 20)
Y = X**2
dset.add_columns(x=X, y=Y, z=X+2, w=Y**3)
dset.add_parameter(name='truc', group='main', value='3.14', unit='eV')
dset.add_parameter(name='machin', group='main', value='abc', unit='')
# Z = [0,1]
#
# for z in Z:
# for x, y in zip(X, Y):
# dset.add_row(x=x, y=y, z=z, random=np.random.rand())
#
#
view = dset.add_view('my view', autoscale=True)
view.select('x', 'y', where="z<10", legend=r"z = 0")
view.select('x', 'y', where="z>10", legend=r"z = 1")
print(dset.get_parameter(group='main'))
constraint = lambda a, b: (a > 10 and a < 15) and b > 0
indices = list(map(constraint, dset.x, dset.w))
print(dset.y[indices])
#data.view()
import sys
data = Data.load(sys.argv[1])
data.view()

2
src/msspec/looper.py
View File

@ -17,7 +17,7 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/looper.py
# Last modified: Thu, 27 Feb 2025 16:33:09 +0100
# Last modified: Wed, 26 Feb 2025 11:15:54 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
from collections import OrderedDict

28
src/msspec/parameters.py
View File

@ -19,8 +19,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/parameters.py
# Last modified: Mon, 16 Jun 2025 14:42:03 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
# Last modified: Tue, 15 Feb 2022 15:37:28 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>
"""
@ -362,8 +362,8 @@ class PhagenParameters(BaseParameters):
Parameter('l2h', types=(int,), fmt='d', default=4),
Parameter('ionicity', types=dict, default={}),
#Parameter('noproto', allowed_values=('.true.', '.false.'),
# types=(str,), fmt='>7s', default='.true.'),
Parameter('noproto', allowed_values=('.true.', '.false.'),
types=(str,), fmt='>7s', default='.true.'),
#Parameter('absorber', types=(int,), limits=(1, None), fmt='d', default=1),
#Parameter('nosym', types=(str,), allowed_values=('.true.', '.false.'), fmt='s', default='.true.'),
#Parameter('outersph', types=(str,), allowed_values=('.true.', '.false.'), fmt='s', default='.false.'),
@ -400,7 +400,7 @@ class SpecParameters(BaseParameters):
fmt='d'),
Parameter('calctype_ipol', types=int, limits=(-1, 2), default=0,
fmt='d'),
Parameter('calctype_iamp', types=int, limits=(0, 1), default=0,
Parameter('calctype_iamp', types=int, limits=(0, 1), default=1,
fmt='d'),
Parameter('ped_li', types=str, default='1s'),
@ -606,7 +606,7 @@ class SpecParameters(BaseParameters):
Parameter('eigval_beta', types=float, default=1., fmt='.2f'),
Parameter('calc_no', types=int, limits=[0, 8], default=1, fmt='d'),
Parameter('calc_ndif', types=int, limits=[1, 18], default=3,
Parameter('calc_ndif', types=int, limits=[1, 10], default=3,
fmt='d'),
Parameter('calc_ispher', types=int, limits=[0, 1], default=1,
fmt='d'),
@ -621,7 +621,7 @@ class SpecParameters(BaseParameters):
fmt='d'),
Parameter('calc_irdia', types=int, limits=[0, 1], default=0,
fmt='d'),
Parameter('calc_itrtl', types=int, limits=[0, 9], default=0,
Parameter('calc_itrtl', types=int, limits=[1, 9], default=7,
fmt='d'),
Parameter('calc_itest', types=int, limits=[0, 2], default=0,
fmt='d'),
@ -640,7 +640,7 @@ class SpecParameters(BaseParameters):
fmt='d'),
Parameter('calc_ira', types=int, limits=[0, 1], default=0, fmt='d'),
Parameter('calc_ipw', types=int, limits=[0, 1], default=0, fmt='d'),
Parameter('calc_ncut', types=int, limits=[0, 18], default=2,
Parameter('calc_ncut', types=int, limits=[0, 10], default=2,
fmt='d'),
Parameter('calc_pctint', types=float, limits=[1e-4, 999.9999],
default=0.01, fmt='.4f'),
@ -1441,19 +1441,19 @@ class CalculationParameters(BaseParameters):
It is only meaningful for the series expansion algorithm.
Its value is limited to 8 but it is rarely necessary to go beyond
2 or 3."""),
Parameter('scattering_order', types=int, limits=(1, 18), default=3,
Parameter('scattering_order', types=int, limits=(1, 10), default=3,
doc="""
The scattering order. Only meaningful for the 'expansion' algorithm.
Its value is limited to 10."""),
Parameter('renormalization_mode', allowed_values=(None, 'G_n', 'Sigma_n',
'Z_n', 'Pi_1', 'L_n'),
'Z_n', 'Pi_1', 'Lowdin'),
types=(type(None), str), default=None,
doc="""
Enable the calculation of the coefficients for the renormalization of
the multiple scattering series.
You can choose to renormalize in terms of the :math:`G_n`, the
:math:`\\Sigma_n`, the :math:`Z_n`, the Löwdin :math:`\\Pi_1` or
:math:`L_n` matrices"""),
:math:`\\Sigma_n`, the :math:`Z_n`, the :math:`\\Pi_1` or the Lowdin
:math:`K^2` matrices"""),
Parameter('renormalization_omega', types=(int,float,complex),
default=1.+0j,
doc="""
@ -1509,7 +1509,7 @@ class CalculationParameters(BaseParameters):
path is rejected and its contribution to the scattering path
operator won’t be computed.
"""),
Parameter('scattering_order_cutoff', types=int, limits=(0, 18),
Parameter('scattering_order_cutoff', types=int, limits=(0, 10),
default=2, doc="""
Used in conjunction with the ‘*plane_wave_normal*’ filter. It states
to activate the plane wave approximation (which is fast but
@ -1589,7 +1589,7 @@ class CalculationParameters(BaseParameters):
'Sigma_n': 2,
'Z_n' : 3,
'Pi_1' : 4,
'L_n' : 5}
'Lowdin' : 5}
# Check that the method is neither 'Z_n' nor 'K^2' for other
# 'spetroscopy' than EIG
try:

74
src/msspec/phagen/fortran/msxas3.inc
View File

@ -1,74 +0,0 @@
c.. dimensions for the program
integer ua_
parameter ( nat_ = 4000,
$ ua_ = 4000,
$ neq_ = 48,
$ thrs = -0.001d0 )
C
C where :
c
c nat_ maximum number of atoms expected in any
c molecule of interest (including an outer
c sphere. an even number is suggested).
c
c ua_ maximum number of nda's (unique, or
c symmetry-distinct atoms) expected in any
c molecule (including an outer sphere).
c
c neq_ maximum number of atoms expected in
c any symmetry-equivalent set (including
c the nda of the set)
c
c thrs threshold within which atoms with the same
c atomic number are considered equivalent.
c If negative, all atoms are considered prototypical.
c
c Warning: This version of msxas3.inc with program
c phagen_scf_2.3_dp.f
c
c...................................................................
c dimensioning cont and cont_sub source program
c...................................................................
c
integer fl_, rdx_
c
parameter ( rdx_ = 1600,
$ lmax_ = 80,
$ npss = lmax_ + 2,
$ fl_ = 2*npss + 1,
$ nef_ = 10,
$ lexp_ = 10,
$ nep_ = 500 )
c
c where :
c
c rdx_ number of points of the linear-log mesh
c
c lmax_ the maximum l-value used on any sphere
c (suggested value 5 or less if running valence dos section of
c phagen, 60 when calculating atomic t_l)
c
c nef_ effective number of atoms used in the transition
c matrix elements of eels. Put = 1 if not doing a eels
c calculation (suggested value 12)
c
c lexp_ lmax in the expansion of coulomb interaction plus one! temporary
c
c nep_ the maximum number of energy points for which phase
c shifts will be computed.
c
c.......................................................................
c multiple scattering paths, xn programs dimensioning
c.......................................................................
c
c
parameter (natoms=nat_)
c
c
c where:
c
c natoms = number of centers in the system
c
c
c...................................................................

1
src/msspec/phagen/fortran/msxas3.inc Symbolic link
View File

@ -0,0 +1 @@
phagen_2.2_dp/msxas3.inc

27
src/msspec/phagen/fortran/msxasc3.inc
View File

@ -1,27 +0,0 @@
logical vinput, nosym, tdl
character*5 potype
character*1 optrsh
character*2 edge,charelx,edge1,edge2,potgen,relc
character*3 calctype,expmode,eikappr,enunit
character*4 coor
character*6 norman
character*7 ionzst
integer absorber,hole,l2h,hole1,hole2
dimension nz(natoms)
dimension c(natoms,3), rad(natoms), redf(natoms)
dimension neqat(natoms)
dimension nk0(0:lmax_)
c.....Warning: when reordering common/options/, reorder also the same common in
c.....subroutine inpot
common/options/rsh,ovlpfac,vc0,rs0,vinput,absorber,hole,mode,
& ionzst,potype,norman,coor,charelx,edge,potgen,lmax_mode,
& lmaxt,relc,eikappr,optrsh,nosym,tdl
common/atoms/c,rad,redf,charge_ion(100),nat,nz,neqat
c common/azimuth/lin,lmax
common/auger/calctype,expmode,edge1,edge2
common/auger1/lin1,lin2,hole1,hole2,l2h
common/funit/idat,iwr,iphas,iedl0,iwf
common/constant/antoau,ev,pi,pi4,pif,zero,thresh,nk0
c....................................................................
c rpot = if real potential is to be used
c.....................................................................

1
src/msspec/phagen/fortran/msxasc3.inc Symbolic link
View File

@ -0,0 +1 @@
phagen_2.2_dp/msxasc3.inc

74
src/msspec/phagen/fortran/msxast3.inc
View File

@ -1,74 +0,0 @@
c.. dimensions for the program
integer ua_
parameter ( nat_ = 4000,
$ ua_ = 4000,
$ neq_ = 48,
$ thrs = -0.001d0 )
C
C where :
c
c nat_ maximum number of atoms expected in any
c molecule of interest (including an outer
c sphere. an even number is suggested).
c
c ua_ maximum number of nda's (unique, or
c symmetry-distinct atoms) expected in any
c molecule (including an outer sphere).
c
c neq_ maximum number of atoms expected in
c any symmetry-equivalent set (including
c the nda of the set)
c
c thrs threshold within which atoms with the same
c atomic number are considered equivalent.
c If negative, all atoms are considered prototypical.
c
c Warning: This version of msxas3.inc with program
c phagen_scf_2.3_dp.f
c
c...................................................................
c dimensioning cont and cont_sub source program
c...................................................................
c
integer fl_, rdx_
c
parameter ( rdx_ = 1600,
$ lmax_ = 80,
$ npss = lmax_ + 2,
$ fl_ = 2*npss + 1,
$ nef_ = 10,
$ lexp_ = 10,
$ nep_ = 500 )
c
c where :
c
c rdx_ number of points of the linear-log mesh
c
c lmax_ the maximum l-value used on any sphere
c (suggested value 5 or less if running valence dos section of
c phagen, 60 when calculating atomic t_l)
c
c nef_ effective number of atoms used in the transition
c matrix elements of eels. Put = 1 if not doing a eels
c calculation (suggested value 12)
c
c lexp_ lmax in the expansion of coulomb interaction plus one! temporary
c
c nep_ the maximum number of energy points for which phase
c shifts will be computed.
c
c.......................................................................
c multiple scattering paths, xn programs dimensioning
c.......................................................................
c
c
parameter (natoms=nat_)
c
c
c where:
c
c natoms = number of centers in the system
c
c
c...................................................................

82
src/msspec/phagen/fortran/phagen_2.3_dp/phagen_scf_2.3_dp.f
View File

@ -1,9 +1,5 @@
C
CST ==> Phagen to python shared object modifications
CST Phagen becomes a subroutine in a library
CST PROGRAM PHAGEN
subroutine phagen()
CST Phagen to python shared object modifications <==
PROGRAM PHAGEN
C
C ....................................
C .. ..
@ -152,14 +148,6 @@ C
C... Starting to write in the check file IWR
C
WRITE(IWR,1000)
CST ==> Phagen to python shared object modifications
CST Create output folders and open input file
CALL SYSTEM('mkdir -p div/wf')
CALL SYSTEM('mkdir -p plot')
CALL SYSTEM('mkdir -p tl')
CALL SYSTEM('mkdir -p clus')
OPEN(idat, FILE='../input/input.ms', STATUS='old')
CST Phagen to python shared object modifications <==
C
C... Opening the Fortran files
C
@ -313,10 +301,6 @@ C
CLOSE(46)
CLOSE(IWF)
CLOSE(IPHAS)
CST ==> Phagen to python shared object modifications <==
CST explicitely close fort.55
CLOSE(55)
CST Phagen to python shared object modifications <==
C
C Formats:
C
@ -345,7 +329,6 @@ C
complex*16 eelsme,p1,p2,p3,ramfsr1,ramfsr2,ramfsr3
complex*16 p3irreg,p2irreg
C
C
C.....................................................................
C
common /continuum/ emin,emax,delta,cip,gamma,eftri,iexcpot,db
@ -1106,11 +1089,7 @@ c
c.....Write out atomic coordinates in symmetry-program order:
c each prototypical atom is followed by its sym-equivalent atoms
c
CST ==> Phagen to python shared object modifications
CST The clus.out file needs to be opened, so I uncommented the
CSR line below
open (10,file='clus/clus.out',status='unknown')
CST Phagen to python shared object modifications <==
c open (10,file='clus/clus.out',status='unknown')
if( coor.eq.'au ') then
ipha=1
coef=1.d0
@ -1119,10 +1098,7 @@ CST Phagen to python shared object modifications <==
ipha=2
coef=0.529177d0
endif
CST ==> Phagen to python shared object modifications <==
CST I uncommented the line below
write(10,888) ipha
CST Phagen to python shared object modifications <==
c write(10,888) ipha
888 format(30x,i1)
write(7,10) (neqat(i),i=1,nat)
10 format (/,16i5,//)
@ -1140,21 +1116,14 @@ c
no = no + 1
write(7,20) no,nsymbl(k),nzeq(k),xv(k)-x0,
& yv(k)-y0,zv(k)-z0,neqat(k-1)
CST ==> Phagen to python shared object modifications
CST I changed the unit to 10
CST write(7,20) no,nsymbl(k),nzeq(k),(xv(k)-x0)*coef,
CST & (yv(k)-y0)*coef,(zv(k)-z0)*coef,neqat(k-1)
write(10,20) no,nsymbl(k),nzeq(k),(xv(k)-x0)*coef,
write(7,20) no,nsymbl(k),nzeq(k),(xv(k)-x0)*coef,
& (yv(k)-y0)*coef,(zv(k)-z0)*coef,neqat(k-1)
endif
continue
enddo
enddo
c
CST ==> Phagen to python shared object modifications
CST I uncommented the line below
close(10)
CST Phagen to python shared object modifications <==
c close(10)
c
20 format (i5,6x,a4,i5,3f10.4,i5)
c
@ -1698,10 +1667,9 @@ c
common/transform/trans
logical shift_cc
c
CST ==> Phagen to python shared object modifications
CST data zero,thrs/0.0d0,0.001d0/ !if thrs is negative, all cluster atoms are considered prototypical
data zero/0.0d0/
CST Phagen to python shared object modifications <==
c data zero,thrs/0.0d0,-0.001d0/ !if thrs is negative, all cluster atoms
c are considered prototypical
data zero/0.0d0/
c
data jtape/21/
data lunout/7/
@ -2169,10 +2137,7 @@ c dgc contains large component radial functions
c common /deux/ dvn(251), dvf(251), d(251), dc(251), dgc(251,30)
c passc and rho0 contain 4*pi*r^2*rho(r)
c
CST ==> Phagen to python shared object modifications
CST dimension r(mp),r_hs(440),rho0_hs(440)
dimension r(mp),r_hs(*),rho0_hs(*)
CST Phagen to python shared object modifications <==
dimension r(mp),r_hs(440),rho0_hs(440)
C
dimension dum1(mp), dum2(mp)
dimension vcoul(mp), rho0(mp), enp(ms)
@ -3360,11 +3325,6 @@ C
5004 FORMAT(8(I5,1X,F7.2))
WRITE(6,5003)
WRITE(6,*)
ELSE
WRITE(6,5003)
WRITE(6,*) ' External radii read in as: '
WRITE(6,*) ' i rs(i) i=1,natoms '
WRITE(6,5004) (I, RS(I), I=1,NATOMSM)
END IF
IF(NWR1.NE.' PCH') GO TO 999
WRITE(7,*)
@ -7588,14 +7548,7 @@ c.....this subroutine calculates the radial matrix elements
c.....necessary for eels cross-section
c.....using a linear-log mesh
c
CST ==> Phagen to Python shared object modifications
CST I replaced the line below
CST common/mtxele/ nstart,nlast
common/mtxele/ nstart,nlast,dmx(2),dmx1(2),qmx(3),qmx1(3),
$ dxdir,dxexc,nfis,nfis1,nfis2
real*8 nfis,nfis2,nfis1
complex*16 dmx,dmx1,qmx,qmx1,dxdir,dxexc
CST Phagen to Python shared object modifications <==
common/mtxele/ nstart,nlast
c
common/mtxelex/ dmxx(2),dmxx1(2),dmxxa(2),dmxxa1(2),
& qmxx(3),qmxx1(3),qmxxa(3),qmxxa1(3),
@ -16918,14 +16871,7 @@ c.....(i=1,2) for lfin=l0i-1 (i=1) and lfin=l0i+1 (i=2) both for
c.....the regular (dmxx) and irregular solution (dmxx1) using a
c.....linear-log mesh
c
CST ==> Phagen to Python shared object modifications
CST I replaced the line below
CST common/mtxele/ nstart,nlast
common/mtxele/ nstart,nlast,dmx(2),dmx1(2),qmx(3),qmx1(3),
$ dxdir,dxexc,nfis,nfis1,nfis2
real*8 nfis,nfis2,nfis1
complex*16 dmx,dmx1,qmx,qmx1,dxdir,dxexc
CST Phagen to Python shared object modifications <==
common/mtxele/ nstart,nlast
c
common/mtxelex/ dmxx(2),dmxx1(2),dmxxa(2),dmxxa1(2),
& qmxx(3),qmxx1(3),qmxxa(3),qmxxa1(3),
@ -20114,11 +20060,7 @@ C MARK 12 RELEASE. NAG COPYRIGHT 1986.
C MARK 15 REVISED. IER-915 (APR 1991).
C .. Scalar Arguments ..
INTEGER INFO
CST ==> Phagen to Python shared object modifications
CST I changed the dimension to automatic
CST CHARACTER*13 SRNAME
CHARACTER(*) SRNAME
CST Phagen to Python shared object modifications <==
CHARACTER*13 SRNAME
C ..
C
C Purpose

23866
src/msspec/phagen/fortran/phagen_scf.f
View File

File diff suppressed because it is too large Load Diff

1
src/msspec/phagen/fortran/phagen_scf.f Symbolic link
View File

@ -0,0 +1 @@
phagen_2.2_dp/phagen_scf_2.2_dp.f

17
src/msspec/spec/fortran/common_sub/read_data.f
View File

@ -83,9 +83,8 @@ C
REAL TEXTE1(10),TEXTE2(10),TEXTE3(10)
REAL TEXTE4(10),TEXTE5(10),TEXTE6(10)
REAL TEXTE6B(10),TEXTE7(10)
REAL THFWD(NATP_M),THBWD(NATP_M)
REAL THFWD(NATP_M),THBWD(NATP_M),GLG(0:N_GAUNT),NJ(0:N_GAUNT)
REAL ALPHAR,BETAR,RACC
REAL*8 GLG(0:N_GAUNT),NJ(0:N_GAUNT)
C
C
C
@ -922,9 +921,6 @@ C
READ(ICOM,34) OUTFILE2,IUO2
READ(ICOM,34) OUTFILE3,IUO3
READ(ICOM,34) OUTFILE4,IUO4
IF(.NOT.(OUTFILE1.EQ.'-')) THEN
OPEN(UNIT=IUO1, FILE=OUTFILE1, STATUS='UNKNOWN')
ENDIF
C
IUSCR=MAX0(ICOM,IUI2,IUI3,IUI4,IUI5,IUI6,IUI7,IUI8,IUI9,IUO1,IUO2,
&IUO3,IUO4)+1
@ -1408,12 +1404,11 @@ C
C
C Computing the renormalization coefficients
C
C IF(I_REN.LE.4) THEN
C CALL COEF_RENORM(NDIF)
C ELSEIF(I_REN.EQ.5) THEN
C CALL COEF_LOEWDIN(NDIF)
C ENDIF
CALL COEF_RENORM(NDIF)
IF(I_REN.LE.4) THEN
CALL COEF_RENORM(NDIF)
ELSEIF(I_REN.EQ.5) THEN
CALL COEF_LOEWDIN(NDIF)
ENDIF
C
C Storage of the logarithm of the Gamma function GLD(N+1,N_INT)
C for integer (N_INT=1) and semi-integer (N_INT=2) values :

1
src/msspec/spec/fortran/eig/common/plotfd.f
View File

@ -8,7 +8,6 @@ C
C INCLUDE 'spec.inc'
C
C
USE DIM_MOD
USE APPROX_MOD
USE FDIF_MOD
USE INIT_L_MOD, L => LI, I2 => INITL, I3 => NNL, I4 => LF1, I5 =>

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

@ -781,7 +781,7 @@ C=======================================================================
C=======================================================================
MODULE EXPFAC_MOD
IMPLICIT NONE
REAL*8, ALLOCATABLE, DIMENSION(:,:) :: EXPF
REAL, ALLOCATABLE, DIMENSION(:,:) :: EXPF
CONTAINS
SUBROUTINE ALLOC_EXPFAC()
USE DIM_MOD
@ -837,7 +837,7 @@ C=======================================================================
C=======================================================================
MODULE EXPFAC2_MOD
IMPLICIT NONE
REAL*8, ALLOCATABLE, DIMENSION(:,:) :: EXPF2
REAL, ALLOCATABLE, DIMENSION(:,:) :: EXPF2
CONTAINS
SUBROUTINE ALLOC_EXPFAC2()
USE DIM_MOD

2
src/msspec/spec/fortran/phd_mi_noso_nosp_nosym/phddif_mi.f
View File

@ -44,7 +44,7 @@ C
USE TYPCAL_MOD
USE TYPEM_MOD
USE TYPEXP_MOD
USE VALIN_MOD, PHLUM => PHILUM
USE VALIN_MOD
USE VALIN_AV_MOD
USE VALFIN_MOD
C

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

@ -1459,7 +1459,7 @@ CST 9 FORMAT(3X,F9.4,1X,F9.4,5X,E12.6,5X,E12.6)
95 FORMAT(////,31X,'AUGER LINE :',A6,//)
97 FORMAT(///,19X,'(PLANE WAVES MULTIPLE SCATTERING - ORDER ',I1,')')
&
98 FORMAT(///,17X,'(SPHERICAL WAVES MULTIPLE SCATTERING - ORDER ',I2,
98 FORMAT(///,17X,'(SPHERICAL WAVES MULTIPLE SCATTERING - ORDER ',I1,
&')')
100 FORMAT(///,8X,'<<<<<<<<<< WRONG NAME FOR THE INITIAL STATE',' >>
&>>>>>>>>')

4
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths1.f
View File

@ -327,10 +327,6 @@ C
IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN
CALL ARCSIN(COMPL1,CTROIS1,PHIJK)
ENDIF
C PRINT *,"KTYP=",KTYP,"(X,Y,Z)=",SYM_AT(1,KTYP),
C &SYM_AT(2,KTYP),SYM_AT(3,KTYP)
C PRINT *,"JTYP=",JTYP,"(X,Y,Z)=",SYM_AT(1,JTYP),
C &SYM_AT(2,JTYP),SYM_AT(3,JTYP)
CALL EULER(THJK,PHIJK,THIJ,PHIIJ,AIJK,BIJK,CIJK,IEULER)
IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2
CALL MATDIF(NO,ND-1,LF2,JTYP,KTYP,JE,I_ABS,ISPEED,ISPHER,A

4
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths5.f
View File

@ -342,8 +342,8 @@ C
ENDIF
ENDIF
IF(ND.EQ.NDIF) GOTO 32
CALL FINDPATHS6(ND,KTYP,KATL,I_CP,R,XR,YR,ZR,RHOJK,
1 THJK,PHIJK,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
c CALL FINDPATHS(ND,KTYP,KATL,I_CP,R,XR,YR,ZR,RHOJK,
c 1 THJK,PHIJK,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
32 DIJ=DIJ-R(ND)
22 IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1
IT(ND-1)=0

367
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths6.f
View File

@ -1,367 +0,0 @@
C
C=======================================================================
C
SUBROUTINE FINDPATHS6(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,
& PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
C
C This routine generates all the paths and filters them according to the
C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH).
C It corresponds to the spin-independent case from a non spin-orbit
C resolved initial core state LI
C
C Last modified : 16 May 2007
C
USE DIM_MOD
C
USE APPROX_MOD , ILE => ILENGTH, RLE => RLENGTH
USE COOR_MOD
USE DEBWAL_MOD
USE INIT_L_MOD
USE PATH_MOD
USE ROT_MOD
USE TESTPA_MOD
USE TESTPB_MOD
USE TRANS_MOD
USE TLDW_MOD
USE VARIA_MOD
C
DIMENSION XR(NDIF_M),YR(NDIF_M),ZR(NDIF_M)
DIMENSION JPOS(NDIF_M,3),R(NDIF_M)
C
C
C
COMPLEX PW1,PWI,FTHETA,RHOMI,RHOIJ,RHOJK
COMPLEX IC,COMPL1,PW(0:NDIF_M)
COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M)
COMPLEX YLM1(0:NL_M,-NL_M:NL_M)
COMPLEX YLM2(0:NL_M,-NL_M:NL_M),CTL,CTL2
C
DATA XCOMP,PI4,SMALL /1.E-10,12.566371,0.0001/
C
IC=(0.,1.)
IEULER=1
C
IF(IFWD.EQ.1) COSFWDI=COS(RTHFWD(ITYP))
IF(IBWD(ITYP).EQ.1) COSBWDI=COS(RTHBWD(ITYP))
C
C I_CP = 0 : all open paths generated
C I_CP = 1 : only closed paths generated
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN
N_TYP=N_PROT
ELSE
N_TYP=1
ENDIF
C
DO JTYP=1,N_TYP
IF(IFWD.EQ.1) COSFWDJ=COS(RTHFWD(JTYP))
IF(IBWD(JTYP).EQ.1) COSBWDJ=COS(RTHBWD(JTYP))
ND=ND+1
C
C I_ABS = 0 : the atom before the scatterer is not the absorber
C I_ABS = 1 : the atom before the scatterer is the absorber
C I_ABS = 2 : the atom after the scatterer is the absorber (XAS only)
C
IF(ND.EQ.1) THEN
I_ABS=1
ELSE
I_ABS=0
ENDIF
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN
NBTYPJ=NATYP(JTYP)
ELSE
NBTYPJ=1
ENDIF
C
DO JNUM=1,NBTYPJ
JATL=NCORR(JNUM,JTYP)
IF(JATL.EQ.IATL) GOTO 12
XR(ND)=SYM_AT(1,JATL)-SYM_AT(1,IATL)
YR(ND)=SYM_AT(2,JATL)-SYM_AT(2,IATL)
ZR(ND)=SYM_AT(3,JATL)-SYM_AT(3,IATL)
R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND))
DIJ=DIJ+R(ND)
IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1
IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42
JPOS(ND,1)=JTYP
JPOS(ND,2)=JNUM
JPOS(ND,3)=JATL
NPATH(ND)=NPATH(ND)+1.
IF(ND.GT.1) THEN
COSTHMIJ=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))/(
&R(ND)*R(ND-1))
IF(IFWD.EQ.1) THEN
IF(IBWD(ITYP).EQ.0) THEN
IF(COSTHMIJ.LT.COSFWDI) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ELSEIF(IBWD(ITYP).EQ.1) THEN
IF((COSTHMIJ.GT.COSBWDI).AND.(COSTHMIJ.LT.-SMALL)) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
IF((COSTHMIJ.LT.COSFWDI).AND.(COSTHMIJ.GE.0.)) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42
RHOIJ=VK(JE)*R(ND)
CTROIS1=ZR(ND)/R(ND)
IF(CTROIS1.GT.1.) THEN
CTROIS1=1.
ELSEIF(CTROIS1.LT.-1.) THEN
CTROIS1=-1.
ENDIF
THIJ=ACOS(CTROIS1)
COMPL1= XR(ND)+IC*YR(ND)
CALL ARCSIN(COMPL1,CTROIS1,PHIIJ)
IF((ND.GT.1).AND.((ND-1).LT.NDIF)) THEN
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 40
ZSURFI=ZSURF-ZR(ND-1)
IF(IDCM.EQ.1) THEN
UJ2(ITYP)=SIG2(R(ND-1),ITYP)
ENDIF
IF(ABS(ZSURFI).LE.SMALL) THEN
IF(ABS(COSTHMIJ-1.).GT.SMALL) THEN
CSKZ2I=(CTROIS1-COS(THMI))*(CTROIS1-COS(THMI))/(2.-2.*CO
&STHMIJ)
ELSE
CSKZ2I=1.
ENDIF
UII=UJ2(ITYP)*(1.+CSKZ2I*(RSJ-1.))
ELSE
UII=UJ2(ITYP)
ENDIF
IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN
XK2UI2=VK2(JE)*UII
CALL DWSPH(ITYP,JE,XK2UI2,TLT,ISPEED)
ENDIF
40 IF(IDWSPH.EQ.1) THEN
DW(ND-1)=1.
ELSE
DW(ND-1)=EXP(-VK2(JE)*UII*(1.-COSTHMIJ))
ENDIF
ENDIF
IF(ND.EQ.1) THEN
RHO01=RHOIJ
TH01=THIJ
PHI01=PHIIJ
CALL DJMN2(TH01,RLM01,LF2,2)
GOTO 30
ENDIF
IF(IPW.EQ.1) THEN
CALL FACDIF(COSTHMIJ,JPOS(ND-1,1),JE,FTHETA)
PWI=FTHETA*DW(ND-1)/R(ND)
PW(ND)=PW(ND-1)*PWI
CTL2=PI4*PW(ND)*CEX(1)/VK(JE)
LMJ=LMAX(ITYP,JE)
IF(ND.GT.NCUT) THEN
IT(ND)=1
ELSE
IT(ND)=0
ENDIF
CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2)
CALL HARSPH2(NL_M,THIJ,PHIIJ,YLM2,LMJ)
XMAXT=0.
DO LJ=0,LMJ
CTL=CTL2*TL(LJ,1,JTYP,JE)*YLM2(LJ,0)
DO LF=LF1,LF2,ISTEP_LF
PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE)
XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO
ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND)=0
ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 42
CALL EULER(THIJ,PHIIJ,THMI,PHIMI,AMIJ,BMIJ,CMIJ,IEULER)
IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2
CALL MATDIF(NO,ND-1,LF2,ITYP,JTYP,JE,I_ABS,ISPEED,ISPHER,AMIJ,
&BMIJ,CMIJ,RHOMI,RHOIJ)
30 CEX(ND)=CEXP(IC*RHOIJ)/R(ND)
CEXDW(ND)=CEX(ND)*DW(ND-1)
IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN
IF((I_CP.EQ.0).OR.(JATL.EQ.1)) THEN
CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOIJ,THIJ,PHIIJ,F
&REF,IJ,DIJ,TAU)
NPATH2(ND)=NPATH2(ND)+1.
ENDIF
ENDIF
IF(ND.EQ.NDIF) GOTO 42
I_ABS=0
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN
N_TYP=N_PROT
ELSE
N_TYP=1
ENDIF
C
DO KTYP=1,N_TYP
ND=ND+1
IF(ND.GT.NDIF) GOTO 20
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN
NBTYPK=NATYP(KTYP)
ELSE
NBTYPK=1
ENDIF
C
DO KNUM=1,NBTYPK
KATL=NCORR(KNUM,KTYP)
IF(KATL.EQ.JATL) GOTO 22
JPOS(ND,1)=KTYP
JPOS(ND,2)=KNUM
JPOS(ND,3)=KATL
XR(ND)=SYM_AT(1,KATL)-SYM_AT(1,JATL)
YR(ND)=SYM_AT(2,KATL)-SYM_AT(2,JATL)
ZR(ND)=SYM_AT(3,KATL)-SYM_AT(3,JATL)
R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND))
DIJ=DIJ+R(ND)
IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1
IF(IT(ND-1).EQ.1) GOTO 32
RHOJK=R(ND)*VK(JE)
NPATH(ND)=NPATH(ND)+1.
COSTHIJK=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))
&/(R(ND)*R(ND-1))
IF(IFWD.EQ.1) THEN
IF(IBWD(JTYP).EQ.0) THEN
IF(COSTHIJK.LT.COSFWDJ) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ELSEIF(IBWD(JTYP).EQ.1) THEN
IF((COSTHIJK.GT.COSBWDJ).AND. (COSTHIJK.LT.-SMALL)) TH
&EN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
IF((COSTHIJK.LT.COSFWDJ).AND.(COSTHIJK.GE.0.))THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ENDIF
ENDIF
IF(IT(ND-1).EQ.1) GOTO 32
CTROIS1=ZR(ND)/R(ND)
IF(CTROIS1.GT.1) THEN
CTROIS1=1.
ELSEIF(CTROIS1.LT.-1.) THEN
CTROIS1=-1.
ENDIF
THJK=ACOS(CTROIS1)
COMPL1= XR(ND)+IC*YR(ND)
IF(ND-1.LT.NDIF) THEN
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 50
ZSURFJ=ZSURF-ZR(ND-1)
IF(IDCM.EQ.1) THEN
UJ2(JTYP)=SIG2(R(ND-1),JTYP)
ENDIF
IF(ABS(ZSURFJ).LE.SMALL) THEN
IF(ABS(COSTHIJK-1.).GT.SMALL) THEN
CSKZ2J=(CTROIS1-COS(THIJ))*(CTROIS1-COS(THIJ))/(2.-2
&.*COSTHIJK)
ELSE
CSKZ2J=1.
ENDIF
UJJ=UJ2(JTYP)*(1.+CSKZ2J*(RSJ-1.))
ELSE
UJJ=UJ2(JTYP)
ENDIF
IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN
XK2UJ2=VK2(JE)*UJJ
CALL DWSPH(JTYP,JE,XK2UJ2,TLT,ISPEED)
ENDIF
50 IF(IDWSPH.EQ.1) THEN
DW(ND-1)=1.
ELSE
DW(ND-1)=EXP(-VK2(JE)*UJJ*(1.-COSTHIJK))
ENDIF
ENDIF
IF(IPW.EQ.1) THEN
CALL FACDIF(COSTHIJK,JPOS(ND-1,1),JE,FTHETA)
PWI=FTHETA*DW(ND-1)/R(ND)
PW(ND)=PW(ND-1)*PWI
CTL2=PI4*PW(ND)*CEX(1)/VK(JE)
CALL ARCSIN(COMPL1,CTROIS1,PHIJK)
LMJ=LMAX(JTYP,JE)
IF(ND.GT.NCUT) THEN
IT(ND)=1
ELSE
IT(ND)=0
ENDIF
CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2)
CALL HARSPH2(NL_M,THJK,PHIJK,YLM2,LMJ)
XMAXT=0.
DO LJ=0,LMJ
CTL=CTL2*TL(LJ,1,KTYP,JE)*YLM2(LJ,0)
DO LF=LF1,LF2,ISTEP_LF
PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE)
XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO
ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))
& IT(ND)=0
ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32
IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN
CALL ARCSIN(COMPL1,CTROIS1,PHIJK)
ENDIF
CALL EULER(THJK,PHIJK,THIJ,PHIIJ,AIJK,BIJK,CIJK,IEULER)
IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2
CALL MATDIF(NO,ND-1,LF2,JTYP,KTYP,JE,I_ABS,ISPEED,ISPHER,A
&IJK,BIJK,CIJK,RHOIJ,RHOJK)
CEX(ND)=CEXP(IC*RHOJK)/R(ND)
CEXDW(ND)=CEX(ND)*DW(ND-1)
IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN
IF((I_CP.EQ.0).OR.(KATL.EQ.1)) THEN
CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOJK,THJK,PHI
&JK,FREF,IJ,DIJ,TAU)
NPATH2(ND)=NPATH2(ND)+1.
ENDIF
ENDIF
IF(ND.EQ.NDIF) GOTO 32
CALL FINDPATHS7(ND,KTYP,KATL,I_CP,R,XR,YR,ZR,RHOJK,
1 THJK,PHIJK,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
32 DIJ=DIJ-R(ND)
22 IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1
IT(ND-1)=0
IN(ND-1)=0
ENDDO
20 CONTINUE
ND=ND-1
ENDDO
42 DIJ=DIJ-R(ND)
12 IF(ND.GT.1) THEN
IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1
IT(ND-1)=0
IN(ND-1)=0
ENDIF
ENDDO
ND=ND-1
ENDDO
C
RETURN
C
END

367
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths7.f
View File

@ -1,367 +0,0 @@
C
C=======================================================================
C
SUBROUTINE FINDPATHS7(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,
& PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
C
C This routine generates all the paths and filters them according to the
C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH).
C It corresponds to the spin-independent case from a non spin-orbit
C resolved initial core state LI
C
C Last modified : 16 May 2007
C
USE DIM_MOD
C
USE APPROX_MOD , ILE => ILENGTH, RLE => RLENGTH
USE COOR_MOD
USE DEBWAL_MOD
USE INIT_L_MOD
USE PATH_MOD
USE ROT_MOD
USE TESTPA_MOD
USE TESTPB_MOD
USE TRANS_MOD
USE TLDW_MOD
USE VARIA_MOD
C
DIMENSION XR(NDIF_M),YR(NDIF_M),ZR(NDIF_M)
DIMENSION JPOS(NDIF_M,3),R(NDIF_M)
C
C
C
COMPLEX PW1,PWI,FTHETA,RHOMI,RHOIJ,RHOJK
COMPLEX IC,COMPL1,PW(0:NDIF_M)
COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M)
COMPLEX YLM1(0:NL_M,-NL_M:NL_M)
COMPLEX YLM2(0:NL_M,-NL_M:NL_M),CTL,CTL2
C
DATA XCOMP,PI4,SMALL /1.E-10,12.566371,0.0001/
C
IC=(0.,1.)
IEULER=1
C
IF(IFWD.EQ.1) COSFWDI=COS(RTHFWD(ITYP))
IF(IBWD(ITYP).EQ.1) COSBWDI=COS(RTHBWD(ITYP))
C
C I_CP = 0 : all open paths generated
C I_CP = 1 : only closed paths generated
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN
N_TYP=N_PROT
ELSE
N_TYP=1
ENDIF
C
DO JTYP=1,N_TYP
IF(IFWD.EQ.1) COSFWDJ=COS(RTHFWD(JTYP))
IF(IBWD(JTYP).EQ.1) COSBWDJ=COS(RTHBWD(JTYP))
ND=ND+1
C
C I_ABS = 0 : the atom before the scatterer is not the absorber
C I_ABS = 1 : the atom before the scatterer is the absorber
C I_ABS = 2 : the atom after the scatterer is the absorber (XAS only)
C
IF(ND.EQ.1) THEN
I_ABS=1
ELSE
I_ABS=0
ENDIF
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN
NBTYPJ=NATYP(JTYP)
ELSE
NBTYPJ=1
ENDIF
C
DO JNUM=1,NBTYPJ
JATL=NCORR(JNUM,JTYP)
IF(JATL.EQ.IATL) GOTO 12
XR(ND)=SYM_AT(1,JATL)-SYM_AT(1,IATL)
YR(ND)=SYM_AT(2,JATL)-SYM_AT(2,IATL)
ZR(ND)=SYM_AT(3,JATL)-SYM_AT(3,IATL)
R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND))
DIJ=DIJ+R(ND)
IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1
IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42
JPOS(ND,1)=JTYP
JPOS(ND,2)=JNUM
JPOS(ND,3)=JATL
NPATH(ND)=NPATH(ND)+1.
IF(ND.GT.1) THEN
COSTHMIJ=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))/(
&R(ND)*R(ND-1))
IF(IFWD.EQ.1) THEN
IF(IBWD(ITYP).EQ.0) THEN
IF(COSTHMIJ.LT.COSFWDI) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ELSEIF(IBWD(ITYP).EQ.1) THEN
IF((COSTHMIJ.GT.COSBWDI).AND.(COSTHMIJ.LT.-SMALL)) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
IF((COSTHMIJ.LT.COSFWDI).AND.(COSTHMIJ.GE.0.)) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42
RHOIJ=VK(JE)*R(ND)
CTROIS1=ZR(ND)/R(ND)
IF(CTROIS1.GT.1.) THEN
CTROIS1=1.
ELSEIF(CTROIS1.LT.-1.) THEN
CTROIS1=-1.
ENDIF
THIJ=ACOS(CTROIS1)
COMPL1= XR(ND)+IC*YR(ND)
CALL ARCSIN(COMPL1,CTROIS1,PHIIJ)
IF((ND.GT.1).AND.((ND-1).LT.NDIF)) THEN
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 40
ZSURFI=ZSURF-ZR(ND-1)
IF(IDCM.EQ.1) THEN
UJ2(ITYP)=SIG2(R(ND-1),ITYP)
ENDIF
IF(ABS(ZSURFI).LE.SMALL) THEN
IF(ABS(COSTHMIJ-1.).GT.SMALL) THEN
CSKZ2I=(CTROIS1-COS(THMI))*(CTROIS1-COS(THMI))/(2.-2.*CO
&STHMIJ)
ELSE
CSKZ2I=1.
ENDIF
UII=UJ2(ITYP)*(1.+CSKZ2I*(RSJ-1.))
ELSE
UII=UJ2(ITYP)
ENDIF
IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN
XK2UI2=VK2(JE)*UII
CALL DWSPH(ITYP,JE,XK2UI2,TLT,ISPEED)
ENDIF
40 IF(IDWSPH.EQ.1) THEN
DW(ND-1)=1.
ELSE
DW(ND-1)=EXP(-VK2(JE)*UII*(1.-COSTHMIJ))
ENDIF
ENDIF
IF(ND.EQ.1) THEN
RHO01=RHOIJ
TH01=THIJ
PHI01=PHIIJ
CALL DJMN2(TH01,RLM01,LF2,2)
GOTO 30
ENDIF
IF(IPW.EQ.1) THEN
CALL FACDIF(COSTHMIJ,JPOS(ND-1,1),JE,FTHETA)
PWI=FTHETA*DW(ND-1)/R(ND)
PW(ND)=PW(ND-1)*PWI
CTL2=PI4*PW(ND)*CEX(1)/VK(JE)
LMJ=LMAX(ITYP,JE)
IF(ND.GT.NCUT) THEN
IT(ND)=1
ELSE
IT(ND)=0
ENDIF
CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2)
CALL HARSPH2(NL_M,THIJ,PHIIJ,YLM2,LMJ)
XMAXT=0.
DO LJ=0,LMJ
CTL=CTL2*TL(LJ,1,JTYP,JE)*YLM2(LJ,0)
DO LF=LF1,LF2,ISTEP_LF
PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE)
XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO
ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND)=0
ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 42
CALL EULER(THIJ,PHIIJ,THMI,PHIMI,AMIJ,BMIJ,CMIJ,IEULER)
IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2
CALL MATDIF(NO,ND-1,LF2,ITYP,JTYP,JE,I_ABS,ISPEED,ISPHER,AMIJ,
&BMIJ,CMIJ,RHOMI,RHOIJ)
30 CEX(ND)=CEXP(IC*RHOIJ)/R(ND)
CEXDW(ND)=CEX(ND)*DW(ND-1)
IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN
IF((I_CP.EQ.0).OR.(JATL.EQ.1)) THEN
CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOIJ,THIJ,PHIIJ,F
&REF,IJ,DIJ,TAU)
NPATH2(ND)=NPATH2(ND)+1.
ENDIF
ENDIF
IF(ND.EQ.NDIF) GOTO 42
I_ABS=0
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN
N_TYP=N_PROT
ELSE
N_TYP=1
ENDIF
C
DO KTYP=1,N_TYP
ND=ND+1
IF(ND.GT.NDIF) GOTO 20
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN
NBTYPK=NATYP(KTYP)
ELSE
NBTYPK=1
ENDIF
C
DO KNUM=1,NBTYPK
KATL=NCORR(KNUM,KTYP)
IF(KATL.EQ.JATL) GOTO 22
JPOS(ND,1)=KTYP
JPOS(ND,2)=KNUM
JPOS(ND,3)=KATL
XR(ND)=SYM_AT(1,KATL)-SYM_AT(1,JATL)
YR(ND)=SYM_AT(2,KATL)-SYM_AT(2,JATL)
ZR(ND)=SYM_AT(3,KATL)-SYM_AT(3,JATL)
R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND))
DIJ=DIJ+R(ND)
IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1
IF(IT(ND-1).EQ.1) GOTO 32
RHOJK=R(ND)*VK(JE)
NPATH(ND)=NPATH(ND)+1.
COSTHIJK=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))
&/(R(ND)*R(ND-1))
IF(IFWD.EQ.1) THEN
IF(IBWD(JTYP).EQ.0) THEN
IF(COSTHIJK.LT.COSFWDJ) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ELSEIF(IBWD(JTYP).EQ.1) THEN
IF((COSTHIJK.GT.COSBWDJ).AND. (COSTHIJK.LT.-SMALL)) TH
&EN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
IF((COSTHIJK.LT.COSFWDJ).AND.(COSTHIJK.GE.0.))THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ENDIF
ENDIF
IF(IT(ND-1).EQ.1) GOTO 32
CTROIS1=ZR(ND)/R(ND)
IF(CTROIS1.GT.1) THEN
CTROIS1=1.
ELSEIF(CTROIS1.LT.-1.) THEN
CTROIS1=-1.
ENDIF
THJK=ACOS(CTROIS1)
COMPL1= XR(ND)+IC*YR(ND)
IF(ND-1.LT.NDIF) THEN
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 50
ZSURFJ=ZSURF-ZR(ND-1)
IF(IDCM.EQ.1) THEN
UJ2(JTYP)=SIG2(R(ND-1),JTYP)
ENDIF
IF(ABS(ZSURFJ).LE.SMALL) THEN
IF(ABS(COSTHIJK-1.).GT.SMALL) THEN
CSKZ2J=(CTROIS1-COS(THIJ))*(CTROIS1-COS(THIJ))/(2.-2
&.*COSTHIJK)
ELSE
CSKZ2J=1.
ENDIF
UJJ=UJ2(JTYP)*(1.+CSKZ2J*(RSJ-1.))
ELSE
UJJ=UJ2(JTYP)
ENDIF
IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN
XK2UJ2=VK2(JE)*UJJ
CALL DWSPH(JTYP,JE,XK2UJ2,TLT,ISPEED)
ENDIF
50 IF(IDWSPH.EQ.1) THEN
DW(ND-1)=1.
ELSE
DW(ND-1)=EXP(-VK2(JE)*UJJ*(1.-COSTHIJK))
ENDIF
ENDIF
IF(IPW.EQ.1) THEN
CALL FACDIF(COSTHIJK,JPOS(ND-1,1),JE,FTHETA)
PWI=FTHETA*DW(ND-1)/R(ND)
PW(ND)=PW(ND-1)*PWI
CTL2=PI4*PW(ND)*CEX(1)/VK(JE)
CALL ARCSIN(COMPL1,CTROIS1,PHIJK)
LMJ=LMAX(JTYP,JE)
IF(ND.GT.NCUT) THEN
IT(ND)=1
ELSE
IT(ND)=0
ENDIF
CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2)
CALL HARSPH2(NL_M,THJK,PHIJK,YLM2,LMJ)
XMAXT=0.
DO LJ=0,LMJ
CTL=CTL2*TL(LJ,1,KTYP,JE)*YLM2(LJ,0)
DO LF=LF1,LF2,ISTEP_LF
PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE)
XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO
ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))
& IT(ND)=0
ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32
IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN
CALL ARCSIN(COMPL1,CTROIS1,PHIJK)
ENDIF
CALL EULER(THJK,PHIJK,THIJ,PHIIJ,AIJK,BIJK,CIJK,IEULER)
IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2
CALL MATDIF(NO,ND-1,LF2,JTYP,KTYP,JE,I_ABS,ISPEED,ISPHER,A
&IJK,BIJK,CIJK,RHOIJ,RHOJK)
CEX(ND)=CEXP(IC*RHOJK)/R(ND)
CEXDW(ND)=CEX(ND)*DW(ND-1)
IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN
IF((I_CP.EQ.0).OR.(KATL.EQ.1)) THEN
CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOJK,THJK,PHI
&JK,FREF,IJ,DIJ,TAU)
NPATH2(ND)=NPATH2(ND)+1.
ENDIF
ENDIF
IF(ND.EQ.NDIF) GOTO 32
CALL FINDPATHS8(ND,KTYP,KATL,I_CP,R,XR,YR,ZR,RHOJK,
1 THJK,PHIJK,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
32 DIJ=DIJ-R(ND)
22 IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1
IT(ND-1)=0
IN(ND-1)=0
ENDDO
20 CONTINUE
ND=ND-1
ENDDO
42 DIJ=DIJ-R(ND)
12 IF(ND.GT.1) THEN
IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1
IT(ND-1)=0
IN(ND-1)=0
ENDIF
ENDDO
ND=ND-1
ENDDO
C
RETURN
C
END

367
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths8.f
View File

@ -1,367 +0,0 @@
C
C=======================================================================
C
SUBROUTINE FINDPATHS8(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,
& PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
C
C This routine generates all the paths and filters them according to the
C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH).
C It corresponds to the spin-independent case from a non spin-orbit
C resolved initial core state LI
C
C Last modified : 16 May 2007
C
USE DIM_MOD
C
USE APPROX_MOD , ILE => ILENGTH, RLE => RLENGTH
USE COOR_MOD
USE DEBWAL_MOD
USE INIT_L_MOD
USE PATH_MOD
USE ROT_MOD
USE TESTPA_MOD
USE TESTPB_MOD
USE TRANS_MOD
USE TLDW_MOD
USE VARIA_MOD
C
DIMENSION XR(NDIF_M),YR(NDIF_M),ZR(NDIF_M)
DIMENSION JPOS(NDIF_M,3),R(NDIF_M)
C
C
C
COMPLEX PW1,PWI,FTHETA,RHOMI,RHOIJ,RHOJK
COMPLEX IC,COMPL1,PW(0:NDIF_M)
COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M)
COMPLEX YLM1(0:NL_M,-NL_M:NL_M)
COMPLEX YLM2(0:NL_M,-NL_M:NL_M),CTL,CTL2
C
DATA XCOMP,PI4,SMALL /1.E-10,12.566371,0.0001/
C
IC=(0.,1.)
IEULER=1
C
IF(IFWD.EQ.1) COSFWDI=COS(RTHFWD(ITYP))
IF(IBWD(ITYP).EQ.1) COSBWDI=COS(RTHBWD(ITYP))
C
C I_CP = 0 : all open paths generated
C I_CP = 1 : only closed paths generated
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN
N_TYP=N_PROT
ELSE
N_TYP=1
ENDIF
C
DO JTYP=1,N_TYP
IF(IFWD.EQ.1) COSFWDJ=COS(RTHFWD(JTYP))
IF(IBWD(JTYP).EQ.1) COSBWDJ=COS(RTHBWD(JTYP))
ND=ND+1
C
C I_ABS = 0 : the atom before the scatterer is not the absorber
C I_ABS = 1 : the atom before the scatterer is the absorber
C I_ABS = 2 : the atom after the scatterer is the absorber (XAS only)
C
IF(ND.EQ.1) THEN
I_ABS=1
ELSE
I_ABS=0
ENDIF
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN
NBTYPJ=NATYP(JTYP)
ELSE
NBTYPJ=1
ENDIF
C
DO JNUM=1,NBTYPJ
JATL=NCORR(JNUM,JTYP)
IF(JATL.EQ.IATL) GOTO 12
XR(ND)=SYM_AT(1,JATL)-SYM_AT(1,IATL)
YR(ND)=SYM_AT(2,JATL)-SYM_AT(2,IATL)
ZR(ND)=SYM_AT(3,JATL)-SYM_AT(3,IATL)
R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND))
DIJ=DIJ+R(ND)
IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1
IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42
JPOS(ND,1)=JTYP
JPOS(ND,2)=JNUM
JPOS(ND,3)=JATL
NPATH(ND)=NPATH(ND)+1.
IF(ND.GT.1) THEN
COSTHMIJ=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))/(
&R(ND)*R(ND-1))
IF(IFWD.EQ.1) THEN
IF(IBWD(ITYP).EQ.0) THEN
IF(COSTHMIJ.LT.COSFWDI) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ELSEIF(IBWD(ITYP).EQ.1) THEN
IF((COSTHMIJ.GT.COSBWDI).AND.(COSTHMIJ.LT.-SMALL)) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
IF((COSTHMIJ.LT.COSFWDI).AND.(COSTHMIJ.GE.0.)) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42
RHOIJ=VK(JE)*R(ND)
CTROIS1=ZR(ND)/R(ND)
IF(CTROIS1.GT.1.) THEN
CTROIS1=1.
ELSEIF(CTROIS1.LT.-1.) THEN
CTROIS1=-1.
ENDIF
THIJ=ACOS(CTROIS1)
COMPL1= XR(ND)+IC*YR(ND)
CALL ARCSIN(COMPL1,CTROIS1,PHIIJ)
IF((ND.GT.1).AND.((ND-1).LT.NDIF)) THEN
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 40
ZSURFI=ZSURF-ZR(ND-1)
IF(IDCM.EQ.1) THEN
UJ2(ITYP)=SIG2(R(ND-1),ITYP)
ENDIF
IF(ABS(ZSURFI).LE.SMALL) THEN
IF(ABS(COSTHMIJ-1.).GT.SMALL) THEN
CSKZ2I=(CTROIS1-COS(THMI))*(CTROIS1-COS(THMI))/(2.-2.*CO
&STHMIJ)
ELSE
CSKZ2I=1.
ENDIF
UII=UJ2(ITYP)*(1.+CSKZ2I*(RSJ-1.))
ELSE
UII=UJ2(ITYP)
ENDIF
IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN
XK2UI2=VK2(JE)*UII
CALL DWSPH(ITYP,JE,XK2UI2,TLT,ISPEED)
ENDIF
40 IF(IDWSPH.EQ.1) THEN
DW(ND-1)=1.
ELSE
DW(ND-1)=EXP(-VK2(JE)*UII*(1.-COSTHMIJ))
ENDIF
ENDIF
IF(ND.EQ.1) THEN
RHO01=RHOIJ
TH01=THIJ
PHI01=PHIIJ
CALL DJMN2(TH01,RLM01,LF2,2)
GOTO 30
ENDIF
IF(IPW.EQ.1) THEN
CALL FACDIF(COSTHMIJ,JPOS(ND-1,1),JE,FTHETA)
PWI=FTHETA*DW(ND-1)/R(ND)
PW(ND)=PW(ND-1)*PWI
CTL2=PI4*PW(ND)*CEX(1)/VK(JE)
LMJ=LMAX(ITYP,JE)
IF(ND.GT.NCUT) THEN
IT(ND)=1
ELSE
IT(ND)=0
ENDIF
CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2)
CALL HARSPH2(NL_M,THIJ,PHIIJ,YLM2,LMJ)
XMAXT=0.
DO LJ=0,LMJ
CTL=CTL2*TL(LJ,1,JTYP,JE)*YLM2(LJ,0)
DO LF=LF1,LF2,ISTEP_LF
PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE)
XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO
ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND)=0
ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 42
CALL EULER(THIJ,PHIIJ,THMI,PHIMI,AMIJ,BMIJ,CMIJ,IEULER)
IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2
CALL MATDIF(NO,ND-1,LF2,ITYP,JTYP,JE,I_ABS,ISPEED,ISPHER,AMIJ,
&BMIJ,CMIJ,RHOMI,RHOIJ)
30 CEX(ND)=CEXP(IC*RHOIJ)/R(ND)
CEXDW(ND)=CEX(ND)*DW(ND-1)
IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN
IF((I_CP.EQ.0).OR.(JATL.EQ.1)) THEN
CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOIJ,THIJ,PHIIJ,F
&REF,IJ,DIJ,TAU)
NPATH2(ND)=NPATH2(ND)+1.
ENDIF
ENDIF
IF(ND.EQ.NDIF) GOTO 42
I_ABS=0
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN
N_TYP=N_PROT
ELSE
N_TYP=1
ENDIF
C
DO KTYP=1,N_TYP
ND=ND+1
IF(ND.GT.NDIF) GOTO 20
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN
NBTYPK=NATYP(KTYP)
ELSE
NBTYPK=1
ENDIF
C
DO KNUM=1,NBTYPK
KATL=NCORR(KNUM,KTYP)
IF(KATL.EQ.JATL) GOTO 22
JPOS(ND,1)=KTYP
JPOS(ND,2)=KNUM
JPOS(ND,3)=KATL
XR(ND)=SYM_AT(1,KATL)-SYM_AT(1,JATL)
YR(ND)=SYM_AT(2,KATL)-SYM_AT(2,JATL)
ZR(ND)=SYM_AT(3,KATL)-SYM_AT(3,JATL)
R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND))
DIJ=DIJ+R(ND)
IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1
IF(IT(ND-1).EQ.1) GOTO 32
RHOJK=R(ND)*VK(JE)
NPATH(ND)=NPATH(ND)+1.
COSTHIJK=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))
&/(R(ND)*R(ND-1))
IF(IFWD.EQ.1) THEN
IF(IBWD(JTYP).EQ.0) THEN
IF(COSTHIJK.LT.COSFWDJ) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ELSEIF(IBWD(JTYP).EQ.1) THEN
IF((COSTHIJK.GT.COSBWDJ).AND. (COSTHIJK.LT.-SMALL)) TH
&EN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
IF((COSTHIJK.LT.COSFWDJ).AND.(COSTHIJK.GE.0.))THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ENDIF
ENDIF
IF(IT(ND-1).EQ.1) GOTO 32
CTROIS1=ZR(ND)/R(ND)
IF(CTROIS1.GT.1) THEN
CTROIS1=1.
ELSEIF(CTROIS1.LT.-1.) THEN
CTROIS1=-1.
ENDIF
THJK=ACOS(CTROIS1)
COMPL1= XR(ND)+IC*YR(ND)
IF(ND-1.LT.NDIF) THEN
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 50
ZSURFJ=ZSURF-ZR(ND-1)
IF(IDCM.EQ.1) THEN
UJ2(JTYP)=SIG2(R(ND-1),JTYP)
ENDIF
IF(ABS(ZSURFJ).LE.SMALL) THEN
IF(ABS(COSTHIJK-1.).GT.SMALL) THEN
CSKZ2J=(CTROIS1-COS(THIJ))*(CTROIS1-COS(THIJ))/(2.-2
&.*COSTHIJK)
ELSE
CSKZ2J=1.
ENDIF
UJJ=UJ2(JTYP)*(1.+CSKZ2J*(RSJ-1.))
ELSE
UJJ=UJ2(JTYP)
ENDIF
IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN
XK2UJ2=VK2(JE)*UJJ
CALL DWSPH(JTYP,JE,XK2UJ2,TLT,ISPEED)
ENDIF
50 IF(IDWSPH.EQ.1) THEN
DW(ND-1)=1.
ELSE
DW(ND-1)=EXP(-VK2(JE)*UJJ*(1.-COSTHIJK))
ENDIF
ENDIF
IF(IPW.EQ.1) THEN
CALL FACDIF(COSTHIJK,JPOS(ND-1,1),JE,FTHETA)
PWI=FTHETA*DW(ND-1)/R(ND)
PW(ND)=PW(ND-1)*PWI
CTL2=PI4*PW(ND)*CEX(1)/VK(JE)
CALL ARCSIN(COMPL1,CTROIS1,PHIJK)
LMJ=LMAX(JTYP,JE)
IF(ND.GT.NCUT) THEN
IT(ND)=1
ELSE
IT(ND)=0
ENDIF
CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2)
CALL HARSPH2(NL_M,THJK,PHIJK,YLM2,LMJ)
XMAXT=0.
DO LJ=0,LMJ
CTL=CTL2*TL(LJ,1,KTYP,JE)*YLM2(LJ,0)
DO LF=LF1,LF2,ISTEP_LF
PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE)
XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO
ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))
& IT(ND)=0
ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32
IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN
CALL ARCSIN(COMPL1,CTROIS1,PHIJK)
ENDIF
CALL EULER(THJK,PHIJK,THIJ,PHIIJ,AIJK,BIJK,CIJK,IEULER)
IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2
CALL MATDIF(NO,ND-1,LF2,JTYP,KTYP,JE,I_ABS,ISPEED,ISPHER,A
&IJK,BIJK,CIJK,RHOIJ,RHOJK)
CEX(ND)=CEXP(IC*RHOJK)/R(ND)
CEXDW(ND)=CEX(ND)*DW(ND-1)
IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN
IF((I_CP.EQ.0).OR.(KATL.EQ.1)) THEN
CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOJK,THJK,PHI
&JK,FREF,IJ,DIJ,TAU)
NPATH2(ND)=NPATH2(ND)+1.
ENDIF
ENDIF
IF(ND.EQ.NDIF) GOTO 32
CALL FINDPATHS9(ND,KTYP,KATL,I_CP,R,XR,YR,ZR,RHOJK,
1 THJK,PHIJK,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
32 DIJ=DIJ-R(ND)
22 IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1
IT(ND-1)=0
IN(ND-1)=0
ENDDO
20 CONTINUE
ND=ND-1
ENDDO
42 DIJ=DIJ-R(ND)
12 IF(ND.GT.1) THEN
IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1
IT(ND-1)=0
IN(ND-1)=0
ENDIF
ENDDO
ND=ND-1
ENDDO
C
RETURN
C
END

367
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths9.f
View File

@ -1,367 +0,0 @@
C
C=======================================================================
C
SUBROUTINE FINDPATHS9(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,
& PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
C
C This routine generates all the paths and filters them according to the
C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH).
C It corresponds to the spin-independent case from a non spin-orbit
C resolved initial core state LI
C
C Last modified : 16 May 2007
C
USE DIM_MOD
C
USE APPROX_MOD , ILE => ILENGTH, RLE => RLENGTH
USE COOR_MOD
USE DEBWAL_MOD
USE INIT_L_MOD
USE PATH_MOD
USE ROT_MOD
USE TESTPA_MOD
USE TESTPB_MOD
USE TRANS_MOD
USE TLDW_MOD
USE VARIA_MOD
C
DIMENSION XR(NDIF_M),YR(NDIF_M),ZR(NDIF_M)
DIMENSION JPOS(NDIF_M,3),R(NDIF_M)
C
C
C
COMPLEX PW1,PWI,FTHETA,RHOMI,RHOIJ,RHOJK
COMPLEX IC,COMPL1,PW(0:NDIF_M)
COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M)
COMPLEX YLM1(0:NL_M,-NL_M:NL_M)
COMPLEX YLM2(0:NL_M,-NL_M:NL_M),CTL,CTL2
C
DATA XCOMP,PI4,SMALL /1.E-10,12.566371,0.0001/
C
IC=(0.,1.)
IEULER=1
C
IF(IFWD.EQ.1) COSFWDI=COS(RTHFWD(ITYP))
IF(IBWD(ITYP).EQ.1) COSBWDI=COS(RTHBWD(ITYP))
C
C I_CP = 0 : all open paths generated
C I_CP = 1 : only closed paths generated
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN
N_TYP=N_PROT
ELSE
N_TYP=1
ENDIF
C
DO JTYP=1,N_TYP
IF(IFWD.EQ.1) COSFWDJ=COS(RTHFWD(JTYP))
IF(IBWD(JTYP).EQ.1) COSBWDJ=COS(RTHBWD(JTYP))
ND=ND+1
C
C I_ABS = 0 : the atom before the scatterer is not the absorber
C I_ABS = 1 : the atom before the scatterer is the absorber
C I_ABS = 2 : the atom after the scatterer is the absorber (XAS only)
C
IF(ND.EQ.1) THEN
I_ABS=1
ELSE
I_ABS=0
ENDIF
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN
NBTYPJ=NATYP(JTYP)
ELSE
NBTYPJ=1
ENDIF
C
DO JNUM=1,NBTYPJ
JATL=NCORR(JNUM,JTYP)
IF(JATL.EQ.IATL) GOTO 12
XR(ND)=SYM_AT(1,JATL)-SYM_AT(1,IATL)
YR(ND)=SYM_AT(2,JATL)-SYM_AT(2,IATL)
ZR(ND)=SYM_AT(3,JATL)-SYM_AT(3,IATL)
R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND))
DIJ=DIJ+R(ND)
IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1
IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42
JPOS(ND,1)=JTYP
JPOS(ND,2)=JNUM
JPOS(ND,3)=JATL
NPATH(ND)=NPATH(ND)+1.
IF(ND.GT.1) THEN
COSTHMIJ=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))/(
&R(ND)*R(ND-1))
IF(IFWD.EQ.1) THEN
IF(IBWD(ITYP).EQ.0) THEN
IF(COSTHMIJ.LT.COSFWDI) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ELSEIF(IBWD(ITYP).EQ.1) THEN
IF((COSTHMIJ.GT.COSBWDI).AND.(COSTHMIJ.LT.-SMALL)) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
IF((COSTHMIJ.LT.COSFWDI).AND.(COSTHMIJ.GE.0.)) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42
RHOIJ=VK(JE)*R(ND)
CTROIS1=ZR(ND)/R(ND)
IF(CTROIS1.GT.1.) THEN
CTROIS1=1.
ELSEIF(CTROIS1.LT.-1.) THEN
CTROIS1=-1.
ENDIF
THIJ=ACOS(CTROIS1)
COMPL1= XR(ND)+IC*YR(ND)
CALL ARCSIN(COMPL1,CTROIS1,PHIIJ)
IF((ND.GT.1).AND.((ND-1).LT.NDIF)) THEN
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 40
ZSURFI=ZSURF-ZR(ND-1)
IF(IDCM.EQ.1) THEN
UJ2(ITYP)=SIG2(R(ND-1),ITYP)
ENDIF
IF(ABS(ZSURFI).LE.SMALL) THEN
IF(ABS(COSTHMIJ-1.).GT.SMALL) THEN
CSKZ2I=(CTROIS1-COS(THMI))*(CTROIS1-COS(THMI))/(2.-2.*CO
&STHMIJ)
ELSE
CSKZ2I=1.
ENDIF
UII=UJ2(ITYP)*(1.+CSKZ2I*(RSJ-1.))
ELSE
UII=UJ2(ITYP)
ENDIF
IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN
XK2UI2=VK2(JE)*UII
CALL DWSPH(ITYP,JE,XK2UI2,TLT,ISPEED)
ENDIF
40 IF(IDWSPH.EQ.1) THEN
DW(ND-1)=1.
ELSE
DW(ND-1)=EXP(-VK2(JE)*UII*(1.-COSTHMIJ))
ENDIF
ENDIF
IF(ND.EQ.1) THEN
RHO01=RHOIJ
TH01=THIJ
PHI01=PHIIJ
CALL DJMN2(TH01,RLM01,LF2,2)
GOTO 30
ENDIF
IF(IPW.EQ.1) THEN
CALL FACDIF(COSTHMIJ,JPOS(ND-1,1),JE,FTHETA)
PWI=FTHETA*DW(ND-1)/R(ND)
PW(ND)=PW(ND-1)*PWI
CTL2=PI4*PW(ND)*CEX(1)/VK(JE)
LMJ=LMAX(ITYP,JE)
IF(ND.GT.NCUT) THEN
IT(ND)=1
ELSE
IT(ND)=0
ENDIF
CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2)
CALL HARSPH2(NL_M,THIJ,PHIIJ,YLM2,LMJ)
XMAXT=0.
DO LJ=0,LMJ
CTL=CTL2*TL(LJ,1,JTYP,JE)*YLM2(LJ,0)
DO LF=LF1,LF2,ISTEP_LF
PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE)
XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO
ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND)=0
ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 42
CALL EULER(THIJ,PHIIJ,THMI,PHIMI,AMIJ,BMIJ,CMIJ,IEULER)
IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2
CALL MATDIF(NO,ND-1,LF2,ITYP,JTYP,JE,I_ABS,ISPEED,ISPHER,AMIJ,
&BMIJ,CMIJ,RHOMI,RHOIJ)
30 CEX(ND)=CEXP(IC*RHOIJ)/R(ND)
CEXDW(ND)=CEX(ND)*DW(ND-1)
IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN
IF((I_CP.EQ.0).OR.(JATL.EQ.1)) THEN
CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOIJ,THIJ,PHIIJ,F
&REF,IJ,DIJ,TAU)
NPATH2(ND)=NPATH2(ND)+1.
ENDIF
ENDIF
IF(ND.EQ.NDIF) GOTO 42
I_ABS=0
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN
N_TYP=N_PROT
ELSE
N_TYP=1
ENDIF
C
DO KTYP=1,N_TYP
ND=ND+1
IF(ND.GT.NDIF) GOTO 20
C
IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN
NBTYPK=NATYP(KTYP)
ELSE
NBTYPK=1
ENDIF
C
DO KNUM=1,NBTYPK
KATL=NCORR(KNUM,KTYP)
IF(KATL.EQ.JATL) GOTO 22
JPOS(ND,1)=KTYP
JPOS(ND,2)=KNUM
JPOS(ND,3)=KATL
XR(ND)=SYM_AT(1,KATL)-SYM_AT(1,JATL)
YR(ND)=SYM_AT(2,KATL)-SYM_AT(2,JATL)
ZR(ND)=SYM_AT(3,KATL)-SYM_AT(3,JATL)
R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND))
DIJ=DIJ+R(ND)
IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1
IF(IT(ND-1).EQ.1) GOTO 32
RHOJK=R(ND)*VK(JE)
NPATH(ND)=NPATH(ND)+1.
COSTHIJK=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))
&/(R(ND)*R(ND-1))
IF(IFWD.EQ.1) THEN
IF(IBWD(JTYP).EQ.0) THEN
IF(COSTHIJK.LT.COSFWDJ) THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ELSEIF(IBWD(JTYP).EQ.1) THEN
IF((COSTHIJK.GT.COSBWDJ).AND. (COSTHIJK.LT.-SMALL)) TH
&EN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
IF((COSTHIJK.LT.COSFWDJ).AND.(COSTHIJK.GE.0.))THEN
NTHOF=NTHOF+1
IN(ND-1)=1
IF(NTHOF.GT.NTHOUT) THEN
IT(ND-1)=1
ENDIF
ENDIF
ENDIF
ENDIF
IF(IT(ND-1).EQ.1) GOTO 32
CTROIS1=ZR(ND)/R(ND)
IF(CTROIS1.GT.1) THEN
CTROIS1=1.
ELSEIF(CTROIS1.LT.-1.) THEN
CTROIS1=-1.
ENDIF
THJK=ACOS(CTROIS1)
COMPL1= XR(ND)+IC*YR(ND)
IF(ND-1.LT.NDIF) THEN
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 50
ZSURFJ=ZSURF-ZR(ND-1)
IF(IDCM.EQ.1) THEN
UJ2(JTYP)=SIG2(R(ND-1),JTYP)
ENDIF
IF(ABS(ZSURFJ).LE.SMALL) THEN
IF(ABS(COSTHIJK-1.).GT.SMALL) THEN
CSKZ2J=(CTROIS1-COS(THIJ))*(CTROIS1-COS(THIJ))/(2.-2
&.*COSTHIJK)
ELSE
CSKZ2J=1.
ENDIF
UJJ=UJ2(JTYP)*(1.+CSKZ2J*(RSJ-1.))
ELSE
UJJ=UJ2(JTYP)
ENDIF
IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN
XK2UJ2=VK2(JE)*UJJ
CALL DWSPH(JTYP,JE,XK2UJ2,TLT,ISPEED)
ENDIF
50 IF(IDWSPH.EQ.1) THEN
DW(ND-1)=1.
ELSE
DW(ND-1)=EXP(-VK2(JE)*UJJ*(1.-COSTHIJK))
ENDIF
ENDIF
IF(IPW.EQ.1) THEN
CALL FACDIF(COSTHIJK,JPOS(ND-1,1),JE,FTHETA)
PWI=FTHETA*DW(ND-1)/R(ND)
PW(ND)=PW(ND-1)*PWI
CTL2=PI4*PW(ND)*CEX(1)/VK(JE)
CALL ARCSIN(COMPL1,CTROIS1,PHIJK)
LMJ=LMAX(JTYP,JE)
IF(ND.GT.NCUT) THEN
IT(ND)=1
ELSE
IT(ND)=0
ENDIF
CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2)
CALL HARSPH2(NL_M,THJK,PHIJK,YLM2,LMJ)
XMAXT=0.
DO LJ=0,LMJ
CTL=CTL2*TL(LJ,1,KTYP,JE)*YLM2(LJ,0)
DO LF=LF1,LF2,ISTEP_LF
PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE)
XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO
ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))
& IT(ND)=0
ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32
IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN
CALL ARCSIN(COMPL1,CTROIS1,PHIJK)
ENDIF
CALL EULER(THJK,PHIJK,THIJ,PHIIJ,AIJK,BIJK,CIJK,IEULER)
IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2
CALL MATDIF(NO,ND-1,LF2,JTYP,KTYP,JE,I_ABS,ISPEED,ISPHER,A
&IJK,BIJK,CIJK,RHOIJ,RHOJK)
CEX(ND)=CEXP(IC*RHOJK)/R(ND)
CEXDW(ND)=CEX(ND)*DW(ND-1)
IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN
IF((I_CP.EQ.0).OR.(KATL.EQ.1)) THEN
CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOJK,THJK,PHI
&JK,FREF,IJ,DIJ,TAU)
NPATH2(ND)=NPATH2(ND)+1.
ENDIF
ENDIF
IF(ND.EQ.NDIF) GOTO 32
c CALL FINDPATHS(ND,KTYP,KATL,I_CP,R,XR,YR,ZR,RHOJK,
c 1 THJK,PHIJK,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
32 DIJ=DIJ-R(ND)
22 IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1
IT(ND-1)=0
IN(ND-1)=0
ENDDO
20 CONTINUE
ND=ND-1
ENDDO
42 DIJ=DIJ-R(ND)
12 IF(ND.GT.1) THEN
IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1
IT(ND-1)=0
IN(ND-1)=0
ENDIF
ENDDO
ND=ND-1
ENDDO
C
RETURN
C
END

323
src/msspec/spec/fortran/renormalization/renormalization.f
View File

@ -1,66 +1,53 @@
SUBROUTINE COEF_RENORM(NDIF)
!
! This subroutine computes the coefficients for the renormalization
! of the multiple scattering series. These coefficients are
! expressed as C_REN(K) where K is the multiple scattering order.
! REN2 is the value of the mixing (or renormalization) parameter.
!
! NDIF is the scattering order at which the series is truncated,
! so that K varies from 0 to NDIF.
!
! COMMON /RENORM/:
!
! I_REN = 1 : renormalization in terms of G_n matrices (n : N_REN)
! = 2 : renormalization in terms of the Sigma_n matrices
! = 3 : renormalization in terms of the Z_n matrices
! = 4 : Löwdin renormalization in terms of the Pi_1 matrices
! = 5 : Löwdin renormalization in terms of the L_n matrices
!
! N_REN = renormalization order n
!
! REN = REN_R + i * REN_I : omega
!
!
! Reference: A. Takatsu, S. Tricot, P. Schieffer, K. Dunseath,
! M. Terao-Dunseath, K. Hatada and D. Sébilleau,
! Phys. Chem. Chem. Phys., 2022, 24, 5658
!
!
! Authors : D. Sébilleau, A. Takatsu, M. Terao-Dunseath, K. Dunseath
!
!
! Last modified (DS,ST): 26 Feb 2025
!
C
C This subroutine computes the coefficients for the renormalization
C of the multiple scattering series. These coefficients are
C expressed as C_REN(K) where K is the multiple scattering order.
C REN2 is the value of the mixing (or renormalization) parameter.
C
C NDIF is the scattering order at which the series is truncated,
C so that K varies from 0 to NDIF.
C
C COMMON /RENORM/:
C
C I_REN = 1 : renormalization in terms of G_n matrices (n : N_REN)
C = 2 : renormalization in terms of the Sigma_n matrices
C = 3 : renormalization in terms of the Z_n matrices
C = 4 : renormalization in terms of the Pi_1 matrix
C
C N_REN = n
C
C REN = REN_R+IC*REN_I : omega
C
C Last modified : 11 Apr 2019
C
USE DIM_MOD
USE C_RENORM_MOD
USE RENORM_MOD
!
INTEGER M_MIN, M_MAX
!
C
REAL C(0:NDIF_M,0:NDIF_M)
!
COMPLEX X(0:NDIF_M,0:NDIF_M),SUM_L,POWER
C
COMPLEX REN,REN2,COEF1,COEF2,ZEROC,ONEC,IC
COMPLEX Y1(0:NDIF_M,0:NDIF_M)
!
!
C
C
ZEROC=(0.,0.)
ONEC=(1.,0.)
IC=(0.,1.)
!
C
REN=REN_R+IC*REN_I ! omega
!
! Initialisation of renormalization coefficients
!
C
C Initialisation of renormalization coefficients
C
DO J=0,NDIF
C_REN(J)=ZEROC
ENDDO
!
! Computing the binomial coefficients C(N,K) = (N) = N! / K! (N-K)!
! (K)
!CCC 2019.06.09 Aika
C
C Computing the binomial coefficients C(N,K) = (N) = N! / K! (N-K)!
C (K)
CCCC 2019.06.09 Aika
c=0.0
!CCC 2019.06.09 Aika
CCCC 2019.06.09 Aika
C(0,0)=1.
C(1,0)=1.
C(1,1)=1.
@ -71,121 +58,153 @@
C(N,K)=C(N-1,K)+C(N-1,K-1)
ENDDO
ENDDO
!
C
IF(I_REN.LE.3) THEN
!
! Computing the modified renormalization parameter REN2 (g_n,s_n,zeta_n)
!
C
C Computing the modified renormalization parameter REN2 (g_n,s_n,zeta_n)
C
IF(I_REN.EQ.1) THEN
!
!.....(g_n,G_n) renormalization
!
C
C.....(g_n,G_n) renormalization
C
REN2=REN**N_REN ! g_n = omega^n
!
C
ELSEIF(I_REN.EQ.2) THEN
!
!.....(s_{n},Sigma_n) renormalization
!
C
C.....(s_{n},Sigma_n) renormalization
C
REN2=(ONEC-REN**(N_REN+1))/(FLOAT(N_REN+1)*(ONEC-REN)) ! s_n
!
C
ELSEIF(I_REN.EQ.3) THEN
!
!.....(zeta_{n},Z_n) renormalization
!
! 2019.04.29
! REN2=(REN-ONEC)**(N_REN+1) ! zeta_n
! 2019.06.09
! REN2=-(REN-ONEC)**(N_REN+1) ! zeta_n
C
C.....(zeta_{n},Z_n) renormalization
C
C 2019.04.29
C REN2=(REN-ONEC)**(N_REN+1) ! zeta_n
C 2019.06.09
C REN2=-(REN-ONEC)**(N_REN+1) ! zeta_n
REN2=-(ONCE-REN)**(N_REN+1) ! zeta_n
!
C
ENDIF
!
! AT & MTD 2019.04.17 - summation over j ?
C
C AT & MTD 2019.04.17 - summation over j ?
DO K=0,NDIF
C_REN(K)=ZEROC
c_ren(k)=zeroc
DO J=K,NDIF
C_REN(K)=C_REN(K)+C(J,K)*(ONEC-REN2)**(J-K)
C_REN(K)=C_REN(K)+c(j,k)*(ONEC-REN2)**(J-K)
ENDDO
C_REN(K)=C_REN(K)*REN2**(K+1)
c_ren(k)=c_ren(k)*ren2**(k+1)
ENDDO
!
! DO K=0,NDIF
! COEF1=REN2**(K+1)
! DO J=K,NDIF
! COEF2=(ONEC-REN2)**(J-K)
! C_REN(J)=C_REN(J)+COEF1*COEF2*C(J,K)
! ENDDO
! ENDDO
!
C
C DO K=0,NDIF
C COEF1=REN2**(K+1)
C DO J=K,NDIF
C COEF2=(ONEC-REN2)**(J-K)
C C_REN(J)=C_REN(J)+COEF1*COEF2*C(J,K)
C ENDDO
C ENDDO
C
ELSEIF(I_REN.EQ.4) THEN
!
! Loewdin (Pi_1) renormalization for n = 1
!
! Notation: Y1(K,M) : [Y_1^k]_m
!
! with k : scattering order
! m : summation index
!
COEF1 = ONEC - REN ! (1 - omega)
!
Y1(0,0) = ONEC !
!
DO K = 1, NDIF !
M_MAX = MIN(K,NDIF) !
M_MIN = INT(K / 2) !
DO M = M_MIN, M_MAX !
COEF2 = (REN**(K-M)) * (COEF1**(2*M-K)) !
Y1(K,M) = COEF2 * ( C(M,K-M) + COEF1 * C(M,K-M-1) ) !
END DO !
END DO
! !
C_REN(0) = ONEC !
C_REN(1) = ONEC !
DO K = 2, NDIF !
C_REN(K) = ZEROC !
DO M = M_MIN, M_MAX !
C_REN(K) = C_REN(K) + Y1(M,K) !
END DO !
END DO !
ELSE IF(I_REN.EQ.5) THEN
!
! Loewdin L_n(omega,NDIF) renormalization
!
! Notation: X(K,N) = X_n(omega,k)
!
!
! Computing the X(N,K) coefficients, with K <= N
!
POWER = ONEC / REN !
DO N = 0, NDIF !
POWER = POWER * REN ! omega^n
IF(N == 0) THEN !
X(N,0) = ONEC !
ELSE !
X(N,0) = ZEROC !
END IF !
DO K = 1, NDIF !
IF(K > N) THEN !
X(N,K) = ZEROC !
ELSE IF(K == N) THEN !
X(N,K) = POWER * X(N-1,K-1) !
ELSE !
X(N,K) = X(N-1,K) * (REN - POWER) + POWER * X(N-1,K-1)!
END IF !
END DO !
END DO !
!
! Calculation of L_n(omega,NDIF)
!
DO N = 0, NDIF !
SUM_L = ZEROC !
DO K = N, NDIF !
SUM_L = SUM_L + X(K,N) !
END DO !
C_REN(N) = SUM_L !
END DO !
!
END IF !
!
C
C Loewdin (Pi_1) renormalization for n = 1
C
C Notation: Y1(M,K) : [Y_1^m]_k
C 2019.06.09
C Notation: Y1(K,M) : [Y_1^k]_m
C
COEF1=ONEC-REN ! (1 - omega)
DO K=0,NDIF
Y1(K,K)=COEF1**K
IF(K.LT.NDIF) THEN
DO M=K+1,NDIF
COEF2=(REN**(M-K))*(COEF1**(2*K-M))
C 2019.04.19 AT & MTD
C Y1(M,K)=COEF2*(C(K,M-K)+COEF1*C(K,M-K-1))
Y1(K,M)=COEF2*(C(K,M-K)+COEF1*C(K,M-K-1))
ENDDO
ENDIF
ENDDO
C
DO K=0,NDIF
IN=INT(K/2)
C_REN(K)=ZEROC
DO M=IN,K
C_REN(K)=C_REN(K)+Y1(M,K)
ENDDO
ENDDO
C
ENDIF
C
END
C
C=======================================================================
C
SUBROUTINE COEF_LOEWDIN(NDIF)
C
C This subroutine computes the coefficients for the Loewdin expansion
C of the multiple scattering series. These coefficients are
C expressed as C_LOW(K) where K is the multiple scattering order.
C REN is the value of the mixing (or renormalization) parameter.
C NDIF is the scattering order at which the series is truncated,
C so that K varies from 0 to NDIF.
C
C Corresponds to parameter I_REN = 5
C
C Notation: X(K,N) = X_n(omega,k)
C
C
C Last modified : 11 Apr 2019
C
USE DIM_MOD
USE C_RENORM_MOD, C_LOW => C_REN
USE RENORM_MOD
C
COMPLEX X(0:NDIF_M,0:NDIF_M),SUM_L,POWER
COMPLEX REN,ZEROC,ONEC,IC
C
C
ZEROC=(0.,0.)
ONEC=(1.,0.)
IC=(0.,1.)
C
REN=REN_R+IC*REN_I ! omega
C
C Initialisation of renormalization coefficients
C
DO J=0,NDIF
C_LOW(J)=ZEROC
ENDDO
C
C Computing the X(N,K) coefficients, with K <= N
C
POWER=ONEC/REN
DO N=0,NDIF
POWER=POWER*REN ! omega^n
IF(N.EQ.0) THEN
X(N,0)=ONEC
ELSE
X(N,0)=ZEROC
ENDIF
DO K=1,NDIF
IF(K.GT.N) THEN
X(N,K)=ZEROC
ELSEIF(K.EQ.N) THEN
X(N,K)=POWER*X(N-1,K-1)
ELSE
X(N,K)=X(N-1,K)*(REN-POWER) + POWER*X(N-1,K-1)
ENDIF
ENDDO
ENDDO
C
C Calculation of L_n(omega,NDIF)
C
DO N=0,NDIF
SUM_L=ZEROC
DO K=N,NDIF
SUM_L=SUM_L+X(K,N)
ENDDO
C_LOW(N)=SUM_L
ENDDO
END

2
src/msspec/utils.py
View File

@ -19,7 +19,7 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/utils.py
# Last modified: Thu, 27 Feb 2025 16:33:09 +0100
# Last modified: Wed, 26 Feb 2025 11:15:03 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>

5
src/options.mk
View File

@ -3,7 +3,7 @@ PYMAJ = 3
PYMIN = 5
FC = gfortran
F2PY = python -m numpy.f2py --f77exec=$(FC) --f90exec=$(FC)
F2PY = f2py --f77exec=$(FC) --f90exec=$(FC)
NO_VENV = 0
DEBUG = 0
@ -18,8 +18,7 @@ INSTALL_PREFIX = $(HOME)/.local
GFORTRAN_FFLAGS = -O2 -ffast-math
GFORTRAN_FFLAGS_DBG = -g -Wall -Wextra -Warray-temporaries -Wconversion
GFORTRAN_FFLAGS_DBG += -fbacktrace -ffree-line-length-0 -fcheck=all
GFORTRAN_FFLAGS_DBG += -ffpe-trap=zero,overflow,underflow,invalid,denormal
GFORTRAN_FFLAGS_DBG += -finit-real=nan
GFORTRAN_FFLAGS_DBG += -ffpe-trap=zero,overflow,underflow -finit-real=nan
################################################################################
################################################################################

2
src/pip.freeze
View File

@ -1,7 +1,6 @@
ase
h5py
ipython
looseversion
lxml
matplotlib
numpy
@ -9,6 +8,5 @@ Pint
pandas
pycairo
scipy
setuptools==73.0.1
setuptools-scm
terminaltables