From 9c10e1aa99b4bd499cfbb5770678e7323a9fa220 Mon Sep 17 00:00:00 2001 From: Felix Moncada Date: Wed, 27 Nov 2024 17:29:00 +0100 Subject: [PATCH 1/3] Adding options to add custom potentials and particles. Also setting the mass from the isotope as the default value for the nuclear mass. Solves issues #31 and #100 Added the option to define "eta" (particles per orbital) directly from the input. As example, one test defining positron in closed shell. Removed the TIP4P particles and potentials from lib/databases, and now they are passed directly from the input (see these tests as examples). Added a test for the isotopes mass. Merged externalPotential and interPotential into a single module GTFpotential Now the script moves the fchk to/from the scratch. Also, when looking for the coefficients from these files, now it does not stop when it fails to find a fchk for a species, instead does a HCore guess. However, the fchks tests are still missing (related to #103) --- bin/lowdin | 40 ++- lib/basis/H2O-1S1P1D | 44 --- lib/basis/PSX-DZ | 52 ++- lib/dataBases/constantsOfCoupling.lib | 16 - lib/dataBases/elementalParticles.lib | 36 +- src/core/ConstantsOfCoupling.f90 | 86 +++-- src/core/ElementalParticle.f90 | 97 +++--- src/core/ExternalPotential.f90 | 378 -------------------- src/core/GTFPotential.f90 | 309 +++++++++++++++++ src/core/InputManager.f90 | 23 +- src/core/InterPotential.f90 | 298 ---------------- src/core/MolecularSystem.f90 | 62 ++-- src/core/Particle.f90 | 53 ++- src/core/Species.f90 | 13 +- src/core/Units.f90 | 1 + src/ints/DirectIntegralManager.f90 | 6 +- src/ints/G12Integrals.f90 | 12 +- src/ints/IntegralManager.f90 | 6 +- src/ints/Libint2Interface.f90 | 38 +- src/scf/DensityMatrixSCFGuess.f90 | 15 +- src/scf/OrbitalLocalizer.f90 | 17 +- src/scf/SingleSCF.f90 | 2 +- test/CHDTHemu.massTest.lowdin | 23 ++ test/CHDTHemu.massTest.py | 84 +++++ test/He2-C60potential-NOCI.lowdin | 4 +- test/He3-C60potential-NOCI.lowdin | 6 +- test/Ps2F2.RHF-customEta.lowdin | 17 + test/Ps2F2.RHF-customEta.py | 68 ++++ test/TIP4P-dimer-singlet-NOCI.lowdin | 459 ++++++++++++++++++++++++- test/TIP4P-dimer-singlet-direct.lowdin | 447 +++++++++++++++++++++++- test/TIP4P-dimer-singlet-direct.py | 30 +- test/TIP4P-dimer-singlet-memory.lowdin | 447 +++++++++++++++++++++++- test/TIP4P-dimer-singlet-memory.py | 30 +- test/TIP4P-dimer-singlet.lowdin | 447 +++++++++++++++++++++++- test/TIP4P-dimer-singlet.py | 30 +- test/TIP4P-dimer-triplet-NOCI.lowdin | 448 +++++++++++++++++++++++- test/TIP4P-dimer-triplet-direct.lowdin | 193 ++++++++++- test/TIP4P-dimer-triplet-direct.py | 24 +- test/TIP4P-dimer-triplet-memory.lowdin | 201 ++++++++++- test/TIP4P-dimer-triplet-memory.py | 24 +- test/TIP4P-dimer-triplet.lowdin | 193 ++++++++++- test/TIP4P-dimer-triplet.py | 24 +- 42 files changed, 3715 insertions(+), 1088 deletions(-) delete mode 100644 lib/basis/H2O-1S1P1D delete mode 100644 src/core/ExternalPotential.f90 create mode 100644 src/core/GTFPotential.f90 delete mode 100644 src/core/InterPotential.f90 create mode 100644 test/CHDTHemu.massTest.lowdin create mode 100644 test/CHDTHemu.massTest.py create mode 100644 test/Ps2F2.RHF-customEta.lowdin create mode 100644 test/Ps2F2.RHF-customEta.py diff --git a/bin/lowdin b/bin/lowdin index 3a709f9c..6377f7d1 100755 --- a/bin/lowdin +++ b/bin/lowdin @@ -204,6 +204,9 @@ if [ $extFile="lowdin" ]; then else if(toupper($i)==toupper("m")){ printf("\tInputParticle_mass = %s\n",toupper($(i+2)) ) } + else if(toupper($i)==toupper("eta")){ + printf("\tInputParticle_eta = %s\n",toupper($(i+2)) ) + } else if(toupper($i)==toupper("omega")){ printf("\tInputParticle_omega = %s\n",toupper($(i+2)) ) } @@ -381,10 +384,10 @@ if [ $extFile="lowdin" ]; then ' | $SED "s/,/./g" >> $nameFile.aux ########################################### - # Check custom basis sets in the input + # Check custom basis sets/potentials in the input ########################################### - BASIS_NAMES=(`gawk '($1~/BASIS/){print toupper($2)}' $nameFile`) + BASIS_NAMES=(`gawk '($1~/^BASIS$/){print toupper($2)}' $nameFile`) if [ ${#BASIS_NAMES[@]} -gt "0" ] then for BASIS_NAME in ${BASIS_NAMES[@]} @@ -401,6 +404,23 @@ if [ $extFile="lowdin" ]; then done fi + POTENTIALS_NAMES=(`gawk '($1~/^POTENTIAL$/){print toupper($2)}' $nameFile`) + if [ ${#POTENTIALS_NAMES[@]} -gt "0" ] + then + for POTENTIALS_NAME in ${POTENTIALS_NAMES[@]} + do + if [ -e $LOWDIN_DATA/basis/$POTENTIALS_NAME ] + then + echo "## ERROR: ## The custom potential file already exists in " $LOWDIN_DATA/potentials/$POTENTIALS_NAME + echo "Modify the POTENTIALS block in your input and select a different name" + exit 1 + fi + gawk '($1~/POTENTIALS/ && toupper($2)~/^'$POTENTIALS_NAME'$/){flag=1; next} + ($0~/END/){flag=0}; + (flag==1){print toupper($0)}' $nameFile > $POTENTIALS_NAME.$PID + done + fi + ########################################### # Exec lowdin.x @@ -421,6 +441,7 @@ if [ $extFile="lowdin" ]; then cp $nameFile*.vec $LOWDIN_SCRATCH/$nameFile &> /dev/null cp $nameFile*.plainvec $LOWDIN_SCRATCH/$nameFile &> /dev/null + cp $nameFile*.fchk $LOWDIN_SCRATCH/$nameFile &> /dev/null cp $nameFile*.val $LOWDIN_SCRATCH/$nameFile &> /dev/null cp $nameFile*.dens $LOWDIN_SCRATCH/$nameFile &> /dev/null cp $nameFile*.sup $LOWDIN_SCRATCH/$nameFile &> /dev/null @@ -432,7 +453,9 @@ if [ $extFile="lowdin" ]; then mv $nameFile*.over $LOWDIN_SCRATCH/$nameFile &> /dev/null mv $nameFile*.kin $LOWDIN_SCRATCH/$nameFile &> /dev/null mv $nameFile*.coeff $LOWDIN_SCRATCH/$nameFile &> /dev/null - #PID to avoid basis duplicates in simultaneous calculations + cp $nameFile*.gms.bs $LOWDIN_SCRATCH/$nameFile &> /dev/null + + #PID to avoid basis/potentials duplicates in simultaneous calculations if [ ${#BASIS_NAMES[@]} -gt "0" ] then for BASIS_NAME in ${BASIS_NAMES[@]} @@ -440,12 +463,14 @@ if [ $extFile="lowdin" ]; then mv $BASIS_NAME.$PID $LOWDIN_SCRATCH/$nameFile/$BASIS_NAME &> /dev/null done fi - - if [ -e $nameFile.gms.bs ] + if [ ${#POTENTIALS_NAMES[@]} -gt "0" ] then - cp $nameFile.gms.bs $LOWDIN_SCRATCH/$nameFile + for POTENTIALS_NAME in ${POTENTIALS_NAMES[@]} + do + mv $POTENTIALS_NAME.$PID $LOWDIN_SCRATCH/$nameFile/$POTENTIALS_NAME &> /dev/null + done fi - + # setting default number of cores for OpenMP if [ -z "$OMP_NUM_THREADS" ]; then @@ -515,6 +540,7 @@ if [ $extFile="lowdin" ]; then mv $LOWDIN_SCRATCH/$nameFile/$nameFile*.47 $currentPath &> 2 mv $LOWDIN_SCRATCH/$nameFile/*.vec $currentPath &> 2 mv $LOWDIN_SCRATCH/$nameFile/*.plainvec $currentPath &> 2 + mv $LOWDIN_SCRATCH/$nameFile/*.fchk $currentPath &> 2 # mv $LOWDIN_SCRATCH/$nameFile/*.val $currentPath &> 2 mv $LOWDIN_SCRATCH/$nameFile/*.NOCI.coords $currentPath &> 2 mv $LOWDIN_SCRATCH/$nameFile/*.NOCI.s* $currentPath &> 2 diff --git a/lib/basis/H2O-1S1P1D b/lib/basis/H2O-1S1P1D deleted file mode 100644 index cc7027ba..00000000 --- a/lib/basis/H2O-1S1P1D +++ /dev/null @@ -1,44 +0,0 @@ -O-HYDROGEN H_1 (1S) BASIS TYPE: 2 -3 -1 0 1 -14.509888498676842 1.0 -2 1 1 -6.885507269761004 1.0 -3 2 1 -9.023681376783887 1.0 - -O-HYDROGEN H-A_1 (1S) BASIS TYPE: 2 -3 -1 0 1 -14.509888498676842 1.0 -2 1 1 -6.885507269761004 1.0 -3 2 1 -9.023681376783887 1.0 - -O-HYDROGEN H-B_1 (1S) BASIS TYPE: 2 -3 -1 0 1 -14.509888498676842 1.0 -2 1 1 -6.885507269761004 1.0 -3 2 1 -9.023681376783887 1.0 - -O-H-TIP X0.5+ (1S) BASIS TYPE: 2 -3 -1 0 1 -14.509888498676842 1.0 -2 1 1 -6.885507269761004 1.0 -3 2 1 -9.023681376783887 1.0 - -O-H-TIP Y0.5+ (1S) BASIS TYPE: 2 -3 -1 0 1 -14.509888498676842 1.0 -2 1 1 -6.885507269761004 1.0 -3 2 1 -9.023681376783887 1.0 diff --git a/lib/basis/PSX-DZ b/lib/basis/PSX-DZ index 35b403de..0d5600b6 100644 --- a/lib/basis/PSX-DZ +++ b/lib/basis/PSX-DZ @@ -1,5 +1,53 @@ -O-POSITRON E+ (5S) BASIS TYPE: 3 -# (5S)-[5S] +O-POSITRON E+ (5S3P2D) BASIS TYPE: 3 +# (5S3P2D)-[5S3P2D] +10 +1 0 1 +.0189693659 1.0 +1 0 1 +.05186863351164733038 1.0 +1 0 1 +.14182630861506996771 1.0 +1 0 1 +.38780088183468771642 1.0 +1 0 1 +1.06037818667291718709 1.0 +1 1 1 +.0590955656 1.0 +1 1 1 +.16127967470846848928 1.0 +1 1 1 +.44015372744092004227 1.0 +1 2 1 +.11654481880000000000 1.0 +1 2 1 +.31287113568838127412 1.0 + +O-POSITRON E+A (5S3P2D) BASIS TYPE: 3 +# (5S3P2D)-[5S3P2D] +10 +1 0 1 +.0189693659 1.0 +1 0 1 +.05186863351164733038 1.0 +1 0 1 +.14182630861506996771 1.0 +1 0 1 +.38780088183468771642 1.0 +1 0 1 +1.06037818667291718709 1.0 +1 1 1 +.0590955656 1.0 +1 1 1 +.16127967470846848928 1.0 +1 1 1 +.44015372744092004227 1.0 +1 2 1 +.11654481880000000000 1.0 +1 2 1 +.31287113568838127412 1.0 + +O-POSITRON E+B (5S3P2D) BASIS TYPE: 3 +# (5S3P2D)-[5S3P2D] 10 1 0 1 .0189693659 1.0 diff --git a/lib/dataBases/constantsOfCoupling.lib b/lib/dataBases/constantsOfCoupling.lib index bda5b291..1f940075 100644 --- a/lib/dataBases/constantsOfCoupling.lib +++ b/lib/dataBases/constantsOfCoupling.lib @@ -4525,19 +4525,3 @@ LAMBDA = 2.0 PARTICLESFRACTION = 0.5 / -&SPECIE - NAME = "HA-TIP" - SYMBOL = "X0.5+" - KAPPA = -1.0 - ETA = 1.0 - LAMBDA = 1.0 - PARTICLESFRACTION = 1 -/ -&SPECIE - NAME = "HB-TIP" - SYMBOL = "Y0.5+" - KAPPA = -1.0 - ETA = 1.0 - LAMBDA = 1.0 - PARTICLESFRACTION = 1 -/ diff --git a/lib/dataBases/elementalParticles.lib b/lib/dataBases/elementalParticles.lib index d6630a9b..54376343 100644 --- a/lib/dataBases/elementalParticles.lib +++ b/lib/dataBases/elementalParticles.lib @@ -171,7 +171,7 @@ SYMBOL = "HEA3" CATEGORY = "LEPTON" CHARGE = 1 - MASS = 5494.892576965 + MASS = 5495.8851 SPIN = 0.5 / &PARTICLE @@ -179,7 +179,7 @@ SYMBOL = "HEB3" CATEGORY = "LEPTON" CHARGE = 1 - MASS = 5494.892576965 + MASS = 5495.8851 SPIN = 0.5 / &PARTICLE @@ -187,7 +187,7 @@ SYMBOL = "HES3" CATEGORY = "LEPTON" CHARGE = 1 - MASS = 5494.892576965 + MASS = 5495.8851 SPIN = 0.5 / &PARTICLE @@ -195,7 +195,7 @@ SYMBOL = "HEA4" CATEGORY = "LEPTON" CHARGE = 1 - MASS = 7292.327967297 + MASS = 7294.2994 SPIN = 0.5 / &PARTICLE @@ -203,7 +203,7 @@ SYMBOL = "HEB4" CATEGORY = "LEPTON" CHARGE = 1 - MASS = 7292.327967297 + MASS = 7294.2994 SPIN = 0.5 / &PARTICLE @@ -211,31 +211,7 @@ SYMBOL = "HES4" CATEGORY = "LEPTON" CHARGE = 1 - MASS = 7292.327967297 - SPIN = 0.5 -/ -&PARTICLE - NAME = "HA-TIP" - SYMBOL = "X0.5+" - CATEGORY = "LEPTON" - CHARGE = 0.5564 - MASS = 1836.15267247 - SPIN = 0.5 -/ -&PARTICLE - NAME = "HB-TIP" - SYMBOL = "Y0.5+" - CATEGORY = "LEPTON" - CHARGE = 0.5564 - MASS = 1836.15267247 - SPIN = -0.5 -/ -&PARTICLE - NAME = "M-TIP" - SYMBOL = "X1.1-" - CATEGORY = "LEPTON" - CHARGE = -1.1128 - MASS = 1836.15267247 + MASS = 7294.2994 SPIN = 0.5 / &PARTICLE diff --git a/src/core/ConstantsOfCoupling.f90 b/src/core/ConstantsOfCoupling.f90 index 3e3bba2c..0c190802 100644 --- a/src/core/ConstantsOfCoupling.f90 +++ b/src/core/ConstantsOfCoupling.f90 @@ -71,59 +71,53 @@ subroutine ConstantsOfCoupling_load( this, symbolSelected ) !! Looking for library inquire(file=trim(CONTROL_instance%DATA_DIRECTORY)//"/dataBases/constantsOfCoupling.lib", exist=existFile) - if ( existFile ) then - - !! Open library - open(unit=10, file=trim(CONTROL_instance%DATA_DIRECTORY)//"/dataBases/constantsOfCoupling.lib", status="old", form="formatted" ) - - !! Read information - symbol = "NONE" - stat = 0 - - do while(trim(symbol) /= trim(symbolSelected)) - - !! Setting defaults - name = "NONE" - kappa = 0 - eta = 0 - particlesFraction = 1 - - if (stat == -1 ) then - - call ConstantsOfCoupling_exception( ERROR, "Elemental particle: "//trim(symbolSelected)//" NOT found!!", "In ConstantsOfCoupling at load function.") - this%isInstanced = .false. + if ( .not. existFile ) call ConstantsOfCoupling_exception( ERROR, "LOWDIN library not found!! please export lowdinvars.sh file.", "In ConstantsOfCoupling at load function.") + + !! Open library + open(unit=10, file=trim(CONTROL_instance%DATA_DIRECTORY)//"/dataBases/constantsOfCoupling.lib", status="old", form="formatted" ) - end if - - read(10,NML=specie, iostat=stat) - - if (stat > 0 ) then - - call ConstantsOfCoupling_exception( ERROR, "Failed reading ConstantsOfCouplings.lib file!! please check this file.", "In ConstantsOfCoupling at load function.") - - end if + !! Read information + symbol = "NONE" + stat = 0 - end do + do while(trim(symbol) /= trim(symbolSelected)) - !! Set object variables - this%name = name - this%symbol = symbol - this%kappa = kappa - this%eta = eta - this%lambda = lambda - this%particlesFraction = particlesFraction - - !! Debug information. - !! call ConstantsOfCoupling_show(this) - - close(10) + !! Setting defaults + name = "NONE" + kappa = 0 + eta = 0 + particlesFraction = 1 - else + if (stat == -1 ) then - call ConstantsOfCoupling_exception( ERROR, "LOWDIN library not found!! please export lowdinvars.sh file.", "In ConstantsOfCoupling at load function.") + ! call ConstantsOfCoupling_exception( WARNING, "Elemental particle: "//trim(symbolSelected)//" NOT found!!", "Setting default values") + this%isInstanced=.false. + exit + end if - end if + read(10,NML=specie, iostat=stat) + if (stat > 0 ) then + + call ConstantsOfCoupling_exception( ERROR, "Failed reading ConstantsOfCouplings.lib file!! please check this file.", "In ConstantsOfCoupling at load function.") + + end if + + end do + + !! Set object variables + this%name = name + this%symbol = symbol + this%kappa = kappa + this%eta = eta + this%lambda = lambda + this%particlesFraction = particlesFraction + + !! Debug information. + ! call ConstantsOfCoupling_show(this) + + close(10) + !! Done end subroutine ConstantsOfCoupling_load diff --git a/src/core/ElementalParticle.f90 b/src/core/ElementalParticle.f90 index ce1b69b9..5b961321 100644 --- a/src/core/ElementalParticle.f90 +++ b/src/core/ElementalParticle.f90 @@ -40,6 +40,7 @@ module ElementalParticle_ real(8) :: mass real(8) :: charge real(8) :: spin + logical :: custom end type ElementalParticle public :: & @@ -59,6 +60,7 @@ subroutine ElementalParticle_load( this, symbolSelected ) character(*) :: symbolSelected logical :: existFile + logical :: custom integer :: stat integer :: i @@ -81,59 +83,64 @@ subroutine ElementalParticle_load( this, symbolSelected ) !! Looking for library inquire(file=trim(CONTROL_instance%DATA_DIRECTORY)//trim(CONTROL_instance%ELEMENTAL_PARTICLES_DATABASE), exist=existFile) - if ( existFile ) then + if ( .not. existFile ) call ElementalParticle_exception( ERROR, "LOWDIN library not found!! please export lowdinvars.sh file.", "In ElementalParticle at load function.") - !! Open library - open(unit=10, file=trim(CONTROL_instance%DATA_DIRECTORY)//trim(CONTROL_instance%ELEMENTAL_PARTICLES_DATABASE), status="old", form="formatted" ) - - !! Read information - symbol = "NONE" - stat = 0 - - do while(trim(symbol) /= trim(symbolSelected)) + !! Open library + open(unit=10, file=trim(CONTROL_instance%DATA_DIRECTORY)//trim(CONTROL_instance%ELEMENTAL_PARTICLES_DATABASE), status="old", form="formatted" ) + + !! Read information + symbol = "NONE" + stat = 0 + + do while(trim(symbol) /= trim(symbolSelected)) + + !! Setting defaults + name = "NONE" + category = "NONE" + mass = -1 + charge = 0 + spin = 0 + custom = .false. - !! Setting defaults - name = "NONE" - category = "NONE" - mass = -1 - charge = 0 - spin = 0 - - if (stat == -1 ) then - - call ElementalParticle_exception( ERROR, "Elemental particle: "//trim(symbolSelected)//" NOT found!!", "In ElementalParticle at load function.") + if (stat == -1 ) then - end if + ! call ElementalParticle_exception( WARNING, "Elemental particle: "//trim(symbolSelected)//" NOT found in ElementalParticles.lib", "Setting default values") + name = trim(symbolSelected) + symbol = trim(symbolSelected) + category = "FERMION" + mass = 1.0 + charge = 1.0 + spin = 0.5 + custom = .true. + + exit - read(10,NML=particle, iostat=stat) - - if (stat > 0 ) then - - call ElementalParticle_exception( ERROR, "Failed reading ElementalParticles.lib file!! please check this file.", "In ElementalParticle at load function.") - - end if + end if - end do - - !! Set object variables - this%name = name - this%symbol = symbol - this%category = category - this%mass = mass - this%charge = charge - this%spin = spin - - !! Debug information. - !! call ElementalParticle_show(this) - - close(10) - - else + read(10,NML=particle, iostat=stat) + + if (stat > 0 ) then + + call ElementalParticle_exception( ERROR, "Failed reading ElementalParticles.lib file!! please check this file.", "In ElementalParticle at load function.") - call ElementalParticle_exception( ERROR, "LOWDIN library not found!! please export lowdinvars.sh file.", "In ElementalParticle at load function.") + end if - end if + end do + !! Set object variables + this%name = name + this%symbol = symbol + this%category = category + this%mass = mass + this%charge = charge + this%spin = spin + this%custom = custom + + !! Debug information. + ! call ElementalParticle_show(this) + + close(10) + !! Done end subroutine ElementalParticle_load diff --git a/src/core/ExternalPotential.f90 b/src/core/ExternalPotential.f90 deleted file mode 100644 index f6459bfa..00000000 --- a/src/core/ExternalPotential.f90 +++ /dev/null @@ -1,378 +0,0 @@ -!!****************************************************************************** -!! This code is part of LOWDIN Quantum chemistry package -!! -!! this program has been developed under direction of: -!! -!! Prof. A REYES' Lab. Universidad Nacional de Colombia -!! http://www.qcc.unal.edu.co -!! Prof. R. FLORES' Lab. Universidad de Guadajara -!! http://www.cucei.udg.mx/~robertof -!! -!! Todos los derechos reservados, 2013 -!! -!!****************************************************************************** - -!> @brief This module contains all the routines to handle external potentials -!! @author E. F. Posada (efposadac@unal.edu.co) -!! @version 2.0 -!! Creation data : 06-08-10 -!! -!! History change: -!! -!! - 06-08-10 : Sergio A. Gonzalez ( sergmonic@gmail.com ) -!! -# Creacioon del modulo y metodos basicos -!! - 2011-02-14 : Fernando Posada ( efposadac@unal.edu.co ) -!! -# Reescribe y adapta el módulo para su inclusion en Lowdin -module ExternalPotential_ - use ContractedGaussian_ - use String_ - use Matrix_ - use Units_ - use Exception_ - implicit none - - type, public :: ExternalPot - character(20) :: name - character(50) :: specie - character(50) :: ttype - character(50) :: units - integer :: numOfComponents - integer :: iter - type(ContractedGaussian), allocatable :: gaussianComponents(:) - end type - - type, public :: ExternalPotential - integer :: ssize - type(ExternalPot), allocatable :: potentials(:) - logical :: isInstanced - end type - - type(ExternalPotential), public, target :: ExternalPotential_instance - -contains - - - !> - !! @brief Constructor by default - !! @param this - subroutine ExternalPotential_constructor(numberOfPotentials) - implicit none - - integer :: numberOfPotentials - - ExternalPotential_instance%ssize = numberOfPotentials - allocate(ExternalPotential_instance%potentials(numberOfPotentials)) - ExternalPotential_instance%isInstanced = .true. - - end subroutine ExternalPotential_constructor - - !> - !! @brief Destroys the object - !! @param this - subroutine ExternalPotential_destructor() - implicit none - - integer :: i - - do i = 1, ExternalPotential_instance%ssize - if (allocated(ExternalPotential_instance%potentials(i)%gaussianComponents)) deallocate(ExternalPotential_instance%potentials(i)%gaussianComponents) - end do - - if (allocated(ExternalPotential_instance%potentials) ) deallocate(ExternalPotential_instance%potentials) - ExternalPotential_instance%isInstanced=.false. - - end subroutine ExternalPotential_destructor - - !> - !! @brief Shows information of the object - !! @param this - subroutine ExternalPotential_show() - implicit none - type(ExternalPot), pointer :: this - integer :: potId, i - - do potId = 1, ExternalPotential_instance%ssize - this => ExternalPotential_instance%potentials(potId) - - print *,"" - print *,"=======" - print *, "External Potential for ", trim(this%specie), " : ", trim(this%name) - print *, "Type : ", trim(this%ttype) - write(6,"(T10,A20,A10,A10,A10,A10,A20)") "Exponent", "l", "R_x", "R_y", "R_z", "Factor" - - do i=1,this%numOfComponents - write(6,"(T10,F20.10,I10,F10.5,F10.5,F10.5,F20.10)") & - this%gaussianComponents(i)%orbitalExponents(1), & - this%gaussianComponents(i)%angularMoment, this%gaussianComponents(i)%origin(:), & - this%gaussianComponents(i)%contractionCoefficients(1) - end do - end do - - end subroutine ExternalPotential_show - - !> - !! @brief loads information from the input file - !! @param this - !! @author E. F. Posada, 2013 - subroutine ExternalPotential_load(potId, name, species) - implicit none - integer :: potId - character(*) :: name - character(*) :: species - - type(ExternalPot), pointer :: this - integer :: status, i, j - character(150) :: fileName - character(20) :: token - character(10) :: symbol - logical :: existFile, found - - this => ExternalPotential_instance%potentials(potId) - - this%name= trim(name) - this%specie= trim(species) - this%ttype="" - this%units="bohr" - this%numOfComponents=0 - this%iter=1 - - fileName = trim( trim( CONTROL_instance%DATA_DIRECTORY ) // & - trim(CONTROL_instance%POTENTIALS_DATABASE)// String_getUppercase(trim(name))) - - inquire(file=trim(fileName), exist = existFile) - if(existFile) then - - !! Open File - open(unit=30, file=trim(fileName), status="old",form="formatted") - rewind(30) - - found = .false. - - !! Open element and Find the proper potential - do while(found .eqv. .false.) - read(30,*, iostat=status) token - symbol = token(3:) - - !! Some debug information in case of error! - if (status > 0 ) then - - call ExternalPotential_exception(ERROR, & - "ERROR reading ExternalPotential file: "//trim(this%name)//& - " Please check that file!","ExternalPotential module at Load function.") - - end if - - if (status == -1 ) then - - call ExternalPotential_exception(ERROR, & - "The ExternalPotential: "//trim(this%name)//& - " for: "//trim(species)//& - " was not found!","ExternalPotential module at Load function.") - - end if - - if(trim(token(1:2)) == "O-") then - if(trim(symbol) == trim(species)) then - found = .true. - - end if - - end if - - end do - - !! Neglect any comment - token = "#" - do while(trim(token(1:1)) == "#") - - read(30,*) token - - end do - - !! Start reading Potential - backspace(30) - - read(30,*, iostat=status) this%numOfComponents - - !! Some debug information in case of error! - if (status > 0 ) then - - call ExternalPotential_exception(ERROR, & - "ERROR reading ExternalPotential file: "//trim(this%name)//& - " Please check that file!","ExternalPotential module at Load function.") - - end if - - allocate(this%gaussianComponents(this%numOfComponents)) - - do i = 1, this%numOfComponents - - read(30,*,iostat=status) this%gaussianComponents(i)%id, & - this%gaussianComponents(i)%angularMoment - this%gaussianComponents(i)%length = 1 - - !! Some debug information in case of error! - if (status > 0 ) then - - call ExternalPotential_exception(ERROR, & - "ERROR reading ExternalPotential file: "//trim(this%name)//& - " Please check that file!","ExternalPotential module at Load function.") - - end if - - allocate(this%gaussianComponents(i)%orbitalExponents(this%gaussianComponents(i)%length)) - allocate(this%gaussianComponents(i)%contractionCoefficients(this%gaussianComponents(i)%length)) - - do j = 1, this%gaussianComponents(i)%length - - read(30,*,iostat=status) this%gaussianComponents(i)%orbitalExponents(j), & - this%gaussianComponents(i)%contractionCoefficients(j) - read(30,*,iostat=status) this%gaussianComponents(i)%origin - - !! Some debug information in case of error! - if (status > 0 ) then - - call ExternalPotential_exception(ERROR, & - "ERROR reading ExternalPotential file: "//trim(this%name)//& - " Please check that file!","ExternalPotential module at Load function.") - - end if - - end do - - - !! Calculates the number of Cartesian orbitals, by dimensionality - select case(CONTROL_instance%DIMENSIONALITY) - case(3) - this%gaussianComponents(i)%numCartesianOrbital = ( ( this%gaussianComponents(i)%angularMoment + 1_8 )*( this%gaussianComponents(i)%angularMoment + 2_8 ) ) / 2_8 - case(2) - this%gaussianComponents(i)%numCartesianOrbital = ( ( this%gaussianComponents(i)%angularMoment + 1_8 ) ) - case(1) - this%gaussianComponents(i)%numCartesianOrbital = 1 - case default - call ExternalPotential_exception( ERROR, & - "Class object ExternalPotential in load function",& - "This Dimensionality is not available") - end select - - !! Normalize - allocate(this%gaussianComponents(i)%contNormalization(this%gaussianComponents(i)%numCartesianOrbital)) - allocate(this%gaussianComponents(i)%primNormalization(this%gaussianComponents(i)%length, & - this%gaussianComponents(i)%length*this%gaussianComponents(i)%numCartesianOrbital)) - - this%gaussianComponents(i)%contNormalization = 1.0_8 - this%gaussianComponents(i)%primNormalization = 1.0_8 - - call ContractedGaussian_normalizePrimitive(this%gaussianComponents(i)) - call ContractedGaussian_normalizeContraction(this%gaussianComponents(i)) - - !! DEBUG - ! call ContractedGaussian_showInCompactForm(ExternalPotential_instance%potentials(potId)%gaussianComponents(i)) - - end do - - close(30) - - !!DONE - - else - - call ExternalPotential_exception(ERROR, & - "The ExternalPotential file: "//trim(name)//& - " was not found!","ExternalPotential module at Load function.") - - end if - - end subroutine ExternalPotential_load - -! !> -! !! @brief -! !! @param this -! function ExternalPotential_getPotential( this, coords ) result(output) -! implicit none -! type(ExternalPotential) :: this -! real(8) :: coords(3) -! real(8) :: output - -! ! integer :: i - -! ! output=0.0 - -! ! do i=1, this%gaussianComponents%length -! ! output = output+( this%gaussianComponents%contractionCoefficients(i)* & -! ! exp(-this%gaussianComponents%primitives(i)%orbitalExponent*( dot_product(coords,coords) ) ) ) -! ! end do - -! end function ExternalPotential_getPotential - - -! ! function ExternalPotential_getInteractionMtx(this, contractions) result(output) -! subroutine ExternalPotential_getInteractionMtx( this, contractions ) -! implicit none -! type(ExternalPotential) :: this -! type(ContractedGaussian) :: contractions(:) -! type(Matrix) :: output - -! ! integer :: i, j, k, l, m, a, b -! ! integer :: numContractions -! ! real(8), allocatable :: auxVal(:) -! ! type(ContractedGaussian) :: auxContract - -! ! do i = 1, size(contractions) -! ! numContractions = numContractions + contractions(i)%numCartesianOrbital -! ! end do - -! ! call Matrix_constructor(output,int(numContractions,8),int(numContractions,8)) - -! ! a = 1 -! ! b = 1 - -! ! do i=1, size(contractions) - -! ! call ContractedGaussian_product(contractions(i), & -! ! this%gaussianComponents, auxContract) - -! ! do j=1, size(contractions) - -! ! call ContractedGaussian_overlapIntegral( auxContract, contractions(j), auxVal) - -! ! m = 0 -! ! do k = a, auxContract%numCartesianOrbital - 1 -! ! do l = b, contractions(j)%numCartesianOrbital - 1 -! ! m = m + 1 -! ! output%values(k,l)= auxVal(m) -! ! end do -! ! end do -! ! b = b + contractions(j)%numCartesianOrbital -! ! end do -! ! a = a + auxContract%numCartesianOrbital -! ! call ContractedGaussian_destructor(auxContract) -! ! end do - -! ! call Matrix_show(output) - -! ! stop "ExternalPotential_getInteractionMtx" - -! ! ! end function ExternalPotential_getInteractionMtx -! end subroutine ExternalPotential_getInteractionMtx - - - !> - !! @brief Maneja excepciones de la clase - subroutine ExternalPotential_exception( typeMessage, description, debugDescription) - implicit none - integer :: typeMessage - character(*) :: description - character(*) :: debugDescription - - type(Exception) :: ex - - call Exception_constructor( ex , typeMessage ) - call Exception_setDebugDescription( ex, debugDescription ) - call Exception_setDescription( ex, description ) - call Exception_show( ex ) - call Exception_destructor( ex ) - - end subroutine ExternalPotential_exception - -end module ExternalPotential_ diff --git a/src/core/GTFPotential.f90 b/src/core/GTFPotential.f90 new file mode 100644 index 00000000..4ac4b93a --- /dev/null +++ b/src/core/GTFPotential.f90 @@ -0,0 +1,309 @@ +!!****************************************************************************** +!! This code is part of LOWDIN Quantum chemistry package +!! +!! this program has been developed under direction of: +!! +!! Prof. A REYES' Lab. Universidad Nacional de Colombia +!! http://www.qcc.unal.edu.co +!! Prof. R. FLORES' Lab. Universidad de Guadajara +!! http://www.cucei.udg.mx/~robertof +!! +!! Todos los derechos reservados, 2013 +!! +!!****************************************************************************** + +!> @brief This module contains all the routines to handle external and interal GTF potentials +!! @author E. F. Posada (efposadac@unal.edu.co) +!! @version 2.0 +!! Creation data : 06-08-10 +!! +!! History change: +!! +!! - 06-08-10 : Sergio A. Gonzalez ( sergmonic@gmail.com ) +!! -# Creacioon del modulo y metodos basicos +!! - 2011-02-14 : Fernando Posada ( efposadac@unal.edu.co ) +!! -# Reescribe y adapta el módulo para su inclusion en Lowdin +!! - 2024-11-26 : Felix +!! -# Merges ExternalPotential and InternalPotentials modules into a single file (GTFPotential) +module GTFPotential_ + use ContractedGaussian_ + use String_ + use Matrix_ + use Units_ + use Exception_ + implicit none + + type, public :: GaussPot + character(20) :: name + character(50) :: species + character(50) :: otherSpecies + character(50) :: ttype + character(50) :: units + integer :: numOfComponents + integer :: iter + type(ContractedGaussian), allocatable :: gaussianComponents(:) + end type GaussPot + + type, public :: GTFPotential + integer :: ssize + type(GaussPot), allocatable :: potentials(:) + character(50) :: type + logical :: isInstanced + end type GTFPotential + + type(GTFPotential), public, target :: ExternalPotential_instance, InterPotential_instance + +contains + + !> + !! @brief Initializes the class + !! @param this, n + !! @author E. F. Posada, 2013 + subroutine GTFPotential_constructor(this,numberOfPotentials,type) + implicit none + type(GTFPotential) :: this + integer :: numberOfPotentials + character(*) :: type + + + this%ssize = numberOfPotentials + allocate(this%potentials(numberOfPotentials)) + this%isInstanced = .true. + this%type = type + + end subroutine GTFPotential_constructor + + !> + !! @brief Destroys the object + !! @param this + subroutine GTFPotential_destructor(this) + implicit none + type(GTFPotential) :: this + + integer :: i + + do i = 1, this%ssize + if (allocated(this%potentials(i)%gaussianComponents)) deallocate(this%potentials(i)%gaussianComponents) + end do + + if (allocated(this%potentials) ) deallocate(this%potentials) + this%isInstanced=.false. + + end subroutine GTFPotential_destructor + + !> + !! @brief loads information from the input file + !! @param this + !! @author E. F. Posada, 2013 + subroutine GTFPotential_load(this, potId, name, species, otherSpecies) + implicit none + type(GTFPotential) :: this + integer :: potId + character(*) :: name + character(*) :: species + character(*), optional :: otherSpecies + + if(present(otherSpecies)) then + call GaussPot_load(this%potentials(potId), potId, name, species, otherSpecies) + else + call GaussPot_load(this%potentials(potId), potId, name, species) + end if + end subroutine GTFPotential_load + + !> + !! @brief Shows information of the object + !! @param this + subroutine GTFPotential_show(this) + implicit none + type(GTFPotential) :: this + integer :: i, j + + do i=1,this%ssize + if( this%potentials(i)%ttype .eq. "INTERNAL") then + print *,"" + print *,"=======" + write(*,"(A30,A)") "GTF Interparticle potential: ", trim(this%potentials(i)%name) + write(*,"(A4,A10,A5,A10)") "for ", trim(this%potentials(i)%species) ," and ", trim(this%potentials(i)%otherSpecies) + write(*,"(T10,A10,A10)") "Units:", trim(this%potentials(i)%units) + write(*,"(T10,A16,A16)") "Exponent", "Factor" + do j=1,this%potentials(i)%numOfComponents + write(*,"(T10,E16.8,E16.8)") this%potentials(i)%gaussianComponents(j)%orbitalExponents, & + this%potentials(i)%gaussianComponents(j)%contractionCoefficients(1) + end do + else if( this%potentials(i)%ttype .eq. "EXTERNAL") then + print *,"" + print *,"=======" + write(*,"(A25,A20,A5,A10)") "GTF External potential: ", trim(this%potentials(i)%name), " for ", trim(this%potentials(i)%species) + write(*,"(T10,A10,A10)") "Units:", trim(this%potentials(i)%units) + write(*,"(T10,A16,A10,A10,A10,A16)") "Exponent", "R_x", "R_y", "R_z", "Factor" + + do j=1,this%potentials(i)%numOfComponents + write(*,"(T10,E16.8,F10.5,F10.5,F10.5,E16.8)") & + this%potentials(i)%gaussianComponents(j)%orbitalExponents(1), & + this%potentials(i)%gaussianComponents(j)%origin(:), & + this%potentials(i)%gaussianComponents(j)%contractionCoefficients(1) + end do + end if + end do + + end subroutine GTFPotential_show + + !> + !! @brief loads information from the input file + !! @param this + !! @author Felix, 2024 + subroutine GaussPot_load(this, potId, name, species, otherSpecies) + type(GaussPot) :: this + integer :: potId + character(*) :: name + character(*) :: species + character(*), optional :: otherSpecies + + integer :: status, i, j + character(150) :: fileName + character(20) :: token + character(50) :: symbol + logical :: existFile, found + + this%name= trim(name) + this%species= trim(species) + this%units="BOHR" + this%numOfComponents=0 + this%iter=1 + + this%ttype="EXTERNAL" + this%otherSpecies="" + if(present(otherSpecies) ) then + this%ttype="INTERNAL" + this%otherSpecies=otherSpecies + end if + + fileName = trim( trim( CONTROL_instance%DATA_DIRECTORY ) // & + trim(CONTROL_instance%POTENTIALS_DATABASE)// String_getUppercase(trim(name))) + + !! Open Potential file from library + inquire(file=trim(fileName), exist = existFile) + if(existFile) then + open(unit=30, file=trim(fileName), status="old",form="formatted") + else + !! Open Potential file from directory + inquire(file=trim(this%name), exist = existFile) + if(existFile) then + open(unit=30, file=trim(this%name), status="old",form="formatted") + else + !! File not found + call Exception_stopError("The GTFPotential file: "//trim(this%name)//& + " was not found!","GTFPotential module at Load function.") + end if + end if + + !! Open File + rewind(30) + + found = .false. + + !! Open element and Find the proper potential + do while(found .eqv. .false.) + read(30,*, iostat=status) token + symbol = token(3:) + + !! Some debug information in case of error! + if (status > 0 ) call Exception_stopError("ERROR reading InterPotential file: "//trim(this%name)//& + " Please check that file!","GTFPotential module at Load function.") + + if (status == -1 ) call Exception_stopError("The InterPotential: "//trim(this%name)//& + " for: "//trim(species)//trim(otherSpecies)//& + " was not found!","GTFPotential module at Load function.") + + if(trim(token(1:2)) == "O-") then + if(this%ttype=="EXTERNAL" .and. trim(symbol) == trim(species)) found = .true. + if(this%ttype=="INTERNAL" .and. trim(symbol) == trim(species)//trim(otherSpecies)) found = .true. + end if + end do + + !! Neglect any comment + token = "#" + do while(trim(token(1:1)) == "#") + read(30,*) token + end do + + !! Start reading Potential + backspace(30) + read(30,*, iostat=status) this%numOfComponents + + !! Some debug information in case of error! + if (status > 0 ) call Exception_stopError("ERROR reading InternalPotential file: "//trim(this%name)//& + " Please check that file!","GTFPotential module at Load function.") + + allocate(this%gaussianComponents(this%numOfComponents)) + + do i = 1, this%numOfComponents + read(30,*,iostat=status) this%gaussianComponents(i)%id, & + this%gaussianComponents(i)%angularMoment + + if(this%gaussianComponents(i)%angularMoment .gt. 0) then + print *, "Warning! you provided a non-zero angular momentum for a GTFpotential ", this%name ,"this feature is not yet implemented, will be ignored and set to zeo" + this%gaussianComponents(i)%angularMoment=0 + end if + + this%gaussianComponents(i)%length = 1 + + !! Some debug information in case of error! + if (status > 0 ) call Exception_stopError("ERROR reading InternalPotential file: "//trim(this%name)//& + " Please check that file!","GTFPotential module at Load function.") + + allocate(this%gaussianComponents(i)%orbitalExponents(this%gaussianComponents(i)%length)) + allocate(this%gaussianComponents(i)%contractionCoefficients(this%gaussianComponents(i)%length)) + + do j = 1, this%gaussianComponents(i)%length + + read(30,*,iostat=status) this%gaussianComponents(i)%orbitalExponents(j), & + this%gaussianComponents(i)%contractionCoefficients(j) + read(30,*,iostat=status) this%gaussianComponents(i)%origin + + !! Some debug information in case of error! + if (status > 0 ) call Exception_stopError("ERROR reading InternalPotential file: "//trim(this%name)//& + " Please check that file!","GTFPotential module at Load function.") + + end do + + if(this%ttype=="INTERNAL" .and. sum(this%gaussianComponents(i)%origin(:)**2) .gt. CONTROL_instance%DOUBLE_ZERO_THRESHOLD) then + print *, "Warning! you provided a non-zero origin for interpotential ", this%name ,"this feature is not yet implemented, will be ignored and set to zeo" + this%gaussianComponents(i)%origin=0.0_8 + end if + + !! Calculates the number of Cartesian orbitals, by dimensionality + select case(CONTROL_instance%DIMENSIONALITY) + case(3) + this%gaussianComponents(i)%numCartesianOrbital = ( ( this%gaussianComponents(i)%angularMoment + 1_8 )*( this%gaussianComponents(i)%angularMoment + 2_8 ) ) / 2_8 + case(2) + this%gaussianComponents(i)%numCartesianOrbital = ( ( this%gaussianComponents(i)%angularMoment + 1_8 ) ) + case(1) + this%gaussianComponents(i)%numCartesianOrbital = 1 + case default + call Exception_stopError("Class object InternalPotential in load function",& + "This Dimensionality is not available") + end select + + !! Normalize + allocate(this%gaussianComponents(i)%contNormalization(this%gaussianComponents(i)%numCartesianOrbital)) + allocate(this%gaussianComponents(i)%primNormalization(this%gaussianComponents(i)%length, & + this%gaussianComponents(i)%length*this%gaussianComponents(i)%numCartesianOrbital)) + + this%gaussianComponents(i)%contNormalization = 1.0_8 + this%gaussianComponents(i)%primNormalization = 1.0_8 + + call ContractedGaussian_normalizePrimitive(this%gaussianComponents(i)) + call ContractedGaussian_normalizeContraction(this%gaussianComponents(i)) + + !! DEBUG + ! call ContractedGaussian_showInCompactForm(InterPotential_instance%potentials(potId)%gaussianComponents(i)) + + end do + + close(30) + + !!DONE + end subroutine GaussPot_load + +end module GTFPotential_ diff --git a/src/core/InputManager.f90 b/src/core/InputManager.f90 index e5a0abde..506a89fb 100644 --- a/src/core/InputManager.f90 +++ b/src/core/InputManager.f90 @@ -23,8 +23,7 @@ module InputManager_ use Exception_ use Particle_ use MolecularSystem_ - use InterPotential_ - use ExternalPotential_ + use GTFPotential_ implicit none @@ -344,6 +343,7 @@ subroutine InputManager_loadGeometry() real(8):: InputParticle_origin(3) real(8) :: InputParticle_charge real(8) :: InputParticle_mass + integer :: InputParticle_eta real(8) :: InputParticle_omega character(15):: InputParticle_qdoCenterOf character(3):: InputParticle_fixedCoordinates @@ -359,6 +359,7 @@ subroutine InputManager_loadGeometry() InputParticle_basisSetName, & InputParticle_charge, & InputParticle_mass, & + InputParticle_eta, & InputParticle_omega, & InputParticle_qdoCenterOf, & InputParticle_origin, & @@ -538,6 +539,7 @@ subroutine InputManager_loadGeometry() InputParticle_basisSetName = "NONE" InputParticle_charge=0.0_8 InputParticle_mass=0.0_8 + InputParticle_eta=0 InputParticle_omega=0.0_8 InputParticle_qdoCenterOf = "NONE" InputParticle_origin=0.0_8 @@ -620,7 +622,8 @@ subroutine InputManager_loadGeometry() spin="ALPHA", & id = particlesID(speciesID), & charge = InputParticle_charge, & - mass = InputParticle_mass, & + mass = InputParticle_mass, & + eta = InputParticle_eta, & omega = InputParticle_omega ) !!BETA SET @@ -644,6 +647,7 @@ subroutine InputManager_loadGeometry() id = particlesID(speciesID), & charge = InputParticle_charge, & mass = InputParticle_mass, & + eta = InputParticle_eta, & omega = InputParticle_omega ) else @@ -671,6 +675,7 @@ subroutine InputManager_loadGeometry() id = particlesID(speciesID), & charge = InputParticle_charge, & mass = InputParticle_mass, & + eta = InputParticle_eta, & omega = InputParticle_omega ) end if @@ -750,7 +755,7 @@ subroutine InputManager_loadPotentials() ! Load interpotentials if(CONTROL_instance%IS_THERE_INTERPARTICLE_POTENTIAL) then - call InterPotential_constructor(Input_instance%numberOfInterPots) + call GTFPotential_constructor(InterPotential_instance, Input_instance%numberOfInterPots,"INTERNAL") !! Reload input file rewind(4) @@ -760,10 +765,10 @@ subroutine InputManager_loadPotentials() read(4,NML=InterPot, iostat=stat) if( stat > 0 ) then - call InputManager_exception( ERROR, "check the TASKS block in your input file", "InputManager loadTask function" ) + call InputManager_exception( ERROR, "check the INTERPOTENTIAL block in your input file", "InputManager loadTask function" ) end if - call InterPotential_load(potId, trim(InterPot_name), trim(InterPot_specie), trim(InterPot_otherSpecie)) + call GTFPotential_load(InterPotential_instance, potId, trim(InterPot_name), trim(InterPot_specie), trim(InterPot_otherSpecie)) end do @@ -772,7 +777,7 @@ subroutine InputManager_loadPotentials() ! Load External Potentials if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) then - call ExternalPotential_constructor(Input_instance%numberOfExternalPots) + call GTFPotential_constructor(ExternalPotential_instance, Input_instance%numberOfExternalPots,"EXTERNAL") !! Reload input file rewind(4) @@ -782,10 +787,10 @@ subroutine InputManager_loadPotentials() read(4,NML=ExternalPot, iostat=stat) if( stat > 0 ) then - call InputManager_exception( ERROR, "check the TASKS block in your input file", "InputManager loadTask function" ) + call InputManager_exception( ERROR, "check the EXTERPOTENTIAL block in your input file", "InputManager loadTask function" ) end if - call ExternalPotential_load(potId, trim(ExternalPot_name), trim(ExternalPot_specie)) + call GTFPotential_load(ExternalPotential_instance, potId, trim(ExternalPot_name), trim(ExternalPot_specie)) end do end if diff --git a/src/core/InterPotential.f90 b/src/core/InterPotential.f90 deleted file mode 100644 index f30008e7..00000000 --- a/src/core/InterPotential.f90 +++ /dev/null @@ -1,298 +0,0 @@ -!!****************************************************************************** -!! This code is part of LOWDIN Quantum chemistry package -!! -!! this program has been developed under direction of: -!! -!! Prof. A REYES' Lab. Universidad Nacional de Colombia -!! http://www.qcc.unal.edu.co -!! Prof. R. FLORES' Lab. Universidad de Guadajara -!! http://www.cucei.udg.mx/~robertof -!! -!! Todos los derechos reservados, 2013 -!! -!!****************************************************************************** - -!> @brief This module contains all the routines to handle inter particle potentials -!! @author E. F. Posada (efposadac@unal.edu.co) -!! @version 2.0 - -module InterPotential_ - use ContractedGaussian_ - use String_ - use Exception_ - implicit none - - type, public :: InterPot - character(20) :: name - character(50) :: specie - character(50) :: otherSpecie - character(50) :: ttype - character(50) :: units - integer :: numOfComponents - integer :: iter - type(ContractedGaussian), allocatable :: gaussianComponents(:) - end type - - type, public :: InterPotential - integer :: ssize - type(InterPot), allocatable :: Potentials(:) - logical :: isInstanced - end type - - type(InterPotential), public, target :: InterPotential_instance - -contains - - !> - !! @brief Initializes the class - !! @param this - !! @author E. F. Posada, 2013 - subroutine InterPotential_constructor(numberOfPotentials) - implicit none - integer :: numberOfPotentials - - InterPotential_instance%ssize = numberOfPotentials - allocate(InterPotential_instance%potentials(numberOfPotentials)) - InterPotential_instance%isInstanced = .true. - - end subroutine InterPotential_constructor - - !> - !! @brief destroy the class - !! @param this - !! @author E. F. Posada, 2013 - subroutine InterPotential_destructor() - implicit none - - integer :: i - - do i = 1, InterPotential_instance%ssize - if (allocated(InterPotential_instance%potentials(i)%gaussianComponents) ) deallocate(InterPotential_instance%potentials(i)%gaussianComponents) - end do - - if (allocated(InterPotential_instance%potentials) )deallocate(InterPotential_instance%potentials) - - end subroutine InterPotential_destructor - - !> - !! @brief Shows information of the object - !! @param this - subroutine InterPotential_show() - implicit none - integer :: i, j - type(InterPot), pointer :: this - - do i=1,InterPotential_instance%ssize - this => InterPotential_instance%potentials(i) - print *,"" - print *,"=======" - print *, "InterParticle Potential for ", trim(this%specie) ," and ", trim(this%otherSpecie), " : ", trim(this%name) - print *, "Type : ", trim(this%ttype) - write(6,"(T10,A20,A10,A10,A10,A20)") "Exponent", "l", "Factor" - do j=1,this%numOfComponents - write(6,"(T10,F20.15,I10,F20.15)") this%gaussianComponents(j)%orbitalExponents, & - this%gaussianComponents(j)%angularMoment, this%gaussianComponents(j)%contractionCoefficients(1) - end do - end do - - end subroutine InterPotential_show - - - !> - !! @brief loads information from the input file - !! @param this - !! @author E. F. Posada, 2015 - subroutine InterPotential_load(potId, name, species, otherSpecies) - implicit none - integer :: potId - character(*) :: name - character(*) :: species - character(*) :: otherSpecies - - type(InterPot), pointer :: this - integer :: status, i, j - character(150) :: fileName - character(20) :: token - character(20) :: symbol - logical :: existFile, found - - this => InterPotential_instance%potentials(potId) - - this%name= trim(name) - this%specie= trim(species) - this%otherSpecie= trim(otherSpecies) - this%ttype="" - this%units="bohr" - this%numOfComponents=0 - this%iter=1 - - fileName = trim( trim( CONTROL_instance%DATA_DIRECTORY ) // & - trim(CONTROL_instance%POTENTIALS_DATABASE)// String_getUppercase(trim(name))) - - - inquire(file=trim(fileName), exist = existFile) - if(existFile) then - - !! Open File - open(unit=30, file=trim(fileName), status="old",form="formatted") - rewind(30) - - found = .false. - - !! Open element and Find the proper potential - do while(found .eqv. .false.) - read(30,*, iostat=status) token - symbol = token(3:) - - !! Some debug information in case of error! - if (status > 0 ) then - - call InterPotential_exception(ERROR, & - "ERROR reading InterPotential file: "//trim(this%name)//& - " Please check that file!","InternalPotential module at Load function.") - - end if - - if (status == -1 ) then - - call InterPotential_exception(ERROR, & - "The InterPotential: "//trim(this%name)//& - " for: "//trim(species)//trim(otherSpecies)//& - " was not found!","InternalPotential module at Load function.") - - end if - - if(trim(token(1:2)) == "O-") then - if(trim(symbol) == trim(species)//trim(otherSpecies)) then - found = .true. - - end if - - end if - - end do - - !! Neglect any comment - token = "#" - do while(trim(token(1:1)) == "#") - - read(30,*) token - - end do - - !! Start reading Potential - backspace(30) - - read(30,*, iostat=status) this%numOfComponents - - !! Some debug information in case of error! - if (status > 0 ) then - - call InterPotential_exception(ERROR, & - "ERROR reading InternalPotential file: "//trim(this%name)//& - " Please check that file!","InternalPotential module at Load function.") - - end if - - allocate(this%gaussianComponents(this%numOfComponents)) - - do i = 1, this%numOfComponents - - read(30,*,iostat=status) this%gaussianComponents(i)%id, & - this%gaussianComponents(i)%angularMoment - this%gaussianComponents(i)%length = 1 - - !! Some debug information in case of error! - if (status > 0 ) then - - call InterPotential_exception(ERROR, & - "ERROR reading InternalPotential file: "//trim(this%name)//& - " Please check that file!","InternalPotential module at Load function.") - - end if - - allocate(this%gaussianComponents(i)%orbitalExponents(this%gaussianComponents(i)%length)) - allocate(this%gaussianComponents(i)%contractionCoefficients(this%gaussianComponents(i)%length)) - - do j = 1, this%gaussianComponents(i)%length - - read(30,*,iostat=status) this%gaussianComponents(i)%orbitalExponents(j), & - this%gaussianComponents(i)%contractionCoefficients(j) - read(30,*,iostat=status) this%gaussianComponents(i)%origin - - !! Some debug information in case of error! - if (status > 0 ) then - - call InterPotential_exception(ERROR, & - "ERROR reading InternalPotential file: "//trim(this%name)//& - " Please check that file!","InternalPotential module at Load function.") - - end if - - end do - - - !! Calculates the number of Cartesian orbitals, by dimensionality - select case(CONTROL_instance%DIMENSIONALITY) - case(3) - this%gaussianComponents(i)%numCartesianOrbital = ( ( this%gaussianComponents(i)%angularMoment + 1_8 )*( this%gaussianComponents(i)%angularMoment + 2_8 ) ) / 2_8 - case(2) - this%gaussianComponents(i)%numCartesianOrbital = ( ( this%gaussianComponents(i)%angularMoment + 1_8 ) ) - case(1) - this%gaussianComponents(i)%numCartesianOrbital = 1 - case default - call InterPotential_exception( ERROR, & - "Class object InternalPotential in load function",& - "This Dimensionality is not available") - end select - - !! Normalize - allocate(this%gaussianComponents(i)%contNormalization(this%gaussianComponents(i)%numCartesianOrbital)) - allocate(this%gaussianComponents(i)%primNormalization(this%gaussianComponents(i)%length, & - this%gaussianComponents(i)%length*this%gaussianComponents(i)%numCartesianOrbital)) - - this%gaussianComponents(i)%contNormalization = 1.0_8 - this%gaussianComponents(i)%primNormalization = 1.0_8 - - call ContractedGaussian_normalizePrimitive(this%gaussianComponents(i)) - call ContractedGaussian_normalizeContraction(this%gaussianComponents(i)) - - !! DEBUG - ! call ContractedGaussian_showInCompactForm(InterPotential_instance%potentials(potId)%gaussianComponents(i)) - - end do - - close(30) - - !!DONE - - else - - call InterPotential_exception(ERROR, & - "The InternalPotential file: "//trim(name)//& - " was not found!","InternalPotential module at Load function.") - - end if - - end subroutine InterPotential_load - - - !> - !! @brief Handles class exceptions - !< - subroutine InterPotential_exception( typeMessage, description, debugDescription) - implicit none - integer :: typeMessage - character(*) :: description - character(*) :: debugDescription - - type(Exception) :: ex - - call Exception_constructor( ex , typeMessage ) - call Exception_setDebugDescription( ex, debugDescription ) - call Exception_setDescription( ex, description ) - call Exception_show( ex ) - call Exception_destructor( ex ) - - end subroutine InterPotential_exception -end module InterPotential_ diff --git a/src/core/MolecularSystem.f90 b/src/core/MolecularSystem.f90 index b5a766fe..32dced1b 100644 --- a/src/core/MolecularSystem.f90 +++ b/src/core/MolecularSystem.f90 @@ -38,8 +38,7 @@ module MolecularSystem_ use Matrix_ use Vector_ use InternalCoordinates_ - use ExternalPotential_ - use InterPotential_ + use GTFPotential_ implicit none type , public :: MolecularSystem @@ -343,9 +342,10 @@ subroutine MolecularSystem_showInformation(this) print *," MOLECULAR SYSTEM: ",trim(system%name) print *,"-----------------" print *,"" - write (6,"(T5,A16,A)") "DESCRIPTION : ", trim( system%description ) - write (6,"(T5,A16,I3)") "CHARGE : ",system%charge - write (6,"(T5,A16,A4)") "PUNTUAL GROUP : ", "NONE" + write (6,"(T5,A16,A)") "DESCRIPTION : ", trim( system%description ) + write (6,"(T5,A16,I3)") "CHARGE : ",system%charge + write (6,"(T5,A16,F12.4)") "MASS (m_e) : ", MolecularSystem_getTotalMass(system) + write (6,"(T5,A16,A4)") "PUNTUAL GROUP : ", "NONE" print *,"" @@ -382,9 +382,9 @@ subroutine MolecularSystem_showParticlesInformation(this) !! print *,"" print *," INFORMATION OF QUANTUM SPECIES " - write (6,"(T5,A70)") "---------------------------------------------------------------------------------------------" + write (6,"(T5,A85)") "------------------------------------------------------------------------------------------------------------" write (6,"(T10,A2,A4,A8,A12,A4,A5,A6,A5,A6,A5,A4,A5,A12)") "ID", " ","Symbol", " ","mass", " ","charge", " ","omega","","spin","","multiplicity" - write (6,"(T5,A70)") "---------------------------------------------------------------------------------------------" + write (6,"(T5,A85)") "------------------------------------------------------------------------------------------------------------" do i = 1, system%numberOfQuantumSpecies write (6,'(T8,I3.0,A5,A10,A5,F10.4,A5,F5.2,A5,F5.2,A5,F5.2,A5,F5.2)') & @@ -415,16 +415,16 @@ subroutine MolecularSystem_showParticlesInformation(this) print *,"" print *," BASIS SET FOR SPECIES " - write (6,"(T7,A60)") "------------------------------------------------------------" - write (6,"(T10,A8,A10,A8,A5,A12,A5,A9)") "Symbol", " ","N. Basis", " ","N. Particles"," ","Basis Set" - write (6,"(T7,A60)") "------------------------------------------------------------" + write (6,"(T7,A70)") "----------------------------------------------------------------------" + write (6,"(T10,A8,A10,A8,A5,A12,A5,A20)") "Symbol", " ","N. Basis", " ","N. Particles"," ","Basis Set" + write (6,"(T7,A70)") "----------------------------------------------------------------------" !! Only shows the basis-set name of the first particle by specie. do i = 1, system%numberOfQuantumSpecies if( system%species(i)%isElectron .and. CONTROL_instance%IS_OPEN_SHELL ) then - write (6,'(T10,A10,A5,I8,A5,I12,A5,A10)') & + write (6,'(T10,A10,A5,I8,A5,I12,A5,A20)') & trim(system%species(i)%symbol)," ",& !MolecularSystem_getTotalNumberOfContractions(i)," ",& system%species(i)%basisSetSize," ",& @@ -434,7 +434,7 @@ subroutine MolecularSystem_showParticlesInformation(this) else - write (6,'(T10,A10,A5,I8,A5,I12,A5,A10)') & + write (6,'(T10,A10,A5,I8,A5,I12,A5,A20)') & trim(system%species(i)%symbol)," ",& !MolecularSystem_getTotalNumberOfContractions(i)," ",& system%species(i)%basisSetSize," ",& @@ -501,7 +501,7 @@ subroutine MolecularSystem_showParticlesInformation(this) if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) then print *,"" print *," INFORMATION OF EXTERNAL POTENTIALS " - call ExternalPotential_show() + call GTFPotential_show(ExternalPotential_instance) print *,"" print *," END INFORMATION OF EXTERNAL POTENTIALS" print *,"" @@ -510,7 +510,7 @@ subroutine MolecularSystem_showParticlesInformation(this) if(CONTROL_instance%IS_THERE_INTERPARTICLE_POTENTIAL) then print *,"" print *," INFORMATION OF INTER-PARTICLE POTENTIALS " - call InterPotential_show() + call GTFPotential_show(InterPotential_instance) print *,"" print *," END INFORMATION OF INTER-PARTICLE POTENTIALS" print *,"" @@ -671,7 +671,7 @@ subroutine MolecularSystem_saveToFile(targetFilePrefix) do i = 1, ExternalPotential_instance%ssize write(40,*) i write(40,*) ExternalPotential_instance%potentials(i)%name - write(40,*) ExternalPotential_instance%potentials(i)%specie + write(40,*) ExternalPotential_instance%potentials(i)%species end do end if @@ -681,8 +681,8 @@ subroutine MolecularSystem_saveToFile(targetFilePrefix) do i = 1, InterPotential_instance%ssize write(40,*) i write(40,*) InterPotential_instance%potentials(i)%name - write(40,*) InterPotential_instance%potentials(i)%specie - write(40,*) InterPotential_instance%potentials(i)%otherSpecie + write(40,*) InterPotential_instance%potentials(i)%species + write(40,*) InterPotential_instance%potentials(i)%otherSpecies end do end if @@ -950,16 +950,14 @@ subroutine MolecularSystem_loadFromFile( form, targetPrefix ) if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) then read(40,*) auxValue - call ExternalPotential_constructor(auxValue) - - !! FELIX TODO: create function to get potential ID + call GTFPotential_constructor(ExternalPotential_instance,auxValue,"EXTERNAL") do j = 1, ExternalPotential_instance%ssize read(40,*) i read(40,*) name read(40,*) species - call ExternalPotential_load(i, name, species) + call GTFPotential_load(ExternalPotential_instance, i, name, species) end do @@ -968,7 +966,7 @@ subroutine MolecularSystem_loadFromFile( form, targetPrefix ) if(CONTROL_instance%IS_THERE_INTERPARTICLE_POTENTIAL) then read(40,*) auxValue - call InterPotential_constructor(auxValue) + call GTFPotential_constructor(InterPotential_instance,auxValue,"INTERNAL") do j = 1, InterPotential_instance%ssize read(40,*) i @@ -976,7 +974,7 @@ subroutine MolecularSystem_loadFromFile( form, targetPrefix ) read(40,*) species read(40,*) otherSpecies - call InterPotential_load(i, name, species, otherSpecies) + call GTFPotential_load(InterPotential_instance, i, name, species, otherSpecies) end do @@ -1979,13 +1977,14 @@ end subroutine MolecularSystem_changeOrbitalOrder !> !! @brief Lee la matriz de densidad y los orbitales de un archivo fchk tipo Gaussian - subroutine MolecularSystem_readFchk( fileName, coefficients, densityMatrix, nameOfSpecies ) + subroutine MolecularSystem_readFchk( fileName, coefficients, densityMatrix, nameOfSpecies, readSuccess ) implicit none character(*), intent(in) :: fileName type(Matrix), intent(inout) :: coefficients type(Matrix), intent(inout) :: densityMatrix character(*) :: nameOfSpecies + logical, optional :: readSuccess integer :: speciesID integer :: numberOfContractions @@ -2006,9 +2005,15 @@ subroutine MolecularSystem_readFchk( fileName, coefficients, densityMatrix, name speciesID=MolecularSystem_getSpecieID(nameOfSpecies) numberOfContractions=MolecularSystem_getTotalnumberOfContractions( speciesID ) inquire(FILE = trim(fileName), EXIST = existFchk ) - if ( .not. existFchk ) call MolecularSystem_exception( ERROR, "I did not find any .fchk coefficients file", "At MolecularSystem_readFchk") + if ( .not. existFchk .and. present(readSuccess)) then + readSuccess=.false. + call MolecularSystem_exception( WARNING, "I did not find the "//trim(filename)//" coefficients file for "//trim(nameOfSpecies), "At MolecularSystem_readFchk") + return + end if + if ( .not. existFchk) then + call MolecularSystem_exception( ERROR, "I did not find the "//trim(filename)//" coefficients file for "//trim(nameOfSpecies), "At MolecularSystem_readFchk") + end if - fchkUnit = 50 open(unit=fchkUnit, file=filename, status="old", form="formatted", access='sequential', action='read') @@ -2097,8 +2102,8 @@ subroutine MolecularSystem_readFchk( fileName, coefficients, densityMatrix, name ! print *, "density matrix from orbitals read" ! call Matrix_show(densityMatrix) - close(fchkUnit) + if(present(readSuccess)) readSuccess=.true. end subroutine MolecularSystem_readFchk @@ -2783,6 +2788,9 @@ function MolecularSystem_getTotalMass( this, unid ) result( output ) case("AMU") output = output * AMU + case("DALTON") + output = output * DALTON + case default end select diff --git a/src/core/Particle.f90 b/src/core/Particle.f90 index 7a845404..f927d006 100644 --- a/src/core/Particle.f90 +++ b/src/core/Particle.f90 @@ -32,6 +32,7 @@ module Particle_ use Exception_ use CONTROL_ + use Units_ use PhysicalConstants_ use AtomicElement_ use ElementalParticle_ @@ -49,6 +50,7 @@ module Particle_ real(8) :: origin(3) !< Posicion espacial real(8) :: charge !< Carga asociada a la particula. real(8) :: mass !< Masa asociada a la particula. + real(8) :: eta !< Particles per orbital real(8) :: omega !< harmonic oscillator frequency real(8) :: spin !< Especifica el espin de la particula real(8) :: totalCharge !< Carga total asociada a la particula. @@ -83,7 +85,7 @@ module Particle_ !! -Re-written and Verified, 2013. E. F. Posada !! @version 2.0 subroutine Particle_load( this, name, baseName, origin, fix, multiplicity, & - addParticles, subsystem, translationCenter, rotationPoint, rotateAround, spin, id, charge, mass, omega, qdoCenterOf ) + addParticles, subsystem, translationCenter, rotationPoint, rotateAround, spin, id, charge, mass, eta, omega, qdoCenterOf ) implicit none type(particle) :: this character(*), intent(in) :: name @@ -101,6 +103,7 @@ subroutine Particle_load( this, name, baseName, origin, fix, multiplicity, & integer, intent(in) :: id real(8), intent(in), optional :: charge real(8), intent(in), optional :: mass + integer, intent(in), optional :: eta real(8), intent(in), optional :: omega type(AtomicElement) :: element @@ -120,6 +123,7 @@ subroutine Particle_load( this, name, baseName, origin, fix, multiplicity, & real(8) :: auxOrigin(3) real(8) :: auxCharge real(8) :: auxMass + integer :: auxEta real(8) :: auxOmega real(8) :: auxMultiplicity logical :: isDummy @@ -137,6 +141,9 @@ subroutine Particle_load( this, name, baseName, origin, fix, multiplicity, & auxMass=0.0_8 if ( present(mass) ) auxMass= mass + auxEta=0 + if ( present(eta) ) auxEta= eta + auxOmega=0.0_8 if ( present(omega) ) auxOmega= omega @@ -219,6 +226,7 @@ subroutine Particle_load( this, name, baseName, origin, fix, multiplicity, & origin=auxOrigin, & charge=auxCharge, & mass=auxMass, & + eta=auxEta, & omega=auxOmega, & basisSetName=trim(baseName), & elementSymbol=trim(elementSymbol), & @@ -309,6 +317,7 @@ subroutine Particle_load( this, name, baseName, origin, fix, multiplicity, & origin = auxOrigin,& charge = auxCharge,& mass = auxMass,& + eta=auxEta, & omega=auxOmega, & basisSetName = trim(baseName), & elementSymbol = trim(elementSymbol), & @@ -330,10 +339,10 @@ subroutine Particle_load( this, name, baseName, origin, fix, multiplicity, & this%charge = element%atomicNumber end if - if ( auxMass == 0.0_8) then - this%mass = element%massicNumber * PhysicalConstants_NEUTRON_MASS & + if ( auxMass .eq. 0.0_8 .and. element%atomicWeight .eq. 0.0_8 ) & + this%mass = element%massicNumber * PhysicalConstants_NEUTRON_MASS & + element%atomicNumber * (PhysicalConstants_PROTON_MASS - PhysicalConstants_NEUTRON_MASS) - end if + if ( auxMass .eq. 0.0_8 ) this%mass = element%atomicWeight*DALTON - element%atomicNumber*PhysicalConstants_ELECTRON_MASS this%totalCharge = element%atomicNumber this%internalSize = 1 + auxAdditionOfParticles @@ -393,8 +402,14 @@ subroutine Particle_load( this, name, baseName, origin, fix, multiplicity, & this%charge = element%atomicNumber this%totalCharge = element%atomicNumber end if - this%mass = element%massicNumber * PhysicalConstants_NEUTRON_MASS & - + element%atomicNumber * (PhysicalConstants_PROTON_MASS - PhysicalConstants_NEUTRON_MASS) + + + if ( element%atomicWeight .eq. 0.0_8 ) then + this%mass = element%massicNumber * PhysicalConstants_NEUTRON_MASS & + + element%atomicNumber * (PhysicalConstants_PROTON_MASS - PhysicalConstants_NEUTRON_MASS) + else + this%mass = element%atomicWeight*DALTON - element%atomicNumber*PhysicalConstants_ELECTRON_MASS + end if this%internalSize = 1 + auxAdditionOfParticles this%id = id @@ -439,8 +454,14 @@ subroutine Particle_load( this, name, baseName, origin, fix, multiplicity, & if ( auxCharge == 0.0_8) then this%charge = element%atomicNumber end if - this%mass = element%massicNumber * PhysicalConstants_NEUTRON_MASS & - + element%atomicNumber * (PhysicalConstants_PROTON_MASS - PhysicalConstants_NEUTRON_MASS) + + if ( element%atomicWeight .eq. 0.0_8 ) then + this%mass = element%massicNumber * PhysicalConstants_NEUTRON_MASS & + + element%atomicNumber * (PhysicalConstants_PROTON_MASS - PhysicalConstants_NEUTRON_MASS) + else + this%mass = element%atomicWeight*DALTON - element%atomicNumber*PhysicalConstants_ELECTRON_MASS + end if + this%totalCharge = element%atomicNumber this%internalSize = 1 + auxAdditionOfParticles this%id = id @@ -468,12 +489,19 @@ subroutine Particle_load( this, name, baseName, origin, fix, multiplicity, & else if ( present(baseName) .and. trim(baseName) /= "DIRAC" .and. trim(baseName) /= "MM") then call ElementalParticle_load( eparticle, trim(name) ) + + if(eparticle%custom .and. auxMass .ne. 0.0_8 .and. auxCharge .ne. 0.0_8) & + print *, "I'm loading a custom particle from the input", trim(name), "with mass", auxMass, "and charge", auxCharge + if(eparticle%custom .and. auxMass .eq. 0.0_8 .and. auxCharge .eq. 0.0_8) & + call Particle_exception( ERROR, "Elemental particle: "//trim(name)//& + " NOT found in ElementalParticles.lib. If you want to use a custom particle, provide at least its mass and charge in the input", "In Particle_load routine") call Particle_build( this = this, & isQuantum = .true., & origin = auxOrigin, & charge = auxCharge, & mass = auxMass, & + eta=auxEta, & omega = auxOmega, & basisSetName = trim(baseName), & elementSymbol = trim(elementSymbol), & @@ -586,7 +614,7 @@ end subroutine Particle_load !! @param this Quantum particle or point charge !! @author S. A. Gonzalez (before known as Particle_constructor) subroutine Particle_build( this, name, symbol, basisSetName, elementSymbol, nickname, & - origin, mass, omega, qdoCenterOf, charge, totalCharge, spin, & + origin, mass, eta, omega, qdoCenterOf, charge, totalCharge, spin, & owner, subsystem, translationCenter, rotationPoint, rotateAround, massNumber, isQuantum, isDummy) implicit none @@ -600,6 +628,7 @@ subroutine Particle_build( this, name, symbol, basisSetName, elementSymbol, nick character(*), optional, intent(in) :: qdoCenterOf real(8), optional, intent(in) :: origin(3) real(8), optional, intent(in) :: mass + integer, optional, intent(in) :: eta real(8), optional, intent(in) :: omega real(8), optional, intent(in) :: charge real(8), optional, intent(in) :: totalCharge @@ -621,6 +650,7 @@ subroutine Particle_build( this, name, symbol, basisSetName, elementSymbol, nick this%isQuantum = .false. this%origin = 0.0_8 this%mass = PhysicalConstants_ELECTRON_MASS + this%eta = 0 this%omega = 0.0_8 this%qdoCenterOf = "NONE" this%charge = PhysicalConstants_ELECTRON_CHARGE @@ -646,6 +676,7 @@ subroutine Particle_build( this, name, symbol, basisSetName, elementSymbol, nick !! Loads optional information if ( present(origin) ) this%origin = origin if ( present( mass ) ) this%mass = mass + if ( present( eta ) ) this%eta = eta if ( present( omega ) ) this%omega = omega if ( present( qdoCenterOf ) ) this%qdoCenterOf = qdoCenterOf if ( present(charge) ) this%charge = charge @@ -722,6 +753,7 @@ subroutine Particle_destroy( this ) this%origin = 0.0_8 this%charge = 0 this%mass = 0.0_8 + this%eta = 0 this%omega = 0.0_8 this%spin = 0.0_8 this%totalCharge = 0 @@ -760,6 +792,7 @@ subroutine Particle_show( this ) write (6,"(T10,A16,I8)") "Owner : ",this%owner write (6,"(T10,A16,F8.2)") "Charge : ",this%charge write (6,"(T10,A16,F8.2)") "Mass : ",this%mass + write (6,"(T10,A16,I8)") "Eta : ",this%eta write (6,"(T10,A16,F8.2)") "Omega : ",this%omega write (6,"(T10,A16,F8.2)") "QDO center of : ",this%qdoCenterOf write (6,"(T10,A16,F8.2)") "Spin : ",this%spin @@ -828,6 +861,7 @@ subroutine Particle_saveToFile( this, unit ) write(unit,*) this%origin write(unit,*) this%charge write(unit,*) this%mass + write(unit,*) this%eta write(unit,*) this%omega write(unit,*) this%qdoCenterOf write(unit,*) this%spin @@ -888,6 +922,7 @@ subroutine Particle_loadFromFile( this, unit ) read(unit,*) this%origin read(unit,*) this%charge read(unit,*) this%mass + read(unit,*) this%eta read(unit,*) this%omega read(unit,*) this%qdoCenterOf read(unit,*) this%spin diff --git a/src/core/Species.f90 b/src/core/Species.f90 index 061113e8..e4d0de4c 100644 --- a/src/core/Species.f90 +++ b/src/core/Species.f90 @@ -93,9 +93,9 @@ subroutine Species_setSpecie( this, speciesID) if(trim(this%statistics) == "FERMION") then this%kappa = -1.0_8 - this%eta = 2.0_8 - this%lambda = 2.0_8 - this%particlesFraction = 0.5_8 + this%eta = 1.0_8 + this%lambda = 1.0_8 + this%particlesFraction = 1.0_8 else @@ -151,6 +151,13 @@ subroutine Species_setSpecie( this, speciesID) !! Adjust multiplicity this%multiplicity = this%multiplicity + 1 + !! Adjust eta + if(this%particles(1)%eta .ne. 0) then + this%eta = this%particles(1)%eta + this%lambda = this%particles(1)%eta + this%particlesFraction = 1.0_8/this%particles(1)%eta + end if + !! Adjust Occupation number this%ocupationNumber = this%ocupationNumber * this%particlesFraction diff --git a/src/core/Units.f90 b/src/core/Units.f90 index 3a51a81d..7f111dde 100644 --- a/src/core/Units.f90 +++ b/src/core/Units.f90 @@ -38,6 +38,7 @@ module Units_ real(8) , parameter :: DEBYES = 2.541764_8 real(8) , parameter :: ELECTRON_REST = 1.0_8 real(8) , parameter :: AMU = 1.0_8/1822.88850065855_8 + real(8) , parameter :: DALTON = 1822.88850065855_8 real(8) , parameter :: kg = 9.109382616D-31 real(8) , parameter :: DEGREES = 57.29577951_8 real(8) , parameter :: CM_NEG1 = 219476.0_8 diff --git a/src/ints/DirectIntegralManager.f90 b/src/ints/DirectIntegralManager.f90 index 4733500f..6d98d24c 100644 --- a/src/ints/DirectIntegralManager.f90 +++ b/src/ints/DirectIntegralManager.f90 @@ -32,7 +32,7 @@ module DirectIntegralManager_ use RysQuadrature_ use Matrix_ use Stopwatch_ - use ExternalPotential_ + use GTFPotential_ use String_ !# use RysQInts_ !! Please do not remove this line @@ -866,8 +866,8 @@ subroutine DirectIntegralManager_getExternalPotentialIntegrals(molSystem,species do i= 1, ExternalPotential_instance%ssize !if( trim(potential(i)%specie)==trim(interactNameSelected) ) then ! This does not work for UHF ! if ( String_findSubstring(trim( molSystem%species(speciesID)%name ), & - ! trim(String_getUpperCase(trim(ExternalPotential_instance%potentials(i)%specie)))) == 1 ) then - if ( trim( molSystem%species(speciesID)%symbol) == trim(String_getUpperCase(trim(ExternalPotential_instance%potentials(i)%specie))) ) then + ! trim(String_getUpperCase(trim(ExternalPotential_instance%potentials(i)%species)))) == 1 ) then + if ( trim( molSystem%species(speciesID)%symbol) == trim(String_getUpperCase(trim(ExternalPotential_instance%potentials(i)%species))) ) then potID=i exit end if diff --git a/src/ints/G12Integrals.f90 b/src/ints/G12Integrals.f90 index 19446252..1bfc7d4b 100644 --- a/src/ints/G12Integrals.f90 +++ b/src/ints/G12Integrals.f90 @@ -29,7 +29,7 @@ module G12Integrals_ use CONTROL_ use MolecularSystem_ use ContractedGaussian_ - use InterPotential_ + use GTFPotential_ implicit none #define contr(n,m) contractions(n)%contractions(m) @@ -253,7 +253,7 @@ subroutine G12Integrals_diskIntraSpecie(specieID) !Get potential ID do i=1, InterPotential_instance%ssize - if ( trim( MolecularSystem_instance%species(specieID)%symbol) == trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. trim( MolecularSystem_instance%species(specieID)%symbol) == trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie))) ) then + if ( trim( MolecularSystem_instance%species(specieID)%symbol) == trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. trim( MolecularSystem_instance%species(specieID)%symbol) == trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies))) ) then potID=i exit end if @@ -918,14 +918,14 @@ subroutine G12Integrals_G12diskInterSpecie(nameOfSpecie, otherNameOfSpecie, spe !Get potential ID do i=1, InterPotential_instance%ssize if ( (trim(MolecularSystem_instance%species(specieID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim(MolecularSystem_instance%species(otherSpecieID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie)) ) & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies)) ) & ) .or. & (trim( MolecularSystem_instance%species(otherSpecieID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim( MolecularSystem_instance%species(specieID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie)) ) & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies)) ) & ) & ) then potID=i diff --git a/src/ints/IntegralManager.f90 b/src/ints/IntegralManager.f90 index 4c669083..ebba8eed 100644 --- a/src/ints/IntegralManager.f90 +++ b/src/ints/IntegralManager.f90 @@ -37,7 +37,7 @@ module IntegralManager_ use Matrix_ use CosmoCore_ use Stopwatch_ - use ExternalPotential_ + use GTFPotential_ use HarmonicIntegrals_ use FirstDerivativeIntegrals_ @@ -929,8 +929,8 @@ subroutine IntegralManager_writeThreeCenterIntegralsByProduct() do i= 1, ExternalPotential_instance%ssize !if( trim(potential(i)%specie)==trim(interactNameSelected) ) then ! This does not work for UHF ! if ( String_findSubstring(trim( MolecularSystem_instance%species(f)%name ), & - ! trim(String_getUpperCase(trim(ExternalPotential_instance%potentials(i)%specie)))) == 1 ) then - if ( trim( MolecularSystem_instance%species(f)%symbol) == trim(String_getUpperCase(trim(ExternalPotential_instance%potentials(i)%specie))) ) then + ! trim(String_getUpperCase(trim(ExternalPotential_instance%potentials(i)%species)))) == 1 ) then + if ( trim( MolecularSystem_instance%species(f)%symbol) == trim(String_getUpperCase(trim(ExternalPotential_instance%potentials(i)%species))) ) then potID=i exit end if diff --git a/src/ints/Libint2Interface.f90 b/src/ints/Libint2Interface.f90 index 23f5ffa1..78cc9d13 100644 --- a/src/ints/Libint2Interface.f90 +++ b/src/ints/Libint2Interface.f90 @@ -26,7 +26,7 @@ module Libint2Interface_ use, intrinsic :: iso_c_binding use MolecularSystem_ - use InterPotential_ + use GTFPotential_ use ContractedGaussian_ ! use Matrix_ @@ -824,9 +824,9 @@ subroutine Libint2Interface_computeG12Intraspecies_disk(speciesID) !Get potential ID do i=1, InterPotential_instance%ssize if ( trim( MolecularSystem_instance%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim( MolecularSystem_instance%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie))) ) then + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies))) ) then potID=i exit end if @@ -883,14 +883,14 @@ subroutine Libint2Interface_computeG12Interspecies_disk(speciesID,otherSpeciesID !Get potential ID do i=1, InterPotential_instance%ssize if ( (trim(MolecularSystem_instance%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim(MolecularSystem_instance%species(otherSpeciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie)) ) & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies)) ) & ) .or. & (trim( MolecularSystem_instance%species(otherSpeciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim( MolecularSystem_instance%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie)) ) & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies)) ) & ) & ) then potID=i @@ -965,9 +965,9 @@ subroutine Libint2Interface_computeG12Intraspecies_direct(speciesID, density, tw !Get potential ID do i=1, InterPotential_instance%ssize if ( trim( molSys%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim( molSys%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie))) ) then + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies))) ) then potID=i exit end if @@ -1027,14 +1027,14 @@ subroutine Libint2Interface_computeG12Interspecies_direct(speciesID,otherSpecies !Get potential ID do i=1, InterPotential_instance%ssize if ( (trim(molSys%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim(molSys%species(otherSpeciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie)) ) & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies)) ) & ) .or. & (trim(molSys%species(otherSpeciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim(molSys%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie)) ) & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies)) ) & ) & ) then potID=i @@ -1088,9 +1088,9 @@ subroutine Libint2Interface_compute2BodyIntraspecies_direct_all(speciesID, densi !Get potential ID do i=1, InterPotential_instance%ssize if ( trim( molSys%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim( molSys%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie))) ) then + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies))) ) then potID=i exit end if @@ -1142,14 +1142,14 @@ subroutine Libint2Interface_compute2BodyInterspecies_direct_all(speciesID, other !Get potential ID do i=1, InterPotential_instance%ssize if ( (trim(molSys%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim(molSys%species(otherSpeciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie)) ) & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies)) ) & ) .or. & (trim( molSys%species(otherSpeciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%specie))) .and. & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%species))) .and. & trim( molSys%species(speciesID)%symbol) == & - trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecie)) ) & + trim(String_getUpperCase(trim(InterPotential_instance%potentials(i)%otherSpecies)) ) & ) & ) then potID=i diff --git a/src/scf/DensityMatrixSCFGuess.f90 b/src/scf/DensityMatrixSCFGuess.f90 index a64a9a3a..f3a279a3 100644 --- a/src/scf/DensityMatrixSCFGuess.f90 +++ b/src/scf/DensityMatrixSCFGuess.f90 @@ -54,7 +54,7 @@ subroutine DensityMatrixSCFGuess_getGuess( speciesID, hcoreMatrix, transformatio type(MolecularSystem), pointer :: molSys type(Matrix) :: auxMatrix - character(30) :: nameOfSpecies + character(30) :: nameOfSpecies, symbolOfSpecies integer(8) :: orderOfMatrix, occupationNumber logical :: existPlain, existBinnary, readSuccess character(50) :: guessType @@ -71,6 +71,7 @@ subroutine DensityMatrixSCFGuess_getGuess( speciesID, hcoreMatrix, transformatio orderOfMatrix = MolecularSystem_getTotalnumberOfContractions(speciesID,molSys) nameOfSpecies = molSys%species(speciesID)%name + symbolOfSpecies = molSys%species(speciesID)%symbol occupationNumber = MolecularSystem_getOcupationNumber(speciesID,molSys) readSuccess=.false. @@ -79,11 +80,12 @@ subroutine DensityMatrixSCFGuess_getGuess( speciesID, hcoreMatrix, transformatio call Matrix_constructor(densityMatrix, int(orderOfMatrix,8), int(orderOfMatrix,8), 0.0_8 ) call Matrix_constructor(orbitals, int(orderOfMatrix,8), int(orderOfMatrix,8), 0.0_8 ) - + + readSuccess=.false. !!Verifica el archivo que contiene los coeficientes para una especie dada if ( CONTROL_instance%READ_FCHK ) then - call MolecularSystem_readFchk(trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".fchk", orbitals, densityMatrix, nameOfSpecies ) - return + wfnFile=trim(CONTROL_instance%INPUT_FILE)//trim(symbolOfSpecies)//".fchk" + call MolecularSystem_readFchk(wfnFile, orbitals, densityMatrix, nameOfSpecies, readSuccess) else if ( CONTROL_instance%READ_COEFFICIENTS ) then wfnUnit = 30 @@ -107,10 +109,7 @@ subroutine DensityMatrixSCFGuess_getGuess( speciesID, hcoreMatrix, transformatio end if end if end if - - !check if the orbitals were read correctly - if(.not. allocated(orbitals%values) ) readSuccess=.false. - + if(readSuccess .and. printInfo ) print *, "Combination coefficients for ", trim(nameOfSpecies), " were read from ", trim(wfnFile) if(.not. readSuccess) then diff --git a/src/scf/OrbitalLocalizer.f90 b/src/scf/OrbitalLocalizer.f90 index 28e4ca88..006ab893 100644 --- a/src/scf/OrbitalLocalizer.f90 +++ b/src/scf/OrbitalLocalizer.f90 @@ -109,7 +109,7 @@ subroutine OrbitalLocalizer_erkaleLocal(speciesID,densityMatrix,fockMatrix,orbit type(Matrix) :: orbitalCoefficients type(Vector) :: orbitalEnergies - character(30) :: nameOfSpecies + character(30) :: nameOfSpecies, symbolOfSpecies integer :: statusSystem integer :: numberOfContractions @@ -118,11 +118,12 @@ subroutine OrbitalLocalizer_erkaleLocal(speciesID,densityMatrix,fockMatrix,orbit !! Convert lowdin fchk files to erkale chk files nameOfSpecies=MolecularSystem_getNameOfSpecies(speciesID) + symbolOfSpecies=MolecularSystem_getSymbolOfSpecies(speciesID) numberOfContractions = MolecularSystem_getTotalNumberOfContractions(speciesID) open(unit=30, file="erkale.read", status="replace", form="formatted") - write(30,*) "LoadFChk ", trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".fchk" - write(30,*) "SaveChk ", trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".chk" + write(30,*) "LoadFChk ", trim(CONTROL_instance%INPUT_FILE)//trim(symbolOfSpecies)//".fchk" + write(30,*) "SaveChk ", trim(CONTROL_instance%INPUT_FILE)//trim(symbolOfSpecies)//".chk" write(30,*) "Reorthonormalize true" close(30) @@ -132,8 +133,8 @@ subroutine OrbitalLocalizer_erkaleLocal(speciesID,densityMatrix,fockMatrix,orbit !! Localize orbitals open(unit=30, file="erkale.local", status="replace", form="formatted") - write(30,*) "LoadChk ", trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".chk" - write(30,*) "SaveChk ", trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".local.chk" + write(30,*) "LoadChk ", trim(CONTROL_instance%INPUT_FILE)//trim(symbolOfSpecies)//".chk" + write(30,*) "SaveChk ", trim(CONTROL_instance%INPUT_FILE)//trim(symbolOfSpecies)//".local.chk" write(30,*) "Method ", trim(CONTROL_instance%ERKALE_LOCALIZATION_METHOD) write(30,*) "Virtual false" write(30,*) "Maxiter 5000" @@ -149,15 +150,15 @@ subroutine OrbitalLocalizer_erkaleLocal(speciesID,densityMatrix,fockMatrix,orbit !!Convert erkale chk files to lowdin fchk files open(unit=30, file="erkale.write", status="replace", form="formatted") - write(30,*) "LoadChk ", trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".local.chk" - write(30,*) "SaveFChk ", trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".local.fchk" + write(30,*) "LoadChk ", trim(CONTROL_instance%INPUT_FILE)//trim(symbolOfSpecies)//".local.chk" + write(30,*) "SaveFChk ", trim(CONTROL_instance%INPUT_FILE)//trim(symbolOfSpecies)//".local.fchk" close(30) call system("erkale_fchkpt erkale.write") !! Read orbital coefficients from fchk files - call MolecularSystem_readFchk(trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".local.fchk", orbitalCoefficients, densityMatrix, nameOfSpecies ) + call MolecularSystem_readFchk(trim(CONTROL_instance%INPUT_FILE)//trim(symbolOfSpecies)//".local.fchk", orbitalCoefficients, densityMatrix, nameOfSpecies ) orbitalEnergies%values=0.0 !! Molecular orbital fock operator expected value diff --git a/src/scf/SingleSCF.f90 b/src/scf/SingleSCF.f90 index 6cb85fd1..16c219fb 100644 --- a/src/scf/SingleSCF.f90 +++ b/src/scf/SingleSCF.f90 @@ -508,7 +508,7 @@ subroutine SingleSCF_readFrozenOrbitals(wfObject) if (CONTROL_instance%READ_FCHK) then call Matrix_constructor (auxiliaryMatrix, numberOfContractions, numberOfContractions) - call MolecularSystem_readFchk( trim(CONTROL_instance%INPUT_FILE)//trim(wfObject%name)//".fchk", & + call MolecularSystem_readFchk( trim(CONTROL_instance%INPUT_FILE)//trim(wfObject%molSys%species(wfObject%species)%symbol)//".fchk", & wfObject%waveFunctionCoefficients, auxiliaryMatrix, wfObject%name ) call Matrix_destructor(auxiliaryMatrix) diff --git a/test/CHDTHemu.massTest.lowdin b/test/CHDTHemu.massTest.lowdin new file mode 100644 index 00000000..85916399 --- /dev/null +++ b/test/CHDTHemu.massTest.lowdin @@ -0,0 +1,23 @@ + +GEOMETRY +e-(C) NAKAI-CC-PVTZ 0.0000000 0.0000000 0.0000000 +e-(H) NAKAI-CC-PVTZ 0.6283310 0.6283310 0.6283310 +e-(H) NAKAI-CC-PVTZ -0.6283310 -0.6283310 0.6283310 +e-(H) NAKAI-CC-PVTZ -0.6283310 0.6283310 -0.6283310 +e-(H) NAKAI-CC-PVTZ 0.6283310 -0.6283310 -0.6283310 +U- HEMU 0.6283310 0.6283310 0.6283310 +C_12 NAKAI-3-SP 0.0000000 0.0000000 0.0000000 +He_4 NAKAI-3-SP 0.6283310 0.6283310 0.6283310 +H_1 DZSPNB -0.6283310 -0.6283310 0.6283310 +H_2 DZSPNB -0.6283310 0.6283310 -0.6283310 +H_3 DZSPNB 0.6283310 -0.6283310 -0.6283310 +END GEOMETRY + +TASKS + method = "RHF" +END TASKS + +CONTROL + readCoefficients=F + removeTranslationalContamination=T +END CONTROL diff --git a/test/CHDTHemu.massTest.py b/test/CHDTHemu.massTest.py new file mode 100644 index 00000000..e1d5f5a4 --- /dev/null +++ b/test/CHDTHemu.massTest.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python +from __future__ import print_function +import os +import sys +from colorstring import * + +if len(sys.argv)==2: + lowdinbin = sys.argv[1] +else: + lowdinbin = "lowdin2" + +testName = sys.argv[0][:-3] +inputName = testName + ".lowdin" +outputName = testName + ".out" +# Reference values and tolerance + +refValues = { + "Total mass" : [40383.2869, 5E-4], + "E- Kinetic energy" : [36.558549208752,1E-6], + "MUON Kinetic energy" : [4.823223401186,1E-6], + "C_12 Kinetic energy" : [0.020362665915,1E-6], + "HE_4 Kinetic energy" : [0.042798597334,1E-6], + "H_1 Kinetic energy" : [0.018810346458,1E-6], + "H_2 Kinetic energy" : [0.013654378944,1E-6], + "H_3 Kinetic energy" : [0.011128664979,1E-6], + "HF energy" : [-74.168958869716,1E-8] +} + +testValues = dict(refValues) #copy +for value in testValues: #reset + testValues[value] = 0 #reset + +# Run calculation + +status = os.system(lowdinbin + " -i " + inputName) + +if status: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + +output = open(outputName, "r") +outputRead = output.readlines() + +# Values +checkArray=[0,0,0,0] +for i in range(0,len(outputRead)): + line = outputRead[i] + if "TOTAL ENERGY =" in line: + testValues["HF energy"] = float(line.split()[3]) + if "MASS (m_e)" in line: + testValues["Total mass"] = float(line.split()[3]) + if "E- Kinetic energy" in line: + testValues["E- Kinetic energy"] = float(line.split()[4]) + if "MUON Kinetic energy" in line: + testValues["MUON Kinetic energy"] = float(line.split()[4]) + if "C_12 Kinetic energy" in line: + testValues["C_12 Kinetic energy"] = float(line.split()[4]) + if "HE_4 Kinetic energy" in line: + testValues["HE_4 Kinetic energy"] = float(line.split()[4]) + if "H_1 Kinetic energy" in line: + testValues["H_1 Kinetic energy"] = float(line.split()[4]) + if "H_2 Kinetic energy" in line: + testValues["H_2 Kinetic energy"] = float(line.split()[4]) + if "H_3 Kinetic energy" in line: + testValues["H_3 Kinetic energy"] = float(line.split()[4]) + +output.close() + +passTest = True + +for value in refValues: + diffValue = abs(refValues[value][0] - testValues[value]) + if ( diffValue <= refValues[value][1] ): + passTest = passTest * True + else : + passTest = passTest * False + print("%s %.8f %.8f %.2e" % ( value, refValues[value][0], testValues[value], diffValue)) + +if passTest : + print(testName + str_green(" ... OK")) +else: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + diff --git a/test/He2-C60potential-NOCI.lowdin b/test/He2-C60potential-NOCI.lowdin index c06d7379..4883911f 100644 --- a/test/He2-C60potential-NOCI.lowdin +++ b/test/He2-C60potential-NOCI.lowdin @@ -2,8 +2,8 @@ SYSTEM_DESCRIPTION='H' GEOMETRY N0 dirac 0.0 0.0 0.0 rotationPoint=1 -HEA3 HE2-1S 1.818207 0.0 -0.347193 rotateAround=1 -HEB3 HE2-1S -1.818207 0.0 0.347193 rotateAround=1 +HEA3 HE2-1S 1.818207 0.0 -0.347193 rotateAround=1 m=5494.8926 +HEB3 HE2-1S -1.818207 0.0 0.347193 rotateAround=1 m=5494.8926 END GEOMETRY TASKS diff --git a/test/He3-C60potential-NOCI.lowdin b/test/He3-C60potential-NOCI.lowdin index 16870c5d..a5efe877 100644 --- a/test/He3-C60potential-NOCI.lowdin +++ b/test/He3-C60potential-NOCI.lowdin @@ -2,9 +2,9 @@ SYSTEM_DESCRIPTION='H' GEOMETRY N0 dirac 0.0 0.0 0.0 rotationPoint=1 -HEA3 HE2-1S 2.159 0.0 0.0 rotateAround=1 -HEA3 HE2-1S 0.0 2.159 0.0 rotateAround=1 -HEB3 HE2-1S 0.0 0.0 2.159 rotateAround=1 +HEA3 HE2-1S 2.159 0.0 0.0 rotateAround=1 m=5494.8926 +HEA3 HE2-1S 0.0 2.159 0.0 rotateAround=1 m=5494.8926 +HEB3 HE2-1S 0.0 0.0 2.159 rotateAround=1 m=5494.8926 END GEOMETRY TASKS diff --git a/test/Ps2F2.RHF-customEta.lowdin b/test/Ps2F2.RHF-customEta.lowdin new file mode 100644 index 00000000..2484865d --- /dev/null +++ b/test/Ps2F2.RHF-customEta.lowdin @@ -0,0 +1,17 @@ +GEOMETRY + e-(F) aug-cc-pvdz 0.00 0.00 -1.33 addParticles=1 + e-(F) aug-cc-pvdz 0.00 0.00 1.33 addParticles=1 + e+ PSX-DZ 0.00 0.00 -1.33 eta=2 + e+ PSX-DZ 0.00 0.00 1.33 + F dirac 0.00 0.00 -1.33 + F dirac 0.00 0.00 1.33 +END GEOMETRY + +TASKS + method = "RHF" +END TASKS + +CONTROL + readCoefficients=F +END CONTROL + diff --git a/test/Ps2F2.RHF-customEta.py b/test/Ps2F2.RHF-customEta.py new file mode 100644 index 00000000..a812ee62 --- /dev/null +++ b/test/Ps2F2.RHF-customEta.py @@ -0,0 +1,68 @@ +#!/usr/bin/env python +from __future__ import print_function +import os +import sys +from colorstring import * + +if len(sys.argv)==2: + lowdinbin = sys.argv[1] +else: + lowdinbin = "lowdin2" + +testName = sys.argv[0][:-3] +inputName = testName + ".lowdin" +outputName = testName + ".out" + +# Reference values and tolerance +refValues = { + "HF energy" : [-199.213740198536,1E-8], + "eta e+" : [2.0,1E-8], + "occupation e+" : [1.0,1E-8] +} + +testValues = dict(refValues) #copy +for value in testValues: #reset + testValues[value] = 0 #reset + +# Run calculation + +status = os.system(lowdinbin + " -i " + inputName) + +if status: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + +output = open(outputName, "r") +outputRead = output.readlines() + +# Values +flagC=0 +for i in range(0,len(outputRead)): + line = outputRead[i] + if "TOTAL ENERGY =" in line: + testValues["HF energy"] = float(line.split()[3]) + if "CONSTANTS OF COUPLING" in line: + flagC=1 + if "E+" in line and flagC==1: + testValues["eta e+"] = float(line.split()[2]) + testValues["occupation e+"] = float(line.split()[4]) + flagC=0 +output.close() + +passTest = True + +for value in refValues: + diffValue = abs(refValues[value][0] - testValues[value]) + if ( diffValue <= refValues[value][1] ): + passTest = passTest * True + else : + passTest = passTest * False + print("%s %.8f %.8f %.2e" % ( value, refValues[value][0], testValues[value], diffValue)) + +if passTest : + print(testName + str_green(" ... OK")) +else: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + +output.close() diff --git a/test/TIP4P-dimer-singlet-NOCI.lowdin b/test/TIP4P-dimer-singlet-NOCI.lowdin index d886158d..213f2093 100644 --- a/test/TIP4P-dimer-singlet-NOCI.lowdin +++ b/test/TIP4P-dimer-singlet-NOCI.lowdin @@ -2,12 +2,12 @@ SYSTEM_DESCRIPTION='H' GEOMETRY N0 dirac 0 0 0 -X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 translationCenter=1 -Y0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 translationCenter=2 +X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 translationCenter=1 q=0.5564 m=1836.15267247 +Y0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 translationCenter=2 q=0.5564 m=1836.15267247 N0 dirac 0.0 0.0 6.0 -X1.1- dirac 0.0 0.0 6.292152 -X0.5+ dirac 0.0 1.430429 7.107157 -X0.5+ dirac 0.0 -1.430429 7.107157 +X1.1- dirac 0.0 0.0 6.292152 q=-1.1128 +X0.5+ dirac 0.0 1.430429 7.107157 q=0.5564 +X0.5+ dirac 0.0 -1.430429 7.107157 q=0.5564 END GEOMETRY TASKS @@ -36,3 +36,452 @@ INTERPOTENTIAL X0.5+ Y0.5+ VHH-CCSDT Y0.5+ Y0.5+ VHH-CCSDT END INTERPOTENTIAL + +EXTERPOTENTIAL + X0.5+ VOH-CCSDT + Y0.5+ VOH-CCSDT +END EXTERPOTENTIAL + +INTERPOTENTIAL + X0.5+ X0.5+ VHH-CCSDT + X0.5+ Y0.5+ VHH-CCSDT + Y0.5+ Y0.5+ VHH-CCSDT +END INTERPOTENTIAL + +BASIS H2O-1S1P1D +#optimized exponents +O-H-TIP X0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 + +O-H-TIP Y0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 +END BASIS + +POTENTIAL VOH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant RHH +O-X0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 + +O-Y0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 +END POTENTIAL + +POTENTIAL VHH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant ROH +O-X0.5+X0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 + +O-X0.5+Y0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 + +O-Y0.5+Y0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 +END POTENTIAL diff --git a/test/TIP4P-dimer-singlet-direct.lowdin b/test/TIP4P-dimer-singlet-direct.lowdin index 9e5e2eaa..3788cc16 100644 --- a/test/TIP4P-dimer-singlet-direct.lowdin +++ b/test/TIP4P-dimer-singlet-direct.lowdin @@ -1,11 +1,11 @@ GEOMETRY N0 dirac 0 0 0 -X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 -Y0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 +X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 q=0.5564 m=1836.15267247 +Y0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 q=0.5564 m=1836.15267247 N0 dirac 0.0 0.0 6.0 -X1.1- dirac 0.0 0.0 6.292152 -X0.5+ dirac 0.0 1.430429 7.107157 -X0.5+ dirac 0.0 -1.430429 7.107157 +X1.1- dirac 0.0 0.0 6.292152 q=-1.1128 +X0.5+ dirac 0.0 1.430429 7.107157 q=0.5564 +X0.5+ dirac 0.0 -1.430429 7.107157 q=0.5564 END GEOMETRY TASKS @@ -31,3 +31,440 @@ INTERPOTENTIAL Y0.5+ Y0.5+ VHH-CCSDT END INTERPOTENTIAL +BASIS H2O-1S1P1D +#optimized exponents +O-H-TIP X0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 + +O-H-TIP Y0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 +END BASIS + +POTENTIAL VOH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant RHH +O-X0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 + +O-Y0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 +END POTENTIAL + +POTENTIAL VHH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant ROH +O-X0.5+X0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 + +O-X0.5+Y0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 + +O-Y0.5+Y0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 +END POTENTIAL diff --git a/test/TIP4P-dimer-singlet-direct.py b/test/TIP4P-dimer-singlet-direct.py index f3b406ea..bc41eaa3 100644 --- a/test/TIP4P-dimer-singlet-direct.py +++ b/test/TIP4P-dimer-singlet-direct.py @@ -17,11 +17,11 @@ refValues = { "HF energy" : [-0.651377267238,1E-8], -"HA-TIP Ext Pot" : [0.005339817047,1E-4], -"HB-TIP Ext Pot" : [0.005265789260,1E-4], -"HA-TIP/HB-TIP Hartree" : [0.003393498016,1E-4], -"HA-TIP/Fixed interact." : [-0.025239485153,1E-4], -"HB-TIP/Fixed interact." : [-0.010156190638,1E-4] +"X0.5+ Ext Pot" : [0.005339817047,1E-4], +"Y0.5+ Ext Pot" : [0.005265789260,1E-4], +"X0.5+/Y0.5+ Hartree" : [0.003393498016,1E-4], +"X0.5+/Fixed interact." : [-0.025239485153,1E-4], +"Y0.5+/Fixed interact." : [-0.010156190638,1E-4] } testValues = dict(refValues) #copy @@ -44,16 +44,16 @@ line = outputRead[i] if "TOTAL ENERGY =" in line: testValues["HF energy"] = float(line.split()[3]) - if "HA-TIP Ext Pot" in line: - testValues["HA-TIP Ext Pot"] = float(line.split()[5]) - if "HB-TIP Ext Pot" in line: - testValues["HB-TIP Ext Pot"] = float(line.split()[5]) - if "HA-TIP/HB-TIP Hartree" in line: - testValues["HA-TIP/HB-TIP Hartree"] = float(line.split()[4]) - if "HA-TIP/Fixed interact." in line: - testValues["HA-TIP/Fixed interact."] = float(line.split()[4]) - if "HB-TIP/Fixed interact." in line: - testValues["HB-TIP/Fixed interact."] = float(line.split()[4]) + if "X0.5+ Ext Pot" in line: + testValues["X0.5+ Ext Pot"] = float(line.split()[5]) + if "Y0.5+ Ext Pot" in line: + testValues["Y0.5+ Ext Pot"] = float(line.split()[5]) + if "X0.5+/Y0.5+ Hartree" in line: + testValues["X0.5+/Y0.5+ Hartree"] = float(line.split()[4]) + if "X0.5+/Fixed interact." in line: + testValues["X0.5+/Fixed interact."] = float(line.split()[4]) + if "Y0.5+/Fixed interact." in line: + testValues["Y0.5+/Fixed interact."] = float(line.split()[4]) passTest = True diff --git a/test/TIP4P-dimer-singlet-memory.lowdin b/test/TIP4P-dimer-singlet-memory.lowdin index bc275bf3..1dc863f5 100644 --- a/test/TIP4P-dimer-singlet-memory.lowdin +++ b/test/TIP4P-dimer-singlet-memory.lowdin @@ -1,11 +1,11 @@ GEOMETRY N0 dirac 0 0 0 -X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 -Y0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 +X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 q=0.5564 m=1836.15267247 +Y0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 q=0.5564 m=1836.15267247 N0 dirac 0.0 0.0 6.0 -X1.1- dirac 0.0 0.0 6.292152 -X0.5+ dirac 0.0 1.430429 7.107157 -X0.5+ dirac 0.0 -1.430429 7.107157 +X1.1- dirac 0.0 0.0 6.292152 q=-1.1128 +X0.5+ dirac 0.0 1.430429 7.107157 q=0.5564 +X0.5+ dirac 0.0 -1.430429 7.107157 q=0.5564 END GEOMETRY TASKS @@ -31,3 +31,440 @@ INTERPOTENTIAL Y0.5+ Y0.5+ VHH-CCSDT END INTERPOTENTIAL +BASIS H2O-1S1P1D +#optimized exponents +O-H-TIP X0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 + +O-H-TIP Y0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 +END BASIS + +POTENTIAL VOH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant RHH +O-X0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 + +O-Y0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 +END POTENTIAL + +POTENTIAL VHH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant ROH +O-X0.5+X0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 + +O-X0.5+Y0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 + +O-Y0.5+Y0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 +END POTENTIAL diff --git a/test/TIP4P-dimer-singlet-memory.py b/test/TIP4P-dimer-singlet-memory.py index f3b406ea..bc41eaa3 100644 --- a/test/TIP4P-dimer-singlet-memory.py +++ b/test/TIP4P-dimer-singlet-memory.py @@ -17,11 +17,11 @@ refValues = { "HF energy" : [-0.651377267238,1E-8], -"HA-TIP Ext Pot" : [0.005339817047,1E-4], -"HB-TIP Ext Pot" : [0.005265789260,1E-4], -"HA-TIP/HB-TIP Hartree" : [0.003393498016,1E-4], -"HA-TIP/Fixed interact." : [-0.025239485153,1E-4], -"HB-TIP/Fixed interact." : [-0.010156190638,1E-4] +"X0.5+ Ext Pot" : [0.005339817047,1E-4], +"Y0.5+ Ext Pot" : [0.005265789260,1E-4], +"X0.5+/Y0.5+ Hartree" : [0.003393498016,1E-4], +"X0.5+/Fixed interact." : [-0.025239485153,1E-4], +"Y0.5+/Fixed interact." : [-0.010156190638,1E-4] } testValues = dict(refValues) #copy @@ -44,16 +44,16 @@ line = outputRead[i] if "TOTAL ENERGY =" in line: testValues["HF energy"] = float(line.split()[3]) - if "HA-TIP Ext Pot" in line: - testValues["HA-TIP Ext Pot"] = float(line.split()[5]) - if "HB-TIP Ext Pot" in line: - testValues["HB-TIP Ext Pot"] = float(line.split()[5]) - if "HA-TIP/HB-TIP Hartree" in line: - testValues["HA-TIP/HB-TIP Hartree"] = float(line.split()[4]) - if "HA-TIP/Fixed interact." in line: - testValues["HA-TIP/Fixed interact."] = float(line.split()[4]) - if "HB-TIP/Fixed interact." in line: - testValues["HB-TIP/Fixed interact."] = float(line.split()[4]) + if "X0.5+ Ext Pot" in line: + testValues["X0.5+ Ext Pot"] = float(line.split()[5]) + if "Y0.5+ Ext Pot" in line: + testValues["Y0.5+ Ext Pot"] = float(line.split()[5]) + if "X0.5+/Y0.5+ Hartree" in line: + testValues["X0.5+/Y0.5+ Hartree"] = float(line.split()[4]) + if "X0.5+/Fixed interact." in line: + testValues["X0.5+/Fixed interact."] = float(line.split()[4]) + if "Y0.5+/Fixed interact." in line: + testValues["Y0.5+/Fixed interact."] = float(line.split()[4]) passTest = True diff --git a/test/TIP4P-dimer-singlet.lowdin b/test/TIP4P-dimer-singlet.lowdin index 9d540025..7ffc1cfe 100644 --- a/test/TIP4P-dimer-singlet.lowdin +++ b/test/TIP4P-dimer-singlet.lowdin @@ -1,11 +1,11 @@ GEOMETRY N0 dirac 0 0 0 -X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 -Y0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 +X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 q=0.5564 m=1836.15267247 +Y0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 q=0.5564 m=1836.15267247 N0 dirac 0.0 0.0 6.0 -X1.1- dirac 0.0 0.0 6.292152 -X0.5+ dirac 0.0 1.430429 7.107157 -X0.5+ dirac 0.0 -1.430429 7.107157 +X1.1- dirac 0.0 0.0 6.292152 q=-1.1128 +X0.5+ dirac 0.0 1.430429 7.107157 q=0.5564 +X0.5+ dirac 0.0 -1.430429 7.107157 q=0.5564 END GEOMETRY TASKS @@ -31,3 +31,440 @@ INTERPOTENTIAL Y0.5+ Y0.5+ VHH-CCSDT END INTERPOTENTIAL +BASIS H2O-1S1P1D +#optimized exponents +O-H-TIP X0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 + +O-H-TIP Y0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 +END BASIS + +POTENTIAL VOH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant RHH +O-X0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 + +O-Y0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 +END POTENTIAL + +POTENTIAL VHH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant ROH +O-X0.5+X0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 + +O-X0.5+Y0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 + +O-Y0.5+Y0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 +END POTENTIAL diff --git a/test/TIP4P-dimer-singlet.py b/test/TIP4P-dimer-singlet.py index f3b406ea..bc41eaa3 100644 --- a/test/TIP4P-dimer-singlet.py +++ b/test/TIP4P-dimer-singlet.py @@ -17,11 +17,11 @@ refValues = { "HF energy" : [-0.651377267238,1E-8], -"HA-TIP Ext Pot" : [0.005339817047,1E-4], -"HB-TIP Ext Pot" : [0.005265789260,1E-4], -"HA-TIP/HB-TIP Hartree" : [0.003393498016,1E-4], -"HA-TIP/Fixed interact." : [-0.025239485153,1E-4], -"HB-TIP/Fixed interact." : [-0.010156190638,1E-4] +"X0.5+ Ext Pot" : [0.005339817047,1E-4], +"Y0.5+ Ext Pot" : [0.005265789260,1E-4], +"X0.5+/Y0.5+ Hartree" : [0.003393498016,1E-4], +"X0.5+/Fixed interact." : [-0.025239485153,1E-4], +"Y0.5+/Fixed interact." : [-0.010156190638,1E-4] } testValues = dict(refValues) #copy @@ -44,16 +44,16 @@ line = outputRead[i] if "TOTAL ENERGY =" in line: testValues["HF energy"] = float(line.split()[3]) - if "HA-TIP Ext Pot" in line: - testValues["HA-TIP Ext Pot"] = float(line.split()[5]) - if "HB-TIP Ext Pot" in line: - testValues["HB-TIP Ext Pot"] = float(line.split()[5]) - if "HA-TIP/HB-TIP Hartree" in line: - testValues["HA-TIP/HB-TIP Hartree"] = float(line.split()[4]) - if "HA-TIP/Fixed interact." in line: - testValues["HA-TIP/Fixed interact."] = float(line.split()[4]) - if "HB-TIP/Fixed interact." in line: - testValues["HB-TIP/Fixed interact."] = float(line.split()[4]) + if "X0.5+ Ext Pot" in line: + testValues["X0.5+ Ext Pot"] = float(line.split()[5]) + if "Y0.5+ Ext Pot" in line: + testValues["Y0.5+ Ext Pot"] = float(line.split()[5]) + if "X0.5+/Y0.5+ Hartree" in line: + testValues["X0.5+/Y0.5+ Hartree"] = float(line.split()[4]) + if "X0.5+/Fixed interact." in line: + testValues["X0.5+/Fixed interact."] = float(line.split()[4]) + if "Y0.5+/Fixed interact." in line: + testValues["Y0.5+/Fixed interact."] = float(line.split()[4]) passTest = True diff --git a/test/TIP4P-dimer-triplet-NOCI.lowdin b/test/TIP4P-dimer-triplet-NOCI.lowdin index b8411434..8a88242a 100644 --- a/test/TIP4P-dimer-triplet-NOCI.lowdin +++ b/test/TIP4P-dimer-triplet-NOCI.lowdin @@ -2,12 +2,12 @@ SYSTEM_DESCRIPTION='H' GEOMETRY N0 dirac 0 0 0 -X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 translationCenter=1 -X0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 translationCenter=2 +X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 translationCenter=1 q=0.5564 m=1836.15267247 +X0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 translationCenter=2 q=0.5564 m=1836.15267247 N0 dirac 0.0 0.0 6.0 -X1.1- dirac 0.0 0.0 6.292152 -X0.5+ dirac 0.0 1.430429 7.107157 -X0.5+ dirac 0.0 -1.430429 7.107157 +X1.1- dirac 0.0 0.0 6.292152 q=-1.1128 +X0.5+ dirac 0.0 1.430429 7.107157 q=0.5564 +X0.5+ dirac 0.0 -1.430429 7.107157 q=0.5564 END GEOMETRY TASKS @@ -33,3 +33,441 @@ END EXTERPOTENTIAL INTERPOTENTIAL X0.5+ X0.5+ VHH-CCSDT END INTERPOTENTIAL + +BASIS H2O-1S1P1D +#optimized exponents +O-H-TIP X0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 + +O-H-TIP Y0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 +END BASIS + +POTENTIAL VOH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant RHH +O-X0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 + +O-Y0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 +END POTENTIAL + +POTENTIAL VHH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant ROH +O-X0.5+X0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 + +O-X0.5+Y0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 + +O-Y0.5+Y0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 +END POTENTIAL diff --git a/test/TIP4P-dimer-triplet-direct.lowdin b/test/TIP4P-dimer-triplet-direct.lowdin index 15ab0abb..2db120e6 100644 --- a/test/TIP4P-dimer-triplet-direct.lowdin +++ b/test/TIP4P-dimer-triplet-direct.lowdin @@ -1,11 +1,11 @@ GEOMETRY N0 dirac 0 0 0 -X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 -X0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 +X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 q=0.5564 m=1836.15267247 +X0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 q=0.5564 m=1836.15267247 N0 dirac 0.0 0.0 6.0 -X1.1- dirac 0.0 0.0 6.292152 -X0.5+ dirac 0.0 1.430429 7.107157 -X0.5+ dirac 0.0 -1.430429 7.107157 +X1.1- dirac 0.0 0.0 6.292152 q=-1.1128 +X0.5+ dirac 0.0 1.430429 7.107157 q=0.5564 +X0.5+ dirac 0.0 -1.430429 7.107157 q=0.5564 END GEOMETRY TASKS @@ -28,3 +28,186 @@ INTERPOTENTIAL X0.5+ X0.5+ VHH-CCSDT END INTERPOTENTIAL + +BASIS H2O-1S1P1D +#optimized exponents +O-H-TIP X0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 +END BASIS + +POTENTIAL VOH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant RHH +O-X0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 +END POTENTIAL + +POTENTIAL VHH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant ROH +O-X0.5+X0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 +END POTENTIAL diff --git a/test/TIP4P-dimer-triplet-direct.py b/test/TIP4P-dimer-triplet-direct.py index 41f2ebc4..d13b500b 100644 --- a/test/TIP4P-dimer-triplet-direct.py +++ b/test/TIP4P-dimer-triplet-direct.py @@ -17,10 +17,10 @@ refValues = { "HF energy" : [-0.651377261423,1E-8], - "HA-TIP Ext Pot" : [0.010606635479,1E-4], - "HA-TIP/HA-TIP Hartree" : [1.217813835967,1E-4], - "HA-TIP Exchange" : [-1.214419735313,1E-4], - "HA-TIP/Fixed interact." : [-0.035396418562,1E-4] + "X0.5+ Ext Pot" : [0.010606635479,1E-4], + "X0.5+/X0.5+ Hartree" : [1.217813835967,1E-4], + "X0.5+ Exchange" : [-1.214419735313,1E-4], + "X0.5+/Fixed interact." : [-0.035396418562,1E-4] } testValues = dict(refValues) #copy @@ -43,14 +43,14 @@ line = outputRead[i] if "TOTAL ENERGY =" in line: testValues["HF energy"] = float(line.split()[3]) - if "HA-TIP Ext Pot" in line: - testValues["HA-TIP Ext Pot"] = float(line.split()[5]) - if "HA-TIP/HA-TIP Hartree" in line: - testValues["HA-TIP/HA-TIP Hartree"] = float(line.split()[4]) - if "HA-TIP Exchange" in line: - testValues["HA-TIP Exchange"] = float(line.split()[4]) - if "HA-TIP/Fixed interact." in line: - testValues["HA-TIP/Fixed interact."] = float(line.split()[4]) + if "X0.5+ Ext Pot" in line: + testValues["X0.5+ Ext Pot"] = float(line.split()[5]) + if "X0.5+/X0.5+ Hartree" in line: + testValues["X0.5+/X0.5+ Hartree"] = float(line.split()[4]) + if "X0.5+ Exchange" in line: + testValues["X0.5+ Exchange"] = float(line.split()[4]) + if "X0.5+/Fixed interact." in line: + testValues["X0.5+/Fixed interact."] = float(line.split()[4]) passTest = True diff --git a/test/TIP4P-dimer-triplet-memory.lowdin b/test/TIP4P-dimer-triplet-memory.lowdin index c37e1e48..5d36490e 100644 --- a/test/TIP4P-dimer-triplet-memory.lowdin +++ b/test/TIP4P-dimer-triplet-memory.lowdin @@ -1,11 +1,11 @@ GEOMETRY N0 dirac 0 0 0 -X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 -X0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 +X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 q=0.5564 m=1836.15267247 +X0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 q=0.5564 m=1836.15267247 N0 dirac 0.0 0.0 6.0 -X1.1- dirac 0.0 0.0 6.292152 -X0.5+ dirac 0.0 1.430429 7.107157 -X0.5+ dirac 0.0 -1.430429 7.107157 +X1.1- dirac 0.0 0.0 6.292152 q=-1.1128 +X0.5+ dirac 0.0 1.430429 7.107157 q=0.5564 +X0.5+ dirac 0.0 -1.430429 7.107157 q=0.5564 END GEOMETRY TASKS @@ -28,3 +28,194 @@ INTERPOTENTIAL X0.5+ X0.5+ VHH-CCSDT END INTERPOTENTIAL +EXTERPOTENTIAL + X0.5+ VOH-CCSDT +END EXTERPOTENTIAL + +INTERPOTENTIAL + X0.5+ X0.5+ VHH-CCSDT +END INTERPOTENTIAL + +BASIS H2O-1S1P1D +#optimized exponents +O-H-TIP X0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 +END BASIS + +POTENTIAL VOH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant RHH +O-X0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 + +END POTENTIAL + +POTENTIAL VHH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant ROH +O-X0.5+X0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 +END POTENTIAL diff --git a/test/TIP4P-dimer-triplet-memory.py b/test/TIP4P-dimer-triplet-memory.py index 41f2ebc4..d13b500b 100644 --- a/test/TIP4P-dimer-triplet-memory.py +++ b/test/TIP4P-dimer-triplet-memory.py @@ -17,10 +17,10 @@ refValues = { "HF energy" : [-0.651377261423,1E-8], - "HA-TIP Ext Pot" : [0.010606635479,1E-4], - "HA-TIP/HA-TIP Hartree" : [1.217813835967,1E-4], - "HA-TIP Exchange" : [-1.214419735313,1E-4], - "HA-TIP/Fixed interact." : [-0.035396418562,1E-4] + "X0.5+ Ext Pot" : [0.010606635479,1E-4], + "X0.5+/X0.5+ Hartree" : [1.217813835967,1E-4], + "X0.5+ Exchange" : [-1.214419735313,1E-4], + "X0.5+/Fixed interact." : [-0.035396418562,1E-4] } testValues = dict(refValues) #copy @@ -43,14 +43,14 @@ line = outputRead[i] if "TOTAL ENERGY =" in line: testValues["HF energy"] = float(line.split()[3]) - if "HA-TIP Ext Pot" in line: - testValues["HA-TIP Ext Pot"] = float(line.split()[5]) - if "HA-TIP/HA-TIP Hartree" in line: - testValues["HA-TIP/HA-TIP Hartree"] = float(line.split()[4]) - if "HA-TIP Exchange" in line: - testValues["HA-TIP Exchange"] = float(line.split()[4]) - if "HA-TIP/Fixed interact." in line: - testValues["HA-TIP/Fixed interact."] = float(line.split()[4]) + if "X0.5+ Ext Pot" in line: + testValues["X0.5+ Ext Pot"] = float(line.split()[5]) + if "X0.5+/X0.5+ Hartree" in line: + testValues["X0.5+/X0.5+ Hartree"] = float(line.split()[4]) + if "X0.5+ Exchange" in line: + testValues["X0.5+ Exchange"] = float(line.split()[4]) + if "X0.5+/Fixed interact." in line: + testValues["X0.5+/Fixed interact."] = float(line.split()[4]) passTest = True diff --git a/test/TIP4P-dimer-triplet.lowdin b/test/TIP4P-dimer-triplet.lowdin index 75996651..18fd85c8 100644 --- a/test/TIP4P-dimer-triplet.lowdin +++ b/test/TIP4P-dimer-triplet.lowdin @@ -1,11 +1,11 @@ GEOMETRY N0 dirac 0 0 0 -X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 -X0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 +X0.5+ H2O-1S1P1D -0.151399 0.000276 1.807200 q=0.5564 m=1836.15267247 +X0.5+ H2O-1S1P1D -1.715032 0.000843 -0.589570 q=0.5564 m=1836.15267247 N0 dirac 0.0 0.0 6.0 -X1.1- dirac 0.0 0.0 6.292152 -X0.5+ dirac 0.0 1.430429 7.107157 -X0.5+ dirac 0.0 -1.430429 7.107157 +X1.1- dirac 0.0 0.0 6.292152 q=-1.1128 +X0.5+ dirac 0.0 1.430429 7.107157 q=0.5564 +X0.5+ dirac 0.0 -1.430429 7.107157 q=0.5564 END GEOMETRY TASKS @@ -28,3 +28,186 @@ INTERPOTENTIAL X0.5+ X0.5+ VHH-CCSDT END INTERPOTENTIAL + +BASIS H2O-1S1P1D +#optimized exponents +O-H-TIP X0.5+ (1S) BASIS TYPE: 2 +3 +1 0 1 +14.509888498676842 1.0 +2 1 1 +6.885507269761004 1.0 +3 2 1 +9.023681376783887 1.0 +END BASIS + +POTENTIAL VOH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant RHH +O-X0.5+ +27 +1 0 +13.892336977 9.839539381 +0.0 0.0 0.0 +2 0 +10.000000000 15.916259785 +0.0 0.0 0.0 +3 0 +7.498942093 -18.199008237 +0.0 0.0 0.0 +4 0 +5.623413252 6.723316202 +0.0 0.0 0.0 +5 0 +4.216965034 4.671467422 +0.0 0.0 0.0 +6 0 +3.162277660 1.703771884 +0.0 0.0 0.0 +7 0 +2.371373706 0.297676262 +0.0 0.0 0.0 +8 0 +1.778279410 0.967600283 +0.0 0.0 0.0 +9 0 +1.333521432 1.034932150 +0.0 0.0 0.0 +10 0 +1.000000000 0.657621305 +0.0 0.0 0.0 +11 0 +0.749894209 -0.327381267 +0.0 0.0 0.0 +12 0 +0.562341325 0.221528399 +0.0 0.0 0.0 +13 0 +0.421696503 -0.255694044 +0.0 0.0 0.0 +14 0 +0.316227766 0.097853405 +0.0 0.0 0.0 +15 0 +0.237137371 0.744489008 +0.0 0.0 0.0 +16 0 +0.177827941 -0.750090202 +0.0 0.0 0.0 +17 0 +0.133352143 -0.309090719 +0.0 0.0 0.0 +18 0 +0.100000000 -0.998800172 +0.0 0.0 0.0 +19 0 +0.074989421 0.718795407 +0.0 0.0 0.0 +20 0 +0.056234133 0.228002288 +0.0 0.0 0.0 +21 0 +0.042169650 0.573744431 +0.0 0.0 0.0 +22 0 +0.031622777 -0.056969092 +0.0 0.0 0.0 +23 0 +0.023713737 -0.392241803 +0.0 0.0 0.0 +24 0 +0.017782794 -0.569609681 +0.0 0.0 0.0 +25 0 +0.013335214 0.759323484 +0.0 0.0 0.0 +26 0 +0.010000000 -0.194525832 +0.0 0.0 0.0 +27 0 +0.000000000 0.138912412 +0.0 0.0 0.0 +END POTENTIAL + +POTENTIAL VHH-CCSDT +#Fitted from CCSD(T)/def2-TZVPPD with constant ROH +O-X0.5+X0.5+ +26 +1 0 +10.000000000 1.313644267 +0.0 0.0 0.0 +2 0 +7.498942093 3.312208762 +0.0 0.0 0.0 +3 0 +5.623413252 -4.290599715 +0.0 0.0 0.0 +4 0 +4.216965034 -2.400490962 +0.0 0.0 0.0 +5 0 +3.162277660 6.160643793 +0.0 0.0 0.0 +6 0 +2.371373706 2.406172066 +0.0 0.0 0.0 +7 0 +1.778279410 -6.852617537 +0.0 0.0 0.0 +8 0 +1.333521432 -1.043562213 +0.0 0.0 0.0 +9 0 +1.000000000 7.556790649 +0.0 0.0 0.0 +10 0 +0.749894209 0.777037995 +0.0 0.0 0.0 +11 0 +0.562341325 -8.976316536 +0.0 0.0 0.0 +12 0 +0.421696503 0.575296040 +0.0 0.0 0.0 +13 0 +0.316227766 8.672468936 +0.0 0.0 0.0 +14 0 +0.237137371 -0.855388714 +0.0 0.0 0.0 +15 0 +0.177827941 -6.271012179 +0.0 0.0 0.0 +16 0 +0.133352143 1.152328020 +0.0 0.0 0.0 +17 0 +0.100000000 3.403807143 +0.0 0.0 0.0 +18 0 +0.074989421 -1.745300340 +0.0 0.0 0.0 +19 0 +0.056234133 -1.349326021 +0.0 0.0 0.0 +20 0 +0.042169650 1.605561812 +0.0 0.0 0.0 +21 0 +0.031622777 -0.267553445 +0.0 0.0 0.0 +22 0 +0.023713737 0.034566428 +0.0 0.0 0.0 +23 0 +0.017782794 -0.185356609 +0.0 0.0 0.0 +24 0 +0.013335214 -0.103128899 +0.0 0.0 0.0 +25 0 +0.010000000 0.128087663 +0.0 0.0 0.0 +26 0 +0.000000000 0.085213897 +0.0 0.0 0.0 +END POTENTIAL diff --git a/test/TIP4P-dimer-triplet.py b/test/TIP4P-dimer-triplet.py index 41f2ebc4..d13b500b 100644 --- a/test/TIP4P-dimer-triplet.py +++ b/test/TIP4P-dimer-triplet.py @@ -17,10 +17,10 @@ refValues = { "HF energy" : [-0.651377261423,1E-8], - "HA-TIP Ext Pot" : [0.010606635479,1E-4], - "HA-TIP/HA-TIP Hartree" : [1.217813835967,1E-4], - "HA-TIP Exchange" : [-1.214419735313,1E-4], - "HA-TIP/Fixed interact." : [-0.035396418562,1E-4] + "X0.5+ Ext Pot" : [0.010606635479,1E-4], + "X0.5+/X0.5+ Hartree" : [1.217813835967,1E-4], + "X0.5+ Exchange" : [-1.214419735313,1E-4], + "X0.5+/Fixed interact." : [-0.035396418562,1E-4] } testValues = dict(refValues) #copy @@ -43,14 +43,14 @@ line = outputRead[i] if "TOTAL ENERGY =" in line: testValues["HF energy"] = float(line.split()[3]) - if "HA-TIP Ext Pot" in line: - testValues["HA-TIP Ext Pot"] = float(line.split()[5]) - if "HA-TIP/HA-TIP Hartree" in line: - testValues["HA-TIP/HA-TIP Hartree"] = float(line.split()[4]) - if "HA-TIP Exchange" in line: - testValues["HA-TIP Exchange"] = float(line.split()[4]) - if "HA-TIP/Fixed interact." in line: - testValues["HA-TIP/Fixed interact."] = float(line.split()[4]) + if "X0.5+ Ext Pot" in line: + testValues["X0.5+ Ext Pot"] = float(line.split()[5]) + if "X0.5+/X0.5+ Hartree" in line: + testValues["X0.5+/X0.5+ Hartree"] = float(line.split()[4]) + if "X0.5+ Exchange" in line: + testValues["X0.5+ Exchange"] = float(line.split()[4]) + if "X0.5+/Fixed interact." in line: + testValues["X0.5+/Fixed interact."] = float(line.split()[4]) passTest = True From 71f5936054555910c6cec9035a296b19cf3f0904 Mon Sep 17 00:00:00 2001 From: Felix Moncada Date: Wed, 27 Nov 2024 18:00:23 +0100 Subject: [PATCH 2/3] Fixed test broken by the default mass change of He_4 --- test/HemuH-CUSTOM_BASIS.lowdin | 6 +++--- test/HemuH-CUSTOM_BASIS.py | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/test/HemuH-CUSTOM_BASIS.lowdin b/test/HemuH-CUSTOM_BASIS.lowdin index 19b21087..27c3692a 100644 --- a/test/HemuH-CUSTOM_BASIS.lowdin +++ b/test/HemuH-CUSTOM_BASIS.lowdin @@ -2,9 +2,9 @@ GEOMETRY e-[H] cc-pvtz 0.0000 0.0000 0.00000 e-[H] CUSTOM_1 0.0000 0.0000 0.74144 - H_1 CUSTOM_2 0.0000 0.0000 0.00000 - U- CUSTOM_3 0.0000 0.0000 0.74144 - He_4 CUSTOM_3 0.0000 0.0000 0.74144 + H_1 CUSTOM_2 0.0000 0.0000 0.00000 m=1836.1527 + U- CUSTOM_3 0.0000 0.0000 0.74144 m=206.7683 + He_4 CUSTOM_3 0.0000 0.0000 0.74144 m=7349.6727 END GEOMETRY TASKS diff --git a/test/HemuH-CUSTOM_BASIS.py b/test/HemuH-CUSTOM_BASIS.py index 7d461f17..44c69e98 100644 --- a/test/HemuH-CUSTOM_BASIS.py +++ b/test/HemuH-CUSTOM_BASIS.py @@ -15,11 +15,11 @@ # Reference values and tolerance refValues = { -"HF energy" : [-343.383191892820,1E-8], -"U-HOMO" : [-371.890049816287,1E-1], -"H_1-HOMO" : [-1.019360160964,1E-4], -"He_4-HOMO" : [-652.366876763392,1E-1], -"e-HOMO" : [-0.585414450602,1E-4], +"HF energy" : [-343.383218426575,1E-8], +"U-HOMO" : [-371.889981903188,1E-1], +"H_1-HOMO" : [-1.019346186407,1E-4], +"He_4-HOMO" : [-652.365841581870,1E-1], +"e-HOMO" : [-0.585408097570,1E-4], } testValues = dict(refValues) #copy From b3d256f100111586a7a03ebec52de14340bc55a1 Mon Sep 17 00:00:00 2001 From: Felix Moncada Date: Wed, 27 Nov 2024 18:15:24 +0100 Subject: [PATCH 3/3] The custom potentials were not working The tests worked because I didnt remove the potential files from the lib/ folder so it read those and got the correct result --- bin/lowdin | 2 +- lib/potentials/VHH-CCSDT | 567 --------------------------------------- lib/potentials/VOH-CCSDT | 420 ----------------------------- 3 files changed, 1 insertion(+), 988 deletions(-) delete mode 100644 lib/potentials/VHH-CCSDT delete mode 100644 lib/potentials/VOH-CCSDT diff --git a/bin/lowdin b/bin/lowdin index 6377f7d1..da82a9d3 100755 --- a/bin/lowdin +++ b/bin/lowdin @@ -415,7 +415,7 @@ if [ $extFile="lowdin" ]; then echo "Modify the POTENTIALS block in your input and select a different name" exit 1 fi - gawk '($1~/POTENTIALS/ && toupper($2)~/^'$POTENTIALS_NAME'$/){flag=1; next} + gawk '($1~/^POTENTIAL$/ && toupper($2)~/^'$POTENTIALS_NAME'$/){flag=1; next} ($0~/END/){flag=0}; (flag==1){print toupper($0)}' $nameFile > $POTENTIALS_NAME.$PID done diff --git a/lib/potentials/VHH-CCSDT b/lib/potentials/VHH-CCSDT deleted file mode 100644 index 138f0597..00000000 --- a/lib/potentials/VHH-CCSDT +++ /dev/null @@ -1,567 +0,0 @@ -#Fitted from CCSD(T)/def2-TZVPPD with constant ROH -O-H_1H_1 -26 -1 0 -10.000000000 1.313644267 -0.0 0.0 0.0 -2 0 -7.498942093 3.312208762 -0.0 0.0 0.0 -3 0 -5.623413252 -4.290599715 -0.0 0.0 0.0 -4 0 -4.216965034 -2.400490962 -0.0 0.0 0.0 -5 0 -3.162277660 6.160643793 -0.0 0.0 0.0 -6 0 -2.371373706 2.406172066 -0.0 0.0 0.0 -7 0 -1.778279410 -6.852617537 -0.0 0.0 0.0 -8 0 -1.333521432 -1.043562213 -0.0 0.0 0.0 -9 0 -1.000000000 7.556790649 -0.0 0.0 0.0 -10 0 -0.749894209 0.777037995 -0.0 0.0 0.0 -11 0 -0.562341325 -8.976316536 -0.0 0.0 0.0 -12 0 -0.421696503 0.575296040 -0.0 0.0 0.0 -13 0 -0.316227766 8.672468936 -0.0 0.0 0.0 -14 0 -0.237137371 -0.855388714 -0.0 0.0 0.0 -15 0 -0.177827941 -6.271012179 -0.0 0.0 0.0 -16 0 -0.133352143 1.152328020 -0.0 0.0 0.0 -17 0 -0.100000000 3.403807143 -0.0 0.0 0.0 -18 0 -0.074989421 -1.745300340 -0.0 0.0 0.0 -19 0 -0.056234133 -1.349326021 -0.0 0.0 0.0 -20 0 -0.042169650 1.605561812 -0.0 0.0 0.0 -21 0 -0.031622777 -0.267553445 -0.0 0.0 0.0 -22 0 -0.023713737 0.034566428 -0.0 0.0 0.0 -23 0 -0.017782794 -0.185356609 -0.0 0.0 0.0 -24 0 -0.013335214 -0.103128899 -0.0 0.0 0.0 -25 0 -0.010000000 0.128087663 -0.0 0.0 0.0 -26 0 -0.000000000 0.085213897 -0.0 0.0 0.0 - -O-H-A_1H-A_1 -26 -1 0 -10.000000000 1.313644267 -0.0 0.0 0.0 -2 0 -7.498942093 3.312208762 -0.0 0.0 0.0 -3 0 -5.623413252 -4.290599715 -0.0 0.0 0.0 -4 0 -4.216965034 -2.400490962 -0.0 0.0 0.0 -5 0 -3.162277660 6.160643793 -0.0 0.0 0.0 -6 0 -2.371373706 2.406172066 -0.0 0.0 0.0 -7 0 -1.778279410 -6.852617537 -0.0 0.0 0.0 -8 0 -1.333521432 -1.043562213 -0.0 0.0 0.0 -9 0 -1.000000000 7.556790649 -0.0 0.0 0.0 -10 0 -0.749894209 0.777037995 -0.0 0.0 0.0 -11 0 -0.562341325 -8.976316536 -0.0 0.0 0.0 -12 0 -0.421696503 0.575296040 -0.0 0.0 0.0 -13 0 -0.316227766 8.672468936 -0.0 0.0 0.0 -14 0 -0.237137371 -0.855388714 -0.0 0.0 0.0 -15 0 -0.177827941 -6.271012179 -0.0 0.0 0.0 -16 0 -0.133352143 1.152328020 -0.0 0.0 0.0 -17 0 -0.100000000 3.403807143 -0.0 0.0 0.0 -18 0 -0.074989421 -1.745300340 -0.0 0.0 0.0 -19 0 -0.056234133 -1.349326021 -0.0 0.0 0.0 -20 0 -0.042169650 1.605561812 -0.0 0.0 0.0 -21 0 -0.031622777 -0.267553445 -0.0 0.0 0.0 -22 0 -0.023713737 0.034566428 -0.0 0.0 0.0 -23 0 -0.017782794 -0.185356609 -0.0 0.0 0.0 -24 0 -0.013335214 -0.103128899 -0.0 0.0 0.0 -25 0 -0.010000000 0.128087663 -0.0 0.0 0.0 -26 0 -0.000000000 0.085213897 -0.0 0.0 0.0 - -O-H-A_1H-B_1 -26 -1 0 -10.000000000 1.313644267 -0.0 0.0 0.0 -2 0 -7.498942093 3.312208762 -0.0 0.0 0.0 -3 0 -5.623413252 -4.290599715 -0.0 0.0 0.0 -4 0 -4.216965034 -2.400490962 -0.0 0.0 0.0 -5 0 -3.162277660 6.160643793 -0.0 0.0 0.0 -6 0 -2.371373706 2.406172066 -0.0 0.0 0.0 -7 0 -1.778279410 -6.852617537 -0.0 0.0 0.0 -8 0 -1.333521432 -1.043562213 -0.0 0.0 0.0 -9 0 -1.000000000 7.556790649 -0.0 0.0 0.0 -10 0 -0.749894209 0.777037995 -0.0 0.0 0.0 -11 0 -0.562341325 -8.976316536 -0.0 0.0 0.0 -12 0 -0.421696503 0.575296040 -0.0 0.0 0.0 -13 0 -0.316227766 8.672468936 -0.0 0.0 0.0 -14 0 -0.237137371 -0.855388714 -0.0 0.0 0.0 -15 0 -0.177827941 -6.271012179 -0.0 0.0 0.0 -16 0 -0.133352143 1.152328020 -0.0 0.0 0.0 -17 0 -0.100000000 3.403807143 -0.0 0.0 0.0 -18 0 -0.074989421 -1.745300340 -0.0 0.0 0.0 -19 0 -0.056234133 -1.349326021 -0.0 0.0 0.0 -20 0 -0.042169650 1.605561812 -0.0 0.0 0.0 -21 0 -0.031622777 -0.267553445 -0.0 0.0 0.0 -22 0 -0.023713737 0.034566428 -0.0 0.0 0.0 -23 0 -0.017782794 -0.185356609 -0.0 0.0 0.0 -24 0 -0.013335214 -0.103128899 -0.0 0.0 0.0 -25 0 -0.010000000 0.128087663 -0.0 0.0 0.0 -26 0 -0.000000000 0.085213897 -0.0 0.0 0.0 - -O-H-B_1H-B_1 -26 -1 0 -10.000000000 1.313644267 -0.0 0.0 0.0 -2 0 -7.498942093 3.312208762 -0.0 0.0 0.0 -3 0 -5.623413252 -4.290599715 -0.0 0.0 0.0 -4 0 -4.216965034 -2.400490962 -0.0 0.0 0.0 -5 0 -3.162277660 6.160643793 -0.0 0.0 0.0 -6 0 -2.371373706 2.406172066 -0.0 0.0 0.0 -7 0 -1.778279410 -6.852617537 -0.0 0.0 0.0 -8 0 -1.333521432 -1.043562213 -0.0 0.0 0.0 -9 0 -1.000000000 7.556790649 -0.0 0.0 0.0 -10 0 -0.749894209 0.777037995 -0.0 0.0 0.0 -11 0 -0.562341325 -8.976316536 -0.0 0.0 0.0 -12 0 -0.421696503 0.575296040 -0.0 0.0 0.0 -13 0 -0.316227766 8.672468936 -0.0 0.0 0.0 -14 0 -0.237137371 -0.855388714 -0.0 0.0 0.0 -15 0 -0.177827941 -6.271012179 -0.0 0.0 0.0 -16 0 -0.133352143 1.152328020 -0.0 0.0 0.0 -17 0 -0.100000000 3.403807143 -0.0 0.0 0.0 -18 0 -0.074989421 -1.745300340 -0.0 0.0 0.0 -19 0 -0.056234133 -1.349326021 -0.0 0.0 0.0 -20 0 -0.042169650 1.605561812 -0.0 0.0 0.0 -21 0 -0.031622777 -0.267553445 -0.0 0.0 0.0 -22 0 -0.023713737 0.034566428 -0.0 0.0 0.0 -23 0 -0.017782794 -0.185356609 -0.0 0.0 0.0 -24 0 -0.013335214 -0.103128899 -0.0 0.0 0.0 -25 0 -0.010000000 0.128087663 -0.0 0.0 0.0 -26 0 -0.000000000 0.085213897 -0.0 0.0 0.0 - -O-X0.5+X0.5+ -26 -1 0 -10.000000000 1.313644267 -0.0 0.0 0.0 -2 0 -7.498942093 3.312208762 -0.0 0.0 0.0 -3 0 -5.623413252 -4.290599715 -0.0 0.0 0.0 -4 0 -4.216965034 -2.400490962 -0.0 0.0 0.0 -5 0 -3.162277660 6.160643793 -0.0 0.0 0.0 -6 0 -2.371373706 2.406172066 -0.0 0.0 0.0 -7 0 -1.778279410 -6.852617537 -0.0 0.0 0.0 -8 0 -1.333521432 -1.043562213 -0.0 0.0 0.0 -9 0 -1.000000000 7.556790649 -0.0 0.0 0.0 -10 0 -0.749894209 0.777037995 -0.0 0.0 0.0 -11 0 -0.562341325 -8.976316536 -0.0 0.0 0.0 -12 0 -0.421696503 0.575296040 -0.0 0.0 0.0 -13 0 -0.316227766 8.672468936 -0.0 0.0 0.0 -14 0 -0.237137371 -0.855388714 -0.0 0.0 0.0 -15 0 -0.177827941 -6.271012179 -0.0 0.0 0.0 -16 0 -0.133352143 1.152328020 -0.0 0.0 0.0 -17 0 -0.100000000 3.403807143 -0.0 0.0 0.0 -18 0 -0.074989421 -1.745300340 -0.0 0.0 0.0 -19 0 -0.056234133 -1.349326021 -0.0 0.0 0.0 -20 0 -0.042169650 1.605561812 -0.0 0.0 0.0 -21 0 -0.031622777 -0.267553445 -0.0 0.0 0.0 -22 0 -0.023713737 0.034566428 -0.0 0.0 0.0 -23 0 -0.017782794 -0.185356609 -0.0 0.0 0.0 -24 0 -0.013335214 -0.103128899 -0.0 0.0 0.0 -25 0 -0.010000000 0.128087663 -0.0 0.0 0.0 -26 0 -0.000000000 0.085213897 -0.0 0.0 0.0 - -O-X0.5+Y0.5+ -26 -1 0 -10.000000000 1.313644267 -0.0 0.0 0.0 -2 0 -7.498942093 3.312208762 -0.0 0.0 0.0 -3 0 -5.623413252 -4.290599715 -0.0 0.0 0.0 -4 0 -4.216965034 -2.400490962 -0.0 0.0 0.0 -5 0 -3.162277660 6.160643793 -0.0 0.0 0.0 -6 0 -2.371373706 2.406172066 -0.0 0.0 0.0 -7 0 -1.778279410 -6.852617537 -0.0 0.0 0.0 -8 0 -1.333521432 -1.043562213 -0.0 0.0 0.0 -9 0 -1.000000000 7.556790649 -0.0 0.0 0.0 -10 0 -0.749894209 0.777037995 -0.0 0.0 0.0 -11 0 -0.562341325 -8.976316536 -0.0 0.0 0.0 -12 0 -0.421696503 0.575296040 -0.0 0.0 0.0 -13 0 -0.316227766 8.672468936 -0.0 0.0 0.0 -14 0 -0.237137371 -0.855388714 -0.0 0.0 0.0 -15 0 -0.177827941 -6.271012179 -0.0 0.0 0.0 -16 0 -0.133352143 1.152328020 -0.0 0.0 0.0 -17 0 -0.100000000 3.403807143 -0.0 0.0 0.0 -18 0 -0.074989421 -1.745300340 -0.0 0.0 0.0 -19 0 -0.056234133 -1.349326021 -0.0 0.0 0.0 -20 0 -0.042169650 1.605561812 -0.0 0.0 0.0 -21 0 -0.031622777 -0.267553445 -0.0 0.0 0.0 -22 0 -0.023713737 0.034566428 -0.0 0.0 0.0 -23 0 -0.017782794 -0.185356609 -0.0 0.0 0.0 -24 0 -0.013335214 -0.103128899 -0.0 0.0 0.0 -25 0 -0.010000000 0.128087663 -0.0 0.0 0.0 -26 0 -0.000000000 0.085213897 -0.0 0.0 0.0 - -O-Y0.5+Y0.5+ -26 -1 0 -10.000000000 1.313644267 -0.0 0.0 0.0 -2 0 -7.498942093 3.312208762 -0.0 0.0 0.0 -3 0 -5.623413252 -4.290599715 -0.0 0.0 0.0 -4 0 -4.216965034 -2.400490962 -0.0 0.0 0.0 -5 0 -3.162277660 6.160643793 -0.0 0.0 0.0 -6 0 -2.371373706 2.406172066 -0.0 0.0 0.0 -7 0 -1.778279410 -6.852617537 -0.0 0.0 0.0 -8 0 -1.333521432 -1.043562213 -0.0 0.0 0.0 -9 0 -1.000000000 7.556790649 -0.0 0.0 0.0 -10 0 -0.749894209 0.777037995 -0.0 0.0 0.0 -11 0 -0.562341325 -8.976316536 -0.0 0.0 0.0 -12 0 -0.421696503 0.575296040 -0.0 0.0 0.0 -13 0 -0.316227766 8.672468936 -0.0 0.0 0.0 -14 0 -0.237137371 -0.855388714 -0.0 0.0 0.0 -15 0 -0.177827941 -6.271012179 -0.0 0.0 0.0 -16 0 -0.133352143 1.152328020 -0.0 0.0 0.0 -17 0 -0.100000000 3.403807143 -0.0 0.0 0.0 -18 0 -0.074989421 -1.745300340 -0.0 0.0 0.0 -19 0 -0.056234133 -1.349326021 -0.0 0.0 0.0 -20 0 -0.042169650 1.605561812 -0.0 0.0 0.0 -21 0 -0.031622777 -0.267553445 -0.0 0.0 0.0 -22 0 -0.023713737 0.034566428 -0.0 0.0 0.0 -23 0 -0.017782794 -0.185356609 -0.0 0.0 0.0 -24 0 -0.013335214 -0.103128899 -0.0 0.0 0.0 -25 0 -0.010000000 0.128087663 -0.0 0.0 0.0 -26 0 -0.000000000 0.085213897 -0.0 0.0 0.0 diff --git a/lib/potentials/VOH-CCSDT b/lib/potentials/VOH-CCSDT deleted file mode 100644 index d3b76b4d..00000000 --- a/lib/potentials/VOH-CCSDT +++ /dev/null @@ -1,420 +0,0 @@ -#Fitted from CCSD(T)/def2-TZVPPD with constant RHH -O-H-A_1 -27 -1 0 -13.892336977 9.839539381 -0.0 0.0 0.0 -2 0 -10.000000000 15.916259785 -0.0 0.0 0.0 -3 0 -7.498942093 -18.199008237 -0.0 0.0 0.0 -4 0 -5.623413252 6.723316202 -0.0 0.0 0.0 -5 0 -4.216965034 4.671467422 -0.0 0.0 0.0 -6 0 -3.162277660 1.703771884 -0.0 0.0 0.0 -7 0 -2.371373706 0.297676262 -0.0 0.0 0.0 -8 0 -1.778279410 0.967600283 -0.0 0.0 0.0 -9 0 -1.333521432 1.034932150 -0.0 0.0 0.0 -10 0 -1.000000000 0.657621305 -0.0 0.0 0.0 -11 0 -0.749894209 -0.327381267 -0.0 0.0 0.0 -12 0 -0.562341325 0.221528399 -0.0 0.0 0.0 -13 0 -0.421696503 -0.255694044 -0.0 0.0 0.0 -14 0 -0.316227766 0.097853405 -0.0 0.0 0.0 -15 0 -0.237137371 0.744489008 -0.0 0.0 0.0 -16 0 -0.177827941 -0.750090202 -0.0 0.0 0.0 -17 0 -0.133352143 -0.309090719 -0.0 0.0 0.0 -18 0 -0.100000000 -0.998800172 -0.0 0.0 0.0 -19 0 -0.074989421 0.718795407 -0.0 0.0 0.0 -20 0 -0.056234133 0.228002288 -0.0 0.0 0.0 -21 0 -0.042169650 0.573744431 -0.0 0.0 0.0 -22 0 -0.031622777 -0.056969092 -0.0 0.0 0.0 -23 0 -0.023713737 -0.392241803 -0.0 0.0 0.0 -24 0 -0.017782794 -0.569609681 -0.0 0.0 0.0 -25 0 -0.013335214 0.759323484 -0.0 0.0 0.0 -26 0 -0.010000000 -0.194525832 -0.0 0.0 0.0 -27 0 -0.000000000 0.138912412 -0.0 0.0 0.0 - -O-H-B_1 -27 -1 0 -13.892336977 9.839539381 -0.0 0.0 0.0 -2 0 -10.000000000 15.916259785 -0.0 0.0 0.0 -3 0 -7.498942093 -18.199008237 -0.0 0.0 0.0 -4 0 -5.623413252 6.723316202 -0.0 0.0 0.0 -5 0 -4.216965034 4.671467422 -0.0 0.0 0.0 -6 0 -3.162277660 1.703771884 -0.0 0.0 0.0 -7 0 -2.371373706 0.297676262 -0.0 0.0 0.0 -8 0 -1.778279410 0.967600283 -0.0 0.0 0.0 -9 0 -1.333521432 1.034932150 -0.0 0.0 0.0 -10 0 -1.000000000 0.657621305 -0.0 0.0 0.0 -11 0 -0.749894209 -0.327381267 -0.0 0.0 0.0 -12 0 -0.562341325 0.221528399 -0.0 0.0 0.0 -13 0 -0.421696503 -0.255694044 -0.0 0.0 0.0 -14 0 -0.316227766 0.097853405 -0.0 0.0 0.0 -15 0 -0.237137371 0.744489008 -0.0 0.0 0.0 -16 0 -0.177827941 -0.750090202 -0.0 0.0 0.0 -17 0 -0.133352143 -0.309090719 -0.0 0.0 0.0 -18 0 -0.100000000 -0.998800172 -0.0 0.0 0.0 -19 0 -0.074989421 0.718795407 -0.0 0.0 0.0 -20 0 -0.056234133 0.228002288 -0.0 0.0 0.0 -21 0 -0.042169650 0.573744431 -0.0 0.0 0.0 -22 0 -0.031622777 -0.056969092 -0.0 0.0 0.0 -23 0 -0.023713737 -0.392241803 -0.0 0.0 0.0 -24 0 -0.017782794 -0.569609681 -0.0 0.0 0.0 -25 0 -0.013335214 0.759323484 -0.0 0.0 0.0 -26 0 -0.010000000 -0.194525832 -0.0 0.0 0.0 -27 0 -0.000000000 0.138912412 -0.0 0.0 0.0 - -O-H_1 -27 -1 0 -13.892336977 9.839539381 -0.0 0.0 0.0 -2 0 -10.000000000 15.916259785 -0.0 0.0 0.0 -3 0 -7.498942093 -18.199008237 -0.0 0.0 0.0 -4 0 -5.623413252 6.723316202 -0.0 0.0 0.0 -5 0 -4.216965034 4.671467422 -0.0 0.0 0.0 -6 0 -3.162277660 1.703771884 -0.0 0.0 0.0 -7 0 -2.371373706 0.297676262 -0.0 0.0 0.0 -8 0 -1.778279410 0.967600283 -0.0 0.0 0.0 -9 0 -1.333521432 1.034932150 -0.0 0.0 0.0 -10 0 -1.000000000 0.657621305 -0.0 0.0 0.0 -11 0 -0.749894209 -0.327381267 -0.0 0.0 0.0 -12 0 -0.562341325 0.221528399 -0.0 0.0 0.0 -13 0 -0.421696503 -0.255694044 -0.0 0.0 0.0 -14 0 -0.316227766 0.097853405 -0.0 0.0 0.0 -15 0 -0.237137371 0.744489008 -0.0 0.0 0.0 -16 0 -0.177827941 -0.750090202 -0.0 0.0 0.0 -17 0 -0.133352143 -0.309090719 -0.0 0.0 0.0 -18 0 -0.100000000 -0.998800172 -0.0 0.0 0.0 -19 0 -0.074989421 0.718795407 -0.0 0.0 0.0 -20 0 -0.056234133 0.228002288 -0.0 0.0 0.0 -21 0 -0.042169650 0.573744431 -0.0 0.0 0.0 -22 0 -0.031622777 -0.056969092 -0.0 0.0 0.0 -23 0 -0.023713737 -0.392241803 -0.0 0.0 0.0 -24 0 -0.017782794 -0.569609681 -0.0 0.0 0.0 -25 0 -0.013335214 0.759323484 -0.0 0.0 0.0 -26 0 -0.010000000 -0.194525832 -0.0 0.0 0.0 -27 0 -0.000000000 0.138912412 -0.0 0.0 0.0 - -O-X0.5+ -27 -1 0 -13.892336977 9.839539381 -0.0 0.0 0.0 -2 0 -10.000000000 15.916259785 -0.0 0.0 0.0 -3 0 -7.498942093 -18.199008237 -0.0 0.0 0.0 -4 0 -5.623413252 6.723316202 -0.0 0.0 0.0 -5 0 -4.216965034 4.671467422 -0.0 0.0 0.0 -6 0 -3.162277660 1.703771884 -0.0 0.0 0.0 -7 0 -2.371373706 0.297676262 -0.0 0.0 0.0 -8 0 -1.778279410 0.967600283 -0.0 0.0 0.0 -9 0 -1.333521432 1.034932150 -0.0 0.0 0.0 -10 0 -1.000000000 0.657621305 -0.0 0.0 0.0 -11 0 -0.749894209 -0.327381267 -0.0 0.0 0.0 -12 0 -0.562341325 0.221528399 -0.0 0.0 0.0 -13 0 -0.421696503 -0.255694044 -0.0 0.0 0.0 -14 0 -0.316227766 0.097853405 -0.0 0.0 0.0 -15 0 -0.237137371 0.744489008 -0.0 0.0 0.0 -16 0 -0.177827941 -0.750090202 -0.0 0.0 0.0 -17 0 -0.133352143 -0.309090719 -0.0 0.0 0.0 -18 0 -0.100000000 -0.998800172 -0.0 0.0 0.0 -19 0 -0.074989421 0.718795407 -0.0 0.0 0.0 -20 0 -0.056234133 0.228002288 -0.0 0.0 0.0 -21 0 -0.042169650 0.573744431 -0.0 0.0 0.0 -22 0 -0.031622777 -0.056969092 -0.0 0.0 0.0 -23 0 -0.023713737 -0.392241803 -0.0 0.0 0.0 -24 0 -0.017782794 -0.569609681 -0.0 0.0 0.0 -25 0 -0.013335214 0.759323484 -0.0 0.0 0.0 -26 0 -0.010000000 -0.194525832 -0.0 0.0 0.0 -27 0 -0.000000000 0.138912412 -0.0 0.0 0.0 - -O-Y0.5+ -27 -1 0 -13.892336977 9.839539381 -0.0 0.0 0.0 -2 0 -10.000000000 15.916259785 -0.0 0.0 0.0 -3 0 -7.498942093 -18.199008237 -0.0 0.0 0.0 -4 0 -5.623413252 6.723316202 -0.0 0.0 0.0 -5 0 -4.216965034 4.671467422 -0.0 0.0 0.0 -6 0 -3.162277660 1.703771884 -0.0 0.0 0.0 -7 0 -2.371373706 0.297676262 -0.0 0.0 0.0 -8 0 -1.778279410 0.967600283 -0.0 0.0 0.0 -9 0 -1.333521432 1.034932150 -0.0 0.0 0.0 -10 0 -1.000000000 0.657621305 -0.0 0.0 0.0 -11 0 -0.749894209 -0.327381267 -0.0 0.0 0.0 -12 0 -0.562341325 0.221528399 -0.0 0.0 0.0 -13 0 -0.421696503 -0.255694044 -0.0 0.0 0.0 -14 0 -0.316227766 0.097853405 -0.0 0.0 0.0 -15 0 -0.237137371 0.744489008 -0.0 0.0 0.0 -16 0 -0.177827941 -0.750090202 -0.0 0.0 0.0 -17 0 -0.133352143 -0.309090719 -0.0 0.0 0.0 -18 0 -0.100000000 -0.998800172 -0.0 0.0 0.0 -19 0 -0.074989421 0.718795407 -0.0 0.0 0.0 -20 0 -0.056234133 0.228002288 -0.0 0.0 0.0 -21 0 -0.042169650 0.573744431 -0.0 0.0 0.0 -22 0 -0.031622777 -0.056969092 -0.0 0.0 0.0 -23 0 -0.023713737 -0.392241803 -0.0 0.0 0.0 -24 0 -0.017782794 -0.569609681 -0.0 0.0 0.0 -25 0 -0.013335214 0.759323484 -0.0 0.0 0.0 -26 0 -0.010000000 -0.194525832 -0.0 0.0 0.0 -27 0 -0.000000000 0.138912412 -0.0 0.0 0.0