mirror of
https://github.com/LCPQ/quantum_package
synced 2024-06-22 05:02:15 +02:00
Merge branch 'Pseudo'
This commit is contained in:
Thomas Applencourt
commit
1806b97cb8
|
|
@ -8,4 +8,4 @@ include Input_determinants_by_hand;;
|
|||
include Input_electrons;;
|
||||
include Input_mo_basis;;
|
||||
include Input_nuclei;;
|
||||
include Input_auto_generated;;
|
||||
include Input_auto_generated;;
|
||||
|
|
|
|||
|
|
@ -105,7 +105,7 @@ def write_ezfio(res, filename):
|
|||
# Transformt H1 into H
|
||||
import re
|
||||
p = re.compile(ur'(\d*)$')
|
||||
label = [p.sub("", x.name) for x in res.geometry]
|
||||
label = [p.sub("", x.name).capitalize() for x in res.geometry]
|
||||
ezfio.set_nuclei_nucl_label(label)
|
||||
|
||||
ezfio.set_nuclei_nucl_coord(coord_x + coord_y + coord_z)
|
||||
|
|
|
|||
27
scripts/ezfio_interface/upgrade_1.0_2.0.sh
Executable file
27
scripts/ezfio_interface/upgrade_1.0_2.0.sh
Executable file
|
|
@ -0,0 +1,27 @@
|
|||
#!/bin/bash
|
||||
# Convert a old ezfio file (with option.irp.f ezfio_default)
|
||||
# into a new EZFIO.cfg type
|
||||
|
||||
# Hartree Fock
|
||||
# Changin the case, don't know if is needed or not
|
||||
mv $1/Hartree_Fock $1/hartree_fock 2> /dev/null
|
||||
|
||||
mv $1/hartree_Fock/thresh_SCF $1/hartree_fock/thresh_scf 2> /dev/null
|
||||
|
||||
# BiInts
|
||||
mv $1/bi_integrals $1/bielect_integrals 2> /dev/null
|
||||
|
||||
if [ -f $1/bielect_integrals/read_ao_integrals ]; then
|
||||
if [ `cat $1/bielect_integrals/read_ao_integrals` -eq "True" ]
|
||||
then
|
||||
echo "Read" > $1/bielect_integrals/disk_access_ao_integrals
|
||||
|
||||
elif [ `cat bielect_integrals/write_ao_integrals` -eq "True" ]
|
||||
then
|
||||
echo "Write" > $1/bielect_integrals/disk_access_ao_integrals
|
||||
|
||||
else
|
||||
echo "None" > $1/bielect_integrals/disk_access_ao_integrals
|
||||
|
||||
fi
|
||||
fi
|
||||
|
|
@ -42,9 +42,9 @@ then
|
|||
echo "ERROR"
|
||||
exit 1
|
||||
fi
|
||||
${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" $atoms
|
||||
|
||||
|
||||
|
||||
|
||||
#${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}"
|
||||
|
||||
${EMSL_API_ROOT}/EMSL_api.py get_basis_data --save --path="${tmpfile}" --basis="${basis}" --db_path="${EMSL_API_ROOT}/db/Pseudo.db"
|
||||
12
scripts/pseudo/create_ez.sh
Executable file
12
scripts/pseudo/create_ez.sh
Executable file
|
|
@ -0,0 +1,12 @@
|
|||
#!/bin/bash
|
||||
# $1 name
|
||||
# $2 mult
|
||||
|
||||
echo "name" $1
|
||||
echo "basis" $2
|
||||
echo "mul" $3
|
||||
echo "\`get_basis.sh\` need to be changed"
|
||||
|
||||
rm -R $1.ezfio
|
||||
qp_create_ezfio_from_xyz $1.xyz -b $2 -m $3
|
||||
~/quantum_package/scripts/pseudo/put_pseudo_in_ezfio.py $1.ezfio
|
||||
118
scripts/pseudo/elts_num_ele.py
Normal file
118
scripts/pseudo/elts_num_ele.py
Normal file
|
|
@ -0,0 +1,118 @@
|
|||
name_to_elec = {"H": 1,
|
||||
"He": 2,
|
||||
"Li": 3,
|
||||
"Be": 4,
|
||||
"B": 5,
|
||||
"C": 6,
|
||||
"N": 7,
|
||||
"O": 8,
|
||||
"F": 9,
|
||||
"Ne": 10,
|
||||
"Na": 11,
|
||||
"Mg": 12,
|
||||
"Al": 13,
|
||||
"Si": 14,
|
||||
"P": 15,
|
||||
"S": 16,
|
||||
"Cl": 17,
|
||||
"Ar": 18,
|
||||
"K": 19,
|
||||
"Ca": 20,
|
||||
"Sc": 21,
|
||||
"Ti": 22,
|
||||
"V": 23,
|
||||
"Cr": 24,
|
||||
"Mn": 25,
|
||||
"Fe": 26,
|
||||
"Co": 27,
|
||||
"Ni": 28,
|
||||
"Cu": 29,
|
||||
"Zn": 30,
|
||||
"Ga": 31,
|
||||
"Ge": 32,
|
||||
"As": 33,
|
||||
"Se": 34,
|
||||
"Br": 35,
|
||||
"Kr": 36,
|
||||
"Rb": 37,
|
||||
"Sr": 38,
|
||||
"Y": 39,
|
||||
"Zr": 40,
|
||||
"Nb": 41,
|
||||
"Mo": 42,
|
||||
"Tc": 43,
|
||||
"Ru": 44,
|
||||
"Rh": 45,
|
||||
"Pd": 46,
|
||||
"Ag": 47,
|
||||
"Cd": 48,
|
||||
"In": 49,
|
||||
"Sn": 50,
|
||||
"Sb": 51,
|
||||
"Te": 52,
|
||||
"I": 53,
|
||||
"Xe": 54,
|
||||
"Cs": 55,
|
||||
"Ba": 56,
|
||||
"La": 57,
|
||||
"Ce": 58,
|
||||
"Pr": 59,
|
||||
"Nd": 60,
|
||||
"Pm": 61,
|
||||
"Sm": 62,
|
||||
"Eu": 63,
|
||||
"Gd": 64,
|
||||
"Tb": 65,
|
||||
"Dy": 66,
|
||||
"Ho": 67,
|
||||
"Er": 68,
|
||||
"Tm": 69,
|
||||
"Yb": 70,
|
||||
"Lu": 71,
|
||||
"Hf": 72,
|
||||
"Ta": 73,
|
||||
"W": 74,
|
||||
"Re": 75,
|
||||
"Os": 76,
|
||||
"Ir": 77,
|
||||
"Pt": 78,
|
||||
"Au": 79,
|
||||
"Hg": 80,
|
||||
"Tl": 81,
|
||||
"Pb": 82,
|
||||
"Bi": 83,
|
||||
"Po": 84,
|
||||
"At": 85,
|
||||
"Rn": 86,
|
||||
"Fr": 87,
|
||||
"Ra": 88,
|
||||
"Ac": 89,
|
||||
"Th": 90,
|
||||
"Pa": 91,
|
||||
"U": 92,
|
||||
"Np": 93,
|
||||
"Pu": 94,
|
||||
"Am": 95,
|
||||
"Cm": 96,
|
||||
"Bk": 97,
|
||||
"Cf": 98,
|
||||
"Es": 99,
|
||||
"Fm": 100,
|
||||
"Md": 101,
|
||||
"No": 102,
|
||||
"Lr": 103,
|
||||
"Rf": 104,
|
||||
"Db": 105,
|
||||
"Sg": 106,
|
||||
"Bh": 107,
|
||||
"Hs": 108,
|
||||
"Mt": 109,
|
||||
"Ds": 110,
|
||||
"Rg": 111,
|
||||
"Cn": 112,
|
||||
"Uut": 113,
|
||||
"Fl": 114,
|
||||
"Uup": 115,
|
||||
"Lv": 116,
|
||||
"Uus": 117,
|
||||
"Uuo": 118}
|
||||
336
scripts/pseudo/put_pseudo_in_ezfio.py
Executable file
336
scripts/pseudo/put_pseudo_in_ezfio.py
Executable file
|
|
@ -0,0 +1,336 @@
|
|||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Create the pseudo potential for a given atom
|
||||
|
||||
Usage:
|
||||
put_pseudo_in_ezfio.py <ezfio_path>
|
||||
|
||||
Help:
|
||||
atom is the Abreviation of the atom
|
||||
"""
|
||||
|
||||
|
||||
import os
|
||||
import sys
|
||||
from docopt import docopt
|
||||
|
||||
from subprocess import Popen, PIPE
|
||||
|
||||
qpackage_root = os.environ['QPACKAGE_ROOT']
|
||||
|
||||
EZFIO = "{0}/EZFIO".format(qpackage_root)
|
||||
sys.path = [EZFIO + "/Python"] + sys.path
|
||||
|
||||
from ezfio import ezfio
|
||||
|
||||
import re
|
||||
p = re.compile(ur'\|(\d+)><\d+\|')
|
||||
|
||||
|
||||
def get_pseudo_str(l_atom):
|
||||
"""
|
||||
Run EMSL_local for geting the str of the speudo potential
|
||||
|
||||
str_ele :
|
||||
Element Symbol: Na
|
||||
Number of replaced protons: 10
|
||||
Number of projectors: 2
|
||||
|
||||
Pseudopotential data:
|
||||
|
||||
Local component:
|
||||
Coeff. r^n Exp.
|
||||
1.00000000 -1 5.35838717
|
||||
5.35838717 1 3.67918975
|
||||
-2.07764789 0 1.60507673
|
||||
|
||||
Non-local component:
|
||||
Coeff. r^n Exp. Proj.
|
||||
10.69640234 0 1.32389367 |0><0|
|
||||
10.11238853 0 1.14052020 |1><1|
|
||||
"""
|
||||
|
||||
EMSL_root = "{0}/EMSL_Basis/".format(qpackage_root)
|
||||
EMSL_path = "{0}/EMSL_api.py".format(EMSL_root)
|
||||
db_path = "{0}/db/Pseudo.db".format(EMSL_root)
|
||||
|
||||
str_ = ""
|
||||
|
||||
for a in l_atom:
|
||||
l_cmd_atom = ["--atom", a]
|
||||
|
||||
l_cmd_head = [EMSL_path, "get_basis_data",
|
||||
"--db_path", db_path,
|
||||
"--basis", "BFD-Pseudo"]
|
||||
|
||||
process = Popen(l_cmd_head + l_cmd_atom, stdout=PIPE, stderr=PIPE)
|
||||
|
||||
stdout, _ = process.communicate()
|
||||
str_ += stdout.strip() + "\n"
|
||||
|
||||
return str_
|
||||
|
||||
|
||||
def get_v_n_dz_local(str_ele):
|
||||
"""
|
||||
From a str_ele of the pseudo (aka only one ele in the str)
|
||||
get the list ussefull for the Local potential : v_k n_k and dz_k
|
||||
"""
|
||||
l_v_k = []
|
||||
l_n_k = []
|
||||
l_dz_k = []
|
||||
|
||||
for l in str_ele.splitlines():
|
||||
try:
|
||||
v, n, dz = l.split()
|
||||
v = float(v)
|
||||
n = int(n)
|
||||
dz = float(dz)
|
||||
except ValueError:
|
||||
pass
|
||||
else:
|
||||
l_v_k.append(v)
|
||||
l_n_k.append(n)
|
||||
l_dz_k.append(dz)
|
||||
|
||||
return l_v_k, l_n_k, l_dz_k
|
||||
|
||||
|
||||
def get_v_n_dz_l_nonlocal(str_ele):
|
||||
"""
|
||||
From a str_ele of the pseudo (aka only one ele in the str)
|
||||
get the list ussefull for the non Local potential
|
||||
v_kl (v, l)
|
||||
n_k (v, l)
|
||||
dz_k (dz ,l)
|
||||
"""
|
||||
l_v_kl = []
|
||||
l_n_kl = []
|
||||
l_dz_kl = []
|
||||
|
||||
for l in str_ele.splitlines():
|
||||
try:
|
||||
v, n, dz, proj = l.split()
|
||||
v = float(v)
|
||||
n = int(n)
|
||||
dz = float(dz)
|
||||
l = int(p.match(proj).group(1))
|
||||
|
||||
except ValueError:
|
||||
pass
|
||||
else:
|
||||
l_v_kl.append([v])
|
||||
l_n_kl.append([n])
|
||||
l_dz_kl.append([dz])
|
||||
|
||||
if not l_v_kl:
|
||||
l_v_kl.append([0.])
|
||||
l_n_kl.append([0])
|
||||
l_dz_kl.append([0.])
|
||||
|
||||
return l_v_kl, l_n_kl, l_dz_kl
|
||||
|
||||
|
||||
def get_zeff_alpha_beta(str_ele):
|
||||
"""
|
||||
Return the the zeff, alpha num elec and beta num elec
|
||||
Assert ezfio_set_file alredy defined
|
||||
"""
|
||||
|
||||
import re
|
||||
|
||||
# ___
|
||||
# | ._ o _|_
|
||||
# _|_ | | | |_
|
||||
#
|
||||
|
||||
# ~#~#~#~#~#~#~ #
|
||||
# s t r _ e l e #
|
||||
# ~#~#~#~#~#~#~ #
|
||||
|
||||
# m = re.search('Element Symbol: ([a-zA-Z]+)', str_ele)
|
||||
# name = m.group(1).capitalize()
|
||||
name = str_ele.split("\n")[0].strip().capitalize()
|
||||
|
||||
m = re.search('Number of replaced protons: (\d+)', str_ele)
|
||||
z_remove = int(m.group(1))
|
||||
|
||||
# _
|
||||
# |_) _. ._ _ _
|
||||
# | (_| | _> (/_
|
||||
#
|
||||
|
||||
from elts_num_ele import name_to_elec
|
||||
z = name_to_elec[name]
|
||||
|
||||
z_eff = z - z_remove
|
||||
|
||||
alpha = (z_remove / 2)
|
||||
beta = (z_remove / 2)
|
||||
|
||||
# _
|
||||
# |_) _ _|_ ._ ._
|
||||
# | \ (/_ |_ |_| | | |
|
||||
#
|
||||
|
||||
return [z_eff, alpha, beta]
|
||||
|
||||
|
||||
def add_zero(array, size, type):
|
||||
for add in xrange(len(array), size):
|
||||
array.append([type(0)])
|
||||
|
||||
return array
|
||||
|
||||
|
||||
def make_it_square(matrix, dim, type=float):
|
||||
"""
|
||||
matix the matrix to squate
|
||||
dim array [lmax, kmax]
|
||||
type the null value you want
|
||||
[[[28.59107316], [19.37583724]], [[50.25646328]]]
|
||||
=>
|
||||
[[[28.59107316], [19.37583724]], [[50.25646328], [0.0]]]
|
||||
"""
|
||||
|
||||
lmax = dim[0]
|
||||
kmax = dim[1]
|
||||
|
||||
for l_list in matrix:
|
||||
|
||||
l_list = add_zero(l_list, lmax, type)
|
||||
|
||||
for k_list in list_:
|
||||
k_list = add_zero(k_list, kmax, type)
|
||||
|
||||
return matrix
|
||||
|
||||
if __name__ == "__main__":
|
||||
arguments = docopt(__doc__)
|
||||
# ___
|
||||
# | ._ o _|_
|
||||
# _|_ | | | |_
|
||||
#
|
||||
|
||||
# ~#~#~#~#~ #
|
||||
# E Z F I O #
|
||||
# ~#~#~#~#~ #
|
||||
|
||||
ezfio_path = arguments["<ezfio_path>"]
|
||||
ezfio_path = os.path.expanduser(ezfio_path)
|
||||
ezfio_path = os.path.expandvars(ezfio_path)
|
||||
ezfio_path = os.path.abspath(ezfio_path)
|
||||
|
||||
ezfio.set_file("{0}".format(ezfio_path))
|
||||
|
||||
# ~#~#~#~#~#~#~#~#~#~#~ #
|
||||
# P s e u d o _ d a t a #
|
||||
# ~#~#~#~#~#~#~#~#~#~#~ #
|
||||
|
||||
l_ele = ezfio.get_nuclei_nucl_label()
|
||||
str_ = get_pseudo_str(l_ele)
|
||||
|
||||
# _
|
||||
# |_) _. ._ _ _
|
||||
# | (_| | _> (/_
|
||||
#
|
||||
|
||||
l_str_ele = [str_ele for str_ele in str_.split("Element Symbol: ")
|
||||
if str_ele]
|
||||
|
||||
for i in "l_zeff v_k n_k dz_k v_kl n_kl dz_kl".split():
|
||||
exec("{0} = []".format(i))
|
||||
|
||||
alpha_tot = 0
|
||||
beta_tot = 0
|
||||
|
||||
for str_ele in l_str_ele:
|
||||
|
||||
# ~#~#~#~#~ #
|
||||
# S p l i t #
|
||||
# ~#~#~#~#~ #
|
||||
|
||||
l = str_ele.find("Local component:")
|
||||
nl = str_ele.find("Non-local component")
|
||||
|
||||
# ~#~#~#~#~ #
|
||||
# L o c a l #
|
||||
# ~#~#~#~#~ #
|
||||
|
||||
l_v, l_n, l_dz = get_v_n_dz_local(str_ele[l:nl])
|
||||
|
||||
v_k.append(l_v)
|
||||
n_k.append(l_n)
|
||||
dz_k.append(l_dz)
|
||||
|
||||
# ~#~#~#~#~#~#~#~#~ #
|
||||
# N o n _ L o c a l #
|
||||
# ~#~#~#~#~#~#~#~#~ #
|
||||
|
||||
l_v_kl, l_n_kl, l_dz_kl = get_v_n_dz_l_nonlocal(str_ele[nl:])
|
||||
|
||||
v_kl.append(l_v_kl)
|
||||
n_kl.append(l_n_kl)
|
||||
dz_kl.append(l_dz_kl)
|
||||
|
||||
# ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
|
||||
# Z _ e f f , a l p h a / b e t a _ e l e c #
|
||||
# ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
|
||||
|
||||
zeff, alpha, beta = get_zeff_alpha_beta(str_ele)
|
||||
|
||||
alpha_tot += alpha
|
||||
beta_tot += beta
|
||||
l_zeff.append(zeff)
|
||||
# _
|
||||
# /\ _| _| _|_ _ _ _ _|_ o _
|
||||
# /--\ (_| (_| |_ (_) (/_ /_ | | (_)
|
||||
#
|
||||
|
||||
# ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
|
||||
# Z _ e f f , a l p h a / b e t a _ e l e c #
|
||||
# ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
|
||||
|
||||
ezfio.nuclei_nucl_charge = l_zeff
|
||||
|
||||
alpha_tot = ezfio.get_electrons_elec_alpha_num() - alpha_tot
|
||||
beta_tot = ezfio.get_electrons_elec_beta_num() - beta_tot
|
||||
|
||||
ezfio.electrons_elec_alpha_num = alpha_tot
|
||||
ezfio.electrons_elec_beta_num = beta_tot
|
||||
|
||||
# Change all the array 'cause EZFIO
|
||||
# v_kl (v, l) => v_kl(l,v)
|
||||
# v_kl => zip(*_v_kl)
|
||||
# [[7.0, 79.74474797, -49.45159098], [1.0, 5.41040609, -4.60151975]]
|
||||
# [(7.0, 1.0), (79.74474797, 5.41040609), (-49.45159098, -4.60151975)]
|
||||
|
||||
# ~#~#~#~#~ #
|
||||
# L o c a l #
|
||||
# ~#~#~#~#~ #
|
||||
|
||||
klocmax = max([len(i) for i in v_k])
|
||||
ezfio.pseudo_klocmax = klocmax
|
||||
|
||||
ezfio.pseudo_v_k = zip(*v_k)
|
||||
ezfio.pseudo_n_k = zip(*n_k)
|
||||
ezfio.pseudo_dz_k = zip(*dz_k)
|
||||
|
||||
# ~#~#~#~#~#~#~#~#~ #
|
||||
# N o n _ L o c a l #
|
||||
# ~#~#~#~#~#~#~#~#~ #
|
||||
|
||||
lmax = max([len(i) for i in v_kl])
|
||||
kmax = max([len(sublist) for list_ in v_kl for sublist in list_])
|
||||
|
||||
ezfio.pseudo_lmaxpo = lmax
|
||||
ezfio.pseudo_kmax = kmax
|
||||
|
||||
v_kl = make_it_square(v_kl, [lmax, kmax])
|
||||
n_kl = make_it_square(n_kl, [lmax, kmax], int)
|
||||
dz_kl = make_it_square(dz_kl, [lmax, kmax])
|
||||
|
||||
ezfio.pseudo_v_kl = zip(*v_kl)
|
||||
ezfio.pseudo_n_kl = zip(*n_kl)
|
||||
ezfio.pseudo_dz_kl = zip(*dz_kl)
|
||||
|
|
@ -1,6 +1,6 @@
|
|||
# Define here all new external source files and objects.Don't forget to prefix the
|
||||
# object files with IRPF90_temp/
|
||||
SRC=
|
||||
OBJ=
|
||||
SRC=int.f90
|
||||
OBJ=IRPF90_temp/int.o
|
||||
|
||||
include $(QPACKAGE_ROOT)/src/Makefile.common
|
||||
include $(QPACKAGE_ROOT)/src/Makefile.common
|
||||
|
|
@ -57,44 +57,80 @@ Documentation
|
|||
array of the mono electronic hamiltonian on the MOs basis
|
||||
: sum of the kinetic and nuclear electronic potential
|
||||
|
||||
`a_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L252>`_
|
||||
Undocumented
|
||||
|
||||
`b_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L257>`_
|
||||
Undocumented
|
||||
|
||||
`ddfact2 <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L243>`_
|
||||
Undocumented
|
||||
|
||||
`erf0 <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L105>`_
|
||||
Undocumented
|
||||
|
||||
`gammln <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L271>`_
|
||||
Undocumented
|
||||
|
||||
`gammp <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L133>`_
|
||||
Undocumented
|
||||
|
||||
`gcf <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L211>`_
|
||||
Undocumented
|
||||
|
||||
`gser <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L167>`_
|
||||
Undocumented
|
||||
|
||||
`rinteg <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L47>`_
|
||||
Undocumented
|
||||
|
||||
`rintgauss <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L31>`_
|
||||
Undocumented
|
||||
|
||||
`sabpartial <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/need.irp.f#L2>`_
|
||||
Undocumented
|
||||
|
||||
`orthonormalize_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/orthonormalize.irp.f#L1>`_
|
||||
Undocumented
|
||||
|
||||
`ao_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L1>`_
|
||||
interaction nuclear electron
|
||||
|
||||
`ao_nucl_elec_integral_per_atom <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L67>`_
|
||||
`ao_nucl_elec_integral_per_atom <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L243>`_
|
||||
ao_nucl_elec_integral_per_atom(i,j,k) = -<AO(i)|1/|r-Rk|AO(j)>
|
||||
where Rk is the geometry of the kth atom
|
||||
|
||||
`give_polynom_mult_center_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L214>`_
|
||||
`ao_nucl_elec_integral_pseudo <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L72>`_
|
||||
interaction nuclear electron
|
||||
|
||||
`give_polynom_mult_center_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L389>`_
|
||||
Undocumented
|
||||
|
||||
`i_x1_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L342>`_
|
||||
`i_x1_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L517>`_
|
||||
Undocumented
|
||||
|
||||
`i_x2_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L413>`_
|
||||
`i_x2_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L588>`_
|
||||
Undocumented
|
||||
|
||||
`int_gaus_pol <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L484>`_
|
||||
`int_gaus_pol <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L659>`_
|
||||
Undocumented
|
||||
|
||||
`nai_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L135>`_
|
||||
`nai_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L310>`_
|
||||
Undocumented
|
||||
|
||||
`v_e_n <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L465>`_
|
||||
`v_e_n <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L640>`_
|
||||
Undocumented
|
||||
|
||||
`v_phi <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L529>`_
|
||||
`v_phi <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L704>`_
|
||||
Undocumented
|
||||
|
||||
`v_r <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L513>`_
|
||||
`v_r <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L688>`_
|
||||
Undocumented
|
||||
|
||||
`v_theta <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L542>`_
|
||||
`v_theta <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L717>`_
|
||||
Undocumented
|
||||
|
||||
`wallis <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L558>`_
|
||||
`wallis <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L733>`_
|
||||
Undocumented
|
||||
|
||||
`mo_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_mo_ints.irp.f#L1>`_
|
||||
|
|
@ -215,5 +251,8 @@ Documentation
|
|||
array of the integrals of MO_i * y^2 MO_j
|
||||
array of the integrals of MO_i * z^2 MO_j
|
||||
|
||||
`compute_integrals_pseudo <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/test_michel.irp.f#L58>`_
|
||||
Undocumented
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
2101
src/MonoInts/int.f90
Normal file
2101
src/MonoInts/int.f90
Normal file
File diff suppressed because it is too large
Load Diff
|
|
@ -4,63 +4,239 @@
|
|||
END_DOC
|
||||
implicit none
|
||||
double precision :: alpha, beta, gama, delta
|
||||
integer :: i_c,num_A,num_B
|
||||
double precision :: A_center(3),B_center(3),C_center(3)
|
||||
integer :: power_A(3),power_B(3)
|
||||
integer :: i,j,k,l,n_pt_in,m
|
||||
double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
|
||||
integer :: nucl_numC
|
||||
! Important for OpenMP
|
||||
integer :: num_A,num_B
|
||||
double precision :: A_center(3),B_center(3),C_center(3)
|
||||
integer :: power_A(3),power_B(3)
|
||||
integer :: i,j,k,l,n_pt_in,m
|
||||
double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
|
||||
|
||||
ao_nucl_elec_integral = 0.d0
|
||||
ao_nucl_elec_integral = ao_nucl_elec_integral_pseudo ! 0.d0
|
||||
|
||||
! _
|
||||
! /| / |_)
|
||||
! | / | \
|
||||
!
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
|
||||
!$OMP num_A,num_B,Z,c,n_pt_in) &
|
||||
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, &
|
||||
!$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge)
|
||||
n_pt_in = n_pt_max_integrals
|
||||
!$OMP DO SCHEDULE (guided)
|
||||
do j = 1, ao_num
|
||||
power_A(1)= ao_power(j,1)
|
||||
power_A(2)= ao_power(j,2)
|
||||
power_A(3)= ao_power(j,3)
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
|
||||
!$OMP num_A,num_B,Z,c,n_pt_in) &
|
||||
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, &
|
||||
!$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge)
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
!$OMP DO SCHEDULE (guided)
|
||||
|
||||
do j = 1, ao_num
|
||||
num_A = ao_nucl(j)
|
||||
A_center(1) = nucl_coord(num_A,1)
|
||||
A_center(2) = nucl_coord(num_A,2)
|
||||
A_center(3) = nucl_coord(num_A,3)
|
||||
power_A(1:3)= ao_power(j,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
|
||||
do i = 1, ao_num
|
||||
power_B(1)= ao_power(i,1)
|
||||
power_B(2)= ao_power(i,2)
|
||||
power_B(3)= ao_power(i,3)
|
||||
|
||||
num_B = ao_nucl(i)
|
||||
B_center(1) = nucl_coord(num_B,1)
|
||||
B_center(2) = nucl_coord(num_B,2)
|
||||
B_center(3) = nucl_coord(num_B,3)
|
||||
do l=1,ao_prim_num(j)
|
||||
alpha = ao_expo_transp(l,j)
|
||||
power_B(1:3)= ao_power(i,1:3)
|
||||
B_center(1:3) = nucl_coord(num_B,1:3)
|
||||
|
||||
do l=1,ao_prim_num(j)
|
||||
alpha = ao_expo_transp(l,j)
|
||||
|
||||
do m=1,ao_prim_num(i)
|
||||
beta = ao_expo_transp(m,i)
|
||||
c = 0.d0
|
||||
do k = 1, nucl_num
|
||||
double precision :: Z,c
|
||||
Z = nucl_charge(k)
|
||||
C_center(1) = nucl_coord(k,1)
|
||||
C_center(2) = nucl_coord(k,2)
|
||||
C_center(3) = nucl_coord(k,3)
|
||||
c = c+Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in)
|
||||
beta = ao_expo_transp(m,i)
|
||||
|
||||
double precision :: c
|
||||
c = 0.d0
|
||||
|
||||
do k = 1, nucl_num
|
||||
double precision :: Z
|
||||
Z = nucl_charge(k)
|
||||
|
||||
C_center(1:3) = nucl_coord(k,1:3)
|
||||
|
||||
c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in)
|
||||
|
||||
enddo
|
||||
ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + &
|
||||
ao_coef_transp(l,j)*ao_coef_transp(m,i)*c
|
||||
enddo
|
||||
enddo
|
||||
ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) - &
|
||||
ao_coef_transp(l,j)*ao_coef_transp(m,i)*c
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo, (ao_num_align,ao_num)]
|
||||
BEGIN_DOC
|
||||
! interaction nuclear electron
|
||||
END_DOC
|
||||
implicit none
|
||||
double precision :: alpha, beta, gama, delta
|
||||
integer :: num_A,num_B
|
||||
double precision :: A_center(3),B_center(3),C_center(3)
|
||||
integer :: power_A(3),power_B(3)
|
||||
integer :: i,j,k,l,n_pt_in,m
|
||||
double precision :: Vloc, Vpseudo
|
||||
|
||||
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0
|
||||
integer :: thread_num
|
||||
|
||||
ao_nucl_elec_integral_pseudo = 0.d0
|
||||
|
||||
!
|
||||
! | _ _ _. |
|
||||
! |_ (_) (_ (_| |
|
||||
!
|
||||
!! Parameters of the local part of pseudo:
|
||||
|
||||
integer klocmax
|
||||
integer, allocatable :: n_k(:,:)
|
||||
double precision, allocatable :: v_k(:,:), dz_k(:,:)
|
||||
|
||||
call ezfio_get_pseudo_klocmax(klocmax)
|
||||
|
||||
allocate(n_k(nucl_num,klocmax),v_k(nucl_num,klocmax), dz_k(nucl_num,klocmax))
|
||||
|
||||
call ezfio_get_pseudo_v_k(v_k)
|
||||
call ezfio_get_pseudo_n_k(n_k)
|
||||
call ezfio_get_pseudo_dz_k(dz_k)
|
||||
|
||||
!! Dump array
|
||||
integer, allocatable :: n_k_dump(:)
|
||||
double precision, allocatable :: v_k_dump(:), dz_k_dump(:)
|
||||
|
||||
allocate(n_k_dump(1:klocmax), v_k_dump(1:klocmax), dz_k_dump(1:klocmax))
|
||||
|
||||
|
||||
!
|
||||
! |\ | _ ._ | _ _ _. |
|
||||
! | \| (_) | | | (_) (_ (_| |
|
||||
!
|
||||
!! Parameters of non local part of pseudo:
|
||||
|
||||
integer :: kmax,lmax
|
||||
integer, allocatable :: n_kl(:,:,:)
|
||||
double precision, allocatable :: v_kl(:,:,:), dz_kl(:,:,:)
|
||||
|
||||
call ezfio_get_pseudo_lmaxpo(lmax)
|
||||
call ezfio_get_pseudo_kmax(kmax)
|
||||
!lmax plus one -> lmax
|
||||
lmax = lmax - 1
|
||||
|
||||
allocate(n_kl(nucl_num,kmax,0:lmax), v_kl(nucl_num,kmax,0:lmax), dz_kl(nucl_num,kmax,0:lmax))
|
||||
|
||||
call ezfio_get_pseudo_n_kl(n_kl)
|
||||
call ezfio_get_pseudo_v_kl(v_kl)
|
||||
call ezfio_get_pseudo_dz_kl(dz_kl)
|
||||
|
||||
|
||||
!! Dump array
|
||||
integer, allocatable :: n_kl_dump(:,:)
|
||||
double precision, allocatable :: v_kl_dump(:,:), dz_kl_dump(:,:)
|
||||
|
||||
allocate(n_kl_dump(kmax,0:lmax), v_kl_dump(kmax,0:lmax), dz_kl_dump(kmax,0:lmax))
|
||||
|
||||
! _
|
||||
! / _. | _ |
|
||||
! \_ (_| | (_ |_| |
|
||||
!
|
||||
|
||||
write(output_monoints,*) 'Providing the nuclear electron pseudo integrals '
|
||||
|
||||
call wall_time(wall_1)
|
||||
call cpu_time(cpu_1)
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
|
||||
!$OMP num_A,num_B,Z,c,n_pt_in, &
|
||||
!$OMP v_k_dump,n_k_dump, dz_k_dump, n_kl_dump, v_kl_dump, dz_kl_dump, &
|
||||
!$OMP wall_0,wall_2,thread_num, output_monoints) &
|
||||
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, &
|
||||
!$OMP n_pt_max_integrals,ao_nucl_elec_integral_pseudo,nucl_num,nucl_charge, &
|
||||
!$OMP klocmax,lmax,kmax,v_k,n_k, dz_k, n_kl, v_kl, dz_kl, &
|
||||
!$OMP wall_1)
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
!$OMP DO SCHEDULE (guided)
|
||||
|
||||
do j = 1, ao_num
|
||||
|
||||
num_A = ao_nucl(j)
|
||||
power_A(1:3)= ao_power(j,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
|
||||
do i = 1, ao_num
|
||||
|
||||
num_B = ao_nucl(i)
|
||||
power_B(1:3)= ao_power(i,1:3)
|
||||
B_center(1:3) = nucl_coord(num_B,1:3)
|
||||
|
||||
do l=1,ao_prim_num(j)
|
||||
alpha = ao_expo_transp(l,j)
|
||||
|
||||
do m=1,ao_prim_num(i)
|
||||
beta = ao_expo_transp(m,i)
|
||||
double precision :: c
|
||||
c = 0.d0
|
||||
|
||||
do k = 1, nucl_num
|
||||
double precision :: Z
|
||||
Z = nucl_charge(k)
|
||||
|
||||
C_center(1:3) = nucl_coord(k,1:3)
|
||||
|
||||
v_k_dump = v_k(k,1:klocmax)
|
||||
n_k_dump = n_k(k,1:klocmax)
|
||||
dz_k_dump = dz_k(k,1:klocmax)
|
||||
|
||||
c = c + Vloc(klocmax, v_k_dump,n_k_dump, dz_k_dump, &
|
||||
A_center,power_A,alpha,B_center,power_B,beta,C_center)
|
||||
|
||||
|
||||
n_kl_dump = n_kl(k,1:kmax,0:lmax)
|
||||
v_kl_dump = v_kl(k,1:kmax,0:lmax)
|
||||
dz_kl_dump = dz_kl(k,1:kmax,0:lmax)
|
||||
|
||||
c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center)
|
||||
|
||||
enddo
|
||||
ao_nucl_elec_integral_pseudo(i,j) = ao_nucl_elec_integral_pseudo(i,j) + &
|
||||
ao_coef_transp(l,j)*ao_coef_transp(m,i)*c
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall_2)
|
||||
if (thread_num == 0) then
|
||||
if (wall_2 - wall_0 > 1.d0) then
|
||||
wall_0 = wall_2
|
||||
write(output_monoints,*) 100.*float(j)/float(ao_num), '% in ', &
|
||||
wall_2-wall_1, 's'
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
|
||||
! _
|
||||
! | \ _ _. | | _ _ _. _|_ _
|
||||
! |_/ (/_ (_| | | (_) (_ (_| |_ (/_
|
||||
!
|
||||
|
||||
deallocate(n_k,v_k, dz_k)
|
||||
deallocate(n_k_dump,v_k_dump, dz_k_dump)
|
||||
|
||||
deallocate(n_kl,v_kl, dz_kl)
|
||||
deallocate(n_kl_dump,v_kl_dump, dz_kl_dump)
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
|
@ -76,7 +252,6 @@ END_PROVIDER
|
|||
integer :: power_A(3),power_B(3)
|
||||
integer :: i,j,k,l,n_pt_in,m
|
||||
double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
|
||||
integer :: nucl_numC
|
||||
! Important for OpenMP
|
||||
|
||||
ao_nucl_elec_integral_per_atom = 0.d0
|
||||
|
|
|
|||
11
src/MonoInts/pseudo.ezfio_config
Normal file
11
src/MonoInts/pseudo.ezfio_config
Normal file
|
|
@ -0,0 +1,11 @@
|
|||
pseudo
|
||||
do_pseudo logical
|
||||
klocmax integer
|
||||
v_k double precision (nuclei_nucl_num,pseudo_klocmax)
|
||||
n_k integer (nuclei_nucl_num,pseudo_klocmax)
|
||||
dz_k double precision (nuclei_nucl_num,pseudo_klocmax)
|
||||
lmaxpo integer
|
||||
kmax integer
|
||||
v_kl double precision (nuclei_nucl_num,pseudo_kmax,pseudo_lmaxpo)
|
||||
n_kl integer (nuclei_nucl_num,pseudo_kmax,pseudo_lmaxpo)
|
||||
dz_kl double precision (nuclei_nucl_num,pseudo_kmax,pseudo_lmaxpo)
|
||||
200
src/MonoInts/test_michel.irp.f
Normal file
200
src/MonoInts/test_michel.irp.f
Normal file
|
|
@ -0,0 +1,200 @@
|
|||
!!
|
||||
!! Computation of Vps, matrix element of the
|
||||
!! pseudo-potential centered at point C
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
!! Vps= < Phi_A | Vloc(C) + Vpp(C) | Phi_B>
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
!! Phi_M (M=A,B) Cartesian gaussian orbital centered at point M :
|
||||
!! Phi_M = (x-M_x)**n^M_x *(y-M_y)**n^M_y *(z-M_z)**n^M_z exp(-g_M rM**2)
|
||||
!! with rM**2=(x-M_x)**2 + (y-M_y)**2 + (z-M_z)**2
|
||||
!!
|
||||
!!** Vloc(C)= \sum_{k=1}^klocmax v_k rC**n_k exp(-dz_k rC**2)
|
||||
!!
|
||||
!!** Vpp(C)= \sum_{l=0}^lmax v_l(rC) \sum_{m=-l}^{m=l} |Y_lm> <Y_lm|
|
||||
!!
|
||||
!! v_l(rC) = \sum_{k=1}^kmax v_kl rC**n_kl exp(-dz_kl rC**2)
|
||||
!!
|
||||
!! Y_lm : real spherical harmonics:
|
||||
!! Y_00 = sqrt(1/4pi)
|
||||
!! Y_11 = sqrt(3/4pi) xchap
|
||||
!! Y_10 = sqrt(3/4pi) ychap
|
||||
!! Y_1-1 = sqrt(3/4pi) zchap
|
||||
!! Y_22 = sqrt(15/16pi) (xchap**2-ychap**2)
|
||||
!! Y_21 = sqrt(15/4pi) xchap*zchap
|
||||
!! Y_20 = sqrt(15/16pi)(-xchap**2-ychap**2+ 2 zchap**2)
|
||||
!! Y_2-1 = sqrt(15/4pi) ychap*zchap
|
||||
!! Y_2-2 = sqrt(15/4pi) xchap*ychap
|
||||
!! etc.
|
||||
!! xchap=x/r ychap=y/r zchap=z/r
|
||||
!!
|
||||
!! Routine computing <Phi_A|Vpp(C)|Phi_B> :
|
||||
!! function Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c)
|
||||
!! lmax of formula above
|
||||
!! kmax of formula above
|
||||
!! v_kl = array v_kl(kmax_max,0:lmax_max)
|
||||
!! n_kl = array n_kl(kmax_max,0:lmax_max)
|
||||
!! dz_kl = array dz_kl(kmax_max,0:lmax_max)
|
||||
!! n_a(1),n_a(2),n_a(3)
|
||||
!! a(1),a(2),a(3)
|
||||
!! g_a
|
||||
!! n_b(1),n_b(2),n_b(3)
|
||||
!! b(1),b(2),b(3)
|
||||
!! g_b
|
||||
!! c(1),c(2),c(3)
|
||||
!!
|
||||
!! Routine computing <Phi_A|Vloc(C)|Phi_B> :
|
||||
!! function Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c)
|
||||
!! klocmax of formula above
|
||||
!! v_k = array v_k(klocmax_max)
|
||||
!! n_k = array n_k(klocmax_max)
|
||||
!! dz_k= array dz_k(klocmax_max)
|
||||
!! Routine total matrix element <Phi_A|Vloc(C)+Vlpp(C)|Phi_B> :
|
||||
!! function Vps(a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl)
|
||||
!!
|
||||
!! Routines Vps_num, Vpseudo_num, and Vloc_num = brute force numerical
|
||||
!! estimations of the same integrals
|
||||
|
||||
|
||||
program compute_integrals_pseudo
|
||||
implicit none
|
||||
integer n_a(3),n_b(3),npts
|
||||
double precision g_a,g_b,a(3),b(3),c(3)
|
||||
double precision Vpseudo,Vpseudo_num,Vloc,Vloc_num
|
||||
double precision v3,v4
|
||||
|
||||
|
||||
double precision vps,vps_num
|
||||
|
||||
! PSEUDOS
|
||||
integer nptsgridmax,nptsgrid
|
||||
double precision coefs_pseudo,ptsgrid
|
||||
|
||||
double precision rmax
|
||||
double precision time_1,time_2,time_3,time_4,time_5
|
||||
integer kga,kgb,na1,na2,na3,nb1,nb2,nb3
|
||||
|
||||
CALL RANDOM_SEED()
|
||||
|
||||
nptsgrid=50
|
||||
call initpseudos(nptsgrid)
|
||||
|
||||
PROVIDE ezfio_filename
|
||||
|
||||
!
|
||||
! | _ _ _. |
|
||||
! |_ (_) (_ (_| |
|
||||
!
|
||||
|
||||
integer klocmax
|
||||
integer, allocatable :: n_k(:)
|
||||
double precision, allocatable :: v_k(:), dz_k(:)
|
||||
|
||||
call ezfio_get_pseudo_klocmax(klocmax)
|
||||
|
||||
allocate(n_k(klocmax),v_k(klocmax), dz_k(klocmax))
|
||||
|
||||
call ezfio_get_pseudo_v_k(v_k)
|
||||
call ezfio_get_pseudo_n_k(n_k)
|
||||
call ezfio_get_pseudo_dz_k(dz_k)
|
||||
|
||||
print*, "klocmax", klocmax
|
||||
|
||||
print*, "n_k_ezfio", n_k
|
||||
print*, "v_k_ezfio",v_k
|
||||
print*, "dz_k_ezfio", dz_k
|
||||
|
||||
!
|
||||
! |\ | _ ._ | _ _ _. |
|
||||
! | \| (_) | | | (_) (_ (_| |
|
||||
!
|
||||
|
||||
!! Parameters of non local part of pseudo:
|
||||
|
||||
integer :: kmax,lmax
|
||||
integer, allocatable :: n_kl(:,:)
|
||||
double precision, allocatable :: v_kl(:,:), dz_kl(:,:)
|
||||
|
||||
call ezfio_get_pseudo_lmaxpo(lmax)
|
||||
call ezfio_get_pseudo_kmax(kmax)
|
||||
lmax = lmax - 1
|
||||
|
||||
allocate(n_kl(kmax,0:lmax), v_kl(kmax,0:lmax), dz_kl(kmax,0:lmax))
|
||||
|
||||
call ezfio_get_pseudo_n_kl(n_kl)
|
||||
call ezfio_get_pseudo_v_kl(v_kl)
|
||||
call ezfio_get_pseudo_dz_kl(dz_kl)
|
||||
|
||||
|
||||
print*, "lmax",lmax
|
||||
print*, "kmax", kmax
|
||||
|
||||
print*,"n_kl_ezfio", n_kl
|
||||
print*,"v_kl_ezfio", v_kl
|
||||
print*,"dz_kl_ezfio", dz_kl
|
||||
|
||||
! _
|
||||
! / _. | _ |
|
||||
! \_ (_| | (_ |_| |
|
||||
!
|
||||
|
||||
! write(*,*)'a?'
|
||||
! read*,a(1),a(2),a(3)
|
||||
!write(*,*)'b?'
|
||||
!read*,b(1),b(2),b(3)
|
||||
! b(1)=-0.1d0
|
||||
! b(2)=-0.2d0
|
||||
! b(3)=0.3d0
|
||||
! !write(*,*)'a?'
|
||||
! !read*,c(1),c(2),c(3)
|
||||
! c(1)=0.1d0
|
||||
! c(2)=0.2d0
|
||||
! c(3)=0.3d0
|
||||
|
||||
a(1)= 0.d0
|
||||
a(2)= 0.d0
|
||||
a(3)= 0.d0
|
||||
|
||||
b(1)= 0.d0
|
||||
b(2)= 0.d0
|
||||
b(3)= 0.d0
|
||||
|
||||
c(1)= 0.d0
|
||||
c(2)= 0.d0
|
||||
c(3)= 0.d0
|
||||
|
||||
print*,'ntps? rmax for brute force integration'
|
||||
read*,npts,rmax
|
||||
|
||||
do kga=0,5
|
||||
g_a=10.d0**kga
|
||||
do kgb=0,5
|
||||
g_b=10.d0**kgb
|
||||
|
||||
do na1=0,1
|
||||
do na2=0,1
|
||||
do na3=0,1
|
||||
do nb1=0,1
|
||||
do nb2=0,1
|
||||
do nb3=0,1
|
||||
n_a(1)=na1
|
||||
n_a(2)=na2
|
||||
n_a(3)=na3
|
||||
n_b(1)=nb1
|
||||
n_b(2)=nb2
|
||||
n_b(3)=nb3
|
||||
|
||||
v3=Vps(a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl)
|
||||
v4=Vps_num(npts,rmax,a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl)
|
||||
print*,'Vps= ',v3,' Vps_num=',v4,' diff=',v4-v3
|
||||
write(33,'(3f10.6)')v3,v4,v4-v3
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end
|
||||
Loading…
Reference in New Issue
Block a user