From d6cb6538e77f4fab7a790e751a919da2544558b7 Mon Sep 17 00:00:00 2001 From: Ravindra Shinde Date: Thu, 4 Mar 2021 12:14:57 +0100 Subject: [PATCH] failsafe mechanism added; new keywords added; beginning of the inline blocks --- src/build.sh | 2 +- src/iochamp.f90 | 164 +++++------------- src/keywords_m.F90 | 10 ++ ...iodic_table_m.f90 => periodic_table_m.F90} | 38 ++-- src/test-champ.inp | 27 ++- src/test-champ.out | 6 +- 6 files changed, 103 insertions(+), 144 deletions(-) rename src/{periodic_table_m.f90 => periodic_table_m.F90} (76%) diff --git a/src/build.sh b/src/build.sh index c41ab31..37771b0 100755 --- a/src/build.sh +++ b/src/build.sh @@ -1,5 +1,5 @@ #!/bin/bash -ifort -c periodic_table_m.f90 +ifort -c periodic_table_m.F90 ifort -c keywords_m.F90 ifort iochamp.f90 keywords_m.o periodic_table_m.o /usr/local/lib/libfdf.a diff --git a/src/iochamp.f90 b/src/iochamp.f90 index 7889abe..5721ceb 100644 --- a/src/iochamp.f90 +++ b/src/iochamp.f90 @@ -5,7 +5,9 @@ PROGRAM iochamp USE fdf USE prec USE parse - use io_fdf + USE io_fdf + USE utils + ! Note the following two modules are being used to store and process the parsed data use keywords @@ -33,6 +35,7 @@ PROGRAM iochamp !type(fdf_file) :: fdffile integer :: max_iteration, max_iter, linecount, argument(5) real(dp) :: float_value + character(len=1) :: char1 character(len=20) :: real_format = '(A, T20, F14.8)' character(len=20) :: int_format = '(A, T20, I8)' character(len=80) :: string_format = '(A, T40, A)' @@ -176,8 +179,18 @@ PROGRAM iochamp write(6,*) '------------------------------------------------------' + +! Additional keywords. check if they clash with existing + + excess_charge = fdf_integer('excess_charge', 0) + write(6,fmt=int_format) 'Excess charges =', excess_charge + + multiplicity = fdf_integer('multiplicity', 1) ! default multiplicity singlet. An assertion is needed + write(6,fmt=int_format) 'multiplicity =', multiplicity + + ! write(6,'(A,4X)') 'optimize_wavefunction using bline', (subblock(i), i = 1, 4) if (fdf_block('general', bfdf)) then @@ -203,7 +216,7 @@ PROGRAM iochamp ia = 1 open (unit=12,file=file_molecule, iostat=iostat, action='read' ) - if (iostat .ne. 0) error stop "Problem in operning the molecule file" + if (iostat .ne. 0) error stop "Problem in opening the molecule file" read(12,*) natoms print*, "natoms ", natoms if (.not. allocated(cent)) allocate(cent(3,natoms)) @@ -237,8 +250,8 @@ PROGRAM iochamp ! write(*,*) "linecount", fdf_block_linecount("molecule") do while((fdf_bline(bfdf, pline))) - - if (pline%ntokens == 1) then +! get the integer from the first line + if ((pline%id(1) .eq. "i") .and. (pline%ntokens .eq. 1)) then ! check if it is the only integer present in a line natoms = fdf_bintegers(pline, 1) write(*,*) "Number of atoms = ", natoms endif @@ -325,123 +338,40 @@ PROGRAM iochamp write(6,*) '------------------------------------------------------' + ! Determinants as a block. read directly from the input file +! under construction + if (fdf_block('determinants', bfdf)) then + ia = 1 + do while(fdf_bline(bfdf, pline)) + symbol(ia) = fdf_bnames(pline, 1) + do i= 1, 3 + xa(i,ia) = fdf_bvalues(pline, i) + enddo + ia = ia + 1 + enddo + na = ia - 1 - ! if (fdf_block('inline_xyz2', bfdf)) then - ! ! Forward reading - ! write(6,*) 'Reading an inline_xyz2 block ' - ! ia = 1 - - ! do while(fdf_bline(bfdf, pline)) - - ! if (pline%ntokens == 1) then - ! number_of_atoms = fdf_bintegers(pline, 1) - ! write(*,*) "Number of atoms", number_of_atoms - ! endif - ! na = number_of_atoms - - ! if (pline%ntokens == 4) then - ! symbol(ia) = fdf_bnames(pline, 1) - ! do i= 1, 3 - ! xa(i,ia) = fdf_bvalues(pline, i) - ! enddo - ! ia = ia + 1 - ! endif - ! enddo - - ! write(6,*) 'Inline XYZ2 Coordinates block:' - ! do ia= 1, na - ! write(6,'(A4,3F10.6)') symbol(ia), (xa(i,ia),i=1,3) - ! enddo - ! endif - - ! write(6,'(A)') + endif - ! write(6,*) '------------------------------------------------------' - - - ! if (fdf_block('molecule2', bfdf)) then - ! ! External file reading - ! write(6,*) 'beginning of external file coordinates block ' - ! ia = 1 - ! ! write(*,*) "linecount", fdf_block_linecount("molecule") - - ! do while((fdf_bline(bfdf, pline))) - - ! if (pline%ntokens == 1) then - ! number_of_atoms = fdf_bintegers(pline, 1) - ! write(*,*) "number of atoms", number_of_atoms - ! endif - ! na = number_of_atoms - - ! if (pline%ntokens == 4) then - ! symbol(ia) = fdf_bnames(pline, 1) - ! do i= 1, 3 - ! xa(i,ia) = fdf_bvalues(pline, i) - ! enddo - ! ia = ia + 1 - ! endif - ! enddo - ! endif - ! write(6,*) 'Coordinates from Molecule2 block: External file' - ! do ia= 1, na - ! write(6,'(A4,3F10.6)') symbol(ia), (xa(i,ia),i=1,3) - ! enddo - - - - ! write(6,'(A)') - - ! write(6,*) '------------------------------------------------------' - - + if (fdf_block('Coordinates', bfdf)) then + write(6,*) 'Coordinates:' + do ia = 1, na + write(6,'(A, 4x, 3F10.6)') symbol(ia), (xa(i,ia),i=1,3) + enddo + endif + write(6,*) '------------------------------------------------------' - ! if ( fdf_block('ListBlock',bfdf) ) then - ! i = 0 - ! do while ( fdf_bline(bfdf,pline) ) - ! i = i + 1 - ! na = fdf_bnlists(pline) - ! write(*,'(2(a,i0),a)') 'Listblock line: ',i,' has ',na,' lists' - ! do ia = 1 , na - ! j = -1 - - ! call fdf_bilists(pline,ia,j,isa) - ! write(*,'(tr5,2(a,i0),a)') 'list ',ia,' has ',j,' entries' - ! call fdf_bilists(pline,ia,j,isa) - ! write(*,'(tr5,a,1000(tr1,i0))') 'list: ',isa(1:j) - - ! end do - ! end do - ! end if - - - ! if ( fdf_islreal('list_floats') .and. fdf_islist('list_floats') & - ! .and. (.not. fdf_islinteger('list_floats')) ) then - ! na = -1 - ! call fdf_list('list_floats',na,listr) - ! write(*,'(tr1,a,i0,a)') 'list_floats has ',na,' entries' - ! if ( na < 2 ) stop 1 - ! call fdf_list('list_floats',na,listr) - ! write(*,'(tr5,a,1000(tr1,f12.8))') 'list_floats: ',listr(1:na) - ! else - ! write(*,*)'list_floats was not recognized' - ! stop 1 - ! end if - - ! write(6,'(A)') - - ! write(6,*) '------------------------------------------------------' - - write(6,'(A)') " Determinants Block" - - write(6,*) '------------------------------------------------------' - - if (.not. fdf_block('determinants', bfdf)) then - ! External file reading + if ( fdf_load_defined('determinants') ) then + ! External file reading + write(6,'(A)') " Determinants Block" + + write(6,*) '------------------------------------------------------' + write(6,*) 'Reading the determinants block from an external file ' ia = 1 ! call io_status() @@ -450,6 +380,7 @@ PROGRAM iochamp ! print*, "pline obtained", (fdf_bline(bfdf, pline)) open (unit=11,file=file_determinants, iostat=iostat, action='read' ) + if (iostat .ne. 0) error stop "Problem in opening the determinant file" read(11,*) temp1, temp2, nelectrons, temp3, nalpha read(11,*) temp1, ndeterminants, nexcitation @@ -479,12 +410,11 @@ PROGRAM iochamp read(11,*) temp1 if (temp1 == "end" ) write(*,*) "Determinant File read successfully " - - - close(11) - endif + endif ! condition if load determinant is present + + endif ! condition determinant block not present write(6,'(A)') diff --git a/src/keywords_m.F90 b/src/keywords_m.F90 index 7f36029..2a32a86 100644 --- a/src/keywords_m.F90 +++ b/src/keywords_m.F90 @@ -136,6 +136,12 @@ MODULE keywords public :: opt_method public :: multiple_adiag +! public :: nvalence + public :: excess_charge + public :: multiplicity + + + ! Following not yet added ! rlobx(y) Lobachevsky parameters for Fock expansion ! ipq,iacc_rej,icross,icuspg,idiv_v @@ -326,4 +332,8 @@ MODULE keywords character(len=20) :: opt_method +! integer :: nvalence + integer :: multiplicity + integer :: excess_charge + end module \ No newline at end of file diff --git a/src/periodic_table_m.f90 b/src/periodic_table_m.F90 similarity index 76% rename from src/periodic_table_m.f90 rename to src/periodic_table_m.F90 index bbfb299..5b27b30 100644 --- a/src/periodic_table_m.f90 +++ b/src/periodic_table_m.F90 @@ -18,6 +18,8 @@ module periodic_table character(len=20) :: name character(len=3) :: symbol integer(kind=i1) :: znuclear + integer(kind=i1) :: core + integer(kind=i1) :: nvalence integer(kind=i1) :: isotope = 1 ! currently supports only the most abundant ! future plans ! covalent_radii @@ -36,75 +38,75 @@ module periodic_table select case(sym) case("hydrogen", "Hydrogen", "H", "1") - atom = atom_t(name="hydrogen", symbol="H", atomic_mass=1.007825, znuclear=1) + atom = atom_t(name="hydrogen", symbol="H", atomic_mass=1.007825, znuclear=1, core=0, nvalence=1) case("helium", "Helium", "He", "2") - atom = atom_t(name="helium", symbol="He", atomic_mass=4.00260, znuclear=2) + atom = atom_t(name="helium", symbol="He", atomic_mass=4.00260, znuclear=2, core=0, nvalence=2) case("lithium", "Lithium", "Li", "3") - atom = atom_t(name="lithium", symbol="Li", atomic_mass=7.016003, znuclear=3) + atom = atom_t(name="lithium", symbol="Li", atomic_mass=7.016003, znuclear=3, core=2, nvalence=1) case("beryllium", "Beryllium", "Be", "4") - atom = atom_t(name="beryllium", symbol="Be", atomic_mass=9.012182, znuclear=4) + atom = atom_t(name="beryllium", symbol="Be", atomic_mass=9.012182, znuclear=4, core=2, nvalence=2) case("boron", "Boron", "B", "5") - atom = atom_t(name="boron", symbol="B", atomic_mass=11.009305, znuclear=5) + atom = atom_t(name="boron", symbol="B", atomic_mass=11.009305, znuclear=5, core=2, nvalence=3) case("carbon", "Carbon", "C", "6") - atom = atom_t(name="carbon", symbol="C", atomic_mass=12.000000, znuclear=6) + atom = atom_t(name="carbon", symbol="C", atomic_mass=12.000000, znuclear=6, core=2, nvalence=4) case("nitrogen", "Nitrogen", "N", "7") - atom = atom_t(name="nitrogen", symbol="N", atomic_mass=14.003074, znuclear=7) + atom = atom_t(name="nitrogen", symbol="N", atomic_mass=14.003074, znuclear=7, core=2, nvalence=5) case("oxygen", "Oxygen", "O", "8") - atom = atom_t(name="oxygen", symbol="O", atomic_mass=15.994915, znuclear=8) + atom = atom_t(name="oxygen", symbol="O", atomic_mass=15.994915, znuclear=8, core=2, nvalence=6) case("fluorine", "Fluorine", "F", "9") - atom = atom_t(name="fluorine", symbol="F", atomic_mass=18.9984032, znuclear=9) + atom = atom_t(name="fluorine", symbol="F", atomic_mass=18.9984032, znuclear=9, core=2, nvalence=7) case("neon", "Neon", "Ne", "10") - atom = atom_t(name="neon", symbol="Ne", atomic_mass=19.992435, znuclear=10) + atom = atom_t(name="neon", symbol="Ne", atomic_mass=19.992435, znuclear=10, core=2, nvalence=8) case("sodium", "Sodium", "Na", "11") - atom = atom_t(name="sodium", symbol="Na", atomic_mass=22.989767, znuclear=11) + atom = atom_t(name="sodium", symbol="Na", atomic_mass=22.989767, znuclear=11, core=10, nvalence=1) case("magnesium", "Magnesium", "Mg", "12") - atom = atom_t(name="magnesium", symbol="Mg", atomic_mass=23.985042, znuclear=12) + atom = atom_t(name="magnesium", symbol="Mg", atomic_mass=23.985042, znuclear=12, core=10, nvalence=2) case("aluminum", "Aluminum", "aluminium", "Aluminium", "Al", "13") - atom = atom_t(name="aluminum", symbol="Al", atomic_mass=26.981540, znuclear=13) + atom = atom_t(name="aluminum", symbol="Al", atomic_mass=26.981540, znuclear=13, core=10, nvalence=3) case("silicon", "Silicon", "Si", "14") - atom = atom_t(name="silicon", symbol="Si", atomic_mass=27.976927, znuclear=14) + atom = atom_t(name="silicon", symbol="Si", atomic_mass=27.976927, znuclear=14, core=10, nvalence=4) case("phosphorus", "Phosphorus", "P", "15") - atom = atom_t(name="phosphorus",symbol="P", atomic_mass=30.973762, znuclear=15) + atom = atom_t(name="phosphorus",symbol="P", atomic_mass=30.973762, znuclear=15, core=10, nvalence=5) case("sulfur", "Sulfur", "sulphur", "Sulphur", "S", "16") - atom = atom_t(name="sulfur", symbol="S", atomic_mass=31.972070, znuclear=16) + atom = atom_t(name="sulfur", symbol="S", atomic_mass=31.972070, znuclear=16, core=10, nvalence=6) case("chlorine", "Chlorine", "Cl", "17") - atom = atom_t(name="chlorine", symbol="Cl", atomic_mass=34.968852, znuclear=17) + atom = atom_t(name="chlorine", symbol="Cl", atomic_mass=34.968852, znuclear=17, core=10, nvalence=7) case("argon", "Argon", "Ar", "18") - atom = atom_t(name="argon", symbol="Ar", atomic_mass=39.962384, znuclear=18) + atom = atom_t(name="argon", symbol="Ar", atomic_mass=39.962384, znuclear=18, core=10, nvalence=8) case default error stop "Unknown element's atomic number or symbol" diff --git a/src/test-champ.inp b/src/test-champ.inp index 144b21f..c981c10 100644 --- a/src/test-champ.inp +++ b/src/test-champ.inp @@ -1,16 +1,29 @@ + title "A Sample QMC input file parsed by libfdf interfaced to CHAMP" -pool ./pool -pseudopot BFD -basis BFD-T-normf0 +pool ./pool +pseudopot BFD +basis BFD-T-normf0 -%include global.inp - -%block molecule < caffeine.xyz +%include global.inp load molecule caffeine.xyz load basis BFD-T-normf0.bas -load determinants TZ_1M_500.det +#load determinants TZ_1M_500.det + + +%block molecule < caffeine.xyz + + + + +#%block pseudopotential + + +#%end block + + + %module optimization optimize_wavefunction 1 diff --git a/src/test-champ.out b/src/test-champ.out index b759668..9cae397 100644 --- a/src/test-champ.out +++ b/src/test-champ.out @@ -3,7 +3,7 @@ pool ./pool molecule caffeine.xyz pseudopot default.psp # default value basis BFD-T-normf0.bas -determinants TZ_1M_500.det +determinants default.det # default value orbitals default.orb # default value jastrow default.jas # default value jastrow_deriv default.jasder # default value @@ -30,6 +30,8 @@ ifock 1 # default value tau 0.3999999911E-01 etrial -408.1744362 eV # above item originally: etrial -15.00000000 Ha +excess_charge 0 # default value +multiplicity 1 # default value %block molecule %block molecule 24 @@ -59,3 +61,5 @@ etrial -408.1744362 eV H 3.2249 1.5791 0.8255 H 2.6431 2.5130 -0.5793 %endblock molecule +#:block? determinants F +#:defined? determinants F