diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py new file mode 100755 index 00000000..792df9de --- /dev/null +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -0,0 +1,183 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +Create the pseudo potential for a given atom + +Usage: + put_pseudo_in_ezfio.py --ezfio= --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 + """ + + 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) + + l_cmd_atom = [] + 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() + return stdout.strip() + + +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)) + l_n_kl.append((n, l)) + l_dz_kl.append((dz, l)) + + return l_v_kl, l_n_kl, l_dz_kl + + +if __name__ == "__main__": + arguments = docopt(__doc__) + + # ___ + # | ._ o _|_ + # _|_ | | | |_ + # + + # ~#~#~#~#~ # + # E Z F I O # + # ~#~#~#~#~ # + + ezfio_path = arguments["--ezfio"] + 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 = arguments["--atom"] + str_ = get_pseudo_str(l_ele) + + # _ + # |_) _. ._ _ _ + # | (_| | _> (/_ + # + + l_str_ele = [str_ele for str_ele in str_.split("Element Symbol: ") + if str_ele] + + 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 # + # ~#~#~#~#~ # + + print "local" + + l_v, l_n, l_dz = get_v_n_dz_local(str_ele[l:nl]) + + print l_v + print l_n + print l_dz + + ezfio.pseudo_klocmax = len(l_v) + ezfio.pseudo_v_k = l_v + ezfio.pseudo_n_k = l_n + ezfio.pseudo_dz_k = l_dz + + # ~#~#~#~#~#~#~#~#~ # + # N o n _ L o c a l # + # ~#~#~#~#~#~#~#~#~ # + + print "non local" + + l_v_kl, l_n_kl, l_dz_kl = get_v_n_dz_l_nonlocal(str_ele[nl:]) + + print l_v_kl + print l_n_kl + print l_dz_kl + + if l_v_kl: + ezfio.pseudo_lmax = max([i[1] for i in l_v_kl]) + 1 + ezfio.pseudo_kmax = len(l_v_kl) + ezfio.pseudo_v_kl = l_v_kl + ezfio.pseudo_n_kl = l_n_kl + ezfio.pseudo_dz_kl = l_dz_kl diff --git a/setup_environment.sh b/setup_environment.sh index c3dc4194..59f3af70 100755 --- a/setup_environment.sh +++ b/setup_environment.sh @@ -17,7 +17,9 @@ export QPACKAGE_ROOT=\$( cd \$(dirname "\${BASH_SOURCE}") ; pwd -P ) export LD_LIBRARY_PATH="\${QPACKAGE_ROOT}"/lib:\${LD_LIBRARY_PATH} export LIBRARY_PATH="\${QPACKAGE_ROOT}"/lib:\${LIBRARY_PATH} export C_INCLUDE_PATH="\${QPACKAGE_ROOT}"/include:\${C_INCLUDE_PATH} -export PYTHONPATH=\${PYTHONPATH}:"\${QPACKAGE_ROOT}"/scripts:"\${QPACKAGE_ROOT}"/scripts/ezfio_interface +export PYTHONPATH=\${PYTHONPATH}:"\${QPACKAGE_ROOT}"/scripts +export PYTHONPATH=\${PYTHONPATH}:"\${QPACKAGE_ROOT}"/scripts/ezfio_interface +export PYTHONPATH=\${PYTHONPATH}:"\${QPACKAGE_ROOT}"/scripts/pseudo export PATH=\${PATH}:"\${QPACKAGE_ROOT}"/scripts:"\${QPACKAGE_ROOT}"/scripts/ezfio_interface export PATH=\${PATH}:"\${QPACKAGE_ROOT}"/bin export PATH=\${PATH}:"\${QPACKAGE_ROOT}"/ocaml diff --git a/src/Pseudo/NEEDED_MODULES b/src/Pseudo/NEEDED_MODULES index e69de29b..542760a9 100644 --- a/src/Pseudo/NEEDED_MODULES +++ b/src/Pseudo/NEEDED_MODULES @@ -0,0 +1 @@ +Ezfio_files Nuclei Output Utils \ No newline at end of file diff --git a/src/Pseudo/int.f90 b/src/Pseudo/int.f90 index 16aca8e0..b67e6102 100644 --- a/src/Pseudo/int.f90 +++ b/src/Pseudo/int.f90 @@ -1,60 +1,3 @@ -!! -!! 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> : -!! 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 : -!! 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 : -!! 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 - - !! Vps= !! !! with: Vloc(C)=\sum_{k=1}^klocmax v_k rC**n_k exp(-dz_k rC**2) @@ -1019,15 +962,6 @@ end ! Y_l^m(theta,phi) = i^(m+|m|) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 ! P_l^|m|(cos(theta)) exp(i m phi) - double precision function fact(n) - implicit double precision(a-h,o-z) - fact=1.d0 - if(n.eq.0)return - do i=1,n - fact=fact*dfloat(i) - enddo - end - subroutine erreur(x,n,rmoy,error) implicit double precision(a-h,o-z) dimension x(n) diff --git a/src/Pseudo/pseudo.ezfio_config b/src/Pseudo/pseudo.ezfio_config new file mode 100644 index 00000000..1ed75773 --- /dev/null +++ b/src/Pseudo/pseudo.ezfio_config @@ -0,0 +1,10 @@ +pseudo + klocmax integer + v_k double precision (pseudo_klocmax) + n_k integer (pseudo_klocmax) + dz_k double precision (pseudo_klocmax) + lmax integer + kmax integer + v_kl double precision (pseudo_kmax,pseudo_lmax) + n_kl integer (pseudo_kmax,pseudo_lmax) + dz_kl double precision (pseudo_kmax,pseudo_lmax) \ No newline at end of file diff --git a/src/Pseudo/test_michel.irp.f b/src/Pseudo/test_michel.irp.f new file mode 100644 index 00000000..e290ea4a --- /dev/null +++ b/src/Pseudo/test_michel.irp.f @@ -0,0 +1,188 @@ +!! +!! 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> : +!! 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 : +!! 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 : +!! 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_lmax(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 + + 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