From 3e97da35f5d4948fae7387320f0172004079bfc5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 12 Oct 2009 17:37:07 +0200 Subject: [PATCH] First working version. --- bin/eplf_input.py | 26 - bin/to_ezfio.py | 180 +++ eplf.config | 40 + src/Makefile | 2 +- src/ao.irp.f | 93 +- src/ao_oneD_point.irp.f | 66 +- src/debug_eplf.irp.f | 8 +- src/det.irp.f | 33 + src/electrons.irp.f | 51 +- src/eplf.irp.f | 2 +- src/eplf_function.irp.f | 150 +-- src/eplf_grid.irp.f | 43 +- src/ezfio_interface.irp.f | 122 +++ src/file.irp.f | 35 +- src/mo.irp.f | 126 +-- src/mpi.irp.f | 2 - src/nuclei.irp.f | 42 +- test/QCIO_File/ao/num_orb_sym.gz | Bin 63 -> 0 bytes test/QCIO_File/ao/transformation.gz | Bin 201 -> 0 bytes test/QCIO_File/basis/atom.gz | Bin 68 -> 0 bytes test/QCIO_File/basis/coefficient.gz | Bin 308 -> 0 bytes test/QCIO_File/basis/exponent.gz | Bin 232 -> 0 bytes test/QCIO_File/basis/is_cartesian | 1 - test/QCIO_File/basis/num_contr | 1 - test/QCIO_File/basis/num_prim.gz | Bin 78 -> 0 bytes test/QCIO_File/basis/power.gz | Bin 92 -> 0 bytes test/QCIO_File/geometry/atomic_number.gz | Bin 69 -> 0 bytes test/QCIO_File/geometry/charge.gz | Bin 70 -> 0 bytes test/QCIO_File/geometry/coord.gz | Bin 125 -> 0 bytes test/QCIO_File/geometry/label.gz | Bin 59 -> 0 bytes test/QCIO_File/geometry/nuc_energy | 1 - test/QCIO_File/geometry/num_atom | 1 - test/QCIO_File/mo/classif.gz | Bin 59 -> 0 bytes test/QCIO_File/mo/eigenvalue.gz | Bin 209 -> 0 bytes test/QCIO_File/mo/label.gz | Bin 62 -> 0 bytes test/QCIO_File/mo/matrix.gz | Bin 2136 -> 0 bytes test/QCIO_File/mo/num_orb_sym.gz | Bin 66 -> 0 bytes test/QCIO_File/mo/symmetry.gz | Bin 71 -> 0 bytes test/QCIO_File/symmetry/num_sym | 1 - test/QCIO_File/system/energy | 1 - test/QCIO_File/system/num_alpha | 1 - test/QCIO_File/system/num_beta | 1 - test/QCIO_File/system/title | 1 - test/QCIO_File/wf/determinant.gz | Bin 64 -> 0 bytes test/c2h.out | 1264 ++++++++++++++++++++++ 45 files changed, 1926 insertions(+), 368 deletions(-) delete mode 100755 bin/eplf_input.py create mode 100755 bin/to_ezfio.py create mode 100644 eplf.config create mode 100644 src/det.irp.f create mode 100644 src/ezfio_interface.irp.f delete mode 100644 test/QCIO_File/ao/num_orb_sym.gz delete mode 100644 test/QCIO_File/ao/transformation.gz delete mode 100644 test/QCIO_File/basis/atom.gz delete mode 100644 test/QCIO_File/basis/coefficient.gz delete mode 100644 test/QCIO_File/basis/exponent.gz delete mode 100644 test/QCIO_File/basis/is_cartesian delete mode 100644 test/QCIO_File/basis/num_contr delete mode 100644 test/QCIO_File/basis/num_prim.gz delete mode 100644 test/QCIO_File/basis/power.gz delete mode 100644 test/QCIO_File/geometry/atomic_number.gz delete mode 100644 test/QCIO_File/geometry/charge.gz delete mode 100644 test/QCIO_File/geometry/coord.gz delete mode 100644 test/QCIO_File/geometry/label.gz delete mode 100644 test/QCIO_File/geometry/nuc_energy delete mode 100644 test/QCIO_File/geometry/num_atom delete mode 100644 test/QCIO_File/mo/classif.gz delete mode 100644 test/QCIO_File/mo/eigenvalue.gz delete mode 100644 test/QCIO_File/mo/label.gz delete mode 100644 test/QCIO_File/mo/matrix.gz delete mode 100644 test/QCIO_File/mo/num_orb_sym.gz delete mode 100644 test/QCIO_File/mo/symmetry.gz delete mode 100644 test/QCIO_File/symmetry/num_sym delete mode 100644 test/QCIO_File/system/energy delete mode 100644 test/QCIO_File/system/num_alpha delete mode 100644 test/QCIO_File/system/num_beta delete mode 100644 test/QCIO_File/system/title delete mode 100644 test/QCIO_File/wf/determinant.gz create mode 100644 test/c2h.out diff --git a/bin/eplf_input.py b/bin/eplf_input.py deleted file mode 100755 index f509785..0000000 --- a/bin/eplf_input.py +++ /dev/null @@ -1,26 +0,0 @@ -#!/usr/bin/env python - -import sys,os,time -from resultsFile import * - -# Check command line - -if len(sys.argv) != 2: - print "usage: "+sys.argv[0]+" file.out" - sys.exit(2) - -firstArg = sys.argv[1] - -file = open("eplf_grid.data","w") -print >>file, "&eplf_grid" -print >>file, " Npoints=80,80,80" -print >>file, "/" -file.close() - -file = getFile(firstArg) -print firstArg, 'recognized as', str(file).split('.')[-1].split()[0] - - -from qcio import qcio -write_qcioFile(file,"QCIO_File") - diff --git a/bin/to_ezfio.py b/bin/to_ezfio.py new file mode 100755 index 0000000..d30edf6 --- /dev/null +++ b/bin/to_ezfio.py @@ -0,0 +1,180 @@ +#!/usr/bin/env python + +import sys,os,time + +wd = os.path.dirname(__file__) +sys.path += [ wd+"../EZFIO/Python" ] +sys.path += [ "/home/scemama/resultsFile" ] + +from resultsFile import * + +# Check command line + +if len(sys.argv) == 2: + State=0 +elif len(sys.argv) == 3: + State=int(sys.argv[2]) +else: + print "usage: "+sys.argv[0]+" file.out" + sys.exit(2) + +firstArg = sys.argv[1] + +file = getFile(firstArg) +print firstArg, 'recognized as', str(file).split('.')[-1].split()[0] + +from ezfio import ezfio + +def write_ezfioFile(res,filename): + ezfio.set_file(filename) + +# Electrons + ezfio.set_electrons_elec_alpha_num(res.num_alpha) + ezfio.set_electrons_elec_beta_num(res.num_beta) + +# Nuclei + ezfio.set_nuclei_nucl_num(len(res.geometry)) + charge = [] + coord = [] + coord_x = [] + coord_y = [] + coord_z = [] + for a in res.geometry: + charge.append(a.charge) + if res.units == 'BOHR': + coord_x.append(a.coord[0]) + coord_y.append(a.coord[1]) + coord_z.append(a.coord[2]) + else: + coord_x.append(a.coord[0]/a0) + coord_y.append(a.coord[1]/a0) + coord_z.append(a.coord[2]/a0) + ezfio.set_nuclei_nucl_charge(charge) + ezfio.set_nuclei_nucl_coord(coord_x+coord_y+coord_z) + +# AO Basis + import string + is_cartesian = True + at = [] + num_prim = [] + magnetic_number = [] + angular_number = [] + power_x = [] + power_y = [] + power_z = [] + coefficient = [] + exponent = [] + for b in res.basis: + if '+' in b.sym or '-' in b.sym: + is_cartesian = False + names = ["s","p","d","f","g","h","i","j"] + for b in res.basis: + c = b.center + for i,atom in enumerate(res.geometry): + if atom.coord == c: + at.append(i+1) + num_prim.append(len(b.prim)) + if is_cartesian: + s = b.sym + power_x.append( string.count(s,"x") ) + power_y.append( string.count(s,"y") ) + power_z.append( string.count(s,"z") ) + else: + magnetic_number.append(names.index(b.sym[0])) + angular_number.append(int(b.sym[1:])) + coefficient.append( b.coef ) + exponent.append( [ p.expo for p in b.prim ] ) + if not is_cartesian: + print 'Only cartesian basis functions work...' + sys.exit(0) + ezfio.set_ao_basis_ao_num(len(res.basis)) + ezfio.set_ao_basis_ao_nucl(at) + ezfio.set_ao_basis_ao_prim_num(num_prim) + ezfio.set_ao_basis_ao_power(power_x+power_y+power_z) + prim_num_max = ezfio.get_ao_basis_ao_prim_num_max() + len_res_basis = len(res.basis) + for i in range(len(res.basis)): + coefficient[i] += [ 0. for j in range(len(coefficient[i]),prim_num_max) ] + exponent[i] += [ 0. for j in range(len(exponent[i]),prim_num_max) ] + coefficient = reduce(lambda x, y: x+y, coefficient, []) + exponent = reduce(lambda x, y: x+y, exponent, []) + coef = [] + expo = [] + for i in range(prim_num_max): + for j in range(i,len(coefficient),prim_num_max): + coef.append ( coefficient[j] ) + expo.append ( exponent[j] ) + ezfio.set_ao_basis_ao_coef(coef) + ezfio.set_ao_basis_ao_expo(expo) + +# MOs + NumOrbSym = [ s[1] for s in res.symmetries ] + mo_tot_num = sum(NumOrbSym) + ezfio.set_mo_basis_mo_tot_num(mo_tot_num) + if res.occ_num is not None: + ezfio.set_mo_basis_mo_occ(res.occ_num) + + MoTag = res.mo_types[-1] + mo = res.mo_sets[MoTag] + if len(mo) < mo_tot_num: + newmo = orbital() + newmo.eigenvalue = 0. + newmo.vector = [0. for i in range(mo_tot_num)] + while len(mo) < mo_tot_num: + mo.append(newmo) + Energies = [ m.eigenvalue for m in mo ] + ezfio.set_mo_basis_mo_energy(Energies) + + if res.occ_num is not None: + OccNum = res.occ_num[MoTag] + while len(OccNum) < mo_tot_num: + OccNum.append(0.) + ezfio.set_mo_basis_mo_occ(OccNum) + + cls = [ "v" for i in mo ] + for i in res.closed_mos: + cls[i] = 'c' + for i in res.active_mos: + cls[i] = 'a' + ezfio.set_mo_basis_mo_classif(cls) + + MoMatrix = [] + for m in mo: + for coef in m.vector: + MoMatrix.append(coef) + while len(MoMatrix) < len(mo[0].vector)**2: + MoMatrix.append(0.) + ezfio.set_mo_basis_mo_coef(MoMatrix) + +# Determinants + closed_mos = res.closed_mos + nactive = ezfio.get_mo_basis_mo_active_num() + dets_a = [] + dets_b = [] + for d in res.determinants: + dnew_a = [] + dnew_b = [] + for x in d['alpha']: + if x not in closed_mos: + dnew_a.append(x+1) + for x in d['beta']: + if x not in closed_mos: + dnew_b.append(x+1) + for x in range(nactive-len(dnew_b)): + dnew_b.append(0) + dets_a.append( dnew_a ) + dets_b.append( dnew_b ) + + coef = reduce(lambda x, y: x+y,res.det_coefficients,[]) + + if len(dets_a[0]) > 0: + ezfio.set_determinants_det_num(len(dets_a)) + ezfio.set_determinants_det_coef(coef) + ezfio.set_determinants_det_occ(dets_a+dets_b) + else: + ezfio.set_determinants_det_num(1) + ezfio.set_determinants_det_coef([1.]) + ezfio.set_determinants_det_occ(dets_a+dets_b) + + +write_ezfioFile(file,firstArg+".ezfio") diff --git a/eplf.config b/eplf.config new file mode 100644 index 0000000..1de93f5 --- /dev/null +++ b/eplf.config @@ -0,0 +1,40 @@ +electrons + elec_alpha_num integer + elec_beta_num integer + elec_num integer = electrons_elec_alpha_num + electrons_elec_beta_num + +nuclei + nucl_num integer + nucl_charge real (nuclei_nucl_num) + nucl_coord real (nuclei_nucl_num,3) + +determinants + det_num integer + det_occ integer (mo_basis_mo_active_num,determinants_det_num,2) + det_coef real (determinants_det_num) + +ao_basis + ao_num integer + ao_prim_num integer (ao_basis_ao_num) + ao_nucl integer (ao_basis_ao_num) + ao_power integer (ao_basis_ao_num,3) + ao_prim_num_max integer = maxval(ao_basis_ao_prim_num) + ao_coef real (ao_basis_ao_num,ao_basis_ao_prim_num_max) + ao_expo real (ao_basis_ao_num,ao_basis_ao_prim_num_max) + +mo_basis + mo_tot_num integer + mo_coef real (ao_basis_ao_num,mo_basis_mo_tot_num) + mo_classif character (mo_basis_mo_tot_num) + mo_closed_num integer =n_count_ch(mo_basis_mo_classif,size(mo_basis_mo_classif),'c') + mo_active_num integer =n_count_ch(mo_basis_mo_classif,size(mo_basis_mo_classif),'a') + mo_virtual_num integer =n_count_ch(mo_basis_mo_classif,size(mo_basis_mo_classif),'v') + mo_energy real (mo_basis_mo_tot_num) + mo_occ real (mo_basis_mo_tot_num) + +grid + point_num integer (3) + step_size real (3) + origin real (3) + opposite real (3) + diff --git a/src/Makefile b/src/Makefile index dfce639..75e4b3d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -15,7 +15,7 @@ FCFLAGS= -O3 SRC= OBJ= -LIB=-lqcio +LIB=../EZFIO/lib/libezfio.a include irpf90.make diff --git a/src/ao.irp.f b/src/ao.irp.f index 4846ba1..8bf774e 100644 --- a/src/ao.irp.f +++ b/src/ao.irp.f @@ -5,11 +5,11 @@ BEGIN_PROVIDER [ integer, ao_num ] ! Number of atomic orbitals END_DOC -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_basis_num_contr(ao_num) -!$OMP END CRITICAL (qcio_critical) - assert (ao_num > 0) + ao_num = -1 + call get_ao_basis_ao_num(ao_num) + if (ao_num <= 0) then + call abrt(irp_here,'Number of contracted gaussians should be > 0') + endif END_PROVIDER @@ -20,10 +20,16 @@ BEGIN_PROVIDER [ integer, ao_prim_num, (ao_num) ] ! Number of primitives per atomic orbital END_DOC -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_basis_num_prim(ao_prim_num) -!$OMP END CRITICAL (qcio_critical) + ao_prim_num = -1 + call get_ao_basis_ao_prim_num(ao_prim_num) + integer :: i + character*(80) :: message + do i=1,ao_num + if (ao_prim_num(i) <= 0) then + write(message,'(A,I6,A)') 'Number of primitives of contraction ',i,' should be > 0' + call abrt(irp_here,message) + endif + enddo END_PROVIDER @@ -34,10 +40,19 @@ BEGIN_PROVIDER [ integer, ao_nucl, (ao_num) ] ! Nucleus on which the atomic orbital is centered END_DOC -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_basis_atom(ao_nucl) -!$OMP END CRITICAL (qcio_critical) + ao_nucl = -1 + call get_ao_basis_ao_nucl(ao_nucl) + + character*(80) :: message + character*(30) :: range + write(range,'(A,I5,A)') '(1,',nucl_num,')' + integer :: i + do i=1,ao_num + if ( (ao_nucl(i) <= 0) .or. (ao_nucl(i) > nucl_num) ) then + write(message,'(A,I6,A)') 'Contraction ',i,' should be centered on a nucleus in the range'//trim(range) + call abrt(irp_here,message) + endif + enddo END_PROVIDER @@ -47,17 +62,17 @@ BEGIN_PROVIDER [ integer, ao_power, (ao_num,3) ] BEGIN_DOC ! x,y,z powers of the atomic orbital END_DOC - integer :: buffer(3,ao_num) + ao_power = 0 + call get_ao_basis_ao_power(ao_power) + + character*(80) :: message integer :: i,j - -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_basis_power(buffer) -!$OMP END CRITICAL (qcio_critical) - do i=1,3 do j=1,ao_num - ao_power(j,i) = buffer(i,j) + if (ao_power(j,i) < 0) then + write(message,'(A,I1,A,I6,A)') 'Power ',i,' of contraction ',j,' should be > 0' + call abrt(irp_here,message) + endif enddo enddo END_PROVIDER @@ -103,38 +118,42 @@ BEGIN_PROVIDER [ integer, ao_prim_num_max ] END_PROVIDER - BEGIN_PROVIDER [ real, ao_expo, (ao_prim_num_max,ao_num) ] -&BEGIN_PROVIDER [ real, ao_coef, (ao_prim_num_max,ao_num) ] + BEGIN_PROVIDER [ real, ao_expo, (ao_num,ao_prim_num_max) ] +&BEGIN_PROVIDER [ real, ao_coef, (ao_num,ao_prim_num_max) ] implicit none BEGIN_DOC ! Exponents and coefficients of the atomic orbitals END_DOC - double precision :: buffer(ao_prim_num_max,ao_num) - integer :: i,j + ao_expo = 0. + call get_ao_basis_ao_expo(ao_expo) -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_basis_exponent(buffer) -!$OMP END CRITICAL (qcio_critical) + integer :: i,j do i=1,ao_num do j=1,ao_prim_num(i) - ao_expo(j,i) = buffer(j,i) + if (ao_expo(i,j) <= 0.) then + character*(80) :: message + write(message,'(A,I6,A,I6,A)') 'Exponent ',j,' of contracted gaussian ',i,' is < 0' + call abrt(irp_here,message) + endif enddo enddo -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_basis_coefficient(buffer) -!$OMP END CRITICAL (qcio_critical) + ao_coef = 0. + call get_ao_basis_ao_coef(ao_coef) +! Normalization of the AO coefficients +! ------------------------------------ double precision :: norm, norm2 double precision :: goverlap + integer :: pow(3), l do i=1,ao_num do j=1,ao_prim_num(i) - norm = goverlap(ao_expo(j,i),ao_expo(j,i),ao_power(i,:)) - norm = sqrt(norm) - ao_coef(j,i) = buffer(j,i)/norm + pow(1) = ao_power(i,1) + pow(2) = ao_power(i,2) + pow(3) = ao_power(i,3) + norm = goverlap(ao_expo(i,j),ao_expo(i,j),pow) + ao_coef(i,j) = ao_coef(i,j)/sqrt(norm) enddo enddo diff --git a/src/ao_oneD_point.irp.f b/src/ao_oneD_point.irp.f index 5897810..f171b90 100644 --- a/src/ao_oneD_point.irp.f +++ b/src/ao_oneD_point.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [ real, ao_oneD_prim_p, (ao_prim_num_max,ao_num) ] +BEGIN_PROVIDER [ real, ao_oneD_prim_p, (ao_num,ao_prim_num_max) ] implicit none include 'types.F' @@ -9,22 +9,24 @@ BEGIN_PROVIDER [ real, ao_oneD_prim_p, (ao_prim_num_max,ao_num) ] real:: r2, rtemp ! Compute alpha*r or alpha*r^2 - do i=1,ao_num - r2 = point_nucl_dist_2(ao_nucl(i)) - do k=1,ao_prim_num(i) - ao_oneD_prim_p(k,i) = r2 + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_prim_p(i,k) = point_nucl_dist_2(ao_nucl(i)) enddo enddo ! Compute exp(-alpha*r) or exp(-alpha*r^2) - do i=1,ao_num - do k=1,ao_prim_num(i) - ao_oneD_prim_p(k,i) = exp(-ao_oneD_prim_p(k,i)*ao_expo(k,i)) + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_prim_p(i,k) = exp(-ao_oneD_prim_p(i,k)*ao_expo(i,k)) enddo - ! Cut below 1.d-12 - do k=1,ao_prim_num(i) - if ( abs(ao_oneD_prim_p(k,i)) < 1.e-12 ) then - ao_oneD_prim_p(k,i) = 0. + enddo + + ! Cut below 1.d-12 + do k=1,ao_prim_num_max + do i=1,ao_num + if ( abs(ao_oneD_prim_p(i,k)) < 1.e-12 ) then + ao_oneD_prim_p(i,k) = 0. endif enddo enddo @@ -43,15 +45,17 @@ BEGIN_PROVIDER [ real, ao_oneD_p, (ao_num) ] do i=1,ao_num ao_oneD_p(i) = 0. - do k=1,ao_prim_num(i) - ao_oneD_p(i) = ao_oneD_p(i) + ao_coef(k,i)*ao_oneD_prim_p(k,i) + enddo + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_p(i) = ao_oneD_p(i) + ao_coef(i,k)*ao_oneD_prim_p(i,k) enddo enddo END_PROVIDER -BEGIN_PROVIDER [ real, ao_oneD_prim_grad_p, (ao_prim_num_max,ao_num,3) ] +BEGIN_PROVIDER [ real, ao_oneD_prim_grad_p, (ao_num,ao_prim_num_max,3) ] implicit none include 'types.F' @@ -61,10 +65,10 @@ BEGIN_PROVIDER [ real, ao_oneD_prim_grad_p, (ao_prim_num_max,ao_num,3) ] integer :: i, k, l real:: factor do l=1,3 - do i=1,ao_num - factor = -2.*point_nucl_dist_vec(ao_nucl(i),l) - do k=1,ao_prim_num(i) - ao_oneD_prim_grad_p(k,i,l) = factor*ao_expo(k,i)*ao_oneD_prim_p(k,i) + do k=1,ao_prim_num_max + do i=1,ao_num + factor = -2.*point_nucl_dist_vec(ao_nucl(i),l) + ao_oneD_prim_grad_p(i,k,l) = factor*ao_expo(i,k)*ao_oneD_prim_p(i,k) enddo enddo enddo @@ -82,15 +86,17 @@ BEGIN_PROVIDER [ real, ao_oneD_grad_p, (ao_num,3) ] do l=1,3 do i=1,ao_num ao_oneD_grad_p(i,l) = 0. - do k=1,ao_prim_num(i) - ao_oneD_grad_p(i,l) = ao_oneD_grad_p(i,l) + ao_coef(k,i)*ao_oneD_prim_grad_p(k,i,l) - enddo + enddo + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_grad_p(i,l) = ao_oneD_grad_p(i,l) + ao_coef(i,k)*ao_oneD_prim_grad_p(i,k,l) + enddo enddo enddo END_PROVIDER -BEGIN_PROVIDER [ real, ao_oneD_prim_lapl_p, (ao_prim_num_max,ao_num) ] +BEGIN_PROVIDER [ real, ao_oneD_prim_lapl_p, (ao_num,ao_prim_num_max) ] implicit none include 'types.F' @@ -98,10 +104,10 @@ BEGIN_PROVIDER [ real, ao_oneD_prim_lapl_p, (ao_prim_num_max,ao_num) ] ! Laplacian of the one-dimensional component of the primitive AOs END_DOC integer :: i, k, l - do i=1,ao_num - do k=1,ao_prim_num(i) - ao_oneD_prim_lapl_p(k,i) = ao_oneD_prim_p(k,i) * ao_expo(k,i) * & - ( 4.*ao_expo(k,i)*point_nucl_dist_2(ao_nucl(i)) - 6. ) + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_prim_lapl_p(i,k) = ao_oneD_prim_p(i,k) * ao_expo(i,k) * & + ( 4.*ao_expo(i,k)*point_nucl_dist_2(ao_nucl(i)) - 6. ) enddo enddo @@ -119,8 +125,10 @@ BEGIN_PROVIDER [ real, ao_oneD_lapl_p, (ao_num) ] do i=1,ao_num ao_oneD_lapl_p(i) = 0. - do k=1,ao_prim_num(i) - ao_oneD_lapl_p(i) = ao_oneD_lapl_p(i) + ao_coef(k,i)*ao_oneD_prim_lapl_p(k,i) + enddo + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_lapl_p(i) = ao_oneD_lapl_p(i) + ao_coef(i,k)*ao_oneD_prim_lapl_p(i,k) enddo enddo diff --git a/src/debug_eplf.irp.f b/src/debug_eplf.irp.f index 4ecd631..1ab649b 100644 --- a/src/debug_eplf.irp.f +++ b/src/debug_eplf.irp.f @@ -4,10 +4,10 @@ program debug integer :: i,j integer :: k print *, '' - print *, 'Occupation numbers' - do k=1,mo_num - print *, k, mo_occ(k) - enddo +!print *, 'Occupation numbers' +!do k=1,mo_num +! print *, k, mo_occ(k) +!enddo read(*,*) i,j print *, '' diff --git a/src/det.irp.f b/src/det.irp.f new file mode 100644 index 0000000..49f658e --- /dev/null +++ b/src/det.irp.f @@ -0,0 +1,33 @@ +BEGIN_PROVIDER [ integer, det_num ] + implicit none + include "types.F" + BEGIN_DOC +! Number of determinants + END_DOC + + det_num = 1 + call get_determinants_det_num(det_num) + if (det_num < 1) then + call abrt(irp_here,'det_num should be > 0') + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ] +&BEGIN_PROVIDER [ real, det_coef, (det_num) ] + implicit none + + BEGIN_DOC +! det : Description of the active orbitals of the determinants +! det_coef : Determinant coefficients + END_DOC + + det = 0 + det_coef = 0. + call get_determinants_det_occ(det) + call get_determinants_det_coef(det_coef) + +END_PROVIDER + + + diff --git a/src/electrons.irp.f b/src/electrons.irp.f index d6cda6d..675c2da 100644 --- a/src/electrons.irp.f +++ b/src/electrons.irp.f @@ -1,48 +1,41 @@ BEGIN_PROVIDER [ integer, elec_alpha_num ] - - BEGIN_DOC + implicit none + BEGIN_DOC ! Number of alpha electrons - END_DOC - - implicit none - elec_alpha_num = -1 -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_system_num_alpha(elec_alpha_num) -!$OMP END CRITICAL (qcio_critical) - ASSERT (elec_alpha_num > 0) + END_DOC + elec_alpha_num = -1 + call get_electrons_elec_alpha_num(elec_alpha_num) + if (elec_alpha_num <= 0) then + call abrt(irp_here,'Number of alpha electrons should be > 0') + endif END_PROVIDER BEGIN_PROVIDER [ integer, elec_beta_num ] - - BEGIN_DOC + implicit none + BEGIN_DOC ! Number of beta electrons - END_DOC - - implicit none - elec_beta_num = -1 -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_system_num_beta(elec_beta_num) -!$OMP END CRITICAL (qcio_critical) - ASSERT (elec_beta_num >= 0) + END_DOC + elec_beta_num = 0 + call get_electrons_elec_beta_num(elec_beta_num) + if (elec_beta_num < 0) then + call abrt(irp_here,'Number of beta electrons should be >= 0') + endif END_PROVIDER BEGIN_PROVIDER [ integer, elec_num ] - - BEGIN_DOC + implicit none + BEGIN_DOC ! Number of electrons - END_DOC - implicit none + END_DOC - elec_num = elec_alpha_num + elec_beta_num - - ASSERT ( elec_num > 0 ) + elec_num = elec_alpha_num + elec_beta_num + ASSERT ( elec_num > 0 ) END_PROVIDER + BEGIN_PROVIDER [ integer, elec_num_2, (2) ] BEGIN_DOC diff --git a/src/eplf.irp.f b/src/eplf.irp.f index 4647dbb..c764287 100644 --- a/src/eplf.irp.f +++ b/src/eplf.irp.f @@ -1,5 +1,5 @@ program eplf - + provide mpi_rank call write_grid_eplf() call finish() end diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f index 630361b..7242946 100644 --- a/src/eplf_function.irp.f +++ b/src/eplf_function.irp.f @@ -5,7 +5,7 @@ BEGIN_PROVIDER [ real, eplf_gamma ] END_DOC real :: eps eps = -real(dlog(tiny(1.d0))) - eplf_gamma = density_p**(2./3.) * 100.*eps + eplf_gamma = density_p**(2./3.) * eps END_PROVIDER BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ] @@ -34,7 +34,7 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ] PROVIDE mo_coef do i=1,mo_num do j=i,mo_num - mo_eplf_integral_matrix(j,i) = 0. + mo_eplf_integral_matrix(j,i) = 0.d0 enddo @@ -72,8 +72,8 @@ END_PROVIDER double precision :: thr thr = 1.d-12 / eplf_gamma - eplf_up_up = 0. - eplf_up_dn = 0. + eplf_up_up = 0.d0 + eplf_up_dn = 0.d0 PROVIDE mo_coef_transp @@ -178,34 +178,34 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center) real :: gmma, center(3), c - ao_eplf_integral_numeric = 0. + ao_eplf_integral_numeric = 0.d0 do q=1,ao_prim_num(j) do p=1,ao_prim_num(i) - c = ao_coef(p,i)*ao_coef(q,j) + c = ao_coef(i,p)*ao_coef(j,q) integral = & ao_eplf_integral_primitive_oneD_numeric( & - ao_expo(p,i), & + ao_expo(i,p), & nucl_coord(ao_nucl(i),1), & ao_power(i,1), & - ao_expo(q,j), & + ao_expo(j,q), & nucl_coord(ao_nucl(j),1), & ao_power(j,1), & gmma, & center(1)) * & ao_eplf_integral_primitive_oneD_numeric( & - ao_expo(p,i), & + ao_expo(i,p), & nucl_coord(ao_nucl(i),2), & ao_power(i,2), & - ao_expo(q,j), & + ao_expo(j,q), & nucl_coord(ao_nucl(j),2), & ao_power(j,2), & gmma, & center(2)) * & ao_eplf_integral_primitive_oneD_numeric( & - ao_expo(p,i), & + ao_expo(i,p), & nucl_coord(ao_nucl(i),3), & ao_power(i,3), & - ao_expo(q,j), & + ao_expo(j,q), & nucl_coord(ao_nucl(j),3), & ao_power(j,3), & gmma, & @@ -216,60 +216,60 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center) end function -double precision function ao_eplf_integral_primitive_oneD2(a,xa,i,b,xb,j,gmma,xr) - implicit none - include 'constants.F' - - real, intent(in) :: a,b,gmma ! Exponents - real, intent(in) :: xa,xb,xr ! Centers - integer, intent(in) :: i,j ! Powers of xa and xb - integer :: ii, jj, kk, ll - real :: xp1,xp - real :: p1,p - double precision :: S(0:i+1,0:j+1) - double precision :: inv_p, di(max(i,j)), dj(j), c, c1 - - ASSERT (a>0.) - ASSERT (b>0.) - ASSERT (i>=0) - ASSERT (j>=0) - - ! Gaussian product - call gaussian_product(a,xa,b,xb,c1,p1,xp1) - call gaussian_product(p1,xp1,gmma,xr,c,p,xp) - inv_p = 1.d0/p - S(0,0) = dsqrt(pi*inv_p)*c*c1 - - ! Obara-Saika recursion - - do ii=1,max(i,j) - di(ii) = 0.5d0*inv_p*dble(ii) - enddo - - S(1,0) = (xp-xa) * S(0,0) - if (i>1) then - do ii=1,i-1 - S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0) - enddo - endif - - S(0,1) = (xp-xb) * S(0,0) - if (j>1) then - do jj=1,j-1 - S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1) - enddo - endif - - do jj=1,j - S(1,jj) = (xp-xa) * S(0,jj) + di(jj) * S(0,jj-1) - do ii=2,i - S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1) - enddo - enddo - - ao_eplf_integral_primitive_oneD2 = S(i,j) - -end function +!double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) +! implicit none +! include 'constants.F' +! +! real, intent(in) :: a,b,gmma ! Exponents +! real, intent(in) :: xa,xb,xr ! Centers +! integer, intent(in) :: i,j ! Powers of xa and xb +! integer :: ii, jj, kk, ll +! real :: xp1,xp +! real :: p1,p +! double precision :: S(0:i+1,0:j+1) +! double precision :: inv_p, di(max(i,j)), dj(j), c, c1 +! +! ASSERT (a>0.) +! ASSERT (b>0.) +! ASSERT (i>=0) +! ASSERT (j>=0) +! +! ! Gaussian product +! call gaussian_product(a,xa,b,xb,c1,p1,xp1) +! call gaussian_product(p1,xp1,gmma,xr,c,p,xp) +! inv_p = 1.d0/p +! S(0,0) = dsqrt(pi*inv_p)*c*c1 +! +! ! Obara-Saika recursion +! +! do ii=1,max(i,j) +! di(ii) = 0.5d0*inv_p*dble(ii) +! enddo +! +! S(1,0) = (xp-xa) * S(0,0) +! if (i>1) then +! do ii=1,i-1 +! S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0) +! enddo +! endif +! +! S(0,1) = (xp-xb) * S(0,0) +! if (j>1) then +! do jj=1,j-1 +! S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1) +! enddo +! endif +! +! do jj=1,j +! S(1,jj) = (xp-xa) * S(0,jj) + di(jj) * S(0,jj-1) +! do ii=2,i +! S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1) +! enddo +! enddo +! +! ao_eplf_integral_primitive_oneD = S(i,j) +! +!end function double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) implicit none @@ -351,39 +351,39 @@ double precision function ao_eplf_integral(i,j,gmma,center) do p=1,ao_prim_num(i) integral = & ao_eplf_integral_primitive_oneD( & - ao_expo(p,i), & + ao_expo(i,p), & nucl_coord(ao_nucl(i),1), & ao_power(i,1), & - ao_expo(q,j), & + ao_expo(j,q), & nucl_coord(ao_nucl(j),1), & ao_power(j,1), & gmma, & center(1)) * & ao_eplf_integral_primitive_oneD( & - ao_expo(p,i), & + ao_expo(i,p), & nucl_coord(ao_nucl(i),2), & ao_power(i,2), & - ao_expo(q,j), & + ao_expo(j,q), & nucl_coord(ao_nucl(j),2), & ao_power(j,2), & gmma, & center(2)) * & ao_eplf_integral_primitive_oneD( & - ao_expo(p,i), & + ao_expo(i,p), & nucl_coord(ao_nucl(i),3), & ao_power(i,3), & - ao_expo(q,j), & + ao_expo(j,q), & nucl_coord(ao_nucl(j),3), & ao_power(j,3), & gmma, & center(3)) - ao_eplf_integral = ao_eplf_integral + integral*ao_coef(p,i)*ao_coef(q,j) + ao_eplf_integral = ao_eplf_integral + integral*ao_coef(i,p)*ao_coef(j,q) enddo enddo end function - double precision function mo_eplf_integral(i,j) +double precision function mo_eplf_integral(i,j) implicit none integer :: i, j, k, l PROVIDE ao_eplf_integral_matrix @@ -399,5 +399,5 @@ end function endif enddo - end function +end function diff --git a/src/eplf_grid.irp.f b/src/eplf_grid.irp.f index 27173fe..db17b3c 100644 --- a/src/eplf_grid.irp.f +++ b/src/eplf_grid.irp.f @@ -1,10 +1,3 @@ -BEGIN_PROVIDER [ character*(100), grid_data_filename ] - BEGIN_DOC -! Name of the file containing the parameters of the grid - END_DOC - grid_data_filename = 'eplf_grid.data' -END_PROVIDER - BEGIN_PROVIDER [ character*(100), grid_cube_filename ] BEGIN_DOC ! Name of the file containing the parameters of the grid @@ -29,23 +22,16 @@ END_PROVIDER real :: origin(3) real :: opposite(3) - namelist /eplf_grid/ Npoints, step_size, origin, opposite - Npoints (:) = 80 origin (:) = UNDEFINED opposite (:) = UNDEFINED step_size(:) = UNDEFINED - - logical :: do_next - inquire(file=grid_data_filename,exist=do_next) - if (do_next) then - open(99,file=grid_data_filename,action='READ',status='OLD') - read(99,eplf_grid,end=90,err=80) - close(99) - endif - - 90 continue + call get_grid_point_num(Npoints) + call get_grid_origin(origin) + call get_grid_opposite(opposite) + call get_grid_step_size(step_size) + if (origin(1) == UNDEFINED) then integer :: i,l do l=1,3 @@ -81,12 +67,6 @@ END_PROVIDER grid_eplf_y_num = Npoints(2) grid_eplf_z_num = Npoints(3) - return - - 80 continue - print *, 'Error in input file' - stop 1 - END_PROVIDER BEGIN_PROVIDER [ real, grid_eplf, (grid_eplf_x_num,grid_eplf_y_num,grid_eplf_z_num) ] @@ -114,7 +94,7 @@ BEGIN_PROVIDER [ real, grid_eplf, (grid_eplf_x_num,grid_eplf_y_num,grid_eplf_z_n icount = mpi_size do iz=1,grid_eplf_z_num if (mpi_master) then - print *, iz + print *, int(100*dble(iz)/dble(grid_eplf_z_num)), '%' endif point(3) = grid_eplf_origin(3)+(iz-1)*grid_eplf_step(3) do iy=1,grid_eplf_y_num @@ -135,20 +115,19 @@ BEGIN_PROVIDER [ real, grid_eplf, (grid_eplf_x_num,grid_eplf_y_num,grid_eplf_z_n IRP_IF MPI integer :: dim, ierr - real :: buffer(grid_eplf_x_num*grid_eplf_y_num*grid_eplf_z_num) - icount = 0 do iz=1,grid_eplf_z_num + real :: buffer(grid_eplf_x_num*grid_eplf_y_num) + icount = 0 do iy=1,grid_eplf_y_num do ix=1,grid_eplf_x_num buffer(icount+ix) = grid_eplf(ix,iy,iz) enddo icount = icount + grid_eplf_x_num enddo - enddo - dim = grid_eplf_x_num * grid_eplf_y_num * grid_eplf_z_num - call MPI_REDUCE(buffer,grid_eplf,dim,mpi_real, & + dim = grid_eplf_x_num * grid_eplf_y_num + call MPI_REDUCE(buffer,grid_eplf(1,1,iz),dim,mpi_real, & mpi_sum,0,MPI_COMM_WORLD,ierr) - call MPI_BARRIER(MPI_COMM_WORLD,ierr) + enddo IRP_ENDIF END_PROVIDER diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f new file mode 100644 index 0000000..d0aa953 --- /dev/null +++ b/src/ezfio_interface.irp.f @@ -0,0 +1,122 @@ +BEGIN_SHELL [ /usr/bin/python ] + +data = [ \ +("nuclei_nucl_num" , "integer" , "" ), +("nuclei_nucl_charge" , "real" , "(nucl_num)" ), +("nuclei_nucl_coord" , "real" , "(nucl_num,3)" ), +("mo_basis_mo_coef" , "real" , "(ao_num,mo_tot_num)" ), +("mo_basis_mo_occ" , "real" , "(mo_tot_num)" ), +("electrons_elec_alpha_num" , "integer" , "" ), +("electrons_elec_beta_num" , "integer" , "" ), +("ao_basis_ao_num" , "integer" , "" ), +("ao_basis_ao_prim_num" , "integer" , "(ao_num)" ), +("ao_basis_ao_nucl" , "integer" , "(ao_num)" ), +("ao_basis_ao_power" , "integer" , "(ao_num,3)" ), +("ao_basis_ao_expo" , "real" , "(ao_num,ao_prim_num_max)" ), +("ao_basis_ao_coef" , "real" , "(ao_num,ao_prim_num_max)" ), +("determinants_det_num" , "integer" , "" ), +("determinants_det_coef" , "real" , "(det_num)" ), +("determinants_det_occ" , "integer" , "(elec_alpha_num-mo_closed_num,det_num,2)" ), +("grid_point_num" , "integer" , "(3)" ), +("grid_step_size" , "real" , "(3)" ), +("grid_origin" , "real" , "(3)" ), +("grid_opposite" , "real" , "(3)" ), +] + +data_no_set = [\ +("mo_basis_mo_tot_num" , "integer" , "" ), +("mo_basis_mo_active_num" , "integer" , "" ), +("mo_basis_mo_closed_num" , "integer" , "" ), +] + +def do_subst(t0,d): + t = t0 + t = t.replace("$X",d[0]) + t = t.replace("$T",d[1]) + t = t.replace("$D",d[2]) + if d[1].startswith("character"): + size = d[1].split("*")[1][1:-1] + u = "character" + elif d[1].startswith("double precision"): + u = d[1].replace(" ","_") + size = "1" + elif "*" in d[1]: + size = "1" + u = d[1].replace("*","") + else: + size = "1" + u = d[1] + t = t.replace("$U",u) + if d[2] == "": + t = t.replace("$S",size) + else: + if size == "1": + t = t.replace("$S","size(res)") + else: + t = t.replace("$S","%s*size(res)"%(size)) + print t +t0 = """ +subroutine get_$X(res) + implicit none + IRP_IF MPI + include 'mpif.h' + IRP_ENDIF + $T :: res$D + integer :: ierr + logical :: exists +!$OMP CRITICAL (ezfio_critical) + PROVIDE ezfio_filename + if (mpi_master) then + call ezfio_has_$X(exists) + if (exists) then + call ezfio_get_$X(res) + else + call ezfio_set_$X(res) + endif + endif + + IRP_IF MPI + call MPI_BCAST(res,$S,MPI_$U,0,MPI_COMM_WORLD,ierr) + if (ierr /= MPI_SUCCESS) then + call abrt(irp_here,'Unable to broadcast $X') + endif + IRP_ENDIF +!$OMP END CRITICAL (ezfio_critical) + +end +""" + + +t1 = """ +subroutine get_$X(res) + implicit none + IRP_IF MPI + include 'mpif.h' + IRP_ENDIF + $T :: res$D + integer :: ierr +!$OMP CRITICAL (ezfio_critical) + PROVIDE ezfio_filename + if (mpi_master) then + call ezfio_get_$X(res) + endif + + IRP_IF MPI + call MPI_BCAST(res,$S,MPI_$U,0,MPI_COMM_WORLD,ierr) + if (ierr /= MPI_SUCCESS) then + call abrt(irp_here,'Unable to broadcast $X') + endif + IRP_ENDIF +!$OMP END CRITICAL (ezfio_critical) + +end +""" + +for d in data: + do_subst(t0,d) + +for d in data_no_set: + do_subst(t1,d) + +END_SHELL + diff --git a/src/file.irp.f b/src/file.irp.f index e0572cc..eb1915a 100644 --- a/src/file.irp.f +++ b/src/file.irp.f @@ -1,24 +1,31 @@ -BEGIN_PROVIDER [ character*(128), qcio_filename ] - implicit none - BEGIN_DOC -! Name of the QCIO file +BEGIN_PROVIDER [ character*(64), ezfio_filename ] + BEGIN_DOC +! Name of the ezfio file END_DOC - integer :: iargc + IRP_IF MPI + include 'mpif.h' + integer :: ierr + IRP_ENDIF + + if (mpi_master) then + call getarg(1,ezfio_filename) + if (ezfio_filename == '') then + call ezfio_get_filename(ezfio_filename) + endif + endif IRP_IF MPI - PROVIDE mpi_rank + call MPI_BCAST(ezfio_filename,64,MPI_character,0,MPI_COMM_WORLD,ierr) + if (ierr /= MPI_SUCCESS) then + call abrt(irp_here,'Unable to broadcast ezfio_filename') + endif IRP_ENDIF - call getarg(1,qcio_filename) - if (qcio_filename /= '') then - call qcio_set_file(qcio_filename) - else - call qcio_get_filename(qcio_filename) - endif + + call ezfio_set_file(ezfio_filename) if (.not.mpi_master) then - call qcio_set_read_only(.True.) + call ezfio_set_read_only(.True.) endif call barrier END_PROVIDER - diff --git a/src/mo.irp.f b/src/mo.irp.f index 07639b8..a83c022 100644 --- a/src/mo.irp.f +++ b/src/mo.irp.f @@ -1,59 +1,39 @@ BEGIN_PROVIDER [ integer, mo_closed_num ] - implicit none - + implicit none BEGIN_DOC ! Number of closed shell Molecular orbitals END_DOC -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_mo_num_closed(mo_closed_num) -!$OMP END CRITICAL (qcio_critical) - ASSERT (mo_closed_num >= 0) - ASSERT (mo_closed_num <= elec_alpha_num) - ASSERT (mo_closed_num <= elec_beta_num) + mo_closed_num = -1 + call get_mo_basis_mo_closed_num(mo_closed_num) + if (mo_closed_num < 0) then + call abrt(irp_here,'Number of closed-shell MOs should be >=0') + endif + if (mo_closed_num > elec_alpha_num) then + call abrt(irp_here,'Number of closed-should MOs should be <= Number of alpha electrons') + endif + if (mo_closed_num > elec_beta_num) then + call abrt(irp_here,'Number of closed-should MOs should be <= Number of beta electrons') + endif END_PROVIDER BEGIN_PROVIDER [ integer, mo_active_num ] - implicit none - + implicit none BEGIN_DOC ! Number of active Molecular orbitals END_DOC -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_mo_num_active(mo_active_num) -!$OMP END CRITICAL (qcio_critical) - ASSERT (mo_active_num >= 0) - -END_PROVIDER - -BEGIN_PROVIDER [ real, mo_occ, (mo_num) ] - implicit none - - BEGIN_DOC -! Occupation numbers of molecular orbitals - END_DOC - - double precision, allocatable :: buffer(:) - allocate ( buffer(mo_tot_num) ) -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_mo_occupation(buffer) -!$OMP END CRITICAL (qcio_critical) - integer :: i - do i=1,mo_num - mo_occ(i) = buffer(i) - enddo - deallocate(buffer) + mo_active_num = -1 + call get_mo_basis_mo_active_num(mo_active_num) + if (mo_active_num < 0) then + call abrt(irp_here,'Number of active MOs should be >=0') + endif END_PROVIDER BEGIN_PROVIDER [ integer, mo_num ] - implicit none - + implicit none BEGIN_DOC ! Number of Molecular orbitals END_DOC @@ -63,42 +43,29 @@ BEGIN_PROVIDER [ integer, mo_num ] END_PROVIDER -BEGIN_PROVIDER [ real, mo_coef, (ao_num,mo_num) ] - implicit none - +BEGIN_PROVIDER [ real, mo_occ, (mo_tot_num) ] + implicit none BEGIN_DOC -! Molecular orbital coefficients +! Occupation numbers of molecular orbitals END_DOC - integer :: i, j - double precision, allocatable :: buffer(:,:) - allocate (buffer(ao_num,mo_tot_num)) - -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_mo_matrix(buffer) -!$OMP END CRITICAL (qcio_critical) - do j=1,mo_num - do i=1,ao_num - mo_coef(i,j) = buffer(i,j) - enddo - enddo - deallocate (buffer) + call get_mo_basis_mo_occ(mo_occ) END_PROVIDER -BEGIN_PROVIDER [ logical, mo_non_zero, (mo_num,ao_num) ] +BEGIN_PROVIDER [ real, mo_coef, (ao_num,mo_num) ] implicit none - BEGIN_DOC - ! Values where the mo coefficients are /= 0. +! Molecular orbital coefficients END_DOC + integer :: i, j + real :: buffer(ao_num,mo_tot_num) - integer :: i, j, k - - do j=1,ao_num - do k=1,mo_num - mo_non_zero(k,j) = mo_coef_transp(j,k) /= 0. + buffer = 0. + call get_mo_basis_mo_coef(buffer) + do i=1,mo_num + do j=1,ao_num + mo_coef(j,i) = buffer(j,i) enddo enddo @@ -130,39 +97,36 @@ BEGIN_PROVIDER [ logical, mo_is_closed, (mo_num) ] ! mo_is_active : True if mo(i) is an active orbital END_DOC - character, allocatable :: buffer(:) - allocate (buffer(mo_num)) -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_mo_classif(buffer) -!$OMP END CRITICAL - + character :: buffer(mo_tot_num) + integer :: i do i=1,mo_num if ( buffer(i) == 'c' ) then - mo_is_closed(i) = .True. + mo_is_closed(i) = .True. mo_is_active(i) = .False. else if ( buffer(i) == 'a' ) then - mo_is_closed(i) = .False. + mo_is_closed(i) = .False. mo_is_active(i) = .True. else mo_is_closed(i) = .False. mo_is_active(i) = .False. endif enddo - deallocate (buffer) END_PROVIDER + BEGIN_PROVIDER [ integer, mo_tot_num ] BEGIN_DOC -! Total number of MOs in the QCIO file +! Total number of MOs in the EZFIO file END_DOC -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_mo_num_orb_tot(mo_tot_num) -!$OMP END CRITICAL (qcio_critical) - ASSERT (mo_tot_num > 0) + + mo_tot_num = -1 + call get_mo_basis_mo_tot_num(mo_tot_num) + if (mo_tot_num <= 0) then + call abrt(irp_here,'Tota number of MOs can''t be <0') + endif END_PROVIDER + diff --git a/src/mpi.irp.f b/src/mpi.irp.f index a7722d2..bc0dca5 100644 --- a/src/mpi.irp.f +++ b/src/mpi.irp.f @@ -62,8 +62,6 @@ BEGIN_PROVIDER [ integer, mpi_size ] IRP_ENDIF - call iinfo(irp_here,'mpi_size',mpi_size) - END_PROVIDER diff --git a/src/nuclei.irp.f b/src/nuclei.irp.f index b47050e..a00dd85 100644 --- a/src/nuclei.irp.f +++ b/src/nuclei.irp.f @@ -5,11 +5,11 @@ BEGIN_PROVIDER [ integer, nucl_num ] ! Number of nuclei END_DOC -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_geometry_num_atom(nucl_num) -!$OMP END CRITICAL (qcio_critical) - assert (nucl_num > 0) + nucl_num = -1 + call get_nuclei_nucl_num(nucl_num) + if (nucl_num <= 0) then + call abrt(irp_here,'Number of nuclei should be > 0') + endif END_PROVIDER @@ -20,19 +20,16 @@ BEGIN_PROVIDER [ real, nucl_charge, (nucl_num) ] ! Nuclear charge END_DOC - double precision,allocatable :: buffer(:) - allocate(buffer(nucl_num)) -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_geometry_charge(buffer) -!$OMP END CRITICAL (qcio_critical) + nucl_charge = -1.d0 + call get_nuclei_nucl_charge(nucl_charge) integer :: i do i=1,nucl_num - nucl_charge(i) = buffer(i) + if (nucl_charge(i) <= 0.) then + call abrt(irp_here,'Nuclear charges should be > 0') + endif enddo - deallocate(buffer) -END_PROVIDER +END_PROVIDER BEGIN_PROVIDER [ real, nucl_coord, (nucl_num,3) ] implicit none @@ -40,24 +37,13 @@ BEGIN_PROVIDER [ real, nucl_coord, (nucl_num,3) ] BEGIN_DOC ! Nuclear coordinates END_DOC - double precision, allocatable :: buffer(:,:) - allocate (buffer(3,nucl_num)) -!$OMP CRITICAL (qcio_critical) - PROVIDE qcio_filename - call qcio_get_geometry_coord(buffer) -!$OMP END CRITICAL (qcio_critical) - - integer :: i,j - do i=1,3 - do j=1,nucl_num - nucl_coord(j,i) = buffer(i,j) - enddo - enddo - deallocate(buffer) + nucl_coord = 0. + call get_nuclei_nucl_coord(nucl_coord) END_PROVIDER + BEGIN_PROVIDER [ real, nucl_dist_2, (nucl_num,nucl_num) ] &BEGIN_PROVIDER [ real, nucl_dist_vec, (nucl_num,nucl_num,3) ] implicit none diff --git a/test/QCIO_File/ao/num_orb_sym.gz b/test/QCIO_File/ao/num_orb_sym.gz deleted file mode 100644 index 21163f4198eba0b8764df34f4926142dc0f109f1..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 63 zcmb2|=HQT-Y3;@IKhW9JKi(}fCsjW&Uq7!jH$J~8DZaQemmxU7;7QN{7B#L2mtNnE Snx+bl0u0RIn$~#?3=9B~RucFC diff --git a/test/QCIO_File/ao/transformation.gz b/test/QCIO_File/ao/transformation.gz deleted file mode 100644 index 1239a3fe295e6ade39cb0888b46ee63f2063f7e2..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 201 zcmb2|=HQT-Y3;@IKhW9JKi(}fCsjW&U%#X%F|Rl+zbH4cBr`vc;qB#(ybT6Ct_SC7 z-ceR`o+5aoe(8x+msKJm&yA+P^45Q3`KWm9@lPA~-3veWecnt~PKF=MCL9fHj29Mo zGF(t#kTT(9kYaIY-p1t6%#d**j3GmSftQH^tOP8<0OT%kWB}<Kd~e~mmxU7;7QN{hGWaPB3w3V$%?9Y YUC0QT5aXCUi-U>b^@_ampG*u4035Ivs{jB1 diff --git a/test/QCIO_File/basis/coefficient.gz b/test/QCIO_File/basis/coefficient.gz deleted file mode 100644 index f09196f57bffc63a14e7115056edb5b71f45eaa4..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 308 zcmV-40n7d$iwFoDnLA1X|4~CpPhUo9Y-KNEVRLD7FJo_IW@c$)X=QG70PU7LZo@DP zK=+)&dtie9&jE7(D^&`Tx^g?TOOXKNsApP_l3{Y&__BSE_+d?Kn~yq_7`adg-s85T zPGBHPhOBY`nSdck+=G-p!F^W(q$-xQBZD3pW|%H+D^M^R>thQ#XJnXLsOAbBIE92M zk5|P2yMd-gjW~c!ATWxPK>wmU>(~h>D#a^(sbc~xtjyEq1Biw)tJi8ZGbk*zo|Rt_ zA4jsr6W~h%1FtAHvRbcJV4OJy6P<1_z)`|ODOS+B@2>l1-)`%+ZtJ#g>;5n30kQ|k zp2%9a|GTZ*x~HmJDeW>0`_@_9LX9%|;;{%L}kB4*5Ft-RC*Zqf8C6#@~vu`C-y+wcnUA`PRmF#il z%oSYB&_e|a8?uy;#?n3uWP^k_8y3mC#bTL3;fW$-wT}))qMH7ZDbfJ6jA}rKdxv;CL-stZdsX-UIJiP%6(B73C7ytmSvT5xA diff --git a/test/QCIO_File/basis/is_cartesian b/test/QCIO_File/basis/is_cartesian deleted file mode 100644 index 62a6e3c..0000000 --- a/test/QCIO_File/basis/is_cartesian +++ /dev/null @@ -1 +0,0 @@ -T diff --git a/test/QCIO_File/basis/num_contr b/test/QCIO_File/basis/num_contr deleted file mode 100644 index 7493501..0000000 --- a/test/QCIO_File/basis/num_contr +++ /dev/null @@ -1 +0,0 @@ - 35 diff --git a/test/QCIO_File/basis/num_prim.gz b/test/QCIO_File/basis/num_prim.gz deleted file mode 100644 index 4921c03f38fbf9457c464f32ac0c4a0d9dbb790d..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 78 zcmV-U0I~lciwFoDnLA1X|4~CpPhUo9Y-KNEVRLD7FK%^hUvP41Z2(hHFyvCe0LG>Y kSVR@fi5A3WsDd$Nf`&s`a5U)17Nm`j0Jfunzxe_H0A9=da|ou0YO2mk;yD|<46gR*!-{+se4@gqD ze!}*;VT74j0Z^&L83t&&<_ZY0t2JOGn%y5x$gBo}%pEI`7IKbdZ&;f>H%eZ3t+MXyW3lte0 SY3tdbz_9(B{qbf71_l75`xc`B diff --git a/test/QCIO_File/mo/matrix.gz b/test/QCIO_File/mo/matrix.gz deleted file mode 100644 index 5c7d316a60f729a50568fe51f55943f70af60962..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 2136 zcmV-e2&eZSiwFoDnLA1X|4~CpPhUo9Y-KNPZ!c|ObaH8U0PP(~aw9qJKCke7LJR^R z7CxZ=|D^#`RbizKdcqRevbs~pLYbMwzUSK;{^#w}kN7S9$M~<(3}(Oi{dMWrGq54_ z4*>cFFiHWv_2cia|MvbRz}n~^j~e$Kg{{x*tuDNH#uDI15f@!>37BmT`tS@005qPM z(Bndls?ZJ%xa_FjQ?z3n_+5a=X0%_Rk)_X&rbV21vDHyQ^TXZWX|9^d;6<8L$wwg) zrOo5d-->@xFxL!gfHc~D!FXNN47)CBBV`6>>MHq=L_YYwdr-3+SG0P@5`e3tO3^MW z^;*T$m(eiS?JiL?EfFGtT*W?&=>^Bi#aRLT=6H~T9W=Lw!nsw`u8H>X z3z7{f(r_IFOOBE@L^Jkr%eivZa)5LDoV)#Yr@0JG?FZ`N@%N*I+&!nw>@ym4fRfh@ zd#*J%|2qI4M~QC6?pJ6RoJD9K!tyF;cj07~?RHQBo<-VI22UA0W$+Fed7Bc-oG8l&YTF&MH-t{rNSoHWFVsiajeC^11PNz zt^pX%-R5LJIucViYBoNUKj57n(((4a?ti$thg=luvj=f$(IV(!Iofl8qh#*(+uctT z>^aS4%58c3ksuB&<*V5UQoeGx7Sgr_$F3TTB5P2@)KN<%xs9>@{DrojKVt<5Ghw?$ zE%NQq(^87&yvu&Uqs4iE;f$>c?e2cuD&F@C%BqmN2xY+0F0j|6Mwwi*&=UL}dYglI zNYudcx*DW><>9Ly&GWleQ^cAAD7MiLO;C@dMck-S8lN-u{DKPrIWyQ^n)}Ugl&r#YuJVfU|~hN=B1*Z8_u zWqi>NgT2V*blnd4I?)6~Uxqvl?POG=GNMwAuTx{0O%56yw)iDeo>g!=CF+!@Q=#61b6)`0y^2R2WzyzE3N@(eL;gEloKpr*89ZgsqoO$^uL9(g7J+t{ zx8;dP!VDcnMeotNb9IV#Wv0|mTmt4$+v@?iCzH?oL^*R+b}}tNdyIw-s%a2VSA8SH zRZQ1pwW1}!6*AlI_KVSgWTfu7b;cs>a-T2LFZx3P(vV-M*ZH4=_5oxk+J8KNolrB$ zGh9}H6NxowXC67)bHKkrWl%a&8MqoYQ(Rd$D(MzMk0zvM3OjnXg=g+9&3vOn4v4oJ zv^j<3*u^Y|9LePg+3Tu1+t^Dokm}Enn$7`i?9qg!fl0bTWC@5dkK^?Fmc(y2Aj#8) zinAt&n^T|FzY9?wD5;WNJQ$@3uc8wpfzA9Z-^b$*keY`_DdlN%clU&3G+;B@qq|Mv zl&DjpPKo+)5`~R=nt^m%>;B5;oXo{qSD03FkKUmq&k!j+IC_2WQxKHSL9PJT=Eym= zQqUbp$n$Vf;WeVHHNss^EgF<$g(@&l;%kfLF1`PmI+Y{7QYu9cXV7&^5L8*`I`xQ=(3ZIwk6>!U$oO znYhl3W&;)R7e$j6kplq2qNDPW%Gl2^@SA|+O4&VA*d%Gy5}G0{MnmaPaeHX2z-ne+ zuhZlRlCdaKHXgAAyj`2?H`lg3+D`Yk$2?BaV0!cv32ioa*aVAk73)1-XYXLNy&=u{ zl!wPMG_C1d3IoZEFJ#rG@pj>qPkYaQFiK7tJZ12d!D~UGX9b>X%RihKLxbFRcy?m8^j^TGd9BS%Q$EOlzN#*C_-jY~NV^*-wNUD9& zVJaCCmw@X`HSCp-*rq#eq;KUhL(RaE%i!!|QM()Yb}_}SsisC0kRqiIR@DP3ly#oP zSTkg;b!$LGkh?_#d>p2wGR&+C8;oT^;)xsCCs*QZk=f?_K__BBSa0ijvU7Vypc73MSd61 zEK>mf7g)3u)NE?S(!=<2R{pU)j=|g9L(Ty)t3N+}KdN5CbriMC6tQWtdHB8$QLyO8 OQ09LmviXJ3cK`q!{UoIT diff --git a/test/QCIO_File/mo/num_orb_sym.gz b/test/QCIO_File/mo/num_orb_sym.gz deleted file mode 100644 index c7d6bcf3612eb66d78bfd896fe59984325fae30d..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 66 zcmb2|=HQT-Y3;@IKhW9JKi(}fCsjW;Uq7!jH$J~8DZaQemmxU7;7QN{7B#L27hc_s Wnx==;GWuc#7-~1jl;<%pFaQACViiCD diff --git a/test/QCIO_File/mo/symmetry.gz b/test/QCIO_File/mo/symmetry.gz deleted file mode 100644 index fa8f58a846124b2af6e734f5ddcf23064f53c068..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 71 zcmb2|=HQT-Y3;@IKhW9JKi(}fCsjW;U%$9AH#fDUsFEQ#z~D*H0fuACxFTFOYRQVq b7&s~~V(n$0zu+VTHxq;Wi=VzfnHU%VT~`= 49.3885 + Charge= 0.0000 electrons + Dipole moment (field-independent basis, Debye): + X= 0.0000 Y= 0.0000 Z= -0.7447 Tot= 0.7447 + Quadrupole moment (field-independent basis, Debye-Ang): + XX= -12.8334 YY= -12.8334 ZZ= -7.9380 + XY= 0.0000 XZ= 0.0000 YZ= 0.0000 + Traceless Quadrupole moment (field-independent basis, Debye-Ang): + XX= -1.6318 YY= -1.6318 ZZ= 3.2636 + XY= 0.0000 XZ= 0.0000 YZ= 0.0000 + Octapole moment (field-independent basis, Debye-Ang**2): + XXX= 0.0000 YYY= 0.0000 ZZZ= -4.5135 XYY= 0.0000 + XXY= 0.0000 XXZ= -0.5099 XZZ= 0.0000 YZZ= 0.0000 + YYZ= -0.5099 XYZ= 0.0000 + Hexadecapole moment (field-independent basis, Debye-Ang**3): + XXXX= -14.6424 YYYY= -14.6424 ZZZZ= -35.5376 XXXY= 0.0000 + XXXZ= 0.0000 YYYX= 0.0000 YYYZ= 0.0000 ZZZX= 0.0000 + ZZZY= 0.0000 XXYY= -4.8808 XXZZ= -9.7271 YYZZ= -9.7271 + XXYZ= 0.0000 YYXZ= 0.0000 ZZXY= 0.0000 + N-N= 2.032531471843D+01 E-N=-2.179436435543D+02 KE= 7.596465821717D+01 + Symmetry A1 KE= 7.166421322913D+01 + Symmetry A2 KE= 1.251847281254D-51 + Symmetry B1 KE= 2.150222494019D+00 + Symmetry B2 KE= 2.150222494019D+00 + Orbital energies and kinetic energies (alpha): + 1 2 + 1 (SG)--O -11.28905 16.01354 + 2 (SG)--O -11.26180 16.01327 + 3 (SG)--O -1.02950 1.83829 + 4 (SG)--O -0.73153 1.22284 + 5 (PI)--O -0.41971 1.07511 + 6 (PI)--O -0.41971 1.07511 + 7 (SG)--O -0.18326 1.48833 + 8 (PI)--V 0.19537 1.00189 + 9 (PI)--V 0.19537 1.00189 + 10 (SG)--V 0.20534 0.59276 + 11 (SG)--V 0.42376 1.07678 + 12 (SG)--V 0.51778 1.13419 + 13 (PI)--V 0.64551 1.66542 + 14 (PI)--V 0.64551 1.66542 + 15 (SG)--V 0.75277 1.63015 + 16 (PI)--V 0.80367 2.26439 + 17 (PI)--V 0.80367 2.26439 + 18 (SG)--V 0.90360 1.65967 + 19 (DLTA)--V 1.23588 1.77278 + 20 (DLTA)--V 1.23588 1.77278 + 21 (SG)--V 1.26247 2.33402 + 22 (PI)--V 1.29170 1.95698 + 23 (PI)--V 1.29170 1.95698 + 24 (SG)--V 1.32557 2.50812 + 25 (DLTA)--V 1.66864 2.17626 + 26 (DLTA)--V 1.66864 2.17626 + 27 (PI)--V 1.81433 2.48590 + 28 (PI)--V 1.81433 2.48590 + 29 (SG)--V 2.00053 3.04043 + 30 (PI)--V 2.46174 3.22061 + 31 (PI)--V 2.46174 3.22061 + 32 (SG)--V 2.57219 5.01846 + 33 (SG)--V 2.63829 5.52949 + 34 (SG)--V 3.03059 6.90499 + 35 (SG)--V 4.61072 7.92694 + Total kinetic energy from orbitals= 7.745298861090D+01 + Isotropic Fermi Contact Couplings + Atom a.u. MegaHertz Gauss 10(-4) cm-1 + 1 C(13) 0.16226 182.41462 65.09007 60.84697 + 2 C(13) 0.83446 938.08691 334.73272 312.91211 + 3 H(1) 0.00270 12.06532 4.30521 4.02456 + -------------------------------------------------------- + Center ---- Spin Dipole Couplings ---- + 3XX-RR 3YY-RR 3ZZ-RR + -------------------------------------------------------- + 1 Atom -0.077447 -0.077447 0.154894 + 2 Atom -0.421656 -0.421656 0.843311 + 3 Atom -0.010356 -0.010356 0.020711 + -------------------------------------------------------- + XY XZ YZ + -------------------------------------------------------- + 1 Atom 0.000000 0.000000 0.000000 + 2 Atom 0.000000 0.000000 0.000000 + 3 Atom 0.000000 0.000000 0.000000 + -------------------------------------------------------- + + + --------------------------------------------------------------------------------- + Anisotropic Spin Dipole Couplings in Principal Axis System + --------------------------------------------------------------------------------- + + Atom a.u. MegaHertz Gauss 10(-4) cm-1 Axes + + Baa -0.0774 -10.393 -3.708 -3.467 1.0000 0.0000 0.0000 + 1 C(13) Bbb -0.0774 -10.393 -3.708 -3.467 0.0000 1.0000 0.0000 + Bcc 0.1549 20.785 7.417 6.933 0.0000 0.0000 1.0000 + + Baa -0.4217 -56.582 -20.190 -18.874 1.0000 0.0000 0.0000 + 2 C(13) Bbb -0.4217 -56.582 -20.190 -18.874 0.0000 1.0000 0.0000 + Bcc 0.8433 113.164 40.380 37.748 0.0000 0.0000 1.0000 + + Baa -0.0104 -5.525 -1.972 -1.843 0.0000 1.0000 0.0000 + 3 H(1) Bbb -0.0104 -5.525 -1.972 -1.843 1.0000 0.0000 0.0000 + Bcc 0.0207 11.051 3.943 3.686 0.0000 0.0000 1.0000 + + + --------------------------------------------------------------------------------- + + No NMR shielding tensors so no spin-rotation constants. + Leave Link 601 at Sat Mar 28 08:02:19 2009, MaxMem= 104857600 cpu: 0.3 + (Enter /usr/local/g03/l9999.exe) + 1\1\GINC-LPQSV11\SP\ROHF\CC-pVDZ\C2H1(2)\SCEMAMA\28-Mar-2009\0\\#P CC- + PVDZ ROHF GFPRINT POP=FULL 6D 10F\\titre\\0,2\C,0,0.,0.,0.59801\C,0,0. + ,0.,-0.59801\H,0,0.,0.,1.65962\\Version=AM64L-G03RevC.02\State=2-SG\HF + =-76.1419687\RMSD=8.827e-05\Dipole=0.,0.,0.2930001\PG=C*V [C*(H1C1C1)] + \\@ + + + IT WAS AN ACT OF DESPARATION. FOR SIX YEARS I HAD STRUGGLED WITH THE + BLACKBODY THEORY. I KNEW THE PROBLEM WAS FUNDAMENTAL, AND I KNEW + THE ANSWER. I HAD TO FIND A THEORETICAL EXPLANATION AT ANY COST, + EXCEPT FOR THE INVIOLABLITY OF THE TWO LAWS OF THERMODYNAMICS. + -- MAX PLANCK, 1931 + Job cpu time: 0 days 0 hours 0 minutes 9.0 seconds. + File lengths (MBytes): RWF= 13 Int= 0 D2E= 0 Chk= 11 Scr= 1 + Normal termination of Gaussian 03 at Sat Mar 28 08:02:20 2009.